# HG changeset patch
# User rv43
# Date 1649183034 0
# Node ID f9c52762c32c7fd69b6b454d6cd76e7cde4956ec
# Parent b8977c98800be5907fb2a814373825066d29e8f0
"planemo upload for repository https://github.com/rolfverberg/galaxytools commit 7dce44d576e4149f31bdc2ee4dce0bb6962badb6"
diff -r b8977c98800b -r f9c52762c32c __pycache__/msnc_tools.cpython-39.pyc
Binary file __pycache__/msnc_tools.cpython-39.pyc has changed
diff -r b8977c98800b -r f9c52762c32c __pycache__/tomo.cpython-39.pyc
Binary file __pycache__/tomo.cpython-39.pyc has changed
diff -r b8977c98800b -r f9c52762c32c msnc_tools.py
--- a/msnc_tools.py Thu Mar 31 20:48:17 2022 +0000
+++ b/msnc_tools.py Tue Apr 05 18:23:54 2022 +0000
@@ -297,7 +297,7 @@
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):
+ clear=True, extent=None, show_grid=False, grid_color='w', grid_linewidth=1, **kwargs):
if title is not None and not isinstance(title, str):
illegal_value('title', title, 'quickImshow')
return
@@ -327,24 +327,32 @@
path = name
else:
path = f'{path}/{name}'
+ if extent is None:
+ extent = (0, a.shape[1], a.shape[0], 0)
if clear:
plt.close(fig=title)
if save_only:
plt.figure(title)
- plt.imshow(a, **kwargs)
+ plt.imshow(a, extent=extent, **kwargs)
+ if show_grid:
+ ax = plt.gca()
+ ax.grid(color=grid_color, linewidth=grid_linewidth)
plt.savefig(path)
plt.close(fig=title)
#plt.imsave(f'{title}.png', a, **kwargs)
else:
plt.ion()
plt.figure(title)
- plt.imshow(a, **kwargs)
+ plt.imshow(a, extent=extent, **kwargs)
+ if show_grid:
+ ax = plt.gca()
+ ax.grid(color=grid_color, linewidth=grid_linewidth)
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):
+ clear=True, show_grid=False, **kwargs):
if title is not None and not isinstance(title, str):
illegal_value('title', title, 'quickPlot')
return
@@ -383,6 +391,9 @@
plt.plot(*y, **kwargs)
else:
plt.plot(*args, **kwargs)
+ if show_grid:
+ ax = plt.gca()
+ ax.grid(color='k')#, linewidth=1)
plt.savefig(path)
plt.close(fig=title)
else:
@@ -393,6 +404,9 @@
plt.plot(*y, **kwargs)
else:
plt.plot(*args, **kwargs)
+ if show_grid:
+ ax = plt.gca()
+ ax.grid(color='k')#, linewidth=1)
if save_fig:
plt.savefig(path)
plt.pause(1)
@@ -536,15 +550,17 @@
illegal_value('upp', upp, 'selectImageBounds')
return None
bounds = [low, upp]
- a_tmp = a
+ a_tmp = np.copy(a)
+ a_tmp_max = a.max()
if axis:
- a_tmp[:,bounds[0]] = a.max()
- a_tmp[:,bounds[1]] = a.max()
+ a_tmp[:,bounds[0]] = a_tmp_max
+ a_tmp[:,bounds[1]-1] = a_tmp_max
else:
- a_tmp[bounds[0],:] = a.max()
- a_tmp[bounds[1],:] = a.max()
+ a_tmp[bounds[0],:] = a_tmp_max
+ a_tmp[bounds[1]-1,:] = a_tmp_max
print(f'lower bound = {low} (inclusive)\nupper bound = {upp} (exclusive)')
quickImshow(a_tmp, title=title)
+ del a_tmp
if pyip.inputYesNo('Accept these bounds ([y]/n)?: ', blank=True) == 'no':
bounds = selectImageBounds(a, axis, low=low_save, upp=upp_save, num_min=num_min_save,
title=title)
diff -r b8977c98800b -r f9c52762c32c run_tomo_find_center
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/run_tomo_find_center Tue Apr 05 18:23:54 2022 +0000
@@ -0,0 +1,3 @@
+#!/bin/bash
+
+python tomo_find_center.py -i red_stacks.npz -c 'config_l_center.yaml' --row_bounds '620 950' --center_rows '670 890' --output_config 'output.yaml' --recon_center_low 'recon_center_low.png' --recon_center_upp 'recon_center_upp.png'
diff -r b8977c98800b -r f9c52762c32c run_tomo_reconstruct
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/run_tomo_reconstruct Tue Apr 05 18:23:54 2022 +0000
@@ -0,0 +1,3 @@
+#!/bin/bash
+
+python tomo_reconstruct.py -i red_stacks.npz -c 'config_l_recon.yaml' --output_config 'output.yaml' --output_data 'output.npz'
diff -r b8977c98800b -r f9c52762c32c run_tomo_setup
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/run_tomo_setup Tue Apr 05 18:23:54 2022 +0000
@@ -0,0 +1,3 @@
+#!/bin/bash
+
+python tomo_setup.py -i inputfiles_l.txt -c 'config_l_setup.yaml' --theta_range '0.0 180.0 354' --output_config 'output.yaml' --output_data 'output.npz' --dark 'dark.png' --bright 'bright.png' --tomo 'tomo.png' --detectorbounds 'detectorbounds.png' 1 20 1 354
diff -r b8977c98800b -r f9c52762c32c tomo.py
--- a/tomo.py Thu Mar 31 20:48:17 2022 +0000
+++ b/tomo.py Tue Apr 05 18:23:54 2022 +0000
@@ -34,13 +34,16 @@
def __init__(self, num_core):
cpu_count = mp.cpu_count()
+ logging.debug(f'start: num_core={num_core} cpu_count={cpu_count}')
if num_core is None or num_core < 1 or num_core > cpu_count:
self.num_core = cpu_count
else:
self.num_core = num_core
+ logging.debug(f'self.num_core={self.num_core}')
def __enter__(self):
self.num_core_org = ne.set_num_threads(self.num_core)
+ logging.debug(f'self.num_core={self.num_core}')
def __exit__(self, exc_type, exc_value, traceback):
ne.set_num_threads(self.num_core_org)
@@ -647,12 +650,12 @@
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:
+ if self.galaxy_flag:
+ msnc.quickImshow(self.tdf, title='dark field', name=dark_field_pngname,
+ save_fig=True, save_only=True)
+ elif not self.test_mode:
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.
@@ -681,12 +684,12 @@
self.tbf -= self.tdf
else:
logging.warning('Dark field unavailable')
- if not self.test_mode and not self.galaxy_flag:
+ if self.galaxy_flag:
+ msnc.quickImshow(self.tbf, title='bright field', name=bright_field_pngname,
+ save_fig=True, save_only=True)
+ elif not self.test_mode:
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.
@@ -743,14 +746,16 @@
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)
+ x_sum_min = x_sum.min()
+ x_sum_max = x_sum.max()
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)
+ save_fig=True, save_only=True, show_grid=True)
msnc.quickPlot((range(x_sum.size), x_sum),
- ([img_x_bounds[0], img_x_bounds[0]], [x_sum.min(), x_sum.max()], 'r-'),
- ([img_x_bounds[1]-1, img_x_bounds[1]-1], [x_sum.min(), x_sum.max()], 'r-'),
+ ([img_x_bounds[0], img_x_bounds[0]], [x_sum_min, x_sum_max], 'r-'),
+ ([img_x_bounds[1]-1, img_x_bounds[1]-1], [x_sum_min, x_sum_max], 'r-'),
title='sum over theta and y', name=detectorbounds_pngname,
- save_fig=True, save_only=True)
+ save_fig=True, save_only=True, show_grid=True)
# Update config and save to file
if preprocess is None:
@@ -758,6 +763,7 @@
else:
preprocess['img_x_bounds'] = img_x_bounds
self.cf.saveFile(self.config_out)
+ del x_sum
return
# For one tomography stack only: load the first image
@@ -777,9 +783,17 @@
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:
+ tmp = np.copy(tomo_stack[0,:,:])
+ tmp_max = tmp.max()
+ tmp[img_x_bounds[0],:] = tmp_max
+ tmp[img_x_bounds[1]-1,:] = tmp_max
+ msnc.quickImshow(tmp, title=title)
+ del tmp
+ x_sum_min = x_sum.min()
+ x_sum_max = x_sum.max()
msnc.quickPlot((range(x_sum.size), x_sum),
- ([img_x_bounds[0], img_x_bounds[0]], [x_sum.min(), x_sum.max()], 'r-'),
- ([img_x_bounds[1]-1, img_x_bounds[1]-1], [x_sum.min(), x_sum.max()], 'r-'),
+ ([img_x_bounds[0], img_x_bounds[0]], [x_sum_min, x_sum_max], 'r-'),
+ ([img_x_bounds[1]-1, img_x_bounds[1]-1], [x_sum_min, x_sum_max], 'r-'),
title='sum over theta and y')
print(f'lower bound = {img_x_bounds[0]} (inclusive)\n'+
f'upper bound = {img_x_bounds[1]} (exclusive)]')
@@ -799,11 +813,20 @@
save_fig=self.save_plots, save_only=True)
else:
x_sum = np.sum(self.tbf, 1)
+ x_sum_min = x_sum.min()
+ x_sum_max = x_sum.max()
use_bounds = 'no'
if img_x_bounds[0] is not None and img_x_bounds[1] is not None:
+ tmp = np.copy(self.tbf)
+ tmp_max = tmp.max()
+ tmp[img_x_bounds[0],:] = tmp_max
+ tmp[img_x_bounds[1]-1,:] = tmp_max
+ title = 'Bright field'
+ msnc.quickImshow(tmp, title=title)
+ del tmp
msnc.quickPlot((range(x_sum.size), x_sum),
- ([img_x_bounds[0], img_x_bounds[0]], [x_sum.min(), x_sum.max()], 'r-'),
- ([img_x_bounds[1]-1, img_x_bounds[1]-1], [x_sum.min(), x_sum.max()], 'r-'),
+ ([img_x_bounds[0], img_x_bounds[0]], [x_sum_min, x_sum_max], 'r-'),
+ ([img_x_bounds[1]-1, img_x_bounds[1]-1], [x_sum_min, x_sum_max], 'r-'),
title='sum over theta and y')
print(f'lower bound = {img_x_bounds[0]} (inclusive)\n'+
f'upper bound = {img_x_bounds[1]} (exclusive)]')
@@ -820,9 +843,16 @@
x_upp = int(x_upp+(x_upp-x_low)/10)
if x_upp >= x_sum.size:
x_upp = x_sum.size
+ tmp = np.copy(self.tbf)
+ tmp_max = tmp.max()
+ tmp[x_low,:] = tmp_max
+ tmp[x_upp-1,:] = tmp_max
+ title = 'Bright field'
+ msnc.quickImshow(tmp, title=title)
+ del tmp
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-'),
+ ([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)
@@ -837,6 +867,7 @@
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)
+ del x_sum
logging.debug(f'img_x_bounds: {img_x_bounds}')
if self.save_plots_only:
@@ -1069,13 +1100,10 @@
stack.pop('reconstructed', 'reconstructed not found')
find_center = self.config.get('find_center')
if find_center:
+ find_center.pop('completed', 'completed not found')
if self.test_mode:
- find_center.pop('completed', 'completed not found')
find_center.pop('lower_center_offset', 'lower_center_offset not found')
find_center.pop('upper_center_offset', 'upper_center_offset not found')
- else:
- if find_center.get('center_stack_index', -1) == index:
- self.config.pop('find_center')
self.cf.saveFile(self.config_out)
if self.tdf.size:
@@ -1083,7 +1111,7 @@
del tbf
def _reconstructOnePlane(self, tomo_plane_T, center, thetas_deg, eff_pixel_size,
- cross_sectional_dim, plot_sinogram=True):
+ cross_sectional_dim, plot_sinogram=True, num_core=1):
"""Invert the sinogram for a single tomography plane.
"""
# tomo_plane_T index order: column,theta
@@ -1117,24 +1145,28 @@
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)
+ recon_clean = tomopy.misc.corr.remove_ring(recon_clean, rwidth=17, ncore=num_core)
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):
+ def _plotEdgesOnePlane(self, recon_plane, title, name=None, 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)
+ if self.galaxy_flag:
+ msnc.quickImshow(edges[0,:,:], title, name=name, save_fig=True, save_only=True,
+ cmap='gray', vmin=vmin, vmax=vmax)
+ else:
+ msnc.quickImshow(edges[0,:,:], f'{title} coolwarm', path=self.output_folder,
+ save_fig=self.save_plots, save_only=self.save_plots_only, cmap='coolwarm')
+ msnc.quickImshow(edges[0,:,:], f'{title} 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):
+ tol=0.1, num_core=1, recon_pngname=None):
"""Find center for a single tomography plane.
"""
# sinogram index order: theta,column
@@ -1143,20 +1175,31 @@
center = sinogram.shape[1]/2
# try automatic center finding routines for initial value
- tomo_center = tomopy.find_center_vo(sinogram)
+ tomo_center = tomopy.find_center_vo(sinogram, ncore=num_core)
center_offset_vo = tomo_center-center
- if not self.test_mode:
+ if self.test_mode or self.galaxy_flag:
+ logging.info(f'Center at row {row} using Nghia Vo’s method = {center_offset_vo:.2f}')
+ if self.test_mode:
+ del sinogram_T
+ return float(center_offset_vo)
+ else:
print(f'Center at row {row} using Nghia Vo’s method = {center_offset_vo:.2f}')
+ if recon_pngname:
+ logging.warning('Ignoring recon_pngname in _findCenterOnePlane (only for Galaxy)')
recon_plane = self._reconstructOnePlane(sinogram_T, tomo_center, thetas_deg,
- eff_pixel_size, cross_sectional_dim, False)
- if not self.test_mode:
- base_name=f'edges row{row} center_offset_vo{center_offset_vo:.2f}'
- self._plotEdgesOnePlane(recon_plane, base_name)
- use_phase_corr = 'no'
- if not self.test_mode:
- use_phase_corr = pyip.inputYesNo('Try finding center using phase correlation '+
- '(y/[n])? ', blank=True)
- if use_phase_corr == 'yes':
+ eff_pixel_size, cross_sectional_dim, False, num_core)
+ if self.galaxy_flag:
+ assert(isinstance(recon_pngname, str))
+ title = os.path.basename(recon_pngname)
+ self._plotEdgesOnePlane(recon_plane, title, name=recon_pngname)
+ del sinogram_T
+ del recon_plane
+ return float(center_offset_vo)
+ else:
+ title = f'edges row{row} center_offset_vo{center_offset_vo:.2f}'
+ self._plotEdgesOnePlane(recon_plane, title)
+ 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.
@@ -1168,24 +1211,18 @@
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)
- accept_center = 'yes'
- if not self.test_mode:
- accept_center = pyip.inputYesNo('Accept a center location ([y]) or continue '+
- 'search (n)? ', blank=True)
- if accept_center != 'no':
+ eff_pixel_size, cross_sectional_dim, False, num_core)
+ title = f'edges row{row} center_offset{center_offset:.2f}'
+ self._plotEdgesOnePlane(recon_plane, title)
+ if (pyip.inputYesNo('Accept a center location ([y]) or continue '+
+ 'search (n)? ', blank=True) != 'no'):
del sinogram_T
del recon_plane
- if self.test_mode:
+ 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
- else:
- 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:
@@ -1204,9 +1241,9 @@
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)
+ thetas_deg, eff_pixel_size, cross_sectional_dim, False, num_core)
+ title = f'edges row{row} center_offset{center_offset}'
+ self._plotEdgesOnePlane(recon_plane, title)
if pyip.inputInt('\nContinue (0) or end the search (1): ', min=0, max=1):
break
@@ -1266,7 +1303,8 @@
init_recon=tomo_recon_stack, options=options, sinogram_order=True,
algorithm=tomopy.astra, ncore=num_core)
if True:
- tomopy.misc.corr.remove_ring(tomo_recon_stack, rwidth=rwidth, out=tomo_recon_stack)
+ tomopy.misc.corr.remove_ring(tomo_recon_stack, rwidth=rwidth, out=tomo_recon_stack,
+ ncore=num_core)
return tomo_recon_stack
def findImageFiles(self):
@@ -1359,6 +1397,22 @@
return dark_files, bright_files, tomo_stack_files
+ def loadTomoStacks(self, input_name):
+ """Load tomography stacks (only for Galaxy).
+ """
+ assert(self.galaxy_flag)
+ t0 = time()
+ logging.info(f'Loading preprocessed tomography stack from {input_name} ...')
+ stack_info = self.config['stack_info']
+ 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}')
+ logging.info(f'... done in {time()-t0:.2f} seconds!')
+
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):
@@ -1459,13 +1513,34 @@
stack['ref_height'] = stack['ref_height']+pixel_size
self.cf.saveFile(self.config_out)
- def findCenters(self):
+ def findCenters(self, row_bounds=None, center_rows=None, recon_low_pngname=None,
+ recon_upp_pngname=None, num_core=None):
"""Find rotation axis centers for the tomography stacks.
"""
+ if num_core is None:
+ num_core = self.num_core
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('Available stacks: {available_stacks}')
+ if self.galaxy_flag:
+ assert(isinstance(row_bounds, list) and len(row_bounds) == 2)
+ assert(isinstance(center_rows, list) and len(center_rows) == 2)
+ assert(isinstance(recon_low_pngname, str))
+ assert(isinstance(recon_upp_pngname, str))
+ else:
+ if row_bounds:
+ logging.warning('Ignoring row_bounds in findCenters (only for Galaxy)')
+ row_bounds = None
+ if center_rows:
+ logging.warning('Ignoring center_rows in findCenters (only for Galaxy)')
+ center_rows = None
+ if recon_low_pngname:
+ logging.warning('Ignoring recon_low_pngname in findCenters (only for Galaxy)')
+ recon_low_pngname = None
+ if recon_upp_pngname:
+ logging.warning('Ignoring recon_upp_pngname in findCenters (only for Galaxy)')
+ recon_upp_pngname = None
# Check for valid available center stack index
find_center = self.config.get('find_center')
@@ -1481,8 +1556,8 @@
else:
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' and
- find_center.get('completed', False) == True):
+ if (pyip.inputYesNo('Do you want to use this again ([y]/n)? ',
+ blank=True) != 'no' and find_center.get('completed', False) == True):
return
# Load the required preprocessed stack
@@ -1503,6 +1578,9 @@
stacks[0]['preprocessed'] = False
raise OSError('Unable to load the required preprocessed tomography stack')
assert(stacks[0].get('preprocessed', False) == True)
+ elif self.galaxy_flag:
+ logging.error('CHECK/FIX FOR GALAXY')
+ #center_stack_index = stacks[int(num_tomo_stacks/2)]['index']
else:
while True:
if not center_stack_index:
@@ -1528,6 +1606,10 @@
find_center = self.config['find_center']
else:
find_center['center_stack_index'] = center_stack_index
+ if not self.galaxy_flag:
+ row_bounds = find_center.get('row_bounds', None)
+ center_rows = [find_center.get('lower_row', None),
+ find_center.get('upper_row', None)]
# Set thetas (in degrees)
theta_range = self.config['theta_range']
@@ -1545,26 +1627,54 @@
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 not self.test_mode and pyip.inputYesNo('\nDo you want to reconstruct the full '+
- f'sample height ({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))
+ accept = 'yes'
+ if not self.test_mode and not self.galaxy_flag:
+ accept = 'no'
+ print('\nSelect bounds for image reconstruction')
+ if msnc.is_index_range(row_bounds, 0, center_stack.shape[0]):
+ a_tmp = np.copy(center_stack[:,0,:])
+ a_tmp_max = a_tmp.max()
+ a_tmp[row_bounds[0],:] = a_tmp_max
+ a_tmp[row_bounds[1]-1,:] = a_tmp_max
+ print(f'lower bound = {row_bounds[0]} (inclusive)')
+ print(f'upper bound = {row_bounds[1]} (exclusive)')
+ msnc.quickImshow(a_tmp, title=f'center stack theta={theta_start}',
+ aspect='auto')
+ del a_tmp
+ accept = pyip.inputYesNo('Accept these bounds ([y]/n)?: ', blank=True)
+ if accept == 'no':
+ (n1, n2) = msnc.selectImageBounds(center_stack[:,0,:], 0,
+ title=f'center stack theta={theta_start}')
+ else:
+ n1 = row_bounds[0]
+ n2 = row_bounds[1]
else:
+ logging.error('CHECK/FIX FOR GALAXY')
+ tomo_ref_heights = [stack['ref_height'] for stack in stacks]
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
+ 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')
- if not self.test_mode:
- 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)
+ if not self.test_mode and not self.galaxy_flag:
+ tmp = center_stack[:,0,:]
+ tmp_max = tmp.max()
+ tmp[n1,:] = tmp_max
+ tmp[n2-1,:] = tmp_max
+ if msnc.is_index_range(center_rows, 0, tmp.shape[0]):
+ tmp[center_rows[0],:] = tmp_max
+ tmp[center_rows[1]-1,:] = tmp_max
+ extent = [0, tmp.shape[1], tmp.shape[0], 0]
+ msnc.quickImshow(tmp, title=f'center stack theta={theta_start}',
+ path=self.output_folder, extent=extent, save_fig=self.save_plots,
+ save_only=self.save_plots_only, aspect='auto')
+ del tmp
+ #extent = [0, center_stack.shape[2], n2, n1]
+ #msnc.quickImshow(center_stack[n1:n2,0,:], title=f'center stack theta={theta_start}',
+ # path=self.output_folder, extent=extent, save_fig=self.save_plots,
+ # save_only=self.save_plots_only, show_grid=True, aspect='auto')
# Get cross sectional diameter in mm
cross_sectional_dim = center_stack.shape[2]*eff_pixel_size
@@ -1576,9 +1686,9 @@
# Lower row center
use_row = 'no'
use_center = 'no'
- row = find_center.get('lower_row')
+ row = center_rows[0]
if msnc.is_int(row, n1, n2-2):
- if self.test_mode:
+ if self.test_mode or self.galaxy_flag:
assert(row is not None)
use_row = 'yes'
else:
@@ -1605,8 +1715,8 @@
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}')
+ eff_pixel_size, cross_sectional_dim, num_core=num_core,
+ recon_pngname=recon_low_pngname)
# Update config and save to file
find_center['row_bounds'] = [n1, n2]
@@ -1618,9 +1728,9 @@
# Upper row center
use_row = 'no'
use_center = 'no'
- row = find_center.get('upper_row')
+ row = center_rows[1]
if msnc.is_int(row, lower_row+1, n2-1):
- if self.test_mode:
+ if self.test_mode or self.galaxy_flag:
assert(row is not None)
use_row = 'yes'
else:
@@ -1647,7 +1757,8 @@
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)
+ eff_pixel_size, cross_sectional_dim, num_core=num_core,
+ recon_pngname=recon_upp_pngname)
logging.info(f'upper_center_offset = {center_offset}')
del center_stack
@@ -1731,14 +1842,22 @@
# Update config file
self.config = msnc.update('config.txt', 'check_centers', True, 'find_centers')
- def reconstructTomoStacks(self):
+ def reconstructTomoStacks(self, output_name=None, num_core=None):
"""Reconstruct tomography stacks.
"""
+ if num_core is None:
+ num_core = self.num_core
logging.debug('Reconstruct tomography stacks')
stacks = self.config['stack_info']['stacks']
assert(len(self.tomo_stacks) == self.config['stack_info']['num'])
assert(len(self.tomo_stacks) == len(stacks))
assert(len(self.tomo_recon_stacks) == len(stacks))
+ if self.galaxy_flag:
+ assert(isinstance(output_name, str))
+ else:
+ if output_name:
+ logging.warning('Ignoring output_name in reconstructTomoStacks '+
+ '(only used in Galaxy)')
# Get rotation axis rows and centers
find_center = self.config['find_center']
@@ -1784,56 +1903,65 @@
# 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('preprocessed', False) == True)
- 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')
- stack[i]['preprocessed'] = False
- load_error = True
+ if not self.galaxy_flag:
+ 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('preprocessed', False) == True)
+ assert(stack.get('reconstructed', False) == 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, num_core=self.num_core,
- 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([])
+ 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')
+ stack[i]['preprocessed'] = False
+ 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, num_core=num_core,
+ 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 and not self.galaxy_flag:
+ 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
- combine_stacks = self.config.get('combine_stacks')
- if combine_stacks and index in combine_stacks.get('stacks', []):
- combine_stacks['stacks'].remove(index)
- self.cf.saveFile(self.config_out)
+ # Update config and save to file
+ stack['reconstructed'] = True
+ combine_stacks = self.config.get('combine_stacks')
+ if combine_stacks and index in combine_stacks.get('stacks', []):
+ combine_stacks['stacks'].remove(index)
+ self.cf.saveFile(self.config_out)
+
+ # Save reconstructed tomography stack to file
+ if self.galaxy_flag:
+ t0 = time()
+ logging.info(f'Saving reconstructed tomography stack to {output_name} ...')
+ save_stacks = {f'set_{stack["index"]}':tomo_stack
+ for stack,tomo_stack in zip(stacks,self.tomo_recon_stacks)}
+ np.savez(output_name, **save_stacks)
+ logging.info(f'... done in {time()-t0:.2f} seconds!')
def combineTomoStacks(self):
"""Combine the reconstructed tomography stacks.
@@ -2095,6 +2223,7 @@
default=False,
help='Test mode flag')
parser.add_argument('--num_core',
+ type=int,
default=-1,
help='Number of cores')
args = parser.parse_args()
diff -r b8977c98800b -r f9c52762c32c tomo_find_center.py
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/tomo_find_center.py Tue Apr 05 18:23:54 2022 +0000
@@ -0,0 +1,74 @@
+#!/usr/bin/env python3
+
+import logging
+
+import sys
+import re
+import argparse
+
+from tomo import Tomo
+
+def __main__():
+
+ # Parse command line arguments
+ parser = argparse.ArgumentParser(
+ description='Find the center axis for a tomography reconstruction')
+ parser.add_argument('-i', '--input_stacks',
+ help='Preprocessed image file stacks')
+ parser.add_argument('-c', '--config',
+ help='Input config')
+ parser.add_argument('--row_bounds',
+ help='Reconstruction row bounds')
+ parser.add_argument('--center_rows',
+ help='Center finding rows')
+ parser.add_argument('--output_config',
+ help='Output config')
+ parser.add_argument('--recon_center_low',
+ help='Lower reconstructed slice center')
+ parser.add_argument('--recon_center_upp',
+ help='Upper reconstructed slice center')
+ parser.add_argument('-l', '--log',
+ type=argparse.FileType('w'),
+ default=sys.stdout,
+ help='Log file')
+ args = parser.parse_args()
+
+ indexRegex = re.compile(r'\d+')
+ row_bounds = [int(i) for i in indexRegex.findall(args.row_bounds)]
+ center_rows = [int(i) for i in indexRegex.findall(args.center_rows)]
+
+ # 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'row_bounds = {args.row_bounds} {row_bounds}')
+ logging.debug(f'center_rows = {args.center_rows} {center_rows}')
+ logging.debug(f'output_config = {args.output_config}')
+ logging.debug(f'recon_center_low = {args.recon_center_low}')
+ logging.debug(f'recon_center_uppm = {args.recon_center_upp}')
+ 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 preprocessed image files
+ tomo.loadTomoStacks(args.input_stacks)
+
+ # Find centers
+ tomo.findCenters(row_bounds, center_rows, args.recon_center_low, args.recon_center_upp)
+
+if __name__ == "__main__":
+ __main__()
+
diff -r b8977c98800b -r f9c52762c32c tomo_find_center.xml
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/tomo_find_center.xml Tue Apr 05 18:23:54 2022 +0000
@@ -0,0 +1,48 @@
+
+ Find the center axis for a tomography reconstruction
+
+ tomopy
+ h5py
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+@misc{githubsum_files,
+ author = {Verberg, Rolf},
+ year = {2022},
+ title = {Tomo Find Center Axis},
+}
+
+
+
diff -r b8977c98800b -r f9c52762c32c tomo_reconstruct.py
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/tomo_reconstruct.py Tue Apr 05 18:23:54 2022 +0000
@@ -0,0 +1,60 @@
+#!/usr/bin/env python3
+
+import logging
+
+import sys
+import argparse
+
+from tomo import Tomo
+
+def __main__():
+
+ # Parse command line arguments
+ parser = argparse.ArgumentParser(
+ description='Perfrom a tomography reconstruction')
+ parser.add_argument('-i', '--input_stacks',
+ help='Preprocessed image file stacks')
+ parser.add_argument('-c', '--config',
+ help='Input config')
+ parser.add_argument('--output_config',
+ help='Output config')
+ parser.add_argument('--output_data',
+ help='Reconstructed tomography data')
+ parser.add_argument('-l', '--log',
+ type=argparse.FileType('w'),
+ default=sys.stdout,
+ help='Log file')
+ 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.debug(f'input_stacks = {args.input_stacks}')
+ logging.debug(f'config = {args.config}')
+ 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 preprocessed image files
+ tomo.loadTomoStacks(args.input_stacks)
+
+ # Reconstruct tomography stacks
+ tomo.reconstructTomoStacks(args.output_data)
+
+if __name__ == "__main__":
+ __main__()
+
diff -r b8977c98800b -r f9c52762c32c tomo_reconstruct.xml
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/tomo_reconstruct.xml Tue Apr 05 18:23:54 2022 +0000
@@ -0,0 +1,36 @@
+
+ Perform a tomography reconstruction
+
+ tomopy
+ h5py
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+@misc{githubsum_files,
+ author = {Verberg, Rolf},
+ year = {2022},
+ title = {Tomo Reconstruction},
+}
+
+
+
diff -r b8977c98800b -r f9c52762c32c tomo_setup.xml
--- a/tomo_setup.xml Thu Mar 31 20:48:17 2022 +0000
+++ b/tomo_setup.xml Tue Apr 05 18:23:54 2022 +0000
@@ -60,9 +60,6 @@
author = {Verberg, Rolf},
year = {2022},
title = {Tomo Setup},
- publisher = {GitHub},
- journal = {GitHub repository},
- url = {somewhere},
}