Line data Source code
1 : /** \file aperturePhotometer.hpp
2 : * \brief Class for conducting aperture photometry on an image.
3 : * \ingroup image_processing_files
4 : * \author Jared R. Males (jaredmales@gmail.com)
5 : *
6 : */
7 :
8 : //***********************************************************************//
9 : // Copyright 2018-2022 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_aperturePhotometer_hpp
28 : #define improc_aperturePhotometer_hpp
29 :
30 : #include <map>
31 :
32 : #include "eigenImage.hpp"
33 : #include "imageMasks.hpp"
34 :
35 : namespace mx
36 : {
37 : namespace improc
38 : {
39 :
40 : /// Class for performing aperture photometry on images.
41 : /** Designed to very efficiently calculate the cumulative flux as a function of radius,
42 : * with minimal cost for performing the measurement on multiple images.
43 : *
44 : * This initial version assumes the source is at 0.5*(rows-1) and 0.5*(cols-1).
45 : *
46 : * First call resize with the desired size. This organizes the radius vector and
47 : * the index image which maps pixels to position in the radius vector.
48 : *
49 : * Then call cumPhot with the image, which will sum the flux contained within each radius.
50 : *
51 : * \ingroup image_processing
52 : *
53 : * \tparam realT the real floating point type of the data
54 : */
55 : template <typename realT>
56 : class aperturePhotometer
57 : {
58 :
59 : protected:
60 : int m_sizeX{ 0 }; ///< The size of the image in X (rows).
61 : int m_sizeY{ 0 }; ///< The size of the image in Y (columns).
62 :
63 : realT m_xcen{ 0 }; ///< The x coordinate of the center of the image. Radii are caculated relative to this point.
64 : ///< Default is 0.5*(sizeX-1).
65 : realT m_ycen{ 0 }; ///< The y coordinate of the center of the image. Radii are caculated relative to this point.
66 : ///< Default is 0.5*(sizeY-1).
67 :
68 : std::vector<realT> m_radius; ///< Holds the ordered unique radii of the pixels in the image.
69 :
70 : eigenImage<size_t> m_indexIm; ///< Maps a pixel in the image to the position in the radius vector.
71 :
72 : bool m_useNanMask{ false };
73 :
74 : eigenImage<realT> *m_nanMask{ nullptr };
75 :
76 : realT m_bgMin{ 0 };
77 : realT m_bgMax{ 0 };
78 :
79 : std::vector<realT> m_bgx;
80 : std::vector<realT> m_bgy;
81 :
82 : public:
83 : /// Resize the photometer, recalculating the radius vector and index image.
84 : /** If size and center are unchanged, nothing is done.
85 : *
86 : * \note Source must be at geometric center of image.
87 : *
88 : * \returns 0 on success
89 : * \returns -1 on error
90 : */
91 : int resize( int sizeX, ///< [in] The new size in rows
92 : int sizeY, ///< [in] The new size in columns
93 : realT xcen, ///< The x coordinate of the center of the image. Radii are caculated relative to this
94 : ///< point. Default is 0.5*(sizeX-1).
95 : realT ycen ///< The y coordinate of the center of the image. Radii are caculated relative to this
96 : ///< point. Default is 0.5*(sizeY-1).
97 : );
98 :
99 : /// Resize the photometer for a centered image, recalculating the radius vector and index image.
100 : /** If size and center are unchanged, nothing is done.
101 : *
102 : * This version uses the geometric center of the image.
103 : *
104 : * \overload
105 : *
106 : * \returns 0 on success
107 : * \returns -1 on error f
108 : */
109 : int resize( int sizeX, ///< [in] The new size in rows
110 : int sizeY ///< [in] The new size in columns
111 : );
112 :
113 : /// Get the cumulative photometry of an image as a function of radius.
114 : /**
115 : *
116 : * \returns 0 on success
117 : * \returns -1 on error
118 : */
119 : int cumPhot(
120 : std::vector<realT> &cumPhot, ///< [out] the cumulative photometry at each point in the radius vector. Resized.
121 : eigenImage<realT> &im, ///< [in] the image on which to perform the calculation
122 : realT maxr =
123 : 0 ///< [in] [optional] the maximum radius to which to calculate. If <= 0 then max possible is used.
124 : );
125 :
126 : /// Get the cumulative photometry of an image as a function of radius.
127 : /**
128 : * \overload
129 : *
130 : * \returns 0 on success
131 : * \returns -1 on error
132 : */
133 : int cumPhot(
134 : std::vector<realT> &cumPhot, ///< [out] the cumulative photometry at each point in the radius vector. Resized.
135 : std::vector<realT> &
136 : deltaPhot, ///< [out] [optional] the photometry at each point in the radius vector, not-cumulative. Resized.
137 : eigenImage<realT> &im, ///< [in] the image on which to perform the calculation
138 : realT maxr =
139 : 0 ///< [in] [optional] the maximum radius to which to calculate. If <= 0 then max possible is used.
140 : );
141 :
142 : protected:
143 : /// Get the cumulative photometry of an image as a function of radius.
144 : /** This is called by the public cumPhot functions.
145 : *
146 : *
147 : * \returns 0 on success
148 : * \returns -1 on error
149 : */
150 : int cumPhotWork(
151 : std::vector<realT> &cumPhot, ///< [out] the cumulative photometry at each point in the radius vector. Resized.
152 : std::vector<realT> *deltaPhot, ///< [out] if not `nullptr`, this vector is filled with the sum at each radius.
153 : eigenImage<realT> &im, ///< [in] the image on which to perform the calculation
154 : realT maxr =
155 : 0 ///< [in] [optional] the maximum radius to which to calculate. If <= 0 then max possible is used.
156 : );
157 :
158 : public:
159 : /// Get the radius value at an index of the vector.
160 0 : realT radius( size_t i /**< [in] the index of the radius vector */ )
161 : {
162 0 : return m_radius[i];
163 : }
164 :
165 0 : void bgMin( realT bgm )
166 : {
167 0 : m_bgMin = bgm;
168 0 : }
169 :
170 0 : void bgMax( realT bgx )
171 : {
172 0 : m_bgMax = bgx;
173 0 : }
174 :
175 0 : void bg( realT bgm, realT bgx )
176 : {
177 0 : m_bgMin = bgm;
178 0 : m_bgMax = bgx;
179 0 : }
180 : };
181 :
182 : template <typename realT>
183 1 : int aperturePhotometer<realT>::resize( int sizeX, int sizeY, realT xcen, realT ycen )
184 : {
185 : // Don't bother if we don't have to.
186 1 : if( sizeX == m_sizeX && sizeY == m_sizeY && xcen == m_xcen && ycen == m_ycen )
187 0 : return 0;
188 :
189 1 : m_sizeX = sizeX;
190 1 : m_sizeY = sizeY;
191 1 : m_xcen = xcen;
192 1 : m_ycen = ycen;
193 :
194 1 : m_indexIm.resize( m_sizeX, m_sizeY );
195 :
196 : // Get radius of each pixel up front, since we need it twice.
197 1 : eigenImage<realT> radIm;
198 1 : radIm.resize( m_sizeX, m_sizeY );
199 1 : radiusImage( radIm, m_xcen, m_ycen ); // mx function to fill in the radius array
200 :
201 1 : std::map<realT, size_t> unrads; // Map of unique radii to their indices in the radius vector.
202 :
203 : // First populate the map
204 4 : for( int i = 0; i < radIm.cols(); ++i )
205 : {
206 12 : for( int j = 0; j < radIm.rows(); ++j )
207 : {
208 : // unrads[radIm(j,i)] = 0;
209 9 : unrads.insert( std::pair<realT, size_t>( radIm( j, i ), 0 ) );
210 : }
211 : }
212 :
213 : // Now fill in the radius vector, and set map values to the indices
214 1 : m_radius.resize( unrads.size() );
215 1 : typename std::map<realT, size_t>::iterator it;
216 1 : int n = 0;
217 4 : for( it = unrads.begin(); it != unrads.end(); ++it )
218 : {
219 3 : m_radius[n] = it->first;
220 3 : it->second = n;
221 3 : ++n;
222 : }
223 : // Finally, fill in the index image with index of the radius vector for that pixel.
224 4 : for( int i = 0; i < radIm.cols(); ++i )
225 : {
226 12 : for( int j = 0; j < radIm.rows(); ++j )
227 : {
228 9 : it = unrads.find( radIm( j, i ) );
229 9 : m_indexIm( j, i ) = it->second;
230 : }
231 : }
232 :
233 1 : if( m_bgMax > m_bgMin && m_bgMin > 0 )
234 : {
235 0 : m_bgx.clear();
236 0 : m_bgy.clear();
237 :
238 0 : for( int cc = 0; cc < radIm.cols(); ++cc )
239 : {
240 0 : for( int rr = 0; rr < radIm.rows(); ++rr )
241 : {
242 0 : if( radIm( rr, cc ) >= m_bgMin && radIm( rr, cc ) <= m_bgMax )
243 : {
244 0 : m_bgx.push_back( rr );
245 0 : m_bgy.push_back( cc );
246 : }
247 : }
248 : }
249 : }
250 :
251 1 : return 0;
252 1 : }
253 :
254 : template <typename realT>
255 1 : int aperturePhotometer<realT>::resize( int sizeX, int sizeY )
256 : {
257 1 : realT xcen = 0.5 * ( 1.0 * sizeX - 1.0 );
258 1 : realT ycen = 0.5 * ( 1.0 * sizeY - 1.0 );
259 :
260 1 : return resize( sizeX, sizeY, xcen, ycen );
261 : }
262 :
263 : template <typename realT>
264 1 : int aperturePhotometer<realT>::cumPhot( std::vector<realT> &cumPhot, eigenImage<realT> &im, realT maxr )
265 : {
266 1 : return cumPhotWork( cumPhot, nullptr, im, maxr );
267 : }
268 :
269 : template <typename realT>
270 0 : int aperturePhotometer<realT>::cumPhot( std::vector<realT> &cumPhot,
271 : std::vector<realT> &deltaPhot,
272 : eigenImage<realT> &im,
273 : realT maxr )
274 : {
275 0 : return cumPhotWork( cumPhot, &deltaPhot, im, maxr );
276 : }
277 :
278 : template <typename realT>
279 1 : int aperturePhotometer<realT>::cumPhotWork( std::vector<realT> &cumPhot,
280 : std::vector<realT> *deltaPhot,
281 : eigenImage<realT> &im,
282 : realT maxr )
283 : {
284 1 : realT xcen = m_xcen; // 0.5*(m_indexIm.rows()-1);
285 1 : realT ycen = m_ycen; // 0.5*(m_indexIm.cols()-1);
286 :
287 1 : if( maxr <= 0 )
288 1 : maxr = m_radius.back();
289 :
290 1 : cumPhot.resize( m_radius.size(), 0 );
291 :
292 : // Make sure we stay in bounds
293 1 : int x0 = xcen - maxr;
294 1 : if( x0 < 0 )
295 0 : x0 = 0;
296 1 : int x1 = xcen + maxr;
297 1 : if( x1 > m_indexIm.rows() - 1 )
298 0 : x1 = m_indexIm.rows() - 1;
299 :
300 1 : int y0 = ycen - maxr;
301 1 : if( y0 < 0 )
302 0 : y0 = 0;
303 1 : int y1 = ycen + maxr;
304 1 : if( y1 > m_indexIm.cols() - 1 )
305 0 : y1 = m_indexIm.rows() - 1;
306 :
307 1 : realT bg = 0;
308 : // If using background annulus, first calc background;
309 1 : if( m_bgx.size() > 0 )
310 : {
311 0 : std::vector<realT> bgann( m_bgx.size() );
312 0 : for( size_t n = 0; n < m_bgx.size(); ++n )
313 : {
314 0 : bgann[n] = im( (int)m_bgx[n], (int)m_bgy[n] );
315 : }
316 :
317 0 : bg = math::vectorMedian( bgann );
318 0 : std::cerr << "background: " << bg << "\n";
319 0 : }
320 :
321 : // First get sum at each unique radius
322 : // Note: switched order of cols/rows indices for speeed
323 1 : if( maxr == m_radius.back() )
324 : {
325 4 : for( int i = y0; i <= y1; ++i )
326 : {
327 12 : for( int j = x0; j <= x1; ++j )
328 : {
329 9 : cumPhot[m_indexIm( j, i )] += im( j, i ) - bg;
330 : }
331 : }
332 : }
333 : else
334 : {
335 0 : for( int i = y0; i <= y1; ++i )
336 : {
337 0 : for( int j = x0; j <= x1; ++j )
338 : {
339 0 : if( m_radius[m_indexIm( j, i )] <= maxr )
340 0 : cumPhot[m_indexIm( j, i )] += im( j, i );
341 : }
342 : }
343 : }
344 :
345 : // If deltaPhot is desired, output it here.
346 1 : if( deltaPhot != nullptr )
347 : {
348 0 : deltaPhot->assign( cumPhot.begin(), cumPhot.end() );
349 : }
350 :
351 : // Now perform cumulative sum
352 3 : for( int i = 1; i < cumPhot.size(); ++i )
353 : {
354 2 : cumPhot[i] += cumPhot[i - 1];
355 : }
356 :
357 1 : return 0;
358 : }
359 :
360 : extern template class aperturePhotometer<float>;
361 :
362 : extern template class aperturePhotometer<double>;
363 :
364 : } // namespace improc
365 : } // namespace mx
366 :
367 : #endif // improc_aperturePhotometer_hpp
|
