mxlib
c++ tools for analyzing astronomical data and other tasks by Jared R. Males. [git repo]
Loading...
Searching...
No Matches
fraunhoferPropagator_test.cpp
Go to the documentation of this file.
1/** \file fraunhoferPropagator_test.cpp
2 * \brief Tests CPU and CUDA Fraunhofer propagation.
3 */
4#include "../../catch2/catch.hpp"
5#ifdef MXLIB_CUDA
6#include "../../cudaTestUtils.hpp"
7#endif
8
9#define MX_NO_ERROR_REPORTS
10
11#include <iostream>
12
13// #define DEBUG
14
18#include "../../../include/math/constants.hpp"
19
20/// Make an Airy pattern and go back to pupil on CPU
21/**
22 * \ingroup fraunhoferPropagator_unit_tests
23 */
24TEST_CASE( "Make an Airy pattern and go back to pupil on CPU", "[wfp]" )
25{
26 typedef float realT;
27 typedef std::complex<realT> complexT;
28 typedef mx::improc::eigenImage<std::complex<realT>> complexFieldT;
29
30 int wfSz = 512;
31 int pupSz = 128;
33 pupil.resize( pupSz, pupSz );
34 pupil.setZero();
35 mx::improc::maskCircle( pupil, 63.5, 63.5, 63.5, 1 );
36
37 complexFieldT complexPupil, complexFocal;
38 mx::improc::eigenImage<realT> realFocal, realPupil;
39
40 complexPupil.resize( wfSz, wfSz );
41 complexFocal.resize( wfSz, wfSz );
42 realFocal.resize( wfSz, wfSz );
43 realPupil.resize( wfSz, wfSz );
44
46 mx::wfp::makeComplexPupil( complexPupil, pupil, wfSz );
47
48 fi.propagatePupilToFocal( complexFocal, complexPupil );
49
50 mx::wfp::extractIntensityImage( realFocal, 0, complexFocal.rows(), 0, complexFocal.cols(), complexFocal, 0, 0 );
51
52 realT expwr = pupil.square().sum();
53 realT pwr = realFocal.sum();
54
55 REQUIRE_THAT( expwr, Catch::Matchers::WithinAbs( pwr, ( pwr ) * ( 1e-4 ) ) );
56
57 complexPupil.setZero();
58 fi.propagateFocalToPupil( complexPupil, complexFocal );
59
60 mx::wfp::extractIntensityImage( realPupil, 0, complexPupil.rows(), 0, complexPupil.cols(), complexPupil, 0, 0 );
61
62 REQUIRE_THAT( realPupil.sum(), Catch::Matchers::WithinAbs( pupil.sum(), pupil.sum() * 1e-4 ) );
63}
64
65#if defined( MXLIB_CUDA ) || defined( __DOXY_ONLY__ )
66/// Make an Airy pattern and go back to pupil on GPU
67/**
68 * \ingroup fraunhoferPropagator_unit_tests
69 */
70TEST_CASE( "Make an Airy pattern and go back to pupil on GPU", "[wfp]" )
71{
72 if( !mxlibTest::cudaDeviceAvailable() )
73 {
74 WARN( "CUDA runtime is available but no CUDA device is present" );
75 return;
76 }
77
78 typedef float realT;
79 typedef std::complex<realT> complexT;
80 typedef mx::improc::eigenImage<std::complex<realT>> complexFieldT;
81
82 int wfSz = 512;
83 int pupSz = 128;
85 pupil.resize( pupSz, pupSz );
86 pupil.setZero();
87 mx::improc::maskCircle( pupil, 63.5, 63.5, 63.5, 1 );
88
89 complexFieldT complexPupil, complexFocal;
90 mx::improc::eigenImage<realT> realFocal, realPupil;
91
92 complexPupil.resize( wfSz, wfSz );
93 complexFocal.resize( wfSz, wfSz );
94 realFocal.resize( wfSz, wfSz );
95 realPupil.resize( wfSz, wfSz );
96
98 mx::wfp::makeComplexPupil( complexPupil, pupil, wfSz );
99
100 mx::cuda::cudaPtr<complexT> dev_complexPupil, dev_complexFocal;
101
102 dev_complexPupil.upload( complexPupil.data(), complexPupil.rows(), complexPupil.cols() );
103 dev_complexFocal.resize( complexPupil.rows(), complexPupil.cols() );
104
105 BREAD_CRUMB;
106
107 fi.propagatePupilToFocal( dev_complexFocal, dev_complexPupil );
108
109 BREAD_CRUMB;
110
111 dev_complexFocal.download( complexFocal.data() );
112
113 BREAD_CRUMB;
114
115 mx::wfp::extractIntensityImage( realFocal, 0, complexFocal.rows(), 0, complexFocal.cols(), complexFocal, 0, 0 );
116
117 BREAD_CRUMB;
118
119 realT expwr = pupil.square().sum();
120 realT pwr = realFocal.sum();
121
122 REQUIRE_THAT( expwr, Catch::Matchers::WithinAbs( pwr, ( pwr ) * ( 1e-4 ) ) );
123
124 BREAD_CRUMB;
125
126 fi.propagateFocalToPupil( dev_complexPupil, dev_complexFocal );
127
128 BREAD_CRUMB;
129
130 complexPupil.setZero(); // make sure
131 dev_complexPupil.download( complexPupil.data() );
132
133 mx::wfp::extractIntensityImage( realPupil, 0, complexPupil.rows(), 0, complexPupil.cols(), complexPupil, 0, 0 );
134
135 REQUIRE_THAT( realPupil.sum(), Catch::Matchers::WithinAbs( pupil.sum(), pupil.sum() * 1e-4 ) );
136}
137#endif // MXLIB_CUDA
Class to perform Fraunhofer propagation between pupil and focal planes.
error_t propagateFocalToPupil(arrayT &complexPupil, arrayT &complexFocal, bool doCenter=true)
Propagate the wavefront from Focal plane to Pupil plane.
error_t propagatePupilToFocal(arrayT &complexFocal, arrayT &complexPupil, bool doCenter=true)
Propagate the wavefront from the pupil plane to the focal plane.
Tools for using the eigen library for image processing.
Declares and defines a class for Fraunhofer propagation of optical wavefronts.
Eigen::Array< scalarT, -1, -1 > eigenImage
Definition of the eigenImage type, which is an alias for Eigen::Array.
TEST_CASE("Make an Airy pattern and go back to pupil on CPU", "[wfp]")
Make an Airy pattern and go back to pupil on CPU.
void maskCircle(arrayT &m, typename arrayT::Scalar xcen, typename arrayT::Scalar ycen, typename arrayT::Scalar rad, typename arrayT::Scalar val, typename arrayT::Scalar pixbuf=0.5)
Mask a circle in an image.
void makeComplexPupil(arrayOutT &complexPupil, const arrayInT &realPupil, int wavefrontSizePixels)
Create a complex pupil plane wavefront from a real amplitude mask.
Declares and defines functions to work with image masks.