savu.plugins.reconstructions.base_recon.BaseRecon.set_centre_of_rotation() - Code Metrics - Inspection of "Merge pull request #924 from DiamondLightSource/ng..." - DiamondLightSource/Savu - Measure and Improve Code Quality continuously with Scrutinizer

Test Failed

Push — master ( f084af...daa7d2 )

unknown

created 2022-05-24 11:56 UTC

BaseRecon.set_centre_of_rotation() A

↳ Parent: savu.plugins.reconstructions.base_recon

Complexity

Conditions

Size

Total Lines	22
Code Lines	16

Duplication

Lines	0
Ratio	0 %

Importance

Changes

Metric	Value
cc	5
eloc	16
nop	4
dl	0
loc	22
rs	9.1333
c	0
b	0
f	0

# Copyright 2014 Diamond Light Source Ltd.
#
# 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.

"""
.. module:: base_recon
   :platform: Unix
   :synopsis: A base class for all reconstruction methods

.. moduleauthor:: Mark Basham <[email protected]>

"""

import math
import copy
import numpy as np
import subprocess as sp
import os
np.seterr(divide='ignore', invalid='ignore')

import savu.core.utils as cu
from savu.plugins.plugin import Plugin

MAX_OUTER_PAD = 2.1


class BaseRecon(Plugin):


    def __init__(self, name='BaseRecon'):
        super(BaseRecon, self).__init__(name)
        self.nOut = 1
        self.nIn = 1
        self.scan_dim = None
        self.rep_dim = None
        self.br_vol_shape = None
        self.frame_angles = None
        self.frame_cors = None
        self.projection_shifts = None
        self.frame_init_data = None
        self.centre = None
        self.base_pad_amount = None
        self.padding_alg = False
        self.cor_shift = 0
        self.init_vol = False
        self.cor_as_dataset = False

    def base_pre_process(self):
        in_data, out_data = self.get_datasets()
        in_pData, out_pData = self.get_plugin_datasets()
        self.pad_dim = \
            in_pData[0].get_data_dimension_by_axis_label('x', contains=True)
        in_meta_data = self.get_in_meta_data()[0]

        self.exp.log(self.name + " End")
        self.br_vol_shape = out_pData[0].get_shape()
        self.set_centre_of_rotation(in_data[0], out_data[0], in_meta_data)
        self.set_projection_shifts(in_data[0], out_data[0], in_meta_data)

        self.main_dir = in_data[0].get_data_patterns()['SINOGRAM']['main_dir']
        self.angles = in_meta_data.get('rotation_angle')

        if len(self.angles.shape) != 1:
            self.scan_dim = in_data[0].get_data_dimension_by_axis_label('scan')
        self.slice_dirs = out_data[0].get_slice_dimensions()

        shape = in_pData[0].get_shape()
        factor = self.__get_outer_pad()
        self.sino_pad = int(math.ceil(factor * shape[self.pad_dim]))

        self.sino_func, self.cor_func = self.set_function(shape)

        self.range = self.parameters['force_zero']
        self.fix_sino = self.get_sino_centre_method()

    def __get_outer_pad(self):
        pad = self.parameters['outer_pad'] if 'outer_pad' in self.parameters \
            else False
        # length of diagonal of square is side*sqrt(2)
        factor = math.sqrt(2) - 1
        if isinstance(pad, bool):
            return factor if pad is True else 0

        factor = float(pad)
        if factor > MAX_OUTER_PAD:
            factor = MAX_OUTER_PAD
            msg = 'Maximum outer_pad value is 2.1, using this instead'
            cu.user_message(msg)
        return factor

    def get_vol_shape(self):
        return self.br_vol_shape

    def set_projection_shifts(self, inData, outData, mData):
        # get experimental metadata of projection_shifts 
        if 'projection_shifts' in list(self.exp.meta_data.dict.keys()):
            self.projection_shifts = self.exp.meta_data.dict['projection_shifts']
        else:
            proj_shifts = np.zeros((inData.get_shape()[0], 2))  # initialise a 2d array of projection shifts
            self.exp.meta_data.set('projection_shifts', proj_shifts)
        outData.meta_data.set("projection_shifts", copy.deepcopy(self.projection_shifts))

    def set_centre_of_rotation(self, inData, outData, mData):
        # if cor has been passed as a dataset then do nothing
        if isinstance(self.parameters['centre_of_rotation'], str):
            return
        if 'centre_of_rotation' in list(mData.get_dictionary().keys()):
            cor = self.__set_param_from_meta_data(mData, inData, 'centre_of_rotation')
        else:
            val = self.parameters['centre_of_rotation']
            if isinstance(val, dict):
                cor = self.__polyfit_cor(val, inData)
            else:
                sdirs = inData.get_slice_dimensions()
                cor = np.ones(np.prod([inData.get_shape()[i] for i in sdirs]))
                # if centre of rotation has not been set then fix it in the
                # centre
                val = val if val != 0 else \
                    (self.get_vol_shape()[self._get_detX_dim()]) / 2.0
                cor *= val
                # mData.set('centre_of_rotation', cor) see Github ticket
        self.cor = cor
        outData.meta_data.set("centre_of_rotation", copy.deepcopy(self.cor))
        self.centre = self.cor[0]

    def populate_metadata_to_output(self, inData, outData, mData, meta_list):
        # writing into the metadata associated with the output (reconstruction)
        for meta_items in meta_list:
            outData.meta_data.set(meta_items, copy.deepcopy(mData.get(meta_items)))

        xDim = inData.get_data_dimension_by_axis_label('x', contains=True)
        det_length = inData.get_shape()[xDim]
        outData.meta_data.set("detector_x_length", copy.deepcopy(det_length))

    def __set_param_from_meta_data(self, mData, inData, meta_string):
        meta_param = mData.get(meta_string)
        sdirs = inData.get_slice_dimensions()
        total_frames = np.prod([inData.get_shape()[i] for i in sdirs])
        if total_frames > len(meta_param):
            meta_param = np.tile(meta_param, total_frames // len(meta_param))
        return meta_param

    def __polyfit_cor(self, cor_dict, inData):
        if 'detector_y' in list(inData.meta_data.get_dictionary().keys()):
            y = inData.meta_data.get('detector_y')
        else:
            yDim = inData.get_data_dimension_by_axis_label('detector_y')
            y = np.arange(inData.get_shape()[yDim])

        z = np.polyfit(list(map(int, list(cor_dict.keys()))), list(cor_dict.values()), 1)
        p = np.poly1d(z)
        cor = p(y)
        return cor

    def set_function(self, pad_shape):
        centre_pad = self.parameters['centre_pad'] if 'centre_pad' in \
            self.parameters else False
        if not centre_pad:
            def cor_func(cor):
                return cor
            if self.parameters['log']:
                sino_func = self.__make_lambda()
            else:
                sino_func = self.__make_lambda(log=False)
        else:
            def cor_func(cor):
                return cor + self.sino_pad
            if self.parameters['log']:
                sino_func = self.__make_lambda(pad=pad_shape)
            else:
                sino_func = self.__make_lambda(pad=pad_shape, log=False)
        return sino_func, cor_func

    def __make_lambda(self, log=True, pad=False):
        log_func = 'np.nan_to_num(sino)' if not log else self.parameters['log_func']
        if pad:
            pad_tuples, mode = self.__get_pad_values(pad)
            log_func = log_func.replace(
                    'sino', 'np.pad(sino, %s, "%s")' % (pad_tuples, mode))
        return eval("lambda sino: " + log_func)

    def __get_pad_values(self, pad_shape):
        mode = 'edge'
        pad_tuples = [(0, 0)] * (len(pad_shape) - 1)
        pad_tuples.insert(self.pad_dim, (self.sino_pad, self.sino_pad))
        pad_tuples = tuple(pad_tuples)
        return pad_tuples, mode

    def base_process_frames_before(self, data):
        """
        Reconstruct a single sinogram with the provided centre of rotation
        """
        sl = self.get_current_slice_list()[0]
        init = data[1] if self.init_vol else None
        angles = \
            self.angles[:, sl[self.scan_dim]] if self.scan_dim else self.angles
        angles = np.squeeze(angles)

        self.frame_angles = angles

        dim_sl = sl[self.main_dir]

        if self.cor_as_dataset:
            self.frame_cors = self.cor_func(data[len(data) - 1])
        else:
            frame_nos = \
                self.get_plugin_in_datasets()[0].get_current_frame_idx()
            a = self.cor[tuple([frame_nos])]
            self.frame_cors = self.cor_func(a)

        # for extra padded frames that make up the numbers
        if not self.frame_cors.shape:
            self.frame_cors = np.array([self.centre])

        len_data = len(np.arange(dim_sl.start, dim_sl.stop, dim_sl.step))

        missing = [self.centre] * (len(self.frame_cors) - len_data)
        self.frame_cors = np.append(self.frame_cors, missing)

        # fix to remove NaNs in the initialised image
        if init is not None:
            init[np.isnan(init)] == 0.0
        self.frame_init_data = init

        data[0] = self.fix_sino(self.sino_func(data[0]), self.frame_cors[0])
        return data

    def base_process_frames_after(self, data):
        lower_range, upper_range = self.range
        data = np.nan_to_num(data)
        if lower_range is not None:
            data[data < lower_range] = 0.0
        if upper_range is not None:
            data[data > upper_range] = 0.0
        return data

    def pad_sino(self, sino, cor):
        """  Pad the sinogram so the centre of rotation is at the centre. """
        detX = self._get_detX_dim()
        pad = self.get_centre_offset(sino, cor, detX)
        self.cor_shift = pad[0]
        pad_tuples = [(0, 0)] * (len(sino.shape) - 1)
        pad_tuples.insert(detX, tuple(pad))
        self.__set_pad_amount(max(pad))
        return np.pad(sino, tuple(pad_tuples), mode='edge')

    def _get_detX_dim(self):
        pData = self.get_plugin_in_datasets()[0]
        return pData.get_data_dimension_by_axis_label('x', contains=True)

    def get_centre_offset(self, sino, cor, detX):
        centre_pad = self.br_array_pad(cor, sino.shape[detX])
        sino_width = sino.shape[detX]
        new_width = sino_width + max(centre_pad)
        sino_pad = int(math.ceil(float(sino_width) / new_width * self.sino_pad) // 2)
        pad = np.array([sino_pad]*2) + centre_pad
        return pad

    def get_centre_shift(self, sino, cor):
        detX = self._get_detX_dim()
        return max(self.get_centre_offset(sino, self.centre, detX))

    def crop_sino(self, sino, cor):
        """  Crop the sinogram so the centre of rotation is at the centre. """
        detX = self._get_detX_dim()
        start, stop = self.br_array_pad(cor, sino.shape[detX])[::-1]
        self.cor_shift = -start
        sl = [slice(None)] * len(sino.shape)
        sl[detX] = slice(start, sino.shape[detX] - stop)
        sino = sino[tuple(sl)]
        self.set_mask(sino.shape)
        return sino

    def br_array_pad(self, ctr, nPixels):
        width = nPixels - 1.0
        alen = ctr
        blen = width - ctr
        mid = (width - 1.0) / 2.0
        shift = round(abs(blen - alen))
        p_low = 0 if (ctr > mid) else shift
        p_high = shift + 0 if (ctr > mid) else 0
        return np.array([int(p_low), int(p_high)])

    def keep_sino(self, sino, cor):
        """ No change to the sinogram """
        return sino

    def get_sino_centre_method(self):
        centre_pad = self.keep_sino
        if 'centre_pad' in list(self.parameters.keys()):
            cpad = self.parameters['centre_pad']
            if not (cpad is True or cpad is False):
                raise Exception('Unknown value for "centre_pad", please choose'
                                ' True or False.')
            centre_pad = self.pad_sino if cpad else self.crop_sino
        return centre_pad

    def __set_pad_amount(self, pad_amount):
        self.base_pad_amount = pad_amount

    def get_pad_amount(self):
        return self.base_pad_amount

    def get_fov_fraction(self, sino, cor):
        """ Get the fraction of the original FOV that can be reconstructed due\
        to offset centre """
        pData = self.get_plugin_in_datasets()[0]
        detX = pData.get_data_dimension_by_axis_label('x', contains=True)
        original_length = sino.shape[detX]
        shift = self.get_centre_shift(sino, cor)
        return (original_length - shift) / float(original_length)

    def get_reconstruction_alg(self):
        return None

    def get_angles(self):
        """ Get the angles associated with the current sinogram(s).

        :returns: Angles of the current frames.
        :rtype: np.ndarray
        """
        return self.frame_angles

    def get_proj_shifts(self):
        """ Get the 2D (X-Y) shifts associated with every projection frame

        :returns: projecton shifts for the current frames.
        :rtype: np.ndarray
        """
        return self.projection_shifts

    def get_cors(self):
        """
        Get the centre of rotations associated with the current sinogram(s).

        :returns: Centre of rotation values for the current frames.
        :rtype: np.ndarray
        """
        return self.frame_cors + self.cor_shift

    def set_mask(self, shape):
        pass

    def get_initial_data(self):
        """
        Get the initial data (if it is exists) associated with the current \
        sinogram(s).

        :returns: The section of the initialisation data associated with the \
            current frames.
        :rtype: np.ndarray or None
        """
        return self.frame_init_data

    def get_frame_params(self):
        params = [self.get_cors(), self.get_angles(), self.get_vol_shape(),
                  self.get_initial_data()]
        return params

    def get_frame_shifts(self):
        return self.get_proj_shifts()

    def setup(self):
        in_dataset, out_dataset = self.get_datasets()
        # reduce the data as per data_subset parameter
        self.preview_flag = \
            self.set_preview(in_dataset[0], self.parameters['preview'])

        self._set_volume_dimensions(in_dataset[0])
        axis_labels = self._get_axis_labels(in_dataset[0])
        shape = self._get_shape(in_dataset[0])

        # output dataset
        out_dataset[0].create_dataset(axis_labels=axis_labels, shape=shape)
        out_dataset[0].add_volume_patterns(*self._get_volume_dimensions())

        # set information relating to the plugin data
        in_pData, out_pData = self.get_plugin_datasets()

        self.init_vol = 1 if 'init_vol' in list(self.parameters.keys()) and\
            self.parameters['init_vol'] else 0
        self.cor_as_dataset = 1 if isinstance(
            self.parameters['centre_of_rotation'], str) else 0

        for i in range(len(in_dataset) - self.init_vol - self.cor_as_dataset):
            in_pData[i].plugin_data_setup('SINOGRAM', self.get_max_frames(),
                                          slice_axis=self.get_slice_axis())
            idx = 1

        # initial volume dataset
        if self.init_vol:
#            from savu.data.data_structures.data_types import Replicate
#            if self.rep_dim:
#                in_dataset[idx].data = Replicate(
#                    in_dataset[idx], in_dataset[0].get_shape(self.rep_dim))
            in_pData[1].plugin_data_setup('VOLUME_XZ', self.get_max_frames())
            idx += 1


        # cor dataset
        if self.cor_as_dataset:
            self.cor_as_dataset = True
            in_pData[idx].plugin_data_setup('METADATA', self.get_max_frames())

        # set pattern_name and nframes to process for all datasets
        out_pData[0].plugin_data_setup('VOLUME_XZ', self.get_max_frames())

        meta_list = ['rotation_angle']  # metadata list to populate
        in_meta_data = self.get_in_meta_data()[0]

        if 'projection_shifts' in list(self.exp.meta_data.dict.keys()):
            self.projection_shifts = self.exp.meta_data.dict['projection_shifts']
        else:
            self.projection_shifts = np.zeros((in_dataset[0].get_shape()[self.volX], 2))  # initialise a 2d array of projection shifts
            self.exp.meta_data.set('projection_shifts', copy.deepcopy(self.projection_shifts))

        out_dataset[0].meta_data.set("projection_shifts", copy.deepcopy(self.projection_shifts))
        self.populate_metadata_to_output(in_dataset[0], out_dataset[0], in_meta_data, meta_list)

    def _get_axis_labels(self, in_dataset):
        """
        Get the new axis labels for the output dataset - this is now a volume.

        Parameters
        ----------
        in_dataset : :class:`savu.data.data_structures.data.Data`
            The input dataset to the plugin.

        Returns
        -------
        labels : dict
            The axis labels for the dataset that is output from the plugin.

        """
        labels = in_dataset.data_info.get('axis_labels')[0]
        volX, volY, volZ = self._get_volume_dimensions()
        labels = [str(volX) + '.voxel_x.voxels', str(volZ) + '.voxel_z.voxels']
        if volY:
            labels.append(str(volY) + '.voxel_y.voxels')
        labels = {in_dataset: labels}
        return labels

    def _set_volume_dimensions(self, data):
        """
        Map the input dimensions to the output volume dimensions

        Parameters
        ----------
        in_dataset : :class:`savu.data.data_structures.data.Data`
            The input dataset to the plugin.
        """
        data._finalise_patterns()
        self.volX = data.get_data_dimension_by_axis_label("rotation_angle")
        self.volZ = data.get_data_dimension_by_axis_label("x", contains=True)
        self.volY = data.get_data_dimension_by_axis_label(
            "y", contains=True, exists=True)

    def _get_volume_dimensions(self):
        return self.volX, self.volY, self.volZ

    def _get_shape(self, in_dataset):
        shape = list(in_dataset.get_shape())
        volX, volY, volZ = self._get_volume_dimensions()

        if self.parameters['vol_shape'] in ('auto', 'fixed'):
            shape[volX] = shape[volZ]
        else:
            shape[volX] = self.parameters['vol_shape']
            shape[volZ] = self.parameters['vol_shape']

        if 'resolution' in self.parameters.keys():
            shape[volX] = int(shape[volX] // self.parameters['resolution'])
            shape[volZ] = int(shape[volZ] // self.parameters['resolution'])
        return tuple(shape)

    def get_max_frames(self):
        """
        Number of data frames to pass to each instance of process_frames func

        Returns
        -------
        str or int
            "single", "multiple" or integer (only to be used if the number of
                                             frames MUST be fixed.)
        """
        if self._max_frames:
            frames_max = self._max_frames
        else:
            frames_max = 'multiple'
        return frames_max
        return 'multiple'

    def get_slice_axis(self):
        """
        Fix the fastest changing slice dimension

        Returns
        -------
        str or None
            str should be the axis_label corresponding to the fastest changing
            dimension

        """
        return None

    def nInput_datasets(self):
        nIn = 1
        if 'init_vol' in self.parameters.keys() and \
                self.parameters['init_vol']:
            if len(self.parameters['init_vol'].split('.')) == 3:
                name, temp, self.rep_dim = self.parameters['init_vol']
                self.parameters['init_vol'] = name
            nIn += 1
            self.parameters['in_datasets'].append(self.parameters['init_vol'])
        if isinstance(self.parameters['centre_of_rotation'], str):
            self.parameters['in_datasets'].append(
                self.parameters['centre_of_rotation'])
            nIn += 1
        return nIn

    def nOutput_datasets(self):
        return self.nOut

    def reconstruct_pre_process(self):
        """
        Should be overridden to perform pre-processing in a child class
        """
        pass

    def executive_summary(self):
        summary = []
        if not self.preview_flag:
            summary.append(("WARNING: Ignoring preview parameters as a preview"
                            " has already been applied to the data."))
        if len(summary) > 0:
            return summary
        return ["Nothing to Report"]

    def get_gpu_memory(self):
        command = "nvidia-smi --query-gpu=memory.free --format=csv"
        memory_free_info = sp.check_output(command.split()).decode('ascii').split('\n')[:-1][1:]
        memory_free_values = [int(x.split()[0]) for i, x in enumerate(memory_free_info)]
        return memory_free_values

    def _handle_case_1(self, in_text):
        final_idx = eval(in_text)
        return final_idx

    def _handle_case_2(self, in_text):
        pos1 = [idx - 2 for idx, char in enumerate(in_text) if char == "."]
        pos2 = [idx for idx, char in enumerate(in_text) if char == ")"]
        list_idx1 = []
        for i, pos in enumerate(pos1):
            cmd = in_text[pos:pos2[i] + 1]
            list_idx1.extend(eval(cmd))
        for i, pos in enumerate(pos1):
            tmp = "-1"
            for j in range(2, 1 + pos2[i] - pos):
                tmp = tmp + " "
            in_text = in_text.replace(in_text[pos:pos2[i] + 1], tmp)
        text = in_text.replace("[", "")
        text = text.replace("]", "")
        out_text = text.split(",")
        list_idx = [int(x) for x in out_text]
        list_idx = np.asarray(list_idx)
        list_idx = list_idx[list_idx != -1]
        final_idx = np.sort(np.concatenate((np.asarray(list_idx1), list_idx)))
        return np.unique(np.int16(final_idx))

    def _handle_case_3(self, in_text):
        text = in_text.replace("[", "")
        text = text.replace("]", "")
        out_text = text.split(",")
        final_idx = []
        for x in out_text:
            try:
                num = int(x)
                final_idx.append(num)
            except ValueError:
                pass
        if final_idx:
            final_idx = np.unique(np.sort(final_idx))
        else:
            final_idx = None
        return final_idx

    def get_skipping_indices(self, in_text):
        if isinstance(in_text, str):
            if "np." in in_text:
                final_idx = self._handle_case_1(in_text)
                if not isinstance(final_idx, np.ndarray):
                    final_idx = self._handle_case_2(in_text)
                else:
                    final_idx = np.unique(np.sort(np.ndarray.flatten(final_idx)))
            else:
                final_idx = self._handle_case_3(in_text)
            if not isinstance(final_idx, np.ndarray):
                raise ValueError("Incorrect syntax!!!")
        elif isinstance(in_text, list):
            final_idx = np.unique(np.sort(np.asarray(in_text)))
        else:
            final_idx = None
        return final_idx


1			# Copyright 2014 Diamond Light Source Ltd.
2			#
3			# Licensed under the Apache License, Version 2.0 (the "License");
4			# you may not use this file except in compliance with the License.
5			# You may obtain a copy of the License at
6			#
7			# http://www.apache.org/licenses/LICENSE-2.0
8			#
9			# Unless required by applicable law or agreed to in writing, software
10			# distributed under the License is distributed on an "AS IS" BASIS,
11			# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
12			# See the License for the specific language governing permissions and
13			# limitations under the License.
14
15			"""
16			.. module:: base_recon
17			:platform: Unix
18			:synopsis: A base class for all reconstruction methods
19
20			.. moduleauthor:: Mark Basham <[email protected]>
21
22			"""
23
24			import math
25			import copy
26			import numpy as np
27			import subprocess as sp
28			import os
29			np.seterr(divide='ignore', invalid='ignore')
30
31			import savu.core.utils as cu
32			from savu.plugins.plugin import Plugin
33
34			MAX_OUTER_PAD = 2.1
35
36
37			class BaseRecon(Plugin):
38
39
40			def __init__(self, name='BaseRecon'):
41			super(BaseRecon, self).__init__(name)
42			self.nOut = 1
43			self.nIn = 1
44			self.scan_dim = None
45			self.rep_dim = None
46			self.br_vol_shape = None
47			self.frame_angles = None
48			self.frame_cors = None
49			self.projection_shifts = None
50			self.frame_init_data = None
51			self.centre = None
52			self.base_pad_amount = None
53			self.padding_alg = False
54			self.cor_shift = 0
55			self.init_vol = False
56			self.cor_as_dataset = False
57
58			def base_pre_process(self):
59			in_data, out_data = self.get_datasets()
60			in_pData, out_pData = self.get_plugin_datasets()
61			self.pad_dim = \
62			in_pData[0].get_data_dimension_by_axis_label('x', contains=True)
63			in_meta_data = self.get_in_meta_data()[0]
64
65			self.exp.log(self.name + " End")
66			self.br_vol_shape = out_pData[0].get_shape()
67			self.set_centre_of_rotation(in_data[0], out_data[0], in_meta_data)
68			self.set_projection_shifts(in_data[0], out_data[0], in_meta_data)
69
70			self.main_dir = in_data[0].get_data_patterns()['SINOGRAM']['main_dir']
71			self.angles = in_meta_data.get('rotation_angle')
72
73			if len(self.angles.shape) != 1:
74			self.scan_dim = in_data[0].get_data_dimension_by_axis_label('scan')
75			self.slice_dirs = out_data[0].get_slice_dimensions()
76
77			shape = in_pData[0].get_shape()
78			factor = self.__get_outer_pad()
79			self.sino_pad = int(math.ceil(factor * shape[self.pad_dim]))
80
81			self.sino_func, self.cor_func = self.set_function(shape)
82
83			self.range = self.parameters['force_zero']
84			self.fix_sino = self.get_sino_centre_method()
85
86			def __get_outer_pad(self):
87			pad = self.parameters['outer_pad'] if 'outer_pad' in self.parameters \
88			else False
89			# length of diagonal of square is side*sqrt(2)
90			factor = math.sqrt(2) - 1
91			if isinstance(pad, bool):
92			return factor if pad is True else 0
93
94			factor = float(pad)
95			if factor > MAX_OUTER_PAD:
96			factor = MAX_OUTER_PAD
97			msg = 'Maximum outer_pad value is 2.1, using this instead'
98			cu.user_message(msg)
99			return factor
100
101			def get_vol_shape(self):
102			return self.br_vol_shape
103
104			def set_projection_shifts(self, inData, outData, mData):
105			# get experimental metadata of projection_shifts
106			if 'projection_shifts' in list(self.exp.meta_data.dict.keys()):
107			self.projection_shifts = self.exp.meta_data.dict['projection_shifts']
108			else:
109			proj_shifts = np.zeros((inData.get_shape()[0], 2)) # initialise a 2d array of projection shifts
110			self.exp.meta_data.set('projection_shifts', proj_shifts)
111			outData.meta_data.set("projection_shifts", copy.deepcopy(self.projection_shifts))
112
113			def set_centre_of_rotation(self, inData, outData, mData):
114			# if cor has been passed as a dataset then do nothing
115			if isinstance(self.parameters['centre_of_rotation'], str):
116			return
117			if 'centre_of_rotation' in list(mData.get_dictionary().keys()):
118			cor = self.__set_param_from_meta_data(mData, inData, 'centre_of_rotation')
119			else:
120			val = self.parameters['centre_of_rotation']
121			if isinstance(val, dict):
122			cor = self.__polyfit_cor(val, inData)
123			else:
124			sdirs = inData.get_slice_dimensions()
125			cor = np.ones(np.prod([inData.get_shape()[i] for i in sdirs]))
126			# if centre of rotation has not been set then fix it in the
127			# centre
128			val = val if val != 0 else \
129			(self.get_vol_shape()[self._get_detX_dim()]) / 2.0
130			cor *= val
131			# mData.set('centre_of_rotation', cor) see Github ticket
132			self.cor = cor
133			outData.meta_data.set("centre_of_rotation", copy.deepcopy(self.cor))
134			self.centre = self.cor[0]
135
136			def populate_metadata_to_output(self, inData, outData, mData, meta_list):
137			# writing into the metadata associated with the output (reconstruction)
138			for meta_items in meta_list:
139			outData.meta_data.set(meta_items, copy.deepcopy(mData.get(meta_items)))
140
141			xDim = inData.get_data_dimension_by_axis_label('x', contains=True)
142			det_length = inData.get_shape()[xDim]
143			outData.meta_data.set("detector_x_length", copy.deepcopy(det_length))
144
145			def __set_param_from_meta_data(self, mData, inData, meta_string):
146			meta_param = mData.get(meta_string)
147			sdirs = inData.get_slice_dimensions()
148			total_frames = np.prod([inData.get_shape()[i] for i in sdirs])
149			if total_frames > len(meta_param):
150			meta_param = np.tile(meta_param, total_frames // len(meta_param))
151			return meta_param
152
153			def __polyfit_cor(self, cor_dict, inData):
154			if 'detector_y' in list(inData.meta_data.get_dictionary().keys()):
155			y = inData.meta_data.get('detector_y')
156			else:
157			yDim = inData.get_data_dimension_by_axis_label('detector_y')
158			y = np.arange(inData.get_shape()[yDim])
159
160			z = np.polyfit(list(map(int, list(cor_dict.keys()))), list(cor_dict.values()), 1)
161			p = np.poly1d(z)
162			cor = p(y)
163			return cor
164
165			def set_function(self, pad_shape):
166			centre_pad = self.parameters['centre_pad'] if 'centre_pad' in \
167			self.parameters else False
168			if not centre_pad:
169			def cor_func(cor):
170			return cor
171			if self.parameters['log']:
172			sino_func = self.__make_lambda()
173			else:
174			sino_func = self.__make_lambda(log=False)
175			else:
176			def cor_func(cor):
177			return cor + self.sino_pad
178			if self.parameters['log']:
179			sino_func = self.__make_lambda(pad=pad_shape)
180			else:
181			sino_func = self.__make_lambda(pad=pad_shape, log=False)
182			return sino_func, cor_func
183
184			def __make_lambda(self, log=True, pad=False):
185			log_func = 'np.nan_to_num(sino)' if not log else self.parameters['log_func']
186			if pad:
187			pad_tuples, mode = self.__get_pad_values(pad)
188			log_func = log_func.replace(
189			'sino', 'np.pad(sino, %s, "%s")' % (pad_tuples, mode))
190			return eval("lambda sino: " + log_func)
191
192			def __get_pad_values(self, pad_shape):
193			mode = 'edge'
194			pad_tuples = [(0, 0)] * (len(pad_shape) - 1)
195			pad_tuples.insert(self.pad_dim, (self.sino_pad, self.sino_pad))
196			pad_tuples = tuple(pad_tuples)
197			return pad_tuples, mode
198
199			def base_process_frames_before(self, data):
200			"""
201			Reconstruct a single sinogram with the provided centre of rotation
202			"""
203			sl = self.get_current_slice_list()[0]
204			init = data[1] if self.init_vol else None
205			angles = \
206			self.angles[:, sl[self.scan_dim]] if self.scan_dim else self.angles
207			angles = np.squeeze(angles)
208
209			self.frame_angles = angles
210
211			dim_sl = sl[self.main_dir]
212
213			if self.cor_as_dataset:
214			self.frame_cors = self.cor_func(data[len(data) - 1])
215			else:
216			frame_nos = \
217			self.get_plugin_in_datasets()[0].get_current_frame_idx()
218			a = self.cor[tuple([frame_nos])]
219			self.frame_cors = self.cor_func(a)
220
221			# for extra padded frames that make up the numbers
222			if not self.frame_cors.shape:
223			self.frame_cors = np.array([self.centre])
224
225			len_data = len(np.arange(dim_sl.start, dim_sl.stop, dim_sl.step))
226
227			missing = [self.centre] * (len(self.frame_cors) - len_data)
228			self.frame_cors = np.append(self.frame_cors, missing)
229
230			# fix to remove NaNs in the initialised image
231			if init is not None:
232			init[np.isnan(init)] == 0.0
233			self.frame_init_data = init
234
235			data[0] = self.fix_sino(self.sino_func(data[0]), self.frame_cors[0])
236			return data
237
238			def base_process_frames_after(self, data):
239			lower_range, upper_range = self.range
240			data = np.nan_to_num(data)
241			if lower_range is not None:
242			data[data < lower_range] = 0.0
243			if upper_range is not None:
244			data[data > upper_range] = 0.0
245			return data
246
247			def pad_sino(self, sino, cor):
248			""" Pad the sinogram so the centre of rotation is at the centre. """
249			detX = self._get_detX_dim()
250			pad = self.get_centre_offset(sino, cor, detX)
251			self.cor_shift = pad[0]
252			pad_tuples = [(0, 0)] * (len(sino.shape) - 1)
253			pad_tuples.insert(detX, tuple(pad))
254			self.__set_pad_amount(max(pad))
255			return np.pad(sino, tuple(pad_tuples), mode='edge')
256
257			def _get_detX_dim(self):
258			pData = self.get_plugin_in_datasets()[0]
259			return pData.get_data_dimension_by_axis_label('x', contains=True)
260
261			def get_centre_offset(self, sino, cor, detX):
262			centre_pad = self.br_array_pad(cor, sino.shape[detX])
263			sino_width = sino.shape[detX]
264			new_width = sino_width + max(centre_pad)
265			sino_pad = int(math.ceil(float(sino_width) / new_width * self.sino_pad) // 2)
266			pad = np.array([sino_pad]*2) + centre_pad
267			return pad
268
269			def get_centre_shift(self, sino, cor):
270			detX = self._get_detX_dim()
271			return max(self.get_centre_offset(sino, self.centre, detX))
272
273			def crop_sino(self, sino, cor):
274			""" Crop the sinogram so the centre of rotation is at the centre. """
275			detX = self._get_detX_dim()
276			start, stop = self.br_array_pad(cor, sino.shape[detX])[::-1]
277			self.cor_shift = -start
278			sl = [slice(None)] * len(sino.shape)
279			sl[detX] = slice(start, sino.shape[detX] - stop)
280			sino = sino[tuple(sl)]
281			self.set_mask(sino.shape)
282			return sino
283
284			def br_array_pad(self, ctr, nPixels):
285			width = nPixels - 1.0
286			alen = ctr
287			blen = width - ctr
288			mid = (width - 1.0) / 2.0
289			shift = round(abs(blen - alen))
290			p_low = 0 if (ctr > mid) else shift
291			p_high = shift + 0 if (ctr > mid) else 0
292			return np.array([int(p_low), int(p_high)])
293
294			def keep_sino(self, sino, cor):
295			""" No change to the sinogram """
296			return sino
297
298			def get_sino_centre_method(self):
299			centre_pad = self.keep_sino
300			if 'centre_pad' in list(self.parameters.keys()):
301			cpad = self.parameters['centre_pad']
302			if not (cpad is True or cpad is False):
303			raise Exception('Unknown value for "centre_pad", please choose'
304			' True or False.')
305			centre_pad = self.pad_sino if cpad else self.crop_sino
306			return centre_pad
307
308			def __set_pad_amount(self, pad_amount):
309			self.base_pad_amount = pad_amount
310
311			def get_pad_amount(self):
312			return self.base_pad_amount
313
314			def get_fov_fraction(self, sino, cor):
315			""" Get the fraction of the original FOV that can be reconstructed due\
316			to offset centre """
317			pData = self.get_plugin_in_datasets()[0]
318			detX = pData.get_data_dimension_by_axis_label('x', contains=True)
319			original_length = sino.shape[detX]
320			shift = self.get_centre_shift(sino, cor)
321			return (original_length - shift) / float(original_length)
322
323			def get_reconstruction_alg(self):
324			return None
325
326			def get_angles(self):
327			""" Get the angles associated with the current sinogram(s).
328
329			:returns: Angles of the current frames.
330			:rtype: np.ndarray
331			"""
332			return self.frame_angles
333
334			def get_proj_shifts(self):
335			""" Get the 2D (X-Y) shifts associated with every projection frame
336
337			:returns: projecton shifts for the current frames.
338			:rtype: np.ndarray
339			"""
340			return self.projection_shifts
341
342			def get_cors(self):
343			"""
344			Get the centre of rotations associated with the current sinogram(s).
345
346			:returns: Centre of rotation values for the current frames.
347			:rtype: np.ndarray
348			"""
349			return self.frame_cors + self.cor_shift
350
351			def set_mask(self, shape):
352			pass
353
354			def get_initial_data(self):
355			"""
356			Get the initial data (if it is exists) associated with the current \
357			sinogram(s).
358
359			:returns: The section of the initialisation data associated with the \
360			current frames.
361			:rtype: np.ndarray or None
362			"""
363			return self.frame_init_data
364
365			def get_frame_params(self):
366			params = [self.get_cors(), self.get_angles(), self.get_vol_shape(),
367			self.get_initial_data()]
368			return params
369
370			def get_frame_shifts(self):
371			return self.get_proj_shifts()
372
373			def setup(self):
374			in_dataset, out_dataset = self.get_datasets()
375			# reduce the data as per data_subset parameter
376			self.preview_flag = \
377			self.set_preview(in_dataset[0], self.parameters['preview'])
378
379			self._set_volume_dimensions(in_dataset[0])
380			axis_labels = self._get_axis_labels(in_dataset[0])
381			shape = self._get_shape(in_dataset[0])
382
383			# output dataset
384			out_dataset[0].create_dataset(axis_labels=axis_labels, shape=shape)
385			out_dataset[0].add_volume_patterns(*self._get_volume_dimensions())
386
387			# set information relating to the plugin data
388			in_pData, out_pData = self.get_plugin_datasets()
389
390			self.init_vol = 1 if 'init_vol' in list(self.parameters.keys()) and\
391			self.parameters['init_vol'] else 0
392			self.cor_as_dataset = 1 if isinstance(
393			self.parameters['centre_of_rotation'], str) else 0
394
395			for i in range(len(in_dataset) - self.init_vol - self.cor_as_dataset):
396			in_pData[i].plugin_data_setup('SINOGRAM', self.get_max_frames(),
397			slice_axis=self.get_slice_axis())
398			idx = 1
399
400			# initial volume dataset
401			if self.init_vol:
402			# from savu.data.data_structures.data_types import Replicate
403			# if self.rep_dim:
404			# in_dataset[idx].data = Replicate(
405			# in_dataset[idx], in_dataset[0].get_shape(self.rep_dim))
406			in_pData[1].plugin_data_setup('VOLUME_XZ', self.get_max_frames())
407			idx += 1
			0 ignored issues – show introduced 2021-02-25 17:14 UTC by Report Bug Copy Issue Report The variable `idx` does not seem to be defined in case the `for` loop on line `395` is not entered. Are you sure this can never be the case? Loading history...
408
409			# cor dataset
410			if self.cor_as_dataset:
411			self.cor_as_dataset = True
412			in_pData[idx].plugin_data_setup('METADATA', self.get_max_frames())
413
414			# set pattern_name and nframes to process for all datasets
415			out_pData[0].plugin_data_setup('VOLUME_XZ', self.get_max_frames())
416
417			meta_list = ['rotation_angle'] # metadata list to populate
418			in_meta_data = self.get_in_meta_data()[0]
419
420			if 'projection_shifts' in list(self.exp.meta_data.dict.keys()):
421			self.projection_shifts = self.exp.meta_data.dict['projection_shifts']
422			else:
423			self.projection_shifts = np.zeros((in_dataset[0].get_shape()[self.volX], 2)) # initialise a 2d array of projection shifts
424			self.exp.meta_data.set('projection_shifts', copy.deepcopy(self.projection_shifts))
425
426			out_dataset[0].meta_data.set("projection_shifts", copy.deepcopy(self.projection_shifts))
427			self.populate_metadata_to_output(in_dataset[0], out_dataset[0], in_meta_data, meta_list)
428
429			def _get_axis_labels(self, in_dataset):
430			"""
431			Get the new axis labels for the output dataset - this is now a volume.
432
433			Parameters
434			----------
435			in_dataset : :class:`savu.data.data_structures.data.Data`
436			The input dataset to the plugin.
437
438			Returns
439			-------
440			labels : dict
441			The axis labels for the dataset that is output from the plugin.
442
443			"""
444			labels = in_dataset.data_info.get('axis_labels')[0]
445			volX, volY, volZ = self._get_volume_dimensions()
446			labels = [str(volX) + '.voxel_x.voxels', str(volZ) + '.voxel_z.voxels']
447			if volY:
448			labels.append(str(volY) + '.voxel_y.voxels')
449			labels = {in_dataset: labels}
450			return labels
451
452			def _set_volume_dimensions(self, data):
453			"""
454			Map the input dimensions to the output volume dimensions
455
456			Parameters
457			----------
458			in_dataset : :class:`savu.data.data_structures.data.Data`
459			The input dataset to the plugin.
460			"""
461			data._finalise_patterns()
462			self.volX = data.get_data_dimension_by_axis_label("rotation_angle")
463			self.volZ = data.get_data_dimension_by_axis_label("x", contains=True)
464			self.volY = data.get_data_dimension_by_axis_label(
465			"y", contains=True, exists=True)
466
467			def _get_volume_dimensions(self):
468			return self.volX, self.volY, self.volZ
469
470			def _get_shape(self, in_dataset):
471			shape = list(in_dataset.get_shape())
472			volX, volY, volZ = self._get_volume_dimensions()
473
474			if self.parameters['vol_shape'] in ('auto', 'fixed'):
475			shape[volX] = shape[volZ]
476			else:
477			shape[volX] = self.parameters['vol_shape']
478			shape[volZ] = self.parameters['vol_shape']
479
480			if 'resolution' in self.parameters.keys():
481			shape[volX] = int(shape[volX] // self.parameters['resolution'])
482			shape[volZ] = int(shape[volZ] // self.parameters['resolution'])
483			return tuple(shape)
484
485			def get_max_frames(self):
486			"""
487			Number of data frames to pass to each instance of process_frames func
488
489			Returns
490			-------
491			str or int
492			"single", "multiple" or integer (only to be used if the number of
493			frames MUST be fixed.)
494			"""
495			if self._max_frames:
496			frames_max = self._max_frames
497			else:
498			frames_max = 'multiple'
499			return frames_max
500			return 'multiple'
501
502			def get_slice_axis(self):
503			"""
504			Fix the fastest changing slice dimension
505
506			Returns
507			-------
508			str or None
509			str should be the axis_label corresponding to the fastest changing
510			dimension
511
512			"""
513			return None
514
515			def nInput_datasets(self):
516			nIn = 1
517			if 'init_vol' in self.parameters.keys() and \
518			self.parameters['init_vol']:
519			if len(self.parameters['init_vol'].split('.')) == 3:
520			name, temp, self.rep_dim = self.parameters['init_vol']
521			self.parameters['init_vol'] = name
522			nIn += 1
523			self.parameters['in_datasets'].append(self.parameters['init_vol'])
524			if isinstance(self.parameters['centre_of_rotation'], str):
525			self.parameters['in_datasets'].append(
526			self.parameters['centre_of_rotation'])
527			nIn += 1
528			return nIn
529
530			def nOutput_datasets(self):
531			return self.nOut
532
533			def reconstruct_pre_process(self):
534			"""
535			Should be overridden to perform pre-processing in a child class
536			"""
537			pass
538
539			def executive_summary(self):
540			summary = []
541			if not self.preview_flag:
542			summary.append(("WARNING: Ignoring preview parameters as a preview"
543			" has already been applied to the data."))
544			if len(summary) > 0:
545			return summary
546			return ["Nothing to Report"]
547
548			def get_gpu_memory(self):
549			command = "nvidia-smi --query-gpu=memory.free --format=csv"
550			memory_free_info = sp.check_output(command.split()).decode('ascii').split('\n')[:-1][1:]
551			memory_free_values = [int(x.split()[0]) for i, x in enumerate(memory_free_info)]
552			return memory_free_values
553
554			def _handle_case_1(self, in_text):
555			final_idx = eval(in_text)
556			return final_idx
557
558			def _handle_case_2(self, in_text):
559			pos1 = [idx - 2 for idx, char in enumerate(in_text) if char == "."]
560			pos2 = [idx for idx, char in enumerate(in_text) if char == ")"]
561			list_idx1 = []
562			for i, pos in enumerate(pos1):
563			cmd = in_text[pos:pos2[i] + 1]
564			list_idx1.extend(eval(cmd))
565			for i, pos in enumerate(pos1):
566			tmp = "-1"
567			for j in range(2, 1 + pos2[i] - pos):
568			tmp = tmp + " "
569			in_text = in_text.replace(in_text[pos:pos2[i] + 1], tmp)
570			text = in_text.replace("[", "")
571			text = text.replace("]", "")
572			out_text = text.split(",")
573			list_idx = [int(x) for x in out_text]
574			list_idx = np.asarray(list_idx)
575			list_idx = list_idx[list_idx != -1]
576			final_idx = np.sort(np.concatenate((np.asarray(list_idx1), list_idx)))
577			return np.unique(np.int16(final_idx))
578
579			def _handle_case_3(self, in_text):
580			text = in_text.replace("[", "")
581			text = text.replace("]", "")
582			out_text = text.split(",")
583			final_idx = []
584			for x in out_text:
585			try:
586			num = int(x)
587			final_idx.append(num)
588			except ValueError:
589			pass
590			if final_idx:
591			final_idx = np.unique(np.sort(final_idx))
592			else:
593			final_idx = None
594			return final_idx
595
596			def get_skipping_indices(self, in_text):
597			if isinstance(in_text, str):
598			if "np." in in_text:
599			final_idx = self._handle_case_1(in_text)
600			if not isinstance(final_idx, np.ndarray):
601			final_idx = self._handle_case_2(in_text)
602			else:
603			final_idx = np.unique(np.sort(np.ndarray.flatten(final_idx)))
604			else:
605			final_idx = self._handle_case_3(in_text)
606			if not isinstance(final_idx, np.ndarray):
607			raise ValueError("Incorrect syntax!!!")
608			elif isinstance(in_text, list):
609			final_idx = np.unique(np.sort(np.asarray(in_text)))
610			else:
611			final_idx = None
612			return final_idx
613

DiamondLightSource / Savu

Push — master ( f084af...daa7d2 )

BaseRecon.set_centre_of_rotation() A

Complexity

Size

Duplication

Importance

Duplication Side-by-Side

Filter issues like