clear all
close all
rand('state',1')
randn('state',1')

% Getting vertical and horizontal disparity distribution

% This version just gives the start and end point; it does not do the smoothing. The idea is, I will then be able
% to play with the range, smoothing spaceconstant etc.


% Ok it's difficult to start from the position of objects in space
% So I will start from the retinae.

torsionlam = 2*pi/180; % for how HH torsion varies - as stated on p. 1987
torsionmu = 0.2;% for how HH torsion varies - as stated on p. 1987
f = 1; % radius of spherical retina (doesn't actually affect results but I like to keep it in so I can check dimensionality)
I = 3; % half interocular distance - in cm (NB this assumes a slightly low itnerocular of 6cm)
% As stated on p. 1989, ratio of I/f=3.

L = [I;0;0]; R = [-I;0;0]; % nodal points of the two eyes
xi  = 15 ; % we look at the central 15 degrees of visision
tanxi =tan(xi* pi/180);

% As stated in legend on Fig 5: "Helmholtz elevation V and mean
% Helmholtz azimuth Hc are picked from a gaussian with mean 0 degree and SD
% either 10 degrees (A) or 20 degrees (B)." This is where this is set:
HVSD = 20; % (so as in Fig 5B)

% As stated in legend to Fig 5: "To obtain vergence D, a random number
% is picked from a gaussian with mean 0 degree and SD either 3 degrees (A) or 8
% degrees (B), and then vergence is set to the absolute value of this number."
% This is where this is set:
DSD = 8; % (so as in Fig 5B)

% radius and aziretmuth of postions on the retina
% tan(angx) = x/f
% angx = [-1:0.1:1]*xi
% x = f*tan(angx)
% radius and aziretmuth of postions on the retina
th = [0.05:0.05:1]*2*pi;
r = f*tan([0:0.05:1]*xi*pi/180);
% As stated on p. 1989: "The starting points (x, y)
% were equally spaced in latitude and longitude on the hemispherical retina:
% x = f tan ? cos ?, y = f tan ? sin ?, where ? covered the range from 0 to 360
% degrees in steps of 18 degrees and ? ran from 0 to 15 degrees in steps of
% 0.75 degrees."



costh=cos(th);
sinth=sin(th);
N=1000;
nth=length(th);
nr=length(r);
dispx1=zeros(1,N*nr*nth*2);
dispx2=zeros(1,N*nr*nth*2);
dispy1=zeros(1,N*nr*nth*2);
dispy2=zeros(1,N*nr*nth*2);
cnt=0;
for j=1:N
    % Pick a new eye posture
    H = randn(1) * HVSD*pi/180;
    D = abs(randn(1) * DSD*pi/180);
    V = randn(1) * HVSD*pi/180;  
    
    HL = H - D/2;
    HR = H + D/2;
    TL = HHTorsion(HL,V,torsionlam + torsionmu*D);
    TR = HHTorsion(HR,V,-torsionlam - torsionmu*D);
    
    RotMatL = RotationMatrix3(V,[1 0 0]) * RotationMatrix3(HL,[0 1 0]) * RotationMatrix3(TL,[0 0 1]); 
    RotMatR = RotationMatrix3(V,[1 0 0]) * RotationMatrix3(HR,[0 1 0]) * RotationMatrix3(TR,[0 0 1]);
    
    %============================================================
    %% i) First way : starting in left retina
    
    % Epipole = point where interocular axis intersects right retina.
    % Independent of left eye position so we can work it out up front
    % using Eq 2.3 in the paper:
    q0 = f/tan(HR) * [-cos(TR) ; sin(TR) ]; 
    q0sq = dot(q0,q0); 
    q0minusasq = (q0sq-(f*tanxi)^2);
    
    % For speeding up computation of qinf:
    A = f* [ -cos(D)*cos(TR)*cos(TL)-sin(TR) *sin(TL) , cos(D)*cos(TR)*sin(TL) - sin(TR)*cos(TL) , -sin(D)*cos(TR); ...
            cos(D)*sin(TR)*cos(TL) - cos(TR) *sin(TL) , - cos(D)*sin(TR)*sin(TL) - cos(TR)*cos(TL), sin(D)*sin(TR)];
    B = [sin(D)*cos(TL); -sin(D)*sin(TL);-cos(D)];
    % (cf eq 2.4)
    
    for jr=1:length(r)
        for jt=1:length(th)
            % Left image position
            xL = r(jr) .* costh(jt);
            yL = r(jr) .* sinth(jt);            
            
            % check that this retinal position doesn't view back into the head:
            chk = RotMatL*[xL;yL;-f];
            if chk(3)>0
                plot(xL,yL,'ro');hold on
            else
                plot(xL,yL,'k.');hold on
                
                
                if HR~=0
                    
                    % Now find gradient and intercept of corresponding R epipolar line
                    % R epi line is yR = aR xR + bR
                    % qinf = f /( sin(D)*cos(TL)*xL - sin(D)*sin(TL)*yL-cos(D)*f) *  ...
                    %    [ -cos(D)*cos(TR)*cos(TL)*xL - sin(TR) *sin(TL) *xL + cos(D)*cos(TR)*sin(TL)*yL - sin(TR)*cos(TL)*yL - sin(D)*cos(TR)*f ; ...
                    %        cos(D)*sin(TR)*cos(TL)*xL - cos(TR) *sin(TL) *xL - cos(D)*sin(TR)*sin(TL)*yL - cos(TR)*cos(TL)*yL + sin(D)*sin(TR)*f ] ;
                    dott = dot(B,[xL;yL;f]);
                    qinf = A*[xL;yL;f] ./ dott;
                    % The epi line goes from q0 to qinf. Write rR = q0 + b(qinf-q0)
                    % Solve for b where it interesects the circle of central vision:
                    epidirection = qinf-q0;
                    epidirsq = dot(epidirection,epidirection);
                    q0dotepidir = dot(q0,epidirection);
                    discrim = q0dotepidir^2 - q0minusasq * epidirsq;
                    % if this is imaginary, the epipolar line does not fall within the circle of central vision; if it is zero, it just touches it.
                    
                    b1 = ( -q0dotepidir - sqrt(discrim) ) / epidirsq;
                    b2 = ( -q0dotepidir + sqrt(discrim) ) / epidirsq;
                    if discrim>0 && ( (sin(HR)/dott<0 && ((0<=b1 && b1<=1) || (0<=b2 && b2<=1) ) ) ||  (sin(HR)/dott>0 && b2>1 ) )
                        if sin(HR)/dott<0 && ((0<=b1 && b1<=1) || (0<=b2 && b2<=1) )
                            % epipolar line goes in between 0 and 1
                            if b1<=0  
                                % set b1=0;
                                xR1 = q0(1);
                                yR1 = q0(2);
                            else
                                % This is an expanded form which will be more accurate when q0 is infinite
                                xR1 = ( q0(1)*epidirection(2) - q0(2) * epidirection(1) ) * epidirection(2) / epidirsq - sqrt(discrim)/ epidirsq * epidirection(1);
                                yR1 = (q0(2)*epidirection(1)-q0(1) * epidirection(2))*epidirection(1) / epidirsq - sqrt(discrim)/ epidirsq * epidirection(2);
                            end
                            if b2>=1
                                % set b2=1;
                                xR2 = qinf(1);
                                yR2 = qinf(2);
                            else
                                % This is an expanded form which will be more accurate when q0 is infinite
                                xR2 = ( q0(1)*epidirection(2) - q0(2) * epidirection(1) ) * epidirection(2) / epidirsq + sqrt(discrim)/ epidirsq * epidirection(1);
                                yR2 = ( q0(2)*epidirection(1) - q0(1) * epidirection(2) ) * epidirection(1) / epidirsq + sqrt(discrim)/ epidirsq * epidirection(2);
                            end                        
                        end % sign HR
                        if sin(HR)/dott>0 & b2>1
                            % epipolar line goes from b=1 to b=infty
                            if b1<1  
                                % set b1=1;
                                xR1 = qinf(1);
                                yR1 = qinf(2);
                            else
                                % This is an expanded form which will be more accurate when q0 is infinite
                                xR1 = ( q0(1)*epidirection(2) - q0(2) * epidirection(1) ) * epidirection(2) / epidirsq - sqrt(discrim)/ epidirsq * epidirection(1);
                                yR1 = ( q0(2)*epidirection(1) - q0(1) * epidirection(2) ) * epidirection(1) / epidirsq - sqrt(discrim)/ epidirsq * epidirection(2);
                            end
                            % This is an expanded form which will be more accurate when q0 is infinite
                            xR2 = ( q0(1)*epidirection(2) - q0(2) * epidirection(1) ) * epidirection(2) / epidirsq + sqrt(discrim)/ epidirsq * epidirection(1);
                            yR2 = ( q0(2)*epidirection(1) - q0(1) * epidirection(2) ) * epidirection(1) / epidirsq + sqrt(discrim)/ epidirsq * epidirection(2);
                        end % sign HR            
                        
                        cnt=cnt+1;
                        dispx1(cnt) = 180/pi*(atan(xL)-atan(xR1));
                        dispy1(cnt) = 180/pi*(atan(yL)-atan(yR1));
                        dispx2(cnt) = 180/pi*(atan(xL)-atan(xR2));
                        dispy2(cnt) = 180/pi*(atan(yL)-atan(yR2));              
                        
                        %%%%%%%%%%%%%%%%%%%%5% check:
                        % This code isn't necessary but is just me checkign I haven't messed up the maths.
                        % Pick a random position on the epipolar line:
                        if sin(HR)/dott>0
                            beta=1+rand(1)*1e2;
                        elseif sin(HR)/dott<0
                            beta=rand(1);
                        end
                        xR = q0(1)+beta*epidirection(1); yR = q0(2)+beta*epidirection(2);
                        % (xR,yR) is now a random position on the epipolar line in the right eye.
                        % Here is where the two image points are in retinal coordinates
                        retR = [xR;yR;-f];                    retL = [xL;yL;-f];      
                        % We need to check they are consitent with a given object.
                        L=[I;0;0];            R=[-I;0;0];                               
                        mL = RotMatL * retL;            mR = RotMatR * retR;                spaceL = L + mL;            spaceR = R + mR;                
                        % spaceL, spaceR are now vectors describing where
                        % the image-points are in head-centered "space"
                        % coordinates.
                        invRotMatR = RotationMatrix3(-TR,[0 0 1]) * RotationMatrix3(-HR,[0 1 0]) * RotationMatrix3(-V,[1 0 0]);
                        fovL = L + RotMatL*[0;0;-f];                fovR = R + RotMatR*[0;0;-f];   
                        retnormL = RotMatL*[0;0;1];            retnormR = RotMatR*[0;0;1];      
                        % Work out where the object corresponding to these
                        % image points is. It is at the intersection of
                        % these two lines: (a) the line through L, the
                        % centre of projection of the left eye, and spaceL,
                        % the  position of the left image in headcentred
                        % coordinates, and (b) ditto for the right eye. 
                        % THat is, it is at the solution of the vector equation
                        % object = L - retlam * mL = R - retgam * mR
                        % where the parameters retlam and retgam have to be
                        % solved for. You can work out that:
                        % retlam = ( dot(L-R,mL) * dot(mR,mR) - dot(L-R,mR) * dot(mL,mR) ) / ( dot(mL,mL) * dot(mR,mR) - (dot(mL,mR))^2);
                        % or equivalently
                        retlam = beta/(beta-1)*2*I*sin(HR)/dott
                        % so
                        object = L - retlam * mL; 
                        % Now project this object back onto the two retinas:
                        [tmp,chkretimageL] = ImagePos(object,HL,V,TL,L,f);   
                        [tmp,chkretimageR] = ImagePos(object,HR,V,TR,R,f);    
                        % and check it agrees with our original starting
                        % point
                        % Stop running if we encounter an error: if the
                        % object is behind the eyes, or the relative and absolute error
                        % in image position both exceed 1e-4:
                        if retmu<0 |retlam<0 | object(3)<0 | (abs((chkretimageL(1) - xL)/xL)>1e-4 & abs(chkretimageL(1)-xL)>1e-4) | (abs((chkretimageL(2) - yL)/yL)>1e-4 & abs(chkretimageL(2)-yL)>1e-4) | (abs((chkretimageR(1) - xR)/xR)>1e-4 & abs(chkretimageR(1)-xR)>1e-4) | (abs((chkretimageR(2) - yR)/yR)>1e-4 & abs(chkretimageR(2)-yR)>1e-4)
                            retmu
                            retlam
                            object
                            chkretimageL'
                            [xL yL]
                            chkretimageR'
                            [xR yR]                    
                            stop
                        end
                        %%%%%%%%%%%%%%%%%%%%5% end check:                                                
                    end % checking whether epipolar line falls within central vision
                else %if HR==0
                    % Now find gradient and intercept of corresponding R epipolar line
                    % R epi line is yR = aR xR + bR
                    % qinf = f /( sin(D)*cos(TL)*xL - sin(D)*sin(TL)*yL-cos(D)*f) *  ...
                    %    [ -cos(D)*cos(TR)*cos(TL)*xL - sin(TR) *sin(TL) *xL + cos(D)*cos(TR)*sin(TL)*yL - sin(TR)*cos(TL)*yL - sin(D)*cos(TR)*f ; ...
                    %        cos(D)*sin(TR)*cos(TL)*xL - cos(TR) *sin(TL) *xL - cos(D)*sin(TR)*sin(TL)*yL - cos(TR)*cos(TL)*yL + sin(D)*sin(TR)*f ] ;
                    dott = dot(B,[xL;yL;f]);
                    qinf = A*[xL;yL;f] ./ dott;
                    % The epi line goes from q0 to qinf. Write rR = q0 + b(qinf-q0)
                    % Solve for b where it interesects the circle of central vision:
                    epidirection = f*2*I*[ cos(HR)*cos(TR)  ;  -sin(TR)] / dott;
                    epidirsq = dot(epidirection,epidirection);
                    qinfdotepidir = dot(qinf,epidirection);
                    qinfminusasq = (dot(qinf,qinf)-(f*tanxi)^2);
                    discrim = qinfdotepidir^2 - qinfminusasq * epidirsq;
                    % if this is imaginary, the epipolar line does not fall within the circle of central vision; if it is zero, it just touches it.
                    b1 = ( -qinfdotepidir - sqrt(discrim) ) / epidirsq;
                    b2 = ( -qinfdotepidir + sqrt(discrim) ) / epidirsq;
                    if discrim>0 
                        % epipolar line goes in between beta=0 and beta=infty
                        if b1>0
                            xR1 = qinf(1) + b1*epidirection(1);
                            yR1 = qinf(2) + b1*epidirection(2);
                        else
                            xR1=qinf(1);
                            yR1=qinf(2);
                        end
                        if b2>0
                            xR2 = qinf(1) + b2*epidirection(1);
                            yR2 = qinf(2) + b2*epidirection(2);
                        else
                            xR2=qinf(1);
                            yR2=qinf(2);
                        end        
                        
                        cnt=cnt+1;
                        dispx1(cnt) = 180/pi*(atan(xL)-atan(xR1));
                        dispy1(cnt) = 180/pi*(atan(yL)-atan(yR1));
                        dispx2(cnt) = 180/pi*(atan(xL)-atan(xR2));
                        dispy2(cnt) = 180/pi*(atan(yL)-atan(yR2));              
                        
                        %%%%%%%%%%%%%%%%%%%%5% check:
                        L=[I;0;0];            R=[-I;0;0];                               
                        invRotMatR = RotationMatrix3(-TR,[0 0 1]) * RotationMatrix3(-HR,[0 1 0]) * RotationMatrix3(-V,[1 0 0]);
                        fovL = L + RotMatL*[0;0;-f];                fovR = R + RotMatR*[0;0;-f];                       
                        if sin(HR)/dott>0
                            beta=1+rand(1)*1e2;
                        elseif sin(HR)/dott<0
                            beta=rand(1);
                        end
                        xR = qinf(1)+beta*epidirection(1); yR = qinf(2)+beta*epidirection(2);                   retR = [xR;yR;-f];                    retL = [xL;yL;-f];      
                        mL = RotMatL * retL;            mR = RotMatR * retR;                spaceL = L + mL;            spaceR = R + mR;                
                        retnormL = RotMatL*[0;0;1];            retnormR = RotMatR*[0;0;1];                 
                        retlam = beta/(beta-1)*2*I*sin(HR)/dott;
                        %%                retmu = f /  (2*sin(HR)*I - sin(HR)*retlam*cos(HL)*cos(TL)*xL + sin(HR)*retlam*cos(HL)*sin(TL)*yL + sin(HR)*retlam*sin(HL)*f + cos(HR)*retlam*xL*sin(HL)*cos(TL) - cos(HR)*retlam*yL*sin(HL)*sin(TL) + cos(HR)*retlam*cos(HL)*f);
                        retmu = f / (2*sin(HR)*I - retlam*dott);
                        object = L - retlam * mL;            [tmp,chkretimageL] = ImagePos(object,HL,V,TL,L,f);                            
                        chkretimageR = -retmu * invRotMatR *(object-R);
                        if retmu<0 |retlam<0 | object(3)<0 | (abs((chkretimageL(1) - xL)/xL)>1e-4 & abs(chkretimageL(1)-xL)>1e-4) | (abs((chkretimageL(2) - yL)/yL)>1e-4 & abs(chkretimageL(2)-yL)>1e-4) | (abs((chkretimageR(1) - xR)/xR)>1e-4 & abs(chkretimageR(1)-xR)>1e-4) | (abs((chkretimageR(2) - yR)/yR)>1e-4 & abs(chkretimageR(2)-yR)>1e-4)
                            retmu
                            retlam
                            object
                            chkretimageL'
                            [xL yL]
                            chkretimageR'
                            [xR yR]                    
                            stop
                        end
                        %%%%%%%%%%%%%%%%%%%%5% end check:
                    end % checking whether epipolar line falls within central vision
                end % checking whether HR=0
            end% checking whether retinal position is acceptable
        end % next th
    end % next r
    
    %============================================================
    %% ii) Second way : starting in right retina      
    q0 = f/tan(HL) * [-cos(TL) ; sin(TL) ];
    q0sq = dot(q0,q0);
    q0minusasq = (q0sq-(f*tanxi)^2);
    
    % For speeding up computation of qinf:
    A = f* [ -cos(D)*cos(TL)*cos(TR)-sin(TL) *sin(TR) , cos(D)*cos(TL)*sin(TR) - sin(TL)*cos(TR) ,sin(D)*cos(TL); ...
            cos(D)*sin(TL)*cos(TR) - cos(TL) *sin(TR) , - cos(D)*sin(TL)*sin(TR) - cos(TL)*cos(TR),-sin(D)*sin(TL)];
    B = [-sin(D)*cos(TR); + sin(D)*sin(TR);-cos(D)];

    for jr=1:length(r)
        for jt=1:length(th)
            % Reft image position
            xR = r(jr) .* costh(jt);
            yR = r(jr) .* sinth(jt);
                        
            % check that this retinal position doesn't view back into the head:
            chk = RotMatR*[xR;yR;-f];
            if chk(3)>0
                plot(xR,yR,'ro');hold on
            else
                plot(xR,yR,'k.');hold on
                if HL~=0
                    % Now find gradient and intercept of corresponding L epipolar line
                    % L epi line is yL = aL xL + bL
                    %qinf = f /( sin(D)*cos(TR)*xR + sin(D)*sin(TR)*yR-cos(D)*f) *  ...
                    %    [ -cos(D)*cos(TL)*cos(TR)*xR - sin(TL) *sin(TR) *xR + cos(D)*cos(TL)*sin(TR)*yR - sin(TL)*cos(TR)*yR - sin(D)*cos(TL)*f ; ...
                    %        cos(D)*sin(TL)*cos(TR)*xR - cos(TL) *sin(TR) *xR - cos(D)*sin(TL)*sin(TR)*yR - cos(TL)*cos(TR)*yR + sin(D)*sin(TL)*f ] ;    
                    % but it is quicker to write
                    dott = dot(B,[xR;yR;f]);
                    qinf = A*[xR;yR;f] ./ dott;
                    % The epi line goes from q0 to qinf. Write rL = q0 + retmu(qinf-q0): retmu goes from 0 to 1
                    % Solve for retmu where it interesects the circle of central vision:
                    epidirection = qinf-q0;
                    epidirsq = dot(epidirection,epidirection);
                    q0dotepidir = dot(q0,epidirection);
                    discrim = q0dotepidir^2 - q0minusasq * epidirsq;
                    % if this is negative, the epipolar line does not fall within the circle of central vision; if it is zero, it just touches it.
                    
                    b1 = ( -q0dotepidir - sqrt(discrim) ) / epidirsq;
                    b2 = ( -q0dotepidir + sqrt(discrim) ) / epidirsq;
                    % need truncated line to fall witihin central vision:
                    if discrim>0 & ( (sin(HL)/dott<0 & b2>=1 ) | ( sin(HL)<0 & ( (0<=b1 & b1<=1) | (0<=b2 & b2<=1) )  ))
                        %
                        if sin(HL)/dott<0 & b2>=1
                            % epipolar line goes from 1 to infinity
                            if b1<1  
                                xL1 = qinf(1);
                                yL1 = qinf(2);
                            else
                                % This is an expanded form which will be more accurate when q0 is infinite
                                xL1 = ( q0(1)*epidirection(2) - q0(2) * epidirection(1) ) * epidirection(2) / epidirsq - sqrt(discrim)/ epidirsq * epidirection(1);
                                yL1 = (q0(2)*epidirection(1)-q0(1) * epidirection(2))*epidirection(1) / epidirsq - sqrt(discrim)/ epidirsq * epidirection(2);
                            end
                            % This is an expanded form which will be more accurate when q0 is infinite
                            xL2 = ( q0(1)*epidirection(2) - q0(2) * epidirection(1) ) * epidirection(2) / epidirsq + sqrt(discrim)/ epidirsq * epidirection(1);
                            yL2 = ( q0(2)*epidirection(1) - q0(1) * epidirection(2) ) * epidirection(1) / epidirsq + sqrt(discrim)/ epidirsq * epidirection(2);
                        end % sign of HL
                        if sin(HL)/dott>0 & ( (0<=b1 & b1<=1) | (0<=b2 & b2<=1) )  
                            % epipolar line goes in between 0 and 1
                            if b1<0  
                                % set b1=0;
                                xL1 = q0(1);
                                yL1 = q0(2);
                            else
                                % This is an expanded form which will be more accurate when q0 is infinite
                                xL1 = ( q0(1)*epidirection(2) - q0(2) * epidirection(1) ) * epidirection(2) / epidirsq - sqrt(discrim)/ epidirsq * epidirection(1);
                                yL1 = (q0(2)*epidirection(1)-q0(1) * epidirection(2))*epidirection(1) / epidirsq - sqrt(discrim)/ epidirsq * epidirection(2);
                            end
                            if b2>1
                                % set retmu2=1;
                                xL2 = qinf(1);
                                yL2 = qinf(2);
                            else
                                % This is an expanded form which will be more accurate when q0 is infinite
                                xL2 = ( q0(1)*epidirection(2) - q0(2) * epidirection(1) ) * epidirection(2) / epidirsq + sqrt(discrim)/ epidirsq * epidirection(1);
                                yL2 = ( q0(2)*epidirection(1) - q0(1) * epidirection(2) ) * epidirection(1) / epidirsq + sqrt(discrim)/ epidirsq * epidirection(2);
                            end  
                        end % sign of HL
                        
                        cnt=cnt+1;
                        dispx1(cnt) = 180/pi*(atan(xL1)-atan(xR));
                        dispy1(cnt) = 180/pi*(atan(yL1)-atan(yR));
                        dispx2(cnt) = 180/pi*(atan(xL2)-atan(xR));
                        dispy2(cnt) = 180/pi*(atan(yL2)-atan(yR));              
                        
                        %%%%%%%%%%%%%%%%%%%%5% check:
                        L=[I;0;0];            R=[-I;0;0];       
                        invRotMatL = RotationMatrix3(-TL,[0 0 1]) * RotationMatrix3(-HL,[0 1 0]) * RotationMatrix3(-V,[1 0 0]);
                        fovL = L + RotMatL*[0;0;-f];                fovR = R + RotMatR*[0;0;-f];            
                        if sin(HL)/dott<0
                            beta=1+rand(1)*1e2;
                        elseif sin(HL)/dott>0
                            beta=rand(1);
                        end
                        xL = q0(1)+beta*epidirection(1); yL = q0(2)+beta*epidirection(2);                      retL = [xL;yL;-f];      retR = [xR;yR;-f];                            
                        mL = RotMatL * retL;            mR = RotMatR * retR;                spaceL = L + mL;            spaceR = R + mR;                
                        retnormL = RotMatL*[0;0;1];            retnormR = RotMatR*[0;0;1];                 
                        retlam = -beta/(beta-1)*2*I*sin(HL)/dott;
                        retmu = f / (-2*sin(HL)*I - retlam*dott);
                        object = R - retlam * mR;            [tmp,chkretimageR] = ImagePos(object,HR,V,TR,R,f);            
                        chkretimageL = - retmu * invRotMatL*(object-L);
                        if retmu<0 |retlam<0 | object(3)<0 | (abs((chkretimageL(1) - xL)/xL)>1e-4 & abs(chkretimageL(1)-xL)>1e-4) | (abs((chkretimageL(2) - yL)/yL)>1e-4 & abs(chkretimageL(2)-yL)>1e-4) | (abs((chkretimageR(1) - xR)/xR)>1e-4 & abs(chkretimageR(1)-xR)>1e-4) | (abs((chkretimageR(2) - yR)/yR)>1e-4 & abs(chkretimageR(2)-yR)>1e-4)
                            retmu
                            retlam
                            object
                            chkretimageL'
                            [xL yL]
                            chkretimageR'
                            [xR yR]                    
                            stop
                        end
                        %%%%%%%%%%%%%%%%%%%%5% end check:
                        
                        
                    end % checking whether epipolar line falls within central vision
                    
                    
                else %if HL==0
                    % Now find gradient and intercept of corresponding R epipolar line
                    dott = dot(B,[xR;yR;f]);
                    qinf = A*[xR;yR;f] ./ dott;
                    % The epi line goes from q0 to qinf. Write rR = q0 + b(qinf-q0)
                    % Solve for b where it interesects the circle of central vision:
                    epidirection = -f*2*I*[ cos(HL)*cos(TL)  ;  -sin(TL)] / dott;
                    epidirsq = dot(epidirection,epidirection);
                    qinfdotepidir = dot(qinf,epidirection);
                    qinfminusasq = (dot(qinf,qinf)-(f*tanxi)^2);
                    discrim = qinfdotepidir^2 - qinfminusasq * epidirsq;
                    % if this is imaginary, the epipolar line does not fall within the circle of central vision; if it is zero, it just touches it.
                    b1 = ( -qinfdotepidir - sqrt(discrim) ) / epidirsq;
                    b2 = ( -qinfdotepidir + sqrt(discrim) ) / epidirsq;
                    if discrim>0 
                        % epipolar line goes in between beta=0 and beta=infty
                        if b1>0
                            xL1 = qinf(1) + b1*epidirection(1);
                            yL1 = qinf(2) + b1*epidirection(2);
                        else
                            xL1=qinf(1);
                            yL1=qinf(2);
                        end
                        if b2>0
                            xL2 = qinf(1) + b2*epidirection(1);
                            yL2 = qinf(2) + b2*epidirection(2);
                        else
                            xL2=qinf(1);
                            yL2=qinf(2);
                        end        
                        
                        cnt=cnt+1;
                        dispx1(cnt) = 180/pi*(atan(xL1)-atan(xR));
                        dispy1(cnt) = 180/pi*(atan(yL1)-atan(yR));
                        dispx2(cnt) = 180/pi*(atan(xL2)-atan(xR));
                        dispy2(cnt) = 180/pi*(atan(yL2)-atan(yR));              
                        
                        %%%%%%%%%%%%%%%%%%%%5% check:
                        L=[I;0;0];            R=[-I;0;0];                               
                        invRotMatR = RotationMatrix3(-TR,[0 0 1]) * RotationMatrix3(-HR,[0 1 0]) * RotationMatrix3(-V,[1 0 0]);
                        fovL = L + RotMatL*[0;0;-f];                fovR = R + RotMatR*[0;0;-f];                       
                        if sin(HR)/dott>0
                            beta=1+rand(1)*1e2;
                        elseif sin(HR)/dott<0
                            beta=rand(1);
                        end
                        xL = qinf(1)+beta*epidirection(1); yL = qinf(2)+beta*epidirection(2);                   retR = [xR;yR;-f];                    retL = [xL;yL;-f];      
                        mL = RotMatL * retL;            mR = RotMatR * retR;                spaceL = L + mL;            spaceR = R + mR;                
                        retnormL = RotMatL*[0;0;1];            retnormR = RotMatR*[0;0;1];                 
                        retlam = beta/(beta-1)*2*I*sin(HR)/dott;
                        %%                retmu = f /  (2*sin(HR)*I - sin(HR)*retlam*cos(HL)*cos(TL)*xL + sin(HR)*retlam*cos(HL)*sin(TL)*yL + sin(HR)*retlam*sin(HL)*f + cos(HR)*retlam*xL*sin(HL)*cos(TL) - cos(HR)*retlam*yL*sin(HL)*sin(TL) + cos(HR)*retlam*cos(HL)*f);
                        retmu = f / (2*sin(HR)*I - retlam*dott);
                        object = L - retlam * mL;            [tmp,chkretimageL] = ImagePos(object,HL,V,TL,L,f);                            
                        chkretimageR = -retmu * invRotMatR *(object-R);
                        if retmu<0 |retlam<0 | object(3)<0 | (abs((chkretimageL(1) - xL)/xL)>1e-4 & abs(chkretimageL(1)-xL)>1e-4) | (abs((chkretimageL(2) - yL)/yL)>1e-4 & abs(chkretimageL(2)-yL)>1e-4) | (abs((chkretimageR(1) - xR)/xR)>1e-4 & abs(chkretimageR(1)-xR)>1e-4) | (abs((chkretimageR(2) - yR)/yR)>1e-4 & abs(chkretimageR(2)-yR)>1e-4)
                            retmu
                            retlam
                            object
                            chkretimageL'
                            [xL yL]
                            chkretimageR'
                            [xR yR]                    
                            stop
                        end
                        %%%%%%%%%%%%%%%%%%%%5% end check:
                    end % checking whether epipolar line falls within central vision
                end % checking whether HL=0
                
            end% checking whether retinal position is acceptable
            
        end % next th
    end % next r
    
    j
end %next fixation point
dispx1 = dispx1(1:cnt);
dispx2 = dispx2(1:cnt);
dispy1 = dispy1(1:cnt);
dispy2 = dispy2(1:cnt);
save(sprintf('DisparityDistrib_skel_rth_xi15_HVSD%d_DSD%d',HVSD,DSD));


