diff --git a/readers/NSIDataReader.py b/readers/NSIDataReader.py new file mode 100644 index 0000000..1db7e23 --- /dev/null +++ b/readers/NSIDataReader.py @@ -0,0 +1,337 @@ +# Copyright 2025 UKRI-STFC +# Copyright 2025 University of Manchester + +# 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. +# +# Authors: +# Evan Kiely (WMG) + +from cil.framework import AcquisitionData, AcquisitionGeometry +from cil.io.TIFF import TIFFStackReader +import warnings +import numpy as np +import os + + +class NSIDataReader(object): + + def __init__(self, + **kwargs): + '''Basic reader for xtekct files + + Parameters + ---------- + + + file_name: str with full path to file. TIFFs will be retrieved from the same path. + + roi: dictionary with roi to load + {'angle': (start, end, step), + 'horizontal': (start, end, step), + 'vertical': (start, end, step)} + Files are stacked along axis_0. axis_1 and axis_2 correspond + to row and column dimensions, respectively. + Files are stacked in alphabetic order. + To skip projections or to change number of projections to load, + adjust 'angle'. For instance, 'angle': (100, 300) + will skip first 100 projections and will load 200 projections. + 'angle': -1 is a shortcut to load all elements along axis. + Start and end can be specified as None which is equivalent + to start = 0 and end = load everything to the end, respectively. + Start and end also can be negative. + + normalise: bool, normalises loaded projections by detector + white level (I_0). Default value is True + + fliplr: bool, default = False, flip projections in the left-right direction + (about vertical axis) + + mode: str, 'bin' (default) or 'slice'. In bin mode, 'step' number + of pixels is binned together, values of resulting binned + pixels are calculated as average. + In 'slice' mode 'step' defines standard numpy slicing. + Note: in general + output array size in bin mode != output array size in slice mode + + Output + ------ + + Acquisition data with corresponding geometry, arranged as ['angle', horizontal'] + if a single slice is loaded and ['vertical, 'angle', horizontal'] + if more than 1 slices are loaded. + + ''' + + self.file_name = kwargs.get('file_name', None) + self.roi = kwargs.get('roi', {'angle': -1, 'horizontal': -1, 'vertical': -1}) + self.normalise = kwargs.get('normalise', True) + self.mode = kwargs.get('mode', 'bin') + self.fliplr = kwargs.get('fliplr', False) + + if self.file_name is not None: + self.set_up(file_name = self.file_name, + roi = self.roi, + normalise = self.normalise, + mode = self.mode, + fliplr = self.fliplr) + + def set_up(self, + file_name = None, + roi = {'angle': -1, 'horizontal': -1, 'vertical': -1}, + normalise = True, + mode = 'bin', + fliplr = False): + + self.file_name = file_name + self.roi = roi + self.normalise = normalise + self.mode = mode + self.fliplr = fliplr + self.tiff_directory_path = os.path.dirname(self.file_name) + + if self.file_name == None: + raise Exception('Path to file is required.') + + # check if data file exists + if not(os.path.isfile(self.file_name)): + raise Exception('File\n {}\n does not exist.'.format(self.file_name)) + + # check labels + for key in self.roi.keys(): + if key not in ['angle', 'horizontal', 'vertical']: + raise Exception("Wrong label. One of the following is expected: angle, horizontal, vertical") + + roi = self.roi.copy() + + if 'angle' not in roi.keys(): + roi['angle'] = -1 + + if 'horizontal' not in roi.keys(): + roi['horizontal'] = -1 + + if 'vertical' not in roi.keys(): + roi['vertical'] = -1 + + # parse data file + with open(self.file_name, 'r', errors='replace') as f: + content = f.readlines() + + content = [x.strip() for x in content] + + #initialise parameters + detector_offset_h = 0 + detector_offset_v = 0 + object_tilt_deg = 0 + object_offset_x = None + object_roll_deg = None + centre_of_rotation_top = 0 + centre_of_rotation_bottom = 0 + angular_step = 0 + + for line in content: + # number of projections + if line.startswith(""): + num_projections = int(line.split('>')[1]) + + # number of pixels along Y axis + elif line.startswith(""): + pixel_num_h_0 = int(line.split('>')[1]) + + # number of pixels along X axis + elif line.startswith(""): + pixel_num_v_0 = int(line.split('>')[1]) + + # pixel size along X and y axis + elif line.startswith(""): + pixel_size_h_0 = float(line.split('>')[1]) + pixel_size_v_0 = pixel_size_h_0 + + elif line.startswith(""): + source_to_det = float(line.split('>')[1]) + + # source to detector distance + elif line.startswith(""): + source_to_origin = float(line.split('>')[1]) + + # angular step + elif line.startswith(""): + if angular_step == 0: + angular_step = float(line.split('>')[1]) + + + self._experiment_name = 'scan' + + self._roi_par = [[0, num_projections, 1] ,[0, pixel_num_v_0, 1], [0, pixel_num_h_0, 1]] + + for key in roi.keys(): + if key == 'angle': + idx = 0 + elif key == 'vertical': + idx = 1 + elif key == 'horizontal': + idx = 2 + if roi[key] != -1: + for i in range(2): + if roi[key][i] != None: + if roi[key][i] >= 0: + self._roi_par[idx][i] = roi[key][i] + else: + self._roi_par[idx][i] = self._roi_par[idx][1]+roi[key][i] + if len(roi[key]) > 2: + if roi[key][2] != None: + if roi[key][2] > 0: + self._roi_par[idx][2] = roi[key][2] + else: + raise Exception("Negative step is not allowed") + + if self.mode == 'bin': + # calculate number of pixels and pixel size + pixel_num_v = (self._roi_par[1][1] - self._roi_par[1][0]) // self._roi_par[1][2] + pixel_num_h = (self._roi_par[2][1] - self._roi_par[2][0]) // self._roi_par[2][2] + pixel_size_v = pixel_size_v_0 * self._roi_par[1][2] + pixel_size_h = pixel_size_h_0 * self._roi_par[2][2] + else: # slice + pixel_num_v = int(np.ceil((self._roi_par[1][1] - self._roi_par[1][0]) / self._roi_par[1][2])) + pixel_num_h = int(np.ceil((self._roi_par[2][1] - self._roi_par[2][0]) / self._roi_par[2][2])) + + pixel_size_v = pixel_size_v_0 + pixel_size_h = pixel_size_h_0 + + det_start_0 = -(pixel_num_h_0 / 2) + det_start = det_start_0 + self._roi_par[2][0] + det_end = det_start + pixel_num_h * self._roi_par[2][2] + det_pos_h = (det_start + det_end) * 0.5 * pixel_size_h_0 + detector_offset_h + + det_start_0 = -(pixel_num_v_0 / 2) + det_start = det_start_0 + self._roi_par[1][0] + det_end = det_start + pixel_num_v * self._roi_par[1][2] + det_pos_v = (det_start + det_end) * 0.5 * pixel_size_v_0 + detector_offset_v + + #angles from data file ignore *.ang and _ctdata.txt as not correct + angles = np.asarray( [ angular_step * proj for proj in range(num_projections) ] , dtype=np.float32) + + if self.mode == 'bin': + n_elem = (self._roi_par[0][1] - self._roi_par[0][0]) // self._roi_par[0][2] + shape = (n_elem, self._roi_par[0][2]) + angles = angles[self._roi_par[0][0]:(self._roi_par[0][0] + n_elem * self._roi_par[0][2])].reshape(shape).mean(1) + else: + angles = angles[slice(self._roi_par[0][0], self._roi_par[0][1], self._roi_par[0][2])] + + + if object_offset_x == None: + object_offset_x = (centre_of_rotation_bottom+centre_of_rotation_top)* 0.5 * source_to_origin /source_to_det + + if object_roll_deg == None: + x = (centre_of_rotation_top-centre_of_rotation_bottom) + y = pixel_num_v_0 * pixel_size_v_0 + object_roll = -np.arctan2(x, y) + else: + object_roll = object_roll_deg * np.pi /180. + + object_tilt = -object_tilt_deg * np.pi /180. + + tilt_matrix = np.eye(3) + tilt_matrix[1][1] = tilt_matrix[2][2] = np.cos(object_tilt) + tilt_matrix[1][2] = -np.sin(object_tilt) + tilt_matrix[2][1] = np.sin(object_tilt) + + roll_matrix = np.eye(3) + roll_matrix[0][0] = roll_matrix[2][2] = np.cos(object_roll) + roll_matrix[0][2] = np.sin(object_roll) + roll_matrix[2][0] = -np.sin(object_roll) + + #order of construction may be reversed, but unlikely to have both in a dataset + rot_matrix = np.matmul(tilt_matrix,roll_matrix) + rotation_axis_direction = rot_matrix.dot([0,0,1]) + + + if self.fliplr: + origin = 'top-left' + else: + origin = 'top-right' + + if pixel_num_v == 1 and (self._roi_par[1][0]+self._roi_par[1][1]) // 2 == pixel_num_v_0 // 2: + self._ag = AcquisitionGeometry.create_Cone2D(source_position=[0, -source_to_origin], + rotation_axis_position=[-object_offset_x, 0], + detector_position=[-det_pos_h, source_to_det-source_to_origin]) + self._ag.set_angles(angles, + angle_unit='degree') + + self._ag.set_panel(pixel_num_h, pixel_size=pixel_size_h, origin=origin) + + self._ag.set_labels(labels=['angle', 'horizontal']) + else: + self._ag = AcquisitionGeometry.create_Cone3D(source_position=[0, -source_to_origin, 0], + rotation_axis_position=[-object_offset_x, 0, 0], + rotation_axis_direction=rotation_axis_direction, + detector_position=[-det_pos_h, source_to_det-source_to_origin, det_pos_v]) + self._ag.set_angles(angles, + angle_unit='degree') + + self._ag.set_panel(( pixel_num_h, pixel_num_v), + pixel_size=( pixel_size_v, pixel_size_h), + origin=origin) + + self._ag.set_labels(labels=['angle', 'vertical', 'horizontal']) + + def get_geometry(self): + + ''' + Return AcquisitionGeometry object + ''' + + return self._ag + def get_roi(self): + '''returns the roi''' + roi = self._roi_par[:] + if self._ag.dimension == '2D': + roi.pop(1) + + roidict = {} + for i,el in enumerate(roi): + # print (i, el) + roidict['axis_{}'.format(i)] = tuple(el) + return roidict + + def read(self): + + ''' + Reads projections and return AcquisitionData container + ''' + + reader = TIFFStackReader() + + roi = self.get_roi() + + reader.set_up(file_name = self.tiff_directory_path, + roi=roi, mode=self.mode) + + ad = reader.read_as_AcquisitionData(self._ag) + + if (self.normalise): + #ad.array[ad.array < 1] = 1 + # cast the data read to float32 + ad = ad / np.float32(float(1)) + + + if self.fliplr: + dim = ad.get_dimension_axis('horizontal') + ad.array = np.flip(ad.array, dim) + + return ad + + def load_projections(self): + '''alias of read for backward compatibility''' + return self.read() + +