star-vx

Structuring voxels for ray-tracing
git clone git://git.meso-star.fr/star-vx.git
Log | Files | Refs | README | LICENSE
commit a895f034092faa6eea8b4ea17d5730c21c927336
parent 3a5dd4bc4993914b918f777b8a69bcf4d8fa5e85
Author: Vincent Forest <vincent.forest@meso-star.com>
Date:   Fri,  2 Mar 2018 15:35:00 +0100

Rename the library in Star-VoXel

Diffstat:

MREADME.md | 4++--


Mcmake/CMakeLists.txt | 60++++++++++++++++++++++++++++++------------------------------


Dsrc/htvox.h | 196-------------------------------------------------------------------------------


Dsrc/htvox_c.h | 75---------------------------------------------------------------------------


Dsrc/htvox_device.c | 124-------------------------------------------------------------------------------


Dsrc/htvox_device.h | 39---------------------------------------


Dsrc/htvox_octree.c | 820-------------------------------------------------------------------------------


Dsrc/htvox_octree.h | 40----------------------------------------


Dsrc/htvox_octree_buffer.c | 198-------------------------------------------------------------------------------


Dsrc/htvox_octree_buffer.h | 239-------------------------------------------------------------------------------


Dsrc/htvox_scene_trace_ray.c | 228-------------------------------------------------------------------------------


Asrc/svx.h | 196+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++


Asrc/svx_c.h | 75+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++


Asrc/svx_device.c | 124+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++


Asrc/svx_device.h | 39+++++++++++++++++++++++++++++++++++++++


Asrc/svx_octree.c | 820+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++


Asrc/svx_octree.h | 40++++++++++++++++++++++++++++++++++++++++


Asrc/svx_octree_buffer.c | 198+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++


Asrc/svx_octree_buffer.h | 239+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++


Asrc/svx_scene_trace_ray.c | 228+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++


Dsrc/test_htvox_device.c | 65-----------------------------------------------------------------


Dsrc/test_htvox_octree.c | 434-------------------------------------------------------------------------------


Dsrc/test_htvox_utils.h | 34----------------------------------


Asrc/test_svx_device.c | 65+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++


Asrc/test_svx_octree.c | 434+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++


Asrc/test_svx_utils.h | 34++++++++++++++++++++++++++++++++++


26 files changed, 2524 insertions(+), 2524 deletions(-)

diff --git a/README.md b/README.md
@@ -1,4 +1,4 @@
-# High-Tune VoXel
+# Star-VoXel
 
 This library manage a 3D scene composed of voxels.
 
@@ -20,7 +20,7 @@ informations on CMake.
 
 ## Licenses
 
-This is a free software copyright (C) CNRS 2018 released under the GPL v3+
+This is a free software copyright (C) 2018 CNRS released under the GPL v3+
 license: GNU GPL version 3 or later. You are welcome to redistribute it under
 certain conditions; refer to the COPYING file for details.
 
diff --git a/cmake/CMakeLists.txt b/cmake/CMakeLists.txt
@@ -1,4 +1,4 @@
-# Copyright (C) CNRS 2018
+# Copyright (C) 2018 CNRS
 #
 # This program is free software: you can redistribute it and/or modify
 # it under the terms of the GNU General Public License as published by
@@ -14,10 +14,10 @@
 # along with this program. If not, see <http://www.gnu.org/licenses/>.
 
 cmake_minimum_required(VERSION 2.8)
-project(htvox C)
+project(svx C)
 enable_testing()
 
-set(HTVOX_SOURCE_DIR ${PROJECT_SOURCE_DIR}/../src)
+set(SVX_SOURCE_DIR ${PROJECT_SOURCE_DIR}/../src)
 option(NO_TEST "Do not build tests" OFF)
 
 ################################################################################
@@ -40,34 +40,34 @@ set(VERSION_MINOR 0)
 set(VERSION_PATCH 0)
 set(VERSION ${VERSION_MAJOR}.${VERSION_MINOR}.${VERSION_PATCH})
 
-set(HTVOX_FILES_SRC
-  htvox_device.c
-  htvox_octree.c
-  htvox_octree_buffer.c)
+set(SVX_FILES_SRC
+  svx_device.c
+  svx_octree.c
+  svx_octree_buffer.c)
 set(SSOL_FILES_INC
-  htvox_device.h
-  htvox_octree.h
-  htvox_octree_buffer.h)
+  svx_device.h
+  svx_octree.h
+  svx_octree_buffer.h)
 set(SSOL_FILES_INC_API
-  htvox.h)
+  svx.h)
 
-set(HTVOX_FILES_DOC COPYING README.md)
+set(SVX_FILES_DOC COPYING README.md)
 
-# Prepend each file in the `HTVOX_FILES_<SRC|INC>' list by `HTVOX_SOURCE_DIR'
-rcmake_prepend_path(HTVOX_FILES_SRC ${HTVOX_SOURCE_DIR})
-rcmake_prepend_path(HTVOX_FILES_INC ${HTVOX_SOURCE_DIR})
-rcmake_prepend_path(HTVOX_FILES_INC_API ${HTVOX_SOURCE_DIR})
-rcmake_prepend_path(HTVOX_FILES_DOC ${PROJECT_SOURCE_DIR}/../)
+# Prepend each file in the `SVX_FILES_<SRC|INC>' list by `SVX_SOURCE_DIR'
+rcmake_prepend_path(SVX_FILES_SRC ${SVX_SOURCE_DIR})
+rcmake_prepend_path(SVX_FILES_INC ${SVX_SOURCE_DIR})
+rcmake_prepend_path(SVX_FILES_INC_API ${SVX_SOURCE_DIR})
+rcmake_prepend_path(SVX_FILES_DOC ${PROJECT_SOURCE_DIR}/../)
 
-add_library(htvox SHARED ${HTVOX_FILES_SRC} ${HTVOX_FILES_INC} ${HTVOX_FILES_INC_API})
-target_link_libraries(htvox RSys)
+add_library(svx SHARED ${SVX_FILES_SRC} ${SVX_FILES_INC} ${SVX_FILES_INC_API})
+target_link_libraries(svx RSys)
 
-set_target_properties(htvox PROPERTIES
-  DEFINE_SYMBOL HTVOX_SHARED_BUILD
+set_target_properties(svx PROPERTIES
+  DEFINE_SYMBOL SVX_SHARED_BUILD
   VERSION ${VERSION}
   SOVERSION ${VERSION_MAJOR})
 
-rcmake_setup_devel(htvox HTVOX ${VERSION} high_tune/htvox.h)
+rcmake_setup_devel(svx SVX ${VERSION} high_tune/svx.h)
 
 ################################################################################
 # Add tests
@@ -75,9 +75,9 @@ rcmake_setup_devel(htvox HTVOX ${VERSION} high_tune/htvox.h)
 if(NOT NO_TEST)
   function(build_test _name)
     add_executable(${_name}
-      ${HTVOX_SOURCE_DIR}/${_name}.c
-      ${HTVOX_SOURCE_DIR}/test_htvox_utils.h)
-    target_link_libraries(${_name} htvox RSys)
+      ${SVX_SOURCE_DIR}/${_name}.c
+      ${SVX_SOURCE_DIR}/test_svx_utils.h)
+    target_link_libraries(${_name} svx RSys)
   endfunction()
 
   function(new_test _name)
@@ -85,17 +85,17 @@ if(NOT NO_TEST)
     add_test(${_name} ${_name})
   endfunction()
 
-  new_test(test_htvox_device)
-  new_test(test_htvox_octree)
+  new_test(test_svx_device)
+  new_test(test_svx_octree)
 endif()
 
 ################################################################################
 # Define output & install directories
 ################################################################################
-install(TARGETS htvox
+install(TARGETS svx
   ARCHIVE DESTINATION bin
   LIBRARY DESTINATION lib
   RUNTIME DESTINATION bin)
-install(FILES ${HTVOX_FILES_INC_API} DESTINATION include/high_tune)
-install(FILES ${HTVOX_FILES_DOC} DESTINATION share/doc/htvox)
+install(FILES ${SVX_FILES_INC_API} DESTINATION include/high_tune)
+install(FILES ${SVX_FILES_DOC} DESTINATION share/doc/svx)
 
diff --git a/src/htvox.h b/src/htvox.h
@@ -1,196 +0,0 @@
-/* Copyright (C) CNRS 2018
- *
- * This program is free software: you can redistribute it and/or modify
- * it under the terms of the GNU General Public License as published by
- * the Free Software Foundation, either version 3 of the License, or
- * (at your option) any later version.
- *
- * This program is distributed in the hope that it will be useful,
- * but WITHOUT ANY WARRANTY; without even the implied warranty of
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
- * GNU General Public License for more details.
- *
- * You should have received a copy of the GNU General Public License
- * along with this program. If not, see <http://www.gnu.org/licenses/>. */
-
-#ifndef HTVOX_H
-#define HTVOX_H
-
-#include <rsys/rsys.h>
-#include <float.h>
-
-/* Library symbol management */
-#if defined(HTVOX_SHARED_BUILD) /* Build shared library */
-  #define HTVOX_API extern EXPORT_SYM
-#elif defined(HTVOX_STATIC) /* Use/build static library */
-  #define HTVOX_API extern LOCAL_SYM
-#else /* Use shared library */
-  #define HTVOX_API extern IMPORT_SYM
-#endif
-
-/* Helper macro that asserts if the invocation of the htvox function `Func'
- * returns an error. One should use this macro on htvox function calls for
- * which no explicit error checking is performed */
-#ifndef NDEBUG
-  #define HTVOX(Func) ASSERT(htvox_ ## Func == RES_OK)
-#else
-  #define HTVOX(Func) htvox_ ## Func
-#endif
-
-/* Maximum memory size of a voxel */
-#define HTVOX_MAX_SIZEOF_VOXEL (sizeof(double)*16)
-
-struct htvox_voxel {
-  const void* data; /* Data of the voxel */
-  size_t id; /* Indentifier of the voxel */
-  size_t depth; /* Depth of the voxel, in [0, htvox_octree_desc.depth[ */
-  int is_leaf; /* Define if the voxel is a leaf into the hierarchy */
-};
-
-#define HTVOX_VOXEL_NULL__ { NULL, SIZE_MAX, SIZE_MAX, 0 }
-static const struct htvox_voxel HTVOX_VOXEL_NULL = HTVOX_VOXEL_NULL__;
-
-#define HTVOX_VOXEL_NONE(Voxel) ((Voxel)->id == HTVOX_VOXEL_NULL.id)
-
-/* Descriptor of a voxel */
-struct htvox_voxel_desc {
-  /* Retrieve the data of the voxels*/
-  void
-  (*get)
-    (const size_t xyz[3], /* Voxel coordinate in voxel space */
-     void* dst, /* Where to store data */
-     void* ctx); /* Pointer toward user data */
-
-  /* Merge the data of N voxels */
-  void
-  (*merge)
-    (void* dst, /* Merged data */
-     const void* voxels[], /* Data to merge */
-     const size_t nvoxels, /* #submitted data */
-     void* ctx); /* Pointer toward user data */
-
-  /* Check if the voxel's data can be merged */
-  int
-  (*challenge_merge)
-    (const void* voxels[], /* Data candidates to the merge */
-     const size_t nvoxels, /* #candidates */
-     void* ctx); /* Pointer toward user data */
-
-  void* context; /* Client side data sent as the last argument of the clbbs */
-  size_t size; /* Size in bytes of a voxel. Must be <= HTVOX_MAX_SIZEOF_VOXE */
-};
-
-#define HTVOX_VOXEL_DESC_NULL__ { NULL, NULL, NULL, NULL, 0 }
-static const struct htvox_voxel_desc HTVOX_VOXEL_DESC_NULL =
-  HTVOX_VOXEL_DESC_NULL__;
-
-struct htvox_octree_desc {
-  /* Submitted Axis Aligned Bounding Box */
-  double lower[3], upper[3];
-
-  size_t nleaves; /* #leaves */
-  size_t nvoxels; /* #voxels, i.e. #leaves + #parents */
-  size_t depth; /* Depth of the octree */
-};
-
-#define HTVOX_OCTREE_DESC_NULL__ \
-  {{DBL_MAX, DBL_MAX, DBL_MAX}, {-DBL_MAX,-DBL_MAX,-DBL_MAX}, 0, 0, 0}
-static const struct htvox_octree_desc HTVOX_OCTREE_DESC_NULL =
-  HTVOX_OCTREE_DESC_NULL__;
-
-struct htvox_hit {
-  double distance; /* Distance from the ray origin to the impacted voxel */
-  struct htvox_voxel voxel; /* Intersected voxel */
-};
-
-#define HTVOX_HIT_NULL__ {DBL_MAX, HTVOX_VOXEL_NULL__}
-static const struct htvox_hit HTVOX_HIT_NULL = HTVOX_HIT_NULL__;
-
-#define HTVOX_HIT_NONE(Hit) ((Hit)->distance >= DBL_MAX)
-
-/* Forward declaration of external data types */
-struct logger;
-struct mem_allocator;
-
-/* Forward declaration of opaque data types */
-struct htvox_device;
-struct htvox_octree;
-
-BEGIN_DECLS
-
-/*******************************************************************************
- * Device
- ******************************************************************************/
-HTVOX_API res_T
-htvox_device_create
-  (struct logger* logger,
-   struct mem_allocator* allocator, /* NULL <=> use default allocator */
-   const int verbose, /* Verbosity level */
-   struct htvox_device** htvox);
-
-HTVOX_API res_T
-htvox_device_ref_get
-  (struct htvox_device* htvox);
-
-HTVOX_API res_T
-htvox_device_ref_put
-  (struct htvox_device* htvox);
-
-/*******************************************************************************
- * Octree
- ******************************************************************************/
-HTVOX_API res_T
-htvox_octree_create
-  (struct htvox_device* dev,
-   const double lower[3], /* Lower bound of the octree */
-   const double upper[3], /* Upper bound of the octree */
-   const size_t nvoxels[3], /* # voxels along the 3 axis */
-   const struct htvox_voxel_desc* desc, /* Descriptor of a voxel */
-   struct htvox_octree** octree);
-
-HTVOX_API res_T
-htvox_octree_ref_get
-  (struct htvox_octree* octree);
-
-HTVOX_API res_T
-htvox_octree_ref_put
-  (struct htvox_octree* octree);
-
-HTVOX_API res_T
-htvox_octree_get_desc
-  (const struct htvox_octree* octree,
-   struct htvox_octree_desc* desc);
-
-HTVOX_API res_T
-htvox_octree_for_each_leaf
-  (struct htvox_octree* octree,
-   void (*functor)
-    (const void* val, /* Value of the voxel */
-     const size_t ileaf, /* Identifier of the leaf in [0, #leafs[ */
-     const double low[3], /* World space lower bound of the voxel */
-     const double upp[3], /* World space upper bound of the voxel */
-     void* ctx), /* Client data */
-   void* context); /* Client data sent as the last argument of the callback */
-
-HTVOX_API res_T
-htvox_octree_trace_ray
-  (struct htvox_octree* octree,
-   const double ray_origin[3],
-   const double ray_direction[3],
-   const double ray_range[2],
-   struct htvox_hit* hit);
-
-HTVOX_API res_T
-htvox_octree_at
-  (struct htvox_octree* octree,
-   const double position[3],
-   int (*filter)  /* Filter function. May be NULL <=> traverse up to leaves */
-    (const struct htvox_voxel* voxel,
-     const double position[3],
-     void* ctx),
-   void* context, /* Client data sent as the last argument of the filter func */
-   struct htvox_voxel* voxel);
-
-#endif /* HTVOX_H */
-
-
diff --git a/src/htvox_c.h b/src/htvox_c.h
@@ -1,75 +0,0 @@
-/* Copyright (C) CNRS 2018
- *
- * This program is free software: you can redistribute it and/or modify
- * it under the terms of the GNU General Public License as published by
- * the Free Software Foundation, either version 3 of the License, or
- * (at your option) any later version.
- *
- * This program is distributed in the hope that it will be useful,
- * but WITHOUT ANY WARRANTY; without even the implied warranty of
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
- * GNU General Public License for more details.
- *
- * You should have received a copy of the GNU General Public License
- * along with this program. If not, see <http://www.gnu.org/licenses/>. */
-
-#ifndef HTVOX_C_H
-#define HTVOX_C_H
-
-#include <rsys/rsys.h>
-
-static INLINE uint64_t
-morton3D_encode_u21(const uint32_t u21)
-{
-  uint64_t u64 = u21 & ((1<<21) - 1);
-  ASSERT(u21 <= ((1 << 21) - 1));
-  u64 = (u64 | (u64 << 32)) & 0xFFFF00000000FFFF;
-  u64 = (u64 | (u64 << 16)) & 0x00FF0000FF0000FF;
-  u64 = (u64 | (u64 << 8))  & 0xF00F00F00F00F00F;
-  u64 = (u64 | (u64 << 4))  & 0x30C30C30C30C30C3;
-  u64 = (u64 | (u64 << 2))  & 0x9249249249249249;
-  return u64;
-}
-
-static INLINE uint32_t
-morton3D_decode_u21(const uint64_t u64)
-{
-  uint64_t tmp = (u64 & 0x9249249249249249);
-  tmp = (tmp | (tmp >> 2))  & 0x30C30C30C30C30C3;
-  tmp = (tmp | (tmp >> 4))  & 0xF00F00F00F00F00F;
-  tmp = (tmp | (tmp >> 8))  & 0x00FF0000FF0000FF;
-  tmp = (tmp | (tmp >> 16)) & 0xFFFF00000000FFFF;
-  tmp = (tmp | (tmp >> 32)) & 0x00000000FFFFFFFF;
-  ASSERT(tmp <= ((1<<21)-1));
-  return (uint32_t)tmp;
-}
-
-/* Count the number of bits set to 1 */
-static FINLINE int
-popcount(const uint8_t x)
-{
-  int n = x - ((x >> 1) & 0x55);
-  n = (n & 0x33) + ((n >> 2) & 0x33);
-  n = (n + (n >> 4)) & 0x0f;
-  return (n * 0x0101) >> 8;
-}
-
-static INLINE uint64_t
-morton_xyz_encode_u21(const uint32_t xyz[3])
-{
-  return (morton3D_encode_u21(xyz[0]) << 2)
-       | (morton3D_encode_u21(xyz[1]) << 1)
-       | (morton3D_encode_u21(xyz[2]) << 0);
-}
-
-static INLINE void
-morton_xyz_decode_u21(const uint64_t code, uint32_t xyz[3])
-{
-  ASSERT(xyz && code < ((1ul << 63)-1));
-  xyz[0] = (uint32_t)morton3D_decode_u21(code >> 2);
-  xyz[1] = (uint32_t)morton3D_decode_u21(code >> 1);
-  xyz[2] = (uint32_t)morton3D_decode_u21(code >> 0);
-}
-
-#endif /* HTVOX_C_H */
-
diff --git a/src/htvox_device.c b/src/htvox_device.c
@@ -1,124 +0,0 @@
-/* Copyright (C) CNRS 2018
- *
- * This program is free software: you can redistribute it and/or modify
- * it under the terms of the GNU General Public License as published by
- * the Free Software Foundation, either version 3 of the License, or
- * (at your option) any later version.
- *
- * This program is distributed in the hope that it will be useful,
- * but WITHOUT ANY WARRANTY; without even the implied warranty of
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
- * GNU General Public License for more details.
- *
- * You should have received a copy of the GNU General Public License
- * along with this program. If not, see <http://www.gnu.org/licenses/>. */
-
-#include "htvox.h"
-#include "htvox_device.h"
-
-#include <rsys/logger.h>
-#include <rsys/mem_allocator.h>
-
-/*******************************************************************************
- * Helper functions
- ******************************************************************************/
-static void
-log_msg
-  (const struct htvox_device* dev,
-   const enum log_type stream,
-   const char* msg,
-   va_list vargs)
-{
-  ASSERT(dev && msg);
-  if(dev->verbose) {
-    res_T res; (void)res;
-    res = logger_vprint(dev->logger, stream, msg, vargs);
-    ASSERT(res == RES_OK);
-  }
-}
-
-static void
-device_release(ref_T* ref)
-{
-  struct htvox_device* dev;
-  ASSERT(ref);
-  dev = CONTAINER_OF(ref, struct htvox_device, ref);
-  MEM_RM(dev->allocator, dev);
-}
-
-/*******************************************************************************
- * Exported functions
- ******************************************************************************/
-res_T
-htvox_device_create
-  (struct logger* log,
-   struct mem_allocator* mem_allocator,
-   const int verbose,
-   struct htvox_device** out_dev)
-{
-  struct htvox_device* dev = NULL;
-  struct mem_allocator* allocator = NULL;
-  struct logger* logger = NULL;
-  res_T res = RES_OK;
-
-  if(!out_dev) {
-    res = RES_BAD_ARG;
-    goto error;
-  }
-
-  allocator = mem_allocator ? mem_allocator : &mem_default_allocator;
-  logger = log ? log : LOGGER_DEFAULT;
-
-  dev = MEM_CALLOC(allocator, 1, sizeof(struct htvox_device));
-  if(!dev) {
-    res = RES_MEM_ERR;
-    goto error;
-  }
-
-
-  ref_init(&dev->ref);
-  dev->allocator = allocator;
-  dev->logger = logger;
-  dev->verbose = verbose;
-
-exit:
-  if(out_dev) *out_dev = dev;
-  return res;
-error:
-  if(dev) {
-    HTVOX(device_ref_put(dev));
-    dev = NULL;
-  }
-  goto exit;
-}
-
-res_T
-htvox_device_ref_get(struct htvox_device* dev)
-{
-  if(!dev) return RES_BAD_ARG;
-  ref_get(&dev->ref);
-  return RES_OK;
-}
-
-res_T
-htvox_device_ref_put(struct htvox_device* dev)
-{
-  if(!dev) return RES_BAD_ARG;
-  ref_put(&dev->ref, device_release);
-  return RES_OK;
-}
-
-/*******************************************************************************
- * Local function
- ******************************************************************************/
-void
-log_err(const struct htvox_device* dev, const char* msg, ...)
-{
-  va_list vargs_list;
-  ASSERT(dev && msg);
-
-  va_start(vargs_list, msg);
-  log_msg(dev, LOG_ERROR, msg, vargs_list);
-  va_end(vargs_list);
-}
-
diff --git a/src/htvox_device.h b/src/htvox_device.h
@@ -1,39 +0,0 @@
-/* Copyright (C) CNRS 2018
- *
- * This program is free software: you can redistribute it and/or modify
- * it under the terms of the GNU General Public License as published by
- * the Free Software Foundation, either version 3 of the License, or
- * (at your option) any later version.
- *
- * This program is distributed in the hope that it will be useful,
- * but WITHOUT ANY WARRANTY; without even the implied warranty of
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
- * GNU General Public License for more details.
- *
- * You should have received a copy of the GNU General Public License
- * along with this program. If not, see <http://www.gnu.org/licenses/>. */
-
-#ifndef HTVOX_DEVICE_H
-#define HTVOX_DEVICE_H
-
-#include <rsys/ref_count.h>
-
-struct htvox_device {
-  int verbose; /* Verbosity level */
-  struct logger* logger;
-  struct mem_allocator* allocator;
-  ref_T ref;
-};
-
-extern LOCAL_SYM void
-log_err
-  (const struct htvox_device* dev,
-   const char* fmt,
-  ...)
-#ifdef COMPILER_GCC
-  __attribute((format(printf, 2, 3)))
-#endif
-  ;
-
-#endif /* HTVOX_DEVICE_H */
-
diff --git a/src/htvox_octree.c b/src/htvox_octree.c
@@ -1,820 +0,0 @@
-/* Copyright (C) CNRS 2018
- *
- * This program is free software: you can redistribute it and/or modify
- * it under the terms of the GNU General Public License as published by
- * the Free Software Foundation, either version 3 of the License, or
- * (at your option) any later version.
- *
- * This program is distributed in the hope that it will be useful,
- * but WITHOUT ANY WARRANTY; without even the implied warranty of
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
- * GNU General Public License for more details.
- *
- * You should have received a copy of the GNU General Public License
- * along with this program. If not, see <http://www.gnu.org/licenses/>. */
-
-#include "htvox.h"
-#include "htvox_c.h"
-#include "htvox_device.h"
-#include "htvox_octree.h"
-
-#include <rsys/double3.h>
-#include <rsys/mem_allocator.h>
-
-#define OCTREE_DEPTH_MAX 16 /* Maximum depth of an octree */
-
-struct voxel {
-  uint64_t mcode; /* Morton code of the voxel */
-  void* data; /* Data of the voxel */
-};
-static const struct voxel VOXEL_NULL = {0, NULL};
-
-struct octree_node {
-  struct octree_index ichild_node; /* Index of the 1st child node */
-  struct octree_index ichild_attr; /* Index of the 1st child attr */
-  uint8_t is_valid; /* Mask defining whether the children are valid or not */
-  uint8_t is_leaf; /* Mask defining whether the children are leaves or not */
-  ALIGN(16) char data[8][HTVOX_MAX_SIZEOF_VOXEL]; /* Data of the leaves */
-};
-
-/* Stacked children of an octree node */
-struct stack {
-  struct octree_node nodes[8]; /* List of registered children nodes */
-  uint8_t mask; /* Mask of valid children nodes (0 = empty) */
-};
-
-struct octree_builder {
-  struct stack stacks[OCTREE_DEPTH_MAX];
-  struct octree_buffer* buffer;
-
-  const struct htvox_voxel_desc* desc;
-  size_t nleaves; /* Number of emitted leaves */
-
-  int octree_depth; /* Maximum octree depth */
-  int non_empty_lvl; /* Index of the 1st non empty level */
-
-  uint64_t mcode; /* Morton code of the last registered voxels */
-};
-
-/*******************************************************************************
- * Helper function
- ******************************************************************************/
-#ifndef NDEBUG
-static void
-check_octree(struct octree_buffer* buf, const struct octree_index root)
-{
-  const struct octree_xnode* node;
-  int ichild;
-  ASSERT(buf);
-
-  node = octree_buffer_get_node(buf, root);
-  FOR_EACH(ichild, 0, 8) {
-    const int ichild_flag = BIT(ichild);
-    if((node->is_valid & ichild_flag) == 0) continue;
-
-    if(node->is_leaf & ichild_flag) {
-      struct octree_index iattr;
-      iattr = octree_buffer_get_child_attr_index(buf, root, ichild);
-      ASSERT(octree_buffer_get_attr(buf, iattr) != NULL);
-    } else {
-      struct octree_index child;
-      child = octree_buffer_get_child_node_index(buf, root, ichild);
-      check_octree(buf, child);
-    }
-  }
-}
-#endif
-
-static FINLINE size_t
-absolute_attr_index
-  (const struct octree_buffer* buf,
-   const struct octree_index index)
-{
-  ASSERT(buf);
-  return index.ipage * buf->pagesize/buf->voxsize + index.inode;
-
-}
-
-static INLINE int
-check_htvox_voxel_desc(const struct htvox_voxel_desc* desc)
-{
-  return desc
-      && desc->get
-      && desc->merge
-      && desc->challenge_merge
-      && desc->size > 0
-      && desc->size <= HTVOX_MAX_SIZEOF_VOXEL;
-}
-
-static INLINE void
-stack_clear(struct stack* stack)
-{
-  int inode;
-  FOR_EACH(inode, 0, 8) {
-    int ichild;
-    stack->nodes[inode].is_leaf = 0;
-    stack->nodes[inode].is_valid = 0;
-    stack->nodes[inode].ichild_node = OCTREE_INDEX_NULL;
-    stack->nodes[inode].ichild_attr = OCTREE_INDEX_NULL;
-    FOR_EACH(ichild, 0, 8) {
-      memset(stack->nodes[inode].data[ichild], 0, HTVOX_MAX_SIZEOF_VOXEL);
-    }
-  }
-  stack->mask = 0;
-}
-
-/* Build a parent octree node from the registered stack nodes */
-static INLINE void
-stack_setup_node
-  (struct stack* stack,
-   struct octree_node* node,
-   const struct htvox_voxel_desc* desc)
-{
-  int ichild;
-  ASSERT(stack && node && check_htvox_voxel_desc(desc));
-
-  node->ichild_node = OCTREE_INDEX_NULL;
-  node->ichild_attr = OCTREE_INDEX_NULL;
-  node->is_valid = stack->mask;
-  node->is_leaf = 0;
-
-  if(stack->mask == 0) return; /* Empty stack */
-
-  /* Try to merge the child's leaves */
-  FOR_EACH(ichild, 0, 8) {
-    const void* data[8];
-    const uint8_t ichild_flag = (uint8_t)BIT(ichild);
-    struct octree_node* child = stack->nodes + ichild;
-    int igrandchild;
-
-    /* Fetch the grandchildren data */
-    FOR_EACH(igrandchild, 0, 8) {
-      data[igrandchild] = child->data[igrandchild];
-    }
-
-    desc->merge(node->data[ichild], data, 8, desc->context);
-
-    if(child->is_leaf==0xFF && desc->challenge_merge(data, 8, desc->context)) {
-      /* The node becomes a leaf : the children does not exist anymore */
-      node->is_leaf |= ichild_flag;
-      stack->mask ^= ichild_flag;
-    }
-  }
-}
-
-static res_T
-stack_write
-  (struct stack* stack, /* Node to write */
-   struct octree_buffer* buf, /* Buffer where nodes are written */
-   struct octree_index* out_index, /* Index of the first written node */
-   size_t* out_nleaves) /* #writen leaves */
-{
-  struct octree_index nodes_id = OCTREE_INDEX_NULL;
-  struct octree_node* node = NULL;
-  size_t nleaves = 0;
-  int inode;
-  res_T res = RES_OK;
-  ASSERT(stack && buf && out_index && out_nleaves);
-
-  /* No registered nodes, this means that the nodes were merged in an higher
-   * level */
-  if(!stack->mask) goto exit;
-
-  /* Write the attrib of the children */
-  FOR_EACH(inode, 0, 8) {
-    char* data = NULL;
-    size_t nattrs = 0;
-    size_t nvoxs = 0;
-    int ichild = 0;
-
-    if((stack->mask & BIT(inode)) == 0) continue; /* Empty node */
-
-    node = stack->nodes + inode;
-
-    nattrs =  (size_t)popcount(node->is_valid);
-    ASSERT(nattrs > 0);
-
-    res = octree_buffer_alloc_attrs(buf, nattrs, &node->ichild_attr);
-    if(res != RES_OK) goto error;
-
-    data = octree_buffer_get_attr(buf, node->ichild_attr);
-    nvoxs = 0;
-    FOR_EACH(ichild, 0, 8) {
-      if(!(node->is_valid & BIT(ichild))) continue;
-      memcpy(data + nvoxs*buf->voxsize, node->data[ichild], buf->voxsize);
-      ++nvoxs;
-    }
-    ASSERT(nvoxs == nattrs);
-    nleaves += (size_t)popcount(node->is_leaf);
-  }
-
-  do {
-    struct octree_index index = OCTREE_INDEX_NULL;
-    struct octree_xnode* xnodes = NULL;
-    const size_t nnodes = (size_t)popcount(stack->mask);
-    size_t ixnode = 0;
-
-    /* Alloc the octree nodes */
-    res = octree_buffer_alloc_nodes(buf, (size_t)nnodes, &nodes_id);
-    if(res != RES_OK) goto error;
-    xnodes = octree_buffer_get_node(buf, nodes_id);
-
-    FOR_EACH(inode, 0, 8) {
-      uint16_t attr_offset = UINT16_MAX;
-      uint16_t node_offset = UINT16_MAX;
-
-      if((stack->mask & BIT(inode)) == 0) continue; /* Empty node */
-      node = stack->nodes + inode;
-
-      /* Setup the offset toward the children and children attribs */
-      if(node->ichild_node.ipage == nodes_id.ipage) {
-        node_offset = node->ichild_node.inode;
-        attr_offset = node->ichild_attr.inode;
-      }
-
-      /* The page id of the children is not the same as that of node */
-      if(node_offset > OCTREE_XNODE_MAX_CHILDREN_OFFSET) {
-        res = octree_buffer_alloc_far_index(buf, &index);
-        if(res != RES_OK) break;
-        *octree_buffer_get_far_index(buf, index) = node->ichild_node;
-        node_offset = OCTREE_XNODE_FLAG_FAR_INDEX | index.inode;
-      }
-
-      /* The page id of the attribs is not tthe same as that of node */
-      if(attr_offset > OCTREE_XNODE_FLAG_FAR_INDEX) {
-        res = octree_buffer_alloc_far_index(buf, &index);
-        if(res != RES_OK) break;
-        *octree_buffer_get_far_index(buf, index) = node->ichild_attr;
-        attr_offset = OCTREE_XNODE_FLAG_FAR_INDEX | index.inode;
-      }
-
-      xnodes[ixnode].node_offset = node_offset;
-      xnodes[ixnode].attr_offset = attr_offset;
-      xnodes[ixnode].is_valid = node->is_valid;
-      xnodes[ixnode].is_leaf = node->is_leaf;
-      ++ixnode;
-    }
-  /* inode < 8 <=> not enough memory in the current page. A far index could not
-   * be stored in the same page of its associated node. The write process was
-   * stoped. Rewrite the whole stacked nodes in a new page. */
-  } while(inode < 8);
-
-
-exit:
-  /* Return the index toward the first writen nodes */
-  *out_index = nodes_id;
-  *out_nleaves = nleaves;
-  return res;
-error:
-  goto exit;
-}
-
-static res_T
-octree_builder_init
-  (struct octree_builder* bldr,
-   const size_t definition,
-   const struct htvox_voxel_desc* desc,
-   struct octree_buffer* buffer)
-{
-  int ilvl;
-  res_T res = RES_OK;
-  ASSERT(bldr && IS_POW2(definition) && check_htvox_voxel_desc(desc));
-  memset(bldr, 0, sizeof(struct octree_builder));
-
-  /* Compute the maximum depth of the octree */
-  bldr->octree_depth = log2i((int)definition);
-  if(bldr->octree_depth > OCTREE_DEPTH_MAX) {
-    res = RES_MEM_ERR;
-    goto error;
-  }
-
-  /* Init the per octree level stack */
-  FOR_EACH(ilvl, 0, bldr->octree_depth) {
-    stack_clear(&bldr->stacks[ilvl]);
-  }
-
-  octree_buffer_clear(buffer);
-  bldr->nleaves = 0;
-  bldr->desc = desc;
-  bldr->buffer = buffer;
-  bldr->non_empty_lvl = bldr->octree_depth - 1;
-
-exit:
-  return res;
-error:
-  goto exit;
-}
-
-static res_T
-octree_builder_add_voxel
-  (struct octree_builder* bldr,
-   const struct voxel* vox)
-{
-  uint64_t mcode_xor;
-  int inode;
-  int ichild;
-  uint8_t ichild_flag;
-  res_T res = RES_OK;
-  ASSERT(bldr && vox && vox->mcode >= bldr->mcode);
-
-  /* Define if the bits in [4 .. 63] of the previous and the next Morton
-   * codes are the same */
-  mcode_xor = bldr->mcode ^ vox->mcode;
-
-  /* The next voxel is not in the current node */
-  if(mcode_xor >= 8) {
-    size_t ilvl;
-
-    inode = (bldr->mcode >> 3) & 7;
-    bldr->stacks[0].mask |= (uint8_t)BIT(inode);
-
-    /* Flush the stack of the octree level that does not contain the next
-     * voxel. The last octree level is actually the level that contains the
-     * root node while the penultimate describes the root node itself. These 2
-     * levels contain all the voxels and can be skipped */
-    FOR_EACH(ilvl, 0, bldr->octree_depth-2/*The 2 last leaves contain all voxels*/) {
-      struct octree_node* stack_node;
-      uint64_t mcode_max_lvl;
-      size_t nleaves;
-
-      /* Compute the maximum morton code value for the current octree level */
-      mcode_max_lvl = 8lu/*#children*/ * (1lu<<(3*(ilvl+1)))/*#voxels per children*/;
-
-      if(mcode_xor < mcode_max_lvl) break;
-
-      /* Retrieve the node index of the next level */
-      inode = (bldr->mcode >> (3*(ilvl+2))) & 7;
-
-      /* The next voxel is not in the ilvl^th stack. Setup the parent node of the
-       * nodes registered into the stack */
-      stack_node = &bldr->stacks[ilvl+1].nodes[inode];
-      stack_setup_node(&bldr->stacks[ilvl], stack_node, bldr->desc);
-      bldr->stacks[ilvl+1].mask |= (uint8_t)BIT(inode);
-
-      /* Write the nodes of the stack of the current octree level into the buf */
-      res = stack_write
-        (&bldr->stacks[ilvl], bldr->buffer, &stack_node->ichild_node, &nleaves);
-      if(res != RES_OK) goto error;
-
-      bldr->nleaves += nleaves;
-      if(nleaves) bldr->non_empty_lvl = MMIN(bldr->non_empty_lvl, (int)ilvl);
-
-      /* Reset the current stack */
-      stack_clear(&bldr->stacks[ilvl]);
-    }
-  }
-
-  /* Retrieve the index of the current voxel and of its parent */
-  ichild = vox->mcode & 7;
-  inode = (vox->mcode >> 3) & 7;
-  ichild_flag = (uint8_t)BIT(ichild);
-
-  /* Register the voxel */
-  memcpy(bldr->stacks[0].nodes[inode].data[ichild], vox->data, bldr->desc->size);
-  bldr->stacks[0].nodes[inode].is_valid |= ichild_flag;
-  bldr->stacks[0].nodes[inode].is_leaf |= ichild_flag;
-
-  /* Update morton code of the last registered voxel */
-  bldr->mcode = vox->mcode;
-
-exit:
-  return res;
-error:
-  goto exit;
-}
-
-static INLINE res_T
-octree_builder_finalize
-  (struct octree_builder* bldr,
-   struct octree_index* root_id,
-   void* root_data)
-{
-  const void* data[8];
-  size_t inode;
-  size_t nleaves;
-  int ilvl;
-  int ichild;
-  res_T res = RES_OK;
-  ASSERT(bldr);
-
-  inode = (bldr->mcode >> 3) & 7;
-  bldr->stacks[0].mask |= (uint8_t)BIT(inode);
-
-  /* Flush the stacked nodes */
-  FOR_EACH(ilvl, 0, bldr->octree_depth-1) {
-    struct octree_node* parent_node;
-
-    if(bldr->stacks[ilvl].mask == 0) continue;
-
-    /* Retrieve the node index of the next level */
-    inode = (bldr->mcode >> (3*(ilvl+2))) & 7;
-
-    /* Setup the parent node of the nodes registered into the current stack */
-    parent_node = &bldr->stacks[ilvl+1].nodes[inode]; /* Fetch the parent node */
-    stack_setup_node(&bldr->stacks[ilvl], parent_node, bldr->desc);
-    bldr->stacks[ilvl+1].mask |= (uint8_t)BIT(inode);
-
-    /* Write the stacked nodes of the current level */
-    res = stack_write
-      (&bldr->stacks[ilvl], bldr->buffer, &parent_node->ichild_node, &nleaves);
-    if(res != RES_OK) goto error;
-    bldr->nleaves += nleaves;
-  }
-
-  ilvl = bldr->octree_depth-1; /* Root level */
-
-  /* Write the root node */
-  res = stack_write(&bldr->stacks[ilvl], bldr->buffer, root_id, &nleaves);
-  if(res != RES_OK) goto error;
-  bldr->nleaves += nleaves;
-
-  /* Setup the root attribs */
-  ASSERT(bldr->stacks[ilvl].mask == 1); /* Only the root node is active */
-  FOR_EACH(ichild, 0, 8) {
-    data[ichild] = bldr->stacks[ilvl].nodes[0].data[ichild];
-  }
-  bldr->desc->merge(root_data, data, 8, bldr->desc->context);
-
-exit:
-  return res;
-error:
-  goto exit;
-}
-
-static void
-octree_release(ref_T* ref)
-{
-  struct htvox_octree* oct;
-  struct htvox_device* dev;
-  ASSERT(ref);
-  oct = CONTAINER_OF(ref, struct htvox_octree, ref);
-  octree_buffer_release(&oct->buffer);
-  dev = oct->dev;
-  MEM_RM(dev->allocator, oct);
-  HTVOX(device_ref_put(dev));
-}
-
-/*******************************************************************************
- * Exported functions
- ******************************************************************************/
-res_T
-htvox_octree_create
-  (struct htvox_device* dev,
-   const double lower[3], /* Lower bound of the octree */
-   const double upper[3], /* Upper bound of the octree */
-   const size_t nvoxels[3], /* # voxels along the 3 axis */
-   const struct htvox_voxel_desc* desc, /* Descriptor of a voxel */
-   struct htvox_octree** out_oct)
-{
-  struct htvox_octree* oct = NULL;
-  double vox_sz[3]; /* World space size of a voxel */
-  struct octree_builder bldr;
-  struct voxel vox = VOXEL_NULL;
-  uint64_t mcode_max;
-  uint64_t mcode;
-  res_T res = RES_OK;
-
-  if(!dev || !lower || !upper || !nvoxels
-  || !check_htvox_voxel_desc(desc) || !out_oct) {
-    res = RES_BAD_ARG;
-    goto error;
-  }
-  if(lower[0] >= upper[0]
-  || lower[1] >= upper[1]
-  || lower[2] >= upper[2]) {
-    log_err(dev,
-      "%s: the submitted AABB is degenerated.\n"
-      "\tlower={%g, %g, %g}, upper={%g, %g, %g}.\n",
-      FUNC_NAME, SPLIT3(lower), SPLIT3(upper));
-    res = RES_BAD_ARG;
-    goto error;
-  }
-  if(!nvoxels[0] || !nvoxels[1] || !nvoxels[2]) {
-    log_err(dev,
-      "%s: the number of voxels along each axis cannot be null.\n"
-      "\t#voxels XYZ = {%lu, %lu, %lu}.\n",
-      FUNC_NAME,
-      (unsigned long)nvoxels[0],
-      (unsigned long)nvoxels[1],
-      (unsigned long)nvoxels[2]);
-    res = RES_BAD_ARG;
-    goto error;
-  }
-  oct = MEM_CALLOC(dev->allocator, 1, sizeof(struct htvox_octree));
-  if(!oct) {
-    res = RES_MEM_ERR;
-    goto error;
-  }
-  ref_init(&oct->ref);
-  HTVOX(device_ref_get(dev));
-  octree_buffer_init(dev->allocator, desc->size, &oct->buffer);
-  oct->dev = dev;
-
-  /* Compute the octree definition */
-  oct->definition = MMAX(nvoxels[0], 2);
-  oct->definition = MMAX(nvoxels[1], oct->definition);
-  oct->definition = MMAX(nvoxels[2], oct->definition);
-  oct->definition = round_up_pow2(oct->definition);
-
-  /* Intialize the octree builder */
-  res = octree_builder_init(&bldr, oct->definition, desc, &oct->buffer);
-  if(res != RES_OK) goto error;
-
-  vox.data = MEM_CALLOC(dev->allocator, 1, desc->size);
-  if(!vox.data) {
-    res = RES_MEM_ERR;
-    goto error;
-  }
-
-  mcode_max = oct->definition * oct->definition * oct->definition;
-  FOR_EACH(mcode, 0, mcode_max) {
-    size_t xyz[3];
-    uint32_t ui3[3];
-
-    morton_xyz_decode_u21(mcode, ui3);
-
-    /* Out of bound voxels */
-    if(ui3[0] >= nvoxels[0]
-    || ui3[1] >= nvoxels[1]
-    || ui3[2] >= nvoxels[2])
-      continue;
-
-    /* Retrieve the voxel data from the caller */
-    xyz[0] = (size_t)ui3[0];
-    xyz[1] = (size_t)ui3[1];
-    xyz[2] = (size_t)ui3[2];
-    desc->get(xyz, vox.data, desc->context);
-    vox.mcode = mcode;
-
-    /* Register the voxel against the octree */
-    res = octree_builder_add_voxel(&bldr, &vox);
-    if(res != RES_OK) goto error;
-  }
-
-  res = octree_builder_finalize(&bldr, &oct->root, oct->root_attr);
-  if(res != RES_OK) goto error;
-
-  oct->nleaves = bldr.nleaves;
-  oct->depth = (size_t)(bldr.octree_depth - bldr.non_empty_lvl)
-    + 1 /* leaf level */;
-  ASSERT(bldr.octree_depth > bldr.non_empty_lvl);
-
-#ifndef NDEBUG
-  check_octree(&oct->buffer, oct->root);
-#endif
-  /* Setup the octree AABB in world space */
-  d3_set(oct->lower, lower);
-  d3_set(oct->upper, upper);
-
-  /* Compute the voxel size in world space */
-  vox_sz[0] = (upper[0] - lower[0]) / (double)nvoxels[0];
-  vox_sz[1] = (upper[1] - lower[1]) / (double)nvoxels[1];
-  vox_sz[2] = (upper[2] - lower[2]) / (double)nvoxels[2];
-
-  /* Compute the octree AABB in world space */
-  d3_set(oct->oclow, lower);
-  oct->ocupp[0] = (double)oct->definition * vox_sz[0];
-  oct->ocupp[1] = (double)oct->definition * vox_sz[1];
-  oct->ocupp[2] = (double)oct->definition * vox_sz[2];
-
-  /* Adjust the octree definition with respect to the octree depth since finest
-   * levels might be removed by the build process due to the merging process */
-  oct->definition = oct->definition / (1u << bldr.non_empty_lvl);
-
-exit:
-  if(vox.data) MEM_RM(dev->allocator, vox.data);
-  if(out_oct) *out_oct = oct;
-  return res;
-error:
-  if(oct) {
-    HTVOX(octree_ref_put(oct));
-    oct = NULL;
-  }
-  goto exit;
-}
-
-res_T
-htvox_octree_ref_get(struct htvox_octree* oct)
-{
-  if(!oct) return RES_BAD_ARG;
-  ref_get(&oct->ref);
-  return RES_OK;
-}
-
-res_T
-htvox_octree_ref_put(struct htvox_octree* oct)
-{
-  if(!oct) return RES_BAD_ARG;
-  ref_put(&oct->ref, octree_release);
-  return RES_OK;
-}
-
-res_T
-htvox_octree_get_desc
-  (const struct htvox_octree* oct, struct htvox_octree_desc* desc)
-{
-  if(!oct || !desc) return RES_BAD_ARG;
-  d3_set(desc->lower, oct->lower);
-  d3_set(desc->upper, oct->upper);
-  desc->nleaves = oct->nleaves;
-  desc->nvoxels =  absolute_attr_index(&oct->buffer, oct->buffer.attr_head)
-    + 1; /* Root node */
-  desc->depth = oct->depth;
-  return RES_OK;
-}
-
-res_T
-htvox_octree_for_each_leaf
-  (struct htvox_octree* oct,
-   void (*func)
-      (const void* val,
-       const size_t ileaf,
-       const double low[3],
-       const double upp[3],
-       void* ctx),
-   void* ctx)
-{
-  struct stack_entry {
-    struct octree_index inode;
-    double low[3];
-    double upp[3];
-  } stack[OCTREE_DEPTH_MAX*8];
-  int istack;
-  size_t ileaf = 0;
-
-  if(!oct || !func) return RES_BAD_ARG;
-
-  stack[0].inode = oct->root;
-  stack[0].low[0] = oct->oclow[0];
-  stack[0].low[1] = oct->oclow[1];
-  stack[0].low[2] = oct->oclow[2];
-  stack[0].upp[0] = oct->ocupp[0];
-  stack[0].upp[1] = oct->ocupp[1];
-  stack[0].upp[2] = oct->ocupp[2];
-  istack = 1;
-
-  do {
-    const struct stack_entry entry = stack[--istack];
-    struct octree_xnode* node;
-    int ichild;
-
-    node = octree_buffer_get_node(&oct->buffer, entry.inode);
-
-    FOR_EACH(ichild, 0, 8) {
-      const uint8_t ichild_flag = (uint8_t)BIT(ichild);
-      double half_sz[3]; /* Half size of the current node */
-      double mid[3]; /* Middle point of the current node */
-      double upp[3];
-      double low[3];
-
-      if((node->is_valid & ichild_flag) == 0) continue; /* Empty node */
-
-      half_sz[0] = (entry.upp[0] - entry.low[0])*0.5;
-      half_sz[1] = (entry.upp[1] - entry.low[1])*0.5;
-      half_sz[2] = (entry.upp[2] - entry.low[2])*0.5;
-      mid[0] = entry.low[0] + half_sz[0];
-      mid[1] = entry.low[1] + half_sz[1];
-      mid[2] = entry.low[2] + half_sz[2];
-      low[0] = ichild&4 ? mid[0] : entry.low[0];
-      low[1] = ichild&2 ? mid[1] : entry.low[1];
-      low[2] = ichild&1 ? mid[2] : entry.low[2];
-      upp[0] = low[0] + half_sz[0];
-      upp[1] = low[1] + half_sz[1];
-      upp[2] = low[2] + half_sz[2];
-
-      if(node->is_leaf & ichild_flag) {
-        struct octree_index iattr;
-        const void* val;
-
-        iattr = octree_buffer_get_child_attr_index
-          (&oct->buffer, entry.inode, ichild);
-        val = octree_buffer_get_attr(&oct->buffer, iattr);
-
-        ASSERT(upp[0] <= oct->upper[0]);
-        ASSERT(upp[1] <= oct->upper[1]);
-        ASSERT(upp[2] <= oct->upper[2]);
-        ASSERT(low[0] >= oct->lower[0]);
-        ASSERT(low[1] >= oct->lower[1]);
-        ASSERT(low[2] >= oct->lower[2]);
-
-        func(val, ileaf, low, upp, ctx);
-        ileaf++;
-      } else {
-        struct stack_entry* top = stack + istack;
-
-        top->inode = octree_buffer_get_child_node_index
-          (&oct->buffer, entry.inode, ichild);
-        top->low[0] = low[0];
-        top->low[1] = low[1];
-        top->low[2] = low[2];
-        top->upp[0] = upp[0];
-        top->upp[1] = upp[1];
-        top->upp[2] = upp[2];
-        ++istack;
-      }
-    }
-  } while(istack);
-
-  return RES_OK;
-}
-
-res_T
-htvox_octree_at
-  (struct htvox_octree* oct,
-   const double position[3],
-   int (*filter)  /* Filter function. May be NULL */
-    (const struct htvox_voxel* voxel,
-     const double position[3],
-     void* ctx),
-   void* context, /* Client data sent as the last argument of the filter func */
-   struct htvox_voxel* voxel)
-{
-  struct octree_index inode;
-  struct octree_index iattr;
-  struct htvox_voxel vox = HTVOX_VOXEL_NULL;
-  double scale_exp2;
-  double low[3];
-  double pos[3];
-  res_T res = RES_OK;
-
-  if(!oct || !position || !voxel) {
-    res = RES_BAD_ARG;
-    goto error;
-  }
-
-  *voxel = HTVOX_VOXEL_NULL;
-
-  /* The position is outside the octree */
-  if(position[0] > oct->upper[0] || position[0] < oct->lower[0]
-  || position[1] > oct->upper[1] || position[1] < oct->lower[1]
-  || position[2] > oct->upper[2] || position[2] < oct->lower[2]) {
-    goto exit;
-  }
-
-  /* Transform the position in the normalized octree space,
-   * i.e. octree lies in [0, 1]^2 */
-  pos[0] = (position[0] - oct->oclow[0]) / (oct->ocupp[0] - oct->oclow[0]);
-  pos[1] = (position[1] - oct->oclow[1]) / (oct->ocupp[1] - oct->oclow[1]);
-  pos[2] = (position[2] - oct->oclow[2]) / (oct->ocupp[2] - oct->oclow[2]);
-
-  /* Initialized the lower left corner of the current node */
-  low[0] = 0;
-  low[1] = 0;
-  low[2] = 0;
-
-  /* Root voxel */
-  vox.depth = 0;
-  vox.is_leaf = 0;
-  if(filter) {
-    vox.data = oct->root_attr;
-    vox.id = absolute_attr_index(&oct->buffer, oct->buffer.attr_head);
-    if(!filter(&vox, position, context)) { *voxel = vox; goto exit; }
-  }
-
-  scale_exp2 = 0.5;
-  inode = oct->root;
-  for(;;) {
-    struct octree_xnode* node = octree_buffer_get_node(&oct->buffer, inode);
-    int ichild;
-    uint8_t ichild_flag;
-    double mid[3];
-
-    ++vox.depth;
-
-    /* Compute the middle point of the node */
-    mid[0] = low[0] + scale_exp2;
-    mid[1] = low[1] + scale_exp2;
-    mid[2] = low[2] + scale_exp2;
-
-    /* Define the child of the current node into which pos `lies' and compute
-     * its lower left corner */
-    ichild = 0;
-    if(pos[0] > mid[0]) { ichild |= 4; low[0] = mid[0]; }
-    if(pos[1] > mid[1]) { ichild |= 2; low[1] = mid[1]; }
-    if(pos[2] > mid[2]) { ichild |= 1; low[2] = mid[2]; }
-
-    ichild_flag = (uint8_t)BIT(ichild);
-    if((node->is_valid & ichild_flag) == 0) break; /* Empty node */
-
-    vox.is_leaf = (node->is_leaf & ichild_flag) != 0;
-    if(filter || vox.is_leaf) {
-      iattr = octree_buffer_get_child_attr_index(&oct->buffer, inode, ichild);
-      vox.data = octree_buffer_get_attr(&oct->buffer, iattr);
-      vox.id = absolute_attr_index(&oct->buffer, iattr);
-      vox.is_leaf = (node->is_leaf & ichild_flag) != 0;
-      if(vox.is_leaf || !filter(&vox, position, context)) {
-        *voxel = vox;
-        break;
-      }
-    }
-
-    inode = octree_buffer_get_child_node_index(&oct->buffer, inode, ichild);
-    scale_exp2 *= 0.5;
-  }
-
-exit:
-  return res;
-error:
-  goto exit;
-}
-
diff --git a/src/htvox_octree.h b/src/htvox_octree.h
@@ -1,40 +0,0 @@
-/* Copyright (C) CNRS 2018
- *
- * This program is free software: you can redistribute it and/or modify
- * it under the terms of the GNU General Public License as published by
- * the Free Software Foundation, either version 3 of the License, or
- * (at your option) any later version.
- *
- * This program is distributed in the hope that it will be useful,
- * but WITHOUT ANY WARRANTY; without even the implied warranty of
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
- * GNU General Public License for more details.
- *
- * You should have received a copy of the GNU General Public License
- * along with this program. If not, see <http://www.gnu.org/licenses/>. */
-
-#ifndef HTVOX_OCTREE_H
-#define HTVOX_OCTREE_H
-
-#include "htvox_octree_buffer.h"
-#include <rsys/ref_count.h>
-
-struct htvox_octree {
-  size_t definition; /* #voxels along the X, Y and Z axis */
-
-  double lower[3], upper[3]; /* Submitted World space AABB of the octree  */
-  double oclow[3], ocupp[3]; /* Adjusted World space AABB of the octree */
-
-  struct octree_buffer buffer;
-  struct octree_index root; /* Index toward the children of the root */
-  ALIGN(16) char root_attr[HTVOX_MAX_SIZEOF_VOXEL]; /* Attribute of the root */
-
-  size_t nleaves; /* #leaves */
-  size_t depth; /* Maximum depth of the octree */
-
-  struct htvox_device* dev;
-  ref_T ref;
-};
-
-#endif /* HTVOX_OCTREE_H */
-
diff --git a/src/htvox_octree_buffer.c b/src/htvox_octree_buffer.c
@@ -1,198 +0,0 @@
-/* Copyright (C) CNRS 2018
- *
- * This program is free software: you can redistribute it and/or modify
- * it under the terms of the GNU General Public License as published by
- * the Free Software Foundation, either version 3 of the License, or
- * (at your option) any later version.
- *
- * This program is distributed in the hope that it will be useful,
- * but WITHOUT ANY WARRANTY; without even the implied warranty of
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
- * GNU General Public License for more details.
- *
- * You should have received a copy of the GNU General Public License
- * along with this program. If not, see <http://www.gnu.org/licenses/>. */
-
-#include "htvox_octree_buffer.h"
-
-#include <rsys/mem_allocator.h>
-
-#include <unistd.h>
-
-/*******************************************************************************
- * Helper functions
- ******************************************************************************/
-static INLINE res_T
-ensure_allocated_nodes(struct octree_buffer* buf, const size_t nnodes)
-{
-  char* node_page = NULL;
-  size_t nnode_pages = 0;
-  res_T res = RES_OK;
-  ASSERT(buf);
-
-  if(buf->node_head.ipage != OCTREE_INDEX_NULL.ipage
-  && buf->node_head.inode + nnodes <= buf->pagesize/sizeof(struct octree_xnode))
-    goto exit;
-
-  nnode_pages = darray_page_size_get(&buf->node_pages);
-  if(nnode_pages > UINT32_MAX) { res = RES_MEM_ERR; goto error; }
-  ASSERT(nnode_pages == buf->node_head.ipage + 1);
-
-  /* Alloc and register a node page containing the node and the far indices */
-  node_page = MEM_ALLOC(buf->allocator, buf->pagesize);
-  if(!node_page) { res = RES_MEM_ERR; goto error; }
-  res = darray_page_push_back(&buf->node_pages, &node_page);
-  if(res != RES_OK) goto error;
-
-  buf->node_head.inode = 0;
-  buf->node_head.ipage = (uint32_t)nnode_pages;
-
-exit:
-  return res;
-error:
-  if(node_page) MEM_RM(buf->allocator, node_page);
-  CHK(darray_page_resize(&buf->node_pages, nnode_pages) == RES_OK);
-  goto exit;
-}
-
-static INLINE res_T
-ensure_allocated_attrs(struct octree_buffer* buf, const size_t nattrs)
-{
-  char* attr_page = NULL;
-  size_t nattr_pages = 0;
-  res_T res = RES_OK;
-  ASSERT(buf);
-
-  if(buf->attr_head.ipage != OCTREE_INDEX_NULL.ipage
-  && buf->attr_head.inode + nattrs <= buf->pagesize/buf->voxsize)
-    goto exit;
-
-  nattr_pages = darray_page_size_get(&buf->attr_pages);
-  if(nattr_pages > UINT32_MAX) { res = RES_MEM_ERR; goto error; }
-  ASSERT(nattr_pages == buf->attr_head.ipage + 1);
-
-  /* Alloc and register a attr page */
-  attr_page = MEM_ALLOC(buf->allocator, buf->pagesize);
-  if(!attr_page) { res = RES_MEM_ERR; goto error; }
-  res = darray_page_push_back(&buf->attr_pages, &attr_page);
-  if(res != RES_OK) goto error;
-
-  buf->attr_head.inode = 0;
-  buf->attr_head.ipage = (uint32_t)nattr_pages;
-
-exit:
-  return res;
-error:
-  if(attr_page) MEM_RM(buf->allocator, attr_page);
-  CHK(darray_page_resize(&buf->attr_pages, nattr_pages) == RES_OK);
-  goto exit;
-}
-
-/*******************************************************************************
- * Local functions
- ******************************************************************************/
-void
-octree_buffer_init
-  (struct mem_allocator* allocator,
-   const size_t voxel_size,
-   struct octree_buffer* buf)
-{
-  ASSERT(buf && allocator);
-  memset(buf, 0, sizeof(struct octree_buffer));
-  buf->pagesize = (size_t)sysconf(_SC_PAGESIZE);
-  buf->voxsize = voxel_size;
-  darray_page_init(allocator, &buf->node_pages);
-  darray_page_init(allocator, &buf->attr_pages);
-  buf->node_head = OCTREE_INDEX_NULL;
-  buf->attr_head = OCTREE_INDEX_NULL;
-  buf->allocator = allocator;
-  CHK(buf->voxsize <= buf->pagesize);
-}
-
-void
-octree_buffer_release(struct octree_buffer* buf)
-{
-  ASSERT(buf);
-  octree_buffer_clear(buf);
-  darray_page_release(&buf->node_pages);
-  darray_page_release(&buf->attr_pages);
-}
-
-res_T
-octree_buffer_alloc_nodes
-  (struct octree_buffer* buf,
-   const size_t nnodes,
-   struct octree_index* first_node)
-{
-  res_T res = RES_OK;
-  ASSERT(buf && first_node);
-
-  if(nnodes > buf->pagesize / sizeof(struct octree_xnode))
-    return RES_MEM_ERR;
-
-  res = ensure_allocated_nodes(buf, nnodes);
-  if(res != RES_OK) return res;
-
-  *first_node = buf->node_head;
-  buf->node_head.inode = (uint16_t)(buf->node_head.inode + nnodes);
-  return RES_OK;
-
-}
-
-res_T
-octree_buffer_alloc_attrs
-  (struct octree_buffer* buf,
-   const size_t nattrs,
-   struct octree_index* first_attr)
-{
-  res_T res = RES_OK;
-  ASSERT(buf && first_attr);
-
-  if(nattrs > buf->pagesize / buf->voxsize) return RES_MEM_ERR;
-
-  res = ensure_allocated_attrs(buf, nattrs);
-  if(res != RES_OK) return res;
-
-  *first_attr = buf->attr_head;
-  buf->attr_head.inode = (uint16_t)(buf->attr_head.inode + nattrs);
-  return RES_OK;
-}
-
-res_T
-octree_buffer_alloc_far_index
-  (struct octree_buffer* buf,
-   struct octree_index* id)
-{
-  size_t remaining_size;
-  size_t skipped_nnodes;
-  STATIC_ASSERT(sizeof(struct octree_index) >= sizeof(struct octree_xnode),
-    Unexpected_type_size);
-
-  remaining_size = buf->pagesize - buf->node_head.inode*sizeof(struct octree_xnode);
-
-  /* Not enough memory in the current page */
-  if(sizeof(struct octree_index) > remaining_size) return RES_MEM_ERR;
-
-  *id = buf->node_head;
-  skipped_nnodes = sizeof(struct octree_index) / sizeof(struct octree_xnode);
-  buf->node_head.inode = (uint16_t)(buf->node_head.inode + skipped_nnodes);
-  return RES_OK;
-}
-
-void
-octree_buffer_clear(struct octree_buffer* buf)
-{
-  size_t i;
-  ASSERT(buf);
-  FOR_EACH(i, 0, darray_page_size_get(&buf->node_pages)) {
-    MEM_RM(buf->allocator, darray_page_data_get(&buf->node_pages)[i]);
-  }
-  FOR_EACH(i, 0, darray_page_size_get(&buf->attr_pages)) {
-    MEM_RM(buf->allocator, darray_page_data_get(&buf->attr_pages)[i]);
-  }
-  darray_page_purge(&buf->node_pages);
-  darray_page_purge(&buf->attr_pages);
-  buf->node_head = OCTREE_INDEX_NULL;
-  buf->attr_head = OCTREE_INDEX_NULL;
-}
-
diff --git a/src/htvox_octree_buffer.h b/src/htvox_octree_buffer.h
@@ -1,239 +0,0 @@
-/* Copyright (C) CNRS 2018
- *
- * This program is free software: you can redistribute it and/or modify
- * it under the terms of the GNU General Public License as published by
- * the Free Software Foundation, either version 3 of the License, or
- * (at your option) any later version.
- *
- * This program is distributed in the hope that it will be useful,
- * but WITHOUT ANY WARRANTY; without even the implied warranty of
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
- * GNU General Public License for more details.
- *
- * You should have received a copy of the GNU General Public License
- * along with this program. If not, see <http://www.gnu.org/licenses/>. */
-
-#ifndef HTVOX_OCTREE_BUFFER_H
-#define HTVOX_OCTREE_BUFFER_H
-
-#include "htvox_c.h"
-#include <rsys/dynamic_array.h>
-
-/*
- * Buffer containing the data of an octree. These data are partitioned in fixed
- * size memory pages whose capacity is defined on buffer initialisation with
- * respect to the page size of the system.
- *
- * The children of a node are stored consecutively into a page. The parent node
- * directly references its first valid children excepted if they lie in a
- * different page. In this case, the node references a `struct octree_index',
- * stored into the same page, that defines the absolute position of its first
- * valid child into the whole list of node pages
- *
- * The data of the nodes are stored in their own memory pages. The attribs of
- * the children of a node are stored consecutively into a page. If the page
- * identifier of the attribs is the same of the page into which their parent
- * node lies, then the node saves the index toward the first valid attrib into
- * the page of attribs. In the other case, the node references a `struct
- * octree_index', stored into the same page of the node, that defines the
- * absolute position of its first valid attrib into the buffer of attribs.
- */
-
-#define OCTREE_XNODE_FLAG_FAR_INDEX (1u<<15)
-#define OCTREE_XNODE_MAX_CHILDREN_OFFSET (OCTREE_XNODE_FLAG_FAR_INDEX-1)
-#define OCTREE_XNODE_MASK OCTREE_XNODE_MAX_CHILDREN_OFFSET
-
-struct octree_xnode {
-  /* Offset to retrieve the children. If the OCTREE_XNODE_FLAG_FAR_INDEX bit is
-   * not set, the children are stored in the same page at the position `offset
-   * & OCTREE_XNODE_MASK'. If OCTREE_XNODE_FLAG_FAR_INDEX is set, `offset &
-   * OCTREE_XNODE_MASK' reference an octree_index toward the node children */
-  uint16_t node_offset;
-  uint16_t attr_offset;
-  uint8_t is_valid; /* Mask defining if the children are valid */
-  uint8_t is_leaf; /* Mask defining if the children are leaves */
-  uint16_t dummy__; /* Ensure a size of 8 Bytes */
-};
-STATIC_ASSERT(sizeof(struct octree_xnode) == 8,
-  Unexpected_sizeof_octree_xnode);
-
-#define OCTREE_INDEX_IPAGE_MAX UINT32_MAX
-#define OCTREE_INDEX_INODE_MAX UINT16_MAX
-
-struct octree_index {
-  uint32_t ipage; /* Identifier of the page */
-  uint16_t inode; /* Identifier of the node in the page */
-  uint16_t dummy__; /* Padding to ensure the octree index is 8 bytes lenght */
-};
-STATIC_ASSERT(sizeof(struct octree_index) == 8,
-  Unexpected_sizeof_octree_index);
-#define OCTREE_INDEX_NULL__ {UINT32_MAX, UINT16_MAX, UINT16_MAX}
-static const struct octree_index OCTREE_INDEX_NULL = OCTREE_INDEX_NULL__;
-#define OCTREE_INDEX_EQ(A, B) ((A)->inode==(B)->inode && (A)->ipage==(B)->ipage)
-
-/* Define the dynamic array of pages */
-#define DARRAY_NAME page
-#define DARRAY_DATA char*
-#include <rsys/dynamic_array.h>
-
-struct octree_buffer {
-  size_t pagesize; /* Memory page size in bytes */
-  size_t voxsize; /* Memory size of a voxel in bytes */
-
-  struct darray_page node_pages; /* List of pages storing nodes */
-  struct darray_page attr_pages; /* List of pages storing node attributes */
-  struct octree_index node_head; /* Index of the next valid node */
-  struct octree_index attr_head; /* Index of the next valid attr */
-
-  struct mem_allocator* allocator;
-};
-
-extern LOCAL_SYM void
-octree_buffer_init
-  (struct mem_allocator* allocator,
-   const size_t sizeof_voxel, /* Size in bytes of a voxel */
-   struct octree_buffer* buf);
-
-extern LOCAL_SYM void
-octree_buffer_release
-  (struct octree_buffer* buf);
-
-extern LOCAL_SYM res_T
-octree_buffer_alloc_nodes
-  (struct octree_buffer* buf,
-   const size_t nnodes,
-   struct octree_index* first_node); /* Index toward the 1st allocated node */
-
-extern LOCAL_SYM res_T
-octree_buffer_alloc_attrs
-  (struct octree_buffer* buf,
-   const size_t nattrs,
-   struct octree_index* first_attr); /* Index toward the 1st allocated attrib */
-
-/* Allocate an octree_index in the current buffer page. Return RES_MEM_ERR if
- * the node index cannot be allocated in the current page. In this case one
- * have to alloc new nodes */
-extern LOCAL_SYM res_T
-octree_buffer_alloc_far_index
-  (struct octree_buffer* buf,
-   struct octree_index* id); /* Index toward the allocated far index */
-
-extern LOCAL_SYM void
-octree_buffer_clear
-  (struct octree_buffer* buf);
-
-static FINLINE int
-octree_buffer_is_empty(const struct octree_buffer* buf)
-{
-  ASSERT(buf);
-  return darray_page_size_get(&buf->node_pages) == 0;
-}
-
-static FINLINE struct octree_xnode*
-octree_buffer_get_node
-  (struct octree_buffer* buf,
-   const struct octree_index id)
-{
-  char* mem;
-  ASSERT(buf && id.inode < buf->pagesize/sizeof(struct octree_xnode));
-  ASSERT(id.ipage < darray_page_size_get(&buf->node_pages));
-  mem = darray_page_data_get(&buf->node_pages)[id.ipage];
-  mem += id.inode * sizeof(struct octree_xnode);
-  return (struct octree_xnode*)mem;
-}
-
-static FINLINE void*
-octree_buffer_get_attr
-  (struct octree_buffer* buf,
-   const struct octree_index id)
-{
-  char* mem;
-  ASSERT(buf && id.inode < buf->pagesize/buf->voxsize);
-  ASSERT(id.ipage < darray_page_size_get(&buf->attr_pages));
-  mem = darray_page_data_get(&buf->attr_pages)[id.ipage];
-  mem += id.inode * buf->voxsize;
-  return mem;
-}
-
-static FINLINE struct octree_index*
-octree_buffer_get_far_index
-  (struct octree_buffer* buf,
-   const struct octree_index id)
-{
-  char* mem;
-  ASSERT(buf && id.inode < buf->pagesize/sizeof(struct octree_xnode));
-  ASSERT(id.ipage < darray_page_size_get(&buf->node_pages));
-  mem = darray_page_data_get(&buf->node_pages)[id.ipage];
-  mem += id.inode * sizeof(struct octree_xnode);
-  return (struct octree_index*)mem;
-}
-
-static FINLINE struct octree_index
-octree_buffer_get_child_node_index
-  (struct octree_buffer* buf,
-   const struct octree_index id,
-   const int ichild) /* in [0, 7] */
-{
-  struct octree_index child_id = OCTREE_INDEX_NULL;
-  struct octree_xnode* node = NULL;
-  uint16_t offset;
-  const int ichild_flag = BIT(ichild);
-  int ichild_off;
-  uint8_t mask;
-
-  ASSERT(ichild >= 0 && ichild < 8 && buf);
-
-  node = octree_buffer_get_node(buf, id);
-  mask = (uint8_t)(node->is_valid & ~node->is_leaf);
-  ASSERT(mask & ichild_flag);
-
-  /* Compute the child offset from the first child node */
-  ichild_off = popcount((uint8_t)((ichild_flag-1) & (int)mask));
-
-  offset = node->node_offset & OCTREE_XNODE_MASK;
-  if(!(node->node_offset & OCTREE_XNODE_FLAG_FAR_INDEX)) {
-    child_id.ipage = id.ipage;
-    child_id.inode = (uint16_t)(offset + ichild_off);
-  } else {
-    char* mem = darray_page_data_get(&buf->node_pages)[id.ipage];
-    child_id = *(struct octree_index*)(mem+offset*(sizeof(struct octree_xnode)));
-    child_id.inode = (uint16_t)(child_id.inode + ichild_off);
-  }
-  return child_id;
-}
-
-static FINLINE struct octree_index
-octree_buffer_get_child_attr_index
-  (struct octree_buffer* buf,
-   const struct octree_index id,
-   const int ichild) /* In [0, 7] */
-{
-  struct octree_index child_id = OCTREE_INDEX_NULL;
-  struct octree_xnode* node = NULL;
-  uint16_t offset;
-  const int ichild_flag = BIT(ichild);
-  int ichild_off;
-  uint8_t mask;
-
-  ASSERT(ichild >= 0 && ichild < 8 && buf);
-
-  node = octree_buffer_get_node(buf, id);
-  mask = node->is_valid;
-  ASSERT(mask & ichild_flag);
-
-  /* Compute the attr offset from the first child node */
-  ichild_off = popcount((uint8_t)((ichild_flag-1) & (int)mask));
-
-  offset = node->attr_offset & OCTREE_XNODE_MASK;
-  if(!(node->attr_offset & OCTREE_XNODE_FLAG_FAR_INDEX)) {
-    child_id.ipage = id.ipage;
-    child_id.inode = (uint16_t)(offset + ichild_off);
-  } else {
-    char* mem = darray_page_data_get(&buf->node_pages)[id.ipage];
-    child_id = *(struct octree_index*)(mem+offset*(sizeof(struct octree_xnode)));
-    child_id.inode = (uint16_t)(child_id.inode + ichild_off);
-  }
-  return child_id;
-}
-
-#endif /* HTVOX_OCTREE_BUFFER_H */
diff --git a/src/htvox_scene_trace_ray.c b/src/htvox_scene_trace_ray.c
@@ -1,228 +0,0 @@
-/* Copyright (C) CNRS 2018
- *
- * This program is free software: you can redistribute it and/or modify
- * it under the terms of the GNU General Public License as published by
- * the Free Software Foundation, either version 3 of the License, or
- * (at your option) any later version.
- *
- * This program is distributed in the hope that it will be useful,
- * but WITHOUT ANY WARRANTY; without even the implied warranty of
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
- * GNU General Public License for more details.
- *
- * You should have received a copy of the GNU General Public License
- * along with this program. If not, see <http://www.gnu.org/licenses/>. */
-
-#include "htvox.h"
-#include "htvox_scene.h"
-
-struct ray {
-  float org[3];
-  float range[2];
-  float ts[3]; /* 1 / -abs(dir) */
-  uint32_t octant_mask;
-};
-
-/*******************************************************************************
- * Helper functions
- ******************************************************************************/
-/* Return RES_BAD_OP if the ray cannot intersect the scene */
-static res_T
-setup_ray
-  (const struct scene* scn,
-   const double org[3],
-   const double dir[3],
-   const double range[2]
-   struct ray* ray)
-{
-  double rcp_ocsz[3]; /* Reciprocal size of the World space octree AABB */
-  ASSERT(scn);
-
-  if(range[0] >= range[1]) return RES_BAD_OP; /* Disabled ray */
-
-  /* Ray paralelle to an axis and that does not intersect the scene AABB */
-  if((!dir[0] && (org[0] < scn->lower[0] || org[0] > scn->upper[0]))
-  || (!dir[1] && (org[1] < scn->lower[1] || org[1] > scn->upper[1]))
-  || (!dir[2] && (org[2] < scn->lower[2] || org[2] > scn->upper[2]))) {
-    return RES_BAD_OP;
-  }
-
-  /* Compute reciprocal size of the world space octree AABB */
-  rcp_ocsz[0] = 1.0 / (scn->ocupp[0] - scn->oclow[0]);
-  rcp_ocsz[1] = 1.0 / (scn->ocupp[1] - scn->oclow[1]);
-  rcp_ocsz[2] = 1.0 / (scn->ocupp[2] - scn->oclow[2]);
-
-  /* Transform the ray origin in the [1, 2] space */
-  ray->org[0] = (float)((org[0] - scn->oclow[0]) * rcp_ocsz[0] + 1.f);
-  ray->org[1] = (float)((org[1] - scn->oclow[1]) * rcp_ocsz[1] + 1.f);
-  ray->org[2] = (float)((org[2] - scn->oclow[2]) * rcp_ocsz[2] + 1.f);
-
-  /* Transform the range in the normalized octree space */
-  ray->range[0] = (float)(range[0] * rcp_ocsz[0]);
-  ray->range[1] = (float)(range[1] * rcp_ocsz[1]);
-
-  /* Transform the direction in the normalized octree space */
-  ray->dir[0] = (float)(dir[0] * rcp_ocsz[0]);
-  ray->dir[1] = (float)(dir[1] * rcp_ocsz[1]);
-  ray->dir[2] = (float)(dir[2] * rcp_ocsz[2]);
-
-  /* Let a ray defines as org + t*dir and X the coordinate of an axis aligned
-   * plane. The ray intersects X at t = (X - org)/dir = (X - org) * ts; with ts
-   * = 1/dir. Note that one assume that dir is always negative. */
-  ray->ts[0] = (float)(1.0 / -fabs(dir[0]));
-  ray->ts[1] = (float)(1.0 / -fabs(dir[1]));
-  ray->ts[2] = (float)(1.0 / -fabs(dir[2]));
-
-  /* Mirror rays with position directions */
-  ray->octant_mask = 0;
-  if(ray_dir[0] > 0) { ray->octant_mask ^= 4; ray->org[0] = 3.f - ray->org[0]; }
-  if(ray_dir[1] > 0) { ray->octant_mask ^= 2; ray->org[0] = 3.f - ray->org[1]; }
-  if(ray_dir[2] > 0) { ray->octant_mask ^= 1; ray->org[0] = 3.f - ray->org[2]; }
-
-  return RES_OK;
-}
-
-static res_T
-trace_ray(const struct scene* scn, const struct* ray, struct htvox_hit* hit)
-{
-  #define SCALE_MAX 23 /* #mantisse bits */
-  struct octree_index inode;
-  float scale_exp2;
-  float t_min, tmax;
-  float low[3], upp[3]; /* AABB of the current node in normalized space */
-  float mid[3]; /* Middle point of the current node in normalized space */
-  float corner[3];
-  uint32_t scale; /* stack index */
-  uint32_t ichild;
-  ASSERT(scn && ray && hit);
-  
-  hit = HTVOX_HIT_NONE; /* Initialise the hit to "no intersection" */
-
-  f3_splat(lower, 1);
-  f3_splat(upper, 2);
-
-  /* Compute the min/max ray intersection with the octree AABB in normalized
-   * space. Note that in this space, the octree AABB is in [1, 2]^3 */
-  t_min =      (2 - ray->org[0]) * ray->ts[0];
-  t_min = MMAX((2 - ray->org[1]) * ray->ts[1], t_min);
-  t_min = MMAX((2 - ray->org[2]) * ray->ts[2], t_min);  
-  t_max =      (1 - ray->org[0]) * ray->ts[0];
-  t_max = MMIN((1 - ray->org[1]) * ray->ts[1], t_max);
-  t_max = MMIN((1 - ray->org[2]) * ray->ts[2], t_max);
-  t_min = MMAX(ray->range[0], t_min);
-  t_max = MMIN(ray->range[1], t_max);
-  if(t_min >= t_max) return RES_OK; /* No intersection */
-
-  /* Traverrsal initialisation */
-  inode = scn->root;
-  scale_exp2 = 0.5f;
-  scale = SCALE_MAX - 1;
-
-  /* Define the first child id and the position of its lower left corner in the
-   * normalized octree space, i.e. in [1, 2]^3 */
-  ichild = 0;
-  corner[0] = 1.f;
-  corner[1] = 1.f;
-  corner[2] = 1.f;
-  if((1.5f - ray->org[0])*ray->ts[0] > t_min) { ichild ^= 4; corner[0] = 1.5f; }
-  if((1.5f - ray->org[1])*ray->ts[1] > t_min) { ichild ^= 2; corner[1] = 1.5f; }
-  if((1.5f - ray->org[2])*ray->ts[2] > t_min) { ichild ^= 1; corner[2] = 1.5f; }
-
-  /* Octree traversal */
-  scale_max = scale + 1;
-  while(scale < scale_max) {
-    const struct octree_xnode* node;
-    float t_corner[3];
-    float t_max_corner;
-    float t_max_child;
-    uint32_t ichild_adjusted = ichild ^ ray->octant_mask;
-    uint32_t istep;
-
-    inode = octree_buffer_get_child_index
-      (&scn->buffer, inode, (int)ichild_adjusted);
-    node = octree_buffer_get_node(&scn->buffer, inode);
-
-    /* Compute the exit point of the ray in the current child node */
-    t_corner[0] = (corner[0] - ray->org[0])*ray->ts[0];
-    t_corner[1] = (corner[1] - ray->org[1])*ray->ts[1];
-    t_corner[2] = (corner[2] - ray->org[2])*ray->ts[2];
-    t_max_corner = MMIN(MMIN(t_corner[0], t_corner[1]), t_corner[2]);
-
-    /* Traverse the current child */
-    if(!OCTREE_XNODE_IS_EMPTY(node)
-    && t_min <= (t_max_child = MMIN(t_max, t_max_corner))) {
-      float scale_exp2_child;
-      float t_max_parent;
-      float t_center[3];
-
-      t_max_parent = t_max;
-      t_max = t_max_child;
-
-      if(OCTREE_XNODE_IS_LEAF(node)) break; /* Leaf node */
-
-      scale_exp2_child = scale_exp2 * 0.5f;
-
-      /* center = corner - scale_exp2_child =>
-       * t_center = ts*(corner + scale_exp2_child - org)
-       * t_center = t_corner + ts*scale_exp2_child 
-       * Anyway we perforrm the whole computation to avoid numerical issues */
-      t_center[0] = (corner[0] - ray->org[0] + scale_exp2_child) * ray->ts[0];
-      t_center[1] = (corner[1] - ray->org[1] + scale_exp2_child) * ray->ts[1];
-      t_center[2] = (corner[2] - ray->org[2] + scale_exp2_child) * ray->ts[2];
-
-      /* Push the parent node */
-      stack[scale].t_max = t_max_parent;
-      stack[scale].inode = inode;
-
-      /* Define the id and the lower left corner of the first child */
-      ichild = 0;
-      if(t_center[0] > t_min) { ichild ^= 4; corner[0] += scale_exp2_child; }
-      if(t_center[1] > t_min) { ichild ^= 2; corner[1] += scale_exp2_child; }
-      if(t_center[2] > t_min) { ichild ^= 1; corner[2] += scale_exp2_child; }
-
-      --scale;
-      scale_exp2 = scale_exp2_child;
-      continue;
-    }
-
-    /* Define the id and the lower left corner of the next child */
-    istep = 0;
-    if(t_corner[0] <= t_max_corner) { istep ^= 4; corner[0] -= scale_exp2; }
-    if(t_corner[1] <= t_max_corner) { istep ^= 2; corner[1] -= scale_exp2; }
-    if(t_corner[2] <= t_max_corner) { istep ^= 1; corner[2] -= scale_exp2; }
-    ichild ^= istep; 
-  }
-}
-
-/*******************************************************************************
- * Exported functions
- ******************************************************************************/
-res_T
-htvox_scene_trace_ray
-  (struct htvox_scene* scn,
-   const double org[3],
-   const double dir[3],
-   const double range[2],
-   struct htvox_hit* hit)
-{
-  struct ray ray;
-  res_T res = RES_OK;
-
-  if(!scn || !org || !range || !hit) {
-    res = RES_BAD_ARG;
-    goto error;
-  }
-
-  *hit = HTVOX_HIT_NULL;
-
-  res = setup_ray(scn, org, dir, range, &ray);
-  if(res == RES_BAD_OP) { /* The ray cannot intersect the scene. */
-    res = RES_OK;
-    goto exit;
-  }
-
-exit:
-  return res;
-error:
-  goto exit;
-}
diff --git a/src/svx.h b/src/svx.h
@@ -0,0 +1,196 @@
+/* Copyright (C) 2018 CNRS
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. */
+
+#ifndef SVX_H
+#define SVX_H
+
+#include <rsys/rsys.h>
+#include <float.h>
+
+/* Library symbol management */
+#if defined(SVX_SHARED_BUILD) /* Build shared library */
+  #define SVX_API extern EXPORT_SYM
+#elif defined(SVX_STATIC) /* Use/build static library */
+  #define SVX_API extern LOCAL_SYM
+#else /* Use shared library */
+  #define SVX_API extern IMPORT_SYM
+#endif
+
+/* Helper macro that asserts if the invocation of the svx function `Func'
+ * returns an error. One should use this macro on svx function calls for
+ * which no explicit error checking is performed */
+#ifndef NDEBUG
+  #define SVX(Func) ASSERT(svx_ ## Func == RES_OK)
+#else
+  #define SVX(Func) svx_ ## Func
+#endif
+
+/* Maximum memory size of a voxel */
+#define SVX_MAX_SIZEOF_VOXEL (sizeof(double)*16)
+
+struct svx_voxel {
+  const void* data; /* Data of the voxel */
+  size_t id; /* Indentifier of the voxel */
+  size_t depth; /* Depth of the voxel, in [0, svx_octree_desc.depth[ */
+  int is_leaf; /* Define if the voxel is a leaf into the hierarchy */
+};
+
+#define SVX_VOXEL_NULL__ { NULL, SIZE_MAX, SIZE_MAX, 0 }
+static const struct svx_voxel SVX_VOXEL_NULL = SVX_VOXEL_NULL__;
+
+#define SVX_VOXEL_NONE(Voxel) ((Voxel)->id == SVX_VOXEL_NULL.id)
+
+/* Descriptor of a voxel */
+struct svx_voxel_desc {
+  /* Retrieve the data of the voxels*/
+  void
+  (*get)
+    (const size_t xyz[3], /* Voxel coordinate in voxel space */
+     void* dst, /* Where to store data */
+     void* ctx); /* Pointer toward user data */
+
+  /* Merge the data of N voxels */
+  void
+  (*merge)
+    (void* dst, /* Merged data */
+     const void* voxels[], /* Data to merge */
+     const size_t nvoxels, /* #submitted data */
+     void* ctx); /* Pointer toward user data */
+
+  /* Check if the voxel's data can be merged */
+  int
+  (*challenge_merge)
+    (const void* voxels[], /* Data candidates to the merge */
+     const size_t nvoxels, /* #candidates */
+     void* ctx); /* Pointer toward user data */
+
+  void* context; /* Client side data sent as the last argument of the clbbs */
+  size_t size; /* Size in bytes of a voxel. Must be <= SVX_MAX_SIZEOF_VOXE */
+};
+
+#define SVX_VOXEL_DESC_NULL__ { NULL, NULL, NULL, NULL, 0 }
+static const struct svx_voxel_desc SVX_VOXEL_DESC_NULL =
+  SVX_VOXEL_DESC_NULL__;
+
+struct svx_octree_desc {
+  /* Submitted Axis Aligned Bounding Box */
+  double lower[3], upper[3];
+
+  size_t nleaves; /* #leaves */
+  size_t nvoxels; /* #voxels, i.e. #leaves + #parents */
+  size_t depth; /* Depth of the octree */
+};
+
+#define SVX_OCTREE_DESC_NULL__ \
+  {{DBL_MAX, DBL_MAX, DBL_MAX}, {-DBL_MAX,-DBL_MAX,-DBL_MAX}, 0, 0, 0}
+static const struct svx_octree_desc SVX_OCTREE_DESC_NULL =
+  SVX_OCTREE_DESC_NULL__;
+
+struct svx_hit {
+  double distance; /* Distance from the ray origin to the impacted voxel */
+  struct svx_voxel voxel; /* Intersected voxel */
+};
+
+#define SVX_HIT_NULL__ {DBL_MAX, SVX_VOXEL_NULL__}
+static const struct svx_hit SVX_HIT_NULL = SVX_HIT_NULL__;
+
+#define SVX_HIT_NONE(Hit) ((Hit)->distance >= DBL_MAX)
+
+/* Forward declaration of external data types */
+struct logger;
+struct mem_allocator;
+
+/* Forward declaration of opaque data types */
+struct svx_device;
+struct svx_octree;
+
+BEGIN_DECLS
+
+/*******************************************************************************
+ * Device
+ ******************************************************************************/
+SVX_API res_T
+svx_device_create
+  (struct logger* logger,
+   struct mem_allocator* allocator, /* NULL <=> use default allocator */
+   const int verbose, /* Verbosity level */
+   struct svx_device** svx);
+
+SVX_API res_T
+svx_device_ref_get
+  (struct svx_device* svx);
+
+SVX_API res_T
+svx_device_ref_put
+  (struct svx_device* svx);
+
+/*******************************************************************************
+ * Octree
+ ******************************************************************************/
+SVX_API res_T
+svx_octree_create
+  (struct svx_device* dev,
+   const double lower[3], /* Lower bound of the octree */
+   const double upper[3], /* Upper bound of the octree */
+   const size_t nvoxels[3], /* # voxels along the 3 axis */
+   const struct svx_voxel_desc* desc, /* Descriptor of a voxel */
+   struct svx_octree** octree);
+
+SVX_API res_T
+svx_octree_ref_get
+  (struct svx_octree* octree);
+
+SVX_API res_T
+svx_octree_ref_put
+  (struct svx_octree* octree);
+
+SVX_API res_T
+svx_octree_get_desc
+  (const struct svx_octree* octree,
+   struct svx_octree_desc* desc);
+
+SVX_API res_T
+svx_octree_for_each_leaf
+  (struct svx_octree* octree,
+   void (*functor)
+    (const void* val, /* Value of the voxel */
+     const size_t ileaf, /* Identifier of the leaf in [0, #leafs[ */
+     const double low[3], /* World space lower bound of the voxel */
+     const double upp[3], /* World space upper bound of the voxel */
+     void* ctx), /* Client data */
+   void* context); /* Client data sent as the last argument of the callback */
+
+SVX_API res_T
+svx_octree_trace_ray
+  (struct svx_octree* octree,
+   const double ray_origin[3],
+   const double ray_direction[3],
+   const double ray_range[2],
+   struct svx_hit* hit);
+
+SVX_API res_T
+svx_octree_at
+  (struct svx_octree* octree,
+   const double position[3],
+   int (*filter)  /* Filter function. May be NULL <=> traverse up to leaves */
+    (const struct svx_voxel* voxel,
+     const double position[3],
+     void* ctx),
+   void* context, /* Client data sent as the last argument of the filter func */
+   struct svx_voxel* voxel);
+
+#endif /* SVX_H */
+
+
diff --git a/src/svx_c.h b/src/svx_c.h
@@ -0,0 +1,75 @@
+/* Copyright (C) 2018 CNRS
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. */
+
+#ifndef SVX_C_H
+#define SVX_C_H
+
+#include <rsys/rsys.h>
+
+static INLINE uint64_t
+morton3D_encode_u21(const uint32_t u21)
+{
+  uint64_t u64 = u21 & ((1<<21) - 1);
+  ASSERT(u21 <= ((1 << 21) - 1));
+  u64 = (u64 | (u64 << 32)) & 0xFFFF00000000FFFF;
+  u64 = (u64 | (u64 << 16)) & 0x00FF0000FF0000FF;
+  u64 = (u64 | (u64 << 8))  & 0xF00F00F00F00F00F;
+  u64 = (u64 | (u64 << 4))  & 0x30C30C30C30C30C3;
+  u64 = (u64 | (u64 << 2))  & 0x9249249249249249;
+  return u64;
+}
+
+static INLINE uint32_t
+morton3D_decode_u21(const uint64_t u64)
+{
+  uint64_t tmp = (u64 & 0x9249249249249249);
+  tmp = (tmp | (tmp >> 2))  & 0x30C30C30C30C30C3;
+  tmp = (tmp | (tmp >> 4))  & 0xF00F00F00F00F00F;
+  tmp = (tmp | (tmp >> 8))  & 0x00FF0000FF0000FF;
+  tmp = (tmp | (tmp >> 16)) & 0xFFFF00000000FFFF;
+  tmp = (tmp | (tmp >> 32)) & 0x00000000FFFFFFFF;
+  ASSERT(tmp <= ((1<<21)-1));
+  return (uint32_t)tmp;
+}
+
+/* Count the number of bits set to 1 */
+static FINLINE int
+popcount(const uint8_t x)
+{
+  int n = x - ((x >> 1) & 0x55);
+  n = (n & 0x33) + ((n >> 2) & 0x33);
+  n = (n + (n >> 4)) & 0x0f;
+  return (n * 0x0101) >> 8;
+}
+
+static INLINE uint64_t
+morton_xyz_encode_u21(const uint32_t xyz[3])
+{
+  return (morton3D_encode_u21(xyz[0]) << 2)
+       | (morton3D_encode_u21(xyz[1]) << 1)
+       | (morton3D_encode_u21(xyz[2]) << 0);
+}
+
+static INLINE void
+morton_xyz_decode_u21(const uint64_t code, uint32_t xyz[3])
+{
+  ASSERT(xyz && code < ((1ul << 63)-1));
+  xyz[0] = (uint32_t)morton3D_decode_u21(code >> 2);
+  xyz[1] = (uint32_t)morton3D_decode_u21(code >> 1);
+  xyz[2] = (uint32_t)morton3D_decode_u21(code >> 0);
+}
+
+#endif /* SVX_C_H */
+
diff --git a/src/svx_device.c b/src/svx_device.c
@@ -0,0 +1,124 @@
+/* Copyright (C) 2018 CNRS
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. */
+
+#include "svx.h"
+#include "svx_device.h"
+
+#include <rsys/logger.h>
+#include <rsys/mem_allocator.h>
+
+/*******************************************************************************
+ * Helper functions
+ ******************************************************************************/
+static void
+log_msg
+  (const struct svx_device* dev,
+   const enum log_type stream,
+   const char* msg,
+   va_list vargs)
+{
+  ASSERT(dev && msg);
+  if(dev->verbose) {
+    res_T res; (void)res;
+    res = logger_vprint(dev->logger, stream, msg, vargs);
+    ASSERT(res == RES_OK);
+  }
+}
+
+static void
+device_release(ref_T* ref)
+{
+  struct svx_device* dev;
+  ASSERT(ref);
+  dev = CONTAINER_OF(ref, struct svx_device, ref);
+  MEM_RM(dev->allocator, dev);
+}
+
+/*******************************************************************************
+ * Exported functions
+ ******************************************************************************/
+res_T
+svx_device_create
+  (struct logger* log,
+   struct mem_allocator* mem_allocator,
+   const int verbose,
+   struct svx_device** out_dev)
+{
+  struct svx_device* dev = NULL;
+  struct mem_allocator* allocator = NULL;
+  struct logger* logger = NULL;
+  res_T res = RES_OK;
+
+  if(!out_dev) {
+    res = RES_BAD_ARG;
+    goto error;
+  }
+
+  allocator = mem_allocator ? mem_allocator : &mem_default_allocator;
+  logger = log ? log : LOGGER_DEFAULT;
+
+  dev = MEM_CALLOC(allocator, 1, sizeof(struct svx_device));
+  if(!dev) {
+    res = RES_MEM_ERR;
+    goto error;
+  }
+
+
+  ref_init(&dev->ref);
+  dev->allocator = allocator;
+  dev->logger = logger;
+  dev->verbose = verbose;
+
+exit:
+  if(out_dev) *out_dev = dev;
+  return res;
+error:
+  if(dev) {
+    SVX(device_ref_put(dev));
+    dev = NULL;
+  }
+  goto exit;
+}
+
+res_T
+svx_device_ref_get(struct svx_device* dev)
+{
+  if(!dev) return RES_BAD_ARG;
+  ref_get(&dev->ref);
+  return RES_OK;
+}
+
+res_T
+svx_device_ref_put(struct svx_device* dev)
+{
+  if(!dev) return RES_BAD_ARG;
+  ref_put(&dev->ref, device_release);
+  return RES_OK;
+}
+
+/*******************************************************************************
+ * Local function
+ ******************************************************************************/
+void
+log_err(const struct svx_device* dev, const char* msg, ...)
+{
+  va_list vargs_list;
+  ASSERT(dev && msg);
+
+  va_start(vargs_list, msg);
+  log_msg(dev, LOG_ERROR, msg, vargs_list);
+  va_end(vargs_list);
+}
+
diff --git a/src/svx_device.h b/src/svx_device.h
@@ -0,0 +1,39 @@
+/* Copyright (C) 2018 CNRS
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. */
+
+#ifndef SVX_DEVICE_H
+#define SVX_DEVICE_H
+
+#include <rsys/ref_count.h>
+
+struct svx_device {
+  int verbose; /* Verbosity level */
+  struct logger* logger;
+  struct mem_allocator* allocator;
+  ref_T ref;
+};
+
+extern LOCAL_SYM void
+log_err
+  (const struct svx_device* dev,
+   const char* fmt,
+  ...)
+#ifdef COMPILER_GCC
+  __attribute((format(printf, 2, 3)))
+#endif
+  ;
+
+#endif /* SVX_DEVICE_H */
+
diff --git a/src/svx_octree.c b/src/svx_octree.c
@@ -0,0 +1,820 @@
+/* Copyright (C) 2018 CNRS
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. */
+
+#include "svx.h"
+#include "svx_c.h"
+#include "svx_device.h"
+#include "svx_octree.h"
+
+#include <rsys/double3.h>
+#include <rsys/mem_allocator.h>
+
+#define OCTREE_DEPTH_MAX 16 /* Maximum depth of an octree */
+
+struct voxel {
+  uint64_t mcode; /* Morton code of the voxel */
+  void* data; /* Data of the voxel */
+};
+static const struct voxel VOXEL_NULL = {0, NULL};
+
+struct octree_node {
+  struct octree_index ichild_node; /* Index of the 1st child node */
+  struct octree_index ichild_attr; /* Index of the 1st child attr */
+  uint8_t is_valid; /* Mask defining whether the children are valid or not */
+  uint8_t is_leaf; /* Mask defining whether the children are leaves or not */
+  ALIGN(16) char data[8][SVX_MAX_SIZEOF_VOXEL]; /* Data of the leaves */
+};
+
+/* Stacked children of an octree node */
+struct stack {
+  struct octree_node nodes[8]; /* List of registered children nodes */
+  uint8_t mask; /* Mask of valid children nodes (0 = empty) */
+};
+
+struct octree_builder {
+  struct stack stacks[OCTREE_DEPTH_MAX];
+  struct octree_buffer* buffer;
+
+  const struct svx_voxel_desc* desc;
+  size_t nleaves; /* Number of emitted leaves */
+
+  int octree_depth; /* Maximum octree depth */
+  int non_empty_lvl; /* Index of the 1st non empty level */
+
+  uint64_t mcode; /* Morton code of the last registered voxels */
+};
+
+/*******************************************************************************
+ * Helper function
+ ******************************************************************************/
+#ifndef NDEBUG
+static void
+check_octree(struct octree_buffer* buf, const struct octree_index root)
+{
+  const struct octree_xnode* node;
+  int ichild;
+  ASSERT(buf);
+
+  node = octree_buffer_get_node(buf, root);
+  FOR_EACH(ichild, 0, 8) {
+    const int ichild_flag = BIT(ichild);
+    if((node->is_valid & ichild_flag) == 0) continue;
+
+    if(node->is_leaf & ichild_flag) {
+      struct octree_index iattr;
+      iattr = octree_buffer_get_child_attr_index(buf, root, ichild);
+      ASSERT(octree_buffer_get_attr(buf, iattr) != NULL);
+    } else {
+      struct octree_index child;
+      child = octree_buffer_get_child_node_index(buf, root, ichild);
+      check_octree(buf, child);
+    }
+  }
+}
+#endif
+
+static FINLINE size_t
+absolute_attr_index
+  (const struct octree_buffer* buf,
+   const struct octree_index index)
+{
+  ASSERT(buf);
+  return index.ipage * buf->pagesize/buf->voxsize + index.inode;
+
+}
+
+static INLINE int
+check_svx_voxel_desc(const struct svx_voxel_desc* desc)
+{
+  return desc
+      && desc->get
+      && desc->merge
+      && desc->challenge_merge
+      && desc->size > 0
+      && desc->size <= SVX_MAX_SIZEOF_VOXEL;
+}
+
+static INLINE void
+stack_clear(struct stack* stack)
+{
+  int inode;
+  FOR_EACH(inode, 0, 8) {
+    int ichild;
+    stack->nodes[inode].is_leaf = 0;
+    stack->nodes[inode].is_valid = 0;
+    stack->nodes[inode].ichild_node = OCTREE_INDEX_NULL;
+    stack->nodes[inode].ichild_attr = OCTREE_INDEX_NULL;
+    FOR_EACH(ichild, 0, 8) {
+      memset(stack->nodes[inode].data[ichild], 0, SVX_MAX_SIZEOF_VOXEL);
+    }
+  }
+  stack->mask = 0;
+}
+
+/* Build a parent octree node from the registered stack nodes */
+static INLINE void
+stack_setup_node
+  (struct stack* stack,
+   struct octree_node* node,
+   const struct svx_voxel_desc* desc)
+{
+  int ichild;
+  ASSERT(stack && node && check_svx_voxel_desc(desc));
+
+  node->ichild_node = OCTREE_INDEX_NULL;
+  node->ichild_attr = OCTREE_INDEX_NULL;
+  node->is_valid = stack->mask;
+  node->is_leaf = 0;
+
+  if(stack->mask == 0) return; /* Empty stack */
+
+  /* Try to merge the child's leaves */
+  FOR_EACH(ichild, 0, 8) {
+    const void* data[8];
+    const uint8_t ichild_flag = (uint8_t)BIT(ichild);
+    struct octree_node* child = stack->nodes + ichild;
+    int igrandchild;
+
+    /* Fetch the grandchildren data */
+    FOR_EACH(igrandchild, 0, 8) {
+      data[igrandchild] = child->data[igrandchild];
+    }
+
+    desc->merge(node->data[ichild], data, 8, desc->context);
+
+    if(child->is_leaf==0xFF && desc->challenge_merge(data, 8, desc->context)) {
+      /* The node becomes a leaf : the children does not exist anymore */
+      node->is_leaf |= ichild_flag;
+      stack->mask ^= ichild_flag;
+    }
+  }
+}
+
+static res_T
+stack_write
+  (struct stack* stack, /* Node to write */
+   struct octree_buffer* buf, /* Buffer where nodes are written */
+   struct octree_index* out_index, /* Index of the first written node */
+   size_t* out_nleaves) /* #writen leaves */
+{
+  struct octree_index nodes_id = OCTREE_INDEX_NULL;
+  struct octree_node* node = NULL;
+  size_t nleaves = 0;
+  int inode;
+  res_T res = RES_OK;
+  ASSERT(stack && buf && out_index && out_nleaves);
+
+  /* No registered nodes, this means that the nodes were merged in an higher
+   * level */
+  if(!stack->mask) goto exit;
+
+  /* Write the attrib of the children */
+  FOR_EACH(inode, 0, 8) {
+    char* data = NULL;
+    size_t nattrs = 0;
+    size_t nvoxs = 0;
+    int ichild = 0;
+
+    if((stack->mask & BIT(inode)) == 0) continue; /* Empty node */
+
+    node = stack->nodes + inode;
+
+    nattrs =  (size_t)popcount(node->is_valid);
+    ASSERT(nattrs > 0);
+
+    res = octree_buffer_alloc_attrs(buf, nattrs, &node->ichild_attr);
+    if(res != RES_OK) goto error;
+
+    data = octree_buffer_get_attr(buf, node->ichild_attr);
+    nvoxs = 0;
+    FOR_EACH(ichild, 0, 8) {
+      if(!(node->is_valid & BIT(ichild))) continue;
+      memcpy(data + nvoxs*buf->voxsize, node->data[ichild], buf->voxsize);
+      ++nvoxs;
+    }
+    ASSERT(nvoxs == nattrs);
+    nleaves += (size_t)popcount(node->is_leaf);
+  }
+
+  do {
+    struct octree_index index = OCTREE_INDEX_NULL;
+    struct octree_xnode* xnodes = NULL;
+    const size_t nnodes = (size_t)popcount(stack->mask);
+    size_t ixnode = 0;
+
+    /* Alloc the octree nodes */
+    res = octree_buffer_alloc_nodes(buf, (size_t)nnodes, &nodes_id);
+    if(res != RES_OK) goto error;
+    xnodes = octree_buffer_get_node(buf, nodes_id);
+
+    FOR_EACH(inode, 0, 8) {
+      uint16_t attr_offset = UINT16_MAX;
+      uint16_t node_offset = UINT16_MAX;
+
+      if((stack->mask & BIT(inode)) == 0) continue; /* Empty node */
+      node = stack->nodes + inode;
+
+      /* Setup the offset toward the children and children attribs */
+      if(node->ichild_node.ipage == nodes_id.ipage) {
+        node_offset = node->ichild_node.inode;
+        attr_offset = node->ichild_attr.inode;
+      }
+
+      /* The page id of the children is not the same as that of node */
+      if(node_offset > OCTREE_XNODE_MAX_CHILDREN_OFFSET) {
+        res = octree_buffer_alloc_far_index(buf, &index);
+        if(res != RES_OK) break;
+        *octree_buffer_get_far_index(buf, index) = node->ichild_node;
+        node_offset = OCTREE_XNODE_FLAG_FAR_INDEX | index.inode;
+      }
+
+      /* The page id of the attribs is not tthe same as that of node */
+      if(attr_offset > OCTREE_XNODE_FLAG_FAR_INDEX) {
+        res = octree_buffer_alloc_far_index(buf, &index);
+        if(res != RES_OK) break;
+        *octree_buffer_get_far_index(buf, index) = node->ichild_attr;
+        attr_offset = OCTREE_XNODE_FLAG_FAR_INDEX | index.inode;
+      }
+
+      xnodes[ixnode].node_offset = node_offset;
+      xnodes[ixnode].attr_offset = attr_offset;
+      xnodes[ixnode].is_valid = node->is_valid;
+      xnodes[ixnode].is_leaf = node->is_leaf;
+      ++ixnode;
+    }
+  /* inode < 8 <=> not enough memory in the current page. A far index could not
+   * be stored in the same page of its associated node. The write process was
+   * stoped. Rewrite the whole stacked nodes in a new page. */
+  } while(inode < 8);
+
+
+exit:
+  /* Return the index toward the first writen nodes */
+  *out_index = nodes_id;
+  *out_nleaves = nleaves;
+  return res;
+error:
+  goto exit;
+}
+
+static res_T
+octree_builder_init
+  (struct octree_builder* bldr,
+   const size_t definition,
+   const struct svx_voxel_desc* desc,
+   struct octree_buffer* buffer)
+{
+  int ilvl;
+  res_T res = RES_OK;
+  ASSERT(bldr && IS_POW2(definition) && check_svx_voxel_desc(desc));
+  memset(bldr, 0, sizeof(struct octree_builder));
+
+  /* Compute the maximum depth of the octree */
+  bldr->octree_depth = log2i((int)definition);
+  if(bldr->octree_depth > OCTREE_DEPTH_MAX) {
+    res = RES_MEM_ERR;
+    goto error;
+  }
+
+  /* Init the per octree level stack */
+  FOR_EACH(ilvl, 0, bldr->octree_depth) {
+    stack_clear(&bldr->stacks[ilvl]);
+  }
+
+  octree_buffer_clear(buffer);
+  bldr->nleaves = 0;
+  bldr->desc = desc;
+  bldr->buffer = buffer;
+  bldr->non_empty_lvl = bldr->octree_depth - 1;
+
+exit:
+  return res;
+error:
+  goto exit;
+}
+
+static res_T
+octree_builder_add_voxel
+  (struct octree_builder* bldr,
+   const struct voxel* vox)
+{
+  uint64_t mcode_xor;
+  int inode;
+  int ichild;
+  uint8_t ichild_flag;
+  res_T res = RES_OK;
+  ASSERT(bldr && vox && vox->mcode >= bldr->mcode);
+
+  /* Define if the bits in [4 .. 63] of the previous and the next Morton
+   * codes are the same */
+  mcode_xor = bldr->mcode ^ vox->mcode;
+
+  /* The next voxel is not in the current node */
+  if(mcode_xor >= 8) {
+    size_t ilvl;
+
+    inode = (bldr->mcode >> 3) & 7;
+    bldr->stacks[0].mask |= (uint8_t)BIT(inode);
+
+    /* Flush the stack of the octree level that does not contain the next
+     * voxel. The last octree level is actually the level that contains the
+     * root node while the penultimate describes the root node itself. These 2
+     * levels contain all the voxels and can be skipped */
+    FOR_EACH(ilvl, 0, bldr->octree_depth-2/*The 2 last leaves contain all voxels*/) {
+      struct octree_node* stack_node;
+      uint64_t mcode_max_lvl;
+      size_t nleaves;
+
+      /* Compute the maximum morton code value for the current octree level */
+      mcode_max_lvl = 8lu/*#children*/ * (1lu<<(3*(ilvl+1)))/*#voxels per children*/;
+
+      if(mcode_xor < mcode_max_lvl) break;
+
+      /* Retrieve the node index of the next level */
+      inode = (bldr->mcode >> (3*(ilvl+2))) & 7;
+
+      /* The next voxel is not in the ilvl^th stack. Setup the parent node of the
+       * nodes registered into the stack */
+      stack_node = &bldr->stacks[ilvl+1].nodes[inode];
+      stack_setup_node(&bldr->stacks[ilvl], stack_node, bldr->desc);
+      bldr->stacks[ilvl+1].mask |= (uint8_t)BIT(inode);
+
+      /* Write the nodes of the stack of the current octree level into the buf */
+      res = stack_write
+        (&bldr->stacks[ilvl], bldr->buffer, &stack_node->ichild_node, &nleaves);
+      if(res != RES_OK) goto error;
+
+      bldr->nleaves += nleaves;
+      if(nleaves) bldr->non_empty_lvl = MMIN(bldr->non_empty_lvl, (int)ilvl);
+
+      /* Reset the current stack */
+      stack_clear(&bldr->stacks[ilvl]);
+    }
+  }
+
+  /* Retrieve the index of the current voxel and of its parent */
+  ichild = vox->mcode & 7;
+  inode = (vox->mcode >> 3) & 7;
+  ichild_flag = (uint8_t)BIT(ichild);
+
+  /* Register the voxel */
+  memcpy(bldr->stacks[0].nodes[inode].data[ichild], vox->data, bldr->desc->size);
+  bldr->stacks[0].nodes[inode].is_valid |= ichild_flag;
+  bldr->stacks[0].nodes[inode].is_leaf |= ichild_flag;
+
+  /* Update morton code of the last registered voxel */
+  bldr->mcode = vox->mcode;
+
+exit:
+  return res;
+error:
+  goto exit;
+}
+
+static INLINE res_T
+octree_builder_finalize
+  (struct octree_builder* bldr,
+   struct octree_index* root_id,
+   void* root_data)
+{
+  const void* data[8];
+  size_t inode;
+  size_t nleaves;
+  int ilvl;
+  int ichild;
+  res_T res = RES_OK;
+  ASSERT(bldr);
+
+  inode = (bldr->mcode >> 3) & 7;
+  bldr->stacks[0].mask |= (uint8_t)BIT(inode);
+
+  /* Flush the stacked nodes */
+  FOR_EACH(ilvl, 0, bldr->octree_depth-1) {
+    struct octree_node* parent_node;
+
+    if(bldr->stacks[ilvl].mask == 0) continue;
+
+    /* Retrieve the node index of the next level */
+    inode = (bldr->mcode >> (3*(ilvl+2))) & 7;
+
+    /* Setup the parent node of the nodes registered into the current stack */
+    parent_node = &bldr->stacks[ilvl+1].nodes[inode]; /* Fetch the parent node */
+    stack_setup_node(&bldr->stacks[ilvl], parent_node, bldr->desc);
+    bldr->stacks[ilvl+1].mask |= (uint8_t)BIT(inode);
+
+    /* Write the stacked nodes of the current level */
+    res = stack_write
+      (&bldr->stacks[ilvl], bldr->buffer, &parent_node->ichild_node, &nleaves);
+    if(res != RES_OK) goto error;
+    bldr->nleaves += nleaves;
+  }
+
+  ilvl = bldr->octree_depth-1; /* Root level */
+
+  /* Write the root node */
+  res = stack_write(&bldr->stacks[ilvl], bldr->buffer, root_id, &nleaves);
+  if(res != RES_OK) goto error;
+  bldr->nleaves += nleaves;
+
+  /* Setup the root attribs */
+  ASSERT(bldr->stacks[ilvl].mask == 1); /* Only the root node is active */
+  FOR_EACH(ichild, 0, 8) {
+    data[ichild] = bldr->stacks[ilvl].nodes[0].data[ichild];
+  }
+  bldr->desc->merge(root_data, data, 8, bldr->desc->context);
+
+exit:
+  return res;
+error:
+  goto exit;
+}
+
+static void
+octree_release(ref_T* ref)
+{
+  struct svx_octree* oct;
+  struct svx_device* dev;
+  ASSERT(ref);
+  oct = CONTAINER_OF(ref, struct svx_octree, ref);
+  octree_buffer_release(&oct->buffer);
+  dev = oct->dev;
+  MEM_RM(dev->allocator, oct);
+  SVX(device_ref_put(dev));
+}
+
+/*******************************************************************************
+ * Exported functions
+ ******************************************************************************/
+res_T
+svx_octree_create
+  (struct svx_device* dev,
+   const double lower[3], /* Lower bound of the octree */
+   const double upper[3], /* Upper bound of the octree */
+   const size_t nvoxels[3], /* # voxels along the 3 axis */
+   const struct svx_voxel_desc* desc, /* Descriptor of a voxel */
+   struct svx_octree** out_oct)
+{
+  struct svx_octree* oct = NULL;
+  double vox_sz[3]; /* World space size of a voxel */
+  struct octree_builder bldr;
+  struct voxel vox = VOXEL_NULL;
+  uint64_t mcode_max;
+  uint64_t mcode;
+  res_T res = RES_OK;
+
+  if(!dev || !lower || !upper || !nvoxels
+  || !check_svx_voxel_desc(desc) || !out_oct) {
+    res = RES_BAD_ARG;
+    goto error;
+  }
+  if(lower[0] >= upper[0]
+  || lower[1] >= upper[1]
+  || lower[2] >= upper[2]) {
+    log_err(dev,
+      "%s: the submitted AABB is degenerated.\n"
+      "\tlower={%g, %g, %g}, upper={%g, %g, %g}.\n",
+      FUNC_NAME, SPLIT3(lower), SPLIT3(upper));
+    res = RES_BAD_ARG;
+    goto error;
+  }
+  if(!nvoxels[0] || !nvoxels[1] || !nvoxels[2]) {
+    log_err(dev,
+      "%s: the number of voxels along each axis cannot be null.\n"
+      "\t#voxels XYZ = {%lu, %lu, %lu}.\n",
+      FUNC_NAME,
+      (unsigned long)nvoxels[0],
+      (unsigned long)nvoxels[1],
+      (unsigned long)nvoxels[2]);
+    res = RES_BAD_ARG;
+    goto error;
+  }
+  oct = MEM_CALLOC(dev->allocator, 1, sizeof(struct svx_octree));
+  if(!oct) {
+    res = RES_MEM_ERR;
+    goto error;
+  }
+  ref_init(&oct->ref);
+  SVX(device_ref_get(dev));
+  octree_buffer_init(dev->allocator, desc->size, &oct->buffer);
+  oct->dev = dev;
+
+  /* Compute the octree definition */
+  oct->definition = MMAX(nvoxels[0], 2);
+  oct->definition = MMAX(nvoxels[1], oct->definition);
+  oct->definition = MMAX(nvoxels[2], oct->definition);
+  oct->definition = round_up_pow2(oct->definition);
+
+  /* Intialize the octree builder */
+  res = octree_builder_init(&bldr, oct->definition, desc, &oct->buffer);
+  if(res != RES_OK) goto error;
+
+  vox.data = MEM_CALLOC(dev->allocator, 1, desc->size);
+  if(!vox.data) {
+    res = RES_MEM_ERR;
+    goto error;
+  }
+
+  mcode_max = oct->definition * oct->definition * oct->definition;
+  FOR_EACH(mcode, 0, mcode_max) {
+    size_t xyz[3];
+    uint32_t ui3[3];
+
+    morton_xyz_decode_u21(mcode, ui3);
+
+    /* Out of bound voxels */
+    if(ui3[0] >= nvoxels[0]
+    || ui3[1] >= nvoxels[1]
+    || ui3[2] >= nvoxels[2])
+      continue;
+
+    /* Retrieve the voxel data from the caller */
+    xyz[0] = (size_t)ui3[0];
+    xyz[1] = (size_t)ui3[1];
+    xyz[2] = (size_t)ui3[2];
+    desc->get(xyz, vox.data, desc->context);
+    vox.mcode = mcode;
+
+    /* Register the voxel against the octree */
+    res = octree_builder_add_voxel(&bldr, &vox);
+    if(res != RES_OK) goto error;
+  }
+
+  res = octree_builder_finalize(&bldr, &oct->root, oct->root_attr);
+  if(res != RES_OK) goto error;
+
+  oct->nleaves = bldr.nleaves;
+  oct->depth = (size_t)(bldr.octree_depth - bldr.non_empty_lvl)
+    + 1 /* leaf level */;
+  ASSERT(bldr.octree_depth > bldr.non_empty_lvl);
+
+#ifndef NDEBUG
+  check_octree(&oct->buffer, oct->root);
+#endif
+  /* Setup the octree AABB in world space */
+  d3_set(oct->lower, lower);
+  d3_set(oct->upper, upper);
+
+  /* Compute the voxel size in world space */
+  vox_sz[0] = (upper[0] - lower[0]) / (double)nvoxels[0];
+  vox_sz[1] = (upper[1] - lower[1]) / (double)nvoxels[1];
+  vox_sz[2] = (upper[2] - lower[2]) / (double)nvoxels[2];
+
+  /* Compute the octree AABB in world space */
+  d3_set(oct->oclow, lower);
+  oct->ocupp[0] = (double)oct->definition * vox_sz[0];
+  oct->ocupp[1] = (double)oct->definition * vox_sz[1];
+  oct->ocupp[2] = (double)oct->definition * vox_sz[2];
+
+  /* Adjust the octree definition with respect to the octree depth since finest
+   * levels might be removed by the build process due to the merging process */
+  oct->definition = oct->definition / (1u << bldr.non_empty_lvl);
+
+exit:
+  if(vox.data) MEM_RM(dev->allocator, vox.data);
+  if(out_oct) *out_oct = oct;
+  return res;
+error:
+  if(oct) {
+    SVX(octree_ref_put(oct));
+    oct = NULL;
+  }
+  goto exit;
+}
+
+res_T
+svx_octree_ref_get(struct svx_octree* oct)
+{
+  if(!oct) return RES_BAD_ARG;
+  ref_get(&oct->ref);
+  return RES_OK;
+}
+
+res_T
+svx_octree_ref_put(struct svx_octree* oct)
+{
+  if(!oct) return RES_BAD_ARG;
+  ref_put(&oct->ref, octree_release);
+  return RES_OK;
+}
+
+res_T
+svx_octree_get_desc
+  (const struct svx_octree* oct, struct svx_octree_desc* desc)
+{
+  if(!oct || !desc) return RES_BAD_ARG;
+  d3_set(desc->lower, oct->lower);
+  d3_set(desc->upper, oct->upper);
+  desc->nleaves = oct->nleaves;
+  desc->nvoxels =  absolute_attr_index(&oct->buffer, oct->buffer.attr_head)
+    + 1; /* Root node */
+  desc->depth = oct->depth;
+  return RES_OK;
+}
+
+res_T
+svx_octree_for_each_leaf
+  (struct svx_octree* oct,
+   void (*func)
+      (const void* val,
+       const size_t ileaf,
+       const double low[3],
+       const double upp[3],
+       void* ctx),
+   void* ctx)
+{
+  struct stack_entry {
+    struct octree_index inode;
+    double low[3];
+    double upp[3];
+  } stack[OCTREE_DEPTH_MAX*8];
+  int istack;
+  size_t ileaf = 0;
+
+  if(!oct || !func) return RES_BAD_ARG;
+
+  stack[0].inode = oct->root;
+  stack[0].low[0] = oct->oclow[0];
+  stack[0].low[1] = oct->oclow[1];
+  stack[0].low[2] = oct->oclow[2];
+  stack[0].upp[0] = oct->ocupp[0];
+  stack[0].upp[1] = oct->ocupp[1];
+  stack[0].upp[2] = oct->ocupp[2];
+  istack = 1;
+
+  do {
+    const struct stack_entry entry = stack[--istack];
+    struct octree_xnode* node;
+    int ichild;
+
+    node = octree_buffer_get_node(&oct->buffer, entry.inode);
+
+    FOR_EACH(ichild, 0, 8) {
+      const uint8_t ichild_flag = (uint8_t)BIT(ichild);
+      double half_sz[3]; /* Half size of the current node */
+      double mid[3]; /* Middle point of the current node */
+      double upp[3];
+      double low[3];
+
+      if((node->is_valid & ichild_flag) == 0) continue; /* Empty node */
+
+      half_sz[0] = (entry.upp[0] - entry.low[0])*0.5;
+      half_sz[1] = (entry.upp[1] - entry.low[1])*0.5;
+      half_sz[2] = (entry.upp[2] - entry.low[2])*0.5;
+      mid[0] = entry.low[0] + half_sz[0];
+      mid[1] = entry.low[1] + half_sz[1];
+      mid[2] = entry.low[2] + half_sz[2];
+      low[0] = ichild&4 ? mid[0] : entry.low[0];
+      low[1] = ichild&2 ? mid[1] : entry.low[1];
+      low[2] = ichild&1 ? mid[2] : entry.low[2];
+      upp[0] = low[0] + half_sz[0];
+      upp[1] = low[1] + half_sz[1];
+      upp[2] = low[2] + half_sz[2];
+
+      if(node->is_leaf & ichild_flag) {
+        struct octree_index iattr;
+        const void* val;
+
+        iattr = octree_buffer_get_child_attr_index
+          (&oct->buffer, entry.inode, ichild);
+        val = octree_buffer_get_attr(&oct->buffer, iattr);
+
+        ASSERT(upp[0] <= oct->upper[0]);
+        ASSERT(upp[1] <= oct->upper[1]);
+        ASSERT(upp[2] <= oct->upper[2]);
+        ASSERT(low[0] >= oct->lower[0]);
+        ASSERT(low[1] >= oct->lower[1]);
+        ASSERT(low[2] >= oct->lower[2]);
+
+        func(val, ileaf, low, upp, ctx);
+        ileaf++;
+      } else {
+        struct stack_entry* top = stack + istack;
+
+        top->inode = octree_buffer_get_child_node_index
+          (&oct->buffer, entry.inode, ichild);
+        top->low[0] = low[0];
+        top->low[1] = low[1];
+        top->low[2] = low[2];
+        top->upp[0] = upp[0];
+        top->upp[1] = upp[1];
+        top->upp[2] = upp[2];
+        ++istack;
+      }
+    }
+  } while(istack);
+
+  return RES_OK;
+}
+
+res_T
+svx_octree_at
+  (struct svx_octree* oct,
+   const double position[3],
+   int (*filter)  /* Filter function. May be NULL */
+    (const struct svx_voxel* voxel,
+     const double position[3],
+     void* ctx),
+   void* context, /* Client data sent as the last argument of the filter func */
+   struct svx_voxel* voxel)
+{
+  struct octree_index inode;
+  struct octree_index iattr;
+  struct svx_voxel vox = SVX_VOXEL_NULL;
+  double scale_exp2;
+  double low[3];
+  double pos[3];
+  res_T res = RES_OK;
+
+  if(!oct || !position || !voxel) {
+    res = RES_BAD_ARG;
+    goto error;
+  }
+
+  *voxel = SVX_VOXEL_NULL;
+
+  /* The position is outside the octree */
+  if(position[0] > oct->upper[0] || position[0] < oct->lower[0]
+  || position[1] > oct->upper[1] || position[1] < oct->lower[1]
+  || position[2] > oct->upper[2] || position[2] < oct->lower[2]) {
+    goto exit;
+  }
+
+  /* Transform the position in the normalized octree space,
+   * i.e. octree lies in [0, 1]^2 */
+  pos[0] = (position[0] - oct->oclow[0]) / (oct->ocupp[0] - oct->oclow[0]);
+  pos[1] = (position[1] - oct->oclow[1]) / (oct->ocupp[1] - oct->oclow[1]);
+  pos[2] = (position[2] - oct->oclow[2]) / (oct->ocupp[2] - oct->oclow[2]);
+
+  /* Initialized the lower left corner of the current node */
+  low[0] = 0;
+  low[1] = 0;
+  low[2] = 0;
+
+  /* Root voxel */
+  vox.depth = 0;
+  vox.is_leaf = 0;
+  if(filter) {
+    vox.data = oct->root_attr;
+    vox.id = absolute_attr_index(&oct->buffer, oct->buffer.attr_head);
+    if(!filter(&vox, position, context)) { *voxel = vox; goto exit; }
+  }
+
+  scale_exp2 = 0.5;
+  inode = oct->root;
+  for(;;) {
+    struct octree_xnode* node = octree_buffer_get_node(&oct->buffer, inode);
+    int ichild;
+    uint8_t ichild_flag;
+    double mid[3];
+
+    ++vox.depth;
+
+    /* Compute the middle point of the node */
+    mid[0] = low[0] + scale_exp2;
+    mid[1] = low[1] + scale_exp2;
+    mid[2] = low[2] + scale_exp2;
+
+    /* Define the child of the current node into which pos `lies' and compute
+     * its lower left corner */
+    ichild = 0;
+    if(pos[0] > mid[0]) { ichild |= 4; low[0] = mid[0]; }
+    if(pos[1] > mid[1]) { ichild |= 2; low[1] = mid[1]; }
+    if(pos[2] > mid[2]) { ichild |= 1; low[2] = mid[2]; }
+
+    ichild_flag = (uint8_t)BIT(ichild);
+    if((node->is_valid & ichild_flag) == 0) break; /* Empty node */
+
+    vox.is_leaf = (node->is_leaf & ichild_flag) != 0;
+    if(filter || vox.is_leaf) {
+      iattr = octree_buffer_get_child_attr_index(&oct->buffer, inode, ichild);
+      vox.data = octree_buffer_get_attr(&oct->buffer, iattr);
+      vox.id = absolute_attr_index(&oct->buffer, iattr);
+      vox.is_leaf = (node->is_leaf & ichild_flag) != 0;
+      if(vox.is_leaf || !filter(&vox, position, context)) {
+        *voxel = vox;
+        break;
+      }
+    }
+
+    inode = octree_buffer_get_child_node_index(&oct->buffer, inode, ichild);
+    scale_exp2 *= 0.5;
+  }
+
+exit:
+  return res;
+error:
+  goto exit;
+}
+
diff --git a/src/svx_octree.h b/src/svx_octree.h
@@ -0,0 +1,40 @@
+/* Copyright (C) 2018 CNRS
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. */
+
+#ifndef SVX_OCTREE_H
+#define SVX_OCTREE_H
+
+#include "svx_octree_buffer.h"
+#include <rsys/ref_count.h>
+
+struct svx_octree {
+  size_t definition; /* #voxels along the X, Y and Z axis */
+
+  double lower[3], upper[3]; /* Submitted World space AABB of the octree  */
+  double oclow[3], ocupp[3]; /* Adjusted World space AABB of the octree */
+
+  struct octree_buffer buffer;
+  struct octree_index root; /* Index toward the children of the root */
+  ALIGN(16) char root_attr[SVX_MAX_SIZEOF_VOXEL]; /* Attribute of the root */
+
+  size_t nleaves; /* #leaves */
+  size_t depth; /* Maximum depth of the octree */
+
+  struct svx_device* dev;
+  ref_T ref;
+};
+
+#endif /* SVX_OCTREE_H */
+
diff --git a/src/svx_octree_buffer.c b/src/svx_octree_buffer.c
@@ -0,0 +1,198 @@
+/* Copyright (C) 2018 CNRS
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. */
+
+#include "svx_octree_buffer.h"
+
+#include <rsys/mem_allocator.h>
+
+#include <unistd.h>
+
+/*******************************************************************************
+ * Helper functions
+ ******************************************************************************/
+static INLINE res_T
+ensure_allocated_nodes(struct octree_buffer* buf, const size_t nnodes)
+{
+  char* node_page = NULL;
+  size_t nnode_pages = 0;
+  res_T res = RES_OK;
+  ASSERT(buf);
+
+  if(buf->node_head.ipage != OCTREE_INDEX_NULL.ipage
+  && buf->node_head.inode + nnodes <= buf->pagesize/sizeof(struct octree_xnode))
+    goto exit;
+
+  nnode_pages = darray_page_size_get(&buf->node_pages);
+  if(nnode_pages > UINT32_MAX) { res = RES_MEM_ERR; goto error; }
+  ASSERT(nnode_pages == buf->node_head.ipage + 1);
+
+  /* Alloc and register a node page containing the node and the far indices */
+  node_page = MEM_ALLOC(buf->allocator, buf->pagesize);
+  if(!node_page) { res = RES_MEM_ERR; goto error; }
+  res = darray_page_push_back(&buf->node_pages, &node_page);
+  if(res != RES_OK) goto error;
+
+  buf->node_head.inode = 0;
+  buf->node_head.ipage = (uint32_t)nnode_pages;
+
+exit:
+  return res;
+error:
+  if(node_page) MEM_RM(buf->allocator, node_page);
+  CHK(darray_page_resize(&buf->node_pages, nnode_pages) == RES_OK);
+  goto exit;
+}
+
+static INLINE res_T
+ensure_allocated_attrs(struct octree_buffer* buf, const size_t nattrs)
+{
+  char* attr_page = NULL;
+  size_t nattr_pages = 0;
+  res_T res = RES_OK;
+  ASSERT(buf);
+
+  if(buf->attr_head.ipage != OCTREE_INDEX_NULL.ipage
+  && buf->attr_head.inode + nattrs <= buf->pagesize/buf->voxsize)
+    goto exit;
+
+  nattr_pages = darray_page_size_get(&buf->attr_pages);
+  if(nattr_pages > UINT32_MAX) { res = RES_MEM_ERR; goto error; }
+  ASSERT(nattr_pages == buf->attr_head.ipage + 1);
+
+  /* Alloc and register a attr page */
+  attr_page = MEM_ALLOC(buf->allocator, buf->pagesize);
+  if(!attr_page) { res = RES_MEM_ERR; goto error; }
+  res = darray_page_push_back(&buf->attr_pages, &attr_page);
+  if(res != RES_OK) goto error;
+
+  buf->attr_head.inode = 0;
+  buf->attr_head.ipage = (uint32_t)nattr_pages;
+
+exit:
+  return res;
+error:
+  if(attr_page) MEM_RM(buf->allocator, attr_page);
+  CHK(darray_page_resize(&buf->attr_pages, nattr_pages) == RES_OK);
+  goto exit;
+}
+
+/*******************************************************************************
+ * Local functions
+ ******************************************************************************/
+void
+octree_buffer_init
+  (struct mem_allocator* allocator,
+   const size_t voxel_size,
+   struct octree_buffer* buf)
+{
+  ASSERT(buf && allocator);
+  memset(buf, 0, sizeof(struct octree_buffer));
+  buf->pagesize = (size_t)sysconf(_SC_PAGESIZE);
+  buf->voxsize = voxel_size;
+  darray_page_init(allocator, &buf->node_pages);
+  darray_page_init(allocator, &buf->attr_pages);
+  buf->node_head = OCTREE_INDEX_NULL;
+  buf->attr_head = OCTREE_INDEX_NULL;
+  buf->allocator = allocator;
+  CHK(buf->voxsize <= buf->pagesize);
+}
+
+void
+octree_buffer_release(struct octree_buffer* buf)
+{
+  ASSERT(buf);
+  octree_buffer_clear(buf);
+  darray_page_release(&buf->node_pages);
+  darray_page_release(&buf->attr_pages);
+}
+
+res_T
+octree_buffer_alloc_nodes
+  (struct octree_buffer* buf,
+   const size_t nnodes,
+   struct octree_index* first_node)
+{
+  res_T res = RES_OK;
+  ASSERT(buf && first_node);
+
+  if(nnodes > buf->pagesize / sizeof(struct octree_xnode))
+    return RES_MEM_ERR;
+
+  res = ensure_allocated_nodes(buf, nnodes);
+  if(res != RES_OK) return res;
+
+  *first_node = buf->node_head;
+  buf->node_head.inode = (uint16_t)(buf->node_head.inode + nnodes);
+  return RES_OK;
+
+}
+
+res_T
+octree_buffer_alloc_attrs
+  (struct octree_buffer* buf,
+   const size_t nattrs,
+   struct octree_index* first_attr)
+{
+  res_T res = RES_OK;
+  ASSERT(buf && first_attr);
+
+  if(nattrs > buf->pagesize / buf->voxsize) return RES_MEM_ERR;
+
+  res = ensure_allocated_attrs(buf, nattrs);
+  if(res != RES_OK) return res;
+
+  *first_attr = buf->attr_head;
+  buf->attr_head.inode = (uint16_t)(buf->attr_head.inode + nattrs);
+  return RES_OK;
+}
+
+res_T
+octree_buffer_alloc_far_index
+  (struct octree_buffer* buf,
+   struct octree_index* id)
+{
+  size_t remaining_size;
+  size_t skipped_nnodes;
+  STATIC_ASSERT(sizeof(struct octree_index) >= sizeof(struct octree_xnode),
+    Unexpected_type_size);
+
+  remaining_size = buf->pagesize - buf->node_head.inode*sizeof(struct octree_xnode);
+
+  /* Not enough memory in the current page */
+  if(sizeof(struct octree_index) > remaining_size) return RES_MEM_ERR;
+
+  *id = buf->node_head;
+  skipped_nnodes = sizeof(struct octree_index) / sizeof(struct octree_xnode);
+  buf->node_head.inode = (uint16_t)(buf->node_head.inode + skipped_nnodes);
+  return RES_OK;
+}
+
+void
+octree_buffer_clear(struct octree_buffer* buf)
+{
+  size_t i;
+  ASSERT(buf);
+  FOR_EACH(i, 0, darray_page_size_get(&buf->node_pages)) {
+    MEM_RM(buf->allocator, darray_page_data_get(&buf->node_pages)[i]);
+  }
+  FOR_EACH(i, 0, darray_page_size_get(&buf->attr_pages)) {
+    MEM_RM(buf->allocator, darray_page_data_get(&buf->attr_pages)[i]);
+  }
+  darray_page_purge(&buf->node_pages);
+  darray_page_purge(&buf->attr_pages);
+  buf->node_head = OCTREE_INDEX_NULL;
+  buf->attr_head = OCTREE_INDEX_NULL;
+}
+
diff --git a/src/svx_octree_buffer.h b/src/svx_octree_buffer.h
@@ -0,0 +1,239 @@
+/* Copyright (C) 2018 CNRS
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. */
+
+#ifndef SVX_OCTREE_BUFFER_H
+#define SVX_OCTREE_BUFFER_H
+
+#include "svx_c.h"
+#include <rsys/dynamic_array.h>
+
+/*
+ * Buffer containing the data of an octree. These data are partitioned in fixed
+ * size memory pages whose capacity is defined on buffer initialisation with
+ * respect to the page size of the system.
+ *
+ * The children of a node are stored consecutively into a page. The parent node
+ * directly references its first valid children excepted if they lie in a
+ * different page. In this case, the node references a `struct octree_index',
+ * stored into the same page, that defines the absolute position of its first
+ * valid child into the whole list of node pages
+ *
+ * The data of the nodes are stored in their own memory pages. The attribs of
+ * the children of a node are stored consecutively into a page. If the page
+ * identifier of the attribs is the same of the page into which their parent
+ * node lies, then the node saves the index toward the first valid attrib into
+ * the page of attribs. In the other case, the node references a `struct
+ * octree_index', stored into the same page of the node, that defines the
+ * absolute position of its first valid attrib into the buffer of attribs.
+ */
+
+#define OCTREE_XNODE_FLAG_FAR_INDEX (1u<<15)
+#define OCTREE_XNODE_MAX_CHILDREN_OFFSET (OCTREE_XNODE_FLAG_FAR_INDEX-1)
+#define OCTREE_XNODE_MASK OCTREE_XNODE_MAX_CHILDREN_OFFSET
+
+struct octree_xnode {
+  /* Offset to retrieve the children. If the OCTREE_XNODE_FLAG_FAR_INDEX bit is
+   * not set, the children are stored in the same page at the position `offset
+   * & OCTREE_XNODE_MASK'. If OCTREE_XNODE_FLAG_FAR_INDEX is set, `offset &
+   * OCTREE_XNODE_MASK' reference an octree_index toward the node children */
+  uint16_t node_offset;
+  uint16_t attr_offset;
+  uint8_t is_valid; /* Mask defining if the children are valid */
+  uint8_t is_leaf; /* Mask defining if the children are leaves */
+  uint16_t dummy__; /* Ensure a size of 8 Bytes */
+};
+STATIC_ASSERT(sizeof(struct octree_xnode) == 8,
+  Unexpected_sizeof_octree_xnode);
+
+#define OCTREE_INDEX_IPAGE_MAX UINT32_MAX
+#define OCTREE_INDEX_INODE_MAX UINT16_MAX
+
+struct octree_index {
+  uint32_t ipage; /* Identifier of the page */
+  uint16_t inode; /* Identifier of the node in the page */
+  uint16_t dummy__; /* Padding to ensure the octree index is 8 bytes lenght */
+};
+STATIC_ASSERT(sizeof(struct octree_index) == 8,
+  Unexpected_sizeof_octree_index);
+#define OCTREE_INDEX_NULL__ {UINT32_MAX, UINT16_MAX, UINT16_MAX}
+static const struct octree_index OCTREE_INDEX_NULL = OCTREE_INDEX_NULL__;
+#define OCTREE_INDEX_EQ(A, B) ((A)->inode==(B)->inode && (A)->ipage==(B)->ipage)
+
+/* Define the dynamic array of pages */
+#define DARRAY_NAME page
+#define DARRAY_DATA char*
+#include <rsys/dynamic_array.h>
+
+struct octree_buffer {
+  size_t pagesize; /* Memory page size in bytes */
+  size_t voxsize; /* Memory size of a voxel in bytes */
+
+  struct darray_page node_pages; /* List of pages storing nodes */
+  struct darray_page attr_pages; /* List of pages storing node attributes */
+  struct octree_index node_head; /* Index of the next valid node */
+  struct octree_index attr_head; /* Index of the next valid attr */
+
+  struct mem_allocator* allocator;
+};
+
+extern LOCAL_SYM void
+octree_buffer_init
+  (struct mem_allocator* allocator,
+   const size_t sizeof_voxel, /* Size in bytes of a voxel */
+   struct octree_buffer* buf);
+
+extern LOCAL_SYM void
+octree_buffer_release
+  (struct octree_buffer* buf);
+
+extern LOCAL_SYM res_T
+octree_buffer_alloc_nodes
+  (struct octree_buffer* buf,
+   const size_t nnodes,
+   struct octree_index* first_node); /* Index toward the 1st allocated node */
+
+extern LOCAL_SYM res_T
+octree_buffer_alloc_attrs
+  (struct octree_buffer* buf,
+   const size_t nattrs,
+   struct octree_index* first_attr); /* Index toward the 1st allocated attrib */
+
+/* Allocate an octree_index in the current buffer page. Return RES_MEM_ERR if
+ * the node index cannot be allocated in the current page. In this case one
+ * have to alloc new nodes */
+extern LOCAL_SYM res_T
+octree_buffer_alloc_far_index
+  (struct octree_buffer* buf,
+   struct octree_index* id); /* Index toward the allocated far index */
+
+extern LOCAL_SYM void
+octree_buffer_clear
+  (struct octree_buffer* buf);
+
+static FINLINE int
+octree_buffer_is_empty(const struct octree_buffer* buf)
+{
+  ASSERT(buf);
+  return darray_page_size_get(&buf->node_pages) == 0;
+}
+
+static FINLINE struct octree_xnode*
+octree_buffer_get_node
+  (struct octree_buffer* buf,
+   const struct octree_index id)
+{
+  char* mem;
+  ASSERT(buf && id.inode < buf->pagesize/sizeof(struct octree_xnode));
+  ASSERT(id.ipage < darray_page_size_get(&buf->node_pages));
+  mem = darray_page_data_get(&buf->node_pages)[id.ipage];
+  mem += id.inode * sizeof(struct octree_xnode);
+  return (struct octree_xnode*)mem;
+}
+
+static FINLINE void*
+octree_buffer_get_attr
+  (struct octree_buffer* buf,
+   const struct octree_index id)
+{
+  char* mem;
+  ASSERT(buf && id.inode < buf->pagesize/buf->voxsize);
+  ASSERT(id.ipage < darray_page_size_get(&buf->attr_pages));
+  mem = darray_page_data_get(&buf->attr_pages)[id.ipage];
+  mem += id.inode * buf->voxsize;
+  return mem;
+}
+
+static FINLINE struct octree_index*
+octree_buffer_get_far_index
+  (struct octree_buffer* buf,
+   const struct octree_index id)
+{
+  char* mem;
+  ASSERT(buf && id.inode < buf->pagesize/sizeof(struct octree_xnode));
+  ASSERT(id.ipage < darray_page_size_get(&buf->node_pages));
+  mem = darray_page_data_get(&buf->node_pages)[id.ipage];
+  mem += id.inode * sizeof(struct octree_xnode);
+  return (struct octree_index*)mem;
+}
+
+static FINLINE struct octree_index
+octree_buffer_get_child_node_index
+  (struct octree_buffer* buf,
+   const struct octree_index id,
+   const int ichild) /* in [0, 7] */
+{
+  struct octree_index child_id = OCTREE_INDEX_NULL;
+  struct octree_xnode* node = NULL;
+  uint16_t offset;
+  const int ichild_flag = BIT(ichild);
+  int ichild_off;
+  uint8_t mask;
+
+  ASSERT(ichild >= 0 && ichild < 8 && buf);
+
+  node = octree_buffer_get_node(buf, id);
+  mask = (uint8_t)(node->is_valid & ~node->is_leaf);
+  ASSERT(mask & ichild_flag);
+
+  /* Compute the child offset from the first child node */
+  ichild_off = popcount((uint8_t)((ichild_flag-1) & (int)mask));
+
+  offset = node->node_offset & OCTREE_XNODE_MASK;
+  if(!(node->node_offset & OCTREE_XNODE_FLAG_FAR_INDEX)) {
+    child_id.ipage = id.ipage;
+    child_id.inode = (uint16_t)(offset + ichild_off);
+  } else {
+    char* mem = darray_page_data_get(&buf->node_pages)[id.ipage];
+    child_id = *(struct octree_index*)(mem+offset*(sizeof(struct octree_xnode)));
+    child_id.inode = (uint16_t)(child_id.inode + ichild_off);
+  }
+  return child_id;
+}
+
+static FINLINE struct octree_index
+octree_buffer_get_child_attr_index
+  (struct octree_buffer* buf,
+   const struct octree_index id,
+   const int ichild) /* In [0, 7] */
+{
+  struct octree_index child_id = OCTREE_INDEX_NULL;
+  struct octree_xnode* node = NULL;
+  uint16_t offset;
+  const int ichild_flag = BIT(ichild);
+  int ichild_off;
+  uint8_t mask;
+
+  ASSERT(ichild >= 0 && ichild < 8 && buf);
+
+  node = octree_buffer_get_node(buf, id);
+  mask = node->is_valid;
+  ASSERT(mask & ichild_flag);
+
+  /* Compute the attr offset from the first child node */
+  ichild_off = popcount((uint8_t)((ichild_flag-1) & (int)mask));
+
+  offset = node->attr_offset & OCTREE_XNODE_MASK;
+  if(!(node->attr_offset & OCTREE_XNODE_FLAG_FAR_INDEX)) {
+    child_id.ipage = id.ipage;
+    child_id.inode = (uint16_t)(offset + ichild_off);
+  } else {
+    char* mem = darray_page_data_get(&buf->node_pages)[id.ipage];
+    child_id = *(struct octree_index*)(mem+offset*(sizeof(struct octree_xnode)));
+    child_id.inode = (uint16_t)(child_id.inode + ichild_off);
+  }
+  return child_id;
+}
+
+#endif /* SVX_OCTREE_BUFFER_H */
diff --git a/src/svx_scene_trace_ray.c b/src/svx_scene_trace_ray.c
@@ -0,0 +1,228 @@
+/* Copyright (C) 2018 CNRS
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. */
+
+#include "svx.h"
+#include "svx_scene.h"
+
+struct ray {
+  float org[3];
+  float range[2];
+  float ts[3]; /* 1 / -abs(dir) */
+  uint32_t octant_mask;
+};
+
+/*******************************************************************************
+ * Helper functions
+ ******************************************************************************/
+/* Return RES_BAD_OP if the ray cannot intersect the scene */
+static res_T
+setup_ray
+  (const struct scene* scn,
+   const double org[3],
+   const double dir[3],
+   const double range[2]
+   struct ray* ray)
+{
+  double rcp_ocsz[3]; /* Reciprocal size of the World space octree AABB */
+  ASSERT(scn);
+
+  if(range[0] >= range[1]) return RES_BAD_OP; /* Disabled ray */
+
+  /* Ray paralelle to an axis and that does not intersect the scene AABB */
+  if((!dir[0] && (org[0] < scn->lower[0] || org[0] > scn->upper[0]))
+  || (!dir[1] && (org[1] < scn->lower[1] || org[1] > scn->upper[1]))
+  || (!dir[2] && (org[2] < scn->lower[2] || org[2] > scn->upper[2]))) {
+    return RES_BAD_OP;
+  }
+
+  /* Compute reciprocal size of the world space octree AABB */
+  rcp_ocsz[0] = 1.0 / (scn->ocupp[0] - scn->oclow[0]);
+  rcp_ocsz[1] = 1.0 / (scn->ocupp[1] - scn->oclow[1]);
+  rcp_ocsz[2] = 1.0 / (scn->ocupp[2] - scn->oclow[2]);
+
+  /* Transform the ray origin in the [1, 2] space */
+  ray->org[0] = (float)((org[0] - scn->oclow[0]) * rcp_ocsz[0] + 1.f);
+  ray->org[1] = (float)((org[1] - scn->oclow[1]) * rcp_ocsz[1] + 1.f);
+  ray->org[2] = (float)((org[2] - scn->oclow[2]) * rcp_ocsz[2] + 1.f);
+
+  /* Transform the range in the normalized octree space */
+  ray->range[0] = (float)(range[0] * rcp_ocsz[0]);
+  ray->range[1] = (float)(range[1] * rcp_ocsz[1]);
+
+  /* Transform the direction in the normalized octree space */
+  ray->dir[0] = (float)(dir[0] * rcp_ocsz[0]);
+  ray->dir[1] = (float)(dir[1] * rcp_ocsz[1]);
+  ray->dir[2] = (float)(dir[2] * rcp_ocsz[2]);
+
+  /* Let a ray defines as org + t*dir and X the coordinate of an axis aligned
+   * plane. The ray intersects X at t = (X - org)/dir = (X - org) * ts; with ts
+   * = 1/dir. Note that one assume that dir is always negative. */
+  ray->ts[0] = (float)(1.0 / -fabs(dir[0]));
+  ray->ts[1] = (float)(1.0 / -fabs(dir[1]));
+  ray->ts[2] = (float)(1.0 / -fabs(dir[2]));
+
+  /* Mirror rays with position directions */
+  ray->octant_mask = 0;
+  if(ray_dir[0] > 0) { ray->octant_mask ^= 4; ray->org[0] = 3.f - ray->org[0]; }
+  if(ray_dir[1] > 0) { ray->octant_mask ^= 2; ray->org[0] = 3.f - ray->org[1]; }
+  if(ray_dir[2] > 0) { ray->octant_mask ^= 1; ray->org[0] = 3.f - ray->org[2]; }
+
+  return RES_OK;
+}
+
+static res_T
+trace_ray(const struct scene* scn, const struct* ray, struct svx_hit* hit)
+{
+  #define SCALE_MAX 23 /* #mantisse bits */
+  struct octree_index inode;
+  float scale_exp2;
+  float t_min, tmax;
+  float low[3], upp[3]; /* AABB of the current node in normalized space */
+  float mid[3]; /* Middle point of the current node in normalized space */
+  float corner[3];
+  uint32_t scale; /* stack index */
+  uint32_t ichild;
+  ASSERT(scn && ray && hit);
+  
+  hit = SVX_HIT_NONE; /* Initialise the hit to "no intersection" */
+
+  f3_splat(lower, 1);
+  f3_splat(upper, 2);
+
+  /* Compute the min/max ray intersection with the octree AABB in normalized
+   * space. Note that in this space, the octree AABB is in [1, 2]^3 */
+  t_min =      (2 - ray->org[0]) * ray->ts[0];
+  t_min = MMAX((2 - ray->org[1]) * ray->ts[1], t_min);
+  t_min = MMAX((2 - ray->org[2]) * ray->ts[2], t_min);  
+  t_max =      (1 - ray->org[0]) * ray->ts[0];
+  t_max = MMIN((1 - ray->org[1]) * ray->ts[1], t_max);
+  t_max = MMIN((1 - ray->org[2]) * ray->ts[2], t_max);
+  t_min = MMAX(ray->range[0], t_min);
+  t_max = MMIN(ray->range[1], t_max);
+  if(t_min >= t_max) return RES_OK; /* No intersection */
+
+  /* Traverrsal initialisation */
+  inode = scn->root;
+  scale_exp2 = 0.5f;
+  scale = SCALE_MAX - 1;
+
+  /* Define the first child id and the position of its lower left corner in the
+   * normalized octree space, i.e. in [1, 2]^3 */
+  ichild = 0;
+  corner[0] = 1.f;
+  corner[1] = 1.f;
+  corner[2] = 1.f;
+  if((1.5f - ray->org[0])*ray->ts[0] > t_min) { ichild ^= 4; corner[0] = 1.5f; }
+  if((1.5f - ray->org[1])*ray->ts[1] > t_min) { ichild ^= 2; corner[1] = 1.5f; }
+  if((1.5f - ray->org[2])*ray->ts[2] > t_min) { ichild ^= 1; corner[2] = 1.5f; }
+
+  /* Octree traversal */
+  scale_max = scale + 1;
+  while(scale < scale_max) {
+    const struct octree_xnode* node;
+    float t_corner[3];
+    float t_max_corner;
+    float t_max_child;
+    uint32_t ichild_adjusted = ichild ^ ray->octant_mask;
+    uint32_t istep;
+
+    inode = octree_buffer_get_child_index
+      (&scn->buffer, inode, (int)ichild_adjusted);
+    node = octree_buffer_get_node(&scn->buffer, inode);
+
+    /* Compute the exit point of the ray in the current child node */
+    t_corner[0] = (corner[0] - ray->org[0])*ray->ts[0];
+    t_corner[1] = (corner[1] - ray->org[1])*ray->ts[1];
+    t_corner[2] = (corner[2] - ray->org[2])*ray->ts[2];
+    t_max_corner = MMIN(MMIN(t_corner[0], t_corner[1]), t_corner[2]);
+
+    /* Traverse the current child */
+    if(!OCTREE_XNODE_IS_EMPTY(node)
+    && t_min <= (t_max_child = MMIN(t_max, t_max_corner))) {
+      float scale_exp2_child;
+      float t_max_parent;
+      float t_center[3];
+
+      t_max_parent = t_max;
+      t_max = t_max_child;
+
+      if(OCTREE_XNODE_IS_LEAF(node)) break; /* Leaf node */
+
+      scale_exp2_child = scale_exp2 * 0.5f;
+
+      /* center = corner - scale_exp2_child =>
+       * t_center = ts*(corner + scale_exp2_child - org)
+       * t_center = t_corner + ts*scale_exp2_child 
+       * Anyway we perforrm the whole computation to avoid numerical issues */
+      t_center[0] = (corner[0] - ray->org[0] + scale_exp2_child) * ray->ts[0];
+      t_center[1] = (corner[1] - ray->org[1] + scale_exp2_child) * ray->ts[1];
+      t_center[2] = (corner[2] - ray->org[2] + scale_exp2_child) * ray->ts[2];
+
+      /* Push the parent node */
+      stack[scale].t_max = t_max_parent;
+      stack[scale].inode = inode;
+
+      /* Define the id and the lower left corner of the first child */
+      ichild = 0;
+      if(t_center[0] > t_min) { ichild ^= 4; corner[0] += scale_exp2_child; }
+      if(t_center[1] > t_min) { ichild ^= 2; corner[1] += scale_exp2_child; }
+      if(t_center[2] > t_min) { ichild ^= 1; corner[2] += scale_exp2_child; }
+
+      --scale;
+      scale_exp2 = scale_exp2_child;
+      continue;
+    }
+
+    /* Define the id and the lower left corner of the next child */
+    istep = 0;
+    if(t_corner[0] <= t_max_corner) { istep ^= 4; corner[0] -= scale_exp2; }
+    if(t_corner[1] <= t_max_corner) { istep ^= 2; corner[1] -= scale_exp2; }
+    if(t_corner[2] <= t_max_corner) { istep ^= 1; corner[2] -= scale_exp2; }
+    ichild ^= istep; 
+  }
+}
+
+/*******************************************************************************
+ * Exported functions
+ ******************************************************************************/
+res_T
+svx_scene_trace_ray
+  (struct svx_scene* scn,
+   const double org[3],
+   const double dir[3],
+   const double range[2],
+   struct svx_hit* hit)
+{
+  struct ray ray;
+  res_T res = RES_OK;
+
+  if(!scn || !org || !range || !hit) {
+    res = RES_BAD_ARG;
+    goto error;
+  }
+
+  *hit = SVX_HIT_NULL;
+
+  res = setup_ray(scn, org, dir, range, &ray);
+  if(res == RES_BAD_OP) { /* The ray cannot intersect the scene. */
+    res = RES_OK;
+    goto exit;
+  }
+
+exit:
+  return res;
+error:
+  goto exit;
+}
diff --git a/src/test_htvox_device.c b/src/test_htvox_device.c
@@ -1,65 +0,0 @@
-/* Copyright (C) CNRS 2018
- *
- * This program is free software: you can redistribute it and/or modify
- * it under the terms of the GNU General Public License as published by
- * the Free Software Foundation, either version 3 of the License, or
- * (at your option) any later version.
- *
- * This program is distributed in the hope that it will be useful,
- * but WITHOUT ANY WARRANTY; without even the implied warranty of
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
- * GNU General Public License for more details.
- *
- * You should have received a copy of the GNU General Public License
- * along with this program. If not, see <http://www.gnu.org/licenses/>. */
-
-#include "htvox.h"
-#include "test_htvox_utils.h"
-
-#include <rsys/logger.h>
-
-static void
-log_stream(const char* msg, void* ctx)
-{
-  ASSERT(msg);
-  (void)msg, (void)ctx;
-  printf("%s\n", msg);
-}
-
-int
-main(int argc, char** argv)
-{
-  struct mem_allocator allocator;
-  struct logger logger;
-  struct htvox_device* htvox;
-  (void)argc, (void)argv;
-
-  CHK(htvox_device_create(NULL, NULL, 0, NULL) == RES_BAD_ARG);
-  CHK(htvox_device_create(NULL, NULL, 0, &htvox) == RES_OK);
-
-  CHK(htvox_device_ref_get(NULL) == RES_BAD_ARG);
-  CHK(htvox_device_ref_get(htvox) == RES_OK);
-  CHK(htvox_device_ref_put(NULL) == RES_BAD_ARG);
-  CHK(htvox_device_ref_put(htvox) == RES_OK);
-  CHK(htvox_device_ref_put(htvox) == RES_OK);
-
-  CHK(mem_init_proxy_allocator(&allocator, &mem_default_allocator) == RES_OK);
-  CHK(htvox_device_create(NULL, &allocator, 1, &htvox) == RES_OK);
-  CHK(htvox_device_ref_put(htvox) == RES_OK);
-
-  CHK(logger_init(&allocator, &logger) == RES_OK);
-  logger_set_stream(&logger, LOG_OUTPUT, log_stream, NULL);
-  logger_set_stream(&logger, LOG_ERROR, log_stream, NULL);
-  logger_set_stream(&logger, LOG_WARNING, log_stream, NULL);
-
-  CHK(htvox_device_create(&logger, &allocator, 0, &htvox) == RES_OK);
-  CHK(htvox_device_ref_put(htvox) == RES_OK);
-  CHK(htvox_device_create(&logger, NULL, 0, &htvox) == RES_OK);
-  CHK(htvox_device_ref_put(htvox) == RES_OK);
-
-  logger_release(&logger);
-  check_memory_allocator(&allocator);
-  mem_shutdown_proxy_allocator(&allocator);
-  CHK(mem_allocated_size() == 0);
-  return 0;
-}
diff --git a/src/test_htvox_octree.c b/src/test_htvox_octree.c
@@ -1,434 +0,0 @@
-/* Copyright (C) CNRS 2018
- *
- * This program is free software: you can redistribute it and/or modify
- * it under the terms of the GNU General Public License as published by
- * the Free Software Foundation, either version 3 of the License, or
- * (at your option) any later version.
- *
- * This program is distributed in the hope that it will be useful,
- * but WITHOUT ANY WARRANTY; without even the implied warranty of
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
- * GNU General Public License for more details.
- *
- * You should have received a copy of the GNU General Public License
- * along with this program. If not, see <http://www.gnu.org/licenses/>. */
-
-#include "htvox.h"
-#include "htvox_c.h" /* For morton_xyz_encode_u21 */
-#include "test_htvox_utils.h"
-
-#include <rsys/double3.h>
-
-struct check_context {
-  double* lower;
-  double* upper;
-  size_t* nvoxels;
-};
-
-static int
-no_merge(const void* voxels[], const size_t nvoxels, void* ctx)
-{
-  CHK(voxels != NULL);
-  CHK(nvoxels != 0);
-  CHK((intptr_t)ctx == 0xDECAFBAD);
-  return 0; /* Merge nothing */
-}
-
-static int
-merge_level0(const void** voxels, const size_t nvoxels, void* ctx)
-{
-  double min_val = DBL_MAX;
-  double max_val =-DBL_MAX;
-  size_t i;
-  CHK(voxels != NULL);
-  CHK(nvoxels != 0);
-  CHK((intptr_t)ctx == 0xDECAFBAD);
-
-  FOR_EACH(i, 0, nvoxels) {
-    const double* val = voxels[i];
-    min_val = MMIN(min_val, *val);
-    max_val = MMAX(max_val, *val);
-  }
-
-  return (max_val - min_val) < 8;
-}
-
-static void
-keep_max(void* dst, const void* voxels[], const size_t nvoxels,  void* ctx)
-{
-  double* vox_dst = dst;
-  double max_val = -DBL_MAX;
-  size_t i;
-
-  CHK(dst != NULL);
-  CHK(voxels != NULL);
-  CHK(nvoxels != 0);
-  CHK(ctx != NULL);
-
-  FOR_EACH(i, 0, nvoxels) {
-    const double* val = voxels[i];
-    max_val = MMAX(max_val, *val);
-  }
-  *vox_dst = max_val;
-}
-
-static void
-get(const size_t xyz[3], void* dst, void* ctx)
-{
-  uint32_t ui3[3];
-  uint64_t mcode;
-  double* val = dst;
-  CHK(xyz != NULL);
-  CHK(val != NULL);
-  CHK((intptr_t)ctx == 0xDECAFBAD);
-
-  ui3[0] = (uint32_t)xyz[0];
-  ui3[1] = (uint32_t)xyz[1];
-  ui3[2] = (uint32_t)xyz[2];
-
-  mcode = morton_xyz_encode_u21(ui3);
-  *val = (double)mcode;
-}
-
-static void
-check_voxel
-  (const void* val,
-   const size_t ivoxel,
-   const double low[3],
-   const double upp[3],
-   void* context)
-{
-  const double* dbl = val;
-  struct check_context* ctx = context;
-  uint64_t mcode;
-  uint32_t xyz[3];
-  double lower[3];
-  double delta[3];
-
-  CHK(val != NULL);
-  CHK(low != NULL);
-  CHK(upp != NULL);
-  CHK(ctx != NULL);
-  CHK(low[0] < upp[0]);
-  CHK(low[1] < upp[1]);
-  CHK(low[2] < upp[2]);
-  CHK(ivoxel < ctx->nvoxels[0]*ctx->nvoxels[1]*ctx->nvoxels[2]);
-  CHK(*dbl >= 0);
-
-  mcode = (uint64_t)(*dbl);
-  CHK(*dbl == (double)mcode);
-
-  delta[0] = (ctx->upper[0] - ctx->lower[0]) / (double)ctx->nvoxels[0];
-  delta[1] = (ctx->upper[1] - ctx->lower[1]) / (double)ctx->nvoxels[1];
-  delta[2] = (ctx->upper[2] - ctx->lower[2]) / (double)ctx->nvoxels[2];
-
-  morton_xyz_decode_u21(mcode, xyz);
-  lower[0] = xyz[0] * delta[0];
-  lower[1] = xyz[1] * delta[1];
-  lower[2] = xyz[2] * delta[2];
-
-  CHK(eq_eps(lower[0], low[0], 1.e-6));
-  CHK(eq_eps(lower[1], low[1], 1.e-6));
-  CHK(eq_eps(lower[2], low[2], 1.e-6));
-  CHK(eq_eps(lower[0] + delta[0], upp[0], 1.e-6));
-  CHK(eq_eps(lower[1] + delta[1], upp[1], 1.e-6));
-  CHK(eq_eps(lower[2] + delta[2], upp[2], 1.e-6));
-}
-
-static void
-write_points
-  (const void* val,
-   const size_t ivoxel,
-   const double low[3],
-   const double upp[3],
-   void* context)
-{
-  FILE* stream = context;
-  (void)val, (void)ivoxel;
-  CHK(stream != NULL);
-  CHK(low != NULL);
-  CHK(upp != NULL);
-  fprintf(stream,
-    "%g %g %g\n%g %g %g\n%g %g %g\n%g %g %g\n"
-    "%g %g %g\n%g %g %g\n%g %g %g\n%g %g %g\n",
-    low[0], low[1], low[2],
-    upp[0], low[1], low[2],
-    low[0], upp[1], low[2],
-    upp[0], upp[1], low[2],
-    low[0], low[1], upp[2],
-    upp[0], low[1], upp[2],
-    low[0], upp[1], upp[2],
-    upp[0], upp[1], upp[2]);
-}
-
-static void
-write_cells
-  (const void* val,
-   const size_t ivoxel,
-   const double low[3],
-   const double upp[3],
-   void* context)
-{
-  FILE* stream = context;
-  (void)ivoxel, (void)val;
-  CHK(stream != NULL);
-  CHK(low != NULL);
-  CHK(upp != NULL);
-  fprintf(stream, "8 %lu %lu %lu %lu %lu %lu %lu %lu\n",
-    (unsigned long)(ivoxel*8 + 0),
-    (unsigned long)(ivoxel*8 + 1),
-    (unsigned long)(ivoxel*8 + 2),
-    (unsigned long)(ivoxel*8 + 3),
-    (unsigned long)(ivoxel*8 + 4),
-    (unsigned long)(ivoxel*8 + 5),
-    (unsigned long)(ivoxel*8 + 6),
-    (unsigned long)(ivoxel*8 + 7));
-}
-
-static void
-write_scalars
-  (const void* val,
-   const size_t ivoxel,
-   const double low[3],
-   const double upp[3],
-   void* context)
-{
-  FILE* stream = context;
-  (void)ivoxel;
-  CHK(stream != NULL);
-  CHK(low != NULL);
-  CHK(upp != NULL);
-  fprintf(stream, "%g\n", *(double*)val);
-}
-
-static void
-dump_data(FILE* stream, struct htvox_octree* oct)
-{
-  struct htvox_octree_desc desc;
-  size_t ileaf;
-
-  CHK(stream != NULL);
-  CHK(oct != NULL);
-
-  fprintf(stream, "# vtk DataFile Version 2.0\n");
-  fprintf(stream, "Volume\n");
-  fprintf(stream, "ASCII\n");
-  fprintf(stream, "DATASET UNSTRUCTURED_GRID\n");
-
-  CHK(htvox_octree_get_desc(oct, &desc) == RES_OK);
-  fprintf(stream, "POINTS %lu float\n", (unsigned long)(desc.nleaves* 8));
-  CHK(htvox_octree_for_each_leaf(oct, write_points, stream) == RES_OK);
-
-  fprintf(stream, "CELLS %lu %lu\n",
-    (unsigned long)desc.nleaves,
-    (unsigned long)(desc.nleaves*(8/*#verts per leaf*/ + 1/*1st field of a cell*/)));
-  CHK(htvox_octree_for_each_leaf(oct, write_cells, stream) == RES_OK);
-
-  fprintf(stream, "CELL_TYPES %lu\n", (unsigned long)desc.nleaves);
-  FOR_EACH(ileaf, 0, desc.nleaves) fprintf(stream, "11\n");
-
-  fprintf(stream, "CELL_DATA %lu\n", (unsigned long)desc.nleaves);
-  fprintf(stream, "SCALARS K float 1\n");
-  fprintf(stream, "LOOKUP_TABLE default\n");
-  CHK(htvox_octree_for_each_leaf(oct, write_scalars, stream) == RES_OK);
-}
-
-static int
-max_lod
-  (const struct htvox_voxel* vox,
-   const double pos[3],
-   void* ctx)
-{
-  size_t* pdepth = ctx;
-
-  CHK(vox != NULL);
-  CHK(pos != NULL);
-  CHK(ctx != NULL);
-  CHK(vox->depth <= *pdepth);
-  return vox->depth < *pdepth;
-}
-
-
-static void
-test_at_accessor(struct htvox_octree* oct, const size_t nvoxels[3])
-{
-  struct htvox_octree_desc octdesc;
-  size_t nvxls;
-  double pos[3];
-  double delta[3];
-  size_t depth;
-  size_t x, y, z;
-
-  CHK(nvoxels != NULL);
-  CHK(htvox_octree_get_desc(oct, &octdesc) == RES_OK);
-
-  delta[0] = (octdesc.upper[0] - octdesc.lower[0])/(double)nvoxels[0];
-  delta[1] = (octdesc.upper[1] - octdesc.lower[1])/(double)nvoxels[1];
-  delta[2] = (octdesc.upper[2] - octdesc.lower[2])/(double)nvoxels[2];
-
-  nvxls = nvoxels[0];
-  nvxls = MMAX(nvxls, nvoxels[1]);
-  nvxls = MMAX(nvxls, nvoxels[2]);
-
-  depth = octdesc.depth;
-
-  while(depth-- != 0) {
-    const size_t iter = octdesc.depth - depth - 1;
-
-    FOR_EACH(x, 0, nvxls) {
-      pos[0] = octdesc.lower[0] + ((double)x+1.0/(1u<<(2+iter)))*delta[0];
-      if(x*(1u<<iter) >= nvoxels[0]) break;
-
-      FOR_EACH(y, 0, nvxls) {
-        pos[1] = octdesc.lower[1] + ((double)y+1.0/(1u<<(2+iter)))*delta[1];
-        if(y*(1u<<iter) >= nvoxels[1]) break;
-
-        FOR_EACH(z, 0, nvxls) {
-          struct htvox_voxel vox;
-          uint32_t ui3[3];
-          uint64_t mcode;
-
-          pos[2] = octdesc.lower[2] + ((double)z+1.0/(1u<<(2+iter)))*delta[2];
-          if(z*(1u<<iter) >= nvoxels[2]) break;
-
-          CHK(htvox_octree_at(oct, pos, max_lod, &depth, &vox) == RES_OK);
-          CHK(!HTVOX_VOXEL_NONE(&vox));
-
-          ui3[0] = (uint32_t)x * (1u << iter) + ((1u << iter) - 1);
-          ui3[1] = (uint32_t)y * (1u << iter) + ((1u << iter) - 1);
-          ui3[2] = (uint32_t)z * (1u << iter) + ((1u << iter) - 1);
-          ui3[0] = MMIN(ui3[0], (uint32_t)(nvoxels[0]-1));
-          ui3[1] = MMIN(ui3[1], (uint32_t)(nvoxels[1]-1));
-          ui3[2] = MMIN(ui3[2], (uint32_t)(nvoxels[2]-1));
-
-          mcode = morton_xyz_encode_u21(ui3);
-          CHK(*((double*)vox.data) == mcode);
-          CHK(vox.is_leaf == (octdesc.depth-1==depth));
-          CHK(vox.depth == depth);
-        }
-      }
-    }
-
-    nvxls = nvxls == 1 ? 0 : (nvxls + 1/*ceil*/) / 2;
-    d3_muld(delta, delta, 2);
-  }
-  CHK(nvxls == 0);
-}
-
-int
-main(int argc, char** argv)
-{
-  struct htvox_device* dev = NULL;
-  struct htvox_octree* oct = NULL;
-  struct mem_allocator allocator;
-  struct htvox_voxel_desc voxdesc = HTVOX_VOXEL_DESC_NULL;
-  struct htvox_octree_desc octdesc = HTVOX_OCTREE_DESC_NULL;
-  double low[3];
-  double upp[3];
-  size_t nvxls[3];
-  struct check_context ctx;
-  void* ptr = (void*)0xDECAFBAD;
-  (void)argc, (void)argv;
-
-  CHK(mem_init_proxy_allocator(&allocator, &mem_default_allocator) == RES_OK);
-
-  CHK(htvox_device_create(NULL, &allocator, 1, &dev) == RES_OK);
-
-  d3_splat(low, 0.0);
-  d3_splat(upp, 1.0);
-  nvxls[0] = nvxls[1] = nvxls[2] = 5;
-
-  ctx.lower = low;
-  ctx.upper = upp;
-  ctx.nvoxels = nvxls;
-
-  #define NEW_SCN htvox_octree_create
-
-  voxdesc.get = get;
-  voxdesc.merge = keep_max;
-  voxdesc.challenge_merge = no_merge;
-  voxdesc.context = ptr;
-  voxdesc.size = sizeof(double);
-  CHK(NEW_SCN(dev, low, upp, nvxls, &voxdesc, &oct) == RES_OK);
-
-  CHK(htvox_octree_ref_get(NULL) == RES_BAD_ARG);
-  CHK(htvox_octree_ref_get(oct) == RES_OK);
-  CHK(htvox_octree_ref_put(NULL) == RES_BAD_ARG);
-  CHK(htvox_octree_ref_put(oct) == RES_OK);
-  CHK(htvox_octree_ref_put(oct) == RES_OK);
-
-  upp[0] = low[0];
-  CHK(NEW_SCN(dev, low, upp, nvxls, &voxdesc, &oct) == RES_BAD_ARG);
-  upp[0] = 1.0;
-
-  nvxls[2] = 0;
-  CHK(NEW_SCN(dev, low, upp, nvxls, &voxdesc, &oct) == RES_BAD_ARG);
-  nvxls[2] = nvxls[0];
-
-  CHK(NEW_SCN(NULL, low, upp, nvxls, &voxdesc, &oct) == RES_BAD_ARG);
-  CHK(NEW_SCN(dev, NULL, upp, nvxls, &voxdesc, &oct) == RES_BAD_ARG);
-  CHK(NEW_SCN(dev, low, NULL, nvxls, &voxdesc, &oct) == RES_BAD_ARG);
-  CHK(NEW_SCN(dev, low, upp, NULL, &voxdesc, &oct) == RES_BAD_ARG);
-  CHK(NEW_SCN(dev, low, upp, nvxls, &voxdesc, NULL) == RES_BAD_ARG);
-
-  voxdesc.get = NULL;
-  CHK(NEW_SCN(dev, low, upp, nvxls, &voxdesc, &oct) == RES_BAD_ARG);
-  voxdesc.get = get;
-
-  voxdesc.merge = NULL;
-  CHK(NEW_SCN(dev, low, upp, nvxls, &voxdesc, &oct) == RES_BAD_ARG);
-  voxdesc.merge = keep_max;
-
-  voxdesc.challenge_merge = NULL;
-  CHK(NEW_SCN(dev, low, upp, nvxls, &voxdesc, &oct) == RES_BAD_ARG);
-  voxdesc.challenge_merge = no_merge;
-
-  voxdesc.size = 0;
-  CHK(NEW_SCN(dev, low, upp, nvxls, &voxdesc, &oct) == RES_BAD_ARG);
-  voxdesc.size = sizeof(double);
-
-  voxdesc.size = HTVOX_MAX_SIZEOF_VOXEL + 1;
-  CHK(NEW_SCN(dev, low, upp, nvxls, &voxdesc, &oct) == RES_BAD_ARG);
-  voxdesc.size = sizeof(double);
-
-  CHK(NEW_SCN(dev, low, upp, nvxls, &voxdesc, &oct) == RES_OK);
-
-  CHK(htvox_octree_for_each_leaf(oct, check_voxel, &ctx) == RES_OK);
-
-  CHK(htvox_octree_get_desc(NULL, &octdesc) == RES_BAD_ARG);
-  CHK(htvox_octree_get_desc(oct, NULL) == RES_BAD_ARG);
-  CHK(htvox_octree_get_desc(oct, &octdesc) == RES_OK);
-  CHK(octdesc.nleaves == 5*5*5);
-  CHK(octdesc.nvoxels == octdesc.nleaves + 36/*#parents*/);
-  CHK(octdesc.depth == 4);
-
-  d3_splat(low, 0);
-  d3_splat(upp, 1);
-  CHK(d3_eq(octdesc.lower, low));
-  CHK(d3_eq(octdesc.upper, upp));
-
-  test_at_accessor(oct, nvxls);
-
-  CHK(htvox_octree_ref_put(oct) == RES_OK);
-
-  nvxls[0] = nvxls[1] = nvxls[2] = 32;
-  voxdesc.challenge_merge = merge_level0;
-  CHK(NEW_SCN(dev, low, upp, nvxls, &voxdesc, &oct) == RES_OK);
-  CHK(htvox_octree_get_desc(oct, &octdesc) == RES_OK);
-  CHK(octdesc.nleaves == nvxls[0]*nvxls[1]*nvxls[2] / 8);
-  CHK(octdesc.nvoxels == (octdesc.nleaves*8 - 1) / 7);
-  CHK(octdesc.depth == (size_t)log2i((int)(nvxls[0]/2))+1);
-
-  dump_data(stdout, oct);
-
-  #undef NEW_SCN
-
-  CHK(htvox_device_ref_put(dev) == RES_OK);
-  CHK(htvox_octree_ref_put(oct) == RES_OK);
-
-  check_memory_allocator(&allocator);
-  mem_shutdown_proxy_allocator(&allocator);
-  CHK(mem_allocated_size() == 0);
-  return 0;
-}
-
diff --git a/src/test_htvox_utils.h b/src/test_htvox_utils.h
@@ -1,34 +0,0 @@
-/* Copyright (C) CNRS 2018
- *
- * This program is free software: you can redistribute it and/or modify
- * it under the terms of the GNU General Public License as published by
- * the Free Software Foundation, either version 3 of the License, or
- * (at your option) any later version.
- *
- * This program is distributed in the hope that it will be useful,
- * but WITHOUT ANY WARRANTY; without even the implied warranty of
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
- * GNU General Public License for more details.
- *
- * You should have received a copy of the GNU General Public License
- * along with this program. If not, see <http://www.gnu.org/licenses/>. */
-
-#ifndef TEST_HTVOX_UTILS_H
-#define TEST_HTVOX_UTILS_H
-
-#include <rsys/mem_allocator.h>
-#include <stdio.h>
-
-static INLINE void
-check_memory_allocator(struct mem_allocator* allocator)
-{
-  if(MEM_ALLOCATED_SIZE(allocator)) {
-    char dump[512];
-    MEM_DUMP(allocator, dump, sizeof(dump)/sizeof(char));
-    fprintf(stderr, "%s\n", dump);
-    FATAL("Memory leaks\n");
-  }
-}
-
-#endif /* TEST_HTVOX_UTILS_H */
-
diff --git a/src/test_svx_device.c b/src/test_svx_device.c
@@ -0,0 +1,65 @@
+/* Copyright (C) 2018 CNRS
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. */
+
+#include "svx.h"
+#include "test_svx_utils.h"
+
+#include <rsys/logger.h>
+
+static void
+log_stream(const char* msg, void* ctx)
+{
+  ASSERT(msg);
+  (void)msg, (void)ctx;
+  printf("%s\n", msg);
+}
+
+int
+main(int argc, char** argv)
+{
+  struct mem_allocator allocator;
+  struct logger logger;
+  struct svx_device* svx;
+  (void)argc, (void)argv;
+
+  CHK(svx_device_create(NULL, NULL, 0, NULL) == RES_BAD_ARG);
+  CHK(svx_device_create(NULL, NULL, 0, &svx) == RES_OK);
+
+  CHK(svx_device_ref_get(NULL) == RES_BAD_ARG);
+  CHK(svx_device_ref_get(svx) == RES_OK);
+  CHK(svx_device_ref_put(NULL) == RES_BAD_ARG);
+  CHK(svx_device_ref_put(svx) == RES_OK);
+  CHK(svx_device_ref_put(svx) == RES_OK);
+
+  CHK(mem_init_proxy_allocator(&allocator, &mem_default_allocator) == RES_OK);
+  CHK(svx_device_create(NULL, &allocator, 1, &svx) == RES_OK);
+  CHK(svx_device_ref_put(svx) == RES_OK);
+
+  CHK(logger_init(&allocator, &logger) == RES_OK);
+  logger_set_stream(&logger, LOG_OUTPUT, log_stream, NULL);
+  logger_set_stream(&logger, LOG_ERROR, log_stream, NULL);
+  logger_set_stream(&logger, LOG_WARNING, log_stream, NULL);
+
+  CHK(svx_device_create(&logger, &allocator, 0, &svx) == RES_OK);
+  CHK(svx_device_ref_put(svx) == RES_OK);
+  CHK(svx_device_create(&logger, NULL, 0, &svx) == RES_OK);
+  CHK(svx_device_ref_put(svx) == RES_OK);
+
+  logger_release(&logger);
+  check_memory_allocator(&allocator);
+  mem_shutdown_proxy_allocator(&allocator);
+  CHK(mem_allocated_size() == 0);
+  return 0;
+}
diff --git a/src/test_svx_octree.c b/src/test_svx_octree.c
@@ -0,0 +1,434 @@
+/* Copyright (C) 2018 CNRS
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. */
+
+#include "svx.h"
+#include "svx_c.h" /* For morton_xyz_encode_u21 */
+#include "test_svx_utils.h"
+
+#include <rsys/double3.h>
+
+struct check_context {
+  double* lower;
+  double* upper;
+  size_t* nvoxels;
+};
+
+static int
+no_merge(const void* voxels[], const size_t nvoxels, void* ctx)
+{
+  CHK(voxels != NULL);
+  CHK(nvoxels != 0);
+  CHK((intptr_t)ctx == 0xDECAFBAD);
+  return 0; /* Merge nothing */
+}
+
+static int
+merge_level0(const void** voxels, const size_t nvoxels, void* ctx)
+{
+  double min_val = DBL_MAX;
+  double max_val =-DBL_MAX;
+  size_t i;
+  CHK(voxels != NULL);
+  CHK(nvoxels != 0);
+  CHK((intptr_t)ctx == 0xDECAFBAD);
+
+  FOR_EACH(i, 0, nvoxels) {
+    const double* val = voxels[i];
+    min_val = MMIN(min_val, *val);
+    max_val = MMAX(max_val, *val);
+  }
+
+  return (max_val - min_val) < 8;
+}
+
+static void
+keep_max(void* dst, const void* voxels[], const size_t nvoxels,  void* ctx)
+{
+  double* vox_dst = dst;
+  double max_val = -DBL_MAX;
+  size_t i;
+
+  CHK(dst != NULL);
+  CHK(voxels != NULL);
+  CHK(nvoxels != 0);
+  CHK(ctx != NULL);
+
+  FOR_EACH(i, 0, nvoxels) {
+    const double* val = voxels[i];
+    max_val = MMAX(max_val, *val);
+  }
+  *vox_dst = max_val;
+}
+
+static void
+get(const size_t xyz[3], void* dst, void* ctx)
+{
+  uint32_t ui3[3];
+  uint64_t mcode;
+  double* val = dst;
+  CHK(xyz != NULL);
+  CHK(val != NULL);
+  CHK((intptr_t)ctx == 0xDECAFBAD);
+
+  ui3[0] = (uint32_t)xyz[0];
+  ui3[1] = (uint32_t)xyz[1];
+  ui3[2] = (uint32_t)xyz[2];
+
+  mcode = morton_xyz_encode_u21(ui3);
+  *val = (double)mcode;
+}
+
+static void
+check_voxel
+  (const void* val,
+   const size_t ivoxel,
+   const double low[3],
+   const double upp[3],
+   void* context)
+{
+  const double* dbl = val;
+  struct check_context* ctx = context;
+  uint64_t mcode;
+  uint32_t xyz[3];
+  double lower[3];
+  double delta[3];
+
+  CHK(val != NULL);
+  CHK(low != NULL);
+  CHK(upp != NULL);
+  CHK(ctx != NULL);
+  CHK(low[0] < upp[0]);
+  CHK(low[1] < upp[1]);
+  CHK(low[2] < upp[2]);
+  CHK(ivoxel < ctx->nvoxels[0]*ctx->nvoxels[1]*ctx->nvoxels[2]);
+  CHK(*dbl >= 0);
+
+  mcode = (uint64_t)(*dbl);
+  CHK(*dbl == (double)mcode);
+
+  delta[0] = (ctx->upper[0] - ctx->lower[0]) / (double)ctx->nvoxels[0];
+  delta[1] = (ctx->upper[1] - ctx->lower[1]) / (double)ctx->nvoxels[1];
+  delta[2] = (ctx->upper[2] - ctx->lower[2]) / (double)ctx->nvoxels[2];
+
+  morton_xyz_decode_u21(mcode, xyz);
+  lower[0] = xyz[0] * delta[0];
+  lower[1] = xyz[1] * delta[1];
+  lower[2] = xyz[2] * delta[2];
+
+  CHK(eq_eps(lower[0], low[0], 1.e-6));
+  CHK(eq_eps(lower[1], low[1], 1.e-6));
+  CHK(eq_eps(lower[2], low[2], 1.e-6));
+  CHK(eq_eps(lower[0] + delta[0], upp[0], 1.e-6));
+  CHK(eq_eps(lower[1] + delta[1], upp[1], 1.e-6));
+  CHK(eq_eps(lower[2] + delta[2], upp[2], 1.e-6));
+}
+
+static void
+write_points
+  (const void* val,
+   const size_t ivoxel,
+   const double low[3],
+   const double upp[3],
+   void* context)
+{
+  FILE* stream = context;
+  (void)val, (void)ivoxel;
+  CHK(stream != NULL);
+  CHK(low != NULL);
+  CHK(upp != NULL);
+  fprintf(stream,
+    "%g %g %g\n%g %g %g\n%g %g %g\n%g %g %g\n"
+    "%g %g %g\n%g %g %g\n%g %g %g\n%g %g %g\n",
+    low[0], low[1], low[2],
+    upp[0], low[1], low[2],
+    low[0], upp[1], low[2],
+    upp[0], upp[1], low[2],
+    low[0], low[1], upp[2],
+    upp[0], low[1], upp[2],
+    low[0], upp[1], upp[2],
+    upp[0], upp[1], upp[2]);
+}
+
+static void
+write_cells
+  (const void* val,
+   const size_t ivoxel,
+   const double low[3],
+   const double upp[3],
+   void* context)
+{
+  FILE* stream = context;
+  (void)ivoxel, (void)val;
+  CHK(stream != NULL);
+  CHK(low != NULL);
+  CHK(upp != NULL);
+  fprintf(stream, "8 %lu %lu %lu %lu %lu %lu %lu %lu\n",
+    (unsigned long)(ivoxel*8 + 0),
+    (unsigned long)(ivoxel*8 + 1),
+    (unsigned long)(ivoxel*8 + 2),
+    (unsigned long)(ivoxel*8 + 3),
+    (unsigned long)(ivoxel*8 + 4),
+    (unsigned long)(ivoxel*8 + 5),
+    (unsigned long)(ivoxel*8 + 6),
+    (unsigned long)(ivoxel*8 + 7));
+}
+
+static void
+write_scalars
+  (const void* val,
+   const size_t ivoxel,
+   const double low[3],
+   const double upp[3],
+   void* context)
+{
+  FILE* stream = context;
+  (void)ivoxel;
+  CHK(stream != NULL);
+  CHK(low != NULL);
+  CHK(upp != NULL);
+  fprintf(stream, "%g\n", *(double*)val);
+}
+
+static void
+dump_data(FILE* stream, struct svx_octree* oct)
+{
+  struct svx_octree_desc desc;
+  size_t ileaf;
+
+  CHK(stream != NULL);
+  CHK(oct != NULL);
+
+  fprintf(stream, "# vtk DataFile Version 2.0\n");
+  fprintf(stream, "Volume\n");
+  fprintf(stream, "ASCII\n");
+  fprintf(stream, "DATASET UNSTRUCTURED_GRID\n");
+
+  CHK(svx_octree_get_desc(oct, &desc) == RES_OK);
+  fprintf(stream, "POINTS %lu float\n", (unsigned long)(desc.nleaves* 8));
+  CHK(svx_octree_for_each_leaf(oct, write_points, stream) == RES_OK);
+
+  fprintf(stream, "CELLS %lu %lu\n",
+    (unsigned long)desc.nleaves,
+    (unsigned long)(desc.nleaves*(8/*#verts per leaf*/ + 1/*1st field of a cell*/)));
+  CHK(svx_octree_for_each_leaf(oct, write_cells, stream) == RES_OK);
+
+  fprintf(stream, "CELL_TYPES %lu\n", (unsigned long)desc.nleaves);
+  FOR_EACH(ileaf, 0, desc.nleaves) fprintf(stream, "11\n");
+
+  fprintf(stream, "CELL_DATA %lu\n", (unsigned long)desc.nleaves);
+  fprintf(stream, "SCALARS K float 1\n");
+  fprintf(stream, "LOOKUP_TABLE default\n");
+  CHK(svx_octree_for_each_leaf(oct, write_scalars, stream) == RES_OK);
+}
+
+static int
+max_lod
+  (const struct svx_voxel* vox,
+   const double pos[3],
+   void* ctx)
+{
+  size_t* pdepth = ctx;
+
+  CHK(vox != NULL);
+  CHK(pos != NULL);
+  CHK(ctx != NULL);
+  CHK(vox->depth <= *pdepth);
+  return vox->depth < *pdepth;
+}
+
+
+static void
+test_at_accessor(struct svx_octree* oct, const size_t nvoxels[3])
+{
+  struct svx_octree_desc octdesc;
+  size_t nvxls;
+  double pos[3];
+  double delta[3];
+  size_t depth;
+  size_t x, y, z;
+
+  CHK(nvoxels != NULL);
+  CHK(svx_octree_get_desc(oct, &octdesc) == RES_OK);
+
+  delta[0] = (octdesc.upper[0] - octdesc.lower[0])/(double)nvoxels[0];
+  delta[1] = (octdesc.upper[1] - octdesc.lower[1])/(double)nvoxels[1];
+  delta[2] = (octdesc.upper[2] - octdesc.lower[2])/(double)nvoxels[2];
+
+  nvxls = nvoxels[0];
+  nvxls = MMAX(nvxls, nvoxels[1]);
+  nvxls = MMAX(nvxls, nvoxels[2]);
+
+  depth = octdesc.depth;
+
+  while(depth-- != 0) {
+    const size_t iter = octdesc.depth - depth - 1;
+
+    FOR_EACH(x, 0, nvxls) {
+      pos[0] = octdesc.lower[0] + ((double)x+1.0/(1u<<(2+iter)))*delta[0];
+      if(x*(1u<<iter) >= nvoxels[0]) break;
+
+      FOR_EACH(y, 0, nvxls) {
+        pos[1] = octdesc.lower[1] + ((double)y+1.0/(1u<<(2+iter)))*delta[1];
+        if(y*(1u<<iter) >= nvoxels[1]) break;
+
+        FOR_EACH(z, 0, nvxls) {
+          struct svx_voxel vox;
+          uint32_t ui3[3];
+          uint64_t mcode;
+
+          pos[2] = octdesc.lower[2] + ((double)z+1.0/(1u<<(2+iter)))*delta[2];
+          if(z*(1u<<iter) >= nvoxels[2]) break;
+
+          CHK(svx_octree_at(oct, pos, max_lod, &depth, &vox) == RES_OK);
+          CHK(!SVX_VOXEL_NONE(&vox));
+
+          ui3[0] = (uint32_t)x * (1u << iter) + ((1u << iter) - 1);
+          ui3[1] = (uint32_t)y * (1u << iter) + ((1u << iter) - 1);
+          ui3[2] = (uint32_t)z * (1u << iter) + ((1u << iter) - 1);
+          ui3[0] = MMIN(ui3[0], (uint32_t)(nvoxels[0]-1));
+          ui3[1] = MMIN(ui3[1], (uint32_t)(nvoxels[1]-1));
+          ui3[2] = MMIN(ui3[2], (uint32_t)(nvoxels[2]-1));
+
+          mcode = morton_xyz_encode_u21(ui3);
+          CHK(*((double*)vox.data) == mcode);
+          CHK(vox.is_leaf == (octdesc.depth-1==depth));
+          CHK(vox.depth == depth);
+        }
+      }
+    }
+
+    nvxls = nvxls == 1 ? 0 : (nvxls + 1/*ceil*/) / 2;
+    d3_muld(delta, delta, 2);
+  }
+  CHK(nvxls == 0);
+}
+
+int
+main(int argc, char** argv)
+{
+  struct svx_device* dev = NULL;
+  struct svx_octree* oct = NULL;
+  struct mem_allocator allocator;
+  struct svx_voxel_desc voxdesc = SVX_VOXEL_DESC_NULL;
+  struct svx_octree_desc octdesc = SVX_OCTREE_DESC_NULL;
+  double low[3];
+  double upp[3];
+  size_t nvxls[3];
+  struct check_context ctx;
+  void* ptr = (void*)0xDECAFBAD;
+  (void)argc, (void)argv;
+
+  CHK(mem_init_proxy_allocator(&allocator, &mem_default_allocator) == RES_OK);
+
+  CHK(svx_device_create(NULL, &allocator, 1, &dev) == RES_OK);
+
+  d3_splat(low, 0.0);
+  d3_splat(upp, 1.0);
+  nvxls[0] = nvxls[1] = nvxls[2] = 5;
+
+  ctx.lower = low;
+  ctx.upper = upp;
+  ctx.nvoxels = nvxls;
+
+  #define NEW_SCN svx_octree_create
+
+  voxdesc.get = get;
+  voxdesc.merge = keep_max;
+  voxdesc.challenge_merge = no_merge;
+  voxdesc.context = ptr;
+  voxdesc.size = sizeof(double);
+  CHK(NEW_SCN(dev, low, upp, nvxls, &voxdesc, &oct) == RES_OK);
+
+  CHK(svx_octree_ref_get(NULL) == RES_BAD_ARG);
+  CHK(svx_octree_ref_get(oct) == RES_OK);
+  CHK(svx_octree_ref_put(NULL) == RES_BAD_ARG);
+  CHK(svx_octree_ref_put(oct) == RES_OK);
+  CHK(svx_octree_ref_put(oct) == RES_OK);
+
+  upp[0] = low[0];
+  CHK(NEW_SCN(dev, low, upp, nvxls, &voxdesc, &oct) == RES_BAD_ARG);
+  upp[0] = 1.0;
+
+  nvxls[2] = 0;
+  CHK(NEW_SCN(dev, low, upp, nvxls, &voxdesc, &oct) == RES_BAD_ARG);
+  nvxls[2] = nvxls[0];
+
+  CHK(NEW_SCN(NULL, low, upp, nvxls, &voxdesc, &oct) == RES_BAD_ARG);
+  CHK(NEW_SCN(dev, NULL, upp, nvxls, &voxdesc, &oct) == RES_BAD_ARG);
+  CHK(NEW_SCN(dev, low, NULL, nvxls, &voxdesc, &oct) == RES_BAD_ARG);
+  CHK(NEW_SCN(dev, low, upp, NULL, &voxdesc, &oct) == RES_BAD_ARG);
+  CHK(NEW_SCN(dev, low, upp, nvxls, &voxdesc, NULL) == RES_BAD_ARG);
+
+  voxdesc.get = NULL;
+  CHK(NEW_SCN(dev, low, upp, nvxls, &voxdesc, &oct) == RES_BAD_ARG);
+  voxdesc.get = get;
+
+  voxdesc.merge = NULL;
+  CHK(NEW_SCN(dev, low, upp, nvxls, &voxdesc, &oct) == RES_BAD_ARG);
+  voxdesc.merge = keep_max;
+
+  voxdesc.challenge_merge = NULL;
+  CHK(NEW_SCN(dev, low, upp, nvxls, &voxdesc, &oct) == RES_BAD_ARG);
+  voxdesc.challenge_merge = no_merge;
+
+  voxdesc.size = 0;
+  CHK(NEW_SCN(dev, low, upp, nvxls, &voxdesc, &oct) == RES_BAD_ARG);
+  voxdesc.size = sizeof(double);
+
+  voxdesc.size = SVX_MAX_SIZEOF_VOXEL + 1;
+  CHK(NEW_SCN(dev, low, upp, nvxls, &voxdesc, &oct) == RES_BAD_ARG);
+  voxdesc.size = sizeof(double);
+
+  CHK(NEW_SCN(dev, low, upp, nvxls, &voxdesc, &oct) == RES_OK);
+
+  CHK(svx_octree_for_each_leaf(oct, check_voxel, &ctx) == RES_OK);
+
+  CHK(svx_octree_get_desc(NULL, &octdesc) == RES_BAD_ARG);
+  CHK(svx_octree_get_desc(oct, NULL) == RES_BAD_ARG);
+  CHK(svx_octree_get_desc(oct, &octdesc) == RES_OK);
+  CHK(octdesc.nleaves == 5*5*5);
+  CHK(octdesc.nvoxels == octdesc.nleaves + 36/*#parents*/);
+  CHK(octdesc.depth == 4);
+
+  d3_splat(low, 0);
+  d3_splat(upp, 1);
+  CHK(d3_eq(octdesc.lower, low));
+  CHK(d3_eq(octdesc.upper, upp));
+
+  test_at_accessor(oct, nvxls);
+
+  CHK(svx_octree_ref_put(oct) == RES_OK);
+
+  nvxls[0] = nvxls[1] = nvxls[2] = 32;
+  voxdesc.challenge_merge = merge_level0;
+  CHK(NEW_SCN(dev, low, upp, nvxls, &voxdesc, &oct) == RES_OK);
+  CHK(svx_octree_get_desc(oct, &octdesc) == RES_OK);
+  CHK(octdesc.nleaves == nvxls[0]*nvxls[1]*nvxls[2] / 8);
+  CHK(octdesc.nvoxels == (octdesc.nleaves*8 - 1) / 7);
+  CHK(octdesc.depth == (size_t)log2i((int)(nvxls[0]/2))+1);
+
+  dump_data(stdout, oct);
+
+  #undef NEW_SCN
+
+  CHK(svx_device_ref_put(dev) == RES_OK);
+  CHK(svx_octree_ref_put(oct) == RES_OK);
+
+  check_memory_allocator(&allocator);
+  mem_shutdown_proxy_allocator(&allocator);
+  CHK(mem_allocated_size() == 0);
+  return 0;
+}
+
diff --git a/src/test_svx_utils.h b/src/test_svx_utils.h
@@ -0,0 +1,34 @@
+/* Copyright (C) 2018 CNRS
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. */
+
+#ifndef TEST_HTVOX_UTILS_H
+#define TEST_HTVOX_UTILS_H
+
+#include <rsys/mem_allocator.h>
+#include <stdio.h>
+
+static INLINE void
+check_memory_allocator(struct mem_allocator* allocator)
+{
+  if(MEM_ALLOCATED_SIZE(allocator)) {
+    char dump[512];
+    MEM_DUMP(allocator, dump, sizeof(dump)/sizeof(char));
+    fprintf(stderr, "%s\n", dump);
+    FATAL("Memory leaks\n");
+  }
+}
+
+#endif /* TEST_HTVOX_UTILS_H */
+