# isotropic non-linear diffusion smoothing (Perona-Malik)

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:

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

(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 close merge delete

2

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

I proposed this program

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

( 2015-09-17 10:36:00 -0500 )edit
1

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

( 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)

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

super!!! thanks again ;-)

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

Sort by » oldest newest most voted

#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 http://image.diku.dk/imagecanon/material/PeronaMalik1990.pdf (IEEE PAMI v12 n 7 1990)
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
imshow("Original",x);

Mat x1,xc;

while (t<20)
{
Mat D; // defined just before equation (5) in text
D = Mat::zeros(x0.size(),CV_32F);
for (int i=0;i<x0.rows;i++)
for (int j = 0; j < x0.cols; 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;
else
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))
D.at<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);
u1++;
for (int j = 1; j < x0.cols-1; j++,u1++)
{
// Value of I at (i+1,j),(i,j+1)...(i,j)
float ip10=x0.at<float>(i+1, j),i0p1=x0.at<float>(i, j+1);
float im10=x0.at<float>(i-1, j),i0m1=x0.at<float>(i, j-1),i00=x0.at<float>(i, j);
// Value of D at at (i+1,j),(i,j+1)...(i,j)
float cp10=D.at<float>(i+1, j),c0p1=D.at<float>(i, j+1);
float cm10=D.at<float>(i-1, j),c0m1=D.at<float>(i, j-1),c00=D.at<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)
}
}
x1.copyTo(x0);
x0.convertTo(xc,CV_8U);
imshow("Perrony x0",xc);
cout << "*";
waitKey(10);
t=t+lambda;
}

waitKey();
return 0;
}

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 :-).

( 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);
dI00.at<float>(i, j) = xx;
if (maxD<fabs(xx))
maxD=fabs(xx);
intxx+=fabs(xx);
// equation (9)
}
}
lambda=100/maxD;
cout <<" lambda = "<< lambda<<"\t Maxd"<<maxD << "\t"<<intxx<<"\n";
for (int i = 1; i < x0.rows-1; i++)
{
float *u1 = (float*)x1.ptr(i);
u1++;
for (int j = 1; j < x0.cols-1; j++,u1++)
{
*u1 = x0.at<float>(i, j) + lambda/4*dI00.at<float>(i, j);
// equation (9)
}

( 2015-09-24 06:29:58 -0500 )edit

It's better to process boundary

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

( 2015-09-26 04:26:40 -0500 )edit
( 2015-10-01 00:34:01 -0500 )edit

@LBerger thanks really nice link ;-)

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

Official site

GitHub

Wiki

Documentation