Line data Source code
1 : /** \file array2FitGaussian2D.hpp
2 : * \author Jared R. Males
3 : * \brief Wrapper for a native array to pass to \ref mx::math::fit::levmarInterface, with 2D Gaussian details.
4 : * \ingroup fitting_files
5 : *
6 : */
7 :
8 : //***********************************************************************//
9 : // Copyright 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 math_fit_array2FitGaussian2D_hpp
28 : #define math_fit_array2FitGaussian2D_hpp
29 :
30 : #include "../../mxlib.hpp"
31 :
32 : namespace mx
33 : {
34 : namespace math
35 : {
36 : namespace fit
37 : {
38 :
39 : /// Wrapper for a native array to pass to \ref levmarInterface, with 2D Gaussian details.
40 : /** Supports fixing G0, G, x0, and y0 independently. The shape and orientation can be fixed, but
41 : * for the general form, sigma_x, sigma_y, and theta can only be fixed together.
42 : * \ingroup gaussian_peak_fit
43 : */
44 : template <typename realT>
45 : struct array2FitGaussian2D
46 : {
47 : realT *data{ nullptr }; ///< Pointer to the array
48 : size_t nx{ 0 }; ///< X dimension of the array
49 : size_t ny{ 0 }; ///< Y dimension of the array
50 :
51 : realT *weights{ nullptr }; /**< Pointer to the (optional) weight array. Normally this should be the
52 : square root of the variance.*/
53 :
54 : realT *mask{ nullptr }; ///< Pointer to the (optional) mask array. Any 0 pixels are excluded from the fit.
55 :
56 : realT m_G0{ 0 };
57 : realT m_G{ 0 };
58 : realT m_x0{ 0 };
59 : realT m_y0{ 0 };
60 : realT m_sigma_x{ 0 };
61 : realT m_sigma_y{ 0 };
62 : realT m_theta{ 0 };
63 :
64 : realT m_a{ 0 };
65 : realT m_b{ 0 };
66 : realT m_c{ 0 };
67 :
68 : realT m_sigma{ 0 };
69 :
70 : int m_G0_idx{ 0 };
71 : int m_G_idx{ 1 };
72 : int m_x0_idx{ 2 };
73 : int m_y0_idx{ 3 };
74 : int m_sigma_x_idx{ 4 }; ///< Index of sigma_x in the parameters. Re-used for a.
75 : int m_sigma_y_idx{ 5 }; ///< Index of sigma_y in the parameters. Re-used for b.
76 : int m_theta_idx{ 6 }; ///< Index of theta in the parameters. Re-used for c.
77 :
78 : int m_sigma_idx{ 4 }; ///< Index of sigma in the symmetric case
79 :
80 : int m_nparams{ 7 };
81 : int m_maxNparams{ 7 };
82 :
83 : void setSymmetric();
84 :
85 : void setGeneral();
86 :
87 : /// Set whether each parameter is fixed.
88 : /** Sets the parameter indices appropriately.
89 : */
90 : void setFixed( bool G0, bool G, bool x0, bool y0, bool sigma_x, bool sigma_y, bool theta );
91 :
92 : realT G0( realT *p );
93 :
94 : void G0( realT *p, realT nG0 );
95 :
96 : realT G( realT *p );
97 :
98 : void G( realT *p, realT nG );
99 :
100 : realT x0( realT *p );
101 :
102 : void x0( realT *p, realT nx0 );
103 :
104 : realT y0( realT *p );
105 :
106 : void y0( realT *p, realT ny0 );
107 :
108 : realT sigma_x( realT *p );
109 :
110 : void sigma_x( realT *p, realT nsigma_x );
111 :
112 : realT sigma_y( realT *p );
113 :
114 : void sigma_y( realT *p, realT nsigma_y );
115 :
116 : realT theta( realT *p );
117 :
118 : void theta( realT *p, realT ntheta );
119 :
120 : realT sigma( realT *p );
121 :
122 : void sigma( realT *p, realT nsigma );
123 :
124 : realT a( realT *p );
125 :
126 : void a( realT *p, realT na );
127 :
128 : realT b( realT *p );
129 :
130 : void b( realT *p, realT nb );
131 :
132 : realT c( realT *p );
133 :
134 : void c( realT *p, realT nc );
135 :
136 : int nparams();
137 : };
138 :
139 : template <typename realT>
140 0 : void array2FitGaussian2D<realT>::setSymmetric()
141 : {
142 0 : m_maxNparams = 5;
143 0 : }
144 :
145 : template <typename realT>
146 21 : void array2FitGaussian2D<realT>::setGeneral()
147 : {
148 21 : m_maxNparams = 7;
149 21 : }
150 :
151 : template <typename realT>
152 0 : void array2FitGaussian2D<realT>::setFixed( bool G0, bool G, bool x0, bool y0, bool sigma_x, bool sigma_y, bool theta )
153 : {
154 0 : int idx = 0;
155 :
156 0 : if( G0 )
157 0 : m_G0_idx = -1;
158 : else
159 0 : m_G0_idx = idx++;
160 :
161 0 : if( G )
162 0 : m_G_idx = -1;
163 : else
164 0 : m_G_idx = idx++;
165 :
166 0 : if( x0 )
167 0 : m_x0_idx = -1;
168 : else
169 0 : m_x0_idx = idx++;
170 :
171 0 : if( y0 )
172 0 : m_y0_idx = -1;
173 : else
174 0 : m_y0_idx = idx++;
175 :
176 0 : if( m_maxNparams == 5 )
177 : {
178 0 : if( sigma_x )
179 0 : m_sigma_idx = -1;
180 : else
181 0 : m_sigma_idx = idx++;
182 : }
183 0 : else if( sigma_x && sigma_y && theta )
184 : {
185 0 : m_sigma_x_idx = -1;
186 0 : m_sigma_y_idx = -1;
187 0 : m_theta_idx = -1;
188 : }
189 : else
190 : {
191 0 : if( sigma_x || sigma_y || theta )
192 : {
193 0 : internal::mxlib_error_report(error_t::notimpl, "cannot fix sigma_x, sigma_y, and theta separately" );
194 : }
195 :
196 0 : m_sigma_x_idx = idx++;
197 0 : m_sigma_y_idx = idx++;
198 0 : m_theta_idx = idx++;
199 : }
200 :
201 0 : m_nparams = idx;
202 0 : }
203 :
204 : template <typename realT>
205 301 : realT array2FitGaussian2D<realT>::G0( realT *p )
206 : {
207 301 : if( m_G0_idx < 0 )
208 : {
209 0 : return m_G0;
210 : }
211 : else
212 : {
213 301 : return p[m_G0_idx];
214 : }
215 : }
216 :
217 : template <typename realT>
218 8 : void array2FitGaussian2D<realT>::G0( realT *p, realT nG0 )
219 : {
220 8 : if( m_G0_idx < 0 )
221 : {
222 0 : m_G0 = nG0;
223 : }
224 : else
225 : {
226 8 : p[m_G0_idx] = nG0;
227 : }
228 8 : }
229 :
230 : template <typename realT>
231 309 : realT array2FitGaussian2D<realT>::G( realT *p )
232 : {
233 309 : if( m_G_idx < 0 )
234 : {
235 0 : return m_G;
236 : }
237 : else
238 : {
239 309 : return p[m_G_idx];
240 : }
241 : }
242 :
243 : template <typename realT>
244 8 : void array2FitGaussian2D<realT>::G( realT *p, realT nG )
245 : {
246 8 : if( m_G_idx < 0 )
247 : {
248 0 : m_G = nG;
249 : }
250 : else
251 : {
252 8 : p[m_G_idx] = nG;
253 : }
254 8 : }
255 :
256 : template <typename realT>
257 309 : realT array2FitGaussian2D<realT>::x0( realT *p )
258 : {
259 309 : if( m_x0_idx < 0 )
260 : {
261 0 : return m_x0;
262 : }
263 : else
264 : {
265 309 : return p[m_x0_idx];
266 : }
267 : }
268 :
269 : template <typename realT>
270 8 : void array2FitGaussian2D<realT>::x0( realT *p, realT nx0 )
271 : {
272 8 : if( m_x0_idx < 0 )
273 : {
274 0 : m_x0 = nx0;
275 : }
276 : else
277 : {
278 8 : p[m_x0_idx] = nx0;
279 : }
280 8 : }
281 :
282 : template <typename realT>
283 309 : realT array2FitGaussian2D<realT>::y0( realT *p )
284 : {
285 309 : if( m_y0_idx < 0 )
286 : {
287 0 : return m_y0;
288 : }
289 : else
290 : {
291 309 : return p[m_y0_idx];
292 : }
293 : }
294 :
295 : template <typename realT>
296 8 : void array2FitGaussian2D<realT>::y0( realT *p, realT ny0 )
297 : {
298 8 : if( m_y0_idx < 0 )
299 : {
300 0 : m_y0 = ny0;
301 : }
302 : else
303 : {
304 8 : p[m_y0_idx] = ny0;
305 : }
306 8 : }
307 :
308 : template <typename realT>
309 8 : realT array2FitGaussian2D<realT>::sigma_x( realT *p )
310 : {
311 8 : if( m_sigma_x_idx < 0 )
312 : {
313 0 : return m_sigma_x;
314 : }
315 : else
316 : {
317 8 : return p[m_sigma_x_idx];
318 : }
319 : }
320 :
321 : template <typename realT>
322 24 : void array2FitGaussian2D<realT>::sigma_x( realT *p, realT nsigma_x )
323 : {
324 24 : if( m_sigma_x_idx < 0 )
325 : {
326 0 : m_sigma_x = nsigma_x;
327 : }
328 : else
329 : {
330 24 : p[m_sigma_x_idx] = nsigma_x;
331 : }
332 24 : }
333 :
334 : template <typename realT>
335 8 : realT array2FitGaussian2D<realT>::sigma_y( realT *p )
336 : {
337 8 : if( m_sigma_y_idx < 0 )
338 : {
339 0 : return m_sigma_y;
340 : }
341 : else
342 : {
343 8 : return p[m_sigma_y_idx];
344 : }
345 : }
346 :
347 : template <typename realT>
348 24 : void array2FitGaussian2D<realT>::sigma_y( realT *p, realT nsigma_y )
349 : {
350 24 : if( m_sigma_y_idx < 0 )
351 : {
352 0 : m_sigma_y = nsigma_y;
353 : }
354 : else
355 : {
356 24 : p[m_sigma_y_idx] = nsigma_y;
357 : }
358 24 : }
359 :
360 : template <typename realT>
361 8 : realT array2FitGaussian2D<realT>::theta( realT *p )
362 : {
363 8 : if( m_theta_idx < 0 )
364 : {
365 0 : return m_theta;
366 : }
367 : else
368 : {
369 8 : return p[m_theta_idx];
370 : }
371 : }
372 :
373 : template <typename realT>
374 24 : void array2FitGaussian2D<realT>::theta( realT *p, realT ntheta )
375 : {
376 24 : if( m_theta_idx < 0 )
377 : {
378 0 : m_theta = ntheta;
379 : }
380 : else
381 : {
382 24 : p[m_theta_idx] = ntheta;
383 : }
384 24 : }
385 :
386 : template <typename realT>
387 0 : realT array2FitGaussian2D<realT>::sigma( realT *p )
388 : {
389 0 : if( m_sigma_idx < 0 )
390 : {
391 0 : return m_sigma;
392 : }
393 : else
394 : {
395 0 : return p[m_sigma_idx];
396 : }
397 : }
398 :
399 : template <typename realT>
400 0 : void array2FitGaussian2D<realT>::sigma( realT *p, realT nsigma )
401 : {
402 0 : if( m_sigma_idx < 0 )
403 : {
404 0 : m_sigma = nsigma;
405 : }
406 : else
407 : {
408 0 : p[m_sigma_idx] = nsigma;
409 : }
410 0 : }
411 :
412 : template <typename realT>
413 317 : realT array2FitGaussian2D<realT>::a( realT *p )
414 : {
415 : // This aliases sigma_x after param normalization
416 317 : if( m_sigma_x_idx < 0 )
417 : {
418 0 : return m_a;
419 : }
420 : else
421 : {
422 317 : return p[m_sigma_x_idx];
423 : }
424 : }
425 :
426 : template <typename realT>
427 8 : void array2FitGaussian2D<realT>::a( realT *p, realT na )
428 : {
429 : // This aliases sigma_x after param normalization
430 8 : if( m_sigma_x_idx < 0 )
431 : {
432 0 : m_a = na;
433 : }
434 : else
435 : {
436 8 : p[m_sigma_x_idx] = na;
437 : }
438 8 : }
439 :
440 : template <typename realT>
441 317 : realT array2FitGaussian2D<realT>::b( realT *p )
442 : {
443 : // This aliases sigma_y after param normalization
444 317 : if( m_sigma_y_idx < 0 )
445 : {
446 0 : return m_b;
447 : }
448 : else
449 : {
450 317 : return p[m_sigma_y_idx];
451 : }
452 : }
453 :
454 : template <typename realT>
455 8 : void array2FitGaussian2D<realT>::b( realT *p, realT nb )
456 : {
457 : // This aliases sigma_y after param normalization
458 8 : if( m_sigma_y_idx < 0 )
459 : {
460 0 : m_b = nb;
461 : }
462 : else
463 : {
464 8 : p[m_sigma_y_idx] = nb;
465 : }
466 8 : }
467 :
468 : template <typename realT>
469 317 : realT array2FitGaussian2D<realT>::c( realT *p )
470 : {
471 : // This aliases theta after param normalization
472 317 : if( m_theta_idx < 0 )
473 : {
474 0 : return m_c;
475 : }
476 : else
477 : {
478 317 : return p[m_theta_idx];
479 : }
480 : }
481 :
482 : template <typename realT>
483 8 : void array2FitGaussian2D<realT>::c( realT *p, realT nc )
484 : {
485 : // This aliases theta after param normalization
486 8 : if( m_theta_idx < 0 )
487 : {
488 0 : m_c = nc;
489 : }
490 : else
491 : {
492 8 : p[m_theta_idx] = nc;
493 : }
494 8 : }
495 :
496 : template <typename realT>
497 21 : int array2FitGaussian2D<realT>::nparams()
498 : {
499 21 : return m_nparams;
500 : }
501 :
502 : } // namespace fit
503 : } // namespace math
504 :
505 : } // namespace mx
506 :
507 : #endif // math_fit_array2FitGaussian2D_hpp
|
