File Exchange

Anisotropic Diffusion (Perona & Malik)

version 1.0.0.0 (8.01 KB) by Daniel Lopes

Daniel Lopes (view profile)

A set of filters that perform 1D, 2D and 3D conventional anisotropic diffusion (gray scale data).

Updated 16 May 2007

Anisotropic diffusion is a powerful image enhancer and restorer based on the PDE of heat transfer.

The implementation details are described in "P. Perona and J. Malik, Scale-Space and Edge Detection Using Anisotropic Diffusion, IEEE Transactions on Pattern Analysis and Machine Intelligence, 12(7):629-639, July 1990" and in "G. Grieg, O. Kubler, R. Kikinis, and F. A. Jolesz, Nonlinear Anisotropic Filtering of MRI Data, IEEE Transactions on Medical Imaging, 11(2):221-232, June 1992".

The algorithm was implemented to perform upon 1D, 2D and 3D gray-scale signals.

Cite As

Daniel Lopes (2019). Anisotropic Diffusion (Perona & Malik) (https://www.mathworks.com/matlabcentral/fileexchange/14995-anisotropic-diffusion-perona-malik), MATLAB Central File Exchange. Retrieved .

nbro

BG

BG (view profile)

Can anyone confirm or disprove this?

anisodiff1D has 2 neighbors.

Max delta_t is 1/(1 + 1) = 1/2.

anisodiff2D has 8 neighbors with two weights

Max delta_t is 1/(1 + 1 + 1 + 1 + 1/2 + 1/2 + 1/2 + 1/2) = 1/6

anisodiff3D has 26 neighbors with variable weights depending on voxel dimensions.

Max delta_t depends on the voxel dimensions. It is not exactly 3/44.

For isometric weights, I get 1/(2+2+2+8+16+24) = 1/54. The dimension dependent range goes from 0 to \inf. Consider the voxel dimensions very large or very small. These can be lower or higher than the algorithm's specified max value. This suggests that for many cases the specified max value will lead to unstable conditions.

I started with Peter Kovesi's anisodiff. His 2D implementation has 4 neighbors.

The max value is 1(1+1+1+1) = 1/4.

His [code](https://www.peterkovesi.com/matlabfns/Spatial/anisodiff.m) confirms this value.

Efe Akmansu

Christian Marås

Christian Marås (view profile)

Sorry I meant since it uses 8 diffusion directions*

Christian Marås

Christian Marås (view profile)

Shouldn't the delta variable be between 0< delta<1/8 since this implementation uses 4 diffusion directions?

Tyou Syouten

Tyou Syouten (view profile)

The edge is not well preserved in this implementation. "imshow(ad,[])" to display the output image.

RAVI KUMAR

RAVI KUMAR (view profile)

how to retrieve the 2D image after diffusion??

phzb

phzb (view profile)

Please help: any idea how the 1/7 maximum delta_t is calculated? Want to extend this to bigger dimension. Thanks!

A-J KHAN

A-J KHAN (view profile)

Hi dears, i used conventional anisotropic diffusion process to reduce isolated papillary muscle and intensity heterogeneity using LEFT VENTRICLE CARDIAC MRI images, plz some one help me in this stage.

just yy

Sebastian Ubalde

Sebastian Ubalde

Sebastian Ubalde (view profile)

Nice contribution.

For those using anisodiff1D, consider replacing imfilter(.,.,'conv') by conv(.,.,'same'). It's quite faster.

Royi Avital

Royi Avital (view profile)

I think you only applied the Divergence operator once when you should apply it twice.

Stefan Karlsson

Stefan Karlsson (view profile)

There is no productive 2D example provided that show the superiority of this method over regular gaussian smoothing. A simple combination of Gaussian and Median filtering provides better denoising than what is proposed:

s = phantom(512) + randn(512);
num_iter = 55;
delta_t = 1/7;
kappa = 30;
option = 2;

It is very likely that the authors generalization of using 8 point neighbourhoods is valid, but if so it changes the valid range of the parameters, and thus one cant easily recreate the results of the original Perona-Malik paper. If a simple example was provided that showed how the 2D method actually outperforms gaussian blurring, then I might be inclined to read the secondary paper referenced.

the interpretation of the kappa parameter is nicely given in weickerts book: http://www.lpi.tel.uva.es/muitic/pim/docus/anisotropic_diffusion.pdf, chapter 1, figure 1.1, b... (its called lambda there)

I cant get that interpretation to work for me here. Whats wrong? How can i strengthen the edges in the image using this, while blurring along isolines!?

Again, if only a positive example was provided for the 2D case, this might all be cleared up.

SHANYUE

jia

jia (view profile)

very clear and have many details about the algorithm

GIRE

GIRE (view profile)

Reena Singh Singh

Reena Singh Singh (view profile)

I ve tried the code it is very useful bt i want to generate diffusion function to calculate flux function hw to do?

Przem

Przem (view profile)

makes more sense. it is the actual 4 point PDE as propsed in the perona-malik paper

Reena Singh Singh

Reena Singh Singh (view profile)

thanx fr the code it is easy to understand bt in the perona-malik paper divergence is taken for the diffusion process nd in dis code diff-image is calculated in some different way. can u plz explain m dis?

Andre

Andre (view profile)

Arif Harmanci, if you input a column vector as required by the code you should get the expected behaviour. (in your example, just changing a_n to a_n' in the function call is sufficient).

George Ntokas

George Ntokas (view profile)

I have a question:at the example in Anisodiff2D you wrote this line: s = phantom(512) + randn(512).The number 512 is about memory wright?But what exactly?And how i can use another image?Thank you in advance.

Jiucheng Nie

Jiucheng Nie (view profile)

for the 2D case: why all the flux from 8 directions are summed up? The paper suggests to use flus(north)-flus(south) and so on. There at least is 'minus' instead 'plus' all the way around. Could you explain me why?

Arif Harmanci

Arif Harmanci (view profile)

I was looking at this code and I am pretty sure that there is a problem with the 1D code. I did not try 2d code. I tried following:

% Generate random regions.
a = [zeros(1, 2000) ones(1, 1000) zeros(1, 3000) 5 * ones(1, 4000) zeros(1, 6000) 3 * ones(1, 5000) zeros(1, 4000)];

a_n = a + rand(1, length(a));

diff_sig = anisodiff1D(a_n, 1000, 1/3, 30, 2);

Then I get all 0's, that is, everything is smoothed.

If you go back to Kovesi's code, he did not use imfilter or conv, he computed the differences by basic vector subtractions. When I replace the imfilter to vector subtractions to compute nablaW/E, then I get smoothing only at the noisy regions, which is the expected result.

Would be great if anyone could comment on this.

lina

lina (view profile)

Do you know how to change this function to be applicable to convert in C code using matlab coder??

Tripp

Tripp (view profile)

Very useful and powerful filter

Xianming Liu

Xianming Liu (view profile)

Many thanks. That's just what I want.

Faraz

Faraz (view profile)

Many thanks, but it gives me a blur image of a mamogram instead of enchancing the contrast?
can you/anyone comment y?

Anand Muglikar

Anand Muglikar (view profile)

This is an awesome file that I used for driving out the noise in blood-slides and preserving the cell boundaries to detect malaria parasite.

matlab user

matlab user (view profile)

thank you !

one question is that I can't duplicate the fig 3 in
G. Grieg, O. Kubler, R. Kikinis, and F. A. Jolesz, Nonlinear Anisotropic Filtering of MRI Data, IEEE Transactions on Medical Imaging, 11(2):221-232, June 1992".

the fig 3 show iterative edge sharpening, but my simulation gives me reverse result... can anyone give me a sample code for results similar to fig 3?

Nitin

Nitin (view profile)

Thanks for sharing

Joshua Groves

Joshua Groves (view profile)

Thanks for this. I translated your code and created a Python version that I will hopefully be using for image analysis.

Forrester Cole

Forrester Cole (view profile)

Nice work. Though I think the examples could be a bit more helpful. For example, setting KAPPA = 0.2 or so, rather than KAPPA = 30, in the examples would show the anisotropic behavior more clearly. It took me a while to figure out that I needed to bring KAPPA way down from 30 before I saw results similar to Fig. 11 in the Perona and Malik paper.

Forrester Cole

Forrester Cole (view profile)

Solaree Shi

Thanks for it!

Thanks for your code. it is very clear and useful.

Daniel San Martin Pascal Filho

Daniel meu chará, você me quebrou um galhão com essa implementação. Obrigado.

Andrew Rein

Thanks for the code! it is very clearly written. For the 2D case in your code it is defined dx,dy and dd. What's the purpose of this? I didn't see it in the paper (have only read Perona). Also delta (lambda in the paper) suggest up to 1/4. Is it because you're using extra neighbours for the implementation?

Daniel Mark

Thank you very much for this excellent code.

albert j.

ivan scardanzan

Anisotropic Diffusion
File Id: 3710 Average rating: 4.62
Size: 142 KB # of reviews: 26
Subscribers: 8
Keywords: NonLinear Diffusion, Anisotropic Diffusion, Coherence-Enhancing Diffus

jianfei ge

Thanks!

Bruno Carozza

Thanks, works perferctly!

Daniel Lopes

Anisotropic diffusion is known to:
- smooth the signal;
- preserve strong edges;
- enhance the contrast of edges.

Contrary to spatially filtering the signal with a Gaussian, anisotropic diffusion doesn't smooth across the edges.

A recommend the following procedures to clarify your doubt relatively to the anisotropic implementation:

- run the example of anisodiff1D(.) [the visual effect of anisotropic filtering is more evident when applied to a 1D signal]

- instead of using a Gaussian with a small sigma try a higher one, e.g.,
hG = fspecial('gaussian',31,9); % sigma = 9;
Who will see that even if you increase the number of diffusion iterations the edges will prevail, contrary to the smoothing provided by the Gaussian filter.

Another advice: try reading and understanding the code with greater attention before rating it. The references to the articles aren´t there for nothing you now! :P

David Gavilan

Is this really anisotropic? In file anisodiff2D.m there's an example using s = phantom(512) + randn(512); ad1 = anisodiff2D(s,15,1/7,30,2); Now try with a simple gaussian: hG = fspecial('gaussian',31,3); ad2 = imfilter(s,hG,'conv','symmetric'); Compare ad1 and ad2 and tell me if there's any difference.

MATLAB Release Compatibility
Created with R14SP2
Compatible with any release
Platform Compatibility
Windows macOS Linux