Code covered by the BSD License

### Highlights from Accurate Fast Marching

4.48148
4.5 | 30 ratings Rate this file 87 Downloads (last 30 days) File Size: 121 KB File ID: #24531 Version: 1.12

# Accurate Fast Marching

### Dirk-Jan Kroon (view profile)

23 Jun 2009 (Updated )

Multistencils second order Fast Marching 2D and 3D including rk4 shortest path and skeletonize

File Information
Description

Descriptions:

- The function MSFM2D/MSFM3D calculates the shortest distance from a list of points to all other pixels in an 2D or 3D image, using the Multistencil Fast Marching Method (MSFM). This method gives more accurate distances by using second order derivatives and cross neighbors.

- The function Skeleton will calculate an accurate skeleton (centerlines) of an object represented by an binary image / volume using the fast-marching distance transform.

- The function Shortestpath traces the shortest path from start point to source point using Euler or Runge Kutta 4 in the 2D or 3D distance map.

Implementation:
The 2D fast marching method is implemented in both Matlab-code and c-code. The c-code uses a custom build unsorted binary tree minimum search which out performs a normal binary sorted tree. The c-code is more than 500 times as fast as the matlab-code (compiled with the Microsoft Visual compiler).

Literature:
We used two papers:
- J. Andreas Baerentzen "On the implementation of fast marching methods for 3D lattices"
- M. Sabry Hassouna et al. "Multistencils Fast Marching Methods: A Highly Accurate Solution to the Eikonal Equation on Cartesian Domains"
- R. van Uitert et al. "Subvoxel precise skeletons of volumetric data based on fast marching methods"

We compared the results of our implementation with the results in the paper:
- Normal fast marching 1th order, exact the same results.
- 2th order, significant smaller errors than in the paper.
- Multistencil 1th order, larger errors than in the paper
- Multistencil 2th order, significant worse results than published in the paper. (Note : Our results are in accordance to other existing implementations )

The last version of our code produces better result than in the paper or other literature. This is achieved by solving the polynomial roots using all the available information, as described by the comment of Olivier Roy below.

Examples:
Compile the c-code with mex msfm2d.c; mex msfm3d.c; mex rk4.c;

Try the examples in the help of msfm2d, shortestpath and skeleton

Acknowledgements

This file inspired Microscopy Image Browser (Mib).

MATLAB release MATLAB 7.8 (R2009a)
06 Jul 2016 Wouter

### Wouter (view profile)

Some more compiling notes for Mac OS X 10.11.5 with Matlab 2016a:
- remove all instances of "__inline" in rk4.c
- change all "/* comment styles */" into "// comment" styles in all *.c files

12 Jun 2016 Sandra Caballé

### Sandra Caballé (view profile)

16 Apr 2016 Andrew Brown

### Andrew Brown (view profile)

15 Mar 2016 Zongbao Wang

### Zongbao Wang (view profile)

Hi All,
who has the c++ code of fast marching algorithm? could you please send me?

Best,

07 Mar 2016 ks

### ks (view profile)

@Ebrahim deleting the "__inline" prefix to interpgrad3d on line 76 fixes the problem.

Comment only
28 Dec 2015 Ebrahim Rasromani

### Ebrahim Rasromani (view profile)

I get an error compiling rk4.c (Error output shown below). Any suggestions on how to fix this?

>> mex('rk4.c')
Building with 'Xcode with Clang'.
Error using mex
Undefined symbols for architecture x86_64:
_RK4STEP_3D in rk4.o
clang: error: linker command failed with exit code 1 (use -v to see invocation)

Comment only
16 Dec 2015 zoza azawi

### zoza azawi (view profile)

hi, thank you for shared your code.
but I'm new mat-lab programmer, can anyone tell me how to run this code, please?

30 Sep 2015 mathcode

### mathcode (view profile)

I like this code since it gives each individual branch with coordinates sorted. However, I recently cannot run the 3d matrix (even the size is very small). It hangs over the first iteration forever. I run it before with much larger matrix, and it worked well. Just suddenly stop working, not even the exampled volume. I re downloaded, but no help. Any thoughts? Thanks!

Comment only
30 Sep 2015 zhang yanhui

### zhang yanhui (view profile)

Hi Dirk,

Just found an error in the matlab function file named 'msfm2d'. It is located in the sub-function 'CalculateDistance'. There is a line where you try to calculate the distance using the cross directions, i.e., 'Coeff=Coeff + [...]'. I think 'Coeff' should be reinitialized with a correction as 'Coeff=[...]'. I guess this may be a translation error since I didn't find it in your corresponding C version .

20 May 2015 Ilya Belevich

### Ilya Belevich (view profile)

17 Feb 2015 Evgeny Pr

### Evgeny Pr (view profile)

16 Jan 2015 Charles Cyiza

### Charles Cyiza (view profile)

This it is a good support and can helps a a lot.
But one question: how can I run this code for my 3D vascular tree already extracted with threshold method now I need to extract the centerline?

Thank you !

16 Dec 2014 William

### William (view profile)

Cannot run skeletization in 3D as msfm3d.m not a function resulting in the following error:

Attempt to execute SCRIPT msfm3d as a function:
C:\Users\Will\Desktop\Matlab Tools\functions\msfm3d.m

Error in msfm (line 105)
T=msfm3d(F, SourcePoints, UseSecond, UseCross);

Error in skeleton>getBoundaryDistance (line 201)
BoundaryDistance = msfm(SpeedImage, SourcePoint, false, true);

Error in skeleton (line 66)
BoundaryDistance=getBoundaryDistance(I,IS3D);

Comment only
28 Nov 2014 jianbo

### jianbo (view profile)

Some compiling notes for clang (mac):

change 'inline' to 'static inline' before compiling with mex

Comment only
11 Nov 2014 Nidhin

### Nidhin (view profile)

Hi Dirk,

Thanks a lot for sharing your code! One thing though, quite basic stuff I am afraid, but hopefully you can help: I cannot compile rk4, I get the error below.

Any help will be more than welcomed!

Pablo
Building with 'Xcode with Clang'.
Error using mex
Undefined symbols for architecture x86_64:
_RK4STEP_3D in rk4.o
clang: error: linker command failed with exit code 1 (use -v to see invocation)

Error in compile_c_files (line 12)
mex('rk4.c');

Comment only
06 Nov 2014 Lowell Liu

### Lowell Liu (view profile)

Seems to have error for constant speed field tracking.

SpeedImage=ones(1000,1000)*0+1500;
SourcePoint=[800;400];
DistanceMap= msfm(SpeedImage, SourcePoint);
figure, imshow(DistanceMap,[0 3400])
StartPoint=[100;100];
ShortestLine=shortestpath(DistanceMap,StartPoint,SourcePoint);
hold on, plot(ShortestLine(:,2),ShortestLine(:,1),'r')

17 Apr 2014 Nicolas Yu

### Nicolas Yu (view profile)

This work is good. But I have one suggestion that why not sort these vessel tree branches with bi-tree, so we are easy to show it consistently ?

05 Jun 2013 Ben

### Ben (view profile)

Thanks for sharing! The skeleton runs quite slow when dealing with 3D data. Is there a fast version (e.g. mex)?

14 Mar 2013 Burcu Guldur

### Burcu Guldur (view profile)

Great submission! Thanks.
I have a quick question, I am using skeleton function on a 'C' shaped section (3 segments) but I end up with 2 segments only (1 long and 1 short). Any ideas?

06 Mar 2013 Are Losnegård

### Are Losnegård (view profile)

Great toolbox! Would it be possible to allow for anisotropic voxels? Thanks

16 Feb 2013 hakan

### hakan (view profile)

skeleton function example doesn't work.somebody help me please

16 Feb 2013 hakan

### hakan (view profile)

skeleton function example doesn't work.somebody help me please

29 Dec 2012 ted p teng

### ted p teng (view profile)

Excellent code!
In response Mark Drakesmith's question in April, one can call the shortestpath function once by specifying the StartPoint and SourcePoint in the skeleton function.

Comment only
28 Dec 2012 ted p teng

### ted p teng (view profile)

08 Jul 2012 svetlana

### svetlana (view profile)

19 Jun 2012 Urs Utzinger

### Urs Utzinger (view profile)

I believe by creating groups of the starting points one could use parfor (on those groups) to compute distance maps and take advantage of multiple CPU cores with each core working on a group of the distance points.

Comment only
25 Apr 2012 Mark Drakesmith

### Mark Drakesmith (view profile)

Excellent. Just what I need for skeletonisation of 3D images. Has saved me lots of time!

Just one thing I am wondering... Is it possible to restrict the skeletonisation to a single, most-probable, branch? All my images are of a single branch but, depending on which threshold I use, the algorithm gives me multiple branches.

Many thanks.

25 Apr 2012 Mark Drakesmith

### Mark Drakesmith (view profile)

17 Apr 2012 svetlana

### svetlana (view profile)

Hi every one.
I wrote FMM code in C++ by myself but for large number of grids, I think it take too much time. Would you please tell me know that typically how long it should take for running the efficient code in a mesh that contains 1 million grids.(100*100*100).

12 Dec 2011 ilkay oksuz

### ilkay oksuz (view profile)

Can anyone help me how to use the compile c-files code in 64-bit computers.

I could not compile it on 64-bit. I end up with the following error.

Error using ==> mex at 208
Unable to complete successfully.

On 32-bit the memory of my computer is not sufficient even for 3-D example..

Thanks for the contribution by the way..

Comment only
25 Oct 2011 olivier cros

### olivier cros (view profile)

Dear Dirk-Jan Kroon,

For a sake a fast computation, I ran the example given in the skeleton file with a reduced version of the vessels3d =>V=V(1:128,1:128,1:128);

After final completion, if you display the skeleton along with an isosurface of the volume of interest (V), and rotate the figure, the skeleton seems to connect the three separated parts into one skeleton.

I think there should be some sort of check on whether the skeleton is passing a background area where there is no foreground voxels representing the blood vessel.

Unless I am doing something wrong, and if so, I will immediately apologize to you.

I also discovered a tiny mistake, easily solvable, in the example of your skeleton code (line 42 - commented):

%V = imfill(Vraw > 30000,'holes');

It should be V instead of Vraw.

Otherwise, thank you for sharing this code with us.

Best regards,

Olivier Cros.

Comment only
27 Jul 2011 Matt

### Matt (view profile)

Hello, going along what Erik Valenti said, is there a .m version of the msmf3d.c code? Thanks

Comment only
27 Jun 2011 Erik Valenti

### Erik Valenti (view profile)

Thank you for this code! It is a great help for my research. Has anyone come across problems with running out of memory? If so, have you found a way to fix this problem.

Comment only
27 Jun 2011 Erik Valenti

### Erik Valenti (view profile)

Is there a way to get a MATLAB .m version of the msmf3d.c code? We are trying to keep everything to a MATLAB program if possible. Thanks

Comment only
08 May 2011 Zara

### Zara (view profile)

I find a large difference in my results between using and not using cross-neighbours, and would like to know what exactly the cross-neighbours are that you use - as I can not find a reference for them anywhere.
Are they the diagonal neighbours calculated via directional derivatives (Hassouna 2007)?

Comment only
06 May 2011 Tieyuan

### Tieyuan (view profile)

Thanks for sharing code. It is faster one than what I used for my research.

One problem I found: when I used constant speed, shortestpath code will crash. Also, shortestpath for 3D case (just like 2D since constant value in y direction) is not consistent with 2D case. I can send your my results if you want to take a look at them.

Comment only
18 Jan 2011 Baran Aydogan

### Baran Aydogan (view profile)

Here is a compiler error and the fix:

My compiler (gcc version "4.4.3-4ubuntu5)") gives the following error when compiling "rk4.c":

------------------
rk4.c: In function ‘RK4STEP_2D’:
rk4.c:133: error: expected expression before ‘/’ token

mex: compile of ' "rk4.c"' failed.

??? Error using ==> mex at 222
Unable to complete successfully.

Error in ==> compile_c_files at 12
mex('rk4.c');
----------------------------------

The referred line says:
//double D[2],dl;

FIX:
Just replace the line with:
/*double D[2],dl;*/

Thanks for the great contribution!

Comment only
13 Jan 2011 Marios Karaoulis

### Marios Karaoulis (view profile)

Can you provide the m files of the c code? Some of us do not have matlab compiler.

Comment only
11 Jan 2011 xiang fiona

### xiang fiona (view profile)

Excellent job!
But I have a problem when I test it on my data, local tumor vascular which is unconnected volume.However,The result skeleton was turned out to be connected and passed through two separated vessels.How can I fix this?

03 Jan 2011 Waiwai Wang

### Waiwai Wang (view profile)

03 Jan 2011 Waiwai Wang

### Waiwai Wang (view profile)

Very useful submission，and I have a question that is how long it will be finished based on 3D image during your experiment?

Comment only
13 Dec 2010 Dirk-Jan Kroon

### Dirk-Jan Kroon (view profile)

*Olivier Roy,
Thanks for your very useful suggestion, to increase the accuracy.

Comment only
02 Dec 2010 Chen Shuo

### Chen Shuo (view profile)

My i ask a question: I prceed the codes as following, but the shortestpath is not finded. What is the reason?
My testing code:
SourcePoint = [26; 100]; %Starting point
F = zeros([101 101]);
F(9:27,10:12) = 1;
F(25:27,12:100) = 1;
SpeedImage = F*1000+0.001;
[X Y] = ndgrid(1:101, 1:101);
T1 = sqrt((X-SourcePoint(1)).^2 + (Y-SourcePoint(2)).^2);
tic; T1_MSFM2 = msfm(SpeedImage, SourcePoint, true, true); toc;
figure, imshow(T1_MSFM2,[ ]); colormap gray;
StartPoint=[10;11];
ShortestLine=shortestpath(T1_MSFM2,StartPoint,SourcePoint);
hold on, plot(ShortestLine(:,2),ShortestLine(:,1),'g')

Comment only
01 Dec 2010 Olivier Roy

### Olivier Roy (view profile)

Great submission thanks. The tree structure you use is very efficient.

I compared, as you did, the accuracy of your implementation with the results reported in M. Sabry Hassouna et al. "Multistencils Fast Marching Methods: A Highly Accurate Solution to the Eikonal Equation on Cartesian Domains". Not sure exactly what they do but I found that the accuracy depends a lot on how the results of the two stencils are combined (which is somewhat arbitrary since we do not know a priori which stencil provides the most accurate result).

To improve the accuracy here is the trick: instead of computing the distance with the first and second stencils separately, simply sum the corresponding second order polynomials and solve the resulting second order equation. In other words, simply replace

if(usecross) {
Coeff[0]=0; Coeff[1]=0; Coeff[2]=-1/(max(pow2(Fij),eps));
...

by

if(usecross) {
Coeff[0]+=0; Coeff[1]+=0; Coeff[2]+=-1/(max(pow2(Fij),eps));

in the CalculateDistance function.

With this little modification, I could obtain better results that the aforementioned paper for their Experiment 1.

Note that a weighted sum of the two polynomials can also be done.

Hope this helps,
Olivier

15 Oct 2010 Florence

### Florence (view profile)

Nice work !
But there are prohibitive memory problems when using the C version. The matlab script stops with the message (OUT OF MEMORY) when calling the MSFM2D several times inside the script.

29 Mar 2010 Dirk-Jan Kroon

### Dirk-Jan Kroon (view profile)

*Daniela
The speed function may contain all values larger then zero.
But in practice you have numerical round offs, and because the next "time front" always depend on the old "time front", this error will grow and can sometimes become very high.

Comment only
05 Mar 2010 Daniela

### Daniela (view profile)

The code is easy to understand, but the speed function must be between [1e-8,1], otherwise you get a complex number to value of T, where T is arrival time - this is not correct. Anybody observe/solve that? *Notice that such a constraint for F is theoretically wrong.*

22 Oct 2009 Matteo

### Matteo (view profile)

Does this work with R2007a? Thank you

Comment only
27 Aug 2009 Nora

### Nora (view profile)

Nice submission, works perfectly, and well commented. Thank you!
One thing I was wondering:
Is it possible to give a maximum distance to compute, and not fill the entire image with correct distances? In my problem that would speed up computation considerably.

Furthermore, I deleted the line
if(itt==652221) { printf("569 \n"); }
in msfm2d.c as it put my command window full with 569s.

07 Aug 2009 Petros Mpogiatzis

### Petros Mpogiatzis (view profile)

Great submition, works perfectly, simple to understand how it works, well done!

24 Jun 2009 Sebastien PARIS

### Sebastien PARIS (view profile)

Fast and work perfectly. Very usefull for Path Planning

23 Jun 2009 1.1

text

28 Sep 2009 1.2

Changed c-code comments // to /* */
In help added, that speed function must be larger than zero

01 Oct 2009 1.3

Linux Ubuntu Tested

27 Oct 2010 1.4

Added skeletonize method for vessel centerline extraction

21 Dec 2010 1.8

Centerline in 3D now works robust. More Accurate see comment of Olivier Roy.

14 Jan 2011 1.12

Fixed Crash-bug, and solved bug when curvature and 2e order are enabled. Shortest-path optimized.