In [1]:

    

import numpy as np
import skimage.measure
# plot cube with based on edge=scale

scale = 500

shape = (np.asarray([3, 3, 3]) * scale).astype(np.int)
data = np.zeros(shape)
data[(1*scale):(2*scale), (1*scale):(2*scale), (1*scale):(2*scale)] = 2
# data[10:20, 10:20, 10:20] = 2

out = skimage.measure.marching_cubes(
    data, level=1)
vertices, faces = out
surface_area_numeric = skimage.measure.mesh_surface_area(vertices,faces)
surface_area_analytic = 6. * scale**2

print("S numeric : ", surface_area_numeric)
print("S analytic: ", surface_area_analytic)
print("Error: ", (surface_area_analytic-surface_area_numeric)/surface_area_analytic)