Merge branch 'development' into entropy_z_crashing_103

laserchicken/compute_neighbors.py

@@ -1,49 +1,95 @@
+import sys
+import math
+import numpy as np
+from scipy.spatial import cKDTree
+from psutil import virtual_memory
+
 from laserchicken import utils, kd_tree
 from laserchicken.volume_specification import Sphere, InfiniteCylinder
+from laserchicken.keys import point
+
+
+def frange(x_value, y_value, jump):
+    while x_value < y_value:
+        yield x_value
+        x_value += jump
+
+
+MEMORY_THRESHOLD = 0.5
+POINTCLOUD_DIST = 10
 
 
 def compute_cylinder_neighborhood(environment_pc, target_pc, radius):
-    """Find the indices of points within a cylindrical neighbourhood (using KD Tree)
-    for a given point of a target point cloud among the points from an environment point cloud.
+    """Find the indices of points within a cylindrical neighbourhood (using KD Tree) for a given point of a target
+    point cloud among the points from an environment point cloud.
 
     :param environment_pc: environment point cloud
     :param target_pc: point cloud that contains the points at which neighborhoods are to be calculated
     :param radius: search radius for neighbors
     :return: indices of neighboring points from the environment point cloud for each target point
              the returned indices also contains the index of the target point.
     """
+    avg_points_cyl = (radius * radius * math.pi) * POINTCLOUD_DIST
+    x = target_pc[point]['x']['data']
+
+    cyl_size = avg_points_cyl * np.size(x) * sys.getsizeof(int)
+    mem_size = virtual_memory().total
 
-    env_tree = kd_tree.get_kdtree_for_pc(environment_pc)
-    target_tree = kd_tree.get_kdtree_for_pc(target_pc)
-    return target_tree.query_ball_tree(env_tree, radius)
+    print("Cylinder size in Bytes: %s" % cyl_size)
+    print("Memory size in Bytes: %s" % mem_size)
+
+    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: %d" % num_points)
+
+        env_tree = kd_tree.get_kdtree_for_pc(environment_pc)
+
+        for i in range(0, np.size(x), num_points):
+            box_points = np.column_stack(
+                (x[i:min(i + num_points, np.size(x))], y[i:min(i + num_points, np.size(x))]))
+            target_box_tree = cKDTree(
+                box_points, compact_nodes=False, balanced_tree=False)
+            yield target_box_tree.query_ball_tree(env_tree, radius)
+
+    else:
+        print("Start tree creation")
+        env_tree = kd_tree.get_kdtree_for_pc(environment_pc)
+        print("Done with env tree creation")
+        target_tree = kd_tree.get_kdtree_for_pc(target_pc)
+        print("Done with target tree creation")
+        yield target_tree.query_ball_tree(env_tree, radius)
 
 
 def compute_sphere_neighborhood(environment_pc, target_pc, radius):
     """
-    Find the indices of points within a spherical neighbourhood for a given point of a target point cloud among
-    the points from an environment point cloud.
+    Find the indices of points within a spherical neighbourhood for a given point of a target point cloud among the
+    points from an environment point cloud.
 
     :param environment_pc: environment point cloud
     :param target_pc: point cloud that contains the points at which neighborhoods are to be calculated
     :param radius: search radius for neighbors
     :return: indices of neighboring points from the environment point cloud for each target point
     """
+    neighbors = compute_cylinder_neighborhood(
+        environment_pc, target_pc, radius)
 
-    neighborhood_indices = compute_cylinder_neighborhood(environment_pc, target_pc, radius)
-
-    result = []
-    for i in range(len(neighborhood_indices)):
-        target_x, target_y, target_z = utils.get_point(target_pc, i)
-        neighbor_indices = neighborhood_indices[i]
-        result_indices = []
-        for j in neighbor_indices:
-            env_x, env_y, env_z = utils.get_point(environment_pc, j)
-            if abs(target_z - env_z) > radius:
-                continue
-            if (env_x - target_x) ** 2 + (env_y - target_y) ** 2 + (env_z - target_z) ** 2 <=  radius ** 2 :
-                result_indices.append(j)
-        result.append(result_indices)
-    return result
+    for neighborhood_indices in neighbors:
+        result = []
+        for i in range(len(neighborhood_indices)):
+            target_x, target_y, target_z = utils.get_point(target_pc, i)
+            neighbor_indices = neighborhood_indices[i]
+            result_indices = []
+            for j in neighbor_indices:
+                env_x, env_y, env_z = utils.get_point(environment_pc, j)
+                if abs(target_z - env_z) > radius:
+                    continue
+                if (env_x - target_x) ** 2 + (env_y - target_y) ** 2 + (env_z - target_z) ** 2 <= radius ** 2:
+                    result_indices.append(j)
+            result.append(result_indices)
+        yield result
 
 
 def compute_neighborhoods(env_pc, target_pc, volume_description):
@@ -56,8 +102,15 @@ 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:
-        return compute_sphere_neighborhood(env_pc, target_pc, volume_description.radius)
+        neighbors = compute_sphere_neighborhood(
+            env_pc, target_pc, volume_description.radius)
     elif volume_type == InfiniteCylinder.TYPE:
-        return compute_cylinder_neighborhood(env_pc, target_pc, volume_description.radius)
-    raise ValueError('Neighborhood computation error because volume type "{}" is unknown.'.format(volume_type))
+        neighbors = compute_cylinder_neighborhood(
+            env_pc, target_pc, volume_description.radius)
+    else:
+        raise ValueError(
+            'Neighborhood computation error because volume type "{}" is unknown.'.format(volume_type))
+    for x in neighbors:
+        yield x

laserchicken/feature_extractor/__init__.py

@@ -40,7 +40,10 @@ def compute_features(env_point_cloud, neighborhoods, target_point_cloud, feature
     >>> point_cloud = read_ply.read('data1.ply')
     >>> target_point_cloud = read_ply.read('data2.ply')
     >>> volume = volume_specification.InfiniteCylinder(4)
-    >>> neighborhoods = compute_neighborhoods(point_cloud, target_point_cloud, volume)
+    >>> neighbors = compute_neighborhoods(point_cloud, target_point_cloud, volume)
+    >>> neighborhoods = []
+    >>> for x in neighbors:
+    >>>   neighborhoods += x
     >>> compute_features(point_cloud, neighborhoods, target_point_cloud, ['eigenv_1', 'kurto_z'], volume)
     >>> eigenv_1 = target_point_cloud[point]['eigenv_1']['data']
 
@@ -69,8 +72,10 @@ def compute_features(env_point_cloud, neighborhoods, target_point_cloud, feature
             start = time.time()
 
         extractor = FEATURES[feature]()
-        _add_or_update_feature(env_point_cloud, neighborhoods, target_point_cloud, extractor, volume, overwrite, kwargs)
-        utils.add_metadata(target_point_cloud, type(extractor).__module__, extractor.get_params())
+        _add_or_update_feature(env_point_cloud, neighborhoods,
+                               target_point_cloud, extractor, volume, overwrite, kwargs)
+        utils.add_metadata(target_point_cloud, type(
+            extractor).__module__, extractor.get_params())
 
         if verbose:
             elapsed = time.time() - start
@@ -91,7 +96,8 @@ def _add_or_update_feature(env_point_cloud, neighborhoods, target_point_cloud, e
         setattr(extractor, k, kwargs[k])
     provided_features = extractor.provides()
     n_features = len(provided_features)
-    feature_values = [np.empty(n_targets, dtype=np.float64) for i in range(n_features)]
+    feature_values = [np.empty(n_targets, dtype=np.float64)
+                      for i in range(n_features)]
     for target_index in range(n_targets):
         point_values = extractor.extract(env_point_cloud, neighborhoods[target_index], target_point_cloud,
                                          target_index, volume)
@@ -103,7 +109,8 @@ def _add_or_update_feature(env_point_cloud, neighborhoods, target_point_cloud, e
     for i in range(n_features):
         feature = provided_features[i]
         if overwrite or (feature not in target_point_cloud[keys.point]):
-            target_point_cloud[keys.point][feature] = {"type": 'float64', "data": feature_values[i]}
+            target_point_cloud[keys.point][feature] = {
+                "type": 'float64', "data": feature_values[i]}
 
 
 def _make_feature_list(feature_names):

laserchicken/feature_extractor/abc.py

@@ -14,7 +14,8 @@ def requires(cls):
 
         :return: List of feature names
         """
-        raise NotImplementedError("Class %s doesn't implement get_requirements()" % (cls.__name__))
+        raise NotImplementedError(
+            "Class %s doesn't implement get_requirements()" % (cls.__name__))
 
     @classmethod
     def provides(cls):
@@ -27,7 +28,8 @@ def provides(cls):
 
         :return: List of feature names
         """
-        raise NotImplementedError("Class %s doesn't implement get_names()" % (cls.__name__))
+        raise NotImplementedError(
+            "Class %s doesn't implement get_names()" % (cls.__name__))
 
     def extract(self, point_cloud, neighborhood, target_point_cloud, target_index, volume_description):
         """
@@ -40,7 +42,8 @@ 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
         """
-        raise NotImplementedError("Class %s doesn't implement extract_features()" % (self.__class__.__name__))
+        raise NotImplementedError(
+            "Class %s doesn't implement extract_features()" % (self.__class__.__name__))
 
     def get_params(self):
         """

laserchicken/feature_extractor/echo_ratio_feature_extractor.py

@@ -62,7 +62,8 @@ def extract(self, point_cloud, neighborhood, target_point_cloud, target_index, v
 
         xyz0 = self.get_target_position(target_point_cloud, target_index)
 
-        n_sphere = np.sum(np.sum((xyz - xyz0) ** 2, 1) <= volume_description.radius ** 2)
+        n_sphere = np.sum(np.sum((xyz - xyz0) ** 2, 1) <=
+                          volume_description.radius ** 2)
         return n_sphere / n_cylinder * 100.
 
     @staticmethod

laserchicken/feature_extractor/entropy_feature_extractor.py

@@ -34,7 +34,8 @@ def extract(self, source_pc, neighborhood, target_pc, target_index, volume_descr
         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/normal_plane_feature_extractor.py

@@ -3,6 +3,7 @@
 from laserchicken.keys import point
 from laserchicken.utils import fit_plane_svd
 
+
 class NormalPlaneFeatureExtractor(AbstractFeatureExtractor):
     """Feature extractor for the normal and slope of a plane."""
 
@@ -29,10 +30,9 @@ def provides(cls):
 
         :return: List of feature names
         """
-        return ['normal_vector_1','normal_vector_2','normal_vector_3','slope']
+        return ['normal_vector_1', 'normal_vector_2', 'normal_vector_3', 'slope']
 
     def extract(self, sourcepc, neighborhood, targetpc, targetindex, volume):
-
         """
         Extract the feature value(s) of the point cloud at location of the target.
 
@@ -48,10 +48,10 @@ def extract(self, sourcepc, neighborhood, targetpc, targetindex, volume):
         y = sourcepc[point]['y']['data'][neighborhood]
         z = sourcepc[point]['z']['data'][neighborhood]
 
-        nvect = fit_plane_svd(x,y,z)
-        slope = np.dot(nvect,np.array([0.,0.,1.]))
+        nvect = fit_plane_svd(x, y, z)
+        slope = np.dot(nvect, np.array([0., 0., 1.]))
 
-        return nvect[0],nvect[1],nvect[2],slope
+        return nvect[0], nvect[1], nvect[2], slope
 
     def get_params(self):
         """

laserchicken/feature_extractor/pulse_penetration_feature_extractor.py

@@ -52,11 +52,15 @@ 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
         """
-        class_neighbors = np.array(point_cloud[point]['raw_classification']["data"])[neighborhood]
-        ground_indices = self._get_ground_indices(class_neighbors, self.ground_tags)
+        class_neighbors = np.array(point_cloud[point]['raw_classification']["data"])[
+            neighborhood]
+        ground_indices = self._get_ground_indices(
+            class_neighbors, self.ground_tags)
 
-        pulse_penetration_ratio = self._get_pulse_penetration_ratio(ground_indices, class_neighbors)
-        density_absolute_mean = self._get_density_absolute_mean(ground_indices, point_cloud)
+        pulse_penetration_ratio = self._get_pulse_penetration_ratio(
+            ground_indices, class_neighbors)
+        density_absolute_mean = self._get_density_absolute_mean(
+            ground_indices, point_cloud)
 
         return pulse_penetration_ratio, density_absolute_mean
 
@@ -81,7 +85,8 @@ def _get_density_absolute_mean(ground_indices, source_point_cloud):
         if n_ground == 0:
             density_absolute_mean = 0.
         else:
-            density_absolute_mean = float(len(z_ground[z_ground > np.mean(z_ground)])) / n_ground * 100.
+            density_absolute_mean = float(
+                len(z_ground[z_ground > np.mean(z_ground)])) / n_ground * 100.
         return density_absolute_mean
 
     def get_params(self):