Milena (Olena)
User documentation 2.0a Id
|
00001 // Copyright (C) 2001, 2002, 2003, 2004, 2007, 2008, 2009, 2010, 2011 00002 // EPITA Research and Development Laboratory (LRDE) 00003 // 00004 // This file is part of Olena. 00005 // 00006 // Olena is free software: you can redistribute it and/or modify it under 00007 // the terms of the GNU General Public License as published by the Free 00008 // Software Foundation, version 2 of the License. 00009 // 00010 // Olena is distributed in the hope that it will be useful, 00011 // but WITHOUT ANY WARRANTY; without even the implied warranty of 00012 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 00013 // General Public License for more details. 00014 // 00015 // You should have received a copy of the GNU General Public License 00016 // along with Olena. If not, see <http://www.gnu.org/licenses/>. 00017 // 00018 // As a special exception, you may use this file as part of a free 00019 // software project without restriction. Specifically, if other files 00020 // instantiate templates or use macros or inline functions from this 00021 // file, or you compile this file and link it with other files to produce 00022 // an executable, this file does not by itself cause the resulting 00023 // executable to be covered by the GNU General Public License. This 00024 // exception does not however invalidate any other reasons why the 00025 // executable file might be covered by the GNU General Public License. 00026 00027 #ifndef MLN_LINEAR_GAUSSIAN_HH 00028 # define MLN_LINEAR_GAUSSIAN_HH 00029 00036 00037 # include <vector> 00038 # include <cmath> 00039 00040 # include <mln/core/concept/image.hh> 00041 # include <mln/core/alias/point2d.hh> 00042 # include <mln/core/alias/dpoint1d.hh> 00043 # include <mln/core/alias/dpoint3d.hh> 00044 # include <mln/extension/adjust_fill.hh> 00045 # include <mln/geom/ncols.hh> 00046 # include <mln/geom/nrows.hh> 00047 # include <mln/geom/min_col.hh> 00048 # include <mln/geom/max_col.hh> 00049 # include <mln/geom/min_row.hh> 00050 # include <mln/geom/max_row.hh> 00051 # include <mln/geom/min_sli.hh> 00052 # include <mln/geom/max_sli.hh> 00053 # include <mln/geom/ninds.hh> 00054 # include <mln/geom/nslis.hh> 00055 # include <mln/data/paste.hh> 00056 # include <mln/data/stretch.hh> 00057 # include <mln/algebra/vec.hh> 00058 00059 00060 namespace mln 00061 { 00062 00063 namespace linear 00064 { 00065 00070 template <typename I> 00071 mln_concrete(I) 00072 gaussian(const Image<I>& input, float sigma); 00073 00074 00075 template <typename I> 00076 mln_concrete(I) 00077 gaussian(const Image<I>& input, float sigma, int dir); 00078 00079 00080 # ifndef MLN_INCLUDE_ONLY 00081 00082 namespace impl 00083 { 00084 00085 typedef float norm_fun(float, float, 00086 float, float, 00087 float, float, 00088 float, float, 00089 float, float, 00090 int&); 00091 00092 struct recursivefilter_coef_ 00093 { 00094 00098 recursivefilter_coef_(float a0, float a1, 00099 float b0, float b1, 00100 float c0, float c1, 00101 float w0, float w1, 00102 float s, norm_fun norm); 00103 std::vector<float> n, d, nm, dm; 00104 }; 00105 00106 inline 00107 recursivefilter_coef_::recursivefilter_coef_(float a0, float a1, 00108 float b0, float b1, 00109 float c0, float c1, 00110 float w0, float w1, 00111 float s, norm_fun norm) 00112 { 00113 n.reserve(5); 00114 d.reserve(5); 00115 nm.reserve(5); 00116 dm.reserve(5); 00117 00118 b0 /= s; 00119 b1 /= s; 00120 w0 /= s; 00121 w1 /= s; 00122 00123 float sin0 = std::sin(w0); 00124 float sin1 = std::sin(w1); 00125 float cos0 = std::cos(w0); 00126 float cos1 = std::cos(w1); 00127 00128 int sign = 1; 00129 float n_ = norm(a0, a1, b0, b1, c0, c1, cos0, sin0, cos1, sin1, sign); 00130 00131 a0 /= n_; 00132 a1 /= n_; 00133 c0 /= n_; 00134 c1 /= n_; 00135 00136 n[3] = 00137 std::exp(-b1 - 2*b0) * (c1 * sin1 - cos1 * c0) + 00138 std::exp(-b0 - 2*b1) * (a1 * sin0 - cos0 * a0); 00139 n[2] = 00140 2 * std::exp(-b0 - b1) * ((a0 + c0) * cos1 * cos0 - 00141 cos1 * a1 * sin0 - 00142 cos0 * c1 * sin1) + 00143 c0 * std::exp(-2*b0) + a0 * std::exp(-2*b1); 00144 n[1] = 00145 std::exp(-b1) * (c1 * sin1 - (c0 + 2 * a0) * cos1) + 00146 std::exp(-b0) * (a1 * sin0 - (2 * c0 + a0) * cos0); 00147 n[0] = 00148 a0 + c0; 00149 00150 d[4] = 00151 std::exp(-2 * b0 - 2 * b1); 00152 d[3] = 00153 -2 * cos0 * std::exp(-b0 - 2*b1) - 00154 2 * cos1 * std::exp(-b1 - 2*b0); 00155 d[2] = 00156 4 * cos1 * cos0 * std::exp(-b0 - b1) + 00157 std::exp(-2*b1) + std::exp(-2*b0); 00158 d[1] = 00159 -2 * std::exp(-b1) * cos1 - 2 * std::exp(-b0) * cos0; 00160 00161 for (unsigned i = 1; i <= 3; ++i) 00162 { 00163 dm[i] = d[i]; 00164 nm[i] = float(sign) * (n[i] - d[i] * n[0]); 00165 } 00166 dm[4] = d[4]; 00167 nm[4] = float(sign) * (-d[4] * n[0]); 00168 } 00169 00170 00171 00172 template <class WorkType, class I> 00173 inline 00174 void 00175 recursivefilter_(I& ima, 00176 const recursivefilter_coef_& c, 00177 const mln_psite(I)& start, 00178 const mln_psite(I)& finish, 00179 int len, 00180 const mln_deduce(I, psite, delta)& d) 00181 { 00182 std::vector<WorkType> tmp1(len); 00183 std::vector<WorkType> tmp2(len); 00184 00185 // The fourth degree approximation implies to have a special 00186 // look on the four first points we consider that there is 00187 // no signal before 0 (to be discussed) 00188 00189 // -- 00190 // Causal part 00191 00192 tmp1[0] = 00193 c.n[0] * ima(start); 00194 00195 tmp1[1] = 00196 c.n[0] * ima(start + d) 00197 + c.n[1] * ima(start) 00198 - c.d[1] * tmp1[0]; 00199 00200 tmp1[2] = 00201 c.n[0] * ima(start + d + d) 00202 + c.n[1] * ima(start + d) 00203 + c.n[2] * ima(start) 00204 - c.d[1] * tmp1[1] 00205 - c.d[2] * tmp1[0]; 00206 00207 tmp1[3] = 00208 c.n[0] * ima(start + d + d + d) 00209 + c.n[1] * ima(start + d + d) 00210 + c.n[2] * ima(start + d) 00211 + c.n[3] * ima(start) 00212 - c.d[1] * tmp1[2] - c.d[2] * tmp1[1] 00213 - c.d[3] * tmp1[0]; 00214 00215 mln_psite(I) current(start + d + d + d + d); 00216 for (mln_deduce(I, site, coord) i = 4; i < len; ++i) 00217 { 00218 tmp1[i] = 00219 c.n[0] * ima(current) 00220 + c.n[1] * ima(current - d) 00221 + c.n[2] * ima(current - d - d) 00222 + c.n[3] * ima(current - d - d - d) 00223 - c.d[1] * tmp1[i - 1] - c.d[2] * tmp1[i - 2] 00224 - c.d[3] * tmp1[i - 3] - c.d[4] * tmp1[i - 4]; 00225 current = current + d; 00226 } 00227 00228 // Non causal part 00229 00230 tmp2[len - 1] = WorkType(); // FIXME : = 0, literal::zero ...? 00231 00232 tmp2[len - 2] = 00233 c.nm[1] * ima(finish); 00234 00235 tmp2[len - 3] = 00236 c.nm[1] * ima(finish - d) 00237 + c.nm[2] * ima(finish) 00238 - c.dm[1] * tmp2[len - 2]; 00239 00240 tmp2[len - 4] = 00241 c.nm[1] * ima(finish - d - d) 00242 + c.nm[2] * ima(finish - d) 00243 + c.nm[3] * ima(finish) 00244 - c.dm[1] * tmp2[len - 3] 00245 - c.dm[2] * tmp2[len - 2]; 00246 00247 current = finish - d - d - d ; 00248 00249 for (int i = len - 5; i >= 0; --i) 00250 { 00251 tmp2[i] = 00252 c.nm[1] * ima(current) 00253 + c.nm[2] * ima(current + d) 00254 + c.nm[3] * ima(current + d + d) 00255 + c.nm[4] * ima(current + d + d + d) 00256 - c.dm[1] * tmp2[i + 1] - c.dm[2] * tmp2[i + 2] 00257 - c.dm[3] * tmp2[i + 3] - c.dm[4] * tmp2[i + 4]; 00258 current = current - d; 00259 } 00260 00261 // Combine results from causal and non-causal parts. 00262 current = start; 00263 for (int i = 0; i < len; ++i) 00264 { 00265 ima(current) = tmp1[i] + tmp2[i]; 00266 current = current + d; 00267 } 00268 } 00269 00270 00271 inline 00272 float gaussian_norm_coef_(float a0, float a1, 00273 float b0, float b1, 00274 float c0, float c1, 00275 float cos0, float sin0, 00276 float cos1, float sin1, 00277 int& sign) 00278 { 00279 float expb0 = std::exp(b0); 00280 float exp2b0 = std::exp(2.f * b0); 00281 00282 float scale0 = 1 + exp2b0 - 2 * cos0 * expb0; 00283 float scaleA = 2 * a1 * sin0 * expb0 - a0 * (1 - exp2b0); 00284 00285 float expb1 = std::exp(b1); 00286 float exp2b1 = std::exp(2.f * b1); 00287 00288 float scale1 = 1 + exp2b1 - 2 * cos1 * expb1; 00289 float scaleC = 2 * c1 * sin1 * expb1 - c0 * (1 - exp2b1); 00290 00291 float sumA = scaleA / scale0; 00292 float sumC = scaleC / scale1; 00293 00294 sign = 1; 00295 00296 return (sumA + sumC); 00297 } 00298 00299 inline 00300 float gaussian_1st_deriv_coef_norm_(float a0, float a1, 00301 float b0, float b1, 00302 float c0, float c1, 00303 float cos0, float sin0, 00304 float cos1, float sin1, 00305 int& sign) 00306 { 00307 float expb0 = std::exp(b0); 00308 float exp2b0 = std::exp(2.f * b0); 00309 00310 float scale0 = 1 + exp2b0 - 2 * cos0 * expb0; 00311 scale0 *= scale0; 00312 float scaleA = - 2 * a1 * sin0 * expb0 * (1 - exp2b0) + 00313 2 * a0 * expb0 * (2 * expb0 - cos0 * (1 + exp2b0)); 00314 00315 float expb1 = std::exp(b1); 00316 float exp2b1 = std::exp(2.f * b1); 00317 00318 float scale1 = 1 + exp2b1 - 2 * cos1 * expb1; 00319 scale1 *= scale1; 00320 float scaleC = - 2 * c1 * sin1 * expb1 * (1 - exp2b1) + 00321 2 * c0 * expb1 * (2 * expb1 - cos1 * (1 + exp2b1)); 00322 00323 float sumA = scaleA / scale0; 00324 float sumC = scaleC / scale1; 00325 00326 sign = -1; 00327 00328 return (sumA + sumC); 00329 } 00330 00331 inline 00332 float gaussian_2nd_deriv_coef_norm_(float a0, float a1, 00333 float b0, float b1, 00334 float c0, float c1, 00335 float cos0, float sin0, 00336 float cos1, float sin1, 00337 int& sign) 00338 { 00339 float expb0 = std::exp(b0); 00340 float exp2b0 = std::exp(2.f * b0); 00341 00342 float scale0 = 1 + exp2b0 - 2 * cos0 * expb0; 00343 scale0 *= scale0 * scale0; 00344 00345 float scaleA = a1 * sin0 * expb0 * 00346 (1 + expb0 * (2 * cos0 * (1 + exp2b0) + exp2b0 - 6)) + 00347 a0 * expb0 * (2 * expb0 * (2 - cos0 * cos0) * 00348 (1 - exp2b0) - cos0 * (1 - exp2b0 * exp2b0)); 00349 00350 float expb1 = std::exp(b1); 00351 float exp2b1 = std::exp(2.f * b1); 00352 00353 float scale1 = 1 + exp2b1 - 2 * cos1 * expb1; 00354 scale1 *= scale1 * scale1; 00355 00356 float scaleC = c1 * sin1 * expb1 * 00357 (1 + expb1 * (2 * cos1 * (1 + exp2b1) + exp2b1 - 6)) + 00358 c0 * expb1 * (2 * expb1 * (2 - cos1 * cos1) * 00359 (1 - exp2b1) - cos1 * (1 - exp2b1 * exp2b1)); 00360 00361 float sumA = scaleA / scale0; 00362 float sumC = scaleC / scale1; 00363 sign = 1; 00364 00365 return (sumA + sumC); 00366 } 00367 00368 00369 template <class I, class F> 00370 inline 00371 void 00372 generic_filter_(trait::image::dimension::one_d, 00373 Image<I>& img_, const F& coef, int dir) 00374 { 00375 I& img = exact(img_); 00376 mln_precondition(dir < I::site::dim); 00377 00378 00379 recursivefilter_<mln_value(I)>(img, coef, 00380 point1d(static_cast<def::coord>(-img.border())), 00381 point1d(static_cast<def::coord>(geom::ninds(img) - 1 + 00382 img.border())), 00383 geom::ninds(img) + 2 * img.border(), 00384 dpoint1d(1)); 00385 } 00386 00387 template <class I, class F> 00388 inline 00389 void 00390 generic_filter_(trait::image::dimension::two_d, 00391 Image<I>& img_, const F& coef, int dir) 00392 { 00393 I& img = exact(img_); 00394 00395 mln_precondition(dir < I::site::dim); 00396 00397 00398 if (dir == 0) 00399 { 00400 // Apply on rows. 00401 for (def::coord j = geom::min_col(img); j <= geom::max_col(img); ++j) 00402 recursivefilter_<mln_value(I)>(img, coef, 00403 point2d(static_cast<def::coord>(geom::min_row(img) - img.border()), 00404 static_cast<def::coord>(j)), 00405 point2d(static_cast<def::coord>(geom::max_row(img) + 00406 img.border()), 00407 static_cast<def::coord>(j)), 00408 geom::nrows(img) + 2 * img.border(), 00409 dpoint2d(1, 0)); 00410 } 00411 00412 if (dir == 1) 00413 { 00414 // Apply on columns. 00415 for (def::coord i = geom::min_row(img); i <= geom::max_row(img); ++i) 00416 recursivefilter_<mln_value(I)>(img, coef, 00417 point2d(static_cast<def::coord>(i), 00418 static_cast<def::coord>(geom::min_col(img) - img.border())), 00419 point2d(static_cast<def::coord>(i), 00420 static_cast<def::coord>(geom::max_col(img) + 00421 img.border())), 00422 geom::ncols(img) + 2 * img.border(), 00423 dpoint2d(0, 1)); 00424 } 00425 } 00426 00427 template <class I, class F> 00428 inline 00429 void 00430 generic_filter_(trait::image::dimension::three_d, 00431 Image<I>& img_, const F& coef, int dir) 00432 { 00433 I& img = exact(img_); 00434 mln_precondition(dir < I::site::dim); 00435 00436 if (dir == 0) 00437 { 00438 // Apply on slices. 00439 for (def::coord j = geom::min_row(img); j <= geom::max_row(img); ++j) 00440 for (def::coord k = geom::min_col(img); k <= geom::max_col(img); ++k) 00441 recursivefilter_<mln_value(I)>(img, coef, 00442 point3d(static_cast<def::coord>(-img.border()), 00443 static_cast<def::coord>(j), 00444 static_cast<def::coord>(k)), 00445 point3d(static_cast<def::coord>(geom::nslis(img) - 1 + 00446 img.border()), 00447 static_cast<def::coord>(j), 00448 static_cast<def::coord>(k)), 00449 geom::nslis(img) + 2 * 00450 img.border(), 00451 dpoint3d(1, 0, 0)); 00452 } 00453 00454 00455 if (dir == 1) 00456 { 00457 // Apply on rows. 00458 for (def::coord i = geom::min_sli(img); i <= geom::max_sli(img); ++i) 00459 for (def::coord k = geom::min_col(img); k <= geom::max_col(img); ++k) 00460 recursivefilter_<mln_value(I)>(img, coef, 00461 point3d(static_cast<def::coord>(i), 00462 static_cast<def::coord>(-img.border()), 00463 static_cast<def::coord>(k)), 00464 point3d(static_cast<def::coord>(i), 00465 static_cast<def::coord>(geom::nrows(img) - 1 + 00466 img.border()), 00467 static_cast<def::coord>(k)), 00468 geom::nrows(img) + 2 * 00469 img.border(), 00470 dpoint3d(0, 1, 0)); 00471 } 00472 00473 if (dir == 2) 00474 { 00475 // Apply on columns. 00476 for (def::coord i = geom::min_sli(img); i <= geom::max_sli(img); ++i) 00477 for (def::coord j = geom::min_row(img); j <= geom::max_row(img); ++j) 00478 recursivefilter_<mln_value(I)>(img, coef, 00479 point3d(static_cast<def::coord>(i), 00480 static_cast<def::coord>(j), 00481 static_cast<def::coord>(-img.border())), 00482 point3d(static_cast<def::coord>(i), 00483 static_cast<def::coord>(j), 00484 static_cast<def::coord>(geom::ncols(img) - 00485 1 + img.border())), 00486 geom::ncols(img) + 2 * 00487 img.border(), 00488 dpoint3d(0, 0, 1)); 00489 } 00490 } 00491 00492 00493 00494 template <class I, class F, class O> 00495 inline 00496 void 00497 generic_filter_common_(trait::value::nature::floating, 00498 const Image<I>& in, 00499 const F& coef, 00500 float sigma, 00501 Image<O>& out) 00502 { 00503 mln_ch_value(O, float) work_img(exact(in).domain()); 00504 data::paste(in, work_img); 00505 extension::adjust_fill(work_img, 4, 0); 00506 00507 // On tiny sigma, Derich algorithm doesn't work. 00508 // It is the same thing that to convolve with a Dirac. 00509 if (sigma > 0.006) 00510 for (int i = 0; i < I::site::dim; ++i) 00511 generic_filter_(mln_trait_image_dimension(I)(), 00512 work_img, coef, i); 00513 00514 // We don't need to convert work_img 00515 data::paste(work_img, out); 00516 } 00517 00518 template <class I, class F, class O> 00519 inline 00520 void 00521 generic_filter_common_(trait::value::nature::floating, 00522 const Image<I>& in, 00523 const F& coef, 00524 float sigma, 00525 Image<O>& out, 00526 int dir) 00527 { 00528 mln_ch_value(O, float) work_img(exact(in).domain()); 00529 data::paste(in, work_img); 00530 extension::adjust_fill(work_img, 4, 0); 00531 00532 // On tiny sigma, Derich algorithm doesn't work. 00533 // It is the same thing that to convolve with a Dirac. 00534 if (sigma > 0.006) 00535 generic_filter_(mln_trait_image_dimension(I)(), 00536 work_img, coef, dir); 00537 00538 // We don't need to convert work_img 00539 data::paste(work_img, out); 00540 } 00541 00542 00543 template <class I, class F, class O> 00544 inline 00545 void 00546 generic_filter_common_(trait::value::nature::scalar, 00547 const Image<I>& in, 00548 const F& coef, 00549 float sigma, 00550 Image<O>& out) 00551 { 00552 mln_ch_value(O, float) work_img(exact(in).domain()); 00553 data::paste(in, work_img); 00554 extension::adjust_fill(work_img, 4, 0); 00555 00556 // On tiny sigma, Derich algorithm doesn't work. 00557 // It is the same thing that to convolve with a Dirac. 00558 if (sigma > 0.006) 00559 for (int i = 0; i < I::site::dim; ++i) 00560 generic_filter_(mln_trait_image_dimension(I)(), 00561 work_img, coef, i); 00562 00563 // Convert work_img into result type 00564 data::paste(data::stretch(mln_value(I)(), work_img), out); 00565 } 00566 00567 template <class I, class F, class O> 00568 inline 00569 void 00570 generic_filter_common_(trait::value::nature::scalar, 00571 const Image<I>& in, 00572 const F& coef, 00573 float sigma, 00574 Image<O>& out, 00575 int dir) 00576 { 00577 mln_ch_value(O, float) work_img(exact(in).domain()); 00578 data::paste(in, work_img); 00579 extension::adjust_fill(work_img, 4, 0); 00580 00581 // On tiny sigma, Derich algorithm doesn't work. 00582 // It is the same thing that to convolve with a Dirac. 00583 if (sigma > 0.006) 00584 generic_filter_(mln_trait_image_dimension(I)(), 00585 work_img, coef, dir); 00586 00587 // Convert work_img into result type 00588 data::paste(data::stretch(mln_value(I)(), work_img), out); 00589 } 00590 00591 00592 00593 template <class I, class F, class O> 00594 inline 00595 void 00596 generic_filter_common_(trait::value::nature::vectorial, 00597 const Image<I>& in, 00598 const F& coef, 00599 float sigma, 00600 Image<O>& out) 00601 { 00602 // typedef algebra::vec<3, float> vec3f; 00603 // mln_ch_value(O, vec3f) work_img(exact(in).domain()); 00604 // FIXME : paste does not work (rgb8 -> vec3f). 00605 data::paste(in, out); 00606 00607 // On tiny sigma, Derich algorithm doesn't work. 00608 // It is the same thing that to convolve with a Dirac. 00609 if (sigma > 0.006) 00610 for (int i = 0; i < I::site::dim; ++i) 00611 generic_filter_(mln_trait_image_dimension(I)(), 00612 out, coef, i); 00613 } 00614 00615 template <class I, class F, class O> 00616 inline 00617 void 00618 generic_filter_common_(trait::value::nature::vectorial, 00619 const Image<I>& in, 00620 const F& coef, 00621 float sigma, 00622 Image<O>& out, 00623 int dir) 00624 { 00625 // typedef algebra::vec<3, float> vec3f; 00626 // mln_ch_value(O, vec3f) work_img(exact(in).domain()); 00627 // FIXME : paste does not work (rgb8 -> vec3f). 00628 data::paste(in, out); 00629 00630 // On tiny sigma, Derich algorithm doesn't work. 00631 // It is the same thing that to convolve with a Dirac. 00632 if (sigma > 0.006) 00633 generic_filter_(mln_trait_image_dimension(I)(), 00634 out, coef, dir); 00635 } 00636 00637 } // end of namespace mln::linear::impl 00638 00639 // Facade. 00640 00650 template <typename I> 00651 inline 00652 mln_concrete(I) 00653 gaussian(const Image<I>& input, float sigma, int dir) 00654 { 00655 mln_precondition(exact(input).is_valid()); 00656 mln_precondition(dir < I::site::dim); 00657 00658 mln_concrete(I) output; 00659 initialize(output, input); 00660 impl::recursivefilter_coef_ coef(1.68f, 3.735f, 00661 1.783f, 1.723f, 00662 -0.6803f, -0.2598f, 00663 0.6318f, 1.997f, 00664 sigma, impl::gaussian_norm_coef_); 00665 00666 impl::generic_filter_common_(mln_trait_value_nature(mln_value(I))(), 00667 input, coef, sigma, output, dir); 00668 return output; 00669 } 00670 00671 00682 template <typename I> 00683 inline 00684 mln_concrete(I) 00685 gaussian_1st_derivative(const Image<I>& input, float sigma, int dir) 00686 { 00687 mln_precondition(exact(input).is_valid()); 00688 mln_precondition(dir < I::site::dim); 00689 00690 mln_concrete(I) output; 00691 initialize(output, input); 00692 00693 impl::recursivefilter_coef_ 00694 coef(-0.6472f, -4.531f, 00695 1.527f, 1.516f, 00696 0.6494f, 0.9557f, 00697 0.6719f, 2.072f, 00698 sigma, impl::gaussian_1st_deriv_coef_norm_); 00699 00700 impl::generic_filter_common_(mln_trait_value_nature(mln_value(I))(), 00701 input, coef, sigma, output, dir); 00702 return output; 00703 } 00704 00715 template <typename I> 00716 inline 00717 mln_concrete(I) 00718 gaussian_2nd_derivative(const Image<I>& input, float sigma, int dir) 00719 { 00720 mln_precondition(exact(input).is_valid()); 00721 mln_precondition(dir < I::site::dim); 00722 00723 mln_concrete(I) output; 00724 initialize(output, input); 00725 00726 impl::recursivefilter_coef_ 00727 coef(-1.331f, 3.661f, 00728 1.24f, 1.314f, 00729 0.3225f, -1.738f, 00730 0.748f, 2.166f, 00731 sigma, impl::gaussian_2nd_deriv_coef_norm_); 00732 00733 impl::generic_filter_common_(mln_trait_value_nature(mln_value(I))(), 00734 input, coef, sigma, output, dir); 00735 return output; 00736 } 00737 00738 00739 00740 00741 00747 template <typename I> 00748 inline 00749 mln_concrete(I) 00750 gaussian(const Image<I>& input, float sigma) 00751 { 00752 mln_precondition(exact(input).is_valid()); 00753 00754 mln_concrete(I) output; 00755 initialize(output, input); 00756 00757 impl::recursivefilter_coef_ 00758 coef(1.68f, 3.735f, 00759 1.783f, 1.723f, 00760 -0.6803f, -0.2598f, 00761 0.6318f, 1.997f, 00762 sigma, impl::gaussian_norm_coef_); 00763 impl::generic_filter_common_(mln_trait_value_nature(mln_value(I))(), 00764 input, coef, sigma, output); 00765 00766 return output; 00767 } 00768 00769 00776 template <typename I> 00777 inline 00778 mln_concrete(I) 00779 gaussian_1st_derivative(const Image<I>& input, float sigma) 00780 { 00781 mln_precondition(exact(input).is_valid()); 00782 00783 mln_concrete(I) output; 00784 initialize(output, input); 00785 00786 impl::recursivefilter_coef_ 00787 coef(-0.6472f, -4.531f, 00788 1.527f, 1.516f, 00789 0.6494f, 0.9557f, 00790 0.6719f, 2.072f, 00791 sigma, impl::gaussian_1st_deriv_coef_norm_); 00792 impl::generic_filter_common_(mln_trait_value_nature(mln_value(I))(), 00793 input, coef, sigma, output); 00794 return output; 00795 } 00796 00797 00804 template <typename I> 00805 inline 00806 mln_concrete(I) 00807 gaussian_2nd_derivative(const Image<I>& input, float sigma) 00808 { 00809 mln_precondition(exact(input).is_valid()); 00810 00811 mln_concrete(I) output; 00812 initialize(output, input); 00813 00814 impl::recursivefilter_coef_ 00815 coef(-1.331f, 3.661f, 00816 1.24f, 1.314f, 00817 0.3225f, -1.738f, 00818 0.748f, 2.166f, 00819 sigma, impl::gaussian_2nd_deriv_coef_norm_); 00820 impl::generic_filter_common_(mln_trait_value_nature(mln_value(I))(), 00821 input, coef, sigma, output); 00822 return output; 00823 } 00824 00825 # endif // ! MLN_INCLUDE_ONLY 00826 00827 } // end of namespace mln::linear 00828 00829 } // end of namespace mln 00830 00831 00832 #endif // ! MLN_LINEAR_GAUSSIAN_HH