C++ 中点是否在 3D 凸多边形内
4.90/5 (28投票s)
一个算法，用于确定点是否在给定多边形顶点 C++ 的 3D 凸多边形内。
引言
本文实现了一种算法，该算法利用平面法向量和点到平面距离向量的方向来确定一个点是否在一个给定多边形顶点的 3D 凸多边形内。
这是原始的 C++ 版本，我已经将该算法移植到 C# 版本、Java 版本、JavaScript 版本、PHP 版本、Python 版本、Perl 版本 和 Fortran。
这些版本涵盖了不同的编程类型，从编译语言到解释语言，所有都支持面向对象编程，这使得这一系列移植具有清晰、稳定和一致的通用程序架构。
现在提供了一个 MASM 汇编 版本，尽管汇编语言没有面向对象特性，但在其 STRUCT 和 invoke 机制的帮助下，我们可以使用基于库的架构使程序与所有上述高级语言具有相同的结构。这应该是本系列移植的最后一个版本。
背景
我有幸在激光扫描 3D 领域中一窥点云处理。该领域的一个问题是如何用特定的几何边界来分离点。例如，给定 8 个立方体顶点，如何获取立方体内的所有点。
可以在网上找到几种解决该问题的方法。例如，光线追踪，通过检查交点数量是否为奇数来确定点是否在内部。在某些特殊情况下，也使用笛卡尔坐标变换，例如，如果边界是一个没有任何变形的完美立方体。
我最终以一种简单的方式解决了这个问题。这个想法主要受到讨论点平面距离的参考文章的启发：https://mathworld.net.cn/Point-PlaneDistance.html。
这是问题和我的解决方案。
提问
对于一个给定 N 个顶点的 3D 凸多边形，确定一个 3D 点 (x, y, z) 是否在多边形内。
解决方案
一个 3D 凸多边形有许多面，一个面有一个面平面，该面位于该平面上。
一个面平面有一个向外法向量，它指向多边形的外部。
点到面平面的距离定义了一个几何向量，如果距离向量与向外法向量方向相反，则该点位于面平面的“内部半空间”中，否则，它位于面平面的“外部半空间”中。
如果一个点位于 3D 凸多边形所有面的“内部半空间”中，则该点被确定为在 3D 多边形内。
这就是该算法的基本思想。
算法几何图
在使用上述方法之前，还需要解决另一个问题，我们需要获取多边形所有面的面平面向外法向量。以下是实现此目的的步骤
- 从 N 个顶点中获取任意 3 个顶点的组合，用于构造一个三角形平面。
- 通过 3 个顶点
P0、P1和P2三角形中 2 个边向量u和v的叉积，得到法向量n = (A, B, C)。顶点P0(x0, y0, z0)是 2 个边向量的共同起点。 - 得到方程为
A * x + B * y + C * z + D = 0的面平面。其中D = -(A * x0 + B * y0 + c * z0)。 - 通过检查所有其他顶点是否都在三角形平面的同一半空间中来确定三角形平面是否是面平面。此步骤需要凸面假设。对于凹 3D 多边形，此规则无法区分多边形的真实面和多边形的内部三角形。
- 通过将与这 3 个顶点位于同一三角形平面上的任何其他顶点附加到这 3 个顶点来获取面顶点。
- 总而言之，我们找到了所有面及其向外法向量和完整的顶点。
这是基本公式
对于一个点 pt(X, Y, Z) 和一个平面 pl(A, B, C, D)，值 (pt.X * pl.A + pt.Y * pl.B + pt.Z * pl.C + pl.D) 的符号能够确定点所在的半空间。
Using the Code
该算法被封装到 VC++ DLL GeoProc 中。主测试程序 LASProc 从 LAS 文件读取点云数据，然后将所有内部点写入新的 LAS 文件。LAS 文件可以使用免费软件 FugroViewer 显示。
要使用算法 DLL，首先像这样构造一个多边形对象
// Create polygon object
GeoPolygon polygonObj = GeoPolygon(verticesVector);
然后，构造主处理对象
// Create main process object
GeoPolygonProc procObj = GeoPolygonProc(polygonObj);
检查点 (x, y, z) 是否在 CubeVertices 内的主程序
// Main procedure to check if a point is inside polygon
procObj.PointInside3DPolygon(x, y, z)
以下是 GeoProc DLL 库中的类
代码带有注释，几乎不言自明。
GeoPoint.h
#pragma once
namespace GeoProc
{
// A 3D Geometry Point
class GeoPoint
{
public:
double x;
double y;
double z;
__declspec(dllexport) GeoPoint(void);
__declspec(dllexport) ~GeoPoint(void);
__declspec(dllexport) GeoPoint(double x, double y, double z)
{
this->x = x;
this->y = y;
this->z = z;
}
__declspec(dllexport) friend GeoPoint operator+(const GeoPoint& p0, const GeoPoint& p1)
{
return GeoPoint(p0.x + p1.x, p0.y + p1.y, p0.z + p1.z);
}
};
}
GeoVector.h
#pragma once
#include "GeoPoint.h"
namespace GeoProc
{
// A 3D Geometry Vector
class GeoVector
{
GeoPoint _p0; // vector begin point
GeoPoint _p1; // vector end point
double _x; // vector x axis projection value
double _y; // vector y axis projection value
double _z; // vector z axis projection value
public:
__declspec(dllexport) GeoVector(void);
__declspec(dllexport) ~GeoVector(void);
__declspec(dllexport) GeoVector(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;
}
__declspec(dllexport) friend GeoVector operator*(const GeoVector& u, const GeoVector& v)
{
double x = u._y * v._z - u._z * v._y;
double y = u._z * v._x - u._x * v._z;
double z = u._x * v._y - u._y * v._x;
GeoPoint p0 = v._p0;
GeoPoint p1 = p0 + GeoPoint(x, y, z);
return GeoVector(p0, p1);
}
__declspec(dllexport) GeoPoint GetP0(){return this->_p0;}
__declspec(dllexport) GeoPoint GetP1(){return this->_p1;}
__declspec(dllexport) double GetX(){return this->_x;}
__declspec(dllexport) double GetY(){return this->_y;}
__declspec(dllexport) double GetZ(){return this->_z;}
};
}
GeoPlane.h
#pragma once
#include "GeoVector.h"
namespace GeoProc
{
// A 3D Geometry Plane
class GeoPlane
{
public:
// Plane Equation: A * x + B * y + C * z + D = 0
double a;
double b;
double c;
double d;
__declspec(dllexport) GeoPlane(void);
__declspec(dllexport) ~GeoPlane(void);
__declspec(dllexport) GeoPlane(double a, double b, double c, double d)
{
this->a = a;
this->b = b;
this->c = c;
this->d = d;
}
__declspec(dllexport) GeoPlane(GeoPoint p0, GeoPoint p1, GeoPoint p2)
{
GeoVector v = GeoVector(p0, p1);
GeoVector u = GeoVector(p0, p2);
GeoVector n = u * v;
// normal vector
this->a = n.GetX();
this->b = n.GetY();
this->c = n.GetZ();
this->d = -(this->a * p0.x + this->b * p0.y + this->c * p0.z);
}
__declspec(dllexport) GeoPlane operator-()
{
return GeoPlane(-this->a, -this->b, -this->c, -this->d);
}
__declspec(dllexport) friend double operator*(const GeoPlane& pl, const GeoPoint& pt)
{
return (pt.x * pl.a + pt.y * pl.b + pt.z * pl.c + pl.d);
}
};
}
GeoFace.h
#pragma once
#include <vector>
#include "GeoPoint.h"
#include "Utility.h"
namespace GeoProc
{
// A Face of a 3D Polygon
class GeoFace
{
// Vertices in one face of the 3D polygon
std::vector<GeoPoint> _pts;
// Vertices Index
std::vector<int> _idx;
// Number of vertices
int _n;
public:
__declspec(dllexport) GeoFace(void);
__declspec(dllexport) ~GeoFace(void)
{
// free memory
Utility::FreeVectorMemory(this->_pts);
Utility::FreeVectorMemory(this->_idx);
}
__declspec(dllexport) GeoFace(std::vector<GeoPoint> pts, std::vector<int> idx)
{
this->_n = pts.size();
for(int i=0;i<_n;i++)
{
this->_pts.push_back(pts[i]);
this->_idx.push_back(idx[i]);
}
}
__declspec(dllexport) int GetN()
{
return this->_n;
}
__declspec(dllexport) std::vector<int> GetI()
{
return this->_idx;
}
__declspec(dllexport) std::vector<GeoPoint> GetV()
{
return this->_pts;
}
};
}
GeoPolygon.h
#pragma once
#include <vector>
#include "GeoPoint.h"
#include "Utility.h"
namespace GeoProc
{
// A 3D Polygon
class GeoPolygon
{
// Vertices of the 3D polygon
std::vector<GeoPoint> _pts;
// Indexes of Vertices
std::vector<int> _idx;
// Number of Vertices
int _n;
public:
__declspec(dllexport) GeoPolygon(void);
__declspec(dllexport) ~GeoPolygon(void)
{
// free memory
Utility::FreeVectorMemory(this->_pts);
Utility::FreeVectorMemory(this->_idx);
}
__declspec(dllexport) GeoPolygon(std::vector<GeoPoint> pts)
{
this->_n = pts.size();
for(int i=0;i<_n;i++)
{
this->_pts.push_back(pts[i]);
this->_idx.push_back(i);
}
}
__declspec(dllexport) int GetN()
{
return this->_n;
}
__declspec(dllexport) std::vector<int> GetI()
{
return this->_idx;
}
__declspec(dllexport) std::vector<GeoPoint> GetV()
{
return this->_pts;
}
};
}
GeoPolygonProc.h
#pragma once
#include <vector>
#include "GeoPoint.h"
#include "GeoVector.h"
#include "GeoFace.h"
#include "GeoPlane.h"
#include "GeoPolygon.h"
#include "Utility.h"
namespace GeoProc
{
// 3D Polygon Process
class GeoPolygonProc
{
// Polygon
GeoPolygon _polygon;
// Polygon Boundary
double _x0, _x1, _y0, _y1, _z0, _z1;
// Polygon faces
std::vector<GeoFace> _Faces;
// Polygon face planes
std::vector<GeoPlane> _FacePlanes;
// Number of faces
int _NumberOfFaces;
// Maximum point to face plane distance error,
// point is considered in the face plane if its distance is less than this error
double _MaxDisError;
__declspec(dllexport) void GeoPolygonProc::Set3DPolygonBoundary();
__declspec(dllexport) void GeoPolygonProc::Set3DPolygonUnitError();
__declspec(dllexport) void GeoPolygonProc::SetConvex3DFaces();
public:
__declspec(dllexport) GeoPolygonProc(void) {}
__declspec(dllexport) ~GeoPolygonProc(void) {}
__declspec(dllexport) GeoPolygonProc(GeoPolygon polygon)
{
this->_polygon = polygon;
Set3DPolygonBoundary();
Set3DPolygonUnitError();
SetConvex3DFaces();
}
__declspec(dllexport) GeoPolygon GetPolygon() { return _polygon; }
__declspec(dllexport) double GetX0() { return this->_x0; }
__declspec(dllexport) double GetX1() { return this->_x1; }
__declspec(dllexport) double GetY0() { return this->_y0; }
__declspec(dllexport) double GetY1() { return this->_y1; }
__declspec(dllexport) double GetZ0() { return this->_z0; }
__declspec(dllexport) double GetZ1() { return this->_z1; }
__declspec(dllexport) std::vector<GeoFace> GetFaces() { return this->_Faces; }
__declspec(dllexport) std::vector<GeoPlane> GetFacePlanes()
{ return this->_FacePlanes; }
__declspec(dllexport) int GetNumberOfFaces() { return this->_NumberOfFaces; }
__declspec(dllexport) double GetMaxDisError() { return this->_MaxDisError; }
__declspec(dllexport) void SetMaxDisError(double value){this->_MaxDisError=value;}
__declspec(dllexport) bool GeoPolygonProc::PointInside3DPolygon(double x, double y,
double z);
};
}
Utility.h
#pragma once
#include <vector>
#include <algorithm>
class Utility
{
public:
Utility(void);
~Utility(void);
template<typename T>
static bool ContainsVector(std::vector<std::vector<T>>
vectorList, std::vector<T> vectorItem)
{
sort(vectorItem.begin(), vectorItem.end());
for (int i = 0; i < vectorList.size(); i++)
{
std::vector<T> temp = vectorList[i];
if (temp.size() == vectorItem.size())
{
sort(temp.begin(), temp.end());
if (equal(temp.begin(), temp.end(), vectorItem.begin()))
{
return true;
}
}
}
return false;
}
template<typename T>
static void FreeVectorMemory(std::vector<T> &obj)
{
obj.clear();
std::vector<T>().swap(obj);
}
template<typename T>
static void FreeVectorListMemory(std::vector<std::vector<T>> &objList)
{
for(int i=0; i<objList.size(); i++)
{
objList[i].clear();
std::vector<T>().swap(objList[i]);
}
objList.clear();
std::vector<std::vector<T>>().swap(objList);
}
};
GeoPolygonProc.cpp
这是引用 GeoProc 库的主类。它能够扩展其功能，基于 GeoProc 库进行更复杂的 3D 多边形几何处理。当前的 GeoPolygonProc 为 3D 多边形边界、面平面和面提供了基本设置。特定的 public 函数 PointInside3DPolygon 是添加其他函数来解决 3D 多边形问题的一个示例。
SetConvex3DFaces 是核心方法，它将多边形表示从完整的顶点集转换为多边形面和面平面。过程如下
- 声明一个 1d
std::vector vertices以从给定的GeoPolygon对象中获取所有顶点。 - 声明一个 2d
std::vector faceVerticeIndex，第一维是面索引，第二维是面顶点索引。 - 声明一个 1d
std:vector fpOutward用于所有面平面。 - 遍历所有给定顶点，用任意 3 个顶点的组合构造一个三角形平面
trianglePlane。 - 通过检查所有其他顶点是否在同一半空间中来确定
trianglePlane是否是面平面。 - 声明一个 1d
std::vector verticeIndexInOneFace用于一个面中完整的顶点索引。 - 查找与三角形平面位于同一平面中的其他顶点，将它们放入 1d 数组
pointInSamePlaneIndex中。 - 通过添加 3 个三角形顶点和在同一平面中的其他顶点，将唯一的面添加到
faceVerticeIndex中。 - 将唯一的面平面添加到
fpOutward - 遍历所有面，获取所有面平面和面。
#include "math.h"
#include "GeoPolygonProc.h"
double MaxUnitMeasureError = 0.001;
namespace GeoProc
{
void GeoPolygonProc::Set3DPolygonBoundary()
{
std::vector<GeoPoint> vertices = _polygon.GetV();
int n = _polygon.GetN();
this->_x0 = vertices[0].x;
this->_x0 = vertices[0].x;
this->_y0 = vertices[0].y;
this->_y1 = vertices[0].y;
this->_z0 = vertices[0].z;
this->_z1 = vertices[0].z;
for(int i=1;i<n;i++)
{
if(vertices[i].x < _x0) this->_x0 = vertices[i].x;
if(vertices[i].y < _y0) this->_y0 = vertices[i].y;
if(vertices[i].z < _z0) this->_z0 = vertices[i].z;
if(vertices[i].x > _x1) this->_x1 = vertices[i].x;
if(vertices[i].y > _y1) this->_y1 = vertices[i].y;
if(vertices[i].z > _z1) this->_z1 = vertices[i].z;
}
}
void GeoPolygonProc::Set3DPolygonUnitError()
{
this->_MaxDisError = ( (fabs(_x0) + fabs(_x1) +
fabs(_y0) + fabs(_y1) +
fabs(_z0) + fabs(_z1) ) / 6 * MaxUnitMeasureError);
}
void GeoPolygonProc::SetConvex3DFaces()
{
// get polygon vertices
std::vector<GeoPoint> vertices = this->_polygon.GetV();
// get number of polygon vertices
int n = this->_polygon.GetN();
// face planes
std::vector<GeoPlane> fpOutward;
// 2d vertices indexes, first dimension is face index,
// second dimension is vertice indexes in one face
std::vector<std::vector<int>> faceVerticeIndex;
for(int i=0; i< n; i++)
{
// triangle point 1
GeoPoint p0 = vertices[i];
for(int j= i+1; j< n; j++)
{
// triangle point 2
GeoPoint p1 = vertices[j];
for(int k = j + 1; k<n; k++)
{
// triangle point 3
GeoPoint p2 = vertices[k];
GeoPlane trianglePlane = GeoPlane(p0, p1, p2);
int onLeftCount = 0;
int onRightCount = 0;
int m = 0;
std::vector<int> pointInSamePlaneIndex;
for(int l = 0; l < n ; l ++)
{
// check any vertices other than the 3 triangle vertices
if(l != i && l != j && l != k)
{
GeoPoint pt = vertices[l];
double dis = trianglePlane * pt;
// add next vertice that is in the triangle plane
if(fabs(dis) < this->_MaxDisError)
{
pointInSamePlaneIndex.push_back(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)
{
std::vector<int> faceVerticeIndexInOneFace;
// add 3 triangle vertices to the triangle plane
faceVerticeIndexInOneFace.push_back(i);
faceVerticeIndexInOneFace.push_back(j);
faceVerticeIndexInOneFace.push_back(k);
// add other same plane vetirces in this triangle plane
for(int p = 0; p < pointInSamePlaneIndex.size(); p ++)
{
faceVerticeIndexInOneFace.push_back(pointInSamePlaneIndex[p]);
}
// check if it is a new face
if(!Utility::ContainsVector
(faceVerticeIndex, faceVerticeIndexInOneFace))
{
// add the new face
faceVerticeIndex.push_back(faceVerticeIndexInOneFace);
// add the new face plane
if(onRightCount == 0)
{
fpOutward.push_back(trianglePlane);
}
else if (onLeftCount == 0)
{
fpOutward.push_back(-trianglePlane);
}
}
// free local memory
Utility::FreeVectorMemory(faceVerticeIndexInOneFace);
}
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
}
// free local memory
Utility::FreeVectorMemory(pointInSamePlaneIndex);
} // k loop
} // j loop
} // i loop
// set number of faces
this->_NumberOfFaces = faceVerticeIndex.size();
for(int i = 0; i<this->_NumberOfFaces; i++)
{
// set face planes
this->_FacePlanes.push_back(GeoPlane(fpOutward[i].a, fpOutward[i].b,
fpOutward[i].c, fpOutward[i].d));
// face vertices
std::vector<GeoPoint> v;
// face vertices indexes
std::vector<int> idx;
// number of vertices of the face
int count = faceVerticeIndex[i].size();
for(int j = 0; j< count; j++)
{
idx.push_back(faceVerticeIndex[i][j]);
v.push_back(GeoPoint(vertices[ idx[j] ].x,
vertices[ idx[j] ].y,
vertices[ idx[j] ].z));
}
// set faces
this->_Faces.push_back(GeoFace(v, idx));
}
// free local memory
Utility::FreeVectorMemory(fpOutward);
Utility::FreeVectorListMemory<int>(faceVerticeIndex);
}
bool GeoPolygonProc::PointInside3DPolygon(double x, double y, double z)
{
GeoPoint pt = GeoPoint(x, y, z);
for (int i=0; i < this->GetNumberOfFaces(); i++)
{
double dis = this->GetFacePlanes()[i] * pt;
// If the point is in the same half space with normal vector
// for any face of the 3D convex polygon, 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 faces,
// then it is inside of the 3D polygon
return true;
}
}
这是 使用 GeoProc DLL 库的 LAS 处理程序
LASProc.cpp
#include <tchar.h>
#include <iostream>
#include <fstream>
#include "GeoPolygonProc.h"
using namespace GeoProc;
GeoPolygonProc GetProcObj()
{
// Initialize cube
GeoPoint CubeVerticesArray[8] =
{
GeoPoint(-27.28046, 37.11775, -39.03485),
GeoPoint(-44.40014, 38.50727, -28.78860),
GeoPoint(-49.63065, 20.24757, -35.05160),
GeoPoint(-32.51096, 18.85805, -45.29785),
GeoPoint(-23.59142, 10.81737, -29.30445),
GeoPoint(-18.36091, 29.07707, -23.04144),
GeoPoint(-35.48060, 30.46659, -12.79519),
GeoPoint(-40.71110, 12.20689, -19.05819)
};
// Create polygon object
std::vector<GeoPoint> verticesVector( CubeVerticesArray,
CubeVerticesArray + sizeof(CubeVerticesArray) / sizeof(GeoPoint) );
GeoPolygon polygonObj = GeoPolygon(verticesVector);
// Create main process object
GeoPolygonProc procObj = GeoPolygonProc(polygonObj);
return procObj;
}
void ProcLAS(GeoPolygonProc procObj)
{
std::ifstream rbFile;
rbFile.open("..\\LASInput\\cube.las",std::ios::binary|std::ios::in);
std::fstream wbFile;
wbFile.open("..\\LASOutput\\cube_inside.las",std::ios::binary|std::ios::out);
std::ofstream wtFile;
wtFile.open("..\\LASOutput\\cube_inside.txt",std::ios::out);
wtFile.precision(4);
// LAS public header variables
unsigned long offset_to_point_data_value;
unsigned long variable_len_num;
unsigned char record_type_c;
unsigned short record_len;
unsigned long record_num;
double x_scale, y_scale, z_scale;
double x_offset, y_offset, z_offset;
// Offset bytes and data types are based on LAS 1.2 Format
// Read LAS public header
rbFile.seekg(96);
rbFile.read ((char *)&offset_to_point_data_value, 4);
rbFile.read ((char *)&variable_len_num, 4);
rbFile.read ((char *)&record_type_c, 1);
rbFile.read ((char *)&record_len, 2);
rbFile.read ((char *)&record_num, 4);
rbFile.seekg(131);
rbFile.read ((char *)&x_scale, 8);
rbFile.read ((char *)&y_scale, 8);
rbFile.read ((char *)&z_scale, 8);
rbFile.read ((char *)&x_offset, 8);
rbFile.read ((char *)&y_offset, 8);
rbFile.read ((char *)&z_offset, 8);
// LAS header buffer
char *bufHeader = (char *)malloc(offset_to_point_data_value);
// Get LAS header
rbFile.seekg(0);
rbFile.read(bufHeader, offset_to_point_data_value);
// Write LAS header
wbFile.write(bufHeader, offset_to_point_data_value);
// LAS point coordinates
double x, y, z; // LAS point actual coordinates
long xi, yi, zi; // LAS point record coordinates
// Number of inside points
int insideCount = 0;
// Point record buffer
char *bufRecord = (char *)malloc(record_len);
// Process point records
for(unsigned int i=0;i<record_num;i++)
{
int record_loc = offset_to_point_data_value + record_len * i;
rbFile.seekg(record_loc);
// Record coordinates
rbFile.read ((char *)&xi, 4);
rbFile.read ((char *)&yi, 4);
rbFile.read ((char *)&zi, 4);
// Actual coordinates
x = (xi * x_scale) + x_offset;
y = (yi * y_scale) + y_offset;
z = (zi * z_scale) + z_offset;
// Filter out the points outside of boundary to reduce computing
if( x>procObj.GetX0() && x<procObj.GetX1() &&
y>procObj.GetY0() && y<procObj.GetY1() &&
z>procObj.GetZ0() && z<procObj.GetZ1())
{
// Main Process to check if the point is inside of the cube
if(procObj.PointInside3DPolygon(x, y, z))
{
// Write the actual coordinates of inside point to text file
wtFile << std::fixed << x << " " <<
std::fixed << y << " " << std::fixed
<< z << std::endl;
// Get point record
rbFile.seekg(record_loc);
rbFile.read(bufRecord, record_len);
// Write point record to binary LAS file
wbFile.write(bufRecord, record_len);
insideCount++;
}
}
}
// Update total record numbers in output binary LAS file
wbFile.seekp(107);
wbFile.write((char *)&insideCount, 4);
rbFile.close();
wbFile.close();
wtFile.close();
free(bufHeader);
free(bufRecord);
}
int _tmain(int argc, _TCHAR* argv[])
{
// Create procInst
GeoPolygonProc procObj = GetProcObj();
// Process LAS
ProcLAS(procObj);
return 0;
}
GeoProcTest.cpp
测试项目旨在确保 GeoProc DLL 库中的所有功能都能正常工作。
#include <tchar.h>
#include <iostream>
#include "GeoPolygonProc.h"
using namespace GeoProc;
GeoPoint p1 = GeoPoint( - 27.28046, 37.11775, - 39.03485);
GeoPoint p2 = GeoPoint( - 44.40014, 38.50727, - 28.78860);
GeoPoint p3 = GeoPoint( - 49.63065, 20.24757, - 35.05160);
GeoPoint p4 = GeoPoint( - 32.51096, 18.85805, - 45.29785);
GeoPoint p5 = GeoPoint( - 23.59142, 10.81737, - 29.30445);
GeoPoint p6 = GeoPoint( - 18.36091, 29.07707, - 23.04144);
GeoPoint p7 = GeoPoint( - 35.48060, 30.46659, - 12.79519);
GeoPoint p8 = GeoPoint( - 40.71110, 12.20689, - 19.05819);
GeoPoint verticesArray[8] = { p1, p2, p3, p4, p5, p6, p7, p8 };
std::vector<GeoPoint> verticesVector( verticesArray,
verticesArray + sizeof(verticesArray) / sizeof(GeoPoint) );
GeoPolygon polygon = GeoPolygon(verticesVector);
void GeoPoint_Test()
{
std::cout << "GeoPoint Test:" << std::endl;
GeoPoint pt = p1 + p2;
std::cout << pt.x << ", " <<
pt.y << ", " << pt.z << std::endl;
}
void GeoVector_Test()
{
std::cout << "GeoVector Test:" << std::endl;
GeoVector v1 = GeoVector(p1, p2);
GeoVector v2 = GeoVector(p1, p3);
GeoVector v3 = v2 * v1;
std::cout << v3.GetX() << ", " <<
v3.GetY() << ", " << v3.GetZ() << std::endl;
}
void GeoPlane_Test()
{
std::cout << "GeoPlane Test:" << std::endl;
GeoPlane pl = GeoPlane(p1, p2, p3);
std::cout << pl.a << ", " << pl.b << ",
" << pl.c << ", " << pl.d << std::endl;
pl = GeoPlane(1.0, 2.0, 3.0, 4.0);
std::cout << pl.a << ", " << pl.b << ",
" << pl.c << ", " << pl.d << std::endl;
pl = -pl;
std::cout << pl.a << ", " << pl.b << ",
" << pl.c << ", " << pl.d << std::endl;
double dis = pl * GeoPoint(1.0, 2.0, 3.0);
std::cout << dis << std::endl;
}
void GeoPolygon_Test()
{
std::cout << "GeoPolygon Test:" << std::endl;
std::vector<int> idx = polygon.GetI();
std::vector<GeoPoint> v = polygon.GetV();
for(int i=0; i<polygon.GetN(); i++)
{
std::cout << idx[i] << ": (" << v[i].x << ",
" << v[i].y << ", "
<< v[i].z << ")" << std::endl;
}
}
void GeoFace_Test()
{
std::cout << "GeoFace Test:" << std::endl;
GeoPoint faceVerticesArray[4] = { p1, p2, p3, p4 };
std::vector<GeoPoint> faceVerticesVector( faceVerticesArray,
faceVerticesArray + sizeof(faceVerticesArray) / sizeof(GeoPoint) );
int faceVerticeIndexArray[4] = { 1, 2, 3, 4 };
std::vector<int> faceVerticeIndexVector( faceVerticeIndexArray,
faceVerticeIndexArray + sizeof(faceVerticeIndexArray) / sizeof(int) );
GeoFace face = GeoFace(faceVerticesVector, faceVerticeIndexVector);
std::vector<int> idx = face.GetI();
std::vector<GeoPoint> v = face.GetV();
for(int i=0; i<face.GetN(); i++)
{
std::cout << idx[i] << ": (" << v[i].x << ",
" << v[i].y << ", "
<< v[i].z << ")" << std::endl;
}
}
void Utility_Test()
{
std::cout << "Utility Test:" << std::endl;
int arr1[4] = { 1, 2, 3, 4 };
std::vector<int> arr1Vector( arr1, arr1 + sizeof(arr1) / sizeof(int) );
int arr2[4] = { 4, 5, 6, 7 };
std::vector<int> arr2Vector( arr2, arr2 + sizeof(arr2) / sizeof(int) );
std::vector<std::vector<int>> vector_2d;
vector_2d.push_back(arr1Vector);
vector_2d.push_back(arr2Vector);
int arr3[4] = { 2, 3, 1, 4 };
std::vector<int> arr3Vector( arr3, arr3 + sizeof(arr3) / sizeof(int) );
int arr4[4] = { 2, 3, 1, 5 };
std::vector<int> arr4Vector( arr4, arr4 + sizeof(arr4) / sizeof(int) );
bool b1 = Utility::ContainsVector(vector_2d, arr3Vector);
bool b2 = Utility::ContainsVector(vector_2d, arr4Vector);
std::cout << b1 << ", " << b2 << std::endl;
std::cout << arr1Vector.size() << std::endl;
Utility::FreeVectorMemory<int>(arr1Vector);
std::cout << arr1Vector.size() << std::endl;
std::cout << vector_2d.size() << std::endl;
Utility::FreeVectorListMemory<int>(vector_2d);
std::cout << vector_2d.size() << std::endl;
}
void GeoPolygonProc_Test()
{
std::cout << "GeoPolygonProc Test:" << std::endl;
GeoPolygonProc procObj = GeoPolygonProc(polygon);
std::cout << procObj.GetX0() << ", " <<
procObj.GetX1() << ", " <<
procObj.GetY0() << ", " <<
procObj.GetY1() << ", " <<
procObj.GetZ0() << ", " <<
procObj.GetZ1() << ", " << std::endl;
std::cout << procObj.GetMaxDisError() << std::endl;
procObj.SetMaxDisError(0.0125);
std::cout << procObj.GetMaxDisError() << std::endl;
int count = procObj.GetNumberOfFaces();
std::vector<GeoPlane> facePlanes = procObj.GetFacePlanes();
std::vector<GeoFace> faces = procObj.GetFaces();
std::cout << count << std::endl;
std::cout << "Face Planes:" << std::endl;
for(int i=0; i<count; i++)
{
std::cout << facePlanes[i].a << ", " <<
facePlanes[i].b << ", " <<
facePlanes[i].a << ", " <<
facePlanes[i].d << std::endl;
}
std::cout << "Faces:" << std::endl;
for(int i=0; i<count; i++)
{
std::cout << "Face #" << i + 1 << " :" <<std::endl;
GeoFace face = faces[i];
std::vector<int> idx = face.GetI();
std::vector<GeoPoint> v = face.GetV();
for(int i=0; i<face.GetN(); i++)
{
std::cout << idx[i] << ": (" << v[i].x << ",
" << v[i].y << ", "
<< v[i].z << ")" << std::endl;
}
}
GeoPoint insidePoint = GeoPoint( - 28.411750, 25.794500, - 37.969000);
GeoPoint outsidePoint = GeoPoint( - 28.411750, 25.794500, - 50.969000);
bool b1 = procObj.PointInside3DPolygon(insidePoint.x, insidePoint.y, insidePoint.z);
bool b2 = procObj.PointInside3DPolygon(outsidePoint.x, outsidePoint.y, outsidePoint.z);
std::cout << b1 << ", " << b2 << std::endl;
}
int _tmain(int argc, _TCHAR* argv[])
{
GeoPoint_Test();
GeoVector_Test();
GeoPlane_Test();
GeoPolygon_Test();
GeoFace_Test();
Utility_Test();
GeoPolygonProc_Test();
}
以下是测试项目输出
GeoPoint Test:
-71.6806, 75.625, -67.8234
GeoVector Test:
-178.391, 160.814, -319.868
GeoPlane Test:
-178.391, 160.814, -319.868, -23321.6
1, 2, 3, 4
-1, -2, -3, -4
-18
GeoPolygon Test:
0: (-27.2805, 37.1178, -39.0348)
1: (-44.4001, 38.5073, -28.7886)
2: (-49.6307, 20.2476, -35.0516)
3: (-32.511, 18.858, -45.2978)
4: (-23.5914, 10.8174, -29.3044)
5: (-18.3609, 29.0771, -23.0414)
6: (-35.4806, 30.4666, -12.7952)
7: (-40.7111, 12.2069, -19.0582)
GeoFace Test:
1: (-27.2805, 37.1178, -39.0348)
2: (-44.4001, 38.5073, -28.7886)
3: (-49.6307, 20.2476, -35.0516)
4: (-32.511, 18.858, -45.2978)
Utility Test:
1, 0
4
0
2
0
GeoPolygonProc Test:
-49.6307, -18.3609, 10.8174, 38.5073, -45.2978, -12.7952,
0.0292349
0.0125
6
Face Planes:
-178.391, 160.814, -178.391, -23321.6
104.61, 365.194, 104.61, -5811.87
342.393, -27.7904, 342.393, 2372.96
-342.394, 27.7906, -342.394, -10373
-104.61, -365.194, -104.61, -2188.14
178.391, -160.814, 178.391, 15321.6
Faces:
Face #1 :
0: (-27.2805, 37.1178, -39.0348)
1: (-44.4001, 38.5073, -28.7886)
2: (-49.6307, 20.2476, -35.0516)
3: (-32.511, 18.858, -45.2978)
Face #2 :
0: (-27.2805, 37.1178, -39.0348)
1: (-44.4001, 38.5073, -28.7886)
5: (-18.3609, 29.0771, -23.0414)
6: (-35.4806, 30.4666, -12.7952)
Face #3 :
0: (-27.2805, 37.1178, -39.0348)
3: (-32.511, 18.858, -45.2978)
4: (-23.5914, 10.8174, -29.3044)
5: (-18.3609, 29.0771, -23.0414)
Face #4 :
1: (-44.4001, 38.5073, -28.7886)
2: (-49.6307, 20.2476, -35.0516)
6: (-35.4806, 30.4666, -12.7952)
7: (-40.7111, 12.2069, -19.0582)
Face #5 :
2: (-49.6307, 20.2476, -35.0516)
3: (-32.511, 18.858, -45.2978)
4: (-23.5914, 10.8174, -29.3044)
7: (-40.7111, 12.2069, -19.0582)
Face #6 :
4: (-23.5914, 10.8174, -29.3044)
5: (-18.3609, 29.0771, -23.0414)
6: (-35.4806, 30.4666, -12.7952)
7: (-40.7111, 12.2069, -19.0582)
1, 0
关注点
该算法不需要耗时的数学计算，即三角形函数、平方根等。其性能良好。虽然该算法不完全适用于 3D 凹多边形，但如果您能够以某种方式（例如手动）为 3D 凹多边形提供所有面的向外法向量，那么该算法仍然可以在一半的程度上可行。
历史
2016 年 1 月 21 日：第三版日期
- 添加了一个新项目
GeoProcTest来测试GeoProcDLL 库中的所有功能 - 所有向量都加上了
std::前缀，并删除了所有using namespace std - 修复了函数
FreeVectorMemory和FreeVectorListMemory中的一个错误，通过传递引用作为参数来修改函数内部的值
2016 年 1 月 19 日：第二版日期
- 使用 STL
std::vector替换new运算符进行动态数组分配。std::vector在操作数组方面具有一些优点。它不需要大小来分配数组，而 new 运算符需要；它能够通过索引访问数组，而std::list不能。 - 删除了 VC++ 预编译头以最小化对 Windows 的依赖
- 在 LAS 文件读写过程中，使用 STL 文件流函数，即
ifstream、seekg、read等，替换 C 库函数，即FILE、fread、fwrite等。