From 68200aa55722b583fbd2934ae81252ee4b6e4cd2 Mon Sep 17 00:00:00 2001 From: Edoardo Pasca Date: Tue, 6 Mar 2018 12:16:53 +0000 Subject: New datasetprocessors (#33) * added simple test for SinogramData * initial version of processors This contains a version of a simple normalization algorithm. --- Wrappers/Python/ccpi/framework.py | 6 +- Wrappers/Python/ccpi/processors.py | 164 +++++++++++++++++++++++++++++++++++++ 2 files changed, 169 insertions(+), 1 deletion(-) create mode 100755 Wrappers/Python/ccpi/processors.py diff --git a/Wrappers/Python/ccpi/framework.py b/Wrappers/Python/ccpi/framework.py index 22324be..2b0ba76 100644 --- a/Wrappers/Python/ccpi/framework.py +++ b/Wrappers/Python/ccpi/framework.py @@ -105,7 +105,7 @@ class DataSet(object): self.shape = numpy.shape(array) self.number_of_dimensions = len (self.shape) self.dimension_labels = {} - + if dimension_labels is not None and \ len (dimension_labels) == self.number_of_dimensions: for i in range(self.number_of_dimensions): @@ -751,5 +751,9 @@ if __name__ == '__main__': print (type(volume3 + 2)) + s = [i for i in range(3 * 4 * 4)] + s = numpy.reshape(numpy.asarray(s), (3,4,4)) + sino = SinogramData( s ) + \ No newline at end of file diff --git a/Wrappers/Python/ccpi/processors.py b/Wrappers/Python/ccpi/processors.py new file mode 100755 index 0000000..f7dcf4e --- /dev/null +++ b/Wrappers/Python/ccpi/processors.py @@ -0,0 +1,164 @@ +# -*- 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 +import numpy +import h5py + +class NormalizationDataSetProcessor(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: SinogramDataSet + ''' + + def __init__(self): + kwargs = { + 'flat_field' :None, + 'dark_field' :None, + # very small number. Used when there is a division by zero + 'tolerance' : 1e-5 + } + + #DataSetProcessor.__init__(self, **kwargs) + super(NormalizationDataSetProcessor, self).__init__(**kwargs) + + 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( + [ NormalizationDataSetProcessor.normalizeProjection( + projection, flat, dark, self.tolerance) \ + for projection in projections.as_array() ] + ) + y = DataSet( a , True, + dimension_labels=projections.dimension_labels ) + return y + + +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') + + sino = SinogramData( proj ) + + + normalizer = NormalizationDataSetProcessor() + 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( + [NormalizationDataSetProcessor.normalizeProjection( p, flat, dark, 1e-5 ) + for p in proj] + ) + + print ("Numpy min {0} max {1}".format(norm1.min(), norm1.max())) \ No newline at end of file -- cgit v1.2.3