tomo: tomo.py comparison

comparison tomo.py @ 0:cb1b0d757704 draft

"planemo upload for repository https://github.com/rolfverberg/galaxytools commit 2da52c7db6def807073a1d437a00e0e2a8e7e72e"

author	rv43
date	Tue, 29 Mar 2022 16:10:16 +0000
parents
children	e4778148df6b

comparison

equal deleted inserted replaced

--1:000000000000
+:cb1b0d757704
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+"""
+Created on Fri Dec 10 09:54:37 2021
+@author: rv43
+"""
+import logging
+import os
+import sys
+import getopt
+import re
+import io
+try:
+import pyinputplus as pyip
+except:
+pass
+import numpy as np
+import numexpr as ne
+import multiprocessing as mp
+import scipy.ndimage as spi
+import tomopy
+from time import time
+from skimage.transform import iradon
+from skimage.restoration import denoise_tv_chambolle
+import msnc_tools as msnc
+class set_numexpr_threads:
+def __init__(self, nthreads):
+cpu_count = mp.cpu_count()
+if nthreads is None or nthreads > cpu_count:
+self.n = cpu_count
+else:
+self.n = nthreads
+def __enter__(self):
+self.oldn = ne.set_num_threads(self.n)
+def __exit__(self, exc_type, exc_value, traceback):
+ne.set_num_threads(self.oldn)
+class ConfigTomo(msnc.Config):
+"""Class for processing a config file.
+"""
+def __init__(self, config_file=None, config_dict=None):
+super().__init__(config_file, config_dict)
+def _validate_txt(self):
+"""Returns False if any required config parameter is illegal or missing.
+"""
+is_valid = True
+# Check for required first-level keys
+pars_required = ['tdf_data_path', 'tbf_data_path', 'detector_id']
+pars_missing = []
+is_valid = super().validate(pars_required, pars_missing)
+if len(pars_missing) > 0:
+logging.error(f'Missing item(s) in config file: {", ".join(pars_missing)}')
+self.detector_id = self.config.get('detector_id')
+# Find tomography dark field images file/folder
+self.tdf_data_path = self.config.get('tdf_data_path')
+# Find tomography bright field images file/folder
+self.tbf_data_path = self.config.get('tbf_data_path')
+# Check number of tomography image stacks
+self.num_tomo_stacks = self.config.get('num_tomo_stacks', 1)
+if not msnc.is_int(self.num_tomo_stacks, 1):
+self.num_tomo_stacks = None
+msnc.illegal_value('num_tomo_stacks', self.num_tomo_stacks, 'config file')
+return False
+logging.info(f'num_tomo_stacks = {self.num_tomo_stacks}')
+# Find tomography images file/folders and stack parameters
+tomo_data_paths_indices = sorted({key:value for key,value in self.config.items()
+if 'tomo_data_path' in key}.items())
+if len(tomo_data_paths_indices) != self.num_tomo_stacks:
+logging.error(f'Incorrect number of tomography data path names in config file')
+is_valid = False
+self.tomo_data_paths = [tomo_data_paths_indices[i][1] for i in range(self.num_tomo_stacks)]
+self.tomo_data_indices = [msnc.get_trailing_int(tomo_data_paths_indices[i][0])
+if msnc.get_trailing_int(tomo_data_paths_indices[i][0]) else None
+for i in range(self.num_tomo_stacks)]
+tomo_ref_height_indices = sorted({key:value for key,value in self.config.items()
+if 'z_pos' in key}.items())
+if self.num_tomo_stacks > 1 and len(tomo_ref_height_indices) != self.num_tomo_stacks:
+logging.error(f'Incorrect number of tomography reference heights in config file')
+is_valid = False
+if len(tomo_ref_height_indices):
+self.tomo_ref_heights = [
+tomo_ref_height_indices[i][1] for i in range(self.num_tomo_stacks)]
+else:
+self.tomo_ref_heights = [0.0]*self.num_tomo_stacks
+# Check tomo angle (theta) range
+self.start_theta = self.config.get('start_theta', 0.)
+if not msnc.is_num(self.start_theta, 0.):
+msnc.illegal_value('start_theta', self.start_theta, 'config file')
+is_valid = False
+logging.debug(f'start_theta = {self.start_theta}')
+self.end_theta = self.config.get('end_theta', 180.)
+if not msnc.is_num(self.end_theta, self.start_theta):
+msnc.illegal_value('end_theta', self.end_theta, 'config file')
+is_valid = False
+logging.debug(f'end_theta = {self.end_theta}')
+self.num_thetas = self.config.get('num_thetas')
+if not (self.num_thetas is None or msnc.is_int(self.num_thetas, 1)):
+msnc.illegal_value('num_thetas', self.num_thetas, 'config file')
+self.num_thetas = None
+is_valid = False
+logging.debug(f'num_thetas = {self.num_thetas}')
+return is_valid
+def _validate_yaml(self):
+"""Returns False if any required config parameter is illegal or missing.
+"""
+is_valid = True
+# Check for required first-level keys
+pars_required = ['dark_field', 'bright_field', 'stack_info', 'detector']
+pars_missing = []
+is_valid = super().validate(pars_required, pars_missing)
+if len(pars_missing) > 0:
+logging.error(f'Missing item(s) in config file: {", ".join(pars_missing)}')
+self.detector_id = self.config['detector'].get('id')
+# Find tomography dark field images file/folder
+self.tdf_data_path = self.config['dark_field'].get('data_path')
+# Find tomography bright field images file/folder
+self.tbf_data_path = self.config['bright_field'].get('data_path')
+# Check number of tomography image stacks
+stack_info = self.config['stack_info']
+self.num_tomo_stacks = stack_info.get('num', 1)
+if not msnc.is_int(self.num_tomo_stacks, 1):
+self.num_tomo_stacks = None
+msnc.illegal_value('stack_info:num', self.num_tomo_stacks, 'config file')
+return False
+logging.info(f'num_tomo_stacks = {self.num_tomo_stacks}')
+# Find tomography images file/folders and stack parameters
+stacks = stack_info.get('stacks')
+if stacks is None or len(stacks) is not self.num_tomo_stacks:
+msnc.illegal_value('stack_info:stacks', stacks, 'config file')
+return False
+self.tomo_data_paths = []
+self.tomo_data_indices = []
+self.tomo_ref_heights = []
+for stack in stacks:
+self.tomo_data_paths.append(stack.get('data_path'))
+self.tomo_data_indices.append(stack.get('index'))
+self.tomo_ref_heights.append(stack.get('ref_height'))
+# Check tomo angle (theta) range
+theta_range = self.config.get('theta_range')
+if theta_range is None:
+self.start_theta = 0.
+self.end_theta = 180.
+self.num_thetas = None
+else:
+self.start_theta = theta_range.get('start', 0.)
+if not msnc.is_num(self.start_theta, 0.):
+msnc.illegal_value('theta_range:start', self.start_theta, 'config file')
+is_valid = False
+logging.debug(f'start_theta = {self.start_theta}')
+self.end_theta = theta_range.get('end', 180.)
+if not msnc.is_num(self.end_theta, self.start_theta):
+msnc.illegal_value('theta_range:end', self.end_theta, 'config file')
+is_valid = False
+logging.debug(f'end_theta = {self.end_theta}')
+self.num_thetas = theta_range.get('num')
+if self.num_thetas and not msnc.is_int(self.num_thetas, 1):
+msnc.illegal_value('theta_range:num', self.num_thetas, 'config file')
+self.num_thetas = None
+is_valid = False
+logging.debug(f'num_thetas = {self.num_thetas}')
+return is_valid
+def validate(self):
+"""Returns False if any required config parameter is illegal or missing.
+"""
+is_valid = True
+# Check work_folder (shared by both file formats)
+work_folder = os.path.abspath(self.config.get('work_folder', ''))
+if not os.path.isdir(work_folder):
+msnc.illegal_value('work_folder', work_folder, 'config file')
+is_valid = False
+logging.info(f'work_folder: {work_folder}')
+# Check data filetype (shared by both file formats)
+self.data_filetype = self.config.get('data_filetype', 'tif')
+if not isinstance(self.data_filetype, str) or (self.data_filetype != 'tif' and
+self.data_filetype != 'h5'):
+msnc.illegal_value('data_filetype', self.data_filetype, 'config file')
+if self.suffix == '.yml' or self.suffix == '.yaml':
+is_valid = self._validate_yaml()
+elif self.suffix == '.txt':
+is_valid = self._validate_txt()
+else:
+logging.error(f'Undefined or illegal config file extension: {self.suffix}')
+# Find tomography bright field images file/folder
+if self.tdf_data_path:
+if self.data_filetype == 'h5':
+if isinstance(self.tdf_data_path, str):
+if not os.path.isabs(self.tdf_data_path):
+self.tdf_data_path = os.path.abspath(
+f'{work_folder}/{self.tdf_data_path}')
+else:
+msnc.illegal_value('tdf_data_path', tdf_data_fil, 'config file')
+is_valid = False
+else:
+if isinstance(self.tdf_data_path, int):
+self.tdf_data_path = os.path.abspath(
+f'{work_folder}/{self.tdf_data_path}/nf')
+elif isinstance(self.tdf_data_path, str):
+if not os.path.isabs(self.tdf_data_path):
+self.tdf_data_path = os.path.abspath(
+f'{work_folder}/{self.tdf_data_path}')
+else:
+msnc.illegal_value('tdf_data_path', self.tdf_data_path, 'config file')
+is_valid = False
+logging.info(f'dark field images path = {self.tdf_data_path}')
+# Find tomography bright field images file/folder
+if self.tbf_data_path:
+if self.data_filetype == 'h5':
+if isinstance(self.tbf_data_path, str):
+if not os.path.isabs(self.tbf_data_path):
+self.tbf_data_path = os.path.abspath(
+f'{work_folder}/{self.tbf_data_path}')
+else:
+msnc.illegal_value('tbf_data_path', tbf_data_fil, 'config file')
+is_valid = False
+else:
+if isinstance(self.tbf_data_path, int):
+self.tbf_data_path = os.path.abspath(
+f'{work_folder}/{self.tbf_data_path}/nf')
+elif isinstance(self.tbf_data_path, str):
+if not os.path.isabs(self.tbf_data_path):
+self.tbf_data_path = os.path.abspath(
+f'{work_folder}/{self.tbf_data_path}')
+else:
+msnc.illegal_value('tbf_data_path', self.tbf_data_path, 'config file')
+is_valid = False
+logging.info(f'bright field images path = {self.tbf_data_path}')
+# Find tomography images file/folders and stack parameters
+tomo_data_paths = []
+tomo_data_indices = []
+tomo_ref_heights = []
+for data_path, index, ref_height in zip(self.tomo_data_paths, self.tomo_data_indices,
+self.tomo_ref_heights):
+if self.data_filetype == 'h5':
+if isinstance(data_path, str):
+if not os.path.isabs(data_path):
+data_path = os.path.abspath(f'{work_folder}/{data_path}')
+else:
+msnc.illegal_value(f'stack_info:stacks:data_path', data_path, 'config file')
+is_valid = False
+data_path = None
+else:
+if isinstance(data_path, int):
+data_path = os.path.abspath(f'{work_folder}/{data_path}/nf')
+elif isinstance(data_path, str):
+if not os.path.isabs(data_path):
+data_path = os.path.abspath(f'{work_folder}/{data_path}')
+else:
+msnc.illegal_value(f'stack_info:stacks:data_path', data_path, 'config file')
+is_valid = False
+data_path = None
+tomo_data_paths.append(data_path)
+if index is None:
+if self.num_tomo_stacks > 1:
+logging.error('Missing stack_info:stacks:index in config file')
+is_valid = False
+index = None
+else:
+index = 1
+elif not isinstance(index, int):
+msnc.illegal_value(f'stack_info:stacks:index', index, 'config file')
+is_valid = False
+index = None
+tomo_data_indices.append(index)
+if ref_height is None:
+if self.num_tomo_stacks > 1:
+logging.error('Missing stack_info:stacks:ref_height in config file')
+is_valid = False
+ref_height = None
+else:
+ref_height = 0.
+elif not msnc.is_num(ref_height):
+msnc.illegal_value(f'stack_info:stacks:ref_height', ref_height, 'config file')
+is_valid = False
+ref_height = None
+# Set reference heights relative to first stack
+if (len(tomo_ref_heights) and msnc.is_num(ref_height) and
+msnc.is_num(tomo_ref_heights[0])):
+ref_height = (round(ref_height-tomo_ref_heights[0], 3))
+tomo_ref_heights.append(ref_height)
+tomo_ref_heights[0] = 0.0
+logging.info('tomography data paths:')
+for i in range(self.num_tomo_stacks):
+logging.info(f'    {tomo_data_paths[i]}')
+logging.info(f'tomography data path indices: {tomo_data_indices}')
+logging.info(f'tomography reference heights: {tomo_ref_heights}')
+# Update config in memory
+if self.suffix == '.txt':
+self.config = {}
+dark_field = self.config.get('dark_field')
+if dark_field is None:
+self.config['dark_field'] = {'data_path' : self.tdf_data_path}
+else:
+self.config['dark_field']['data_path'] = self.tdf_data_path
+bright_field = self.config.get('bright_field')
+if bright_field is None:
+self.config['bright_field'] = {'data_path' : self.tbf_data_path}
+else:
+self.config['bright_field']['data_path'] = self.tbf_data_path
+detector = self.config.get('detector')
+if detector is None:
+self.config['detector'] = {'id' : self.detector_id}
+else:
+detector['id'] = self.detector_id
+self.config['work_folder'] = work_folder
+self.config['data_filetype'] = self.data_filetype
+stack_info = self.config.get('stack_info')
+if stack_info is None:
+stacks = []
+for i in range(self.num_tomo_stacks):
+stacks.append({'data_path' : tomo_data_paths[i], 'index' : tomo_data_indices[i],
+'ref_height' : tomo_ref_heights[i]})
+self.config['stack_info'] = {'num' : self.num_tomo_stacks, 'stacks' : stacks}
+else:
+stack_info['num'] = self.num_tomo_stacks
+stacks = stack_info.get('stacks')
+for i,stack in enumerate(stacks):
+stack['data_path'] = tomo_data_paths[i]
+stack['index'] = tomo_data_indices[i]
+stack['ref_height'] = tomo_ref_heights[i]
+if self.num_thetas:
+theta_range = {'start' : self.start_theta, 'end' : self.end_theta,
+'num' : self.num_thetas}
+else:
+theta_range = {'start' : self.start_theta, 'end' : self.end_theta}
+self.config['theta_range'] = theta_range
+# Cleanup temporary validation variables
+del self.tdf_data_path
+del self.tbf_data_path
+del self.detector_id
+del self.data_filetype
+del self.num_tomo_stacks
+del self.tomo_data_paths
+del self.tomo_data_indices
+del self.tomo_ref_heights
+del self.start_theta
+del self.end_theta
+del self.num_thetas
+return is_valid
+class Tomo:
+"""Processing tomography data with misalignment.
+"""
+def __init__(self, config_file=None, config_dict=None, config_out=None, output_folder='.',
+log_level='INFO', log_stream='tomo.log', galaxy_flag=False, test_mode=False):
+"""Initialize with optional config input file or dictionary
+"""
+self.ncore = mp.cpu_count()
+self.config_out = config_out
+self.output_folder = output_folder
+self.galaxy_flag = galaxy_flag
+self.test_mode = test_mode
+self.save_plots = True # Make input argument?
+self.save_plots_only = True # Make input argument?
+self.cf = None
+self.config = None
+self.is_valid = True
+self.tdf = np.array([])
+self.tbf = np.array([])
+self.tomo_stacks = []
+self.tomo_recon_stacks = []
+# Validate input parameters
+if config_file is not None and not os.path.isfile(config_file):
+raise OSError(f'Invalid config_file input {config_file} {type(config_file)}')
+if config_dict is not None and not isinstance(config_dict, dict):
+raise ValueError(f'Invalid config_dict input {config_dict} {type(config_dict)}')
+if config_out is not None:
+if isinstance(config_out, str):
+if isinstance(log_stream, str):
+path = os.path.split(log_stream)[0]
+if path and not os.path.isdir(path):
+raise OSError(f'Invalid log_stream path')
+else:
+raise OSError(f'Invalid config_out input {config_out} {type(config_out)}')
+if not os.path.isdir(output_folder):
+raise OSError(f'Invalid output_folder input {output_folder} {type(output_folder)}')
+if isinstance(log_stream, str):
+path = os.path.split(log_stream)[0]
+if path and not os.path.isdir(path):
+raise OSError(f'Invalid log_stream path')
+if not os.path.isabs(path):
+log_stream = f'{output_folder}/{log_stream}'
+if not isinstance(galaxy_flag, bool):
+raise ValueError(f'Invalid galaxy_flag input {galaxy_flag} {type(galaxy_flag)}')
+if not isinstance(test_mode, bool):
+raise ValueError(f'Invalid test_mode input {test_mode} {type(test_mode)}')
+# Set log configuration
+logging_format = '%(asctime)s : %(levelname)s - %(module)s : %(funcName)s - %(message)s'
+if self.test_mode:
+self.save_plots_only = True
+if isinstance(log_stream, str):
+logging.basicConfig(filename=f'{log_stream}', filemode='w',
+format=logging_format, level=logging.WARNING, force=True)
+elif isinstance(log_stream, io.TextIOWrapper):
+logging.basicConfig(filemode='w', format=logging_format, level=logging.WARNING,
+stream=log_stream, force=True)
+else:
+raise ValueError(f'Invalid log_stream: {log_stream}')
+logging.warning('Ignoring log_level argument in test mode')
+else:
+level = getattr(logging, log_level.upper(), None)
+if not isinstance(level, int):
+raise ValueError(f'Invalid log_level: {log_level}')
+if log_stream is sys.stdout:
+logging.basicConfig(format=logging_format, level=level, force=True,
+handlers=[logging.StreamHandler()])
+else:
+if isinstance(log_stream, str):
+logging.basicConfig(filename=f'{log_stream}', filemode='w',
+format=logging_format, level=level, force=True)
+elif isinstance(log_stream, io.TextIOWrapper):
+logging.basicConfig(filemode='w', format=logging_format, level=level,
+stream=log_stream, force=True)
+else:
+raise ValueError(f'Invalid log_stream: {log_stream}')
+stream_handler = logging.StreamHandler()
+logging.getLogger().addHandler(stream_handler)
+stream_handler.setLevel(logging.WARNING)
+stream_handler.setFormatter(logging.Formatter(logging_format))
+# Check/set output config file name
+if self.config_out is None:
+self.config_out = f'{self.output_folder}/config.yaml'
+elif (self.config_out is os.path.basename(self.config_out) and
+not os.path.isabs(self.config_out)):
+self.config_out = f'{self.output_folder}/{self.config_out}'
+logging.info(f'ncore = {self.ncore}')
+logging.debug(f'config_file = {config_file}')
+logging.debug(f'config_dict = {config_dict}')
+logging.debug(f'config_out = {self.config_out}')
+logging.debug(f'output_folder = {self.output_folder}')
+logging.debug(f'log_stream = {log_stream}')
+logging.debug(f'log_level = {log_level}')
+logging.debug(f'galaxy_flag = {self.galaxy_flag}')
+logging.debug(f'test_mode = {self.test_mode}')
+# Create config object and load config file
+self.cf = ConfigTomo(config_file, config_dict)
+if not self.cf.load_flag:
+self.is_valid = False
+return
+if self.galaxy_flag:
+self.ncore = 1 #RV can I set this? mp.cpu_count()
+assert(self.output_folder == '.')
+assert(self.test_mode is False)
+self.save_plots = True
+self.save_plots_only = True
+else:
+# Input validation is already performed during link_data_to_galaxy
+# Check config file parameters
+self.is_valid =  self.cf.validate()
+# Load detector info file
+df = msnc.Detector(self.cf.config['detector']['id'])
+# Check detector info file parameters
+if df.validate():
+pixel_size = df.getPixelSize()
+num_rows, num_columns = df.getDimensions()
+if not pixel_size or not num_rows or not num_columns:
+self.is_valid = False
+else:
+pixel_size = None
+num_rows = None
+num_columns = None
+self.is_valid = False
+# Update config
+self.cf.config['detector']['pixel_size'] = pixel_size
+self.cf.config['detector']['rows'] = num_rows
+self.cf.config['detector']['columns'] = num_columns
+logging.debug(f'pixel_size = self.cf.config["detector"]["pixel_size"]')
+logging.debug(f'num_rows: {self.cf.config["detector"]["rows"]}')
+logging.debug(f'num_columns: {self.cf.config["detector"]["columns"]}')
+# Safe config to file
+if self.is_valid:
+self.cf.saveFile(self.config_out)
+# Initialize shortcut to config
+self.config = self.cf.config
+# Initialize tomography stack
+num_tomo_stacks = self.config['stack_info']['num']
+if num_tomo_stacks:
+self.tomo_stacks = [np.array([]) for _ in range(num_tomo_stacks)]
+self.tomo_recon_stacks = [np.array([]) for _ in range(num_tomo_stacks)]
+logging.debug(f'save_plots = {self.save_plots}')
+logging.debug(f'save_plots_only = {self.save_plots_only}')
+def _selectImageRanges(self, available_stacks=None):
+"""Select image files to be included in analysis.
+"""
+self.is_valid = True
+stack_info = self.config['stack_info']
+if available_stacks is None:
+available_stacks = [False]*stack_info['num']
+elif len(available_stacks) != stack_info['num']:
+logging.warning('Illegal dimension of available_stacks in getImageFiles '+
+f'({len(available_stacks)}');
+available_stacks = [False]*stack_info['num']
+# Check number of tomography angles/thetas
+num_thetas = self.config['theta_range'].get('num')
+if num_thetas is None:
+num_thetas = pyip.inputInt('\nEnter the number of thetas (>0): ', greaterThan=0)
+elif not msnc.is_int(num_thetas, 0):
+msnc.illegal_value('num_thetas', num_thetas, 'config file')
+self.is_valid = False
+return
+self.config['theta_range']['num'] = num_thetas
+logging.debug(f'num_thetas = {self.config["theta_range"]["num"]}')
+# Find tomography dark field images
+dark_field = self.config['dark_field']
+img_start = dark_field.get('img_start', -1)
+img_offset = dark_field.get('img_offset', -1)
+num_imgs = dark_field.get('num', 0)
+if not self.test_mode:
+img_start, img_offset, num_imgs = msnc.selectImageRange(img_start, img_offset,
+num_imgs, 'dark field')
+if img_start < 0 or num_imgs < 1:
+logging.error('Unable to find suitable dark field images')
+if dark_field['data_path']:
+self.is_valid = False
+dark_field['img_start'] = img_start
+dark_field['img_offset'] = img_offset
+dark_field['num'] = num_imgs
+logging.debug(f'Dark field image start index: {dark_field["img_start"]}')
+logging.debug(f'Dark field image offset: {dark_field["img_offset"]}')
+logging.debug(f'Number of dark field images: {dark_field["num"]}')
+# Find tomography bright field images
+bright_field = self.config['bright_field']
+img_start = bright_field.get('img_start', -1)
+img_offset = bright_field.get('img_offset', -1)
+num_imgs = bright_field.get('num', 0)
+if not self.test_mode:
+img_start, img_offset, num_imgs = msnc.selectImageRange(img_start, img_offset,
+num_imgs, 'bright field')
+if img_start < 0 or num_imgs < 1:
+logging.error('Unable to find suitable bright field images')
+if bright_field['data_path']:
+self.is_valid = False
+bright_field['img_start'] = img_start
+bright_field['img_offset'] = img_offset
+bright_field['num'] = num_imgs
+logging.debug(f'Bright field image start index: {bright_field["img_start"]}')
+logging.debug(f'Bright field image offset: {bright_field["img_offset"]}')
+logging.debug(f'Number of bright field images: {bright_field["num"]}')
+# Find tomography images
+for i,stack in enumerate(stack_info['stacks']):
+# Check if stack is already loaded or available
+if self.tomo_stacks[i].size or available_stacks[i]:
+continue
+index = stack['index']
+img_start = stack.get('img_start', -1)
+img_offset = stack.get('img_offset', -1)
+num_imgs = stack.get('num', 0)
+if not self.test_mode:
+img_start, img_offset, num_imgs = msnc.selectImageRange(img_start, img_offset,
+num_imgs, f'tomography stack {index}', num_thetas)
+if img_start < 0 or num_imgs != num_thetas:
+logging.error('Unable to find suitable tomography images')
+self.is_valid = False
+stack['img_start'] = img_start
+stack['img_offset'] = img_offset
+stack['num'] = num_imgs
+logging.debug(f'Tomography stack {index} image start index: {stack["img_start"]}')
+logging.debug(f'Tomography stack {index} image offset: {stack["img_offset"]}')
+logging.debug(f'Number of tomography images for stack {index}: {stack["num"]}')
+available_stacks[i] = True
+# Safe updated config to file
+if self.is_valid:
+self.cf.saveFile(self.config_out)
+return
+def _genDark(self, tdf_files, dark_field_pngname):
+"""Generate dark field.
+"""
+# Load the dark field images
+logging.debug('Loading dark field...')
+dark_field = self.config['dark_field']
+tdf_stack = msnc.loadImageStack(tdf_files, self.config['data_filetype'],
+dark_field['img_offset'], dark_field['num'])
+# Take median
+self.tdf = np.median(tdf_stack, axis=0)
+del tdf_stack
+# RV make input of some kind (not always needed)
+tdf_cutoff = 21
+self.tdf[self.tdf > tdf_cutoff] = np.nan
+tdf_mean = np.nanmean(self.tdf)
+logging.debug(f'tdf_cutoff = {tdf_cutoff}')
+logging.debug(f'tdf_mean = {tdf_mean}')
+np.nan_to_num(self.tdf, copy=False, nan=tdf_mean, posinf=tdf_mean, neginf=0.)
+if not self.test_mode and not self.galaxy_flag:
+msnc.quickImshow(self.tdf, title='dark field', path=self.output_folder,
+save_fig=self.save_plots, save_only=self.save_plots_only)
+elif self.galaxy_flag:
+msnc.quickImshow(self.tdf, title='dark field', name=dark_field_pngname,
+save_fig=True, save_only=True)
+def _genBright(self, tbf_files, bright_field_pngname):
+"""Generate bright field.
+"""
+# Load the bright field images
+logging.debug('Loading bright field...')
+bright_field = self.config['bright_field']
+tbf_stack = msnc.loadImageStack(tbf_files, self.config['data_filetype'],
+bright_field['img_offset'], bright_field['num'])
+# Take median
+"""Median or mean: It may be best to try the median because of some image
+artifacts that arise due to crinkles in the upstream kapton tape windows
+causing some phase contrast images to appear on the detector.
+One thing that also may be useful in a future implementation is to do a
+brightfield adjustment on EACH frame of the tomo based on a ROI in the
+corner of the frame where there is no sample but there is the direct X-ray
+beam because there is frame to frame fluctuations from the incoming beam.
+We don’t typically account for them but potentially could.
+"""
+self.tbf = np.median(tbf_stack, axis=0)
+del tbf_stack
+# Subtract dark field
+if self.tdf.size:
+self.tbf -= self.tdf
+else:
+logging.warning('Dark field unavailable')
+if not self.test_mode and not self.galaxy_flag:
+msnc.quickImshow(self.tbf, title='bright field', path=self.output_folder,
+save_fig=self.save_plots, save_only=self.save_plots_only)
+elif self.galaxy_flag:
+msnc.quickImshow(self.tbf, title='bright field', name=bright_field_pngname,
+save_fig=True, save_only=True)
+def _setDetectorBounds(self, tomo_stack_files, tomo_field_pngname, detectorbounds_pngname):
+"""Set vertical detector bounds for image stack.
+"""
+preprocess = self.config.get('preprocess')
+if preprocess is None:
+img_x_bounds = [None, None]
+else:
+img_x_bounds = preprocess.get('img_x_bounds', [0, self.tbf.shape[0]])
+if self.test_mode:
+# Update config and save to file
+if preprocess is None:
+self.cf.config['preprocess'] = {'img_x_bounds' : [0, self.tbf.shape[0]]}
+else:
+preprocess['img_x_bounds'] = img_x_bounds
+self.cf.saveFile(self.config_out)
+return
+# Check reference heights
+pixel_size = self.config['detector']['pixel_size']
+if pixel_size is None:
+raise ValueError('Detector pixel size unavailable')
+if not self.tbf.size:
+raise ValueError('Bright field unavailable')
+num_x_min = None
+num_tomo_stacks = self.config['stack_info']['num']
+stacks = self.config['stack_info']['stacks']
+if num_tomo_stacks > 1:
+delta_z = stacks[1]['ref_height']-stacks[0]['ref_height']
+for i in range(2, num_tomo_stacks):
+delta_z = min(delta_z, stacks[i]['ref_height']-stacks[i-1]['ref_height'])
+logging.debug(f'delta_z = {delta_z}')
+num_x_min = int(delta_z/pixel_size)+1
+logging.debug(f'num_x_min = {num_x_min}')
+if num_x_min > self.tbf.shape[0]:
+logging.warning('Image bounds and pixel size prevent seamless stacking')
+num_x_min = self.tbf.shape[0]
+# Select image bounds
+if self.galaxy_flag:
+if num_x_min is None or num_x_min >= self.tbf.shape[0]:
+img_x_bounds = [0, self.tbf.shape[0]]
+else:
+margin = int(num_x_min/10)
+offset = min(0, int((self.tbf.shape[0]-num_x_min)/2-margin))
+img_x_bounds = [offset, num_x_min+offset+2*margin]
+tomo_stack = msnc.loadImageStack(tomo_stack_files[0], self.config['data_filetype'],
+stacks[0]['img_offset'], 1)
+x_sum = np.sum(tomo_stack[0,:,:], 1)
+title = f'tomography image at theta={self.config["theta_range"]["start"]}'
+msnc.quickImshow(tomo_stack[0,:,:], title=title, name=tomo_field_pngname,
+save_fig=True, save_only=True)
+msnc.quickPlot((range(x_sum.size), x_sum),
+([img_x_bounds[0], img_x_bounds[0]], [x_sum.min(), x_sum.max()], 'r-'),
+([img_x_bounds[1]-1, img_x_bounds[1]-1], [x_sum.min(), x_sum.max()], 'r-'),
+title='sum over theta and y', name=detectorbounds_pngname,
+save_fig=True, save_only=True)
+# Update config and save to file
+if preprocess is None:
+self.cf.config['preprocess'] = {'img_x_bounds' : img_x_bounds}
+else:
+preprocess['img_x_bounds'] = img_x_bounds
+self.cf.saveFile(self.config_out)
+return
+# For one tomography stack only: load the first image
+title = None
+msnc.quickImshow(self.tbf, title='bright field')
+if num_tomo_stacks == 1:
+tomo_stack = msnc.loadImageStack(tomo_stack_files[0], self.config['data_filetype'],
+stacks[0]['img_offset'], 1)
+title = f'tomography image at theta={self.config["theta_range"]["start"]}'
+msnc.quickImshow(tomo_stack[0,:,:], title=title)
+tomo_or_bright = pyip.inputNum('\nSelect image bounds from bright field (0) or '+
+'from first tomography image (1): ', min=0, max=1)
+else:
+print('\nSelect image bounds from bright field')
+tomo_or_bright = 0
+if tomo_or_bright:
+x_sum = np.sum(tomo_stack[0,:,:], 1)
+use_bounds = 'no'
+if img_x_bounds[0] is not None and img_x_bounds[1] is not None:
+if img_x_bounds[0] < 0:
+msnc.illegal_value('img_x_bounds[0]', img_x_bounds[0], 'config file')
+img_x_bounds[0] = 0
+if not img_x_bounds[0] < img_x_bounds[1] <= x_sum.size:
+msnc.illegal_value('img_x_bounds[1]', img_x_bounds[1], 'config file')
+img_x_bounds[1] = x_sum.size
+tomo_tmp = tomo_stack[0,:,:]
+tomo_tmp[img_x_bounds[0],:] = tomo_stack[0,:,:].max()
+tomo_tmp[img_x_bounds[1]-1,:] = tomo_stack[0,:,:].max()
+title = f'tomography image at theta={self.config["theta_range"]["start"]}'
+msnc.quickImshow(tomo_stack[0,:,:], title=title)
+msnc.quickPlot((range(x_sum.size), x_sum),
+([img_x_bounds[0], img_x_bounds[0]], [x_sum.min(), x_sum.max()], 'r-'),
+([img_x_bounds[1]-1, img_x_bounds[1]-1], [x_sum.min(), x_sum.max()], 'r-'),
+title='sum over theta and y')
+print(f'lower bound = {img_x_bounds[0]} (inclusive)\n'+
+f'upper bound = {img_x_bounds[1]} (exclusive)]')
+use_bounds =  pyip.inputYesNo('Accept these bounds ([y]/n)?: ', blank=True)
+if use_bounds == 'no':
+img_x_bounds = msnc.selectImageBounds(tomo_stack[0,:,:], 0,
+img_x_bounds[0], img_x_bounds[1], num_x_min,
+f'tomography image at theta={self.config["theta_range"]["start"]}')
+if num_x_min is not None and img_x_bounds[1]-img_x_bounds[0]+1 < num_x_min:
+logging.warning('Image bounds and pixel size prevent seamless stacking')
+tomo_tmp = tomo_stack[0,:,:]
+tomo_tmp[img_x_bounds[0],:] = tomo_stack[0,:,:].max()
+tomo_tmp[img_x_bounds[1]-1,:] = tomo_stack[0,:,:].max()
+title = f'tomography image at theta={self.config["theta_range"]["start"]}'
+msnc.quickImshow(tomo_stack[0,:,:], title=title, path=self.output_folder,
+save_fig=self.save_plots, save_only=True)
+msnc.quickPlot(range(img_x_bounds[0], img_x_bounds[1]),
+x_sum[img_x_bounds[0]:img_x_bounds[1]],
+title='sum over theta and y', path=self.output_folder,
+save_fig=self.save_plots, save_only=True)
+else:
+x_sum = np.sum(self.tbf, 1)
+use_bounds = 'no'
+if img_x_bounds[0] is not None and img_x_bounds[1] is not None:
+if img_x_bounds[0] < 0:
+msnc.illegal_value('img_x_bounds[0]', img_x_bounds[0], 'config file')
+img_x_bounds[0] = 0
+if not img_x_bounds[0] < img_x_bounds[1] <= x_sum.size:
+msnc.illegal_value('img_x_bounds[1]', img_x_bounds[1], 'config file')
+img_x_bounds[1] = x_sum.size
+msnc.quickPlot((range(x_sum.size), x_sum),
+([img_x_bounds[0], img_x_bounds[0]], [x_sum.min(), x_sum.max()], 'r-'),
+([img_x_bounds[1]-1, img_x_bounds[1]-1], [x_sum.min(), x_sum.max()], 'r-'),
+title='sum over theta and y')
+print(f'lower bound = {img_x_bounds[0]} (inclusive)\n'+
+f'upper bound = {img_x_bounds[1]} (exclusive)]')
+use_bounds =  pyip.inputYesNo('Accept these bounds ([y]/n)?: ', blank=True)
+if use_bounds == 'no':
+fit = msnc.fitStep(y=x_sum, model='rectangle', form='atan')
+x_low = fit.get('center1', None)
+x_upp = fit.get('center2', None)
+if (x_low is not None and x_low >= 0 and x_upp is not None and
+x_low < x_upp < x_sum.size):
+x_low = int(x_low-(x_upp-x_low)/10)
+if x_low < 0:
+x_low = 0
+x_upp = int(x_upp+(x_upp-x_low)/10)
+if x_upp >= x_sum.size:
+x_upp = x_sum.size
+msnc.quickPlot((range(x_sum.size), x_sum),
+([x_low, x_low], [x_sum.min(), x_sum.max()], 'r-'),
+([x_upp, x_upp], [x_sum.min(), x_sum.max()], 'r-'),
+title='sum over theta and y')
+print(f'lower bound = {x_low} (inclusive)\nupper bound = {x_upp} (exclusive)]')
+use_fit =  pyip.inputYesNo('Accept these bounds ([y]/n)?: ', blank=True)
+if use_fit == 'no':
+img_x_bounds = msnc.selectArrayBounds(x_sum, img_x_bounds[0], img_x_bounds[1],
+num_x_min, 'sum over theta and y')
+else:
+img_x_bounds = [x_low, x_upp]
+if num_x_min is not None and img_x_bounds[1]-img_x_bounds[0]+1 < num_x_min:
+logging.warning('Image bounds and pixel size prevent seamless stacking')
+msnc.quickPlot(range(img_x_bounds[0], img_x_bounds[1]),
+x_sum[img_x_bounds[0]:img_x_bounds[1]],
+title='sum over theta and y', path=self.output_folder,
+save_fig=self.save_plots, save_only=True)
+logging.debug(f'img_x_bounds: {img_x_bounds}')
+if self.save_plots_only:
+msnc.clearFig('bright field')
+msnc.clearFig('sum over theta and y')
+if title:
+msnc.clearFig(title)
+# Update config and save to file
+if preprocess is None:
+self.cf.config['preprocess'] = {'img_x_bounds' : img_x_bounds}
+else:
+preprocess['img_x_bounds'] = img_x_bounds
+self.cf.saveFile(self.config_out)
+def _setZoomOrSkip(self):
+"""Set zoom and/or theta skip to reduce memory the requirement for the analysis.
+"""
+preprocess = self.config.get('preprocess')
+zoom_perc = 100
+if not self.galaxy_flag:
+if preprocess is None or 'zoom_perc' not in preprocess:
+if pyip.inputYesNo(
+'\nDo you want to zoom in to reduce memory requirement (y/[n])? ',
+blank=True) == 'yes':
+zoom_perc = pyip.inputInt('    Enter zoom percentage [1, 100]: ',
+min=1, max=100)
+else:
+zoom_perc = preprocess['zoom_perc']
+if msnc.is_num(zoom_perc, 1., 100.):
+zoom_perc = int(zoom_perc)
+else:
+msnc.illegal_value('zoom_perc', zoom_perc, 'config file')
+zoom_perc = 100
+num_theta_skip = 0
+if not self.galaxy_flag:
+if preprocess is None or 'num_theta_skip' not in preprocess:
+if pyip.inputYesNo(
+'Do you want to skip thetas to reduce memory requirement (y/[n])? ',
+blank=True) == 'yes':
+num_theta_skip = pyip.inputInt('    Enter the number skip theta interval'+
+f' [0, {self.num_thetas-1}]: ', min=0, max=self.num_thetas-1)
+else:
+num_theta_skip = preprocess['num_theta_skip']
+if not msnc.is_int(num_theta_skip, 0):
+msnc.illegal_value('num_theta_skip', num_theta_skip, 'config file')
+num_theta_skip = 0
+logging.info(f'zoom_perc = {zoom_perc}')
+logging.info(f'num_theta_skip = {num_theta_skip}')
+# Update config and save to file
+if preprocess is None:
+self.cf.config['preprocess'] = {'zoom_perc' : zoom_perc,
+'num_theta_skip' : num_theta_skip}
+else:
+preprocess['zoom_perc'] = zoom_perc
+preprocess['num_theta_skip'] = num_theta_skip
+self.cf.saveFile(self.config_out)
+def _loadTomo(self, base_name, index, required=False):
+"""Load a tomography stack.
+"""
+# stack order: row,theta,column
+zoom_perc = None
+preprocess = self.config.get('preprocess')
+if preprocess:
+zoom_perc = preprocess.get('zoom_perc')
+if zoom_perc is None or zoom_perc == 100:
+title = f'{base_name} fullres'
+else:
+title = f'{base_name} {zoom_perc}p'
+title += f'_{index}'
+tomo_file = re.sub(r"\s+", '_', f'{self.output_folder}/{title}.npy')
+load_flag = 'no'
+available = False
+if os.path.isfile(tomo_file):
+available = True
+if required:
+load_flag = 'yes'
+else:
+load_flag = pyip.inputYesNo(f'\nDo you want to load {tomo_file} (y/n)? ')
+stack = np.array([])
+if load_flag == 'yes':
+t0 = time()
+logging.info(f'Loading {tomo_file} ...')
+try:
+stack = np.load(tomo_file)
+except IOError or ValueError:
+stack = np.array([])
+logging.error(f'Error loading {tomo_file}')
+logging.info(f'... done in {time()-t0:.2f} seconds!')
+if stack.size:
+msnc.quickImshow(stack[:,0,:], title=title, path=self.output_folder,
+save_fig=self.save_plots, save_only=self.save_plots_only)
+return stack, available
+def _saveTomo(self, base_name, stack, index=None):
+"""Save a tomography stack.
+"""
+zoom_perc = None
+preprocess = self.config.get('preprocess')
+if preprocess:
+zoom_perc = preprocess.get('zoom_perc')
+if zoom_perc is None or zoom_perc == 100:
+title = f'{base_name} fullres'
+else:
+title = f'{base_name} {zoom_perc}p'
+if index:
+title += f'_{index}'
+tomo_file = re.sub(r"\s+", '_', f'{self.output_folder}/{title}.npy')
+t0 = time()
+logging.info(f'Saving {tomo_file} ...')
+np.save(tomo_file, stack)
+logging.info(f'... done in {time()-t0:.2f} seconds!')
+def _genTomo(self, tomo_stack_files, available_stacks):
+"""Generate tomography fields.
+"""
+stacks = self.config['stack_info']['stacks']
+assert(len(self.tomo_stacks) == self.config['stack_info']['num'])
+assert(len(self.tomo_stacks) == len(stacks))
+if len(available_stacks) != len(stacks):
+logging.warning('Illegal dimension of available_stacks in _genTomo'+
+f'({len(available_stacks)}');
+available_stacks = [False]*self.num_tomo_stacks
+preprocess = self.config.get('preprocess')
+if preprocess is None:
+img_x_bounds = [0, self.tbf.shape[0]]
+img_y_bounds = [0, self.tbf.shape[1]]
+zoom_perc = preprocess['zoom_perc']
+num_theta_skip = preprocess['num_theta_skip']
+else:
+img_x_bounds = preprocess.get('img_x_bounds', [0, self.tbf.shape[0]])
+img_y_bounds = preprocess.get('img_y_bounds', [0, self.tbf.shape[1]])
+zoom_perc = 100
+num_theta_skip = 0
+if self.tdf.size:
+tdf = self.tdf[img_x_bounds[0]:img_x_bounds[1],img_y_bounds[0]:img_y_bounds[1]]
+else:
+logging.warning('Dark field unavailable')
+if not self.tbf.size:
+raise ValueError('Bright field unavailable')
+tbf = self.tbf[img_x_bounds[0]:img_x_bounds[1],img_y_bounds[0]:img_y_bounds[1]]
+for i,stack in enumerate(stacks):
+# Check if stack is already loaded or available
+if self.tomo_stacks[i].size or available_stacks[i]:
+continue
+# Load a stack of tomography images
+t0 = time()
+tomo_stack = msnc.loadImageStack(tomo_stack_files[i], self.config['data_filetype'],
+stack['img_offset'], self.config['theta_range']['num'], num_theta_skip,
+img_x_bounds, img_y_bounds)
+tomo_stack = tomo_stack.astype('float64')
+logging.debug(f'loading took {time()-t0:.2f} seconds!')
+# Subtract dark field
+if self.tdf.size:
+t0 = time()
+with set_numexpr_threads(self.ncore):
+ne.evaluate('tomo_stack-tdf', out=tomo_stack)
+logging.debug(f'subtracting dark field took {time()-t0:.2f} seconds!')
+# Normalize
+t0 = time()
+with set_numexpr_threads(self.ncore):
+ne.evaluate('tomo_stack/tbf', out=tomo_stack, truediv=True)
+logging.debug(f'normalizing took {time()-t0:.2f} seconds!')
+# Remove non-positive values and linearize data
+t0 = time()
+cutoff = 1.e-6
+with set_numexpr_threads(self.ncore):
+ne.evaluate('where(tomo_stack<cutoff, cutoff, tomo_stack)', out=tomo_stack)
+with set_numexpr_threads(self.ncore):
+ne.evaluate('-log(tomo_stack)', out=tomo_stack)
+logging.debug('removing non-positive values and linearizing data took '+
+f'{time()-t0:.2f} seconds!')
+# Get rid of nans/infs that may be introduced by normalization
+t0 = time()
+np.where(np.isfinite(tomo_stack), tomo_stack, 0.)
+logging.debug(f'remove nans/infs took {time()-t0:.2f} seconds!')
+# Downsize tomography stack to smaller size
+tomo_stack = tomo_stack.astype('float32')
+if not self.galaxy_flag:
+index = stack['index']
+title = f'red stack fullres {index}'
+if not self.test_mode:
+msnc.quickImshow(tomo_stack[0,:,:], title=title, path=self.output_folder,
+save_fig=self.save_plots, save_only=self.save_plots_only)
+if zoom_perc != 100:
+t0 = time()
+logging.info(f'Zooming in ...')
+tomo_zoom_list = []
+for j in range(tomo_stack.shape[0]):
+tomo_zoom = spi.zoom(tomo_stack[j,:,:], 0.01*zoom_perc)
+tomo_zoom_list.append(tomo_zoom)
+tomo_stack = np.stack([tomo_zoom for tomo_zoom in tomo_zoom_list])
+logging.info(f'... done in {time()-t0:.2f} seconds!')
+del tomo_zoom_list
+if not self.galaxy_flag:
+title = f'red stack {zoom_perc}p {index}'
+if not self.test_mode:
+msnc.quickImshow(tomo_stack[0,:,:], title=title, path=self.output_folder,
+save_fig=self.save_plots, save_only=self.save_plots_only)
+# Convert tomography stack from theta,row,column to row,theta,column
+tomo_stack = np.swapaxes(tomo_stack, 0, 1)
+# Save tomography stack to file
+if not self.galaxy_flag:
+if not self.test_mode:
+self._saveTomo('red stack', tomo_stack, index)
+else:
+np.savetxt(self.output_folder+f'red_stack_{index}.txt',
+tomo_stack[0,:,:], fmt='%.6e')
+# Combine stacks
+t0 = time()
+self.tomo_stacks[i] = tomo_stack
+logging.debug(f'combining nstack took {time()-t0:.2f} seconds!')
+# Update config and save to file
+stack['preprocessed'] = True
+self.cf.saveFile(self.config_out)
+if self.tdf.size:
+del tdf
+del tbf
+def _reconstructOnePlane(self, tomo_plane_T, center, thetas_deg, eff_pixel_size,
+cross_sectional_dim, plot_sinogram=True):
+"""Invert the sinogram for a single tomography plane.
+"""
+# tomo_plane_T index order: column,theta
+assert(0 <= center < tomo_plane_T.shape[0])
+center_offset = center-tomo_plane_T.shape[0]/2
+two_offset = 2*int(np.round(center_offset))
+two_offset_abs = np.abs(two_offset)
+max_rad = int(0.5*(cross_sectional_dim/eff_pixel_size)*1.1) # 10% slack to avoid edge effects
+dist_from_edge = max(1, int(np.floor((tomo_plane_T.shape[0]-two_offset_abs)/2.)-max_rad))
+if two_offset >= 0:
+logging.debug(f'sinogram range = [{two_offset+dist_from_edge}, {-dist_from_edge}]')
+sinogram = tomo_plane_T[two_offset+dist_from_edge:-dist_from_edge,:]
+else:
+logging.debug(f'sinogram range = [{dist_from_edge}, {two_offset-dist_from_edge}]')
+sinogram = tomo_plane_T[dist_from_edge:two_offset-dist_from_edge,:]
+if plot_sinogram:
+msnc.quickImshow(sinogram.T, f'sinogram center offset{center_offset:.2f}',
+path=self.output_folder, save_fig=self.save_plots,
+save_only=self.save_plots_only, aspect='auto')
+# Inverting sinogram
+t0 = time()
+recon_sinogram = iradon(sinogram, theta=thetas_deg, circle=True)
+logging.debug(f'inverting sinogram took {time()-t0:.2f} seconds!')
+del sinogram
+# Removing ring artifacts
+# RV parameters for the denoise, gaussian, and ring removal will be different for different feature sizes
+t0 = time()
+#        recon_sinogram = filters.gaussian(recon_sinogram, 3.0)
+recon_sinogram = spi.gaussian_filter(recon_sinogram, 0.5)
+recon_clean = np.expand_dims(recon_sinogram, axis=0)
+del recon_sinogram
+recon_clean = tomopy.misc.corr.remove_ring(recon_clean, rwidth=17)
+logging.debug(f'filtering and removing ring artifact took {time()-t0:.2f} seconds!')
+return recon_clean
+def _plotEdgesOnePlane(self, recon_plane, base_name, weight=0.001):
+# RV parameters for the denoise, gaussian, and ring removal will be different for different feature sizes
+edges = denoise_tv_chambolle(recon_plane, weight = weight)
+vmax = np.max(edges[0,:,:])
+vmin = -vmax
+msnc.quickImshow(edges[0,:,:], f'{base_name} coolwarm', path=self.output_folder,
+save_fig=self.save_plots, save_only=self.save_plots_only, cmap='coolwarm')
+msnc.quickImshow(edges[0,:,:], f'{base_name} gray', path=self.output_folder,
+save_fig=self.save_plots, save_only=self.save_plots_only, cmap='gray',
+vmin=vmin, vmax=vmax)
+del edges
+def _findCenterOnePlane(self, sinogram, row, thetas_deg, eff_pixel_size, cross_sectional_dim,
+tol=0.1):
+"""Find center for a single tomography plane.
+"""
+# sinogram index order: theta,column
+# need index order column,theta for iradon, so take transpose
+sinogram_T = sinogram.T
+center = sinogram.shape[1]/2
+# try automatic center finding routines for initial value
+tomo_center = tomopy.find_center_vo(sinogram)
+center_offset_vo = tomo_center-center
+print(f'Center at row {row} using Nghia Vo’s method = {center_offset_vo:.2f}')
+recon_plane = self._reconstructOnePlane(sinogram_T, tomo_center, thetas_deg,
+eff_pixel_size, cross_sectional_dim, False)
+base_name=f'edges row{row} center_offset_vo{center_offset_vo:.2f}'
+self._plotEdgesOnePlane(recon_plane, base_name)
+if pyip.inputYesNo('Try finding center using phase correlation (y/[n])? ',
+blank=True) == 'yes':
+tomo_center = tomopy.find_center_pc(sinogram, sinogram, tol=0.1,
+rotc_guess=tomo_center)
+error = 1.
+while error > tol:
+prev = tomo_center
+tomo_center = tomopy.find_center_pc(sinogram, sinogram, tol=tol,
+rotc_guess=tomo_center)
+error = np.abs(tomo_center-prev)
+center_offset = tomo_center-center
+print(f'Center at row {row} using phase correlation = {center_offset:.2f}')
+recon_plane = self._reconstructOnePlane(sinogram_T, tomo_center, thetas_deg,
+eff_pixel_size, cross_sectional_dim, False)
+base_name=f'edges row{row} center_offset{center_offset:.2f}'
+self._plotEdgesOnePlane(recon_plane, base_name)
+if pyip.inputYesNo('Accept a center location ([y]) or continue search (n)? ',
+blank=True) != 'no':
+del sinogram_T
+del recon_plane
+center_offset = pyip.inputNum(
+f'    Enter chosen center offset [{-int(center)}, {int(center)}] '+
+f'([{center_offset_vo}])): ', blank=True)
+if center_offset == '':
+center_offset = center_offset_vo
+return float(center_offset)
+while True:
+center_offset_low = pyip.inputInt('\nEnter lower bound for center offset '+
+f'[{-int(center)}, {int(center)}]: ', min=-int(center), max=int(center))
+center_offset_upp = pyip.inputInt('Enter upper bound for center offset '+
+f'[{center_offset_low}, {int(center)}]: ',
+min=center_offset_low, max=int(center))
+if center_offset_upp == center_offset_low:
+center_offset_step = 1
+else:
+center_offset_step = pyip.inputInt('Enter step size for center offset search '+
+f'[1, {center_offset_upp-center_offset_low}]: ',
+min=1, max=center_offset_upp-center_offset_low)
+for center_offset in range(center_offset_low, center_offset_upp+center_offset_step,
+center_offset_step):
+logging.info(f'center_offset = {center_offset}')
+recon_plane = self._reconstructOnePlane(sinogram_T, center_offset+center,
+thetas_deg, eff_pixel_size, cross_sectional_dim, False)
+base_name=f'edges row{row} center_offset{center_offset}'
+self._plotEdgesOnePlane(recon_plane, base_name)
+if pyip.inputInt('\nContinue (0) or end the search (1): ', min=0, max=1):
+break
+del sinogram_T
+del recon_plane
+center_offset = pyip.inputNum(f'    Enter chosen center offset '+
+f'[{-int(center)}, {int(center)}]: ', min=-int(center), max=int(center))
+return float(center_offset)
+def _reconstructOneTomoStack(self, tomo_stack, thetas, row_bounds=None,
+center_offsets=[], sigma=0.1, rwidth=30, ncore=1, algorithm='gridrec',
+run_secondary_sirt=False, secondary_iter=100):
+"""reconstruct a single tomography stack.
+"""
+# stack order: row,theta,column
+# thetas must be in radians
+# centers_offset: tomography axis shift in pixels relative to column center
+# RV should we remove stripes?
+# https://tomopy.readthedocs.io/en/latest/api/tomopy.prep.stripe.html
+# RV should we remove rings?
+# https://tomopy.readthedocs.io/en/latest/api/tomopy.misc.corr.html
+# RV: Add an option to do (extra) secondary iterations later or to do some sort of convergence test?
+if row_bounds is None:
+row_bounds = [0, tomo_stack.shape[0]]
+else:
+if not (0 <= row_bounds[0] <= row_bounds[1] <= tomo_stack.shape[0]):
+raise ValueError('Illegal row bounds in reconstructOneTomoStack')
+if thetas.size != tomo_stack.shape[1]:
+raise ValueError('theta dimension mismatch in reconstructOneTomoStack')
+if not len(center_offsets):
+centers = np.zeros((row_bounds[1]-row_bounds[0]))
+elif len(center_offsets) == 2:
+centers = np.linspace(center_offsets[0], center_offsets[1],
+row_bounds[1]-row_bounds[0])
+else:
+if center_offsets.size != row_bounds[1]-row_bounds[0]:
+raise ValueError('center_offsets dimension mismatch in reconstructOneTomoStack')
+centers = center_offsets
+centers += tomo_stack.shape[2]/2
+if True:
+tomo_stack = tomopy.prep.stripe.remove_stripe_fw(tomo_stack[row_bounds[0]:row_bounds[1]],
+sigma=sigma, ncore=ncore)
+else:
+tomo_stack = tomo_stack[row_bounds[0]:row_bounds[1]]
+tomo_recon_stack = tomopy.recon(tomo_stack, thetas, centers, sinogram_order=True,
+algorithm=algorithm, ncore=ncore)
+if run_secondary_sirt and secondary_iter > 0:
+#options = {'method':'SIRT_CUDA', 'proj_type':'cuda', 'num_iter':secondary_iter}
+#RV: doesn't work for me: "Error: CUDA error 803: system has unsupported display driver /
+#                          cuda driver combination."
+#options = {'method':'SIRT', 'proj_type':'linear', 'MinConstraint': 0, 'num_iter':secondary_iter}
+#SIRT did not finish while running overnight
+#options = {'method':'SART', 'proj_type':'linear', 'num_iter':secondary_iter}
+options = {'method':'SART', 'proj_type':'linear', 'MinConstraint': 0, 'num_iter':secondary_iter}
+tomo_recon_stack  = tomopy.recon(tomo_stack, thetas, centers, init_recon=tomo_recon_stack,
+options=options, sinogram_order=True, algorithm=tomopy.astra, ncore=ncore)
+if True:
+tomopy.misc.corr.remove_ring(tomo_recon_stack, rwidth=rwidth, out=tomo_recon_stack)
+return tomo_recon_stack
+def findImageFiles(self):
+"""Find all available image files.
+"""
+self.is_valid = True
+# Find dark field images
+dark_field = self.config['dark_field']
+img_start, num_imgs, dark_files = msnc.findImageFiles(
+dark_field['data_path'], self.config['data_filetype'], 'dark field')
+if img_start < 0 or num_imgs < 1:
+logging.error('Unable to find suitable dark field images')
+if dark_field['data_path']:
+self.is_valid = False
+dark_field['num'] = num_imgs
+dark_field['img_start'] = img_start
+logging.info(f'Number of dark field images = {dark_field["num"]}')
+logging.info(f'Dark field image start index = {dark_field["img_start"]}')
+# Find bright field images
+bright_field = self.config['bright_field']
+img_start, num_imgs, bright_files  = msnc.findImageFiles(
+bright_field['data_path'], self.config['data_filetype'], 'bright field')
+if img_start < 0 or num_imgs < 1:
+logging.error('Unable to find suitable bright field images')
+self.is_valid = False
+bright_field['num'] = num_imgs
+bright_field['img_start'] = img_start
+logging.info(f'Number of bright field images = {bright_field["num"]}')
+logging.info(f'Bright field image start index = {bright_field["img_start"]}')
+# Find tomography images
+tomo_stack_files = []
+for stack in self.config['stack_info']['stacks']:
+index = stack['index']
+img_start, num_imgs, tomo_files = msnc.findImageFiles(
+stack['data_path'], self.config['data_filetype'], f'tomography set {index}')
+if img_start < 0 or num_imgs < 1:
+logging.error('Unable to find suitable tomography images')
+self.is_valid = False
+stack['num'] = num_imgs
+stack['img_start'] = img_start
+logging.info(f'Number of tomography images for set {index} = {stack["num"]}')
+logging.info(f'Tomography set {index} image start index = {stack["img_start"]}')
+tomo_stack_files.append(tomo_files)
+del tomo_files
+# Safe updated config
+if self.is_valid:
+self.cf.saveFile(self.config_out)
+return dark_files, bright_files, tomo_stack_files
+def genTomoStacks(self, tdf_files=None, tbf_files=None, tomo_stack_files=None,
+dark_field_pngname=None, bright_field_pngname=None, tomo_field_pngname=None,
+detectorbounds_pngname=None, output_name=None):
+"""Preprocess tomography images.
+"""
+# Try loading any already preprocessed stacks (skip in Galaxy)
+# preprocessed stack order for each one in stack: row,theta,column
+stack_info = self.config['stack_info']
+stacks = stack_info['stacks']
+num_tomo_stacks = stack_info['num']
+assert(num_tomo_stacks == len(self.tomo_stacks))
+available_stacks = [False]*num_tomo_stacks
+if self.galaxy_flag:
+assert(tdf_files is None or isinstance(tdf_files, list))
+assert(isinstance(tbf_files, list))
+assert(isinstance(tomo_stack_files, list))
+assert(num_tomo_stacks == len(tomo_stack_files))
+assert(isinstance(dark_field_pngname, str))
+assert(isinstance(bright_field_pngname, str))
+assert(isinstance(tomo_field_pngname, str))
+assert(isinstance(detectorbounds_pngname, str))
+assert(isinstance(output_name, str))
+else:
+if tdf_files:
+logging.warning('Ignoring tdf_files in genTomoStacks (only for Galaxy)')
+if tbf_files:
+logging.warning('Ignoring tbf_files in genTomoStacks (only for Galaxy)')
+if tomo_stack_files:
+logging.warning('Ignoring tomo_stack_files in genTomoStacks (only for Galaxy)')
+tdf_files, tbf_files, tomo_stack_files = self.findImageFiles()
+if not self.is_valid:
+return
+for i,stack in enumerate(stacks):
+if not self.tomo_stacks[i].size and stack.get('preprocessed', False):
+self.tomo_stacks[i], available_stacks[i] = \
+self._loadTomo('red stack', stack['index'])
+if dark_field_pngname:
+logging.warning('Ignoring dark_field_pngname in genTomoStacks (only for Galaxy)')
+if bright_field_pngname:
+logging.warning('Ignoring bright_field_pngname in genTomoStacks (only for Galaxy)')
+if tomo_field_pngname:
+logging.warning('Ignoring tomo_field_pngname in genTomoStacks (only for Galaxy)')
+if detectorbounds_pngname:
+logging.warning('Ignoring detectorbounds_pngname in genTomoStacks '+
+'(only used in Galaxy)')
+if output_name:
+logging.warning('Ignoring output_name in genTomoStacks '+
+'(only used in Galaxy)')
+# Preprocess any unloaded stacks
+if False in available_stacks:
+logging.debug('Preprocessing tomography images')
+# Check required image files (skip in Galaxy)
+if not self.galaxy_flag:
+self._selectImageRanges(available_stacks)
+if not self.is_valid:
+return
+# Generate dark field
+if tdf_files:
+self._genDark(tdf_files, dark_field_pngname)
+# Generate bright field
+self._genBright(tbf_files, bright_field_pngname)
+# Set vertical detector bounds for image stack
+self._setDetectorBounds(tomo_stack_files, tomo_field_pngname, detectorbounds_pngname)
+# Set zoom and/or theta skip to reduce memory the requirement
+self._setZoomOrSkip()
+# Generate tomography fields
+self._genTomo(tomo_stack_files, available_stacks)
+# Save tomography stack to file
+if self.galaxy_flag:
+t0 = time()
+logging.info(f'Saving preprocessed tomography stack to file ...')
+save_stacks = {f'set_{stack["index"]}':tomo_stack
+for stack,tomo_stack in zip(stacks,self.tomo_stacks)}
+np.savez(output_name, **save_stacks)
+logging.info(f'... done in {time()-t0:.2f} seconds!')
+del available_stacks
+# Adjust sample reference height, update config and save to file
+preprocess = self.config.get('preprocess')
+if preprocess is None:
+img_x_bounds = [0, self.tbf.shape[0]]
+else:
+img_x_bounds = preprocess.get('img_x_bounds', [0, self.tbf.shape[0]])
+pixel_size = self.config['detector']['pixel_size']
+if pixel_size is None:
+raise ValueError('Detector pixel size unavailable')
+pixel_size *= img_x_bounds[0]
+for stack in stacks:
+stack['ref_height'] = stack['ref_height']+pixel_size
+self.cf.saveFile(self.config_out)
+def findCenters(self):
+"""Find rotation axis centers for the tomography stacks.
+"""
+logging.debug('Find centers for tomography stacks')
+stacks = self.config['stack_info']['stacks']
+available_stacks = [stack['index'] for stack in stacks if stack.get('preprocessed', False)]
+logging.debug('Avaliable stacks: {available_stacks}')
+# Check for valid available center stack index
+find_center = self.config.get('find_center')
+if find_center and 'center_stack_index' in find_center:
+center_stack_index = find_center['center_stack_index']
+if (not isinstance(center_stack_index, int) or
+center_stack_index not in available_stacks):
+msnc.illegal_value('center_stack_index', center_stack_index, 'config file')
+else:
+if self.test_mode:
+find_center['completed'] = True
+self.cf.saveFile(self.config_out)
+return
+print('\nFound calibration center offset info for stack '+
+f'{center_stack_index}')
+if pyip.inputYesNo('Do you want to use this again (y/n)? ') == 'yes':
+find_center['completed'] = True
+self.cf.saveFile(self.config_out)
+return
+# Load the required preprocessed stack
+# preprocessed stack order: row,theta,column
+num_tomo_stacks = self.config['stack_info']['num']
+assert(len(stacks) == num_tomo_stacks)
+assert(len(self.tomo_stacks) == num_tomo_stacks)
+if num_tomo_stacks == 1:
+center_stack_index = stacks[0]['index']
+if not self.tomo_stacks[0].size:
+self.tomo_stacks[0], available = self._loadTomo('red stack', center_stack_index,
+required=True)
+center_stack = self.tomo_stacks[0]
+if not center_stack.size:
+logging.error('Unable to load the required preprocessed tomography stack')
+return
+assert(stacks[0].get('preprocessed', False) == True)
+else:
+while True:
+center_stack_index = pyip.inputInt('\nEnter tomography stack index to get '
+f'rotation axis centers {available_stacks}: ')
+while center_stack_index not in available_stacks:
+center_stack_index = pyip.inputInt('\nEnter tomography stack index to get '
+f'rotation axis centers {available_stacks}: ')
+tomo_stack_index = available_stacks.index(center_stack_index)
+if not self.tomo_stacks[tomo_stack_index].size:
+self.tomo_stacks[tomo_stack_index], available = self._loadTomo(
+'red stack', center_stack_index, required=True)
+center_stack = self.tomo_stacks[tomo_stack_index]
+if not center_stack.size:
+logging.error(f'Unable to load the {center_stack_index}th '+
+'preprocessed tomography stack, pick another one')
+else:
+break
+assert(stacks[tomo_stack_index].get('preprocessed', False) == True)
+if find_center is None:
+self.config['find_center'] = {'center_stack_index' : center_stack_index}
+find_center = self.config['find_center']
+# Set thetas (in degrees)
+theta_range = self.config['theta_range']
+theta_start = theta_range['start']
+theta_end = theta_range['end']
+num_theta = theta_range['num']
+num_theta_skip = self.config['preprocess'].get('num_theta_skip', 0)
+thetas_deg = np.linspace(theta_start, theta_end, int(num_theta/(num_theta_skip+1)),
+endpoint=False)
+# Get non-overlapping sample row boundaries
+zoom_perc = self.config['preprocess'].get('zoom_perc', 100)
+pixel_size = self.config['detector']['pixel_size']
+if pixel_size is None:
+raise ValueError('Detector pixel size unavailable')
+eff_pixel_size = 100.*pixel_size/zoom_perc
+logging.debug(f'eff_pixel_size = {eff_pixel_size}')
+tomo_ref_heights = [stack['ref_height'] for stack in stacks]
+if num_tomo_stacks == 1:
+n1 = 0
+height = center_stack.shape[0]*eff_pixel_size
+if pyip.inputYesNo('\nDo you want to reconstruct the full samply height '+
+f'({height:.3f} mm) (y/n)? ') == 'no':
+height = pyip.inputNum('\nEnter the desired reconstructed sample height '+
+f'in mm [0, {height:.3f}]: ', min=0, max=height)
+n1 = int(0.5*(center_stack.shape[0]-height/eff_pixel_size))
+else:
+n1 = int((1.+(tomo_ref_heights[0]+center_stack.shape[0]*eff_pixel_size-
+tomo_ref_heights[1])/eff_pixel_size)/2)
+n2 = center_stack.shape[0]-n1
+logging.info(f'n1 = {n1}, n2 = {n2} (n2-n1) = {(n2-n1)*eff_pixel_size:.3f} mm')
+if not center_stack.size:
+RuntimeError('Center stack not loaded')
+msnc.quickImshow(center_stack[:,0,:], title=f'center stack theta={theta_start}',
+path=self.output_folder, save_fig=self.save_plots, save_only=self.save_plots_only)
+# Get cross sectional diameter in mm
+cross_sectional_dim = center_stack.shape[2]*eff_pixel_size
+logging.debug(f'cross_sectional_dim = {cross_sectional_dim}')
+# Determine center offset at sample row boundaries
+logging.info('Determine center offset at sample row boundaries')
+# Lower row center
+use_row = False
+use_center = False
+row = find_center.get('lower_row')
+if msnc.is_int(row, n1, n2-2):
+msnc.quickImshow(center_stack[:,0,:], title=f'theta={theta_start}', aspect='auto')
+use_row = pyip.inputYesNo('\nCurrent row index for lower center = '
+f'{row}, use this value (y/n)? ')
+if self.save_plots_only:
+msnc.clearFig(f'theta={theta_start}')
+if use_row:
+center_offset = find_center.get('lower_center_offset')
+if msnc.is_num(center_offset):
+use_center = pyip.inputYesNo('Current lower center offset = '+
+f'{center_offset}, use this value (y/n)? ')
+if not use_center:
+if not use_row:
+msnc.quickImshow(center_stack[:,0,:], title=f'theta={theta_start}', aspect='auto')
+row = pyip.inputInt('\nEnter row index to find lower center '+
+f'[[{n1}], {n2-2}]: ', min=n1, max=n2-2, blank=True)
+if row == '':
+row = n1
+if self.save_plots_only:
+msnc.clearFig(f'theta={theta_start}')
+# center_stack order: row,theta,column
+center_offset = self._findCenterOnePlane(center_stack[row,:,:], row, thetas_deg,
+eff_pixel_size, cross_sectional_dim)
+logging.info(f'Lower center offset = {center_offset}')
+# Update config and save to file
+find_center['row_bounds'] = [n1, n2]
+find_center['lower_row'] = row
+find_center['lower_center_offset'] = center_offset
+self.cf.saveFile(self.config_out)
+lower_row = row
+# Upper row center
+use_row = False
+use_center = False
+row = find_center.get('upper_row')
+if msnc.is_int(row, lower_row+1, n2-1):
+msnc.quickImshow(center_stack[:,0,:], title=f'theta={theta_start}', aspect='auto')
+use_row = pyip.inputYesNo('\nCurrent row index for upper center = '
+f'{row}, use this value (y/n)? ')
+if self.save_plots_only:
+msnc.clearFig(f'theta={theta_start}')
+if use_row:
+center_offset = find_center.get('upper_center_offset')
+if msnc.is_num(center_offset):
+use_center = pyip.inputYesNo('Current upper center offset = '+
+f'{center_offset}, use this value (y/n)? ')
+if not use_center:
+if not use_row:
+msnc.quickImshow(center_stack[:,0,:], title=f'theta={theta_start}', aspect='auto')
+row = pyip.inputInt('\nEnter row index to find upper center '+
+f'[{lower_row+1}, [{n2-1}]]: ', min=lower_row+1, max=n2-1, blank=True)
+if row == '':
+row = n2-1
+if self.save_plots_only:
+msnc.clearFig(f'theta={theta_start}')
+# center_stack order: row,theta,column
+center_offset = self._findCenterOnePlane(center_stack[row,:,:], row, thetas_deg,
+eff_pixel_size, cross_sectional_dim)
+logging.info(f'upper_center_offset = {center_offset}')
+del center_stack
+# Update config and save to file
+find_center['upper_row'] = row
+find_center['upper_center_offset'] = center_offset
+find_center['completed'] = True
+self.cf.saveFile(self.config_out)
+def checkCenters(self):
+"""Check centers for the tomography stacks.
+"""
+#RV TODO load all stacks and check at all stack boundaries
+return
+logging.debug('Check centers for tomography stacks')
+center_stack_index = self.config.get('center_stack_index')
+if center_stack_index is None:
+raise ValueError('Unable to read center_stack_index from config')
+center_stack_index = self.tomo_stacks[self.tomo_data_indices.index(center_stack_index)]
+lower_row = self.config.get('lower_row')
+if lower_row is None:
+raise ValueError('Unable to read lower_row from config')
+lower_center_offset = self.config.get('lower_center_offset')
+if lower_center_offset is None:
+raise ValueError('Unable to read lower_center_offset from config')
+upper_row = self.config.get('upper_row')
+if upper_row is None:
+raise ValueError('Unable to read upper_row from config')
+upper_center_offset = self.config.get('upper_center_offset')
+if upper_center_offset is None:
+raise ValueError('Unable to read upper_center_offset from config')
+center_slope = (upper_center_offset-lower_center_offset)/(upper_row-lower_row)
+shift = upper_center_offset-lower_center_offset
+if lower_row == 0:
+logging.warning(f'lower_row == 0: one row offset between both planes')
+else:
+lower_row -= 1
+lower_center_offset -= center_slope
+# stack order: stack,row,theta,column
+if center_stack_index:
+stack1 = self.tomo_stacks[center_stack_index-1]
+stack2 = self.tomo_stacks[center_stack_index]
+if not stack1.size:
+logging.error(f'Unable to load required tomography stack {stack1}')
+elif not stack2.size:
+logging.error(f'Unable to load required tomography stack {stack1}')
+else:
+assert(0 <= lower_row < stack2.shape[0])
+assert(0 <= upper_row < stack1.shape[0])
+plane1 = stack1[upper_row,:]
+plane2 = stack2[lower_row,:]
+for i in range(-2, 3):
+shift_i = shift+2*i
+plane1_shifted = spi.shift(plane2, [0, shift_i])
+msnc.quickPlot((plane1[0,:],), (plane1_shifted[0,:],),
+title=f'stacks {stack1} {stack2} shifted {2*i} theta={self.start_theta}',
+path=self.output_folder, save_fig=self.save_plots,
+save_only=self.save_plots_only)
+if center_stack_index < self.num_tomo_stacks-1:
+stack1 = self.tomo_stacks[center_stack_index]
+stack2 = self.tomo_stacks[center_stack_index+1]
+if not stack1.size:
+logging.error(f'Unable to load required tomography stack {stack1}')
+elif not stack2.size:
+logging.error(f'Unable to load required tomography stack {stack1}')
+else:
+assert(0 <= lower_row < stack2.shape[0])
+assert(0 <= upper_row < stack1.shape[0])
+plane1 = stack1[upper_row,:]
+plane2 = stack2[lower_row,:]
+for i in range(-2, 3):
+shift_i = -shift+2*i
+plane1_shifted = spi.shift(plane2, [0, shift_i])
+msnc.quickPlot((plane1[0,:],), (plane1_shifted[0,:],),
+title=f'stacks {stack1} {stack2} shifted {2*i} theta={start_theta}',
+path=self.output_folder, save_fig=self.save_plots,
+save_only=self.save_plots_only)
+del plane1, plane2, plane1_shifted
+# Update config file
+self.config = msnc.update('config.txt', 'check_centers', True, 'find_centers')
+def reconstructTomoStacks(self):
+"""Reconstruct tomography stacks.
+"""
+logging.debug('Reconstruct tomography stacks')
+# Get rotation axis rows and centers
+find_center = self.config['find_center']
+lower_row = find_center.get('lower_row')
+if lower_row is None:
+logging.error('Unable to read lower_row from config')
+return
+lower_center_offset = find_center.get('lower_center_offset')
+if lower_center_offset is None:
+logging.error('Unable to read lower_center_offset from config')
+return
+upper_row = find_center.get('upper_row')
+if upper_row is None:
+logging.error('Unable to read upper_row from config')
+return
+upper_center_offset = find_center.get('upper_center_offset')
+if upper_center_offset is None:
+logging.error('Unable to read upper_center_offset from config')
+return
+logging.debug(f'lower_row = {lower_row} upper_row = {upper_row}')
+assert(lower_row < upper_row)
+center_slope = (upper_center_offset-lower_center_offset)/(upper_row-lower_row)
+# Set thetas (in radians)
+theta_range = self.config['theta_range']
+theta_start = theta_range['start']
+theta_end = theta_range['end']
+num_theta = theta_range['num']
+num_theta_skip = self.config['preprocess'].get('num_theta_skip', 0)
+thetas = np.radians(np.linspace(theta_start, theta_end,
+int(num_theta/(num_theta_skip+1)), endpoint=False))
+# Reconstruct tomo stacks
+zoom_perc = self.config['preprocess'].get('zoom_perc', 100)
+if zoom_perc == 100:
+basetitle = 'recon stack full'
+else:
+basetitle = f'recon stack {zoom_perc}p'
+load_error = False
+stacks = self.config['stack_info']['stacks']
+for i,stack in enumerate(stacks):
+# Check if stack can be loaded
+# reconstructed stack order for each one in stack : row/z,x,y
+# preprocessed stack order for each one in stack: row,theta,column
+index = stack['index']
+available = False
+if stack.get('reconstructed', False):
+self.tomo_recon_stacks[i], available = self._loadTomo('recon stack', index)
+if self.tomo_recon_stacks[i].size or available:
+if self.tomo_stacks[i].size:
+self.tomo_stacks[i] = np.array([])
+assert(stack.get('reconstructed', False) == True)
+continue
+else:
+stack['reconstructed'] = False
+if not self.tomo_stacks[i].size:
+self.tomo_stacks[i], available = self._loadTomo('red stack', index,
+required=True)
+if not self.tomo_stacks[i].size:
+logging.error(f'Unable to load tomography stack {index} for reconstruction')
+load_error = True
+continue
+assert(0 <= lower_row < upper_row < self.tomo_stacks[i].shape[0])
+center_offsets = [lower_center_offset-lower_row*center_slope,
+upper_center_offset+(self.tomo_stacks[i].shape[0]-1-upper_row)*center_slope]
+t0 = time()
+self.tomo_recon_stacks[i]= self._reconstructOneTomoStack(self.tomo_stacks[i],
+thetas, center_offsets=center_offsets, sigma=0.1, ncore=self.ncore,
+algorithm='gridrec', run_secondary_sirt=True, secondary_iter=25)
+logging.info(f'Reconstruction of stack {index} took {time()-t0:.2f} seconds!')
+if not self.test_mode:
+row_slice = int(self.tomo_stacks[i].shape[0]/2)
+title = f'{basetitle} {index} slice{row_slice}'
+msnc.quickImshow(self.tomo_recon_stacks[i][row_slice,:,:], title=title,
+path=self.output_folder, save_fig=self.save_plots,
+save_only=self.save_plots_only)
+msnc.quickPlot(self.tomo_recon_stacks[i]
+[row_slice,int(self.tomo_recon_stacks[i].shape[2]/2),:],
+title=f'{title} cut{int(self.tomo_recon_stacks[i].shape[2]/2)}',
+path=self.output_folder, save_fig=self.save_plots,
+save_only=self.save_plots_only)
+self._saveTomo('recon stack', self.tomo_recon_stacks[i], index)
+#                else:
+#                    np.savetxt(self.output_folder+f'recon_stack_{index}.txt',
+#                            self.tomo_recon_stacks[i][row_slice,:,:], fmt='%.6e')
+self.tomo_stacks[i] = np.array([])
+# Update config and save to file
+stack['reconstructed'] = True
+self.cf.saveFile(self.config_out)
+def combineTomoStacks(self):
+"""Combine the reconstructed tomography stacks.
+"""
+# stack order: stack,row(z),x,y
+logging.debug('Combine reconstructed tomography stacks')
+# Load any unloaded reconstructed stacks
+stacks = self.config['stack_info']['stacks']
+for i,stack in enumerate(stacks):
+if not self.tomo_recon_stacks[i].size and stack.get('reconstructed', False):
+self.tomo_recon_stacks[i], available = self._loadTomo('recon stack',
+stack['index'], required=True)
+if not self.tomo_recon_stacks[i].size or not available:
+logging.error(f'Unable to load reconstructed stack {stack["index"]}')
+return
+if i:
+if (self.tomo_recon_stacks[i].shape[1] != self.tomo_recon_stacks[0].shape[1] or
+self.tomo_recon_stacks[i].shape[1] != self.tomo_recon_stacks[0].shape[1]):
+logging.error('Incompatible reconstructed tomography stack dimensions')
+return
+# Get center stack boundaries
+row_bounds = self.config['find_center']['row_bounds']
+if not msnc.is_index_range(row_bounds, 0, self.tomo_recon_stacks[0].shape[0]):
+msnc.illegal_value('row_bounds', row_bounds, 'config file')
+return
+# Selecting xy bounds
+tomosum = 0
+#RV FIX :=
+#[tomosum := tomosum+np.sum(tomo_recon_stack, axis=(0,2)) for tomo_recon_stack in
+#        self.tomo_recon_stacks]
+combine_stacks =self.config.get('combine_stacks')
+if combine_stacks and 'x_bounds' in combine_stacks:
+x_bounds = combine_stacks['x_bounds']
+if not msnc.is_index_range(x_bounds, 0, self.tomo_recon_stacks[0].shape[1]):
+msnc.illegal_value('x_bounds', x_bounds, 'config file')
+elif not self.test_mode:
+msnc.quickPlot(tomosum, title='recon stack sum yz')
+if pyip.inputYesNo('\nCurrent image x-bounds: '+
+f'[{x_bounds[0]}, {x_bounds[1]}], use these values ([y]/n)? ',
+blank=True) == 'no':
+if pyip.inputYesNo(
+'Do you want to change the image x-bounds ([y]/n)? ',
+blank=True) == 'no':
+x_bounds = [0, self.tomo_recon_stacks[0].shape[1]]
+else:
+x_bounds = msnc.selectArrayBounds(tomosum, title='recon stack sum yz')
+else:
+msnc.quickPlot(tomosum, title='recon stack sum yz')
+if pyip.inputYesNo('\nDo you want to change the image x-bounds (y/n)? ') == 'no':
+x_bounds = [0, self.tomo_recon_stacks[0].shape[1]]
+else:
+x_bounds = msnc.selectArrayBounds(tomosum, title='recon stack sum yz')
+if False and self.test_mode:
+np.savetxt(self.output_folder+'recon_stack_sum_yz.txt',
+tomosum[x_bounds[0]:x_bounds[1]], fmt='%.6e')
+if self.save_plots_only:
+msnc.clearFig('recon stack sum yz')
+tomosum = 0
+#RV FIX :=
+#[tomosum := tomosum+np.sum(tomo_recon_stack, axis=(0,1)) for tomo_recon_stack in
+#        self.tomo_recon_stacks]
+if combine_stacks and 'y_bounds' in combine_stacks:
+y_bounds = combine_stacks['y_bounds']
+if not msnc.is_index_range(x_bounds, 0, self.tomo_recon_stacks[0].shape[1]):
+msnc.illegal_value('y_bounds', y_bounds, 'config file')
+elif not self.test_mode:
+msnc.quickPlot(tomosum, title='recon stack sum xz')
+if pyip.inputYesNo('\nCurrent image y-bounds: '+
+f'[{y_bounds[0]}, {y_bounds[1]}], use these values ([y]/n)? ',
+blank=True) == 'no':
+if pyip.inputYesNo(
+'Do you want to change the image y-bounds ([y]/n)? ',
+blank=True) == 'no':
+y_bounds = [0, self.tomo_recon_stacks[0].shape[1]]
+else:
+y_bounds = msnc.selectArrayBounds(tomosum, title='recon stack sum yz')
+else:
+msnc.quickPlot(tomosum, title='recon stack sum xz')
+if pyip.inputYesNo('\nDo you want to change the image y-bounds (y/n)? ') == 'no':
+y_bounds = [0, self.tomo_recon_stacks[0].shape[1]]
+else:
+y_bounds = msnc.selectArrayBounds(tomosum, title='recon stack sum xz')
+if False and self.test_mode:
+np.savetxt(self.output_folder+'recon_stack_sum_xz.txt',
+tomosum[y_bounds[0]:y_bounds[1]], fmt='%.6e')
+if self.save_plots_only:
+msnc.clearFig('recon stack sum xz')
+# Combine reconstructed tomography stacks
+logging.info(f'Combining reconstructed stacks ...')
+t0 = time()
+num_tomo_stacks = self.config['stack_info']['num']
+if num_tomo_stacks == 1:
+low_bound = row_bounds[0]
+else:
+low_bound = 0
+tomo_recon_combined = self.tomo_recon_stacks[0][low_bound:row_bounds[1]:,
+x_bounds[0]:x_bounds[1],y_bounds[0]:y_bounds[1]]
+if num_tomo_stacks > 2:
+tomo_recon_combined = np.concatenate([tomo_recon_combined]+
+[self.tomo_recon_stacks[i][row_bounds[0]:row_bounds[1],
+x_bounds[0]:x_bounds[1],y_bounds[0]:y_bounds[1]]
+for i in range(1, num_tomo_stacks-1)])
+if num_tomo_stacks > 1:
+tomo_recon_combined = np.concatenate([tomo_recon_combined]+
+[self.tomo_recon_stacks[num_tomo_stacks-1][row_bounds[0]:,
+x_bounds[0]:x_bounds[1],y_bounds[0]:y_bounds[1]]])
+logging.info(f'... done in {time()-t0:.2f} seconds!')
+tomosum = np.sum(tomo_recon_combined, axis=(1,2))
+if self.test_mode:
+zoom_perc = self.config['preprocess'].get('zoom_perc', 100)
+filename = 'recon combined sum xy'
+if zoom_perc is None or zoom_perc == 100:
+filename += ' fullres.dat'
+else:
+filename += f' {zoom_perc}p.dat'
+msnc.quickPlot(tomosum, title='recon combined sum xy',
+path=self.output_folder, save_fig=self.save_plots,
+save_only=self.save_plots_only)
+if False:
+np.savetxt(self.output_folder+'recon_combined_sum_xy.txt',
+tomosum, fmt='%.6e')
+np.savetxt(self.output_folder+'recon_combined.txt',
+tomo_recon_combined[int(tomo_recon_combined.shape[0]/2),:,:], fmt='%.6e')
+combine_stacks =self.config.get('combine_stacks')
+# Update config and save to file
+if combine_stacks:
+combine_stacks['x_bounds'] = x_bounds
+combine_stacks['y_bounds'] = y_bounds
+else:
+self.config['combine_stacks'] = {'x_bounds' : x_bounds, 'y_bounds' : y_bounds}
+self.cf.saveFile(self.config_out)
+return
+msnc.quickPlot(tomosum, title='recon combined sum xy')
+if pyip.inputYesNo(
+'\nDo you want to change the image z-bounds (y/[n])? ',
+blank=True) != 'yes':
+z_bounds = [0, tomo_recon_combined.shape[0]]
+else:
+z_bounds = msnc.selectArrayBounds(tomosum, title='recon combined sum xy')
+if z_bounds[0] != 0 or z_bounds[1] != tomo_recon_combined.shape[0]:
+tomo_recon_combined = tomo_recon_combined[z_bounds[0]:z_bounds[1],:,:]
+logging.info(f'tomo_recon_combined.shape = {tomo_recon_combined.shape}')
+if self.save_plots_only:
+msnc.clearFig('recon combined sum xy')
+# Plot center slices
+msnc.quickImshow(tomo_recon_combined[int(tomo_recon_combined.shape[0]/2),:,:],
+title=f'recon combined xslice{int(tomo_recon_combined.shape[0]/2)}',
+path=self.output_folder, save_fig=self.save_plots,
+save_only=self.save_plots_only)
+msnc.quickImshow(tomo_recon_combined[:,int(tomo_recon_combined.shape[1]/2),:],
+title=f'recon combined yslice{int(tomo_recon_combined.shape[1]/2)}',
+path=self.output_folder, save_fig=self.save_plots,
+save_only=self.save_plots_only)
+msnc.quickImshow(tomo_recon_combined[:,:,int(tomo_recon_combined.shape[2]/2)],
+title=f'recon combined zslice{int(tomo_recon_combined.shape[2]/2)}',
+path=self.output_folder, save_fig=self.save_plots,
+save_only=self.save_plots_only)
+# Save combined reconstructed tomo stacks
+base_name = 'recon combined'
+combined_stacks = []
+for stack in stacks:
+base_name += f' {stack["index"]}'
+combined_stacks.append(stack['index'])
+self._saveTomo(base_name, tomo_recon_combined)
+# Update config and save to file
+if combine_stacks:
+combine_stacks['x_bounds'] = x_bounds
+combine_stacks['y_bounds'] = y_bounds
+combine_stacks['stacks'] = combined_stacks
+else:
+self.config['combine_stacks'] = {'x_bounds' : x_bounds, 'y_bounds' : y_bounds,
+'stacks' : combined_stacks}
+self.cf.saveFile(self.config_out)
+def runTomo(config_file=None, config_dict=None, output_folder='.', log_level='INFO',
+test_mode=False):
+"""Run a tomography analysis.
+"""
+# Instantiate Tomo object
+tomo = Tomo(config_file=config_file, output_folder=output_folder, log_level=log_level,
+test_mode=test_mode)
+if not tomo.is_valid:
+raise ValueError('Invalid config and/or detector file provided.')
+# Preprocess the image files
+num_tomo_stacks = tomo.config['stack_info']['num']
+assert(num_tomo_stacks == len(tomo.tomo_stacks))
+preprocess = tomo.config.get('preprocess', None)
+preprocessed_stacks = []
+if preprocess:
+preprocessed_stacks = [stack['index'] for stack in tomo.config['stack_info']['stacks']
+if stack.get('preprocessed', False)]
+if not len(preprocessed_stacks):
+tomo.genTomoStacks()
+if not tomo.is_valid:
+IOError('Unable to load all required image files.')
+find_center = tomo.config.get('find_center')
+if find_center and find_center.get('completed', False):
+center_stack_index = find_center['center_stack_index']
+if not center_stack_index in preprocessed_stacks:
+find_center['completed'] = False
+#RV FIX
+#        tomo.config.pop('check_center', 'check_center not found')
+#        combined_stacks = tomo.config.get('combined_stacks')
+#        if combined_stacks:
+#            combined_stacks['completed'] = False
+tomo.cf.saveFile(tomo.config_out)
+# Find centers
+find_center = tomo.config.get('find_center')
+if find_center is None or not find_center.get('completed', False):
+tomo.findCenters()
+# Check centers
+#if num_tomo_stacks > 1 and not tomo.config.get('check_centers', False):
+#    tomo.checkCenters()
+# Reconstruct tomography stacks
+if len(tomo.config.get('reconstructed_stacks', [])) != tomo.config['stack_info']['num']:
+tomo.reconstructTomoStacks()
+# Combine reconstructed tomography stacks
+combined_stacks = tomo.config.get('combined_stacks')
+if combined_stacks is None or not combined_stacks.get('completed', False):
+tomo.combineTomoStacks()
+#%%============================================================================
+if __name__ == '__main__':
+# Parse command line arguments
+arguments = sys.argv[1:]
+config_file = None
+output_folder = '.'
+log_level = 'INFO'
+test_mode = False
+try:
+opts, args = getopt.getopt(arguments,"hc:o:l:t")
+except getopt.GetoptError:
+print('usage: tomo.py -c <config_file> -o <output_folder> -l <log_level> -t')
+sys.exit(2)
+for opt, arg in opts:
+if opt == '-h':
+print('usage: tomo.py -c <config_file> -o <output_folder> -l <log_level> -t')
+sys.exit()
+elif opt in ("-c"):
+config_file = arg
+elif opt in ("-o"):
+output_folder = arg
+elif opt in ("-l"):
+log_level = arg
+elif opt in ("-t"):
+test_mode = True
+if config_file is None:
+if os.path.isfile('config.yaml'):
+config_file = 'config.yaml'
+else:
+config_file = 'config.txt'
+# Set basic log configuration
+logging_format = '%(asctime)s : %(levelname)s - %(module)s : %(funcName)s - %(message)s'
+if not test_mode:
+level = getattr(logging, log_level.upper(), None)
+if not isinstance(level, int):
+raise ValueError(f'Invalid log_level: {log_level}')
+logging.basicConfig(format=logging_format, level=level, force=True,
+handlers=[logging.StreamHandler()])
+logging.debug(f'config_file = {config_file}')
+logging.debug(f'output_folder = {output_folder}')
+logging.debug(f'log_level = {log_level}')
+logging.debug(f'test_mode = {test_mode}')
+# Run tomography analysis
+runTomo(config_file=config_file, output_folder=output_folder, log_level=log_level,
+test_mode=test_mode)
+#%%============================================================================
+input('Press any key to continue')
+#%%============================================================================

Mercurial > repos > rv43 > tomo

comparison tomo.py @ 0:cb1b0d757704 draft