Important: Please read the installation page for details about how to install the toolboxes. $\newcommand{\dotp}[2]{\langle #1, #2 \rangle}$ $\newcommand{\enscond}[2]{\lbrace #1, #2 \rbrace}$ $\newcommand{\pd}[2]{ \frac{ \partial #1}{\partial #2} }$ $\newcommand{\umin}[1]{\underset{#1}{\min}\;}$ $\newcommand{\umax}[1]{\underset{#1}{\max}\;}$ $\newcommand{\umin}[1]{\underset{#1}{\min}\;}$ $\newcommand{\uargmin}[1]{\underset{#1}{argmin}\;}$ $\newcommand{\norm}[1]{\|#1\|}$ $\newcommand{\abs}[1]{\left|#1\right|}$ $\newcommand{\choice}[1]{ \left\{ \begin{array}{l} #1 \end{array} \right. }$ $\newcommand{\pa}[1]{\left(#1\right)}$ $\newcommand{\diag}[1]{{diag}\left( #1 \right)}$ $\newcommand{\qandq}{\quad\text{and}\quad}$ $\newcommand{\qwhereq}{\quad\text{where}\quad}$ $\newcommand{\qifq}{ \quad \text{if} \quad }$ $\newcommand{\qarrq}{ \quad \Longrightarrow \quad }$ $\newcommand{\ZZ}{\mathbb{Z}}$ $\newcommand{\CC}{\mathbb{C}}$ $\newcommand{\RR}{\mathbb{R}}$ $\newcommand{\EE}{\mathbb{E}}$ $\newcommand{\Zz}{\mathcal{Z}}$ $\newcommand{\Ww}{\mathcal{W}}$ $\newcommand{\Vv}{\mathcal{V}}$ $\newcommand{\Nn}{\mathcal{N}}$ $\newcommand{\NN}{\mathcal{N}}$ $\newcommand{\Hh}{\mathcal{H}}$ $\newcommand{\Bb}{\mathcal{B}}$ $\newcommand{\Ee}{\mathcal{E}}$ $\newcommand{\Cc}{\mathcal{C}}$ $\newcommand{\Gg}{\mathcal{G}}$ $\newcommand{\Ss}{\mathcal{S}}$ $\newcommand{\Pp}{\mathcal{P}}$ $\newcommand{\Ff}{\mathcal{F}}$ $\newcommand{\Xx}{\mathcal{X}}$ $\newcommand{\Mm}{\mathcal{M}}$ $\newcommand{\Ii}{\mathcal{I}}$ $\newcommand{\Dd}{\mathcal{D}}$ $\newcommand{\Ll}{\mathcal{L}}$ $\newcommand{\Tt}{\mathcal{T}}$ $\newcommand{\si}{\sigma}$ $\newcommand{\al}{\alpha}$ $\newcommand{\la}{\lambda}$ $\newcommand{\ga}{\gamma}$ $\newcommand{\Ga}{\Gamma}$ $\newcommand{\La}{\Lambda}$ $\newcommand{\si}{\sigma}$ $\newcommand{\Si}{\Sigma}$ $\newcommand{\be}{\beta}$ $\newcommand{\de}{\delta}$ $\newcommand{\De}{\Delta}$ $\newcommand{\phi}{\varphi}$ $\newcommand{\th}{\theta}$ $\newcommand{\om}{\omega}$ $\newcommand{\Om}{\Omega}$
This tour studies the computation of the medial axis using the Fast Marching.
from __future__ import division
import nt_toolbox as nt
from nt_solutions import shapes_6_medialaxis as solutions
%matplotlib inline
%load_ext autoreload
%autoreload 2
The Voronoi diagram is the segmentation of the image given by the region of influence of the set of starting points.
Load a distance map.
n = 200
W = load_image('mountain', n)
W = rescale(W, .25, 1)
Select seed points.
pstart = [[20; 20] [120; 100] [180; 30] [60; 160]]
nbound = size(pstart, 2)
Display the map and the points.
ms = 20
clf; hold on
imageplot(W)
h = plot(pstart(2, : ), pstart(1, : ), '.r'); set(h, 'MarkerSize', ms)
Compute the geodesic distant to the whole set of points.
[D, S, Q] = perform_fast_marching(W, pstart)
Display the geodesic distance.
clf; hold on
imageplot(convert_distance_color(D, W))
h = plot(pstart(2, : ), pstart(1, : ), '.r'); set(h, 'MarkerSize', ms)
Display the Voronoi Segmentation.
clf; hold on
imageplot(Q)
h = plot(pstart(2, : ), pstart(1, : ), '.r'); set(h, 'MarkerSize', ms)
colormap jet(256)
The medial axis is difficult to extract from the singularity of the distance map. It is much more robust to extract it from the discontinuities in the Voronoi index map |Q|.
Compute the derivative, the gradient.
G = grad(Q)
Take it modulo |nbound|.
G(G <-nbound/ 2) = G(G <-nbound/ 2) + nbound
G(G >nbound/ 2) = G(G >nbound/ 2) - nbound
Compute the norm of the gadient.
G = sqrt(sum(G.^2, 3))
Compute the medial axis by thresholding the gradient magnitude.
B = 1 - (G >.1)
Display.
clf; hold on
imageplot(B)
h = plot(pstart(2, : ), pstart(1, : ), '.r'); set(h, 'MarkerSize', ms)
The sekeleton, also called Medial Axis, is the set of points where the geodesic distance is singular.
A binary shape is represented as a binary image.
n = 200
name = 'chicken'
M = load_image(name, n)
M = perform_blurring(M, 5)
M = double(rescale(M) >.5)
if M(1) = =1
M = 1-M
Compute its boundary, that is going to be the set of starting points.
pstart = compute_shape_boundary(M)
nbound = size(pstart, 2)
Display the metric.
lw = 2
clf; hold on
imageplot(-M)
h = plot(pstart(2, : ), pstart(1, : ), 'r'); set(h, 'LineWidth', lw); axis ij
Parameters for the Fast Marching: constant speed |W|, but retricted using |L| to the inside of the shape.
W = ones(n)
L = zeros(n)-Inf; L(M = =1) = + Inf
Compute the fast marching, from the boundary points.
options.constraint_map = L
[D, S, Q] = perform_fast_marching(W, pstart, options)
D(M = =0) = Inf
Display the distance function to the boundary.
hold on
display_shape_function(D)
h = plot(pstart(2, : ), pstart(1, : ), 'r'); set(h, 'LineWidth', lw); axis ij
Display the index of the closest boundary point.
hold on
display_shape_function(Q)
h = plot(pstart(2, : ), pstart(1, : ), 'r'); set(h, 'LineWidth', lw); axis ij
Exercise 1
Compute the norm of the gradient |G| modulo |nbound|. Be careful to remove the boundary of the shape from this indicator. Display the thresholded gradient map. radient ompute the norm of the gadient. emove the boundary to the skeletton.
solutions.exo1()
## Insert your code here.
Exercise 2
Display the Skeleton obtained for different threshold values.
solutions.exo2()
## Insert your code here.