BernsteinMulti.h

//===========================================================================
// GoTools - SINTEF Geometry Tools
//
// GoTools module: Implicitization, version 1.0
//
// 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 _BERNSTEINMULTI_H
#define _BERNSTEINMULTI_H

#include "Array.h"
#include <vector>
#include <iostream>


namespace Go {



class BernsteinPoly;


class BernsteinMulti {
public:
    BernsteinMulti() : degu_(-1), degv_(-1) { }
    BernsteinMulti(int m, int n,
                   const std::vector<double>& coefs)
        : degu_(m), degv_(n), coefs_(coefs)
    { 
        ALWAYS_ERROR_IF((degu_+1) * (degv_+1) != int(coefs_.size()),
                        "Number of coefficients does not match degrees");
    }
    template <typename ForwardIterator>
    BernsteinMulti(int m, int n,
                   ForwardIterator begin, ForwardIterator end)
        : degu_(m), degv_(n), coefs_(begin, end) { }
    explicit BernsteinMulti(double coef)
        : degu_(0), degv_(0), coefs_(std::vector<double>(1, coef)) { }

    int degreeU() const
    { return degu_; }
    int degreeV() const
    { return degv_; }

    const double& operator[] (int num) const
    { return coefs_[num]; }
    double& operator[] (int num)
    { return coefs_[num]; }
    std::vector<double>::iterator coefsBegin()
    { return coefs_.begin(); }
    std::vector<double>::const_iterator coefsBegin() const
    { return coefs_.begin(); }
    std::vector<double>::iterator coefsEnd()
    { return coefs_.end(); }
    std::vector<double>::const_iterator coefsEnd() const
    { return coefs_.end(); }

    double operator() (double u, double v) const;
    double operator() (const Vector2D& pt) const
    { return (*this)(pt[0], pt[1]); }
    double operator() (const std::vector<double>& pt) const
    { return (*this)(pt[0], pt[1]); }

    bool isZero(const double eps = 0.0) const;
    bool isStrictlyPositive(const double eps = 0.0) const;
    bool isStrictlyNegative(const double eps = 0.0) const;
    bool isNonNegative(const double eps = 0.0) const;

    double norm() const;
    void normalize();

    double mean() const;

    BernsteinMulti deriv(int der1, int der2) const;
    BernsteinMulti detHess() const;
    BernsteinMulti traceHess() const;

    double blossom(const std::vector<double>& uvec,
                   const std::vector<double>& vvec) const;
    BernsteinMulti pickDomain(double u0, double u1,
                              double v0, double v1) const;
    BernsteinPoly pickLine(Array<double,2> a, 
                           Array<double,2> b) const;

    void degreeElevate(int du, int dv);

    BernsteinMulti& operator*= (const BernsteinMulti& multi);
    BernsteinMulti& operator*= (double c);
    
    BernsteinMulti& operator+= (const BernsteinMulti& multi);
    BernsteinMulti& operator+= (double c);
    
    BernsteinMulti& operator-= (const BernsteinMulti& multi);
    BernsteinMulti& operator-= (double c);
    
    BernsteinMulti& operator/= (double c);

    BernsteinPoly bindU(double u) const;
    BernsteinPoly bindV(double v) const;

    void read(std::istream& is);
    void write(std::ostream& os) const;

private:
    // Degrees in the u'th and v'th direction.
    int degu_, degv_;
    // Coefficients are stored in a vector with (degu_+1)(degv_+1)
    // entries. The "u-index" runs fastest; i.e. for control points
    // c_ij, the sequence is
    // c00, c10, c20, ..., c01, c11, c21, ..., cmn
    std::vector<double> coefs_;
};


inline BernsteinMulti operator* (const BernsteinMulti& m1,
                                 const BernsteinMulti& m2)
{
    BernsteinMulti tmp = m1;
    tmp *= m2;
    return tmp;
}


inline BernsteinMulti operator* (const BernsteinMulti& m1,
                                 double c)
{
    BernsteinMulti tmp = m1;
    tmp *= c;
    return tmp;
}


inline BernsteinMulti operator* (double c,
                                 const BernsteinMulti& m1)
{
    BernsteinMulti tmp = m1;
    tmp *= c;
    return tmp;
}


inline BernsteinMulti operator+ (const BernsteinMulti& m1,
                                 const BernsteinMulti& m2)
{
    BernsteinMulti tmp = m1;
    tmp += m2;
    return tmp;
}


inline BernsteinMulti operator+ (const BernsteinMulti& m,
                                 double c)
{
    BernsteinMulti tmp = m;
    tmp += c;
    return tmp;
}


inline BernsteinMulti operator+ (double c,
                                 const BernsteinMulti& m)
{
    BernsteinMulti tmp = m;
    tmp += c;
    return tmp;
}


inline BernsteinMulti operator- (const BernsteinMulti& m1,
                                 const BernsteinMulti& m2)
{
    BernsteinMulti tmp = m1;
    tmp -= m2;
    return tmp;
}


inline BernsteinMulti operator- (const BernsteinMulti& m,
                                 double c)
{
    BernsteinMulti tmp = m;
    tmp -= c;
    return tmp;
}


inline BernsteinMulti operator- (double c,
                                 const BernsteinMulti& m)
{
    BernsteinMulti tmp(c);
    tmp -= m;
    return tmp;
}


inline BernsteinMulti operator/ (const BernsteinMulti& m,
                                 double c)
{
    BernsteinMulti tmp = m;
    tmp /= c;
    return tmp;
}


} // namespace Go


namespace std {



inline std::istream& operator >> (std::istream& is,
                                  Go::BernsteinMulti& m)
{
    m.read(is);
    return is;
}


inline std::ostream& operator << (std::ostream& os,
                                  const Go::BernsteinMulti& m)
{
    m.write(os);
    return os;
}


} // namespace std


#endif // _BERNSTEINMULTI_H

---