There exists a typo in my Siggraph99 paper.
At three lines below Eq. (2), B_k(...) in the denominator of beta_ij
should be squared: B_k(...) should be replaced by B_k^2(...).
Chang Hwan Seol at KAIST pointed out my mistake.
The below is my single-level B-spline approximation codes.


// ---------------  QmApproximate.h  --------------------------//

class QmScatteredData
{
  public:
	int	position;
	vector	value;
};

class QmApproximate
{
  private:
	int	num;
	vector*	cp;

  public:
	QmApproximate() { num=0; cp=NULL; }

	int	getSize() const { return num; }
	void	setSize( int n );

	vector	evaluate( m_real ) const;
	void	discretize( int, vector* );
	void	approximate( int, int, QmScatteredData *data );

	static m_real B0( m_real t ) { return (1.0-t)*(1.0-t)*(1.0-t) / 6.0; }
	static m_real B1( m_real t ) { return (3.0*t*t*t - 6.0*t*t + 4.0) / 6.0; }
	static m_real B2( m_real t ) { return (-3.0*t*t*t + 3.0*t*t + 3.0*t + 1.0) / 6.0; }
	static m_real B3( m_real t ) { return t*t*t / 6.0; }
};




// ---------------  QmApproximate.cc  --------------------------//

#include "MATHCLASS/mathclass.h" // it defines the "vector" class
#include "qmApproximate.h"

void
QmApproximate::setSize( int n )
{
	if ( n==num ) return;

	delete[] cp; cp = NULL;

	num = n;
	if ( n>0 ) cp = new vector[n+3];
}

vector
QmApproximate::evaluate( m_real t ) const
{
	int i = t;
	if ( i < 0     ) i = 0;
	if ( i > num-1 ) i = num-1;

	m_real f = t - i;

	return B0(f)*cp[i] + B1(f)*cp[i+1] + B2(f)*cp[i+2] + B3(f)*cp[i+3];
}

void
QmApproximate::discretize( int num_frames, vector *data )
{
	for( int i=0; i<num_frames; i++ )
	{
		m_real t = num * i;
		data[i] = evaluate( t/num_frames );
	}
}

void
QmApproximate::approximate( int num_frames, int num_data, QmScatteredData *data )
{
	static vectorN weight; weight.setSize( num+3 );

	for( int i=0; i<num+3; i++ )
	{
		cp[i] = vector(0,0,0);
		weight[i] = 0.0;
	}

	for( int j=0; j<num_data; j++ )
	{
        	m_real t = num * data[j].position;
        	t /= num_frames;

        	i = t;
        	if ( i < 0     ) i = 0;
        	if ( i > num-1 ) i = num-1;

		m_real f = t - i;

		m_real r0 = B0(f); m_real s0 = r0*r0;
		m_real r1 = B1(f); m_real s1 = r1*r1;
		m_real r2 = B2(f); m_real s2 = r2*r2;
		m_real r3 = B3(f); m_real s3 = r3*r3;

		m_real sum = s0 + s1 + s2 + s3;

		cp[i  ] += s0 * data[j].value * r0 / sum;
		cp[i+1] += s1 * data[j].value * r1 / sum;
		cp[i+2] += s2 * data[j].value * r2 / sum;
		cp[i+3] += s3 * data[j].value * r3 / sum;

		weight[i  ] += s0;
		weight[i+1] += s1;
		weight[i+2] += s2;
		weight[i+3] += s3;
	}

	for( i=0; i<num+3; i++ )
		if ( weight[i]>EPS )
			cp[i] /= weight[i];
		else
			cp[i] = vector(0,0,0);
}