From 74ad0344ce64f650d9d710ca35e250349f4f141c Mon Sep 17 00:00:00 2001 From: Julien Staub Date: Sat, 19 Jul 2025 16:17:22 +0200 Subject: [PATCH 1/8] don't use pow which is slow --- zernikeprocess.cpp | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/zernikeprocess.cpp b/zernikeprocess.cpp index 1ad0bfff..f43401a3 100644 --- a/zernikeprocess.cpp +++ b/zernikeprocess.cpp @@ -275,12 +275,12 @@ zernikePolar *zernikePolar::get_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.); + rho3 = rho2 * rho; + rho4 = rho3 * rho; + rho5 = rho4 * rho; + rho6 = rho5 * rho; + rho8 = rho6 * rho; + rho10 = rho8 * rho2; costheta = cos(theta); sintheta = sin(theta); cos2theta = cos(2. * theta); From e463a8ff6fbeb8011349b3c56bf55b3291a27332 Mon Sep 17 00:00:00 2001 From: Julien Staub Date: Sat, 19 Jul 2025 17:11:46 +0200 Subject: [PATCH 2/8] optimize zernike computation --- dftarea.cpp | 4 +- surfacemanager.cpp | 12 +-- zernikeeditdlg.cpp | 4 +- zernikeprocess.cpp | 239 +++++++++++++++++++-------------------------- zernikeprocess.h | 11 ++- 5 files changed, 118 insertions(+), 152 deletions(-) diff --git a/dftarea.cpp b/dftarea.cpp index ad151f09..49e427d7 100644 --- a/dftarea.cpp +++ b/dftarea.cpp @@ -1475,11 +1475,11 @@ qDebug() << "rec" << left << top << width << height; zernikePolar &zpolar = *zernikePolar::get_Instance(); - zpolar.init(r0,t0); + zpolar.init(r0,t0,100); // debug print the difference // for (int i = 0; i < 100; ++i){ -// qDebug() << "psi Zernike" <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); } } } @@ -753,11 +753,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); + zpolar.init(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..233e52c0 100644 --- a/zernikeeditdlg.cpp +++ b/zernikeeditdlg.cpp @@ -65,12 +65,12 @@ void zernikeEditDlg::on_createSurface_clicked() if (rho <= 1.) { double theta = atan2(uy,ux); - zpolar.init(rho,theta); + zpolar.init(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/zernikeprocess.cpp b/zernikeprocess.cpp index f43401a3..9f704f03 100644 --- a/zernikeprocess.cpp +++ b/zernikeprocess.cpp @@ -273,153 +273,114 @@ zernikePolar *zernikePolar::get_Instance(){ return m_instance; } -void zernikePolar::init(double rho, double theta){ - rho2 = rho * rho; +// Having all terms computed at once here let's compiler optimize the code better +void zernikePolar::init(double rho, double theta, size_t nbTerms){ + if(nbTerms > 49){ + qWarning() << "zernikePolar::init: maxorder is limited to 49, setting to 49"; + nbTerms = 49; + } + + m_nbTermsComputed = nbTerms; + + zernTerms[0] = 1.; + + rho2 = rho * rho; + costheta = cos(theta); + sintheta = sin(theta); + cos2theta = cos(2. * theta); + sin2theta = 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); + + // 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 * rho; - rho10 = rho8 * rho2; - 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); + + 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 = cos(5. * theta); sin5theta = sin(5. * theta); - } -double zernikePolar::zernike(int n, double rho, double 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; + } - 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: + if(nbTerms > 35) { - return 1. + rho2 * (-6 + 6. * rho2); - break; + zernTerms[36] = rho6 * cos(6. * theta); + zernTerms[37] = rho6 * 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. * rho8 * rho2; } - 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; +} +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.; } - return 0.; - - } 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.; @@ -676,21 +637,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); + zpolar.init(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; @@ -782,7 +743,7 @@ cv::Mat zernikeProcess::null_unwrapped(wavefront&wf, std::vector zerns, rho = m_rhoTheta.row(0)(i); theta = m_rhoTheta.row(1)(i); if (!md->m_useAnnular){ - zpolar.init(rho,theta); + zpolar.init(rho,theta, Z_TERMS); } sz = unwrapped.at(y,x); @@ -792,7 +753,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,14 +763,14 @@ 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); + nz += defocus * zpolar.zernike(z); + nz -= zerns[z] * zpolar.zernike(z); } else if (!enables[z]){ if (!md->m_useAnnular) - nz -= zerns[z] * zpolar.zernike(z,rho, theta); + nz -= zerns[z] * zpolar.zernike(z); else { nz -= zerns[z] * m_zerns(i,z) ; @@ -885,11 +846,11 @@ void zernikeProcess::fillVoid(wavefront &wf){ theY.push_back(y); } else { - zpolar.init(rho,theta); + zpolar.init(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 +930,11 @@ void zernikeProcess::fillVoid(wavefront &wf){ theY.push_back(y); } else{ - zpolar.init(rho,theta); + zpolar.init(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 5dadb91e..b57e1302 100644 --- a/zernikeprocess.h +++ b/zernikeprocess.h @@ -25,6 +25,8 @@ #include "mainwindow.h" #include "armadillo" #include + +//TODO remove this if possible and use a class extern std::vector zernEnables; extern int Zw[]; extern double BestSC; @@ -137,14 +139,15 @@ public slots: }; +//TODO separate class into different files 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); + static zernikePolar *get_Instance(); //TODO check if we really need a singleton here + void init(double rho, double theta, size_t nbTerms = 48); + double zernike(size_t z); private: static zernikePolar *m_instance; double rho2; @@ -164,6 +167,8 @@ class zernikePolar : public QObject double sin4theta; double cos5theta; double sin5theta; + size_t m_nbTermsComputed; + double zernTerms[49]; // TODO if no singleton, use correct length }; void debugZernRoutines(); From 24b2c78e17ca680c95c19820f427f52e763d18cc Mon Sep 17 00:00:00 2001 From: Julien Staub Date: Sat, 19 Jul 2025 17:57:21 +0200 Subject: [PATCH 3/8] separate differnt class in different files --- DFTFringe.pro | 2 + DFTFringe_Dale.pro | 2 + surfacemanager.cpp | 1 + zernikeeditdlg.cpp | 1 + zernikepolar.cpp | 137 +++++++++++++++++++++++++++++++++++++++++++++ zernikepolar.h | 56 ++++++++++++++++++ zernikeprocess.cpp | 116 +------------------------------------- zernikeprocess.h | 31 +--------- 8 files changed, 201 insertions(+), 145 deletions(-) create mode 100644 zernikepolar.cpp create mode 100644 zernikepolar.h diff --git a/DFTFringe.pro b/DFTFringe.pro index cf7ec6e0..f85f7a71 100644 --- a/DFTFringe.pro +++ b/DFTFringe.pro @@ -262,6 +262,7 @@ SOURCES += SingleApplication/singleapplication.cpp \ zapm.cpp \ zernikedlg.cpp \ zernikeeditdlg.cpp \ + zernikepolar.cpp \ zernikeprocess.cpp \ zernikes.cpp \ zernikesmoothingdlg.cpp @@ -389,6 +390,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 c4307ac2..25c9280f 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/surfacemanager.cpp b/surfacemanager.cpp index ccb4cd65..3b94899f 100644 --- a/surfacemanager.cpp +++ b/surfacemanager.cpp @@ -41,6 +41,7 @@ #include #include #include "zernikeprocess.h" +#include "zernikepolar.h" #include #include #include "rotationdlg.h" diff --git a/zernikeeditdlg.cpp b/zernikeeditdlg.cpp index 233e52c0..5b9c359c 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) : diff --git a/zernikepolar.cpp b/zernikepolar.cpp new file mode 100644 index 00000000..f4ba4d9d --- /dev/null +++ b/zernikepolar.cpp @@ -0,0 +1,137 @@ +/****************************************************************************** +** +** 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 + +/* +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; +} + +// Having all terms computed at once here let's compiler optimize the code better +void zernikePolar::init(double rho, double theta, size_t nbTerms){ + if(nbTerms > 49){ + qWarning() << "zernikePolar::init: maxorder is limited to 49, setting to 49"; + nbTerms = 49; + } + + m_nbTermsComputed = nbTerms; + + zernTerms[0] = 1.; + + rho2 = rho * rho; + costheta = cos(theta); + sintheta = sin(theta); + cos2theta = cos(2. * theta); + sin2theta = 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); + + // 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 * rho; + cos3theta = cos(3. * theta); + sin3theta = sin(3. * theta); + cos4theta = cos(4. * theta); + sin4theta = 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 = cos(5. * theta); + sin5theta = 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 * cos(6. * theta); + zernTerms[37] = rho6 * 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. * rho8 * 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..b859ef36 --- /dev/null +++ b/zernikepolar.h @@ -0,0 +1,56 @@ +/****************************************************************************** +** +** 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 QObject +{ + Q_OBJECT +public: + explicit zernikePolar(){}; + static zernikePolar *get_Instance(); //TODO check if we really need a singleton here + void init(double rho, double theta, size_t nbTerms = 48); + double zernike(size_t z); +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; + size_t m_nbTermsComputed; + double zernTerms[49]; // TODO if no singleton, use correct length +}; + + +#endif \ No newline at end of file diff --git a/zernikeprocess.cpp b/zernikeprocess.cpp index 9f704f03..85811a39 100644 --- a/zernikeprocess.cpp +++ b/zernikeprocess.cpp @@ -16,6 +16,7 @@ ****************************************************************************/ #include "zernikeprocess.h" +#include "zernikepolar.h" #include #include #include "mainwindow.h" @@ -261,122 +262,7 @@ Dim big As Double, dum As Double, pivinv As Double, temp As Double 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; -} -// Having all terms computed at once here let's compiler optimize the code better -void zernikePolar::init(double rho, double theta, size_t nbTerms){ - if(nbTerms > 49){ - qWarning() << "zernikePolar::init: maxorder is limited to 49, setting to 49"; - nbTerms = 49; - } - - m_nbTermsComputed = nbTerms; - - zernTerms[0] = 1.; - - rho2 = rho * rho; - costheta = cos(theta); - sintheta = sin(theta); - cos2theta = cos(2. * theta); - sin2theta = 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); - - // 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 * rho; - cos3theta = cos(3. * theta); - sin3theta = sin(3. * theta); - cos4theta = cos(4. * theta); - sin4theta = 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 = cos(5. * theta); - sin5theta = 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 * cos(6. * theta); - zernTerms[37] = rho6 * 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. * rho8 * 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.; - } -} double Zernike(int n, double X, double Y) diff --git a/zernikeprocess.h b/zernikeprocess.h index b57e1302..625de234 100644 --- a/zernikeprocess.h +++ b/zernikeprocess.h @@ -140,36 +140,7 @@ public slots: }; //TODO separate class into different files -class zernikePolar : public QObject -{ - Q_OBJECT -public: - explicit zernikePolar(){}; - static zernikePolar *get_Instance(); //TODO check if we really need a singleton here - void init(double rho, double theta, size_t nbTerms = 48); - double zernike(size_t z); -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; - size_t m_nbTermsComputed; - double zernTerms[49]; // TODO if no singleton, use correct length -}; + void debugZernRoutines(); From 4e6be683022168e28578344f293d69f22156271e Mon Sep 17 00:00:00 2001 From: Julien Staub Date: Sat, 19 Jul 2025 18:11:46 +0200 Subject: [PATCH 4/8] some cleanup --- mainwindow.cpp | 1 + zernikepolar.cpp | 5 +--- zernikeprocess.cpp | 45 +++++++++++-------------------- zernikeprocess.h | 66 +++------------------------------------------- zernikes.cpp | 11 +++++++- 5 files changed, 31 insertions(+), 97 deletions(-) diff --git a/mainwindow.cpp b/mainwindow.cpp index 4794b430..4acdd422 100644 --- a/mainwindow.cpp +++ b/mainwindow.cpp @@ -568,6 +568,7 @@ void MainWindow::updateMetrics(wavefront& wf){ } metrics->setZernTitle(ztitle); double z8 = zernTablemodel->values[8]; + double BestSC; if (m_mirrorDlg->doNull && wf.useSANull){ BestSC = z8/m_mirrorDlg->z8; } diff --git a/zernikepolar.cpp b/zernikepolar.cpp index f4ba4d9d..2189ae40 100644 --- a/zernikepolar.cpp +++ b/zernikepolar.cpp @@ -19,10 +19,7 @@ #include "zernikepolar.h" #include -/* -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){ diff --git a/zernikeprocess.cpp b/zernikeprocess.cpp index 85811a39..5b26eba8 100644 --- a/zernikeprocess.cpp +++ b/zernikeprocess.cpp @@ -27,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 -}; - // @@ -132,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) { @@ -152,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 @@ -168,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]; @@ -261,11 +247,14 @@ Dim big As Double, dum As Double, pivinv As Double, temp As Double delete[] ipiv; delete[] indxr; delete[] indxc; -} - +}*/ -double Zernike(int n, double X, double Y) +/* +Public Function Zernike(n As Integer, X As Double, Y As Double) As Double +' N is Zernike number as given by James Wyant +*/ +/*double Zernike(int n, double X, double Y) { static double X2 = 0., X3 = 0., X4 = 0.; static double Y2 = 0., Y3 = 0., Y4 = 0.; @@ -431,9 +420,7 @@ double Zernike(int n, double X, double Y) return(0.); } return d; - - -} +}*/ diff --git a/zernikeprocess.h b/zernikeprocess.h index 625de234..75128cc3 100644 --- a/zernikeprocess.h +++ b/zernikeprocess.h @@ -26,67 +26,11 @@ #include "armadillo" #include -//TODO remove this if possible and use a class +//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(Mat wf, 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; - -class zernProcess : public QObject -{ - Q_OBJECT -private: - static zernProcess *m_Instance; - bool m_needsInit; +//double zernike(int n, double x, double y); +//void gauss_jordan(int n, double* Am, double* Bm); -public: - zernikeData m_zData; - - explicit zernProcess(QObject *parent = 0); - static zernProcess *get_Instance(); - void unwrap_to_zernikes(wavefront &wf, int zterms = Z_TERMS); - cv::Mat null_unwrapped(wavefront&wf, std::vector zerns, std::vector enables,int start_term =0, int last_term = Z_TERMS); - void ZernFitWavefront( wavefront &wf); - void initGrid(wavefront &wf, int maxOrder); - void initGrid(int width, double radius, double cx, double cy, int maxOrder, double inside = 0); - void unwrap_to_zernikes(zern_generator *zg, cv::Mat wf, cv::Mat mask); - cv::Mat makeSurfaceFromZerns(int border, bool doColor); - //arma::mat rhotheta( int width, double radius, double cx, double cy, - //double insideRad, const wavefront *wf = 0); - //arma::mat zpmC(arma::rowvec rho, arma::rowvec theta, int maxorder); - //arma::mat zapmC(const arma::rowvec& rho, const arma::rowvec& theta, const int& maxorder=12); - void fillVoid(wavefront &wf); - void setMaxOrder(int n); - int getNumberOfTerms(); - int m_abc; - mirrorDlg *md; - MainWindow *mw; - //arma::mat m_zerns; - QVector m_norms; -signals: -void statusBarUpdate(QString, int); -public slots: - -}; class zernikeProcess : public QObject @@ -139,9 +83,5 @@ public slots: }; -//TODO separate class into different files - - -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 } From 7e11ad03d71b33bebdaede33fbfd32c824521dd8 Mon Sep 17 00:00:00 2001 From: Julien Staub Date: Sat, 19 Jul 2025 18:35:36 +0200 Subject: [PATCH 5/8] fix build --- zernikepolar.cpp | 25 +++++++++++++------------ 1 file changed, 13 insertions(+), 12 deletions(-) diff --git a/zernikepolar.cpp b/zernikepolar.cpp index 2189ae40..28ce14ce 100644 --- a/zernikepolar.cpp +++ b/zernikepolar.cpp @@ -17,6 +17,7 @@ ****************************************************************************/ #include "zernikepolar.h" +#include #include @@ -40,10 +41,10 @@ void zernikePolar::init(double rho, double theta, size_t nbTerms){ zernTerms[0] = 1.; rho2 = rho * rho; - costheta = cos(theta); - sintheta = sin(theta); - cos2theta = cos(2. * theta); - sin2theta = sin(2. * theta); + costheta = std::cos(theta); + sintheta = std::sin(theta); + cos2theta = std::cos(2. * theta); + sin2theta = std::sin(2. * theta); zernTerms[1] = rho * costheta; zernTerms[2] = rho * sintheta; zernTerms[3] = -1. + 2. * rho2; @@ -63,10 +64,10 @@ void zernikePolar::init(double rho, double theta, size_t nbTerms){ rho5 = rho4 * rho; rho6 = rho5 * rho; rho8 = rho6 * rho; - cos3theta = cos(3. * theta); - sin3theta = sin(3. * theta); - cos4theta = cos(4. * theta); - sin4theta = sin(4. * theta); + 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; @@ -88,8 +89,8 @@ void zernikePolar::init(double rho, double theta, size_t nbTerms){ if(nbTerms > 24) { rho10 = rho8 * rho2; - cos5theta = cos(5. * theta); - sin5theta = sin(5. * theta); + cos5theta = std::cos(5. * theta); + sin5theta = std::sin(5. * theta); zernTerms[25] = rho5 * cos5theta; zernTerms[26] = rho5 * sin5theta; @@ -106,8 +107,8 @@ void zernikePolar::init(double rho, double theta, size_t nbTerms){ if(nbTerms > 35) { - zernTerms[36] = rho6 * cos(6. * theta); - zernTerms[37] = rho6 * sin(6. * theta); + 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; From 84bee4641e635d2f653e51be8985b569e131c4a4 Mon Sep 17 00:00:00 2001 From: Julien Staub Date: Sat, 19 Jul 2025 19:35:41 +0200 Subject: [PATCH 6/8] fix typo in zerns --- zernikepolar.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/zernikepolar.cpp b/zernikepolar.cpp index 28ce14ce..9f43ab3a 100644 --- a/zernikepolar.cpp +++ b/zernikepolar.cpp @@ -63,7 +63,7 @@ void zernikePolar::init(double rho, double theta, size_t nbTerms){ rho4 = rho3 * rho; rho5 = rho4 * rho; rho6 = rho5 * rho; - rho8 = rho6 * rho; + rho8 = rho6 * rho2; cos3theta = std::cos(3. * theta); sin3theta = std::sin(3. * theta); cos4theta = std::cos(4. * theta); From 385380a86741a6a16e8ed9a141e344579ca2cdaa Mon Sep 17 00:00:00 2001 From: Julien Staub Date: Sat, 19 Jul 2025 20:19:24 +0200 Subject: [PATCH 7/8] typo --- zernikepolar.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/zernikepolar.cpp b/zernikepolar.cpp index 9f43ab3a..a6e67945 100644 --- a/zernikepolar.cpp +++ b/zernikepolar.cpp @@ -119,7 +119,7 @@ void zernikePolar::init(double rho, double theta, size_t nbTerms){ 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. * rho8 * rho2; + zernTerms[48] = 1. - 42. * rho2 + 420. * rho4 - 1680. * rho6 + 3150. * rho8 -2772. * rho10 + 924. * rho10 * rho2; } } From 155dfa99afa687897f8bd1931b66d8199a75239b Mon Sep 17 00:00:00 2001 From: Julien Staub Date: Sun, 20 Jul 2025 09:06:32 +0200 Subject: [PATCH 8/8] made a normal C++ class. no singleton, no QT --- dftarea.cpp | 4 +--- surfacemanager.cpp | 6 ++---- zernikeeditdlg.cpp | 3 +-- zernikepolar.cpp | 36 +++++++++++++++++++--------------- zernikepolar.h | 31 +++++------------------------ zernikeprocess.cpp | 49 +++++++++++++++++++++++++--------------------- 6 files changed, 56 insertions(+), 73 deletions(-) diff --git a/dftarea.cpp b/dftarea.cpp index 49e427d7..71521a7a 100644 --- a/dftarea.cpp +++ b/dftarea.cpp @@ -1473,9 +1473,7 @@ qDebug() << "rec" << left << top << width << height; qDebug() << "Rho " << r0 << "Theta" << t0; - - zernikePolar &zpolar = *zernikePolar::get_Instance(); - zpolar.init(r0,t0,100); + zernikePolar zpolar(r0,t0,100); // debug print the difference // for (int i = 0; i < 100; ++i){ diff --git a/surfacemanager.cpp b/surfacemanager.cpp index 3b94899f..fc16d20a 100644 --- a/surfacemanager.cpp +++ b/surfacemanager.cpp @@ -435,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; @@ -448,7 +447,7 @@ cv::Mat SurfaceManager::computeWaveFrontFromZernikes(int wx, int wy, std::vector { double S1 = 0; double theta = atan2(y1,x1); - zpolar.init(rho,theta,zernsToUse.size()); + zernikePolar zpolar(rho, theta, zernsToUse.size()); for (int ii = 0; ii < zernsToUse.size(); ++ii) { int z = zernsToUse[ii]; @@ -745,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; @@ -754,7 +752,7 @@ void SurfaceManager::useDemoWaveFront(){ double y1 = (double)(j - (ycen )) /rad; rho = sqrt(x1 * x1 + y1 * y1); double theta = atan2(y1,x1); - zpolar.init(rho,theta,10); + zernikePolar zpolar(rho, theta, 10); if (rho <= 1.) { diff --git a/zernikeeditdlg.cpp b/zernikeeditdlg.cpp index 5b9c359c..21e4c3c8 100644 --- a/zernikeeditdlg.cpp +++ b/zernikeeditdlg.cpp @@ -54,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; @@ -66,7 +65,7 @@ void zernikeEditDlg::on_createSurface_clicked() if (rho <= 1.) { double theta = atan2(uy,ux); - zpolar.init(rho, theta, tableModel->rowCount()); + zernikePolar zpolar(rho, theta, tableModel->rowCount()); double s1 = 0; for (int z = 0; z < tableModel->rowCount(); ++z){ double v = tableModel->values[z]; diff --git a/zernikepolar.cpp b/zernikepolar.cpp index a6e67945..633484e0 100644 --- a/zernikepolar.cpp +++ b/zernikepolar.cpp @@ -20,17 +20,9 @@ #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 -zernikePolar *zernikePolar::m_instance = 0; -zernikePolar *zernikePolar::get_Instance(){ - if (m_instance == 0){ - m_instance = new zernikePolar; - } - return m_instance; -} - -// Having all terms computed at once here let's compiler optimize the code better -void zernikePolar::init(double rho, double theta, size_t nbTerms){ if(nbTerms > 49){ qWarning() << "zernikePolar::init: maxorder is limited to 49, setting to 49"; nbTerms = 49; @@ -40,11 +32,11 @@ void zernikePolar::init(double rho, double theta, size_t nbTerms){ zernTerms[0] = 1.; - rho2 = rho * rho; - costheta = std::cos(theta); - sintheta = std::sin(theta); - cos2theta = std::cos(2. * theta); - sin2theta = std::sin(2. * theta); + 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; @@ -54,11 +46,23 @@ void zernikePolar::init(double rho, double theta, size_t nbTerms){ 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; diff --git a/zernikepolar.h b/zernikepolar.h index b859ef36..c832c03e 100644 --- a/zernikepolar.h +++ b/zernikepolar.h @@ -19,37 +19,16 @@ #ifndef ZERNIKEPOLAR_H #define ZERNIKEPOLAR_H -#include +#include -class zernikePolar : public QObject +class zernikePolar { - Q_OBJECT public: - explicit zernikePolar(){}; - static zernikePolar *get_Instance(); //TODO check if we really need a singleton here - void init(double rho, double theta, size_t nbTerms = 48); + zernikePolar(double rho, double theta, size_t nbTerms = 48); double zernike(size_t z); 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; - size_t m_nbTermsComputed; - double zernTerms[49]; // TODO if no singleton, use correct length + size_t m_nbTermsComputed = 0; + double zernTerms[49]; }; diff --git a/zernikeprocess.cpp b/zernikeprocess.cpp index 5b26eba8..427fd19b 100644 --- a/zernikeprocess.cpp +++ b/zernikeprocess.cpp @@ -497,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 { @@ -510,7 +509,7 @@ 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, zterms); + zernikePolar zpolar(rho, theta, zterms); for ( int i = 0; i < zterms; ++i) { @@ -597,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. @@ -613,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, Z_TERMS); - } + // 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; @@ -626,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); + nz -= scz8 * zpolar->zernike(8); else { nz -= scz8 * m_zerns(i, 8); } @@ -636,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); - nz -= zerns[z] * zpolar.zernike(z); - + 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); + 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; @@ -667,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; @@ -719,7 +724,7 @@ void zernikeProcess::fillVoid(wavefront &wf){ theY.push_back(y); } else { - zpolar.init(rho, theta, wf.InputZerns.size()); + zernikePolar zpolar(rho, theta, wf.InputZerns.size()); double v = 0.; for (size_t z = 0; z < wf.InputZerns.size(); ++z){ @@ -803,7 +808,7 @@ void zernikeProcess::fillVoid(wavefront &wf){ theY.push_back(y); } else{ - zpolar.init(rho, theta, wf.InputZerns.size()); + zernikePolar zpolar(rho, theta, wf.InputZerns.size()); double v = 0.; for (size_t z = 0; z < wf.InputZerns.size(); ++z){