Mercurial > repos > rv43 > tomo

comparison tomo.py @ 13:40395e60d2be 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 e61a62f5b5950b3142fedad2c0e49a00f66f9980"
author rv43
date Fri, 08 Apr 2022 17:49:28 +0000
parents ec7c3e84d611
children 1bcca1f2adb4
comparison
equal deleted inserted replaced
12:27d4a36a8c4a 13:40395e60d2be
1147 del recon_sinogram 1147 del recon_sinogram
1148 recon_clean = tomopy.misc.corr.remove_ring(recon_clean, rwidth=17, ncore=num_core) 1148 recon_clean = tomopy.misc.corr.remove_ring(recon_clean, rwidth=17, ncore=num_core)
1149 logging.debug(f'filtering and removing ring artifact took {time()-t0:.2f} seconds!') 1149 logging.debug(f'filtering and removing ring artifact took {time()-t0:.2f} seconds!')
1150 return recon_clean 1150 return recon_clean
1151 1151
1152 def _plotEdgesOnePlane(self, recon_plane, title, path=None, name=None, weight=0.001): 1152 def _plotEdgesOnePlane(self, recon_plane, title, path=None, weight=0.001):
1153 # RV parameters for the denoise, gaussian, and ring removal will be different for different feature sizes 1153 # RV parameters for the denoise, gaussian, and ring removal will be different for different feature sizes
1154 edges = denoise_tv_chambolle(recon_plane, weight = weight) 1154 edges = denoise_tv_chambolle(recon_plane, weight = weight)
1155 vmax = np.max(edges[0,:,:]) 1155 vmax = np.max(edges[0,:,:])
1156 vmin = -vmax 1156 vmin = -vmax
1157 if self.galaxy_flag: 1157 if self.galaxy_flag:
1158 msnc.quickImshow(edges[0,:,:], title, path=path, name=name, save_fig=True, 1158 msnc.quickImshow(edges[0,:,:], title, path=path, save_fig=True, save_only=True,
1159 save_only=True, cmap='gray', vmin=vmin, vmax=vmax) 1159 cmap='gray', vmin=vmin, vmax=vmax)
1160 else: 1160 else:
1161 if path is None: 1161 if path is None:
1162 path=self.output_folder 1162 path=self.output_folder
1163 msnc.quickImshow(edges[0,:,:], f'{title} coolwarm', path=path, 1163 msnc.quickImshow(edges[0,:,:], f'{title} coolwarm', path=path,
1164 save_fig=self.save_plots, save_only=self.save_plots_only, cmap='coolwarm') 1164 save_fig=self.save_plots, save_only=self.save_plots_only, cmap='coolwarm')
1166 save_fig=self.save_plots, save_only=self.save_plots_only, cmap='gray', 1166 save_fig=self.save_plots, save_only=self.save_plots_only, cmap='gray',
1167 vmin=vmin, vmax=vmax) 1167 vmin=vmin, vmax=vmax)
1168 del edges 1168 del edges
1169 1169
1170 def _findCenterOnePlane(self, sinogram, row, thetas_deg, eff_pixel_size, cross_sectional_dim, 1170 def _findCenterOnePlane(self, sinogram, row, thetas_deg, eff_pixel_size, cross_sectional_dim,
1171 tol=0.1, num_core=1, recon_pngname=None, galaxy_param=None): 1171 tol=0.1, num_core=1, galaxy_param=None):
1172 """Find center for a single tomography plane. 1172 """Find center for a single tomography plane.
1173 """ 1173 """
1174 if self.galaxy_flag:
1175 assert(galaxy_param)
1176 if not os.path.exists('find_center_pngs'):
1177 os.mkdir('find_center_pngs')
1174 # sinogram index order: theta,column 1178 # sinogram index order: theta,column
1175 # need index order column,theta for iradon, so take transpose 1179 # need index order column,theta for iradon, so take transpose
1176 sinogram_T = sinogram.T 1180 sinogram_T = sinogram.T
1177 center = sinogram.shape[1]/2 1181 center = sinogram.shape[1]/2
1178 1182
1183 logging.info(f'Center at row {row} using Nghia Vo’s method = {center_offset_vo:.2f}') 1187 logging.info(f'Center at row {row} using Nghia Vo’s method = {center_offset_vo:.2f}')
1184 del sinogram_T 1188 del sinogram_T
1185 return float(center_offset_vo) 1189 return float(center_offset_vo)
1186 elif self.galaxy_flag: 1190 elif self.galaxy_flag:
1187 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}')
1188 if recon_pngname: 1192 recon_plane = self._reconstructOnePlane(sinogram_T, tomo_center, thetas_deg,
1189 assert(isinstance(recon_pngname, str)) 1193 eff_pixel_size, cross_sectional_dim, False, num_core)
1190 recon_plane = self._reconstructOnePlane(sinogram_T, tomo_center, thetas_deg, 1194 title = f'edges row{row} center offset{center_offset_vo:.2f} Vo'
1191 eff_pixel_size, cross_sectional_dim, False, num_core) 1195 self._plotEdgesOnePlane(recon_plane, title, path='find_center_pngs')
1192 title = os.path.basename(recon_pngname) 1196 del recon_plane
1193 self._plotEdgesOnePlane(recon_plane, title, name=recon_pngname) 1197 if not galaxy_param['center_type_selector']:
1194 del recon_plane
1195 if not galaxy_param or not galaxy_param['center_type_selector']:
1196 del sinogram_T 1198 del sinogram_T
1197 return float(center_offset_vo) 1199 return float(center_offset_vo)
1198 else: 1200 else:
1199 print(f'Center at row {row} using Nghia Vo’s method = {center_offset_vo:.2f}') 1201 print(f'Center at row {row} using Nghia Vo’s method = {center_offset_vo:.2f}')
1200 if recon_pngname:
1201 logging.warning('Ignoring recon_pngname in _findCenterOnePlane (only for Galaxy)')
1202 recon_plane = self._reconstructOnePlane(sinogram_T, tomo_center, thetas_deg, 1202 recon_plane = self._reconstructOnePlane(sinogram_T, tomo_center, thetas_deg,
1203 eff_pixel_size, cross_sectional_dim, False, num_core) 1203 eff_pixel_size, cross_sectional_dim, False, num_core)
1204 title = f'edges row{row} center_offset_vo{center_offset_vo:.2f}' 1204 title = f'edges row{row} center offset{center_offset_vo:.2f} Vo'
1205 self._plotEdgesOnePlane(recon_plane, title) 1205 self._plotEdgesOnePlane(recon_plane, title)
1206 if not self.galaxy_flag: 1206 if not self.galaxy_flag:
1207 if (pyip.inputYesNo('Try finding center using phase correlation '+ 1207 if (pyip.inputYesNo('Try finding center using phase correlation '+
1208 '(y/[n])? ', blank=True) == 'yes'): 1208 '(y/[n])? ', blank=True) == 'yes'):
1209 tomo_center = tomopy.find_center_pc(sinogram, sinogram, tol=0.1, 1209 tomo_center = tomopy.find_center_pc(sinogram, sinogram, tol=0.1,
1216 error = np.abs(tomo_center-prev) 1216 error = np.abs(tomo_center-prev)
1217 center_offset = tomo_center-center 1217 center_offset = tomo_center-center
1218 print(f'Center at row {row} using phase correlation = {center_offset:.2f}') 1218 print(f'Center at row {row} using phase correlation = {center_offset:.2f}')
1219 recon_plane = self._reconstructOnePlane(sinogram_T, tomo_center, thetas_deg, 1219 recon_plane = self._reconstructOnePlane(sinogram_T, tomo_center, thetas_deg,
1220 eff_pixel_size, cross_sectional_dim, False, num_core) 1220 eff_pixel_size, cross_sectional_dim, False, num_core)
1221 title = f'edges row{row} center_offset{center_offset:.2f}' 1221 title = f'edges row{row} center_offset{center_offset:.2f} PC'
1222 self._plotEdgesOnePlane(recon_plane, title) 1222 self._plotEdgesOnePlane(recon_plane, title)
1223 if (pyip.inputYesNo('Accept a center location ([y]) or continue '+ 1223 if (pyip.inputYesNo('Accept a center location ([y]) or continue '+
1224 'search (n)? ', blank=True) != 'no'): 1224 'search (n)? ', blank=True) != 'no'):
1225 center_offset = pyip.inputNum( 1225 center_offset = pyip.inputNum(
1226 f' Enter chosen center offset [{-int(center)}, {int(center)}] '+ 1226 f' Enter chosen center offset [{-int(center)}, {int(center)}] '+
1230 del sinogram_T 1230 del sinogram_T
1231 del recon_plane 1231 del recon_plane
1232 return float(center_offset) 1232 return float(center_offset)
1233 1233
1234 # perform center finding search 1234 # perform center finding search
1235 if self.galaxy_flag and not os.path.exists('png_files'):
1236 os.mkdir('png_files')
1237 while True: 1235 while True:
1238 if self.galaxy_flag and galaxy_param and galaxy_param['center_type_selector']: 1236 if self.galaxy_flag and galaxy_param and galaxy_param['center_type_selector']:
1239 set_center = center_offset_vo 1237 set_center = center_offset_vo
1240 if galaxy_param['center_type_selector'] == 'user': 1238 if galaxy_param['center_type_selector'] == 'user':
1241 set_center = galaxy_param['set_center'] 1239 set_center = galaxy_param['set_center']
1269 logging.info(f'center_offset = {center_offset:.2f}') 1267 logging.info(f'center_offset = {center_offset:.2f}')
1270 recon_plane = self._reconstructOnePlane(sinogram_T, center_offset+center, 1268 recon_plane = self._reconstructOnePlane(sinogram_T, center_offset+center,
1271 thetas_deg, eff_pixel_size, cross_sectional_dim, False, num_core) 1269 thetas_deg, eff_pixel_size, cross_sectional_dim, False, num_core)
1272 title = f'edges row{row} center_offset{center_offset:.2f}' 1270 title = f'edges row{row} center_offset{center_offset:.2f}'
1273 if self.galaxy_flag: 1271 if self.galaxy_flag:
1274 self._plotEdgesOnePlane(recon_plane, title, path='png_files') 1272 self._plotEdgesOnePlane(recon_plane, title, path='find_center_pngs')
1275 else: 1273 else:
1276 self._plotEdgesOnePlane(recon_plane, title) 1274 self._plotEdgesOnePlane(recon_plane, title)
1277 if self.galaxy_flag or pyip.inputInt('\nContinue (0) or end the search (1): ', 1275 if self.galaxy_flag or pyip.inputInt('\nContinue (0) or end the search (1): ',
1278 min=0, max=1): 1276 min=0, max=1):
1279 break 1277 break
1758 if row == '': 1756 if row == '':
1759 row = n1 1757 row = n1
1760 if self.save_plots_only: 1758 if self.save_plots_only:
1761 msnc.clearFig(f'theta={theta_start}') 1759 msnc.clearFig(f'theta={theta_start}')
1762 # center_stack order: row,theta,column 1760 # center_stack order: row,theta,column
1763 if galaxy_param:
1764 recon_pngname = galaxy_param['recon_center_low']
1765 else:
1766 recon_pngname = None
1767 center_offset = self._findCenterOnePlane(center_stack[row,:,:], row, thetas_deg, 1761 center_offset = self._findCenterOnePlane(center_stack[row,:,:], row, thetas_deg,
1768 eff_pixel_size, cross_sectional_dim, num_core=num_core, 1762 eff_pixel_size, cross_sectional_dim, num_core=num_core,
1769 recon_pngname=recon_pngname, galaxy_param=galaxy_param) 1763 galaxy_param=galaxy_param)
1770 logging.info(f'lower_center_offset = {center_offset:.2f} {type(center_offset)}') 1764 logging.info(f'lower_center_offset = {center_offset:.2f} {type(center_offset)}')
1771 print(center_offset) 1765 print(center_offset)
1772 1766
1773 # Update config and save to file 1767 # Update config and save to file
1774 find_center['row_bounds'] = [n1, n2] 1768 find_center['row_bounds'] = [n1, n2]
1807 if row == '': 1801 if row == '':
1808 row = n2-1 1802 row = n2-1
1809 if self.save_plots_only: 1803 if self.save_plots_only:
1810 msnc.clearFig(f'theta={theta_start}') 1804 msnc.clearFig(f'theta={theta_start}')
1811 # center_stack order: row,theta,column 1805 # center_stack order: row,theta,column
1812 if galaxy_param:
1813 recon_pngname = galaxy_param['recon_center_upp']
1814 else:
1815 recon_pngname = None
1816 center_offset = self._findCenterOnePlane(center_stack[row,:,:], row, thetas_deg, 1806 center_offset = self._findCenterOnePlane(center_stack[row,:,:], row, thetas_deg,
1817 eff_pixel_size, cross_sectional_dim, num_core=num_core, 1807 eff_pixel_size, cross_sectional_dim, num_core=num_core,
1818 recon_pngname=recon_pngname, galaxy_param=galaxy_param) 1808 galaxy_param=galaxy_param)
1819 logging.info(f'upper_center_offset = {center_offset:.2f}') 1809 logging.info(f'upper_center_offset = {center_offset:.2f}')
1820 del center_stack 1810 del center_stack
1821 1811
1822 # Update config and save to file 1812 # Update config and save to file
1823 find_center['upper_row'] = row 1813 find_center['upper_row'] = row