From 03021b2c22e86821c1df1f89d5774062120bbc3c Mon Sep 17 00:00:00 2001 From: Srikanth Nagella Date: Fri, 16 Mar 2018 05:57:13 +0000 Subject: Diamond and Xtek Readers --- Wrappers/Python/ccpi/io/__init__.py | 18 +++ Wrappers/Python/ccpi/io/reader.py | 301 ++++++++++++++++++++++++++++++++++++ Wrappers/Python/test/TestReader.py | 134 ++++++++++++++++ 3 files changed, 453 insertions(+) create mode 100644 Wrappers/Python/ccpi/io/__init__.py create mode 100644 Wrappers/Python/ccpi/io/reader.py create mode 100644 Wrappers/Python/test/TestReader.py 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 -- cgit v1.2.3