% 9th Janb 06. This differs from FitPentagon.m in that it evaluates the integral in real space for high freqs (when the physical RF is small) and in Fourier space for 
% low freqs (when the Fourier-space RF is small). For each SF, it works out which is going to be faster and uses that method.

clear all
close all
imL = double(flipud(imread('left2.png'))) - 127;
imR = double(flipud(imread('right2.png'))) - 127;

[ny,nx]=size(imL);
x= ([1:nx]/nx-0.5)*10;
y= ([1:ny]/ny-0.5)*10*ny/nx;

FT=1;
PHYS=2;
SFBW = 1.5;
SFs = 0.5;%[4 2 1 0.5];%[0.25 0.5 1 2 ];
RFphase = 0;
oris = [30:30:180]*pi/180;
xprefs = x(round([nx/8  : 7*nx/8]));
yprefs = y(round([nx/8  : 7*nx/8]));;

SDs = sqrt(log(2))./(2*pi*SFs).*(2^SFBW+1)/(2^SFBW-1);
[x2,y2]=meshgrid(x,y);
immx=max([max(max(imL)) max(max(imR))]);
immn=min([min(min(imL)) min(min(imR))]);

% Take FTs of images
[FTimL,freqx,freqy]=jfft2(x,y,imL);
[FTimR,freqx,freqy]=jfft2(x,y,imR);
[freqx2,freqy2]=meshgrid(freqx,freqy);

dfx = min(diff(freqx))
dfy = min(diff(freqy))
dx = min(diff(x))
dy = min(diff(y))

imfig=figure('pos',[3         253        1275         447]);
subplot(1,2,1)
imagesc(x,y,imL)
caxis([immn immx])
subplot(1,2,2)
imagesc(x,y,imR)
caxis([immn immx])
set(gca,'ydir','normal')
for j=1:2
    subplot(1,2,j)
    hold on
    set(gca,'ydir','normal')
    axis equal tight
    colormap gray
    % Mark on places where disparity is evaluated
    for jx=1:length(xprefs)
        plot(xprefs(jx)*ones(size(yprefs)),yprefs,'m.')
    end
    % Mark on biggest and smallest channel
    if j==1
        chanSF = min(SFs);
    else
        chanSF = max(SFs);
    end
    chanSD = sqrt(log(2))/(2*pi*chanSF)*(2^SFBW+1)/(2^SFBW-1);
    RF = exp( - 0.5*(x2.^2+y2.^2)/chanSD^2) .* cos(2*pi*chanSF.*x2 - RFphase);
    contour(x,y,RF,[0.2:0.2:1]*max(max(abs(RF))),'y');   
    contour(x,y,RF,-[0.2:0.2:1]*max(max(abs(RF))),'y--');   
end

chk = figure;
res = figure('pos',[ 152         338        1061         351]);
Dx_phase = zeros(length(yprefs),length(xprefs),length(SFs),length(oris));
C_phase = Dx_phase;
M_phase = Dx_phase;
minphase = Dx_phase;
C_maxcorr = Dx_phase;
Dx_maxcorr = Dx_phase;
M_maxcorr = Dx_phase;
Dx_maxE = Dx_phase;
methused = Dx_phase;
noris = length(oris);
nSFs = length(SFs);
method=PHYS;0
for jSF=1:nSFs
    chanSF = SFs(jSF);
    chanSD = sqrt(log(2))/(2*pi*chanSF)*(2^SFBW+1)/(2^SFBW-1);
    halfSDsq=0.5/chanSD^2;
    Dxprefs = [-0.2 : min([0.1 1/SFs(jSF)/5]) : 0.2];
    
    for jy=1:length(yprefs)
        ypref = yprefs(jy);
        for jx=1:length(xprefs)
            xpref = xprefs(jx);
            
            % For physical-space method:
            % To save time, we are going to begin by taking out a chunk of the image and just looking at that.
            jj = find( abs(x2-xpref) < 5*chanSD+max(abs(Dxprefs))/2 & abs(y2-ypref) < 5*chanSD);
            x2red = x2(jj);
            y2red = y2(jj);
            imLred = imL(jj);
            imRred = imR(jj);
            
            
            % For checking:%%%%%%%%%%%%%%%%%%%%%%%%
            xx = unique(x2red);
            yy = unique(y2red);
            figure(chk)
            subplot(2,2,1)
            hold off
            imagesc(x,y,imL);
            subplot(2,2,2)
            hold off
            imagesc(x,y,imR);
            for k=1:2
                subplot(2,2,k)
                hold on
                plot([xx(1) xx(end)],yy(1)*[1 1],'m-')
                plot([xx(1) xx(end)],yy(end)*[1 1],'m-')
                plot(xx(1)*[1 1],[yy(1) yy(end)],'m-')
                plot(xx(end)*[1 1],[yy(1) yy(end)],'m-')
            end
            subplot(2,2,3)
            hold off
            imagesc(xx,yy,reshape(imLred,length(yy),length(xx)))
            caxis([immn immx])
            hold on
            plot(xprefs(jx),yprefs(jy),'m.')
            subplot(2,2,4)
            hold off
            imagesc(xx,yy,reshape(imLred,length(yy),length(xx)))
            caxis([immn immx])
            hold on
            RF = exp( - 0.5*((x2red-xpref).^2+(y2red-ypref).^2)/chanSD^2) .* cos(2*pi*chanSF.*(x2red-xpref) - RFphase);
            RF = reshape(RF,length(yy),length(xx));
            contour(xx,yy,RF,[0.2:0.2:1]*max(max(abs(RF))),'y');   
            contour(xx,yy,RF,-[0.2:0.2:1]*max(max(abs(RF))),'y--');    
            colormap gray
            drawnow
            %%%%%%%%%%%%%%%%%%%%%%%%%
            for jori=1:noris
                % For the 1st orientation in each SF, look at doing the integrals in Fourier space
                % If this 1st orientation looked as if it would be faster to do the integrals in Fourier space, 
                % then keep on in Fourier space for the rest of this SF channel:
                if jori==1 | method == FT
                    % For Fourier-space method
                    x2rot = cos(oris(jori)).*x2 + sin(oris(jori)).*y2;
                    RFE = exp( - halfSDsq*(x2.^2+y2.^2) ) .* cos(2*pi*chanSF.*x2rot - RFphase);
                    RFO = exp( - halfSDsq*(x2.^2+y2.^2) ) .* sin(2*pi*chanSF.*x2rot - RFphase);
                    % RF = RFE.*cos(Dphipref/2) + RFO.*sin(Dphipref/2)
                    [FTRFE,freqx,freqy] = jfft2(x,y,RFE);
                    [FTRFO,freqx,freqy] = jfft2(x,y,RFO);
                    jFT = find(abs(FTRFE)>3.7e-6.*max2d(abs(FTRFE)) | abs(FTRFO)>3.7e-6.*max2d(abs(FTRFO)));
                    freqx=freqx2(jFT);
                    freqy=freqy2(jFT);
                    FTimLred = FTimL(jFT);
                    FTimRred = FTimR(jFT);
                    FTRFE = FTRFE(jFT);
                    FTRFO = FTRFO(jFT);
                    if length(jFT)<length(jj)
                        method=FT;
                    else
                        method=PHYS;
                    end
                    
                end
                if method==FT
                    % then do the integral in Fourier space
                    [Dx_phase(jy,jx,jSF,jori),C_phase(jy,jx,jSF,jori),M_phase(jy,jx,jSF,jori),minphase(jy,jx,jSF,jori),...
                            Dx_maxcorr(jy,jx,jSF,jori),C_maxcorr(jy,jx,jSF,jori),M_maxcorr(jy,jx,jSF,jori),...
                            Dx_maxE(jy,jx,jSF,jori)] ...
                        = FindDisparityFT(i*2*pi*freqx,i*2*pi*freqy,FTimLred,FTimRred,FTRFE,FTRFO,xpref,ypref,Dxprefs);  
                    methused(jy,jx,jSF,jori) = FT;
                else
                    % do thte integral in physical space
                    [Dx_phase(jy,jx,jSF,jori),C_phase(jy,jx,jSF,jori),M_phase(jy,jx,jSF,jori),minphase(jy,jx,jSF,jori),...
                            Dx_maxcorr(jy,jx,jSF,jori),C_maxcorr(jy,jx,jSF,jori),M_maxcorr(jy,jx,jSF,jori),...
                            Dx_maxE(jy,jx,jSF,jori)] ...
                        = FindDisparity(x2red,y2red,imLred,imRred,chanSF,chanSD,oris(jori),RFphase,xpref,ypref,Dxprefs);    
                    methused(jy,jx,jSF,jori) = PHYS;
                end
            end
            %%%%%%%%%%%%%%%%%%%%%%%
            figure(res);
            for jori=1:noris
                subplot(nSFs,noris,jori+(jSF-1)*noris)
                jimagesc(xprefs,yprefs,Dx_phase(:,:,jSF,jori),colormap(jet),min(Dxprefs),max(Dxprefs))
                set(gca,'ydir','norm')
                title(sprintf('%4.1cpd, %d^o',SFs(jSF),round(180/pi*oris(jori))))
            end
            drawnow
            %%%%%%%%%%%%%%%%%%%%%%%%%
            
            
            
        end
        
        save FitPentagon_SFBW1.5fineSF16.mat
    end
end

