Fitting circles to crescent shaped clusters in raster data automatically

Question

I have mobile lidar scanner (MLS) that was collected in a forest with big, mature trees.

My research involves fitting circles to tree stems to estimate their diameters at different heights.

Once I normalized the cloud and segmented out the trees, I sliced the point cloud into horizontal sections (like floors in a building), then converted them into geotiffs. The results look something like this in a good case:

As it stands, I'm using the 'draw circle from 3 points tool' in qgis shape digitizing toolbar to manually fit circles to each cluster. The result looks something like this:

Is there any intuitive way to automate this process? I find that doing it manually it is difficult to be consistent in addition to being extremely tedious. It would save me tons of time and effort if this could be avoided!

I think there may be a way to accomplish this via something like 'contour detection' in opencv, but I dont know how to go about integrating it with qgis.

Answer 1

Yes, there is a way to accomplish this via something like 'contour detection' in opencv, however, I used this tutorial for adapting following code in QGIS.

import cv2
import numpy as np

image=cv2.imread('circles.tif')
cv2.imshow('input image',image)

orig_image = image.copy()

gray=cv2.cvtColor(image,cv2.COLOR_BGR2GRAY)
ret, thresh=cv2.threshold(gray,127,255,cv2.THRESH_BINARY_INV)

_, contours, hierarchy=cv2.findContours(thresh.copy(),cv2.RETR_LIST,cv2.CHAIN_APPROX_NONE)

areas = []

for c in contours:
    area = cv2.contourArea(c)
    areas.append(area)

max_area = max(areas)

for c in contours:
    area = cv2.contourArea(c)
    x,y,w,h=cv2.boundingRect(c)
    cv2.rectangle(orig_image,(x,y),(x+w,y+h),(0,0,255),2)
    #calculate accuracy as a percent of contour perimeter
    accuracy=0.01*cv2.arcLength(c,True)
    approx=cv2.approxPolyDP(c,accuracy,True)
    if area == max_area:
        pol = c
        cv2.drawContours(image,[approx],0,(0,255,0),2)

cv2.imshow('Approx polyDP', image)

points = []
points2 = []

xmin = 523321.0347018922
ymax = 4795980.40210969 
xmax = 523322.2847018922 
ymin = 4795979.15210969

for point in pol:
    points2.append(QgsPointXY(xmin + 0.003*point[0][0], ymax - 0.003*point[0][1]))

polygon = [points2]

epsg = 32610

uri = "Polygon?crs=epsg:" + str(epsg) + "&field=id:integer""&index=yes"

mem_layer = QgsVectorLayer(uri,
                           'polygon',
                           'memory')

prov = mem_layer.dataProvider()

feats = [ QgsFeature() for i in range(len(polygon)) ]

for i, feat in enumerate(feats):
    feat.setAttributes([i])
    feat.setGeometry(QgsGeometry.fromPolygonXY([polygon[i]]))

prov.addFeatures(feats)

QgsProject.instance().addMapLayer(mem_layer)

Above code uses your image (I named it as circle.tif) that I projected with EPSG:32610 for an arbitrary point in Oregon and for producing a diameter about 0.7 m. It looks in QGIS as follows:

After running above code, it was produced polygon of following image:

After fixing the geometry and producing a Concave Hull with 20-nearest neighbors (this can be also automated in above script), I used a plugin developed by me for creating a circle in point (523321.79227764974, 4795979.818539615) with radius 0.35 (diameter 0.7 m). Unfortunately, my plugin uses center and radius for producing a circle but, it can be adapted for drawing a circle with 3 points and automatically indicates radius and diameter. These points are the key to fully automate complete process.

Editing Note:

Two points of circle can be obtained from intersection with oriented minimum bounding box and a third point from intersection of a perpendicular to this side with concave hull geometry.

Code for drawing circle with three points can be produced easily.