1 /*=============================================================================
2 File: rec_filter.cpp
3 Purpose:
4 Revision: $Id: rec_filter.cpp,v 1.2 2002/05/13 21:07:45 philosophil Exp $
5 Created by: Philippe Lavoie (18 February 1999)
6 Modified by:
7
8 Copyright notice:
9 Copyright (C) 1996-1997 Philippe Lavoie
10
11 This library is free software; you can redistribute it and/or
12 modify it under the terms of the GNU Library General Public
13 License as published by the Free Software Foundation; either
14 version 2 of the License, or (at your option) any later version.
15
16 This library is distributed in the hope that it will be useful,
17 but WITHOUT ANY WARRANTY; without even the implied warranty of
18 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
19 Library General Public License for more details.
20
21 You should have received a copy of the GNU Library General Public
22 License along with this library; if not, write to the Free
23 Software Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
24 =============================================================================*/
25
26 #ifndef REC_FILTER_SOURCE
27 #define REC_FILTER_SOURCE
28
29 #include "rec_filter.h"
30
31 /*!
32 */
33 namespace PLib{
34
toParams(double a1,double a2,double a3,double a4,double a5,double a6,double a7,double a8,double b1,double b2,double c1,double c2,double k,Params & params)35 void toParams(double a1, double a2, double a3, double a4, double a5, double a6, double a7, double a8, double b1, double b2, double c1, double c2, double k, Params& params){
36 params.a1 = a1 ;
37 params.a2 = a2 ;
38 params.a3 = a3 ;
39 params.a4 = a4 ;
40 params.a5 = a5 ;
41 params.a6 = a6 ;
42 params.a7 = a7 ;
43 params.a8 = a8 ;
44 params.b1 = b1 ;
45 params.b2 = b2 ;
46 params.c1 = c1 ;
47 params.c2 = c2 ;
48 params.k = k ;
49 }
50
fromParams(const Params & params,double & a1,double & a2,double & a3,double & a4,double & a5,double & a6,double & a7,double & a8,double & b1,double & b2,double & c1,double & c2,double & k)51 void fromParams(const Params& params, double &a1, double &a2, double &a3, double &a4, double &a5, double &a6, double &a7, double &a8, double &b1, double &b2, double &c1, double &c2, double &k){
52 a1 = params.a1 ;
53 a2 = params.a2 ;
54 a3 = params.a3 ;
55 a4 = params.a4 ;
56 a5 = params.a5 ;
57 a6 = params.a6 ;
58 a7 = params.a7 ;
59 a8 = params.a8 ;
60 b1 = params.b1 ;
61 b2 = params.b2 ;
62 c1 = params.c1 ;
63 c2 = params.c2 ;
64 k = params.k ;
65 }
66
generalRFx(const Params & params,const Basic2DArray<double> & x,Basic2DArray<double> & y)67 void generalRFx(const Params& params, const Basic2DArray<double> &x, Basic2DArray<double> &y){
68 Basic2DArray<double> y1,y2 ;
69 y1.resize(x);
70 y2.resize(x);
71 y.resize(x);
72
73 double a1,a2,a3,a4,a5,a6,a7,a8 ;
74 double b1,b2 ;
75 double c1,c2 ;
76 double k ;
77
78 fromParams(params,a1,a2,a3,a4,a5,a6,a7,a8,b1,b2,c1,c2,k) ;
79 int i,j ;
80 const int rows = x.rows() ;
81 const int cols = x.cols() ;
82
83 for(i=0;i<rows-1;++i){
84 y1(i,0) = a1*x(i,0) ;
85 y1(i,1) = a1*x(i,1) + a2*x(i,0) + b1*y1(i,0) ;
86 for(j=2;j<cols-1;++j){
87 y1(i,j) = a1*x(i,j) + a2*x(i,j-1) + b1*y1(i,j-1) + b2*y1(i,j-2) ;
88 }
89 }
90
91 for(i=rows-1;i>=0;--i){
92 j = cols-1 ;
93 y2(i,j) = 0 ;
94 y(i,j) = c1*(y1(i,j)+y2(i,j));
95 j = cols-2 ;
96 y2(i,j) = a3*x(i,j+1) + b1*y2(i,j+1) ;
97 y(i,j) = c1*(y1(i,j)+y2(i,j));
98 for(j=cols-3;j>=0;--j){
99 y2(i,j) = a3*x(i,j+1) + a4*x(i,j+2) + b1*y2(i,j+1) + b2*y2(i,j+2) ;
100 y(i,j) = c1*(y1(i,j)+y2(i,j));
101 }
102 }
103
104 }
105
106
generalRFy(const Params & params,const Basic2DArray<double> & r,Basic2DArray<double> & y)107 void generalRFy(const Params& params, const Basic2DArray<double> &r, Basic2DArray<double> &y){
108 Basic2DArray<double> y1,y2 ;
109 y1.resize(r);
110 y2.resize(r);
111 y.resize(r);
112
113 double a1,a2,a3,a4,a5,a6,a7,a8 ;
114 double b1,b2 ;
115 double c1,c2 ;
116 double k ;
117
118 fromParams(params,a1,a2,a3,a4,a5,a6,a7,a8,b1,b2,c1,c2,k) ;
119 int i,j ;
120 const int rows = r.rows() ;
121 const int cols = r.cols() ;
122
123 for(j=0;j<cols-1;++j){
124 y1(0,j) = a5*r(0,j) ;
125 y1(1,j) = a5*r(1,j) + a6*r(0,j) + b1*y1(0,j);
126 for(i=2;i<rows-1;++i){
127 y1(i,j) = a5*r(i,j) + a6*r(i-1,j) + b1*y1(i-1,j) + b2*y1(i-2,j) ;
128 }
129 }
130
131 for(j=cols-1;j>=0;--j){
132 i = rows-1 ;
133 y2(i,j) = 0 ;
134 y(i,j) = c2*(y1(i,j)+y2(i,j));
135 i = rows-2 ;
136 y2(i,j) = a7*r(i+1,j) + b1*y2(i+1,j) ;
137 y(i,j) = c2*(y1(i,j)+y2(i,j));
138 for(i=rows-3;i>=0;--i){
139 y2(i,j) = a7*r(i+1,j) + a8*r(i+2,j) + b1*y2(i+1,j) + b2*y2(i+2,j) ;
140 y(i,j) = c2*(y1(i,j)+y2(i,j));
141 }
142 }
143 }
144
145
146
generalRF(const Params & params,const Basic2DArray<double> & x,Basic2DArray<double> & y)147 void generalRF(const Params& params, const Basic2DArray<double> &x, Basic2DArray<double> &y){
148 Basic2DArray<double> y1,y2,r ;
149 y1.resize(x);
150 y2.resize(x);
151 r.resize(x) ;
152 y.resize(x);
153
154 double a1,a2,a3,a4,a5,a6,a7,a8 ;
155 double b1,b2 ;
156 double c1,c2 ;
157 double k ;
158
159 fromParams(params,a1,a2,a3,a4,a5,a6,a7,a8,b1,b2,c1,c2,k) ;
160
161 int i,j ;
162 const int rows = x.rows() ;
163 const int cols = x.cols() ;
164
165 for(i=0;i<rows-1;++i){
166 y1(i,0) = a1*x(i,0) ;
167 y1(i,1) = a1*x(i,1) + a2*x(i,0) + b1*y1(i,0) ;
168 for(j=2;j<cols-1;++j){
169 y1(i,j) = a1*x(i,j) + a2*x(i,j-1) + b1*y1(i,j-1) + b2*y1(i,j-2) ;
170 }
171 }
172
173 for(i=rows-1;i>=0;--i){
174 j = cols-1 ;
175 y2(i,j) = 0 ;
176 r(i,j) = c1*(y1(i,j)+y2(i,j));
177 j = cols-2 ;
178 y2(i,j) = a3*x(i,j+1) + b1*y2(i,j+1) ;
179 r(i,j) = c1*(y1(i,j)+y2(i,j));
180 for(j=cols-3;j>=0;--j){
181 y2(i,j) = a3*x(i,j+1) + a4*x(i,j+2) + b1*y2(i,j+1) + b2*y2(i,j+2) ;
182 r(i,j) = c1*(y1(i,j)+y2(i,j));
183 }
184 }
185
186 for(j=0;j<cols-1;++j){
187 y1(0,j) = a5*r(0,j) ;
188 y1(1,j) = a5*r(1,j) + a6*r(0,j) + b1*y1(0,j);
189 for(i=2;i<rows-1;++i){
190 y1(i,j) = a5*r(i,j) + a6*r(i-1,j) + b1*y1(i-1,j) + b2*y1(i-2,j) ;
191 }
192 }
193
194 for(j=cols-1;j>=0;--j){
195 i = rows-1 ;
196 y2(i,j) = 0 ;
197 y(i,j) = c2*(y1(i,j)+y2(i,j));
198 i = rows-2 ;
199 y2(i,j) = a7*r(i+1,j) + b1*y2(i+1,j) ;
200 y(i,j) = c2*(y1(i,j)+y2(i,j));
201 for(i=rows-3;i>=0;--i){
202 y2(i,j) = a7*r(i+1,j) + a8*r(i+2,j) + b1*y2(i+1,j) + b2*y2(i+2,j) ;
203 y(i,j) = c2*(y1(i,j)+y2(i,j));
204 }
205 }
206
207 }
208
smooth2ndOrder(Params & params,double alpha)209 void smooth2ndOrder(Params& params, double alpha){
210 const double expa = exp(-alpha) ;
211 const double expa2 = exp(-2*alpha) ;
212 params.k = (1-expa)*(1-expa)/(1+2*alpha*expa-expa2) ;
213 params.b1 = 2*expa ;
214 params.b2 = -expa2 ;
215 params.a1 = params.a5 = params.k ;
216 params.a2 = params.a6 = params.k*expa*(alpha-1) ;
217 params.a3 = params.a7 = params.k*expa*(alpha+1) ;
218 params.a4 = params.a8 = -params.k*expa2;
219 params.c1 = params.c2 = 1 ;
220 }
221
xderiv2ndOrderSmooth(Params & params,double alpha)222 void xderiv2ndOrderSmooth(Params& params, double alpha){
223 const double expa = exp(-alpha) ;
224 const double expa2 = exp(-2*alpha) ;
225 params.k = (1-expa)*(1-expa)/(1+2*alpha*expa-expa2) ;
226 params.b1 = 2*expa ;
227 params.b2 = -expa2 ;
228 params.a1 = 0 ;
229 params.a2 = 1 ;
230 params.a3 = -1 ;
231 params.a4 = 0 ;
232 params.a5 = params.k ;
233 params.a6 = params.k*expa*(alpha-1) ;
234 params.a7 = params.k*expa*(alpha+1) ;
235 params.a8 = -params.k*expa2 ;
236 params.c1 = -1*(1-expa)*(1-expa);
237 params.c2 = 1 ;
238 }
239
yderiv2ndOrderSmooth(Params & params,double alpha)240 void yderiv2ndOrderSmooth(Params& params, double alpha){
241 const double expa = exp(-alpha) ;
242 const double expa2 = exp(-2*alpha) ;
243 params.k = (1-expa)*(1-expa)/(1+2*alpha*expa-expa2) ;
244 params.b1 = 2*expa ;
245 params.b2 = -expa2 ;
246 params.a1 = params.k ;
247 params.a2 = params.k*expa*(alpha-1) ;
248 params.a3 = params.k*expa*(alpha+1) ;
249 params.a4 = -params.k*expa2 ;
250 params.a5 = 0 ;
251 params.a6 = 1 ;
252 params.a7 = -1 ;
253 params.a8 = 0 ;
254 params.c1 = 1 ;
255 params.c2 = -1*(1-expa)*(1-expa);
256 }
257
xderiv2ndOrder(Params & params,double alpha)258 void xderiv2ndOrder(Params& params, double alpha){
259 const double expa = exp(-alpha) ;
260 const double expa2 = exp(-2*alpha) ;
261 params.k = (1-expa)*(1-expa)/(expa) ;
262 params.b1 = 2*expa ;
263 params.b2 = -2*expa2 ;
264 params.a1 = 0 ;
265 params.a2 = 1 ;
266 params.a3 = -1 ;
267 params.a4 = 0 ;
268 params.c1 = params.k*expa ;
269 params.c2 = params.k*expa ;
270 }
271
yderiv2ndOrder(Params & params,double alpha)272 void yderiv2ndOrder(Params& params, double alpha){
273 const double expa = exp(-alpha) ;
274 const double expa2 = exp(-2*alpha) ;
275 params.k = (1-expa)*(1-expa)/(expa) ;
276 params.b1 = 2*expa ;
277 params.b2 = -2*expa2 ;
278 params.a5 = 0 ;
279 params.a6 = 1 ;
280 params.a7 = -1 ;
281 params.a8 = 0 ;
282 params.c1 = params.k*expa ;
283 params.c2 = params.k*expa ;
284 }
285
smooth1stOrder(Params & params,double alpha,double k0)286 void smooth1stOrder(Params& params, double alpha, double k0){
287 const double expa = exp(-alpha) ;
288 params.k = (1-expa)/(1+expa) ; //(1-exp(-2*alpha))/(2*alpha*expa) ;
289 params.b1 = expa ;
290 params.b2 = 0;
291 params.a1 = 1 ;
292 params.a2 = 0 ;
293 params.a3 = expa ;
294 params.a4 = 0 ;
295 params.a5 = 1 ;
296 params.a6 = 0 ;
297 params.a7 = expa ;
298 params.a8 = 0 ;
299 params.c1 = params.k ;
300 params.c2 = params.k ;
301 }
302
LL1stOrder(Params & params,double alpha)303 void LL1stOrder(Params& params, double alpha){
304 const double expa = exp(-alpha) ;
305 const double expa2 = exp(-2*alpha);
306 params.k = (1-expa2)/(2*alpha*expa) ; //(1-expa)/(1+expa) ;
307 params.b1 = 2*expa ;
308 params.b2 = -expa2;
309 params.a1 = 0 ;
310 params.a2 = 1 ;
311 params.a3 = 1 ;
312 params.a4 = 0 ;
313 params.a5 = 0 ;
314 params.a6 = 1 ;
315 params.a7 = 1 ;
316 params.a8 = 0 ;
317 params.c1 = params.c2 = (1-expa2)/2 ;
318 }
319
ubyteToDouble(const Basic2DArray<unsigned char> & img,Basic2DArray<double> & mat)320 void ubyteToDouble(const Basic2DArray<unsigned char>& img, Basic2DArray<double>& mat){
321 mat.resize(img.rows(),img.cols());
322 for(int i=img.rows()-1;i>=0;--i)
323 for(int j=img.cols()-1;j>=0;--j)
324 mat(i,j) = img(i,j) ;
325 }
326
doubleToUbyte(const Basic2DArray<double> & mat,Basic2DArray<unsigned char> & img)327 void doubleToUbyte(const Basic2DArray<double>& mat, Basic2DArray<unsigned char>& img){
328 img.resize(mat.rows(),mat.cols());
329 for(int i=img.rows()-1;i>=0;--i)
330 for(int j=img.cols()-1;j>=0;--j)
331 img(i,j) = (unsigned char)mat(i,j) ;
332 }
333
quadInterp(double x,double x0,double f0,double x1,double f1,double x2,double f2)334 double quadInterp(double x, double x0, double f0, double x1, double f1, double x2, double f2){
335 const double f01 = firstDivDiff(x0,f0,x1,f1) ;
336 const double f012 = secondDivDiff(x0,f0,x1,f1,x2,f2);
337 double fx = f0 + (x-x0)*f01 + (x-x0)*(x-x1)*f012 ;
338 return fx;
339 }
340
findEdge(const Basic2DArray<double> & dx,const Basic2DArray<double> & dy,Basic2DArray<double> & edge,Basic2DArray<double> & gradN,double thresh)341 int findEdge(const Basic2DArray<double>& dx, const Basic2DArray<double>& dy, Basic2DArray<double> &edge, Basic2DArray<double>& gradN, double thresh){
342 if(dx.rows() != dy.rows() || dx.cols() != dy.cols())
343 return 0;
344 edge.resize(dx) ;
345 gradN.resize(dx);
346 //Basic2DArray<double> gradN(dx.rows(),dx.cols());
347 int i,j ;
348
349 for(i=0;i<dx.rows();++i)
350 for(j=0;j<dx.cols();++j){
351 gradN(i,j) = sqrt(dx(i,j)*dx(i,j) + dy(i,j)*dy(i,j)) ;
352 }
353
354
355 for(i=1;i<dx.rows()-1;++i)
356 for(j=1;j<dx.cols()-1;++j){
357 if(absolute(dx(i,j)) > absolute(dy(i,j))){
358 double d = 1/dx(i,j) ;
359 double y = dy(i,j)/dx(i,j) ;
360 double b = quadInterp(i+y,i-1,gradN(i-1,j+1),i,gradN(i,j+1),i+1,gradN(i+1,j+1));
361 double c = quadInterp(i-y,i-1,gradN(i-1,j-1),i,gradN(i,j-1),i+1,gradN(i+1,j-1));
362 if(gradN(i,j) >= b && gradN(i,j) >= c && gradN(i,j) > thresh)
363 edge(i,j) = 200;
364 else
365 edge(i,j) = 0 ;
366
367 }
368 else{
369 double d = 1/dy(i,j) ;
370 double x = dx(i,j)/dy(i,j) ;
371 double b = quadInterp(j-x,j-1,gradN(i-1,j-1),j,gradN(i-1,j),j+1,gradN(i-1,j+1));
372 double c = quadInterp(j+x,j-1,gradN(i+1,j-1),j,gradN(i+1,j),j+1,gradN(i+1,j+1));
373 if(gradN(i,j) >= b && gradN(i,j) >= c && gradN(i,j) > thresh)
374 edge(i,j) = 200;
375 else
376 edge(i,j) = 0 ;
377 }
378 }
379 return 1;
380 }
381
findSubEdge(const Basic2DArray<double> & dx,const Basic2DArray<double> & dy,Basic2DArray<double> & edge,double thresh,const Basic2DArray<double> & in)382 int findSubEdge(const Basic2DArray<double>& dx, const Basic2DArray<double>& dy, Basic2DArray<double> &edge, double thresh, const Basic2DArray<double>& in){
383 if(dx.rows() != dy.rows() || dx.cols() != dy.cols())
384 return 0;
385 const int zoom = 7 ;
386 const int zoom2 = 4 ;
387 edge.resize(dx.rows()*zoom,dx.cols()*zoom) ;
388 Basic2DArray<double> gradN(dx.rows(),dx.cols());
389 int i,j ;
390
391 for(i=0;i<dx.rows();++i)
392 for(j=0;j<dx.cols();++j){
393 gradN(i,j) = sqrt(dx(i,j)*dx(i,j) + dy(i,j)*dy(i,j)) ;
394 }
395
396
397 for(i=1;i<dx.rows()-1;++i)
398 for(j=1;j<dx.cols()-1;++j)
399 for(int k=0;k<zoom;++k)
400 for(int l=0;l<zoom;l++)
401 edge(i*zoom+k-zoom2,j*zoom+l-zoom2) = (in(i,j)>250) ? 250 : in(i,j);
402
403
404 for(i=1;i<dx.rows()-1;++i)
405 for(j=1;j<dx.cols()-1;++j){
406 if(absolute(dx(i,j)) > absolute(dy(i,j))){
407 double d = 1/dx(i,j) ;
408 double y = dy(i,j)/dx(i,j) ;
409 double c = quadInterp(i+y,i-1,gradN(i-1,j+1),i,gradN(i,j+1),i+1,gradN(i+1,j+1));
410 double a = quadInterp(i-y,i-1,gradN(i-1,j-1),i,gradN(i,j-1),i+1,gradN(i+1,j-1));
411 if(gradN(i,j) >= a && gradN(i,j) >= c && gradN(i,j) > thresh){
412 double m = (a-c)/( 2*(a-2*gradN(i,j)+c)) ;
413 double di, dj ;
414 if(absolute(m)>0.5) cerr << " m = " << m << " " << a << " " << gradN(i,j) << " " << c << " at " << i << j << endl ;
415 dj = zoom*m*absolute(dx(i,j))/gradN(i,j) ;
416 di = zoom*m*absolute(dy(i,j))/gradN(i,j) ;
417 if(m<0 && dx(i,j)*dy(i,j)<0)
418 di *= -1 ;
419 if(m>0 && dx(i,j)*dy(i,j)<0)
420 di *= -1 ;
421 for(int k=0;k<zoom;++k)
422 for(int l=0;l<zoom;l++)
423 edge(i*zoom+k-zoom2,j*zoom+l-zoom2) = 255;
424 edge(i*zoom+di,j*zoom+dj) = 254 ;
425 }
426 }
427 else{
428 double d = 1/dy(i,j) ;
429 double x = dx(i,j)/dy(i,j) ;
430 double a = quadInterp(j-x,j-1,gradN(i-1,j-1),j,gradN(i-1,j),j+1,gradN(i-1,j+1));
431 double c = quadInterp(j+x,j-1,gradN(i+1,j-1),j,gradN(i+1,j),j+1,gradN(i+1,j+1));
432 if(gradN(i,j) >= a && gradN(i,j) >= c && gradN(i,j) > thresh){
433 double m = (a-c)/( 2*(a-2*gradN(i,j)+c)) ;
434 double di, dj ;
435 if(absolute(m)>0.5) cerr << " m = " << m << " " << a << " " << gradN(i,j) << " " << c << " at " << i << j << endl ;
436 dj = zoom*m*absolute(dx(i,j))/gradN(i,j);
437 if(m<0 && dx(i,j)*dy(i,j)<0)
438 dj *= -1 ;
439 if(m>0 && dx(i,j)*dy(i,j)<0)
440 dj *= -1 ;
441 di = zoom*m*absolute(dy(i,j))/gradN(i,j) ;
442 for(int k=0;k<zoom;++k)
443 for(int l=0;l<zoom;l++)
444 edge(i*zoom+k-zoom2,j*zoom+l-zoom2) = 255 ;
445 edge(i*zoom+di,j*zoom+dj) = 254 ;
446 }
447 }
448 }
449 return 1;
450 }
451
452
453 /*
454 int main(){
455 IM_Image img;
456 if(!img.read("test.png")){
457 cerr << "Can't open file test.png\n" ;
458 }
459 Matrix<double> in,dx,dy,e,e2,g ;
460 ubyteToDouble(img,in);
461 Params params ;
462
463 const double alpha = 1.0 ;
464
465 xderiv2ndOrder(params,alpha);
466 generalRFx(params,in,dx);
467
468 yderiv2ndOrder(params,alpha);
469 generalRFy(params,in,dy);
470 cout << "Testing the quad interpolation\n" ;
471
472 cout << quadInterp(2.5,0,10,5,5,10,0) << endl ;
473 cout << quadInterp(-0.5,-1,4,0,5,1,1) << endl ;
474 cout << quadInterp(-0.4,-1,4,0,5,1,1) << endl ;
475 cout << quadInterp(-0.3,-1,4,0,5,1,1) << endl ;
476 cout << quadInterp(-0.2,-1,4,0,5,1,1) << endl ;
477 cout << quadInterp(9.2,8.0,2.0794,9.0,2.1972,9.5,2.2513) << endl ;
478
479 findEdge(dx,dy,e,g,20);
480 findSubEdge(dx,dy,e2,20,in) ;
481
482 dx += 128 ;
483 doubleToUbyte(dx,img);
484 img.write("outputA.png");
485
486 dy += 128 ;
487 doubleToUbyte(dy,img);
488 img.write("outputB.png");
489
490 doubleToUbyte(e,img);
491 img.write("edge.png");
492
493 doubleToUbyte(e2,img);
494 img.write("edge2.png");
495
496 doubleToUbyte(g,img);
497 img.write("grad.png");
498
499
500
501 }
502 */
503
504 template <class T>
RecursiveFilter(const Basic2DArray<T> & in,Basic2DArray<T> & out)505 RecursiveFilter<T>::RecursiveFilter(const Basic2DArray<T>& in, Basic2DArray<T>& out): input(in), output(out) {
506 input_ = new Basic2DArray<double>(in.rows(),in.cols()) ;
507 output_ = new Basic2DArray<double>(out.rows(),out.cols()) ;
508 output.resize(out.rows(),out.cols());
509 toDouble(input,*input_);
510 }
511
512 template <>
RecursiveFilter(const Basic2DArray<double> & in,Basic2DArray<double> & out)513 RecursiveFilter<double>::RecursiveFilter(const Basic2DArray<double>& in, Basic2DArray<double>& out): input(in), output(out) {
514 input_ = const_cast<Basic2DArray<double>*>(&in) ;
515 //input_ = &in ;
516 output_ = &out ;
517 output.resize(out.rows(),out.cols());
518 }
519
520 template <class T>
compute_xderiv2ndOrderSmooth(double alpha)521 void RecursiveFilter<T>::compute_xderiv2ndOrderSmooth(double alpha){
522 xderiv2ndOrderSmooth(params,alpha) ;
523 generalRFx(params,*input_,*output_);
524 cerr << "size = " << output_->rows() << "," << output_->cols() << endl ;
525 fromDouble(*output_,output) ;
526 }
527
528 template <class T>
compute_yderiv2ndOrderSmooth(double alpha)529 void RecursiveFilter<T>::compute_yderiv2ndOrderSmooth(double alpha){
530 yderiv2ndOrderSmooth(params,alpha);
531 generalRFy(params,*input_,*output_);
532 fromDouble(*output_,output) ;
533 }
534
535 template <class T>
compute_xderiv2ndOrder(double alpha)536 void RecursiveFilter<T>::compute_xderiv2ndOrder(double alpha){
537 xderiv2ndOrder(params,alpha) ;
538 generalRFx(params,*input_,*output_);
539 fromDouble(*output_,output) ;
540 }
541
542 template <class T>
compute_yderiv2ndOrder(double alpha)543 void RecursiveFilter<T>::compute_yderiv2ndOrder(double alpha){
544 yderiv2ndOrder(params,alpha);
545 generalRFy(params,*input_,*output_);
546 fromDouble(*output_,output) ;
547 }
548
549 template <class T>
compute_smooth1stOrder(double alpha,double k0)550 void RecursiveFilter<T>::compute_smooth1stOrder(double alpha, double k0){
551 smooth1stOrder(params,alpha,k0) ;
552 generalRF(params,*input_,*output_);
553 fromDouble(*output_,output);
554 }
555
556
557 template <class T>
compute_smooth1stOrder_x(double alpha,double k0)558 void RecursiveFilter<T>::compute_smooth1stOrder_x(double alpha, double k0){
559 smooth1stOrder(params,alpha,k0) ;
560 generalRFx(params,*input_,*output_);
561 fromDouble(*output_,output) ;
562 }
563
564
565 template <class T>
compute_smooth1stOrder_y(double alpha,double k0)566 void RecursiveFilter<T>::compute_smooth1stOrder_y(double alpha, double k0){
567 smooth1stOrder(params,alpha,k0) ;
568 generalRFy(params,*input_,*output_);
569 fromDouble(*output_,output) ;
570 }
571
572
573 template <class T>
compute_smooth2ndOrder(double alpha)574 void RecursiveFilter<T>::compute_smooth2ndOrder(double alpha){
575 smooth2ndOrder(params,alpha) ;
576 generalRF(params,*input_,*output_);
577 fromDouble(*output_,output) ;
578 }
579
580
581 template <class T>
compute_smooth2ndOrder_x(double alpha)582 void RecursiveFilter<T>::compute_smooth2ndOrder_x(double alpha){
583 smooth2ndOrder(params,alpha) ;
584 generalRFx(params,*input_,*output_);
585 fromDouble(*output_,output) ;
586 }
587
588
589 template <class T>
compute_smooth2ndOrder_y(double alpha)590 void RecursiveFilter<T>::compute_smooth2ndOrder_y(double alpha){
591 smooth2ndOrder(params,alpha) ;
592 generalRFy(params,*input_,*output_);
593 fromDouble(*output_,output) ;
594 }
595
596
597 template <class T>
compute_LL1stOrder(double alpha)598 void RecursiveFilter<T>::compute_LL1stOrder(double alpha){
599 LL1stOrder(params,alpha);
600 generalRF(params,*input_,*output_);
601 fromDouble(*output_,output) ;
602 }
603
604
605 template <class T>
compute_LL1stOrder_x(double alpha)606 void RecursiveFilter<T>::compute_LL1stOrder_x(double alpha){
607 LL1stOrder(params,alpha);
608 generalRFx(params,*input_,*output_);
609 fromDouble(*output_,output) ;
610 }
611
612
613 template <class T>
compute_LL1stOrder_y(double alpha)614 void RecursiveFilter<T>::compute_LL1stOrder_y(double alpha){
615 LL1stOrder(params,alpha);
616 generalRFy(params,*input_,*output_);
617 fromDouble(*output_,output) ;
618 }
619
620 }
621
622
623
624 #endif
625