Line data Source code
1 : /** \file zernike.hpp
2 : * \author Jared R. Males (jaredmales@gmail.com)
3 : * \brief Working with the Zernike polynomials.
4 : *
5 : * \todo the basic zernike polys should be in math::func.
6 : *
7 : * \ingroup signal_processing_files
8 : *
9 : */
10 :
11 : //***********************************************************************//
12 : // Copyright 2015-2020 Jared R. Males (jaredmales@gmail.com)
13 : //
14 : // This file is part of mxlib.
15 : //
16 : // mxlib is free software: you can redistribute it and/or modify
17 : // it under the terms of the GNU General Public License as published by
18 : // the Free Software Foundation, either version 3 of the License, or
19 : // (at your option) any later version.
20 : //
21 : // mxlib is distributed in the hope that it will be useful,
22 : // but WITHOUT ANY WARRANTY; without even the implied warranty of
23 : // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
24 : // GNU General Public License for more details.
25 : //
26 : // You should have received a copy of the GNU General Public License
27 : // along with mxlib. If not, see <http://www.gnu.org/licenses/>.
28 : //***********************************************************************//
29 :
30 : #ifndef math_zernike_hpp
31 : #define math_zernike_hpp
32 :
33 : #include <cmath>
34 : #include <complex>
35 : #include <vector>
36 :
37 : #include "../mxlib.hpp"
38 : #include "../math/func/bessel.hpp"
39 : #include "../math/func/jinc.hpp"
40 : #include "../math/func/factorial.hpp"
41 : #include "../math/func/sign.hpp"
42 : #include "../math/constants.hpp"
43 :
44 : namespace mx
45 : {
46 : namespace sigproc
47 : {
48 :
49 : /**
50 : * \ingroup zernike_basis
51 : * @{
52 : */
53 :
54 : /// Get the Zernike coefficients n,m corrresponding the Noll index j.
55 : /** Calculates the values of (n,m) for an index j following Noll (1976) \cite noll_1976
56 : * See also: http://en.wikipedia.org/wiki/Zernike_polynomials
57 : *
58 : * If j is odd, this returns m <= 0.
59 : *
60 : *
61 : * \retval 0 on success
62 : * \retval -1 on error (j < 1)
63 : *
64 : */
65 : int noll_nm( int &n, ///< [out] n the radial index of the Zernike polynomial
66 : int &m, ///< [out] m the azimuthal index of the Zernnike polynomial. m < 0 if j odd.
67 : int j ///< [in] j the Noll index, j > 0.
68 : );
69 :
70 : /// Get the Noll index j corresponding to Zernike coefficients n,m
71 : /** Calculates the value j for(n,m) following Noll (1976) \cite noll_1976
72 : * See also: http://en.wikipedia.org/wiki/Zernike_polynomials
73 : *
74 : * \retval >= 0 on success
75 : * \retval -1 on error (n-m odd)
76 : *
77 : */
78 : int noll_j( unsigned n, ///< [in] n the radial index of the Zernike polynomial
79 : int m ///< [in] m the azimuthal index of the Zernnike polynomial.
80 : );
81 :
82 : /// Get the number of Zernikes up to and including a radial order.
83 : /** Calculates the total number of Zernike polynomials through radial order \p n. See Noll (1976) \cite noll_1976
84 : * See also: http://en.wikipedia.org/wiki/Zernike_polynomials
85 : *
86 : * \retval the number of
87 : * \retval -1 on error (n-m odd)
88 : *
89 : */
90 : int nZernRadOrd( unsigned n /**< [n] the radial order */ );
91 :
92 : /// Calculate the coefficients of a Zernike radial polynomial
93 : /**
94 : * \retval 0 on success
95 : * \retval -1 on error
96 : *
97 : * \tparam realT is a real floating type
98 : */
99 : template <typename realT>
100 0 : int zernikeRCoeffs(
101 : std::vector<realT> &c, ///< [out] allocated to length \f$ 0.5(n-m)+1\f$ and filled with the coefficients.
102 : int n, ///< [in] the radial index of the Zernike polynomial.
103 : int m ///< [in] the azimuthal index of the Zernike polynomial.
104 : )
105 : {
106 0 : m = abs( m );
107 :
108 0 : if( n < m )
109 : {
110 0 : internal::mxlib_error_report( error_t::invalidarg, "n cannot be less than m in the Zernike polynomials" );
111 0 : return -1;
112 : }
113 :
114 : // If odd, it's 0.
115 0 : if( ( n - m ) % 2 > 0 )
116 : {
117 0 : c.resize( 1, 0 );
118 0 : return 0;
119 : }
120 :
121 0 : int ul = 0.5 * ( n - m ) + 1;
122 :
123 0 : c.resize( ul );
124 :
125 0 : for( int k = 0; k < ul; ++k )
126 : {
127 0 : c[k] = pow( -1.0, k ) * math::func::factorial<realT>( n - k ) /
128 0 : ( math::func::factorial<realT>( k ) * math::func::factorial<realT>( 0.5 * ( n + m ) - k ) *
129 0 : math::func::factorial<realT>( 0.5 * ( n - m ) - k ) );
130 : }
131 :
132 0 : return 0;
133 : }
134 :
135 : // Explicit instantiations:
136 : extern template int zernikeRCoeffs<float>( std::vector<float> &c, int n, int m );
137 :
138 : extern template int zernikeRCoeffs<double>( std::vector<double> &c, int n, int m );
139 :
140 : extern template int zernikeRCoeffs<long double>( std::vector<long double> &c, int n, int m );
141 :
142 : #ifdef HASQUAD
143 : extern template int zernikeRCoeffs<__float128>( std::vector<__float128> &c, int n, int m );
144 : #endif
145 :
146 : /// Calculate the value of a Zernike radial polynomial at a given separation.
147 : /**
148 : *
149 : * \retval -9999 indicates a possible error
150 : * \retval R the value of the Zernike radial polynomial otherwise
151 : *
152 : * \tparam realT is a real floating type
153 : * \tparam calcRealT is a real floating type used for internal calcs, should be at least double.
154 : */
155 : template <typename realT, typename calcRealT>
156 0 : realT zernikeR( realT rho, ///< [in] the radial coordinate, \f$ 0 \le \rho \le 1 \f$.
157 : int n, ///< [in] the radial index of the Zernike polynomial.
158 : int m, ///< [in] the azimuthal index of the Zernike polynomial.
159 : std::vector<calcRealT> &c /**< [in] contains the radial polynomial coeeficients,
160 : and must be of length \f$ 0.5(n-m)+1\f$. */
161 : )
162 : {
163 0 : m = abs( m );
164 :
165 : // If odd, it's 0.
166 0 : if( ( n - m ) % 2 > 0 )
167 : {
168 0 : c.resize( 1, 0 );
169 0 : return 0.0;
170 : }
171 :
172 0 : if( c.size() != 0.5 * ( n - m ) + 1 )
173 : {
174 0 : internal::mxlib_error_report( error_t::invalidarg, "c vector has incorrect length for n and m." );
175 0 : return -9999;
176 : }
177 :
178 0 : realT R = 0.0;
179 0 : for( size_t k = 0; k < c.size(); ++k )
180 : {
181 0 : R += c[k] * pow( rho, n - 2 * k );
182 : }
183 :
184 0 : return R;
185 : }
186 :
187 : extern template float zernikeR<float, double>( float rho, int n, int m, std::vector<double> &c );
188 :
189 : extern template double zernikeR<double, double>( double rho, int n, int m, std::vector<double> &c );
190 :
191 : extern template long double
192 : zernikeR<long double, long double>( long double rho, int n, int m, std::vector<long double> &c );
193 :
194 : #ifdef HASQUAD
195 : extern template __float128 zernikeR<__float128, __float128>( __float128 rho, int n, int m, std::vector<__float128> &c );
196 : #endif
197 :
198 : /// Calculate the value of a Zernike radial polynomial at a given separation.
199 : /**
200 : * \retval -9999 indicates a possible error
201 : * \retval R the value of the Zernike radial polynomial otherwise
202 : *
203 : * \tparam realT is a real floating type
204 : * \tparam calcRealT is a real floating type used for internal calculations, should be at least double
205 : */
206 : template <typename realT, typename calcRealT>
207 0 : realT zernikeR( realT rho, ///< [in] the radial coordinate, \f$ 0 \le \rho \le 1 \f$.
208 : int n, ///< [in] the radial index of the Zernike polynomial.
209 : int m ///< [in] the azimuthal index of the Zernike polynomial.
210 : )
211 : {
212 0 : m = abs( m );
213 :
214 : // If odd, it's 0.
215 0 : if( ( n - m ) % 2 > 0 )
216 : {
217 0 : return 0.0;
218 : }
219 :
220 0 : std::vector<calcRealT> c;
221 :
222 0 : if( zernikeRCoeffs( c, n, m ) < 0 )
223 0 : return -9999;
224 :
225 0 : return zernikeR<realT, calcRealT>( rho, n, m, c );
226 0 : }
227 :
228 : extern template float zernikeR<float, double>( float rho, int n, int m );
229 :
230 : extern template double zernikeR<double, double>( double rho, int n, int m );
231 :
232 : extern template long double zernikeR<long double, long double>( long double rho, int n, int m );
233 :
234 : #ifdef HASQUAD
235 : extern template __float128 zernikeR<__float128, __float128>( __float128 rho, int n, int m );
236 : #endif
237 :
238 : /// Calculate the value of a Zernike radial polynomial at a given radius and angle.
239 : /**
240 : * \retval -9999 indicates a possible error
241 : * \retval R the value of the Zernike radial polynomial otherwise
242 : *
243 : * \tparam realT is a real floating type
244 : * \tparam calcRealT is a real floating type used for internal calculations, should be at least double
245 : */
246 : template <typename realT, typename calcRealT>
247 0 : realT zernike( realT rho, /**< [in] the radial coordinate, \f$ 0 \le \rho \le 1 \f$.*/
248 : realT phi, /**< [in] the azimuthal angle (in radians)*/
249 : int n, /**< [in] the radial index of the Zernike polynomial.*/
250 : int m, /**< [in] the azimuthal index of the Zernike polynomial.*/
251 : std::vector<calcRealT> &c /**< [in] contains the radial polynomial coeeficients, and
252 : must be of length \f$ 0.5(n-m)+1\f$.*/
253 : )
254 : {
255 : realT azt;
256 :
257 0 : if( n == 0 && m == 0 )
258 : {
259 0 : return 1.0;
260 : }
261 :
262 0 : if( m < 0 )
263 : {
264 0 : azt = math::root_two<realT>() * sin( -m * phi );
265 : }
266 0 : else if( m > 0 )
267 : {
268 0 : azt = math::root_two<realT>() * cos( m * phi );
269 : }
270 : else
271 : {
272 0 : azt = 1.0;
273 : }
274 :
275 0 : return sqrt( (realT)n + 1 ) * zernikeR<realT, calcRealT>( rho, n, m, c ) * azt;
276 : }
277 :
278 : extern template float zernike<float, double>( float rho, float phi, int n, int m, std::vector<double> &c );
279 :
280 : extern template double zernike<double, double>( double rho, double phi, int n, int m, std::vector<double> &c );
281 :
282 : extern template long double
283 : zernike<long double, long double>( long double rho, long double phi, int n, int m, std::vector<long double> &c );
284 :
285 : #ifdef HASQUAD
286 : extern template __float128
287 : zernike<__float128, __float128>( __float128 rho, __float128 phi, int n, int m, std::vector<__float128> &c );
288 : #endif
289 :
290 : /// Calculate the value of a Zernike radial polynomial at a given radius and angle.
291 : /**
292 : *
293 : * \retval -9999 indicates a possible error
294 : * \retval R the value of the Zernike radial polynomial otherwise
295 : *
296 : * \tparam realT is a real floating type
297 : * \tparam calcRealT is a real floating type used for internal calculations, should be at least double
298 : */
299 : template <typename realT, typename calcRealT>
300 0 : realT zernike( realT rho, ///< [in] the radial coordinate, \f$ 0 \le \rho \le 1 \f$.
301 : realT phi, ///< [in] the azimuthal angle (in radians)
302 : int n, ///< [in] the radial index of the Zernike polynomial.
303 : int m ///< [in] the azimuthal index of the Zernike polynomial.
304 : )
305 : {
306 :
307 0 : std::vector<calcRealT> c;
308 :
309 0 : if( zernikeRCoeffs<calcRealT>( c, n, m ) < 0 )
310 0 : return -9999;
311 :
312 0 : return zernike<realT, calcRealT>( rho, phi, n, m, c );
313 0 : }
314 :
315 : extern template float zernike<float, double>( float rho, float phi, int n, int m );
316 :
317 : extern template double zernike<double, double>( double rho, double phi, int n, int m );
318 :
319 : extern template long double zernike<long double, long double>( long double rho, long double phi, int n, int m );
320 :
321 : #ifdef HASQUAD
322 : extern template __float128 zernike<__float128, __float128>( __float128 rho, __float128 phi, int n, int m );
323 : #endif
324 :
325 : /// Calculate the value of a Zernike radial polynomial at a given radius and angle.
326 : /**
327 : * \retval -9999 indicates a possible error
328 : * \retval R the value of the Zernike radial polynomial otherwise
329 : *
330 : * \tparam realT is a real floating type
331 : * \tparam calcRealT is a real floating type used for internal calculations, should be at least double
332 : */
333 : template <typename realT, typename calcRealT>
334 0 : realT zernike( realT rho, ///< [in] the radial coordinate, \f$ 0 \le \rho \le 1 \f$.
335 : realT phi, ///< [in] the azimuthal angle (in radians)
336 : int j ///< [in] the Noll index of the Zernike polynomial.
337 : )
338 : {
339 : int n, m;
340 :
341 : // Get n and m from j
342 0 : if( noll_nm( n, m, j ) < 0 )
343 0 : return -9999;
344 :
345 0 : return zernike<realT, calcRealT>( rho, phi, n, m );
346 : }
347 :
348 : extern template float zernike<float, double>( float rho, float phi, int j );
349 :
350 : extern template double zernike<double, double>( double rho, double phi, int j );
351 :
352 : extern template long double zernike<long double, long double>( long double rho, long double phi, int j );
353 :
354 : #ifdef HASQUAD
355 : extern template __float128 zernike<__float128, __float128>( __float128 rho, __float128 phi, int j );
356 : #endif
357 :
358 : /// Fill in an Eigen-like array with a Zernike polynomial
359 : /** Sets any pixel which is at rad \<= r \< rad+0.5 pixels to rho = 1, to be consistent with mx::circularPupil
360 : *
361 : * \tparam realT is a real floating type
362 : * \tparam calcRealT is a real floating type used for internal calculations, should be at least double
363 : */
364 : template <typename arrayT, typename calcRealT, int overscan = 2>
365 : int zernike( arrayT &arr, /**< [out] allocated array with an Eigen-like interface. The rows()
366 : and cols() members are used to size the polynomial. */
367 : int n, /**< [in] the radial index of the polynomial*/
368 : int m, /**< [in] the azimuthal index of the polynomial*/
369 : typename arrayT::Scalar xcen, /**< [in] the x coordinate of the desired center of the polynomial,
370 : in pixels*/
371 : typename arrayT::Scalar ycen, /**< [in] the y coordinate of the desired center of the polynomial,
372 : in pixels*/
373 : typename arrayT::Scalar rad = -1 /**< [in] the desired radius. If rad \<= 0, then the maximum radius
374 : based on dimensions of m is used.*/
375 : )
376 : {
377 : typedef typename arrayT::Scalar realT;
378 : realT x;
379 : realT y;
380 : realT r, rho;
381 : realT phi;
382 :
383 : std::vector<calcRealT> c;
384 :
385 : if( zernikeRCoeffs( c, n, m ) < 0 )
386 : return -1;
387 :
388 : size_t l0 = arr.rows();
389 : size_t l1 = arr.cols();
390 :
391 : if( rad <= 0 )
392 : rad = 0.5 * std::min( l0 - 1, l1 - 1 );
393 :
394 : for( size_t i = 0; i < l0; ++i )
395 : {
396 : for( size_t j = 0; j < l1; ++j )
397 : {
398 : x = i - xcen;
399 : y = j - ycen;
400 :
401 : r = std::sqrt( x * x + y * y );
402 :
403 : // This is to be consistent with mx::circularPupil while still respecting the Zernike rules
404 : if( r > rad && r <= rad + ( 1.0 / overscan ) )
405 : r = rad;
406 :
407 : rho = r / rad;
408 :
409 : if( rho <= 1.0 )
410 : {
411 : phi = std::atan2( y, x );
412 : arr( i, j ) = zernike( rho, phi, n, m, c );
413 : }
414 : else
415 : {
416 : arr( i, j ) = 0.0;
417 : }
418 : }
419 : }
420 : return 0;
421 : }
422 :
423 : /// Fill in an Eigen-like array with a Zernike polynomial
424 : /** Sets any pixel which is at rad \<= r \<= rad+0.5 pixels to rho = 1, to be consistent with mx::circularPupil
425 : *
426 : *
427 : * \tparam realT is a real floating type
428 : * \tparam calcRealT is a real floating type used for internal calculations, should be at least double
429 : */
430 : template <typename arrayT, typename calcRealT>
431 : int zernike(
432 : arrayT &arr, /**< [out] is the allocated array with an Eigen-like interface. The rows() and
433 : cols() members are used to size the polynomial.*/
434 : int j, /**< [in] is the Noll index of the polynomial*/
435 : typename arrayT::Scalar xcen, /**< [in] is the x coordinate of the desired center of the polynomial, in pixels */
436 : typename arrayT::Scalar ycen, /**< [in] is the y coordinate of the desired center of the polynomial, in pixels*/
437 : typename arrayT::Scalar rad = -1 /**< [in] is the desired radius. If rad \<= 0, then the maximum radius
438 : based on dimensions of m is used.*/
439 : )
440 : {
441 : typedef typename arrayT::Scalar realT;
442 :
443 : int n, m;
444 :
445 : if( noll_nm( n, m, j ) < 0 )
446 : return -1;
447 :
448 : return zernike<arrayT, calcRealT>( arr, n, m, xcen, ycen, rad );
449 : }
450 :
451 : /// Fill in an Eigen-like array with a Zernike polynomial
452 : /** The geometric center of the array, 0.5*(arr.rows()-1), 0.5*(arr.cols()-1), is used as the center.
453 : * Sets any pixel which is at rad \<= r \< rad+0.5 pixels to rho = 1, to be consistent with mx::circularPupil
454 : *
455 : * \tparam realT is a real floating type
456 : * \tparam calcRealT is a real floating type used for internal calculations, should be at least double
457 : */
458 : template <typename arrayT, typename calcRealT>
459 : int zernike( arrayT &arr, /**< [out] allocated array with an Eigen-like interface.
460 : The rows() and cols() members are used to size
461 : the polynomial.*/
462 : int n, /**< [in] the radial index of the polynomial*/
463 : int m, /**< [in] the azimuthal index of the polynomial*/
464 : typename arrayT::Scalar rad = -1 /**< [in] [opt] the desired radius. If rad \<= 0, then the
465 : maximum radius based on dimensions of m is used.*/
466 : )
467 : {
468 : typename arrayT::Scalar xcen = 0.5 * ( arr.rows() - 1.0 );
469 : typename arrayT::Scalar ycen = 0.5 * ( arr.cols() - 1.0 );
470 :
471 : return zernike<arrayT, calcRealT>( arr, n, m, xcen, ycen, rad );
472 : }
473 :
474 : /// Fill in an Eigen-like array with a Zernike polynomial
475 : /** The geometric center of the array, 0.5*(arr.rows()-1), 0.5*(arr.cols()-1), is used as the center.
476 : * Sets any pixel which is at rad \<= r \< rad+0.5 pixels to rho = 1, to be consistent with mx::circularPupil
477 : *
478 : * \tparam arrayT is an Eigen-like array of real floating type
479 : * \tparam calcRealT is a real floating type used for internal calculations, should be at least double
480 : */
481 : template <typename arrayT, typename calcRealT>
482 : int zernike( arrayT &arr, /**< [out] the allocated array with an Eigen-like interface.
483 : The rows() and cols() members are used to size
484 : the polynomial.*/
485 : int j, /**< [in] the Noll index of the polynomial*/
486 : typename arrayT::Scalar rad = -1 /**< [in] [opt] the desired radius. If rad \<= 0, then the maximum
487 : radius based on dimensions of m is used.*/
488 : )
489 : {
490 : typename arrayT::Scalar xcen = 0.5 * ( arr.rows() - 1.0 );
491 : typename arrayT::Scalar ycen = 0.5 * ( arr.cols() - 1.0 );
492 :
493 : return zernike<arrayT, calcRealT>( arr, j, xcen, ycen, rad );
494 : }
495 :
496 : /// Fill in an Eigencube-like array with Zernike polynomials in Noll order
497 : /** The cube is pre-allocated to set the image size and the number of modes.
498 : *
499 : * \returns 0 on success
500 : * \returns -1 on error
501 : *
502 : * \tparam cubeT is an Eigencube-like array with real floating point type
503 : * \tparam calcRealT is a real floating type used for internal calculations, should be at least double
504 : */
505 : template <typename cubeT, typename calcRealT>
506 : int zernikeBasis( cubeT &cube, /**< [in/out] the pre-allocated cube which will be filled
507 : with the Zernike basis*/
508 : typename cubeT::Scalar rad = -1, /**< [in] [opt] the radius of the aperture. If -1 then the full
509 : image size is used.*/
510 : int minj = 2 /**< [in] [opt] the minimum j value to include. The default is j=2,
511 : which skips piston (j=1).*/
512 : )
513 : {
514 : typedef typename cubeT::imageT arrayT;
515 :
516 : typename cubeT::imageT im;
517 :
518 : im.resize( cube.rows(), cube.cols() );
519 :
520 : int rv;
521 : for( int i = 0; i < cube.planes(); ++i )
522 : {
523 : rv = zernike<arrayT, calcRealT>( im, minj + i, rad );
524 :
525 : if( rv < 0 )
526 : {
527 : return rv;
528 : }
529 : cube.image( i ) = im;
530 : }
531 :
532 : return 0;
533 : }
534 :
535 : /// Calculate the square-normed Fourier transform of a Zernike polynomial at position (k,phi)
536 : /** Implements Equation (8) of Noll (1976) \cite noll_1976.
537 : *
538 : * \todo need a more robust jinc_n function for n > 1
539 : *
540 : *
541 : * \returns the value of |Q(k,phi)|^2
542 : *
543 : * \tparam realT is the floating point type used for arithmetic
544 : */
545 : template <typename realT>
546 : std::complex<realT> zernikeQ( realT k, /**< [in] the radial coordinate of normalized spatial frequency. This is in the
547 : \cite noll_1976 convention of cycles-per-radius.*/
548 : realT phi, ///< [in] the azimuthal coordinate of normalized spatial frequency
549 : int n, ///< [in] the Zernike polynomial n
550 : int m ///< [in] the Zernike polynomial m
551 : )
552 : {
553 :
554 : std::complex<realT> Q;
555 :
556 : // sloppy implementation of jinc_n for k ~ 0
557 : if( k < 1e-12 )
558 : {
559 : if( n == 0 )
560 : Q = 1.0;
561 : else
562 : Q = 0.0;
563 : }
564 : else
565 : {
566 : Q = math::func::bessel_j( n + 1, math::two_pi<realT>() * k ) / ( math::pi<realT>() * k );
567 : }
568 :
569 : Q = sqrt( n + 1 ) * Q;
570 :
571 : if( m > 0 ) // Even j (see Noll 1976)
572 : {
573 : Q = Q * pow( -1, 0.5 * ( n - m ) ) * pow( std::complex<realT>( { 0, 1 } ), m ) * sqrt( 2 ) * cos( m * phi );
574 : }
575 : else if( m < 0 ) // Odd j (see Noll 1976) , but m can't really be neg
576 : {
577 : Q = Q * pow( -1, 0.5 * ( n + m ) ) * pow( std::complex<realT>( { 0, 1 } ), -m ) * sqrt( 2 ) * sin( -m * phi );
578 : }
579 : else
580 : {
581 : Q = Q * pow( -1, 0.5 * n );
582 : }
583 :
584 : return Q;
585 : }
586 :
587 : /// Calculate the square-normed Fourier transform of a Zernike polynomial at position (k,phi)
588 : /** Implements Equation (8) of Noll (1976) \cite noll_1976.
589 : *
590 : * \todo need a more robust jinc_n function for n > 1
591 : *
592 : *
593 : * \returns the value of |Q(k,phi)|^2
594 : *
595 : * \tparam realT is the floating point type used for arithmetic
596 : */
597 : template <typename realT>
598 1024 : realT zernikeQNorm( realT k, /**< [in] the radial coordinate of normalized spatial frequency. This is in the
599 : \cite noll_1976 convention of cycles-per-radius.*/
600 : realT phi, ///< [in] the azimuthal coordinate of normalized spatial frequency
601 : int n, ///< [in] the Zernike polynomial n
602 : int m ///< [in] the Zernike polynomial m
603 : )
604 : {
605 :
606 : realT B;
607 :
608 : // sloppy implementation of jinc_n for k ~ 0
609 1024 : if( k < 0.00001 )
610 : {
611 1 : if( n == 0 )
612 : {
613 1 : B = 1.0;
614 : }
615 : else
616 : {
617 0 : B = 0.0;
618 : }
619 : }
620 : else
621 : {
622 1023 : B = math::func::bessel_j( n + 1, math::two_pi<realT>() * k ) / ( math::pi<realT>() * k );
623 : }
624 :
625 1024 : realT Q2 = ( n + 1 ) * ( B * B );
626 :
627 1024 : if( m > 0 ) // Even j (see Noll 1976)
628 : {
629 0 : Q2 = 2 * Q2 * pow( cos( m * phi ), 2 );
630 : }
631 1024 : else if( m < 0 ) // Odd j (see Noll 1976)
632 : {
633 0 : Q2 = 2 * Q2 * pow( sin( -m * phi ), 2 );
634 : }
635 :
636 1024 : return Q2;
637 : }
638 :
639 : extern template float zernikeQNorm<float>( float k, float phi, int n, int m );
640 :
641 : extern template double zernikeQNorm<double>( double k, double phi, int n, int m );
642 :
643 : extern template long double zernikeQNorm<long double>( long double k, long double phi, int n, int m );
644 :
645 : #ifdef HASQUAD
646 : extern template __float128 zernikeQNorm<__float128>( __float128 k, __float128 phi, int n, int m );
647 : #endif
648 :
649 : /// Calculate the square-normed Fourier transform of a Zernike polynomial at position (k,phi)
650 : /** Implements Equation (8) of Noll (1976) \cite noll_1976.
651 : *
652 : * \returns the value of |Q(k,phi)|^2
653 : *
654 : * \tparam realT is the floating point type used for arithmetic
655 : *
656 : */
657 : template <typename realT>
658 : realT zernikeQNorm( realT k, /**< [in] the radial coordinate of normalized spatial frequency. This is in the
659 : \cite noll_1976 convention of cycles-per-radius.*/
660 : realT phi, ///< [in] the azimuthal coordinate of normalized spatial frequency
661 : int j ///< [in] the Zernike polynomial index j (Noll convention)
662 : )
663 : {
664 : int n, m;
665 :
666 : noll_nm( n, m, j );
667 :
668 : return zernikeQNorm( k, phi, n, m );
669 : }
670 :
671 : /// Fill in an Eigen-like array with the square-normed Fourier transform of a Zernike polynomial
672 : /** The array is filled in with \f$\lvert Q(k,\phi)\rvert^2\f$ according to Equation (8) of Noll (1976).
673 : *
674 : *
675 : * \returns 0 on success
676 : * \returns -1 on error
677 : *
678 : * \tparam arrayT is the Eigen-like array type. Arithmetic will be done in arrayT::Scalar.
679 : */
680 : template <typename arrayT>
681 1 : int zernikeQNorm( arrayT &arr, /**< [out] the allocated array. The rows() and cols() members are used to size
682 : the transform.*/
683 : arrayT &k, /**< [in] the normalized spatial frequency magnitude at each pixel, in the Noll (1976)
684 : convention of cycles-per-radius.*/
685 : arrayT &phi, ///< [in] the spatial frequency angle at each pixel
686 : int j ///< [in] the polynomial index in the Noll convention \cite noll_1976
687 : )
688 : {
689 1 : if( arr.rows() != k.rows() || arr.cols() != k.cols() )
690 : {
691 0 : internal::mxlib_error_report( error_t::invalidarg, "output array and input k are not the same size" );
692 0 : return -1;
693 : }
694 :
695 1 : if( arr.rows() != phi.rows() || arr.cols() != phi.cols() )
696 : {
697 0 : internal::mxlib_error_report( error_t::invalidarg, "output array and input phi are not the same size" );
698 0 : return -1;
699 : }
700 :
701 : int n, m;
702 1 : if( noll_nm( n, m, j ) < 0 )
703 0 : return -1; // noll_nm will explain error
704 :
705 33 : for( size_t i = 0; i < arr.rows(); ++i )
706 : {
707 1056 : for( size_t j = 0; j < arr.cols(); ++j )
708 : {
709 1024 : arr( i, j ) = zernikeQNorm( k( i, j ), phi( i, j ), n, m );
710 : }
711 : }
712 1 : return 0;
713 : }
714 :
715 : /// Calculate the spatial power spectrum of Piston
716 : template <typename realT>
717 : realT zernikePPiston( const realT &kD /**< [in] Spatial frequency in diameter units, i.e. cycles per aperture.*/ )
718 : {
719 : return 4 * pow( math::func::jinc( math::pi<realT>() * kD ), 2 );
720 : }
721 :
722 : /// Calculate the spatial power spectrum of Tip \& Tilt
723 : template <typename realT>
724 : realT zernikePTipTilt( const realT &kD /**< [in] Spatial frequency in diameter units, i.e. cycles per aperture.*/ )
725 : {
726 : return 16 * pow( math::func::jincN( 2, math::pi<realT>() * kD ), 2 );
727 : }
728 :
729 : /// Calculate the spatial power spectrum of Defocus
730 : template <typename realT>
731 : realT zernikePDefocus( const realT &kD /**< [in] Spatial frequency in diameter units, i.e. cycles per aperture.*/ )
732 : {
733 : return 12 * pow( math::func::jincN( 3, math::pi<realT>() * kD ), 2 );
734 : }
735 :
736 : /// Calculate the spatial power spectrum of Astigmatism
737 : template <typename realT>
738 : realT zernikePAstig( const realT &kD /**< [in] Spatial frequency in diameter units, i.e. cycles per aperture.*/ )
739 : {
740 : return 24 * pow( math::func::jincN( 3, math::pi<realT>() * kD ), 2 );
741 : }
742 :
743 : /// Calculate the spatial power spectrum of Coma
744 : template <typename realT>
745 : realT zernikePComa( const realT &kD /**< [in] Spatial frequency in diameter units, i.e. cycles per aperture.*/ )
746 : {
747 : return 32 * pow( math::func::jincN( 4, math::pi<realT>() * kD ), 2 );
748 : }
749 :
750 : /// Calculate the spatial power spectrum of Trefoil
751 : template <typename realT>
752 : realT zernikePTrefoil( const realT &kD /**< [in] Spatial frequency in diameter units, i.e. cycles per aperture.*/ )
753 : {
754 : return 32 * pow( math::func::jincN( 4, math::pi<realT>() * kD ), 2 );
755 : }
756 :
757 : /// Get the degrees of correction coefficient for Zernike polynomials in Kolmogorov turbulence.
758 : /** Returns the coefficient from Table IV of Noll (1976) \cite noll_1976 for the given index.
759 : * Given this, the total variance in radians at wavelength \f$\lambda\f$ for a diameter \f$D\f$
760 : * aperture after correcting \f$j\f$ modes can be calculated from
761 : * \f[
762 : * realT c = zernikeModeDOCKolmogorov( j );
763 : * var_lambda = c * pow(D/r_0_lambda, math::five_thirds<realT>)
764 : * \f]
765 : * \returns 0 if noll_j is 0, indicating an error
766 : * \returns the coefficient
767 : */
768 : template <typename realT>
769 : realT zernikeModeDOCKolmogorov( unsigned noll_j /**< [in] the mode index, must be greater than 0 */ )
770 : {
771 : realT c;
772 :
773 : switch( noll_j )
774 : {
775 : case 1:
776 : c = 1.0299;
777 : break;
778 : case 2:
779 : c = 0.582;
780 : break;
781 : case 3:
782 : c = 0.134;
783 : break;
784 : case 4:
785 : c = 0.111;
786 : break;
787 : case 5:
788 : c = 0.0880;
789 : break;
790 : case 6:
791 : c = 0.0648;
792 : break;
793 : case 7:
794 : c = 0.0587;
795 : break;
796 : case 8:
797 : c = 0.0525;
798 : break;
799 : case 9:
800 : c = 0.0463;
801 : break;
802 : case 10:
803 : c = 0.0401;
804 : break;
805 : case 11:
806 : c = 0.0377;
807 : break;
808 : case 12:
809 : c = 0.0352;
810 : break;
811 : case 13:
812 : c = 0.0328;
813 : break;
814 : case 14:
815 : c = 0.0304;
816 : break;
817 : case 15:
818 : c = 0.0279;
819 : break;
820 : case 16:
821 : c = 0.0267;
822 : break;
823 : case 17:
824 : c = 0.0255;
825 : break;
826 : case 18:
827 : c = 0.0243;
828 : break;
829 : case 19:
830 : c = 0.0232;
831 : break;
832 : case 20:
833 : c = 0.0220;
834 : break;
835 : case 21:
836 : c = 0.0208;
837 : break;
838 : default:
839 : if( noll_j == 0 )
840 : {
841 : return 0;
842 : }
843 :
844 : c = 0.2944 * pow( static_cast<realT>( noll_j ), -1 * math::half_root_three<realT>() );
845 : }
846 :
847 : return c;
848 : }
849 :
850 : /// Get the degrees of correction for Zernike polynomials in Kolmogorov turbulence.
851 : /** Returns the degree of correction from Table IV of Noll (1976) \cite noll_1976 for the given index.
852 : * This is the total variance in radians for Fried parameter $r_0$ for a diameter $D$.
853 : * Equivalent to:
854 : * \f[
855 : * realT c = zernikeModeDOCKolmogorov( j );
856 : * var_lambda = c * pow(D/r_0_lambda, math::five_thirds<realT>)
857 : * \f]
858 : * \returns 0 if noll_j is 0, indicating an error
859 : * \returns the coefficient
860 : */
861 : template <typename realT>
862 : realT zernikeModeDOCKolmogorov( unsigned noll_j, /**< [in] the mode index, must be greater than 0 */
863 : realT D, /**< [in] the aperture diameter, in same units as r_0 */
864 : realT r_0 /** <[in] the Fried parameter, in same units as D */
865 : )
866 : {
867 : if( noll_j == 0 )
868 : {
869 : return 0;
870 : }
871 :
872 : return zernikeModeDOCKolmogorov<realT>( noll_j ) * pow( D / r_0, math::five_thirds<realT>() );
873 : }
874 :
875 : /// Get the difference in degrees of correction coefficient for Zernike polynomials in Kolmogorov turbulence.
876 : /** Returns the difference in coefficients from Table IV of Noll (1976) \cite noll_1976 for the given index
877 : * from the previous mode. Given this, the variance in radians-squared at wavelength \f$\lambda\f$ for a
878 : * diameter \f$D\f$ aperture in the \f$j\f$-th mode can be calculated from
879 : * \f[
880 : * realT c = zernikeModeDOCDiffKolmogorov( j );
881 : * var_j_lambda = c * pow(D/r_0_lambda, math::five_thirds<realT>)
882 : * \f]
883 : *
884 : * \returns 0 if noll_j is 0 or 1, indicating an error
885 : * \returns the coefficient difference
886 : */
887 : template <typename realT>
888 : realT zernikeModeDOCDiffKolmogorov( unsigned noll_j /**< [in] the mode index, must be greater than 1 */ )
889 : {
890 : realT c;
891 :
892 : switch( noll_j )
893 : {
894 : case 2:
895 : c = 1.0299 - 0.582;
896 : break;
897 : case 3:
898 : c = 0.582 - 0.134;
899 : break;
900 : case 4:
901 : c = 0.134 - 0.111;
902 : break;
903 : case 5:
904 : c = 0.111 - 0.0880;
905 : break;
906 : case 6:
907 : c = 0.0880 - 0.0648;
908 : break;
909 : case 7:
910 : c = 0.0648 - 0.0587;
911 : break;
912 : case 8:
913 : c = 0.0587 - 0.0525;
914 : break;
915 : case 9:
916 : c = 0.0525 - 0.0463;
917 : break;
918 : case 10:
919 : c = 0.0463 - 0.0401;
920 : break;
921 : case 11:
922 : c = 0.0401 - 0.0377;
923 : break;
924 : case 12:
925 : c = 0.0377 - 0.0352;
926 : break;
927 : case 13:
928 : c = 0.0352 - 0.0328;
929 : break;
930 : case 14:
931 : c = 0.0328 - 0.0304;
932 : break;
933 : case 15:
934 : c = 0.0304 - 0.0279;
935 : break;
936 : case 16:
937 : c = 0.0279 - 0.0267;
938 : break;
939 : case 17:
940 : c = 0.0267 - 0.0255;
941 : break;
942 : case 18:
943 : c = 0.0255 - 0.0243;
944 : break;
945 : case 19:
946 : c = 0.0243 - 0.0232;
947 : break;
948 : case 20:
949 : c = 0.0232 - 0.0220;
950 : break;
951 : case 21:
952 : c = 0.0220 - 0.0208;
953 : break;
954 : case 22:
955 : c = 0.0208 - 0.2944 * pow( static_cast<realT>( 22 ), -1 * math::half_root_three<realT>() );
956 : break;
957 : default:
958 : if( noll_j < 2 )
959 : {
960 : return 0;
961 : }
962 :
963 : c = 0.2944 * ( pow( static_cast<realT>( noll_j - 1 ), -1 * math::half_root_three<realT>() ) -
964 : pow( static_cast<realT>( noll_j ), -1 * math::half_root_three<realT>() ) );
965 : }
966 :
967 : return c;
968 : }
969 :
970 : /// Get the variance for a single Zernike polynomial in Kolmogorov turbulence.
971 : /** Returns the variance from Table IV of Noll (1976) \cite noll_1976 for the given mode index.
972 : * For the $j$-th mode on a diameter $D$ aperture and Fried parameter $r_0$ this is equivalent to calculating
973 : * \f[
974 : * realT c = zernikeModeDOCDiffKolmogorov( j );
975 : * var_j_lambda = c * pow(D/r_0_lambda, math::five_thirds<realT>)
976 : * \f]
977 : *
978 : * \returns 0 if noll_j is 0 or 1, indicating an error
979 : * \returns the variance for the \param noll_j mode in radians squared.
980 : */
981 : template <typename realT>
982 : realT zernikeModeDOCDiffKolmogorov( unsigned noll_j, /**< [in] the mode index, must be greater than 0 */
983 : realT D, /**< [in] the aperture diameter, in same units as r_0 */
984 : realT r_0 /** <[in] the Fried parameter, in same units as D */ )
985 : {
986 : if( noll_j < 2 )
987 : {
988 : return 0;
989 : }
990 :
991 : return zernikeModeDOCDiffKolmogorov<realT>( noll_j ) * pow( D / r_0, math::five_thirds<realT>() );
992 : }
993 :
994 : ///@} signal_processing
995 :
996 : } // namespace sigproc
997 : } // namespace mx
998 :
999 : #endif // math_zernike_hpp
|
