Click here to Skip to main content
14,971,209 members
Articles / Programming Languages / Java / Java SE
Article
Posted 12 Jan 2016

Stats

15.9K views
291 downloads
3 bookmarked

Point Inside 3D Convex Polygon in Java

Rate me:
Please Sign up or sign in to vote.
5.00/5 (2 votes)
12 Jan 2016CPOL
An algorithm to determine if a point is inside a 3D convex polygon for a given polygon vertices in Java

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:

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

Then, construct the main process instance:

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

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

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

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

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

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

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

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

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

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.

Comments and Discussions

 
Praisevote of 5 Pin
Beginner Luck25-Jan-16 18:43
professionalBeginner Luck25-Jan-16 18:43 

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.