1
$\begingroup$

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:

  1. Convolution of a gaussian profile with a rectangular function results in a gaussian profile.
  2. The PSF and the observed image have the same size (440 x 440).
  3. 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:

enter image description here

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:

Division in fourier domain

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:

Van Cittert Implementation

Richardson-Lucy implementation

What we do not understand is:

  1. In the wiener implementation, where does the noise come from, provided that both the input images did not have any noise?
  2. 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)?
  3. Is there another approach to this problem?
  4. Is something amiss in my code?

Any help is greatly appreciated!

$\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.