package algorithmrepository; import java.io.FileOutputStream; import java.io.IOException; import java.io.PrintStream; import java.util.ArrayList; import java.util.Collections; import java.util.Comparator; import java.util.Random; import oneLinersForAlgoTests.algoAsciiMatrixFile; import oneLinersForAlgoTests.algoBinaryMatrixFile; import oneLinersForAlgoTests.algoOneLiners; /** * Creates a Delaunay triangulation of set of points given in the constructor and * provides linear interpolation from those fixed triangles. The function values can be modified * on-the-fly with setF(). * * Core originally ported from p bourke's triangulate.c * http://astronomy.swin.edu.au/~pbourke/terrain/triangulate/triangulate.c * * The core code here is numerically glitchy - on large numbers of systematically * varying points I have seen it generate extra triangles that cross several others. * * @author oliford, Jakob Svensson */ public class DelaunayTriangulationInterpolation { /** points are provided in constructor */ private double points[][]; private int triangles[][]; private double triangleCoords[][][]; private double f[]; int nPoints; private double outsideValue = Double.NaN; PolygonGridLookup triangleLookup; /** Build Delaunay trianglulation interpolator on given points, using given triangles * function values must be set later by setF() at points * will return NaN when evaluated outside any triangle * @param points points[index][x/y] */ public DelaunayTriangulationInterpolation(int triangles[][], double points[][]) { this(triangles, points, null, Double.NaN); } /** Build Delaunay trianglulation interpolator on given points, using given triangles * @param points points[index][x/y] * @param f function values at points * @param outsideValue value when evaluated outside any triangle */ public DelaunayTriangulationInterpolation(int triangles[][], double points[][], double f[]) { this(triangles, points, f, Double.NaN); } /** Build Delaunay trianglulation interpolator on given points, using given triangles * @param points points[index][x/y] * @param f function values at points * @param outsideValue value when evaluated outside any triangle */ public DelaunayTriangulationInterpolation(int triangles[][], double points[][], double f[], double outsideValue){ this.points = points; this.f = f; this.outsideValue = outsideValue; this.triangles = triangles; createTriangleCoords(); triangleLookup = new PolygonGridLookup(triangleCoords, 100, 100); } /** Build Delaunay trianglulation interpolator on given points, using internal Pure Java routine * function values must be set later by setF() at points * will return NaN when evaluated outside any triangle * @param points points[index][x/y] */ public DelaunayTriangulationInterpolation(double points[][]) { this(points, null, Double.NaN); } /** Build Delaunay trianglulation interpolator on given points, using internal Pure Java routine * function values must be set later by setF() at points * will return NaN when evaluated outside any triangle * @param points points[index][x/y] */ public DelaunayTriangulationInterpolation(double points[][], double f[]) { this(points, f, Double.NaN); } /** Build Delaunay trianglulation interpolator on given points, using internal Pure Java routine * @param points points[index][x/y] * @param f function values at points * @param outsideValue value when evaluated outside any triangle */ public DelaunayTriangulationInterpolation(double points[][], double f[], double outsideValue) { this.f = f; this.outsideValue = outsideValue; this.points = points; nPoints = points.length; //our internal points need to be 3 elements bigger for storage for supertriangle verticies double sortedPoints[][] = new double[nPoints+3][]; System.arraycopy(points, 0, sortedPoints, 0, nPoints); sortedPoints[nPoints] = new double[]{ Double.NaN, Double.NaN }; sortedPoints[nPoints+1] = new double[]{ Double.NaN, Double.NaN }; sortedPoints[nPoints+2] = new double[]{ Double.NaN, Double.NaN }; //and needs to be in order int unsortedPointIndex[] = sortPoints(sortedPoints); int trianglesOfSortedIndices[][] = createTriangles(sortedPoints); //now create the list of triangles, pointing at the original point indices triangles = new int[trianglesOfSortedIndices.length][3]; for(int i=0; i < trianglesOfSortedIndices.length; i++){ triangles[i][0] = unsortedPointIndex[trianglesOfSortedIndices[i][0]]; triangles[i][1] = unsortedPointIndex[trianglesOfSortedIndices[i][1]]; triangles[i][2] = unsortedPointIndex[trianglesOfSortedIndices[i][2]]; } createTriangleCoords(); triangleLookup = new PolygonGridLookup(triangleCoords, 100, 100); } /** Sets the function value at each of the points given in the constructor */ public void setF(double[] f) { this.f = f; } /** Performs a linear interpolation of the relevant triangle at a given point */ public double[] eval(double x[], double y[]){ double ret[] = new double[x.length]; for(int i=0; i < x.length; i++){ int j = triangleLookup.getPolygonAtPoint(x[i], y[i]); //int j = triangleLookup.getPolygonAtPointBruteForce(x[i], y[i]); if(j < 0) { ret[i] = outsideValue; continue; } double tx[] = new double[]{ points[ triangles[j][0] ][0], points[ triangles[j][1] ][0], points[ triangles[j][2] ][0] }; double ty[] = new double[]{ points[ triangles[j][0] ][1], points[ triangles[j][1] ][1], points[ triangles[j][2] ][1] }; double tz[] = new double[]{ f[ triangles[j][0] ], f[ triangles[j][1] ], f[ triangles[j][2] ] }; ret[i] = triangleInterp(tx, ty, tz, x[i], y[i]); } return ret; } /** Returns the triangles as indices into the points array * @return triangles[triangleIndex][i] - indices into points array of i'th index of triangle (i=0,1,2) */ public int[][] getTriangleVertexIndices(){ return triangles; } /** Returns the coordinates on the triangles * @return triangles[triangleIndex][x, y] */ public double[][] getTriangleVertexCoords(){ double vertices[][] = new double[triangles.length][6]; for (int i=0; i < triangles.length; i++) { vertices[i][0] = points[ triangles[i][0] ][0]; vertices[i][1] = points[ triangles[i][0] ][1]; vertices[i][2] = points[ triangles[i][1] ][0]; vertices[i][3] = points[ triangles[i][1] ][1]; vertices[i][4] = points[ triangles[i][2] ][0]; vertices[i][5] = points[ triangles[i][2] ][1]; } return vertices; } public double[][] getPoints(){ return points; } public PolygonGridLookup getPolygonLookup(){ return triangleLookup; } private void createTriangleCoords(){ triangleCoords = new double[triangles.length][3][]; for(int i=0; i < triangles.length; i++){ for(int j=0; j < 3; j++) triangleCoords[i][j] = points[ triangles[i][j] ]; } } /** *************** STATIC SUPPORT MEMBERS - (the guts) ******************* */ public static final double numericalTolerance = 1e-10; /** * Returns the circle that circumscribes a given triangle. * * @param x The x-coordinates of the triangle points * @param y The y-coordinates of the triangle points * @return The centre and radius of the circumcircle. * ret[0] is the centre x-coordinate, ret[1] centre y-coordinate, * ret[1] is the square radius. */ public static final double[] circumCircle(double[] x, double[] y) { double x1=x[0],x2=x[1],x3=x[2]; double y1=y[0],y2=y[1],y3=y[2]; double m1,m2,mx1,mx2,my1,my2; double dx,dy,rsqr; double xc, yc, r; // Check for coincident points if ( Math.abs(y1-y2) < numericalTolerance && Math.abs(y2-y3) < numericalTolerance ) { //throw new RuntimeException("CircumCircle: Points are coincident."); System.err.println("ERROR: CircumCircle: Points are coincident."); } if ( Math.abs(y2-y1) < numericalTolerance ) { m2 = - (x3-x2) / (y3-y2); mx2 = (x2 + x3) / 2.0; my2 = (y2 + y3) / 2.0; xc = (x2 + x1) / 2.0; yc = m2 * (xc - mx2) + my2; } else if ( Math.abs(y3-y2) < numericalTolerance ) { m1 = - (x2-x1) / (y2-y1); mx1 = (x1 + x2) / 2.0; my1 = (y1 + y2) / 2.0; xc = (x3 + x2) / 2.0; yc = m1 * (xc - mx1) + my1; } else { m1 = - (x2-x1) / (y2-y1); m2 = - (x3-x2) / (y3-y2); mx1 = (x1 + x2) / 2.0; mx2 = (x2 + x3) / 2.0; my1 = (y1 + y2) / 2.0; my2 = (y2 + y3) / 2.0; xc = (m1 * mx1 - m2 * mx2 + my2 - my1) / (m1 - m2); yc = m1 * (xc - mx1) + my1; } dx = x2 - xc; dy = y2 - yc; rsqr = dx*dx + dy*dy; return new double[] { xc, yc, rsqr }; } /** Given points, create Delaunay triganulation. * * Originally: * --------------------------------------------------------------------------------- * ported from p bourke's triangulate.c * http://astronomy.swin.edu.au/~pbourke/terrain/triangulate/triangulate.c * * fjenett, 20th february 2005, offenbach-germany. * contact: http://www.florianjenett.de/ * --------------------------------------------------------------------------------- * * @param points Array of points to create triangles from. The array should be 3 elements bigger * than required as the last three elements are used for the temporary super-triangle. * Points must be in order in one coordinate direction. It works if it is increasing in points[][0] * but may also work for other orders. * @return triangles[triangleIndex][i] - indices into points array of i'th index of triangle (i=0,1,2) */ public static final int[][] createTriangles(double points[][]) { boolean complete[] = null; int edges[][] = null; int nEdges = 0; int maxTriangles; int maxEdges = 200; int nVertices = points.length - 3; int triangles[][] = null; double xp, yp, x1, y1, x2, y2, x3, y3, xc, yc, r; double xmin, xmax, ymin, ymax, xmid, ymid; double dx, dy, dmax; int nTriangles = 0; /* Allocate memory for the completeness list, flag for each triangle */ maxTriangles = 4*nVertices; complete = new boolean[maxTriangles]; triangles = new int[maxTriangles][3]; /* Allocate edge list */ edges = new int[maxEdges][2]; for (int ie=0; ie xmax) xmax = points[i][0]; if (points[i][1] < ymin) ymin = points[i][1]; if (points[i][1] > ymax) ymax = points[i][1]; } dx = xmax - xmin; dy = ymax - ymin; dmax = (dx > dy) ? dx : dy; xmid = (xmax + xmin) / 2.0; ymid = (ymax + ymin) / 2.0; /* Set up the supertriangle This is a triangle which encompasses all the sample points. The supertriangle coordinates are added to the end of the vertex list. The supertriangle is the first triangle in the triangle list. */ points[nVertices+0][0] = xmid - 2.0 * dmax; points[nVertices+0][1] = ymid - dmax; points[nVertices+1][0] = xmid; points[nVertices+1][1] = ymid + 2.0 * dmax; points[nVertices+2][0] = xmid + 2.0 * dmax; points[nVertices+2][1] = ymid - dmax; triangles[0][0] = nVertices; triangles[0][1] = nVertices+1; triangles[0][2] = nVertices+2; complete[0] = false; nTriangles = 1; /* Include each point one at a time into the existing mesh */ for (int i=0; i < nVertices; i++) { xp = points[i][0]; yp = points[i][1]; nEdges = 0; /* Set up the edge buffer. If the point (xp,yp) lies inside the circumcircle then the three edges of that triangle are added to the edge buffer and that triangle is removed. */ for (int j=0; j < nTriangles; j++) { if (complete[j]) continue; x1 = points[ triangles[j][0] ][0]; y1 = points[ triangles[j][0] ][1]; x2 = points[ triangles[j][1] ][0]; y2 = points[ triangles[j][1] ][1]; x3 = points[ triangles[j][2] ][0]; y3 = points[ triangles[j][2] ][1]; double circ[] = circumCircle(new double[]{ x1, x2, x3 }, new double[]{ y1, y2, y3}); double pointRadius2 = (xp - circ[0])*(xp - circ[0]) + (yp - circ[1])*(yp - circ[1]); xc = circ[0]; yc = circ[1]; r = Math.sqrt(circ[2]); if (xc + r < xp) complete[j] = true; if (pointRadius2 <= circ[2]) { /* Check that we haven't exceeded the edge list size */ if (nEdges+3 >= maxEdges) { maxEdges += 100; int[][] edges_n = new int[maxEdges][]; System.arraycopy(edges, 0, edges_n, 0, edges.length); for (int ie=edges.length; ie= maxTriangles) throw new RuntimeException("Run out of space for triangles. Not sure if this should be possible."); triangles[nTriangles][0] = edges[j][0]; triangles[nTriangles][1] = edges[j][1]; triangles[nTriangles][2] = i; complete[nTriangles] = false; nTriangles++; } } /* Remove triangles with supertriangle vertices These are triangles which have a vertex number greater than nv */ for (int i=0;i= nVertices || triangles[i][1] >= nVertices || triangles[i][2] >= nVertices) { triangles[i] = triangles[nTriangles-1]; nTriangles--; i--; } } //create triangles list only as long as the ones we have int ret[][] = new int[nTriangles][]; System.arraycopy(triangles, 0, ret, 0, nTriangles); return ret; } /** Sorts the points in points[idx][x, y] into ascending order of x * @return the indices into the original array, of each element of the new one * */ public static final int[] sortPoints(double points[][]){ class CompareByIndex implements Comparator { double points[][]; public CompareByIndex(double points[][]) { this.points = points; } @Override public int compare(Integer idx1, Integer idx2){ return Double.compare(points[idx1][0], points[idx2][0]); } } //create list of indices ArrayList pointIdxList = new ArrayList(points.length); for (int i=0; i -numericalTolerance) && (v > -numericalTolerance) && (u + v < (1 + numericalTolerance)); } /** Curtesy of biki_, * http://www.gamedev.net/community/forums/topic.asp?topic_id=380061 * @return Value of z at (px, py) on the plane that passes through * the three points in (x[0-2], y[0-2], z[0-2]) * */ private static final double triangleInterp(double x[], double y[], double z[], double px, double py){ double P = x[0]-x[1], Q = y[0]-y[1]; double R = x[0]-x[2], S = y[0]-y[2]; double T = z[0]-z[1], U = z[0]-z[2]; double det = P*S - R*Q; double A = (T*S - U*Q) / det; double B = (P*U - R*T) / det; double C = z[0] - A*x[0] - B*y[0]; return A*px + B*py + C; } }