basis functions for certain NU(R)B knot vectors

basis functions for certain NU(R)B knot vectors

Post by Dr. Thomas Kreb » Sat, 01 Apr 2000 04:00:00



Hello

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

 
 
 

basis functions for certain NU(R)B knot vectors

Post by Norbert Irme » Sun, 02 Apr 2000 04:00:00


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

Hi,

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:}                                  

/*----------------------------------------------------------------------
  Extracts a part of a NURBS curve (ranging from 'u1' to 'u2')
----------------------------------------------------------------------*/  
template< class T, class HP >
void ExtractCurve(
         const T u1, const T u2,
         const int n, const int p,
         Ptr< T > U, Ptr< HP > Pw,
         int *retn, Ptr< T > *retU, Ptr< HP > *retPw )
{
  int a,ra,sa,b,rb,sb,N,i,r;
  HP v1,v2;
  T alpha;

  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:}        

The above source code and more functions from the NURBS book can be found at

  http://www.veryComputer.com/

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

Regards,
Norbert

 
 
 

1. nu(r) bs question.

Hello there,

I am using the openGL nurbs functions to check a nonuniform b-spline drawing
function I programmed, and it turns out that when drawing the curve with my
own routine, it doesn't draw the last point (so also the line to that
point).

Ik think the problem lies in the following. When drawing a segment, I let a
variable u vary from 0.0 to 1.0 with a certain delta. Now when I reach u=1.0
I skip to the next segment and calculate the point at u=0.0 (which should be
the same as u=1.0 (or 0.999???) of the previous segment.

So the problem lies in the last point of the last segment. How do I
calculate this?

 tried using u=0.99999 and it works (I think),  but is there a better
mathematical solution, instead of this quick hack???

Thanks
  Roger

2. particle question

3. radial basis function

4. SOS

5. Laws' basis functions

6. Importing GIF into Corel Depth

7. NURBS basis function on SGI's

8. Creative fonts

9. Bartel's basis evaluation function

10. Derivative of B-Spline basis function..........

11. radial basis function

12. Fast interpolation using radial basis functions

13. Fitting radial basis function to points cloud