diff --git a/doc/index.md b/doc/index.md
index 65c1ac3..48dc204 100644
--- a/doc/index.md
+++ b/doc/index.md
@@ -10,7 +10,7 @@ Welcome to Idefix's documentation!
    :maxdepth: 2
    :caption: Contents:
 
-
+.. _voxels
 
 Indices and tables
 ==================
diff --git a/doc/voxels.md b/doc/voxels.md
new file mode 100644
index 0000000..0f4230f
--- /dev/null
+++ b/doc/voxels.md
@@ -0,0 +1,40 @@
+Documentation and developer's notes on voxels in IDEFIX package.
+
+Voxels data
+---
+
+- Density
+- Distribution
+    + Mean
+    + Variance
+    + Mode
+    + Entropy
+    + Min
+    + Max
+    + Quantil
+
+Voxels glossary
+---
+
+- Step, resolution, bins...
+
+Algorithm
+---
+
+### Mode
+
+Insight and warn when use with high entropic data.
+
+- Slow (usual mode)
+- Quick (versus boost)
+
+Ideas
+---
+
+- Unified voxel format
+    + Do not store complete matrix
+        * Point cloud like to benefit of PC IO for import/export.
+        * Matrices always start to 0,0,0 image like, need metadata (start,
+          steps, dtypes, field names...)
+    + Allow matrice aumgentation
+
diff --git a/idefix/__init__.py b/idefix/__init__.py
index 3a51d50..96b59e1 100644
--- a/idefix/__init__.py
+++ b/idefix/__init__.py
@@ -9,7 +9,6 @@ Utils and production pipelines for processing LiDAR point clouds.
 
 """
 
-__all__ = ['utils', 'io']
+__all__ = ['utils', 'io', 'vxl']
 
-from . import utils
-from . import io
+from . import utils, io, vxl
diff --git a/idefix/io.py b/idefix/io.py
index 099f1d8..0841f0f 100644
--- a/idefix/io.py
+++ b/idefix/io.py
@@ -3,9 +3,9 @@
 # author Florent Guiotte <florent.guiotte@uhb.fr>
 # version 0.0
 # date 04 mars 2019
-"""Abstract
+""" IO IDEFIX sub-package
 
-doc.
+General functions to load and dump data in various format.
 """
 
 import logging
@@ -14,7 +14,7 @@ import numpy as np
 from numpy.lib import recfunctions as rfn
 import laspy
 
-log = logging.getLogger(__name__) 
+log = logging.getLogger(__name__)
 
 def load_las(fname):
     '''Load a LAS file into idefix point cloud format.
@@ -51,7 +51,7 @@ def load_las(fname):
         raise IOError(msg)
 
     log.debug('Extract spatial data')
-    spatial = np.core.records.fromarrays([np.array((infile.x, infile.y, infile.z)).T], 
+    spatial = np.core.records.fromarrays([np.array((infile.x, infile.y, infile.z)).T],
                                          dtype=[('spatial', np.float, 3)])
 
     log.debug('Extract feature data')
@@ -95,7 +95,7 @@ def load_txt(fname, header, delimiter=' '):
         Names of the columns contained in the text point cloud file.
     delimiter : str, optional
         String used to separate values. The default is whitespace.
-        
+
     Returns
     -------
     pcloud : recarray
@@ -128,7 +128,7 @@ def load_txt(fname, header, delimiter=' '):
     raw_txt = np.loadtxt(fname, delimiter=delimiter, dtype=dtype)
 
     log.debug('Extract spatial data')
-    spatial = np.core.records.fromarrays([np.array([raw_txt[x] for x in ('x', 'y', 'z')]).T], 
+    spatial = np.core.records.fromarrays([np.array([raw_txt[x] for x in ('x', 'y', 'z')]).T],
                                          dtype=[('spatial', np.float, 3)])
 
     log.debug('Extract feature data')
diff --git a/idefix/vxl.py b/idefix/vxl.py
new file mode 100644
index 0000000..3d99c72
--- /dev/null
+++ b/idefix/vxl.py
@@ -0,0 +1,60 @@
+#!/usr/bin/env python
+# file vxl.py
+# author Florent Guiotte <florent.guiotte@uhb.fr>
+# version 0.0
+# date 19 mars 2019
+""" Voxel IDEFIX sub-package
+
+General functions to transform point clouds to voxels compatible with numpy.
+"""
+
+import logging
+import numpy as np
+from .utils import bbox
+import ipdb
+
+log = logging.getLogger(__name__)
+
+def _ui_step(step):
+    '''User input management for step (number or array)
+    '''
+    try:
+        iter(step)
+        if len(step) != 3:
+            msg = 'Wrong steps input, 3 steps expected in step = \'{}\''.format(step)
+            log.error(msg)
+            raise IOError(msg)
+    except TypeError:
+        step = [step] * 3
+    return step
+
+def get_grid(spatial, step):
+    '''Return grid bins.
+
+    Compute the grid bins of a spatial point cloud or corresponding bounding
+    box according to given step (or steps for anisotropic grid).
+
+    Parameters
+    ----------
+    spatial : array (n, 3)
+        The spatial point cloud or the corresponding bounding box to grid.
+    step : number or array or tuple
+        The step of the grid, can be a number to get an isotropic grid, or an
+        iterable of size 3 (required) to get an anisotropic grid.
+
+    Returns
+    -------
+    grid : array of array (3,)
+        Grid of spatial given step. Return three arrays (not necessarily of the
+        same size) defining the bins of axis `x`, `y` and `z`.
+    '''
+    bb = bbox(spatial)
+    step = _ui_step(step)
+
+    #ipdb.set_trace()
+    grid = []
+    for start, stop, s in zip(bb[0], bb[1], step):
+        grid += [np.arange(start, stop + 2*s, s)]
+
+    return grid
+
diff --git a/test/test_vxl.py b/test/test_vxl.py
new file mode 100644
index 0000000..e3a6ed6
--- /dev/null
+++ b/test/test_vxl.py
@@ -0,0 +1,75 @@
+#!/usr/bin/env python
+# file test_vxl.py
+# author Florent Guiotte <florent.guiotte@uhb.fr>
+# version 0.0
+# date 18 mars 2019
+""" Test functions for the voxels `vxl` subpackage.
+
+"""
+
+from dataclasses import dataclass
+import pytest
+import numpy as np
+from idefix import vxl
+
+@dataclass
+class Pcloud:
+    spatial: np.ndarray
+    feature: np.ndarray
+
+def data_pc(datadir, set_id):
+    path = datadir.join('pc{}.txt'.format(set_id))
+    data = np.loadtxt(path)
+    return Pcloud(data[:,:3], data[:,3])
+
+#def data_vxl(datadir, set_id, step, method):
+#    pass
+
+@pytest.fixture
+def data_0_vxl():
+    def _data_0_vxl(method, resolution):
+        if method == 'mean':
+            pass
+
+def data_grid(datadir, set_id, step_id):
+    def _read(fname):
+        with open(fname, 'r') as f:
+            grid = []
+            for arr in f.readlines():
+                grid += [np.array([float(x) for x in arr.split(' ')])]
+        return grid
+
+    path = datadir.join('pc{}_grid_s{}.txt'.format(set_id, step_id))
+    return _read(path)
+    
+@pytest.mark.parametrize('set_id, step, grid_id', [
+    ('0', 1., '1'),
+    ('0', .1, '0_1'),
+])
+def test_get_grid(datadir, set_id, step, grid_id):
+    spatial = data_pc(datadir, set_id).spatial
+    res = data_grid(datadir, set_id, grid_id)
+
+    assert spatial is not None, 'Test data empty, test function is broken!'
+    assert spatial.shape[-1] == 3, 'Test data malformed, test function is broken!'
+    assert res is not None, 'Test data empty, test function is broken!'
+    assert len(res) == 3, 'Test data malformed, test function is broken!'
+
+    test = vxl.get_grid(spatial, step)
+
+    assert test is not None, 'Function did not return anything :('
+    assert len(test) == 3, 'Function doesn\'t give right number of axis'
+
+    for axis_test, axis_truth in zip(test, res):
+        assert axis_test.size == axis_truth.size, 'Wrong size for axis'
+        assert (axis_test == axis_truth).all(), 'Axis inequality between truth and test'
+                         
+def test_grid():
+    """
+    - dtype
+    - method
+    - mask
+    - data
+    """
+    pass
+
diff --git a/test/test_vxl/pc0.txt b/test/test_vxl/pc0.txt
new file mode 100644
index 0000000..b5d64b8
--- /dev/null
+++ b/test/test_vxl/pc0.txt
@@ -0,0 +1,8 @@
+# x y z feature
+1 1 1 2
+1 3 2 5
+1 3 2 5
+1 3 2 10
+1 3 2 20
+10 10 10 1
+5 5 5 0
diff --git a/test/test_vxl/pc0_grid_s0_1.txt b/test/test_vxl/pc0_grid_s0_1.txt
new file mode 100644
index 0000000..1dcf6e5
--- /dev/null
+++ b/test/test_vxl/pc0_grid_s0_1.txt
@@ -0,0 +1,3 @@
+1.0 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 1.9 2.0 2.1 2.2 2.3 2.4 2.5 2.6 2.7 2.8 2.9 3.0 3.1 3.2 3.3 3.4 3.5 3.6 3.7 3.8 3.9 4.0 4.1 4.2 4.3 4.4 4.5 4.6 4.7 4.8 4.9 5.0 5.1 5.2 5.3 5.4 5.5 5.6 5.7 5.8 5.9 6.0 6.1 6.2 6.3 6.4 6.5 6.6 6.7 6.8 6.9 7.0 7.1 7.2 7.3 7.4 7.5 7.6 7.7 7.8 7.9 8.0 8.1 8.2 8.3 8.4 8.5 8.6 8.7 8.8 8.9 9.0 9.1 9.2 9.3 9.4 9.5 9.6 9.7 9.8 9.9 10.0 10.1
+1.0 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 1.9 2.0 2.1 2.2 2.3 2.4 2.5 2.6 2.7 2.8 2.9 3.0 3.1 3.2 3.3 3.4 3.5 3.6 3.7 3.8 3.9 4.0 4.1 4.2 4.3 4.4 4.5 4.6 4.7 4.8 4.9 5.0 5.1 5.2 5.3 5.4 5.5 5.6 5.7 5.8 5.9 6.0 6.1 6.2 6.3 6.4 6.5 6.6 6.7 6.8 6.9 7.0 7.1 7.2 7.3 7.4 7.5 7.6 7.7 7.8 7.9 8.0 8.1 8.2 8.3 8.4 8.5 8.6 8.7 8.8 8.9 9.0 9.1 9.2 9.3 9.4 9.5 9.6 9.7 9.8 9.9 10.0 10.1
+1.0 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 1.9 2.0 2.1 2.2 2.3 2.4 2.5 2.6 2.7 2.8 2.9 3.0 3.1 3.2 3.3 3.4 3.5 3.6 3.7 3.8 3.9 4.0 4.1 4.2 4.3 4.4 4.5 4.6 4.7 4.8 4.9 5.0 5.1 5.2 5.3 5.4 5.5 5.6 5.7 5.8 5.9 6.0 6.1 6.2 6.3 6.4 6.5 6.6 6.7 6.8 6.9 7.0 7.1 7.2 7.3 7.4 7.5 7.6 7.7 7.8 7.9 8.0 8.1 8.2 8.3 8.4 8.5 8.6 8.7 8.8 8.9 9.0 9.1 9.2 9.3 9.4 9.5 9.6 9.7 9.8 9.9 10.0 10.1
diff --git a/test/test_vxl/pc0_grid_s0_15.txt b/test/test_vxl/pc0_grid_s0_15.txt
new file mode 100644
index 0000000..1ad8874
--- /dev/null
+++ b/test/test_vxl/pc0_grid_s0_15.txt
@@ -0,0 +1,3 @@
+1.0 1.15 1.3 1.45 1.6 1.75 1.9 2.05 2.2 2.35 2.5 2.65 2.8 2.95 3.1 3.25 3.4 3.55 3.7 3.85 4.0 4.15 4.3 4.45 4.6 4.75 4.9 5.05 5.2 5.35 5.5 5.65 5.8 5.95 6.1 6.25 6.4 6.55 6.7 6.85 7.0 7.15 7.3 7.45 7.6 7.75 7.9 8.05 8.2 8.35 8.5 8.65 8.8 8.95 9.1 9.25 9.4 9.55 9.7 9.85 10.0 10.15
+1.0 1.15 1.3 1.45 1.6 1.75 1.9 2.05 2.2 2.35 2.5 2.65 2.8 2.95 3.1 3.25 3.4 3.55 3.7 3.85 4.0 4.15 4.3 4.45 4.6 4.75 4.9 5.05 5.2 5.35 5.5 5.65 5.8 5.95 6.1 6.25 6.4 6.55 6.7 6.85 7.0 7.15 7.3 7.45 7.6 7.75 7.9 8.05 8.2 8.35 8.5 8.65 8.8 8.95 9.1 9.25 9.4 9.55 9.7 9.85 10.0 10.15
+1.0 1.15 1.3 1.45 1.6 1.75 1.9 2.05 2.2 2.35 2.5 2.65 2.8 2.95 3.1 3.25 3.4 3.55 3.7 3.85 4.0 4.15 4.3 4.45 4.6 4.75 4.9 5.05 5.2 5.35 5.5 5.65 5.8 5.95 6.1 6.25 6.4 6.55 6.7 6.85 7.0 7.15 7.3 7.45 7.6 7.75 7.9 8.05 8.2 8.35 8.5 8.65 8.8 8.95 9.1 9.25 9.4 9.55 9.7 9.85 10.0 10.15
diff --git a/test/test_vxl/pc0_grid_s1-1-2.txt b/test/test_vxl/pc0_grid_s1-1-2.txt
new file mode 100644
index 0000000..d4260bc
--- /dev/null
+++ b/test/test_vxl/pc0_grid_s1-1-2.txt
@@ -0,0 +1,3 @@
+1.0 2.0 3.0 4.0 5.0 6.0 7.0 8.0 9.0 10.0 11.0
+1.0 2.0 3.0 4.0 5.0 6.0 7.0 8.0 9.0 10.0 11.0
+1.0 3.0 5.0 7.0 9.0 11.0
diff --git a/test/test_vxl/pc0_grid_s1.txt b/test/test_vxl/pc0_grid_s1.txt
new file mode 100644
index 0000000..ae019d4
--- /dev/null
+++ b/test/test_vxl/pc0_grid_s1.txt
@@ -0,0 +1,3 @@
+1.0 2.0 3.0 4.0 5.0 6.0 7.0 8.0 9.0 10.0 11.0
+1.0 2.0 3.0 4.0 5.0 6.0 7.0 8.0 9.0 10.0 11.0
+1.0 2.0 3.0 4.0 5.0 6.0 7.0 8.0 9.0 10.0 11.0
diff --git a/test/test_vxl/pc0_grid_s2.txt b/test/test_vxl/pc0_grid_s2.txt
new file mode 100644
index 0000000..fd10de6
--- /dev/null
+++ b/test/test_vxl/pc0_grid_s2.txt
@@ -0,0 +1,3 @@
+1.0 3.0 5.0 7.0 9.0 11.0
+1.0 3.0 5.0 7.0 9.0 11.0
+1.0 3.0 5.0 7.0 9.0 11.0
diff --git a/test/test_vxl/pc0_vxl_s1.txt b/test/test_vxl/pc0_vxl_s1.txt
new file mode 100644
index 0000000..b5b4dc0
--- /dev/null
+++ b/test/test_vxl/pc0_vxl_s1.txt
@@ -0,0 +1,5 @@
+# x y z density mean mode 
+0 0 0 1 2 2
+0 2 1 4 10 5
+9 9 9 1 1 1
+4 4 4 1 0 0