Improve DFTFringe speed on zernike related computation

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

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 \

dftarea.cpp

@@ -140,7 +140,7 @@
     QSettings set;
     m_center_filter = set.value("DFT Center Filter", 10).toDouble();
     qDebug() << "init center" << m_center_filter;
     emit updateFilterSize(m_center_filter);
     installEventFilter(this);
 
     /*
@@ -1471,13 +1471,11 @@
     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" <<i <<  coords(0,i) << zpolar.zernike(i,r0, t0);
+//        qDebug() << "psi Zernike" <<i <<  coords(0,i) << zpolar.zernike(i);
 //    }
 
     //cv::Mat modul = modu.reshape(0,rows);

mainwindow.cpp

@@ -567,6 +567,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;
     }

surfacemanager.cpp

@@ -42,6 +42,7 @@
 #include <algorithm>
 #include <vector>
 #include "zernikeprocess.h"
+#include "zernikepolar.h"
 #include <QTimer>
 #include <qprinter.h>
 #include "rotationdlg.h"
@@ -434,7 +435,6 @@
 
     std::vector<bool> &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 @@
             {
                 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);
                         }
                     }
                 }
@@ -549,7 +549,7 @@
 
 SurfaceManager::~SurfaceManager(){
     spdlog::get("logger")->trace("SurfaceManager::~SurfaceManager");
     for(wavefront* wf : m_wavefronts){
         delete wf;
     }
     if(m_standAstigWizard != nullptr){
@@ -744,7 +744,6 @@
     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 @@
             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<double>(j,i) = S1;
             }
@@ -1549,7 +1548,7 @@
     QString maxkey;
     qDebug() << "sizes" << sizes;
 
     foreach(QString v, sizes.keys()){
         int a = sizes[v].length();
         if (a > max) {
             max = a;

zernikeeditdlg.cpp

@@ -3,6 +3,7 @@
 #include <QFileDialog>
 #include <fstream>
 #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<double>(y,x) = s1;

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 <http://www.gnu.org/licenses/>.
+
+****************************************************************************/
+
+#include "zernikepolar.h"
+#include <cmath>
+#include <QDebug>
+
+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.;
+    }
+}

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 <http://www.gnu.org/licenses/>.
+
+****************************************************************************/
+
+#ifndef ZERNIKEPOLAR_H
+#define ZERNIKEPOLAR_H
+
+#include <cstddef>
+
+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