Ground Extraction Using Differential Morphological Profiles

a.md

The algorithm is described in [1]. I've started a Jupyter Notebook that outlines the steps, but for those who just want to run with it, here is a pipeline that will apply the filter. Doing it this way, it's a bit slow (think 30 minutes for about 1.4M points on my MBP), but results look pretty good - perhaps good enough to justify coding up in C++ as a new stage.

Requires PDAL built with Python support, and Python packages astropy, numpy, pandas, and scipy. The parameters S, k, n, and b could be exposed as pdalargs to be user-defined. These values were appropriate for some of the ISPRS datasets.

[1] D. Mongus, N. Lukac, and B. Zalik, “Ground and building extraction from LiDAR data based on differential morphological profiles and locally fitted surfaces,” ISPRS J. Photogramm. Remote Sens., vol. 93, no. January, pp. 145–156, 2014.

dmp.json

{
    "pipeline":[
        {
            "type":"filters.python",
            "module":"anything",
            "function":"filter",
            "source":"from astropy.convolution import Gaussian2DKernel, convolve
from scipy import spatial, ndimage
from scipy.ndimage import morphology

import numpy as np
import pandas as pd

S = 30
k = 0.2
n = 0.3
b = 0.2

def idw(data):
    valid = np.argwhere(~np.isnan(data))
    nans = np.argwhere(np.isnan(data))
    dfg = pd.DataFrame(valid, columns=['X','Y'])
    dfng = pd.DataFrame(nans, columns=['X','Y'])
    tree = spatial.cKDTree(dfg)
    for _, row in dfng.iterrows():
        d, idx = tree.query([row.X, row.Y], k=3)
        d = np.power(d,-2)
        v0 = data[dfg.iloc[idx[0]][0], dfg.iloc[idx[0]][1]]
        v1 = data[dfg.iloc[idx[1]][0], dfg.iloc[idx[1]][1]]
        v2 = data[dfg.iloc[idx[2]][0], dfg.iloc[idx[2]][1]]
        data[row.X, row.Y] = (v0*d[0] + v1*d[1] + v2*d[2])/(d[0]+d[1]+d[2])
    return data

def filter(ins, outs):
    cls = ins['Classification']
    cls.fill(1)
    
    df3D = pd.DataFrame({'X': ins['X'], 'Y': ins['Y'], 'Z': ins['Z']})
    density = len(df3D) / (df3D.Y.ptp() * df3D.X.ptp())
    hres = 1. / density
    
    xi = np.ogrid[df3D.X.min():df3D.X.max():hres]
    yi = np.ogrid[df3D.Y.min():df3D.Y.max():hres]

    groups = df3D.groupby([np.digitize(df3D.X, xi), np.digitize(df3D.Y, yi)])

    cz = np.empty((yi.size, xi.size))
    cz.fill(np.nan)
    for name, group in groups:
        idx = group.Z.idxmin()
        cz[name[1]-1, name[0]-1] = df3D.Z.loc[idx]

    cz = idw(cz)

    struct = ndimage.iterate_structure(ndimage.generate_binary_structure(2, 1), 11).astype(int)
    opened = morphology.grey_opening(cz, structure=struct)

    struct = ndimage.iterate_structure(ndimage.generate_binary_structure(2, 1), 9).astype(int)
    closed = morphology.grey_closing(opened, structure=struct)

    lowx, lowy = np.where((closed - cz) >= 1.0)

    cz[lowx, lowy] = closed[lowx, lowy]
    stdev = 14
    G = Gaussian2DKernel(stdev)
    low = convolve(np.pad(cz,2*stdev,'edge'), G)[2*stdev:-2*stdev,2*stdev:-2*stdev]
    high = cz - low
    
    granulometry = []
    for scale in xrange(1, S+1):
        struct = ndimage.iterate_structure(ndimage.generate_binary_structure(2, 1), scale).astype(int)
        o = morphology.grey_opening(high, structure=struct)
        granulometry.append(o)
    
    out = []
    for i in xrange(1, len(granulometry)):
        out.append(granulometry[i-1]-granulometry[i])
    
    xs, ys = out[0].shape
    gprime = np.maximum.reduce(out)
    gstar = np.zeros((xs,ys))
    gplus = np.zeros((xs,ys))
    for ii in xrange(0,xs):
        for jj in xrange(0,ys):
            for kk in xrange(0,len(out)):
                if out[kk][ii,jj] < gprime[ii,jj]:
                    gplus[ii,jj] += out[kk][ii,jj]
                if out[kk][ii,jj] == gprime[ii,jj]:
                    gplus[ii,jj] += out[kk][ii,jj]
                    gstar[ii,jj] = kk
                    break
    
    T = k * gstar + n
    Sg = gprime < T
    F=cz.copy()
    F[np.where(Sg==0)] = np.nan
    G=idw(F)
    struct = ndimage.iterate_structure(ndimage.generate_binary_structure(2, 1), 1).astype(int)
    gradDTM = morphology.grey_dilation(G, structure=struct)
    xi = np.ogrid[df3D.X.min():df3D.X.max():hres]
    yi = np.ogrid[df3D.Y.min():df3D.Y.max():hres]
    ground = []
    nonground = []
    for i in range(0,len(df3D)):
        xbin = np.digitize(df3D.X[i], xi)
        ybin = np.digitize(df3D.Y[i], yi)
        if df3D.Z[i] < gradDTM[ybin-1,xbin-1]+b:
            ground.append(i)
        
    cls[ground] = 2

    outs['Classification'] = cls
    return True"
        }
    ]
}