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