root / trunk / extensions / extGraph / src / org / gvsig / fmap / algorithm / triangulation / visad / DelaunayOverlap.java @ 39203
History | View | Annotate | Download (83.7 KB)
| 1 |
//
|
|---|---|
| 2 |
// DelaunayOverlap.java
|
| 3 |
//
|
| 4 |
|
| 5 |
/*
|
| 6 |
VisAD system for interactive analysis and visualization of numerical
|
| 7 |
data. Copyright (C) 1996 - 2002 Bill Hibbard, Curtis Rueden, Tom
|
| 8 |
Rink, Dave Glowacki, Steve Emmerson, Tom Whittaker, Don Murray, and
|
| 9 |
Tommy Jasmin.
|
| 10 |
|
| 11 |
This library is free software; you can redistribute it and/or
|
| 12 |
modify it under the terms of the GNU Library General Public
|
| 13 |
License as published by the Free Software Foundation; either
|
| 14 |
version 2 of the License, or (at your option) any later version.
|
| 15 |
|
| 16 |
This library is distributed in the hope that it will be useful,
|
| 17 |
but WITHOUT ANY WARRANTY; without even the implied warranty of
|
| 18 |
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
|
| 19 |
Library General Public License for more details.
|
| 20 |
|
| 21 |
You should have received a copy of the GNU Library General Public
|
| 22 |
License along with this library; if not, write to the Free
|
| 23 |
Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
|
| 24 |
MA 02111-1307, USA
|
| 25 |
*/
|
| 26 |
|
| 27 |
package org.gvsig.fmap.algorithm.triangulation.visad; |
| 28 |
|
| 29 |
// AWT is used in main method, for testing
|
| 30 |
import java.awt.Canvas; |
| 31 |
import java.awt.Color; |
| 32 |
import java.awt.Dimension; |
| 33 |
import java.awt.Frame; |
| 34 |
import java.awt.Graphics; |
| 35 |
import java.awt.GridLayout; |
| 36 |
import java.awt.Panel; |
| 37 |
import java.awt.Toolkit; |
| 38 |
import java.awt.event.WindowAdapter; |
| 39 |
import java.awt.event.WindowEvent; |
| 40 |
import java.awt.event.WindowListener; |
| 41 |
|
| 42 |
/**
|
| 43 |
DelaunayOverlap quickly constructs a triangulation of a series
|
| 44 |
of partially overlapping two-dimensional gridded sets.<P>
|
| 45 |
|
| 46 |
Grids should be aligned so that they overlap in a roughly
|
| 47 |
vertical column (i.e., as samples[*][1] increases or decreases,
|
| 48 |
the grid number increases or decreases respectively).
|
| 49 |
DelaunayOverlap does not handle grids which overlap across a
|
| 50 |
horizontal row.<P>
|
| 51 |
*/
|
| 52 |
public class DelaunayOverlap extends Delaunay { |
| 53 |
|
| 54 |
/* This method could result in overlapping triangles in the
|
| 55 |
triangulation in the intersecting section of two grids.
|
| 56 |
For "bow-tie" shaped data, this problem is unlikely to
|
| 57 |
occur, but for data with the left and right bowed inwards,
|
| 58 |
overlapping triangles are more likely. The latter type of
|
| 59 |
data is also more likely to have a greater number of
|
| 60 |
Type 3 lost points. */
|
| 61 |
|
| 62 |
// topology signature of grid boxes
|
| 63 |
private boolean sig; |
| 64 |
|
| 65 |
/** constructor--lenx and leny are the dimensions of each gridded set. */
|
| 66 |
public DelaunayOverlap(double[][] samples, int lenx, int leny) |
| 67 |
throws VisADException {
|
| 68 |
if (samples.length < 2) { |
| 69 |
throw new SetException("DelaunayOverlap: bad dimension"); |
| 70 |
} |
| 71 |
if (lenx < 2 || leny < 2) { |
| 72 |
throw new SetException("DelaunayOverlap: must have at least " |
| 73 |
+"two points per dimension");
|
| 74 |
} |
| 75 |
if (samples[0].length < lenx*leny) { |
| 76 |
throw new SetException("DelaunayOverlap: not enough samples"); |
| 77 |
} |
| 78 |
int lenp = lenx*leny;
|
| 79 |
int leng = (lenx-1)*(leny-1); |
| 80 |
int numgrids = samples[0].length/lenp; |
| 81 |
if (numgrids < 2) { |
| 82 |
throw new SetException("DelaunayOverlap: not enough grids " |
| 83 |
+"(try Gridded2DSet instead)");
|
| 84 |
} |
| 85 |
|
| 86 |
// samples consistency check for each grid
|
| 87 |
sig = ( (samples[0][1]-samples[0][0]) |
| 88 |
*(samples[1][lenx+1]-samples[1][1]) |
| 89 |
- (samples[1][1]-samples[1][0]) |
| 90 |
*(samples[0][lenx+1]-samples[0][1]) > 0); |
| 91 |
for (int curgrid=0; curgrid<numgrids; curgrid++) { |
| 92 |
for (int gy=0; gy<leny-1; gy++) { |
| 93 |
for (int gx=0; gx<lenx-1; gx++) { |
| 94 |
double v00x, v00y, v10x, v10y, v01x, v01y, v11x, v11y;
|
| 95 |
int q = curgrid*lenp+gy*lenx+gx;
|
| 96 |
v00x = samples[0][q];
|
| 97 |
v00y = samples[1][q];
|
| 98 |
v10x = samples[0][q+1]; |
| 99 |
v10y = samples[1][q+1]; |
| 100 |
v01x = samples[0][q+lenx];
|
| 101 |
v01y = samples[1][q+lenx];
|
| 102 |
v11x = samples[0][q+lenx+1]; |
| 103 |
v11y = samples[1][q+lenx+1]; |
| 104 |
if ( ( (v10x-v00x)*(v11y-v10y)
|
| 105 |
- (v10y-v00y)*(v11x-v10x) > 0 != sig)
|
| 106 |
|| ( (v11x-v10x)*(v01y-v11y) |
| 107 |
- (v11y-v10y)*(v01x-v11x) > 0 != sig)
|
| 108 |
|| ( (v01x-v11x)*(v00y-v01y) |
| 109 |
- (v01y-v11y)*(v00x-v01x) > 0 != sig)
|
| 110 |
|| ( (v00x-v01x)*(v10y-v00y) |
| 111 |
- (v00y-v01y)*(v10x-v00x) > 0 != sig) ) {
|
| 112 |
throw new SetException("DelaunayOverlap: samples from grid " |
| 113 |
+curgrid+" do not form a valid grid "
|
| 114 |
+"("+gx+","+gy+")"); |
| 115 |
} |
| 116 |
} |
| 117 |
} |
| 118 |
} |
| 119 |
|
| 120 |
// arrays for keeping track of which points fall in which grid boxes
|
| 121 |
int ngrid = (numgrids-1)*(lenx-1)*(leny-1); |
| 122 |
int[][] grid = new int[ngrid][5]; |
| 123 |
int[] gsize = new int[ngrid]; |
| 124 |
int[] glen = new int[ngrid]; |
| 125 |
for (int g=0; g<numgrids-1; g++) { |
| 126 |
for (int gy=0; gy<leny-1; gy++) { |
| 127 |
for (int gx=0; gx<lenx-1; gx++) { |
| 128 |
int qg = leng*g + (lenx-1)*gy + gx; |
| 129 |
gsize[qg] = 4;
|
| 130 |
glen[qg] = 0;
|
| 131 |
} |
| 132 |
} |
| 133 |
} |
| 134 |
|
| 135 |
// Broken grid boxes (severed from a previous grid by a line)
|
| 136 |
// broken[*][0] = left edge break point
|
| 137 |
// broken[*][X] = break point of X'th grid box column
|
| 138 |
// broken[*][lenx] = right edge break point
|
| 139 |
int[][] broken = new int[numgrids][lenx+1]; |
| 140 |
for (int g=0; g<numgrids; g++) { |
| 141 |
for (int pix=0; pix<lenx+1; pix++) { |
| 142 |
broken[g][pix] = -1;
|
| 143 |
} |
| 144 |
} |
| 145 |
|
| 146 |
// left and right outer hulls
|
| 147 |
int[] leftedge = new int[numgrids*leny]; |
| 148 |
int ledges = 0; |
| 149 |
int[] rightedge = new int[numgrids*leny]; |
| 150 |
int redges = 0; |
| 151 |
|
| 152 |
// the trinum arrays keep track of the triangle
|
| 153 |
// numbers that border the edge arrays;
|
| 154 |
// the triedge arrays keep track of the edge numbers
|
| 155 |
// of the triangles that border the edge arrays
|
| 156 |
int[] ltrinum = new int[numgrids*leny]; |
| 157 |
int[] ltriedge = new int[numgrids*leny]; |
| 158 |
int[] rtrinum = new int[numgrids*leny]; |
| 159 |
int[] rtriedge = new int[numgrids*leny]; |
| 160 |
|
| 161 |
// arrays for saving previous values of edge information
|
| 162 |
int[] old_leftedge; |
| 163 |
int[] old_rightedge; |
| 164 |
int old_ledges;
|
| 165 |
int old_redges;
|
| 166 |
int[] old_ltrinum; |
| 167 |
int[] old_ltriedge; |
| 168 |
int[] old_rtrinum; |
| 169 |
int[] old_rtriedge; |
| 170 |
|
| 171 |
// boolean matrix of whether each point fell into a grid
|
| 172 |
boolean[] ptfell = new boolean[numgrids*lenx*leny]; |
| 173 |
for (int g=0; g<numgrids; g++) { |
| 174 |
for (int line=0; line<leny; line++) { |
| 175 |
for (int pix=0; pix<lenx; pix++) { |
| 176 |
ptfell[lenp*g+lenx*line+pix] = false;
|
| 177 |
} |
| 178 |
} |
| 179 |
} |
| 180 |
|
| 181 |
// walk construction helper arrays, for use during box triangulation
|
| 182 |
int[] bottom = new int[lenx-1]; |
| 183 |
for (int i=0; i<lenx-1; i++) bottom[i] = -1; |
| 184 |
// array for saving previous values of bottom
|
| 185 |
int[] old_bottom; |
| 186 |
|
| 187 |
// more walk helper arrays, for use during the zipping algorithm.
|
| 188 |
// tri[*][0] = right line of grid box
|
| 189 |
// tri[*][1] = bottom line of grid box
|
| 190 |
// tri[*][2] = left line of grid box
|
| 191 |
// tri[*][istwo[boxnum] ? 2 : 0] = top line of grid box
|
| 192 |
// tlr[*][0] = triangle containing top line
|
| 193 |
// tlr[*][1] = triangle containing left line
|
| 194 |
// tlr[*][2] = triangle containing right line
|
| 195 |
int[][] tlr = new int[leng][3]; |
| 196 |
boolean[] istwo = new boolean[leng]; |
| 197 |
|
| 198 |
// expandable triangle storage arrays
|
| 199 |
int[][] tri; |
| 200 |
int[][] walk; |
| 201 |
int trisize = 0; |
| 202 |
int trilen = 0; |
| 203 |
|
| 204 |
int curgrid;
|
| 205 |
int prevgrid = 0; |
| 206 |
|
| 207 |
|
| 208 |
// ******************** MAIN ALGORITHM ******************** //
|
| 209 |
|
| 210 |
// *** PHASE 1: *** Pre-triangulation preparation
|
| 211 |
|
| 212 |
// *** PHASE 1a: *** Figure out where all the points lie
|
| 213 |
// (i.e., locate containing grid boxes).
|
| 214 |
|
| 215 |
// array of pointers to last correct grid box of each grid
|
| 216 |
int[] foundx = new int[numgrids]; |
| 217 |
int[] foundy = new int[numgrids]; |
| 218 |
|
| 219 |
// center grid box of a grid
|
| 220 |
int cenx = lenx/2 - 1; |
| 221 |
int ceny = leny/2 - 1; |
| 222 |
|
| 223 |
// current point being examined in 1-D rasterization
|
| 224 |
int curpt;
|
| 225 |
|
| 226 |
// examine all grids except the first one
|
| 227 |
for (curgrid=1; curgrid<numgrids; curgrid++) { |
| 228 |
|
| 229 |
// initialize found arrays to center of each grid
|
| 230 |
for (int g=0; g<numgrids; g++) { |
| 231 |
foundx[g] = cenx; |
| 232 |
foundy[g] = ceny; |
| 233 |
} |
| 234 |
|
| 235 |
// keeps track of whether entire scan line is outside previous grid
|
| 236 |
int lfound;
|
| 237 |
// keeps track of whether entire grid is outside previous grid
|
| 238 |
int gfound;
|
| 239 |
|
| 240 |
// examine every point in the current grid
|
| 241 |
gfound = curgrid; |
| 242 |
for (int line=0; line<leny; line++) { |
| 243 |
lfound = curgrid; |
| 244 |
for (int pix=0; pix<lenx; pix++) { |
| 245 |
|
| 246 |
// the point to be located
|
| 247 |
curpt = lenp*curgrid + lenx*line + pix; |
| 248 |
double x = samples[0][curpt]; |
| 249 |
double y = samples[1][curpt]; |
| 250 |
|
| 251 |
// look through all applicable previous grids until box found
|
| 252 |
boolean found = false; |
| 253 |
for (int g=prevgrid; g<curgrid; g++) { |
| 254 |
|
| 255 |
// search grid #g for containing grid box of current point
|
| 256 |
int gx = foundx[g];
|
| 257 |
int gy = foundy[g];
|
| 258 |
int ogx, ogy;
|
| 259 |
int q, qg;
|
| 260 |
double ax, ay, bx, by, cx, cy, dx, dy;
|
| 261 |
|
| 262 |
// marching boxes--find correct grid box
|
| 263 |
while (!found) {
|
| 264 |
q = lenp*g + lenx*gy + gx; // index into samples
|
| 265 |
ax = samples[0][q];
|
| 266 |
ay = samples[1][q];
|
| 267 |
bx = samples[0][q+1]; |
| 268 |
by = samples[1][q+1]; |
| 269 |
cx = samples[0][q+lenx];
|
| 270 |
cy = samples[1][q+lenx];
|
| 271 |
dx = samples[0][q+lenx+1]; |
| 272 |
dy = samples[1][q+lenx+1]; |
| 273 |
ogx = gx; |
| 274 |
ogy = gy; |
| 275 |
|
| 276 |
// determine marching direction
|
| 277 |
switch (((ax - bx)*(y - by)
|
| 278 |
- (ay - by)*(x - bx) < 0 == sig ? 0 : 1) |
| 279 |
+ ((bx - dx)*(y - dy) |
| 280 |
- (by - dy)*(x - dx) < 0 == sig ? 0 : 2) |
| 281 |
+ ((dx - cx)*(y - cy) |
| 282 |
- (dy - cy)*(x - cx) < 0 == sig ? 0 : 4) |
| 283 |
+ ((cx - ax)*(y - ay) |
| 284 |
- (cy - ay)*(x - ax) < 0 == sig ? 0 : 8)) { |
| 285 |
|
| 286 |
case 0: // inside this box |
| 287 |
qg = leng*g + (lenx-1)*gy + gx; // index into grid |
| 288 |
found = true;
|
| 289 |
if (lfound > g) lfound = g;
|
| 290 |
if (gfound > g) gfound = g;
|
| 291 |
foundx[g] = gx; |
| 292 |
foundy[g] = gy; |
| 293 |
|
| 294 |
// mark grid box as containing current point
|
| 295 |
grid[qg][glen[qg]++] = curpt; |
| 296 |
|
| 297 |
// expand grid array as necessary
|
| 298 |
if (glen[qg] > gsize[qg]) {
|
| 299 |
int[] newg = new int[gsize[q]+gsize[q]+1]; |
| 300 |
System.arraycopy(grid[qg], 0, newg, 0, gsize[qg]); |
| 301 |
grid[qg] = newg; |
| 302 |
gsize[qg] += gsize[qg]; |
| 303 |
} |
| 304 |
|
| 305 |
// mark point as being found
|
| 306 |
ptfell[curpt] = true;
|
| 307 |
|
| 308 |
// break box to the lower left of point
|
| 309 |
if (broken[curgrid][pix] < line) {
|
| 310 |
broken[curgrid][pix] = line; |
| 311 |
} |
| 312 |
|
| 313 |
// break box to the lower right of point
|
| 314 |
if (broken[curgrid][pix+1] < line) { |
| 315 |
broken[curgrid][pix+1] = line;
|
| 316 |
} |
| 317 |
break;
|
| 318 |
|
| 319 |
case 1: case 11: // up from this box |
| 320 |
gy--; |
| 321 |
break;
|
| 322 |
case 2: case 7: // right from this box |
| 323 |
gx++; |
| 324 |
break;
|
| 325 |
case 3: // up and right from this box |
| 326 |
gx++; |
| 327 |
gy--; |
| 328 |
break;
|
| 329 |
case 4: case 14: // down from this box |
| 330 |
gy++; |
| 331 |
break;
|
| 332 |
case 6: // down and right from this box |
| 333 |
gx++; |
| 334 |
gy++; |
| 335 |
break;
|
| 336 |
case 8: case 13: // left from this box |
| 337 |
gx--; |
| 338 |
break;
|
| 339 |
case 9: // up and left from this box |
| 340 |
gx--; |
| 341 |
gy--; |
| 342 |
break;
|
| 343 |
case 12: // down and left from this box |
| 344 |
gx--; |
| 345 |
gy++; |
| 346 |
break;
|
| 347 |
|
| 348 |
default: // theoretically, should never occur |
| 349 |
throw new SetException("DelaunayOverlap: " |
| 350 |
+"pathological grid");
|
| 351 |
} |
| 352 |
|
| 353 |
if (gx > lenx-2) gx = lenx-2; |
| 354 |
if (gx < 0) gx = 0; |
| 355 |
if (gy > leny-2) gy = leny-2; |
| 356 |
if (gy < 0) gy = 0; |
| 357 |
|
| 358 |
// check if point is outside the grid
|
| 359 |
if (ogx == gx && ogy == gy && !found) break; |
| 360 |
} |
| 361 |
|
| 362 |
// no need to check remaining grids if a grid box was located
|
| 363 |
if (found) break; |
| 364 |
} |
| 365 |
} |
| 366 |
|
| 367 |
// if entire scan line wasn't in a previous grid, advance prevgrid
|
| 368 |
prevgrid = lfound; |
| 369 |
} |
| 370 |
// if entire grid wasn't in a previous grid, advance prevgrid
|
| 371 |
prevgrid = gfound; |
| 372 |
} |
| 373 |
|
| 374 |
// *** PHASE 1b: *** Break boxes that contain a point from
|
| 375 |
// the bottom line of the previous grid.
|
| 376 |
// Search logic from above is duplicated
|
| 377 |
// to maximize efficiency of the search.
|
| 378 |
|
| 379 |
// offset of the bottom scan line of a grid
|
| 380 |
int z = lenp-lenx;
|
| 381 |
|
| 382 |
// re-initialize arrays to center of each grid (old values are
|
| 383 |
// irrelevant now; center is actually a better starting guess)
|
| 384 |
for (int g=0; g<numgrids; g++) { |
| 385 |
foundx[g] = cenx; |
| 386 |
foundy[g] = ceny; |
| 387 |
} |
| 388 |
|
| 389 |
// search all grids but the last one
|
| 390 |
for (curgrid=0; curgrid<numgrids-1; curgrid++) { |
| 391 |
|
| 392 |
// search every point in the bottom scan line of curgrid
|
| 393 |
for (int pix=0; pix<lenx; pix++) { |
| 394 |
|
| 395 |
// the point to be located
|
| 396 |
curpt = lenp*curgrid + z + pix; |
| 397 |
double x = samples[0][curpt]; |
| 398 |
double y = samples[1][curpt]; |
| 399 |
int g = curgrid+1; |
| 400 |
|
| 401 |
// search grid #g for containing grid box of current point
|
| 402 |
int gx = foundx[g];
|
| 403 |
int gy = foundy[g];
|
| 404 |
int ogx, ogy, q;
|
| 405 |
boolean found = false; |
| 406 |
double ax, ay, bx, by, cx, cy, dx, dy;
|
| 407 |
|
| 408 |
// marching boxes--find correct grid box
|
| 409 |
while (!found) {
|
| 410 |
q = lenp*g + lenx*gy + gx; |
| 411 |
ax = samples[0][q];
|
| 412 |
ay = samples[1][q];
|
| 413 |
bx = samples[0][q+1]; |
| 414 |
by = samples[1][q+1]; |
| 415 |
cx = samples[0][q+lenx];
|
| 416 |
cy = samples[1][q+lenx];
|
| 417 |
dx = samples[0][q+lenx+1]; |
| 418 |
dy = samples[1][q+lenx+1]; |
| 419 |
ogx = gx; |
| 420 |
ogy = gy; |
| 421 |
|
| 422 |
// determine marching direction
|
| 423 |
switch (((ax - bx)*(y - by)
|
| 424 |
- (ay - by)*(x - bx) < 0 == sig ? 0 : 1) |
| 425 |
+ ((bx - dx)*(y - dy) |
| 426 |
- (by - dy)*(x - dx) < 0 == sig ? 0 : 2) |
| 427 |
+ ((dx - cx)*(y - cy) |
| 428 |
- (dy - cy)*(x - cx) < 0 == sig ? 0 : 4) |
| 429 |
+ ((cx - ax)*(y - ay) |
| 430 |
- (cy - ay)*(x - ax) < 0 == sig ? 0 : 8)) { |
| 431 |
|
| 432 |
case 0: // inside this box |
| 433 |
found = true;
|
| 434 |
foundx[g] = gx; |
| 435 |
foundy[g] = gy; |
| 436 |
|
| 437 |
// break containing box
|
| 438 |
if (broken[g][gx+1] < gy) broken[g][gx+1] = gy; |
| 439 |
break;
|
| 440 |
|
| 441 |
case 1: case 11: // up from this box |
| 442 |
gy--; |
| 443 |
break;
|
| 444 |
case 2: case 7: // right from this box |
| 445 |
gx++; |
| 446 |
break;
|
| 447 |
case 3: // up and right from this box |
| 448 |
gx++; |
| 449 |
gy--; |
| 450 |
break;
|
| 451 |
case 4: case 14: // down from this box |
| 452 |
gy++; |
| 453 |
break;
|
| 454 |
case 6: // down and right from this box |
| 455 |
gx++; |
| 456 |
gy++; |
| 457 |
break;
|
| 458 |
case 8: case 13: // left from this box |
| 459 |
gx--; |
| 460 |
break;
|
| 461 |
case 9: // up and left from this box |
| 462 |
gx--; |
| 463 |
gy--; |
| 464 |
break;
|
| 465 |
case 12: // down and left from this box |
| 466 |
gx--; |
| 467 |
gy++; |
| 468 |
break;
|
| 469 |
|
| 470 |
default: // theoretically, should never occur |
| 471 |
throw new SetException("DelaunayOverlap: " |
| 472 |
+"pathological grid");
|
| 473 |
} |
| 474 |
if (gx > lenx-2) gx = lenx-2; |
| 475 |
if (gx < 0) gx = 0; |
| 476 |
if (gy > leny-2) gy = leny-2; |
| 477 |
if (gy < 0) gy = 0; |
| 478 |
|
| 479 |
// check if point is outside the grid
|
| 480 |
if (ogx == gx && ogy == gy && !found) break; |
| 481 |
} |
| 482 |
} |
| 483 |
} |
| 484 |
|
| 485 |
// *** PHASE 2: *** triangulate points inside grid boxes
|
| 486 |
|
| 487 |
// set leftedge to the first column of grid 0
|
| 488 |
leftedge[0] = 0; |
| 489 |
for (int i=1; i<leny; i++) { |
| 490 |
leftedge[i] = leftedge[i-1] + lenx;
|
| 491 |
} |
| 492 |
ledges += leny; |
| 493 |
|
| 494 |
// set rightedge to the last column of grid 0
|
| 495 |
rightedge[0] = lenx-1; |
| 496 |
for (int i=1; i<leny; i++) { |
| 497 |
rightedge[i] = rightedge[i-1] + lenx;
|
| 498 |
} |
| 499 |
redges += leny; |
| 500 |
|
| 501 |
// determine necessary tri array size
|
| 502 |
// Note: This memory allocation strategy could be more precise.
|
| 503 |
// It should loop through only unbroken boxes, and figure
|
| 504 |
// out the number of triangles in the zipped sections and
|
| 505 |
// the number of lost points. If these steps were
|
| 506 |
// implemented, then the tri and walk arrays would not be
|
| 507 |
// needed, and triangulation data could be stored directly
|
| 508 |
// in the Tri and Walk arrays.
|
| 509 |
for (curgrid=0; curgrid<numgrids; curgrid++) { |
| 510 |
for (int gy=0; gy<leny-1; gy++) { |
| 511 |
for (int gx=0; gx<lenx-1; gx++) { |
| 512 |
int npts = (curgrid == numgrids-1) |
| 513 |
? 0 : glen[leng*curgrid + (lenx-1)*gy + gx]; |
| 514 |
trisize += 2*npts+2; |
| 515 |
} |
| 516 |
} |
| 517 |
} |
| 518 |
// probably won't be more than 2*leny Type 2 lost points per grid
|
| 519 |
trisize += 2*leny*numgrids;
|
| 520 |
|
| 521 |
// allocate memory
|
| 522 |
tri = new int[trisize+2][3]; |
| 523 |
walk = new int[trisize+2][3]; |
| 524 |
for (int i=0; i<trisize+2; i++) { |
| 525 |
walk[i][0] = -1; |
| 526 |
walk[i][1] = -1; |
| 527 |
walk[i][2] = -1; |
| 528 |
} |
| 529 |
|
| 530 |
// examine every grid
|
| 531 |
for (curgrid=0; curgrid<numgrids; curgrid++) { |
| 532 |
|
| 533 |
// save old bottom values
|
| 534 |
old_bottom = bottom; |
| 535 |
bottom = new int[lenx-1]; |
| 536 |
for (int i=0; i<lenx-1; i++) bottom[i] = -1; |
| 537 |
|
| 538 |
// triangulate all unbroken boxes in this grid
|
| 539 |
for (int gx=0; gx<lenx-1; gx++) { |
| 540 |
for (int gy=broken[curgrid][gx+1]+1; gy<leny-1; gy++) { |
| 541 |
|
| 542 |
int q = lenp*curgrid + lenx*gy + gx;
|
| 543 |
int qmg = (lenx-1)*gy + gx; |
| 544 |
int qg = leng*curgrid + qmg;
|
| 545 |
int A = q;
|
| 546 |
int B = q + 1; |
| 547 |
int C = q + lenx;
|
| 548 |
int D = q + lenx + 1; |
| 549 |
|
| 550 |
double ax, ay, bx, by, cx, cy, dx, dy;
|
| 551 |
int G1, G2, G3;
|
| 552 |
|
| 553 |
switch (curgrid == numgrids-1 ? 0 : glen[qg]) { |
| 554 |
case 0: // find the best diagonal of the box |
| 555 |
ax = samples[0][A];
|
| 556 |
ay = samples[1][A];
|
| 557 |
bx = samples[0][B];
|
| 558 |
by = samples[1][B];
|
| 559 |
cx = samples[0][C];
|
| 560 |
cy = samples[1][C];
|
| 561 |
dx = samples[0][D];
|
| 562 |
dy = samples[1][D];
|
| 563 |
double abx = ax - bx;
|
| 564 |
double aby = ay - by;
|
| 565 |
double acx = ax - cx;
|
| 566 |
double acy = ay - cy;
|
| 567 |
double dbx = dx - bx;
|
| 568 |
double dby = dy - by;
|
| 569 |
double dcx = dx - cx;
|
| 570 |
double dcy = dy - cy;
|
| 571 |
double Q = abx*acx + aby*acy;
|
| 572 |
double R = dbx*abx + dby*aby;
|
| 573 |
double S = acx*dcx + acy*dcy;
|
| 574 |
double T = dbx*dcx + dby*dcy;
|
| 575 |
boolean diag;
|
| 576 |
if (Q < 0 && T < 0 || R > 0 && S > 0) diag = true; |
| 577 |
else if (R < 0 && S < 0 || Q > 0 && T > 0) diag = false; |
| 578 |
else if ((Q < 0 ? Q : T) < (R < 0 ? R : S)) diag = true; |
| 579 |
else diag = false; |
| 580 |
|
| 581 |
if (diag) {
|
| 582 |
tri[trilen][0] = B;
|
| 583 |
tri[trilen][1] = D;
|
| 584 |
tri[trilen][2] = A;
|
| 585 |
walk[trilen][2] = bottom[gx];
|
| 586 |
if (bottom[gx] >= 0) walk[bottom[gx]][1] = trilen; |
| 587 |
walk[trilen][1] = trilen+1; |
| 588 |
trilen++; |
| 589 |
tri[trilen][0] = A;
|
| 590 |
tri[trilen][1] = D;
|
| 591 |
tri[trilen][2] = C;
|
| 592 |
if (gx > 0 && broken[curgrid][gx] < gy) { |
| 593 |
walk[trilen][2] = tlr[qmg-1][2]; |
| 594 |
walk[tlr[qmg-1][2]][0] = trilen; |
| 595 |
} |
| 596 |
else walk[trilen][2] = -1; |
| 597 |
walk[trilen][0] = trilen-1; |
| 598 |
trilen++; |
| 599 |
bottom[gx] = trilen-1;
|
| 600 |
tlr[qmg][0] = trilen-2; |
| 601 |
tlr[qmg][1] = trilen-1; |
| 602 |
tlr[qmg][2] = trilen-2; |
| 603 |
istwo[qmg] = true;
|
| 604 |
} |
| 605 |
else {
|
| 606 |
tri[trilen][0] = A;
|
| 607 |
tri[trilen][1] = B;
|
| 608 |
tri[trilen][2] = C;
|
| 609 |
walk[trilen][0] = bottom[gx];
|
| 610 |
if (bottom[gx] >= 0) walk[bottom[gx]][1] = trilen; |
| 611 |
if (gx > 0 && broken[curgrid][gx] < gy) { |
| 612 |
walk[trilen][2] = tlr[qmg-1][2]; |
| 613 |
walk[tlr[qmg-1][2]][0] = trilen; |
| 614 |
} |
| 615 |
else walk[trilen][2] = -1; |
| 616 |
walk[trilen][1] = trilen+1; |
| 617 |
trilen++; |
| 618 |
tri[trilen][0] = B;
|
| 619 |
tri[trilen][1] = D;
|
| 620 |
tri[trilen][2] = C;
|
| 621 |
walk[trilen][2] = trilen-1; |
| 622 |
trilen++; |
| 623 |
bottom[gx] = trilen-1;
|
| 624 |
tlr[qmg][0] = trilen-2; |
| 625 |
tlr[qmg][1] = trilen-2; |
| 626 |
tlr[qmg][2] = trilen-1; |
| 627 |
istwo[qmg] = false;
|
| 628 |
} |
| 629 |
break;
|
| 630 |
|
| 631 |
case 1: // draw the four triangles of the point |
| 632 |
G1 = grid[qg][0];
|
| 633 |
tri[trilen][0] = B;
|
| 634 |
tri[trilen][1] = G1;
|
| 635 |
tri[trilen][2] = A;
|
| 636 |
walk[trilen][2] = bottom[gx];
|
| 637 |
if (bottom[gx] >= 0) walk[bottom[gx]][1] = trilen; |
| 638 |
walk[trilen][0] = trilen+1; |
| 639 |
walk[trilen][1] = trilen+3; |
| 640 |
trilen++; |
| 641 |
tri[trilen][0] = B;
|
| 642 |
tri[trilen][1] = D;
|
| 643 |
tri[trilen][2] = G1;
|
| 644 |
walk[trilen][1] = trilen+1; |
| 645 |
walk[trilen][2] = trilen-1; |
| 646 |
trilen++; |
| 647 |
tri[trilen][0] = G1;
|
| 648 |
tri[trilen][1] = D;
|
| 649 |
tri[trilen][2] = C;
|
| 650 |
walk[trilen][0] = trilen-1; |
| 651 |
walk[trilen][2] = trilen+1; |
| 652 |
trilen++; |
| 653 |
tri[trilen][0] = A;
|
| 654 |
tri[trilen][1] = G1;
|
| 655 |
tri[trilen][2] = C;
|
| 656 |
if (gx > 0 && broken[curgrid][gx] < gy) { |
| 657 |
walk[trilen][2] = tlr[qmg-1][2]; |
| 658 |
walk[tlr[qmg-1][2]][0] = trilen; |
| 659 |
} |
| 660 |
else walk[trilen][2] = -1; |
| 661 |
walk[trilen][0] = trilen-3; |
| 662 |
walk[trilen][1] = trilen-1; |
| 663 |
trilen++; |
| 664 |
bottom[gx] = trilen-2;
|
| 665 |
tlr[qmg][0] = trilen-4; |
| 666 |
tlr[qmg][1] = trilen-1; |
| 667 |
tlr[qmg][2] = trilen-3; |
| 668 |
istwo[qmg] = true;
|
| 669 |
break;
|
| 670 |
|
| 671 |
default: // triangulate first two points |
| 672 |
// use a diagonal comparison
|
| 673 |
G1 = grid[qg][0];
|
| 674 |
G2 = grid[qg][1];
|
| 675 |
double Gdx = samples[0][G2] - samples[0][G1]; |
| 676 |
double Gdy = samples[1][G2] - samples[1][G1]; |
| 677 |
ax = samples[0][A];
|
| 678 |
ay = samples[1][A];
|
| 679 |
bx = samples[0][B];
|
| 680 |
by = samples[1][B];
|
| 681 |
cx = samples[0][C];
|
| 682 |
cy = samples[1][C];
|
| 683 |
dx = samples[0][D];
|
| 684 |
dy = samples[1][D];
|
| 685 |
double val1 = Gdx*(ax-dx) + Gdy*(ay-dy);
|
| 686 |
if (val1 < 0) val1 = -val1; |
| 687 |
double val2 = Gdx*(bx-cx) + Gdy*(by-cy);
|
| 688 |
if (val2 < 0) val2 = -val2; |
| 689 |
|
| 690 |
if (val1 > val2) {
|
| 691 |
double Qx1 = ax - samples[0][G1]; |
| 692 |
double Qy1 = ay - samples[1][G1]; |
| 693 |
double Qx2 = ax - samples[0][G2]; |
| 694 |
double Qy2 = ay - samples[1][G2]; |
| 695 |
if (Qx1*Qx1 + Qy1*Qy1 < Qx2*Qx2 + Qy2*Qy2) {
|
| 696 |
// G1 is closer to A than G2
|
| 697 |
tri[trilen][0] = B;
|
| 698 |
tri[trilen][1] = G1;
|
| 699 |
tri[trilen][2] = A;
|
| 700 |
walk[trilen][2] = bottom[gx];
|
| 701 |
if (bottom[gx] >= 0) walk[bottom[gx]][1] = trilen; |
| 702 |
walk[trilen][0] = trilen+5; |
| 703 |
walk[trilen][1] = trilen+1; |
| 704 |
trilen++; |
| 705 |
tri[trilen][0] = A;
|
| 706 |
tri[trilen][1] = G1;
|
| 707 |
tri[trilen][2] = C;
|
| 708 |
if (gx > 0 && broken[curgrid][gx] < gy) { |
| 709 |
walk[trilen][2] = tlr[qmg-1][2]; |
| 710 |
walk[tlr[qmg-1][2]][0] = trilen; |
| 711 |
} |
| 712 |
else walk[trilen][2] = -1; |
| 713 |
walk[trilen][0] = trilen-1; |
| 714 |
walk[trilen][1] = trilen+3; |
| 715 |
trilen++; |
| 716 |
tri[trilen][0] = G2;
|
| 717 |
tri[trilen][1] = D;
|
| 718 |
tri[trilen][2] = C;
|
| 719 |
walk[trilen][0] = trilen+1; |
| 720 |
walk[trilen][2] = trilen+2; |
| 721 |
trilen++; |
| 722 |
tri[trilen][0] = B;
|
| 723 |
tri[trilen][1] = D;
|
| 724 |
tri[trilen][2] = G2;
|
| 725 |
walk[trilen][1] = trilen-1; |
| 726 |
walk[trilen][2] = trilen+2; |
| 727 |
trilen++; |
| 728 |
tri[trilen][0] = C;
|
| 729 |
tri[trilen][1] = G1;
|
| 730 |
tri[trilen][2] = G2;
|
| 731 |
walk[trilen][0] = trilen-3; |
| 732 |
walk[trilen][1] = trilen+1; |
| 733 |
walk[trilen][2] = trilen-2; |
| 734 |
trilen++; |
| 735 |
tri[trilen][0] = B;
|
| 736 |
tri[trilen][1] = G2;
|
| 737 |
tri[trilen][2] = G1;
|
| 738 |
walk[trilen][0] = trilen-2; |
| 739 |
walk[trilen][1] = trilen-1; |
| 740 |
walk[trilen][2] = trilen-5; |
| 741 |
trilen++; |
| 742 |
bottom[gx] = trilen-4;
|
| 743 |
tlr[qmg][0] = trilen-6; |
| 744 |
tlr[qmg][1] = trilen-5; |
| 745 |
tlr[qmg][2] = trilen-3; |
| 746 |
istwo[qmg] = true;
|
| 747 |
} |
| 748 |
else {
|
| 749 |
// G2 is closer to A than G1
|
| 750 |
tri[trilen][0] = B;
|
| 751 |
tri[trilen][1] = G2;
|
| 752 |
tri[trilen][2] = A;
|
| 753 |
walk[trilen][2] = bottom[gx];
|
| 754 |
if (bottom[gx] >= 0) walk[bottom[gx]][1] = trilen; |
| 755 |
walk[trilen][0] = trilen+5; |
| 756 |
walk[trilen][1] = trilen+1; |
| 757 |
trilen++; |
| 758 |
tri[trilen][0] = A;
|
| 759 |
tri[trilen][1] = G2;
|
| 760 |
tri[trilen][2] = C;
|
| 761 |
if (gx > 0 && broken[curgrid][gx] < gy) { |
| 762 |
walk[trilen][2] = tlr[qmg-1][2]; |
| 763 |
walk[tlr[qmg-1][2]][0] = trilen; |
| 764 |
} |
| 765 |
else walk[trilen][2] = -1; |
| 766 |
walk[trilen][0] = trilen-1; |
| 767 |
walk[trilen][1] = trilen+3; |
| 768 |
trilen++; |
| 769 |
tri[trilen][0] = G1;
|
| 770 |
tri[trilen][1] = D;
|
| 771 |
tri[trilen][2] = C;
|
| 772 |
walk[trilen][0] = trilen+1; |
| 773 |
walk[trilen][2] = trilen+2; |
| 774 |
trilen++; |
| 775 |
tri[trilen][0] = B;
|
| 776 |
tri[trilen][1] = D;
|
| 777 |
tri[trilen][2] = G1;
|
| 778 |
walk[trilen][1] = trilen-1; |
| 779 |
walk[trilen][2] = trilen+2; |
| 780 |
trilen++; |
| 781 |
tri[trilen][0] = C;
|
| 782 |
tri[trilen][1] = G2;
|
| 783 |
tri[trilen][2] = G1;
|
| 784 |
walk[trilen][0] = trilen-3; |
| 785 |
walk[trilen][1] = trilen+1; |
| 786 |
walk[trilen][2] = trilen-2; |
| 787 |
trilen++; |
| 788 |
tri[trilen][0] = B;
|
| 789 |
tri[trilen][1] = G1;
|
| 790 |
tri[trilen][2] = G2;
|
| 791 |
walk[trilen][0] = trilen-2; |
| 792 |
walk[trilen][1] = trilen-1; |
| 793 |
walk[trilen][2] = trilen-5; |
| 794 |
trilen++; |
| 795 |
bottom[gx] = trilen-4;
|
| 796 |
tlr[qmg][0] = trilen-6; |
| 797 |
tlr[qmg][1] = trilen-5; |
| 798 |
tlr[qmg][2] = trilen-3; |
| 799 |
istwo[qmg] = true;
|
| 800 |
} |
| 801 |
} |
| 802 |
else {
|
| 803 |
double Qx1 = bx - samples[0][G1]; |
| 804 |
double Qy1 = by - samples[1][G1]; |
| 805 |
double Qx2 = bx - samples[0][G2]; |
| 806 |
double Qy2 = by - samples[1][G2]; |
| 807 |
if (Qx1*Qx1 + Qy1*Qy1 < Qx2*Qx2 + Qy2*Qy2) {
|
| 808 |
// G1 is closer to B than G2
|
| 809 |
tri[trilen][0] = B;
|
| 810 |
tri[trilen][1] = G1;
|
| 811 |
tri[trilen][2] = A;
|
| 812 |
walk[trilen][2] = bottom[gx];
|
| 813 |
if (bottom[gx] >= 0) walk[bottom[gx]][1] = trilen; |
| 814 |
walk[trilen][0] = trilen+3; |
| 815 |
walk[trilen][1] = trilen+4; |
| 816 |
trilen++; |
| 817 |
tri[trilen][0] = A;
|
| 818 |
tri[trilen][1] = G2;
|
| 819 |
tri[trilen][2] = C;
|
| 820 |
if (gx > 0 && broken[curgrid][gx] < gy) { |
| 821 |
walk[trilen][2] = tlr[qmg-1][2]; |
| 822 |
walk[tlr[qmg-1][2]][0] = trilen; |
| 823 |
} |
| 824 |
else walk[trilen][2] = -1; |
| 825 |
walk[trilen][0] = trilen+3; |
| 826 |
walk[trilen][1] = trilen+1; |
| 827 |
trilen++; |
| 828 |
tri[trilen][0] = G2;
|
| 829 |
tri[trilen][1] = D;
|
| 830 |
tri[trilen][2] = C;
|
| 831 |
walk[trilen][0] = trilen+3; |
| 832 |
walk[trilen][2] = trilen-1; |
| 833 |
trilen++; |
| 834 |
tri[trilen][0] = B;
|
| 835 |
tri[trilen][1] = D;
|
| 836 |
tri[trilen][2] = G1;
|
| 837 |
walk[trilen][1] = trilen+2; |
| 838 |
walk[trilen][2] = trilen-3; |
| 839 |
trilen++; |
| 840 |
tri[trilen][0] = A;
|
| 841 |
tri[trilen][1] = G1;
|
| 842 |
tri[trilen][2] = G2;
|
| 843 |
walk[trilen][0] = trilen-4; |
| 844 |
walk[trilen][1] = trilen+1; |
| 845 |
walk[trilen][2] = trilen-3; |
| 846 |
trilen++; |
| 847 |
tri[trilen][0] = D;
|
| 848 |
tri[trilen][1] = G2;
|
| 849 |
tri[trilen][2] = G1;
|
| 850 |
walk[trilen][0] = trilen-3; |
| 851 |
walk[trilen][1] = trilen-1; |
| 852 |
walk[trilen][2] = trilen-2; |
| 853 |
trilen++; |
| 854 |
bottom[gx] = trilen-4;
|
| 855 |
tlr[qmg][0] = trilen-6; |
| 856 |
tlr[qmg][1] = trilen-5; |
| 857 |
tlr[qmg][2] = trilen-3; |
| 858 |
istwo[qmg] = true;
|
| 859 |
} |
| 860 |
else {
|
| 861 |
// G2 is closer to B than G1
|
| 862 |
tri[trilen][0] = B;
|
| 863 |
tri[trilen][1] = G2;
|
| 864 |
tri[trilen][2] = A;
|
| 865 |
walk[trilen][2] = bottom[gx];
|
| 866 |
if (bottom[gx] >= 0) walk[bottom[gx]][1] = trilen; |
| 867 |
walk[trilen][0] = trilen+3; |
| 868 |
walk[trilen][1] = trilen+4; |
| 869 |
trilen++; |
| 870 |
tri[trilen][0] = A;
|
| 871 |
tri[trilen][1] = G1;
|
| 872 |
tri[trilen][2] = C;
|
| 873 |
if (gx > 0 && broken[curgrid][gx] < gy) { |
| 874 |
walk[trilen][2] = tlr[qmg-1][2]; |
| 875 |
walk[tlr[qmg-1][2]][0] = trilen; |
| 876 |
} |
| 877 |
else walk[trilen][2] = -1; |
| 878 |
walk[trilen][0] = trilen+3; |
| 879 |
walk[trilen][1] = trilen+1; |
| 880 |
trilen++; |
| 881 |
tri[trilen][0] = G1;
|
| 882 |
tri[trilen][1] = D;
|
| 883 |
tri[trilen][2] = C;
|
| 884 |
walk[trilen][0] = trilen+3; |
| 885 |
walk[trilen][2] = trilen-1; |
| 886 |
trilen++; |
| 887 |
tri[trilen][0] = B;
|
| 888 |
tri[trilen][1] = D;
|
| 889 |
tri[trilen][2] = G2;
|
| 890 |
walk[trilen][1] = trilen+2; |
| 891 |
walk[trilen][2] = trilen-3; |
| 892 |
trilen++; |
| 893 |
tri[trilen][0] = A;
|
| 894 |
tri[trilen][1] = G2;
|
| 895 |
tri[trilen][2] = G1;
|
| 896 |
walk[trilen][0] = trilen-4; |
| 897 |
walk[trilen][1] = trilen+1; |
| 898 |
walk[trilen][2] = trilen-3; |
| 899 |
trilen++; |
| 900 |
tri[trilen][0] = D;
|
| 901 |
tri[trilen][1] = G1;
|
| 902 |
tri[trilen][2] = G2;
|
| 903 |
walk[trilen][0] = trilen-3; |
| 904 |
walk[trilen][1] = trilen-1; |
| 905 |
walk[trilen][2] = trilen-2; |
| 906 |
trilen++; |
| 907 |
bottom[gx] = trilen-4;
|
| 908 |
tlr[qmg][0] = trilen-6; |
| 909 |
tlr[qmg][1] = trilen-5; |
| 910 |
tlr[qmg][2] = trilen-3; |
| 911 |
istwo[qmg] = true;
|
| 912 |
} |
| 913 |
} |
| 914 |
|
| 915 |
// subtriangulate remaining points
|
| 916 |
int starttri = trilen-1; |
| 917 |
int maxit = 2*glen[qg]+2; |
| 918 |
for (int i=2; i<glen[qg]; i++) { |
| 919 |
int pt = grid[qg][i];
|
| 920 |
|
| 921 |
// find containing triangle of point pt
|
| 922 |
int curtri = starttri;
|
| 923 |
int p0 = -1; |
| 924 |
int p1 = -1; |
| 925 |
int p2 = -1; |
| 926 |
int w0 = -1; |
| 927 |
int w1 = -1; |
| 928 |
int w2 = -1; |
| 929 |
boolean found = false; |
| 930 |
int itnum;
|
| 931 |
for (itnum=0; itnum<maxit && !found; itnum++) { |
| 932 |
p0 = tri[curtri][0];
|
| 933 |
p1 = tri[curtri][1];
|
| 934 |
p2 = tri[curtri][2];
|
| 935 |
w0 = walk[curtri][0];
|
| 936 |
w1 = walk[curtri][1];
|
| 937 |
w2 = walk[curtri][2];
|
| 938 |
double ptx = samples[0][pt]; |
| 939 |
double pty = samples[1][pt]; |
| 940 |
double p0x = samples[0][p0]; |
| 941 |
double p0y = samples[1][p0]; |
| 942 |
double p1x = samples[0][p1]; |
| 943 |
double p1y = samples[1][p1]; |
| 944 |
double p2x = samples[0][p2]; |
| 945 |
double p2y = samples[1][p2]; |
| 946 |
double p01x = p0x-p1x;
|
| 947 |
double p01y = p0y-p1y;
|
| 948 |
double p02x = p0x-p2x;
|
| 949 |
double p02y = p0y-p2y;
|
| 950 |
double p12x = p1x-p2x;
|
| 951 |
double p12y = p1y-p2y;
|
| 952 |
double c0 = p01x*(p0y-pty) - p01y*(p0x-ptx);
|
| 953 |
double c1 = p12x*(p1y-pty) - p12y*(p1x-ptx);
|
| 954 |
double c2 = p02x*(pty-p2y) - p02y*(ptx-p2x);
|
| 955 |
boolean t0 = ( c0 == 0 || |
| 956 |
c0 > 0 == p01x*p02y - p01y*p02x > 0 ); |
| 957 |
boolean t1 = ( c1 == 0 || |
| 958 |
c1 > 0 == p12x*p01y - p12y*p01x > 0 ); |
| 959 |
boolean t2 = ( c2 == 0 || |
| 960 |
c2 > 0 == p02x*p12y - p02y*p12x > 0 ); |
| 961 |
|
| 962 |
if (!t0 && !t1 && !t2) {
|
| 963 |
throw new SetException("DelaunayOverlap: " |
| 964 |
+"subtriangulation error");
|
| 965 |
} |
| 966 |
else if (!t0) { |
| 967 |
if (curtri != tlr[qmg][2]) curtri = w0; |
| 968 |
else throw new SetException("DelaunayOverlap: " |
| 969 |
+"subtriangulation error");
|
| 970 |
} |
| 971 |
else if (!t1) { |
| 972 |
if (curtri != bottom[gx]) curtri = w1;
|
| 973 |
else throw new SetException("DelaunayOverlap: " |
| 974 |
+"subtriangulation error");
|
| 975 |
} |
| 976 |
else if (!t2) { |
| 977 |
if (curtri != tlr[qmg][0] && curtri != tlr[qmg][1]) { |
| 978 |
curtri = w2; |
| 979 |
} |
| 980 |
else throw new SetException("DelaunayOverlap: " |
| 981 |
+"subtriangulation error");
|
| 982 |
} |
| 983 |
else found = true; |
| 984 |
} |
| 985 |
|
| 986 |
// should never happen
|
| 987 |
if (!found) throw new SetException("DelaunayOverlap: " |
| 988 |
+"subtriangulation error");
|
| 989 |
else if (curtri == bottom[gx]) { |
| 990 |
// curtri is on bottom edge of grid box
|
| 991 |
|
| 992 |
// find adjoining triangles' connecting edges
|
| 993 |
int we0 = -1; |
| 994 |
int we2 = -1; |
| 995 |
for (int w=0; w<3; w++) { |
| 996 |
if (walk[w0][w] == curtri) we0 = w;
|
| 997 |
if (walk[w2][w] == curtri) we2 = w;
|
| 998 |
} |
| 999 |
if (we0 < 0 || we2 < 0) { |
| 1000 |
throw new SetException("DelaunayOverlap: " |
| 1001 |
+"subtriangulation error");
|
| 1002 |
} |
| 1003 |
|
| 1004 |
// define three new triangles
|
| 1005 |
// the first new triangle overwrites the old triangle
|
| 1006 |
tri[curtri][0] = pt;
|
| 1007 |
// tri[curtri][1] = p1; <-- already true
|
| 1008 |
// tri[curtri][2] = p2; <-- already true
|
| 1009 |
walk[curtri][0] = trilen;
|
| 1010 |
// walk[curtri][1] = w1; <-- already true
|
| 1011 |
// walk[w1][we1] = curtri; <-- already true
|
| 1012 |
walk[curtri][2] = trilen+1; |
| 1013 |
|
| 1014 |
tri[trilen][0] = p1;
|
| 1015 |
tri[trilen][1] = pt;
|
| 1016 |
tri[trilen][2] = p0;
|
| 1017 |
walk[trilen][0] = curtri;
|
| 1018 |
walk[trilen][1] = trilen+1; |
| 1019 |
walk[trilen][2] = w0;
|
| 1020 |
walk[w0][we0] = trilen; |
| 1021 |
trilen++; |
| 1022 |
|
| 1023 |
tri[trilen][0] = p0;
|
| 1024 |
tri[trilen][1] = pt;
|
| 1025 |
tri[trilen][2] = p2;
|
| 1026 |
walk[trilen][0] = trilen-1; |
| 1027 |
walk[trilen][1] = curtri;
|
| 1028 |
walk[trilen][2] = w2;
|
| 1029 |
walk[w2][we2] = trilen; |
| 1030 |
trilen++; |
| 1031 |
} |
| 1032 |
else if (curtri == tlr[qmg][0] || curtri == tlr[qmg][1]) { |
| 1033 |
// tlr[qmg][0]: curtri is on top edge of grid box
|
| 1034 |
// tlr[qmg][1]: curtri is on left edge of grid box
|
| 1035 |
|
| 1036 |
// find adjoining triangles' connecting edges
|
| 1037 |
int we0 = -1; |
| 1038 |
int we1 = -1; |
| 1039 |
for (int w=0; w<3; w++) { |
| 1040 |
if (walk[w0][w] == curtri) we0 = w;
|
| 1041 |
if (walk[w1][w] == curtri) we1 = w;
|
| 1042 |
} |
| 1043 |
if (we0 < 0 || we1 < 0) { |
| 1044 |
throw new SetException("DelaunayOverlap: " |
| 1045 |
+"subtriangulation error");
|
| 1046 |
} |
| 1047 |
|
| 1048 |
// define three new triangles
|
| 1049 |
// the first new triangle overwrites the old triangle
|
| 1050 |
// tri[curtri][0] = p0; <-- already true
|
| 1051 |
tri[curtri][1] = pt;
|
| 1052 |
// tri[curtri][2] = p2; <-- already true
|
| 1053 |
walk[curtri][0] = trilen;
|
| 1054 |
walk[curtri][1] = trilen+1; |
| 1055 |
// walk[curtri][2] = w2; <-- already true
|
| 1056 |
// walk[w2][we2] = curtri; <-- already true
|
| 1057 |
|
| 1058 |
tri[trilen][0] = p0;
|
| 1059 |
tri[trilen][1] = p1;
|
| 1060 |
tri[trilen][2] = pt;
|
| 1061 |
walk[trilen][0] = w0;
|
| 1062 |
walk[w0][we0] = trilen; |
| 1063 |
walk[trilen][1] = trilen+1; |
| 1064 |
walk[trilen][2] = curtri;
|
| 1065 |
trilen++; |
| 1066 |
|
| 1067 |
tri[trilen][0] = pt;
|
| 1068 |
tri[trilen][1] = p1;
|
| 1069 |
tri[trilen][2] = p2;
|
| 1070 |
walk[trilen][0] = trilen-1; |
| 1071 |
walk[trilen][1] = w1;
|
| 1072 |
walk[w1][we1] = trilen; |
| 1073 |
walk[trilen][2] = curtri;
|
| 1074 |
trilen++; |
| 1075 |
} |
| 1076 |
else {
|
| 1077 |
// if curtri == tlr[qmg][2],
|
| 1078 |
// then curtri is on right edge of grid box,
|
| 1079 |
// otherwise, curtri is not on border of grid box
|
| 1080 |
|
| 1081 |
// find adjoining triangles' connecting edges
|
| 1082 |
int we1 = -1; |
| 1083 |
int we2 = -1; |
| 1084 |
for (int w=0; w<3; w++) { |
| 1085 |
if (walk[w1][w] == curtri) we1 = w;
|
| 1086 |
if (walk[w2][w] == curtri) we2 = w;
|
| 1087 |
} |
| 1088 |
if (we1 < 0 || we2 < 0) { |
| 1089 |
throw new SetException("DelaunayOverlap: " |
| 1090 |
+"subtriangulation error");
|
| 1091 |
} |
| 1092 |
|
| 1093 |
// define three new triangles
|
| 1094 |
// the first new triangle overwrites the old triangle
|
| 1095 |
// tri[curtri][0] = p0; <-- already true
|
| 1096 |
// tri[curtri][1] = p1; <-- already true
|
| 1097 |
tri[curtri][2] = pt;
|
| 1098 |
// walk[curtri][0] = w0; <-- already true
|
| 1099 |
// walk[w0][we0] = curtri; <-- already true
|
| 1100 |
walk[curtri][1] = trilen;
|
| 1101 |
walk[curtri][2] = trilen+1; |
| 1102 |
|
| 1103 |
tri[trilen][0] = pt;
|
| 1104 |
tri[trilen][1] = p1;
|
| 1105 |
tri[trilen][2] = p2;
|
| 1106 |
walk[trilen][0] = curtri;
|
| 1107 |
walk[trilen][1] = w1;
|
| 1108 |
walk[w1][we1] = trilen; |
| 1109 |
walk[trilen][2] = trilen+1; |
| 1110 |
trilen++; |
| 1111 |
|
| 1112 |
tri[trilen][0] = p0;
|
| 1113 |
tri[trilen][1] = pt;
|
| 1114 |
tri[trilen][2] = p2;
|
| 1115 |
walk[trilen][0] = curtri;
|
| 1116 |
walk[trilen][1] = trilen-1; |
| 1117 |
walk[trilen][2] = w2;
|
| 1118 |
walk[w2][we2] = trilen; |
| 1119 |
trilen++; |
| 1120 |
} |
| 1121 |
} |
| 1122 |
break;
|
| 1123 |
} |
| 1124 |
|
| 1125 |
} |
| 1126 |
} |
| 1127 |
|
| 1128 |
// *** PHASE 3: *** triangulate connection between grids
|
| 1129 |
|
| 1130 |
if (curgrid == 0) { |
| 1131 |
// construct ltrinum, ltriedge, rtrinum, and rtriedge
|
| 1132 |
for (int i=0; i<ledges-1; i++) { |
| 1133 |
ltrinum[i] = tlr[i*(lenx-1)][1]; |
| 1134 |
ltriedge[i] = 2;
|
| 1135 |
} |
| 1136 |
for (int i=0; i<redges-1; i++) { |
| 1137 |
rtrinum[i] = tlr[(i+1)*(lenx-1)-1][2]; |
| 1138 |
rtriedge[i] = 0;
|
| 1139 |
} |
| 1140 |
} |
| 1141 |
else {
|
| 1142 |
// starting and ending indices of grid connection
|
| 1143 |
int start1; // index into samples (left edge of curgrid) |
| 1144 |
int start2; // index into leftedge |
| 1145 |
int end1; // index into samples (right edge of curgrid) |
| 1146 |
int end2; // index into rightedge |
| 1147 |
|
| 1148 |
// find start1 and start2
|
| 1149 |
int tmup = curgrid*lenp + lenx*(broken[curgrid][0]+1); |
| 1150 |
double tx = samples[0][tmup]; |
| 1151 |
double ty = samples[1][tmup]; |
| 1152 |
boolean sig = samples[1][leftedge[0]] |
| 1153 |
< samples[1][leftedge[ledges-1]]; |
| 1154 |
|
| 1155 |
if (samples[1][leftedge[ledges-1]] < ty != sig |
| 1156 |
&& samples[1][leftedge[ledges-1]] != ty) { |
| 1157 |
// use a binary search to find prospective start2
|
| 1158 |
int min = 0; |
| 1159 |
int max = ledges-1; |
| 1160 |
int sg = (min+max)/2; |
| 1161 |
int itnum = 0; |
| 1162 |
while (samples[1][leftedge[sg+1]] < ty |
| 1163 |
== samples[1][leftedge[sg]] < ty && itnum < ledges) {
|
| 1164 |
if (samples[1][leftedge[sg]] < ty == sig) { |
| 1165 |
// point is closer to ledges-1
|
| 1166 |
min = sg; |
| 1167 |
} |
| 1168 |
else {
|
| 1169 |
// point is closer to 0
|
| 1170 |
if (sg == 0) { |
| 1171 |
throw new SetException("DelaunayOverlap: " |
| 1172 |
+"illegal grid overlap");
|
| 1173 |
} |
| 1174 |
max = sg; |
| 1175 |
} |
| 1176 |
sg = (min+max)/2;
|
| 1177 |
if (sg == ledges-1) { |
| 1178 |
throw new SetException("DelaunayOverlap: " |
| 1179 |
+"pathological grid overlap");
|
| 1180 |
} |
| 1181 |
itnum++; |
| 1182 |
} |
| 1183 |
if (itnum >= ledges) {
|
| 1184 |
throw new SetException("DelaunayOverlap: corrupt " |
| 1185 |
+"left edge structure");
|
| 1186 |
} |
| 1187 |
int clep = leftedge[sg];
|
| 1188 |
double cx = samples[0][clep]; |
| 1189 |
double cy = samples[1][clep]; |
| 1190 |
int clep1 = leftedge[sg+1]; |
| 1191 |
double c1x = samples[0][clep1]; |
| 1192 |
double c1y = samples[1][clep1]; |
| 1193 |
|
| 1194 |
// perform cross-product test to verify that points are correct
|
| 1195 |
int pt = curgrid*lenp - lenx + 1; |
| 1196 |
double px = samples[0][pt]; |
| 1197 |
double py = samples[1][pt]; |
| 1198 |
|
| 1199 |
if ( (tx-cx)*(ty-c1y) - (ty-cy)*(tx-c1x) > 0 |
| 1200 |
== (px-cx)*(py-c1y) - (py-cy)*(px-c1x) > 0 ) {
|
| 1201 |
// inverted overlap, must adopt slightly different strategy
|
| 1202 |
start2 = ledges-1;
|
| 1203 |
|
| 1204 |
// use a binary search to find start1
|
| 1205 |
min = 0;
|
| 1206 |
max = leny-1;
|
| 1207 |
sg = (min+max)/2;
|
| 1208 |
itnum = 0;
|
| 1209 |
double ll = samples[1][leftedge[ledges-1]]; |
| 1210 |
int offst = curgrid*lenp;
|
| 1211 |
int offst2 = offst+lenx;
|
| 1212 |
while (samples[1][lenx*sg+offst] < ll |
| 1213 |
== samples[1][lenx*sg+offst2] < ll && itnum < leny) {
|
| 1214 |
if (ll > samples[1][lenx*sg+offst]) { |
| 1215 |
// point is closer to leny-1
|
| 1216 |
min = sg; |
| 1217 |
} |
| 1218 |
else {
|
| 1219 |
// point is closer to 0
|
| 1220 |
if (sg == 0) { |
| 1221 |
throw new SetException("DelaunayOverlap: " |
| 1222 |
+"illegal grid overlap");
|
| 1223 |
} |
| 1224 |
max = sg; |
| 1225 |
} |
| 1226 |
sg = (min+max)/2;
|
| 1227 |
if (sg == leny-1) { |
| 1228 |
throw new SetException("DelaunayOverlap: " |
| 1229 |
+"pathological grid overlap");
|
| 1230 |
} |
| 1231 |
itnum++; |
| 1232 |
} |
| 1233 |
if (itnum >= leny) {
|
| 1234 |
throw new SetException("DelaunayOverlap: corrupt " |
| 1235 |
+"grid structure");
|
| 1236 |
} |
| 1237 |
start1 = lenx*sg+offst2; |
| 1238 |
} |
| 1239 |
else {
|
| 1240 |
// normal overlap
|
| 1241 |
start1 = tmup; |
| 1242 |
start2 = sg; |
| 1243 |
} |
| 1244 |
} |
| 1245 |
else {
|
| 1246 |
// tmup is outside range of leftedge, use last leftedge pt
|
| 1247 |
start1 = tmup; |
| 1248 |
start2 = ledges-1;
|
| 1249 |
} |
| 1250 |
|
| 1251 |
// find end1 and end2
|
| 1252 |
tmup = curgrid*lenp + lenx*(broken[curgrid][lenx]+2) - 1; |
| 1253 |
tx = samples[0][tmup];
|
| 1254 |
ty = samples[1][tmup];
|
| 1255 |
sig = samples[1][rightedge[0]] |
| 1256 |
< samples[1][rightedge[redges-1]]; |
| 1257 |
|
| 1258 |
if (samples[1][rightedge[redges-1]] < ty != sig |
| 1259 |
&& samples[1][rightedge[redges-1]] != ty) { |
| 1260 |
// use a binary search to find prospective end2
|
| 1261 |
int min = 0; |
| 1262 |
int max = redges-1; |
| 1263 |
int sg = (min+max)/2; |
| 1264 |
int itnum = 0; |
| 1265 |
while (samples[1][rightedge[sg+1]] < ty |
| 1266 |
== samples[1][rightedge[sg]] < ty && itnum < redges) {
|
| 1267 |
if (samples[1][rightedge[sg]] < ty == sig) { |
| 1268 |
// point is closer to redges-1
|
| 1269 |
min = sg; |
| 1270 |
} |
| 1271 |
else {
|
| 1272 |
// point is closer to 0
|
| 1273 |
if (sg == 0) { |
| 1274 |
throw new SetException("DelaunayOverlap: " |
| 1275 |
+"illegal grid overlap");
|
| 1276 |
} |
| 1277 |
max = sg; |
| 1278 |
} |
| 1279 |
sg = (min+max)/2;
|
| 1280 |
if (sg == redges-1) { |
| 1281 |
throw new SetException("DelaunayOverlap: " |
| 1282 |
+"pathological grid overlap");
|
| 1283 |
} |
| 1284 |
itnum++; |
| 1285 |
} |
| 1286 |
if (itnum >= redges) {
|
| 1287 |
throw new SetException("DelaunayOverlap: corrupt " |
| 1288 |
+"right edge structure");
|
| 1289 |
} |
| 1290 |
int crep = rightedge[sg];
|
| 1291 |
double cx = samples[0][crep]; |
| 1292 |
double cy = samples[1][crep]; |
| 1293 |
int crep1 = rightedge[sg+1]; |
| 1294 |
double c1x = samples[0][crep1]; |
| 1295 |
double c1y = samples[1][crep1]; |
| 1296 |
|
| 1297 |
// perform cross-product test to verify that points are correct
|
| 1298 |
int pt = curgrid*lenp - 2; |
| 1299 |
double px = samples[0][pt]; |
| 1300 |
double py = samples[1][pt]; |
| 1301 |
|
| 1302 |
if ( (tx-cx)*(ty-c1y) - (ty-cy)*(tx-c1x) > 0 |
| 1303 |
== (px-cx)*(py-c1y) - (py-cy)*(px-c1x) > 0 ) {
|
| 1304 |
// inverted overlap, must adopt slightly different strategy
|
| 1305 |
end2 = redges-1;
|
| 1306 |
|
| 1307 |
// use a binary search to find end1
|
| 1308 |
min = 0;
|
| 1309 |
max = leny-1;
|
| 1310 |
sg = (min+max)/2;
|
| 1311 |
itnum = 0;
|
| 1312 |
double ll = samples[1][rightedge[redges-1]]; |
| 1313 |
int offst = curgrid*lenp+lenx-1; |
| 1314 |
int offst2 = offst+lenx;
|
| 1315 |
while (samples[1][lenx*sg+offst] < ll |
| 1316 |
== samples[1][lenx*sg+offst2] < ll && itnum < leny) {
|
| 1317 |
if (ll > samples[1][lenx*sg+offst]) { |
| 1318 |
// point is closer to leny-1
|
| 1319 |
min = sg; |
| 1320 |
} |
| 1321 |
else {
|
| 1322 |
// point is closer to 0
|
| 1323 |
if (sg == 0) { |
| 1324 |
throw new SetException("DelaunayOverlap: " |
| 1325 |
+"illegal grid overlap");
|
| 1326 |
} |
| 1327 |
max = sg; |
| 1328 |
} |
| 1329 |
sg = (min+max)/2;
|
| 1330 |
if (sg == leny-1) { |
| 1331 |
throw new SetException("DelaunayOverlap: " |
| 1332 |
+"pathological grid overlap");
|
| 1333 |
} |
| 1334 |
itnum++; |
| 1335 |
} |
| 1336 |
if (itnum >= leny) {
|
| 1337 |
throw new SetException("DelaunayOverlap: corrupt " |
| 1338 |
+"grid structure");
|
| 1339 |
} |
| 1340 |
end1 = lenx*sg+offst2; |
| 1341 |
} |
| 1342 |
else {
|
| 1343 |
// normal overlap
|
| 1344 |
end1 = tmup; |
| 1345 |
end2 = sg; |
| 1346 |
} |
| 1347 |
} |
| 1348 |
else {
|
| 1349 |
// tmup is outside range of rightedge, use last rightedge pt
|
| 1350 |
end1 = tmup; |
| 1351 |
end2 = redges-1;
|
| 1352 |
} |
| 1353 |
|
| 1354 |
// save old edge values
|
| 1355 |
old_leftedge = leftedge; |
| 1356 |
old_rightedge = rightedge; |
| 1357 |
old_ledges = ledges; |
| 1358 |
old_redges = redges; |
| 1359 |
old_ltrinum = ltrinum; |
| 1360 |
old_ltriedge = ltriedge; |
| 1361 |
old_rtrinum = rtrinum; |
| 1362 |
old_rtriedge = rtriedge; |
| 1363 |
|
| 1364 |
// update leftedge and rightedge arrays
|
| 1365 |
leftedge = new int[numgrids*leny]; |
| 1366 |
ltrinum = new int[numgrids*leny]; |
| 1367 |
ltriedge = new int[numgrids*leny]; |
| 1368 |
ledges = 0;
|
| 1369 |
rightedge = new int[numgrids*leny]; |
| 1370 |
rtrinum = new int[numgrids*leny]; |
| 1371 |
rtriedge = new int[numgrids*leny]; |
| 1372 |
redges = 0;
|
| 1373 |
System.arraycopy(old_leftedge, 0, leftedge, 0, start2+1); |
| 1374 |
System.arraycopy(old_ltrinum, 0, ltrinum, 0, start2); |
| 1375 |
System.arraycopy(old_ltriedge, 0, ltriedge, 0, start2); |
| 1376 |
ledges += start2+1;
|
| 1377 |
System.arraycopy(old_rightedge, 0, rightedge, 0, end2+1); |
| 1378 |
System.arraycopy(old_rtrinum, 0, rtrinum, 0, end2); |
| 1379 |
System.arraycopy(old_rtriedge, 0, rtriedge, 0, end2); |
| 1380 |
redges += end2+1;
|
| 1381 |
for (int i=start1; i<=(curgrid+1)*lenp-lenx; i+=lenx) { |
| 1382 |
leftedge[ledges++] = i; |
| 1383 |
} |
| 1384 |
for (int i=end1; i<(curgrid+1)*lenp; i+=lenx) { |
| 1385 |
rightedge[redges++] = i; |
| 1386 |
} |
| 1387 |
int curledge = start2;
|
| 1388 |
int curredge = redges - leny + broken[curgrid][lenx-1]; |
| 1389 |
|
| 1390 |
// indices into current quadrilateral's points
|
| 1391 |
int base1, oneUp1;
|
| 1392 |
int base2, oneUp2;
|
| 1393 |
|
| 1394 |
// x and y coordinates of base1 and oneUp1
|
| 1395 |
int base1x, base1y, oneUp1x, oneUp1y;
|
| 1396 |
|
| 1397 |
// flag indicating which array base2 & oneUp2 reference
|
| 1398 |
// b2lbr = 0 --> base2 & oneUp2 reference old_leftedge array
|
| 1399 |
// b2lbr = 1 --> base2 & oneUp2 reference samples array
|
| 1400 |
// (bottom of curgrid-1)
|
| 1401 |
// b2lbr = 2 --> base2 & oneUp2 reference old_rightedge array
|
| 1402 |
int b2lbr = 0; |
| 1403 |
|
| 1404 |
// initialize base1 and base2
|
| 1405 |
base1 = start1; |
| 1406 |
base2 = start2; |
| 1407 |
if (base2 == old_ledges-1) { |
| 1408 |
base2 = curgrid*lenp - lenx; |
| 1409 |
b2lbr = 1;
|
| 1410 |
} |
| 1411 |
base1x = 0;
|
| 1412 |
base1y = base1 % lenp / lenx; |
| 1413 |
boolean down = false; |
| 1414 |
|
| 1415 |
// initialize oneUp1 and oneUp2
|
| 1416 |
if (base1y > 0 && broken[curgrid][0] < base1y - 1) { |
| 1417 |
// move up
|
| 1418 |
oneUp1x = 0;
|
| 1419 |
oneUp1y = base1y - 1;
|
| 1420 |
} |
| 1421 |
else if (base1y == leny-1 || broken[curgrid][1] < base1y) { |
| 1422 |
// move right
|
| 1423 |
oneUp1x = 1;
|
| 1424 |
oneUp1y = base1y; |
| 1425 |
} |
| 1426 |
else {
|
| 1427 |
// move down
|
| 1428 |
oneUp1x = 0;
|
| 1429 |
oneUp1y = base1y + 1;
|
| 1430 |
down = true;
|
| 1431 |
} |
| 1432 |
oneUp1 = curgrid*lenp + oneUp1y*lenx + oneUp1x; |
| 1433 |
oneUp2 = base2 + 1;
|
| 1434 |
|
| 1435 |
boolean diag = false; |
| 1436 |
boolean firsttri = true; |
| 1437 |
while (base1 != end1 || base2 != end2 || b2lbr != 2) { |
| 1438 |
|
| 1439 |
// determine which diagonal to take
|
| 1440 |
if (base1 == end1) diag = true; |
| 1441 |
else if (base2 == end2 && b2lbr == 2) diag = false; |
| 1442 |
else {
|
| 1443 |
// set up variables
|
| 1444 |
double ax = samples[0][base1]; |
| 1445 |
double ay = samples[1][base1]; |
| 1446 |
double cx = samples[0][oneUp1]; |
| 1447 |
double cy = samples[1][oneUp1]; |
| 1448 |
double bx, by, dx, dy;
|
| 1449 |
if (b2lbr == 0) { |
| 1450 |
bx = samples[0][old_leftedge[base2]];
|
| 1451 |
by = samples[1][old_leftedge[base2]];
|
| 1452 |
dx = samples[0][old_leftedge[oneUp2]];
|
| 1453 |
dy = samples[1][old_leftedge[oneUp2]];
|
| 1454 |
} |
| 1455 |
else if (b2lbr == 1) { |
| 1456 |
bx = samples[0][base2];
|
| 1457 |
by = samples[1][base2];
|
| 1458 |
dx = samples[0][oneUp2];
|
| 1459 |
dy = samples[1][oneUp2];
|
| 1460 |
} |
| 1461 |
else {
|
| 1462 |
bx = samples[0][old_rightedge[base2]];
|
| 1463 |
by = samples[1][old_rightedge[base2]];
|
| 1464 |
dx = samples[0][old_rightedge[oneUp2]];
|
| 1465 |
dy = samples[1][old_rightedge[oneUp2]];
|
| 1466 |
} |
| 1467 |
|
| 1468 |
// perform diagonal test
|
| 1469 |
double abx = ax - bx;
|
| 1470 |
double aby = ay - by;
|
| 1471 |
double acx = ax - cx;
|
| 1472 |
double acy = ay - cy;
|
| 1473 |
double dbx = dx - bx;
|
| 1474 |
double dby = dy - by;
|
| 1475 |
double dcx = dx - cx;
|
| 1476 |
double dcy = dy - cy;
|
| 1477 |
double Q = abx*acx + aby*acy;
|
| 1478 |
double R = dbx*abx + dby*aby;
|
| 1479 |
double S = acx*dcx + acy*dcy;
|
| 1480 |
double T = dbx*dcx + dby*dcy;
|
| 1481 |
boolean QD = abx*acy - aby*acx > 0; |
| 1482 |
boolean RD = dbx*aby - dby*abx > 0; |
| 1483 |
boolean SD = acx*dcy - acy*dcx > 0; |
| 1484 |
boolean TD = dcx*dby - dcy*dbx > 0; |
| 1485 |
boolean sigD = (QD ? 1 : 0) + (RD ? 1 : 0) |
| 1486 |
+ (SD ? 1 : 0) + (TD ? 1 : 0) < 2; |
| 1487 |
if (QD == sigD) diag = true; |
| 1488 |
else if (RD == sigD || SD == sigD) diag = false; |
| 1489 |
else if (TD == sigD) diag = true; |
| 1490 |
else if (Q < 0 && T < 0 || R > 0 && S > 0) diag = true; |
| 1491 |
else if (R < 0 && S < 0 || Q > 0 && T > 0) diag = false; |
| 1492 |
else if ((Q < 0 ? Q : T) < (R < 0 ? R : S)) diag = true; |
| 1493 |
else diag = false; |
| 1494 |
} |
| 1495 |
|
| 1496 |
// walk[*][0] --> NEXT
|
| 1497 |
// walk[*][diag ? 1 : 2] --> PREVIOUS
|
| 1498 |
// walk[*][diag ? 2 : 1] --> OUTSIDE
|
| 1499 |
|
| 1500 |
if (diag) { // diagonal goes from base1 to oneUp2 |
| 1501 |
|
| 1502 |
// update tri array
|
| 1503 |
if (b2lbr == 0) { |
| 1504 |
tri[trilen][0] = old_leftedge[oneUp2];
|
| 1505 |
tri[trilen][1] = base1;
|
| 1506 |
tri[trilen][2] = old_leftedge[base2];
|
| 1507 |
} |
| 1508 |
else if (b2lbr == 1) { |
| 1509 |
tri[trilen][0] = oneUp2;
|
| 1510 |
tri[trilen][1] = base1;
|
| 1511 |
tri[trilen][2] = base2;
|
| 1512 |
} |
| 1513 |
else {
|
| 1514 |
tri[trilen][0] = old_rightedge[oneUp2];
|
| 1515 |
tri[trilen][1] = base1;
|
| 1516 |
tri[trilen][2] = old_rightedge[base2];
|
| 1517 |
} |
| 1518 |
|
| 1519 |
// update walk array
|
| 1520 |
|
| 1521 |
// current --> next
|
| 1522 |
walk[trilen][0] = -1; |
| 1523 |
|
| 1524 |
// current <--> previous
|
| 1525 |
if (firsttri) {
|
| 1526 |
ltrinum[curledge] = trilen; |
| 1527 |
ltriedge[curledge] = 2;
|
| 1528 |
curledge++; |
| 1529 |
firsttri = false;
|
| 1530 |
} |
| 1531 |
else {
|
| 1532 |
walk[trilen][1] = trilen-1; |
| 1533 |
walk[trilen-1][0] = trilen; |
| 1534 |
} |
| 1535 |
|
| 1536 |
// current <--> outside
|
| 1537 |
if (b2lbr == 0) { |
| 1538 |
int x = old_ltrinum[base2];
|
| 1539 |
walk[trilen][2] = x;
|
| 1540 |
walk[x][old_ltriedge[base2]] = trilen; |
| 1541 |
} |
| 1542 |
else if (b2lbr == 1) { |
| 1543 |
int x = old_bottom[base2 % lenx];
|
| 1544 |
walk[trilen][2] = x;
|
| 1545 |
walk[x][1] = trilen;
|
| 1546 |
} |
| 1547 |
else {
|
| 1548 |
int x = old_rtrinum[oneUp2];
|
| 1549 |
walk[trilen][2] = x;
|
| 1550 |
walk[x][old_rtriedge[oneUp2]] = trilen; |
| 1551 |
} |
| 1552 |
|
| 1553 |
// update base2
|
| 1554 |
base2 = oneUp2; |
| 1555 |
if (b2lbr == 0 && base2 == old_ledges-1) { |
| 1556 |
b2lbr = 1;
|
| 1557 |
base2 = curgrid*lenp - lenx; |
| 1558 |
} |
| 1559 |
if (b2lbr == 1 && base2 == curgrid*lenp - 1) { |
| 1560 |
b2lbr = 2;
|
| 1561 |
base2 = old_redges-1;
|
| 1562 |
} |
| 1563 |
|
| 1564 |
// update oneUp2
|
| 1565 |
oneUp2 = base2 + (b2lbr == 2 ? -1 : 1); |
| 1566 |
} |
| 1567 |
else { // diagonal goes from base2 to oneUp1 |
| 1568 |
|
| 1569 |
// update tri array
|
| 1570 |
if (b2lbr == 0) { |
| 1571 |
tri[trilen][0] = old_leftedge[base2];
|
| 1572 |
tri[trilen][1] = oneUp1;
|
| 1573 |
tri[trilen][2] = base1;
|
| 1574 |
} |
| 1575 |
else if (b2lbr == 1) { |
| 1576 |
tri[trilen][0] = base2;
|
| 1577 |
tri[trilen][1] = oneUp1;
|
| 1578 |
tri[trilen][2] = base1;
|
| 1579 |
} |
| 1580 |
else {
|
| 1581 |
tri[trilen][0] = old_rightedge[base2];
|
| 1582 |
tri[trilen][1] = oneUp1;
|
| 1583 |
tri[trilen][2] = base1;
|
| 1584 |
} |
| 1585 |
|
| 1586 |
// update walk array
|
| 1587 |
|
| 1588 |
// current --> next
|
| 1589 |
walk[trilen][0] = -1; |
| 1590 |
|
| 1591 |
// current <--> previous
|
| 1592 |
if (firsttri) {
|
| 1593 |
ltrinum[curledge] = trilen; |
| 1594 |
ltriedge[curledge] = 1;
|
| 1595 |
curledge++; |
| 1596 |
firsttri = false;
|
| 1597 |
} |
| 1598 |
else {
|
| 1599 |
walk[trilen][2] = trilen-1; |
| 1600 |
walk[trilen-1][0] = trilen; |
| 1601 |
} |
| 1602 |
|
| 1603 |
// current <--> outside
|
| 1604 |
if (oneUp1 - base1 == -lenx) {
|
| 1605 |
// base1 -> oneUp1 = up
|
| 1606 |
if (base1x < lenx-1) { |
| 1607 |
int x = tlr[(lenx-1)*oneUp1y+oneUp1x][1]; |
| 1608 |
walk[trilen][1] = x;
|
| 1609 |
walk[x][2] = trilen;
|
| 1610 |
} |
| 1611 |
else {
|
| 1612 |
// base1 is on the right edge of current grid
|
| 1613 |
walk[trilen][1] = -1; |
| 1614 |
|
| 1615 |
// update right edge triangle structure
|
| 1616 |
rtrinum[curredge] = trilen; |
| 1617 |
rtriedge[curredge] = 1;
|
| 1618 |
curredge--; |
| 1619 |
} |
| 1620 |
} |
| 1621 |
else if (oneUp1 - base1 == 1) { |
| 1622 |
// base1 -> oneUp1 = right
|
| 1623 |
if (base1y < leny-1) { |
| 1624 |
int inx = (lenx-1)*base1y+base1x; |
| 1625 |
int x = tlr[inx][0]; |
| 1626 |
walk[trilen][1] = x;
|
| 1627 |
walk[x][istwo[inx] ? 2 : 0] = trilen; |
| 1628 |
} |
| 1629 |
else {
|
| 1630 |
// base1 is on the bottom edge of current grid
|
| 1631 |
walk[trilen][1] = -1; |
| 1632 |
bottom[base1x] = trilen; |
| 1633 |
} |
| 1634 |
} |
| 1635 |
else {
|
| 1636 |
// base1 -> oneUp1 = down
|
| 1637 |
if (base1x > 0) { |
| 1638 |
int x = tlr[(lenx-1)*base1y+base1x-1][2]; |
| 1639 |
walk[trilen][1] = x;
|
| 1640 |
walk[x][0] = trilen;
|
| 1641 |
} |
| 1642 |
else {
|
| 1643 |
// base1 is on the left edge of current grid
|
| 1644 |
walk[trilen][1] = -1; |
| 1645 |
ltrinum[curledge] = trilen; |
| 1646 |
ltriedge[curledge] = 1;
|
| 1647 |
curledge++; |
| 1648 |
} |
| 1649 |
} |
| 1650 |
|
| 1651 |
// update base1
|
| 1652 |
base1 = oneUp1; |
| 1653 |
base1x = oneUp1x; |
| 1654 |
base1y = oneUp1y; |
| 1655 |
|
| 1656 |
// update oneUp1
|
| 1657 |
if (broken[curgrid][base1x == 0 ? 0 : base1x+1] < base1y - 1 |
| 1658 |
&& base1y > 0 && !down) {
|
| 1659 |
// move up
|
| 1660 |
oneUp1x = base1x; |
| 1661 |
oneUp1y = base1y - 1;
|
| 1662 |
} |
| 1663 |
else if ( base1y == leny-1 || (base1x < lenx-1 |
| 1664 |
&& broken[curgrid][base1x+1] < base1y) ) {
|
| 1665 |
// move right
|
| 1666 |
oneUp1x = base1x+1;
|
| 1667 |
oneUp1y = base1y; |
| 1668 |
down = false;
|
| 1669 |
} |
| 1670 |
else {
|
| 1671 |
// move down
|
| 1672 |
oneUp1x = base1x; |
| 1673 |
oneUp1y = base1y + 1;
|
| 1674 |
down = true;
|
| 1675 |
} |
| 1676 |
oneUp1 = curgrid*lenp + oneUp1y*lenx + oneUp1x; |
| 1677 |
} |
| 1678 |
|
| 1679 |
trilen++; |
| 1680 |
|
| 1681 |
// expand tri array if necessary
|
| 1682 |
if (trilen > trisize) {
|
| 1683 |
trisize += trisize; |
| 1684 |
int[][] newtri = new int[trisize+2][3]; |
| 1685 |
int[][] newwalk = new int[trisize+2][3]; |
| 1686 |
for (int i=0; i<trilen; i++) { |
| 1687 |
newtri[i][0] = tri[i][0]; |
| 1688 |
newtri[i][1] = tri[i][1]; |
| 1689 |
newtri[i][2] = tri[i][2]; |
| 1690 |
newwalk[i][0] = walk[i][0]; |
| 1691 |
newwalk[i][1] = walk[i][1]; |
| 1692 |
newwalk[i][2] = walk[i][2]; |
| 1693 |
} |
| 1694 |
for (int i=trilen; i<trisize+2; i++) { |
| 1695 |
newwalk[i][0] = -1; |
| 1696 |
newwalk[i][1] = -1; |
| 1697 |
newwalk[i][2] = -1; |
| 1698 |
} |
| 1699 |
tri = newtri; |
| 1700 |
walk = newwalk; |
| 1701 |
} |
| 1702 |
} |
| 1703 |
|
| 1704 |
// add triangle number (trilen-1) to rtrinum and rtriedge
|
| 1705 |
rtrinum[curredge] = trilen-1;
|
| 1706 |
rtriedge[curredge] = (diag ? 2 : 1); |
| 1707 |
curredge = redges - leny + broken[curgrid][lenx-1] + 1; |
| 1708 |
|
| 1709 |
// finish updating ltrinum, ltriedge, rtrinum, and rtriedge
|
| 1710 |
int x = broken[curgrid][1]+1; |
| 1711 |
for (int i=x; i<leny-1; i++) { |
| 1712 |
ltrinum[curledge] = tlr[(lenx-1)*i][1]; |
| 1713 |
ltriedge[curledge] = 2;
|
| 1714 |
curledge++; |
| 1715 |
} |
| 1716 |
x = broken[curgrid][lenx-1]+1; |
| 1717 |
for (int i=x; i<leny-1; i++) { |
| 1718 |
rtrinum[curredge] = tlr[(lenx-1)*(i+1)-1][2]; |
| 1719 |
rtriedge[curredge] = 0;
|
| 1720 |
curredge++; |
| 1721 |
} |
| 1722 |
} |
| 1723 |
} |
| 1724 |
|
| 1725 |
// *** PHASE 4: *** sub-triangulate lost grid points
|
| 1726 |
|
| 1727 |
for (curgrid=0; curgrid<numgrids; curgrid++) { |
| 1728 |
|
| 1729 |
// *** PHASE 4a: *** sub-triangulate Type 1 lost grid points
|
| 1730 |
|
| 1731 |
// find all Type 1 lost grid points and deal with them.
|
| 1732 |
// Type 1 lost grid points are those that fall into a broken box
|
| 1733 |
|
| 1734 |
// The final grid's boxes cannot contain any points, and
|
| 1735 |
// thus cannot contain any Type 1 lost points
|
| 1736 |
if (curgrid < numgrids-1) { |
| 1737 |
int curtri = trilen/2; // a reasonable starting guess |
| 1738 |
for (int gx=0; gx<lenx-1; gx++) { |
| 1739 |
for (int gy=0; gy<=broken[curgrid][gx+1]; gy++) { |
| 1740 |
|
| 1741 |
// subtriangulate all points within broken grid boxes
|
| 1742 |
int qg = leng*curgrid + (lenx-1)*gy + gx; |
| 1743 |
|
| 1744 |
if (curtri < 0) curtri = trilen/2; |
| 1745 |
for (int pt=0; pt<glen[qg]; pt++) { |
| 1746 |
|
| 1747 |
// point (px, py) is broken
|
| 1748 |
double Px = samples[0][grid[qg][pt]]; |
| 1749 |
double Py = samples[1][grid[qg][pt]]; |
| 1750 |
|
| 1751 |
// locate containing triangle of (px, py)
|
| 1752 |
boolean located = false; |
| 1753 |
int itnum;
|
| 1754 |
for (itnum=0; itnum<trilen && !located; itnum++) { |
| 1755 |
// define data
|
| 1756 |
int t0 = tri[curtri][0]; |
| 1757 |
int t1 = tri[curtri][1]; |
| 1758 |
int t2 = tri[curtri][2]; |
| 1759 |
double Ax = samples[0][t0]; |
| 1760 |
double Ay = samples[1][t0]; |
| 1761 |
double Bx = samples[0][t1]; |
| 1762 |
double By = samples[1][t1]; |
| 1763 |
double Cx = samples[0][t2]; |
| 1764 |
double Cy = samples[1][t2]; |
| 1765 |
|
| 1766 |
// test whether point is contained in current triangle
|
| 1767 |
double tval0 = (Bx-Ax)*(Py-Ay) - (By-Ay)*(Px-Ax);
|
| 1768 |
double tval1 = (Cx-Bx)*(Py-By) - (Cy-By)*(Px-Bx);
|
| 1769 |
double tval2 = (Ax-Cx)*(Py-Cy) - (Ay-Cy)*(Px-Cx);
|
| 1770 |
boolean test0 = (tval0 == 0) || ( (tval0 > 0) == ( |
| 1771 |
(Bx-Ax)*(Cy-Ay) - (By-Ay)*(Cx-Ax) > 0) );
|
| 1772 |
boolean test1 = (tval1 == 0) || ( (tval1 > 0) == ( |
| 1773 |
(Cx-Bx)*(Ay-By) - (Cy-By)*(Ax-Bx) > 0) );
|
| 1774 |
boolean test2 = (tval2 == 0) || ( (tval2 > 0) == ( |
| 1775 |
(Ax-Cx)*(By-Cy) - (Ay-Cy)*(Bx-Cx) > 0) );
|
| 1776 |
|
| 1777 |
// figure out which triangle to go to next
|
| 1778 |
if (!test0 && !test1 && !test2) {
|
| 1779 |
throw new SetException("DelaunayOverlap: corrupt " |
| 1780 |
+"triangle structure");
|
| 1781 |
} |
| 1782 |
if (!test0 && !test1) {
|
| 1783 |
int tri0 = walk[curtri][0]; |
| 1784 |
int tri1 = walk[curtri][1]; |
| 1785 |
if (tri0 == -1) curtri = tri1; |
| 1786 |
else curtri = tri0;
|
| 1787 |
} |
| 1788 |
else if (!test0 && !test2) { |
| 1789 |
int tri0 = walk[curtri][0]; |
| 1790 |
int tri2 = walk[curtri][2]; |
| 1791 |
if (tri0 == -1) curtri = tri2; |
| 1792 |
else curtri = tri0;
|
| 1793 |
} |
| 1794 |
else if (!test1 && !test2) { |
| 1795 |
int tri1 = walk[curtri][1]; |
| 1796 |
int tri2 = walk[curtri][2]; |
| 1797 |
if (tri1 == -1) curtri = tri2; |
| 1798 |
else curtri = tri1;
|
| 1799 |
} |
| 1800 |
else if (!test0) curtri = walk[curtri][0]; |
| 1801 |
else if (!test1) curtri = walk[curtri][1]; |
| 1802 |
else if (!test2) curtri = walk[curtri][2]; |
| 1803 |
else located = true; |
| 1804 |
|
| 1805 |
// point is outside current triangulation
|
| 1806 |
if (curtri < 0) itnum = trilen; |
| 1807 |
} |
| 1808 |
|
| 1809 |
// if itnum == trilen, the point is permanently lost
|
| 1810 |
if (itnum < trilen) {
|
| 1811 |
// subtriangulate point grid[qg][pt]
|
| 1812 |
// (curtri is the containing triangle)
|
| 1813 |
// get curtri's point information
|
| 1814 |
int q = grid[qg][pt];
|
| 1815 |
int ct0 = tri[curtri][0]; |
| 1816 |
int ct1 = tri[curtri][1]; |
| 1817 |
int ct2 = tri[curtri][2]; |
| 1818 |
|
| 1819 |
// get curtri's walk information
|
| 1820 |
int T0 = walk[curtri][0]; |
| 1821 |
int T1 = walk[curtri][1]; |
| 1822 |
int T2 = walk[curtri][2]; |
| 1823 |
int TE0, TE1, TE2;
|
| 1824 |
if (T0 == -1) TE0 = -1; |
| 1825 |
else if (walk[T0][0] == curtri) TE0 = 0; |
| 1826 |
else if (walk[T0][1] == curtri) TE0 = 1; |
| 1827 |
else if (walk[T0][2] == curtri) TE0 = 2; |
| 1828 |
else throw new SetException("DelaunayOverlap: corrupt " |
| 1829 |
+"walk structure");
|
| 1830 |
if (T1 == -1) TE1 = -1; |
| 1831 |
else if (walk[T1][0] == curtri) TE1 = 0; |
| 1832 |
else if (walk[T1][1] == curtri) TE1 = 1; |
| 1833 |
else if (walk[T1][2] == curtri) TE1 = 2; |
| 1834 |
else throw new SetException("DelaunayOverlap: corrupt " |
| 1835 |
+"walk structure");
|
| 1836 |
if (T2 == -1) TE2 = -1; |
| 1837 |
else if (walk[T2][0] == curtri) TE2 = 0; |
| 1838 |
else if (walk[T2][1] == curtri) TE2 = 1; |
| 1839 |
else if (walk[T2][2] == curtri) TE2 = 2; |
| 1840 |
else throw new SetException("DelaunayOverlap: corrupt " |
| 1841 |
+"walk structure");
|
| 1842 |
|
| 1843 |
// define three new triangles
|
| 1844 |
// the first new triangle overwrites the old triangle
|
| 1845 |
// tri[curtri][0] = ct0; <-- already true
|
| 1846 |
// tri[curtri][1] = ct1; <-- already true
|
| 1847 |
tri[curtri][2] = q;
|
| 1848 |
// walk[curtri][0] = T0; <-- already true
|
| 1849 |
// if (TE0 >= 0) walk[T0][TE0] = curtri; <-- already true
|
| 1850 |
walk[curtri][1] = trilen;
|
| 1851 |
walk[curtri][2] = trilen+1; |
| 1852 |
|
| 1853 |
tri[trilen][0] = q;
|
| 1854 |
tri[trilen][1] = ct1;
|
| 1855 |
tri[trilen][2] = ct2;
|
| 1856 |
walk[trilen][0] = curtri;
|
| 1857 |
walk[trilen][1] = T1;
|
| 1858 |
if (TE1 >= 0) walk[T1][TE1] = trilen; |
| 1859 |
walk[trilen][2] = trilen+1; |
| 1860 |
trilen++; |
| 1861 |
|
| 1862 |
tri[trilen][0] = ct0;
|
| 1863 |
tri[trilen][1] = q;
|
| 1864 |
tri[trilen][2] = ct2;
|
| 1865 |
walk[trilen][0] = curtri;
|
| 1866 |
walk[trilen][1] = trilen-1; |
| 1867 |
walk[trilen][2] = T2;
|
| 1868 |
if (TE2 >= 0) walk[T2][TE2] = trilen; |
| 1869 |
trilen++; |
| 1870 |
|
| 1871 |
// expand tri array if necessary
|
| 1872 |
if (trilen > trisize) {
|
| 1873 |
trisize += trisize; |
| 1874 |
int[][] newtri = new int[trisize+2][3]; |
| 1875 |
int[][] newwalk = new int[trisize+2][3]; |
| 1876 |
for (int i=0; i<trilen; i++) { |
| 1877 |
newtri[i][0] = tri[i][0]; |
| 1878 |
newtri[i][1] = tri[i][1]; |
| 1879 |
newtri[i][2] = tri[i][2]; |
| 1880 |
newwalk[i][0] = walk[i][0]; |
| 1881 |
newwalk[i][1] = walk[i][1]; |
| 1882 |
newwalk[i][2] = walk[i][2]; |
| 1883 |
} |
| 1884 |
for (int i=trilen; i<trisize+2; i++) { |
| 1885 |
newwalk[i][0] = -1; |
| 1886 |
newwalk[i][1] = -1; |
| 1887 |
newwalk[i][2] = -1; |
| 1888 |
} |
| 1889 |
tri = newtri; |
| 1890 |
walk = newwalk; |
| 1891 |
} |
| 1892 |
|
| 1893 |
} |
| 1894 |
} |
| 1895 |
} |
| 1896 |
} |
| 1897 |
} |
| 1898 |
|
| 1899 |
// *** PHASE 4b: *** sub-triangulate Type 2 lost grid points
|
| 1900 |
|
| 1901 |
// find all Type 2 lost grid points and deal with them.
|
| 1902 |
// Type 2 lost grid points are those that did not fall into any
|
| 1903 |
// box, but which have both grid boxes below them broken.
|
| 1904 |
int curtri = trilen/2; |
| 1905 |
for (int line=0; line<leny-1; line++) { |
| 1906 |
for (int pix=1; pix<lenx-1; pix++) { |
| 1907 |
|
| 1908 |
// find out if point is lost
|
| 1909 |
int q = lenp*curgrid + lenx*line + pix;
|
| 1910 |
if (!ptfell[q] && broken[curgrid][pix] >= line
|
| 1911 |
&& broken[curgrid][pix+1] >= line) {
|
| 1912 |
|
| 1913 |
// point (px, py) is lost
|
| 1914 |
double Px = samples[0][q]; |
| 1915 |
double Py = samples[1][q]; |
| 1916 |
|
| 1917 |
// locate containing triangle of (px, py)
|
| 1918 |
boolean located = false; |
| 1919 |
int itnum;
|
| 1920 |
for (itnum=0; itnum<trilen && !located; itnum++) { |
| 1921 |
// define data
|
| 1922 |
int t0 = tri[curtri][0]; |
| 1923 |
int t1 = tri[curtri][1]; |
| 1924 |
int t2 = tri[curtri][2]; |
| 1925 |
double Ax = samples[0][t0]; |
| 1926 |
double Ay = samples[1][t0]; |
| 1927 |
double Bx = samples[0][t1]; |
| 1928 |
double By = samples[1][t1]; |
| 1929 |
double Cx = samples[0][t2]; |
| 1930 |
double Cy = samples[1][t2]; |
| 1931 |
|
| 1932 |
// test whether point is contained in current triangle
|
| 1933 |
double tval0 = (Bx-Ax)*(Py-Ay) - (By-Ay)*(Px-Ax);
|
| 1934 |
double tval1 = (Cx-Bx)*(Py-By) - (Cy-By)*(Px-Bx);
|
| 1935 |
double tval2 = (Ax-Cx)*(Py-Cy) - (Ay-Cy)*(Px-Cx);
|
| 1936 |
boolean test0 = (tval0 == 0) || ( (tval0 > 0) == ( |
| 1937 |
(Bx-Ax)*(Cy-Ay) - (By-Ay)*(Cx-Ax) > 0) );
|
| 1938 |
boolean test1 = (tval1 == 0) || ( (tval1 > 0) == ( |
| 1939 |
(Cx-Bx)*(Ay-By) - (Cy-By)*(Ax-Bx) > 0) );
|
| 1940 |
boolean test2 = (tval2 == 0) || ( (tval2 > 0) == ( |
| 1941 |
(Ax-Cx)*(By-Cy) - (Ay-Cy)*(Bx-Cx) > 0) );
|
| 1942 |
|
| 1943 |
// figure out which triangle to go to next
|
| 1944 |
if (!test0 && !test1 && !test2) {
|
| 1945 |
throw new SetException("DelaunayOverlap: corrupt " |
| 1946 |
+"triangle structure");
|
| 1947 |
} |
| 1948 |
if (!test0 && !test1) {
|
| 1949 |
int tri0 = walk[curtri][0]; |
| 1950 |
int tri1 = walk[curtri][1]; |
| 1951 |
if (tri0 == -1) curtri = tri1; |
| 1952 |
else curtri = tri0;
|
| 1953 |
} |
| 1954 |
else if (!test0 && !test2) { |
| 1955 |
int tri0 = walk[curtri][0]; |
| 1956 |
int tri2 = walk[curtri][2]; |
| 1957 |
if (tri0 == -1) curtri = tri2; |
| 1958 |
else curtri = tri0;
|
| 1959 |
} |
| 1960 |
else if (!test1 && !test2) { |
| 1961 |
int tri1 = walk[curtri][1]; |
| 1962 |
int tri2 = walk[curtri][2]; |
| 1963 |
if (tri1 == -1) curtri = tri2; |
| 1964 |
else curtri = tri1;
|
| 1965 |
} |
| 1966 |
else if (!test0) curtri = walk[curtri][0]; |
| 1967 |
else if (!test1) curtri = walk[curtri][1]; |
| 1968 |
else if (!test2) curtri = walk[curtri][2]; |
| 1969 |
else located = true; |
| 1970 |
|
| 1971 |
// point is outside current triangulation
|
| 1972 |
// (this point is permanently lost)
|
| 1973 |
if (curtri < 0) itnum = trilen; |
| 1974 |
} |
| 1975 |
|
| 1976 |
// if itnum == trilen, the point is permanently lost
|
| 1977 |
if (itnum < trilen) {
|
| 1978 |
// subtriangulate point q
|
| 1979 |
// (curtri is the containing triangle)
|
| 1980 |
// get curtri's point information
|
| 1981 |
int ct0 = tri[curtri][0]; |
| 1982 |
int ct1 = tri[curtri][1]; |
| 1983 |
int ct2 = tri[curtri][2]; |
| 1984 |
|
| 1985 |
// get curtri's walk information
|
| 1986 |
int T0 = walk[curtri][0]; |
| 1987 |
int T1 = walk[curtri][1]; |
| 1988 |
int T2 = walk[curtri][2]; |
| 1989 |
int TE0, TE1, TE2;
|
| 1990 |
if (T0 == -1) TE0 = -1; |
| 1991 |
else if (walk[T0][0] == curtri) TE0 = 0; |
| 1992 |
else if (walk[T0][1] == curtri) TE0 = 1; |
| 1993 |
else if (walk[T0][2] == curtri) TE0 = 2; |
| 1994 |
else throw new SetException("DelaunayOverlap: corrupt " |
| 1995 |
+"walk structure");
|
| 1996 |
if (T1 == -1) TE1 = -1; |
| 1997 |
else if (walk[T1][0] == curtri) TE1 = 0; |
| 1998 |
else if (walk[T1][1] == curtri) TE1 = 1; |
| 1999 |
else if (walk[T1][2] == curtri) TE1 = 2; |
| 2000 |
else throw new SetException("DelaunayOverlap: corrupt " |
| 2001 |
+"walk structure");
|
| 2002 |
if (T2 == -1) TE2 = -1; |
| 2003 |
else if (walk[T2][0] == curtri) TE2 = 0; |
| 2004 |
else if (walk[T2][1] == curtri) TE2 = 1; |
| 2005 |
else if (walk[T2][2] == curtri) TE2 = 2; |
| 2006 |
else throw new SetException("DelaunayOverlap: corrupt " |
| 2007 |
+"walk structure");
|
| 2008 |
|
| 2009 |
// define three new triangles
|
| 2010 |
// the first new triangle overwrites the old triangle
|
| 2011 |
// tri[curtri][0] = ct0; <-- already true
|
| 2012 |
// tri[curtri][1] = ct1; <-- already true
|
| 2013 |
tri[curtri][2] = q;
|
| 2014 |
// walk[curtri][0] = T0; <-- already true
|
| 2015 |
// if (TE0 >= 0) walk[T0][TE0] = curtri; <-- already true
|
| 2016 |
walk[curtri][1] = trilen;
|
| 2017 |
walk[curtri][2] = trilen+1; |
| 2018 |
|
| 2019 |
tri[trilen][0] = q;
|
| 2020 |
tri[trilen][1] = ct1;
|
| 2021 |
tri[trilen][2] = ct2;
|
| 2022 |
walk[trilen][0] = curtri;
|
| 2023 |
walk[trilen][1] = T1;
|
| 2024 |
if (TE1 >= 0) walk[T1][TE1] = trilen; |
| 2025 |
walk[trilen][2] = trilen+1; |
| 2026 |
trilen++; |
| 2027 |
|
| 2028 |
tri[trilen][0] = ct0;
|
| 2029 |
tri[trilen][1] = q;
|
| 2030 |
tri[trilen][2] = ct2;
|
| 2031 |
walk[trilen][0] = curtri;
|
| 2032 |
walk[trilen][1] = trilen-1; |
| 2033 |
walk[trilen][2] = T2;
|
| 2034 |
if (TE2 >= 0) walk[T2][TE2] = trilen; |
| 2035 |
trilen++; |
| 2036 |
|
| 2037 |
// expand tri array if necessary
|
| 2038 |
if (trilen > trisize) {
|
| 2039 |
trisize += trisize; |
| 2040 |
int[][] newtri = new int[trisize+2][3]; |
| 2041 |
int[][] newwalk = new int[trisize+2][3]; |
| 2042 |
for (int i=0; i<trilen; i++) { |
| 2043 |
newtri[i][0] = tri[i][0]; |
| 2044 |
newtri[i][1] = tri[i][1]; |
| 2045 |
newtri[i][2] = tri[i][2]; |
| 2046 |
newwalk[i][0] = walk[i][0]; |
| 2047 |
newwalk[i][1] = walk[i][1]; |
| 2048 |
newwalk[i][2] = walk[i][2]; |
| 2049 |
} |
| 2050 |
for (int i=trilen; i<trisize+2; i++) { |
| 2051 |
newwalk[i][0] = -1; |
| 2052 |
newwalk[i][1] = -1; |
| 2053 |
newwalk[i][2] = -1; |
| 2054 |
} |
| 2055 |
tri = newtri; |
| 2056 |
walk = newwalk; |
| 2057 |
} |
| 2058 |
|
| 2059 |
} |
| 2060 |
} |
| 2061 |
} |
| 2062 |
} |
| 2063 |
|
| 2064 |
// There is a third type of lost grid point. If a point on
|
| 2065 |
// the left or right edge of a grid is cut off from the rest
|
| 2066 |
// of its grid by a zigzagging grid edge from a previous grid,
|
| 2067 |
// that point is Type 3 lost and will not be incorporated into
|
| 2068 |
// the triangulation. This case is not yet handled by
|
| 2069 |
// DelaunayOverlap's algorithm, but will be added to the
|
| 2070 |
// constructor if data sets exhibit properties which cause
|
| 2071 |
// the loss of many points.
|
| 2072 |
} |
| 2073 |
|
| 2074 |
// *** PHASE 5: *** finish the triangulation
|
| 2075 |
// copy information into Tri and Walk arrays
|
| 2076 |
Tri = new int[trilen][3]; |
| 2077 |
Walk = new int[trilen][3]; |
| 2078 |
for (int i=0; i<trilen; i++) { |
| 2079 |
Tri[i][0] = tri[i][0]; |
| 2080 |
Tri[i][1] = tri[i][1]; |
| 2081 |
Tri[i][2] = tri[i][2]; |
| 2082 |
Walk[i][0] = walk[i][0]; |
| 2083 |
Walk[i][1] = walk[i][1]; |
| 2084 |
Walk[i][2] = walk[i][2]; |
| 2085 |
} |
| 2086 |
|
| 2087 |
// call more generic method for constructing Vertices and Edges arrays
|
| 2088 |
finish_triang(samples); |
| 2089 |
} |
| 2090 |
|
| 2091 |
static final double[][] m_samples = { // grid 1 x-coordinates |
| 2092 |
{ 65, 142, 215, 315, 373, 435, 39, 118,
|
| 2093 |
202, 320, 373, 455, 40, 114, 206, 307, |
| 2094 |
384, 457, 66, 128, 208, 308, 380, 436, |
| 2095 |
// grid 2 x-coordinates
|
| 2096 |
83, 144, 201, 293, 354, 433, 60, 135, |
| 2097 |
202, 285, 355, 456, 59, 136, 204, 284, |
| 2098 |
362, 456, 75, 138, 207, 285, 363, 438, |
| 2099 |
// grid 3 x-coordinates
|
| 2100 |
90, 153, 217, 292, 358, 441, 61, 145, |
| 2101 |
216, 292, 366, 452, 55, 143, 213, 295, |
| 2102 |
373, 463, 80, 148, 217, 295, 375, 444}, |
| 2103 |
// grid 1 y-coordinates
|
| 2104 |
{ 67, 87, 103, 104, 72, 42, 122, 148,
|
| 2105 |
160, 165, 156, 109, 212, 211, 203, 200, |
| 2106 |
204, 211, 282, 263, 248, 243, 256, 287, |
| 2107 |
// grid 2 y-coordinates
|
| 2108 |
166, 187, 201, 207, 185, 155, 230, 235, |
| 2109 |
240, 241, 235, 210, 283, 275, 270, 267, |
| 2110 |
273, 280, 338, 310, 299, 297, 305, 331, |
| 2111 |
// grid 3 y-coordinates
|
| 2112 |
247, 270, 277, 276, 265, 233, 295, 319, |
| 2113 |
325, 322, 306, 281, 368, 371, 368, 362, |
| 2114 |
372, 376, 464, 431, 417, 418, 424, 455} }; |
| 2115 |
|
| 2116 |
static final int m_lenx = 6; |
| 2117 |
static final int m_leny = 4; |
| 2118 |
static final int m_numgrids = 3; |
| 2119 |
|
| 2120 |
static final Color[] m_col = {Color.black, Color.gray, Color.blue}; |
| 2121 |
|
| 2122 |
static DelaunayOverlap delO = null; |
| 2123 |
|
| 2124 |
/* run 'java visad.DelaunayOverlap' to test the DelaunayOverlap class */
|
| 2125 |
public static void main(String[] argv) throws VisADException { |
| 2126 |
|
| 2127 |
Frame frame = new Frame("DelaunayOverlap"); |
| 2128 |
WindowListener l = new WindowAdapter() { |
| 2129 |
public void windowClosing(WindowEvent e) { |
| 2130 |
System.exit(0); |
| 2131 |
} |
| 2132 |
}; |
| 2133 |
frame.addWindowListener(l); |
| 2134 |
frame.setSize(500, 600); |
| 2135 |
Dimension screenSize = Toolkit.getDefaultToolkit().getScreenSize(); |
| 2136 |
frame.setLocation(screenSize.width/2 - 250, |
| 2137 |
screenSize.height/2 - 300); |
| 2138 |
|
| 2139 |
Panel big_panel = new Panel(); |
| 2140 |
big_panel.setLayout(new GridLayout(1, 1)); |
| 2141 |
frame.add(big_panel); |
| 2142 |
|
| 2143 |
Canvas gcan = new Canvas() { |
| 2144 |
public void paint(Graphics gr) { |
| 2145 |
// draw triangulation
|
| 2146 |
if (delO != null) { |
| 2147 |
gr.setColor(Color.red);
|
| 2148 |
for (int t=0; t<delO.Tri.length; t++) { |
| 2149 |
for (int e=0; e<3; e++) { |
| 2150 |
int gx1 = (int)m_samples[0][delO.Tri[t][e]]; |
| 2151 |
int gy1 = (int)m_samples[1][delO.Tri[t][e]]; |
| 2152 |
int gx2 = (int)m_samples[0][delO.Tri[t][(e+1)%3]]; |
| 2153 |
int gy2 = (int)m_samples[1][delO.Tri[t][(e+1)%3]]; |
| 2154 |
gr.drawLine(gx1, gy1, gx2, gy2); |
| 2155 |
} |
| 2156 |
} |
| 2157 |
} |
| 2158 |
|
| 2159 |
// plot points
|
| 2160 |
int colnum = 0; |
| 2161 |
for (int p=0; p<m_numgrids*m_lenx*m_leny; p++) { |
| 2162 |
int x = (int)m_samples[0][p]; |
| 2163 |
int y = (int)m_samples[1][p]; |
| 2164 |
if (p % (m_lenx*m_leny) == 0) { |
| 2165 |
gr.setColor(m_col[colnum++]); |
| 2166 |
} |
| 2167 |
gr.fillRect(x-2, y-2, 4, 4); |
| 2168 |
} |
| 2169 |
} |
| 2170 |
}; |
| 2171 |
big_panel.add(gcan); |
| 2172 |
|
| 2173 |
System.out.println("Constructing triangulation..."); |
| 2174 |
long startTime = System.currentTimeMillis(); |
| 2175 |
delO = new DelaunayOverlap(m_samples, m_lenx, m_leny);
|
| 2176 |
long endTime = System.currentTimeMillis(); |
| 2177 |
frame.setVisible(true);
|
| 2178 |
System.out.println("Algorithm completed successfully."); |
| 2179 |
double algTime = (double)(endTime-startTime) / 1000; |
| 2180 |
System.out.println("Triangulation took "+algTime+" seconds."); |
| 2181 |
|
| 2182 |
// Test the triangulation for errors
|
| 2183 |
System.out.println("Testing triangulation..."); |
| 2184 |
boolean good = delO.test(m_samples);
|
| 2185 |
if (!good) System.out.println("Errors in triangulation!"); |
| 2186 |
else System.out.println("Triangulation was constructed correctly."); |
| 2187 |
} |
| 2188 |
|
| 2189 |
} |
| 2190 |
|