Add squash function.

doc/voxels.md

@@ -39,11 +39,11 @@ Voxels data
     + Orientation PCA
     + Normals
 - Feature distribution
-    + **Mean**
-    + Variance
-    + **Mode**
+    + **Mean**: Reduce noise in the cell
+    + Variance: Characterize the distribution of the data
+    + **Mode**: Return the majority value (e.g. labels)
     + Entropy
-    + Min
+    + *Min*: To create last echoes map
     + Max
     + Quantil
 

idefix/vxl.py

@@ -262,10 +262,81 @@ def _human_to_bytes(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)
 
+def _squash_position(voxel_grid, method, axis):
+    squash_mask = np.zeros_like(voxel_grid, dtype=np.int)
+    mask_idx = (~voxel_grid.mask).nonzero()
+    squash_mask[mask_idx] = mask_idx[axis]
+
+    if method == 'top':
+        squash_id = squash_mask.max(axis=axis).astype(np.uint)
+    elif method == 'center':
+        squash_id = np.ma.median(squash_mask, axis=axis).astype(np.uint)
+    elif method == 'bottom':
+        squash_id = squash_mask.min(axis=axis).astype(np.uint)
+
+    xy_where = np.nonzero(~squash_id.mask)
+    voxel_grid_where = list(xy_where)
+    voxel_grid_where.insert(axis%(len(voxel_grid_where)+1), squash_id.compressed())
+
+    raster = np.zeros_like(squash_id)
+    raster[xy_where] = voxel_grid[tuple(voxel_grid_where)]
+
+    return raster
+
+def squash(voxel_grid, method='top', axis=-1):
+    """Flatten a voxel grid.
+
+    Squash the voxel grid along `axis` according to `method` into a raster.
+
+    The squash methods proposed are :
+
+    - Position based in the "column" (i.e. along axis).
+        + 'top': The first non empty cells (from top) is returned.
+        + 'center': The most centered cell is returned.
+        + 'bottom': The last
+    - Cell description in the "column".
+        + 'count': The number of non empty cells.
+        + 'mean': The mean value of the non empty cells.
+        + 'median': The median value...
+        + 'std': ...
+        + 'min': ...
+        + 'max': ...
+
+    Parameters
+    ----------
+    voxel_grid : masked array (3D)
+        The voxel grid (binned point cloud) to squash.
+    method : str
+        The squash method. It can be 'top', 'center', 'bottom', 'count', 'min',
+        'mean', 'max', 'std' or 'median'. Default is 'top'.
+    axis : number
+        The axis to squash along. Default is last (i.e. 2 for 3D voxel grid).
+
+    Return
+    ------
+    raster_grid : masked array (2D)
+        The squashed raster.
+    """
+    if method in ('top', 'center', 'bottom'):
+        return _squash_position(voxel_grid, method, axis)
+    elif method == 'count':
+        return ~voxel_grid.mask.sum(axis=axis)
+    elif method == 'mean':
+        return voxel_grid.mean(axis=axis)
+    elif method == 'median':
+        return np.ma.median(voxel_grid, axis=axis)
+    elif method == 'min':
+        return voxel_grid.min(axis=axis)
+    elif method == 'max':
+        return voxel_grid.max(axis=axis)
+    elif method == 'std':
+        return voxel_grid.std(axis=axis)
+    
+    raise NotImplementedError('Method \'{}\' does not exist.'.format(method))
 

test/test_vxl.py

@@ -39,7 +39,7 @@ def data_vxl(datadir, set_id, grid_id, method):
 
         i = fields.index(feature_name)
 
-        data = np.loadtxt('test/test_vxl/pc0_vxl_s1.txt')
+        data = np.loadtxt(fname)
         spatial = data[:,:3].astype(np.intp)
         feature = data[:,i]
 
@@ -48,6 +48,14 @@ def data_vxl(datadir, set_id, grid_id, method):
     path = datadir.join('pc{}_vxl_s{}.txt'.format(set_id, grid_id))
     return _load_vxl(path, method)
 
+def data_raster(datadir, set_id, grid_id, axis, method):
+    def _load_raster(fname):
+        data = np.loadtxt(fname)
+        return np.ma.masked_array(data, data==0)
+
+    path = datadir.join('pc{}_vxl_s{}_raster_{}_{}.txt'.format(set_id, grid_id, axis, method))
+    return _load_raster(path)
+
 def data_grid(datadir, set_id, step_id):
     def _read(fname):
         with open(fname, 'r') as f:
@@ -183,4 +191,26 @@ def test__geo_to_np_coordinate():
 
     assert (raster_truth == vxl._geo_to_np_coordinate(raster)).all(), 'Missmatch between 3D raters' 
 
+@pytest.mark.parametrize('set_id, grid_id, axis, method', [
+    (1, 1, 2, 'top'),
+    (1, 1, 2, 'center'),
+    (1, 1, 2, 'bottom'),
+    (1, 1, 2, 'mean'),
+    (1, 1, 2, 'max'),
+    (1, 1, 2, 'min'),
+    (1, 1, 2, 'median'),
+    (1, 1, 0, 'top'),
+    (1, 1, 0, 'center'),
+    (1, 1, 0, 'bottom'),
+    (1, 1, 1, 'top'),
+    (1, 1, 1, 'center'),
+    (1, 1, 1, 'bottom'),
+])
+def test_squash(datadir, set_id, grid_id, axis, method):
+    vxld = data_vxl(datadir, set_id, grid_id, 'density' )
+    truth = data_raster(datadir, set_id, grid_id, axis, method)
+    res = vxl.squash(vxld, method, axis)
 
+    assert res is not None, 'Tested function did not return anything :('
+    assert res.shape == truth.shape, 'Missmatch between truth and tested shape'
+    assert np.allclose(res, truth), 'Missmatch between truth and tested raster'

test/test_vxl/pc1_vxl_s1.txt

@@ -0,0 +1,9 @@
+# x y z density mean mode
+1 2 0 2 0 0
+1 3 0 8 0 0
+3 0 1 7 0 0
+3 4 1 42 0 0
+0 0 2 1 0 0
+1 2 2 4 0 0
+1 2 3 5 0 0
+2 2 3 6 0 0

test/test_vxl/pc1_vxl_s1_raster_0_bottom.txt

@@ -0,0 +1,5 @@
+0 7 1 0
+0 0 0 0
+2 0 4 5
+8 0 0 0
+0 42 0 0

test/test_vxl/pc1_vxl_s1_raster_0_center.txt

@@ -0,0 +1,5 @@
+0 7 1 0
+0 0 0 0
+2 0 4 5
+8 0 0 0
+0 42 0 0

test/test_vxl/pc1_vxl_s1_raster_0_top.txt

@@ -0,0 +1,5 @@
+0 7 1 0
+0 0 0 0
+2 0 4 6
+8 0 0 0
+0 42 0 0

test/test_vxl/pc1_vxl_s1_raster_1_bottom.txt

@@ -0,0 +1,4 @@
+0 0 1 0
+0 0 4 5
+0 0 0 6
+0 7 0 0

test/test_vxl/pc1_vxl_s1_raster_1_center.txt

@@ -0,0 +1,4 @@
+0 0 1 0
+0 0 4 5
+0 0 0 6
+0 7 0 0

test/test_vxl/pc1_vxl_s1_raster_1_top.txt

@@ -0,0 +1,4 @@
+0 0 1 0
+0 0 4 5
+0 0 0 6
+0 42 0 0

test/test_vxl/pc1_vxl_s1_raster_2_bottom.txt

@@ -0,0 +1,4 @@
+1 0 0 0 0
+0 0 2 8 0
+0 0 6 0 0
+7 0 0 0 42

test/test_vxl/pc1_vxl_s1_raster_2_center.txt

@@ -0,0 +1,4 @@
+1 0 0 0 0
+0 0 4 8 0
+0 0 6 0 0
+7 0 0 0 42

test/test_vxl/pc1_vxl_s1_raster_2_max.txt

@@ -0,0 +1,4 @@
+1 0 0 0 0
+0 0 5 8 0
+0 0 6 0 0
+7 0 0 0 42

test/test_vxl/pc1_vxl_s1_raster_2_mean.txt

@@ -0,0 +1,4 @@
+1 0 0 0 0
+0 0 3.66667 8 0
+0 0 6 0 0
+7 0 0 0 42

test/test_vxl/pc1_vxl_s1_raster_2_median.txt

@@ -0,0 +1,4 @@
+1 0 0 0 0
+0 0 4 8 0
+0 0 6 0 0
+7 0 0 0 42

test/test_vxl/pc1_vxl_s1_raster_2_min.txt

@@ -0,0 +1,4 @@
+1 0 0 0 0
+0 0 2 8 0
+0 0 6 0 0
+7 0 0 0 42

test/test_vxl/pc1_vxl_s1_raster_2_top.txt

@@ -0,0 +1,4 @@
+1 0 0 0 0
+0 0 5 8 0
+0 0 6 0 0
+7 0 0 0 42