De convolution linéaire de deux images sous Matlab à l'aide de fft2

Je voudrais prendre deux images et de convolution ensemble dans Matlab à l'aide de la FFT 2D sans avoir recours à la conv2 fonction. Cependant, je suis incertain de la façon dont les matrices doivent être matelassés et préparé pour le produit de convolution.

L'opération mathématique est la suivante:

UN * B = C

Ci-dessus, * est le produit de convolution de l'opérateur (Lien wikipédia).

La suite de programme Matlab montre la différence entre le rembourrage et de ne pas remplir les matrices. Je les soupçonne de ne pas remplir les matrices de résultats dans une circulaire de convolution, mais je voudrais effectuer une convolution linéaire sans aliasing.

Si je ne pad les deux matrices, alors comment puis-je tronquer la sortie de la convolution de sorte que C est de la même taille que Un et 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')

Ici est la sortie du programme:

De convolution linéaire de deux images sous Matlab à l'aide de fft2
De convolution linéaire de deux images sous Matlab à l'aide de fft2
De convolution linéaire de deux images sous Matlab à l'aide de fft2
De convolution linéaire de deux images sous Matlab à l'aide de fft2

OriginalL'auteur Nicholas Kinar | 2012-09-03

  1. 11

    Sans rembourrage, le résultat sera équivalent à convolution circulaire comme vous le soulignez. Pour de convolution linéaire, en réalisant le produit de convolution de 2 images (2D signaux) A*B la sortie complète sera de taille Ma+Mb-1 x Na+Nb-1, où Ma x Na, Mb x Nb la taille des images A et B resp.

    Après un rembourrage à la taille attendue, la multiplication et la transformation de retour, via ifft2, vous pouvez garder la partie centrale de l'image résultante (correspondant le plus souvent à la plus grande de A et 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,[]);

    Maintenant, comparez cela à faire spatiale-domaine de convolution, à l'aide de conv2D

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

    Résultats sont visuellement le même, et l'erreur totale entre les deux (en raison de l'arrondissement) de l'ordre de e-10.

    C'est tout simplement merveilleux, gevang. Merci beaucoup pour une réponse complète! Une question: Depuis l'opération de convolution est liée à la corrélation, je ferais le même pour le 2D corrélation des deux matrices? Serait ifft2(fft2(M1,mm,nn).*fft2(fliplr(flipud(M2))),mm,nn) être le code Matlab pour la corrélation, et puis j'aimerais garder la partie centrale de l'image, de la même manière comme indiqué dans votre réponse ci-dessus?
    En outre, si A et B ont été complexes, les matrices, serait le produit de convolution (ou de corrélation) code de l'opération fonctionne exactement de la même manière?
    Depuis la convolution (et à transformée de Fourier) sont linéaires des opérations et de distribution avec plus, l'équivalence tiendra pour les signaux de la forme A + Aj, c'est à dire que vous aurez une somme de circonvolutions entre les combinaisons des parties réelles et imaginaires des images de la taille d'origine. Vous pouvez le vérifier en remplaçant A et B avec des matrices complexes ci-dessus. Si les deux sont complexes, de visualiser le abs() de la sortie ou de comparer la différence entre le real() et imag() parties de F et D.
    Merci encore, gevang. Oui, après quelques essais, je peux confirmer que le code de convolution travaille dans Matlab pour les signaux de la forme A + Aj.
    En outre, en utilisant le logiciel Matlab xcorr2 fonctionner comme un plus lents test, il semble que la 2D corrélation de la A et B les matrices peuvent être calculées par ifft2(fft2(A,mm,nn).*fft2(fliplr(flipud(B)),mm,nn)); ou, de manière équivalente ifft2(fft2(A,mm,nn).*fft2(rot90(B,2),mm,nn)). Une référence pour ce qui est de la DSP livre de Steven Smith: dspguide.com/ch24/6.htm.

    OriginalL'auteur gevang