JavaScript中点在3D凸多边形内部的判断
5.00/5 (2投票s)
在JavaScript中判断一个点是否在一个给定多边形顶点构成的3D凸多边形内部的算法。
引言
这是一个原始 C++ 算法的 JavaScript 版本，可以在这里找到。
背景
请参考原始 C++ 算法 这里。
Using the Code
该算法被封装到一个 JavaScript 库文件夹 GeoProc 中。 主测试程序 JSLASProc 从 LAS 文件读取点云数据，然后将所有内部点写入一个新的 LAS 文件。 可以使用免费软件 FugroViewer 显示 LAS 文件。
GeoProc 库文件夹中的每个 JavaScript 文件都有用法和测试部分。 要测试每个函数，请取消注释测试部分并在浏览器中打开 GeoProc.htm。
GeoPoint.js
Javascript 没有运算符重载，Add 方法用于模拟运算符重载。 有关其用法，请查看 GeoVector.js 中 Add 方法的引用。
// Constructor
function GeoPoint(x, y, z) // double x, y, z
{
this.x = x;
this.y = y;
this.z = z;
}
// public Instance method to simulate overloading binary operator add (GeoPoint + GeoPoint)
GeoPoint.prototype.Add = function(p) // GeoPoint p
{
return new GeoPoint(this.x + p.x, this.y + p.y, this.z + p.z);
};
GeoVector.js
Multiple 方法用于模拟 GeoVector 叉积运算符。 有关其用法，请查看 GeoPlane.js 中 Multiple 方法的引用。
// Constructor
function GeoVector(p0, p1) // GeoPoint p0, p1
{
// vector begin point
this.p0 = p0;
// vector end point
this.p1 = p1;
// vector x axis projection value
this.x = p1.x - p0.x;
// vector y axis projection value
this.y = p1.y - p0.y;
// vector z axis projection value
this.z = p1.z - p0.z;
}
// public Instance method to simulate overloading
// binary operator multiple (GeoVector * GeoVector)
GeoVector.prototype.Multiple = function(v) // GeoVector v
{
var x = this.y * v.z - this.z * v.y;
var y = this.z * v.x - this.x * v.z;
var z = this.x * v.y - this.y * v.x;
var p0 = v.p0; // GeoPoint
var p1 = p0.Add(new GeoPoint(x, y, z)); // GeoPoint
return new GeoVector(p0, p1);
};
GeoPlane.js
Negative 方法用于模拟一元负运算符，Multiple 方法用于模拟 GeoPlane 和 GeoPoint 的点积，它们的引用在 GeoPolygonProc.js 中。
// Constructor
function GeoPlane(a, b, c, d) // double a, b, c, d
{
this.a = a;
this.b = b;
this.c = c;
this.d = d;
};
// public Static method to simulate the second constructor
GeoPlane.Create = function(p0, p1, p2) // GeoPoint p0, p1, p2
{
var v = new GeoVector(p0, p1);
var u = new GeoVector(p0, p2);
// normal vector.
var n = u.Multiple(v); // GeoVector
var a = n.x; // double
var b = n.y; // double
var c = n.z; // double
var d = - (a * p0.x + b * p0.y + c * p0.z); // double
return new GeoPlane(a, b, c, d);
}
// public Instance method to simulate overloading unary operator negative - GeoPlane
GeoPlane.prototype.Negative = function()
{
return new GeoPlane( - this.a, - this.b, - this.c, - this.d);
};
// public Instance method to simulate overloading binary operator multiple
// (GeoPlane * GeoPoint)
GeoPlane.prototype.Multiple = function(p) // GeoPoint p
{
return (this.a * p.x + this.b * p.y + this.c * p.z + this.d); // double
};
GeoFace.js
如果使用 [] 创建 Array，则默认 Array 长度为 1，因此将长度设置为 0 是必要的，以帮助其引用使用默认的 0 length 假设。
// Constructor
function GeoFace(p, idx) // List<GeoPoint> p, List<int> idx
{
// new List<GeoPoint>()
this.v = [];
this.v.length = 0;
// new List<GeoPoint>()
this.idx = [];
this.idx.length = 0;
this.n = p.length;
for(var i = 0; i < this.n; i ++ )
{
this.v.push(p[i]);
this.idx.push(idx[i]);
}
}
GeoPolygon.js
// Constructor
function GeoPolygon(p) // List<GeoPoint> p
{
// new List<GeoPoint> ()
this.v = [];
this.v.length = 0;
// new List<int> ()
this.idx = [];
this.idx.length = 0;
this.n = p.length;
for(var i = 0; i < this.n; i ++ )
{
this.v.push(p[i]);
this.idx.push(i);
}
}
Utility.js
ContainsList 方法是如何在 JavaScript 中声明一个 public static 方法的示例。参数列表是一个 2D Array，listItem 是一个 1D Array，该函数的唯一引用在 GeoPolygonProc.js 中，它用于通过检查 2D Array faceVerticeIndex 是否已包含新找到的面顶点 1D Array verticeIndexInOneFace 来避免添加重复面。
// Constructor
function Utility()
{
}
// public Static method
Utility.ContainsList = function(list, listItem) // List<List<T>> list, List<T> listItem
{
if(list == null || listItem == null)
{
return false;
}
listItem.sort();
for (var i = 0; i < list.length; i ++ )
{
var temp = list[i]; // List<T>
if (temp.length == listItem.length)
{
temp.sort();
var itemEqual = true;
for(var j = 0; j < listItem.length; j ++ )
{
if(temp[j] != listItem[j])
{
itemEqual = false;
break;
}
}
if(itemEqual)
{
return true;
}
}
else
{
return false;
}
}
return false;
}
GeoPolygonProc.js
构造函数使用 "call" 方法来引用外部 private 方法。 所有 private 方法都用于设置 GeoPolygonProc public 实例成员。
this.Faces: 所有面的 1DGeoFaceArraythis.FacePlanes: 所有面平面的 1DGeoPlaneArraythis.NumberOfFaces: 面的数量this.x0,this.x1,this.y0,this.y1,this.z0,this.z1: 多边形边界
唯一的 public 方法是 PointInside3DPolygon。 所有 public 实例成员和方法都在 JSLASProc.js 中引用。
faceVerticeIndex: 2D Array，每个面中的索引。面索引是第一主要维度，一个面中的顶点索引是第二次要维度。 原始 3D 多边形顶点索引根据不同的面拆分为不同的组。maxError: 此值是点到面三角形平面距离阈值，用于确定另一个顶点是否在面三角形平面中。
var MaxUnitMeasureError = 0.001;
// Constructor
function GeoPolygonProc(polygonInst) // GeoPolygon polygonInst
{
// Set boundary
Set3DPolygonBoundary.call(this, polygonInst);
// Set maximum point to face plane distance error,
Set3DPolygonUnitError.call(this);
// Set faces and face planes
SetConvex3DFaces.call(this, polygonInst);
}
// private Instance method
function Set3DPolygonBoundary(polygon) // GeoPolygon polygon
{
// List<GeoPoint>
var vertices = polygon.v;
var n = polygon.n;
var xmin, xmax, ymin, ymax, zmin, zmax;
xmin = xmax = vertices[0].x;
ymin = ymax = vertices[0].y;
zmin = zmax = vertices[0].z;
for (var i = 1; i < n; i ++ )
{
if (vertices[i].x < xmin) xmin = vertices[i].x;
if (vertices[i].y < ymin) ymin = vertices[i].y;
if (vertices[i].z < zmin) zmin = vertices[i].z;
if (vertices[i].x > xmax) xmax = vertices[i].x;
if (vertices[i].y > ymax) ymax = vertices[i].y;
if (vertices[i].z > zmax) zmax = vertices[i].z;
}
// double
this.x0 = xmin;
this.x1 = xmax;
this.y0 = ymin;
this.y1 = ymax;
this.z0 = zmin;
this.z1 = zmax;
}
// private Instance method
function Set3DPolygonUnitError()
{
// double
this.MaxDisError = ((Math.abs(this.x0) + Math.abs(this.x1) +
Math.abs(this.y0) + Math.abs(this.y1) +
Math.abs(this.z0) + Math.abs(this.z1)) / 6 * MaxUnitMeasureError);
}
// private Instance method
function SetConvex3DFaces(polygon) // GeoPolygon polygon
{
// new List<GeoFace>()
var faces = [];
faces.length = 0;
// new List<GeoPlane>()
var facePlanes = [];
facePlanes.length = 0;
var numberOfFaces;
var maxError = this.MaxDisError;
// vertices of 3D polygon
var vertices = polygon.v; // new List<GeoPoint>()
var n = polygon.n;
// vertices indexes for all faces
// vertices index is the original index value in the input polygon
// new List<List<int>>()
var faceVerticeIndex = [[]];
faceVerticeIndex.length = 0;
// face planes for all faces
// new List<GeoPlane> ()
var fpOutward = [];
fpOutward.length = 0;
for(var i = 0; i < n; i ++ )
{
// triangle point 1
// GeoPoint
var p0 = vertices[i];
for(var j = i + 1; j < n; j ++ )
{
// triangle point 2
// GeoPoint
var p1 = vertices[j];
for(var k = j + 1; k < n; k ++ )
{
// triangle point 3
// GeoPoint
var p2 = vertices[k];
var trianglePlane = GeoPlane.Create(p0, p1, p2);
var onLeftCount = 0;
var onRightCount = 0;
// indexes of points that lie in same plane with face triangle plane
// new List<int>()
var pointInSamePlaneIndex = [];
pointInSamePlaneIndex.length = 0;
for(var l = 0; l < n ; l ++ )
{
// any point other than the 3 triangle points
if(l != i && l != j && l != k)
{
// GeoPoint
var p = vertices[l];
// double
var dis = trianglePlane.Multiple(p);
// next point is in the triangle plane
if(Math.abs(dis) < maxError)
{
pointInSamePlaneIndex.push(l);
}
else
{
if(dis < 0)
{
onLeftCount ++ ;
}
else
{
onRightCount ++ ;
}
}
}
}
// This is a face for a CONVEX 3d polygon.
// For a CONCAVE 3d polygon, this maybe not a face.
if(onLeftCount == 0 || onRightCount == 0)
{
// new List<int>()
var verticeIndexInOneFace = [];
verticeIndexInOneFace.length = 0;
// triangle plane
verticeIndexInOneFace.push(i);
verticeIndexInOneFace.push(j);
verticeIndexInOneFace.push(k);
var m = pointInSamePlaneIndex.length;
if(m > 0) // there are other vertices in this triangle plane
{
for(var p = 0; p < m; p ++ )
{
verticeIndexInOneFace.push(pointInSamePlaneIndex[p]);
}
}
// if verticeIndexInOneFace is a new face
if ( ! Utility.ContainsList(faceVerticeIndex, verticeIndexInOneFace))
{
// add it in the faceVerticeIndex list
faceVerticeIndex.push(verticeIndexInOneFace);
// add the trianglePlane in the face plane list fpOutward.
if (onRightCount == 0)
{
fpOutward.push(trianglePlane);
}
else if (onLeftCount == 0)
{
fpOutward.push(trianglePlane.Negative());
}
}
}
else
{
// possible reasons :
// 1. the plane is not a face of a convex 3d polygon,
// it is a plane crossing the convex 3d polygon.
// 2. the plane is a face of a concave 3d polygon
}
} // k loop
} // j loop
} // i loop
numberOfFaces = faceVerticeIndex.length;
for (var i = 0; i < numberOfFaces; i ++ )
{
facePlanes.push(new GeoPlane(fpOutward[i].a, fpOutward[i].b,
fpOutward[i].c, fpOutward[i].d));
// new List<GeoPoint>()
var gp = [];
gp.length = 0;
// new List<int>()
var vi = [];
vi.length = 0;
var count = faceVerticeIndex[i].length;
for (var j = 0; j < count; j ++ )
{
vi.push(faceVerticeIndex[i][j]);
gp.push( new GeoPoint(vertices[ vi[j] ].x,
vertices[ vi[j] ].y,
vertices[ vi[j] ].z));
}
faces.push(new GeoFace(gp, vi));
}
// List<GeoFace>
this.Faces = faces;
// List<GeoPlane>
this.FacePlanes = facePlanes;
this.NumberOfFaces = numberOfFaces;
}
// public Instance method
GeoPolygonProc.prototype.PointInside3DPolygon = function(x, y, z) // GeoPoint x, y, z
{
var P = new GeoPoint(x, y, z);
for (var i = 0; i < this.NumberOfFaces; i ++ )
{
var dis = this.FacePlanes[i].Multiple(P);
// If the point is in the same half space with normal vector for any face of the cube,
// then it is outside of the 3D polygon
if (dis > 0)
{
return false;
}
}
// If the point is in the opposite half space with normal vector for all 6 faces,
// then it is inside of the 3D polygon
return true;
}
JSLASProc.js
要读取 LAS 文件，使用 html5 FileReader、DataView 和 ArrayBuffer，"slice" 方法用于模拟文件指针移动的 Seek 函数。
headerBuffer 被延后添加到输出 Array lasInsideBytes 的前面，因为输出 LAS 标头中的记录数需要使用实际写入计数 insideCount 进行更新。"unshift" 方法用于在 Array 前面添加缓冲区。
FileReader 是异步的，readAsArrayBuffer 的结果仅在其 "onloadend" 事件中准备就绪。
目前，Firefox 和 IE 仍然不支持 html5 FileSystem & FileWriter API。 因此，未选择这种最新技术进行文件输出。 另一个选择是 ActiveXObject，同样的原因被丢弃，它对浏览器不友好。
Blob 和 createObjectURL 最终用于创建输出文本文件和二进制 LAS 文件，然后将它们附加到下载链接，通过单击下载链接输出文件。
这里有趣的一件事是 Array 的长度是它的项目计数，如果使用 "push" 方法 10 次，那么 Array 的长度是 10，无论每次推送多少字节。
另一个小技巧是在文本文件输出字符串变量 lasTextStr 中使用 "\r\n" 作为回车符。 单独使用 "\r" 和 "\n" 都不起作用。
var offset_to_point_data_value, record_num, record_len;
var x_scale, y_scale, z_scale, x_offset, y_offset, z_offset;
function GetLASFile()
{
var files = document.getElementById("LASFiles").files;
if ( ! files.length)
{
alert('Please select a LAS file!');
return;
}
var file = files[0];
return file;
}
function AppendInfoInPage(infoStr)
{
var node = document.getElementById('infoID');
var newNode = document.createElement('div');
newNode.appendChild(document.createTextNode(infoStr));
node.appendChild(newNode);
}
function MakeTextFile(str)
{
var blobData = new Blob([str],
{
type : 'text/plain'
}
);
textFile = window.URL.createObjectURL(blobData);
return textFile;
}
function MakeBinaryFile(bytesArray)
{
var blobData = new Blob(bytesArray,
{
type : "octet/stream"
}
);
binFile = window.URL.createObjectURL(blobData);
return binFile;
}
function ActiveTextDownloadLink(downloadID, downloadStr)
{
var downloadLink = document.getElementById(downloadID);
downloadLink.href = MakeTextFile(downloadStr);
downloadLink.style.display = 'block';
}
function ActiveBinaryDownloadLink(downloadID, downloadBytesArray)
{
var downloadLink = document.getElementById(downloadID);
downloadLink.href = MakeBinaryFile(downloadBytesArray);
downloadLink.style.display = 'block';
}
function WriteInsidePointsFile(procInst, lasBuffer)
{
AppendInfoInPage("First 10 Inside Points:");
var insideCount = 0;
var lasTextStr = "";
var lasInsideBytes = [];
lasInsideBytes.length = 0;
for(var i = 0; i < record_num; i ++ )
{
var record_loc = offset_to_point_data_value + record_len * i;
var coorBuffer = lasBuffer.slice(record_loc, record_loc + 12);
var dv = new DataView(coorBuffer);
var littleEndian = true;
// Record coordinates
var pos = 0;
var x = dv.getInt32(pos, littleEndian);
pos += 4;
var y = dv.getInt32(pos, littleEndian);
pos += 4;
var z = dv.getInt32(pos, littleEndian);
// Actual coordinates
var ax = (x * x_scale) + x_offset;
var ay = (y * y_scale) + y_offset;
var az = (z * z_scale) + z_offset;
// Filter out the points outside of boundary to reduce computing
if(ax > procInst.x0 && ax < procInst.x1 &&
ay > procInst.y0 && ay < procInst.y1 &&
az > procInst.z0 && az < procInst.z1)
{
// Main Process to check if the point is inside of the cube
if(procInst.PointInside3DPolygon(ax, ay, az))
{
var coorStr = ax.toFixed(6) + ', ' + ay.toFixed(6) + ', ' + az.toFixed(6);
if(insideCount < 10)
{
AppendInfoInPage(coorStr);
}
// Write the actual coordinates of inside point to text file
lasTextStr += coorStr + "\r\n";
// Get point record
var recordBuffer = lasBuffer.slice(record_loc, record_loc + record_len);
// Write inside point LAS record
lasInsideBytes.push(recordBuffer);
insideCount ++ ;
}
}
}
AppendInfoInPage("Total Inside Points Count: " + insideCount);
AppendInfoInPage("Total Points Count: " + record_num);
// Get LAS header
var headerBuffer = lasBuffer.slice(0, offset_to_point_data_value);
// Update total point number with actual writting count
var dv = new DataView(headerBuffer);
dv.setUint32(107, insideCount, true);
// Write LAS header in front
lasInsideBytes.unshift(headerBuffer);
ActiveTextDownloadLink('insideTextLink', lasTextStr);
ActiveBinaryDownloadLink('insideLASLink', lasInsideBytes);
}
function ProcLASFileHeader(lasBuffer)
{
// LAS Header part 1
var partBuffer = lasBuffer.slice(96, 96 + 15);
// html5 DataView
var dv = new DataView(partBuffer);
var littleEndian = true;
var pos = 0;
offset_to_point_data_value = dv.getUint32(pos, littleEndian);
pos += 4;
var variable_len_num = dv.getUint32(pos, littleEndian);
pos += 4;
var record_type_c = dv.getUint8(pos, littleEndian);
pos += 1;
record_len = dv.getUint16(pos, littleEndian);
pos += 2;
record_num = dv.getUint32(pos, littleEndian);
// LAS Header part 2
partBuffer = lasBuffer.slice(131, 131 + 8 * 6);
dv = new DataView(partBuffer);
pos = 0;
x_scale = dv.getFloat64(pos, littleEndian);
pos += 8;
y_scale = dv.getFloat64(pos, littleEndian);
pos += 8;
z_scale = dv.getFloat64(pos, littleEndian);
pos += 8;
x_offset = dv.getFloat64(pos, littleEndian);
pos += 8;
y_offset = dv.getFloat64(pos, littleEndian);
pos += 8;
z_offset = dv.getFloat64(pos, littleEndian);
pos += 8;
// Verify the result via displaying them in page
AppendInfoInPage("offset_to_point_data_value: " + offset_to_point_data_value +
", record_len: " + record_len +
", record_num: " + record_num);
AppendInfoInPage("x_scale: " + x_scale + ", y_scale: " + y_scale + ", z_scale: " + z_scale +
", x_offset: " + x_offset + ", y_offset: " + y_offset +
", z_offset: " + z_offset);
}
function LASProc()
{
var p1 = new GeoPoint( - 27.28046, 37.11775, - 39.03485);
var p2 = new GeoPoint( - 44.40014, 38.50727, - 28.78860);
var p3 = new GeoPoint( - 49.63065, 20.24757, - 35.05160);
var p4 = new GeoPoint( - 32.51096, 18.85805, - 45.29785);
var p5 = new GeoPoint( - 23.59142, 10.81737, - 29.30445);
var p6 = new GeoPoint( - 18.36091, 29.07707, - 23.04144);
var p7 = new GeoPoint( - 35.48060, 30.46659, - 12.79519);
var p8 = new GeoPoint( - 40.71110, 12.20689, - 19.05819);
var gp = [p1, p2, p3, p4, p5, p6, p7, p8];
var gpInst = new GeoPolygon(gp);
var procInst = new GeoPolygonProc(gpInst);
// Get LAS file object via html file selector
// Javascript does not be allowed to open a file from file system as default
var lasFile = GetLASFile();
// html5 FileReader is Asynchronous
var lasFileReader = new FileReader();
// Asynchronous read, process the result till the read is done
lasFileReader.readAsArrayBuffer(lasFile);
// process the result buffer till Asynchronous readAsArrayBuffer is done
lasFileReader.onloadend = function(evt)
{
var lasBuffer = lasFileReader.result;
ProcLASFileHeader(lasBuffer);
WriteInsidePointsFile(procInst, lasBuffer);
}
}
JSLASProc.htm
通过文件选择器选择 LAS 文件 cube.las 后，单击按钮 "LAS Process"。 这两个下载链接将出现在页面中，并且还会显示一些简短的 LAS 信息。
<html>
<!-- include all JavaScript in GeoProc JS library
with the same order of GeoProc\GeoProc.htm -->
<script src="..\GeoProc\Utility.js"> </script>
<script src="..\GeoProc\GeoPoint.js"> </script>
<script src="..\GeoProc\GeoVector.js"> </script>
<script src="..\GeoProc\GeoPlane.js"> </script>
<script src="..\GeoProc\GeoPolygon.js"> </script>
<script src="..\GeoProc\GeoFace.js"> </script>
<script src="..\GeoProc\GeoPolygonProc.js"> </script>
<head></head>
<body>
<!-- include LAS Proc JavaScript -->
<script src="JSLASProc.js"> </script>
<div>
<form>
<input type="file" id="LASFiles" name="lasfile" />
<input type="button" id="btnLASProc" value="LAS Process" onclick="LASProc();"/>
<a download="cube_inside.txt" id="insideTextLink" style="display: none">
Download Inside Text</a>
<a download="cube_inside.las" id="insideLASLink" style="display: none">
Download Inside LAS</a>
</form>
</div>
<div id=infoID></div>
</body>
</html>
关注点
JavaScript 版本在 C#、Java 和 C++ 版本中是最快的。
历史
- 2016 年 1 月 15 日：初始版本