conjugate tutorial seems to be broken, made some changes but not fixed yet

Former-commit-id: ad4d28f4987eebe06425e341c9fa798a8a865091

include/igl/ConjugateFFSolverData.h

@@ -10,7 +10,9 @@
 #include "igl_inline.h"
 #include <Eigen/Core>
 #include <Eigen/Sparse>
-
+#include <igl/matlab_format.h>
+#include <iostream>
+using namespace std;
 namespace igl 
 {
   // Data class for the Conjugate Frame Field Solver
@@ -122,6 +124,8 @@ IGL_INLINE void igl::ConjugateFFSolverData<DerivedV, DerivedF>::computeCurvature
   kmax = kmax.bottomRows(numF);
   kmin = kmin.bottomRows(numF);
 
+  cerr<<igl::matlab_format(kmax,"kmax")<<endl;
+  cerr<<igl::matlab_format(kmin,"kmin")<<endl;
   //  kmax = dmax3.rowwise().norm();
   //  kmin = dmin3.rowwise().norm();
 

include/igl/conjugate_frame_fields.cpp

@@ -20,46 +20,45 @@ namespace igl {
   public:
     IGL_INLINE ConjugateFFSolver(const ConjugateFFSolverData<DerivedV, DerivedF> &_data,
                                  int _maxIter = 50,
-                                 const typename DerivedV::Scalar &_lambdaOrtho = .1,
-                                 const typename DerivedV::Scalar &_lambdaInit = 100,
-                                 const typename DerivedV::Scalar &_lambdaMultFactor = 1.01,
+                                 const typename DerivedV::Scalar _lambdaOrtho = .1,
+                                 const typename DerivedV::Scalar _lambdaInit = 100,
+                                 const typename DerivedV::Scalar _lambdaMultFactor = 1.01,
                                  bool _doHardConstraints = true);
-    IGL_INLINE bool solve(const Eigen::VectorXi &isConstrained,
-                          const Eigen::PlainObjectBase<DerivedO> &initialSolution,
-                          Eigen::PlainObjectBase<DerivedO> &output,
-                          typename DerivedV::Scalar *lambdaOut = NULL);
-
+    IGL_INLINE typename DerivedV::Scalar solve(const Eigen::VectorXi &isConstrained,
+                                               const Eigen::PlainObjectBase<DerivedO> &initialSolution,
+                                               Eigen::PlainObjectBase<DerivedO> &output);
+    
   private:
-
+    
     const ConjugateFFSolverData<DerivedV, DerivedF> &data;
-
+    
     //polyVF data
     Eigen::Matrix<std::complex<typename DerivedV::Scalar>, Eigen::Dynamic, 1> Acoeff, Bcoeff;
     Eigen::Matrix<typename DerivedV::Scalar, Eigen::Dynamic, 2> pvU, pvV;
     typename DerivedV::Scalar lambda;
-
+    
     //parameters
     typename DerivedV::Scalar lambdaOrtho;
     typename DerivedV::Scalar lambdaInit,lambdaMultFactor;
     int maxIter;
     bool doHardConstraints;
-
+    
     IGL_INLINE void localStep();
     IGL_INLINE void getPolyCoeffsForLocalSolve(const Eigen::Matrix<typename DerivedV::Scalar, 4, 1> &s,
                                                const Eigen::Matrix<typename DerivedV::Scalar, 4, 1> &z,
                                                Eigen::Matrix<typename DerivedV::Scalar, Eigen::Dynamic, 1> &polyCoeff);
-
+    
     IGL_INLINE void globalStep(const Eigen::Matrix<int, Eigen::Dynamic, 1>  &isConstrained,
                                const Eigen::Matrix<std::complex<typename DerivedV::Scalar>, Eigen::Dynamic, 1>  &Ak,
                                const Eigen::Matrix<std::complex<typename DerivedV::Scalar>, Eigen::Dynamic, 1>  &Bk);
     IGL_INLINE void minQuadWithKnownMini(const Eigen::SparseMatrix<std::complex<typename DerivedV::Scalar> > &Q,
-                         const Eigen::SparseMatrix<std::complex<typename DerivedV::Scalar> > &f,
-                         const Eigen::VectorXi isConstrained,
-                         const Eigen::Matrix<std::complex<typename DerivedV::Scalar>, Eigen::Dynamic, 1> &xknown,
+                                         const Eigen::SparseMatrix<std::complex<typename DerivedV::Scalar> > &f,
+                                         const Eigen::VectorXi isConstrained,
+                                         const Eigen::Matrix<std::complex<typename DerivedV::Scalar>, Eigen::Dynamic, 1> &xknown,
                                          Eigen::Matrix<std::complex<typename DerivedV::Scalar>, Eigen::Dynamic, 1> &x);
     IGL_INLINE void setFieldFromCoefficients();
     IGL_INLINE void setCoefficientsFromField();
-
+    
   };
 }
 
@@ -69,9 +68,9 @@ template <typename DerivedV, typename DerivedF, typename DerivedO>
 IGL_INLINE igl::ConjugateFFSolver<DerivedV, DerivedF, DerivedO>::
 ConjugateFFSolver(const ConjugateFFSolverData<DerivedV, DerivedF> &_data,
                   int _maxIter,
-                  const typename DerivedV::Scalar &_lambdaOrtho,
-                  const typename DerivedV::Scalar &_lambdaInit,
-                  const typename DerivedV::Scalar &_lambdaMultFactor,
+                  const typename DerivedV::Scalar _lambdaOrtho,
+                  const typename DerivedV::Scalar _lambdaInit,
+                  const typename DerivedV::Scalar _lambdaMultFactor,
                   bool _doHardConstraints):
 data(_data),
 lambdaOrtho(_lambdaOrtho),
@@ -102,7 +101,7 @@ getPolyCoeffsForLocalSolve(const Eigen::Matrix<typename DerivedV::Scalar, 4, 1>
   typename DerivedV::Scalar z1 = z(1);
   typename DerivedV::Scalar z2 = z(2);
   typename DerivedV::Scalar z3 = z(3);
-
+  
   polyCoeff.resize(7,1);
   polyCoeff(0) =  s0*s0* s1*s1* s2*s2* s3* z3*z3 +  s0*s0* s1*s1* s2* s3*s3* z2*z2 +  s0*s0* s1* s2*s2* s3*s3* z1*z1 +  s0* s1*s1* s2*s2* s3*s3* z0*z0 ;
   polyCoeff(1) = 2* s0*s0* s1*s1* s2* s3* z2*z2 + 2* s0*s0* s1*s1* s2* s3* z3*z3 + 2* s0*s0* s1* s2*s2* s3* z1*z1 + 2* s0*s0* s1* s2*s2* s3* z3*z3 + 2* s0*s0* s1* s2* s3*s3* z1*z1 + 2* s0*s0* s1* s2* s3*s3* z2*z2 + 2* s0* s1*s1* s2*s2* s3* z0*z0 + 2* s0* s1*s1* s2*s2* s3* z3*z3 + 2* s0* s1*s1* s2* s3*s3* z0*z0 + 2* s0* s1*s1* s2* s3*s3* z2*z2 + 2* s0* s1* s2*s2* s3*s3* z0*z0 + 2* s0* s1* s2*s2* s3*s3* z1*z1 ;
@@ -111,7 +110,7 @@ getPolyCoeffsForLocalSolve(const Eigen::Matrix<typename DerivedV::Scalar, 4, 1>
   polyCoeff(4) =  s0*s0* s1* z1*z1 +  s0*s0* s2* z2*z2 +  s0*s0* s3* z3*z3 +  s0* s1*s1* z0*z0 + 4* s0* s1* s2* z0*z0 + 4* s0* s1* s2* z1*z1 + 4* s0* s1* s2* z2*z2 + 4* s0* s1* s3* z0*z0 + 4* s0* s1* s3* z1*z1 + 4* s0* s1* s3* z3*z3 +  s0* s2*s2* z0*z0 + 4* s0* s2* s3* z0*z0 + 4* s0* s2* s3* z2*z2 + 4* s0* s2* s3* z3*z3 +  s0* s3*s3* z0*z0 +  s1*s1* s2* z2*z2 +  s1*s1* s3* z3*z3 +  s1* s2*s2* z1*z1 + 4* s1* s2* s3* z1*z1 + 4* s1* s2* s3* z2*z2 + 4* s1* s2* s3* z3*z3 +  s1* s3*s3* z1*z1 +  s2*s2* s3* z3*z3 +  s2* s3*s3* z2*z2;
   polyCoeff(5) = 2* s0* s1* z0*z0 + 2* s0* s1* z1*z1 + 2* s0* s2* z0*z0 + 2* s0* s2* z2*z2 + 2* s0* s3* z0*z0 + 2* s0* s3* z3*z3 + 2* s1* s2* z1*z1 + 2* s1* s2* z2*z2 + 2* s1* s3* z1*z1 + 2* s1* s3* z3*z3 + 2* s2* s3* z2*z2 + 2* s2* s3* z3*z3 ;
   polyCoeff(6) =  s0* z0*z0 +  s1* z1*z1 +  s2* z2*z2 +  s3* z3*z3;
-
+  
 }
 
 
@@ -124,12 +123,12 @@ localStep()
     Eigen::Matrix<typename DerivedV::Scalar, 4, 1> xproj; xproj << pvU.row(j).transpose(),pvV.row(j).transpose();
     Eigen::Matrix<typename DerivedV::Scalar, 4, 1> z = data.UH[j].transpose()*xproj;
     Eigen::Matrix<typename DerivedV::Scalar, 4, 1> x;
-
+    
     Eigen::Matrix<typename DerivedV::Scalar, Eigen::Dynamic, 1> polyCoeff;
     getPolyCoeffsForLocalSolve(data.s[j], z, polyCoeff);
     Eigen::Matrix<std::complex<typename DerivedV::Scalar>, Eigen::Dynamic, 1> roots;
     igl::polyRoots<typename DerivedV::Scalar, typename DerivedV::Scalar> (polyCoeff, roots );
-
+    
     //  find closest real root to xproj
     typename DerivedV::Scalar minDist = 1e10;
     for (int i =0; i< 6; ++i)
@@ -144,9 +143,9 @@ localStep()
         minDist = dist;
         x = candidate;
       }
-
+      
     }
-
+    
     pvU.row(j) << x(0),x(1);
     pvV.row(j) << x(2),x(3);
   }
@@ -174,15 +173,15 @@ globalStep(const Eigen::Matrix<int, Eigen::Dynamic, 1>  &isConstrained,
            const Eigen::Matrix<std::complex<typename DerivedV::Scalar>, Eigen::Dynamic, 1>  &Bk)
 {
   setCoefficientsFromField();
-
+  
   Eigen::SparseMatrix<std::complex<typename DerivedV::Scalar> > I;
   igl::speye(data.numF, data.numF, I);
   Eigen::SparseMatrix<std::complex<typename DerivedV::Scalar> > QA = data.DDA+lambda*data.planarityWeight+lambdaOrtho*I;
   Eigen::SparseMatrix<std::complex<typename DerivedV::Scalar> > fA = (-2*lambda*data.planarityWeight*Acoeff).sparseView();
-
+  
   Eigen::SparseMatrix<std::complex<typename DerivedV::Scalar> > QB = data.DDB+lambda*data.planarityWeight;
   Eigen::SparseMatrix<std::complex<typename DerivedV::Scalar> > fB = (-2*lambda*data.planarityWeight*Bcoeff).sparseView();
-
+  
   if(doHardConstraints)
   {
     minQuadWithKnownMini(QA, fA, isConstrained, Ak, Acoeff);
@@ -196,7 +195,7 @@ globalStep(const Eigen::Matrix<int, Eigen::Dynamic, 1>  &isConstrained,
     minQuadWithKnownMini(QB, fB, isknown_, xknown_, Bcoeff);
   }
   setFieldFromCoefficients();
-
+  
 }
 
 
@@ -210,10 +209,10 @@ setFieldFromCoefficients()
     //    matlab code from roots (given there are no trailing zeros in the polynomial coefficients)
     Eigen::Matrix<std::complex<typename DerivedV::Scalar>, Eigen::Dynamic, 1> polyCoeff(5,1);
     polyCoeff<<1., 0., -Acoeff(i), 0., Bcoeff(i);
-
+    
     Eigen::Matrix<std::complex<typename DerivedV::Scalar>, Eigen::Dynamic, 1> roots;
     polyRoots<std::complex<typename DerivedV::Scalar>>(polyCoeff,roots);
-
+    
     std::complex<typename DerivedV::Scalar> u = roots[0];
     int maxi = -1;
     float maxd = -1;
@@ -230,7 +229,7 @@ setFieldFromCoefficients()
     pvU(i,0) = real(u); pvU(i,1) = imag(u);
     pvV(i,0) = real(v); pvV(i,1) = imag(v);
   }
-
+  
 }
 
 template<typename DerivedV, typename DerivedF, typename DerivedO>
@@ -242,11 +241,11 @@ minQuadWithKnownMini(const Eigen::SparseMatrix<std::complex<typename DerivedV::S
                      Eigen::Matrix<std::complex<typename DerivedV::Scalar>, Eigen::Dynamic, 1> &x)
 {
   int N = Q.rows();
-
+  
   int nc = xknown.rows();
   Eigen::VectorXi known; known.setZero(nc,1);
   Eigen::VectorXi unknown; unknown.setZero(N-nc,1);
-
+  
   int indk = 0, indu = 0;
   for (int i = 0; i<N; ++i)
     if (isConstrained[i])
@@ -259,21 +258,21 @@ minQuadWithKnownMini(const Eigen::SparseMatrix<std::complex<typename DerivedV::S
       unknown[indu] = i;
       indu++;
     }
-
+  
   Eigen::SparseMatrix<std::complex<typename DerivedV::Scalar>> Quu, Quk;
-
+  
   igl::slice(Q,unknown, unknown, Quu);
   igl::slice(Q,unknown, known, Quk);
-
-
+  
+  
   std::vector<typename Eigen::Triplet<std::complex<typename DerivedV::Scalar> > > tripletList;
-
+  
   Eigen::SparseMatrix<std::complex<typename DerivedV::Scalar> > fu(N-nc,1);
-
+  
   igl::slice(f,unknown, Eigen::VectorXi::Zero(1,1), fu);
-
+  
   Eigen::SparseMatrix<std::complex<typename DerivedV::Scalar> > rhs = (Quk*xknown).sparseView()+.5*fu;
-
+  
   Eigen::SparseLU< Eigen::SparseMatrix<std::complex<typename DerivedV::Scalar>>> solver;
   solver.compute(-Quu);
   if(solver.info()!=Eigen::Success)
@@ -287,7 +286,7 @@ minQuadWithKnownMini(const Eigen::SparseMatrix<std::complex<typename DerivedV::S
     std::cerr<<"Solving failed!"<<std::endl;
     return;
   }
-
+  
   indk = 0, indu = 0;
   x.setZero(N,1);
   for (int i = 0; i<N; ++i)
@@ -295,21 +294,20 @@ minQuadWithKnownMini(const Eigen::SparseMatrix<std::complex<typename DerivedV::S
       x[i] = xknown[indk++];
     else
       x[i] = b.coeff(indu++,0);
-
+  
 }
 
 
 template<typename DerivedV, typename DerivedF, typename DerivedO>
-IGL_INLINE bool igl::ConjugateFFSolver<DerivedV, DerivedF, DerivedO>::
+IGL_INLINE typename DerivedV::Scalar igl::ConjugateFFSolver<DerivedV, DerivedF, DerivedO>::
 solve(const Eigen::VectorXi &isConstrained,
       const Eigen::PlainObjectBase<DerivedO> &initialSolution,
-      Eigen::PlainObjectBase<DerivedO> &output,
-      typename DerivedV::Scalar *lambdaOut)
+      Eigen::PlainObjectBase<DerivedO> &output)
 {
   int numConstrained = isConstrained.sum();
   // coefficient values
   Eigen::Matrix<std::complex<typename DerivedV::Scalar>, Eigen::Dynamic, 1> Ak, Bk;
-
+  
   pvU.resize(data.numF,2);
   pvV.resize(data.numF,2);
   for (int fi = 0; fi <data.numF; ++fi)
@@ -334,54 +332,54 @@ solve(const Eigen::VectorXi &isConstrained,
       ind ++;
     }
   }
-
-
-
+  
+  
+  
   typename DerivedV::Scalar smoothnessValue;
   Eigen::Matrix<typename DerivedV::Scalar, Eigen::Dynamic, 1> conjValues;
   typename DerivedV::Scalar meanConj;
   typename DerivedV::Scalar maxConj;
-
+  
   data.evaluateConjugacy(pvU, pvV, conjValues);
   meanConj = conjValues.cwiseAbs().mean();
   maxConj = conjValues.cwiseAbs().maxCoeff();
   printf("Initial max non-conjugacy: %.5g\n",maxConj);
-
+  
   smoothnessValue = (Acoeff.adjoint()*data.DDA*Acoeff + Bcoeff.adjoint()*data.DDB*Bcoeff).real()[0];
   printf("\n\nInitial smoothness: %.5g\n",smoothnessValue);
-
+  
   lambda = lambdaInit;
-
+  
   bool doit = false;
   for (int iter = 0; iter<maxIter; ++iter)
   {
     printf("\n\n--- Iteration %d ---\n",iter);
-
+    
     typename DerivedV::Scalar oldMeanConj = meanConj;
-
+    
     localStep();
     globalStep(isConstrained, Ak, Bk);
-
-
+    
+    
     smoothnessValue = (Acoeff.adjoint()*data.DDA*Acoeff + Bcoeff.adjoint()*data.DDB*Bcoeff).real()[0];
-
+    
     printf("Smoothness: %.5g\n",smoothnessValue);
-
+    
     data.evaluateConjugacy(pvU, pvV, conjValues);
     meanConj = conjValues.cwiseAbs().mean();
     maxConj = conjValues.cwiseAbs().maxCoeff();
     printf("Mean/Max non-conjugacy: %.5g, %.5g\n",meanConj,maxConj);
     typename DerivedV::Scalar diffMeanConj = fabs(oldMeanConj-meanConj);
-
+    
     if (diffMeanConj<1e-4)
       doit = true;
-
+    
     if (doit)
       lambda = lambda*lambdaMultFactor;
     printf(" %d %.5g %.5g\n",iter, smoothnessValue,maxConj);
-
+    
   }
-
+  
   output.setZero(data.numF,6);
   for (int fi=0; fi<data.numF; ++fi)
   {
@@ -390,26 +388,23 @@ solve(const Eigen::VectorXi &isConstrained,
     output.block(fi,0, 1, 3) = pvU(fi,0)*b1 + pvU(fi,1)*b2;
     output.block(fi,3, 1, 3) = pvV(fi,0)*b1 + pvV(fi,1)*b2;
   }
-
-  if (lambdaOut)
-    *lambdaOut = lambda;
-
-  return true;
+  
+  return lambda;
 }
 
 
 
 template <typename DerivedV, typename DerivedF, typename DerivedO>
 IGL_INLINE void igl::conjugate_frame_fields(const Eigen::PlainObjectBase<DerivedV> &V,
-                                            const Eigen::PlainObjectBase<DerivedF> &F,
-                                            const Eigen::VectorXi &isConstrained,
-                                            const Eigen::PlainObjectBase<DerivedO> &initialSolution,
-                                            Eigen::PlainObjectBase<DerivedO> &output,
-                                            int maxIter,
-                                            const typename DerivedV::Scalar &lambdaOrtho,
-                                            const typename DerivedV::Scalar &lambdaInit,
-                                            const typename DerivedV::Scalar &lambdaMultFactor,
-                                            bool doHardConstraints)
+                                       const Eigen::PlainObjectBase<DerivedF> &F,
+                                       const Eigen::VectorXi &isConstrained,
+                                       const Eigen::PlainObjectBase<DerivedO> &initialSolution,
+                                       Eigen::PlainObjectBase<DerivedO> &output,
+                                       int maxIter,
+                                       const typename DerivedV::Scalar lambdaOrtho,
+                                       const typename DerivedV::Scalar lambdaInit,
+                                       const typename DerivedV::Scalar lambdaMultFactor,
+                                       bool doHardConstraints)
 {
   igl::ConjugateFFSolverData<DerivedV, DerivedF> csdata(V, F);
   igl::ConjugateFFSolver<DerivedV, DerivedF, DerivedO> cs(csdata, maxIter, lambdaOrtho, lambdaInit, lambdaMultFactor, doHardConstraints);
@@ -417,22 +412,22 @@ IGL_INLINE void igl::conjugate_frame_fields(const Eigen::PlainObjectBase<Derived
 }
 
 template <typename DerivedV, typename DerivedF, typename DerivedO>
-IGL_INLINE void igl::conjugate_frame_fields(const igl::ConjugateFFSolverData<DerivedV, DerivedF> &csdata,
-                                            const Eigen::VectorXi &isConstrained,
-                                            const Eigen::PlainObjectBase<DerivedO> &initialSolution,
-                                            Eigen::PlainObjectBase<DerivedO> &output,
-                                            int maxIter,
-                                            const typename DerivedV::Scalar &lambdaOrtho,
-                                            const typename DerivedV::Scalar &lambdaInit,
-                                            const typename DerivedV::Scalar &lambdaMultFactor,
-                                            bool doHardConstraints,
-                                            typename DerivedV::Scalar *lambdaOut)
+IGL_INLINE typename DerivedV::Scalar igl::conjugate_frame_fields(const igl::ConjugateFFSolverData<DerivedV, DerivedF> &csdata,
+                                                                 const Eigen::VectorXi &isConstrained,
+                                                                 const Eigen::PlainObjectBase<DerivedO> &initialSolution,
+                                                                 Eigen::PlainObjectBase<DerivedO> &output,
+                                                                 int maxIter,
+                                                                 const typename DerivedV::Scalar lambdaOrtho,
+                                                                 const typename DerivedV::Scalar lambdaInit,
+                                                                 const typename DerivedV::Scalar lambdaMultFactor,
+                                                                 bool doHardConstraints)
 {
   igl::ConjugateFFSolver<DerivedV, DerivedF, DerivedO> cs(csdata, maxIter, lambdaOrtho, lambdaInit, lambdaMultFactor, doHardConstraints);
-  cs.solve(isConstrained, initialSolution, output, lambdaOut);
+  typename DerivedV::Scalar lambdaOut = cs.solve(isConstrained, initialSolution, output);
+  return lambdaOut;
 }
 
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template specialization
-template void igl::conjugate_frame_fields<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1> >(igl::ConjugateFFSolverData<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::Matrix<int, -1, 1, 0, -1, 1> const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&, int, Eigen::Matrix<double, -1, -1, 0, -1, -1>::Scalar const&, Eigen::Matrix<double, -1, -1, 0, -1, -1>::Scalar const&, Eigen::Matrix<double, -1, -1, 0, -1, -1>::Scalar const&, bool, Eigen::Matrix<double, -1, -1, 0, -1, -1>::Scalar*);
+template Eigen::Matrix<double, -1, -1, 0, -1, -1>::Scalar igl::conjugate_frame_fields<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1> >(igl::ConjugateFFSolverData<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::Matrix<int, -1, 1, 0, -1, 1> const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&, int, Eigen::Matrix<double, -1, -1, 0, -1, -1>::Scalar, Eigen::Matrix<double, -1, -1, 0, -1, -1>::Scalar, Eigen::Matrix<double, -1, -1, 0, -1, -1>::Scalar, bool);
 #endif

include/igl/conjugate_frame_fields.h

@@ -17,17 +17,7 @@ namespace igl {
   //todo
   // TODO: isConstrained should become a list of indices for consistency with
   //       n_polyvector
-  /// Given 2 vectors centered on origin calculate the rotation matrix from first to the second
-
-  // Inputs:
-  //   v0, v1         the two #3 by 1 vectors
-  //   normalized     boolean, if false, then the vectors are normalized prior to the calculation
-  // Output:
-  //                  3 by 3 rotation matrix that takes v0 to v1
-  //
-  template <typename DerivedV, typename DerivedF>
-  class ConjugateFFSolverData;
-
+  
   template <typename DerivedV, typename DerivedF, typename DerivedO>
   IGL_INLINE void conjugate_frame_fields(const Eigen::PlainObjectBase<DerivedV> &V,
                                          const Eigen::PlainObjectBase<DerivedF> &F,
@@ -35,23 +25,22 @@ namespace igl {
                                          const Eigen::PlainObjectBase<DerivedO> &initialSolution,
                                          Eigen::PlainObjectBase<DerivedO> &output,
                                          int _maxIter = 50,
-                                         const typename DerivedV::Scalar &_lambdaOrtho = .1,
-                                         const typename DerivedV::Scalar &_lambdaInit = 100,
-                                         const typename DerivedV::Scalar &_lambdaMultFactor = 1.01,
+                                         const typename DerivedV::Scalar _lambdaOrtho = .1,
+                                         const typename DerivedV::Scalar _lambdaInit = 100,
+                                         const typename DerivedV::Scalar _lambdaMultFactor = 1.01,
                                          bool _doHardConstraints = true);
-
+  
   template <typename DerivedV, typename DerivedF, typename DerivedO>
-  IGL_INLINE void conjugate_frame_fields(const ConjugateFFSolverData<DerivedV, DerivedF> &csdata,
-                                         const Eigen::VectorXi &isConstrained,
-                                         const Eigen::PlainObjectBase<DerivedO> &initialSolution,
-                                         Eigen::PlainObjectBase<DerivedO> &output,
-                                         int _maxIter = 50,
-                                         const typename DerivedV::Scalar &_lambdaOrtho = .1,
-                                         const typename DerivedV::Scalar &_lambdaInit = 100,
-                                         const typename DerivedV::Scalar &_lambdaMultFactor = 1.01,
-                                         bool _doHardConstraints = true,
-                                         typename DerivedV::Scalar *lambdaOut = NULL);
-
+  IGL_INLINE typename DerivedV::Scalar conjugate_frame_fields(const ConjugateFFSolverData<DerivedV, DerivedF> &csdata,
+                                                              const Eigen::VectorXi &isConstrained,
+                                                              const Eigen::PlainObjectBase<DerivedO> &initialSolution,
+                                                              Eigen::PlainObjectBase<DerivedO> &output,
+                                                              int _maxIter = 50,
+                                                              const typename DerivedV::Scalar _lambdaOrtho = .1,
+                                                              const typename DerivedV::Scalar _lambdaInit = 100,
+                                                              const typename DerivedV::Scalar _lambdaMultFactor = 1.01,
+                                                              bool _doHardConstraints = true);
+  
 };
 
 

tutorial/508_ConjugateField/main.cpp

@@ -33,8 +33,6 @@ Eigen::MatrixXd conjugate_pvf;
 Eigen::VectorXd conjugacy_s;
 Eigen::VectorXd conjugacy_c;
 
-igl::ConjugateFFSolverData<Eigen::MatrixXd, Eigen::MatrixXi> *csdata;
-
 
 bool key_down(igl::viewer::Viewer& viewer, unsigned char key, int modifier)
 {
@@ -143,7 +141,7 @@ int main(int argc, char *argv[])
   igl::n_polyvector(V, F, b, bc, smooth_pvf);
 
   // Initialize conjugate field with smooth field
-  csdata = new igl::ConjugateFFSolverData<Eigen::MatrixXd,Eigen::MatrixXi>(V,F);
+  igl::ConjugateFFSolverData<Eigen::MatrixXd, Eigen::MatrixXi> csdata(V,F);
   conjugate_pvf = smooth_pvf;
 
 
@@ -153,12 +151,11 @@ int main(int argc, char *argv[])
   double lambdaInit = 100;
   double lambdaMultFactor = 1.01;
   bool doHardConstraints = true;
-  double lambdaOut;
   VectorXi isConstrained = VectorXi::Constant(F.rows(),0);
   for (unsigned i=0; i<b.size(); ++i)
     isConstrained(b(i)) = 1;
 
-  igl::conjugate_frame_fields(*csdata, isConstrained, conjugate_pvf, conjugate_pvf, conjIter, lambdaOrtho, lambdaInit, lambdaMultFactor, doHardConstraints, &lambdaOut);
+  double lambdaOut = igl::conjugate_frame_fields(csdata, isConstrained, conjugate_pvf, conjugate_pvf, conjIter, lambdaOrtho, lambdaInit, lambdaMultFactor, doHardConstraints);
 
   // local representations of field vectors
   Eigen::Matrix<double, Eigen::Dynamic, 2> pvU, pvV;
@@ -168,13 +165,13 @@ int main(int argc, char *argv[])
   const Eigen::MatrixXd &Vs = smooth_pvf.rightCols(3);
   pvU << igl::dot_row(Us,B1), igl::dot_row(Us,B2);
   pvV << igl::dot_row(Vs,B1), igl::dot_row(Vs,B2);
-  csdata->evaluateConjugacy(pvU, pvV, conjugacy_s);
+  csdata.evaluateConjugacy(pvU, pvV, conjugacy_s);
   //conjugate
   const Eigen::MatrixXd &Uc = conjugate_pvf.leftCols(3);
   const Eigen::MatrixXd &Vc = conjugate_pvf.rightCols(3);
   pvU << igl::dot_row(Uc,B1), igl::dot_row(Uc,B2);
   pvV << igl::dot_row(Vc,B1), igl::dot_row(Vc,B2);
-  csdata->evaluateConjugacy(pvU, pvV, conjugacy_c);
+  csdata.evaluateConjugacy(pvU, pvV, conjugacy_c);
   // Launch the viewer
   igl::viewer::Viewer viewer;
   viewer.core.invert_normals = true;