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