From 5e72985e2983e9be8117696fe0ec02388e7153f1 Mon Sep 17 00:00:00 2001
From: "Suren A. Chilingaryan" <csa@ipecompute4.ands.kit.edu>
Date: Wed, 27 Jul 2022 00:28:05 +0200
Subject: Latest stuff I found from Harish
---
dummy.py | 132 ++++++++++++++++++++++++++
src/__init__.py | 2 +
src/operators/ProjectionOperator.py | 67 +++++++++++++
src/operators/UfoStackedProjector.py | 60 ++++++++++++
src/operators/UfoStandardProjector.py | 58 +++++++++++
src/operators/__init__.py | 3 +
src/processors/FBP.py | 32 +++++++
src/processors/FBP_Stacked.py | 85 +++++++++++++++++
src/processors/FBP_Standard.py | 80 ++++++++++++++++
src/processors/UfoStackedBackProjector.py | 101 ++++++++++++++++++++
src/processors/UfoStackedForwardProjector.py | 100 +++++++++++++++++++
src/processors/UfoStandardBackProjector.py | 93 ++++++++++++++++++
src/processors/UfoStandardForwardProjector.py | 89 +++++++++++++++++
src/processors/__init__.py | 5 +
14 files changed, 907 insertions(+)
create mode 100644 dummy.py
create mode 100644 src/__init__.py
create mode 100644 src/operators/ProjectionOperator.py
create mode 100644 src/operators/UfoStackedProjector.py
create mode 100644 src/operators/UfoStandardProjector.py
create mode 100644 src/operators/__init__.py
create mode 100644 src/processors/FBP.py
create mode 100644 src/processors/FBP_Stacked.py
create mode 100644 src/processors/FBP_Standard.py
create mode 100644 src/processors/UfoStackedBackProjector.py
create mode 100644 src/processors/UfoStackedForwardProjector.py
create mode 100644 src/processors/UfoStandardBackProjector.py
create mode 100644 src/processors/UfoStandardForwardProjector.py
create mode 100644 src/processors/__init__.py
diff --git a/dummy.py b/dummy.py
new file mode 100644
index 0000000..047c021
--- /dev/null
+++ b/dummy.py
@@ -0,0 +1,132 @@
+# cil imports
+from cil.framework import ImageData, ImageGeometry
+from cil.framework import AcquisitionGeometry, AcquisitionData
+from cil.optimisation.algorithms import CGLS, SIRT
+
+import src.operators.ProjectionOperator as UfoProjectionOperator
+import cil.plugins.astra.operators.ProjectionOperator as AstraProjectionOperator
+import cil.plugins.astra.processors.FBP
+
+
+# External imports
+import tifffile
+import os
+import numpy as np
+import matplotlib.pyplot as plt
+import logging
+import time
+
+logging.basicConfig(level = logging.INFO)
+
+# set up default colour map for visualisation
+#cmap = "gray"
+
+# set the backend for FBP and the ProjectionOperator
+#device = 'gpu'
+
+def reconstruct():
+ # number of pixels
+ n_pixels = 4096
+ z_size = 1024
+ #slices=2
+ #rounding='int8'
+ #rounding='half'
+ #rounding='single'
+ #n_pixels = 1024
+
+ # Angles
+ angles = np.linspace(0, 180, n_pixels, endpoint=False, dtype=np.float32)
+
+
+ # Setup acquisition geometry
+ # with sufficient number of projections
+ ag = AcquisitionGeometry.create_Parallel3D()\
+ .set_angles(angles)\
+ .set_panel([n_pixels,z_size], pixel_size=1)
+ #ag = AcquisitionGeometry.create_Parallel2D()\
+ # .set_angles(angles)\
+ # .set_panel(n_pixels,pixel_size=1)
+
+ ag.set_labels(['vertical', 'angle', 'horizontal'])
+
+ # Setup image geometry
+ ig = ImageGeometry(voxel_num_x=n_pixels,
+ voxel_num_y=n_pixels,
+ voxel_num_z=z_size,
+ voxel_size_x=1,
+ voxel_size_y=1,
+ voxel_size_z=1)
+
+
+ # Get phantom
+ #home='/ccpi/data/new/phantom/'
+ home='/ccpi/data/new/1024/phantom/'
+ #home='/ccpi/data/test/1024/slices/'
+# phantom = []
+# for i in sorted(os.listdir(home)):
+# for i in range(0, z_size):
+ #file = os.path.join(home,i)
+ #img = tifffile.imread(file)
+# img = np.random.random_sample((n_pixels,n_pixels)).astype(np.float32)
+# if len(img.shape)==3:
+# img = np.squeeze(img, axis=0)
+# phantom.append(img)
+# phantom = np.asarray(phantom)
+ phantom = np.random.random_sample((z_size,n_pixels,n_pixels)).astype(np.float32)
+
+ print(phantom.shape)
+ phantom_reshaped = np.transpose(phantom, (0, 1, 2))
+ print(phantom_reshaped.shape)
+ phantom = ImageData(phantom_reshaped, deep_copy=False, geometry=ig)
+ phantom_arr = phantom.as_array()
+ tifffile.imsave('phantom_3d.tif', phantom_arr[:,:,0])
+
+ # Create projection operator using Astra-Toolbox.
+ out_file = "ufo_single.tif"
+ sino_file= "sino_single.tif"
+
+# Ufo-harish
+ A = UfoProjectionOperator(ig, ag, True, rounding, slices)
+# Ufo-farago
+# A = UfoProjectionOperator(ig, ag, False, rounding, slices)
+# Astra
+# A = AstraProjectionOperator(ig, ag, 'gpu')
+
+ #A.dot_test(A)
+
+ # Create an acqusition data (numerically)
+ sino = A.direct(phantom)
+ sino_arr = sino.as_array()
+ tifffile.imsave(sino_file,sino_arr[0,:,:])
+
+ # back_projection
+ back_projection = A.adjoint(sino)
+
+ # initial estimate - zero array in this case
+ initial = ig.allocate(0)
+
+ # setup CGLS
+ f = open("single.txt", "a")
+ for i in range(20,21):
+ total_time=0
+ for j in range(0,1):
+ cgls = CGLS(initial=initial,
+ operator=A,
+ data=sino,
+ tolerance=1e-16,
+ max_iteration = 50,
+ update_objective_interval = 1 )
+
+ # run N interations
+ start = time.time()
+ cgls.run(i, verbose=True)
+ end = time.time()
+ if(j>0):
+ total_time += (end-start)
+ f.write('Iteration {} Time {} \n'.format(i,str(total_time/3)))
+ f.close()
+ print("finished reconstruction in single precision")
+
+if __name__ == '__main__':
+ reconstruct()
+
diff --git a/src/__init__.py b/src/__init__.py
new file mode 100644
index 0000000..900d738
--- /dev/null
+++ b/src/__init__.py
@@ -0,0 +1,2 @@
+from .operators import ProjectionOperator
+from .processors import FBP
diff --git a/src/operators/ProjectionOperator.py b/src/operators/ProjectionOperator.py
new file mode 100644
index 0000000..76d5c0c
--- /dev/null
+++ b/src/operators/ProjectionOperator.py
@@ -0,0 +1,67 @@
+# ========================================================================
+# Copyright 2019 Science Technology Facilities Council
+# Copyright 2019 University of Manchester
+#
+# This work is part of the Core Imaging Library developed by Science Technology
+# Facilities Council and University of Manchester
+#
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+#
+# http://www.apache.org/licenses/LICENSE-2.0.txt
+#
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+#
+# =========================================================================
+
+from cil.framework import DataOrder
+from cil.optimisation.operators import LinearOperator, ChannelwiseOperator
+from src.operators.UfoStandardProjector import UfoStandardProjector
+from src.operators.UfoStackedProjector import UfoStackedProjector
+
+class ProjectionOperator(LinearOperator):
+ """Ufo projector modified to use DataSet and geometry."""
+
+ def __init__(self, geomv, geomp, stacked=True, precision_mode='single', stack_num=2):
+
+ super(ProjectionOperator, self).__init__(domain_geometry=geomv, range_geometry=geomp)
+
+ DataOrder.check_order_for_engine('astra', geomv)
+ DataOrder.check_order_for_engine('astra', geomp)
+
+ self.volume_geometry = geomv
+ self.sinogram_geometry = geomp
+
+ sinogram_geometry_sc = geomp.subset(channel=0)
+ volume_geometry_sc = geomv.subset(channel=0)
+
+ if stacked==True:
+ operator = UfoStackedProjector(volume_geometry_sc, sinogram_geometry_sc, precision_mode, stack_num)
+ else:
+ operator = UfoStandardProjector(volume_geometry_sc, sinogram_geometry_sc)
+
+ if geomp.channels > 1:
+ operator_full = ChannelwiseOperator(operator, self.sinogram_geometry.channels, dimension='prepend')
+ self.operator = operator_full
+ else:
+ self.operator = operator
+
+ def direct(self, IM, out=None):
+ return self.operator.direct(IM, out=out)
+
+ def adjoint(self, DATA, out=None):
+ return self.operator.adjoint(DATA, out=out)
+
+ def calculate_norm(self):
+ return self.operator.norm()
+
+ def domain_geometry(self):
+ return self.volume_geometry
+
+ def range_geometry(self):
+ return self.sinogram_geometry
diff --git a/src/operators/UfoStackedProjector.py b/src/operators/UfoStackedProjector.py
new file mode 100644
index 0000000..e2e65fd
--- /dev/null
+++ b/src/operators/UfoStackedProjector.py
@@ -0,0 +1,60 @@
+# ========================================================================
+# Copyright 2019 Science Technology Facilities Council
+# Copyright 2019 University of Manchester
+#
+# This work is part of the Core Imaging Library developed by Science Technology
+# Facilities Council and University of Manchester
+#
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+#
+# http://www.apache.org/licenses/LICENSE-2.0.txt
+#
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+#
+# =========================================================================
+
+from cil.optimisation.operators import LinearOperator
+from src.processors import UfoStackedForwardProjector, UfoStackedBackProjector
+
+class UfoStackedProjector(LinearOperator):
+ """UFO projector modified to use DataSet and geometry."""
+
+ def __init__(self, geomv, geomp, precision_mode='single', stack_num=2):
+
+ super(UfoStackedProjector, self).__init__(domain_geometry=geomv, range_geometry=geomp)
+
+ self.sinogram_geometry = geomp
+ self.volume_geometry = geomv
+
+ self.fp = UfoStackedForwardProjector(volume_geometry=geomv, sinogram_geometry=geomp,
+ precision_mode=precision_mode, stack_num=stack_num)
+ self.bp = UfoStackedBackProjector(volume_geometry=geomv, sinogram_geometry=geomp,
+ precision_mode=precision_mode, stack_num=stack_num)
+
+ def direct(self, IM, out=None):
+ self.fp.set_input(IM)
+
+ if out is None:
+ return self.fp.get_output()
+ else:
+ out.fill(self.fp.get_output())
+
+ def adjoint(self, DATA, out=None):
+ self.bp.set_input(DATA)
+
+ if out is None:
+ return self.bp.get_output()
+ else:
+ out.fill(self.bp.get_output())
+
+ def domain_geometry(self):
+ return self.volume_geometry
+
+ def range_geometry(self):
+ return self.sinogram_geometry
diff --git a/src/operators/UfoStandardProjector.py b/src/operators/UfoStandardProjector.py
new file mode 100644
index 0000000..cef1a36
--- /dev/null
+++ b/src/operators/UfoStandardProjector.py
@@ -0,0 +1,58 @@
+# ========================================================================
+# Copyright 2019 Science Technology Facilities Council
+# Copyright 2019 University of Manchester
+#
+# This work is part of the Core Imaging Library developed by Science Technology
+# Facilities Council and University of Manchester
+#
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+#
+# http://www.apache.org/licenses/LICENSE-2.0.txt
+#
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+#
+# =========================================================================
+
+from cil.optimisation.operators import LinearOperator
+from src.processors import UfoStandardForwardProjector, UfoStandardBackProjector
+
+class UfoStandardProjector(LinearOperator):
+ """UFO projector modified to use DataSet and geometry."""
+
+ def __init__(self, geomv, geomp):
+
+ super(UfoStandardProjector, self).__init__(domain_geometry=geomv, range_geometry=geomp)
+
+ self.sinogram_geometry = geomp
+ self.volume_geometry = geomv
+
+ self.fp = UfoStandardForwardProjector(volume_geometry=geomv, sinogram_geometry=geomp)
+ self.bp = UfoStandardBackProjector(volume_geometry=geomv, sinogram_geometry=geomp)
+
+ def direct(self, IM, out=None):
+ self.fp.set_input(IM)
+
+ if out is None:
+ return self.fp.get_output()
+ else:
+ out.fill(self.fp.get_output())
+
+ def adjoint(self, DATA, out=None):
+ self.bp.set_input(DATA)
+
+ if out is None:
+ return self.bp.get_output()
+ else:
+ out.fill(self.bp.get_output())
+
+ def domain_geometry(self):
+ return self.volume_geometry
+
+ def range_geometry(self):
+ return self.sinogram_geometry
diff --git a/src/operators/__init__.py b/src/operators/__init__.py
new file mode 100644
index 0000000..b5e6d8c
--- /dev/null
+++ b/src/operators/__init__.py
@@ -0,0 +1,3 @@
+from .ProjectionOperator import ProjectionOperator
+from .UfoStandardProjector import UfoStandardProjector
+from .UfoStackedProjector import UfoStackedProjector
diff --git a/src/processors/FBP.py b/src/processors/FBP.py
new file mode 100644
index 0000000..c48c590
--- /dev/null
+++ b/src/processors/FBP.py
@@ -0,0 +1,32 @@
+from cil.framework import DataProcessor
+from cil.framework import DataOrder
+from src.processors.FBP_Standard import FBP_Standard
+from src.processors.FBP_Stacked import FBP_Stacked
+
+class FBP(DataProcessor):
+ def __init__(self, volume_geometry, sinogram_geometry, stacked=False, precision_mode='single', stack_num=2):
+ if stacked==True:
+ processor = FBP_Stacked(volume_geometry, sinogram_geometry, precision_mode, stack_num)
+ else:
+ processor = FBP_Standard(volume_geometry, sinogram_geometry)
+
+ super(FBP, self).__init__(volume_geometry=volume_geometry, sinogram_geometry=sinogram_geometry,
+ stacked=stacked,
+ precision_mode=precision_mode, stack_num=stack_num, processor=processor)
+
+ self.processor = processor
+
+ def set_input(self, dataset):
+ return self.processor.set_input(dataset)
+
+ def get_input(self):
+ return self.processor.get_input()
+
+ def get_output(self, out=None):
+ return self.processor.get_output(out=None)
+
+ def check_input(self, dataset):
+ return self.processor.check_input(dataset)
+
+ def process(self, out=None):
+ return self.processor.process(out=None)
diff --git a/src/processors/FBP_Stacked.py b/src/processors/FBP_Stacked.py
new file mode 100644
index 0000000..b6310a1
--- /dev/null
+++ b/src/processors/FBP_Stacked.py
@@ -0,0 +1,85 @@
+from cil.framework import DataProcessor, ImageData
+import numpy as np
+
+import gi
+gi.require_version('Ufo','0.0')
+from gi.repository import Ufo
+
+class FBP_Stacked(DataProcessor):
+ def __init__(self, volume_geometry,
+ sinogram_geometry,
+ precision_mode='single',
+ stack_num=2):
+
+ kwargs = {
+ 'volume_geometry': volume_geometry,
+ 'sinogram_geometry': sinogram_geometry,
+ 'precision_mode': precision_mode,
+ 'stack_num': stack_num}
+
+ super(FBP_Stacked, self).__init__(**kwargs)
+
+ self.precision_mode = precision_mode
+ self.stack_num = stack_num
+ self.set_ImageGeometry(volume_geometry)
+ self.set_AcquisitionGeometry(sinogram_geometry)
+
+ def check_input(self, dataset):
+ if self.sinogram_geometry.dimension == '2D':
+ raise ValueError("Expected input dimensions is 3, got {0}" \
+ .format(dataset.number_of_dimensions))
+ return True
+
+ def set_ImageGeometry(self, volume_geometry):
+ self.volume_geometry = volume_geometry
+
+ def set_AcquisitionGeometry(self, sinogram_geometry):
+ self.sinogram_geometry = sinogram_geometry
+
+ def process(self, **kwargs):
+ # Get DATA
+ DATA = self.get_input()
+ DATA = DATA.as_array()
+
+ pm = Ufo.PluginManager()
+ graph = Ufo.TaskGraph()
+ scheduler = Ufo.Scheduler()
+
+ read = pm.get_task('memory-in')
+ read.props.pointer = DATA.__array_interface__['data'][0]
+ read.props.width = DATA.shape[1]
+ read.props.height = DATA.shape[0]
+ read.props.number = DATA.shape[2]
+ read.props.bitdepth = 32
+
+ fft = pm.get_task('fft')
+ fft.props.dimensions = 1
+
+ filter = pm.get_task('filter')
+
+ ifft = pm.get_task('ifft')
+ ifft.props.dimensions = 1
+
+ stack = pm.get_task('stack')
+ stack.props.number = self.stack_num
+
+ stacked_bp = pm.get_task('stacked-backproject')
+ stacked_bp.props.axis_pos = DATA.shape[1] / 2
+ stacked_bp.props.precision_mode = self.precision_mode
+
+ vol_arr = np.zeros(self.volume_geometry.shape, dtype=np.float32)
+
+ write = pm.get_task('memory-out')
+ write.props.pointer = vol_arr.__array_interface__['data'][0]
+ write.props.max_size = vol_arr.nbytes
+
+ graph.connect_nodes(read, fft)
+ graph.connect_nodes(fft, filter)
+ graph.connect_nodes(filter, ifft)
+ graph.connect_nodes(ifft, stack)
+ graph.connect_nodes(stack, stacked_bp)
+ graph.connect_nodes(stacked_bp, write)
+
+ scheduler.run(graph)
+ vol_arr = ImageData(vol_arr, deep_copy=False, geometry=self.volume_geometry.copy(), suppress_warning=True)
+ return vol_arr
diff --git a/src/processors/FBP_Standard.py b/src/processors/FBP_Standard.py
new file mode 100644
index 0000000..3fdb91c
--- /dev/null
+++ b/src/processors/FBP_Standard.py
@@ -0,0 +1,80 @@
+from cil.framework import DataProcessor, ImageData
+import numpy as np
+
+import gi
+gi.require_version('Ufo','0.0')
+from gi.repository import Ufo
+
+class FBP_Standard(DataProcessor):
+ def __init__(self, volume_geometry,
+ sinogram_geometry,
+ z_dim = None):
+
+ kwargs = {
+ 'volume_geometry': volume_geometry,
+ 'sinogram_geometry': sinogram_geometry,
+ 'z_dim': z_dim}
+
+ super(FBP_Standard, self).__init__(**kwargs)
+
+ self.z_dim = 0
+ self.set_ImageGeometry(volume_geometry)
+ self.set_AcquisitionGeometry(sinogram_geometry)
+
+ def check_input(self, dataset):
+ if self.sinogram_geometry.dimension == '2D':
+ self.z_dim = 1
+ elif self.sinogram_geometry.dimension == '3D':
+ self.z_dim = self.sinogram_geometry.shape[2]
+
+ return True
+
+ def set_ImageGeometry(self, volume_geometry):
+ self.volume_geometry = volume_geometry
+
+ def set_AcquisitionGeometry(self, sinogram_geometry):
+ self.sinogram_geometry = sinogram_geometry
+
+ def process(self, **kwargs):
+ # Get DATA
+ DATA = self.get_input()
+ DATA = DATA.as_array()
+
+ pm = Ufo.PluginManager()
+ graph = Ufo.TaskGraph()
+ scheduler = Ufo.Scheduler()
+
+ read = pm.get_task('memory-in')
+ read.props.pointer = DATA.__array_interface__['data'][0]
+ read.props.width = DATA.shape[1]
+ read.props.height = DATA.shape[0]
+ read.props.number = self.z_dim
+ read.props.bitdepth = 32
+
+ fft = pm.get_task('fft')
+ fft.props.dimensions = 1
+
+ filter = pm.get_task('filter')
+
+ ifft = pm.get_task('ifft')
+ ifft.props.dimensions = 1
+
+ bp = pm.get_task('backproject')
+ bp.props.axis_pos = DATA.shape[1] / 2
+
+ vol_arr = np.zeros(self.volume_geometry.shape, dtype=np.float32)
+
+ write = pm.get_task('memory-out')
+ write.props.pointer = vol_arr.__array_interface__['data'][0]
+ write.props.max_size = vol_arr.nbytes
+
+ graph.connect_nodes(read, fft)
+ graph.connect_nodes(fft, filter)
+ graph.connect_nodes(filter, ifft)
+ graph.connect_nodes(ifft, bp)
+ graph.connect_nodes(bp, write)
+
+ scheduler.run(graph)
+ vol_arr = ImageData(vol_arr, deep_copy=False, geometry=self.volume_geometry.copy(), suppress_warning=True)
+
+ return vol_arr
diff --git a/src/processors/UfoStackedBackProjector.py b/src/processors/UfoStackedBackProjector.py
new file mode 100644
index 0000000..a08b0c4
--- /dev/null
+++ b/src/processors/UfoStackedBackProjector.py
@@ -0,0 +1,101 @@
+import tifffile
+from cil.framework import DataProcessor, ImageData, DataOrder
+import gi
+
+gi.require_version('Ufo', '0.0')
+from gi.repository import Ufo
+import numpy as np
+
+
+class UfoStackedBackProjector(DataProcessor):
+
+ def __init__(self,
+ volume_geometry=None,
+ sinogram_geometry=None,
+ precision_mode=None,
+ stack_num=None):
+ kwargs = {
+ 'volume_geometry': volume_geometry,
+ 'sinogram_geometry': sinogram_geometry,
+ 'precision_mode': precision_mode,
+ 'stack_num': stack_num}
+
+ super(UfoStackedBackProjector, self).__init__(**kwargs)
+
+ self.precision_mode = precision_mode
+ self.stack_num = stack_num
+ self.set_ImageGeometry(volume_geometry)
+ self.set_AcquisitionGeometry(sinogram_geometry)
+
+ def check_input(self, dataset):
+
+ if self.sinogram_geometry.shape != dataset.geometry.shape:
+ raise ValueError("Dataset not compatible with geometry used to create the projector")
+
+ return True
+
+ def set_ImageGeometry(self, volume_geometry):
+
+ DataOrder.check_order_for_engine('astra', volume_geometry)
+
+ if len(volume_geometry.dimension_labels) > 3:
+ raise ValueError("Supports 2D and 3D data only, got {0}".format(volume_geometry.number_of_dimensions))
+
+ self.volume_geometry = volume_geometry.copy()
+
+ def set_AcquisitionGeometry(self, sinogram_geometry):
+
+ DataOrder.check_order_for_engine('astra', sinogram_geometry)
+
+ if len(sinogram_geometry.dimension_labels) > 3:
+ raise ValueError("Supports 2D and 3D data only, got {0}".format(self.volume_geometry.number_of_dimensions))
+
+ self.sinogram_geometry = sinogram_geometry.copy()
+
+ def process(self, out=None):
+
+ DATA = self.get_input()
+
+ data_temp = DATA.as_array()
+
+ arr_out = self.create_backprojection(data_temp, out)
+
+ if out is None:
+ out = ImageData(arr_out, deep_copy=False, geometry=self.volume_geometry.copy(), suppress_warning=True)
+ return out
+ else:
+ out.fill(arr_out)
+
+ def create_backprojection(self, sino, vol):
+
+ pm = Ufo.PluginManager()
+ graph = Ufo.TaskGraph()
+ scheduler = Ufo.Scheduler()
+
+ read = pm.get_task('memory-in')
+ read.props.pointer = sino.__array_interface__['data'][0]
+ read.props.width = sino.shape[0]
+ read.props.height = sino.shape[1]
+ read.props.number = sino.shape[2]
+ read.props.bitdepth = 32
+
+ stack = pm.get_task('stack')
+ stack.props.number = self.stack_num
+
+ stacked_bp = pm.get_task('stacked-backproject')
+ stacked_bp.props.axis_pos = sino.shape[1] / 2
+ stacked_bp.props.precision_mode = self.precision_mode
+
+ vol_arr = np.zeros(self.volume_geometry.shape, dtype=np.float32)
+
+ write = pm.get_task('memory-out')
+ write.props.pointer = vol_arr.__array_interface__['data'][0]
+ write.props.max_size = vol_arr.nbytes
+
+ graph.connect_nodes(read, stack)
+ graph.connect_nodes(stack, stacked_bp)
+ graph.connect_nodes(stacked_bp, write)
+
+ scheduler.run(graph)
+
+ return vol_arr
diff --git a/src/processors/UfoStackedForwardProjector.py b/src/processors/UfoStackedForwardProjector.py
new file mode 100644
index 0000000..86513c6
--- /dev/null
+++ b/src/processors/UfoStackedForwardProjector.py
@@ -0,0 +1,100 @@
+from cil.framework import DataProcessor, AcquisitionData, DataOrder
+
+import gi
+
+gi.require_version('Ufo', '0.0')
+from gi.repository import Ufo
+
+import numpy as np
+
+
+class UfoStackedForwardProjector(DataProcessor):
+ def __init__(self,
+ volume_geometry=None,
+ sinogram_geometry=None,
+ precision_mode=None,
+ stack_num=None):
+
+ kwargs = {
+ 'volume_geometry': volume_geometry,
+ 'sinogram_geometry': sinogram_geometry,
+ 'precision_mode': precision_mode,
+ 'stack_num': stack_num
+ }
+
+ super(UfoStackedForwardProjector, self).__init__(**kwargs)
+
+ self.precision_mode = precision_mode
+ self.stack_num = stack_num
+ self.set_ImageGeometry(volume_geometry)
+ self.set_AcquisitionGeometry(sinogram_geometry)
+
+ def check_input(self, dataset):
+
+ if self.volume_geometry.shape != dataset.geometry.shape:
+ raise ValueError("Dataset not compatible with geometry used to create the projector")
+
+ return True
+
+ def set_ImageGeometry(self, volume_geometry):
+
+ DataOrder.check_order_for_engine('astra', volume_geometry)
+
+ if len(volume_geometry.dimension_labels) > 3:
+ raise ValueError("Supports 2D and 3D data only, got {0}".format(volume_geometry.number_of_dimensions))
+
+ self.volume_geometry = volume_geometry.copy()
+
+ def set_AcquisitionGeometry(self, sinogram_geometry):
+
+ DataOrder.check_order_for_engine('astra', sinogram_geometry)
+
+ if len(sinogram_geometry.dimension_labels) > 3:
+ raise ValueError("Supports 2D and 3D data only, got {0}".format(self.volume_geometry.number_of_dimensions))
+
+ self.sinogram_geometry = sinogram_geometry.copy()
+
+ def process(self, out=None):
+
+ IM = self.get_input()
+ data_temp = IM.as_array()
+
+ arr_out = self.create_sinogram(data_temp, out)
+
+ if out is None:
+ out = AcquisitionData(arr_out, deep_copy=False, geometry=self.sinogram_geometry.copy(),suppress_warning=True)
+ return out
+ else:
+ out.fill(arr_out)
+
+ def create_sinogram(self, vol, sino):
+ pm = Ufo.PluginManager()
+ graph = Ufo.TaskGraph()
+ scheduler = Ufo.Scheduler()
+
+ read = pm.get_task('memory-in')
+ read.props.pointer = vol.__array_interface__['data'][0]
+ read.props.width = vol.shape[0]
+ read.props.height = vol.shape[1]
+ read.props.number = vol.shape[2]
+ read.props.bitdepth = 32
+
+ stack = pm.get_task('stack')
+ stack.props.number = self.stack_num
+
+ stacked_fp = pm.get_task('stacked-forwardproject')
+ stacked_fp.props.number = vol.shape[0]
+ stacked_fp.props.precision_mode = self.precision_mode
+
+ sino_arr = np.zeros(self.sinogram_geometry.shape, dtype=np.float32)
+
+ write = pm.get_task('memory-out')
+ write.props.pointer = sino_arr.__array_interface__['data'][0]
+ write.props.max_size = sino_arr.nbytes
+
+ graph.connect_nodes(read, stack)
+ graph.connect_nodes(stack, stacked_fp)
+ graph.connect_nodes(stacked_fp, write)
+
+ scheduler.run(graph)
+ return sino_arr
diff --git a/src/processors/UfoStandardBackProjector.py b/src/processors/UfoStandardBackProjector.py
new file mode 100644
index 0000000..ee27753
--- /dev/null
+++ b/src/processors/UfoStandardBackProjector.py
@@ -0,0 +1,93 @@
+#import tifffile
+from cil.framework import DataProcessor, ImageData, DataOrder
+import gi
+
+gi.require_version('Ufo', '0.0')
+from gi.repository import Ufo
+import numpy as np
+
+
+class UfoStandardBackProjector(DataProcessor):
+
+ def __init__(self,
+ volume_geometry=None,
+ sinogram_geometry=None):
+ kwargs = {
+ 'volume_geometry': volume_geometry,
+ 'sinogram_geometry': sinogram_geometry}
+
+ super(UfoStandardBackProjector, self).__init__(**kwargs)
+
+ self.set_ImageGeometry(volume_geometry)
+ self.set_AcquisitionGeometry(sinogram_geometry)
+
+ def check_input(self, dataset):
+
+ if self.sinogram_geometry.shape != dataset.geometry.shape:
+ raise ValueError("Dataset not compatible with geometry used to create the projector")
+
+ return True
+
+ def set_ImageGeometry(self, volume_geometry):
+
+ DataOrder.check_order_for_engine('astra', volume_geometry)
+
+ if len(volume_geometry.dimension_labels) > 3:
+ raise ValueError("Supports 2D and 3D data only, got {0}".format(volume_geometry.number_of_dimensions))
+
+ self.volume_geometry = volume_geometry.copy()
+
+ def set_AcquisitionGeometry(self, sinogram_geometry):
+
+ DataOrder.check_order_for_engine('astra', sinogram_geometry)
+
+ if len(sinogram_geometry.dimension_labels) > 3:
+ raise ValueError("Supports 2D and 3D data only, got {0}".format(self.volume_geometry.number_of_dimensions))
+
+ self.sinogram_geometry = sinogram_geometry.copy()
+
+ def process(self, out=None):
+
+ DATA = self.get_input()
+
+ data_temp = DATA.as_array()
+
+ arr_out = self.create_backprojection(data_temp, out)
+
+ if out is None:
+ out = ImageData(arr_out, deep_copy=False, geometry=self.volume_geometry.copy(), suppress_warning=True)
+ return out
+ else:
+ out.fill(arr_out)
+
+ def create_backprojection(self, sino, vol):
+
+ pm = Ufo.PluginManager()
+ graph = Ufo.TaskGraph()
+ scheduler = Ufo.Scheduler()
+
+ read = pm.get_task('memory-in')
+ read.props.pointer = sino.__array_interface__['data'][0]
+ read.props.width = sino.shape[1]
+ read.props.height = sino.shape[0]
+ if(len(sino.shape)==2):
+ read.props.number = 1
+ else:
+ read.props.number = sino.shape[2]
+ read.props.bitdepth = 32
+
+ backproject = pm.get_task('backproject')
+ backproject.props.axis_pos = sino.shape[1] / 2
+
+ vol_arr = np.zeros(self.volume_geometry.shape, dtype=np.float32)
+
+ write = pm.get_task('memory-out')
+ write.props.pointer = vol_arr.__array_interface__['data'][0]
+ write.props.max_size = vol_arr.nbytes
+
+ graph.connect_nodes(read, backproject)
+ graph.connect_nodes(backproject, write)
+
+ scheduler.run(graph)
+
+ return vol_arr
diff --git a/src/processors/UfoStandardForwardProjector.py b/src/processors/UfoStandardForwardProjector.py
new file mode 100644
index 0000000..bb6e0b4
--- /dev/null
+++ b/src/processors/UfoStandardForwardProjector.py
@@ -0,0 +1,89 @@
+from cil.framework import DataProcessor, AcquisitionData, DataOrder
+
+import gi
+gi.require_version('Ufo', '0.0')
+from gi.repository import Ufo
+
+import numpy as np
+
+class UfoStandardForwardProjector(DataProcessor):
+ def __init__(self,
+ volume_geometry=None,
+ sinogram_geometry=None):
+ kwargs = {
+ 'volume_geometry': volume_geometry,
+ 'sinogram_geometry': sinogram_geometry
+ }
+
+ super(UfoStandardForwardProjector, self).__init__(**kwargs)
+
+ self.set_ImageGeometry(volume_geometry)
+ self.set_AcquisitionGeometry(sinogram_geometry)
+
+ def check_input(self, dataset):
+
+ if self.volume_geometry.shape != dataset.geometry.shape:
+ raise ValueError("Dataset not compatible with geometry used to create the projector")
+
+ return True
+
+ def set_ImageGeometry(self, volume_geometry):
+
+ DataOrder.check_order_for_engine('astra', volume_geometry)
+
+ if len(volume_geometry.dimension_labels) > 3:
+ raise ValueError("Supports 2D and 3D data only, got {0}".format(volume_geometry.number_of_dimensions))
+
+ self.volume_geometry = volume_geometry.copy()
+
+ def set_AcquisitionGeometry(self, sinogram_geometry):
+
+ DataOrder.check_order_for_engine('astra', sinogram_geometry)
+
+ if len(sinogram_geometry.dimension_labels) > 3:
+ raise ValueError("Supports 2D and 3D data only, got {0}".format(self.volume_geometry.number_of_dimensions))
+
+ self.sinogram_geometry = sinogram_geometry.copy()
+
+ def process(self, out=None):
+
+ IM = self.get_input()
+ data_temp = IM.as_array()
+
+ arr_out = self.create_sinogram(data_temp, out)
+
+ if out is None:
+ out = AcquisitionData(arr_out, deep_copy=False, geometry=self.sinogram_geometry.copy(),suppress_warning=True)
+ return out
+ else:
+ out.fill(arr_out)
+
+ def create_sinogram(self, vol, sino):
+ pm = Ufo.PluginManager()
+ graph = Ufo.TaskGraph()
+ scheduler = Ufo.Scheduler()
+
+ read = pm.get_task('memory-in')
+ read.props.pointer = vol.__array_interface__['data'][0]
+ read.props.width = vol.shape[1]
+ read.props.height = vol.shape[0]
+ if(len(vol.shape) == 2):
+ read.props.number = 1
+ else:
+ read.props.number = vol.shape[2]
+ read.props.bitdepth = 32
+
+ forwardproject = pm.get_task('forwardproject')
+ forwardproject.props.number = vol.shape[0]
+
+ sino_arr = np.zeros(self.sinogram_geometry.shape, dtype=np.float32)
+
+ write = pm.get_task('memory-out')
+ write.props.pointer = sino_arr.__array_interface__['data'][0]
+ write.props.max_size = sino_arr.nbytes
+
+ graph.connect_nodes(read, forwardproject)
+ graph.connect_nodes(forwardproject, write)
+
+ scheduler.run(graph)
+ return sino_arr
diff --git a/src/processors/__init__.py b/src/processors/__init__.py
new file mode 100644
index 0000000..c913057
--- /dev/null
+++ b/src/processors/__init__.py
@@ -0,0 +1,5 @@
+from .UfoStandardForwardProjector import UfoStandardForwardProjector
+from .UfoStandardBackProjector import UfoStandardBackProjector
+from .UfoStackedForwardProjector import UfoStackedForwardProjector
+from .UfoStackedBackProjector import UfoStackedBackProjector
+from .FBP import FBP
--
cgit v1.2.3