Concave hull with missing edges

Question

I have a binary mask from which I want to extract a contour. All "outside" edges of the binary mask are "true edges" with very high probability. Keeping these edges fixed, my goal is now to interpolate the "missing" edges (example image and desired results below). I have tried using Delaunay triangulation, without much success (see code below). However, I'm not even sure that's the best approach, as I will lose some of the detail from "true edges" in the process.

Is Delaunay triangulation appropriate here? If so, what's wrong with the code below? Is there a better algorithm to solve this kind of problem? How could this be implemented in Python?

Current code (does not work)

import numpy as np
from scipy.spatial import Delaunay
from descartes import PolygonPatch
from shapely.geometry import MultiLineString
from shapely.ops import cascaded_union, polygonize
from skimage.draw import polygon
from skimage import segmentation, morphology

def triangulate(mask, alpha=1000):

  contours = measure.find_contours(mask, 0.8)

  points = []

  for n, contour in enumerate(contours):

    for m in xrange(0, len(contour[:, 0])):
      y = contour[:, 0][m]
      x = contour[:, 1][m]
      points.append([y, x])

  points = np.asarray(points)

  tri = Delaunay(points)
  edges = set()
  edge_points = []

  def add_edge(i, j):
    if (i, j) in edges or (j, i) in edges: return
    edges.add( (i, j) )
    edge_points.append(points[ [i, j] ])

  for ia, ib, ic in tri.vertices:
    pa = points[ia]
    pb = points[ib]
    pc = points[ic]

    # Lengths of sides of triangle
    a = math.sqrt((pa[0]-pb[0])**2 + (pa[1]-pb[1])**2)
    b = math.sqrt((pb[0]-pc[0])**2 + (pb[1]-pc[1])**2)
    c = math.sqrt((pc[0]-pa[0])**2 + (pc[1]-pa[1])**2)

    # Semiperimeter of triangle
    s = (a + b + c)/2.0

    # Area of triangle by Heron's formula
    try:
      area = math.sqrt(s*(s-a)*(s-b)*(s-c))
      circum_r = a*b*c/(4.0*area)

      # Here's the radius filter.
      #if circum_r < 1.0/alpha:
      if circum_r < alpha:
          add_edge(ia, ib)
          add_edge(ib, ic)
          add_edge(ic, ia)
    except:
      print('Triangulation error')

  m = MultiLineString(edge_points)
  triangles = list(polygonize(m))

  poly = PolygonPatch(cascaded_union(triangles), alpha=0.5)

  vertices = poly.get_path().vertices

  rr, cc = polygon(vertices[:,0], vertices[:,1])

  img = np.zeros(im.shape)
  img[rr, cc] = 1

  return img

Answer 1

As nikie pointed out, the morphological snakes will solve your problem if performance is not a requirement.

When you increase the number of iterations, the snake tends to cover up perfectly the shape, but that's not what you want, so you need to figure out a combination of parameters that work.

I played around a little bit and found a good approximation:

For this shape I used the following code:

def test_concave():
    # Load the image.
    imgcolor = imread("testimages/concave_hull.jpg")/255.0
    img = rgb2gray(imgcolor)

    # g(I)
    gI = morphsnakes.gborders(img, alpha=1000, sigma=7)

    # Morphological GAC. Initialization of the level-set.
    mgac = morphsnakes.MorphGAC(gI, smoothing=3, threshold=0.55, balloon=-1)
    mgac.levelset = circle_levelset(img.shape, (50, 199), 200)

    # Visual evolution.
    ppl.figure()
    morphsnakes.evolve_visual(mgac, num_iters=62, background=imgcolor)


if __name__ == '__main__':
    print """"""
    test_concave()
    ppl.show()