function [nn,img,xax,yax] = makePinwheelMap(nVox,voxSize,nAngles,pixpermm)

if ~exist('pixpermm','var')
    pixpermm = 25;
end

if ~exist('nAngles','var')
    nAngles = 8;
end

%initial parameters:
FOV = nVox*voxSize; %Field of view (mm)
pixpermm = 25;  %Number of pixels per millimeter

sz = pixpermm*FOV;  %Total size of image (in pixels)

angleList = linspace(-180,180,nAngles+1);   %this determines the histogram bins
%9 numbers define 8 bins


%Rojer and Schwartz' method of bandpassing random noise:

%Rojer, A.S. and E.L. Schwartz, Cat and monkey cortical columnar patterns
%modeled by bandpass-filtered 2D white noise. Biol Cybern, 1990. 62(5): p. 381-91.

%Make random noise: complex numbers where the angle is the orientation
z = exp(sqrt(-1)*rand(sz)*pi*2);

%Make a gabor filter
filtSz = 3; % 3mm
sig = .8; %  .8 mm
freq = 1; %cycles/mm (Try zero for big columns)
filtPix = ceil(filtSz*pixpermm);
[x,y] = meshgrid(linspace(-filtSz/2,filtSz/2,filtPix),linspace(-filtSz/2,filtSz/2,filtPix));
r = sqrt(x.^2+y.^2);
filt = exp(-r.^2/sig.^2).*cos(2*pi*freq*r);  %Gabor

%Convolve z with the filter
bandpassz = conv2(z,filt,'same');

%scale it from 0 to 64
img = 64*(angle(bandpassz)+pi)/(pi*2);

%These are vectors for x and y axes, in units of millemters:
xax = linspace(.5/pixpermm,FOV(2)-.5/pixpermm,FOV(2)*pixpermm);
yax = linspace(.5/pixpermm,FOV(1)-.5/pixpermm,FOV(1)*pixpermm);


nn = zeros(nAngles,nVox(1),nVox(2));
%loop through the voxels
for i=1:nVox(1)  %columns
    for j=1:nVox(2) %rows
        %define out the pixels inside the voxel
        xx = (j-1)*voxSize;
        yy = (i-1)*voxSize;
        inVoxel = bandpassz(yy*pixpermm+1:(yy+voxSize)*pixpermm,xx*pixpermm+1:(xx+voxSize)*pixpermm);
        a = 180*angle(inVoxel(:))/pi;
        %calculate the pixels for each bin
        for k=1:nAngles
            nn(k,i,j) = sum(a>=angleList(k) & a<angleList(k+1));
        end
    end
end
%scale
nn = nn*(length(angleList)-1)/length(inVoxel(:));

nn = reshape(nn,nAngles,prod(nVox));