8

What is the most efficient way to find the minimal wall-thickness (value and location) of a complex, non convex polygon area including holes? See example of a polygon in blue, with the minimum wall-thickness in red, although in this case the location is ambiguous, if the two adjacent lines are parallel.

enter image description here

So far, we have tried:

  • Subdividing polygon lines and finding minimum point-point lines inside polygon (brute force, not efficient for complex polygons with >10'000 points)
  • Delaunay triangulation and finding minimal edges inside polygon. Not precise enough, only feasible if combined with subdivision of polygon lines first. Here's example (Nr 3), where the Delaunay triangulation would not find simplex edges in red but would miss out on the minimal wall-thickness in the green box:
    enter image description here
  • Iteratively increasing erosion buffer to find the minimal inset, where the erosion polygon breaks into multiple parts = half of minimal wall-thickness. Problem is to find the location of the minimum wall-thickness with this approach afterwards. Furthermore, erosion doesn't always break into multiple parts and misses "dead ends". Here's an example (Nr 2) that erodes to a line and does give the wrong minimal wall-thickness:
    enter image description here
  • Finding medial axis first, then search for minimum circle on the medial axis that is covering but not overlapping the polygon area. Edit: Problematic are the many "wrong candidates" on the medial axis: Eg. (Nr 1) circle A would be wrong, circle B would denote the correct minimum wall-thickness:
    enter image description here
  • Get the distance between all the line pairs to find the closest ones. – bugmenot123 Feb 5 '18 at 11:46
  • So what was wrong with the medial axis approach? – Hornbydd Feb 5 '18 at 11:47
  • 1
    @Hornbydd: Problem was, that there are many circles on the medial axis that are touching corners but are not defining wall-thickness. See second example: circle A would be wrong, circle B would be the correct location of the minimal wall-thickness. So medial axis looks like a compuational costly detour and not providing the right answer... – Oliver Staubli Feb 5 '18 at 13:51
  • 1
    If you do erosion until the polygon degenerates to two polygons touching at a point, then the location will be where a circle of radius the same as the buffer centred at the point touches the original polygon. That's a hypothesis presented without proof but I can't see a counterexample... – Spacedman Feb 5 '18 at 14:16
  • 1
    @OliverStaubli My suggestion is to check not only edges of delaunay triangles, but also heights of those triangles that have one edge on the boundary and the other two in the interior of the polygon. In example Nr.3 the height of the triangle under the green square is what you're looking for. (depending on the constraints of the triangulation you might need to also filter out some of the candidates in obtuse triangles) – mkadunc Feb 6 '18 at 17:14
2

Two more ideas to throw in the pot:

  1. Rasterize your polygon and use a distance transform (returns image of the shortest distance from each nonzero pixel to a zero-pixel). Skeletonize your rasterized image, then take the values of the distance transformed image along the skeleton. From that set you will have some minimums that should correspond to your narrowest width. This set can be used as initial search points to then implement your brute force approach. I should note that the skeleton will bisect the corners of objects, and at those locations the distance transform will approach zero (as you get closer to the object's boundary). This may be problematic, but represents an issue with your problem--why can't the smallest width be at a corner (and be essentially zero)? You could address this problem by setting a threshold on the shortest distance around the perimeter between the two points (if they lie on the same object). You can use a geodesic distance transform on the set of perimeter pixels to quickly find that value.

    This method requires you to make a decision about the resolution of the rasterized polygon, which introduces some scale-dependency. And if you choose too high a resolution, the distance transform can be time-consuming. But generally they're pretty fast. This method might not give you the precision you want, but it could at least give you a much smaller set of locations that needed to be checked.

  2. Your brute force method isn't a bad place to start. I had a similar problem where I had to find all intersections of a (long) line with itself, and I was able to vastly speed up search time using a kd-tree search algorithm (I used rangesearch in Matlab at the time) to find points within a neighborhood first. That way you're only brute-forcing a small subset of the total number of points.

  • Thanks @jon. Both promising approaches. I was considering kdTree already but hoped the problem described has a well known "copy-paste" solution :-) Will have to dig deeper... – Oliver Staubli Feb 5 '18 at 16:27
1

One of most efficient methods to find the minimal wall-thickness (value and location) of a complex, non convex polygon area including holes, it could be by using a regularly spaced layer (or random) of points for determining, first, closest segment with context for each point and, next, the point for intersection between incremental segment and opposite side polygon; based in directors cosines.

Incremental distances can be used until first segment reaches and intersects some side polygon (the minimal wall-thickness).

To try out my approach I cloned your polygon with holes and created a random points layer inside polygon with 100 points; as it can be observed at following image:

enter image description here

PyQGIS used code looks as follows:

import math

def azimuth(point1, point2):
   return point1.azimuth(point2) #in degrees

def cosdir_azim(azim):
   azim = math.radians(azim)
   cosa = math.sin(azim)
   cosb = math.cos(azim)
   return cosa,cosb

registry = QgsMapLayerRegistry.instance()    
polygon = registry.mapLayersByName('polygon_with_holes')
point_layer = registry.mapLayersByName('Random_points')
points = [ feat.geometry().asPoint() for feat in point_layer[0].getFeatures() ]
feat_polygon = polygon[0].getFeatures().next()

#producing rings polygons
rings_polygon = feat_polygon.geometry().asPolygon()
segments = []
epsg = point_layer[0].crs().authid()
uri = "LineString?crs=" + epsg + "&field=id:integer""&index=yes"
mem_layer = QgsVectorLayer(uri,
                           'increasing_segments',
                           'memory')
prov = mem_layer.dataProvider()
length = 10
pt2 = 0 
k = 0
while pt2 == 0:
    for i, point in enumerate(points):
        #determining closest distance to vertex or side polygon
        dist1 = feat_polygon.geometry().closestSegmentWithContext(point)[0]
        #determining point with closest distance to vertex or side polygon
        pt = feat_polygon.geometry().closestSegmentWithContext(point)[1]
        cosa, cosb = cosdir_azim(azimuth(pt, point))
        #extending segment in opposite direction based in director cosine and length 
        op_pt  = QgsPoint(point.x() + (length*cosa), point.y() + (length*cosb))
        segments.append([pt,op_pt])
        geom = QgsGeometry.fromPolyline([point,op_pt])
        for ring in rings_polygon:
            geom_ring = QgsGeometry.fromPolyline(ring)
            if geom.intersects(geom_ring):
                pt3 = geom.intersection(geom_ring)
                pt2 = pt3.distance(QgsGeometry.fromPoint(point))
                ms = [pt3.asPoint(), pt]
    length += 100
    k += 1
new_segments = segments[len(segments) -1 - len(segments)/k: len(segments) - 1]

feats = [ QgsFeature() for i in range(len(new_segments)) ]
for i,feat in enumerate(feats):
    feat.setAttributes([i])
    geom = QgsGeometry.fromPolyline(new_segments[i])
    feat.setGeometry(geom)

prov.addFeatures(feats)
QgsMapLayerRegistry.instance().addMapLayer(mem_layer)
minimum_segment = QgsGeometry().fromPolyline(ms).exportToWkt()
print minimum_segment, k

and it produces a memory layer of incremental distances (only for visualization purposes) and print minimal wall-thickness in WKT format.

After running the code at Python Console of QGIS I got result of following image:

enter image description here

It can be observed that only one incremental distance reached the opposite side first in expected area.

Printed WKT format (for minimal wall-thickness) is used with QuickWKT plugin of QGIS for visualizing that segment at following image:

enter image description here

The slight inclination was produced because "closest segment with context" was joined to a vertex; instead side polygon. However, it can be avoided with a code exception or more points.

0

The medial axis approach is correct, you just need a criterion to ignore bad circles: Each circle at the medial axis touches the surface in two (or more) points. Imagine vectors from the circle's center to these points on the surface and observe that the angle between is 180° for the good circle B and only 90° for the bad circle A.

enter image description here

0

A generic polygon drawing algorithm works by sorting line segments top to bottom, and working through them to draw orthogonal lines of pixels. As this way of breaking down polygons (without curvatures) is very quick, it can be used as a basis. Instead of only going top-to-bottom, you can go 0, 30, 60, and 90 degrees, and find your shortest orthogonal section (=minimal wall-thickness!), you only have to calculate that once per point and not for any kind of 'pixel resolution'.

See computer-graphics-scan-line-polygon-fill-algorithm

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.