sourceFinder.hpp

Go to the documentation of this file.

/** \file sourceFinder.hpp
 * \brief Declares and defines a class for finding sources in images
 * \ingroup image_processing_files
 * \author Jared R. Males (jaredmales@gmail.com)
 *
 */
 
//***********************************************************************//
// Copyright 2021 Jared R. Males (jaredmales@gmail.com)
//
// This file is part of mxlib.
//
// mxlib is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation, either version 3 of the License, or
// (at your option) any later version.
//
// mxlib is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License
// along with mxlib.  If not, see <http://www.gnu.org/licenses/>.
//***********************************************************************//
 
#ifndef improc_sourceFinder_hpp
#define improc_sourceFinder_hpp
 
#include "eigenImage.hpp"
 
#include "../math/vectorUtils.hpp"
 
namespace mx
{
namespace improc
{
 
/// Find stars in an image by SNR thresholding
/**
 * \ingroup image_utils
 */
template <typename _realT>
class sourceFinder
{
 
  public:
    typedef _realT realT; ///< The real valued type in which to do calculations
 
    typedef std::pair<int, int> pixelT; ///< Pixel coordinates
 
  protected:
    realT m_threshold{ 5 }; ///< The SNR threshold to use
 
    realT m_minsep{ 2 }; ///< The minimum separation between stars.  Closer SNR detections are treated as the same star.
 
    bool m_useMedian{
        true }; ///< If true, then the median is used for standard deviation calculation. The mean is used otherwise.
 
    // Working memory, retained to help with repeated calls.
    std::vector<float> m_vals;
    std::vector<pixelT> m_pixs;
    std::vector<size_t> m_snr;
 
  public:
    /// Default c'tor
    sourceFinder();
 
    /// Constructor to set up the algorithm
    sourceFinder( realT thresh, ///< [in]
                  realT minsep  ///< [in]
    );
 
    /// Clear the working memory
    void clear();
 
    /// Set the SNR threshold
    void treshold( const realT &thresh /**< [in] the new threshold*/ );
 
    /// Get the SNR threshold
    /** \returns the current value of m_threshold;
     */
    realT threshold();
 
    /// Set the minimum separation
    void minsep( const realT &msep /**< [in] the new minimum separation*/ );
 
    /// Get the minimum separation between stars
    /** \returns the current value of m_minsep;
     */
    realT minsep();
 
    /// Set the useMedian flag
    void useMedian( const bool &useMed /**< [in] the new useMedian flag value*/ );
 
    /// Get the useMedian flag
    /** \returns the current value of m_useMedian;
     */
    bool useMedian();
 
    /// Find stars in an image
    /** Returns the coordinates as (row, col) of the highest SNR pixel in each unique star.
     *
     * \returns a vector of pixel coordinates
     */
    std::vector<pixelT> operator()( const eigenImage<realT> &im /**< [in] the image to search for stars*/ );
};
 
template <typename realT>
sourceFinder<realT>::sourceFinder()
{
}
 
template <typename realT>
sourceFinder<realT>::sourceFinder( realT thresh, realT minsep ) : m_threshold( thresh ), m_minsep( minsep )
{
}
 
template <typename realT>
void sourceFinder<realT>::treshold( const realT &thresh )
{
    m_threshold = thresh;
}
 
template <typename realT>
void sourceFinder<realT>::clear()
{
    m_vals.clear();
    m_pixs.clear();
    m_snr.clear();
}
 
template <typename realT>
realT sourceFinder<realT>::threshold()
{
    return m_threshold;
}
 
template <typename realT>
void sourceFinder<realT>::minsep( const realT &msep )
{
    m_minsep = msep;
}
 
template <typename realT>
realT sourceFinder<realT>::minsep()
{
    return m_minsep;
}
 
template <typename realT>
void sourceFinder<realT>::useMedian( const bool &useMed )
{
    m_useMedian = useMed;
}
 
template <typename realT>
bool sourceFinder<realT>::useMedian()
{
    return m_useMedian;
}
 
template <typename realT>
std::vector<typename sourceFinder<realT>::pixelT> sourceFinder<realT>::operator()( const eigenImage<realT> &im )
{
    m_vals.resize( im.rows() * im.cols() );
    m_pixs.resize( im.rows() * im.cols() );
 
    int p = 0;
    for( int cc = 0; cc < im.cols(); ++cc )
    {
        for( int rr = 0; rr < im.rows(); ++rr )
        {
            m_vals[p] = im( rr, cc );
            m_pixs[p] = { rr, cc };
            ++p;
        }
    }
 
    realT mm;
 
    if( m_useMedian )
        mm = math::vectorMedian( m_vals );
    else
        mm = math::vectorMean( m_vals );
    realT std = sqrt( math::vectorVariance( m_vals, mm ) );
 
    // Get the SNR for any pixels above the threshold
    m_snr.clear();
    for( size_t n = 0; n < m_vals.size(); ++n )
    {
        if( ( m_vals[n] - mm ) / std >= m_threshold )
        {
            m_snr.push_back( n );
        }
    }
 
    // Delete the pixels above threshold which are with minsep of a higher pixel
    realT mms2 = m_minsep * m_minsep; // save an op
    for( long n = 0; n < m_snr.size(); ++n )
    {
        for( long m = 0; m < m_snr.size(); ++m )
        {
            if( m == n )
                continue;
 
            // Squared distance:
            realT d2 = pow( m_pixs[m_snr[n]].first - m_pixs[m_snr[m]].first, 2 ) +
                       pow( m_pixs[m_snr[n]].second - m_pixs[m_snr[m]].second, 2 );
 
            if( d2 < mms2 )
            {
                if( m_vals[m_snr[n]] < m_vals[m_snr[m]] ) // If this is the higher pixel, delete the lower
                {
                    m_snr.erase( m_snr.begin() + n );
                    --n;
                    break; // So we're done comparing to the one we just deleted.
                }
                else // This is lower, so delete it.
                {
                    m_snr.erase( m_snr.begin() + m );
                    if( m < n )
                        --n; // changing size behind
                    --m;
                    continue;
                }
            }
        }
    }
 
    std::vector<pixelT> found( m_snr.size() );
 
    for( size_t n = 0; n < m_snr.size(); ++n )
    {
        found[n] = m_pixs[m_snr[n]];
    }
 
    return found;
}
 
} // namespace improc
} // namespace mx
 
#endif // improc_sourceFinder_hpp