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_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
00179
00180
00181
00182
00183
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
00222
00223 tmp2[len - 1] = WorkType();
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
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
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
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
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
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
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
00501
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
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
00526
00527 if (sigma > 0.006)
00528 generic_filter_(mln_trait_image_dimension(I)(),
00529 work_img, coef, dir);
00530
00531
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
00550
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
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
00575
00576 if (sigma > 0.006)
00577 generic_filter_(mln_trait_image_dimension(I)(),
00578 work_img, coef, dir);
00579
00580
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
00596
00597
00598 data::paste(in, out);
00599
00600
00601
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
00619
00620
00621 data::paste(in, out);
00622
00623
00624
00625 if (sigma > 0.006)
00626 generic_filter_(mln_trait_image_dimension(I)(),
00627 out, coef, dir);
00628 }
00629
00630 }
00631
00632
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 }
00821
00822 }
00823
00824
00825 #endif // ! MLN_LINEAR_GAUSSIAN_HH