I've designed a filter based on this equation

$$ E(u) = \frac{1}{2}\int_{\Omega} \left[(u-f)^2 + g(\lVert \nabla u \rVert)\right]dxdy $$

Euler-Lagrange equations in this case are (should be) given by

$$ u - f - \frac{1}{2}\nabla^T \cdot \left( \frac{\nabla u}{\lVert \nabla u\rVert}g'(\lVert \nabla u\rVert)\right) = 0 $$

For given $\lambda > 0$ if I chose $g(x) = e^{-x/\lambda}$ the EL equations are $$u - f + \frac{1}{2\lambda}\nabla^T \cdot \left( \frac{\nabla u}{\lVert \nabla u\rVert}e^{-\lVert \nabla u\rVert/\lambda}\right) = 0$$

Continuous gradient descent in this case is given by

$$ u_t = f - u - \frac{1}{2\lambda}\nabla^T \cdot \left( \frac{\nabla u}{\lVert \nabla u\rVert}e^{-\lVert \nabla u\rVert/\lambda}\right) $$

Basically my idea would be... if the gradient magnitude is very high then $u = f$ (i.e. I'd like to preserve strong edges) otherwise (small variations) $u$ should be tending to some constant, or anyway smaller variations.

Now... This is my implementation in opencv

void g(InputArray input, OutputArray output) {
    float lambda = 0.0000001;
    Mat rval, arg = input.getMat();
    //cv::exp(-arg / lambda2, tmp);

void dg(InputArray input, OutputArray output) {
    Mat rval, tmp, arg = input.getMat();
    float lambda = 0.0000001;
    rval = (-1.0/lambda)*tmp;

void evolveImage(const cv::Mat & input, cv::Mat & output, void (*g)(cv::InputArray, cv::OutputArray), void(*dg)(cv::InputArray, cv::OutputArray), const int & k_max) {
    int X = input.cols;
    int Y = input.rows;

    Mat imfloat;
    input.convertTo(imfloat, CV_32FC1);
    Mat u;
    cv::Mat & f = imfloat;

    double minval, maxval;

    minMaxLoc(u, &minval, &maxval);
    Mat gradientMagnitude;

    Mat dx, dy;
    Mat gprime, a, b;
    for (auto k = 0; k < k_max; ++k) {

        Sobel(u, dx, CV_32F, 1, 0, 3); //Estimation differences in x
        Sobel(u, dy, CV_32F, 0, 1, 3); //Estimation differences in y
        cv::sqrt(dx.mul(dx) + dy.mul(dy),gradientMagnitude);
        dg(gradientMagnitude, gprime);
        minMaxLoc(dx, &minval, &maxval);
        minMaxLoc(dy, &minval, &maxval);
        a = gprime.mul(dx);
        b = gprime.mul(dy);
        Sobel(a, dx, CV_32F, 1, 0, 3);
        Sobel(b, dy, CV_32F, 0, 1, 3);

        u = f + 0.5f*(dx + dy);


    minMaxLoc(u, &minval, &maxval);
    u.convertTo(output, CV_8UC1, 255 / (maxval - minval), -minval * 255 / (maxval - minval));

I have this messy function that can be called by a main if you do cut and paste

int continuousGradientDescentExample(int argc, char **argv) {
    if (argc != 2) {
        std::cout << "Need filename..." << std::endl;
        return 1;

    Mat im = imread(argv[1], CV_LOAD_IMAGE_GRAYSCALE);

    namedWindow("Test", WINDOW_AUTOSIZE);
    imshow("Test", im);

    Mat imfloat;
    im.convertTo(imfloat, CV_32FC1);
    Mat u(imfloat.size(), imfloat.type());

    int & X = im.cols;
    int & Y = im.rows;

    evolveImage(imfloat, u, g, dg, 50);

    std::cout << type2str(u.type()) << std::endl;

    Mat output(Size(2 * X, Y), CV_8UC1);
    im.copyTo(output(Rect(0, 0, X, Y)));
    u.convertTo(u, CV_8UC1);
    u.copyTo(output(Rect(X, 0, X, Y)));

    namedWindow("Original vs Anisotropic diffusion", WINDOW_AUTOSIZE);
    imshow("Original vs Anisotropic diffusion", output);


    return 1;

The code is working (i.e. doesn't crash, you can try it if you want). But the thing is no matter how I choose $\lambda$ in the evolveImage the matrices $a,b$ are always $0$ which implies the update of $u$ doesn't produce any effect. To me the equation makes sense, but I can't understand why there's no difference at all.

Can anyone give me an insight?


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.