demoSegmentation.py

## @package python.demoSegmentation
## @brief unkown python script prefix=autobuild.swdvlp64.opals.distro.demo
# 
#   example script showing how to use opalsSegmentation
#   and access its results based on the segment manager.
#   If the parameter alphaRadius is used, segment alpha shapes can be accessed and visualized.
#   For both segmentation methods (conditional clustering and plane extraction), example parameters are set.
# 

from __future__ import print_function  # print function syntax as in python 3

import matplotlib
import scipy
import numpy as np

import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

from opals.tools.flags import Debug_Mode

if Debug_Mode == False:
    from opals import pyDM, Segmentation, Import, Normals  # for pyDM in release mode
else:
    from opals import pyDM_d as pyDM  # for pyDM in debug mode
    from opals import Segmentation_d as Segmentation
    from opals import Import_d as Import
    from opals import Normals_d as Normals

import opals

def set_axes_radius(ax, origin, radius):
    ax.set_xlim3d([origin[0] - radius, origin[0] + radius])
    ax.set_ylim3d([origin[1] - radius, origin[1] + radius])
    ax.set_zlim3d([origin[2] - radius, origin[2] + radius])

##
# Make axes of 3D plot have equal scale so that spheres appear as spheres,
#     cubes as cubes, etc..  This is one possible solution to Matplotlib's
#     ax.set_aspect('equal') and ax.axis('equal') not working for 3D.
# 
#     Input
#       ax: a matplotlib axis, e.g., as output from plt.gca().
#     
def set_axes_equal(ax):

    limits = np.array([
        ax.get_xlim3d(),
        ax.get_ylim3d(),
        ax.get_zlim3d(),
    ])

    origin = np.mean(limits, axis=1)
    radius = 0.5 * np.max(np.abs(limits[:, 1] - limits[:, 0]))
    set_axes_radius(ax, origin, radius)


##
# output all attribute values within an addinfo object
def output_addinfo(info):
    print(" attributes: ",end="")
    if info.columns() == 0:
        print("NONE")
    else:
        for i in range(info.columns()):
            if i > 0:
                print(", ", end="")

            print(info.name(i)+"=", end="")
            if info.isNull(i):
                print("NULL", end="")
            else:
                print(info.get(i), end="")

        print() #line feed


##
# Run necessary preprocessing steps for the point cloud
#         e.g. Import point cloud into OPALS Datamanager (odm) with opalsImport,
#              Compute normal vectors with opalsNormals
#     
def preprocessing_steps(name):
    imp = Import.Import(inFile=name+".laz",tileSize=100,filter="generic[z>265]")
    imp.run()

    normals = Normals.Normals(inFile=name+".odm", normalsAlg=opals.Types.NormalsAlgorithm.simplePlane, neighbours=8, searchMode=opals.Types.SearchMode.d2, searchRadius=1)
    normals.run()


# ATTENTION: programm needs to be excecuted within the opals demo directory or
#            copy strip21.laz in your current directory

# necessary preprocessing steps
filename_stem = "strip21"
preprocessing_steps(filename_stem)

# Select SegmentationMethod here:
# mode = opals.Types.SegmentationMethod.planeExtraction
# or
mode = opals.Types.SegmentationMethod.condClustering

# perform segmentation using plane extraction mode and alpha shape computation
mod = Segmentation.Segmentation()

filename = filename_stem+".odm"
if mode == opals.Types.SegmentationMethod.planeExtraction:
    mod.inFile = filename
    mod.searchRadius = 1
    mod.minSegSize = 50
    mod.method = opals.Types.SegmentationMethod.planeExtraction
    mod.planeExtraction.maxDist = 0.2
    mod.planeExtraction.maxSigma = 0.15
    mod.planeExtraction.seedCalculator = "NormalSigma0<0.02 AND Z > 275 ? NormalSigma0 : invalid"
    mod.filter = ["Region[529568.8 5338773.5 529863 5338773.9 529862.6 5338642.4 529706 5338685.4 529569.1 5338753.8]"]
    mod.alphaRadius = 1
    mod.sort = opals.Types.SegmentSort.descendingPSize
elif mode == opals.Types.SegmentationMethod.condClustering:
    mod.inFile = filename
    mod.searchRadius = 1
    mod.minSegSize = 50
    mod.method = opals.Types.SegmentationMethod.condClustering
    mod.condClustering.criterion = "normalX*n[0].normalX+normalY*n[0].normalY+normalZ*n[0].normalZ>0.996"
    # mod.condClustering.alphaShapeRefPlane = [0, 0, 1]
    mod.alphaRadius = 1
    mod.sort = opals.Types.SegmentSort.descendingPSize
mod.run()

# get segmentation manager (gives access to the segment objects)
segManager = mod.segments

# pyDM.Datamanager.load parameter: filename(string), readOnly(bool) threadSafety(bool)
odm = pyDM.Datamanager.load(filename, True, False)
pi = odm.getPointIndex()  # get point index object of odm

print("The segmentation process has found "+str(segManager.sizeSegments())+" segment(s)")

max_pt_output = 5
ax = Axes3D(plt.figure())
ax.set_aspect('equal')
# iterate over all segments and output relevant segments properties
for seg in segManager.getSegments():
    # print basic segment information
    print("Segment-Nr %d:\t %d points; COG (%.3f / %.3f / %.3f)" % (seg.getID(), seg.sizePoints(), seg.getCoG().x, seg.getCoG().y, seg.getCoG().z ) )

    # get plane parameters
    try:
        planeParams = seg.getPlaneParams()
        assert len(planeParams) == 5
        print("\tplane parameters: a=%.4f b=%.4f c=%.4f d=%.4f (sigma0=%.3f)" % tuple(planeParams) )
    except:
        print("No plane params")

    # get alpha shape (can be empty in case the computation failed)
    alphaShape = seg.getAlphaShape()
    if alphaShape:
        print("\talpha shape consists of %d parts:" % alphaShape.sizePart())
        # get random color for segment alpha shape
        col = matplotlib.colors.rgb2hex(scipy.rand(3))
        for idx, part in enumerate(alphaShape.parts()):
            print("\t\t%d part has %3d points" % (idx,part.sizePoint()), end="")
            vtx_x = []
            vtx_y = []
            vtx_z = []
            for pt in part.points():
                vtx_x.append(pt.x)
                vtx_y.append(pt.y)
                vtx_z.append(pt.z)
            if len(vtx_x) < 3:
                print("Len < 3")
            # plot alpha shape
            tri = ax.plot3D(vtx_x, vtx_y, vtx_z, color=col)

            output_addinfo(part.info())

    # get segment point ids and print the first "max_pt_output" points
    print("\touput the first %d points: " % max_pt_output)
    ptIds = seg.getPoints()
    for idx, id in enumerate(ptIds):
        pt = pi.getPoint(id)
        print("\t\t%2d.point (%.3f / %.3f / %.3f)" % (idx, pt.x, pt.y, pt.z) )
        if idx+1 >= max_pt_output:
            break

set_axes_equal(ax)
plt.axis('off')

#plt.show()                   # use this command to display 3d graphic view
plt.savefig(filename_stem+"_alphaShapes.svg") # write figure as svg graph