astra-toolbox / astra-toolbox

ASTRA Tomography Toolbox
GNU General Public License v3.0
455 stars 208 forks source link

Cuda memory error #347

Closed vlibX closed 1 year ago

vlibX commented 2 years ago

Hello ! I am using a script to reconstruct an unusual geometry (vertical fan beam on a linear detector, translated step by step). I have only 60 projections of 1792980 pixels. The reconstructed volume is 400250*250 voxels. I use a nvidia rtx 3090 with 24Go RAM.

I was nicely provided a reconstruction script that worked very well when the sinogram was downsampled from a factor 4. Given I have much than enough memory, I decided to try and reconstruct with the full sinogram. I then obtained this error message:

Error: CUDA error 11: invalid argument. Error: Failed to allocate 1792x58800x1 GPU array Error: CUDA error 11: invalid argument. Assertion failed: err == cudaSuccess, file E:/wjp/git/astra-toolbox/cuda/3d/util3d.cu, line 383

Could anyone please help me with this issue ?

V.

askorikov commented 1 year ago

Hi! Could you maybe post the full script (well the part related to reconstruction) if you couldn't resolve this problem in the end?

vlibX commented 1 year ago

Hi !

Sure thing, here is the script below.

Thank you for having a look !

def reconstruction(dirname, listImages, i0, sid, sdd, spacing, sf):     sinogram = None     n_proj = len(listImages)     mult = n_proj/sflistImages[0].shape[1]     n_lines = mult//n_proj     newWidth = n_linessf     image_dims = (int(listImages[0].shape[0]), int(newWidth))     for imarray in listImages:         imarray = imarray.astype(np.float32)         if sinogram is None:             if imarray.shape != image_dims:                 imarray = imarray[:, (imarray.shape[1]-image_dims[1]) // 2:image_dims[1]+(imarray.shape[1]-image_dims[1])//2]             sinogram = imarray         else:             if imarray.shape != image_dims:                 imarray = imarray[:, (imarray.shape[1]-image_dims[1]) // 2:image_dims[1]+(imarray.shape[1]-image_dims[1])//2]             sinogram = np.concatenate((sinogram, imarray), axis=1)

    sinogram[sinogram > i0] = i0     #n_proj = sinogram.shape[1]//image_dims[1]     sinogram = np.log(i0/sinogram)     sinogram = sinogram[:, :, np.newaxis]     sinogram = sinogram[:, ::sf]

    n_lines = (sinogram.shape[1])//(n_proj)     n_pixels = sinogram.shape[0]     SOD = sid  # 1895 #mm     SDD = sdd  # 2324 #mm     len_pixel = spacing  # 0.8 #mm     len_lineshift = spacing*sf  # mm - takes into account the downsample, equals to 4

    vol_geom = astra.create_vol_geom(         400, 250, 250, -400, 400, 0, 1280, -400, 400)

    angles = np.linspace(0, 2np.pi, n_proj, False)     vectors = np.zeros((n_projn_lines, 12))

    for angle_idx in range(n_proj):         for line_idx in range(n_lines):             i = line_idx + angle_idx*n_lines             a = -angles[angle_idx]

            s = np.array([-SOD, (line_idx-n_lines/2)*len_lineshift])             # source             vectors[i, 0] = np.dot(np.array([np.cos(a), -np.sin(a)]), s)             vectors[i, 1] = 0  # fixed             vectors[i, 2] = np.dot(np.array([np.sin(a), np.cos(a)]), s)

            d = np.array([SDD-SOD, (line_idx-n_lines/2)len_lineshift])             # center of detector             vectors[i, 3] = np.dot(np.array([np.cos(a), -np.sin(a)]), d)             vectors[i, 4] = len_pixeln_pixels/2  # fixed             vectors[i, 5] = np.dot(np.array([np.sin(a), np.cos(a)]), d)

            # vector from detector pixel (0,0) to (0,1)             vectors[i, 6] = 0  # fixed             vectors[i, 7] = len_pixel  # fixed             vectors[i, 8] = 0  # fixed

            # vector from detector pixel (0,0) to (1,0)             vectors[i, 9] = -np.sin(a)len_lineshift             vectors[i, 10] = 0  # fixed             vectors[i, 11] = np.cos(a)len_lineshift

    # invertire n_pixels e 1 se le cose vanno male     proj_geom = astra.create_proj_geom('cone_vec', 1, n_pixels, vectors)     sinogram = np.transpose(sinogram, (2, 1, 0))     proj_id = astra.data3d.create('-proj3d', proj_geom, sinogram)     rec_id = astra.data3d.create('-vol', vol_geom)

    cfg = astra.astra_dict('SIRT3D_CUDA')     cfg['ReconstructionDataId'] = rec_id     cfg['ProjectionDataId'] = proj_id     cfg['option'] = {}

    astra.astra.set_gpu_index(0, 0)     alg_id = astra.algorithm.create(cfg)     astra.algorithm.run(alg_id, allParams["nbIter"])

    rec = copy.deepcopy(astra.data3d.get(rec_id))

    sp = 2

    minRec = np.min(rec)     maxRec = np.max(rec)     minTarget = 0     maxTarget = 65000

    c1 = (maxTarget-minTarget) / (maxRec - minRec)     c2 = maxTarget-c1*maxRec

    rec = rec*c1+c2     rec = rec.astype(np.uint16)

    vol = sitk.GetImageFromArray(rec)     vol.SetSpacing(spnp.ones(3))     dims = np.array(list(rec.shape))     vol.SetOrigin((1-dims)0.5*sp)

    sitk.WriteImage(vol, os.path.join(dirname, "reconstruction.mhd"))

    return vol

On 12.12.22 10:46, Alexander Skorikov wrote:

Hi! Could you maybe post the full script (well the part related to reconstruction) if you couldn't resolve this problem in the end?

— Reply to this email directly, view it on GitHub https://github.com/astra-toolbox/astra-toolbox/issues/347#issuecomment-1346178806, or unsubscribe https://github.com/notifications/unsubscribe-auth/A3PJLSZY5XNOQSVVAHSEOCTWM3X75ANCNFSM6AAAAAAQ7UCV4M. You are receiving this because you authored the thread.Message ID: @.***>

askorikov commented 1 year ago

I'm not sure if I understood correctly what type of geometry you have, but there is a dedicated fanflat geometry type in Astra (https://www.astra-toolbox.com/apiref/creators.html#astra.creators.create_proj_geom). I think using it (if applicable) could solve this CUDA error.

vlibX commented 1 year ago

Hi Alexander,

I attached a picture that describes our geometry.  Basically, we have a 1D detector, placed vertically. The barrel to be imaged is sequentially translated in front of the detector.  Hence, the geometry is neither a cone beam nor a fan beam CT.

Domenico Iuso wrote the script for us. Here is his take at the memory issue:

"

I confirm my assumption: on the CUDA programming guide ( https://docs.nvidia.com/cuda/cuda-c-programming-guide/index.html#features-and-technical-specifications__technical-specifications-per-compute-capability ) it is written that there is an hard-limit on the size for the allocated 3D arrays, which is the one you are hitting on dimension 2.

It can, though, be solved, but the solution involves some programming: The usual geometries of Astra have a 2D detector and a certain number of views (the memory representation is [first detector line, views, second detector line]), while in this case there is a much higher number of views (to compensate the 1D detector) - so both the [number of rotation

"

On 12.12.22 12:55, Alexander Skorikov wrote:

I'm not sure if I understood correctly what type of geometry you have, but there is a dedicated fanflat geometry type in Astra (https://www.astra-toolbox.com/apiref/creators.html#astra.creators.create_proj_geom). I think using it (if applicable) could solve this CUDA error.

— Reply to this email directly, view it on GitHub https://github.com/astra-toolbox/astra-toolbox/issues/347#issuecomment-1346344356, or unsubscribe https://github.com/notifications/unsubscribe-auth/A3PJLS5G4G4L6TEQRULTRY3WM4HBLANCNFSM6AAAAAAQ7UCV4M. You are receiving this because you authored the thread.Message ID: @.***>

askorikov commented 1 year ago

The picture didn't get attached I guess, but I agree with Domenico (at least the first point, the second one will depend on geometry).

askorikov commented 1 year ago

Closing the issue. Let us know if the problem is still relevant.