Point Inside 3D Convex Polygon in Java

• Download FugroViewerSetup22 - 2.9 MB
• Download JGeoProc - 5.1 MB

Introduction

This is a Java version of the original C++ algorithm which can be found here.

Background

Please refer to the original C++ algorithm here.

Using the Code

The algorithm is wrapped into a Java class library package GeoProc. The main test program CSLASProc 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 class library, construct a polygon instance first like this:

// Create polygon instance
GeoPolygon polygonInst = new GeoPolygon(CubeVertices);

Then, construct the main process instance:

// Create main process instance
GeoPolygonProc procInst = new GeoPolygonProc(polygonInst);

Main procedure to check if a point (ax, ay, az) is inside the CubeVertices:

// Main Process to check if the point is inside of the cube                
procInst.PointInside3DPolygon(ax, ay, az)

Here are the classes in the GeoProc package:

GeoPoint.java

public class GeoPoint {
    
    private double x;
    private double y;
    private double z;
    
    public double getX() { return x; }    
    public void setX(double x) { this.x = x;}
    
    public double getY() { return y; }    
    public void setY(double y) { this.y = y;}
    
    public double getZ() { return z; }    
    public void setZ(double z) { this.z = z;}
    
    public GeoPoint(){}
    
    public GeoPoint(double x, double y, double z)        
    {
        this.x=x;
        this.y=y;
        this.z=z;    
    }    

    public static GeoPoint Add(GeoPoint p0, GeoPoint p1)
    {
        return new GeoPoint(p0.x + p1.x, p0.y + p1.y, p0.z + p1.z);
    }
}

GeoVector.java

public class GeoVector {
    
    private GeoPoint p0; // vector begin point
    private GeoPoint p1; // vector end point
    private double x; // vector x axis projection value
    private double y; // vector y axis projection value
    private double z; // vector z axis projection value
    
    public GeoPoint getP0() {return this.p0;}
    public GeoPoint getP1() {return this.p1;}    
    public double getX() {return this.x;}
    public double getY() {return this.y;}
    public double getZ() {return this.z;}  

    public GeoVector() {}
        
    public GeoVector(GeoPoint p0, GeoPoint p1)
    {
        this.p0 = p0;
        this.p1 = p1;            
        this.x = p1.getX() - p0.getX();
        this.y = p1.getY() - p0.getY();
        this.z = p1.getZ() - p0.getZ();
    }        

    public static GeoVector Multiple(GeoVector u, GeoVector v)
    {
        double x = u.getY() * v.getZ() - u.getZ() * v.getY();
        double y = u.getZ() * v.getX() - u.getX() * v.getZ();
        double z = u.getX() * v.getY() - u.getY() * v.getX();
        
        GeoPoint p0 = v.getP0();
        GeoPoint p1 = GeoPoint.Add(p0, new GeoPoint(x, y, z));

        return new GeoVector(p0, p1);
    }
}

GeoPlane.java

public class GeoPlane {
    
    // Plane Equation: a * x + b * y + c * z + d = 0

    private double a;
    private double b;
    private double c;
    private double d;
    
    public double getA() { return this.a; }
    public double getB() { return this.b; }
    public double getC() { return this.c; }
    public double getD() { return this.d; }

    public GeoPlane() {}
    
    public GeoPlane(double a, double b, double c, double d)
    {
        this.a = a;
        this.b = b;
        this.c = c;
        this.d = d;
    }

    public GeoPlane(GeoPoint p0, GeoPoint p1, GeoPoint p2)
    {        
        GeoVector v = new GeoVector(p0, p1);

        GeoVector u = new GeoVector(p0, p2);

        GeoVector n = GeoVector.Multiple(u, v);

        // normal vector        
        double a = n.getX();
        double b = n.getY();
        double c = n.getZ();                
        double d = -(a * p0.getX() + b * p0.getY() + c * p0.getZ());
        
        this.a = a;
        this.b = b;
        this.c = c;
        this.d = d;        
    }

    public static GeoPlane Negative(GeoPlane pl)
    {
        return new GeoPlane(-pl.getA(), -pl.getB(), -pl.getC(), -pl.getD());
    }

    public static double Multiple(GeoPoint pt, GeoPlane pl)
    {
        return (pt.getX() * pl.getA() + pt.getY() * pl.getB() + 
                pt.getZ() * pl.getC() + pl.getD());
    }
}

GeoFace.java

public class GeoFace {
    
    // Vertices in one face of the 3D polygon
    private ArrayList<GeoPoint> v;
    
    // Vertices Index
    private ArrayList<Integer> idx;

    // Number of vertices    
    private int n;        

    public ArrayList<GeoPoint> getV() { return this.v; }

    public ArrayList<Integer> getI() { return this.idx; }
    
    public int getN() { return this.n; }

    public GeoFace(){}
    
    public GeoFace(ArrayList<GeoPoint> p, ArrayList<Integer> idx)
    {            
        this.v = new ArrayList<GeoPoint>();

        this.idx = new ArrayList<Integer>();

        this.n = p.size();
        
        for(int i=0;i<n;i++)
        {
            this.v.add(p.get(i));
            this.idx.add(idx.get(i));
        }        
    }   
}

GeoPolygon.java

public class GeoPolygon {
    
    // Vertices of the 3D polygon
    private ArrayList<GeoPoint> v;
    
    // Vertices Index
    private ArrayList<Integer> idx;

    // Number of vertices    
    private int n;        

    public ArrayList<GeoPoint> getV() { return this.v; }

    public ArrayList<Integer> getI() { return this.idx; }
    
    public int getN() { return this.n; }

    public GeoPolygon(){}
    
    public GeoPolygon(ArrayList<GeoPoint> p)
    {            
        this.v = new ArrayList<GeoPoint>();

        this.idx = new ArrayList<Integer>();

        this.n = p.size();
        
        for(int i=0;i<n;i++)
        {
            this.v.add(p.get(i));
            this.idx.add(i);
        }        
    }   
}

GeoPolygonProc.java

public class GeoPolygonProc {
    
    private double MaxUnitMeasureError = 0.001;

    // Polygon Boundary
    private double X0, X1, Y0, Y1, Z0, Z1;    

    // Polygon faces
    private ArrayList<GeoFace> Faces;

    // Polygon face planes
    private ArrayList<GeoPlane> FacePlanes;

    // Number of faces
    private int NumberOfFaces;    

    // Maximum point to face plane distance error,
    // point is considered in the face plane if its distance is less than this error
    private double MaxDisError;

    public double getX0() { return this.X0; }
    public double getX1() { return this.X1; }
    public double getY0() { return this.Y0; }
    public double getY1() { return this.Y1; }
    public double getZ0() { return this.Z0; }
    public double getZ1() { return this.Z1; }
    public ArrayList<GeoFace> getFaces() { return this.Faces; }
    public ArrayList<GeoPlane> GetFacePlanes() { return this.FacePlanes; }
    public int getNumberOfFaces() { return this.NumberOfFaces; }
    
    public GeoPolygonProc(){}        
        
    public GeoPolygonProc(GeoPolygon polygonInst)
    {
        
        // Set boundary
        this.Set3DPolygonBoundary(polygonInst);

        // Set maximum point to face plane distance error,         
        this.Set3DPolygonUnitError(polygonInst);

        // Set faces and face planes        
        this.SetConvex3DFaces(polygonInst);      
    }    

    public boolean PointInside3DPolygon(double x, double y, double z)
    {
        GeoPoint P = new GeoPoint(x, y, z);
        
        for (int i = 0; i < this.NumberOfFaces; i++)
        {

            double dis = GeoPlane.Multiple(P, this.FacePlanes.get(i));

            // If the point is in the same half space with normal vector for any face of the cube,
            // then it is outside of the 3D polygon        
            if (dis > 0)
            {
                return false;
            }
        }

        // If the point is in the opposite half space with normal vector for all 6 faces,
        // then it is inside of the 3D polygon
        return true;            
    }

    private void Set3DPolygonUnitError(GeoPolygon polygon)
    {        
        this.MaxDisError = ((Math.abs(this.X0) + Math.abs(this.X1) +
                Math.abs(this.Y0) + Math.abs(this.Y1) +
                Math.abs(this.Z0) + Math.abs(this.Z1)) / 6 * MaxUnitMeasureError);        
    }

    private void Set3DPolygonBoundary(GeoPolygon polygon)
    {
        ArrayList<GeoPoint> vertices = polygon.getV();

        int n = polygon.getN();
        
        double xmin, xmax, ymin, ymax, zmin, zmax;

        xmin = xmax = vertices.get(0).getX();
        ymin = ymax = vertices.get(0).getY();
        zmin = zmax = vertices.get(0).getZ();

        for (int i = 1; i < n; i++)
        {
            if (vertices.get(i).getX() < xmin) xmin = vertices.get(i).getX();
            if (vertices.get(i).getY() < ymin) ymin = vertices.get(i).getY();
            if (vertices.get(i).getZ() < zmin) zmin = vertices.get(i).getZ();
            if (vertices.get(i).getX() > xmax) xmax = vertices.get(i).getX();
            if (vertices.get(i).getY() > ymax) ymax = vertices.get(i).getY();
            if (vertices.get(i).getZ() > zmax) zmax = vertices.get(i).getZ();
        }        
        
        this.X0 = xmin;
        this.X1 = xmax;
        this.Y0 = ymin;
        this.Y1 = ymax;
        this.Z0 = zmin;
        this.Z1 = zmax;
    }
    
    private void SetConvex3DFaces(GeoPolygon polygon)
    {        
        ArrayList<GeoFace> faces = new ArrayList<GeoFace>();

        ArrayList<GeoPlane> facePlanes = new ArrayList<GeoPlane>();
        
        int numberOfFaces;
        
        double maxError = this.MaxDisError;
        
        // vertices of 3D polygon
        ArrayList<GeoPoint> vertices = polygon.getV();
      
        int n = polygon.getN();

        // vertices indexes for all faces
        // vertices index is the original index value in the input polygon
        ArrayList<ArrayList<Integer>> faceVerticeIndex = new ArrayList<ArrayList<Integer>>();
        
        // face planes for all faces
        ArrayList<GeoPlane> fpOutward = new ArrayList<GeoPlane>();    
            
        for(int i=0; i< n; i++)
        {
            // triangle point 1
            GeoPoint p0 = vertices.get(i);

            for(int j= i+1; j< n; j++)
            {
                // triangle point 2
                GeoPoint p1 = vertices.get(j);

                for(int k = j + 1; k<n; k++)
                {
                    // triangle point 3
                    GeoPoint p2 = vertices.get(k);
    
                    GeoPlane trianglePlane = new GeoPlane(p0, p1, p2);
                
                    int onLeftCount = 0;
                    int onRightCount = 0;

                    // indexes of points that lie in same plane with face triangle plane
                    ArrayList<Integer> pointInSamePlaneIndex = new ArrayList<Integer>();
        
                    for(int l = 0; l < n ; l ++)
                    {                        
                        // any point other than the 3 triangle points
                        if(l != i && l != j && l != k)
                        {
                            GeoPoint p = vertices.get(l);

                            double dis = GeoPlane.Multiple(p, trianglePlane);
                            
                            // next point is in the triangle plane
                            if(Math.abs(dis) < maxError)
                            {                            
                                pointInSamePlaneIndex.add(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)
                    {
                        ArrayList<Integer> verticeIndexInOneFace = new ArrayList<Integer>();
                       
                        // triangle plane
                        verticeIndexInOneFace.add(i);
                        verticeIndexInOneFace.add(j);
                        verticeIndexInOneFace.add(k);
                        
                        int m = pointInSamePlaneIndex.size();

                        if(m > 0) // there are other vertices in this triangle plane
                        {
                            for(int p = 0; p < m; p ++)
                            {
                                verticeIndexInOneFace.add(pointInSamePlaneIndex.get(p));
                            }                        
                        }

                        // if verticeIndexInOneFace is a new face,
                        // add it in the faceVerticeIndex list,
                        // add the trianglePlane in the face plane list fpOutward
                        if (!Utility.ContainsList(faceVerticeIndex, verticeIndexInOneFace))
                        {
                            faceVerticeIndex.add(verticeIndexInOneFace);

                            if (onRightCount == 0)
                            {
                                fpOutward.add(trianglePlane);
                            }
                            else if (onLeftCount == 0)
                            {
                                fpOutward.add(GeoPlane.Negative(trianglePlane));
                            }
                        }
                    }
                    else
                    {                    
                        // possible reasons:
                        // 1. the plane is not a face of a convex 3d polygon,
                        //    it is a plane crossing the convex 3d polygon.
                        // 2. the plane is a face of a concave 3d polygon
                    }

                } // k loop
            } // j loop        
        } // i loop                        

        // return number of faces
        numberOfFaces = faceVerticeIndex.size();
        
        for (int i = 0; i < numberOfFaces; i++)
        {                
            // return face planes
            facePlanes.add(new GeoPlane(fpOutward.get(i).getA(), fpOutward.get(i).getB(), 
                                        fpOutward.get(i).getC(), fpOutward.get(i).getD()));

            ArrayList<GeoPoint> gp = new ArrayList<GeoPoint>();

            ArrayList<Integer> vi = new ArrayList<Integer>();
            
            int count = faceVerticeIndex.get(i).size();
            for (int j = 0; j < count; j++)
            {
                vi.add(faceVerticeIndex.get(i).get(j));
                gp.add( new GeoPoint(vertices.get(vi.get(j)).getX(),
                                     vertices.get(vi.get(j)).getY(),
                                     vertices.get(vi.get(j)).getZ()));
            }

            // return faces
            faces.add(new GeoFace(gp, vi));
        }
        
        this.Faces = faces;
        this.FacePlanes = facePlanes;
        this.NumberOfFaces = numberOfFaces;        
    }
}

Points of Interest

Unlike C# or C++, Java has no operator overloading and class member property. Java also doesn't support returning values from function arguments.

Since Java is cross-platform network language, it uses Big-Endian for its byte order, this introduces some additional works to process the LAS file, which is Little-Endian data.

History

• January 12, 2016: Initial version