## 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: Algorithms, Boundaries, EKL, Geometry & Topology
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:

Thank you!

Thanks Renaud! That’s looks more elegant ?