Filling a cavity with a level set methodΒΆ
This example shows how to fill a concave cavity using a level set method.

Edge of segmented region shown as a red outline on top of the original image. The animation shows evolution of the edge as a function of iteration number.
def levelset_fill_cavity():
"""
Demonstrates how to fill a non-convex cavity using a level-set method.
"""
# Create an image that contains a rectangle with a cavity in its edge.
# For this example we make a 2D image for easier visualization,
# but everything should work the same for a 3D volume image, too.
c = 30 # Radius of the image
geom = pi.newimage(ImageDataType.UINT8, 2 * c, 2 * c)
pi.box(geom, [c - 15, c - 15, 0], 30, 128) # The box
pi.sphere(geom, [c, c, 0], 8, 0) # Part of cavity
pi.sphere(geom, [c, c - 6, 0], 6, 0) # Part of cavity
pi.sphere(geom, [c, c - 12, 0], 6, 0) # Part of cavity
# Noise to make the image more realistic
pi.noise(geom, 0, 20)
# Save the geometry
pi.writetif(geom, output_file('levelset_geom'))
# Initialize an image that holds the level set function.
# After iteration below, the segmentation is the part of phi
# that has positive value.
phi = pi.newlike(geom, ImageDataType.FLOAT32)
# Update phi iteratively
for n in np.arange(0, 200):
print(f"Iteration {n}" )
# Construct force field that acts on the surface defined by phi
# The force will be the sum of three terms.
F = pi.newlike(phi)
# First term: the force is positive inside the object and negative everywhere else.
# This makes the surface take the shape of the object.
pi.copy(geom, F)
pi.divide(F, 128)
pi.subtract(F, 0.5)
# Second term: Penalize high curvature by making curvy
# regions less curved.
kappa = pi.newlike(phi)
pi.meancurvature(phi, kappa, 0.5)
# Multiply the curvature values to scale them correctly.
pi.multiply(kappa, -5)
# Remove negative curvature values. They correspond to
# convex shapes, and we want to zero those so that
# they don't have any effect on the surface.
pi.max(kappa, 0)
# Add kappa to the force term
pi.add(F, kappa)
# Third term: Normal force
# This term makes the surface move towards its normal.
# This term is not required in this example.
#L = pi.newlike(phi)
#pi.gradientmagnitude(phi, L, 0.75, 0)
#pi.multiply(L, 0.01)
#pi.add(F, L)
# Smooth the total force a little bit
# This is not strictly by the book, but smoothing seems to
# make phi converge faster to a smoother result.
tmp = pi.newlike(F)
pi.gaussfilter(F, tmp, 0.5, BoundaryCondition.NEAREST)
pi.set(F, tmp)
# Multiply by time step
dt = 1
pi.multiply(F, dt)
# Add force*dt to the surface
pi.add(phi, F)
# Convert phi to segmentation and save it for visualization purposes
vis = pi.newlike(phi)
pi.copy(phi, vis)
pi.threshold(vis, 0)
pi.convert(vis, ImageDataType.UINT8)
pi.writetif(vis, output_file(f"./levelset_result/{n}"))