I am currently working on a project, where we record an electron beam profile using a target. The obtained image is a result of convolution of the actual beam profile and the aperture wherein the aperture is respresented by an image of the same size as the beam profile with all pixels lying inside the circle (aperture) having high intensity values( 60000 for 16bit images) and those lying outside the circle having 0 intensity values.
I have been trying for the past two weeks to implement an appropriate deconvolution to obtain the actual beam profile from the observed beam profile signal.
Our Assumptions:
- Convolution of a gaussian profile with a rectangular function results in a gaussian profile.
- The PSF and the observed image have the same size (440 x 440).
- The observed image has a slightly larger diameter than the expected original image (due to convolution).
I have tried using Wiener Deconvolution implementation mentioned here: https://stackoverflow.com/questions/40713929/weiner-deconvolution-using-opencv
with the following results on simulated images where it seems to work fine:
However on an actual beam profile the out put is a rather noisy image.
Our expectations were in fact to retrieve a gauss profile with a smaller diameter. All the images (except the last two) used were simulated images which do not have any noise or background values.
The code used for Wiener Deconvolution is as follows:
bool wienerDeconvolution(const cv::Mat *src1, const cv::Mat *src2, cv::Mat *dst)
{
// create necessary "real" input images
cv::Mat realImage1(src1->size(), CV_64F);
src1->convertTo(realImage1, CV_64F);
cv::Mat realImage2(src2->size(), CV_64F);
src2->convertTo(realImage2, CV_64F);
// DFT of the source image
cv::Mat planesSrc[] = {cv::Mat_<double>(realImage1), cv::Mat::zeros(realImage1.size(), CV_64F)};
cv::Mat complexSrc;
merge(planesSrc, 2, complexSrc); // Add to the expanded another plane with zeros
dft(complexSrc, complexSrc, cv::DFT_COMPLEX_OUTPUT);
cv::split(complexSrc, planesSrc);
cv::Mat planesSrcOrg[] = {cv::Mat_<double>(realImage1), cv::Mat::zeros(realImage1.size(), CV_64F)};
planesSrcOrg[0] = planesSrc[0](cv::Rect(0, 0, realImage1.rows, realImage1.cols));
planesSrcOrg[1] = planesSrc[1](cv::Rect(0, 0, realImage1.rows, realImage1.cols));
// DFT of the PSF
cv::Mat planesPsf[] = {cv::Mat_<double>(realImage2), cv::Mat::zeros(realImage2.size(), CV_64F)};
cv::Mat complexPsf;
merge(planesPsf, 2, complexPsf); // Add to the expanded another plane with zeros
dft(complexPsf, complexPsf, cv::DFT_COMPLEX_OUTPUT);
cv::split(complexPsf, planesPsf);
cv::Mat planesPsfOrg[] = {cv::Mat_<double>(realImage2), cv::Mat::zeros(realImage2.size(), CV_64F)};
planesPsfOrg[0] = planesPsf[0](cv::Rect(0, 0, realImage2.rows, realImage2.cols));
planesPsfOrg[1] = planesPsf[1](cv::Rect(0, 0, realImage2.rows, realImage2.cols));
//--------------------------------------------------
cv::Mat Fu[2];
Fu[0] = cv::Mat::zeros(realImage1.rows, realImage1.cols, CV_64FC1);
Fu[1] = cv::Mat::zeros(realImage1.rows, realImage1.cols, CV_64FC1);
std::complex<double> a;
std::complex<double> b;
std::complex<double> c;
double Hf_Re;
double Hf_Im;
double Phf;
double hfz;
double hz;
double A;
const double eps = 1E-8;
for (int i = 0; i < realImage1.rows; i++)
{
for (int j = 0; j < realImage1.cols; j++)
{
Hf_Re = planesPsf[0].at<double>(i, j);
Hf_Im = planesPsf[1].at<double>(i, j);
Phf = (Hf_Re * Hf_Re + Hf_Im * Hf_Im);
hfz = (Phf < eps) * eps;
hz = (realImage2.at<double>(i, j) > 0);
A = Phf / (Phf + hz + 0.01);
a = std::complex<double>(planesSrc[0].at<double>(i, j), planesSrc[1].at<double>(i, j));
b = std::complex<double>(Hf_Re + hfz, Hf_Im + hfz);
c = a / b; // Deconvolution :) other work to avoid division by zero
Fu[0].at<double>(i, j) = (c.real() * A);
Fu[1].at<double>(i, j) = (c.imag() * A);
}
}
cv::Mat outputFft;
cv::merge(Fu, 2, outputFft);
dft(outputFft, outputFft, cv::DFT_INVERSE + cv::DFT_SCALE);
split(outputFft, Fu);
cv::Mat dstFinal = Fu[0];
reOrderMagnitudeImage(dstFinal);
const double rangeMax {1.0 * 60000/65535.0};
cv::normalize(dstFinal, dstFinal, 0, rangeMax, CV_MINMAX);
dstFinal.convertTo(*dst, CV_16U, 65535);
return true;
}
Next I tried the solution mentioned in link: 2D deconvolution of recorded electron beam data wherein as per my understanding they perform an element wise division in the fourier domain (since deconvolution is nothing but division in fourier domain). The code is as below:
bool deconvolution2(const cv::Mat *src1, const cv::Mat *src2, cv::Mat *dst)
{
// Convert the input images into CV_64F forma
cv::Mat srcImage(src1->size(), CV_64F);
src1->convertTo(srcImage, CV_64F);
cv::Mat psfImage(src2->size(), CV_64F);
src1->convertTo(psfImage, CV_64F);
// Perform forward fourier transform
cv::Mat planesSrc[] = {cv::Mat_<double>(srcImage), cv::Mat::zeros(srcImage.rows, srcImage.cols, CV_64F)};
cv::Mat planesPsf[] = {cv::Mat_<double>(psfImage), cv::Mat::zeros(psfImage.rows, psfImage.cols, CV_64F)};
cv::Mat complexSrc, complexPsf;
cv::merge(planesSrc, 2, complexSrc);
cv::merge(planesPsf, 2, complexPsf);
cv::dft(complexSrc, complexSrc, cv::DFT_COMPLEX_OUTPUT);
cv::dft(complexPsf, complexPsf, cv::DFT_COMPLEX_OUTPUT);
// Separate real and imaginary parts
cv::split(complexSrc, planesSrc);
cv::split(complexPsf, planesPsf);
double epsilon = 1E-8;
// Perform elementwise division
cv::Mat dest[] = {cv::Mat::zeros(complexSrc.rows, complexSrc.cols, CV_64F),
cv::Mat::zeros(complexSrc.rows, complexSrc.cols, CV_64F)};
for (int row = 0; row < complexSrc.rows; ++row)
{
for (int col = 0; col < complexSrc.cols; ++col)
{
double denom = planesPsf[0].at<double>(row, col) * planesPsf[0].at<double>(row, col) +
planesPsf[1].at<double>(row, col) * planesPsf[1].at<double>(row, col) + epsilon;
double real = planesSrc[0].at<double>(row, col) * planesPsf[0].at<double>(row, col) +
planesSrc[1].at<double>(row, col) * planesPsf[1].at<double>(row, col);
double imaginary = planesSrc[1].at<double>(row, col) * planesPsf[0].at<double>(row, col) -
planesSrc[0].at<double>(row, col) * planesPsf[1].at<double>(row, col);
dest[0].at<double>(row, col) = real / denom;
dest[1].at<double>(row, col) = imaginary / denom;
}
}
cv::Mat destOut;
cv::merge(dest, 2, destOut);
// Perform idft
cv::dft(destOut, destOut, cv::DFT_INVERSE + cv::DFT_SCALE);
cv::split(destOut, dest);
cv::Mat destReordered = dest[0];
reOrderMagnitudeImage(destReordered);
cv::normalize(destReordered, destReordered, 0, 1 , CV_MINMAX);
destReordered.convertTo(*dst, CV_16UC1, 65535);
return true;
}
with the following output:
Next I implemented Van-Cittert Deconvolution: https://arxiv.org/pdf/1707.09177.pdf and https://de.wikipedia.org/wiki/Van-Cittert-Dekonvolution and Richardson-Lucy deconvolution: https://github.com/chrrrisw/RL_deconv/blob/master/rl_deconv.cpp code as follows:
bool rlDeconvolution(const cv::Mat *src1, const cv::Mat *src2, cv::Mat *dst, int numIteration)
{
// create necessary "real" input images
cv::Mat observedImage(src1->size(), CV_64F);
src1->convertTo(observedImage, CV_64F);
cv::Mat psfImage(src2->size(), CV_64F);
src2->convertTo(psfImage, CV_64F);
cv::Mat latent_est = cv::Mat(src1->size(), CV_64F, cv::Scalar(0.5));
// Flip the point spread function (NOT the inverse)
cv::Mat psf_hat = cv::Mat(psfImage.size(), CV_64FC1);
int psf_row_max = psfImage.rows - 1;
int psf_col_max = psfImage.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) =
psfImage.at<double>(row, col);
}
}
cv::Mat est_conv;
cv::Mat relative_blur;
cv::Mat error_est;
// Iterate
for (int i=0; i<numIteration; i++) {
filter2D(latent_est, est_conv, -1, psfImage);
// Element-wise division
relative_blur = observedImage.mul(1.0/est_conv);
filter2D(relative_blur, error_est, -1, psf_hat);
// Element-wise multiplication
latent_est = latent_est.mul(error_est);
}
latent_est.convertTo(*dst, CV_16UC1, 65535);
return true;
}
bool FilterFunctions::vanCittertDeconvolution(const cv::Mat *src1, const cv::Mat *src2, cv::Mat *dst, int numIteration)
{
// create necessary "real" input images
cv::Mat observedImage(src1->size(), CV_32S);
src1->convertTo(observedImage, CV_32S);
cv::Mat psfImage(src2->size(), CV_32S);
src2->convertTo(psfImage, CV_32S);
// convolve the observed image with the gauss kernel -> apperture represented by src2 in our case
cv::Mat latentEst = observedImage.clone();
cv::Mat estConv, diffImage, sumImage;
for (int count = 0; count < numIteration; ++count)
{
// cv::filter2D(latentEst, estConv, -1, psfImage);
convolution(&latentEst, &psfImage, &estConv);
// TypeConversion
cv::Mat estConv32S;
estConv.convertTo(estConv32S, CV_32S);
// create a difference image
cv::subtract(observedImage, estConv32S, diffImage);
// sum of the diff image and the original image
cv::add(diffImage, observedImage, latentEst);
}
// normalize the estimated image
cv::normalize(latentEst, latentEst, 0, 60000./65535., CV_MINMAX);
latentEst.convertTo(*dst, CV_16UC1, 65535);
return true;
}
with following outputs:
What we do not understand is:
- In the wiener implementation, where does the noise come from, provided that both the input images did not have any noise?
- Are these methods the correct approach for deconvolution between a rectangular and a gaussian function (in almost all links on deconvolution, they speak about a gaussian PSF while our PSF is a rectangular function)?
- Is there another approach to this problem?
- Is something amiss in my code?
Any help is greatly appreciated!
