Line data Source code
1 : /** \file sourceFinder.hpp
2 : * \brief Declares and defines a class for finding sources in images
3 : * \ingroup image_processing_files
4 : * \author Jared R. Males (jaredmales@gmail.com)
5 : *
6 : */
7 :
8 : //***********************************************************************//
9 : // Copyright 2021 Jared R. Males (jaredmales@gmail.com)
10 : //
11 : // This file is part of mxlib.
12 : //
13 : // mxlib is free software: you can redistribute it and/or modify
14 : // it under the terms of the GNU General Public License as published by
15 : // the Free Software Foundation, either version 3 of the License, or
16 : // (at your option) any later version.
17 : //
18 : // mxlib is distributed in the hope that it will be useful,
19 : // but WITHOUT ANY WARRANTY; without even the implied warranty of
20 : // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
21 : // GNU General Public License for more details.
22 : //
23 : // You should have received a copy of the GNU General Public License
24 : // along with mxlib. If not, see <http://www.gnu.org/licenses/>.
25 : //***********************************************************************//
26 :
27 : #ifndef improc_sourceFinder_hpp
28 : #define improc_sourceFinder_hpp
29 :
30 : #include "eigenImage.hpp"
31 :
32 : #include "../math/vectorUtils.hpp"
33 :
34 : namespace mx
35 : {
36 : namespace improc
37 : {
38 :
39 : /// Find stars in an image by SNR thresholding
40 : /**
41 : * \ingroup image_utils
42 : */
43 : template <typename _realT>
44 : class sourceFinder
45 : {
46 :
47 : public:
48 : typedef _realT realT; ///< The real valued type in which to do calculations
49 :
50 : typedef std::pair<int, int> pixelT; ///< Pixel coordinates
51 :
52 : protected:
53 : realT m_threshold{ 5 }; ///< The SNR threshold to use
54 :
55 : realT m_minsep{ 2 }; ///< The minimum separation between stars. Closer SNR detections are treated as the same star.
56 :
57 : bool m_useMedian{
58 : true }; ///< If true, then the median is used for standard deviation calculation. The mean is used otherwise.
59 :
60 : // Working memory, retained to help with repeated calls.
61 : std::vector<float> m_vals;
62 : std::vector<pixelT> m_pixs;
63 : std::vector<size_t> m_snr;
64 :
65 : public:
66 : /// Default c'tor
67 : sourceFinder();
68 :
69 : /// Constructor to set up the algorithm
70 : sourceFinder( realT thresh, ///< [in]
71 : realT minsep ///< [in]
72 : );
73 :
74 : /// Clear the working memory
75 : void clear();
76 :
77 : /// Set the SNR threshold
78 : void treshold( const realT &thresh /**< [in] the new threshold*/ );
79 :
80 : /// Get the SNR threshold
81 : /** \returns the current value of m_threshold;
82 : */
83 : realT threshold();
84 :
85 : /// Set the minimum separation
86 : void minsep( const realT &msep /**< [in] the new minimum separation*/ );
87 :
88 : /// Get the minimum separation between stars
89 : /** \returns the current value of m_minsep;
90 : */
91 : realT minsep();
92 :
93 : /// Set the useMedian flag
94 : void useMedian( const bool &useMed /**< [in] the new useMedian flag value*/ );
95 :
96 : /// Get the useMedian flag
97 : /** \returns the current value of m_useMedian;
98 : */
99 : bool useMedian();
100 :
101 : /// Find stars in an image
102 : /** Returns the coordinates as (row, col) of the highest SNR pixel in each unique star.
103 : *
104 : * \returns a vector of pixel coordinates
105 : */
106 : std::vector<pixelT> operator()( const eigenImage<realT> &im /**< [in] the image to search for stars*/ );
107 : };
108 :
109 : template <typename realT>
110 0 : sourceFinder<realT>::sourceFinder()
111 : {
112 0 : }
113 :
114 : template <typename realT>
115 0 : sourceFinder<realT>::sourceFinder( realT thresh, realT minsep ) : m_threshold( thresh ), m_minsep( minsep )
116 : {
117 0 : }
118 :
119 : template <typename realT>
120 0 : void sourceFinder<realT>::treshold( const realT &thresh )
121 : {
122 0 : m_threshold = thresh;
123 0 : }
124 :
125 : template <typename realT>
126 0 : void sourceFinder<realT>::clear()
127 : {
128 0 : m_vals.clear();
129 0 : m_pixs.clear();
130 0 : m_snr.clear();
131 0 : }
132 :
133 : template <typename realT>
134 0 : realT sourceFinder<realT>::threshold()
135 : {
136 0 : return m_threshold;
137 : }
138 :
139 : template <typename realT>
140 0 : void sourceFinder<realT>::minsep( const realT &msep )
141 : {
142 0 : m_minsep = msep;
143 0 : }
144 :
145 : template <typename realT>
146 0 : realT sourceFinder<realT>::minsep()
147 : {
148 0 : return m_minsep;
149 : }
150 :
151 : template <typename realT>
152 0 : void sourceFinder<realT>::useMedian( const bool &useMed )
153 : {
154 0 : m_useMedian = useMed;
155 0 : }
156 :
157 : template <typename realT>
158 0 : bool sourceFinder<realT>::useMedian()
159 : {
160 0 : return m_useMedian;
161 : }
162 :
163 : template <typename realT>
164 0 : std::vector<typename sourceFinder<realT>::pixelT> sourceFinder<realT>::operator()( const eigenImage<realT> &im )
165 : {
166 0 : m_vals.resize( im.rows() * im.cols() );
167 0 : m_pixs.resize( im.rows() * im.cols() );
168 :
169 0 : int p = 0;
170 0 : for( int cc = 0; cc < im.cols(); ++cc )
171 : {
172 0 : for( int rr = 0; rr < im.rows(); ++rr )
173 : {
174 0 : m_vals[p] = im( rr, cc );
175 0 : m_pixs[p] = { rr, cc };
176 0 : ++p;
177 : }
178 : }
179 :
180 : realT mm;
181 :
182 0 : if( m_useMedian )
183 0 : mm = math::vectorMedian( m_vals );
184 : else
185 0 : mm = math::vectorMean( m_vals );
186 0 : realT std = sqrt( math::vectorVariance( m_vals, mm ) );
187 :
188 : // Get the SNR for any pixels above the threshold
189 0 : m_snr.clear();
190 0 : for( size_t n = 0; n < m_vals.size(); ++n )
191 : {
192 0 : if( ( m_vals[n] - mm ) / std >= m_threshold )
193 : {
194 0 : m_snr.push_back( n );
195 : }
196 : }
197 :
198 : // Delete the pixels above threshold which are with minsep of a higher pixel
199 0 : realT mms2 = m_minsep * m_minsep; // save an op
200 0 : for( long n = 0; n < m_snr.size(); ++n )
201 : {
202 0 : for( long m = 0; m < m_snr.size(); ++m )
203 : {
204 0 : if( m == n )
205 0 : continue;
206 :
207 : // Squared distance:
208 0 : realT d2 = pow( m_pixs[m_snr[n]].first - m_pixs[m_snr[m]].first, 2 ) +
209 0 : pow( m_pixs[m_snr[n]].second - m_pixs[m_snr[m]].second, 2 );
210 :
211 0 : if( d2 < mms2 )
212 : {
213 0 : if( m_vals[m_snr[n]] < m_vals[m_snr[m]] ) // If this is the higher pixel, delete the lower
214 : {
215 0 : m_snr.erase( m_snr.begin() + n );
216 0 : --n;
217 0 : break; // So we're done comparing to the one we just deleted.
218 : }
219 : else // This is lower, so delete it.
220 : {
221 0 : m_snr.erase( m_snr.begin() + m );
222 0 : if( m < n )
223 0 : --n; // changing size behind
224 0 : --m;
225 0 : continue;
226 : }
227 : }
228 : }
229 : }
230 :
231 0 : std::vector<pixelT> found( m_snr.size() );
232 :
233 0 : for( size_t n = 0; n < m_snr.size(); ++n )
234 : {
235 0 : found[n] = m_pixs[m_snr[n]];
236 : }
237 :
238 0 : return found;
239 : }
240 :
241 : } // namespace improc
242 : } // namespace mx
243 :
244 : #endif // improc_sourceFinder_hpp
|
