--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/.shed.yml	Thu Mar 24 16:59:59 2022 +0000
@@ -0,0 +1,6 @@
+categories:
+- Structural Materials Analysis
+description: Preprocess tomography images
+name: tomo_setup
+owner: ximgchess
+remote_repository_url: https://github.com/ximg-chess/galaxytools/tools/tomo

--- a/msnc_tools.py	Thu Mar 24 16:57:45 2022 +0000
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,853 +0,0 @@
-#!/usr/bin/env python3
-
-# -*- coding: utf-8 -*-
-"""
-Created on Mon Dec  6 15:36:22 2021
-
-@author: rv43
-"""
-
-import logging
-
-import os
-import sys
-import re
-import yaml
-import h5py
-import pyinputplus as pyip
-import numpy as np
-import imageio as img
-import matplotlib.pyplot as plt
-from time import time
-from ast import literal_eval
-from lmfit.models import StepModel, RectangleModel
-
-def depth_list(L): return isinstance(L, list) and max(map(depth_list, L))+1
-def depth_tuple(T): return isinstance(T, tuple) and max(map(depth_tuple, T))+1
-
-def is_int(v, v_min=None, v_max=None):
-    """Value is an integer in range v_min <= v <= v_max.
-    """
-    if not isinstance(v, int):
-        return False
-    if (v_min != None and v < v_min) or (v_max != None and v > v_max):
-        return False
-    return True
-
-def is_num(v, v_min=None, v_max=None):
-    """Value is a number in range v_min <= v <= v_max.
-    """
-    if not isinstance(v, (int,float)):
-        return False
-    if (v_min != None and v < v_min) or (v_max != None and v > v_max):
-        return False
-    return True
-
-def is_index(v, v_min=0, v_max=None):
-    """Value is an array index in range v_min <= v < v_max.
-    """
-    if not isinstance(v, int):
-        return False
-    if v < v_min or (v_max != None and v >= v_max):
-        return False
-    return True
-
-def is_index_range(v, v_min=0, v_max=None):
-    """Value is an array index range in range v_min <= v[0] <= v[1] < v_max.
-    """
-    if not (isinstance(v, list) and len(v) == 2 and isinstance(v[0], int) and
-            isinstance(v[1], int)):
-        return False
-    if not 0 <= v[0] < v[1] or (v_max != None and v[1] >= v_max):
-        return False
-    return True
-
-def illegal_value(name, value, location=None, exit_flag=False):
-    if not isinstance(location, str):
-        location = ''
-    else:
-        location = f'in {location} '
-    if isinstance(name, str):
-        logging.error(f'Illegal value for {name} {location}({value}, {type(value)})')
-    else:
-        logging.error(f'Illegal value {location}({value}, {type(value)})')
-    if exit_flag:
-        exit(1)
-
-def get_trailing_int(string):
-    indexRegex = re.compile(r'\d+$')
-    mo = indexRegex.search(string)
-    if mo == None:
-        return None
-    else:
-        return int(mo.group())
-
-def findImageFiles(path, filetype, name=None, select_range=False, num_required=None):
-    if isinstance(name, str):
-        name = f' {name} '
-    else:
-        name = ' '
-    # Find available index range
-    if filetype == 'tif':
-        if not isinstance(path, str) and not os.path.isdir(path):
-            illegal_value('path', path, 'findImageRange')
-            return -1, 0, []
-        indexRegex = re.compile(r'\d+')
-        # At this point only tiffs
-        files = sorted([f for f in os.listdir(path) if os.path.isfile(os.path.join(path, f)) and
-                f.endswith('.tif') and indexRegex.search(f)])
-        num_imgs = len(files)
-        if num_imgs < 1:
-            logging.warning('No available'+name+'files')
-            return -1, 0, []
-        first_index = indexRegex.search(files[0]).group()
-        last_index = indexRegex.search(files[-1]).group()
-        if first_index == None or last_index == None:
-            logging.error('Unable to find correctly indexed'+name+'images')
-            return -1, 0, []
-        first_index = int(first_index)
-        last_index = int(last_index)
-        if num_imgs != last_index-first_index+1:
-            logging.error('Non-consecutive set of indices for'+name+'images')
-            return -1, 0, []
-        paths = [os.path.join(path, f) for f in files]
-    elif filetype == 'h5':
-        if not isinstance(path, str) or not os.path.isfile(path):
-            illegal_value('path', path, 'findImageRange')
-            return -1, 0, []
-        # At this point only h5 in alamo2 detector style
-        first_index = 0
-        with h5py.File(path, 'r') as f:
-            num_imgs = f['entry/instrument/detector/data'].shape[0]
-            last_index = num_imgs-1
-        paths = [path]
-    else:
-        illegal_value('filetype', filetype, 'findImageRange')
-        return -1, 0, []
-    logging.debug('\nNumber of available'+name+f'images: {num_imgs}')
-    logging.debug('Index range of available'+name+f'images: [{first_index}, '+
-            f'{last_index}]')
-
-    return first_index, num_imgs, paths
-
-def selectImageRange(first_index, offset, num_imgs, name=None, num_required=None):
-    if isinstance(name, str):
-        name = f' {name} '
-    else:
-        name = ' '
-    # Check existing values
-    use_input = 'no'
-    if (is_int(first_index, 0) and is_int(offset, 0) and is_int(num_imgs, 1)):
-        if offset < 0:
-            use_input = pyip.inputYesNo('\nCurrent'+name+f'first index = {first_index}, '+
-                    'use this value ([y]/n)? ', blank=True)
-        else:
-            use_input = pyip.inputYesNo('\nCurrent'+name+'first index/offset = '+
-                    f'{first_index}/{offset}, use these values ([y]/n)? ',
-                    blank=True)
-        if use_input != 'no':
-            use_input = pyip.inputYesNo('Current number of'+name+'images = '+
-                    f'{num_imgs}, use this value ([y]/n)? ',
-                    blank=True)
-    if use_input == 'yes':
-        return first_index, offset, num_imgs
-
-    # Check range against requirements
-    if num_imgs < 1:
-        logging.warning('No available'+name+'images')
-        return -1, -1, 0
-    if num_required == None:
-        if num_imgs == 1:
-            return first_index, 0, 1
-    else:
-        if not is_int(num_required, 1):
-            illegal_value('num_required', num_required, 'selectImageRange')
-            return -1, -1, 0
-        if num_imgs < num_required:
-            logging.error('Unable to find the required'+name+
-                    f'images ({num_imgs} out of {num_required})')
-            return -1, -1, 0
-
-    # Select index range
-    if num_required == None:
-        last_index = first_index+num_imgs
-        use_all = f'Use all ([{first_index}, {last_index}])'
-        pick_offset = 'Pick a first index offset and a number of images'
-        pick_bounds = 'Pick the first and last index'
-        menuchoice = pyip.inputMenu([use_all, pick_offset, pick_bounds], numbered=True)
-        if menuchoice == use_all:
-            offset = 0
-        elif menuchoice == pick_offset:
-            offset = pyip.inputInt('Enter the first index offset'+
-                    f' [0, {last_index-first_index}]: ', min=0, max=last_index-first_index)
-            first_index += offset
-            if first_index == last_index:
-                num_imgs = 1
-            else:
-                num_imgs = pyip.inputInt(f'Enter the number of images [1, {num_imgs}]: ',
-                        min=1, max=num_imgs)
-        else:
-            offset = pyip.inputInt(f'Enter the first index [{first_index}, {last_index}]: ',
-                    min=first_index, max=last_index)-first_index
-            first_index += offset
-            num_imgs = pyip.inputInt(f'Enter the last index [{first_index}, {last_index}]: ',
-                    min=first_index, max=last_index)-first_index+1
-    else:
-        use_all = f'Use ([{first_index}, {first_index+num_required-1}])'
-        pick_offset = 'Pick the first index offset'
-        menuchoice = pyip.inputMenu([use_all, pick_offset], numbered=True)
-        offset = 0
-        if menuchoice == pick_offset:
-            offset = pyip.inputInt('Enter the first index offset'+
-                    f'[0, {num_imgs-num_required}]: ', min=0, max=num_imgs-num_required)
-            first_index += offset
-        num_imgs = num_required
-
-    return first_index, offset, num_imgs
-
-def loadImage(f, img_x_bounds=None, img_y_bounds=None):
-    """Load a single image from file.
-    """
-    if not os.path.isfile(f):
-        logging.error(f'Unable to load {f}')
-        return None
-    img_read = img.imread(f)
-    if not img_x_bounds:
-        img_x_bounds = [0, img_read.shape[0]]
-    else:
-        if (not isinstance(img_x_bounds, list) or len(img_x_bounds) != 2 or 
-                not (0 <= img_x_bounds[0] < img_x_bounds[1] <= img_read.shape[0])):
-            logging.error(f'inconsistent row dimension in {f}')
-            return None
-    if not img_y_bounds:
-        img_y_bounds = [0, img_read.shape[1]]
-    else:
-        if (not isinstance(img_y_bounds, list) or len(img_y_bounds) != 2 or 
-                not (0 <= img_y_bounds[0] < img_y_bounds[1] <= img_read.shape[0])):
-            logging.error(f'inconsistent column dimension in {f}')
-            return None
-    return img_read[img_x_bounds[0]:img_x_bounds[1],img_y_bounds[0]:img_y_bounds[1]]
-
-def loadImageStack(files, filetype, img_offset, num_imgs, num_img_skip=0,
-        img_x_bounds=None, img_y_bounds=None):
-    """Load a set of images and return them as a stack.
-    """
-    logging.debug(f'img_offset = {img_offset}')
-    logging.debug(f'num_imgs = {num_imgs}')
-    logging.debug(f'num_img_skip = {num_img_skip}')
-    logging.debug(f'\nfiles:\n{files}\n')
-    img_stack = np.array([])
-    if filetype == 'tif':
-        img_read_stack = []
-        i = 1
-        t0 = time()
-        for f in files[img_offset:img_offset+num_imgs:num_img_skip+1]:
-            if not i%20:
-                logging.info(f'    loading {i}/{num_imgs}: {f}')
-            else:
-                logging.debug(f'    loading {i}/{num_imgs}: {f}')
-            img_read = loadImage(f, img_x_bounds, img_y_bounds)
-            img_read_stack.append(img_read)
-            i += num_img_skip+1
-        img_stack = np.stack([img_read for img_read in img_read_stack])
-        logging.info(f'... done in {time()-t0:.2f} seconds!')
-        logging.debug(f'img_stack shape = {np.shape(img_stack)}')
-        del img_read_stack, img_read
-    elif filetype == 'h5':
-        if not isinstance(files[0], str) and not os.path.isfile(files[0]):
-            illegal_value('files[0]', files[0], 'loadImageStack')
-            return img_stack
-        t0 = time()
-        with h5py.File(files[0], 'r') as f:
-            shape = f['entry/instrument/detector/data'].shape
-            if len(shape) != 3:
-                logging.error(f'inconsistent dimensions in {files[0]}')
-            if not img_x_bounds:
-                img_x_bounds = [0, shape[1]]
-            else:
-                if (not isinstance(img_x_bounds, list) or len(img_x_bounds) != 2 or 
-                        not (0 <= img_x_bounds[0] < img_x_bounds[1] <= shape[1])):
-                    logging.error(f'inconsistent row dimension in {files[0]}')
-            if not img_y_bounds:
-                img_y_bounds = [0, shape[2]]
-            else:
-                if (not isinstance(img_y_bounds, list) or len(img_y_bounds) != 2 or 
-                        not (0 <= img_y_bounds[0] < img_y_bounds[1] <= shape[2])):
-                    logging.error(f'inconsistent column dimension in {files[0]}')
-            img_stack = f.get('entry/instrument/detector/data')[
-                    img_offset:img_offset+num_imgs:num_img_skip+1,
-                    img_x_bounds[0]:img_x_bounds[1],img_y_bounds[0]:img_y_bounds[1]]
-        logging.info(f'... done in {time()-t0:.2f} seconds!')
-    else:
-        illegal_value('filetype', filetype, 'findImageRange')
-    return img_stack
-
-def clearFig(title):
-    if not isinstance(title, str):
-        illegal_value('title', title, 'clearFig')
-        return
-    plt.close(fig=re.sub(r"\s+", '_', title))
-
-def quickImshow(a, title=None, path=None, name=None, save_fig=False, save_only=False,
-            clear=True, **kwargs):
-    if title != None and not isinstance(title, str):
-        illegal_value('title', title, 'quickImshow')
-        return
-    if path is not None and not isinstance(path, str):
-        illegal_value('path', path, 'quickImshow')
-        return
-    if not isinstance(save_fig, bool):
-        illegal_value('save_fig', save_fig, 'quickImshow')
-        return
-    if not isinstance(save_only, bool):
-        illegal_value('save_only', save_only, 'quickImshow')
-        return
-    if not isinstance(clear, bool):
-        illegal_value('clear', clear, 'quickImshow')
-        return
-    if not title:
-        title='quick_imshow'
-    else:
-        title = re.sub(r"\s+", '_', title)
-    if name is None:
-        if path is None:
-            path = f'{title}.png'
-        else:
-            path = f'{path}/{title}.png'
-    else:
-        if path is None:
-            path = name
-        else:
-            path = f'{path}/{name}'
-    if clear:
-        plt.close(fig=title)
-    if save_only:
-        plt.figure(title)
-        plt.imshow(a, **kwargs)
-        plt.savefig(path)
-        plt.close(fig=title)
-        #plt.imsave(f'{title}.png', a, **kwargs)
-    else:
-        plt.ion()
-        plt.figure(title)
-        plt.imshow(a, **kwargs)
-        if save_fig:
-            plt.savefig(path)
-        plt.pause(1)
-
-def quickPlot(*args, title=None, path=None, name=None, save_fig=False, save_only=False,
-        clear=True, **kwargs):
-    if title != None and not isinstance(title, str):
-        illegal_value('title', title, 'quickPlot')
-        return
-    if path is not None and not isinstance(path, str):
-        illegal_value('path', path, 'quickPlot')
-        return
-    if not isinstance(save_fig, bool):
-        illegal_value('save_fig', save_fig, 'quickPlot')
-        return
-    if not isinstance(save_only, bool):
-        illegal_value('save_only', save_only, 'quickPlot')
-        return
-    if not isinstance(clear, bool):
-        illegal_value('clear', clear, 'quickPlot')
-        return
-    if not title:
-        title = 'quick_plot'
-    else:
-        title = re.sub(r"\s+", '_', title)
-    if name is None:
-        if path is None:
-            path = f'{title}.png'
-        else:
-            path = f'{path}/{title}.png'
-    else:
-        if path is None:
-            path = name
-        else:
-            path = f'{path}/{name}'
-    if clear:
-        plt.close(fig=title)
-    if save_only:
-        plt.figure(title)
-        if depth_tuple(args) > 1:
-           for y in args:
-               plt.plot(*y, **kwargs)
-        else:
-            plt.plot(*args, **kwargs)
-        plt.savefig(path)
-        plt.close(fig=title)
-    else:
-        plt.ion()
-        plt.figure(title)
-        if depth_tuple(args) > 1:
-           for y in args:
-               plt.plot(*y, **kwargs)
-        else:
-            plt.plot(*args, **kwargs)
-        if save_fig:
-            plt.savefig(path)
-        plt.pause(1)
-
-def selectArrayBounds(a, x_low=None, x_upp=None, num_x_min=None,
-        title='select array bounds'):
-    """Interactively select the lower and upper data bounds for a numpy array.
-    """
-    if not isinstance(a, np.ndarray) or a.ndim != 1:
-        logging.error('Illegal array type or dimension in selectArrayBounds')
-        return None
-    if num_x_min == None:
-        num_x_min = 1
-    else:
-        if num_x_min < 2 or num_x_min > a.size:
-            logging.warning('Illegal input for num_x_min in selectArrayBounds, input ignored')
-            num_x_min = 1
-    if x_low == None:
-        x_min = 0
-        x_max = a.size
-        x_low_max = a.size-num_x_min
-        while True:
-            quickPlot(range(x_min, x_max), a[x_min:x_max], title=title)
-            zoom_flag = pyip.inputInt('Set lower data bound (0) or zoom in (1)?: ',
-                    min=0, max=1)
-            if zoom_flag:
-                x_min = pyip.inputInt(f'    Set lower zoom index [0, {x_low_max}]: ',
-                        min=0, max=x_low_max)
-                x_max = pyip.inputInt(f'    Set upper zoom index [{x_min+1}, {x_low_max+1}]: ',
-                        min=x_min+1, max=x_low_max+1)
-            else:
-                x_low = pyip.inputInt(f'    Set lower data bound [0, {x_low_max}]: ',
-                        min=0, max=x_low_max)
-                break
-    else:
-        if not is_int(x_low, 0, a.size-num_x_min):
-            illegal_value('x_low', x_low, 'selectArrayBounds')
-            return None
-    if x_upp == None:
-        x_min = x_low+num_x_min
-        x_max = a.size
-        x_upp_min = x_min
-        while True:
-            quickPlot(range(x_min, x_max), a[x_min:x_max], title=title)
-            zoom_flag = pyip.inputInt('Set upper data bound (0) or zoom in (1)?: ',
-                    min=0, max=1)
-            if zoom_flag:
-                x_min = pyip.inputInt(f'    Set upper zoom index [{x_upp_min}, {a.size-1}]: ',
-                        min=x_upp_min, max=a.size-1)
-                x_max = pyip.inputInt(f'    Set upper zoom index [{x_min+1}, {a.size}]: ',
-                        min=x_min+1, max=a.size)
-            else:
-                x_upp = pyip.inputInt(f'    Set upper data bound [{x_upp_min}, {a.size}]: ',
-                        min=x_upp_min, max=a.size)
-                break
-    else:
-        if not is_int(x_upp, x_low+num_x_min, a.size):
-            illegal_value('x_upp', x_upp, 'selectArrayBounds')
-            return None
-    print(f'lower bound = {x_low} (inclusive)\nupper bound = {x_upp} (exclusive)]')
-    bounds = [x_low, x_upp]
-    #quickPlot(range(bounds[0], bounds[1]), a[bounds[0]:bounds[1]], title=title)
-    quickPlot((range(a.size), a), ([bounds[0], bounds[0]], [a.min(), a.max()], 'r-'),
-            ([bounds[1], bounds[1]], [a.min(), a.max()], 'r-'), title=title)
-    if pyip.inputYesNo('Accept these bounds ([y]/n)?: ', blank=True) == 'no':
-        bounds = selectArrayBounds(a, title=title)
-    return bounds
-
-def selectImageBounds(a, axis, low=None, upp=None, num_min=None,
-        title='select array bounds'):
-    """Interactively select the lower and upper data bounds for a 2D numpy array.
-    """
-    if not isinstance(a, np.ndarray) or a.ndim != 2:
-        logging.error('Illegal array type or dimension in selectImageBounds')
-        return None
-    if axis < 0 or axis >= a.ndim:
-        illegal_value('axis', axis, 'selectImageBounds')
-        return None
-    if num_min == None:
-        num_min = 1
-    else:
-        if num_min < 2 or num_min > a.shape[axis]:
-            logging.warning('Illegal input for num_min in selectImageBounds, input ignored')
-            num_min = 1
-    if low == None:
-        min_ = 0
-        max_ = a.shape[axis]
-        low_max = a.shape[axis]-num_min
-        while True:
-            if axis:
-                quickImshow(a[:,min_:max_], title=title, aspect='auto',
-                        extent=[min_,max_,a.shape[0],0])
-            else:
-                quickImshow(a[min_:max_,:], title=title, aspect='auto',
-                        extent=[0,a.shape[1], max_,min_])
-            zoom_flag = pyip.inputInt('Set lower data bound (0) or zoom in (1)?: ',
-                    min=0, max=1)
-            if zoom_flag:
-                min_ = pyip.inputInt(f'    Set lower zoom index [0, {low_max}]: ',
-                        min=0, max=low_max)
-                max_ = pyip.inputInt(f'    Set upper zoom index [{min_+1}, {low_max+1}]: ',
-                        min=min_+1, max=low_max+1)
-            else:
-                low = pyip.inputInt(f'    Set lower data bound [0, {low_max}]: ',
-                        min=0, max=low_max)
-                break
-    else:
-        if not is_int(low, 0, a.shape[axis]-num_min):
-            illegal_value('low', low, 'selectImageBounds')
-            return None
-    if upp == None:
-        min_ = low+num_min
-        max_ = a.shape[axis]
-        upp_min = min_
-        while True:
-            if axis:
-                quickImshow(a[:,min_:max_], title=title, aspect='auto',
-                        extent=[min_,max_,a.shape[0],0])
-            else:
-                quickImshow(a[min_:max_,:], title=title, aspect='auto',
-                        extent=[0,a.shape[1], max_,min_])
-            zoom_flag = pyip.inputInt('Set upper data bound (0) or zoom in (1)?: ',
-                    min=0, max=1)
-            if zoom_flag:
-                min_ = pyip.inputInt(f'    Set upper zoom index [{upp_min}, {a.shape[axis]-1}]: ',
-                        min=upp_min, max=a.shape[axis]-1)
-                max_ = pyip.inputInt(f'    Set upper zoom index [{min_+1}, {a.shape[axis]}]: ',
-                        min=min_+1, max=a.shape[axis])
-            else:
-                upp = pyip.inputInt(f'    Set upper data bound [{upp_min}, {a.shape[axis]}]: ',
-                        min=upp_min, max=a.shape[axis])
-                break
-    else:
-        if not is_int(upp, low+num_min, a.shape[axis]):
-            illegal_value('upp', upp, 'selectImageBounds')
-            return None
-    print(f'lower bound = {low} (inclusive)\nupper bound = {upp} (exclusive)')
-    bounds = [low, upp]
-    a_tmp = a
-    if axis:
-        a_tmp[:,bounds[0]] = a.max()
-        a_tmp[:,bounds[1]] = a.max()
-    else:
-        a_tmp[bounds[0],:] = a.max()
-        a_tmp[bounds[1],:] = a.max()
-    quickImshow(a_tmp, title=title)
-    if pyip.inputYesNo('Accept these bounds ([y]/n)?: ', blank=True) == 'no':
-        bounds = selectImageBounds(a, title=title)
-    return bounds
-
-def fitStep(x=None, y=None, model='step', form='arctan'):
-    if not isinstance(y, np.ndarray) or y.ndim != 1:
-        logging.error('Illegal array type or dimension for y in fitStep')
-        return
-    if isinstance(x, type(None)):
-        x = np.array(range(y.size))
-    elif not isinstance(x, np.ndarray) or x.ndim != 1 or x.size != y.size:
-        logging.error('Illegal array type or dimension for x in fitStep')
-        return
-    if not isinstance(model, str) or not model in ('step', 'rectangle'):
-        illegal_value('model', model, 'fitStepModel')
-        return
-    if not isinstance(form, str) or not form in ('linear', 'atan', 'arctan', 'erf', 'logistic'):
-        illegal_value('form', form, 'fitStepModel')
-        return
-
-    if model == 'step':
-        mod = StepModel(form=form)
-    else:
-        mod = RectangleModel(form=form)
-    pars = mod.guess(y, x=x)
-    out  = mod.fit(y, pars, x=x)
-    #print(out.fit_report())
-    #quickPlot((x,y),(x,out.best_fit))
-    return out.best_values
-
-class Config:
-    """Base class for processing a config file or dictionary.
-    """
-    def __init__(self, config_file=None, config_dict=None):
-        self.config = {}
-        self.load_flag = False
-        self.suffix = None
-
-        # Load config file 
-        if config_file is not None and config_dict is not None:
-            logging.warning('Ignoring config_dict (both config_file and config_dict are specified)')
-        if config_file is not None:
-           self.loadFile(config_file)
-        elif config_dict is not None:
-           self.loadDict(config_dict)
-
-    def loadFile(self, config_file):
-        """Load a config file.
-        """
-        if self.load_flag:
-            logging.warning('Overwriting any previously loaded config file')
-        self.config = {}
-
-        # Ensure config file exists
-        if not os.path.isfile(config_file):
-            logging.error(f'Unable to load {config_file}')
-            return
-
-        # Load config file
-        self.suffix = os.path.splitext(config_file)[1]
-        if self.suffix == '.yml' or self.suffix == '.yaml':
-            with open(config_file, 'r') as f:
-                self.config = yaml.safe_load(f)
-        elif self.suffix == '.txt':
-            with open(config_file, 'r') as f:
-                lines = f.read().splitlines()
-            self.config = {item[0].strip():literal_eval(item[1].strip()) for item in
-                    [line.split('#')[0].split('=') for line in lines if '=' in line.split('#')[0]]}
-        else:
-            logging.error(f'Illegal config file extension: {self.suffix}')
-
-        # Make sure config file was correctly loaded
-        if isinstance(self.config, dict):
-            self.load_flag = True
-        else:
-            logging.error(f'Unable to load dictionary from config file: {config_file}')
-            self.config = {}
-
-    def loadDict(self, config_dict):
-        """Takes a dictionary and places it into self.config.
-        """
-        exit('loadDict not tested yet, what format do we follow: txt or yaml?')
-        if self.load_flag:
-            logging.warning('Overwriting the previously loaded config file')
-
-        if isinstance(config_dict, dict):
-            self.config = config_dict
-            self.load_flag = True
-        else:
-            logging.error(f'Illegal dictionary config object: {config_dict}')
-            self.config = {}
-
-    def saveFile(self, config_file):
-        """Save the config file (as a yaml file only right now).
-        """
-        suffix = os.path.splitext(config_file)[1]
-        if suffix != '.yml' and suffix != '.yaml':
-            logging.error(f'Illegal config file extension: {suffix}')
-
-        # Check if config file exists
-        if os.path.isfile(config_file):
-            logging.info(f'Updating {config_file}')
-        else:
-            logging.info(f'Saving {config_file}')
-
-        # Save config file
-        with open(config_file, 'w') as f:
-            yaml.dump(self.config, f)
-
-    def validate(self, pars_required, pars_missing=None):
-        """Returns False if any required first level keys are missing.
-        """
-        if not self.load_flag:
-            logging.error('Load a config file prior to calling Config.validate')
-        pars = [p for p in pars_required if p not in self.config]
-        if isinstance(pars_missing, list):
-            pars_missing.extend(pars)
-        elif pars_missing is not None:
-            illegal_value('pars_missing', pars_missing, 'Config.validate')
-        if len(pars) > 0:
-            return False
-        return True
-
-#RV FIX this is for a txt file, obsolete?
-#    def update_txt(self, config_file, key, value, search_string=None, header=None):
-#        if not self.load_flag:
-#            logging.error('Load a config file prior to calling Config.update')
-#
-#        if not os.path.isfile(config_file):
-#            logging.error(f'Unable to load {config_file}')
-#            lines = []
-#        else:
-#            with open(config_file, 'r') as f:
-#                lines = f.read().splitlines()
-#        config = {item[0].strip():literal_eval(item[1].strip()) for item in
-#                [line.split('#')[0].split('=') for line in lines if '=' in line.split('#')[0]]}
-#        if not isinstance(key, str):
-#            illegal_value('key', key, 'Config.update')
-#            return config
-#        if isinstance(value, str):
-#            newline = f"{key} = '{value}'"
-#        else:
-#            newline = f'{key} = {value}'
-#        if key in config.keys():
-#            # Update key with value
-#            for index,line in enumerate(lines):
-#                if '=' in line:
-#                    item = line.split('#')[0].split('=')
-#                    if item[0].strip() == key:
-#                        lines[index] = newline
-#                        break
-#        else:
-#            # Insert new key/value pair
-#            if search_string != None:
-#                if isinstance(search_string, str):
-#                    search_string = [search_string]
-#                elif not isinstance(search_string, (tuple, list)):
-#                    illegal_value('search_string', search_string, 'Config.update')
-#                    search_string = None
-#            update_flag = False
-#            if search_string != None:
-#                indices = [[index for index,line in enumerate(lines) if item in line]
-#                        for item in search_string]
-#                for i,index in enumerate(indices):
-#                    if index:
-#                        if len(search_string) > 1 and key < search_string[i]:
-#                            lines.insert(index[0], newline)
-#                        else:
-#                            lines.insert(index[0]+1, newline)
-#                        update_flag = True
-#                        break
-#            if not update_flag:
-#                if isinstance(header, str):
-#                    lines += ['', header, newline]
-#                else:
-#                    lines += ['', newline]
-#        # Write updated config file
-#        with open(config_file, 'w') as f:
-#            for line in lines:
-#                f.write(f'{line}\n')
-#        # Update loaded config
-#        config['key'] = value
-#    
-#RV update and bring into Config if needed again
-#def search(config_file, search_string):
-#    if not os.path.isfile(config_file):
-#        logging.error(f'Unable to load {config_file}')
-#        return False
-#    with open(config_file, 'r') as f:
-#        lines = f.read()
-#        if search_string in lines:
-#            return True
-#    return False
-
-class Detector:
-    """Class for processing a detector info file or dictionary.
-    """
-    def __init__(self, detector_id):
-        self.detector = {}
-        self.load_flag = False
-        self.validate_flag = False
-
-        # Load detector file 
-        self.loadFile(detector_id)
-
-    def loadFile(self, detector_id):
-        """Load a detector file.
-        """
-        if self.load_flag:
-            logging.warning('Overwriting the previously loaded detector file')
-        self.detector = {}
-
-        # Ensure detector file exists
-        if not isinstance(detector_id, str):
-            illegal_value('detector_id', detector_id, 'Detector.loadFile')
-            return
-        detector_file = f'{detector_id}.yaml'
-        if not os.path.isfile(detector_file):
-            detector_file = self.config['detector_id']+'.yaml'
-            if not os.path.isfile(detector_file):
-                logging.error(f'Unable to load detector info file for {detector_id}')
-                return
-
-        # Load detector file
-        with open(detector_file, 'r') as f:
-            self.detector = yaml.safe_load(f)
-
-        # Make sure detector file was correctly loaded
-        if isinstance(self.detector, dict):
-            self.load_flag = True
-        else:
-            logging.error(f'Unable to load dictionary from detector file: {detector_file}')
-            self.detector = {}
-
-    def validate(self):
-        """Returns False if any config parameters is illegal or missing.
-        """
-        if not self.load_flag:
-            logging.error('Load a detector file prior to calling Detector.validate')
-
-        # Check for required first-level keys
-        pars_required = ['detector', 'lens_magnification']
-        pars_missing = [p for p in pars_required if p not in self.detector]
-        if len(pars_missing) > 0:
-            logging.error(f'Missing item(s) in detector file: {", ".join(pars_missing)}')
-            return False
-
-        is_valid = True
-
-        # Check detector pixel config keys
-        pixels = self.detector['detector'].get('pixels')
-        if not pixels:
-            pars_missing.append('detector:pixels')
-        else:
-            rows = pixels.get('rows')
-            if not rows:
-                pars_missing.append('detector:pixels:rows')
-            columns = pixels.get('columns')
-            if not columns:
-                pars_missing.append('detector:pixels:columns')
-            size = pixels.get('size')
-            if not size:
-                pars_missing.append('detector:pixels:size')
-
-        if not len(pars_missing):
-            self.validate_flag = True
-        else:
-            is_valid = False
-
-        return is_valid
-
-    def getPixelSize(self):
-        """Returns the detector pixel size.
-        """
-        if not self.validate_flag:
-            logging.error('Validate detector file info prior to calling Detector.getPixelSize')
-
-        lens_magnification = self.detector.get('lens_magnification')
-        if not isinstance(lens_magnification, (int,float)) or lens_magnification <= 0.:
-            illegal_value('lens_magnification', lens_magnification, 'detector file')
-            return 0
-        pixel_size = self.detector['detector'].get('pixels').get('size')
-        if isinstance(pixel_size, (int,float)):
-            if pixel_size <= 0.:
-                illegal_value('pixel_size', pixel_size, 'detector file')
-                return 0
-            pixel_size /= lens_magnification
-        elif isinstance(pixel_size, list):
-            if ((len(pixel_size) > 2) or
-                    (len(pixel_size) == 2 and pixel_size[0] != pixel_size[1])):
-                illegal_value('pixel size', pixel_size, 'detector file')
-                return 0
-            elif not is_num(pixel_size[0], 0.):
-                illegal_value('pixel size', pixel_size, 'detector file')
-                return 0
-            else:
-                pixel_size = pixel_size[0]/lens_magnification
-        else:
-            illegal_value('pixel size', pixel_size, 'detector file')
-            return 0
-
-        return pixel_size
-
-    def getDimensions(self):
-        """Returns the detector pixel dimensions.
-        """
-        if not self.validate_flag:
-            logging.error('Validate detector file info prior to calling Detector.getDimensions')
-
-        pixels = self.detector['detector'].get('pixels')
-        num_rows = pixels.get('rows')
-        if not is_int(num_rows, 1):
-            illegal_value('rows', num_rows, 'detector file')
-            return (0, 0)
-        num_columns = pixels.get('columns')
-        if not is_int(num_columns, 1):
-            illegal_value('columns', num_columns, 'detector file')
-            return (0, 0)
-
-        return num_rows, num_columns

--- a/tomo.py	Thu Mar 24 16:57:45 2022 +0000
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,2007 +0,0 @@
-#!/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
-import pyinputplus as pyip
-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 = []
-
-        # 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))
-
-        # Set output config file name
-        if self.config_out is None:
-            if self.config is None:
-                self.config_out = 'config.yaml'
-            else:
-                self.config_out = config_file
-
-        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 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 selectImageRanges(self, available_stacks=None):
-        """Find and check all required image files.
-        """
-        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 = num_thetas
-            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 = [0, self.tbf.shape[0]]
-        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' : img_x_bounds}
-            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], img_x_bounds[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],:] = 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], img_x_bounds[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],:] = 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], img_x_bounds[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 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 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(self.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')
-#%%============================================================================

--- a/tomo_setup.py	Thu Mar 24 16:57:45 2022 +0000
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,169 +0,0 @@
-#!/usr/bin/env python3
-
-import logging
-
-import os
-import sys
-import re
-import yaml
-import argparse
-import numpy as np
-
-from tomo import Tomo
-import msnc_tools as msnc
-
-def __main__():
-
-    # Parse command line arguments
-    parser = argparse.ArgumentParser(
-            description='Setup tomography reconstruction')
-    parser.add_argument('-i', '--inputfiles',
-            default='inputfiles.txt',
-            help='Input file collections')
-    parser.add_argument('-c', '--config',
-            help='Input config')
-    parser.add_argument('--theta_range',
-            help='Theta range (lower bound, upper bound, number of angles)')
-    parser.add_argument('--dark',
-            help='Dark field')
-    parser.add_argument('--bright',
-            help='Bright field')
-    parser.add_argument('--tomo',
-            help='First tomography image')
-    parser.add_argument('--detectorbounds',
-            help='Detector bounds')
-    parser.add_argument('--output_config',
-            help='Output config')
-    parser.add_argument('--output_data',
-            help='Preprocessed tomography data')
-    parser.add_argument('-l', '--log', 
-            type=argparse.FileType('w'),
-            default=sys.stdout,
-            help='Log file')
-    parser.add_argument('tomo_ranges', metavar='N', type=int, nargs='+')
-    args = parser.parse_args()
-
-    # Set basic log configuration
-    logging_format = '%(asctime)s : %(levelname)s - %(module)s : %(funcName)s - %(message)s'
-    log_level = 'INFO'
-    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.info(f'config = {args.config}\n')
-    logging.info(f'theta_range = {args.theta_range.split()}\n')
-    logging.info(f'dark = {args.dark}\n')
-    logging.info(f'bright = {args.bright}\n')
-    logging.info(f'tomo = {args.tomo}\n')
-    logging.info(f'detectorbounds = {args.detectorbounds}\n')
-    logging.info(f'output_config = {args.output_config}\n')
-    logging.info(f'output_data = {args.output_data}\n')
-    logging.info(f'log = {args.log}')
-    logging.info(f'is log stdout? {args.log == sys.stdout}')
-    logging.info(f'tomoranges = {args.tomo_ranges}\n\n')
-
-    # Read input files and collect data files info
-    logging.debug('data files:\n')
-    datasets = []
-    with open(args.inputfiles) as cf:
-        for line in cf:
-            if not line.strip() or line.startswith('#'):
-                continue
-            fields = [x.strip() for x in line.split('\t')]
-            filepath = fields[0]
-            element_identifier = fields[1] if len(fields) > 1 else fields[0].split('/')[-1]
-            datasets.append({'element_identifier' : fields[1], 'filepath' : filepath})
-    logging.debug(f'\ndatasets: {datasets}\n')
-
-    # Read and sort data files
-    collections = []
-    for dataset in datasets:
-        element_identifier = [x.strip() for x in dataset['element_identifier'].split('_')]
-        if len(element_identifier) > 1:
-            name = element_identifier[0]
-        else:
-            name = 'other'
-        filepath = dataset['filepath']
-        if not len(collections):
-            collections = [{'name' : name, 'filepaths' : [filepath]}]
-        else:
-            collection = [c for c in collections if c['name'] == name]
-            if len(collection):
-                collection[0]['filepaths'].append(filepath)
-            else:
-                collection = {'name' : name, 'filepaths' : [filepath]}
-                collections.append(collection)
-    logging.debug(f'\ncollections:\n{collections}\n\n')
-    if len(args.tomo_ranges) != 2*len(collections):
-        raise ValueError('Inconsistent tomo ranges size.')
-
-    # Instantiate Tomo object
-    tomo = Tomo(config_file=args.config, config_out=args.output_config, log_level=log_level,
-            log_stream=args.log, galaxy_flag=True)
-    if not tomo.is_valid:
-        raise ValueError('Invalid config file provided.')
-    logging.debug(f'config:\n{tomo.config}\n\n')
-
-    # Set theta inputs
-    theta_range = args.theta_range.split()
-    config_theta_range = tomo.config.get('theta_range')
-    if config_theta_range is None:
-        config_tomo.config['theta_range'] = {'start' : float(theta_range[0]),
-            'end' : float(theta_range[1]), 'num' : int(theta_range[2])}
-    else:
-        config_theta_range['start'] = float(theta_range[0])
-        config_theta_range['end'] = float(theta_range[1])
-        config_theta_range['num'] = int(theta_range[2])
-    tomo.cf.saveFile(args.output_config)
-
-    # Find dark field files
-    dark_field = tomo.config['dark_field']
-    tdf_files = [c['filepaths'] for c in collections if c['name'] == 'tdf']
-    if len(tdf_files) != 1 or len(tdf_files[0]) < 1:
-        logging.warning('Unable to obtain dark field files')
-        assert(dark_field['data_path'] is None)
-        assert(dark_field['img_start'] == -1)
-        assert(not dark_field['num'])
-        tdf_files = [None]
-        num_collections = 0
-    else:
-        dark_field['img_offset'] = args.tomo_ranges[0]
-        dark_field['num'] = args.tomo_ranges[1]
-        num_collections = 1
-
-    # Find bright field files
-    bright_field = tomo.config['bright_field']
-    bright_field['img_offset'] = args.tomo_ranges[2*num_collections]
-    bright_field['num'] = args.tomo_ranges[2*num_collections+1]
-    tbf_files = [c['filepaths'] for c in collections if c['name'] == 'tbf']
-    if len(tbf_files) != 1 or len(tbf_files[0]) < 1:
-        exit('Unable to obtain bright field files')
-    num_collections += 1
-
-    # Find tomography files
-    stack_info = tomo.config['stack_info']
-    if stack_info['num'] != len(collections) - num_collections:
-        raise ValueError('Inconsistent number of tomography data image sets')
-    tomo_stack_files = []
-    for stack in stack_info['stacks']:
-        stack['img_offset'] = args.tomo_ranges[2*num_collections]
-        stack['num'] = args.tomo_ranges[2*num_collections+1]
-        tomo_files = [c['filepaths'] for c in collections if c['name'] == f'set{stack["index"]}']
-        if len(tomo_files) != 1 or len(tomo_files[0]) < 1:
-            exit(f'Unable to obtain tomography images for set {stack["index"]}\n')
-        tomo_stack_files.append(tomo_files[0])
-        num_collections += 1
-
-    # Preprocess the image files
-    tomo.genTomoStacks(tdf_files[0], tbf_files[0], tomo_stack_files, args.dark, args.bright,
-        args.tomo, args.detectorbounds, args.output_data)
-    if not tomo.is_valid:
-        IOError('Unable to load all required image files.')
-
-#RV make start_theta, end_theta and num_theta inputs
-
-if __name__ == "__main__":
-    __main__()
-