00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026 #ifndef MLN_LINEAR_GAUSSIAN_DIRECTIONAL_2D_HH
00027 # define MLN_LINEAR_GAUSSIAN_DIRECTIONAL_2D_HH
00028
00034
00035 #include <mln/core/image/image2d.hh>
00036 #include <mln/extension/adjust_fill.hh>
00037 #include <mln/extension/adjust_duplicate.hh>
00038
00039
00040
00041 namespace mln
00042 {
00043
00044 namespace linear
00045 {
00046
00047
00048 template <typename I>
00049 mln_concrete(I)
00050 gaussian_directional_2d(const Image<I>& input,
00051 unsigned dir, double sigma,
00052 const mln_value(I)& bdr);
00053
00054
00055
00056 # ifndef MLN_INCLUDE_ONLY
00057
00058 namespace my
00059 {
00060
00061 struct recursivefilter_coef_
00062 {
00063 enum FilterType { DericheGaussian,
00064 DericheGaussianFirstDerivative,
00065 DericheGaussianSecondDerivative };
00066
00067 std::vector<double> n, d, nm, dm;
00068 double sumA, sumC;
00069
00070 recursivefilter_coef_(double a0, double a1,
00071 double b0, double b1,
00072 double c0, double c1,
00073 double w0, double w1,
00074 double s, FilterType filter_type)
00075 {
00076 n.reserve(5);
00077 d.reserve(5);
00078 nm.reserve(5);
00079 dm.reserve(5);
00080
00081 b0 /= s;
00082 b1 /= s;
00083 w0 /= s;
00084 w1 /= s;
00085
00086 double sin0 = sin(w0);
00087 double sin1 = sin(w1);
00088 double cos0 = cos(w0);
00089 double cos1 = cos(w1);
00090
00091 switch (filter_type) {
00092
00093 case DericheGaussian :
00094 {
00095 sumA =
00096 (2.0 * a1 * exp( b0 ) * cos0 * cos0 - a0 * sin0 * exp( 2.0 * b0 )
00097 + a0 * sin0 - 2.0 * a1 * exp( b0 )) /
00098 (( 2.0 * cos0 * exp( b0 ) - exp( 2.0 * b0 ) - 1 ) * sin0);
00099
00100 sumC =
00101 (2.0 * c1 * exp( b1 ) * cos1 * cos1 - c0 * sin1 * exp( 2.0 * b1 )
00102 + c0 * sin1 - 2.0 * c1 * exp( b1 ))
00103 / (( 2.0 * cos1 * exp( b1 ) - exp( 2.0 * b1 ) - 1 ) * sin1);
00104 break;
00105 }
00106
00107 case DericheGaussianFirstDerivative :
00108 {
00109 sumA = -2.f *
00110 (a0 * cos0 - a1 * sin0 + a1 * sin0 * exp( 2.0 * b0 )
00111 + a0 * cos0 * exp( 2.0 * b0 ) - 2.0 * a0 * exp( b0 ))
00112 * exp( b0 )
00113 / (exp( 4.0 * b0 ) - 4.0 * cos0 * exp( 3.0 * b0 )
00114 + 2.0 * exp( 2.0 * b0 ) + 4.0 * cos0 * cos0 * exp( 2.0 * b0 )
00115 + 1.0 - 4.0 * cos0 * exp( b0 ));
00116 sumC = - 2.f *
00117 (c0 * cos1 - c1 * sin1 + c1 * sin1 * exp( 2.0 * b1 )
00118 + c0 * cos1 * exp( 2.0 * b1 ) - 2.0 * c0 * exp( b1 ))
00119 * exp( b1 ) /
00120 (exp( 4.0 * b1 ) - 4.0 * cos1 * exp( 3.0 * b1 )
00121 + 2.0 * exp( 2.0 * b1 ) + 4.0 * cos1 * cos1 * exp( 2.0 * b1 )
00122 + 1.0 - 4.0 * cos1 * exp( b1 ));
00123 break;
00124 }
00125
00126 case DericheGaussianSecondDerivative :
00127 {
00128 double aux;
00129 aux =
00130 12.0 * cos0 * exp( 3.0 * b0 ) - 3.0 * exp( 2.0 * b0 )
00131 + 8.0 * cos0 * cos0 * cos0 * exp( 3.0 * b0 ) - 12.0 * cos0 * cos0 *
00132 exp( 4.0 * b0 )
00133 - (3.0 * exp( 4.0 * b0 ))
00134 + 6.0 * cos0 * exp( 5.0 * b0 ) - exp( 6.0 * b0 ) + 6.0 * cos0 * exp
00135 ( b0 )
00136 - ( 1.0 + 12.0 * cos0 * cos0 * exp( 2.0 * b0 ) );
00137 sumA =
00138 4.0 * a0 * sin0 * exp( 3.0 * b0 ) + a1 * cos0 * cos0 * exp( 4.0 * b0 )
00139 - ( 4.0 * a0 * sin0 * exp( b0 ) + 6.0 * a1 * cos0 * cos0 * exp( 2.0 * b0 ) )
00140 + 2.0 * a1 * cos0 * cos0 * cos0 * exp( b0 ) - 2.0 * a1 * cos0 * exp(b0)
00141 + 2.0 * a1 * cos0 * cos0 * cos0 * exp( 3.0 * b0 ) - 2.0 * a1 * cos0
00142 * exp( 3.0 * b0 )
00143 + a1 * cos0 * cos0 - a1 * exp( 4.0 * b0 )
00144 + 2.0 * a0 * sin0 * cos0 * cos0 * exp( b0 ) - 2.0 * a0 * sin0 * cos0
00145 * cos0 * exp( 3.0 * b0 )
00146 - ( a0 * sin0 * cos0 * exp( 4.0 * b0 ) + a1 )
00147 + 6.0 * a1 * exp( 2.0 * b0 ) + a0 * cos0 * sin0
00148 * 2.0 * exp( b0 ) / ( aux * sin0 );
00149 aux =
00150 12.0 * cos1 * exp( 3.0 * b1 ) - 3.0 * exp( 2.0 * b1 )
00151 + 8.0 * cos1 * cos1 * cos1 * exp( 3.0 * b1 ) - 12.0 * cos1 * cos1 *
00152 exp( 4.0 * b1 )
00153 - 3.0 * exp( 4.0 * b1 )
00154 + 6.0 * cos1 * exp( 5.0 * b1 ) - exp( 6.0 * b1 ) + 6.0 * cos1 * exp
00155 ( b1 )
00156 - ( 1.0 + 12.0 * cos1 * cos1 * exp( 2.0 * b1 ) );
00157 sumC = 4.0 * c0 * sin1 * exp( 3.0 * b1 ) + c1 * cos1 * cos1 * exp( 4.0 * b1 )
00158 - ( 4.0 * c0 * sin1 * exp( b1 ) + 6.0 * c1 * cos1 * cos1 * exp( 2.0 * b1 ) )
00159 + 2.0 * c1 * cos1 * cos1 * cos1 * exp( b1 ) - 2.0 * c1 * cos1 * exp( b1 )
00160 + 2.0 * c1 * cos1 * cos1 * cos1 * exp( 3.0 * b1 ) - 2.0 * c1 * cos1
00161 * exp( 3.0 * b1 )
00162 + c1 * cos1 * cos1 - c1 * exp( 4.0 * b1 )
00163 + 2.0 * c0 * sin1 * cos1 * cos1 * exp( b1 ) - 2.0 * c0 * sin1 * cos1
00164 * cos1 * exp( 3.0 * b1 )
00165 - ( c0 * sin1 * cos1 * exp( 4.0 * b1 ) + c1 )
00166 + 6.0 * c1 * exp( 2.0 * b1 ) + c0 * cos1 * sin1
00167 * 2.0 * exp( b1 ) / ( aux * sin1 );
00168 sumA /= 2;
00169 sumC /= 2;
00170 break;
00171 }
00172 }
00173
00174 a0 /= (sumA + sumC);
00175 a1 /= (sumA + sumC);
00176 c0 /= (sumA + sumC);
00177 c1 /= (sumA + sumC);
00178
00179 n[3] =
00180 exp( -b1 - 2*b0 ) * (c1 * sin1 - cos1 * c0) +
00181 exp( -b0 - 2*b1 ) * (a1 * sin0 - cos0 * a0);
00182 n[2] =
00183 2 * exp(-b0 - b1) * ((a0 + c0) * cos1 * cos0 -
00184 cos1 * a1 * sin0 -
00185 cos0 * c1 * sin1) +
00186 c0 * exp(-2 * b0) + a0 * exp(-2 * b1);
00187 n[1] =
00188 exp(-b1) * (c1 * sin1 - (c0 + 2*a0) * cos1) +
00189 exp(-b0) * (a1 * sin0 - (2*c0 + a0) * cos0);
00190 n[0] =
00191 a0 + c0;
00192
00193 d[4] = exp(-2 * b0 - 2 * b1);
00194 d[3] =
00195 -2 * cos0 * exp(-b0 - 2*b1) -
00196 2 * cos1 * exp(-b1 - 2*b0);
00197 d[2] =
00198 4 * cos1 * cos0 * exp(-b0 - b1) +
00199 exp(-2*b1) + exp(-2*b0);
00200 d[1] =
00201 -2*exp(-b1) * cos1 - 2 * exp(-b0) * cos0;
00202
00203 switch (filter_type) {
00204 case DericheGaussian :
00205 case DericheGaussianSecondDerivative :
00206 {
00207 for (unsigned i = 1; i <= 3; ++i)
00208 {
00209 dm[i] = d[i];
00210 nm[i] = n[i] - d[i] * n[0];
00211 }
00212 dm[4] = d[4];
00213 nm[4] = -d[4] * n[0];
00214 break;
00215 }
00216 case DericheGaussianFirstDerivative :
00217 {
00218 for (unsigned i = 1; i <= 3; ++i)
00219 {
00220 dm[i] = d[i];
00221 nm[i] = - (n[i] - d[i] * n[0]);
00222 }
00223 dm[4] = d[4];
00224 nm[4] = d[4] * n[0];
00225 break;
00226 }
00227 }
00228 }
00229 };
00230
00231
00232 }
00233
00234
00235
00236
00237
00238
00239
00240
00241 template <typename I, typename C>
00242 inline
00243 void
00244 recursivefilter_directional_generic(I& ima,
00245 const C& c,
00246 const mln_psite(I)& start,
00247 const mln_psite(I)& finish,
00248 int len,
00249 const mln_deduce(I, psite, delta)& d)
00250 {
00251 std::vector<double> tmp1(len);
00252 std::vector<double> tmp2(len);
00253
00254 tmp1[0] =
00255 c.n[0]*ima(start);
00256
00257 tmp1[1] =
00258 c.n[0]*ima(start + d)
00259 + c.n[1]*ima(start)
00260 - c.d[1]*tmp1[0];
00261
00262 tmp1[2] =
00263 c.n[0]*ima(start + d + d)
00264 + c.n[1]*ima(start + d)
00265 + c.n[2]*ima(start)
00266 - c.d[1]*tmp1[1]
00267 - c.d[2]*tmp1[0];
00268
00269 tmp1[3] =
00270 c.n[0]*ima(start + d + d + d)
00271 + c.n[1]*ima(start + d + d)
00272 + c.n[2]*ima(start + d)
00273 + c.n[3]*ima(start)
00274 - c.d[1]*tmp1[2] - c.d[2]*tmp1[1]
00275 - c.d[3]*tmp1[0];
00276
00277 mln_site(I) current = start + d + d + d + d;
00278 for (int i = 4; i < len; ++i)
00279 {
00280 tmp1[i] =
00281 c.n[0]*ima(current)
00282 + c.n[1]*ima(current - d)
00283 + c.n[2]*ima(current - d - d)
00284 + c.n[3]*ima(current - d - d - d)
00285 - c.d[1]*tmp1[i - 1] - c.d[2]*tmp1[i - 2]
00286 - c.d[3]*tmp1[i - 3] - c.d[4]*tmp1[i - 4];
00287 current += d;
00288 }
00289
00290
00291
00292 tmp2[len - 1] = 0;
00293
00294 tmp2[len - 2] =
00295 c.nm[1]*ima(finish);
00296
00297 tmp2[len - 3] =
00298 c.nm[1]*ima(finish - d)
00299 + c.nm[2]*ima(finish)
00300 - c.dm[1]*tmp2[len-2];
00301
00302 tmp2[len - 4] =
00303 c.nm[1]*ima(finish - d - d)
00304 + c.nm[2]*ima(finish - d)
00305 + c.nm[3]*ima(finish)
00306 - c.dm[1]*tmp2[len-3]
00307 - c.dm[2]*tmp2[len-2];
00308
00309 current = finish - d - d - d ;
00310
00311 for (int i = len - 5; i >= 0; --i)
00312 {
00313 tmp2[i] =
00314 c.nm[1]*ima(current)
00315 + c.nm[2]*ima(current + d)
00316 + c.nm[3]*ima(current + d + d)
00317 + c.nm[4]*ima(current + d + d + d)
00318 - c.dm[1]*tmp2[i+1] - c.dm[2]*tmp2[i+2]
00319 - c.dm[3]*tmp2[i+3] - c.dm[4]*tmp2[i+4];
00320 current -= d;
00321 }
00322
00323
00324
00325 current = start;
00326 for (int i = 0; i < len; ++i)
00327 {
00328 ima(current) = (tmp1[i] + tmp2[i]);
00329 current += d;
00330 }
00331 }
00332
00333
00334
00335
00336 template <typename I, typename C>
00337 inline
00338 void
00339 recursivefilter_directional_fastest(I& ima,
00340 const C& c,
00341 const mln_psite(I)& start,
00342 const mln_psite(I)& finish,
00343 int len,
00344 const mln_deduce(I, psite, delta)& d,
00345 const mln_value(I)& bdr)
00346 {
00347
00348
00349
00350 std::vector<double> tmp1(len);
00351 std::vector<double> tmp2(len);
00352
00353 unsigned delta_offset = ima.delta_index(d);
00354 unsigned
00355 o_start = ima.index_of_point(start),
00356 o_start_d = o_start + delta_offset,
00357 o_start_dd = o_start + 2 * delta_offset,
00358 o_start_ddd = o_start + 3 * delta_offset;
00359
00360 tmp1[0] =
00361 c.n[0] * ima.element(o_start);
00362
00363 tmp1[1] = 0
00364 + c.n[0] * ima.element(o_start_d)
00365 + c.n[1] * ima.element(o_start)
00366 - c.d[1] * tmp1[0];
00367
00368 tmp1[2] = 0
00369 + c.n[0] * ima.element(o_start_dd)
00370 + c.n[1] * ima.element(o_start_d)
00371 + c.n[2] * ima.element(o_start)
00372 - c.d[1] * tmp1[1]
00373 - c.d[2] * tmp1[0];
00374
00375 tmp1[3] = 0
00376 + c.n[0] * ima.element(o_start_ddd)
00377 + c.n[1] * ima.element(o_start_dd)
00378 + c.n[2] * ima.element(o_start_d)
00379 + c.n[3] * ima.element(o_start)
00380 - c.d[1] * tmp1[2] - c.d[2] * tmp1[1]
00381 - c.d[3] * tmp1[0];
00382
00383 unsigned
00384 o_current = o_start + 4 * delta_offset,
00385 o_current_d = o_current - delta_offset,
00386 o_current_dd = o_current - 2 * delta_offset,
00387 o_current_ddd = o_current - 3 * delta_offset;
00388
00389 for (int i = 4; i < len; ++i)
00390 {
00391 tmp1[i] = 0
00392 + c.n[0] * ima.element(o_current)
00393 + c.n[1] * ima.element(o_current_d)
00394 + c.n[2] * ima.element(o_current_dd)
00395 + c.n[3] * ima.element(o_current_ddd)
00396 - c.d[1] * tmp1[i - 1] - c.d[2] * tmp1[i - 2]
00397 - c.d[3] * tmp1[i - 3] - c.d[4] * tmp1[i - 4];
00398 o_current += delta_offset;
00399 o_current_d += delta_offset;
00400 o_current_dd += delta_offset;
00401 o_current_ddd += delta_offset;
00402 }
00403
00404
00405
00406 (void) bdr;
00407
00408
00409 unsigned
00410 o_finish = ima.index_of_point(finish),
00411 o_finish_d = o_finish - delta_offset,
00412 o_finish_dd = o_finish - 2 * delta_offset;
00413
00414 tmp2[len - 1] = 0;
00415
00416 tmp2[len - 2] =
00417 c.nm[1] * ima.element(o_finish);
00418
00419 tmp2[len - 3] = 0
00420 + c.nm[1] * ima.element(o_finish_d)
00421 + c.nm[2] * ima.element(o_finish)
00422 - c.dm[1] * tmp2[len-2];
00423
00424 tmp2[len - 4] = 0
00425 + c.nm[1] * ima.element(o_finish_dd)
00426 + c.nm[2] * ima.element(o_finish_d)
00427 + c.nm[3] * ima.element(o_finish)
00428 - c.dm[1] * tmp2[len-3]
00429 - c.dm[2] * tmp2[len-2];
00430
00431 o_current = o_finish - 3 * delta_offset;
00432 o_current_d = o_current + delta_offset;
00433 o_current_dd = o_current + 2 * delta_offset;
00434 o_current_ddd = o_current + 3 * delta_offset;
00435
00436 for (int i = len - 5; i >= 0; --i)
00437 {
00438 tmp2[i] = 0
00439 + c.nm[1] * ima.element(o_current)
00440 + c.nm[2] * ima.element(o_current_d)
00441 + c.nm[3] * ima.element(o_current_dd)
00442 + c.nm[4] * ima.element(o_current_ddd)
00443 - c.dm[1] * tmp2[i+1] - c.dm[2] * tmp2[i+2]
00444 - c.dm[3] * tmp2[i+3] - c.dm[4] * tmp2[i+4];
00445 o_current -= delta_offset;
00446 o_current_d -= delta_offset;
00447 o_current_dd -= delta_offset;
00448 o_current_ddd -= delta_offset;
00449 }
00450
00451
00452
00453 o_current = o_start;
00454 for (int i = 0; i < len; ++i)
00455 {
00456 ima.element(o_current) = static_cast<mln_value(I)>(tmp1[i] + tmp2[i]);
00457 o_current += delta_offset;
00458 }
00459 }
00460
00461
00462 template <typename I>
00463 inline
00464 mln_concrete(I)
00465 gaussian_directional_2d(const Image<I>& input_,
00466 unsigned dir, double sigma,
00467 const mln_value(I)& bdr)
00468 {
00469 trace::entering("linear::gaussian_directional_2d");
00470
00471 typedef mln_site(I) P;
00472 mlc_bool(P::dim == 2)::check();
00473
00474 const I& input = exact(input_);
00475
00476 mln_precondition(dir == 0 || dir == 1);
00477 mln_precondition(input.is_valid());
00478
00479 my::recursivefilter_coef_ coef(1.68f, 3.735f,
00480 1.783f, 1.723f,
00481 -0.6803f, -0.2598f,
00482 0.6318f, 1.997f,
00483 sigma,
00484 my::recursivefilter_coef_::DericheGaussian);
00485
00486 extension::adjust_fill(input, 5 * int(sigma + .50001) + 1, bdr);
00487 mln_concrete(I) output = duplicate(input);
00488
00489 if (sigma < 0.006)
00490 return output;
00491
00492 int
00493 nrows = geom::nrows(input),
00494 ncols = geom::ncols(input),
00495 b = input.border();
00496
00497 if (dir == 0)
00498 {
00499 for (int j = 0; j < ncols; ++j)
00500 recursivefilter_directional_fastest(output, coef,
00501 point2d(- b, j),
00502 point2d(nrows - 1 + b, j),
00503 nrows + 2 * b,
00504 dpoint2d(1, 0),
00505 bdr);
00506 }
00507
00508 if (dir == 1)
00509 {
00510 for (int i = 0; i < nrows; ++i)
00511 recursivefilter_directional_fastest(output, coef,
00512 point2d(i, - b),
00513 point2d(i, ncols - 1 + b),
00514 ncols + 2 * b,
00515 dpoint2d(0, 1),
00516 bdr);
00517 }
00518
00519 trace::exiting("linear::gaussian_directional_2d");
00520 return output;
00521 }
00522
00523 # endif // ! MLN_INCLUDE_ONLY
00524
00525 }
00526
00527 }
00528
00529
00530 #endif // ! MLN_LINEAR_GAUSSIAN_DIRECTIONAL_2D_HH