Merge pull request #138 from eEcoLiDAR/development

Development

Author: cwmeijer

.gitignore

@@ -105,3 +105,5 @@ ENV/
 .idea/tasks.xml
 .idea/eEcoLiDAR.iml
 .idea/misc.xml
+.idea/dictionaries/
+.idea/inspectionProfiles/

.travis.yml

@@ -3,7 +3,6 @@ dist: trusty
 sudo: false
 
 python:
-- "2.7"
 - "3.5"
 - "3.6"
 

CHANGELOG.md

@@ -4,14 +4,23 @@ All notable changes to this project will be documented in this file.
 The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/)
 and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html).
 
-## [Unreleased]
+## 0.3.0 - 2019-02-27
 ### Added
+- Normalization module
 - General tests that all current and future feature extractors will be checked against.
+- Possibility to have a randomly subsampled (fixed) number of neighbors (eg for faster feature calculation)
+- Added feature extractors based on normalized height
 
 ## Changed
+- Many feature calculations are done in a vectorized way
+- z_entropy feature renamed to entropy_z for consistency
+- density_absolute_mean feature renamed to density_absolute_mean_z for consistency
+- perc_10 features renamed to perc_10_z for consistency
 
 ## Fixed
-- Fixed many feature extractors for corner cases (e.g. zero points)
+- Corner cases for many feature extractors (e.g. zero points)
+- Bug in slope feature calculation
+- Bug in fix density absolute mean feature calculation
 
 ## Removed
 

laserchicken/_version.py

@@ -1 +1 @@
-__version__ = '0.2.0'
+__version__ = '0.3.0'

laserchicken/compute_neighbors.py

@@ -80,10 +80,12 @@ def compute_sphere_neighborhood(environment_pc, target_pc, radius):
     neighborhoods = compute_cylinder_neighborhood(
         environment_pc, target_pc, radius)
 
+    counter = 0  # Target and neighborhood indices are going to be out of sync in loop below.
     for neighborhood_indices in neighborhoods:
         result = []
         for i, _ in enumerate(neighborhood_indices):
-            target_x, target_y, target_z = utils.get_point(target_pc, i)
+            target_x, target_y, target_z = utils.get_point(target_pc, counter)
+            counter += 1
             result_indices = []
             for j in neighborhood_indices[i]:
                 env_x, env_y, env_z = utils.get_point(environment_pc, j)
@@ -106,20 +108,22 @@ def compute_cell_neighborhood(environment_pc, target_pc, side_length):
     :return: indices of neighboring points from the environment point cloud for each target point
     """
 
-    max_radius = 0.5*math.sqrt((side_length ** 2) + (side_length ** 2))
+    max_radius = 0.5 * math.sqrt((side_length ** 2) + (side_length ** 2))
 
     neighbors = compute_cylinder_neighborhood(
         environment_pc, target_pc, max_radius)
 
+    counter = 0  # Target and neighborhood indices are going to be out of sync in loop below.
     for neighborhood_indices in neighbors:
         result = []
         for i, _ in enumerate(neighborhood_indices):
-            target_x, target_y, _ = utils.get_point(target_pc, i)
+            target_x, target_y, _ = utils.get_point(target_pc, counter)
+            counter += 1
             neighbor_indices = neighborhood_indices[i]
             result_indices = []
             for j in neighbor_indices:
                 env_x, env_y, _ = utils.get_point(environment_pc, j)
-                if ((abs(target_x - env_x)) > 0.5*side_length) or ((abs(target_y - env_y)) > 0.5*side_length):
+                if ((abs(target_x - env_x)) > 0.5 * side_length) or ((abs(target_y - env_y)) > 0.5 * side_length):
                     continue
                 else:
                     result_indices.append(j)
@@ -141,10 +145,12 @@ def compute_cube_neighborhood(environment_pc, target_pc, side_length):
     neighbors = compute_cell_neighborhood(
         environment_pc, target_pc, side_length)
 
+    counter = 0
     for neighborhood_indices in neighbors:
         result = []
         for i, _ in enumerate(neighborhood_indices):
-            _, _, target_z = utils.get_point(target_pc, i)
+            _, _, target_z = utils.get_point(target_pc, counter)
+            counter += 1
             neighbor_indices = neighborhood_indices[i]
             result_indices = []
             for j in neighbor_indices:
@@ -157,33 +163,41 @@ def compute_cube_neighborhood(environment_pc, target_pc, side_length):
         yield result
 
 
-def compute_neighborhoods(env_pc, target_pc, volume_description):
+def compute_neighborhoods(env_pc, target_pc, volume_description, sample_size=None):
     """
     Find a subset of points in a neighbourhood in the environment point cloud for each point in a target point cloud.
 
     :param env_pc: environment point cloud
     :param target_pc: point cloud that contains the points at which neighborhoods are to be calculated
     :param volume_description: volume object that describes the shape and size of the search volume
+    :param sample_size: maximum number of neighbors returned per target point; if None (default), all are returned
     :return: indices of neighboring points from the environment point cloud for each target point
     """
     volume_type = volume_description.get_type()
-    
+
     if volume_type == Cell.TYPE:
-        neighbors1 = compute_cell_neighborhood(env_pc, target_pc, volume_description.side_length)
-        for x in neighbors1:
-            yield x
+        neighbors = compute_cell_neighborhood(env_pc, target_pc, volume_description.side_length)
     elif volume_type == Cube.TYPE:
         neighbors = compute_cube_neighborhood(env_pc, target_pc, volume_description.side_length)
-        for x in neighbors:
-            yield x
     elif volume_type == Sphere.TYPE:
         neighbors = compute_sphere_neighborhood(env_pc, target_pc, volume_description.radius)
-        for x in neighbors:
-            yield x
     elif volume_type == InfiniteCylinder.TYPE:
         neighbors = compute_cylinder_neighborhood(env_pc, target_pc, volume_description.radius)
-        for x in neighbors:
-            yield x
     else:
         raise ValueError(
             'Neighborhood computation error because volume type "{}" is unknown.'.format(volume_type))
+
+    for x in neighbors:
+        yield _subsample_if_necessary(x, sample_size)
+
+
+def _subsample_if_necessary(x, sample_size):
+    if sample_size:
+        sampled = [_subsample(elements, sample_size) if len(elements) > sample_size else elements for elements in x]
+        return sampled
+    else:
+        return x
+
+
+def _subsample(neighborhood, sample_size=None):
+    return list(np.random.choice(neighborhood, sample_size, replace=False))

laserchicken/feature_extractor/__init__.py

@@ -10,28 +10,49 @@
 from laserchicken import keys, utils
 from .density_feature_extractor import PointDensityFeatureExtractor
 from .echo_ratio_feature_extractor import EchoRatioFeatureExtractor
-from .eigenvals_feature_extractor import EigenValueFeatureExtractor
-from .entropy_feature_extractor import EntropyFeatureExtractor
-from .height_statistics_feature_extractor import HeightStatisticsFeatureExtractor
-from .normal_plane_feature_extractor import NormalPlaneFeatureExtractor
-from .percentile_feature_extractor import PercentileFeatureExtractor
+from .eigenvals_feature_extractor import EigenValueVectorizeFeatureExtractor
+from .entropy_z_feature_extractor import EntropyZFeatureExtractor
+from .percentile_z_feature_extractor import PercentileZFeatureExtractor
 from .pulse_penetration_feature_extractor import PulsePenetrationFeatureExtractor
 from .sigma_z_feature_extractor import SigmaZFeatureExtractor
+from .median_z_feature_extractor import MedianZFeatureExtractor
+from .range_z_feature_extractor import RangeZFeatureExtractor
+from .var_z_feature_extractor import VarianceZFeatureExtractor
+from .mean_std_coeff_z_feature_extractor import MeanStdCoeffZFeatureExtractor
+from .skew_z_feature_extractor import SkewZFeatureExtractor
+from .kurtosis_z_feature_extractor import KurtosisZFeatureExtractor
+from .skew_norm_z_feature_extractor import SkewNormZFeatureExtractor
+from .mean_std_coeff_norm_z_feature_extractor import MeanStdCoeffNormZFeatureExtractor
+from .var_norm_z_feature_extractor import VarianceNormZFeatureExtractor
+from .range_norm_z_feature_extractor import RangeNormZFeatureExtractor
+from .kurtosis_norm_z_feature_extractor import KurtosisNormZFeatureExtractor
+from .entropy_norm_z_feature_extractor import EntropyNormZFeatureExtractor
+from .median_norm_z_feature_extractor import MedianNormZFeatureExtractor
+from .percentile_norm_z_feature_extractor import PercentileNormZFeatureExtractor
+from .density_absolute_mean_z_feature_extractor import DensityAbsoluteMeanZFeatureExtractor
+from .density_absolute_mean_norm_z_feature_extractor import DensityAbsoluteMeanNormZFeatureExtractor
+
+
+def _create_feature_map(module_name=__name__):
+    """Construct a mapping from feature names to feature extractor classes."""
+    name_extractor_pairs = _find_name_extractor_pairs(module_name)
+    return {feature_name: extractor for feature_name, extractor in name_extractor_pairs}
 
 
-def _feature_map(module_name=__name__):
-    """Construct a mapping from feature names to feature extractor classes."""
+def _find_name_extractor_pairs(module_name):
     module = importlib.import_module(module_name)
-    return {feature_name: extractor
-            for name, extractor in vars(module).items() if re.match('^[A-Z][a-zA-Z0-9_]*FeatureExtractor$', name)
-            for feature_name in extractor.provides()
-            }
+    name_extractor_pairs = [(feature_name, extractor)
+                            for name, extractor in vars(module).items() if
+                            re.match('^[A-Z][a-zA-Z0-9_]*FeatureExtractor$', name)
+                            for feature_name in extractor.provides()]
+    return name_extractor_pairs
 
 
-FEATURES = _feature_map()
+FEATURES = _create_feature_map()
 
 
-def compute_features(env_point_cloud, neighborhoods, target_idx_base, target_point_cloud, feature_names, volume, overwrite=False,
+def compute_features(env_point_cloud, neighborhoods, target_idx_base, target_point_cloud, feature_names, volume,
+                     overwrite=False,
                      verbose=True, **kwargs):
     """
     Compute features for each target and store result as point attributes in target point cloud.
@@ -60,18 +81,23 @@ def compute_features(env_point_cloud, neighborhoods, target_idx_base, target_poi
     :return: None, results are stored in attributes of the target point cloud
     """
     _verify_feature_names(feature_names)
-    ordered_features = _make_feature_list(feature_names)
+    wanted_feature_names = feature_names + [existing_feature for existing_feature in target_point_cloud[keys.point]]
+    extended_features = _make_extended_feature_list(feature_names)
+    features_to_do = extended_features
 
-    for feature in ordered_features:
-        if (target_idx_base == 0) and (not overwrite) and (feature in target_point_cloud[keys.point]):
+    while features_to_do:
+        feature_name = features_to_do[0]
+
+        if (target_idx_base == 0) and (not overwrite) and (feature_name in target_point_cloud[keys.point]):
             continue  # Skip feature calc if it is already there and we do not overwrite
 
+        extractor = FEATURES[feature_name]()
+
         if verbose:
-            sys.stdout.write('Feature "{}"\n'.format(feature))
+            sys.stdout.write('Feature(s) "{}"'.format(extractor.provides()))
             sys.stdout.flush()
             start = time.time()
 
-        extractor = FEATURES[feature]()
         _add_or_update_feature(env_point_cloud, neighborhoods, target_idx_base,
                                target_point_cloud, extractor, volume, overwrite, kwargs)
         utils.add_metadata(target_point_cloud, type(
@@ -82,6 +108,21 @@ def compute_features(env_point_cloud, neighborhoods, target_idx_base, target_poi
             sys.stdout.write(' took {:.2f} seconds\n'.format(elapsed))
             sys.stdout.flush()
 
+        for provided_feature in extractor.provides():
+            if provided_feature in features_to_do:
+                features_to_do.remove(provided_feature)
+
+    _keep_only_wanted_features(target_point_cloud, wanted_feature_names)
+
+
+def _keep_only_wanted_features(target_point_cloud, wanted_feature_names):
+    redundant_features = [f for f in target_point_cloud[keys.point] if f not in wanted_feature_names]
+    if redundant_features:
+        print('The following unrequested features were calculated as a side effect, but will not be returned:',
+              redundant_features)
+    for f in redundant_features:
+        target_point_cloud[keys.point].pop(f)
+
 
 def _verify_feature_names(feature_names):
     unknown_features = [f for f in feature_names if f not in FEATURES]
@@ -90,8 +131,8 @@ def _verify_feature_names(feature_names):
                          .format(', '.join(unknown_features), ', '.join(FEATURES.keys())))
 
 
-def _add_or_update_feature(env_point_cloud, neighborhoods, target_idx_base, target_point_cloud, extractor, volume, overwrite, kwargs):
-    #n_targets = len(target_point_cloud[keys.point]["x"]["data"])
+def _add_or_update_feature(env_point_cloud, neighborhoods, target_idx_base, target_point_cloud, extractor, volume,
+                           overwrite, kwargs):
     n_targets = len(neighborhoods)
 
     for k in kwargs:
@@ -101,35 +142,73 @@ def _add_or_update_feature(env_point_cloud, neighborhoods, target_idx_base, targ
     feature_values = [np.empty(n_targets, dtype=np.float64)
                       for i in range(n_features)]
 
-    print("Number of targets: %d, number of features: %d" % (n_targets, n_features))
+    if hasattr(extractor, 'is_vectorized'):
+        _add_or_update_feature_in_chunks(env_point_cloud, extractor, feature_values, n_features, n_targets,
+                                         neighborhoods,
+                                         target_idx_base, target_point_cloud, volume)
+    else:
+        _add_or_update_feature_one_by_one(env_point_cloud, extractor, feature_values, n_features, n_targets,
+                                          neighborhoods, target_idx_base, target_point_cloud, volume)
+
+    for i in range(n_features):
+        feature = provided_features[i]
+        if target_idx_base != 0:
+            if feature not in target_point_cloud[keys.point]:
+                continue
+
+            target_point_cloud[keys.point][feature]["data"] = np.hstack(
+                [target_point_cloud[keys.point][feature]["data"], feature_values[i]])
+
+        elif overwrite or (feature not in target_point_cloud[keys.point]):
+            target_point_cloud[keys.point][feature] = {
+                "type": 'float64', "data": feature_values[i]}
+
+
+def _add_or_update_feature_one_by_one(env_point_cloud, extractor, feature_values, n_features, n_targets, neighborhoods,
+                                      target_idx_base, target_point_cloud, volume):
     for target_index in range(n_targets):
         point_values = extractor.extract(env_point_cloud, neighborhoods[target_index], target_point_cloud,
-                                         target_index+target_idx_base, volume)
+                                         target_index + target_idx_base, volume)
         if n_features > 1:
             for i in range(n_features):
                 feature_values[i][target_index] = point_values[i]
         else:
             feature_values[0][target_index] = point_values
-    for i in range(n_features):
-        feature = provided_features[i]
-        if (target_idx_base != 0):
-            target_point_cloud[keys.point][feature]["data"] = np.append(target_point_cloud[keys.point][feature]["data"], feature_values[i])
-        elif (overwrite or (feature not in target_point_cloud[keys.point])) and (target_idx_base == 0):
-            target_point_cloud[keys.point][feature] = {
-                "type": 'float64', "data": feature_values[i]}
 
 
-def _make_feature_list(feature_names):
-    feature_list = reversed(_make_feature_list_helper(feature_names))
+def _add_or_update_feature_in_chunks(env_point_cloud, extractor, feature_values, n_features, n_targets, neighborhoods,
+                                     target_idx_base, target_point_cloud, volume):
+    chunk_size = 100000
+    print('calculating {} in chunks'.format(extractor.provides()))
+    for chunk_no in range(int(np.math.ceil(n_targets / chunk_size))):
+        i_start = chunk_no * chunk_size
+        i_end = min((chunk_no + 1) * chunk_size, n_targets)
+        target_indices = np.arange(i_start, i_end)
+        point_values = extractor.extract(env_point_cloud, neighborhoods[i_start:i_end], target_point_cloud,
+                                         target_indices + target_idx_base, volume)
+
+        if n_features > 1:
+            for i in range(n_features):
+                feature_values[i][target_indices] = point_values[i]
+        else:
+            feature_values[0][target_indices] = point_values
+
+
+def _make_extended_feature_list(feature_names):
+    feature_list = reversed(_make_extended_feature_list_helper(feature_names))
+    return _remove_duplicates(feature_list)
+
+
+def _remove_duplicates(feature_list):
     seen = set()
     return [f for f in feature_list if not (f in seen or seen.add(f))]
 
 
-def _make_feature_list_helper(feature_names):
+def _make_extended_feature_list_helper(feature_names):
     feature_list = feature_names
     for feature_name in feature_names:
         extractor = FEATURES[feature_name]()
         dependencies = extractor.requires()
         feature_list.extend(dependencies)
-        feature_list.extend(_make_feature_list_helper(dependencies))
+        feature_list.extend(_make_extended_feature_list_helper(dependencies))
     return feature_list

laserchicken/feature_extractor/density_absolute_mean_norm_z_feature_extractor.py

@@ -0,0 +1,33 @@
+"""Pulse penetration ratio and density absolute mean calculations.
+
+See https://github.com/eEcoLiDAR/eEcoLiDAR/issues/23.
+"""
+
+import numpy as np
+
+from laserchicken.feature_extractor.abc import AbstractFeatureExtractor
+from laserchicken.feature_extractor.density_absolute_mean_z_feature_extractor import \
+    DensityAbsoluteMeanZFeatureExtractor
+from laserchicken.keys import point, normalized_height
+
+# classification according to
+# http://www.asprs.org/wp-content/uploads/2010/12/LAS_1-4_R6.pdf
+GROUND_TAGS = [2]
+
+
+class DensityAbsoluteMeanNormZFeatureExtractor(DensityAbsoluteMeanZFeatureExtractor):
+    """Feature extractor for the point density."""
+    DATA_KEY = normalized_height
+
+    @classmethod
+    def provides(cls):
+        """
+        Get a list of names of the feature values.
+
+        This will return as many names as the number feature values that will be returned.
+        For instance, if a feature extractor returns the first 3 eigen values, this method
+        should return 3 names, for instance 'eigen_value_1', 'eigen_value_2' and 'eigen_value_3'.
+
+        :return: List of feature names
+        """
+        return ['density_absolute_mean_norm_z']