Mercurial > repos > rv43 > tomo

comparison tomo.py @ 31:063e1e847d95 draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
"planemo upload for repository https://github.com/rolfverberg/galaxytools commit bd7bbc4facb864b77e56b8c22bcbb497d1815e95"
author rv43
date Mon, 18 Apr 2022 21:54:56 +0000
parents c4a90d5b1ed2
children b1f92b63f84c
comparison
equal deleted inserted replaced
30:c4a90d5b1ed2 31:063e1e847d95
1181 # need index order column,theta for iradon, so take transpose 1181 # need index order column,theta for iradon, so take transpose
1182 sinogram_T = sinogram.T 1182 sinogram_T = sinogram.T
1183 center = sinogram.shape[1]/2 1183 center = sinogram.shape[1]/2
1184 1184
1185 # try automatic center finding routines for initial value 1185 # try automatic center finding routines for initial value
1186 #num_core=1
1187 tomo_center = tomopy.find_center_vo(sinogram, ncore=num_core) 1186 tomo_center = tomopy.find_center_vo(sinogram, ncore=num_core)
1188 center_offset_vo = tomo_center-center 1187 center_offset_vo = tomo_center-center
1189 if self.test_mode: 1188 if self.test_mode:
1190 logging.info(f'Center at row {row} using Nghia Vo’s method = {center_offset_vo:.2f}') 1189 logging.info(f'Center at row {row} using Nghia Vo’s method = {center_offset_vo:.2f}')
1191 del sinogram_T 1190 del sinogram_T
1299 # RV should we remove stripes? 1298 # RV should we remove stripes?
1300 # https://tomopy.readthedocs.io/en/latest/api/tomopy.prep.stripe.html 1299 # https://tomopy.readthedocs.io/en/latest/api/tomopy.prep.stripe.html
1301 # RV should we remove rings? 1300 # RV should we remove rings?
1302 # https://tomopy.readthedocs.io/en/latest/api/tomopy.misc.corr.html 1301 # https://tomopy.readthedocs.io/en/latest/api/tomopy.misc.corr.html
1303 # RV: Add an option to do (extra) secondary iterations later or to do some sort of convergence test? 1302 # RV: Add an option to do (extra) secondary iterations later or to do some sort of convergence test?
1303 print(f'OK AA row_bounds = {row_bounds}')
1304 if row_bounds is None: 1304 if row_bounds is None:
1305 row_bounds = [0, tomo_stack.shape[0]] 1305 row_bounds = [0, tomo_stack.shape[0]]
1306 else: 1306 else:
1307 if not msnc.is_index_range(row_bounds, 0, tomo_stack.shape[0]): 1307 if not msnc.is_index_range(row_bounds, 0, tomo_stack.shape[0]):
1308 raise ValueError('Illegal row bounds in reconstructOneTomoStack') 1308 raise ValueError('Illegal row bounds in reconstructOneTomoStack')
1316 else: 1316 else:
1317 if center_offsets.size != row_bounds[1]-row_bounds[0]: 1317 if center_offsets.size != row_bounds[1]-row_bounds[0]:
1318 raise ValueError('center_offsets dimension mismatch in reconstructOneTomoStack') 1318 raise ValueError('center_offsets dimension mismatch in reconstructOneTomoStack')
1319 centers = center_offsets 1319 centers = center_offsets
1320 centers += tomo_stack.shape[2]/2 1320 centers += tomo_stack.shape[2]/2
1321 print(f'num_core = {num_core}')
1322 print(f'thetas = {thetas}')
1323 print(f'centers = {centers}')
1321 if True: 1324 if True:
1322 tomo_stack = tomopy.prep.stripe.remove_stripe_fw( 1325 tomo_stack = tomopy.prep.stripe.remove_stripe_fw(
1323 tomo_stack[row_bounds[0]:row_bounds[1]], sigma=sigma, ncore=num_core) 1326 tomo_stack[row_bounds[0]:row_bounds[1]], sigma=sigma, ncore=num_core)
1324 else: 1327 else:
1325 tomo_stack = tomo_stack[row_bounds[0]:row_bounds[1]] 1328 tomo_stack = tomo_stack[row_bounds[0]:row_bounds[1]]
1326 print('OK AA') 1329 print('OK BB')
1327 tomo_recon_stack = tomopy.recon(tomo_stack, thetas, centers, sinogram_order=True, 1330 tomo_recon_stack = tomopy.recon(tomo_stack, thetas, centers, sinogram_order=True,
1328 algorithm=algorithm, ncore=num_core) 1331 algorithm=algorithm, ncore=num_core)
1329 print('OK BB') 1332 print('OK CC')
1330 if run_secondary_sirt and secondary_iter > 0: 1333 if run_secondary_sirt and secondary_iter > 0:
1331 #options = {'method':'SIRT_CUDA', 'proj_type':'cuda', 'num_iter':secondary_iter} 1334 #options = {'method':'SIRT_CUDA', 'proj_type':'cuda', 'num_iter':secondary_iter}
1332 #RV: doesn't work for me: "Error: CUDA error 803: system has unsupported display driver / 1335 #RV: doesn't work for me: "Error: CUDA error 803: system has unsupported display driver /
1333 # cuda driver combination." 1336 # cuda driver combination."
1334 #options = {'method':'SIRT', 'proj_type':'linear', 'MinConstraint': 0, 'num_iter':secondary_iter} 1337 #options = {'method':'SIRT', 'proj_type':'linear', 'MinConstraint': 0, 'num_iter':secondary_iter}
1337 options = {'method':'SART', 'proj_type':'linear', 'MinConstraint': 0, 1340 options = {'method':'SART', 'proj_type':'linear', 'MinConstraint': 0,
1338 'num_iter':secondary_iter} 1341 'num_iter':secondary_iter}
1339 tomo_recon_stack = tomopy.recon(tomo_stack, thetas, centers, 1342 tomo_recon_stack = tomopy.recon(tomo_stack, thetas, centers,
1340 init_recon=tomo_recon_stack, options=options, sinogram_order=True, 1343 init_recon=tomo_recon_stack, options=options, sinogram_order=True,
1341 algorithm=tomopy.astra, ncore=num_core) 1344 algorithm=tomopy.astra, ncore=num_core)
1342 print('OK CC') 1345 print('OK DD')
1343 if True: 1346 if True:
1344 tomopy.misc.corr.remove_ring(tomo_recon_stack, rwidth=rwidth, out=tomo_recon_stack, 1347 tomopy.misc.corr.remove_ring(tomo_recon_stack, rwidth=rwidth, out=tomo_recon_stack,
1345 ncore=num_core) 1348 ncore=num_core)
1346 return tomo_recon_stack 1349 return tomo_recon_stack
1347 1350
1454 def genTomoStacks(self, galaxy_param=None, num_core=None): 1457 def genTomoStacks(self, galaxy_param=None, num_core=None):
1455 """Preprocess tomography images. 1458 """Preprocess tomography images.
1456 """ 1459 """
1457 if num_core is None: 1460 if num_core is None:
1458 num_core = self.num_core 1461 num_core = self.num_core
1462 logging.info(f'num_core = {self.num_core}')
1459 # Try loading any already preprocessed stacks (skip in Galaxy) 1463 # Try loading any already preprocessed stacks (skip in Galaxy)
1460 # preprocessed stack order for each one in stack: row,theta,column 1464 # preprocessed stack order for each one in stack: row,theta,column
1461 stack_info = self.config['stack_info'] 1465 stack_info = self.config['stack_info']
1462 stacks = stack_info['stacks'] 1466 stacks = stack_info['stacks']
1463 num_tomo_stacks = stack_info['num'] 1467 num_tomo_stacks = stack_info['num']
1537 def findCenters(self, galaxy_param=None, num_core=None): 1541 def findCenters(self, galaxy_param=None, num_core=None):
1538 """Find rotation axis centers for the tomography stacks. 1542 """Find rotation axis centers for the tomography stacks.
1539 """ 1543 """
1540 if num_core is None: 1544 if num_core is None:
1541 num_core = self.num_core 1545 num_core = self.num_core
1546 logging.info(f'num_core = {self.num_core}')
1542 logging.debug('Find centers for tomography stacks') 1547 logging.debug('Find centers for tomography stacks')
1543 stacks = self.config['stack_info']['stacks'] 1548 stacks = self.config['stack_info']['stacks']
1544 available_stacks = [stack['index'] for stack in stacks if stack.get('preprocessed', False)] 1549 available_stacks = [stack['index'] for stack in stacks if stack.get('preprocessed', False)]
1545 logging.debug('Available stacks: {available_stacks}') 1550 logging.debug('Available stacks: {available_stacks}')
1546 if self.galaxy_flag: 1551 if self.galaxy_flag:
1888 def reconstructTomoStacks(self, galaxy_param=None, num_core=None): 1893 def reconstructTomoStacks(self, galaxy_param=None, num_core=None):
1889 """Reconstruct tomography stacks. 1894 """Reconstruct tomography stacks.
1890 """ 1895 """
1891 if num_core is None: 1896 if num_core is None:
1892 num_core = self.num_core 1897 num_core = self.num_core
1898 logging.info(f'num_core = {self.num_core}')
1893 if self.galaxy_flag: 1899 if self.galaxy_flag:
1894 assert(galaxy_param) 1900 assert(galaxy_param)
1895 if not os.path.exists('center_slice_pngs'): 1901 if not os.path.exists('center_slice_pngs'):
1896 os.mkdir('center_slice_pngs') 1902 os.mkdir('center_slice_pngs')
1897 logging.debug('Reconstruct tomography stacks') 1903 logging.debug('Reconstruct tomography stacks')
1988 print(f'lower_row = {lower_row} upper_row = {upper_row} self.tomo_stacks[i].shape[0] = {self.tomo_stacks[i].shape[0]}') 1994 print(f'lower_row = {lower_row} upper_row = {upper_row} self.tomo_stacks[i].shape[0] = {self.tomo_stacks[i].shape[0]}')
1989 assert(0 <= lower_row < upper_row < self.tomo_stacks[i].shape[0]) 1995 assert(0 <= lower_row < upper_row < self.tomo_stacks[i].shape[0])
1990 center_offsets = [lower_center_offset-lower_row*center_slope, 1996 center_offsets = [lower_center_offset-lower_row*center_slope,
1991 upper_center_offset+(self.tomo_stacks[i].shape[0]-1-upper_row)*center_slope] 1997 upper_center_offset+(self.tomo_stacks[i].shape[0]-1-upper_row)*center_slope]
1992 t0 = time() 1998 t0 = time()
1993 num_core = 1
1994 print(f'OK4 {i} c') 1999 print(f'OK4 {i} c')
1995 self.tomo_recon_stacks[i]= self._reconstructOneTomoStack(self.tomo_stacks[i], 2000 self.tomo_recon_stacks[i]= self._reconstructOneTomoStack(self.tomo_stacks[i],
1996 thetas, center_offsets=center_offsets, sigma=0.1, num_core=num_core, 2001 thetas, center_offsets=center_offsets, sigma=0.1, num_core=num_core,
1997 algorithm='gridrec', run_secondary_sirt=True, secondary_iter=25) 2002 algorithm='gridrec', run_secondary_sirt=True, secondary_iter=25)
1998 print(f'OK4 {i} d') 2003 print(f'OK4 {i} d')
1999 exit('Done')
2000 logging.info(f'Reconstruction of stack {index} took {time()-t0:.2f} seconds!') 2004 logging.info(f'Reconstruction of stack {index} took {time()-t0:.2f} seconds!')
2001 if self.galaxy_flag: 2005 if self.galaxy_flag:
2002 x_slice = int(self.tomo_stacks[i].shape[0]/2) 2006 x_slice = int(self.tomo_stacks[i].shape[0]/2)
2003 title = f'{basetitle} {index} xslice{x_slice}' 2007 title = f'{basetitle} {index} xslice{x_slice}'
2004 msnc.quickImshow(self.tomo_recon_stacks[i][x_slice,:,:], title=title, 2008 msnc.quickImshow(self.tomo_recon_stacks[i][x_slice,:,:], title=title,