Milena (Olena): coefficients.hh Source File

00001 // Copyright (C) 2001, 2002, 2003, 2004, 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_INTERNAL_COEFFICIENTS_HH
00027 # define MLN_LINEAR_GAUSSIAN_INTERNAL_COEFFICIENTS_HH
00028 
00035 
00036 # include <vector>
00037 # include <cmath>
00038 
00039 
00040 
00041 namespace mln
00042 {
00043 
00044   namespace linear
00045   {
00046 
00047     namespace gaussian
00048     {
00049 
00050       namespace internal
00051       {
00052 
00053 
00054         typedef double norm_fun(double, double,
00055                                 double, double,
00056                                 double, double,
00057                                 double, double,
00058                                 double, double,
00059                                 int&);
00060 
00061         struct coefficients
00062         {
00064           //
00065           coefficients(double a0, double a1,
00066                       double b0, double b1,
00067                       double c0, double c1,
00068                       double w0, double w1,
00069                       double s, norm_fun norm);
00070           std::vector<double> n, d, nm, dm;
00071         };
00072 
00073 
00074         coefficients coefficients_not_derivative(double sigma);
00075         coefficients coefficients_1st_derivative(double sigma);
00076         coefficients coefficients_2nd_derivative(double sigma);
00077 
00078 
00079 
00080 # ifndef MLN_INCLUDE_ONLY
00081 
00082 
00083         inline
00084         coefficients::coefficients(double a0, double a1,
00085                                  double b0, double b1,
00086                                  double c0, double c1,
00087                                  double w0, double w1,
00088                                  double s, norm_fun norm)
00089         {
00090           n.reserve(5);
00091           d.reserve(5);
00092           nm.reserve(5);
00093           dm.reserve(5);
00094 
00095           b0 /= s;
00096           b1 /= s;
00097           w0 /= s;
00098           w1 /= s;
00099 
00100           double sin0 = std::sin(w0);
00101           double sin1 = std::sin(w1);
00102           double cos0 = std::cos(w0);
00103           double cos1 = std::cos(w1);
00104 
00105           int sign = 1;
00106           double n_ = norm(a0, a1, b0, b1, c0, c1, cos0, sin0, cos1, sin1, sign);
00107 
00108           a0 /= n_;
00109           a1 /= n_;
00110           c0 /= n_;
00111           c1 /= n_;
00112 
00113           n[3] =
00114             std::exp(-b1 - 2*b0) * (c1 * sin1 - cos1 * c0) +
00115             std::exp(-b0 - 2*b1) * (a1 * sin0 - cos0 * a0);
00116           n[2] =
00117             2 * std::exp(-b0 - b1) * ((a0 + c0) * cos1 * cos0 -
00118                                       cos1 * a1 * sin0 -
00119                                       cos0 * c1 * sin1) +
00120             c0 * std::exp(-2*b0) + a0 * std::exp(-2*b1);
00121           n[1] =
00122             std::exp(-b1) * (c1 * sin1 - (c0 + 2 * a0) * cos1) +
00123             std::exp(-b0) * (a1 * sin0 - (2 * c0 + a0) * cos0);
00124           n[0] =
00125             a0 + c0;
00126 
00127           d[4] =
00128             std::exp(-2 * b0 - 2 * b1);
00129           d[3] =
00130             -2 * cos0 * std::exp(-b0 - 2*b1) -
00131             2 * cos1 * std::exp(-b1 - 2*b0);
00132           d[2] =
00133             4 * cos1 * cos0 * std::exp(-b0 - b1) +
00134             std::exp(-2*b1) + std::exp(-2*b0);
00135           d[1] =
00136             -2 * std::exp(-b1) * cos1 - 2 * std::exp(-b0) * cos0;
00137 
00138           for (unsigned i = 1; i <= 3; ++i)
00139             {
00140               dm[i] = d[i];
00141               nm[i] = double(sign) * (n[i] - d[i] * n[0]);
00142             }
00143           dm[4] = d[4];
00144           nm[4] = double(sign) * (-d[4] * n[0]);
00145         }
00146 
00147 
00148 
00149         inline
00150         double norm_not_derivative(double a0, double a1,
00151                                    double b0, double b1,
00152                                    double c0, double c1,
00153                                    double cos0, double sin0,
00154                                    double cos1, double sin1,
00155                                    int& sign)
00156         {
00157           double expb0 = std::exp(b0);
00158           double exp2b0 = std::exp(2.f * b0);
00159 
00160           double scale0 = 1 + exp2b0 - 2 * cos0 * expb0;
00161           double scaleA = 2 * a1 * sin0 * expb0 - a0 * (1 - exp2b0);
00162 
00163           double expb1 = std::exp(b1);
00164           double exp2b1 = std::exp(2.f * b1);
00165 
00166           double scale1 = 1 + exp2b1 - 2 * cos1 * expb1;
00167           double scaleC = 2 * c1 * sin1 * expb1 - c0 * (1 - exp2b1);
00168 
00169           double sumA = scaleA / scale0;
00170           double sumC = scaleC / scale1;
00171 
00172           sign = 1;
00173 
00174           return sumA + sumC;
00175         }
00176 
00177 
00178         inline
00179         double norm_1st_derivative(double a0, double a1,
00180                                    double b0, double b1,
00181                                    double c0, double c1,
00182                                    double cos0, double sin0,
00183                                    double cos1, double sin1,
00184                                    int& sign)
00185         {
00186           double expb0 = std::exp(b0);
00187           double exp2b0 = std::exp(2.f * b0);
00188 
00189           double scale0 = 1 + exp2b0 - 2 * cos0 * expb0;
00190           scale0 *= scale0;
00191           double scaleA = - 2 * a1 * sin0 * expb0 * (1 - exp2b0) +
00192             2 * a0 * expb0 * (2 * expb0 - cos0 * (1 + exp2b0));
00193 
00194           double expb1 = std::exp(b1);
00195           double exp2b1 = std::exp(2.f * b1);
00196 
00197           double scale1 = 1 + exp2b1 - 2 * cos1 * expb1;
00198           scale1 *= scale1;
00199           double scaleC = - 2 * c1 * sin1 * expb1 * (1 - exp2b1) +
00200             2 * c0 * expb1 * (2 * expb1 - cos1 * (1 + exp2b1));
00201 
00202           double sumA = scaleA / scale0;
00203           double sumC = scaleC / scale1;
00204 
00205           sign = -1;
00206 
00207           return sumA + sumC;
00208         }
00209 
00210 
00211         inline
00212         double norm_2nd_derivative(double a0, double a1,
00213                                    double b0, double b1,
00214                                    double c0, double c1,
00215                                    double cos0, double sin0,
00216                                    double cos1, double sin1,
00217                                    int& sign)
00218         {
00219           double expb0 = std::exp(b0);
00220           double exp2b0 = std::exp(2.f * b0);
00221 
00222           double scale0 = 1 + exp2b0 - 2 * cos0 * expb0;
00223           scale0 *= scale0 * scale0;
00224 
00225           double scaleA = a1 * sin0 * expb0 *
00226             (1 + expb0 * (2 * cos0 * (1 + exp2b0) + exp2b0 - 6)) +
00227             a0 * expb0 * (2 * expb0 * (2 - cos0 * cos0) *
00228                           (1 - exp2b0) - cos0 * (1 - exp2b0 * exp2b0));
00229 
00230           double expb1 = std::exp(b1);
00231           double exp2b1 = std::exp(2.f * b1);
00232 
00233           double scale1 = 1 + exp2b1 - 2 * cos1 * expb1;
00234           scale1 *= scale1 * scale1;
00235 
00236           double scaleC = c1 * sin1 * expb1 *
00237             (1 + expb1 * (2 * cos1 * (1 + exp2b1) + exp2b1 - 6)) +
00238             c0 * expb1 * (2 * expb1 * (2 - cos1 * cos1) *
00239                           (1 - exp2b1) - cos1 * (1 - exp2b1 * exp2b1));
00240 
00241           double sumA = scaleA / scale0;
00242           double sumC = scaleC / scale1;
00243           sign = 1;
00244 
00245           return sumA + sumC;
00246         }
00247 
00248 
00249         inline
00250         coefficients coefficients_not_derivative(double sigma)
00251         {
00252           coefficients tmp(+1.6800, +3.7350,
00253                           +1.7830, +1.7230,
00254                           -0.6803, -0.2598,
00255                           +0.6318, +1.9970,
00256                           sigma,
00257                           norm_not_derivative);
00258           return tmp;
00259         }
00260 
00261 
00262         inline
00263         coefficients coefficients_1st_derivative(double sigma)
00264         {
00265           coefficients tmp(-0.6472, -4.5310,
00266                           +1.5270, +1.5160,
00267                           +0.6494, +0.9557,
00268                           +0.6719, +2.0720,
00269                           sigma,
00270                           norm_1st_derivative);
00271           return tmp;
00272         }
00273 
00274 
00275         inline
00276         coefficients coefficients_2nd_derivative(double sigma)
00277         {
00278           coefficients tmp(-1.3310, +3.661,
00279                           +1.2400, +1.314,
00280                           +0.3225, -1.738,
00281                           +0.7480, +2.166,
00282                           sigma,
00283                           norm_2nd_derivative);
00284           return tmp;
00285         }
00286 
00287 # endif // ! MLN_INCLUDE_ONLY
00288 
00289       } // end of namespace mln::linear::gaussian::internal
00290 
00291     } // end of namespace mln::linear::gaussian
00292 
00293   } // end of namespace mln::linear
00294 
00295 } // end of namespace mln
00296 
00297 
00298 #endif // ! MLN_LINEAR_GAUSSIAN_INTERNAL_COEFFICIENTS_HH