📄 content.html.svn-base
字号:
<!DOCTYPE html PUBLIC "-//W3C//DTD XHTML 1.0 Strict//EN"><html xmlns:mwsh="http://www.mathworks.com/namespace/mcode/v1/syntaxhighlight.dtd"> <head> <meta http-equiv="Content-Type" content="text/html; charset=utf-8"> <!--This HTML is auto-generated from an M-file.To make changes, update the M-file and republish this document. --> <title>Toolbox Fast Marching - A toolbox for Fast Marching and level sets computations</title> <meta name="generator" content="MATLAB 7.4"> <meta name="date" content="2008-08-16"> <meta name="m-file" content="content"> <LINK REL="stylesheet" HREF="style.css" TYPE="text/css"> </head> <body> <div class="content"> <h1>Toolbox Fast Marching - A toolbox for Fast Marching and level sets computations</h1> <introduction> <p>Copyright (c) 2008 Gabriel Peyre</p> </introduction> <h2>Contents</h2> <div> <ul> <li><a href="#2">2D Fast Marching Computations</a></li> <li><a href="#7">3D volumetric dataset</a></li> <li><a href="#10">3D Fast Marching</a></li> <li><a href="#13">Fast Marching on 3D meshes</a></li> <li><a href="#17">Fast Marching inside a 2D shape</a></li> <li><a href="#23">Farthest point sampling</a></li> <li><a href="#24">Anisotropic Fast Marching</a></li> </ul> </div> <p>First add some directories to the path</p><pre class="codeinput">path(path, <span class="string">'toolbox/'</span>);path(path, <span class="string">'data/'</span>);</pre><h2>2D Fast Marching Computations<a name="2"></a></h2> <p>The toolbox allows to compute a distance map from a set of starting point using an arbitrary isotropic metric. This distance map is then used to compute geodesic paths from any point that joins the closest starting point, according to the geodesic distance. </p> <p>We first load a 2D map, and ask to the user for the starting points.</p><pre class="codeinput">n = 400;[M,W] = load_potential_map(<span class="string">'road2'</span>, n);start_point = [16;219];end_point = [394;192];<span class="comment">% You can use instead the function</span><span class="comment">% [start_point,end_point] = pick_start_end_point(M);</span></pre><p>We then perform the front propagation</p><pre class="codeinput">clear <span class="string">options</span>;options.nb_iter_max = Inf;options.end_points = end_point; <span class="comment">% stop propagation when end point reached</span>[D,S] = perform_fast_marching(W, start_point, options);<span class="comment">% nicde color for display</span>A = convert_distance_color(D);imageplot({W A}, {<span class="string">'Potential map'</span> <span class="string">'Distance to starting point'</span>});colormap <span class="string">gray(256)</span>;</pre><img vspace="5" hspace="5" src="content_01.png"> <p>At last we can extract the geodesic curve, and display the result.</p><pre class="codeinput">gpath = compute_geodesic(D,end_point);plot_fast_marching_2d(W,S,gpath,start_point,end_point, options);title(<span class="string">'Shortest path'</span>);</pre><img vspace="5" hspace="5" src="content_02.png"> <p>One can also use several starting and ending points and perform progressive propagation</p><pre class="codeinput">[M,W] = load_potential_map(<span class="string">'mountain'</span>, n);<span class="comment">% random seeding of the points</span>nstart = 8; nend = 30;start_points = round( rand(2,nstart)*(n-1)+1 );end_points = round( rand(2,nend)*(n-1)+1 );<span class="comment">% FM computation</span>options.end_points = []; A = {};nlist = round( linspace(0.05,1,8)*n^2 );<span class="keyword">for</span> i=1:length(nlist) options.nb_iter_max = nlist(i); [D,S] = perform_fast_marching(W, start_points, options); mask = repmat(S==-1,[1 1 3]); A{i} = convert_distance_color(D); A{i} = mask.*A{i} + (1-mask) .* repmat(rescale(W),[1 1 3]);<span class="keyword">end</span>clf;imageplot(A, <span class="string">''</span>, 2,4);<span class="comment">% paths extractions</span>gpath = compute_geodesic(D,end_points);<span class="comment">% display</span>clf;subplot(1,2,1);imageplot(W,<span class="string">'Potential map'</span>);subplot(1,2,2);plot_fast_marching_2d(A{end},[],gpath,start_points,end_points, options);title(<span class="string">'Shortest paths'</span>);</pre><img vspace="5" hspace="5" src="content_03.png"> <h2>3D volumetric dataset<a name="7"></a></h2> <p>You can load a 3D volume of data</p><pre class="codeinput"><span class="comment">% load the whole volume</span>load <span class="string">../toolbox_fast_marching_data/brain1-crop-256.mat</span><span class="comment">% crop to retain only the central part</span>n = 100;M = rescale( crop(M,n) );<span class="comment">% display some horizontal slices</span>slices = round(linspace(10,n-10,6));Mlist = mat2cell( M(:,:,slices), n, n, ones(6,1));clf; imageplot(Mlist);</pre><img vspace="5" hspace="5" src="content_04.png"> <p>You can also perform a volumetric display. By changing the value of <tt>options.center</tt>, you can display various level sets. </p><pre class="codeinput">clf;options.center = .35;imageplot(M,options);title(<span class="string">'center=.35'</span>);</pre><img vspace="5" hspace="5" src="content_05.png"> <h2>3D Fast Marching<a name="10"></a></h2> <p>The procedure for volumetric FM is the same as for 2D FM.</p> <p>We ask to the user for some input points.</p><pre class="codeinput">delta = 5;start_point = [91;15;delta];<span class="comment">% You can use instead the function pick_start_end_point</span></pre><p>We compute a potential that is high only very close to the value of <tt>M</tt> at the selected point. </p><pre class="codeinput">W = abs(M-M(start_point(1),start_point(2),start_point(3)));W = rescale(-W,.001,1);<span class="comment">% perform the front propagation</span>options.nb_iter_max = Inf;[D,S] = perform_fast_marching(W, start_point, options);<span class="comment">% display the distance map</span>D1 = rescale(D); D1(D>.7) = 0;clf; imageplot(D1,options);alphamap(<span class="string">'rampup'</span>);colormap <span class="string">jet(256)</span>;</pre><img vspace="5" hspace="5" src="content_06.png"> <h2>Fast Marching on 3D meshes<a name="13"></a></h2> <p>The FM algorithm also works on triangulated meshes. The metric is the metric induced by the embedding of the 3D surface in R^3, but one can also modulate this metric by a scalar metric map. </p> <p>The triangulated mesh is loaded from a <tt>.off</tt> file. You can see the toolbox mesh for more function for mesh processing. </p><pre class="codeinput"><span class="comment">% load the mesh</span>[vertex,faces] = read_mesh(<span class="string">'elephant-50kv'</span>);<span class="comment">% display the mesh</span>clf;plot_mesh(vertex, faces);shading <span class="string">interp</span>;</pre><img vspace="5" hspace="5" src="content_07.png"> <p>Select some starting point and do the propagation</p><pre class="codeinput">start_points = 20361;[D,S,Q] = perform_fast_marching_mesh(vertex, faces, start_points, options);</pre><p>Extract some geodesics and display the result.</p><pre class="codeinput">npaths = 30; nverts = size(vertex,2);<span class="comment">% select some points that are far enough from the starting point</span>[tmp,I] = sort( D(:) ); I = I(end:-1:1); I = I(1:round(nverts*1));end_points = floor( rand(npaths,1)*(length(I)-1) )+1;end_points = I(end_points);<span class="comment">% precompute some usefull information about the mesh</span>options.v2v = compute_vertex_ring(faces);options.e2f = compute_edge_face_ring(faces);<span class="comment">% extract the geodesics</span>options.method = <span class="string">'continuous'</span>;options.verb = 0;paths = compute_geodesic_mesh(D,vertex,faces, end_points, options);<span class="comment">% display</span>options.colorfx = <span class="string">'equalize'</span>;plot_fast_marching_mesh(vertex,faces, D, paths, options);shading <span class="string">interp</span>;</pre><img vspace="5" hspace="5" src="content_08.png"> <h2>Fast Marching inside a 2D shape<a name="17"></a></h2> <p>It is possible to retrict the propagation inside a 2D shape.</p> <p>First we load a 2D shape, which is a BW image</p><pre class="codeinput">n = 128;M = rescale( load_image(<span class="string">'chicken'</span>,n), 0,1 );M = double(M>0.5);</pre><p>Compute geodesic distance to a point inside the shape</p><pre class="codeinput">start_points = [65;65];<span class="comment">% use constant metric</span>W = ones(n);<span class="comment">% create a mask to restrict propagation</span>L = zeros(n)-Inf; L(M==1) = +Inf;options.constraint_map = L;<span class="comment">% do the FM computation</span>[D,S,Q] = perform_fast_marching(W, start_points, options);</pre><p>Select a set of end points, here we locate them on the boundary of the shape.</p><pre class="codeinput">bound = compute_shape_boundary(M);nbound = size(bound,1); npaths = 40;sel = round(linspace(1,nbound+1,npaths+1)); sel(end) = [];end_points = bound(sel,:)';</pre><p>Extract the geodesics from these ending points</p><pre class="codeinput">D1 = D; D1(M==0) = 1e9;paths = compute_geodesic(D1,end_points);</pre><p>Do a nice display</p><pre class="codeinput">ms = 30; lw = 3; <span class="comment">% size for plot</span><span class="comment">% display</span>A = convert_distance_color(D);clf; hold <span class="string">on</span>;imageplot(A); axis <span class="string">image</span>; axis <span class="string">off</span>;<span class="keyword">for</span> i=1:npaths h = plot( paths{i}(2,:), paths{i}(1,:), <span class="string">'k'</span> ); set(h, <span class="string">'LineWidth'</span>, lw); h = plot(end_points(2,i),end_points(1,i), <span class="string">'.b'</span>); set(h, <span class="string">'MarkerSize'</span>, ms);<span class="keyword">end</span>h = plot(start_points(2),start_points(1), <span class="string">'.r'</span>);set(h, <span class="string">'MarkerSize'</span>, ms);hold <span class="string">off</span>;colormap <span class="string">jet(256)</span>;axis <span class="string">ij</span>;</pre><img vspace="5" hspace="5" src="content_09.png"> <h2>Farthest point sampling<a name="23"></a></h2> <p>One can sample an image using a greedy algorithm that distribute, at each step, a new point that is located the farthest away from the previously sampled set of points. </p><pre class="codeinput">clear <span class="string">options</span>;n = 400;[M,W] = load_potential_map(<span class="string">'mountain'</span>, n);npoints_list = round(linspace(20,200,6));landmark = [];options.verb = 0;ms = 15;clf;<span class="keyword">for</span> i=1:length(npoints_list) nbr_landmarks = npoints_list(i); landmark = perform_farthest_point_sampling( W, landmark, nbr_landmarks-size(landmark,2), options ); <span class="comment">% compute the associated triangulation</span> [D,Z,Q] = perform_fast_marching(W, landmark); <span class="comment">% display sampling and distance function</span> D = perform_histogram_equalization(D,linspace(0,1,n^2)); subplot(2,3,i); hold <span class="string">on</span>; imageplot(D'); plot(landmark(1,:), landmark(2,:), <span class="string">'r.'</span>, <span class="string">'MarkerSize'</span>, ms); title([num2str(nbr_landmarks) <span class="string">' points'</span>]); hold <span class="string">off</span>; colormap <span class="string">jet(256)</span>;<span class="keyword">end</span></pre><img vspace="5" hspace="5" src="content_10.png"> <h2>Anisotropic Fast Marching<a name="24"></a></h2> <p>The toolbox also implements anisotropic FM, where at each location, the metric is given by a tensor field.</p> <p>We first create a 2D vector field. The tensor field will be aligned with this tensor field.</p><pre class="codeinput">n = 200;U = randn(n,n,2);options.bound = <span class="string">'per'</span>;U = perform_vf_normalization( perform_blurring(U, 30,options) );</pre><p>Create a set of ending points, located on the boundary of the image.</p><pre class="codeinput">start_points = [n;n]/2;<span class="comment">% end points on the boundary</span>t = 1:n;x = [t t*0+n t(end-1:-1:1) t*0+1];y = [t*0+1 t t*0+n t(end-1:-1:1)];npaths = 50;s = round(linspace( 1,length(x), npaths+1) ); s(end) = [];end_points = cat(1, x(s),y(s));</pre><p>We build a metric with an increasing level of anisotropy, and perform the FM each time.</p><pre class="codeinput"><span class="comment">% test for various degree of anisotropy</span>aniso_list = [.01 .05 .1 .2 .5 1];ms = 30; lw = 3; <span class="comment">% display params</span>clf;<span class="keyword">for</span> ianiso = 1:length(aniso_list) <span class="comment">% build the tensor field</span> aniso = aniso_list(ianiso); V = cat(3, -U(:,:,2), U(:,:,1)); <span class="comment">% orthogonal vector</span> T = perform_tensor_recomp(U,V, ones(n),ones(n)*aniso ); <span class="comment">% propagation</span> [D,S,Q] = perform_fast_marching(T, start_points); <span class="comment">% for sexy display</span> D1 = perform_histogram_equalization(D, linspace(0,1,n^2)); <span class="comment">% extract tons of geodesics</span> paths = compute_geodesic(D,end_points, options); <span class="comment">% display</span> subplot(2,3,ianiso); hold <span class="string">on</span>; imageplot(D1); axis <span class="string">image</span>; axis <span class="string">off</span>; colormap <span class="string">jet(256)</span>; title([<span class="string">'Anisotropy='</span> num2str(1-aniso)]); <span class="keyword">for</span> i=1:npaths end_point = end_points(:,i); h = plot( paths{i}(2,:), paths{i}(1,:), <span class="string">'k'</span> ); set(h, <span class="string">'LineWidth'</span>, lw); h = plot(end_point(2),end_point(1), <span class="string">'.b'</span>); set(h, <span class="string">'MarkerSize'</span>, ms); <span class="keyword">end</span> h = plot(start_points(2),start_points(1), <span class="string">'.r'</span>); set(h, <span class="string">'MarkerSize'</span>, ms); hold <span class="string">off</span>; colormap <span class="string">jet(256)</span>; axis <span class="string">ij</span>;<span class="keyword">end</span></pre><img vspace="5" hspace="5" src="content_11.png"> <p class="footer"><br> Copyright ® 2008 Gabriel Peyre<br></p> </div> <!--##### SOURCE BEGIN #####%% Toolbox Fast Marching - A toolbox for Fast Marching and level sets computations%% Copyright (c) 2008 Gabriel Peyre%%%% First add some directories to the pathpath(path, 'toolbox/');path(path, 'data/');%% 2D Fast Marching Computations% The toolbox allows to compute a distance map from a set of starting point% using an arbitrary isotropic metric. This distance map is then used to% compute geodesic paths from any point that joins the closest starting point,% according to the geodesic distance.%%% We first load a 2D map, and ask to the user for the starting points.n = 400;[M,W] = load_potential_map('road2', n);start_point = [16;219];end_point = [394;192];% You can use instead the function% [start_point,end_point] = pick_start_end_point(M);%%% We then perform the front propagationclear options;options.nb_iter_max = Inf;options.end_points = end_point; % stop propagation when end point reached[D,S] = perform_fast_marching(W, start_point, options);% nicde color for display⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -