decon() sometimes changes the size of image.

[ajan890]

Hi @tlambert03,

While testing my code, I noticed that sometimes the image (in this case, ndarrays loaded from .nrrd files) would change shape seemingly randomly. I am unsure about the cause of this and believe it is some sort of bug present in the library. I would appreciate it if you are able to give some pointers on how I can work around this issue.

Additionally, I noticed that the edges of the image are not being deconvolved. Correct me if I am wrong, but I think this is due to the nature of the algorithm, hence the padding before applying the deconvolution.

I have provided a snippet of my code, as well as some sample data in case you want to reproduce the error. You will have to extract the zip folder and provide the file paths as displayed in the main function.

Thank you!

sample_data.zip

import nrrd
from pathlib import Path
from numpy import pad
from pycudadecon import decon, make_otf

def process_cube(
        input_file: Path,   # path to nrrd file
        temp_folder: Path,  # path to a temporary folder (folder must exist before function runs)
        psf: str,           # string representing path to psf file
        pad_amount: int,    # amount to pad to each side of image before applying deconvolution
        n_iters: int = 9,
        dz_data: float = 0.7,
        dx_data: float = 0.7,
        dz_psf: float = 1.8,
        dxy_psf: float = 0.2,
        background: int = 0,
        wavelength_em: int = 525,
        na: float = 0.4,
        nimm: float = 1.42
): 
    # read input file
    img, header = nrrd.read(input_file.__str__())

    # generate otf file from psf file
    otf_file = make_otf(
            psf,
            outpath=(temp_folder / 'otf.tif').__str__(),
            dzpsf=dz_psf,
            dxpsf=dxy_psf,
            wavelength=wavelength_em,
            na=na,
            nimm=nimm,
            fixorigin=0,
        )

    print("Image shape before pad: ", img.shape)

    # pad to ensure edges of image are deconvolved
    img = pad(img, pad_amount, 'reflect')
    
    print("Image shape (after pad), before deconvolution: ", img.shape)

    # run deconvolution    
    img = decon(
        img,
        otf_file,
        fpattern=None,
        n_iters=n_iters,
        dzdata=dz_data,
        dxdata=dx_data,
        dzpsf=dz_psf,
        dxpsf=dxy_psf,
        background=background,
        wavelength=wavelength_em,
        na=na,
        nimm=nimm
    )
    print("Image shape after deconvolution: ", img.shape)
    
    # ... rest of code is irrelevant

if __name__ == '__main__':
    # change pad_amount to different integers.  sometimes the deconvolution causes image to change size.
    # .nrrd file along with psf.tif are provided as sample data.

    process_cube(
        Path("C:/path/to/file/y00006400.x00019200.nrrd"),
        Path("C:/path/to/temp/folder"),
        "C:/path/to/psf.tif",
        pad_amount=1
    )

[tlambert03]

Hi @ajan890 and thanks for your note, and the example code. Could you save me a bit of time and also paste the output that you're getting when you run that script? I'm specifically curious to see exactly what the shapes are.

I suspect the reason is that cudadecon pads the image to an efficient size for FFTs (which is one that can be factorized into 2s, 3s, 5s, and 7s as per CUFFT recs). It does that here:

https://github.com/scopetools/cudadecon/blob/a679bf4c2a7e2c7579ca47a7719d0e99ea72974a/src/RL-Biggs-Andrews.cpp#L18-L43

If it's really important for your application that the output be exactly the same size as the input, I think i would probably create one additional small function that simply pads or crops the result to the size of your input image (and if you wanted to contribute that here, i would be happy to discuss a PR)

Additionally, I noticed that the edges of the image are not being deconvolved. Correct me if I am wrong, but I think this is due to the nature of the algorithm, hence the padding before applying the deconvolution.

Indeed, deconvolution will always have some edge artifacts, and your strategy of padding a bit with reflection (and preferably chopping that off afterwards) is about as good of a strategy as any to try to get everything you can out of the edges.

[he113-code]

Hi, I have noticed that my images are also cropped, and as I am slicing the image for prediction, I would like to return images of the same size.Can you tell me what action I should take to avoid this?
code:

import numpy as np
from pycudadecon import decon

psf_array = np.ones((100,100,100))
image = np.zeros(((440, 512, 512)))
print(1,image.shape)
result22 = decon(image, psf_array,n_iters = 5)
print(2,result22.shape)

return：
C:\Users\Peng\AppData\Local\Temp\tmpbepdgztj.tif
callednx=100, ny=100, nz=100
no peak: a1=5022346917897180067849241100288.000000, a2=5022346917897180067849241100288.000000, a3=5022346917897180067849241100288.000000, slope=0.000000, curvature=0.000000
Center of mass is (74.000, 60.000, 0.600)
Background is -0.037
In radialft()
output nz=432
X_k allocated. 432MB 21541MB free Pinning raw.data's Host RAM. Copy raw.data to X_k HostToDevice. Done.
CC allocated. 432MB 21109MB free
rawGPUbuf allocated. 432MB 20677MB free
X_kminus1 allocated. 432MB 20245MB free
Y_k allocated. 432MB 19813MB free
G_kminus1 allocated. 432MB 19381MB free
G_kminus2 allocated. 432MB 18949MB free
fftGPUbuf allocated. 433MB 18515MB free

The final image obtained, z = 432

What rules should zxy conform to in order to avoid cropping? @tlambert03

[tlambert03]

What rules should zxy conform to in order to avoid cropping? @tlambert03

Hi @he113-code, please see my previous comment. It is not a simple explanation. I linked to the logic that determines whether cropping occurs or not. and also made a suggestion:

If it's really important for your application that the output be exactly the same size as the input, I think i would probably create one additional small function that simply pads or crops the result to the size of your input image

[he113-code]

Thanks for the tip, it took me some time to understand the C++ function, for the follow up tip I am providing the python code below:

def not_good_dimension(num):
    if num in (2, 3, 5, 7):
        return False
    elif num % 2 == 0:
        return not_good_dimension(num // 2)
    elif num % 3 == 0:
        return not_good_dimension(num // 3)
    elif num % 5 == 0:
        return not_good_dimension(num // 5)
    elif num % 7 == 0:
        return not_good_dimension(num // 7)
    else:
        return True

def find_optimal_dimension(in_size, step=-1):
    out_size = in_size
    while not_good_dimension(out_size):
        out_size += step
    return out_size


initial_size = 440
optimal_size = find_optimal_dimension(initial_size)
print(f"Optimal dimension for FFT: {optimal_size}")