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;

% oisH.timeStep

%% 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));