From ee0a3bf634e92351ad82c09b4e94bcb64ee52d73 Mon Sep 17 00:00:00 2001 From: Edoardo Pasca Date: Wed, 11 Apr 2018 16:18:20 +0100 Subject: initial working release with simple_demo_ccpi.py --- Wrappers/Python/ccpi/plugins/__init__.py | 0 Wrappers/Python/ccpi/plugins/ops.py | 102 ++++ Wrappers/Python/ccpi/plugins/processors.py | 792 +++++++++++++++++++++++++++++ Wrappers/Python/conda-recipe/bld.bat | 10 + Wrappers/Python/conda-recipe/build.sh | 9 + Wrappers/Python/conda-recipe/meta.yaml | 27 + Wrappers/Python/setup.py | 58 +++ Wrappers/Python/wip/simple_demo_ccpi.py | 272 ++++++++++ 8 files changed, 1270 insertions(+) create mode 100644 Wrappers/Python/ccpi/plugins/__init__.py create mode 100755 Wrappers/Python/ccpi/plugins/ops.py create mode 100755 Wrappers/Python/ccpi/plugins/processors.py create mode 100644 Wrappers/Python/conda-recipe/bld.bat create mode 100644 Wrappers/Python/conda-recipe/build.sh create mode 100644 Wrappers/Python/conda-recipe/meta.yaml create mode 100644 Wrappers/Python/setup.py create mode 100755 Wrappers/Python/wip/simple_demo_ccpi.py (limited to 'Wrappers/Python') diff --git a/Wrappers/Python/ccpi/plugins/__init__.py b/Wrappers/Python/ccpi/plugins/__init__.py new file mode 100644 index 0000000..e69de29 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(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 -- cgit v1.2.3