构造Delaunay三角形网格

Delaunay是一种在离散点序列中快速构造三角形网格的方法，本代码依据的Delaunay三角形的性质：在已知的Dalaunay三角化的网格上加入一点P，只需要删除所有外接圆包含此点的三角形，并连接P与所有可见的点(即连接后不会与其他边相交)，则形成的网格仍然满足Delaunay三角剖分的条件。

方法：
1、构造超大三角形，使得所有离散点均落在该三角形的内部；
2、以该超大三角形作为Delaunay三角形集D的首个成员；
3、对所有离散点集里的每个点，搜索D中满足外接圆包含该点的三角形集T；
4、新点与T构成三角形集N，在D中删除T，并加入N；
5、重复第3、4步；
6、删除D中所有与超大三角形有关的三角形。

 

---

Delaunay.java

---

001  /**
002     Delaunay 方法构造三角形网格
003 
004         PACKAGE: cma.common.isoline
005        FILENAME: Delaunay.java
006        LANGUAGE: Java2 v1.4
007        ORIGINAL: none
008     DESCRIPTION:
009          CREATE: 2007-07-06 17:43:49
010           UPDATE 2007-07-08
011        MODIFIER: 刘泽军 (BJ0773@gmail.com)
012                  广西气象减灾研究所
013                  Guangxi Institude of Meteorology and Disaster-reducing Research(GIMDR)
014 
015       REFERENCE: 参考了  <b> Delaunay.pas </b>
016 
017         COMPILE: javac TriangleVertex.java Delaunay.java
018 
019     由于研究不透其Triangulate方法，故按以下条件重写：
020     1、搜索满足外接圆包含新点的三角形序列；
021     2、在此序列中搜索非共有的边，即此三角形序列构成的边界；
022     3、判断新点与某三角形中作为边界的边是否可以建立新三角形，条件是：新点与该边对应的顶点同处该边界边的同侧。
023 
024  */
025 
026  /*
027  用法：
028       int     points  = 10000;
029       Point2D.Double[]    pts = new Point2D.Double[points];
030       for(int i=0;i<points;i++) {
031           pts[i]  = new Point2D.Double(//在中国区域(70-140E, 10-60N)产生随机点
032                           70.0 + 69.0 * Math.random() + Math.random(),
033                           10.0 + 49.0 * Math.random() + Math.random());
034       }
035 
036       Delaunay    delaunayT   = new Delaunay(pts);//创建Delaunay对象
037       int         nTriangle   = delaunayT.build();//创建Delaunay网络
038 
039       int             width   = 1600;
040       int             height  = 1024;
041       BufferedImage   image   = new BufferedImage(width, height, BufferedImage.TYPE_3BYTE_BGR);
042       Graphics2D      g       = image.createGraphics();
043       g.setColor(new Color(192,240,255));//Micaps1.0的背景色
044       g.fillRect(0, 0, width, height);//背景色填充
045 
046       Linear      coordinate  = new Linear(110.0, 35.0,width/2, height/2, 15, 15);//线性投影
047       Diamond09.drawBorderline(g, Color.gray, coordinate, "/path/to/ProvinceMap.dat");//画省界、国界、洲界
048       coordinate.drawGridLine(g, null, Color.green, 10, 10);//画经纬线
049 
050       delaunayT.draw(g, Color.blue, coordinate);//显示Delaunay网格
051 
052       //输出到浏览器
053       ServletOutputStream sos     = response.getOutputStream();
054       JPEGImageEncoder    encoder = JPEGCodec.createJPEGEncoder(sos);
055 
056       //以下三行改进图片质量
057       JPEGEncodeParam param = encoder.getDefaultJPEGEncodeParam(image);
058       param.setQuality(1.0f, false);
059       encoder.setJPEGEncodeParam(param);
060 
061       encoder.encode(image);
062  */
063 
064  //Credit to Paul Bourke (pbourke@swin.edu.au) for the original Fortran 77 Program :))
065  //Conversion to Visual Basic by EluZioN (EluZioN@casesladder.com)
066  //Conversion from VB to Delphi6 by Dr Steve Evans (steve@lociuk.com)
067  //08Jul2007 Conversion from Delphi6 to Java 1.4 by LIU Zejun (BJ0773@gmail.com)
068  ///
069  //June 2002 Update by Dr Steve Evans (steve@lociuk.com): Heap memory allocation
070  //added to prevent stack overflow when MaxVertices and MaxTriangles are very large.
071  //Additional Updates in June 2002:
072  //Bug in InCircle function fixed. Radius r := Sqrt(rsqr).
073  //Check for duplicate points added when inserting new point.
074  //For speed, all points pre-sorted in x direction using quicksort algorithm and
075  //triangles flagged when no longer needed. The circumcircle centre and radius of
076  //the triangles are now stored to improve calculation time.
077  ///
078  //08 Jul 2007 Update by LIU Zejun (BJ0773@gmail.com)
079  //Rewrite the FUNCTION InCircle(rename to circumcircle) to improve calculation time.
080  //Rewrite the FUNCTION Triangulate(rename to build) by used another algorithm.
081  //Convert all data dimension to DYNAMIC(used java.util.Vector).
082  ///
083  //You can use this code however you like providing the above credits remain in tact
084 
085 
086  package  cma.common.isoline; //要用到cma.common.isoline.TriangleVertex
087 
088  import  java.io.*;
089  import  java.lang.*;
090  import  java.util.*;
091  import  java.awt.*;
092  import  java.awt.geom.*;
093  import  java.text.DecimalFormat;
094  import  java.awt.image.*;
095  import  com.sun.image.codec.jpeg.*;
096 
097  import  cma.common.projection.*; //仅在draw函数中用到
098 
099  public class  Delaunay  {
100       public       Point2D.Double []     points;      //离散点序列
101       //triangle circle radius 的长度一致
102       public       Vector              triangle;    //三角形序列，存储TriangleVertex对象(各顶点在points中的索引)
103       public       Vector              circle;      //三角形的外接圆圆心，存储Point2D.Double对象
104       public       Vector              radius;      //三角形的外接圆半径，存储Double对象
105       private      int                  nt;          //三角形总数
106       private      int                  np;          //离散点总数
107 
108  /**
109    * 功能：
110    *      构造函数
111    * 参数：
112    *      pts     - 离散点序列，一般为经纬度值，长度必须>=3
113    * 返回值：
114    *      是否成功
115    */
116       public  Delaunay ( Point2D.Double []  pts ) {
117           np      = pts.length;
118           points  =  new  Point2D.Double [ np+ 3 ] ; //最后三个数据存放超大三角形的三个顶点
119           for ( int  i= 0 ;i<np;i++ ) {
120               points [ i ]    =  new  Point2D.Double ( pts [ i ] .x, pts [ i ] .y ) ;
121           }
122           triangle    =  new  Vector () ;
123           circle      =  new  Vector () ;
124           radius      =  new  Vector () ;
125           nt          =  0 ; //
126       }
127 
128  /*
129  超大三角形(SupperTriangle)示意图
130  假设为经纬度坐标(若为屏幕坐标，则为倒三角形，不影响计算)
131                   E
132                   ^
133                  /|/
134                 / | /
135                /  |  /
136               /30 | 30/
137              /    |    /
138             /     |     /
139            /      |H     /
140          D+-------+-------+C
141          /|       |       |/
142         / |       |       | /
143        /  |       |       |  /
144       /   |  xmid o ymid  |   /
145      /    |       |       |    /
146     /     |       |       |     /
147    / 60   |       |       |   60 /
148  /-------+-------+-------+-------/
149  F      A               B       G
150  等腰三角形△EFG，o(xmid,ymid)为源数据的中心
151  xmid=(xmin+xmax)/2
152  ymid=(ymin+ymax)/2
153 
154  AB=BC=CD=DA=EH=dmax=Math.max(xmax-xmin, ymax-ymin)
155  FA=DH=CH=BG=AB/2
156  ∠EFG=60度
157 
158  F=(xmid-dmax, ymid-dmax/2)
159  G=(xmid+dmax, ymid-dmax/2)
160  E=(xmid, ymid+dmax/2+dmax)
161  */
162 
163  /**
164    * 功能：
165    *      构造一个SuperTriangle，使得所有离散点均落在该三角形内
166    * 参数：
167    *      无
168    * 返回值：
169    *      无
170    */
171       private  boolean  createEdge () {
172           double       dmax, xmin, ymin, xmax, ymax, xmid, ymid;
173           xmin        = points [ 0 ] .x;
174           ymin        = points [ 0 ] .y;
175           xmax        = points [ 0 ] .x;
176           ymax        = points [ 0 ] .y;
177           for ( int  i= 1 ;i<np;i++ ) {
178               xmin    = Math.min ( xmin, points [ i ] .x ) ;
179               ymin    = Math.min ( ymin, points [ i ] .y ) ;
180               xmax    = Math.max ( xmax, points [ i ] .x ) ;
181               ymax    = Math.max ( ymax, points [ i ] .y ) ;
182           }
183           xmin    = Math.floor ( xmin ) ;
184           xmax    = Math.ceil ( xmax ) ;
185           ymin    = Math.floor ( ymin ) ;
186           ymax    = Math.ceil ( ymax ) ;
187 
188           //源数据区域的中心
189           xmid    =  ( xmin+xmax ) / 2.0 ;
190           ymid    =  ( ymin+ymax ) / 2.0 ;
191 
192           //经向和纬向最大距离的极值dmax，以中心为(xmid,ymid)、边长为dmax的正方形作为源数据的范围
193           dmax    = Math.max ( xmax-xmin, ymax-ymin ) ;
194           dmax    = Math.ceil ( 1.1 *dmax ) ; //增加0.1*dmax可确保离散点不在超大三角形的边上
195 
196           //超大三角形的顶点
197           points [ np+ 0 ]     =  new  Point2D.Double ( xmid-dmax, ymid-dmax/ 2 ) ;
198           points [ np+ 1 ]     =  new  Point2D.Double ( xmid, ymid+dmax/ 2 +dmax ) ;
199           points [ np+ 2 ]     =  new  Point2D.Double ( xmid+dmax, ymid-dmax/ 2 ) ;
200 
201           //创建超大三角形
202           boolean  enabled = add ( np+ 0 , np+ 1 , np+ 2 ) ;
203           return ( enabled ) ;
204       }
205 
206  /**
207    * 功能：
208    *      删除与SupperTriangle各顶点相关的三角形，包括SupperTriangle
209    * 参数：
210    *      无
211    * 返回值：
212    *      无
213    */
214       private  void  deleteEdge () {
215           int  sA  = np +  0 ;
216           int  sB  = np +  1 ;
217           int  sC  = np +  2 ;
218           TriangleVertex  tv;
219           for ( int  i=nt- 1 ;i>= 0 ;i-- ) {
220               tv  =  ( TriangleVertex ) triangle.get ( i ) ;
221               if (  tv.exists ( sA )  || tv.exists ( sB )  || tv.exists ( sC ) ) { //某个顶点与超大三角形相同
222                   this .delete ( i ) ;
223               }
224           }
225       }
226 
227  /**
228    * 功能：求三角形的外接圆圆心位置
229    * 参数：
230    *      x1,y1   - 顶点A的位置
231    *      x2,y2   - 顶点B的位置
232    *      x3,y3   - 顶点C的位置
233    * 返回值：
234    *      外接圆的圆心位置
235    * 公式：
236    *      根据“外接圆心到三顶点的距离均相等”推导得到
237    * (1) r*r = (x-x1)*(x-x1)+(y-y1)*(y-y1) = x*x-2*x*x1+x1*x1+y*y-x*y*y1+y1*y1
238    * (2) r*r = (x-x2)*(x-x2)+(y-y2)*(y-y2) = x*x-2*x*x2+x2*x2+y*y-x*y*y2+y2*y2
239    * (3) r*r = (x-x3)*(x-x3)+(y-y3)*(y-y3) = x*x-2*x*x3+x3*x3+y*y-x*y*y3+y3*y3
240    * 分别相减，略掉x*x和y*y
241    * x=[(y2-y1)*(y3*y3-y1*y1+x3*x3-x1*x1)-(y3-y1)*(y2*y2-y1*y1+x2*x2-x1*x1)] / {2*(x3-x1)*(y2-y1)-2*((x2-x1)*(y3-y1)]}
242    * y=[(x2-x1)*(x3*x3-x1*x1+y3*y3-y1*y1)-(x3-x1)*(x2*x2-x1*x1+y2*y2-y1*y1)] / {2*(y3-y1)*(x2-x1)-2*((y2-y1)*(x3-x1)]}
243    */
244       public  Point2D.Double circumcircle (
245               double  x1,  double  y1,  double  x2,  double  y2,  double  x3,  double  y3 ) {
246           double   x1x1    = x1 * x1;
247           double   x2x2    = x2 * x2;
248           double   x3x3    = x3 * x3;
249           double   y1y1    = y1 * y1;
250           double   y2y2    = y2 * y2;
251           double   y3y3    = y3 * y3;
252           double   x2_x1   = x2 - x1;
253           double   x3_x1   = x3 - x1;
254           double   y2_y1   = y2 - y1;
255           double   y3_y1   = y3 - y1;
256           double   x = (( y2_y1 ) * ( y3y3 -y1y1 +x3x3 -x1x1 )  - ( y3_y1 ) * ( y2y2 -y1y1 +x2x2 -x1x1  ))  /  ( 2 * ( x3_x1 ) * ( y2_y1 ) - 2 * ( x2_x1 ) * ( y3_y1 )) ;
257           double   y = (( x2_x1 ) * ( x3x3 -x1x1 +y3y3 -y1y1 )  - ( x3_x1 ) * ( x2x2 -x1x1 +y2y2 -y1y1  ))  /  ( 2 * ( y3_y1 ) * ( x2_x1 ) - 2 * ( y2_y1 ) * ( x3_x1 )) ;
258           return (
259               new  Point2D.Double ( x, y )
260           ) ;
261       }
262 
263  /**
264    * 功能：
265    *      判断点相对于矢线a->b的位置
266    * 参数：
267    *      x0,y0   - 要判断的点
268    *      x1,y1   - 线段端点a的位置
269    *      x2,y2   - 线段端点b的位置
270    * 返回值：
271    *      -1  - 在矢线a->b的左侧
272    *      0   - 在矢线a->b上(或延长线上)
273    *      1   - 在矢线a->b的右侧
274    */
275       public  int  position (
276               double  x0,  double  y0,
277               double  x1,  double  y1,
278               double  x2,  double  y2 ) {
279           double   xab   =    ( x0-x1 ) * ( y2-y1 ) - ( y0-y1 ) * ( x2-x1 ) ;
280           return ( xab <  0.0  ? - 1  : xab ==  0.0  ?  0  :  1 ) ;
281       }
282 
283  /**
284    * 功能：
285    *      判断AB与CD是否相交
286    * 参数：
287    *      xa,ya   - 线段AB的端点a的位置
288    *      xb,yb   - 线段AB的端点b的位置
289    *      xc,yc   - 线段CD的端点c的位置
290    *      xd,yd   - 线段CD的端点d的位置
291    * 返回值：
292    *      true    - 相交
293    *      false   - 相交
294    */
295       public  boolean  cross (
296               double  xa,  double  ya,
297               double  xb,  double  yb,
298               double  xc,  double  yc,
299               double  xd,  double  yd ) {
300           int  a_cd    =  this .position ( xa, ya, xc, yc, xd, yd ) ;
301           int  b_cd    =  this .position ( xb, yb, xc, yc, xd, yd ) ;
302           return ( 1 !=a_cd*b_cd ) ; //A与B在CD的同侧，则a_cd与b_cd同号，异侧则异号，某点在CD加上，则乘积为0
303       }
304 
305  /**
306    * 功能：
307    *      获得外接圆包含指定点的三角形序列中的索引值列表
308    * 参数：
309    *      p       - 待检查的点
310    * 返回值：
311    *      序号列表
312    */
313       private  int []  contains ( Point2D.Double p ) {
314           if (  nt ==  0  ) {
315               return ( null ) ;
316           }
317           Vector          lists   =  new  Vector () ;
318           Point2D.Double  pc;  //三角形的外接圆圆心位置
319           double           r;   //三角形的外接圆半径
320           double           pr2;
321           int []            res     =  new  int [ nt ] ;
322           int              count   =  0 ;
323           Arrays.fill ( res, - 1 ) ;
324           for ( int  i= 1 ;i<nt;i++ ) { //[0]是超大三角形，
325               r   =  (( Double ) radius.get ( i )) .doubleValue () ;
326               pc  =  ( Point2D.Double ) circle.get ( i ) ;
327               pr2 =  ( p.x - pc.x ) * ( p.x - pc.x ) + ( p.y - pc.y ) * ( p.y - pc.y ) ;
328               //与外接圆圆心的距离不超过外接圆的半径
329               if (  pr2 <= r*r  ) {
330                   res [ count++ ]     = i;
331               }
332           }
333           if (  count ==  0  ) {
334               return ( null ) ;
335           }
336           int []    results =  new  int [ count ] ;
337           for ( int  i= 0 ;i<count;i++ ) {
338               results [ i ]   = res [ i ] ;
339           }
340           return ( results ) ;
341       }
342 
343  /**
344    * 功能：
345    *      清空创建的三角形
346    * 参数：
347    *      无
348    * 返回值：
349    *      无
350    */
351       public  void  clear () {
352           triangle.clear () ;
353           circle.clear () ;
354           radius.clear () ;
355           nt  =  0 ;
356       }
357 
358  /**
359    * 功能：
360    *      删除指定的三角形
361    * 参数：
362    *      index   - 三角形索引值
363    * 返回值：
364    *      无
365    */
366       public  void  remove ( int  index ) {
367           if (  index >  0  && index < nt  ) { //不允许从类外删除SuperTriangle
368               delete ( index ) ;
369           }
370       }
371 
372  /**
373    * 功能：
374    *      删除指定的三角形
375    * 参数：
376    *      index   - 三角形索引值
377    * 返回值：
378    *      无
379    */
380       private  void  delete ( int  index ) {
381           if (  index >=  0  && index < nt  ) {
382               TriangleVertex  tv  =  ( TriangleVertex ) triangle.get ( index ) ;
383  //          System.out.println("Delaunay.java # 234 : remove T" + String.valueOf(index) + "=" + String.valueOf(tv.A) + "," + String.valueOf(tv.B) + "," + String.valueOf(tv.C));
384               triangle.remove ( index ) ;
385               circle.remove ( index ) ;
386               radius.remove ( index ) ;
387               nt--;
388           }
389       }
390 
391  /**
392    * 功能：
393    *      删除指定的三角形
394    * 参数：
395    *      index   - 三角形索引值
396    * 返回值：
397    *      无
398    */
399       public  void  remove ( int []  index ) {
400           if (  null  == index || index.length ==  0  ) {
401               return ;
402           }
403           int []    idx =  new  int [ index.length ] ;
404           for ( int  i= 0 ;i<index.length;i++ ) {
405               idx [ i ]   = index [ i ] ;
406           }
407           Arrays.sort ( idx ) ;
408           for ( int  i=idx.length- 1 ;i>= 0 ;i-- ) {
409               this .remove ( idx [ i ]) ;
410           }
411       }
412 
413  /** 注：此方法已废弃
414    * 功能：
415    *      根据index点与△abc的位置关系确定增加的三角形()
416    * 参数：
417    *      index   - 新点的索引
418    *      abc     - 已知的三角形
419    * 返回值：
420    *      能否成功创建
421    */
422       public  boolean  add ( int  index, TriangleVertex abc ) {
423           //index点相对于△abc三边的位置
424           int  posAB   =  this .position (
425                           points [ index ] .x, points [ index ] .y,
426                           points [ abc.A ] .x, points [ abc.A ] .y,
427                           points [ abc.B ] .x, points [ abc.B ] .y
428                       ) ;
429           int  posBC   =  this .position (
430                           points [ index ] .x, points [ index ] .y,
431                           points [ abc.B ] .x, points [ abc.B ] .y,
432                           points [ abc.C ] .x, points [ abc.C ] .y
433                       ) ;
434           int  posCA   =  this .position (
435                           points [ index ] .x, points [ index ] .y,
436                           points [ abc.C ] .x, points [ abc.C ] .y,
437                           points [ abc.A ] .x, points [ abc.A ] .y
438                       ) ;
439           boolean  inTriangle  = //全部在A->B->C->A的同一侧，说明在三角形内
440                               ( posAB == - 1  && posBC == - 1  && posCA == - 1 )  ||
441                               ( posAB ==   1  && posBC ==   1  && posCA ==   1 ) ;
442           boolean  enabled     =  false ;
443           if (  inTriangle  ) { //在△abc内
444  //          System.out.println("Delaunay.java # 295 : P" + String.valueOf(index) + " in T(" + String.valueOf(abc.A) + "," + String.valueOf(abc.B) + "," + String.valueOf(abc.C) + ")");
445               enabled =
446                           this .add ( index, abc.A, abc.B )  &&
447                           this .add ( index, abc.B, abc.C )  &&
448                           this .add ( index, abc.C, abc.A ) ;
449           }
450           else if (  posBC == posCA  ) { //在AB边上或AB边外
451  //          System.out.println("Delaunay.java # 302 : " + String.valueOf(nt) + "=" + String.valueOf(index) + "," + String.valueOf(abc.A) + "," + String.valueOf(abc.B));
452               enabled =  this .add ( index, abc.A, abc.B ) ;
453           }
454           else if (  posAB == posCA  ) { //在BC边上或BC边外
455  //          System.out.println("Delaunay.java # 306 : " + String.valueOf(nt) + "=" + String.valueOf(index) + "," + String.valueOf(abc.B) + "," + String.valueOf(abc.C));
456               enabled =  this .add ( index, abc.B, abc.C ) ;
457           }
458           else if (  posAB == posBC  ) { //在CA边上或CA边外
459  //          System.out.println("Delaunay.java # 310 : T" + String.valueOf(nt) + "=" + String.valueOf(index) + "," + String.valueOf(abc.C) + "," + String.valueOf(abc.A));
460               enabled =  this .add ( index, abc.C, abc.A ) ;
461           }
462           else  {
463  //          System.out.println("Delaunay.java # 314 : none");
464               enabled =  false ;
465           }
466           return ( enabled ) ;
467       }
468 
469  /**
470    * 功能：
471    *      增加一个三角形，并计算其外接圆圆心和半径
472    * 参数：
473    *      a,b,c   - 顶点的索引(必须保证是一个有效的三角形)
474    * 返回值：
475    *      能否成功创建
476    */
477       public  boolean  add ( int  a,  int  b,  int  c ) {
478           boolean  enabled     =  false ;
479           //创建一个三角形
480           TriangleVertex  tv  =  new  TriangleVertex () ;
481           enabled = tv.reset ( a, b, c ) ;
482           if (  !enabled  ) {
483               return ( false ) ;
484           }
485           int  ntSave  = nt; //记录现有的项数
486           try  {
487               triangle.add ( tv ) ;
488 
489               //三角形的外接圆圆心
490               Point2D.Double  p   =  this .circumcircle (
491                                           points [ a ] .x, points [ a ] .y,
492                                           points [ b ] .x, points [ b ] .y,
493                                           points [ c ] .x, points [ c ] .y ) ;
494               circle.add ( p ) ;
495 
496               //三角形的外接圆半径
497               Double          d   =  new  Double (
498                                       Math.sqrt (
499                                           ( points [ a ] .x-p.x ) * ( points [ a ] .x-p.x )  +
500                                           ( points [ a ] .y-p.y ) * ( points [ a ] .y-p.y )
501                                       )
502                                   ) ;
503               radius.add ( d ) ;
504  //          System.out.println("Delaunay.java # 355 : T" + String.valueOf(nt) + "=" + String.valueOf(tv.A) + "," + String.valueOf(tv.B) + "," + String.valueOf(tv.C));
505               nt++;
506               return ( true ) ;
507           }
508           catch ( Exception ex ) { //出错时删除已经增加的项
509               for ( int  i=triangle.size () - 1 ;i>ntSave;i-- ) {
510                   triangle.remove ( i ) ;
511               }
512               for ( int  i=circle.size () - 1 ;i>ntSave;i-- ) {
513                   circle.remove ( i ) ;
514               }
515               for ( int  i=radius.size () - 1 ;i>ntSave;i-- ) {
516                   radius.remove ( i ) ;
517               }
518               nt  = ntSave;
519               System.out.println ( ex.getMessage ()) ;
520               ex.printStackTrace () ;
521               return ( false ) ;
522           }
523       }
524 
525  /**
526    * 功能：
527    *      增加一系列三角形，并计算外接圆圆心和半径
528    * 参数：
529    *      index   - 离散点索引
530    *      lists   - 外接圆包含该离散点的三角形序列
531    * 返回值：
532    *      无
533    */
534       public  void  add ( int  index,  int []  lists ) {
535           if (  null  == lists  ) {
536               return ;
537           }
538  /*      if( lists.length == 1 ) {
539               this.add(index, (TriangleVertex)triangle.get(lists[0]));
540               return;
541           }
542  */       int []    idx     =  new  int [ lists.length ] ;
543           for ( int  i= 0 ;i<lists.length;i++ ) {
544               idx [ i ]   = lists [ i ] ;
545           }
546           Arrays.sort ( idx ) ;
547           TriangleVertex []     tvs         =  new  TriangleVertex [ idx.length ] ;
548           boolean [][]          existsSize  =  new  boolean [ 3 ][ idx.length ] ; //判断是否为共有边
549           Arrays.fill ( existsSize [ 0 ] ,  false ) ; //AB边列表
550           Arrays.fill ( existsSize [ 1 ] ,  false ) ; //BC边列表
551           Arrays.fill ( existsSize [ 2 ] ,  false ) ; //CA边列表
552           for ( int  i= 0 ;i<idx.length;i++ ) {
553               tvs [ i ]   =  ( TriangleVertex ) triangle.get ( idx [ i ]) ;
554           }
555           for ( int  i= 0 ;i<tvs.length;i++ ) { //在三角形序列中查找非共有的边，即序列构成的区域的边界
556               for ( int  j= 0 ;j<tvs.length;j++ ) {
557                   if (  i != j  ) {
558                       existsSize [ 0 ][ i ]     = existsSize [ 0 ][ i ]  || tvs [ j ] .exists ( tvs [ i ] .A, tvs [ i ] .B ) ; //其它三角形是否也存在AB边
559                       existsSize [ 1 ][ i ]     = existsSize [ 1 ][ i ]  || tvs [ j ] .exists ( tvs [ i ] .B, tvs [ i ] .C ) ; //其它三角形是否也存在BC边
560                       existsSize [ 2 ][ i ]     = existsSize [ 2 ][ i ]  || tvs [ j ] .exists ( tvs [ i ] .C, tvs [ i ] .A ) ; //其它三角形是否也存在CA边
561                   }
562               }
563           }
564           int  pos1, pos2;
565           for ( int  i= 0 ;i<tvs.length;i++ ) { //非共有的边，与P(index)点构成新三角形
566               if (      !existsSize [ 0 ][ i ]  && //非共有的边
567                       ! this .cross ( //PC与AB不相交
568                           points [ index ] .x, points [ index ] .y,
569                           points [ tvs [ i ] .C ] .x, points [ tvs [ i ] .C ] .y,
570                           points [ tvs [ i ] .A ] .x, points [ tvs [ i ] .A ] .y,
571                           points [ tvs [ i ] .B ] .x, points [ tvs [ i ] .B ] .y
572                       )
573                   ) {
574                   this .add ( index, tvs [ i ] .A, tvs [ i ] .B ) ; //构成△PAB
575               }
576               if (  !existsSize [ 1 ][ i ]   && //非共有的边
577                       ! this .cross ( //PA与BC不相交
578                           points [ index ] .x, points [ index ] .y,
579                           points [ tvs [ i ] .A ] .x, points [ tvs [ i ] .A ] .y,
580                           points [ tvs [ i ] .B ] .x, points [ tvs [ i ] .B ] .y,
581                           points [ tvs [ i ] .C ] .x, points [ tvs [ i ] .C ] .y
582                       )
583                   ) {
584                   this .add ( index, tvs [ i ] .B, tvs [ i ] .C ) ; //构成△PBC
585               }
586               if (  !existsSize [ 2 ][ i ]   && //非共有的边
587                       ! this .cross ( //PB与CA不相交
588                           points [ index ] .x, points [ index ] .y,
589                           points [ tvs [ i ] .B ] .x, points [ tvs [ i ] .B ] .y,
590                           points [ tvs [ i ] .C ] .x, points [ tvs [ i ] .C ] .y,
591                           points [ tvs [ i ] .A ] .x, points [ tvs [ i ] .A ] .y
592                       )
593                   ) {
594                   this .add ( index, tvs [ i ] .C, tvs [ i ] .A ) ; //构成△PCA
595               }
596           }
597       }
598 
599  /**
600    * 功能：
601    *      建立三角形网格
602    * 参数：
603    *      无
604    * 返回值：
605    *      成功创建的三角形个数
606    */
607       public  int  build () {
608           //清空已经存在的三角形网格
609           this .clear () ;
610 
611           if (  !createEdge ()    || //创建超大三角形
612               ! this .add ( 0 , np+ 0 , np+ 1 )     || //第一个离散点与超大三角形构成三个小三角形
613               ! this .add ( 0 , np+ 1 , np+ 2 )     ||
614               ! this .add ( 0 , np+ 0 , np+ 2 )    ) {
615               return ( 0 ) ;
616           }
617           for ( int  i= 1 ;i<np;i++ ) { //对所有离散点循环([0]已经使用)
618               int []    lists   =  this .contains ( points [ i ]) ; //获得外接圆包含该点的三角形序列
619               this .add ( i, lists ) ; //添加三角形
620               this .remove ( lists ) ;
621           }
622           deleteEdge () ;
623           return ( nt ) ;
624       }
625 
626  /**
627    * 功能：
628    *      显示三角形网格
629    * 参数：
630    *      g       - 图形设备
631    *      c       - 颜色
632    *      crd     - 坐标投影对象
633    * 返回值：
634    *      无
635    */
636       public  void  draw ( Graphics2D g, Color c, Coordinate crd ) {
637           Color   saveColor   = g.getColor () ;
638           g.setColor ( c ) ;
639           TriangleVertex  tv;
640           Point   p1, p2, p3 /*, p0*/ ;
641           int      r;
642           for ( int  i= 0 ;i<nt;i++ ) {
643               tv  =  ( TriangleVertex ) triangle.get ( i ) ;
644               p1  = crd.getPosition ( points [ tv.A ] .x, points [ tv.A ] .y ) ;
645               p2  = crd.getPosition ( points [ tv.B ] .x, points [ tv.B ] .y ) ;
646               p3  = crd.getPosition ( points [ tv.C ] .x, points [ tv.C ] .y ) ;
647               g.drawLine ( p1.x, p1.y, p2.x, p2.y ) ;
648               g.drawLine ( p2.x, p2.y, p3.x, p3.y ) ;
649               g.drawLine ( p3.x, p3.y, p1.x, p1.y ) ;
650           }
651           for ( int  i= 0 ;i<np;i++ ) {
652               p1  = crd.getPosition ( points [ i ] .x, points [ i ] .y ) ;
653               g.setColor ( Color.red ) ;
654               g.fillRect ( p1.x- 1 , p1.y- 1 ,  3 ,  3 ) ;
655               g.setColor ( Color.blue ) ;
656               g.drawString ( String.valueOf ( i ) , p1.x, p1.y ) ;
657           }
658           g.setColor ( saveColor ) ;
659       }
660 
661  }

 

---

TriangleVertex.java

---

001  /*
002 
003     三角形顶点
004 
005         PACKAGE: cma.common.isoline
006        FILENAME: TriangleVertex.java
007        LANGUAGE: Java2 v1.4
008        ORIGINAL: none
009     DESCRIPTION: Triangle Vertex Index
010          CREATE: 2007-07-06 23:09:14
011          UPDATE: 2007-07-06
012          AUTHOR: 刘泽军 (BJ0773@gmail.com)
013                  广西气象减灾研究所
014                  Guangxi Institude of Meteorology and Disaster-reducing Research(GIMDR)
015 
016         COMPILE: javac TriangleVertex.java
017 
018  */
019 
020  package  cma.common.isoline;
021 
022  import  java.lang.Math.*;
023 
024  public class  TriangleVertex  {
025 
026       //以从小到大的顺序存放 A < B < C
027       public  int       A;       //顶点A的在离散点序列中的索引
028       public  int       B;       //顶点B的在离散点序列中的索引
029       public  int       C;       //顶点C的在离散点序列中的索引
030 
031  /**
032    * 功能：
033    *      重新设定
034    * 参数：
035    *      i,j,k   - 三顶点在离散点序列中的索引
036    * 返回值：
037    *      无
038    */
039       public  boolean  reset ( int  i,  int  j,  int  k ) {
040           if (      i <   0  || j <   0  || k <   0  || //顶点必须是正确的索引值
041                   i == j || j == k || k == i  ) { //顶点不能相同
042               A   = - 1 ;
043               B   = - 1 ;
044               C   = - 1 ;
045               return ( false ) ;
046           }
047           else  {
048               A   = Math.min ( Math.min ( i, j ) , k ) ;
049               C   = Math.max ( Math.max ( i, j ) , k ) ;
050               B   =
051                       i != A && i != C ? i :
052                       j != A && j != C ? j :
053                       k; //k != A && k != C ? k : -1;
054  /*          if( B == -1 ) {
055                   A   = -1;
056                   C   = -1;
057               }
058  */
059               return ( true ) ;
060           }
061       }
062 
063  /**
064    * 功能：
065    *      构造函数
066    * 参数：
067    *      无
068    * 返回值：
069    *      无
070    */
071       public  TriangleVertex () {
072           A   = - 1 ;
073           B   = - 1 ;
074           C   = - 1 ;
075       }
076 
077  /**
078    * 功能：
079    *      构造函数
080    * 参数：
081    *      i,j,k   - 三顶点的索引
082    * 返回值：
083    *      无
084    */
085       public  TriangleVertex ( int  i,  int  j,  int  k ) {
086           reset ( i, j, k ) ;
087       }
088 
089  /**
090    * 功能：
091    *      判断三角形是否等价于指定的三角形
092    * 备注：
093    *      两个均由(-1, -1, -1)构成的三角形是不相等的，因为不是一个有效的三角形
094    * 参数：
095    *      a,b     - 两个顶点的索引值
096    * 返回值：
097    *      true    - 存在
098    *      false   - 不存在
099    */
100       public  boolean  equals ( int  a,  int  b,  int  c ) {
101           return ( exists ( a, b, c )) ;
102       }
103 
104  /**
105    * 功能：
106    *      判断三角形是否有三个顶点与指定的参数相同，即两个三角形相等
107    * 参数：
108    *      a,b,c   - 三个顶点的索引值
109    * 返回值：
110    *      true    - 存在
111    *      false   - 不存在
112    */
113       public  boolean  exists ( int  a,  int  b,  int  c ) {
114           return (
115               a != b && b != c && c != a &&
116               exists ( a )  &&
117               exists ( b )  &&
118               exists ( c )
119           ) ;
120       }
121 
122  /**
123    * 功能：
124    *      判断三角形是否有两个顶点与指定的参数相同
125    * 参数：
126    *      a,b     - 两个顶点的索引值
127    * 返回值：
128    *      true    - 存在
129    *      false   - 不存在
130    */
131       public  boolean  exists ( int  a,  int  b ) {
132           return ( a != b && exists ( a )  && exists ( b )) ;
133       }
134 
135  /**
136    * 功能：
137    *      判断三角形是否存在指定的顶点
138    * 参数：
139    *      a       - 顶点的索引值
140    * 返回值：
141    *      true    - 存在
142    *      false   - 不存在
143    */
144       public  boolean  exists ( int  a ) {
145           return (
146               a >=  0  &&  ( a == A || a == B || a == C )
147           ) ;
148       }
149 
150  }