15,662,970 members
Articles / Programming Languages / C++
Article
Posted 19 Dec 2015

82.3K views
58 bookmarked

# Point Inside 3D Convex Polygon in C++

Rate me:
An algorithm to determine if a point is inside a 3D convex polygon for a given polygon vertices in C++

## Introduction

This article implements an algorithm to utilize plane normal vector and direction of point to plane distance vector to determine if a point is inside a 3D convex polygon for a given polygon vertices.

This is the original C++ version, I already ported the algorithm to C# version, Java version, JavaScript version, PHP version, Python version, Perl version and Fortran.

These versions cover different programming types, from compiled language to interpreted language, all support Object Oriented programming, which makes these series portings have a clear, stable and consistent common program architecture.

A MASM Assembly version is now available, although Assembly language has no OO feature, but with the help of its STRUCT and invoke mechanism, we can use a library based architecture to make the program to have same structure with all above advanced languages. This should be the last one for this series portings.

## Background

I had a chance to get a glimpse for the point cloud process in laser scanning 3D field. One problem in the field is how to separate the points with a certain geometry boundary. For example, with a given 8 cube vertices, how to get all points inside the cube.

A few approaches can be found online to solve the problem. For example, Ray Tracing, which checks if the intersects number is odd to determine if the point is inside. Cartesian coordinate transform is also used in some special cases, for example, if the boundary is a perfect cube without any distortion.

I ended up this investigation with a simple way to solve it. The idea is mainly inspired by the reference article to discuss point-plane distance: http://mathworld.wolfram.com/Point-PlaneDistance.html.

Here is the question and my solution.

### Question

For a given 3D convex polygon with N vertices, determine if a 3D point (x, y, z) is inside the polygon.

### Solution

A 3D convex polygon has many faces, a face has a face plane where the face lies in.

A face plane has an outward normal vector, which directs to outside of the polygon.

A point to face plane distance defines a geometry vector, if the distance vector has an opposite direction with the outward normal vector, then the point is in "inside half space" of the face plane, otherwise, it is in "outside half space" of the face plane.

A point is determined to be inside of the 3D polygon if the point is in "inside half space" for all faces of the 3D convex polygon.

That is the basic idea of this algorithm.

Algorithm geometry diagram:

There is another question to be solved prior to utilizing the above method, we need to get face plane outward normal vectors for all faces of the polygon. Here are steps to accomplish this:

1. Get any 3 vertices combinations from N vertices, which is prepared to construct a triangle plane.
2. Get normal vector `n = (A, B, C)` with the cross production of 2 edge vectors `u` and `v` from the 3 vertices `P0`, `P1` and `P2` triangle. The vertice `P0(x0, y0, z0)` is the common begin point of the 2 edge vectors.
3. Get the face plane with equation `A * x + B * y + C * z + D = 0`. Here `D = -(A * x0 + B * y0 + c * z0)`.
4. Determine if the triangle plane is a face plane with checking if all other vertices points are in the same half space of the triangle plane. This step requires the convex assumption. For a concave 3D polygon, it is not distinguished between a polygon real face and a polygon inner triangle with this rule.
5. Get the face vertices with appending any other vertices that lie in the same triangle plane to these 3 vertices.
6. After all, we find all faces with their outward normal vectors and their complete vertices.

Here is the essential formula:

For a point `pt(X, Y, Z)` and a plane `pl(A, B, C, D)`, the sign of value `(pt.X * pl.A + pt.Y * pl.B + pt.Z * pl.C + pl.D)` is able to determine which half space the point is in.

## Using the Code

The algorithm is wrapped into a VC++ DLL `GeoProc`. The main test program `LASProc` 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 algorithm DLL, construct a polygon object first like this:

C++
```// Create polygon object
GeoPolygon polygonObj = GeoPolygon(verticesVector);```

Then, construct the main process object:

C++
```// Create main process object
GeoPolygonProc procObj = GeoPolygonProc(polygonObj);```

Main procedure to check if a point (`x, y, z`) is inside the `CubeVertices`:

C++
```// Main procedure to check if a point is inside polygon
procObj.PointInside3DPolygon(x, y, z)```

Here are the classes in the GeoProc DLL library:

The codes are almost self-evident with comments.

### GeoPoint.h

C++
```#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

C++
```#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

C++
```#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

C++
```#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

C++
```#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

C++
```#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

C++
```#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

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 an example for adding other functions to solve 3D polygon problems.

`SetConvex3DFaces` is the essential method, it transforms the polygon presentation from a complete polygon vertice set to polygon faces and faces planes. The procedure is as given below:

1. Declare a 1d `std::vector vertices `to get all vertices from the given `GeoPolygon` object.
2. Declare a 2d `std::vector faceVerticeIndex`, first dimension is face index, second dimension is face vertice index.
3. Declare a 1d `std:vector fpOutward `for all face planes.
4. Go through all given vertices, construct a triangle plane `trianglePlane `with any 3 vertices combination.
5. Determine if `trianglePlane `is a face plane through checking if all other vertices are in the same half space.
6. Declare a 1d `std::vector verticeIndexInOneFace `for a complete vertices indexes in one face.
7. Find other vertices in the same plane with the triangle plane, put them in 1d array `pointInSamePlaneIndex`.
8. Add unique face to `faceVerticeIndex `with adding the 3 triangle vertices and the other same plane vertices.
9. Add unique face plane to `fpOutward`
10. Go through all faces and get all face planes and faces.
C++
```#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))
{
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;
}
}```

Here is the LAS Process program with using GeoProc DLL library:

### LASProc.cpp

C++
```#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);

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

rbFile.seekg(96);

rbFile.seekg(131);

rbFile.seekg(0);

// 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

// 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);

// 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(bufRecord);
}

int _tmain(int argc, _TCHAR* argv[])
{
// Create procInst
GeoPolygonProc procObj = GetProcObj();

// Process LAS
ProcLAS(procObj);

return 0;
}```

### GeoProcTest.cpp

The test project is to ensure all functions in `GeoProc` DLL library works.

C++
```#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();
}```

Here is the test project output:

```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```

## Points of Interest

The algorithm does not need time consuming math computing, i.e., triangle function, square root, etc. Its performance is good. Although the algorithm is not suitable for a 3D concave polygon in a whole way, but if you can provide all face outward normal vector for a 3D concave polygon in some way, i.e., manually, then the algorithm is still feasible in a half way.

## History

##### January 21st, 2016: Third version date
1. Added a new project `GeoProcTest` to test all functions in `GeoProc` DLL library
2. Prefixed all vector with `std::` and removed all using `namespace std`
3. Fixed a bug in function `FreeVectorMemory` and `FreeVectorListMemory` with passing reference as argument to modify them inside the functions
##### January 19th, 2016: Second version date
1. Used STL `std::vector` to replace the `new` operator for dynamic array allocation. `std::vector` has a few advantages regarding how to manipulate array. It does not need a size to allocate array, while new operator does; It is capable of accessing array by index, while `std::list` cannot.
2. Removed the VC++ Precompiled Header to minimize the Windows dependency
3. Used STL file stream functions, i.e., `ifstream`, `seekg`, `read`, etc. to replace C library functions, i.e., `FILE`, `fread`, `fwrite`, etc. in the LAS file read-write process

## Reference

Written By
Software Developer
My name is Jiyang Hou (or John Hou). I was born in HeiLongJiang province in north east of China. I got all my educations in China. My university major is Geophysics, but my main professional role is software developer. My biggest accomplishment so far is quit smoking about 5 years ago after almost 20 years smoking history. I am still interested on programming beside making living with it like many other developers. I immigrated to Canada in 2003 and became a permanent resident till now. I live in Calgary, Alberta, Canada. You can reach me by jyhou69@gmail.com regarding to any questions, comments, advice, etc.

 First Prev Next
 Is it typo? Koon Tae Ki21-May-18 20:43 Koon Tae Ki 21-May-18 20:43
 btw JustWatchLittle 15-Dec-17 12:17 JustWatchLittle 15-Dec-17 12:17
 Non-exhaustive approach YvesDaoust7-Mar-16 23:58 YvesDaoust 7-Mar-16 23:58
 please avoid prefix '_' in C/C++ symbols Stefan_Lang28-Jan-16 21:36 Stefan_Lang 28-Jan-16 21:36
 As a DLL The_Inventor15-Jan-16 18:27 The_Inventor 15-Jan-16 18:27
 Re: As a DLL sisira22-Jan-16 13:27 sisira 22-Jan-16 13:27
 Re: As a DLL The_Inventor22-Jan-16 16:09 The_Inventor 22-Jan-16 16:09
 I believe that your reply to my Praise to the author of this particular article was meant for a different question / comment. Once the plane and normal of the surface has been defined then +1 , 0 , -1 , depending on how one writes the logic algorithm. The World as we think we know it Has a lot more to it than meets the eye. A Mad Scientist who has seen it for himself....
 Re: As a DLL Rick York20-Apr-17 10:07 Rick York 20-Apr-17 10:07
 Re: As a DLL The_Inventor20-Apr-17 17:10 The_Inventor 20-Apr-17 17:10
 getting the direction of the normal Member 73165214-Jan-16 10:51 Member 731652 14-Jan-16 10:51
 Polyhedron, not polygon Stefan_Lang4-Jan-16 3:35 Stefan_Lang 4-Jan-16 3:35
 What kind of task could be solved using this algorithm? Alexander Voronin21-Dec-15 22:15 Alexander Voronin 21-Dec-15 22:15
 Re: What kind of task could be solved using this algorithm? `Randor` 27-Sep-17 13:18 `Randor` 27-Sep-17 13:18
 Alternative way LucaBarbucci21-Dec-15 0:42 LucaBarbucci 21-Dec-15 0:42
 Re: Alternative way Stefan_Lang5-Jan-16 1:38 Stefan_Lang 5-Jan-16 1:38
 Last Visit: 31-Dec-99 18:00     Last Update: 28-May-23 20:02 Refresh 1