构造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: 刘泽军 ([e-mail])
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 ([e-mail]) for the original Fortran 77 Program :))
065 //Conversion to Visual Basic by EluZioN ([e-mail])
066 //Conversion from VB to Delphi6 by Dr Steve Evans ([e-mail])
067 //08Jul2007 Conversion from Delphi6 to Java 1.4 by LIU Zejun ([e-mail])
068 ///////////////////////////////////////////////////////////////////////////////
069 //June 2002 Update by Dr Steve Evans ([e-mail]): 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 ([e-mail])
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: 刘泽军 ([e-mail])
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 }