function R = x_cntr(im,sz,method, output,th);

% X_CNTR computes the maximum of a correlation function IM with subpix accuracy
% 2*sz+1 is the size of the window over which the center is computed
% The threshold is set automatically set so that the external frame is
% always below the thhreshold.
% Output can be X or Y.
% Xavier Trepat 2007

try
    
    [mi,mj]	= find(im == max(max(im)));
    mi = mi(1);
    mj = mj(1); 
    N = size(im,1);

    switch (method),
        case '2DCentroid';
            th = max([im(max(mi-sz-1,1):min(mi+sz+1,N), max(mj-sz-1,1))'...
                    im(max(mi-sz-1,1):min(mi+sz+1,N), min(mj+sz+1,N))'...
                    im(max(mi-sz-1,1), max(mj-sz-1,1):min(mj+sz+1,N))...
                    im(min(mi+sz+1,N), max(mj-sz-1,1):min(mj+sz+1,N))]); % Set the th automatically.
            
            im1 = im(mi-sz-1:mi+sz+1, mj-sz-1:mj+sz+1); 
            
            im1 = im1-th; 
            im1(find(im1<0))=0;
            grdi = meshgrid(mi-sz-1 : mi+sz+1)'; 
            grdj = meshgrid(mj-sz-1 : mj+sz+1);
            
            xc = sum(sum(im1.*grdj))/sum(sum(im1));
            yc = sum(sum(im1.*grdi))/sum(sum(im1));
            
        case '1DCentroid';
            
            im1=im;
             switch output
                case 'X'
                    th=max(im1(mi,mj-sz-1),im1(mi,mj+sz+1));   
                case 'Y'
                    th=max(im1(mi-sz-1,mj),im1(mi+sz+1,mj));   
            end
            im1 = im1-th;
            im1(find(im1<0))=0;
            grdi = meshgrid(mi-sz : mi+sz,1)'; 
            grdj = meshgrid(mj-sz : mj+sz,1);
            xc = sum(im1(mi,mj-sz:mj+sz).*grdj)/sum(im1(mi,mj-sz:mj+sz));
            yc = sum(im1(mi-sz:mi+sz,mj).*grdi)/sum(im1(mi-sz:mi+sz,mj));
            
        case '2DPoly';
            
            im1 = im(mi-sz+1:mi+sz-1, mj-sz+1:mj+sz-1);     
            im1(find(im1<0))=0;
            grdi = meshgrid(mi-sz+1 : mi+sz-1)'; 
            grdj = meshgrid(mj-sz+1 : mj+sz-1);
            P = polyfitweighted2(mi-sz+1:mi+sz-1, mj-sz+1:mj+sz-1, im1, 2, ones(size(im1)));
            
            i_vals = meshgrid(min2(grdi):0.01:max2(grdi))';
            j_vals = meshgrid(min2(grdj):0.01:max2(grdj));
            im1_vals = polyval2(P,min2(grdi):0.01:max2(grdi),min2(grdj):0.01:max2(grdj));
            
            xc = j_vals(find (im1_vals == max2(im1_vals)));
            yc = i_vals(find (im1_vals == max2(im1_vals)));
            
        case '1DPoly';
            
            im1x = im(mi, mj-sz:mj+sz); 
            im1y = im(mi-sz:mi+sz,mj); 
            grdi = meshgrid(mi-sz : mi+sz,1)'; 
            grdj = meshgrid(mj-sz : mj+sz,1);
            switch output
                case 'X'
                    [a,b] = polyfit(grdj,im1x,2);
                    j_vals = min(grdj):0.01:max(grdj);
                    im1_vals = polyval(a, j_vals);
                    xc = j_vals(find (im1_vals == max(im1_vals)));       
                case 'Y'                    
                    [a,b] = polyfit(grdi,im1y,2);
                    i_vals = min(grdi):0.01:max(grdi);
                    im1_vals = polyval(a, i_vals);
                    yc = i_vals(find (im1_vals == max(im1_vals)));    
            end
    end
    
    
    if(output=='X')
        R = xc-size(im,2)/2-1;    
    elseif(output=='Y')
        R = yc-size(im,2)/2-1;    
    else
        R = NaN;
    end
    
    
catch
    R=NaN;
end
