diff --git a/DFTFringe.pro b/DFTFringe.pro index d3b3dc7c..e4af46a4 100644 --- a/DFTFringe.pro +++ b/DFTFringe.pro @@ -143,10 +143,11 @@ RESOURCES += DFTResources.qrc TRANSLATIONS += dftfringe_fr.ts -INCLUDEPATH += ./bezier ./SingleApplication +INCLUDEPATH += ./bezier ./SingleApplication ./zernike SOURCES += SingleApplication/singleapplication.cpp \ SingleApplication/singleapplication_p.cpp \ + zernike/zapm.cpp \ annulushelpdlg.cpp \ arbitrarywavefronthelp.cpp \ arbitrarywavwidget.cpp \ @@ -261,7 +262,6 @@ SOURCES += SingleApplication/singleapplication.cpp \ wavefrontloader.cpp \ wftexaminer.cpp \ wftstats.cpp \ - zapm.cpp \ zernikedlg.cpp \ zernikeeditdlg.cpp \ zernikepolar.cpp \ @@ -272,6 +272,7 @@ SOURCES += SingleApplication/singleapplication.cpp \ HEADERS += bezier/bezier.h \ SingleApplication/singleapplication_p.h \ SingleApplication/singleapplication.h \ + zernike/zapm_interface.h \ annulushelpdlg.h \ arbitrarywavefronthelp.h \ astigpolargraph.h \ diff --git a/DFTFringe_Dale.pro b/DFTFringe_Dale.pro index de04bf23..a88e4fe1 100644 --- a/DFTFringe_Dale.pro +++ b/DFTFringe_Dale.pro @@ -59,7 +59,6 @@ SOURCES += main.cpp \ dftcolormap.cpp \ surfaceanalysistools.cpp \ surfacemanager.cpp \ - zapm.cpp \ zernikedlg.cpp \ zernikepolar.cpp \ zernikeprocess.cpp \ @@ -147,6 +146,7 @@ SOURCES += main.cpp \ psitiltoptions.cpp \ contourrulerparams.cpp \ zernikesmoothingdlg.cpp \ + zernike/zapm.cpp \ SingleApplication/singleapplication.cpp \ SingleApplication/singleapplication_p.cpp @@ -275,10 +275,11 @@ HEADERS += mainwindow.h \ contourrulerparams.h \ zernikesmoothingdlg.h \ bezier/bezier.h \ + zernike/zapm_interface.h \ SingleApplication/singleapplication.h \ SingleApplication/singleapplication_p.h -INCLUDEPATH += ./bezier ./SingleApplication +INCLUDEPATH += ./bezier ./SingleApplication ./zernike FORMS += mainwindow.ui \ annulushelpdlg.ui \ diff --git a/DFTFringe_QT5.pro b/DFTFringe_QT5.pro index 0271d3fc..72741196 100644 --- a/DFTFringe_QT5.pro +++ b/DFTFringe_QT5.pro @@ -142,10 +142,11 @@ RESOURCES += DFTResources.qrc TRANSLATIONS += dftfringe_fr.ts -INCLUDEPATH += ./bezier ./SingleApplication +INCLUDEPATH += ./bezier ./SingleApplication ./zernike SOURCES += SingleApplication/singleapplication.cpp \ SingleApplication/singleapplication_p.cpp \ + zernike/zapm.cpp \ annulushelpdlg.cpp \ arbitrarywavefronthelp.cpp \ arbitrarywavwidget.cpp \ @@ -260,7 +261,6 @@ SOURCES += SingleApplication/singleapplication.cpp \ wavefrontloader.cpp \ wftexaminer.cpp \ wftstats.cpp \ - zapm.cpp \ zernikedlg.cpp \ zernikeeditdlg.cpp \ zernikepolar.cpp \ @@ -271,6 +271,7 @@ SOURCES += SingleApplication/singleapplication.cpp \ HEADERS += bezier/bezier.h \ SingleApplication/singleapplication_p.h \ SingleApplication/singleapplication.h \ + zernike/zapm_interface.h \ annulushelpdlg.h \ arbitrarywavefronthelp.h \ astigpolargraph.h \ diff --git a/dftarea.cpp b/dftarea.cpp index 757ffe1d..bd024da9 100644 --- a/dftarea.cpp +++ b/dftarea.cpp @@ -1271,7 +1271,9 @@ cv::Mat atan2Mat(cv::Mat y, cv::Mat x){ } #include - +// zpmCxx commented out, not used, might be needed in doPSIstep4 +// I (atsju) recommend deleting this part and add original or modified file from Mike when/if needed +/* arma::mat zpmCxx(double rho, double theta, int maxorder) { int m, n, n0, mmax = maxorder/2; @@ -1345,6 +1347,7 @@ arma::mat zpmCxx(double rho, double theta, int maxorder) { } return zm; } +*/ #include "outlinedialog.h" void DFTArea::doPSIstep4(cv::Mat images, QVector phases){ diff --git a/mainwindow.cpp b/mainwindow.cpp index 8f267eda..d80c1c55 100644 --- a/mainwindow.cpp +++ b/mainwindow.cpp @@ -49,6 +49,8 @@ #include "colorchannel.h" #include "opencv2/opencv.hpp" #include +#include "zapm_interface.h" +#include "zernikeprocess.h" using namespace QtConcurrent; std::vector g_wavefronts; @@ -1845,11 +1847,7 @@ void MainWindow::on_actionCreate_Movie_of_wavefronts_triggered() this->setCursor(Qt::ArrowCursor); QApplication::restoreOverrideCursor(); } -arma::mat zapm(const arma::vec& rho, const arma::vec& theta, - const double& eps, const int& maxorder=12) ; -#include "armadillo" -void dumpArma(arma::mat mm, QString title = "", QVector colHeading = QVector(0), - QVector RowLable = QVector(0)); + void MainWindow::on_actionDebugStuff_triggered() { zernikeProcess *zp = zernikeProcess::get_Instance(); diff --git a/mikespsiinterface.h_ b/mikespsiinterface.h_ deleted file mode 100644 index fcc07e46..00000000 --- a/mikespsiinterface.h_ +++ /dev/null @@ -1,10 +0,0 @@ -#ifndef MIKESPSIINTERFACE_H -#define MIKESPSIINTERFACE_H -#include "armadillo" -arma::mat zpmC(arma::vec rho, arma::vec theta, int maxorder); - -typedef struct { arma::mat phi; arma::mat mod; arma::mat phase; arma::mat zcx; int iter; arma::vec sse;} psitiltReturn; -arma::mat pxls(const arma::mat& im, const arma::rowvec& phases, const arma::mat& zcs, const arma::mat& coords); -arma::mat pwrap(const arma::mat& phase); -arma::mat rhotheta(const double rows, const double cols, double obspercent=0); -#endif // MIKESPSIINTERFACE_H diff --git a/mikespsirinterface.cpp_ b/mikespsirinterface.cpp_ deleted file mode 100644 index 83ede20d..00000000 --- a/mikespsirinterface.cpp_ +++ /dev/null @@ -1,260 +0,0 @@ -// portions of this code are -// Copyright (C) 2016 Michael Peck -// -// License: MIT -// -#include -#include -#include "mikespsiinterface.h" -#include -//#include -#include -#include "qmath.h" -#include -#include -// create two row vectors containing rho and theta values of the inclosed circular pupil -// return them as two columns of a mat. - - - - -// [[Rcpp::export]] - - -typedef struct {arma::vec x; bool convergence; int iter;} lmreturn; - -/************** - * - * A stupidly simple implementation of a Levenberg-Marquardt - * algorithm for nonlinear least squares. - * - * Adapted directly from the algorithm description on p. 27 - * of the booklet "Methods for Non-Linear Least Squares - * Problems" by Madsen, Nielsen & Tingleff (2004) - * download from http://www2.imm.dtu.dk/pubdb/views/edoc_download.php/3215/pdf/imm3215.pdf - * -*/ - -lmreturn mylevmar(arma::vec (*res_fun)(const arma::vec&, const QMap&, double abar, double bbar), - arma::mat (*jac_fun)(const arma::vec&, const QMap&, double bbar), - unsigned int m, int n, arma::vec x0, QMap adata, double abar, double bbar, - int itmax = 100, arma::vec control= {1.e-3, 2., sqrt(arma::datum::eps), sqrt(arma::datum::eps)}) { - - int k=0; - double tau = control(0); - double nu = control(1); - double eps1 = control(2); - double eps2 = control(3); - - arma::mat A(n, n); - arma::vec g(n); - arma::vec feval(m); - arma::vec fnew(m); - arma::mat Jac(m, n); - arma::mat In = arma::eye(n, n); - arma::vec x = x0; - arma::vec hx(n), xnew(n); - bool found; - double mu; - double F, Fnew, dL, rho; - lmreturn retvals; - - Jac = jac_fun(x, adata,bbar); - feval = res_fun(x, adata, abar, bbar); - A = Jac.t() * Jac; - g = Jac.t() * feval; - mu = tau * max(A.diag()); - found = (norm(g, "Inf") <= eps1); - - while (!found && k 0.) { - x = xnew; - Jac = jac_fun(x, adata,bbar); - A = Jac.t() * Jac; - g = Jac.t() * fnew; - feval = fnew; - mu *= fmax(1./3., 1 - pow(2.*rho -1., 3)); - nu = control(1); - found = (arma::norm(g, "Inf") <= eps1); - } else { - mu *= nu; - nu *= 2.; - } - } - - retvals.x = x; - retvals.convergence = found; - retvals.iter = k; - - return retvals; - -} - -// pixel-wise least squares with known phases, tilts, etc. - -arma::mat pxls(const arma::mat& im, const arma::rowvec& phases, const arma::mat& zcs, const arma::mat& coords) { - unsigned int nr = im.n_rows; - unsigned int nf = im.n_cols; - - arma::rowvec ph(nf); - arma::mat A(3, nf); - arma::rowvec b(3); - arma::mat B(nr, 3); - - for (unsigned int n=0; n(nf), arma::cos(ph)), arma::sin(ph)); - b = im.row(n) * pinv(A); - B.row(n) = b; - } - return B; -} - -// residuals and analytic Jacobian for mylevmar. Exported just in case I want to use minpack.lm - -// [[Rcpp::export]] - -arma::vec res_frame(const arma::vec& pars, const QMap & adata, double abar, double bbar) { - arma::vec img = adata["img"]; - arma::mat coords = adata["coords"]; - - arma::vec phi = adata["phi"]; - int np = pars.n_elem; - arma::vec ph(img.n_elem); - ph = pars[0] + coords * pars.tail(np-1); - return img - abar - bbar * cos(phi + ph); -} - -// [[Rcpp::export]] - -arma::mat jac_frame(const arma::vec& pars, const QMap & adata, double bbar) { - arma::vec img = adata["img"]; - arma::mat coords = adata["coords"]; - - arma::vec phi = adata["phi"]; - int np = pars.n_elem; - unsigned int m = img.n_elem; - arma::vec ph(m); - arma::vec sp(m); - arma::mat jac(m, np); - ph = phi + pars[0] + coords * pars.tail(np-1); - sp = sin(ph); - jac = arma::join_rows(arma::ones(m), coords); - jac.each_col() %= sp; - return bbar*jac; -} - -/************ - * - * wrap phase into [-pi, pi) - * This is slower than R version - * but export it anyway with a different name. - * -*/ - -//[[Rcpp::export]] - -arma::mat pwrap(const arma::mat& phase) { - unsigned int nr = phase.n_rows; - unsigned int nc = phase.n_cols; - unsigned int ne = phase.n_elem; - arma::mat wphase(nr, nc); - - for (unsigned int i=0; i(np-1, N); - arma::vec t0(np-1); - arma::vec sse = arma::zeros(maxiter); - double abar, bbar; - QMap adata; - adata["coords"] = coords; - arma::vec pars(np); - arma::mat pt_last = join_cols(phases, zcs); - arma::mat pt = pt_last; - double dpt; - int i; - lmreturn lmrets; - - S = arma::join_cols(arma::join_cols(arma::ones(N), arma::cos(phases)), arma::sin(phases)); - Phi = images * pinv(S); - - for (i=0; i -# include "armadillo" +#include "zapm_interface.h" +#include +#include "armadillo" #include #include @@ -296,7 +297,7 @@ mat rzernike_ann(const vec& rho, const double& eps, const int& n, const int& m, //' //' @md // [[Rcpp::export]] -mat zapm(const vec& rho, const vec& theta, const double& eps, const int& maxorder=12) { +mat zapm(const vec& rho, const vec& theta, const double& eps, const int& maxorder) { int j, k, nmax, nz, mmax = maxorder/2; uword nr = rho.size(); @@ -337,7 +338,8 @@ mat zapm(const vec& rho, const vec& theta, const double& eps, const int& maxorde RZ = rzernike_ann(rho, eps, maxorder, 0, xq, qwts); for (int n=0; n<=maxorder; n += 2) { k = (n*n)/4 + n; - zm.col(k) = RZ.col(n/2); + // Dale: normalizers aka m_norms[] are not used here in DFTFringe + zm.col(k) = /* std::sqrt(n+1.)* */RZ.col(n/2); } for (int m=1; m<=mmax; m++) { @@ -348,9 +350,11 @@ mat zapm(const vec& rho, const vec& theta, const double& eps, const int& maxorde j = 0; for (int n=m; n<= nmax; n += 2) { k = ((n+m)*(n+m))/4 + n - m; - zm.col(k) = RZ.col(j) % cosmtheta.col(m-1); + // Dale: normalizers aka m_norms[] are not used here in DFTFringe + zm.col(k) = /* std::sqrt(2.*(n+1.)) * */ RZ.col(j) % cosmtheta.col(m-1); k++; - zm.col(k) = RZ.col(j) % sinmtheta.col(m-1); + // Dale: normalizers aka m_norms[] are not used here in DFTFringe + zm.col(k) = /* std::sqrt(2.*(n+1.)) * */ RZ.col(j) % sinmtheta.col(m-1); j++; } } diff --git a/zernike/zapm_interface.h b/zernike/zapm_interface.h new file mode 100644 index 00000000..7bd39382 --- /dev/null +++ b/zernike/zapm_interface.h @@ -0,0 +1,4 @@ +#include "armadillo" + +arma::mat zapm(const arma::vec &rho, const arma::vec &theta, const double &eps, + const int &maxorder = 12); diff --git a/zernikeprocess.cpp b/zernikeprocess.cpp index 538cae5a..f2f558ef 100644 --- a/zernikeprocess.cpp +++ b/zernikeprocess.cpp @@ -27,6 +27,7 @@ #include "myutils.h" //#include "arbitrarywavefrontdlg.h" #include "userdrawnprofiledlg.h" +#include "zapm_interface.h" std::vector zernEnables; @@ -47,82 +48,6 @@ std::vector zNulls; -arma::mat zapm(const arma::vec& rho, const arma::vec& theta, - const double& eps, const int& maxorder=12) ; -cv::Mat zpmCx(QVector rho, QVector theta, int maxorder) { - - int m, n, n0, mmax = maxorder/2; - int i, imax = rho.size(); - int order, nm, nm1mm1, nm1mp1, nm2m; - int ncol = (mmax+1)*(mmax+1); - double a0; - double cosmtheta[mmax], sinmtheta[mmax]; - cv::Mat_ zm(imax, ncol); - - //do some rudimentary error checking - - //if (rho.length() != theta.length()) stop("Numeric vectors must be same length"); - //if ((maxorder % 2) != 0) stop("maxorder must be even"); - - //good enough - - - for (i=0; i0 - - for (m=1; m <= mmax; m++) { - zm(i, m*m) = rho[i] * zm(i, (m-1)*(m-1)); - } - - // non-symmetric terms - - for (order=4; order<=maxorder; order+=2) { - for (m=order/2-1; m>0; m--) { - n=order-m; - nm = order*order/4 + n - m; - nm1mm1 = (order-2)*(order-2)/4 + n - m; - nm1mp1 = nm - 2; - nm2m = nm1mm1 - 2; - zm(i, nm) = rho[i]*(zm(i, nm1mm1) + zm(i, nm1mp1)) - zm(i, nm2m); - } - - // m=0 (symmetric) term - nm = order*order/4 + order; - nm1mp1 = nm-2; - nm2m = (order-2)*(order-2)/4+order-2; - zm(i, nm) = 2.*rho[i]*zm(i, nm1mp1) - zm(i, nm2m); - } - - // now multiply each column by normalizing factor and cos, sin - - n0 = 1; - for (order=2; order <= maxorder; order+=2) { - for(m=order/2; m>0; m--) { - n=order-m; - a0 = sqrt(2.*(n+1)); - zm(i, n0+1) = a0*sinmtheta[m-1]*zm(i, n0); - zm(i, n0) *= a0*cosmtheta[m-1]; - n0 += 2; - } - n = order; - zm(i, n0) *= sqrt(n+1.); - n0++; - } - } - return zm; -} long fact( int n1, int n2) @@ -1145,63 +1070,6 @@ arma::mat zernikeProcess::zpmC(arma::rowvec rho, arma::rowvec theta, int maxorde return zm; } -arma::mat zernikeProcess::zapmC(const arma::rowvec& rho, const arma::rowvec& theta, const int& maxorder) { - - unsigned int nrow = rho.size(); - int mmax = maxorder/2; - int ncol = (mmax+1)*(mmax+1); - int i, j, m, nj; - double zpnorm = std::sqrt((double) nrow); - arma::mat zm(nrow, ncol), annzm(nrow, ncol); - - //do some rudimentary error checking - -// if (rho.size() != theta.size()) stop("Numeric vectors must be same length"); -// if ((maxorder % 2) != 0) stop("maxorder must be even"); - - //good enough - - zm = zpmC(rho, theta, maxorder); - - // for each azimuthal index m find the column indexes - // of the zp matrix having that value of m. That - // subset is what we need to orthogonalize. - // Note this is done separately for "sine" and "cosine" components. - - nj = maxorder/2 + 1; - for (m=0; m 0) { - qr_econ(Q, R, zm.cols(jsin)); - sgn = sign(R.diag()); - annzm.cols(jsin) = zpnorm * Q * diagmat(sgn); - } - --nj; - } - - // highest order only has one term - - j = mmax*mmax; - - annzm.col(j) = zm.col(j) * zpnorm / norm(zm.col(j)); - annzm.col(j+1) = zm.col(j+1) * zpnorm / norm(zm.col(j+1)); - - return annzm; -} - void zernikeProcess::initGrid(int width, double radius, double cx, double cy, int maxOrder, @@ -1344,8 +1212,7 @@ std::vector zernikeProcess::ZernFitWavefront(wavefront &wf){ } //debug routine to see what is matrix values. -void dumpArma(arma::mat mm, QString title = "", QVector colHeading = QVector(0), - QVector RowLable = QVector(0)){ +void dumpArma(arma::mat mm, QString title, QVector colHeading, QVector RowLable){ arma::mat theMat = mm; QString transposed(" "); diff --git a/zernikeprocess.h b/zernikeprocess.h index 75128cc3..44433a96 100644 --- a/zernikeprocess.h +++ b/zernikeprocess.h @@ -32,6 +32,9 @@ extern std::vector zernEnables; //void gauss_jordan(int n, double* Am, double* Bm); +void dumpArma(arma::mat mm, QString title = "", + QVector colHeading = QVector(0), + QVector RowLable = QVector(0)); class zernikeProcess : public QObject { @@ -62,7 +65,6 @@ class zernikeProcess : public QObject arma::mat rhotheta( int width, double radius, double cx, double cy, 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();