gsBSpline.h

Go to the documentation of this file.

#pragma once

#include <gsCore/gsLinearAlgebra.h>
#include <gsCore/gsGeometry.h>

#include <gsNurbs/gsBSplineBasis.h>
#include <gsNurbs/gsBoehm.h>
#include <gsNurbs/gsBSplineSolver.h>

#include <gsUtils/gsPointGrid.h>

namespace gismo {
namespace internal {

template<class T>
struct gsCurveIntersectionResult;

} // namespace internal
} // namespace gismo

namespace gismo
{
    
template<class T>
class gsBSpline : public gsGeoTraits<1,T>::GeometryBase
{
public: 
    typedef gsKnotVector<T> KnotVectorType;

    typedef typename gsGeoTraits<1,T>::GeometryBase Base;

    typedef gsBSplineBasis<T> Basis;

    typedef memory::shared_ptr< gsBSpline > Ptr;

    typedef memory::unique_ptr< gsBSpline > uPtr;
    
public:
    
    gsBSpline() { }

    // enable swap
    using gsGeometry<T>::swap;
    
    gsBSpline( const Basis & basis, gsMatrix<T> coefs )
    : Base( basis, give(coefs) )
    {
        if( this->basis().isPeriodic() )
            this->basis().expandCoefs(m_coefs);
    }
    
    gsBSpline(KnotVectorType KV, gsMatrix<T> coefs, bool periodic = false )
    {
        this->m_basis = new Basis(give(KV));
        m_coefs.swap(coefs);

        if( periodic )
        {
            const index_t sz = this->basis().size();

            this->basis().setPeriodic();

            if ( m_coefs.rows() == sz )
            {
                this->basis().trimCoefs(m_coefs);
            }
            else if  ( m_coefs.rows() == this->basis().size() )
            {
                this->basis().expandCoefs(m_coefs);
            }
            else
            {
                GISMO_ERROR("Wrong number of coefficients for periodic basis.");
            }
        }
        else // non-periodic
        {
            if( this->m_coefs.rows() + KV.degree() + 1 != static_cast<int>( this->knots().size() ) )
                gsWarn << "gsBSpline Warning: #Knots="<< this->knots().size()<< ", #coefs="<< this->m_coefs.rows() <<"\n";
        }
    }
    

    gsBSpline(T u0, T u1, unsigned interior, int degree,
              gsMatrix<T> coefs, unsigned mult_interior=1, bool periodic = false)
    {
        this->m_basis = new Basis(u0, u1, interior, degree, mult_interior, periodic );
        if( periodic )
        {
            GISMO_ASSERT( this->basis().numCrossingFunctions() == 1,
                          "Knot vector with clamped knots should lead to m_periodic == 1.");

            this->m_coefs = this->basis().perCoefs(coefs);
            GISMO_ASSERT( this->m_coefs.rows() >= this->basis().size(),
                          "You should give at least as many coefficients as the size of the basis after converting to periodic.");
            if( this->m_coefs.rows() > this->basis().size() )
                gsWarn << "Some of the last coefficients are going to be lost during conversion to periodic basis.\n";
        }
        else // non-periodic
        {
            this->m_coefs.swap(coefs);
            GISMO_ASSERT( this->m_coefs.rows() == this->basis().size(),
                          "Number of coefficients does not match the size of the basis.");
        }
    }

    GISMO_CLONE_FUNCTION(gsBSpline)
    
    GISMO_BASIS_ACCESSORS    
    
public:
    
    std::ostream &print(std::ostream &os) const
    {
        os << "BSpline curve "<< "of degree "<< this->basis().degree()<< ", "<<  this->basis().knots() <<".\n";
        os << "with control points "<< this->m_coefs.row(0)<<" ... "<<this->m_coefs.bottomRows(1) << ".\n";
        if( this->basis().isPeriodic() )
            os << "Periodic with overlay " << this->basis().numCrossingFunctions() << ".\n";
        return os;
    }
    
    T domainStart() const { return *this->basis().knots().domainBegin(); }
    
    T domainEnd() const { return *this->basis().knots().domainEnd(); }
    
    KnotVectorType & knots(const int i = 0) 
    {
        GISMO_UNUSED(i);
        GISMO_ASSERT( i==0, "Requested knots of invalid direction "<< i );
        return this->basis().knots();
    }

    const KnotVectorType & knots(const int i = 0) const 
    {
        GISMO_UNUSED(i);
        GISMO_ASSERT( i==0, "Requested knots of invalid direction "<< i );
        return this->basis().knots(); 
    }

    short_t degree(short_t i = 0) const
    {
        GISMO_UNUSED(i);
        GISMO_ASSERT( i==0, "Requested knots of invalid direction "<< i );
        return this->basis().degree(); 
    }
    
    // compatible curves: same degree, same first/last p+1 knots
    void isCompatible( gsGeometry<T> * )
    { GISMO_NO_IMPLEMENTATION }
    
    // compatible curves: same degree, same first/last p+1 knots
    void makeCompatible( gsGeometry<T> * )
    { GISMO_NO_IMPLEMENTATION }


    void merge( gsGeometry<T> * otherG );

    gsBSpline<T> segmentFromTo(T u0, T u1, T tolerance=1e-15) const;

    void splitAt(T u0, gsBSpline<T>& left,  gsBSpline<T>& right, T tolerance=1e-15) const;

    gsMultiPatch<T> toBezier(T tolerance=1e-15) const;

    std::vector<internal::gsCurveIntersectionResult<T>> intersect(const gsBSpline<T>& other, T tolerance = 1e-5, T curvatureTolerance=1+1e-6) const;

    T pseudoCurvature() const;

    void insertKnot( T knot, index_t i = 1);
private:
    // Resolve hidden overload w.r.t. gsGeometry
    virtual void insertKnot( T knot, index_t dir, index_t i = 1);
public:

    template <class It>
    void insertKnots( It inBegin, It inEnd)
    {
        if( this->basis().isPeriodic() )
        {
            // We assume we got valid (i.e., non-NULL) iterators; I don't think we have a reasonable way to test it in GISMO_ASSERT.

            GISMO_ASSERT( (*inBegin > *(this->knots().begin()))
                          && (*(inEnd-1) < *(this->knots().end()-1)),
                          "Please, ask me to insert knots inside the knot interval only." );

            std::vector<T> newKnots(inBegin,inEnd);

            // It can happen that user would ask us to insert knots outside the active range.
            // Should it happen, we shift the values and then sort the knot vector to remain non-decreasing.

            T activeLength = this->basis()._activeLength();
            T blue1 = *(this->basis().knots().begin()  + this->degree() );
            T blue2 = *(this->basis().knots().end()    - this->degree() );
            typename std::vector<T>::iterator begin = newKnots.begin();
            typename std::vector<T>::iterator end   = newKnots.end();
            for( typename std::vector<T>::iterator it = begin; it != end; ++it )
            {
                if( *it < blue1 )
                    *it += activeLength;
                else if( *it > blue2 )
                    *it -= activeLength;
                else if( (*it == blue1) || (*it == blue2) )
                {
                    // I would like to erase it and go on but for that the iterator is not enough.
                    // ANGELOS, wouldn't it be better to put a GISMO_ASSERT here instead? Or throw an exception?
                    gsWarn << "Interrupting insertKnots(): Cannot insert value equal to the p+1st knot from the beginning or end into a periodic curve.\n";
                    return;
                }
            }
            /* // This way it would be nicer but it does not cover all the dubious cases.
            for( It it = begin; *it < blue1; ++it )
                *it += activeLength;
            for( It it = end-1; *it > blue2; --it )
                *it -= activeLength;*/

            std::sort( begin, end );

            // We copy some of the control points to pretend for a while that the basis is not periodic.
            //gsMatrix<T> trueCoefs = this->basis().perCoefs( this->coefs() );
            gsBoehmRefine( this->basis().knots(), this->coefs(), this->degree(), begin, end );
            //this->basis().enforceOuterKnotsPeriodic();
            //this->coefs() = trueCoefs; // Isn't this memory-leaking?
            // Periodic basis is restored using the so-called Student's Algorithm: forget as much as possible =).
            //this->coefs().conservativeResize( this->basis().size(), this->coefs().cols() );
        }
        else // non-periodic
            gsBoehmRefine( this->basis().knots(), this->coefs(), this->degree(), inBegin, inEnd);
    }

    // Look at gsGeometry class for a description
    void degreeElevate(short_t const i = 1, short_t const dir = -1);

    // Under Construction..( TO DO )
    bool contains( gsMatrix<T> const & p, T const & tol = 1e-6 )
        {
            assert( p.cols()==1 );
            gsBSplineSolver<T> slv;
            std::vector<T> roots;
            short_t dim = this->geoDim();
            gsMatrix<T> ev, tmp(1,1);
            int i(1);

            slv.allRoots( *this, roots, 0, p(0,0), 1e-6, 100);
            for ( typename std::vector<T>::const_iterator it=roots.begin();
                  it != roots.end(); ++it)
            {
                tmp(0,0)= *it;
                this->eval_into(tmp,ev);
                for (i=1; i!=dim; ++i )
                    if ( math::abs( ev(i,0) - p(i,0) ) > tol ) 
                        break;
                if (i==dim) 
                    return true;
            }
            return false;
        }

    gsMatrix<T> sample(int npoints = 50) const
    {      
        gsMatrix<T> images;
        gsMatrix<T> interval = this->parameterRange();
        const gsMatrix<T> pts = gsPointGrid(interval, npoints );
        this->eval_into( pts, images );
        return images;
    }

    void setPeriodic(bool flag = true)
    {
        this->basis().setPeriodic(flag);
        this->m_coefs = this->basis().perCoefs( this->m_coefs );
    }

    inline int numCoefs() const { return this->m_coefs.rows() - this->basis().numCrossingFunctions(); }

    inline bool isPeriodic() const { return this->basis().isPeriodic(); }


    bool isClosed(T tol = 1e-10) const 
    { 
        return ( this->basis().isPeriodic() || 
                 (m_coefs.row(0) - m_coefs.row(m_coefs.rows()-1)).squaredNorm()<tol );  
    }

    bool isOn(gsMatrix<T> const &u, T tol = 1e-3) const;

    bool isPatchCorner(gsMatrix<T> const &v, T tol = 1e-3) const;

    void findCorner(const gsMatrix<T>   & v, 
                    gsVector<index_t,1> & curr,
                    T tol = 1e-3);

    void setOriginCorner(gsMatrix<T> const &v);

    void setFurthestCorner(gsMatrix<T> const &v);

    void swapDirections(const unsigned i, const unsigned j);

protected:
    
    using Base::m_coefs;
    
    // TODO Check function
    // check function: check the coefficient number, degree, knot vector ...
}; // class gsBSpline


#ifdef GISMO_WITH_PYBIND11

  void pybind11_init_gsBSpline(pybind11::module &m);

#endif // GISMO_WITH_PYBIND11


/*
// Product of spline functions
template<class T>
gsBSpline<T> operator*(const gsBSpline<T> & lhs, const gsBSpline<T> & rhs) 
{
    // Note: quick-fix implementation. Better algorithms  exist
    gsKnotVector<T> kv = lhs.knots().knotUnion( lhs.knots() );
    kv.degreeElevate( lhs.degree() + rhs.degree() - kv.degree() );

    gsMatrix<T> pts;
    kv.greville_into(pts);
    const gsMatrix<T> ev = (lhs.eval(pts)).array() * (rhs.eval(pts)).array();
    
    // fixme: avoid temporaries here
    return *(static_cast<gsBSpline<T>*>(gsBSplineBasis<T>(kv).interpolateData(ev,pts)));
}
*/

} // namespace gismo


#ifndef GISMO_BUILD_LIB
#include GISMO_HPP_HEADER(gsBSpline.hpp)
#endif