eddy/zerndeco
1
0
Fork 0
mirror of https://github.com/eddyem/zerndeco.git synced 2025-12-06 02:25:10 +03:00
Code Issues Packages Projects Releases Wiki Activity
zerndeco/zernike_annular.c
eddyem abee079d26 copy
2016-07-14 09:47:49 +03:00

389 lines
12 KiB
C
Raw Permalink Blame History

/*
* zernike_annular.c - annular Zernike polynomials
*
* Copyright 2013 Edward V. Emelianoff <eddy@sao.ru>
*
* This program is free software; you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation; either version 2 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program; if not, write to the Free Software
* Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston,
* MA 02110-1301, USA.
*/
/*
* These polynomials realized according to formulae in:
@ARTICLE{1981JOSA...71...75M,
author = {{Mahajan}, V.~N.},
title = "{Zernike annular polynomials for imaging systems with annular pupils.}",
journal = {Journal of the Optical Society of America (1917-1983)},
keywords = {Astronomical Optics},
year = 1981,
volume = 71,
pages = {75-85},
adsurl = {http://adsabs.harvard.edu/abs/1981JOSA...71...75M},
adsnote = {Provided by the SAO/NASA Astrophysics Data System}
}
@ARTICLE{2006ApOpt..45.3442H,
author = {{Hou}, X. and {Wu}, F. and {Yang}, L. and {Wu}, S. and {Chen}, Q.
},
title = "{Full-aperture wavefront reconstruction from annular subaperture interferometric data by use of Zernike annular polynomials and a matrix method for testing large aspheric surfaces}",
journal = {\ao},
keywords = {Mathematics, Optical data processing, Metrology, Optical testing},
year = 2006,
month = may,
volume = 45,
pages = {3442-3455},
doi = {10.1364/AO.45.003442},
adsurl = {http://adsabs.harvard.edu/abs/2006ApOpt..45.3442H},
adsnote = {Provided by the SAO/NASA Astrophysics Data System}
}
*/
#include "zernike.h"
#include "zern_private.h"
#include "usefull_macros.h"
// epsilon (==rmin)
double eps = -1., eps2 = -1.; // this global variable changes with each conv_r call!
/**
* convert coordinates of points on ring into normalized coordinates in [0, 1]
* for further calculations of annular Zernike (expression for R_k^0)
* @param r0 (i) - original coordinates array
* @param Sz - size of r0
* @return dinamically allocated array with new coordinates, which zero element is
* (0,0)
*/
polar *conv_r(polar *r0, int Sz){
int x;
polar *Pn = MALLOC(polar, Sz);
double rmax = -1., rmin = 2.;
polar *p = r0, *p1 = Pn;
for(x = 0; x < Sz; x++, p++){
if(p->r < rmin) rmin = p->r;
if(p->r > rmax) rmax = p->r;
}
if(rmin > rmax || fabs(rmin - rmax) <= DBL_EPSILON || rmax > 1. || rmin < 0.)
ERRX(_("polar coordinates should be in [0,1]!"));
eps = rmin;
eps2 = eps*eps;
p = r0;
//rmax = 1. - eps2; // 1 - eps^2 -- denominator
rmax = rmax*rmax - eps2;
for(x = 0; x < Sz; x++, p++, p1++){
p1->theta = p->theta;
p1->r = sqrt((p->r*p->r - eps2) / rmax); // r = (r^2 - eps^2) / (1 - eps^2)
}
return Pn;
}
/**
* Allocate 2D array (H arrays of length W)
* @param len - array's width (number of elements in each 1D array
* @param num - array's height (number of 1D arrays)
* @return dinamically allocated 2D array filled with zeros
*/
_2D *_2ALLOC(int len, int num){
_2D *r = MALLOC(_2D, 1);
r->data = MALLOC(double*, num);
int x;
for(x = 0; x < num; x++)
r->data[x] = MALLOC(double, len);
r->len = (size_t)len;
r->num = (size_t)num;
return r;
}
void _2FREE(_2D **a){
int x, W = (*a)->num;
for(x = 0; x < W; x++) free((*a)->data[x]);
free((*a)->data);
FREE(*a);
}
/**
* Zernike radial function R_n^0
* @param n - order of polynomial
* @param Sz - number of points
* @param P (i) - modified polar coordinates: sqrt((r^2-eps^2)/(1-eps^2))
* @param R (o) - array with size Sz to which polynomial will be stored
*/
void Rn0(int n, int Sz, polar *P, double *R){
if(Sz < 1 || !P)
ERRX(_("Size of matrix must be > 0!"));
if(n < 0)
ERRX(_("Order of Zernike polynomial must be > 0!"));
if(!FK) build_factorial();
int j, k, iup = n/2;
double **Rpow = build_rpowR(n, Sz, P);
// now fill output matrix
double *Zptr = R;
polar *p = P;
for(j = 0; j < Sz; j++, p++, Zptr++){
double Z = 0.;
if(p->r > 1.) continue; // throw out points with R>1
// calculate R_n^m
for(k = 0; k <= iup; k++){ // Sum
double F = FK[iup-k];
double z = (1. - 2. * (k % 2)) * FK[n - k] // (-1)^k * (n-k)!
/(//------------------------------------- ---------------------
FK[k]*F*F // k!((n/2-k)!)^2
);
Z += z * Rpow[n-2*k][j]; // *R^{n-2k}
}
*Zptr = Z;
//DBG("order %d, point (%g, %g): %g\n", n, p->r, p->theta, Z);
}
// free unneeded memory
free_rpow(&Rpow, n);
}
/**
* Convert parameters n, m into index in Noll notation, p
* @param n, m - Zernike orders
* @return order in Noll notation
*/
int p_from_nm(int n, int m){
return (n*n + 2*n + m)/2;
}
/**
* Create array [p][pt] of annular Zernike polynomials
* @param pmax - max number of polynomial in Noll notation
* @param Sz - size of points array
* @param P (i) - array with points' coordinates
* @param Norm(o)- array [0..pmax] of sum^2 for each polynomial (dynamically allocated) or NULL
* @return
*/
//#undef DBG
//#define DBG(...)
_2D *ann_Z(int pmax, int Sz, polar *P, double **Norm){
int m, n, j, nmax, mmax, jmax, mW, jm;
convert_Zidx(pmax, &nmax, &mmax);
mmax = nmax; // for a case when p doesn't include a whole n powers
jmax = nmax / 2;
mW = mmax + 1; // width of row
jm = jmax + mW; // height of Q0 and h
polar *Pn = conv_r(P, Sz);
if(eps < 0. || eps2 < 0. || eps > 1. || eps2 > 1.)
ERRX(_("Wrong epsilon value!!!"));
// !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
_2D *h = _2ALLOC(mW, jm); // constants -> array [(jmax+1)+((mmax+1)-1)][mmax+1]
_2D *Q0= _2ALLOC(mW, jm); // Q(0) -> array [(jmax+1)+((mmax+1)-1)][mmax+1]
_2D *Q = _2ALLOC(Sz, (jmax + 1)*mW); // functions -> array [(jmax+1)*(mmax+1)][Sz]
DBG("Allocated Q&h, size of Q: %zdx%zd, size of H: %zdx%zd\n", Q->num, Q->len, h->num, h->len);
DBG("pmax = %d, nmax = mmax = %d, jmax = %d, jm = %d\n", pmax, nmax, jmax, jm);
DBG("eps=%g, eps2=%g\n", eps, eps2);
// in Q index of [j][m][point] is [j*mmax+m][point]
double **hd = h->data, **Qd = Q->data, **Q0d = Q0->data;
// fill h_j^0 & Q_j^0
for(j = 0; j < jm; j++){ // main matrix + triangle
DBG("j = %d, m = 0\n", j);
polar Zero = {0., 0.};
// h_j^0 = (1-eps^2)/(2(2j+1))
hd[j][0] = (1.-eps2)/(2.*j+1.)/2.;
// Q_j^0(rho, eps) = R_{2j}^0(r)
DBG("Q,h [%d][0], h=%g\n", j, hd[j][0]);
if(j <= jmax) Rn0(2*j, Sz, Pn, Qd[mW*j]);
Rn0(2*j, 1, &Zero, &Q0d[j][0]); // fill also Q(0)
}
FREE(Pn);
// fill h_j^m & Q_j&m
int JMAX = jm-1; // max number of j
for(m = 1; m <= mmax; m++, JMAX--){ // m in outern cycle because of Q_{j+1}!
DBG("\n\nm=%d\n",m);
int mminusone = m - 1;
int idx = mminusone;// index Q_j^{m-1}
for(j = 0; j < JMAX; j++, idx+=mW){
DBG("\nj=%d\n",j);
// 2(2j+2m-1) * h_j^{m-1}
// H_j^m = ------------------------------
// (j+m)(1-eps^2) * Q_j^{m-1}(0)
double H = 2*(2*(j+m)-1)/(j+m)/(1-eps2) * hd[j][mminusone]/Q0d[j][mminusone];
// h_j^m = -H_j^m * Q_{j+1}^{m-1}(0)
hd[j][m] = -H * Q0d[j+1][mminusone];
DBG("hd[%d][%d] = %g (H = %g)\n", j, m, hd[j][m], H);
// j Q_i^{m-1}(0) * Q_i^{m-1}(r)
// Q_j^m(r) = H_j^m * Sum -----------------------------
// i=0 h_i^{m-1}
if(j <= jmax){
int pt;
for(pt = 0; pt < Sz; pt++){ // cycle by points
int i, idx1 = mminusone;
double S = 0;
for(i = 0; i <=j; i++, idx1+=mW){ // cycle Sum
S += Q0d[i][mminusone] * Qd[idx1][pt] / hd[i][mminusone];
}
Qd[idx+1][pt] = H * S; // Q_j^m(rho^2)
//DBG("Q[%d][%d](%d) = %g\n", j, m, pt, H * S);
}
}
// and Q(0):
double S = 0.;
int i;
for(i = 0; i <=j; i++){
double qim = Q0d[i][mminusone];
S += qim*qim / hd[i][mminusone];
}
Q0d[j][m] = H * S; // Q_j^m(0)
DBG("Q[%d][%d](0) = %g\n", j, m, H * S);
}
}
_2FREE(&Q0);
// allocate memory for Z_p
_2D *Zarr = _2ALLOC(Sz, pmax + 1); // pmax is max number _including_
double **Zd = Zarr->data;
// now we should calculate R_n^m
DBG("R_{2j}^0\n");
// R_{2j}^0(rho) = Q_j^0(rho^2)
size_t S = Sz*sizeof(double);
for(j = 0; j <= jmax; j++){
int pidx = p_from_nm(2*j, 0);
if(pidx > pmax) break;
DBG("R[%d]\n", pidx);
memcpy(Zd[pidx], Qd[j*mW], S);
}
DBG("R_n^n\n");
// R_n^n(rho) = rho^n / sqrt(Sum_{i=0}^n[eps^{2i}])
double **Rpow = build_rpowR(nmax, Sz, P); // powers of rho
double epssum = 1., epspow = 1.;
for(n = 1; n <= nmax; n++){
int pidx = p_from_nm(n, -n); // R_n^{-n}
if(pidx > pmax) break;
DBG("R[%d] (%d, %d)\n", pidx, n, -n);
epspow *= eps2; // eps^{2n}
epssum += epspow;// sum_0^n eps^{2n}
double Sq = sqrt(epssum); // sqrt(sum...)
double *rptr = Zd[pidx], *pptr = Rpow[n];
int pt;
for(pt = 0; pt < Sz; pt++, rptr++, pptr++)
*rptr = *pptr / Sq; // rho^n / sqrt(...)
int pidx1 = p_from_nm(n, n); // R_n^{n}
if(pidx1 > pmax) break;
DBG("R[%d] (%d, %d)\n", pidx1, n, n);
memcpy(Zd[pidx1], Zd[pidx], S);
}
DBG("R_n^m\n");
/* / 1 - eps^2 \
R_{2j+|m|}^{|m|}(rho, eps) = sqrt| ------------------ | * rho^m * Q_j^m(rho^2)
\ 2(2j+|m|+1)*h_j^m / */
for(j = 1; j <= jmax; j++){
for(m = 1; m <= mmax; m++){
int N = 2*j + m, pt, pidx = p_from_nm(N, -m);
if(pidx > pmax) continue;
DBG("R[%d]\n", pidx);
double *rptr = Zd[pidx], *pptr = Rpow[m], *qptr = Qd[j*mW+m];
double hjm = hd[j][m];
for(pt = 0; pt < Sz; pt++, rptr++, pptr++, qptr++)
*rptr = sqrt((1-eps2)/2/(N+1)/hjm) * (*pptr) * (*qptr);
int pidx1 = p_from_nm(N, m);
if(pidx1 > pmax) continue;
DBG("R[%d]\n", pidx1);
memcpy(Zd[pidx1], Zd[pidx], S);
}
}
free_rpow(&Rpow, nmax);
_2FREE(&h);
_2FREE(&Q);
// last step: fill Z_n^m and find sum^2
int p;
double *sum2 = MALLOC(double, pmax+1);
for(p = 0; p <= pmax; p++){
convert_Zidx(p, &n, &m);
DBG("p=%d (n=%d, m=%d)\n", p, n, m);
int pt;
double (*fn)(double val);
double empty(double _U_ val){return 1.;}
if(m == 0) fn = empty;
else{if(m > 0)
fn = cos;
else
fn = sin;}
double m_abs = fabs((double)m), *rptr = Zd[p], *sptr = &sum2[p];
polar *Pptr = P;
for(pt = 0; pt < Sz; pt++, rptr++, Pptr++){
*rptr *= fn(m_abs * Pptr->theta);
//DBG("\tpt %d, val = %g\n", pt, *rptr);
*sptr += (*rptr) * (*rptr);
}
DBG("SUM p^2 = %g\n", *sptr);
}
if(Norm) *Norm = sum2;
else FREE(sum2);
return Zarr;
}
//#undef DBG
//#define DBG(...) do{fprintf(stderr, __VA_ARGS__); }while(0)
/**
* Compose wavefront by koefficients of annular Zernike polynomials
* @params like for ZcomposeR
* @return composed wavefront (dynamically allocated)
*/
double *ann_Zcompose(int Zsz, double *Zidxs, int Sz, polar *P){
int i;
double *image = MALLOC(double, Sz);
_2D *Zannular = ann_Z(Zsz-1, Sz, P, NULL);
double **Zcoeff = Zannular->data;
for(i = 0; i < Zsz; i++){ // now we fill array
double K = Zidxs[i];
if(fabs(K) < Z_prec) continue;
int j;
double *iptr = image, *zptr = Zcoeff[i];
for(j = 0; j < Sz; j++, iptr++, zptr++){
*iptr += K * (*zptr); // add next Zernike polynomial
}
}
_2FREE(&Zannular);
return image;
}
/**
* Decompose wavefront by annular Zernike
* @params like ZdecomposeR
* @return array of Zernike coefficients
*/
double *ann_Zdecompose(int Nmax, int Sz, polar *P, double *heights, int *Zsz, int *lastIdx){
int i, pmax, maxIdx = 0;
double *Zidxs = NULL, *icopy = NULL;
pmax = (Nmax + 1) * (Nmax + 2) / 2; // max Z number in Noll notation
Zidxs = MALLOC(double, pmax);
icopy = MALLOC(double, Sz);
memcpy(icopy, heights, Sz*sizeof(double)); // make image copy to leave it unchanged
*Zsz = pmax;
double *norm;
_2D *Zannular = ann_Z(pmax-1, Sz, P, &norm);
double **Zcoeff = Zannular->data;
for(i = 0; i < pmax; i++){ // now we fill array
int j;
double *iptr = icopy, *zptr = Zcoeff[i], K = 0., norm_i = norm[i];
for(j = 0; j < Sz; j++, iptr++, zptr++)
K += (*zptr) * (*iptr) / norm_i; // multiply matrixes to get coefficient
if(fabs(K) < Z_prec)
continue; // there's no need to substract values that are less than our precision
Zidxs[i] = K;
maxIdx = i;
iptr = icopy; zptr = Zcoeff[i];
for(j = 0; j < Sz; j++, iptr++, zptr++)
*iptr -= K * (*zptr); // subtract composed coefficient to reduce high orders values
}
if(lastIdx) *lastIdx = maxIdx;
_2FREE(&Zannular);
FREE(norm);
FREE(icopy);
return Zidxs;
}
Reference in New Issue View Git Blame Copy Permalink