%%%%% Compare Martin's with Qian's model

%%%%% Martin: convolves phase-shifted flicker filters with image in left and right eye (simple cell)

%%%%% integrates over space and add squared input (complex cell)

%%%%% Qian : convolves phase-shifted motion filters with image in left and right eye (simple cell)

%%%%% integrate over space and add squared input (complex cell)

%%%%% by A Dolia and M Lages

%%%%% TO DO:

%%%%% make activation of flicker direction-selective

%%%%% smaller receptive fields tuned to left and right motion

%%%%% shift receptive fields according to stimulus phase shift of -pi/4...+pi/4

%%%%% adjust filter size and position

%%%%% 5 December 2003

%function [ComplexCell, maxv, maxi ] = StereoGaborWaveletV2(LeftImage, RightImage, Gabor_options)

% --------------- general settings ---------------------

clear all;

close all;

SpatialFrequency = 0.75; % was 1.0 cycle per unit

SpatialSize = 1.0; % SD of spatial envelope

SpatialStep = 0.1;

LeftLimit = -5.0;

RightLimit = 5.0;

TemporalFrequency = 0.75;

TemporalSize = 1.0; % SD of temporal envelope

TemporalStep = 0.1;

LowLimit = -5.0;

TopLimit = 5.0;

PhaseX = 0; % stimulus phase shift varied between -pi/2 and +pi/2

% Spatial and temporal extent

[filter_x_location,filter_y_location] = meshgrid(LowLimit : TemporalStep: TopLimit, LeftLimit: SpatialStep: RightLimit);

x = filter_x_location;

y = filter_y_location;

filter_bufL = (x).^2/(2*SpatialSize^2) + y.^2/(2*TemporalSize^2); % +0.25 equivalent to +pi/4??

filter_bufR = (x).^2/(2*SpatialSize^2) + y.^2/(2*TemporalSize^2); % -0.25 equivalent to -pi/4??

filter_envelopeL = exp(-filter_bufL);

filter_envelopeR = exp(-filter_bufR);

% Stimulus input to left and right eye

LeftImage = cos((2*pi*SpatialFrequency ) .* x - (2*pi*TemporalFrequency ) .* y + PhaseX/2);

RightImage = cos((2*pi*SpatialFrequency ) .* x - (2*pi*TemporalFrequency ) .* y - PhaseX/2);

%LeftImageSpectra = fft2(LeftImage);

%RightImageSpectra = fft2(RightImage);

%j = sqrt(-1);

ndepth = 24; % number of depth planes

delta_phase = pi/(ndepth-1); % pi/(ndepth-1) % extent in depth (rad)

for jj=1: ndepth

phase_diff = -0.5 * pi + (jj-1) * delta_phase; % -pi/2 to + pi/2

%PhaseL = -pi/2.0 + phase_diff;

%PhaseR = -pi/2.0 - phase_diff;

PhaseL = phase_diff/2;

PhaseR = -phase_diff/2;

%%%%%%%%%%%%%%% Martin's flicker filters (no motion)

%%%% positive velocity cos A + cos B = 2 cos( (A+B)/2 )cos( (A-B)/2)

MP_Left_cos_filter = 2 * filter_envelopeL .* cos((2*pi*SpatialFrequency ) .* x + PhaseL) .* ...

cos((2*pi*TemporalFrequency ) .* y );

MP_Right_cos_filter = -2 * filter_envelopeR .* sin((2*pi*SpatialFrequency ) .* x + PhaseR) .* ...

sin((2*pi*TemporalFrequency ) .* y );

%%%% positive velocity sin A + sin B = 2 sin( (A+B)/2 )cos( (A-B)/2)

MP_Left_sin_filter = 2 * filter_envelopeL .* sin((2*pi*SpatialFrequency ) .* x + PhaseL) .* ...

cos((2*pi*TemporalFrequency ) .* y );

MP_Right_sin_filter = 2 * filter_envelopeR .* cos((2*pi*SpatialFrequency ) .* x + PhaseR) .* ...

sin((2*pi*TemporalFrequency ) .* y );

MP_cos_simple_cell = zeros(201,1);

MP_sin_simple_cell = zeros(201,1);

% convolve image with filter over time then add over space

for i=1:101

MP_cos_simple_cell = MP_cos_simple_cell + conv(LeftImage(:,i),MP_Left_cos_filter(:,i)) + conv(RightImage(:,i),MP_Right_cos_filter(:,i));

MP_sin_simple_cell = MP_sin_simple_cell + conv(LeftImage(:,i),MP_Left_sin_filter(:,i)) + conv(RightImage(:,i),MP_Right_sin_filter(:,i));

end

% complex cell - squared input

MP_ComplexCell{jj} = MP_cos_simple_cell .* MP_cos_simple_cell + MP_sin_simple_cell .* MP_sin_simple_cell;

%%%%%%%%%%%%%%%%%%% Qian's motion filters (no flicker)

QP_Left_cos_filter = filter_envelopeL .* cos( 2*pi*SpatialFrequency .* x + PhaseL + ...

2*pi*TemporalFrequency .* y );

QP_Right_cos_filter= filter_envelopeR .* cos( 2*pi*SpatialFrequency .* x + PhaseR + ...

2*pi*TemporalFrequency .* y );

QP_Left_sin_filter = filter_envelopeL .* sin( 2*pi * SpatialFrequency .* x + PhaseL + ...

2*pi * TemporalFrequency .* y );

QP_Right_sin_filter= filter_envelopeR .* sin( 2*pi*SpatialFrequency .* x + PhaseR + ...

2*pi*TemporalFrequency .* y );

QP_cos_simple_cell = zeros(201,1);

QP_sin_simple_cell = zeros(201,1);

for i=1:101 % add over space

% convolve image with filter over time (columns)

QP_cos_simple_cell = QP_cos_simple_cell + conv(LeftImage(:,i),QP_Left_cos_filter(:,i)) + conv(RightImage(:,i),QP_Right_cos_filter(:,i));

QP_sin_simple_cell = QP_sin_simple_cell + conv(LeftImage(:,i),QP_Left_sin_filter(:,i)) + conv(RightImage(:,i),QP_Right_sin_filter(:,i));

end

% complex cell with squared input

QP_ComplexCell{jj} = (QP_cos_simple_cell) .* (QP_cos_simple_cell) + QP_sin_simple_cell .* QP_sin_simple_cell;

end % for jj

%%%%%% Determine max activity and depth plane for Martin's model

MP_maxv = MP_ComplexCell{1};

MP_minv = MP_ComplexCell{1};

MP_maxi = ones(201,1);

MP_maxv_depth = 1;

%ones(size(MP_cos_simple_cell,1));

for jj=2:ndepth

MP_max_value=max(MP_maxv);

for i_row=1 : size(MP_cos_simple_cell,1)

% for i_col=1 : size(LeftImage,2)

if MP_ComplexCell{jj}(i_row) > MP_maxv(i_row)

MP_maxv(i_row) = MP_ComplexCell{jj}(i_row); % sample max over depths at given time (i_row)

MP_maxi(i_row) = jj;

end % if

if MP_ComplexCell{jj}(i_row) > MP_max_value % find highest overall value

MP_max_value = MP_ComplexCell{jj}(i_row);

MP_maxv_depth = jj;

end % if

if MP_ComplexCell{jj}(i_row) < MP_minv(i_row)

MP_minv(i_row) = MP_ComplexCell{jj}(i_row);

end % if

% end % i_col

end % i_row

end; % jj

%%%%%% Determine max activity and depth plane for Qian's model

QP_maxv = QP_ComplexCell{1};

QP_minv = QP_ComplexCell{1};

QP_maxi = ones(201,1);

QP_maxv_depth = 1;

for jj=2:ndepth

QP_max_value=max(QP_maxv);

for i_row=1 : size(QP_cos_simple_cell,1)

% for i_col=1 : size(LeftImage,2)

if QP_ComplexCell{jj}(i_row) > QP_maxv(i_row)

QP_maxv(i_row) = QP_ComplexCell{jj}(i_row);

QP_maxi(i_row) = jj;

end % if

if QP_ComplexCell{jj}(i_row) > QP_max_value % find highest overall value

QP_max_value = QP_ComplexCell{jj}(i_row);

QP_maxv_depth = jj;

end % if

if QP_ComplexCell{jj}(i_row) < QP_minv(i_row)

QP_minv(i_row) = QP_ComplexCell{jj}(i_row);

end % if

% end % i_col

end % i_row

end; % jj

% QP_max_value = max(QP_maxv);

% QP_min_value = min(QP_minv);

%

% figure;

% subplot(2,2,1); imshow(mat2gray(Left_cos_filter ), [0 1]);

% subplot(2,2,2); imshow(mat2gray(Right_cos_filter), [0 1]);

% subplot(2,2,3); imshow(mat2gray(Left_sin_filter ), [0 1]);

% subplot(2,2,4); imshow(mat2gray(Right_sin_filter), [0 1]);

%figure;

%for jj=1:ndepth

% subplot(5,5, jj); imshow( mat2gray(ComplexCell{jj}/max_value), [0 1]);

%end

row = size(MP_ComplexCell{jj},1);

col = size(MP_ComplexCell{jj},2);

max_value = max([MP_max_value, QP_max_value]); %MN_max_value]) %; QP_max_value; QN_max_value])

% figure;

% for jj=1:ndepth

% subplot(5,5, jj); plot(MP_ComplexCell{jj}(:, round(col/2) ) ); axis([1 row 0 max_value]);

% end

% xlabel('space'); ylabel('time'); title('MP, Com. cell');

%

% figure;

% for jj=1:ndepth

% subplot(5,5, jj); plot(MN_ComplexCell{jj}(:, round(col/2) ) ); axis([1 row 0 max_value]);

% end

% xlabel('space'); ylabel('time'); title('MN, Com. cell');

%

% figure;

% for jj=1:ndepth

% subplot(5,5, jj); plot(QP_ComplexCell{jj}(:, round(col/2) ) ); axis([1 row 0 max_value]);

% end

% xlabel('space'); ylabel('time'); title('QP, Com. cell');

%

% figure;

% for jj=1:ndepth

% subplot(5,5, jj); plot(QN_ComplexCell{jj}(:, round(col/2) ) ); axis([1 row 0 max_value]);

% end

% xlabel('space'); ylabel('time'); title('QN, Com. cell');

%subplot(5,5, 25); imshow(mat2gray(maxv), [0 1]);

figure;

subplot(5,2,1) ; imshow(mat2gray(LeftImage ), [0 1]);

subplot(5,2,2) ; imshow(mat2gray(RightImage), [0 1]);

subplot(5,2,3); imshow(mat2gray(MP_ComplexCell{MP_maxv_depth}), [0 1]);

subplot(5,2,4); imshow(mat2gray(QP_ComplexCell{QP_maxv_depth}), [0 1]);

subplot(5,2,5) ; plot(MP_maxv); axis([1 row 0 MP_max_value]); xlabel('time'); ylabel('activation'); title('Martin, Complex Cells, MPos');

subplot(5,2,6) ; plot(MP_maxi); axis([1 row 1 ndepth ]); xlabel('time'); ylabel('depth (index)'); title('Martin, Index, MPos');

subplot(5,2,7) ; plot(QP_maxv); axis([1 row 0 QP_max_value]); xlabel('time'); ylabel('activation'); title('Qian, Com. Cells, MP');

subplot(5,2,8) ; plot(QP_maxi); axis([1 row 1 ndepth ]); xlabel('time'); ylabel('depth (index)'); title('Qian, Index, MP');

% subplot(5,2,7) ; plot(MN_maxv(round(col/2),:)); axis([1 row 0 max_value]); xlabel('time'); %ylabel('activation');title('Martin, Complex Cells, MNeg');

% subplot(5,2,8) ; plot(MN_maxi(round(col/2),:)); axis([1 row 1 ndepth ]); xlabel('time'); ylabel('depth %(index)');title('Martin, Index, MNeg');

%

% subplot(5,2,9) ; plot(QN_maxv(round(col/2),:)); axis([1 row 0 max_value]); xlabel('time'); ylabel('activation'); %title('Qian, Com. Cells, QN');

% subplot(5,2,10); plot(QN_maxi(round(col/2),:)); axis([1 row 1 ndepth ]); xlabel('time'); ylabel('depth %(index)'); title('Qian, Index, QN');