function [pixel_radius,pixel_center] = findcenterbypartialring(S, threepoints)
%SAXS/FINDCENTERBYPARTIALRING Compute q/pixel calibration and center for SAXS
%   image with partial ring.
%
%   [pixel] = FINDCENTERBYPARTIALRING(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_radius is a two-element vector returning the horizontal and vertical radius 
%   of the calibration circle and the


[xcen,ycen,expecrad]=solvecircle(threepoints);
threepointsangle=atan2((threepoints(:,2)-ycen),(threepoints(:,1)-xcen))*360/2/pi;
threepointsangle=threepointsangle+(threepointsangle<0)*360;

%The strategy is now to find improved positions of the 3 clicked points and
% redo the solve center routine

%To do this we convert the image into cylindrical coordinates, (using the approach
% from polar.m - including the use of polcart and interpolate and then
%find the massaverage radius of the peak corresponding to a particular
%angle (for each point)

%To get adequate statistics we take a 3 degree slice. i.e. the cylindrical
%image will have a 120 angular pixels and the the number of radius pixels
%corresponding to the largest number of pixels in the original image.

delth=1;  
threepointanglepixel=round(threepointsangle/delth);

%pixelmax=floor(expecrad*1.2);
%num_radius_pixels=pixelmax;


[ysize xsize] = size(S.pixels);

% Get maximum radius from center.
corners = [1 1 xsize xsize; 1 ysize 1 ysize];
%corners = [xcen-pixelmax xcen-pixelmax xcen+pixelmax xcen+pixelmax; ycen-pixelmax ycen+pixelmax ycen-pixelmax ycen+pixelmax];
rmax = max(sqrt(sum(((corners - repmat([xcen;ycen], 1, 4))).^2)));

% Create arrays containing the x and y pixel coordinates.
[x y] = meshgrid(1:xsize, 1:ysize);

% Create a set of polar coordinates to compute.

[xip yip] = meshgrid(1:rmax,(0:delth:360)*2*pi/360);
[xi yi] = pol2cart(yip, xip);

% Interpolate our polar coordinate values from the image array.
imgrth = interp2((x - xcen), (y - ycen), S.pixels, xi, yi,'*nearest');

s1=1/8;             % this values is the how much "extra" to include around the expected peak. 
% 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.


s=3;
try
    peakdata1=sum(imgrth(max(0,round(threepointanglepixel(1)-s)):min(round(threepointanglepixel(1)+s),360),floor((1-s1)*expecrad):floor((1+s1)*expecrad)));
    peak1=massaverage1D(peakdata1)+floor((1-s1)*expecrad);
    [xpos1,ypos1]=pol2cart(threepointanglepixel(1)/180*pi,peak1);

    peakdata2=sum(imgrth(max(0,round(threepointanglepixel(2)-s)):min(round(threepointanglepixel(2)+s),360),floor((1-s1)*expecrad):floor((1+s1)*expecrad)));
    peak2=massaverage1D(peakdata2)+floor((1-s1)*expecrad);
    [xpos2,ypos2]=pol2cart(threepointanglepixel(2)/180*pi,peak2);
    
    peakdata3=sum(imgrth(max(0,round(threepointanglepixel(3)-s)):min(round(threepointanglepixel(3)+s),360),floor((1-s1)*expecrad):floor((1+s1)*expecrad)));
    peak3=massaverage1D(peakdata3)+floor((1-s1)*expecrad);
    [xpos3,ypos3]=pol2cart(threepointanglepixel(3)/180*pi,peak3);
    
    newpoints=[xpos1+xcen ypos1+ycen; xpos2+xcen ypos2+ycen; xpos3+xcen ypos3+ycen];
    [xcennew,ycennew,expecradnew]=solvecircle(newpoints);

    pixel_radius(1)=round(expecradnew*10)/10;
    pixel_radius(2)=round(expecradnew*10)/10;
    pixel_center(1)=round(xcennew*10)/10;
    pixel_center(2)=round(ycennew*10)/10;
catch
    errordlg('Auto-routine was unable to find a fitting circle based on your inputs. Try again or use manual input ','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);
