diff options
| -rw-r--r-- | dummy.py | 132 | ||||
| -rw-r--r-- | src/__init__.py | 2 | ||||
| -rw-r--r-- | src/operators/ProjectionOperator.py | 67 | ||||
| -rw-r--r-- | src/operators/UfoStackedProjector.py | 60 | ||||
| -rw-r--r-- | src/operators/UfoStandardProjector.py | 58 | ||||
| -rw-r--r-- | src/operators/__init__.py | 3 | ||||
| -rw-r--r-- | src/processors/FBP.py | 32 | ||||
| -rw-r--r-- | src/processors/FBP_Stacked.py | 85 | ||||
| -rw-r--r-- | src/processors/FBP_Standard.py | 80 | ||||
| -rw-r--r-- | src/processors/UfoStackedBackProjector.py | 101 | ||||
| -rw-r--r-- | src/processors/UfoStackedForwardProjector.py | 100 | ||||
| -rw-r--r-- | src/processors/UfoStandardBackProjector.py | 93 | ||||
| -rw-r--r-- | src/processors/UfoStandardForwardProjector.py | 89 | ||||
| -rw-r--r-- | src/processors/__init__.py | 5 | 
14 files changed, 907 insertions, 0 deletions
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  | 
