Milena (Olena): gaussian.hh Source File

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