tomo: tomo.py comparison

comparison tomo.py @ 3:f9c52762c32c draft

"planemo upload for repository https://github.com/rolfverberg/galaxytools commit 7dce44d576e4149f31bdc2ee4dce0bb6962badb6"

author	rv43
date	Tue, 05 Apr 2022 18:23:54 +0000
parents	b8977c98800b
children	845270a96464

comparison

equal deleted inserted replaced

-:b8977c98800b
+:f9c52762c32c
 class set_numexpr_threads:
 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)
 class ConfigTomo(msnc.Config):
 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:
+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.
 """
 # Load the bright field images
 # 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:
+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.
 """
 preprocess = self.config.get('preprocess')
 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)
+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[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[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:
 self.cf.config['preprocess'] = {'img_x_bounds' : img_x_bounds}
 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
 title = None
 msnc.quickImshow(self.tbf, title='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:
+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[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[1]-1, img_x_bounds[1]-1], [x_sum_min, x_sum_max], 'r-'),
 title='sum over theta and y')
 print(f'lower bound = {img_x_bounds[0]} (inclusive)\n'+
 f'upper bound = {img_x_bounds[1]} (exclusive)]')
 use_bounds =  pyip.inputYesNo('Accept these bounds ([y]/n)?: ', blank=True)
 if use_bounds == 'no':
 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)
+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[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[1]-1, img_x_bounds[1]-1], [x_sum_min, x_sum_max], 'r-'),
 title='sum over theta and y')
 print(f'lower bound = {img_x_bounds[0]} (inclusive)\n'+
 f'upper bound = {img_x_bounds[1]} (exclusive)]')
 use_bounds =  pyip.inputYesNo('Accept these bounds ([y]/n)?: ', blank=True)
 if use_bounds == 'no':
 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
+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_low, x_low], [x_sum_min, x_sum_max], 'r-'),
-([x_upp, x_upp], [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],
 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)
+del x_sum
 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')
 # Update config and save to file
 stack['preprocessed'] = True
 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:
 del tdf
 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
 assert(0 <= center < tomo_plane_T.shape[0])
 center_offset = center-tomo_plane_T.shape[0]/2
 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)
+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,
+if self.galaxy_flag:
-save_fig=self.save_plots, save_only=self.save_plots_only, cmap='coolwarm')
+msnc.quickImshow(edges[0,:,:], title, name=name, save_fig=True, save_only=True,
-msnc.quickImshow(edges[0,:,:], f'{base_name} gray', path=self.output_folder,
+cmap='gray', vmin=vmin, vmax=vmax)
-save_fig=self.save_plots, save_only=self.save_plots_only, cmap='gray',
+else:
-vmin=vmin, vmax=vmax)
+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
 # 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)
+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)
+eff_pixel_size, cross_sectional_dim, False, num_core)
-if not self.test_mode:
+if self.galaxy_flag:
-base_name=f'edges row{row} center_offset_vo{center_offset_vo:.2f}'
+assert(isinstance(recon_pngname, str))
-self._plotEdgesOnePlane(recon_plane, base_name)
+title = os.path.basename(recon_pngname)
-use_phase_corr = 'no'
+self._plotEdgesOnePlane(recon_plane, title, name=recon_pngname)
-if not self.test_mode:
+del sinogram_T
-use_phase_corr = pyip.inputYesNo('Try finding center using phase correlation '+
+del recon_plane
-'(y/[n])? ', blank=True)
+return float(center_offset_vo)
-if use_phase_corr == 'yes':
+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.
 while error > tol:
 prev = tomo_center
 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)
+eff_pixel_size, cross_sectional_dim, False, num_core)
-base_name=f'edges row{row} center_offset{center_offset:.2f}'
+title = f'edges row{row} center_offset{center_offset:.2f}'
-self._plotEdgesOnePlane(recon_plane, base_name)
+self._plotEdgesOnePlane(recon_plane, title)
-accept_center = 'yes'
+if (pyip.inputYesNo('Accept a center location ([y]) or continue '+
-if not self.test_mode:
+'search (n)? ', blank=True) != 'no'):
-accept_center = pyip.inputYesNo('Accept a center location ([y]) or continue '+
-'search (n)? ', blank=True)
-if accept_center != '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:
 center_offset_low = pyip.inputInt('\nEnter lower bound for center offset '+
 f'[{-int(center)}, {int(center)}]: ', min=-int(center), max=int(center))
 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)
+thetas_deg, eff_pixel_size, cross_sectional_dim, False, num_core)
-base_name=f'edges row{row} center_offset{center_offset}'
+title = f'edges row{row} center_offset{center_offset}'
-self._plotEdgesOnePlane(recon_plane, base_name)
+self._plotEdgesOnePlane(recon_plane, title)
 if pyip.inputInt('\nContinue (0) or end the search (1): ', min=0, max=1):
 break
 del sinogram_T
 del recon_plane
 '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=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):
 """Find all available image files.
 """
 # Safe updated config
 if self.is_valid:
 self.cf.saveFile(self.config_out)
 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):
 """Preprocess tomography images.
 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):
+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')
 center_stack_index = None
 if find_center and 'center_stack_index' in find_center:
 if self.test_mode:
 assert(find_center.get('completed', False) == False)
 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
+if (pyip.inputYesNo('Do you want to use this again ([y]/n)? ',
-find_center.get('completed', False) == True):
+blank=True) != 'no' and find_center.get('completed', False) == True):
 return
 # Load the required preprocessed stack
 # preprocessed stack order: row,theta,column
 num_tomo_stacks = self.config['stack_info']['num']
 center_stack = self.tomo_stacks[0]
 if not center_stack.size:
 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:
 center_stack_index = pyip.inputInt('\nEnter tomography stack index to get '
 f'rotation axis centers {available_stacks}: ')
 if find_center is None:
 self.config['find_center'] = {'center_stack_index' : center_stack_index}
 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']
 theta_start = theta_range['start']
 theta_end = theta_range['end']
 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
+accept = 'yes'
-height = center_stack.shape[0]*eff_pixel_size
+if not self.test_mode and not self.galaxy_flag:
-if not self.test_mode and pyip.inputYesNo('\nDo you want to reconstruct the full '+
+accept = 'no'
-f'sample height ({height:.3f} mm) (y/n)? ') == 'no':
+print('\nSelect bounds for image reconstruction')
-height = pyip.inputNum('\nEnter the desired reconstructed sample height '+
+if msnc.is_index_range(row_bounds, 0, center_stack.shape[0]):
-f'in mm [0, {height:.3f}]: ', min=0, max=height)
+a_tmp = np.copy(center_stack[:,0,:])
-n1 = int(0.5*(center_stack.shape[0]-height/eff_pixel_size))
+a_tmp_max = a_tmp.max()
-else:
+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
 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:
+if not self.test_mode and not self.galaxy_flag:
-msnc.quickImshow(center_stack[:,0,:], title=f'center stack theta={theta_start}',
+tmp = center_stack[:,0,:]
-path=self.output_folder, save_fig=self.save_plots,
+tmp_max = tmp.max()
-save_only=self.save_plots_only)
+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
 logging.debug(f'cross_sectional_dim = {cross_sectional_dim}')
 logging.info('Determine center offset at sample row boundaries')
 # 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:
 msnc.quickImshow(center_stack[:,0,:], title=f'theta={theta_start}', aspect='auto')
 use_row = pyip.inputYesNo('\nCurrent row index for lower center = '
 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)
+eff_pixel_size, cross_sectional_dim, num_core=num_core,
-logging.info(f'Lower center offset = {center_offset}')
+recon_pngname=recon_low_pngname)
 # Update config and save to file
 find_center['row_bounds'] = [n1, n2]
 find_center['lower_row'] = row
 find_center['lower_center_offset'] = center_offset
 lower_row = row
 # 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:
 msnc.quickImshow(center_stack[:,0,:], title=f'theta={theta_start}', aspect='auto')
 use_row = pyip.inputYesNo('\nCurrent row index for upper center = '
 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)
+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
 # Update config and save to file
 find_center['upper_row'] = row
 del plane1, plane2, plane1_shifted
 # 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']
 lower_row = find_center.get('lower_row')
 if lower_row is None:
 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 not self.galaxy_flag:
-if stack.get('reconstructed', False):
+available = False
-self.tomo_recon_stacks[i], available = self._loadTomo('recon stack', index)
+if stack.get('reconstructed', False):
-if self.tomo_recon_stacks[i].size or available:
+self.tomo_recon_stacks[i], available = self._loadTomo('recon stack', index)
-if self.tomo_stacks[i].size:
+if self.tomo_recon_stacks[i].size or available:
-self.tomo_stacks[i] = np.array([])
+if self.tomo_stacks[i].size:
-assert(stack.get('preprocessed', False) == True)
+self.tomo_stacks[i] = np.array([])
-assert(stack.get('reconstructed', False) == True)
+assert(stack.get('preprocessed', False) == True)
+assert(stack.get('reconstructed', False) == True)
+continue
+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
-else:
+assert(0 <= lower_row < upper_row < self.tomo_stacks[i].shape[0])
-stack['reconstructed'] = False
+center_offsets = [lower_center_offset-lower_row*center_slope,
-if not self.tomo_stacks[i].size:
+upper_center_offset+(self.tomo_stacks[i].shape[0]-1-upper_row)*center_slope]
-self.tomo_stacks[i], available = self._loadTomo('red stack', index,
+t0 = time()
-required=True)
+self.tomo_recon_stacks[i]= self._reconstructOneTomoStack(self.tomo_stacks[i],
-if not self.tomo_stacks[i].size:
+thetas, center_offsets=center_offsets, sigma=0.1, num_core=num_core,
-logging.error(f'Unable to load tomography stack {index} for reconstruction')
+algorithm='gridrec', run_secondary_sirt=True, secondary_iter=25)
-stack[i]['preprocessed'] = False
+logging.info(f'Reconstruction of stack {index} took {time()-t0:.2f} seconds!')
-load_error = True
+if not self.test_mode and not self.galaxy_flag:
-continue
+row_slice = int(self.tomo_stacks[i].shape[0]/2)
-assert(0 <= lower_row < upper_row < self.tomo_stacks[i].shape[0])
+title = f'{basetitle} {index} slice{row_slice}'
-center_offsets = [lower_center_offset-lower_row*center_slope,
+msnc.quickImshow(self.tomo_recon_stacks[i][row_slice,:,:], title=title,
-upper_center_offset+(self.tomo_stacks[i].shape[0]-1-upper_row)*center_slope]
+path=self.output_folder, save_fig=self.save_plots,
-t0 = time()
+save_only=self.save_plots_only)
-self.tomo_recon_stacks[i]= self._reconstructOneTomoStack(self.tomo_stacks[i],
+msnc.quickPlot(self.tomo_recon_stacks[i]
-thetas, center_offsets=center_offsets, sigma=0.1, num_core=self.num_core,
+[row_slice,int(self.tomo_recon_stacks[i].shape[2]/2),:],
-algorithm='gridrec', run_secondary_sirt=True, secondary_iter=25)
+title=f'{title} cut{int(self.tomo_recon_stacks[i].shape[2]/2)}',
-logging.info(f'Reconstruction of stack {index} took {time()-t0:.2f} seconds!')
+path=self.output_folder, save_fig=self.save_plots,
-if not self.test_mode:
+save_only=self.save_plots_only)
-row_slice = int(self.tomo_stacks[i].shape[0]/2)
+self._saveTomo('recon stack', self.tomo_recon_stacks[i], index)
-title = f'{basetitle} {index} slice{row_slice}'
+#            else:
-msnc.quickImshow(self.tomo_recon_stacks[i][row_slice,:,:], title=title,
+#                np.savetxt(self.output_folder+f'recon_stack_{index}.txt',
-path=self.output_folder, save_fig=self.save_plots,
+#                        self.tomo_recon_stacks[i][row_slice,:,:], fmt='%.6e')
-save_only=self.save_plots_only)
+self.tomo_stacks[i] = np.array([])
-msnc.quickPlot(self.tomo_recon_stacks[i]
-[row_slice,int(self.tomo_recon_stacks[i].shape[2]/2),:],
+# Update config and save to file
-title=f'{title} cut{int(self.tomo_recon_stacks[i].shape[2]/2)}',
+stack['reconstructed'] = True
-path=self.output_folder, save_fig=self.save_plots,
+combine_stacks = self.config.get('combine_stacks')
-save_only=self.save_plots_only)
+if combine_stacks and index in combine_stacks.get('stacks', []):
-self._saveTomo('recon stack', self.tomo_recon_stacks[i], index)
+combine_stacks['stacks'].remove(index)
-#                else:
+self.cf.saveFile(self.config_out)
-#                    np.savetxt(self.output_folder+f'recon_stack_{index}.txt',
-#                            self.tomo_recon_stacks[i][row_slice,:,:], fmt='%.6e')
+# Save reconstructed tomography stack to file
-self.tomo_stacks[i] = np.array([])
+if self.galaxy_flag:
+t0 = time()
-# Update config and save to file
+logging.info(f'Saving reconstructed tomography stack to {output_name} ...')
-stack['reconstructed'] = True
+save_stacks = {f'set_{stack["index"]}':tomo_stack
-combine_stacks = self.config.get('combine_stacks')
+for stack,tomo_stack in zip(stacks,self.tomo_recon_stacks)}
-if combine_stacks and index in combine_stacks.get('stacks', []):
+np.savez(output_name, **save_stacks)
-combine_stacks['stacks'].remove(index)
+logging.info(f'... done in {time()-t0:.2f} seconds!')
-self.cf.saveFile(self.config_out)
 def combineTomoStacks(self):
 """Combine the reconstructed tomography stacks.
 """
 # stack order: stack,row(z),x,y
 parser.add_argument('-t', '--test_mode',
 action='store_true',
 default=False,
 help='Test mode flag')
 parser.add_argument('--num_core',
+type=int,
 default=-1,
 help='Number of cores')
 args = parser.parse_args()
 if args.config is None:

Mercurial > repos > rv43 > tomo

comparison tomo.py @ 3:f9c52762c32c draft