0

Goal

I want to plot a large number of cubes (arranged in a 3D grid) with different colors and opacities.

Current State and question

I have come up with a solution using vispy, but the performance is very poor - drawing takes very long and the window is very unresponsive. Also, there seem to be some glitches in the visualization, but I could live with those.

Is there a more efficient/elegant way to implement that? I am open to using other packages (I have tried open3d but found it difficult to specify colors and opacities - the documentation is not very verbose). However, I need to use python.

What I did so far

The first problem I had to solve with vispy was that I was unable to create cubes at custom positions. I therefore wrote a subclass that can do that:

import vispy.visuals
from vispy.geometry import create_box

class PositionedCubeVisual(vispy.visuals.BoxVisual):
    
    def __init__(self, size=1, position=(0, 0, 0), width_segments=1,
                 height_segments=1, depth_segments=1, planes=None,
                 vertex_colors=None, face_colors=None,
                 color=(0.5, 0.5, 1, 1), edge_color=None, **kwargs):
        
        vertices, filled_indices, outline_indices = create_box(
            size, size, size, width_segments, height_segments,
            depth_segments, planes)
        
        for column, pos in zip(vertices['position'].T, position):
            column += pos
        
        self._mesh = vispy.visuals.MeshVisual(vertices['position'], filled_indices,
                                vertex_colors, face_colors, color)
        if edge_color:
            self._border = vispy.visuals.MeshVisual(vertices['position'], outline_indices,
                                      color=edge_color, mode='lines')
        else:
            self._border = vispy.visuals.MeshVisual()

        vispy.visuals.CompoundVisual.__init__(self, [self._mesh, self._border], **kwargs)
        self.mesh.set_gl_state(polygon_offset_fill=True,
                               polygon_offset=(1, 1), depth_test=True)
        
PositionedCube = vispy.scene.visuals.create_visual_node(PositionedCubeVisual)

I then plot the cubes as follows:

import numpy as np
import vispy.scene

def plot_grid_cubes(x, y, z, c=None, size=1, alpha=0.1, edge_color="black", 
                    cmap="viridis", bgcolor="#FFFFFF"):
    
    canvas = vispy.scene.SceneCanvas(keys='interactive', show=True)
    
    view = canvas.central_widget.add_view()
    view.bgcolor = bgcolor
    view.camera = 'turntable'
    
    c = get_color_array(c, alpha, min(len(x), len(y), len(z)), cmap)
    
    for xx, yy, zz, cc in zip(x, y, z, c):
        
        cube = PositionedCube(size, (xx, yy, zz), color=cc, edge_color=edge_color, parent=view.scene)
    
    canvas.app.run()


def get_color_array(c, alpha, size, cmap):
    if c is not None:
        cmap = cm.get_cmap(cmap)
        if hasattr(c, "__iter__"):
            c = np.array(c, copy=True, dtype=float)
            c -= c.min() 
            c *= 255/c.max()
        return cmap(c.astype(int), alpha)
    else:
        color = np.ones((size, 4))
        color[:, 3] = alpha
        return color

This can then be applied as follows:

plot_grid_cubes([0, 1], [0, 1], [0, 1], c=[0.3, 0.5], alpha=[0.3, 0.8])

The example above works great, but it becomes poor if I plot thousands of cubes.

Samufi
  • 2,465
  • 3
  • 19
  • 43

2 Answers2

1

Regarding performance on vispy, you may want to read this:

Each Visual object in VisPy is an OpenGL Program consisting of at least a vertex shader and a fragment shader (see Modern OpenGL). In general, except for some very specific cases, OpenGL Programs can only be executed one at a time by a single OpenGL context. This means that in your VisPy visualization each Visual object you tell VisPy to draw will extend how long each update (draw) takes. When frames per second (FPS) or responsiveness are a concern, this means each Visual you add reduces the performance of your visualization.

While VisPy is constantly striving to improve performance, there are things that you can do in the mean time (depending on your particular case). The most important change that you can make is to lower the number of Visual objects you have. For a lot of Visuals it is possible to combine them into one by putting a little extra work into the data you provide them. For example, instead of creating 10 separate LineVisuals, create 1 LineVisual that draws 10 lines. While this is a simple example, the same concept applies to other types of Visuals in VisPy and more complex use cases. As a last resort for the most complex cases, a custom Visual (custom shader code) may be necessary. Before writing a custom Visual, check with VisPy maintainers by asking a question on gitter or creating a question as a GitHub issue.

Now for the BoxVisual this is a little difficult because as far as I can tell this "convenience" class doesn't allow you to make many boxes with a single BoxVisual instance. Since you are already comfortable making a Visual subclass I would recommend making the MeshVisuals yourself and providing the vertices for each box as one single position array.

As for not being able to specify position, this won't apply to your custom Visual class that will use the all-in-one array since you'll be providing each position at the beginning, but I thought I should still describe it. It is unfortunate that the BoxVisual is trying to be so convenient that it isn't helpful in this case since other Visuals allow you to pass your vertex positions on creation. In other cases or when you only want to make small modifications, typically what is done in VisPy is to use one or more "transforms" added to the Visual to shift (transform) the positions passed to the Visual. For example:

from vispy.visuals.transforms import STTransform

cube = ... create a cube visual ...
cube.transform = STTransform(scale=(1.0, 1.0, 1.0), translate=(0.0, 0.0, 0.0))

where you change the scale and translate values as needed to effect (X, Y, Z) coordinate values. After this, if you modify the cube.transform.translate = (new_x, new_y, new_z) property (or .scale or use another transform class) directly this has the benefit of only modifying that property on the GPU and not needing to recompute and resend the vertex positions (better performance).

djhoese
  • 3,567
  • 1
  • 27
  • 45
0

Inspired by @djhoese's answer, I have adjusted my initial solution attempt to build a single mesh only. The idea is to add the meshes of the individual cubes to one large mesh, which is then shown. The performance of the interactive visualization is as good as may be expected (some glitches remain); more efficient mesh creation algorithms are definitely possible.

First define a visual for multiple cubes:

from vispy.geometry import create_box
from vispy.scene import visuals
import vispy.scene
import vispy.visuals
from itertools import repeat

class MultiCubeVisual(vispy.visuals.BoxVisual):
    
    def __init__(self, size=1, positions=[(0, 0, 0)], width_segments=1,
                 height_segments=1, depth_segments=1, planes=None,
                 vertex_colors=None, face_colors=None,
                 color=(0.5, 0.5, 1, 0.5), edge_color=None, **kwargs):
        
        vertices = []
        filled_indices = []
        outline_indices = []
        all_vertex_colors = [] 
        all_face_colors = []
        
        if vertex_colors is None:
            vertex_colors = repeat(vertex_colors)
            all_vertex_colors = None
        
        if face_colors is None:
            face_colors = repeat(face_colors)
            all_face_colors = None
        
        maxIndex = 0
        for position, vertex_color, face_color in zip(positions, vertex_colors, face_colors):
            box_vertices, box_filled_indices, box_outline_indices = create_box(
                size, size, size, width_segments, height_segments,
                depth_segments, planes)
            
            for column, pos in zip(box_vertices['position'].T, position):
                column += pos

            box_filled_indices += maxIndex
            box_outline_indices += maxIndex
            vertexNumber = len(box_vertices)
            maxIndex += vertexNumber
            vertices.append(box_vertices)
            filled_indices.append(box_filled_indices)
            outline_indices.append(box_outline_indices)
            
            if vertex_color is not None:
                all_vertex_colors.extend([vertex_color] * vertexNumber)
            if face_color is not None:
                all_face_colors.extend([face_color] * len(box_filled_indices))
        
        vertices = np.concatenate(vertices)
        filled_indices = np.vstack(filled_indices)
        outline_indices = np.vstack(outline_indices)
        
        self._mesh = vispy.visuals.MeshVisual(vertices['position'], filled_indices,
                                all_vertex_colors, all_face_colors, color)
        if edge_color:
            self._border = vispy.visuals.MeshVisual(vertices['position'], outline_indices,
                                      color=edge_color, mode='lines')
        else:
            self._border = vispy.visuals.MeshVisual()

        vispy.visuals.CompoundVisual.__init__(self, [self._mesh, self._border], **kwargs)
        self.mesh.set_gl_state(polygon_offset_fill=True,
                               polygon_offset=(1, 1), depth_test=True)


MultiCube = vispy.scene.visuals.create_visual_node(MultiCubeVisual)

Now plotting works like this:

def get_color_array(c, alpha, size, cmap="viridis"):
    if c is not None:
        cmap = cm.get_cmap(cmap, 256)
        if hasattr(c, "__iter__"):
            c = np.array(c, copy=True, dtype=float)
            c -= c.min() 
            c *= 255/c.max()
        return cmap(c.astype(int), alpha)
    else:
        color = np.ones((size, 4))
        color[:, 3] = alpha
        return color
    

def plot_grid_cubes(x, y, z, c=None, size=1, alpha=0.1, edge_color="black", 
                    cmap="viridis", bgcolor="#FFFFFF", figure=None, show=True):
    
    if figure is None:
        figure = VispyFigure(bgcolor)
        
    mergedData = np.vstack((x, y, z)).T
    
    c = get_color_array(c, alpha, mergedData.shape[1], cmap)
    cubes = MultiCube(size, mergedData, edge_color=edge_color, face_colors=c, 
                      parent=figure.view.scene)

This can be applied as suggested in the question:

plot_grid_cubes([0, 1], [0, 1], [0, 1], c=[0.3, 0.5], alpha=[0.3, 0.8])
Samufi
  • 2,465
  • 3
  • 19
  • 43