Mercurial > repos > rv43 > tomo

comparison tomo.py @ 27:9d631a5e118f 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 aa8c5847f708ec917089859955382dfe3e82567e"
author rv43
date Mon, 18 Apr 2022 21:18:05 +0000
parents 7791b6e7d4e6
children c7211931fc99
comparison
equal deleted inserted replaced
26:7791b6e7d4e6 27:9d631a5e118f
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 1186 num_core=1
1187 print(f'OK2 numcore = {num_core}') 1187 #tomo_center = tomopy.find_center_vo(sinogram, ncore=1)
1188 tomo_center = tomopy.find_center_vo(sinogram, ncore=num_core) 1188 tomo_center = tomopy.find_center_vo(sinogram, ncore=num_core)
1189 center_offset_vo = tomo_center-center 1189 center_offset_vo = tomo_center-center
1190 print(f'center_offset_vo = {center_offset_vo}')
1191 if self.test_mode: 1190 if self.test_mode:
1192 logging.info(f'Center at row {row} using Nghia Vo’s method = {center_offset_vo:.2f}') 1191 logging.info(f'Center at row {row} using Nghia Vo’s method = {center_offset_vo:.2f}')
1193 del sinogram_T 1192 del sinogram_T
1194 return float(center_offset_vo) 1193 return float(center_offset_vo)
1195 elif self.galaxy_flag: 1194 elif self.galaxy_flag:
1197 recon_plane = self._reconstructOnePlane(sinogram_T, tomo_center, thetas_deg, 1196 recon_plane = self._reconstructOnePlane(sinogram_T, tomo_center, thetas_deg,
1198 eff_pixel_size, cross_sectional_dim, False, num_core) 1197 eff_pixel_size, cross_sectional_dim, False, num_core)
1199 title = f'edges row{row} center offset{center_offset_vo:.2f} Vo' 1198 title = f'edges row{row} center offset{center_offset_vo:.2f} Vo'
1200 self._plotEdgesOnePlane(recon_plane, title, path='find_center_pngs') 1199 self._plotEdgesOnePlane(recon_plane, title, path='find_center_pngs')
1201 del recon_plane 1200 del recon_plane
1202 print(f'center_type_selector = {galaxy_param["center_type_selector"]}')
1203 if not galaxy_param['center_type_selector']: 1201 if not galaxy_param['center_type_selector']:
1204 del sinogram_T 1202 del sinogram_T
1205 return float(center_offset_vo) 1203 return float(center_offset_vo)
1206 else: 1204 else:
1207 print(f'Center at row {row} using Nghia Vo’s method = {center_offset_vo:.2f}') 1205 print(f'Center at row {row} using Nghia Vo’s method = {center_offset_vo:.2f}')
1265 center_offset_step = pyip.inputInt('Enter step size for center offset search '+ 1263 center_offset_step = pyip.inputInt('Enter step size for center offset search '+
1266 f'[1, {center_offset_upp-center_offset_low}]: ', 1264 f'[1, {center_offset_upp-center_offset_low}]: ',
1267 min=1, max=center_offset_upp-center_offset_low) 1265 min=1, max=center_offset_upp-center_offset_low)
1268 num_center_offset = 1+int((center_offset_upp-center_offset_low)/center_offset_step) 1266 num_center_offset = 1+int((center_offset_upp-center_offset_low)/center_offset_step)
1269 center_offsets = np.linspace(center_offset_low, center_offset_upp, num_center_offset) 1267 center_offsets = np.linspace(center_offset_low, center_offset_upp, num_center_offset)
1270 print(f'center_offsets = {center_offsets}')
1271 for center_offset in center_offsets: 1268 for center_offset in center_offsets:
1272 if center_offset == center_offset_vo: 1269 if center_offset == center_offset_vo:
1273 continue 1270 continue
1274 logging.info(f'center_offset = {center_offset:.2f}') 1271 logging.info(f'center_offset = {center_offset:.2f}')
1275 recon_plane = self._reconstructOnePlane(sinogram_T, center_offset+center, 1272 recon_plane = self._reconstructOnePlane(sinogram_T, center_offset+center,
1277 title = f'edges row{row} center_offset{center_offset:.2f}' 1274 title = f'edges row{row} center_offset{center_offset:.2f}'
1278 if self.galaxy_flag: 1275 if self.galaxy_flag:
1279 self._plotEdgesOnePlane(recon_plane, title, path='find_center_pngs') 1276 self._plotEdgesOnePlane(recon_plane, title, path='find_center_pngs')
1280 else: 1277 else:
1281 self._plotEdgesOnePlane(recon_plane, title) 1278 self._plotEdgesOnePlane(recon_plane, title)
1282 print('OK3')
1283 if self.galaxy_flag or pyip.inputInt('\nContinue (0) or end the search (1): ', 1279 if self.galaxy_flag or pyip.inputInt('\nContinue (0) or end the search (1): ',
1284 min=0, max=1): 1280 min=0, max=1):
1285 break 1281 break
1286 1282
1287 del sinogram_T 1283 del sinogram_T
1738 if use_row != 'no': 1734 if use_row != 'no':
1739 center_offset = find_center.get('lower_center_offset') 1735 center_offset = find_center.get('lower_center_offset')
1740 if msnc.is_num(center_offset): 1736 if msnc.is_num(center_offset):
1741 use_center = pyip.inputYesNo('Current lower center offset = '+ 1737 use_center = pyip.inputYesNo('Current lower center offset = '+
1742 f'{center_offset}, use this value ([y]/n)? ', blank=True) 1738 f'{center_offset}, use this value ([y]/n)? ', blank=True)
1743 logging.info(f'use_center = {use_center}') 1739 #logging.info(f'use_center = {use_center}')
1744 logging.info(f'use_row = {use_row}') 1740 #logging.info(f'use_row = {use_row}')
1745 if use_center == 'no': 1741 if use_center == 'no':
1746 if use_row == 'no': 1742 if use_row == 'no':
1747 if not self.test_mode: 1743 if not self.test_mode:
1748 msnc.quickImshow(center_stack[:,0,:], title=f'theta={theta_start}', 1744 msnc.quickImshow(center_stack[:,0,:], title=f'theta={theta_start}',
1749 aspect='auto') 1745 aspect='auto')
1785 if use_row != 'no': 1781 if use_row != 'no':
1786 center_offset = find_center.get('upper_center_offset') 1782 center_offset = find_center.get('upper_center_offset')
1787 if msnc.is_num(center_offset): 1783 if msnc.is_num(center_offset):
1788 use_center = pyip.inputYesNo('Current upper center offset = '+ 1784 use_center = pyip.inputYesNo('Current upper center offset = '+
1789 f'{center_offset}, use this value ([y]/n)? ', blank=True) 1785 f'{center_offset}, use this value ([y]/n)? ', blank=True)
1790 logging.info(f'use_center = {use_center}') 1786 #logging.info(f'use_center = {use_center}')
1791 logging.info(f'use_row = {use_row}') 1787 #logging.info(f'use_row = {use_row}')
1792 if use_center == 'no': 1788 if use_center == 'no':
1793 if use_row == 'no': 1789 if use_row == 'no':
1794 if not self.test_mode: 1790 if not self.test_mode:
1795 msnc.quickImshow(center_stack[:,0,:], title=f'theta={theta_start}', 1791 msnc.quickImshow(center_stack[:,0,:], title=f'theta={theta_start}',
1796 aspect='auto') 1792 aspect='auto')
1899 logging.debug('Reconstruct tomography stacks') 1895 logging.debug('Reconstruct tomography stacks')
1900 stacks = self.config['stack_info']['stacks'] 1896 stacks = self.config['stack_info']['stacks']
1901 assert(len(self.tomo_stacks) == self.config['stack_info']['num']) 1897 assert(len(self.tomo_stacks) == self.config['stack_info']['num'])
1902 assert(len(self.tomo_stacks) == len(stacks)) 1898 assert(len(self.tomo_stacks) == len(stacks))
1903 assert(len(self.tomo_recon_stacks) == len(stacks)) 1899 assert(len(self.tomo_recon_stacks) == len(stacks))
1900 print('OK1')
1904 if self.galaxy_flag: 1901 if self.galaxy_flag:
1905 assert(isinstance(galaxy_param, dict)) 1902 assert(isinstance(galaxy_param, dict))
1906 # Get rotation axis centers 1903 # Get rotation axis centers
1907 center_offsets = galaxy_param['center_offsets'] 1904 center_offsets = galaxy_param['center_offsets']
1908 assert(isinstance(center_offsets, list) and len(center_offsets) == 2) 1905 assert(isinstance(center_offsets, list) and len(center_offsets) == 2)
1916 galaxy_param = None 1913 galaxy_param = None
1917 lower_center_offset = None 1914 lower_center_offset = None
1918 upper_center_offset = None 1915 upper_center_offset = None
1919 1916
1920 # Get rotation axis rows and centers 1917 # Get rotation axis rows and centers
1918 print('OK2')
1921 find_center = self.config['find_center'] 1919 find_center = self.config['find_center']
1922 lower_row = find_center.get('lower_row') 1920 lower_row = find_center.get('lower_row')
1923 if lower_row is None: 1921 if lower_row is None:
1924 logging.error('Unable to read lower_row from config') 1922 logging.error('Unable to read lower_row from config')
1925 return 1923 return
1949 num_theta_skip = self.config['preprocess'].get('num_theta_skip', 0) 1947 num_theta_skip = self.config['preprocess'].get('num_theta_skip', 0)
1950 thetas = np.radians(np.linspace(theta_start, theta_end, 1948 thetas = np.radians(np.linspace(theta_start, theta_end,
1951 int(num_theta/(num_theta_skip+1)), endpoint=False)) 1949 int(num_theta/(num_theta_skip+1)), endpoint=False))
1952 1950
1953 # Reconstruct tomo stacks 1951 # Reconstruct tomo stacks
1952 print('OK3')
1954 zoom_perc = self.config['preprocess'].get('zoom_perc', 100) 1953 zoom_perc = self.config['preprocess'].get('zoom_perc', 100)
1955 if zoom_perc == 100: 1954 if zoom_perc == 100:
1956 basetitle = 'recon stack fullres' 1955 basetitle = 'recon stack fullres'
1957 else: 1956 else:
1958 basetitle = f'recon stack {zoom_perc}p' 1957 basetitle = f'recon stack {zoom_perc}p'
1961 for i,stack in enumerate(stacks): 1960 for i,stack in enumerate(stacks):
1962 # Check if stack can be loaded 1961 # Check if stack can be loaded
1963 # reconstructed stack order for each one in stack : row/z,x,y 1962 # reconstructed stack order for each one in stack : row/z,x,y
1964 # preprocessed stack order for each one in stack: row,theta,column 1963 # preprocessed stack order for each one in stack: row,theta,column
1965 index = stack['index'] 1964 index = stack['index']
1965 print(f'OK4 {i} a')
1966 if not self.galaxy_flag: 1966 if not self.galaxy_flag:
1967 available = False 1967 available = False
1968 if stack.get('reconstructed', False): 1968 if stack.get('reconstructed', False):
1969 self.tomo_recon_stacks[i], available = self._loadTomo('recon stack', index) 1969 self.tomo_recon_stacks[i], available = self._loadTomo('recon stack', index)
1970 if self.tomo_recon_stacks[i].size or available: 1970 if self.tomo_recon_stacks[i].size or available:
1971 if self.tomo_stacks[i].size: 1971 if self.tomo_stacks[i].size:
1972 self.tomo_stacks[i] = np.array([]) 1972 self.tomo_stacks[i] = np.array([])
1973 assert(stack.get('preprocessed', False) == True) 1973 assert(stack.get('preprocessed', False) == True)
1974 assert(stack.get('reconstructed', False) == True) 1974 assert(stack.get('reconstructed', False) == True)
1975 continue 1975 continue
1976 print(f'OK4 {i} b')
1976 stack['reconstructed'] = False 1977 stack['reconstructed'] = False
1977 if not self.tomo_stacks[i].size: 1978 if not self.tomo_stacks[i].size:
1978 self.tomo_stacks[i], available = self._loadTomo('red stack', index, 1979 self.tomo_stacks[i], available = self._loadTomo('red stack', index,
1979 required=True) 1980 required=True)
1980 if not self.tomo_stacks[i].size: 1981 if not self.tomo_stacks[i].size:
1988 t0 = time() 1989 t0 = time()
1989 self.tomo_recon_stacks[i]= self._reconstructOneTomoStack(self.tomo_stacks[i], 1990 self.tomo_recon_stacks[i]= self._reconstructOneTomoStack(self.tomo_stacks[i],
1990 thetas, center_offsets=center_offsets, sigma=0.1, num_core=num_core, 1991 thetas, center_offsets=center_offsets, sigma=0.1, num_core=num_core,
1991 algorithm='gridrec', run_secondary_sirt=True, secondary_iter=25) 1992 algorithm='gridrec', run_secondary_sirt=True, secondary_iter=25)
1992 logging.info(f'Reconstruction of stack {index} took {time()-t0:.2f} seconds!') 1993 logging.info(f'Reconstruction of stack {index} took {time()-t0:.2f} seconds!')
1994 print(f'OK4 {i} c')
1993 if self.galaxy_flag: 1995 if self.galaxy_flag:
1994 x_slice = int(self.tomo_stacks[i].shape[0]/2) 1996 x_slice = int(self.tomo_stacks[i].shape[0]/2)
1995 title = f'{basetitle} {index} xslice{x_slice}' 1997 title = f'{basetitle} {index} xslice{x_slice}'
1996 msnc.quickImshow(self.tomo_recon_stacks[i][x_slice,:,:], title=title, 1998 msnc.quickImshow(self.tomo_recon_stacks[i][x_slice,:,:], title=title,
1997 path='center_slice_pngs', save_fig=True, save_only=True) 1999 path='center_slice_pngs', save_fig=True, save_only=True)
2021 stack['reconstructed'] = True 2023 stack['reconstructed'] = True
2022 combine_stacks = self.config.get('combine_stacks') 2024 combine_stacks = self.config.get('combine_stacks')
2023 if combine_stacks and index in combine_stacks.get('stacks', []): 2025 if combine_stacks and index in combine_stacks.get('stacks', []):
2024 combine_stacks['stacks'].remove(index) 2026 combine_stacks['stacks'].remove(index)
2025 self.cf.saveFile(self.config_out) 2027 self.cf.saveFile(self.config_out)
2028 print(f'OK4 {i} d')
2026 2029
2027 # Save reconstructed tomography stack to file 2030 # Save reconstructed tomography stack to file
2031 print('OK5')
2028 if self.galaxy_flag: 2032 if self.galaxy_flag:
2029 t0 = time() 2033 t0 = time()
2030 output_name = galaxy_param['output_name'] 2034 output_name = galaxy_param['output_name']
2031 logging.info(f'Saving reconstructed tomography stack to {output_name} ...') 2035 logging.info(f'Saving reconstructed tomography stack to {output_name} ...')
2032 save_stacks = {f'set_{stack["index"]}':tomo_stack 2036 save_stacks = {f'set_{stack["index"]}':tomo_stack
2033 for stack,tomo_stack in zip(stacks,self.tomo_recon_stacks)} 2037 for stack,tomo_stack in zip(stacks,self.tomo_recon_stacks)}
2034 np.savez(output_name, **save_stacks) 2038 np.savez(output_name, **save_stacks)
2035 logging.info(f'... done in {time()-t0:.2f} seconds!') 2039 logging.info(f'... done in {time()-t0:.2f} seconds!')
2040 print('OK6')
2036 2041
2037 def combineTomoStacks(self): 2042 def combineTomoStacks(self):
2038 """Combine the reconstructed tomography stacks. 2043 """Combine the reconstructed tomography stacks.
2039 """ 2044 """
2040 # stack order: stack,row(z),x,y 2045 # stack order: stack,row(z),x,y