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). */