/********************************************************************* Simple b-spline curve algorithm Copyright 1994 by Keith Vertanen ([email protected]) Released to the public domain (your mileage may vary)**********************************************************************/#ifndef __BSP_H__ #define __BSP_H__struct point { double x; double y; double z; };void compute_intervals(int *u, int n, int t); double blend(int k, int t, int *u, double v); void compute_point(int *u, int n, int t, double v, point *control, point *output); void bspline(int n, int t, point *control, point *output, int num_output); void makebspline(point *output, int num_output);#endif //#include "stdafx.h" #include "Bsp.h"/********************************************************************* Parameters: n - the number of control points minus 1 t - the degree of the polynomial plus 1 control - control point array made up of point stucture output - array in which the calculate spline points are to be put num_output - how many points on the spline are to be calculatedPre-conditions: n+2>t (no curve results if n+2<=t) control array contains the number of points specified by n output array is the proper size to hold num_output point structures **********************************************************************/ void bspline(int n, int t, point *control, point *output, int num_output) { int *u; double increment,interval; point calcxyz; int output_index;
u=new int[n+t+1];//7+4+1=12 compute_intervals(u, n, t); //==int u[12]={0,0,0,0,1,2,3,4,5,5,5,5}; // 7-4+2=5 smooth=32 increment=(double) (n-t+2)/(num_output-1);//=0.16129 how much parameter goes up each time interval=0;
for (output_index=0; output_index<num_output-1; output_index++) { compute_point(u, n, t, interval, control, &calcxyz); output[output_index].x = calcxyz.x; output[output_index].y = calcxyz.y; output[output_index].z = calcxyz.z; interval=interval+increment; // increment our parameter 0.16129*31=4.999 } output[num_output-1].x=control[n].x; // put in the last point output[num_output-1].y=control[n].y; output[num_output-1].z=control[n].z; // delete u; } // compute_point() calls it double blend(int k, int t, int *u, double interval) // calculate the blending value { double value;
if (t==1) // base case for the recursion { if ((u[k]<=interval) && (interval<u[k+1])) value=1; else value=0; } else {// t <> 1 if ((u[k+t-1]==u[k]) && (u[k+t]==u[k+1])) // check for divide by zero value = 0; else if (u[k+t-1]==u[k]) // if a term's denominator is zero,use just the other value = (u[k+t] - interval) / (u[k+t] - u[k+1]) * blend(k+1, t-1, u, interval);//recursive else if (u[k+t]==u[k+1]) value = (interval - u[k]) / (u[k+t-1] - u[k]) * blend(k, t-1, u, interval); else value = (interval - u[k]) / (u[k+t-1] - u[k]) * blend(k, t-1, u, interval) + (u[k+t] - interval) / (u[k+t] - u[k+1]) * blend(k+1, t-1, u, interval); } return value; } //1 k t u v //blend(3, 1, int * 0x00345678, double 0.0) //blend(2, 2, int * 0x00345678, double 0.0) //blend(1, 3, int * 0x00345678, double 0.0) //blend(0, 4, int * 0x00345678, double 0.0) //2 //blend(2, 2, int * 0x00345678, double 0.0) //blend(1, 3, int * 0x00345678, double 0.0) //blend(0, 4, int * 0x00345678, double 0.0) //3 //blend(1, 3, int * 0x00345678, double 0.0) //blend(0, 4, int * 0x00345678, double 0.0) //4 //blend(0, 4, int * 0x00345678, double 0.0) // void compute_intervals(int *u, int n, int t)// figure out the knots {//n=7,t=4 int j;
for (j=0; j<=n+t; j++)//7+4 { if (j<t) u[j]=0;//<4 =[0]=[1]=[2]=[3]=0 else if ((t<=j) && (j<=n)) u[j]=j-t+1;//>=4 ;<=7 =[4]=1 [5]=2 [6]=3 [7]=4 else if (j>n) u[j]=n-t+2;//>7 =[8]=5 [9]=5 [10]=5 [11]=5 // if n-t=-2 then we're screwed, everything goes to 0 } } // void compute_point(int *u, int n, int t, double interval, point *control, point *output) {// 7 4 0->31*0.16129 int k; double temp;
// initialize the variables that will hold our outputted point output->x=0; output->y=0; output->z=0;
for (k=0; k<=n; k++) { temp = blend(k,t,u,interval); // same blend is used for each dimension coordinate output->x = output->x + (control[k]).x * temp; output->y = output->y + (control[k]).y * temp; output->z = output->z + (control[k]).z * temp; } }
#define __BSP_H__struct point
{
double x;
double y;
double z;
};void compute_intervals(int *u, int n, int t);
double blend(int k, int t, int *u, double v);
void compute_point(int *u, int n, int t, double v, point *control, point *output);
void bspline(int n, int t, point *control, point *output, int num_output);
void makebspline(point *output, int num_output);#endif
//#include "stdafx.h"
#include "Bsp.h"/*********************************************************************
Parameters:
n - the number of control points minus 1
t - the degree of the polynomial plus 1
control - control point array made up of point stucture
output - array in which the calculate spline points are to be put
num_output - how many points on the spline are to be calculatedPre-conditions:
n+2>t (no curve results if n+2<=t)
control array contains the number of points specified by n
output array is the proper size to hold num_output point structures
**********************************************************************/
void bspline(int n, int t, point *control, point *output, int num_output)
{
int *u;
double increment,interval;
point calcxyz;
int output_index;
u=new int[n+t+1];//7+4+1=12
compute_intervals(u, n, t);
//==int u[12]={0,0,0,0,1,2,3,4,5,5,5,5};
// 7-4+2=5 smooth=32
increment=(double) (n-t+2)/(num_output-1);//=0.16129 how much parameter goes up each time
interval=0;
for (output_index=0; output_index<num_output-1; output_index++)
{
compute_point(u, n, t, interval, control, &calcxyz);
output[output_index].x = calcxyz.x;
output[output_index].y = calcxyz.y;
output[output_index].z = calcxyz.z;
interval=interval+increment; // increment our parameter 0.16129*31=4.999
}
output[num_output-1].x=control[n].x; // put in the last point
output[num_output-1].y=control[n].y;
output[num_output-1].z=control[n].z;
//
delete u;
}
// compute_point() calls it
double blend(int k, int t, int *u, double interval) // calculate the blending value
{
double value;
if (t==1) // base case for the recursion
{
if ((u[k]<=interval) && (interval<u[k+1])) value=1;
else value=0;
}
else
{// t <> 1
if ((u[k+t-1]==u[k]) && (u[k+t]==u[k+1])) // check for divide by zero
value = 0;
else if (u[k+t-1]==u[k]) // if a term's denominator is zero,use just the other
value = (u[k+t] - interval) / (u[k+t] - u[k+1]) * blend(k+1, t-1, u, interval);//recursive
else if (u[k+t]==u[k+1])
value = (interval - u[k]) / (u[k+t-1] - u[k]) * blend(k, t-1, u, interval);
else
value = (interval - u[k]) / (u[k+t-1] - u[k]) * blend(k, t-1, u, interval) +
(u[k+t] - interval) / (u[k+t] - u[k+1]) * blend(k+1, t-1, u, interval);
}
return value;
}
//1 k t u v
//blend(3, 1, int * 0x00345678, double 0.0)
//blend(2, 2, int * 0x00345678, double 0.0)
//blend(1, 3, int * 0x00345678, double 0.0)
//blend(0, 4, int * 0x00345678, double 0.0)
//2
//blend(2, 2, int * 0x00345678, double 0.0)
//blend(1, 3, int * 0x00345678, double 0.0)
//blend(0, 4, int * 0x00345678, double 0.0)
//3
//blend(1, 3, int * 0x00345678, double 0.0)
//blend(0, 4, int * 0x00345678, double 0.0)
//4
//blend(0, 4, int * 0x00345678, double 0.0)
//
void compute_intervals(int *u, int n, int t)// figure out the knots
{//n=7,t=4
int j;
for (j=0; j<=n+t; j++)//7+4
{
if (j<t) u[j]=0;//<4 =[0]=[1]=[2]=[3]=0
else if ((t<=j) && (j<=n)) u[j]=j-t+1;//>=4 ;<=7 =[4]=1 [5]=2 [6]=3 [7]=4
else if (j>n) u[j]=n-t+2;//>7 =[8]=5 [9]=5 [10]=5 [11]=5
// if n-t=-2 then we're screwed, everything goes to 0
}
}
//
void compute_point(int *u, int n, int t, double interval, point *control, point *output)
{// 7 4 0->31*0.16129
int k;
double temp;
// initialize the variables that will hold our outputted point
output->x=0;
output->y=0;
output->z=0;
for (k=0; k<=n; k++)
{
temp = blend(k,t,u,interval); // same blend is used for each dimension coordinate
output->x = output->x + (control[k]).x * temp;
output->y = output->y + (control[k]).y * temp;
output->z = output->z + (control[k]).z * temp;
}
}