Line data Source code
1 : /** \file fitEmpirical.hpp
2 : * \author Jared R. Males
3 : * \brief Tools for fitting empirical functions to data.
4 : * \ingroup fitting_files
5 : *
6 : */
7 :
8 : //***********************************************************************//
9 : // Copyright 2024 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 fitEmpirical_hpp
28 : #define fitEmpirical_hpp
29 :
30 : #include "levmarInterface.hpp"
31 :
32 : #include "../../improc/eigenImage.hpp"
33 : #include "../../improc/imageTransforms.hpp"
34 :
35 : namespace mx
36 : {
37 : namespace math
38 : {
39 : namespace fit
40 : {
41 :
42 : // forward
43 : template <typename realT>
44 : struct array2FitEmpirical;
45 :
46 : // forward
47 : template <typename _realT>
48 : struct empirical2D_sym_fitter;
49 :
50 : /// Class to manage fitting a 2D Moffat to data via the \ref levmarInterface
51 : /** Fits \ref gen_math_moffats to a 2-dimensional array of data.
52 : *
53 : * This class allows for treating any of the parameters as fixed.
54 : *
55 : * \tparam fitterT a type meeting the requirements specified in \ref levmarInterface.
56 : *
57 : * \ingroup moffat_peak_fit
58 : *
59 : */
60 : template <typename fitterT>
61 : class fitEmpirical2DGen : public levmarInterface<fitterT>
62 : {
63 :
64 : public:
65 : typedef typename fitterT::realT realT;
66 :
67 : public:
68 : array2FitEmpirical<realT> arr;
69 :
70 0 : void initialize()
71 : {
72 0 : this->allocate_params( arr.nparams() );
73 0 : this->adata = &arr;
74 0 : }
75 :
76 : public:
77 0 : fitEmpirical2DGen()
78 0 : {
79 0 : initialize();
80 0 : }
81 :
82 0 : ~fitEmpirical2DGen()
83 : {
84 0 : }
85 :
86 : /// Set whether each parameter is fixed.
87 : /** Sets the parameter indices appropriately.
88 : */
89 0 : void setFixed( bool scale, ///< [in] if true, then scale will be not be part of the fit
90 : bool dx, ///< [in] if true, then dx will be not be part of the fit
91 : bool dy ///< [in] if true, then dy will be not be part of the fit
92 : )
93 : {
94 0 : arr.setFixed( scale, dx, dy );
95 0 : this->allocate_params( arr.nparams() );
96 0 : }
97 :
98 : /// Set the initial guess for the empirical fit
99 0 : void setGuess( realT scale, ///< [in] the scale factr
100 : realT dx, ///< [in] the x shift
101 : realT dy ///< [in] the y shift
102 : )
103 : {
104 0 : arr.scale( this->p, scale );
105 0 : arr.dx( this->p, dx );
106 0 : arr.dy( this->p, dy );
107 0 : }
108 :
109 : /// Set the data aray.
110 0 : void setArray( const improc::eigenImage<realT> *data, ///< [in] pointer to an nx X ny array of data to be fit
111 : const improc::eigenImage<realT> *ref, ///< [in] pointer to the empirical function to fit
112 : const improc::eigenImage<realT> *weights )
113 : {
114 0 : arr.setup( data, ref, weights );
115 :
116 0 : this->n = arr.m_nx * arr.m_ny;
117 0 : }
118 :
119 : /// Set the data aray.
120 0 : void setArray( const improc::eigenImage<realT> *data, ///< [in] pointer to an nx X ny array of data to be fit
121 : const improc::eigenImage<realT> *ref ///< [in] pointer to the empirical function to fit
122 : )
123 : {
124 0 : setArray( data, ref, nullptr );
125 0 : }
126 :
127 : /// Get the current value of the scale factor.
128 : /**
129 : * \returns the current value of scale
130 : */
131 0 : realT scale()
132 : {
133 0 : return arr.scale( this->p );
134 : }
135 :
136 : /// Get the current value of dx, the x shift
137 : /**
138 : * \returns the current value of dx
139 : */
140 0 : realT dx()
141 : {
142 0 : return arr.dx( this->p );
143 : }
144 :
145 : /// Get the current value of dy, the y shift
146 : /**
147 : * \returns the current value of dy
148 : */
149 0 : realT dy()
150 : {
151 0 : return arr.dy( this->p );
152 : }
153 : };
154 :
155 : /// Wrapper for a native array to pass to \ref levmarInterface, with empirical function fit details.
156 : /** \ingroup moffat_peak_fit
157 : */
158 : template <typename realT>
159 : struct array2FitEmpirical
160 : {
161 : const improc::eigenImage<realT> *m_data{ nullptr }; ///< Pointer to the data array
162 : const improc::eigenImage<realT> *m_ref{ nullptr }; ///< Pointer to the reference image to fit to the data
163 : const improc::eigenImage<realT> *m_weights{ nullptr }; ///< Pointer to the weight image
164 : improc::eigenImage<realT> m_refShifted; ///< Working memory for the shifted reference image
165 :
166 : size_t m_nx{ 0 }; ///< X dimension of the array
167 : size_t m_ny{ 0 }; ///< Y dimension of the array
168 :
169 : realT m_scale{ 0 };
170 : realT m_dx{ 0 };
171 : realT m_dy{ 0 };
172 :
173 : int m_scale_idx{ 0 };
174 : int m_dx_idx{ 1 };
175 : int m_dy_idx{ 2 };
176 :
177 : int m_nparams = 3;
178 :
179 0 : void setup( const improc::eigenImage<realT> *data,
180 : const improc::eigenImage<realT> *ref,
181 : const improc::eigenImage<realT> *weights )
182 : {
183 0 : m_nx = data->rows();
184 0 : m_ny = data->cols();
185 :
186 0 : if( ref->rows() != m_nx || ref->cols() != m_ny )
187 : {
188 0 : std::cerr << "ref and data not same size\n";
189 0 : exit( -1 );
190 : }
191 :
192 0 : if( weights )
193 : {
194 0 : if( weights->rows() != m_nx || weights->cols() != m_ny )
195 : {
196 0 : std::cerr << "ref and data not same size\n";
197 0 : exit( -1 );
198 : }
199 : }
200 :
201 0 : m_data = data;
202 0 : m_ref = ref;
203 0 : m_weights = weights;
204 0 : m_refShifted.resize( m_nx, m_ny );
205 0 : }
206 :
207 0 : void setup( const improc::eigenImage<realT> *data, const improc::eigenImage<realT> *ref )
208 : {
209 0 : setup( data, ref, nullptr );
210 0 : }
211 :
212 : /// Set whether each parameter is fixed.
213 : /** Sets the parameter indices appropriately.
214 : */
215 0 : void setFixed( bool scale, bool dx, bool dy )
216 : {
217 0 : int idx = 0;
218 :
219 0 : if( scale )
220 0 : m_scale_idx = -1;
221 : else
222 0 : m_scale_idx = idx++;
223 :
224 0 : if( dx )
225 0 : m_dx_idx = -1;
226 : else
227 0 : m_dx_idx = idx++;
228 :
229 0 : if( dy )
230 0 : m_dy_idx = -1;
231 : else
232 0 : m_dy_idx = idx++;
233 :
234 0 : m_nparams = idx;
235 0 : }
236 :
237 0 : realT scale( realT *p )
238 : {
239 0 : if( m_scale_idx < 0 )
240 : {
241 0 : return m_scale;
242 : }
243 : else
244 : {
245 0 : return p[m_scale_idx];
246 : }
247 : }
248 :
249 0 : void scale( realT *p, realT nscale )
250 : {
251 0 : if( m_scale_idx < 0 )
252 : {
253 0 : m_scale = nscale;
254 : }
255 : else
256 : {
257 0 : p[m_scale_idx] = nscale;
258 : }
259 0 : }
260 :
261 0 : realT dx( realT *p )
262 : {
263 0 : if( m_dx_idx < 0 )
264 : {
265 0 : return m_dx;
266 : }
267 : else
268 : {
269 0 : return p[m_dx_idx];
270 : }
271 : }
272 :
273 0 : void dx( realT *p, realT ndx )
274 : {
275 0 : if( m_dx_idx < 0 )
276 : {
277 0 : m_dx = ndx;
278 : }
279 : else
280 : {
281 0 : p[m_dx_idx] = ndx;
282 : }
283 0 : }
284 :
285 0 : realT dy( realT *p )
286 : {
287 0 : if( m_dy_idx < 0 )
288 : {
289 0 : return m_dy;
290 : }
291 : else
292 : {
293 0 : return p[m_dy_idx];
294 : }
295 : }
296 :
297 0 : void dy( realT *p, realT ndy )
298 : {
299 0 : if( m_dy_idx < 0 )
300 : {
301 0 : m_dy = ndy;
302 : }
303 : else
304 : {
305 0 : p[m_dy_idx] = ndy;
306 : }
307 0 : }
308 :
309 0 : int nparams()
310 : {
311 0 : return m_nparams;
312 : }
313 : };
314 :
315 : ///\ref levmarInterface fitter structure for 2D empirical functions.
316 : /** \ingroup moffat_peak_fit
317 : *
318 : */
319 : template <typename _realT>
320 : struct empirical2D_fitter
321 : {
322 : typedef _realT realT;
323 :
324 : static const int nparams = 3;
325 :
326 0 : static void func( realT *p, realT *hx, int m, int n, void *adata )
327 : {
328 0 : array2FitEmpirical<realT> *arr = (array2FitEmpirical<realT> *)adata;
329 :
330 : size_t idx_dat;
331 :
332 0 : idx_dat = 0;
333 :
334 0 : realT scale = arr->scale( p );
335 0 : realT dx = arr->dx( p );
336 0 : realT dy = arr->dy( p );
337 :
338 0 : improc::imageShift( arr->m_refShifted, *( arr->m_ref ), dx, dy, improc::cubicConvolTransform<realT>() );
339 :
340 0 : arr->m_refShifted *= scale;
341 :
342 0 : if( arr->m_weights )
343 : {
344 0 : for( int cc = 0; cc < arr->m_ny; ++cc )
345 : {
346 0 : for( int rr = 0; rr < arr->m_nx; ++rr )
347 : {
348 0 : hx[idx_dat] =
349 0 : ( *arr->m_weights )( rr, cc ) * ( ( *arr->m_data )( rr, cc ) - arr->m_refShifted( rr, cc ) );
350 :
351 0 : ++idx_dat;
352 : }
353 : }
354 : }
355 : else
356 : {
357 0 : for( int cc = 0; cc < arr->m_ny; ++cc )
358 : {
359 0 : for( int rr = 0; rr < arr->m_nx; ++rr )
360 : {
361 0 : hx[idx_dat] = ( ( *( arr->m_data ) )( rr, cc ) - arr->m_refShifted( rr, cc ) );
362 :
363 0 : ++idx_dat;
364 : }
365 : }
366 : }
367 0 : }
368 : };
369 :
370 : /// Alias for the fitEmpirical2D type fitting the symmetric Moffat profile.
371 : /** \ingroup moffat_peak_fit
372 : */
373 : template <typename realT>
374 : using fitEmpirical2D = mx::math::fit::fitEmpirical2DGen<mx::math::fit::empirical2D_fitter<realT>>;
375 :
376 : } // namespace fit
377 : } // namespace math
378 :
379 : } // namespace mx
380 :
381 : #endif // fitEmpirical_hpp
|
