标签：
剧情提要：[工程师阿伟]正在和[机器小伟]一起研究[星球战争 BC2799 至 BC2700（公元前28世纪）]。
这一百年，照样没发生什么有营养的战争，阿伟耐不住寂寞，就去看看神农氏在干嘛。
神农氏正对着一张草图发呆。
这是神马？阿伟问他。
“噢，我们不是占领了斧燧氏的地盘嘛，但是原来的土地也不能放弃啊，你来看，在地图上是这个样子。“
“嗯哼”。
”我在想怎样在这些地方之间走动更有效率一点，毕竟到处都是洪水和沼泽，能走的道也不是太多。“
”然后你就把它画成一团乱麻了？“阿伟问他。
”画着画着不就多了么。“
”这个简单，你看，我先画个简单点的给你看。“阿伟拿出了一张图。
然后像你那个，只要念一声咒语，啊呜。
这不久成了。
神农高兴起来了。
“不过，我可要批评你了，你怎么还不去打仗啊，你不打仗我写什么呀。”阿伟不客气地指出了主要问题。
“这个，我还要改进一样我们部落的生产工具，还有战争的武器啊，前不久刚得到一张设计图，就是不知道怎么样制造啊。“
这个，阿伟觉得应该先试一下，看能不能帮帮他。
一阵昏天黑地地运算过后，呯地一声响，从[机器小伟]的肚子里飘出来一张草图。
阿伟运起火眼金睛一看，这图缺陷颇多啊，看来算法还是不够威武，只能无奈地对神农说，看来你们的工具还得慢慢地来，战争的事情看来还得往后再安排啊。
好了，不管怎么说，阿伟还是决定给神农留下全套设计思路，希望他早日造出武器，去发动战争。
彩蛋：
<span style="font-size:18px;">/**
* @usage 分块治图取轮廓
* @author mw
* @date 2015年12月20日 星期日 10:46:21
* @param
* @return
*
*/
function Nexus() {
//方格大小范围，比如2*2, 10*10, 100*100个象素
this.range = 100;
//方格序数
this.xth = 0;
this.yth = 0;
//方格开始和结束的x,y象素点坐标
this.xBeg = 0;
this.xEnd = 0;
this.yBeg = 0;
this.yEnd = 0;
//候选点的x, y极值
this.xMax = 0;
this.xMin = 0;
this.yMax = 0;
this.yMin = 0;
//候选点个数和比例
this.count = 0;
this.percent = 0;
//候选点的x, y总值
this.xTotal = 0;
this.yTotal = 0;
//候选点族的中心
this.xCenter = 0;
this.yCenter = 0;
//候选点族的半径区域
this.r = 0;
//需要传入参数xth, yth, range
this.init = function(xth, yth, range) {
this.xth = xth;
this.yth = yth;
this.range = range;
//方格开始和结束的x,y象素点坐标
this.xBeg = 0 + this.xth * this.range;
this.xEnd = 0 + (this.xth+1) * this.range;
this.yBeg = 0 + this.yth * this.range;
this.yEnd = 0 + (this.yth+1) * this.range;
};
//传入7个参数
this.calc = function(count, xTotal, yTotal, xMax, yMax, xMin, yMin) {
this.count = count;
this.xTotal = xTotal;
this.yTotal = yTotal;
this.xMax = xMax;
this.xMin = xMin;
this.yMax = yMax;
this.yMin = yMin;
this.percent = this.count / (this.range * this.range);
//候选点族的中心
this.xCenter = this.count==0 ? (this.xBeg+this.xEnd)/2 : this.xTotal / this.count;
this.yCenter = this.count==0 ? (this.yBeg+this.yEnd)/2 : this.yTotal / this.count;
//候选点族的半径区域
this.r = Math.min(
Math.abs(this.xMax - this.xCenter),
Math.abs(this.xMin - this.xCenter),
Math.abs(this.yMax - this.yCenter),
Math.abs(this.yMin - this.yCenter),
this.range
);
};
}
/**
* @usage 压缩中心点数量
* @author mw
* @date 2015年12月23日 星期三 11:45:43
* @param
* @return
*
*/
//压缩中心点数量
function centerPoint(arr, tolerance) {
tolerance = tolerance ? tolerance : 50;
var array = new Array();
var indexArray = new Array();
var len = arr.length;
var xTotal=0, yTotal=0, count=0;
for (var i =0 ; i < len; i++) {
if (indexArray[i]) continue;
//第i号已处理
indexArray[i]=1;
xTotal = arr[i][2];
yTotal = arr[i][3];
count=1;
for (var j =i+1; j < len; j++) {
if (indexArray[j]) continue;
if (Math.abs(arr[i][0] - arr[j][0])<=1 &&
Math.abs(arr[i][1] - arr[j][1])<=1) {
if (Math.abs(arr[i][2]-arr[j][2])+
Math.abs(arr[i][3] - arr[j][3]) < tolerance) {
indexArray[j]=1;
xTotal += arr[j][2];
yTotal += arr[j][3];
count++;
}
}
}
array.push([arr[i][0], arr[i][1], Math.round(xTotal/count), Math.round(yTotal/count)]);
}
array.sort(function(a, b) {
if (a[1] < b[1]) return -1;
else if (a[1] > b[1]) {
return 1;
}
else {
if (a[0] < b[0]) return -1;
else if (a[0] > b[0]) return 1;
else return 0;
}
return 1;
});
var pointInfo = "$picDataArray = [";
len = array.length;
for (var i = 0; i < len; i++) {
pointInfo += "["+array[i][0]+", "+array[i][1]+", "+array[i][2]+", "
+array[i][3]+"], ";
}
pointInfo += "];";
document.body.appendChild(document.createTextNode(pointInfo));
return array;
}
/**
* @usage 按区格压缩分块点，依赖Nexus节点类
* @author mw
* @date 2015年12月23日 星期三 11:45:43
* @param
* @return
*
*/
//按区格压缩分块点
function step1() {
var arr = new Array();
var pointInfo="$picDataArray = [";
//图片
var image = new Image();
image.src = "./1.jpg";
//只处理这100*100个象素
var width = 700;
var height = 600;
//格子大小，行和列共有多少格
var range = 20;
var rows = height / range;
var cols = width / range;
//确定范围
var xBeg, xEnd, yBeg, yEnd;
//待计算参数
var count, xTotal, yTotal, xMin, yMin, xMax, yMax;
var gap = 50;
image.onload = function() {
plot.drawImage(image);
var imagedata = plot.getImageData(0, 0, width, height);
//计算
for (var i = 0; i < cols; i++) {
for (var j = 0; j < rows; j++) {
var nexus = new Nexus();
nexus.init(i, j, range);
//确定范围
xBeg = nexus.xBeg;
xEnd = nexus.xEnd;
yBeg = nexus.yBeg;
yEnd = nexus.yEnd;
//待计算参数
xMin = nexus.xMin;
xMax = nexus.xMax;
yMin = nexus.yMin;
yMax = nexus.yMax;
count = 0;
xTotal = 0;
yTotal = 0;
//水平方向找差异
for (var col = xBeg+1; col < xEnd; col++) {
for (var row = yBeg+2; row<yEnd-2; row++) {
//pos最小为1
pos =row * width + col;
R0 = imagedata.data[4 * (pos-1)];
R1 = imagedata.data[4 * pos];
G0 = imagedata.data[4 * (pos-1)+1];
G1 = imagedata.data[4 * pos+1];
B0 = imagedata.data[4 * (pos-1)+2]
B1 = imagedata.data[4 * pos + 2]
//简单容差判断
if (Math.abs(R1-R0) > gap ||
Math.abs(G1-G0)>gap ||
Math.abs(B1-B0)>gap){
count++;
xTotal += col;
yTotal += row;
if (xMin > col) xMin = col;
if (xMax < col) xMax = col;
if (yMin > row) yMin = row;
if (yMax < row) yMax = row;
}
}
}
//垂直方向找差异
for (var col = xBeg+2; col < xEnd-2; col++) {
for (var row = yBeg+1; row<yEnd; row++) {
//pos最小为第二行
pos =row * width + col;
R0 = imagedata.data[4 * (pos-width)];
R1 = imagedata.data[4 * pos];
G0 = imagedata.data[4 * (pos-width)+1];
G1 = imagedata.data[4 * pos+1];
B0 = imagedata.data[4 * (pos-width)+2];
B1 = imagedata.data[4 * pos + 2];
//简单容差判断
if (Math.abs(R1-R0) > gap ||
Math.abs(G1-G0)>gap ||
Math.abs(B1-B0)>gap) {
count++;
xTotal += col;
yTotal += row;
if (xMin > col) xMin = col;
if (xMax < col) xMax = col;
if (yMin > row) yMin = row;
if (yMax < row) yMax = row;
}
}
}
nexus.calc(count, xTotal, yTotal, xMax, yMax, xMin, yMin);
arr.push(nexus);
}
}
arr.sort(function(a, b) {
if (a[1] < b[1]) return -1;
else if (a[1] > b[1]) {
return 1;
}
else {
if (a[0] < b[0]) return -1;
else if (a[0] > b[0]) return 1;
else return 0;
}
return 1;
});
var nexus = new Nexus();
for (var i = 0; i < arr.length; i++) {
nexus = arr[i];
if (nexus.count > 10 && nexus.percent < 0.9) {
pointInfo += '[' + nexus.yth + ', ' + nexus.xth + ',' +
nexus.xCenter.toFixed(0) + ', ' + nexus.yCenter.toFixed(0)+'], ';
}
if (nexus.count > 10 && nexus.percent < 0.9) {
fillCircle(nexus.xCenter, nexus.yCenter, 4);
}
else {
fillCircle(nexus.xCenter, nexus.yCenter, 1);
}
}
pointInfo += '];';
var pointInfoNode = document.createTextNode(pointInfo);
document.body.appendChild(pointInfoNode);
}
}
//方便函数
//绘点并产生正确的点个数
function step2() {
var arr = new Array();
arr = centerPoint($picDataArray, 100);
var image = new Image();
image.src = "./1.jpg";
image.onload = function() {
plot.drawImage(image);
var len = arr.length;
plot.save()
.setFillStyle('blue');
for (var i = 0; i<len; i++) {
//去取图片边界上的分隔点
if (Math.abs(arr[i][2] - 350) <=320 && Math.abs(arr[i][3]-300)<=270) {
fillCircle(arr[i][2], arr[i][3], 10);
}
}
plot.restore();
}
} </span>
<span style="font-size:18px;"> //delaunay
if (1) {
var r = 20;
config.setSector(1,1,1,1);
config.graphPaper2D(0, 0, r);
config.axis2D(0, 0, 180);
var style = 'blue';
var scale = 1;
var showLable = 1;
var mapWidth = 700, mapHeight = 600;
var transform = new Transform();
$Verts = transform.translate(transform.flipX($Verts), -mapWidth/2, mapHeight/2);
var triTrans = [];
var triNum = $Triangles.length;
for (var i = 0; i < triNum; i++) {
triTrans.push(transform.translate(transform.flipX($Triangles[i]), -mapWidth/2, mapHeight/2));
}
for (var i = 0; i < triNum; i++) {
shape.strokeDraw([].concat(triTrans[i]), style, scale);
}
shape.pointDraw([].concat($Verts), 'red', scale, showLable);
}
if (0) {
step2();
}</span><span style="font-size:18px;">'''
#伪代码思路
subroutine triangulate
input : vertex list
output : triangle list
initialize the triangle list
determine the supertriangle
add supertriangle vertices to the end of the vertex list
add the supertriangle to the triangle list
for each sample point in the vertex list
initialize the edge buffer
for each triangle currently in the triangle list
calculate the triangle circumcircle center and radius
if the point lies in the triangle circumcircle then
add the three triangle edges to the edge buffer
remove the triangle from the triangle list
endif
endfor
delete all doubly specified edges from the edge buffer
this leaves the edges of the enclosing polygon only
add to the triangle list all triangles formed between the point
and the edges of the enclosing polygon
endfor
remove any triangles from the triangle list that use the supertriangle vertices
remove the supertriangle vertices from the vertex list
end
算法伪代码
input: 顶点列表(vertices) //vertices为外部生成的随机或乱序顶点列表
output:已确定的三角形列表(triangles)
初始化顶点列表
创建索引列表(indices = new Array(vertices.length)) //indices数组中的值为0,1,2,3,......,vertices.length-1
基于vertices中的顶点x坐标对indices进行sort //sort后的indices值顺序为顶点坐标x从小到大排序（也可对y坐标，本例中针对x坐标）
确定超级三角形
将超级三角形保存至未确定三角形列表（temp triangles）
将超级三角形push到triangles列表
遍历基于indices顺序的vertices中每一个点 //基于indices后，则顶点则是由x从小到大出现
初始化边缓存数组（edge buffer）
遍历temp triangles中的每一个三角形
计算该三角形的圆心和半径
如果该点在外接圆的右侧
则该三角形为Delaunay三角形，保存到triangles
并在temp里去除掉
跳过
如果该点在外接圆外（即也不是外接圆右侧）
则该三角形为不确定 //后面会在问题中讨论
跳过
如果该点在外接圆内
则该三角形不为Delaunay三角形
将三边保存至edge buffer
在temp中去除掉该三角形
对edge buffer进行去重
将edge buffer中的边与当前的点进行组合成若干三角形并保存至temp triangles中
将triangles与temp triangles进行合并
除去与超级三角形有关的三角形
end
Watson算法的基本步骤是：
1、构造一个超级三角形，包含所有散点，放入三角形链表。
2、将点集中的散点依次插入，在三角形链表中找出外接圆包含插入点的三角形（称为该点的影响三角形），删除影响三角形的公共边，将插入点同影响三角形的全部顶点连接起来，完成一个点在Delaunay三角形链表中的插入。
3、根据优化准则对局部新形成的三角形优化。将形成的三角形放入Delaunay三角形链表。
4、循环执行上述第2步，直到所有散点插入完毕。
'''
import geo;
###
# @usage Bowyer-Watson算法进行Delaunay三角剖分
# @author mw
# @date 2016年07月15日 星期五 10:31:36
# @param
# @return
#
###
class Delaunay():
#设置顶点
#vertices是[[x_0, y_0], [x_1, y_1], ...]顶点对格式
def setVertice(self, vertices):
return vertices;
def sortVerticebyX(self, vertices):
v_x = sorted(vertices, key = lambda a : (a[0], a[1]));
return v_x;
def sortVerticebyY(self, vertices):
v_y = sorted(vertices, key = lambda a : (a[1], a[0]));
return v_y;
#去除重复点
def removeDup(self, vertices):
v_new = [];
len_1 = len(vertices);
len_2 = 0;
for i in range(len_1):
len_2 = len(v_new);
if (len_2 < 1):
v_new.append(vertices[i]);
else:
for j in range(len_2):
if v_new[j] == vertices[i]:
break;
if (j >= len_2-1):
v_new.append(vertices[i]);
return v_new;
#计算边界
def calcBound(self, vertices):
len_ = len(vertices)
v_x = self.sortVerticebyX(vertices);
xMin = v_x[0][0];
xMax = v_x[len_-1][0];
v_y = self.sortVerticebyY(vertices);
yMin = v_y[0][1];
yMax = v_y[len_-1][1];
return [xMin, xMax, yMin, yMax];
#超级三角形
def superTri(self, vertices):
bound = self.calcBound(vertices);
xMin = bound[0]-10;
xMax = bound[1]+10;
yMin = bound[2]-10;
yMax = bound[3]+10;
xCenter = (xMin+xMax)/2;
yCenter = (yMin+yMax)/2;
xR = (xMax-xMin)/2;
yR = (yMax-yMin)/2;
xMin_ = xCenter-2*xR;
xMax_ = xCenter+2*xR;
yMin_ = yMin;
yMax_ = yMin + 4*yR;
return [[xMin_, yMin_], [xMax_, yMin_], [xCenter, yMax_]];
#计算剖分三角形
def calcDelaunayTri(self, vertices, mode = 1):
#移除重复点
vertices = self.removeDup(vertices);
#按X坐标由小到大排序
vertices = self.sortVerticebyX(vertices);
#顶点数量
vertNum = len(vertices);
#临时三角形存放处
tempTriArray = [];
#三角形存放处
triArray = [];
#边存放处
edgeArray = [];
supertri = self.superTri(vertices);
tempTriArray.append(supertri);
triArray.append(supertri);
for i in range(vertNum):
P0 = vertices[i];
tmpTriNum = len(tempTriArray);
print('顶点{0} --> {1}个临时三角形'.format(i, tmpTriNum));
edgeArray = [];
tmpTri = [];
for j in range(tmpTriNum):
P1, P2, P3 = tempTriArray[j][0], tempTriArray[j][1], tempTriArray[j][2];
#调用geo的circle方法
circleProp = geo.circle(P1, P2, P3);
#取得圆心和半径
P_center, R = circleProp[0], circleProp[1];
#print(P_center, 'R = ', R);
d = geo.distance2D(P0, P_center);
if (P0[0] > P_center[0]+R):
#对比点是在圆外右侧的三角形，已经得到晋级确认
triArray.append([P1, P2, P3]);
elif (d > R):
#不确定的三角形，不理它
tmpTri.append([P1, P2, P3]);
else:
#对比点在圆内，这个三角形被淘汰
edgeArray.append([P1, P2]);
edgeArray.append([P2, P3]);
edgeArray.append([P3, P1]);
edgeArray = self.removeDupEdge(edgeArray);
edges = len(edgeArray);
print('顶点{0} --> {1}条边'.format(i, edges));
for k in range(edges):
P1, P2 = edgeArray[k][0], edgeArray[k][1];
if (geo.pointInLine(P1, P2, P0) == False):
tmpTri.append([P0, P1, P2]);
#临时数组已经重新安排
tempTriArray = [];
tempTriArray = self.removeDupTri(tmpTri);
triArray += tempTriArray;
triArray = self.removeDupTri(triArray);
if (mode == 0):
return triArray;
else:
newTriArray = [];
triNum = len(triArray);
for i in range(triNum):
tri_ = triArray[i];
relate = False;
for j in range(3):
if relate == True:
break;
for k in range(3):
if tri_[j] == supertri[k]:
relate = True;
break;
if relate == False:
newTriArray.append(tri_);
return newTriArray;
#移除相同的三角形
def removeDupTri(self, triArray):
newTriArray = [];
for i in range(len(triArray)):
len_ = len(newTriArray);
if (len_ <= 0):
newTriArray.append(triArray[i]);
else:
for j in range(len_):
if self.judgeSameTri(newTriArray[j], triArray[i]) == True:
break;
if (j >= len_ -1):
newTriArray.append(triArray[i]);
return newTriArray;
#判断两个三角形相同
#三角形格式[P1, P2, P3], P是顶点
def judgeSameTri(self, tri_1, tri_2):
P_11, P_12, P_13, P_21, P_22, P_23 = tri_1[0], tri_1[1], tri_1[2], tri_2[0], tri_2[1], tri_2[2];
tri_1 = sorted(tri_1, key = lambda a:(a[0], a[1]));
tri_2 = sorted(tri_2, key = lambda a:(a[0], a[1]));
if (tri_1 == tri_2):
return True;
else:
return False;
#移除相同的边，本算法的去重指的是要成对的却除相同的边
#而不是只允许出现一次那种
def removeDupEdge(self, edgeArray):
newEdgeArray = [];
for i in range(len(edgeArray)):
len_ = len(newEdgeArray);
if (len_ <= 0):
newEdgeArray.append(edgeArray[i]);
else:
for j in range(len_):
if self.judgeSameEdge(newEdgeArray[j], edgeArray[i]) == True:
newEdgeArray = newEdgeArray[:j]+newEdgeArray[j+1:];
break;
if (j >= len_ -1):
newEdgeArray.append(edgeArray[i]);
return newEdgeArray;
#判断两条边相同
#边格式[P1, P2], P是顶点
def judgeSameEdge(self, edge_1, edge_2):
P_11, P_12, P_21, P_22 = edge_1[0], edge_1[1], edge_2[0], edge_2[1];
if (P_11 == P_21 and P_12 == P_22) or (P_11 == P_22 and P_12 == P_21):
return True;
else:
return False;
if __name__ == '__main__':
delaunay = Delaunay();
picDataArray = [[11, 0, 22, 241], [13, 1, 29, 279], [15, 1, 43, 317], [10, 2, 65, 203], [14, 3, 80, 289], [9, 4, 98, 186], [12, 5, 117, 257], [14, 5, 106, 285], [7, 6, 144, 149], [10, 7, 160, 219], [7, 8, 183, 140], [9, 9, 195, 206], [6, 10, 219, 141], [8, 11, 241, 162], [11, 11, 242, 237], [9, 13, 282, 195], [12, 13, 279, 252], [9, 15, 305, 203], [11, 15, 327, 234], [13, 16, 338, 278], [11, 17, 350, 232], [13, 18, 383, 275], [15, 18, 374, 314], [13, 20, 409, 284], [17, 21, 441, 353], [15, 23, 470, 324], [19, 24, 499, 393], [17, 25, 519, 356], [21, 26, 541, 436], [19, 27, 565, 392], [22, 28, 580, 459], [24, 29, 601, 480], [21, 31, 639, 437], [24, 31, 642, 500], [22, 33, 679, 464], [24, 33, 670, 499]];
verts = [];
mapWidth, mapHeight = 700/2, 600/2;
for i in range(len(picDataArray)):
x = picDataArray[i][2];
y = picDataArray[i][3];
if (abs(x-mapWidth)<=350 and abs(y-mapHeight)<=270):
verts.append([x, y]);
print('顶点集：', verts);
if (len(verts) >= 3):
a = delaunay.calcDelaunayTri(verts);
print('写入文件开始。>>>');
fout = open('output.txt', 'w');
print('共有顶点{0}个，共有三角形{1}个'.format(len(verts), len(a)));
s = '$Verts = ';
s += str(verts);
fout.write(s + '\n\n\n\n');
s = '$Triangles = ';
s += str(a);
fout.write(s+'\n');
fout.close();
print('写入文件完毕。');</span>
本节到此结束，欲知后事如何，请看下回分解。
[从头读历史] 第306节 星球战争 BC2799 至 BC2700（公元前28世纪）
标签：
原文地址：http://blog.csdn.net/mwsister/article/details/51924396
