Add load point cloud in IO

idefix/io.py

@@ -11,6 +11,7 @@ General functions to load and dump data in various format.
 import logging
 from pathlib import Path
 import numpy as np
+import numpy.core.records as rcd
 from numpy.lib import recfunctions as rfn
 import laspy
 
@@ -36,11 +37,7 @@ def load_las(fname):
 
         [(spatial), (feature, [f1, f2, ..., fn])]
     '''
-    fname = Path(fname)
-    if not fname.is_file():
-        msg = 'No such file: \'{}\''.format(fname)
-        log.error(msg)
-        raise IOError(msg)
+    fname = _get_verify_path(fname)
 
     log.info('Loading LAS file \'{}\'...'.format(fname))
     try:
@@ -106,11 +103,7 @@ def load_txt(fname, header, delimiter=' ', dtype=None):
 
         [(spatial), (feature, [f1, f2, ..., fn])]
     '''
-    fname = Path(fname)
-    if not fname.is_file():
-        msg = 'No such file: \'{}\''.format(fname)
-        log.error(msg)
-        raise IOError(msg)
+    fname = _get_verify_path(fname)
 
     if dtype is not None:
         assert len(dtype) == len(header), 'dtype and header must be the same size'
@@ -153,3 +146,50 @@ def load_txt(fname, header, delimiter=' ', dtype=None):
     pcloud = rfn.append_fields(spatial, 'feature', feature, usemask=False, asrecarray=True)
 
     return pcloud
+
+def _get_verify_path(fname):
+    fname = Path(fname)
+    if not fname.is_file():
+        msg = 'No such file: \'{}\''.format(fname)
+        log.error(msg)
+        raise IOError(msg)
+    return fname
+
+def _arr_to_rec(arr):
+    """Array to record array.
+
+    Used for point clouds, should work for everything else tho...
+    """
+    arrays = []; dtypes = []
+    for k in arr.dtype.fields.keys():
+        arrays += [arr[k]]
+        dtypes += [(k, arr.dtype[k])]
+    return np.core.records.fromarrays(arrays, dtypes)
+
+def load_pc(fname):
+    """Load point cloud from file.
+    
+    Loader for point clouds containted in '.pc', '.pcz' or compatible '.npy',
+    '.npz' files. This "point cloud" format is based on NumPy files, with small
+    overhead to manage record array and multispectral point clouds.
+
+    Parameters
+    ----------
+    fname : str, Path
+        Path to the point cloud file to load.
+
+    Returns
+    -------
+    point_cloud : recarray or tuple of recarray
+        The point cloud respecting or tuple of point clouds (for multispectral
+        point cloud files).
+    """
+    log.info('Loading point cloud file \'{}\')'.format(fname))
+
+    fname = _get_verify_path(fname)
+
+    archive = np.load(fname)
+    if len(archive.files) == 1:
+        return _arr_to_rec(archive[archive.files[0]])
+    else:
+        return tuple(_arr_to_rec(archive[arr]) for arr in archive.files)

idefix/vxl.py

@@ -12,7 +12,7 @@ import logging
 import numpy as np
 import humanize
 from .utils import bbox
-import ipdb
+import mayavi.mlab as mlab
 
 log = logging.getLogger(__name__)
 
@@ -270,6 +270,9 @@ def _geo_to_np_coordinate(raster):
     '''
     return np.flip(np.swapaxes(raster, 0, 1), 0)
 
+def _np_to_geo_coordinate(raster):
+    return np.swapaxes(np.flip(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()
@@ -342,3 +345,37 @@ def squash(voxel_grid, method='top', axis=-1):
     
     raise NotImplementedError('Method \'{}\' does not exist.'.format(method))
 
+def plot(voxel_grid, vmin=None, vmax=None):
+    """Plot voxel grid with Mayavi.
+
+    Parameters
+    ----------
+    voxel_grid : masked array (3D)
+        The voxel grid to plot.
+    vmin, vmax : scalar, optional
+        Define the data range that the colormap cover.
+
+    Returns
+    -------
+    figure : mlab figure
+        The figure instance.
+
+    Examples
+    --------
+    >>> a = np.random.random((10,10,10))
+    >>> view = {}
+    >>> mlab.clf()
+    >>> vxl.plot(a)
+    >>> mlab.view(**view)
+    >>> mlab.savefig(fname, magnification=4)
+    >>> mlab.show()
+    """
+    points = np.where(~voxel_grid.mask)
+
+    if vmin or vmax:
+        disp_value = np.clip(voxel_grid[points], vmin, vmax)
+    else:
+        disp_value = voxel_grid[points]
+
+    voxels_display = mlab.points3d(*points, disp_value, mode='cube', scale_factor=1, scale_mode='none', opacity=1., colormap='viridis')
+    return voxels_display

mlab_test.py

@@ -0,0 +1,22 @@
+#!/usr/bin/env python
+# file mlab_test.py
+# author Florent Guiotte <florent.guiotte@uhb.fr>
+# version 0.0
+# date 11 avril 2019
+"""Abstract
+
+doc.
+"""
+
+from idefix import vxl
+import mayavi.mlab as mlab
+import numpy as np
+
+spatial = np.random.random((10000, 3))
+feature = np.random.random(10000)
+
+grid = vxl.get_grid(spatial, .1)
+vg = vxl.bin(grid, spatial, feature, 'mean')
+
+vxl.plot(vg)
+mlab.show()

test/pytest.ini

@@ -0,0 +1,7 @@
+[pytest]
+filterwarnings = 
+    ignore::DeprecationWarning:sqlalchemy.*:
+    ignore::DeprecationWarning:apptools.*:
+    ignore::DeprecationWarning:pyface.*:
+    ignore::DeprecationWarning:traits.*:
+    ignore:.*escape sequence.*:DeprecationWarning

test/test_io.py

@@ -85,3 +85,30 @@ def test_load_txt(datadir, fname, head, separator, exp_point_count, exp_field_co
     if dtype is not None:
         for feature_name, feature_dtype in zip(head[3:], dtype[3:]):
             assert result.feature[feature_name].dtype == feature_dtype, "Missmatch between specified dtype and returned feature dtype"
+
+@pytest.mark.parametrize('fname, exp_point_count, exp_field_count', [
+    ('test.npz', 58629, 2, ),
+    ('test_compressed.npz', 58629, 2,),
+])
+def test_load_pc(datadir, fname, exp_point_count, exp_field_count):
+    fname = datadir.join(fname)
+    
+    # Raise "No such file"
+    with pytest.raises(IOError) as e_info:
+        io.load_pc('not_as_file.las')
+
+    # Open file without exception
+    try:
+        result = io.load_pc(fname)
+    except IOError:
+        pytest.fail('Opening legit file without exception')
+
+    assert result.size == exp_point_count, "Return correct point count"
+
+    assert result['spatial'].shape[-1] == 3, "Return ndarray with spatial field"
+
+    assert (result['spatial'] == result.spatial).all(), "Quick access with records array"
+
+    assert len(result['feature'].dtype) == exp_field_count, "Return ndarray with attribute fields"
+
+    assert result.spatial.dtype == np.float, "Dtype of spatial is np.float"