Line data Source code
1 : /** \file imCenterCircleSym.hpp
2 : * \brief A class to find PSF centers using circular symmetry.
3 : * \ingroup image_processing_files
4 : * \author Jared R. Males (jaredmales@gmail.com)
5 : *
6 : */
7 :
8 : //***********************************************************************//
9 : // Copyright 2015, 2016, 2017 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 __imCenterCircleSym_hpp__
28 : #define __imCenterCircleSym_hpp__
29 :
30 : #include "../math/vectorUtils.hpp"
31 : #include "../math/fit/fitGaussian.hpp"
32 : #include "imageFilters.hpp"
33 : #include "eigenImage.hpp"
34 :
35 : // #define ICCS_OMP
36 : namespace mx
37 : {
38 : namespace improc
39 : {
40 :
41 : /// Find the center of the image of a point source using circular symmetry of the PSF.
42 : /** Conducts a grid search through possible center positions. Finds the point with the minimum
43 : * total variance around concentric circles about that point. The cost grid is fit with a Gaussian.
44 : *
45 : * \tparam realT is a real floating point type.
46 : *
47 : * \ingroup image_reg
48 : */
49 : template <typename realT>
50 : struct imCenterCircleSym
51 : {
52 :
53 : eigenImage<realT> *mask; ///< An optional pointer to a mask. If not NULL and the same size as the input image, any
54 : ///< 0 pixels in this image are ignored.
55 :
56 : eigenImage<realT> grid; ///< The cost grid.
57 : eigenImage<realT> smGrid; ///< The smoothed cost grid.
58 :
59 : realT maxPix; ///< The maximum range of the grid search, the grid will be +/- maxPix. Default: 5.0.
60 : realT dPix; ///< The spacing of the points in the grid search. Default: 0.1.
61 :
62 : realT minRad; ///< The minimum radius to include in the cost function. Default: 1.0.
63 : realT maxRad; ///< The maximum radius to include in the cost function Default: 20.0.
64 : realT dRad; ///< The spacing of the radii. Default: 1.0.
65 :
66 : realT smWidth; ///< FWHM of the Guassian smoothing kernel. Default is 10 (seems to work well regardless of dPix).
67 : ///< Set to 0 for no smoothing.
68 : realT guessWidth; ///< Guess in original pixels for the FWHM of the cost grid. Default is 0.5, which seems to work
69 : ///< well in most cases.
70 :
71 : /// The fitter for finding the centroid of the grid. You can use this to inspect the details of the fit if desired.
72 : // mx::math::fit::fitGaussian2D<mx::gaussian2D_gen_fitter<realT>> fit;
73 : mx::math::fit::fitGaussian2Dgen<realT> fit;
74 :
75 : realT _x0; ///< Internal storage of the guess center position
76 : realT _y0; ///< Internal storage of the guess center position.
77 :
78 0 : imCenterCircleSym()
79 0 : {
80 :
81 0 : mask = 0;
82 :
83 0 : init( 5, 0.1, 1.0, 20.0, 1.0 );
84 :
85 0 : smWidth = 10;
86 0 : guessWidth = 0.5;
87 0 : }
88 :
89 : /// Initialize the grid.
90 0 : int init( realT maxP, ///< [in] the new value of maxPix
91 : realT dP, ///< [in] the new value of dPix
92 : realT minR, ///< [in] the new value of minRad
93 : realT maxR, ///< [in] the new value of maxRad
94 : realT dR ///< [in] the new value of dRad
95 : )
96 : {
97 0 : maxPix = maxP;
98 0 : dPix = dP;
99 0 : minRad = minR;
100 0 : maxRad = maxR;
101 0 : dRad = dR;
102 0 : return 0;
103 : }
104 :
105 : /// Peform the grid search and fit.
106 0 : int center( realT &x, ///< [out] the x-coordinate of the estimated center of rotational symmetry
107 : realT &y, ///< [out] the y-coordinate estimated center of rotational symmetry
108 : eigenImage<realT> &im, ///< [in] the image to centroid
109 : realT x0, ///< [in] an initial guess at the x-coordinate, the grid is centered about this point.
110 : realT y0 ///< [in] an initial guess at the y-coordinate, the grid is centered about this point.
111 : )
112 : {
113 0 : _x0 = x0;
114 0 : _y0 = y0;
115 :
116 0 : bool doMask = false;
117 0 : if( mask )
118 : {
119 0 : if( mask->rows() == im.rows() && mask->cols() == im.cols() )
120 0 : doMask = true;
121 : }
122 :
123 0 : grid.resize( 2 * ( maxPix / dPix ) + 1, 2 * ( maxPix / dPix ) + 1 );
124 :
125 : #ifdef ICCS_OMP
126 : #pragma omp parallel
127 : {
128 : #endif
129 0 : int nRad = ( maxRad - minRad ) / dRad + 1;
130 0 : std::vector<std::vector<realT>> rads;
131 0 : rads.resize( nRad );
132 :
133 0 : int N = 2 * ( maxPix / dPix ) + 1;
134 :
135 : realT startX, startY;
136 : #ifdef ICCS_OMP
137 : #pragma omp for
138 : #endif
139 0 : for( int i = 0; i < N; ++i )
140 : {
141 0 : startX = x0 - maxPix + i * dPix;
142 0 : for( int j = 0; j < N; ++j )
143 : {
144 0 : startY = y0 - maxPix + j * dPix;
145 :
146 0 : for( int k = 0; k < rads.size(); ++k )
147 : {
148 0 : rads[k].clear();
149 : }
150 :
151 0 : for( int ii = startX - maxRad; ii <= startX + maxRad; ++ii )
152 : {
153 0 : for( int jj = startY - maxRad; jj <= startY + maxRad; ++jj )
154 : {
155 :
156 0 : if( doMask )
157 : {
158 0 : if( ( *mask )( ii, jj ) == 0 )
159 0 : continue;
160 : }
161 :
162 : realT rad;
163 0 : rad = sqrt( pow( ii - startX, 2 ) + pow( jj - startY, 2 ) );
164 :
165 0 : if( rad < minRad || rad >= maxRad )
166 0 : continue;
167 :
168 0 : rads[( rad - minRad ) / dRad].push_back( im( ii, jj ) );
169 : } // jj
170 : } // ii
171 :
172 0 : grid( i, j ) = 0;
173 :
174 : realT mean, var;
175 0 : for( int k = 0; k < rads.size(); ++k )
176 : {
177 0 : if( rads[k].size() <= 1 )
178 0 : continue;
179 0 : mean = math::vectorMean( rads[k] );
180 0 : var = math::vectorVariance( rads[k], mean );
181 :
182 0 : grid( i, j ) += var / pow( mean, 2 );
183 : }
184 : } // j
185 : } // i
186 : #ifdef ICCS_OMP
187 : } // omp parallel
188 : #endif
189 : //-------- Now Smooth the grid, and multiply by -1 for fitting -------------
190 :
191 0 : if( smWidth > 0 )
192 : {
193 0 : filterImage( smGrid, grid, gaussKernel<eigenImage<realT>, 2>( 10 ), 0 );
194 : }
195 : else
196 : {
197 0 : smGrid = grid;
198 : }
199 :
200 0 : smGrid *= -1;
201 :
202 : // Get initial guess for the fit.
203 : realT mx, mn, gx, gy;
204 :
205 0 : mx = smGrid.maxCoeff( &gx, &gy );
206 0 : mn = smGrid.minCoeff();
207 :
208 : //------------- Now fit the smoothed cost grid to a 2D elliptical Gaussian ------------
209 : /** \todo This needs to be able to handle highly non-symmetric peaks much better.
210 : */
211 0 : fit.setArray( smGrid.data(), smGrid.rows(), smGrid.cols() );
212 0 : fit.setGuess( mn, mx - mn, gx, gy, guessWidth / dPix, guessWidth / dPix, 0 );
213 :
214 0 : fit.fit();
215 :
216 : // The estimated enter of circular symmetry from the fit.
217 0 : x = x0 - maxPix + fit.x0() * dPix;
218 0 : y = y0 - maxPix + fit.y0() * dPix;
219 :
220 0 : return 0;
221 0 : }
222 :
223 : /// Output the results to a stream
224 : /** Prints a result summary to the input stream.
225 : *
226 : * \tparam iosT is a std::ostream-like type.
227 : * \tparam comment is a comment character to start each line. Can be '\0'.
228 : */
229 : template <typename iosT, char comment = '#'>
230 : iosT &dumpResults( iosT &ios )
231 : {
232 : char c[] = { comment, '\0' };
233 :
234 : ios << c << "--------------------------------------\n";
235 : ios << c << "mx::improc::imCenterCircleSym Results \n";
236 : ios << c << "--------------------------------------\n";
237 : ios << c << "Estimated x-center: " << _x0 - maxPix + fit.x0() * dPix << "\n";
238 : ios << c << "Estimated y-center: " << _y0 - maxPix + fit.y0() * dPix << "\n";
239 : ios << c << "Cost half-width: "
240 : << 0.5 * sigma2fwhm( sqrt( pow( fit.sigma_x(), 2 ) + pow( fit.sigma_y(), 2 ) ) * dPix ) << " pixels\n";
241 : ios << c << "--------------------------------------\n";
242 : ios << c << "Setup:\n";
243 : ios << c << "--------------------------------------\n";
244 : ios << c << "maxPix: " << maxPix << "\n";
245 : ios << c << "dPix: " << dPix << "\n";
246 : ios << c << "minRad: " << minRad << "\n";
247 : ios << c << "maxRad: " << maxRad << "\n";
248 : ios << c << "dRad: " << dRad << "\n";
249 : ios << c << "smWidth: " << smWidth << "\n";
250 : ios << c << "guessWidth: " << guessWidth << "\n";
251 : ios << c << "--------------------------------------\n";
252 : ios << c << "Fit results:\n";
253 : fit.dumpReport( ios );
254 : ios << c << "--------------------------------------\n";
255 :
256 : return ios;
257 : }
258 :
259 : /// Output the results to std::cout
260 : /** Prints a result summary to std::cout
261 : *
262 : * \tparam comment is a comment character to start each line. Can be '\0'.
263 : */
264 : template <char comment = '#'>
265 : std::ostream &dumpResults()
266 : {
267 : return dumpResults<std::ostream, comment>( std::cout );
268 : }
269 : };
270 :
271 : } // namespace improc
272 : } // namespace mx
273 :
274 : #endif //__imCenterCircleSym_hpp__
|
