--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/tomo.py	Tue Mar 29 16:10:16 2022 +0000
@@ -0,0 +1,2037 @@
+#!/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')
+#%%============================================================================