function [pixel_radius,pixel_center] = findkcal(calibdata, threepoints)
%SAXS/FINDKCAL Compute q/pixel calibration for SAXS image.
%
%   [pixel_radius,pixel_center] = FINDKCAL(S, THREEPOINTS) returns the horizontal and vertical 
%   radius (pixel_radius(1) and pixel_radius(2)) of a SAXS image, as well as the center
%   (pixel_center(1) and pixel_center(2))using data from a reference file 
%   (for example taken with Silver Behenate)and 3 approximate points 
%   (THREEPOINTS) on the partial 2D calibration circle.
%
%   pixel_center and pixel_radius are two-element vector returning 
%   the horizontal and vertical radius of the calibration circle and the


% First transfer to non-centered data 
threepoints(:,1)=threepoints(:,1);
threepoints(:,2)=threepoints(:,2);

% First get the approximate circle parameters using the 3 points
[xcen,ycen,expecrad]=solvecircle(threepoints);

s=floor(expecrad/10); % this insures that the center determination is not scewed by the projection of a ring, 
% in that it only takes a projection over 1/10 of a
% radian or approximately 0.5%. this will skew the mass average by less
% than 0.1%
s1=1/8;             % this values is the how much "extra" to include around the expected peak. I
% if the expected detector distance is not accurate
% enough or one is using a very high order reflection (with close
% peaks) this number may have to be adjusted.

try
    xpeakplusdata=sum(calibdata.pixels(round(ycen-s):round(ycen+s),floor(xcen+(1-s1)*expecrad):floor(xcen+(1+s1)*expecrad)));
    xpeakminusdata=sum(calibdata.pixels(round(ycen-s):round(ycen+s),floor(xcen-(1+s1)*expecrad):floor(xcen-(1-s1)*expecrad)));
    calibxpeakplus=massaverage1D(xpeakplusdata);
    calibxpeakminus=massaverage1D(xpeakminusdata);
    calibxpeak=((xcen+(1-s1)*expecrad+calibxpeakplus)-(xcen-(1+s1)*expecrad+calibxpeakminus))/2+1; %the +1 -1 is found empirically and is not understood presently.;
    xcenpeak=  ((xcen+(1-s1)*expecrad+calibxpeakplus)+(xcen-(1+s1)*expecrad+calibxpeakminus))/2-1;
    %disp(['calculated horizontal (x) center is ',num2str(xcenpeak)])

    ypeakplusdata=sum((calibdata.pixels(floor(ycen+(1-s1)*expecrad):floor(ycen+(1+s1)*expecrad),round(xcen-s):round(xcen+s)))');
    ypeakminusdata=sum((calibdata.pixels(floor(ycen-(1+s1)*expecrad):floor(ycen-(1-s1)*expecrad),round(xcen-s):round(xcen+s)))');
    calibypeakplus=massaverage1D(ypeakplusdata);
    calibypeakminus=massaverage1D(ypeakminusdata);
    calibypeak=((ycen+(1-s1)*expecrad+calibypeakplus)-(ycen-(1+s1)*expecrad+calibypeakminus))/2+1;
    ycenpeak=((ycen+(1-s1)*expecrad+calibypeakplus)+(ycen-(1+s1)*expecrad+calibypeakminus))/2-1;
    %disp(['calculated horizontal (y) center is ',num2str(ycenpeak)])

    pixel_radius(1)=round(calibxpeak*10)/10;
    pixel_radius(2)=round(calibypeak*10)/10;
    pixel_center(1)=round(xcenpeak*10)/10;
    pixel_center(2)=round(ycenpeak*10)/10;

catch
    errordlg('Auto-routine was unable to find a calibration based on your inputs. Try again or use manual calibration','Auto-Calibration Error','modal')
    pixel_radius=[ ];
    pixel_center=[ ];
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5
function [xcen]=massaverage1D(Y)
%This function calculates the mass-average xcen of any distribution Y provided to it
%It will fail for complicated curves but for simple curves it should do a
%good solid job
%removing a baseline
Ymax=length(Y);
Ybase=(Y(1)+Y(2)+Y(Ymax)+Y(Ymax-1))/4;
Y=Y-Ybase;
%establishing the cumulative distribution
Ysum=sum(Y);
Ycumsumx=cumsum(Y);
%finding the mean
xcen=find(abs((Ycumsumx)-0.5*Ysum) == min(abs((Ycumsumx)-0.5*Ysum)));

% -------------------------------------------------------------
function [xc,yc,R]=solvecircle(P)
%Finds the center and radius of a circle given 3 points on the circle
x1=P(1,1);
x2=P(2,1);
x3=P(3,1);
y1=P(1,2);
y2=P(2,2);
y3=P(3,2);

yc=((x1*x1-x2*x2)/(x1-x2)-(x2*x2-x3*x3)/(x2-x3));
yc=yc+(y1*y1-y2*y2)/(x1-x2)-(y2*y2-y3*y3)/(x2-x3);
yc=yc/((y1-y2)/(x1-x2)-(y2-y3)/(x2-x3))/2;

xc=-yc*(y2-y3)/(x2-x3)+1/2*((x2*x2-x3*x3)/(x2-x3)+(y2*y2-y3*y3)/(x2-x3));

R=sqrt((x1-xc)^2+(y1-yc)^2);
