Point Inside 3D Convex Polygon in JavaScript
5.00/5 (2 votes)
An algorithm to determine if a point is inside a 3D convex polygon for a given polygon vertices in JavaScript
Introduction
This is a JavaScript version of the original C++ algorithm which can be found here.
Background
Please refer to the original C++ algorithm here.
Using the Code
The algorithm is wrapped into a JavaScript library folder GeoProc. The main test program JSLASProc reads point cloud data from a LAS file, then writes all inside points into a new LAS file. A LAS file is able to be displayed with a freeware FugroViewer.
Each JavaScript in GeoProc library folder has an usage and test part. To test each function, uncomment the test part and open GeoProc.htm in browser.
GeoPoint.js
Javascript has no operator overloading, the Add method is to simulate operator overloading. Check the reference of Add method in GeoVector.js for its usage.
// 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
The method Multiple is to simulate GeoVector cross product operator. Check the reference of Multiple method in GeoPlane.js for its usage.
// 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
The method Negative is to simulate unary negative operator, the method Multiple is to simulate dot product of GeoPlane and GeoPoint, their references are in 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
The default Array length is 1 if using [] to create Array, so set length to 0 is neccessary to help their references to use default 0 length assumption.
// 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
The method ContainsList is an example for how to declare a public static method in JavaScript. The argument list is a 2D Array, the listItem is a 1D Array, the only reference of the function is in GeoPolygonProc.js, it is used to avoid adding duplicate face via checking if the 2D Array faceVerticeIndex already contains the newly found face vertices 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
The constructor uses "call" method to reference outside private method. All private methods are used to set GeoPolygonProc public instance members:
this.Faces: 1DGeoFaceArray for all facesthis.FacePlanes: 1DGeoPlaneArray for all faces planesthis.NumberOfFaces: Number of facesthis.x0,this.x1,this.y0,this.y1,this.z0,this.z1: Polygon boundary
The only public method is PointInside3DPolygon. All public instance member and methods are referenced in JSLASProc.js.
faceVerticeIndex: 2D Array, indexes in each face for all faces. Face index is first major dimension, vertice indexes in one face is second minor dimension. The original 3D polygon vertices indexes are split to difference groups according to different faces.maxError: This value is a point to face triangle plane distance threshold to determine if the other vertice is in the face triangle plane.
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
To read LAS file, html5 FileReader, DataView and ArrayBuffer are used, "slice" method is used to simulate Seek function for file pointer movement.
The headerBuffer is postponded to be added in front of the output Array lasInsideBytes since record number in the output LAS header needs to update with actual writting count insideCount. The "unshift" method is used to add buffer in front of Array.
FileReader is Asynchronous, the result of readAsArrayBuffer is only ready inside its "onloadend" event.
At the time, html5 FileSystem & FileWriter API are still not supported by Firefox and IE. So this latest technique is not chosen for file output. Another choice is ActiveXObject, same reason it is discarded, it is not browser friendly.
Blob and createObjectURL are finally used to create the outout text file and binary LAS file, then attach them to download links, output the files via clicking the download links.
An interesting thing here is an Array length is its items count, if "push" methods are used 10 times, then the Array length is 10, no matter how many bytes are pushed in for each time.
Another small trick is using "\r\n" as carrage return in text file output string variable lasTextStr. Both "\r" and "\n" do not work.
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
After select the LAS file cube.las via the file selector, click button "LAS Process". The two download links will appear in the page and some brief LAS info will show up as well.
<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>
Points of Interest
JavaScript version is fastest in the C#, Java and C++ versions.
History
- 15th January, 2016: Initial version