package esurfer;
import java.awt.geom.Point2D;
import java.io.BufferedReader;
import java.io.BufferedWriter;
import java.io.FileWriter;
import java.io.IOException;
import java.io.Reader;
import java.io.UnsupportedEncodingException;
import java.util.*;
/*
* ported from p bourke's triangulate.c
* http://astronomy.swin.edu.au/~pbourke/terrain/triangulate/triangulate.c
*
* fjenett, 20th february 2005, offenbach-germany. contact:
* http://www.florianjenett.de/
*
* run like this: javac *.java java triangulate
*
* to view the output: http://processing.org/
*
*/
class Triangle {
TQEPoint vertex1, vertex2, vertex3;
Triangle(TQEPoint vertex1, TQEPoint vertex2, TQEPoint vertex3) {
this.vertex1 = vertex1;
this.vertex2 = vertex2;
this.vertex3 = vertex3;
}
}
public class Triangulate {
//public static double EPSILON = 0.000001;
private static boolean DEBUG = false;
private static float dMax = 0;
static ArrayList Triangulate(int nv, TPoint2D points[], ArrayList borderList) {
ArrayList triangleList = new ArrayList();
return triangleList;
}
/**
* 追踪离散点边界
* @param points 离散点数组
* @return 边界点集链表，每项为GeoPoint
*/
public static ArrayList traceBoundary(TPoint2D points[]) {
ArrayList borderList = new ArrayList(); // 边界链表
int np = points.length;
Point2D[] pts2d = new Point2D.Double[np];
for (int i = 0; i < np; i++) {
pts2d[i] = new Point2D.Double(points[i].x, points[i].y);
}
Point2D[] brdVertex = QuickHull.computeConvexHull(pts2d);
np = brdVertex.length;
for (int i = 0; i < np; i++) {
borderList.add(new GeoPoint(brdVertex[i].getX(), brdVertex[i].getY()));
}
return borderList;
}
/**
* 追踪三角网及边界
* @param points 离散点数组
* @return 存储TIN三角形数组及边界数组的哈希表，取得方法:
* 三角形： (ArrayList)map.get("triangle")，链表元素为 Triangle 类型
* 边界： (ArrayList)map.get("border")，链表元素为 GeoPoint 类型
**/
public static HashMap createTIN(ArrayList pointsList, int borderParam) {
// 找出x最大值及y最大值
double xMax = Double.MIN_VALUE;
double yMax = Double.MIN_VALUE;
int numOfPoints = pointsList.size();
for (int i = 0; i < numOfPoints; i++) {
TPoint2D p = (TPoint2D)pointsList.get(i);
if (p.x > xMax) {
xMax = p.x;
}
if (p.y > yMax) {
yMax = p.y;
}
}
HashMap mapTriBrd = new HashMap(); // 用于存储三角表链表triangleList及边界链表borderList
TPoint2D[] tempPoints = new TPoint2D[pointsList.size()];
tempPoints = (TPoint2D[])pointsList.toArray(tempPoints);
TPoint2D[] points = new TPoint2D[pointsList.size()];
points = (TPoint2D[])pointsList.toArray(tempPoints);
// 追踪边界
ArrayList borderList = traceBoundary(tempPoints); // 初始边界链表(凸壳)
GeoPoint[] brdPoints = new GeoPoint[borderList.size()];
brdPoints = (GeoPoint[])borderList.toArray(brdPoints);
// 计算凸壳面积
float hullArea = (float)calculatePolyArea(brdPoints);
// 计算离散点密度
float density = points.length / hullArea;
// 估算相邻离散点平均距离
float avgDistance = (float)(1 / Math.sqrt(density));
// 边界线段最大长度
dMax = borderParam * avgDistance;
// 边界收缩
boolean hasNewBrdPnt;
GeoPoint startPnt = null, endPnt = null;
do {
hasNewBrdPnt = false; // 判断是否还有新的边界点
for (int i = 0; i < borderList.size(); i++) {
startPnt = (GeoPoint)borderList.get(i);
if ( i == (borderList.size() - 1) ) {
endPnt = (GeoPoint)borderList.get(0);
} else {
endPnt = (GeoPoint)borderList.get(i + 1);
}
if (pnt2PntDistance(startPnt, endPnt) > dMax) { // 边界线段长度大于 dMax
// 寻找左侧的一个最近点，加入边界链表
float temp = 0;
int nearPntIdx = -1;
for (int j = 0; j < points.length; j++) {
// 判断是否边界点
boolean isBrdPoint = false;
for (int k = 0; k < borderList.size(); k++) {
GeoPoint brdPnt = (GeoPoint)borderList.get(k);
if (points[j].x == brdPnt.x && points[j].y == brdPnt.y) {
isBrdPoint = true;
break;
}
}
if (isBrdPoint) continue;
// 在边界线左侧找出第三个点的索引，使角v1v3v2最大(用余弦定理求角v1v3v2)。
float pnt2Line = posOfLine(points[j], startPnt, endPnt);
if (pnt2Line < 0) { // 在左侧
float squareEdgeV1V2 = (float)( (startPnt.x - endPnt.x) * (startPnt.x - endPnt.x) + (startPnt.y - endPnt.y) * (startPnt.y - endPnt.y) );
float squareEdgeV2V3 = (float)( (endPnt.x - points[j].x) * (endPnt.x - points[j].x) + (endPnt.y - points[j].y) * (endPnt.y - points[j].y) );
float squareEdgeV3V1 = (float)( (startPnt.x - points[j].x) * (startPnt.x - points[j].x) + (startPnt.y - points[j].y) * (startPnt.y - points[j].y) );
float cosV132 = (float)( 1 - (squareEdgeV2V3 + squareEdgeV3V1 - squareEdgeV1V2) / (2 * Math.sqrt(squareEdgeV2V3 * squareEdgeV3V1)) );
if (cosV132 > temp + 0.00001) {
temp = cosV132;
nearPntIdx = j;
}
}
}
if (nearPntIdx != -1) {
hasNewBrdPnt = true;
// 将最近点加入边界链表
borderList.add(i + 1, new GeoPoint(points[nearPntIdx].x, points[nearPntIdx].y));
i += 2;
}
}
}
if (!hasNewBrdPnt) break;
} while (true);
// 从离散点中分离出边界点和边界内其他离散点
int i = 0;
ArrayList brdZList = new ArrayList();
while (borderList.size() > 0 && i < borderList.size()) {
GeoPoint p = (GeoPoint)borderList.get(i);
int j = 0;
boolean find = false;
while (pointsList.size() > 0 && j < pointsList.size()) {
TPoint2D p2d = (TPoint2D)pointsList.get(j);
if ( Math.abs(p.x - p2d.x) < Geometry.EPS && Math.abs(p.y - p2d.y) < Geometry.EPS) {
// 将带有Z值的点加入边界链表
brdZList.add(new TPoint2D(p2d.x, p2d.y, p2d.z));
// 将已经找到的边界点从原边界链表及离散点集中删除
borderList.remove(i);
borderList.trimToSize();
pointsList.remove(j);
pointsList.trimToSize();
find = true;
break;
}
j++;
}
if (!find) {
i++;
}
}
TPoint2D[] brdPnts = new TPoint2D[brdZList.size()];
brdPnts = (TPoint2D[])brdZList.toArray(brdPnts);
borderList.clear();
for (int k = 0; k < brdPnts.length; k++) {
borderList.add(new GeoPoint(brdPnts[k].x, brdPnts[k].y));
}
QuadMesh mesh = new QuadMesh();
评论0