# 3DFilter for PSF image Richard-Lucy Deconvolution

Hello everyone,

the code posted below shows an implementation of a Richard-Lucy deconvolution, which is supposed to recover a blurred image using the point spread function. I found this on GitHub (RL_deconv/rl_deconv.cpp at master · chrrrisw/RL_deconv · GitHub) and it was a good starting point for my problem. The input given is a 3D image and a 2D PSF.

As I said before, this program only works with 2D psf images, but I would to do this with a 3D psf. My research for a 3D filter was unfortunately without result.

My question now is whether there is a 3D filter in OpenCV or is there another approach for this task?

I would appreciate any help.

``````#define _USE_MATH_DEFINES

#include "opencv2/opencv.hpp"
#include <iostream>
#include <cmath>

using namespace cv;
using namespace std;

// From wikipedia:
//
// def RL_deconvolution(observed, psf, iterations):
//     # initial estimate is arbitrary - uniform 50% grey works fine
//     latent_est = 0.5*np.ones(observed.shape)
//     # create an inverse psf
//     psf_hat = psf[::-1,::-1]
//     # iterate towards ML estimate for the latent image
//     for i in np.arange(iterations):
//         est_conv      = cv2.filter2D(latent_est,-1,psf)
//         relative_blur = observed/est_conv;
//         error_est     = cv2.filter2D(relative_blur,-1,psf_hat)
//         latent_est    = latent_est * error_est
//     return latent_est

static int image_type;

Mat RL_deconvolution(Mat observed, Mat psf, int iterations) {

Scalar grey;

// Uniform grey starting estimation
switch (image_type) {
case CV_64FC1:
grey = Scalar(0.5);
case CV_64FC3:
grey = Scalar(0.5, 0.5, 0.5);
}
Mat latent_est = Mat(observed.size(), image_type, grey);

// Flip the point spread function (NOT the inverse)
Mat psf_hat = Mat(psf.size(), CV_64FC1);
int psf_row_max = psf.rows - 1;
int psf_col_max = psf.cols - 1;
for (int row = 0; row <= psf_row_max; row++) {
for (int col = 0; col <= psf_col_max; col++) {
psf_hat.at<double>(psf_row_max - row, psf_col_max - col) =
psf.at<double>(row, col);
}
}

Mat est_conv;
Mat relative_blur;
Mat error_est;

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

filter2D(latent_est, est_conv, -1, psf);

// Element-wise division
relative_blur = observed.mul(1.0 / est_conv);

filter2D(relative_blur, error_est, -1, psf_hat);

// Element-wise multiplication
latent_est = latent_est.mul(error_est);
}

return latent_est;
}

int main(/*int argc, const char** argv*/)
{

//if (argc != 3) {
//	cout << "Usage: " << argv[0] << " image iterations" << "\n";
//	return -1;
//}

//int iterations = atoi(argv[2]);
int iterations = 10;

Mat original_image;

int num_channels = original_image.channels();
switch (num_channels) {
case 1:
image_type = CV_64FC1;
break;
case 3:
image_type = CV_64FC3;
break;
default:
return -2;
}

// This is a hack, assumes too much
int divisor;
switch (original_image.elemSize() / num_channels) {
case 1:
divisor = 255;
break;
case 2:
divisor = 65535;
break;
default:
return -3;
}

// From here on, use 64-bit floats
// Convert original_image to float
Mat float_image;
original_image.convertTo(float_image, image_type);
float_image *= 1.0 / divisor;
namedWindow("Float", WINDOW_AUTOSIZE); // CV_WINDOW_AUTOSIZE replaced with WINDOW_AUTOSIZE
imshow("Float", float_image);

// Calculate a gaussian blur psf.
double sigma_row = 9.0;
double sigma_col = 5.0;
int psf_size = 5;
double mean_row = 0.0;
double mean_col = psf_size / 2.0;
double sum = 0.0;
double temp;
Mat psf = Mat(Size(psf_size, psf_size), CV_64FC1, 0.0);

for (int j = 0; j < psf.rows; j++) {
for (int k = 0; k < psf.cols; k++) {
temp = exp(
-0.5 * (
pow((j - mean_row) / sigma_row, 2.0) +
pow((k - mean_col) / sigma_col, 2.0))) /
(2 * M_PI * sigma_row * sigma_col);
sum += temp;
psf.at<double>(j, k) = temp;
}
}

// Normalise the psf.
for (int row = 0; row < psf.rows; row++) {
// cout << row << " ";
for (int col = 0; col < psf.cols; col++) {
psf.at<double>(row, col) /= sum;
// cout << psf.at<double>(row, col) << " ";
}
// cout << "\n";
}

//Mat psf;

// Blur the float_image with the psf.
Mat blurred_float;
blurred_float = float_image.clone();
filter2D(float_image, blurred_float, -1, psf);
namedWindow("BlurredFloat", WINDOW_AUTOSIZE); // CV_WINDOW_AUTOSIZE replaced with WINDOW_AUTOSIZE
imshow("BlurredFloat", blurred_float);

Mat estimation = RL_deconvolution(blurred_float, psf, iterations);
namedWindow("Estimation", WINDOW_AUTOSIZE); // CV_WINDOW_AUTOSIZE replaced with WINDOW_AUTOSIZE
imshow("Estimation", estimation);

waitKey(0); //wait infinite time for a keypress

destroyWindow("Float");
destroyWindow("BlurredFloat");
destroyWindow("Estimation");

return 0;
}
``````

does that mean you just need to “lift” that code from 2D to 3D?

so you’ll need to replace `filter2D()` with something like a “filter3D”, which opencv doesn’t have, but you could build. you could use fourier transforms to make that fast. IIRC, filter2D uses those too, unless the direct calculation is faster (small kernels)

also you’ll need to lift all the 2D data to 3D. OpenCV’s Mat class can do that (arbitrary dimensions) but it’s not terribly convenient. I’d recommend working in python. then you have numpy, scipy, scikit-image, … and OpenCV can chip in too, if the array shape isn’t too wild.

all the “heavy” calculations appear to happen in the iteration that uses `filter2D`. everything else looks like constant effort, so it’s OK to be written in plain python (even if not, there’s `numba`)

does that mean you just need to “lift” that code from 2D to 3D?

Yes, I would like to “lift” the 2D psf into a 3D psf or rather load a 3D psf via OpenCV.

I’d recommend working in python.

Python is of course better, but have to implement this in C++ and I think OpenCV is for that case the best package.

so you’ll need to replace `filter2D()` with something like a “filter3D”, which opencv doesn’t have, but you could build. you could use fourier transforms to make that fast. IIRC, filter2D uses those too, unless the direct calculation is faster (small kernels)

Yeah the problem is, that there isn’ t a 3D filter in OpenCV. The question is: Is there already an approach of a 3D filter or do I simply have to orientate myself on the 2D filter and build a 3D filter based on it?