package edu.stanford.rsl.conrad.geometry.splines;
import java.io.BufferedReader;
import java.io.FileReader;
import java.io.IOException;
import java.util.ArrayList;
import java.util.List;
import edu.stanford.rsl.conrad.geometry.AbstractCurve;
import edu.stanford.rsl.conrad.geometry.AbstractShape;
import edu.stanford.rsl.conrad.geometry.AbstractSurface;
import edu.stanford.rsl.conrad.geometry.General;
import edu.stanford.rsl.conrad.geometry.shapes.compound.CompoundShape;
import edu.stanford.rsl.conrad.geometry.shapes.compound.LinearOctree;
import edu.stanford.rsl.conrad.geometry.shapes.compound.NestedOctree;
import edu.stanford.rsl.conrad.geometry.shapes.simple.Edge;
import edu.stanford.rsl.conrad.geometry.shapes.simple.PointND;
import edu.stanford.rsl.conrad.geometry.shapes.simple.StraightLine;
import edu.stanford.rsl.conrad.geometry.shapes.simple.Triangle;
import edu.stanford.rsl.conrad.geometry.transforms.Transform;
import edu.stanford.rsl.conrad.numerics.SimpleOperators;
import edu.stanford.rsl.conrad.numerics.SimpleVector;
import edu.stanford.rsl.jpop.FunctionOptimizer;
import edu.stanford.rsl.jpop.OptimizableFunction;
/**
* Class to model a surface made of BSplines.
*
*
* @author akmaier
*
*/
public class SurfaceBSpline extends AbstractSurface {
/**
*
*/
private static final long serialVersionUID = 4942618494585588413L;
protected SimpleVector uKnots;
protected SimpleVector vKnots;
protected int numberOfUPoints;
protected int numberOfVPoints;
protected ArrayList<PointND> points;
protected BSpline [] vSplines;
protected BSpline uSpline;
protected int dimension;
protected boolean clockwise;
public static final int TESSELATE_COMPOUND_SHAPE = 0x1;
public static final int TESSELATE_COMPOUND_OF_COMPOUND_SHAPES = 0x2;
public static final int TESSELATE_COMPOUND_OF_OCTREES = 0x3;
public static final int TESSELATE_LINEAR_OCTREE = 0x4;
public static final int TESSELATE_NESTED_OCTREE = 0x5;
private String title;
public ArrayList<PointND> getControlPoints(){
return points;
}
public SimpleVector getUKnots(){
return uKnots;
}
public SimpleVector getVKnots(){
return vKnots;
}
public void setUKnot(SimpleVector uknots){
uKnots = uknots;
}
public void setVKnot(SimpleVector vknots){
vKnots = vknots;
}
public boolean isClockwise(){
return clockwise;
}
public void setClockwise(boolean in_clockwise){
clockwise = in_clockwise;
}
public BSpline[] getVSplines(){
return vSplines;
}
/**
* @return the title
*/
public String getTitle() {
return title;
}
/**
* @param title the title to set
*/
public void setTitle(String title) {
this.title = title;
}
protected synchronized void init(){
dimension = points.get(0).getDimension();
int degreeU = 0;
while(uKnots.getElement(degreeU) == 0) degreeU++;
degreeU --;
int degreeV = 0;
while(uKnots.getElement(degreeV) == 0) degreeV++;
degreeV --;
numberOfUPoints = uKnots.getLen() - ((degreeU+1));
numberOfVPoints = vKnots.getLen() - ((degreeV+1));
assert(numberOfUPoints * numberOfVPoints == points.size());
ArrayList<PointND> temp = new ArrayList<PointND>();
for(int i = 0; i < numberOfUPoints; i++){
temp.add(points.get(0));
}
uSpline = new BSpline(temp, uKnots);
vSplines = new BSpline[numberOfUPoints];
for(int i = 0; i < numberOfUPoints; i++){
temp = new ArrayList<PointND>();
for(int j = 0; j < numberOfVPoints; j++){
temp.add(points.get((i*numberOfVPoints)+j));
}
vSplines[i] = new BSpline(temp, vKnots);
}
max = new PointND(-Double.MAX_VALUE, -Double.MAX_VALUE, -Double.MAX_VALUE);
min = new PointND(Double.MAX_VALUE, Double.MAX_VALUE, Double.MAX_VALUE);
for (PointND p: points) {
max.updateIfHigher(p);
min.updateIfLower(p);
}
generateBoundingPlanes();
}
protected boolean checkPointsClockwise(PointND a, PointND b, PointND center)
{
/*
* after: http://stackoverflow.com/questions/6989100/sort-points-in-clockwise-order
*/
if (a.get(0) >= 0 && b.get(0) < 0)
return true;
if (a.get(0) == 0 && b.get(0) == 0)
return a.get(1) > b.get(1);
// compute the cross product of vectors (center -> a) x (center -> b)
double det = ((a.get(0)-center.get(0)) * (b.get(1)-center.get(1)) - (b.get(0) - center.get(0)) * (a.get(1) - center.get(1)));
if (det < 0) {
return true;
} else if (det > 0) {
return false;
}
// points a and b are on the same line from the center
// check which point is closer to the center
double d1 = (a.get(0)-center.get(0)) * (a.get(0)-center.get(0)) + (a.get(1)-center.get(1)) * (a.get(1)-center.get(1));
double d2 = (b.get(0)-center.get(0)) * (b.get(0)-center.get(0)) + (b.get(1)-center.get(1)) * (b.get(1)-center.get(1));
return d1 > d2;
}
/**
* Constructor for a surface BSpline.
*
* @param controlPoints the control points
* @param uKnots the knot vector in u direction
* @param vKnots the knot vector in v direction
*/
public SurfaceBSpline(ArrayList<PointND> controlPoints, SimpleVector uKnots, SimpleVector vKnots) {
points = controlPoints;
this.uKnots = uKnots;
this.vKnots = vKnots;
init();
/* // Determination of spline orientation (increasing u,v clockwise or counterclockwise)
// assumption: splines are slice by slice
PointND center = new PointND(0,0,0);
for (PointND p: points) {
for (int d=0; d < 3; d++) center.set(d, p.get(d) / (double) points.size() + center.get(d));
}
int clockwisePoints=0;
for (int i=0; i < points.size()-1;i++) {
// actually incorrect calculation, checkPointsClockwise() only compares x,y coordinates
// works for heart splines (by chance?)
if (checkPointsClockwise(points.get(i), points.get(i+1), center))
clockwisePoints++;
else
clockwisePoints--;
}
clockwise = clockwisePoints > 0; // TODO: does not work for artery trees, nonconvex shapes
*/
// Determination of spline orientation (increasing u,v clockwise or counterclockwise)
// assumption: splines are slice by slice with increasing z
int clockwisePoints=0;
for (int i=0; i < numberOfUPoints; i++) {
PointND center = new PointND(0,0,0);
int point_index = 0;
for (int j=0; j < numberOfVPoints; j++) {
point_index = i*numberOfVPoints+j;
for (int d=0; d < 3; d++) center.set(d, points.get(point_index).get(d) / (double) numberOfVPoints + center.get(d));
}
for (int j=0; j < numberOfVPoints-1; j++) {
point_index = i*numberOfVPoints+j;
if (checkPointsClockwise(points.get(point_index), points.get(point_index+1), center))
clockwisePoints++;
else
clockwisePoints--;
}
}
clockwise = clockwisePoints > 0;
//System.out.println("Clockwisepoints " +clockwisePoints+ " for a total of "+points.size());
}
/**
* Constructor for a surface BSpline.
*
* @param list the control points
* @param uKnots the knot vector in u direction
* @param vKnots the knot vector in v direction
*/
public SurfaceBSpline(ArrayList<PointND> list, double[] uKnots,
double[] vKnots) {
this(list, new PointND(uKnots).getAbstractVector(), new PointND(vKnots).getAbstractVector());
}
/**
* Constructor for a surface BSpline.
*
* @param title the title of the BSpline
* @param list the control points
* @param uKnots the knot vector in u direction
* @param vKnots the knot vector in v direction
*/
public SurfaceBSpline(String title, ArrayList<PointND> list, double[] uKnots,
double[] vKnots) {
this(list, new PointND(uKnots).getAbstractVector(), new PointND(vKnots).getAbstractVector());
this.title = title;
}
/**
* Constructor for a surface BSpline.
*
* @param title the title of the BSpline
* @param list the control points
* @param uKnots the knot vector in u direction
* @param vKnots the knot vector in v direction
*/
public SurfaceBSpline(String title, ArrayList<PointND> list, SimpleVector uKnots,
SimpleVector vKnots) {
this(list, uKnots, vKnots);
this.title = title;
}
public static ArrayList<PointND> cloneList(List<PointND> list) {
ArrayList<PointND> clone = new ArrayList<PointND>(list.size());
for(PointND item: list) clone.add(item.clone());
return clone;
}
public SurfaceBSpline(SurfaceBSpline spline) {
this(cloneList(spline.points), new SimpleVector(spline.uKnots), new SimpleVector(spline.vKnots));
min = (spline.min!=null) ? spline.min.clone() : null;
max = (spline.max!=null) ? spline.max.clone() : null;
title = spline.title;
}
@Override
public PointND evaluate(double u, double v){
SimpleVector point = new SimpleVector(dimension);
for (int i = 0; i < numberOfUPoints; i++){
//System.out.println("Upoints: " + uPoints + " " + i);
double weight = uSpline.getWeight(u, i);
point.add(vSplines[i].evaluate(v).getAbstractVector().multipliedBy(weight));
}
PointND revan = new PointND(point);
return revan;
}
@Override
public int getDimension() {
return dimension;
}
/**
* Reads a BSpline from a BufferedReader and returns it.
* @param reader the reader
* @return the surface BSpline
* @throws IOException
*/
public static SurfaceBSpline readBSpline(BufferedReader reader) throws IOException{
boolean done = false;
int uPoints = -1;
int vPoints = -1;
double [] vKnots = null;
double [] uKnots = null;
String last = "No title";
String line = "No title";
ArrayList<PointND> list = new ArrayList<PointND>();
while (!done){
if (!line.equals("")) last = line ;
line = reader.readLine();
if (line == null){
throw new IOException("End of File reached");
}
if (line.contains(":M")){
String [] elements = line.split("\\s+");
uPoints = Integer.parseInt(elements[0]);
//System.out.println("M = " + uPoints);
break;
}
}
while (!done){
line = reader.readLine();
if (line.contains(":M")){
String [] elements = line.split("\\s+");
uPoints = Integer.parseInt(elements[0]);
//System.out.println("M = " + uPoints);
}
if (line.contains(":N")){
String [] elements = line.split("\\s+");
vPoints = Integer.parseInt(elements[0]);
//System.out.println("N = " + vPoints);
}
if (line.contains("U Knot Vector")){
vKnots = readKnotVector(reader, vPoints);
}
if (line.contains("V Knot Vector")){
uKnots = readKnotVector(reader, uPoints);
}
if (line.contains("Control Points")){
for (int i = 0; i < vPoints; i++){
list.add(new PointND(read3DPoint(reader)));
}
for (int j = 1; j < uPoints; j++) {
reader.readLine();
for (int i = 0; i < vPoints; i++){
list.add(new PointND(read3DPoint(reader)));
}
}
//System.out.println("Last Point " + list.get(list.size()-1).toString());
done = true;
break;
}
}
//System.out.println("SurfaceBSpline read: " + last);
//System.out.println("#CP = " + list.size());
//
//return new SurfaceBSpline(last, list, uKnots, vKnots);
if ((uKnots[0] == uKnots[1]) && (uKnots[2] == uKnots[3])) return new SurfaceUniformCubicBSpline(last, list, uKnots, vKnots);
else return new SurfaceBSpline(last, list, uKnots, vKnots);
}
public static ArrayList<SurfaceBSpline> readSplinesFromFile(String filename) throws IOException{
FileReader file = new FileReader(filename);
BufferedReader bf = new BufferedReader(file);
ArrayList<SurfaceBSpline> list = new ArrayList<SurfaceBSpline>();
// read all surface splines in the file.
boolean reading = true;
while (reading) {
try {
SurfaceBSpline spline = SurfaceBSpline.readBSpline(bf);
list.add(spline);
} catch (IOException e){
reading = false;
}
}
return list;
}
private static double [] readKnotVector(BufferedReader reader, int knotPoints) throws IOException{
double [] uKnots = null;
int degree = 0;
double number = 0;
while (number == 0){
number = readDouble(reader);
degree++;
}
degree -= 2;
//System.out.println("Degree: " + degree);
uKnots = new double [degree + knotPoints +1];
for (int i = 0; i < degree + 1; i++){
uKnots[i] = 0;
//System.out.println(uKnots[i]);
}
uKnots[degree + 1] = number;
//System.out.println(uKnots[degree+1]);
for (int i = degree +2; i< knotPoints +1 +degree; i++){
uKnots[i] = readDouble(reader);
//System.out.println(uKnots[i]);
}
return uKnots;
}
private static double readDouble(BufferedReader reader) throws IOException{
String read = reader.readLine();
String [] elements = read.split("\\s+");
if (elements.length == 2) {
return Double.parseDouble(elements[1]);
} else {
return Double.parseDouble(elements[0]);
}
}
private static double [] read3DPoint(BufferedReader reader) throws IOException {
double [] point = new double [3];
String read = reader.readLine();
String [] elements = read.split("\\s+");
if (elements.length == 4) {
point[0] = Double.parseDouble(elements[1]);
point[1] = Double.parseDouble(elements[2]);
point[2] = Double.parseDouble(elements[3]);
} else {
point[0] = Double.parseDouble(elements[0]);
point[1] = Double.parseDouble(elements[1]);
point[2] = Double.parseDouble(elements[2]);
}
return point;
}
/**
* Computes approximate u, v coordinates to of the closest control point to the given point.
* @param point the point
* @return the approximate (u, v) coordinates of the control point.
*/
public double [] computeInitialUV(PointND point){
int closeI = 0;
int closeJ = 0;
double minDistance = Double.MAX_VALUE;
for (int i = 0; i < numberOfUPoints; i++){
for (int j = 0; j < numberOfVPoints; j++){
double distance = vSplines[i].getControlPoint(j).euclideanDistance(point);
if (distance < minDistance){
minDistance = distance;
closeI = i;
closeJ = j;
}
}
}
int offsetU = (uSpline.getDegree() + 1) /2;
int offsetV = (vSplines[0].getDegree() + 1) /2;
return new double[] {Math.log(uSpline.getKnotVectorEntry(closeI+offsetU)), Math.log(vSplines[closeI].getKnotVectorEntry(closeJ + offsetV))};
}
@Override
public synchronized ArrayList<PointND> intersect(AbstractCurve other) {
if (other instanceof StraightLine){
StraightLine line = (StraightLine) other;
// These two calls are not thread safe. Hence, synchronization is required.
ArrayList<PointND> list = getHitsOnBoundingBox(line);
//System.out.println(list.size());
if (list.size() > 1) {
boolean simple = false;
ArrayList<PointND> intersectionPoints = new ArrayList<PointND>();
if (!simple) {
Edge connection = new Edge(list.get(0), list.get(list.size()-1));
//System.out.println(connection);
// points on one third to the end and one third from the start.
list.add(connection.evaluate(connection.getLastInternalIndex() / 3.0));
list.add(connection.evaluate((2.0 * connection.getLastInternalIndex()) / 3.0));
for (PointND p : list) {
BSplineIntersector intersection = new BSplineIntersector(line, this);
FunctionOptimizer optimizer = new FunctionOptimizer(2);
optimizer.setInitialX(computeInitialUV(p));
optimizer.optimizeFunction(intersection);
double [] uvvector = optimizer.getOptimum();
PointND hit = evaluate(Math.exp(uvvector[0]), Math.exp(uvvector[1]));
if (!intersectionPoints.contains(hit)){
intersectionPoints.add(hit);
}
}
} else {
intersectionPoints = list;
}
return intersectionPoints;
} else {
return null;
}
} else
throw new RuntimeException("Intersection between BSplineSurfaces and other AbstractCurves are not yet implemented.");
}
public PointND [] intersectDeCasteljau(StraightLine line){
BSplineIntersector intersection = new BSplineIntersector(line, this);
// TODO
throw new UnsupportedOperationException("This method is not implemented yet. See Nishita, Sederberg, Kakimoto. Ray Tracing Rational Surface Patches. Computer Graphics, 1990. 24(4):337-45.");
}
/**
* Tesselates the BSplineSurface into a mesh of Triangles. The parameters samplingU and samplingV define the number of points in each of the internal dimensions. Based on these points a triangle surface mesh is build. The resulting triangles are stored in a CompoundShape for each neighboring point pair in u direction. The u-triangle rings are then put into another CompoundShape. This configuration was optimal when ray tracing in x direction using XCAT BSplines.
* @param samplingU number of points in u direction
* @param samplingV number of points in v direction
* @return the tesselated mesh
*/
public AbstractShape tessellateMesh(double samplingU, double samplingV){
AbstractShape s = tessellateMesh(samplingU, samplingV, TESSELATE_COMPOUND_OF_COMPOUND_SHAPES);
s.setName(getTitle());
return s;
}
/**
* Tesselates the BSplineSurface into a mesh of Triangles. The parameters samplingU and samplingV define the number of points in each of the internal dimensions. Based on these points a triangle surface mesh is build. The resulting triangles are stored in a CompoundShape for each neighboring point pair in u direction. The mode determines the internal structure of the tesselated mesh.
* @param samplingU number of points in u direction
* @param samplingV number of points in v direction
* @param mode the internal representation of the mesh. Possible modes are
* <li>TESSELATE_COMPOUND_SHAPE</li>
* <li>TESSELATE_COMPOUND_OF_COMPOUND_SHAPES (default)</li>
* <li>TESSELATE_COMPOUND_OF_OCTREES</li>
* <li>TESSELATE_LINEAR_OCTREE</li>
* <li>TESSELATE_NESTED_OCTREE</li>
* @return the tesselated mesh
*/
public AbstractShape tessellateMesh(double samplingU, double samplingV, int mode){
switch (mode){
case TESSELATE_COMPOUND_SHAPE: return tessellateMeshCompoundShape(samplingU, samplingV);
case TESSELATE_COMPOUND_OF_COMPOUND_SHAPES: return tessellateMeshNestedCompoundShapes(samplingU, samplingV);
case TESSELATE_COMPOUND_OF_OCTREES: return tessellateMeshWithCompundShapesAndLinearOctrees(samplingU, samplingV);
case TESSELATE_LINEAR_OCTREE: return tessellateMeshLinearOctree(samplingU, samplingV);
case TESSELATE_NESTED_OCTREE: return tessellateMeshNestedOctree(samplingU, samplingV);
default: return tessellateMesh(samplingU, samplingV);
}
}
private AbstractShape tessellateMeshCompoundShape(double samplingU, double samplingV){
PointND [] pts = getRasterPoints(samplingU, samplingV);
CompoundShape superShape = new CompoundShape();
CompoundShape mesh = superShape;
ArrayList<PointND> lastSlice = new ArrayList<PointND>();
for(double i =1 ; i < samplingU; i++){
double lasti = i - 1;
if (lasti < 0) lasti = samplingU - 1;
for (double j = 0; j < samplingV; j++){
double lastj = j - 1;
if (lastj < 0) lastj = samplingV - 1;
if (i == samplingU-1) lastSlice.add(pts[(int)(i*samplingV+j)]);
try{
Triangle t1 = new Triangle(pts[(int)(lasti*samplingV+lastj)] ,
pts[(int)(i*samplingV+lastj)] ,
pts[(int)(lasti*samplingV+j)]);
mesh.add(t1);
} catch (Exception e){
if (!e.getLocalizedMessage().contains("direction vector")){
System.out.println(e.getLocalizedMessage() +" "+ i + " " + j + " " + lasti + " " + lastj);
}
}
try{
Triangle t2 = new Triangle(pts[(int)(i*samplingV+lastj)] ,
pts[(int)(lasti*samplingV+j)] ,
pts[(int)(i*samplingV+j)]);
mesh.add(t2);
} catch (Exception e){
if (!e.getLocalizedMessage().contains("direction vector")){
System.out.println(e.getLocalizedMessage() +" "+ i + " " + j + " " + lasti + " " + lastj);
}
}
}
}
if (lastSlice.size() > 0) {
PointND center = General.getGeometricCenter(lastSlice);
for (int i = 1; i < lastSlice.size(); i++){
try{
mesh.add(new Triangle(center, lastSlice.get(i), lastSlice.get(i-1)));
} catch (Exception e){
if (!e.getLocalizedMessage().contains("direction vector")){
System.out.println(e.getLocalizedMessage());
}
}
}
try{
mesh.add(new Triangle(center, lastSlice.get(0), lastSlice.get(lastSlice.size()-1)));
} catch (Exception e){
if (!e.getLocalizedMessage().contains("direction vector")){
System.out.println(e.getLocalizedMessage());
}
}
}
System.out.println("Triangles: " + superShape + " " + superShape.getInternalDimension() + " " + superShape.size());
return superShape;
}
private AbstractShape tessellateMeshLinearOctree(double samplingU, double samplingV){
PointND [] pts = getRasterPoints(samplingU, samplingV);
CompoundShape superShape = new LinearOctree(min, max);
CompoundShape mesh = superShape;
ArrayList<PointND> lastSlice = new ArrayList<PointND>();
for(double i =1 ; i < samplingU; i++){
double lasti = i - 1;
if (lasti < 0) lasti = samplingU - 1;
for (double j = 0; j < samplingV; j++){
double lastj = j - 1;
if (lastj < 0) lastj = samplingV - 1;
if (i == samplingU-1) lastSlice.add(pts[(int)(i*samplingV+j)]);
try{
Triangle t1 = new Triangle(pts[(int)(lasti*samplingV+lastj)] ,
pts[(int)(i*samplingV+lastj)] ,
pts[(int)(lasti*samplingV+j)]);
mesh.add(t1);
} catch (Exception e){
if (!e.getLocalizedMessage().equals("Second given direction vector must not be null!")){
System.out.println(e.getLocalizedMessage() +" "+ i + " " + j + " " + lasti + " " + lastj);
}
}
try{
Triangle t2 = new Triangle(pts[(int)(i*samplingV+lastj)] ,
pts[(int)(lasti*samplingV+j)] ,
pts[(int)(i*samplingV+j)]);
mesh.add(t2);
} catch (Exception e){
if (!e.getLocalizedMessage().equals("Second given direction vector must not be null!")){
System.out.println(e.getLocalizedMessage() +" "+ i + " " + j + " " + lasti + " " + lastj);
}
}
}
}
if (lastSlice.size() > 0) {
PointND center = General.getGeometricCenter(lastSlice);
for (int i = 1; i < lastSlice.size(); i++){
try{
mesh.add(new Triangle(center, lastSlice.get(i), lastSlice.get(i-1)));
} catch (Exception e){
if (!e.getLocalizedMessage().equals("Second given direction vector must not be null!")){
System.out.println(e.getLocalizedMessage());
}
}
}
try{
mesh.add(new Triangle(center, lastSlice.get(0), lastSlice.get(lastSlice.size()-1)));
} catch (Exception e){
if (!e.getLocalizedMessage().equals("Second given direction vector must not be null!")){
System.out.println(e.getLocalizedMessage());
}
}
}
System.out.println("Triangles: " + superShape + " " + superShape.getInternalDimension() + " " + superShape.size());
return superShape;
}
private AbstractShape tessellateMeshNestedOctree(double samplingU, double samplingV){
PointND [] pts = getRasterPoints(samplingU, samplingV);
CompoundShape superShape = new NestedOctree(min, max);
CompoundShape mesh = superShape;
ArrayList<PointND> lastSlice = new ArrayList<PointND>();
for(double i =1 ; i < samplingU; i++){
double lasti = i - 1;
if (lasti < 0) lasti = samplingU - 1;
for (double j = 0; j < samplingV; j++){
double lastj = j - 1;
if (lastj < 0) lastj = samplingV - 1;
if (i == samplingU-1) lastSlice.add(pts[(int)(i*samplingV+j)]);
try{
Triangle t1 = new Triangle(pts[(int)(lasti*samplingV+lastj)] ,
pts[(int)(i*samplingV+lastj)] ,
pts[(int)(lasti*samplingV+j)]);
mesh.add(t1);
} catch (Exception e){
if (!e.getLocalizedMessage().equals("Second given direction vector must not be null!")){
System.out.println(e.getLocalizedMessage() +" "+ i + " " + j + " " + lasti + " " + lastj);
}
}
try{
Triangle t2 = new Triangle(pts[(int)(i*samplingV+lastj)] ,
pts[(int)(lasti*samplingV+j)] ,
pts[(int)(i*samplingV+j)]);
mesh.add(t2);
} catch (Exception e){
if (!e.getLocalizedMessage().equals("Second given direction vector must not be null!")){
System.out.println(e.getLocalizedMessage() +" "+ i + " " + j + " " + lasti + " " + lastj);
}
}
}
}
if (lastSlice.size() > 0) {
PointND center = General.getGeometricCenter(lastSlice);
for (int i = 1; i < lastSlice.size(); i++){
try{
mesh.add(new Triangle(center, lastSlice.get(i), lastSlice.get(i-1)));
} catch (Exception e){
if (!e.getLocalizedMessage().equals("Second given direction vector must not be null!")){
System.out.println(e.getLocalizedMessage());
}
}
}
try{
mesh.add(new Triangle(center, lastSlice.get(0), lastSlice.get(lastSlice.size()-1)));
} catch (Exception e){
if (!e.getLocalizedMessage().equals("Second given direction vector must not be null!")){
System.out.println(e.getLocalizedMessage());
}
}
}
mesh.getMin();
System.out.println("Triangles: " + superShape + " " + superShape.getInternalDimension() + " " + superShape.size());
return superShape;
}
private AbstractShape tessellateMeshNestedCompoundShapes(double samplingU, double samplingV){
PointND [] pts = getRasterPoints(samplingU, samplingV);
CompoundShape superShape = new CompoundShape();
ArrayList<PointND> lastSlice = new ArrayList<PointND>();
ArrayList<PointND> firstSlice = new ArrayList<PointND>();
for(double i =1 ; i < samplingU; i++){
CompoundShape mesh = new CompoundShape();
double lasti = i - 1;
if (lasti < 0) lasti = samplingU - 1;
for (double j = 0; j < samplingV; j++){
double lastj = j - 1;
if (lastj < 0) lastj = samplingV - 1;
if (i == samplingU-1) lastSlice.add(pts[(int)(i*samplingV+j)]);
if (i == 1) firstSlice.add(pts[(int)((i-1)*samplingV+j)]);
try{
Triangle t1 = new Triangle(pts[(int)(lasti*samplingV+lastj)] ,
pts[(int)(i*samplingV+lastj)] ,
pts[(int)(lasti*samplingV+j)]);
mesh.add(t1);
} catch (Exception e){
if (!e.getLocalizedMessage().contains("direction vector")){
System.out.println(e.getLocalizedMessage() +" "+ i + " " + j + " " + lasti + " " + lastj);
}
}
try{
Triangle t2 = new Triangle(pts[(int)(i*samplingV+lastj)] ,
pts[(int)(lasti*samplingV+j)] ,
pts[(int)(i*samplingV+j)]);
mesh.add(t2);
} catch (Exception e){
if (!e.getLocalizedMessage().contains("direction vector")){
System.out.println(e.getLocalizedMessage() +" "+ i + " " + j + " " + lasti + " " + lastj);
}
}
}
superShape.add(mesh);
}
superShape.addAll(General.createTrianglesFromPlanarPointSet(lastSlice));
superShape.addAll(General.createTrianglesFromPlanarPointSet(firstSlice));
//System.out.println("Triangles: " + superShape + " " + superShape.getInternalDimension() + " " + superShape.size());
return superShape;
}
private AbstractShape tessellateMeshWithCompundShapesAndLinearOctrees(double samplingU, double samplingV){
PointND [] pts = getRasterPoints(samplingU, samplingV);
CompoundShape superShape = new CompoundShape();
ArrayList<PointND> lastSlice = new ArrayList<PointND>();
for(double i =1 ; i < samplingU; i++){
CompoundShape mesh = new CompoundShape();
double lasti = i - 1;
if (lasti < 0) lasti = samplingU - 1;
for (double j = 0; j < samplingV; j++){
double lastj = j - 1;
if (lastj < 0) lastj = samplingV - 1;
if (i == samplingU-1) lastSlice.add(pts[(int)(i*samplingV+j)]);
try{
Triangle t1 = new Triangle(pts[(int)(lasti*samplingV+lastj)] ,
pts[(int)(i*samplingV+lastj)] ,
pts[(int)(lasti*samplingV+j)]);
mesh.add(t1);
} catch (Exception e){
if (!e.getLocalizedMessage().equals("Second given direction vector must not be null!")){
System.out.println(e.getLocalizedMessage() +" "+ i + " " + j + " " + lasti + " " + lastj);
}
}
try{
Triangle t2 = new Triangle(pts[(int)(i*samplingV+lastj)] ,
pts[(int)(lasti*samplingV+j)] ,
pts[(int)(i*samplingV+j)]);
mesh.add(t2);
} catch (Exception e){
if (!e.getLocalizedMessage().equals("Second given direction vector must not be null!")){
System.out.println(e.getLocalizedMessage() +" "+ i + " " + j + " " + lasti + " " + lastj);
}
}
}
LinearOctree tree = new LinearOctree(mesh.getMin(), mesh.getMax());
tree.addAll(mesh);
System.out.println(tree);
superShape.add(tree);
}
if (lastSlice.size() > 0) {
CompoundShape mesh = new CompoundShape();
PointND center = General.getGeometricCenter(lastSlice);
for (int i = 1; i < lastSlice.size(); i++){
try{
mesh.add(new Triangle(center, lastSlice.get(i), lastSlice.get(i-1)));
} catch (Exception e){
if (!e.getLocalizedMessage().equals("Second given direction vector must not be null!")){
System.out.println(e.getLocalizedMessage());
}
}
}
try{
mesh.add(new Triangle(center, lastSlice.get(0), lastSlice.get(lastSlice.size()-1)));
} catch (Exception e){
if (!e.getLocalizedMessage().equals("Second given direction vector must not be null!")){
System.out.println(e.getLocalizedMessage());
}
}
LinearOctree tree = new LinearOctree(mesh.getMin(), mesh.getMax());
tree.addAll(mesh);
System.out.println(tree);
superShape.add(tree);
}
System.out.println("Triangles: " + superShape + " " + superShape.getInternalDimension() + " " + superShape.size());
return superShape;
}
@Override
public void applyTransform(Transform t) {
ArrayList<PointND> newPoints = new ArrayList<PointND>();
for (int i =0; i< points.size(); i++){
newPoints.add(t.transform(points.get(i)));
}
points = newPoints;
init();
}
@Override
public boolean isBounded() {
return true;
}
/**
* Class to intersect bsplines with a StraightLine. Uses a procedure similar to Toth if used as OptimizableFunction. Optimization after Nishita will be implemented next.
* @author akmaier
*
*/
private class BSplineIntersector implements OptimizableFunction {
double [] planeCoeff1, planeCoeff2;
StraightLine line;
SurfaceBSpline spline;
public BSplineIntersector(StraightLine line, SurfaceBSpline spline) {
this.line = line;
this.spline = spline;
SimpleVector e1 = new SimpleVector(line.getDirection().getLen());
e1.setElementValue(0, 1);
SimpleVector planeNormal1 = General.crossProduct(line.getDirection(), e1);
if (planeNormal1.normL2() < 0.1) {
e1.setElementValue(0, 0);
e1.setElementValue(1, 1);
planeNormal1 = General.crossProduct(line.getDirection(), e1);
}
planeNormal1.normalizeL2();
SimpleVector planeNormal2 = General.crossProduct(line.getDirection(),planeNormal1);
planeNormal2.normalizeL2();
SimpleVector point = line.getPoint().getAbstractVector();
planeCoeff1 = new double [4];
planeCoeff2 = new double [4];
planeNormal1.copyTo(planeCoeff1);
planeNormal2.copyTo(planeCoeff2);
planeCoeff1[3] = - SimpleOperators.multiplyInnerProd(planeNormal1, point);
planeCoeff2[3] = - SimpleOperators.multiplyInnerProd(planeNormal2, point);
}
/**
* computes d1 (after Nishita)
* @param x
* @return
*/
public double computeD1 (double [] x){
double d1 = 0;
for (int i = 0; i < numberOfUPoints; i++){
double weight = uSpline.getWeight(x[0], i);
double d1Part = 0;
for (int j = 0; j < numberOfVPoints; j++){
double weight2 = vSplines[i].getWeight(x[1], j);
d1Part += weight2 * vSplines[i].distance(planeCoeff1, j);
}
d1 += weight * d1Part;
}
return d1;
}
/**
* computes d2 (after Nishita)
* @param x
* @return
*/
public double computeD2 (double [] x){
double d2 = 0;
for (int i = 0; i < numberOfUPoints; i++){
double weight = uSpline.getWeight(x[0], i);
double d2Part = 0;
for (int j = 0; j < numberOfVPoints; j++){
double weight2 = vSplines[i].getWeight(x[1], j);
d2Part += weight2 * vSplines[i].distance(planeCoeff2, j);
}
d2 += weight * d2Part;
}
return d2;
}
public double evaluate_planebased(double[] x, int block) {
double d1 = 0;
double d2 = 0;
for (int i = 0; i < numberOfUPoints; i++){
double weight = uSpline.getWeight(Math.exp(x[0]), i);
double d1Part = 0;
double d2Part = 0;
for (int j = 0; j < numberOfVPoints; j++){
double weight2 = vSplines[i].getWeight(Math.exp(x[1]), j);
d1Part += weight2 * vSplines[i].distance(planeCoeff1, j);
d2Part += weight2 * vSplines[i].distance(planeCoeff2, j);
}
d1 += weight * d1Part;
d2 += weight * d2Part;
}
double revan = Math.pow(d1, 2) + Math.pow(d2, 2);
return revan;
}
@Override
public double evaluate(double[] x, int block) {
PointND p = spline.evaluate(Math.exp(x[0]), Math.exp(x[1]));
double revan = line.computeDistanceTo(p);
return revan;
}
@Override
public int getNumberOfProcessingBlocks() {
return 1;
}
@Override
public void setNumberOfProcessingBlocks(int number) {
}
}
@Override
public PointND[] getRasterPoints(int number) {
int subNumber = (int) Math.sqrt(number);
PointND [] array = new PointND[subNumber * subNumber];
for (double i = 0; i <subNumber ; i++){
for (double j = 0; j < subNumber; j++){
array[(int)(i+(j*subNumber))] = evaluate(i/(subNumber-1), j/(subNumber-1));
}
}
return array;
}
/**
* Binary Representation of a Surface BSpline:
* <pre>
* type
* total size in float values
* ID
* # of u knots
* # of v knots
* # of u points
* # of v points
* u knot vector
* v knot vector
* u BSpline
* v BSplines
* </pre>
* @return the binary representation for use with OpenCL
*/
public float[] getBinaryRepresentation() {
float type = BSpline.SPLINE2DTO3D;
int totalsize;
float uknot = uKnots.getLen();
float vknot = vKnots.getLen();
float uPoints = this.numberOfUPoints;
float vPoints = this.numberOfVPoints;
float [] uSpline = this.uSpline.getBinaryRepresentation();
float [][] vSplines = new float [this.vSplines.length][];
int vSplinesLength = 0;
for (int i = 0; i< vSplines.length; i++){
vSplines[i] = this.vSplines[i].getBinaryRepresentation();
vSplinesLength += vSplines[i].length;
}
totalsize = 7 + uSpline.length + vSplines.length + uKnots.getLen() + vKnots.getLen();
float [] binary = new float [totalsize];
binary[0] = type;
binary[1] = totalsize;
binary[2] = BSpline.getID(title);
binary[3] = uknot;
binary[4] = vknot;
binary[5] = uPoints;
binary[6] = vPoints;
int index =7;
for (int i=0; i<uknot; i++){
binary[index] = (float)uKnots.getElement(i);
index ++;
}
for (int i=0; i<vknot; i++){
binary[index] = (float)vKnots.getElement(i);
index ++;
}
for (int i=0; i<uSpline.length; i++){
binary[index] = uSpline[i];
index ++;
}
for (int i=0; i<vSplines.length; i++){
for (int j = 0; j< vSplines[i].length; j++){
binary[index] = vSplines[i][j];
index ++;
}
}
return binary;
}
@Override
public AbstractShape tessellate(double accuracy) {
return this.tessellateMesh(100.0 / accuracy, 100.0 / accuracy);
}
/**
* @return the numberOfUPoints
*/
public int getNumberOfUPoints() {
return numberOfUPoints;
}
/**
* @param numberOfUPoints the numberOfUPoints to set
*/
public void setNumberOfUPoints(int numberOfUPoints) {
this.numberOfUPoints = numberOfUPoints;
}
/**
* @return the numberOfVPoints
*/
public int getNumberOfVPoints() {
return numberOfVPoints;
}
/**
* @param numberOfVPoints the numberOfVPoints to set
*/
public void setNumberOfVPoints(int numberOfVPoints) {
this.numberOfVPoints = numberOfVPoints;
}
public PointND[] getRasterPoints(double samplingU, double samplingV) {
PointND [] pts = new PointND[((int)samplingU) * ((int)samplingV)];
for(double i =0 ; i < samplingU; i++){
for (double j = 0; j < samplingV; j++){
PointND p = evaluate(i/samplingU, j/ samplingV);
pts[(int)(i*samplingV+j)] = p;
}
}
return pts;
}
@Override
public AbstractShape clone() {
return new SurfaceBSpline(this);
}
}
/*
* Copyright (C) 2010-2014 Andreas Maier
* CONRAD is developed as an Open Source project under the GNU General Public License (GPL).
*/