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 : * \test Scenario: testing noll_nm \ref tests_sigproc_zernike_noll_nm "[test doc]"
65 : */
66 : int noll_nm( int &n, ///< [out] n the radial index of the Zernike polynomial
67 : int &m, ///< [out] m the azimuthal index of the Zernnike polynomial. m < 0 if j odd.
68 : int j ///< [in] j the Noll index, j > 0.
69 : );
70 :
71 : /// Get the Noll index j corresponding to Zernike coefficients n,m
72 : /** Calculates the value j for(n,m) following Noll (1976) \cite noll_1976
73 : * See also: http://en.wikipedia.org/wiki/Zernike_polynomials
74 : *
75 : * \retval >= 0 on success
76 : * \retval -1 on error (n-m odd)
77 : *
78 : */
79 : int noll_j( unsigned n, ///< [in] n the radial index of the Zernike polynomial
80 : int m ///< [in] m the azimuthal index of the Zernnike polynomial.
81 : );
82 :
83 : /// Get the number of Zernikes up to and including a radial order.
84 : /** Calculates the total number of Zernike polynomials through radial order \p n. See Noll (1976) \cite noll_1976
85 : * See also: http://en.wikipedia.org/wiki/Zernike_polynomials
86 : *
87 : * \retval the number of
88 : * \retval -1 on error (n-m odd)
89 : *
90 : */
91 : int nZernRadOrd( unsigned n /**< [n] the radial order */ );
92 :
93 : /// Calculate the coefficients of a Zernike radial polynomial
94 : /**
95 : * \retval 0 on success
96 : * \retval -1 on error
97 : *
98 : * \tparam realT is a real floating type
99 : */
100 : template <typename realT>
101 : int zernikeRCoeffs(
102 : std::vector<realT> &c, ///< [out] allocated to length \f$ 0.5(n-m)+1\f$ and filled with the coefficients.
103 : int n, ///< [in] the radial index of the Zernike polynomial.
104 : int m ///< [in] the azimuthal index of the Zernike polynomial.
105 : )
106 : {
107 : m = abs( m );
108 :
109 : if( n < m )
110 : {
111 : internal::mxlib_error_report( error_t::invalidarg, "n cannot be less than m in the Zernike polynomials" );
112 : return -1;
113 : }
114 :
115 : // If odd, it's 0.
116 : if( ( n - m ) % 2 > 0 )
117 : {
118 : c.resize( 1, 0 );
119 : return 0;
120 : }
121 :
122 : int ul = 0.5 * ( n - m ) + 1;
123 :
124 : c.resize( ul );
125 :
126 : for( int k = 0; k < ul; ++k )
127 : {
128 : c[k] = pow( -1.0, k ) * math::func::factorial<realT>( n - k ) /
129 : ( math::func::factorial<realT>( k ) * math::func::factorial<realT>( 0.5 * ( n + m ) - k ) *
130 : math::func::factorial<realT>( 0.5 * ( n - m ) - k ) );
131 : }
132 :
133 : return 0;
134 : }
135 :
136 : // Explicit instantiations:
137 : extern template int zernikeRCoeffs<float>( std::vector<float> &c, int n, int m );
138 :
139 : extern template int zernikeRCoeffs<double>( std::vector<double> &c, int n, int m );
140 :
141 : extern template int zernikeRCoeffs<long double>( std::vector<long double> &c, int n, int m );
142 :
143 : #ifdef HASQUAD
144 : extern template int zernikeRCoeffs<__float128>( std::vector<__float128> &c, int n, int m );
145 : #endif
146 :
147 : /// Calculate the value of a Zernike radial polynomial at a given separation.
148 : /**
149 : *
150 : * \retval -9999 indicates a possible error
151 : * \retval R the value of the Zernike radial polynomial otherwise
152 : *
153 : * \tparam realT is a real floating type
154 : * \tparam calcRealT is a real floating type used for internal calcs, should be at least double.
155 : */
156 : template <typename realT, typename calcRealT>
157 : realT zernikeR( realT rho, ///< [in] the radial coordinate, \f$ 0 \le \rho \le 1 \f$.
158 : int n, ///< [in] the radial index of the Zernike polynomial.
159 : int m, ///< [in] the azimuthal index of the Zernike polynomial.
160 : std::vector<calcRealT> &c /**< [in] contains the radial polynomial coeeficients,
161 : and must be of length \f$ 0.5(n-m)+1\f$. */
162 : )
163 : {
164 : m = abs( m );
165 :
166 : // If odd, it's 0.
167 : if( ( n - m ) % 2 > 0 )
168 : {
169 : c.resize( 1, 0 );
170 : return 0.0;
171 : }
172 :
173 : if( c.size() != 0.5 * ( n - m ) + 1 )
174 : {
175 : internal::mxlib_error_report( error_t::invalidarg, "c vector has incorrect length for n and m." );
176 : return -9999;
177 : }
178 :
179 : realT R = 0.0;
180 : for( size_t k = 0; k < c.size(); ++k )
181 : {
182 : R += c[k] * pow( rho, n - 2 * k );
183 : }
184 :
185 : return R;
186 : }
187 :
188 : extern template float zernikeR<float, double>( float rho, int n, int m, std::vector<double> &c );
189 :
190 : extern template double zernikeR<double, double>( double rho, int n, int m, std::vector<double> &c );
191 :
192 : extern template long double
193 : zernikeR<long double, long double>( long double rho, int n, int m, std::vector<long double> &c );
194 :
195 : #ifdef HASQUAD
196 : extern template __float128 zernikeR<__float128, __float128>( __float128 rho, int n, int m, std::vector<__float128> &c );
197 : #endif
198 :
199 : /// Calculate the value of a Zernike radial polynomial at a given separation.
200 : /**
201 : * \retval -9999 indicates a possible error
202 : * \retval R the value of the Zernike radial polynomial otherwise
203 : *
204 : * \tparam realT is a real floating type
205 : * \tparam calcRealT is a real floating type used for internal calculations, should be at least double
206 : */
207 : template <typename realT, typename calcRealT>
208 : realT zernikeR( realT rho, ///< [in] the radial coordinate, \f$ 0 \le \rho \le 1 \f$.
209 : int n, ///< [in] the radial index of the Zernike polynomial.
210 : int m ///< [in] the azimuthal index of the Zernike polynomial.
211 : )
212 : {
213 : m = abs( m );
214 :
215 : // If odd, it's 0.
216 : if( ( n - m ) % 2 > 0 )
217 : {
218 : return 0.0;
219 : }
220 :
221 : std::vector<calcRealT> c;
222 :
223 : if( zernikeRCoeffs( c, n, m ) < 0 )
224 : return -9999;
225 :
226 : return zernikeR<realT, calcRealT>( rho, n, m, c );
227 : }
228 :
229 : extern template float zernikeR<float, double>( float rho, int n, int m );
230 :
231 : extern template double zernikeR<double, double>( double rho, int n, int m );
232 :
233 : extern template long double zernikeR<long double, long double>( long double rho, int n, int m );
234 :
235 : #ifdef HASQUAD
236 : extern template __float128 zernikeR<__float128, __float128>( __float128 rho, int n, int m );
237 : #endif
238 :
239 : /// Calculate the value of a Zernike radial polynomial at a given radius and angle.
240 : /**
241 : * \retval -9999 indicates a possible error
242 : * \retval R the value of the Zernike radial polynomial otherwise
243 : *
244 : * \tparam realT is a real floating type
245 : * \tparam calcRealT is a real floating type used for internal calculations, should be at least double
246 : */
247 : template <typename realT, typename calcRealT>
248 : realT zernike( realT rho, /**< [in] the radial coordinate, \f$ 0 \le \rho \le 1 \f$.*/
249 : realT phi, /**< [in] the azimuthal angle (in radians)*/
250 : int n, /**< [in] the radial index of the Zernike polynomial.*/
251 : int m, /**< [in] the azimuthal index of the Zernike polynomial.*/
252 : std::vector<calcRealT> &c /**< [in] contains the radial polynomial coeeficients, and
253 : must be of length \f$ 0.5(n-m)+1\f$.*/
254 : )
255 : {
256 : realT azt;
257 :
258 : if( n == 0 && m == 0 )
259 : {
260 : return 1.0;
261 : }
262 :
263 : if( m < 0 )
264 : {
265 : azt = math::root_two<realT>() * sin( -m * phi );
266 : }
267 : else if( m > 0 )
268 : {
269 : azt = math::root_two<realT>() * cos( m * phi );
270 : }
271 : else
272 : {
273 : azt = 1.0;
274 : }
275 :
276 : return sqrt( (realT)n + 1 ) * zernikeR<realT, calcRealT>( rho, n, m, c ) * azt;
277 : }
278 :
279 : extern template float zernike<float, double>( float rho, float phi, int n, int m, std::vector<double> &c );
280 :
281 : extern template double zernike<double, double>( double rho, double phi, int n, int m, std::vector<double> &c );
282 :
283 : extern template long double
284 : zernike<long double, long double>( long double rho, long double phi, int n, int m, std::vector<long double> &c );
285 :
286 : #ifdef HASQUAD
287 : extern template __float128
288 : zernike<__float128, __float128>( __float128 rho, __float128 phi, int n, int m, std::vector<__float128> &c );
289 : #endif
290 :
291 : /// Calculate the value of a Zernike radial polynomial at a given radius and angle.
292 : /**
293 : *
294 : * \retval -9999 indicates a possible error
295 : * \retval R the value of the Zernike radial polynomial otherwise
296 : *
297 : * \tparam realT is a real floating type
298 : * \tparam calcRealT is a real floating type used for internal calculations, should be at least double
299 : */
300 : template <typename realT, typename calcRealT>
301 : realT zernike( realT rho, ///< [in] the radial coordinate, \f$ 0 \le \rho \le 1 \f$.
302 : realT phi, ///< [in] the azimuthal angle (in radians)
303 : int n, ///< [in] the radial index of the Zernike polynomial.
304 : int m ///< [in] the azimuthal index of the Zernike polynomial.
305 : )
306 : {
307 :
308 : std::vector<calcRealT> c;
309 :
310 : if( zernikeRCoeffs<calcRealT>( c, n, m ) < 0 )
311 : return -9999;
312 :
313 : return zernike<realT, calcRealT>( rho, phi, n, m, c );
314 : }
315 :
316 : extern template float zernike<float, double>( float rho, float phi, int n, int m );
317 :
318 : extern template double zernike<double, double>( double rho, double phi, int n, int m );
319 :
320 : extern template long double zernike<long double, long double>( long double rho, long double phi, int n, int m );
321 :
322 : #ifdef HASQUAD
323 : extern template __float128 zernike<__float128, __float128>( __float128 rho, __float128 phi, int n, int m );
324 : #endif
325 :
326 : /// Calculate the value of a Zernike radial polynomial at a given radius and angle.
327 : /**
328 : * \retval -9999 indicates a possible error
329 : * \retval R the value of the Zernike radial polynomial otherwise
330 : *
331 : * \tparam realT is a real floating type
332 : * \tparam calcRealT is a real floating type used for internal calculations, should be at least double
333 : */
334 : template <typename realT, typename calcRealT>
335 : realT zernike( realT rho, ///< [in] the radial coordinate, \f$ 0 \le \rho \le 1 \f$.
336 : realT phi, ///< [in] the azimuthal angle (in radians)
337 : int j ///< [in] the Noll index of the Zernike polynomial.
338 : )
339 : {
340 : int n, m;
341 :
342 : // Get n and m from j
343 : if( noll_nm( n, m, j ) < 0 )
344 : return -9999;
345 :
346 : return zernike<realT, calcRealT>( rho, phi, n, m );
347 : }
348 :
349 : extern template float zernike<float, double>( float rho, float phi, int j );
350 :
351 : extern template double zernike<double, double>( double rho, double phi, int j );
352 :
353 : extern template long double zernike<long double, long double>( long double rho, long double phi, int j );
354 :
355 : #ifdef HASQUAD
356 : extern template __float128 zernike<__float128, __float128>( __float128 rho, __float128 phi, int j );
357 : #endif
358 :
359 : /// Fill in an Eigen-like array with a Zernike polynomial
360 : /** Sets any pixel which is at rad \<= r \< rad+0.5 pixels to rho = 1, to be consistent with mx::circularPupil
361 : *
362 : * \tparam realT is a real floating type
363 : * \tparam calcRealT is a real floating type used for internal calculations, should be at least double
364 : */
365 : template <typename arrayT, typename calcRealT, int overscan = 2>
366 : int zernike( arrayT &arr, /**< [out] allocated array with an Eigen-like interface. The rows()
367 : and cols() members are used to size the polynomial. */
368 : int n, /**< [in] the radial index of the polynomial*/
369 : int m, /**< [in] the azimuthal index of the polynomial*/
370 : typename arrayT::Scalar xcen, /**< [in] the x coordinate of the desired center of the polynomial,
371 : in pixels*/
372 : typename arrayT::Scalar ycen, /**< [in] the y coordinate of the desired center of the polynomial,
373 : in pixels*/
374 : typename arrayT::Scalar rad = -1 /**< [in] the desired radius. If rad \<= 0, then the maximum radius
375 : based on dimensions of m is used.*/
376 : )
377 : {
378 : typedef typename arrayT::Scalar realT;
379 : realT x;
380 : realT y;
381 : realT r, rho;
382 : realT phi;
383 :
384 : std::vector<calcRealT> c;
385 :
386 : if( zernikeRCoeffs( c, n, m ) < 0 )
387 : return -1;
388 :
389 : size_t l0 = arr.rows();
390 : size_t l1 = arr.cols();
391 :
392 : if( rad <= 0 )
393 : rad = 0.5 * std::min( l0 - 1, l1 - 1 );
394 :
395 : for( size_t i = 0; i < l0; ++i )
396 : {
397 : for( size_t j = 0; j < l1; ++j )
398 : {
399 : x = i - xcen;
400 : y = j - ycen;
401 :
402 : r = std::sqrt( x * x + y * y );
403 :
404 : // This is to be consistent with mx::circularPupil while still respecting the Zernike rules
405 : if( r > rad && r <= rad + ( 1.0 / overscan ) )
406 : r = rad;
407 :
408 : rho = r / rad;
409 :
410 : if( rho <= 1.0 )
411 : {
412 : phi = std::atan2( y, x );
413 : arr( i, j ) = zernike( rho, phi, n, m, c );
414 : }
415 : else
416 : {
417 : arr( i, j ) = 0.0;
418 : }
419 : }
420 : }
421 : return 0;
422 : }
423 :
424 : /// Fill in an Eigen-like array with a Zernike polynomial
425 : /** Sets any pixel which is at rad \<= r \<= rad+0.5 pixels to rho = 1, to be consistent with mx::circularPupil
426 : *
427 : *
428 : * \tparam realT is a real floating type
429 : * \tparam calcRealT is a real floating type used for internal calculations, should be at least double
430 : */
431 : template <typename arrayT, typename calcRealT>
432 : int zernike(
433 : arrayT &arr, /**< [out] is the allocated array with an Eigen-like interface. The rows() and
434 : cols() members are used to size the polynomial.*/
435 : int j, /**< [in] is the Noll index of the polynomial*/
436 : typename arrayT::Scalar xcen, /**< [in] is the x coordinate of the desired center of the polynomial, in pixels */
437 : typename arrayT::Scalar ycen, /**< [in] is the y coordinate of the desired center of the polynomial, in pixels*/
438 : typename arrayT::Scalar rad = -1 /**< [in] is the desired radius. If rad \<= 0, then the maximum radius
439 : based on dimensions of m is used.*/
440 : )
441 : {
442 : typedef typename arrayT::Scalar realT;
443 :
444 : int n, m;
445 :
446 : if( noll_nm( n, m, j ) < 0 )
447 : return -1;
448 :
449 : return zernike<arrayT, calcRealT>( arr, n, m, xcen, ycen, rad );
450 : }
451 :
452 : /// Fill in an Eigen-like array with a Zernike polynomial
453 : /** The geometric center of the array, 0.5*(arr.rows()-1), 0.5*(arr.cols()-1), is used as the center.
454 : * Sets any pixel which is at rad \<= r \< rad+0.5 pixels to rho = 1, to be consistent with mx::circularPupil
455 : *
456 : * \tparam realT is a real floating type
457 : * \tparam calcRealT is a real floating type used for internal calculations, should be at least double
458 : */
459 : template <typename arrayT, typename calcRealT>
460 : int zernike( arrayT &arr, /**< [out] allocated array with an Eigen-like interface.
461 : The rows() and cols() members are used to size
462 : the polynomial.*/
463 : int n, /**< [in] the radial index of the polynomial*/
464 : int m, /**< [in] the azimuthal index of the polynomial*/
465 : typename arrayT::Scalar rad = -1 /**< [in] [opt] the desired radius. If rad \<= 0, then the
466 : maximum radius based on dimensions of m is used.*/
467 : )
468 : {
469 : typename arrayT::Scalar xcen = 0.5 * ( arr.rows() - 1.0 );
470 : typename arrayT::Scalar ycen = 0.5 * ( arr.cols() - 1.0 );
471 :
472 : return zernike<arrayT, calcRealT>( arr, n, m, xcen, ycen, rad );
473 : }
474 :
475 : /// Fill in an Eigen-like array with a Zernike polynomial
476 : /** The geometric center of the array, 0.5*(arr.rows()-1), 0.5*(arr.cols()-1), is used as the center.
477 : * Sets any pixel which is at rad \<= r \< rad+0.5 pixels to rho = 1, to be consistent with mx::circularPupil
478 : *
479 : * \tparam arrayT is an Eigen-like array of real floating type
480 : * \tparam calcRealT is a real floating type used for internal calculations, should be at least double
481 : */
482 : template <typename arrayT, typename calcRealT>
483 : int zernike( arrayT &arr, /**< [out] the allocated array with an Eigen-like interface.
484 : The rows() and cols() members are used to size
485 : the polynomial.*/
486 : int j, /**< [in] the Noll index of the polynomial*/
487 : typename arrayT::Scalar rad = -1 /**< [in] [opt] the desired radius. If rad \<= 0, then the maximum
488 : radius based on dimensions of m is used.*/
489 : )
490 : {
491 : typename arrayT::Scalar xcen = 0.5 * ( arr.rows() - 1.0 );
492 : typename arrayT::Scalar ycen = 0.5 * ( arr.cols() - 1.0 );
493 :
494 : return zernike<arrayT, calcRealT>( arr, j, xcen, ycen, rad );
495 : }
496 :
497 : /// Fill in an Eigencube-like array with Zernike polynomials in Noll order
498 : /** The cube is pre-allocated to set the image size and the number of modes.
499 : *
500 : * \returns 0 on success
501 : * \returns -1 on error
502 : *
503 : * \tparam cubeT is an Eigencube-like array with real floating point type
504 : * \tparam calcRealT is a real floating type used for internal calculations, should be at least double
505 : */
506 : template <typename cubeT, typename calcRealT>
507 : int zernikeBasis( cubeT &cube, /**< [in/out] the pre-allocated cube which will be filled
508 : with the Zernike basis*/
509 : typename cubeT::Scalar rad = -1, /**< [in] [opt] the radius of the aperture. If -1 then the full
510 : image size is used.*/
511 : int minj = 2 /**< [in] [opt] the minimum j value to include. The default is j=2,
512 : which skips piston (j=1).*/
513 : )
514 : {
515 : typedef typename cubeT::imageT arrayT;
516 :
517 : typename cubeT::imageT im;
518 :
519 : im.resize( cube.rows(), cube.cols() );
520 :
521 : int rv;
522 : for( int i = 0; i < cube.planes(); ++i )
523 : {
524 : rv = zernike<arrayT, calcRealT>( im, minj + i, rad );
525 :
526 : if( rv < 0 )
527 : {
528 : return rv;
529 : }
530 : cube.image( i ) = im;
531 : }
532 :
533 : return 0;
534 : }
535 :
536 : /// Calculate the square-normed Fourier transform of a Zernike polynomial at position (k,phi)
537 : /** Implements Equation (8) of Noll (1976) \cite noll_1976.
538 : *
539 : * \todo need a more robust jinc_n function for n > 1
540 : *
541 : * \test Scenario: testing zernikeQNorm \ref tests_sigproc_zernike_zernikeQNorm "[test doc]"
542 : *
543 : * \returns the value of |Q(k,phi)|^2
544 : *
545 : * \tparam realT is the floating point type used for arithmetic
546 : */
547 : template <typename realT>
548 : std::complex<realT> zernikeQ( realT k, /**< [in] the radial coordinate of normalized spatial frequency. This is in the
549 : \cite noll_1976 convention of cycles-per-radius.*/
550 : realT phi, ///< [in] the azimuthal coordinate of normalized spatial frequency
551 : int n, ///< [in] the Zernike polynomial n
552 : int m ///< [in] the Zernike polynomial m
553 : )
554 : {
555 :
556 : std::complex<realT> Q;
557 :
558 : // sloppy implementation of jinc_n for k ~ 0
559 : if( k < 1e-12 )
560 : {
561 : if( n == 0 )
562 : Q = 1.0;
563 : else
564 : Q = 0.0;
565 : }
566 : else
567 : {
568 : Q = math::func::bessel_j( n + 1, math::two_pi<realT>() * k ) / ( math::pi<realT>() * k );
569 : }
570 :
571 : Q = sqrt( n + 1 ) * Q;
572 :
573 : if( m > 0 ) // Even j (see Noll 1976)
574 : {
575 : Q = Q * pow( -1, 0.5 * ( n - m ) ) * pow( std::complex<realT>( { 0, 1 } ), m ) * sqrt( 2 ) * cos( m * phi );
576 : }
577 : else if( m < 0 ) // Odd j (see Noll 1976) , but m can't really be neg
578 : {
579 : Q = Q * pow( -1, 0.5 * ( n + m ) ) * pow( std::complex<realT>( { 0, 1 } ), -m ) * sqrt( 2 ) * sin( -m * phi );
580 : }
581 : else
582 : {
583 : Q = Q * pow( -1, 0.5 * n );
584 : }
585 :
586 : return Q;
587 : }
588 :
589 : /// Calculate the square-normed Fourier transform of a Zernike polynomial at position (k,phi)
590 : /** Implements Equation (8) of Noll (1976) \cite noll_1976.
591 : *
592 : * \todo need a more robust jinc_n function for n > 1
593 : *
594 : * \test Scenario: testing zernikeQNorm \ref tests_sigproc_zernike_zernikeQNorm "[test doc]"
595 : *
596 : * \returns the value of |Q(k,phi)|^2
597 : *
598 : * \tparam realT is the floating point type used for arithmetic
599 : */
600 : template <typename realT>
601 : realT zernikeQNorm( realT k, /**< [in] the radial coordinate of normalized spatial frequency. This is in the
602 : \cite noll_1976 convention of cycles-per-radius.*/
603 : realT phi, ///< [in] the azimuthal coordinate of normalized spatial frequency
604 : int n, ///< [in] the Zernike polynomial n
605 : int m ///< [in] the Zernike polynomial m
606 : )
607 : {
608 :
609 : realT B;
610 :
611 : // sloppy implementation of jinc_n for k ~ 0
612 : if( k < 0.00001 )
613 : {
614 : if( n == 0 )
615 : {
616 : B = 1.0;
617 : }
618 : else
619 : {
620 : B = 0.0;
621 : }
622 : }
623 : else
624 : {
625 : B = math::func::bessel_j( n + 1, math::two_pi<realT>() * k ) / ( math::pi<realT>() * k );
626 : }
627 :
628 : realT Q2 = ( n + 1 ) * ( B * B );
629 :
630 : if( m > 0 ) // Even j (see Noll 1976)
631 : {
632 : Q2 = 2 * Q2 * pow( cos( m * phi ), 2 );
633 : }
634 : else if( m < 0 ) // Odd j (see Noll 1976)
635 : {
636 : Q2 = 2 * Q2 * pow( sin( -m * phi ), 2 );
637 : }
638 :
639 : return Q2;
640 : }
641 :
642 : extern template float zernikeQNorm<float>( float k, float phi, int n, int m );
643 :
644 : extern template double zernikeQNorm<double>( double k, double phi, int n, int m );
645 :
646 : extern template long double zernikeQNorm<long double>( long double k, long double phi, int n, int m );
647 :
648 : #ifdef HASQUAD
649 : extern template __float128 zernikeQNorm<__float128>( __float128 k, __float128 phi, int n, int m );
650 : #endif
651 :
652 : /// Calculate the square-normed Fourier transform of a Zernike polynomial at position (k,phi)
653 : /** Implements Equation (8) of Noll (1976) \cite noll_1976.
654 : *
655 : * \returns the value of |Q(k,phi)|^2
656 : *
657 : * \tparam realT is the floating point type used for arithmetic
658 : *
659 : */
660 : template <typename realT>
661 : realT zernikeQNorm( realT k, /**< [in] the radial coordinate of normalized spatial frequency. This is in the
662 : \cite noll_1976 convention of cycles-per-radius.*/
663 : realT phi, ///< [in] the azimuthal coordinate of normalized spatial frequency
664 : int j ///< [in] the Zernike polynomial index j (Noll convention)
665 : )
666 : {
667 : int n, m;
668 :
669 : noll_nm( n, m, j );
670 :
671 : return zernikeQNorm( k, phi, n, m );
672 : }
673 :
674 : /// Fill in an Eigen-like array with the square-normed Fourier transform of a Zernike polynomial
675 : /** The array is filled in with the values of |Q(k,phi)|^2 according to Equation (8) of Noll (1976) \cite noll_1976.
676 : *
677 : * \test Scenario: testing zernikeQNorm \ref tests_sigproc_zernike_zernikeQNorm "[test doc]"
678 : *
679 : * \returns 0 on success
680 : * \returns -1 on error
681 : *
682 : * \tparam arrayT is the Eigen-like array type. Arithmetic will be done in arrayT::Scalar.
683 : */
684 : template <typename arrayT>
685 1 : int zernikeQNorm( arrayT &arr, /**< [out] the allocated array. The rows() and cols() members are used to size
686 : the transform.*/
687 : arrayT &k, /**< [in] the normalized spatial frequency magnitude at each pixel. This is in the \cite
688 : noll_1976 convention of cycles-per-radius.*/
689 : arrayT &phi, ///< [in] the spatial frequency angle at each pixel
690 : int j ///< [in] the polynomial index in the Noll convention \cite noll_1976
691 : )
692 : {
693 1 : if( arr.rows() != k.rows() || arr.cols() != k.cols() )
694 : {
695 0 : internal::mxlib_error_report( error_t::invalidarg, "output array and input k are not the same size" );
696 0 : return -1;
697 : }
698 :
699 1 : if( arr.rows() != phi.rows() || arr.cols() != phi.cols() )
700 : {
701 0 : internal::mxlib_error_report( error_t::invalidarg, "output array and input phi are not the same size" );
702 0 : return -1;
703 : }
704 :
705 : int n, m;
706 1 : if( noll_nm( n, m, j ) < 0 )
707 0 : return -1; // noll_nm will explain error
708 :
709 33 : for( size_t i = 0; i < arr.rows(); ++i )
710 : {
711 1056 : for( size_t j = 0; j < arr.cols(); ++j )
712 : {
713 1024 : arr( i, j ) = zernikeQNorm( k( i, j ), phi( i, j ), n, m );
714 : }
715 : }
716 1 : return 0;
717 : }
718 :
719 : /// Calculate the spatial power spectrum of Piston
720 : template <typename realT>
721 : realT zernikePPiston( const realT &kD /**< [in] Spatial frequency in diameter units, i.e. cycles per aperture.*/ )
722 : {
723 : return 4 * pow( math::func::jinc( math::pi<realT>() * kD ), 2 );
724 : }
725 :
726 : /// Calculate the spatial power spectrum of Tip \& Tilt
727 : template <typename realT>
728 : realT zernikePTipTilt( const realT &kD /**< [in] Spatial frequency in diameter units, i.e. cycles per aperture.*/ )
729 : {
730 : return 16 * pow( math::func::jincN( 2, math::pi<realT>() * kD ), 2 );
731 : }
732 :
733 : /// Calculate the spatial power spectrum of Defocus
734 : template <typename realT>
735 : realT zernikePDefocus( const realT &kD /**< [in] Spatial frequency in diameter units, i.e. cycles per aperture.*/ )
736 : {
737 : return 12 * pow( math::func::jincN( 3, math::pi<realT>() * kD ), 2 );
738 : }
739 :
740 : /// Calculate the spatial power spectrum of Astigmatism
741 : template <typename realT>
742 : realT zernikePAstig( const realT &kD /**< [in] Spatial frequency in diameter units, i.e. cycles per aperture.*/ )
743 : {
744 : return 24 * pow( math::func::jincN( 3, math::pi<realT>() * kD ), 2 );
745 : }
746 :
747 : /// Calculate the spatial power spectrum of Coma
748 : template <typename realT>
749 : realT zernikePComa( const realT &kD /**< [in] Spatial frequency in diameter units, i.e. cycles per aperture.*/ )
750 : {
751 : return 32 * pow( math::func::jincN( 4, math::pi<realT>() * kD ), 2 );
752 : }
753 :
754 : /// Calculate the spatial power spectrum of Trefoil
755 : template <typename realT>
756 : realT zernikePTrefoil( const realT &kD /**< [in] Spatial frequency in diameter units, i.e. cycles per aperture.*/ )
757 : {
758 : return 32 * pow( math::func::jincN( 4, math::pi<realT>() * kD ), 2 );
759 : }
760 :
761 : ///@} signal_processing
762 :
763 : } // namespace sigproc
764 : } // namespace mx
765 :
766 : #endif // math_zernike_hpp
|
