Ask Your Question

isotropic non-linear diffusion smoothing (Perona-Malik)

asked 2015-09-16 09:41:58 -0500

theodore gravatar image

updated 2015-09-16 09:51:17 -0500

Now in addition to the my previous thread regarding isotropic linear diffusion smoothing, I want to solve the non-linear version of it based on the Perona-Malik approach. Again we have the following formula:

image description(1)

, where D is not a constant but varies across the image domain. A popular choice is the Perona-Malik diffusivity which can be given as:

image description(2)

where λ is the contrast parameter. ∇u(x, y) is the gradient of the image at pixel (x, y). It seems that KAZE features embed this functionality, therefore I have a look at the source code. Specifically the formula (2) is implemented with the following function:

/* ************************************************************************* */
 * @brief This function computes the Perona and Malik conductivity coefficient g2
 * g2 = 1 / (1 + dL^2 / k^2)
 * @Param Lx First order image derivative in X-direction (horizontal)
 * @Param Ly First order image derivative in Y-direction (vertical)
 * @Param dst Output image
 * @Param k Contrast factor parameter
void pm_g2(const cv::Mat &Lx, const cv::Mat& Ly, cv::Mat& dst, float k) {

    Size sz = Lx.size();
    dst.create(sz, Lx.type());
    float k2inv = 1.0f / (k * k);

    for(int y = 0; y < sz.height; y++) {
        const float *Lx_row = Lx.ptr<float>(y);
        const float *Ly_row = Ly.ptr<float>(y);
        float* dst_row = dst.ptr<float>(y);
        for(int x = 0; x < sz.width; x++) {
            dst_row[x] = 1.0f / (1.0f + ((Lx_row[x] * Lx_row[x] + Ly_row[x] * Ly_row[x]) * k2inv));

However, while I am trying to use it regarding to what @LBerger came up in the other thread I cannot get the correct output. What do I miss again?

I see that the author apply some scalar non-linear diffusion step functionality which I do not really understand what it is about it, here is the function. I tried what @Guanta suggested in the other thread

I think in the evolution_ - Vector (KAZEFeatures.h) are the evolutions over time, so if you'd take the last element and then from that element's struct (TEvolution.h) the Lsmooth Mat should be the image which has been smoothed. To create the evolution_-vector you need to call Create_Nonlinear_Scale_Space() .

but it did not give a good output :-(.

edit retag flag offensive close merge delete



So I haven't studied link given. I have read perona malik paper

I proposed this program

LBerger gravatar imageLBerger ( 2015-09-17 07:26:23 -0500 )edit

@LBerger man I really appreciate your time to help me. I tested you code and it seems to provide the expected result :-). However, if it is not a harsh I would like you to provide some comments in your code in order to understand what you are doing in each step, I guess this would be helpful for other users as well. Moreover, again make it an answer in order to accept it as correct. Again many thanks ;-)

theodore gravatar imagetheodore ( 2015-09-17 10:36:00 -0500 )edit

That's an interesting question and I like giving my time for tah's sort of problem

I don't think that code is ready for an answer. In this paper equation 8 (p 632) it's OK but equation (9) and (10) I m' not sure. I have inserted some comment in code

LBerger gravatar imageLBerger ( 2015-09-17 11:56:58 -0500 )edit

I think my code is good because i have found this post ( sorry it's in french)

LBerger gravatar imageLBerger ( 2015-09-18 06:00:05 -0500 )edit

super!!! thanks again ;-)

theodore gravatar imagetheodore ( 2015-09-18 08:10:57 -0500 )edit

1 answer

Sort by » oldest newest most voted

answered 2015-09-18 06:03:09 -0500

LBerger gravatar image

updated 2015-09-23 11:46:09 -0500

theodore gravatar image

it is my answer :

#include <opencv2/opencv.hpp> 
#include <opencv2/core/core.hpp>
#include <opencv2/imgproc/imgproc.hpp> 
#include "opencv2/highgui/highgui.hpp"
#include <iostream>
#include <fstream>
#include <string>
#include <sstream>
#include "opencv2/core/ocl.hpp"

using namespace cv; 
using namespace std;

int main (int argc,char **argv)

// Refernce (IEEE PAMI v12 n 7 1990)
Mat x = imread("f:/lib/opencv/samples/data/lena.jpg",CV_LOAD_IMAGE_GRAYSCALE);
Mat x0;
x.convertTo(x0, CV_32FC1);

double t=0;
double lambda=0.25; // Defined in equation (7)
double K=10; // defined after equation(13) in text

Mat x1,xc;

 while (t<20)
Mat D; // defined just before equation (5) in text
Mat gradxX,gradyX; // Image Gradient t time 
D = Mat::zeros(x0.size(),CV_32F);
for (int i=0;i<x0.rows;i++)
    for (int j = 0; j < x0.cols; j++)
        float gx =<float>(i, j), gy =<float>(i,j);
        float d;
        if (i==0 || i== x0.rows-1 || j==0 || j==x0.cols-1) // conduction coefficient set to 1 p633 after equation 13
            d =1.0/(1.0+abs((gx*gx + gy*gy))/(K*K)); // expression of g(gradient(I))
            //d =-exp(-(gx*gx + gy*gy)/(K*K)); // expression of g(gradient(I))<float>(i, j) = d;
x1 = Mat::zeros(x0.size(),CV_32F);
for (int i = 1; i < x0.rows-1; i++)
    float *u1 = (float*)x1.ptr(i);
    for (int j = 1; j < x0.cols-1; j++,u1++)
        // Value of I at (i+1,j),(i,j+1)...(i,j)
        float<float>(i+1, j),<float>(i, j+1);
        float<float>(i-1, j),<float>(i, j-1),<float>(i, j);
        // Value of D at at (i+1,j),(i,j+1)...(i,j)
        float<float>(i+1, j),<float>(i, j+1);
        float<float>(i-1, j),<float>(i, j-1),<float>(i, j);
        // Equation (7) p632
        *u1 = i00 + lambda/4*( (cp10+c00)*(ip10-i00) + (c0p1+c00)*(i0p1-i00) + (cm10+c00)*(im10-i00)+ (c0m1+c00)*(i0m1-i00));
        // equation (9)
imshow("Perrony x0",xc);
cout << "*";

return 0;
edit flag offensive delete link more


@LBerger, while the above code works perfectly, in my task I should use K = 0.5 (i.e. contrast value) over t = {1, 5, 10, 30, 100} and obtain a similar results as it gives now with K in higher values e.g. >10. So my guess is that the loop that simulates the diffusion time and the way that t is handled (i.e. t = t + lambda;) somehow needs some optimization, do you have an idea about it? Also I made some corrections on the code, since there were typos :-).

theodore gravatar imagetheodore ( 2015-09-23 12:00:28 -0500 )edit

I'm not sure that result are similar for t=30 and 100. I have changed lambda in source code. i cannot see convergence

        // Equation (7) p632
        double xx=(cp10+c00)*(ip10-i00) + (c0p1+c00)*(i0p1-i00) + (cm10+c00)*(im10-i00)+ (c0m1+c00)*(i0m1-i00);<float>(i, j) = xx;
        if (maxD<fabs(xx))
        // equation (9)
cout <<" lambda = "<< lambda<<"\t Maxd"<<maxD << "\t"<<intxx<<"\n";
for (int i = 1; i < x0.rows-1; i++)
    float *u1 = (float*)x1.ptr(i);
    for (int j = 1; j < x0.cols-1; j++,u1++)
        *u1 =<float>(i, j) + lambda/4*<float>(i, j);
        // equation (9)
LBerger gravatar imageLBerger ( 2015-09-24 06:29:58 -0500 )edit

It's better to process boundary

LBerger gravatar imageLBerger ( 2015-09-25 11:35:47 -0500 )edit

@LBerger thanks again, I'll have a look and I will tell you if I have any comment.

theodore gravatar imagetheodore ( 2015-09-26 04:26:40 -0500 )edit
LBerger gravatar imageLBerger ( 2015-10-01 00:34:01 -0500 )edit

@LBerger thanks really nice link ;-)

theodore gravatar imagetheodore ( 2015-10-01 03:54:48 -0500 )edit

Question Tools



Asked: 2015-09-16 09:41:58 -0500

Seen: 3,719 times

Last updated: Sep 23 '15