//
// DelaunayOverlap.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;
// AWT is used in main method, for testing
import java.awt.*;
import java.awt.event.*;
/**
DelaunayOverlap quickly constructs a triangulation of a series
of partially overlapping two-dimensional gridded sets.
Grids should be aligned so that they overlap in a roughly
vertical column (i.e., as samples[*][1] increases or decreases,
the grid number increases or decreases respectively).
DelaunayOverlap does not handle grids which overlap across a
horizontal row.
*/
public class DelaunayOverlap extends Delaunay {
/* This method could result in overlapping triangles in the
triangulation in the intersecting section of two grids.
For "bow-tie" shaped data, this problem is unlikely to
occur, but for data with the left and right bowed inwards,
overlapping triangles are more likely. The latter type of
data is also more likely to have a greater number of
Type 3 lost points. */
// topology signature of grid boxes
private boolean sig;
/**
* Construct a Delaunay triangulation of the points in the
* samples array, which are a sequence of 2-D grids of size
* lenx * leny, and which may overlap with each other. This
* situation can arise when each grid is a scan from a polar
* orbiting satellite and the scans are perpendicular to the
* direction of travel of the satellite. The scan grids widen
* away from the satellite sub-point and can overlap.
* @param samples locations of points for topology - dimensioned
* float[dimension][number_of_points]
* @param lenx x size of scan grid
* @param leny y size of scan grid
* @throws VisADException a VisAD error occurred
*/
public DelaunayOverlap(float[][] samples, int lenx, int leny)
throws VisADException {
if (samples.length < 2) {
throw new SetException("DelaunayOverlap: bad dimension");
}
if (lenx < 2 || leny < 2) {
throw new SetException("DelaunayOverlap: must have at least "
+"two points per dimension");
}
if (samples[0].length < lenx*leny) {
throw new SetException("DelaunayOverlap: not enough samples");
}
int lenp = lenx*leny;
int leng = (lenx-1)*(leny-1);
int numgrids = samples[0].length/lenp;
if (numgrids < 2) {
throw new SetException("DelaunayOverlap: not enough grids "
+"(try Gridded2DSet instead)");
}
// samples consistency check for each grid
sig = ( (samples[0][1]-samples[0][0])
*(samples[1][lenx+1]-samples[1][1])
- (samples[1][1]-samples[1][0])
*(samples[0][lenx+1]-samples[0][1]) > 0);
for (int curgrid=0; curgrid 0 != sig)
|| ( (v11x-v10x)*(v01y-v11y)
- (v11y-v10y)*(v01x-v11x) > 0 != sig)
|| ( (v01x-v11x)*(v00y-v01y)
- (v01y-v11y)*(v00x-v01x) > 0 != sig)
|| ( (v00x-v01x)*(v10y-v00y)
- (v00y-v01y)*(v10x-v00x) > 0 != sig) ) {
throw new SetException("DelaunayOverlap: samples from grid "
+curgrid+" do not form a valid grid "
+"("+gx+","+gy+")");
}
}
}
}
// arrays for keeping track of which points fall in which grid boxes
int ngrid = (numgrids-1)*(lenx-1)*(leny-1);
int[][] grid = new int[ngrid][5];
int[] gsize = new int[ngrid];
int[] glen = new int[ngrid];
for (int g=0; g g) lfound = g;
if (gfound > g) gfound = g;
foundx[g] = gx;
foundy[g] = gy;
// mark grid box as containing current point
grid[qg][glen[qg]++] = curpt;
// expand grid array as necessary
if (glen[qg] > gsize[qg]) {
int[] newg = new int[gsize[q]+gsize[q]+1];
System.arraycopy(grid[qg], 0, newg, 0, gsize[qg]);
grid[qg] = newg;
gsize[qg] += gsize[qg];
}
// mark point as being found
ptfell[curpt] = true;
// break box to the lower left of point
if (broken[curgrid][pix] < line) {
broken[curgrid][pix] = line;
}
// break box to the lower right of point
if (broken[curgrid][pix+1] < line) {
broken[curgrid][pix+1] = line;
}
break;
case 1: case 11: // up from this box
gy--;
break;
case 2: case 7: // right from this box
gx++;
break;
case 3: // up and right from this box
gx++;
gy--;
break;
case 4: case 14: // down from this box
gy++;
break;
case 6: // down and right from this box
gx++;
gy++;
break;
case 8: case 13: // left from this box
gx--;
break;
case 9: // up and left from this box
gx--;
gy--;
break;
case 12: // down and left from this box
gx--;
gy++;
break;
default: // theoretically, should never occur
throw new SetException("DelaunayOverlap: "
+"pathological grid");
}
if (gx > lenx-2) gx = lenx-2;
if (gx < 0) gx = 0;
if (gy > leny-2) gy = leny-2;
if (gy < 0) gy = 0;
// check if point is outside the grid
if (ogx == gx && ogy == gy && !found) break;
}
// no need to check remaining grids if a grid box was located
if (found) break;
}
}
// if entire scan line wasn't in a previous grid, advance prevgrid
prevgrid = lfound;
}
// if entire grid wasn't in a previous grid, advance prevgrid
prevgrid = gfound;
}
// *** PHASE 1b: *** Break boxes that contain a point from
// the bottom line of the previous grid.
// Search logic from above is duplicated
// to maximize efficiency of the search.
// offset of the bottom scan line of a grid
int z = lenp-lenx;
// re-initialize arrays to center of each grid (old values are
// irrelevant now; center is actually a better starting guess)
for (int g=0; g lenx-2) gx = lenx-2;
if (gx < 0) gx = 0;
if (gy > leny-2) gy = leny-2;
if (gy < 0) gy = 0;
// check if point is outside the grid
if (ogx == gx && ogy == gy && !found) break;
}
}
}
// *** PHASE 2: *** triangulate points inside grid boxes
// set leftedge to the first column of grid 0
leftedge[0] = 0;
for (int i=1; i 0 && S > 0) diag = true;
else if (R < 0 && S < 0 || Q > 0 && T > 0) diag = false;
else if ((Q < 0 ? Q : T) < (R < 0 ? R : S)) diag = true;
else diag = false;
if (diag) {
tri[trilen][0] = B;
tri[trilen][1] = D;
tri[trilen][2] = A;
walk[trilen][2] = bottom[gx];
if (bottom[gx] >= 0) walk[bottom[gx]][1] = trilen;
walk[trilen][1] = trilen+1;
trilen++;
tri[trilen][0] = A;
tri[trilen][1] = D;
tri[trilen][2] = C;
if (gx > 0 && broken[curgrid][gx] < gy) {
walk[trilen][2] = tlr[qmg-1][2];
walk[tlr[qmg-1][2]][0] = trilen;
}
else walk[trilen][2] = -1;
walk[trilen][0] = trilen-1;
trilen++;
bottom[gx] = trilen-1;
tlr[qmg][0] = trilen-2;
tlr[qmg][1] = trilen-1;
tlr[qmg][2] = trilen-2;
istwo[qmg] = true;
}
else {
tri[trilen][0] = A;
tri[trilen][1] = B;
tri[trilen][2] = C;
walk[trilen][0] = bottom[gx];
if (bottom[gx] >= 0) walk[bottom[gx]][1] = trilen;
if (gx > 0 && broken[curgrid][gx] < gy) {
walk[trilen][2] = tlr[qmg-1][2];
walk[tlr[qmg-1][2]][0] = trilen;
}
else walk[trilen][2] = -1;
walk[trilen][1] = trilen+1;
trilen++;
tri[trilen][0] = B;
tri[trilen][1] = D;
tri[trilen][2] = C;
walk[trilen][2] = trilen-1;
trilen++;
bottom[gx] = trilen-1;
tlr[qmg][0] = trilen-2;
tlr[qmg][1] = trilen-2;
tlr[qmg][2] = trilen-1;
istwo[qmg] = false;
}
break;
case 1: // draw the four triangles of the point
G1 = grid[qg][0];
tri[trilen][0] = B;
tri[trilen][1] = G1;
tri[trilen][2] = A;
walk[trilen][2] = bottom[gx];
if (bottom[gx] >= 0) walk[bottom[gx]][1] = trilen;
walk[trilen][0] = trilen+1;
walk[trilen][1] = trilen+3;
trilen++;
tri[trilen][0] = B;
tri[trilen][1] = D;
tri[trilen][2] = G1;
walk[trilen][1] = trilen+1;
walk[trilen][2] = trilen-1;
trilen++;
tri[trilen][0] = G1;
tri[trilen][1] = D;
tri[trilen][2] = C;
walk[trilen][0] = trilen-1;
walk[trilen][2] = trilen+1;
trilen++;
tri[trilen][0] = A;
tri[trilen][1] = G1;
tri[trilen][2] = C;
if (gx > 0 && broken[curgrid][gx] < gy) {
walk[trilen][2] = tlr[qmg-1][2];
walk[tlr[qmg-1][2]][0] = trilen;
}
else walk[trilen][2] = -1;
walk[trilen][0] = trilen-3;
walk[trilen][1] = trilen-1;
trilen++;
bottom[gx] = trilen-2;
tlr[qmg][0] = trilen-4;
tlr[qmg][1] = trilen-1;
tlr[qmg][2] = trilen-3;
istwo[qmg] = true;
break;
default: // triangulate first two points
// use a diagonal comparison
G1 = grid[qg][0];
G2 = grid[qg][1];
float Gdx = samples[0][G2] - samples[0][G1];
float Gdy = samples[1][G2] - samples[1][G1];
ax = samples[0][A];
ay = samples[1][A];
bx = samples[0][B];
by = samples[1][B];
cx = samples[0][C];
cy = samples[1][C];
dx = samples[0][D];
dy = samples[1][D];
float val1 = Gdx*(ax-dx) + Gdy*(ay-dy);
if (val1 < 0) val1 = -val1;
float val2 = Gdx*(bx-cx) + Gdy*(by-cy);
if (val2 < 0) val2 = -val2;
if (val1 > val2) {
float Qx1 = ax - samples[0][G1];
float Qy1 = ay - samples[1][G1];
float Qx2 = ax - samples[0][G2];
float Qy2 = ay - samples[1][G2];
if (Qx1*Qx1 + Qy1*Qy1 < Qx2*Qx2 + Qy2*Qy2) {
// G1 is closer to A than G2
tri[trilen][0] = B;
tri[trilen][1] = G1;
tri[trilen][2] = A;
walk[trilen][2] = bottom[gx];
if (bottom[gx] >= 0) walk[bottom[gx]][1] = trilen;
walk[trilen][0] = trilen+5;
walk[trilen][1] = trilen+1;
trilen++;
tri[trilen][0] = A;
tri[trilen][1] = G1;
tri[trilen][2] = C;
if (gx > 0 && broken[curgrid][gx] < gy) {
walk[trilen][2] = tlr[qmg-1][2];
walk[tlr[qmg-1][2]][0] = trilen;
}
else walk[trilen][2] = -1;
walk[trilen][0] = trilen-1;
walk[trilen][1] = trilen+3;
trilen++;
tri[trilen][0] = G2;
tri[trilen][1] = D;
tri[trilen][2] = C;
walk[trilen][0] = trilen+1;
walk[trilen][2] = trilen+2;
trilen++;
tri[trilen][0] = B;
tri[trilen][1] = D;
tri[trilen][2] = G2;
walk[trilen][1] = trilen-1;
walk[trilen][2] = trilen+2;
trilen++;
tri[trilen][0] = C;
tri[trilen][1] = G1;
tri[trilen][2] = G2;
walk[trilen][0] = trilen-3;
walk[trilen][1] = trilen+1;
walk[trilen][2] = trilen-2;
trilen++;
tri[trilen][0] = B;
tri[trilen][1] = G2;
tri[trilen][2] = G1;
walk[trilen][0] = trilen-2;
walk[trilen][1] = trilen-1;
walk[trilen][2] = trilen-5;
trilen++;
bottom[gx] = trilen-4;
tlr[qmg][0] = trilen-6;
tlr[qmg][1] = trilen-5;
tlr[qmg][2] = trilen-3;
istwo[qmg] = true;
}
else {
// G2 is closer to A than G1
tri[trilen][0] = B;
tri[trilen][1] = G2;
tri[trilen][2] = A;
walk[trilen][2] = bottom[gx];
if (bottom[gx] >= 0) walk[bottom[gx]][1] = trilen;
walk[trilen][0] = trilen+5;
walk[trilen][1] = trilen+1;
trilen++;
tri[trilen][0] = A;
tri[trilen][1] = G2;
tri[trilen][2] = C;
if (gx > 0 && broken[curgrid][gx] < gy) {
walk[trilen][2] = tlr[qmg-1][2];
walk[tlr[qmg-1][2]][0] = trilen;
}
else walk[trilen][2] = -1;
walk[trilen][0] = trilen-1;
walk[trilen][1] = trilen+3;
trilen++;
tri[trilen][0] = G1;
tri[trilen][1] = D;
tri[trilen][2] = C;
walk[trilen][0] = trilen+1;
walk[trilen][2] = trilen+2;
trilen++;
tri[trilen][0] = B;
tri[trilen][1] = D;
tri[trilen][2] = G1;
walk[trilen][1] = trilen-1;
walk[trilen][2] = trilen+2;
trilen++;
tri[trilen][0] = C;
tri[trilen][1] = G2;
tri[trilen][2] = G1;
walk[trilen][0] = trilen-3;
walk[trilen][1] = trilen+1;
walk[trilen][2] = trilen-2;
trilen++;
tri[trilen][0] = B;
tri[trilen][1] = G1;
tri[trilen][2] = G2;
walk[trilen][0] = trilen-2;
walk[trilen][1] = trilen-1;
walk[trilen][2] = trilen-5;
trilen++;
bottom[gx] = trilen-4;
tlr[qmg][0] = trilen-6;
tlr[qmg][1] = trilen-5;
tlr[qmg][2] = trilen-3;
istwo[qmg] = true;
}
}
else {
float Qx1 = bx - samples[0][G1];
float Qy1 = by - samples[1][G1];
float Qx2 = bx - samples[0][G2];
float Qy2 = by - samples[1][G2];
if (Qx1*Qx1 + Qy1*Qy1 < Qx2*Qx2 + Qy2*Qy2) {
// G1 is closer to B than G2
tri[trilen][0] = B;
tri[trilen][1] = G1;
tri[trilen][2] = A;
walk[trilen][2] = bottom[gx];
if (bottom[gx] >= 0) walk[bottom[gx]][1] = trilen;
walk[trilen][0] = trilen+3;
walk[trilen][1] = trilen+4;
trilen++;
tri[trilen][0] = A;
tri[trilen][1] = G2;
tri[trilen][2] = C;
if (gx > 0 && broken[curgrid][gx] < gy) {
walk[trilen][2] = tlr[qmg-1][2];
walk[tlr[qmg-1][2]][0] = trilen;
}
else walk[trilen][2] = -1;
walk[trilen][0] = trilen+3;
walk[trilen][1] = trilen+1;
trilen++;
tri[trilen][0] = G2;
tri[trilen][1] = D;
tri[trilen][2] = C;
walk[trilen][0] = trilen+3;
walk[trilen][2] = trilen-1;
trilen++;
tri[trilen][0] = B;
tri[trilen][1] = D;
tri[trilen][2] = G1;
walk[trilen][1] = trilen+2;
walk[trilen][2] = trilen-3;
trilen++;
tri[trilen][0] = A;
tri[trilen][1] = G1;
tri[trilen][2] = G2;
walk[trilen][0] = trilen-4;
walk[trilen][1] = trilen+1;
walk[trilen][2] = trilen-3;
trilen++;
tri[trilen][0] = D;
tri[trilen][1] = G2;
tri[trilen][2] = G1;
walk[trilen][0] = trilen-3;
walk[trilen][1] = trilen-1;
walk[trilen][2] = trilen-2;
trilen++;
bottom[gx] = trilen-4;
tlr[qmg][0] = trilen-6;
tlr[qmg][1] = trilen-5;
tlr[qmg][2] = trilen-3;
istwo[qmg] = true;
}
else {
// G2 is closer to B than G1
tri[trilen][0] = B;
tri[trilen][1] = G2;
tri[trilen][2] = A;
walk[trilen][2] = bottom[gx];
if (bottom[gx] >= 0) walk[bottom[gx]][1] = trilen;
walk[trilen][0] = trilen+3;
walk[trilen][1] = trilen+4;
trilen++;
tri[trilen][0] = A;
tri[trilen][1] = G1;
tri[trilen][2] = C;
if (gx > 0 && broken[curgrid][gx] < gy) {
walk[trilen][2] = tlr[qmg-1][2];
walk[tlr[qmg-1][2]][0] = trilen;
}
else walk[trilen][2] = -1;
walk[trilen][0] = trilen+3;
walk[trilen][1] = trilen+1;
trilen++;
tri[trilen][0] = G1;
tri[trilen][1] = D;
tri[trilen][2] = C;
walk[trilen][0] = trilen+3;
walk[trilen][2] = trilen-1;
trilen++;
tri[trilen][0] = B;
tri[trilen][1] = D;
tri[trilen][2] = G2;
walk[trilen][1] = trilen+2;
walk[trilen][2] = trilen-3;
trilen++;
tri[trilen][0] = A;
tri[trilen][1] = G2;
tri[trilen][2] = G1;
walk[trilen][0] = trilen-4;
walk[trilen][1] = trilen+1;
walk[trilen][2] = trilen-3;
trilen++;
tri[trilen][0] = D;
tri[trilen][1] = G1;
tri[trilen][2] = G2;
walk[trilen][0] = trilen-3;
walk[trilen][1] = trilen-1;
walk[trilen][2] = trilen-2;
trilen++;
bottom[gx] = trilen-4;
tlr[qmg][0] = trilen-6;
tlr[qmg][1] = trilen-5;
tlr[qmg][2] = trilen-3;
istwo[qmg] = true;
}
}
// subtriangulate remaining points
int starttri = trilen-1;
int maxit = 2*glen[qg]+2;
for (int i=2; i 0 == p01x*p02y - p01y*p02x > 0 );
boolean t1 = ( c1 == 0 ||
c1 > 0 == p12x*p01y - p12y*p01x > 0 );
boolean t2 = ( c2 == 0 ||
c2 > 0 == p02x*p12y - p02y*p12x > 0 );
if (!t0 && !t1 && !t2) {
throw new SetException("DelaunayOverlap: "
+"subtriangulation error");
}
else if (!t0) {
if (curtri != tlr[qmg][2]) curtri = w0;
else throw new SetException("DelaunayOverlap: "
+"subtriangulation error");
}
else if (!t1) {
if (curtri != bottom[gx]) curtri = w1;
else throw new SetException("DelaunayOverlap: "
+"subtriangulation error");
}
else if (!t2) {
if (curtri != tlr[qmg][0] && curtri != tlr[qmg][1]) {
curtri = w2;
}
else throw new SetException("DelaunayOverlap: "
+"subtriangulation error");
}
else found = true;
}
// should never happen
if (!found) throw new SetException("DelaunayOverlap: "
+"subtriangulation error");
else if (curtri == bottom[gx]) {
// curtri is on bottom edge of grid box
// find adjoining triangles' connecting edges
int we0 = -1;
int we2 = -1;
for (int w=0; w<3; w++) {
if (walk[w0][w] == curtri) we0 = w;
if (walk[w2][w] == curtri) we2 = w;
}
if (we0 < 0 || we2 < 0) {
throw new SetException("DelaunayOverlap: "
+"subtriangulation error");
}
// define three new triangles
// the first new triangle overwrites the old triangle
tri[curtri][0] = pt;
// tri[curtri][1] = p1; <-- already true
// tri[curtri][2] = p2; <-- already true
walk[curtri][0] = trilen;
// walk[curtri][1] = w1; <-- already true
// walk[w1][we1] = curtri; <-- already true
walk[curtri][2] = trilen+1;
tri[trilen][0] = p1;
tri[trilen][1] = pt;
tri[trilen][2] = p0;
walk[trilen][0] = curtri;
walk[trilen][1] = trilen+1;
walk[trilen][2] = w0;
walk[w0][we0] = trilen;
trilen++;
tri[trilen][0] = p0;
tri[trilen][1] = pt;
tri[trilen][2] = p2;
walk[trilen][0] = trilen-1;
walk[trilen][1] = curtri;
walk[trilen][2] = w2;
walk[w2][we2] = trilen;
trilen++;
}
else if (curtri == tlr[qmg][0] || curtri == tlr[qmg][1]) {
// tlr[qmg][0]: curtri is on top edge of grid box
// tlr[qmg][1]: curtri is on left edge of grid box
// find adjoining triangles' connecting edges
int we0 = -1;
int we1 = -1;
for (int w=0; w<3; w++) {
if (walk[w0][w] == curtri) we0 = w;
if (walk[w1][w] == curtri) we1 = w;
}
if (we0 < 0 || we1 < 0) {
throw new SetException("DelaunayOverlap: "
+"subtriangulation error");
}
// define three new triangles
// the first new triangle overwrites the old triangle
// tri[curtri][0] = p0; <-- already true
tri[curtri][1] = pt;
// tri[curtri][2] = p2; <-- already true
walk[curtri][0] = trilen;
walk[curtri][1] = trilen+1;
// walk[curtri][2] = w2; <-- already true
// walk[w2][we2] = curtri; <-- already true
tri[trilen][0] = p0;
tri[trilen][1] = p1;
tri[trilen][2] = pt;
walk[trilen][0] = w0;
walk[w0][we0] = trilen;
walk[trilen][1] = trilen+1;
walk[trilen][2] = curtri;
trilen++;
tri[trilen][0] = pt;
tri[trilen][1] = p1;
tri[trilen][2] = p2;
walk[trilen][0] = trilen-1;
walk[trilen][1] = w1;
walk[w1][we1] = trilen;
walk[trilen][2] = curtri;
trilen++;
}
else {
// if curtri == tlr[qmg][2],
// then curtri is on right edge of grid box,
// otherwise, curtri is not on border of grid box
// find adjoining triangles' connecting edges
int we1 = -1;
int we2 = -1;
for (int w=0; w<3; w++) {
if (walk[w1][w] == curtri) we1 = w;
if (walk[w2][w] == curtri) we2 = w;
}
if (we1 < 0 || we2 < 0) {
throw new SetException("DelaunayOverlap: "
+"subtriangulation error");
}
// define three new triangles
// the first new triangle overwrites the old triangle
// tri[curtri][0] = p0; <-- already true
// tri[curtri][1] = p1; <-- already true
tri[curtri][2] = pt;
// walk[curtri][0] = w0; <-- already true
// walk[w0][we0] = curtri; <-- already true
walk[curtri][1] = trilen;
walk[curtri][2] = trilen+1;
tri[trilen][0] = pt;
tri[trilen][1] = p1;
tri[trilen][2] = p2;
walk[trilen][0] = curtri;
walk[trilen][1] = w1;
walk[w1][we1] = trilen;
walk[trilen][2] = trilen+1;
trilen++;
tri[trilen][0] = p0;
tri[trilen][1] = pt;
tri[trilen][2] = p2;
walk[trilen][0] = curtri;
walk[trilen][1] = trilen-1;
walk[trilen][2] = w2;
walk[w2][we2] = trilen;
trilen++;
}
}
break;
}
}
}
// *** PHASE 3: *** triangulate connection between grids
if (curgrid == 0) {
// construct ltrinum, ltriedge, rtrinum, and rtriedge
for (int i=0; i= ledges) {
throw new SetException("DelaunayOverlap: corrupt "
+"left edge structure");
}
int clep = leftedge[sg];
float cx = samples[0][clep];
float cy = samples[1][clep];
int clep1 = leftedge[sg+1];
float c1x = samples[0][clep1];
float c1y = samples[1][clep1];
// perform cross-product test to verify that points are correct
int pt = curgrid*lenp - lenx + 1;
float px = samples[0][pt];
float py = samples[1][pt];
if ( (tx-cx)*(ty-c1y) - (ty-cy)*(tx-c1x) > 0
== (px-cx)*(py-c1y) - (py-cy)*(px-c1x) > 0 ) {
// inverted overlap, must adopt slightly different strategy
start2 = ledges-1;
// use a binary search to find start1
min = 0;
max = leny-1;
sg = (min+max)/2;
itnum = 0;
float ll = samples[1][leftedge[ledges-1]];
int offst = curgrid*lenp;
int offst2 = offst+lenx;
while (samples[1][lenx*sg+offst] < ll
== samples[1][lenx*sg+offst2] < ll && itnum < leny) {
if (ll > samples[1][lenx*sg+offst]) {
// point is closer to leny-1
min = sg;
}
else {
// point is closer to 0
if (sg == 0) {
throw new SetException("DelaunayOverlap: "
+"illegal grid overlap");
}
max = sg;
}
sg = (min+max)/2;
if (sg == leny-1) {
throw new SetException("DelaunayOverlap: "
+"pathological grid overlap");
}
itnum++;
}
if (itnum >= leny) {
throw new SetException("DelaunayOverlap: corrupt "
+"grid structure");
}
start1 = lenx*sg+offst2;
}
else {
// normal overlap
start1 = tmup;
start2 = sg;
}
}
else {
// tmup is outside range of leftedge, use last leftedge pt
start1 = tmup;
start2 = ledges-1;
}
// find end1 and end2
tmup = curgrid*lenp + lenx*(broken[curgrid][lenx]+2) - 1;
tx = samples[0][tmup];
ty = samples[1][tmup];
sig = samples[1][rightedge[0]]
< samples[1][rightedge[redges-1]];
if (samples[1][rightedge[redges-1]] < ty != sig
&& samples[1][rightedge[redges-1]] != ty) {
// use a binary search to find prospective end2
int min = 0;
int max = redges-1;
int sg = (min+max)/2;
int itnum = 0;
while (samples[1][rightedge[sg+1]] < ty
== samples[1][rightedge[sg]] < ty && itnum < redges) {
if (samples[1][rightedge[sg]] < ty == sig) {
// point is closer to redges-1
min = sg;
}
else {
// point is closer to 0
if (sg == 0) {
throw new SetException("DelaunayOverlap: "
+"illegal grid overlap");
}
max = sg;
}
sg = (min+max)/2;
if (sg == redges-1) {
throw new SetException("DelaunayOverlap: "
+"pathological grid overlap");
}
itnum++;
}
if (itnum >= redges) {
throw new SetException("DelaunayOverlap: corrupt "
+"right edge structure");
}
int crep = rightedge[sg];
float cx = samples[0][crep];
float cy = samples[1][crep];
int crep1 = rightedge[sg+1];
float c1x = samples[0][crep1];
float c1y = samples[1][crep1];
// perform cross-product test to verify that points are correct
int pt = curgrid*lenp - 2;
float px = samples[0][pt];
float py = samples[1][pt];
if ( (tx-cx)*(ty-c1y) - (ty-cy)*(tx-c1x) > 0
== (px-cx)*(py-c1y) - (py-cy)*(px-c1x) > 0 ) {
// inverted overlap, must adopt slightly different strategy
end2 = redges-1;
// use a binary search to find end1
min = 0;
max = leny-1;
sg = (min+max)/2;
itnum = 0;
float ll = samples[1][rightedge[redges-1]];
int offst = curgrid*lenp+lenx-1;
int offst2 = offst+lenx;
while (samples[1][lenx*sg+offst] < ll
== samples[1][lenx*sg+offst2] < ll && itnum < leny) {
if (ll > samples[1][lenx*sg+offst]) {
// point is closer to leny-1
min = sg;
}
else {
// point is closer to 0
if (sg == 0) {
throw new SetException("DelaunayOverlap: "
+"illegal grid overlap");
}
max = sg;
}
sg = (min+max)/2;
if (sg == leny-1) {
throw new SetException("DelaunayOverlap: "
+"pathological grid overlap");
}
itnum++;
}
if (itnum >= leny) {
throw new SetException("DelaunayOverlap: corrupt "
+"grid structure");
}
end1 = lenx*sg+offst2;
}
else {
// normal overlap
end1 = tmup;
end2 = sg;
}
}
else {
// tmup is outside range of rightedge, use last rightedge pt
end1 = tmup;
end2 = redges-1;
}
// save old edge values
old_leftedge = leftedge;
old_rightedge = rightedge;
old_ledges = ledges;
old_redges = redges;
old_ltrinum = ltrinum;
old_ltriedge = ltriedge;
old_rtrinum = rtrinum;
old_rtriedge = rtriedge;
// update leftedge and rightedge arrays
leftedge = new int[numgrids*leny];
ltrinum = new int[numgrids*leny];
ltriedge = new int[numgrids*leny];
ledges = 0;
rightedge = new int[numgrids*leny];
rtrinum = new int[numgrids*leny];
rtriedge = new int[numgrids*leny];
redges = 0;
System.arraycopy(old_leftedge, 0, leftedge, 0, start2+1);
System.arraycopy(old_ltrinum, 0, ltrinum, 0, start2);
System.arraycopy(old_ltriedge, 0, ltriedge, 0, start2);
ledges += start2+1;
System.arraycopy(old_rightedge, 0, rightedge, 0, end2+1);
System.arraycopy(old_rtrinum, 0, rtrinum, 0, end2);
System.arraycopy(old_rtriedge, 0, rtriedge, 0, end2);
redges += end2+1;
for (int i=start1; i<=(curgrid+1)*lenp-lenx; i+=lenx) {
leftedge[ledges++] = i;
}
for (int i=end1; i<(curgrid+1)*lenp; i+=lenx) {
rightedge[redges++] = i;
}
int curledge = start2;
int curredge = redges - leny + broken[curgrid][lenx-1];
// indices into current quadrilateral's points
int base1, oneUp1;
int base2, oneUp2;
// x and y coordinates of base1 and oneUp1
int base1x, base1y, oneUp1x, oneUp1y;
// flag indicating which array base2 & oneUp2 reference
// b2lbr = 0 --> base2 & oneUp2 reference old_leftedge array
// b2lbr = 1 --> base2 & oneUp2 reference samples array
// (bottom of curgrid-1)
// b2lbr = 2 --> base2 & oneUp2 reference old_rightedge array
int b2lbr = 0;
// initialize base1 and base2
base1 = start1;
base2 = start2;
if (base2 == old_ledges-1) {
base2 = curgrid*lenp - lenx;
b2lbr = 1;
}
base1x = 0;
base1y = base1 % lenp / lenx;
boolean down = false;
// initialize oneUp1 and oneUp2
if (base1y > 0 && broken[curgrid][0] < base1y - 1) {
// move up
oneUp1x = 0;
oneUp1y = base1y - 1;
}
else if (base1y == leny-1 || broken[curgrid][1] < base1y) {
// move right
oneUp1x = 1;
oneUp1y = base1y;
}
else {
// move down
oneUp1x = 0;
oneUp1y = base1y + 1;
down = true;
}
oneUp1 = curgrid*lenp + oneUp1y*lenx + oneUp1x;
oneUp2 = base2 + 1;
boolean diag = false;
boolean firsttri = true;
while (base1 != end1 || base2 != end2 || b2lbr != 2) {
// determine which diagonal to take
if (base1 == end1) diag = true;
else if (base2 == end2 && b2lbr == 2) diag = false;
else {
// set up variables
float ax = samples[0][base1];
float ay = samples[1][base1];
float cx = samples[0][oneUp1];
float cy = samples[1][oneUp1];
float bx, by, dx, dy;
if (b2lbr == 0) {
bx = samples[0][old_leftedge[base2]];
by = samples[1][old_leftedge[base2]];
dx = samples[0][old_leftedge[oneUp2]];
dy = samples[1][old_leftedge[oneUp2]];
}
else if (b2lbr == 1) {
bx = samples[0][base2];
by = samples[1][base2];
dx = samples[0][oneUp2];
dy = samples[1][oneUp2];
}
else {
bx = samples[0][old_rightedge[base2]];
by = samples[1][old_rightedge[base2]];
dx = samples[0][old_rightedge[oneUp2]];
dy = samples[1][old_rightedge[oneUp2]];
}
// perform diagonal test
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 sigD = (QD ? 1 : 0) + (RD ? 1 : 0)
+ (SD ? 1 : 0) + (TD ? 1 : 0) < 2;
if (QD == sigD) diag = true;
else if (RD == sigD || SD == sigD) diag = false;
else if (TD == sigD) diag = true;
else if (Q < 0 && T < 0 || R > 0 && S > 0) diag = true;
else if (R < 0 && S < 0 || Q > 0 && T > 0) diag = false;
else if ((Q < 0 ? Q : T) < (R < 0 ? R : S)) diag = true;
else diag = false;
}
// walk[*][0] --> NEXT
// walk[*][diag ? 1 : 2] --> PREVIOUS
// walk[*][diag ? 2 : 1] --> OUTSIDE
if (diag) { // diagonal goes from base1 to oneUp2
// update tri array
if (b2lbr == 0) {
tri[trilen][0] = old_leftedge[oneUp2];
tri[trilen][1] = base1;
tri[trilen][2] = old_leftedge[base2];
}
else if (b2lbr == 1) {
tri[trilen][0] = oneUp2;
tri[trilen][1] = base1;
tri[trilen][2] = base2;
}
else {
tri[trilen][0] = old_rightedge[oneUp2];
tri[trilen][1] = base1;
tri[trilen][2] = old_rightedge[base2];
}
// update walk array
// current --> next
walk[trilen][0] = -1;
// current <--> previous
if (firsttri) {
ltrinum[curledge] = trilen;
ltriedge[curledge] = 2;
curledge++;
firsttri = false;
}
else {
walk[trilen][1] = trilen-1;
walk[trilen-1][0] = trilen;
}
// current <--> outside
if (b2lbr == 0) {
int x = old_ltrinum[base2];
walk[trilen][2] = x;
walk[x][old_ltriedge[base2]] = trilen;
}
else if (b2lbr == 1) {
int x = old_bottom[base2 % lenx];
walk[trilen][2] = x;
walk[x][1] = trilen;
}
else {
int x = old_rtrinum[oneUp2];
walk[trilen][2] = x;
walk[x][old_rtriedge[oneUp2]] = trilen;
}
// update base2
base2 = oneUp2;
if (b2lbr == 0 && base2 == old_ledges-1) {
b2lbr = 1;
base2 = curgrid*lenp - lenx;
}
if (b2lbr == 1 && base2 == curgrid*lenp - 1) {
b2lbr = 2;
base2 = old_redges-1;
}
// update oneUp2
oneUp2 = base2 + (b2lbr == 2 ? -1 : 1);
}
else { // diagonal goes from base2 to oneUp1
// update tri array
if (b2lbr == 0) {
tri[trilen][0] = old_leftedge[base2];
tri[trilen][1] = oneUp1;
tri[trilen][2] = base1;
}
else if (b2lbr == 1) {
tri[trilen][0] = base2;
tri[trilen][1] = oneUp1;
tri[trilen][2] = base1;
}
else {
tri[trilen][0] = old_rightedge[base2];
tri[trilen][1] = oneUp1;
tri[trilen][2] = base1;
}
// update walk array
// current --> next
walk[trilen][0] = -1;
// current <--> previous
if (firsttri) {
ltrinum[curledge] = trilen;
ltriedge[curledge] = 1;
curledge++;
firsttri = false;
}
else {
walk[trilen][2] = trilen-1;
walk[trilen-1][0] = trilen;
}
// current <--> outside
if (oneUp1 - base1 == -lenx) {
// base1 -> oneUp1 = up
if (base1x < lenx-1) {
int x = tlr[(lenx-1)*oneUp1y+oneUp1x][1];
walk[trilen][1] = x;
walk[x][2] = trilen;
}
else {
// base1 is on the right edge of current grid
walk[trilen][1] = -1;
// update right edge triangle structure
rtrinum[curredge] = trilen;
rtriedge[curredge] = 1;
curredge--;
}
}
else if (oneUp1 - base1 == 1) {
// base1 -> oneUp1 = right
if (base1y < leny-1) {
int inx = (lenx-1)*base1y+base1x;
int x = tlr[inx][0];
walk[trilen][1] = x;
walk[x][istwo[inx] ? 2 : 0] = trilen;
}
else {
// base1 is on the bottom edge of current grid
walk[trilen][1] = -1;
bottom[base1x] = trilen;
}
}
else {
// base1 -> oneUp1 = down
if (base1x > 0) {
int x = tlr[(lenx-1)*base1y+base1x-1][2];
walk[trilen][1] = x;
walk[x][0] = trilen;
}
else {
// base1 is on the left edge of current grid
walk[trilen][1] = -1;
ltrinum[curledge] = trilen;
ltriedge[curledge] = 1;
curledge++;
}
}
// update base1
base1 = oneUp1;
base1x = oneUp1x;
base1y = oneUp1y;
// update oneUp1
if (broken[curgrid][base1x == 0 ? 0 : base1x+1] < base1y - 1
&& base1y > 0 && !down) {
// move up
oneUp1x = base1x;
oneUp1y = base1y - 1;
}
else if ( base1y == leny-1 || (base1x < lenx-1
&& broken[curgrid][base1x+1] < base1y) ) {
// move right
oneUp1x = base1x+1;
oneUp1y = base1y;
down = false;
}
else {
// move down
oneUp1x = base1x;
oneUp1y = base1y + 1;
down = true;
}
oneUp1 = curgrid*lenp + oneUp1y*lenx + oneUp1x;
}
trilen++;
// expand tri array if necessary
if (trilen > trisize) {
trisize += trisize;
int[][] newtri = new int[trisize+2][3];
int[][] newwalk = new int[trisize+2][3];
for (int i=0; i 0) == (
(Bx-Ax)*(Cy-Ay) - (By-Ay)*(Cx-Ax) > 0) );
boolean test1 = (tval1 == 0) || ( (tval1 > 0) == (
(Cx-Bx)*(Ay-By) - (Cy-By)*(Ax-Bx) > 0) );
boolean test2 = (tval2 == 0) || ( (tval2 > 0) == (
(Ax-Cx)*(By-Cy) - (Ay-Cy)*(Bx-Cx) > 0) );
// figure out which triangle to go to next
if (!test0 && !test1 && !test2) {
throw new SetException("DelaunayOverlap: corrupt "
+"triangle structure");
}
if (!test0 && !test1) {
int tri0 = walk[curtri][0];
int tri1 = walk[curtri][1];
if (tri0 == -1) curtri = tri1;
else curtri = tri0;
}
else if (!test0 && !test2) {
int tri0 = walk[curtri][0];
int tri2 = walk[curtri][2];
if (tri0 == -1) curtri = tri2;
else curtri = tri0;
}
else if (!test1 && !test2) {
int tri1 = walk[curtri][1];
int tri2 = walk[curtri][2];
if (tri1 == -1) curtri = tri2;
else curtri = tri1;
}
else if (!test0) curtri = walk[curtri][0];
else if (!test1) curtri = walk[curtri][1];
else if (!test2) curtri = walk[curtri][2];
else located = true;
// point is outside current triangulation
if (curtri < 0) itnum = trilen;
}
// if itnum == trilen, the point is permanently lost
if (itnum < trilen) {
// subtriangulate point grid[qg][pt]
// (curtri is the containing triangle)
// get curtri's point information
int q = grid[qg][pt];
int ct0 = tri[curtri][0];
int ct1 = tri[curtri][1];
int ct2 = tri[curtri][2];
// get curtri's walk information
int T0 = walk[curtri][0];
int T1 = walk[curtri][1];
int T2 = walk[curtri][2];
int TE0, TE1, TE2;
if (T0 == -1) TE0 = -1;
else if (walk[T0][0] == curtri) TE0 = 0;
else if (walk[T0][1] == curtri) TE0 = 1;
else if (walk[T0][2] == curtri) TE0 = 2;
else throw new SetException("DelaunayOverlap: corrupt "
+"walk structure");
if (T1 == -1) TE1 = -1;
else if (walk[T1][0] == curtri) TE1 = 0;
else if (walk[T1][1] == curtri) TE1 = 1;
else if (walk[T1][2] == curtri) TE1 = 2;
else throw new SetException("DelaunayOverlap: corrupt "
+"walk structure");
if (T2 == -1) TE2 = -1;
else if (walk[T2][0] == curtri) TE2 = 0;
else if (walk[T2][1] == curtri) TE2 = 1;
else if (walk[T2][2] == curtri) TE2 = 2;
else throw new SetException("DelaunayOverlap: corrupt "
+"walk structure");
// define three new triangles
// the first new triangle overwrites the old triangle
// tri[curtri][0] = ct0; <-- already true
// tri[curtri][1] = ct1; <-- already true
tri[curtri][2] = q;
// walk[curtri][0] = T0; <-- already true
// if (TE0 >= 0) walk[T0][TE0] = curtri; <-- already true
walk[curtri][1] = trilen;
walk[curtri][2] = trilen+1;
tri[trilen][0] = q;
tri[trilen][1] = ct1;
tri[trilen][2] = ct2;
walk[trilen][0] = curtri;
walk[trilen][1] = T1;
if (TE1 >= 0) walk[T1][TE1] = trilen;
walk[trilen][2] = trilen+1;
trilen++;
tri[trilen][0] = ct0;
tri[trilen][1] = q;
tri[trilen][2] = ct2;
walk[trilen][0] = curtri;
walk[trilen][1] = trilen-1;
walk[trilen][2] = T2;
if (TE2 >= 0) walk[T2][TE2] = trilen;
trilen++;
// expand tri array if necessary
if (trilen > trisize) {
trisize += trisize;
int[][] newtri = new int[trisize+2][3];
int[][] newwalk = new int[trisize+2][3];
for (int i=0; i= line
&& broken[curgrid][pix+1] >= line) {
// point (px, py) is lost
float Px = samples[0][q];
float Py = samples[1][q];
// locate containing triangle of (px, py)
boolean located = false;
int itnum;
for (itnum=0; itnum 0) == (
(Bx-Ax)*(Cy-Ay) - (By-Ay)*(Cx-Ax) > 0) );
boolean test1 = (tval1 == 0) || ( (tval1 > 0) == (
(Cx-Bx)*(Ay-By) - (Cy-By)*(Ax-Bx) > 0) );
boolean test2 = (tval2 == 0) || ( (tval2 > 0) == (
(Ax-Cx)*(By-Cy) - (Ay-Cy)*(Bx-Cx) > 0) );
// figure out which triangle to go to next
if (!test0 && !test1 && !test2) {
throw new SetException("DelaunayOverlap: corrupt "
+"triangle structure");
}
if (!test0 && !test1) {
int tri0 = walk[curtri][0];
int tri1 = walk[curtri][1];
if (tri0 == -1) curtri = tri1;
else curtri = tri0;
}
else if (!test0 && !test2) {
int tri0 = walk[curtri][0];
int tri2 = walk[curtri][2];
if (tri0 == -1) curtri = tri2;
else curtri = tri0;
}
else if (!test1 && !test2) {
int tri1 = walk[curtri][1];
int tri2 = walk[curtri][2];
if (tri1 == -1) curtri = tri2;
else curtri = tri1;
}
else if (!test0) curtri = walk[curtri][0];
else if (!test1) curtri = walk[curtri][1];
else if (!test2) curtri = walk[curtri][2];
else located = true;
// point is outside current triangulation
// (this point is permanently lost)
if (curtri < 0) itnum = trilen;
}
// if itnum == trilen, the point is permanently lost
if (itnum < trilen) {
// subtriangulate point q
// (curtri is the containing triangle)
// get curtri's point information
int ct0 = tri[curtri][0];
int ct1 = tri[curtri][1];
int ct2 = tri[curtri][2];
// get curtri's walk information
int T0 = walk[curtri][0];
int T1 = walk[curtri][1];
int T2 = walk[curtri][2];
int TE0, TE1, TE2;
if (T0 == -1) TE0 = -1;
else if (walk[T0][0] == curtri) TE0 = 0;
else if (walk[T0][1] == curtri) TE0 = 1;
else if (walk[T0][2] == curtri) TE0 = 2;
else throw new SetException("DelaunayOverlap: corrupt "
+"walk structure");
if (T1 == -1) TE1 = -1;
else if (walk[T1][0] == curtri) TE1 = 0;
else if (walk[T1][1] == curtri) TE1 = 1;
else if (walk[T1][2] == curtri) TE1 = 2;
else throw new SetException("DelaunayOverlap: corrupt "
+"walk structure");
if (T2 == -1) TE2 = -1;
else if (walk[T2][0] == curtri) TE2 = 0;
else if (walk[T2][1] == curtri) TE2 = 1;
else if (walk[T2][2] == curtri) TE2 = 2;
else throw new SetException("DelaunayOverlap: corrupt "
+"walk structure");
// define three new triangles
// the first new triangle overwrites the old triangle
// tri[curtri][0] = ct0; <-- already true
// tri[curtri][1] = ct1; <-- already true
tri[curtri][2] = q;
// walk[curtri][0] = T0; <-- already true
// if (TE0 >= 0) walk[T0][TE0] = curtri; <-- already true
walk[curtri][1] = trilen;
walk[curtri][2] = trilen+1;
tri[trilen][0] = q;
tri[trilen][1] = ct1;
tri[trilen][2] = ct2;
walk[trilen][0] = curtri;
walk[trilen][1] = T1;
if (TE1 >= 0) walk[T1][TE1] = trilen;
walk[trilen][2] = trilen+1;
trilen++;
tri[trilen][0] = ct0;
tri[trilen][1] = q;
tri[trilen][2] = ct2;
walk[trilen][0] = curtri;
walk[trilen][1] = trilen-1;
walk[trilen][2] = T2;
if (TE2 >= 0) walk[T2][TE2] = trilen;
trilen++;
// expand tri array if necessary
if (trilen > trisize) {
trisize += trisize;
int[][] newtri = new int[trisize+2][3];
int[][] newwalk = new int[trisize+2][3];
for (int i=0; i
Ads help maintain this website.