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 |
---|