Fix coordinate system

doc/voxels.md

@@ -1,18 +1,55 @@
 Documentation and developer's notes on voxels in IDEFIX package.
 
+Main concepts
+=============
+
+3D Images
+---------
+
+### Image referential system
+
+Similarly to 2D images of go referenced data, the 3D images coordinates system
+does not directly correspond to the one contained in the LiDAR point clouds.
+
+When comparing NumPy and geographic coordinate systems for 2D raters we can
+notice these orientations:
+
+```
+NumPy       Geophic
+  y      
+ +-->         ^
+x|           y|
+ V            +-->
+                x
+ Coordinate systems
+```
+
+For 3D data, the NumPy coordinate respect the right hand rule while the LiDAR
+*z* axis will directly correspond to altitude, thus breaking the right hand
+rule.
+
+IDEFIX silently project the data from one system to another while providing
+utility function to get back any projection or axis labels.
+
 Voxels data
 ---
 
-- Density
+- Point cloud arrangement
-- Distribution
+    + **Density**
-    + Mean
+    + Orientation PCA
+    + Normals
+- Feature distribution
+    + **Mean**
     + Variance
-    + Mode
+    + **Mode**
     + Entropy
     + Min
     + Max
     + Quantil
 
+Annexes
+=======
+
 Voxels glossary
 ---
 

idefix/vxl.py

@@ -102,20 +102,21 @@ def bin(grid, spatial, feature=None, method='density'):
     log.info('Bining point cloud in grid...')
 
     if method == 'density':
-        return _bin_density(grid, spatial)
+        geo_rst = _bin_density(grid, spatial)
     else:
         if feature is None:
             msg = 'Missing required argument : \'feature\'.'
             log.error(msg)
             raise ValueError(msg)
         if method == 'mean':
-        return _bin_mean(grid, spatial, feature)
+            geo_rst = _bin_mean(grid, spatial, feature)
-    if method == 'mode':
+        elif method == 'mode':
-        return _bin_mode(grid, spatial, feature)
+            geo_rst = _bin_mode(grid, spatial, feature)
-
+        else:
             msg = 'Method \'{}\' does not exist.'.format(method)
             log.error(msg)
             raise NotImplementedError(msg)
+    return _geo_to_np_coordinate(geo_rst)
 
 def _bin_density(grid, spatial):
     '''Bin spatial in a grid, density method.
@@ -258,3 +259,13 @@ def _human_to_bytes(human_size):
         if human_size.endswith(k):
             return float(human_size.strip(k)) * 1024 ** v
     raise IOError('Did not understand size: \'{}\''.format(human_size))
+
+def _geo_to_np_coordinate(raster):
+    '''Geographic to numpy coordinate system.
+    
+    Transfer the raster (2D and 3D) from a geographic coordinate system to the
+    numpy coordinate system.
+    '''
+    return np.flip(np.swapaxes(raster, 0, 1), 0)
+
+

test/test_vxl.py

@@ -157,3 +157,30 @@ def test_insight(method):
         vxl.insight(grid, method=method, mem_limit='3KB')
     with pytest.raises(MemoryError) as e_info:
         vxl.insight(grid, method=method, mem_limit=3000)
+
+def test__geo_to_np_coordinate():
+    raster = np.zeros((5, 5), dtype=np.uint8)
+    raster[0, 0] = 42
+    raster[4, 4] = 25
+    raster[2, 0] = 7
+
+    raster_truth = np.zeros_like(raster)
+    raster_truth[-1, 0] = 42
+    raster_truth[0, -1] = 25
+    raster_truth[-1, 2] = 7
+
+    assert (raster_truth == vxl._geo_to_np_coordinate(raster)).all(), 'Missmatch between 2D raters' 
+
+    raster = np.zeros((5, 5, 3), dtype=np.uint8)
+    raster[0, 0, 0] = 42
+    raster[4, 4, 1] = 25
+    raster[2, 0, 2] = 7
+
+    raster_truth = np.zeros_like(raster)
+    raster_truth[-1, 0, 0] = 42
+    raster_truth[0, -1, 1] = 25
+    raster_truth[-1, 2, 2] = 7
+
+    assert (raster_truth == vxl._geo_to_np_coordinate(raster)).all(), 'Missmatch between 3D raters' 
+
+

test/test_vxl/pc0_vxl_s1.txt

@@ -1,5 +1,5 @@
 # x y z density mean mode
-0 0 0 1 2 2
+9 0 0 1 2 2
-0 2 1 4 10 5
+7 0 1 4 10 5
-9 9 9 1 1 1
+0 9 9 1 1 1
-4 4 4 1 0 0
+5 4 4 1 0 0