001: /*
002: * $RCSfile: Numerics.java,v $
003: *
004: * Copyright (c) 2007 Sun Microsystems, Inc. All rights reserved.
005: *
006: * Redistribution and use in source and binary forms, with or without
007: * modification, are permitted provided that the following conditions
008: * are met:
009: *
010: * - Redistribution of source code must retain the above copyright
011: * notice, this list of conditions and the following disclaimer.
012: *
013: * - Redistribution in binary form must reproduce the above copyright
014: * notice, this list of conditions and the following disclaimer in
015: * the documentation and/or other materials provided with the
016: * distribution.
017: *
018: * Neither the name of Sun Microsystems, Inc. or the names of
019: * contributors may be used to endorse or promote products derived
020: * from this software without specific prior written permission.
021: *
022: * This software is provided "AS IS," without a warranty of any
023: * kind. ALL EXPRESS OR IMPLIED CONDITIONS, REPRESENTATIONS AND
024: * WARRANTIES, INCLUDING ANY IMPLIED WARRANTY OF MERCHANTABILITY,
025: * FITNESS FOR A PARTICULAR PURPOSE OR NON-INFRINGEMENT, ARE HEREBY
026: * EXCLUDED. SUN MICROSYSTEMS, INC. ("SUN") AND ITS LICENSORS SHALL
027: * NOT BE LIABLE FOR ANY DAMAGES SUFFERED BY LICENSEE AS A RESULT OF
028: * USING, MODIFYING OR DISTRIBUTING THIS SOFTWARE OR ITS
029: * DERIVATIVES. IN NO EVENT WILL SUN OR ITS LICENSORS BE LIABLE FOR
030: * ANY LOST REVENUE, PROFIT OR DATA, OR FOR DIRECT, INDIRECT, SPECIAL,
031: * CONSEQUENTIAL, INCIDENTAL OR PUNITIVE DAMAGES, HOWEVER CAUSED AND
032: * REGARDLESS OF THE THEORY OF LIABILITY, ARISING OUT OF THE USE OF OR
033: * INABILITY TO USE THIS SOFTWARE, EVEN IF SUN HAS BEEN ADVISED OF THE
034: * POSSIBILITY OF SUCH DAMAGES.
035: *
036: * You acknowledge that this software is not designed, licensed or
037: * intended for use in the design, construction, operation or
038: * maintenance of any nuclear facility.
039: *
040: * $Revision: 1.4 $
041: * $Date: 2007/02/09 17:20:20 $
042: * $State: Exp $
043: */
044:
045: // ----------------------------------------------------------------------
046: //
047: // The reference to Fast Industrial Strength Triangulation (FIST) code
048: // in this release by Sun Microsystems is related to Sun's rewrite of
049: // an early version of FIST. FIST was originally created by Martin
050: // Held and Joseph Mitchell at Stony Brook University and is
051: // incorporated by Sun under an agreement with The Research Foundation
052: // of SUNY (RFSUNY). The current version of FIST is available for
053: // commercial use under a license agreement with RFSUNY on behalf of
054: // the authors and Stony Brook University. Please contact the Office
055: // of Technology Licensing at Stony Brook, phone 631-632-9009, for
056: // licensing information.
057: //
058: // ----------------------------------------------------------------------
059: package com.sun.j3d.utils.geometry;
060:
061: import java.io.*;
062: import java.util.*;
063: import javax.vecmath.*;
064:
065: class Numerics {
066:
067: static double max3(double a, double b, double c) {
068: return (((a) > (b)) ? (((a) > (c)) ? (a) : (c))
069: : (((b) > (c)) ? (b) : (c)));
070: }
071:
072: static double min3(double a, double b, double c) {
073: return (((a) < (b)) ? (((a) < (c)) ? (a) : (c))
074: : (((b) < (c)) ? (b) : (c)));
075: }
076:
077: static boolean lt(double a, double eps) {
078: return ((a) < -eps);
079: }
080:
081: static boolean le(double a, double eps) {
082: return (a <= eps);
083: }
084:
085: static boolean ge(double a, double eps) {
086: return (!((a) <= -eps));
087: }
088:
089: static boolean eq(double a, double eps) {
090: return (((a) <= eps) && !((a) < -eps));
091: }
092:
093: static boolean gt(double a, double eps) {
094: return !((a) <= eps);
095: }
096:
097: static double baseLength(Tuple2f u, Tuple2f v) {
098: double x, y;
099: x = (v).x - (u).x;
100: y = (v).y - (u).y;
101: return Math.abs(x) + Math.abs(y);
102: }
103:
104: static double sideLength(Tuple2f u, Tuple2f v) {
105: double x, y;
106: x = (v).x - (u).x;
107: y = (v).y - (u).y;
108: return x * x + y * y;
109: }
110:
111: /**
112: * This checks whether i3, which is collinear with i1, i2, is
113: * between i1, i2. note that we rely on the lexicographic sorting of the
114: * points!
115: */
116: static boolean inBetween(int i1, int i2, int i3) {
117: return ((i1 <= i3) && (i3 <= i2));
118: }
119:
120: static boolean strictlyInBetween(int i1, int i2, int i3) {
121: return ((i1 < i3) && (i3 < i2));
122: }
123:
124: /**
125: * this method computes the determinant det(points[i],points[j],points[k])
126: * in a consistent way.
127: */
128: static double stableDet2D(Triangulator triRef, int i, int j, int k) {
129: double det;
130: Point2f numericsHP, numericsHQ, numericsHR;
131:
132: // if((triRef.inPointsList(i)==false)||(triRef.inPointsList(j)==false)||
133: // (triRef.inPointsList(k)==false))
134: // System.out.println("Numerics.stableDet2D Not inPointsList " + i + " " + j
135: // + " " + k);
136:
137: if ((i == j) || (i == k) || (j == k)) {
138: det = 0.0;
139: } else {
140: numericsHP = triRef.points[i];
141: numericsHQ = triRef.points[j];
142: numericsHR = triRef.points[k];
143:
144: if (i < j) {
145: if (j < k) /* i < j < k */
146: det = Basic.det2D(numericsHP, numericsHQ,
147: numericsHR);
148: else if (i < k) /* i < k < j */
149: det = -Basic.det2D(numericsHP, numericsHR,
150: numericsHQ);
151: else
152: /* k < i < j */
153: det = Basic.det2D(numericsHR, numericsHP,
154: numericsHQ);
155: } else {
156: if (i < k) /* j < i < k */
157: det = -Basic.det2D(numericsHQ, numericsHP,
158: numericsHR);
159: else if (j < k) /* j < k < i */
160: det = Basic.det2D(numericsHQ, numericsHR,
161: numericsHP);
162: else
163: /* k < j < i */
164: det = -Basic.det2D(numericsHR, numericsHQ,
165: numericsHP);
166: }
167: }
168:
169: return det;
170: }
171:
172: /**
173: * Returns the orientation of the triangle.
174: * @return +1 if the points i, j, k are given in CCW order;
175: * -1 if the points i, j, k are given in CW order;
176: * 0 if the points i, j, k are collinear.
177: */
178: static int orientation(Triangulator triRef, int i, int j, int k) {
179: int ori;
180: double numericsHDet;
181: numericsHDet = stableDet2D(triRef, i, j, k);
182: // System.out.println("orientation : numericsHDet " + numericsHDet);
183: if (lt(numericsHDet, triRef.epsilon))
184: ori = -1;
185: else if (gt(numericsHDet, triRef.epsilon))
186: ori = 1;
187: else
188: ori = 0;
189: return ori;
190: }
191:
192: /**
193: * This method checks whether l is in the cone defined by i, j and j, k
194: */
195: static boolean isInCone(Triangulator triRef, int i, int j, int k,
196: int l, boolean convex) {
197: boolean flag;
198: int numericsHOri1, numericsHOri2;
199:
200: // if((triRef.inPointsList(i)==false)||(triRef.inPointsList(j)==false)||
201: // (triRef.inPointsList(k)==false)||(triRef.inPointsList(l)==false))
202: // System.out.println("Numerics.isInCone Not inPointsList " + i + " " + j
203: // + " " + k + " " + l);
204:
205: flag = true;
206: if (convex) {
207: if (i != j) {
208: numericsHOri1 = orientation(triRef, i, j, l);
209: // System.out.println("isInCone : i != j, numericsHOri1 = " + numericsHOri1);
210: if (numericsHOri1 < 0)
211: flag = false;
212: else if (numericsHOri1 == 0) {
213: if (i < j) {
214: if (!inBetween(i, j, l))
215: flag = false;
216: } else {
217: if (!inBetween(j, i, l))
218: flag = false;
219: }
220: }
221: }
222: if ((j != k) && (flag == true)) {
223: numericsHOri2 = orientation(triRef, j, k, l);
224: // System.out.println("isInCone : ((j != k) && (flag == true)), numericsHOri2 = " +
225: // numericsHOri2);
226: if (numericsHOri2 < 0)
227: flag = false;
228: else if (numericsHOri2 == 0) {
229: if (j < k) {
230: if (!inBetween(j, k, l))
231: flag = false;
232: } else {
233: if (!inBetween(k, j, l))
234: flag = false;
235: }
236: }
237: }
238: } else {
239: numericsHOri1 = orientation(triRef, i, j, l);
240: if (numericsHOri1 <= 0) {
241: numericsHOri2 = orientation(triRef, j, k, l);
242: if (numericsHOri2 < 0)
243: flag = false;
244: }
245: }
246: return flag;
247: }
248:
249: /**
250: * Returns convex angle flag.
251: * @return 0 ... if angle is 180 degrees <br>
252: * 1 ... if angle between 0 and 180 degrees <br>
253: * 2 ... if angle is 0 degrees <br>
254: * -1 ... if angle between 180 and 360 degrees <br>
255: * -2 ... if angle is 360 degrees <br>
256: */
257: static int isConvexAngle(Triangulator triRef, int i, int j, int k,
258: int ind) {
259: int angle;
260: double numericsHDot;
261: int numericsHOri1;
262: Point2f numericsHP, numericsHQ;
263:
264: // if((triRef.inPointsList(i)==false)||(triRef.inPointsList(j)==false)||
265: // (triRef.inPointsList(k)==false))
266: // System.out.println("Numerics.isConvexAngle: Not inPointsList " + i + " " + j
267: // + " " + k);
268:
269: if (i == j) {
270: if (j == k) {
271: // all three vertices are identical; we set the angle to 1 in
272: // order to enable clipping of j.
273: return 1;
274: } else {
275: // two of the three vertices are identical; we set the angle to 1
276: // in order to enable clipping of j.
277: return 1;
278: }
279: } else if (j == k) {
280: // two vertices are identical. we could either determine the angle
281: // by means of yet another lengthy analysis, or simply set the
282: // angle to -1. using -1 means to err on the safe side, as all the
283: // incarnations of this vertex will be clipped right at the start
284: // of the ear-clipping algorithm. thus, eventually there will be no
285: // other duplicates at this vertex position, and the regular
286: // classification of angles will yield the correct answer for j.
287: return -1;
288: } else {
289: numericsHOri1 = orientation(triRef, i, j, k);
290: // System.out.println("i " + i + " j " + j + " k " + k + " ind " + ind +
291: // ". In IsConvexAngle numericsHOri1 is " +
292: // numericsHOri1);
293: if (numericsHOri1 > 0) {
294: angle = 1;
295: } else if (numericsHOri1 < 0) {
296: angle = -1;
297: } else {
298: // 0, 180, or 360 degrees.
299: numericsHP = new Point2f();
300: numericsHQ = new Point2f();
301: Basic.vectorSub2D(triRef.points[i], triRef.points[j],
302: numericsHP);
303: Basic.vectorSub2D(triRef.points[k], triRef.points[j],
304: numericsHQ);
305: numericsHDot = Basic.dotProduct2D(numericsHP,
306: numericsHQ);
307: if (numericsHDot < 0.0) {
308: // 180 degrees.
309: angle = 0;
310: } else {
311: // 0 or 360 degrees? this cannot be judged locally, and more
312: // work is needed.
313:
314: angle = spikeAngle(triRef, i, j, k, ind);
315: // System.out.println("SpikeAngle return is "+ angle);
316: }
317: }
318: }
319: return angle;
320: }
321:
322: /**
323: * This method checks whether point i4 is inside of or on the boundary
324: * of the triangle i1, i2, i3.
325: */
326: static boolean pntInTriangle(Triangulator triRef, int i1, int i2,
327: int i3, int i4) {
328: boolean inside;
329: int numericsHOri1;
330:
331: inside = false;
332: numericsHOri1 = orientation(triRef, i2, i3, i4);
333: if (numericsHOri1 >= 0) {
334: numericsHOri1 = orientation(triRef, i1, i2, i4);
335: if (numericsHOri1 >= 0) {
336: numericsHOri1 = orientation(triRef, i3, i1, i4);
337: if (numericsHOri1 >= 0)
338: inside = true;
339: }
340: }
341: return inside;
342: }
343:
344: /**
345: * This method checks whether point i4 is inside of or on the boundary
346: * of the triangle i1, i2, i3. it also returns a classification if i4 is
347: * on the boundary of the triangle (except for the edge i2, i3).
348: */
349: static boolean vtxInTriangle(Triangulator triRef, int i1, int i2,
350: int i3, int i4, int[] type) {
351: boolean inside;
352: int numericsHOri1;
353:
354: inside = false;
355: numericsHOri1 = orientation(triRef, i2, i3, i4);
356: if (numericsHOri1 >= 0) {
357: numericsHOri1 = orientation(triRef, i1, i2, i4);
358: if (numericsHOri1 > 0) {
359: numericsHOri1 = orientation(triRef, i3, i1, i4);
360: if (numericsHOri1 > 0) {
361: inside = true;
362: type[0] = 0;
363: } else if (numericsHOri1 == 0) {
364: inside = true;
365: type[0] = 1;
366: }
367: } else if (numericsHOri1 == 0) {
368: numericsHOri1 = orientation(triRef, i3, i1, i4);
369: if (numericsHOri1 > 0) {
370: inside = true;
371: type[0] = 2;
372: } else if (numericsHOri1 == 0) {
373: inside = true;
374: type[0] = 3;
375: }
376: }
377: }
378: return inside;
379: }
380:
381: /**
382: * Checks whether the line segments i1, i2 and i3, i4 intersect. no
383: * intersection is reported if they intersect at a common vertex.
384: * the function assumes that i1 <= i2 and i3 <= i4. if i3 or i4 lies
385: * on i1, i2 then an intersection is reported, but no intersection is
386: * reported if i1 or i2 lies on i3, i4. this function is not symmetric!
387: */
388: static boolean segIntersect(Triangulator triRef, int i1, int i2,
389: int i3, int i4, int i5) {
390: int ori1, ori2, ori3, ori4;
391:
392: // if((triRef.inPointsList(i1)==false)||(triRef.inPointsList(i2)==false)||
393: // (triRef.inPointsList(i3)==false)||(triRef.inPointsList(i4)==false))
394: // System.out.println("Numerics.segIntersect Not inPointsList " + i1 + " " + i2
395: // + " " + i3 + " " + i4);
396: //
397: // if((i1 > i2) || (i3 > i4))
398: // System.out.println("Numerics.segIntersect i1>i2 or i3>i4 " + i1 + " " + i2
399: // + " " + i3 + " " + i4);
400:
401: if ((i1 == i2) || (i3 == i4))
402: return false;
403: if ((i1 == i3) && (i2 == i4))
404: return true;
405:
406: if ((i3 == i5) || (i4 == i5))
407: ++(triRef.identCntr);
408:
409: ori3 = orientation(triRef, i1, i2, i3);
410: ori4 = orientation(triRef, i1, i2, i4);
411: if (((ori3 == 1) && (ori4 == 1))
412: || ((ori3 == -1) && (ori4 == -1)))
413: return false;
414:
415: if (ori3 == 0) {
416: if (strictlyInBetween(i1, i2, i3))
417: return true;
418: if (ori4 == 0) {
419: if (strictlyInBetween(i1, i2, i4))
420: return true;
421: } else
422: return false;
423: } else if (ori4 == 0) {
424: if (strictlyInBetween(i1, i2, i4))
425: return true;
426: else
427: return false;
428: }
429:
430: ori1 = orientation(triRef, i3, i4, i1);
431: ori2 = orientation(triRef, i3, i4, i2);
432: if (((ori1 <= 0) && (ori2 <= 0))
433: || ((ori1 >= 0) && (ori2 >= 0)))
434: return false;
435:
436: return true;
437: }
438:
439: /**
440: * this function computes a quality measure of a triangle i, j, k.
441: * it returns the ratio `base / height', where base is the length of the
442: * longest side of the triangle, and height is the normal distance
443: * between the vertex opposite of the base side and the base side. (as
444: * usual, we again use the l1-norm for distances.)
445: */
446: static double getRatio(Triangulator triRef, int i, int j, int k) {
447: double area, a, b, c, base, ratio;
448: Point2f p, q, r;
449:
450: // if((triRef.inPointsList(i)==false)||(triRef.inPointsList(j)==false)||
451: // (triRef.inPointsList(k)==false))
452: // System.out.println("Numerics.getRatio: Not inPointsList " + i + " " + j
453: // + " " + k);
454:
455: p = triRef.points[i];
456: q = triRef.points[j];
457: r = triRef.points[k];
458:
459: a = baseLength(p, q);
460: b = baseLength(p, r);
461: c = baseLength(r, q);
462: base = max3(a, b, c);
463:
464: if ((10.0 * a) < Math.min(b, c))
465: return 0.1;
466:
467: area = stableDet2D(triRef, i, j, k);
468: if (lt(area, triRef.epsilon)) {
469: area = -area;
470: } else if (!gt(area, triRef.epsilon)) {
471: if (base > a)
472: return 0.1;
473: else
474: return Double.MAX_VALUE;
475: }
476:
477: ratio = base * base / area;
478:
479: if (ratio < 10.0)
480: return ratio;
481: else {
482: if (a < base)
483: return 0.1;
484: else
485: return ratio;
486: }
487: }
488:
489: static int spikeAngle(Triangulator triRef, int i, int j, int k,
490: int ind) {
491: int ind1, ind2, ind3;
492: int i1, i2, i3;
493:
494: // if((triRef.inPointsList(i)==false)||(triRef.inPointsList(j)==false)||
495: // (triRef.inPointsList(k)==false))
496: // System.out.println("Numerics.spikeAngle: Not inPointsList " + i + " " + j
497: // + " " + k);
498:
499: ind2 = ind;
500: i2 = triRef.fetchData(ind2);
501:
502: // if(i2 != j)
503: // System.out.println("Numerics.spikeAngle: i2 != j " + i2 + " " + j );
504:
505: ind1 = triRef.fetchPrevData(ind2);
506: i1 = triRef.fetchData(ind1);
507:
508: // if(i1 != i)
509: // System.out.println("Numerics.spikeAngle: i1 != i " + i1 + " " + i );
510:
511: ind3 = triRef.fetchNextData(ind2);
512: i3 = triRef.fetchData(ind3);
513:
514: // if(i3 != k)
515: // System.out.println("Numerics.spikeAngle: i3 != k " + i3 + " " + k );
516:
517: return recSpikeAngle(triRef, i, j, k, ind1, ind3);
518: }
519:
520: static int recSpikeAngle(Triangulator triRef, int i1, int i2,
521: int i3, int ind1, int ind3) {
522: int ori, ori1, ori2, i0, ii1, ii2;
523: Point2f pq, pr;
524: double dot;
525:
526: if (ind1 == ind3) {
527: // all points are collinear??? well, then it does not really matter
528: // which angle is returned. perhaps, -2 is the best bet as my code
529: // likely regards this contour as a hole.
530: return -2;
531: }
532:
533: if (i1 != i3) {
534: if (i1 < i2) {
535: ii1 = i1;
536: ii2 = i2;
537: } else {
538: ii1 = i2;
539: ii2 = i1;
540: }
541: if (inBetween(ii1, ii2, i3)) {
542: i2 = i3;
543: ind3 = triRef.fetchNextData(ind3);
544: i3 = triRef.fetchData(ind3);
545:
546: if (ind1 == ind3)
547: return 2;
548: ori = orientation(triRef, i1, i2, i3);
549: if (ori > 0)
550: return 2;
551: else if (ori < 0)
552: return -2;
553: else
554: return recSpikeAngle(triRef, i1, i2, i3, ind1, ind3);
555: } else {
556: i2 = i1;
557: ind1 = triRef.fetchPrevData(ind1);
558: i1 = triRef.fetchData(ind1);
559: if (ind1 == ind3)
560: return 2;
561: ori = orientation(triRef, i1, i2, i3);
562: if (ori > 0)
563: return 2;
564: else if (ori < 0)
565: return -2;
566: else
567: return recSpikeAngle(triRef, i1, i2, i3, ind1, ind3);
568: }
569: } else {
570: i0 = i2;
571: i2 = i1;
572: ind1 = triRef.fetchPrevData(ind1);
573: i1 = triRef.fetchData(ind1);
574:
575: if (ind1 == ind3)
576: return 2;
577: ind3 = triRef.fetchNextData(ind3);
578: i3 = triRef.fetchData(ind3);
579: if (ind1 == ind3)
580: return 2;
581: ori = orientation(triRef, i1, i2, i3);
582: if (ori > 0) {
583: ori1 = orientation(triRef, i1, i2, i0);
584: if (ori1 > 0) {
585: ori2 = orientation(triRef, i2, i3, i0);
586: if (ori2 > 0)
587: return -2;
588: }
589: return 2;
590: } else if (ori < 0) {
591: ori1 = orientation(triRef, i2, i1, i0);
592: if (ori1 > 0) {
593: ori2 = orientation(triRef, i3, i2, i0);
594: if (ori2 > 0)
595: return 2;
596: }
597: return -2;
598: } else {
599: pq = new Point2f();
600: Basic.vectorSub2D(triRef.points[i1], triRef.points[i2],
601: pq);
602: pr = new Point2f();
603: Basic.vectorSub2D(triRef.points[i3], triRef.points[i2],
604: pr);
605: dot = Basic.dotProduct2D(pq, pr);
606: if (dot < 0.0) {
607: ori = orientation(triRef, i2, i1, i0);
608: if (ori > 0)
609: return 2;
610: else
611: return -2;
612: } else {
613: return recSpikeAngle(triRef, i1, i2, i3, ind1, ind3);
614: }
615: }
616: }
617: }
618:
619: /**
620: * computes the signed angle between p, p1 and p, p2.
621: *
622: * warning: this function does not handle a 180-degree angle correctly!
623: * (this is no issue in our application, as we will always compute
624: * the angle centered at the mid-point of a valid diagonal.)
625: */
626: static double angle(Triangulator triRef, Point2f p, Point2f p1,
627: Point2f p2) {
628: int sign;
629: double angle1, angle2, angle;
630: Point2f v1, v2;
631:
632: sign = Basic.signEps(Basic.det2D(p2, p, p1), triRef.epsilon);
633:
634: if (sign == 0)
635: return 0.0;
636:
637: v1 = new Point2f();
638: v2 = new Point2f();
639: Basic.vectorSub2D(p1, p, v1);
640: Basic.vectorSub2D(p2, p, v2);
641:
642: angle1 = Math.atan2(v1.y, v1.x);
643: angle2 = Math.atan2(v2.y, v2.x);
644:
645: if (angle1 < 0.0)
646: angle1 += 2.0 * Math.PI;
647: if (angle2 < 0.0)
648: angle2 += 2.0 * Math.PI;
649:
650: angle = angle1 - angle2;
651: if (angle > Math.PI)
652: angle = 2.0 * Math.PI - angle;
653: else if (angle < -Math.PI)
654: angle = 2.0 * Math.PI + angle;
655:
656: if (sign == 1) {
657: if (angle < 0.0)
658: return -angle;
659: else
660: return angle;
661: } else {
662: if (angle > 0.0)
663: return -angle;
664: else
665: return angle;
666: }
667: }
668:
669: }
|