The projected image is too large, causing the reconstruction to be unresponsive
1372484434 opened this issue · 19 comments
Expected Behavior
The projected image is too large, causing the reconstruction to be unresponsive
Actual Behavior
Code to reproduce the problem (If applicable)
Specifications
- MATLAB version:
- OS:
- CUDA version:
The projected image is too large, specifically (50,1024,1024), using fdk sart and other reconstruction does not respond, and it has been stuck, but other data sets can be used.
Hi @1372484434 , The projected image you suggest is certainly not very large, and also not very large for TIGRE at all.
Can you please give me information on your system, and also code that would reproduce the effect?
Below is part of my code. When executed here, it stops moving.“ recon = algs.sart(projections, geo, data['train']['angles'], niter=100, gpuids=gpuids)”.When I run other datasets, such as the projected image is (50, 256, 256), etc., it can run normally.
code：
if data['randomAngle'] is False: data['train'] = {'angles': np.linspace(0, data['totalAngle'] / 180 * np.pi, data['numTrain']+1)[:-1] + data['startAngle']/ 180 * np.pi} else: data['train'] = {'angles': np.sort(np.random.rand(data['numTrain']) * data['totalAngle'] / 180 * np.pi) + data['startAngle']/ 180 * np.pi} projections = tigre.Ax(np.transpose(img, (2, 1, 0)).copy(), geo, data['train']['angles'])[:, ::-1, :] plt.figure() plt.imshow(projections[0, :, :]) plt.show() gpuids = gpu.getGpuIds(listGpuNames[0]) print(gpuids) # geo['angles'] = data['train']['angles'] #recon = algs.fdk(projections, geo, data['train']['angles'], gpuids=gpuids) #recon = algs.ossart(projections, geo, data['train']['angles'], 100, blocksize=20, gpuids=gpuids) recon = algs.sart(projections, geo, data['train']['angles'], niter=100, gpuids=gpuids)
The following is the configuration information：
Parameters to setup CT simulation.
dataset
numTrain: 50
numVal: 50
System configuration
DSD: 1500.0 # Distance Source Detector (mm)
DSO: 1000.0 # Distance Source Origin (mm)
Detector parameters
nDetector: # Number of pixels (px)
- 1024
- 1024
dDetector: # Size of each pixel (mm) - 1.0
- 1.0
Image parameters
nVoxel: # Number of voxels (vx)
- 512
- 512
- 463
dVoxel: # size of each voxel (mm) - 0.625
- 0.625
- 1.0
Offsets
offOrigin: # Offset of image from origin (mm)
- 0 # x direction
- 0 # y direction
- 0 # z direction
offDetector: # Offset of Detector (only in two direction) (mm) - 0 # u direction
- 0 # v direction
Auxiliary
accuracy: 0.5 # Accuracy of FWD proj (vx/sample)
Mode
mode: cone # X-ray source mode parallel/cone
filter: null
Angles
totalAngle: 180.0 # Total angle (degree)
startAngle: 0.0 # Start angle (degree)
randomAngle: False
CT
convert: True
rescale_slope: 1.0
rescale_intercept: 0.0
normalize: True
Noise
noise: 0
Please produce a minimal example that reproduces the code please, something that I can copy-paste and it would intermediately run and reproduce your issue.
@1372484434 Sorry, if I were to copy paste this code, it wold not work because it is not complete. Can yo uplease make a complete example?
It is complete, you can try to run it, can you download the attachment I sent?
@1372484434 It is not properly formatted, does not have imports, etc. The first error would be data not defined. It is not complete.
Can you see this?
No, that file has not reached my email, but you are responding in github, so maybe it doesn't allow you to add files.
Can you provide me with an email to send to you? Github does not allow sending this 128MB email
I tried mine. When the projection angle is (19, 1024, 1024), SART can iterate normally. When the projection angle is greater than 19, it will not run. If I don’t provide the data for you to run it, , you may not experience my failure
@1372484434 have you tried by simulating projections with TIGRE?
Are you saying that its your data values, and not the size of the data, that makes TIGRE fail?
I doubt that is the case, but if it is, then you know where to look for errors: your data.
import tigre
from tigre.utilities.geometry import Geometry
import yaml
from tigre.algorithms import sart
import scipy.io
import scipy.ndimage.interpolation
import numpy as np
import matplotlib.pyplot as plt
import tigre.utilities.gpu as gpu
#os.environ["CUDA_VISIBLE_DEVICES"] = "0"
listGpuNames = gpu.getGpuNames()
if len(listGpuNames) == 0:
print("Error: No gpu found")
else:
for id in range(len(listGpuNames)):
print("{}: {}".format(id, listGpuNames[id]))
gpuids = gpu.getGpuIds(listGpuNames[0])
print(gpuids)
I did use TIGRE to simulate projection.'
def main():
matPath = f'./img.mat'#./dataGenerator/raw/chest/img.mat
configPath = f'./config.yml'#./dataGenerator/raw/chest/config.yml
generator(matPath, configPath, True)
class ConeGeometry_special(Geometry):
def __init__(self, data):
Geometry.__init__(self)
# VARIABLE DESCRIPTION UNITS
# -------------------------------------------------------------------------------------
self.DSD = data['DSD'] / 1000 # Distance Source Detector (m)1.5米
self.DSO = data['DSO'] / 1000 # Distance Source Origin (m)1米
# Detector parameters
self.nDetector = np.array(data['nDetector']) # number of pixels (px)[256,256]
self.dDetector = np.array(data['dDetector']) / 1000 # size of each pixel (m)[0.0015,0.0015]
self.sDetector = self.nDetector * self.dDetector # total size of the detector (m)[0.384,0.384]
# Image parameters
self.nVoxel = np.array(data['nVoxel'][::-1]) # number of voxels (vx)[128,128,128]
self.dVoxel = np.array(data['dVoxel'][::-1]) / 1000 # size of each voxel (m)[0.001,0.001,0.001]
self.sVoxel = self.nVoxel * self.dVoxel # total size of the image (m)[0.128,0.128,0.128]
# Offsets
self.offOrigin = np.array(data['offOrigin'][::-1]) / 1000 # Offset of image from origin (m)图像与原点的偏移[0.,0.,0.]
self.offDetector = np.array(
[data['offDetector'][1], data['offDetector'][0], 0]) / 1000 # Offset of Detector (m)[0.,0.,0.]
# Auxiliary
self.accuracy = data['accuracy'] # Accuracy of FWD proj (vx/sample) # noqa: E5010.5
# Mode
self.mode = data['mode'] # parallel, cone ...
self.filter = data['filter']
def convert_to_attenuation(data: np.array, rescale_slope: float, rescale_intercept: float):
"""
CT scan is measured using Hounsfield units (HU). We need to convert it to attenuation.
The HU is first computed with rescaling parameters:
HU = slope * data + intercept
Then HU is converted to attenuation:
mu = mu_water + HU/1000x(mu_water-mu_air)
mu_water = 0.206
mu_air=0.0004
Args:
data (np.array(X, Y, Z)): CT data.
rescale_slope (float): rescale slope.
rescale_intercept (float): rescale intercept.
Returns:
mu (np.array(X, Y, Z)): attenuation map.
"""
HU = data * rescale_slope + rescale_intercept
mu_water = 0.206
mu_air = 0.0004
mu = mu_water + (mu_water - mu_air) / 1000 * HU
# mu = mu * 100
return mu
def loadImage(dirname, nVoxels, convert, rescale_slope, rescale_intercept, normalize=True):
"""
Load CT image.
"""
test_data = scipy.io.loadmat(dirname)
# Loads data in F_CONTIGUOUS MODE (column major), convert to Row major
image_ori = test_data["img"].astype(np.float32)
if convert:
print('Convert from HU to attenuation') ##在这里不需要将HU转为衰减系数
image = convert_to_attenuation(image_ori, rescale_slope, rescale_intercept)
else:
image = image_ori
image_max = np.max(image)
image_min = np.min(image)
image_mean = np.mean(image)
print('Range of CT image is [%f, %f], mean: %f' % (image_min, image_max, image_mean))
if normalize and image_min !=0 and image_max != 1:
print('Normalize range to [0, 1]')
image = (image - image_min) / (image_max - image_min)
return image
def generator(matPath, configPath, show=False):
"""
Generate projections given CT image and configuration.
"""
# Load configuration
with open(configPath, 'r') as handle:
data = yaml.safe_load(handle)
# Load CT image
geo = ConeGeometry_special(data)
img = loadImage(matPath, data['nVoxel'], data['convert'],
data['rescale_slope'], data['rescale_intercept'], data['normalize'])
plt.figure()
plt.imshow(img[:,:,0])
plt.show()
# Generate training images
data['train'] = {'angles': np.linspace(0, data['totalAngle'] / 180 * np.pi, data['numTrain']+1)[:-1] + data['startAngle']/ 180 * np.pi}
projections = tigre.Ax(np.transpose(img, (2, 1, 0)).copy(), geo, data['train']['angles'])[:, ::-1, :]
#projections.tofile('./abdomen.raw')
plt.figure()
plt.imshow(projections[0, :, :])
plt.show()
#recon = algs.fdk(projections, geo, data['train']['angles'], gpuids=gpuids)
#recon = algs.ossart(projections, geo, data['train']['angles'], 100, blocksize=20, gpuids=gpuids)
recon = sart(projections, geo, data['train']['angles'], niter=100, gpuids=gpuids)
recon.tofile('./abdomen.raw')
print(recon.shape)
plt.figure()
plt.imshow(recon[0, :, :])
plt.show()
if __name__ == '__main__':
main()
Apologies, I don't want to sound rude, but I am very very busy. I am willing to take time to help you debug your code, as I am 99.9% sure this is a problem in your side, not in TIGRE, as I do routinely reconstruct images bigger than what you suggest.
However, as I am very busy, I need you to take a couple of hours, and reduce your code the to minimum. Eg. Hounsfield units are irrelevat to the error. Or this code does not work without the variable data.
If you take your time to make this code a self-suficient shortest possible code that I can copy paste and run, I will be able to help you, otherwise, you may need to wait a month or two until I have time to rewrite it to figure out what could be wrong.
Please also do provide information about your system. e.g. RAM size, OS, etc.
ubuntu、Linux、11GB、
@1372484434 I think you can put your data in the Google Drive link so that many people can help you find the problem.
@1372484434 this may be related to #496