content: package edu.stanford.rsl.conrad.geometry.splines;
  
  import java.util.ArrayList;
  import java.util.Arrays;
  import java.util.Collections;
  import java.util.HashMap;
  import java.util.Iterator;
  import java.util.Map;
  
  import edu.stanford.rsl.conrad.geometry.AbstractCurve;
  import edu.stanford.rsl.conrad.geometry.AbstractShape;
  import edu.stanford.rsl.conrad.geometry.shapes.simple.PointND;
  import edu.stanford.rsl.conrad.geometry.transforms.Transform;
  import edu.stanford.rsl.conrad.numerics.SimpleVector;
  
  public class BSpline extends AbstractCurve {
  	
  	boolean IS_CLAMPED = true;
  	private double lowerBound = 0;
  	private double delta = 1;
  
  	/**
  	 * 
  	 */
  	public static final float BSPLINECOLLECTION = 0;
  	public static final float SPLINE1DTO3D = 1;
  	public static final float SPLINE2DTO3D = 2;
  	public static final float SPLINE3DTO3D = 3;
  	private static Map<String, Float> idMap = Collections.synchronizedMap(new HashMap<String, Float>());
  	
  	protected static final long serialVersionUID = -1510774879340694844L;
  	
  	protected int degree;
  	private double [] knots = null;
  	private ArrayList<PointND> controlPoints;
  	protected int dimension;
  	
  	public BSpline (ArrayList<PointND> controlPoints, double ... uVector){
  		this.setControlPoints(controlPoints);
  		setKnots(uVector);
  		degree = ((getKnots().length - controlPoints.size()));
  		checkIfClamped();
  	}
  	
  	public BSpline (ArrayList<PointND> controlPoints, int degree, double ... uVector){
  		this.setControlPoints(controlPoints);
  		setKnots(uVector);
  		this.degree = degree;
  		checkIfClamped();
  	}
  	
  	public BSpline(BSpline bs){
  		super(bs);
  		degree = bs.degree;
  		dimension = bs.dimension;
  		if (bs.knots != null){
  			knots = new double[bs.knots.length];
  			System.arraycopy(bs.knots, 0, knots, 0, bs.knots.length);
  		}
  		else
  			knots = null;
  		
  		if (bs.controlPoints != null){
  			Iterator<PointND> it = bs.controlPoints.iterator();
  			controlPoints = new ArrayList<PointND>();
  			while (it.hasNext()) {
  				PointND p = it.next();
  				controlPoints.add((p!=null) ? p.clone() : null);
  			}
  		}
  		else{
  			controlPoints = null;
  		}
  		checkIfClamped();
  	}
  	
  	private void checkIfClamped(){
  		for(int i = 0; i < degree; i++){
  			if( (knots[i] != 0) && (knots[knots.length-1-i] != 1) ){
  				this.IS_CLAMPED = false;
  				lowerBound =   knots[degree];
  				delta = (knots[knots.length-degree]-lowerBound);
  				break;
  			}				
  		}
  	}
  	
  	/**
  	Non-recursive implementation of the N-function.
  	*/
  	protected double N(double internalCoordinate, int index) {
  		//double [] knot = getKnots();
  		double d = 0;
  		int [] coefficientArrayA = new int[2 * degree];
  		int [] coefficientArrayC = new int[2 * degree];
  		int degreeMinus2 = degree - 2;
  		for (int j = 0; j < degree; j++) {
  			double firstKnot = knots[index+j];
  			double secondKnot = knots[index+j+1];
  			if (internalCoordinate >= firstKnot && internalCoordinate <= secondKnot && firstKnot != secondKnot) {
  				// set required coefficients to 0
  				for (int k = degree - j - 1; k >= 0; k--){
  					coefficientArrayA[k] = 0;
  				}
  				if (j > 0) {
  					// set required coefficients to their index
  					for (int k = 0; k < j; k++){
  						coefficientArrayC[k] = k;
  					}
  					coefficientArrayC[j] = Integer.MAX_VALUE;
  				}
  				else {
  					coefficientArrayC[0] = degreeMinus2;
  					coefficientArrayC[1] = degree;
  				}
  				int z = 0;
  				while (true) {
  					if (coefficientArrayC[z] < coefficientArrayC[z+1] - 1) {
  						double e = 1.0;
  						int bCounter = 0;
  						int y = degreeMinus2 - j;
  						int p = j - 1;
  
  						for (int m = degreeMinus2, n = degree; m >= 0; m--, n--) {
  							if (p >= 0 && coefficientArrayC[p] == m) {
  								int w = index + bCounter;
  								double kd = knots[w+n];
  								e *= (kd - internalCoordinate) / (kd - knots[w+1]);
  								bCounter++;
  								p--;
  							}
  							else {
  								int w = index + coefficientArrayA[y];
  								double kw = knots[w];
  								e *= (internalCoordinate - kw) / (knots[w+n-1] - kw);
  								y--;
  							}
  						}
  						// this code updates the a-counter
  						if (j > 0) {
  							int g = 0;
  							boolean reset = false;
  
  							while (true) {
  								coefficientArrayA[g]++;
  
  								if (coefficientArrayA[g] > j) {
  									g++;
  									reset = true;
  								}
  								else {
  									if (reset) {
  										for (int h = g - 1; h >= 0; h--)
  											coefficientArrayA[h] = coefficientArrayA[g];
  									}
  									break;
  								}
  							}
  						}
  
  						d += e;
  
  						// this code updates the bit-counter
  						coefficientArrayC[z]++;
  						if (coefficientArrayC[z] > degreeMinus2) break;
  
  						for (int k = 0; k < z; k++)
  							coefficientArrayC[k] = k;
  						z = 0;
  					}
  					else {
  						z++;
  					}
  				}
  				break; // required to prevent spikes
  			}
  		}
  
  		return d;
  	}
  	
  	public double [] evalFast(double t) {
  		if(!IS_CLAMPED){
  			t = lowerBound + t*delta;
  		}
  		double [] p = new double [getControlPoints().get(0).getDimension()];
  		int context = degree+1;
  		int index = Arrays.binarySearch(getKnots(), t);
  		int numPts = getControlPoints().size();
  		if (index < 0) index *= -1; 
  		int start = index - context;
  		int end = index + context;
  		if (start < 0) start = 0;
  		if (end > numPts) end = numPts;
  		//System.out.print("t = " + t +" ");
  		for (int i = start; i < end; i++) {
  			double w = N(t, i);
  			//double w = N(t, i, degree);
  			double[] loc = getControlPoints().get(i).getCoordinates();
  			for (int j = 0; j < loc.length; j++){
  				p[j] += (loc[j] * w); //pt[i][j] * w);
  				//if (j==0) System.out.print("weight " + w + " ");
  			}
  		}
  		//System.out.println(p [0] + " " + p[1]);
  		return p;
  	}
  	
  	/**
  	 * Constructor for a BSpline using ArbitraryPoints and a weight vector as SimpleVector
  	 * @param controlPoints the control points
  	 * @param knotVector the weight vector
  	 */
  	public BSpline (ArrayList<PointND> controlPoints, SimpleVector knotVector){
  		this(controlPoints, new PointND(knotVector).getCoordinates());
  	}
  	
  	@Override
  	public PointND evaluate(double u) {
  		double [] point = evalFast(u);
  		return new PointND(point);
  	}
  	
  	
  
  	@Override
  	public int getDimension() {
  		return dimension;
  	}
  	
  	/**
  	 * Computes N iteratively.
  	 * @param u the point on the curve
  	 * @param i the control point
  	 * @return the weight
  	 */
  	public double getWeight(double u, int i){
  		return N(u, i);
  	}
  
  	/**
  	 * Computes the bounding box for the curve.
  	 */
  	public void computeBounds(){
  		min = getControlPoints().get(0);
  		max = getControlPoints().get(0);
  		for (int i = 1; i < getControlPoints().size(); i++){
  			min.updateIfLower(getControlPoints().get(i));
  			max.updateIfHigher(getControlPoints().get(i));
  		}
  	}
  	
  	/**
  	 * Returns the i-th control point
  	 * @param i the index
  	 * @return the control point
  	 */
  	public PointND getControlPoint(int i){
  		return getControlPoints().get(i);
  	}
  	
  	/**
  	 * Returns the i-th knot vector entry.
  	 * 
  	 * @param i
  	 * @return the knot vector entry
  	 */
  	public double getKnotVectorEntry(int i){
  		return getKnots()[i];
  	}
  	
  	/**
  	 * Returns the degree of the spline.
  	 * @return the degree
  	 */
  	public int getDegree(){
  		return degree;
  	}
  	
  	@Override
  	public ArrayList<PointND> intersect(AbstractCurve other) {
  		throw new RuntimeException("intersection between BSplines and other curves are not yet implemented");
  	}
  
  	@Override
  	public boolean isBounded() {
  		return true;
  	}
  
  	public double distance(double [] planeCoefficients, double u) {
  		return distanceFast(planeCoefficients, u);
  	}
  	
  	public double distance(double [] planeCoefficients, int i) {
  		return distanceFast(planeCoefficients, i);
  	}
  	
  	public double distanceFast(double[] plane, double u) {
  		int numPts = getControlPoints().size();
  		//System.out.println("GroupSize: " + numPts);
  		int context = degree + 2;
  		int index = Arrays.binarySearch(getKnots(), u);
  		if (index < 0) index *= -1; 
  		//System.out.println(gi.getGroupSize() + " " + index + " " + t);
  		int start = index - context;
  		int end = index + context;
  		if (start < 0) start = 0;
  		if (end > numPts) end = numPts;
  		double revan = 0;
  		for (int i = start; i < end; i++) {
  			double w = N(u, i);
  			//double w = N(t, i, degree);
  			double[] loc = getControlPoints().get(i).getCoordinates();
  			revan += (loc[0] * plane[0] +
  					loc[1] * plane[1] +
  					loc[2] * plane[2] +
  					plane[3]) * w;
  		}
  		return revan;
  	}
  
  	@Override
  	public void applyTransform(Transform t) {
  		ArrayList<PointND> newPoints = new ArrayList<PointND>();
  		for (int i = 0; i < getControlPoints().size(); i++){
  			newPoints.add(t.transform(getControlPoints().get(i)));
  		}
  		setControlPoints(newPoints);
  	}
  
  	@Override
  	public PointND[] getRasterPoints(int number) {
  		PointND [] array = new PointND[number];
  		for (double i = 0; i <number ; i++){
  			array[(int)(i)] = evaluate(i/(number-1));
  		}
  		return array;
  	}
  
  	public static float getID(String title) {
  		Float result = idMap.get(title);
  		if (result == null){
  			result = new Float(idMap.keySet().size());
  			idMap.put(title, result);
  		}
  		return result;
  	}
  
  	/**
  	 * Rewrites the BSpline into a float representation:
  	 * <pre>
  	 * type
  	 * total size
  	 * degree
  	 * # of entries in knot vector
  	 * # of control points
  	 * knovector entries (1 float value)
  	 * control points (3 float values: x,y,z)
  	 * </pre>
  	 * @return the binary represenation for use with openCL
  	 */
  	public float[] getBinaryRepresentation() {
  		int totalsize = this.getKnots().length + (this.getControlPoints().size()*3) + 5;
  		float [] binary = new float [totalsize];
  		int index = 0;
  		binary[index] = SPLINE1DTO3D;
  		index ++;
  		binary[index] = totalsize;
  		index ++;
  		binary[index] = this.degree;
  		index ++;
  		binary[index] = this.getKnots().length;
  		index ++;
  		binary[index] = this.getControlPoints().size();
  		index ++;
  		for (int i = 0; i < getKnots().length; i++){
  			binary[index] = (float)getKnots()[i];
  			index++;
  		}
  		for (int i = 0; i < getControlPoints().size(); i++){
  			binary[index] = (float) getControlPoints().get(i).get(0);
  			index++;
  			binary[index] = (float) getControlPoints().get(i).get(1);
  			index++;
  			binary[index] = (float) getControlPoints().get(i).get(2);
  			index++;
  		}
  		return binary;
  	}
  	
  	/**
  	 * Calculates the derivative of a spline
  	 * @return
  	 */
  	public BSpline getDerivative(){
  		int dim = controlPoints.get(0).getDimension();
  		// implementation following James D Emery: "B-Splines And Divided Differences"
  		double[] knots = this.knots;
  		
  		ArrayList<PointND> q = new ArrayList<PointND>();
  		
  		for(int i = 0; i < this.controlPoints.size()+1; i++){
  			double factor = degree-1;
  			factor /= (knots[i+degree-1]-knots[i]);
  			
  			SimpleVector pHigh = (i==controlPoints.size()) ? 
  					new PointND(new double[dim]).getAbstractVector() :
  						new SimpleVector(controlPoints.get(i).getAbstractVector());
  			SimpleVector pLow = (i==0) ? new PointND(new double[dim]).getAbstractVector() :
  					controlPoints.get(i-1).getAbstractVector();
  			pHigh.subtract(pLow);
  			q.add(new PointND(pHigh.multipliedBy(factor).copyAsDoubleArray()));
  		}
  		q.remove(q.size()-1);
  		q.remove(0);
  		knots = new double[this.knots.length-2];
  		for(int i = 0; i < knots.length; i++){
  			knots[i] = (IS_CLAMPED) ? this.knots[i+1]:(i/(knots.length-1));
  		}
  		
  		return new BSpline(q,knots);
  	}
  
  	public void setKnots(double [] knots) {
  		this.knots = knots;
  	}
  
  	public double [] getKnots() {
  		return knots;
  	}
  
  	/**
  	 * @param controlPoints the controlPoints to set
  	 */
  	public void setControlPoints(ArrayList<PointND> controlPoints) {
  		this.controlPoints = controlPoints;
  	}
  
  	/**
  	 * @return the controlPoints
  	 */
  	public ArrayList<PointND> getControlPoints() {
  		return controlPoints;
  	}
  
  	@Override
  	public AbstractShape clone() {
  		return new BSpline(this);
  	}
  	
  }
  /*
   * Copyright (C) 2010-2014 Andreas Maier
   * CONRAD is developed as an Open Source project under the GNU General Public License (GPL).
  */