7

我想拍摄两张图像并在 Matlab 中使用 2D FFT 将它们卷积在一起而不求助于该conv2函数。但是，我不确定应该如何正确填充矩阵并为卷积做准备。

数学运算如下：

A * B = C

在上面，* 是卷积算子（维基百科链接）。

下面的 Matlab 程序显示了填充和不填充矩阵之间的区别。我怀疑不填充矩阵会导致循环卷积，但我想执行没有混叠的线性卷积。

如果我确实填充了这两个矩阵，那么如何截断卷积的输出以使C与A和B的大小相同？

A = rgb2gray(im2double(imread('1.png'))); % input A
B = rgb2gray(im2double(imread('2.png'))); % kernel B

figure;
imagesc(A); colormap gray;
title ('A')

figure;
imagesc(B); colormap gray;
title ('B')

[m,n] = size(A);
mm = 2*m - 1;
nn = 2*n - 1;

C = (ifft2(fft2(A,mm,nn).* fft2(B,mm,nn)));

figure;
imagesc(C); colormap gray;
title ('C with padding')

C0 = (ifft2(fft2(A).* fft2(B)));

figure;
imagesc(C0); colormap gray;
title ('C without padding')

这是程序的输出：

一个 乙 C C0

4

2 回答 2

12

正如您所指出的，如果没有填充，结果将等同于循环卷积。对于线性卷积，在卷积 2 个图像（2D 信号）A*B 时，完整输出的大小为Ma+Mb-1 x Na+Nb-1，其中Ma x Na, Mb x Nb图像 A 和 B 的大小分别为 。

在填充到预期大小、相乘和变换回来后，通过ifft2，您可以保留结果图像的中心部分（通常对应于 A 和 B 中最大的一个）。 

A = double(imread('cameraman.tif'))./255; % image
B = fspecial('gaussian', [15 15], 2); % some 2D filter function

[m,n] = size(A);
[mb,nb] = size(B); 
% output size 
mm = m + mb - 1;
nn = n + nb - 1;

% pad, multiply and transform back
C = ifft2(fft2(A,mm,nn).* fft2(B,mm,nn));

% padding constants (for output of size == size(A))
padC_m = ceil((mb-1)./2);
padC_n = ceil((nb-1)./2);

% frequency-domain convolution result
D = C(padC_m+1:m+padC_m, padC_n+1:n+padC_n); 
figure; imshow(D,[]);

现在，将上述方法与进行空间域卷积进行比较，使用conv2D

 % space-domain convolution result
 F = conv2(A,B,'same');
 figure; imshow(F,[]);

结果在视觉上是相同的，两者之间的总误差（由于四舍五入）大约为e-10.

于 2012-09-03T22:40:07.927 回答
1

我创建了一个基本上conv2()在频域中的 MATLAB 函数：

function [ mO ] = ImageConvFrequencyDomain( mI, mH, convShape )
% ----------------------------------------------------------------------------------------------- %
% [ mO ] = ImageConvFrequencyDomain( mI, mH, convShape )
% Applies Image Convolution in the Frequency Domain.
% Input:
%   - mI                -   Input Image.
%                           Structure: Matrix.
%                           Type: 'Single' / 'Double' (Single Channel).
%                           Range: (-inf, inf).
%   - mH                -   Filtering Kernel.
%                           Structure: Matrix.
%                           Type: 'Single' / 'Double'.
%                           Range: (-inf, inf).
%   - convShape         -   Convolution Shape.
%                           Sets the convolution shape.
%                           Structure: Scalar.
%                           Type: 'Single' / 'Double'.
%                           Range: {1, 2, 3}.
% Output:
%   - mI                -   Output Image.
%                           Structure: Matrix (Single Channel).
%                           Type: 'Single' / 'Double'.
%                           Range: (-inf, inf).
% References:
%   1.  MATLAB's 'conv2()' - https://www.mathworks.com/help/matlab/ref/conv2.html.
% Remarks:
%   1.  A
% TODO:
%   1.  
%   Release Notes:
%   -   1.0.000     29/04/2021  Royi Avital     RoyiAvital@yahoo.com
%       *   First release version.
% ----------------------------------------------------------------------------------------------- %

CONV_SHAPE_FULL     = 1;
CONV_SHAPE_SAME     = 2;
CONV_SHAPE_VALID    = 3;

numRows     = size(mI, 1);
numCols     = size(mI, 2);

numRowsKernel = size(mH, 1);
numColsKernel = size(mH, 2);

switch(convShape)
    case(CONV_SHAPE_FULL)
        numRowsFft  = numRows + numRowsKernel - 1;
        numColsFft  = numCols + numColsKernel - 1;
        firstRowIdx = 1;
        firstColIdx = 1;
        lastRowIdx  = numRowsFft;
        lastColdIdx = numColsFft;
    case(CONV_SHAPE_SAME)
        numRowsFft  = numRows + numRowsKernel;
        numColsFft  = numCols + numColsKernel;
        firstRowIdx = ceil((numRowsKernel + 1) / 2);
        firstColIdx = ceil((numColsKernel + 1) / 2);
        lastRowIdx  = firstRowIdx + numRows - 1;
        lastColdIdx = firstColIdx + numCols - 1;
    case(CONV_SHAPE_VALID)
        numRowsFft = numRows;
        numColsFft = numCols;
        firstRowIdx = numRowsKernel;
        firstColIdx = numColsKernel;
        % The Kernel when transformed is shifted (Namely its (0, 0) is top
        % left not middle).
        lastRowIdx  = numRowsFft;
        lastColdIdx = numColsFft;
end

mO = ifft2(fft2(mI, numRowsFft, numColsFft) .* fft2(mH, numRowsFft, numColsFft), 'symmetric');
mO = mO(firstRowIdx:lastRowIdx, firstColIdx:lastColdIdx);


end

它完全兼容并经过验证。
完整代码可在我的StackExchange Signal Processing Q74803 GitHub 存储库中找到。

于 2021-04-29T12:17:05.670 回答