svn-gvsig-desktop / trunk / extensions / extGraph / src / org / gvsig / fmap / algorithm / triangulation / visad / DelaunayWatson.java @ 39203
History | View | Annotate | Download (11.1 KB)
| 1 |
//
|
|---|---|
| 2 |
// DelaunayWatson.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 |
import java.util.Random; |
| 30 |
|
| 31 |
/*
|
| 32 |
* The Delaunay triangulation/tetrahedralization algorithm in this class is
|
| 33 |
* originally from nnsort.c by David F. Watson:
|
| 34 |
*
|
| 35 |
* nnsort() finds the Delaunay triangulation of the two- or three-component
|
| 36 |
* vectors in 'data_list' and returns a list of simplex vertices in 'vertices'
|
| 37 |
* with the corresponding circumcentre and squared radius in the rows of
|
| 38 |
* 'circentres'. nnsort() also can be used to find the ordered convex hull of
|
| 39 |
* the two- or three-component vectors in 'data_list' and returns a list of
|
| 40 |
* (d-1)-facet vertices in 'vertices' (dummy filename for 'circentres' must be
|
| 41 |
* used). nnsort() was written by Dave Watson and uses the algorithm described
|
| 42 |
* in - Watson, D.F., 1981, Computing the n-dimensional Delaunay tessellation
|
| 43 |
* with application to Voronoi polytopes: The Computer J., 24(2), p. 167-172.
|
| 44 |
*
|
| 45 |
* additional information about this algorithm can be found in - CONTOURING: A
|
| 46 |
* guide to the analysis and display of spatial data, by David F. Watson,
|
| 47 |
* Pergamon Press, 1992, ISBN 0 08 040286 0
|
| 48 |
*/
|
| 49 |
|
| 50 |
/**
|
| 51 |
* DelaunayWatson represents an O(N^2) method with low overhead to find the
|
| 52 |
* Delaunay triangulation or tetrahedralization of a set of samples of R^2 or
|
| 53 |
* R^3.
|
| 54 |
* <P>
|
| 55 |
*/
|
| 56 |
public class DelaunayWatson extends Delaunay { |
| 57 |
|
| 58 |
private static final float BIGNUM = (float) 1E37; |
| 59 |
private static final float EPSILON = 0.00001f; |
| 60 |
// temporary storage size factor
|
| 61 |
private static final int TSIZE = 75; |
| 62 |
// factor (>=1) for radius of control points
|
| 63 |
private static final float RANGE = 10.0f; |
| 64 |
|
| 65 |
public DelaunayWatson(double[][] samples) throws VisADException { |
| 66 |
int dim = samples.length;
|
| 67 |
int nrs = samples[0].length; |
| 68 |
|
| 69 |
double xx, yy, bgs;
|
| 70 |
int i0, i1, i2, i3, i4, i5, i6, i7, i8, i9, i11;
|
| 71 |
int[] ii = new int[3]; |
| 72 |
int dm, dim1, nts, tsz;
|
| 73 |
|
| 74 |
double[][] mxy = new double[2][dim]; |
| 75 |
for (i0 = 0; i0 < dim; i0++) |
| 76 |
mxy[0][i0] = -(mxy[1][i0] = BIGNUM); |
| 77 |
dim1 = dim + 1;
|
| 78 |
double[][] wrk = new double[dim][dim1]; |
| 79 |
for (i0 = 0; i0 < dim; i0++) |
| 80 |
for (i1 = 0; i1 < dim1; i1++) |
| 81 |
wrk[i0][i1] = -RANGE; |
| 82 |
for (i0 = 0; i0 < dim; i0++) |
| 83 |
wrk[i0][i0] = RANGE * (3 * dim - 1); |
| 84 |
|
| 85 |
double[][] pts = new double[nrs + dim1][dim]; |
| 86 |
for (i0 = 0; i0 < nrs; i0++) { |
| 87 |
if (dim < 3) { |
| 88 |
pts[i0][0] = samples[0][i0]; |
| 89 |
pts[i0][1] = samples[1][i0]; |
| 90 |
} else {
|
| 91 |
pts[i0][0] = samples[0][i0]; |
| 92 |
pts[i0][1] = samples[1][i0]; |
| 93 |
pts[i0][2] = samples[2][i0]; |
| 94 |
} |
| 95 |
// compute bounding box
|
| 96 |
for (i1 = 0; i1 < dim; i1++) { |
| 97 |
if (mxy[0][i1] < pts[i0][i1]) |
| 98 |
mxy[0][i1] = pts[i0][i1]; // max |
| 99 |
if (mxy[1][i1] > pts[i0][i1]) |
| 100 |
mxy[1][i1] = pts[i0][i1]; // min |
| 101 |
} |
| 102 |
} |
| 103 |
|
| 104 |
for (bgs = 0, i0 = 0; i0 < dim; i0++) { |
| 105 |
mxy[0][i0] -= mxy[1][i0]; |
| 106 |
if (bgs < mxy[0][i0]) |
| 107 |
bgs = mxy[0][i0];
|
| 108 |
} |
| 109 |
// now bgs = largest range
|
| 110 |
// add random perturbations to points
|
| 111 |
bgs *= EPSILON; |
| 112 |
|
| 113 |
Random rand = new Random(367); |
| 114 |
for (i0 = 0; i0 < nrs; i0++) |
| 115 |
for (i1 = 0; i1 < dim; i1++) { |
| 116 |
// random numbers [0, 1]
|
| 117 |
pts[i0][i1] += bgs * (0.5 - rand.nextDouble() / 0x7fffffff); |
| 118 |
} |
| 119 |
for (i0 = 0; i0 < dim1; i0++) |
| 120 |
for (i1 = 0; i1 < dim; i1++) { |
| 121 |
pts[nrs + i0][i1] = mxy[1][i1] + wrk[i1][i0] * mxy[0][i1]; |
| 122 |
} |
| 123 |
for (i1 = 1, i0 = 2; i0 < dim1; i0++) |
| 124 |
i1 *= i0; |
| 125 |
tsz = TSIZE * i1; |
| 126 |
int[][] tmp = new int[tsz + 1][dim]; |
| 127 |
// storage allocation - increase value of `i1' for 3D if necessary
|
| 128 |
i1 *= (nrs + 50 * i1);
|
| 129 |
/* WLH 4 Nov 97 */
|
| 130 |
if (dim == 3) |
| 131 |
i1 *= 10;
|
| 132 |
/* end WLH 4 Nov 97 */
|
| 133 |
int[] id = new int[i1]; |
| 134 |
for (i0 = 0; i0 < i1; i0++) |
| 135 |
id[i0] = i0; |
| 136 |
int[][] a3s = new int[i1][dim1]; |
| 137 |
double[][] ccr = new double[i1][dim1]; |
| 138 |
for (a3s[0][0] = nrs, i0 = 1; i0 < dim1; i0++) |
| 139 |
a3s[0][i0] = a3s[0][i0 - 1] + 1; |
| 140 |
for (ccr[0][dim] = BIGNUM, i0 = 0; i0 < dim; i0++) |
| 141 |
ccr[0][i0] = 0; |
| 142 |
nts = i4 = 1;
|
| 143 |
dm = dim - 1;
|
| 144 |
for (i0 = 0; i0 < nrs; i0++) { |
| 145 |
i1 = i7 = -1;
|
| 146 |
i9 = 0;
|
| 147 |
Loop3: for (i11 = 0; i11 < nts; i11++) { |
| 148 |
i1++; |
| 149 |
while (a3s[i1][0] < 0) |
| 150 |
i1++; |
| 151 |
xx = ccr[i1][dim]; |
| 152 |
for (i2 = 0; i2 < dim; i2++) { |
| 153 |
xx -= (pts[i0][i2] - ccr[i1][i2]) |
| 154 |
* (pts[i0][i2] - ccr[i1][i2]); |
| 155 |
if (xx < 0) |
| 156 |
continue Loop3;
|
| 157 |
} |
| 158 |
i9--; |
| 159 |
i4--; |
| 160 |
id[i4] = i1; |
| 161 |
Loop2: for (i2 = 0; i2 < dim1; i2++) { |
| 162 |
ii[0] = 0; |
| 163 |
if (ii[0] == i2) |
| 164 |
ii[0]++;
|
| 165 |
for (i3 = 1; i3 < dim; i3++) { |
| 166 |
ii[i3] = ii[i3 - 1] + 1; |
| 167 |
if (ii[i3] == i2)
|
| 168 |
ii[i3]++; |
| 169 |
} |
| 170 |
if (i7 > dm) {
|
| 171 |
i8 = i7; |
| 172 |
Loop1: for (i3 = 0; i3 <= i8; i3++) { |
| 173 |
for (i5 = 0; i5 < dim; i5++) { |
| 174 |
if (a3s[i1][ii[i5]] != tmp[i3][i5])
|
| 175 |
continue Loop1;
|
| 176 |
} |
| 177 |
for (i6 = 0; i6 < dim; i6++) |
| 178 |
tmp[i3][i6] = tmp[i8][i6]; |
| 179 |
i7--; |
| 180 |
continue Loop2;
|
| 181 |
} |
| 182 |
} |
| 183 |
if (++i7 > tsz) {
|
| 184 |
int newtsz = 2 * tsz; |
| 185 |
int[][] newtmp = new int[newtsz + 1][dim]; |
| 186 |
System.arraycopy(tmp, 0, newtmp, 0, tsz); |
| 187 |
tsz = newtsz; |
| 188 |
tmp = newtmp; |
| 189 |
// WLH 23 july 97
|
| 190 |
// throw new VisADException(
|
| 191 |
// "DelaunayWatson: Temporary storage exceeded");
|
| 192 |
} |
| 193 |
for (i3 = 0; i3 < dim; i3++) |
| 194 |
tmp[i7][i3] = a3s[i1][ii[i3]]; |
| 195 |
} |
| 196 |
a3s[i1][0] = -1; |
| 197 |
} |
| 198 |
for (i1 = 0; i1 <= i7; i1++) { |
| 199 |
for (i2 = 0; i2 < dim; i2++) { |
| 200 |
for (wrk[i2][dim] = 0, i3 = 0; i3 < dim; i3++) { |
| 201 |
wrk[i2][i3] = pts[tmp[i1][i2]][i3] - pts[i0][i3]; |
| 202 |
wrk[i2][dim] += wrk[i2][i3] |
| 203 |
* (pts[tmp[i1][i2]][i3] + pts[i0][i3]) / 2;
|
| 204 |
} |
| 205 |
} |
| 206 |
if (dim < 3) { |
| 207 |
xx = wrk[0][0] * wrk[1][1] - wrk[1][0] * wrk[0][1]; |
| 208 |
ccr[id[i4]][0] = (wrk[0][2] * wrk[1][1] - wrk[1][2] |
| 209 |
* wrk[0][1]) |
| 210 |
/ xx; |
| 211 |
ccr[id[i4]][1] = (wrk[0][0] * wrk[1][2] - wrk[1][0] |
| 212 |
* wrk[0][2]) |
| 213 |
/ xx; |
| 214 |
} else {
|
| 215 |
xx = (wrk[0][0] * (wrk[1][1] * wrk[2][2] - wrk[2][1] |
| 216 |
* wrk[1][2])) |
| 217 |
- (wrk[0][1] * (wrk[1][0] * wrk[2][2] - wrk[2][0] |
| 218 |
* wrk[1][2])) |
| 219 |
+ (wrk[0][2] * (wrk[1][0] * wrk[2][1] - wrk[2][0] |
| 220 |
* wrk[1][1])); |
| 221 |
ccr[id[i4]][0] = ((wrk[0][3] * (wrk[1][1] * wrk[2][2] - wrk[2][1] |
| 222 |
* wrk[1][2])) |
| 223 |
- (wrk[0][1] * (wrk[1][3] * wrk[2][2] - wrk[2][3] |
| 224 |
* wrk[1][2])) + (wrk[0][2] * (wrk[1][3] |
| 225 |
* wrk[2][1] - wrk[2][3] * wrk[1][1]))) |
| 226 |
/ xx; |
| 227 |
ccr[id[i4]][1] = ((wrk[0][0] * (wrk[1][3] * wrk[2][2] - wrk[2][3] |
| 228 |
* wrk[1][2])) |
| 229 |
- (wrk[0][3] * (wrk[1][0] * wrk[2][2] - wrk[2][0] |
| 230 |
* wrk[1][2])) + (wrk[0][2] * (wrk[1][0] |
| 231 |
* wrk[2][3] - wrk[2][0] * wrk[1][3]))) |
| 232 |
/ xx; |
| 233 |
ccr[id[i4]][2] = ((wrk[0][0] * (wrk[1][1] * wrk[2][3] - wrk[2][1] |
| 234 |
* wrk[1][3])) |
| 235 |
- (wrk[0][1] * (wrk[1][0] * wrk[2][3] - wrk[2][0] |
| 236 |
* wrk[1][3])) + (wrk[0][3] * (wrk[1][0] |
| 237 |
* wrk[2][1] - wrk[2][0] * wrk[1][1]))) |
| 238 |
/ xx; |
| 239 |
} |
| 240 |
for (ccr[id[i4]][dim] = 0, i2 = 0; i2 < dim; i2++) { |
| 241 |
ccr[id[i4]][dim] += (pts[i0][i2] - ccr[id[i4]][i2]) |
| 242 |
* (pts[i0][i2] - ccr[id[i4]][i2]); |
| 243 |
a3s[id[i4]][i2] = tmp[i1][i2]; |
| 244 |
} |
| 245 |
a3s[id[i4]][dim] = i0; |
| 246 |
i4++; |
| 247 |
i9++; |
| 248 |
} |
| 249 |
nts += i9; |
| 250 |
} |
| 251 |
|
| 252 |
/*
|
| 253 |
* OUTPUT is in a3s ARRAY needed output is: Tri - array of pointers from
|
| 254 |
* triangles or tetrahedra to their corresponding vertices Vertices -
|
| 255 |
* array of pointers from vertices to their corresponding triangles or
|
| 256 |
* tetrahedra Walk - array of pointers from triangles or tetrahedra to
|
| 257 |
* neighboring triangles or tetrahedra Edges - array of pointers from
|
| 258 |
* each triangle or tetrahedron's edges to their corresponding triangles
|
| 259 |
* or tetrahedra
|
| 260 |
*
|
| 261 |
* helpers: nverts - number of triangles or tetrahedra per vertex
|
| 262 |
*/
|
| 263 |
|
| 264 |
// compute number of triangles or tetrahedra
|
| 265 |
int[] nverts = new int[nrs]; |
| 266 |
for (int i = 0; i < nrs; i++) |
| 267 |
nverts[i] = 0;
|
| 268 |
int ntris = 0; |
| 269 |
i0 = -1;
|
| 270 |
for (i11 = 0; i11 < nts; i11++) { |
| 271 |
i0++; |
| 272 |
while (a3s[i0][0] < 0) |
| 273 |
i0++; |
| 274 |
if (a3s[i0][0] < nrs) { |
| 275 |
ntris++; |
| 276 |
if (dim < 3) { |
| 277 |
nverts[a3s[i0][0]]++;
|
| 278 |
nverts[a3s[i0][1]]++;
|
| 279 |
nverts[a3s[i0][2]]++;
|
| 280 |
} else {
|
| 281 |
nverts[a3s[i0][0]]++;
|
| 282 |
nverts[a3s[i0][1]]++;
|
| 283 |
nverts[a3s[i0][2]]++;
|
| 284 |
nverts[a3s[i0][3]]++;
|
| 285 |
} |
| 286 |
} |
| 287 |
} |
| 288 |
Vertices = new int[nrs][]; |
| 289 |
for (int i = 0; i < nrs; i++) |
| 290 |
Vertices[i] = new int[nverts[i]]; |
| 291 |
for (int i = 0; i < nrs; i++) |
| 292 |
nverts[i] = 0;
|
| 293 |
|
| 294 |
// build Tri & Vertices components
|
| 295 |
Tri = new int[ntris][dim1]; |
| 296 |
int a, b, c, d;
|
| 297 |
int itri = 0; |
| 298 |
i0 = -1;
|
| 299 |
for (i11 = 0; i11 < nts; i11++) { |
| 300 |
i0++; |
| 301 |
while (a3s[i0][0] < 0) |
| 302 |
i0++; |
| 303 |
if (a3s[i0][0] < nrs) { |
| 304 |
if (dim < 3) { |
| 305 |
a = a3s[i0][0];
|
| 306 |
b = a3s[i0][1];
|
| 307 |
c = a3s[i0][2];
|
| 308 |
Vertices[a][nverts[a]] = itri; |
| 309 |
nverts[a]++; |
| 310 |
Vertices[b][nverts[b]] = itri; |
| 311 |
nverts[b]++; |
| 312 |
Vertices[c][nverts[c]] = itri; |
| 313 |
nverts[c]++; |
| 314 |
Tri[itri][0] = a;
|
| 315 |
Tri[itri][1] = b;
|
| 316 |
Tri[itri][2] = c;
|
| 317 |
} else {
|
| 318 |
a = a3s[i0][0];
|
| 319 |
b = a3s[i0][1];
|
| 320 |
c = a3s[i0][2];
|
| 321 |
d = a3s[i0][3];
|
| 322 |
Vertices[a][nverts[a]] = itri; |
| 323 |
nverts[a]++; |
| 324 |
Vertices[b][nverts[b]] = itri; |
| 325 |
nverts[b]++; |
| 326 |
Vertices[c][nverts[c]] = itri; |
| 327 |
nverts[c]++; |
| 328 |
Vertices[d][nverts[d]] = itri; |
| 329 |
nverts[d]++; |
| 330 |
Tri[itri][0] = a;
|
| 331 |
Tri[itri][1] = b;
|
| 332 |
Tri[itri][2] = c;
|
| 333 |
Tri[itri][3] = d;
|
| 334 |
} |
| 335 |
itri++; |
| 336 |
} |
| 337 |
} |
| 338 |
|
| 339 |
// call more generic method for constructing Walk and Edges arrays
|
| 340 |
finish_triang(samples); |
| 341 |
} |
| 342 |
|
| 343 |
/*
|
| 344 |
* DO NOT DELETE THIS COMMENTED CODE IT CONTAINS ALGORITHM DETAILS NOT CAST
|
| 345 |
* INTO JAVA (YET?) i0 = -1; for (i11=0; i11<nts; i11++) { i0++; while
|
| 346 |
* (a3s[i0][0] < 0) i0++; if (a3s[i0][0] < nrs) { for (i1=0; i1<dim; i1++)
|
| 347 |
* for (i2=0; i2<dim; i2++) { wrk[i1][i2] = pts[a3s[i0][i1]][i2] -
|
| 348 |
* pts[a3s[i0][dim]][i2]; } if (dim < 3) { xx = wrk[0][0] * wrk[1][1] -
|
| 349 |
* wrk[0][1] * wrk[1][0]; if (fabs(xx) > EPSILON) { if (xx < 0)
|
| 350 |
* fprintf(afile,"%d %d %d\n",a3s[i0][0],a3s[i0][2],a3s[i0][1]); else
|
| 351 |
* fprintf(afile,"%d %d %d\n",a3s[i0][0],a3s[i0][1],a3s[i0][2]);
|
| 352 |
* fprintf(bfile,"%e %e %e\n",ccr[i0][0],ccr[i0][1],ccr[i0][2]); } } else {
|
| 353 |
* xx = ((wrk[0][0] * (wrk[1][1] * wrk[2][2] - wrk[2][1] * wrk[1][2])) -
|
| 354 |
* (wrk[0][1] * (wrk[1][0] * wrk[2][2] - wrk[2][0] * wrk[1][2])) +
|
| 355 |
* (wrk[0][2] * (wrk[1][0] * wrk[2][1] - wrk[2][0] * wrk[1][1]))); if
|
| 356 |
* (fabs(xx) > EPSILON) { if (xx < 0) fprintf(afile,"%d %d %d %d\n",
|
| 357 |
* a3s[i0][0],a3s[i0][1],a3s[i0][3],a3s[i0][2]); else fprintf(afile,"%d %d
|
| 358 |
* %d %d\n", a3s[i0][0],a3s[i0][1],a3s[i0][2],a3s[i0][3]); fprintf(bfile,"%e
|
| 359 |
* %e %e %e\n", ccr[i0][0],ccr[i0][1],ccr[i0][2],ccr[i0][3]); } } } }
|
| 360 |
*/
|
| 361 |
|
| 362 |
} |