Ask Your Question
1

Get a 2D gaussian kernel for deconvolution

asked 2017-12-07 04:46:13 -0600

lock042 gravatar image

Hello,

I am trying to apply the Lucy-Richardson deconvolution algorithm with the following code:

static Mat RLTikh_deconvolution(Mat observed, Mat psf, int iterations) {

    Mat deconv = observed.clone();
    double mu = 0.01;

    // Iterate
    for (int i = 0; i < iterations; i++) {

            // Temporary matrix
            Mat ratio;
            filter2D(deconv, ratio, deconv.depth(), psf, Point(-1, -1), 0,
            BORDER_REFLECT);

            divide(observed, ratio, ratio);

            filter2D(ratio, ratio, ratio.depth(), psf, Point(-1, -1), 0,
            BORDER_REFLECT);

            // TV Regularization
            Mat denom;
            Laplacian(deconv, denom, deconv.depth(), 1, 1, 0, BORDER_REFLECT);
            denom = 1.0 - 2.0 * mu * denom;
            divide(ratio, denom, ratio);

            // Apply iteration on the estimate
            multiply(deconv, ratio, deconv);
    }

    return deconv;
}

I get the psf mat with Mat psf = getGaussianKernel(13, -1);

The problem is that this kernel is a 1D kernel but I would like to apply a 2D one. I believe I could apply the kernel two times (in X direction and Y direction) to get the result I want. But how I can do it ?

edit retag flag offensive close merge delete

1 answer

Sort by ยป oldest newest most voted
1

answered 2017-12-07 05:21:39 -0600

berak gravatar image

updated 2017-12-07 05:27:59 -0600

you can use sepFilter2D with your 1d kernel:

I believe I could apply the kernel two times (in X direction and Y direction)

that's exactly, what it does.

static Mat RLTikh_deconvolution(Mat observed, Mat psf, int iterations) {

    Mat deconv = observed.clone();
    double mu = 0.01;

    // Iterate
    for (int i = 0; i < iterations; i++) {

            // Temporary matrix
            Mat ratio;
            sepFilter2D(deconv, ratio, deconv.depth(), psf, psf, Point(-1, -1), 0,
            BORDER_REFLECT);

            divide(observed, ratio, ratio);

            sepFilter2D(ratio, ratio, ratio.depth(), psf, psf, Point(-1, -1), 0,
            BORDER_REFLECT);

            // TV Regularization
            Mat denom;
            Laplacian(deconv, denom, deconv.depth(), 1, 1, 0, BORDER_REFLECT);
            denom = 1.0 - 2.0 * mu * denom;
            divide(ratio, denom, ratio);

            // Apply iteration on the estimate
            multiply(deconv, ratio, deconv);
    }

    return deconv;
}
edit flag offensive delete link more

Comments

1

So easy. Sorry for this dummy question :). And of course, thank you for your answer.

lock042 gravatar imagelock042 ( 2017-12-07 05:30:43 -0600 )edit
1

well, you have to know , that it exists, then it's easy ...

berak gravatar imageberak ( 2017-12-07 05:39:08 -0600 )edit

Question Tools

1 follower

Stats

Asked: 2017-12-07 04:46:13 -0600

Seen: 2,126 times

Last updated: Dec 07 '17