Line data Source code
1 : /** \file fitGaussian.hpp
2 : * \author Jared R. Males
3 : * \brief Tools for fitting Gaussians to data.
4 : * \ingroup fitting_files
5 : *
6 : */
7 :
8 : //***********************************************************************//
9 : // Copyright 2015, 2016, 2017 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 fitGaussian_hpp
28 : #define fitGaussian_hpp
29 :
30 : #include "array2FitGaussian1D.hpp"
31 : #include "array2FitGaussian2D.hpp"
32 :
33 : #include "../../mxlib.hpp"
34 : #include "levmarInterface.hpp"
35 : #include "../func/gaussian.hpp"
36 : #include "../constants.hpp"
37 :
38 : #include "../../improc/eigenImage.hpp"
39 :
40 : namespace mx
41 : {
42 : namespace math
43 : {
44 : namespace fit
45 : {
46 :
47 : /** \defgroup gaussian_peak_fit Gaussians
48 : * \brief Fitting Gaussians to data.
49 : *
50 : * The Gaussian function is fit to data.
51 : *
52 : * \ingroup peak_fit
53 : */
54 :
55 : // forward
56 : template <typename _realT>
57 : struct gaussian1D_fitter;
58 :
59 : ///\ref levmarInterface fitter structure for the symmetric Gaussian.
60 : /** \ingroup gaussian_peak_fit
61 : *
62 : */
63 : template <typename _realT>
64 : struct gaussian1D_fitter
65 : {
66 : typedef _realT realT;
67 :
68 : static const int nparams = 4;
69 :
70 0 : static void func( realT *p, realT *hx, int m, int n, void *adata )
71 : {
72 0 : array2FitGaussian1D<realT> *arr = (array2FitGaussian1D<realT> *)adata;
73 :
74 0 : realT G0 = arr->G0( p );
75 0 : realT G = arr->G( p );
76 0 : realT x0 = arr->x0( p );
77 0 : realT sigma = arr->sigma( p );
78 :
79 0 : if( arr->m_mask )
80 : {
81 0 : if( arr->m_coords )
82 : {
83 0 : for( int i = 0; i < arr->m_nx; ++i )
84 : {
85 0 : if( arr->m_mask[i] == 0 )
86 0 : continue;
87 0 : hx[i] = func::gaussian<realT>( arr->m_coords[i], G0, G, x0, sigma ) - arr->m_data[i];
88 : }
89 : }
90 : else
91 : {
92 0 : for( int i = 0; i < arr->m_nx; ++i )
93 : {
94 0 : if( arr->m_mask[i] == 0 )
95 0 : continue;
96 0 : hx[i] = func::gaussian<realT>( i, G0, G, x0, sigma ) - arr->m_data[i];
97 : }
98 : }
99 : }
100 : else
101 : {
102 0 : if( arr->m_coords )
103 : {
104 0 : for( int i = 0; i < arr->m_nx; ++i )
105 : {
106 0 : hx[i] = func::gaussian<realT>( arr->m_coords[i], G0, G, x0, sigma ) - arr->m_data[i];
107 : }
108 : }
109 : else
110 : {
111 0 : for( int i = 0; i < arr->m_nx; ++i )
112 : {
113 0 : hx[i] = func::gaussian<realT>( i, G0, G, x0, sigma ) - arr->m_data[i];
114 : }
115 : }
116 : }
117 0 : }
118 : };
119 :
120 : extern template struct gaussian1D_fitter<float>;
121 : extern template struct gaussian1D_fitter<double>;
122 :
123 : /// Class to manage fitting a 1D Gaussian to data via the \ref levmarInterface
124 : /** Fits the following function to the data:
125 : * \f$ G(x) = G_0 + G\exp[-(0.5/\sigma^2)((x-x_0)^2)]\f$
126 : *
127 : * Can use a vector of x-coordinates, or use the array index as the position corresponding to each y-value.
128 : * A mask can be provided, where any 0 entries cause that position to be ignored.
129 : *
130 : * Any of the parameters can be fixed, and not included in the fit.
131 : *
132 : * \code
133 : *
134 : * \endcode
135 : *
136 : *
137 : * \tparam fitterT a type meeting the above requirements.
138 : *
139 : * \ingroup gaussian_peak_fit
140 : *
141 : */
142 : // template<typename fitterT>
143 : template <typename _realT>
144 : class fitGaussian1D : public levmarInterface<gaussian1D_fitter<_realT>> // fitterT>
145 : {
146 :
147 : public:
148 : typedef _realT realT;
149 : typedef gaussian1D_fitter<realT> fitterT;
150 :
151 : protected:
152 : array2FitGaussian1D<realT> arr;
153 :
154 0 : void initialize()
155 : {
156 0 : this->allocate_params( arr.nparams() );
157 0 : this->adata = &arr;
158 0 : }
159 :
160 : public:
161 0 : fitGaussian1D()
162 0 : {
163 0 : this->initialize();
164 0 : }
165 :
166 0 : ~fitGaussian1D()
167 : {
168 0 : }
169 :
170 : /// Set whether each parameter is fixed.
171 : /** Sets the parameter indices appropriately.
172 : */
173 0 : void setFixed( bool G0, ///< [in] if true, then G0 will be not be part of the fit
174 : bool G, ///< [in] if true, then G will be not be part of the fit
175 : bool x0, ///< [in] if true, then x0 will be not be part of the fit
176 : bool sigma ///< [in] if true, then sigma will be not be part of the fit
177 : )
178 : {
179 0 : arr.setFixed( G0, G, x0, sigma );
180 0 : this->allocate_params( arr.nparams() );
181 0 : }
182 :
183 : /// Set the initial guess for a symmetric Gaussian.
184 : /** Also works for the general case, setting the same width in both directions.
185 : */
186 0 : void setGuess( realT G0, ///< [in] the constant background level
187 : realT G, ///< [in] the peak scaling
188 : realT x0, ///< [in] the center x-coordinate
189 : realT sigma ///< [in] the width parameter
190 : )
191 : {
192 0 : arr.G0( this->p, G0 );
193 0 : arr.G( this->p, G );
194 0 : arr.x0( this->p, x0 );
195 0 : arr.sigma( this->p, sigma );
196 0 : }
197 :
198 : /// Set the data aray.
199 0 : void setArray( realT *data, int nx )
200 : {
201 0 : arr.m_data = data;
202 0 : arr.m_nx = nx;
203 :
204 0 : this->n = nx;
205 0 : }
206 :
207 : /// Set the data aray.
208 0 : void setArray( realT *data, realT *coords, int nx )
209 : {
210 0 : arr.m_data = data;
211 0 : arr.m_coords = coords;
212 0 : arr.m_nx = nx;
213 :
214 0 : this->n = nx;
215 0 : }
216 :
217 : /// Do the fit.
218 0 : int fit()
219 : {
220 : fitterT fitter;
221 :
222 0 : return levmarInterface<fitterT>::fit();
223 : }
224 :
225 : /// Get the current value of G0, the constant.
226 : /**
227 : * \returns the current value of G0, which is p[0].
228 : */
229 0 : realT G0()
230 : {
231 0 : return arr.G0( this->p );
232 : }
233 :
234 : /// Get the peak scaling.
235 0 : realT G()
236 : {
237 0 : return arr.G( this->p );
238 : }
239 :
240 : /// Get the center x-coordinate
241 0 : realT x0()
242 : {
243 0 : return arr.x0( this->p );
244 : }
245 :
246 : /// Return sigma
247 : /**
248 : */
249 0 : realT sigma()
250 : {
251 0 : return arr.sigma( this->p );
252 : }
253 : };
254 :
255 : extern template class fitGaussian1D<float>;
256 : extern template class fitGaussian1D<double>;
257 :
258 : /// Alias for the fitGaussian1D type fitting the gaussian.
259 : /** \ingroup gaussian_peak_fit
260 : */
261 : // template<typename realT>
262 : // using fitGaussian1D = mx::math::fit::fitGaussian1D<mx::math::fit::gaussian1D_fitter<realT>>;
263 :
264 : // forward
265 : template <typename _realT>
266 : struct gaussian2D_sym_fitter;
267 :
268 : // forward
269 : template <typename _realT>
270 : struct gaussian2D_gen_fitter;
271 :
272 : template <typename _realT>
273 : struct gaussian2D_gen_fitter_bgfixed;
274 :
275 : /// Class to manage fitting a 2D Gaussian to data via the \ref levmarInterface
276 : /** Can fit either the symmetric Gaussian, or the rotated asymmetric Gaussian. Individual parameters can be fixed,
277 : * except in the asymmetric case \f$\sigma_x, \sigma_y, \theta\f$ must currently be fixed together.
278 : *
279 : * In addition to the requirements on fitterT specified by \ref levmarInterface
280 : * this class also requires the following in fitterT
281 : * \code
282 : * static const int nparams = 7; //where the number 7 is replaced by the number of parameters that fitterT expects to
283 : * fit.
284 : *
285 : * void paramNormalizer( array2FitGaussian2D * arr,
286 : * realT *,
287 : * int); //A function which converts from input parameters to fitting parameters. May do nothing.
288 : *
289 : * \endcode
290 : *
291 : *
292 : * \tparam fitterT a type meeting the above requirements.
293 : *
294 : * \ingroup gaussian_peak_fit
295 : *
296 : * "[test doc]"
297 : *
298 : */
299 : template <typename fitterT>
300 : class fitGaussian2D : public levmarInterface<fitterT>
301 : {
302 :
303 : public:
304 : typedef typename fitterT::realT realT;
305 :
306 : array2FitGaussian2D<realT>
307 : arr; /**< Data array to pass to the levmar library.
308 : Contains the actual data plus the parameters if fixed.*/
309 :
310 21 : void initialize()
311 : {
312 : if( fitterT::maxNparams == 5 )
313 : {
314 0 : arr.setSymmetric();
315 : }
316 : else
317 : {
318 21 : arr.setGeneral();
319 : }
320 21 : this->allocate_params( arr.nparams() );
321 21 : this->adata = &arr;
322 21 : }
323 :
324 21 : fitGaussian2D()
325 21 : {
326 21 : initialize();
327 21 : }
328 :
329 21 : ~fitGaussian2D()
330 : {
331 21 : }
332 :
333 : /// Set whether each parameter is fixed.
334 : /** Sets the parameter indices appropriately.
335 : */
336 0 : void setFixed( bool G0, ///< [in] if true, then G0 will be not be part of the fit
337 : bool G, ///< [in] if true, then G will be not be part of the fit
338 : bool x0, ///< [in] if true, then x0 will be not be part of the fit
339 : bool y0, ///< [in] if true, then y0 will be not be part of the fit
340 : bool sigma_x, ///< [in] if true, then sigma_x will be not be part of the fit
341 : bool sigma_y, ///< [in] if true, then sigma_y will be not be part of the fit
342 : bool theta ///< [in] if true, then theta will be not be part of the fit
343 : )
344 : {
345 0 : arr.setFixed( G0, G, x0, y0, sigma_x, sigma_y, theta );
346 0 : this->allocate_params( arr.nparams() );
347 0 : }
348 :
349 : /// Set whether each parameter is fixed.
350 : /** Sets the parameter indices appropriately.
351 : */
352 0 : void setFixed( bool G0, ///< [in] if true, then G0 will be not be part of the fit
353 : bool G, ///< [in] if true, then G will be not be part of the fit
354 : bool x0, ///< [in] if true, then x0 will be not be part of the fit
355 : bool y0, ///< [in] if true, then y0 will be not be part of the fit
356 : bool sigma ///< [in] if true, then sigma will be not be part of the fit
357 : )
358 : {
359 0 : arr.setFixed( G0, G, x0, y0, sigma, sigma, sigma );
360 0 : this->allocate_params( arr.nparams() );
361 0 : }
362 :
363 : /// Set the initial guess for a symmetric Gaussian.
364 : /** Also works for the general case, setting the same width in both directions.
365 : */
366 0 : void setGuess( realT G0, ///< [in] the constant background level
367 : realT G, ///< [in] the peak scaling
368 : realT x0, ///< [in] the center x-coordinate
369 : realT y0, ///< [in] the center y-coordinate
370 : realT sigma ///< [in] the width parameter
371 : )
372 : {
373 0 : arr.G0( this->p, G0 );
374 0 : arr.G( this->p, G );
375 0 : arr.x0( this->p, x0 );
376 0 : arr.y0( this->p, y0 );
377 0 : arr.sigma( this->p, sigma );
378 0 : }
379 :
380 : /// Set the initial guess for the general Gaussian.
381 8 : void setGuess( realT G0, ///< [in] the constant background level
382 : realT G, ///< [in] the peak scaling
383 : realT x0, ///< [in] the center x-coordinate
384 : realT y0, ///< [in] the center y-coordinate
385 : realT sigma_x, ///< [in] the width parameter in the rotated x-direction (the long axis)
386 : realT sigma_y, ///< [in] the width parameter in the rotated y-direction
387 : realT theta ///< [in] the angle of the long axis (always sigma_x)
388 : )
389 : {
390 :
391 8 : arr.G0( this->p, G0 );
392 8 : arr.G( this->p, G );
393 8 : arr.x0( this->p, x0 );
394 8 : arr.y0( this->p, y0 );
395 8 : arr.sigma_x( this->p, sigma_x );
396 8 : arr.sigma_y( this->p, sigma_y );
397 8 : arr.theta( this->p, theta );
398 8 : }
399 :
400 : /// Set the data aray.
401 8 : void setArray( realT *data, ///< [in] The 2D array of data to fit.
402 : int nx, ///< [in] the x size of the data
403 : int ny ///< [in] the y size of the data
404 : )
405 : {
406 8 : arr.data = data;
407 8 : arr.nx = nx;
408 8 : arr.ny = ny;
409 8 : arr.mask = nullptr;
410 8 : this->n = nx * ny;
411 8 : }
412 :
413 : /// Set the data aray, with a mask.
414 0 : void setArray(
415 : realT *data, ///< [in] The 2D array of data to fit.
416 : int nx, ///< [in] the x size of the data
417 : int ny, ///< [in] the y size of the data
418 : realT *mask ///< [in] Array of same size as data. Any 0 pixels in this array will be excluded from the fit.
419 : )
420 : {
421 0 : arr.data = data;
422 0 : arr.nx = nx;
423 0 : arr.ny = ny;
424 0 : arr.mask = mask;
425 :
426 0 : this->n = 0;
427 : int idx_mat;
428 0 : for( int j = 0; j < nx; ++j )
429 : {
430 0 : for( int i = 0; i < ny; ++i )
431 : {
432 0 : idx_mat = i + j * nx;
433 :
434 0 : this->n += arr.mask[idx_mat];
435 : }
436 : }
437 0 : }
438 :
439 : /// Do the fit.
440 8 : int fit()
441 : {
442 : fitterT fitter;
443 8 : fitter.paramNormalizer( &arr, this->p, 1 );
444 :
445 8 : levmarInterface<fitterT>::fit();
446 :
447 8 : fitter.paramNormalizer( &arr, this->p, -1 );
448 8 : fitter.paramNormalizer(
449 8 : &arr, this->init_p, -1 ); // The normalized version is stored, so fix it before possible output.
450 :
451 8 : return 0;
452 : }
453 :
454 : /// Get the current value of G0, the constant.
455 : /**
456 : * \returns the current value of G0.
457 : */
458 0 : realT G0()
459 : {
460 0 : return arr.G0( this->p );
461 : }
462 :
463 : /// Get the peak scaling.
464 8 : realT G()
465 : {
466 8 : return arr.G( this->p );
467 : }
468 :
469 : /// Get the center x-coordinate
470 8 : realT x0()
471 : {
472 8 : return arr.x0( this->p );
473 : }
474 :
475 : /// Get the center y-coordinate
476 8 : realT y0()
477 : {
478 8 : return arr.y0( this->p );
479 : }
480 :
481 : /// Return the width parameter
482 : /** As described for the symmetric Gaussian.
483 : *
484 : * For the general Gaussian, this returns \f$ \sigma = \sqrt{ \sigma_x^2 + \sigma_y^2} \f$.
485 : */
486 0 : realT sigma()
487 : {
488 0 : return arr.sigma( this->p );
489 : }
490 :
491 : /// Return the full-width at half maximum
492 : /** This is a simple scaling of the sigma() result
493 : */
494 0 : realT fwhm()
495 : {
496 0 : return func::sigma2fwhm( arr.sigma( this->p ) );
497 : }
498 :
499 : /// Return the width parameter on the long axis.
500 0 : realT sigma_x()
501 : {
502 0 : return arr.sigma_x( this->p );
503 : }
504 :
505 : /// Return the width parameter on the short axis.
506 0 : realT sigma_y()
507 : {
508 0 : return arr.sigma_y( this->p );
509 : }
510 :
511 : /// Return the orientation of the long axis.
512 0 : realT theta()
513 : {
514 0 : return arr.theta( this->p );
515 : }
516 : };
517 :
518 : ///\ref levmarInterface fitter structure for the symmetric Gaussian.
519 : /** \ingroup gaussian_peak_fit
520 : *
521 : * "[test doc]"
522 : */
523 : template <typename _realT>
524 : struct gaussian2D_sym_fitter
525 : {
526 : typedef _realT realT;
527 :
528 : static const int maxNparams = 5;
529 :
530 0 : static void func( realT *p, realT *hx, int m, int n, void *adata )
531 : {
532 0 : array2Fit<realT> *arr = (array2Fit<realT> *)adata;
533 :
534 : size_t idx_mat, idx_dat;
535 :
536 0 : idx_dat = 0;
537 :
538 0 : for( int j = 0; j < arr->ny; j++ )
539 : {
540 0 : for( int i = 0; i < arr->nx; i++ )
541 : {
542 0 : idx_mat = i + j * arr->nx;
543 :
544 0 : hx[idx_dat] = func::gaussian2D<realT>( i, j, p[0], p[1], p[2], p[3], p[4] ) - arr->data[idx_mat];
545 :
546 0 : idx_dat++;
547 : }
548 : }
549 0 : }
550 :
551 : /// Does nothing in this case.
552 0 : void paramNormalizer( array2FitGaussian2D<realT> *arr, realT *p, int dir )
553 : {
554 0 : return;
555 : }
556 : };
557 :
558 : ///\ref levmarInterface fitter structure for the general elliptical Gaussian.
559 : /** \ingroup gaussian_peak_fit
560 : *
561 : */
562 : template <typename _realT>
563 : struct gaussian2D_gen_fitter
564 : {
565 : typedef _realT realT;
566 :
567 : static const int maxNparams = 7;
568 :
569 301 : static void func( realT *p, realT *hx, int m, int n, void *adata )
570 : {
571 301 : array2FitGaussian2D<realT> *arr = (array2FitGaussian2D<realT> *)adata;
572 :
573 : size_t idx_mat, idx_dat;
574 :
575 301 : realT G0 = arr->G0( p ); // 0
576 301 : realT G = arr->G( p ); // 1
577 301 : realT x0 = arr->x0( p ); // 2
578 301 : realT y0 = arr->y0( p ); // 3
579 301 : realT a = arr->a( p ); // 4
580 301 : realT b = arr->b( p ); // 5
581 301 : realT c = arr->c( p ); // 6
582 :
583 : // Check for positive-definiteness of {{a b}{b c}}
584 301 : if( a * c - b * b <= 0 || a <= 0 || c <= 0 || a + c <= 2 * fabs( b ) )
585 : {
586 0 : idx_dat = 0;
587 : // If it's not positive-definite, then we just fill in with the value of the image itself.
588 0 : if( arr->mask == nullptr )
589 : {
590 0 : if( arr->weights == nullptr)
591 : {
592 0 : for( int j = 0; j < arr->ny; ++j )
593 : {
594 0 : for( int i = 0; i < arr->ny; ++i )
595 : {
596 0 : idx_mat = i + j * arr->nx;
597 0 : hx[idx_dat] = arr->data[idx_mat];
598 0 : ++idx_dat;
599 : }
600 : }
601 : }
602 : else
603 : {
604 0 : for( int j = 0; j < arr->ny; ++j )
605 : {
606 0 : for( int i = 0; i < arr->ny; ++i )
607 : {
608 0 : idx_mat = i + j * arr->nx;
609 0 : hx[idx_dat] = arr->data[idx_mat] * arr->weights[idx_mat];
610 0 : ++idx_dat;
611 : }
612 : }
613 : }
614 : }
615 : else
616 : {
617 0 : for( int j = 0; j < arr->ny; ++j )
618 : {
619 0 : for( int i = 0; i < arr->ny; ++i )
620 : {
621 0 : idx_mat = i + j * arr->nx;
622 0 : if( arr->mask[idx_mat] == 0 )
623 : {
624 0 : continue;
625 : }
626 0 : hx[idx_dat] = arr->data[idx_mat];
627 0 : ++idx_dat;
628 : }
629 : }
630 : }
631 0 : return;
632 : }
633 :
634 : // If positive-definite, now actually calculate
635 301 : idx_dat = 0;
636 :
637 301 : if( arr->mask == nullptr )
638 : {
639 301 : if(arr->weights == nullptr)
640 : {
641 9764 : for( int j = 0; j < arr->ny; ++j )
642 : {
643 307700 : for( int i = 0; i < arr->nx; ++i )
644 : {
645 298237 : idx_mat = i + j * arr->nx;
646 :
647 298237 : hx[idx_dat] = func::gaussian2D<realT>( i, j, G0, G, x0, y0, a, b, c ) - arr->data[idx_mat];
648 :
649 298237 : ++idx_dat;
650 : }
651 : }
652 : }
653 : else
654 : {
655 0 : for( int j = 0; j < arr->ny; ++j )
656 : {
657 0 : for( int i = 0; i < arr->nx; ++i )
658 : {
659 0 : idx_mat = i + j * arr->nx;
660 :
661 0 : hx[idx_dat] = (func::gaussian2D<realT>( i, j, G0, G, x0, y0, a, b, c ) - arr->data[idx_mat])*arr->weights[idx_mat];
662 :
663 0 : ++idx_dat;
664 : }
665 : }
666 : }
667 : }
668 : else
669 : {
670 0 : for( int j = 0; j < arr->ny; ++j )
671 : {
672 0 : for( int i = 0; i < arr->nx; ++i )
673 : {
674 0 : idx_mat = i + j * arr->nx;
675 :
676 0 : if( arr->mask[idx_mat] == 0 )
677 : {
678 0 : continue;
679 : }
680 : else
681 : {
682 0 : hx[idx_dat] = func::gaussian2D<realT>( i, j, G0, G, x0, y0, a, b, c ) - arr->data[idx_mat];
683 :
684 0 : ++idx_dat;
685 : }
686 : }
687 : }
688 : }
689 : }
690 :
691 24 : void paramNormalizer( array2FitGaussian2D<realT> *arr, realT *p, int dir )
692 : {
693 : // Prepare for fit
694 24 : if( dir == 1 )
695 : {
696 : realT na, nb, nc;
697 8 : realT sx = arr->sigma_x( p );
698 8 : realT sy = arr->sigma_y( p );
699 8 : realT th = arr->theta( p );
700 :
701 8 : func::gaussian2D_rot2gen( na, nb, nc, sx, sy, th );
702 8 : arr->a( p, na );
703 8 : arr->b( p, nb );
704 8 : arr->c( p, nc );
705 8 : return;
706 : }
707 :
708 : // Convert after fit
709 16 : if( dir == -1 )
710 : {
711 : realT sx, sy, th;
712 16 : realT na = arr->a( p );
713 16 : realT nb = arr->b( p );
714 16 : realT nc = arr->c( p );
715 16 : func::gaussian2D_gen2rot( sx, sy, th, na, nb, nc );
716 :
717 16 : arr->sigma_x( p, sx );
718 16 : arr->sigma_y( p, sy );
719 16 : arr->theta( p, th );
720 16 : return;
721 : }
722 : }
723 : };
724 :
725 : extern template class fitGaussian2D<mx::math::fit::gaussian2D_sym_fitter<float>>;
726 : extern template class fitGaussian2D<mx::math::fit::gaussian2D_sym_fitter<double>>;
727 : extern template class fitGaussian2D<mx::math::fit::gaussian2D_gen_fitter<float>>;
728 : extern template class fitGaussian2D<mx::math::fit::gaussian2D_gen_fitter<double>>;
729 :
730 : /// Alias for the fitGaussian2D type fitting the symmetric gaussian.
731 : /** \ingroup gaussian_peak_fit
732 : */
733 : template <typename realT>
734 : using fitGaussian2Dsym = fitGaussian2D<mx::math::fit::gaussian2D_sym_fitter<realT>>;
735 :
736 : /// Alias for the fitGaussian2D type fitting the general elliptical gaussian.
737 : /** \ingroup gaussian_peak_fit
738 : */
739 : template <typename realT>
740 : using fitGaussian2Dgen = fitGaussian2D<mx::math::fit::gaussian2D_gen_fitter<realT>>;
741 :
742 : /// Form an estimate of the parameters of an elliptical Gaussian from a 2D image.
743 : /** Note that this assumes that there is no constant value (i.e. it is zero-ed).
744 : *
745 : * \ingroup gaussian_peak_fit
746 : */
747 : template <typename realT>
748 0 : int guessGauss2D_ang( realT &Ag, ///< [out] estimate of the peak
749 : realT &xg, ///< [out] estimate of the x-coordinate of the peak
750 : realT &yg, ///< [out] estimate of the y-coordinate of the peak
751 : realT &xFWHM, ///< [out] estimate of the x-FWHM
752 : realT &yFWHM, ///< [out] estimate of the y-FWHM
753 : realT &angG, ///< [out] estimate of the angle of the ellipse
754 : mx::improc::eigenImage<realT> &im, ///< [in] the image with an elliptical gaussian
755 : realT maxWidth, ///< [in] the width of the box to search for the maximum
756 : realT widthWidth, ///< [in] the radius of the circle to search for the widths
757 : realT nAngs, ///< [in] the number of angles at which to search for the widths
758 : realT xg0, ///< [in] an initial guess at the x-coordinate of the peak
759 : realT yg0 ///< [in] an initial guess at the y-coordinate of the peak
760 : )
761 : {
762 0 : xg = xg0;
763 0 : yg = yg0;
764 0 : Ag = im( (int)xg, (int)yg );
765 0 : for( int i = 0; i < 2 * maxWidth + 1; ++i )
766 : {
767 0 : for( int j = 0; j < 2 * maxWidth + 1; ++j )
768 : {
769 0 : if( im( (int)( xg0 - maxWidth + i ), (int)( yg0 - maxWidth + j ) ) > Ag )
770 : {
771 0 : Ag = im( (int)( xg0 - maxWidth + i ), (int)( yg0 - maxWidth + j ) );
772 0 : xg = xg0 - maxWidth + i;
773 0 : yg = yg0 - maxWidth + j;
774 : }
775 : }
776 : }
777 :
778 0 : realT dAng = math::two_pi<realT>() / nAngs;
779 : // std::vector<realT> dist(nAngs);
780 :
781 : realT c, s;
782 :
783 0 : realT maxD = 0;
784 0 : int maxDidx = 0;
785 0 : realT minD = widthWidth;
786 0 : int minDidx = 0;
787 :
788 0 : for( int i = 0; i < nAngs; ++i )
789 : {
790 0 : c = cos( i * dAng );
791 0 : s = sin( i * dAng );
792 :
793 0 : for( int j = 0; j < widthWidth; ++j )
794 : {
795 0 : if( im( (int)( xg + j * c ), (int)( yg + j * s ) ) <= 0.5 * Ag )
796 : {
797 : // dist[i] = j;
798 :
799 0 : if( j > maxD )
800 : {
801 0 : maxD = j;
802 0 : maxDidx = i;
803 : }
804 :
805 0 : if( j < minD )
806 : {
807 0 : minD = j;
808 0 : minDidx = i;
809 : }
810 0 : break;
811 : }
812 : }
813 : }
814 :
815 : // Take minang and move it by 90 degrees
816 0 : realT minang = fmod( minDidx * dAng - 0.5 * math::pi<realT>(), math::pi<realT>() );
817 0 : if( minang < 0 )
818 0 : minang = fmod( minang + math::two_pi<realT>(), math::pi<realT>() );
819 :
820 0 : realT maxang = fmod( maxDidx * dAng, math::pi<realT>() );
821 :
822 : // Now average
823 0 : angG = 0.5 * ( minang + maxang );
824 :
825 0 : xFWHM = 2 * maxD;
826 0 : yFWHM = 2 * minD;
827 :
828 0 : return 0; ///\returns 0 if successful
829 : }
830 :
831 : extern template int guessGauss2D_ang<float>( float &Ag,
832 : float &xg,
833 : float &yg,
834 : float &xFWHM,
835 : float &yFWHM,
836 : float &angG,
837 : mx::improc::eigenImage<float> &im,
838 : float maxWidth,
839 : float widthWidth,
840 : float nAngs,
841 : float xg0,
842 : float yg0 );
843 :
844 : extern template int guessGauss2D_ang<double>( double &Ag,
845 : double &xg,
846 : double &yg,
847 : double &xFWHM,
848 : double &yFWHM,
849 : double &angG,
850 : mx::improc::eigenImage<double> &im,
851 : double maxWidth,
852 : double widthWidth,
853 : double nAngs,
854 : double xg0,
855 : double yg0 );
856 :
857 : } // namespace fit
858 : } // namespace math
859 :
860 : } // namespace mx
861 :
862 : #endif // fitGaussian_hpp
|
