Skip to content

Commit

Permalink
remove the duplicate zernike function
Browse files Browse the repository at this point in the history
  • Loading branch information
@rockfool
rockfool committed Sep 21, 2019
1 parent de0ad46 commit 88cceb1
Show file tree
Hide file tree
Showing 2 changed files with 3 additions and 200 deletions.
6 changes: 3 additions & 3 deletions src/material.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1515,7 +1515,7 @@ PolyProperty::evaluate_zernike1d(Position r) {
rho = std::sqrt(r.x * r.x + r.y * r.y) / coeffs_[0];
c_index = 2;
for( i=0; i <= order_; i=i+2){
poly_val = calc_zn_one_d(i,rho,0.0);
calc_zn_rad(i, rho, &poly_val);
property += poly_val * coeffs_[c_index-1];
c_index = c_index + 1;
}
Expand All @@ -1530,9 +1530,9 @@ PolyProperty::evaluate_zernike(Position r) {
double property {0.0};
//! Get normalized positions
rho = sqrt( r.x * r.x + r.y * r.y) / coeffs_[0];
phi = std::atan2(r.y,r.x);
phi = std::atan2(r.y, r.x);
k = 2; //! Keeps tracking of index in this % coeffs
calc_zn_old(order_, rho, phi, poly_norm_, poly_results_);
calc_zn(order_, rho, phi, poly_results_);
for( i=0; i<=order_; i++){
for( j=1; j<=i+1; j++){
property += poly_results_[k-1-1] * coeffs_[k-1];
Expand Down
197 changes: 0 additions & 197 deletions src/math_functions.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -887,202 +887,5 @@ std::complex<double> w_derivative(std::complex<double> z, int order)
}
}

// cvmt
double calc_zn_one_d(int n, double rho, double phi) {
double zn {0.0}; //! The resultant Z_n(uvw)
//! Right now this is only in 2D.
//! n == radial degree
//! m == azimuthal frequency
//!
switch(n){
case 0:
//! n = 0, m = 0
zn = ( ( 1.00 ) ) \
* ( 1.000000000000 );
break;
case 2:
//! n = 2, m = 0
zn = ( ( -1.00 ) + \
( 2.00 ) *rho * rho ) \
* ( 1.732050807569 );
break;
case 4:
zn = ( ( 1.00 ) + \
( -6.00 ) *rho * rho + \
( 6.00 ) *rho *rho *rho * rho ) \
* ( 2.236067977500 );
break;
case 6:
//! n = 6, m = 0
zn = ( ( -1.00 ) + \
( 12.00 ) *rho * rho + \
( -30.00 ) *rho *rho *rho * rho + \
( 20.00 ) *rho *rho *rho *rho *rho * rho ) \
* ( 2.645751311065 );
break;
case 8:
//! n = 8, m = 0
zn = ( ( 1.00 ) + \
( -20.00 ) *rho * rho + \
( 90.00 ) *rho *rho *rho * rho + \
( -140.00 ) *rho *rho *rho *rho *rho * rho + \
( 70.00 ) *rho *rho *rho *rho *rho *rho *rho * rho ) \
* ( 3.000000000000 );
break;
case 10:
//! n = 10, m = 0
zn = ( ( -1.00 ) + \
( 30.00 ) *rho * rho + \
( -210.00 ) *rho *rho *rho * rho + \
( 560.00 ) *rho *rho *rho *rho *rho * rho + \
( -630.00 ) *rho *rho *rho *rho *rho *rho *rho * rho + \
( 252.00 ) *rho *rho *rho *rho *rho *rho *rho *rho *rho * rho ) \
* ( 3.316624790355 );
break;
case 12:
//! n = 12, m = 0
zn = ( ( 1.00 ) + \
( -42.00 ) *rho * rho + \
( 420.00 ) *rho *rho *rho * rho + \
( -1680.00 ) *rho *rho *rho *rho *rho * rho + \
( 3150.00 ) *rho *rho *rho *rho *rho *rho *rho * rho + \
( -2772.00 ) *rho *rho *rho *rho *rho *rho *rho *rho *rho * rho + \
( 924.00 ) *rho *rho *rho *rho *rho *rho *rho *rho *rho *rho *rho * rho ) \
* ( 3.605551275464 );
break;
case 14:
//! n = 14, m = 0
zn = ( ( -1.00 ) + \
( 56.00 ) *rho * rho + \
( -756.00 ) *rho *rho *rho * rho + \
( 4200.00 ) *rho *rho *rho *rho *rho * rho + \
( -11550.00 ) *rho *rho *rho *rho *rho *rho *rho * rho + \
( 16632.00 ) *rho *rho *rho *rho *rho *rho *rho *rho *rho * rho + \
( -12012.00 ) *rho *rho *rho *rho *rho *rho *rho *rho *rho *rho *rho * rho + \
( 3432.00 ) *rho *rho *rho *rho *rho *rho *rho *rho *rho *rho *rho *rho *rho * rho ) \
* ( 3.872983346207 );
break;
case 16:
//! n = 16, m = 0
zn = ( ( 1.00 ) + \
( -72.00 ) *rho * rho + \
( 1260.00 ) *rho *rho *rho * rho + \
( -9240.00 ) *rho *rho *rho *rho *rho * rho + \
( 34650.00 ) *rho *rho *rho *rho *rho *rho *rho * rho + \
( -72072.00 ) *rho *rho *rho *rho *rho *rho *rho *rho *rho * rho + \
( 84084.00 ) *rho *rho *rho *rho *rho *rho *rho *rho *rho *rho *rho * rho + \
( -51480.00 ) *rho *rho *rho *rho *rho *rho *rho *rho *rho *rho *rho *rho *rho * rho + \
( 12870.00 ) *rho *rho *rho *rho *rho *rho *rho *rho *rho *rho *rho *rho *rho *rho *rho * rho ) \
* ( 4.123105625618 );
break;
case 18:
//! n = 18, m = 0
zn = ( ( -1.00 ) + \
( 90.00 ) *rho * rho + \
( -1980.00 ) *rho *rho *rho * rho + \
( 18480.00 ) *rho *rho *rho *rho *rho * rho + \
( -90090.00 ) *rho *rho *rho *rho *rho *rho *rho * rho + \
( 252252.00 ) *rho *rho *rho *rho *rho *rho *rho *rho *rho * rho + \
( -420420.00 ) *rho *rho *rho *rho *rho *rho *rho *rho *rho *rho *rho * rho + \
( 411840.00 ) *rho *rho *rho *rho *rho *rho *rho *rho *rho *rho *rho *rho *rho * rho + \
( -218790.00 ) *rho *rho *rho *rho *rho *rho *rho *rho *rho *rho *rho *rho *rho *rho *rho * rho + \
( 48620.00 ) *rho *rho *rho *rho *rho *rho *rho *rho *rho *rho *rho *rho *rho *rho *rho *rho *rho * rho ) \
* ( 4.358898943541 );
break;
default:
zn = 0.0;
break;
}
return zn;
}

void calc_zn_old(int n, double rho, double phi, double sqrt_norm[], double zn[]){
//! This efficient method for calculation R(m,n) is taken from
//! Chong, C. W., Raveendran, P., & Mukundan, R. (2003). A comparative
//! analysis of algorithms for fast computation of Zernike moments.
//! Pattern Recognition, 36(3), 731-742.
//!
//! n ! Maximum order
//! rho ! Radial location in the unit disk
//! phi ! Theta (radians) location in the unit disk
//! sqrt_norm(:)! Sqrt normalization for each order n
//! zn(:) ! The resulting list of coefficients
double sin_phi, cos_phi; //! Sine and Cosine of phi
std::vector<double> sin_phi_vec; //! Contains sin(n*phi)
std::vector<double> cos_phi_vec; //! Contains cos(n*phi)
std::vector<std::vector<double>> zn_mat;//! Matrix form of the coefficients which is
//! easier to work with
int i,p,q; //! Loop counters
double k1, k2, k3, k4; //! Variables for R_m_n calculation
//
//! n == radial degree
//! m == azimuthal frequency
sin_phi_vec.resize(n+1);
cos_phi_vec.resize(n+1);
zn_mat.resize(n+1);
for(i=0; i<=n+1;i++){
zn_mat[i].resize(n+1);
}
//! Deterine vector of sin(n*phi) and cos(n*phi)
sin_phi = std::sin(phi);
cos_phi = std::cos(phi);

sin_phi_vec[1-1] = 1.0;
cos_phi_vec[1-1] = 1.0;

sin_phi_vec[2-1] = 2.0 * cos_phi;
cos_phi_vec[2-1] = cos_phi;

for( i=3; i<=n+1; i++){
sin_phi_vec[i-1] = 2.0 * cos_phi * sin_phi_vec[i-1-1] - sin_phi_vec[i-2-1];
cos_phi_vec[i-1] = 2.0 * cos_phi * cos_phi_vec[i-1-1] - cos_phi_vec[i-2-1];
}

for( i=1; i<=n+1;i++){
sin_phi_vec[i-1] = sin_phi_vec[i-1] * sin_phi;
}

//! Calculate R_m_n(rho)
//! Fill the main diagonal first
for( p=0; p<=n; p++){
zn_mat[p+1-1][p+1-1] = std::pow(rho,p);
}
//! Fill in the second diagonal
for( q=0; q<=n-2; q++){
zn_mat[q+2+1-1][q+1-1] = (q+2) * zn_mat[q+2+1-1][q+2+1-1] - (q+1) * zn_mat[q+1-1][q+1-1];
}
//! Fill in the rest of the values using the original results
for( p=4; p<=n; p++){
k2 = 2 * p * (p - 1) * (p - 2);
for( q=p-4; q>=0; q=q-2){
k1 = (p + q) * (p - q) * (p - 2) / 2;
k3 = -pow(q,2)*(p - 1) - p * (p - 1) * (p - 2);
k4 = -p * (p + q - 2) * (p - q - 2) / 2;
zn_mat[p+1-1][q+1-1] = ((k2 * rho * rho + k3) * zn_mat[p-2+1-1][q+1-1]
+ k4 * zn_mat[p-4+1-1][q+1-1]) / k1;
}
}
//
//! Roll into a single vector for easier computation later
//! The vector is ordered (0,0), (1,-1), (1,1), (2,-2), (2,0),
//! (2, 2), .... in (n,m) indices
//! Note that the cos and sin vectors are offest by one
//! sin_phi_vec = [sin(x), sin(2x), sin(3x) ...]
//! cos_phi_vec = [1.0, cos(x), cos(2x)... ]
i = 1;
for( p=0; p<=n; p++){
for( q=-p; q<=p; q=q+2){
if (q < 0) {
zn[i-1] = zn_mat[p+1-1][abs(q)+1-1] * sin_phi_vec[p-1] * sqrt_norm[i-1];
} else if ( q == 0) {
zn[i-1] = zn_mat[p+1-1][q+1-1] * sqrt_norm[i-1];
} else {
zn[i-1] = zn_mat[p+1-1][q+1-1] * cos_phi_vec[p+1-1] * sqrt_norm[i-1];
}
i = i + 1;
}
}

}

} // namespace openmc

0 comments on commit 88cceb1

Please sign in to comment.