Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

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

Open
1372484434 opened this issue Sep 17, 2023 · 19 comments
Open

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

1372484434 opened this issue Sep 17, 2023 · 19 comments

Comments

@1372484434
Copy link

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:
@1372484434
Copy link
Author

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
Copy link
Member

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
Copy link
Author

1372484434 commented Sep 17, 2023 •
edited
Loading

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)

  • 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

@AnderBiguri
Copy link
Member

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
Copy link
Author

1372484434 commented Sep 18, 2023 via email

@AnderBiguri
Copy link
Member

@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
Copy link
Author

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

@AnderBiguri
Copy link
Member

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

@1372484434
Copy link
Author

Can you see this?问题

@AnderBiguri
Copy link
Member

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
Copy link
Author

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

@AnderBiguri
Copy link
Member

AnderBiguri commented Sep 18, 2023 via email

@1372484434
Copy link
Author

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
Copy link
Member

@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
Copy link
Author

1372484434 commented Sep 18, 2023 •
edited by AnderBiguri
Loading 

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
Copy link
Member

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
Copy link
Author

ubuntu、Linux、11GB、image

@TianSong1991
Copy link

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

@AnderBiguri
Copy link
Member

@1372484434 this may be related to #496

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Assignees
No one assigned
Labels
None yet
Projects
None yet
Milestone
No milestone
Development

No branches or pull requests
3 participants
@AnderBiguri @TianSong1991 @1372484434