// cubic B spline curve fitting
// B-splines
//An order   B-spline is formed by joining several pieces of polynomials of degree   with at most   continuity at the breakpoints. A set of non-descending breaking points   defines a knot vector 
 //T=(to,t1,...,tm)(1) which determines the parametrization of the basis functions.
//Given a knot vector   , the associated B-spline basis functions,   , are defined as: 
// 	 	 	             (1.58)
//
//for   , and 
// 	 	 	             (1.59)
//
//for   and   . These equations have the following properties:
//"	Positivity:   , for   .
//"	Local support:   , for   , and   .
//"	Partition of unity:   , for   .
//"	Recursion: Given by (1.59).
//"	Continuity:   has   continuity at each simple knot.
//The concept of nodes or Greville abscissae [130,92], which are the averages of the knots, are important in B-spline approximations [130,452] and defined as follows: 
// 	 	 	(1.60)
//
//The node   generally lies near the parameter value which corresponds to a maximum of the basis function  [345,314].
//The derivative of the B-spline basis function is given by [314] 
// 	 	 	(1.61)
//B-spline curve
//A B-spline curve is defined as a linear combination of control points   and B-spline basis functions   given by 
// 	 	 	(1.62)
//
//In this context the control points are called de Boor points. The basis function   is defined on a knot vector 
 //	 	 	(1.63)
//
//where there are   elements, i.e. the number of control points   plus the order of the curve   . Each knotspan   is mapped onto a polynomial curve between two successive joints   and   . Normalization of the knot vector, so it covers the interval [0,1], is helpful in improving numerical accuracy in floating point arithmetic computation due to the higher density of floating point numbers in this interval [133,300].
// 
//Figure 1.10: An order four B-spline basis functions with uniform knot vector
// 
//Figure 1.11: A clamped cubic B-spline curve
//A B-spline curve has the following properties:
//"	Geometry invariance property: Partition of unity property of the B-spline assures the invariance of the shape of the B-spline curve under translation and rotation.
//"	End points geometric property:
//o	Unlike Bézier curves, B-spline curves do not in general pass through the two end control points. Increasing the multiplicity of a knot reduces the continuity of the curve at that knot. Specifically, the curve is   times continuously differentiable at a knot with multiplicity     , and thus has   continuity. Therefore, the control polygon will coincide with the curve at a knot of multiplicity   , and a knot with multiplicity  indicates   continuity, or a discontinuous curve. Repeating the knots at the end   times will force the endpoints to coincide with the control polygon. Thus the first and the last control points of a curve with a knot vector described by 
// 	 	 	(1.64)
//o	
//coincide with the endpoints of the curve. Such knot vectors and curves are known as clamped [314]. In other words, clamped/unclamped refers to whether both ends of the knot vector have multiplicity equal to   or not. Figure 1.10 shows cubic B-spline basis functions defined on a knot vector  . A clamped cubic B-spline curve based on this knot vector is illustrated in Fig. 1.11 with its control polygon.
//o	B-spline curves with a knot vector (1.64) are tangent to the control polygon at their endpoints. This is derived from the fact that the first derivative of a B-spline curve is given by [175] 
 // 	 	(1.65)
//o	
//where the knot vector is obtained by dropping the first and last knots from (1.64), i.e. 
// 	 	 	(1.66)
//o	
//and 
// 	 	 	(1.67)
// 	 	 	(1.68)
//o	
//
//"	Convex hull property: The convex hull property for B-splines applies locally, so that a span lies within the convex hull of the control points that affect it. This provides a tighter convex hull property than that of a Bézier curve, as can be seen in Fig. 1.11. The   -th span of the cubic B-spline curve in Fig. 1.11 lies within the convex hull formed by control points  ,   ,   ,   . In other words, a B-spline curve must lie within the union of all such convex hulls formed by   successive control points [130].
//"	Local support property: A single span of a B-spline curve is controlled only by   control points, and any control point affects   spans. Specifically, changing   affects the curve in the parameter range   and the curve at a point   where   is determined completely by the control points   as shown in Fig.1.11.
//"	Variation diminishing property:
//o	2-D: The number of intersections of a straight line with a planar B-spline curve is no greater than the number of intersections of the line with the control polygon. A line intersecting the convex hull of a planar B-spline curve may intersect the curve transversally, be tangent to the curve, or not intersect the curve at all. It may not, however, intersect the curve more times than it intersects the control polygon.
//o	3-D: The same relation holds true for a plane with a 3-D space B-spline curve.
//"	B-spline to Bézier property: From the discussion of end points geometric property, it can be seen that a Bézier curve of order   (degree   ) is a B-spline curve with no internal knots and the end knots repeated   times. The knot vector is thus 
 //	 	 	(1.69)
//"	
//where   or   .
 
import java.io.*;
import javax.swing.*;

class TNAcs
{
	public static double [] hi (double ti[])
{
		int n=ti.length-1;
		int nn=n+2;
		double h[]=new double [nn];
	           for (int i=0;i<n;i++)
               {h[i]=ti[i+1]-ti[i];}
               h[0]=h[1]; h[n]=h[n-1];
               return h;
}
public static double [] B_Spline(double ti[],double yi[])
{
        int n=ti.length-1;
		int nn=n+2;
        int i;
        double delta,gamma,p,q,r;
        double a[]=new double[nn];
        double h[]=new double[nn];
        h=hi(ti);
        delta=-1;
        gamma=3.0*yi[0];
        p=delta*gamma;
        q=3.0;
        for(i=0;i<=n;i++)
        {r=h[i+1]/h[i];
  delta=-r*delta;
  gamma=-r*gamma+(r+1)*yi[i];
  p=p+gamma*delta;
  q=q+delta*delta;
  }
  a[0]=Math.sqrt(p/q);
   for(i=0;i<(n+1);i++)
  {a[i+1]=((h[i]+h[i+1])*yi[i]-h[i+1]*a[i])/h[i];}
  return a;
}     
public static double funcB_Spline(double ti[],double a[],double x)
{
  double h[]=hi(ti);
  int i;
  int n=ti.length-1;
  double nn=n+2;
  double d,e;
  for(i=n-2;i>=0;i--)
    {if((x-ti[i])>=0) break;}
  i++;
  double bs= a[i+2]*(((ti[i+1]-x+h[i+2])*(ti[i+1]-x+h[i+2]))/(h[i+1]+h[i+2])*h[i+2])+a[i+1]*((((x-ti[i]+h[i])*(ti[i+1]-x))/((h[i]+h[i+1])*h[i+1])+(ti[i+1]-x+h[i+2])*(x-ti[i])/((h[i+1]+h[i+2])*h[i+1])))+a[i]*((x-ti[i]+h[i])*(x-ti[i]+h[i]))/((h[i]+h[i+1])*h[i]);
  return bs;
  }
  public static double[][]  funcB_Spline(double ti[],double yi[],int aradegersayisi)
 {
	int n=ti.length;
	int nn=(n-1)*(aradegersayisi+1);
	double z[][]=new double[2][nn];
    double A[]=B_Spline(ti,yi);	
    double dx=0;
    int k=0;
    int i;
    for(i=0;i<(n-1);i++)
       {  z[0][k]=ti[i];z[1][k]=funcB_Spline(ti,A,z[0][k]);k++;
          for(int j=0;j<aradegersayisi;j++) 
             {dx=(ti[i+1]-ti[i])/((double)aradegersayisi+1.0);
          z[0][k]=z[0][k-1]+dx;z[1][k]=funcB_Spline(ti,A,z[0][k]);k++;}
       }
       
    return z;
 } 
 public static void main(String args[]) throws IOException
  {  
  double x[]={0.0,  0.6 ,1.5 ,1.7 ,1.9 ,2.1,2.3 ,2.6 ,2.8,3.0 ,3.6 ,4.7  ,5.2  ,5.7  ,5.8  , 6.0 , 6.4 ,6.9,7.6,8.0};
  double y[]={-0.8,-0.34,0.29,0.59,0.23,0.1,0.28,1.03,1.5,1.44,0.74,-0.82,-1.27,-0.92,-0.92,-1.04,-0.79,-0.06,1.0,0.0};
  int n=x.length;
  double z1[][]=funcB_Spline(x,y,5	); 
  Plot pp=new Plot(z1[0],z1[1]);
  pp.setPlabel("B-spline curve fitting");
  pp.setXlabel("x");
  pp.setYlabel("y");
  pp.setPlotType(0,0); 
  pp.addData(x,y);
  pp.setPlotType(1,24);
  pp.plot();
  }
}
