function center = findcenter(s, quad, status, show)
%SAXS/FINDCENTER Compute symmetric center of SAXS image.
%   This is a back end for CENTER, and for use from other functions, such
%   as inside GUIs.  We don't validate arguments here.
%
%   C = FINDCENTER(S, DIST, STATUS, PLOT) returns a computed center (in raw
%   pixel coordinates) of SAXS image S, searching DIST pixels away from
%   each prospective center for symmetry.  If DIST = 0 FINDCENTER guesses a
%   reasonable value depending on the image size.  If STATUS is nonzero
%   FINDCENTER reports stepwise operations in the command window.  If PLOT
%   is nonzero FINDCENTER plots the image and search points.
%
%   See also SAXS, SAXS/CENTER, SAXS/SETCENTER.

% Parameters.
% KDJ 1/1-2011  By Default uses maximum gradients to locate center for
% PILATUS files
beam_ignore = s.cmask_radius;
img = s.pixels;
[ysize, xsize] = size(s.pixels);
algotype=s.detectortype;

% Starting point.
% If axis location 0,0 is inside the image start the search from the
% nearest pixel.  Otherwise start with the central pixel in the image.
pixorigin = pixel(s, [0, 0]);  % always returns nearest pixel
if any(pixorigin == 1) || any(pixorigin == [xsize, ysize])
    % Origin on the edge of the image.  Move to the middle.
    xs = round(ysize/2);  % starting x coordinate
    ys = round(xsize/2);  % starting y coordinate
else
    xs = pixorigin(2);
    ys = pixorigin(1);
end

if ~isempty(strfind(algotype,'PILATUS'))
     d2nearedge=2*floor(min([abs(ysize-xs) xs abs(xsize-ys) ys])/2);
     qperpix=s.kcal./s.pixelcal;
     beam_ignore=10;
     npix=d2nearedge/2-beam_ignore;
    % first generate the radial average from 350 to 10 degree.
      i350=average_raw(s, 'x',[beam_ignore*qperpix(2) 350 d2nearedge/2*qperpix(2) 10],npix);
      i0=average_raw(s, 'x',[beam_ignore*qperpix(2) 0 d2nearedge/2*qperpix(2) 10],npix);
    % then generate the radial average from 180-190 
      i170=average_raw(s, 'x',[beam_ignore*qperpix(2) 170 d2nearedge/2*qperpix(2) 10],npix);

      % then generate the radial average from 60-80 
      i60=average_raw(s, 'x',[beam_ignore*qperpix(2) 60 d2nearedge/2*qperpix(2) 10],npix);
    % then generate the radial average from 240-260
      i240=average_raw(s, 'x',[beam_ignore*qperpix(2) 240 d2nearedge/2*qperpix(2) 10],npix);
    % now smooth
      i350g=gradient(conv((i350.pixels+i0.pixels)/2,[1 1 1 1 1]/5,'valid'));
      i170g=gradient(conv(i170.pixels,[1 1 1 1 1]/5,'valid'));
      i60g=gradient(conv(i60.pixels,[1 1 1 1 1]/5,'valid'));
      i240g=gradient(conv(i240.pixels,[1 1 1 1 1]/5,'valid'));
      xaxisfullg=conv(i240.xaxisfull,[1 1 1 1 1]/5,'valid');
      
      %figure 
      %plot(i350.xaxisfullg,i350.pixels,i170.xaxisfull,i170.pixels,xaxisfull,i350s,xaxisfull,i170s)
      %figure 
      %plot(i80.xaxisfull,i80.pixels,i260.xaxisfull,i260.pixels,xaxisfull,i260s,xaxisfull,i260s)
      
      
      %figure 
      %plot(xaxisfullg,i350g,xaxisfullg,i170g)
      %figure 
      %plot(xaxisfullg,i60g,xaxisfullg,i240g)
      
      %now find a minimum of the gradient for each       
      %and perform a polymomial fit to get the best value for the minimum
      mintemp=max(2,find(i350g==min(i350g),1));
      P=polyfit(xaxisfullg(mintemp-1:mintemp+1),i350g(mintemp-1:mintemp+1),2);
      min350=-P(2)/2/P(1);
      
      mintemp=max(2,find(i170g==min(i170g),1)); %
      P=polyfit(xaxisfullg(mintemp-1:mintemp+1),i170g(mintemp-1:mintemp+1),2);
      min170=-P(2)/2/P(1);
      
      mintemp=max(2,find(i60g==min(i60g),1));
      P=polyfit(xaxisfullg(mintemp-1:mintemp+1),i60g(mintemp-1:mintemp+1),2);
      min60=-P(2)/2/P(1);
      
      mintemp=max(2,find(i240g==min(i240g),1));
      P=polyfit(xaxisfullg(mintemp-1:mintemp+1),i240g(mintemp-1:mintemp+1),2);
      min240=-P(2)/2/P(1);
      

      if (mintemp==2)
         %polynomial fit did not do the job since there was no "peak" in
         %the data. Will insted compar intensities.
         minval=min([i350.pixels(2) i170.pixels(2) i60.pixels(2) i240.pixels(2)]);
         
         mintemp=find(i350.pixels==min(abs(i350.pixels-minval)).*sign(i350.pixels-minval)+minval,1);
         xx=i350.xaxisfull;yy=i350.pixels-minval;
         min350 = fzero(@(x) interp1(xx,yy,x),i350.xaxisfull(mintemp));
         
         mintemp=find(i170.pixels==min(abs(i170.pixels-minval)).*sign(i170.pixels-minval)+minval,1)
         xx=i170.xaxisfull;yy=i170.pixels-minval;
         min170 = fzero(@(x) interp1(xx,yy,x),i170.xaxisfull(mintemp));
         
         mintemp=find(i60.pixels==min(abs(i60.pixels-minval)).*sign(i60.pixels-minval)+minval,1)
         xx=i60.xaxisfull;yy=i60.pixels-minval;
         min60 = fzero(@(x) interp1(xx,yy,x),i60.xaxisfull(mintemp));
         
         mintemp=find(i240.pixels==min(abs(i240.pixels-minval)).*sign(i240.pixels-minval)+minval,1)
         xx=i240.xaxisfull;yy=i240.pixels-minval;
         min240 = fzero(@(x) interp1(xx,yy,x),i240.xaxisfull(mintemp));
          
      end
      ycorr=(min350-min170)/2/qperpix(2);
      xcorr=(min60-min240)/2/qperpix(1);
      
      
      a=xs+xcorr;
      b=ys+ycorr;

      
      
    
    
else   
    if status
        disp(['Search ', num2str(xsize), ' x ', num2str(ysize), ' image.'])
        disp(['Begin at x = ', num2str(xs), ', y = ', num2str(ys), '.'])
        if show
            disp('Plot image and search points.')
        end
        disp(['Ignore ', num2str(2 * beam_ignore), ' x ', ...
            num2str(2 * beam_ignore), ' box around candidate centers.'])
        disp(['Use ', num2str(quad), ' x ', ...
            num2str(quad), ' quadrants to test symmetry.'])
    end
    
    if show, imagexrd(img, 'log'), hold on, end
    
    if status
        disp('Simplexes...')
    end
    
    objf = zeros(3,1);
    
    % Define first simplex
    curp = [xs ys
        xs+1 ys
        xs   ys+1];
    
    % Calcuate objective function for first simplex
    for i = 1:3
        objf(i) = imcnobjf(img,curp(i,1),curp(i,2),quad,beam_ignore,algotype);
    end
    
    if show, plot(curp(:,2)+.5,curp(:,1)+.5,'+m'), end  % plot simplex on image
    if status, fprintf(1,'.'), end
    
    lastx = max(objf);
    xcur = min(objf);
    
    % Get the simplex going as long as the new point is better than the old
    while xcur < lastx
        % Determine the max point
        [maxo,ind] = max(objf);
        % Calculate a new point by flipping the simplex
        if ind == 1
            newp = curp(1,:) - 2*(curp(1,:) - mean(curp([2 3],:)));
        elseif ind == 2
            newp = curp(2,:) - 2*(curp(2,:) - mean(curp([1 3],:)));
        elseif ind == 3
            newp = curp(3,:) - 2*(curp(3,:) - mean(curp([1 2],:)));
        end
        % Replace the old max with the new point
        curp(ind,:) = newp;
        % Plot the new point on the image
        if show, plot(curp(ind,2)+.5,curp(ind,1)+.5,'+m'), end
        if status, fprintf(1,'.'), end
        % Evaluate the objective function on the new point
        objf(ind) = imcnobjf(img,curp(ind,1),curp(ind,2),quad,beam_ignore,algotype);
        % Set xcur to be the new point
        xcur = objf(ind);
        % Set lastx to be the former maximum
        lastx = maxo;
    end
    if status, fprintf(1,'\n'), end
    
    % Find the minimum point
    [mino,ind] = min(objf);
    a = curp(ind,1);
    b = curp(ind,2);
    
    % Now search around the point
    if status, disp('Grid search...'), end
    objf = zeros(25,1);
    xcp = zeros(25,1);
    ycp = zeros(25,1);
    eflag = 1;
    while eflag == 1
        k = 0;
        for i = 1:5
            for j = 1:5
                k = k+1;
                objf(k) = imcnobjf(img,a-3+i,b-3+j,quad,beam_ignore,algotype);
                if show, plot(b-3+j+.5,a-3+i+.5,'og'), end
                if status, fprintf(1,'.'), end
                xcp(k) = a-3+i;
                ycp(k) = b-3+j;
            end
        end
        objfm = reshape(objf,5,5);
        %[xcp ycp objf/1000000000]
        [mx,my] = find(objfm==min(objf));
        [t,mo] = min(objf);
        b = ycp(mo);
        a = xcp(mo);
        if show, plot(b+.5,a+.5,'xy'), end
        % Check to make sure minimum point isn't on an edge
        if ((mx~=1)&(mx~=5))&((my~=1)&(my~=5))
            %disp('not on an edge')
            if status, disp(' '), disp('Quadratic fit...'), end
            dat = objfm(mx-1:mx+1,my-1:my+1);
            xdat = [-1 -1 -1 0 0 0 1 1 1; -1 0 1 -1 0 1 -1 0 1]';
            coeff = [ones(9,1) xdat xdat.^2]\dat(:);
            xs = (-coeff(2)/coeff(4));
            ys = (-coeff(3)/coeff(5));
            eflag = 0;
        end
    end
    if status, fprintf(1,'\n'), end
    if abs(xs)<1 & abs(ys)<1
        a = a + .5 + xs;
        b = b + .5 + ys;
    else
        a = a + .5;
        b = b + .5;
        %disp(objfm), disp([xs ys])
    end
end
if show
    hold off
    title('Image with Estimated Center Shown along with Test Points')
    hline(a), vline(b)
end

if status
    disp(['Calculated center column = ', num2str(b)])
    disp(['Calculated center row    = ', num2str(a)])
end

% % Return center coordinates in axes units.
% center = axiscoords(s, [b, a]);
% No, just return raw pixel coordinates.
center = [b, a];

% -------------------------------------------------------------
% New parameterized objective function.
function objf = imcnobjf(img,xc,yc,quad,ignore,algotype)
%IMCNOBJF Calculates objective function for XRD center finder
ul = img(xc-(quad-1):xc,yc-(quad-1):yc);
lr = img(xc+1:xc+quad,yc+1:yc+quad);
ur = img(xc+1:xc+quad,yc-(quad-1):yc);
ll = img(xc-(quad-1):xc,yc:yc+(quad-1));

if strcmp(algotype,'PILATUS 300K')
    dif1 = double(ul) - double(flipud(ur));
    dif1((quad-(ignore-1)):quad,(quad-(ignore-1)):quad) = 0;
    %dif1(isnan(dif1))=0;
    dif2 = double(ul) - double(fliplr(ll));
    dif2(1:ignore,(quad-(ignore-1)):quad) = 0;
    objf = sum(dif1(~isnan(dif1)).^2) + sum(dif2(~isnan(dif2)).^2);
else
    dif1 = double(ul) - double(rot90(lr,2));
    dif1((quad-(ignore-1)):quad,(quad-(ignore-1)):quad) = 0;
    dif2 = double(ur) - double(rot90(ll,2));
    dif2(1:ignore,(quad-(ignore-1)):quad) = 0;
    objf = sum(dif1(:).^2) + sum(dif2(:).^2);
end

