BernsteinTriangularPoly.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 _BERNSTEINTRIANGULARPOLY_H
#define _BERNSTEINTRIANGULARPOLY_H


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


namespace Go {



class BernsteinTriangularPoly {
public:
    BernsteinTriangularPoly() : deg_(-1) { }
    BernsteinTriangularPoly(int deg, const std::vector<double>& coefs)
        : deg_(deg), coefs_(coefs) { }
    template <typename ForwardIterator>
    BernsteinTriangularPoly(int deg,
                            ForwardIterator begin, ForwardIterator end)
        : deg_(deg), coefs_(begin, end) { }
    explicit BernsteinTriangularPoly(double coef)
        : deg_(0), coefs_(1, coef) { }

    int degree() const
    { return deg_; }

    const double operator[] (int num) const
    { return coefs_[num]; }
    double& operator[] (int num)
    { return coefs_[num]; }

    template <typename T>
    T operator() (const Array<T, 3>& u) const
    {
        // Adapted from the de Casteljau algorithm 'tri_decas' in Farin,
        // "Curves and Surfaces for CAGD".

        ASSERT(deg_ >= 0);

        int sz = coefs_.size();
        std::vector<T> tmp(sz);
        for (int i = 0; i < sz; ++i) {
            tmp[i] = T(coefs_[i]);
        }

        for (int r = 1; r <= deg_; ++r) {
            int m = -1;
            for (int i = 0; i <= deg_-r; ++i) {
                for (int l = 0; l <= i; ++l) {
                    m++;
                    tmp[m] = u[0] * tmp[m]
                        + u[1] * tmp[m + 1 + i]
                        + u[2] * tmp[m + 2 + i];
                }
            }
        }

        return tmp[0];
    }

    double norm() const;
    void normalize();

    void deriv(int der, const Vector3D& d,
               BernsteinTriangularPoly& btp) const;

    template <typename T>
    T blossom(const std::vector<Array<T, 3> >& uvec) const
    {
        // This routine is a natural extension of operator().

        ASSERT(deg_ >= 0);

        int sz = coefs_.size();
        std::vector<T> tmp(sz);
        for (int i = 0; i < sz; ++i)
            tmp[i] = T(coefs_[i]);

        for (int r = 1; r <= deg_; ++r) {
            const T& u = uvec[r-1][0];
            const T& v = uvec[r-1][1];
            const T& w = uvec[r-1][2];
            int m = -1;
            for (int i = 0; i <= deg_-r; ++i) {
                for (int l = 0; l <= i; ++l) {
                    m++;
                    tmp[m] = u * tmp[m]
                        + v * tmp[m + 1 + i]
                        + w * tmp[m + 2 + i];
                }
            }
        }

        return tmp[0];
    }

    BernsteinTriangularPoly& operator*= (const BernsteinTriangularPoly& poly);
    BernsteinTriangularPoly& operator*= (double c);

    BernsteinTriangularPoly& operator+= (const BernsteinTriangularPoly& poly);
    BernsteinTriangularPoly& operator+= (double c);

    BernsteinTriangularPoly& operator-= (const BernsteinTriangularPoly& poly);
    BernsteinTriangularPoly& operator-= (double c);

    BernsteinTriangularPoly& operator/= (double c);

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

private:
    int deg_;
    //typedef ScratchVect<double, 21> VecType;
    typedef std::vector<double> VecType;
    VecType coefs_;
};


inline BernsteinTriangularPoly operator* (const BernsteinTriangularPoly& p1,
                                          const BernsteinTriangularPoly& p2)
{
    BernsteinTriangularPoly tmp = p1;
    tmp *= p2;
    return tmp;
}


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


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


inline BernsteinTriangularPoly operator+ (const BernsteinTriangularPoly& p1,
                                          const BernsteinTriangularPoly& p2)
{
    BernsteinTriangularPoly tmp = p1;
    tmp += p2;
    return tmp;
}


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


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


inline BernsteinTriangularPoly operator- (const BernsteinTriangularPoly& p1,
                                          const BernsteinTriangularPoly& p2)
{
    BernsteinTriangularPoly tmp = p1;
    tmp -= p2;
    return tmp;
}


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


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


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


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

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


} // namespace Go



#endif // _BERNSTEINTRIANGULARPOLY_H

---