[从头读历史] 第306节 星球战争 BC2799 至 BC2700（公元前28世纪）

[工程师阿伟]正在和[机器小伟]一起研究[星球战争 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;">	//delaunayif (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 listinitialize the triangle listdetermine the supertriangleadd supertriangle vertices to the end of the vertex listadd the supertriangle to the triangle listfor each sample point in the vertex listinitialize the edge bufferfor each triangle currently in the triangle listcalculate the triangle circumcircle center and radiusif the point lies in the triangle circumcircle thenadd the three triangle edges to the edge bufferremove the triangle from the triangle listendifendfordelete all doubly specified edges from the edge bufferthis leaves the edges of the enclosing polygon onlyadd to the triangle list all triangles formed between the point and the edges of the enclosing polygonendforremove any triangles from the triangle list that use the supertriangle verticesremove 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进行合并除去与超级三角形有关的三角形
endWatson算法的基本步骤是：
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>

本节到此结束，欲知后事如何，请看下回分解。