Very Computer

I have a problem calculating the basis functions for knot vectors that
are not clamped.
I use the algorithm from the Piegl/Tiller NURBS book (ALGORITHM A2.2, page 70)
Piegl/Tiller always assume that the knot vector is clamped, i.e. knot values
are of multiplicity order (degree + 1) at the start and the end.
Unfortunately the data I have to deal with doesn't always looks like that,
e.g. I have a knot vector for a curve of degree 3 that begins like:
-7.62939453147204E-06,-7.62939453147204E-06,0.0,0.0,0.00390624999999359 ...
and even ends with:
1.,1.,1.00390624999999,1.00781249999999
i.e. multiplities look like:
2,2,1...2,1,1 instead of 4,1...1,4

I already changed the algorithm to find the span for a parameter u to return
the real span [1..m-1] instead of Piegl/Tillers version that returns only spans
between [p..m-p].

But now I have problems (I think) getting the basis functions and evaluating the
curve. My methods GetBasisFunctions (what implements Piegl/Tillers BasisFuns)
and Evaluate (what implements Piegl/Tillers CurvePoint) look like follows:

  DoubleVector
  B_Spline_Curve::GetBasisFunctions(double u) const
  {
    DoubleVector U = GetKnotVector();
    const int i = FindSpan(u);
    DoubleVector N(m_degree + 1);
    unsigned int p = m_degree;
    N[0] = 1.;
    DoubleVector left(p + 1);
    DoubleVector right(p + 1);
    for (unsigned int j = 1; j <= p; ++j) {
      left[j] = u - U[i + 1 - j];
      right[j] = U[i + j] - u;
      double saved = 0.;
      for (unsigned int r = 0; r < j; ++r) {
        double temp = N[r] / (right[r + 1] + left[j - r]);
        N[r] = saved + right[r + 1] * temp;
        saved = left[j - r] * temp;
      }
      N[j] = saved;
    }
    return N;
  }

  Cartesian_Point*
  B_Spline_Curve::Evaluate(ParameterValue value) const
  {
    Cartesian_Point* result = NULL;
    int span = FindSpan(value);
    DoubleVector knotvector = GetKnotVector();
    int ksize = knotvector.size();
    assert(span > 0 && span < ksize);
    if ((span > 0) && (span < (ksize))) {
      result = new Cartesian_Point(0., 0., 0.);
      DoubleVector basisfunction = GetBasisFunctions(value);
      int ctrlpoint_index = 0;
      int ctrlpoint_count = m_control_points_list.size();
      for (unsigned int i = 0; i <= m_degree; ++i) {
        ctrlpoint_index = span - m_degree + i;
        if (ctrlpoint_index >= 0 && ctrlpoint_index < ctrlpoint_count - 1) {
          Cartesian_Point* pnt = m_control_points_list[ctrlpoint_index];
          Cartesian_Point* epnt = basisfunction[i] * *pnt;
          Cartesian_Point* old = result;
          result = *old + *epnt;
        }
      }
    }
    return result;
  }

Can somebody enlighten me here?

Thanks in advance
Thomas
--
Dr. Thomas Krebs
Zuken Ltd.
Roonstr. 21
90429 Nrnberg
Germany
Tel.:  ++49 911 9269 111
Fax.: ++49 911 9269 200

Quote:> I have a problem calculating the basis functions for knot vectors that
> are not clamped.
> I use the algorithm from the Piegl/Tiller NURBS book (ALGORITHM A2.2, page 70)
> Piegl/Tiller always assume that the knot vector is clamped, i.e. knot values
> are of multiplicity order (degree + 1) at the start and the end.
> Unfortunately the data I have to deal with doesn't always looks like that,
> e.g. I have a knot vector for a curve of degree 3 that begins like:
> -7.62939453147204E-06,-7.62939453147204E-06,0.0,0.0,0.00390624999999359 ...
> and even ends with:
> 1.,1.,1.00390624999999,1.00781249999999

I suppose that this are rounding errors:
-7.62939453147204E-06 should be 0.0, and 1.00390624999999 should be 1.0.

It would be better to remove them before trying to find the knot span.

If you have a real unclamped knot vector ( like U=(1,2,3,4,5,6,7,8) )
you can use the the knot insertion algorithm given in the NURBS book (with a slight modification) to extract the parameter interval for which the curve points
are defined.
In the avove example the new knot vector would be U=(4,4,4,4,5,5,5,5).
After this you can normalize the knot vector to get U=(0,0,0,0,1,1,1,1).

My implementation of ExtractCurve looks like this:

/*-------------------------------------------------------------------------
  Finds the knotspan containing 'u' in a knot vector 'U' (degree
  'p','n'+1 knots (see "The NURBS book")), and counts the multiplicity
  of 'u'
-------------------------------------------------------------------------*/
template< class T >
inline int FindSpanMultip( const T u, const int n,  const int p,
                           Ptr< T > U, int *s )
{
  int l,low,mid,high;

/* --- Knotenspanne suchen (von rechts):   */

  if( u == U[n+1] )
  {
    if( u == U[n+p+1] ) // clamped
    {
      *s = p + 1;
      return(n);
    }

    l = n;           // if unclamped
    while( l >= 0 )  //   count multiplicity
    {
      if( U[l] != u ) break;
      l--;
    }
    *s = n+1 - l;
    return( n+1 );   // knot span = n+1
  }

  low = p;
  high = n+1;
  mid = (low+high) / 2 ;  

  while( (u < U[mid]) || (u >= U[mid+1]) )
  {
    if( u < U[mid] )
      high = mid;
    else
      low = mid;

    mid = (low+high) / 2;
  }

  l = mid;
  while( l >= 0 )
  {
    if( U[l] != u ) break;
    l--;
  }
  *s = mid - l;

  return(mid);    

Quote:}                                  

  Ptr< HP > Cw;
  Ptr< T > Uh;

  a = FindSpanMultip( u1, n, p, U, &sa );
  if( sa > p )          /* Klammerknoten am Anfang */
    ra = 0;
  else
    ra = p - sa;

  b = FindSpanMultip( u2, n, p, U, &sb );
  if( sb > p )         /* Klammerknoten am Ende */
  {
   sb = p;
   rb = 0;
   b += p;
  }
  else
    rb = p - sb;

  N = (b-sb) - a + p;

  Cw.reserve_pool( N+1 );
  Uh.reserve_pool( N+p+2 );

/* ---- Linken Knoten einfuegen, bis Vielfachheit = p ------------------- */

  for( i = 0; i <= p; i++ )
  {
    Cw[i] = Pw[a-p+i];
    Uh[i] = u1;
  }

  for( r = 1; r <= ra; r++ )
  {
    for( i = 0; i <= ra-r; i++ )
    {
      alpha = (u1 - U[a-p+r+i]) / (U[a+r+i] - U[a-p+r+i]);
      v1 = ((T)1.0-alpha) * Cw[i];
      v2 = alpha * Cw[i+1];
      Cw[i] =  v1 + v2;
    }
  }    

/* ---- Unveraenderten Teil uebernehmen ----------------------------------- */

  for( i = 0; i < (b-sb)-a; i++ )
  {
    Uh[p+1+i] = U[a+1+i];
    Cw[p+1+i] = Pw[a+1+i];
  }

/* ---- Rechten Knoten einfuegen, bis Vielfachheit = p --------------------- */  

  for( r = 1; r <= rb; r++ )
  {
    for( i = 0; i <= rb-r; i++ )
    {
      alpha = (u2 - U[b-sb-i]) / (U[b-sb-i+p] - U[b-sb-i]);
      v1 = alpha * Cw[N-i];
      v2 = ((T)1.0-alpha) * Cw[N-i-1];
      Cw[N-i] = v1 + v2;      
    }
  }

  for( i = 0; i <= p; i++ )
    Uh[N+1+i] = u2;

  NormalizeKnotVector( N, p, Uh );

  *retn = N;
  *retPw = Cw;
  *retU = Uh;  

Quote:}        

  http://www.veryComputer.com/

(i couldn't resist to do this piece of shameless adverti*t, sorry :)))

Regards,
Norbert