Fresco/Berlin/include/Berlin/nurbs Vertex.hh,NONE,1.1 array.hh,NONE,1.1 domain.hh,NONE,1.1 nurbs.hh,NONE,1.1 point.hh,NONE,1.1

Update of /cvs/fresco/Fresco/Berlin/include/Berlin/nurbs
In directory purcel:/tmp/cvs-serv7158/nurbs

Added Files:
	Vertex.hh array.hh domain.hh nurbs.hh point.hh 
Log Message:
some nurbs code

--- NEW FILE: Vertex.hh ---
/*$Id: Vertex.hh,v 1.1 2003/04/25 16:46:49 stefan Exp $
 *
 * This source file is a part of the Fresco Project.
 * Copyright (C) 2003 Stefan Seefeld <[email protected]> 
 * http://www.fresco.org
 *
 * This library is free software; you can redistribute it and/or
 * modify it under the terms of the GNU Library General Public
 * License as published by the Free Software Foundation; either
 * version 2 of the License, or (at your option) any later version.
 *
 * This library 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
 * Library General Public License for more details.
 *
 * You should have received a copy of the GNU Library General Public
 * License along with this library; if not, write to the
 * Free Software Foundation, Inc., 675 Mass Ave, Cambridge,
 * MA 02139, USA.
 */
#ifndef _Berlin_nurbs_Vertex_hh
#define _Berlin_nurbs_Vertex_hh

namespace Berlin
{
namespace nurbs
{

inline Vertex make_vertex(double x, double y, double z)
{
  Vertex v;
  v.x = x;
  v.y = y;
  v.z = z;
  return v;
}

inline Vertex &assign(Vertex &v, double s) { v.x = v.y = v.z = s; return v;}

inline Vertex &operator+=(Vertex &p, const Vertex &q)
{
  p.x += q.x;
  p.y += q.y;
  p.z += q.z;
  return p;
}

inline Vertex &operator-=(Vertex &p, const Vertex &q)
{
  p.x -= q.x;
  p.y -= q.y;
  p.z -= q.z;
  return p;
}

inline Vertex &operator*=(Vertex &v, double s)
{
  v.x *= s;
  v.y *= s;
  v.z *= s;
  return v;
}

inline Vertex &operator/=(Vertex &v, double s)
{
  v.x /= s;
  v.y /= s;
  v.z /= s;
  return v;
}

inline Vertex operator+(const Vertex &p, const Vertex &q)
{
  Vertex result(p);
  return result += q;
}

inline Vertex operator-(const Vertex &p, const Vertex &q)
{
  Vertex result(p);
  return result -= q;
}

inline Vertex operator*(const Vertex &p, double scalar)
{
  Vertex result(p);
  return result *= scalar;
}

inline Vertex operator/(const Vertex &p, double scalar)
{
  Vertex result(p);
  return result /= scalar;
}

inline double scalar(const Vertex &p, const Vertex &q)
{
  return p.x * q.x + p.y * q.y + p.z * q.z;
}

inline double norm(const Vertex &p)
{
  return std::sqrt(scalar(p, p));
}

}
}

#endif

--- NEW FILE: array.hh ---
/*$Id: array.hh,v 1.1 2003/04/25 16:46:49 stefan Exp $
 *
 * This source file is a part of the Fresco Project.
 * Copyright (C) 2003 Stefan Seefeld <[email protected]> 
 * http://www.fresco.org
 *
 * This library is free software; you can redistribute it and/or
 * modify it under the terms of the GNU Library General Public
 * License as published by the Free Software Foundation; either
 * version 2 of the License, or (at your option) any later version.
 *
 * This library 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
 * Library General Public License for more details.
 *
 * You should have received a copy of the GNU Library General Public
 * License along with this library; if not, write to the
 * Free Software Foundation, Inc., 675 Mass Ave, Cambridge,
 * MA 02139, USA.
 */
#ifndef _Berlin_nurbs_array_hh
#define _Berlin_nurbs_array_hh

namespace Berlin
{
namespace nurbs
{

template<class T, std::size_t N>
class array 
{
public:
  typedef T              value_type;
  typedef T             *iterator;
  typedef const T       *const_iterator;
  typedef std::reverse_iterator<iterator> reverse_iterator;
  typedef std::reverse_iterator<const_iterator> const_reverse_iterator;
  typedef T             &reference;
  typedef const T       &const_reference;
  typedef std::size_t    size_type;
  typedef std::ptrdiff_t difference_type;
    
  array() {}
  array(const array<T, N> &other) { std::copy(other.begin(), other.end(), begin());}
  array(T value) { std::fill_n(begin(), N, value);}
  template <typename T2>
  array<T,N> &operator= (const array<T2,N> &rhs) 
  {
    std::copy(rhs.begin(),rhs.end(), begin());
    return *this;
  }
  void assign(const T &value) { std::fill_n(begin(), N, value);}

  iterator begin() { return my_data;}
  const_iterator begin() const { return my_data;}
  iterator end() { return my_data+N;}
  const_iterator end() const { return my_data+N;}
  
  reverse_iterator rbegin() { return reverse_iterator(end());}
  const_reverse_iterator rbegin() const { return const_reverse_iterator(end());}
  reverse_iterator rend() { return reverse_iterator(begin());}
  const_reverse_iterator rend() const { return const_reverse_iterator(begin());}

  reference operator[](size_type i) { return my_data[i];}
  const_reference operator[](size_type i) const { return my_data[i];}

  reference front() { return my_data[0];}
  const_reference front() const { return my_data[0];}
  reference back() { return my_data[N-1];}
  const_reference back() const { return my_data[N-1];}
  
  void swap (array<T,N> &y) { std::swap_ranges(begin(),end(),y.begin());}
  const T *data() const { return my_data;}
private:
  T my_data[N];
};

template<class T, std::size_t N>
bool operator== (const array<T,N> &x, const array<T,N> &y) 
{
  return std::equal(x.begin(), x.end(), y.begin());
}

template<class T, std::size_t N>
bool operator< (const array<T,N> &x, const array<T,N> &y) 
{
  return std::lexicographical_compare(x.begin(),x.end(),y.begin(),y.end());
}

template<class T, std::size_t N>
bool operator!= (const array<T,N> &x, const array<T,N> &y)
{
  return !(x==y);
}

template<class T, std::size_t N>
bool operator> (const array<T,N> &x, const array<T,N> &y)
{
  return y<x;
}

template<class T, std::size_t N>
bool operator<= (const array<T,N> &x, const array<T,N> &y)
{
  return !(y<x);
}

template<class T, std::size_t N>
bool operator>= (const array<T,N> &x, const array<T,N> &y)
{
  return !(x<y);
}

template<class T, std::size_t N>
inline void swap (array<T,N> &x, array<T,N> &y)
{
  x.swap(y);
}

}
}

#endif

--- NEW FILE: domain.hh ---
/*$Id: domain.hh,v 1.1 2003/04/25 16:46:49 stefan Exp $
 *
 * This source file is a part of the Fresco Project.
 * Copyright (C) 2003 Stefan Seefeld <[email protected]> 
 * http://www.fresco.org
 *
 * This library is free software; you can redistribute it and/or
 * modify it under the terms of the GNU Library General Public
 * License as published by the Free Software Foundation; either
 * version 2 of the License, or (at your option) any later version.
 *
 * This library 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
 * Library General Public License for more details.
 *
 * You should have received a copy of the GNU Library General Public
 * License along with this library; if not, write to the
 * Free Software Foundation, Inc., 675 Mass Ave, Cambridge,
 * MA 02139, USA.
 */
#ifndef _Berlin_nurbs_domain_hh
#define _Berlin_nurbs_domain_hh

#include <ostream>
#include <algorithm>
#include <cstdarg>

namespace Berlin
{
namespace nurbs
{

template <size_t D>
struct index
{
  index() {}
  index(const size_t *i) { std::copy(i, i + D, values);}
  index(index<D-1> rhs, size_t lhs)
  {
    std::copy(rhs.values, rhs.values + D - 1, values);
    values[D-1] = lhs;
  }
  index(const index<D> &idx) { std::copy(idx.values, idx.values + D, values);}
  index<D+1> operator[](size_t i) const { return index<D+1>(*this, i);}
  size_t values[D+1];
};

template <typename T, size_t D> 
class domain
{
public:
  typedef size_t vector_index;
  static const nurbs::index<0> index;
  static const size_t dimensions = D;

  domain(const nurbs::index<D> &);
  domain(const size_t *);
  domain(const domain<T, D> &);
  domain<T,D> &operator=(const domain<T, D> &);
  ~domain();

  size_t size(size_t i) const { return my_sizes.values[i];}
  const size_t *sizes() const { return my_sizes.values;}
  const T *data() const { return my_data;}
  size_t length() const { return my_length;}

  const T &operator() (nurbs::index<D>) const;
  T &operator() (nurbs::index<D>);

  const T &operator[](vector_index i) const { return my_data[i];}
  T &operator[](vector_index i) { return my_data[i];}
private:
  nurbs::index<D> my_sizes;
  size_t          my_length;
  T              *my_data;
};

template <typename T, size_t D>
const nurbs::index<0> domain<T, D>::index;

template <typename T, size_t D>
domain<T, D>::domain(const nurbs::index<D> &idx)
  : my_sizes(idx)
{
  my_length = 1;
  for(size_t i = 0; i < D; ++i) my_length *= idx.values[i];
  my_data = new T[my_length];
}

template <typename T, size_t D>
domain<T, D>::domain(const size_t *idx)
  : my_sizes(idx)
{
  my_length = 1;
  for(size_t i = 0; i < D; ++i) my_length *= idx[i];
  my_data = new T[my_length];
}

template <typename T, size_t D>
domain<T, D>::domain(const domain<T, D> &m)
  : my_sizes(m.my_sizes),
    my_length(m.my_length),
    my_data(new T[my_length])
{
  std::copy(m.my_data, m.my_data + my_length, my_data);
}

template <typename T, size_t D>
domain<T, D>::~domain()
{
  delete [] my_data;
}

template <typename T, unsigned int D>
domain<T, D> &domain<T, D>::operator=(const domain<T, D> &m)
{
  my_length = m.my_length;
  my_data = m.my_data;
  my_sizes = m.my_sizes;
  return *this;
}

template<typename T, size_t D>
const T &domain<T,D>::operator()(nurbs::index<D> idx) const
{
  size_t sum = 0;
  for(size_t j = 0; j < D; ++j)
  {
    size_t prod = 1;
    for(size_t k = 0; k < j; ++k) prod *= my_sizes.values[k];
    sum += idx.values[j] * prod;  
    prod = 1;
  }
  return my_data[sum];
}

template<typename T, size_t D>
T &domain<T,D>::operator()(nurbs::index<D> idx)
{
  size_t sum = 0;
  for(size_t j = 0; j < D; ++j)
  {
    size_t prod = 1;
    for(size_t k = 0; k < j; ++k) prod *= my_sizes.values[k];
    sum += idx.values[j] * prod;  
    prod = 1;
  }
  return my_data[sum];
}

}
}

#endif

--- NEW FILE: nurbs.hh ---
/*$Id: nurbs.hh,v 1.1 2003/04/25 16:46:49 stefan Exp $
 *
 * This source file is a part of the Fresco Project.
 * Copyright (C) 2002 Massimo Ricci
 * Copyright (C) 2003 Stefan Seefeld <[email protected]> 
 * http://www.fresco.org
 *
 * This library is free software; you can redistribute it and/or
 * modify it under the terms of the GNU Library General Public
 * License as published by the Free Software Foundation; either
 * version 2 of the License, or (at your option) any later version.
 *
 * This library 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
 * Library General Public License for more details.
 *
 * You should have received a copy of the GNU Library General Public
 * License along with this library; if not, write to the
 * Free Software Foundation, Inc., 675 Mass Ave, Cambridge,
 * MA 02139, USA.
 */
#ifndef _Berlin_nurbs_basic_nurbs_hh
#define _Berlin_nurbs_basic_nurbs_hh

#include <Berlin/nurbs/domain.hh>
#include <Berlin/nurbs/array.hh>
#include <vector>
#include <string>

namespace Berlin
{
namespace nurbs
{

template <typename T>
struct Weighted
{
  T numerator;
  double denominator;
};

template <typename K>
double basis_function(size_t i, size_t d, double p, const K &k)
{
  if(d == 0)
  {
    if(k[i] <= p && p <= k[i + 1]) return 1.;
    else return 0.;
  }
  else
  {
    if(k[i] <= p && p <= k[i + d + 1])
    {
      double e1 = k[i + d] - k[i];
      double e2 = k[i + d + 1] - k[i + 1];

      if (e1 == 0. && e2 == 0.) return 0.;
	  
      if (e1 == 0. && e2 != 0.)
	return (k[i + d + 1] - p) / e2 * basis_function(i + 1, d - 1, p, k);
	  
      if (e1 != 0. && e2 == 0.)
	return (p - k[i]) / e1 * basis_function(i, d - 1, p, k);

      if (e1 != 0. && e2 != 0.)
	return (k[i + d + 1] - p) / e2 * basis_function(i + 1, d - 1, p, k) +
	  (p - k[i]) / e1 * basis_function(i, d - 1, p, k);
    }
    else return 0.;
  }
  return 0.;
}

template <typename K>
double basis_function_derivative(size_t i, size_t d, double p, const K &k)
{
  if (d == 0) return 0.;
  else
  {
    double d1 = k[i + d] - k[i];
    double d2 = k[i + d + 1] - k[i + 1];
      
    double e1, e2, e3, e4;
      
    e1 = e2 = e3 = e4 = 0.0;
      
    if(d1 != 0.)
    {
      e1 = (1. / d1) * basis_function(i, d - 1, p, k);
      e2 = (p - k[i]) / d1 * basis_function_derivative(i, d - 1, p, k);
    }
      
    if(d2 != 0.)
    {
      e3 = (-1. / d2) * basis_function(i + 1, d - 1, p, k);
      e4 = (k[i + d + 1] - p) / d2 * basis_function_derivative(i + 1, d - 1, p, k);
    }
    return e1 + e2 + e3 + e4;
  }
  return 0.;
}

template <typename T, size_t P>
Weighted<T> eval_at(const domain<T, P> &ctrls,
                    const domain<double, P> weights,
                    const array<size_t, P> &degrees,
                    const array<std::vector<double>, P> &knots,
                    const array<size_t, P> &steps,
                    array<double, P> params,
                    array<size_t, P> idxs,
                    size_t r)
{
  Weighted<T> sum;
  sum.numerator = 0.;
  sum.denominator = 0.;

  if(r > 1)
  {
    for(idxs[r - 1] = 0; idxs[r - 1] < ctrls.size(r - 1); ++idxs[r - 1])
    {
      double basis = basis_function(idxs[r - 1], degrees[r - 1],
                                    params[r - 1], knots[r - 1]);
      if(basis != 0.)
      {
        Weighted<T> point = eval_at(ctrls, weights, degrees, knots, steps,
                                    params, idxs, r - 1);
        sum.numerator += point.numerator * basis;
        sum.denominator += point.denominator * basis;
      }
    }
  }
  else
  {
    for(idxs[r - 1] = 0; idxs[r - 1] < ctrls.size(r - 1); ++idxs[r - 1])
    {
      double basis = basis_function(idxs[r - 1], degrees[r - 1],
				    params[r - 1], knots[r - 1]);
      if(basis != 0.)
      {
        sum.numerator += ctrls(idxs.data()) * basis * weights(idxs.data());
        sum.denominator += basis * weights(idxs.data());
      }
    }
  }
  return sum;
}

template <typename T, size_t P>
T eval_at(const domain<T, P> &ctrls,
          const domain<double, P> &weights,
          const array<size_t, P> &degrees,
          const array<std::vector<double>, P> &knots,
          const array<size_t, P> &steps,
          const array<double, P> &params)
{
  array<size_t, P> idx(0);
  Weighted<T> point = eval_at(ctrls, weights, degrees, knots, steps,
                              params, idx, P);
  return point.numerator / point.denominator;
}

template <typename T, size_t P>
void eval(const domain<T, P> &ctrls, 
          const domain<double, P> &weights,
          const array<size_t, P> &degrees,
          const array<std::vector<double>, P> &knots,
          const array<size_t, P> &steps,
          array<double, P> params,
          size_t r,
          std::vector<T> &points)
{
  size_t step = steps[r - 1];
  size_t length = ctrls.size(r - 1) + 1;
  size_t plength = length - degrees[r - 1] +
    (length - degrees[r - 1] - 1) * (step - 1);

  std::vector<double> pvalues(plength);

  size_t index = degrees[r - 1];

  if(step == 1)
  {
    pvalues = knots[r - 1];
    plength = length;
  }
  else
  {
    size_t h = 0;
    for(size_t i = index; i < length - 1; ++i)
    {
      pvalues[h] = knots[r - 1][i];
      double d = (knots[r - 1][i + 1] - knots[r - 1][i]) / step;
      for(size_t j = 0; j < step - 1; ++j)
      {
        ++h;
        pvalues[h] = pvalues[h - 1] + d;
      }
      ++h;
    }
    pvalues[h] = knots[r - 1][length - 1];
    index = 0;
  }
   
  for(size_t i = index; i < plength; ++i)
  {
    params[r - 1] = pvalues[i];
    if(r == 1)
      points.push_back(eval_at(ctrls, weights, degrees, knots, steps,
                               params));
    else eval(ctrls, weights, degrees, knots, steps,
              params, r - 1, points);
  }
}

template <typename T, size_t P>
domain<T, P> *eval(const domain<T, P> &ctrls,
		   const domain<double, P> &weights,
		   const array<size_t, P> &degrees,
		   const array<std::vector<double>, P> &knots,
		   const array<size_t, P> &steps)
{
  std::vector<T> points;
  array<size_t, P> length;
  array<double, P> params;

  for(size_t i = 0; i < P; ++i)
  {
    params[i] = knots[i][0];
    length[i] = ctrls.size(i) - degrees[i] + 1 +
      (ctrls.size(i) - degrees[i]) * (steps[i] - 1);
  }
  
  eval(ctrls, weights, degrees, knots, steps, params, P, points);

  domain<T, P> *retn = new domain<T, P>(length.data());  
  for(size_t i = 0; i < points.size(); ++i) (*retn)[i] = points[i];
  return retn;
}

template <typename T, size_t P>
array<T, P + 1>
eval_with_derivations_at(const domain<T, P> &ctrls,
                         const domain<double, P> &weights,
                         const array<size_t, P> &degrees,
                         const array<std::vector<double>, P> &knots,
                         const array<size_t, P> &steps,
                         const array<double, P> &params)
{
  array<size_t, P> idx(0);
  array<Weighted<T>, P + 1> points = eval_with_derivations_at(ctrls, weights,
                                                              degrees, knots, steps,
                                                              params, idx, P);
  T numerator = points[0].numerator;
  double denominator = points[0].denominator;

  array<T, P + 1> derivations;

  for(size_t k = 1; k < P + 1; ++k)
  {
    derivations[k] = points[k].numerator * denominator;
    derivations[k] -= numerator * points[k].denominator;
    if(denominator != 0) derivations[k] /= denominator * denominator;
    else assign(derivations[k], 0.);
  }

  if(denominator != 0) derivations[0] = numerator / denominator;
  else assign(derivations[0], 0.);
  return derivations;
}

template <typename T, size_t P>
domain<array<T, P + 1>, P> *
eval_with_derivations(const domain<T, P> &ctrls,
                      const domain<double, P> &weights,
                      const array<size_t, P> &degrees,
                      const array<std::vector<double>, P> &knots,
                      const array<size_t, P> &steps)
{
  array<size_t, P> length;
  array<double, P> params;

  for(size_t i = 0; i < P; ++i)
  {
    params[i] = knots[i][0];
    length[i] = ctrls.size(i) - degrees[i] + 1 +
      (ctrls.size(i) - degrees[i]) * (steps[i] - 1);
  }
   
  std::vector<array<T, P + 1> > points;
  eval_with_derivations(ctrls, weights, degrees, knots, steps,
                        params, P, points);

  domain<array<T, P + 1>, P> *retn = new domain<array<T, P + 1>, P>(length.data());
  for(size_t i = 0; i < points.size(); ++i) (*retn)[i] = points[i];
  return retn;
}

template <typename T, size_t P>
array<Weighted<T>, P + 1>
eval_with_derivations_at(const domain<T, P> &ctrls,
                         const domain<double, P> &weights,
                         const array<size_t, P> &degrees,
                         const array<std::vector<double>, P> &knots,
                         const array<size_t, P> &steps,
                         array<double, P> params,
                         array<size_t, P> idxs,
                         size_t r)
{
  array<Weighted<T>, P + 1> sum;
  for(size_t i = 0; i < P + 1; ++i)
  {
    assign(sum[i].numerator, 0.);
    sum[i].denominator = 0.;
  }
  if(r > 1)
  {
    array<Weighted<T>, P + 1> points;
    for(idxs[r - 1] = 0; idxs[r - 1] < ctrls.size(r - 1); ++idxs[r - 1])
    {
      double derivative = basis_function_derivative(idxs[r - 1], degrees[r - 1],
                                                    params[r - 1], knots[r - 1]);
      double basis = basis_function(idxs[r - 1], degrees[r - 1],
                                    params[r -  1], knots[r - 1]);
      if(basis != 0 || derivative != 0)
        points = eval_with_derivations_at(ctrls, weights, degrees, knots, steps,
                                          params, idxs, r - 1);
      if(basis != 0)
      {
        sum[0].numerator += points[0].numerator * basis;
        sum[0].denominator += points[0].denominator * basis;
      }  
      for(size_t k = 1; k < P + 1; ++k)
        if(k == r)
        {
          if(derivative != 0)
          {
            sum[k].numerator += points[k].numerator * derivative;
            sum[k ].denominator += points[k].denominator * derivative;
          }
        }
        else
        {
          if(basis != 0)
          {
            sum[k].numerator += points[k].numerator * basis;
            sum[k].denominator += points[k].denominator * basis;
          }
        }
    }
  }
  else
  {
    for(idxs[r - 1] = 0; idxs[r - 1] < ctrls.size(r - 1); ++idxs[r - 1])
    {
      double basis = basis_function(idxs[r - 1], degrees[r - 1],
                                    params[r - 1], knots[r - 1]);
      double derivative = basis_function_derivative(idxs[r - 1], degrees[r - 1],
                                                    params[r - 1], knots[r - 1]);
      if(basis != 0)
      {
        sum[0].numerator += ctrls(idxs.data()) * basis * weights(idxs.data());
        sum[0].denominator += basis * weights(idxs.data());
      }
      for(size_t k = 1; k < P + 1; ++k)
        if(k == r && derivative != 0)
        {
          sum[k].numerator += ctrls(idxs.data()) * derivative * weights(idxs.data());
          sum[k].denominator += weights(idxs.data()) * derivative;
        }
        else
        {
          if(basis != 0)
          {
            sum[k].numerator += ctrls(idxs.data()) * basis * weights(idxs.data());
            sum[k].denominator += weights(idxs.data()) * basis;
          }
        }
    }
  }
  return sum;
}

template <typename T, size_t P>
void eval_with_derivations(const domain<T, P> &ctrls,
                           const domain<double, P> &weights,
                           const array<size_t, P> &degrees,
                           const array<std::vector<double>, P> knots,
                           const array<size_t, P> &steps,
                           array<double, P> params,
                           size_t r,
                           std::vector<array<T, P + 1> > &points)
{
  size_t step = steps[r - 1];
  size_t length = ctrls.size(r - 1) + 1;
  size_t plength = length - degrees[r - 1] +
    (length - degrees[r - 1] - 1) * (step - 1);

  std::vector<double> pvalues(plength);
  size_t index = degrees[r - 1];

  if(step == 1)
  {
    std::copy(knots[r-1].begin(), knots[r-1].end(), pvalues.begin());
    plength = length;
  }
  else
  {
    size_t h = 0;
    for(size_t i = index; i < length - 1; ++i)
    {
      pvalues[h] = knots[r - 1][i];
      double d = (knots[r - 1][i + 1] - knots[r - 1][i]) / step;
      for(size_t j = 0; j < step - 1; ++j)
      {
        ++h;
        pvalues[h] = pvalues[h - 1] + d;
      }
      ++h;
    }
    pvalues[h] = knots[r - 1][length - 1];
    index = 0;
  }
  for(size_t i = index; i < plength; ++i)
  {
    params[r - 1] = pvalues[i];
    if(r == 1)
    {
      points.push_back(eval_with_derivations_at(ctrls, weights, degrees,
                                                knots, steps, params)); 
    }
    else eval_with_derivations(ctrls, weights, degrees, knots, steps,
                               params, r - 1, points);
  }
}

template <typename T, size_t P>
domain<T, P> *evaluate(const domain<T, P> &ctrls,
                       const domain<double, P> &weights,
                       const array<size_t, P> &degrees,
                       const array<size_t, P> &steps)
{
  array<std::vector<double>, P> knots;
  for(size_t i = 0; i < P; ++i)
  {
    knots[i].resize(ctrls.size(i) + degrees[i] + 1);
    double delta = 1. / (ctrls.size(i) + degrees[i]);
    for (size_t j = 0; j != ctrls.size(i) + degrees[i] + 1; ++j)
      knots[i][j] = j * delta;
  }
  return eval(ctrls, weights, degrees, knots, steps);
}

template <typename T, size_t P>
domain<T, P> *evaluate(const domain<T, P> &ctrls,
                       const domain<double, P> &weights,
                       const array<size_t, P> &degrees,
                       const array<std::vector<double>, P> &knots,
                       const array<size_t, P> &steps)
{
  return eval(ctrls, weights, degrees, knots, steps);
}

template <typename T, size_t P>
domain<array<T, P + 1>, P> *
evaluate_with_derivations(const domain<T, P> &ctrls,
                          const domain<double, P> &weights,
                          const array<size_t, P> &degrees,
                          const array<size_t, P> &steps)
{
  array<std::vector<double>, P> knots;
  for(size_t i = 0; i < P; ++i)
  {
    knots[i].resize(ctrls.size(i) + degrees[i] + 1);
    double delta = 1. / (ctrls.size(i) + degrees[i]);
    for (size_t j = 0; j != ctrls.size(i) + degrees[i] + 1; ++j)
      knots[i][j] = j * delta;
  }
  return eval_with_derivations(ctrls, weights, degrees, knots, steps);
}

template <typename T, size_t P>
domain<array<T, P + 1>, P> *
evaluate_with_derivations(const domain<T, P> &ctrls,
                          const domain<double, P> &weights,
                          const array<size_t, P> &degrees,
                          const array<std::vector<double>, P> &knots,
                          const array<size_t, P> &steps)
{
  return eval_with_derivations(ctrls, weights, degrees, knots, steps);
}

}
}

#endif





--- NEW FILE: point.hh ---
/*$Id: point.hh,v 1.1 2003/04/25 16:46:49 stefan Exp $
 *
 * This source file is a part of the Fresco Project.
 * Copyright (C) 2003 Stefan Seefeld <[email protected]> 
 * http://www.fresco.org
 *
 * This library is free software; you can redistribute it and/or
 * modify it under the terms of the GNU Library General Public
 * License as published by the Free Software Foundation; either
 * version 2 of the License, or (at your option) any later version.
 *
 * This library 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
 * Library General Public License for more details.
 *
 * You should have received a copy of the GNU Library General Public
 * License along with this library; if not, write to the
 * Free Software Foundation, Inc., 675 Mass Ave, Cambridge,
 * MA 02139, USA.
 */
#ifndef _Berlin_nurbs_point_hh
#define _Berlin_nurbs_point_hh

namespace Berlin
{
namespace nurbs
{

template <typename T, size_t D>
class point
{
public:
  static const size_t dimensions = D;
  point() { for(size_t i = 0; i < D; ++i) my_coords[i] = 0;}

  T &operator[](size_t i) { return my_coords[i];}
  const T &operator[](size_t i) const { return my_coords[i];}

  point<T,D> &operator+=(const point<T,D> &);
  point<T,D> &operator-=(const point<T,D> &);
  point<T,D> &operator*=(T);
  point<T,D> &operator/=(T);
private:
  T my_coords[D];
};

template <typename T, size_t D>
point<T, D> &assign(point<T, D> &p, T s) { p[0] = p[1] = p[2] = s; return p;}

template <typename T, size_t D>
point<T, D> &point<T, D>::operator+=(const point<T,D> &p)
{
  T *tmp1 = my_coords;
  const T *tmp2 = p.my_coords;
  for(size_t i = 0; i < D; ++i) *tmp1++ += *tmp2++;
  return *this;
}

template <typename T, size_t D>
point <T, D> &point<T, D>::operator-=(const point<T,D> &p)
{
  T *tmp1 = my_coords;
  const T *tmp2 = p.my_coords;
  for(size_t i = 0; i < D; ++i) *tmp1++ -= *tmp2++;
  return *this;
}

template <typename T, size_t D>
point <T, D> &point<T, D>::operator*=(T s)
{
  T *tmp1 = my_coords;
  for(size_t i = 0; i < D; ++i) *tmp1++ *= s;
  return *this;
}

template <typename T, size_t D>
point <T, D> &point<T, D>::operator/=(T s)
{
  T *tmp1 = my_coords;
  for(size_t i = 0; i < D; ++i) *tmp1++ /= s;
  return *this;
}

template <typename T, size_t D>
point<T, D> operator+(const point<T, D> &p, const point<T, D> &q)
{
  point<T, D> result(p);
  return result += q;
}

template <typename T, size_t D>
point<T, D> operator-(const point<T, D> &p, const point<T, D> &q)
{
  point<T, D> result(p);
  return result -= q;
}

template <typename T, size_t D>
point<T, D> operator*(const point<T, D> &p, T scalar)
{
  point<T, D> result(p);
  return result *= scalar;
}

template <typename T, size_t D>
point <T, D> operator/(const point<T, D> &p, T scalar)
{
  point<T, D> result(p);
  return result /= scalar;
}

template <typename T, size_t D>
T scalar(const point<T, D> &p, const point<T, D> &q)
{
  T sum = 0.;
  for(size_t i = 0; i < D; ++i) sum += p[i] * q[i];
  return sum;
}

template <typename T, size_t D>
T norm(const point<T, D> &p)
{
  return std::sqrt(scalar(p, p));
}

}
}

#endif