00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026 #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 }
00290
00291 }
00292
00293 }
00294
00295 }
00296
00297
00298 #endif // ! MLN_LINEAR_GAUSSIAN_INTERNAL_COEFFICIENTS_HH