SplineSurface.h

//===========================================================================
// GoTools - SINTEF Geometry Tools version 1.1
//
// GoTools module: CORE
//
// Copyright (C) 2000-2007 SINTEF ICT, Applied Mathematics, Norway.
//
// This program is free software; you can redistribute it and/or          
// modify it under the terms of the GNU General Public License            
// as published by the Free Software Foundation version 2 of the License. 
//
// This program is distributed in the hope that it will be useful,        
// but WITHOUT ANY WARRANTY; without even the implied warranty of         
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the          
// GNU General Public License for more details.                           
//
// You should have received a copy of the GNU General Public License      
// along with this program; if not, write to the Free Software            
// Foundation, Inc.,                                                      
// 59 Temple Place - Suite 330,                                           
// Boston, MA  02111-1307, USA.                                           
//
// Contact information: E-mail: tor.dokken@sintef.no                      
// SINTEF ICT, Department of Applied Mathematics,                         
// P.O. Box 124 Blindern,                                                 
// 0314 Oslo, Norway.                                                     
//
// Other licenses are also available for this software, notably licenses
// for:
// - Building commercial software.                                        
// - Building software whose source code you wish to keep private.        
//===========================================================================
#ifndef _GOSPLINESURFACE_H
#define _GOSPLINESURFACE_H


#include "ParamSurface.h"
#include "BsplineBasis.h"
#include "RectDomain.h"
#include "ScratchVect.h"


namespace Go
{


class Interpolator;
class SplineCurve;
class DirectionCone;

class SplineSurface : public ParamSurface
{
public:
    SplineSurface()
        : dim_(-1), rational_(false)
    {
    }

    template <typename RandomIterator1,
              typename RandomIterator2,
              typename RandomIterator3>
    SplineSurface(int number1,
                    int number2,
                    int order1,
                    int order2,
                    RandomIterator1 knot1start,
                    RandomIterator2 knot2start,
                    RandomIterator3 coefsstart,
                    int dim,
                    bool rational = false)
        : dim_(dim), rational_(rational),
          basis_u_(number1, order1, knot1start),
          basis_v_(number2, order2, knot2start)
    {
        if (rational) {
            int n = (dim+1)*number1*number2;
            rcoefs_.resize(n);
            std::copy(coefsstart, coefsstart + n, rcoefs_.begin());
            coefs_.resize(dim*number1*number2);
            updateCoefsFromRcoefs();
        } else {
            int n = dim*number1*number2;
            coefs_.resize(n);
            std::copy(coefsstart, coefsstart + n, coefs_.begin());
        }
    }

    template <typename RandomIterator>
    SplineSurface(const BsplineBasis& basis_u,
                    const BsplineBasis& basis_v,
                    RandomIterator coefsstart,
                    int dim,
                    bool rational = false)
        : dim_(dim), rational_(rational),
          basis_u_(basis_u),
          basis_v_(basis_v)
    {
        int number1 = basis_u.numCoefs();
        int number2 = basis_v.numCoefs();
        if (rational) {
            int n = (dim+1)*number1*number2;
            rcoefs_.resize(n);
            std::copy(coefsstart, coefsstart + n, rcoefs_.begin());
            coefs_.resize(dim*number1*number2);
            updateCoefsFromRcoefs();
        } else {
            int n = dim*number1*number2;
            coefs_.resize(n);
            std::copy(coefsstart, coefsstart + n, coefs_.begin());
        }
    }

    virtual ~SplineSurface();

    // inherited from Streamable
    virtual void read (std::istream& is);

    // inherited from Streamable
    virtual void write (std::ostream& os) const;

    // inherited from GeomObject
    virtual BoundingBox boundingBox() const;

    // inherited from GeomObject
    virtual int dimension() const;
    
    void swap(SplineSurface& other);

    // inherited from GeomObject
    virtual ClassType instanceType() const;

    // inherited from GeomObject
    static ClassType classType()
    { return Class_SplineSurface; }
// #if ((_MSC_VER > 0) && (_MSC_VER < 1300))
//     virtual GeomObject* clone() const
//     { return new SplineSurface(*this); }
// #else
    virtual SplineSurface* clone() const
    { return new SplineSurface(*this); }
// #endif

    // inherited from ParamSurface
// #if ((_MSC_VER > 0) && (_MSC_VER < 1300))
//     virtual const Domain& parameterDomain() const;
// #else
    virtual const RectDomain& parameterDomain() const;
// #endif

    // inherited from ParamSurface
    virtual RectDomain containingDomain() const;

    // inherited from ParamSurface
    virtual CurveLoop outerBoundaryLoop(double degenerate_epsilon
                                          = DEFAULT_SPACE_EPSILON) const;
    // inherited from ParamSurface
    virtual std::vector<CurveLoop> allBoundaryLoops(double degenerate_epsilon
                                                      = DEFAULT_SPACE_EPSILON) const;

    // inherited from ParamSurface
    virtual void point(Point& pt, double upar, double vpar) const;

    // Output: Partial derivatives up to order derivs (pts[0]=S(u,v),
    // pts[1]=dS/du=S_u, pts[2]=S_v, pts[3]=S_uu, pts[4]=S_uv, pts[5]=S_vv, ...)
    // inherited from ParamSurface
    virtual void point(std::vector<Point>& pts, 
                       double upar, double vpar,
                       int derivs,
                       bool u_from_right = true,
                       bool v_from_right = true,
                       double resolution = 1.0e-12) const;

    virtual double startparam_u() const;

    virtual double startparam_v() const;

    virtual double endparam_u() const;

    virtual double endparam_v() const;

    // inherited from ParamSurface
    virtual void normal(Point& n, double upar, double vpar) const;

    enum NormalConeMethod { 
        SederbergMeyers = 0,
        SMCornersFirst = 1
    };
    DirectionCone normalCone(NormalConeMethod method) const;

    virtual DirectionCone normalCone() const;

    // inherited from ParamSurface
    virtual DirectionCone tangentCone(bool pardir_is_u) const;

    // Creates a composite box enclosing the surface. The composite box
    // consists of an inner and an edge box. The inner box is
    // made from the interior control points of the surface, while the
    // edge box is made from the boundary control points.
    // Inherited from ParamSurface
    virtual CompositeBox compositeBox() const;

    SplineSurface* normal() const;

    SplineSurface* normalSurface() const;

    SplineSurface* derivSurface(int ider1, int ider2) const;

    SplineSurface* subSurface(double from_upar,
                              double from_vpar,
                              double to_upar,
                              double to_vpar,
                              double fuzzy =
                              DEFAULT_PARAMETER_EPSILON) const;

    // inherited from ParamSurface
    virtual std::vector<boost::shared_ptr<ParamSurface> >
    subSurfaces(double from_upar, double from_vpar,
                double to_upar, double to_vpar,
                double fuzzy = DEFAULT_PARAMETER_EPSILON) const;

    // inherited from ParamSurface
    virtual void closestPoint(const Point& pt,
                              double&        clo_u,
                              double&        clo_v, 
                              Point&         clo_pt,
                              double&        clo_dist,
                              double         epsilon,
                              const RectDomain* domain_of_interest = NULL,
                              double   *seed = 0) const;

    // inherited from ParamSurface
    virtual void closestBoundaryPoint(const Point& pt,
                                      double&        clo_u,
                                      double&        clo_v, 
                                      Point&         clo_pt,
                                      double&        clo_dist,
                                      double         epsilon,
                                      const RectDomain* rd = NULL,
                                      double *seed = 0) const;

    void appendSurface(ParamSurface* sf, int join_dir,
                       int continuity, double& dist, bool repar=true);

    void appendSurface(SplineSurface* sf, int join_dir, int continuity, double& dist);

    void appendSurface(ParamSurface* sf, int join_dir, bool repar=true);

    // inherited from ParamSurface
    virtual void getBoundaryInfo(Point& pt1, Point& pt2, 
                                 double epsilon, SplineCurve*& cv,
                                 SplineCurve*& crosscv, double knot_tol = 1e-05) const;

    void getBoundaryInfo(double par1, double par2,
                         int bdindex, SplineCurve*& cv,
                         SplineCurve*& crosscv, double knot_tol = 1e-05) const;

    void getBoundaryIdx(Point& pt1, Point& pt2, 
                        double epsilon, int& bdindex,
                        double& par1, double& par2, double knot_tol = 1e-05) const;

    // corresponding boundary.  The Boundaries are indexed as in \getBoundaryIdx().
    void getBoundaryIdx(Point& pt1, double epsilon, int& bdindex,
                        double knot_tol = 1e-05) const;

    // inherited from ParamSurface
    virtual void turnOrientation();

    // inherited from ParamSurface
    virtual void swapParameterDirection();

    // inherited from ParamSurface
    virtual void reverseParameterDirection(bool direction_is_u);







    int boundaryIndex(Point& param_pt1, Point& param_pt2) const;

    const BsplineBasis& basis_u() const
    { return basis_u_; }

    const BsplineBasis& basis_v() const
    { return basis_v_; }

    const BsplineBasis& basis(int i) const
    { return (i==0) ? basis_u_ : basis_v_; }

    int numCoefs_u() const
    { return basis_u_.numCoefs(); }

    int numCoefs_v() const
    { return basis_v_.numCoefs(); }

    int order_u() const
    { return basis_u_.order(); }

    int order_v() const
    { return basis_v_.order(); }

    bool rational() const
    { return rational_; }
    
    std::vector<double>::iterator coefs_begin()
    { return coefs_.begin(); }

    std::vector<double>::iterator coefs_end()
    { return coefs_.end(); }

    std::vector<double>::const_iterator coefs_begin() const
    { return coefs_.begin(); }

    std::vector<double>::const_iterator coefs_end() const
    { return coefs_.end(); }

    std::vector<double>::iterator rcoefs_begin()
    { return rcoefs_.begin(); }

    std::vector<double>::iterator rcoefs_end()
    { return rcoefs_.end(); }

    std::vector<double>::const_iterator rcoefs_begin() const
    { return rcoefs_.begin(); }

    std::vector<double>::const_iterator rcoefs_end() const
    { return rcoefs_.end(); }

    // inherited from ParamSurface
    virtual bool isDegenerate(bool& b, bool& r,
                      bool& t, bool& l, double tolerance) const;

    void setParameterDomain(double u1, double u2, double v1, double v2);

    void insertKnot_u(double apar);
    
    void insertKnot_u(const std::vector<double>& new_knots);

    void insertKnot_v(double apar);
    
    void insertKnot_v(const std::vector<double>& new_knots);

    void removeKnot_u(double tpar);

    void removeKnot_v(double tpar);

    void makeBernsteinKnotsU();

    void makeBernsteinKnotsV();

    int numberOfPatches_u() const;

    int numberOfPatches_v() const;

    void raiseOrder(int raise_u, int raise_v);

    SplineCurve* constParamCurve(double parameter,
                                 bool pardir_is_u) const;

    void constParamCurve(double parameter, 
                         bool pardir_is_u, 
                         SplineCurve*& cv, 
                         SplineCurve*& crosscv) const;

    // inherited from ParamSurface
    virtual std::vector< boost::shared_ptr<ParamCurve> >
    constParamCurves(double parameter, bool pardir_is_u) const;



  // @@sbr What about the orientation? Seems it follows that of the sf.
    SplineCurve* edgeCurve(int ccw_edge_number) const;

    void interpolate(Interpolator& interpolator1,
                     Interpolator& interpolator2,
                     int num_points1,
                     int num_points2,
                     int dim,
                     const double* param1_start,
                     const double* param2_start,
                     const double* data_start);

    void gridEvaluator(int num_u, int num_v,
                       std::vector<double>& points,
                       std::vector<double>& normals,
                       std::vector<double>& param_u,
                       std::vector<double>& param_v) const;

    void gridEvaluator(int num_u, int num_v,
                       std::vector<double>& points,
                       std::vector<double>& param_u,
                       std::vector<double>& param_v) const;

    // inherited from ParamSurface
    virtual double nextSegmentVal(int dir, double par, bool forward, double tol) const;


private:

    struct degenerate_info
    {
        bool is_set_;
        bool b_, r_, t_, l_;
        
        degenerate_info()
        { is_set_ = false; }
    };

    // Canonical data
    int dim_;
    bool rational_;
    BsplineBasis basis_u_;
    BsplineBasis basis_v_;
    std::vector<double> coefs_;   
    std::vector<double> rcoefs_;  

    // Generated data
    mutable RectDomain domain_;
    mutable CurveLoop spatial_boundary_;
    mutable degenerate_info degen_;


    // Helper functions
    void updateCoefsFromRcoefs();
    std::vector<double>& activeCoefs() { return rational_ ? rcoefs_ : coefs_; }
    bool normal_not_failsafe(Point& n, double upar, double vpar) const;
    bool search_for_normal(bool interval_in_u,
                           double fixed_parameter,
                           double interval_start, // normal is not defined here
                           double interval_end,
                           Point& normal) const;

    // Members new to this class
    void pointsGrid(int m1, int m2, int derivs,
                    const double* basisvals1,
                    const double* basisvals2,
                    const int* knotint1,
                    const int* knotint2,
                    double* result,
                    double* normals = 0) const;

    // Rewritten pointsGrid, to avoid reformatting results.
    // Does not return the derivatives, works only for a 3D non-rational
    // spline.
    void pointsGridNoDerivs(int m1, int m2,
                            const double* basisvals1,
                            const double* basisvals2,
                            const int* knotint1,
                            const int* knotint2,
                            double* result,
                            double* normals = 0) const;

    // Actually computes the closest point. Only difference is the explicit
    // robust_seedfind-parameter, which is always true in the virtual
    // closestPoint() function, when that function calls closestPointImpl().
    void closestPointImpl(const Point& pt,
                          double&        clo_u,
                          double&        clo_v, 
                          Point&         clo_pt,
                          double&        clo_dist,
                          double         epsilon,
                          const RectDomain* domain_of_interest = NULL,
                          bool robust_seedfind = true,
                          double   *seed = 0) const;

    void s1773(const double ppoint[],double aepsge, double estart[],double eend[],double enext[],
               double gpos[],int *jstat) const;

    void s1773_s9corr(double gd[],double acoef1,double acoef2,
                      double astart1,double aend1,double astart2,double aend2) const;

    void s1773_s9dir(double *cdist,double *cdiff1,double *cdiff2,
                     double PS[],const double *eval1,std::vector<Point> eval2,
                     double aepsge, int idim,int *jstat) const;


};


} // namespace Go




#endif // _GOSPLINESURFACE_H

---