Skip to content

Commit

Permalink
Merge pull request #14 from dkazanc/typecheck
Browse files Browse the repository at this point in the history
new interface to create Objects and Sinograms.
  • Loading branch information
dkazanc committed Apr 1, 2018
2 parents 4f03d08 + efa3fd0 commit fae3849
Show file tree
Hide file tree
Showing 20 changed files with 1,576 additions and 1,312 deletions.
2 changes: 1 addition & 1 deletion functions/TomoP2DModel.c
Expand Up @@ -301,4 +301,4 @@ void mexFunction(
printf("%s %i %s \n", "Model no. ", ModelSelected, "is not found!");
mexErrMsgTxt("No object found, check models file");
}
}
}
6 changes: 3 additions & 3 deletions functions/TomoP2DModelSino.c
Expand Up @@ -173,13 +173,13 @@ void mexFunction(
printf("%s %f\n", "b (object size) must be positive in [0,2] range, the given value is", b);
mexErrMsgTxt("b (object size) must be positive in [0,2] range");
break; }
printf("\nObject : %s \nC0 : %f \nx0 : %f \ny0 : %f \na : %f \nb : %f \n", tmpstr2, C0, x0, y0, a, b);
printf("\nObject : %s \nC0 : %f \nx0 : %f \ny0 : %f \na : %f \nb : %f \n", tmpstr2, C0, x0, y0, a, b);

TomoP2DObjectSino_core(A, N, P, Th, (int)NStructElems, CenTypeIn, tmpstr2, C0, x0, y0, b, a, -psi_gr1, 0); /* Matlab */
}
}
else {
/**************************************************/
/**************************************************/
printf("\n %s %i %s \n", "Temporal 2D+time model", ModelSelected, " is selected");
/* temporal phantom 2D + time (3D) */
const mwSize N_dims[3] = {NStructElems, P, steps}; /*format: X-detectors, Y-angles dim, Time-Frames*/
Expand Down Expand Up @@ -288,7 +288,7 @@ void mexFunction(
/*loop over time frames*/
for(tt=0; tt < steps; tt++) {

TomoP2DObjectSino_core(A, N, P, Th, (int)NStructElems, CenTypeIn, tmpstr2, C_t, -y_t, x_t, b_t, a_t, -phi_t, tt);
TomoP2DObjectSino_core(A, N, P, Th, (int)NStructElems, CenTypeIn, tmpstr2, C_t, x_t, y_t, b_t, a_t, -phi_t, tt);

/* calculating new coordinates of an object */
if (distance != 0.0f) {
Expand Down
951 changes: 475 additions & 476 deletions functions/TomoP2DModelSino_core.c

Large diffs are not rendered by default.

657 changes: 330 additions & 327 deletions functions/TomoP2DModel_core.c

Large diffs are not rendered by default.

2 changes: 1 addition & 1 deletion functions/TomoP3DModel_core.c
Expand Up @@ -148,7 +148,7 @@ float TomoP3DObject_core(float *A, int N, char *Object,
}}}
}
if (strcmp("cuboid",Object) == 0) {
/* the object is a cube */
/* the object is a cuboid */
float x0r, y0r, HX, HY;
a2 = 0.5f*a;
b2 = 0.5f*b;
Expand Down
48 changes: 24 additions & 24 deletions functions/models/Phantom2DLibrary.dat
Expand Up @@ -148,30 +148,30 @@ Model : 11;
Components : 25;
TimeSteps : 1;
Object : ellipse .3e00 0e00 0.0e00 .9500 .95e00 00.e00;
Object : ellipse .5e00 0.0e00 0.7e00 0.2e00 0.2e00 00.e00;
Object : ellipse .500e00 -.6e00 .35e00 0.2e00 0.2e00 00.e00;
Object : ellipse .500e00 -.6e00 -.35e00 0.2e00 0.2e00 00.e00;
Object : ellipse .500e00 0.0e00 -.7e00 0.2e00 0.2e00 00.e00;
Object : ellipse .500e00 .6e00 -.35e00 0.2e00 0.2e00 00.e00;
Object : ellipse .500e00 .6e00 .35e00 0.2e00 0.2e00 00.e00;
Object : ellipse .500e00 -.2e00 .35e00 0.1e00 0.1e00 00.e00;
Object : ellipse .500e00 -.4e00 .0e00 0.1e00 0.1e00 00.e00;
Object : ellipse .500e00 -.2e00 -.35e00 0.1e00 0.1e00 00.e00;
Object : ellipse .500e00 .2e00 -.35e00 0.1e00 0.1e00 00.e00;
Object : ellipse .500e00 .4e00 .0e00 0.1e00 0.1e00 00.e00;
Object : ellipse .500e00 0.2e00 .35e00 0.1e00 0.1e00 00.e00;
Object : rectangle .500e00 0.0e00 0.085e00 0.1e00 0.1e00 00.e00;
Object : rectangle .500e00 0.085e00 0.00e00 0.1e00 0.1e00 00.e00;
Object : rectangle .500e00 -0.085e00 0.0e00 0.1e00 0.1e00 00.e00;
Object : rectangle .500e00 0.0e00 -0.085e00 0.1e00 0.1e00 00.e00;
Object : ellipse .5e00 0.0e00 0.0e00 0.02e00 0.11e00 45.e00;
Object : ellipse .5e00 0.0e00 0.0e00 0.02e00 0.11e00 -45.e00;
Object : ellipse .800e00 .7e00 .0e00 0.01e00 0.01e00 00.e00;
Object : ellipse .80e00 -.38e00 -.65e00 0.01e00 0.01e00 00.e00;
Object : ellipse .90e00 -.38e00 .65e00 0.01e00 0.01e00 00.e00;
Object : ellipse .90e00 .0e00 -.35e00 0.01e00 0.01e00 00.e00;
Object : ellipse .900e00 .32e00 .18e00 0.01e00 0.01e00 00.e00;
Object : ellipse .800e00 -.32e00 .18e00 0.01e00 0.01e00 00.e00;
Object : ellipse .5e00 0.7 0.0 0.2e00 0.2e00 00.e00;
Object : ellipse .500e00 0.35 -0.6 0.2e00 0.2e00 00.e00;
Object : ellipse .500e00 -0.35 -0.6 0.2e00 0.2e00 00.e00;
Object : ellipse .500e00 -0.7 0.0 0.2e00 0.2e00 00.e00;
Object : ellipse .500e00 -.35e00 .6e00 0.2e00 0.2e00 00.e00;
Object : ellipse .500e00 .35e00 .6e00 0.2e00 0.2e00 00.e00;
Object : ellipse .500e00 .35e00 -.2e00 0.1e00 0.1e00 00.e00;
Object : ellipse .500e00 .0e00 -.4e00 0.1e00 0.1e00 00.e00;
Object : ellipse .500e00 -.35e00 -.2e00 0.1e00 0.1e00 00.e00;
Object : ellipse .500e00 -.35e00 .2e00 0.1e00 0.1e00 00.e00;
Object : ellipse .500e00 .0e00 .4e00 0.1e00 0.1e00 00.e00;
Object : ellipse .500e00 .35e00 0.2 0.1e00 0.1e00 00.e00;
Object : rectangle .500e00 0.085e00 0.0 0.1e00 0.1e00 00.e00;
Object : rectangle .500e00 0.0 0.085 0.1e00 0.1e00 00.e00;
Object : rectangle .500e00 0.0 -0.085 0.1e00 0.1e00 00.e00;
Object : rectangle .500e00 -0.085 0.0 0.1e00 0.1e00 00.e00;
Object : ellipse .5e00 0.0 0.0 0.02e00 0.11e00 45.e00;
Object : ellipse .5e00 0.0 0.0 0.02e00 0.11e00 -45.e00;
Object : ellipse .800e00 0.0 .7e00 0.01e00 0.01e00 00.e00;
Object : ellipse .80e00 -.65e00 -.38e00 0.01e00 0.01e00 00.e00;
Object : ellipse .90e00 .65e00 -.38e00 0.01e00 0.01e00 00.e00;
Object : ellipse .90e00 -.35e00 .0e00 0.01e00 0.01e00 00.e00;
Object : ellipse .900e00 .18e00 .32e00 0.01e00 0.01e00 00.e00;
Object : ellipse .800e00 .18e00 -.32e00 0.01e00 0.01e00 00.e00;
#----------------------------------------------------
# composite - rectangles and ellipses
Model : 12;
Expand Down
6 changes: 3 additions & 3 deletions functions/models/Phantom3DLibrary.dat
Expand Up @@ -296,12 +296,12 @@ Object : gaussian 1.00 0.55 0.55 0.0 0.25 0.25 0.25 0.0 0.0 0.0;
Object : gaussian 1.00 0.55 -0.55 0.0 0.25 0.25 0.25 0.0 0.0 0.0;
Object : gaussian 1.00 -0.55 -0.55 0.0 0.25 0.25 0.25 0.0 0.0 0.0;
#----------------------------------------------------
# 1 volumetric Gaussian + 1 paraboloid
# 1 gaussian + 1 cuboid
Model : 10;
Components : 2;
TimeSteps : 1;
Object : gaussian 1. -0.3 0.1 0. 0.3 0.2 0.5 30. 0. 0.;
Object : paraboloid 1. 0.4 -0.1 0. 0.2 0.2 0.3 0. 0. 0.;
Object : gaussian 1.00 -0.25 -0.15 0.0 0.3 0.2 0.3 35.0 0.0 0.0;
Object : cuboid 1.00 0.1 0.2 0.0 0.15 0.35 0.6 -60.0 0.0 0.0;
#----------------------------------------------------
#****************************************************
#---------TEMPORAL (4D) Phantoms [no . >=100]--------
Expand Down
2 changes: 1 addition & 1 deletion matlab/Model2DGeneratorDemo.m
Expand Up @@ -11,7 +11,7 @@
% adding paths
addpath('../functions/models/'); addpath('compiled/'); addpath('supplem/');

ModelNo = 12; % Select a model from Phantom2DLibrary.dat
ModelNo = 4; % Select a model from Phantom2DLibrary.dat
% Define phantom dimensions
N = 512; % x-y size (squared image)

Expand Down
2 changes: 1 addition & 1 deletion matlab/Model3DGeneratorDemo.m
Expand Up @@ -9,7 +9,7 @@
% adding paths
addpath('../functions/models/'); addpath('compiled/'); addpath('supplem/');

ModelNo = 1; % Select a model
ModelNo = 10; % Select a model
% Define phantom dimensions
N = 256; % x-y-z size (cubic image)

Expand Down
2 changes: 1 addition & 1 deletion matlab/Object3DGeneratorDemo.m
Expand Up @@ -9,7 +9,7 @@
addpath('../functions/models/'); addpath('compiled/'); addpath('supplem/');

% Define phantom dimensions
N = 1024; % x-y-z size (cubic image)
N = 256; % x-y-z size (cubic image)

% define parameters
paramsObject.Ob = 'gaussian';
Expand Down
3 changes: 2 additions & 1 deletion matlab/supplem/rec2Dastra.m
Expand Up @@ -7,7 +7,6 @@
cfg = astra_struct('FBP_CUDA');
cfg.FilterType = 'Ram-Lak';


if (strcmp(device, 'cpu') == 1)
proj_id = astra_create_projector('strip', proj_geom, vol_geom);
cfg.ProjectorId = proj_id;
Expand All @@ -19,6 +18,7 @@

rec_id = astra_mex_data2d('create', '-vol', vol_geom);
sinogram_id = astra_mex_data2d('create', '-sino', proj_geom, sino);

cfg.ReconstructionDataId = rec_id;
cfg.ProjectionDataId = sinogram_id;

Expand All @@ -33,6 +33,7 @@
astra_mex_algorithm('delete', alg_id);
astra_mex_data2d('delete', sinogram_id);
astra_mex_data2d('delete', rec_id);

if (strcmp(device, 'cpu') == 1)
astra_mex_data2d('delete', proj_id);
end
10 changes: 7 additions & 3 deletions python/Demos/DemoModel.py
Expand Up @@ -18,10 +18,11 @@
from tomophantom import TomoP2D
from astraOP import AstraTools
#%%
model = 4
model = 1
N_size = 512
#specify a full path to the parameters file
pathTP = '../../functions/models/Phantom2DLibrary.dat'
#objlist = modelfile2Dtolist(pathTP, model) # one can extract parameters
#This will generate a N_size x N_size phantom (2D)
phantom_2D = TomoP2D.Model(model, N_size, pathTP)

Expand All @@ -46,7 +47,7 @@
plt.colorbar(ticks=[0, 150, 250], orientation='vertical')
plt.title('{}''{}'.format('Analytical sinogram of model no.',model))
#%%
Atools = AstraTools(P, angles_rad - 0.5*np.pi, N_size, 'cpu') # initiate a class object
Atools = AstraTools(P, angles_rad, N_size, 'cpu') # initiate a class object
sino_num_ASTRA = Atools.forwproj(phantom_2D) # generate numerical sino (Ax)
#x = Atools.backproj(sino_an) # generate backprojection (A'b)

Expand All @@ -58,10 +59,11 @@
plt.imshow(sino_num_ASTRA,cmap="BuPu")
plt.title('Numerical sinogram')
plt.show()
#calculate norm
rmse1 = np.linalg.norm(sino_an - sino_num_ASTRA)/np.linalg.norm(sino_num_ASTRA)
#%%
print ("Reconstructing analytical sinogram using FBP (astra TB)...")
FBPrec1 = Atools.fbp2D(sino_an)

plt.figure(3)
plt.imshow(FBPrec1, vmin=0, vmax=1, cmap="BuPu")
plt.colorbar(ticks=[0, 0.5, 1], orientation='vertical')
Expand All @@ -79,6 +81,7 @@
plt.imshow(abs(FBPrec1-FBPrec2), vmin=0, vmax=0.05, cmap="BuPu")
plt.colorbar(ticks=[0, 0.02, 0.05], orientation='vertical')
plt.title('FBP rec differences')
rmse2 = np.linalg.norm(FBPrec1 - FBPrec2)/np.linalg.norm(FBPrec2)
#%%
"""
print ("Reconstructing using SIRT...")
Expand Down Expand Up @@ -122,3 +125,4 @@
plt.title('3D Phantom, sagittal view')
plt.show()
"""
#%%
87 changes: 87 additions & 0 deletions python/Demos/DemoModel2.py
@@ -0,0 +1,87 @@
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
GPLv3 license (ASTRA toolbox)
Note that the TomoPhantom package is released under Apache License, Version 2.0
Script to generate 2D/3D analytical phantoms and their sinograms
If one needs to modify/add phantoms, please edit Phantom2DLibrary.dat or
Phantom3DLibrary.dat
>>>> Prerequisites: ASTRA toolbox, if one needs to do reconstruction <<<<<
Main difference from DemoModel.py is that we extract all parameters from the
library file using Python and then pass it to the Object function instead of model
Run demo from the folder "Demos"
@author: Daniil Kazantsev
"""
import numpy as np
import matplotlib.pyplot as plt
from tomophantom import TomoP2D
from astraOP import AstraTools
from libraryToDict import modelfile2Dtolist

#%%
model = 11
N_size = 512
#specify a full path to the parameters file
pathTP = '../../functions/models/Phantom2DLibrary.dat'
objlist = modelfile2Dtolist(pathTP, model) # extract parameters using Python
#This will generate a N_size x N_size phantom (2D)
phantom_2D = TomoP2D.Object(N_size, objlist)

plt.close('all')
plt.figure(1)
plt.rcParams.update({'font.size': 21})
plt.imshow(phantom_2D, vmin=0, vmax=1, cmap="BuPu")
plt.colorbar(ticks=[0, 0.5, 1], orientation='vertical')
plt.title('{}''{}'.format('2D Phantom using model no.',model))
#%%
# create sinogram analytically
angles_num = int(0.5*np.pi*N_size); # angles number
angles = np.linspace(0,180,angles_num,dtype='float32')
angles_rad = angles*(np.pi/180)
P = int(np.sqrt(2)*N_size) #detectors

sino_an = TomoP2D.ObjectSino(N_size, P, angles, objlist)

plt.figure(2)
plt.rcParams.update({'font.size': 21})
plt.imshow(sino_an, cmap="BuPu")
plt.colorbar(ticks=[0, 150, 250], orientation='vertical')
plt.title('{}''{}'.format('Analytical sinogram of model no.',model))
#%%
Atools = AstraTools(P, angles_rad, N_size, 'cpu') # initiate a class object
sino_num_ASTRA = Atools.forwproj(phantom_2D) # generate numerical sino (Ax)
#x = Atools.backproj(sino_an) # generate backprojection (A'b)

plt.figure(3)
plt.subplot(121)
plt.imshow(sino_an,cmap="BuPu")
plt.title('Analytical sinogram')
plt.subplot(122)
plt.imshow(sino_num_ASTRA,cmap="BuPu")
plt.title('Numerical sinogram')
plt.show()
#%%
print ("Reconstructing analytical sinogram using FBP (astra TB)...")
FBPrec1 = Atools.fbp2D(sino_an)
plt.figure(4)
plt.imshow(FBPrec1, vmin=0, vmax=1, cmap="BuPu")
plt.colorbar(ticks=[0, 0.5, 1], orientation='vertical')
plt.title('FBP Reconstructed Phantom (analyt)')
#%%
print ("Reconstructing numerical sinogram using FBP (astra TB)...")
FBPrec2 = Atools.fbp2D(sino_num_ASTRA)
plt.figure(5)
plt.imshow(FBPrec2, vmin=0, vmax=1, cmap="BuPu")
plt.colorbar(ticks=[0, 0.5, 1], orientation='vertical')
plt.title('FBP Reconstructed Phantom (numeric)')
#%%
plt.figure(6)
plt.imshow(abs(FBPrec1-FBPrec2), vmin=0, vmax=0.05, cmap="BuPu")
plt.colorbar(ticks=[0, 0.02, 0.05], orientation='vertical')
plt.title('FBP rec differences')
rmse2 = np.linalg.norm(FBPrec1 - FBPrec2)/np.linalg.norm(FBPrec2)
#%%
5 changes: 3 additions & 2 deletions python/Demos/DemoModel_temporal.py
Expand Up @@ -20,12 +20,13 @@
#%%
model = 102 # note that the selected model is temporal (2D + time)
N_size = 512
timeframes = 25
#specify a full path to the parameters file
pathTP = '../../functions/models/Phantom2DLibrary.dat'
#This will generate a N_size x N_size x Time frames phantom (2D + time)
phantom_2Dt = TomoP2D.ModelTemporal(model, N_size, pathTP)

timeframes = 25
plt.close('all')
plt.figure(1)
plt.rcParams.update({'font.size': 21})
plt.title('{}''{}'.format('2D+t phantom using model no.',model))
Expand Down Expand Up @@ -53,7 +54,7 @@
plt.draw
#%%
# reconstruct
Atools = AstraTools(P, angles_rad - 0.5*np.pi, N_size, 'cpu') # initiate a class object
Atools = AstraTools(P, angles_rad, N_size, 'cpu') # initiate a class object
FBPrec = Atools.fbp2D(sino[15,:,:].transpose())

plt.figure(3)
Expand Down

0 comments on commit fae3849

Please sign in to comment.