//
// Delaunay.java
//
/*
VisAD system for interactive analysis and visualization of numerical
data. Copyright (C) 1996 - 2017 Bill Hibbard, Curtis Rueden, Tom
Rink, Dave Glowacki, Steve Emmerson, Tom Whittaker, Don Murray, and
Tommy Jasmin.
This library is free software; you can redistribute it and/or
modify it under the terms of the GNU Library General Public
License as published by the Free Software Foundation; either
version 2 of the License, or (at your option) any later version.
This library is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
Library General Public License for more details.
You should have received a copy of the GNU Library General Public
License along with this library; if not, write to the Free
Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
MA 02111-1307, USA
*/
package visad;
/**
Delaunay represents an abstract class for calculating an
N-dimensional Delaunay triangulation, that can be extended
to allow for various triangulation algorithms.
*/
public abstract class Delaunay implements java.io.Serializable {
/**
* triangles/tetrahedra --> vertices.
* Tri = new int[ntris][dim + 1]
*
*
This is the key output, a list of triangles (in two dimensions,
* tetrahedra in three dimensions, etc). ntris is the number of triangles.
*
*
In 2-D, Tri[i] is an array of 3 integers, which are three indices into
* the samples[0] and samples[1] arrays to get the x and y values of the
* three vertices of the triangle.
*
*
In 3-D, Tri[i] is an array of 4 integers, which are four indices into
* the samples[0], samples[1] and samples[2] arrays to get the x, y and z
* values of the four vertices of the tetrahedron.
*
*
This pattern continues for higher dimensionalities.
*/
public int[][] Tri;
/**
* vertices --> triangles/tetrahedra.
* Vertices = new int[nrs][nverts[i]]
*
*
nrs is the number of samples (the length of the samples[0] and
* samples[1] arrays. For sample i, Vertices[i] is a (variable length) list
* of indices into the Tri array above, giving the indices of the triangles
* that include vertex i.
*
*
nverts is an array as the second index of the Vertices array since
* different vertices may be part of different numbers of triangles.
*
*
You can use Tri and Vertices together to traverse the triangulation.
* If you don't need to traverse, then you can probably ignore all arrays
* except Tri.
*/
public int[][] Vertices;
/**
* triangles/tetrahedra --> triangles/tetrahedra.
* Walk = new int[ntris][dim + 1]
*
*
Also useful for traversing the triangulation, in this case giving the
* indices of triangles that share edges with the current triangle.
*/
public int[][] Walk;
/**
* tri/tetra edges --> global edge number.
* Edges = new int[ntris][3 * (dim - 1)];
*
*
'global edge number' is the number of an edge that is unique among the
* whole triangulation. This number is not an index into any array, but will
* match for a shared edge between two triangles.
*/
public int[][] Edges;
/** number of unique global edge numbers */
public int NumEdges;
private boolean nonConvex = false;
/**
* The abstract constructor initializes the class's data arrays.
* @throws VisADException a VisAD error occurred
*/
public Delaunay() throws VisADException {
Tri = null;
Vertices = null;
Walk = null;
Edges = null;
NumEdges = 0;
}
/**
* set flag indicating this Delaunay topology is non-convex
*/
public void setNonConvex() {
nonConvex = true;
}
/**
* @return flag indicating whether this Delaunay topology is non-convex
*/
public boolean getNonConvex() {
return nonConvex;
}
/**
* @return clone of this Delaunay as a DelaunayCustom
*/
public Object clone() {
try {
return new DelaunayCustom(null, Tri, Vertices, Walk, Edges, NumEdges);
}
catch (VisADException e) {
throw new VisADError("Delaunay.clone: " + e.toString());
}
}
/**
* The factory class method heuristically decides which extension
* to the Delaunay abstract class to use in order to construct the
* fastest triangulation, and calls that extension, returning the
* finished triangulation. The method chooses from among the Fast,
* Clarkson, and Watson methods.
* @param samples locations of points for topology - dimensioned
* float[dimension][number_of_points]
* @param exact flag indicating need for exact Delaunay triangulation
* @return a topology using an appropriate sub-class of Delaunay
* @throws VisADException a VisAD error occurred
*/
public static Delaunay factory(float[][] samples, boolean exact)
throws VisADException {
/* Note: Clarkson doesn't work well for very closely clumped site values,
since the algorithm rounds each value to the nearest integer
before computing the triangulation. This fact should probably
be taken into account in this factory algorithm, but as of yet
is not. In other words, if you need an exact triangulation
and have more than 3000 data sites, and they have closely
clumped values, be sure to scale them up before calling the
factory method. */
/* Note: The factory method will not take new Delaunay extensions into
account unless it is extended as well. */
int choice;
int FAST = 0;
int CLARKSON = 1;
int WATSON = 2;
int dim = samples.length;
if (dim < 2) throw new VisADException("Delaunay.factory: "
+"dimension must be 2 or higher");
// only Clarkson can handle triangulations in high dimensions
if (dim > 3) {
choice = CLARKSON;
}
else {
int nrs = samples[0].length;
for (int i=1; i 10000) {
// use fast in 2-D with a very large set and exact not required
choice = FAST;
}
else if (nrs > 3000) {
// use Clarkson for large sets
choice = CLARKSON;
}
else {
choice = WATSON;
}
}
try {
if (choice == FAST) {
// triangulate with the Fast method and one improvement pass
DelaunayFast delan = new DelaunayFast(samples);
delan.improve(samples, 1);
return (Delaunay) delan;
}
if (choice == CLARKSON) {
// triangulate with the Clarkson method
DelaunayClarkson delan = new DelaunayClarkson(samples);
return (Delaunay) delan;
}
if (choice == WATSON) {
// triangulate with the Watson method
DelaunayWatson delan = new DelaunayWatson(samples);
return (Delaunay) delan;
}
}
catch (Exception e) {
if (choice != CLARKSON) {
try {
// triangulate with the Clarkson method
DelaunayClarkson delan = new DelaunayClarkson(samples);
return (Delaunay) delan;
}
catch (Exception ee) {
}
}
}
return null;
}
/**
* alters the values of the samples by multiplying them by
* the mult factor
* @param samples locations of points for topology - dimensioned
* float[dimension][number_of_points]
* @param mult multiplication factor
* @param copy specifies whether scale should modify and return the
* argument samples array or a copy
* @return array of scaled values
*/
public static float[][] scale(float[][] samples, float mult,
boolean copy) {
int dim = samples.length;
int nrs = samples[0].length;
for (int i=1; i= nrs) {
if (printErrors) {
System.err.println("Delaunay.test: illegal triangle vertex (" +
"Tri[" + i + "][" + j + "]=" + Tri[i][j] + "; nrs=" + nrs + ")");
}
return false;
}
}
}
// verify that all points are in at least one triangle
int[] nverts = new int[nrs];
for (int i=0; i 2) {
throw new UnimplementedException("Delaunay.improve: dimension " +
"must be 2!");
}
int ntris = Tri.length;
int nrs = samples[0].length;
for (int i=1; i= 0) {
int[] trit2 = Tri[t2];
int[] walkt2 = Walk[t2];
int[] edgest2 = Edges[t2];
// check if the diagonal needs to be flipped
int f = (walkt2[0] == t) ? 0 :
(walkt2[1] == t) ? 1 : 2;
int A = (e + 2) % 3;
int B = (A + 1) % 3;
int C = (B + 1) % 3;
int D = (f + 2) % 3;
float ax = samp0[trit[A]];
float ay = samp1[trit[A]];
float bx = samp0[trit[B]];
float by = samp1[trit[B]];
float cx = samp0[trit[C]];
float cy = samp1[trit[C]];
float dx = samp0[trit2[D]];
float dy = samp1[trit2[D]];
float abx = ax - bx;
float aby = ay - by;
float acx = ax - cx;
float acy = ay - cy;
float dbx = dx - bx;
float dby = dy - by;
float dcx = dx - cx;
float dcy = dy - cy;
float Q = abx*acx + aby*acy;
float R = dbx*abx + dby*aby;
float S = acx*dcx + acy*dcy;
float T = dbx*dcx + dby*dcy;
boolean QD = abx*acy - aby*acx >= 0;
boolean RD = dbx*aby - dby*abx >= 0;
boolean SD = acx*dcy - acy*dcx >= 0;
boolean TD = dcx*dby - dcy*dbx >= 0;
boolean sig = (QD ? 1 : 0) + (RD ? 1 : 0)
+ (SD ? 1 : 0) + (TD ? 1 : 0) < 2;
boolean d;
if (QD == sig) d = true;
else if (RD == sig) d = false;
else if (SD == sig) d = false;
else if (TD == sig) d = true;
else if (Q < 0 && T < 0 || R > 0 && S > 0) d = true;
else if (R < 0 && S < 0 || Q > 0 && T > 0) d = false;
else if ((Q < 0 ? Q : T) < (R < 0 ? R : S)) d = true;
else d = false;
if (d) {
// diagonal needs to be swapped
eflipped = true;
int n1 = trit[A];
int n2 = trit[B];
int n3 = trit[C];
int n4 = trit2[D];
int w1 = walkt[A];
int w2 = walkt[C];
int e1 = edgest[A];
int e2 = edgest[C];
int w3, w4, e3, e4;
if (trit2[(D+1)%3] == trit[C]) {
w3 = walkt2[D];
w4 = walkt2[(D+2)%3];
e3 = edgest2[D];
e4 = edgest2[(D+2)%3];
}
else {
w3 = walkt2[(D+2)%3];
w4 = walkt2[D];
e3 = edgest2[(D+2)%3];
e4 = edgest2[D];
}
// update Tri array
trit[0] = n1;
trit[1] = n2;
trit[2] = n4;
trit2[0] = n1;
trit2[1] = n4;
trit2[2] = n3;
// update Walk array
walkt[0] = w1;
walkt[1] = w4;
walkt[2] = t2;
walkt2[0] = t;
walkt2[1] = w3;
walkt2[2] = w2;
if (w2 >= 0) {
int val = (Walk[w2][0] == t) ? 0
: (Walk[w2][1] == t) ? 1 : 2;
Walk[w2][val] = t2;
}
if (w4 >= 0) {
int val = (Walk[w4][0] == t2) ? 0
: (Walk[w4][1] == t2) ? 1 : 2;
Walk[w4][val] = t;
}
// update Edges array
edgest[0] = e1;
edgest[1] = e4;
// Edges[t][2] and Edges[t2][0] stay the same
edgest2[1] = e3;
edgest2[2] = e2;
// update Vertices array
int[] vertn1 = Vertices[n1];
int[] vertn2 = Vertices[n2];
int[] vertn3 = Vertices[n3];
int[] vertn4 = Vertices[n4];
int ln1 = vertn1.length;
int ln2 = vertn2.length;
int ln3 = vertn3.length;
int ln4 = vertn4.length;
int[] tn1 = new int[ln1 + 1]; // Vertices[n1] adds t2
int[] tn2 = new int[ln2 - 1]; // Vertices[n2] loses t2
int[] tn3 = new int[ln3 - 1]; // Vertices[n3] loses t
int[] tn4 = new int[ln4 + 1]; // Vertices[n4] adds t
System.arraycopy(vertn1, 0, tn1, 0, ln1);
tn1[ln1] = t2;
int c = 0;
for (int i=0; i= 0) {
int cside = -1;
for (int k=0; k<3; k++) {
if (Walk[othtri][k] == i) cside = k;
}
if (cside != -1) {
Edges[othtri][cside] = NumEdges;
}
else {
throw new SetException("Delaunay.finish_triang: " +
"error in triangulation!");
}
}
Edges[i][j] = NumEdges++;
}
}
}
}
else { // mdim == 3
int[] ptlook1 = {0, 0, 0, 1, 1, 2};
int[] ptlook2 = {1, 2, 3, 2, 3, 3};
for (int i=0; i= 0) {
int cside = -1;
for (int k=0; k<3; k++) {
if (Walk[othtri][k] == i) cside = k;
}
if (cside != -1) {
Edges[othtri][cside] = NumEdges;
}
else {
throw new SetException("Delaunay.finish_triang: " +
"error in triangulation!");
}
}
Edges[i][j] = NumEdges++;
}
}
}
}
else { // dim == 3
int[] ptlook1 = {0, 0, 0, 1, 1, 2};
int[] ptlook2 = {1, 2, 3, 2, 3, 3};
for (int i=0; i " + samples[0][i] + " " +
samples[1][i] + " " + samples[2][i] + "\n");
}
s.append("\n");
}
s.append("\nTri (triangles -> vertices) " + Tri.length + "\n");
for (int i=0; i ");
for (int j=0; j triangles) " + Vertices.length + "\n");
for (int i=0; i ");
for (int j=0; j triangles) " + Walk.length + "\n");
for (int i=0; i ");
for (int j=0; j global edges) " + Edges.length + "\n");
for (int i=0; i ");
for (int j=0; j
Ads help maintain this website.