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