dkazanc / tomophantom Goto Github PK

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

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

import tomophantom
import os
path = os.path.dirname(tomophantom.__file__)

import tomophantom
import tomophantom.supp.artifacts
from tomophantom.supp.artifacts import ArtifactClass

artifact = ArtifactClass(sino)

from tomophantom.supp.artifacts import ArtifactClass

artifact = ArtifactClass(sino)