PCL StarDetection

jmurphy

Well-known member
NormalizeScaleGradient needs a fast way to find all the stars in the reference and target images.
The javascript version uses a slightly modified version of StarDetector.jsh

For each star, I need the approximate size (for example, the star's bounding rectangle, or a radius), and the (x,y) position.

I notice that the SubframeSelector-StarDetector.cpp provides most of what I require. I would just need to modify it to return the star's bounding rectangle.

Is this likely to be the best option, or is there a better way of doing this?
Thanks, John Murphy
 
Hi John,

The StarAlignment tool uses an improved version of our star detection algorithm with optional multithreaded PSF fitting. Allow me some time to extract the relevant source code and I'll share it with you here.
 
Hi John,

The StarAlignment tool uses an improved version of our star detection algorithm with optional multithreaded PSF fitting. Allow me some time to extract the relevant source code and I'll share it with you here.
Fantastic! Thanks. :)
 
Here we go. The Star structure defines the centroid position (in image coordinates), the total flux (with local background subtracted), and a goodness-of-fit estimate for stars where a PSF has been fitted. You can extend it easily as needed:

Star.h
C++:
#ifndef __PCL_Star_h
#define __PCL_Star_h

#include <pcl/Point.h>

namespace pcl
{

// ----------------------------------------------------------------------------

struct Star
{
   typedef DPoint::component  component;

   DPoint pos = 0.0;     // Barycenter position in pixels, image coordinates.
   float  flux = 0.0F;   // Total flux minus local background.
   float  mad = 0.0F;    // Average absolute deviation of measured PSF.

   Star() = default;

   Star( const DPoint& p, float f, float m = 0 )
      : pos( p ), flux( f ), mad( m )
   {
   }

   Star( double x, double y, float f, float m = 0 )
      : pos( x, y ), flux( f ), mad( m )
   {
   }

   Star( const Star& ) = default;
   Star& operator =( const Star& ) = default;

   // Two stars are equal iif they are located at the same coordinates
   bool operator ==( const Star& s ) const
   {
      return pos == s.pos;
   }

   // Stars are sorted by brightness in descending order (brighter stars come first)
   bool operator <( const Star& s ) const
   {
      return flux > s.flux;
   }

   // Coordinate selection - for compatibility with pcl::QuadTree
   component operator []( int i ) const
   {
      return pos[i];
   }
};

// ----------------------------------------------------------------------------

} // pcl

#endif   // __PCL_Star_h


The StarDetector class implements our star detection algorithm. The header file comes first:

StarDetector.h
C++:
#ifndef __PCL_StarDetector_h
#define __PCL_StarDetector_h

#include <pcl/Array.h>
#include <pcl/ImageVariant.h>

#include "Star.h"

namespace pcl
{

// ----------------------------------------------------------------------------

class PCL_CLASS StarDetector
{
public:

   typedef Array<Star>  star_list;

   /*
    * Detected stars, sorted by brightness in descending order.
    */
   star_list S;

   StarDetector( const ImageVariant& img,
                 int   structureLayers = 5,
                 int   noiseLayers = 1,
                 int   hotPixelFilterRadius = 1,
                 int   noiseReductionFilterRadius = 0,
                 int   minStructureSize = 0,
                 float sensitivity = 0.1F,
                 float peakResponse = 0.75F,
                 float maxDistortion = 0.5F,
                 float upperLimit = 1.0F,
                 bool  invert = false,
                 bool  fitPSF = false,
                 float psfTolerance = 1.5,
                 const UInt8Image* mask = nullptr );

   virtual ~StarDetector()
   {
   }

   static void GetStructureMap( ImageVariant& image, int structureLayers = 5, int noiseLayers = 1,
                                int hotPixelFilterRadius = 1, bool invert = false, const UInt8Image* mask = nullptr );

   static void GetStructures( ImageVariant& image, int structureLayers = 5, int noiseLayers = 1,
                              int hotPixelFilterRadius = 1, bool invert = false, const UInt8Image* mask = nullptr );
};

// ----------------------------------------------------------------------------

} // pcl

#endif   // __PCL_StarDetector_h


And now the implementation of StarDetector:

StarDetector.cpp
C++:
#include "StarDetector.h"

#include <pcl/ATrousWaveletTransform.h>
#include <pcl/Console.h>
#include <pcl/Convolution.h>
#include <pcl/GaussianFilter.h>
#include <pcl/MetaModule.h>
#include <pcl/MorphologicalTransformation.h>
#include <pcl/PSFFit.h>
#include <pcl/Selection.h>
#include <pcl/SeparableConvolution.h>

#include <float.h>

namespace pcl
{

// ----------------------------------------------------------------------------

/*
 * The local background is evaluated for each star on an inflated rectangular
 * region around the star detection structure. bkgDelta is the inflation
 * distance in pixels.
 */
#define BKG_DELTA    3

/*
 * Stretch factor for the barycenter search algorithm, in sigma units.
 * Increase it to make the algorithm more robust to nearby structures, such as
 * multiple/crowded stars and small nebular features. However, too large of a
 * stretch factor will make the algorithm less accurate.
 */
#define XY_STRETCH   1.5

/*
 * Computational cost of the structure detection algorithm.
 */
#define SD_COST      11

// ----------------------------------------------------------------------------

static void HotPixelFilter( Image& image, int hotPixelFilterRadius )
{
   if ( hotPixelFilterRadius > 0 ) // N
   {
      if ( hotPixelFilterRadius > 1 )
         MorphologicalTransformation( MedianFilter(), CircularStructure( 2*hotPixelFilterRadius + 1 ) ) >> image;
      else
         MorphologicalTransformation( MedianFilter(), BoxStructure( 3 ) ) >> image;
   }
   else
      image.Status() += image.NumberOfPixels();
}

// ----------------------------------------------------------------------------

static void NoiseLayersFilter( Image& image, int noiseLayers )
{
   if ( noiseLayers > 0 ) // N
   {
      GaussianFilter G( 1 + (1 << noiseLayers) );
      if ( G.Size() < SeparableConvolution::FasterThanNonseparableFilterSize( Thread::NumberOfThreads( PCL_MAX_PROCESSORS ) ) )
         Convolution( G ) >> image;
      else
         SeparableConvolution( G.AsSeparableFilter() ) >> image;
   }
   else
      image.Status() += image.NumberOfPixels();
}

// ----------------------------------------------------------------------------

static void GetStructureMap( Image& map, int structureLayers, int noiseLayers, int hotPixelFilterRadius,
                             const UInt8Image* mask = nullptr )
{
   bool initializeStatus = map.Status().IsInitializationEnabled();
   if ( initializeStatus )
   {
      map.Status().Initialize( "Structure map", map.NumberOfPixels()*SD_COST );
      map.Status().DisableInitialization();
   }

   /*
    * Optional hot pixel removal with a median filter.
    */
   HotPixelFilter( map, hotPixelFilterRadius ); // N

   /*
    * Optional noise reduction with a low-pass filter.
    */
   NoiseLayersFilter( map, noiseLayers ); // N

   /*
    * Flatten the image with a high-pass filter.
    */
   {
      Image s( map );
      GaussianFilter G( 1 + (1 << structureLayers) );
      if ( G.Size() < SeparableConvolution::FasterThanNonseparableFilterSize( Thread::NumberOfThreads( PCL_MAX_PROCESSORS ) ) )
         Convolution( G ) >> s;  // N
      else
         SeparableConvolution( G.AsSeparableFilter() ) >> s;
      map.Status() = s.Status();
      map.Sub( s );     // N
      map.Truncate();   // N
      map.Rescale();    // N
   }

   /*
    * Strength the smallest structures with a dilation filter.
    */
   MorphologicalTransformation( DilationFilter(), BoxStructure( 3 ) ) >> map; // N

   /*
    * Adaptive binarization based on noise evaluation.
    */
   double median = map.Median();
   if ( 1 + median == 1 )
   {
      // An empty image, or a noise-free synthetic star field.
      map.SetRangeClipping( 0, 1 );
      median = map.Median();
      map.Binarize( median + map.MAD( median ) );  // N
      map.DisableRangeClipping();
   }
   else
   {
      // A "natural image" or something similar.
      static const float B3[] =
      {
         0.003906F, 0.015625F, 0.023438F, 0.015625F, 0.003906F,
         0.015625F, 0.062500F, 0.093750F, 0.062500F, 0.015625F,
         0.023438F, 0.093750F, 0.140625F, 0.093750F, 0.023438F,
         0.015625F, 0.062500F, 0.093750F, 0.062500F, 0.015625F,
         0.003906F, 0.015625F, 0.023438F, 0.015625F, 0.003906F
      };
      static const float B3v[] = { 0.0625F, 0.25F, 0.375F, 0.25F, 0.0625F };
      static const float B3k[] = { 0.8907F, 0.2007F, 0.0856F };

      ATrousWaveletTransform::WaveletScalingFunction H;
      if ( SeparableConvolution::FasterThanNonseparableFilterSize( Thread::NumberOfThreads( PCL_MAX_PROCESSORS ) ) > 5 )
         H.Set( KernelFilter( B3, 5 ) );
      else
         H.Set( SeparableFilter( B3v, B3v, 5 ) );
      ATrousWaveletTransform W( H, 2 );
      W << map; // 2*N
      float noise = W.NoiseKSigma( 1/*j*/, 3/*k*/, 0.01F/*eps*/, 10/*nit*/ )/B3k[1];
      map.Binarize( median + 3*noise ); // N
   }

   /*
    * Optional detection mask.
    */
   if ( mask != nullptr )
      map.Mul( *mask );          // N
   else
      map.Status() += map.NumberOfPixels();

   if ( initializeStatus )
      map.Status().EnableInitialization();
}

// ----------------------------------------------------------------------------

static void GetStructures( Image& image, int structureLayers, int noiseLayers, int hotPixelFilterRadius, const UInt8Image* mask )
{
   bool initializeStatus = image.Status().IsInitializationEnabled();
   if ( initializeStatus )
   {
      image.Status().Initialize( "Extracting structures", image.NumberOfPixels()*(SD_COST + 1) );
      image.Status().DisableInitialization();
   }

   Image map( image );
   GetStructureMap( map, structureLayers, noiseLayers, hotPixelFilterRadius, mask );
   image.Status() = map.Status();
   image.Mul( map );               // N

   if ( initializeStatus )
      image.Status().EnableInitialization();
}

// ----------------------------------------------------------------------------

struct StarData
{
   DPoint pos;  // barycenter image coordinates
   Rect   rect;
   float  bkg;  // local background
   float  flux; // total flux
   float  peak; // peak value
   int    size; // structure size in square pixels
};

static void
GetStarParameters( StarData& data, const Image& image, const Rect& rect, const Array<Point>& sp )
{
   // Calculate a mean local background as the median of four selected regions
   // of background pixels.
   {
      Rect r = rect.InflatedBy( BKG_DELTA );
      FMatrix bt = FMatrix::FromImage( image, Rect(    r.x0,    r.y0,    r.x1, rect.y0 ) );
      FMatrix bl = FMatrix::FromImage( image, Rect(    r.x0, rect.y0, rect.x0, rect.y1 ) );
      FMatrix bb = FMatrix::FromImage( image, Rect(    r.x0, rect.y1,    r.x1,    r.y1 ) );
      FMatrix br = FMatrix::FromImage( image, Rect( rect.x1, rect.y0,    r.x1, rect.y1 ) );
      Array<float> b;
      b.Add( bt.Begin(), bt.End() );
      b.Add( bl.Begin(), bl.End() );
      b.Add( bb.Begin(), bb.End() );
      b.Add( br.Begin(), br.End() );
      data.rect = r;
      data.bkg = Median( b.Begin(), b.End() );
   }

   // Compute barycenter coordinates.
   {
      FMatrix Z = FMatrix::FromImage( image, rect );
      for ( float& z : Z )
         if ( z < 1.0e-07 ) // zero out insignificant values
            z = 0;
      Z.Truncate( Range( FMatrix( Z ).Median() + XY_STRETCH*Z.StdDev(), 0.0, 1.0 ), 1.0 );
      Z.Rescale();
      const float* z = Z.Begin();
      double sx = 0, sy = 0, sz = 0;
      for ( int y = rect.y0; y < rect.y1; ++y )
         for ( int x = rect.x0; x < rect.x1; ++x, ++z )
            if ( *z > 0 )
            {
               sx += *z * x;
               sy += *z * y;
               sz += *z;
            }
      if ( 1 + sz == 1 )
      {
         data.flux = 0;
         return;
      }
      data.pos.x = sx/sz;
      data.pos.y = sy/sz;
   }

   // Total flux, peak value, and structure size.
   data.flux = data.peak = 0;
   for ( Array<Point>::const_iterator i = sp.Begin(); i != sp.End(); ++i )
   {
      float f = image( *i );
      data.flux += f;
      if ( f > data.peak )
         data.peak = f;
   }
   data.size = int( sp.Length() );
};

// ----------------------------------------------------------------------------

struct PSFFitData
{
   DPoint pos;
   Rect   rect;
   float  flux;
};

class PSFFitThread : public Thread
{
public:

   StarDetector::star_list stars;

   PSFFitThread( const AbstractImage::ThreadData& data,
                 Image& image, float tolerance, const Array<PSFFitData>& psfData,
                 int start, int end )
      : m_data( data )
      , m_image( image )
      , m_tolerance( tolerance )
      , m_psfData( psfData )
      , m_start( start )
      , m_end( end )
   {
   }

   PCL_HOT_FUNCTION void Run() override
   {
      INIT_THREAD_MONITOR()

      for ( int i = m_start; i < m_end; ++i )
      {
         const PSFFitData& d = m_psfData[i];

         // Adaptive PSF sampling region
         int size = Max( d.rect.Width(), d.rect.Height() );
         Rect rect = DRect( d.pos + (0.5-size), d.pos + (0.5+size) ).TruncatedToInt();
         for ( double m0 = 1; ; )
         {
            double m = FMatrix::FromImage( m_image, rect ).Median();
            if ( m0 <= m || (m0 - m)/m0 < 0.01 )
               break;
            m0 = m;
            rect.InflateBy( 1, 1 );
         }

         PSFFit fit( ImageVariant( &m_image ), d.pos + 0.5, rect, PSFunction::Gaussian, true/*circular*/ );
         if ( fit )
            if ( DRect( rect ).DeflatedBy( rect.Width()*0.15 ).Includes( fit.psf.c0 ) )
               if ( fit.psf.c0.DistanceTo( d.pos + 0.5 ) < m_tolerance )
                  stars << Star( fit.psf.c0 - 0.5, d.flux/*fit.psf.A*/, fit.psf.mad ); // flux is more robust than amplitude here

         UPDATE_THREAD_MONITOR( 16 )
      }
   }

private:

   const AbstractImage::ThreadData& m_data;
         Image&                     m_image;
         float                      m_tolerance;
   const Array<PSFFitData>&         m_psfData;
         int                        m_start, m_end;
};

static void GetStars( StarDetector::star_list& S, Image& image,
                      int structureLayers, int noiseLayers,
                      int hotPixelFilterRadius, int noiseReductionFilterRadius, int minStructureSize,
                      float sensitivity, float peakResponse, float maxDistortion,
                      float upperLimit, bool invert, bool fitPSF, float psfTolerance,
                      const UInt8Image* mask )
{
   /*
    * If the invert flag is set, then we are looking for dark structures on a
    * bright background.
    */
   if ( invert )
      image.Invert();

   bool initializeStatus = image.Status().IsInitializationEnabled();

   /*
    * Optional noise reduction
    */
   if ( noiseReductionFilterRadius > 0 )
   {
      if ( initializeStatus )
      {
         image.Status().Initialize( "Noise reduction", 2*image.NumberOfPixels() );
         image.Status().DisableInitialization();
      }

      // If noise reduction is enabled, remove hot pixels first - Otherwise
      // they would be "promoted to stars".
      HotPixelFilter( image, hotPixelFilterRadius ); // N
      hotPixelFilterRadius = 0;

      GaussianFilter G( (noiseReductionFilterRadius << 1)|1 );
      if ( G.Size() < SeparableConvolution::FasterThanNonseparableFilterSize( Thread::NumberOfThreads( PCL_MAX_PROCESSORS ) ) )
         Convolution( G ) >> image;  // N
      else
         SeparableConvolution( G.AsSeparableFilter() ) >> image;

      if ( initializeStatus )
         image.Status().EnableInitialization();
   }

   /*
    * Structure map
    */
   Image map( image );
   GetStructureMap( map, structureLayers, noiseLayers, hotPixelFilterRadius, nullptr/*mask*/ );

   if ( initializeStatus )
   {
      image.Status().Initialize( "Detecting stars", size_type( image.Width()-1 )*size_type( image.Height()-1 ) );
      image.Status().DisableInitialization();
   }

   Array<PSFFitData> psfData;

   /*
    * Star detection
    */
   for ( int y0 = 0, x1 = image.Width()-1, y1 = image.Height()-1; y0 < y1; ++y0 )
      for ( int x0 = 0; x0 < x1; ++x0, ++image.Status() )
      {
         // Exclude background pixels and already visited pixels
         if ( map( x0, y0 ) == 0 )
            continue;

         // Star bounding rectangle
         Rect r( x0, y0, x0+1, y0+1 );

         // Current star points
         Array<Point> sp;

         // Grow star region downward
         for ( int y = y0, x = x0, xa, xb; ; )
         {
            // Add this pixel to the current star
            sp << Point( x, y );

            // Explore left segment of this row
            for ( xa = x; xa > 0; )
            {
               if ( map( xa-1, y ) == 0 )
                  break;
               --xa;
               sp << Point( xa, y );
            }

            // Explore right segment of this row
            for ( xb = x; xb < x1; )
            {
               if ( map( xb+1, y ) == 0 )
                  break;
               ++xb;
               sp << Point( xb, y );
            }

            // xa and xb are now the left and right boundary limits,
            // respectively, of this row in the current star.

            if ( xa < r.x0 )  // update left boundary
               r.x0 = xa;

            if ( xb >= r.x1 ) // update right boundary
               r.x1 = xb + 1; // bottom-right corner excluded (PCL-specific)

            // Prepare for next row
            ++y;

            // Decide whether we are done with this star now, or if
            // there is at least one more row that must be explored.

            bool nextRow = false;

            // We explore the next row from left to right. We'll continue
            // gathering pixels if we find at least one nonzero pixel.
            for ( x = xa; x <= xb; ++x )
               if ( map( x, y ) != 0 )
               {
                  nextRow = true;
                  break;
               }

            if ( !nextRow )
               break;

            // Update bottom boundary
            r.y1 = y + 1;  // Rectangle<int> *excludes* the bottom-right corner

            // Terminate if we reach the last row of the image
            if ( y == y1 )
               break;
         }

         /*
          * If this is a reliable star, compute its barycenter coordinates
          * and add it to the star list.
          *
          * Rejection criteria:
          *
          * * If this structure is touching a border of the image, reject
          *   it. We cannot compute an accurate position for a clipped star.
          *
          * * Too small structures are rejected. This mainly prevents
          *   inclusion of hot (or cold) pixels. This condition is enforced
          *   by the hot pixel removal and noise reduction steps performed
          *   during the structure detection phase.
          *
          * * Too large structures are rejected. This prevents inclusion of
          *   extended nonstellar objects and saturated bright stars. This
          *   is also part of the structure detection algorithm.
          *
          * * Too elongated stars are rejected. The maxDistortion parameter
          *   determines the maximum distortion allowed. A perfect square
          *   has distortion = 1. The distortion of a perfect circle is
          *   pi/4, or about 0.8.
          *
          * * We don't trust stars whose centroids are too misplaced with
          *   respect to their peak positions. This prevents detection of
          *   multiple stars, where an accurate position cannot be computed.
          *
          * * Too flat structures are rejected. The peakResponse parameter
          *   defines a peakedness threshold necessary for a structure to be
          *   idenfified as a valid star.
          */
         if ( r.Width() > 1 && r.Height() > 1 )
            if ( r.y0 > 0 && r.y1 <= y1 && r.x0 > 0 && r.x1 <= x1 )
               if ( sp.Length() >= size_type( minStructureSize ) )
               {
                  int d = Max( r.Width(), r.Height() );
                  if ( float( sp.Length() )/d/d > maxDistortion )
                  {
                     StarData data;
                     GetStarParameters( data, image, r, sp );
                     if ( data.flux > 0 )
                        if ( data.peak <= upperLimit )
                        {
                           Point p = (data.pos + 0.5).TruncatedToInt();
                           if ( mask == nullptr || (*mask)( p ) != 0 )
                           {
                              // Detection level corrected for peak response
                              float f = data.peak - (1 - peakResponse)*data.flux/data.size;
                              if ( data.bkg == 0 || (f - data.bkg)/data.bkg > sensitivity )
                                 if ( image( p ) > 0.85*data.peak )
                                    if ( image.Median( r ) < peakResponse*data.peak )
                                       if ( fitPSF )
                                          psfData << PSFFitData{ data.pos, data.rect, data.flux };
                                       else
                                          S << Star( data.pos, data.flux );
                           }
                        }
                  }
               }

         // Erase this structure
         for ( const Point& p : sp )
            map( p ) = 0;
      }

   /*
    * PSF fitting
    */
   if ( !psfData.IsEmpty() )
   {
      if ( initializeStatus )
      {
         image.Status().EnableInitialization();
         image.Status().Initialize( String().Format( "Fitting %u stars", psfData.Length() ), psfData.Length() );
         image.Status().DisableInitialization();
      }

      Array<size_type> L = Thread::OptimalThreadLoads( psfData.Length() );
      ReferenceArray<PSFFitThread> threads;
      AbstractImage::ThreadData data( image, psfData.Length() );
      for ( int i = 0, n = 0; i < int( L.Length() ); n += int( L[i++] ) )
         threads.Add( new PSFFitThread( data, image, psfTolerance, psfData, n, n + int( L[i] ) ) );
      AbstractImage::RunThreads( threads, data );

      for ( const PSFFitThread& thread : threads )
         S << thread.stars;

      threads.Destroy();
   }

   if ( initializeStatus )
      image.Status().EnableInitialization();
}

// ----------------------------------------------------------------------------

template <class P>
static Image GetIntensity( const GenericImage<P>& image )
{
   StatusMonitor status = image.Status();
   image.Status() = StatusMonitor();
   Image I;
   image.GetIntensity( I );
   image.Status() = I.Status() = status;
   return I;
}

// ----------------------------------------------------------------------------

template <class P>
static void GetStars( StarDetector::star_list& S, const GenericImage<P>& image,
                      int structureLayers, int noiseLayers,
                      int hotPixelFilterRadius, int noiseReductionFilterRadius, int minStructureSize,
                      float sensitivity, float peakResponse, float maxDistortion,
                      float upperLimit, bool invert, bool fitPSF, float psfTolerance,
                      const UInt8Image* mask )
{
   image.PushSelections();
   image.ResetSelection();

   Image I = GetIntensity( image );
   GetStars( S, I, structureLayers, noiseLayers,
             hotPixelFilterRadius, noiseReductionFilterRadius, minStructureSize,
             sensitivity, peakResponse, maxDistortion, upperLimit, invert, fitPSF, psfTolerance, mask );

   image.PopSelections();
}

// ----------------------------------------------------------------------------

StarDetector::StarDetector( const ImageVariant& image, int structureLayers, int noiseLayers,
                            int hotPixelFilterRadius, int noiseReductionFilterRadius, int minStructureSize,
                            float sensitivity, float peakResponse, float maxDistortion,
                            float upperLimit, bool invert, bool fitPSF, float psfTolerance,
                            const UInt8Image* mask )
{
   if ( mask != nullptr )
      if ( mask->IsEmpty() )
         mask = nullptr;

   if ( !image.IsComplexSample() )
      if ( image.IsFloatSample() )
         switch ( image.BitsPerSample() )
         {
         case 32:
            GetStars( S, static_cast<const Image&>( *image ), structureLayers, noiseLayers,
                      hotPixelFilterRadius, noiseReductionFilterRadius, minStructureSize,
                      sensitivity, peakResponse, maxDistortion, upperLimit, invert, fitPSF, psfTolerance, mask );
            break;
         case 64:
            GetStars( S, static_cast<const DImage&>( *image ), structureLayers, noiseLayers,
                      hotPixelFilterRadius, noiseReductionFilterRadius, minStructureSize,
                      sensitivity, peakResponse, maxDistortion, upperLimit, invert, fitPSF, psfTolerance, mask );
            break;
         }
      else
         switch ( image.BitsPerSample() )
         {
         case 8:
            GetStars( S, static_cast<const UInt8Image&>( *image ), structureLayers, noiseLayers,
                      hotPixelFilterRadius, noiseReductionFilterRadius, minStructureSize,
                      sensitivity, peakResponse, maxDistortion, upperLimit, invert, fitPSF, psfTolerance, mask );
            break;
         case 16:
            GetStars( S, static_cast<const UInt16Image&>( *image ), structureLayers, noiseLayers,
                      hotPixelFilterRadius, noiseReductionFilterRadius, minStructureSize,
                      sensitivity, peakResponse, maxDistortion, upperLimit, invert, fitPSF, psfTolerance, mask );
            break;
         case 32:
            GetStars( S, static_cast<const UInt32Image&>( *image ), structureLayers, noiseLayers,
                      hotPixelFilterRadius, noiseReductionFilterRadius, minStructureSize,
                      sensitivity, peakResponse, maxDistortion, upperLimit, invert, fitPSF, psfTolerance, mask );
            break;
         }

   // Sort stars in decreasing brightness order
   S.Sort();
}

// ----------------------------------------------------------------------------

static void GetStructureMapImp( Image& image, int structureLayers, int noiseLayers,
                                int hotPixelFilterRadius, bool invert,
                                const UInt8Image* mask )
{
   if ( image.IsColor() )
   {
      Image I = GetIntensity( image );
      if ( invert )
         I.Invert();
      GetStructureMap( I, structureLayers, noiseLayers, hotPixelFilterRadius, mask );
      image = I;
   }
   else
   {
      if ( invert )
         image.Invert();
      GetStructureMap( image, structureLayers, noiseLayers, hotPixelFilterRadius, mask );
   }
}

template <class P>
static void GetStructureMapImp( GenericImage<P>& image, int structureLayers, int noiseLayers,
                                int hotPixelFilterRadius, bool invert,
                                const UInt8Image* mask )
{
   if ( !image.IsColor() )
      image.SelectChannel( 0 );
   Image I = GetIntensity( image );
   if ( invert )
      I.Invert();
   GetStructureMap( I, structureLayers, noiseLayers, hotPixelFilterRadius, mask );
   image = I;
}

void StarDetector::GetStructureMap( ImageVariant& image, int structureLayers, int noiseLayers,
                                    int hotPixelFilterRadius, bool invert,
                                    const UInt8Image* mask )
{
   if ( mask != nullptr )
      if ( mask->IsEmpty() )
         mask = nullptr;

   if ( !image.IsComplexSample() )
      if ( image.IsFloatSample() )
         switch ( image.BitsPerSample() )
         {
         case 32:
            GetStructureMapImp( static_cast<Image&>( *image ), structureLayers, noiseLayers,
                                hotPixelFilterRadius, invert, mask );
            break;
         case 64:
            GetStructureMapImp( static_cast<DImage&>( *image ), structureLayers, noiseLayers,
                                hotPixelFilterRadius, invert, mask );
            break;
         }
      else
         switch ( image.BitsPerSample() )
         {
         case 8:
            GetStructureMapImp( static_cast<UInt8Image&>( *image ), structureLayers, noiseLayers,
                                hotPixelFilterRadius, invert, mask );
            break;
         case 16:
            GetStructureMapImp( static_cast<UInt16Image&>( *image ), structureLayers, noiseLayers,
                                hotPixelFilterRadius, invert, mask );
            break;
         case 32:
            GetStructureMapImp( static_cast<UInt32Image&>( *image ), structureLayers, noiseLayers,
                                hotPixelFilterRadius, invert, mask );
            break;
         }
}

// ----------------------------------------------------------------------------

static void GetStructuresImp( Image& image, int structureLayers, int noiseLayers,
                              int hotPixelFilterRadius, bool invert,
                              const UInt8Image* mask )
{
   if ( image.IsColor() )
   {
      Image I = GetIntensity( image );
      if ( invert )
         I.Invert();
      GetStructures( I, structureLayers, noiseLayers, hotPixelFilterRadius, mask );
      image = I;
   }
   else
   {
      if ( invert )
         image.Invert();
      GetStructures( image, structureLayers, noiseLayers, hotPixelFilterRadius, mask );
   }
}

template <class P>
static void GetStructuresImp( GenericImage<P>& image, int structureLayers, int noiseLayers,
                              int hotPixelFilterRadius, bool invert,
                              const UInt8Image* mask )
{
   if ( !image.IsColor() )
      image.SelectChannel( 0 );
   Image I = GetIntensity( image );
   I.Status() = image.Status();
   if ( invert )
      I.Invert();
   GetStructures( I, structureLayers, noiseLayers, hotPixelFilterRadius, mask );
   image = I;
}

void StarDetector::GetStructures( ImageVariant& image, int structureLayers, int noiseLayers,
                                  int hotPixelFilterRadius, bool invert,
                                  const UInt8Image* mask )
{
   if ( mask != nullptr )
      if ( mask->IsEmpty() )
         mask = nullptr;

   if ( !image.IsComplexSample() )
      if ( image.IsFloatSample() )
         switch ( image.BitsPerSample() )
         {
         case 32:
               GetStructuresImp( static_cast<Image&>( *image ), structureLayers, noiseLayers,
                                 hotPixelFilterRadius, invert, mask );
               break;
         case 64:
            GetStructuresImp( static_cast<DImage&>( *image ), structureLayers, noiseLayers,
                              hotPixelFilterRadius, invert, mask );
            break;
         }
      else
         switch ( image.BitsPerSample() )
         {
         case 8:
            GetStructuresImp( static_cast<UInt8Image&>( *image ), structureLayers, noiseLayers,
                              hotPixelFilterRadius, invert, mask );
            break;
         case 16:
            GetStructuresImp( static_cast<UInt16Image&>( *image ), structureLayers, noiseLayers,
                              hotPixelFilterRadius, invert, mask );
            break;
         case 32:
            GetStructuresImp( static_cast<UInt32Image&>( *image ), structureLayers, noiseLayers,
                              hotPixelFilterRadius, invert, mask );
            break;
         }
}

// ----------------------------------------------------------------------------

} // pcl

Note that the detection routine has been designed for increased robustness of star positions, which is an important requirement for image registration. See the cascade of rejection conditions in lines 536 to 560 of StarDetector.cpp. For your specific application, maybe you'll prefer to relax some of these conditions, or to add additional ones better tailored to your requirements.

As for including the detection rectangle in the Star structure, you can get rectangle coordinates in two ways:

- If you don't fit PSFs, the 'r' variable (declared at line 439 of StarDetector.cpp) contains the bounding rectangle for each detected star.

- If you fit PSFs, the 'rect' variable (declared at line 343 of StarDetector.cpp) is the PSF sampling region for each star, computed adaptively by a robust median stabilization algorithm.

I think the code is sufficiently self-documented, but of course don't hesitate to ask if you have doubts or need further clarification.
 
Here we go. The Star structure defines the centroid position (in image coordinates), the total flux (with local background subtracted), and a goodness-of-fit estimate for stars where a PSF has been fitted. You can extend it easily as needed:
...
Note that the detection routine has been designed for increased robustness of star positions, which is an important requirement for image registration. See the cascade of rejection conditions in lines 536 to 560 of StarDetector.cpp. For your specific application, maybe you'll prefer to relax some of these conditions, or to add additional ones better tailored to your requirements.

As for including the detection rectangle in the Star structure, you can get rectangle coordinates in two ways:

- If you don't fit PSFs, the 'r' variable (declared at line 439 of StarDetector.cpp) contains the bounding rectangle for each detected star.
- If you fit PSFs, the 'rect' variable (declared at line 343 of StarDetector.cpp) is the PSF sampling region for each star, computed adaptively by a robust median stabilization algorithm.

I think the code is sufficiently self-documented, but of course don't hesitate to ask if you have doubts or need further clarification.
I have successfully added the code, and am now detecting stars :)
It works great! Thanks!
 
Back
Top