/*
* Copyright (C) 2010-2014 - Andreas Maier
* CONRAD is developed as an Open Source project under the GNU General Public License (GPL).
*/package edu.stanford.rsl.conrad.geometry.splines.fitting;
import java.util.ArrayList;
import java.util.Arrays;
import edu.stanford.rsl.apps.gui.opengl.BSplineVolumeRenderer;
import edu.stanford.rsl.conrad.geometry.shapes.simple.PointND;
import edu.stanford.rsl.conrad.geometry.splines.SurfaceBSpline;
import edu.stanford.rsl.conrad.utils.Configuration;
/**
* Class to fit a grid of data points with an arbitrary degree tensor product B-Spline surface.
* See L Piegl. "The NURBS book". 2012, page 376
* @author Xinyun Li
*
*/
public class BSplineSurfaceInterpolation {
public static final int UNIFORM = 0;
public static final int CHORD = 1;
public static final int CENTRIPETAL = 2;
private ArrayList<PointND> dataPoints;
protected int dimension;
protected int degreeU;//degree
protected int degreeV;
protected int m; //number of columns along u direction
protected int n; //number of rows along v direction
protected double[] uParameters;
protected double[] vParameters;
/**
* Creates a new instance of BSplineSurfaceInterpolation
* @param dataPoints a list of points to be fitted to a BSpline surface,
* actually is reshaped from a grid of points row by row into one list
* @param degreeU degree in the u direction
* @param degreeV degree in the v direction
* @param m number of columns of the grid along u direction
* @param n number of rows of the grid along v direction
*/
public BSplineSurfaceInterpolation(ArrayList<PointND> dataPoints, int degreeU, int degreeV, int m, int n){
this.dataPoints = dataPoints;
this.degreeU = degreeU;
this.degreeV = degreeV;
this.dimension = dataPoints.get(0).getDimension();
this.m = m;
this.n = n;
}
/**
* core method, apply the interpolation, which fits the given data grid to a B-Spline surface of the degree given in the Constructor
* Smooth close is not supported, in order to achieve this, duplicate the first point in the end to get a clamped close
* @param parametrization three parameterization options : uniform, chord length and centripetal
* @return the fitted B-Spline surface
*/
public SurfaceBSpline applyInterpolation(int parametrization){
//only support clamped knots
int endpointCondition = 0;
uParameters = buildUParameterVector(parametrization, endpointCondition);
vParameters = buildVParameterVector(parametrization, endpointCondition);
double[] uKnots = buildKnotVector(uParameters, degreeU, endpointCondition);
System.out.println("uKnot:" + Arrays.toString(uKnots));
double[] vKnots = buildKnotVector(vParameters, degreeV, endpointCondition);
System.out.println("vKnot:" + Arrays.toString(vKnots));
ArrayList<PointND> uControlPoints = new ArrayList<PointND>();
for(int l = 0; l < m; l ++){
ArrayList<PointND> uPoints = new ArrayList<PointND>();
for(int k = 0; k < n; k ++)
uPoints.add(dataPoints.get(k * m + l));
//System.out.println("Udata"+uPoints);
BSplineCurveInterpolation uSplineInterpolation = new BSplineCurveInterpolation(uPoints, degreeU);
ArrayList<PointND> temp = uSplineInterpolation.applyInterpolation(uParameters, uKnots, endpointCondition).getControlPoints();
//System.out.println("UCP_inter"+temp);
uControlPoints.addAll(temp);
}
//uCP should be saved as columns, but as they were saved as arraylist ,they should be read via columns
ArrayList<PointND> ControlPoints = new ArrayList<PointND>();
for(int k = 0; k < n; k ++){
ArrayList<PointND> vPoints = new ArrayList<PointND>();
for(int l = 0; l < m; l ++)
vPoints.add(uControlPoints.get(l * n + k));
// System.out.println("VCP_inter"+vPoints);
BSplineCurveInterpolation vSplineInterpolation = new BSplineCurveInterpolation(vPoints, degreeV);
ArrayList<PointND> temp = vSplineInterpolation.applyInterpolation(vParameters, vKnots, endpointCondition).getControlPoints();
//System.out.println("VCP_final"+temp);
ControlPoints.addAll(temp);
}
return new SurfaceBSpline(ControlPoints, uKnots, vKnots);
}
protected double[] buildUParameterVector(int parametrization, int endpointCondition){
double[] uParameters = new double[n];
uParameters[n - 1] = 1d;
int num = m;
for(int l = 0; l < m; l++){
double total = 0d;
double[] cds = new double[n];
for (int k = 1; k < n; k ++){
cds[k] = dataPoints.get(k * m + l).euclideanDistance(dataPoints.get((k - 1) * m + l));
total += cds[k];
}
if (total == 0) num--;
else {
double d = 0d;
for(int k = 1; k < n - 1; k ++){
d += cds[k];
uParameters[k] += d / total;
}
}
}
if(num == 0) System.out.println("error");
for(int k = 1; k < n - 1; k ++) uParameters[k] /= num;
return uParameters;
}
protected double[] buildVParameterVector(int parametrization, int endpointCondition){
double[] vParameters = new double[m];
vParameters[m - 1] = 1d;
int num = n;
for (int k = 0; k < n; k++){
double total = 0d;
double[] cds = new double[m];
for (int l = 1; l < m; l ++){
cds[l] = dataPoints.get(k * m + l).euclideanDistance(dataPoints.get(k * m + l - 1));
total += cds[l];
}
if (total == 0) num--;
else{
double d = 0d;
for (int l = 1; l < m - 1; l ++){
d += cds[l];
vParameters[l] += d / total;
}
}
}
if(num == 0) System.out.println("error");
for(int l = 1; l < m - 1; l ++) vParameters[l] /= num;
return vParameters;
}
protected double[] buildKnotVector(double[] parameters, int degree, int endpointCondition){
int num = parameters.length;
if (endpointCondition == 1){
int k = num + degree + 1; //number of knots
double[] uKnots = new double[k];
uKnots[num] = 1;
for (int j = 1; j <= num - degree - 1; j++){
double sum = 0;
for (int i = j; i < j + degree; i++)
sum += parameters[i];
uKnots[j + degree] = sum / degree;
}
for(int i = 0; i < degree; i++){
uKnots[degree - i - 1] = uKnots[degree - i] - (uKnots[num - i] - uKnots[num - i - 1]);
uKnots[num + i + 1] = uKnots[num + i] + (uKnots[degree + i + 1] - uKnots[degree + i]);
}
return uKnots;
} else if (endpointCondition == 2){
//uniform knots
int k = num + degree * 2 + 1; //number of knots
double[] knots = new double[k];
for(int j = 0; j < k; j++)
knots[j] = ((double)(j - degree)) / num;
return knots;
} else {
int k = num + degree + 1; //number of knots
double[] Knots = new double[k];
// compute the knot vector
for(int i = 0; i <= degree; i++){
Knots[i] = 0;
Knots[k-i-1] = 1;
}
for(int j = 1; j <= num - degree - 1; j++){
double sum = 0;
for(int i = j; i < j + degree; i++)
sum += parameters[i];
Knots[j + degree] = sum / degree;
}
return Knots;
}
}
public double[] getUParameters() {return uParameters;}
public double[] getVParameters() {return vParameters;}
public static void main(String[] args) {
// TODO Auto-generated method stub
Configuration.loadConfiguration();
//generate a semi-sphere samples
ArrayList<PointND> list = new ArrayList<PointND>();
int r = 1;
list.add(new PointND (0, r, 1));
list.add(new PointND (0.707 * r, 0.707 * r, 1));
list.add(new PointND (r, 0, 1));
list.add(new PointND (0.707 * r, -0.707 * r, 1));
list.add(new PointND (0, -r, 1));
list.add(new PointND (-0.707 * r, -0.707 * r, 1));
list.add(new PointND (-r, 0, 1));
list.add(new PointND (-0.707 * r, 0.707 * r, 1));
list.add(new PointND (0, r, 1));
r = 2;
list.add(new PointND (0, r, 2));
list.add(new PointND (0.707 * r, 0.707 * r, 2));
list.add(new PointND (r, 0, 2));
list.add(new PointND (0.707 * r, -0.707 * r, 2));
list.add(new PointND (0, -r, 2));
list.add(new PointND (-0.707 * r, -0.707 * r, 2));
list.add(new PointND (-r, 0, 2));
list.add(new PointND (-0.707 * r, 0.707 * r,2));
list.add(new PointND (0, r, 2));
r = 3;
list.add(new PointND (0, r, 3));
list.add(new PointND (0.707 * r, 0.707 * r, 3));
list.add(new PointND (r, 0, 3));
list.add(new PointND (0.707 * r, -0.707 * r, 3));
list.add(new PointND (0, -r, 3));
list.add(new PointND (-0.707 * r, -0.707 * r, 3));
list.add(new PointND (-r, 0, 3));
list.add(new PointND (-0.707 * r, 0.707 * r, 3));
list.add(new PointND (0, r, 3));
r = 5;
list.add(new PointND (0, r, 5));
list.add(new PointND (0.707 * r, 0.707 * r, 5));
list.add(new PointND (r, 0, 5));
list.add(new PointND (0.707 * r, -0.707 * r, 5));
list.add(new PointND (0, -r, 5));
list.add(new PointND (-0.707 * r, -0.707 * r, 5));
list.add(new PointND (-r, 0, 5));
list.add(new PointND (-0.707 * r, 0.707 * r, 5));
list.add(new PointND (0, r, 5));
//B-Spline surface interpolation
int degreeU = 2, degreeV = 2;
int m = 9, n = 4;
BSplineSurfaceInterpolation test = new BSplineSurfaceInterpolation(list, degreeU, degreeV, m, n);
SurfaceBSpline Sbsp = test.applyInterpolation(1);
for(int k = 0; k < test.uParameters.length; k++) {
for(int j = 0; j < test.vParameters.length; j++) {
PointND res = Sbsp.evaluate(test.uParameters[k], test.vParameters[j]);
System.out.println("RES at ("+test.uParameters[k]+ ", " + test.vParameters[j] + ") is "+res);
}
}
//Visualization
BSplineVolumeRenderer test1 = new BSplineVolumeRenderer(Sbsp);
}
}