%%%%%%%%%%%%
%%%%%%%%%%%%
%%% Generate a multiresolution variable-contrast grating in order to test contrast sensitivity function
%%%%%%%%%%%%
%%%%%%%%%%%%

%%% (c) Pietro Perona, California Institute of Technology. January 21, 2002


%%% Parameters of display being used
DISPLAY.COLS = 1024;
DISPLAY.ROWS = 768;
GRATING.PERIOD = [3 200];
GRATING.CONTRAST = [0.002 1];

%% Create the grating

%% Vertical component
M = log(GRATING.CONTRAST(2));
m = log(GRATING.CONTRAST(1));
DC = (M-m)/(DISPLAY.ROWS-1); %% Contrast increment
IY = exp([m:DC:M]');
figure(1); plot(IY);

%% Horizontal component
m = log(1/GRATING.PERIOD(2));
M = log(1/GRATING.PERIOD(1));
DP = (M-m)/(DISPLAY.COLS-1); %% Period increment
P = exp(m:DP:M);
IX = sin([1:DISPLAY.COLS].*P);
figure(2); plot(IX);
figure(3); plot(P)

%% Final image
I = 1/2 + (1-IY*ones(1,DISPLAY.COLS))/2 + IY * IX /2;
figure(4); image(256*I); colormap(gray(256));
imwrite(I,'grating.jpg','jpg');

gmap = gray(256);
rgmap = [gmap(:,1) 0.5*(1-gmap(:,1)) 0*gmap(:,1)]; %% Dim the green a bit to make it equiluminant with red

Irg = 1/2 + IY * IX /2;
figure(5); image(256*Irg); colormap(rgmap);
imwrite(round(Irg*255+1),rgmap,'rg_grating.jpg','jpg');
imwrite(round(Irg*255+1),gmap,'g_grating.jpg','jpg');  %% Another version of the graylevel image