
clear all
close all

% Getting vertical and horizontal disparity distribution

% This reads in a skeletn file and smoothes it.
skeletonfile = 'DisparityDistrib_skel_rth_xi15_HVSD20_DSD8.mat';
load(skeletonfile)
fleshfile = strrep(skeletonfile,'skel','flesh');

% Set SD of smoothing:
SD=0.0005;
% Set up disparity "bins". Notice I have chosen a wider range for dx than dy.
dx=[-1:SD:1];
dy=[-0.05:SD:0.05];
% Write these chosen limits into the filename:
fleshfile = strrep(fleshfile,'.mat','_xlim1_ylim0.05.mat');

[dx2,dy2] = meshgrid(dx,dy);

prob = zeros(size(dx2));
recip2SDsq = 1/2/SD^2;
recipSDsqrt2 = 1/SD/sqrt(2);

length(dispx1)
halfdisp = 0.5*sqrt((dispx1-dispx2).^2 + (dispy1-dispy2).^2);
a = atan2([dispy1-dispy2],[dispx1-dispx2]);
cosa = cos(a); sina = sin(a);
dispmidx = 0.5*(dispx1+dispx2);
dispmidy = 0.5*(dispy1+dispy2);
figure
jj=fix(rand(1,length(dispx1))*length(dispx1));
vol=0;

for j=1:length(dispx1)
    j
    dxp = (dx2-dispmidx(j))*cosa(j) + (dy2-dispmidy(j))*sina(j);
    dyp = -(dx2-dispmidx(j))*sina(j) + (dy2-dispmidy(j))*cosa(j);
    prob = prob +  0.5/sqrt(2*pi)/SD * exp(-dyp.^2.*recip2SDsq) .* ( erf((halfdisp(j)-dxp)*recipSDsqrt2) + erf((halfdisp(j)+dxp)*recipSDsqrt2) );
    % calculatse contrib to total volume under this "rod" of disparity
    vol = vol + 2*halfdisp(j) ;
    if fix(j/500)==j/500
        imagesc(dx,dy,prob);
        drawnow
        % What percentage of the distribution is displayed?
        sum(sum(prob))*min(diff(dx))*min(diff(dy))/vol
		save DisparityDistrib_flesh_rth_xi15_HVSD10_DSD3_xlim1_ylim0.05.mat
    end
end %next fixation point
imagesc(dx,dy,prob);
hold on
mn = 1/sqrt(2*pi)/SD *exp(-1);
contour(dx,dy,prob,[0 mn],'w');
z=prctile2dPDF(prob,min(diff(dx))*min(diff(dy)),vol,90)
contour(dx,dy,prob,[0 z],'w');
sum(sum(prob(prob>z)))*min(diff(dx))*min(diff(dy))/vol

ylabel('vertical disparity  \deltay','fontsize',14)
xlabel('horizontal disparity \deltax','fontsize',14)
