4

I am trying to use skimage.measure.marching_cubes_lewiner to resolve some isosurface f(x,y,z)=0. In my case f is strongly nonlinear, and is best mapped when coordinates are given with logarithmic spacing. Because the marching cubes wants a regular grid, to build the voxels, I am working on a meshgrid of coordinates X,Y,Z which correspond to the log10 of my original coordinates, so that my isosurface is equivalently given by f(10**X,10**Y,10**Z)=0. Everything would be fine, if it were not for the fact that, say I am working with X,Y,Z in [-1.5,2]^3 (equivalent to x,y,z in [0.03,100.]^3), the vertex coordinates of the solution given by skimage.measure.marching_cubes_lewiner, are not in this cube.

Following the answer to another related question on SO, I thought it could be due to the fact that, probably the algorithm works thinking of a unitary volume, so that I need to set the right spacing input argument in my call of skimage.measure.marching_cubes_lewiner. In this fashion, say I am mapping my function f on a grid of N points per coordinate, so that I am increasing exponents by numpy.diff([-1.5,2])/N per coordinate, I accordingly call:

import numpy as np
from skimage import measure as msr

def f(x,y,z):
    val = ... # some lengthy code to define my implicit function
    return val 

# Define ranges of my coordinates
xRange = [0.03,100.]
yRange = [0.03,100.]
zRange = [0.03,100.]

XRange = np.log10(xRange)
YRange = np.log10(yRange)
ZRange = np.log10(zRange)   

# Create regular grid
N = 50  # number of points per coordinate
X,Y,Z = np.mesh[XRange[0]:XRange[1]:N*1j,
                YRange[0]:YRange[1]:N*1j,
                ZRange[0]:ZRange[1]:N*1j]
F = f(10**X,10**Y,10**Z)
sol,_,_,_ = skimage.measure.marching_cubes_lewiner(F,0.0,spacing(np.diff(XRange)/N,np.diff(YRange)/N,np.diff(ZRange)/N))

yet, unexpectedly, the coordinates of the solution points generally seem in [0,Vx]*[0,Vy]*[0,Vz] with Vx>XRange[-1], Vy>YRange[-1] and Vz>ZRange[-1]. I have no clue of why this happens and how I could properly rescale the coordinates of my isosurface solution, to the real units of my problem.

maurizio
  • 745
  • 1
  • 7
  • 25
  • Anyone can help on this? Is it perhaps that I simply need to shift the vertex coordinates by `XRange[0]`, `YRange[0]` and `ZRange[0]`. Doing this trick results seems consistent, but I cannot cross-check unless I outsource a different implementation of MCs in python (other than the `skimage.measure` one). – maurizio Mar 06 '19 at 07:55
  • In case someone gets stucked on this same issue, it seems that my above comment is the answer to my question. But before I post it separately, unless someone confirms it meanwhile, I want first to compare the solutions for my isosurface solution `skimage.measure` 's MCs with those independently estimated by global solvers. – maurizio Mar 06 '19 at 14:12
  • Maybe you could open an issue in the scikit-image GitHub repo, linking this question. https://github.com/scikit-image/scikit-image/issues Maybe it is just a documentation issue or maybe it is an implementation issue. In any case, I think you should let them know :) – Eskapp Mar 08 '19 at 19:24

0 Answers0