gsFitting.h

Go to the documentation of this file.

 
#pragma once
 
#include <gsCore/gsForwardDeclarations.h>
#include <gsMSplines/gsMappedBasis.h>
#include <gsMSplines/gsMappedSpline.h>
 
namespace gismo
{
 
template<class T>
class gsFitting
{
public:
    gsFitting()
    {
        m_basis = nullptr;
        m_result= nullptr;
    }
 
    gsFitting(gsMatrix<T> const & param_values,
              gsMatrix<T> const & points,
              gsBasis<T>  & basis);
 
    gsFitting(gsMatrix<T> const & param_values,
              gsMatrix<T> const & points,
              gsVector<index_t>  offset,
              gsMappedBasis<2,T>  & mbasis) ;
 
    virtual ~gsFitting();
 
public:
 
    void compute(T lambda = 0);
 
    void updateGeometry(gsMatrix<T> coefficients, gsMatrix<T> parameters);
 
 
    void initializeGeometry(const gsMatrix<T> & coefficients, const gsMatrix<T> & parameters);
 
    enum tdm_method
    {
        tdm_boundary_pdm, // TDM
        pdm, // PDM
        hybrid_pdm_tdm_boundary_pdm, // HDM with constant weight
        hybrid_error_pdm_tdm_boundary_pdm, // HDM with weight based on point-wise error
        hybrid_curvature_pdm_tdm_boundary_pdm, // HDM with weight based on curvature
    };
 
    void compute_tdm(T lambda, T mu, T sigma, const std::vector<index_t> & interpIdx,
                     tdm_method method = hybrid_curvature_pdm_tdm_boundary_pdm);
 
 
    bool is_corner(gsMatrix<T> & parametric_domain, gsVector<T> & parameter);
 
 
 
    bool is_point_within_cell(const gsMatrix<T>& parameter, const gsMatrix<T>& element);
    bool is_point_within_cell(const T x, const T y, const gsMatrix<T>& element);
 
    bool is_point_inside_support(const gsMatrix<T>& parameter, const gsMatrix<T>& support);
    bool is_point_inside_support(const T x, const T y, const gsMatrix<T>& support);
 
 
    void parameterCorrection(T accuracy = 1e-8,
                             index_t maxIter = 10,
                             T tolOrth = 1e-6);
 
    void parameterProjectionSepBoundary(T accuracy,const std::vector<index_t>& interpIdx);
 
    void parameterCorrectionSepBoundary_pdm(T accuracy, index_t maxIter, const std::vector<index_t>& sepIndex);
 
    void parameterCorrectionSepBoundary_tdm(T accuracy, index_t maxIter, T mu, T sigma, const std::vector<index_t>& sepIndex, tdm_method method = hybrid_curvature_pdm_tdm_boundary_pdm);
 
    void computeErrors();
 
    gsMatrix<T> pointWiseErrors(const gsMatrix<> & parameters,const gsMatrix<> & points);
 
    std::vector<T> computeErrors(const gsMatrix<> & param_values,const gsMatrix<> & points);
 
    void computeMaxNormErrors();
 
    void computeApproxError(T & error, int type = 0) const;
 
    void get_Error(std::vector<T>& errors, int type = 0) const;
 
    T minPointError() const { return m_min_error; }
 
    T maxPointError() const { return m_max_error; }
 
    const std::vector<T> & pointWiseErrors() const
    {
        return m_pointErrors;
    }
 
    size_t numPointsBelow(T threshold) const
    {
        const size_t result=
            std::count_if(m_pointErrors.begin(), m_pointErrors.end(),
                          GS_BIND2ND(std::less<T>(), threshold));
        return result;
    }
 
    void iterativeCompute( T const & tolerance, unsigned const & num_iters = 10);
 
    void applySmoothing(T lambda, gsSparseMatrix<T> & A_mat);
    gsSparseMatrix<T> smoothingMatrix(T lambda) const;
    void assembleSystem(gsSparseMatrix<T>& A_mat, gsMatrix<T>& B);
 
    void compute_normals(const index_t & num_int, const gsMatrix<T> & params_int, gsSparseMatrix<T> & N_int);
 
    gsMatrix<T> fill_pointWiseErrors(const index_t & num_int, T & max_err_int);
 
    gsMatrix<T> principal_curvatures(const gsMatrix<T> & params);
 
    gsMatrix<T> inverse_principal_curvatures(const index_t & num_int, const gsMatrix<T> & params_int);
 
    void blending_weights(const gsSparseMatrix<T> & N_int, const index_t & num_int, const T & mu, const T & sigma,
                          const gsMatrix<T> & params_int, tdm_method method, gsSparseMatrix<T> & NNT);
 
    void assembleSystem(const gsMatrix<T> & points_int, const gsMatrix<T> & params_int,
                        const gsMatrix<T> & points_bdy, const gsMatrix<T> & params_bdy,
                        const index_t & num_basis, const gsSparseMatrix<T> & NNT,
                        gsSparseMatrix<T> & A_tilde, gsMatrix<T> & rhs);
 
public:
 
    gsGeometry<T> * result() const { return m_result; }
 
    const gsMappedSpline<2,T> & mresult() const { return m_mresult; }
 
    const gsBasis<T> & getBasis() const {return *static_cast<const gsBasis<T>*>(m_basis);}
 
    void setBasis(gsBasis<T> & basis) {m_basis=&basis;}
 
    gsMatrix<T> & returnParamValues() {return m_param_values;}
 
    gsMatrix<T> returnPoints() const {return m_points;}
 
    void setConstraints(const gsSparseMatrix<T>& lhs, const gsMatrix<T>& rhs)
    {
        m_constraintsLHS = lhs;
        m_constraintsRHS = rhs;
    }
 
    void setConstraints(const std::vector<index_t>& indices,
                        const std::vector<gsMatrix<T> >& coefs);
 
    void setConstraints(const std::vector<boxSide>& fixedSides);
 
    void setConstraints(const std::vector<boxSide>& fixedSides,
            const std::vector<gsBSpline<T> >& fixedCurves);
    void setConstraints(const std::vector<boxSide>& fixedSides,
            const std::vector<gsGeometry<T> * >& fixedCurves);
 
    void initParametricDomain()
    {
        m_uMin = m_param_values.row(0).minCoeff();
        m_uMax = m_param_values.row(0).maxCoeff();
        m_vMin = m_param_values.row(1).minCoeff();
        m_vMax = m_param_values.row(1).maxCoeff();
 
        gsInfo << "Parametric domain: ["
               << m_uMin << ", " << m_uMax << "] x ["
               << m_vMin << ", " << m_vMax << "]" << std::endl;
    }
 
    T lambda() const {return m_last_lambda;}
 
 
private:
    void extendSystem(gsSparseMatrix<T>& A_mat, gsMatrix<T>& m_B);
 
protected:
 
    void assembleBlockB(const gsMatrix<T>& points,
                        const gsMatrix<T>& params,
                        index_t num_basis,
                        gsSparseMatrix<T>& result) const
    {
        GISMO_UNUSED(points);
        GISMO_UNUSED(num_basis);
        //index_t num_pts = points.rows();
        gsSparseMatrix<T> sparseColloc = m_result->basis().collocationMatrix(params);
        threeOnDiag(sparseColloc, result);
    }
 
    void assembleBlockX(const gsMatrix<T>& points,
                        gsMatrix<T>& result) const
    {
        result.resize(points.rows() * 3, 1);
        result << points.col(0), points.col(1), points.col(2);
    }
 
protected:
 
    //gsOptionList
 
    gsMatrix<T> m_param_values;
 
    gsMatrix<T> m_points;
 
    // Patch offsets
    gsVector<index_t> m_offset;
 
    gsFunctionSet<T> * m_basis;
 
    gsGeometry<T> * m_result;
 
    gsMappedSpline<2,T>  m_mresult;
 
    // All point-wise errors
    std::vector<T> m_pointErrors;
 
    // Interior c1 and c2 curvature, as column of the matrix
    gsMatrix<T> m_pointCurvature;
 
    mutable T m_last_lambda;
 
    T m_max_error;
 
    T m_min_error;
 
    gsSparseMatrix<T> m_constraintsLHS;
 
    gsMatrix<T>       m_constraintsRHS;
 
    T m_uMin, m_uMax, m_vMin, m_vMax;
 
private:
    //void applySmoothing(T lambda, gsMatrix<T> & A_mat);
 
}; // class gsFitting
 
 
#ifdef GISMO_WITH_PYBIND11
 
  void pybind11_init_gsFitting(pybind11::module &m);
 
#endif // GISMO_WITH_PYBIND11
 
 
}// namespace gismo
 
 
#ifndef GISMO_BUILD_LIB
#include GISMO_HPP_HEADER(gsFitting.hpp)
#endif