Add Perona–Malik anisotropic diffusion algorithm

[simmplecoder]

Description

Implement anisotropic diffusion as defined by Perona-Malik equation.

There are quite a few variations of the solution for the PDE, so I would like to implement the algorithm in the following way:

template <typename NablaStrategy, typename DiffusivityStrategy, typename ReductionStrategy>
struct diffuse_strategy { 
    NablaStrategy nabla_strategy;
    DiffusivityStrategy diffusivity_strategy;
    ReductionStrategy reduction_strategy;
    diffuse_strategy(<each strategy here>);
    bool needs_padding;
    template <typename View>
    typename View::value_type operator()(View view, point_t current_point) {
        auto nabla = nabla_strategy(view, current_point);
        auto diffusivity = diffusivity_strategy(nabla);
        return reduction_strategy(nabla, diffusivity);
   }
};

With that, we can have defines for the known ones:

using opencv_strategy = ...;
using classic_strategy = ...; // the one in the paper

Then just have

template <typename InputView, typename OutputView, typename ComputationStrategy = opencv_strategy>
void anisotropic_diffusion(InputView input_view, OutputView output_view, ComputationStrategy strategy, uint64_t iteration_count)
{
    // pad if the strategy needs padding
    // initialize two buffers
    // copy input into the first buffer
   <for iteration count>
   {
       <for every x y>
            <second_buffer(x, y) = strategy(first_buffer(x, y))>
       swap(first_buffer, second_buffer);
   }
   copy_pixels(first_buffer, output_view);
}

Note that kappa and delta_t (diffusion speed) will be encoded inside the strategy, thus there no need for more arguments. I believe it will be possible to call it like the following:

anisotropic_diffusion(input, output, {kappa, delta_t}, iter_count);

In the worst case I will just provide an overload.

References

• Scale-Space and Edge Detection Using Anisotropic Diffusion, Pietro Perona, Jitendra Malik, 1990
• http://www.mat.unimi.it/users/naldi/perona-malik.pdf
• http://image.diku.dk/imagecanon/material/PeronaMalik1990.pdf

Tasklist

• Implement the algorithm
• Add docs
• Add example
• Add test case(s)
• Ensure all CI builds pass
• Review and approve

[mloskot]

@simmplecoder Is this draft or ready for review?

For review, I'd suggest you to add tests which reviewers can run. The examples are useful too, but they don't verify results as well as tests.

Any particular reason you added the implementation to image_processing/numeric.hpp?
I'd suggest to add new image_processing/diffusion.hpp file.

[simmplecoder]

@mloskot , the algorithm itself is ready, e.g. the interface is the one I want. I will deal with the CI failing tonight.

About testing: it is a bit complicated. We can do simple sanity check plus some precomputed examples. The source can be matlab scripts presented by the authors, but I'm not sure about the licensing issues. Will it be okay to use the output of it given we do not own matlab commercial license (I have student one) and our users will not have the license? I could search for other sources but the one from paper would be ideal.

Also, the performance is faster than what I thought. I believe it is reasonable for a single threaded CPU implementation.

About the location: I'm not aware of any other diffusion algorithms that are widely used. Though I believe it is reasonable to put there, will move it tonight.

[mloskot]

About testing: it is a bit complicated.
We can do simple sanity check plus some precomputed examples.

Right, so we are getting to the point when our tests need to be based on comparing whole or part of output image against reference image.

I knew I should have implemented the PNM ASCII support done #359 long time ago, actually, for #353 during GSoC 2019. Having this, we could simply write algorithms output to PAM file, then we could do simple diff against corresponding reference file also in PAM.
I do believe test images in human-friendly text form, without using any intermediate image format libraries, would help a lot for testing. I'll see if I can work it out over the weekend.

About the location: I'm not aware of any other diffusion algorithms that are widely used.

Even for single one, I think it's a good idea reflect some categories in structure of the source files.
Second, it helps to keep files of reasonable size, instead of just few gigantic headers which are difficult to maintain.

[mloskot]

First nit-picking review :-)

[simmplecoder]

@mloskot, after working out the equations for tests, I found that not handling the border and corner cases makes this implementation broken for sufficiently high kappa and iteration counts. I will need to handle those cases, but the code will be qutie a bit more brittle, and possibly way larger.

My plan is to introduce a function that will handle the edge cases in a transparent way, but that might not be possible. I will report back once I figure it out (no more than one day should be needed).

[mloskot]

@simmplecoder Your last comment suggest this is still work in progress or is it ready for review?

[simmplecoder]

@mloskot , I've been on a search for a bug, and I believe I found it. The issue is that I'm using completely wrong formula there. I thought that the version that I wrote is actually an improved version of the original, but I was wrong. First of all the summation step should be sum between differences of opposing directions. Second, the flux (which is the product in the summation) step is not diffusivity by difference, but rather diffusivity by double difference. I will check this out and report back.

[simmplecoder]

@mloskot , I'm a bit lost on what people call anisotropic diffusion. I implemented the one that is clearly mentioned in the book, but it still leaks luminosity heavily. I will check out opencv and if it behaves the same way, I will just drop the second test. If my implementation is wrong, I will probably put this on a backburner until the later parts of GSoC.

Book I'm using: link. The formula is on page 81 (page 9 of the pdf). I will deal with the CI failures after figure out if my formula is correct.

[mloskot]

@simmplecoder I think OpenCV implementation is based this original paper, http://image.diku.dk/imagecanon/material/PeronaMalik1990.pdf, which probably is the same as the one you followed (https://github.com/opencv/opencv/issues/8361). Different implementations seem to solve the PDE using also different flux functions, e.g. OpenCV allows two alternatives.

I will try to dig this a bit further and look at your code.
Meanwhile, I think it is a good idea to move to other topic for a break, then get back to it with fresh view.

[simmplecoder]

@mloskot , I believe the current version is the correct one according to the paper. I have updated the description though, please have a look.

[mloskot]

@simmplecoder The PR is getting close to be merged, would it be possible to ensure the CI builds pass?

[simmplecoder]

@mloskot , I almost finished everything. The only thing left is to cover laplace functions with tests, and adapt demo to include command line argument to choose between default, classic and matlab versions of anisotropic diffusion. Code is well tested. Even the diffusivity functions that are not part of any configuration are tested so that in case somebody uses them they are not broken.

Overall this is quite a bit of code, please let me know if some parts don't make sense. I guess I'll need to write up good docs pages for this, though we will need to start from the beginning (topic - finite differences and PDEs). I guess we could do docs at the end of GSoC.

[mloskot]

@simmplecoder First, I'm terribly sorry, for the long silence again.
I am catching up with the reviews now, but first I need to work out the merge of develop to master for Boost 1.74 release before the master closes today.

[mloskot]

@simmplecoder Excellent work!

This is a complex topic and, I think, you've managed to very well, despite quite a number of hurdles encountered on your way. The interface for the anisotropic diffusion is as simple as possible, yet extensible via pluggable strategies.

I'm happy to accept this PR.

The only aspects that could be improved is adding some code comments that could help a reader trace the calculation process and find references between the implementation and the Perona-Malik paper.
A piece of documentation would be awesome too.

---

I think, we can assume there are no objections to merging it.
Please, let me know if/when this is ready to merge.

[simmplecoder]

I had paranoia about current version not being correct. Texts mention that conductivity in south direction is negative of north. What they don't mention is that they take different cells, the ones shifted by one in south direction.

Lets say diff is the vertical difference between elements, upper cell minus lower one. If we will shift cells by one and take negative, it will be I[x][y + 1] - I[x][y] which is exactly what current code does. I'm just documenting it here in case somebody comes by with questions and I forget I had a look at this already.

[mloskot]

@simmplecoder Good finding. I think I've seen different approaches to orientation/origin of kernel in some papers.
The freedom of choice seems similarly to the Hough transform where the angles and radius are anchored differently for the parameter space.

[mloskot]

@simmplecoder Finally, I'm happy to approve this PR

I'm still hoping we will be able to add some piece of documentation for this feature,
but that has to happen after the merge.

[simmplecoder]

@mloskot Thank you! I will write down the documentation after we settle #517 . I'll ask somebody who knows and document it. That seems to be the most critical piece.