Mercurial > repos > rv43 > tomo
changeset 48:059819ea1f0e draft
"planemo upload for repository https://github.com/rolfverberg/galaxytools commit b97c6a8f181dea6d2f9cfa2069b86d30dcc47b4d"
| author | rv43 |
|---|---|
| date | Wed, 27 Apr 2022 17:28:26 +0000 |
| parents | c693117d01e9 |
| children | 26f99fdd8d61 |
| files | tomo.py tomo_combine.py tomo_combine.xml tomo_find_center.xml tomo_reconstruct.xml |
| diffstat | 5 files changed, 288 insertions(+), 81 deletions(-) [+] |
line wrap: on
line diff
--- a/tomo.py Mon Apr 25 20:02:15 2022 +0000 +++ b/tomo.py Wed Apr 27 17:28:26 2022 +0000 @@ -1564,7 +1564,7 @@ return dark_files, bright_files, tomo_stack_files - def loadTomoStacks(self, input_name): + def loadTomoStacks(self, input_name, recon_flag=False): """Load tomography stacks (only for Galaxy). """ assert(self.galaxy_flag) @@ -1574,10 +1574,16 @@ stacks = stack_info.get('stacks') assert(len(self.tomo_stacks) == stack_info['num']) with np.load(input_name) as f: - for i,stack in enumerate(stacks): - self.tomo_stacks[i] = f[f'set_{stack["index"]}'] - logging.info(f'loaded stack {i}: index = {stack["index"]}, shape = '+ - f'{self.tomo_stacks[i].shape}') + if recon_flag: + for i,stack in enumerate(stacks): + self.tomo_recon_stacks[i] = f[f'set_{stack["index"]}'] + logging.info(f'loaded stack {i}: index = {stack["index"]}, shape = '+ + f'{self.tomo_recon_stacks[i].shape}') + else: + for i,stack in enumerate(stacks): + self.tomo_stacks[i] = f[f'set_{stack["index"]}'] + logging.info(f'loaded stack {i}: index = {stack["index"]}, shape = '+ + f'{self.tomo_stacks[i].shape}') logging.info(f'... done in {time()-t0:.2f} seconds!') def genTomoStacks(self, galaxy_param=None, num_core=None): @@ -2192,23 +2198,41 @@ tomosum = np.concatenate([tomosum]+ [np.sum(self.tomo_recon_stacks[i][row_bounds[0]:row_bounds[1],:,:], axis=(1,2)) for i in range(1, num_tomo_stacks-1)]) - print(f'tomosum.shape = {tomosum.shape}') if num_tomo_stacks > 1: tomosum = np.concatenate([tomosum, np.sum(self.tomo_recon_stacks[num_tomo_stacks-1][row_bounds[0]:,:,:], axis=(1,2))]) - print(f'tomosum.shape = {tomosum.shape}') msnc.quickPlot(tomosum, title='recon stack sum xy', path='center_slice_pngs', save_fig=True, save_only=True) - def combineTomoStacks(self): + def combineTomoStacks(self, galaxy_param=None): """Combine the reconstructed tomography stacks. """ # stack order: stack,row(z),x,y + if self.galaxy_flag: + assert(galaxy_param) + if not os.path.exists('combine_pngs'): + os.mkdir('combine_pngs') logging.debug('Combine reconstructed tomography stacks') - # Load any unloaded reconstructed stacks stack_info = self.config['stack_info'] stacks = stack_info['stacks'] + assert(len(self.tomo_recon_stacks) == stack_info['num']) + assert(len(self.tomo_recon_stacks) == len(stacks)) + if self.galaxy_flag: + assert(isinstance(galaxy_param, dict)) + # Get image bounds + x_bounds = galaxy_param['x_bounds'] + assert(isinstance(x_bounds, list) and len(x_bounds) == 2) + y_bounds = galaxy_param['y_bounds'] + assert(isinstance(y_bounds, list) and len(y_bounds) == 2) + z_bounds = galaxy_param['z_bounds'] + assert(isinstance(z_bounds, list) and len(z_bounds) == 2) + else: + if galaxy_param: + logging.warning('Ignoring galaxy_param in reconstructTomoStacks (only for Galaxy)') + galaxy_param = None + + # Load any unloaded reconstructed stacks for i,stack in enumerate(stacks): available = False if not self.tomo_recon_stacks[i].size and stack.get('reconstructed', False): @@ -2236,69 +2260,97 @@ [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)? ', + if self.galaxy_flag: + if x_bounds[0] == -1: + x_bounds[0] = 0 + if x_bounds[1] == -1: + x_bounds[1] = tomosum.size + if not msnc.is_index_range(x_bounds, 0, tomosum.size): + msnc.illegal_value('x_bounds', x_bounds, 'galaxy input') + tomosum_min = tomosum.min() + tomosum_max = tomosum.max() + msnc.quickPlot((range(tomosum.size), tomosum), + ([x_bounds[0], x_bounds[0]], [tomosum_min, tomosum_max], 'r-'), + ([x_bounds[1]-1, x_bounds[1]-1], [tomosum_min, tomosum_max], 'r-'), + title=f'recon stack sum yz', path='combine_pngs', save_fig=True, save_only=True) + else: + if combine_stacks and 'x_bounds' in combine_stacks: + x_bounds = combine_stacks['x_bounds'] + if not msnc.is_index_range(x_bounds, 0, tomosum.size): + 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': - 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]] + if pyip.inputYesNo( + 'Do you want to change the image x-bounds ([y]/n)? ', + blank=True) == 'no': + x_bounds = [0, tomosum.size] + else: + x_bounds = msnc.selectArrayBounds(tomosum, title='recon stack sum yz') else: - x_bounds = msnc.selectArrayBounds(tomosum, title='recon stack sum yz') - if False and self.test_mode: - np.savetxt(f'{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') + 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, tomosum.size] + else: + x_bounds = msnc.selectArrayBounds(tomosum, title='recon stack sum yz') + if False and self.test_mode: + np.savetxt(f'{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') # Selecting y bounds (in xz-plane) 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(y_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)? ', + if self.galaxy_flag: + if y_bounds[0] == -1: + y_bounds[0] = 0 + if y_bounds[1] == -1: + y_bounds[1] = tomosum.size + if not msnc.is_index_range(y_bounds, 0, tomosum.size): + msnc.illegal_value('y_bounds', y_bounds, 'galaxy input') + tomosum_min = tomosum.min() + tomosum_max = tomosum.max() + msnc.quickPlot((range(tomosum.size), tomosum), + ([y_bounds[0], y_bounds[0]], [tomosum_min, tomosum_max], 'r-'), + ([y_bounds[1]-1, y_bounds[1]-1], [tomosum_min, tomosum_max], 'r-'), + title=f'recon stack sum xz', path='combine_pngs', save_fig=True, save_only=True) + else: + if combine_stacks and 'y_bounds' in combine_stacks: + y_bounds = combine_stacks['y_bounds'] + if not msnc.is_index_range(y_bounds, 0, tomosum.size): + 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': - 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]] + if pyip.inputYesNo( + 'Do you want to change the image y-bounds ([y]/n)? ', + blank=True) == 'no': + y_bounds = [0, tomosum.size] + else: + y_bounds = msnc.selectArrayBounds(tomosum, title='recon stack sum yz') else: - y_bounds = msnc.selectArrayBounds(tomosum, title='recon stack sum xz') - if False and self.test_mode: - np.savetxt(f'{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') + 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, tomosum.size] + else: + y_bounds = msnc.selectArrayBounds(tomosum, title='recon stack sum xz') + if False and self.test_mode: + np.savetxt(f'{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 = stack_info['num'] + num_tomo_stacks = len(stacks) if num_tomo_stacks == 1: low_bound = row_bounds[0] else: @@ -2333,7 +2385,7 @@ np.savetxt(f'{self.output_folder}/recon_combined_sum_xy.txt', tomosum, fmt='%.6e') np.savetxt(f'{self.output_folder}/recon_combined.txt', - tomo_recon_combined[int(tomo_recon_combined.shape[0]/2),:,:], fmt='%.6e') + tomo_recon_combined[int(tomosum.size/2),:,:], fmt='%.6e') # Update config and save to file if combine_stacks: @@ -2347,37 +2399,63 @@ return # Selecting z bounds (in xy-plane) - 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]] + if self.galaxy_flag: + if z_bounds[0] == -1: + z_bounds[0] = 0 + if z_bounds[1] == -1: + z_bounds[1] = tomosum.size + if not msnc.is_index_range(z_bounds, 0, tomosum.size): + msnc.illegal_value('z_bounds', z_bounds, 'galaxy input') + tomosum_min = tomosum.min() + tomosum_max = tomosum.max() + msnc.quickPlot((range(tomosum.size), tomosum), + ([z_bounds[0], z_bounds[0]], [tomosum_min, tomosum_max], 'r-'), + ([z_bounds[1]-1, z_bounds[1]-1], [tomosum_min, tomosum_max], 'r-'), + title=f'recon stack sum xy', path='combine_pngs', save_fig=True, save_only=True) 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]: + 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, tomosum.size] + else: + z_bounds = msnc.selectArrayBounds(tomosum, title='recon combined sum xy') + if self.save_plots_only: + msnc.clearFig('recon combined sum xy') + if z_bounds[0] != 0 or z_bounds[1] != tomosum.size: tomo_recon_combined = tomo_recon_combined[z_bounds[0]:z_bounds[1],:,:] - if self.save_plots_only: - msnc.clearFig('recon combined sum xy') # Plot center slices + if self.galaxy_flag: + path = 'combine_pngs' + save_fig = True + save_only = True + else: + 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[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) + path=path, save_fig=save_fig, save_only=save_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) + path=path, save_fig=save_fig, save_only=save_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) + path=path, save_fig=save_fig, save_only=save_only) - # Save combined reconstructed tomo stacks - base_name = 'recon combined' - for stack in stacks: - base_name += f' {stack["index"]}' - self._saveTomo(base_name, tomo_recon_combined) + # Save combined reconstructed tomography stack to file + if self.galaxy_flag: + t0 = time() + output_name = galaxy_param['output_name'] + logging.info(f'Saving combined reconstructed tomography stack to {output_name} ...') + np.save(output_name, tomo_recon_combined) + logging.info(f'... done in {time()-t0:.2f} seconds!') + else: + base_name = 'recon combined' + for stack in stacks: + base_name += f' {stack["index"]}' + self._saveTomo(base_name, tomo_recon_combined) # Update config and save to file if combine_stacks:
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/tomo_combine.py Wed Apr 27 17:28:26 2022 +0000 @@ -0,0 +1,82 @@ +#!/usr/bin/env python3 + +import logging + +import sys +import argparse +import tracemalloc + +from tomo import Tomo + +def __main__(): + + # Parse command line arguments + parser = argparse.ArgumentParser( + description='Combine reconstructed tomography stacks') + parser.add_argument('-i', '--input_stacks', + help='Reconstructed image file stacks') + parser.add_argument('-c', '--config', + help='Input config') + parser.add_argument('--x_bounds', + required=True, nargs=2, type=int, help='Reconstructed range in x direction') + parser.add_argument('--y_bounds', + required=True, nargs=2, type=int, help='Reconstructed range in y direction') + parser.add_argument('--z_bounds', + required=True, nargs=2, type=int, help='Reconstructed range in z direction') + parser.add_argument('--output_config', + help='Output config') + parser.add_argument('--output_data', + help='Combined tomography stacks') + parser.add_argument('-l', '--log', + type=argparse.FileType('w'), + default=sys.stdout, + help='Log file') + args = parser.parse_args() + + # Starting memory monitoring + tracemalloc.start() + + # 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.debug(f'input_stacks = {args.input_stacks}') + logging.debug(f'config = {args.config}') + logging.debug(f'x_bounds = {args.x_bounds} {type(args.x_bounds)}') + logging.debug(f'y_bounds = {args.y_bounds} {type(args.y_bounds)}') + logging.debug(f'z_bounds = {args.z_bounds} {type(args.z_bounds)}') + logging.debug(f'output_config = {args.output_config}') + logging.debug(f'output_data = {args.output_data}') + logging.debug(f'log = {args.log}') + logging.debug(f'is log stdout? {args.log is sys.stdout}') + + # 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}') + + # Load reconstructed image files + tomo.loadTomoStacks(args.input_stacks, recon_flag=True) + + # Combined reconstructed tomography stacks + galaxy_param = {'x_bounds' : args.x_bounds, 'y_bounds' : args.y_bounds, + 'z_bounds' : args.z_bounds, 'output_name' : args.output_data} + logging.info(f'galaxy_param = {galaxy_param}') + tomo.combineTomoStacks(galaxy_param) + + # Displaying memory usage + logging.info(f'Memory usage: {tracemalloc.get_traced_memory()}') + + # stopping memory monitoring + tracemalloc.stop() + +if __name__ == "__main__": + __main__() +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/tomo_combine.xml Wed Apr 27 17:28:26 2022 +0000 @@ -0,0 +1,47 @@ +<tool id="tomo_combine" name="Tomo Combine Reconstructed Stacks" version="0.1.0" python_template_version="3.9"> + <description>Combine reconstructed tomography stacks</description> + <macros> + <import>tomo_macros.xml</import> + </macros> + <expand macro="requirements" /> + <command detect_errors="exit_code"><![CDATA[ + mkdir combine_pngs; + $__tool_directory__/tomo_combine.py + -i '$recon_stacks' + -c '$config' + --x_bounds $x_bounds.low $x_bounds.upp + --y_bounds $y_bounds.low $y_bounds.upp + --z_bounds $z_bounds.low $z_bounds.upp + --output_data 'output_data.npy' + --output_config 'output_config.yaml' + -l '$log' + ]]></command> + <inputs> + <expand macro="common_inputs"/> + <param name="recon_stacks" type='data' format='npz' optional='false' label="Reconstructed stacks"/> + <section name="x_bounds" title="Reconstructed range in x direction"> + <param name="low" type="integer" value = "-1" label="Lower x-bound"/> + <param name="upp" type="integer" value = "-1" label="Upper x-bound"/> + </section> + <section name="y_bounds" title="Reconstructed range in y direction"> + <param name="low" type="integer" value = "-1" label="Lower y-bound"/> + <param name="upp" type="integer" value = "-1" label="Upper y-bound"/> + </section> + <section name="z_bounds" title="Reconstructed range in z direction"> + <param name="low" type="integer" value = "-1" label="Lower z-bound"/> + <param name="upp" type="integer" value = "-1" label="Upper z-bound"/> + </section> + </inputs> + <outputs> + <expand macro="common_outputs"/> + <data name="output_data" format="npy" label="Combined tomography stacks" from_work_dir="output_data.npy"/> + <collection name="combine_pngs" type="list" label="Recontructed slices midway in each combined dimension"> + <discover_datasets pattern="__name_and_ext__" directory="combine_pngs"/> + </collection> + <data name="output_config" format="yaml" label="Output config reconstruction" from_work_dir="output_config.yaml"/> + </outputs> + <help><![CDATA[ + Combine reconstructed tomography images. + ]]></help> + <expand macro="citations"/> +</tool>
--- a/tomo_find_center.xml Mon Apr 25 20:02:15 2022 +0000 +++ b/tomo_find_center.xml Wed Apr 27 17:28:26 2022 +0000 @@ -55,7 +55,7 @@ </inputs> <outputs> <expand macro="common_outputs"/> - <collection name="recon_center_set" type="list" label="Recontructed slices at various centers"> + <collection name="find_center_pngs" type="list" label="Recontructed slices at various centers"> <discover_datasets pattern="__name_and_ext__" directory="find_center_pngs"/> </collection> <data name="output_config" format="yaml" label="Output config find center" from_work_dir="output_config.yaml"/>
--- a/tomo_reconstruct.xml Mon Apr 25 20:02:15 2022 +0000 +++ b/tomo_reconstruct.xml Wed Apr 27 17:28:26 2022 +0000 @@ -25,7 +25,7 @@ <outputs> <expand macro="common_outputs"/> <data name="output_data" format="npz" label="Reconstructed tomography stacks" from_work_dir="output_data.npz"/> - <collection name="recon_slices" type="list" label="Recontructed slices midway in each dimension"> + <collection name="center_slice_pngs" type="list" label="Recontructed slices midway in each dimension"> <discover_datasets pattern="__name_and_ext__" directory="center_slice_pngs"/> </collection> <data name="output_config" format="yaml" label="Output config reconstruction" from_work_dir="output_config.yaml"/>