PHP 中的点在 3D 凸多边形内
5.00/5 (1投票)
一个算法，用于确定一个点是否在给定多边形顶点的 3D 凸多边形内，使用 PHP 实现。
引言
这是原始 C++ 算法的 PHP 版本，可以在 这里 找到。
背景
请参考原始 C++ 算法 这里。
Using the Code
该算法被封装到 PHP 类库文件夹 GeoProc 中。主测试程序 PHPLASProc 从 LAS 文件中读取点云数据，然后将所有内部点写入一个新的 LAS 文件。可以使用免费软件 FugroViewer 显示 LAS 文件。
要使用 GeoProc 库，首先准备一个 GeoPolygonProc 实例，如下所示
$gp = [$p1, $p2, $p3, $p4, $p5, $p6, $p7, $p8];
$gpInst = new GeoPolygon($gp);
$procInst = new GeoPolygonProc($gpInst);
在上面的代码中，$gp 是一个 GeoPoint 数组；$gpInst 是一个从 GeoPoint 数组创建的 GeoPolygon 实例，作为其顶点；$procInst 是一个从 GeoPolygon 实例创建的 GeoPolygonProc 实例。
在主测试程序 PHPLASProc.php 中，首先使用 $pronInst 检查 3D 多边形的边界以进行粗略筛选，然后使用其主要的 public 方法 PointInside3DPolygon 获取实际的内部点，如下所示
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))
以下是所有 GeoProc 类
GeoPoint.php
PHP 没有运算符重载，使用一个 static 函数 Add 来模拟两个 GeoPoint 的 Add 运算符。它在 GeoVector.php 的 Multiple 函数中被引用。
三个坐标成员 $x、$y、$z 声明为 public 以简化其引用。
<?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
static 函数 Multiple 用于模拟两个 GeoVector 对象的乘法运算符。它在 GeoPlane.php 的 Create 函数中被引用。
GeoVector 成员 $p0、$p1、$x、$y、$z 声明为 private，并使用五个 public get 方法来模拟其 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
这是几何平面的通用声明，用于描述 3D 多边形的一个面所在的平面。它在 GeoPolygonProc.php 中的 GeoPlane 数组 $FacePlanes 中被引用。
<?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
该类使用其所有顶点声明了 3D 多边形。它用于构造 GeoPolygonProc 实例，并初始化其边界、面平面和面。它的引用在 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
该类声明了一个具有面顶点的 3D 多边形面。面顶点索引来自 GeoPolygon 的顶点索引，该面属于该顶点索引。
<?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
该类只有一个函数来检查 2d 数组是否包含 1d 数组。该函数用于防止 2d 面索引数组 $faceVerticeIndex 从 GeoPolygonProc.php 中添加重复的新面 $verticeIndexInOneFace。
<?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
这是引用 GeoProc 库的主类。它能够扩展其功能，以便基于 GeoProc 库进行更复杂的 3D 多边形几何处理。当前的 GeoPolygonProc 提供了 3D 多边形边界、面平面和面的基本设置。特定的 public 函数 PointInside3DPolygon 是添加其他函数以解决 3D 多边形问题的模型。
SetConvex3DFaces 是基本方法，它从简单顶点提取所有面信息。该过程是
- 声明一个 2d 数组
$faceVerticeIndex，第一个维度是面索引，第二个维度是面顶点索引。 - 遍历所有给定的顶点，使用任意 3 个顶点组合构造一个三角形平面。
- 通过检查所有其他顶点是否在同一半空间中来确定三角形平面是否为一个面。
- 声明一个 1d 数组
$verticeIndexInOneFace，用于一个面中的完整顶点索引。 - 找到与三角形平面在同一平面中的其他顶点，将它们放入 1d 数组
$pointInSamePlaneIndex中。 - 通过添加 3 个三角形顶点和同一平面的顶点，将唯一的面添加到
$faceVerticeIndex。 - 获取所有面平面和面。
<?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
这是使用 GeoProc 库的主测试程序，它读取原始 LAS 文件 cube.las，然后获取所有内部点，针对给定的 8 个 3D 多边形顶点，同时将它们写入一个新的 LAS 文件和一个坐标纯文本文件。
<?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);
?>
关注点
PHP 二进制文件处理提供了 fopen、fread、fseek、fwrite 和 fclose 方法，这看起来与 C 风格非常相似，但是必须使用其独特的 "pack" 和 "unpack" 方法，并且必须考虑大端或小端字节顺序，这些特性也在 Python 中使用。
历史
- 2016 年 1 月 17 日：初始版本