可以在图像的非矩形部分进行傅里叶变换

信息处理 matlab 傅里叶变换 噪音 图片
2022-02-02 11:20:03

亲爱的信号处理读者，

我想在图像的某些部分中引入“噪声”。

到目前为止，我出于类似目的使用矩形图像，并在 MATLAB 中使用（逆）傅里叶变换 (2D) 进行了以下操作：

fft = fft2(image);    % run fourier transform on image
amplitude = abs(fft); % extract power spectrum>
newPhase = ( (rand(size(amplitude)) * 2) - 1) * pi;          %generate random phase
newImage = real(ifft2(amplitude. *  exp(1i * newPhase)));    %use both to generate noised image

这使整个图像的相位随机化，保留了功率谱，并且应该保证输入和输出中的空间频谱相同。

我现在想做的是将这种效果仅应用于图像的一部分，而不是矩形（它是脸部的“内部”部分，适合椭圆形）。所以我想从图像中“剪掉”这个椭圆，把它变成上面的“相位噪声”（功率谱/空间频率恒定），但它又回到了我之前从源图像中剪下来的“洞” .

我的“主要解决方案”当然是 FFT/iFFT 也可以在椭圆上运行，所以我可以继续使用上面显示的方法。经过一些谷歌搜索后，我对该选项感到非常悲观 - 有没有希望？

我发现了以下可能有希望的方法的踪迹，但由于我对 FFT 背后的数学不是很熟悉，我会很感激任何提示/解释（为什么）这些想法中的哪些朝着好的方向发展：

  • 以某种方式将椭圆“拉伸”为矩形，按原样应用我的方法，将生成的矩形拉伸回输入椭圆形状。在这里，我担心矩形/椭圆变换可能会破坏匹配的功率谱。您知道在不丢失这些信息的情况下进行转换的好方法吗？（这让我想到了这个想法）

这篇文章还有两个提示：

  • 使用 3D 傅立叶变换，其中（据我所知）第 3 维提供有关应使用图像的哪一部分的信息。任何想法/经验如何实现这一点？

  • 使用“椭圆离散傅里叶变换” - 但我不明白它做什么（以及如何）以及如何实现它。

  • 在这个论坛中，有人向相关问题提出了某种窗口逐渐变细的建议（据我所知，这是为了“过滤”只进入图像的影响，因为感兴趣的椭圆位于黑色背景上）我也不明白，另外我不确定它是否可以帮助我，因为听起来这个问题不需要椭圆作为输出。

提前感谢您的任何建议！

干杯。简

1个回答

即使图像部分是矩形的，FFT 结果也会受到边界效应的污染。零件的 FFT 是矩形零件无限平铺的傅里叶变换。（椭圆形或圆形部分不能平铺。不过，六边形平铺会非常接近。）如果从一个瓷砖到另一个瓷砖之间存在不连续性，这将产生无法描述零件内部的高频分量。随机化阶段将在空间上转换这些组件，以便它们也出现在零件内部，这可能不是您想要的。在边界处镜像并不能消除导数的不连续性，并且需要将 FFT 维度加倍。

这就是平滑渐缩的二维窗函数可以提供帮助的地方。对于圆形零件，您可能喜欢Tukey 窗口，它以与零件中心的距离作为参数。

一维图基窗

对于椭圆，相应地拉伸窗口函数。将数据乘以窗口函数以选择要分析的区域，并对结果进行 FFT。然后进行相位随机化。通过将噪声与椭圆盘状掩码相乘来掩盖噪声，该掩码在圆盘外的值为 0，在圆盘内的值为 1。您还可以使用窗口函数作为掩码来获得模糊边界。我还没有找到一种完美的方法来标准化噪声，但是匹配均值和标准偏差都由掩码加权可以很好地工作：

在此处输入图像描述

将输入图像乘以互补掩码：1 - 掩码，并将两者相加。八度音源：

pkg load signal;
r = 120; # Circle radius
n = 512; # Image width and height
tukeyr = 0.5; # Window transition relative size
lena = double(imread("lena.gif"))/255;
tukeyn = round(r/(1-tukeyr));
x = repmat([1:n], n, 1);
y = repmat([1:n]', 1, n);
dist = sqrt((x-n/2).^2 + (y-n/2).^2);
tukey = tukeywin(tukeyn*2, tukeyr)(tukeyn+1:tukeyn*2);
win = reshape(interp1(1:tukeyn, tukey, reshape(dist, n*n, 1), 0), n, n);
randphase = 2*pi*rand(n, n);
noisefft = abs(fft2(lena.*win)).*exp(j*randphase);
noise = real(ifft2(noisefft));

#Circular mask
weight = r-dist;
weight = (weight >= 1)*1 + (weight > 0 & weight < 1).*weight;
#Or use window as mask
#weight = win;

weightsum = sum(sum(weight));
noiseweightedmean = sum(sum(noise.*weight))/weightsum;
noiseweightedstdev = sqrt(sum(sum((noise-noiseweightedmean).^2.*weight))/weightsum);
lenaweightedmean = sum(sum(lena.*weight))/weightsum;
lenaweightedstdev = sqrt(sum(sum((lena-lenaweightedmean).^2.*weight))/weightsum);
normalizednoise = ((noise - noiseweightedmean) / noiseweightedstdev * lenaweightedstdev) + lenaweightedmean;
imshow(lena.*(1-weight).+normalizednoise.*weight);

PS我发现高通过滤用于噪声生成的Lena图像的副本，通过将其一些最低频率的bin设置为零，与源图像的fft相比，将减少最终图像的fft幅度的误差. 我认为原因是低频幅度波动会减少。在圆（椭圆）内，这些在实践中表现为零频率，并且在没有它们的情况下需要较少的归一化。所以这是要调查的一件事。

其它你可能感兴趣的问题
上一篇OFDM发射机带宽 下一篇将 1D 窗口函数扩展到 3D（或更高）