diff options
| -rw-r--r-- | Wrappers/Python/ccpi/reconstruction/geoms.py | 26 | ||||
| -rwxr-xr-x[-rw-r--r--] | Wrappers/Python/wip/simple_mc_demo.py (renamed from Wrappers/Python/test/simple_mc_demo.py) | 55 |
2 files changed, 50 insertions, 31 deletions
diff --git a/Wrappers/Python/ccpi/reconstruction/geoms.py b/Wrappers/Python/ccpi/reconstruction/geoms.py index 9f81fa0..a887c0c 100644 --- a/Wrappers/Python/ccpi/reconstruction/geoms.py +++ b/Wrappers/Python/ccpi/reconstruction/geoms.py @@ -2,12 +2,12 @@ class VolumeGeometry: def __init__(self, - voxel_num_x=None, - voxel_num_y=None, - voxel_num_z=None, - voxel_size_x=None, - voxel_size_y=None, - voxel_size_z=None, + voxel_num_x=0, + voxel_num_y=0, + voxel_num_z=0, + voxel_size_x=1, + voxel_size_y=1, + voxel_size_z=1, center_x=0, center_y=0, center_z=0, @@ -37,10 +37,16 @@ class VolumeGeometry: return self.center_y + 0.5*self.voxel_num_y*self.voxel_size_y def getMinZ(self): - return self.center_z - 0.5*self.voxel_num_z*self.voxel_size_z + if not voxel_num_z == 0: + return self.center_z - 0.5*self.voxel_num_z*self.voxel_size_z + else: + return 0 def getMaxZ(self): - return self.center_z + 0.5*self.voxel_num_z*self.voxel_size_z + if not voxel_num_z == 0: + return self.center_z + 0.5*self.voxel_num_z*self.voxel_size_z + else: + return 0 class SinogramGeometry: @@ -49,9 +55,9 @@ class SinogramGeometry: geom_type, dimension, angles, - pixel_num_h=None, + pixel_num_h=0, pixel_size_h=1, - pixel_num_v=None, + pixel_num_v=0, pixel_size_v=1, dist_source_center=None, dist_center_detector=None, diff --git a/Wrappers/Python/test/simple_mc_demo.py b/Wrappers/Python/wip/simple_mc_demo.py index a6959f9..54c9134 100644..100755 --- a/Wrappers/Python/test/simple_mc_demo.py +++ b/Wrappers/Python/wip/simple_mc_demo.py @@ -11,44 +11,57 @@ from ccpi.reconstruction.ops import * from ccpi.reconstruction.astra_ops import * from ccpi.reconstruction.geoms import * -import numpy as np +import numpy import matplotlib.pyplot as plt -test_case = 2 # 1=parallel2D, 2=cone2D +test_case = 1 # 1=parallel2D, 2=cone2D # Set up phantom N = 128 - +M = 100 numchannels = 3 -x = np.zeros((N,N,1,numchannels)) +vg = VolumeGeometry(voxel_num_x=N,voxel_num_y=N,channels=numchannels) +Phantom = VolumeData(geometry=vg) + +x = Phantom.as_array() +x[0 , round(N/4):round(3*N/4) , round(N/4):round(3*N/4) ] = 1.0 +x[0 , round(N/8):round(7*N/8) , round(3*N/8):round(5*N/8)] = 2.0 + +x[1 , round(N/4):round(3*N/4) , round(N/4):round(3*N/4) ] = 0.7 +x[1 , round(N/8):round(7*N/8) , round(3*N/8):round(5*N/8)] = 1.2 + +x[2 , round(N/4):round(3*N/4) , round(N/4):round(3*N/4) ] = 1.5 +x[2 , round(N/8):round(7*N/8) , round(3*N/8):round(5*N/8)] = 2.2 + +#x = numpy.zeros((N,N,1,numchannels)) -x[round(N/4):round(3*N/4),round(N/4):round(3*N/4),:,0] = 1.0 -x[round(N/8):round(7*N/8),round(3*N/8):round(5*N/8),:,0] = 2.0 +#x[round(N/4):round(3*N/4),round(N/4):round(3*N/4),:,0] = 1.0 +#x[round(N/8):round(7*N/8),round(3*N/8):round(5*N/8),:,0] = 2.0 -x[round(N/4):round(3*N/4),round(N/4):round(3*N/4),:,1] = 0.7 -x[round(N/8):round(7*N/8),round(3*N/8):round(5*N/8),:,1] = 1.2 +#x[round(N/4):round(3*N/4),round(N/4):round(3*N/4),:,1] = 0.7 +#x[round(N/8):round(7*N/8),round(3*N/8):round(5*N/8),:,1] = 1.2 -x[round(N/4):round(3*N/4),round(N/4):round(3*N/4),:,2] = 1.5 -x[round(N/8):round(7*N/8),round(3*N/8):round(5*N/8),:,2] = 2.2 +#x[round(N/4):round(3*N/4),round(N/4):round(3*N/4),:,2] = 1.5 +#x[round(N/8):round(7*N/8),round(3*N/8):round(5*N/8),:,2] = 2.2 f, axarr = plt.subplots(1,numchannels) for k in numpy.arange(3): - axarr[k].imshow(x[:,:,0,k],vmin=0,vmax=2.5) + axarr[k].imshow(x[k],vmin=0,vmax=2.5) plt.show() -vg = VolumeGeometry(N,N,None, 1,1,None,channels=numchannels) +#vg = VolumeGeometry(N,N,None, 1,1,None,channels=numchannels) -Phantom = VolumeData(x,geometry=vg) +#Phantom = VolumeData(x,geometry=vg) # Set up measurement geometry angles_num = 20; # angles number if test_case==1: - angles = np.linspace(0,np.pi,angles_num,endpoint=False) + angles = numpy.linspace(0,numpy.pi,angles_num,endpoint=False) elif test_case==2: - angles = np.linspace(0,2*np.pi,angles_num,endpoint=False) + angles = numpy.linspace(0,2*numpy.pi,angles_num,endpoint=False) else: NotImplemented @@ -76,7 +89,7 @@ elif test_case==2: channels=numchannels) # ASTRA operator using volume and sinogram geometries -Aop = AstraProjectorMC(vg, pg, 'gpu') +Aop = AstraProjectorMC(vg, pg, 'cpu') # Try forward and backprojection @@ -84,21 +97,21 @@ b = Aop.direct(Phantom) fb, axarrb = plt.subplots(1,numchannels) for k in numpy.arange(3): - axarrb[k].imshow(b.as_array()[:,:,0,k],vmin=0,vmax=250) + axarrb[k].imshow(b.as_array()[k],vmin=0,vmax=250) plt.show() out2 = Aop.adjoint(b) fo, axarro = plt.subplots(1,numchannels) for k in range(3): - axarro[k].imshow(out2.as_array()[:,:,0,k],vmin=0,vmax=3500) + axarro[k].imshow(out2.as_array()[k],vmin=0,vmax=3500) plt.show() # Create least squares object instance with projector and data. f = Norm2sq(Aop,b,c=0.5) # Initial guess -x_init = VolumeData(np.zeros(x.shape),geometry=vg) +x_init = VolumeData(numpy.zeros(x.shape),geometry=vg) # FISTA options opt = {'tol': 1e-4, 'iter': 200} @@ -109,7 +122,7 @@ x_fista0, it0, timing0, criter0 = FISTA(x_init, f, None, opt) ff0, axarrf0 = plt.subplots(1,numchannels) for k in numpy.arange(3): - axarrf0[k].imshow(x_fista0.as_array()[:,:,0,k],vmin=0,vmax=2.5) + axarrf0[k].imshow(x_fista0.as_array()[k],vmin=0,vmax=2.5) plt.show() plt.semilogy(criter0) @@ -126,7 +139,7 @@ x_fista1, it1, timing1, criter1 = FISTA(x_init, f, g0, opt) ff1, axarrf1 = plt.subplots(1,numchannels) for k in numpy.arange(3): - axarrf1[k].imshow(x_fista1.as_array()[:,:,0,k],vmin=0,vmax=2.5) + axarrf1[k].imshow(x_fista1.as_array()[k],vmin=0,vmax=2.5) plt.show() plt.semilogy(criter1) |