Merge pull request #116 from eEcoLiDAR/development

Development

Author: cwmeijer

.idea/eEcoLiDAR.iml

@@ -0,0 +1,18 @@
+# Changelog
+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]
+### Added
+- General tests that all current and future feature extractors will be checked against.
+
+## Changed
+
+## Fixed
+- Fixed many feature extractors for corner cases (e.g. zero points)
+
+## Removed
+
+## [0.1.0] - 2018-04-17

.idea/misc.xml

@@ -41,9 +41,8 @@ def compute_cylinder_neighborhood(environment_pc, target_pc, radius):
     if cyl_size > mem_size * MEMORY_THRESHOLD:
         y = target_pc[point]['y']['data']
 
-        num_points = math.floor(mem_size * MEMORY_THRESHOLD / \
-            (avg_points_cyl * sys.getsizeof(int)))
-        print("Number of points: %f" % num_points)
+        num_points = int(math.floor(mem_size * MEMORY_THRESHOLD / (avg_points_cyl * sys.getsizeof(int))))
+        print("Number of points: %d" % num_points)
 
         env_tree = kd_tree.get_kdtree_for_pc(environment_pc)
 
@@ -102,7 +101,6 @@ def compute_neighborhoods(env_pc, target_pc, volume_description):
     :return: indices of neighboring points from the environment point cloud for each target point
     """
     volume_type = volume_description.get_type()
-    neighbors = []
     if volume_type == Sphere.TYPE:
         neighbors = compute_sphere_neighborhood(
             env_pc, target_pc, volume_description.radius)

CHANGELOG.md

@@ -31,29 +31,20 @@ def provides(cls):
         """
         return ['point_density']
 
-    def extract(self, sourcepc, neighborhood, targetpc, targetindex, volume):
+    def extract(self, source_point_cloud, neighborhood, target_point_cloud, target_index, volume):
         """
-        Extract the feature value(s) of the point cloud at location of the target.
+        Extract either the surface density or volume density depending on the volume type.
 
-        :param point_cloud: environment (search space) point cloud
+        :param source_point_cloud: environment (search space) point cloud
         :param neighborhood: array of indices of points within the point_cloud argument
         :param target_point_cloud: point cloud that contains target point
         :param target_index: index of the target point in the target pointcloud
-        :param volume_description: volume object that describes the shape and size of the search volume
+        :param volume: volume object that describes the shape and size of the search volume
         :return: feature value
         """
 
-        if sourcepc is not None and isinstance(neighborhood, list):
-            npts = float(len(sourcepc[point]['x']['data'][neighborhood]))
+        npts = float(len(neighborhood))
 
-        elif targetpc is not None:
-            npts = float(len(targetpc[point]['x']['data']))
-        else:
-            raise ValueError("You can either specify a sourcepc and a neighborhhod or a targetpc\n\
-                              example\nextractror.extract(sourcepc,index,None,None,volume)\n\
-                              extractror.extract(None,None,targetpc,None,volume)")
-
-        # sphere/cylinder cases
         if volume.get_type() == Sphere.TYPE:
             vol = volume.calculate_volume()
             return npts / vol

laserchicken/compute_neighbors.py

@@ -28,14 +28,15 @@ def get_params(self):
         return p
 
     def extract(self, source_pc, neighborhood, target_pc, target_index, volume_description):
+        if len(neighborhood) == 0:
+            return 0
         z = source_pc[keys.point]["z"]["data"][neighborhood]
         _z_min = np.min(z) if self.z_min is None else self.z_min
         _z_max = np.max(z) if self.z_max is None else self.z_max
-        if (_z_min == _z_max):
-          return 0
+        if _z_min == _z_max:
+            return 0
         n_bins = int(np.ceil((_z_max - _z_min) / self.layer_thickness))
-        data = np.histogram(z, bins=n_bins, range=(
-            _z_min, _z_max), density=True)[0]
+        data = np.histogram(z, bins=n_bins, range=(_z_min, _z_max), density=True)[0]
         entropy_func = np.vectorize(_x_log_2x)
         norm = np.sum(data)
         return -(entropy_func(data / norm)).sum()

laserchicken/feature_extractor/density_feature_extractor.py

@@ -6,6 +6,9 @@
 
 
 class HeightStatisticsFeatureExtractor(AbstractFeatureExtractor):
+    DEFAULT_MAX = float('NaN')
+    DEFAULT_MIN = float('NaN')
+
     @classmethod
     def requires(cls):
         return []
@@ -15,15 +18,18 @@ def provides(cls):
         return ['max_z', 'min_z', 'mean_z', 'median_z', 'std_z', 'var_z', 'range', 'coeff_var_z', 'skew_z', 'kurto_z']
 
     def extract(self, sourcepc, neighborhood, targetpc, targetindex, volume_description):
-        z = sourcepc[point]['z']['data'][neighborhood]
-        max_z = np.max(z)
-        min_z = np.min(z)
-        mean_z = np.mean(z)
-        median_z = np.median(z)
-        std_z = np.std(z)
-        var_z = np.var(z)
-        range_z = max_z - min_z
-        coeff_var_z = np.std(z) / np.mean(z)
-        skew_z = stat.skew(z)
-        kurto_z = stat.kurtosis(z)
+        if neighborhood:
+            z = sourcepc[point]['z']['data'][neighborhood]
+            max_z = np.max(z) if len(z) > 0 else self.DEFAULT_MAX
+            min_z = np.min(z) if len(z) > 0 else self.DEFAULT_MIN
+            mean_z = np.mean(z)
+            median_z = np.median(z)
+            std_z = np.std(z)
+            var_z = np.var(z)
+            range_z = max_z - min_z
+            coeff_var_z = np.std(z) / np.mean(z)
+            skew_z = stat.skew(z)
+            kurto_z = stat.kurtosis(z)
+        else:
+            max_z = min_z = mean_z = median_z = std_z = var_z = range_z = coeff_var_z = skew_z = kurto_z = np.NaN
         return max_z, min_z, mean_z, median_z, std_z, var_z, range_z, coeff_var_z, skew_z, kurto_z

laserchicken/feature_extractor/entropy_feature_extractor.py

@@ -32,26 +32,30 @@ def provides(cls):
         """
         return ['normal_vector_1', 'normal_vector_2', 'normal_vector_3', 'slope']
 
-    def extract(self, sourcepc, neighborhood, targetpc, targetindex, volume):
+    def extract(self, point_cloud, neighborhood, target_point_cloud, target_index, volume_description):
         """
         Extract the feature value(s) of the point cloud at location of the target.
 
         :param point_cloud: environment (search space) point cloud
         :param neighborhood: array of indices of points within the point_cloud argument
         :param target_point_cloud: point cloud that contains target point
-        :param target_index: index of the target point in the target pointcloud
+        :param target_index: index of the target point in the target point cloud
         :param volume_description: volume object that describes the shape and size of the search volume
         :return: feature value
         """
 
-        x = sourcepc[point]['x']['data'][neighborhood]
-        y = sourcepc[point]['y']['data'][neighborhood]
-        z = sourcepc[point]['z']['data'][neighborhood]
+        x = _to_float64(point_cloud[point]['x']['data'][neighborhood])
+        y = _to_float64(point_cloud[point]['y']['data'][neighborhood])
+        z = _to_float64(point_cloud[point]['z']['data'][neighborhood])
 
-        nvect = fit_plane_svd(x, y, z)
-        slope = np.dot(nvect, np.array([0., 0., 1.]))
+        try:
+            normal_vector = fit_plane_svd(x, y, z)
+            slope = np.dot(normal_vector, np.array([0., 0., 1.]))
+        except np.linalg.linalg.LinAlgError:
+            normal_vector = [0.0, 0.0, 0.0]
+            slope = 0.0
 
-        return nvect[0], nvect[1], nvect[2], slope
+        return normal_vector[0], normal_vector[1], normal_vector[2], slope
 
     def get_params(self):
         """
@@ -60,3 +64,13 @@ def get_params(self):
         Needed for provenance.
         """
         return ()
+
+
+def _to_float64(x_org):
+    """
+    Returns the input if dtype is float64 or a copy of it with this dtype if it isn't.
+    NB. Returns reference to input array or copy.
+    :param x_org:
+    :return: array with float64 dtype
+    """
+    return x_org if x_org.dtype == np.float64 else np.array(x_org, dtype=np.float64)

laserchicken/feature_extractor/height_statistics_feature_extractor.py

@@ -8,14 +8,14 @@
 from laserchicken.feature_extractor.abc import AbstractFeatureExtractor
 from laserchicken.keys import point
 
+# classification according to
+# http://www.asprs.org/wp-content/uploads/2010/12/LAS_1-4_R6.pdf
+GROUND_TAGS = [2]
+
 
 class PulsePenetrationFeatureExtractor(AbstractFeatureExtractor):
     """Feature extractor for the point density."""
 
-    # classification according to
-    # http://www.asprs.org/wp-content/uploads/2010/12/LAS_1-4_R6.pdf
-    ground_tags = [2]
-
     @classmethod
     def requires(cls):
         """
@@ -52,10 +52,14 @@ def extract(self, point_cloud, neighborhood, target_point_cloud, target_index, v
         :param volume_description: volume object that describes the shape and size of the search volume
         :return: feature value
         """
+        if 'raw_classification' not in point_cloud[point]:
+            raise ValueError(
+                'Missing raw_classification attribute which is necessary for calculating pulse_penetratio and '
+                'density_absolute_mean features.')
+
         class_neighbors = [point_cloud[point]['raw_classification']["data"][n] for n in neighborhood]
 
-        ground_indices = self._get_ground_indices(
-            class_neighbors, self.ground_tags)
+        ground_indices = self._get_ground_indices(class_neighbors, GROUND_TAGS)
 
         pulse_penetration_ratio = self._get_pulse_penetration_ratio(
             ground_indices, class_neighbors)
@@ -74,7 +78,7 @@ def _get_ground_indices(point_cloud, ground_tags):
 
     @staticmethod
     def _get_pulse_penetration_ratio(ground_indices, class_neighbors):
-        n_total = len(class_neighbors)
+        n_total = np.max((len(class_neighbors), 1))
         n_ground = len(ground_indices)
         return float(n_ground) / n_total
 

laserchicken/feature_extractor/normal_plane_feature_extractor.py

@@ -4,6 +4,7 @@
 See https://github.com/eEcoLiDAR/eEcoLiDAR/issues/20
 """
 import numpy as np
+from numpy.linalg import LinAlgError
 
 from laserchicken.feature_extractor.abc import AbstractFeatureExtractor
 from laserchicken.utils import get_point, fit_plane
@@ -49,9 +50,12 @@ def extract(self, source_point_cloud, neighborhood, target_point_cloud, target_i
         :return:
         """
         x, y, z = get_point(source_point_cloud, neighborhood)
-        plane_estimator = fit_plane(x, y, z)
-        normalized = z - plane_estimator(x, y)
-        return np.std(normalized)
+        try:
+            plane_estimator = fit_plane(x, y, z)
+            normalized = z - plane_estimator(x, y)
+            return np.std(normalized)
+        except LinAlgError:
+            return 0
 
     def get_params(self):
         """