| 查看: 1279 | 回复: 0 | ||
xiashaoyan新虫 (正式写手)
|
[求助]
基于相位相关和重采样的亚像素图像配准算法问题
|
|
在这里，亚像素图像配准算法是采用的方法是相位相关和重采样，实现的程序参考网址为：http://www.mathworks.com/matlabc ... t/dftregistration.m，所有的程序将在下面贴出来，程序中有几个地方看了很长时间看不明白，所以写了这个帖子，希望研究过亚像素图像配准算法或者看懂程序的前辈能给小弟指教一下，疑惑的问题将在附件中的程序疑惑1,2,3，具体的疑惑为： 疑惑1： 在讨论中，有多个地方出现了(出现的位置在程序疑惑1和2用红色方框标记出来了) if rloc > md2 row_shift = rloc - m - 1; else row_shift = rloc - 1; end if cloc > nd2 col_shift = cloc - n - 1; else col_shift = cloc - 1; end 这里的代码的意思，好像是平移的意思，但是具体为什么这样操作思索了很长时间还是不明白？？？？ 疑惑2：(出现的位置在程序疑惑3用红色方框标记出来了) 程序中 buf1ft是基准图像,buf2ft待匹配的图像，usfac是图像的采样倍数，程序中讨论了usfac=1时，用FFT计算，而usfac不等于1的计算方法不同于usfac==1的情况，另外usfac>2时是用DFT计算的，所以这里我就疑惑了，为什么讨论usfac不同的值时，计算的方法不一样？？？ 为了方便大家理解，我截出了亚像素图像配准的计算过程图，见附件 亚像素图像匹配算法的计算过程.PNG 程序的全部代码： function [output Greg] = dftregistration(buf1ft,buf2ft,usfac) % function [output Greg] = dftregistration(buf1ft,buf2ft,usfac); % Efficient subpixel image registration by crosscorrelation. This code % gives the same precision as the FFT upsampled cross correlation in a % small fraction of the computation time and with reduced memory % requirements. It obtains an initial estimate of the crosscorrelation peak % by an FFT and then refines the shift estimation by upsampling the DFT % only in a small neighborhood of that estimate by means of a % matrix-multiply DFT. With this procedure all the image points are used to % compute the upsampled crosscorrelation. % Manuel Guizar - Dec 13, 2007 % Portions of this code were taken from code written by Ann M. Kowalczyk % and James R. Fienup. % J.R. Fienup and A.M. Kowalczyk, "Phase retrieval for a complex-valued % object by using a low-resolution image," J. Opt. Soc. Am. A 7, 450-458 % (1990). % Citation for this algorithm: % Manuel Guizar-Sicairos, Samuel T. Thurman, and James R. Fienup, % "Efficient subpixel image registration algorithms," Opt. Lett. 33, % 156-158 (2008). % Inputs % buf1ft Fourier transform of reference image, % DC in (1,1) [DO NOT FFTSHIFT] % buf2ft Fourier transform of image to register, % DC in (1,1) [DO NOT FFTSHIFT] % usfac Upsampling factor (integer). Images will be registered to % within 1/usfac of a pixel. For example usfac = 20 means the % images will be registered within 1/20 of a pixel. (default = 1) % Outputs % output = [error,diffphase,net_row_shift,net_col_shift] % error Translation invariant normalized RMS error between f and g % diffphase Global phase difference between the two images (should be % zero if images are non-negative). % net_row_shift net_col_shift Pixel shifts between images % Greg (Optional) Fourier transform of registered version of buf2ft, % the global phase difference is compensated for. % Default usfac to 1 if exist('usfac')~=1, usfac=1; end % Compute error for no pixel shift if usfac == 0, CCmax = sum(sum(buf1ft.*conj(buf2ft))); rfzero = sum(abs(buf1ft(  ).^2);rgzero = sum(abs(buf2ft( ).^2); error = 1.0 - CCmax.*conj(CCmax)/(rgzero*rfzero); error = sqrt(abs(error)); diffphase=atan2(imag(CCmax),real(CCmax)); output=[error,diffphase]; % Whole-pixel shift - Compute crosscorrelation by an IFFT and locate the % peak elseif usfac == 1, [m,n]=size(buf1ft); CC = ifft2(buf1ft.*conj(buf2ft)); [max1,loc1] = max(CC); [max2,loc2] = max(max1); rloc=loc1(loc2); cloc=loc2; CCmax=CC(rloc,cloc); rfzero = sum(abs(buf1ft( ).^2)/(m*n);rgzero = sum(abs(buf2ft( ).^2)/(m*n); error = 1.0 - CCmax.*conj(CCmax)/(rgzero(1,1)*rfzero(1,1)); error = sqrt(abs(error)); diffphase=atan2(imag(CCmax),real(CCmax)); md2 = fix(m/2); nd2 = fix(n/2); if rloc > md2 row_shift = rloc - m - 1; else row_shift = rloc - 1; end if cloc > nd2 col_shift = cloc - n - 1; else col_shift = cloc - 1; end output=[error,diffphase,row_shift,col_shift]; % Partial-pixel shift else % First upsample by a factor of 2 to obtain initial estimate % Embed Fourier data in a 2x larger array [m,n]=size(buf1ft); mlarge=m*2; nlarge=n*2; CC=zeros(mlarge,nlarge); CC(m+1-fix(m/2):m+1+fix((m-1)/2),n+1-fix(n/2):n+1+fix((n-1)/2)) = ... fftshift(buf1ft).*conj(fftshift(buf2ft)); % Compute crosscorrelation and locate the peak CC = ifft2(ifftshift(CC)); % Calculate cross-correlation [max1,loc1] = max(CC); [max2,loc2] = max(max1); rloc=loc1(loc2);cloc=loc2; CCmax=CC(rloc,cloc); % Obtain shift in original pixel grid from the position of the % crosscorrelation peak [m,n] = size(CC); md2 = fix(m/2); nd2 = fix(n/2); if rloc > md2 row_shift = rloc - m - 1; else row_shift = rloc - 1; end if cloc > nd2 col_shift = cloc - n - 1; else col_shift = cloc - 1; end row_shift=row_shift/2; col_shift=col_shift/2; % If upsampling > 2, then refine estimate with matrix multiply DFT if usfac > 2, %%% DFT computation %%% % Initial shift estimate in upsampled grid row_shift = round(row_shift*usfac)/usfac; col_shift = round(col_shift*usfac)/usfac; dftshift = fix(ceil(usfac*1.5)/2); %% Center of output array at dftshift+1 % Matrix multiply DFT around the current shift estimate CC = conj(dftups(buf2ft.*conj(buf1ft),ceil(usfac*1.5),ceil(usfac*1.5),usfac,... dftshift-row_shift*usfac,dftshift-col_shift*usfac))/(md2*nd2*usfac^2); % Locate maximum and map back to original pixel grid [max1,loc1] = max(CC); [max2,loc2] = max(max1); rloc = loc1(loc2); cloc = loc2; CCmax = CC(rloc,cloc); rg00 = dftups(buf1ft.*conj(buf1ft),1,1,usfac)/(md2*nd2*usfac^2); rf00 = dftups(buf2ft.*conj(buf2ft),1,1,usfac)/(md2*nd2*usfac^2); rloc = rloc - dftshift - 1; cloc = cloc - dftshift - 1; row_shift = row_shift + rloc/usfac; col_shift = col_shift + cloc/usfac; % If upsampling = 2, no additional pixel shift refinement else rg00 = sum(sum( buf1ft.*conj(buf1ft) ))/m/n; rf00 = sum(sum( buf2ft.*conj(buf2ft) ))/m/n; end error = 1.0 - CCmax.*conj(CCmax)/(rg00*rf00); error = sqrt(abs(error)); 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, row_shift = 0; end if nd2 == 1, col_shift = 0; end output=[error,diffphase,row_shift,col_shift]; end % Compute registered version of buf2ft if (nargout > 1)&&(usfac > 0), [nr,nc]=size(buf2ft); Nr = ifftshift([-fix(nr/2):ceil(nr/2)-1]); Nc = ifftshift([-fix(nc/2):ceil(nc/2)-1]); [Nc,Nr] = meshgrid(Nc,Nr); Greg = buf2ft.*exp(i*2*pi*(-row_shift*Nr/nr-col_shift*Nc/nc)); Greg = Greg*exp(i*diffphase); elseif (nargout > 1)&&(usfac == 0) Greg = buf2ft*exp(i*diffphase); end return function out=dftups(in,nor,noc,usfac,roff,coff) % function out=dftups(in,nor,noc,usfac,roff,coff); % Upsampled DFT by matrix multiplies, can compute an upsampled DFT in just % a small region. % usfac Upsampling factor (default usfac = 1) % [nor,noc] Number of pixels in the output upsampled DFT, in % units of upsampled pixels (default = size(in)) % roff, coff Row and column offsets, allow to shift the output array to % a region of interest on the DFT (default = 0) % Recieves DC in upper left corner, image center must be in (1,1) % Manuel Guizar - Dec 13, 2007 % Modified from dftus, by J.R. Fienup 7/31/06 % This code is intended to provide the same result as if the following % operations were performed % - Embed the array "in" in an array that is usfac times larger in each % dimension. ifftshift to bring the center of the image to (1,1). % - Take the FFT of the larger array % - Extract an [nor, noc] region of the result. Starting with the % [roff+1 coff+1] element. % It achieves this result by computing the DFT in the output array without % the need to zeropad. Much faster and memory efficient than the % zero-padded FFT approach if [nor noc] are much smaller than [nr*usfac nc*usfac] [nr,nc]=size(in); % Set defaults if exist('roff')~=1, roff=0; end if exist('coff')~=1, coff=0; end if exist('usfac')~=1, usfac=1; end if exist('noc')~=1, noc=nc; end if exist('nor')~=1, nor=nr; end % Compute kernels and obtain DFT by matrix products kernc=exp((-i*2*pi/(nc*usfac))*( ifftshift([0:nc-1]).' - floor(nc/2) )*( [0:noc-1] - coff )); kernr=exp((-i*2*pi/(nr*usfac))*( [0:nor-1].' - roff )*( ifftshift([0:nr-1]) - floor(nr/2) )); out=kernr*in*kernc; return  亚像素图像匹配算法的计算过程.PNG  程序疑惑1.PNG  程序疑惑2.PNG  程序疑惑3.PNG |
回复此楼» 猜你喜欢
请问有评职称，把科研教学业绩算分排序的高校吗
已经有6人回复
-
2025冷门绝学什么时候出结果
已经有6人回复
-
Bioresource Technology期刊，第一次返修的时候被退回好几次了
已经有7人回复
-
真诚求助：手里的省社科项目结项要求主持人一篇中文核心，有什么渠道能发核心吗
已经有8人回复
-
寻求一种能扛住强氧化性腐蚀性的容器密封件
已经有5人回复
-
请问哪里可以有青B申请的本子可以借鉴一下。
已经有4人回复
-
请问下大家为什么这个铃木偶联几乎不反应呢
已经有5人回复
-
天津工业大学郑柳春团队欢迎化学化工、高分子化学或有机合成方向的博士生和硕士生加入
已经有4人回复
-
康复大学泰山学者周祺惠团队招收博士研究生
已经有6人回复
-
AI论文写作工具：是科研加速器还是学术作弊器？
已经有3人回复
