PCL
PSFFit.h
Go to the documentation of this file.
1 // ____ ______ __
2 // / __ \ / ____// /
3 // / /_/ // / / /
4 // / ____// /___ / /___ PixInsight Class Library
5 // /_/ \____//_____/ PCL 2.6.5
6 // ----------------------------------------------------------------------------
7 // pcl/PSFFit.h - Released 2024-01-13T15:47:58Z
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-2024 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 (https://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 = 0, // 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  NumberOfFunctions,
103  Auto = 999, // Automatic selection for optimization, process-dependent.
104  Default = Gaussian
105  };
106 }
107 
108 // ----------------------------------------------------------------------------
109 
125 namespace PSFFitStatus
126 {
127  enum value_type
128  {
129  Invalid = -1,
130  NotFitted = 0,
131  FittedOk,
132  BadParameters,
133  NoSolution,
134  NoConvergence,
135  InaccurateSolution,
136  UnknownError
137  };
138 }
139 
140 // ----------------------------------------------------------------------------
141 
146 struct PSFData
147 {
151  using psf_function = PSFunction::value_type;
152 
156  using psf_fit_status = PSFFitStatus::value_type;
157 
163 
164  psf_function function = PSFunction::Invalid;
165  bool circular = false;
166  psf_fit_status status = PSFFitStatus::NotFitted;
167  bool celestial = false;
168  double B = 0;
169  double A = 0;
170  DPoint b0 = 0.0;
171  DPoint c0 = 0.0;
172  DPoint q0 = 0.0;
173  double sx = 0;
174  double sy = 0;
175  double theta = 0;
176  double beta = 0;
177  double flux = 0;
178  double signal = 0;
179  unsigned signalCount = 0;
180  double mad = 0;
185  PSFData() = default;
186 
190  PSFData( const PSFData& ) = default;
191 
195  PSFData( PSFData&& ) = default;
196 
200  PSFData& operator =( const PSFData& ) = default;
201 
205  PSFData& operator =( PSFData&& ) = default;
206 
212  operator bool() const
213  {
214  return status == PSFFitStatus::FittedOk;
215  }
216 
222  operator double() const
223  {
224  return (signal != 0) ? signal : flux;
225  }
226 
233  double operator []( int i ) const noexcept
234  {
235  return (i == 0) ? c0.x : c0.y;
236  }
237 
242  String FunctionName() const;
243 
247  String StatusText() const;
248 
257  double FWHMx() const
258  {
259  return FWHM( function, sx, beta );
260  }
261 
270  double FWHMy() const
271  {
272  return FWHM( function, sy, beta );
273  }
274 
283  double FWTMx() const
284  {
285  return FWTM( function, sx, beta );
286  }
287 
296  double FWTMy() const
297  {
298  return FWTM( function, sy, beta );
299  }
300 
305  double Volume() const
306  {
307  return A*Volume( function, sx, sy, beta );
308  }
309 
324  DRect Bounds() const
325  {
326  double d = (circular ? FWHMx() : Max( FWHMx(), FWHMy() ))/2;
327  return DRect( c0.x-d, c0.y-d, c0.x+d, c0.y+d );
328  }
329 
346  {
347  double d = (circular ? FWHMx() : Max( FWHMx(), FWHMy() ))/2;
348  return DRect( c0.x-d, c0.y-d, c0.x+d, c0.y+d );
349  }
350 
357  {
358  return celestial;
359  }
360 
365  void ToImage( Image& ) const;
366 
382  static double FWHM( psf_function function, double sigma, double beta = 2 )
383  {
384  PCL_PRECONDITION( beta > 0 )
385  switch ( function )
386  {
387  case PSFunction::Gaussian: return 2.3548200450309493 * sigma;
388  case PSFunction::Moffat: return 2 * sigma * Sqrt( Pow2( 1/beta ) - 1 );
389  case PSFunction::MoffatA: return 0.5358113941912513 * sigma;
390  case PSFunction::Moffat8: return 0.6016900619596693 * sigma;
391  case PSFunction::Moffat6: return 0.6998915581984769 * sigma;
392  case PSFunction::Moffat4: return 0.8699588840921645 * sigma;
393  case PSFunction::Moffat25: return 1.1305006161394060 * sigma;
394  case PSFunction::Moffat15: return 1.5328418730817597 * sigma;
395  case PSFunction::Lorentzian: return 2 * sigma;
396  case PSFunction::VariableShape: return 2 * sigma * Pow( beta*0.6931471805599453, 1/beta );
397  default: return 0; // ?!
398  }
399  }
400 
416  static double FWTM( psf_function function, double sigma, double beta = 2 )
417  {
418  PCL_PRECONDITION( beta > 0 )
419  switch ( function )
420  {
421  case PSFunction::Gaussian: return 4.291932052578694 * sigma;
422  case PSFunction::Moffat: return 2 * sigma * Sqrt( Pow( 10.0, 1/beta ) - 1 );
423  case PSFunction::MoffatA: return 1.017694279819175 * sigma;
424  case PSFunction::Moffat8: return 1.155026289161115 * sigma;
425  case PSFunction::Moffat6: return 1.367917055412454 * sigma;
426  case PSFunction::Moffat4: return 1.764402913213331 * sigma;
427  case PSFunction::Moffat25: return 2.459175822514185 * sigma;
428  case PSFunction::Moffat15: return 3.816589489904712 * sigma;
429  case PSFunction::Lorentzian: return 6 * sigma;
430  case PSFunction::VariableShape: return 2 * sigma * Pow( beta*2.302585092994045, 1/beta );
431  default: return 0; // ?!
432  }
433  }
434 
456  static double Volume( psf_function function, double sigma_x, double sigma_y, double beta = 2 )
457  {
458  PCL_PRECONDITION( beta > 0 )
459  PCL_PRECONDITION( function != PSFunction::Lorentzian && (function != PSFunction::Moffat || beta > 1) )
460  switch ( function )
461  {
462  case PSFunction::Gaussian: return 6.2831853071795862 * sigma_x*sigma_y;
463  case PSFunction::Moffat: return 3.1415926535897931 * sigma_x*sigma_y/(beta - 1);
464  case PSFunction::MoffatA: return 0.3490658503988659 * sigma_x*sigma_y;
465  case PSFunction::Moffat8: return 0.4487989505128276 * sigma_x*sigma_y;
466  case PSFunction::Moffat6: return 0.6283185307179586 * sigma_x*sigma_y;
467  case PSFunction::Moffat4: return 1.0471975511965976 * sigma_x*sigma_y;
468  case PSFunction::Moffat25: return 2.0943951023931953 * sigma_x*sigma_y;
469  case PSFunction::Moffat15: return 6.2831853071795862 * sigma_x*sigma_y;
470  default: return 0; // ?!
471  }
472  }
473 };
474 
475 // ----------------------------------------------------------------------------
476 
482 class PSFFit
483 {
484 public:
485 
489  using psf_function = PSFData::psf_function;
490 
494  using psf_fit_status = PSFData::psf_fit_status;
495 
500 
567  PSFFit( const ImageVariant& image,
568  const DPoint& center, const DRect& rect,
569  psf_function function = PSFunction::Gaussian, bool circular = false,
570  float betaMin = 1.0F, float betaMax = 4.0F,
571  double tolerance = 1.0e-08, float bkgMaxVar = 0.1F, float growthForFlux = 1.0F );
572 
576  PSFFit( const PSFFit& ) = default;
577 
581  PSFFit( PSFFit&& ) = default;
582 
586  PSFFit& operator =( const PSFFit& x ) = default;
587 
591  PSFFit& operator =( PSFFit&& x ) = default;
592 
597  operator bool() const
598  {
599  return psf;
600  }
601 
602 private:
603 
604  Matrix S; // matrix of sampled pixel data
605  Vector P; // vector of function parameters
606 
607  // The initial local background measured on the sampling region and the
608  // maximum allowed relative difference, for stabilization of local
609  // background PSF parameters.
610  double m_bkg;
611  float m_bkgMaxVar;
612 
613  // Growing factor in units of the Full Width at Tenth Maximum (FWTM) for
614  // extension/contraction of the PSF flux measurement region.
615  float m_growthForFlux = 1.0F;
616 
617  // Keep track of successive beta values in L-M iterations for stabilization
618  // of shape parameters. This guarantees convergence for Moffat functions.
619  mutable float m_beta;
620 
621  Vector GoodnessOfFit( psf_function, bool circular, bool test = false ) const;
622 
623  friend class PSFFitEngine;
624 };
625 
626 // ----------------------------------------------------------------------------
627 
628 } // pcl
629 
630 #endif // __PCL_PSFFit_h
631 
632 // ----------------------------------------------------------------------------
633 // EOF pcl/PSFFit.h - Released 2024-01-13T15:47:58Z
pcl::PSFData::Volume
double Volume() const
Definition: PSFFit.h:305
pcl
PCL root namespace.
Definition: AbstractImage.h:76
pcl::PSFData::A
double A
Function amplitude (or estimated maximum) in pixel value units.
Definition: PSFFit.h:169
pcl::PSFData::PSFData
PSFData()=default
pcl::PSFData
PSF fit parameters.
Definition: PSFFit.h:146
pcl::PSFData::mad
double mad
Definition: PSFFit.h:180
pcl::PSFData::StatusText
String StatusText() const
pcl::PSFFit::PSFFit
PSFFit(const ImageVariant &image, const DPoint &center, const DRect &rect, psf_function function=PSFunction::Gaussian, bool circular=false, float betaMin=1.0F, float betaMax=4.0F, double tolerance=1.0e-08, float bkgMaxVar=0.1F, float growthForFlux=1.0F)
pcl::PSFData::FWTM
static double FWTM(psf_function function, double sigma, double beta=2)
Definition: PSFFit.h:416
pcl::PSFData::theta
double theta
Rotation angle of the sx axis in degrees, in the [0,180) range.
Definition: PSFFit.h:175
pcl::String
Unicode (UTF-16) string.
Definition: String.h:8112
pcl::Max
constexpr const T & Max(const T &a, const T &b) noexcept
Definition: Utility.h:119
pcl::GenericPoint
A generic point in the two-dimensional space.
Definition: Point.h:99
pcl::GenericMatrix
Generic dynamic matrix of arbitrary dimensions.
Definition: Matrix.h:122
Matrix.h
pcl::PSFData::FWTMx
double FWTMx() const
Definition: PSFFit.h:283
pcl::PSFData::Volume
static double Volume(psf_function function, double sigma_x, double sigma_y, double beta=2)
Definition: PSFFit.h:456
pcl::PSFData::sx
double sx
Function width in pixels on the X axis, sx >= sy.
Definition: PSFFit.h:173
pcl::PSFFit::operator=
PSFFit & operator=(const PSFFit &x)=default
pcl::PSFData::operator[]
double operator[](int i) const noexcept
Definition: PSFFit.h:233
pcl::PSFData::psf_function
PSFunction::value_type psf_function
Definition: PSFFit.h:151
pcl::PSFData::FWHMy
double FWHMy() const
Definition: PSFFit.h:270
pcl::PSFData::FWHMx
double FWHMx() const
Definition: PSFFit.h:257
pcl::PSFData::q0
DPoint q0
Centroid equatorial coordinates, when celestial=true.
Definition: PSFFit.h:172
pcl::PSFData::ToImage
void ToImage(Image &) const
pcl::Pow
Complex< T1 > Pow(const Complex< T1 > &c, T2 x) noexcept
Definition: Complex.h:747
pcl::PSFData::component
DPoint::component component
Definition: PSFFit.h:162
pcl::PSFData::flux
double flux
Total flux above the local background, measured over the rectangular fitting region.
Definition: PSFFit.h:177
pcl::PSFData::FWTMy
double FWTMy() const
Definition: PSFFit.h:296
pcl::PSFData::celestial
bool celestial
True iff equatorial coordinates are available.
Definition: PSFFit.h:167
pcl::PSFData::beta
double beta
Moffat beta or shape parameter (dimensionless).
Definition: PSFFit.h:176
pcl::PSFData::circular
bool circular
Circular or elliptical PSF.
Definition: PSFFit.h:165
pcl::PSFData::c0
DPoint c0
Fitted centroid position in image coordinates.
Definition: PSFFit.h:171
pcl::GenericRectangle
A generic rectangle in the two-dimensional space.
Definition: Rectangle.h:313
pcl::PSFData::psf_fit_status
PSFFitStatus::value_type psf_fit_status
Definition: PSFFit.h:156
pcl::PSFData::HasCelestialCoordinates
bool HasCelestialCoordinates() const
Definition: PSFFit.h:356
pcl::Sqrt
Complex< T > Sqrt(const Complex< T > &c) noexcept
Definition: Complex.h:674
pcl::PSFData::FWTMBounds
DRect FWTMBounds() const
Definition: PSFFit.h:345
pcl::PSFData::FWHM
static double FWHM(psf_function function, double sigma, double beta=2)
Definition: PSFFit.h:382
pcl::PSFData::B
double B
Local background estimate in pixel value units.
Definition: PSFFit.h:168
pcl::GenericPoint::x
component x
Abscissa (horizontal, or X-axis coordinate).
Definition: Point.h:111
pcl::PSFData::signal
double signal
Total flux above the local background, measured over the elliptical PSF region.
Definition: PSFFit.h:178
pcl::PSFData::b0
DPoint b0
Barycenter position (initial star detection position) in image coordinates.
Definition: PSFFit.h:170
pcl::PSFData::sy
double sy
Function width in pixels on the Y axis, sx >= sy.
Definition: PSFFit.h:174
pcl::PSFData::Bounds
DRect Bounds() const
Definition: PSFFit.h:324
pcl::PSFData::signalCount
unsigned signalCount
Number of pixels used for signal evaluation.
Definition: PSFFit.h:179
pcl::PSFData::operator=
PSFData & operator=(const PSFData &)=default
String.h
pcl::PSFFit
Numerical Point Spread Function (PSF) fit to a source in an image.
Definition: PSFFit.h:482
pcl::GenericPoint::y
component y
Ordinate (vertical, or Y-axis coordinate).
Definition: Point.h:112
pcl::Pow2
constexpr T Pow2(T x) noexcept
Definition: Math.h:1717
pcl::PSFData::FunctionName
String FunctionName() const
pcl::ImageVariant
Acts like a union for all types of images in PCL, with optional class-wide ownership of transported i...
Definition: ImageVariant.h:321
Defs.h
pcl::PSFFit::psf_function
PSFData::psf_function psf_function
Definition: PSFFit.h:489
pcl::PSFFit::psf
PSFData psf
Definition: PSFFit.h:499
pcl::GenericVector< double >
ImageVariant.h
pcl::GenericPoint::component
T component
Definition: Point.h:106
pcl::GenericImage
Implements a generic, two-dimensional, shared or local image.
Definition: Image.h:277
pcl::PSFData::status
psf_fit_status status
Status code (PSFFitStatus namespace).
Definition: PSFFit.h:166