From 1355d4be61c26230fb86abb03fadf1e2e6dfe5d9 Mon Sep 17 00:00:00 2001 From: jakobsj Date: Tue, 6 Mar 2018 11:16:32 +0000 Subject: Geometry (#31) * started geometry * Simple geometry implemented #2 * sino geometry simplified to standard * Added get min/max methods vol geom * Basic vol and sino geoms working in astra_op and simple_demo * Fixed style * Added center/offset to vol geom. --- Wrappers/Python/ccpi/reconstruction/astra_ops.py | 63 ++++++++++++++++ Wrappers/Python/ccpi/reconstruction/geoms.py | 96 ++++++++++++++++++++++++ Wrappers/Python/test/simple_demo.py | 12 ++- 3 files changed, 168 insertions(+), 3 deletions(-) create mode 100644 Wrappers/Python/ccpi/reconstruction/geoms.py diff --git a/Wrappers/Python/ccpi/reconstruction/astra_ops.py b/Wrappers/Python/ccpi/reconstruction/astra_ops.py index 0c95b73..bc6f5e1 100644 --- a/Wrappers/Python/ccpi/reconstruction/astra_ops.py +++ b/Wrappers/Python/ccpi/reconstruction/astra_ops.py @@ -21,6 +21,69 @@ import numpy from ccpi.framework import SinogramData, VolumeData from ccpi.reconstruction.ops import PowerMethodNonsquare +class AstraProjectorSimple(Operator): + """ASTRA projector modified to use DataSet and geometry.""" + def __init__(self, geomv, geomp, device): + super(AstraProjectorSimple, self).__init__() + + # ASTRA Volume geometry + self.vol_geom = astra.create_vol_geom(geomv.voxel_num_x, \ + geomv.voxel_num_y, \ + geomv.getMinX(), \ + geomv.getMaxX(), \ + geomv.getMinY(), \ + geomv.getMaxY()) + + # ASTRA Projections geometry + if geomp.dimension == '2D': + if geomp.geom_type == 'parallel': + self.proj_geom = astra.create_proj_geom('parallel', \ + geomp.pixel_size_h, \ + geomp.pixel_num_h, \ + geomp.angles) + elif geomp.geom_type == 'cone': + NotImplemented + else: + NotImplemented + + # ASTRA projector + if device == 'cpu': + # Note that 'line' is only for parallel (2D) and only one option + self.proj_id = astra.create_projector('line', self.proj_geom, + self.vol_geom) # for CPU + elif device == 'gpu': + self.proj_id = astra.create_projector('cuda', self.proj_geom, + self.vol_geom) # for GPU + else: + NotImplemented + + self.s1 = None + + def direct(self, IM): + + sinogram_id, DATA = astra.create_sino(IM.as_array(), self.proj_id) + astra.data2d.delete(sinogram_id) + return SinogramData(DATA) + + def adjoint(self, DATA): + rec_id, IM = astra.create_backprojection(DATA.as_array(), self.proj_id) + astra.data2d.delete(rec_id) + return VolumeData(IM) + + def delete(self): + astra.data2d.delete(self.proj_id) + + def get_max_sing_val(self): + self.s1, sall, svec = PowerMethodNonsquare(self,10) + return self.s1 + + def size(self): + return ( (self.proj_geom['ProjectionAngles'].size, \ + self.proj_geom['DetectorCount']), \ + (self.vol_geom['GridColCount'], \ + self.vol_geom['GridRowCount']) ) + + class AstraProjector(Operator): """A simple 2D/3D parallel/fan beam projection/backprojection class based on ASTRA toolbox""" diff --git a/Wrappers/Python/ccpi/reconstruction/geoms.py b/Wrappers/Python/ccpi/reconstruction/geoms.py new file mode 100644 index 0000000..f7f3329 --- /dev/null +++ b/Wrappers/Python/ccpi/reconstruction/geoms.py @@ -0,0 +1,96 @@ + +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, \ + center_x=0, \ + center_y=0, \ + center_z=0): + + self.voxel_num_x = voxel_num_x + self.voxel_num_y = voxel_num_y + self.voxel_num_z = voxel_num_z + self.voxel_size_x = voxel_size_x + self.voxel_size_y = voxel_size_y + self.voxel_size_z = voxel_size_z + self.center_x = center_x + self.center_y = center_y + self.center_z = center_z + + def getMinX(self): + return self.center_x - 0.5*self.voxel_num_x*self.voxel_size_x + + def getMaxX(self): + return self.center_x + 0.5*self.voxel_num_x*self.voxel_size_x + + def getMinY(self): + return self.center_y - 0.5*self.voxel_num_y*self.voxel_size_y + + def getMaxY(self): + 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 + + def getMaxZ(self): + return self.center_z + 0.5*self.voxel_num_z*self.voxel_size_z + + +class SinogramGeometry: + + def __init__(self, \ + geom_type, \ + dimension, \ + angles, \ + pixel_num_h=None, \ + pixel_size_h=None, \ + pixel_num_v=None, \ + pixel_size_v=None, \ + dist_source_center=None, \ + dist_center_detector=None, \ + ): + """ + General inputs for standard type projection geometries + detectorDomain or detectorpixelSize: + If 2D + If scalar: Width of detector or single detector pixel + If 2-vec: Error + If 3D + If scalar: Width in both dimensions + If 2-vec: Vertical then horizontal size + grid + If 2D + If scalar: number of detectors + If 2-vec: error + If 3D + If scalar: Square grid that size + If 2-vec vertical then horizontal size + cone or parallel + 2D or 3D + parallel_parameters: ? + cone_parameters: + source_to_center_dist (if parallel: NaN) + center_to_detector_dist (if parallel: NaN) + standard or nonstandard (vec) geometry + angles + angles_format radians or degrees + """ + self.geom_type = geom_type # 'parallel' or 'cone' + self.dimension = dimension # 2D or 3D + self.angles = angles + + self.dist_source_center = dist_source_center + self.dist_center_detector = dist_center_detector + + self.pixel_num_h = pixel_num_h + self.pixel_size_h = pixel_size_h + self.pixel_num_v = pixel_num_v + self.pixel_size_v = pixel_size_v + + + diff --git a/Wrappers/Python/test/simple_demo.py b/Wrappers/Python/test/simple_demo.py index 4b7e89e..7a6b049 100644 --- a/Wrappers/Python/test/simple_demo.py +++ b/Wrappers/Python/test/simple_demo.py @@ -8,6 +8,7 @@ from ccpi.reconstruction.algs import * from ccpi.reconstruction.funcs import * from ccpi.reconstruction.ops import * from ccpi.reconstruction.astra_ops import * +from ccpi.reconstruction.geoms import * import numpy as np import matplotlib.pyplot as plt @@ -22,10 +23,10 @@ x[round(N/8):round(7*N/8),round(3*N/8):round(5*N/8)] = 2.0 plt.imshow(x) plt.show() -#vg = VolumeGeometry(grid=(N,N), domain=((-N/2,N/2),(-N/2,N/2))) +vg = VolumeGeometry(N,N,None, 1,1,None) -#Phantom = VolumeData(x,geometry=vg) Phantom = VolumeData(x) +#Phantom = VolumeData(x) # Set up measurement geometry angles_num = 20; # angles number @@ -37,14 +38,19 @@ det_num = N SourceOrig = 500 OrigDetec = 0 +# Parallelbeam geometry test +pg = SinogramGeometry('parallel','2D',angles,det_num,det_w) + # Set up ASTRA projector #Aop = AstraProjector(vg, angles, N,'gpu') #Aop = AstraProjector(det_w, det_num, angles, projtype='parallel',device='gpu') -Aop = AstraProjector(det_w, det_num, SourceOrig, +Aop_old = AstraProjector(det_w, det_num, SourceOrig, OrigDetec, angles, N,'fanbeam','gpu') +Aop = AstraProjectorSimple(vg, pg, 'gpu') + # Try forward and backprojection b = Aop.direct(Phantom) -- cgit v1.2.3