3D Projection of phantoms with shapes of elliptical cylinder and cuboid don't have tilting rectangular shapes whatever Euler angles set

[lokeaichirou]

Code that randomly generates 3D phantoms including elliptical cylinders only

# Import the necessary libraries.
# ~~~python
import os
import numpy as np
import matplotlib.pyplot as plt
import tifffile
import random
from copy import deepcopy
import numpy as np

from tomophantom import TomoP3D

# Define the path and name of the `.dat` file that will contain the phantoms model description. If the file exists, the new models will be appended.python

path_1 = '/home/mli/tomograms/pycharm_demos/single_density/without_Gaussian_shapes/normal_test/without_noise/1024*1024_with_oversampling=2/Generate_Training_Data/Phantom3DLibrary_For_att_without_Gaussian_3.dat'
path_2 = '/home/mli/tomograms/pycharm_demos/single_density/without_Gaussian_shapes/normal_test/without_noise/1024*1024_with_oversampling=2/Generate_Training_Data/Phantom3DLibrary_For_phase_without_Gaussian_3.dat'

# Generate the models:

number_of_phantoms = 1000

# components = ['gaussian','paraboloid', 'ellipsoid', 'cone', 'cuboid', 'elliptical_cylinder']
# the unused components are not uniform in density.

components = ['elliptical_cylinder']
mu = [0.328]
delta = [481]

for n in range(1, number_of_phantoms + 1):
    num_of_components = random.randint(1, 5)

    objects_1 = []
    objects_2 = []
    component_subset = random.choices(components, k=random.randint(1, 3))
    for i in range(num_of_components):
        obj_1 = random.choice(component_subset)
        obj_2 = deepcopy(obj_1)
        upper_lim_half_width = random.choice([0.5])

        c0 = 0
        m0 = mu[c0]
        d0 = delta[c0]

        if obj_1 == 'elliptical_cylinder' or 'cuboid':

            p_or_n_x = random.choices([-1, 1], k=1)
            p_or_n_y = random.choices([-1, 1], k=1)
            p_or_n_z = random.choices([-1, 1], k=1)

            x0 = round(random.uniform(-1.0*p_or_n_x[0], -0.75*p_or_n_x[0]), 2)
            y0 = round(random.uniform(-1.0*p_or_n_y[0], -0.75*p_or_n_y[0]), 2)
            z0 = round(random.uniform(-1.0*p_or_n_z[0], -0.75*p_or_n_z[0]), 2)
            a = round(random.uniform(0.3, upper_lim_half_width), 3)
            b = round(random.uniform(0.3, upper_lim_half_width), 3)
            c = round(random.uniform(0.3, upper_lim_half_width), 3)
            # alpha = round(random.uniform(-180, 180), 2)
            alpha = 125
            # beta = round(random.uniform(-180, 180), 2)
            beta = 70
            # gamma = round(random.uniform(-180, 180), 2)
            gamma = 50
            print('ellip')

        else:
            x0 = round(random.uniform(-0.5, 0.5), 2)
            y0 = round(random.uniform(-0.5, 0.5), 2)
            z0 = round(random.uniform(-0.5, 0.5), 2)
            a = round(random.uniform(0.01, upper_lim_half_width), 3)
            b = round(random.uniform(0.01, upper_lim_half_width), 3)
            c = round(random.uniform(0.01, upper_lim_half_width), 3)
            alpha = round(random.uniform(-180, 180), 2)
            beta = round(random.uniform(-180, 180), 2)
            gamma = round(random.uniform(-180, 180), 2)

        obj_1 += ' {c0} {x0} {y0} {x0} {a} {b} {c} {alpha} {beta} {gamma}'.format(c0=m0,
                                                                                  # round(random.uniform(0.1,1),2),
                                                                                  x0=x0,
                                                                                  y0=y0,
                                                                                  z0=z0,
                                                                                  a=a,
                                                                                  b=b,
                                                                                  c=c,
                                                                                  alpha=alpha,
                                                                                  beta=beta,
                                                                                  gamma=gamma)

        obj_2 += ' {c0} {x0} {y0} {x0} {a} {b} {c} {alpha} {beta} {gamma}'.format(c0=d0,
                                                                                  # round(random.uniform(0.1,1),2),
                                                                                  x0=x0,
                                                                                  y0=y0,
                                                                                  z0=z0,
                                                                                  a=a,
                                                                                  b=b,
                                                                                  c=c,
                                                                                  alpha=alpha,
                                                                                  beta=beta,
                                                                                  gamma=gamma)

        objects_1.append(obj_1)
        objects_2.append(obj_2)

    phantom_string_1 = '''#----------------------------------------------------
# random phantom number {num}
Model : {num};
Components : {comp};
TimeSteps : 1;
'''.format(num=n, comp=num_of_components)

    phantom_string_2 = '''#----------------------------------------------------
# random phantom number {num}
Model : {num};
Components : {comp};
TimeSteps : 1;
'''.format(num=n, comp=num_of_components)

    for obj in objects_1:
        phantom_string_1 += 'Object : ' + obj + '\n'

    with open(path_1, 'a') as file:
        file.write(phantom_string_1)

    for obj in objects_2:
        phantom_string_2 += 'Object : ' + obj + '\n'

    with open(path_2, 'a') as file:
        file.write(phantom_string_2)

Code that make projections of 3D randomly generated phantoms

import os
import numpy as np
import matplotlib.pyplot as plt
import tifffile
from sklearn import preprocessing
from tomophantom import TomoP3D

# Define the path and name of the `.dat` file that will contain the phantoms model description. If the file exists, the new models will be appended.python

path_1 = '/home/mli/tomograms/pycharm_demos/single_density/without_Gaussian_shapes/normal_test/without_noise/1024*1024_with_oversampling=2/Generate_Training_Data/Phantom3DLibrary_For_att_without_Gaussian_3.dat'
path_2 = '/home/mli/tomograms/pycharm_demos/single_density/without_Gaussian_shapes/normal_test/without_noise/1024*1024_with_oversampling=2/Generate_Training_Data/Phantom3DLibrary_For_phase_without_Gaussian_3.dat'

# Define the image size.

N_size = 1024

#Horiz_det = int(np.sqrt(2)*N_size) # detector column count (horizontal)
Horiz_det = N_size # detector column count (horizontal)
print("Horiz_det: ", Horiz_det)

Vert_det = N_size # detector row count (vertical) (no reason for it to be > N)
print("Vert_det: ", Vert_det)

# angles_num = int(0.5*np.pi*N_size); # angles number
angles_num = 5
print("angles_num: ", angles_num)

angles = np.linspace(0.0, 179.9, angles_num, dtype='float32') # in degrees
print("angles: ", angles)

# Created files
os.makedirs('../projection_2', exist_ok=True)
os.makedirs('../projection_2/tif_format', exist_ok=True)
os.makedirs('../projection_2/npy_format', exist_ok=True)

# Create the directories to where the files will be saved in tif format. For attenuation
os.makedirs('../projection_2/tif_format/attenuation_projection_in_tif_format', exist_ok=True)

# Create the directories to where the files will be saved in tif format. For phase
os.makedirs('../projection_2/tif_format/phase_projection_in_tif_format', exist_ok=True)

# Create the directories to where the files will be saved in npy format. For attenuation
os.makedirs('../projection_2/npy_format/attenuation_projection_in_npy_format', exist_ok=True)

# Create the directories to where the files will be saved in npy format. For phase
os.makedirs('../projection_2/npy_format/phase_projection_in_npy_format', exist_ok=True)

# Create the directories to where the files contain the stacks of attenuation and phase together and save in .npy format
os.makedirs('../projection_2/npy_format/projection_stack_in_npy_format', exist_ok=True)

# Generate the phantoms from models *a* to *b* from the *Phantom3DLibrary3.dat* library using TomoP3D.ModelSino.
pixel_size = 0.7e-6
scaling_factor = pixel_size * 100

# Projection
N = 1024
# For attenuation
projData3D_analyt_list_For_attenuation = []
for model in range(901, 1001):  # For phantom No.1 ~ 1000
    projData3D_analyt_attenuation = TomoP3D.ModelSino(model, N_size, Horiz_det, Vert_det, angles,
                                                      path_1) * scaling_factor
    print('projection shape: ', projData3D_analyt_attenuation.shape)
    projData3D_analyt_list_For_attenuation.append(projData3D_analyt_attenuation)

# For phase
projData3D_analyt_list_For_phase = []
for model in range(901, 1001):  # For phantom No.1 ~ 99
    projData3D_analyt_phase = TomoP3D.ModelSino(model, N_size, Horiz_det, Vert_det, angles, path_2) * scaling_factor
    projData3D_analyt_list_For_phase.append(projData3D_analyt_phase)

# At this step you may: extract one slice for each projection sinogram, extract one angle
# Save in tiff picure format

# For attenuation
for model in range(1, 101):      # For phantom No.1 ~ 1000
    projected_3D = projData3D_analyt_list_For_attenuation[model-1]
    projected_slice = projected_3D[::, 0, ::]
    tifffile.imsave('../projection_2/tif_format/attenuation_projection_in_tif_format/{:05d}.tif'.format(model+900), projected_slice.astype(np.float32))
    np.save('../projection_2/npy_format/attenuation_projection_in_npy_format/{:05d}.npy'.format(model+900), projected_slice.astype(np.float32))

# For phase
for model in range(1, 101):      # For phantom No.1 ~ 1000
    projected_3D = projData3D_analyt_list_For_phase[model-1]
    projected_slice = projected_3D[::, 0, ::]
    tifffile.imsave('../projection_2/tif_format/phase_projection_in_tif_format/{:05d}.tif'.format(model+900), projected_slice.astype(np.float32))
    np.save('../projection_2/npy_format/phase_projection_in_npy_format/{:05d}.npy'.format(model+900), projected_slice.astype(np.float32))

# At this step you may: extract one slice for each projection sinogram (both attenuation and phase), extract one angle
# Stack them together and save in .npy format

for model in range(1, 101):      # For phantom No.1 ~ 1000
    projected_slice = [0, 0]
    projected_3D_attenuation = projData3D_analyt_list_For_attenuation[model-1]
    projected_3D_phase = projData3D_analyt_list_For_phase[model-1]
    projected_slice[0] = projected_3D_attenuation[::, 0, ::]
    projected_slice[1] = projected_3D_phase[::, 0, ::]
    projected_slice = np.array(projected_slice)
    np.save('../projection_2/npy_format/projection_stack_in_npy_format/{:05d}.npy'.format(model+900), projected_slice.astype(np.float32))

[paskino]

I see what you are trying to do.

The python interface is not very nice, but I guess what you want to achieve I would do with TomoP3D.ObjectSino

Notice that the paramenter objlist is a list with dictionaries describing each object in the phantom along these lines:

obj1 = {}
obj1['C0'] = c0
obj1['x0'] = x0
...
obj2 = {} ...

objlist = [obj1, obj2,]

I also see that the angles phi2 and phi3 will not be used, probably a bug.

[lokeaichirou]

I see what you are trying to do.

The python interface is not very nice, but I guess what you want to achieve I would do with TomoP3D.ObjectSino

Notice that the paramenter objlist is a list with dictionaries describing each object in the phantom along these lines:

obj1 = {}
obj1['C0'] = c0
obj1['x0'] = x0
...
obj2 = {} ...

objlist = [obj1, obj2,]

I also see that the angles phi2 and phi3 will not be used, probably a bug.

Thanks for suggestion. I will do it. Since ph12 and phi3 will not be used (be zero in practice), it means the 3D projection images will still remain like the rectangular shapes without any tilting even I set phi2 and phi3 that makes phantom has a certain tilting angle? right?

[paskino]

Apparently yes. But it should be an easy fix.

[paskino]

I made a PR #88 which should fix the bug.

[paskino]

@lokeaichirou I just merged #88 which should fix the bug. It's being built now. In a few hours you should be able to find the updated version.

Please install as conda install -c ccpi/label/dev tomophantom to get it.

[dkazanc]

Hi @lokeaichirou , I removed two of Euler angles on purpose. The reason is that the 3D analytical projection data generator is not finished so providing these angles to the function will not give you the right answer. It is a limitation at the moment but I might find some time to complete it.

The 3D phantom, however can be generated by providing all 3 Euler angles. The objects can be controlled fully. So to avoid errors with someone using all 3 angles for the phantom and for the data I have removed them. So it is a hack rather than a bug, sorry. If you really need a full 3D control (tilting, rolling, panning), then it is possible to pass all angles while generating a phantom and then use a numerical forward projector (e.g. Astra toolbox) to generate data. thanks

[lokeaichirou]

Hi @lokeaichirou , I removed two of Euler angles on purpose. The reason is that the 3D analytical projection data generator is not finished so providing these angles to the function will not give you the right answer. It is a limitation at the moment but I might find some time to complete it.

The 3D phantom, however can be generated by providing all 3 Euler angles. The objects can be controlled fully. So to avoid errors with someone using all 3 angles for the phantom and for the data I have removed them. So it is a hack rather than a bug, sorry. If you really need a full 3D control (tilting, rolling, panning), then it is possible to pass all angles while generating a phantom and then use a numerical forward projector (e.g. Astra toolbox) to generate data. thanks Hi @dkazanc . Many thanks to your answer, it;s helpful. May I ask is it possible to use the tomophantom to generate random 3D phantoms first and then use astra to make their projections?because my project requires me to generate random 3D phantoms, and tomophantom is a tool can do that.

[dkazanc]

Hi @lokeaichirou, sure you can do it with this piece of code This is to generate a 'foam' type phantoms i.e. with spheres, ellipses and Gaussians. One probably can make a similar random generation with cubes but it is not implemented.

You probably can modify it in a way that you generate objects using all 3 angles and then use Astra to generate projections.