From 0541b735ffc64e350f33dfd56f276efaf933baf0 Mon Sep 17 00:00:00 2001 From: Romain Franconville Date: Thu, 4 Dec 2014 11:19:23 -0500 Subject: [PATCH 1/6] Added DFT registration functions --- src/dftRegistration.jl | 176 +++++++++++++++++++++++++++++++++++++++++ 1 file changed, 176 insertions(+) create mode 100644 src/dftRegistration.jl diff --git a/src/dftRegistration.jl b/src/dftRegistration.jl new file mode 100644 index 00000000..b66081a9 --- /dev/null +++ b/src/dftRegistration.jl @@ -0,0 +1,176 @@ +### DFT registration +## Algorithm from (and code very inspired by the matlab code provided by the authors): +## Manuel Guizar-Sicairos, Samuel T. Thurman, and James R. Fienup, +## "Efficient subpixel image registration algorithms," Opt. Lett. 33, 156-158 (2008). + +### Apply dft registration to an image/array (or each frame of such an image if the input is 3D) compared to a reference image (by default the first image of the image stack inputed). If align==false, returns an array with one row per frame in the image, and columns corresponding to the error, the phase difference, the estimated row shift and the estimated column shift. Otherwise returns the registered image. +function dftReg(imgser;ref::AbstractArray=imgser[:,:,1],ufac::Int64=10,align::Bool=false) + ref = fft(ref) + imgF = fft(imgser,(1,2)) + if ndims(imgF) == 2 + results = dftRegfft(data(ref),imgF,ufac).' + else + results = zeros(size(imgF)[3],4) + for i=1:size(imgF)[3] + results[i,:] = dftRegfft(data(ref),imgF[:,:,i],ufac) + end + end + if align + return(alignFromDft(imgser,results)) + else + return(results) + end +end + + +### Align an image using the results of the dftReg function. +function alignFromDft(img2reg::AbstractArray,dftRegRes) + imRes = similar(img2reg,Float64) + img2regF = fft(img2reg,(1,2)) + if ndims(img2reg)==2 + return(subpixelshift(img2regF,dftRegRes[3],dftRegRes[4],dftRegRes[2])) + end + for i=1:length(dftRegRes) + imRes[:,:,i] = subpixelshift(img2regF[:,:,i],dftRegRes[i,3],dftRegRes[i,4],dftRegRes[i,2]) + end + imRes +end + + +### The DFT algorithm : computes misalignment error, phase difference and row and column shifts between two images. The input is the fft transform of the two images, and the oversampling factor usfac. +function dftRegfft(reffft,imgfft,usfac) + if usfac==0 + ## Compute error for no pixel shift + CCmax = sum(reffft.*conj(imgfft)) + rfzero = sumabs2(reffft) + rgzero = sumabs2(imgfft) + error = 1 - CCmax*conj(CCmax)/(rgzero*rfzero) + diffphase = atan2(imag(CCmax),real(CCmax)) + output = ["error" => error, "diffphase" => diffphase] + elseif usfac==1 + ## Whole-pixel shift - Compute crosscorrelation by an IFFT and locate the peak + L = length(reffft) + CC = ifft(reffft.*conj(imgfft)) + loc = indmax(abs(CC)) + CCmax=CC[loc] + rfzero = sumabs2(reffft)/L + rgzero = sumabs2(imgfft)/L + error = abs(1 - CCmax*conj(CCmax)/(rgzero*rfzero)) + diffphase = atan2(imag(CCmax),real(CCmax)) + + (m,n) = size(reffft) + md2 = div(m,2) + nd2 = div(n,2) + locI = ind2sub((m,n),loc) + + if locI[1]>md2 + rowShift = locI[1]-m-1 + else rowShift = locI[1]-1 + end + if locI[2]>nd2 + colShift = locI[2]-n-1 + else colShift = locI[2]-1 + end + output = ["error"=>error,"diffphase"=>diffphase,"rowShift"=>rowShift,"colShift"=>colShift] + else + ## Partial pixel shift + + ##First upsample by a factor of 2 to obtain initial estimate + ##Embed Fourier data in a 2x larger array + dimR = size(reffft) + dimRL = map((x) -> x*2, dimR) + CC = zeros(Complex,dimRL) + CC[(dimR[1]+1-div(dimR[1],2)):(dimR[1]+1+div((dimR[1]-1),2)),(dimR[2]+1-div(dimR[2],2)):(dimR[2]+1+div((dimR[2]-1),2))] = fftshift(reffft).*conj(fftshift(imgfft)) + + ## Compute crosscorrelation and locate the peak + CC = ifft(ifftshift(CC)) + loc = indmax(abs(CC)) + (m,n) = size(CC) + locI = ind2sub((m,n),loc) + CCmax = CC[loc] + + ## Obtain shift in original pixel grid from the position of the crosscorrelation peak + md2 = div(m,2) + nd2 = div(n,2) + + if locI[1] > md2 + rowShift = locI[1]-m-1 + else rowShift = locI[1]-1 + end + + if locI[2] > nd2 + colShift = locI[2]-n-1 + else colShift = locI[2]-1 + end + rowShift = rowShift/2 + colShift = colShift/2 + + ## If upsampling > 2, then refine estimate with matrix multiply DFT + if usfac > 2 + ### DFT Computation ### + # Initial shift estimate in upsampled grid + rowShift = iround(rowShift*usfac)/usfac + colShift = iround(colShift*usfac)/usfac + dftShift = div(ceil(usfac*1.5),2) ## center of output array at dftshift+1 + ## Matrix multiply DFT around the current shift estimate + CC = conj(dftups(imgfft.*conj(reffft),ceil(Int64,usfac*1.5),ceil(Int64,usfac*1.5),usfac,dftShift-rowShift*usfac,dftShift-colShift*usfac))/(md2*nd2*usfac^2) + ## Locate maximum and map back to original pixel grid + loc = indmax(abs(CC)) + locI = ind2sub(size(CC),loc) + CCmax = CC[loc] + rg00 = dftups(reffft.*conj(reffft),1,1,usfac)[1]/(md2*nd2*usfac^2) + rf00 = dftups(imgfft.*conj(imgfft),1,1,usfac)[1]/(md2*nd2*usfac^2) + locI = map((x) -> x - dftShift - 1,locI) + rowShift = rowShift + locI[1]/usfac + colShift = colShift + locI[2]/usfac + else + rg00 = sum(reffft.*conj(reffft))/m/n + rf00 = sum(imgfft.*conj(imgfft))/m/n + end + error = abs(1 - CCmax*conj(CCmax)/(rg00*rf00)) + diffphase = atan2(imag(CCmax),real(CCmax)) + ## If its only one row or column the shift along that dimension has no effect. We set to zero. + if md2 == 1 + rowShift = 0 + end + if nd2==1 + colShift = 0 + end + output = [error,diffphase,rowShift,colShift] + end + output +end + +### Upsampled DFT by matrix multiplies, can compute an upsampled DFT in just a small region +function dftups(inp,nor,noc,usfac=1,roff=0,coff=0) + (nr,nc) = size(inp) + kernc = exp((-1im*2*pi/(nc*usfac))*((ifftshift(0:(nc-1)))-floor(nc/2))*((0:(noc-1))-coff).') + kernr = exp((-1im*2*pi/(nr*usfac))*((0:(nor-1))-roff)*(ifftshift(0:(nr-1))-floor(nr/2)).') + kernr*data(inp)*kernc +end + +### Translate a 2D image/array at subpixel resolution. Outputs the original images translated. If the input is of complex type, it is assumed to be the Fourier transform of the image to shift. +function subpixelshift(img::AbstractArray,rowShift,colShift,diffphase) + img = fft(img) + (nr,nc) = size(img) + Nr = ifftshift((-div(nr,2)):(ceil(Int64,nr/2)-1)) + Nc = ifftshift((-div(nc,2)):(ceil(Int64,nc/2)-1)) + Greg = data(img).* exp(2im*pi*((-rowShift*Nr/nr).-(colShift*Nc/nc).')) + Greg = Greg .* exp(1im*diffphase) + copyproperties(img,real(ifft(Greg))) +end + +function subpixelshift(img::AbstractArray{Complex{Float64},2},rowShift,colShift,diffphase) + (nr,nc) = size(img) + Nr = ifftshift((-div(nr,2)):(ceil(Int64,nr/2)-1)) + Nc = ifftshift((-div(nc,2)):(ceil(Int64,nc/2)-1)) + Greg = data(img) .* exp(2im*pi*((-rowShift*Nr/nr).-(colShift*Nc/nc).')) + Greg = Greg .* exp(1im*diffphase) + copyproperties(img,real(ifft(Greg))) +end + + + + + + From 6900249ba8d7c7bdaec35516c2bf18df6f28c18e Mon Sep 17 00:00:00 2001 From: Romain Franconville Date: Thu, 4 Dec 2014 11:19:23 -0500 Subject: [PATCH 2/6] Added DFT registration functions Adding DFT registration functions --- src/Images.jl | 9 ++- src/dftRegistration.jl | 176 +++++++++++++++++++++++++++++++++++++++++ 2 files changed, 184 insertions(+), 1 deletion(-) create mode 100644 src/dftRegistration.jl diff --git a/src/Images.jl b/src/Images.jl index 413454b9..cd0ed111 100644 --- a/src/Images.jl +++ b/src/Images.jl @@ -46,6 +46,7 @@ include("labeledarrays.jl") include("algorithms.jl") include("connected.jl") include("edge.jl") +include("dftRegistration.jl") __init__() = LibMagick.init() @@ -249,7 +250,13 @@ export # types thin_edges_nonmaxsup_subpix, # phantoms - shepp_logan + shepp_logan, + + # DFT registration + dftReg, + alignFromDft, + subpixelshift + export # Deprecated exports ClipMin, diff --git a/src/dftRegistration.jl b/src/dftRegistration.jl new file mode 100644 index 00000000..b66081a9 --- /dev/null +++ b/src/dftRegistration.jl @@ -0,0 +1,176 @@ +### DFT registration +## Algorithm from (and code very inspired by the matlab code provided by the authors): +## Manuel Guizar-Sicairos, Samuel T. Thurman, and James R. Fienup, +## "Efficient subpixel image registration algorithms," Opt. Lett. 33, 156-158 (2008). + +### Apply dft registration to an image/array (or each frame of such an image if the input is 3D) compared to a reference image (by default the first image of the image stack inputed). If align==false, returns an array with one row per frame in the image, and columns corresponding to the error, the phase difference, the estimated row shift and the estimated column shift. Otherwise returns the registered image. +function dftReg(imgser;ref::AbstractArray=imgser[:,:,1],ufac::Int64=10,align::Bool=false) + ref = fft(ref) + imgF = fft(imgser,(1,2)) + if ndims(imgF) == 2 + results = dftRegfft(data(ref),imgF,ufac).' + else + results = zeros(size(imgF)[3],4) + for i=1:size(imgF)[3] + results[i,:] = dftRegfft(data(ref),imgF[:,:,i],ufac) + end + end + if align + return(alignFromDft(imgser,results)) + else + return(results) + end +end + + +### Align an image using the results of the dftReg function. +function alignFromDft(img2reg::AbstractArray,dftRegRes) + imRes = similar(img2reg,Float64) + img2regF = fft(img2reg,(1,2)) + if ndims(img2reg)==2 + return(subpixelshift(img2regF,dftRegRes[3],dftRegRes[4],dftRegRes[2])) + end + for i=1:length(dftRegRes) + imRes[:,:,i] = subpixelshift(img2regF[:,:,i],dftRegRes[i,3],dftRegRes[i,4],dftRegRes[i,2]) + end + imRes +end + + +### The DFT algorithm : computes misalignment error, phase difference and row and column shifts between two images. The input is the fft transform of the two images, and the oversampling factor usfac. +function dftRegfft(reffft,imgfft,usfac) + if usfac==0 + ## Compute error for no pixel shift + CCmax = sum(reffft.*conj(imgfft)) + rfzero = sumabs2(reffft) + rgzero = sumabs2(imgfft) + error = 1 - CCmax*conj(CCmax)/(rgzero*rfzero) + diffphase = atan2(imag(CCmax),real(CCmax)) + output = ["error" => error, "diffphase" => diffphase] + elseif usfac==1 + ## Whole-pixel shift - Compute crosscorrelation by an IFFT and locate the peak + L = length(reffft) + CC = ifft(reffft.*conj(imgfft)) + loc = indmax(abs(CC)) + CCmax=CC[loc] + rfzero = sumabs2(reffft)/L + rgzero = sumabs2(imgfft)/L + error = abs(1 - CCmax*conj(CCmax)/(rgzero*rfzero)) + diffphase = atan2(imag(CCmax),real(CCmax)) + + (m,n) = size(reffft) + md2 = div(m,2) + nd2 = div(n,2) + locI = ind2sub((m,n),loc) + + if locI[1]>md2 + rowShift = locI[1]-m-1 + else rowShift = locI[1]-1 + end + if locI[2]>nd2 + colShift = locI[2]-n-1 + else colShift = locI[2]-1 + end + output = ["error"=>error,"diffphase"=>diffphase,"rowShift"=>rowShift,"colShift"=>colShift] + else + ## Partial pixel shift + + ##First upsample by a factor of 2 to obtain initial estimate + ##Embed Fourier data in a 2x larger array + dimR = size(reffft) + dimRL = map((x) -> x*2, dimR) + CC = zeros(Complex,dimRL) + CC[(dimR[1]+1-div(dimR[1],2)):(dimR[1]+1+div((dimR[1]-1),2)),(dimR[2]+1-div(dimR[2],2)):(dimR[2]+1+div((dimR[2]-1),2))] = fftshift(reffft).*conj(fftshift(imgfft)) + + ## Compute crosscorrelation and locate the peak + CC = ifft(ifftshift(CC)) + loc = indmax(abs(CC)) + (m,n) = size(CC) + locI = ind2sub((m,n),loc) + CCmax = CC[loc] + + ## Obtain shift in original pixel grid from the position of the crosscorrelation peak + md2 = div(m,2) + nd2 = div(n,2) + + if locI[1] > md2 + rowShift = locI[1]-m-1 + else rowShift = locI[1]-1 + end + + if locI[2] > nd2 + colShift = locI[2]-n-1 + else colShift = locI[2]-1 + end + rowShift = rowShift/2 + colShift = colShift/2 + + ## If upsampling > 2, then refine estimate with matrix multiply DFT + if usfac > 2 + ### DFT Computation ### + # Initial shift estimate in upsampled grid + rowShift = iround(rowShift*usfac)/usfac + colShift = iround(colShift*usfac)/usfac + dftShift = div(ceil(usfac*1.5),2) ## center of output array at dftshift+1 + ## Matrix multiply DFT around the current shift estimate + CC = conj(dftups(imgfft.*conj(reffft),ceil(Int64,usfac*1.5),ceil(Int64,usfac*1.5),usfac,dftShift-rowShift*usfac,dftShift-colShift*usfac))/(md2*nd2*usfac^2) + ## Locate maximum and map back to original pixel grid + loc = indmax(abs(CC)) + locI = ind2sub(size(CC),loc) + CCmax = CC[loc] + rg00 = dftups(reffft.*conj(reffft),1,1,usfac)[1]/(md2*nd2*usfac^2) + rf00 = dftups(imgfft.*conj(imgfft),1,1,usfac)[1]/(md2*nd2*usfac^2) + locI = map((x) -> x - dftShift - 1,locI) + rowShift = rowShift + locI[1]/usfac + colShift = colShift + locI[2]/usfac + else + rg00 = sum(reffft.*conj(reffft))/m/n + rf00 = sum(imgfft.*conj(imgfft))/m/n + end + error = abs(1 - CCmax*conj(CCmax)/(rg00*rf00)) + diffphase = atan2(imag(CCmax),real(CCmax)) + ## If its only one row or column the shift along that dimension has no effect. We set to zero. + if md2 == 1 + rowShift = 0 + end + if nd2==1 + colShift = 0 + end + output = [error,diffphase,rowShift,colShift] + end + output +end + +### Upsampled DFT by matrix multiplies, can compute an upsampled DFT in just a small region +function dftups(inp,nor,noc,usfac=1,roff=0,coff=0) + (nr,nc) = size(inp) + kernc = exp((-1im*2*pi/(nc*usfac))*((ifftshift(0:(nc-1)))-floor(nc/2))*((0:(noc-1))-coff).') + kernr = exp((-1im*2*pi/(nr*usfac))*((0:(nor-1))-roff)*(ifftshift(0:(nr-1))-floor(nr/2)).') + kernr*data(inp)*kernc +end + +### Translate a 2D image/array at subpixel resolution. Outputs the original images translated. If the input is of complex type, it is assumed to be the Fourier transform of the image to shift. +function subpixelshift(img::AbstractArray,rowShift,colShift,diffphase) + img = fft(img) + (nr,nc) = size(img) + Nr = ifftshift((-div(nr,2)):(ceil(Int64,nr/2)-1)) + Nc = ifftshift((-div(nc,2)):(ceil(Int64,nc/2)-1)) + Greg = data(img).* exp(2im*pi*((-rowShift*Nr/nr).-(colShift*Nc/nc).')) + Greg = Greg .* exp(1im*diffphase) + copyproperties(img,real(ifft(Greg))) +end + +function subpixelshift(img::AbstractArray{Complex{Float64},2},rowShift,colShift,diffphase) + (nr,nc) = size(img) + Nr = ifftshift((-div(nr,2)):(ceil(Int64,nr/2)-1)) + Nc = ifftshift((-div(nc,2)):(ceil(Int64,nc/2)-1)) + Greg = data(img) .* exp(2im*pi*((-rowShift*Nr/nr).-(colShift*Nc/nc).')) + Greg = Greg .* exp(1im*diffphase) + copyproperties(img,real(ifft(Greg))) +end + + + + + + From 3fdebea44c2ebbeab1a96f1c6099f57d4a930ac8 Mon Sep 17 00:00:00 2001 From: Romain Franconville Date: Thu, 4 Dec 2014 13:03:41 -0500 Subject: [PATCH 3/6] Typo fixed --- src/dftRegistration.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/dftRegistration.jl b/src/dftRegistration.jl index b66081a9..07dcae4b 100644 --- a/src/dftRegistration.jl +++ b/src/dftRegistration.jl @@ -30,7 +30,7 @@ function alignFromDft(img2reg::AbstractArray,dftRegRes) if ndims(img2reg)==2 return(subpixelshift(img2regF,dftRegRes[3],dftRegRes[4],dftRegRes[2])) end - for i=1:length(dftRegRes) + for i=1:size(dftRegRes)[1] imRes[:,:,i] = subpixelshift(img2regF[:,:,i],dftRegRes[i,3],dftRegRes[i,4],dftRegRes[i,2]) end imRes From 0cd7bcbb3e6f5b402ae49e415e4f1a2bce68bd63 Mon Sep 17 00:00:00 2001 From: Romain Franconville Date: Thu, 4 Dec 2014 17:28:35 -0500 Subject: [PATCH 4/6] Correcting an error in the usfac=1 case. --- src/dftRegistration.jl | 23 ++++++++++++----------- 1 file changed, 12 insertions(+), 11 deletions(-) diff --git a/src/dftRegistration.jl b/src/dftRegistration.jl index 07dcae4b..a9062bcb 100644 --- a/src/dftRegistration.jl +++ b/src/dftRegistration.jl @@ -44,18 +44,19 @@ function dftRegfft(reffft,imgfft,usfac) CCmax = sum(reffft.*conj(imgfft)) rfzero = sumabs2(reffft) rgzero = sumabs2(imgfft) - error = 1 - CCmax*conj(CCmax)/(rgzero*rfzero) + error = sqrt(abs(1 - CCmax*conj(CCmax)/(rgzero*rfzero))) diffphase = atan2(imag(CCmax),real(CCmax)) output = ["error" => error, "diffphase" => diffphase] elseif usfac==1 ## Whole-pixel shift - Compute crosscorrelation by an IFFT and locate the peak L = length(reffft) - CC = ifft(reffft.*conj(imgfft)) + CC = fftshift(reffft).*conj(fftshift(imgfft)) + CC = ifft(ifftshift(CC)) loc = indmax(abs(CC)) CCmax=CC[loc] rfzero = sumabs2(reffft)/L rgzero = sumabs2(imgfft)/L - error = abs(1 - CCmax*conj(CCmax)/(rgzero*rfzero)) + error = sqrt(abs(1 - CCmax*conj(CCmax)/(rgzero*rfzero))) diffphase = atan2(imag(CCmax),real(CCmax)) (m,n) = size(reffft) @@ -71,7 +72,7 @@ function dftRegfft(reffft,imgfft,usfac) colShift = locI[2]-n-1 else colShift = locI[2]-1 end - output = ["error"=>error,"diffphase"=>diffphase,"rowShift"=>rowShift,"colShift"=>colShift] + output = [error,diffphase,rowShift,colShift] else ## Partial pixel shift @@ -118,16 +119,16 @@ function dftRegfft(reffft,imgfft,usfac) loc = indmax(abs(CC)) locI = ind2sub(size(CC),loc) CCmax = CC[loc] - rg00 = dftups(reffft.*conj(reffft),1,1,usfac)[1]/(md2*nd2*usfac^2) - rf00 = dftups(imgfft.*conj(imgfft),1,1,usfac)[1]/(md2*nd2*usfac^2) + rgzero = dftups(reffft.*conj(reffft),1,1,usfac)[1]/(md2*nd2*usfac^2) + rfzero = dftups(imgfft.*conj(imgfft),1,1,usfac)[1]/(md2*nd2*usfac^2) locI = map((x) -> x - dftShift - 1,locI) rowShift = rowShift + locI[1]/usfac colShift = colShift + locI[2]/usfac else - rg00 = sum(reffft.*conj(reffft))/m/n - rf00 = sum(imgfft.*conj(imgfft))/m/n + rgzero = sum(reffft.*conj(reffft))/m/n + rfzero = sum(imgfft.*conj(imgfft))/m/n end - error = abs(1 - CCmax*conj(CCmax)/(rg00*rf00)) + error = sqrt(abs(1 - CCmax*conj(CCmax)/(rgzero*rfzero))) diffphase = atan2(imag(CCmax),real(CCmax)) ## If its only one row or column the shift along that dimension has no effect. We set to zero. if md2 == 1 @@ -156,7 +157,7 @@ function subpixelshift(img::AbstractArray,rowShift,colShift,diffphase) Nr = ifftshift((-div(nr,2)):(ceil(Int64,nr/2)-1)) Nc = ifftshift((-div(nc,2)):(ceil(Int64,nc/2)-1)) Greg = data(img).* exp(2im*pi*((-rowShift*Nr/nr).-(colShift*Nc/nc).')) - Greg = Greg .* exp(1im*diffphase) + Greg = Greg * exp(1im*diffphase) copyproperties(img,real(ifft(Greg))) end @@ -165,7 +166,7 @@ function subpixelshift(img::AbstractArray{Complex{Float64},2},rowShift,colShift, Nr = ifftshift((-div(nr,2)):(ceil(Int64,nr/2)-1)) Nc = ifftshift((-div(nc,2)):(ceil(Int64,nc/2)-1)) Greg = data(img) .* exp(2im*pi*((-rowShift*Nr/nr).-(colShift*Nc/nc).')) - Greg = Greg .* exp(1im*diffphase) + Greg = Greg * exp(1im*diffphase) copyproperties(img,real(ifft(Greg))) end From 21d6467c805f99ae21a69f9c8c9248997e869cae Mon Sep 17 00:00:00 2001 From: Romain Franconville Date: Thu, 4 Dec 2014 17:43:13 -0500 Subject: [PATCH 5/6] Removed the phase difference information, as it is irrelevant for positive valued images. --- src/dftRegistration.jl | 16 ++++++---------- 1 file changed, 6 insertions(+), 10 deletions(-) diff --git a/src/dftRegistration.jl b/src/dftRegistration.jl index a9062bcb..b7b9681f 100644 --- a/src/dftRegistration.jl +++ b/src/dftRegistration.jl @@ -28,10 +28,10 @@ function alignFromDft(img2reg::AbstractArray,dftRegRes) imRes = similar(img2reg,Float64) img2regF = fft(img2reg,(1,2)) if ndims(img2reg)==2 - return(subpixelshift(img2regF,dftRegRes[3],dftRegRes[4],dftRegRes[2])) + return(subpixelshift(img2regF,dftRegRes[2],dftRegRes[3])) end for i=1:size(dftRegRes)[1] - imRes[:,:,i] = subpixelshift(img2regF[:,:,i],dftRegRes[i,3],dftRegRes[i,4],dftRegRes[i,2]) + imRes[:,:,i] = subpixelshift(img2regF[:,:,i],dftRegRes[i,2],dftRegRes[i,3]) end imRes end @@ -57,7 +57,6 @@ function dftRegfft(reffft,imgfft,usfac) rfzero = sumabs2(reffft)/L rgzero = sumabs2(imgfft)/L error = sqrt(abs(1 - CCmax*conj(CCmax)/(rgzero*rfzero))) - diffphase = atan2(imag(CCmax),real(CCmax)) (m,n) = size(reffft) md2 = div(m,2) @@ -72,7 +71,7 @@ function dftRegfft(reffft,imgfft,usfac) colShift = locI[2]-n-1 else colShift = locI[2]-1 end - output = [error,diffphase,rowShift,colShift] + output = [error,rowShift,colShift] else ## Partial pixel shift @@ -129,7 +128,6 @@ function dftRegfft(reffft,imgfft,usfac) rfzero = sum(imgfft.*conj(imgfft))/m/n end error = sqrt(abs(1 - CCmax*conj(CCmax)/(rgzero*rfzero))) - diffphase = atan2(imag(CCmax),real(CCmax)) ## If its only one row or column the shift along that dimension has no effect. We set to zero. if md2 == 1 rowShift = 0 @@ -137,7 +135,7 @@ function dftRegfft(reffft,imgfft,usfac) if nd2==1 colShift = 0 end - output = [error,diffphase,rowShift,colShift] + output = [error,rowShift,colShift] end output end @@ -151,22 +149,20 @@ function dftups(inp,nor,noc,usfac=1,roff=0,coff=0) end ### Translate a 2D image/array at subpixel resolution. Outputs the original images translated. If the input is of complex type, it is assumed to be the Fourier transform of the image to shift. -function subpixelshift(img::AbstractArray,rowShift,colShift,diffphase) +function subpixelshift(img::AbstractArray,rowShift,colShift) img = fft(img) (nr,nc) = size(img) Nr = ifftshift((-div(nr,2)):(ceil(Int64,nr/2)-1)) Nc = ifftshift((-div(nc,2)):(ceil(Int64,nc/2)-1)) Greg = data(img).* exp(2im*pi*((-rowShift*Nr/nr).-(colShift*Nc/nc).')) - Greg = Greg * exp(1im*diffphase) copyproperties(img,real(ifft(Greg))) end -function subpixelshift(img::AbstractArray{Complex{Float64},2},rowShift,colShift,diffphase) +function subpixelshift(img::AbstractArray{Complex{Float64},2},rowShift,colShift) (nr,nc) = size(img) Nr = ifftshift((-div(nr,2)):(ceil(Int64,nr/2)-1)) Nc = ifftshift((-div(nc,2)):(ceil(Int64,nc/2)-1)) Greg = data(img) .* exp(2im*pi*((-rowShift*Nr/nr).-(colShift*Nc/nc).')) - Greg = Greg * exp(1im*diffphase) copyproperties(img,real(ifft(Greg))) end From 82264477146dbe8bc3b9d86fc9322abde2d8e6f3 Mon Sep 17 00:00:00 2001 From: Romain Franconville Date: Thu, 4 Dec 2014 18:10:09 -0500 Subject: [PATCH 6/6] Removed last mention to diffphase. --- src/dftRegistration.jl | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/src/dftRegistration.jl b/src/dftRegistration.jl index b7b9681f..66db7b71 100644 --- a/src/dftRegistration.jl +++ b/src/dftRegistration.jl @@ -45,8 +45,7 @@ function dftRegfft(reffft,imgfft,usfac) rfzero = sumabs2(reffft) rgzero = sumabs2(imgfft) error = sqrt(abs(1 - CCmax*conj(CCmax)/(rgzero*rfzero))) - diffphase = atan2(imag(CCmax),real(CCmax)) - output = ["error" => error, "diffphase" => diffphase] + output = error elseif usfac==1 ## Whole-pixel shift - Compute crosscorrelation by an IFFT and locate the peak L = length(reffft)