CODEWITHSUNDEEP
pythonimage-processinginterpolationtriangulationdelaunay
user2398029
user2398029

Reputation: 6937

Concave hull with missing edges

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?

enter image description here

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

Upvotes: 4

Views: 2900

Answers (2)

mescarra
mescarra

Reputation: 725

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:

Convex Hull

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()

Upvotes: 6

Martin Beckett
Martin Beckett

Reputation: 96109

Convex hull to get the outer boundary then convexity defect to see how far away the nearest real edge is from the hull.

If the distance is small then use that contour as a real outer otherwise if the defect is large, then it is a real gap and use the convex boundary between the adjacent points

Upvotes: 2

Related Questions