Merge pull request #66 from eEcoLiDAR/feature_eigenvals

Feature eigenvals

laserchicken/compute_neighbors.py

@@ -15,17 +15,9 @@ def compute_cylinder_neighbourhood_indicies(env_pc, target_pc, radius):
 
     return neighb_points_indices
 
-def compute_cylinder_neighbourhoods(env_pc, target_pc, radius):
-    ''' Function to create a cylindrincal neighbourhood for a given point of
-    a target point cloud among the points from an environment point cloud'''
-
-    neighb_points_indices = compute_cylinder_neighbourhood_indicies(env_pc, target_pc, radius)
-
-    return [utils.copy_pointcloud(env_pc,iarray) for iarray in neighb_points_indices]
-
-def compute_sphere_neighbourhoods(env_pc, target_pc, radius):
-    ''' Function to create a spherical neighbourhood for a given point of
-    a target point cloud among the points from an environment point cloud'''
+def compute_sphere_neighbourhood_indicies(env_pc, target_pc, radius):
+    ''' Function to find the indicies of points within a spherical neighbourhood
+    for a given point of a target point cloud among the points from an environment point cloud'''
 
     neighb_points_indices = compute_cylinder_neighbourhood_indicies(env_pc, target_pc, radius)
 
@@ -39,9 +31,21 @@ def compute_sphere_neighbourhoods(env_pc, target_pc, radius):
             if(abs(targetz - envz) > radius): continue
             if((envx - targetx)**2 + (envy - targety)**2 + (envz - targetz)**2 <= radius**2):
                 resultindices.append(j)
-        result.append(utils.copy_pointcloud(env_pc,resultindices))
+        result.append(resultindices)
     return result
 
+def compute_cylinder_neighbourhoods(env_pc, target_pc, radius):
+    ''' Function to create a cylindrincal neighbourhood for a given point of
+    a target point cloud among the points from an environment point cloud'''
+
+    neighb_points_indices = compute_cylinder_neighbourhood_indicies(env_pc, target_pc, radius)
+
+    return [utils.copy_pointcloud(env_pc,iarray) for iarray in neighb_points_indices]
+
+def compute_sphere_neighbourhoods(env_pc, target_pc, radius):
+    ''' Function to create a spherical neighbourhood for a given point of
+    a target point cloud among the points from an environment point cloud'''
 
+    neighb_points_indices = compute_sphere_neighbourhood_indicies(env_pc, target_pc, radius)
 
-    return #[utils.copy_pointcloud(env_pc,iarray) for iarray in neighb_points_indices]
+    return [utils.copy_pointcloud(env_pc,iarray) for iarray in neighb_points_indices]

laserchicken/feature_extractor/__init__.py

@@ -2,11 +2,11 @@
 import importlib
 import re
 import numpy as np
-from laserchicken import keys,utils
-#from .gijs_elena_feature import MyFeatureExtractor
+from laserchicken import keys, utils
+from .eigenvals_feature_extractor import EigenValueFeatureExtractor
 
 
-def _feature_map(module_name = __name__):
+def _feature_map(module_name=__name__):
     """Construct a mapping from feature names to feature extractor classes."""
     module = importlib.import_module(module_name)
     return {
@@ -19,29 +19,30 @@ def _feature_map(module_name = __name__):
 FEATURES = _feature_map()
 
 
-def compute_features(env_point_cloud, neighborhoods, target_point_cloud, feature_names, overwrite = False):
+def compute_features(env_point_cloud, neighborhoods, target_point_cloud, feature_names, overwrite=False):
     ordered_features = _make_feature_list(feature_names)
     targetsize = len(target_point_cloud[keys.point]["x"]["data"])
     for feature in ordered_features:
-        if((not overwrite) and (feature in target_point_cloud[keys.point])):
-            continue # Skip feature calc if it is already there and we do not overwrite
+        if ((not overwrite) and (feature in target_point_cloud[keys.point])):
+            continue  # Skip feature calc if it is already there and we do not overwrite
         extractor = FEATURES[feature]()
         providedfeatures = extractor.provides()
         numfeatures = len(providedfeatures)
-        featurevalues = [np.empty(targetsize,dtype = np.float64) for i in range(numfeatures)]
+        featurevalues = [np.empty(targetsize, dtype=np.float64) for i in range(numfeatures)]
         for target_index in range(targetsize):
-            pointvalues = extractor.extract(env_point_cloud, neighborhoods[target_index], target_point_cloud, target_index)
-            if(numfeatures > 1):
+            pointvalues = extractor.extract(env_point_cloud, neighborhoods[target_index], target_point_cloud,
+                                            target_index)
+            if numfeatures > 1:
                 for i in range(numfeatures):
                     featurevalues[i][target_index] = pointvalues[i]
             else:
                 featurevalues[0][target_index] = pointvalues
         for i in range(numfeatures):
             fname = providedfeatures[i]
-            if(overwrite or (not fname in target_point_cloud[keys.point])):
+            if (overwrite or (fname not in target_point_cloud[keys.point])):
                 # Set feature values if it is not there (or we want to overwrite)
-                target_point_cloud[keys.point][fname] = {"type" : np.float64, "data" : featurevalues[i]}
-        utils.add_metadata(target_point_cloud,type(extractor).__module__,extractor.get_params())
+                target_point_cloud[keys.point][fname] = {"type": np.float64, "data": featurevalues[i]}
+        utils.add_metadata(target_point_cloud, type(extractor).__module__, extractor.get_params())
 
 
 def _make_feature_list(feature_names):

laserchicken/feature_extractor/eigenvals_feature_extractor.py

@@ -0,0 +1,47 @@
+import numpy as np
+
+from laserchicken import utils
+from laserchicken.feature_extractor.abc import AbstractFeatureExtractor
+
+
+def structure_tensor(points):
+    """
+    Computes the structure tensor of points by computing the eigenvalues
+    and eigenvectors of the covariance matrix of a point cloud.
+    Parameters
+    ----------
+    points : (Mx3) array
+        X, Y and Z coordinates of points.
+    Returns
+    -------
+    eigenvalues : (1x3) array
+        The eigenvalues corresponding to the eigenvectors of the covariance
+        matrix.
+    eigenvectors : (3,3) array
+        The eigenvectors of the covariance matrix.
+    """
+    if points.shape[0] > 3:
+        cov_mat = np.cov(points, rowvar=False)
+        eigenvalues, eigenvectors = np.linalg.eig(cov_mat)
+        order = np.argsort(-eigenvalues)
+        eigenvalues = eigenvalues[order]
+        eigenvectors = eigenvectors[:, order]
+        return eigenvalues, eigenvectors
+    else:
+        raise ValueError('Not enough points to compute eigenvalues/vectors.')
+
+
+class EigenValueFeatureExtractor(AbstractFeatureExtractor):
+    @classmethod
+    def requires(cls):
+        return []
+
+    @classmethod
+    def provides(cls):
+        return ['eigenv_1', 'eigenv_2', 'eigenv_3']
+
+    def extract(self, sourcepc, neighborhood, targetpc, targetindex):
+        nbptsX, nbptsY, nbptsZ = utils.get_point(sourcepc, neighborhood)
+        matrix = np.column_stack((nbptsX, nbptsY, nbptsZ))
+        eigenvals, eigenvecs = structure_tensor(matrix)
+        return [eigenvals[0], eigenvals[1], eigenvals[2]]

laserchicken/test_feature_extractor/test_eigenvals_feature_extractor.py

@@ -0,0 +1,36 @@
+import os
+import random
+import unittest
+
+from laserchicken import compute_neighbors, feature_extractor, keys, read_las, utils
+
+
+class TestExtractEigenValues(unittest.TestCase):
+
+    _test_file_name = 'AHN3.las'
+    _test_data_source = 'testdata'
+    pointcloud = None
+
+    def test_eigenvalues_in_cylinders(self):
+        """Test computing of eigenvalues in cylinder."""
+        num_all_pc_points = len(self.pointcloud[keys.point]["x"]["data"])
+        rand_indices = [random.randint(0, num_all_pc_points) for p in range(20)]
+        target_pointcloud = utils.copy_pointcloud(self.pointcloud, rand_indices)
+        numtargets = len(target_pointcloud[keys.point]["x"]["data"])
+        radius = 2.5
+        result_index_lists = compute_neighbors.compute_cylinder_neighbourhood_indicies(
+            self.pointcloud, target_pointcloud, radius)
+        feature_extractor.compute_features(self.pointcloud, result_index_lists, target_pointcloud,
+                                           ["eigenv_1", "eigenv_2", "eigenv_3"])
+        for i in range(numtargets):
+            lambda1, lambda2, lambda3 = utils.get_features(target_pointcloud, i, ["eigenv_1", "eigenv_2", "eigenv_3"])
+            self.assertTrue(lambda1 >= lambda2 and lambda2 >= lambda3)
+        self.assertEqual("laserchicken.feature_extractor.eigenvals_feature_extractor",
+                         target_pointcloud[keys.provenance][0]["module"])
+
+    def setUp(self):
+        self.pointcloud = read_las.read(os.path.join(self._test_data_source, self._test_file_name))
+        random.seed(102938482634)
+
+    def tearDown(self):
+        pass

laserchicken/test_feature_extractor/test_extract_features.py

@@ -1,42 +1,51 @@
 """Test feature extraction."""
-import pytest
 import numpy as np
-from laserchicken import feature_extractor,keys,test_tools
+import pytest
+
+from laserchicken import feature_extractor, keys, test_tools
 
 from . import __name__ as test_module_name
 
-# Overwrite the available feature extractors with test feature extractors
-feature_extractor.FEATURES = feature_extractor._feature_map(test_module_name)
+
+@pytest.fixture(scope='module', autouse=True)
+def override_features():
+    """Overwrite the available feature extractors with test feature extractors."""
+    feature_extractor.FEATURES = feature_extractor._feature_map(test_module_name)
+    yield
+    feature_extractor.FEATURES = feature_extractor._feature_map(feature_extractor.__name__)
 
 
-def _compute_features(target,feature_names,overwrite = False):
+def _compute_features(target, feature_names, overwrite=False):
     neighborhoods = [[] for i in range(len(target["vertex"]["x"]["data"]))]
     feature_extractor.compute_features({}, neighborhoods, target, feature_names, overwrite)
     return target
 
 
 def test_extract_single_feature():
     target = test_tools.ComplexTestData.get_point_cloud()
-    _compute_features(target,['test3_a'])
-    assert ('test1_b' in target[keys.point])
+    _compute_features(target, ['test3_a'])
+    assert 'test1_b' in target[keys.point]
     assert all(target[keys.point]['test3_a']['data'] == target[keys.point]['z']['data'])
 
+
 def test_extract_multiple_features():
     target = test_tools.ComplexTestData.get_point_cloud()
-    feature_names = ['test3_a','test2_b']
-    target = _compute_features(target,feature_names)
+    feature_names = ['test3_a', 'test2_b']
+    target = _compute_features(target, feature_names)
     assert ('test3_a' in target[keys.point] and 'test2_b' in target[keys.point])
 
+
 def test_extract_does_not_overwrite():
     target = test_tools.ComplexTestData.get_point_cloud()
-    target[keys.point]['test2_b'] = {"type":np.float64,"data":[0.9,0.99,0.999,0.9999]}
-    feature_names = ['test3_a','test2_b']
-    target = _compute_features(target,feature_names)
-    assert (target[keys.point]['test2_b']['data'][2] == 0.999)
+    target[keys.point]['test2_b'] = {"type": np.float64, "data": [0.9, 0.99, 0.999, 0.9999]}
+    feature_names = ['test3_a', 'test2_b']
+    target = _compute_features(target, feature_names)
+    assert target[keys.point]['test2_b']['data'][2] == 0.999
+
 
 def test_extract_can_overwrite():
     target = test_tools.ComplexTestData.get_point_cloud()
-    target[keys.point]['test2_b'] = {"type":np.float64,"data":[0.9,0.99,0.999,0.9999]}
-    feature_names = ['test3_a','test2_b']
-    target = _compute_features(target,feature_names,overwrite = True)
-    assert (target[keys.point]['test2_b']['data'][2] == 11.5)
+    target[keys.point]['test2_b'] = {"type": np.float64, "data": [0.9, 0.99, 0.999, 0.9999]}
+    feature_names = ['test3_a', 'test2_b']
+    target = _compute_features(target, feature_names, overwrite=True)
+    assert target[keys.point]['test2_b']['data'][2] == 11.5

laserchicken/test_feature_extractor/test_feature_map.py

@@ -1,10 +1,20 @@
 """Test that the map from feature names to extractor classes is correct."""
-from laserchicken.feature_extractor import _feature_map
+import pytest
+
+from laserchicken import feature_extractor
 
 from . import __name__ as test_module_name
 from . import Test1FeatureExtractor, Test2FeatureExtractor, Test3FeatureExtractor, TestBrokenFeatureExtractor
 
 
+@pytest.fixture(scope='module', autouse=True)
+def override_features():
+    """Overwrite the available feature extractors with test feature extractors."""
+    feature_extractor.FEATURES = feature_extractor._feature_map(test_module_name)
+    yield
+    feature_extractor.FEATURES = feature_extractor._feature_map(feature_extractor.__name__)
+
+
 def test__feature_map():
     feature_map = {
         'test1_a': Test1FeatureExtractor,
@@ -15,4 +25,4 @@ def test__feature_map():
         'test3_a': Test3FeatureExtractor,
         'test_broken': TestBrokenFeatureExtractor,
     }
-    assert feature_map == _feature_map(test_module_name)
+    assert feature_map == feature_extractor._feature_map(test_module_name)