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 : constexpr floatT twosqrt2log2()
50 : {
51 : 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 : floatT sigma2fwhm( floatT sig /**< [in] the Gaussian width parameter */ )
80 : {
81 : 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 : 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 : 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 Scenario: Verify direction and accuracy of various image shifts \ref tests_improc_imageTransforms_imageShift
123 : "[test doc]"
124 : *
125 : */
126 : template <typename realT>
127 566796 : realT gaussian2D( const realT x, ///< [in] the x-position at which to evaluate the Gaussian
128 : const realT y, ///< [in] the y-positoin at which to evaluate the Gaussian
129 : const realT G0, ///< [in] the constant to add to the Gaussian
130 : const realT G, ///< [in] the scaling factor (peak height is G-G0)
131 : const realT x0, ///< [in] the x-coordinate of the center
132 : const realT y0, ///< [in] the y-coordinate of the center
133 : const realT sigma ///< [in] the width of the Gaussian.
134 : )
135 : {
136 : return G0 +
137 566796 : G * std::exp( -( 0.5 / ( sigma * sigma ) * ( ( ( x - x0 ) * ( x - x0 ) ) + ( ( y - y0 ) * ( y - y0 ) ) ) ) );
138 : }
139 :
140 : /// Fill in an array with the 2D arbitrarily-centered symmetric Gaussian
141 : /**
142 : * At each pixel (x,y) of the array this computes:
143 : *
144 : * \f$ f(x,y) = G_0 + G\exp[-(0.5/\sigma^2)((x-x_0)^2+(y-y_0)^2)] \f$
145 : *
146 : * \tparam realT is the type to use for arithmetic
147 : *
148 : * \ingroup gen_math_gaussians
149 : */
150 : template <typename realT>
151 33 : void gaussian2D( realT *arr, ///< [out] is the allocated array to fill in
152 : size_t nx, ///< [in] is the size of the x dimension of the array
153 : size_t ny, ///< [in] is the size of the y dimension of the array
154 : const realT G0, ///< [in] is the constant to add to the Gaussian
155 : const realT G, ///< [in] is the scaling factor (peak height = G-G0)
156 : const realT x0, ///< [in] is the x-coordinate of the center
157 : const realT y0, ///< [in] is the y-coordinate of the center
158 : const realT sigma ///< [in] is the third rotation and scaling factor
159 : )
160 : {
161 : size_t idx;
162 :
163 3501 : for( size_t j = 0; j < ny; ++j )
164 : {
165 570264 : for( size_t i = 0; i < nx; ++i )
166 : {
167 566796 : idx = i + j * nx;
168 :
169 566796 : arr[idx] = gaussian2D( (realT)i, (realT)j, G0, G, x0, y0, sigma );
170 : }
171 : }
172 33 : }
173 :
174 : /// Find value at position (x,y) of the 2D general elliptical Gaussian
175 : /**
176 : * Computes:
177 : *
178 : * \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$
179 : *
180 : * where, for a counter-clockwise rotation \f$ \theta \f$, we have
181 : *
182 : * \f$ a = \frac{\cos^2 \theta}{\sigma_x^2} + \frac{\sin^2\theta}{\sigma_y^2}\f$
183 : *
184 : * \f$ b = \frac{\sin 2\theta}{2} \left(\frac{1}{\sigma_x^2} - \frac{1}{\sigma_y^2} \right)\f$
185 : *
186 : * \f$ c = \frac{\sin^2 \theta}{\sigma_x^2} + \frac{\cos^2\theta}{\sigma_y^2}\f$
187 : *
188 : * In this version the parameters are specified directly as (a,b,c), in particular avoiding the trig function calls.
189 : * This should be much more efficient, and so this version should be used inside fitting routines, etc. However
190 : * note that the matrix {a b}{b c} must be <a
191 : * href="http://mathworld.wolfram.com/PositiveDefiniteMatrix.html">positive-definite</a> otherwise infinities can result
192 : * from the argument of the exponent being positive.
193 : *
194 : * The functions gaussian2D_gen2rot() and gaussian2D_rot2gen() provide conversions from
195 : * (a,b,c) to (\f$\sigma_x\f$,\f$\sigma_y\f$,\f$\theta\f$) and back. The function gaussian2D_ang() is a
196 : * wrapper for this function, which instead accepts (\f$\sigma_x\f$,\f$\sigma_y\f$,\f$\theta\f$) as inputs.
197 : *
198 : * \returns the value of the 2D elliptical Gaussian at (x,y)
199 : *
200 : * \tparam realT is the type to use for arithmetic
201 : *
202 : * \ingroup gen_math_gaussians
203 : */
204 : template <typename realT>
205 : realT gaussian2D( const realT x, ///< [in] the x-position at which to evaluate the Gaussian
206 : const realT y, ///< [in] the y-positoin at which to evaluate the Gaussian
207 : const realT G0, ///< [in] the constant to add to the Gaussian
208 : const realT G, ///< [in] the scaling factor (peak = G-G0)
209 : const realT x0, ///< [in] the x-coordinate of the center
210 : const realT y0, ///< [in] the y-coordinate of the center
211 : const realT a, ///< [in] the first rotation and scaling factor
212 : const realT b, ///< [in] the second rotation and scaling factor
213 : const realT c ///< [in] the third rotation and scaling factor
214 : )
215 : {
216 : realT dx = x - x0;
217 : realT dy = y - y0;
218 :
219 : return G0 + G * std::exp( -0.5 * ( a * dx * dx + 2. * b * dx * dy + c * dy * dy ) );
220 : }
221 :
222 : /// Convert from (a,b,c) to (\f$\sigma_x\f$,\f$\sigma_y\f$,\f$\theta\f$) for the elliptical Gaussian.
223 : /** The general 2D elliptical Gaussian
224 : *
225 : * \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$
226 : *
227 : * can be expressed in terms of widths and a rotation angle using:
228 : *
229 : * \f$ a = \frac{\cos^2 \theta}{\sigma_x^2} + \frac{\sin^2\theta}{\sigma_y^2}\f$
230 : *
231 : * \f$ b = \frac{\sin 2\theta}{2} \left(\frac{1}{\sigma_x^2} - \frac{1}{\sigma_y^2} \right)\f$
232 : *
233 : * \f$ c = \frac{\sin^2 \theta}{\sigma_x^2} + \frac{\cos^2\theta}{\sigma_y^2}\f$
234 : *
235 : * where \f$ \theta \f$ specifies a counter-clockwise rotation.
236 : *
237 : * This function calculates (\f$\sigma_x\f$,\f$\sigma_y\f$,\f$\theta\f$) from inputs (a,b,c). It adopts
238 : * the convention that the long axis of the ellipse is \f$\sigma_x\f$, and \f$\theta\f$ is chosen
239 : * appropriately. Note that \f$-\frac{\pi}{2} < \theta < \frac{\pi}{2}\f$.
240 : *
241 : * \tparam realT is the type to use for arithmetic
242 : *
243 : * \ingroup gen_math_gaussians
244 : */
245 : template <typename realT>
246 : void gaussian2D_gen2rot( realT &sigma_x, ///< [out] the width parameter in the rotated x-direction
247 : realT &sigma_y, ///< [out] the width parameter in the rotated y-direction
248 : realT &theta, ///< [out] the c.c.w rotation
249 : const realT a, ///< [in] the first rotation and scaling parameter
250 : const realT b, ///< [in] the second rotation and scaling parameter
251 : const realT c ///< [in] the the third rotation and scaling parameter
252 : )
253 : {
254 : realT x1, x2, s, s2, theta0, theta1;
255 :
256 : realT arg = a * a - 2 * a * c + 4 * b * b + c * c;
257 : if( arg < 0 )
258 : {
259 : x2 = 0.5 * ( a + c );
260 : x1 = x2;
261 : }
262 : else
263 : {
264 : // There are two roots, one is 1/sigma_x^2 and one is 1/sigma_y^2
265 : x2 = 0.5 * ( a + c ) + 0.5 * sqrt( arg );
266 : x1 = 0.5 * ( a + c ) - 0.5 * sqrt( arg );
267 : }
268 :
269 : // Our convention is that the larger width is sigma_x
270 : sigma_x = sqrt( 1. / x1 );
271 : sigma_y = sqrt( 1. / x2 );
272 :
273 : //----------------------------------------------------------------------------//
274 : // Choosing theta:
275 : // There is an ambiguity in theta and which direction (x,y) is the long axis
276 : // Here we choose theta so that it specifies the long direction, sigma_x
277 : //----------------------------------------------------------------------------//
278 :
279 : // s is (sin(theta))^2
280 : // s2 = (a-x2)*(a-x2)/(b*b + (a-x2)*(a-x2));
281 : s = ( a - x1 ) * ( a - x1 ) / ( b * b + ( a - x1 ) * ( a - x1 ) );
282 :
283 : // First check if x1-x2 will be close to zero
284 :
285 : if( fabs( x1 - x2 ) < 1e-12 )
286 : {
287 : theta = 0;
288 : return;
289 : }
290 :
291 : theta0 = 0.5 * asin( 2 * b / ( x1 - x2 ) ); // This always gives the correct sign
292 : theta1 = asin( sqrt( s ) ); // This always gives the correct magnitude, and seems to be more accurate
293 :
294 : // Compare signs. If they match, then we use theta1
295 : if( std::signbit( theta0 ) == std::signbit( theta1 ) )
296 : {
297 : theta = theta1;
298 : }
299 : else
300 : {
301 : // Otherwise, we switch quadrants
302 : theta = asin( -sqrt( s ) );
303 : }
304 : }
305 :
306 : /// Convert from (\f$\sigma_x\f$,\f$\sigma_y\f$,\f$\theta\f$) to (a,b,c) for the elliptical Gaussian.
307 : /** The general 2D elliptical Gaussian
308 : *
309 : * \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$
310 : *
311 : * can be expressed in terms of widths and a rotation angle using:
312 : *
313 : * \f$ a = \frac{\cos^2 \theta}{\sigma_x^2} + \frac{\sin^2\theta}{\sigma_y^2}\f$
314 : *
315 : * \f$ b = \frac{\sin 2\theta}{2} \left(\frac{1}{\sigma_x^2} - \frac{1}{\sigma_y^2} \right)\f$
316 : *
317 : * \f$ c = \frac{\sin^2 \theta}{\sigma_x^2} + \frac{\cos^2\theta}{\sigma_y^2}\f$
318 : *
319 : * where \f$ \theta \f$ specifies a counter-clockwise rotation.
320 : *
321 : * This function calculates (a,b,c) from inputs (\f$\sigma_x\f$,\f$\sigma_y\f$,\f$\theta\f$).
322 : *
323 : * \tparam realT is the type to use for arithmetic
324 : *
325 : * \ingroup gen_math_gaussians
326 : */
327 : template <typename realT>
328 : void gaussian2D_rot2gen( realT &a, ///< [out] the first rotation and scaling parameter
329 : realT &b, ///< [out] the second rotation and scaling parameter
330 : realT &c, ///< [out] the the third rotation and scaling parameter
331 : const realT sigma_x, ///< [in] the width parameter in the rotated x-direction
332 : const realT sigma_y, ///< [in] the width parameter in the rotated y-direction
333 : const realT theta ///< [in] the c.c.w rotation
334 : )
335 : {
336 : realT sn, cs, sx2, sy2;
337 :
338 : sn = sin( theta );
339 : cs = cos( theta );
340 : sx2 = sigma_x * sigma_x;
341 : sy2 = sigma_y * sigma_y;
342 :
343 : a = cs * cs / sx2 + sn * sn / sy2;
344 : b = sn * cs * ( 1. / sx2 - 1. / sy2 );
345 : c = sn * sn / sx2 + cs * cs / sy2;
346 : }
347 :
348 : /// Find value at position (x,y) of the 2D rotated elliptical Gaussian
349 : /**
350 : * Computes:
351 : * \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$
352 : *
353 : * where, for a counter-clockwise rotation \f$ \theta \f$, we have
354 : *
355 : * \f$ a = \frac{\cos^2 \theta}{\sigma_x^2} + \frac{\sin^2\theta}{\sigma_y^2}\f$
356 : *
357 : * \f$ b = \frac{\sin 2\theta}{2} \left(\frac{1}{\sigma_x^2} - \frac{1}{\sigma_y^2} \right)\f$
358 : *
359 : * \f$ c = \frac{\sin^2 \theta}{\sigma_x^2} + \frac{\cos^2\theta}{\sigma_y^2}\f$
360 : *
361 : * This is a convenience wrapper for the general elliptical Gaussian function gaussian2D(), where here
362 : * (a,b,c) are first calculated from (\f$\sigma_x\f$,\f$\sigma_y\f$,\f$\theta\f$). This will in general be
363 : * slower due to the trig function calls, so the (a,b,c) version should be used most of the time.
364 : *
365 : * The functions gaussian2D_gen2rot() and gaussian2D_rot2gen() provide conversions from
366 : * (a,b,c) to (\f$\sigma_x\f$,\f$\sigma_y\f$,\f$\theta\f$) and back. The function gaussian2D_rot() is a
367 : * wrapper for this function, which instead accepts (\f$\sigma_x\f$,\f$\sigma_y\f$,\f$\theta\f$) as inputs.
368 : *
369 : * \returns the value of the 2D elliptical Gaussian at (x,y)
370 : *
371 : * \tparam realT is the type to use for arithmetic
372 : *
373 : * \ingroup gen_math_gaussians
374 : */
375 : template <typename realT>
376 : realT gaussian2D_ang( const realT x, ///< [in] the x-position at which to evaluate the Gaussian
377 : const realT y, ///< [in] the y-positoin at which to evaluate the Gaussian
378 : const realT G0, ///< [in] the constant to add to the Gaussian
379 : const realT G, ///< [in] the scaling factor (peak height = G-G0)
380 : const realT x0, ///< [in] the x-coordinate of the center
381 : const realT y0, ///< [in] the y-coordinate of the center
382 : const realT sigma_x, ///< [in] the width in the rotated x direction
383 : const realT sigma_y, ///< [in] the width in the rotated y direction
384 : const realT theta ///< [in] the counter-clockwise rotation angle
385 : )
386 : {
387 : realT a, b, c;
388 :
389 : gaussian2D_rot2gen( a, b, c, sigma_x, sigma_y, theta );
390 :
391 : return gaussian2D( x, y, G0, G, x0, y0, a, b, c );
392 : }
393 :
394 : /// Fill in an array with the 2D general elliptical Gaussian
395 : /**
396 : * At each pixel (x,y) of the array this computes:
397 : *
398 : * \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$
399 : *
400 : * where, for a counter-clockwise rotation \f$ \theta \f$, we have
401 : *
402 : * \f$ a = \frac{\cos^2 \theta}{\sigma_x^2} + \frac{\sin^2\theta}{\sigma_y^2}\f$
403 : *
404 : * \f$ b = \frac{\sin 2\theta}{2} \left(\frac{1}{\sigma_x^2} - \frac{1}{\sigma_y^2} \right)\f$
405 : *
406 : * \f$ c = \frac{\sin^2 \theta}{\sigma_x^2} + \frac{\cos^2\theta}{\sigma_y^2}\f$
407 : *
408 : * In this version the parameters are specified directly as (a,b,c), in particular avoiding the trig function calls.
409 : * This should be much more efficient, and so this version should be used inside fitting routines, etc.
410 : *
411 : * The functions gaussian2D_gen2rot() and gaussian2D_rot2gen() provide conversions from
412 : * (a,b,c) to (\f$\sigma_x\f$,\f$\sigma_y\f$,\f$\theta\f$) and back. The function gaussian2D_ang() is a
413 : * wrapper for this function, which instead accepts (\f$\sigma_x\f$,\f$\sigma_y\f$,\f$\theta\f$) as inputs.
414 : *
415 : * \tparam realT is the type to use for arithmetic
416 : *
417 : * \ingroup gen_math_gaussians
418 : */
419 : template <typename realT>
420 : void gaussian2D( realT *arr, ///< [out] the native array to populate with the Gaussian function
421 : size_t nx, ///< [in] the x-size of the array, in pixels
422 : size_t ny, ///< [in] the y-size of the array, in pixels
423 : const realT G0, ///< [in] the constant to add to the Gaussian
424 : const realT G, ///< [in] the scaling factor (peak = G)
425 : const realT x0, ///< [in] the x-coordinate of the center, in pixels
426 : const realT y0, ///< [in] the y-coordinate of the center, in pixels
427 : const realT a, ///< [in] is the first rotation and scaling factor
428 : const realT b, ///< [in] is the second rotation and scaling factor
429 : const realT c ///< [in] is the third rotation and scaling factor
430 : )
431 : {
432 : size_t idx;
433 :
434 : for( size_t j = 0; j < ny; ++j )
435 : {
436 : for( size_t i = 0; i < nx; ++i )
437 : {
438 : idx = i + j * nx;
439 :
440 : arr[idx] = gaussian2D( (realT)i, (realT)j, G0, G, x0, y0, a, b, c );
441 : }
442 : }
443 : }
444 :
445 : /// Fill in an array with the 2D general elliptical Gaussian
446 : /**
447 : * At each pixel (x,y) of the array this computes:
448 : *
449 : * \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$
450 : *
451 : * where, for a counter-clockwise rotation \f$ \theta \f$, we have
452 : *
453 : * \f$ a = \frac{\cos^2 \theta}{\sigma_x^2} + \frac{\sin^2\theta}{\sigma_y^2}\f$
454 : *
455 : * \f$ b = \frac{\sin 2\theta}{2} \left(\frac{1}{\sigma_x^2} - \frac{1}{\sigma_y^2} \right)\f$
456 : *
457 : * \f$ c = \frac{\sin^2 \theta}{\sigma_x^2} + \frac{\cos^2\theta}{\sigma_y^2}\f$
458 : *
459 : * This is a convenience wrapper for the general elliptical Gaussian function
460 : * gaussian2D( realT *, size_t, size_t, const realT, const realT, const realT, const realT, const realT, const realT,
461 : * const realT) , where here (a,b,c) are first calculated from (\f$\sigma_x\f$,\f$\sigma_y\f$,\f$\theta\f$).
462 : *
463 : * The functions gaussian2D_gen2rot() and gaussian2D_rot2gen() provide conversions from
464 : * (a,b,c) to (\f$\sigma_x\f$,\f$\sigma_y\f$,\f$\theta\f$) and back.
465 : *
466 : *
467 : * \tparam realT is the type to use for arithmetic
468 : *
469 : * \ingroup gen_math_gaussians
470 : */
471 : template <typename realT>
472 : void gaussian2D_ang( realT *arr, ///< [out] the native array to populate with the Gaussian function
473 : size_t nx, ///< [in] the x-size of the array, in pixels
474 : size_t ny, ///< [in] the y-size of the array, in pixels
475 : const realT G0, ///< [in] the constant to add to the Gaussian
476 : const realT G, ///< [in] the scaling factor (peak = G)
477 : const realT x0, ///< [in] the x-coordinate of the center, in pixels
478 : const realT y0, ///< [in] the y-coordinate of the center, in pixels
479 : const realT sigma_x, ///< [in] the width in the rotated x direction, in pixels
480 : const realT sigma_y, ///< [in] the width in the rotated y direction, in pixels
481 : const realT theta ///< [in] the counter-clockwise rotation angle, in radians
482 : )
483 : {
484 : realT a, b, c;
485 : size_t idx;
486 :
487 : // Get the (a,b,c) parameters so the trig only happens once
488 : gaussian2D_rot2gen( a, b, c, sigma_x, sigma_y, theta );
489 :
490 : gaussian2D( arr, nx, ny, G0, G, x0, y0, a, b, c );
491 : }
492 :
493 : /// Calculate the Jacobian at position (x,y) for the 2D general elliptical Gaussian
494 : /** \note this has not been verified and may be incorrect.
495 : *
496 : * Given:
497 : *
498 : * \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$
499 : *
500 : * we calculate
501 : *
502 : * \f$ \frac{\partial G}{\partial G_0} = 1 \f$
503 : *
504 : * \f$ \frac{\partial G}{\partial G} = (G - G_0)/A \f$
505 : *
506 : * \f$ \frac{\partial G}{\partial x_0} = 0.5(G-G_0)( 2a(x-x_0) + b(y-y_0)) \f$
507 : *
508 : * \f$ \frac{\partial G}{\partial y_0} = 0.5(G-G_0) ( b(x-x_0) + 2c(y-y_0)) \f$
509 : *
510 : * \f$ \frac{\partial G}{\partial a} = -0.5(G-G_0) ( (x-x_0)^2 ) \f$
511 : *
512 : * \f$ \frac{\partial G}{\partial b} = -0.5(G-G_0) ( (x-x_0)(y-y_0) ) \f$
513 : *
514 : * \f$ \frac{\partial G}{\partial c} = -0.5(G-G_0) ( (y-y_0)^2 )\f$
515 : *
516 : * \param j is a 7 element vector which is populated with the derivatives
517 : * \param x is the x-position at which to evaluate the Gaussian
518 : * \param y is the y-positoin at which to evaluate the Gaussian
519 : * \param G0 is the constant to add to the Gaussian
520 : * \param G is the scaling factor (peak = A)
521 : * \param x0 is the x-coordinate of the center
522 : * \param y0 is the y-coordinate of the center
523 : * \param a is the first rotation and scaling factor
524 : * \param b is the second rotation and scaling factor
525 : * \param c is the third rotation and scaling factor
526 : *
527 : *
528 : * \tparam realT is the type to use for arithmetic
529 : *
530 : * \ingroup gen_math_gaussians
531 : */
532 : template <typename realT>
533 : void gaussian2D_jacobian( realT *j,
534 : const realT x,
535 : const realT y,
536 : const realT G0,
537 : const realT G,
538 : const realT x0,
539 : const realT y0,
540 : const realT a,
541 : const realT b,
542 : const realT c )
543 : {
544 : realT G_G0 = gaussian2D<realT>( x, y, 0, G, x0, y0, a, b, c );
545 :
546 : j[0] = (realT)1;
547 :
548 : j[1] = ( G_G0 ) / G;
549 :
550 : j[2] = 0.5 * G_G0 * ( 2 * a * ( x - x0 ) + b * ( y - y0 ) );
551 :
552 : j[3] = 0.5 * G_G0 * ( b * ( x - x0 ) + 2 * c * ( y - y0 ) );
553 :
554 : j[4] = -0.5 * G_G0 * ( ( x - x0 ) * ( x - x0 ) );
555 :
556 : j[5] = -0.5 * G_G0 * ( ( x - x0 ) * ( y - y0 ) );
557 :
558 : j[6] = -0.5 * G_G0 * ( ( y - y0 ) * ( y - y0 ) );
559 : }
560 :
561 : } // namespace func
562 : } // namespace math
563 : } // namespace mx
564 :
565 : #endif // gaussian_hpp
|
