{"worksheets": [{"cells": [{"cell_type": "markdown", "metadata": {}, "source": ["Fast Marching in 3D", "===================", "", "*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 of Fast Marching methods in 2D."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["from __future__ import division", "import nt_toolbox as nt", "from nt_solutions import fastmarching_2_3d as solutions", "%matplotlib inline", "%load_ext autoreload", "%autoreload 2"]}, {"cell_type": "markdown", "metadata": {}, "source": ["3D Volumetric Datasets", "----------------------", "A volumetric data is simply a 3D array.", "", "", "We load a volumetric data."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["name = 'vessels'", "options.nbdims = 3", "M = read_bin(name, options)", "M = rescale(M)"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Size of the image (here it is a cube)."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["n = size(M, 1)"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Such a volumetric dataset is more difficult to visualize than a standard 2D image.", "You can render slices along each X/Y/Z direction."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["imageplot(M(: , : , 50), 'X/ Y slice', 1, 3, 1)", "imageplot(squeeze(M(: , 50, : )), 'X/ Z slice', 1, 3, 2)", "imageplot(squeeze(M(50, : , : )), 'Y/ Z slice', 1, 3, 3)"]}, {"cell_type": "markdown", "metadata": {}, "source": ["We can display some horizontal slices."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["slices = round(linspace(10, n-10, 4))", "", "for i in 1: length(slices):", "    s = slices(i)", "    imageplot(M(: , : , s), strcat(['Z = ' num2str(s)]), 2, 2, i)"]}, {"cell_type": "markdown", "metadata": {}, "source": ["You can also perform a volumetric rendering.", "In order to do so, you need to set up a correct alpha mapping to make transparent some parts of the volume.", "Here, each time the options.center value is increased."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["h = vol3d('cdata', M, 'texture', '2D')", "view(3); axis off"]}, {"cell_type": "markdown", "metadata": {}, "source": ["set up a colormap"]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["colormap bone(256)"]}, {"cell_type": "markdown", "metadata": {}, "source": ["set up an alpha map"]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["options.sigma = .08; % control the width of the non-transparent region", "options.center = .4; % here a value in [0, 1]", "a = compute_alpha_map('gaussian', options); % you can plot(a) to see the alphamap"]}, {"cell_type": "markdown", "metadata": {}, "source": ["refresh the rendering"]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["vol3d(h)"]}, {"cell_type": "markdown", "metadata": {}, "source": ["__Exercise 1__", "", "Try with other alphamapping and colormapping", "et up a colormap", "et up an alpha map", "efresh the rendering"]}, {"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": ["We can display an isosurface of the dataset (here we sub-sample to speed", "up the computation)."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["sel = 1: 2: n", "", "isosurface(M(sel, sel, sel), .5)", "axis('off')"]}, {"cell_type": "markdown", "metadata": {}, "source": ["3D Shortest Paths", "-----------------", "The definition of shortest path extend to any dimension, for instance 3D.", "", "", "Geodesic distances can be computed on a 3D volume using the", "Fast Marching. The important point here is to define the correct", "potential field W that should be large in the region where you want the front to move fast.", "It means that geodesic will follow these regions.", "", "", "", "Select the starting points."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["delta = 5", "start_point = [107; 15; delta]"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Compute the (inverse of the) potential that is small  close", "to the value of |M| at the selected point"]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["W = abs(M - M(start_point(1), start_point(2), start_point(3)))"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Rescale above some threshold to avoid too small potentials."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["W = rescale(W, 1e-2, 1)"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Perform the front propagation."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["options.nb_iter_max = Inf", "[D, S] = perform_fast_marching(1./ W, start_point, options)"]}, {"cell_type": "markdown", "metadata": {}, "source": ["In order to extract a geodesic, we need to select an ending point", "and perform a descent of the distance function D from this starting point.", "The selection is done by choosing a point of low distance value in the slice D(:,:,n-delta)."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["imageplot(D(: , : , delta), '', 1, 2, 1)", "imageplot(D(: , : , n-delta), '', 1, 2, 2)", "colormap(jet(256))"]}, {"cell_type": "markdown", "metadata": {}, "source": ["__Exercise 2__", "", "select the point (x,y) of minimum value in the slice |D(:,:,n-delta)|.", "hint: use functions 'min' and 'ind2sub'"]}, {"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": ["Extract the geodesic by discrete descent."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["options.method = 'discrete'", "minpath = compute_geodesic(D, end_point, options)"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Draw the path."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["Dend = D(end_point(1), end_point(2), end_point(3))", "D1 = double(D <= Dend)", "", "plot_fast_marching_3d(M, D1, minpath, start_point, end_point)"]}, {"cell_type": "markdown", "metadata": {}, "source": ["__Exercise 3__", "", "Select other starting points. In order to do so, ask the user to", "click on a starting point in a given horizontal slice |W(:,:,delta)|.", "You can do this by using |ginput(1)|", "on the plane |Z=delta|."]}, {"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."]}]}], "nbformat": 3, "metadata": {"name": ""}, "nbformat_minor": 0}