version 1.12.0.0 (121 KB) by
Dirk-Jan Kroon

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

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

Dirk-Jan Kroon (2021). Accurate Fast Marching (https://www.mathworks.com/matlabcentral/fileexchange/24531-accurate-fast-marching), MATLAB Central File Exchange. Retrieved .

Created with
R2009a

Compatible with any release

**Inspired:**
Microscopy Image Browser (MIB), Microscopy Image Browser 2 (MIB2)

More Files in the Power Electronics Control Community

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!Create scripts with code, output, and formatted text in a single executable document.

Peshala ThibbotuwawaLandon ChadThe paper by Hassouna (2007) discusses the angles between the derivatives within the stencils. I understand that your 2D code doesn't need to consider the angles because both of the stencils contain orthogonal derivative directions. However, the last 2 stencils in your 3D code do NOT contain orthogonal derivatives, so shouldn't you have to account for the angles between the derivatives?

Ali Özgür ArgunsahDeleting all __inline and iszero function in common.c fixed issues with ubuntu 18.04 gcc 7.5

ROY EL ZEGHONDIHi ,

Im using the skeleton function to get the centerline of vessels i have with their be=ranches , it works really well however it is missing some branches.

I am getting this message (The shortest path trace did not finish at the source point) and I think it is the cause of the missing branches .

Can someone confirm if this is the issue and how can i solve it to get my missing branches and maybe elaborate on the meaning of the message ?

Thank you

brskaFollowing up on my last comment, I was able to resolve the issue. It turned out to just be a problem in the coordinates. Be careful if using the mex compiled files, any "internal errors" might be your actually be your own mistakes but you won't get a useful error message. I was able to resolve the issue by running it on a clean install without the mex compiled files which yielded the proper error message.

brskaHaving a similar issue to some people below but with Matlab 2020. I get the error: "MATLAB has encountered an internal problem and needs to close." I have already tried to build all the c files with "mex -compatibleArrayDims" but the error persists.

PseudomonasAnyone compiled the Mex files with Matlab 2019a on Ubuntu 18.04 (gcc 7.4.0) by chance?

mr janemcdull hallIt might be necessary to compile the c functions with mex in compatible mode, for example:

mex -compatibleArrayDims rk4.c

Rockson RocksonProf. Kroon,

I am a researcher at the Academy of Mathematics and Systems Science, Chinese Academy of Sciences.

I reimplemented numerical experiences of the paper, Multistencils Fast Marching Methods:A Highly Accurate Solution to the Eikonal Equation on Cartesian Domains, through your code. When I used your program to test the convergence order of the MSFM2 method, I found that it did not reach the second-order convergence and only reached the first-order convergence. In other words, I used your code to reimplement the result of TABLE 10 of the reference paper, and there occured a problem.

I don't know where the problem is. In the attached file, the m file is my test code, and the pdf file is the reference paper for your program.

Thanks. I am looking forward to your reply eagerly. I can hardly find your email address on the internet. You can reply to my E-mail, my E-mail address is 'lusong@lsec.cc.ac.cn' .

Best wishes,

Song Lu,

Institution: LSEC, AMSS, CAS

Javier E. SantosHello all,

I'm using matlab 2017a and whenever I run:

SourcePoint=rand(2,100)*255+1;

SpeedImage = ones([256 256]);

tic; T1_MSFM2 = msfm(SpeedImage, SourcePoint, true, true); toc;

figure, imshow(T1_MSFM2,[]); colormap(hot(256));

I get back a binary image with -1s and 0s.

Anyone knows how to fix this?

Noah MitchellWhen running shortestpath.m, MATLAB crashes: "MATLAB has encountered an internal problem and needs to close." I tried deleting "__inline" throughout rk4.c, and I also tried replacing "__inline" with "static inline". No change in outcome. Any suggestions?

Eric ChadwickI've gotten some very nice results with this method, however I have noticed that many smaller features of my image are ignored and do not receive skeleton lines. I have a scan of a full mouse lung so I go from vessels of ~1000 um in diameter to about 15 um. These smaller ones seem to not have centrelines. Is there any parameter I could change to improve this?

Michel KocherHello, I have a problem compiling rk4.c on a Mac architecture. I got this message :Error using mex

Undefined symbols for architecture x86_64:

"_interpgrad3d", referenced from:

_RK4STEP_3D in rk4.o

ld: symbol(s) not found for architecture x86_64

clang: error: linker command failed with exit code 1 (use -v to see invocation)

Error in compile_c_files (line 15)

mex -compatibleArrayDims rk4.c

Do you have any advice ?

renaissanceMy bad. I should have compiled the code first.

renaissanceHi, I tried to run the example inside skeleton.m, but the script failed with the following errors. Please advise how to fix it. Thank you.

***********************************************

I=im2double(rgb2gray(imread('images/vessels2d.png')));

% Convert double image to logical

Ibin=I<0.5;

% Use fastmarching to find the skeleton

S=skeleton(Ibin);

% Display the skeleton

figure, imshow(Ibin); hold on;

for i=1:length(S)

L=S{i};

plot(L(:,2),L(:,1),'-','Color',rand(1,3));

end

Distance Map Constructed

Find Branches Iterations : 0

Attempt to execute SCRIPT rk4 as a function:

C:\Matlab_work\FastMarching_version3b\shortestpath\rk4.m

Error in shortestpath (line 64)

EndPoint=rk4(StartPoint, GradientVolume, Stepsize);

Error in skeleton (line 97)

ShortestLine=shortestpath(T,StartPoint,SourcePoint,1,'rk4');

Ilaria Olivetisorry, I haven't read the dmi pr's answer. Now the example in shortestpath.m is running incredibly slow. Has anyone else had the same problem?

Jinlong FuHi, please read the comments from 'carter griest'. It may solve your problem about the code.

carter griest@hone Greenhall, take dmi pr's advice and recompile all of the .c files with -compatibleArrayDims like this:

mex -compatibleArrayDims msfm2d.c

mex -compatibleArrayDims msfm3d.c

cd('shortestpath');

mex -compatibleArrayDims rk4.c

yimin zhangFu-Chen Chenjohn greenhallWhen running the example in shortestpath.m, I keep getting the "MATLAB has encountered an internal problem and needs to close." error message. Looking into the details, it looks like there's an access violation in the rk4.mexw64 file. I am running MATLAB 2017b on windows 10. Has anyone else run into this issue?

carlos niñoExcellent implementation

I have a question the dx and the dy, they can be varied.

William WarrinerTested on 201 cubed image with uniform speed. Seed was center, i.e. [ 101, 101, 101 ]. The resulting contour at distance=50 was as posted in this image. https://imgur.com/a/59yh2 Expected something closer to a sphere, or at least something with two-fold symmetry over each axial plane. Please investigate so I can rate you up!

William WarrinerTested on 201 cubed image with uniform speed. Seed was center, i.e. [ 101, 101, 101 ]. The result was as posted in this image. https://imgur.com/a/59yh2 Expected something closer to a sphere, or at least something with two-fold symmetry over each axial plane. Please investigate so I can rate you up!

mingzhi chenmingzhi chendmi prTo use with new versions of Matlab just add -compatibleArrayDims when compile a mex files, like this:

mex -compatibleArrayDims msfm2d.c

Martin DenkMartin DenkTomasz NogaAnyone knows how to de-bug the toolbox for Matlab 2017a (as it only returns '-1' with this version of matlab)?

hao liuI tried this for all the night, it executes very slow, the problem is it can't find all the centerline of my vessel, feel so frustrated...

hao liuthanks for the code i have tried to find a code like this in github for a week...what a surprise!!!

Karena CaiIs there a way to refine the resolution of the mesh used to calculate the target distances for the msfm.m function for the 2d image?

Karena CaiKevin Gobron@Enoch Pae I didn't use the functions you tried but all I can tell you is theses scritps doesn't work on Matlab R2017a (I tried). It works properly on Matlab 2015 though (with an appropriate compiler).

Kevin GobronPelegThank you Dirk-Jan Kroon.

I used the ShortestPath on 3D binary data and the script worked successfully.

Enoch PaeI am currently using mac os Sierra 10.12.4 and Matlab R2017a.

Whenever I try to use the skeleton.m on the examples, it gets trapped in the while loop in shortestpath.m

Turns out it gets problem with the rk4 function in

line 64 EndPoint=rk4(StartPoint, GradientVolume, Stepsize);

the EndPoint keeps returning NaN's.

I've erased all the __inline in rk4.c and compiled it.

Anybody else encountered such problem?

Also, is it right that the getBoundaryDistance function is supposed to return a matrix with only -1, 0 values?

line 73 [SourcePoint,maxD]=maxDistancePoint(BoundaryDistance,I,IS3D);

in skeleton.m

Then why bother getting maxD, if it's always 0? just curious.

Thanks

Enoch PaeThanks for the code, msfm works great.

Enoch PaeThanks for the code, by the way, msfm works great.

S SSThank you very much for your submission. It is very useful.

I use skeleton to find the centerline of an aorta. The problem when I have big matrix (512*512*559) it take a lot of time and I have no result. How can I resolve this problem? Thank you

AlnsGRFantastic, clean skeleton in 3D, albeit a bit slow (average 5 min) but worth the time.

Benas SEDIT:

have you tried to find skeleton of the square, which height and width are equal to some odd number, let's say 41? In that case algorithm throws an error:

Output argument "S" (and maybe others) not assigned during call to "skeleton>OrganizeSkeleton".

Temporary solution:

change line

if(linelength<maxD*2), break; end;

to

if(linelength<maxD*1.2), break; end;

int the skeleton.m file

Jim McElwaineYes it only returns -1 with recent versions of matlab

There are a lot of compiler warnings which may be the problem

chang shuwhen I run Example 2D

it retuens T1_FMM1 T1_MSFM1 T1_FMM2 T1_MSFM2, but these value are all -1

what is the problem?

Shengnan LiuJulien SmetsEli YorkWith the Windows 10 updates, there are a lot of issues using mex to compile. Could you please just compile these files into matlab functions and make them available? Thank you!

WouterSome 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

Sandra CaballéAndrew BrownZongbao WangHi All,

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

Best,

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

Ebrahim RasromaniI 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:

"_interpgrad3d", referenced from:

_RK4STEP_3D in rk4.o

ld: symbol(s) not found for architecture x86_64

clang: error: linker command failed with exit code 1 (use -v to see invocation)

zoza azawihi, thank you for shared your code.

but I'm new mat-lab programmer, can anyone tell me how to run this code, please?

mathcodeI 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!

zhang yanhuiHi 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 .

Ilya BelevichEvgeny PrCharles CyizaThis 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 !

WilliamCannot 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);

Please advise, many thanks

jianboSome compiling notes for clang (mac):

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

NidhinHi 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:

"_interpgrad3d", referenced from:

_RK4STEP_3D in rk4.o

ld: symbol(s) not found for architecture x86_64

clang: error: linker command failed with exit code 1 (use -v to see invocation)

Error in compile_c_files (line 12)

mex('rk4.c');

Lowell LiuSeems 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')

Nicolas YuThis 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 ?

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

Burcu GuldurGreat 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?

Are LosnegårdGreat toolbox! Would it be possible to allow for anisotropic voxels? Thanks

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

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

pangyutengExcellent 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.

pangyutengsvetlanaUrs UtzingerI 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.

Mark DrakesmithExcellent. 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.

Mark DrakesmithsvetlanaMohammadHi 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).

Thanks for your help in advance.

ilkay oksuzCan 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..

olivier crosDear 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.

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

Erik ValentiThank 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.

Erik ValentiIs 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

ZaraI 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)?

TieyuanThanks 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.

Baran AydoganHere 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!

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

xiang fionaExcellent 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?

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

Dirk-Jan Kroon*Olivier Roy,

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

Chen ShuoMy 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')

Olivier RoyGreat 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

FlorenceNice 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.

Dirk-Jan Kroon*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.

DanielaThe 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.*

MatteoDoes this work with R2007a? Thank you

NoraNice 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.

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

Sebastien PARISFast and work perfectly. Very usefull for Path Planning