[37] | 1 | function [out]=fftolamopt2(a,b,opt,shape) |
---|
| 2 | % [out]=fftolamopt2(a,b,siz1,siz2,shape) |
---|
| 3 | % |
---|
| 4 | % Overlap-add method FFT-based 2D convolution |
---|
| 5 | % Example: |
---|
| 6 | % load fftexecutiontimes; % load FFTrv, FFTiv and IFFTiv in workspace |
---|
| 7 | % a = rand(500,500); % first image |
---|
| 8 | % b = rand(340,220); % second image |
---|
| 9 | % opt = detbestlength2(FFTrv,FFTiv,IFFTiv,size(a),size(b),isreal(a),isreal(b)); % optimized parameters |
---|
| 10 | % y0 = fftolamopt2(a,b,opt); % equivalent to y0 = conv2(a,b); |
---|
| 11 | % |
---|
| 12 | % INPUT |
---|
| 13 | % a: first image (2D double matrix) |
---|
| 14 | % b: second image (2D double matrix) |
---|
| 15 | % opt: the optimized parameters calculated by detbestlength.m function |
---|
| 16 | % opt = detbestlength(FFTrv,FFTiv,IFFTiv,size(a),size(b)); |
---|
| 17 | % shape: returns a subsection of the 2D convolution with size specified by |
---|
| 18 | % 'shape': |
---|
| 19 | % 'full' - (default) returns the full 2-D convolution, |
---|
| 20 | % 'same' - returns the central part of the convolution |
---|
| 21 | % that is the same size as A. |
---|
| 22 | % 'valid' - returns only those parts of the convolution |
---|
| 23 | % that are computed without the zero-padded |
---|
| 24 | % edges. size(C) = [ma-mb+1,na-nb+1] when |
---|
| 25 | % all(size(A) >= size(B)), otherwise C is empty. |
---|
| 26 | % See also conv2. |
---|
| 27 | % OUTPUT |
---|
| 28 | % out: 2D convolution of a and b matrices: out = conv2(a,b); |
---|
| 29 | |
---|
| 30 | |
---|
| 31 | % Original size |
---|
| 32 | [z1x,z1y] = size(a); |
---|
| 33 | [z2x,z2y] = size(b); |
---|
| 34 | |
---|
| 35 | % Reverse a and b if necessary |
---|
| 36 | if opt.inverse |
---|
| 37 | atemp = a; |
---|
| 38 | a = b; |
---|
| 39 | b = atemp; |
---|
| 40 | end |
---|
| 41 | |
---|
| 42 | fftorder = zeros(2,1); |
---|
| 43 | ifftorder = zeros(2,1); |
---|
| 44 | fftsize = zeros(2,1); |
---|
| 45 | filterord = zeros(2,1); |
---|
| 46 | filtersiz = zeros(2,1); |
---|
| 47 | |
---|
| 48 | if (opt.fftxfirst == 1) |
---|
| 49 | fftorder(1) = 1; |
---|
| 50 | fftorder(2) = 2; |
---|
| 51 | fftsize(1) = opt.nfftx; |
---|
| 52 | fftsize(2) = opt.nffty; |
---|
| 53 | else |
---|
| 54 | fftorder(1) = 2; |
---|
| 55 | fftorder(2) = 1; |
---|
| 56 | fftsize(1) = opt.nffty; |
---|
| 57 | fftsize(2) = opt.nfftx; |
---|
| 58 | end |
---|
| 59 | |
---|
| 60 | |
---|
| 61 | if (opt.ifftxfirst == 1) |
---|
| 62 | ifftorder(1) = 1; |
---|
| 63 | ifftorder(2) = 2; |
---|
| 64 | else |
---|
| 65 | ifftorder(1) = 2; |
---|
| 66 | ifftorder(2) = 1; |
---|
| 67 | end |
---|
| 68 | |
---|
| 69 | if opt.filterxfirst==1 |
---|
| 70 | filterord(1) = 1; |
---|
| 71 | filterord(2) = 2; |
---|
| 72 | |
---|
| 73 | filtersiz(1) = opt.nfftx; |
---|
| 74 | filtersiz(2) = opt.nffty; |
---|
| 75 | else |
---|
| 76 | filterord(1) = 2; |
---|
| 77 | filterord(2) = 1; |
---|
| 78 | |
---|
| 79 | filtersiz(1) = opt.nffty; |
---|
| 80 | filtersiz(2) = opt.nfftx; |
---|
| 81 | end |
---|
| 82 | |
---|
| 83 | siz1 = opt.nfftx; |
---|
| 84 | siz2 = opt.nffty; |
---|
| 85 | |
---|
| 86 | [ax,ay] = size(a); |
---|
| 87 | [bx,by] = size(b); |
---|
| 88 | dimx = ax+bx-1; |
---|
| 89 | dimy = ay+by-1; |
---|
| 90 | nfftx = siz1; |
---|
| 91 | nffty = siz2; |
---|
| 92 | Lx = nfftx-bx+1; |
---|
| 93 | Ly = nffty-by+1; |
---|
| 94 | B = fft(fft(b,filtersiz(1),filterord(1)),filtersiz(2),filterord(2)); |
---|
| 95 | out = zeros(dimx,dimy); |
---|
| 96 | x0 = 1; |
---|
| 97 | while x0 <= ax |
---|
| 98 | x1 = min(x0+Lx-1,ax); |
---|
| 99 | y0 = 1; |
---|
| 100 | endx = min(dimx,x0+nfftx-1); |
---|
| 101 | while y0 <= ay |
---|
| 102 | y1 = min(y0+Ly-1,ay); |
---|
| 103 | endy = min(dimy,y0+nffty-1); |
---|
| 104 | X = fft(fft(a(x0:x1,y0:y1),fftsize(1),fftorder(1)),fftsize(2),fftorder(2)); |
---|
| 105 | Y = ifft(ifft(X.*B,[],ifftorder(1)),[],ifftorder(2)); |
---|
| 106 | out(x0:endx,y0:endy) = out(x0:endx,y0:endy)+Y(1:(endx-x0+1),1:(endy-y0+1)); |
---|
| 107 | y0 = y0+Ly; |
---|
| 108 | end |
---|
| 109 | x0 = x0+Lx; |
---|
| 110 | end |
---|
| 111 | if isreal(a) && isreal(b) |
---|
| 112 | out=real(out); |
---|
| 113 | end |
---|
| 114 | if nargin<4 || strcmp(shape,'full') |
---|
| 115 | return; |
---|
| 116 | end |
---|
| 117 | if strcmp(shape,'valid') |
---|
| 118 | if ((z1x<z2x)||(z1y<z2y)) |
---|
| 119 | out = []; |
---|
| 120 | else |
---|
| 121 | px = z2x; |
---|
| 122 | py = z2y; |
---|
| 123 | out = out(px:px+z1x-z2x,py:py+z1y-z2y); |
---|
| 124 | end |
---|
| 125 | return; |
---|
| 126 | end |
---|
| 127 | if strcmp(shape,'same') |
---|
| 128 | px = ((z2x-1)+mod((z2x-1),2))/2; |
---|
| 129 | py = ((z2y-1)+mod((z2y-1),2))/2; |
---|
| 130 | out = out(px+1:px+z1x,py+1:py+z1y); |
---|
| 131 | return; |
---|
| 132 | end |
---|