前面讲到psf的傅里叶变换尺度与图像不统一而带来的计算上的问题。后面我就根据matlab维纳滤波的源代码进行分析,找出计算流程。
首先从deconvwnr.m开始。函数的编写提供了比较广泛的接口输入,我们抛去这些只关注核心计算步骤。
在deconvwnr.m里还定义了其他两个函数,parse_inputs和CreateNDfrom1D.第一个函数式用来检验输入参数的类型及判别是否正确;第二个函数是针对输入参数有噪声自相关和图像自相关而设计的,如果采用K值来设定滤波器的话,该函数用不到,这里只研究输入为k值的情况。
[I, PSF, ncorr, icorr, sizeI, classI, sizePSF, numNSdim] = parse_inputs(varargin{:});
首先用parse_input函数将输入的参数进行处理,得到一些新的相关变量。假设我们使用的是J = DECONVWNR(I,PSF,NSR)这种形式,NSR就是我们使用的k值。经过parse_input函数对三个输入参数的处理我们也会得到一些新参数。其中I是待处理的图像,PSF是点扩散函数,ncorr为我们输入的k值,icorr为空,sizeI为图像的尺寸,classI图像的数据类型,numNSdim暂时不会用到。
otf = psf2otf(PSF,sizeI);
将psf根据图像的大小变为光学传递函数,在后面的使用中就当做是psf的傅里叶变换。
if isempty(icorr),
检验icorr是否为空,根据我们的使用这一定是空的。
K = ncorr;
nojunk = [1:prod(sizeI)].';% all points are included in the algorithm这两句是icorr为空下执行的。首先给k赋值。举例说明第二句在执行什么,假设sizeI为[300,200],prod(sizeI)=60000,[1:prod(sizeI)]=[1,2,3...,60000],nojunk就是将其转置一下。实际意义就是把图像的每个像素都进行编号,然后把这些编号排列成列状。
跳过一大段涉及自相关函数输入的程序段。进入图像重建。
Denom = abs(otf(nojunk(:))).^2 + K;
这句就是分母计算参考前文所述,分母为:
otf(nojunk(:))找出与图像像素点相对应的在psf傅里叶变换上的点,abs求出其模,然后.^2就是将整个矩阵的元素挨个进行平方,再每个元素分别加上k。这样就按照公式将分母矩阵构建完成了。
Nomin = conj(otf).*fftn(I);
分子是稍加变形过后的结果,||H(u,v)||^2/H(u,v)=H*(u,v),所以分子采用该计算方法。
然后跳过一些检验性质的代码。
JFT(prod(sizeI)) = 0;
创建一个具有和原图像一样像素数量的行向量JFT,初始化为0.
JFT(nojunk(:)) = Nomin(nojunk(:))./Denom;
将JFT的每个位置都赋上值。
J = real(ifftn(reshape(JFT,sizeI)));
将JFT这样的行向量根据图像的size转换形状,再做逆傅里叶变换恢复出原图像。
以上就是维纳滤波根据教材上给出的频率域公式的计算方法。