Introduction:
This toolbox will perform Anisotropic NonLinear Diffusion filtering on a 2D gray/color or 3D image. This filtering will reduce the image noise while preserving the region edges, and also enhancing the edges by smoothing along them.
This is one of the more advanced image enhancement methods available, and also contains HDCS from october 2009. The result looks like an artist painted the image, with clear brush strokes along the image edges and ridges, see screenshot.
My papers about the code are included:
 "Coherence Filtering to Enhance the Mandibular Canal in ConeBeam CT Data", IEEEEMBS Benelux Chapter Symposium, 2009.
 "Optimized Anisotropic Rotational Invariant Diffusion Scheme on ConeBeam CT", MICCAI, 2010
Method:
The basis of the method used is the one introduced by Weickert.
1, Calculate Hessian from every pixel of the Gaussian smoothed input image
2, Gaussian Smooth the Hessian, and calculate its eigenvectors and values (Image edges give large eigenvalues, and the eigenvectors corresponding to those eigenvalues describe the direction of the edge)
3, The eigenvectors are used as diffusion tensor directions. The amplitude of the diffusion in those 3 directions is based on the eigen values and determined by Weickerts equation.
4, An Finite Difference scheme is used to do the diffusion
5, Back to step 1, till a certain diffusion time is reached.
Diffusion schemes:
There are several diffusion schemes available: standard, implicit, nonegative discretization, and also a rotation invariant scheme, and a novel diffusion scheme with new optimized derivatives.
Mex Files:
All 3D files are not only available as Matlab but also as Ccode /MEX files, to increase speed and reduce the amount of memory used. Compile the ccode by executing compile_c_files.m.
Literature (Full list in the included paper):
 Weickert : "A Scheme for CoherenceEnhancing Diffusion Filtering with Optimized Rotation Invariance"
 Mendrik et al, "Noise Reduction in Computed Tomography Scans Using 3D Anisotropic Hybrid Diffusion With Continuous Switch", October 2009
 Weickert : "Anisotropic Diffusion in Image Processing", Thesis 1996
 Laura Fritz : "DiffusionBased Applications for Interactive Medical Image Segmentation"
 Siham Tabik, et al. : "Multiprocessing of Anisotropic Nonlinear Diffusion for filtering 3D image"
Usage:
Read the help of CoherenceFilter, compile the ccode and try the examples in the help.
Please report bugs, successes and questions.
DirkJan Kroon (2020). Image Edge Enhancing Coherence Filter Toolbox (https://www.mathworks.com/matlabcentral/fileexchange/25449imageedgeenhancingcoherencefiltertoolbox), MATLAB Central File Exchange. Retrieved .
1.7.0.0  Added paper, some minor changes in values 

1.6.0.0  Added new Diffusion scheme which uses optimized 5x5 second order derivatives... 

1.5.0.0  LCC compiler fix 

1.4.0.0  Added Adriënne M. Mendrik et al. "Hybrid Diffusion with Continuous Switch (HDCS)", october 2009 IEEE transactions on medical imaging. 

1.3.0.0  Literature 

1.2.0.0  Fixed boundary check bug in backward derivatives. 

1.1.0.0  Tested Linux Ubuntu 
Inspired: Microscopy Image Browser (MIB), Anisotropic Diffusion with Memory based on Speckle Statistics for Ultrasound Images, CSKMorphometrics
Create scripts with code, output, and formatted text in a single executable document.
Hello, I'm experiencing the same problem as Christian Sanders and rui gao both in Matlab 2018b AND 2019B. I already tried to change the C functions as Christian Sanders but the error remains. I don't understand C and i would deeply appreciate if you can help me solve this problem as right now it's impossible to use the filter.
Error using CoherenceFilterStep2D
Requested 2379411882152x1 (17728.0GB) array exceeds maximum array size preference. Creation of arrays greater
than this limit may take a long time and cause MATLAB to become unresponsive. See array size limit or
preference panel for more information.
Error in CoherenceFilter (line 249)
u=CoherenceFilterStep2D(u,Options);
Error in Untitled (line 22)
filter_coherence=CoherenceFilter(image_original);
(My MATLAB version is MATLAB R2018b)
Every time when I call the function which relates to C files, I will always meet the error which said something about memory and something else.
The functions are 'imgaussian.m', derivatives.m' and 'CoherenceFilterStep2D.m'.
If I changed the name of those files, like 'imgaussian_works.m', I can run the example successfully. Only the example  Rotation Invariant  cannot be gained since 'CoherenceFilterStep2D.m' contains nothing.
I just started learning Matlab and no backgrounds about C, any helpful information for me would be appreciated!
I have also met a similar problem with Christian Sanders, and the hint from my MATLAB R2018b is:
Error using CoherenceFilterStep2D
Requested 1103806595329x1 (8224.0GB) array exceeds maximum array size preference. Creation of arrays greater than this limit may take a long time and cause MATLAB to become unresponsive. See array size limit or preference panel for more information.
Some default values in 'Options' are not the same as the description in the comments. Does that mean that the values in the 'Options' result in this undesirable feedback?
Another question is the thing which relates to compiling. If I compile the file, 'compile_c_files.m', directly, I got the hint:
"Error in compile_c_files (line 32)
mex diffusion_scheme_3D_novel_getUpdate v;"
But if I modify the 32nd line with "mex diffusion_scheme_3D_novel_getUpdate.c v", it would be compilied successfully. Just like Christopher's comment on 30 Jan 2014.
For any of the examples (this is the last one of the list), the program is requesting an extensively large amount of memory. I'm using the 2018b version.
Error using imgaussian
Requested 279172874305x11547402067903512641x2483235742355800016 (17179869184.0GB) array
exceeds maximum array size preference. Creation of arrays greater than this limit may
take a long time and cause MATLAB to become unresponsive. See array size limit or
preference panel for more information.
Error in CoherenceFilter>Anisotropic_step3D (line 314)
usigma=imgaussian(u,Options.sigma,4*Options.sigma);
Error in CoherenceFilter (line 258)
u=Anisotropic_step3D(u,Options);
Error in examples (line 43)
JS = CoherenceFilter(V,struct('T',5,'dt',0.15,'Scheme','S','eigenmode',4));
Fixed by changing the following in each c file:
mxCreateNumericArray(3, dims, mxDOUBLE_CLASS, mxREAL) TO> mxCreateNumericArray(3, dimsu_const, mxDOUBLE_CLASS, mxREAL)
I assumed the interface has changed after all these years
Thank you for sharing the codes.
sir i am beginner ,plz tell me how to run this code for 2d fundus image for image enhancement
I keep having an access violation exception whenever I try to run the 2d example. Has anyone else had this problem? I've run the compile script and I don't see any errors.
now the compiler can run how can I apply it on random images ???
good job
Nice toolbox but I have two issues:
First your default T value is 2 not 5 as described in the help text.
Second in the main while loop for 3D you compare the Scheme strings and look for R:
if(strcmpi(Options.Scheme,'R'))
u=CoherenceFilterStep3D(u,Options);
else
u=Anisotropic_step3D(u,Options);
end
in Anisotropic_step3D you then look again for R:
switch lower(Options.Scheme)
case 'r'
u=diffusion_scheme_3D_rotation_invariant(u,Dxx,Dxy,Dxz,Dyy,Dyz,Dzz,Options.dt);
So diffusion_scheme_3D_rotation_invariant will actually never be called.
I don't know if this intentional, since for 2D your first comparison looks very different:
if(strcmpi(Options.Scheme,'R')&&(Options.eigenmode==0)&&(exist('CoherenceFilterStep2D')==3))
u=CoherenceFilterStep2D(u,Options);
else
u=Anisotropic_step2D(u,Options);
end
I'd like to agree with / confirm Michelle's comment. compile_c_files gives an error " mex: no file name given." unless you change "diffusion_scheme_3D_novel_getUpdate" to "diffusion_scheme_3D_novel_getUpdate.c"
Thanks!
Hi,
there's a little typo in compile_c_files.m.
Last file to be compiled has not the .c extension.
This: "diffusion_scheme_3D_novel_getUpdate.c"
Thanks.
solved! it was some mac os problem with gcc!
Hi
I am getting the allowing error…. any help is greatly appreciated. I am using matlab 2012b on Mac os. thanks
Sourin
derivatives.c: In function 'gradient2Dx':
derivatives.c:24: warning: incompatible implicit declaration of builtin function 'malloc'
derivatives.c: In function 'gradient2Dy':
derivatives.c:57: warning: incompatible implicit declaration of builtin function 'malloc'
derivatives.c: In function 'gradient3Dx_double':
derivatives.c:110: warning: incompatible implicit declaration of builtin function 'malloc'
derivatives.c: In function 'gradient3Dx_float':
derivatives.c:195: warning: incompatible implicit declaration of builtin function 'malloc'
derivatives.c: In function 'gradient3Dy_float':
derivatives.c:280: warning: incompatible implicit declaration of builtin function 'malloc'
derivatives.c: In function 'gradient3Dz_float':
derivatives.c:364: warning: incompatible implicit declaration of builtin function 'malloc'
derivatives.c: In function 'gradient3Dy_double':
derivatives.c:449: warning: incompatible implicit declaration of builtin function 'malloc'
derivatives.c: In function 'gradient3Dz_double':
derivatives.c:533: warning: incompatible implicit declaration of builtin function 'malloc'
mex: compile of ' "derivatives.c"' failed.
Error using mex (line 206)
Unable to complete successfully.
Error in compile_c_files (line 8)
mex(filename,'v');
When filtering 3D images I noticed that there is not option to specify the voxel resolution. To me this implies that in creating this function the author assumed that the voxels of the input image have isotropic resolution. Unfortunately this is usually not the case and must be accounted for when discretizing the diffusion equation.
Aside from this, I think this is a great submission. Many thanks!.
The original stencel can be found here:
http://www.scribd.com/doc/52387266/JahneBHandbookofComputerVisionandApplicationsVol2SignalProcessingandPatternRecognition
Hi,
I have found an error in the code.
The 2D "standard" diffusion algorithm may be incorrect. The stencil shown in Weickertt (2002) is not the same as that referenced (from the handbook volume 2). That stencil has both terms in the topleft and bottom left corners negative. I beleive it is a typo in the Weickertt paper unless I am missing something.
Hi,
I have found an error in the code.
The 2D "standard" diffusion algorithm may be incorrect. The stencil shown in Weickertt (2002) is not the same as that referenced (from the handbook volume 2). That stencil has both terms in the topleft and bottom left corners negative. I beleive it is a typo in the Weickertt paper unless I am missing something.
I think it is useful to explain how to tune the parameters and constants. The default configuration oversmooth my image. I have a very thin object in a noisy image that disappeared by the filter.
thank you.
hi, can I request everyone to report runtime of their experiment and size of image volume.
thanks.
where can imresize2 (that's called in hessian) be found?
it is a useful tool. Thanks the author very much!
it is a elaborate tool. it is very usefull to my research. Thanks the author a lot!
To solve the mex compilation error :
on OSX, add the following as line 3 of diffusion_scheme_3D_novel_getUpdate.c :
#include "string.h"
and rerun the mex script
Compliments for this excellent code
I experience the same problem as Anthony Kilburg:
diffusion_scheme_3D_novel_getUpdate.c:126: warning: incompatible implicit declaration of builtin function 'memcpy'
Probably due to OS X or my compiler? (gcc4.2)
Google 'says' I have to include <String.h>, but this does not seems to solve it.
Dirk,
Thanks for sharing, your profession in matlab is very impressive. the package is really nonlinear anisotropic diffusion, not only coherence diffusion, so I strongly suggest changing the name. I found a problem in functions2D/EigenVectors2D.m the eigenvectors v1 v2 corresponds to mu2 and mu1, please check it, it can be a serious bug.
Hello,
I reached an error while attempting to compile the Ccode files. Error message is as follows:
??? Error using ==> mex at 221
Unable to complete successfully.
Error in ==> compile_c_files at 32
mex diffusion_scheme_3D_novel_getUpdate v;
Upon review of the code I downloaded, I found that all mex funtions call for the function name .c v for all the except the function I am having trouble with. Is the call in the function meant to be:
mex diffusion_scheme_3D_novel_getUpdate.c v;
Thank you
in diffusion_scheme_3D_standard.c , /* Compute tensordriven diffusion (as in [1] pp. 8082) */ which Literature is quoted ? What's [1] ?
hi sir,
i just downloaded ur work. i find an error on line 218 . it says
Error in ==> CoherenceFilter at 218
if(size(u,3)<4), u=double(u); else u=single(u);
A great code for studying diffusion filtering implementation. Thanks a lot for sharing.
I have two questions about the scheme
1) implicit scheme : does it means using semiimplicit AOS scheme to implementation?
2) diffusion_scheme_3D_implicit.m , you said "!! Scheme is unstable, and not ready to use yet." what is the problem?
I uploaded a new version with derivative boundary check bug fixed, probably tomorrow online. thnx Elmar.
*Elmar
Thank you for sharing the solution for your NaN problem.
Now, I know it's probably caused by a boundary check problem.
DirkJan
I permuted the directions of the volume so that it is 147x227x227, now it works perfectly
* Elmar,
Thank you for your report, I never experienced NaN's my self.
Maybe you can add a small value to rule out division by zero.
If you still have problems, maybe you can share a part of the volume data with me?
Thanx, DirkJan
Thanks for sharing this tool! works very well and fast, but i'm experiencing some problems:
I'm trying to filter a 3D MRI volume of size 227x227x147. Depending on the value of 'rho' I use I get a variable number of slices in the third dimension that are all 'NaN'. The bigger rho, the more NaN slices are in the output image volume. For every value of T used the number of NaNslices accumulates, so that after several timesteps the whole image volume consists of NaNs.
What could be the reason for the NaNs in the output?
My original image is single and normalized to a range of 0 to 1 and does not contain any NaNs.
thx, Elmar