1
$\begingroup$

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) {
    //cv::exp(-im,out);
    float lambda = 0.0000001;
    Mat rval, arg = input.getMat();
    //cv::exp(-arg / lambda2, tmp);
    cv::exp(-arg/lambda,rval);
    output.assign(rval);
}

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

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;
    imfloat.copyTo(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);
    waitKey(0);

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

    waitKey(0);

    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?

$\endgroup$

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.