PCL
PSFFit.h
Go to the documentation of this file.
1 // ____ ______ __
2 // / __ \ / ____// /
3 // / /_/ // / / /
4 // / ____// /___ / /___ PixInsight Class Library
5 // /_/ \____//_____/ PCL 2.1.19
6 // ----------------------------------------------------------------------------
7 // pcl/PSFFit.h - Released 2019-11-07T10:59:34Z
8 // ----------------------------------------------------------------------------
9 // This file is part of the PixInsight Class Library (PCL).
10 // PCL is a multiplatform C++ framework for development of PixInsight modules.
11 //
12 // Copyright (c) 2003-2019 Pleiades Astrophoto S.L. All Rights Reserved.
13 //
14 // Redistribution and use in both source and binary forms, with or without
15 // modification, is permitted provided that the following conditions are met:
16 //
17 // 1. All redistributions of source code must retain the above copyright
18 // notice, this list of conditions and the following disclaimer.
19 //
20 // 2. All redistributions in binary form must reproduce the above copyright
21 // notice, this list of conditions and the following disclaimer in the
22 // documentation and/or other materials provided with the distribution.
23 //
24 // 3. Neither the names "PixInsight" and "Pleiades Astrophoto", nor the names
25 // of their contributors, may be used to endorse or promote products derived
26 // from this software without specific prior written permission. For written
27 // permission, please contact info@pixinsight.com.
28 //
29 // 4. All products derived from this software, in any form whatsoever, must
30 // reproduce the following acknowledgment in the end-user documentation
31 // and/or other materials provided with the product:
32 //
33 // "This product is based on software from the PixInsight project, developed
34 // by Pleiades Astrophoto and its contributors (http://pixinsight.com/)."
35 //
36 // Alternatively, if that is where third-party acknowledgments normally
37 // appear, this acknowledgment must be reproduced in the product itself.
38 //
39 // THIS SOFTWARE IS PROVIDED BY PLEIADES ASTROPHOTO AND ITS CONTRIBUTORS
40 // "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
41 // TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
42 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL PLEIADES ASTROPHOTO OR ITS
43 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
44 // EXEMPLARY OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, BUSINESS
45 // INTERRUPTION; PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; AND LOSS OF USE,
46 // DATA OR PROFITS) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
47 // CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
48 // ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
49 // POSSIBILITY OF SUCH DAMAGE.
50 // ----------------------------------------------------------------------------
51 
52 #ifndef __PCL_PSFFit_h
53 #define __PCL_PSFFit_h
54 
56 
57 #include <pcl/Defs.h>
58 #include <pcl/Diagnostics.h>
59 
60 #include <pcl/ImageVariant.h>
61 #include <pcl/Matrix.h>
62 #include <pcl/String.h>
63 
64 namespace pcl
65 {
66 
67 // ----------------------------------------------------------------------------
68 
87 namespace PSFunction
88 {
89  enum value_type
90  {
91  Invalid = -1, // Represents an invalid or unsupported PSF type
92  Gaussian, // Gaussian PSF
93  Moffat, // Moffat PSF with a fitted beta parameter
94  MoffatA, // Moffat PSF with fixed beta=10
95  Moffat8, // Moffat PSF with fixed beta=8
96  Moffat6, // Moffat PSF with fixed beta=6
97  Moffat4, // Moffat PSF with fixed beta=4
98  Moffat25, // Moffat PSF with fixed beta=2.5
99  Moffat15, // Moffat PSF with fixed beta=1.5
100  Lorentzian, // Lorentzian PSF, equivalent to a Moffat PSF with fixed beta=1
101  VariableShape, // Variable shape PSF
102 
103  NumberOfFunctions,
104 
105  Default = Gaussian
106  };
107 }
108 
109 // ----------------------------------------------------------------------------
110 
126 namespace PSFFitStatus
127 {
128  enum value_type
129  {
130  Invalid = -1,
131  NotFitted,
132  FittedOk,
133  BadParameters,
134  NoSolution,
135  NoConvergence,
136  InaccurateSolution,
137  UnknownError
138  };
139 }
140 
141 // ----------------------------------------------------------------------------
142 
147 struct PSFData
148 {
152  typedef PSFunction::value_type psf_function;
153 
157  typedef PSFFitStatus::value_type psf_fit_status;
158 
159  psf_function function = PSFunction::Invalid;
160  bool circular = false;
161  psf_fit_status status = PSFFitStatus::NotFitted;
162  double B = 0;
163  double A = 0;
164  DPoint c0 = 0.0;
165  DPoint q0 = 0.0;
166  bool celestial = false;
167  double sx = 0;
168  double sy = 0;
169  double theta = 0;
170  double beta = 0;
171  double mad = 0;
172 
176  PSFData() = default;
177 
181  PSFData( const PSFData& ) = default;
182 
186  PSFData( PSFData&& ) = default;
187 
191  PSFData& operator =( const PSFData& ) = default;
192 
196  PSFData& operator =( PSFData&& ) = default;
197 
203  operator bool() const
204  {
205  return status == PSFFitStatus::FittedOk;
206  }
207 
213  {
214  switch ( function )
215  {
216  case PSFunction::Gaussian: return "Gaussian";
217  case PSFunction::Moffat: return "Moffat";
218  case PSFunction::MoffatA: return "Moffat10";
219  case PSFunction::Moffat8: return "Moffat8";
220  case PSFunction::Moffat6: return "Moffat6";
221  case PSFunction::Moffat4: return "Moffat4";
222  case PSFunction::Moffat25: return "Moffat25";
223  case PSFunction::Moffat15: return "Moffat15";
224  case PSFunction::Lorentzian: return "Lorentzian";
225  case PSFunction::VariableShape: return "VariableShape";
226  default: return "Unknown"; // ?!
227  }
228  }
229 
234  {
235  switch ( status )
236  {
237  case PSFFitStatus::NotFitted: return "Not fitted";
238  case PSFFitStatus::FittedOk: return "Fitted Ok";
239  case PSFFitStatus::BadParameters: return "Bad parameters";
240  case PSFFitStatus::NoSolution: return "No solution";
241  case PSFFitStatus::NoConvergence: return "No convergence";
242  case PSFFitStatus::InaccurateSolution: return "Inaccurate solution";
243  default:
244  case PSFFitStatus::UnknownError: return "Unknown error";
245  }
246  }
247 
256  double FWHMx() const
257  {
258  return FWHM( function, sx, beta );
259  }
260 
269  double FWHMy() const
270  {
271  return FWHM( function, sy, beta );
272  }
273 
288  DRect Bounds() const
289  {
290  double d = (circular ? FWHMx() : Max( FWHMx(), FWHMy() ))/2;
291  return DRect( c0.x-d, c0.y-d, c0.x+d, c0.y+d );
292  }
293 
300  {
301  return celestial;
302  }
303 
308  void ToImage( Image& ) const;
309 
325  static double FWHM( psf_function function, double sigma, double beta = 2 )
326  {
327  PCL_PRECONDITION( beta > 0 )
328  switch ( function )
329  {
330  case PSFunction::Gaussian: return 2.3548200450309493 * sigma;
331  case PSFunction::Moffat: return 2 * sigma * Sqrt( Pow2( 1/beta ) - 1 );
332  case PSFunction::MoffatA: return 0.5358113941912513 * sigma;
333  case PSFunction::Moffat8: return 0.6016900619596693 * sigma;
334  case PSFunction::Moffat6: return 0.6998915581984769 * sigma;
335  case PSFunction::Moffat4: return 0.8699588840921645 * sigma;
336  case PSFunction::Moffat25: return 1.1305006161394060 * sigma;
337  case PSFunction::Moffat15: return 1.5328418730817597 * sigma;
338  case PSFunction::Lorentzian: return 2 * sigma;
339  case PSFunction::VariableShape: return 2 * sigma * Pow( beta*0.6931471805599453, 1/beta );
340  default: return 0; // ?!
341  }
342  }
343 };
344 
345 // ----------------------------------------------------------------------------
346 
351 class PSFFit
352 {
353 public:
354 
358  typedef PSFData::psf_function psf_function;
359 
363  typedef PSFData::psf_fit_status psf_fit_status;
364 
369 
400  PSFFit( const ImageVariant& image,
401  const DPoint& center, const DRect& rect,
402  psf_function function = PSFunction::Gaussian, bool circular = false );
403 
407  PSFFit( const PSFFit& ) = default;
408 
412  PSFFit( PSFFit&& ) = default;
413 
417  PSFFit& operator =( const PSFFit& x ) = default;
418 
422  PSFFit& operator =( PSFFit&& x ) = default;
423 
428  operator bool() const
429  {
430  return psf;
431  }
432 
433 private:
434 
435  Matrix S; // matrix of sampled data
436  Vector P; // vector of function parameters
437  mutable double beta0;
438 
439  Vector GoodnessOfFit( psf_function, bool circular ) const;
440 
441  friend class PSFFitEngine;
442 };
443 
444 // ----------------------------------------------------------------------------
445 
446 } // pcl
447 
448 #endif // __PCL_PSFFit_h
449 
450 // ----------------------------------------------------------------------------
451 // EOF pcl/PSFFit.h - Released 2019-11-07T10:59:34Z
PSF fit parameters.
Definition: PSFFit.h:147
String FunctionToString() const
Definition: PSFFit.h:212
PSFData psf
Definition: PSFFit.h:368
PCL root namespace.
Definition: AbstractImage.h:76
constexpr T Pow2(T x)
Definition: Math.h:1607
Complex< T > Sqrt(const Complex< T > &c)
Definition: Complex.h:665
String StatusToString() const
Definition: PSFFit.h:233
Numerical Point Spread Function (PSF) fit to a source in an image.
Definition: PSFFit.h:351
PSFData::psf_function psf_function
Definition: PSFFit.h:358
64-bit floating-point rectangle in the R^2 space.
Acts like a union for all types of images in PCL, with optional class-wide ownership of transported i...
Definition: ImageVariant.h:317
64-bit floating point real matrix.
Unicode (UTF-16) string.
Definition: String.h:7911
constexpr const T & Max(const T &a, const T &b)
Definition: Utility.h:119
PSFFitStatus::value_type psf_fit_status
Definition: PSFFit.h:157
bool HasCelestialCoordinates() const
Definition: PSFFit.h:299
PSFData::psf_fit_status psf_fit_status
Definition: PSFFit.h:363
64-bit floating point real vector.
64-bit floating-point point in the R^2 space.
PSFunction::value_type psf_function
Definition: PSFFit.h:152
double FWHMx() const
Definition: PSFFit.h:256
DRect Bounds() const
Definition: PSFFit.h:288
double FWHMy() const
Definition: PSFFit.h:269
static double FWHM(psf_function function, double sigma, double beta=2)
Definition: PSFFit.h:325
32-bit floating point real image.
Complex< T1 > Pow(const Complex< T1 > &c, T2 x)
Definition: Complex.h:738