Milena (Olena): gaussian_directional

00001 // Copyright (C) 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_DIRECTIONAL_2D_HH
00027 # define MLN_LINEAR_GAUSSIAN_DIRECTIONAL_2D_HH
00028 
00034 
00035 #include <mln/core/image/image2d.hh>
00036 #include <mln/extension/adjust_fill.hh>
00037 #include <mln/extension/adjust_duplicate.hh>
00038 
00039 
00040 
00041 namespace mln
00042 {
00043 
00044   namespace linear
00045   {
00046 
00047 
00048     template <typename I>
00049     mln_concrete(I)
00050     gaussian_directional_2d(const Image<I>& input,
00051                             unsigned dir, double sigma,
00052                             const mln_value(I)& bdr);
00053 
00054 
00055 
00056 # ifndef MLN_INCLUDE_ONLY
00057 
00058     namespace my
00059     {
00060 
00061       struct recursivefilter_coef_
00062       {
00063         enum FilterType { DericheGaussian,
00064                           DericheGaussianFirstDerivative,
00065                           DericheGaussianSecondDerivative };
00066 
00067         std::vector<double> n, d, nm, dm;
00068         double sumA, sumC;
00069 
00070         recursivefilter_coef_(double a0, double a1,
00071                               double b0, double b1,
00072                               double c0, double c1,
00073                               double w0, double w1,
00074                               double s, FilterType filter_type)
00075         {
00076           n.reserve(5);
00077           d.reserve(5);
00078           nm.reserve(5);
00079           dm.reserve(5);
00080 
00081           b0 /= s;
00082           b1 /= s;
00083           w0 /= s;
00084           w1 /= s;
00085 
00086           double sin0   = sin(w0);
00087           double sin1   = sin(w1);
00088           double cos0   = cos(w0);
00089           double cos1   = cos(w1);
00090 
00091           switch (filter_type) {
00092 
00093           case DericheGaussian :
00094             {
00095               sumA  =
00096                 (2.0 * a1 * exp( b0 ) * cos0 * cos0 - a0 * sin0 * exp( 2.0 * b0 )
00097                  + a0 * sin0 - 2.0 * a1 * exp( b0 )) /
00098                 (( 2.0 * cos0 * exp( b0 ) - exp( 2.0 * b0 ) - 1 ) * sin0);
00099 
00100               sumC  =
00101                 (2.0 * c1 * exp( b1 ) * cos1 * cos1 - c0 * sin1 * exp( 2.0 * b1 )
00102                  + c0 * sin1 - 2.0 * c1 * exp( b1 ))
00103                 / (( 2.0 * cos1 * exp( b1 ) - exp( 2.0 * b1 ) - 1 ) * sin1);
00104               break;
00105             }
00106 
00107           case DericheGaussianFirstDerivative :
00108             {
00109               sumA  = -2.f *
00110                 (a0 * cos0 - a1 * sin0 + a1 * sin0 * exp( 2.0 * b0 )
00111                  + a0 * cos0 * exp( 2.0 * b0 ) - 2.0 * a0 * exp( b0 ))
00112                 * exp( b0 )
00113                 / (exp( 4.0 * b0 ) - 4.0 * cos0 * exp( 3.0 * b0 )
00114                    + 2.0 * exp( 2.0 * b0 ) + 4.0 * cos0 * cos0 * exp( 2.0 * b0 )
00115                    + 1.0 - 4.0 * cos0 * exp( b0 ));
00116               sumC  = - 2.f *
00117                 (c0 * cos1 - c1 * sin1 + c1 * sin1 * exp( 2.0 * b1 )
00118                  + c0 * cos1 * exp( 2.0 * b1 ) - 2.0 * c0 * exp( b1 ))
00119                 * exp( b1 ) /
00120                 (exp( 4.0 * b1 ) - 4.0 * cos1 * exp( 3.0 * b1 )
00121                  + 2.0 * exp( 2.0 * b1 ) + 4.0 * cos1 * cos1 * exp( 2.0 * b1 )
00122                  + 1.0 - 4.0 * cos1 * exp( b1 ));
00123               break;
00124             }
00125 
00126           case DericheGaussianSecondDerivative :
00127             {
00128               double aux;
00129               aux   =
00130                 12.0 * cos0 * exp( 3.0 * b0 ) - 3.0 * exp( 2.0 * b0 )
00131                 + 8.0 * cos0 * cos0 * cos0 * exp( 3.0 * b0 ) - 12.0 * cos0 * cos0 *
00132                 exp( 4.0 * b0 )
00133                 - (3.0 * exp( 4.0 * b0 ))
00134                 + 6.0 * cos0 * exp( 5.0 * b0 ) -  exp( 6.0 * b0 ) + 6.0 * cos0 * exp
00135                 ( b0 )
00136                 - ( 1.0 + 12.0 * cos0 * cos0 * exp( 2.0 * b0 ) );
00137               sumA  =
00138                 4.0 * a0 * sin0 * exp( 3.0 * b0 ) + a1 * cos0 * cos0 * exp( 4.0 * b0 )
00139                 - ( 4.0 * a0 * sin0 * exp( b0 ) + 6.0 * a1 * cos0 * cos0 * exp( 2.0 * b0 ) )
00140                 + 2.0 * a1 * cos0 * cos0 * cos0 * exp( b0 ) - 2.0 * a1 * cos0 * exp(b0)
00141                 + 2.0 * a1 * cos0 * cos0 * cos0 * exp( 3.0 * b0 ) - 2.0 * a1 * cos0
00142                 * exp( 3.0 * b0 )
00143                 + a1 * cos0 * cos0 - a1 * exp( 4.0 * b0 )
00144                 + 2.0 * a0 * sin0 * cos0 * cos0 * exp( b0 ) - 2.0 * a0 * sin0 * cos0
00145                 * cos0 * exp( 3.0 * b0 )
00146                 - ( a0 * sin0 * cos0 * exp( 4.0 * b0 ) + a1 )
00147                 + 6.0 * a1 * exp( 2.0 * b0 ) + a0 * cos0 * sin0
00148                 * 2.0 * exp( b0 ) / ( aux * sin0 );
00149               aux   =
00150                 12.0 * cos1 * exp( 3.0 * b1 ) - 3.0 * exp( 2.0 * b1 )
00151                 + 8.0 * cos1 * cos1 * cos1 * exp( 3.0 * b1 ) - 12.0 * cos1 * cos1 *
00152                 exp( 4.0 * b1 )
00153                 - 3.0 * exp( 4.0 * b1 )
00154                 + 6.0 * cos1 * exp( 5.0 * b1 ) -  exp( 6.0 * b1 ) + 6.0 * cos1 * exp
00155                 ( b1 )
00156                 - ( 1.0 + 12.0 * cos1 * cos1 * exp( 2.0 * b1 ) );
00157               sumC  = 4.0 * c0 * sin1 * exp( 3.0 * b1 ) + c1 * cos1 * cos1 * exp( 4.0 * b1 )
00158                 - ( 4.0 * c0 * sin1 * exp( b1 ) + 6.0 * c1 * cos1 * cos1 * exp( 2.0 * b1 ) )
00159                 + 2.0 * c1 * cos1 * cos1 * cos1 * exp( b1 ) - 2.0 * c1 * cos1 * exp( b1 )
00160                 + 2.0 * c1 * cos1 * cos1 * cos1 * exp( 3.0 * b1 ) - 2.0 * c1 * cos1
00161                 * exp( 3.0 * b1 )
00162                 + c1 * cos1 * cos1 - c1 * exp( 4.0 * b1 )
00163                 + 2.0 * c0 * sin1 * cos1 * cos1 * exp( b1 ) - 2.0 * c0 * sin1 * cos1
00164                 * cos1 * exp( 3.0 * b1 )
00165                 - ( c0 * sin1 * cos1 * exp( 4.0 * b1 ) + c1 )
00166                 + 6.0 * c1 * exp( 2.0 * b1 ) + c0 * cos1 * sin1
00167                 * 2.0 * exp( b1 ) / ( aux * sin1 );
00168               sumA /= 2;
00169               sumC /= 2;
00170               break;
00171             }
00172           }
00173 
00174           a0 /= (sumA + sumC);
00175           a1 /= (sumA + sumC);
00176           c0 /= (sumA + sumC);
00177           c1 /= (sumA + sumC);
00178 
00179           n[3] =
00180             exp( -b1 - 2*b0 ) * (c1 * sin1 - cos1 * c0) +
00181             exp( -b0 - 2*b1 ) * (a1 * sin0 - cos0 * a0);
00182           n[2] =
00183             2 * exp(-b0 - b1) * ((a0 + c0) * cos1 * cos0 -
00184                                  cos1 * a1 * sin0 -
00185                                  cos0 * c1 * sin1) +
00186             c0 * exp(-2 * b0) + a0 * exp(-2 * b1);
00187           n[1] =
00188             exp(-b1) * (c1 * sin1 - (c0 + 2*a0) * cos1) +
00189             exp(-b0) * (a1 * sin0 - (2*c0 + a0) * cos0);
00190           n[0] =
00191             a0 + c0;
00192 
00193           d[4] = exp(-2 * b0 - 2 * b1);
00194           d[3] =
00195             -2 * cos0 * exp(-b0 - 2*b1) -
00196             2 * cos1 * exp(-b1 - 2*b0);
00197           d[2] =
00198             4 * cos1 * cos0 * exp(-b0 - b1) +
00199             exp(-2*b1) + exp(-2*b0);
00200           d[1] =
00201             -2*exp(-b1) * cos1 - 2 * exp(-b0) * cos0;
00202 
00203           switch (filter_type) {
00204           case DericheGaussian :
00205           case DericheGaussianSecondDerivative :
00206             {
00207               for (unsigned i = 1; i <= 3; ++i)
00208                 {
00209                   dm[i] = d[i];
00210                   nm[i] = n[i] - d[i] * n[0];
00211                 }
00212               dm[4] = d[4];
00213               nm[4] = -d[4] * n[0];
00214               break;
00215             }
00216           case DericheGaussianFirstDerivative :
00217             {
00218               for (unsigned i = 1; i <= 3; ++i)
00219                 {
00220                   dm[i] = d[i];
00221                   nm[i] = - (n[i] - d[i] * n[0]);
00222                 }
00223               dm[4] = d[4];
00224               nm[4] = d[4] * n[0];
00225               break;
00226             }
00227           }
00228         }
00229       };
00230 
00231 
00232     } // end of namespace mln::linear::my
00233 
00234 
00235 
00236 
00237     // FIXME: in the "generic" code below there is no test that we
00238     // actually stay in the image domain...
00239 
00240 
00241     template <typename I, typename C>
00242     inline
00243     void
00244     recursivefilter_directional_generic(I& ima,
00245                                         const C& c,
00246                                         const mln_psite(I)& start,
00247                                         const mln_psite(I)& finish,
00248                                         int len,
00249                                         const mln_deduce(I, psite, delta)& d)
00250     {
00251       std::vector<double> tmp1(len);
00252       std::vector<double> tmp2(len);
00253 
00254       tmp1[0] =
00255         c.n[0]*ima(start);
00256 
00257       tmp1[1] =
00258         c.n[0]*ima(start + d)
00259         + c.n[1]*ima(start)
00260         - c.d[1]*tmp1[0];
00261 
00262       tmp1[2] =
00263         c.n[0]*ima(start + d + d)
00264         + c.n[1]*ima(start + d)
00265         + c.n[2]*ima(start)
00266         - c.d[1]*tmp1[1]
00267         - c.d[2]*tmp1[0];
00268 
00269       tmp1[3] =
00270         c.n[0]*ima(start + d + d + d)
00271         + c.n[1]*ima(start + d + d)
00272         + c.n[2]*ima(start + d)
00273         + c.n[3]*ima(start)
00274         - c.d[1]*tmp1[2] - c.d[2]*tmp1[1]
00275         - c.d[3]*tmp1[0];
00276 
00277       mln_site(I) current = start + d + d + d + d;
00278       for (int i = 4; i < len; ++i)
00279         {
00280           tmp1[i] =
00281             c.n[0]*ima(current)
00282             + c.n[1]*ima(current - d)
00283             + c.n[2]*ima(current - d - d)
00284             + c.n[3]*ima(current - d - d - d)
00285             - c.d[1]*tmp1[i - 1] - c.d[2]*tmp1[i - 2]
00286             - c.d[3]*tmp1[i - 3] - c.d[4]*tmp1[i - 4];
00287           current += d;
00288         }
00289 
00290       // Non causal part
00291 
00292       tmp2[len - 1] = 0;
00293 
00294       tmp2[len - 2] =
00295         c.nm[1]*ima(finish);
00296 
00297       tmp2[len - 3] =
00298         c.nm[1]*ima(finish - d)
00299         + c.nm[2]*ima(finish)
00300         - c.dm[1]*tmp2[len-2];
00301 
00302       tmp2[len - 4] =
00303         c.nm[1]*ima(finish - d - d)
00304         + c.nm[2]*ima(finish - d)
00305         + c.nm[3]*ima(finish)
00306         - c.dm[1]*tmp2[len-3]
00307         - c.dm[2]*tmp2[len-2];
00308 
00309       current = finish - d - d - d ;
00310 
00311       for (int i = len - 5; i >= 0; --i)
00312         {
00313           tmp2[i] =
00314             c.nm[1]*ima(current)
00315             + c.nm[2]*ima(current + d)
00316             + c.nm[3]*ima(current + d + d)
00317             + c.nm[4]*ima(current + d + d + d)
00318             - c.dm[1]*tmp2[i+1] - c.dm[2]*tmp2[i+2]
00319             - c.dm[3]*tmp2[i+3] - c.dm[4]*tmp2[i+4];
00320           current -= d;
00321         }
00322 
00323       // Combine results from causal and non-causal parts.
00324 
00325       current = start;
00326       for (int i = 0; i < len; ++i)
00327         {
00328           ima(current) = (tmp1[i] + tmp2[i]);
00329           current += d;
00330         }
00331     }
00332 
00333 
00334 
00335 
00336     template <typename I, typename C>
00337     inline
00338     void
00339     recursivefilter_directional_fastest(I& ima,
00340                                         const C& c,
00341                                         const mln_psite(I)& start,
00342                                         const mln_psite(I)& finish,
00343                                         int len,
00344                                         const mln_deduce(I, psite, delta)& d,
00345                                         const mln_value(I)& bdr)
00346     {
00347       //       extension::adjust_fill(ima, 5 * int(151 + .50001) + 1, bdr);
00348       // extension::fill(ima, bdr);
00349 
00350       std::vector<double> tmp1(len);
00351       std::vector<double> tmp2(len);
00352 
00353       unsigned delta_offset = ima.delta_index(d);
00354       unsigned
00355         o_start = ima.index_of_point(start),
00356         o_start_d   = o_start +     delta_offset,
00357         o_start_dd  = o_start + 2 * delta_offset,
00358         o_start_ddd = o_start + 3 * delta_offset;
00359 
00360       tmp1[0] =
00361         c.n[0] * ima.element(o_start);
00362 
00363       tmp1[1] = 0
00364         + c.n[0] * ima.element(o_start_d)
00365         + c.n[1] * ima.element(o_start)
00366         - c.d[1] * tmp1[0];
00367 
00368       tmp1[2] = 0
00369         + c.n[0] * ima.element(o_start_dd)
00370         + c.n[1] * ima.element(o_start_d)
00371         + c.n[2] * ima.element(o_start)
00372         - c.d[1] * tmp1[1]
00373         - c.d[2] * tmp1[0];
00374 
00375       tmp1[3] = 0
00376         + c.n[0] * ima.element(o_start_ddd)
00377         + c.n[1] * ima.element(o_start_dd)
00378         + c.n[2] * ima.element(o_start_d)
00379         + c.n[3] * ima.element(o_start)
00380         - c.d[1] * tmp1[2] - c.d[2] * tmp1[1]
00381         - c.d[3] * tmp1[0];
00382 
00383       unsigned
00384         o_current = o_start + 4 * delta_offset,
00385         o_current_d   = o_current -     delta_offset,
00386         o_current_dd  = o_current - 2 * delta_offset,
00387         o_current_ddd = o_current - 3 * delta_offset;
00388 
00389       for (int i = 4; i < len; ++i)
00390         {
00391           tmp1[i] = 0
00392             + c.n[0] * ima.element(o_current)
00393             + c.n[1] * ima.element(o_current_d)
00394             + c.n[2] * ima.element(o_current_dd)
00395             + c.n[3] * ima.element(o_current_ddd)
00396             - c.d[1] * tmp1[i - 1] - c.d[2] * tmp1[i - 2]
00397             - c.d[3] * tmp1[i - 3] - c.d[4] * tmp1[i - 4];
00398           o_current     += delta_offset;
00399           o_current_d   += delta_offset;
00400           o_current_dd  += delta_offset;
00401           o_current_ddd += delta_offset;
00402         }
00403 
00404       // Non causal part
00405 
00406       (void) bdr;
00407       // extension::fill(ima, bdr);
00408 
00409       unsigned
00410         o_finish   = ima.index_of_point(finish),
00411         o_finish_d  = o_finish -     delta_offset,
00412         o_finish_dd = o_finish - 2 * delta_offset;
00413 
00414       tmp2[len - 1] = 0;
00415 
00416       tmp2[len - 2] =
00417         c.nm[1] * ima.element(o_finish);
00418 
00419       tmp2[len - 3] = 0
00420         + c.nm[1] * ima.element(o_finish_d)
00421         + c.nm[2] * ima.element(o_finish)
00422         - c.dm[1] * tmp2[len-2];
00423 
00424       tmp2[len - 4] = 0
00425         + c.nm[1] * ima.element(o_finish_dd)
00426         + c.nm[2] * ima.element(o_finish_d)
00427         + c.nm[3] * ima.element(o_finish)
00428         - c.dm[1] * tmp2[len-3]
00429         - c.dm[2] * tmp2[len-2];
00430 
00431       o_current = o_finish - 3 * delta_offset;
00432       o_current_d   = o_current +     delta_offset;
00433       o_current_dd  = o_current + 2 * delta_offset;
00434       o_current_ddd = o_current + 3 * delta_offset;
00435 
00436       for (int i = len - 5; i >= 0; --i)
00437         {
00438           tmp2[i] = 0
00439             + c.nm[1] * ima.element(o_current)
00440             + c.nm[2] * ima.element(o_current_d)
00441             + c.nm[3] * ima.element(o_current_dd)
00442             + c.nm[4] * ima.element(o_current_ddd)
00443             - c.dm[1] * tmp2[i+1] - c.dm[2] * tmp2[i+2]
00444             - c.dm[3] * tmp2[i+3] - c.dm[4] * tmp2[i+4];
00445           o_current     -= delta_offset;
00446           o_current_d   -= delta_offset;
00447           o_current_dd  -= delta_offset;
00448           o_current_ddd -= delta_offset;
00449         }
00450 
00451       // Combine results from causal and non-causal parts.
00452 
00453       o_current = o_start;
00454       for (int i = 0; i < len; ++i)
00455         {
00456           ima.element(o_current) = static_cast<mln_value(I)>(tmp1[i] + tmp2[i]);
00457           o_current += delta_offset;
00458         }
00459     }
00460 
00461 
00462     template <typename I>
00463     inline
00464     mln_concrete(I)
00465     gaussian_directional_2d(const Image<I>& input_,
00466                             unsigned dir, double sigma,
00467                             const mln_value(I)& bdr)
00468     {
00469       trace::entering("linear::gaussian_directional_2d");
00470 
00471       typedef mln_site(I) P;
00472       mlc_bool(P::dim == 2)::check();
00473 
00474       const I& input = exact(input_);
00475 
00476       mln_precondition(dir == 0 || dir == 1);
00477       mln_precondition(input.is_valid());
00478 
00479       my::recursivefilter_coef_ coef(1.68f, 3.735f,
00480                                      1.783f, 1.723f,
00481                                      -0.6803f, -0.2598f,
00482                                      0.6318f, 1.997f,
00483                                      sigma,
00484                                      my::recursivefilter_coef_::DericheGaussian);
00485 
00486       extension::adjust_fill(input, 5 * int(sigma + .50001) + 1, bdr);
00487       mln_concrete(I) output = duplicate(input);
00488 
00489       if (sigma < 0.006)
00490         return output;
00491 
00492       int
00493         nrows = geom::nrows(input),
00494         ncols = geom::ncols(input),
00495         b     = input.border();
00496 
00497       if (dir == 0)
00498         {
00499           for (int j = 0; j < ncols; ++j)
00500             recursivefilter_directional_fastest(output, coef,
00501                                                 point2d(- b, j),
00502                                                 point2d(nrows - 1 + b, j),
00503                                                 nrows + 2 * b,
00504                                                 dpoint2d(1, 0),
00505                                                 bdr);
00506         }
00507 
00508       if (dir == 1)
00509         {
00510           for (int i = 0; i < nrows; ++i)
00511             recursivefilter_directional_fastest(output, coef,
00512                                                 point2d(i, - b),
00513                                                 point2d(i, ncols - 1 + b),
00514                                                 ncols + 2 * b,
00515                                                 dpoint2d(0, 1),
00516                                                 bdr);
00517         }
00518 
00519       trace::exiting("linear::gaussian_directional_2d");
00520       return output;
00521     }
00522 
00523 # endif // ! MLN_INCLUDE_ONLY
00524 
00525   } // end of namespace mln::linear
00526 
00527 } // end of namespace mln
00528 
00529 
00530 #endif // ! MLN_LINEAR_GAUSSIAN_DIRECTIONAL_2D_HH
gaussian_directional_2d.hh
