Matlab psf2otf与fft2函数的关系

Matlab 专栏收录该内容
23 篇文章 0 订阅

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]

l=010141010

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

 F = fft2(l,3,3)

F=01.50002.5981i1.5000+2.5981i1.50002.5981i3.0000+5.1962i6.00001.5000+2.5981i6.00003.00005.1962i

如果使用psf2otf函数呢

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

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))

S=411100100

然后再进行 fft 变换

FS = fft2(S)

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

更多阅读

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

  • 13
    点赞
  • 2
    评论
  • 18
    收藏
  • 一键三连
    一键三连
  • 扫一扫,分享海报

©️2021 CSDN 皮肤主题: 编程工作室 设计师:CSDN官方博客 返回首页
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、C币套餐、付费专栏及课程。

余额充值