Point Inside 3D Convex Polygon in PHP
5.00/5 (1 vote)
An algorithm to determine if a point is inside a 3D convex polygon for a given polygon vertices in PHP
Introduction
This is a PHP 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 PHP class library folder GeoProc. The main test program PHPLASProc 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.
To consume the GeoProc library, prepare a GeoPolygonProc instance first like this:
$gp = [$p1, $p2, $p3, $p4, $p5, $p6, $p7, $p8];
$gpInst = new GeoPolygon($gp);
$procInst = new GeoPolygonProc($gpInst);
In the above code, $gp is a GeoPoint array; $gpInst is a GeoPolygon instance created from a GeoPoint array as its vertices; $procInst is a GeoPolygonProc instance created from a GeoPolygon instance.
In the main test program PHPLASProc.php, the $pronInst is used to check the 3D polygon boundary first to pre-filter the inside points for coarse screening, then use its main public method PointInside3DPolygon to get the actual inside points as below:
if($x > $procInst->getX0() && $x < $procInst->getX1() &&
$y > $procInst->getY0() && $x < $procInst->getY1() &&
$z > $procInst->getZ0() && $x < $procInst->getZ1() )
{
// Main Process to check if the point is inside of the cube
if($procInst->PointInside3DPolygon($x, $y, $z))
Here are all the GeoProc classes:
GeoPoint.php
PHP has no operator overloading, a static function Add is used to simulate Add operator for two GeoPoint. It is referenced in the Multiple function in GeoVector.php.
The three coordinates member $x, $y, $z are declared as public to simply its referencing.
<?php
class GeoPoint
{
public $x;
public $y;
public $z;
public function __construct($x, $y, $z) // double x, y, z
{
$this->x = $x;
$this->y = $y;
$this->z = $z;
}
// public static function to simulate
// binary operator add overloading: (GeoPoint + GeoPoint)
public static function Add(GeoPoint $p0, GeoPoint $p1)
{
return new GeoPoint($p0->x + $p1->x, $p0->y + $p1->y, $p0->z + $p1->z);
}
}
?>
GeoVector.php
The static function Multiple is to simulate multiple operator for two GeoVector object. It is referenced in the Create function in GeoPlane.php.
The GeoVector member $p0, $p1, $x, $y, $z are declared as private and use five public get methods to simulate its getter.
<?php
require_once ('GeoPoint.php');
class GeoVector
{
private $p0; // vector begin point
private $p1; // vector end point
private $x; // vector x axis projection value
private $y; // vector y axis projection value
private $z; // vector z axis projection value
public function getP0() {return $this->p0;}
public function getP1() {return $this->p1;}
public function getX() {return $this->x;}
public function getY() {return $this->y;}
public function getZ() {return $this->z;}
public function __construct(GeoPoint $p0, GeoPoint $p1)
{
$this->p0 = $p0;
$this->p1 = $p1;
$this->x = $p1->x - $p0->x;
$this->y = $p1->y - $p0->y;
$this->z = $p1->z - $p0->z;
}
// public Instance method to simulate
// binary operator multiple overloading: (GeoVector * GeoVector)
public static function Multiple(GeoVector $u, GeoVector $v)
{
$x = $u->y * $v->z - $u->z * $v->y;
$y = $u->z * $v->getX() - $u->x * $v->getZ();
$z = $u->x * $v->getY() - $u->y * $v->getX();
$p0 = $v->getP0(); // GeoPoint
$p1 = GeoPoint::Add($p0, new GeoPoint($x, $y, $z)); // GeoPoint
return new GeoVector($p0, $p1);
}
}
?>
GeoPlane.php
This is a generic declaration for a geometry plane, it is used to describe a face plane where a face of the 3D polygon lies in. It is referenced in a GeoPlane array $FacePlanes in GeoPolygonProc.php.
<?php
require_once ('GeoPoint.php');
require_once ('GeoVector.php');
class GeoPlane
{
// Plane Equation: a * x + b * y + c * z + d = 0
private $a;
private $b;
private $c;
private $d;
// public instance function to simulate class property getter
public function getA() {return $this->a;}
public function getB() {return $this->b;}
public function getC() {return $this->c;}
public function getD() {return $this->d;}
// Constructor
public function __construct($a, $b, $c, $d) // double a, b, c, d
{
$this->a = $a;
$this->b = $b;
$this->c = $c;
$this->d = $d;
}
// public static function to simulate constructor overloading
public static function Create
(GeoPoint $p0, GeoPoint $p1, GeoPoint $p2) // GeoPoint p0, p1, p2
{
$v = new GeoVector($p0, $p1);
$u = new GeoVector($p0, $p2);
// normal vector.
$n = GeoVector::Multiple($u, $v); // GeoVector
$a = $n->getX(); // double
$b = $n->getY(); // double
$c = $n->getZ(); // double
$d = - ($a * $p0->x + $b * $p0->y + $c * $p0->z); // double
return new GeoPlane($a, $b, $c, $d);
}
// Simulate uary operator negative overloading: -GeoPlane
public static function Negative(GeoPlane $pl)
{
return new GeoPlane( - $pl->getA(), - $pl->getB(), - $pl->getC(), - $pl->getD());
}
// Simulate binary operator multiple overloading binary operator multiple:
// GeoPlane * GeoPoint
public static function Multiple(GeoPlane $pl, GeoPoint $pt)
{
return ($pl->getA() * $pt->x + $pl->getB() *
$pt->y + $pl->getC() * $pt->z + $pl->getD()); // double
}
}
?>
GeoPolygon.php
The class declared the 3D polgyon with its all vertices. It is used to construct GeoPolygonProc instance with initializing its boundary, face planes and faces. Its references are in GeoPolygonProc.
<?php
class GeoPolygon
{
private $v; // 3D polygon vertices coordinates: GeoPoint Array
private $idx; // 3D polygon vertices indexes: Integer Array
private $n; // number of 3D polygon vertices: Integer
public function getV(){return $this->v;}
public function getIdx(){return $this->idx;}
public function getN(){return $this->n;}
public function __construct($v) // $v: polygon vertices coordinates: GeoPoint Array
{
// alloc instance array memory
$this->v = [];
$this->idx = [];
$this->n = count($v);
for($i=0; $i< $this->n; $i++)
{
array_push($this->v, $v[$i]);
array_push($this->idx, $i);
}
}
public function __destruct()
{
// free instance array memory
unset($this->v);
unset($this->idx);
}
}
?>
GeoFace.php
The class declares a 3D polygon face with the face vertices. The face vertices indexes are from a GeoPolygon vertice indexes which the face belongs to.
<?php
class GeoFace
{
private $v; // 3D polygon face vertices coordinates: GeoPoint Array
private $idx; // 3D polygon face vertices indexes: Integer Array
private $n; // number of 3D polygon face vertices: Integer
public function getV(){return $this->v;}
public function getIdx(){return $this->idx;}
public function getN(){return $this->n;}
// $v: 3D polygon face vertices coordinates: GeoPoint Array
// $idx: 3D polygon face vertices index
public function __construct(array $v, array $idx)
{
// alloc instance array memory
$this->v = [];
$this->idx = [];
$this->n = count($v);
for($i=0; $i< $this->n; $i++)
{
array_push($this->v, $v[$i]);
array_push($this->idx, $idx[$i]);
}
}
public function __destruct()
{
// free instance array memory
unset($this->v);
unset($this->idx);
}
}
?>
Utility.php
The class has only one function to check if a 2d array contains a 1d array. The function is used to prevent the 2d face index array $faceVerticeIndex from adding the duplicate new face $verticeIndexInOneFace in GeoPolygonProc.php.
<?php
class Utility
{
// public Static method
public static function ContainsList
($listList, $listItem) // listList: [[]] list, listItem: []
{
if($listList == null || $listItem == null)
{
return false;
}
sort($listItem);
for ($i = 0; $i < count($listList); $i ++ )
{
$temp = $listList[$i];
if (count($temp) == count($listItem))
{
sort($temp);
$itemEqual = true;
for($j = 0; $j < count($listItem); $j ++ )
{
if($temp[$j] != $listItem[$j])
{
$itemEqual = false;
break;
}
}
if($itemEqual)
{
return true;
}
}
else
{
return false;
}
}
return false;
}
}
?>
GeoPolygonProc.php
This is the main class for referencing the GeoProc library. It is capable of extending its functions to do more complex 3D polygon geometry processing based on the GeoProc library. The current GeoPolygonProc provides basic settings for 3D polygon boundary, face planes and faces. The particular public function PointInside3DPolygon is a role model for adding other functions to solve 3D polygon problems.
SetConvex3DFaces is the essential method, it extracts all faces information from the simple vertices. The procedure is:
- Declare a 2d array
$faceVerticeIndex, first dimension is face index, second dimension is face vertice index. - Go through all given vertices, construct a triangle plane with any 3 vertices combination.
- Determine if the triangle plane is a face through checking if all other vertices in same half space.
- Declare a 1d array
$verticeIndexInOneFacefor a complete vertices indexes in one face. - Find other vertices in the same plane with the triangle plane, put them in 1d array
$pointInSamePlaneIndex. - Add unique face to
$faceVerticeIndexwith adding the 3 triangle vertices and the same plane vertices. - Get all face planes and faces.
<?php
require_once ('GeoPoint.php');
require_once ('GeoVector.php');
require_once ('GeoPlane.php');
require_once ('GeoPolygon.php');
require_once ('GeoFace.php');
require_once ('Utility.php');
class GeoPolygonProc
{
private $MaxUnitMeasureError = 0.001;
private $Faces; // GeoFace Array
private $FacePlanes; // GeoPlane Array
private $NumberOfFaces; // Integer
private $x0, $y0, $x1, $y1, $z0, $z1; // 3D polygon boundary: Double
public $MaxDisError; // Allow to set and get
public function getFaces(){return $this->Faces;}
public function getFacePlanes(){return $this->FacePlanes;}
public function getNumberOfFaces(){return $this->NumberOfFaces;}
public function getX0(){return $this->x0;}
public function getX1(){return $this->x1;}
public function getY0(){return $this->y0;}
public function getY1(){return $this->y1;}
public function getZ0(){return $this->z0;}
public function getZ1(){return $this->z1;}
// Constructor
public function __construct($polygonInst) // $polygonInst: GeoPolygon
{
// Set boundary
$this->Set3DPolygonBoundary($polygonInst);
// Set maximum point to face plane distance error,
$this->Set3DPolygonUnitError();
// Set faces and face planes
$this->SetConvex3DFaces($polygonInst);
}
public function __destruct()
{
// free instance array memory
unset($this->Faces);
unset($this->FaceFacePlaness);
}
// private Instance method
private function Set3DPolygonBoundary($polygon) // $polygon: GeoPolygon polygon
{
$vertices = $polygon->getV(); // GeoPoint Array
$n = $polygon->getN();
$xmin = $vertices[0]->x;
$xmax = $vertices[0]->x;
$ymin = $vertices[0]->y;
$ymax = $vertices[0]->y;
$zmin = $vertices[0]->z;
$zmax = $vertices[0]->z;
for ($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
private function Set3DPolygonUnitError()
{
$avarageError = (abs($this->x0) + abs($this->x1) +
abs($this->y0) + abs($this->y1) +
abs($this->z0) + abs($this->z1)) / 6;
$this->MaxDisError = $avarageError * $this->MaxUnitMeasureError;
}
// private Instance method
private function SetConvex3DFaces($polygon) // $polygon: GeoPolygon
{
// alloc instance array memory
$this->Faces = []; // GeoFace Array
$this->FacePlanes = []; // GeoPlane Array
// local 2d integer array memory, only declare as 1d array
// with face index as major array index
// push minor face indexes 1d array as its element later
// vertices indexes for all faces
// vertices index is the original index value in the input polygon
$faceVerticeIndex = [];
// local GeoPlane array memory
// face planes for all faces
$fpOutward = [];
// vertices of polygon
$vertices = $polygon->getV(); // GeoPoint Array
// number of polygon vertices
$n = $polygon->getN();
for($i = 0; $i < $n; $i ++ )
{
// triangle point 1
$p0 = $vertices[$i];
for($j = $i + 1; $j < $n; $j ++ )
{
// triangle point 2
$p1 = $vertices[$j];
for($k = $j + 1; $k < $n; $k ++ )
{
// triangle point 3
$p2 = $vertices[$k];
$trianglePlane = GeoPlane::Create($p0, $p1, $p2);
$onLeftCount = 0;
$onRightCount = 0;
// alloc local Integer array memory
// indexes of points that lie in same plane with face triangle plane
// new List<int>()
$pointInSamePlaneIndex = [];
for($l = 0; $l < $n ; $l ++ )
{
// any point other than the 3 triangle points
if($l != $i && $l != $j && $l != $k)
{
// GeoPoint
$pt = new GeoPoint($vertices[$l]->x,
$vertices[$l]->y, $vertices[$l]->z);
// double
$dis = GeoPlane::Multiple($trianglePlane, $pt);
// next point is in the triangle plane
if(abs($dis) < $this->MaxDisError)
{
array_push($pointInSamePlaneIndex, $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)
{
// alloc local Integer array memory
$verticeIndexInOneFace = [];
// triangle plane
array_push($verticeIndexInOneFace, $i);
array_push($verticeIndexInOneFace, $j);
array_push($verticeIndexInOneFace, $k);
$m = count($pointInSamePlaneIndex);
if($m > 0) // there are other vertices in this triangle plane
{
for($p = 0; $p < $m; $p ++ )
{
array_push($verticeIndexInOneFace, $pointInSamePlaneIndex[$p]);
}
}
// if verticeIndexInOneFace is a new face
if ( ! Utility::ContainsList($faceVerticeIndex, $verticeIndexInOneFace))
{
// add it in the faceVerticeIndex list
array_push($faceVerticeIndex, $verticeIndexInOneFace);
// add the trianglePlane in the face plane list fpOutward.
if ($onRightCount == 0)
{
array_push($fpOutward, $trianglePlane);
}
else if ($onLeftCount == 0)
{
array_push($fpOutward, GeoPlane::Negative($trianglePlane));
}
}
}
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
$this->NumberOfFaces = count($faceVerticeIndex);
for ($i = 0; $i < $this->NumberOfFaces; $i ++ )
{
array_push($this->FacePlanes, new GeoPlane(
$fpOutward[$i]->getA(), $fpOutward[$i]->getB(),
$fpOutward[$i]->getC(), $fpOutward[$i]->getD() ));
// alloc local GeoPoint array memory
$gp = [];
// alloc local Integer array memory
$vi = [];
$count = count($faceVerticeIndex[$i]);
for ($j = 0; $j < $count; $j ++ )
{
array_push($vi, $faceVerticeIndex[$i][$j]);
array_push($gp, new GeoPoint($vertices[ $vi[$j] ]->x,
$vertices[ $vi[$j] ]->y,
$vertices[ $vi[$j] ]->z));
}
array_push($this->Faces, new GeoFace($gp, $vi));
}
// free local array memory
unset($faceVerticeIndex);
unset($fpOutward);
unset($pointInSamePlaneIndex);
unset($verticeIndexInOneFace);
unset($gp);
unset($vi);
}
// public Instance method
public function PointInside3DPolygon($x, $y, $z) // $x, $y, $z: GeoPoint
{
$pt = new GeoPoint($x, $y, $z);
for ($i = 0; $i < $this->NumberOfFaces; $i ++ )
{
$dis = GeoPlane::Multiple($this->FacePlanes[$i], $pt);
// 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;
}
}
?>
PHPLASProc.php
This is main test program to consume the GeoProc library, it reads original LAS file cube.las, then gets all inside points against the given 8 3D polygon vertices while writing them to a new LAS file and a coordinate plain text file.
<?php
// include files
require_once ('../GeoProc/GeoPoint.php');
require_once ('../GeoProc/GeoVector.php');
require_once ('../GeoProc/GeoPlane.php');
require_once ('../GeoProc/GeoPolygon.php');
require_once ('../GeoProc/GeoFace.php');
require_once ('../GeoProc/GeoPolygonProc.php');
// prepare GeoPolygonProc instance $procInst
$p1 = new GeoPoint( - 27.28046, 37.11775, - 39.03485);
$p2 = new GeoPoint( - 44.40014, 38.50727, - 28.78860);
$p3 = new GeoPoint( - 49.63065, 20.24757, - 35.05160);
$p4 = new GeoPoint( - 32.51096, 18.85805, - 45.29785);
$p5 = new GeoPoint( - 23.59142, 10.81737, - 29.30445);
$p6 = new GeoPoint( - 18.36091, 29.07707, - 23.04144);
$p7 = new GeoPoint( - 35.48060, 30.46659, - 12.79519);
$p8 = new GeoPoint( - 40.71110, 12.20689, - 19.05819);
$gp = [$p1, $p2, $p3, $p4, $p5, $p6, $p7, $p8];
$gpInst = new GeoPolygon($gp);
$procInst = new GeoPolygonProc($gpInst);
// process LAS file with $procInst
// open LAS file to read
$lasFile = "../LASInput/cube.las";
$rbFile = fopen($lasFile, "rb");
// open LAS file to write
$wbFile = fopen("../LASOutput/cube_inside.las", "wb");
// open Text file to write
$wtFile = fopen("../LASOutput/cube_inside.txt", "w");
// read LAS header parameters in Little-Endian
fseek($rbFile, 96);
$offset_to_point_data_value = unpack('V', fread($rbFile, 4))[1];
$variable_len_num = unpack('V', fread($rbFile, 4))[1];
$record_type = unpack('C', fread($rbFile, 1))[1];
$record_len = unpack('v', fread($rbFile, 2))[1];
$record_num = unpack('V', fread($rbFile, 4))[1];
fseek($rbFile, 131);
$x_scale = unpack('d', fread($rbFile, 8))[1];
$y_scale = unpack('d', fread($rbFile, 8))[1];
$z_scale = unpack('d', fread($rbFile, 8))[1];
$x_offset = unpack('d', fread($rbFile, 8))[1];
$y_offset = unpack('d', fread($rbFile, 8))[1];
$z_offset = unpack('d', fread($rbFile, 8))[1];
// read LAS header
fseek($rbFile, 0);
$headerBuffer = fread($rbFile, $offset_to_point_data_value);
// write LAS public header to LAS
fwrite($wbFile, $headerBuffer);
$insideCount = 0;
// write Point record
for($i = 0; $i < $record_num; $i++)
{
$record_loc = $offset_to_point_data_value + $record_len * $i;
fseek($rbFile, $record_loc);
$xi = unpack('l', fread($rbFile, 4))[1];
$yi = unpack('l', fread($rbFile, 4))[1];
$zi = unpack('l', fread($rbFile, 4))[1];
$x = $xi * $x_scale + $x_offset;
$y = $yi * $y_scale + $y_offset;
$z = $zi * $z_scale + $z_offset;
if($x > $procInst->getX0() && $x < $procInst->getX1() &&
$y > $procInst->getY0() && $x < $procInst->getY1() &&
$z > $procInst->getZ0() && $x < $procInst->getZ1() )
{
// Main Process to check if the point is inside of the cube
if($procInst->PointInside3DPolygon($x, $y, $z))
{
// Write the actual coordinates of inside point to text file
fprintf($wtFile, "%15.6f %15.6f %15.6f \n", $x, $y, $z);
fseek($rbFile, $record_loc);
// Read LAS point record
$recordBuffer = fread($rbFile, $record_len);
// Write LAS point record
fwrite($wbFile, $recordBuffer);
$insideCount++;
}
}
}
// Update total record numbers with actual writing record number
// in output LAS file header
fseek($wbFile, 107);
fwrite($wbFile, pack('V', $insideCount), 4);
// close files
fclose($rbFile);
fclose($wbFile);
fclose($wtFile);
?>
Points of Interest
The PHP binary file process provides fopen, fread, fseek, fwrite and fclose methods, this looks like very similar with C style, but its distinctive "pack" and "unpack" methods must be used, and the Big Endian or Little Endian bytes order must be considered, these features are also used in Python.
History
- 17th January, 2016: Initial version