Isosurface Extraction with Dual Contouring

#1

Hi together :slight_smile: ,

for my current project i try to solve the following problem:

My starting point is an arbitrary point cloud consisting of 3D-Coordinates that is to be voxelized to form a regular scalar field.

–> From this regular scalar field the extraction of an isosurface should be performed within Blender.

The hurdle to overcome is that i am not going to use the Marching Cubes algorithm, which is perfectly working as
it is implemented in the skimage.measures module of python and can be achieved by running the following code.
(The package skimage has to be copied to the following path ‘Blender Foundation\Blender\2.77\python\lib\site-packages’).
Maybe also the package scipy.

import numpy as np
import bpy
from skimage import measure
from skimage.draw import ellipsoid

def createMeshFromData(name, origin, verts, faces):
    # Create mesh and object
    me = bpy.data.meshes.new(name+'Mesh')
    ob = bpy.data.objects.new(name, me)
    ob.location = origin
    ob.show_name = True
 
    # Link object to scene and make active
    scn = bpy.context.scene
    scn.objects.link(ob)
    scn.objects.active = ob
    ob.select = True
 
    # Create mesh from given verts, faces.
    me.from_pydata(verts, [], faces)
    # Update mesh with new data
    me.update()    
    return ob

#http://scikit-image.org/docs/dev/auto_examples/plot_marching_cubes.html
# Generate a level set about zero of two identical ellipsoids in 3D
ellip_base = ellipsoid(6, 10, 16, levelset=True)
ellip_double = np.concatenate((ellip_base[:-1, ...],
                               ellip_base[2:, ...]), axis=0)
# Use marching cubes to obtain the surface mesh of these ellipsoids
isovalue = 0.5
# MC
verts, faces = measure.marching_cubes(ellip_double, isovalue)
# faces has to be of type int not int32numpy
ob = createMeshFromData('marching_cubes',(0,-30,0),verts,(faces.astype(int)).tolist())

# workaround to REMESH the MC-Mesh
ob2 = createMeshFromData('dual_contouring',(0,0,0),verts,(faces.astype(int)).tolist())

bpy.ops.object.modifier_add(type='REMESH')
bpy.context.object.modifiers["Remesh"].octree_depth = 5
bpy.context.object.modifiers["Remesh"].mode = 'SMOOTH'

My goal is now to substitute the Marching Cubes Algorithm by the Dual Contouring Algorithm by Ju, 2002.
source: ‘http://www.frankpetterson.com/publications/dualcontour/dualcontour.pdf’
and extract an isosurface based on a certain isovalue from the voxelgrid.

The dual contouring source code is already available as c++ code written by the author of the paper here (https://github.com/aewallin/dualcontouring).

An extended version of the algorithm (Polymender - ‘http://www.cs.wustl.edu/~taoju/code/polymender.htm’) ) is also already implemented in the source code in form of the remesh modifier by Nicholas Bishop here (C:\blender-2.77a\intern\dualcon\intern) / https://wiki.blender.org/index.php/Dev:Ref/Release_Notes/2.62/Remesh_Modifier.

The problem is that as the modifiers name connotes it needs a polygonal mesh as input and remeshes it using the Dual Contouring algorithm.
It works basically in this three steps:

1. Scan-conversion: Embed the input model (figure 2 (a)) in a uniformly spaced grid, and mark edges on the grid that intersect the polygons as intersection edges. For efficiency, cells containing intersection edges are stored in an octree (figure 2 (b)).
2. Sign generation: At the grid points, generate signs that are consistent with the intersection edges, so that each cell edge intersecting the model should exhibit a sign change (figure 2 (c)).
3. Surface reconstruction: Reconstruct a closed surface on the signed grid by contouring. Dual contouring can be used to reproduce sharp features when Hermite data is stored on the intersection edges (figure 2 (d)).

The readme file from the Polymender software says the pure dual contouring c++ code files can take two kinds of input:

.dcf (Dual Contouring Format) and
.sog (Signed Octree with Geometry).

Both formats store an octree grid with
inside/outside signs. DCF contains intersection points and normals on
grid edges, whereas SOG contains a single point location within each
non-empty grid cell. Both formats can be produced from a polygonal
model, via scan-conversion, using the Polymender software (Ju).

Is there any possibility to use at least parts of the already implemented source code in order to only extract an isosurface with dual contouring from a point cloud?

That means to jump in at for example step 3, the surface reconstruction, by for example providing the intersection points or the signed grid as .dcf or .sog file.

My initial thoughts were to calculate the interpolated intersections based on the isovalue from the voxelgrid and write them as Ju did in the octree grid as .dcf or .sog.
After that maybe only the DC-code can be used to contour the grid (https://github.com/aewallin/dualcontouring).

Another general question regarding the way to achieve a project like this one is:
What do you think is the most sophisticated way to access already existing c++ code within Blender? May i try to reuse/rewrite/abbreviate some of the
already existing source code? or access the c++ library via ctypes in a python addon? I would like to hear your suggestions.

I hope anyone of you maybe has already done something with dual contouring in Blender too and is able to give me some advide.
Thank you in advance!

Regards, Lukas