Nexus read slice (#144)

* replaces in with explicit test

closes #139

* add read of single slice or subset

* get_image_keys

* reader reads single slice

* adds get_acquisition_data_slice and get_acquisition_data_subset

it is possible to read in only a subset of the dataset.

* added multifile_nexus.py

* cleaned example

* fix normaliser and reader

* wip

* wip for reader

* added support for non standard?? nexus files

updated example with normalisation.

* fix camelcase

* Added initial unittest for reader

bugfix for python 2.7

* fixes for python 2.7

* remove duplicate test requirement

* minimal unittest for get_acquisition_data_subset

7 files changed, 762 insertions, 71 deletions

diff --git a/Wrappers/Python/ccpi/framework.py b/Wrappers/Python/ccpi/framework.py
index 4d74d2b..9938fb7 100644
--- a/Wrappers/Python/ccpi/framework.py
+++ b/Wrappers/Python/ccpi/framework.py
@@ -17,18 +17,16 @@
 #   See the License for the specific language governing permissions and
 #   limitations under the License.
 
+from __future__ import absolute_import
 from __future__ import division
-import abc
+from __future__ import print_function
+from __future__ import unicode_literals
+
 import numpy
 import sys
 from datetime import timedelta, datetime
 import warnings
 
-if sys.version_info[0] >= 3 and sys.version_info[1] >= 4:
-    ABC = abc.ABC
-else:
-    ABC = abc.ABCMeta('ABC', (), {})
-
 def find_key(dic, val):
     """return the key of dictionary dic given the value"""
     return [k for k, v in dic.items() if v == val][0]
@@ -376,7 +374,6 @@ class DataContainer(object):
         return self.shape == other.shape
     
     ## algebra 
-    
     def __add__(self, other , out=None, *args, **kwargs):
         if issubclass(type(other), DataContainer):    
             if self.check_dimensions(other):
@@ -525,9 +522,9 @@ class DataContainer(object):
     # must return self
     
     
+    
     def __iadd__(self, other):
         if isinstance(other, (int, float)) :
-            #print ("__iadd__", self.array.shape)
             numpy.add(self.array, other, out=self.array)
         elif issubclass(type(other), DataContainer):
             if self.check_dimensions(other):
@@ -539,11 +536,7 @@ class DataContainer(object):
     
     def __imul__(self, other):
         if isinstance(other, (int, float)) :
-            #print ("__imul__", self.array.shape)
-            #print ("type(self)", type(self))
-            #print ("self.array", self.array, other)
             arr = self.as_array()
-            #print ("arr", arr)
             numpy.multiply(arr, other, out=arr)
         elif issubclass(type(other), DataContainer):
             if self.check_dimensions(other):
@@ -627,13 +620,14 @@ class DataContainer(object):
                                  .format(len(self.shape),len(new_order)))
         
                 
-    def copy(self): 
+    def copy(self):
         '''alias of clone'''    
         return self.clone()
     
     ## binary operations
             
     def pixel_wise_binary(self,pwop, x2 , out=None, *args,  **kwargs):    
+
         if out is None:
             if isinstance(x2, (int, float, complex)):
                 out = pwop(self.as_array() , x2 , *args, **kwargs )
@@ -657,11 +651,8 @@ class DataContainer(object):
                 raise ValueError(message(type(self),"Wrong size for data memory: ", out.shape,self.shape))
         elif issubclass(type(out), DataContainer) and isinstance(x2, (int,float,complex)):
             if self.check_dimensions(out):
+
                 pwop(self.as_array(), x2, out=out.as_array(), *args, **kwargs )
-                #return type(self)(out.as_array(),
-                #       deep_copy=False, 
-                #       dimension_labels=self.dimension_labels,
-                #       geometry=self.geometry)
                 return out
             else:
                 raise ValueError(message(type(self),"Wrong size for data memory: ", out.shape,self.shape))
@@ -694,7 +685,6 @@ class DataContainer(object):
         return self.pixel_wise_binary(numpy.maximum, x2=x2, out=out, *args, **kwargs)
     
     ## unary operations
-    
     def pixel_wise_unary(self,pwop, out=None, *args,  **kwargs):
         if out is None:
             out = pwop(self.as_array() , *args, **kwargs )
@@ -705,19 +695,11 @@ class DataContainer(object):
         elif issubclass(type(out), DataContainer):
             if self.check_dimensions(out):
                 pwop(self.as_array(), out=out.as_array(), *args, **kwargs )
-                #return type(self)(out.as_array(),
-                #       deep_copy=False, 
-                #       dimension_labels=self.dimension_labels,
-                #       geometry=self.geometry)
             else:
                 raise ValueError(message(type(self),"Wrong size for data memory: ", out.shape,self.shape))
         elif issubclass(type(out), numpy.ndarray):
             if self.array.shape == out.shape and self.array.dtype == out.dtype:
                 pwop(self.as_array(), out=out, *args, **kwargs)
-                #return type(self)(out,
-                #       deep_copy=False, 
-                #       dimension_labels=self.dimension_labels,
-                #       geometry=self.geometry)
         else:
             raise ValueError (message(type(self),  "incompatible class:" , pwop.__name__, type(out)))
     
@@ -725,11 +707,6 @@ class DataContainer(object):
         return self.pixel_wise_unary(numpy.abs, out=out, *args,  **kwargs)
     
     def sign(self, out=None, *args,  **kwargs):
-        #out = numpy.sign(self.as_array() )
-        #return type(self)(out,
-        #               deep_copy=False, 
-        #               dimension_labels=self.dimension_labels,
-        #               geometry=self.geometry)
         return self.pixel_wise_unary(numpy.sign , out=out, *args,  **kwargs)
     
     def sqrt(self, out=None, *args,  **kwargs):
@@ -743,7 +720,6 @@ class DataContainer(object):
     ## reductions
     def sum(self, out=None, *args, **kwargs):
         return self.as_array().sum(*args, **kwargs)
-        #return self.pixel_wise_unary(numpy.sum, out=out, *args, **kwargs)
     
 class ImageData(DataContainer):
     '''DataContainer for holding 2D or 3D DataContainer'''
diff --git a/Wrappers/Python/ccpi/io/reader.py b/Wrappers/Python/ccpi/io/reader.py
index b0b5414..856f5e0 100644
--- a/Wrappers/Python/ccpi/io/reader.py
+++ b/Wrappers/Python/ccpi/io/reader.py
@@ -22,6 +22,11 @@ This is a reader module with classes for loading 3D datasets.
 

 @author: Mr. Srikanth Nagella

 '''

+from __future__ import absolute_import

+from __future__ import division

+from __future__ import print_function

+from __future__ import unicode_literals

+

 from ccpi.framework import AcquisitionGeometry

 from ccpi.framework import AcquisitionData

 import numpy as np

@@ -53,7 +58,18 @@ class NexusReader(object):
         self.angles = None

         self.geometry = None

         self.filename = nexus_filename

-        

+        self.key_path = 'entry1/tomo_entry/instrument/detector/image_key'

+        self.data_path = 'entry1/tomo_entry/data/data'

+        self.angle_path = 'entry1/tomo_entry/data/rotation_angle'

+    

+    def get_image_keys(self):

+        try:

+            with NexusFile(self.filename,'r') as file:    

+                return np.array(file[self.key_path])

+        except KeyError as ke:

+            raise KeyError("get_image_keys: " , ke.args[0] , self.key_path)

+            

+    

     def load(self, dimensions=None, image_key_id=0):  

         '''

         This is generic loading function of flat field, dark field and projection data.

@@ -64,10 +80,10 @@ class NexusReader(object):
             return        

         try:

             with NexusFile(self.filename,'r') as file:                

-                image_keys = np.array(file['entry1/tomo_entry/instrument/detector/image_key'])                

+                image_keys = np.array(file[self.key_path])                

                 projections = None

                 if dimensions == None:

-                    projections = np.array(file['entry1/tomo_entry/data/data'])

+                    projections = np.array(file[self.data_path])

                     result = projections[image_keys==image_key_id]

                     return result

                 else:

@@ -76,8 +92,8 @@ class NexusReader(object):
                     projection_indexes = index_array[0][dimensions[0]]

                     new_dimensions = list(dimensions)

                     new_dimensions[0]= projection_indexes

-                    new_dimensions = tuple(new_dimensions)              

-                    result = np.array(file['entry1/tomo_entry/data/data'][new_dimensions])

+                    new_dimensions = tuple(new_dimensions)

+                    result = np.array(file[self.data_path][new_dimensions])

                     return result

         except:

             print("Error reading nexus file")

@@ -88,6 +104,12 @@ class NexusReader(object):
         Loads the projection data from the nexus file.

         returns: numpy array with projection data

         '''

+        try:

+            if 0 not in self.get_image_keys():

+                raise ValueError("Projections are not in the data. Data Path " , 

+                                 self.data_path)

+        except KeyError as ke:

+            raise KeyError(ke.args[0] , self.data_path)

         return self.load(dimensions, 0)

     

     def load_flat(self, dimensions=None):

@@ -95,6 +117,12 @@ class NexusReader(object):
         Loads the flat field data from the nexus file.

         returns: numpy array with flat field data

         '''        

+        try:

+            if 1 not in self.get_image_keys():

+                raise ValueError("Flats are not in the data. Data Path " , 

+                                 self.data_path)

+        except KeyError as ke:

+            raise KeyError(ke.args[0] , self.data_path)

         return self.load(dimensions, 1)

     

     def load_dark(self, dimensions=None):

@@ -102,6 +130,12 @@ class NexusReader(object):
         Loads the Dark field data from the nexus file.

         returns: numpy array with dark field data

         '''        

+        try:

+            if 2 not in self.get_image_keys():

+                raise ValueError("Darks are not in the data. Data Path " , 

+                             self.data_path)

+        except KeyError as ke:

+            raise KeyError(ke.args[0] , self.data_path)

         return self.load(dimensions, 2)

     

     def get_projection_angles(self):

@@ -114,11 +148,11 @@ class NexusReader(object):
             return        

         try:

             with NexusFile(self.filename,'r') as file:                

-                angles = np.array(file['entry1/tomo_entry/data/rotation_angle'],np.float32)

-                image_keys = np.array(file['entry1/tomo_entry/instrument/detector/image_key'])                

+                angles = np.array(file[self.angle_path],np.float32)

+                image_keys = np.array(file[self.key_path])                

                 return angles[image_keys==0]

         except:

-            print("Error reading nexus file")

+            print("get_projection_angles Error reading nexus file")

             raise        

 

     

@@ -132,8 +166,8 @@ class NexusReader(object):
             return

         try:

             with NexusFile(self.filename,'r') as file:                

-                projections = file['entry1/tomo_entry/data/data']

-                image_keys = np.array(file['entry1/tomo_entry/instrument/detector/image_key'])

+                projections = file[self.data_path]

+                image_keys = np.array(file[self.key_path])

                 dims = list(projections.shape)

                 dims[0] = dims[1]

                 dims[1] = np.sum(image_keys==0)

@@ -151,15 +185,25 @@ class NexusReader(object):
         if self.filename is None:

             return

         try:

-            with NexusFile(self.filename,'r') as file:                

-                projections = file['entry1/tomo_entry/data/data']

-                image_keys = np.array(file['entry1/tomo_entry/instrument/detector/image_key'])

+            with NexusFile(self.filename,'r') as file:

+                try:                

+                    projections = file[self.data_path]

+                except KeyError as ke:

+                    raise KeyError('Error: data path {0} not found\n{1}'\

+                                   .format(self.data_path, 

+                                           ke.args[0]))

+                #image_keys = np.array(file[self.key_path])

+                image_keys = self.get_image_keys()

                 dims = list(projections.shape)

                 dims[0] = np.sum(image_keys==0)

                 return tuple(dims)

         except:

-            print("Error reading nexus file")

-            raise  

+            print("Warning: Error reading image_keys trying accessing data on " , self.data_path)

+            with NexusFile(self.filename,'r') as file:

+                dims = file[self.data_path].shape

+                return tuple(dims)

+            

+              

         

     def get_acquisition_data(self, dimensions=None):

         '''

@@ -179,6 +223,160 @@ class NexusReader(object):
         return AcquisitionData(data, geometry=geometry, 

                                dimension_labels=['angle','vertical','horizontal'])   

     

+    def get_acquisition_data_subset(self, ymin=None, ymax=None):

+        '''

+        This method load the acquisition data and given dimension and returns an AcquisitionData Object

+        '''

+        if not h5pyAvailable:

+            raise Exception("Error: h5py is not installed")

+        if self.filename is None:

+            return        

+        try:

+            

+                

+            with NexusFile(self.filename,'r') as file:    

+                try:

+                    dims = self.get_projection_dimensions()

+                except KeyError:

+                    pass

+                dims = file[self.data_path].shape

+                if ymin is None and ymax is None:

+                    data = np.array(file[self.data_path])

+                else:

+                    if ymin is None:

+                        ymin = 0

+                        if ymax > dims[1]:

+                            raise ValueError('ymax out of range')

+                        data = np.array(file[self.data_path][:,:ymax,:])

+                    elif ymax is None:        

+                        ymax = dims[1]

+                        if ymin < 0:

+                            raise ValueError('ymin out of range')

+                        data = np.array(file[self.data_path][:,ymin:,:])

+                    else:

+                        if ymax > dims[1]:

+                            raise ValueError('ymax out of range')

+                        if ymin < 0:

+                            raise ValueError('ymin out of range')

+                        

+                        data = np.array(file[self.data_path]

+                            [: , ymin:ymax , :] )

+                

+        except:

+            print("Error reading nexus file")

+            raise

+        

+        

+        try:

+            angles = self.get_projection_angles()

+        except KeyError as ke:

+            n = data.shape[0]

+            angles = np.linspace(0, n, n+1, dtype=np.float32)

+        

+        if ymax-ymin > 1:        

+            

+            geometry = AcquisitionGeometry('parallel', '3D', 

+                                       angles,

+                                       pixel_num_h          = dims[2],

+                                       pixel_size_h         = 1 ,

+                                       pixel_num_v          = ymax-ymin,

+                                       pixel_size_v         = 1,

+                                       dist_source_center   = None, 

+                                       dist_center_detector = None, 

+                                       channels             = 1)

+            return AcquisitionData(data, False, geometry=geometry, 

+                               dimension_labels=['angle','vertical','horizontal']) 

+        elif ymax-ymin == 1:

+            geometry = AcquisitionGeometry('parallel', '2D', 

+                                       angles,

+                                       pixel_num_h          = dims[2],

+                                       pixel_size_h         = 1 ,

+                                       dist_source_center   = None, 

+                                       dist_center_detector = None, 

+                                       channels             = 1)

+            return AcquisitionData(data.squeeze(), False, geometry=geometry, 

+                               dimension_labels=['angle','horizontal']) 

+    def get_acquisition_data_slice(self, y_slice=0):

+        return self.get_acquisition_data_subset(ymin=y_slice , ymax=y_slice+1)

+    def get_acquisition_data_whole(self):

+        with NexusFile(self.filename,'r') as file:    

+            try:

+                dims = self.get_projection_dimensions()

+            except KeyError:

+                print ("Warning: ")

+                dims = file[self.data_path].shape

+                

+            ymin = 0 

+            ymax = dims[1] - 1

+            

+            return self.get_acquisition_data_subset(ymin=ymin, ymax=ymax)

+            

+        

+    

+    def list_file_content(self):

+        try:

+            with NexusFile(self.filename,'r') as file:                

+                file.visit(print)

+        except:

+            print("Error reading nexus file")

+            raise  

+    def get_acquisition_data_batch(self, bmin=None, bmax=None):

+        if not h5pyAvailable:

+            raise Exception("Error: h5py is not installed")

+        if self.filename is None:

+            return        

+        try:

+            

+                

+            with NexusFile(self.filename,'r') as file:    

+                try:

+                    dims = self.get_projection_dimensions()

+                except KeyError:

+                    dims = file[self.data_path].shape

+                if bmin is None or bmax is None:

+                    raise ValueError('get_acquisition_data_batch: please specify fastest index batch limits')

+                    

+                if bmin >= 0 and bmin < bmax and bmax <= dims[0]:

+                    data = np.array(file[self.data_path][bmin:bmax])

+                else:

+                    raise ValueError('get_acquisition_data_batch: bmin {0}>0 bmax {1}<{2}'.format(bmin, bmax, dims[0]))

+                    

+        except:

+            print("Error reading nexus file")

+            raise

+        

+        

+        try:

+            angles = self.get_projection_angles()[bmin:bmax]

+        except KeyError as ke:

+            n = data.shape[0]

+            angles = np.linspace(0, n, n+1, dtype=np.float32)[bmin:bmax]

+        

+        if bmax-bmin > 1:        

+            

+            geometry = AcquisitionGeometry('parallel', '3D', 

+                                       angles,

+                                       pixel_num_h          = dims[2],

+                                       pixel_size_h         = 1 ,

+                                       pixel_num_v          = bmax-bmin,

+                                       pixel_size_v         = 1,

+                                       dist_source_center   = None, 

+                                       dist_center_detector = None, 

+                                       channels             = 1)

+            return AcquisitionData(data, False, geometry=geometry, 

+                               dimension_labels=['angle','vertical','horizontal']) 

+        elif bmax-bmin == 1:

+            geometry = AcquisitionGeometry('parallel', '2D', 

+                                       angles,

+                                       pixel_num_h          = dims[2],

+                                       pixel_size_h         = 1 ,

+                                       dist_source_center   = None, 

+                                       dist_center_detector = None, 

+                                       channels             = 1)

+            return AcquisitionData(data.squeeze(), False, geometry=geometry, 

+                               dimension_labels=['angle','horizontal']) 

+        

+    

           

 class XTEKReader(object):

     '''

diff --git a/Wrappers/Python/ccpi/processors.py b/Wrappers/Python/ccpi/processors.py
index a2d0446..6a9057a 100755
--- a/Wrappers/Python/ccpi/processors.py
+++ b/Wrappers/Python/ccpi/processors.py
@@ -19,6 +19,7 @@
 

 from ccpi.framework import DataProcessor, DataContainer, AcquisitionData,\

  AcquisitionGeometry, ImageGeometry, ImageData

+from ccpi.reconstruction.parallelbeam import alg as pbalg

 import numpy

 from scipy import ndimage

 

@@ -39,8 +40,8 @@ class Normalizer(DataProcessor):
     

     def __init__(self, flat_field = None, dark_field = None, tolerance = 1e-5):

         kwargs = {

-                  'flat_field'  : None, 

-                  'dark_field'  : None,

+                  'flat_field'  : flat_field, 

+                  'dark_field'  : dark_field,

                   # very small number. Used when there is a division by zero

                   'tolerance'   : tolerance

                   }

@@ -53,7 +54,8 @@ class Normalizer(DataProcessor):
             self.set_dark_field(dark_field)

     

     def check_input(self, dataset):

-        if dataset.number_of_dimensions == 3:

+        if dataset.number_of_dimensions == 3 or\

+           dataset.number_of_dimensions == 2:

                return True

         else:

             raise ValueError("Expected input dimensions is 2 or 3, got {0}"\

@@ -65,7 +67,7 @@ class Normalizer(DataProcessor):
                 raise ValueError('Dark Field should be 2D')

             elif len(numpy.shape(df)) == 2:

                 self.dark_field = df

-        elif issubclass(type(df), DataSet):

+        elif issubclass(type(df), DataContainer):

             self.dark_field = self.set_dark_field(df.as_array())

     

     def set_flat_field(self, df):

@@ -74,7 +76,7 @@ class Normalizer(DataProcessor):
                 raise ValueError('Flat Field should be 2D')

             elif len(numpy.shape(df)) == 2:

                 self.flat_field = df

-        elif issubclass(type(df), DataSet):

+        elif issubclass(type(df), DataContainer):

             self.flat_field = self.set_flat_field(df.as_array())

     

     @staticmethod

@@ -86,22 +88,40 @@ class Normalizer(DataProcessor):
             c[ ~ numpy.isfinite( c )] = tolerance # set to not zero if 0/0 

         return c

     

+    @staticmethod

+    def estimate_normalised_error(projection, flat, dark, delta_flat, delta_dark):

+        '''returns the estimated relative error of the normalised projection

+        

+        n = (projection - dark) / (flat - dark)

+        Dn/n = (flat-dark + projection-dark)/((flat-dark)*(projection-dark))*(Df/f + Dd/d)

+        ''' 

+        a = (projection - dark)

+        b = (flat-dark)

+        df = delta_flat / flat

+        dd = delta_dark / dark

+        rel_norm_error = (b + a) / (b * a) * (df + dd)

+        return rel_norm_error

+        

     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() ]

-                )

+        if projections.number_of_dimensions == 3:

+            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() ]

+                    )

+        elif projections.number_of_dimensions == 2:

+            a = Normalizer.normalize_projection(projections.as_array(), 

+                                                flat, dark, self.tolerance)

         y = type(projections)( a , True, 

                     dimension_labels=projections.dimension_labels,

                     geometry=projections.geometry)

@@ -388,3 +408,107 @@ class CenterOfRotationFinder(DataProcessor):
         

         return cor

 

+            

+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
\ No newline at end of file
diff --git a/Wrappers/Python/conda-recipe/meta.yaml b/Wrappers/Python/conda-recipe/meta.yaml
index 45d25f9..5de01c5 100644
--- a/Wrappers/Python/conda-recipe/meta.yaml
+++ b/Wrappers/Python/conda-recipe/meta.yaml
@@ -1,6 +1,6 @@
 package:
   name: ccpi-framework
-  version: 0.11.2
+  version: 0.11.3
 
 
 build:
@@ -9,6 +9,12 @@ build:
 #    - CIL_VERSION   
   number: 0
   
+  
+test:
+  requires:
+    - python-wget
+    - cvxpy # [not win]
+    
 requirements:
   build:
     - python
@@ -19,13 +25,8 @@ requirements:
     - numpy
     - scipy
     - matplotlib
-  run:
-    - python
-    - numpy
-    - cvxpy # [not win]
-    - scipy
-    - matplotlib
-	
+    - h5py
+  
 about:
   home: http://www.ccpi.ac.uk
   license:  Apache 2.0 License
diff --git a/Wrappers/Python/conda-recipe/run_test.py b/Wrappers/Python/conda-recipe/run_test.py
index d6bc340..5bf6538 100755
--- a/Wrappers/Python/conda-recipe/run_test.py
+++ b/Wrappers/Python/conda-recipe/run_test.py
@@ -20,6 +20,13 @@ from ccpi.optimisation.ops import LinearOperatorMatrix
 from ccpi.optimisation.ops import TomoIdentity

 from ccpi.optimisation.ops import Identity

 

+

+import numpy.testing

+import wget

+import os

+from ccpi.io.reader import NexusReader

+

+

 try:

     from cvxpy import *

     cvx_not_installable = False

@@ -879,6 +886,84 @@ class TestFunction(unittest.TestCase):
 #             s.split(2)

 # =============================================================================

 

+class TestNexusReader(unittest.TestCase):

+

+    def setUp(self):

+        wget.download('https://github.com/DiamondLightSource/Savu/raw/master/test_data/data/24737_fd.nxs')

+        self.filename = '24737_fd.nxs'

+

+    def tearDown(self):

+        os.remove(self.filename)

+

+

+    def testGetDimensions(self):

+        nr = NexusReader(self.filename)

+        self.assertEqual(nr.get_sinogram_dimensions(), (135, 91, 160), "Sinogram dimensions are not correct")

+        

+    def testGetProjectionDimensions(self):

+        nr = NexusReader(self.filename)

+        self.assertEqual(nr.get_projection_dimensions(), (91,135,160), "Projection dimensions are not correct")        

+        

+    def testLoadProjectionWithoutDimensions(self):

+        nr = NexusReader(self.filename)

+        projections = nr.load_projection()        

+        self.assertEqual(projections.shape, (91,135,160), "Loaded projection data dimensions are not correct")        

+

+    def testLoadProjectionWithDimensions(self):

+        nr = NexusReader(self.filename)

+        projections = nr.load_projection((slice(0,1), slice(0,135), slice(0,160)))        

+        self.assertEqual(projections.shape, (1,135,160), "Loaded projection data dimensions are not correct")        

+            

+    def testLoadProjectionCompareSingle(self):

+        nr = NexusReader(self.filename)

+        projections_full = nr.load_projection()

+        projections_part = nr.load_projection((slice(0,1), slice(0,135), slice(0,160))) 

+        numpy.testing.assert_array_equal(projections_part, projections_full[0:1,:,:])

+        

+    def testLoadProjectionCompareMulti(self):

+        nr = NexusReader(self.filename)

+        projections_full = nr.load_projection()

+        projections_part = nr.load_projection((slice(0,3), slice(0,135), slice(0,160))) 

+        numpy.testing.assert_array_equal(projections_part, projections_full[0:3,:,:])

+        

+    def testLoadProjectionCompareRandom(self):

+        nr = NexusReader(self.filename)

+        projections_full = nr.load_projection()

+        projections_part = nr.load_projection((slice(1,8), slice(5,10), slice(8,20))) 

+        numpy.testing.assert_array_equal(projections_part, projections_full[1:8,5:10,8:20])                

+        

+    def testLoadProjectionCompareFull(self):

+        nr = NexusReader(self.filename)

+        projections_full = nr.load_projection()

+        projections_part = nr.load_projection((slice(None,None), slice(None,None), slice(None,None))) 

+        numpy.testing.assert_array_equal(projections_part, projections_full[:,:,:])

+        

+    def testLoadFlatCompareFull(self):

+        nr = NexusReader(self.filename)

+        flats_full = nr.load_flat()

+        flats_part = nr.load_flat((slice(None,None), slice(None,None), slice(None,None))) 

+        numpy.testing.assert_array_equal(flats_part, flats_full[:,:,:])

+        

+    def testLoadDarkCompareFull(self):

+        nr = NexusReader(self.filename)

+        darks_full = nr.load_dark()

+        darks_part = nr.load_dark((slice(None,None), slice(None,None), slice(None,None))) 

+        numpy.testing.assert_array_equal(darks_part, darks_full[:,:,:])

+        

+    def testProjectionAngles(self):

+        nr = NexusReader(self.filename)

+        angles = nr.get_projection_angles()

+        self.assertEqual(angles.shape, (91,), "Loaded projection number of angles are not correct")        

+        

+    def test_get_acquisition_data_subset(self):

+        nr = NexusReader(self.filename)

+        key = nr.get_image_keys()

+        sl = nr.get_acquisition_data_subset(0,10)

+        data = nr.get_acquisition_data().subset(['vertical','horizontal'])

+        

+        self.assertTrue(sl.shape , (10,data.shape[1]))

+        

 

 if __name__ == '__main__':

     unittest.main()

+    

diff --git a/Wrappers/Python/setup.py b/Wrappers/Python/setup.py
index ef579e9..b584344 100644
--- a/Wrappers/Python/setup.py
+++ b/Wrappers/Python/setup.py
@@ -23,7 +23,7 @@ import os
 import sys
 
 
-cil_version='0.11.2'
+cil_version='0.11.3'
 
 setup(
     name="ccpi-framework",
diff --git a/Wrappers/Python/wip/multifile_nexus.py b/Wrappers/Python/wip/multifile_nexus.py
new file mode 100755
index 0000000..d1ad438
--- /dev/null
+++ b/Wrappers/Python/wip/multifile_nexus.py
@@ -0,0 +1,307 @@
+# -*- coding: utf-8 -*-

+"""

+Created on Wed Aug 15 16:00:53 2018

+

+@author: ofn77899

+"""

+

+import os

+from ccpi.io.reader import NexusReader

+

+from sys import getsizeof

+

+import matplotlib.pyplot as plt

+

+from ccpi.framework import DataProcessor, DataContainer

+from ccpi.processors import Normalizer

+from ccpi.processors import CenterOfRotationFinder

+import numpy

+

+

+class averager(object):

+    def __init__(self):

+        self.reset()

+        

+    def reset(self):

+        self.N = 0

+        self.avg = 0

+        self.min = 0

+        self.max = 0

+        self.var = 0

+        self.ske = 0

+        self.kur = 0

+    

+    def add_reading(self, val):

+        

+        if (self.N == 0):

+           self.avg = val;

+           self.min = val;

+           self.max = val;

+        elif (self.N == 1):

+        #//set min/max

+            self.max = val if val > self.max else self.max

+            self.min = val if val < self.min else self.min

+            

+      

+            thisavg = (self.avg + val)/2

+            #// initial value is in avg

+            self.var = (self.avg - thisavg)*(self.avg-thisavg) + (val - thisavg) * (val-thisavg)

+            self.ske = self.var * (self.avg - thisavg)

+            self.kur = self.var * self.var

+            self.avg = thisavg

+        else:

+            self.max = val if val > self.max else self.max

+            self.min = val if val < self.min else self.min

+            

+            M = self.N

+

+            #// b-factor =(<v>_N + v_(N+1)) / (N+1)

+            #float b = (val -avg) / (M+1);

+            b = (val -self.avg) / (M+1)

+            

+            self.kur = self.kur + (M *(b*b*b*b) * (1+M*M*M))- (4* b * self.ske) + (6 * b *b * self.var * (M-1))

+    

+            self.ske = self.ske + (M * b*b*b *(M-1)*(M+1)) - (3 * b * self.var * (M-1))

+    

+            #//var = var * ((M-1)/M) + ((val - avg)*(val - avg)/(M+1)) ;

+            self.var = self.var * ((M-1)/M) + (b * b * (M+1))

+            

+            self.avg = self.avg * (M/(M+1)) + val / (M+1)

+    

+        self.N += 1

+    

+    def stats(self, vector):

+        i = 0

+        while i < vector.size:

+            self.add_reading(vector[i])

+            i+=1

+

+avg = averager()

+a = numpy.linspace(0,39,40)

+avg.stats(a)

+print ("average" , avg.avg, a.mean())

+print ("variance" , avg.var, a.var())

+b = a - a.mean()

+b *= b

+b = numpy.sqrt(sum(b)/(a.size-1))

+print ("std" , numpy.sqrt(avg.var), b)

+#%%         

+

+class DataStatMoments(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, axis, skewness=False, kurtosis=False, offset=0):

+        kwargs = {

+                  'axis'     : axis,

+                  'skewness' : skewness,

+                  'kurtosis' : kurtosis,

+                  'offset'   : offset,

+                  }

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

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

+        

+    

+    def check_input(self, dataset):

+        #self.N = dataset.get_dimension_size(self.axis)

+        return True

+    

+    @staticmethod

+    def add_sample(dataset, N, axis, stats=None, offset=0):

+        #dataset = self.get_input()

+        if (N == 0):

+            # get the axis number along to calculate the stats

+            

+            

+            #axs = dataset.get_dimension_size(self.axis)

+            # create a placeholder for the output

+            if stats is None:

+                ax = dataset.get_dimension_axis(axis)

+                shape = [dataset.shape[i] for i in range(len(dataset.shape)) if i != ax]

+                # output has 4 components avg, std, skewness and kurtosis + last avg+ (val-thisavg)

+                shape.insert(0, 4+2)

+                stats = numpy.zeros(shape)

+                

+            stats[0] = dataset.subset(**{axis:N-offset}).array[:]

+            

+            #avg = val

+        elif (N == 1):

+            # val

+            stats[5] = dataset.subset(**{axis:N-offset}).array

+            stats[4] = stats[0] + stats[5]

+            stats[4] /= 2         # thisavg

+            stats[5] -= stats[4]  # (val - thisavg) 

+            

+            #float thisavg = (avg + val)/2;

+            

+            #// initial value is in avg

+            #var = (avg - thisavg)*(avg-thisavg) + (val - thisavg) * (val-thisavg);

+            stats[1] = stats[5] * stats[5] + stats[5] * stats[5]

+            #skewness = var * (avg - thisavg);

+            stats[2] = stats[1] * stats[5]

+            #kurtosis = var * var;

+            stats[3] = stats[1] * stats[1]

+            #avg = thisavg;

+            stats[0] = stats[4]

+        else:

+        

+            #float M = (float)N;

+            M = N

+            #val

+            stats[4] = dataset.subset(**{axis:N-offset}).array

+            #// b-factor =(<v>_N + v_(N+1)) / (N+1)

+            stats[5] = stats[4] - stats[0]

+            stats[5] /= (M+1)    # b factor

+            #float b = (val -avg) / (M+1);

+    

+            #kurtosis = kurtosis + (M *(b*b*b*b) * (1+M*M*M))- (4* b * skewness) + (6 * b *b * var * (M-1));

+            #if self.kurtosis:

+            #    stats[3] += (M * stats[5] * stats[5] * stats[5] * stats[5]) - \

+            #                (4 * stats[5] * stats[2]) + \

+            #                (6 * stats[5] * stats[5] * stats[1] * (M-1))

+                        

+            #skewness = skewness + (M * b*b*b *(M-1)*(M+1)) - (3 * b * var * (M-1));

+            #if self.skewness:

+            #    stats[2] = stats[2] +  (M * stats[5]* stats[5] * stats[5] * (M-1)*(M-1) ) -\

+            #                3 * stats[5] * stats[1] * (M-1) 

+            #//var = var * ((M-1)/M) + ((val - avg)*(val - avg)/(M+1)) ;

+            #var = var * ((M-1)/M) + (b * b * (M+1));

+            stats[1] = ((M-1)/M) * stats[1] + (stats[5] * stats[5] * (M+1))

+           

+            #avg = avg * (M/(M+1)) + val / (M+1)

+            stats[0] = stats[0] * (M/(M+1)) + stats[4] / (M+1)

+    

+        N += 1

+        return stats , N

+    

+           

+    def process(self):

+        

+        data = self.get_input()

+        

+        #stats, i = add_sample(0)

+        N = data.get_dimension_size(self.axis)

+        ax = data.get_dimension_axis(self.axis)

+        stats = None

+        i = 0

+        while i < N:

+            stats , i = DataStatMoments.add_sample(data, i, self.axis, stats, offset=self.offset)

+        self.offset += N

+        labels = ['StatisticalMoments']

+        

+        labels += [data.dimension_labels[i] \

+                   for i in range(len(data.dimension_labels)) if i != ax]

+        y = DataContainer( stats[:4] , False, 

+                    dimension_labels=labels)

+        return y

+

+directory = r'E:\Documents\Dataset\CCPi\Nexus_test'

+data_path="entry1/instrument/pco1_hw_hdf_nochunking/data"

+

+reader = NexusReader(os.path.join( os.path.abspath(directory) , '74331.nxs'))

+

+print ("read flat")

+read_flat = NexusReader(os.path.join( os.path.abspath(directory) , '74240.nxs'))

+read_flat.data_path = data_path

+flatsslice = read_flat.get_acquisition_data_whole()

+avg = DataStatMoments('angle')

+

+avg.set_input(flatsslice)

+flats = avg.get_output()

+

+ave = averager()

+ave.stats(flatsslice.array[:,0,0])

+

+print ("avg" , ave.avg, flats.array[0][0][0])

+print ("var" , ave.var, flats.array[1][0][0])

+

+print ("read dark")

+read_dark = NexusReader(os.path.join( os.path.abspath(directory) , '74243.nxs'))

+read_dark.data_path = data_path

+

+## darks are very many so we proceed in batches

+total_size = read_dark.get_projection_dimensions()[0]

+print ("total_size", total_size)

+

+batchsize = 40

+if batchsize > total_size:

+    batchlimits = [batchsize * (i+1) for i in range(int(total_size/batchsize))] + [total_size-1]

+else:

+    batchlimits = [total_size]

+#avg.N = 0

+avg.offset = 0

+N = 0

+for batch in range(len(batchlimits)):  

+    print ("running batch " , batch)

+    bmax = batchlimits[batch]

+    bmin = bmax-batchsize

+    

+    darksslice = read_dark.get_acquisition_data_batch(bmin,bmax)

+    if batch == 0:

+        #create stats

+        ax = darksslice.get_dimension_axis('angle')

+        shape = [darksslice.shape[i] for i in range(len(darksslice.shape)) if i != ax]

+        # output has 4 components avg, std, skewness and kurtosis + last avg+ (val-thisavg)

+        shape.insert(0, 4+2)

+        print ("create stats shape ", shape)

+        stats = numpy.zeros(shape)

+    print ("N" , N)        

+    #avg.set_input(darksslice)

+    i = bmin

+    while i < bmax:

+        stats , i = DataStatMoments.add_sample(darksslice, i, 'angle', stats, bmin)

+        print ("{0}-{1}-{2}".format(bmin, i , bmax ) )

+    

+darks = stats

+#%%

+

+fig = plt.subplot(2,2,1)

+fig.imshow(flats.subset(StatisticalMoments=0).array)

+fig = plt.subplot(2,2,2)

+fig.imshow(numpy.sqrt(flats.subset(StatisticalMoments=1).array))

+fig = plt.subplot(2,2,3)

+fig.imshow(darks[0])

+fig = plt.subplot(2,2,4)

+fig.imshow(numpy.sqrt(darks[1]))

+plt.show()

+

+

+#%%

+norm = Normalizer(flat_field=flats.array[0,200,:], dark_field=darks[0,200,:])

+#norm.set_flat_field(flats.array[0,200,:])

+#norm.set_dark_field(darks.array[0,200,:])

+norm.set_input(reader.get_acquisition_data_slice(200))

+

+n = Normalizer.normalize_projection(norm.get_input().as_array(), flats.array[0,200,:], darks[0,200,:], 1e-5)

+#dn_n= Normalizer.estimate_normalised_error(norm.get_input().as_array(), flats.array[0,200,:], darks[0,200,:],

+#                                           numpy.sqrt(flats.array[1,200,:]), numpy.sqrt(darks[1,200,:]))

+#%%

+

+

+#%%

+fig = plt.subplot(2,1,1)

+

+

+fig.imshow(norm.get_input().as_array())

+fig = plt.subplot(2,1,2)

+fig.imshow(n)

+

+#fig = plt.subplot(3,1,3)

+#fig.imshow(dn_n)

+

+

+plt.show()

+

+

+

+

+

+