Merge pull request #119 from eEcoLiDAR/development

Development

laserchicken/_version.py

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

laserchicken/compute_neighbors.py

@@ -5,7 +5,7 @@
 from psutil import virtual_memory
 
 from laserchicken import utils, kd_tree
-from laserchicken.volume_specification import Sphere, InfiniteCylinder
+from laserchicken.volume_specification import Sphere, InfiniteCylinder, Cell, Cube
 from laserchicken.keys import point
 
 
@@ -16,7 +16,7 @@ def frange(x_value, y_value, jump):
 
 
 MEMORY_THRESHOLD = 0.5
-POINTCLOUD_DIST = 10
+POINT_CLOUD_DIST = 10
 
 
 def compute_cylinder_neighborhood(environment_pc, target_pc, radius):
@@ -29,9 +29,14 @@ def compute_cylinder_neighborhood(environment_pc, target_pc, radius):
     :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
+    avg_points_cyl = (radius * radius * math.pi) * POINT_CLOUD_DIST
     x = target_pc[point]['x']['data']
 
+    if len(environment_pc[point]['x']['data']) == 0:
+        for _ in x:
+            yield [[] for _ in x]
+        return
+
     cyl_size = avg_points_cyl * np.size(x) * sys.getsizeof(int)
     mem_size = virtual_memory().total
 
@@ -72,16 +77,15 @@ def compute_sphere_neighborhood(environment_pc, target_pc, radius):
     :param radius: search radius for neighbors
     :return: indices of neighboring points from the environment point cloud for each target point
     """
-    neighbors = compute_cylinder_neighborhood(
+    neighborhoods = compute_cylinder_neighborhood(
         environment_pc, target_pc, radius)
 
-    for neighborhood_indices in neighbors:
+    for neighborhood_indices in neighborhoods:
         result = []
-        for i in range(len(neighborhood_indices)):
+        for i, _ in enumerate(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:
+            for j in neighborhood_indices[i]:
                 env_x, env_y, env_z = utils.get_point(environment_pc, j)
                 if abs(target_z - env_z) > radius:
                     continue
@@ -91,6 +95,68 @@ def compute_sphere_neighborhood(environment_pc, target_pc, radius):
         yield result
 
 
+def compute_cell_neighborhood(environment_pc, target_pc, side_length):
+    """
+    Find the indices of points within a square 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 side_length: search radius for neighbors
+    :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))
+
+    neighbors = compute_cylinder_neighborhood(
+        environment_pc, target_pc, max_radius)
+
+    for neighborhood_indices in neighbors:
+        result = []
+        for i, _ in enumerate(neighborhood_indices):
+            target_x, target_y, _ = utils.get_point(target_pc, i)
+            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):
+                    continue
+                else:
+                    result_indices.append(j)
+            result.append(result_indices)
+        yield result
+
+
+def compute_cube_neighborhood(environment_pc, target_pc, side_length):
+    """
+    Find the indices of points within a square 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 side_length: search radius for neighbors
+    :return: indices of neighboring points from the environment point cloud for each target point
+    """
+
+    neighbors = compute_cell_neighborhood(
+        environment_pc, target_pc, side_length)
+
+    for neighborhood_indices in neighbors:
+        result = []
+        for i, _ in enumerate(neighborhood_indices):
+            _, _, target_z = utils.get_point(target_pc, i)
+            neighbor_indices = neighborhood_indices[i]
+            result_indices = []
+            for j in neighbor_indices:
+                _, _, env_z = utils.get_point(environment_pc, j)
+                if abs(target_z - env_z) > side_length:
+                    continue
+                else:
+                    result_indices.append(j)
+            result.append(result_indices)
+        yield result
+
+
 def compute_neighborhoods(env_pc, target_pc, volume_description):
     """
     Find a subset of points in a neighbourhood in the environment point cloud for each point in a target point cloud.
@@ -101,14 +167,23 @@ 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()
-    if volume_type == Sphere.TYPE:
-        neighbors = compute_sphere_neighborhood(
-            env_pc, target_pc, volume_description.radius)
+
+    if volume_type == Cell.TYPE:
+        neighbors1 = compute_cell_neighborhood(env_pc, target_pc, volume_description.side_length)
+        for x in neighbors1:
+            yield x
+    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)
+        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 x

laserchicken/feature_extractor/density_feature_extractor.py

@@ -31,27 +31,32 @@ def provides(cls):
         """
         return ['point_density']
 
-    def extract(self, source_point_cloud, neighborhood, target_point_cloud, target_index, volume):
+    def extract(self, source_pc, neighborhood, target_pc, target_index, volume_description):
         """
-        Extract either the surface density or volume density depending on the volume type.
+        Extract the feature value(s) of the point cloud at location of the target.
 
-        :param source_point_cloud: environment (search space) point cloud
+        :param source_pc: 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_pc: point cloud that contains target point
         :param target_index: index of the target point in the target pointcloud
-        :param volume: volume object that describes the shape and size of the search volume
+        :param volume_description: volume object that describes the shape and size of the search volume
         :return: feature value
         """
 
-        npts = float(len(neighborhood))
+        if source_pc is not None and isinstance(neighborhood, list):
+            n_points = float(len(source_pc[point]['x']['data'][neighborhood]))
+
+        elif target_pc is not None:
+            n_points = float(len(target_pc[point]['x']['data']))
+        else:
+            raise ValueError("You can either specify a sourcepc and a neighborhood or a targetpc\n\
+                              example\nextractror.extract(sourcepc,index,None,None,volume)\n\
+                              extractror.extract(None,None,targetpc,None,volume)")
+
+        area_or_volume = volume_description.calculate_area_or_volume()
+        return n_points / area_or_volume
 
-        if volume.get_type() == Sphere.TYPE:
-            vol = volume.calculate_volume()
-            return npts / vol
 
-        elif volume.get_type() == InfiniteCylinder.TYPE:
-            area = volume.calculate_base_area()
-            return npts / area
 
     def get_params(self):
         """

laserchicken/feature_extractor/test_all_features.py

@@ -8,7 +8,7 @@
 from laserchicken import read_las
 from laserchicken.feature_extractor.pulse_penetration_feature_extractor import GROUND_TAGS
 from laserchicken.keys import point
-from laserchicken.utils import copy_pointcloud
+from laserchicken.utils import copy_point_cloud
 from laserchicken.volume_specification import InfiniteCylinder
 
 from . import _feature_map
@@ -22,9 +22,9 @@
 
 _CYLINDER = InfiniteCylinder(4)
 _PC_260807 = read_las.read(os.path.join(_TEST_DATA_SOURCE, _TEST_FILE_NAME))
-_PC_1000 = copy_pointcloud(_PC_260807, array_mask=(
+_PC_1000 = copy_point_cloud(_PC_260807, array_mask=(
     np.random.choice(range(len(_PC_260807[keys.point]['x']['data'])), size=1000, replace=False)))
-_PC_10 = copy_pointcloud(_PC_260807, array_mask=(
+_PC_10 = copy_point_cloud(_PC_260807, array_mask=(
     np.random.choice(range(len(_PC_260807[keys.point]['x']['data'])), size=10, replace=False)))
 _1000_NEIGHBORHOODS_IN_260807 = next(compute_neighbors.compute_neighborhoods(_PC_260807, _PC_1000, _CYLINDER))
 _10_NEIGHBORHOODS_IN_260807 = next(compute_neighbors.compute_neighborhoods(_PC_260807, _PC_10, _CYLINDER))
@@ -35,16 +35,16 @@
 
 @pytest.mark.parametrize("feature", feature_names)
 def test_completeTile_consistentOutput(feature):
-    target_point_cloud = copy_pointcloud(_PC_1000)
-    compute_features(copy_pointcloud(_PC_260807), _1000_NEIGHBORHOODS_IN_260807, 0, target_point_cloud,
+    target_point_cloud = copy_point_cloud(_PC_1000)
+    compute_features(copy_point_cloud(_PC_260807), _1000_NEIGHBORHOODS_IN_260807, 0, target_point_cloud,
                      [feature], volume=_CYLINDER)
     _assert_consistent_attribute_length(target_point_cloud)
 
 
 @pytest.mark.parametrize("feature", feature_names)
 def test_manyTargets_consistentOutput(feature):
-    target_point_cloud = copy_pointcloud(_PC_260807)
-    compute_features(copy_pointcloud(_PC_10), _260807_NEIGHBORHOODS_IN_10, 0, target_point_cloud,
+    target_point_cloud = copy_point_cloud(_PC_260807)
+    compute_features(copy_point_cloud(_PC_10), _260807_NEIGHBORHOODS_IN_10, 0, target_point_cloud,
                      [feature], volume=_CYLINDER)
     _assert_consistent_attribute_length(target_point_cloud)
 
@@ -106,9 +106,9 @@ def _test_consistent_output_with_n_neighbors(feature, n_neighbors):
 @pytest.mark.parametrize("feature", feature_names)
 def test_inputNotChanged(feature):
     original_environment = _PC_260807
-    environment = copy_pointcloud(original_environment)
+    environment = copy_point_cloud(original_environment)
     original_targets = _PC_10
-    targets = copy_pointcloud(original_targets)
+    targets = copy_point_cloud(original_targets)
     original_neighborhoods = _10_NEIGHBORHOODS_IN_260807
     neighborhoods = [[e for e in l] for l in original_neighborhoods]
 

laserchicken/feature_extractor/test_density_extractor.py

@@ -6,8 +6,9 @@
 
 from laserchicken import keys, read_las, utils
 from laserchicken.compute_neighbors import compute_neighborhoods
-from laserchicken.volume_specification import Sphere, InfiniteCylinder
+from laserchicken.volume_specification import Sphere, InfiniteCylinder, Cell, Cube
 from laserchicken.feature_extractor.density_feature_extractor import PointDensityFeatureExtractor
+from laserchicken.test_tools import create_point_cloud
 
 
 class TestDensityFeatureExtractorSphere(unittest.TestCase):
@@ -17,12 +18,8 @@ class TestDensityFeatureExtractorSphere(unittest.TestCase):
 
     def test_sphere(self):
         """Compute the density for a sphere given as index of the source pc."""
-        neighbors = compute_neighborhoods(self.point_cloud,
-                                          self.targetpc,
-                                          self.sphere)
-        neighbors_index = []
-        for x in neighbors:
-            neighbors_index += x
+        neighbors_index = list(compute_neighborhoods(self.point_cloud, self.targetpc, self.sphere))
+
         extractor = PointDensityFeatureExtractor()
         for index in neighbors_index:
             d = extractor.extract(self.point_cloud, index,
@@ -31,7 +28,7 @@ def test_sphere(self):
 
     def _get_central_point(self):
         """Get the central point."""
-        return utils.copy_pointcloud(self.point_cloud, [0])
+        return utils.copy_point_cloud(self.point_cloud, [0])
 
     def _set_sphere_data(self):
         """Generate a pc of points regularly positionned on a two spheres of radius 1 and 2."""
@@ -81,12 +78,8 @@ class TestDensityFeatureExtractorCylinder(unittest.TestCase):
 
     def test_cylinder(self):
         """Compute the density for a cylinder given as index of source pc."""
-        neighbors = compute_neighborhoods(self.point_cloud,
-                                          self.targetpc,
-                                          self.cyl)
-        neighbors_index = []
-        for x in neighbors:
-            neighbors_index += x
+        neighbors_index = compute_neighborhoods(self.point_cloud, self.targetpc, self.cyl)
+
         extractor = PointDensityFeatureExtractor()
         for index in neighbors_index:
             d = extractor.extract(self.point_cloud, index,
@@ -95,7 +88,7 @@ def test_cylinder(self):
 
     def _get_central_point(self):
         """Get the central point."""
-        return utils.copy_pointcloud(self.point_cloud, [0])
+        return utils.copy_point_cloud(self.point_cloud, [0])
 
     def _set_cylinder_data(self):
 
@@ -141,6 +134,40 @@ def tearDowm(self):
         pass
 
 
+class TestDensityFeatureForCell(unittest.TestCase):
+    def test_cell(self):
+        n_included = 123
+        n_excluded = 456
+        x = np.append(np.zeros(n_included), np.ones(n_excluded))
+        environment = create_point_cloud(x, x, x)
+        target = create_point_cloud(np.zeros(1), np.zeros(1), np.zeros(1))
+
+        cell = Cell(1)  # area = 1.0
+        neighbors_index = compute_neighborhoods(environment, target, cell)
+
+        extractor = PointDensityFeatureExtractor()
+        for index in neighbors_index:
+            d = extractor.extract(environment, index, target, [0], cell)
+            self.assertEqual(d, n_included)
+
+
+class TestDensityFeatureForCube(unittest.TestCase):
+    def test_cell(self):
+        n_included = 123
+        n_excluded = 456
+        x = np.append(np.zeros(n_included), np.ones(n_excluded))
+        environment = create_point_cloud(x, x, x)
+        target = create_point_cloud(np.zeros(1), np.zeros(1), np.zeros(1))
+
+        cube = Cube(1)  # volume = 1.0
+        neighbors_index = compute_neighborhoods(environment, target, cube)
+
+        extractor = PointDensityFeatureExtractor()
+        for index in neighbors_index:
+            d = extractor.extract(environment, index, target, [0], cube)
+            self.assertEqual(d, n_included)
+
+
 class TestDensityFeatureOnRealData(unittest.TestCase):
     """Test density extractor on real data and make sure it doesn't crash."""
 
@@ -178,7 +205,7 @@ def _get_random_targets(self):
         num_all_pc_points = len(self.point_cloud[keys.point]["x"]["data"])
         rand_indices = [random.randint(0, num_all_pc_points)
                         for p in range(20)]
-        return utils.copy_pointcloud(self.point_cloud, rand_indices)
+        return utils.copy_point_cloud(self.point_cloud, rand_indices)
 
     def setUp(self):
         """Set up the test."""

laserchicken/feature_extractor/test_echo_ratio_extractor.py

@@ -83,7 +83,7 @@ def _set_cylinder_data(self):
 
     def _get_central_point(self, index):
         """Get the central point."""
-        return utils.copy_pointcloud(self.point_cloud, [index])
+        return utils.copy_point_cloud(self.point_cloud, [index])
 
     def setUp(self):
         """
@@ -165,7 +165,7 @@ def _get_random_targets(self):
         num_all_pc_points = len(self.point_cloud[keys.point]["x"]["data"])
         rand_indices = [random.randint(0, num_all_pc_points)
                         for p in range(20)]
-        return utils.copy_pointcloud(self.point_cloud, rand_indices)
+        return utils.copy_point_cloud(self.point_cloud, rand_indices)
 
 
 if __name__ == '__main__':

laserchicken/feature_extractor/test_eigenvals_feature_extractor.py

@@ -20,7 +20,7 @@ def test_eigenvalues_in_cylinders(self):
         num_all_pc_points = len(point_cloud[keys.point]["x"]["data"])
         rand_indices = [random.randint(0, num_all_pc_points)
                         for _ in range(20)]
-        target_point_cloud = utils.copy_pointcloud(point_cloud, rand_indices)
+        target_point_cloud = utils.copy_point_cloud(point_cloud, rand_indices)
         n_targets = len(target_point_cloud[keys.point]["x"]["data"])
         radius = 2.5
         neighbors = compute_neighbors.compute_cylinder_neighborhood(