Import Upstream version 0.5.0+dfsg

18 files changed, 7808 insertions, 0 deletions

diff --git a/silx/gui/plot3d/scene/__init__.py b/silx/gui/plot3d/scene/__init__.py
new file mode 100644
index 0000000..25a7171
--- /dev/null
+++ b/silx/gui/plot3d/scene/__init__.py
@@ -0,0 +1,34 @@
+# coding: utf-8
+# /*##########################################################################
+#
+# Copyright (c) 2015-2017 European Synchrotron Radiation Facility
+#
+# Permission is hereby granted, free of charge, to any person obtaining a copy
+# of this software and associated documentation files (the "Software"), to deal
+# in the Software without restriction, including without limitation the rights
+# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+# copies of the Software, and to permit persons to whom the Software is
+# furnished to do so, subject to the following conditions:
+#
+# The above copyright notice and this permission notice shall be included in
+# all copies or substantial portions of the Software.
+#
+# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
+# THE SOFTWARE.
+#
+# ###########################################################################*/
+"""This module provides a 3D graphics scene graph structure."""
+
+__authors__ = ["T. Vincent"]
+__license__ = "MIT"
+__date__ = "08/11/2016"
+
+
+from .core import Elem, Group, PrivateGroup  # noqa
+from .viewport import Viewport  # noqa
+from .window import Window  # noqa
diff --git a/silx/gui/plot3d/scene/axes.py b/silx/gui/plot3d/scene/axes.py
new file mode 100644
index 0000000..528e4f7
--- /dev/null
+++ b/silx/gui/plot3d/scene/axes.py
@@ -0,0 +1,224 @@
+# coding: utf-8
+# /*##########################################################################
+#
+# Copyright (c) 2016-2017 European Synchrotron Radiation Facility
+#
+# Permission is hereby granted, free of charge, to any person obtaining a copy
+# of this software and associated documentation files (the "Software"), to deal
+# in the Software without restriction, including without limitation the rights
+# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+# copies of the Software, and to permit persons to whom the Software is
+# furnished to do so, subject to the following conditions:
+#
+# The above copyright notice and this permission notice shall be included in
+# all copies or substantial portions of the Software.
+#
+# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
+# THE SOFTWARE.
+#
+# ###########################################################################*/
+"""Primitive displaying a text field in the scene."""
+
+from __future__ import absolute_import, division, unicode_literals
+
+__authors__ = ["T. Vincent"]
+__license__ = "MIT"
+__date__ = "17/10/2016"
+
+
+import logging
+import numpy
+
+from ...plot._utils import ticklayout
+
+from . import core, primitives, text, transform
+
+
+_logger = logging.getLogger(__name__)
+
+
+class LabelledAxes(primitives.GroupBBox):
+    """A group displaying a bounding box with axes labels around its children.
+    """
+
+    def __init__(self):
+        super(LabelledAxes, self).__init__()
+        self._ticksForBounds = None
+
+        self._font = text.Font()
+
+        # TODO offset labels from anchor in pixels
+
+        self._xlabel = text.Text2D(font=self._font)
+        self._xlabel.align = 'center'
+        self._xlabel.transforms = [self._boxTransforms,
+                                   transform.Translate(tx=0.5)]
+        self._children.insert(-1, self._xlabel)
+
+        self._ylabel = text.Text2D(font=self._font)
+        self._ylabel.align = 'center'
+        self._ylabel.transforms = [self._boxTransforms,
+                                   transform.Translate(ty=0.5)]
+        self._children.insert(-1, self._ylabel)
+
+        self._zlabel = text.Text2D(font=self._font)
+        self._zlabel.align = 'center'
+        self._zlabel.transforms = [self._boxTransforms,
+                                   transform.Translate(tz=0.5)]
+        self._children.insert(-1, self._zlabel)
+
+        # Init tick lines with dummy pos
+        self._tickLines = primitives.DashedLines(
+            positions=((0., 0., 0.), (0., 0., 0.)))
+        self._tickLines.dash = 5, 10
+        self._tickLines.visible = False
+        self._children.insert(-1, self._tickLines)
+
+        self._tickLabels = core.Group()
+        self._children.insert(-1, self._tickLabels)
+
+        # Sync color
+        self.tickColor = 1., 1., 1., 1.
+
+    @property
+    def tickColor(self):
+        """Color of ticks and text labels.
+
+        This does NOT set bounding box color.
+        Use :attr:`color` for the bounding box.
+        """
+        return self._xlabel.foreground
+
+    @tickColor.setter
+    def tickColor(self, color):
+        self._xlabel.foreground = color
+        self._ylabel.foreground = color
+        self._zlabel.foreground = color
+        transparentColor = color[0], color[1], color[2], color[3] * 0.6
+        self._tickLines.setAttribute('color', transparentColor)
+        for label in self._tickLabels.children:
+            label.foreground = color
+
+    @property
+    def font(self):
+        """Font of axes text labels (Font)"""
+        return self._font
+
+    @font.setter
+    def font(self, font):
+        self._font = font
+        self._xlabel.font = font
+        self._ylabel.font = font
+        self._zlabel.font = font
+        for label in self._tickLabels.children:
+            label.font = font
+
+    @property
+    def xlabel(self):
+        """Text label of the X axis (str)"""
+        return self._xlabel.text
+
+    @xlabel.setter
+    def xlabel(self, text):
+        self._xlabel.text = text
+
+    @property
+    def ylabel(self):
+        """Text label of the Y axis (str)"""
+        return self._ylabel.text
+
+    @ylabel.setter
+    def ylabel(self, text):
+        self._ylabel.text = text
+
+    @property
+    def zlabel(self):
+        """Text label of the Z axis (str)"""
+        return self._zlabel.text
+
+    @zlabel.setter
+    def zlabel(self, text):
+        self._zlabel.text = text
+
+    def _updateTicks(self):
+        """Check if ticks need update and update them if needed."""
+        bounds = self._group.bounds(transformed=False, dataBounds=True)
+        if bounds is None:  # No content
+            if self._ticksForBounds is not None:
+                self._ticksForBounds = None
+                self._tickLines.visible = False
+                self._tickLabels.children = []  # Reset previous labels
+
+        elif (self._ticksForBounds is None or
+                not numpy.all(numpy.equal(bounds, self._ticksForBounds))):
+            self._ticksForBounds = bounds
+
+            # Update ticks
+            ticklength = numpy.abs(bounds[1] - bounds[0])
+
+            xticks, xlabels = ticklayout.ticks(*bounds[:, 0])
+            yticks, ylabels = ticklayout.ticks(*bounds[:, 1])
+            zticks, zlabels = ticklayout.ticks(*bounds[:, 2])
+
+            # Update tick lines
+            coords = numpy.empty(
+                ((len(xticks) + len(yticks) + len(zticks)), 4, 3),
+                dtype=numpy.float32)
+            coords[:, :, :] = bounds[0, :]  # account for offset from origin
+
+            xcoords = coords[:len(xticks)]
+            xcoords[:, :, 0] = numpy.asarray(xticks)[:, numpy.newaxis]
+            xcoords[:, 1, 1] += ticklength[1]  # X ticks on XY plane
+            xcoords[:, 3, 2] += ticklength[2]  # X ticks on XZ plane
+
+            ycoords = coords[len(xticks):len(xticks) + len(yticks)]
+            ycoords[:, :, 1] = numpy.asarray(yticks)[:, numpy.newaxis]
+            ycoords[:, 1, 0] += ticklength[0]  # Y ticks on XY plane
+            ycoords[:, 3, 2] += ticklength[2]  # Y ticks on YZ plane
+
+            zcoords = coords[len(xticks) + len(yticks):]
+            zcoords[:, :, 2] = numpy.asarray(zticks)[:, numpy.newaxis]
+            zcoords[:, 1, 0] += ticklength[0]  # Z ticks on XZ plane
+            zcoords[:, 3, 1] += ticklength[1]  # Z ticks on YZ plane
+
+            self._tickLines.setPositions(coords.reshape(-1, 3))
+            self._tickLines.visible = True
+
+            # Update labels
+            color = self.tickColor
+            offsets = bounds[0] - ticklength / 20.
+            labels = []
+            for tick, label in zip(xticks, xlabels):
+                text2d = text.Text2D(text=label, font=self.font)
+                text2d.align = 'center'
+                text2d.foreground = color
+                text2d.transforms = [transform.Translate(
+                    tx=tick, ty=offsets[1], tz=offsets[2])]
+                labels.append(text2d)
+
+            for tick, label in zip(yticks, ylabels):
+                text2d = text.Text2D(text=label, font=self.font)
+                text2d.align = 'center'
+                text2d.foreground = color
+                text2d.transforms = [transform.Translate(
+                    tx=offsets[0], ty=tick, tz=offsets[2])]
+                labels.append(text2d)
+
+            for tick, label in zip(zticks, zlabels):
+                text2d = text.Text2D(text=label, font=self.font)
+                text2d.align = 'center'
+                text2d.foreground = color
+                text2d.transforms = [transform.Translate(
+                    tx=offsets[0], ty=offsets[1], tz=tick)]
+                labels.append(text2d)
+
+            self._tickLabels.children = labels  # Reset previous labels
+
+    def prepareGL2(self, context):
+        self._updateTicks()
+        super(LabelledAxes, self).prepareGL2(context)
diff --git a/silx/gui/plot3d/scene/camera.py b/silx/gui/plot3d/scene/camera.py
new file mode 100644
index 0000000..8cc279d
--- /dev/null
+++ b/silx/gui/plot3d/scene/camera.py
@@ -0,0 +1,350 @@
+# coding: utf-8
+# /*##########################################################################
+#
+# Copyright (c) 2015-2017 European Synchrotron Radiation Facility
+#
+# Permission is hereby granted, free of charge, to any person obtaining a copy
+# of this software and associated documentation files (the "Software"), to deal
+# in the Software without restriction, including without limitation the rights
+# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+# copies of the Software, and to permit persons to whom the Software is
+# furnished to do so, subject to the following conditions:
+#
+# The above copyright notice and this permission notice shall be included in
+# all copies or substantial portions of the Software.
+#
+# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
+# THE SOFTWARE.
+#
+# ###########################################################################*/
+"""This module provides classes to handle a perspective projection in 3D."""
+
+from __future__ import absolute_import, division, unicode_literals
+
+__authors__ = ["T. Vincent"]
+__license__ = "MIT"
+__date__ = "25/07/2016"
+
+
+import numpy
+
+from . import transform
+
+
+# CameraExtrinsic #############################################################
+
+class CameraExtrinsic(transform.Transform):
+    """Transform matrix to handle camera position and orientation.
+
+    :param position: Coordinates of the point of view.
+    :type position: numpy.ndarray-like of 3 float32.
+    :param direction: Sight direction vector.
+    :type direction: numpy.ndarray-like of 3 float32.
+    :param up: Vector pointing upward in the image plane.
+    :type up: numpy.ndarray-like of 3 float32.
+    """
+
+    def __init__(self, position=(0., 0., 0.),
+                 direction=(0., 0., -1.),
+                 up=(0., 1., 0.)):
+
+        super(CameraExtrinsic, self).__init__()
+        self._position = None
+        self.position = position  # set _position
+        self._side = 1., 0., 0.
+        self._up = 0., 1., 0.
+        self._direction = 0., 0., -1.
+        self.setOrientation(direction=direction, up=up)  # set _direction, _up
+
+    def _makeMatrix(self):
+        return transform.mat4LookAtDir(self._position,
+                                       self._direction, self._up)
+
+    def copy(self):
+        """Return an independent copy"""
+        return CameraExtrinsic(self.position, self.direction, self.up)
+
+    def setOrientation(self, direction=None, up=None):
+        """Set the rotation of the point of view.
+
+        :param direction: Sight direction vector or
+                          None to keep the current one.
+        :type direction: numpy.ndarray-like of 3 float32 or None.
+        :param up: Vector pointing upward in the image plane or
+                   None to keep the current one.
+        :type up: numpy.ndarray-like of 3 float32 or None.
+        :raises RuntimeError: if the direction and up are parallel.
+        """
+        if direction is None:  # Use current direction
+            direction = self.direction
+        else:
+            assert len(direction) == 3
+            direction = numpy.array(direction, copy=True, dtype=numpy.float32)
+            direction /= numpy.linalg.norm(direction)
+
+        if up is None:  # Use current up
+            up = self.up
+        else:
+            assert len(up) == 3
+            up = numpy.array(up, copy=True, dtype=numpy.float32)
+
+        # Update side and up to make sure they are perpendicular and normalized
+        side = numpy.cross(direction, up)
+        sidenormal = numpy.linalg.norm(side)
+        if sidenormal == 0.:
+            raise RuntimeError('direction and up vectors are parallel.')
+            # Alternative: when one of the input parameter is None, it is
+            # possible to guess correct vectors using previous direction and up
+        side /= sidenormal
+        up = numpy.cross(side, direction)
+        up /= numpy.linalg.norm(up)
+
+        self._side = side
+        self._up = up
+        self._direction = direction
+        self.notify()
+
+    @property
+    def position(self):
+        """Coordinates of the point of view as a numpy.ndarray of 3 float32."""
+        return self._position.copy()
+
+    @position.setter
+    def position(self, position):
+        assert len(position) == 3
+        self._position = numpy.array(position, copy=True, dtype=numpy.float32)
+        self.notify()
+
+    @property
+    def direction(self):
+        """Sight direction (ndarray of 3 float32)."""
+        return self._direction.copy()
+
+    @direction.setter
+    def direction(self, direction):
+        self.setOrientation(direction=direction)
+
+    @property
+    def up(self):
+        """Vector pointing upward in the image plane (ndarray of 3 float32).
+        """
+        return self._up.copy()
+
+    @up.setter
+    def up(self, up):
+        self.setOrientation(up=up)
+
+    @property
+    def side(self):
+        """Vector pointing towards the side of the image plane.
+
+        ndarray of 3 float32"""
+        return self._side.copy()
+
+    def move(self, direction, step=1.):
+        """Move the camera relative to the image plane.
+
+        :param str direction: Direction relative to image plane.
+                              One of: 'up', 'down', 'left', 'right',
+                              'forward', 'backward'.
+        :param float step: The step of the pan to perform in the coordinate
+                           in which the camera position is defined.
+        """
+        if direction in ('up', 'down'):
+            vector = self.up * (1. if direction == 'up' else -1.)
+        elif direction in ('left', 'right'):
+            vector = self.side * (1. if direction == 'right' else -1.)
+        elif direction in ('forward', 'backward'):
+            vector = self.direction * (1. if direction == 'forward' else -1.)
+        else:
+            raise ValueError('Unsupported direction: %s' % direction)
+
+        self.position += step * vector
+
+    def rotate(self, direction, angle=1.):
+        """First-person rotation of the camera towards the direction.
+
+        :param str direction: Direction of movement relative to image plane.
+                              In: 'up', 'down', 'left', 'right'.
+        :param float angle: The angle in degrees of the rotation.
+        """
+        if direction in ('up', 'down'):
+            axis = self.side * (1. if direction == 'up' else -1.)
+        elif direction in ('left', 'right'):
+            axis = self.up * (1. if direction == 'left' else -1.)
+        else:
+            raise ValueError('Unsupported direction: %s' % direction)
+
+        matrix = transform.mat4RotateFromAngleAxis(numpy.radians(angle), *axis)
+        newdir = numpy.dot(matrix[:3, :3], self.direction)
+
+        if direction in ('up', 'down'):
+            # Rotate up to avoid up and new direction to be (almost) co-linear
+            newup = numpy.dot(matrix[:3, :3], self.up)
+            self.setOrientation(newdir, newup)
+        else:
+            # No need to rotate up here as it is the rotation axis
+            self.direction = newdir
+
+    def orbit(self, direction, center=(0., 0., 0.), angle=1.):
+        """Rotate the camera around a point.
+
+        :param str direction: Direction of movement relative to image plane.
+                              In: 'up', 'down', 'left', 'right'.
+        :param center: Position around which to rotate the point of view.
+        :type center: numpy.ndarray-like of 3 float32.
+        :param float angle: he angle in degrees of the rotation.
+        """
+        if direction in ('up', 'down'):
+            axis = self.side * (1. if direction == 'down' else -1.)
+        elif direction in ('left', 'right'):
+            axis = self.up * (1. if direction == 'right' else -1.)
+        else:
+            raise ValueError('Unsupported direction: %s' % direction)
+
+        # Rotate viewing direction
+        rotmatrix = transform.mat4RotateFromAngleAxis(
+            numpy.radians(angle), *axis)
+        self.direction = numpy.dot(rotmatrix[:3, :3], self.direction)
+
+        # Rotate position around center
+        center = numpy.array(center, copy=False, dtype=numpy.float32)
+        matrix = numpy.dot(transform.mat4Translate(*center), rotmatrix)
+        matrix = numpy.dot(matrix, transform.mat4Translate(*(-center)))
+        position = numpy.append(self.position, 1.)
+        self.position = numpy.dot(matrix, position)[:3]
+
+    _RESET_CAMERA_ORIENTATIONS = {
+        'side': ((-1., -1., -1.), (0., 1., 0.)),
+        'front': ((0., 0., -1.), (0., 1., 0.)),
+        'back': ((0., 0., 1.), (0., 1., 0.)),
+        'top': ((0., -1., 0.), (0., 0., -1.)),
+        'bottom': ((0., 1., 0.), (0., 0., 1.)),
+        'right': ((-1., 0., 0.), (0., 1., 0.)),
+        'left': ((1., 0., 0.), (0., 1., 0.))
+    }
+
+    def reset(self, face=None):
+        """Reset the camera position to pre-defined orientations.
+
+        :param str face: The direction of the camera in:
+                         side, front, back, top, bottom, right, left.
+        """
+        if face not in self._RESET_CAMERA_ORIENTATIONS:
+            raise ValueError('Unsupported face: %s' % face)
+
+        distance = numpy.linalg.norm(self.position)
+        direction, up = self._RESET_CAMERA_ORIENTATIONS[face]
+        self.setOrientation(direction, up)
+        self.position = - self.direction * distance
+
+
+class Camera(transform.Transform):
+    """Combination of camera projection and position.
+
+    See :class:`Perspective` and :class:`CameraExtrinsic`.
+
+    :param float fovy: Vertical field-of-view in degrees.
+    :param float near: The near clipping plane Z coord (strictly positive).
+    :param float far: The far clipping plane Z coord (> near).
+    :param size: Viewport's size used to compute the aspect ratio.
+    :type size: 2-tuple of float (width, height).
+    :param position: Coordinates of the point of view.
+    :type position: numpy.ndarray-like of 3 float32.
+    :param direction: Sight direction vector.
+    :type direction: numpy.ndarray-like of 3 float32.
+    :param up: Vector pointing upward in the image plane.
+    :type up: numpy.ndarray-like of 3 float32.
+    """
+
+    def __init__(self, fovy=30., near=0.1, far=1., size=(1., 1.),
+                 position=(0., 0., 0.),
+                 direction=(0., 0., -1.), up=(0., 1., 0.)):
+        super(Camera, self).__init__()
+        self._intrinsic = transform.Perspective(fovy, near, far, size)
+        self._intrinsic.addListener(self._transformChanged)
+        self._extrinsic = CameraExtrinsic(position, direction, up)
+        self._extrinsic.addListener(self._transformChanged)
+
+    def _makeMatrix(self):
+        return numpy.dot(self.intrinsic.matrix, self.extrinsic.matrix)
+
+    def _transformChanged(self, source):
+        """Listener of intrinsic and extrinsic camera parameters instances."""
+        if source is not self:
+            self.notify()
+
+    def resetCamera(self, bounds):
+        """Change camera to have the bounds in the viewing frustum.
+
+        It updates the camera position and depth extent.
+        Camera sight direction and up are not affected.
+
+        :param bounds: The axes-aligned bounds to include.
+        :type bounds: numpy.ndarray: ((xMin, yMin, zMin), (xMax, yMax, zMax))
+        """
+
+        center = 0.5 * (bounds[0] + bounds[1])
+        radius = numpy.linalg.norm(0.5 * (bounds[1] - bounds[0]))
+
+        if isinstance(self.intrinsic, transform.Perspective):
+            # Get the viewpoint distance from the bounds center
+            minfov = numpy.radians(self.intrinsic.fovy)
+            width, height = self.intrinsic.size
+            if width < height:
+                minfov *= width / height
+
+            offset = radius / numpy.sin(0.5 * minfov)
+
+            # Update camera
+            self.extrinsic.position = \
+                center - offset * self.extrinsic.direction
+            self.intrinsic.setDepthExtent(offset - radius, offset + radius)
+
+        elif isinstance(self.intrinsic, transform.Orthographic):
+            # Y goes up
+            self.intrinsic.setClipping(
+                left=center[0] - radius,
+                right=center[0] + radius,
+                bottom=center[1] - radius,
+                top=center[1] + radius)
+
+            # Update camera
+            self.extrinsic.position = 0, 0, 0
+            self.intrinsic.setDepthExtent(center[2] - radius,
+                                          center[2] + radius)
+        else:
+            raise RuntimeError('Unsupported camera: %s' % self.intrinsic)
+
+    @property
+    def intrinsic(self):
+        """Intrinsic camera parameters, i.e., projection matrix."""
+        return self._intrinsic
+
+    @intrinsic.setter
+    def intrinsic(self, intrinsic):
+        self._intrinsic.removeListener(self._transformChanged)
+        self._intrinsic = intrinsic
+        self._intrinsic.addListener(self._transformChanged)
+
+    @property
+    def extrinsic(self):
+        """Extrinsic camera parameters, i.e., position and orientation."""
+        return self._extrinsic
+
+    def move(self, *args, **kwargs):
+        """See :meth:`CameraExtrinsic.move`."""
+        self.extrinsic.move(*args, **kwargs)
+
+    def rotate(self, *args, **kwargs):
+        """See :meth:`CameraExtrinsic.rotate`."""
+        self.extrinsic.rotate(*args, **kwargs)
+
+    def orbit(self, *args, **kwargs):
+        """See :meth:`CameraExtrinsic.orbit`."""
+        self.extrinsic.orbit(*args, **kwargs)
diff --git a/silx/gui/plot3d/scene/core.py b/silx/gui/plot3d/scene/core.py
new file mode 100644
index 0000000..a293f28
--- /dev/null
+++ b/silx/gui/plot3d/scene/core.py
@@ -0,0 +1,334 @@
+# coding: utf-8
+# /*##########################################################################
+#
+# Copyright (c) 2015-2017 European Synchrotron Radiation Facility
+#
+# Permission is hereby granted, free of charge, to any person obtaining a copy
+# of this software and associated documentation files (the "Software"), to deal
+# in the Software without restriction, including without limitation the rights
+# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+# copies of the Software, and to permit persons to whom the Software is
+# furnished to do so, subject to the following conditions:
+#
+# The above copyright notice and this permission notice shall be included in
+# all copies or substantial portions of the Software.
+#
+# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
+# THE SOFTWARE.
+#
+# ###########################################################################*/
+"""This module provides the base scene structure.
+
+This module provides the classes for describing a tree structure with
+rendering and picking API.
+All nodes inherit from :class:`Base`.
+Nodes with children are provided with :class:`PrivateGroup` and
+:class:`Group` classes.
+Leaf rendering nodes should inherit from :class:`Elem`.
+"""
+
+from __future__ import absolute_import, division, unicode_literals
+
+__authors__ = ["T. Vincent"]
+__license__ = "MIT"
+__date__ = "25/07/2016"
+
+
+import itertools
+import weakref
+
+import numpy
+
+from . import event
+from . import transform
+
+from .viewport import Viewport
+
+
+# Nodes #######################################################################
+
+class Base(event.Notifier):
+    """A scene node with common features."""
+
+    def __init__(self):
+        super(Base, self).__init__()
+        self._visible = True
+        self._pickable = False
+
+        self._parentRef = None
+
+        self._transforms = transform.TransformList()
+        self._transforms.addListener(self._transformChanged)
+
+    # notifying properties
+
+    visible = event.notifyProperty('_visible',
+                                   doc="Visibility flag of the node")
+    pickable = event.notifyProperty('_pickable',
+                                    doc="True to make node pickable")
+
+    # Access to tree path
+
+    @property
+    def parent(self):
+        """Parent or None if no parent"""
+        return None if self._parentRef is None else self._parentRef()
+
+    def _setParent(self, parent):
+        """Set the parent of this node.
+
+        For internal use.
+
+        :param Base parent: The parent.
+        """
+        if parent is not None and self._parentRef is not None:
+            raise RuntimeError('Trying to add a node at two places.')
+            # Alternative: remove it from previous children list
+        self._parentRef = None if parent is None else weakref.ref(parent)
+
+    @property
+    def path(self):
+        """Tuple of scene nodes, from the tip of the tree down to this node.
+
+        If this tree is attached to a :class:`Viewport`,
+        then the :class:`Viewport` is the first element of path.
+        """
+        if self.parent is None:
+            return self,
+        elif isinstance(self.parent, Viewport):
+            return self.parent, self
+        else:
+            return self.parent.path + (self, )
+
+    @property
+    def viewport(self):
+        """The viewport this node is attached to or None."""
+        root = self.path[0]
+        return root if isinstance(root, Viewport) else None
+
+    @property
+    def objectToNDCTransform(self):
+        """Transform from object to normalized device coordinates.
+
+        Do not forget perspective divide.
+        """
+        # Using the Viewport's transforms property to proxy the camera
+        path = self.path
+        assert isinstance(path[0], Viewport)
+        return transform.StaticTransformList(elem.transforms for elem in path)
+
+    @property
+    def objectToSceneTransform(self):
+        """Transform from object to scene.
+
+        Combine transforms up to the Viewport (not including it).
+        """
+        path = self.path
+        if isinstance(path[0], Viewport):
+            path = path[1:]  # Remove viewport to remove camera transforms
+        return transform.StaticTransformList(elem.transforms for elem in path)
+
+    # transform
+
+    @property
+    def transforms(self):
+        """List of transforms defining the frame of this node relative
+        to its parent."""
+        return self._transforms
+
+    @transforms.setter
+    def transforms(self, iterable):
+        self._transforms.removeListener(self._transformChanged)
+        if isinstance(iterable, transform.TransformList):
+            # If it is a TransformList, do not create one to enable sharing.
+            self._transforms = iterable
+        else:
+            assert hasattr(iterable, '__iter__')
+            self._transforms = transform.TransformList(iterable)
+        self._transforms.addListener(self._transformChanged)
+
+    def _transformChanged(self, source):
+        self.notify()  # Broadcast transform notification
+
+    # Bounds
+
+    _CUBE_CORNERS = numpy.array(list(itertools.product((0., 1.), repeat=3)),
+                                dtype=numpy.float32)
+    """Unit cube corners used to transform bounds"""
+
+    def _bounds(self, dataBounds=False):
+        """Override in subclass to provide bounds in object coordinates"""
+        return None
+
+    def bounds(self, transformed=False, dataBounds=False):
+        """Returns the bounds of this node aligned with the axis,
+        with or without transform applied.
+
+        :param bool transformed: False to give bounds in object coordinates
+                                 (the default), True to apply this object's
+                                 transforms.
+        :param bool dataBounds: False to give bounds of vertices (the default),
+                                True to give bounds of the represented data.
+        :return: The bounds: ((xMin, yMin, zMin), (xMax, yMax, zMax)) or None
+                 if no bounds.
+        :rtype: numpy.ndarray of float
+        """
+        bounds = self._bounds(dataBounds)
+
+        if transformed and bounds is not None:
+            bounds = self.transforms.transformBounds(bounds)
+
+        return bounds
+
+    # Rendering
+
+    def prepareGL2(self, ctx):
+        """Called before the rendering to prepare OpenGL resources.
+
+        Override in subclass.
+        """
+        pass
+
+    def renderGL2(self, ctx):
+        """Called to perform the OpenGL rendering.
+
+        Override in subclass.
+        """
+        pass
+
+    def render(self, ctx):
+        """Called internally to perform rendering."""
+        if self.visible:
+            ctx.pushTransform(self.transforms)
+            self.prepareGL2(ctx)
+            self.renderGL2(ctx)
+            ctx.popTransform()
+
+    def postRender(self, ctx):
+        """Hook called when parent's node render is finished.
+
+        Called in the reverse of rendering order (i.e., last child first).
+
+        Meant for nodes that modify the :class:`RenderContext` ctx to
+        reset their modifications.
+        """
+        pass
+
+    def pick(self, ctx, x, y, depth=None):
+        """True/False picking, should be fast"""
+        if self.pickable:
+            pass
+
+    def pickRay(self, ctx, ray):
+        """Picking returning list of ray intersections."""
+        if self.pickable:
+            pass
+
+
+class Elem(Base):
+    """A scene node that does some rendering."""
+
+    def __init__(self):
+        super(Elem, self).__init__()
+        # self.showBBox = False  # Here or outside scene graph?
+        # self.clipPlane = None  # This needs to be handled in the shader
+
+
+class PrivateGroup(Base):
+    """A scene node that renders its (private) childern.
+
+    :param iterable children: :class:`Base` nodes to add as children
+    """
+
+    class ChildrenList(event.NotifierList):
+        """List of children with notification and children's parent update."""
+
+        def _listWillChangeHook(self, methodName, *args, **kwargs):
+            super(PrivateGroup.ChildrenList, self)._listWillChangeHook(
+                methodName, *args, **kwargs)
+            for item in self:
+                item._setParent(None)
+
+        def _listWasChangedHook(self, methodName, *args, **kwargs):
+            for item in self:
+                item._setParent(self._parentRef())
+            super(PrivateGroup.ChildrenList, self)._listWasChangedHook(
+                methodName, *args, **kwargs)
+
+        def __init__(self, parent, children):
+            self._parentRef = weakref.ref(parent)
+            super(PrivateGroup.ChildrenList, self).__init__(children)
+
+    def __init__(self, children=()):
+        super(PrivateGroup, self).__init__()
+        self.__children = PrivateGroup.ChildrenList(self, children)
+        self.__children.addListener(self._updated)
+
+    @property
+    def _children(self):
+        """List of children to be rendered.
+
+        This private attribute is meant to be used by subclass.
+        """
+        return self.__children
+
+    @_children.setter
+    def _children(self, iterable):
+        self.__children.removeListener(self._updated)
+        for item in self.__children:
+            item._setParent(None)
+        del self.__children  # This is needed
+        self.__children = PrivateGroup.ChildrenList(self, iterable)
+        self.__children.addListener(self._updated)
+        self.notify()
+
+    def _updated(self, source, *args, **kwargs):
+        """Listen for updates"""
+        if source is not self:  # Avoid infinite recursion
+            self.notify(*args, **kwargs)
+
+    def _bounds(self, dataBounds=False):
+        """Compute the bounds from transformed children bounds"""
+        bounds = []
+        for child in self._children:
+            if child.visible:
+                childBounds = child.bounds(
+                    transformed=True, dataBounds=dataBounds)
+                if childBounds is not None:
+                    bounds.append(childBounds)
+
+        if len(bounds) == 0:
+            return None
+        else:
+            bounds = numpy.array(bounds, dtype=numpy.float32)
+            return numpy.array((bounds[:, 0, :].min(axis=0),
+                                bounds[:, 1, :].max(axis=0)),
+                               dtype=numpy.float32)
+
+    def prepareGL2(self, ctx):
+        pass
+
+    def renderGL2(self, ctx):
+        """Render all children"""
+        for child in self._children:
+            child.render(ctx)
+        for child in reversed(self._children):
+            child.postRender(ctx)
+
+
+class Group(PrivateGroup):
+    """A scene node that renders its (public) children."""
+
+    @property
+    def children(self):
+        """List of children to be rendered."""
+        return self._children
+
+    @children.setter
+    def children(self, iterable):
+        self._children = iterable
diff --git a/silx/gui/plot3d/scene/cutplane.py b/silx/gui/plot3d/scene/cutplane.py
new file mode 100644
index 0000000..79b4168
--- /dev/null
+++ b/silx/gui/plot3d/scene/cutplane.py
@@ -0,0 +1,374 @@
+# coding: utf-8
+# /*##########################################################################
+#
+# Copyright (c) 2016-2017 European Synchrotron Radiation Facility
+#
+# Permission is hereby granted, free of charge, to any person obtaining a copy
+# of this software and associated documentation files (the "Software"), to deal
+# in the Software without restriction, including without limitation the rights
+# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+# copies of the Software, and to permit persons to whom the Software is
+# furnished to do so, subject to the following conditions:
+#
+# The above copyright notice and this permission notice shall be included in
+# all copies or substantial portions of the Software.
+#
+# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
+# THE SOFTWARE.
+#
+# ###########################################################################*/
+"""A cut plane in a 3D texture: hackish implementation...
+"""
+
+from __future__ import absolute_import, division, unicode_literals
+
+__authors__ = ["T. Vincent"]
+__license__ = "MIT"
+__date__ = "05/10/2016"
+
+import string
+import numpy
+
+from ... import _glutils
+from ..._glutils import gl
+
+from .function import Colormap
+from .primitives import Box, Geometry, PlaneInGroup
+from . import transform, utils
+
+
+class ColormapMesh3D(Geometry):
+    """A 3D mesh with color from a 3D texture."""
+
+    _shaders = ("""
+    attribute vec3 position;
+    attribute vec3 normal;
+
+    uniform mat4 matrix;
+    uniform mat4 transformMat;
+    //uniform mat3 matrixInvTranspose;
+    uniform vec3 dataScale;
+
+    varying vec4 vCameraPosition;
+    varying vec3 vPosition;
+    varying vec3 vNormal;
+    varying vec3 vTexCoords;
+
+    void main(void)
+    {
+        vCameraPosition = transformMat * vec4(position, 1.0);
+        //vNormal = matrixInvTranspose * normalize(normal);
+        vPosition = position;
+        vTexCoords = dataScale * position;
+        vNormal = normal;
+        gl_Position = matrix * vec4(position, 1.0);
+    }
+    """,
+                string.Template("""
+    varying vec4 vCameraPosition;
+    varying vec3 vPosition;
+    varying vec3 vNormal;
+    varying vec3 vTexCoords;
+    uniform sampler3D data;
+    uniform float alpha;
+
+    $colormapDecl
+
+    $clippingDecl
+    $lightingFunction
+
+    void main(void)
+    {
+        float value = texture3D(data, vTexCoords).r;
+        vec4 color = $colormapCall(value);
+        color.a = alpha;
+
+        $clippingCall(vCameraPosition);
+
+        gl_FragColor = $lightingCall(color, vPosition, vNormal);
+    }
+    """))
+
+    def __init__(self, position, normal, data, copy=True,
+                 mode='triangles', indices=None, colormap=None):
+        assert mode in self._TRIANGLE_MODES
+        data = numpy.array(data, copy=copy, order='C')
+        assert data.ndim == 3
+        self._data = data
+        self._texture = None
+        self._update_texture = True
+        self._update_texture_filter = False
+        self._alpha = 1.
+        self._colormap = colormap or Colormap()  # Default colormap
+        self._colormap.addListener(self._cmapChanged)
+        self._interpolation = 'linear'
+        super(ColormapMesh3D, self).__init__(mode,
+                                             indices,
+                                             position=position,
+                                             normal=normal)
+
+        self.isBackfaceVisible = True
+
+    def setData(self, data, copy=True):
+        data = numpy.array(data, copy=copy, order='C')
+        assert data.ndim == 3
+        self._data = data
+        self._update_texture = True
+
+    def getData(self, copy=True):
+        return numpy.array(self._data, copy=copy)
+
+    @property
+    def interpolation(self):
+        """The texture interpolation mode: 'linear' or 'nearest'"""
+        return self._interpolation
+
+    @interpolation.setter
+    def interpolation(self, interpolation):
+        assert interpolation in ('linear', 'nearest')
+        self._interpolation = interpolation
+        self._update_texture_filter = True
+        self.notify()
+
+    @property
+    def alpha(self):
+        """Transparency of the plane, float in [0, 1]"""
+        return self._alpha
+
+    @alpha.setter
+    def alpha(self, alpha):
+        self._alpha = float(alpha)
+
+    @property
+    def colormap(self):
+        """The colormap used by this primitive"""
+        return self._colormap
+
+    def _cmapChanged(self, source, *args, **kwargs):
+        """Broadcast colormap changes"""
+        self.notify(*args, **kwargs)
+
+    def prepareGL2(self, ctx):
+        if self._texture is None or self._update_texture:
+            if self._texture is not None:
+                self._texture.discard()
+
+            if self.interpolation == 'nearest':
+                filter_ = gl.GL_NEAREST
+            else:
+                filter_ = gl.GL_LINEAR
+            self._update_texture = False
+            self._update_texture_filter = False
+            self._texture = _glutils.Texture(
+                gl.GL_R32F, self._data, gl.GL_RED,
+                minFilter=filter_,
+                magFilter=filter_,
+                wrap=gl.GL_CLAMP_TO_EDGE)
+
+        if self._update_texture_filter:
+            self._update_texture_filter = False
+            if self.interpolation == 'nearest':
+                filter_ = gl.GL_NEAREST
+            else:
+                filter_ = gl.GL_LINEAR
+            self._texture.minFilter = filter_
+            self._texture.magFilter = filter_
+
+        super(ColormapMesh3D, self).prepareGL2(ctx)
+
+    def renderGL2(self, ctx):
+        fragment = self._shaders[1].substitute(
+            clippingDecl=ctx.clipper.fragDecl,
+            clippingCall=ctx.clipper.fragCall,
+            lightingFunction=ctx.viewport.light.fragmentDef,
+            lightingCall=ctx.viewport.light.fragmentCall,
+            colormapDecl=self.colormap.decl,
+            colormapCall=self.colormap.call
+            )
+        program = ctx.glCtx.prog(self._shaders[0], fragment)
+        program.use()
+
+        ctx.viewport.light.setupProgram(ctx, program)
+        self.colormap.setupProgram(ctx, program)
+
+        if not self.isBackfaceVisible:
+            gl.glCullFace(gl.GL_BACK)
+            gl.glEnable(gl.GL_CULL_FACE)
+
+        program.setUniformMatrix('matrix', ctx.objectToNDC.matrix)
+        program.setUniformMatrix('transformMat',
+                                 ctx.objectToCamera.matrix,
+                                 safe=True)
+        gl.glUniform1f(program.uniforms['alpha'], self._alpha)
+
+        shape = self._data.shape
+        scales = 1./shape[2], 1./shape[1], 1./shape[0]
+        gl.glUniform3f(program.uniforms['dataScale'], *scales)
+
+        gl.glUniform1i(program.uniforms['data'], self._texture.texUnit)
+
+        ctx.clipper.setupProgram(ctx, program)
+
+        self._texture.bind()
+        self._draw(program)
+
+        if not self.isBackfaceVisible:
+            gl.glDisable(gl.GL_CULL_FACE)
+
+
+class CutPlane(PlaneInGroup):
+    """A cutting plane in a 3D texture"""
+
+    def __init__(self, point=(0., 0., 0.), normal=(0., 0., 1.)):
+        self._data = None
+        self._mesh = None
+        self._alpha = 1.
+        self._interpolation = 'linear'
+        self._colormap = Colormap()
+        super(CutPlane, self).__init__(point, normal)
+
+    def setData(self, data, copy=True):
+        if data is None:
+            self._data = None
+            if self._mesh is not None:
+                self._children.remove(self._mesh)
+            self._mesh = None
+
+        else:
+            data = numpy.array(data, copy=copy, order='C')
+            assert data.ndim == 3
+            self._data = data
+            if self._mesh is not None:
+                self._mesh.setData(data, copy=False)
+
+    def getData(self, copy=True):
+        return None if self._mesh is None else self._mesh.getData(copy=copy)
+
+    @property
+    def alpha(self):
+        return self._alpha
+
+    @alpha.setter
+    def alpha(self, alpha):
+        self._alpha = float(alpha)
+        if self._mesh is not None:
+            self._mesh.alpha = alpha
+
+    @property
+    def colormap(self):
+        return self._colormap
+
+    @property
+    def interpolation(self):
+        """The texture interpolation mode: 'linear' (default) or 'nearest'"""
+        return self._interpolation
+
+    @interpolation.setter
+    def interpolation(self, interpolation):
+        assert interpolation in ('nearest', 'linear')
+        if interpolation != self.interpolation:
+            self._interpolation = interpolation
+            if self._mesh is not None:
+                self._mesh.interpolation = interpolation
+
+    def prepareGL2(self, ctx):
+        if self.isValid:
+
+            contourVertices = self.contourVertices
+
+            if (self.interpolation == 'nearest' and
+                    contourVertices is not None and len(contourVertices)):
+                # Avoid cut plane co-linear with array bin edges
+                for index, normal in enumerate(((1., 0., 0.), (0., 1., 0.), (0., 0., 1.))):
+                    if (numpy.all(numpy.equal(self.plane.normal, normal)) and
+                            int(self.plane.point[index]) == self.plane.point[index]):
+                        contourVertices += self.plane.normal * 0.01  # Add an offset
+                        break
+
+            if self._mesh is None and self._data is not None:
+                self._mesh = ColormapMesh3D(contourVertices,
+                                            normal=self.plane.normal,
+                                            data=self._data,
+                                            copy=False,
+                                            mode='fan',
+                                            colormap=self.colormap)
+                self._mesh.alpha = self._alpha
+                self._interpolation = self.interpolation
+                self._children.insert(0, self._mesh)
+
+            if self._mesh is not None:
+                if (contourVertices is None or
+                        len(contourVertices) == 0):
+                    self._mesh.visible = False
+                else:
+                    self._mesh.visible = True
+                    self._mesh.setAttribute('normal', self.plane.normal)
+                    self._mesh.setAttribute('position', contourVertices)
+
+        super(CutPlane, self).prepareGL2(ctx)
+
+    def renderGL2(self, ctx):
+        with self.viewport.light.turnOff():
+            super(CutPlane, self).renderGL2(ctx)
+
+    def _bounds(self, dataBounds=False):
+        if not dataBounds:
+            vertices = self.contourVertices
+            if vertices is not None:
+                return numpy.array(
+                    (vertices.min(axis=0), vertices.max(axis=0)),
+                    dtype=numpy.float32)
+            else:
+                return None  # Plane in not slicing the data volume
+        else:
+            if self._data is None:
+                return None
+            else:
+                depth, height, width = self._data.shape
+                return numpy.array(((0., 0., 0.),
+                                    (width, height, depth)),
+                                   dtype=numpy.float32)
+
+    @property
+    def contourVertices(self):
+        """The vertices of the contour of the plane/bounds intersection."""
+        # TODO copy from PlaneInGroup, refactor all that!
+        bounds = self.bounds(dataBounds=True)
+        if bounds is None:
+            return None  # No bounds: no vertices
+
+        # Check if cache is valid and return it
+        cachebounds, cachevertices = self._cache
+        if numpy.all(numpy.equal(bounds, cachebounds)):
+            return cachevertices
+
+        # Cache is not OK, rebuild it
+        boxvertices = bounds[0] + Box._vertices.copy()*(bounds[1] - bounds[0])
+        lineindices = Box._lineIndices
+        vertices = utils.boxPlaneIntersect(
+            boxvertices, lineindices, self.plane.normal, self.plane.point)
+
+        self._cache = bounds, vertices if len(vertices) != 0 else None
+
+        return self._cache[1]
+
+    # Render transforms RW, TODO refactor this!
+    @property
+    def transforms(self):
+        return self._transforms
+
+    @transforms.setter
+    def transforms(self, iterable):
+        self._transforms.removeListener(self._transformChanged)
+        if isinstance(iterable, transform.TransformList):
+            # If it is a TransformList, do not create one to enable sharing.
+            self._transforms = iterable
+        else:
+            assert hasattr(iterable, '__iter__')
+            self._transforms = transform.TransformList(iterable)
+        self._transforms.addListener(self._transformChanged)
diff --git a/silx/gui/plot3d/scene/event.py b/silx/gui/plot3d/scene/event.py
new file mode 100644
index 0000000..7b85434
--- /dev/null
+++ b/silx/gui/plot3d/scene/event.py
@@ -0,0 +1,225 @@
+# coding: utf-8
+# /*##########################################################################
+#
+# Copyright (c) 2015-2017 European Synchrotron Radiation Facility
+#
+# Permission is hereby granted, free of charge, to any person obtaining a copy
+# of this software and associated documentation files (the "Software"), to deal
+# in the Software without restriction, including without limitation the rights
+# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+# copies of the Software, and to permit persons to whom the Software is
+# furnished to do so, subject to the following conditions:
+#
+# The above copyright notice and this permission notice shall be included in
+# all copies or substantial portions of the Software.
+#
+# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
+# THE SOFTWARE.
+#
+# ###########################################################################*/
+"""This module provides a simple generic notification system."""
+
+from __future__ import absolute_import, division, unicode_literals
+
+__authors__ = ["T. Vincent"]
+__license__ = "MIT"
+__date__ = "25/07/2016"
+
+
+import logging
+
+from silx.utils.weakref import WeakList
+
+_logger = logging.getLogger(__name__)
+
+
+# Notifier ####################################################################
+
+class Notifier(object):
+    """Base class for object with notification mechanism."""
+
+    def __init__(self):
+        self._listeners = WeakList()
+
+    def addListener(self, listener):
+        """Register a listener.
+
+        Adding an already registered listener has no effect.
+
+        :param callable listener: The function or method to register.
+        """
+        if listener not in self._listeners:
+            self._listeners.append(listener)
+        else:
+            _logger.warning('Ignoring addition of an already registered listener')
+
+    def removeListener(self, listener):
+        """Remove a previously registered listener.
+
+        :param callable listener: The function or method to unregister.
+        """
+        try:
+            self._listeners.remove(listener)
+        except ValueError:
+            _logger.warn('Trying to remove a listener that is not registered')
+
+    def notify(self, *args, **kwargs):
+        """Notify all registered listeners with the given parameters.
+
+        Listeners are called directly in this method.
+        Listeners are called in the order they were registered.
+        """
+        for listener in self._listeners:
+            listener(self, *args, **kwargs)
+
+
+def notifyProperty(attrName, copy=False, converter=None, doc=None):
+    """Create a property that adds notification to an attribute.
+
+    :param str attrName: The name of the attribute to wrap.
+    :param bool copy: Whether to return a copy of the attribute
+                      or not (the default).
+    :param converter: Function converting input value to appropriate type
+                      This function takes a single argument and return the
+                      converted value.
+                      It can be used to perform some asserts.
+    :param str doc: The docstring of the property
+    :return: A property with getter and setter
+    """
+    if copy:
+        def getter(self):
+            return getattr(self, attrName).copy()
+    else:
+        def getter(self):
+            return getattr(self, attrName)
+
+    if converter is None:
+        def setter(self, value):
+            if getattr(self, attrName) != value:
+                setattr(self, attrName, value)
+                self.notify()
+
+    else:
+        def setter(self, value):
+            value = converter(value)
+            if getattr(self, attrName) != value:
+                setattr(self, attrName, value)
+                self.notify()
+
+    return property(getter, setter, doc=doc)
+
+
+class HookList(list):
+    """List with hooks before and after modification."""
+
+    def __init__(self, iterable):
+        super(HookList, self).__init__(iterable)
+
+        self._listWasChangedHook('__init__', iterable)
+
+    def _listWillChangeHook(self, methodName, *args, **kwargs):
+        """To override. Called before modifying the list.
+
+        This method is called with the name of the method called to
+        modify the list and its parameters.
+        """
+        pass
+
+    def _listWasChangedHook(self, methodName, *args, **kwargs):
+        """To override. Called after modifying the list.
+
+        This method is called with the name of the method called to
+        modify the list and its parameters.
+        """
+        pass
+
+    # Wrapping methods that modify the list
+
+    def _wrapper(self, methodName, *args, **kwargs):
+        """Generic wrapper of list methods calling the hooks."""
+        self._listWillChangeHook(methodName, *args, **kwargs)
+        result = getattr(super(HookList, self),
+                         methodName)(*args, **kwargs)
+        self._listWasChangedHook(methodName, *args, **kwargs)
+        return result
+
+    # Add methods
+
+    def __iadd__(self, *args, **kwargs):
+        return self._wrapper('__iadd__', *args, **kwargs)
+
+    def __imul__(self, *args, **kwargs):
+        return self._wrapper('__imul__', *args, **kwargs)
+
+    def append(self, *args, **kwargs):
+        return self._wrapper('append', *args, **kwargs)
+
+    def extend(self, *args, **kwargs):
+        return self._wrapper('extend', *args, **kwargs)
+
+    def insert(self, *args, **kwargs):
+        return self._wrapper('insert', *args, **kwargs)
+
+    # Remove methods
+
+    def __delitem__(self, *args, **kwargs):
+        return self._wrapper('__delitem__', *args, **kwargs)
+
+    def __delslice__(self, *args, **kwargs):
+        return self._wrapper('__delslice__', *args, **kwargs)
+
+    def remove(self, *args, **kwargs):
+        return self._wrapper('remove', *args, **kwargs)
+
+    def pop(self, *args, **kwargs):
+        return self._wrapper('pop', *args, **kwargs)
+
+    # Set methods
+
+    def __setitem__(self, *args, **kwargs):
+        return self._wrapper('__setitem__', *args, **kwargs)
+
+    def __setslice__(self, *args, **kwargs):
+        return self._wrapper('__setslice__', *args, **kwargs)
+
+    # In place methods
+
+    def sort(self, *args, **kwargs):
+        return self._wrapper('sort', *args, **kwargs)
+
+    def reverse(self, *args, **kwargs):
+        return self._wrapper('reverse', *args, **kwargs)
+
+
+class NotifierList(HookList, Notifier):
+    """List of Notifiers with notification mechanism.
+
+    This class registers itself as a listener of the list items.
+
+    The default listener method forward notification from list items
+    to the listeners of the list.
+    """
+
+    def __init__(self, iterable=()):
+        Notifier.__init__(self)
+        HookList.__init__(self, iterable)
+
+    def _listWillChangeHook(self, methodName, *args, **kwargs):
+        for item in self:
+            item.removeListener(self._notified)
+
+    def _listWasChangedHook(self, methodName, *args, **kwargs):
+        for item in self:
+            item.addListener(self._notified)
+        self.notify()
+
+    def _notified(self, source, *args, **kwargs):
+        """Default listener forwarding list item changes to its listeners."""
+        # Avoid infinite recursion if the list is listening itself
+        if source is not self:
+            self.notify(*args, **kwargs)
diff --git a/silx/gui/plot3d/scene/function.py b/silx/gui/plot3d/scene/function.py
new file mode 100644
index 0000000..80ac820
--- /dev/null
+++ b/silx/gui/plot3d/scene/function.py
@@ -0,0 +1,471 @@
+# coding: utf-8
+# /*##########################################################################
+#
+# Copyright (c) 2015-2017 European Synchrotron Radiation Facility
+#
+# Permission is hereby granted, free of charge, to any person obtaining a copy
+# of this software and associated documentation files (the "Software"), to deal
+# in the Software without restriction, including without limitation the rights
+# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+# copies of the Software, and to permit persons to whom the Software is
+# furnished to do so, subject to the following conditions:
+#
+# The above copyright notice and this permission notice shall be included in
+# all copies or substantial portions of the Software.
+#
+# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
+# THE SOFTWARE.
+#
+# ###########################################################################*/
+"""This module provides functions to add to shaders."""
+
+from __future__ import absolute_import, division, unicode_literals
+
+__authors__ = ["T. Vincent"]
+__license__ = "MIT"
+__date__ = "08/11/2016"
+
+
+import contextlib
+import logging
+import numpy
+
+from ..._glutils import gl
+
+from . import event
+from . import utils
+
+
+_logger = logging.getLogger(__name__)
+
+
+class ProgramFunction(object):
+    """Class providing a function to add to a GLProgram shaders.
+    """
+
+    def setupProgram(self, context, program):
+        """Sets-up uniforms of a program using this shader function.
+
+        :param RenderContext context: The current rendering context
+        :param GLProgram program: The program to set-up.
+                                  It MUST be in use and using this function.
+        """
+        pass
+
+
+class ClippingPlane(ProgramFunction):
+    """Description of a clipping plane and rendering.
+
+    Convention: Clipping is performed in camera/eye space.
+
+    :param point: Local coordinates of a point on the plane.
+    :type point: numpy.ndarray-like of 3 float32
+    :param normal: Local coordinates of the plane normal.
+    :type normal: numpy.ndarray-like of 3 float32
+    """
+
+    _fragDecl = """
+    /* Clipping plane */
+    /* as rx + gy + bz + a > 0, clipping all positive */
+    uniform vec4 planeEq;
+
+    /* Position is in camera/eye coordinates */
+
+    bool isClipped(vec4 position) {
+        vec4 tmp = planeEq * position;
+        float value = tmp.x + tmp.y + tmp.z + planeEq.a;
+        return (value < 0.0001);
+    }
+
+    void clipping(vec4 position) {
+        if (isClipped(position)) {
+            discard;
+        }
+    }
+    /* End of clipping */
+    """
+
+    _fragDeclNoop = """
+    bool isClipped(vec4 position)
+    {
+        return false;
+    }
+
+    void clipping(vec4 position) {}
+    """
+
+    def __init__(self, point=(0., 0., 0.), normal=(0., 0., 0.)):
+        self._plane = utils.Plane(point, normal)
+
+    @property
+    def plane(self):
+        """Plane parameters in camera space."""
+        return self._plane
+
+    # GL2
+
+    @property
+    def fragDecl(self):
+        return self._fragDecl if self.plane.isPlane else self._fragDeclNoop
+
+    @property
+    def fragCall(self):
+        return "clipping"
+
+    def setupProgram(self, context, program):
+        """Sets-up uniforms of a program using this shader function.
+
+        :param RenderContext context: The current rendering context
+        :param GLProgram program: The program to set-up.
+                                  It MUST be in use and using this function.
+        """
+        if self.plane.isPlane:
+            gl.glUniform4f(program.uniforms['planeEq'], *self.plane.parameters)
+
+
+class DirectionalLight(event.Notifier, ProgramFunction):
+    """Description of a directional Phong light.
+
+    :param direction: The direction of the light or None to disable light
+    :type direction: ndarray of 3 floats or None
+    :param ambient: RGB ambient light
+    :type ambient: ndarray of 3 floats in [0, 1], default: (1., 1., 1.)
+    :param diffuse: RGB diffuse light parameter
+    :type diffuse: ndarray of 3 floats in [0, 1], default: (0., 0., 0.)
+    :param specular: RGB specular light parameter
+    :type specular: ndarray of 3 floats in [0, 1], default: (1., 1., 1.)
+    :param int shininess: The shininess of the material for specular term,
+                          default: 0 which disables specular component.
+    """
+
+    fragmentShaderFunction = """
+    /* Lighting */
+    struct DLight {
+        vec3 lightDir; // Direction of light in object space
+        vec3 ambient;
+        vec3 diffuse;
+        vec3 specular;
+        float shininess;
+        vec3 viewPos; // Camera position in object space
+    };
+
+    uniform DLight dLight;
+
+    vec4 lighting(vec4 color, vec3 position, vec3 normal)
+    {
+        normal = normalize(normal);
+        // 1-sided
+        float nDotL = max(0.0, dot(normal, - dLight.lightDir));
+
+        // 2-sided
+        //float nDotL = dot(normal, - dLight.lightDir);
+        //if (nDotL < 0.) {
+        //    nDotL = - nDotL;
+        //    normal = - normal;
+        //}
+
+        float specFactor = 0.;
+        if (dLight.shininess > 0. && nDotL > 0.) {
+            vec3 reflection = reflect(dLight.lightDir, normal);
+            vec3 viewDir = normalize(dLight.viewPos - position);
+            specFactor = max(0.0, dot(reflection, viewDir));
+            if (specFactor > 0.) {
+                specFactor = pow(specFactor, dLight.shininess);
+            }
+        }
+
+        vec3 enlightedColor = color.rgb * (dLight.ambient +
+                                           dLight.diffuse * nDotL) +
+                                           dLight.specular * specFactor;
+
+        return vec4(enlightedColor.rgb, color.a);
+    }
+    /* End of Lighting */
+    """
+
+    fragmentShaderFunctionNoop = """
+    vec4 lighting(vec4 color, vec3 position, vec3 normal)
+    {
+        return color;
+    }
+    """
+
+    def __init__(self, direction=None,
+                 ambient=(1., 1., 1.), diffuse=(0., 0., 0.),
+                 specular=(1., 1., 1.), shininess=0):
+        super(DirectionalLight, self).__init__()
+        self._direction = None
+        self.direction = direction  # Set _direction
+        self._isOn = True
+        self._ambient = ambient
+        self._diffuse = diffuse
+        self._specular = specular
+        self._shininess = shininess
+
+    ambient = event.notifyProperty('_ambient')
+    diffuse = event.notifyProperty('_diffuse')
+    specular = event.notifyProperty('_specular')
+    shininess = event.notifyProperty('_shininess')
+
+    @property
+    def isOn(self):
+        """True if light is on, False otherwise."""
+        return self._isOn and self._direction is not None
+
+    @isOn.setter
+    def isOn(self, isOn):
+        self._isOn = bool(isOn)
+
+    @contextlib.contextmanager
+    def turnOff(self):
+        """Context manager to temporary turn off lighting during rendering.
+
+        >>> with light.turnOff():
+        ...     # Do some rendering without lighting
+        """
+        wason = self._isOn
+        self._isOn = False
+        yield
+        self._isOn = wason
+
+    @property
+    def direction(self):
+        """The direction of the light, or None if light is not on."""
+        return self._direction
+
+    @direction.setter
+    def direction(self, direction):
+        if direction is None:
+            self._direction = None
+        else:
+            assert len(direction) == 3
+            direction = numpy.array(direction, dtype=numpy.float32, copy=True)
+            norm = numpy.linalg.norm(direction)
+            assert norm != 0
+            self._direction = direction / norm
+        self.notify()
+
+    # GL2
+
+    @property
+    def fragmentDef(self):
+        """Definition to add to fragment shader"""
+        if self.isOn:
+            return self.fragmentShaderFunction
+        else:
+            return self.fragmentShaderFunctionNoop
+
+    @property
+    def fragmentCall(self):
+        """Function name to call in fragment shader"""
+        return "lighting"
+
+    def setupProgram(self, context, program):
+        """Sets-up uniforms of a program using this shader function.
+
+        :param RenderContext context: The current rendering context
+        :param GLProgram program: The program to set-up.
+                                  It MUST be in use and using this function.
+        """
+        if self.isOn and self._direction is not None:
+            # Transform light direction from camera space to object coords
+            lightdir = context.objectToCamera.transformDir(
+                self._direction, direct=False)
+            lightdir /= numpy.linalg.norm(lightdir)
+
+            gl.glUniform3f(program.uniforms['dLight.lightDir'], *lightdir)
+
+            # Convert view position to object coords
+            viewpos = context.objectToCamera.transformPoint(
+                numpy.array((0., 0., 0., 1.), dtype=numpy.float32),
+                direct=False,
+                perspectiveDivide=True)[:3]
+            gl.glUniform3f(program.uniforms['dLight.viewPos'], *viewpos)
+
+            gl.glUniform3f(program.uniforms['dLight.ambient'], *self.ambient)
+            gl.glUniform3f(program.uniforms['dLight.diffuse'], *self.diffuse)
+            gl.glUniform3f(program.uniforms['dLight.specular'], *self.specular)
+            gl.glUniform1f(program.uniforms['dLight.shininess'],
+                           self.shininess)
+
+
+class Colormap(event.Notifier, ProgramFunction):
+    # TODO use colors for out-of-bound values, for <=0 with log, for nan
+    # TODO texture-based colormap
+
+    decl = """
+    #define CMAP_GRAY   0
+    #define CMAP_R_GRAY 1
+    #define CMAP_RED    2
+    #define CMAP_GREEN  3
+    #define CMAP_BLUE   4
+    #define CMAP_TEMP   5
+
+    uniform struct {
+        int id;
+        bool isLog;
+        float min;
+        float oneOverRange;
+    } cmap;
+
+    const float oneOverLog10 = 0.43429448190325176;
+
+    vec4 colormap(float value) {
+        if (cmap.isLog) { /* Log10 mapping */
+            if (value > 0.0) {
+                value = clamp(cmap.oneOverRange *
+                              (oneOverLog10 * log(value) - cmap.min),
+                              0.0, 1.0);
+            } else {
+                value = 0.0;
+            }
+        } else { /* Linear mapping */
+            value = clamp(cmap.oneOverRange * (value - cmap.min), 0.0, 1.0);
+        }
+
+        if (cmap.id == CMAP_GRAY) {
+            return vec4(value, value, value, 1.0);
+        }
+        else if (cmap.id == CMAP_R_GRAY) {
+            float invValue = 1.0 - value;
+            return vec4(invValue, invValue, invValue, 1.0);
+        }
+        else if (cmap.id == CMAP_RED) {
+            return vec4(value, 0.0, 0.0, 1.0);
+        }
+        else if (cmap.id == CMAP_GREEN) {
+            return vec4(0.0, value, 0.0, 1.0);
+        }
+        else if (cmap.id == CMAP_BLUE) {
+            return vec4(0.0, 0.0, value, 1.0);
+        }
+        else if (cmap.id == CMAP_TEMP) {
+            //red: 0.5->0.75: 0->1
+            //green: 0.->0.25: 0->1; 0.75->1.: 1->0
+            //blue: 0.25->0.5: 1->0
+            return vec4(
+                clamp(4.0 * value - 2.0, 0.0, 1.0),
+                1.0 - clamp(4.0 * abs(value - 0.5) - 1.0, 0.0, 1.0),
+                1.0 - clamp(4.0 * value - 1.0, 0.0, 1.0),
+                1.0);
+        }
+        else {
+            /* Unknown colormap */
+            return vec4(0.0, 0.0, 0.0, 1.0);
+        }
+    }
+    """
+
+    call = "colormap"
+
+    _COLORMAPS = {
+        'gray': 0,
+        'reversed gray': 1,
+        'red': 2,
+        'green': 3,
+        'blue': 4,
+        'temperature': 5
+    }
+
+    COLORMAPS = tuple(_COLORMAPS.keys())
+    """Tuple of supported colormap names."""
+
+    NORMS = 'linear', 'log'
+    """Tuple of supported normalizations."""
+
+    def __init__(self, name='gray', norm='linear', range_=(1., 10.)):
+        """Shader function to apply a colormap to a value.
+
+        :param str name: Name of the colormap.
+        :param str norm: Normalization to apply: 'linear' (default) or 'log'.
+        :param range_: Range of value to map to the colormap.
+        :type range_: 2-tuple of float (begin, end).
+        """
+        super(Colormap, self).__init__()
+
+        # Init privates to default
+        self._name, self._norm, self._range = 'gray', 'linear', (1., 10.)
+
+        # Set to param values through properties to go through asserts
+        self.name = name
+        self.norm = norm
+        self.range_ = range_
+
+    @property
+    def name(self):
+        """Name of the colormap in use."""
+        return self._name
+
+    @name.setter
+    def name(self, name):
+        if name != self._name:
+            assert name in self.COLORMAPS
+            self._name = name
+            self.notify()
+
+    @property
+    def norm(self):
+        """Normalization to use for colormap mapping.
+
+        Either 'linear' (the default) or 'log' for log10 mapping.
+        With 'log' normalization, values <= 0. are set to 1. (i.e. log == 0)
+        """
+        return self._norm
+
+    @norm.setter
+    def norm(self, norm):
+        if norm != self._norm:
+            assert norm in self.NORMS
+            self._norm = norm
+            if norm == 'log':
+                self.range_ = self.range_  # To test for positive range_
+            self.notify()
+
+    @property
+    def range_(self):
+        """Range of values to map to the colormap.
+
+        2-tuple of floats: (begin, end).
+        The begin value is mapped to the origin of the colormap and the
+        end value is mapped to the other end of the colormap.
+        The colormap is reversed if begin > end.
+        """
+        return self._range
+
+    @range_.setter
+    def range_(self, range_):
+        assert len(range_) == 2
+        range_ = float(range_[0]), float(range_[1])
+
+        if self.norm == 'log' and (range_[0] <= 0. or range_[1] <= 0.):
+            _logger.warn(
+                "Log normalization and negative range: updating range.")
+            minPos = numpy.finfo(numpy.float32).tiny
+            range_ = max(range_[0], minPos), max(range_[1], minPos)
+
+        if range_ != self._range:
+            self._range = range_
+            self.notify()
+
+    def setupProgram(self, context, program):
+        """Sets-up uniforms of a program using this shader function.
+
+        :param RenderContext context: The current rendering context
+        :param GLProgram program: The program to set-up.
+                                  It MUST be in use and using this function.
+        """
+        gl.glUniform1i(program.uniforms['cmap.id'], self._COLORMAPS[self.name])
+        gl.glUniform1i(program.uniforms['cmap.isLog'], self._norm == 'log')
+
+        min_, max_ = self.range_
+        if self._norm == 'log':
+            min_, max_ = numpy.log10(min_), numpy.log10(max_)
+
+        gl.glUniform1f(program.uniforms['cmap.min'], min_)
+        gl.glUniform1f(program.uniforms['cmap.oneOverRange'],
+                       (1. / (max_ - min_)) if max_ != min_ else 0.)
diff --git a/silx/gui/plot3d/scene/interaction.py b/silx/gui/plot3d/scene/interaction.py
new file mode 100644
index 0000000..68bfc13
--- /dev/null
+++ b/silx/gui/plot3d/scene/interaction.py
@@ -0,0 +1,652 @@
+# coding: utf-8
+# /*##########################################################################
+#
+# Copyright (c) 2015-2017 European Synchrotron Radiation Facility
+#
+# Permission is hereby granted, free of charge, to any person obtaining a copy
+# of this software and associated documentation files (the "Software"), to deal
+# in the Software without restriction, including without limitation the rights
+# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+# copies of the Software, and to permit persons to whom the Software is
+# furnished to do so, subject to the following conditions:
+#
+# The above copyright notice and this permission notice shall be included in
+# all copies or substantial portions of the Software.
+#
+# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
+# THE SOFTWARE.
+#
+# ###########################################################################*/
+"""This module provides interaction to plug on the scene graph."""
+
+from __future__ import absolute_import
+
+__authors__ = ["T. Vincent"]
+__license__ = "MIT"
+__date__ = "25/07/2016"
+
+import logging
+import numpy
+
+from silx.gui.plot.Interaction import \
+    StateMachine, State, LEFT_BTN, RIGHT_BTN  # , MIDDLE_BTN
+
+from . import transform
+
+
+_logger = logging.getLogger(__name__)
+
+
+# ClickOrDrag #################################################################
+
+# TODO merge with silx.gui.plot.Interaction.ClickOrDrag
+class ClickOrDrag(StateMachine):
+    """Click or drag interaction for a given button."""
+
+    DRAG_THRESHOLD_SQUARE_DIST = 5 ** 2
+
+    class Idle(State):
+        def onPress(self, x, y, btn):
+            if btn == self.machine.button:
+                self.goto('clickOrDrag', x, y)
+                return True
+
+    class ClickOrDrag(State):
+        def enterState(self, x, y):
+            self.initPos = x, y
+
+        enter = enterState  # silx v.0.3 support, remove when 0.4 out
+
+        def onMove(self, x, y):
+            dx = (x - self.initPos[0]) ** 2
+            dy = (y - self.initPos[1]) ** 2
+            if (dx ** 2 + dy ** 2) >= self.machine.DRAG_THRESHOLD_SQUARE_DIST:
+                self.goto('drag', self.initPos, (x, y))
+
+        def onRelease(self, x, y, btn):
+            if btn == self.machine.button:
+                self.machine.click(x, y)
+                self.goto('idle')
+
+    class Drag(State):
+        def enterState(self, initPos, curPos):
+            self.initPos = initPos
+            self.machine.beginDrag(*initPos)
+            self.machine.drag(*curPos)
+
+        enter = enterState  # silx v.0.3 support, remove when 0.4 out
+
+        def onMove(self, x, y):
+            self.machine.drag(x, y)
+
+        def onRelease(self, x, y, btn):
+            if btn == self.machine.button:
+                self.machine.endDrag(self.initPos, (x, y))
+                self.goto('idle')
+
+    def __init__(self, button=LEFT_BTN):
+        self.button = button
+        states = {
+            'idle': ClickOrDrag.Idle,
+            'clickOrDrag': ClickOrDrag.ClickOrDrag,
+            'drag': ClickOrDrag.Drag
+        }
+        super(ClickOrDrag, self).__init__(states, 'idle')
+
+    def click(self, x, y):
+        """Called upon a left or right button click.
+        To override in a subclass.
+        """
+        pass
+
+    def beginDrag(self, x, y):
+        """Called at the beginning of a drag gesture with left button
+        pressed.
+        To override in a subclass.
+        """
+        pass
+
+    def drag(self, x, y):
+        """Called on mouse moved during a drag gesture.
+        To override in a subclass.
+        """
+        pass
+
+    def endDrag(self, x, y):
+        """Called at the end of a drag gesture when the left button is
+        released.
+        To override in a subclass.
+        """
+        pass
+
+
+# CameraRotate ################################################################
+
+class CameraRotate(ClickOrDrag):
+    """Camera rotation using an arcball-like interaction."""
+
+    def __init__(self, viewport, orbitAroundCenter=True, button=RIGHT_BTN):
+        self._viewport = viewport
+        self._orbitAroundCenter = orbitAroundCenter
+        self._reset()
+        super(CameraRotate, self).__init__(button)
+
+    def _reset(self):
+        self._origin, self._center = None, None
+        self._startExtrinsic = None
+
+    def click(self, x, y):
+        pass  # No interaction yet
+
+    def beginDrag(self, x, y):
+        centerPos = None
+        if not self._orbitAroundCenter:
+            # Try to use picked object position as center of rotation
+            ndcZ = self._viewport._pickNdcZGL(x, y)
+            if ndcZ != 1.:
+                # Hit an object, use picked point as center
+                centerPos = self._viewport._getXZYGL(x, y)  # Can return None
+
+        if centerPos is None:
+            # Not using picked position, use scene center
+            bounds = self._viewport.scene.bounds(transformed=True)
+            centerPos = 0.5 * (bounds[0] + bounds[1])
+
+        self._center = transform.Translate(*centerPos)
+        self._origin = x, y
+        self._startExtrinsic = self._viewport.camera.extrinsic.copy()
+
+    def drag(self, x, y):
+        if self._center is None:
+            return
+
+        dx, dy = self._origin[0] - x, self._origin[1] - y
+
+        if dx == 0 and dy == 0:
+            direction = self._startExtrinsic.direction
+            up = self._startExtrinsic.up
+            position = self._startExtrinsic.position
+        else:
+            minsize = min(self._viewport.size)
+            distance = numpy.sqrt(dx ** 2 + dy ** 2)
+            angle = distance / minsize * numpy.pi
+
+            # Take care of y inversion
+            direction = dx * self._startExtrinsic.side - \
+                dy * self._startExtrinsic.up
+            direction /= numpy.linalg.norm(direction)
+            axis = numpy.cross(direction, self._startExtrinsic.direction)
+            axis /= numpy.linalg.norm(axis)
+
+            # Orbit start camera with current angle and axis
+            # Rotate viewing direction
+            rotation = transform.Rotate(numpy.degrees(angle), *axis)
+            direction = rotation.transformDir(self._startExtrinsic.direction)
+            up = rotation.transformDir(self._startExtrinsic.up)
+
+            # Rotate position around center
+            trlist = transform.StaticTransformList((
+                self._center,
+                rotation,
+                self._center.inverse()))
+            position = trlist.transformPoint(self._startExtrinsic.position)
+
+        camerapos = self._viewport.camera.extrinsic
+        camerapos.setOrientation(direction, up)
+        camerapos.position = position
+
+    def endDrag(self, x, y):
+        self._reset()
+
+
+# CameraSelectPan #############################################################
+
+class CameraSelectPan(ClickOrDrag):
+    """Picking on click and pan camera on drag."""
+
+    def __init__(self, viewport, button=LEFT_BTN, selectCB=None):
+        self._viewport = viewport
+        self._selectCB = selectCB
+        self._lastPosNdc = None
+        super(CameraSelectPan, self).__init__(button)
+
+    def click(self, x, y):
+        if self._selectCB is not None:
+            ndcZ = self._viewport._pickNdcZGL(x, y)
+            position = self._viewport._getXZYGL(x, y)
+            # This assume no object lie on the far plane
+            # Alternative, change the depth range so that far is < 1
+            if ndcZ != 1. and position is not None:
+                self._selectCB((x, y, ndcZ), position)
+
+    def beginDrag(self, x, y):
+        ndc = self._viewport.windowToNdc(x, y)
+        ndcZ = self._viewport._pickNdcZGL(x, y)
+        # ndcZ is the panning plane
+        if ndc is not None and ndcZ is not None:
+            self._lastPosNdc = numpy.array((ndc[0], ndc[1], ndcZ, 1.),
+                                           dtype=numpy.float32)
+        else:
+            self._lastPosNdc = None
+
+    def drag(self, x, y):
+        if self._lastPosNdc is not None:
+            ndc = self._viewport.windowToNdc(x, y)
+            if ndc is not None:
+                ndcPos = numpy.array((ndc[0], ndc[1], self._lastPosNdc[2], 1.),
+                                     dtype=numpy.float32)
+
+                # Convert last and current NDC positions to scene coords
+                scenePos = self._viewport.camera.transformPoint(
+                    ndcPos, direct=False, perspectiveDivide=True)
+                lastScenePos = self._viewport.camera.transformPoint(
+                    self._lastPosNdc, direct=False, perspectiveDivide=True)
+
+                # Get translation in scene coords
+                translation = scenePos[:3] - lastScenePos[:3]
+                self._viewport.camera.extrinsic.position -= translation
+
+                # Store for next drag
+                self._lastPosNdc = ndcPos
+
+    def endDrag(self, x, y):
+        self._lastPosNdc = None
+
+
+# CameraWheel #################################################################
+
+class CameraWheel(object):
+    """StateMachine like class, just handling wheel events."""
+
+    # TODO choose scale of motion? Translation or Scale?
+    def __init__(self, viewport, mode='center', scaleTransform=None):
+        assert mode in ('center', 'position', 'scale')
+        self._viewport = viewport
+        if mode == 'center':
+            self._zoomTo = self._zoomToCenter
+        elif mode == 'position':
+            self._zoomTo = self._zoomToPosition
+        elif mode == 'scale':
+            self._zoomTo = self._zoomByScale
+            self._scale = scaleTransform
+        else:
+            raise ValueError('Unsupported mode: %s' % mode)
+
+    def handleEvent(self, eventName, *args, **kwargs):
+        if eventName == 'wheel':
+            return self._zoomTo(*args, **kwargs)
+
+    def _zoomToCenter(self, x, y, angleInDegrees):
+        """Zoom to center of display.
+
+        Only works with perspective camera.
+        """
+        direction = 'forward' if angleInDegrees > 0 else 'backward'
+        self._viewport.camera.move(direction)
+        return True
+
+    def _zoomToPositionAbsolute(self, x, y, angleInDegrees):
+        """Zoom while keeping pixel under mouse invariant.
+
+        Only works with perspective camera.
+        """
+        ndc = self._viewport.windowToNdc(x, y)
+        if ndc is not None:
+            near = numpy.array((ndc[0], ndc[1], -1., 1.), dtype=numpy.float32)
+
+            nearscene = self._viewport.camera.transformPoint(
+                near, direct=False, perspectiveDivide=True)
+
+            far = numpy.array((ndc[0], ndc[1], 1., 1.), dtype=numpy.float32)
+            farscene = self._viewport.camera.transformPoint(
+                far, direct=False, perspectiveDivide=True)
+
+            dirscene = farscene[:3] - nearscene[:3]
+            dirscene /= numpy.linalg.norm(dirscene)
+
+            if angleInDegrees < 0:
+                dirscene *= -1.
+
+            # TODO which scale
+            self._viewport.camera.extrinsic.position += dirscene
+        return True
+
+    def _zoomToPosition(self, x, y, angleInDegrees):
+        """Zoom while keeping pixel under mouse invariant."""
+        projection = self._viewport.camera.intrinsic
+        extrinsic = self._viewport.camera.extrinsic
+
+        if isinstance(projection, transform.Perspective):
+            # For perspective projection, move camera
+            ndc = self._viewport.windowToNdc(x, y)
+            if ndc is not None:
+                ndcz = self._viewport._pickNdcZGL(x, y)
+
+                position = numpy.array((ndc[0], ndc[1], ndcz),
+                                       dtype=numpy.float32)
+                positionscene = self._viewport.camera.transformPoint(
+                    position, direct=False, perspectiveDivide=True)
+
+                camtopos = extrinsic.position - positionscene
+
+                step = 0.2 * (1. if angleInDegrees < 0 else -1.)
+                extrinsic.position += step * camtopos
+
+        elif isinstance(projection, transform.Orthographic):
+            # For orthographic projection, change projection borders
+            ndcx, ndcy = self._viewport.windowToNdc(x, y, checkInside=False)
+
+            step = 0.2 * (1. if angleInDegrees < 0 else -1.)
+
+            dx = (ndcx + 1) / 2.
+            stepwidth = step * (projection.right - projection.left)
+            left = projection.left - dx * stepwidth
+            right = projection.right + (1. - dx) * stepwidth
+
+            dy = (ndcy + 1) / 2.
+            stepheight = step * (projection.top - projection.bottom)
+            bottom = projection.bottom - dy * stepheight
+            top = projection.top + (1. - dy) * stepheight
+
+            projection.setClipping(left, right, bottom, top)
+
+        else:
+            raise RuntimeError('Unsupported camera', projection)
+        return True
+
+    def _zoomByScale(self, x, y, angleInDegrees):
+        """Zoom by scaling scene (do not keep pixel under mouse invariant)."""
+        scalefactor = 1.1
+        if angleInDegrees < 0.:
+            scalefactor = 1. / scalefactor
+        self._scale.scale = scalefactor * self._scale.scale
+
+        self._viewport.adjustCameraDepthExtent()
+        return True
+
+
+# FocusManager ################################################################
+
+class FocusManager(StateMachine):
+    """Manages focus across multiple event handlers
+
+    On press an event handler can acquire focus.
+    By default it looses focus when all buttons are released.
+    """
+    class Idle(State):
+        def onPress(self, x, y, btn):
+            for eventHandler in self.machine.eventHandlers:
+                requestfocus = eventHandler.handleEvent('press', x, y, btn)
+                if requestfocus:
+                    self.goto('focus', eventHandler, btn)
+                    break
+
+        def _processEvent(self, *args):
+            for eventHandler in self.machine.eventHandlers:
+                consumeevent = eventHandler.handleEvent(*args)
+                if consumeevent:
+                    break
+
+        def onMove(self, x, y):
+            self._processEvent('move', x, y)
+
+        def onRelease(self, x, y, btn):
+            self._processEvent('release', x, y, btn)
+
+        def onWheel(self, x, y, angle):
+            self._processEvent('wheel', x, y, angle)
+
+    class Focus(State):
+        def enterState(self, eventHandler, btn):
+            self.eventHandler = eventHandler
+            self.focusBtns = {btn}  # Set
+
+        enter = enterState  # silx v.0.3 support, remove when 0.4 out
+
+        def onPress(self, x, y, btn):
+            self.focusBtns.add(btn)
+            self.eventHandler.handleEvent('press', x, y, btn)
+
+        def onMove(self, x, y):
+            self.eventHandler.handleEvent('move', x, y)
+
+        def onRelease(self, x, y, btn):
+            self.focusBtns.discard(btn)
+            requestfocus = self.eventHandler.handleEvent('release', x, y, btn)
+            if len(self.focusBtns) == 0 and not requestfocus:
+                self.goto('idle')
+
+        def onWheel(self, x, y, angleInDegrees):
+            self.eventHandler.handleEvent('wheel', x, y, angleInDegrees)
+
+    def __init__(self, eventHandlers=()):
+        self.eventHandlers = list(eventHandlers)
+
+        states = {
+            'idle': FocusManager.Idle,
+            'focus': FocusManager.Focus
+        }
+        super(FocusManager, self).__init__(states, 'idle')
+
+    def cancel(self):
+        for handler in self.eventHandlers:
+            handler.cancel()
+
+
+# CameraControl ###############################################################
+
+class CameraControl(FocusManager):
+    """Combine wheel, selectPan and rotate state machine."""
+    def __init__(self, viewport,
+                 orbitAroundCenter=False,
+                 mode='center', scaleTransform=None,
+                 selectCB=None):
+        handlers = (CameraWheel(viewport, mode, scaleTransform),
+                    CameraSelectPan(viewport, LEFT_BTN, selectCB),
+                    CameraRotate(viewport, orbitAroundCenter, RIGHT_BTN))
+        super(CameraControl, self).__init__(handlers)
+
+
+# PlaneRotate #################################################################
+
+class PlaneRotate(ClickOrDrag):
+    """Plane rotation using arcball interaction.
+
+    Arcball ref.:
+    Ken Shoemake. ARCBALL: A user interface for specifying three-dimensional
+    orientation using a mouse. In Proc. GI '92. (1992). pp. 151-156.
+    """
+
+    def __init__(self, viewport, plane, button=RIGHT_BTN):
+        self._viewport = viewport
+        self._plane = plane
+        self._reset()
+        super(PlaneRotate, self).__init__(button)
+
+    def _reset(self):
+        self._beginNormal, self._beginCenter = None, None
+
+    def click(self, x, y):
+        pass  # No interaction
+
+    @staticmethod
+    def _sphereUnitVector(radius, center, position):
+        """Returns the unit vector of the projection of position on a sphere.
+
+        It assumes an orthographic projection.
+        For perspective projection, it gives an approximation, but it
+        simplifies computations and results in consistent arcball control
+        in control space.
+
+        All parameters must be in screen coordinate system
+        (either pixels or normalized coordinates).
+
+        :param float radius: The radius of the sphere.
+        :param center: (x, y) coordinates of the center.
+        :param position: (x, y) coordinates of the cursor position.
+        :return: Unit vector.
+        :rtype: numpy.ndarray of 3 floats.
+        """
+        center, position = numpy.array(center), numpy.array(position)
+
+        # Normalize x and y on a unit circle
+        spherecoords = (position - center) / float(radius)
+        squarelength = numpy.sum(spherecoords ** 2)
+
+        # Project on the unit sphere and compute z coordinates
+        if squarelength > 1.0:  # Outside sphere: project
+            spherecoords /= numpy.sqrt(squarelength)
+            zsphere = 0.0
+        else:  # In sphere: compute z
+            zsphere = numpy.sqrt(1. - squarelength)
+
+        spherecoords = numpy.append(spherecoords, zsphere)
+        return spherecoords
+
+    def beginDrag(self, x, y):
+        # Makes sure the point defining the plane is at the center as
+        # it will be the center of rotation (as rotation is applied to normal)
+        self._plane.plane.point = self._plane.center
+
+        # Store the plane normal
+        self._beginNormal = self._plane.plane.normal
+
+        _logger.debug(
+            'Begin arcball, plane center %s', str(self._plane.center))
+
+        # Do the arcball on the screen
+        radius = min(self._viewport.size)
+        if self._plane.center is None:
+            self._beginCenter = None
+
+        else:
+            center = self._plane.objectToNDCTransform.transformPoint(
+                self._plane.center, perspectiveDivide=True)
+            self._beginCenter = self._viewport.ndcToWindow(
+                center[0], center[1], checkInside=False)
+
+            self._startVector = self._sphereUnitVector(
+                radius, self._beginCenter, (x, y))
+
+    def drag(self, x, y):
+        if self._beginCenter is None:
+            return
+
+        # Compute rotation: this is twice the rotation of the arcball
+        radius = min(self._viewport.size)
+        currentvector = self._sphereUnitVector(
+            radius, self._beginCenter, (x, y))
+        crossprod = numpy.cross(self._startVector, currentvector)
+        dotprod = numpy.dot(self._startVector, currentvector)
+
+        quaternion = numpy.append(crossprod, dotprod)
+        # Rotation was computed with Y downward, but apply in NDC, invert Y
+        quaternion[1] *= -1.
+
+        rotation = transform.Rotate()
+        rotation.quaternion = quaternion
+
+        # Convert to NDC, rotate, convert back to object
+        normal = self._plane.objectToNDCTransform.transformNormal(
+            self._beginNormal)
+        normal = rotation.transformNormal(normal)
+        normal = self._plane.objectToNDCTransform.transformNormal(
+            normal, direct=False)
+        self._plane.plane.normal = normal
+
+    def endDrag(self, x, y):
+        self._reset()
+
+
+# PlanePan ###################################################################
+
+class PlanePan(ClickOrDrag):
+    """Pan a plane along its normal on drag."""
+
+    def __init__(self, viewport, plane, button=LEFT_BTN):
+        self._plane = plane
+        self._viewport = viewport
+        self._beginPlanePoint = None
+        self._beginPos = None
+        self._dragNdcZ = 0.
+        super(PlanePan, self).__init__(button)
+
+    def click(self, x, y):
+        pass
+
+    def beginDrag(self, x, y):
+        ndc = self._viewport.windowToNdc(x, y)
+        ndcZ = self._viewport._pickNdcZGL(x, y)
+        # ndcZ is the panning plane
+        if ndc is not None and ndcZ is not None:
+            ndcPos = numpy.array((ndc[0], ndc[1], ndcZ, 1.),
+                                 dtype=numpy.float32)
+            scenePos = self._viewport.camera.transformPoint(
+                ndcPos, direct=False, perspectiveDivide=True)
+            self._beginPos = self._plane.objectToSceneTransform.transformPoint(
+                scenePos, direct=False)
+            self._dragNdcZ = ndcZ
+        else:
+            self._beginPos = None
+            self._dragNdcZ = 0.
+
+        self._beginPlanePoint = self._plane.plane.point
+
+    def drag(self, x, y):
+        if self._beginPos is not None:
+            ndc = self._viewport.windowToNdc(x, y)
+            if ndc is not None:
+                ndcPos = numpy.array((ndc[0], ndc[1], self._dragNdcZ, 1.),
+                                     dtype=numpy.float32)
+
+                # Convert last and current NDC positions to scene coords
+                scenePos = self._viewport.camera.transformPoint(
+                    ndcPos, direct=False, perspectiveDivide=True)
+                curPos = self._plane.objectToSceneTransform.transformPoint(
+                    scenePos, direct=False)
+
+                # Get translation in scene coords
+                translation = curPos[:3] - self._beginPos[:3]
+
+                newPoint = self._beginPlanePoint + translation
+
+                # Keep plane point in bounds
+                bounds = self._plane.parent.bounds(dataBounds=True)
+                if bounds is not None:
+                    newPoint = numpy.clip(
+                        newPoint, a_min=bounds[0], a_max=bounds[1])
+
+                    # Only update plane if it is in some bounds
+                    self._plane.plane.point = newPoint
+
+    def endDrag(self, x, y):
+        self._beginPlanePoint = None
+
+
+# PlaneControl ################################################################
+
+class PlaneControl(FocusManager):
+    """Combine wheel, selectPan and rotate state machine for plane control."""
+    def __init__(self, viewport, plane,
+                 mode='center', scaleTransform=None):
+        handlers = (CameraWheel(viewport, mode, scaleTransform),
+                    PlanePan(viewport, plane, LEFT_BTN),
+                    PlaneRotate(viewport, plane, RIGHT_BTN))
+        super(PlaneControl, self).__init__(handlers)
+
+
+class PanPlaneRotateCameraControl(FocusManager):
+    """Combine wheel, pan plane and camera rotate state machine."""
+    def __init__(self, viewport, plane,
+                 mode='center', scaleTransform=None):
+        handlers = (CameraWheel(viewport, mode, scaleTransform),
+                    PlanePan(viewport, plane, LEFT_BTN),
+                    CameraRotate(viewport,
+                                 orbitAroundCenter=False,
+                                 button=RIGHT_BTN))
+        super(PanPlaneRotateCameraControl, self).__init__(handlers)
diff --git a/silx/gui/plot3d/scene/primitives.py b/silx/gui/plot3d/scene/primitives.py
new file mode 100644
index 0000000..ca2616a
--- /dev/null
+++ b/silx/gui/plot3d/scene/primitives.py
@@ -0,0 +1,1764 @@
+# coding: utf-8
+# /*##########################################################################
+#
+# Copyright (c) 2015-2017 European Synchrotron Radiation Facility
+#
+# Permission is hereby granted, free of charge, to any person obtaining a copy
+# of this software and associated documentation files (the "Software"), to deal
+# in the Software without restriction, including without limitation the rights
+# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+# copies of the Software, and to permit persons to whom the Software is
+# furnished to do so, subject to the following conditions:
+#
+# The above copyright notice and this permission notice shall be included in
+# all copies or substantial portions of the Software.
+#
+# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
+# THE SOFTWARE.
+#
+# ###########################################################################*/
+
+from __future__ import absolute_import, division, unicode_literals
+
+__authors__ = ["T. Vincent"]
+__license__ = "MIT"
+__date__ = "25/07/2016"
+
+
+import collections
+import ctypes
+from functools import reduce
+import logging
+import string
+
+import numpy
+
+from silx.gui.plot.Colors import rgba
+
+from ... import _glutils
+from ..._glutils import gl
+
+from . import event
+from . import core
+from . import transform
+from . import utils
+
+_logger = logging.getLogger(__name__)
+
+
+# Geometry ####################################################################
+
+class Geometry(core.Elem):
+    """Set of vertices with normals and colors.
+
+    :param str mode: OpenGL drawing mode:
+                     lines, line_strip, loop, triangles, triangle_strip, fan
+    :param indices: Array of vertex indices or None
+    :param bool copy: True (default) to copy the data, False to use as is.
+    :param attributes: Provide list of attributes as extra parameters.
+    """
+
+    _ATTR_INFO = {
+        'position': {'dims': (1, 2), 'lastDim': (2, 3, 4)},
+        'normal': {'dims': (1, 2), 'lastDim': (3,)},
+        'color': {'dims': (1, 2), 'lastDim': (3, 4)},
+    }
+
+    _MODE_CHECKS = {  # Min, Modulo
+        'lines': (2, 2), 'line_strip': (2, 0), 'loop': (2, 0),
+        'points': (1, 0),
+        'triangles': (3, 3), 'triangle_strip': (3, 0), 'fan': (3, 0)
+    }
+
+    _MODES = {
+        'lines': gl.GL_LINES,
+        'line_strip': gl.GL_LINE_STRIP,
+        'loop': gl.GL_LINE_LOOP,
+
+        'points': gl.GL_POINTS,
+
+        'triangles': gl.GL_TRIANGLES,
+        'triangle_strip': gl.GL_TRIANGLE_STRIP,
+        'fan': gl.GL_TRIANGLE_FAN
+    }
+
+    _LINE_MODES = 'lines', 'line_strip', 'loop'
+
+    _TRIANGLE_MODES = 'triangles', 'triangle_strip', 'fan'
+
+    def __init__(self, mode, indices=None, copy=True, **attributes):
+        super(Geometry, self).__init__()
+
+        self._vbos = {}  # Store current vbos
+        self._unsyncAttributes = []  # Store attributes to copy to vbos
+        self.__bounds = None  # Cache object's bounds
+
+        assert mode in self._MODES
+        self._mode = mode
+
+        # Set attributes
+        self._attributes = {}
+        for name, data in attributes.items():
+            self.setAttribute(name, data, copy=copy)
+
+        # Set indices
+        self._indices = None
+        self.setIndices(indices, copy=copy)
+
+        # More consistency checks
+        mincheck, modulocheck = self._MODE_CHECKS[self._mode]
+        if self._indices is not None:
+            nbvertices = len(self._indices)
+        else:
+            nbvertices = self.nbVertices
+        assert nbvertices >= mincheck
+        if modulocheck != 0:
+            assert (nbvertices % modulocheck) == 0
+
+    @staticmethod
+    def _glReadyArray(array, copy=True):
+        """Making a contiguous array, checking float types.
+
+        :param iterable array: array-like data to prepare for attribute
+        :param bool copy: True to make a copy of the array, False to use as is
+        """
+        # Convert single value (int, float, numpy types) to tuple
+        if not isinstance(array, collections.Iterable):
+            array = (array, )
+
+        # Makes sure it is an array
+        array = numpy.array(array, copy=False)
+
+        # Cast all float to float32
+        dtype = None
+        if numpy.dtype(array.dtype).kind == 'f':
+            dtype = numpy.float32
+
+        return numpy.array(array, dtype=dtype, order='C', copy=copy)
+
+    @property
+    def nbVertices(self):
+        """Returns the number of vertices of current attributes.
+
+        It returns None if there is no attributes.
+        """
+        for array in self._attributes.values():
+            if len(array.shape) == 2:
+                return len(array)
+        return None
+
+    def setAttribute(self, name, array, copy=True):
+        """Set attribute with provided array.
+
+        :param str name: The name of the attribute
+        :param array: Array-like attribute data or None to remove attribute
+        :param bool copy: True (default) to copy the data, False to use as is
+        """
+        # This triggers associated GL resources to be garbage collected
+        self._vbos.pop(name, None)
+
+        if array is None:
+            self._attributes.pop(name, None)
+
+        else:
+            array = self._glReadyArray(array, copy=copy)
+
+            if name not in self._ATTR_INFO:
+                _logger.info('Not checking attibute %s dimensions', name)
+            else:
+                checks = self._ATTR_INFO[name]
+
+                if (len(array.shape) == 1 and checks['lastDim'] == (1,) and
+                        len(array) > 1):
+                    array = array.reshape((len(array), 1))
+
+                # Checks
+                assert len(array.shape) in checks['dims'], "Attr %s" % name
+                assert array.shape[-1] in checks['lastDim'], "Attr %s" % name
+
+            # Check length against another attribute array
+            # Causes problems when updating
+            # nbVertices = self.nbVertices
+            # if len(array.shape) == 2 and nbVertices is not None:
+            #     assert len(array) == nbVertices
+
+            self._attributes[name] = array
+            if len(array.shape) == 2:  # Store this in a VBO
+                self._unsyncAttributes.append(name)
+
+            if name == 'position':  # Reset bounds
+                self.__bounds = None
+
+        self.notify()
+
+    def getAttribute(self, name, copy=True):
+        """Returns the numpy.ndarray corresponding to the name attribute.
+
+        :param str name: The name of the attribute to get.
+        :param bool copy: True to get a copy (default),
+                          False to get internal array (DO NOT MODIFY)
+        :return: The corresponding array or None if no corresponding attribute.
+        :rtype: numpy.ndarray
+        """
+        attr = self._attributes.get(name, None)
+        return None if attr is None else numpy.array(attr, copy=copy)
+
+    def useAttribute(self, program, name=None):
+        """Enable and bind attribute(s) for a specific program.
+
+        This MUST be called with OpenGL context active and after prepareGL2
+        has been called.
+
+        :param GLProgram program: The program for which to set the attributes
+        :param str name: The attribute name to set or None to set then all
+        """
+        if name is None:
+            for name in program.attributes:
+                self.useAttribute(program, name)
+
+        else:
+            attribute = program.attributes.get(name)
+            if attribute is None:
+                return
+
+            vboattrib = self._vbos.get(name)
+            if vboattrib is not None:
+                gl.glEnableVertexAttribArray(attribute)
+                vboattrib.setVertexAttrib(attribute)
+
+            elif name not in self._attributes:
+                gl.glDisableVertexAttribArray(attribute)
+
+            else:
+                array = self._attributes[name]
+                assert array is not None
+
+                if len(array.shape) == 1:
+                    assert len(array) in (1, 2, 3, 4)
+                    gl.glDisableVertexAttribArray(attribute)
+                    _glVertexAttribFunc = getattr(
+                        _glutils.gl, 'glVertexAttrib{}f'.format(len(array)))
+                    _glVertexAttribFunc(attribute, *array)
+                else:
+                    # TODO As is this is a never event, remove?
+                    gl.glEnableVertexAttribArray(attribute)
+                    gl.glVertexAttribPointer(
+                        attribute,
+                        array.shape[-1],
+                        _glutils.numpyToGLType(array.dtype),
+                        gl.GL_FALSE,
+                        0,
+                        array)
+
+    def setIndices(self, indices, copy=True):
+        """Set the primitive indices to use.
+
+        :param indices: Array-like of uint primitive indices or None to unset
+        :param bool copy: True (default) to copy the data, False to use as is
+        """
+        # Trigger garbage collection of previous indices VBO if any
+        self._vbos.pop('__indices__', None)
+
+        if indices is None:
+            self._indices = None
+        else:
+            indices = self._glReadyArray(indices, copy=copy).ravel()
+            assert indices.dtype.name in ('uint8', 'uint16', 'uint32')
+            if _logger.getEffectiveLevel() <= logging.DEBUG:
+                # This might be a costy check
+                assert indices.max() < self.nbVertices
+            self._indices = indices
+
+    def getIndices(self, copy=True):
+        """Returns the numpy.ndarray corresponding to the indices.
+
+        :param bool copy: True to get a copy (default),
+                          False to get internal array (DO NOT MODIFY)
+        :return: The primitive indices array or None if not set.
+        :rtype: numpy.ndarray or None
+        """
+        if self._indices is None:
+            return None
+        else:
+            return numpy.array(self._indices, copy=copy)
+
+    def _bounds(self, dataBounds=False):
+        if self.__bounds is None:
+            self.__bounds = numpy.zeros((2, 3), dtype=numpy.float32)
+            # Support vertex with to 2 to 4 coordinates
+            positions = self._attributes['position']
+            self.__bounds[0, :positions.shape[1]] = positions.min(axis=0)[:3]
+            self.__bounds[1, :positions.shape[1]] = positions.max(axis=0)[:3]
+        return self.__bounds.copy()
+
+    def prepareGL2(self, ctx):
+        # TODO manage _vbo and multiple GL context + allow to share them !
+        # TODO make one or multiple VBO depending on len(vertices),
+        # TODO use a general common VBO for small amount of data
+        for name in self._unsyncAttributes:
+            array = self._attributes[name]
+            self._vbos[name] = ctx.glCtx.makeVboAttrib(array)
+        self._unsyncAttributes = []
+
+        if self._indices is not None and '__indices__' not in self._vbos:
+            vbo = ctx.glCtx.makeVbo(self._indices,
+                                    usage=gl.GL_STATIC_DRAW,
+                                    target=gl.GL_ELEMENT_ARRAY_BUFFER)
+            self._vbos['__indices__'] = vbo
+
+    def _draw(self, program=None, nbVertices=None):
+        """Perform OpenGL draw calls.
+
+        :param GLProgram program:
+            If not None, call :meth:`useAttribute` for this program.
+        :param int nbVertices:
+            The number of vertices to render or None to render all vertices.
+        """
+        if program is not None:
+            self.useAttribute(program)
+
+        if self._indices is None:
+            if nbVertices is None:
+                nbVertices = self.nbVertices
+            gl.glDrawArrays(self._MODES[self._mode], 0, nbVertices)
+        else:
+            if nbVertices is None:
+                nbVertices = self._indices.size
+            with self._vbos['__indices__']:
+                gl.glDrawElements(self._MODES[self._mode],
+                                  nbVertices,
+                                  _glutils.numpyToGLType(self._indices.dtype),
+                                  ctypes.c_void_p(0))
+
+
+# Lines #######################################################################
+
+class Lines(Geometry):
+    """A set of segments"""
+    _shaders = ("""
+    attribute vec3 position;
+    attribute vec3 normal;
+    attribute vec4 color;
+
+    uniform mat4 matrix;
+    uniform mat4 transformMat;
+
+    varying vec4 vCameraPosition;
+    varying vec3 vPosition;
+    varying vec3 vNormal;
+    varying vec4 vColor;
+
+    void main(void)
+    {
+        gl_Position = matrix * vec4(position, 1.0);
+        vCameraPosition = transformMat * vec4(position, 1.0);
+        vPosition = position;
+        vNormal = normal;
+        vColor = color;
+    }
+    """,
+                string.Template("""
+    varying vec4 vCameraPosition;
+    varying vec3 vPosition;
+    varying vec3 vNormal;
+    varying vec4 vColor;
+
+    $clippingDecl
+    $lightingFunction
+
+    void main(void)
+    {
+        $clippingCall(vCameraPosition);
+        gl_FragColor = $lightingCall(vColor, vPosition, vNormal);
+    }
+    """))
+
+    def __init__(self, positions, normals=None, colors=(1., 1., 1., 1.),
+                 indices=None, mode='lines', width=1.):
+        if mode == 'strip':
+            mode = 'line_strip'
+        assert mode in self._LINE_MODES
+
+        self._width = width
+        self._smooth = True
+
+        super(Lines, self).__init__(mode, indices,
+                                    position=positions,
+                                    normal=normals,
+                                    color=colors)
+
+    width = event.notifyProperty('_width', converter=float,
+                                 doc="Width of the line in pixels.")
+
+    smooth = event.notifyProperty(
+        '_smooth',
+        converter=bool,
+        doc="Smooth line rendering enabled (bool, default: True)")
+
+    def renderGL2(self, ctx):
+        # Prepare program
+        isnormals = 'normal' in self._attributes
+        if isnormals:
+            fraglightfunction = ctx.viewport.light.fragmentDef
+        else:
+            fraglightfunction = ctx.viewport.light.fragmentShaderFunctionNoop
+
+        fragment = self._shaders[1].substitute(
+            clippingDecl=ctx.clipper.fragDecl,
+            clippingCall=ctx.clipper.fragCall,
+            lightingFunction=fraglightfunction,
+            lightingCall=ctx.viewport.light.fragmentCall)
+        prog = ctx.glCtx.prog(self._shaders[0], fragment)
+        prog.use()
+
+        if isnormals:
+            ctx.viewport.light.setupProgram(ctx, prog)
+
+        gl.glLineWidth(self.width)
+
+        prog.setUniformMatrix('matrix', ctx.objectToNDC.matrix)
+        prog.setUniformMatrix('transformMat',
+                              ctx.objectToCamera.matrix,
+                              safe=True)
+
+        ctx.clipper.setupProgram(ctx, prog)
+
+        with gl.enabled(gl.GL_LINE_SMOOTH, self._smooth):
+            self._draw(prog)
+
+
+class DashedLines(Lines):
+    """Set of dashed lines
+
+    This MUST be defined as a set of lines (no strip or loop).
+    """
+
+    _shaders = ("""
+    attribute vec3 position;
+    attribute vec3 origin;
+    attribute vec3 normal;
+    attribute vec4 color;
+
+    uniform mat4 matrix;
+    uniform mat4 transformMat;
+    uniform vec2 viewportSize;  /* Width, height of the viewport */
+
+    varying vec4 vCameraPosition;
+    varying vec3 vPosition;
+    varying vec3 vNormal;
+    varying vec4 vColor;
+    varying vec2 vOriginFragCoord;
+
+    void main(void)
+    {
+        gl_Position = matrix * vec4(position, 1.0);
+        vCameraPosition = transformMat * vec4(position, 1.0);
+        vPosition = position;
+        vNormal = normal;
+        vColor = color;
+
+        vec4 clipOrigin = matrix * vec4(origin, 1.0);
+        vec4 ndcOrigin = clipOrigin / clipOrigin.w;  /* Perspective divide */
+        /* Convert to same frame as gl_FragCoord: lower-left, pixel center at 0.5, 0.5 */
+        vOriginFragCoord = (ndcOrigin.xy + vec2(1.0, 1.0)) * 0.5 * viewportSize + vec2(0.5, 0.5);
+    }
+    """,  # noqa
+                string.Template("""
+    varying vec4 vCameraPosition;
+    varying vec3 vPosition;
+    varying vec3 vNormal;
+    varying vec4 vColor;
+    varying vec2 vOriginFragCoord;
+
+    uniform vec2 dash;
+
+    $clippingDecl
+    $lightingFunction
+
+    void main(void)
+    {
+        /* Discard off dash fragments */
+        float lineDist = distance(vOriginFragCoord, gl_FragCoord.xy);
+        if (mod(lineDist, dash.x + dash.y) > dash.x) {
+            discard;
+        }
+        $clippingCall(vCameraPosition);
+        gl_FragColor = $lightingCall(vColor, vPosition, vNormal);
+    }
+    """))
+
+    def __init__(self, positions, colors=(1., 1., 1., 1.),
+                 indices=None, width=1.):
+        self._dash = 1, 0
+        super(DashedLines, self).__init__(positions=positions,
+                                          colors=colors,
+                                          indices=indices,
+                                          mode='lines',
+                                          width=width)
+
+    @property
+    def dash(self):
+        """Dash of the line as a 2-tuple of lengths in pixels: (on, off)"""
+        return self._dash
+
+    @dash.setter
+    def dash(self, dash):
+        dash = float(dash[0]), float(dash[1])
+        if dash != self._dash:
+            self._dash = dash
+            self.notify()
+
+    def getPositions(self, copy=True):
+        """Get coordinates of lines.
+
+        :param bool copy: True to get a copy, False otherwise
+        :returns: Coordinates of lines
+        :rtype: numpy.ndarray of float32 of shape (N, 2, Ndim)
+        """
+        return self.getAttribute('position', copy=copy)
+
+    def setPositions(self, positions, copy=True):
+        """Set line coordinates.
+
+        :param positions: Array of line coordinates
+        :param bool copy: True to copy input array, False to use as is
+        """
+        self.setAttribute('position', positions, copy=copy)
+        # Update line origins from given positions
+        origins = numpy.array(positions, copy=True, order='C')
+        origins[1::2] = origins[::2]
+        self.setAttribute('origin', origins, copy=False)
+
+    def renderGL2(self, context):
+        # Prepare program
+        isnormals = 'normal' in self._attributes
+        if isnormals:
+            fraglightfunction = context.viewport.light.fragmentDef
+        else:
+            fraglightfunction = \
+                context.viewport.light.fragmentShaderFunctionNoop
+
+        fragment = self._shaders[1].substitute(
+            clippingDecl=context.clipper.fragDecl,
+            clippingCall=context.clipper.fragCall,
+            lightingFunction=fraglightfunction,
+            lightingCall=context.viewport.light.fragmentCall)
+        program = context.glCtx.prog(self._shaders[0], fragment)
+        program.use()
+
+        if isnormals:
+            context.viewport.light.setupProgram(context, program)
+
+        gl.glLineWidth(self.width)
+
+        program.setUniformMatrix('matrix', context.objectToNDC.matrix)
+        program.setUniformMatrix('transformMat',
+                                 context.objectToCamera.matrix,
+                                 safe=True)
+
+        gl.glUniform2f(
+            program.uniforms['viewportSize'], *context.viewport.size)
+        gl.glUniform2f(program.uniforms['dash'], *self.dash)
+
+        context.clipper.setupProgram(context, program)
+
+        self._draw(program)
+
+
+class Box(core.PrivateGroup):
+    """Rectangular box"""
+
+    _lineIndices = numpy.array((
+        (0, 1), (1, 2), (2, 3), (3, 0),  # Lines with z=0
+        (0, 4), (1, 5), (2, 6), (3, 7),  # Lines from z=0 to z=1
+        (4, 5), (5, 6), (6, 7), (7, 4)),  # Lines with z=1
+        dtype=numpy.uint8)
+
+    _faceIndices = numpy.array(
+        (0, 3, 1, 2, 5, 6, 4, 7, 7, 6, 6, 2, 7, 3, 4, 0, 5, 1),
+        dtype=numpy.uint8)
+
+    _vertices = numpy.array((
+        # Corners with z=0
+        (0., 0., 0.), (1., 0., 0.), (1., 1., 0.), (0., 1., 0.),
+        # Corners with z=1
+        (0., 0., 1.), (1., 0., 1.), (1., 1., 1.), (0., 1., 1.)),
+        dtype=numpy.float32)
+
+    def __init__(self, size=(1., 1., 1.),
+                 stroke=(1., 1., 1., 1.),
+                 fill=(1., 1., 1., 0.)):
+        super(Box, self).__init__()
+
+        self._fill = Mesh3D(self._vertices,
+                            colors=rgba(fill),
+                            mode='triangle_strip',
+                            indices=self._faceIndices)
+        self._fill.visible = self.fillColor[-1] != 0.
+
+        self._stroke = Lines(self._vertices,
+                             indices=self._lineIndices,
+                             colors=rgba(stroke),
+                             mode='lines')
+        self._stroke.visible = self.strokeColor[-1] != 0.
+        self.strokeWidth = 1.
+
+        self._children = [self._stroke, self._fill]
+
+        self._size = None
+        self.size = size
+
+    @property
+    def size(self):
+        """Size of the box (sx, sy, sz)"""
+        return self._size
+
+    @size.setter
+    def size(self, size):
+        assert len(size) == 3
+        size = tuple(size)
+        if size != self.size:
+            self._size = size
+            self._fill.setAttribute(
+                'position',
+                self._vertices * numpy.array(size, dtype=numpy.float32))
+            self._stroke.setAttribute(
+                'position',
+                self._vertices * numpy.array(size, dtype=numpy.float32))
+            self.notify()
+
+    @property
+    def strokeSmooth(self):
+        """True to draw smooth stroke, False otherwise"""
+        return self._stroke.smooth
+
+    @strokeSmooth.setter
+    def strokeSmooth(self, smooth):
+        smooth = bool(smooth)
+        if smooth != self._stroke.smooth:
+            self._stroke.smooth = smooth
+            self.notify()
+
+    @property
+    def strokeWidth(self):
+        """Width of the stroke (float)"""
+        return self._stroke.width
+
+    @strokeWidth.setter
+    def strokeWidth(self, width):
+        width = float(width)
+        if width != self.strokeWidth:
+            self._stroke.width = width
+            self.notify()
+
+    @property
+    def strokeColor(self):
+        """RGBA color of the box lines (4-tuple of float in [0, 1])"""
+        return tuple(self._stroke.getAttribute('color', copy=False))
+
+    @strokeColor.setter
+    def strokeColor(self, color):
+        color = rgba(color)
+        if color != self.strokeColor:
+            self._stroke.setAttribute('color', color)
+            # Fully transparent = hidden
+            self._stroke.visible = color[-1] != 0.
+            self.notify()
+
+    @property
+    def fillColor(self):
+        """RGBA color of the box faces (4-tuple of float in [0, 1])"""
+        return tuple(self._fill.getAttribute('color', copy=False))
+
+    @fillColor.setter
+    def fillColor(self, color):
+        color = rgba(color)
+        if color != self.fillColor:
+            self._fill.setAttribute('color', color)
+            # Fully transparent = hidden
+            self._fill.visible = color[-1] != 0.
+            self.notify()
+
+    @property
+    def fillCulling(self):
+        return self._fill.culling
+
+    @fillCulling.setter
+    def fillCulling(self, culling):
+        self._fill.culling = culling
+
+
+class Axes(Lines):
+    """3D RGB orthogonal axes"""
+    _vertices = numpy.array(((0., 0., 0.), (1., 0., 0.),
+                             (0., 0., 0.), (0., 1., 0.),
+                             (0., 0., 0.), (0., 0., 1.)),
+                            dtype=numpy.float32)
+
+    _colors = numpy.array(((255, 0, 0, 255), (255, 0, 0, 255),
+                           (0, 255, 0, 255), (0, 255, 0, 255),
+                           (0, 0, 255, 255), (0, 0, 255, 255)),
+                          dtype=numpy.uint8)
+
+    def __init__(self):
+        super(Axes, self).__init__(self._vertices,
+                                   colors=self._colors,
+                                   width=3.)
+
+
+class BoxWithAxes(Lines):
+    """Rectangular box with RGB OX, OY, OZ axes
+
+    :param color: RGBA color of the box
+    """
+
+    _vertices = numpy.array((
+        # Axes corners
+        (0., 0., 0.), (1., 0., 0.),
+        (0., 0., 0.), (0., 1., 0.),
+        (0., 0., 0.), (0., 0., 1.),
+        # Box corners with z=0
+        (1., 0., 0.), (1., 1., 0.), (0., 1., 0.),
+        # Box corners with z=1
+        (0., 0., 1.), (1., 0., 1.), (1., 1., 1.), (0., 1., 1.)),
+        dtype=numpy.float32)
+
+    _axesColors = numpy.array(((1., 0., 0., 1.), (1., 0., 0., 1.),
+                               (0., 1., 0., 1.), (0., 1., 0., 1.),
+                               (0., 0., 1., 1.), (0., 0., 1., 1.)),
+                              dtype=numpy.float32)
+
+    _lineIndices = numpy.array((
+        (0, 1), (2, 3), (4, 5),  # Axes lines
+        (6, 7), (7, 8),  # Box lines with z=0
+        (6, 10), (7, 11), (8, 12),  # Box lines from z=0 to z=1
+        (9, 10), (10, 11), (11, 12), (12, 9)),  # Box lines with z=1
+        dtype=numpy.uint8)
+
+    def __init__(self, color=(1., 1., 1., 1.)):
+        self._color = (1., 1., 1., 1.)
+        colors = numpy.ones((len(self._vertices), 4), dtype=numpy.float32)
+        colors[:len(self._axesColors), :] = self._axesColors
+
+        super(BoxWithAxes, self).__init__(self._vertices,
+                                          indices=self._lineIndices,
+                                          colors=colors,
+                                          width=2.)
+        self.color = color
+
+    @property
+    def color(self):
+        """The RGBA color to use for the box: 4 float in [0, 1]"""
+        return self._color
+
+    @color.setter
+    def color(self, color):
+        color = rgba(color)
+        if color != self._color:
+            self._color = color
+            colors = numpy.empty((len(self._vertices), 4), dtype=numpy.float32)
+            colors[:len(self._axesColors), :] = self._axesColors
+            colors[len(self._axesColors):, :] = self._color
+            self.setAttribute('color', colors)  # Do the notification
+
+
+class PlaneInGroup(core.PrivateGroup):
+    """A plane using its parent bounds to display a contour.
+
+    If plane is outside the bounds of its parent, it is not visible.
+
+    Cannot set the transform attribute of this primitive.
+    This primitive never has any bounds.
+    """
+    # TODO inherit from Lines directly?, make sure the plane remains visible?
+
+    def __init__(self, point=(0., 0., 0.), normal=(0., 0., 1.)):
+        super(PlaneInGroup, self).__init__()
+        self._cache = None, None  # Store bounds, vertices
+        self._outline = None
+
+        self._color = None
+        self.color = 1., 1., 1., 1.  # Set _color
+        self._width = 2.
+
+        self._plane = utils.Plane(point, normal)
+        self._plane.addListener(self._planeChanged)
+
+    def moveToCenter(self):
+        """Place the plane at the center of the data, not changing orientation.
+        """
+        if self.parent is not None:
+            bounds = self.parent.bounds(dataBounds=True)
+            if bounds is not None:
+                center = (bounds[0] + bounds[1]) / 2.
+                _logger.debug('Moving plane to center: %s', str(center))
+                self.plane.point = center
+
+    @property
+    def color(self):
+        """Plane outline color (array of 4 float in [0, 1])."""
+        return self._color.copy()
+
+    @color.setter
+    def color(self, color):
+        self._color = numpy.array(color, copy=True, dtype=numpy.float32)
+        if self._outline is not None:
+            self._outline.setAttribute('color', self._color)
+        self.notify()  # This is OK as Lines are rebuild for each rendering
+
+    @property
+    def width(self):
+        """Width of the plane stroke in pixels"""
+        return self._width
+
+    @width.setter
+    def width(self, width):
+        self._width = float(width)
+        if self._outline is not None:
+            self._outline.width = self._width  # Sync width
+
+    # Plane access
+
+    @property
+    def plane(self):
+        """The plane parameters in the frame of the object."""
+        return self._plane
+
+    def _planeChanged(self, source):
+        """Listener of plane changes: clear cache and notify listeners."""
+        self._cache = None, None
+        self.notify()
+
+    # Disable some scene features
+
+    @property
+    def transforms(self):
+        # Ready-only transforms to prevent using it
+        return self._transforms
+
+    def _bounds(self, dataBounds=False):
+        # This is bound less as it uses the bounds of its parent.
+        return None
+
+    @property
+    def contourVertices(self):
+        """The vertices of the contour of the plane/bounds intersection."""
+        parent = self.parent
+        if parent is None:
+            return None  # No parent: no vertices
+
+        bounds = parent.bounds(dataBounds=True)
+        if bounds is None:
+            return None  # No bounds: no vertices
+
+        # Check if cache is valid and return it
+        cachebounds, cachevertices = self._cache
+        if numpy.all(numpy.equal(bounds, cachebounds)):
+            return cachevertices
+
+        # Cache is not OK, rebuild it
+        boxvertices = bounds[0] + Box._vertices.copy()*(bounds[1] - bounds[0])
+        lineindices = Box._lineIndices
+        vertices = utils.boxPlaneIntersect(
+            boxvertices, lineindices, self.plane.normal, self.plane.point)
+
+        self._cache = bounds, vertices if len(vertices) != 0 else None
+
+        return self._cache[1]
+
+    @property
+    def center(self):
+        """The center of the plane/bounds intersection points."""
+        if not self.isValid:
+            return None
+        else:
+            return numpy.mean(self.contourVertices, axis=0)
+
+    @property
+    def isValid(self):
+        """True if a contour is defined, False otherwise."""
+        return self.plane.isPlane and self.contourVertices is not None
+
+    def prepareGL2(self, ctx):
+        if self.isValid:
+            if self._outline is None:  # Init outline
+                self._outline = Lines(self.contourVertices,
+                                      mode='loop',
+                                      colors=self.color)
+                self._outline.width = self._width
+                self._children.append(self._outline)
+
+            # Update vertices, TODO only when necessary
+            self._outline.setAttribute('position', self.contourVertices)
+
+            super(PlaneInGroup, self).prepareGL2(ctx)
+
+    def renderGL2(self, ctx):
+        if self.isValid:
+            super(PlaneInGroup, self).renderGL2(ctx)
+
+
+# Points ######################################################################
+
+_POINTS_ATTR_INFO = Geometry._ATTR_INFO.copy()
+_POINTS_ATTR_INFO.update(value={'dims': (1, 2), 'lastDim': (1,)},
+                         size={'dims': (1, 2), 'lastDim': (1,)},
+                         symbol={'dims': (1, 2), 'lastDim': (1,)})
+
+
+class Points(Geometry):
+    """A set of data points with an associated value and size."""
+    _shaders = ("""
+    #version 120
+
+    attribute vec3 position;
+    attribute float symbol;
+    attribute float value;
+    attribute float size;
+
+    uniform mat4 matrix;
+    uniform mat4 transformMat;
+
+    uniform vec2 valRange;
+
+    varying vec4 vCameraPosition;
+    varying float vSymbol;
+    varying float vNormValue;
+    varying float vSize;
+
+    void main(void)
+    {
+        vSymbol = symbol;
+
+        vNormValue = clamp((value - valRange.x) / (valRange.y - valRange.x),
+                           0.0, 1.0);
+
+        bool isValueInRange = value >= valRange.x && value <= valRange.y;
+        if (isValueInRange) {
+            gl_Position = matrix * vec4(position, 1.0);
+        } else {
+            gl_Position = vec4(2.0, 0.0, 0.0, 1.0); /* Get clipped */
+        }
+        vCameraPosition = transformMat * vec4(position, 1.0);
+
+        gl_PointSize = size;
+        vSize = size;
+    }
+    """,
+                string.Template("""
+    #version 120
+
+    varying vec4 vCameraPosition;
+    varying float vSize;
+    varying float vSymbol;
+    varying float vNormValue;
+
+    $clippinDecl
+
+    /* Circle */
+    #define SYMBOL_CIRCLE 1.0
+
+    float alphaCircle(vec2 coord, float size) {
+        float radius = 0.5;
+        float r = distance(coord, vec2(0.5, 0.5));
+        return clamp(size * (radius - r), 0.0, 1.0);
+    }
+
+    /* Half lines */
+    #define SYMBOL_H_LINE 2.0
+    #define LEFT 1.0
+    #define RIGHT 2.0
+    #define SYMBOL_V_LINE 3.0
+    #define UP 1.0
+    #define DOWN 2.0
+
+    float alphaLine(vec2 coord, float size, float direction)
+    {
+        vec2 delta = abs(size * (coord - 0.5));
+
+        if (direction == SYMBOL_H_LINE) {
+            return (delta.y < 0.5) ? 1.0 : 0.0;
+        }
+        else if (direction == SYMBOL_H_LINE + LEFT) {
+            return (coord.x <= 0.5 && delta.y < 0.5) ? 1.0 : 0.0;
+        }
+        else if (direction == SYMBOL_H_LINE + RIGHT) {
+            return (coord.x >= 0.5 && delta.y < 0.5) ? 1.0 : 0.0;
+        }
+        else if (direction == SYMBOL_V_LINE) {
+            return (delta.x < 0.5) ? 1.0 : 0.0;
+        }
+        else if (direction == SYMBOL_V_LINE + UP) {
+            return (coord.y <= 0.5 && delta.x < 0.5) ? 1.0 : 0.0;
+        }
+        else if (direction == SYMBOL_V_LINE + DOWN) {
+             return (coord.y >= 0.5 && delta.x < 0.5) ? 1.0 : 0.0;
+        }
+        return 1.0;
+    }
+
+    void main(void)
+    {
+        $clippingCall(vCameraPosition);
+
+        gl_FragColor = vec4(0.5 * vNormValue + 0.5, 0.0, 0.0, 1.0);
+
+        float alpha = 1.0;
+        float symbol = floor(vSymbol);
+        if (1 == 1) { //symbol == SYMBOL_CIRCLE) {
+            alpha = alphaCircle(gl_PointCoord, vSize);
+        }
+        else if (symbol >= SYMBOL_H_LINE &&
+                 symbol <= (SYMBOL_V_LINE + DOWN)) {
+            alpha = alphaLine(gl_PointCoord, vSize, symbol);
+        }
+        if (alpha == 0.0) {
+            discard;
+        }
+        gl_FragColor.a *= alpha;
+    }
+    """))
+
+    _ATTR_INFO = _POINTS_ATTR_INFO
+
+    # TODO Add colormap, light?
+
+    def __init__(self, vertices, values=0., sizes=1., indices=None,
+                 symbols=0.,
+                 minValue=None, maxValue=None):
+        super(Points, self).__init__('points', indices,
+                                     position=vertices,
+                                     value=values,
+                                     size=sizes,
+                                     symbol=symbols)
+
+        values = self._attributes['value']
+        self._minValue = values.min() if minValue is None else minValue
+        self._maxValue = values.max() if maxValue is None else maxValue
+
+    minValue = event.notifyProperty('_minValue')
+    maxValue = event.notifyProperty('_maxValue')
+
+    def renderGL2(self, ctx):
+        fragment = self._shaders[1].substitute(
+            clippingDecl=ctx.clipper.fragDecl,
+            clippingCall=ctx.clipper.fragCall)
+        prog = ctx.glCtx.prog(self._shaders[0], fragment)
+        prog.use()
+
+        gl.glEnable(gl.GL_VERTEX_PROGRAM_POINT_SIZE)  # OpenGL 2
+        gl.glEnable(gl.GL_POINT_SPRITE)  # OpenGL 2
+        # gl.glEnable(gl.GL_PROGRAM_POINT_SIZE)
+
+        prog.setUniformMatrix('matrix', ctx.objectToNDC.matrix)
+        prog.setUniformMatrix('transformMat',
+                              ctx.objectToCamera.matrix,
+                              safe=True)
+
+        ctx.clipper.setupProgram(ctx, prog)
+
+        gl.glUniform2f(prog.uniforms['valRange'], self.minValue, self.maxValue)
+
+        self._draw(prog)
+
+
+class ColorPoints(Geometry):
+    """A set of points with an associated color and size."""
+
+    _shaders = ("""
+    #version 120
+
+    attribute vec3 position;
+    attribute float symbol;
+    attribute vec4 color;
+    attribute float size;
+
+    uniform mat4 matrix;
+    uniform mat4 transformMat;
+
+    varying vec4 vCameraPosition;
+    varying float vSymbol;
+    varying vec4 vColor;
+    varying float vSize;
+
+    void main(void)
+    {
+        vCameraPosition = transformMat * vec4(position, 1.0);
+        vSymbol = symbol;
+        vColor = color;
+        gl_Position = matrix * vec4(position, 1.0);
+        gl_PointSize = size;
+        vSize = size;
+    }
+    """,
+                string.Template("""
+    #version 120
+
+    varying vec4 vCameraPosition;
+    varying float vSize;
+    varying float vSymbol;
+    varying vec4 vColor;
+
+    $clippingDecl;
+
+    /* Circle */
+    #define SYMBOL_CIRCLE 1.0
+
+    float alphaCircle(vec2 coord, float size) {
+        float radius = 0.5;
+        float r = distance(coord, vec2(0.5, 0.5));
+        return clamp(size * (radius - r), 0.0, 1.0);
+    }
+
+    /* Half lines */
+    #define SYMBOL_H_LINE 2.0
+    #define LEFT 1.0
+    #define RIGHT 2.0
+    #define SYMBOL_V_LINE 3.0
+    #define UP 1.0
+    #define DOWN 2.0
+
+    float alphaLine(vec2 coord, float size, float direction)
+    {
+        vec2 delta = abs(size * (coord - 0.5));
+
+        if (direction == SYMBOL_H_LINE) {
+            return (delta.y < 0.5) ? 1.0 : 0.0;
+        }
+        else if (direction == SYMBOL_H_LINE + LEFT) {
+            return (coord.x <= 0.5 && delta.y < 0.5) ? 1.0 : 0.0;
+        }
+        else if (direction == SYMBOL_H_LINE + RIGHT) {
+            return (coord.x >= 0.5 && delta.y < 0.5) ? 1.0 : 0.0;
+        }
+        else if (direction == SYMBOL_V_LINE) {
+            return (delta.x < 0.5) ? 1.0 : 0.0;
+        }
+        else if (direction == SYMBOL_V_LINE + UP) {
+            return (coord.y <= 0.5 && delta.x < 0.5) ? 1.0 : 0.0;
+        }
+        else if (direction == SYMBOL_V_LINE + DOWN) {
+             return (coord.y >= 0.5 && delta.x < 0.5) ? 1.0 : 0.0;
+        }
+        return 1.0;
+    }
+
+    void main(void)
+    {
+        $clippingCall(vCameraPosition);
+
+        gl_FragColor = vColor;
+
+        float alpha = 1.0;
+        float symbol = floor(vSymbol);
+        if (1 == 1) { //symbol == SYMBOL_CIRCLE) {
+            alpha = alphaCircle(gl_PointCoord, vSize);
+        }
+        else if (symbol >= SYMBOL_H_LINE &&
+                 symbol <= (SYMBOL_V_LINE + DOWN)) {
+            alpha = alphaLine(gl_PointCoord, vSize, symbol);
+        }
+        if (alpha == 0.0) {
+            discard;
+        }
+        gl_FragColor.a *= alpha;
+    }
+    """))
+
+    _ATTR_INFO = _POINTS_ATTR_INFO
+
+    def __init__(self, vertices, colors=(1., 1., 1., 1.), sizes=1.,
+                 indices=None, symbols=0.,
+                 minValue=None, maxValue=None):
+        super(ColorPoints, self).__init__('points', indices,
+                                          position=vertices,
+                                          color=colors,
+                                          size=sizes,
+                                          symbol=symbols)
+
+    def renderGL2(self, ctx):
+        fragment = self._shaders[1].substitute(
+            clippingDecl=ctx.clipper.fragDecl,
+            clippingCall=ctx.clipper.fragCall)
+        prog = ctx.glCtx.prog(self._shaders[0], fragment)
+        prog.use()
+
+        gl.glEnable(gl.GL_VERTEX_PROGRAM_POINT_SIZE)  # OpenGL 2
+        gl.glEnable(gl.GL_POINT_SPRITE)  # OpenGL 2
+        # gl.glEnable(gl.GL_PROGRAM_POINT_SIZE)
+
+        prog.setUniformMatrix('matrix', ctx.objectToNDC.matrix)
+        prog.setUniformMatrix('transformMat',
+                              ctx.objectToCamera.matrix,
+                              safe=True)
+
+        ctx.clipper.setupProgram(ctx, prog)
+
+        self._draw(prog)
+
+
+class GridPoints(Geometry):
+    # GLSL 1.30 !
+    """Data points on a regular grid with an associated value and size."""
+    _shaders = ("""
+    #version 130
+
+    in float value;
+    in float size;
+
+    uniform ivec3 gridDims;
+    uniform mat4 matrix;
+    uniform mat4 transformMat;
+    uniform vec2 valRange;
+
+    out vec4 vCameraPosition;
+    out float vNormValue;
+
+    //ivec3 coordsFromIndex(int index, ivec3 shape)
+    //{
+        /*Assumes that data is stored as z-major, then y, contiguous on x
+        */
+   //     int yxPlaneSize = shape.y * shape.x; /* nb of elem in 2d yx plane */
+   //     int z = index / yxPlaneSize;
+   //     int yxIndex = index - z * yxPlaneSize; /* index in 2d yx plane */
+   //     int y = yxIndex / shape.x;
+   //     int x = yxIndex - y * shape.x;
+   //     return ivec3(x, y, z);
+   // }
+
+    ivec3 coordsFromIndex(int index, ivec3 shape)
+    {
+        /*Assumes that data is stored as x-major, then y, contiguous on z
+        */
+        int yzPlaneSize = shape.y * shape.z; /* nb of elem in 2d yz plane */
+        int x = index / yzPlaneSize;
+        int yzIndex = index - x * yzPlaneSize; /* index in 2d yz plane */
+        int y = yzIndex / shape.z;
+        int z = yzIndex - y * shape.z;
+        return ivec3(x, y, z);
+    }
+
+    void main(void)
+    {
+        vNormValue = clamp((value - valRange.x) / (valRange.y - valRange.x),
+                           0.0, 1.0);
+
+        bool isValueInRange = value >= valRange.x && value <= valRange.y;
+        if (isValueInRange) {
+            /* Retrieve 3D position from gridIndex */
+            vec3 coords = vec3(coordsFromIndex(gl_VertexID, gridDims));
+            vec3 position = coords / max(vec3(gridDims) - 1.0, 1.0);
+            gl_Position = matrix * vec4(position, 1.0);
+            vCameraPosition = transformMat * vec4(position, 1.0);
+        } else {
+            gl_Position = vec4(2.0, 0.0, 0.0, 1.0); /* Get clipped */
+            vCameraPosition = vec4(0.0, 0.0, 0.0, 0.0);
+        }
+
+        gl_PointSize = size;
+    }
+    """,
+                string.Template("""
+    #version 130
+
+    in vec4 vCameraPosition;
+    in float vNormValue;
+    out vec4 fragColor;
+
+    $clippingDecl
+
+    void main(void)
+    {
+        $clippingCall(vCameraPosition);
+
+        fragColor = vec4(0.5 * vNormValue + 0.5, 0.0, 0.0, 1.0);
+    }
+    """))
+
+    _ATTR_INFO = {
+        'value': {'dims': (1, 2), 'lastDim': (1,)},
+        'size': {'dims': (1, 2), 'lastDim': (1,)}
+    }
+
+    # TODO Add colormap, shape?
+    # TODO could also use a texture to store values
+
+    def __init__(self, values=0., shape=None, sizes=1., indices=None,
+                 minValue=None, maxValue=None):
+        if isinstance(values, collections.Iterable):
+            values = numpy.array(values, copy=False)
+
+            # Test if gl_VertexID will overflow
+            assert values.size < numpy.iinfo(numpy.int32).max
+
+            self._shape = values.shape
+            values = values.ravel()  # 1D to add as a 1D vertex attribute
+
+        else:
+            assert shape is not None
+            self._shape = tuple(shape)
+
+        assert len(self._shape) in (1, 2, 3)
+
+        super(GridPoints, self).__init__('points', indices,
+                                         value=values,
+                                         size=sizes)
+
+        data = self.getAttribute('value', copy=False)
+        self._minValue = data.min() if minValue is None else minValue
+        self._maxValue = data.max() if maxValue is None else maxValue
+
+    minValue = event.notifyProperty('_minValue')
+    maxValue = event.notifyProperty('_maxValue')
+
+    def _bounds(self, dataBounds=False):
+        # Get bounds from values shape
+        bounds = numpy.zeros((2, 3), dtype=numpy.float32)
+        bounds[1, :] = self._shape
+        bounds[1, :] -= 1
+        return bounds
+
+    def renderGL2(self, ctx):
+        fragment = self._shaders[1].substitute(
+            clippingDecl=ctx.clipper.fragDecl,
+            clippingCall=ctx.clipper.fragCall)
+        prog = ctx.glCtx.prog(self._shaders[0], fragment)
+        prog.use()
+
+        gl.glEnable(gl.GL_VERTEX_PROGRAM_POINT_SIZE)  # OpenGL 2
+        gl.glEnable(gl.GL_POINT_SPRITE)  # OpenGL 2
+        # gl.glEnable(gl.GL_PROGRAM_POINT_SIZE)
+
+        prog.setUniformMatrix('matrix', ctx.objectToNDC.matrix)
+        prog.setUniformMatrix('transformMat',
+                              ctx.objectToCamera.matrix,
+                              safe=True)
+
+        ctx.clipper.setupProgram(ctx, prog)
+
+        gl.glUniform3i(prog.uniforms['gridDims'],
+                       self._shape[2] if len(self._shape) == 3 else 1,
+                       self._shape[1] if len(self._shape) >= 2 else 1,
+                       self._shape[0])
+
+        gl.glUniform2f(prog.uniforms['valRange'], self.minValue, self.maxValue)
+
+        self._draw(prog, nbVertices=reduce(lambda a, b: a * b, self._shape))
+
+
+# Spheres #####################################################################
+
+class Spheres(Geometry):
+    """A set of spheres.
+
+    Spheres are rendered as circles using points.
+    This brings some limitations:
+    - Do not support non-uniform scaling.
+    - Assume the projection keeps ratio.
+    - Do not render distorion by perspective projection.
+    - If the sphere center is clipped, the whole sphere is not displayed.
+    """
+    # TODO check those links
+    # Accounting for perspective projection
+    # http://iquilezles.org/www/articles/sphereproj/sphereproj.htm
+
+    # Michael Mara and Morgan McGuire.
+    # 2D Polyhedral Bounds of a Clipped, Perspective-Projected 3D Sphere
+    # Journal of Computer Graphics Techniques, Vol. 2, No. 2, 2013.
+    # http://jcgt.org/published/0002/02/05/paper.pdf
+    # https://research.nvidia.com/publication/2d-polyhedral-bounds-clipped-perspective-projected-3d-sphere
+
+    # TODO some issues with small scaling and regular grid or due to sampling
+
+    _shaders = ("""
+    #version 120
+
+    attribute vec3 position;
+    attribute vec4 color;
+    attribute float radius;
+
+    uniform mat4 transformMat;
+    uniform mat4 projMat;
+    uniform vec2 screenSize;
+
+    varying vec4 vCameraPosition;
+    varying vec3 vPosition;
+    varying vec4 vColor;
+    varying float vViewDepth;
+    varying float vViewRadius;
+
+    void main(void)
+    {
+        vCameraPosition = transformMat * vec4(position, 1.0);
+        gl_Position = projMat * vCameraPosition;
+
+        vPosition = gl_Position.xyz / gl_Position.w;
+
+        /* From object space radius to view space diameter.
+         * Do not support non-uniform scaling */
+        vec4 viewSizeVector = transformMat * vec4(2.0 * radius, 0.0, 0.0, 0.0);
+        float viewSize = length(viewSizeVector.xyz);
+
+        /* Convert to pixel size at the xy center of the view space */
+        vec4 projSize = projMat * vec4(0.5 * viewSize, 0.0,
+                                       vCameraPosition.z, vCameraPosition.w);
+        gl_PointSize = max(1.0, screenSize[0] * projSize.x / projSize.w);
+
+        vColor = color;
+        vViewRadius = 0.5 * viewSize;
+        vViewDepth = vCameraPosition.z;
+    }
+    """,
+                string.Template("""
+    # version 120
+
+    uniform mat4 projMat;
+
+    varying vec4 vCameraPosition;
+    varying vec3 vPosition;
+    varying vec4 vColor;
+    varying float vViewDepth;
+    varying float vViewRadius;
+
+    $clippingDecl
+    $lightingFunction
+
+    void main(void)
+    {
+        $clippingCall(vCameraPosition);
+
+        /* Get normal from point coords */
+        vec3 normal;
+        normal.xy = 2.0 * gl_PointCoord - vec2(1.0);
+        normal.y *= -1.0; /*Invert y to match NDC orientation*/
+        float sqLength = dot(normal.xy, normal.xy);
+        if (sqLength > 1.0) { /* Length -> out of sphere */
+            discard;
+        }
+        normal.z = sqrt(1.0 - sqLength);
+
+        /*Lighting performed in NDC*/
+        /*TODO update this when lighting changed*/
+        //XXX vec3 position = vPosition + vViewRadius * normal;
+        gl_FragColor = $lightingCall(vColor, vPosition, normal);
+
+        /*Offset depth*/
+        float viewDepth = vViewDepth + vViewRadius * normal.z;
+        vec2 clipZW = viewDepth * projMat[2].zw + projMat[3].zw;
+        gl_FragDepth = 0.5 * (clipZW.x / clipZW.y) + 0.5;
+    }
+    """))
+
+    _ATTR_INFO = {
+        'position': {'dims': (2, ), 'lastDim': (2, 3, 4)},
+        'radius': {'dims': (1, 2), 'lastDim': (1, )},
+        'color': {'dims': (1, 2), 'lastDim': (3, 4)},
+    }
+
+    def __init__(self, positions, radius=1., colors=(1., 1., 1., 1.)):
+        self.__bounds = None
+        super(Spheres, self).__init__('points', None,
+                                      position=positions,
+                                      radius=radius,
+                                      color=colors)
+
+    def renderGL2(self, ctx):
+        fragment = self._shaders[1].substitute(
+            clippingDecl=ctx.clipper.fragDecl,
+            clippingCall=ctx.clipper.fragCall,
+            lightingFunction=ctx.viewport.light.fragmentDef,
+            lightingCall=ctx.viewport.light.fragmentCall)
+        prog = ctx.glCtx.prog(self._shaders[0], fragment)
+        prog.use()
+
+        ctx.viewport.light.setupProgram(ctx, prog)
+
+        gl.glEnable(gl.GL_VERTEX_PROGRAM_POINT_SIZE)  # OpenGL 2
+        gl.glEnable(gl.GL_POINT_SPRITE)  # OpenGL 2
+        # gl.glEnable(gl.GL_PROGRAM_POINT_SIZE)
+
+        prog.setUniformMatrix('projMat', ctx.projection.matrix)
+        prog.setUniformMatrix('transformMat',
+                              ctx.objectToCamera.matrix,
+                              safe=True)
+
+        ctx.clipper.setupProgram(ctx, prog)
+
+        gl.glUniform2f(prog.uniforms['screenSize'], *ctx.viewport.size)
+
+        self._draw(prog)
+
+    def _bounds(self, dataBounds=False):
+        if self.__bounds is None:
+            self.__bounds = numpy.zeros((2, 3), dtype=numpy.float32)
+            # Support vertex with to 2 to 4 coordinates
+            positions = self._attributes['position']
+            radius = self._attributes['radius']
+            self.__bounds[0, :positions.shape[1]] = \
+                (positions - radius).min(axis=0)[:3]
+            self.__bounds[1, :positions.shape[1]] = \
+                (positions + radius).max(axis=0)[:3]
+        return self.__bounds.copy()
+
+
+# Meshes ######################################################################
+
+class Mesh3D(Geometry):
+    """A conventional 3D mesh"""
+
+    _shaders = ("""
+    attribute vec3 position;
+    attribute vec3 normal;
+    attribute vec4 color;
+
+    uniform mat4 matrix;
+    uniform mat4 transformMat;
+    //uniform mat3 matrixInvTranspose;
+
+    varying vec4 vCameraPosition;
+    varying vec3 vPosition;
+    varying vec3 vNormal;
+    varying vec4 vColor;
+
+    void main(void)
+    {
+        vCameraPosition = transformMat * vec4(position, 1.0);
+        //vNormal = matrixInvTranspose * normalize(normal);
+        vPosition = position;
+        vNormal = normal;
+        vColor = color;
+        gl_Position = matrix * vec4(position, 1.0);
+    }
+    """,
+                string.Template("""
+    varying vec4 vCameraPosition;
+    varying vec3 vPosition;
+    varying vec3 vNormal;
+    varying vec4 vColor;
+
+    $clippingDecl
+    $lightingFunction
+
+    void main(void)
+    {
+        $clippingCall(vCameraPosition);
+
+        gl_FragColor = $lightingCall(vColor, vPosition, vNormal);
+    }
+    """))
+
+    def __init__(self,
+                 positions,
+                 colors,
+                 normals=None,
+                 mode='triangles',
+                 indices=None):
+        assert mode in self._TRIANGLE_MODES
+        super(Mesh3D, self).__init__(mode, indices,
+                                     position=positions,
+                                     normal=normals,
+                                     color=colors)
+
+        self._culling = None
+
+    @property
+    def culling(self):
+        """Face culling (str)
+
+        One of 'back', 'front' or None.
+        """
+        return self._culling
+
+    @culling.setter
+    def culling(self, culling):
+        assert culling in ('back', 'front', None)
+        if culling != self._culling:
+            self._culling = culling
+            self.notify()
+
+    def renderGL2(self, ctx):
+        isnormals = 'normal' in self._attributes
+        if isnormals:
+            fragLightFunction = ctx.viewport.light.fragmentDef
+        else:
+            fragLightFunction = ctx.viewport.light.fragmentShaderFunctionNoop
+
+        fragment = self._shaders[1].substitute(
+            clippingDecl=ctx.clipper.fragDecl,
+            clippingCall=ctx.clipper.fragCall,
+            lightingFunction=fragLightFunction,
+            lightingCall=ctx.viewport.light.fragmentCall)
+        prog = ctx.glCtx.prog(self._shaders[0], fragment)
+        prog.use()
+
+        if isnormals:
+            ctx.viewport.light.setupProgram(ctx, prog)
+
+        if self.culling is not None:
+            cullFace = gl.GL_FRONT if self.culling == 'front' else gl.GL_BACK
+            gl.glCullFace(cullFace)
+            gl.glEnable(gl.GL_CULL_FACE)
+
+        prog.setUniformMatrix('matrix', ctx.objectToNDC.matrix)
+        prog.setUniformMatrix('transformMat',
+                              ctx.objectToCamera.matrix,
+                              safe=True)
+
+        ctx.clipper.setupProgram(ctx, prog)
+
+        self._draw(prog)
+
+        if self.culling is not None:
+            gl.glDisable(gl.GL_CULL_FACE)
+
+
+# Group ######################################################################
+
+# TODO lighting, clipping as groups?
+# group composition?
+
+class GroupDepthOffset(core.Group):
+    """A group using 2-pass rendering and glDepthRange to avoid Z-fighting"""
+
+    def __init__(self, children=(), epsilon=None):
+        super(GroupDepthOffset, self).__init__(children)
+        self._epsilon = epsilon
+        self.isDepthRangeOn = True
+
+    def prepareGL2(self, ctx):
+        if self._epsilon is None:
+            depthbits = gl.glGetInteger(gl.GL_DEPTH_BITS)
+            self._epsilon = 1. / (1 << (depthbits - 1))
+
+    def renderGL2(self, ctx):
+        if self.isDepthRangeOn:
+            self._renderGL2WithDepthRange(ctx)
+        else:
+            super(GroupDepthOffset, self).renderGL2(ctx)
+
+    def _renderGL2WithDepthRange(self, ctx):
+        # gl.glDepthFunc(gl.GL_LESS)
+        with gl.enabled(gl.GL_CULL_FACE):
+            gl.glCullFace(gl.GL_BACK)
+            for child in self.children:
+                gl.glColorMask(
+                    gl.GL_FALSE, gl.GL_FALSE, gl.GL_FALSE, gl.GL_FALSE)
+                gl.glDepthMask(gl.GL_TRUE)
+                gl.glDepthRange(self._epsilon, 1.)
+
+                child.render(ctx)
+
+                gl.glColorMask(
+                    gl.GL_TRUE, gl.GL_TRUE, gl.GL_TRUE, gl.GL_TRUE)
+                gl.glDepthMask(gl.GL_FALSE)
+                gl.glDepthRange(0., 1. - self._epsilon)
+
+                child.render(ctx)
+
+            gl.glCullFace(gl.GL_FRONT)
+            for child in reversed(self.children):
+                gl.glColorMask(
+                    gl.GL_FALSE, gl.GL_FALSE, gl.GL_FALSE, gl.GL_FALSE)
+                gl.glDepthMask(gl.GL_TRUE)
+                gl.glDepthRange(self._epsilon, 1.)
+
+                child.render(ctx)
+
+                gl.glColorMask(
+                    gl.GL_TRUE, gl.GL_TRUE, gl.GL_TRUE, gl.GL_TRUE)
+                gl.glDepthMask(gl.GL_FALSE)
+                gl.glDepthRange(0., 1. - self._epsilon)
+
+                child.render(ctx)
+
+        gl.glDepthMask(gl.GL_TRUE)
+        gl.glDepthRange(0., 1.)
+        # gl.glDepthFunc(gl.GL_LEQUAL)
+        # TODO use epsilon for all rendering?
+        # TODO issue with picking in depth buffer!
+
+
+class GroupBBox(core.PrivateGroup):
+    """A group displaying a bounding box around the children."""
+
+    def __init__(self, children=(), color=(1., 1., 1., 1.)):
+        super(GroupBBox, self).__init__()
+        self._group = core.Group(children)
+
+        self._boxTransforms = transform.TransformList(
+            (transform.Translate(), transform.Scale()))
+
+        self._boxWithAxes = BoxWithAxes(color)
+        self._boxWithAxes.smooth = False
+        self._boxWithAxes.transforms = self._boxTransforms
+
+        self._children = [self._boxWithAxes, self._group]
+
+    def _updateBoxAndAxes(self):
+        """Update bbox and axes position and size according to children."""
+        bounds = self._group.bounds(dataBounds=True)
+        if bounds is not None:
+            origin = bounds[0]
+            scale = [(d if d != 0. else 1.) for d in bounds[1] - bounds[0]]
+        else:
+            origin, scale = (0., 0., 0.), (1., 1., 1.)
+
+        self._boxTransforms[0].translation = origin
+        self._boxTransforms[1].scale = scale
+
+    def _bounds(self, dataBounds=False):
+        self._updateBoxAndAxes()
+        return super(GroupBBox, self)._bounds(dataBounds)
+
+    def prepareGL2(self, ctx):
+        self._updateBoxAndAxes()
+        super(GroupBBox, self).prepareGL2(ctx)
+
+    # Give access to _group children
+
+    @property
+    def children(self):
+        return self._group.children
+
+    @children.setter
+    def children(self, iterable):
+        self._group.children = iterable
+
+    # Give access to box color
+
+    @property
+    def color(self):
+        """The RGBA color to use for the box: 4 float in [0, 1]"""
+        return self._boxWithAxes.color
+
+    @color.setter
+    def color(self, color):
+        self._boxWithAxes.color = color
+
+
+# Clipping Plane ##############################################################
+
+class ClipPlane(PlaneInGroup):
+    """A clipping plane attached to a box"""
+
+    def renderGL2(self, ctx):
+        super(ClipPlane, self).renderGL2(ctx)
+
+        if self.visible:
+            # Set-up clipping plane for following brothers
+
+            # No need of perspective divide, no projection
+            point = ctx.objectToCamera.transformPoint(self.plane.point,
+                                                      perspectiveDivide=False)
+            normal = ctx.objectToCamera.transformNormal(self.plane.normal)
+            ctx.setClipPlane(point, normal)
+
+    def postRender(self, ctx):
+        if self.visible:
+            # Disable clip planes
+            ctx.setClipPlane()
diff --git a/silx/gui/plot3d/scene/setup.py b/silx/gui/plot3d/scene/setup.py
new file mode 100644
index 0000000..ff4c0a6
--- /dev/null
+++ b/silx/gui/plot3d/scene/setup.py
@@ -0,0 +1,41 @@
+# coding: utf-8
+# /*##########################################################################
+#
+# Copyright (c) 2015-2017 European Synchrotron Radiation Facility
+#
+# Permission is hereby granted, free of charge, to any person obtaining a copy
+# of this software and associated documentation files (the "Software"), to deal
+# in the Software without restriction, including without limitation the rights
+# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+# copies of the Software, and to permit persons to whom the Software is
+# furnished to do so, subject to the following conditions:
+#
+# The above copyright notice and this permission notice shall be included in
+# all copies or substantial portions of the Software.
+#
+# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
+# THE SOFTWARE.
+#
+# ###########################################################################*/
+__authors__ = ["T. Vincent"]
+__license__ = "MIT"
+__date__ = "25/07/2016"
+
+from numpy.distutils.misc_util import Configuration
+
+
+def configuration(parent_package='', top_path=None):
+    config = Configuration('scene', parent_package, top_path)
+    config.add_subpackage('test')
+    return config
+
+
+if __name__ == "__main__":
+    from numpy.distutils.core import setup
+
+    setup(configuration=configuration)
diff --git a/silx/gui/plot3d/scene/test/__init__.py b/silx/gui/plot3d/scene/test/__init__.py
new file mode 100644
index 0000000..fc4621e
--- /dev/null
+++ b/silx/gui/plot3d/scene/test/__init__.py
@@ -0,0 +1,43 @@
+# coding: utf-8
+# /*##########################################################################
+#
+# Copyright (c) 2015-2017 European Synchrotron Radiation Facility
+#
+# Permission is hereby granted, free of charge, to any person obtaining a copy
+# of this software and associated documentation files (the "Software"), to deal
+# in the Software without restriction, including without limitation the rights
+# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+# copies of the Software, and to permit persons to whom the Software is
+# furnished to do so, subject to the following conditions:
+#
+# The above copyright notice and this permission notice shall be included in
+# all copies or substantial portions of the Software.
+#
+# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
+# THE SOFTWARE.
+#
+# ###########################################################################*/
+
+from __future__ import absolute_import, division, unicode_literals
+
+__authors__ = ["T. Vincent"]
+__license__ = "MIT"
+__date__ = "25/07/2016"
+
+
+import unittest
+
+from .test_transform import suite as test_transform_suite
+from .test_utils import suite as test_utils_suite
+
+
+def suite():
+    testsuite = unittest.TestSuite()
+    testsuite.addTest(test_transform_suite())
+    testsuite.addTest(test_utils_suite())
+    return testsuite
diff --git a/silx/gui/plot3d/scene/test/test_transform.py b/silx/gui/plot3d/scene/test/test_transform.py
new file mode 100644
index 0000000..9ea0af1
--- /dev/null
+++ b/silx/gui/plot3d/scene/test/test_transform.py
@@ -0,0 +1,91 @@
+# coding: utf-8
+# /*##########################################################################
+#
+# Copyright (c) 2015-2017 European Synchrotron Radiation Facility
+#
+# Permission is hereby granted, free of charge, to any person obtaining a copy
+# of this software and associated documentation files (the "Software"), to deal
+# in the Software without restriction, including without limitation the rights
+# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+# copies of the Software, and to permit persons to whom the Software is
+# furnished to do so, subject to the following conditions:
+#
+# The above copyright notice and this permission notice shall be included in
+# all copies or substantial portions of the Software.
+#
+# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
+# THE SOFTWARE.
+#
+# ###########################################################################*/
+
+from __future__ import absolute_import, division, unicode_literals
+
+__authors__ = ["T. Vincent"]
+__license__ = "MIT"
+__date__ = "05/01/2017"
+
+
+import numpy
+import unittest
+
+from silx.gui.plot3d.scene import transform
+
+
+class TestTransformList(unittest.TestCase):
+
+    def assertSameArrays(self, a, b):
+        return self.assertTrue(numpy.allclose(a, b, atol=1e-06))
+
+    def testTransformList(self):
+        """Minimalistic test of TransformList"""
+        transforms = transform.TransformList()
+        refmatrix = numpy.identity(4, dtype=numpy.float32)
+        self.assertSameArrays(refmatrix, transforms.matrix)
+
+        # Append translate
+        transforms.append(transform.Translate(1., 1., 1.))
+        refmatrix = numpy.array(((1., 0., 0., 1.),
+                                 (0., 1., 0., 1.),
+                                 (0., 0., 1., 1.),
+                                 (0., 0., 0., 1.)), dtype=numpy.float32)
+        self.assertSameArrays(refmatrix, transforms.matrix)
+
+        # Extend scale
+        transforms.extend([transform.Scale(0.1, 2., 1.)])
+        refmatrix = numpy.dot(refmatrix,
+                              numpy.array(((0.1, 0., 0., 0.),
+                                           (0., 2., 0., 0.),
+                                           (0., 0., 1., 0.),
+                                           (0., 0., 0., 1.)),
+                                          dtype=numpy.float32))
+        self.assertSameArrays(refmatrix, transforms.matrix)
+
+        # Insert rotate
+        transforms.insert(0, transform.Rotate(360.))
+        self.assertSameArrays(refmatrix, transforms.matrix)
+
+        # Update translate and check for listener called
+        self._callCount = 0
+
+        def listener(source):
+            self._callCount += 1
+        transforms.addListener(listener)
+
+        transforms[1].tx += 1
+        self.assertEqual(self._callCount, 1)
+
+
+def suite():
+    testsuite = unittest.TestSuite()
+    testsuite.addTest(
+        unittest.defaultTestLoader.loadTestsFromTestCase(TestTransformList))
+    return testsuite
+
+
+if __name__ == '__main__':
+    unittest.main(defaultTest='suite')
diff --git a/silx/gui/plot3d/scene/test/test_utils.py b/silx/gui/plot3d/scene/test/test_utils.py
new file mode 100644
index 0000000..65c2407
--- /dev/null
+++ b/silx/gui/plot3d/scene/test/test_utils.py
@@ -0,0 +1,275 @@
+# coding: utf-8
+# /*##########################################################################
+#
+# Copyright (c) 2015-2017 European Synchrotron Radiation Facility
+#
+# Permission is hereby granted, free of charge, to any person obtaining a copy
+# of this software and associated documentation files (the "Software"), to deal
+# in the Software without restriction, including without limitation the rights
+# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+# copies of the Software, and to permit persons to whom the Software is
+# furnished to do so, subject to the following conditions:
+#
+# The above copyright notice and this permission notice shall be included in
+# all copies or substantial portions of the Software.
+#
+# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
+# THE SOFTWARE.
+#
+# ###########################################################################*/
+
+from __future__ import absolute_import, division, unicode_literals
+
+__authors__ = ["T. Vincent"]
+__license__ = "MIT"
+__date__ = "25/07/2016"
+
+
+import unittest
+from silx.test.utils import ParametricTestCase
+
+import numpy
+
+from silx.gui.plot3d.scene import utils
+
+
+# angleBetweenVectors #########################################################
+
+class TestAngleBetweenVectors(ParametricTestCase):
+
+    TESTS = {  # name: (refvector, vectors, norm, refangles)
+        'single vector':
+            ((1., 0., 0.), (1., 0., 0.), (0., 0., 1.), 0.),
+        'single vector, no norm':
+            ((1., 0., 0.), (1., 0., 0.), None, 0.),
+
+        'with orthogonal norm':
+            ((1., 0., 0.),
+             ((1., 0., 0.), (0., 1., 0.), (-1., 0., 0.), (0., -1., 0.)),
+             (0., 0., 1.),
+             (0., 90., 180., 270.)),
+
+        'with coplanar norm':  # = similar to no norm
+            ((1., 0., 0.),
+             ((1., 0., 0.), (0., 1., 0.), (-1., 0., 0.), (0., -1., 0.)),
+             (1., 0., 0.),
+             (0., 90., 180., 90.)),
+
+        'without norm':
+            ((1., 0., 0.),
+             ((1., 0., 0.), (0., 1., 0.), (-1., 0., 0.), (0., -1., 0.)),
+             None,
+             (0., 90., 180., 90.)),
+
+        'not unit vectors':
+            ((2., 2., 0.), ((1., 1., 0.), (1., -1., 0.)), None, (0., 90.)),
+    }
+
+    def testAngleBetweenVectorsFunction(self):
+        for name, params in self.TESTS.items():
+            refvector, vectors, norm, refangles = params
+            with self.subTest(name):
+                refangles = numpy.radians(refangles)
+
+                refvector = numpy.array(refvector)
+                vectors = numpy.array(vectors)
+                if norm is not None:
+                    norm = numpy.array(norm)
+
+                testangles = utils.angleBetweenVectors(
+                    refvector, vectors, norm)
+
+                self.assertTrue(
+                    numpy.allclose(testangles, refangles, atol=1e-5))
+
+
+# Plane #######################################################################
+
+class AssertNotificationContext(object):
+    """Context that checks if an event.Notifier is sending events."""
+
+    def __init__(self, notifier, count=1):
+        """Initializer.
+
+        :param event.Notifier notifier: The notifier to test.
+        :param int count: The expected number of calls.
+        """
+        self._notifier = notifier
+        self._callCount = None
+        self._count = count
+
+    def __enter__(self):
+        self._callCount = 0
+        self._notifier.addListener(self._callback)
+
+    def __exit__(self, exc_type, exc_value, traceback):
+        # Do not return True so exceptions are propagated
+        self._notifier.removeListener(self._callback)
+        assert self._callCount == self._count
+        self._callCount = None
+
+    def _callback(self, *args, **kwargs):
+        self._callCount += 1
+
+
+class TestPlaneParameters(ParametricTestCase):
+    """Test Plane.parameters read/write and notifications."""
+
+    PARAMETERS = {
+        'unit normal': (1., 0., 0., 1.),
+        'not unit normal': (1., 1., 0., 1.),
+        'd = 0': (1., 0., 0., 0.)
+    }
+
+    def testParameters(self):
+        """Check parameters read/write and notification."""
+        plane = utils.Plane()
+
+        for name, parameters in self.PARAMETERS.items():
+            with self.subTest(name, parameters=parameters):
+                with AssertNotificationContext(plane):
+                    plane.parameters = parameters
+
+                # Plane parameters are converted to have a unit normal
+                normparams = parameters / numpy.linalg.norm(parameters[:3])
+                self.assertTrue(numpy.allclose(plane.parameters, normparams))
+
+    ZEROS_PARAMETERS = (
+        (0., 0., 0., 0.),
+        (0., 0., 0., 1.)
+    )
+
+    ZEROS = 0., 0., 0., 0.
+
+    def testParametersNoPlane(self):
+        """Test Plane.parameters with ||normal|| == 0 ."""
+        plane = utils.Plane()
+        plane.parameters = self.ZEROS
+
+        for parameters in self.ZEROS_PARAMETERS:
+            with self.subTest(parameters=parameters):
+                with AssertNotificationContext(plane, count=0):
+                    plane.parameters = parameters
+                self.assertTrue(
+                    numpy.allclose(plane.parameters, self.ZEROS, 0., 0.))
+
+
+# unindexArrays ###############################################################
+
+class TestUnindexArrays(ParametricTestCase):
+    """Test unindexArrays function."""
+
+    def testBasicModes(self):
+        """Test for modes: points, lines and triangles"""
+        indices = numpy.array((1, 2, 0))
+        arrays = (numpy.array((0., 1., 2.)),
+                  numpy.array(((0, 0), (1, 1), (2, 2))))
+        refresults = (numpy.array((1., 2., 0.)),
+                      numpy.array(((1, 1), (2, 2), (0, 0))))
+
+        for mode in ('points', 'lines', 'triangles'):
+            with self.subTest(mode=mode):
+                testresults = utils.unindexArrays(mode, indices, *arrays)
+                for ref, test in zip(refresults, testresults):
+                    self.assertTrue(numpy.equal(ref, test).all())
+
+    def testPackedLines(self):
+        """Test for modes: line_strip, loop"""
+        indices = numpy.array((1, 2, 0))
+        arrays = (numpy.array((0., 1., 2.)),
+                  numpy.array(((0, 0), (1, 1), (2, 2))))
+        results = {
+            'line_strip': (
+                numpy.array((1., 2., 2., 0.)),
+                numpy.array(((1, 1), (2, 2), (2, 2), (0, 0)))),
+            'loop': (
+                numpy.array((1., 2., 2., 0., 0., 1.)),
+                numpy.array(((1, 1), (2, 2), (2, 2), (0, 0), (0, 0), (1, 1)))),
+        }
+
+        for mode, refresults in results.items():
+            with self.subTest(mode=mode):
+                testresults = utils.unindexArrays(mode, indices, *arrays)
+                for ref, test in zip(refresults, testresults):
+                    self.assertTrue(numpy.equal(ref, test).all())
+
+    def testPackedTriangles(self):
+        """Test for modes: triangle_strip, fan"""
+        indices = numpy.array((1, 2, 0, 3))
+        arrays = (numpy.array((0., 1., 2., 3.)),
+                  numpy.array(((0, 0), (1, 1), (2, 2), (3, 3))))
+        results = {
+            'triangle_strip': (
+                numpy.array((1., 2., 0., 2., 0., 3.)),
+                numpy.array(((1, 1), (2, 2), (0, 0), (2, 2), (0, 0), (3, 3)))),
+            'fan': (
+                numpy.array((1., 2., 0., 1., 0., 3.)),
+                numpy.array(((1, 1), (2, 2), (0, 0), (1, 1), (0, 0), (3, 3)))),
+        }
+
+        for mode, refresults in results.items():
+            with self.subTest(mode=mode):
+                testresults = utils.unindexArrays(mode, indices, *arrays)
+                for ref, test in zip(refresults, testresults):
+                    self.assertTrue(numpy.equal(ref, test).all())
+
+    def testBadIndices(self):
+        """Test with negative indices and indices higher than array length"""
+        arrays = numpy.array((0, 1)), numpy.array((0, 1, 2))
+
+        # negative indices
+        with self.assertRaises(AssertionError):
+            utils.unindexArrays('points', (-1, 0), *arrays)
+
+        # Too high indices
+        with self.assertRaises(AssertionError):
+            utils.unindexArrays('points', (0, 10), *arrays)
+
+
+# triangleNormals #############################################################
+
+class TestTriangleNormals(ParametricTestCase):
+    """Test triangleNormals function."""
+
+    def test(self):
+        """Test for modes: points, lines and triangles"""
+        positions = numpy.array(
+            ((0., 0., 0.), (1., 0., 0.), (0., 1., 0.),  # normal = Z
+             (1., 1., 1.), (1., 2., 3.), (4., 5., 6.),  # Random triangle
+             # Degenerated triangles:
+             (0., 0., 0.), (1., 0., 0.), (2., 0., 0.),  # Colinear points
+             (1., 1., 1.), (1., 1., 1.), (1., 1., 1.),  # All same point
+             ),
+            dtype='float32')
+
+        normals = numpy.array(
+            ((0., 0., 1.),
+             (-0.40824829,  0.81649658, -0.40824829),
+             (0., 0., 0.),
+             (0., 0., 0.)),
+            dtype='float32')
+
+        testnormals = utils.trianglesNormal(positions)
+        self.assertTrue(numpy.allclose(testnormals, normals))
+
+
+# suite #######################################################################
+
+def suite():
+    testsuite = unittest.TestSuite()
+    for test in (TestAngleBetweenVectors,
+                 TestPlaneParameters,
+                 TestUnindexArrays,
+                 TestTriangleNormals):
+        testsuite.addTest(
+            unittest.defaultTestLoader.loadTestsFromTestCase(test))
+    return testsuite
+
+
+if __name__ == '__main__':
+    unittest.main(defaultTest='suite')
diff --git a/silx/gui/plot3d/scene/text.py b/silx/gui/plot3d/scene/text.py
new file mode 100644
index 0000000..903fc21
--- /dev/null
+++ b/silx/gui/plot3d/scene/text.py
@@ -0,0 +1,534 @@
+# coding: utf-8
+# /*##########################################################################
+#
+# Copyright (c) 2016-2017 European Synchrotron Radiation Facility
+#
+# Permission is hereby granted, free of charge, to any person obtaining a copy
+# of this software and associated documentation files (the "Software"), to deal
+# in the Software without restriction, including without limitation the rights
+# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+# copies of the Software, and to permit persons to whom the Software is
+# furnished to do so, subject to the following conditions:
+#
+# The above copyright notice and this permission notice shall be included in
+# all copies or substantial portions of the Software.
+#
+# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
+# THE SOFTWARE.
+#
+# ###########################################################################*/
+"""Primitive displaying a text field in the scene."""
+
+from __future__ import absolute_import, division, unicode_literals
+
+__authors__ = ["T. Vincent"]
+__license__ = "MIT"
+__date__ = "17/10/2016"
+
+
+import logging
+import numpy
+
+from silx.gui.plot.Colors import rgba
+
+from ... import _glutils
+from ..._glutils import gl
+
+from ..._glutils import font as _font
+from ...plot._utils import ticklayout
+
+from . import event, primitives, core, transform
+
+
+_logger = logging.getLogger(__name__)
+
+
+class Font(event.Notifier):
+    """Description of a font.
+
+    :param str name: Family of the font
+    :param int size: Size of the font in points
+    :param int weight: Font weight
+    :param bool italic: True for italic font, False (default) otherwise
+    """
+
+    def __init__(self, name=None, size=-1, weight=-1, italic=False):
+        self._name = name if name is not None else _font.getDefaultFontFamily()
+        self._size = size
+        self._weight = weight
+        self._italic = italic
+        super(Font, self).__init__()
+
+    name = event.notifyProperty(
+        '_name',
+        doc="""Name of the font (str)""",
+        converter=str)
+
+    size = event.notifyProperty(
+        '_size',
+        doc="""Font size in points (int)""",
+        converter=int)
+
+    weight = event.notifyProperty(
+        '_weight',
+        doc="""Font size in points (int)""",
+        converter=int)
+
+    italic = event.notifyProperty(
+        '_italic',
+        doc="""True for italic (bool)""",
+        converter=bool)
+
+
+class Text2D(primitives.Geometry):
+    """Text field as a 2D texture displayed with bill-boarding
+
+    :param str text: Text to display
+    :param Font font: The font to use
+    """
+
+    # Text anchor values
+    CENTER = 'center'
+
+    LEFT = 'left'
+    RIGHT = 'right'
+
+    TOP = 'top'
+    BASELINE = 'baseline'
+    BOTTOM = 'bottom'
+
+    _ALIGN = LEFT, CENTER, RIGHT
+    _VALIGN = TOP, BASELINE, CENTER, BOTTOM
+
+    _rasterTextCache = {}
+    """Internal cache storing already rasterized text"""
+    # TODO limit cache size and discard least recent used
+
+    def __init__(self, text='', font=None):
+        self._dirtyTexture = True
+        self._dirtyAlign = True
+        self._baselineOffset = 0
+        self._text = text
+        self._font = font if font is not None else Font()
+        self._foreground = 1., 1., 1., 1.
+        self._background = 0., 0., 0., 0.
+        self._overlay = False
+        self._align = 'left'
+        self._valign = 'baseline'
+        self._devicePixelRatio = 1.0  # Store it to check for changes
+
+        self._texture = None
+        self._textureDirty = True
+
+        super(Text2D, self).__init__(
+            'triangle_strip',
+            copy=False,
+            # Keep an array for position as it is bound to attr 0 and MUST
+            # be active and an array at least on Mac OS X
+            position=numpy.zeros((4, 3), dtype=numpy.float32),
+            vertexID=numpy.arange(4., dtype=numpy.float32).reshape(4, 1),
+            offsetInViewportCoords=(0., 0.))
+
+    @property
+    def text(self):
+        """Text displayed by this primitive (str)"""
+        return self._text
+
+    @text.setter
+    def text(self, text):
+        text = str(text)
+        if text != self._text:
+            self._dirtyTexture = True
+            self._text = text
+            self.notify()
+
+    @property
+    def font(self):
+        """Font to use to raster text (Font)"""
+        return self._font
+
+    @font.setter
+    def font(self, font):
+        self._font.removeListener(self._fontChanged)
+        self._font = font
+        self._font.addListener(self._fontChanged)
+        self._fontChanged(self)  # Which calls notify and primitive as dirty
+
+    def _fontChanged(self, source):
+        """Listen for font change"""
+        self._dirtyTexture = True
+        self.notify()
+
+    foreground = event.notifyProperty(
+        '_foreground', doc="""RGBA color of the text: 4 float in [0, 1]""",
+        converter=rgba)
+
+    background = event.notifyProperty(
+        '_background',
+        doc="RGBA background color of the text field: 4 float in [0, 1]",
+        converter=rgba)
+
+    overlay = event.notifyProperty(
+        '_overlay',
+        doc="True to always display text on top of the scene (default: False)",
+        converter=bool)
+
+    def _setAlign(self, align):
+        assert align in self._ALIGN
+        self._align = align
+        self._dirtyAlign = True
+        self.notify()
+
+    align = property(
+        lambda self: self._align,
+        _setAlign,
+        doc="""Horizontal anchor position of the text field (str).
+
+        Either 'left' (default), 'center' or 'right'.""")
+
+    def _setVAlign(self, valign):
+        assert valign in self._VALIGN
+        self._valign = valign
+        self._dirtyAlign = True
+        self.notify()
+
+    valign = property(
+        lambda self: self._valign,
+        _setVAlign,
+        doc="""Vertical anchor position of the text field (str).
+
+        Either 'top', 'baseline' (default), 'center' or 'bottom'""")
+
+    def _raster(self, devicePixelRatio):
+        """Raster current primitive to a bitmap
+
+        :param float devicePixelRatio:
+            The ratio between device and device-independent pixels
+        :return: Corresponding image in grayscale and baseline offset from top
+        :rtype: (HxW numpy.ndarray of uint8, int)
+        """
+        params = (self.text,
+                  self.font.name,
+                  self.font.size,
+                  self.font.weight,
+                  self.font.italic,
+                  devicePixelRatio)
+
+        if params not in self._rasterTextCache:  # Add to cache
+            self._rasterTextCache[params] = _font.rasterText(*params)
+
+        array, offset = self._rasterTextCache[params]
+        return array.copy(), offset
+
+    def _bounds(self, dataBounds=False):
+        return None
+
+    def prepareGL2(self, context):
+        # Check if devicePixelRatio has changed since last rendering
+        devicePixelRatio = context.glCtx.devicePixelRatio
+        if self._devicePixelRatio != devicePixelRatio:
+            self._devicePixelRatio = devicePixelRatio
+            self._dirtyTexture = True
+
+        if self._dirtyTexture:
+            self._dirtyTexture = False
+
+            if self._texture is not None:
+                self._texture.discard()
+            self._texture = None
+            self._baselineOffset = 0
+
+            if self.text:
+                image, self._baselineOffset = self._raster(
+                    self._devicePixelRatio)
+                self._texture = _glutils.Texture(
+                    gl.GL_R8, image, gl.GL_RED,
+                    minFilter=gl.GL_NEAREST,
+                    magFilter=gl.GL_NEAREST,
+                    wrap=gl.GL_CLAMP_TO_EDGE)
+                self._dirtyAlign = True  # To force update of offset
+
+        if self._dirtyAlign:
+            self._dirtyAlign = False
+
+            if self._texture is not None:
+                height, width = self._texture.shape
+
+                if self._align == 'left':
+                    ox = 0.
+                elif self._align == 'center':
+                    ox = - width // 2
+                elif self._align == 'right':
+                    ox = - width
+                else:
+                    _logger.error("Unsupported align: %s", self._align)
+                    ox = 0.
+
+                if self._valign == 'top':
+                    oy = 0.
+                elif self._valign == 'baseline':
+                    oy = self._baselineOffset
+                elif self._valign == 'center':
+                    oy = height // 2
+                elif self._valign == 'bottom':
+                    oy = height
+                else:
+                    _logger.error("Unsupported valign: %s", self._valign)
+                    oy = 0.
+
+                offsets = (ox, oy) + numpy.array(
+                    ((0., 0.), (width, 0.), (0., -height), (width, -height)),
+                    dtype=numpy.float32)
+                self.setAttribute('offsetInViewportCoords', offsets)
+
+        super(Text2D, self).prepareGL2(context)
+
+    def renderGL2(self, context):
+        if not self.text:
+            return  # Nothing to render
+
+        program = context.glCtx.prog(*self._shaders)
+        program.use()
+
+        program.setUniformMatrix('matrix', context.objectToNDC.matrix)
+        gl.glUniform2f(
+            program.uniforms['viewportSize'], *context.viewport.size)
+        gl.glUniform4f(program.uniforms['foreground'], *self.foreground)
+        gl.glUniform4f(program.uniforms['background'], *self.background)
+        gl.glUniform1i(program.uniforms['texture'], self._texture.texUnit)
+        gl.glUniform1i(program.uniforms['isOverlay'],
+                       1 if self._overlay else 0)
+
+        self._texture.bind()
+
+        if not self._overlay or not gl.glGetBoolean(gl.GL_DEPTH_TEST):
+            self._draw(program)
+        else:  # overlay and depth test currently enabled
+            gl.glDisable(gl.GL_DEPTH_TEST)
+            self._draw(program)
+            gl.glEnable(gl.GL_DEPTH_TEST)
+
+    # TODO texture atlas + viewportSize as attribute to chain text rendering
+
+    _shaders = (
+        """
+    attribute vec3 position;
+    attribute vec2 offsetInViewportCoords;  /* Offset in pixels (y upward) */
+    attribute float vertexID;  /* Index of rectangle corner */
+
+    uniform mat4 matrix;
+    uniform vec2 viewportSize;  /* Width, height of the viewport */
+    uniform int isOverlay;
+
+    varying vec2 texCoords;
+
+    void main(void)
+    {
+        vec4 clipPos = matrix * vec4(position, 1.0);
+        vec4 ndcPos = clipPos / clipPos.w;  /* Perspective divide */
+
+        /* Align ndcPos with pixels in viewport-like coords (origin useless) */
+        vec2 viewportPos = floor((ndcPos.xy + vec2(1.0, 1.0)) * 0.5 * viewportSize);
+
+        /* Apply offset in viewport coords */
+        viewportPos += offsetInViewportCoords;
+
+        /* Convert back to NDC */
+        vec2 pointPos = 2.0 * viewportPos / viewportSize - vec2(1.0, 1.0);
+        float z = (isOverlay != 0) ? -1.0 : ndcPos.z;
+        gl_Position = vec4(pointPos, z, 1.0);
+
+        /* Index : texCoords:
+         * 0: (0., 0.)
+         * 1: (1., 0.)
+         * 2: (0., 1.)
+         * 3: (1., 1.)
+         */
+        texCoords = vec2(vertexID == 0.0 || vertexID == 2.0 ? 0.0 : 1.0,
+                         vertexID < 1.5 ? 0.0 : 1.0);
+    }
+    """,  # noqa
+
+        """
+    varying vec2 texCoords;
+
+    uniform vec4 foreground;
+    uniform vec4 background;
+    uniform sampler2D texture;
+
+    void main(void)
+    {
+        float value = texture2D(texture, texCoords).r;
+
+        if (background.a != 0.0) {
+            gl_FragColor = mix(background, foreground, value);
+        } else {
+            gl_FragColor = foreground;
+            gl_FragColor.a *= value;
+            if (gl_FragColor.a <= 0.01) {
+                discard;
+            }
+        }
+    }
+    """)
+
+
+class LabelledAxes(primitives.GroupBBox):
+    """A group displaying a bounding box with axes labels around its children.
+    """
+
+    def __init__(self):
+        super(LabelledAxes, self).__init__()
+        self._ticksForBounds = None
+
+        self._font = Font()
+
+        # TODO offset labels from anchor in pixels
+
+        self._xlabel = Text2D(font=self._font)
+        self._xlabel.align = 'center'
+        self._xlabel.transforms = [self._boxTransforms,
+                                   transform.Translate(tx=0.5)]
+        self._children.append(self._xlabel)
+
+        self._ylabel = Text2D(font=self._font)
+        self._ylabel.align = 'center'
+        self._ylabel.transforms = [self._boxTransforms,
+                                   transform.Translate(ty=0.5)]
+        self._children.append(self._ylabel)
+
+        self._zlabel = Text2D(font=self._font)
+        self._zlabel.align = 'center'
+        self._zlabel.transforms = [self._boxTransforms,
+                                   transform.Translate(tz=0.5)]
+        self._children.append(self._zlabel)
+
+        self._tickLines = primitives.Lines(  # Init tick lines with dummy pos
+            positions=((0., 0., 0.), (0., 0., 0.)),
+            mode='lines')
+        self._tickLines.visible = False
+        self._children.append(self._tickLines)
+
+        self._tickLabels = core.Group()
+        self._children.append(self._tickLabels)
+
+    @property
+    def font(self):
+        """Font of axes text labels (Font)"""
+        return self._font
+
+    @font.setter
+    def font(self, font):
+        self._font = font
+        self._xlabel.font = font
+        self._ylabel.font = font
+        self._zlabel.font = font
+        for label in self._tickLabels.children:
+            label.font = font
+
+    @property
+    def xlabel(self):
+        """Text label of the X axis (str)"""
+        return self._xlabel.text
+
+    @xlabel.setter
+    def xlabel(self, text):
+        self._xlabel.text = text
+
+    @property
+    def ylabel(self):
+        """Text label of the Y axis (str)"""
+        return self._ylabel.text
+
+    @ylabel.setter
+    def ylabel(self, text):
+        self._ylabel.text = text
+
+    @property
+    def zlabel(self):
+        """Text label of the Z axis (str)"""
+        return self._zlabel.text
+
+    @zlabel.setter
+    def zlabel(self, text):
+        self._zlabel.text = text
+
+    def _updateTicks(self):
+        """Check if ticks need update and update them if needed."""
+        bounds = self._group.bounds(transformed=False, dataBounds=True)
+        if bounds is None:  # No content
+            if self._ticksForBounds is not None:
+                self._ticksForBounds = None
+                self._tickLines.visible = False
+                self._tickLabels.children = []  # Reset previous labels
+
+        elif (self._ticksForBounds is None or
+                not numpy.all(numpy.equal(bounds, self._ticksForBounds))):
+            self._ticksForBounds = bounds
+
+            # Update ticks
+            # TODO make ticks having a constant length on the screen
+            ticklength = numpy.abs(bounds[1] - bounds[0]) / 20.
+
+            xticks, xlabels = ticklayout.ticks(*bounds[:, 0])
+            yticks, ylabels = ticklayout.ticks(*bounds[:, 1])
+            zticks, zlabels = ticklayout.ticks(*bounds[:, 2])
+
+            # Update tick lines
+            coords = numpy.empty(
+                ((len(xticks) + len(yticks) + len(zticks)), 4, 3),
+                dtype=numpy.float32)
+            coords[:, :, :] = bounds[0, :]  # account for offset from origin
+
+            xcoords = coords[:len(xticks)]
+            xcoords[:, :, 0] = numpy.asarray(xticks)[:, numpy.newaxis]
+            xcoords[:, 1, 1] += ticklength[1]  # X ticks on XY plane
+            xcoords[:, 3, 2] += ticklength[2]  # X ticks on XZ plane
+
+            ycoords = coords[len(xticks):len(xticks) + len(yticks)]
+            ycoords[:, :, 1] = numpy.asarray(yticks)[:, numpy.newaxis]
+            ycoords[:, 1, 0] += ticklength[0]  # Y ticks on XY plane
+            ycoords[:, 3, 2] += ticklength[2]  # Y ticks on YZ plane
+
+            zcoords = coords[len(xticks) + len(yticks):]
+            zcoords[:, :, 2] = numpy.asarray(zticks)[:, numpy.newaxis]
+            zcoords[:, 1, 0] += ticklength[0]  # Z ticks on XZ plane
+            zcoords[:, 3, 1] += ticklength[1]  # Z ticks on YZ plane
+
+            self._tickLines.setAttribute('position', coords.reshape(-1, 3))
+            self._tickLines.visible = True
+
+            # Update labels
+            offsets = bounds[0] - ticklength
+            labels = []
+            for tick, label in zip(xticks, xlabels):
+                text = Text2D(text=label, font=self.font)
+                text.align = 'center'
+                text.transforms = [transform.Translate(
+                    tx=tick, ty=offsets[1], tz=offsets[2])]
+                labels.append(text)
+
+            for tick, label in zip(yticks, ylabels):
+                text = Text2D(text=label, font=self.font)
+                text.align = 'center'
+                text.transforms = [transform.Translate(
+                    tx=offsets[0], ty=tick, tz=offsets[2])]
+                labels.append(text)
+
+            for tick, label in zip(zticks, zlabels):
+                text = Text2D(text=label, font=self.font)
+                text.align = 'center'
+                text.transforms = [transform.Translate(
+                    tx=offsets[0], ty=offsets[1], tz=tick)]
+                labels.append(text)
+
+            self._tickLabels.children = labels  # Reset previous labels
+
+    def prepareGL2(self, context):
+        self._updateTicks()
+        super(LabelledAxes, self).prepareGL2(context)
diff --git a/silx/gui/plot3d/scene/transform.py b/silx/gui/plot3d/scene/transform.py
new file mode 100644
index 0000000..71a6b74
--- /dev/null
+++ b/silx/gui/plot3d/scene/transform.py
@@ -0,0 +1,968 @@
+# coding: utf-8
+# /*##########################################################################
+#
+# Copyright (c) 2015-2017 European Synchrotron Radiation Facility
+#
+# Permission is hereby granted, free of charge, to any person obtaining a copy
+# of this software and associated documentation files (the "Software"), to deal
+# in the Software without restriction, including without limitation the rights
+# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+# copies of the Software, and to permit persons to whom the Software is
+# furnished to do so, subject to the following conditions:
+#
+# The above copyright notice and this permission notice shall be included in
+# all copies or substantial portions of the Software.
+#
+# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
+# THE SOFTWARE.
+#
+# ###########################################################################*/
+"""This module provides 4x4 matrix operation and classes to handle them."""
+
+from __future__ import absolute_import, division, unicode_literals
+
+__authors__ = ["T. Vincent"]
+__license__ = "MIT"
+__date__ = "25/07/2016"
+
+
+import itertools
+import numpy
+
+from . import event
+
+
+# Functions ###################################################################
+
+# Projections
+
+def mat4LookAtDir(position, direction, up):
+    """Creates matrix to look in direction from position.
+
+    :param position: Array-like 3 coordinates of the point of view position.
+    :param direction: Array-like 3 coordinates of the sight direction vector.
+    :param up: Array-like 3 coordinates of the upward direction
+               in the image plane.
+    :returns: Corresponding matrix.
+    :rtype: numpy.ndarray of shape (4, 4)
+    """
+    assert len(position) == 3
+    assert len(direction) == 3
+    assert len(up) == 3
+
+    direction = numpy.array(direction, copy=True, dtype=numpy.float32)
+    dirnorm = numpy.linalg.norm(direction)
+    assert dirnorm != 0.
+    direction /= dirnorm
+
+    side = numpy.cross(direction,
+                       numpy.array(up, copy=False, dtype=numpy.float32))
+    sidenorm = numpy.linalg.norm(side)
+    assert sidenorm != 0.
+    up = numpy.cross(side / sidenorm, direction)
+    upnorm = numpy.linalg.norm(up)
+    assert upnorm != 0.
+    up /= upnorm
+
+    matrix = numpy.identity(4, dtype=numpy.float32)
+    matrix[0, :3] = side
+    matrix[1, :3] = up
+    matrix[2, :3] = -direction
+    return numpy.dot(matrix,
+                     mat4Translate(-position[0], -position[1], -position[2]))
+
+
+def mat4LookAt(position, center, up):
+    """Creates matrix to look at center from position.
+
+    See gluLookAt.
+
+    :param position: Array-like 3 coordinates of the point of view position.
+    :param center: Array-like 3 coordinates of the center of the scene.
+    :param up: Array-like 3 coordinates of the upward direction
+               in the image plane.
+    :returns: Corresponding matrix.
+    :rtype: numpy.ndarray of shape (4, 4)
+    """
+    position = numpy.array(position, copy=False, dtype=numpy.float32)
+    center = numpy.array(center, copy=False, dtype=numpy.float32)
+    direction = center - position
+    return mat4LookAtDir(position, direction, up)
+
+
+def mat4Frustum(left, right, bottom, top, near, far):
+    """Creates a frustum projection matrix.
+
+    See glFrustum.
+    """
+    return numpy.array((
+        (2.*near / (right-left), 0., (right+left) / (right-left), 0.),
+        (0., 2.*near / (top-bottom), (top+bottom) / (top-bottom), 0.),
+        (0., 0., -(far+near) / (far-near), -2.*far*near / (far-near)),
+        (0., 0., -1., 0.)), dtype=numpy.float32)
+
+
+def mat4Perspective(fovy, width, height, near, far):
+    """Creates a perspective projection matrix.
+
+    Similar to gluPerspective.
+
+    :param float fovy: Field of view angle in degrees in the y direction.
+    :param float width: Width of the viewport.
+    :param float height: Height of the viewport.
+    :param float near: Distance to the near plane (strictly positive).
+    :param float far: Distance to the far plane (strictly positive).
+    :return: Corresponding matrix.
+    :rtype: numpy.ndarray of shape (4, 4)
+    """
+    assert fovy != 0
+    assert height != 0
+    assert width != 0
+    assert near > 0.
+    assert far > near
+    aspectratio = width / height
+    f = 1. / numpy.tan(numpy.radians(fovy) / 2.)
+    return numpy.array((
+        (f / aspectratio, 0., 0., 0.),
+        (0., f, 0., 0.),
+        (0., 0., (far + near) / (near - far), 2. * far * near / (near - far)),
+        (0., 0., -1., 0.)), dtype=numpy.float32)
+
+
+def mat4Orthographic(left, right, bottom, top, near, far):
+    """Creates an orthographic (i.e., parallel) projection matrix.
+
+    See glOrtho.
+    """
+    return numpy.array((
+        (2. / (right - left), 0., 0., - (right + left) / (right - left)),
+        (0., 2. / (top - bottom), 0., - (top + bottom) / (top - bottom)),
+        (0., 0., -2. / (far - near), - (far + near) / (far - near)),
+        (0., 0., 0., 1.)), dtype=numpy.float32)
+
+
+# Affine
+
+def mat4Translate(tx, ty, tz):
+    """4x4 translation matrix."""
+    return numpy.array((
+        (1., 0., 0., tx),
+        (0., 1., 0., ty),
+        (0., 0., 1., tz),
+        (0., 0., 0., 1.)), dtype=numpy.float32)
+
+
+def mat4Scale(sx, sy, sz):
+    """4x4 scale matrix."""
+    return numpy.array((
+        (sx, 0., 0., 0.),
+        (0., sy, 0., 0.),
+        (0., 0., sz, 0.),
+        (0., 0., 0., 1.)), dtype=numpy.float32)
+
+
+def mat4RotateFromAngleAxis(angle, x=0., y=0., z=1.):
+    """4x4 rotation matrix from angle and axis.
+
+    :param float angle: The rotation angle in radians.
+    :param float x: The rotation vector x coordinate.
+    :param float y: The rotation vector y coordinate.
+    :param float z: The rotation vector z coordinate.
+    """
+    ca = numpy.cos(angle)
+    sa = numpy.sin(angle)
+    return numpy.array((
+        ((1.-ca) * x*x + ca,   (1.-ca) * x*y - sa*z, (1.-ca) * x*z + sa*y, 0.),
+        ((1.-ca) * x*y + sa*z, (1.-ca) * y*y + ca,   (1.-ca) * y*z - sa*x, 0.),
+        ((1.-ca) * x*z - sa*y, (1.-ca) * y*z + sa*x, (1.-ca) * z*z + ca, 0.),
+        (0., 0., 0., 1.)), dtype=numpy.float32)
+
+
+def mat4RotateFromQuaternion(quaternion):
+    """4x4 rotation matrix from quaternion.
+
+    :param quaternion: Array-like unit quaternion stored as (x, y, z, w)
+    """
+    quaternion = numpy.array(quaternion, copy=True)
+    quaternion /= numpy.linalg.norm(quaternion)
+
+    qx, qy, qz, qw = quaternion
+    return numpy.array((
+        (1. - 2.*(qy**2 + qz**2), 2.*(qx*qy - qw*qz), 2.*(qx*qz + qw*qy), 0.),
+        (2.*(qx*qy + qw*qz), 1. - 2.*(qx**2 + qz**2), 2.*(qy*qz - qw*qx), 0.),
+        (2.*(qx*qz - qw*qy), 2.*(qy*qz + qw*qx), 1. - 2.*(qx**2 + qy**2), 0.),
+        (0., 0., 0., 1.)), dtype=numpy.float32)
+
+
+def mat4Shear(axis, sx=0., sy=0., sz=0.):
+    """4x4 shear matrix: Skew two axes relative to a third fixed one.
+
+    shearFactor = tan(shearAngle)
+
+    :param str axis: The axis to keep constant and shear against.
+                     In 'x', 'y', 'z'.
+    :param float sx: The shear factor for the X axis relative to axis.
+    :param float sy: The shear factor for the Y axis relative to axis.
+    :param float sz: The shear factor for the Z axis relative to axis.
+    """
+    assert axis in ('x', 'y', 'z')
+
+    matrix = numpy.identity(4, dtype=numpy.float32)
+
+    # Make the shear column
+    index = 'xyz'.find(axis)
+    shearcolumn = numpy.array((sx, sy, sz, 0.), dtype=numpy.float32)
+    shearcolumn[index] = 1.
+    matrix[:, index] = shearcolumn
+    return matrix
+
+
+# Transforms ##################################################################
+
+class Transform(event.Notifier):
+
+    def __init__(self, static=False):
+        """Base class for (row-major) 4x4 matrix transforms.
+
+        :param bool static: False (default) to reset cache when changed,
+                            True for static matrices.
+        """
+        super(Transform, self).__init__()
+        self._matrix = None
+        self._inverse = None
+        if not static:
+            self.addListener(self._changed)  # Listening self for changes
+
+    def __repr__(self):
+        return '%s(%s)' % (self.__class__.__init__,
+                           repr(self.getMatrix(copy=False)))
+
+    def inverse(self):
+        """Return the Transform of the inverse.
+
+        The returned Transform is static, it is not updated when this
+        Transform is modified.
+
+        :return: A Transform which is the inverse of this Transform.
+        """
+        return Inverse(self)
+
+    # Matrix
+
+    def _makeMatrix(self):
+        """Override to build matrix"""
+        return numpy.identity(4, dtype=numpy.float32)
+
+    def _makeInverse(self):
+        """Override to build inverse matrix."""
+        return numpy.linalg.inv(self.getMatrix(copy=False))
+
+    def getMatrix(self, copy=True):
+        """The 4x4 matrix of this transform.
+
+        :param bool copy: True (the default) to get a copy of the matrix,
+                          False to get the internal matrix, do not modify!
+        :return: 4x4 matrix of this transform.
+        """
+        if self._matrix is None:
+            self._matrix = self._makeMatrix()
+        if copy:
+            return self._matrix.copy()
+        else:
+            return self._matrix
+
+    matrix = property(getMatrix, doc="The 4x4 matrix of this transform.")
+
+    def getInverseMatrix(self, copy=False):
+        """The 4x4 matrix of the inverse of this transform.
+
+        :param bool copy: True (the default) to get a copy of the matrix,
+                          False to get the internal matrix, do not modify!
+        :return: 4x4 matrix of the inverse of this transform.
+        """
+        if self._inverse is None:
+            self._inverse = self._makeInverse()
+        if copy:
+            return self._inverse.copy()
+        else:
+            return self._inverse
+
+    inverseMatrix = property(
+        getInverseMatrix,
+        doc="The 4x4 matrix of the inverse of this transform.")
+
+    # Listener
+
+    def _changed(self, source):
+        """Default self listener reseting matrix cache."""
+        self._matrix = None
+        self._inverse = None
+
+    # Multiplication with vectors
+
+    @staticmethod
+    def _prepareVector(vector, w):
+        """Add 4th coordinate (w) to vector if missing."""
+        assert len(vector) in (3, 4)
+        vector = numpy.array(vector, copy=False, dtype=numpy.float32)
+        if len(vector) == 3:
+            vector = numpy.append(vector, w)
+        return vector
+
+    def transformPoint(self, point, direct=True, perspectiveDivide=False):
+        """Apply the transform to a point.
+
+        If len(point) == 3, apply persective divide if possible.
+
+        :param point: Array-like vector of 3 or 4 coordinates.
+        :param bool direct: Whether to apply the direct (True, the default)
+                            or inverse (False) transform.
+        :param bool perspectiveDivide: Whether to apply the perspective divide
+                                       (True) or not (False, the default).
+        :return: The transformed point.
+        :rtype: numpy.ndarray of same length as point.
+        """
+        if direct:
+            matrix = self.getMatrix(copy=False)
+        else:
+            matrix = self.getInverseMatrix(copy=False)
+        result = numpy.dot(matrix, self._prepareVector(point, 1.))
+
+        if perspectiveDivide and result[3] != 0.:
+            result /= result[3]
+
+        if len(point) == 3:
+            return result[:3]
+        else:
+            return result
+
+    def transformDir(self, direction, direct=True):
+        """Apply the transform to a direction.
+
+        :param direction: Array-like vector of 3 coordinates.
+        :param bool direct: Whether to apply the direct (True, the default)
+                            or inverse (False) transform.
+        :return: The transformed direction.
+        :rtype: numpy.ndarray of length 3.
+        """
+        if direct:
+            matrix = self.getMatrix(copy=False)
+        else:
+            matrix = self.getInverseMatrix(copy=False)
+        return numpy.dot(matrix[:3, :3], direction[:3])
+
+    def transformNormal(self, normal, direct=True):
+        """Apply the transform to a normal: R = (M-1)t * V.
+
+        :param normal: Array-like vector of 3 coordinates.
+        :param bool direct: Whether to apply the direct (True, the default)
+                            or inverse (False) transform.
+        :return: The transformed normal.
+        :rtype: numpy.ndarray of length 3.
+        """
+        if direct:
+            matrix = self.getInverseMatrix(copy=False).T
+        else:
+            matrix = self.getMatrix(copy=False).T
+        return numpy.dot(matrix[:3, :3], normal[:3])
+
+    _CUBE_CORNERS = numpy.array(list(itertools.product((0., 1.), repeat=3)),
+                                dtype=numpy.float32)
+    """Unit cube corners used by :meth:`transformRectangularBox`"""
+
+    def transformBounds(self, bounds, direct=True):
+        """Apply the transform to an axes-aligned rectangular box.
+
+        :param bounds: Min and max coords of the box for each axes.
+        :type bounds: 2x3 numpy.ndarray
+        :param bool direct: Whether to apply the direct (True, the default)
+                            or inverse (False) transform.
+        :return: Axes-aligned rectangular box including the transformed box.
+        :rtype: 2x3 numpy.ndarray of float32
+        """
+        corners = numpy.ones((8, 4), dtype=numpy.float32)
+        corners[:, :3] = bounds[0] + \
+            self._CUBE_CORNERS * (bounds[1] - bounds[0])
+
+        if direct:
+            matrix = self.getMatrix(copy=False)
+        else:
+            matrix = self.getInverseMatrix(copy=False)
+
+        # Transform corners
+        cornerstransposed = numpy.dot(matrix, corners.T)
+        cornerstransposed = cornerstransposed / cornerstransposed[3]
+
+        # Get min/max for each axis
+        transformedbounds = numpy.empty((2, 3), dtype=numpy.float32)
+        transformedbounds[0] = cornerstransposed.T[:, :3].min(axis=0)
+        transformedbounds[1] = cornerstransposed.T[:, :3].max(axis=0)
+
+        return transformedbounds
+
+
+class Inverse(Transform):
+    """Transform which is the inverse of another one.
+
+    Static: It never gets updated.
+    """
+
+    def __init__(self, transform):
+        """Initializer.
+
+        :param Transform transform: The transform to invert.
+        """
+
+        super(Inverse, self).__init__(static=True)
+        self._matrix = transform.getInverseMatrix(copy=True)
+        self._inverse = transform.getMatrix(copy=True)
+
+
+class TransformList(Transform, event.HookList):
+    """List of transforms."""
+
+    def __init__(self, iterable=()):
+        Transform.__init__(self)
+        event.HookList.__init__(self, iterable)
+
+    def _listWillChangeHook(self, methodName, *args, **kwargs):
+        for item in self:
+            item.removeListener(self._transformChanged)
+
+    def _listWasChangedHook(self, methodName, *args, **kwargs):
+        for item in self:
+            item.addListener(self._transformChanged)
+        self.notify()
+
+    def _transformChanged(self, source):
+        """Listen to transform changes of the list and its items."""
+        if source is not self:  # Avoid infinite recursion
+            self.notify()
+
+    def _makeMatrix(self):
+        matrix = numpy.identity(4, dtype=numpy.float32)
+        for transform in self:
+            matrix = numpy.dot(matrix, transform.getMatrix(copy=False))
+        return matrix
+
+
+class StaticTransformList(Transform):
+    """Transform that is a snapshot of a list of Transforms
+
+    It does not keep reference to the list of Transforms.
+
+    :param iterable: Iterable of Transform used for initialization
+    """
+
+    def __init__(self, iterable=()):
+        super(StaticTransformList, self).__init__(static=True)
+        matrix = numpy.identity(4, dtype=numpy.float32)
+        for transform in iterable:
+            matrix = numpy.dot(matrix, transform.getMatrix(copy=False))
+        self._matrix = matrix  # Init matrix once
+
+
+# Affine ######################################################################
+
+class Matrix(Transform):
+
+    def __init__(self, matrix=None):
+        """4x4 Matrix.
+
+        :param matrix: 4x4 array-like matrix or None for identity matrix.
+        """
+        super(Matrix, self).__init__(static=True)
+        self.setMatrix(matrix)
+
+    def setMatrix(self, matrix=None):
+        """Update the 4x4 Matrix.
+
+        :param matrix: 4x4 array-like matrix or None for identity matrix.
+        """
+        if matrix is None:
+            self._matrix = numpy.identity(4, dtype=numpy.float32)
+        else:
+            matrix = numpy.array(matrix, copy=True, dtype=numpy.float32)
+            assert matrix.shape == (4, 4)
+            self._matrix = matrix
+        # Reset cached inverse as Transform is declared static
+        self._inverse = None
+        self.notify()
+
+    # Redefined here to add a setter
+    matrix = property(Transform.getMatrix, setMatrix,
+                      doc="The 4x4 matrix of this transform.")
+
+
+class Translate(Transform):
+    """4x4 translation matrix."""
+
+    def __init__(self, tx=0., ty=0., tz=0.):
+        super(Translate, self).__init__()
+        self._tx, self._ty, self._tz = 0., 0., 0.
+        self.setTranslate(tx, ty, tz)
+
+    def _makeMatrix(self):
+        return mat4Translate(self.tx, self.ty, self.tz)
+
+    def _makeInverse(self):
+        return mat4Translate(-self.tx, -self.ty, -self.tz)
+
+    @property
+    def tx(self):
+        return self._tx
+
+    @tx.setter
+    def tx(self, tx):
+        self.setTranslate(tx=tx)
+
+    @property
+    def ty(self):
+        return self._ty
+
+    @ty.setter
+    def ty(self, ty):
+        self.setTranslate(ty=ty)
+
+    @property
+    def tz(self):
+        return self._tz
+
+    @tz.setter
+    def tz(self, tz):
+        self.setTranslate(tz=tz)
+
+    @property
+    def translation(self):
+        return numpy.array((self.tx, self.ty, self.tz), dtype=numpy.float32)
+
+    @translation.setter
+    def translation(self, translations):
+        tx, ty, tz = translations
+        self.setTranslate(tx, ty, tz)
+
+    def setTranslate(self, tx=None, ty=None, tz=None):
+        if tx is not None:
+            self._tx = tx
+        if ty is not None:
+            self._ty = ty
+        if tz is not None:
+            self._tz = tz
+        self.notify()
+
+
+class Scale(Transform):
+    """4x4 scale matrix."""
+
+    def __init__(self, sx=1., sy=1., sz=1.):
+        super(Scale, self).__init__()
+        self._sx, self._sy, self._sz = 0., 0., 0.
+        self.setScale(sx, sy, sz)
+
+    def _makeMatrix(self):
+        return mat4Scale(self.sx, self.sy, self.sz)
+
+    def _makeInverse(self):
+        return mat4Scale(1. / self.sx, 1. / self.sy, 1. / self.sz)
+
+    @property
+    def sx(self):
+        return self._sx
+
+    @sx.setter
+    def sx(self, sx):
+        self.setScale(sx=sx)
+
+    @property
+    def sy(self):
+        return self._sy
+
+    @sy.setter
+    def sy(self, sy):
+        self.setScale(sy=sy)
+
+    @property
+    def sz(self):
+        return self._sz
+
+    @sz.setter
+    def sz(self, sz):
+        self.setScale(sz=sz)
+
+    @property
+    def scale(self):
+        return numpy.array((self._sx, self._sy, self._sz), dtype=numpy.float32)
+
+    @scale.setter
+    def scale(self, scales):
+        sx, sy, sz = scales
+        self.setScale(sx, sy, sz)
+
+    def setScale(self, sx=None, sy=None, sz=None):
+        if sx is not None:
+            assert sx != 0.
+            self._sx = sx
+        if sy is not None:
+            assert sy != 0.
+            self._sy = sy
+        if sz is not None:
+            assert sz != 0.
+            self._sz = sz
+        self.notify()
+
+
+class Rotate(Transform):
+
+    def __init__(self, angle=0., ax=0., ay=0., az=1.):
+        """4x4 rotation matrix.
+
+        :param float angle: The rotation angle in degrees.
+        :param float ax: The x coordinate of the rotation axis.
+        :param float ay: The y coordinate of the rotation axis.
+        :param float az: The z coordinate of the rotation axis.
+        """
+        super(Rotate, self).__init__()
+        self._angle = 0.
+        self._axis = None
+        self.setAngleAxis(angle, (ax, ay, az))
+
+    @property
+    def angle(self):
+        """The rotation angle in degrees."""
+        return self._angle
+
+    @angle.setter
+    def angle(self, angle):
+        self.setAngleAxis(angle=angle)
+
+    @property
+    def axis(self):
+        """The normalized rotation axis as a numpy.ndarray."""
+        return self._axis.copy()
+
+    @axis.setter
+    def axis(self, axis):
+        self.setAngleAxis(axis=axis)
+
+    def setAngleAxis(self, angle=None, axis=None):
+        """Update the angle and/or axis of the rotation.
+
+        :param float angle: The rotation angle in degrees.
+        :param axis: Array-like axis vector (3 coordinates).
+        """
+        if angle is not None:
+            self._angle = angle
+        if axis is not None:
+            assert len(axis) == 3
+            axis = numpy.array(axis, copy=True, dtype=numpy.float32)
+            assert axis.size == 3
+            norm = numpy.linalg.norm(axis)
+            if norm == 0.:  # No axis, set rotation angle to 0.
+                self._angle = 0.
+                self._axis = numpy.array((0., 0., 1.), dtype=numpy.float32)
+            else:
+                self._axis = axis / norm
+
+        if angle is not None or axis is not None:
+            self.notify()
+
+    @property
+    def quaternion(self):
+        """Rotation unit quaternion as (x, y, z, w).
+
+        Where: ||(x, y, z)|| = sin(angle/2),  w = cos(angle/2).
+        """
+        if numpy.linalg.norm(self._axis) == 0.:
+            return numpy.array((0., 0., 0., 1.), dtype=numpy.float32)
+
+        else:
+            quaternion = numpy.empty((4,), dtype=numpy.float32)
+            halfangle = 0.5 * numpy.radians(self.angle)
+            quaternion[0:3] = numpy.sin(halfangle) * self._axis
+            quaternion[3] = numpy.cos(halfangle)
+            return quaternion
+
+    @quaternion.setter
+    def quaternion(self, quaternion):
+        assert len(quaternion) == 4
+
+        # Normalize quaternion
+        quaternion = numpy.array(quaternion, copy=True)
+        quaternion /= numpy.linalg.norm(quaternion)
+
+        # Get angle
+        sinhalfangle = numpy.linalg.norm(quaternion[0:3])
+        coshalfangle = quaternion[3]
+        angle = 2. * numpy.arctan2(sinhalfangle, coshalfangle)
+
+        # Axis will be normalized in setAngleAxis
+        self.setAngleAxis(numpy.degrees(angle), quaternion[0:3])
+
+    def _makeMatrix(self):
+        angle = numpy.radians(self.angle, dtype=numpy.float32)
+        return mat4RotateFromAngleAxis(angle, *self.axis)
+
+    def _makeInverse(self):
+        return numpy.array(self.getMatrix(copy=False).transpose(),
+                           copy=True, order='C',
+                           dtype=numpy.float32)
+
+
+class Shear(Transform):
+
+    def __init__(self, axis, sx=0., sy=0., sz=0.):
+        """4x4 shear/skew matrix of 2 axes relative to the third one.
+
+        :param str axis: The axis to keep fixed, in 'x', 'y', 'z'
+        :param float sx: The shear factor for the x axis.
+        :param float sy: The shear factor for the y axis.
+        :param float sz: The shear factor for the z axis.
+        """
+        assert axis in ('x', 'y', 'z')
+        super(Shear, self).__init__()
+        self._axis = axis
+        self._factors = sx, sy, sz
+
+    @property
+    def axis(self):
+        """The axis against which other axes are skewed."""
+        return self._axis
+
+    @property
+    def factors(self):
+        """The shear factors: shearFactor = tan(shearAngle)"""
+        return self._factors
+
+    def _makeMatrix(self):
+        return mat4Shear(self.axis, *self.factors)
+
+    def _makeInverse(self):
+        sx, sy, sz = self.factors
+        return mat4Shear(self.axis, -sx, -sy, -sz)
+
+
+# Projection ##################################################################
+
+class _Projection(Transform):
+    """Base class for projection matrix.
+
+    Handles near and far clipping plane values.
+    Subclasses must implement :meth:`_makeMatrix`.
+
+    :param float near: Distance to the near plane.
+    :param float far: Distance to the far plane.
+    :param bool checkDepthExtent: Toggle checks near > 0 and far > near.
+    :param size: Viewport's size used to compute the aspect ratio.
+    :type size: 2-tuple of float (width, height).
+    """
+
+    def __init__(self, near, far, checkDepthExtent=False, size=(1., 1.)):
+        super(_Projection, self).__init__()
+        self._checkDepthExtent = checkDepthExtent
+        self._depthExtent = 1, 10
+        self.setDepthExtent(near, far)  # set _depthExtent
+        self._size = 1., 1.
+        self.size = size  # set _size
+
+    def setDepthExtent(self, near=None, far=None):
+        """Set the extent of the visible area along the viewing direction.
+
+        :param float near: The near clipping plane Z coord.
+        :param float far: The far clipping plane Z coord.
+        """
+        near = float(near) if near is not None else self._depthExtent[0]
+        far = float(far) if far is not None else self._depthExtent[1]
+
+        if self._checkDepthExtent:
+            assert near > 0.
+            assert far > near
+
+        self._depthExtent = near, far
+        self.notify()
+
+    @property
+    def near(self):
+        """Distance to the near plane."""
+        return self._depthExtent[0]
+
+    @near.setter
+    def near(self, near):
+        if near != self.near:
+            self.setDepthExtent(near=near)
+
+    @property
+    def far(self):
+        """Distance to the far plane."""
+        return self._depthExtent[1]
+
+    @far.setter
+    def far(self, far):
+        if far != self.far:
+            self.setDepthExtent(far=far)
+
+    @property
+    def size(self):
+        """Viewport size as a 2-tuple of float (width, height)."""
+        return self._size
+
+    @size.setter
+    def size(self, size):
+        assert len(size) == 2
+        self._size = tuple(size)
+        self.notify()
+
+
+class Orthographic(_Projection):
+    """Orthographic (i.e., parallel) projection which keeps aspect ratio.
+
+    Clipping planes are adjusted to match the aspect ratio of
+    the :attr:`size` attribute.
+
+    The left, right, bottom and top parameters defines the area which must
+    always remain visible.
+    Effective clipping planes are adjusted to keep the aspect ratio.
+
+    :param float left: Coord of the left clipping plane.
+    :param float right: Coord of the right clipping plane.
+    :param float bottom: Coord of the bottom clipping plane.
+    :param float top: Coord of the top clipping plane.
+    :param float near: Distance to the near plane.
+    :param float far: Distance to the far plane.
+    :param size: Viewport's size used to compute the aspect ratio.
+    :type size: 2-tuple of float (width, height).
+    """
+
+    def __init__(self, left=0., right=1., bottom=1., top=0., near=-1., far=1.,
+                 size=(1., 1.)):
+        self._left, self._right = left, right
+        self._bottom, self._top = bottom, top
+        super(Orthographic, self).__init__(near, far, checkDepthExtent=False,
+                                           size=size)
+        # _update called when setting size
+
+    def _makeMatrix(self):
+        return mat4Orthographic(
+            self.left, self.right, self.bottom, self.top, self.near, self.far)
+
+    def _update(self, left, right, bottom, top):
+        width, height = self.size
+        aspect = width / height
+
+        orthoaspect = abs(left - right) / abs(bottom - top)
+
+        if orthoaspect >= aspect:  # Keep width, enlarge height
+            newheight = \
+                numpy.sign(top - bottom) * abs(left - right) / aspect
+            bottom = 0.5 * (bottom + top) - 0.5 * newheight
+            top = bottom + newheight
+
+        else:  # Keep height, enlarge width
+            newwidth = \
+                numpy.sign(right - left) * abs(bottom - top) * aspect
+            left = 0.5 * (left + right) - 0.5 * newwidth
+            right = left + newwidth
+
+        # Store values
+        self._left, self._right = left, right
+        self._bottom, self._top = bottom, top
+
+    def setClipping(self, left=None, right=None, bottom=None, top=None):
+        """Set the clipping planes of the projection.
+
+        Parameters are adjusted to keep aspect ratio.
+        If a clipping plane coord is not provided, it uses its current value
+
+        :param float left: Coord of the left clipping plane.
+        :param float right: Coord of the right clipping plane.
+        :param float bottom: Coord of the bottom clipping plane.
+        :param float top: Coord of the top clipping plane.
+        """
+        left = float(left) if left is not None else self.left
+        right = float(right) if right is not None else self.right
+        bottom = float(bottom) if bottom is not None else self.bottom
+        top = float(top) if top is not None else self.top
+
+        self._update(left, right, bottom, top)
+        self.notify()
+
+    left = property(lambda self: self._left,
+                    doc="Coord of the left clipping plane.")
+
+    right = property(lambda self: self._right,
+                     doc="Coord of the right clipping plane.")
+
+    bottom = property(lambda self: self._bottom,
+                      doc="Coord of the bottom clipping plane.")
+
+    top = property(lambda self: self._top,
+                   doc="Coord of the top clipping plane.")
+
+    @property
+    def size(self):
+        """Viewport size as a 2-tuple of float (width, height) or None."""
+        return self._size
+
+    @size.setter
+    def size(self, size):
+        assert len(size) == 2
+        self._size = float(size[0]), float(size[1])
+        self._update(self.left, self.right, self.bottom, self.top)
+        self.notify()
+
+
+class Ortho2DWidget(_Projection):
+    """Orthographic projection with pixel as unit.
+
+    Provides same coordinates as widgets:
+    origin: top left, X axis goes left, Y axis goes down.
+
+    :param float near: Z coordinate of the near clipping plane.
+    :param float far: Z coordinante of the far clipping plane.
+    :param size: Viewport's size used to compute the aspect ratio.
+    :type size: 2-tuple of float (width, height).
+    """
+
+    def __init__(self, near=-1., far=1., size=(1., 1.)):
+
+        super(Ortho2DWidget, self).__init__(near, far, size)
+
+    def _makeMatrix(self):
+        width, height = self.size
+        return mat4Orthographic(0., width, height, 0., self.near, self.far)
+
+
+class Perspective(_Projection):
+    """Perspective projection matrix defined by FOV and aspect ratio.
+
+    :param float fovy: Vertical field-of-view in degrees.
+    :param float near: The near clipping plane Z coord (stricly positive).
+    :param float far: The far clipping plane Z coord (> near).
+    :param size: Viewport's size used to compute the aspect ratio.
+    :type size: 2-tuple of float (width, height).
+    """
+
+    def __init__(self, fovy=90., near=0.1, far=1., size=(1., 1.)):
+
+        super(Perspective, self).__init__(near, far, checkDepthExtent=True)
+        self._fovy = 90.
+        self.fovy = fovy  # Set _fovy
+        self.size = size  # Set _ size
+
+    def _makeMatrix(self):
+        width, height = self.size
+        return mat4Perspective(self.fovy, width, height, self.near, self.far)
+
+    @property
+    def fovy(self):
+        """Vertical field-of-view in degrees."""
+        return self._fovy
+
+    @fovy.setter
+    def fovy(self, fovy):
+        self._fovy = float(fovy)
+        self.notify()
diff --git a/silx/gui/plot3d/scene/utils.py b/silx/gui/plot3d/scene/utils.py
new file mode 100644
index 0000000..930a087
--- /dev/null
+++ b/silx/gui/plot3d/scene/utils.py
@@ -0,0 +1,516 @@
+# coding: utf-8
+# /*##########################################################################
+#
+# Copyright (c) 2015-2017 European Synchrotron Radiation Facility
+#
+# Permission is hereby granted, free of charge, to any person obtaining a copy
+# of this software and associated documentation files (the "Software"), to deal
+# in the Software without restriction, including without limitation the rights
+# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+# copies of the Software, and to permit persons to whom the Software is
+# furnished to do so, subject to the following conditions:
+#
+# The above copyright notice and this permission notice shall be included in
+# all copies or substantial portions of the Software.
+#
+# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
+# THE SOFTWARE.
+#
+# ###########################################################################*/
+"""
+This module provides functions to generate indices, to check intersection
+and to handle planes.
+"""
+
+from __future__ import absolute_import, division, unicode_literals
+
+__authors__ = ["T. Vincent"]
+__license__ = "MIT"
+__date__ = "25/07/2016"
+
+
+import logging
+import numpy
+
+from . import event
+
+
+_logger = logging.getLogger(__name__)
+
+
+# numpy #######################################################################
+
+def _uniqueAlongLastAxis(a):
+    """Numpy unique on the last axis of a 2D array
+
+    Implemented here as not in numpy as of writing.
+
+    See adding axis parameter to numpy.unique:
+    https://github.com/numpy/numpy/pull/3584/files#r6225452
+
+    :param array_like a: Input array.
+    :return: Unique elements along the last axis.
+    :rtype: numpy.ndarray
+    """
+    assert len(a.shape) == 2
+
+    # Construct a type over last array dimension to run unique on a 1D array
+    if a.dtype.char in numpy.typecodes['AllInteger']:
+        # Bit-wise comparison of the 2 indices of a line at once
+        # Expect a C contiguous array of shape N, 2
+        uniquedt = numpy.dtype((numpy.void, a.itemsize * a.shape[-1]))
+    elif a.dtype.char in numpy.typecodes['Float']:
+        uniquedt = [('f{i}'.format(i=i), a.dtype) for i in range(a.shape[-1])]
+    else:
+        raise TypeError("Unsupported type {dtype}".format(dtype=a.dtype))
+
+    uniquearray = numpy.unique(numpy.ascontiguousarray(a).view(uniquedt))
+    return uniquearray.view(a.dtype).reshape((-1, a.shape[-1]))
+
+
+# conversions #################################################################
+
+def triangleToLineIndices(triangleIndices, unicity=False):
+    """Generates lines indices from triangle indices.
+
+    This is generating lines indices for the edges of the triangles.
+
+    :param triangleIndices: The indices to draw a set of vertices as triangles.
+    :type triangleIndices: numpy.ndarray
+    :param bool unicity: If True remove duplicated lines,
+                         else (the default) returns all lines.
+    :return: The indices to draw the edges of the triangles as lines.
+    :rtype: 1D numpy.ndarray of uint16 or uint32.
+    """
+    # Makes sure indices ar packed by triangle
+    triangleIndices = triangleIndices.reshape(-1, 3)
+
+    # Pack line indices by triangle and by edge
+    lineindices = numpy.empty((len(triangleIndices), 3, 2),
+                              dtype=triangleIndices.dtype)
+    lineindices[:, 0] = triangleIndices[:, :2]  # edge = t0, t1
+    lineindices[:, 1] = triangleIndices[:, 1:]  # edge =t1, t2
+    lineindices[:, 2] = triangleIndices[:, ::2]  # edge = t0, t2
+
+    if unicity:
+        lineindices = _uniqueAlongLastAxis(lineindices.reshape(-1, 2))
+
+    # Make sure it is 1D
+    lineindices.shape = -1
+
+    return lineindices
+
+
+def verticesNormalsToLines(vertices, normals, scale=1.):
+    """Return vertices of lines representing normals at given positions.
+
+    :param vertices: Positions of the points.
+    :type vertices: numpy.ndarray with shape: (nbPoints, 3)
+    :param normals: Corresponding normals at the points.
+    :type normals: numpy.ndarray with shape: (nbPoints, 3)
+    :param float scale: The scale factor to apply to normals.
+    :returns: Array of vertices to draw corresponding lines.
+    :rtype: numpy.ndarray with shape: (nbPoints * 2, 3)
+    """
+    linevertices = numpy.empty((len(vertices) * 2, 3), dtype=vertices.dtype)
+    linevertices[0::2] = vertices
+    linevertices[1::2] = vertices + scale * normals
+    return linevertices
+
+
+def unindexArrays(mode, indices, *arrays):
+    """Convert indexed GL primitives to unindexed ones.
+
+    Given indices in arrays and the OpenGL primitive they represent,
+    return the unindexed equivalent.
+
+    :param str mode:
+       Kind of primitive represented by indices.
+       In: points, lines, line_strip, loop, triangles, triangle_strip, fan.
+    :param indices: Indices in other arrays
+    :type indices: numpy.ndarray of dimension 1.
+    :param arrays: Remaining arguments are arrays to convert
+    :return: Converted arrays
+    :rtype: tuple of numpy.ndarray
+    """
+    indices = numpy.array(indices, copy=False)
+
+    assert mode in ('points',
+                    'lines', 'line_strip', 'loop',
+                    'triangles', 'triangle_strip', 'fan')
+
+    if mode in ('lines', 'line_strip', 'loop'):
+        assert len(indices) >= 2
+    elif mode in ('triangles', 'triangle_strip', 'fan'):
+        assert len(indices) >= 3
+
+    assert indices.min() >= 0
+    max_index = indices.max()
+    for data in arrays:
+        assert len(data) >= max_index
+
+    if mode == 'line_strip':
+        unpacked = numpy.empty((2 * (len(indices) - 1),), dtype=indices.dtype)
+        unpacked[0::2] = indices[:-1]
+        unpacked[1::2] = indices[1:]
+        indices = unpacked
+
+    elif mode == 'loop':
+        unpacked = numpy.empty((2 * len(indices),), dtype=indices.dtype)
+        unpacked[0::2] = indices
+        unpacked[1:-1:2] = indices[1:]
+        unpacked[-1] = indices[0]
+        indices = unpacked
+
+    elif mode == 'triangle_strip':
+        unpacked = numpy.empty((3 * (len(indices) - 2),), dtype=indices.dtype)
+        unpacked[0::3] = indices[:-2]
+        unpacked[1::3] = indices[1:-1]
+        unpacked[2::3] = indices[2:]
+        indices = unpacked
+
+    elif mode == 'fan':
+        unpacked = numpy.empty((3 * (len(indices) - 2),), dtype=indices.dtype)
+        unpacked[0::3] = indices[0]
+        unpacked[1::3] = indices[1:-1]
+        unpacked[2::3] = indices[2:]
+        indices = unpacked
+
+    return tuple(numpy.ascontiguousarray(data[indices]) for data in arrays)
+
+
+def trianglesNormal(positions):
+    """Return normal for each triangle.
+
+    :param positions: Serie of triangle's corners
+    :type positions: numpy.ndarray of shape (NbTriangles*3, 3)
+    :return: Normals corresponding to each position.
+    :rtype: numpy.ndarray of shape (NbTriangles, 3)
+    """
+    assert positions.ndim == 2
+    assert positions.shape[1] == 3
+
+    positions = numpy.array(positions, copy=False).reshape(-1, 3, 3)
+
+    normals = numpy.cross(positions[:, 1] - positions[:, 0],
+                          positions[:, 2] - positions[:, 0])
+
+    # Normalize normals
+    if numpy.version.version < '1.8.0':
+        # Debian 7 support: numpy.linalg.norm has no axis argument
+        norms = numpy.array(tuple(numpy.linalg.norm(vec) for vec in normals),
+                            dtype=normals.dtype)
+    else:
+        norms = numpy.linalg.norm(normals, axis=1)
+    norms[norms == 0] = 1
+
+    return normals / norms.reshape(-1, 1)
+
+
+# grid ########################################################################
+
+def gridVertices(dim0Array, dim1Array, dtype):
+    """Generate an array of 2D positions from 2 arrays of 1D coordinates.
+
+    :param dim0Array: 1D array-like of coordinates along the first dimension.
+    :param dim1Array: 1D array-like of coordinates along the second dimension.
+    :param numpy.dtype dtype: Data type of the output array.
+    :return: Array of grid coordinates.
+    :rtype: numpy.ndarray with shape: (len(dim0Array), len(dim1Array), 2)
+    """
+    grid = numpy.empty((len(dim0Array), len(dim1Array), 2), dtype=dtype)
+    grid.T[0, :, :] = dim0Array
+    grid.T[1, :, :] = numpy.array(dim1Array, copy=False)[:, None]
+    return grid
+
+
+def triangleStripGridIndices(dim0, dim1):
+    """Generate indices to draw a grid of vertices as a triangle strip.
+
+    Vertices are expected to be stored as row-major (i.e., C contiguous).
+
+    :param int dim0: The number of rows of vertices.
+    :param int dim1: The number of columns of vertices.
+    :return: The vertex indices
+    :rtype: 1D numpy.ndarray of uint32
+    """
+    assert dim0 >= 2
+    assert dim1 >= 2
+
+    # Filling a row of squares +
+    # an index before and one after for degenerated triangles
+    indices = numpy.empty((dim0 - 1, 2 * (dim1 + 1)), dtype=numpy.uint32)
+
+    # Init indices with minimum indices for each row of squares
+    indices[:] = (dim1 * numpy.arange(dim0 - 1, dtype=numpy.uint32))[:, None]
+
+    # Update indices with offset per row of squares
+    offset = numpy.arange(dim1, dtype=numpy.uint32)
+    indices[:, 1:-1:2] += offset
+    offset += dim1
+    indices[:, 2::2] += offset
+    indices[:, -1] += offset[-1]
+
+    # Remove extra indices for degenerated triangles before returning
+    return indices.ravel()[1:-1]
+
+    # Alternative:
+    # indices = numpy.zeros(2 * dim1 * (dim0 - 1) + 2 * (dim0 - 2),
+    #                      dtype=numpy.uint32)
+    #
+    # offset = numpy.arange(dim1, dtype=numpy.uint32)
+    # for d0Index in range(dim0 - 1):
+    #    start = 2 * d0Index * (dim1 + 1)
+    #    end = start + 2 * dim1
+    #    if d0Index != 0:
+    #        indices[start - 2] = offset[-1]
+    #        indices[start - 1] = offset[0]
+    #    indices[start:end:2] = offset
+    #    offset += dim1
+    #    indices[start + 1:end:2] = offset
+    # return indices
+
+
+def linesGridIndices(dim0, dim1):
+    """Generate indices to draw a grid of vertices as lines.
+
+    Vertices are expected to be stored as row-major (i.e., C contiguous).
+
+    :param int dim0: The number of rows of vertices.
+    :param int dim1: The number of columns of vertices.
+    :return: The vertex indices.
+    :rtype: 1D numpy.ndarray of uint32
+    """
+    # Horizontal and vertical lines
+    nbsegmentalongdim1 = 2 * (dim1 - 1)
+    nbsegmentalongdim0 = 2 * (dim0 - 1)
+
+    indices = numpy.empty(nbsegmentalongdim1 * dim0 +
+                          nbsegmentalongdim0 * dim1,
+                          dtype=numpy.uint32)
+
+    # Line indices over dim0
+    onedim1line = (numpy.arange(nbsegmentalongdim1,
+                                dtype=numpy.uint32) + 1) // 2
+    indices[:dim0 * nbsegmentalongdim1] = \
+        (dim1 * numpy.arange(dim0, dtype=numpy.uint32)[:, None] +
+         onedim1line[None, :]).ravel()
+
+    # Line indices over dim1
+    onedim0line = (numpy.arange(nbsegmentalongdim0,
+                                dtype=numpy.uint32) + 1) // 2
+    indices[dim0 * nbsegmentalongdim1:] = \
+        (numpy.arange(dim1, dtype=numpy.uint32)[:, None] +
+         dim1 * onedim0line[None, :]).ravel()
+
+    return indices
+
+
+# intersection ################################################################
+
+def angleBetweenVectors(refVector, vectors, norm=None):
+    """Return the angle between 2 vectors.
+
+    :param refVector: Coordinates of the reference vector.
+    :type refVector: numpy.ndarray of shape: (NCoords,)
+    :param vectors: Coordinates of the vector(s) to get angle from reference.
+    :type vectors: numpy.ndarray of shape: (NCoords,) or (NbVector, NCoords)
+    :param norm: A direction vector giving an orientation to the angles
+                 or None.
+    :returns: The angles in radians in [0, pi] if norm is None
+              else in [0, 2pi].
+    :rtype: float or numpy.ndarray of shape (NbVectors,)
+    """
+    singlevector = len(vectors.shape) == 1
+    if singlevector:  # Make it a 2D array for the computation
+        vectors = vectors.reshape(1, -1)
+
+    assert len(refVector.shape) == 1
+    assert len(vectors.shape) == 2
+    assert len(refVector) == vectors.shape[1]
+
+    # Normalize vectors
+    refVector /= numpy.linalg.norm(refVector)
+    vectors = numpy.array([v / numpy.linalg.norm(v) for v in vectors])
+
+    dots = numpy.sum(refVector * vectors, axis=-1)
+    angles = numpy.arccos(numpy.clip(dots, -1., 1.))
+    if norm is not None:
+        signs = numpy.sum(norm * numpy.cross(refVector, vectors), axis=-1) < 0.
+        angles[signs] = numpy.pi * 2. - angles[signs]
+
+    return angles[0] if singlevector else angles
+
+
+def segmentPlaneIntersect(s0, s1, planeNorm, planePt):
+    """Compute the intersection of a segment with a plane.
+
+    :param s0: First end of the segment
+    :type s0: 1D numpy.ndarray-like of length 3
+    :param s1: Second end of the segment
+    :type s1: 1D numpy.ndarray-like of length 3
+    :param planeNorm: Normal vector of the plane.
+    :type planeNorm: numpy.ndarray of shape: (3,)
+    :param planePt: A point of the plane.
+    :type planePt: numpy.ndarray of shape: (3,)
+    :return: The intersection points. The number of points goes
+             from 0 (no intersection) to 2 (segment in the plane)
+    :rtype: list of numpy.ndarray
+    """
+    s0, s1 = numpy.asarray(s0), numpy.asarray(s1)
+
+    segdir = s1 - s0
+    dotnormseg = numpy.dot(planeNorm, segdir)
+    if dotnormseg == 0:
+        # line and plane are parallels
+        if numpy.dot(planeNorm, planePt - s0) == 0:  # segment is in plane
+            return [s0, s1]
+        else:  # No intersection
+            return []
+
+    alpha = - numpy.dot(planeNorm, s0 - planePt) / dotnormseg
+    if 0. <= alpha <= 1.:  # Intersection with segment
+        return [s0 + alpha * segdir]
+    else:  # intersection outside segment
+        return []
+
+
+def boxPlaneIntersect(boxVertices, boxLineIndices, planeNorm, planePt):
+    """Return intersection points between a box and a plane.
+
+    :param boxVertices: Position of the corners of the box.
+    :type boxVertices: numpy.ndarray with shape: (8, 3)
+    :param boxLineIndices: Indices of the box edges.
+    :type boxLineIndices: numpy.ndarray-like with shape: (12, 2)
+    :param planeNorm: Normal vector of the plane.
+    :type planeNorm: numpy.ndarray of shape: (3,)
+    :param planePt: A point of the plane.
+    :type planePt: numpy.ndarray of shape: (3,)
+    :return: The found intersection points
+    :rtype: numpy.ndarray with 2 dimensions
+    """
+    segments = numpy.take(boxVertices, boxLineIndices, axis=0)
+
+    points = set()  # Gather unique intersection points
+    for seg in segments:
+        for point in segmentPlaneIntersect(seg[0], seg[1], planeNorm, planePt):
+            points.add(tuple(point))
+    points = numpy.array(list(points))
+
+    if len(points) <= 2:
+        return numpy.array(())
+    elif len(points) == 3:
+        return points
+    else:  # len(points) > 3
+        # Order point to have a polyline lying on the unit cube's faces
+        vectors = points - numpy.mean(points, axis=0)
+        angles = angleBetweenVectors(vectors[0], vectors, planeNorm)
+        points = numpy.take(points, numpy.argsort(angles), axis=0)
+        return points
+
+
+# Plane #######################################################################
+
+class Plane(event.Notifier):
+    """Object handling a plane and notifying plane changes.
+
+    :param point: A point on the plane.
+    :type point: 3-tuple of float.
+    :param normal: Normal of the plane.
+    :type normal: 3-tuple of float.
+    """
+
+    def __init__(self, point=(0., 0., 0.), normal=(0., 0., 1.)):
+        super(Plane, self).__init__()
+
+        assert len(point) == 3
+        self._point = numpy.array(point, copy=True, dtype=numpy.float32)
+        assert len(normal) == 3
+        self._normal = numpy.array(normal, copy=True, dtype=numpy.float32)
+        self.notify()
+
+    def setPlane(self, point=None, normal=None):
+        """Set plane point and normal and notify.
+
+        :param point: A point on the plane.
+        :type point: 3-tuple of float or None.
+        :param normal: Normal of the plane.
+        :type normal: 3-tuple of float or None.
+        """
+        planechanged = False
+
+        if point is not None:
+            assert len(point) == 3
+            point = numpy.array(point, copy=True, dtype=numpy.float32)
+            if not numpy.all(numpy.equal(self._point, point)):
+                self._point = point
+                planechanged = True
+
+        if normal is not None:
+            assert len(normal) == 3
+            normal = numpy.array(normal, copy=True, dtype=numpy.float32)
+
+            norm = numpy.linalg.norm(normal)
+            if norm != 0.:
+                normal /= norm
+
+            if not numpy.all(numpy.equal(self._normal, normal)):
+                self._normal = normal
+                planechanged = True
+
+        if planechanged:
+            _logger.debug('Plane updated:\n\tpoint: %s\n\tnormal: %s',
+                          str(self._point), str(self._normal))
+            self.notify()
+
+    @property
+    def point(self):
+        """A point on the plane."""
+        return self._point.copy()
+
+    @point.setter
+    def point(self, point):
+        self.setPlane(point=point)
+
+    @property
+    def normal(self):
+        """The (normalized) normal of the plane."""
+        return self._normal.copy()
+
+    @normal.setter
+    def normal(self, normal):
+        self.setPlane(normal=normal)
+
+    @property
+    def parameters(self):
+        """Plane equation parameters: a*x + b*y + c*z + d = 0."""
+        return numpy.append(self._normal,
+                            - numpy.dot(self._point, self._normal))
+
+    @parameters.setter
+    def parameters(self, parameters):
+        assert len(parameters) == 4
+        parameters = numpy.array(parameters, dtype=numpy.float32)
+
+        # Normalize normal
+        norm = numpy.linalg.norm(parameters[:3])
+        if norm != 0:
+            parameters /= norm
+
+        normal = parameters[:3]
+        point = - parameters[3] * normal
+        self.setPlane(point, normal)
+
+    @property
+    def isPlane(self):
+        """True if a plane is defined (i.e., ||normal|| != 0)."""
+        return numpy.any(self.normal != 0.)
+
+    def move(self, step):
+        """Move the plane of step along the normal."""
+        self.point += step * self.normal
diff --git a/silx/gui/plot3d/scene/viewport.py b/silx/gui/plot3d/scene/viewport.py
new file mode 100644
index 0000000..83cda43
--- /dev/null
+++ b/silx/gui/plot3d/scene/viewport.py
@@ -0,0 +1,492 @@
+# coding: utf-8
+# /*##########################################################################
+#
+# Copyright (c) 2015-2017 European Synchrotron Radiation Facility
+#
+# Permission is hereby granted, free of charge, to any person obtaining a copy
+# of this software and associated documentation files (the "Software"), to deal
+# in the Software without restriction, including without limitation the rights
+# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+# copies of the Software, and to permit persons to whom the Software is
+# furnished to do so, subject to the following conditions:
+#
+# The above copyright notice and this permission notice shall be included in
+# all copies or substantial portions of the Software.
+#
+# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
+# THE SOFTWARE.
+#
+# ###########################################################################*/
+"""This module provides a class to control a viewport on the rendering window.
+
+The :class:`Viewport` describes a Viewport rendering a scene.
+The attribute :attr:`scene` is the root group of the scene tree.
+:class:`RenderContext` handles the current state during rendering.
+"""
+
+from __future__ import absolute_import, division, unicode_literals
+
+__authors__ = ["T. Vincent"]
+__license__ = "MIT"
+__date__ = "25/07/2016"
+
+
+import numpy
+
+from silx.gui.plot.Colors import rgba
+
+from ..._glutils import gl
+
+from . import camera
+from . import event
+from . import transform
+from .function import DirectionalLight, ClippingPlane
+
+
+class RenderContext(object):
+    """Handle a current rendering context.
+
+    An instance of this class is passed to rendering method through
+    the scene during render.
+
+    User should NEVER use an instance of this class beyond the method
+    it is passed to as an argument (i.e., do not keep a reference to it).
+
+    :param Viewport viewport: The viewport doing the rendering.
+    :param Context glContext: The operating system OpenGL context in use.
+    """
+
+    def __init__(self, viewport, glContext):
+        self._viewport = viewport
+        self._glContext = glContext
+        self._transformStack = [viewport.camera.extrinsic]
+        self._clipPlane = ClippingPlane(normal=(0., 0., 0.))
+
+    @property
+    def viewport(self):
+        """Viewport doing the current rendering"""
+        return self._viewport
+
+    @property
+    def glCtx(self):
+        """The OpenGL context in use"""
+        return self._glContext
+
+    @property
+    def objectToCamera(self):
+        """The current transform from object to camera coords.
+
+        Do not modify.
+        """
+        return self._transformStack[-1]
+
+    @property
+    def projection(self):
+        """Projection transform.
+
+        Do not modify.
+        """
+        return self.viewport.camera.intrinsic
+
+    @property
+    def objectToNDC(self):
+        """The transform from object to NDC (this includes projection).
+
+        Do not modify.
+        """
+        return transform.StaticTransformList(
+            (self.projection, self.objectToCamera))
+
+    def pushTransform(self, transform_, multiply=True):
+        """Push a :class:`Transform` on the transform stack.
+
+        :param Transform transform_: The transform to add to the stack.
+        :param bool multiply:
+            True (the default) to multiply with the top of the stack,
+            False to push the transform as is without multiplication.
+        """
+        if multiply:
+            assert len(self._transformStack) >= 1
+            transform_ = transform.StaticTransformList(
+                (self._transformStack[-1], transform_))
+
+        self._transformStack.append(transform_)
+
+    def popTransform(self):
+        """Pop the transform on top of the stack.
+
+        :return: The Transform that is popped from the stack.
+        """
+        assert len(self._transformStack) > 1
+        return self._transformStack.pop()
+
+    @property
+    def clipper(self):
+        """The current clipping plane
+        """
+        return self._clipPlane
+
+    def setClipPlane(self, point=(0., 0., 0.), normal=(0., 0., 0.)):
+        """Set the clipping plane to use
+
+        For now only handles a single clipping plane.
+
+        :param point: A point of the plane
+        :type point: 3-tuple of float
+        :param normal: Normal vector of the plane or (0, 0, 0) for no clipping
+        :type normal: 3-tuple of float
+        """
+        self._clipPlane = ClippingPlane(point, normal)
+
+
+class Viewport(event.Notifier):
+    """Rendering a single scene through a camera in part of a framebuffer.
+
+    :param int framebuffer: The framebuffer ID this viewport is rendering into
+    """
+
+    def __init__(self, framebuffer=0):
+        from . import Group  # Here to avoid cyclic import
+        super(Viewport, self).__init__()
+        self._dirty = True
+        self._origin = 0, 0
+        self._size = 1, 1
+        self._framebuffer = int(framebuffer)
+        self.scene = Group()  # The stuff to render, add overlaid scenes?
+        self.scene._setParent(self)
+        self.scene.addListener(self._changed)
+        self._background = 0., 0., 0., 1.
+        self._camera = camera.Camera(fovy=30., near=1., far=100.,
+                                     position=(0., 0., 12.))
+        self._camera.addListener(self._changed)
+        self._transforms = transform.TransformList([self._camera])
+
+        self._light = DirectionalLight(direction=(0., 0., -1.),
+                                       ambient=(0.3, 0.3, 0.3),
+                                       diffuse=(0.7, 0.7, 0.7))
+        self._light.addListener(self._changed)
+
+    @property
+    def transforms(self):
+        """Proxy of camera transforms.
+
+        Do not modify the list.
+        """
+        return self._transforms
+
+    def _changed(self, *args, **kwargs):
+        """Callback handling scene updates"""
+        self._dirty = True
+        self.notify()
+
+    @property
+    def dirty(self):
+        """True if scene is dirty and needs redisplay."""
+        return self._dirty
+
+    def resetDirty(self):
+        """Mark the scene as not being dirty.
+
+        To call after rendering.
+        """
+        self._dirty = False
+
+    @property
+    def background(self):
+        """Background color of the viewport (4-tuple of float in [0, 1]"""
+        return self._background
+
+    @background.setter
+    def background(self, color):
+        color = rgba(color)
+        if self._background != color:
+            self._background = color
+            self._changed()
+
+    @property
+    def camera(self):
+        """The camera used to render the scene."""
+        return self._camera
+
+    @property
+    def light(self):
+        """The light used to render the scene."""
+        return self._light
+
+    @property
+    def origin(self):
+        """Origin (ox, oy) of the viewport in pixels"""
+        return self._origin
+
+    @origin.setter
+    def origin(self, origin):
+        ox, oy = origin
+        origin = int(ox), int(oy)
+        if origin != self._origin:
+            self._origin = origin
+            self._changed()
+
+    @property
+    def size(self):
+        """Size (width, height) of the viewport in pixels"""
+        return self._size
+
+    @size.setter
+    def size(self, size):
+        w, h = size
+        size = int(w), int(h)
+        if size != self._size:
+            self._size = size
+
+            self.camera.intrinsic.size = size
+            self._changed()
+
+    @property
+    def shape(self):
+        """Shape (height, width) of the viewport in pixels.
+
+        This is a convenient wrapper to the inverse of size.
+        """
+        return self._size[1], self._size[0]
+
+    @shape.setter
+    def shape(self, shape):
+        self.size = shape[1], shape[0]
+
+    @property
+    def framebuffer(self):
+        """The framebuffer ID this viewport is rendering into (int)"""
+        return self._framebuffer
+
+    @framebuffer.setter
+    def framebuffer(self, framebuffer):
+        self._framebuffer = int(framebuffer)
+
+    def render(self, glContext):
+        """Perform the rendering of the viewport
+
+        :param Context glContext: The context used for rendering"""
+        # Get a chance to run deferred delete
+        glContext.cleanGLGarbage()
+
+        # OpenGL set-up: really need to be done once
+        ox, oy = self.origin
+        w, h = self.size
+        gl.glViewport(ox, oy, w, h)
+
+        gl.glEnable(gl.GL_SCISSOR_TEST)
+        gl.glScissor(ox, oy, w, h)
+
+        gl.glEnable(gl.GL_BLEND)
+        gl.glBlendFunc(gl.GL_SRC_ALPHA, gl.GL_ONE_MINUS_SRC_ALPHA)
+
+        gl.glEnable(gl.GL_DEPTH_TEST)
+        gl.glDepthFunc(gl.GL_LEQUAL)
+        gl.glDepthRange(0., 1.)
+
+        # gl.glEnable(gl.GL_POLYGON_OFFSET_FILL)
+        # gl.glPolygonOffset(1., 1.)
+
+        gl.glHint(gl.GL_LINE_SMOOTH_HINT, gl.GL_NICEST)
+        gl.glEnable(gl.GL_LINE_SMOOTH)
+
+        gl.glClearColor(*self.background)
+
+        # Prepare OpenGL
+        gl.glClear(gl.GL_COLOR_BUFFER_BIT |
+                   gl.GL_STENCIL_BUFFER_BIT |
+                   gl.GL_DEPTH_BUFFER_BIT)
+
+        ctx = RenderContext(self, glContext)
+        self.scene.render(ctx)
+        self.scene.postRender(ctx)
+
+    def adjustCameraDepthExtent(self):
+        """Update camera depth extent to fit the scene bounds.
+
+        Only near and far planes are updated.
+        The scene might still not be fully visible
+        (e.g., if spanning behind the viewpoint with perspective projection).
+        """
+        bounds = self.scene.bounds(transformed=True)
+        bounds = self.camera.extrinsic.transformBounds(bounds)
+
+        if isinstance(self.camera.intrinsic, transform.Perspective):
+            # This needs to be reworked
+            zbounds = - bounds[:, 2]
+            zextent = max(numpy.fabs(zbounds[0] - zbounds[1]), 0.0001)
+            near = max(zextent / 1000., 0.95 * zbounds[1])
+            far = max(near + 0.1, 1.05 * zbounds[0])
+
+            self.camera.intrinsic.setDepthExtent(near, far)
+        elif isinstance(self.camera.intrinsic, transform.Orthographic):
+            # Makes sure z bounds are included
+            border = max(abs(bounds[:, 2]))
+            self.camera.intrinsic.setDepthExtent(-border, border)
+        else:
+            raise RuntimeError('Unsupported camera', self.camera.intrinsic)
+
+    def resetCamera(self):
+        """Change camera to have the whole scene in the viewing frustum.
+
+        It updates the camera position and depth extent.
+        Camera sight direction and up are not affected.
+        """
+        self.camera.resetCamera(self.scene.bounds(transformed=True))
+
+    def orbitCamera(self, direction, angle=1.):
+        """Rotate the camera around center of the scene.
+
+        :param str direction: Direction of movement relative to image plane.
+                              In: 'up', 'down', 'left', 'right'.
+        :param float angle: he angle in degrees of the rotation.
+        """
+        bounds = self.scene.bounds(transformed=True)
+        center = 0.5 * (bounds[0] + bounds[1])
+        self.camera.orbit(direction, center, angle)
+
+    def moveCamera(self, direction, step=0.1):
+        """Move the camera relative to the image plane.
+
+        :param str direction: Direction relative to image plane.
+                              One of: 'up', 'down', 'left', 'right',
+                              'forward', 'backward'.
+        :param float step: The ratio of data to step for each pan.
+        """
+        bounds = self.scene.bounds(transformed=True)
+        bounds = self.camera.extrinsic.transformBounds(bounds)
+        center = 0.5 * (bounds[0] + bounds[1])
+        ndcCenter = self.camera.intrinsic.transformPoint(
+            center, perspectiveDivide=True)
+
+        step *= 2.  # NDC has size 2
+
+        if direction == 'up':
+            ndcCenter[1] -= step
+        elif direction == 'down':
+            ndcCenter[1] += step
+
+        elif direction == 'right':
+            ndcCenter[0] -= step
+        elif direction == 'left':
+            ndcCenter[0] += step
+
+        elif direction == 'forward':
+            ndcCenter[2] += step
+        elif direction == 'backward':
+            ndcCenter[2] -= step
+
+        else:
+            raise ValueError('Unsupported direction: %s' % direction)
+
+        newCenter = self.camera.intrinsic.transformPoint(
+            ndcCenter, direct=False, perspectiveDivide=True)
+
+        self.camera.move(direction, numpy.linalg.norm(newCenter - center))
+
+    def windowToNdc(self, winX, winY, checkInside=True):
+        """Convert position from window to normalized device coordinates.
+
+        If window coordinates are int, they are moved half a pixel
+        to be positioned at the center of pixel.
+
+        :param winX: X window coord, origin left.
+        :param winY: Y window coord, origin top.
+        :param bool checkInside: If True, returns None if position is
+                                 outside viewport.
+        :return: (x, y) Normalize device coordinates in [-1, 1] or None.
+                 Origin center, x to the right, y goes upward.
+        """
+        ox, oy = self._origin
+        width, height = self.size
+
+        # If int, move it to the center of pixel
+        if isinstance(winX, int):
+            winX += 0.5
+        if isinstance(winY, int):
+            winY += 0.5
+
+        x, y = winX - ox, winY - oy
+
+        if checkInside and (x < 0. or x > width or y < 0. or y > height):
+            return None  # Out of viewport
+
+        ndcx = 2. * x / float(width) - 1.
+        ndcy = 1. - 2. * y / float(height)
+        return ndcx, ndcy
+
+    def ndcToWindow(self, ndcX, ndcY, checkInside=True):
+        """Convert position from normalized device coordinates (NDC) to window.
+
+        :param float ndcX: X NDC coord.
+        :param float ndcY: Y NDC coord.
+        :param bool checkInside: If True, returns None if position is
+                                 outside viewport.
+        :return: (x, y) window coordinates or None.
+                 Origin top-left, x to the right, y goes downward.
+        """
+        if (checkInside and
+                (ndcX < -1. or ndcX > 1. or ndcY < -1. or ndcY > 1.)):
+            return None  # Outside viewport
+
+        ox, oy = self._origin
+        width, height = self.size
+
+        winx = ox + width * 0.5 * (ndcX + 1.)
+        winy = oy + height * 0.5 * (1. - ndcY)
+        return winx, winy
+
+    def _pickNdcZGL(self, x, y):
+        """Retrieve depth from depth buffer and return corresponding NDC Z.
+
+        :param int x: In pixels in window coordinates, origin left.
+        :param int y: In pixels in window coordinates, origin top.
+        :return: Normalize device Z coordinate of depth in [-1, 1]
+                 or None if outside viewport.
+        :rtype: float or None
+        """
+        ox, oy = self._origin
+        width, height = self.size
+
+        x = int(x)
+        y = height - int(y)  # Invert y coord
+
+        if x < ox or x > ox + width or y < oy or y > oy + height:
+            # Outside viewport
+            return None
+
+        # Get depth from depth buffer in [0., 1.]
+        # Bind used framebuffer to get depth
+        gl.glBindFramebuffer(gl.GL_FRAMEBUFFER, self.framebuffer)
+        depth = gl.glReadPixels(
+            x, y, 1, 1, gl.GL_DEPTH_COMPONENT, gl.GL_FLOAT)[0]
+        gl.glBindFramebuffer(gl.GL_FRAMEBUFFER, 0)
+        # This is not GL|ES friendly
+
+        # Z in NDC in [-1., 1.]
+        return float(depth) * 2. - 1.
+
+    def _getXZYGL(self, x, y):
+        ndc = self.windowToNdc(x, y)
+        if ndc is None:
+            return None  # Outside viewport
+        ndcz = self._pickNdcZGL(x, y)
+        ndcpos = numpy.array((ndc[0], ndc[1], ndcz, 1.), dtype=numpy.float32)
+
+        camerapos = self.camera.intrinsic.transformPoint(
+            ndcpos, direct=False, perspectiveDivide=True)
+
+        scenepos = self.camera.extrinsic.transformPoint(camerapos,
+                                                        direct=False)
+        return scenepos[:3]
+
+    def pick(self, x, y):
+        pass
+        # ndcX, ndcY = self.windowToNdc(x, y)
+        # ndcNearPt = ndcX, ndcY, -1.
+        # ndcFarPT = ndcX, ndcY, 1.
diff --git a/silx/gui/plot3d/scene/window.py b/silx/gui/plot3d/scene/window.py
new file mode 100644
index 0000000..ad7e6e5
--- /dev/null
+++ b/silx/gui/plot3d/scene/window.py
@@ -0,0 +1,420 @@
+# coding: utf-8
+# /*##########################################################################
+#
+# Copyright (c) 2015-2017 European Synchrotron Radiation Facility
+#
+# Permission is hereby granted, free of charge, to any person obtaining a copy
+# of this software and associated documentation files (the "Software"), to deal
+# in the Software without restriction, including without limitation the rights
+# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+# copies of the Software, and to permit persons to whom the Software is
+# furnished to do so, subject to the following conditions:
+#
+# The above copyright notice and this permission notice shall be included in
+# all copies or substantial portions of the Software.
+#
+# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
+# THE SOFTWARE.
+#
+# ###########################################################################*/
+"""This module provides a class for Viewports rendering on the screen.
+
+The :class:`Window` renders a list of Viewports in the current framebuffer.
+The rendering can be performed in an off-screen framebuffer that is only
+updated when the scene has changed and not each time Qt is requiring a repaint.
+
+The :class:`Context` and :class:`ContextGL2` represent the operating system
+OpenGL context and handle OpenGL resources.
+"""
+
+from __future__ import absolute_import, division, unicode_literals
+
+__authors__ = ["T. Vincent"]
+__license__ = "MIT"
+__date__ = "10/01/2017"
+
+
+import weakref
+import numpy
+
+from ..._glutils import gl
+from ... import _glutils
+
+from . import event
+
+
+class Context(object):
+    """Correspond to an operating system OpenGL context.
+
+    User should NEVER use an instance of this class beyond the method
+    it is passed to as an argument (i.e., do not keep a reference to it).
+
+    :param glContextHandle: System specific OpenGL context handle.
+    """
+
+    def __init__(self, glContextHandle):
+        self._context = glContextHandle
+        self._isCurrent = False
+        self._devicePixelRatio = 1.0
+
+    @property
+    def isCurrent(self):
+        """Whether this OpenGL context is the current one or not."""
+        return self._isCurrent
+
+    def setCurrent(self, isCurrent=True):
+        """Set the state of the OpenGL context to reflect OpenGL state.
+
+        This should not be called from the scene graph, only in the
+        wrapper that handle the OpenGL context to reflect its state.
+
+        :param bool isCurrent: The state of the system OpenGL context.
+        """
+        self._isCurrent = bool(isCurrent)
+
+    @property
+    def devicePixelRatio(self):
+        """Ratio between device and device independent pixels (float)
+
+        This is useful for font rendering.
+        """
+        return self._devicePixelRatio
+
+    @devicePixelRatio.setter
+    def devicePixelRatio(self, ratio):
+        assert ratio > 0
+        self._devicePixelRatio = float(ratio)
+
+    def __enter__(self):
+        self.setCurrent(True)
+        return self
+
+    def __exit__(self, exc_type, exc_val, exc_tb):
+        self.setCurrent(False)
+
+    @property
+    def glContext(self):
+        """The handle to the OpenGL context provided by the system."""
+        return self._context
+
+    def cleanGLGarbage(self):
+        """This is releasing OpenGL resource that are no longer used."""
+        pass
+
+
+class ContextGL2(Context):
+    """Handle a system GL2 context.
+
+    User should NEVER use an instance of this class beyond the method
+    it is passed to as an argument (i.e., do not keep a reference to it).
+
+    :param glContextHandle: System specific OpenGL context handle.
+    """
+    def __init__(self, glContextHandle):
+        super(ContextGL2, self).__init__(glContextHandle)
+
+        self._programs = {}  # GL programs already compiled
+        self._vbos = {}  # GL Vbos already set
+        self._vboGarbage = []  # Vbos waiting to be discarded
+
+    # programs
+
+    def prog(self, vertexShaderSrc, fragmentShaderSrc):
+        """Cache program within context.
+
+        WARNING: No clean-up.
+        """
+        assert self.isCurrent
+        key = vertexShaderSrc, fragmentShaderSrc
+        prog = self._programs.get(key, None)
+        if prog is None:
+            prog = _glutils.Program(vertexShaderSrc, fragmentShaderSrc)
+            self._programs[key] = prog
+        return prog
+
+    # VBOs
+
+    def makeVbo(self, data=None, sizeInBytes=None,
+                usage=None, target=None):
+        """Create a VBO in this context with the data.
+
+        Current limitations:
+
+        - One array per VBO
+        - Do not support sharing VertexBuffer across VboAttrib
+
+        Automatically discards the VBO when the returned
+        :class:`VertexBuffer` istance is deleted.
+
+        :param numpy.ndarray data: 2D array of data to store in VBO or None.
+        :param int sizeInBytes: Size of the VBO or None.
+                                It should be <= data.nbytes if both are given.
+        :param usage: OpenGL usage define in VertexBuffer._USAGES.
+        :param target: OpenGL target in VertexBuffer._TARGETS.
+        :return: The VertexBuffer created in this context.
+        """
+        assert self.isCurrent
+        vbo = _glutils.VertexBuffer(data, sizeInBytes, usage, target)
+        vboref = weakref.ref(vbo, self._deadVbo)
+        # weakref is hashable as far as target is
+        self._vbos[vboref] = vbo.name
+        return vbo
+
+    def makeVboAttrib(self, data, usage=None, target=None):
+        """Create a VBO from data and returns the associated VBOAttrib.
+
+        Automatically discards the VBO when the returned
+        :class:`VBOAttrib` istance is deleted.
+
+        :param numpy.ndarray data: 2D array of data to store in VBO or None.
+        :param usage: OpenGL usage define in VertexBuffer._USAGES.
+        :param target: OpenGL target in VertexBuffer._TARGETS.
+        :returns: A VBOAttrib instance created in this context.
+        """
+        assert self.isCurrent
+        vbo = self.makeVbo(data, usage=usage, target=target)
+
+        assert len(data.shape) <= 2
+        dimension = 1 if len(data.shape) == 1 else data.shape[1]
+
+        return _glutils.VertexBufferAttrib(
+            vbo,
+            type_=_glutils.numpyToGLType(data.dtype),
+            size=data.shape[0],
+            dimension=dimension,
+            offset=0,
+            stride=0)
+
+    def _deadVbo(self, vboRef):
+        """Callback handling dead VBOAttribs."""
+        vboid = self._vbos.pop(vboRef)
+        if self.isCurrent:
+            # Direct delete if context is active
+            gl.glDeleteBuffers(vboid)
+        else:
+            # Deferred VBO delete if context is not active
+            self._vboGarbage.append(vboid)
+
+    def cleanGLGarbage(self):
+        """Delete OpenGL resources that are pending for destruction.
+
+        This requires the associated OpenGL context to be active.
+        This is meant to be called before rendering.
+        """
+        assert self.isCurrent
+        if self._vboGarbage:
+            vboids = self._vboGarbage
+            gl.glDeleteBuffers(vboids)
+            self._vboGarbage = []
+
+
+class Window(event.Notifier):
+    """OpenGL Framebuffer where to render viewports
+
+    :param str mode: Rendering mode to use:
+
+        - 'direct' to render everything for each render call
+        - 'framebuffer' to cache viewport rendering in a texture and
+          update the texture only when needed.
+    """
+
+    _position = numpy.array(((-1., -1., 0., 0.),
+                             (1., -1., 1., 0.),
+                             (-1., 1., 0., 1.),
+                             (1., 1., 1., 1.)),
+                            dtype=numpy.float32)
+
+    _shaders = ("""
+        attribute vec4 position;
+        varying vec2 textureCoord;
+
+        void main(void) {
+            gl_Position = vec4(position.x, position.y, 0., 1.);
+            textureCoord = position.zw;
+        }
+        """,
+                """
+        uniform sampler2D texture;
+        varying vec2 textureCoord;
+
+        void main(void) {
+            gl_FragColor = texture2D(texture, textureCoord);
+        }
+        """)
+
+    def __init__(self, mode='framebuffer'):
+        super(Window, self).__init__()
+        self._dirty = True
+        self._size = 0, 0
+        self._contexts = {}  # To map system GL context id to Context objects
+        self._viewports = event.NotifierList()
+        self._viewports.addListener(self._updated)
+        self._framebufferid = 0
+        self._framebuffers = {}  # Cache of framebuffers
+
+        assert mode in ('direct', 'framebuffer')
+        self._isframebuffer = mode == 'framebuffer'
+
+    @property
+    def dirty(self):
+        """True if this object or any attached viewports is dirty."""
+        for viewport in self._viewports:
+            if viewport.dirty:
+                return True
+        return self._dirty
+
+    @property
+    def size(self):
+        """Size (width, height) of the window in pixels"""
+        return self._size
+
+    @size.setter
+    def size(self, size):
+        w, h = size
+        size = int(w), int(h)
+        if size != self._size:
+            self._size = size
+            self._dirty = True
+            self.notify()
+
+    @property
+    def shape(self):
+        """Shape (height, width) of the window in pixels.
+
+        This is a convenient wrapper to the reverse of size.
+        """
+        return self._size[1], self._size[0]
+
+    @shape.setter
+    def shape(self, shape):
+        self.size = shape[1], shape[0]
+
+    @property
+    def viewports(self):
+        """List of viewports to render in the corresponding framebuffer"""
+        return self._viewports
+
+    @viewports.setter
+    def viewports(self, iterable):
+        self._viewports.removeListener(self._updated)
+        self._viewports = event.NotifierList(iterable)
+        self._viewports.addListener(self._updated)
+        self._dirty = True
+
+    def _updated(self, source, *args, **kwargs):
+        if source is not self:
+            self._dirty = True
+            self.notify(*args, **kwargs)
+
+    framebufferid = property(lambda self: self._framebufferid,
+                             doc="Framebuffer ID used to perform rendering")
+
+    def grab(self, glcontext):
+        """Returns the raster of the scene as an RGB numpy array
+
+        :returns: OpenGL scene RGB bitmap
+        :rtype: numpy.ndarray of uint8 of dimension (height, width, 3)
+        """
+        height, width = self.shape
+        image = numpy.empty((height, width, 3), dtype=numpy.uint8)
+
+        gl.glBindFramebuffer(gl.GL_FRAMEBUFFER, self.framebufferid)
+        gl.glPixelStorei(gl.GL_PACK_ALIGNMENT, 1)
+        gl.glReadPixels(
+            0, 0, width, height, gl.GL_RGB, gl.GL_UNSIGNED_BYTE, image)
+        gl.glBindFramebuffer(gl.GL_FRAMEBUFFER, 0)
+
+        # glReadPixels gives bottom to top,
+        # while images are stored as top to bottom
+        image = numpy.flipud(image)
+
+        return numpy.array(image, copy=False, order='C')
+
+    def render(self, glcontext, devicePixelRatio):
+        """Perform the rendering of attached viewports
+
+        :param glcontext: System identifier of the OpenGL context
+        :param float devicePixelRatio:
+            Ratio between device and device-independent pixels
+        """
+        if glcontext not in self._contexts:
+            self._contexts[glcontext] = ContextGL2(glcontext)  # New context
+
+        with self._contexts[glcontext] as context:
+            context.devicePixelRatio = devicePixelRatio
+            if self._isframebuffer:
+                self._renderWithOffscreenFramebuffer(context)
+            else:
+                self._renderDirect(context)
+
+        self._dirty = False
+
+    def _renderDirect(self, context):
+        """Perform the direct rendering of attached viewports
+
+        :param Context context: Object wrapping OpenGL context
+        """
+        for viewport in self._viewports:
+            viewport.framebuffer = self.framebufferid
+            viewport.render(context)
+            viewport.resetDirty()
+
+    def _renderWithOffscreenFramebuffer(self, context):
+        """Renders viewports in a texture and render this texture on screen.
+
+        The texture is updated only if viewport or size has changed.
+
+        :param ContextGL2 context: Object wrappign OpenGL context
+        """
+        if self.dirty or context not in self._framebuffers:
+            # Need to redraw framebuffer content
+
+            if (context not in self._framebuffers or
+                    self._framebuffers[context].shape != self.shape):
+                # Need to rebuild framebuffer
+
+                if context in self._framebuffers:
+                    self._framebuffers[context].discard()
+
+                fbo = _glutils.FramebufferTexture(gl.GL_RGBA,
+                                                  shape=self.shape,
+                                                  minFilter=gl.GL_NEAREST,
+                                                  magFilter=gl.GL_NEAREST,
+                                                  wrap=gl.GL_CLAMP_TO_EDGE)
+                self._framebuffers[context] = fbo
+                self._framebufferid = fbo.name
+
+            # Render in framebuffer
+            with self._framebuffers[context]:
+                self._renderDirect(context)
+
+        # Render framebuffer texture to screen
+        fbo = self._framebuffers[context]
+        height, width = fbo.shape
+
+        program = context.prog(*self._shaders)
+        program.use()
+
+        gl.glViewport(0, 0, width, height)
+        gl.glDisable(gl.GL_BLEND)
+        gl.glDisable(gl.GL_DEPTH_TEST)
+        gl.glDisable(gl.GL_SCISSOR_TEST)
+        # gl.glScissor(0, 0, width, height)
+        gl.glClearColor(0., 0., 0., 0.)
+        gl.glClear(gl.GL_COLOR_BUFFER_BIT)
+        gl.glUniform1i(program.uniforms['texture'], fbo.texture.texUnit)
+        gl.glEnableVertexAttribArray(program.attributes['position'])
+        gl.glVertexAttribPointer(program.attributes['position'],
+                                 4,
+                                 gl.GL_FLOAT,
+                                 gl.GL_FALSE,
+                                 0,
+                                 self._position)
+        fbo.texture.bind()
+        gl.glDrawArrays(gl.GL_TRIANGLE_STRIP, 0, len(self._position))
+        gl.glBindTexture(gl.GL_TEXTURE_2D, 0)