Line data Source code
1 : /** \file pupilFitter.hpp
2 : * \brief The MagAO-X Pyramid Pupil Fitter class header
3 : *
4 : * \ingroup pupilFit_files
5 : */
6 :
7 : #ifndef pupilFitter_hpp
8 : #define pupilFitter_hpp
9 :
10 : #include <mx/improc/eigenImage.hpp>
11 : #include <mx/improc/imageUtils.hpp>
12 : #include <mx/improc/circleOuterpix.hpp>
13 : #include <mx/improc/imageTransforms.hpp>
14 :
15 : namespace MagAOX
16 : {
17 : namespace app
18 : {
19 :
20 : /// Struct to perform centration and measure diameter of Pyramid pupils
21 : template <typename realT>
22 : struct pupilFitter
23 : {
24 : mx::improc::eigenImage<realT> m_quad;
25 : mx::improc::eigenImage<realT> m_quadMag;
26 : mx::improc::eigenImage<realT> m_circ;
27 : mx::improc::eigenImage<realT> m_circMag;
28 :
29 : unsigned m_rows{ 0 }; ///< [in] the size of a quad, in rows
30 : unsigned m_cols{ 0 }; ///< [in] the size of a quad, in cols
31 :
32 : int m_numPupils{ 4 }; ///< the number of pupils. Default is 4, 3 is also supported.
33 : float m_bgMedIndex{ 0.1867 };
34 : float m_pupMedIndex{ 0.6867 }; ///< the index of the illuminated pupil median in a sorted array.
35 :
36 : float m_bg[4] = { 0, 0, 0, 0 };
37 : float m_med[4] = { 0, 0, 0, 0 };
38 : float m_avgx[4] = { 0, 0, 0, 0 };
39 : float m_avgy[4] = { 0, 0, 0, 0 };
40 : float m_avgr[4] = { 0, 0, 0, 0 };
41 :
42 : float m_thresh{ 0.5 };
43 :
44 : std::vector<realT> m_pixs;
45 :
46 : pupilFitter()
47 : {
48 : }
49 :
50 : int setSize( unsigned rows, ///< [in] the new size of a quad, in rows
51 : unsigned cols ///< [in] the new size of a quad, in cols
52 : );
53 :
54 : /// Returns the quadrant starting coordinates
55 : int quadCoords( size_t &i0, ///< [out] the i coordinate of the lower-left corner of quad
56 : size_t &j0, ///< [out] the j coordinate of the lower-left corner of auad
57 : int quadNo ///< [in] the quadrant number
58 : );
59 :
60 : int threshold( mx::improc::eigenImage<realT> &im );
61 :
62 : int getQuad( mx::improc::eigenImage<realT> &quad, mx::improc::eigenImage<realT> &im, int quadNo );
63 :
64 : int putQuad( mx::improc::eigenImage<realT> &im, mx::improc::eigenImage<realT> &quad, int quadNo );
65 :
66 : int outerpix( float &avgx, float &avgy, float &avgr, int quadNo );
67 :
68 : int fit( mx::improc::eigenImage<realT> &im, mx::improc::eigenImage<realT> &edged );
69 :
70 : // int emitRegion( const std::string fname );
71 : };
72 :
73 : template <typename realT>
74 0 : int pupilFitter<realT>::setSize( unsigned rows, unsigned cols )
75 : {
76 0 : m_rows = rows;
77 0 : m_cols = cols;
78 :
79 0 : m_quad.resize( m_rows, m_cols );
80 0 : m_circ.resize( m_rows, m_cols );
81 :
82 : // Set up for the magnified version
83 0 : m_quadMag.resize( m_quad.rows() * 10, m_quad.cols() * 10 );
84 0 : m_circMag.resize( m_circ.rows() * 10, m_circ.cols() * 10 );
85 :
86 0 : m_pixs.resize( m_quadMag.rows() * m_quadMag.cols() );
87 :
88 0 : return 0;
89 : }
90 :
91 : template <typename realT>
92 0 : int pupilFitter<realT>::quadCoords( size_t &i0, size_t &j0, int quadNo )
93 : {
94 0 : if( m_numPupils == 3 )
95 : {
96 0 : i0 = 0;
97 0 : j0 = 0;
98 :
99 0 : if( quadNo == 1 )
100 : {
101 0 : i0 = m_rows;
102 : }
103 0 : if( quadNo == 2 )
104 : {
105 0 : i0 = 0.5 * m_rows;
106 0 : j0 = m_cols;
107 : }
108 : }
109 : else
110 : {
111 0 : i0 = 0;
112 0 : j0 = 0;
113 :
114 0 : if( quadNo == 1 )
115 0 : i0 = m_rows;
116 0 : if( quadNo == 2 )
117 0 : j0 = m_cols;
118 0 : if( quadNo == 3 )
119 : {
120 0 : i0 = m_rows;
121 0 : j0 = m_cols;
122 : }
123 : }
124 :
125 0 : return 0;
126 : }
127 :
128 : template <typename realT>
129 0 : int pupilFitter<realT>::threshold( mx::improc::eigenImage<realT> &im )
130 : {
131 0 : for( int i = 0; i < im.rows(); ++i )
132 : {
133 0 : for( int j = 0; j < im.cols(); ++j )
134 : {
135 0 : if( im( i, j ) >= m_thresh )
136 0 : im( i, j ) = 1;
137 : else
138 0 : im( i, j ) = 0;
139 : }
140 : }
141 :
142 0 : return 0;
143 : }
144 :
145 : template <typename realT>
146 0 : int pupilFitter<realT>::getQuad( mx::improc::eigenImage<realT> &quad, mx::improc::eigenImage<realT> &im, int quadNo )
147 : {
148 :
149 0 : if( im.rows() != 2 * m_rows || im.cols() != 2 * m_cols )
150 : {
151 0 : return -1;
152 : }
153 :
154 0 : size_t i0 = 0, j0 = 0;
155 0 : quadCoords( i0, j0, quadNo );
156 :
157 0 : for( size_t i = i0; i < i0 + m_rows; ++i )
158 : {
159 0 : for( size_t j = j0; j < j0 + m_cols; ++j )
160 : {
161 0 : quad( i - i0, j - j0 ) = im( i, j );
162 : }
163 : }
164 :
165 0 : return 0;
166 : }
167 :
168 : template <typename realT>
169 0 : int pupilFitter<realT>::putQuad( mx::improc::eigenImage<realT> &im, mx::improc::eigenImage<realT> &quad, int quadNo )
170 : {
171 0 : if( im.rows() != 2 * m_rows || im.cols() != 2 * m_cols )
172 : {
173 0 : return -1;
174 : }
175 :
176 0 : size_t i0 = 0, j0 = 0;
177 0 : quadCoords( i0, j0, quadNo );
178 :
179 0 : for( size_t i = i0; i < i0 + m_rows; ++i )
180 : {
181 0 : for( size_t j = j0; j < j0 + m_cols; ++j )
182 : {
183 0 : im( i, j ) = quad( i - i0, j - j0 );
184 : ;
185 : }
186 : }
187 :
188 0 : return 0;
189 : }
190 :
191 : template <typename realT>
192 0 : int pupilFitter<realT>::outerpix( float &avgx, float &avgy, float &avgr, int quadNo )
193 : {
194 :
195 : realT x0, y0, avgr0;
196 :
197 0 : mx::improc::circleOuterpix( x0, y0, avgr0, avgx, avgy, avgr, m_circMag, m_quadMag );
198 :
199 0 : avgx /= 10.0;
200 0 : avgy /= 10.0;
201 0 : avgr /= 10.0;
202 :
203 0 : size_t i0 = 0, j0 = 0;
204 0 : quadCoords( i0, j0, quadNo );
205 :
206 0 : avgx += i0;
207 0 : avgy += j0;
208 :
209 0 : return 0;
210 : }
211 :
212 : template <typename realT>
213 0 : int pupilFitter<realT>::fit( mx::improc::eigenImage<realT> &im, mx::improc::eigenImage<realT> &edged )
214 : {
215 : try
216 : {
217 :
218 0 : mx::improc::eigenImage<realT> imin = im;
219 0 : im.setZero();
220 0 : for( int i = 0; i < m_numPupils; ++i )
221 : {
222 : // 0) magnify the image
223 0 : getQuad( m_quad, imin, i );
224 :
225 0 : mx::improc::imageMagnify( m_quadMag, m_quad, mx::improc::bilinearTransform<realT>() );
226 :
227 : // 1) normalize by median of pupil
228 :
229 0 : for( size_t j = 0; j < m_pixs.size(); ++j )
230 : {
231 0 : m_pixs[j] = m_quadMag.data()[j];
232 : }
233 :
234 0 : std::sort( m_pixs.begin(), m_pixs.end() );
235 :
236 0 : m_bg[i] = m_pixs[m_bgMedIndex * m_pixs.size()];
237 :
238 0 : m_med[i] = m_pixs[m_pupMedIndex * m_pixs.size()]; // This should be the median of the pupils if the right
239 : // size and contained in the quad
240 :
241 0 : m_quadMag = ( m_quadMag - m_bg[i] ) / ( m_med[i] - m_bg[i] );
242 :
243 : // 2) Threshold the normalized quad
244 0 : threshold( m_quadMag );
245 :
246 : // 3) Find outer pixels and the radius
247 0 : outerpix( m_avgx[i], m_avgy[i], m_avgr[i], i );
248 :
249 : // 4) De-magnify and prepare for putting on the streams
250 0 : mx::improc::imageRebinSum( m_quad, m_quadMag, true );
251 0 : mx::improc::imageRebinSum( m_circ, m_circMag );
252 0 : putQuad( im, m_quad, i );
253 0 : putQuad( edged, m_circ, i );
254 : }
255 :
256 0 : return 0;
257 0 : }
258 0 : catch( const std::exception &e )
259 : {
260 0 : std::cerr << __FILE__ << ' ' << __LINE__ << ' ' << e.what() << '\n';
261 0 : return -1;
262 : }
263 : }
264 : } // namespace app
265 : } // namespace MagAOX
266 :
267 : #endif // pupilFitter_hpp
|
