Milena (Olena): gaussian.hh Source File

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