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 explores the use of Fast Marching methods in 2D.
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
name = 'vessels'
options.nbdims = 3
M = read_bin(name, options)
M = rescale(M)
Size of the image (here it is a cube).
n = size(M, 1)
Such a volumetric dataset is more difficult to visualize than a standard 2D image. You can render slices along each X/Y/Z direction.
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)
We can display some horizontal slices.
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)
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.
h = vol3d('cdata', M, 'texture', '2D')
view(3); axis off
set up a colormap
colormap bone(256)
set up an alpha map
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
refresh the rendering
vol3d(h)
Exercise 1
Try with other alphamapping and colormapping et up a colormap et up an alpha map efresh the rendering
solutions.exo1()
## Insert your code here.
We can display an isosurface of the dataset (here we sub-sample to speed up the computation).
sel = 1: 2: n
isosurface(M(sel, sel, sel), .5)
axis('off')
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.
delta = 5
start_point = [107; 15; delta]
Compute the (inverse of the) potential that is small close to the value of |M| at the selected point
W = abs(M - M(start_point(1), start_point(2), start_point(3)))
Rescale above some threshold to avoid too small potentials.
W = rescale(W, 1e-2, 1)
Perform the front propagation.
options.nb_iter_max = Inf
[D, S] = perform_fast_marching(1./ W, start_point, options)
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).
imageplot(D(: , : , delta), '', 1, 2, 1)
imageplot(D(: , : , n-delta), '', 1, 2, 2)
colormap(jet(256))
Exercise 2
select the point (x,y) of minimum value in the slice |D(:,:,n-delta)|. hint: use functions 'min' and 'ind2sub'
solutions.exo2()
## Insert your code here.
Extract the geodesic by discrete descent.
options.method = 'discrete'
minpath = compute_geodesic(D, end_point, options)
Draw the path.
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)
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|.
solutions.exo3()
## Insert your code here.