Point Inside 3D Convex Polygon in PHP

• Download PHPLASProc - 5.5 MB
• Download FugroViewerSetup22 - 2.9 MB

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:

1. Declare a 2d array $faceVerticeIndex, first dimension is face index, second dimension is face vertice index.
2. Go through all given vertices, construct a triangle plane with any 3 vertices combination.
3. Determine if the triangle plane is a face through checking if all other vertices in same half space.
4. Declare a 1d array $verticeIndexInOneFace for a complete vertices indexes in one face.
5. Find other vertices in the same plane with the triangle plane, put them in 1d array $pointInSamePlaneIndex.
6. Add unique face to $faceVerticeIndex with adding the 3 triangle vertices and the same plane vertices.
7. 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

• 17^th January, 2016: Initial version