% pinwheel_synthesis.m:  generate synthetic orientation map
%
%  This script generates synthetic orientation map data by band-pass
%  filtering uniformly distributed "orientations" represented by unit
%  vectors with random angles between 0 and 180 degrees, following the model
%  of Rojer and Schwartz [1,2]. The constructed map is similar to the
%  patterns found in visual cortex via optical recording [3,4] and
%  electrophysiology [5].
%
%  The following code is a simplified version of the code used to generate
%  the figures of Polimeni et al. [6], and can be found reproduced in Figure
%  9 of the Supporting Text document.
%
%  A electronic version of [6], along with the accompanying Supporting Materials
%  and an electronic version of this script, can be found at
%  http://eslab.bu.edu/publications/polimeni2005physical/
%
%
%   References:
%   ===========
%
%  [1] A. Rojer and E. L. Schwartz. Cat and monkey cortical columnar
%      patterns modeled by bandpass-filtered 2D white noise. Biological
%      Cybernetics, 62:381-391, 1990.
%
%  [2] E. L. Schwartz and A. S. Rojer. Cortical hypercolumns and the
%      topology of random orientation maps. Proceedings of 12th
%      International Conference on Pattern Recognition, pages 150-5
%      vol.2, 1994.
%
%  [3] G. G. Blasdel and G. Salama. Voltage-sensitive dyes reveal a
%      modular organization in monkey striate cortex. Nature,
%      321(6070):579-585, 5-11 June 1986.
%
%  [4] T. Bonhoeffer and A. Grinvald. The layout of iso-orientation
%      domains in area 18 of cat visual cortex: optical imaging reveals
%      a pinwheel-like organization. Journal of Neuroscience,
%      13(10):4157-4180, October 1993.
%
%  [5] N. V. Swindale, J. A. Matsubara, and M. S. Cynader. Surface
%      organization of orientation and direction selectivity in cat area
%      18. Journal of Neuroscience, 7(5):1414-1427, May 1987.
%
%  [6] J. R. Polimeni, D. Granquist-Fraser, R. J. Wood, and E. L. Schwartz.
%      Physical limits to spatial resolution of optical recording:
%      Clarifying the spatial structure of cortical hypercolumns.
%      Proceedings of the National Academy of Sciences of the United States
%      of America, 102(11):4158-4163, 15 March 2005.

% Copyright (C) 2003, 2004, 2005 Jonathan Polimeni <jonnyreb@cns.bu.edu>
%   Computer Vision and Computational Neuroscience Lab
%   Department of Cognitive and Neural Systems
%   Boston University
%   Boston, MA  02215
%
% This program is free software; you can redistribute it and/or modify it
% under the terms of the GNU General Public License as published by the Free
% Software Foundation; either version 2 of the License, or (at your option)
% any later version.
%
% This program is distributed in the hope that it will be useful, but
% WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
% or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
% for more details.
%
% You should have received a copy of the GNU General Public License along
% with this program; if not, write to the Free Software Foundation, Inc., 59
% Temple Place - Suite 330, Boston, MA 02111-1307, USA.
%
% $Author: jonnyreb $
% $Date: 2005/03/13 03:09:11 $
%========================================================================%


% image size (assumed square)
N = 256;

% generate uniform random angle noise
v = rand(N,N).*(2*pi);

% bandpass filter parameters
fwidth = 0.025;
fpass1 = 0.070;
fpass2 = fpass1 + fwidth;
forder = 200;

% use the MATLAB FIR machinery to make a nice bandpass filter
H = fftshift(freqz2(ftrans2(fir1(forder, [fpass1 fpass2])), [N N]));

% bandpass filtering using FFT method for convolution
x = real(ifft2( fft2(cos(v)).*H ));
y = real(ifft2( fft2(sin(v)).*H ));

% compute angle map from real and imaginary parts
phi = atan2(y, x);

% perform phase unwrapping such that 0 <= phi < 2pi
neg = find(phi < 0);
phi(neg) = phi(neg) + 2*pi;

% divide by 2 to relabel "angles" as "orientations", i.e., 0 <= phi < pi
phi = phi / 2;

% open new figure and display
figure; imagesc(phi); axis image;