diff --git a/DFTFringe.pro b/DFTFringe.pro index 2b760168..0271d3fc 100644 --- a/DFTFringe.pro +++ b/DFTFringe.pro @@ -263,6 +263,7 @@ SOURCES += SingleApplication/singleapplication.cpp \ zapm.cpp \ zernikedlg.cpp \ zernikeeditdlg.cpp \ + zernikepolar.cpp \ zernikeprocess.cpp \ zernikes.cpp \ zernikesmoothingdlg.cpp @@ -390,6 +391,7 @@ HEADERS += bezier/bezier.h \ wftstats.h \ zernikedlg.h \ zernikeeditdlg.h \ + zernikepolar.h \ zernikeprocess.h \ zernikes.h \ zernikesmoothingdlg.h diff --git a/DFTFringe_Dale.pro b/DFTFringe_Dale.pro index a35afefb..de04bf23 100644 --- a/DFTFringe_Dale.pro +++ b/DFTFringe_Dale.pro @@ -61,6 +61,7 @@ SOURCES += main.cpp \ surfacemanager.cpp \ zapm.cpp \ zernikedlg.cpp \ + zernikepolar.cpp \ zernikeprocess.cpp \ mirrordlg.cpp \ zernikes.cpp \ @@ -183,6 +184,7 @@ HEADERS += mainwindow.h \ surfaceanalysistools.h \ surfacemanager.h \ zernikedlg.h \ + zernikepolar.h \ zernikeprocess.h \ mirrordlg.h \ zernikes.h \ diff --git a/dftarea.cpp b/dftarea.cpp index f8895b45..757ffe1d 100644 --- a/dftarea.cpp +++ b/dftarea.cpp @@ -1471,13 +1471,11 @@ qDebug() << "rec" << left << top << width << height; qDebug() << "Rho " << r0 << "Theta" << t0; - - zernikePolar &zpolar = *zernikePolar::get_Instance(); - zpolar.init(r0,t0); + zernikePolar zpolar(r0,t0,100); // debug print the difference // for (int i = 0; i < 100; ++i){ -// qDebug() << "psi Zernike" <setZernTitle(ztitle); double z8 = zernTablemodel->values[8]; + double BestSC; if (m_mirrorDlg->doNull && wf.useSANull){ BestSC = z8/m_mirrorDlg->z8; } diff --git a/surfacemanager.cpp b/surfacemanager.cpp index 5b3dcec8..bd01eeb2 100644 --- a/surfacemanager.cpp +++ b/surfacemanager.cpp @@ -42,6 +42,7 @@ #include #include #include "zernikeprocess.h" +#include "zernikepolar.h" #include #include #include "rotationdlg.h" @@ -434,7 +435,6 @@ cv::Mat SurfaceManager::computeWaveFrontFromZernikes(int wx, int wy, std::vector std::vector &en = zernEnables; mirrorDlg *md = mirrorDlg::get_Instance(); - zernikePolar &zpolar = *zernikePolar::get_Instance(); for (int i = 0; i < wx; ++i) { double x1 = (double)(i - (xcen)) / rad; @@ -447,19 +447,19 @@ cv::Mat SurfaceManager::computeWaveFrontFromZernikes(int wx, int wy, std::vector { double S1 = 0; double theta = atan2(y1,x1); - zpolar.init(rho,theta); + zernikePolar zpolar(rho, theta, zernsToUse.size()); for (int ii = 0; ii < zernsToUse.size(); ++ii) { int z = zernsToUse[ii]; if ( z == 3 && m_surfaceTools->m_useDefocus){ - S1 += m_surfaceTools->m_defocus * zpolar.zernike(z,rho,theta); + S1 += m_surfaceTools->m_defocus * zpolar.zernike(z); } else { if (en[z]){ if (z == 8 && md->doNull) - S1 += md->z8 * zpolar.zernike(z,rho, theta); + S1 += md->z8 * zpolar.zernike(z); - S1 += zerns[z] * zpolar.zernike(z,rho, theta); + S1 += zerns[z] * zpolar.zernike(z); } } } @@ -744,7 +744,6 @@ void SurfaceManager::useDemoWaveFront(){ double rho; mirrorDlg *md = mirrorDlg::get_Instance(); cv::Mat result = cv::Mat::zeros(wx,wx, numType); - zernikePolar &zpolar = *zernikePolar::get_Instance(); for (int i = 0; i < wx; ++i) { double x1 = (double)(i - (xcen)) / rad; @@ -753,11 +752,11 @@ void SurfaceManager::useDemoWaveFront(){ double y1 = (double)(j - (ycen )) /rad; rho = sqrt(x1 * x1 + y1 * y1); double theta = atan2(y1,x1); - zpolar.init(rho,theta); + zernikePolar zpolar(rho, theta, 10); if (rho <= 1.) { - double S1 = md->z8 * -.9 * zpolar.zernike(8,rho,theta) + .02* zpolar.zernike(9, rho,theta); + double S1 = md->z8 * -.9 * zpolar.zernike(8) + .02* zpolar.zernike(9); result.at(j,i) = S1; } diff --git a/zernikeeditdlg.cpp b/zernikeeditdlg.cpp index f1a56c41..21e4c3c8 100644 --- a/zernikeeditdlg.cpp +++ b/zernikeeditdlg.cpp @@ -3,6 +3,7 @@ #include #include #include "zernikeprocess.h" +#include "zernikepolar.h" #include "mirrordlg.h" #include "myutils.h" zernikeEditDlg::zernikeEditDlg(SurfaceManager * sfm, QWidget *parent) : @@ -53,7 +54,6 @@ void zernikeEditDlg::on_createSurface_clicked() double xcen = (size -1)/2.; double ycen = xcen; double rad = xcen - 1; - zernikePolar &zpolar = *zernikePolar::get_Instance(); for (int y = 0; y < size; ++y) { double uy = (double)(y - (ycen)) / rad; @@ -65,12 +65,12 @@ void zernikeEditDlg::on_createSurface_clicked() if (rho <= 1.) { double theta = atan2(uy,ux); - zpolar.init(rho,theta); + zernikePolar zpolar(rho, theta, tableModel->rowCount()); double s1 = 0; for (int z = 0; z < tableModel->rowCount(); ++z){ double v = tableModel->values[z]; if (m_zernEnables[z]){ - s1 += v * zpolar.zernike(z, rho, theta); + s1 += v * zpolar.zernike(z); } } result.at(y,x) = s1; diff --git a/zernikepolar.cpp b/zernikepolar.cpp new file mode 100644 index 00000000..633484e0 --- /dev/null +++ b/zernikepolar.cpp @@ -0,0 +1,139 @@ +/****************************************************************************** +** +** Copyright 2016 Dale Eason +** This file is part of DFTFringe +** 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 version 3 of the License + +** DFTFringe 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 DFTFringe. If not, see . + +****************************************************************************/ + +#include "zernikepolar.h" +#include +#include + +zernikePolar::zernikePolar(double rho, double theta, size_t nbTerms) { + // Having all terms computed at once here let's compiler optimize the code better + + if(nbTerms > 49){ + qWarning() << "zernikePolar::init: maxorder is limited to 49, setting to 49"; + nbTerms = 49; + } + + m_nbTermsComputed = nbTerms; + + zernTerms[0] = 1.; + + const double rho2 = rho * rho; + const double costheta = std::cos(theta); + const double sintheta = std::sin(theta); + const double cos2theta = std::cos(2. * theta); + const double sin2theta = std::sin(2. * theta); + zernTerms[1] = rho * costheta; + zernTerms[2] = rho * sintheta; + zernTerms[3] = -1. + 2. * rho2; + zernTerms[4] = rho2 * cos2theta; + zernTerms[5] = rho2 * sin2theta; + zernTerms[6] = rho * (-2. + 3. * rho2) * costheta; + zernTerms[7] = rho * (-2. + 3. * rho2) * sintheta; + zernTerms[8] = 1. + rho2 * (-6 + 6. * rho2); + + double rho3; + double rho4; + double rho5; + double rho6; + double rho8; + double cos3theta; + double sin3theta; + double cos4theta; + double sin4theta; + double rho10; + double cos5theta; + double sin5theta; + + // only compute what is actually needed + // but to avoid complex code I use only 4 ranges + if(nbTerms > 8) + { + rho3 = rho2 * rho; + rho4 = rho3 * rho; + rho5 = rho4 * rho; + rho6 = rho5 * rho; + rho8 = rho6 * rho2; + cos3theta = std::cos(3. * theta); + sin3theta = std::sin(3. * theta); + cos4theta = std::cos(4. * theta); + sin4theta = std::sin(4. * theta); + + zernTerms[9] = rho3 * cos3theta; + zernTerms[10] = rho3 * sin3theta; + zernTerms[11] = rho2 * (-3 + 4 * rho2) * cos2theta; + zernTerms[12] = rho2 * (-3 + 4 * rho2) * sin2theta ; + zernTerms[13] = rho * (3. - 12. * rho2 + 10. * rho4) * costheta; + zernTerms[14] = rho * (3. - 12. * rho2 + 10. * rho4) * sintheta; + zernTerms[15] = -1 + 12 * rho2 - 30. * rho4 + 20. * rho6; + zernTerms[16] = rho4 * cos4theta; + zernTerms[17] = rho4 * sin4theta; + zernTerms[18] = rho3 *( -4. + 5. * rho2) * cos3theta; + zernTerms[19] = rho3 *( -4. + 5. * rho2) * sin3theta; + zernTerms[20] = rho2 * (6. - 20. * rho2 + 15 * rho4)* cos2theta; + zernTerms[21] = rho2 * (6. - 20. * rho2 + 15 * rho4)* sin2theta; + zernTerms[22] = rho * (-4. + 30. * rho2 - 60. * rho4 + 35 * rho6)* costheta; + zernTerms[23] = rho * (-4. + 30. * rho2 - 60. * rho4 + 35 * rho6)* sintheta; + zernTerms[24] = 1. - 20. * rho2 + 90. * rho4 - 140. * rho6 + 70. * rho8; + } + + if(nbTerms > 24) { + rho10 = rho8 * rho2; + cos5theta = std::cos(5. * theta); + sin5theta = std::sin(5. * theta); + + zernTerms[25] = rho5 * cos5theta; + zernTerms[26] = rho5 * sin5theta; + zernTerms[27] = rho4 * (-5. + 6. * rho2) * cos4theta; + zernTerms[28] = rho4 * (-5. + 6. * rho2) * sin4theta; + zernTerms[29] = rho3 * (10. - 30. * rho2 + 21. * rho4) * cos3theta; + zernTerms[30] = rho3 * (10. - 30. * rho2 + 21. * rho4) * sin3theta; + zernTerms[31] = rho2 *(-10. + 60. * rho2 - 105. * rho4 + 56. * rho6) * cos2theta; + zernTerms[32] = rho2 *(-10. + 60. * rho2 - 105. * rho4 + 56. * rho6) * sin2theta; + zernTerms[33] = rho * (5. - 60. * rho2 + 210 * rho4 -280. * rho6 + 126. * rho8) * costheta; + zernTerms[34] = rho * (5. - 60. * rho2 + 210 * rho4 -280. * rho6 + 126. * rho8) * sintheta; + zernTerms[35] = -1 + 30. * rho2 -210 * rho4 + 560. * rho6 - 630 * rho8 + 252. * rho10; + } + + if(nbTerms > 35) + { + zernTerms[36] = rho6 * std::cos(6. * theta); + zernTerms[37] = rho6 * std::sin(6. * theta); + zernTerms[38] = rho5 * (-6. + 7 * rho2) * cos5theta; + zernTerms[39] = rho5 * (-6. + 7 * rho2) * sin5theta; + zernTerms[40] = rho4 * (15. -42. * rho2 + 28. * rho4) * cos4theta; + zernTerms[41] = rho4 * (15. -42. * rho2 + 28. * rho4) * sin4theta; + zernTerms[42] = rho3 * (-20 + 105. * rho2 - 168. * rho4 + 84 * rho6) * cos3theta; + zernTerms[43] = rho3 * (-20. + 105. * rho2 - 168. * rho4 + 84. * rho6) * sin3theta; + zernTerms[44] = rho2 * (15. - 140. * rho2 + 420. * rho4 - 504. * rho6 + 210. * rho8) * cos2theta; + zernTerms[45] = rho2 * (15. - 140. * rho2 + 420. * rho4 - 504. * rho6 + 210. * rho8) * sin2theta; + zernTerms[46] = rho *(-6. + 105 * rho2 - 560. * rho4 + 1260. * rho6 -1260. * rho8 +462. * rho10) * costheta; + zernTerms[47] = rho *(-6. + 105 * rho2 - 560. * rho4 + 1260. * rho6 -1260. * rho8 +462. * rho10) * sintheta; + zernTerms[48] = 1. - 42. * rho2 + 420. * rho4 - 1680. * rho6 + 3150. * rho8 -2772. * rho10 + 924. * rho10 * rho2; + } +} + +double zernikePolar::zernike(size_t n){ + if(n < m_nbTermsComputed) { + return zernTerms[n]; + } + else + { + throw std::out_of_range("Zernike order exceeds maximum computed order"); + return 0.; + } +} \ No newline at end of file diff --git a/zernikepolar.h b/zernikepolar.h new file mode 100644 index 00000000..c832c03e --- /dev/null +++ b/zernikepolar.h @@ -0,0 +1,35 @@ +/****************************************************************************** +** +** Copyright 2016 Dale Eason +** This file is part of DFTFringe +** 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 version 3 of the License + +** DFTFringe 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 DFTFringe. If not, see . + +****************************************************************************/ + +#ifndef ZERNIKEPOLAR_H +#define ZERNIKEPOLAR_H + +#include + +class zernikePolar +{ +public: + zernikePolar(double rho, double theta, size_t nbTerms = 48); + double zernike(size_t z); +private: + size_t m_nbTermsComputed = 0; + double zernTerms[49]; +}; + + +#endif \ No newline at end of file diff --git a/zernikeprocess.cpp b/zernikeprocess.cpp index 2e885cf6..538cae5a 100644 --- a/zernikeprocess.cpp +++ b/zernikeprocess.cpp @@ -16,6 +16,7 @@ ****************************************************************************/ #include "zernikeprocess.h" +#include "zernikepolar.h" #include #include #include "mainwindow.h" @@ -26,19 +27,10 @@ #include "myutils.h" //#include "arbitrarywavefrontdlg.h" #include "userdrawnprofiledlg.h" + + std::vector zernEnables; std::vector zNulls; -double BestSC = -1.; -int Zw[] = { /* n */ - 1, // 0 - 4,4,3, // 1 - 6,6,8,8,5, //8 // 2 - 8,8,10,10,12,12,7, // 15 3 - 10,10,12,12,14,14,16,16,9, //24 4 - 12,12,14,14,16,16,18,18,20,20,11, //35 5 - 14,14,16,16,18,18,20,20,22,22,24,24,13 // 48 -}; - // @@ -131,15 +123,7 @@ cv::Mat zpmCx(QVector rho, QVector theta, int maxorder) { } return zm; } -/* -Public Function ZernikeW(n As Integer) As Double -' N is Zernike number as given by James Wyant -' routine calculates the inverse of the weight in the variance computation -*/ -int ZernikeW(int n) -{ - return(Zw[n]); -} + long fact( int n1, int n2) { @@ -151,9 +135,8 @@ long fact( int n1, int n2) return f; } -void gauss_jordan(int n, double* Am, double* Bm) -{ - /* + +/* Private Sub GJ(n As Integer) Dim indxc(50) As Integer, indxr(50) As Integer, ipiv(50) As Integer Dim i As Integer, icol As Integer, irow As Integer, j As Integer @@ -167,6 +150,10 @@ Dim big As Double, dum As Double, pivinv As Double, temp As Double 'on entry, Am(1 to N, 1 to N) is the matrix coefficients 'on exit, Am is the inverse matrix */ +/* +void gauss_jordan(int n, double* Am, double* Bm) +{ + int* ipiv = new int[n]; int* indxr = new int[n]; int* indxc = new int[n]; @@ -260,166 +247,15 @@ Dim big As Double, dum As Double, pivinv As Double, temp As Double delete[] ipiv; delete[] indxr; delete[] indxc; -} +}*/ + + /* Public Function Zernike(n As Integer, X As Double, Y As Double) As Double ' N is Zernike number as given by James Wyant */ -zernikePolar *zernikePolar::m_instance = 0; -zernikePolar *zernikePolar::get_Instance(){ - if (m_instance == 0){ - m_instance = new zernikePolar; - } - return m_instance; -} - -void zernikePolar::init(double rho, double theta){ - rho2 = rho * rho; - rho3 = pow(rho,3.); - rho4 = pow(rho,4.); - rho5 = pow(rho,5.); - rho6 = pow(rho,6.); - rho8 = pow(rho,8.); - rho10 = pow(rho,10.); - costheta = cos(theta); - sintheta = sin(theta); - cos2theta = cos(2. * theta); - sin2theta = sin(2. * theta); - cos3theta = cos(3. * theta); - sin3theta = sin(3. * theta); - cos4theta = cos(4. * theta); - sin4theta = sin(4. * theta); - cos5theta = cos(5. * theta); - sin5theta = sin(5. * theta); - } - -double zernikePolar::zernike(int n, double rho, double theta){ - - switch(n){ - case 0: return 1.; - break; - case 1: return rho * costheta; - break; - case 2: return rho * sintheta; - break; - case 3: return -1. + 2. * rho * rho; - break; - case 4: return rho2 * cos2theta; - break; - case 5: return rho2 * sin2theta; - break; - case 6: return rho * (-2. + 3. * rho2) * costheta; - break; - case 7: return rho * (-2. + 3. * rho2) * sintheta; - break; - case 8: - { - return 1. + rho2 * (-6 + 6. * rho2); - break; - } - case 9: return rho * rho * rho * cos3theta; - break; - case 10: return rho * rho * rho * sin3theta; - break; - - case 11: - - return rho2 * (-3 + 4 * rho2) * cos2theta; - - break; - case 12: - - return rho2 * (-3 + 4 * rho2) * sin2theta ; - - break; - - case 13: - return rho * (3. - 12. * rho2 + 10. * rho4) * costheta; - break; - case 14: - return rho * (3. - 12. * rho2 + 10. * rho4) * sintheta; - break; - case 15: - return -1 + 12 * rho2 - 30. * rho4 + 20. * rho6; - break; - case 16: return rho4 * cos4theta; - break; - case 17: return rho4 * sin4theta; - break; - case 18: return rho3 *( -4. + 5. * rho2) * cos3theta; - break; - case 19: return rho3 *( -4. + 5. * rho2) * sin3theta; - break; - case 20: return rho2 * (6. - 20. * rho2 + 15 * rho4)* cos2theta; - break; - case 21: return rho2 * (6. - 20. * rho2 + 15 * rho4)* sin2theta; - break; - case 22: return rho * (-4. + 30. * rho2 - 60. * rho4 + 35 * rho6)* costheta; - break; - case 23: return rho * (-4. + 30. * rho2 - 60. * rho4 + 35 * rho6)* sintheta; - break; - case 24: return 1. - 20. * rho2 + 90. * rho4 - 140. * rho6 + 70. * rho8; - break; - case 25: return rho5 * cos5theta; - break; - case 26: return rho5 * sin5theta; - break; - case 27: return rho4 * (-5. + 6. * rho2) * cos4theta; - break; - case 28: return rho4 * (-5. + 6. * rho2) * sin4theta; - break; - case 29: return rho3 * (10. - 30. * rho2 + 21. * rho4) * cos3theta; - break; - case 30: return rho3 * (10. - 30. * rho2 + 21. * rho4) * sin3theta; - break; - case 31: return rho2 *(-10. + 60. * rho2 - 105. * rho4 + 56. * rho6) * cos2theta; - break; - case 32: return rho2 *(-10. + 60. * rho2 - 105. * rho4 + 56. * rho6) * sin2theta; - break; - case 33: return rho * (5. - 60. * rho2 + 210 * rho4 -280. * rho6 + 126. * rho8) * costheta; - break; - case 34: return rho * (5. - 60. * rho2 + 210 * rho4 -280. * rho6 + 126. * rho8) * sintheta; - break; - case 35: return -1 + 30. * rho2 -210 * rho4 + 560. * rho6 - 630 * rho8 + 252. * rho10; - break; - case 36: return rho6 * cos(6. * theta); - break; - case 37: return rho6 * sin(6. * theta); - break; - case 38: return rho5 * (-6. + 7 * rho2) * cos5theta; - break; - case 39: return rho5 * (-6. + 7 * rho2) * sin5theta; - break; - case 40: return rho4 * (15. -42. * rho2 + 28. * rho4) * cos4theta; - break; - case 41: return rho4 * (15. -42. * rho2 + 28. * rho4) * sin4theta; - break; - case 42: return rho3 * (-20 + 105. * rho2 - 168. * rho4 + 84 * rho6) * cos3theta; - break; - case 43: return rho3 * (-20. + 105. * rho2 - 168. * rho4 + 84. * rho6) * sin3theta; - break; - case 44: return rho2 * (15. - 140. * rho2 + 420. * rho4 - 504. * rho6 + 210. * rho8) * cos2theta; - break; - case 45: return rho2 * (15. - 140. * rho2 + 420. * rho4 - 504. * rho6 + 210. * rho8) * sin2theta; - break; - case 46: return rho *(-6. + 105 * rho2 - 560. * rho4 + 1260. * rho6 -1260. * rho8 +462. * rho10) * costheta; - break; - case 47: return rho *(-6. + 105 * rho2 - 560. * rho4 + 1260. * rho6 -1260. * rho8 +462. * rho10) * sintheta; - break; - case 48: return 1. - 42. * rho2 + 420. * rho4 - 1680. * rho6 + 3150. * rho8 -2772. * rho10 + 924. * pow(rho,12.); - break; - - } - return 0.; - - -} - - -double Zernike(int n, double X, double Y) +/*double Zernike(int n, double X, double Y) { - - static double X2 = 0., X3 = 0., X4 = 0.; static double Y2 = 0., Y3 = 0., Y4 = 0.; static double R2 = 0.; @@ -584,9 +420,7 @@ double Zernike(int n, double X, double Y) return(0.); } return d; - - -} +}*/ @@ -663,7 +497,6 @@ void zernikeProcess::unwrap_to_zernikes(wavefront &wf, int zterms){ } double delta = 1./(wf.m_outside.m_radius); - zernikePolar &zpolar = *zernikePolar::get_Instance(); int sampleCnt = 0; for(int y = 0; y < ny; y += step) //for each point on the surface { @@ -676,21 +509,21 @@ void zernikeProcess::unwrap_to_zernikes(wavefront &wf, int zterms){ if ( rho <= 1. && (wf.mask.at(y,x) != 0) && wf.data.at(y,x) != 0.0){ double theta = atan2(uy,ux); - zpolar.init(rho, theta); + zernikePolar zpolar(rho, theta, zterms); for ( int i = 0; i < zterms; ++i) { if (useSvd){ - double t = zpolar.zernike(i, rho, theta); + double t = zpolar.zernike(i); A(sampleCnt,i) = t; } else { - double t = zpolar.zernike(i, rho, theta); + double t = zpolar.zernike(i); for (int j = 0; j < zterms; ++j) { - //Am[ndx] = Am[ndx] + t * zpolar.zernike(j, rho, theta); - A(i,j) += t * zpolar.zernike(j, rho, theta); + //Am[ndx] = Am[ndx] + t * zpolar.zernike(j); + A(i,j) += t * zpolar.zernike(j); } // FN is the OPD at (Xn,Yn) //Bm[i] = Bm[i] + surface.at(y,x) * t; @@ -763,7 +596,6 @@ cv::Mat zernikeProcess::null_unwrapped(wavefront&wf, std::vector zerns, double sz,nz; double rho,theta; std::vector rows, cols; - zernikePolar &zpolar = *zernikePolar::get_Instance(); // make a list of points on the surface containing their rho and theta values as well as their // row column indexes in the matix that contanis the wave front. // annular wave fronts already have this made elsewhere. @@ -779,11 +611,14 @@ cv::Mat zernikeProcess::null_unwrapped(wavefront&wf, std::vector zerns, // this mask test may not be needed any longer but don't have time to check that. if (mask.at(y,x) != 0 && wf.data.at(y,x) != 0.0){ - rho = m_rhoTheta.row(0)(i); - theta = m_rhoTheta.row(1)(i); - if (!md->m_useAnnular){ - zpolar.init(rho,theta); - } + // Declare zernikePolar pointer, construct later if needed + zernikePolar* zpolar = nullptr; + + rho = m_rhoTheta.row(0)(i); + theta = m_rhoTheta.row(1)(i); + if (!md->m_useAnnular){ + zpolar = new zernikePolar(rho,theta, Z_TERMS); + } sz = unwrapped.at(y,x); nz = 0; @@ -792,7 +627,7 @@ cv::Mat zernikeProcess::null_unwrapped(wavefront&wf, std::vector zerns, { if (md->doNull && enables[8]){ if (!md->m_useAnnular) - nz -= scz8 * zpolar.zernike(8,rho, theta); + nz -= scz8 * zpolar->zernike(8); else { nz -= scz8 * m_zerns(i, 8); } @@ -802,21 +637,27 @@ cv::Mat zernikeProcess::null_unwrapped(wavefront&wf, std::vector zerns, for (int z = start_term; z < Z_TERMS; ++z) { if ((z == 3) && doDefocus){ - nz += defocus * zpolar.zernike(z,rho, theta); - nz -= zerns[z] * zpolar.zernike(z,rho,theta); - + if (!md->m_useAnnular) { + nz += defocus * zpolar->zernike(z); + nz -= zerns[z] * zpolar->zernike(z); + } + else { + nz += defocus * m_zerns(i,z); + nz -= zerns[z] * m_zerns(i,z); + } } - else if (!enables[z]){ - if (!md->m_useAnnular) - nz -= zerns[z] * zpolar.zernike(z,rho, theta); + if (!md->m_useAnnular) { + nz -= zerns[z] * zpolar->zernike(z); + } else { - - nz -= zerns[z] * m_zerns(i,z) ; - + nz -= zerns[z] * m_zerns(i,z); } } + } + if (!md->m_useAnnular){ + delete zpolar; } nulled.at(y,x) = sz +nz; @@ -833,8 +674,6 @@ void zernikeProcess::fillVoid(wavefront &wf){ mirrorDlg *md = mirrorDlg::get_Instance(); bool useannular = md->m_useAnnular; - zernikePolar &zpolar = *zernikePolar::get_Instance(); - if (wf.regions.size() > 0){ int x = wf.regions[0][0].x; int y = wf.mask.rows - wf.regions[0][0].y; @@ -885,11 +724,11 @@ void zernikeProcess::fillVoid(wavefront &wf){ theY.push_back(y); } else { - zpolar.init(rho,theta); + zernikePolar zpolar(rho, theta, wf.InputZerns.size()); double v = 0.; for (size_t z = 0; z < wf.InputZerns.size(); ++z){ - v += wf.InputZerns[z] * zpolar.zernike(z,rho, theta); + v += wf.InputZerns[z] * zpolar.zernike(z); } wf.data.at(y,x) = v; } @@ -969,11 +808,11 @@ void zernikeProcess::fillVoid(wavefront &wf){ theY.push_back(y); } else{ - zpolar.init(rho,theta); + zernikePolar zpolar(rho, theta, wf.InputZerns.size()); double v = 0.; for (size_t z = 0; z < wf.InputZerns.size(); ++z){ - v += wf.InputZerns[z] * zpolar.zernike(z,rho, theta); + v += wf.InputZerns[z] * zpolar.zernike(z); } wf.data.at(y,x) = v; } diff --git a/zernikeprocess.h b/zernikeprocess.h index 5b679965..75128cc3 100644 --- a/zernikeprocess.h +++ b/zernikeprocess.h @@ -25,29 +25,12 @@ #include "mainwindow.h" #include "armadillo" #include + +//TODO remove this if possible and use a class. It's undocumented, difficult to undesrstand and maintain extern std::vector zernEnables; -extern int Zw[]; -extern double BestSC; -double zernike(int n, double x, double y); -void gauss_jordan(int n, double* Am, double* Bm); -void ZernikeSmooth(cv::Mat wf, cv::Mat mask); - -typedef struct { - std::vector enables; - std::vector norms; - int maxOrder; - int noOfTerms; - std::vector rowIndex; - std::vector colIndex; - arma::mat zernsArSample; - arma::mat rhoTheta; - int widthOfMatrix; - double radius; - double centerX; - double centerY; - double insideRadius; - bool valid; -} zernikeData; +//double zernike(int n, double x, double y); +//void gauss_jordan(int n, double* Am, double* Bm); + class zernikeProcess : public QObject @@ -100,35 +83,5 @@ public slots: }; -class zernikePolar : public QObject -{ - Q_OBJECT -public: - explicit zernikePolar(){}; - static zernikePolar *get_Instance(); - void init(double rho, double theta); - double zernike(int z, double rho, double theta); -private: - static zernikePolar *m_instance; - double rho2; - double rho3; - double rho4; - double rho5; - double rho6; - double rho8; - double rho10; - double costheta; - double sintheta; - double cos2theta; - double sin2theta; - double cos3theta; - double sin3theta; - double cos4theta; - double sin4theta; - double cos5theta; - double sin5theta; -}; - -void debugZernRoutines(); #endif // ZERNIKEPROCESS_H diff --git a/zernikes.cpp b/zernikes.cpp index c0bbc141..fbe873ba 100644 --- a/zernikes.cpp +++ b/zernikes.cpp @@ -38,6 +38,15 @@ void dump_matrix (double *a, int nrows, int ncols,const char *desc) } } +static const int Zw[] = { /* n */ + 1, // 0 + 4,4,3, // 1 + 6,6,8,8,5, //8 // 2 + 8,8,10,10,12,12,7, // 15 3 + 10,10,12,12,14,14,16,16,9, //24 4 + 12,12,14,14,16,16,18,18,20,20,11, //35 5 + 14,14,16,16,18,18,20,20,22,22,24,24,13 // 48 +}; double computeRMS(int term, double val){ @@ -251,7 +260,7 @@ double *rand_zern ( zern_spec *s) void zern_generator::set_zcoefs(double *data) { - m_zcoef.assign(data,data + get_terms_cnt()); // temp. chanage to param after debug + m_zcoef.assign(data,data + get_terms_cnt()); // temp. chanage to param after debug }