diff options
| author | Edoardo Pasca <edo.paskino@gmail.com> | 2018-04-11 16:18:20 +0100 |
|---|---|---|
| committer | Edoardo Pasca <edo.paskino@gmail.com> | 2018-04-11 16:18:20 +0100 |
| commit | ee0a3bf634e92351ad82c09b4e94bcb64ee52d73 (patch) | |
| tree | 246466b22b14f20e46df2d50104f72e23cb2b782 | |
| parent | a371b63586e0d0b0b4b13b7b3038e0d0cacbe0ca (diff) | |
| download | framework-plugins-ee0a3bf634e92351ad82c09b4e94bcb64ee52d73.tar.gz framework-plugins-ee0a3bf634e92351ad82c09b4e94bcb64ee52d73.tar.bz2 framework-plugins-ee0a3bf634e92351ad82c09b4e94bcb64ee52d73.tar.xz framework-plugins-ee0a3bf634e92351ad82c09b4e94bcb64ee52d73.zip | |
initial working release with simple_demo_ccpi.py
| -rw-r--r-- | Wrappers/Python/ccpi/plugins/__init__.py | 0 | ||||
| -rwxr-xr-x | Wrappers/Python/ccpi/plugins/ops.py | 102 | ||||
| -rwxr-xr-x | Wrappers/Python/ccpi/plugins/processors.py | 792 | ||||
| -rw-r--r-- | Wrappers/Python/conda-recipe/bld.bat | 10 | ||||
| -rw-r--r-- | Wrappers/Python/conda-recipe/build.sh | 9 | ||||
| -rw-r--r-- | Wrappers/Python/conda-recipe/meta.yaml | 27 | ||||
| -rw-r--r-- | Wrappers/Python/setup.py | 58 | ||||
| -rwxr-xr-x | Wrappers/Python/wip/simple_demo_ccpi.py | 272 |
8 files changed, 1270 insertions, 0 deletions
diff --git a/Wrappers/Python/ccpi/plugins/__init__.py b/Wrappers/Python/ccpi/plugins/__init__.py new file mode 100644 index 0000000..e69de29 --- /dev/null +++ b/Wrappers/Python/ccpi/plugins/__init__.py diff --git a/Wrappers/Python/ccpi/plugins/ops.py b/Wrappers/Python/ccpi/plugins/ops.py new file mode 100755 index 0000000..0028cf1 --- /dev/null +++ b/Wrappers/Python/ccpi/plugins/ops.py @@ -0,0 +1,102 @@ +# -*- coding: utf-8 -*-
+# This work is part of the Core Imaging Library developed by
+# Visual Analytics and Imaging System Group of the Science Technology
+# Facilities Council, STFC
+
+# Copyright 2018 Jakob Jorgensen, Daniil Kazantsev and Edoardo Pasca
+
+# 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
+
+# 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.
+
+import numpy
+from ccpi.optimisation.ops import Operator, PowerMethodNonsquare
+from ccpi.framework import ImageData, DataContainer
+from ccpi.plugins.processors import CCPiBackwardProjector, CCPiForwardProjector
+
+class LinearOperatorMatrix(Operator):
+ def __init__(self,A):
+ self.A = A
+ self.s1 = None # Largest singular value, initially unknown
+ super(LinearOperatorMatrix, self).__init__()
+
+ def direct(self,x):
+ return DataContainer(numpy.dot(self.A,x.as_array()))
+
+ def adjoint(self,x):
+ return DataContainer(numpy.dot(self.A.transpose(),x.as_array()))
+
+ def size(self):
+ return self.A.shape
+
+ def get_max_sing_val(self):
+ # If unknown, compute and store. If known, simply return it.
+ if self.s1 is None:
+ self.s1 = svds(self.A,1,return_singular_vectors=False)[0]
+ return self.s1
+ else:
+ return self.s1
+
+
+
+class CCPiProjectorSimple(Operator):
+ """ASTRA projector modified to use DataSet and geometry."""
+ def __init__(self, geomv, geomp):
+ super(CCPiProjectorSimple, self).__init__()
+
+ # Store volume and sinogram geometries.
+ self.acquisition_geometry = geomp
+ self.volume_geometry = geomv
+
+ self.fp = CCPiForwardProjector(image_geometry=geomv,
+ acquisition_geometry=geomp,
+ output_axes_order=['angle','vertical','horizontal'])
+
+ self.bp = CCPiBackwardProjector(image_geometry=geomv,
+ acquisition_geometry=geomp,
+ output_axes_order=['horizontal_x','horizontal_y','vertical'])
+
+ # Initialise empty for singular value.
+ self.s1 = None
+
+ def direct(self, image_data):
+ self.fp.set_input(image_data)
+ out = self.fp.get_output()
+ return out
+
+ def adjoint(self, acquisition_data):
+ self.bp.set_input(acquisition_data)
+ out = self.bp.get_output()
+ return out
+
+ #def delete(self):
+ # astra.data2d.delete(self.proj_id)
+
+ def get_max_sing_val(self):
+ a = PowerMethodNonsquare(self,10)
+ self.s1 = a[0]
+ return self.s1
+
+ def size(self):
+ # Only implemented for 3D
+ return ( (self.acquisition_geometry.angles.size, \
+ self.acquisition_geometry.pixel_num_v,
+ self.acquisition_geometry.pixel_num_h), \
+ (self.volume_geometry.voxel_num_x, \
+ self.volume_geometry.voxel_num_y,
+ self.volume_geometry.voxel_num_z) )
+ def create_image_data(self):
+ x0 = ImageData(geometry = self.volume_geometry,
+ dimension_labels=self.bp.output_axes_order)#\
+ #.subset(['horizontal_x','horizontal_y','vertical'])
+ print (x0)
+ x0.fill(numpy.random.randn(*x0.shape))
+ return x0
\ No newline at end of file diff --git a/Wrappers/Python/ccpi/plugins/processors.py b/Wrappers/Python/ccpi/plugins/processors.py new file mode 100755 index 0000000..2fd3bf1 --- /dev/null +++ b/Wrappers/Python/ccpi/plugins/processors.py @@ -0,0 +1,792 @@ +# -*- coding: utf-8 -*-
+# This work is part of the Core Imaging Library developed by
+# Visual Analytics and Imaging System Group of the Science Technology
+# Facilities Council, STFC
+
+# Copyright 2018 Edoardo Pasca
+
+# 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
+
+# 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 ccpi.framework import DataProcessor, DataContainer, AcquisitionData,\
+ AcquisitionGeometry, ImageGeometry, ImageData
+from ccpi.reconstruction.parallelbeam import alg as pbalg
+import numpy
+import h5py
+from scipy import ndimage
+
+import matplotlib.pyplot as plt
+
+
+class Normalizer(DataProcessor):
+ '''Normalization based on flat and dark
+
+ This processor read in a AcquisitionData and normalises it based on
+ the instrument reading with and without incident photons or neutrons.
+
+ Input: AcquisitionData
+ Parameter: 2D projection with flat field (or stack)
+ 2D projection with dark field (or stack)
+ Output: AcquisitionDataSetn
+ '''
+
+ def __init__(self, flat_field = None, dark_field = None, tolerance = 1e-5):
+ kwargs = {
+ 'flat_field' : None,
+ 'dark_field' : None,
+ # very small number. Used when there is a division by zero
+ 'tolerance' : tolerance
+ }
+
+ #DataProcessor.__init__(self, **kwargs)
+ super(Normalizer, self).__init__(**kwargs)
+ if not flat_field is None:
+ self.set_flat_field(flat_field)
+ if not dark_field is None:
+ self.set_dark_field(dark_field)
+
+ def check_input(self, dataset):
+ if dataset.number_of_dimensions == 3:
+ return True
+ else:
+ raise ValueError("Expected input dimensions is 2 or 3, got {0}"\
+ .format(dataset.number_of_dimensions))
+
+ def set_dark_field(self, df):
+ if type(df) is numpy.ndarray:
+ if len(numpy.shape(df)) == 3:
+ raise ValueError('Dark Field should be 2D')
+ elif len(numpy.shape(df)) == 2:
+ self.dark_field = df
+ elif issubclass(type(df), DataSet):
+ self.dark_field = self.set_dark_field(df.as_array())
+
+ def set_flat_field(self, df):
+ if type(df) is numpy.ndarray:
+ if len(numpy.shape(df)) == 3:
+ raise ValueError('Flat Field should be 2D')
+ elif len(numpy.shape(df)) == 2:
+ self.flat_field = df
+ elif issubclass(type(df), DataSet):
+ self.flat_field = self.set_flat_field(df.as_array())
+
+ @staticmethod
+ def normalize_projection(projection, flat, dark, tolerance):
+ a = (projection - dark)
+ b = (flat-dark)
+ with numpy.errstate(divide='ignore', invalid='ignore'):
+ c = numpy.true_divide( a, b )
+ c[ ~ numpy.isfinite( c )] = tolerance # set to not zero if 0/0
+ return c
+
+ def process(self):
+
+ projections = self.get_input()
+ dark = self.dark_field
+ flat = self.flat_field
+
+ if not (projections.shape[1:] == dark.shape and \
+ projections.shape[1:] == flat.shape):
+ raise ValueError('Flats/Dark and projections size do not match.')
+
+
+ a = numpy.asarray(
+ [ Normalizer.normalize_projection(
+ projection, flat, dark, self.tolerance) \
+ for projection in projections.as_array() ]
+ )
+ y = type(projections)( a , True,
+ dimension_labels=projections.dimension_labels,
+ geometry=projections.geometry)
+ return y
+
+
+class CenterOfRotationFinder(DataProcessor):
+ '''Processor to find the center of rotation in a parallel beam experiment
+
+ This processor read in a AcquisitionDataSet and finds the center of rotation
+ based on Nghia Vo's method. https://doi.org/10.1364/OE.22.019078
+
+ Input: AcquisitionDataSet
+
+ Output: float. center of rotation in pixel coordinate
+ '''
+
+ def __init__(self):
+ kwargs = {
+
+ }
+
+ #DataProcessor.__init__(self, **kwargs)
+ super(CenterOfRotationFinder, self).__init__(**kwargs)
+
+ def check_input(self, dataset):
+ if dataset.number_of_dimensions == 3:
+ if dataset.geometry.geom_type == 'parallel':
+ return True
+ else:
+ raise ValueError('{0} is suitable only for parallel beam geometry'\
+ .format(self.__class__.__name__))
+ else:
+ raise ValueError("Expected input dimensions is 3, got {0}"\
+ .format(dataset.number_of_dimensions))
+
+
+ # #########################################################################
+ # Copyright (c) 2015, UChicago Argonne, LLC. All rights reserved. #
+ # #
+ # Copyright 2015. UChicago Argonne, LLC. This software was produced #
+ # under U.S. Government contract DE-AC02-06CH11357 for Argonne National #
+ # Laboratory (ANL), which is operated by UChicago Argonne, LLC for the #
+ # U.S. Department of Energy. The U.S. Government has rights to use, #
+ # reproduce, and distribute this software. NEITHER THE GOVERNMENT NOR #
+ # UChicago Argonne, LLC MAKES ANY WARRANTY, EXPRESS OR IMPLIED, OR #
+ # ASSUMES ANY LIABILITY FOR THE USE OF THIS SOFTWARE. If software is #
+ # modified to produce derivative works, such modified software should #
+ # be clearly marked, so as not to confuse it with the version available #
+ # from ANL. #
+ # #
+ # Additionally, redistribution and use in source and binary forms, with #
+ # or without modification, are permitted provided that the following #
+ # conditions are met: #
+ # #
+ # * Redistributions of source code must retain the above copyright #
+ # notice, this list of conditions and the following disclaimer. #
+ # #
+ # * Redistributions in binary form must reproduce the above copyright #
+ # notice, this list of conditions and the following disclaimer in #
+ # the documentation and/or other materials provided with the #
+ # distribution. #
+ # #
+ # * Neither the name of UChicago Argonne, LLC, Argonne National #
+ # Laboratory, ANL, the U.S. Government, nor the names of its #
+ # contributors may be used to endorse or promote products derived #
+ # from this software without specific prior written permission. #
+ # #
+ # THIS SOFTWARE IS PROVIDED BY UChicago Argonne, LLC AND CONTRIBUTORS #
+ # "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT #
+ # LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS #
+ # FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL UChicago #
+ # Argonne, LLC OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, #
+ # INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, #
+ # BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; #
+ # LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER #
+ # CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT #
+ # LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN #
+ # ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE #
+ # POSSIBILITY OF SUCH DAMAGE. #
+ # #########################################################################
+
+ @staticmethod
+ def as_ndarray(arr, dtype=None, copy=False):
+ if not isinstance(arr, numpy.ndarray):
+ arr = numpy.array(arr, dtype=dtype, copy=copy)
+ return arr
+
+ @staticmethod
+ def as_dtype(arr, dtype, copy=False):
+ if not arr.dtype == dtype:
+ arr = numpy.array(arr, dtype=dtype, copy=copy)
+ return arr
+
+ @staticmethod
+ def as_float32(arr):
+ arr = CenterOfRotationFinder.as_ndarray(arr, numpy.float32)
+ return CenterOfRotationFinder.as_dtype(arr, numpy.float32)
+
+
+
+
+ @staticmethod
+ def find_center_vo(tomo, ind=None, smin=-40, smax=40, srad=10, step=0.5,
+ ratio=2., drop=20):
+ """
+ Find rotation axis location using Nghia Vo's method. :cite:`Vo:14`.
+
+ Parameters
+ ----------
+ tomo : ndarray
+ 3D tomographic data.
+ ind : int, optional
+ Index of the slice to be used for reconstruction.
+ smin, smax : int, optional
+ Reference to the horizontal center of the sinogram.
+ srad : float, optional
+ Fine search radius.
+ step : float, optional
+ Step of fine searching.
+ ratio : float, optional
+ The ratio between the FOV of the camera and the size of object.
+ It's used to generate the mask.
+ drop : int, optional
+ Drop lines around vertical center of the mask.
+
+ Returns
+ -------
+ float
+ Rotation axis location.
+
+ Notes
+ -----
+ The function may not yield a correct estimate, if:
+
+ - the sample size is bigger than the field of view of the camera.
+ In this case the ``ratio`` argument need to be set larger
+ than the default of 2.0.
+
+ - there is distortion in the imaging hardware. If there's
+ no correction applied, the center of the projection image may
+ yield a better estimate.
+
+ - the sample contrast is weak. Paganin's filter need to be applied
+ to overcome this.
+
+ - the sample was changed during the scan.
+ """
+ tomo = CenterOfRotationFinder.as_float32(tomo)
+
+ if ind is None:
+ ind = tomo.shape[1] // 2
+ _tomo = tomo[:, ind, :]
+
+
+
+ # Reduce noise by smooth filters. Use different filters for coarse and fine search
+ _tomo_cs = ndimage.filters.gaussian_filter(_tomo, (3, 1))
+ _tomo_fs = ndimage.filters.median_filter(_tomo, (2, 2))
+
+ # Coarse and fine searches for finding the rotation center.
+ if _tomo.shape[0] * _tomo.shape[1] > 4e6: # If data is large (>2kx2k)
+ #_tomo_coarse = downsample(numpy.expand_dims(_tomo_cs,1), level=2)[:, 0, :]
+ #init_cen = _search_coarse(_tomo_coarse, smin, smax, ratio, drop)
+ #fine_cen = _search_fine(_tomo_fs, srad, step, init_cen*4, ratio, drop)
+ init_cen = CenterOfRotationFinder._search_coarse(_tomo_cs, smin,
+ smax, ratio, drop)
+ fine_cen = CenterOfRotationFinder._search_fine(_tomo_fs, srad,
+ step, init_cen,
+ ratio, drop)
+ else:
+ init_cen = CenterOfRotationFinder._search_coarse(_tomo_cs,
+ smin, smax,
+ ratio, drop)
+ fine_cen = CenterOfRotationFinder._search_fine(_tomo_fs, srad,
+ step, init_cen,
+ ratio, drop)
+
+ #logger.debug('Rotation center search finished: %i', fine_cen)
+ return fine_cen
+
+
+ @staticmethod
+ def _search_coarse(sino, smin, smax, ratio, drop):
+ """
+ Coarse search for finding the rotation center.
+ """
+ (Nrow, Ncol) = sino.shape
+ centerfliplr = (Ncol - 1.0) / 2.0
+
+ # Copy the sinogram and flip left right, the purpose is to
+ # make a full [0;2Pi] sinogram
+ _copy_sino = numpy.fliplr(sino[1:])
+
+ # This image is used for compensating the shift of sinogram 2
+ temp_img = numpy.zeros((Nrow - 1, Ncol), dtype='float32')
+ temp_img[:] = sino[-1]
+
+ # Start coarse search in which the shift step is 1
+ listshift = numpy.arange(smin, smax + 1)
+ listmetric = numpy.zeros(len(listshift), dtype='float32')
+ mask = CenterOfRotationFinder._create_mask(2 * Nrow - 1, Ncol,
+ 0.5 * ratio * Ncol, drop)
+ for i in listshift:
+ _sino = numpy.roll(_copy_sino, i, axis=1)
+ if i >= 0:
+ _sino[:, 0:i] = temp_img[:, 0:i]
+ else:
+ _sino[:, i:] = temp_img[:, i:]
+ listmetric[i - smin] = numpy.sum(numpy.abs(numpy.fft.fftshift(
+ #pyfftw.interfaces.numpy_fft.fft2(
+ # numpy.vstack((sino, _sino)))
+ numpy.fft.fft2(numpy.vstack((sino, _sino)))
+ )) * mask)
+ minpos = numpy.argmin(listmetric)
+ return centerfliplr + listshift[minpos] / 2.0
+
+ @staticmethod
+ def _search_fine(sino, srad, step, init_cen, ratio, drop):
+ """
+ Fine search for finding the rotation center.
+ """
+ Nrow, Ncol = sino.shape
+ centerfliplr = (Ncol + 1.0) / 2.0 - 1.0
+ # Use to shift the sinogram 2 to the raw CoR.
+ shiftsino = numpy.int16(2 * (init_cen - centerfliplr))
+ _copy_sino = numpy.roll(numpy.fliplr(sino[1:]), shiftsino, axis=1)
+ if init_cen <= centerfliplr:
+ lefttake = numpy.int16(numpy.ceil(srad + 1))
+ righttake = numpy.int16(numpy.floor(2 * init_cen - srad - 1))
+ else:
+ lefttake = numpy.int16(numpy.ceil(
+ init_cen - (Ncol - 1 - init_cen) + srad + 1))
+ righttake = numpy.int16(numpy.floor(Ncol - 1 - srad - 1))
+ Ncol1 = righttake - lefttake + 1
+ mask = CenterOfRotationFinder._create_mask(2 * Nrow - 1, Ncol1,
+ 0.5 * ratio * Ncol, drop)
+ numshift = numpy.int16((2 * srad) / step) + 1
+ listshift = numpy.linspace(-srad, srad, num=numshift)
+ listmetric = numpy.zeros(len(listshift), dtype='float32')
+ factor1 = numpy.mean(sino[-1, lefttake:righttake])
+ num1 = 0
+ for i in listshift:
+ _sino = ndimage.interpolation.shift(
+ _copy_sino, (0, i), prefilter=False)
+ factor2 = numpy.mean(_sino[0,lefttake:righttake])
+ _sino = _sino * factor1 / factor2
+ sinojoin = numpy.vstack((sino, _sino))
+ listmetric[num1] = numpy.sum(numpy.abs(numpy.fft.fftshift(
+ #pyfftw.interfaces.numpy_fft.fft2(
+ # sinojoin[:, lefttake:righttake + 1])
+ numpy.fft.fft2(sinojoin[:, lefttake:righttake + 1])
+ )) * mask)
+ num1 = num1 + 1
+ minpos = numpy.argmin(listmetric)
+ return init_cen + listshift[minpos] / 2.0
+
+ @staticmethod
+ def _create_mask(nrow, ncol, radius, drop):
+ du = 1.0 / ncol
+ dv = (nrow - 1.0) / (nrow * 2.0 * numpy.pi)
+ centerrow = numpy.ceil(nrow / 2) - 1
+ centercol = numpy.ceil(ncol / 2) - 1
+ # added by Edoardo Pasca
+ centerrow = int(centerrow)
+ centercol = int(centercol)
+ mask = numpy.zeros((nrow, ncol), dtype='float32')
+ for i in range(nrow):
+ num1 = numpy.round(((i - centerrow) * dv / radius) / du)
+ (p1, p2) = numpy.int16(numpy.clip(numpy.sort(
+ (-num1 + centercol, num1 + centercol)), 0, ncol - 1))
+ mask[i, p1:p2 + 1] = numpy.ones(p2 - p1 + 1, dtype='float32')
+ if drop < centerrow:
+ mask[centerrow - drop:centerrow + drop + 1,
+ :] = numpy.zeros((2 * drop + 1, ncol), dtype='float32')
+ mask[:,centercol-1:centercol+2] = numpy.zeros((nrow, 3), dtype='float32')
+ return mask
+
+ def process(self):
+
+ projections = self.get_input()
+
+ cor = CenterOfRotationFinder.find_center_vo(projections.as_array())
+
+ return cor
+
+class CCPiForwardProjector(DataProcessor):
+ '''Normalization based on flat and dark
+
+ This processor read in a AcquisitionData and normalises it based on
+ the instrument reading with and without incident photons or neutrons.
+
+ Input: AcquisitionData
+ Parameter: 2D projection with flat field (or stack)
+ 2D projection with dark field (or stack)
+ Output: AcquisitionDataSetn
+ '''
+
+ def __init__(self,
+ image_geometry = None,
+ acquisition_geometry = None,
+ output_axes_order = None):
+ if output_axes_order is None:
+ # default ccpi projector image storing order
+ output_axes_order = ['angle','vertical','horizontal']
+
+ kwargs = {
+ 'image_geometry' : image_geometry,
+ 'acquisition_geometry' : acquisition_geometry,
+ 'output_axes_order' : output_axes_order,
+ 'default_image_axes_order' : ['horizontal_x','horizontal_y','vertical'],
+ 'default_acquisition_axes_order' : ['angle','vertical','horizontal']
+ }
+
+ super(CCPiForwardProjector, self).__init__(**kwargs)
+
+ def check_input(self, dataset):
+ if dataset.number_of_dimensions == 3 or dataset.number_of_dimensions == 2:
+ # sort in the order that this projector needs it
+ return True
+ else:
+ raise ValueError("Expected input dimensions is 2 or 3, got {0}"\
+ .format(dataset.number_of_dimensions))
+
+ def process(self):
+
+ volume = self.get_input()
+ volume_axes = volume.get_data_axes_order(new_order=self.default_image_axes_order)
+ if not volume_axes == [0,1,2]:
+ volume.array = numpy.transpose(volume.array, volume_axes)
+ pixel_per_voxel = 1 # should be estimated from image_geometry and
+ # acquisition_geometry
+ if self.acquisition_geometry.geom_type == 'parallel':
+ #int msize = ndarray_volume.shape(0) > ndarray_volume.shape(1) ? ndarray_volume.shape(0) : ndarray_volume.shape(1);
+ #int detector_width = msize;
+ # detector_width is the max between the shape[0] and shape[1]
+
+
+ #double rotation_center = (double)detector_width/2.;
+ #int detector_height = ndarray_volume.shape(2);
+
+ #int number_of_projections = ndarray_angles.shape(0);
+
+ ##numpy_3d pixels(reinterpret_cast<float*>(ndarray_volume.get_data()),
+ #boost::extents[number_of_projections][detector_height][detector_width]);
+
+ pixels = pbalg.pb_forward_project(volume.as_array(),
+ self.acquisition_geometry.angles,
+ pixel_per_voxel)
+ out = AcquisitionData(geometry=self.acquisition_geometry,
+ label_dimensions=self.default_acquisition_axes_order)
+ out.fill(pixels)
+ out_axes = out.get_data_axes_order(new_order=self.output_axes_order)
+ if not out_axes == [0,1,2]:
+ out.array = numpy.transpose(out.array, out_axes)
+ return out
+ else:
+ raise ValueError('Cannot process cone beam')
+
+class CCPiBackwardProjector(DataProcessor):
+ '''Backward projector
+
+ This processor reads in a AcquisitionData and performs a backward projection,
+ i.e. project to reconstruction space.
+ Notice that it assumes that the center of rotation is in the middle
+ of the horizontal axis: in case when that's not the case it can be chained
+ with the AcquisitionDataPadder.
+
+ Input: AcquisitionData
+ Parameter: 2D projection with flat field (or stack)
+ 2D projection with dark field (or stack)
+ Output: AcquisitionDataSetn
+ '''
+
+ def __init__(self,
+ image_geometry = None,
+ acquisition_geometry = None,
+ output_axes_order=None):
+ if output_axes_order is None:
+ # default ccpi projector image storing order
+ output_axes_order = ['horizontal_x','horizontal_y','vertical']
+ kwargs = {
+ 'image_geometry' : image_geometry,
+ 'acquisition_geometry' : acquisition_geometry,
+ 'output_axes_order' : output_axes_order,
+ 'default_image_axes_order' : ['horizontal_x','horizontal_y','vertical'],
+ 'default_acquisition_axes_order' : ['angle','vertical','horizontal']
+ }
+
+ super(CCPiBackwardProjector, self).__init__(**kwargs)
+
+ def check_input(self, dataset):
+ if dataset.number_of_dimensions == 3 or dataset.number_of_dimensions == 2:
+ #number_of_projections][detector_height][detector_width
+
+ return True
+ else:
+ raise ValueError("Expected input dimensions is 2 or 3, got {0}"\
+ .format(dataset.number_of_dimensions))
+
+ def process(self):
+ projections = self.get_input()
+ projections_axes = projections.get_data_axes_order(new_order=self.default_acquisition_axes_order)
+ if not projections_axes == [0,1,2]:
+ projections.array = numpy.transpose(projections.array, projections_axes)
+
+ pixel_per_voxel = 1 # should be estimated from image_geometry and acquisition_geometry
+ image_geometry = ImageGeometry(voxel_num_x = self.acquisition_geometry.pixel_num_h,
+ voxel_num_y = self.acquisition_geometry.pixel_num_h,
+ voxel_num_z = self.acquisition_geometry.pixel_num_v)
+ # input centered/padded acquisitiondata
+ center_of_rotation = projections.get_dimension_size('horizontal') / 2
+ #print (center_of_rotation)
+ if self.acquisition_geometry.geom_type == 'parallel':
+ back = pbalg.pb_backward_project(
+ projections.as_array(),
+ self.acquisition_geometry.angles,
+ center_of_rotation,
+ pixel_per_voxel
+ )
+ out = ImageData(geometry=self.image_geometry,
+ dimension_labels=self.default_image_axes_order)
+
+ out_axes = out.get_data_axes_order(new_order=self.output_axes_order)
+ if not out_axes == [0,1,2]:
+ back = numpy.transpose(back, out_axes)
+ out.fill(back)
+
+ return out
+
+ else:
+ raise ValueError('Cannot process cone beam')
+
+class AcquisitionDataPadder(DataProcessor):
+ '''Normalization based on flat and dark
+
+ This processor read in a AcquisitionData and normalises it based on
+ the instrument reading with and without incident photons or neutrons.
+
+ Input: AcquisitionData
+ Parameter: 2D projection with flat field (or stack)
+ 2D projection with dark field (or stack)
+ Output: AcquisitionDataSetn
+ '''
+
+ def __init__(self,
+ center_of_rotation = None,
+ acquisition_geometry = None,
+ pad_value = 1e-5):
+ kwargs = {
+ 'acquisition_geometry' : acquisition_geometry,
+ 'center_of_rotation' : center_of_rotation,
+ 'pad_value' : pad_value
+ }
+
+ super(AcquisitionDataPadder, self).__init__(**kwargs)
+
+ def check_input(self, dataset):
+ if self.acquisition_geometry is None:
+ self.acquisition_geometry = dataset.geometry
+ if dataset.number_of_dimensions == 3:
+ return True
+ else:
+ raise ValueError("Expected input dimensions is 2 or 3, got {0}"\
+ .format(dataset.number_of_dimensions))
+
+ def process(self):
+ projections = self.get_input()
+ w = projections.get_dimension_size('horizontal')
+ delta = w - 2 * self.center_of_rotation
+
+ padded_width = int (
+ numpy.ceil(abs(delta)) + w
+ )
+ delta_pix = padded_width - w
+
+ voxel_per_pixel = 1
+ geom = pbalg.pb_setup_geometry_from_acquisition(projections.as_array(),
+ self.acquisition_geometry.angles,
+ self.center_of_rotation,
+ voxel_per_pixel )
+
+ padded_geometry = self.acquisition_geometry.clone()
+
+ padded_geometry.pixel_num_h = geom['n_h']
+ padded_geometry.pixel_num_v = geom['n_v']
+
+ delta_pix_h = padded_geometry.pixel_num_h - self.acquisition_geometry.pixel_num_h
+ delta_pix_v = padded_geometry.pixel_num_v - self.acquisition_geometry.pixel_num_v
+
+ if delta_pix_h == 0:
+ delta_pix_h = delta_pix
+ padded_geometry.pixel_num_h = padded_width
+ #initialize a new AcquisitionData with values close to 0
+ out = AcquisitionData(geometry=padded_geometry)
+ out = out + self.pad_value
+
+
+ #pad in the horizontal-vertical plane -> slice on angles
+ if delta > 0:
+ #pad left of middle
+ command = "out.array["
+ for i in range(out.number_of_dimensions):
+ if out.dimension_labels[i] == 'horizontal':
+ value = '{0}:{1}'.format(delta_pix_h, delta_pix_h+w)
+ command = command + str(value)
+ else:
+ if out.dimension_labels[i] == 'vertical' :
+ value = '{0}:'.format(delta_pix_v)
+ command = command + str(value)
+ else:
+ command = command + ":"
+ if i < out.number_of_dimensions -1:
+ command = command + ','
+ command = command + '] = projections.array'
+ #print (command)
+ else:
+ #pad right of middle
+ command = "out.array["
+ for i in range(out.number_of_dimensions):
+ if out.dimension_labels[i] == 'horizontal':
+ value = '{0}:{1}'.format(0, w)
+ command = command + str(value)
+ else:
+ if out.dimension_labels[i] == 'vertical' :
+ value = '{0}:'.format(delta_pix_v)
+ command = command + str(value)
+ else:
+ command = command + ":"
+ if i < out.number_of_dimensions -1:
+ command = command + ','
+ command = command + '] = projections.array'
+ #print (command)
+ #cleaned = eval(command)
+ exec(command)
+ return out
+
+#class FiniteDifferentiator(DataProcessor):
+# def __init__(self):
+# kwargs = {
+#
+# }
+#
+# super(FiniteDifferentiator, self).__init__(**kwargs)
+#
+# def check_input(self, dataset):
+# return True
+#
+# def process(self):
+# axis = 0
+# d1 = numpy.diff(x,n=1,axis=axis)
+# d1 = numpy.resize(d1, numpy.shape(x))
+
+
+
+def loadNexus(filename):
+ '''Load a dataset stored in a NeXuS file (HDF5)'''
+ ###########################################################################
+ ## Load a dataset
+ nx = h5py.File(filename, "r")
+
+ data = nx.get('entry1/tomo_entry/data/rotation_angle')
+ angles = numpy.zeros(data.shape)
+ data.read_direct(angles)
+
+ data = nx.get('entry1/tomo_entry/data/data')
+ stack = numpy.zeros(data.shape)
+ data.read_direct(stack)
+ data = nx.get('entry1/tomo_entry/instrument/detector/image_key')
+
+ itype = numpy.zeros(data.shape)
+ data.read_direct(itype)
+ # 2 is dark field
+ darks = [stack[i] for i in range(len(itype)) if itype[i] == 2 ]
+ dark = darks[0]
+ for i in range(1, len(darks)):
+ dark += darks[i]
+ dark = dark / len(darks)
+ #dark[0][0] = dark[0][1]
+
+ # 1 is flat field
+ flats = [stack[i] for i in range(len(itype)) if itype[i] == 1 ]
+ flat = flats[0]
+ for i in range(1, len(flats)):
+ flat += flats[i]
+ flat = flat / len(flats)
+ #flat[0][0] = dark[0][1]
+
+
+ # 0 is projection data
+ proj = [stack[i] for i in range(len(itype)) if itype[i] == 0 ]
+ angle_proj = [angles[i] for i in range(len(itype)) if itype[i] == 0 ]
+ angle_proj = numpy.asarray (angle_proj)
+ angle_proj = angle_proj.astype(numpy.float32)
+
+ return angle_proj , numpy.asarray(proj) , dark, flat
+
+
+
+if __name__ == '__main__':
+ angles, proj, dark, flat = loadNexus('../../../data/24737_fd.nxs')
+
+ parallelbeam = AcquisitionGeometry('parallel', '3D' ,
+ angles=angles,
+ pixel_num_h=numpy.shape(proj)[2],
+ pixel_num_v=numpy.shape(proj)[1],
+ )
+
+ dim_labels = ['angles' , 'vertical' , 'horizontal']
+ sino = AcquisitionData( proj , geometry=parallelbeam,
+ dimension_labels=dim_labels)
+
+ normalizer = Normalizer()
+ normalizer.set_input(sino)
+ normalizer.set_flat_field(flat)
+ normalizer.set_dark_field(dark)
+ norm = normalizer.get_output()
+ #print ("Processor min {0} max {1}".format(norm.as_array().min(), norm.as_array().max()))
+
+ #norm1 = numpy.asarray(
+ # [Normalizer.normalize_projection( p, flat, dark, 1e-5 )
+ # for p in proj]
+ # )
+
+ #print ("Numpy min {0} max {1}".format(norm1.min(), norm1.max()))
+
+ cor_finder = CenterOfRotationFinder()
+ cor_finder.set_input(sino)
+ cor = cor_finder.get_output()
+ print ("center of rotation {0} == 86.25?".format(cor))
+
+
+ padder = AcquisitionDataPadder(center_of_rotation=cor,
+ pad_value=0.75,
+ acquisition_geometry=normalizer.get_output().geometry)
+ #padder.set_input(normalizer.get_output())
+ padder.set_input_processor(normalizer)
+
+ #print ("padder ", padder.get_output())
+
+ volume_geometry = ImageGeometry()
+
+ back = CCPiBackwardProjector(acquisition_geometry=parallelbeam,
+ image_geometry=volume_geometry)
+ back.set_input_processor(padder)
+ #back.set_input(padder.get_output())
+ #print (back.image_geometry)
+ forw = CCPiForwardProjector(acquisition_geometry=parallelbeam,
+ image_geometry=volume_geometry)
+ forw.set_input_processor(back)
+
+
+ out = padder.get_output()
+ fig = plt.figure()
+ # projections row
+ a = fig.add_subplot(1,4,1)
+ a.imshow(norm.array[80])
+ a.set_title('orig')
+ a = fig.add_subplot(1,4,2)
+ a.imshow(out.array[80])
+ a.set_title('padded')
+
+ a = fig.add_subplot(1,4,3)
+ a.imshow(back.get_output().as_array()[15])
+ a.set_title('back')
+
+ a = fig.add_subplot(1,4,4)
+ a.imshow(forw.get_output().as_array()[80])
+ a.set_title('forw')
+
+
+ plt.show()
+
+
+ conebeam = AcquisitionGeometry('cone', '3D' ,
+ angles=angles,
+ pixel_num_h=numpy.shape(proj)[2],
+ pixel_num_v=numpy.shape(proj)[1],
+ )
+
+ try:
+ sino2 = AcquisitionData( proj , geometry=conebeam)
+ cor_finder.set_input(sino2)
+ cor = cor_finder.get_output()
+ except ValueError as err:
+ print (err)
\ No newline at end of file diff --git a/Wrappers/Python/conda-recipe/bld.bat b/Wrappers/Python/conda-recipe/bld.bat new file mode 100644 index 0000000..4b4c7f5 --- /dev/null +++ b/Wrappers/Python/conda-recipe/bld.bat @@ -0,0 +1,10 @@ +IF NOT DEFINED CIL_VERSION ( +ECHO CIL_VERSION Not Defined. +exit 1 +) + +ROBOCOPY /E "%RECIPE_DIR%\.." "%SRC_DIR%" + +%PYTHON% setup.py build_py +%PYTHON% setup.py install +if errorlevel 1 exit 1 diff --git a/Wrappers/Python/conda-recipe/build.sh b/Wrappers/Python/conda-recipe/build.sh new file mode 100644 index 0000000..5dd97b0 --- /dev/null +++ b/Wrappers/Python/conda-recipe/build.sh @@ -0,0 +1,9 @@ +if [ -z "$CIL_VERSION" ]; then + echo "Need to set CIL_VERSION" + exit 1 +fi +mkdir ${SRC_DIR}/ccpi +cp -r "${RECIPE_DIR}/../../../" ${SRC_DIR}/ccpi + +cd ${SRC_DIR}/ccpi/Wrappers/Python +$PYTHON setup.py install diff --git a/Wrappers/Python/conda-recipe/meta.yaml b/Wrappers/Python/conda-recipe/meta.yaml new file mode 100644 index 0000000..f97de5e --- /dev/null +++ b/Wrappers/Python/conda-recipe/meta.yaml @@ -0,0 +1,27 @@ +package: + name: ccpi-plugins + version: {{ environ['CIL_VERSION'] }} + + +build: + preserve_egg_dir: False + script_env: + - CIL_VERSION +# number: 0 + +requirements: + build: + - python + - setuptools + + run: + - python + - numpy + - ccpi-framework + - ccpi-reconstruction + - matplotlib + +about: + home: http://www.ccpi.ac.uk + license: Apache 2.0 License + summary: 'CCPi Framework Plugins' diff --git a/Wrappers/Python/setup.py b/Wrappers/Python/setup.py new file mode 100644 index 0000000..f4c42e8 --- /dev/null +++ b/Wrappers/Python/setup.py @@ -0,0 +1,58 @@ +#!/usr/bin/env python +# -*- coding: utf-8 -*- +# This work is part of the Core Imaging Library developed by +# Visual Analytics and Imaging System Group of the Science Technology +# Facilities Council, STFC + +# Copyright 2018 Edoardo Pasca + +# 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 + +# 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 distutils.core import setup +import os +import sys + + + +cil_version=os.environ['CIL_VERSION'] +if cil_version == '': + print("Please set the environmental variable CIL_VERSION") + sys.exit(1) + +setup( + name="ccpi-plugins", + version=cil_version, + packages=['ccpi' , 'ccpi.plugins'], + + # Project uses reStructuredText, so ensure that the docutils get + # installed or upgraded on the target machine + #install_requires=['docutils>=0.3'], + +# package_data={ +# # If any package contains *.txt or *.rst files, include them: +# '': ['*.txt', '*.rst'], +# # And include any *.msg files found in the 'hello' package, too: +# 'hello': ['*.msg'], +# }, + # zip_safe = False, + + # metadata for upload to PyPI + author="Edoardo Pasca", + author_email="edoardo.pasca@stfc.ac.uk", + description='CCPi Core Imaging Library - Python Framework Plugins Module', + license="Apache v2.0", + keywords="Python Framework Plugins", + url="http://www.ccpi.ac.uk", # project home page, if any + + # could also include long_description, download_url, classifiers, etc. +) diff --git a/Wrappers/Python/wip/simple_demo_ccpi.py b/Wrappers/Python/wip/simple_demo_ccpi.py new file mode 100755 index 0000000..10bb75f --- /dev/null +++ b/Wrappers/Python/wip/simple_demo_ccpi.py @@ -0,0 +1,272 @@ +#import sys +#sys.path.append("..") + +from ccpi.framework import ImageData , AcquisitionData, ImageGeometry, AcquisitionGeometry + +from ccpi.optimisation.algs import FISTA, FBPD, CGLS +from ccpi.optimisation.funcs import Norm2sq, Norm1 , TV2D + +from ccpi.plugins.ops import CCPiProjectorSimple +from ccpi.plugins.processors import CCPiForwardProjector, CCPiBackwardProjector + +from ccpi.reconstruction.parallelbeam import alg as pbalg + +import numpy as np +import matplotlib.pyplot as plt + +test_case = 1 # 1=parallel2D, 2=cone2D, 3=parallel3D + +# Set up phantom +N = 128 +vert = 1 +# Set up measurement geometry +angles_num = 20; # angles number +det_w = 1.0 +det_num = N +SourceOrig = 200 +OrigDetec = 0 + +if test_case==1: + angles = np.linspace(0,np.pi,angles_num,endpoint=False,dtype=np.float32)*180/np.pi + #nangles = angles_num + #angles = np.linspace(0,360, nangles, dtype=np.float32) + +elif test_case==2: + angles = np.linspace(0,2*np.pi,angles_num,endpoint=False) +elif test_case == 3: + angles = np.linspace(0,np.pi,angles_num,endpoint=False) +else: + NotImplemented + + +def setupCCPiGeometries(voxel_num_x, voxel_num_y, voxel_num_z, angles, counter): + vg = ImageGeometry(voxel_num_x=voxel_num_x,voxel_num_y=voxel_num_y, voxel_num_z=voxel_num_z) + Phantom_ccpi = ImageData(geometry=vg, + dimension_labels=['horizontal_x','horizontal_y','vertical']) + #.subset(['horizontal_x','horizontal_y','vertical']) + # ask the ccpi code what dimensions it would like + + voxel_per_pixel = 1 + geoms = pbalg.pb_setup_geometry_from_image(Phantom_ccpi.as_array(), + angles, + voxel_per_pixel ) + + pg = AcquisitionGeometry('parallel', + '3D', + angles, + geoms['n_h'],det_w, + geoms['n_v'], det_w #2D in 3D is a slice 1 pixel thick + ) + + center_of_rotation = Phantom_ccpi.get_dimension_size('horizontal_x') / 2 + ad = AcquisitionData(geometry=pg,dimension_labels=['angle','vertical','horizontal']) + geoms_i = pbalg.pb_setup_geometry_from_acquisition(ad.as_array(), + angles, + center_of_rotation, + voxel_per_pixel ) + + #print (counter) + counter+=1 + #print (geoms , geoms_i) + if counter < 4: + if (not ( geoms_i == geoms )): + print ("not equal and {0}".format(counter)) + X = max(geoms['output_volume_x'], geoms_i['output_volume_x']) + Y = max(geoms['output_volume_y'], geoms_i['output_volume_y']) + Z = max(geoms['output_volume_z'], geoms_i['output_volume_z']) + return setupCCPiGeometries(X,Y,Z,angles, counter) + else: + print ("return geoms {0}".format(geoms)) + return geoms + else: + print ("return geoms_i {0}".format(geoms_i)) + return geoms_i + +geoms = setupCCPiGeometries(N,N,vert,angles,0) +#%% +#geoms = {'n_v': 12, 'output_volume_y': 128, 'n_angles': 20, +# 'output_volume_x': 128, 'output_volume_z': 12, 'n_h': 128} +vg = ImageGeometry(voxel_num_x=geoms['output_volume_x'], + voxel_num_y=geoms['output_volume_y'], + voxel_num_z=geoms['output_volume_z']) +Phantom = ImageData(geometry=vg,dimension_labels=['horizontal_x','horizontal_y','vertical']) + 0.1 + + +#x = Phantom.as_array() +i = 0 +while i < geoms['n_v']: + #x = Phantom.subset(vertical=i, dimensions=['horizontal_x','horizontal_y']).array + x = Phantom.subset(vertical=i).array + x[round(N/4):round(3*N/4),round(N/4):round(3*N/4)] = 0.5 + x[round(N/8):round(7*N/8),round(3*N/8):round(5*N/8)] = 0.98 + Phantom.fill(x, vertical=i) + i += 1 + +plt.imshow(Phantom.subset(vertical=0).as_array()) +plt.show() + + + +# Parallelbeam geometry test +if test_case==1: + #Phantom_ccpi = Phantom.subset(['horizontal_x','horizontal_y','vertical']) + #Phantom_ccpi.geometry = vg.clone() + center_of_rotation = Phantom.get_dimension_size('horizontal_x') / 2 + + pg = AcquisitionGeometry('parallel', + '3D', + angles, + geoms['n_h'],det_w, + geoms['n_v'], det_w #2D in 3D is a slice 1 pixel thick + ) +elif test_case==2: + raise NotImplemented('cone beam projector not yet available') + pg = AcquisitionGeometry('cone', + '2D', + angles, + det_num, + det_w, + vert, det_w, #2D in 3D is a slice 1 pixel thick + dist_source_center=SourceOrig, + dist_center_detector=OrigDetec) + +# ASTRA operator using volume and sinogram geometries +#Aop = AstraProjectorSimple(vg, pg, 'cpu') +Cop = CCPiProjectorSimple(vg, pg) + +# Try forward and backprojection +b = Cop.direct(Phantom) +out2 = Cop.adjoint(b) + +#%% +for i in range(b.get_dimension_size('vertical')): + plt.imshow(b.subset(vertical=i).array) + #plt.imshow(Phantom.subset( vertical=i).array) + #plt.imshow(b.array[:,i,:]) + plt.show() +#%% + +plt.imshow(out2.subset( vertical=0).array) +plt.show() + +# Create least squares object instance with projector and data. +f = Norm2sq(Cop,b,c=0.5) + +# Initial guess +x_init = ImageData(geometry=vg, dimension_labels=['horizontal_x','horizontal_y','vertical']) +#invL = 0.5 +#g = f.grad(x_init) +#print (g) +#u = x_init - invL*f.grad(x_init) + +#%% +# Run FISTA for least squares without regularization +opt = {'tol': 1e-4, 'iter': 100} +x_fista0, it0, timing0, criter0 = FISTA(x_init, f, None, opt=opt) + +plt.imshow(x_fista0.subset(vertical=0).array) +plt.title('FISTA0') +plt.show() + +# Now least squares plus 1-norm regularization +lam = 0.1 +g0 = Norm1(lam) + +# Run FISTA for least squares plus 1-norm function. +x_fista1, it1, timing1, criter1 = FISTA(x_init, f, g0,opt=opt) + +plt.imshow(x_fista0.subset(vertical=0).array) +plt.title('FISTA1') +plt.show() + +plt.semilogy(criter1) +plt.show() + +# Run FBPD=Forward Backward Primal Dual method on least squares plus 1-norm +x_fbpd1, it_fbpd1, timing_fbpd1, criter_fbpd1 = FBPD(x_init,None,f,g0,opt=opt) + +plt.imshow(x_fbpd1.subset(vertical=0).array) +plt.title('FBPD1') +plt.show() + +plt.semilogy(criter_fbpd1) +plt.show() + +# Now FBPD for least squares plus TV +#lamtv = 1 +#gtv = TV2D(lamtv) + +#x_fbpdtv, it_fbpdtv, timing_fbpdtv, criter_fbpdtv = FBPD(x_init,None,f,gtv,opt=opt) + +#plt.imshow(x_fbpdtv.subset(vertical=0).array) +#plt.show() + +#plt.semilogy(criter_fbpdtv) +#plt.show() + + +# Run CGLS, which should agree with the FISTA0 +x_CGLS, it_CGLS, timing_CGLS, criter_CGLS = CGLS(x_init, Cop, b, opt=opt) + +plt.imshow(x_CGLS.subset(vertical=0).array) +plt.title('CGLS') +plt.title('CGLS recon, compare FISTA0') +plt.show() + +plt.semilogy(criter_CGLS) +plt.title('CGLS criterion') +plt.show() + + +#%% + +clims = (0,1) +cols = 3 +rows = 2 +current = 1 +fig = plt.figure() +# projections row +a=fig.add_subplot(rows,cols,current) +a.set_title('phantom {0}'.format(np.shape(Phantom.as_array()))) + +imgplot = plt.imshow(Phantom.subset(vertical=0).as_array(),vmin=clims[0],vmax=clims[1]) + +current = current + 1 +a=fig.add_subplot(rows,cols,current) +a.set_title('FISTA0') +imgplot = plt.imshow(x_fista0.subset(vertical=0).as_array(),vmin=clims[0],vmax=clims[1]) + +current = current + 1 +a=fig.add_subplot(rows,cols,current) +a.set_title('FISTA1') +imgplot = plt.imshow(x_fista1.subset(vertical=0).as_array(),vmin=clims[0],vmax=clims[1]) + +current = current + 1 +a=fig.add_subplot(rows,cols,current) +a.set_title('FBPD1') +imgplot = plt.imshow(x_fbpd1.subset(vertical=0).as_array(),vmin=clims[0],vmax=clims[1]) + +current = current + 1 +a=fig.add_subplot(rows,cols,current) +a.set_title('CGLS') +imgplot = plt.imshow(x_CGLS.subset(vertical=0).as_array(),vmin=clims[0],vmax=clims[1]) + +plt.show() +#%% +#current = current + 1 +#a=fig.add_subplot(rows,cols,current) +#a.set_title('FBPD TV') +#imgplot = plt.imshow(x_fbpdtv.subset(vertical=0).as_array(),vmin=clims[0],vmax=clims[1]) + +fig = plt.figure() +# projections row +b=fig.add_subplot(1,1,1) +b.set_title('criteria') +imgplot = plt.loglog(criter0 , label='FISTA0') +imgplot = plt.loglog(criter1 , label='FISTA1') +imgplot = plt.loglog(criter_fbpd1, label='FBPD1') +imgplot = plt.loglog(criter_CGLS, label='CGLS') +#imgplot = plt.loglog(criter_fbpdtv, label='FBPD TV') +b.legend(loc='right') +plt.show() +#%%
\ No newline at end of file |