BernsteinTetrahedralPoly.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 _BERNSTEINTETRAHEDRALPOLY_H
#define _BERNSTEINTETRAHEDRALPOLY_H


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


namespace Go {



class BernsteinPoly;


class BernsteinTetrahedralPoly {
public:
    BernsteinTetrahedralPoly() : deg_(-1) { }
    BernsteinTetrahedralPoly(int deg, const std::vector<double>& coefs)
        : deg_(deg), coefs_(coefs) { }
    template <typename ForwardIterator>
    BernsteinTetrahedralPoly(int deg,
                             ForwardIterator begin, ForwardIterator end)
        : deg_(deg), coefs_(begin, end) { }
    explicit BernsteinTetrahedralPoly(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, 4>& u) const
    {
        // Generalized 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) {
                    int k = (i+1) * (i+2) / 2;
                    for (int j = 0; j <= i; ++j) {
                        for (int l = 0; l <= j; ++l) {
                            m++;
                            tmp[m] = u[0] * tmp[m]
                                + u[1] * tmp[m + k]
                                + u[2] * tmp[m + 1 + j + k]
                                + u[3] * tmp[m + 2 + j + k];
                        }
                    }
                }
        }
        
        return tmp[0];
    }

    double norm() const;
    void normalize();

    void deriv(int der, const Vector4D& d,
               BernsteinTetrahedralPoly& btp) const;

    template <typename T>
    T blossom(const std::vector<Array<T, 4> >& 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];
            const T& x = uvec[r-1][3];
            int m = -1;
            for (int i = 0; i <= deg_-r; ++i) {
                int k = (i+1) * (i+2) / 2;
                for (int j = 0; j <= i; ++j) {
                    for (int l = 0; l <= j; ++l) {
                        m++;
                        tmp[m] = u * tmp[m]
                            + v * tmp[m + k]
                            + w * tmp[m + 1 + j + k]
                            + x * tmp[m + 2 + j + k];
                    }
                }
            }
        }

        return tmp[0];
    }

    BernsteinPoly pickLine(const Array<double, 4>& a,
                           const Array<double, 4>& b) const;

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

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

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

    BernsteinTetrahedralPoly& operator/= (double c);

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

private:
    int deg_;
    std::vector<double> coefs_;
};


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


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


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


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


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


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


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


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


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


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


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


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


} // namespace Go


#endif // _BERNSTEINTETRAHEDRALPOLY_H

---