Matlab psf2otf与fft2函数的关系

本文为转载，原博客地址：http://blog.csdn.net/bluecol/article/details/48288739

psf2otf

---

Convert point-spread function to optical transfer function

Syntax

OTF = psf2otf(PSF) 
OTF = psf2otf(PSF,OUTSIZE)

Description

OTF = psf2otf(PSF) computes the fast Fourier transform (FFT) of the 
point-spread function (PSF) array and creates the optical transfer 
function array, OTF, that is not influenced by the PSF off-centering. 
By default, the OTF array is the same size as the PSF array.

OTF = psf2otf(PSF,OUTSIZE) converts the PSF array into an OTF array, 
where OUTSIZE specifies the size of the OTF array. OUTSIZE cannot be 
smaller than the PSF array size in any dimension.

To ensure that the OTF is not altered because of PSF off-centering, 
psf2otf postpads the PSF array (down or to the right) with 0’s to 
match dimensions specified in OUTSIZE, then circularly shifts the 
values of the PSF array up (or to the left) until the central pixel 
reaches (1,1) position.

Note that this function is used in image convolution/deconvolution 
when the operations involve the FFT.

fft2

---

2-D fast Fourier transform

Syntax

Y = fft2(X) 
Y = fft2(X,m,n)

两者之间的官方解释也说了，psf2otf，其作用是将一个空间点扩散函数转换为频谱面的光学传递函数，执行的也是对PSF的FFT变换。真的只是这样吗？

---

大胆假设 小心求证 
我们用matlab验证一下，首先我们生成PSF函数（其实就是空间卷积嘛，囧），随便来一个拉普拉斯模板吧，So easy。

l = [0 -1 0;-1 4 -1;0 -1 0]
 
 
1

l=⎡⎣⎢0−10−14−10−10⎤⎦⎥

接下来，对 l 进行傅里叶变换，这里为了显示方便，我简单选择m，n都是3。

 F = fft2(l,3,3)
 
 
1

F=⎡⎣⎢0−1.5000−2.5981i−1.5000+2.5981i−1.5000−2.5981i−3.0000+5.1962i6.0000−1.5000+2.5981i6.0000−3.0000−5.1962i⎤⎦⎥

如果使用psf2otf函数呢

 P = psf2otf(l,[3,3])
 
 
1

P=⎡⎣⎢033366666⎤⎦⎥

这里可以看出两者结果是不同的，看下面解释吧，我就不翻译了。 
官方解释如下：

To ensure that the OTF is not altered because of PSF off-centering, 
psf2otf postpads the PSF array (down or to the right) with 0’s to 
match dimensions specified in OUTSIZE, then circularly shifts the 
values of the PSF array up (or to the left) until the central pixel 
reaches (1,1) position.

我们试一下，将 l 的中心元素移动到矩阵的左上角。

S = circshift(l,-floor(size(l)/2))
 
 
1

S=⎡⎣⎢4−1−1−100−100⎤⎦⎥

然后再进行 fft 变换

FS = fft2(S)
 
 
1

FS=⎡⎣⎢033366666⎤⎦⎥

这样两者就一样了。 
实际上psf2otf中最后还舍弃由于浮点数运算带来的微小误差。 
Matlab中psf2otf源码如下：

[psf, psfSize, outSize] = ParseInputs(varargin{:});

if  ~all(psf(:)==0),

   % Pad the PSF to outSize
   padSize = outSize - psfSize;
   psf     = padarray(psf, padSize, 'post');

   % Circularly shift otf so that the "center" of the PSF is at the
   % (1,1) element of the array.
   psf    = circshift(psf,-floor(psfSize/2));

   % Compute the OTF
   otf = fftn(psf);

   % Estimate the rough number of operations involved in the 
   % computation of the FFT.
   nElem = prod(psfSize);
   nOps  = 0;
   for k=1:ndims(psf)
      nffts = nElem/psfSize(k);
      nOps  = nOps + psfSize(k)*log2(psfSize(k))*nffts; 
   end

   % Discard the imaginary part of the psf if it's within roundoff error.
   if max(abs(imag(otf(:))))/max(abs(otf(:))) <= nOps*eps
      otf = real(otf);
   end
else
   otf = zeros(outSize);
end
 
 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31

更多阅读

http://www.imageprocessingplace.com/forumphpBB3/phpBB3/viewtopic.php?f=1&t=9