Point Inside 3D Convex Polygon in JavaScript

• Download JSLASProc - 5.4 MB
• Download FugroViewerSetup22 - 2.9 MB

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: 1D GeoFace Array for all faces
• this.FacePlanes: 1D GeoPlane Array for all faces planes
• this.NumberOfFaces: Number of faces
• this.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

• 15^th January, 2016: Initial version