Mercurial > repos > rv43 > tomo

comparison tomo.py @ 34:63de912323e5 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 c81dfe083e164b51b3688e512f2a0b00b953851c"
author rv43
date Tue, 19 Apr 2022 16:40:34 +0000
parents afb612b64f26
children 20829cd7aaef
comparison
equal deleted inserted replaced
33:afb612b64f26 34:63de912323e5
1298 # RV should we remove stripes? 1298 # RV should we remove stripes?
1299 # https://tomopy.readthedocs.io/en/latest/api/tomopy.prep.stripe.html 1299 # https://tomopy.readthedocs.io/en/latest/api/tomopy.prep.stripe.html
1300 # RV should we remove rings? 1300 # RV should we remove rings?
1301 # https://tomopy.readthedocs.io/en/latest/api/tomopy.misc.corr.html 1301 # https://tomopy.readthedocs.io/en/latest/api/tomopy.misc.corr.html
1302 # 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: 1303 if row_bounds is None:
1305 row_bounds = [0, tomo_stack.shape[0]] 1304 row_bounds = [0, tomo_stack.shape[0]]
1306 else: 1305 else:
1307 if not msnc.is_index_range(row_bounds, 0, tomo_stack.shape[0]): 1306 if not msnc.is_index_range(row_bounds, 0, tomo_stack.shape[0]):
1308 raise ValueError('Illegal row bounds in reconstructOneTomoStack') 1307 raise ValueError('Illegal row bounds in reconstructOneTomoStack')
1316 else: 1315 else:
1317 if center_offsets.size != row_bounds[1]-row_bounds[0]: 1316 if center_offsets.size != row_bounds[1]-row_bounds[0]:
1318 raise ValueError('center_offsets dimension mismatch in reconstructOneTomoStack') 1317 raise ValueError('center_offsets dimension mismatch in reconstructOneTomoStack')
1319 centers = center_offsets 1318 centers = center_offsets
1320 centers += tomo_stack.shape[2]/2 1319 centers += tomo_stack.shape[2]/2
1321 print(f'num_core = {num_core}') 1320 # RV hangs here with more than 24 cores and sge_64G_4
1322 print(f'thetas = {thetas}')
1323 print(f'centers = {centers}')
1324 if True: 1321 if True:
1325 tomo_stack = tomopy.prep.stripe.remove_stripe_fw( 1322 tomo_stack = tomopy.prep.stripe.remove_stripe_fw(
1326 tomo_stack[row_bounds[0]:row_bounds[1]], sigma=sigma, ncore=num_core) 1323 tomo_stack[row_bounds[0]:row_bounds[1]], sigma=sigma, ncore=num_core)
1327 else: 1324 else:
1328 tomo_stack = tomo_stack[row_bounds[0]:row_bounds[1]] 1325 tomo_stack = tomo_stack[row_bounds[0]:row_bounds[1]]
1329 print('OK BB') 1326 logging.info('performing initial reconstruction')
1330 tomo_recon_stack = tomopy.recon(tomo_stack, thetas, centers, sinogram_order=True, 1327 tomo_recon_stack = tomopy.recon(tomo_stack, thetas, centers, sinogram_order=True,
1331 algorithm=algorithm, ncore=num_core) 1328 algorithm=algorithm, ncore=num_core)
1332 print('OK CC')
1333 if run_secondary_sirt and secondary_iter > 0: 1329 if run_secondary_sirt and secondary_iter > 0:
1330 logging.info(f'running {secondary_iter} secondary iterations')
1334 #options = {'method':'SIRT_CUDA', 'proj_type':'cuda', 'num_iter':secondary_iter} 1331 #options = {'method':'SIRT_CUDA', 'proj_type':'cuda', 'num_iter':secondary_iter}
1335 #RV: doesn't work for me: "Error: CUDA error 803: system has unsupported display driver / 1332 #RV: doesn't work for me: "Error: CUDA error 803: system has unsupported display driver /
1336 # cuda driver combination." 1333 # cuda driver combination."
1337 #options = {'method':'SIRT', 'proj_type':'linear', 'MinConstraint': 0, 'num_iter':secondary_iter} 1334 #options = {'method':'SIRT', 'proj_type':'linear', 'MinConstraint': 0, 'num_iter':secondary_iter}
1338 #SIRT did not finish while running overnight 1335 #SIRT did not finish while running overnight
1340 options = {'method':'SART', 'proj_type':'linear', 'MinConstraint': 0, 1337 options = {'method':'SART', 'proj_type':'linear', 'MinConstraint': 0,
1341 'num_iter':secondary_iter} 1338 'num_iter':secondary_iter}
1342 tomo_recon_stack = tomopy.recon(tomo_stack, thetas, centers, 1339 tomo_recon_stack = tomopy.recon(tomo_stack, thetas, centers,
1343 init_recon=tomo_recon_stack, options=options, sinogram_order=True, 1340 init_recon=tomo_recon_stack, options=options, sinogram_order=True,
1344 algorithm=tomopy.astra, ncore=num_core) 1341 algorithm=tomopy.astra, ncore=num_core)
1345 print('OK DD')
1346 if True: 1342 if True:
1347 tomopy.misc.corr.remove_ring(tomo_recon_stack, rwidth=rwidth, out=tomo_recon_stack, 1343 tomopy.misc.corr.remove_ring(tomo_recon_stack, rwidth=rwidth, out=tomo_recon_stack,
1348 ncore=num_core) 1344 ncore=num_core)
1349 return tomo_recon_stack 1345 return tomo_recon_stack
1350 1346
1893 def reconstructTomoStacks(self, galaxy_param=None, num_core=None): 1889 def reconstructTomoStacks(self, galaxy_param=None, num_core=None):
1894 """Reconstruct tomography stacks. 1890 """Reconstruct tomography stacks.
1895 """ 1891 """
1896 if num_core is None: 1892 if num_core is None:
1897 num_core = self.num_core 1893 num_core = self.num_core
1898 logging.info(f'num_core = {num_core}') 1894 logging.info(f'num_core available = {num_core}')
1895 if num_core > 24:
1896 num_core = 24
1897 logging.info(f'num_core used = {num_core}')
1899 if self.galaxy_flag: 1898 if self.galaxy_flag:
1900 assert(galaxy_param) 1899 assert(galaxy_param)
1901 if not os.path.exists('center_slice_pngs'): 1900 if not os.path.exists('center_slice_pngs'):
1902 os.mkdir('center_slice_pngs') 1901 os.mkdir('center_slice_pngs')
1903 logging.debug('Reconstruct tomography stacks') 1902 logging.debug('Reconstruct tomography stacks')
1904 stacks = self.config['stack_info']['stacks'] 1903 stacks = self.config['stack_info']['stacks']
1905 assert(len(self.tomo_stacks) == self.config['stack_info']['num']) 1904 assert(len(self.tomo_stacks) == self.config['stack_info']['num'])
1906 assert(len(self.tomo_stacks) == len(stacks)) 1905 assert(len(self.tomo_stacks) == len(stacks))
1907 assert(len(self.tomo_recon_stacks) == len(stacks)) 1906 assert(len(self.tomo_recon_stacks) == len(stacks))
1908 print('OK1')
1909 if self.galaxy_flag: 1907 if self.galaxy_flag:
1910 assert(isinstance(galaxy_param, dict)) 1908 assert(isinstance(galaxy_param, dict))
1911 # Get rotation axis centers 1909 # Get rotation axis centers
1912 center_offsets = galaxy_param['center_offsets'] 1910 center_offsets = galaxy_param['center_offsets']
1913 assert(isinstance(center_offsets, list) and len(center_offsets) == 2) 1911 assert(isinstance(center_offsets, list) and len(center_offsets) == 2)
1921 galaxy_param = None 1919 galaxy_param = None
1922 lower_center_offset = None 1920 lower_center_offset = None
1923 upper_center_offset = None 1921 upper_center_offset = None
1924 1922
1925 # Get rotation axis rows and centers 1923 # Get rotation axis rows and centers
1926 print('OK2')
1927 find_center = self.config['find_center'] 1924 find_center = self.config['find_center']
1928 lower_row = find_center.get('lower_row') 1925 lower_row = find_center.get('lower_row')
1929 if lower_row is None: 1926 if lower_row is None:
1930 logging.error('Unable to read lower_row from config') 1927 logging.error('Unable to read lower_row from config')
1931 return 1928 return
1955 num_theta_skip = self.config['preprocess'].get('num_theta_skip', 0) 1952 num_theta_skip = self.config['preprocess'].get('num_theta_skip', 0)
1956 thetas = np.radians(np.linspace(theta_start, theta_end, 1953 thetas = np.radians(np.linspace(theta_start, theta_end,
1957 int(num_theta/(num_theta_skip+1)), endpoint=False)) 1954 int(num_theta/(num_theta_skip+1)), endpoint=False))
1958 1955
1959 # Reconstruct tomo stacks 1956 # Reconstruct tomo stacks
1960 print('OK3')
1961 zoom_perc = self.config['preprocess'].get('zoom_perc', 100) 1957 zoom_perc = self.config['preprocess'].get('zoom_perc', 100)
1962 if zoom_perc == 100: 1958 if zoom_perc == 100:
1963 basetitle = 'recon stack fullres' 1959 basetitle = 'recon stack fullres'
1964 else: 1960 else:
1965 basetitle = f'recon stack {zoom_perc}p' 1961 basetitle = f'recon stack {zoom_perc}p'
1968 for i,stack in enumerate(stacks): 1964 for i,stack in enumerate(stacks):
1969 # Check if stack can be loaded 1965 # Check if stack can be loaded
1970 # reconstructed stack order for each one in stack : row/z,x,y 1966 # reconstructed stack order for each one in stack : row/z,x,y
1971 # preprocessed stack order for each one in stack: row,theta,column 1967 # preprocessed stack order for each one in stack: row,theta,column
1972 index = stack['index'] 1968 index = stack['index']
1973 print(f'OK4 {i} a')
1974 if not self.galaxy_flag: 1969 if not self.galaxy_flag:
1975 available = False 1970 available = False
1976 if stack.get('reconstructed', False): 1971 if stack.get('reconstructed', False):
1977 self.tomo_recon_stacks[i], available = self._loadTomo('recon stack', index) 1972 self.tomo_recon_stacks[i], available = self._loadTomo('recon stack', index)
1978 if self.tomo_recon_stacks[i].size or available: 1973 if self.tomo_recon_stacks[i].size or available:
1979 if self.tomo_stacks[i].size: 1974 if self.tomo_stacks[i].size:
1980 self.tomo_stacks[i] = np.array([]) 1975 self.tomo_stacks[i] = np.array([])
1981 assert(stack.get('preprocessed', False) == True) 1976 assert(stack.get('preprocessed', False) == True)
1982 assert(stack.get('reconstructed', False) == True) 1977 assert(stack.get('reconstructed', False) == True)
1983 continue 1978 continue
1984 print(f'OK4 {i} b size = {self.tomo_stacks[i].size}')
1985 stack['reconstructed'] = False 1979 stack['reconstructed'] = False
1986 if not self.tomo_stacks[i].size: 1980 if not self.tomo_stacks[i].size:
1987 self.tomo_stacks[i], available = self._loadTomo('red stack', index, 1981 self.tomo_stacks[i], available = self._loadTomo('red stack', index,
1988 required=True) 1982 required=True)
1989 if not self.tomo_stacks[i].size: 1983 if not self.tomo_stacks[i].size:
1990 logging.error(f'Unable to load tomography stack {index} for reconstruction') 1984 logging.error(f'Unable to load tomography stack {index} for reconstruction')
1991 stack[i]['preprocessed'] = False 1985 stack[i]['preprocessed'] = False
1992 load_error = True 1986 load_error = True
1993 continue 1987 continue
1994 print(f'lower_row = {lower_row} upper_row = {upper_row} self.tomo_stacks[i].shape[0] = {self.tomo_stacks[i].shape[0]}')
1995 assert(0 <= lower_row < upper_row < self.tomo_stacks[i].shape[0]) 1988 assert(0 <= lower_row < upper_row < self.tomo_stacks[i].shape[0])
1996 center_offsets = [lower_center_offset-lower_row*center_slope, 1989 center_offsets = [lower_center_offset-lower_row*center_slope,
1997 upper_center_offset+(self.tomo_stacks[i].shape[0]-1-upper_row)*center_slope] 1990 upper_center_offset+(self.tomo_stacks[i].shape[0]-1-upper_row)*center_slope]
1998 t0 = time() 1991 t0 = time()
1999 print(f'OK4 {i} c')
2000 self.tomo_recon_stacks[i]= self._reconstructOneTomoStack(self.tomo_stacks[i], 1992 self.tomo_recon_stacks[i]= self._reconstructOneTomoStack(self.tomo_stacks[i],
2001 thetas, center_offsets=center_offsets, sigma=0.1, num_core=num_core, 1993 thetas, center_offsets=center_offsets, sigma=0.1, num_core=num_core,
2002 algorithm='gridrec', run_secondary_sirt=True, secondary_iter=25) 1994 algorithm='gridrec', run_secondary_sirt=True, secondary_iter=25)
2003 print(f'OK4 {i} d')
2004 logging.info(f'Reconstruction of stack {index} took {time()-t0:.2f} seconds!') 1995 logging.info(f'Reconstruction of stack {index} took {time()-t0:.2f} seconds!')
2005 if self.galaxy_flag: 1996 if self.galaxy_flag:
2006 x_slice = int(self.tomo_stacks[i].shape[0]/2) 1997 x_slice = int(self.tomo_stacks[i].shape[0]/2)
2007 title = f'{basetitle} {index} xslice{x_slice}' 1998 title = f'{basetitle} {index} xslice{x_slice}'
2008 msnc.quickImshow(self.tomo_recon_stacks[i][x_slice,:,:], title=title, 1999 msnc.quickImshow(self.tomo_recon_stacks[i][x_slice,:,:], title=title,
2033 stack['reconstructed'] = True 2024 stack['reconstructed'] = True
2034 combine_stacks = self.config.get('combine_stacks') 2025 combine_stacks = self.config.get('combine_stacks')
2035 if combine_stacks and index in combine_stacks.get('stacks', []): 2026 if combine_stacks and index in combine_stacks.get('stacks', []):
2036 combine_stacks['stacks'].remove(index) 2027 combine_stacks['stacks'].remove(index)
2037 self.cf.saveFile(self.config_out) 2028 self.cf.saveFile(self.config_out)
2038 print(f'OK4 {i} e')
2039 2029
2040 # Save reconstructed tomography stack to file 2030 # Save reconstructed tomography stack to file
2041 print('OK5')
2042 if self.galaxy_flag: 2031 if self.galaxy_flag:
2043 t0 = time() 2032 t0 = time()
2044 output_name = galaxy_param['output_name'] 2033 output_name = galaxy_param['output_name']
2045 logging.info(f'Saving reconstructed tomography stack to {output_name} ...') 2034 logging.info(f'Saving reconstructed tomography stack to {output_name} ...')
2046 save_stacks = {f'set_{stack["index"]}':tomo_stack 2035 save_stacks = {f'set_{stack["index"]}':tomo_stack
2047 for stack,tomo_stack in zip(stacks,self.tomo_recon_stacks)} 2036 for stack,tomo_stack in zip(stacks,self.tomo_recon_stacks)}
2048 np.savez(output_name, **save_stacks) 2037 np.savez(output_name, **save_stacks)
2049 logging.info(f'... done in {time()-t0:.2f} seconds!') 2038 logging.info(f'... done in {time()-t0:.2f} seconds!')
2050 print('OK6')
2051 2039
2052 def combineTomoStacks(self): 2040 def combineTomoStacks(self):
2053 """Combine the reconstructed tomography stacks. 2041 """Combine the reconstructed tomography stacks.
2054 """ 2042 """
2055 # stack order: stack,row(z),x,y 2043 # stack order: stack,row(z),x,y