Line data Source code
1 : /** \file imagingUtils.hpp
2 : * \brief Utilities for modeling image formation
3 : * \ingroup imaging
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 imagingUtils_hpp
28 : #define imagingUtils_hpp
29 :
30 : #include <cmath>
31 :
32 : #include "../mxlib.hpp"
33 : #include "../improc/eigenImage.hpp"
34 : #include "../math/constants.hpp"
35 :
36 : #include "imagingArray.hpp"
37 :
38 : namespace mx
39 : {
40 :
41 : namespace wfp
42 : {
43 :
44 : /// Calculate the angular plate scale (radians per pixel) of an image after propagation by FFT.
45 : /**
46 : *
47 : * \returns the platescale of the wavefront after propagation by FFT, in radians per pixel.
48 : *
49 : * \tparam realT a real floating point type
50 : *
51 : * \ingroup imaging
52 : */
53 : template <typename realT>
54 : realT fftPlateScale( size_t pixels, ///< [in] the linear dimension of the FFT (including 0 pad, etc.)
55 : realT metersPerPixel, ///< [in] the scale of the input wavefront [m/pix]
56 : realT lambda ///< [in] the wavelength of the wavefront [m]
57 : )
58 : {
59 : return ( lambda / metersPerPixel ) * ( 1. / pixels );
60 : }
61 :
62 : /// Fill in an Eigen-like array with a circular pupil mask.
63 : /** Sets any pixel which is at rad \<= r \< rad+(1.0/overscan) pixels to rho = 1,
64 : *
65 : *
66 : * \retval 0 on success
67 : * \retval -1 on error
68 : *
69 : * \ingroup imaging
70 : */
71 : template <class arrayT>
72 0 : int circularPupil( arrayT &m, ///< [in.out] is the allocated Array. Dimensions are used to create the pupil.
73 : typename arrayT::Scalar eps = 0, ///< [in] [optional] is the central obscuration. 0-1, default is 0.
74 : typename arrayT::Scalar rad = 0, ///< [in] [optional] is the desired radius. If rad \<= 0, then the
75 : ///< maximum radius based on dimensions of m is used. Default is 0.
76 : typename arrayT::Scalar overscan = 0 ///< [in] [optional] overscan in fractional pixels, to include
77 : ///< partial pixels on the edge. Default is 0.
78 : )
79 : {
80 :
81 0 : if( eps < 0 )
82 : {
83 0 : internal::mxlib_error_report( error_t::invalidarg, "Central obscuration can not be < 0." );
84 0 : return -1;
85 : }
86 :
87 0 : if( eps > 1 )
88 : {
89 0 : internal::mxlib_error_report( error_t::invalidarg, "Central obscuration can not be > 1." );
90 0 : return -1;
91 : }
92 :
93 0 : size_t l0 = m.rows();
94 0 : size_t l1 = m.cols();
95 :
96 : typename arrayT::Scalar r;
97 0 : typename arrayT::Scalar xc = 0.5 * ( l0 - 1 );
98 0 : typename arrayT::Scalar yc = 0.5 * ( l1 - 1 );
99 :
100 0 : if( rad <= 0 )
101 0 : rad = 0.5 * std::min( l0 - 1, l1 - 1 );
102 :
103 0 : for( size_t i = 0; i < l0; i++ )
104 : {
105 0 : for( size_t j = 0; j < l1; j++ )
106 : {
107 0 : r = std::sqrt( std::pow( i - xc, 2 ) + std::pow( j - yc, 2 ) );
108 :
109 0 : if( r <= rad + ( overscan ) && r >= eps * rad )
110 0 : m( i, j ) = 1;
111 : else
112 0 : m( i, j ) = 0;
113 : }
114 : }
115 :
116 0 : return 0;
117 : }
118 :
119 : /// Draw a line in an image
120 : /** \todo should handle width much more intelligently, this only works for ~45 degree lines.
121 : *
122 : * \tparam arrayT is an Eigen-like array with public typedef Scalar
123 : */
124 : template <class arrayT>
125 : void drawLine( arrayT &im, ///< [in.out] The input image, modified.
126 : typename arrayT::Scalar x0, ///< [in] the x value, relative to image center, of the starting point
127 : typename arrayT::Scalar y0, ///< [in] the y value, relative to image center, of the starting point
128 : typename arrayT::Scalar x1, ///< [in] the x value, relative to image center, of the end point
129 : typename arrayT::Scalar y1, ///< [in] the y value, relative to image center, of the end point
130 : typename arrayT::Scalar halfWidth ///< [in] the half-width of the line.
131 : )
132 : {
133 :
134 : int d1 = im.rows();
135 : int d2 = im.cols();
136 :
137 : typename arrayT::Scalar xc, yc;
138 :
139 : xc = 0.5 * ( im.rows() - 1 );
140 : yc = 0.5 * ( im.cols() - 1 );
141 :
142 : typename arrayT::Scalar m = ( y1 - y0 ) / ( x1 - x0 );
143 : int y;
144 :
145 : if( x1 > x0 )
146 : {
147 : for( int x = x0; x <= x1; ++x )
148 : {
149 : y = y0 + ( x - x0 ) * m;
150 :
151 : for( int i = 0; i <= halfWidth; ++i )
152 : {
153 : if( x + xc >= 0 && x + xc < d1 && y + yc + i >= 0 && y + yc + i < d2 )
154 : im( x + xc, y + yc + i ) = 0;
155 : if( x + xc >= 0 && x + xc < d1 && y + yc - i >= 0 && y + yc - i < d2 )
156 : im( x + xc, y + yc - i ) = 0;
157 : }
158 : }
159 : }
160 : else
161 : {
162 : for( int x = x1; x <= x0; ++x )
163 : {
164 : y = y0 + ( x - x0 ) * m;
165 :
166 : for( int i = 0; i <= halfWidth; ++i )
167 : {
168 : if( x + xc >= 0 && x + xc < d1 && y + yc + i >= 0 && y + yc + i < d2 )
169 : im( x + xc, y + yc + i ) = 0;
170 : if( x + xc >= 0 && x + xc < d1 && y + yc - i >= 0 && y + yc - i < d2 )
171 : im( x + xc, y + yc - i ) = 0;
172 : }
173 : }
174 : }
175 : }
176 :
177 : /// Create a complex pupil plane wavefront from a real amplitude mask.
178 : /** The input real amplitude mask is placed in the center of a 0-padded complex array.
179 : *
180 : * \param [out] complexPupil the complex pupil plane wavefront
181 : * \param [in] realPupil a real amplitude mask.
182 : * \param [in] wavefrontSizePixels the desired size of the ouput wavefront, should be at least as big as realPupil
183 : *
184 : * \ingroup imaging
185 : */
186 : template <typename arrayOutT, typename arrayInT>
187 1 : void makeComplexPupil( arrayOutT &complexPupil, const arrayInT &realPupil, int wavefrontSizePixels )
188 : {
189 :
190 1 : complexPupil.resize( wavefrontSizePixels, wavefrontSizePixels );
191 1 : complexPupil.setZero();
192 :
193 : // Lower-left corner of insertion region
194 1 : int bl = 0.5 * ( complexPupil.rows() - 1 ) - 0.5 * ( realPupil.rows() - 1. );
195 :
196 129 : for( int i = 0; i < realPupil.cols(); ++i )
197 : {
198 16512 : for( int j = 0; j < realPupil.rows(); ++j )
199 : {
200 16384 : complexPupil( bl + i, bl + j ) =
201 32768 : typename arrayOutT::Scalar( realPupil( i, j ), 0 ); //*exp( typename arrayOutT::Scalar(0,1));
202 : }
203 : }
204 : // complexPupil.block(bl, bl, realPupil.rows(), realPupil.rows()) = realPupil*std::complex<realT>(1,0);
205 1 : }
206 :
207 : /// Create a complex wavefront from a real amplitude and a real phase.
208 : /** The wavefront is placed in the center of a 0-padded complex array.
209 : *
210 : * \param [out] complexWavefront the complex pupil plane wavefront
211 : * \param [in] realAmplitude is the real-valued amplitude.
212 : * \param [in] realPhase is the real-valued phase in radians, same size as realAmplitude
213 : * \param [in] wavefrontSizePixels the desired size of the ouput wavefront, should be at least as big as the real arrays
214 : *
215 : * \ingroup imaging
216 : */
217 : template <typename arrayOutT, typename arrayInT>
218 0 : void makeComplexPupil( arrayOutT &complexWavefront,
219 : const arrayInT &realAmplitude,
220 : const arrayInT &realPhase,
221 : int wavefrontSizePixels )
222 : {
223 :
224 0 : complexWavefront.resize( wavefrontSizePixels, wavefrontSizePixels );
225 0 : complexWavefront.setConstant( typename arrayOutT::Scalar( 0, 0 ) );
226 :
227 : // Lower-left corner of insertion region
228 0 : int bl = 0.5 * ( complexWavefront.rows() - 1 ) - 0.5 * ( realAmplitude.rows() - 1. );
229 :
230 0 : for( int i = 0; i < realAmplitude.cols(); ++i )
231 : {
232 0 : for( int j = 0; j < realAmplitude.rows(); ++j )
233 : {
234 0 : complexWavefront( bl + j, bl + i ) =
235 0 : realAmplitude( j, i ) * exp( ( typename arrayOutT::Scalar( 0, 1 ) ) * realPhase( j, i ) );
236 : }
237 : }
238 0 : }
239 :
240 : /// Apply a tilt to a wavefront
241 : /**
242 : * \param complexWavefront [in.out] the complex wavefront to tilt, will be modified on output
243 : * \param xTilt [input] the amount of tilt in the x direction, in pixels
244 : * \param yTilt [input] the amount of tilt in the y direction, in pixels
245 : *
246 : * \ingroup imaging
247 : */
248 : template <typename wavefrontT>
249 : void tiltWavefront( wavefrontT &complexWavefront,
250 : typename wavefrontT::Scalar::value_type xTilt,
251 : typename wavefrontT::Scalar::value_type yTilt )
252 : {
253 : typedef typename wavefrontT::Scalar complexT;
254 : typedef typename wavefrontT::Scalar::value_type realT;
255 :
256 : realT pi = math::pi<realT>();
257 :
258 : int wfsSizeX = complexWavefront.cols();
259 : int wfsSizeY = complexWavefront.rows();
260 :
261 : realT xCen = 0.5 * ( wfsSizeX - 1 );
262 : realT yCen = 0.5 * ( wfsSizeY - 1 );
263 :
264 : realT argX = 2.0 * pi / ( wfsSizeX - 1.0 );
265 : realT argY = 2.0 * pi / ( wfsSizeY - 1.0 );
266 :
267 : for( int ii = 0; ii < wfsSizeX; ++ii )
268 : {
269 : for( int jj = 0; jj < wfsSizeY; ++jj )
270 : {
271 : complexWavefront( ii, jj ) =
272 : complexWavefront( ii, jj ) *
273 : exp( complexT( (realT)0., argX * xTilt * ( ii - xCen ) + argY * yTilt * ( jj - yCen ) ) );
274 : }
275 : }
276 : }
277 :
278 : /// Extract a block from one image and insert it into a second
279 : template <typename imageT1,
280 : typename imageT2>
281 0 : void extractBlock( imageT1 &dest, ///< [in/out] the image in which to place the extracted block. Must be pre-allocated.
282 : int imX0, ///< [in] the x/row-coord of the lower left pixel of the destination region.
283 : int imXsz, ///< [in] the x size (number of rows) of the block
284 : int imY0, ///< [in] the y/col-coord of the lower left pixel of the desination region
285 : int imYsz, ///< [in] the y=size (number of cols) ov the block
286 : imageT2 &src, ///< [in] the source of the block. Must be large enough.
287 : int wfX0, ///< [in] the x/row-coord of the lower left pixel of the source region
288 : int wfY0 ) ///< [in] the y/col-coord of the lower left pixel of the source region
289 : {
290 0 : int dest_cols = dest.cols();
291 :
292 0 : int src_cols = src.cols();
293 :
294 : typedef typename imageT1::Scalar dataT;
295 :
296 : dataT *dest_data;
297 : dataT *src_data;
298 :
299 0 : for( int j = 0; j < imYsz; ++j )
300 : {
301 0 : dest_data = &dest.data()[imX0 + ( imY0 + j ) * dest_cols];
302 0 : src_data = &src.data()[wfX0 + ( wfY0 + j ) * src_cols];
303 :
304 0 : memcpy( dest_data, src_data, sizeof( dataT ) * imXsz );
305 : }
306 0 : }
307 :
308 : /// Extract a pixels from one image and insert them into a second based on a mask
309 : /** Only pixels with a non-zero value in mask are changed in dest to have the value in src. Other pixels are not
310 : * modified.
311 : */
312 : template <typename imageT1,
313 : typename imageT2,
314 : typename imageT3>
315 : void extractMaskedPixels( imageT1 &dest, ///< [in/out] the image in which to place the extracted pixels. Must be the
316 : ///< same size as src and mask.
317 : const imageT2 &src, ///< [in] the source of the pixels. Must be the same size as mask.
318 : const imageT3 &mask ///< [in] the mask image, where any value other than 0 indicates a pixel
319 : ///< to extract. Must be the same size as src.
320 : )
321 : {
322 : if( dest.rows() != src.rows() )
323 : {
324 : throw mx::exception( error_t::sizeerr, "dest and src do not have same size (rows)" );
325 : }
326 :
327 : if( dest.cols() != src.cols() )
328 : {
329 : throw mx::exception( error_t::sizeerr, "dest and src do not have same size (cols)" );
330 : }
331 :
332 : if( src.rows() != mask.rows() )
333 : {
334 : throw mx::exception( error_t::sizeerr, "src and mask do not have same size (rows)" );
335 : }
336 :
337 : if( src.cols() != mask.cols() )
338 : {
339 : throw mx::exception( error_t::sizeerr, "src and mask do not have same size (cols)" );
340 : }
341 :
342 : for( int cc = 0; cc < dest.cols(); ++cc )
343 : {
344 : for( int rr = 0; rr < dest.rows(); ++rr )
345 : {
346 : if( mask( rr, cc ) != 0 )
347 : {
348 : dest( rr, cc ) = src( rr, cc );
349 : }
350 : }
351 : }
352 : }
353 :
354 : template <typename realImageT, typename complexImageT>
355 2 : void extractIntensityImage(
356 : realImageT &im, int imX0, int imXsz, int imY0, int imYsz, complexImageT &wf, int wfX0, int wfY0 )
357 : {
358 2 : int im_rows = im.cols();
359 :
360 2 : int wf_rows = wf.cols();
361 :
362 : typename realImageT::Scalar *im_data;
363 : typename complexImageT::Scalar *wf_data;
364 :
365 1026 : for( int j = 0; j < imXsz; ++j )
366 : {
367 1024 : im_data = &im.data()[imX0 + ( imY0 + j ) * im_rows];
368 1024 : wf_data = &wf.data()[wfX0 + ( wfY0 + j ) * wf_rows];
369 525312 : for( int i = 0; i < imYsz; ++i )
370 : {
371 524288 : im_data[i] = norm( wf_data[i] );
372 : }
373 : }
374 2 : }
375 :
376 : /// Extract the intensity image from a complex wavefront and accumulate the result
377 : /** Can work on a subset of either or both of the image and the wavefront. This way only
378 : * a small part of the image can be extracted for an oversampled wavefront.
379 : */
380 : template <typename realImageT, typename complexImageT, int cudaGPU = 0>
381 : void extractIntensityImageAccum( realImageT &im, /**<[out] The real image to populate with intensity */
382 : int imX0, /**<[in] The x-coord of the starting image pixel*/
383 : int imXsz, /**<[in] The number of x image pixels*/
384 : int imY0, /**<[in] The y-coord of the starting image pixel */
385 : int imYsz, /**<[in] The number of y image pixels */
386 : complexImageT &wf, /**<[in] The complex wavefront*/
387 : int wfX0, /**<[in] The x-coord of the starting wavefront pixel*/
388 : int wfY0 /**<[in] The y-coord of the starting wavefront pixel*/
389 : );
390 :
391 :
392 :
393 : // CPU versions
394 : template <>
395 : void extractIntensityImageAccum<improc::eigenImage<float>, improc::eigenImage<std::complex<float>>, 0>(
396 : improc::eigenImage<float> &im, int imX0, int imXsz, int imY0, int imYsz, improc::eigenImage<std::complex<float>> &wf, int wfX0, int wfY0 );
397 :
398 : template <>
399 : void extractIntensityImageAccum<improc::eigenImage<double>, improc::eigenImage<std::complex<double>>, 0>(
400 : improc::eigenImage<double> &im, int imX0, int imXsz, int imY0, int imYsz, improc::eigenImage<std::complex<double>> &wf, int wfX0, int wfY0 );
401 :
402 :
403 : #ifdef MXLIB_CUDA
404 :
405 :
406 : // GPU versions
407 : template <>
408 : void extractIntensityImageAccum<cuda::cudaPtr<float>, cuda::cudaPtr<std::complex<float>>, 1>(
409 : cuda::cudaPtr<float> &im, int imX0, int imXsz, int imY0, int imYsz, cuda::cudaPtr<std::complex<float>> &wf, int wfX0, int wfY0 );
410 :
411 : template <>
412 : void extractIntensityImageAccum<cuda::cudaPtr<double>, cuda::cudaPtr<std::complex<double>>, 1>(
413 : cuda::cudaPtr<double> &im, int imX0, int imXsz, int imY0, int imYsz, cuda::cudaPtr<std::complex<double>> &wf, int wfX0, int wfY0 );
414 :
415 : #endif
416 :
417 : } // namespace wfp
418 : } // namespace mx
419 :
420 : #endif //__imagingUtils_hpp__
|
