4D Object By Change - Customizing the algorithm

The 4D-OBC implementation provided by py4dgeo can be customized similarly to how M3C2 is customized. The fundamental idea is to create your own algorithm class derived from py4dgeo’s RegionGrowingAlgorithm class and override only those parts of the algorithm that you want to change. This notebook introduces the currently existing customization points.

Implemented by

Dominic Kempf (@dokempf, Heidelberg University)

Author(s) of the method Katharina Anders, Lukas Winiwarter, Bernhard Höfle (Heidelberg University)

Original publication of the method Anders, K., Winiwarter, L., Mara, H., Lindenbergh, R., Vos, S. E., & Höfle, B. (2021). Fully automatic spatiotemporal segmentation of 3D LiDAR time series for the extraction of natural surface changes. ISPRS Journal of Photogrammetry and Remote Sensing, 173, pp. 297-308. doi: 10.1016/j.isprsjprs.2021.01.015.

import py4dgeo
import _py4dgeo  # The C++ bindings module for py4dgeo

import numpy as np
import tempfile
import shutil

from pathlib import Path
from py4dgeo.util import find_file

Numba is not installed.
Use `pip install numba` to install it to enable faster computations in Vapc.
You can run py4dgeo.Vapc without Numba, but it may be slower.

Now we download data from external zenodo link

# Fetch original test data
test_filename = "synthetic.zip"
original_file = find_file(test_filename)
print(f"File downloaded to: {original_file}")

# Create temporary directory
tmpdir = tempfile.mkdtemp()
zip_file = Path(tmpdir) / test_filename

# Copy original file into temporary directory
shutil.copy(find_file(test_filename), zip_file)
print(f"Working copy: {zip_file}")

File downloaded to: /home/docs/.cache/py4dgeo/synthetic.zip
Working copy: /tmp/tmpqyde7zw4/synthetic.zip

# Load test data
analysis = py4dgeo.SpatiotemporalAnalysis(zip_file)

# We are disabling log messages on this tutorial to increase the readability of the output
import logging

logging.disable()

Distance Measure

By default, 4D-OBC uses a normalized Dynamic Time Warping distance measure. You can provide your own Python function, although all the same warnings as with M3C2 apply: Python code will be significantly slower compared to C++ implementations and will be run sequentially even if you are using OpenMP.

def custom_distance(params: _py4dgeo.TimeseriesDistanceFunctionData):
    mask = ~np.isnan(params.ts1) & ~np.isnan(params.ts2)
    if not np.any(mask):
        return np.nan

    # Mask the two input arrays
    masked_ts1 = params.ts1[mask]
    masked_ts2 = params.ts2[mask]

    return np.sum(np.abs(masked_ts1 - masked_ts2)) / np.sum(
        np.abs(np.maximum(masked_ts1, masked_ts2))
    )

class CustomDistance4DOBC(py4dgeo.RegionGrowingAlgorithm):
    def distance_measure(self):
        return custom_distance

analysis.invalidate_results()
objects = CustomDistance4DOBC(
    neighborhood_radius=2.0,
    seed_subsampling=20,
).run(analysis)

The params data structure passed to the distance function contains the following fields:

• ts1 and ts2 are the time series to compare which may include NaN values

• norm1 and norm2 which are normalization factors

Prioritizing seeds

The 4D-OBC algorithm finds a number of seed locations for region growing and then prioritizes these seeds by sorting them according to a criterion. You can pass your own criterium like this:

def sorting_criterion(seed):
    # Choose a random score, resulting in random seed order
    return np.random.rand()

class CustomSeedSorting(py4dgeo.RegionGrowingAlgorithm):
    def seed_sorting_scorefunction(self):
        return sorting_criterion

analysis.invalidate_results()
objects = CustomSeedSorting(
    neighborhood_radius=2.0,
    seed_subsampling=20,
).run(analysis)

Rejecting grown objects

After the region growing is done, the algorithm class calls it method filter_objects to check whether the object should be used or discarded. The method must return True (keep) or False (discard):

class RejectSmallObjects(py4dgeo.RegionGrowingAlgorithm):
    def filter_objects(self, obj):
        return len(obj.indices) > 10

analysis.invalidate_results()
objects = RejectSmallObjects(
    neighborhood_radius=2.0,
    seed_subsampling=20,
).run(analysis)

Seed point detection

If you would like to run an entirely different algorithm to determine the seeds for region growing, you can do so by overriding find_seedpoints:

from py4dgeo.segmentation import RegionGrowingSeed

class DifferentSeeds(py4dgeo.RegionGrowingAlgorithm):
    def find_seedpoints(self):
        # Use one seed for corepoint 0 and the entire time interval
        return [RegionGrowingSeed(0, 0, self.analysis.distances.shape[1] - 1)]

analysis.invalidate_results()
objects = DifferentSeeds(
    neighborhood_radius=2.0,
    seed_subsampling=20,
).run(analysis)