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_IMPL_HH
00027 # define MLN_LINEAR_GAUSSIAN_IMPL_HH
00028
00035
00036 # include <vector>
00037 # include <cmath>
00038
00039 # include <mln/core/concept/image.hh>
00040 # include <mln/core/alias/point2d.hh>
00041 # include <mln/core/alias/dpoint1d.hh>
00042 # include <mln/core/alias/dpoint3d.hh>
00043 # include <mln/extension/adjust_fill.hh>
00044 # include <mln/geom/ncols.hh>
00045 # include <mln/geom/nrows.hh>
00046 # include <mln/geom/ninds.hh>
00047 # include <mln/geom/nslis.hh>
00048 # include <mln/data/paste.hh>
00049 # include <mln/data/stretch.hh>
00050 # include <mln/algebra/vec.hh>
00051
00052 # include <mln/linear/gaussian/internal/coefficients.hh>
00053
00054
00055 namespace mln
00056 {
00057
00058 namespace linear
00059 {
00060
00061 namespace gaussian
00062 {
00063
00064 namespace impl
00065 {
00066
00067
00068 # ifndef MLN_INCLUDE_ONLY
00069
00070 template <class WorkType, class I>
00071 inline
00072 void
00073 recursivefilter_(I& ima,
00074 const internal::coefficients& c,
00075 const mln_psite(I)& start,
00076 const mln_psite(I)& finish,
00077 int len,
00078 const mln_deduce(I, psite, delta)& d)
00079 {
00080 std::vector<WorkType> tmp1(len);
00081 std::vector<WorkType> tmp2(len);
00082
00083
00084
00085
00086
00087
00088
00089
00090 tmp1[0] =
00091 c.n[0] * ima(start);
00092
00093 tmp1[1] =
00094 c.n[0] * ima(start + d)
00095 + c.n[1] * ima(start)
00096 - c.d[1] * tmp1[0];
00097
00098 tmp1[2] =
00099 c.n[0] * ima(start + d + d)
00100 + c.n[1] * ima(start + d)
00101 + c.n[2] * ima(start)
00102 - c.d[1] * tmp1[1]
00103 - c.d[2] * tmp1[0];
00104
00105 tmp1[3] =
00106 c.n[0] * ima(start + d + d + d)
00107 + c.n[1] * ima(start + d + d)
00108 + c.n[2] * ima(start + d)
00109 + c.n[3] * ima(start)
00110 - c.d[1] * tmp1[2] - c.d[2] * tmp1[1]
00111 - c.d[3] * tmp1[0];
00112
00113 mln_psite(I) current(start + d + d + d + d);
00114 for (mln_deduce(I, site, coord) i = 4; i < len; ++i)
00115 {
00116 tmp1[i] =
00117 c.n[0] * ima(current)
00118 + c.n[1] * ima(current - d)
00119 + c.n[2] * ima(current - d - d)
00120 + c.n[3] * ima(current - d - d - d)
00121 - c.d[1] * tmp1[i - 1] - c.d[2] * tmp1[i - 2]
00122 - c.d[3] * tmp1[i - 3] - c.d[4] * tmp1[i - 4];
00123 current = current + d;
00124 }
00125
00126
00127
00128 tmp2[len - 1] = WorkType();
00129
00130 tmp2[len - 2] =
00131 c.nm[1] * ima(finish);
00132
00133 tmp2[len - 3] =
00134 c.nm[1] * ima(finish - d)
00135 + c.nm[2] * ima(finish)
00136 - c.dm[1] * tmp2[len - 2];
00137
00138 tmp2[len - 4] =
00139 c.nm[1] * ima(finish - d - d)
00140 + c.nm[2] * ima(finish - d)
00141 + c.nm[3] * ima(finish)
00142 - c.dm[1] * tmp2[len - 3]
00143 - c.dm[2] * tmp2[len - 2];
00144
00145 current = finish - d - d - d ;
00146
00147 for (int i = len - 5; i >= 0; --i)
00148 {
00149 tmp2[i] =
00150 c.nm[1] * ima(current)
00151 + c.nm[2] * ima(current + d)
00152 + c.nm[3] * ima(current + d + d)
00153 + c.nm[4] * ima(current + d + d + d)
00154 - c.dm[1] * tmp2[i + 1] - c.dm[2] * tmp2[i + 2]
00155 - c.dm[3] * tmp2[i + 3] - c.dm[4] * tmp2[i + 4];
00156 current = current - d;
00157 }
00158
00159
00160 current = start;
00161 for (int i = 0; i < len; ++i)
00162 {
00163 ima(current) = tmp1[i] + tmp2[i];
00164 current = current + d;
00165 }
00166 }
00167
00168
00169 template <class I, class F>
00170 inline
00171 void
00172 generic_filter_(trait::image::dimension::one_d,
00173 Image<I>& img_, const F& coef, int dir)
00174 {
00175 I& img = exact(img_);
00176 typedef mln_site(I) S;
00177
00178 mln_precondition(dir < S::dim);
00179
00180 recursivefilter_<mln_value(I)>(img, coef,
00181 point1d(static_cast<def::coord>(-img.border())),
00182 point1d(static_cast<def::coord>(geom::ninds(img) - 1 +
00183 img.border())),
00184 geom::ninds(img) + 2 * img.border(),
00185 dpoint1d(1));
00186 }
00187
00188 template <class I, class F>
00189 inline
00190 void
00191 generic_filter_(trait::image::dimension::two_d,
00192 Image<I>& img_, const F& coef, int dir)
00193 {
00194 I& img = exact(img_);
00195 typedef mln_site(I) S;
00196
00197 mln_precondition(dir < S::dim);
00198
00199 if (dir == 0)
00200 {
00201
00202 for (unsigned j = 0; j < geom::ncols(img); ++j)
00203 recursivefilter_<mln_value(I)>(img, coef,
00204 point2d(static_cast<def::coord>(-img.border()),
00205 static_cast<def::coord>(j)),
00206 point2d(static_cast<def::coord>(geom::nrows(img) - 1 +
00207 img.border()),
00208 static_cast<def::coord>(j)),
00209 geom::nrows(img) + 2 * img.border(),
00210 dpoint2d(1, 0));
00211 }
00212
00213 if (dir == 1)
00214 {
00215
00216 for (unsigned i = 0; i < geom::nrows(img); ++i)
00217 recursivefilter_<mln_value(I)>(img, coef,
00218 point2d(static_cast<def::coord>(i),
00219 static_cast<def::coord>(-img.border())),
00220 point2d(static_cast<def::coord>(i),
00221 static_cast<def::coord>(geom::ncols(img) - 1 +
00222 img.border())),
00223 geom::ncols(img) + 2 * img.border(),
00224 dpoint2d(0, 1));
00225 }
00226 }
00227
00228 template <class I, class F>
00229 inline
00230 void
00231 generic_filter_(trait::image::dimension::three_d,
00232 Image<I>& img_, const F& coef, int dir)
00233 {
00234 I& img = exact(img_);
00235 typedef mln_site(I) S;
00236
00237 mln_precondition(dir < S::dim);
00238
00239 if (dir == 0)
00240 {
00241
00242 for (unsigned j = 0; j < geom::nrows(img); ++j)
00243 for (unsigned k = 0; k < geom::ncols(img); ++k)
00244 recursivefilter_<mln_value(I)>(img, coef,
00245 point3d(static_cast<def::coord>(-img.border()),
00246 static_cast<def::coord>(j),
00247 static_cast<def::coord>(k)),
00248 point3d(static_cast<def::coord>(geom::nslis(img) - 1 +
00249 img.border()),
00250 static_cast<def::coord>(j),
00251 static_cast<def::coord>(k)),
00252 geom::nslis(img) + 2 *
00253 img.border(),
00254 dpoint3d(1, 0, 0));
00255 }
00256
00257
00258 if (dir == 1)
00259 {
00260
00261 for (unsigned i = 0; i < geom::nslis(img); ++i)
00262 for (unsigned k = 0; k < geom::ncols(img); ++k)
00263 recursivefilter_<mln_value(I)>(img, coef,
00264 point3d(static_cast<def::coord>(i),
00265 static_cast<def::coord>(-img.border()),
00266 static_cast<def::coord>(k)),
00267 point3d(static_cast<def::coord>(i),
00268 static_cast<def::coord>(geom::nrows(img) - 1 +
00269 img.border()),
00270 static_cast<def::coord>(k)),
00271 geom::nrows(img) + 2 *
00272 img.border(),
00273 dpoint3d(0, 1, 0));
00274 }
00275
00276 if (dir == 2)
00277 {
00278
00279 for (unsigned i = 0; i < geom::nslis(img); ++i)
00280 for (unsigned j = 0; j < geom::nrows(img); ++i)
00281 recursivefilter_<mln_value(I)>(img, coef,
00282 point3d(static_cast<def::coord>(i),
00283 static_cast<def::coord>(j),
00284 static_cast<def::coord>(-img.border())),
00285 point3d(static_cast<def::coord>(i),
00286 static_cast<def::coord>(j),
00287 static_cast<def::coord>(geom::ncols(img) -
00288 1 + img.border())),
00289 geom::ncols(img) + 2 *
00290 img.border(),
00291 dpoint3d(0, 0, 1));
00292 }
00293 }
00294
00295
00296
00297 template <class I, class F, class O>
00298 inline
00299 void
00300 generic_filter_common_(trait::value::nature::floating,
00301 const Image<I>& in,
00302 const F& coef,
00303 double sigma,
00304 Image<O>& out)
00305 {
00306 typedef mln_site(I) S;
00307
00308 mln_ch_value(O, double) work_img(exact(in).domain());
00309 data::paste(in, work_img);
00310 extension::adjust_fill(work_img, 4, 0);
00311
00312
00313
00314 if (sigma > 0.006)
00315 for (int i = 0; i < S::dim; ++i)
00316 generic_filter_(mln_trait_image_dimension(I)(),
00317 work_img, coef, i);
00318
00319
00320 data::paste(work_img, out);
00321 }
00322
00323 template <class I, class F, class O>
00324 inline
00325 void
00326 generic_filter_common_(trait::value::nature::floating,
00327 const Image<I>& in,
00328 const F& coef,
00329 double sigma,
00330 Image<O>& out,
00331 int dir)
00332 {
00333 mln_ch_value(O, double) work_img(exact(in).domain());
00334 data::paste(in, work_img);
00335 extension::adjust_fill(work_img, 4, 0);
00336
00337
00338
00339 if (sigma > 0.006)
00340 generic_filter_(mln_trait_image_dimension(I)(),
00341 work_img, coef, dir);
00342
00343
00344 data::paste(work_img, out);
00345 }
00346
00347
00348 template <class I, class F, class O>
00349 inline
00350 void
00351 generic_filter_common_(trait::value::nature::scalar,
00352 const Image<I>& in,
00353 const F& coef,
00354 double sigma,
00355 Image<O>& out)
00356 {
00357 typedef mln_site(I) S;
00358
00359 mln_ch_value(O, double) work_img(exact(in).domain());
00360 data::paste(in, work_img);
00361 extension::adjust_fill(work_img, 4, 0);
00362
00363
00364
00365 if (sigma > 0.006)
00366 for (int i = 0; i < S::dim; ++i)
00367 generic_filter_(mln_trait_image_dimension(I)(),
00368 work_img, coef, i);
00369
00370
00371 data::paste(data::stretch(mln_value(I)(), work_img), out);
00372 }
00373
00374 template <class I, class F, class O>
00375 inline
00376 void
00377 generic_filter_common_(trait::value::nature::scalar,
00378 const Image<I>& in,
00379 const F& coef,
00380 double sigma,
00381 Image<O>& out,
00382 int dir)
00383 {
00384 mln_ch_value(O, double) work_img(exact(in).domain());
00385 data::paste(in, work_img);
00386 extension::adjust_fill(work_img, 4, 0);
00387
00388
00389
00390 if (sigma > 0.006)
00391 generic_filter_(mln_trait_image_dimension(I)(),
00392 work_img, coef, dir);
00393
00394
00395 data::paste(data::stretch(mln_value(I)(), work_img), out);
00396 }
00397
00398
00399
00400 template <class I, class F, class O>
00401 inline
00402 void
00403 generic_filter_common_(trait::value::nature::vectorial,
00404 const Image<I>& in,
00405 const F& coef,
00406 double sigma,
00407 Image<O>& out)
00408 {
00409 typedef mln_site(I) S;
00410
00411
00412
00413
00414 data::paste(in, out);
00415
00416
00417
00418 if (sigma > 0.006)
00419 for (int i = 0; i < S::dim; ++i)
00420 generic_filter_(mln_trait_image_dimension(I)(),
00421 out, coef, i);
00422 }
00423
00424 template <class I, class F, class O>
00425 inline
00426 void
00427 generic_filter_common_(trait::value::nature::vectorial,
00428 const Image<I>& in,
00429 const F& coef,
00430 double sigma,
00431 Image<O>& out,
00432 int dir)
00433 {
00434
00435
00436
00437 data::paste(in, out);
00438
00439
00440
00441 if (sigma > 0.006)
00442 generic_filter_(mln_trait_image_dimension(I)(),
00443 out, coef, dir);
00444 }
00445
00446
00447 # endif // ! MLN_INCLUDE_ONLY
00448
00449
00450 }
00451
00452 }
00453
00454 }
00455
00456 }
00457
00458
00459 #endif // ! MLN_LINEAR_GAUSSIAN_IMPL_HH