14 Aug, 2019

Bounding Geometry With EKL: 2D Convex Hull

We often need to build simplified representations of geometry to facilitate solving problems related to packing and stacking; area and volume (weight) estimation; and more generally, space reservation. One such simplified representation is finding a 2D convex hull of a point set. Let’s take a look.

Pseudo Code

A convex hull is an unambiguous boundary of edges that encompasses a set of points. The speed of generating it is very driven by the algorithm and the point-set. There are many algorithms to get the hull. The one I implemented here is quite simple. I look for edges on the boundary, which I identify as follows:

    - Pick a starting point, then pick any second point. Make an edge.
    - Test if all other points are to the left (or right) of this edge. If so, store the two points. 
    - Repeat for the next two points in the list. 
    - Once hull points are found, remove collinear ones to build a polygon with minimum number of edges (this is not part of the definition of the hull).

To check of all other points are to the left (or right) of an edge, I computing an area of a triangle made by the first two points, and a third point. If positive (or negative), I increment a variable (for example triangle-count) by 1. At the end of the loop, If the list of triangle count of a direction( positive or negative) is equal to the overall point count -2, then the edge defined by the first two points is a boundary edge.

Here are a couple point-set examples:


Code

Here is the code. Note that I first clean the point set to remove any overlapping points. I also only take points directly under the input geoset. If the intent is to get points nested inside, then use the Query method.


/* Action created by melkhaldi 8/13/2019 */
//INPUTS: PointGeoSet: OpenBodyFeature: Geoset that includes points.

//here, we're taking the points directly inside the geometrical set.
//if the intent is to get all points nested (for example inside other geosets), then use 
//PointGeoSet.Query("Point", "")
let  pointCloud(List)
pointCloud = PointGeoSet.Children.Filter("Point", "")

//let's make sure there are no points with the same coordinates--if found, only take one.

let testPoint, otherPoint(Point)
let t=1

for t while t<= pointCloud.Size(){
    testPoint =pointCloud.GetItem(t)

    let tt=1
    for tt while tt<= pointCloud.Size() {
        if(t==tt) continue
        otherPoint =pointCloud.GetItem(tt)
        if(distance( otherPoint, testPoint)==0mm){
            pointCloud.RemoveItem(t)
            break
        }
    }
}

//let's find the points that are on the boundary
let hullPoints(List)
let A, B, C(Point)
let triangleArea(Area)
let triangles(Integer)
let Ax, Ay, Bx, By, Cx, Cy(LENGTH)
for A inside pointCloud{
    //for every point in B
    for B inside pointCloud{

        //reset triangle count
        triangles=0

        //for every point in C
        for C inside pointCloud {
            //skip if this any of the distances between ABC is zero
          //--to avoid created triangles with zero-length sizes.
            if distance(A, B)==0mm or distance(B, C)==0mm or distance(A,C)==0mm continue
            //compute triangle
            Ax = A.coord(1)
            Ay = A.coord(2)
            Bx  = B.coord(1)
            By  = B.coord(2)
            Cx  = C.coord(1)
            Cy  = C.coord(2)
            triangleArea = (Ax*( By-Cy) + Bx*(Cy-Ay) + Cx*(Ay-By))/2
            //if value is negative (or positive)
          //either one works, but stay consistent
          //zero area is included to account 
          //for points that are collinear.
            if triangleArea<=0m2 {
                //increment triangle count by 1
                triangles = triangles+1
            }
        }

        if triangles ==pointCloud.Size()-2{ 
            //if index is zero, then it is not in the list. add it.
            if ( hullPoints.IndexOf(A, 1)==0 )
                hullPoints.Append(A)
            if ( hullPoints.IndexOf(B, 1)==0 )
                hullPoints.Append(B)
        }

    }
}

//build reference information
let avgPoint(Point)
avgPoint = centerofgravity( hullPoints)

//define any reference segment using the center point (average point) and any other point.
let refAngleStartDir(Direction)
refAngleStartDir =direction( line(avgPoint, hullPoints.GetItem(1) ))    

//define a vertical direction to help identify negative and positive angle directions
//for that, we'll get a mean plane that passes through the point set
let refAngleNormalDir(Direction)
refAngleNormalDir = direction(  planemean(hullPoints))

//measure angles
let ang(Angle)
let pointAnglePairs(List)
let angleEndLine(Line)
//measure angles and make a list of point-angle pairs
let somePoint(Point)
for somePoint inside hullPoints {
    angleEndLine =line(avgPoint, somePoint)
    ang = angleoriented(refAngleStartDir , direction( angleEndLine), refAngleNormalDir)
    pointAnglePairs.Append(  List( somePoint, ang) )
}

//now sort the list of points by angle, and extract the points from the list.
let sortedHullPoints(List)
let sortedPaird(List)
//sort the pairs
sortedPaird =pointAnglePairs.Sort("<", "List", "Angle", "y=x.GetItem(2)")

//get the points for the pairs list
sortedHullPoints = sortedPaird.Extract("List", "Point", "y=x.GetItem(1)")

//now that we have points sorted radially, we can start making line segments 
//making sure that they will only connect head to tail
//so when a polygon is made, it won't self intersect.
//but first, lets remove any points that fall between other two points. 
//this is because we might have picked multiple collinear points.
//so we want to minimize them in order have fewer polygon edges

//lets check on end of hull points
A = sortedHullPoints.GetItem(sortedHullPoints.Size()-1)
B = sortedHullPoints.GetItem(sortedHullPoints.Size())
C = sortedHullPoints.GetItem(1)
//if this is true, the B falls between A and C --that is, ABC are collinear
if( distance(A, B) + distance( B,C) == distance(A,C)){
    sortedHullPoints.RemoveItem(sortedHullPoints.Size())
}
//let's loop
let m=2
for m while m< sortedHullPoints.Size(){
    A = sortedHullPoints.GetItem(m-1)
    B = sortedHullPoints.GetItem(m)
    C = sortedHullPoints.GetItem(m+1)
    if( distance(A, B) + distance( B,C) == distance(A,C)){
        sortedHullPoints.RemoveItem(m)
        m=2
    }
} 
//capture the start of hull points
A = sortedHullPoints.GetItem(sortedHullPoints.Size())
B = sortedHullPoints.GetItem(1)
C = sortedHullPoints.GetItem(2)
if( distance(A, B) + distance( B,C) == distance(A,C)){
    sortedHullPoints.RemoveItem(1)
}

//last , lets check on point count
if sortedHullPoints.Size()<3{
    Notify("Point count is less than 3")
    exit
}

//make a polygon
let aLine(Line)
let lines(List)
let i=1
//loop through the points and make lines
//notice we are stopping before the end by one
//before we are using point at i and i+1
for i while i< sortedHullPoints.Size() {
    aLine=line( sortedHullPoints.GetItem(i), sortedHullPoints.GetItem(i+1 ))
    lines.Append( aLine )
}
//build the last line
aLine = line( sortedHullPoints.GetItem( sortedHullPoints.Size()), sortedHullPoints.GetItem(1))
lines.Append(aLine)

//build the hull
let hull(Curve)
//add it to the tree
hull = new("Curve", "HULL",PointGeoSet: OpenBodyFeature)
//give it the definition as the result of an assemble function
hull = assemble(lines): Curve

//update the geoset
PointGeoSet.Update()

Handling Polygons

We can extend the above code to extract points from a polygon. For example:

let PointGeoSet(OpenBodyFeature)
set PointGeoSet = poly.Owner 

let  pointCloud(List)
pointCloud =  poly.GetSubElements(1, false).Extract("CATEdge", "Curve", "set y=extract(x)").Extract("Curve", "Point", "y= pointoncurveRatio(x,NULL,0,true)")

Take a look:


Since CATIA will only allow creating correct geometry, we can omit the step of removing overlapping points since we are extracting points from a curve (polygon).

oh, I’ll let you explain that one line of code to get points from a polygon 😉

Tags: , , ,

About : Maher Elkhaldi

Maher Elkhaldi is a senior applications engineer at Tesla Motors. He founded the 3DXAutomation blog to help make knowledge of programming CATIA easier to find, and contribute to the open-source community.

2 thoughts on “Bounding Geometry With EKL: 2D Convex Hull”

  1. Very interesting post! Huge fan of this blog.

    When I remove items from a list I always go backward in order not to skip an item when the counter gets incremented over a removed item. Also, I limit the nested loop to the items not scanned yet, this avoid duplicate tests and also removes the if t == tt condition.

    Here’s an optimized version of your code that filters out duplicates points:


    let t (Integer) t = pointCloud.Size() for t while t>0{ testPoint = pointCloud.GetItem(t) let tt = 1 for tt while tt<t{ otherPoint = pointCloud.GetItem(tt) if(distance(otherPoint, testPoint)==0mm){ pointCloud.RemoveItem(t) break } } t = t-1 }

    Thank you!

Leave a Reply

This site uses Akismet to reduce spam. Learn how your comment data is processed.