Merge remote-tracking branch 'origin' into rename_as_ccppetmr

3 files changed, 505 insertions, 578 deletions

diff --git a/Wrappers/Python/ccpi/processors.py b/Wrappers/Python/ccpi/processors.py
index 9138a27..87acf92 100755
--- a/Wrappers/Python/ccpi/processors.py
+++ b/Wrappers/Python/ccpi/processors.py
@@ -1,470 +1,470 @@
-# -*- 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 DataSetProcessor, DataSet, VolumeData, SinogramData, ImageGeometry, AcquisitionGeometry
-import numpy
-import h5py
-from scipy import ndimage
-
-class Normalizer(DataSetProcessor):
-    '''Normalization based on flat and dark
-    
-    This processor read in a SinogramDataSet and normalises it based on 
-    the instrument reading with and without incident photons or neutrons.
-    
-    Input: SinogramDataSet
-    Parameter: 2D projection with flat field (or stack)
-               2D projection with dark field (or stack)
-    Output: SinogramDataSetn
-    '''
-    
-    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
-                  }
-        
-        #DataSetProcessor.__init__(self, **kwargs)
-        super(Normalizer, self).__init__(**kwargs)
-        if not flat_field is None:
-            self.setFlatField(flat_field)
-        if not dark_field is None:
-            self.setDarkField(dark_field)
-    
-    def checkInput(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 setDarkField(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.setDarkField(df.as_array())
-    
-    def setFlatField(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.setDarkField(df.as_array())
-    
-    @staticmethod
-    def normalizeProjection(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.getInput()
-        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.normalizeProjection(
-                        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(DataSetProcessor):
-    '''Processor to find the center of rotation in a parallel beam experiment
-    
-    This processor read in a SinogramDataSet and finds the center of rotation 
-    based on Nghia Vo's method. https://doi.org/10.1364/OE.22.019078
-    
-    Input: SinogramDataSet
-    
-    Output: float. center of rotation in pixel coordinate
-    '''
-    
-    def __init__(self):
-        kwargs = {
-                  
-                  }
-        
-        #DataSetProcessor.__init__(self, **kwargs)
-        super(CenterOfRotationFinder, self).__init__(**kwargs)
-    
-    def checkInput(self, dataset):
-        if dataset.number_of_dimensions == 3:
-            if dataset.geometry.geom_type == 'parallel':
-                return True
-            else:
-                raise ValueError('This algorithm is suitable only for parallel beam geometry')
-        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.getInput()
-        
-        cor = CenterOfRotationFinder.find_center_vo(projections.as_array())
-        
-        return cor
-
-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],
-                                 )    
-    sino = SinogramData( proj , geometry=parallelbeam)
-    
-    normalizer = Normalizer()
-    normalizer.setInput(sino)
-    normalizer.setFlatField(flat)
-    normalizer.setDarkField(dark)
-    norm = normalizer.getOutput()
-    print ("Processor min {0} max {1}".format(norm.as_array().min(), norm.as_array().max()))
-    
-    norm1 = numpy.asarray(
-            [Normalizer.normalizeProjection( 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.setInput(sino)
-    cor = cor_finder.getOutput()
-    print ("center of rotation {0} == 86.25?".format(cor))
-    
-    conebeam = AcquisitionGeometry('cone', '3D' , 
-                                 angles=angles, 
-                                 pixel_num_h=numpy.shape(proj)[2],
-                                 pixel_num_v=numpy.shape(proj)[1],
-                                 )    
-    sino = SinogramData( proj , geometry=conebeam)
-    try:
-        cor_finder.setInput(sino)
-        cor = cor_finder.getOutput()
-    except ValueError as err:
+# -*- 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 DataSetProcessor, DataSet, VolumeData, SinogramData, ImageGeometry, AcquisitionGeometry

+import numpy

+import h5py

+from scipy import ndimage

+

+class Normalizer(DataSetProcessor):

+    '''Normalization based on flat and dark

+    

+    This processor read in a SinogramDataSet and normalises it based on 

+    the instrument reading with and without incident photons or neutrons.

+    

+    Input: SinogramDataSet

+    Parameter: 2D projection with flat field (or stack)

+               2D projection with dark field (or stack)

+    Output: SinogramDataSetn

+    '''

+    

+    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

+                  }

+        

+        #DataSetProcessor.__init__(self, **kwargs)

+        super(Normalizer, self).__init__(**kwargs)

+        if not flat_field is None:

+            self.setFlatField(flat_field)

+        if not dark_field is None:

+            self.setDarkField(dark_field)

+    

+    def checkInput(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 setDarkField(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.setDarkField(df.as_array())

+    

+    def setFlatField(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.setDarkField(df.as_array())

+    

+    @staticmethod

+    def normalizeProjection(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.getInput()

+        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.normalizeProjection(

+                        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(DataSetProcessor):

+    '''Processor to find the center of rotation in a parallel beam experiment

+    

+    This processor read in a SinogramDataSet and finds the center of rotation 

+    based on Nghia Vo's method. https://doi.org/10.1364/OE.22.019078

+    

+    Input: SinogramDataSet

+    

+    Output: float. center of rotation in pixel coordinate

+    '''

+    

+    def __init__(self):

+        kwargs = {

+                  

+                  }

+        

+        #DataSetProcessor.__init__(self, **kwargs)

+        super(CenterOfRotationFinder, self).__init__(**kwargs)

+    

+    def checkInput(self, dataset):

+        if dataset.number_of_dimensions == 3:

+            if dataset.geometry.geom_type == 'parallel':

+                return True

+            else:

+                raise ValueError('This algorithm is suitable only for parallel beam geometry')

+        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.getInput()

+        

+        cor = CenterOfRotationFinder.find_center_vo(projections.as_array())

+        

+        return cor

+

+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],

+                                 )    

+    sino = SinogramData( proj , geometry=parallelbeam)

+    

+    normalizer = Normalizer()

+    normalizer.setInput(sino)

+    normalizer.setFlatField(flat)

+    normalizer.setDarkField(dark)

+    norm = normalizer.getOutput()

+    print ("Processor min {0} max {1}".format(norm.as_array().min(), norm.as_array().max()))

+    

+    norm1 = numpy.asarray(

+            [Normalizer.normalizeProjection( 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.setInput(sino)

+    cor = cor_finder.getOutput()

+    print ("center of rotation {0} == 86.25?".format(cor))

+    

+    conebeam = AcquisitionGeometry('cone', '3D' , 

+                                 angles=angles, 

+                                 pixel_num_h=numpy.shape(proj)[2],

+                                 pixel_num_v=numpy.shape(proj)[1],

+                                 )    

+    sino = SinogramData( proj , geometry=conebeam)

+    try:

+        cor_finder.setInput(sino)

+        cor = cor_finder.getOutput()

+    except ValueError as err:

         print (err)
\ No newline at end of file
diff --git a/Wrappers/Python/wip/simple_demo.py b/Wrappers/Python/wip/simple_demo.py
index 629d9b5..655b68d 100644
--- a/Wrappers/Python/wip/simple_demo.py
+++ b/Wrappers/Python/wip/simple_demo.py
@@ -95,7 +95,7 @@ g0 = Norm1(lam)
 x_fista1, it1, timing1, criter1 = FISTA(x_init, f, g0)
 
 plt.imshow(x_fista1.array)
-plt.title('FISTA')
+plt.title('FISTA1')
 plt.show()
 
 plt.semilogy(criter1)
@@ -106,12 +106,25 @@ opt = {'tol': 1e-4, 'iter': 10000}
 x_fbpd1, it_fbpd1, timing_fbpd1, criter_fbpd1 = FBPD(x_init,None,f,g0,opt=opt)
 
 plt.imshow(x_fbpd1.array)
-plt.title('FBPD')
+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.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(Aop, b, 1000, x_init)
 
@@ -126,6 +139,8 @@ plt.show()
 
 
 #%%
+
+clims = (-0.5,2.5)
 cols = 3
 rows = 2
 current = 1
@@ -133,35 +148,46 @@ fig = plt.figure()
 # projections row
 a=fig.add_subplot(rows,cols,current)
 a.set_title('phantom {0}'.format(np.shape(Phantom.as_array())))
+<<<<<<< HEAD
 imgplot = plt.imshow(Phantom.as_array())
+=======
+imgplot = plt.imshow(Phantom.as_array(),vmin=clims[0],vmax=clims[1])
+>>>>>>> origin
 
 current = current + 1
 a=fig.add_subplot(rows,cols,current)
 a.set_title('FISTA0')
-imgplot = plt.imshow(x_fista0.as_array())
+imgplot = plt.imshow(x_fista0.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.as_array())
+imgplot = plt.imshow(x_fista1.as_array(),vmin=clims[0],vmax=clims[1])
 
 current = current + 1
 a=fig.add_subplot(rows,cols,current)
-a.set_title('FBPD')
-imgplot = plt.imshow(x_fbpd1.as_array())
+a.set_title('FBPD1')
+imgplot = plt.imshow(x_fbpd1.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.as_array())
+imgplot = plt.imshow(x_CGLS.as_array(),vmin=clims[0],vmax=clims[1])
 
 current = current + 1
 a=fig.add_subplot(rows,cols,current)
-a.set_title('criteria')
+a.set_title('FBPD TV')
+imgplot = plt.imshow(x_fbpdtv.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='FBPD')
+imgplot = plt.loglog(criter_fbpd1, label='FBPD1')
 imgplot = plt.loglog(criter_CGLS, label='CGLS')
-a.legend(loc='right')
+imgplot = plt.loglog(criter_fbpdtv, label='FBPD TV')
+b.legend(loc='right')
 plt.show()
 #%%
\ No newline at end of file
diff --git a/Wrappers/Python/wip/simple_demo_tv.py b/Wrappers/Python/wip/simple_demo_tv.py
deleted file mode 100644
index 29e2349..0000000
--- a/Wrappers/Python/wip/simple_demo_tv.py
+++ /dev/null
@@ -1,99 +0,0 @@
-
-
-import sys
-
-sys.path.append("..")
-
-from ccpi.framework import *
-from ccpi.reconstruction.algs import *
-from ccpi.reconstruction.funcs import *
-from ccpi.reconstruction.ops import *
-from ccpi.reconstruction.astra_ops import *
-
-import numpy as np
-import matplotlib.pyplot as plt
-
-test_case = 2   # 1=parallel2D, 2=cone2D
-
-# Set up phantom
-N = 128
-
-x = np.zeros((N,N))
-x[round(N/4):round(3*N/4),round(N/4):round(3*N/4)] = 1.0
-x[round(N/8):round(7*N/8),round(3*N/8):round(5*N/8)] = 2.0
-
-plt.imshow(x)
-plt.show()
-
-vg = ImageGeometry(N,N,None, 1,1,None)
-
-Phantom = VolumeData(x,geometry=vg)
-
-# Set up measurement geometry
-angles_num = 20; # angles number
-
-if test_case==1:
-    angles = np.linspace(0,np.pi,angles_num,endpoint=False)
-elif test_case==2:
-    angles = np.linspace(0,2*np.pi,angles_num,endpoint=False)
-else:
-    NotImplemented
-
-det_w = 1.0
-det_num = N
-SourceOrig = 200
-OrigDetec = 0
-
-# Parallelbeam geometry test
-if test_case==1:
-    pg = AcquisitionGeometry('parallel',
-                          '2D',
-                          angles,
-                          det_num,det_w)
-elif test_case==2:
-    pg = AcquisitionGeometry('cone',
-                          '2D',
-                          angles,
-                          det_num,
-                          det_w,
-                          dist_source_center=SourceOrig, 
-                          dist_center_detector=OrigDetec)
-
-# ASTRA operator using volume and sinogram geometries
-Aop = AstraProjectorSimple(vg, pg, 'gpu')
-
-# Unused old astra projector without geometry
-# Aop_old = AstraProjector(det_w, det_num, SourceOrig, 
-#                      OrigDetec, angles, 
-#                      N,'fanbeam','gpu') 
-
-# Try forward and backprojection
-b = Aop.direct(Phantom)
-out2 = Aop.adjoint(b)
-
-plt.imshow(b.array)
-plt.show()
-
-plt.imshow(out2.array)
-plt.show()
-
-# Create least squares object instance with projector and data.
-f = Norm2sq(Aop,b,c=0.5)
-
-# Initial guess
-x_init = VolumeData(np.zeros(x.shape),geometry=vg)
-
-# Now least squares plus 1-norm regularization
-lam = 1
-g0 = TV2D(lam)
-
-
-# Run FBPD=Forward Backward Primal Dual method on least squares plus 1-norm
-opt = {'tol': 1e-4, 'iter': 10000}
-x_fbpd1, it_fbpd1, timing_fbpd1, criter_fbpd1 = FBPD(x_init,None,f,g0,opt=opt)
-
-plt.imshow(x_fbpd1.array)
-plt.show()
-
-plt.semilogy(criter_fbpd1)
-plt.show()
\ No newline at end of file