3

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:

points

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:

fit

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.

5
  • What you are looking for is called Hough transform. But I can't help you with open cv. – JRR Mar 27 '20 at 23:38
  • Does lidR have this capability? I tried using the package TreeLS (which you showed me) and it definitely uses Hough transform, but the diameters that result seem really unreliable. – Kartograaf Mar 28 '20 at 0:33
  • Not in lidR. And I don't know what TreeLS is capable of. But there are plenty of image processing software with hough tranform – JRR Mar 28 '20 at 0:37
  • side question. How did you export your lidar to a geotiff that look like that? i can export to DEM but nothing like you. – B-C B. Jan 12 at 20:20
  • 1
    First I used package lidR to normalize the cloud, then load in the cloud in 1m vertical increments such as slice <- readLAS(lasfile, select = "*" , filter = "-drop_z_below 1 -drop_z_above 2"). Then, convert the new cloud to a raster and save as tiff file. – Kartograaf Jan 12 at 20:24
3

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:

enter image description here

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

enter image description here

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.

enter image description here

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.

enter image description here

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

enter image description here

4
  • Wow @xunilk this is great! I will spend some time implementing it and if I can make it work I will award you the correct answer. Thank you! – Kartograaf Mar 31 '20 at 17:10
  • 1
    Thank you. This is the link to video youtube.com/watch?v=UtZLgann2ks where is showed fully automate complete process. Basis is the code and different tips about tools of Processing Tool Box posted in my answer. I'm open to any question. – xunilk Apr 1 '20 at 0:50
  • is this suppose to work if there is multiple circular shape. my output to the code is just a rectangle polygon with and hole in one of the positions representing a found shape (which is not actually at the right position). – B-C B. Feb 3 at 4:56
  • I've skipped the first step and created the polygon using polygonize tool. I then got the concave hulls and the oriented extent polygons. I am not sure what you mean in your editing notes about getting the 3 points to make the circle. how does it know which point to get from the interesection, there will be more than 2. and is the third point suppose to be the middle of the circle? What is the perpendicular to the hull? – B-C B. Feb 4 at 0:35

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.