Out of Memory Error for 4000x4000x2096 VOLUME

[elliotaldarson1717]

Hello Kyle, Thanks for such a wonderful toolkit I have 3 questions

1. my actual voxel size is 8.6um here it shows ~34um than i read LEAP docs.. and found that leapct.set_default_volume(4.0) is kind of 4x4 binning .. i made it 1.0 and got following error

numpy.core._exceptions._ArrayMemoryError: Unable to allocate 125. GiB for an array with shape (2096, 4000, 4000) and data type float32

My Setup is Core i9 11th Gen 32GB RAM GPU: RTX4090 .... Do i need to install 128GB ram or other workaround is there.....

2 Is it possible to Preview Selected Slice ???

3. I attached the screenshot of slice ... any suggestion to remove the noise .. I mean should i use pre processing filters or its Due to FBP .. On LEAP project frontpage there is comparision between ASTA FDK and LEAP FDK LEAP FDK is cleaner .. so asking your suggestions.....

[kylechampley]

Hello. Thanks for posting this issue. The answers to your questions are listed below.

1) Yes, I know those voxels are big. I just used bigger voxels so that I could quickly reconstruct the whole volume without running out of CPU memory. I am currently working on a wrapper library for LEAP that will automatically process the data in chunks so that you don't run out of CPU RAM. I'll let you know when this is ready. In the meantime, you can follow this demo script: d13_cropping_subchunking.py which shows you how to divide the volume in smaller chunks and process each chunk sequentially.

2) Yes, you can preview a selected slice. To reconstruct one slice, just do: leapct.set_numZ(1) and then leapct.set_offsetZ(z_position), where z_position is the location of the slice you are interested in. There should be an example of this in the sample script I sent you. You can also use the parameter_sweep algorithm to sample a single slice. This algorithm is a part of the leap preprocessing algorithms.

3) For noise reduction, I recommend regularized iterative reconstruction or applying postprocessing noise filters. I do not recommend denoising the projection data. In general, regularized iterative reconstruction is the most effective method at suppressing noise, but it is very time consuming and requires a lot of RAM. The best denoising iterative algorithm is RWLS. For postprocessing, my favorite denoising method is the guided filter, but Total Variation works really well too. You should be able to find examples of all of these things in the demo_ctypes folder.

Good luck and let me know if you have any further questions!

[elliotaldarson1717]

Hello Kyle I used cropping subchunking ... so now i am able to use FBP for 150 slices.. but result is not the same in 150 slices few shows this type of artefact which is not problem in regular FBP

here is the code

Set scanner geometry

numAngles = 361 numCols = 4000 numRows = 2096 pixelWidth = 11.62/1000.0 pixelHeight = pixelWidth sod = 255.8 sdd = 345.1 tiltAngle = -0.3 leapct.set_conebeam(numAngles, numRows, numCols, pixelHeight, pixelWidth, 0.5(numRows-1), 0.5(numCols-1), leapct.setAngleArray(numAngles, 1.0*numAngles), sod, sdd)

Set reconstruction volume

leapct.set_default_volume()

leapct.set_default_volume()

""" z = leapct.z_samples() z_0 = z[469] leapct.set_default_volume() leapct.set_numZ(1) leapct.set_offsetZ(z_0)

"""

leapct.print_parameters()

Read data

dataPath = os.path.abspath(os.path.dirname('D:\XRAY\6DIMLEAP\'))

dataPath = os.path.abspath(os.path.dirname('D:\XRAY\5STONE\')) files = glob.glob(os.path.join(dataPath, '*.tif')) g = leapct.allocate_projections() for n in range(len(files)): file = files[n] anImage = np.array(imageio.imread(files[n])) g[n,:,:] = anImage[:,:]

Convert to attenuation

ROI = np.array([5, 55, 5, 55]) makeAttenuationRadiographs(leapct, g, None, None, ROI) leapct.find_centerCol(g)

"""Estimate detector tilt alphas = (np.array(range(21))-10)/10.0 f_stack = parameter_sweep(leapct, g, alphas, param='tilt', iz=469, algorithmName='FBP') leapct.display(f_stack) quit()

"""

Perform detector tilt

leapct.convert_to_modularbeam() leapct.rotate_detector(-0.3)

Reconstruct with FBP

f = leapct.allocate_volume()

g=leapct.cropProjections([1000, 1150], None, g) leapct.set_default_volume() leapct.print_parameters()

"Simulate" projection data

startTime = time.time()

Reset the volume array to zero, otherwise iterative reconstruction algorithm will start their iterations

with the true result which is cheating

x,y,z = leapct.voxelSamples(True)

f = leapct.FBP(g)

leapct.display(f)

[elliotaldarson1717]

I tried to use this for denoising https://github.com/LLNL/LEAP/blob/main/demo_leapctype/d33_reducingConeBeamArtifacts.py but i can do only 20 slices at a time BTW i have RTX4090 so i think there is some other problem with Memory .... and in those 20 slices i didnt see any Difference .. may be i am doing wrong ....

[kylechampley]

Oh, sorry, you were using the wrong part of that script. Please read the comments in that demo script file. See the section under the statement elif whichDemo == 2? That is the part you need because it steps through the volume, one chunk at a time. You will have to insert text that saves each chunk to a file to view it at the end.

For the denoising stuff, do not use the d33 script; this is for something else. Look at this script d29_filter_sequence.py for a demonstration of doing iterative reconstruction for denoising purposes. Of course, you probably don't have enough memory for this anyway.

If I were you, I'd do this. First reconstruct your whole volume as a described above and save it to disk. Then load your volume, one chunk at a time and perform denoising on it. I recommend Total Variation or Guided Filter. See the filterSequece part of the docs for an explanation.

Your GPU card is not an issue, LEAP will automatically split up large data to fit into your GPU. Your real issue is that you don't have enough CPU memory for these operations. Iterative reconstruction requires a HUGE amount of memory.

[elliotaldarson1717]

OK .. Got your point.. I read few LEAP docs and working on it ... about RAM i already ordered 128GB .. Few more question .. Does Linux Makes difference .. cant we use Virtual memory? .. once full process is done than i will focus on optimization.. for Denoising there is no enough memory i mean 32GB is not enough but cant we do it in Chunk like reconstruction .. ... Can you tell me optimal Configuration for my Data??? I will order that Config

[kylechampley]

One copy of your volume is about 128 GB, so I don't think 128 GB of RAM will be sufficient. I'd recommend at least 256 GB of RAM.

I'm not sure about the whole virtual memory thing. I don't even use this because it is likely horribly inefficient to use this in tomography.

[elliotaldarson1717]

Kyle, Thanks.. I got the whole idea and Clean result too...Somehow ASDPOCS example gave me better result than RWLS may be due to filter or filter params .. i will check with different settings

RWLS Code

filters = filterSequence(1.0e0) filters.append(TV(leapct, delta=0.01/100.0, p=1.0))

leapct.RWLS(g,f,10,filters,None,'SQS')

ASDPOCS

filters = filterSequence() filters.append(MedianFilter(leapct, 0.0, 5))

filters.append(TV(leapct, delta=0.01/100.0, p=1.0))

leapct.ASDPOCS(g,f,3,10,10,filters)

Now i have very big problem ... i modified that chunk thing but with every iteration RAM usage increases than finally CudaMalloc error than other memory error

First Iteration:

Second Iteration:

Third Iteration

......

Final Error :

Here is the Code : import sys import os import time import glob import imageio import numpy as np from leapctype import leapct = tomographicModels() from leap_preprocessing_algorithms import

Set scanner geometry

numAngles = 361 numCols = 4000 numRows = 2096 pixelWidth = 11.62/1000.0 pixelHeight = pixelWidth sod = 255.8 sdd = 345.1 tiltAngle = -0.3 leapct.set_conebeam(numAngles, numRows, numCols, pixelHeight, pixelWidth, 0.5(numRows-1), 0.5(numCols-1), leapct.setAngleArray(numAngles, 1.0*numAngles), sod, sdd)

Set reconstruction volume

leapct.set_default_volume()

leapct.set_default_volume()

leapct.print_parameters()

Read data

dataPath = os.path.abspath(os.path.dirname('D:\XRAY\6DIMLEAP\')) files = glob.glob(os.path.join(dataPath, '*.tif')) g = leapct.allocate_projections() for n in range(len(files)): file = files[n] anImage = np.array(imageio.imread(files[n])) g[n,:,:] = anImage[:,:]

Convert to attenuation

ROI = np.array([5, 55, 5, 55]) makeAttenuationRadiographs(leapct, g, None, None, ROI) leapct.find_centerCol(g)

Perform detector tilt

leapct.convert_to_modularbeam() leapct.rotate_detector(-0.3)

z = leapct.z_samples() chunkSize = leapct.get_numZ()//8 numChunks = leapct.get_numZ()//chunkSize

leapct_chunk = tomographicModels()

for n in range(numChunks): print('*****') leapct_chunk.copy_parameters(leapct)

sliceStart = n*chunkSize
sliceEnd = min(z.size-1, sliceStart + chunkSize - 1)
numZ = sliceEnd - sliceStart + 1
if numZ <= 0:
    break

leapct_chunk.set_numZ(numZ)
leapct_chunk.set_offsetZ(leapct_chunk.get_offsetZ() + z[sliceStart]-leapct_chunk.get_z0())

rowRange = leapct_chunk.rowRangeNeededForBackprojection()

g_chunk=leapct_chunk.cropProjections(rowRange, None, g)
leapct_chunk.print_parameters()
    #quit()

f_chunk = leapct_chunk.FBP(g_chunk)
leapct_chunk.display(f_chunk)

[elliotaldarson1717]

One strange thing.. that sample works well with Simulated phantom .. i tried to increate numCol=512 to 1024 ... memory usage is the same after every iteration .. somehow its not working with tiff files or my tiff files

[kylechampley]

RWLS and ASDPOCS have regularization parameters that need to be tuned for your specific data. Thus, one cannot say that one algorithm works and another does not until these parameters are tuned. FBP is great because it gives predictable results and is an excellent algorithm for a non-expert, but iterative reconstruction requires one to understand how the algorithms work in order to get the best performance.

In addition, one cannot simply reconstruct a subset of cone-beam slices without correcting for the out of bound slices. I am building automated routines for this, but they are not ready yet. Until these are ready or you figure out how to do this yourself, I would not do iterative reconstruction if I were you.

Instead, I would try postprocess denoising. Some good options are: GuidedFilter 3D Median Filter Total Variation

[elliotaldarson1717]

Hello Kyle...that problem soved i mean in FBP no double edges.. now i have to run that chunk script 8 times... for loop gives memory error in 3rd iteration FBP only ... the sample works great.. i uploaded Memory photos .. is it the problem for my side??? if you see the source code i just replaced Phantom code with tiff code....Here i dont want any correction just For loop not working... regarding Denoising i will try those filters with Different parameters

Result of FPB only ... my main problem of Double Edges solved .... Thank you Very much ...

[elliotaldarson1717]

@kylechampley

Everything sorted out .. Total Variation worked ... i have one weird problem asusal.. Intensity changes in slices.... here is the screenshot .. you can see few slices are dark few are bright .. that creates problem in segmentation ...

it happens randomly here is the code .. am i doing anything wrong??

import sys import os import time import glob import imageio import numpy as np from leapctype import leapct = tomographicModels() from leap_preprocessing_algorithms import import cv2

Set scanner geometry

numAngles = 361 numCols = 2808 numRows = 1902 pixelWidth = 11.62/1000.0 pixelHeight = pixelWidth sod = 275.8 sdd = 365.1 tiltAngle = -0.3 leapct.set_conebeam(numAngles, numRows, numCols, pixelHeight, pixelWidth, 0.5(numRows-1), 0.5(numCols-1), leapct.setAngleArray(numAngles, 1.0*numAngles), sod, sdd)

leapct.set_default_volume()

leapct.set_default_volume(1.0)

""" z = leapct.z_samples() z_0 = z[469] leapct.set_default_volume() leapct.set_numZ(1) leapct.set_offsetZ(z_0)

"""

leapct.print_parameters()

Read data

dataPath = os.path.abspath(os.path.dirname('D:\XData\6DIMLEAP\')) files = glob.glob(os.path.join(dataPath, '*.tif')) g = leapct.allocate_projections() for n in range(len(files)): file = files[n] anImage = np.array(imageio.imread(files[n])) g[n,:,:] = anImage[:,:]

Convert to attenuation

ROI = np.array([5, 55, 5, 55])

makeAttenuationRadiographs(leapct, g, None, None, ROI)

Perform detector tilt

leapct.convert_to_modularbeam() leapct.rotate_detector(-0.3)

z = leapct.z_samples() chunkSize = leapct.get_numZ()//4 numChunks = leapct.get_numZ()//chunkSize

leapct_chunk = tomographicModels() print('*****') leapct_chunk.copy_parameters(leapct)

n=0

output_dir="I:\TEMP\" output_file=os.path.join(output_dir, 'Volume' + str(n) +'.tiff')

sliceStart = n*chunkSize sliceEnd = min(z.size-1, sliceStart + chunkSize - 1) numZ = sliceEnd - sliceStart + 1 if numZ <= 0: quit()

leapct_chunk.set_numZ(numZ) leapct_chunk.set_offsetZ(leapct_chunk.get_offsetZ() + z[sliceStart]-leapct_chunk.get_z0()) rowRange = leapct_chunk.rowRangeNeededForBackprojection()

g_chunk=leapct_chunk.cropProjections(rowRange, None, g)

leapct_chunk.print_parameters()

f_chunk = leapct_chunk.FBP(g_chunk)

cv2.imwrite("C:\asdfasdf.tif",normalized_image)

print("BEFORE FILTER") leapct_chunk.set_numTVneighbors(26) leapct_chunk.MedianFilter(f_chunk, 0.5,5) print("MEDIAN DONE") leapct_chunk.diffuse(f_chunk, delta=0.02/20.0,numIter=20) print("DIFFUSE DONE")

print("SAVING VOLUME")

leapct_chunk.save_volume(output_file, f_chunk,sliceStart)

[kylechampley]

Are you sure the intensity is changing? This may just be an artifact of how Windows shows these images in their file explorer. Often these images are displayed from their min to max value and this can certainly change, slice-by-slice, but the average intensity should be very continuous from slice-to-slice. Can you bring these slice images into a 3D viewer to verify?

[elliotaldarson1717]

Kyle here is the snapshot

I am using ORS Dragonfly for the 3d Mesh Creation .. you can see the missing layer in 3d view those are not missing but low intensity ...

[kylechampley]

That's really weird; I have no idea why that would happen. I am traveling and have limited computer access, so I will try this when I get back. I did try these filters with the d13 demo script and it worked fine.

A few comments. It looks like you're using the default value for "delta" in the TV denoising (diffuse). Did you optimize this value for your data?

[elliotaldarson1717]

Thanks, Tonight I planned to optimize denoising. If you can suggest delta range with steps... It will be helpful

[kylechampley]

There is a good description of the effect of the delta parameter in the documentation of the TV filter sequence.

[kylechampley]

@elliotaldarson1717 so sorry for taking so long to get back to you. I looked at your script and it looks like your data has been cropped or something because the data dimensions are different. You could be getting variable brightness because you are cropping off some of the data that contains attenuation values. CT requires untruncated projections. Anyway, I took the code you have above and made some modifications and it seems to work now. I do feel that you over-smoothed your data, but that's for you to decide. Anyway, the code below should work.

import sys
import os
import time
import glob
import imageio
import numpy as np
from leapctype import *
leapct = tomographicModels()
from leap_preprocessing_algorithms import *

# Set scanner geometry
numAngles = 361
numCols = 4000
numRows = 2096
pixelWidth = 11.62/1000.0
pixelHeight = pixelWidth
sod = 255.8
sdd = 345.1
tiltAngle = -0.3
leapct.set_conebeam(numAngles, numRows, numCols, pixelHeight, pixelWidth, 0.5*(numRows-1), 0.5*(numCols-1), leapct.setAngleArray(numAngles, 1.0*numAngles), sod, sdd)

# Set reconstruction volume
leapct.set_default_volume()
#leapct.set_default_volume(4.0)
leapct.print_parameters()

# Read data
dataPath = os.path.abspath(os.path.dirname(__file__))
files = glob.glob(os.path.join(dataPath, '*.tif'))
g = leapct.allocate_projections()
for n in range(len(files)):
    file = files[n]
    anImage = np.array(imageio.imread(files[n]))
    g[n,:,:] = anImage[:,:]

# Convert to attenuation
ROI = np.array([5, 55, 5, 55])
makeAttenuationRadiographs(leapct, g, None, None, ROI)
leapct.find_centerCol(g)

# Perform detector tilt
leapct.convert_to_modularbeam()
leapct.rotate_detector(-0.3)

z = leapct.z_samples()
chunkSize = leapct.get_numZ()//16
numChunks = leapct.get_numZ()//chunkSize

output_dir = os.path.join(dataPath, 'recon')
if not os.path.exists(output_dir):
    os.makedirs(output_dir)
output_file=os.path.join(output_dir, 'Volume' + '.tiff')

leapct_chunk = tomographicModels()
for n in range(numChunks):
    print('*************************************************************')
    leapct_chunk.copy_parameters(leapct)

    sliceStart = n*chunkSize
    sliceEnd = min(z.size-1, sliceStart + chunkSize - 1)
    numZ = sliceEnd - sliceStart + 1
    if numZ <= 0:
        quit()

    leapct_chunk.set_numZ(numZ)
    leapct_chunk.set_offsetZ(leapct_chunk.get_offsetZ() + z[sliceStart]-leapct_chunk.get_z0())
    rowRange = leapct_chunk.rowRangeNeededForBackprojection()

    g_chunk=leapct_chunk.cropProjections(rowRange, None, g)

    leapct_chunk.print_parameters()

    f_chunk = leapct_chunk.FBP(g_chunk)

    print("BEFORE FILTER")
    leapct_chunk.set_numTVneighbors(26)
    leapct_chunk.MedianFilter(f_chunk, 0.5,5)
    print("MEDIAN DONE")
    leapct_chunk.diffuse(f_chunk, delta=0.02/20.0,numIter=20)
    print("DIFFUSE DONE")

    print("SAVING VOLUME")

    leapct_chunk.save_volume(output_file, f_chunk,sliceStart)

[elliotaldarson1717]

Thanks Kyle .. somehow that out of memory issue solved.. and images looks much better in 3d view... Just a question you wrote "you over-smoothed your data" you mean My projections are oversmooth or you are talking about Median and TV Filters i used in that script??? that filters was for trial to see the effect... if there is any problem in Projection please let me know... i will correct it...

[kylechampley]

Yes, I was just referring to the parameters used in the Median and TV filters.

[hws203]

For using the large size projection image for tomography, I think that binning(2x2) option is best for saving time and memory based on my experience.

[elliotaldarson1717]

@hws203

For using the large size projection image for tomography, I think that binning(2x2) option is best for saving time and memory based on my experience. @hws203 , My application is to find Crack in Gemstone so for Good Spatial Resolution i thought its better not to bin... actually pixel size is ~8.6um i need actual voxel resolution of around 15-20um .. but i read somewhere that if pixel resolution should be higher than what you actual want .. May be we loose details in filtering etc.. .. BTW in which application you are working ???

[elliotaldarson1717]

@hws203 , @kylechampley ... Any suggestion on Projection Image improvement for below kind of noise(Like LOW DOSE HIGH DOSE , More frame averaging, this data is taken on 12 bit OLD detector .. ) ... it prevent me from auto segmentation

[elliotaldarson1717]

@kylechampley

Its not good to flood that Discuttion thread with problem solution so asking here.. in My Application(Single material with some inclusions and cracks only) 8 bit image/slice output is more than enough and desirable.. right now output images are in 32bit float i tried to modify the save_data method of leapctype.py as following .. but both resulted in inferior results here is the modified code of save_data method... check comment out .. i was trying diffrent things...

def save_data(self, fileName, x, T=1.0, offset_0=0.0, offset_1=0.0, offset_2=0.0, sequence_offset=0): """Save 3D data to file (tif sequence, nrrd, or npy)""" volFilePath, dontCare = os.path.split(fileName) if os.path.isdir(volFilePath) == False or os.access(volFilePath, os.W_OK) == False: print('Folder to save data either does not exist or not accessible!') return False

    if has_torch == True and type(x) is torch.Tensor:
        x = x.cpu().detach().numpy()

    if fileName.endswith('.npy'):
        np.save(fileName, x)
        return True
    elif fileName.endswith('.nrrd'):
        try:
            import nrrd

            # https://pynrrd.readthedocs.io/en/latest/examples.html
            header = {'units': ['mm', 'mm', 'mm'], 'spacings': [T, T, T], 'axismins': [offset_0, offset_1, offset_2], 'thicknesses': [T, T, T],}
            nrrd.write(fileName, x, header)
            return True
        except:
            print('Error: Failed to load nrrd library!')
            print('To install this package do: pip install pynrrd')
            return False
    elif fileName.endswith('.tif') or fileName.endswith('.tiff'):
        try:
            #from PIL import Image
            import imageio
            baseName, fileExtension = os.path.splitext(fileName)
            if len(x.shape) <= 2:
                im = x
                #im.save(baseName + '_' + str(int(i)) + fileExtension)
                imageio.imwrite(baseName + fileExtension, im)
                print("IFFFFF CALLEDDDD"+baseName + fileExtension)
            else:
                for i in range(x.shape[0]):
                    import cv2
                    #trialslice = np.squeeze(x[i,:,:])
                    #mat11 =cv2.Mat(trialslice)

                    #normalized_image = cv2.normalize(mat11, None,255,0, cv2.NORM_MINMAX,cv2.CV_8U) 
                    #normalized_image = cv2.equalizeHist(normalized_image)  
                    #cv2.imwrite(baseName + '_' + str(int(i)+sequence_offset) + fileExtension,normalized_image)

                    trialslice = np.squeeze(x[i,:,:])
                    trialsliceAstra2 =trialslice.copy()
                    #if np.min(trialsliceAstra2) <0:
                    #    trialsliceAstra2 += (-np.min(trialsliceAstra2))

                    #trialsliceAstra2[trialsliceAstra2 < 0] = 0
                    #trialsliceAstra2 /= np.max(trialsliceAstra2)
                    #trialsliceAstra2 = np.round(trialsliceAstra2 * 65535).astype(np.uint16)
                    #imageio.imwrite(baseName + '_' + str(int(i)+sequence_offset) + fileExtension, trialsliceAstra2)

                    im = x[i,:,:]
                    #im.save(baseName + '_' + str(int(i)) + fileExtension)
                    imageio.imwrite(baseName + '_' + str(int(i)+sequence_offset) + fileExtension, im)

            return True

        except:
            #print('Error: Failed to load PIL library!')
            #print('To install this package do: pip install Pillow')
            print('Error: Failed to load imageio library!')
            print('To install PIL do: pip install imageio')
            return False
    else:
        print('Error: must be a tif, npy, or nrrd file!')
        return False