{"worksheets": [{"cells": [{"cell_type": "markdown", "metadata": {}, "source": ["Geodesic Farthest Point Sampling", "================================", "", "*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 explores the use geodesic computation to perform image", "sampling."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["from __future__ import division", "import nt_toolbox as nt", "from nt_solutions import fastmarching_5_sampling_2d as solutions", "%matplotlib inline", "%load_ext autoreload", "%autoreload 2"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Voronoi Segmentation", "--------------------", "A geodesic Voronoi segementation is obtained by computing a geodesic distance from", "multiple starting points.", "", "", "Compute an image with bumps."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["n = 256; % size of the image", "sigma = n/ 8; % width of the bumps", "[Y, X] = meshgrid(1: n, 1: n)", "x = n/ 4; y = n/ 4", "M = exp(-((X-x).^2 + (Y-y).^2)/ (2*sigma^2))", "x = 3*n/ 4; y = 3*n/ 4", "M = M + exp(-((X-x).^2 + (Y-y).^2)/ (2*sigma^2))"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Compute a metric by rescaling |M|."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["W = rescale(M, 1e-2, 1)"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Create random starting points."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["m = 20", "pstart = floor(rand(2, m)*(n-1)) + 1"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Perform the propagation using the Fast Marching."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["[D, Z, Q] = perform_fast_marching(1./ W, pstart)"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Display the distance map."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["subplot(1, 2, 1)", "hold on", "imageplot(perform_hist_eq(D, 'linear'));  title('Geodesic distance')", "plot(pstart(2, : ), pstart(1, : ), 'r.')", "subplot(1, 2, 2)", "hold on", "imageplot(Q); title('Voronoi')", "plot(pstart(2, : ), pstart(1, : ), 'r.')", "colormap jet(256)"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Geodesic Delaunay Triangulation", "-------------------------------", "A geodesic Delaunay triangulation is obtained by linking starting points", "whose Voronoi cells touch. This is the dual of the original Voronoi segmentation."]}, {"cell_type": "markdown", "metadata": {}, "source": ["__Exercise 1__", "", "Using |Q|, compute the faces |faces| of the Delaunay triangulation.", "To that end,", "extract each quad of values |Q(i,j),Q(i+1,j),Q(i+1,j+1),Q(i,j+1)|,", "and add a new face when three of these four values are different (this", "corresponds to a Voronoi point).", "Display the obtained triangulation.", "ompute quad of values"]}, {"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": ["Display the obtained triangulation."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["subplot(1, 2, 1)", "hold on", "imageplot(Q, 'Voronoi'); axis ij", "plot(pstart(2, : ), pstart(1, : ), 'r.')", "subplot(1, 2, 2)", "plot_triangulation(pstart, faces, M)", "colormap jet(256)"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Compare this Geodesic Delaunay triangulation with the Euclidean", "triangulation."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["faces_euc = compute_delaunay(pstart)"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Display."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["subplot(1, 2, 1)", "plot_triangulation(pstart, faces, M); title('Geodesic')", "subplot(1, 2, 2)", "plot_triangulation(pstart, faces_euc, M); title('Euclidean')", "colormap jet(256)"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Farthest Point Sampling", "-----------------------", "To sample point uniformly according to the geodesic distance, one can use", "an iterative farthest point sampling scheme.", "", "", "Construct a metric that is large in area where we want to sample more", "points. The front should move slowly in high density sampling region."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["W = rescale(M, 3*1e-1, 1)"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Choose the initial point."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["vertex = [1; 1]"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Compute the geodesic distance."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["[D, Z, Q] = perform_fast_marching(1./ W, vertex)"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Choose the second point as the farthest point."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["[tmp, i] = max(D(: ))", "[x, y] = ind2sub([n n], i)", "vertex(: , end + 1) = [x; y]"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Display distance and points."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["subplot(1, 2, 1)", "hold on", "imageplot(W, 'Metric'); axis ij", "plot(vertex(2, 1), vertex(1, 1), 'r.')", "plot(vertex(2, 2), vertex(1, 2), 'b.')", "subplot(1, 2, 2)", "hold on", "imageplot(perform_hist_eq(D, 'linear'), 'Distance'); axis ij", "plot(vertex(2, 1), vertex(1, 1), 'r.')", "plot(vertex(2, 2), vertex(1, 2), 'b.')", "colormap jet(256)"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Update the value of the distance map with a partial propagation from the last added point."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["options.constraint_map = D", "[D1, Z, Q] = perform_fast_marching(1./ W, vertex(: , end), options)"]}, {"cell_type": "markdown", "metadata": {}, "source": ["display old/new"]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["imageplot(convert_distance_color(D, M, 1), 'Update distance', 1, 3, 1)", "imageplot(convert_distance_color(D1, M, 1), 'Update distance', 1, 3, 2)", "imageplot(convert_distance_color(min(D, D1), M, 1), 'New distance', 1, 3, 3)"]}, {"cell_type": "markdown", "metadata": {}, "source": ["update"]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["D = min(D, D1)"]}, {"cell_type": "markdown", "metadata": {}, "source": ["__Exercise 2__", "", "Iterate the sampling process to add more and more points."]}, {"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": ["__Exercise 3__", "", "Display the geodesic Delaunay triangulation corresponding to the", "sampling"]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["solutions.exo3()"]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["## Insert your code here."]}, {"cell_type": "markdown", "metadata": {}, "source": ["Lloyd Relaxation", "----------------", "The farthest point sampling strategy is greedy and does not move the", "positions of the points once they are seeded.", "", "", "To enhance the sampling, it is possible to relocate iteratively the", "points at the center of the Voronoi cells. This corresponds to the Lloyd", "algorithm, first developped for vector quantization.", "", "", "First with a constant metric."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["n = 512", "W = ones(n)"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Seed random points."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["p = 40", "vertex = floor(rand(2, p)*n-1) + 1"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Compute Voronoi partition."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["[D, Z, Q] = perform_fast_marching(1./ W, vertex)"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Display Vornoi cells."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["clf; hold on", "imageplot(Q')", "h = plot(vertex(1, : ), vertex(2, : ), 'k.')", "set(h, 'MarkerSize', 15)", "colormap(jet(256))"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Re-center each point at the barycenter of its cell."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["for i in 1: p:", "    [x, y] = ind2sub(size(W), find(Q = =i))", "    vertex(: , i) = [mean(x); mean(y)]"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Display updated partitions."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["[D, Z, Q] = perform_fast_marching(1./ W, vertex)", "clf; hold on", "imageplot(Q')", "h = plot(vertex(1, : ), vertex(2, : ), 'k.')", "set(h, 'MarkerSize', 15)", "colormap(jet(256))"]}, {"cell_type": "markdown", "metadata": {}, "source": ["__Exercise 4__", "", "Perform the Lloyd iterative algorithm.", "etric", "eed random points."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["solutions.exo4()"]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["## Insert your code here."]}, {"cell_type": "markdown", "metadata": {}, "source": ["Now we define a non-constent metric."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["n = 256", "x = linspace(-1, 1, n)", "[Y, X] = meshgrid(x, x)", "sigma = .3", "W = exp(-(X.^2 + Y.^2)/ (2*sigma^2))", "W = rescale(W, .02, 1)"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Display."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["imageplot(W)"]}, {"cell_type": "markdown", "metadata": {}, "source": ["__Exercise 5__", "", "Perform the Lloyd iterative algorithm.", "eed random points."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["solutions.exo5()"]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["## Insert your code here."]}]}], "nbformat": 3, "metadata": {"name": ""}, "nbformat_minor": 0}