PlanarSplineImplicitization.h

Go to the documentation of this file.

/* This file is part of MoDELib, the Mechanics Of Defects Evolution Library.
 *
 *
 * MoDELib is distributed without any warranty under the
 * GNU General Public License (GPL) v2 <http://www.gnu.org/licenses/>.
 */
 
#ifndef model_PLANARSPLINEIMPLICITIZATION_H_
#define model_PLANARSPLINEIMPLICITIZATION_H_
 
#include <set>
#include <map>
#include <utility> // for std::pair
#include <cfloat> // for FLT_EPSILON
#include <Eigen/Dense>
#include <Eigen/Eigenvalues>
#include <Eigen/Jacobi>
//#include <GeneralizedEigenSolver.h>
 
namespace model
{
    /**********************************************************************************************************************/
    /* PlanarSplineImplicitization<polyDegree> class template *************************************************************/
    /**********************************************************************************************************************/
    template<short unsigned int polyDegree>
    class PlanarSplineImplicitization
    {
        
        enum {polyCoeff=polyDegree+1};
        typedef Eigen::Matrix<double,polyDegree,polyDegree> MatrixPolyDeg;
        
        /*********************************************************************/
        MatrixPolyDeg fill_Mx() const {
            MatrixPolyDeg temp;
            temp << coeffs(1,1), coeffs(1,2), coeffs(1,3),
            /*****/ coeffs(1,2), coeffs(1,3), 0.0,
            /*****/ coeffs(1,3), 0.0,         0.0;
            return temp;
        }
        
        /*********************************************************************/
        MatrixPolyDeg fill_My() const {
            MatrixPolyDeg temp;
            temp << -coeffs(0,1), -coeffs(0,2), -coeffs(0,3),
            /*****/ -coeffs(0,2), -coeffs(0,3),  0.0,
            /*****/ -coeffs(0,3),  0.0,          0.0;
            return temp;
        }
        
        /*********************************************************************/
        MatrixPolyDeg fill_Mc() const {
            MatrixPolyDeg temp;
            temp << coeffs(0,1)*coeffs(1,0)-coeffs(0,0)*coeffs(1,1),  coeffs(0,2)*coeffs(1,0)-coeffs(0,0)*coeffs(1,2),                              coeffs(0,3)*coeffs(1,0)-coeffs(0,0)*coeffs(1,3),
            /*****/ coeffs(0,2)*coeffs(1,0)-coeffs(0,0)*coeffs(1,2),  coeffs(0,3)*coeffs(1,0)-coeffs(0,0)*coeffs(1,3)+coeffs(0,2)*coeffs(1,1)-coeffs(0,1)*coeffs(1,2),  coeffs(0,3)*coeffs(1,1)-coeffs(0,1)*coeffs(1,3),
            /*****/ coeffs(0,3)*coeffs(1,0)-coeffs(0,0)*coeffs(1,3),  coeffs(0,3)*coeffs(1,1)-coeffs(0,1)*coeffs(1,3),                              coeffs(0,3)*coeffs(1,2)-coeffs(0,2)*coeffs(1,3);
            return temp;
        }
        
        
    public:
        
        const Eigen::Matrix<double,2,polyCoeff> coeffs;
        const MatrixPolyDeg Mx;
        const MatrixPolyDeg My;
        const MatrixPolyDeg Mc;
        
        static double compTol;
        static double physTol;
        
        
        /*********************************************************************/
        PlanarSplineImplicitization(const Eigen::Matrix<double,2,polyCoeff>& coeffs_in) : coeffs(coeffs_in),
        /*                                                                           */   Mx(fill_Mx()),
        /*                                                                           */   My(fill_My()),
        /*                                                                           */   Mc(fill_Mc())
        {
            
            
            assert(polyDegree==3 && "GENERALIZE HERE IN THE WAY WE FILL Mx My M0!!!!");
            
        }
        
        //              /*********************************************************************/
        //              std::pair<bool,double> isOnSpline(const Eigen::Matrix<double,2,1>& P, const double& tol=FLT_EPSILON) const {
        //
        //                      std::pair<bool,double> temp=std::make_pair(false,-1.0);
        //                      MatrixPolyDeg M=Mx*P(0)+My*P(1)+M0;
        //                      if (std::fabs(M.determinant())<tol){
        //                              temp.first=true;
        //                              //Eigen::SelfAdjointEigenSolver<MatrixPolyDeg> eigensolver(M);
        //                              //temp.second=...
        //                              assert(0 && "FINISH HERE!!!!");
        //                      }
        //                      return temp;
        //              }
        
        /*********************************************************************/
        template<short unsigned int otherPolyDegree>
        std::set<std::pair<double,double> > intersectWith(const Eigen::Matrix<double,2,otherPolyDegree+1>& otherCoeffs) const
        {
            
            //                  assert((!SplineDegeneracy<dim,polyDegree>::isLine(coeffs,tol)) && "Spline1 is degenerate.");      // make sure that spline1 is not a line
            //                  assert((!SplineDegeneracy<dim,polyDegree>::isLine(otherCoeffs,tol)) && "Spline2 is degenerate."); // make sure that spline2 is not a line
            
            enum {eigenSize=polyDegree*otherPolyDegree};
            typedef Eigen::Matrix< double, eigenSize,eigenSize> EigenSizeMatrixType;
            
            
            std::set<std::pair<double,double> > intersectionParameters;
 
            
            const MatrixPolyDeg Mn(otherCoeffs(0,otherPolyDegree)*Mx+otherCoeffs(1,otherPolyDegree)*My);
            const Eigen::JacobiSVD<MatrixPolyDeg> jsvd(Mn, Eigen::ComputeFullU | Eigen::ComputeFullV);
            const double condNumber(jsvd.singularValues().maxCoeff()/jsvd.singularValues().minCoeff()); // singular values are positive
            
            
            //            std::cout<<"conditionNumber"<<condNumber<<std::endl;
            
            if (condNumber<1.0/FLT_EPSILON)
            {
                // Assemble the matrix C=inv(B)*A
                EigenSizeMatrixType C(Eigen::Matrix<double, eigenSize,eigenSize>::Zero());
                if (otherPolyDegree>1)
                {
                    C.template block<polyDegree*(otherPolyDegree-1),polyDegree*(otherPolyDegree-1)>(0,polyDegree).setIdentity();
                }
                C.template block<polyDegree,polyDegree>(polyDegree*(otherPolyDegree-1),0*polyDegree)= -jsvd.solve(otherCoeffs(0,0)*Mx+otherCoeffs(1,0)*My+Mc);
                for (int p=1; p<otherPolyDegree;++p)
                {
                    C.template block<polyDegree,polyDegree>(polyDegree*(otherPolyDegree-1),p*polyDegree)= -jsvd.solve(otherCoeffs(0,p)*Mx+otherCoeffs(1,p)*My);
                }
                
                //                std::cout<<"C="<<C<<std::endl;
                
                // Compute the eigenvalues and eigenvectors of C
                const Eigen::EigenSolver<EigenSizeMatrixType> es(C);
                
                //               std::cout<<std::setprecision(15)<<std::scientific<<es.eigenvalues()<<std::endl;
                //               std::cout<<std::setprecision(15)<<std::scientific<<es.eigenvectors()<<std::endl;
                //                std::cout<<es.eigenvalues()<<std::endl;
               
#ifdef _MODEL_BENCH_SPLINEINTERSECTION_
 
                EigenSizeMatrixType A(EigenSizeMatrixType::Identity());
                A.template block<polyDegree,polyDegree>(0,0)=-(otherCoeffs(0,0)*Mx+otherCoeffs(1,0)*My+Mc);
                
                EigenSizeMatrixType B(EigenSizeMatrixType::Zero());
                if (otherPolyDegree>1)
                {
                    B.template block<polyDegree*(otherPolyDegree-1),polyDegree*(otherPolyDegree-1)>(polyDegree,0).setIdentity();
                }
                for (int p=0; p<otherPolyDegree;++p)
                {
                    B.template block<polyDegree,polyDegree>(0,p*polyDegree)=otherCoeffs(0,p+1)*Mx+otherCoeffs(1,p+1)*My;
                }
                
                Eigen::GeneralizedEigenSolver<EigenSizeMatrixType> ges(A,B);
                //            ges.compute(A,B);
                
                std::cout<<"geav"<<ges.eigenvalues()<<std::endl<<std::endl;
                
                std::ofstream ABfile;
                ABfile.open ("AB.txt");
                //ABfile << "Writing this to a file.\n";
                ABfile<< std::scientific<<std::setprecision(15)<<A<<std::endl<<B<<std::endl;
                ABfile.close();
                
                std::ofstream coeffFile;
                coeffFile.open ("coeffFile.txt");
                //ABfile << "Writing this to a file.\n";
                coeffFile<< std::scientific<<std::setprecision(15)<<coeffs<<std::endl<<otherCoeffs<<std::endl;
                coeffFile.close();
                
                //            model::GeneralizedEigenSolver<double,eigenSize> GES;
                //            GES.solve(A,B);
                //
                //            //! 3- Clear and fill intersectionParameters if each real eigenvalue s in [0,1] has eigenvator ratio t=v(1)/v(0) also in [0,1]
                //            for (int i=0; i<eigenSize; ++i) {
                //                if (std::fabs(GES.D(i).imag())<tol) {                                 // consider only eigenvectors with 0 imaginary part
                //                    const double t=GES.V(1,i).real()/GES.V(0,i).real();
                //                    const double s=GES.D(i).real();
                //                    if ((s>=0.0 && s<=1.0) && (t>=0.0 && t<=1.0)) { // t and s must be in [0,1]
                //                        Eigen::Matrix<double,2,1> Pt(coeffs.col(0));
                //                        for(int k=1;k<polyDegree+1;++k){
                //                            Pt+=coeffs.col(k)*std::pow(t,k);
                //                        }
                //                        Eigen::Matrix<double,2,1> Ps(otherCoeffs.col(0));
                //                        for(int k=1;k<otherPolyDegree+1;++k){
                //                            Ps+=otherCoeffs.col(k)*std::pow(s,k);
                //                        }
                //                        if((Pt-Ps).norm()<physTol){
                //                            intersectionParameters.insert(std::make_pair(t,s));
                //                        }
                //                    }
                //                }
                //            }
                
                for(int i=0;i<eigenSize;++i)
                {
                    const std::complex<double> s(es.eigenvalues()(i));
                    const std::complex<double> t( es.eigenvectors().col(i)(1)/es.eigenvectors().col(i)(0) ); // the root for the spline defined by "coeffs"
                    std::cout<<t<<", "<<s<<std::endl;
                }
#endif
                
                if (es.info() == Eigen::Success)
                {
                    //                    std::cout<<"eigenvalues="<<std::endl<<std::setprecision(15)<<es.eigenvalues()<<std::endl;
                    //                    std::cout<<"eigenvectors="<<std::endl<<std::setprecision(15)<<es.eigenvectors()<<std::endl;
                    
                    // Sort real eigenvalues
                    std::map<double,const int> sMap;
                    for (int i=0; i<eigenSize; ++i)
                    {
                        if (std::fabs(es.eigenvalues()(i).imag())<compTol)
                        { // consider only eigenvectors with 0 imaginary part
                            //                            const double t( (es.eigenvectors().col(i)(1)/es.eigenvectors().col(i)(0)).real() ); // the root for the spline defined by "coeffs"
                            //                            const double s( es.eigenvalues()(i).real() ); // the root for the spline defined by "otherCoeffs"
                            //                            if ((s>=0.0 && s<=1.0) && (t>=0.0 && t<=1.0))
                            //                            if ((s>=0.0-compTol && s<=1.0+compTol) && (t>=0.0-compTol && t<=1.0+compTol))
                            //                            { // t and s must be in [0,1]
                            sMap.emplace(es.eigenvalues()(i).real(),i);
                            //                            const double s(  ); // the root for the spline defined by "otherCoeffs"
                            
                            //                                Eigen::Matrix<double,2,1> Pt(coeffs.col(0));
                            //                                for(int k=1;k<polyDegree+1;++k)
                            //                                {
                            //                                    Pt+=coeffs.col(k)*std::pow(t,k);
                            //                                }
                            //                                Eigen::Matrix<double,2,1> Ps(otherCoeffs.col(0));
                            //                                for(int k=1;k<otherPolyDegree+1;++k)
                            //                                {
                            //                                    Ps+=otherCoeffs.col(k)*std::pow(s,k);
                            //                                }
                            //                                if((Pt-Ps).norm()<physTol)
                            //                                { // TO DO: COMPARE THE squaredNorm to avoid computing the root
                            //                                    intersectionParameters.insert(std::make_pair(t,s));
                            //                                }
                            //intersectionParameters.insert(std::make_pair(t,s));
                            //                            }
                        }
                    }
                    
                    // Remove multiple eigenvalues
                    for(std::map<double,const int>::const_iterator last=sMap.begin();last!=sMap.end();)
                    {
                        std::map<double,const int>::const_iterator first = last;
                        std::map<double,const int>::const_iterator next = last;
                        next++;
                        last = sMap.upper_bound(first->first+FLT_EPSILON); // jump to the first element grater than current + FLT_EPSILON
                        //                        std::cout<<"first="<<&first<<std::endl;
                        //                        std::cout<<"next="<<&next<<std::endl;
                        //                        std::cout<<"last="<<&last<<std::endl;
                        //std::cout<<"end="<<&sMap.end()<<std::endl;
                        if(last!=next) // more than one equal elements exist
                        {
                            sMap.erase(first,last); // remove all ements between first (included) and last (excluded)
                        }
                    }
                    
                    // Check that simple eigenvlaues are in range
                    for (auto pair : sMap)
                    {
                        const double& s=pair.first;
                        const int& i=pair.second;
                        const double t( (es.eigenvectors().col(i)(1)/es.eigenvectors().col(i)(0)).real() ); // the root for the spline defined by "coeffs"
                        
                        if ((s>=0.0-compTol && s<=1.0+compTol) && (t>=0.0-compTol && t<=1.0+compTol))
                        { // t and s must be in [0,1]
                            Eigen::Matrix<double,2,1> Pt(coeffs.col(0));
                            for(int k=1;k<polyDegree+1;++k)
                            {
                                Pt+=coeffs.col(k)*std::pow(t,k);
                            }
                            Eigen::Matrix<double,2,1> Ps(otherCoeffs.col(0));
                            for(int k=1;k<otherPolyDegree+1;++k)
                            {
                                Ps+=otherCoeffs.col(k)*std::pow(s,k);
                            }
                            if((Pt-Ps).norm()<physTol)
                            { // TO DO: COMPARE THE squaredNorm to avoid computing the root
                                intersectionParameters.insert(std::make_pair(t,s));
                            }
                        }
                    }
                }
            }
            else{
                //    std::cout<<"ILL_CONDITIONED INTERSECTION"<<std::endl;
            }
            
            return intersectionParameters;
        }
        
    };
    
    
    
    // Declare static data members;
    template<short unsigned int polyDegree>
    double PlanarSplineImplicitization<polyDegree>::compTol=FLT_EPSILON;
    
    // Declare static data members;
    template<short unsigned int polyDegree>
    double PlanarSplineImplicitization<polyDegree>::physTol=DBL_MAX;
    
    
    
    /**************************************************************************/
    /**************************************************************************/
}// namespace model
#endif