Line data Source code
1 : /** \file fraunhoferPropagator.hpp
2 : * \brief Declares and defines a class for Fraunhofer propagation of optical wavefronts
3 : * \ingroup imaging
4 : * \author Jared R. Males (jaredmales@gmail.com)
5 : *
6 : */
7 :
8 : //***********************************************************************//
9 : // Copyright 2015-2020 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 wfp_fraunhoferPropagator_hpp
28 : #define wfp_fraunhoferPropagator_hpp
29 :
30 : #include "../mxlib.hpp"
31 :
32 : #include "../math/constants.hpp"
33 : #include "imagingArray.hpp"
34 : #include "imagingUtils.hpp"
35 :
36 : #include "../math/ft/fftT.hpp"
37 :
38 : #ifdef MXLIB_CUDA
39 : #include "../math/cuda/templateCuda.hpp"
40 : #include "../math/cuda/cudaPtr.hpp"
41 : #include "../math/cuda/templateCufft.hpp"
42 : #include "../math/cuda/templateCublas.hpp"
43 : #endif
44 :
45 : #ifdef DEBUG
46 : #define BREAD_CRUMB std::cout << "DEBUG: " << __FILE__ << " " << __LINE__ << "\n";
47 : #else
48 : #define BREAD_CRUMB
49 : #endif
50 :
51 : namespace mx
52 : {
53 :
54 : namespace wfp
55 : {
56 :
57 : template <typename wavefrontT, int cudaGPU>
58 : struct fraunhoferPropagatorArrayT;
59 :
60 : template <typename wavefrontT>
61 : struct fraunhoferPropagatorArrayT<wavefrontT, 0>
62 : {
63 : typedef wavefrontT arrayT;
64 : };
65 :
66 : #ifdef MXLIB_CUDA
67 : template <typename wavefrontT>
68 : struct fraunhoferPropagatorArrayT<wavefrontT, 1>
69 : {
70 : typedef mx::cuda::cudaPtr<typename wavefrontT::Scalar> arrayT;
71 : };
72 : #endif
73 :
74 : /// Class to perform Fraunhofer propagation between pupil and focal planes
75 : /** This class uses the FFT to propagate between planes, and normalizes so that flux
76 : * is conserved. For propagation from pupil to focal plane, the pupil wavefront is tilted so that
77 : * the focal-plane image is centered at the geometric center of the array. After propagation from
78 : * focal plane to pupil plane, the pupil plane wavefront is un-tilted to restore the
79 : * pupil to its original position.
80 : *
81 : * \tparam _wavefrontT is an Eigen::Array-like type, with std::complex values.
82 : *
83 : * \todo check if we should reverse the FFT orders to be proper
84 : *
85 : * \ingroup imaging
86 : */
87 : template <typename _wavefrontT, int _cudaGPU = 0>
88 : class fraunhoferPropagator
89 : {
90 :
91 : public:
92 : /// The wavefront data type
93 : typedef _wavefrontT wavefrontT;
94 :
95 : static constexpr int cudaGPU = _cudaGPU;
96 :
97 : /// The complex data type
98 : typedef typename wavefrontT::Scalar complexT;
99 :
100 : /// The real data type
101 : typedef typename wavefrontT::Scalar::value_type realT;
102 :
103 : /// The array data type
104 : typedef typename fraunhoferPropagatorArrayT<wavefrontT, cudaGPU>::arrayT arrayT;
105 :
106 : typedef verbose::d verboseT;
107 :
108 : protected:
109 : /// The size of the wavefront in pixels
110 : int m_wavefrontSizePixels{ 0 };
111 :
112 : realT m_xcen{ 0 }; ///< x-coordinate of focal plane center, in pixels
113 : realT m_ycen{ 0 }; ///< x-coordinate of focal plane center, in pixels
114 :
115 : /// Determines how the image is centered.
116 : /** If 0 (default) it is at 0.5*(wfSz-1), if true it is shifted by 0.5*m_wholePixel in each axis.
117 : */
118 : realT m_wholePixel{ 0 };
119 :
120 : /// Phase screen for tilting the pupil plane so that the focal plane image is centered.
121 : fraunhoferPropagatorArrayT<wavefrontT, cudaGPU>::arrayT m_centerFocal;
122 :
123 : /// Phase screen for un-tilting the pupil plane after propagating from a centered focal plane.
124 : fraunhoferPropagatorArrayT<wavefrontT, cudaGPU>::arrayT m_centerPupil;
125 :
126 : /// FFT object for forward FFTs
127 : math::ft::fftT<complexT, complexT, 2, cudaGPU> m_fft_fwd;
128 :
129 : /// FFT object for backward FFTs
130 : math::ft::fftT<complexT, complexT, 2, cudaGPU> m_fft_back;
131 :
132 : public:
133 : /// Constructor
134 : fraunhoferPropagator();
135 :
136 : /// Destructor
137 : ~fraunhoferPropagator();
138 :
139 : /// Get the value of the wholePixel parameter
140 : /** The wholePixel parameter determines how the image is centered. If 0 (default) it is at 0.5*(wfSz-1), if true it
141 : * is shifted by 0.5*m_wholePixel in each axis.
142 : *
143 : * \returns the value of m_wholePixel
144 : */
145 : int wholePixel();
146 :
147 : /// Set the value of the wholePixel parameter
148 : /** The wholePixel parameter determines how the image is centered. If 0 (default) it is at 0.5*(wfSz-1), if true it
149 : * is shifted by 0.5*m_wholePixel in each axis.
150 : */
151 : void wholePixel( realT wp /**< [in] the new wholePixel value */ );
152 :
153 : /// Apply the shift to a pupil wavefront and apply the normalization.
154 : /** This will center the resultant focal plane image.
155 : * You must have allocated the shift screens first, by calling propagatePupilToFocal, propagateFocalToPupil, or
156 : * setWavefrontSizePixels.
157 : */
158 : template <int ccudaGPU = cudaGPU>
159 : error_t shiftPupil( arrayT &complexPupil, /**< [in.out] the complex pupil plane wavefront to shift*/
160 : typename std::enable_if<ccudaGPU == 0>::type * = 0 );
161 :
162 : /// Apply the shift to a pupil wavefront and apply the normalization.
163 : /** This will center the resultant focal plane image.
164 : * You must have allocated the shift screens first, by calling propagatePupilToFocal, propagateFocalToPupil, or
165 : * setWavefrontSizePixels.
166 : */
167 : template <int ccudaGPU = cudaGPU>
168 : error_t shiftPupil( arrayT &complexPupil, /**< [in.out] the complex pupil plane wavefront to shift*/
169 : typename std::enable_if<ccudaGPU == 1>::type * = 0 );
170 :
171 : /// Apply the shift to a pupil wavefront which will restore it to a centered pupil image, with correct flux.
172 : /** You must have allocated the shift screens first, by calling propagatePupilToFocal, propagateFocalToPupil, or
173 : * setWavefrontSizePixels.
174 : */
175 : template <int ccudaGPU = cudaGPU>
176 : error_t unshiftPupil( arrayT &complexPupil, /**< [in.out] the complex pupil plane wavefront to shift*/
177 : typename std::enable_if<ccudaGPU == 0>::type * = 0 );
178 :
179 : /// Apply the shift to a pupil wavefront which will restore it to a centered pupil image, with correct flux.
180 : /** You must have allocated the shift screens first, by calling propagatePupilToFocal, propagateFocalToPupil, or
181 : * setWavefrontSizePixels.
182 : */
183 : template <int ccudaGPU = cudaGPU>
184 : error_t unshiftPupil( arrayT &complexPupil, /**< [in.out] the complex pupil plane wavefront to shift*/
185 : typename std::enable_if<ccudaGPU == 1>::type * = 0 );
186 :
187 : public:
188 : /// Propagate the wavefront from the pupil plane to the focal plane
189 : /** The pupil plane wavefront (complexPupil) is multiplied by a tilt to place the
190 : * image in the geometric center of the focal plane. This can be prevented by
191 : * setting doCenter to false.
192 : *
193 : */
194 : error_t propagatePupilToFocal( arrayT &complexFocal, /**< [out] the focal plane wavefront. Must be
195 : pre-allocated to same size as complexPupil.*/
196 : arrayT &complexPupil, /**< [in] the pupil plane wavefront. Modified due
197 : to application of centering tilt.*/
198 : bool doCenter = true /**< [in] [optional] set to false to not apply
199 : the centering shift*/
200 : );
201 :
202 : /// Propagate the wavefront from Focal plane to Pupil plane
203 : /** After the fourier transform, the output pupil plane wavefront is de-tilted, restoring it
204 : * to the state prior to calling \ref propagatePupilToFocal. This can be prevented by
205 : * setting doCenter to false.
206 : *
207 : */
208 : error_t propagateFocalToPupil( arrayT &complexPupil, /**< [out] the pupil plane wavefront. Must be
209 : pre-allocated to same size as complexFocal.*/
210 : arrayT &complexFocal, /**< [in] the focal plane wavefront. */
211 : bool doCenter = true /**< [in] [optional] set to false to not apply
212 : the centering shift*/
213 : );
214 :
215 : /// Set the size of the wavefront, in pixels
216 : /** Checks if the size changes, does nothing if no change. Otherwise, calls
217 : * \ref makeShiftPhase to pre-calculate the tilt arrays and plans the FFTs.
218 : *
219 : */
220 : void setWavefrontSizePixels( int wfsPix /**< [in] the desired new size of the wavefront */ );
221 :
222 : protected:
223 : template <int ccudaGPU = cudaGPU>
224 : void
225 : setShiftPhase( complexT *shiftFocal, complexT *shiftPupil, typename std::enable_if<ccudaGPU == 0>::type * = 0 );
226 :
227 : template <int ccudaGPU = cudaGPU>
228 : void
229 : setShiftPhase( complexT *shiftFocal, complexT *shiftPupil, typename std::enable_if<ccudaGPU == 1>::type * = 0 );
230 :
231 : /// Calculate the complex tilt arrays for centering and normalizing the wavefronts
232 : /**
233 : */
234 : void makeShiftPhase();
235 :
236 : }; // class fraunhoferPropagator
237 :
238 : template <typename wavefrontT, int cudaGPU>
239 1 : fraunhoferPropagator<wavefrontT, cudaGPU>::fraunhoferPropagator()
240 : {
241 1 : }
242 :
243 : template <typename wavefrontT, int cudaGPU>
244 1 : fraunhoferPropagator<wavefrontT, cudaGPU>::~fraunhoferPropagator()
245 : {
246 1 : }
247 :
248 : template <typename wavefrontT, int cudaGPU>
249 : int fraunhoferPropagator<wavefrontT, cudaGPU>::wholePixel()
250 : {
251 : return m_wholePixel;
252 : }
253 :
254 : template <typename wavefrontT, int cudaGPU>
255 0 : void fraunhoferPropagator<wavefrontT, cudaGPU>::wholePixel( realT wp )
256 : {
257 0 : m_wholePixel = wp;
258 :
259 : // Re-make the shift phase if size already set.
260 0 : if( m_wavefrontSizePixels > 0 )
261 : {
262 0 : makeShiftPhase();
263 : }
264 0 : }
265 :
266 : template <typename wavefrontT, int cudaGPU>
267 : template <int ccudaGPU>
268 1 : error_t fraunhoferPropagator<wavefrontT, cudaGPU>::shiftPupil( arrayT &complexPupil,
269 : typename std::enable_if<ccudaGPU == 0>::type * )
270 : {
271 1 : complexPupil *= m_centerFocal;
272 :
273 1 : return error_t::noerror;
274 : }
275 :
276 : template <typename wavefrontT, int cudaGPU>
277 : template <int ccudaGPU>
278 : error_t fraunhoferPropagator<wavefrontT, cudaGPU>::shiftPupil( arrayT &complexPupil,
279 : typename std::enable_if<ccudaGPU == 1>::type * )
280 : {
281 : BREAD_CRUMB;
282 : #ifdef MXLIB_CUDA
283 :
284 : cudaError_t ce = mx::cuda::elementwiseXxY( reinterpret_cast<cuComplex *>( complexPupil.m_devicePtr ),
285 : reinterpret_cast<cuComplex *>( m_centerFocal.m_devicePtr ),
286 : pow( m_wavefrontSizePixels, 2 ) );
287 :
288 : if( ce != cudaSuccess )
289 : {
290 : return internal::mxlib_error_report<verboseT>( cudaError2error_t( ce ), "from mx::cuda::elementwiseXxY" );
291 : }
292 :
293 : #endif
294 :
295 : BREAD_CRUMB;
296 :
297 : return error_t::noerror;
298 : }
299 :
300 : template <typename wavefrontT, int cudaGPU>
301 : template <int ccudaGPU>
302 1 : error_t fraunhoferPropagator<wavefrontT, cudaGPU>::unshiftPupil( arrayT &complexPupil,
303 : typename std::enable_if<ccudaGPU == 0>::type * )
304 : {
305 1 : complexPupil *= m_centerPupil;
306 :
307 1 : return error_t::noerror;
308 : }
309 :
310 : template <typename wavefrontT, int cudaGPU>
311 : template <int ccudaGPU>
312 : error_t fraunhoferPropagator<wavefrontT, cudaGPU>::unshiftPupil( arrayT &complexPupil,
313 : typename std::enable_if<ccudaGPU == 1>::type * )
314 : {
315 : #ifdef MXLIB_CUDA
316 : cudaError_t ce = mx::cuda::elementwiseXxY( reinterpret_cast<cuComplex *>( complexPupil.m_devicePtr ),
317 : reinterpret_cast<cuComplex *>( m_centerPupil.m_devicePtr ),
318 : pow( m_wavefrontSizePixels, 2 ) );
319 :
320 : if( ce != cudaSuccess )
321 : {
322 : return internal::mxlib_error_report<verboseT>( cudaError2error_t( ce ), "from mx::cuda::elementwiseXxY" );
323 : }
324 : #endif
325 :
326 : return error_t::noerror;
327 : }
328 :
329 : template <typename wavefrontT, int cudaGPU>
330 1 : error_t fraunhoferPropagator<wavefrontT, cudaGPU>::propagatePupilToFocal( arrayT &complexFocal,
331 : arrayT &complexPupil,
332 : bool doCenter /*default = true*/
333 : )
334 : {
335 :
336 : BREAD_CRUMB;
337 :
338 : // First setup the tilt screens (does nothing if there's no change in size)
339 1 : setWavefrontSizePixels( complexPupil.rows() );
340 :
341 : BREAD_CRUMB;
342 :
343 : // Apply the centering shift -- this adjusts by 0.5 pixels and normalizes
344 1 : if( doCenter )
345 : {
346 1 : error_t ec = shiftPupil( complexPupil );
347 :
348 1 : if(ec != error_t::noerror)
349 : {
350 0 : return internal::mxlib_error_report<verboseT>(ec, "from shiftPupil");
351 : }
352 : }
353 :
354 : BREAD_CRUMB;
355 :
356 : // fft_fwd.fft(complexPupil.data(), complexFocal.data() );
357 1 : m_fft_fwd( complexFocal.data(), complexPupil.data() );
358 :
359 : BREAD_CRUMB;
360 :
361 1 : return error_t::noerror;
362 : }
363 :
364 : template <typename wavefrontT, int cudaGPU>
365 1 : error_t fraunhoferPropagator<wavefrontT, cudaGPU>::propagateFocalToPupil( arrayT &complexPupil,
366 : arrayT &complexFocal,
367 : bool doCenter /*default = true*/
368 : )
369 : {
370 : BREAD_CRUMB;
371 :
372 : // First setup the tilt screens (does nothing if there's no change in size)
373 1 : setWavefrontSizePixels( complexPupil.rows() );
374 :
375 : BREAD_CRUMB;
376 : // fft_back.fft( complexFocal.data(), complexPupil.data());
377 1 : m_fft_back( complexPupil.data(), complexFocal.data() );
378 :
379 : BREAD_CRUMB;
380 :
381 : // Unshift the wavefront and normalize
382 1 : if( doCenter )
383 : {
384 1 : error_t ec = unshiftPupil( complexPupil );
385 :
386 1 : if(ec != error_t::noerror)
387 : {
388 0 : return internal::mxlib_error_report<verboseT>(ec, "from unshiftPupil");
389 : }
390 :
391 : }
392 :
393 : BREAD_CRUMB;
394 :
395 1 : return error_t::noerror;
396 : }
397 :
398 : template <typename wavefrontT, int cudaGPU>
399 2 : void fraunhoferPropagator<wavefrontT, cudaGPU>::setWavefrontSizePixels( int wfsPix )
400 : {
401 : BREAD_CRUMB;
402 :
403 : // If no change in size, do nothing
404 2 : if( wfsPix == m_centerFocal.rows() )
405 : {
406 1 : return;
407 : }
408 :
409 : BREAD_CRUMB;
410 :
411 1 : m_wavefrontSizePixels = wfsPix;
412 :
413 : BREAD_CRUMB;
414 :
415 1 : m_xcen = 0.5 * ( wfsPix - 1.0 );
416 1 : m_ycen = 0.5 * ( wfsPix - 1.0 );
417 :
418 : BREAD_CRUMB;
419 :
420 1 : makeShiftPhase();
421 :
422 : BREAD_CRUMB;
423 :
424 1 : m_fft_fwd.plan( wfsPix, wfsPix );
425 :
426 : BREAD_CRUMB;
427 :
428 1 : m_fft_back.plan( wfsPix, wfsPix, math::ft::dir::backward );
429 :
430 : BREAD_CRUMB;
431 : }
432 :
433 : template <typename wavefrontT, int cudaGPU>
434 : template <int ccudaGPU>
435 1 : void fraunhoferPropagator<wavefrontT, cudaGPU>::setShiftPhase( complexT *centerFocal,
436 : complexT *centerPupil,
437 : typename std::enable_if<ccudaGPU == 0>::type * )
438 : {
439 1 : m_centerFocal.resize( m_wavefrontSizePixels, m_wavefrontSizePixels );
440 1 : m_centerPupil.resize( m_wavefrontSizePixels, m_wavefrontSizePixels );
441 :
442 : // Have to do this manually b/c "Eigen-like" might not support copy from a map
443 513 : for( int cc = 0; cc < m_wavefrontSizePixels; ++cc )
444 : {
445 262656 : for( int rr = 0; rr < m_wavefrontSizePixels; ++rr )
446 : {
447 262144 : m_centerFocal( rr, cc ) = centerFocal[cc * m_wavefrontSizePixels + rr];
448 262144 : m_centerPupil( rr, cc ) = centerPupil[cc * m_wavefrontSizePixels + rr];
449 : }
450 : }
451 1 : }
452 :
453 : template <typename wavefrontT, int cudaGPU>
454 : template <int ccudaGPU>
455 : void fraunhoferPropagator<wavefrontT, cudaGPU>::setShiftPhase( complexT *centerFocal,
456 : complexT *centerPupil,
457 : typename std::enable_if<ccudaGPU == 1>::type * )
458 : {
459 : #ifdef MXLIB_CUDA
460 :
461 : BREAD_CRUMB;
462 :
463 : m_centerFocal.upload( centerFocal, m_wavefrontSizePixels, m_wavefrontSizePixels );
464 :
465 : BREAD_CRUMB;
466 :
467 : m_centerPupil.upload( centerPupil, m_wavefrontSizePixels, m_wavefrontSizePixels );
468 :
469 : BREAD_CRUMB;
470 :
471 : #endif
472 : }
473 :
474 : template <typename wavefrontT, int cudaGPU>
475 1 : void fraunhoferPropagator<wavefrontT, cudaGPU>::makeShiftPhase()
476 : {
477 1 : constexpr realT pi = math::pi<realT>();
478 :
479 : // The normalization is included in the tilt.
480 1 : realT norm = 1. / ( m_wavefrontSizePixels * sqrt( 2 ) );
481 1 : complexT cnorm = complexT( norm, norm );
482 :
483 : BREAD_CRUMB;
484 :
485 : /// Host memory to build the shift screens
486 262145 : complexT *centerFocal = new complexT[m_wavefrontSizePixels * m_wavefrontSizePixels];
487 :
488 262145 : complexT *centerPupil = new complexT[m_wavefrontSizePixels * m_wavefrontSizePixels];
489 :
490 : // Shift by 0.5 pixels
491 1 : realT arg = -2.0 * pi * 0.5 * ( m_wavefrontSizePixels - m_wholePixel ) / ( m_wavefrontSizePixels - 1 );
492 :
493 513 : for( int cc = 0; cc < m_wavefrontSizePixels; ++cc )
494 : {
495 262656 : for( int rr = 0; rr < m_wavefrontSizePixels; ++rr )
496 : {
497 262144 : centerFocal[cc * m_wavefrontSizePixels + rr] =
498 262144 : cnorm * exp( complexT( 0., arg * ( ( rr - m_xcen ) + ( cc - m_ycen ) ) ) );
499 262144 : centerPupil[cc * m_wavefrontSizePixels + rr] =
500 262144 : cnorm * exp( complexT( 0., 0.5 * pi - arg * ( ( rr - m_xcen ) + ( cc - m_ycen ) ) ) );
501 : }
502 : }
503 :
504 : BREAD_CRUMB;
505 :
506 1 : setShiftPhase( centerFocal, centerPupil );
507 :
508 : BREAD_CRUMB;
509 :
510 1 : delete[] centerFocal;
511 1 : delete[] centerPupil;
512 1 : }
513 :
514 : } // namespace wfp
515 : } // namespace mx
516 :
517 : #endif // wfp_fraunhoferPropagator_hpp
|
