Matlab中一组图像的频谱

信息处理 图像处理 matlab fft 频谱
2022-02-03 23:03:00

我希望计算一组图像的频率分布。我基本上热衷于找出图像中频率的变化方式。

我知道我可以使用“fft2”函数来计算二维傅里叶变换。我可以为每个图像执行此操作。但我不确定如何获得所有图像的各向同性频率分布。我想找到一种比较频率分布的方法。

解决这个问题的最佳方法是什么?谢谢你的帮助。

3个回答

如果我正确解释了这个问题,那么您正在寻找一种方法来找到每个图像中的频率分布。听起来您可能需要某种 3D 散点图。

我不知道这是什么情节,但我只是在展示这个想法

如果你有 20 张 256x256 的图像,我想这是如何工作的,

  1. 为每张图像获取 2D FFT 或 2D DCT(或分块 FFT/DCT)
  2. 标准化所有值,以便您可以在每个图像之间进行良好的比较
  3. 将所有频域图像连接成 256 x 256 x 20 矩阵
  4. 过滤掉低于某个阈值的值,因为它们不代表图像的重要组成部分。例如,如果角落中代表高频信息的值的幅度非常低,那么删除它们不会对您的图像产生巨大影响,因为大部分信息都在其他地方。但是删除这些数据点会在您的 3D 绘图中创建空洞,因此您可以看到您感兴趣的信息
  5. 以 3D 形式绘制数据。如果您使用的是 Matlab,则可以使用scatter3()如果您应用热图,该图将更易于阅读。

否则,如果您正在寻找单个图像的频率信息,根据您的问题,有几种方法可以表示数据。如果您澄清了您的问题及其应用程序,我可以编辑此答案。

直方图

如果您获取灰度图像的直方图,您将看到每个像素的强度值分布。

最简单的查看方法是在照片编辑软件(如 Photoshop、GIMP、Picasa)中打开图像,然后从菜单中选择直方图。

如果你想使用 Matlab,hist()函数会绘制一个类似于下面的直方图。您可以直接将图像变量传递给 hist() 函数。但是选择一个好的 bin 大小很重要,因为不同的 bin 大小会极大地改变数据的表示。有多种方法可以计算最佳箱尺寸,但您始终可以猜测和检查。

在此处输入图像描述

二维快速傅里叶变换

在此处输入图像描述

按照 Wikipedia 文章进行理论解释,但基本上 2D FFT 为您提供了低频到高频信息强度的图像。图像的中心是低频信息(如平面、天空、皮肤、墙壁),角落是高频信息(边缘、颗粒状噪声、复杂图案)。

如果您使用的是 Matlab,则需要移动 2D FFT 以将低频信息放在图像的中心。标准化 2D FFT 也是一个好主意,这样您就可以更好地看到细节;这类似于增强照片的对比度以查看细节。

2D DCT

在此处输入图像描述

按照维基百科的文章进行理论解释。2D DCT 与 2D FFT 非常相似;它还为您提供了低频到高频信息强度的图像,但低频信息在左上角高频信息在右下角这使得信息更容易操作(至少对我而言)。

通常,2D DCT 用于图像块(8x8 像素)。这为您提供了仍然与其在空间图像中的位置相关的频率信息。

2D DCT 可用于编码图像。JPEG 编码的一种变体使用 2D DCT 变换。

投影到 1D

这与我使用的类似,它试图将 2D 光谱投影到 1D 上。也就是说,它试图执行

f^pol(r)=02πf^(r,θ)dθ

其中 $[u,v] = [r \cos\theta, r \sin\theta]$ 和 $u,v$ 是光谱坐标。代码中有一部分是频谱乘以 $r$,这是为了尝试调整由于离散化导致的靠近中心的插值误差。然后可以在图像之间比较一维投影。[u,v]=[rcosθ,rsinθ] and u,v are the spectrum coordinates. There is a section in the code where the spectrum is multiplied by r, this is to try to adjust for the interpolation errors close to the centre due to discretisation. The 1D projections can then be compared between images.

close all

% Image
I = imread('cameraman.tif');
imagesc(I); colormap gray; title('Original image'); pause;

% Fourier transform
F = ifftshift(fft2(I))./rows./cols;

% Show spectrum (log)
imagesc(log(abs(F))); title('Fourier transform (abs log)'); pause;

% Grid of FFT coordinates
[rows, cols] = size(F);
[ux, uy] = meshgrid(([1:cols]-(fix(cols/2)+1))/(cols-mod(cols,2)), ...
    ([1:rows]-(fix(rows/2)+1))/(rows-mod(rows,2)));

% Convert to polar coordinates
th = atan2(uy,ux);
r = sqrt(ux.^2 + uy.^2);

% Convert to polar coordinates
Fr = F .* r;
imagesc(abs(Fr)); title('Fourier transform x radius'); pause;
rcoords = linspace(0,sqrt(ux(1,1)^2 + uy(1,1)^2),rows);
thcoords = linspace(0,2*pi,cols);
[ri,thi] = meshgrid(rcoords,thcoords);
[x,y] = pol2cart(thi,ri);
Fp = interp2(ux,uy,abs(Fr),x,y);
imagesc(Fp); title('Fourier transform in polar coordinates'); pause;

% Sum columns to give 1D projection
F1D = sum(Fp);
plot(rcoords,F1D); title('Projection onto 1D'); xlim([0 0.5]);

虽然我没有在图像上使用此特定功能的经验,但有时在通信中我们使用频谱图来查看信号的频率分布如何随时间变化,请查看此 wiki 页面http://en.wikipedia.org/wiki /频谱图

本质上,在大小为 2048 的一维数组上,我们可以采用一组较小的 FFT,例如大小为 128 的大小为 16 的 FFT,并随着时间的推移绘制 fft。其他变化包括在采取 fft 时使用重叠。

由于第四维,这将难以在图像处理中使用。

但是,随着时间的推移,您可以减少某些东西的维度以达到这种效果的一种方法是,对整个图像集的单行图像进行 FFT 并将它们绘制为三维表面,然后对每个相关的图像中的行。