Absorptions Code

%% Create harmonic sequences % % Make the images of harmonics % Compute the cone mosaic responses and photocurrent responses % Show how to control eye movements % %

%% Init Parameters

   ieInit;

%% Here is a new function: oisCreate

% We are going to mix a grating (harmonic) with a uniform background. % These are the relative weights of the grating and background. % Set up the timing

% Gaussian onset and offset of the grating

   tSeries = ieScale(fspecial('gaussian',[1,100],10),0,1);

% % Have a look at the stimulus amplitude time series % vcNewGraphWin; plot(stimWeights);

% This is the field of view of the scene.

   sparams.fov = 0.5;
   sparams.meanluminance = 200;

% Initialize the harmonic parameters structure with default % Change entries that are common to uniform and harmonic

   clear params
   for ii=2:-1:1
       params(ii) = harmonicP; 
       params(ii).GaborFlag = 0.15;
       params(ii).freq      = 30;
   end

% params(1) is for the uniform field

   params(1).contrast  = 0.5;  % contrast of the two frequencies

% params(2) is matched and describes the grating

   params(2).contrast = 0.4;
   len = length(tSeries);

% The call to create the optical image sequence

   oisH = oisCreate('harmonic','blend',tSeries(len/2:len/2+1),'testParameters',params,'sceneParameters',sparams);

% Have a look, though the code here is not that great yet so the look is % only approximate.

   oisH.visualize;

%% Now, make the cone mosaic and compute absorptions and current

   fov = oiGet(oisH.oiFixed,'fov');
   emSamples = oisH.length;
   cMosaic = coneMosaic;
   cMosaic.noiseFlag = 'none';
   cMosaic.integrationTime = 0.001;
   cMosaic.setSizeToFOV(0.5*fov);

% The number of eye movement samples should extend as long as the oisH % So these should be equal % % tSamples * cMosaic.integrationTime = oisH.length*oisH.timeStep % %

   tSamples = floor(oisH.length*oisH.timeStep/cMosaic.integrationTime);
   cMosaic.emGenSequence(tSamples);

% Compute and then look

   cMosaic.compute(oisH);
   cMosaic.window;
   fprintf('Spatial frequency %.1f cpd\n',params(ii).freq/sparams.fov);

%% Create the current with and without noise

% cMosaic.os.noiseFlag = true;

   cMosaic.computeCurrent;
   cMosaic.window;

%% Compute and look at the current % cMosaic.os.noiseFlag = false;

   cMosaic.computeCurrent;
   cMosaic.window;

%% How to adjust the eye movement parameters % create simple sinusoidal eye movement

   emF = 3; emA = 3;
   x = round(emA*sin(2*pi*emF*(1:tSamples)/tSamples));
   y = zeros(size(x(:)));
   cMosaic.emPositions = [x(:),y(:)];
   [abs_move, curr_move] = cMosaic.compute(oisH,'currentFlag',true);
   abs_move = squeeze(abs_move);
   curr_move = squeeze(curr_move);
   movement_pos = cMosaic.emPositions;
   cMosaic.plot('eye movement path');  % How did we do?
   set(gca,'xlim',[-15 15],'ylim',[-15 15]);

% no tremor

   em_noMove = emCreate;     % Create an eye movement object
   em_noMove.emFlag = [0 0 0];  % Make sure tremor, draft and saccade are all on
   em_noMove.tremor.amplitude = 0.02;  % Set the big amplitude
   cMosaic.emGenSequence(tSamples,'em',em_noMove);  % Generate the sequence
   [abs_noMove, curr_noMove] = cMosaic.compute(oisH,'currentFlag',true);
   abs_noMove = squeeze(abs_noMove);
   curr_noMove = squeeze(curr_noMove);
   noMove_pos = cMosaic.emPositions;
   cMosaic.plot('eye movement path');  % How did we do?
   set(gca,'xlim',[-15 15],'ylim',[-15 15]);

%% plot movement

   figure;
   hold on;
   plot(1:length(movement_pos),movement_pos,'b');
   plot(1:length(noMove_pos),noMove_pos,'r');
   hold off;
   title('eye movement over time');
   xlabel('time (ms)');
   ylabel('distance (cone width)');

%% check to see if absorptions were actually computed differently % t = 1:size(abs_move,3); % plot(t,squeeze(abs_noMove(45,45,:)),t,squeeze(abs_move(45,45,:)));

%% absorption level comparison

   H = sum(abs_move,3)./sum(abs_noMove,3);
   figure();
   imagesc(H);
   colorbar;
   title('Ratio of total absorptions (with mvmnt / without mvmnt)');

% figure(); % imagesc(sum(abs_noMove,3)); avg_ratio = sum(sum(H))/numel(H)

%% Separate cone types % calculate absorption sums for R,G,B cones

   red_indices = find(cMosaic.pattern == 2);
   green_indices = find(cMosaic.pattern == 3);
   blue_indices = find(cMosaic.pattern == 4);
   avg_ratio_red = sum(sum(H(red_indices)))/numel(H(red_indices));
   avg_ratio_blue = sum(sum(H(blue_indices)))/numel(H(blue_indices));
   avg_ratio_green = sum(sum(H(green_indices)))/numel(H(green_indices));
   ratios = [avg_ratio_red avg_ratio_green avg_ratio_blue];
   figure();
   C = categorical(ratios,ratios,{'L','M','S'});
   h = histogram(C,'BarWidth',0.5);
   set(gca,'ylim',[0 1.2]);
   title('Ratio of absorption total absorptions by cone type');

% figure(); % subplot(3,1,1); % plot(1:length(red_indices),H(red_indices)); % subplot(3,1,2); % plot(1:length(green_indices),H(green_indices)); % subplot(3,1,3); % plot(1:length(blue_indices),H(blue_indices));