diff options
| -rw-r--r-- | Wrappers/Python/ccpi/io/__init__.py | 18 | ||||
| -rw-r--r-- | Wrappers/Python/ccpi/io/reader.py | 301 | ||||
| -rw-r--r-- | Wrappers/Python/test/TestReader.py | 134 |
3 files changed, 453 insertions, 0 deletions
diff --git a/Wrappers/Python/ccpi/io/__init__.py b/Wrappers/Python/ccpi/io/__init__.py new file mode 100644 index 0000000..9233d7a --- /dev/null +++ b/Wrappers/Python/ccpi/io/__init__.py @@ -0,0 +1,18 @@ +# -*- 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.
\ No newline at end of file diff --git a/Wrappers/Python/ccpi/io/reader.py b/Wrappers/Python/ccpi/io/reader.py new file mode 100644 index 0000000..4f7e288 --- /dev/null +++ b/Wrappers/Python/ccpi/io/reader.py @@ -0,0 +1,301 @@ +# -*- coding: utf-8 -*-
+# This work is part of the Core Imaging Library developed by
+# Visual Analytics and Imaging System Group of the Science Technology
+# Facilities Council, STFC
+
+# Copyright 2018 Jakob Jorgensen, Daniil Kazantsev, Edoardo Pasca and Srikanth Nagella
+
+# 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.
+
+'''
+This is a reader module with classes for loading 3D datasets.
+
+@author: Mr. Srikanth Nagella
+'''
+from ccpi.framework import AcquisitionGeometry
+from ccpi.framework import AcquisitionData
+import numpy as np
+import os
+
+h5pyAvailable = True
+try:
+ from h5py import File as NexusFile
+except:
+ h5pyAvailable = False
+
+pilAvailable = True
+try:
+ from PIL import Image
+except:
+ pilAvailable = False
+
+class NexusReader(object):
+ '''
+ Reader class for loading Nexus files.
+ '''
+
+ def __init__(self, nexusFilename=None):
+ '''
+ This takes in input as filename and loads the data dataset.
+ '''
+ self.flat = None
+ self.dark = None
+ self.angles = None
+ self.geometry = None
+ self.filename = nexusFilename
+
+ def load(self, dimensions=None, image_key_id=0):
+ '''
+ This is generic loading function of flat field, dark field and projection data.
+ '''
+ if not h5pyAvailable:
+ raise Exception("Error: h5py is not installed")
+ if self.filename is None:
+ return
+ try:
+ with NexusFile(self.filename,'r') as file:
+ image_keys = np.array(file['entry1/tomo_entry/instrument/detector/image_key'])
+ projections = None
+ if dimensions == None:
+ projections = np.array(file['entry1/tomo_entry/data/data'])
+ result = projections[image_keys==image_key_id]
+ return result
+ else:
+ #When dimensions are specified they need to be mapped to image_keys
+ index_array = np.where(image_keys==image_key_id)
+ 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])
+ return result
+ except:
+ print("Error reading nexus file")
+ raise
+
+ def loadProjection(self, dimensions=None):
+ '''
+ Loads the projection data from the nexus file.
+ returns: numpy array with projection data
+ '''
+ return self.load(dimensions, 0)
+
+ def loadFlat(self, dimensions=None):
+ '''
+ Loads the flat field data from the nexus file.
+ returns: numpy array with flat field data
+ '''
+ return self.load(dimensions, 1)
+
+ def loadDark(self, dimensions=None):
+ '''
+ Loads the Dark field data from the nexus file.
+ returns: numpy array with dark field data
+ '''
+ return self.load(dimensions, 2)
+
+ def getProjectionAngles(self):
+ '''
+ This function returns the projection angles
+ '''
+ if not h5pyAvailable:
+ raise Exception("Error: h5py is not installed")
+ if self.filename is None:
+ 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'])
+ return angles[image_keys==0]
+ except:
+ print("Error reading nexus file")
+ raise
+
+
+ def getSinogramDimensions(self):
+ '''
+ Return the dimensions of the dataset
+ '''
+ if not h5pyAvailable:
+ raise Exception("Error: h5py is not installed")
+ 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'])
+ dims = list(projections.shape)
+ dims[0] = dims[1]
+ dims[1] = np.sum(image_keys==0)
+ return tuple(dims)
+ except:
+ print("Error reading nexus file")
+ raise
+
+ def getProjectionDimensions(self):
+ '''
+ Return the dimensions of the dataset
+ '''
+ if not h5pyAvailable:
+ raise Exception("Error: h5py is not installed")
+ 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'])
+ dims = list(projections.shape)
+ dims[0] = np.sum(image_keys==0)
+ return tuple(dims)
+ except:
+ print("Error reading nexus file")
+ raise
+
+ def getAcquisitionData(self, dimensions=None):
+ '''
+ This method load the acquisition data and given dimension and returns an AcquisitionData Object
+ '''
+ data = self.loadProjection(dimensions)
+ return AcquisitionData(data)
+
+ def getAcquisitionGeometry(self, dimensions=None):
+ '''
+ This method creates the acquisition geometry from a given dimension and return the created object
+ '''
+ return AcquisitionGeometry('parallel', '3D', self.getProjectionAngles())
+
+class XTEKReader(object):
+ '''
+ Reader class for loading XTEK files
+ '''
+
+ def __init__(self, xtekConfigFilename=None):
+ '''
+ This takes in the xtek config filename and loads the dataset and the
+ required geometry parameters
+ '''
+ self.projections = None
+ self.geometry = {}
+ self.filename = xtekConfigFilename
+ self.load()
+
+ def load(self):
+ pixel_num_h = 0
+ pixel_num_v = 0
+ xpixel_size = 0
+ ypixel_size = 0
+ source_x = 0
+ detector_x = 0
+ with open(self.filename) as f:
+ content = f.readlines()
+ content = [x.strip() for x in content]
+ for line in content:
+ if line.startswith("SrcToObject"):
+ source_x = float(line.split('=')[1])
+ elif line.startswith("SrcToDetector"):
+ detector_x = float(line.split('=')[1])
+ elif line.startswith("DetectorPixelsY"):
+ pixel_num_v = int(line.split('=')[1])
+ #self.num_of_vertical_pixels = self.calc_v_alighment(self.num_of_vertical_pixels, self.pixels_per_voxel)
+ elif line.startswith("DetectorPixelsX"):
+ pixel_num_h = int(line.split('=')[1])
+ elif line.startswith("DetectorPixelSizeX"):
+ xpixel_size = float(line.split('=')[1])
+ elif line.startswith("DetectorPixelSizeY"):
+ ypixel_size = float(line.split('=')[1])
+ elif line.startswith("Projections"):
+ self.num_projections = int(line.split('=')[1])
+ elif line.startswith("InitialAngle"):
+ self.initial_angle = float(line.split('=')[1])
+ elif line.startswith("Name"):
+ self.experiment_name = line.split('=')[1]
+ elif line.startswith("Scattering"):
+ self.scattering = float(line.split('=')[1])
+ elif line.startswith("WhiteLevel"):
+ self.white_level = float(line.split('=')[1])
+ elif line.startswith("MaskRadius"):
+ self.mask_radius = float(line.split('=')[1])
+
+ #Read Angles
+ angles = self.readAngles()
+ self.geometry = AcquisitionGeometry('cone', '3D', angles, pixel_num_h, xpixel_size, pixel_num_v, ypixel_size, -1 * source_x,
+ detector_x - source_x,
+ )
+
+ def readAngles(self):
+ """
+ Read the angles file .ang or _ctdata.txt file and returns the angles
+ as an numpy array.
+ """
+ input_path = os.path.dirname(self.filename)
+ angles_ctdata_file = os.path.join(input_path, '_ctdata.txt')
+ angles_named_file = os.path.join(input_path, self.experiment_name+'.ang')
+ angles = np.zeros(self.num_projections,dtype='f')
+ #look for _ctdata.txt
+ if os.path.exists(angles_ctdata_file):
+ #read txt file with angles
+ with open(angles_ctdata_file) as f:
+ content = f.readlines()
+ #skip firt three lines
+ #read the middle value of 3 values in each line as angles in degrees
+ index = 0
+ for line in content[3:]:
+ self.angles[index]=float(line.split(' ')[1])
+ index+=1
+ angles = np.deg2rad(self.angles+self.initial_angle);
+ elif os.path.exists(angles_named_file):
+ #read the angles file which is text with first line as header
+ with open(angles_named_file) as f:
+ content = f.readlines()
+ #skip first line
+ index = 0
+ for line in content[1:]:
+ angles[index] = float(line.split(':')[1])
+ index+=1
+ angles = np.flipud(angles+self.initial_angle) #angles are in the reverse order
+ else:
+ raise RuntimeError("Can't find angles file")
+ return angles
+
+ def loadProjection(self, dimensions=None):
+ '''
+ This method reads the projection images from the directory and returns a numpy array
+ '''
+ if not pilAvailable:
+ raise('Image library pillow is not installed')
+ if dimensions != None:
+ raise('Extracting subset of data is not implemented')
+ input_path = os.path.dirname(self.filename)
+ pixels = np.zeros((self.num_projections, self.geometry.pixel_num_h, self.geometry.pixel_num_v), dtype='float32')
+ for i in range(1, self.num_projections+1):
+ im = Image.open(os.path.join(input_path,self.experiment_name+"_%04d"%i+".tif"))
+ pixels[i-1,:,:] = np.fliplr(np.transpose(np.array(im))) ##Not sure this is the correct way to populate the image
+
+ #normalising the data
+ #TODO: Move this to a processor
+ pixels = pixels - (self.white_level*self.scattering)/100.0
+ pixels[pixels < 0.0] = 0.000001 # all negative values to approximately 0 as the std log of zero and non negative number is not defined
+ return pixels
+
+ def getAcquisitionData(self, dimensions=None):
+ '''
+ This method load the acquisition data and given dimension and returns an AcquisitionData Object
+ '''
+ data = self.loadProjection(dimensions)
+ return AcquisitionData(data)
+
+ def getAcquisitionGeometry(self, dimensions=None):
+ '''
+ This method creates the acquisition geometry from a given dimension and return the created object
+ '''
+ return self.geometry
\ No newline at end of file diff --git a/Wrappers/Python/test/TestReader.py b/Wrappers/Python/test/TestReader.py new file mode 100644 index 0000000..51db052 --- /dev/null +++ b/Wrappers/Python/test/TestReader.py @@ -0,0 +1,134 @@ +# -*- coding: utf-8 -*-
+# This work is part of the Core Imaging Library developed by
+# Visual Analytics and Imaging System Group of the Science Technology
+# Facilities Council, STFC
+
+# Copyright 2018 Jakob Jorgensen, Daniil Kazantsev, Edoardo Pasca and Srikanth Nagella
+
+# 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.
+
+'''
+Unit tests for Readers
+
+@author: Mr. Srikanth Nagella
+'''
+import unittest
+
+import numpy.testing
+import wget
+import os
+from ccpi.io.reader import NexusReader
+from ccpi.io.reader import XTEKReader
+#@unittest.skip
+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.getSinogramDimensions(), (135, 91, 160), "Sinogram dimensions are not correct")
+
+ def testGetProjectionDimensions(self):
+ nr = NexusReader(self.filename)
+ self.assertEqual(nr.getProjectionDimensions(), (91,135,160), "Projection dimensions are not correct")
+
+ def testLoadProjectionWithoutDimensions(self):
+ nr = NexusReader(self.filename)
+ projections = nr.loadProjection()
+ self.assertEqual(projections.shape, (91,135,160), "Loaded projection data dimensions are not correct")
+
+ def testLoadProjectionWithDimensions(self):
+ nr = NexusReader(self.filename)
+ projections = nr.loadProjection((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.loadProjection()
+ projections_part = nr.loadProjection((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.loadProjection()
+ projections_part = nr.loadProjection((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.loadProjection()
+ projections_part = nr.loadProjection((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.loadProjection()
+ projections_part = nr.loadProjection((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.loadFlat()
+ flats_part = nr.loadFlat((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.loadDark()
+ darks_part = nr.loadDark((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.getProjectionAngles()
+ self.assertEqual(angles.shape, (91,), "Loaded projection number of angles are not correct")
+
+class TestXTEKReader(unittest.TestCase):
+
+ def setUp(self):
+ testpath, filename = os.path.split(os.path.realpath(__file__))
+ testpath = os.path.normpath(os.path.join(testpath, '..','..','..'))
+ self.filename = os.path.join(testpath,'data','SophiaBeads','SophiaBeads_64_averaged.xtekct')
+
+ def tearDown(self):
+ pass
+
+ def testLoad(self):
+ xtekReader = XTEKReader(self.filename)
+ self.assertEqual(xtekReader.geometry.pixel_num_h, 500, "Detector pixel X size is not correct")
+ self.assertEqual(xtekReader.geometry.pixel_num_v, 500, "Detector pixel Y size is not correct")
+ self.assertEqual(xtekReader.geometry.dist_source_center, -80.6392412185669, "Distance from source to center is not correct")
+ self.assertEqual(xtekReader.geometry.dist_center_detector, (1007.006 - 80.6392412185669), "Distance from center to detector is not correct")
+
+ def testReadAngles(self):
+ xtekReader = XTEKReader(self.filename)
+ angles = xtekReader.readAngles()
+ self.assertEqual(angles.shape, (63,), "Angles doesn't match")
+ self.assertAlmostEqual(angles[46], -085.717, 3, "46th Angles doesn't match")
+
+ def testLoadProjection(self):
+ xtekReader = XTEKReader(self.filename)
+ pixels = xtekReader.loadProjection()
+ self.assertEqual(pixels.shape, (63, 500, 500), "projections shape doesn't match")
+
+
+
+if __name__ == "__main__":
+ #import sys;sys.argv = ['', 'TestNexusReader.testLoad']
+ unittest.main()
\ No newline at end of file |