PSFFit.h

Go to the documentation of this file.

//     ____   ______ __
//    / __ \ / ____// /
//   / /_/ // /    / /
//  / ____// /___ / /___   PixInsight Class Library
// /_/     \____//_____/   PCL 2.7.0
// ----------------------------------------------------------------------------
// pcl/PSFFit.h - Released 2024-06-18T15:48:54Z
// ----------------------------------------------------------------------------
// This file is part of the PixInsight Class Library (PCL).
// PCL is a multiplatform C++ framework for development of PixInsight modules.
//
// Copyright (c) 2003-2024 Pleiades Astrophoto S.L. All Rights Reserved.
//
// Redistribution and use in both source and binary forms, with or without
// modification, is permitted provided that the following conditions are met:
//
// 1. All redistributions of source code must retain the above copyright
//    notice, this list of conditions and the following disclaimer.
//
// 2. All redistributions in binary form must reproduce the above copyright
//    notice, this list of conditions and the following disclaimer in the
//    documentation and/or other materials provided with the distribution.
//
// 3. Neither the names "PixInsight" and "Pleiades Astrophoto", nor the names
//    of their contributors, may be used to endorse or promote products derived
//    from this software without specific prior written permission. For written
//    permission, please contact info@pixinsight.com.
//
// 4. All products derived from this software, in any form whatsoever, must
//    reproduce the following acknowledgment in the end-user documentation
//    and/or other materials provided with the product:
//
//    "This product is based on software from the PixInsight project, developed
//    by Pleiades Astrophoto and its contributors (https://pixinsight.com/)."
//
//    Alternatively, if that is where third-party acknowledgments normally
//    appear, this acknowledgment must be reproduced in the product itself.
//
// THIS SOFTWARE IS PROVIDED BY PLEIADES ASTROPHOTO AND ITS CONTRIBUTORS
// "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
// TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL PLEIADES ASTROPHOTO OR ITS
// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
// EXEMPLARY OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, BUSINESS
// INTERRUPTION; PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; AND LOSS OF USE,
// DATA OR PROFITS) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
// CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
// ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
// POSSIBILITY OF SUCH DAMAGE.
// ----------------------------------------------------------------------------
 
#ifndef __PCL_PSFFit_h
#define __PCL_PSFFit_h
 
 
#include <pcl/Defs.h>
#include <pcl/Diagnostics.h>
 
#include <pcl/ImageVariant.h>
#include <pcl/Matrix.h>
#include <pcl/String.h>
 
namespace pcl
{
 
// ----------------------------------------------------------------------------
 
namespace PSFunction
{
   enum value_type
   {
      Invalid  = -1, // Represents an invalid or unsupported PSF type
      Gaussian = 0,  // Gaussian PSF
      Moffat,        // Moffat PSF with a fitted beta parameter
      MoffatA,       // Moffat PSF with fixed beta=10
      Moffat8,       // Moffat PSF with fixed beta=8
      Moffat6,       // Moffat PSF with fixed beta=6
      Moffat4,       // Moffat PSF with fixed beta=4
      Moffat25,      // Moffat PSF with fixed beta=2.5
      Moffat15,      // Moffat PSF with fixed beta=1.5
      Lorentzian,    // Lorentzian PSF, equivalent to a Moffat PSF with fixed beta=1
      VariableShape, // Variable shape PSF
      NumberOfFunctions,
      Auto = 999,    // Automatic selection for optimization, process-dependent.
      Default = Gaussian
   };
}
 
// ----------------------------------------------------------------------------
 
namespace PSFFitStatus
{
   enum value_type
   {
      Invalid   = -1,
      NotFitted = 0,
      FittedOk,
      BadParameters,
      NoSolution,
      NoConvergence,
      InaccurateSolution,
      UnknownError
   };
}
 
// ----------------------------------------------------------------------------
 
struct PSFData
{
   using psf_function = PSFunction::value_type;
 
   using psf_fit_status = PSFFitStatus::value_type;
 
   using component = DPoint::component;
 
   psf_function   function = PSFunction::Invalid;   
   bool           circular = false;                 
   psf_fit_status status = PSFFitStatus::NotFitted; 
   bool           celestial = false;   
   double         B = 0;      
   double         A = 0;      
   DPoint         b0 = 0.0;   
   DPoint         c0 = 0.0;   
   DPoint         q0 = 0.0;   
   double         sx = 0;     
   double         sy = 0;     
   double         theta = 0;  
   double         beta = 0;   
   double         flux = 0;   
   double         signal = 0; 
   unsigned       signalCount = 0; 
   double         mad = 0;    
   PSFData() = default;
 
   PSFData( const PSFData& ) = default;
 
   PSFData( PSFData&& ) = default;
 
   PSFData& operator =( const PSFData& ) = default;
 
   PSFData& operator =( PSFData&& ) = default;
 
   operator bool() const
   {
      return status == PSFFitStatus::FittedOk;
   }
 
   operator double() const
   {
      return (signal != 0) ? signal : flux;
   }
 
   double operator []( int i ) const noexcept
   {
      return (i == 0) ? c0.x : c0.y;
   }
 
   String FunctionName() const;
 
   String StatusText() const;
 
   double FWHMx() const
   {
      return FWHM( function, sx, beta );
   }
 
   double FWHMy() const
   {
      return FWHM( function, sy, beta );
   }
 
   double FWTMx() const
   {
      return FWTM( function, sx, beta );
   }
 
   double FWTMy() const
   {
      return FWTM( function, sy, beta );
   }
 
   double Volume() const
   {
      return A*Volume( function, sx, sy, beta );
   }
 
   DRect Bounds( double k = 1.0 ) const
   {
      double dx = k*FWHMx()/2;
      double dy = k*FWHMy()/2;
      return DRect( c0.x-dx, c0.y-dy, c0.x+dx, c0.y+dy );
   }
 
   DRect FWTMBounds( double k = 1.0 ) const
   {
      double dx = k*FWTMx()/2;
      double dy = k*FWTMy()/2;
      return DRect( c0.x-dx, c0.y-dy, c0.x+dx, c0.y+dy );
   }
 
   bool HasCelestialCoordinates() const
   {
      return celestial;
   }
 
   void ToImage( Image& ) const;
 
   static double FWHM( psf_function function, double sigma, double beta = 2 )
   {
      PCL_PRECONDITION( beta > 0 )
      switch ( function )
      {
      case PSFunction::Gaussian:      return 2.3548200450309493 * sigma;
      case PSFunction::Moffat:        return 2 * sigma * Sqrt( Pow2( 1/beta ) - 1 );
      case PSFunction::MoffatA:       return 0.5358113941912513 * sigma;
      case PSFunction::Moffat8:       return 0.6016900619596693 * sigma;
      case PSFunction::Moffat6:       return 0.6998915581984769 * sigma;
      case PSFunction::Moffat4:       return 0.8699588840921645 * sigma;
      case PSFunction::Moffat25:      return 1.1305006161394060 * sigma;
      case PSFunction::Moffat15:      return 1.5328418730817597 * sigma;
      case PSFunction::Lorentzian:    return 2 * sigma;
      case PSFunction::VariableShape: return 2 * sigma * Pow( beta*0.6931471805599453, 1/beta );
      default:                        return 0; // ?!
      }
   }
 
   static double FWTM( psf_function function, double sigma, double beta = 2 )
   {
      PCL_PRECONDITION( beta > 0 )
      switch ( function )
      {
      case PSFunction::Gaussian:      return 4.291932052578694 * sigma;
      case PSFunction::Moffat:        return 2 * sigma * Sqrt( Pow( 10.0, 1/beta ) - 1 );
      case PSFunction::MoffatA:       return 1.017694279819175 * sigma;
      case PSFunction::Moffat8:       return 1.155026289161115 * sigma;
      case PSFunction::Moffat6:       return 1.367917055412454 * sigma;
      case PSFunction::Moffat4:       return 1.764402913213331 * sigma;
      case PSFunction::Moffat25:      return 2.459175822514185 * sigma;
      case PSFunction::Moffat15:      return 3.816589489904712 * sigma;
      case PSFunction::Lorentzian:    return 6 * sigma;
      case PSFunction::VariableShape: return 2 * sigma * Pow( beta*2.302585092994045, 1/beta );
      default:                        return 0; // ?!
      }
   }
 
   static double Volume( psf_function function, double sigma_x, double sigma_y, double beta = 2 )
   {
      PCL_PRECONDITION( beta > 0 )
      PCL_PRECONDITION( function != PSFunction::Lorentzian && (function != PSFunction::Moffat || beta > 1) )
      switch ( function )
      {
      case PSFunction::Gaussian: return 6.2831853071795862 * sigma_x*sigma_y;
      case PSFunction::Moffat:   return 3.1415926535897931 * sigma_x*sigma_y/(beta - 1);
      case PSFunction::MoffatA:  return 0.3490658503988659 * sigma_x*sigma_y;
      case PSFunction::Moffat8:  return 0.4487989505128276 * sigma_x*sigma_y;
      case PSFunction::Moffat6:  return 0.6283185307179586 * sigma_x*sigma_y;
      case PSFunction::Moffat4:  return 1.0471975511965976 * sigma_x*sigma_y;
      case PSFunction::Moffat25: return 2.0943951023931953 * sigma_x*sigma_y;
      case PSFunction::Moffat15: return 6.2831853071795862 * sigma_x*sigma_y;
      default:                   return 0; // ?!
      }
   }
};
 
// ----------------------------------------------------------------------------
 
class PSFFit
{
public:
 
   using psf_function = PSFData::psf_function;
 
   using psf_fit_status = PSFData::psf_fit_status;
 
   PSFData psf;
 
   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 );
 
   PSFFit( const PSFFit& ) = default;
 
   PSFFit( PSFFit&& ) = default;
 
   PSFFit& operator =( const PSFFit& x ) = default;
 
   PSFFit& operator =( PSFFit&& x ) = default;
 
   operator bool() const
   {
      return psf;
   }
 
private:
 
   Matrix S; // matrix of sampled pixel data
   Vector P; // vector of function parameters
 
   // The initial local background measured on the sampling region and the
   // maximum allowed relative difference, for stabilization of local
   // background PSF parameters.
   double m_bkg;
   float  m_bkgMaxVar;
 
   // Growing factor in units of the Full Width at Tenth Maximum (FWTM) for
   // extension/contraction of the PSF flux measurement region.
   float  m_growthForFlux = 1.0F;
 
   // Keep track of successive beta values in L-M iterations for stabilization
   // of shape parameters. This guarantees convergence for Moffat functions.
   mutable float m_beta;
 
   Vector GoodnessOfFit( psf_function, bool circular, bool test = false ) const;
 
   friend class PSFFitEngine;
};
 
// ----------------------------------------------------------------------------
 
} // pcl
 
#endif   // __PCL_PSFFit_h
 
// ----------------------------------------------------------------------------
// EOF pcl/PSFFit.h - Released 2024-06-18T15:48:54Z