Line data Source code
1 : /** \file gaussian.hpp
2 : * \author Jared R. Males
3 : * \brief Declarations for utilities related to the Gaussian function.
4 : * \ingroup gen_math_files
5 : *
6 : */
7 :
8 : //***********************************************************************//
9 : // Copyright 2015, 2016, 2017, 2018, 2019, 2020 Jared R. Males (jaredmales@gmail.com)
10 : //
11 : // This file is part of mxlib.
12 : //
13 : // mxlib is free software: you can redistribute it and/or modify
14 : // it under the terms of the GNU General Public License as published by
15 : // the Free Software Foundation, either version 3 of the License, or
16 : // (at your option) any later version.
17 : //
18 : // mxlib is distributed in the hope that it will be useful,
19 : // but WITHOUT ANY WARRANTY; without even the implied warranty of
20 : // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
21 : // GNU General Public License for more details.
22 : //
23 : // You should have received a copy of the GNU General Public License
24 : // along with mxlib. If not, see <http://www.gnu.org/licenses/>.
25 : //***********************************************************************//
26 :
27 : #ifndef gaussian_hpp
28 : #define gaussian_hpp
29 :
30 : #include <cmath>
31 :
32 : namespace mx
33 : {
34 : namespace math
35 : {
36 : namespace func
37 : {
38 :
39 : /// Constant to convert between the Gaussian width parameter and FWHM
40 : /** Used for
41 : *
42 : * \f$ FWHM = 2\sqrt{2\log2}\sigma \f$
43 : *
44 : * This was calculated in long double precision.
45 : *
46 : * \ingroup gen_math_gaussians
47 : */
48 : template <typename floatT>
49 0 : constexpr floatT twosqrt2log2()
50 : {
51 0 : return static_cast<floatT>( 2.3548200450309493820231386529193992754947713787716 );
52 : }
53 :
54 : /// Convert from FWHM to the Gaussian width parameter
55 : /** Performs the conversion:
56 : *
57 : * \f$ \sigma = 2\sqrt{2\log2}FWHM \f$
58 : *
59 : * \returns the converted value of the Gaussian width parameter
60 : *
61 : * \ingroup gen_math_gaussians
62 : */
63 : template <typename floatT>
64 : floatT fwhm2sigma( floatT fw /**< [in] the full-width at half maximum */ )
65 : {
66 : return fw / twosqrt2log2<floatT>();
67 : }
68 :
69 : /// Convert from Gaussian width parameter to FWHM
70 : /** Performs the conversion:
71 : *
72 : * \f$ FWHM = 2\sqrt{2\log2}\sigma \f$
73 : *
74 : * \returns the converted value of the full-width at half maximum
75 : *
76 : * \ingroup gen_math_gaussians
77 : */
78 : template <typename floatT>
79 0 : floatT sigma2fwhm( floatT sig /**< [in] the Gaussian width parameter */ )
80 : {
81 0 : return sig * twosqrt2log2<floatT>();
82 : }
83 :
84 : /// Find value at position (x) of the 1D arbitrarily-centered symmetric Gaussian
85 : /**
86 : * Computes:
87 : * \f$ G(x) = G_0 + G\exp[-(0.5/\sigma^2)((x-x_0)^2)]\f$
88 : *
89 : *
90 : * \returns the value of the 1D arbitrarily-centered symmetric Gaussian at (x)
91 : *
92 : * \tparam realT is the type to use for arithmetic
93 : *
94 : * \ingroup gen_math_gaussians
95 : *
96 : *
97 : */
98 : template <typename realT>
99 0 : realT gaussian( const realT x, ///< [in] is the x-position at which to evaluate the Gaussian
100 : const realT G0, ///< [in] is the constant to add to the Gaussian
101 : const realT G, ///< [in] is the scaling factor (peak height = G-G0)
102 : const realT x0, ///< [in] is the x-coordinate of the center
103 : const realT sigma ///< [in] is the width of the Gaussian.
104 : )
105 : {
106 0 : return G0 + G * exp( -( static_cast<realT>( 0.5 ) / ( sigma * sigma ) * ( ( ( x - x0 ) * ( x - x0 ) ) ) ) );
107 : }
108 :
109 : /// Find value at position (x,y) of the 2D arbitrarily-centered symmetric Gaussian
110 : /**
111 : * Computes:
112 : * \f[
113 : f(x,y) = G_0 + G \exp[-(0.5/\sigma^2)((x-x_0)^2+(y-y_0)^2)]
114 : * \f]
115 : *
116 : * \returns the value of the 2D arbitrarily-centered symmetric Gaussian at (x,y)
117 : *
118 : * \tparam realT is the type to use for arithmetic
119 : *
120 : * \ingroup gen_math_gaussians
121 : *
122 : "[test doc]"
123 : *
124 : */
125 : template <typename realT>
126 620134 : realT gaussian2D( const realT x, ///< [in] the x-position at which to evaluate the Gaussian
127 : const realT y, ///< [in] the y-positoin at which to evaluate the Gaussian
128 : const realT G0, ///< [in] the constant to add to the Gaussian
129 : const realT G, ///< [in] the scaling factor (peak height is G-G0)
130 : const realT x0, ///< [in] the x-coordinate of the center
131 : const realT y0, ///< [in] the y-coordinate of the center
132 : const realT sigma ///< [in] the width of the Gaussian.
133 : )
134 : {
135 0 : return G0 +
136 620134 : G * std::exp( -( 0.5 / ( sigma * sigma ) * ( ( ( x - x0 ) * ( x - x0 ) ) + ( ( y - y0 ) * ( y - y0 ) ) ) ) );
137 : }
138 :
139 : /// Fill in an array with the 2D arbitrarily-centered symmetric Gaussian
140 : /**
141 : * At each pixel (x,y) of the array this computes:
142 : *
143 : * \f$ f(x,y) = G_0 + G\exp[-(0.5/\sigma^2)((x-x_0)^2+(y-y_0)^2)] \f$
144 : *
145 : * \tparam realT is the type to use for arithmetic
146 : *
147 : * \ingroup gen_math_gaussians
148 : */
149 : template <typename realT>
150 47 : void gaussian2D( realT *arr, ///< [out] is the allocated array to fill in
151 : size_t nx, ///< [in] is the size of the x dimension of the array
152 : size_t ny, ///< [in] is the size of the y dimension of the array
153 : const realT G0, ///< [in] is the constant to add to the Gaussian
154 : const realT G, ///< [in] is the scaling factor (peak height = G-G0)
155 : const realT x0, ///< [in] is the x-coordinate of the center
156 : const realT y0, ///< [in] is the y-coordinate of the center
157 : const realT sigma ///< [in] is the third rotation and scaling factor
158 : )
159 : {
160 : size_t idx;
161 :
162 4357 : for( size_t j = 0; j < ny; ++j )
163 : {
164 624444 : for( size_t i = 0; i < nx; ++i )
165 : {
166 620134 : idx = i + j * nx;
167 :
168 620134 : arr[idx] = gaussian2D( (realT)i, (realT)j, G0, G, x0, y0, sigma );
169 : }
170 : }
171 47 : }
172 :
173 : /// Find value at position (x,y) of the 2D general elliptical Gaussian
174 : /**
175 : * Computes:
176 : *
177 : * \f$ f(x,y) = G_0 + G\exp [-0.5( a(x-x_0)^2 + b(x-x_0)(y-y_0) + c(y-y_0)^2 ]\f$
178 : *
179 : * where, for a counter-clockwise rotation \f$ \theta \f$, we have
180 : *
181 : * \f$ a = \frac{\cos^2 \theta}{\sigma_x^2} + \frac{\sin^2\theta}{\sigma_y^2}\f$
182 : *
183 : * \f$ b = \frac{\sin 2\theta}{2} \left(\frac{1}{\sigma_x^2} - \frac{1}{\sigma_y^2} \right)\f$
184 : *
185 : * \f$ c = \frac{\sin^2 \theta}{\sigma_x^2} + \frac{\cos^2\theta}{\sigma_y^2}\f$
186 : *
187 : * In this version the parameters are specified directly as (a,b,c), in particular avoiding the trig function calls.
188 : * This should be much more efficient, and so this version should be used inside fitting routines, etc. However
189 : * note that the matrix {a b}{b c} must be <a
190 : * href="http://mathworld.wolfram.com/PositiveDefiniteMatrix.html">positive-definite</a> otherwise infinities can result
191 : * from the argument of the exponent being positive.
192 : *
193 : * The functions gaussian2D_gen2rot() and gaussian2D_rot2gen() provide conversions from
194 : * (a,b,c) to (\f$\sigma_x\f$,\f$\sigma_y\f$,\f$\theta\f$) and back. The function gaussian2D_ang() is a
195 : * wrapper for this function, which instead accepts (\f$\sigma_x\f$,\f$\sigma_y\f$,\f$\theta\f$) as inputs.
196 : *
197 : * \returns the value of the 2D elliptical Gaussian at (x,y)
198 : *
199 : * \tparam realT is the type to use for arithmetic
200 : *
201 : * \ingroup gen_math_gaussians
202 : */
203 : template <typename realT>
204 298237 : realT gaussian2D( const realT x, ///< [in] the x-position at which to evaluate the Gaussian
205 : const realT y, ///< [in] the y-positoin at which to evaluate the Gaussian
206 : const realT G0, ///< [in] the constant to add to the Gaussian
207 : const realT G, ///< [in] the scaling factor (peak = G-G0)
208 : const realT x0, ///< [in] the x-coordinate of the center
209 : const realT y0, ///< [in] the y-coordinate of the center
210 : const realT a, ///< [in] the first rotation and scaling factor
211 : const realT b, ///< [in] the second rotation and scaling factor
212 : const realT c ///< [in] the third rotation and scaling factor
213 : )
214 : {
215 298237 : realT dx = x - x0;
216 298237 : realT dy = y - y0;
217 :
218 298237 : return G0 + G * std::exp( -0.5 * ( a * dx * dx + 2. * b * dx * dy + c * dy * dy ) );
219 : }
220 :
221 : /// Convert from (a,b,c) to (\f$\sigma_x\f$,\f$\sigma_y\f$,\f$\theta\f$) for the elliptical Gaussian.
222 : /** The general 2D elliptical Gaussian
223 : *
224 : * \f$ f(x,y) = G_0 + G\exp [-0.5( a(x-x_0)^2 + b(x-x_0)(y-y_0) + c(y-y_0)^2 ]\f$
225 : *
226 : * can be expressed in terms of widths and a rotation angle using:
227 : *
228 : * \f$ a = \frac{\cos^2 \theta}{\sigma_x^2} + \frac{\sin^2\theta}{\sigma_y^2}\f$
229 : *
230 : * \f$ b = \frac{\sin 2\theta}{2} \left(\frac{1}{\sigma_x^2} - \frac{1}{\sigma_y^2} \right)\f$
231 : *
232 : * \f$ c = \frac{\sin^2 \theta}{\sigma_x^2} + \frac{\cos^2\theta}{\sigma_y^2}\f$
233 : *
234 : * where \f$ \theta \f$ specifies a counter-clockwise rotation.
235 : *
236 : * This function calculates (\f$\sigma_x\f$,\f$\sigma_y\f$,\f$\theta\f$) from inputs (a,b,c). It adopts
237 : * the convention that the long axis of the ellipse is \f$\sigma_x\f$, and \f$\theta\f$ is chosen
238 : * appropriately. Note that \f$-\frac{\pi}{2} < \theta < \frac{\pi}{2}\f$.
239 : *
240 : * \tparam realT is the type to use for arithmetic
241 : *
242 : * \ingroup gen_math_gaussians
243 : */
244 : template <typename realT>
245 16 : void gaussian2D_gen2rot( realT &sigma_x, ///< [out] the width parameter in the rotated x-direction
246 : realT &sigma_y, ///< [out] the width parameter in the rotated y-direction
247 : realT &theta, ///< [out] the c.c.w rotation
248 : const realT a, ///< [in] the first rotation and scaling parameter
249 : const realT b, ///< [in] the second rotation and scaling parameter
250 : const realT c ///< [in] the the third rotation and scaling parameter
251 : )
252 : {
253 : realT x1, x2, s, s2, theta0, theta1;
254 :
255 16 : realT arg = a * a - 2 * a * c + 4 * b * b + c * c;
256 16 : if( arg < 0 )
257 : {
258 0 : x2 = 0.5 * ( a + c );
259 0 : x1 = x2;
260 : }
261 : else
262 : {
263 : // There are two roots, one is 1/sigma_x^2 and one is 1/sigma_y^2
264 16 : x2 = 0.5 * ( a + c ) + 0.5 * sqrt( arg );
265 16 : x1 = 0.5 * ( a + c ) - 0.5 * sqrt( arg );
266 : }
267 :
268 : // Our convention is that the larger width is sigma_x
269 16 : sigma_x = sqrt( 1. / x1 );
270 16 : sigma_y = sqrt( 1. / x2 );
271 :
272 : //----------------------------------------------------------------------------//
273 : // Choosing theta:
274 : // There is an ambiguity in theta and which direction (x,y) is the long axis
275 : // Here we choose theta so that it specifies the long direction, sigma_x
276 : //----------------------------------------------------------------------------//
277 :
278 : // s is (sin(theta))^2
279 : // s2 = (a-x2)*(a-x2)/(b*b + (a-x2)*(a-x2));
280 16 : s = ( a - x1 ) * ( a - x1 ) / ( b * b + ( a - x1 ) * ( a - x1 ) );
281 :
282 : // First check if x1-x2 will be close to zero
283 :
284 16 : if( fabs( x1 - x2 ) < 1e-12 )
285 : {
286 16 : theta = 0;
287 16 : return;
288 : }
289 :
290 0 : theta0 = 0.5 * asin( 2 * b / ( x1 - x2 ) ); // This always gives the correct sign
291 0 : theta1 = asin( sqrt( s ) ); // This always gives the correct magnitude, and seems to be more accurate
292 :
293 : // Compare signs. If they match, then we use theta1
294 0 : if( std::signbit( theta0 ) == std::signbit( theta1 ) )
295 : {
296 0 : theta = theta1;
297 : }
298 : else
299 : {
300 : // Otherwise, we switch quadrants
301 0 : theta = asin( -sqrt( s ) );
302 : }
303 : }
304 :
305 : /// Convert from (\f$\sigma_x\f$,\f$\sigma_y\f$,\f$\theta\f$) to (a,b,c) for the elliptical Gaussian.
306 : /** The general 2D elliptical Gaussian
307 : *
308 : * \f$ f(x,y) = G_0 + G\exp [-0.5( a(x-x_0)^2 + b(x-x_0)(y-y_0) + c(y-y_0)^2 ]\f$
309 : *
310 : * can be expressed in terms of widths and a rotation angle using:
311 : *
312 : * \f$ a = \frac{\cos^2 \theta}{\sigma_x^2} + \frac{\sin^2\theta}{\sigma_y^2}\f$
313 : *
314 : * \f$ b = \frac{\sin 2\theta}{2} \left(\frac{1}{\sigma_x^2} - \frac{1}{\sigma_y^2} \right)\f$
315 : *
316 : * \f$ c = \frac{\sin^2 \theta}{\sigma_x^2} + \frac{\cos^2\theta}{\sigma_y^2}\f$
317 : *
318 : * where \f$ \theta \f$ specifies a counter-clockwise rotation.
319 : *
320 : * This function calculates (a,b,c) from inputs (\f$\sigma_x\f$,\f$\sigma_y\f$,\f$\theta\f$).
321 : *
322 : * \tparam realT is the type to use for arithmetic
323 : *
324 : * \ingroup gen_math_gaussians
325 : */
326 : template <typename realT>
327 8 : void gaussian2D_rot2gen( realT &a, ///< [out] the first rotation and scaling parameter
328 : realT &b, ///< [out] the second rotation and scaling parameter
329 : realT &c, ///< [out] the the third rotation and scaling parameter
330 : const realT sigma_x, ///< [in] the width parameter in the rotated x-direction
331 : const realT sigma_y, ///< [in] the width parameter in the rotated y-direction
332 : const realT theta ///< [in] the c.c.w rotation
333 : )
334 : {
335 : realT sn, cs, sx2, sy2;
336 :
337 8 : sn = sin( theta );
338 8 : cs = cos( theta );
339 8 : sx2 = sigma_x * sigma_x;
340 8 : sy2 = sigma_y * sigma_y;
341 :
342 8 : a = cs * cs / sx2 + sn * sn / sy2;
343 8 : b = sn * cs * ( 1. / sx2 - 1. / sy2 );
344 8 : c = sn * sn / sx2 + cs * cs / sy2;
345 8 : }
346 :
347 : /// Find value at position (x,y) of the 2D rotated elliptical Gaussian
348 : /**
349 : * Computes:
350 : * \f$ f(x,y) = G_0 + G\exp [-0.5( a(x-x_0)^2 + b(x-x_0)(y-y_0) + c(y-y_0)^2 ]\f$
351 : *
352 : * where, for a counter-clockwise rotation \f$ \theta \f$, we have
353 : *
354 : * \f$ a = \frac{\cos^2 \theta}{\sigma_x^2} + \frac{\sin^2\theta}{\sigma_y^2}\f$
355 : *
356 : * \f$ b = \frac{\sin 2\theta}{2} \left(\frac{1}{\sigma_x^2} - \frac{1}{\sigma_y^2} \right)\f$
357 : *
358 : * \f$ c = \frac{\sin^2 \theta}{\sigma_x^2} + \frac{\cos^2\theta}{\sigma_y^2}\f$
359 : *
360 : * This is a convenience wrapper for the general elliptical Gaussian function gaussian2D(), where here
361 : * (a,b,c) are first calculated from (\f$\sigma_x\f$,\f$\sigma_y\f$,\f$\theta\f$). This will in general be
362 : * slower due to the trig function calls, so the (a,b,c) version should be used most of the time.
363 : *
364 : * The functions gaussian2D_gen2rot() and gaussian2D_rot2gen() provide conversions from
365 : * (a,b,c) to (\f$\sigma_x\f$,\f$\sigma_y\f$,\f$\theta\f$) and back. The function gaussian2D_rot() is a
366 : * wrapper for this function, which instead accepts (\f$\sigma_x\f$,\f$\sigma_y\f$,\f$\theta\f$) as inputs.
367 : *
368 : * \returns the value of the 2D elliptical Gaussian at (x,y)
369 : *
370 : * \tparam realT is the type to use for arithmetic
371 : *
372 : * \ingroup gen_math_gaussians
373 : */
374 : template <typename realT>
375 : realT gaussian2D_ang( const realT x, ///< [in] the x-position at which to evaluate the Gaussian
376 : const realT y, ///< [in] the y-positoin at which to evaluate the Gaussian
377 : const realT G0, ///< [in] the constant to add to the Gaussian
378 : const realT G, ///< [in] the scaling factor (peak height = G-G0)
379 : const realT x0, ///< [in] the x-coordinate of the center
380 : const realT y0, ///< [in] the y-coordinate of the center
381 : const realT sigma_x, ///< [in] the width in the rotated x direction
382 : const realT sigma_y, ///< [in] the width in the rotated y direction
383 : const realT theta ///< [in] the counter-clockwise rotation angle
384 : )
385 : {
386 : realT a, b, c;
387 :
388 : gaussian2D_rot2gen( a, b, c, sigma_x, sigma_y, theta );
389 :
390 : return gaussian2D( x, y, G0, G, x0, y0, a, b, c );
391 : }
392 :
393 : /// Fill in an array with the 2D general elliptical Gaussian
394 : /**
395 : * At each pixel (x,y) of the array this computes:
396 : *
397 : * \f$ f(x,y) = G_0 + G\exp [-0.5( a(x-x_0)^2 + b(x-x_0)(y-y_0) + c(y-y_0)^2 ]\f$
398 : *
399 : * where, for a counter-clockwise rotation \f$ \theta \f$, we have
400 : *
401 : * \f$ a = \frac{\cos^2 \theta}{\sigma_x^2} + \frac{\sin^2\theta}{\sigma_y^2}\f$
402 : *
403 : * \f$ b = \frac{\sin 2\theta}{2} \left(\frac{1}{\sigma_x^2} - \frac{1}{\sigma_y^2} \right)\f$
404 : *
405 : * \f$ c = \frac{\sin^2 \theta}{\sigma_x^2} + \frac{\cos^2\theta}{\sigma_y^2}\f$
406 : *
407 : * In this version the parameters are specified directly as (a,b,c), in particular avoiding the trig function calls.
408 : * This should be much more efficient, and so this version should be used inside fitting routines, etc.
409 : *
410 : * The functions gaussian2D_gen2rot() and gaussian2D_rot2gen() provide conversions from
411 : * (a,b,c) to (\f$\sigma_x\f$,\f$\sigma_y\f$,\f$\theta\f$) and back. The function gaussian2D_ang() is a
412 : * wrapper for this function, which instead accepts (\f$\sigma_x\f$,\f$\sigma_y\f$,\f$\theta\f$) as inputs.
413 : *
414 : * \tparam realT is the type to use for arithmetic
415 : *
416 : * \ingroup gen_math_gaussians
417 : */
418 : template <typename realT>
419 : void gaussian2D( realT *arr, ///< [out] the native array to populate with the Gaussian function
420 : size_t nx, ///< [in] the x-size of the array, in pixels
421 : size_t ny, ///< [in] the y-size of the array, in pixels
422 : const realT G0, ///< [in] the constant to add to the Gaussian
423 : const realT G, ///< [in] the scaling factor (peak = G)
424 : const realT x0, ///< [in] the x-coordinate of the center, in pixels
425 : const realT y0, ///< [in] the y-coordinate of the center, in pixels
426 : const realT a, ///< [in] is the first rotation and scaling factor
427 : const realT b, ///< [in] is the second rotation and scaling factor
428 : const realT c ///< [in] is the third rotation and scaling factor
429 : )
430 : {
431 : size_t idx;
432 :
433 : for( size_t j = 0; j < ny; ++j )
434 : {
435 : for( size_t i = 0; i < nx; ++i )
436 : {
437 : idx = i + j * nx;
438 :
439 : arr[idx] = gaussian2D( (realT)i, (realT)j, G0, G, x0, y0, a, b, c );
440 : }
441 : }
442 : }
443 :
444 : /// Fill in an array with the 2D general elliptical Gaussian
445 : /**
446 : * At each pixel (x,y) of the array this computes:
447 : *
448 : * \f$ f(x,y) = G_0 + G\exp [-0.5( a(x-x_0)^2 + b(x-x_0)(y-y_0) + c(y-y_0)^2 ]\f$
449 : *
450 : * where, for a counter-clockwise rotation \f$ \theta \f$, we have
451 : *
452 : * \f$ a = \frac{\cos^2 \theta}{\sigma_x^2} + \frac{\sin^2\theta}{\sigma_y^2}\f$
453 : *
454 : * \f$ b = \frac{\sin 2\theta}{2} \left(\frac{1}{\sigma_x^2} - \frac{1}{\sigma_y^2} \right)\f$
455 : *
456 : * \f$ c = \frac{\sin^2 \theta}{\sigma_x^2} + \frac{\cos^2\theta}{\sigma_y^2}\f$
457 : *
458 : * This is a convenience wrapper for the general elliptical Gaussian function
459 : * gaussian2D( realT *, size_t, size_t, const realT, const realT, const realT, const realT, const realT, const realT,
460 : * const realT) , where here (a,b,c) are first calculated from (\f$\sigma_x\f$,\f$\sigma_y\f$,\f$\theta\f$).
461 : *
462 : * The functions gaussian2D_gen2rot() and gaussian2D_rot2gen() provide conversions from
463 : * (a,b,c) to (\f$\sigma_x\f$,\f$\sigma_y\f$,\f$\theta\f$) and back.
464 : *
465 : *
466 : * \tparam realT is the type to use for arithmetic
467 : *
468 : * \ingroup gen_math_gaussians
469 : */
470 : template <typename realT>
471 : void gaussian2D_ang( realT *arr, ///< [out] the native array to populate with the Gaussian function
472 : size_t nx, ///< [in] the x-size of the array, in pixels
473 : size_t ny, ///< [in] the y-size of the array, in pixels
474 : const realT G0, ///< [in] the constant to add to the Gaussian
475 : const realT G, ///< [in] the scaling factor (peak = G)
476 : const realT x0, ///< [in] the x-coordinate of the center, in pixels
477 : const realT y0, ///< [in] the y-coordinate of the center, in pixels
478 : const realT sigma_x, ///< [in] the width in the rotated x direction, in pixels
479 : const realT sigma_y, ///< [in] the width in the rotated y direction, in pixels
480 : const realT theta ///< [in] the counter-clockwise rotation angle, in radians
481 : )
482 : {
483 : realT a, b, c;
484 : size_t idx;
485 :
486 : // Get the (a,b,c) parameters so the trig only happens once
487 : gaussian2D_rot2gen( a, b, c, sigma_x, sigma_y, theta );
488 :
489 : gaussian2D( arr, nx, ny, G0, G, x0, y0, a, b, c );
490 : }
491 :
492 : /// Calculate the Jacobian at position (x,y) for the 2D general elliptical Gaussian
493 : /** \note this has not been verified and may be incorrect.
494 : *
495 : * Given:
496 : *
497 : * \f$ f(x,y) = G_0 + G\exp [-0.5( a(x-x_0)^2 + b(x-x_0)(y-y_0) + c(y-y_0)^2 ]\f$
498 : *
499 : * we calculate
500 : *
501 : * \f$ \frac{\partial G}{\partial G_0} = 1 \f$
502 : *
503 : * \f$ \frac{\partial G}{\partial G} = (G - G_0)/A \f$
504 : *
505 : * \f$ \frac{\partial G}{\partial x_0} = 0.5(G-G_0)( 2a(x-x_0) + b(y-y_0)) \f$
506 : *
507 : * \f$ \frac{\partial G}{\partial y_0} = 0.5(G-G_0) ( b(x-x_0) + 2c(y-y_0)) \f$
508 : *
509 : * \f$ \frac{\partial G}{\partial a} = -0.5(G-G_0) ( (x-x_0)^2 ) \f$
510 : *
511 : * \f$ \frac{\partial G}{\partial b} = -0.5(G-G_0) ( (x-x_0)(y-y_0) ) \f$
512 : *
513 : * \f$ \frac{\partial G}{\partial c} = -0.5(G-G_0) ( (y-y_0)^2 )\f$
514 : *
515 : * \param j is a 7 element vector which is populated with the derivatives
516 : * \param x is the x-position at which to evaluate the Gaussian
517 : * \param y is the y-positoin at which to evaluate the Gaussian
518 : * \param G0 is the constant to add to the Gaussian
519 : * \param G is the scaling factor (peak = A)
520 : * \param x0 is the x-coordinate of the center
521 : * \param y0 is the y-coordinate of the center
522 : * \param a is the first rotation and scaling factor
523 : * \param b is the second rotation and scaling factor
524 : * \param c is the third rotation and scaling factor
525 : *
526 : *
527 : * \tparam realT is the type to use for arithmetic
528 : *
529 : * \ingroup gen_math_gaussians
530 : */
531 : template <typename realT>
532 : void gaussian2D_jacobian( realT *j,
533 : const realT x,
534 : const realT y,
535 : const realT G0,
536 : const realT G,
537 : const realT x0,
538 : const realT y0,
539 : const realT a,
540 : const realT b,
541 : const realT c )
542 : {
543 : realT G_G0 = gaussian2D<realT>( x, y, 0, G, x0, y0, a, b, c );
544 :
545 : j[0] = (realT)1;
546 :
547 : j[1] = ( G_G0 ) / G;
548 :
549 : j[2] = 0.5 * G_G0 * ( 2 * a * ( x - x0 ) + b * ( y - y0 ) );
550 :
551 : j[3] = 0.5 * G_G0 * ( b * ( x - x0 ) + 2 * c * ( y - y0 ) );
552 :
553 : j[4] = -0.5 * G_G0 * ( ( x - x0 ) * ( x - x0 ) );
554 :
555 : j[5] = -0.5 * G_G0 * ( ( x - x0 ) * ( y - y0 ) );
556 :
557 : j[6] = -0.5 * G_G0 * ( ( y - y0 ) * ( y - y0 ) );
558 : }
559 :
560 : } // namespace func
561 : } // namespace math
562 : } // namespace mx
563 :
564 : #endif // gaussian_hpp
|
