Reconstrution with AAPM Low Dose Challage Dataset using helical code

[lixing0810]

Hello, thanks for you amazing job. I want to recontruct CT images from AAPM Low Dose Challage Dataset. But it didn't work. Here is my code: `import time import sys sys.path.append('/home/abc/data_project/LEAP-main/src') from leapctype import leapct = tomographicModels() import scipy.io as scio import h5py numCols = 736 numTurns = 10 numAngles = 35680 numRows = 64 pixelHeight = 1.0947 pixelWidth = 1.2858 centerRow = 0.5(numRows+1) centerCol = 369.625 phis = leapct.setAngleArray(numAngles, 360.0*numTurns) sod = 595 sdd = 1085.6 pitch = 0.8 leapct.set_conebeam(numAngles, numRows, numCols, pixelHeight, pixelWidth, centerRow, centerCol, phis, sod, sdd) leapct.set_normalizedHelicalPitch(pitch) leapct.set_curvedDetector()

numX = 512 numY = 512 numZ = 211 voxelWidth = 0.6934 voxelHeight = 0.6934 offsetX = 0 offsetY = 0 offsetZ = 0 leapct.set_volume(numX, numY, numZ, voxelWidth, voxelHeight, offsetX, offsetY, offsetZ) leapct.set_diameterFOV(512) leapct.print_parameters()

g = np.ascontiguousarray(np.transpose(np.asarray(h5py.File('/home/abc/data_project/quater_dose_506.mat','r')['dicom_matrix']).astype(np.float32),[2,1,0]))

f = leapct.allocateVolume() startTime = time.time() leapct.FBP(g,f) print('Reconstruction Elapsed Time: ' + str(time.time()-startTime))` The reconstruction CT images like below:

Here is the dataset link: 'https://aapm.app.box.com/s/eaw4jddb53keg1bptavvvd1sf4x3pe9h/folder/145431584872' Thanks a lot!

[kylechampley]

Can you send me a link to the CT geometry description?

[lixing0810]

Can you send me a link to the CT geometry description?

@kylechampley of course, its offical provided geomestry is below: L506_4L_100kv_quarterdose1.txt

[lixing0810]

There are 526 CT images after reconstruction and still provided in the offical website. And I found that the voxelWidth and voxelHeight = 0.7422 from pixelspacing key of CT DICOM image files. Thanks a lot!

[kylechampley]

I figured it out. You were quite close. You had the angles specified incorrectly and I had to change the normalized helical pitch to -0.8. The number of reconstruction slices you specified will only reconstruct part of the patient; you may want to increase numZ. Below is my script.

import time
import sys
import glob
from leapctype import *
leapct = tomographicModels()
import scipy.io as scio
import h5py
numCols = 736
numAngles = 51326
numAngles_perRotation = 2304
numRows = 64
pixelHeight = 1.0947
pixelWidth = 1.2858
centerRow = 0.5*(numRows+1)
centerCol = 369.625
phis = leapct.setAngleArray(numAngles, 360.0*numAngles/numAngles_perRotation)
sod = 595
sdd = 1085.6
pitch = -0.8
leapct.set_conebeam(numAngles, numRows, numCols, pixelHeight, pixelWidth, centerRow, centerCol, phis, sod, sdd)
leapct.set_normalizedHelicalPitch(pitch)
leapct.set_curvedDetector()

numX = 512
numY = 512
numZ = 211
voxelWidth = 0.6934
voxelHeight = 0.6934
offsetX = 0
offsetY = 0
offsetZ = 0
leapct.set_volume(numX, numY, numZ, voxelWidth, voxelHeight, offsetX, offsetY, offsetZ)
leapct.set_diameterFOV(numX*voxelWidth)
#leapct.set_default_volume()
leapct.print_parameters()

dataPath = r'C:\Users\champley\Downloads\DICOM-CT-PD_QD\DICOM-CT-PD_QD'
fileName = 'L506_4L_100kv_quarterdose1.*.dcm'

files = glob.glob(os.path.join(dataPath, fileName))
g = leapct.allocate_projections()
for n in range(numAngles):
    file = files[n]
    print('reading ', file)
    anImage = np.array(imageio.imread(file)).T
    g[n,:,:] = anImage[:,:]

f = leapct.allocateVolume()
startTime = time.time()
leapct.FBP(g,f)
print('Reconstruction Elapsed Time: ' + str(time.time()-startTime))
leapct.display(f)

[hws203]

Below is my running result from above script. And I can see that I have to reduce the numAngles = 51326//2, because leapct.allocateVolume() returns error message which comes from my GPU memory shortage for full scan images. I wonder that normalize pitch should be always negative sign or there is some other rule. Is it related to CW or CCW scan direction? And I hope that leapct package can handle such kind of large memory usage case as like pushing step by step image into GPU memory from reading file or host memory.

[lixing0810]

I figured it out. You were quite close. You had the angles specified incorrectly and I had to change the normalized helical pitch to -0.8. The number of reconstruction slices you specified will only reconstruct part of the patient; you may want to increase numZ. Below is my script.

import time
import sys
import glob
from leapctype import *
leapct = tomographicModels()
import scipy.io as scio
import h5py
numCols = 736
numAngles = 51326
numAngles_perRotation = 2304
numRows = 64
pixelHeight = 1.0947
pixelWidth = 1.2858
centerRow = 0.5*(numRows+1)
centerCol = 369.625
phis = leapct.setAngleArray(numAngles, 360.0*numAngles/numAngles_perRotation)
sod = 595
sdd = 1085.6
pitch = -0.8
leapct.set_conebeam(numAngles, numRows, numCols, pixelHeight, pixelWidth, centerRow, centerCol, phis, sod, sdd)
leapct.set_normalizedHelicalPitch(pitch)
leapct.set_curvedDetector()

numX = 512
numY = 512
numZ = 211
voxelWidth = 0.6934
voxelHeight = 0.6934
offsetX = 0
offsetY = 0
offsetZ = 0
leapct.set_volume(numX, numY, numZ, voxelWidth, voxelHeight, offsetX, offsetY, offsetZ)
leapct.set_diameterFOV(numX*voxelWidth)
#leapct.set_default_volume()
leapct.print_parameters()

dataPath = r'C:\Users\champley\Downloads\DICOM-CT-PD_QD\DICOM-CT-PD_QD'
fileName = 'L506_4L_100kv_quarterdose1.*.dcm'

files = glob.glob(os.path.join(dataPath, fileName))
g = leapct.allocate_projections()
for n in range(numAngles):
    file = files[n]
    print('reading ', file)
    anImage = np.array(imageio.imread(file)).T
    g[n,:,:] = anImage[:,:]

f = leapct.allocateVolume()
startTime = time.time()
leapct.FBP(g,f)
print('Reconstruction Elapsed Time: ' + str(time.time()-startTime))
leapct.display(f)

Thanks for your scripts and my reconstruction just like @hws203 showed. And I also wanna to know why the helical pitch to nagative value and if leapct package can deal with the these parameters like StartReconPos and StartEndPos ... to ensure the CT images displaying at normal angles, just like below: not ours like this:

Many thanks for your great job!

[kylechampley]

The LEAP image is rotated because the correct starting angle was not given. I think if you add about 120 to the array of angles, it will be rotated similarly.

The negative value on the helical pitch just signifies the translation of the bed position with respect to the rows of the detector. If you flip the rows of each projection, you can flip the sign of the helical pitch.

I can't exactly be sure what they mean by StartReconPos, but it likely has to do with what slices are reconstructed. To change this with LEAP change the offsetZ and numZ parameters in the set_volume function.

[lixing0810]

@kylechampley Many thx!