Finding the proper image shift from phase correlation peak location

I have developed a script in MATLAB which performs a phase correlation of two images, and successfully produces a delta peak. I'm confused on a foolproof method to extracting the peak location as a shift vector. I have tried the following check, however it breaks in some situations:

if X>=1 && X<N/2 || X<1
    dx=-X;
    else
        dx=-X+N+1;
end
 

if Y>=1 && Y<N/2 || Y<1
    dy=-Y+1;
    else
        dy=-Y+N+1;
end

This website describes a bit more background information, but it is non-specific about how to get the shift vector.

Here is the script I use to perform the correlation so you could check:

clear;clc;close all

%WHEN PHASE CORR IMAGE SHOWS USE BRACKETS TO ZOOM, AFTER CLICKING PEAK HIT
%ENTER

a = imread('cameraman.tif'); %load image
b = imtranslate(a,[33 30]); %shift image for example
imcell=cell(1,2);
imcell{1}=a;
imcell{2}=b;
[m,n] = size(a);
xcenter = n/2;
ycenter = m/2;
sigma = n*.65; %choose window termination point

%create window function

for i=1:m
    for j=1:n
        r = sqrt( (j-xcenter)^2 + (i-ycenter)^2);
        if r<=sigma
            g(i,j) = .5 + .5*cos((pi*r)/sigma); %hanning
        else
            g(i,j) = 0;
        end 
    end
end

g=(g/max(max(g)));

fwa = im2uint8(g.*im2double(a));
fwb = im2uint8(g.*im2double(b));

N = round(1.5*length(a)); %choose a multiple of 16 greater than m,n.
m0=round((N-m)/2);
n0=round((N-n)/2);
fca = zeros(N,N);
fcb = zeros(N,N); %make big blank square image 

% center images in blank square

for i=1:N
    for j=1:N 
        if (j>n0 && j<(n0+n-1) && i>m0 && i<(m0+m-1))
             fca(i,j) = fwa(i-m0,j-n0);
             fcb(i,j) = fwb(i-m0,j-n0);
        else
            fca(i,j) = 0;
            fcb(i,j) = 0;
        end
    end
end

imcell=cell(1,2);
imcell{1}=fca;
imcell{2}=fcb;

%% Correlate
Ga = fft2(fca);
Gb = fft2(fcb);
Aga = abs(Ga);
Agb = abs(Gb);
p=.0001*max(Aga,Agb);
q=.0001*max(Aga,Agb);

 
%Rs = (Ga.*conj(Gb))./sqrt((Ga)^2 + (conj(Gb))^2); %normalized cross power spectrum

Rs = (Ga.*conj(Gb))./((Aga+p) + (Agb+q)); %seminormalized cross power spectrum
[etamax,numax] = size(Rs);

%parameters of weighted highpasslowpass filter

r1 = .4;
sigma1 = .2;
r2 = 0.9;
sigma2 = .25;
 
%high pass weight function
for eta=1:etamax
    for nu=1:numax
       p = ((eta-N/2)^2 + (nu-N/2)^2);      
       if 4/N^2 * p < (r1-sigma1)^2
           Hr1s1(eta,nu) = 0;     
       elseif  (r1-sigma1)^2 <= 4/N^2 * p && 4/N^2 * p<r1^2
           Hr1s1(eta,nu) = .5+.5*cos(pi*(r1-(2/N)*sqrt(p))/sigma1);
       else 
            Hr1s1(eta,nu) = 1;
       end    
    end
end
    
%low pass weight function
for eta=1:etamax
    for nu=1:numax
       p = ((eta-N/2)^2 + (nu-N/2)^2);      
       if 4/N^2 * p < r2^2
           Hr2s2(eta,nu) = 1;     
       elseif  r2^2 <= 4/N^2 * p && 4/N^2 * p<(r2+sigma2)^2
           Hr2s2(eta,nu) = .5+.5*cos(pi*(r2-(2/N)*sqrt(p))/sigma2);
       else 
            Hr2s2(eta,nu) = 0;
       end    
    end
end
 

Hr2s2r1s1 = Hr2s2.*Hr1s1;
Rsfilt = Hr2s2r1s1.*Rs;
PCfilt = ifft2(Rsfilt);
PC = ifft2(Rs);

B = imshow(PCfilt,[]);
X = []; Y = [];

while 0<1
    [x,y,b] = ginput(1); 
    if isempty(b)
        break;
    elseif b==91
        ax = axis; width=ax(2)-ax(1); height=ax(4)-ax(3);
        axis([x-width/2 x+width/2 y-height/2 y+height/2]);
        zoom(1/2);
    elseif b==93
        ax = axis; width=ax(2)-ax(1); height=ax(4)-ax(3);
        axis([x-width/2 x+width/2 y-height/2 y+height/2]);
        zoom(2);    
    else
        X=[X;x];
        Y=[Y;y];
    end
end
 
close all

X;
Y;

if X>=1 && X<N/2 || X<1
    dx=-X;
    else
        dx=-X+N+1;
end
 

if Y>=1 && Y<N/2 || Y<1
    dy=-Y+1;
    else
        dy=-Y+N+1;
end 

dx %output shifts
dy

[Francois]

i dont think there us a full prood method to get the peak everytime, but you can fit the shape with a 2D gaussian, there are are garanty that there is only on peak, you can also use the gaussian to interpolate the peak and get sub pixel accuracy. I am also workign on solar eclipse pictures, the interpolation miss the 0,0 peak everytime.

here is my code for the gaussian fit, i hope it helps

% Define the Gaussian function to interpolate Cross power spectrum
gaussianFunction = @(params, xy) params(1) * exp(-((xy(:,1)-params(2)).^2/(2*params(4)^2) + (xy(:,2)-params(3)).^2/(2*params(5)^2))) + params(6); % Define the Gaussian function
[x, y] = meshgrid(1:size(C_real_Crop,2), 1:size(C_real_Crop,1));           % Define the x and y coordinates for interpolation
xy = [x(:), y(:)];
% do the fit
blurKernel = fspecial('gaussian', [12, 12], 3);                            % ADJUSTABLE PARAMETER
Imageblurred{i} = imfilter(C_real_Crop, blurKernel, 'conv', 'replicate');
initialGuess = [1, size(C_real_Crop,2)/2, size(C_real_Crop,1)/2, 10, 10, 0]; % Initial guess for parameters [A, x0, y0, sigma_x, sigma_y, offset]
fitResult = lsqcurvefit(gaussianFunction, initialGuess, xy, double(C_real_Crop(:)));% Fit the Gaussian to the image
C_real_Crop_fit = reshape(gaussianFunction(fitResult, xy), size(C_real_Crop));


% interpolation for sub pixel level accuracy
[maxValue, maxIndex] = max(C_real_Crop_fit(:));                            % Get coordinates around the maximum
[maxRow, maxCol] = ind2sub(size(C_real_Crop_fit), maxIndex);
[x, y] = meshgrid(maxCol-1:0.1:maxCol+1, maxRow-1:0.1:maxRow+1);           % Sub-pixel refinement using interpolation
xy = [x(:), y(:)];
interpolatedValues = gaussianFunction(fitResult, xy);

[~, subPixelIndex] = max(interpolatedValues); % Find the sub-pixel maximum
maxSubPixelIndices = xy(subPixelIndex, :);