DM_spatial_processing.py

## @package python.demo.DM_spatial_processing
# "
# example script showing how to spatially process an ODM.
# based on a generic neighborhood definition statistical features
# are extracted and attached as attributes to each point
# (script works in python 2 and 3)
# 
 
from __future__ import print_function  # print function syntax as in python 3
import sys
import  datetime
 
from opals import pyDM
 
 
##
# callback object computing the user-defined attributes '_pcount' and '_zmean'
class Kernel(pyDM.KernelPointEx):
 
    def __init__(self):
        # initialize the base class (mandatory!)
        if sys.version_info >= (3, 0):
            super().__init__()
        else:
            super(Kernel, self).__init__()
        return
 
    ##
    # callback notifies about a changed leaf
    def leafChanged(self, leaf, localTree):
        print("process tile with id = %d" % leaf.id() )
        return
 
    ##
    # callback notifies that all leaf handels should be released
    def releaseLeaf(self):
        return
 
    ##
    # callback for processing a point
    def process(self, pt, neighbours):
 
        # add neighbour count (_pcount) and z mean of neighbours (_zmean) to the current point
        zsum = 0
        for n in neighbours.points():
            zsum += n.z
 
        # numpyDict = pyDM.NumpyConverter.get(neighbours)  # get neighbours as numpy dict
 
        pt.info().set(0, neighbours.sizePoint())
        if neighbours.sizePoint() > 0:
            zmean = zsum / float(neighbours.sizePoint())
            pt.info().set(1, zmean)
 
        # return True if the point was changed otherwise False
        # this internally marks the current point leaf as changed (=needs to be written to disk)
        return True
 
 
def DM_spatial_processing(filename, maxOutput=1000):
    # pyDM.Datamanager.load parameter: filename(string), readOnly(bool) threadSafety(bool)
    dm = pyDM.Datamanager.load(filename, False, False)
    if not dm:
        print("Unable to open ODM '" + filename + "'")
        sys.exit(1)
 
    # creates the appropriate attribute layout for processing
    lf = pyDM.AddInfoLayoutFactory()
    lf.addColumn(pyDM.ColumnType.int32, "_pcount")  # column 0
    lf.addColumn(pyDM.ColumnType.float_, "_zmean")  # column 1
    layout = lf.getLayout()
    kernelLayout = lf.getLayout()
 
 
    # create spatial selection
    query = pyDM.QueryDescriptor("sphere(r=1)")  # use a sphere with radius 1m as neighbourhood definition (textual description)
    
    # sphere = pyDM.QuerySphere(r=1)        # sphere with radius 1m (class based description)
    # knn = pyDM.QueryKnn(k=12,dim=3)       # 3d knn with 12 neighbours (class based description)
    # query = pyDM.QueryDescriptor(sphere)           # create actual neighbourhood definition: fixed distance
    # query = pyDM.QueryDescriptor(knn & sphere)     # create actual neighbourhood definition: knn with a maximum search radius
 
    print("Ready to append attribute '_pcount' and '_zmean' to '%s'" % filename)
 
    print("Number of tiles to process = %d" % dm.getPointIndex().sizeLeaf())
 
    # create kernel object
    k = Kernel()
 
    start = datetime.datetime.now()
    # create processor
    # pyDM.ProcessorEx parameter:  dm(datamnager), query(QueryDescriptor),
    #                              processFilter, processLayout, processLayoutReadOnly,
    #                              neighbourFilter, neighbourLayout, neighbourLayoutReadOnly
    processor = pyDM.ProcessorEx(dm, query, None, layout, False, None, None, False)
    processor.run(k)  # perform computation
 
    diff = datetime.datetime.now() - start
    print("Processing took %.2f [s]" % diff.total_seconds() )
 
    # save manager object
    print("Save manager")
    dm.save()
 
    print("done")
 
if len(sys.argv) == 1:
    print("ODM parameter missing")
    sys.exit(-1)
 
DM_spatial_processing(sys.argv[1])