CERN / TIGRE

TIGRE: Tomographic Iterative GPU-based Reconstruction Toolbox
BSD 3-Clause "New" or "Revised" License
582 stars 191 forks source link

The projected image is too large, causing the reconstruction to be unresponsive #486

Open 1372484434 opened 1 year ago

1372484434 commented 1 year ago

Expected Behavior

The projected image is too large, causing the reconstruction to be unresponsive

Actual Behavior

Code to reproduce the problem (If applicable)

Specifications

1372484434 commented 1 year ago

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.

AnderBiguri commented 1 year ago

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?

1372484434 commented 1 year ago

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)`

image 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)

Image parameters

nVoxel: # Number of voxels (vx)

Offsets

offOrigin: # Offset of image from origin (mm)

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

AnderBiguri commented 1 year ago

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 commented 1 year ago

This is the smallest code I can run, 128m because the data is a little big, the code is very concise

发自我的iPhone

------------------ Original ------------------ From: Biguri @.> Date: Sun,Sep 17,2023 8:32 PM To: CERN/TIGRE @.> Cc: 1372484434 @.>, Mention @.> Subject: Re: [CERN/TIGRE] The projected image is too large, causing thereconstruction to be unresponsive (Issue #486)

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.

— Reply to this email directly, view it on GitHub, or unsubscribe. You are receiving this because you were mentioned.Message ID: @.***>

                    从QQ邮箱发来的超大附件    

                                                                    demo.zip                             (128.7MB, 2023年10月18日 14:37)                                                                           进入下载页面                          :https://wx.mail.qq.com/ftn/download?func=3&k=c79b6d38fce80b7ff5e81a3835623833d14d3f30376238331d4b4c4c075b0e005d07190000000a1e5b0207011a540d0a5e49555e0e510d525957050c015230335c015957191851432cc97fd6ba967012d64d4903d69f0ebc394b9a5a71&key=c79b6d38fce80b7ff5e81a3835623833d14d3f30376238331d4b4c4c075b0e005d07190000000a1e5b0207011a540d0a5e49555e0e510d525957050c015230335c015957191851432cc97fd6ba967012d64d4903d69f0ebc394b9a5a71&code=8d487b83&from=
AnderBiguri commented 1 year ago

@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?

1372484434 commented 1 year ago

It is complete, you can try to run it, can you download the attachment I sent?

AnderBiguri commented 1 year ago

@1372484434 It is not properly formatted, does not have imports, etc. The first error would be data not defined. It is not complete.

1372484434 commented 1 year ago

Can you see this?问题

AnderBiguri commented 1 year ago

No, that file has not reached my email, but you are responding in github, so maybe it doesn't allow you to add files.

1372484434 commented 1 year ago

Can you provide me with an email to send to you? Github does not allow sending this 128MB email

AnderBiguri commented 1 year ago

I would much prefer if you could add a simple example in code here. I don't need your data or your code, just take the TIGRE example, and make it fail like it happens to you.

Ander

On Mon, 18 Sept 2023 at 09:30, 1372484434 @.***> wrote:

Can you provide me with an email to send to you? Github does not allow sending this 128MB email

— Reply to this email directly, view it on GitHub https://github.com/CERN/TIGRE/issues/486#issuecomment-1722966055, or unsubscribe https://github.com/notifications/unsubscribe-auth/AC2OENEWMV5ZA3RX7UYGRK3X3ABC5ANCNFSM6AAAAAA43MBBXY . You are receiving this because you commented.Message ID: @.***>

1372484434 commented 1 year ago

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

AnderBiguri commented 1 year ago

@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.

1372484434 commented 1 year ago
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()
AnderBiguri commented 1 year ago

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.

1372484434 commented 1 year ago

ubuntu、Linux、11GB、image

TianSong1991 commented 1 year ago

@1372484434 I think you can put your data in the Google Drive link so that many people can help you find the problem.

AnderBiguri commented 1 year ago

@1372484434 this may be related to #496