{"worksheets": [{"cells": [{"cell_type": "markdown", "metadata": {}, "source": ["Geodesic Medial Axsis", "=====================", "", "*Important:* Please read the [installation page](http://gpeyre.github.io/numerical-tours/installation_python/) 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}$"]}, {"cell_type": "markdown", "metadata": {}, "source": ["This tour studies the computation of the medial axis using the Fast", "Marching."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["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"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Voronoi Diagram", "---------------", "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."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["n = 200", "W = load_image('mountain', n)", "W = rescale(W, .25, 1)"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Select seed points."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["pstart = [[20; 20] [120; 100] [180; 30] [60; 160]]", "nbound = size(pstart, 2)"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Display the map and the points."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["ms = 20", "clf; hold on", "imageplot(W)", "h = plot(pstart(2, : ), pstart(1, : ), '.r'); set(h, 'MarkerSize', ms)"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Compute the geodesic distant to the whole set of points."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["[D, S, Q] = perform_fast_marching(W, pstart)"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Display the geodesic distance."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["clf; hold on", "imageplot(convert_distance_color(D, W))", "h = plot(pstart(2, : ), pstart(1, : ), '.r'); set(h, 'MarkerSize', ms)"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Display the Voronoi Segmentation."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["clf; hold on", "imageplot(Q)", "h = plot(pstart(2, : ), pstart(1, : ), '.r'); set(h, 'MarkerSize', ms)", "colormap jet(256)"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Medial Axis from the Voronoi Map", "--------------------------------", "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."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["G = grad(Q)"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Take it modulo |nbound|."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["G(G <-nbound/ 2) = G(G <-nbound/ 2) + nbound", "G(G >nbound/ 2) = G(G >nbound/ 2) - nbound"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Compute the norm of the gadient."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["G = sqrt(sum(G.^2, 3))"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Compute the medial axis by thresholding the gradient magnitude."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["B = 1 - (G >.1)"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Display."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["clf; hold on", "imageplot(B)", "h = plot(pstart(2, : ), pstart(1, : ), '.r'); set(h, 'MarkerSize', ms)"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Skeleton of a Shape ", "--------------------", "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."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["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"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Compute its boundary, that is going to be the set of starting points."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["pstart = compute_shape_boundary(M)", "nbound = size(pstart, 2)"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Display the metric."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["lw = 2", "clf; hold on", "imageplot(-M)", "h = plot(pstart(2, : ), pstart(1, : ), 'r'); set(h, 'LineWidth', lw); axis ij"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Parameters for the Fast Marching: constant speed |W|, but retricted using |L| to the", "inside of the shape."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["W = ones(n)", "L = zeros(n)-Inf; L(M = =1) = + Inf"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Compute the fast marching, from the boundary points."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["options.constraint_map = L", "[D, S, Q] = perform_fast_marching(W, pstart, options)", "D(M = =0) = Inf"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Display the distance function to the boundary."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["hold on", "display_shape_function(D)", "h = plot(pstart(2, : ), pstart(1, : ), 'r'); set(h, 'LineWidth', lw); axis ij"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Display the index of the closest boundary point."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["hold on", "display_shape_function(Q)", "h = plot(pstart(2, : ), pstart(1, : ), 'r'); set(h, 'LineWidth', lw); axis ij"]}, {"cell_type": "markdown", "metadata": {}, "source": ["__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."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["solutions.exo1()"]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["## Insert your code here."]}, {"cell_type": "markdown", "metadata": {}, "source": ["__Exercise 2__", "", "Display the Skeleton obtained for different threshold values."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["solutions.exo2()"]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["## Insert your code here."]}, {"cell_type": "markdown", "metadata": {}, "source": ["Skeletton", "---------"]}]}], "nbformat": 3, "metadata": {"name": ""}, "nbformat_minor": 0}