Click here to Skip to main content
13,248,107 members (35,239 online)
Click here to Skip to main content
Add your own
alternative version

Stats

41.7K views
2K downloads
53 bookmarked
Posted 19 Dec 2015

Point Inside 3D Convex Polygon in C++

, 2 Mar 2016
Rate this:
Please Sign up or sign in to vote.
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:

// Create polygon object
GeoPolygon polygonObj = GeoPolygon(verticesVector);

Then, construct the main process object:

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

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

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

#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

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.
#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;
    }
}

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

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

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

#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.
December 19th, 2015: Initial version date

Reference

License

This article, along with any associated source code and files, is licensed under The Code Project Open License (CPOL)

Share

About the Author

John Jiyang Hou
Software Developer
Canada Canada
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.

You may also be interested in...

Pro

Comments and Discussions

 
QuestionNon-exhaustive approach Pin
YvesDaoust8-Mar-16 0:58
memberYvesDaoust8-Mar-16 0:58 
Suggestionplease avoid prefix '_' in C/C++ symbols Pin
Stefan_Lang28-Jan-16 22:36
memberStefan_Lang28-Jan-16 22:36 
PraiseAs a DLL Pin
The_Inventor15-Jan-16 19:27
memberThe_Inventor15-Jan-16 19:27 
GeneralRe: As a DLL Pin
sisira22-Jan-16 14:27
membersisira22-Jan-16 14:27 
GeneralRe: As a DLL Pin
The_Inventor22-Jan-16 17:09
memberThe_Inventor22-Jan-16 17:09 
GeneralRe: As a DLL Pin
Rick York20-Apr-17 11:07
memberRick York20-Apr-17 11:07 
GeneralRe: As a DLL Pin
The_Inventor20-Apr-17 18:10
memberThe_Inventor20-Apr-17 18:10 
Questiongetting the direction of the normal Pin
Member 73165214-Jan-16 11:51
memberMember 73165214-Jan-16 11:51 
QuestionPolyhedron, not polygon Pin
Stefan_Lang4-Jan-16 4:35
memberStefan_Lang4-Jan-16 4:35 
QuestionWhat kind of task could be solved using this algorithm? Pin
Alexander Voronin21-Dec-15 23:15
memberAlexander Voronin21-Dec-15 23:15 
AnswerRe: What kind of task could be solved using this algorithm? Pin
Randor 27-Sep-17 14:18
professional Randor 27-Sep-17 14:18 
AnswerAlternative way Pin
LucaBarbucci21-Dec-15 1:42
memberLucaBarbucci21-Dec-15 1:42 
GeneralRe: Alternative way Pin
Stefan_Lang5-Jan-16 2:38
memberStefan_Lang5-Jan-16 2:38 

General General    News News    Suggestion Suggestion    Question Question    Bug Bug    Answer Answer    Joke Joke    Praise Praise    Rant Rant    Admin Admin   

Use Ctrl+Left/Right to switch messages, Ctrl+Up/Down to switch threads, Ctrl+Shift+Left/Right to switch pages.

Permalink | Advertise | Privacy | Terms of Use | Mobile
Web01 | 2.8.171114.1 | Last Updated 3 Mar 2016
Article Copyright 2015 by John Jiyang Hou
Everything else Copyright © CodeProject, 1999-2017
Layout: fixed | fluid