使用FFT属性查找2D函数的导数

Question

我正在尝试使用FFT属性来获得2D函数的第i个导数 - 特别是2D高斯函数.我也想在MATLAB中用数字做这个.

实际上,我不知道我在做什么,但我在互联网上阅读了很多内容并且所有MATLAB都有帮助,似乎没什么可以帮助我的,所以我会在这里问一下.

这是我到目前为止所拥有的:

h = fspecial('gaussian', size(imr), 10);
imshow(h, []);
G = fft2(h);

这是大小的高斯的0阶导数imr,或512 x 512.我知道我必须做反向和乘法,但我不知道是什么.

任何人都可以给我一个关于使用FFT属性获得2D高斯函数的第i个导数的线索吗？

提前致谢

Answer 1

免责声明 - 2015年9月16日编辑

由于我对频域中微分运算的理解存在一个小但基本的缺陷,这篇文章发生了重大变化.自上次迭代以来,很多帖子都发生了变化.

我要感谢Luis Mendo帮助我调试为什么这不适用于除了这个帖子的原始版本中提供的其他图像之外的其他图像.

当涉及离散时间数据时,没有衍生物这样的东西.衍生物是一种仅在连续时间内存在的分析工具.在离散时间,我们只能通过使用差异来近似.因此,我将假设您的意思是差异操作而不是衍生操作.导数的最常见近似之一是使用前向差分运算.具体地说,对于1D离散序列x[n],y[n]应用前向差分运算的输出序列定义为:

y[n] = x[n+1] - x[n], for n = 1, 2, ..., M-1

M是离散序列中的样本总数.当然,由于前向差分操作的性质,将会有一个较少的价值.

要在频域中完成同样的事情,必须使用离散傅里叶变换属性的移位定理.具体地说,给定信号的DFTed版本x[n],即X(k),以下属性将时域和频域联系在一起:

这里,i表示复数或-1的平方根,并且k是频域中的角频率.F^{-1}是信号的逆傅里叶变换X(k),它是时域信号的傅里叶变换x[n]. M也是信号的长度x[n]和m是要通过移动信号的移动量.因此,这是说,如果你想通过移动计算的换档操作m的位置,你只需要采取傅立叶变换,元素的元素乘以每个部件exp(-i*2*pi*m*k/M),其中k的变换,并采取指数点傅立叶该中间结果的逆傅里叶变换.

请特别注意傅里叶变换的最后一个元素是没有意义的,因为差分运算会产生一个元素少一个的信号,因此从该结果中删除最后一个元素非常重要.

如上所述,为了计算导数的近似值,我们需要找到y[n]或的傅里叶变换Y(k),并且可以这样计算:

请注意,这种转变是-1这样的x[n+1] = x[n-(-1)].因此,您可以计算傅里叶变换,将每个对应的值乘以exp(i*2*pi*k/M) - 1并取反.

进入2D并不是一件容易的事.在这里,我假设我们在两个方向上做差异操作 - 水平和垂直.请记住,我们有两个表示水平和垂直空间坐标的变量.我们也称之为空间域而不是时域,因为变量与我们在2D空间中的位置一致.按照上面的公式,进入2D非常简单:

这里,u和v表示所述2D离散差分操作的空间坐标y[u,v]. U并且V是二维频率分量,并且M和N是2D信号的列和行的数.特别注意您实际上正在进行水平差异操作,然后进行垂直差分操作.具体来说,您将独立地为每一行执行水平差异操作,然后分别对每列执行垂直差异操作.实际上,顺序并不重要.你可以按照你选择的顺序做一个跟随另一个.需要注意的是,信号的二维傅立叶变换保持不变,您只需将每个值乘以一些复数常数即可.但是,您需要确保删除结果的最后一行和最后一列,因为每次操作都会产生比每行和每列的原始长度小一个点的信号,这是由于我们所讨论的差异操作观察.

假设您已经拥有该函数的FFT,您只需要获取每个2D空间坐标并乘以(exp(i*2*pi*U/M) - 1)*(exp(i*2*pi*V/N) - 1).对于存储在其中的FFT的每个2D分量(U,V),我们使用这些相同的频率坐标,并将该位置乘以前面这两个事物的乘积.之后,您将采用逆FFT,我们将其与实际的空间域导数进行比较.我们应该看到它们是一样的.

请注意,当您执行2D FFT时,频率会被标准化,使得它们[-1,1]在两个维度之间跨越.在显示两个域中衍生操作的等效性时,我们需要考虑到这一点.此外,为了使事情变得更加简单,如果您移动频率分量以使它们出现在图像的中心,这也是谨慎的做法.当执行2D FFT时,组件被定位成使得DC分量或原点位于图像的左上角.让我们通过使用将所有内容转移到中心fftshift.

现在,让我们从一些小事做起.为了这个过程,我们假设每个维度的大小是奇数,以允许频率的移动fftshift更容易和明确.假设我们有一个101 x 101高斯滤波器,标准偏差为10,然后我们采用FFT的FFT fftshift.所以:

h = fspecial('gaussian', [101 101], 10);
hf = fftshift(fft2(h));

如果我们想要找到导数,我们定义一个频率点网格,它跨越频域中图像的大小,从0作为原点开始,现在是中心.我们可以通过meshgrid这样做:

[U,V] = meshgrid(-50:50, -50:50);

通常,这个水平和垂直坐标跨越[-floor(cols/2),floor(cols/2)]水平和[-floor(rows/2),floor(rows/2)]垂直.请注意,上面的网格与图像具有相同的尺寸,但中心位于(0,0),并且我们在两个维度上都从-50到50.

现在,在两个方向上进行频域差分运算.让我们做这为双方 x和y,所以在两个方向上的一阶导数:

hf2 = (exp(1i*2*pi*U/101) - 1).*(exp(1i*2*pi*V/101) - 1).*hf;

请注意,由于FFT的频率被归一化,我们必须将频率归一化除以101.101 x 101是由于图像的大小.如果我们想要回到时域,只需采用逆FFT.但是,我们需要撤消ifftshift在进行逆FFT之前我们所做的转变.我们还用于real消除由于数值不准确而导致的任何残余复值分量.你会看到输出是复数值的,但是虚数组件的大小非常小,我们可以忽略它们.因此:

out_freq = real(ifft2(ifftshift(hf2)));

还记得删除最后一行和一列:

out_freq = out_freq(1:end-1,1:end-1);

如果要将此输出与过滤器的时域差异操作进行比较h,我们可以通过使用diff行来完全执行此操作,然后使用此结果,我们将遍历列.因此:

out_spatial = diff(diff(h, 1, 1), 1, 2);

第二个参数表示差异操作的顺序.这是一阶差异,因此水平和垂直差异操作都设置为1.第三个参数告诉您正在进行差分的维度.我们首先在列上执行此操作,然后在行上应用中间结果.然后我们可以展示它们的样子:

figure;
subplot(1,2,1);
imshow(out_spatial, []);
title('Derivative from spatial domain');
subplot(1,2,2);
imshow(out_freq, []);
title('Derivative from frequency domain');

我们得到:

凉!这与我所知的高斯导数一致.

奖金

作为一点奖励,我想向您指出康奈尔大学的一些精彩幻灯片:http://www.cs.cornell.edu/courses/cs6670/2011sp/lectures/lec02_filter.pdf

具体来说,如果您在3D中查看它,高斯看起来就是这样,我们有水平坐标和垂直坐标,并且Z值是3D中高斯的高度:

所以,如果我们独立地采用两个方向的导数就会发生这种情况.顶部显示了空间发生的情况,底部显示了如果我们直接在鸟瞰图上看到这一点会发生什么:

表面的最负面部分是用黑色绘制的,表面的最正面部分是用白色绘制的,其他所有部分都在强度方面进行缩放......所以,如果我们采用一个方向的空间导数,然后使用这个中间结果,并采取另一个方向的空间导数,你会看到我们会得到我们上面看到的.玩弄不同的m和n价值观,看看你得到了什么.我在实践中经常使用的这个属性,但它确实很有用.