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
2:b8977c98800b 3:f9c52762c32c
32 32
33 class set_numexpr_threads: 33 class set_numexpr_threads:
34 34
35 def __init__(self, num_core): 35 def __init__(self, num_core):
36 cpu_count = mp.cpu_count() 36 cpu_count = mp.cpu_count()
37 logging.debug(f'start: num_core={num_core} cpu_count={cpu_count}')
37 if num_core is None or num_core < 1 or num_core > cpu_count: 38 if num_core is None or num_core < 1 or num_core > cpu_count:
38 self.num_core = cpu_count 39 self.num_core = cpu_count
39 else: 40 else:
40 self.num_core = num_core 41 self.num_core = num_core
42 logging.debug(f'self.num_core={self.num_core}')
41 43
42 def __enter__(self): 44 def __enter__(self):
43 self.num_core_org = ne.set_num_threads(self.num_core) 45 self.num_core_org = ne.set_num_threads(self.num_core)
46 logging.debug(f'self.num_core={self.num_core}')
44 47
45 def __exit__(self, exc_type, exc_value, traceback): 48 def __exit__(self, exc_type, exc_value, traceback):
46 ne.set_num_threads(self.num_core_org) 49 ne.set_num_threads(self.num_core_org)
47 50
48 class ConfigTomo(msnc.Config): 51 class ConfigTomo(msnc.Config):
645 self.tdf[self.tdf > tdf_cutoff] = np.nan 648 self.tdf[self.tdf > tdf_cutoff] = np.nan
646 tdf_mean = np.nanmean(self.tdf) 649 tdf_mean = np.nanmean(self.tdf)
647 logging.debug(f'tdf_cutoff = {tdf_cutoff}') 650 logging.debug(f'tdf_cutoff = {tdf_cutoff}')
648 logging.debug(f'tdf_mean = {tdf_mean}') 651 logging.debug(f'tdf_mean = {tdf_mean}')
649 np.nan_to_num(self.tdf, copy=False, nan=tdf_mean, posinf=tdf_mean, neginf=0.) 652 np.nan_to_num(self.tdf, copy=False, nan=tdf_mean, posinf=tdf_mean, neginf=0.)
650 if not self.test_mode and not self.galaxy_flag: 653 if self.galaxy_flag:
654 msnc.quickImshow(self.tdf, title='dark field', name=dark_field_pngname,
655 save_fig=True, save_only=True)
656 elif not self.test_mode:
651 msnc.quickImshow(self.tdf, title='dark field', path=self.output_folder, 657 msnc.quickImshow(self.tdf, title='dark field', path=self.output_folder,
652 save_fig=self.save_plots, save_only=self.save_plots_only) 658 save_fig=self.save_plots, save_only=self.save_plots_only)
653 elif self.galaxy_flag:
654 msnc.quickImshow(self.tdf, title='dark field', name=dark_field_pngname,
655 save_fig=True, save_only=True)
656 659
657 def _genBright(self, tbf_files, bright_field_pngname): 660 def _genBright(self, tbf_files, bright_field_pngname):
658 """Generate bright field. 661 """Generate bright field.
659 """ 662 """
660 # Load the bright field images 663 # Load the bright field images
679 # Subtract dark field 682 # Subtract dark field
680 if self.tdf.size: 683 if self.tdf.size:
681 self.tbf -= self.tdf 684 self.tbf -= self.tdf
682 else: 685 else:
683 logging.warning('Dark field unavailable') 686 logging.warning('Dark field unavailable')
684 if not self.test_mode and not self.galaxy_flag: 687 if self.galaxy_flag:
688 msnc.quickImshow(self.tbf, title='bright field', name=bright_field_pngname,
689 save_fig=True, save_only=True)
690 elif not self.test_mode:
685 msnc.quickImshow(self.tbf, title='bright field', path=self.output_folder, 691 msnc.quickImshow(self.tbf, title='bright field', path=self.output_folder,
686 save_fig=self.save_plots, save_only=self.save_plots_only) 692 save_fig=self.save_plots, save_only=self.save_plots_only)
687 elif self.galaxy_flag:
688 msnc.quickImshow(self.tbf, title='bright field', name=bright_field_pngname,
689 save_fig=True, save_only=True)
690 693
691 def _setDetectorBounds(self, tomo_stack_files, tomo_field_pngname, detectorbounds_pngname): 694 def _setDetectorBounds(self, tomo_stack_files, tomo_field_pngname, detectorbounds_pngname):
692 """Set vertical detector bounds for image stack. 695 """Set vertical detector bounds for image stack.
693 """ 696 """
694 preprocess = self.config.get('preprocess') 697 preprocess = self.config.get('preprocess')
741 offset = min(0, int((self.tbf.shape[0]-num_x_min)/2-margin)) 744 offset = min(0, int((self.tbf.shape[0]-num_x_min)/2-margin))
742 img_x_bounds = [offset, num_x_min+offset+2*margin] 745 img_x_bounds = [offset, num_x_min+offset+2*margin]
743 tomo_stack = msnc.loadImageStack(tomo_stack_files[0], self.config['data_filetype'], 746 tomo_stack = msnc.loadImageStack(tomo_stack_files[0], self.config['data_filetype'],
744 stacks[0]['img_offset'], 1) 747 stacks[0]['img_offset'], 1)
745 x_sum = np.sum(tomo_stack[0,:,:], 1) 748 x_sum = np.sum(tomo_stack[0,:,:], 1)
749 x_sum_min = x_sum.min()
750 x_sum_max = x_sum.max()
746 title = f'tomography image at theta={self.config["theta_range"]["start"]}' 751 title = f'tomography image at theta={self.config["theta_range"]["start"]}'
747 msnc.quickImshow(tomo_stack[0,:,:], title=title, name=tomo_field_pngname, 752 msnc.quickImshow(tomo_stack[0,:,:], title=title, name=tomo_field_pngname,
748 save_fig=True, save_only=True) 753 save_fig=True, save_only=True, show_grid=True)
749 msnc.quickPlot((range(x_sum.size), x_sum), 754 msnc.quickPlot((range(x_sum.size), x_sum),
750 ([img_x_bounds[0], img_x_bounds[0]], [x_sum.min(), x_sum.max()], 'r-'), 755 ([img_x_bounds[0], img_x_bounds[0]], [x_sum_min, x_sum_max], 'r-'),
751 ([img_x_bounds[1]-1, img_x_bounds[1]-1], [x_sum.min(), x_sum.max()], 'r-'), 756 ([img_x_bounds[1]-1, img_x_bounds[1]-1], [x_sum_min, x_sum_max], 'r-'),
752 title='sum over theta and y', name=detectorbounds_pngname, 757 title='sum over theta and y', name=detectorbounds_pngname,
753 save_fig=True, save_only=True) 758 save_fig=True, save_only=True, show_grid=True)
754 759
755 # Update config and save to file 760 # Update config and save to file
756 if preprocess is None: 761 if preprocess is None:
757 self.cf.config['preprocess'] = {'img_x_bounds' : img_x_bounds} 762 self.cf.config['preprocess'] = {'img_x_bounds' : img_x_bounds}
758 else: 763 else:
759 preprocess['img_x_bounds'] = img_x_bounds 764 preprocess['img_x_bounds'] = img_x_bounds
760 self.cf.saveFile(self.config_out) 765 self.cf.saveFile(self.config_out)
766 del x_sum
761 return 767 return
762 768
763 # For one tomography stack only: load the first image 769 # For one tomography stack only: load the first image
764 title = None 770 title = None
765 msnc.quickImshow(self.tbf, title='bright field') 771 msnc.quickImshow(self.tbf, title='bright field')
775 tomo_or_bright = 0 781 tomo_or_bright = 0
776 if tomo_or_bright: 782 if tomo_or_bright:
777 x_sum = np.sum(tomo_stack[0,:,:], 1) 783 x_sum = np.sum(tomo_stack[0,:,:], 1)
778 use_bounds = 'no' 784 use_bounds = 'no'
779 if img_x_bounds[0] is not None and img_x_bounds[1] is not None: 785 if img_x_bounds[0] is not None and img_x_bounds[1] is not None:
786 tmp = np.copy(tomo_stack[0,:,:])
787 tmp_max = tmp.max()
788 tmp[img_x_bounds[0],:] = tmp_max
789 tmp[img_x_bounds[1]-1,:] = tmp_max
790 msnc.quickImshow(tmp, title=title)
791 del tmp
792 x_sum_min = x_sum.min()
793 x_sum_max = x_sum.max()
780 msnc.quickPlot((range(x_sum.size), x_sum), 794 msnc.quickPlot((range(x_sum.size), x_sum),
781 ([img_x_bounds[0], img_x_bounds[0]], [x_sum.min(), x_sum.max()], 'r-'), 795 ([img_x_bounds[0], img_x_bounds[0]], [x_sum_min, x_sum_max], 'r-'),
782 ([img_x_bounds[1]-1, img_x_bounds[1]-1], [x_sum.min(), x_sum.max()], 'r-'), 796 ([img_x_bounds[1]-1, img_x_bounds[1]-1], [x_sum_min, x_sum_max], 'r-'),
783 title='sum over theta and y') 797 title='sum over theta and y')
784 print(f'lower bound = {img_x_bounds[0]} (inclusive)\n'+ 798 print(f'lower bound = {img_x_bounds[0]} (inclusive)\n'+
785 f'upper bound = {img_x_bounds[1]} (exclusive)]') 799 f'upper bound = {img_x_bounds[1]} (exclusive)]')
786 use_bounds = pyip.inputYesNo('Accept these bounds ([y]/n)?: ', blank=True) 800 use_bounds = pyip.inputYesNo('Accept these bounds ([y]/n)?: ', blank=True)
787 if use_bounds == 'no': 801 if use_bounds == 'no':
797 x_sum[img_x_bounds[0]:img_x_bounds[1]], 811 x_sum[img_x_bounds[0]:img_x_bounds[1]],
798 title='sum over theta and y', path=self.output_folder, 812 title='sum over theta and y', path=self.output_folder,
799 save_fig=self.save_plots, save_only=True) 813 save_fig=self.save_plots, save_only=True)
800 else: 814 else:
801 x_sum = np.sum(self.tbf, 1) 815 x_sum = np.sum(self.tbf, 1)
816 x_sum_min = x_sum.min()
817 x_sum_max = x_sum.max()
802 use_bounds = 'no' 818 use_bounds = 'no'
803 if img_x_bounds[0] is not None and img_x_bounds[1] is not None: 819 if img_x_bounds[0] is not None and img_x_bounds[1] is not None:
820 tmp = np.copy(self.tbf)
821 tmp_max = tmp.max()
822 tmp[img_x_bounds[0],:] = tmp_max
823 tmp[img_x_bounds[1]-1,:] = tmp_max
824 title = 'Bright field'
825 msnc.quickImshow(tmp, title=title)
826 del tmp
804 msnc.quickPlot((range(x_sum.size), x_sum), 827 msnc.quickPlot((range(x_sum.size), x_sum),
805 ([img_x_bounds[0], img_x_bounds[0]], [x_sum.min(), x_sum.max()], 'r-'), 828 ([img_x_bounds[0], img_x_bounds[0]], [x_sum_min, x_sum_max], 'r-'),
806 ([img_x_bounds[1]-1, img_x_bounds[1]-1], [x_sum.min(), x_sum.max()], 'r-'), 829 ([img_x_bounds[1]-1, img_x_bounds[1]-1], [x_sum_min, x_sum_max], 'r-'),
807 title='sum over theta and y') 830 title='sum over theta and y')
808 print(f'lower bound = {img_x_bounds[0]} (inclusive)\n'+ 831 print(f'lower bound = {img_x_bounds[0]} (inclusive)\n'+
809 f'upper bound = {img_x_bounds[1]} (exclusive)]') 832 f'upper bound = {img_x_bounds[1]} (exclusive)]')
810 use_bounds = pyip.inputYesNo('Accept these bounds ([y]/n)?: ', blank=True) 833 use_bounds = pyip.inputYesNo('Accept these bounds ([y]/n)?: ', blank=True)
811 if use_bounds == 'no': 834 if use_bounds == 'no':
818 if x_low < 0: 841 if x_low < 0:
819 x_low = 0 842 x_low = 0
820 x_upp = int(x_upp+(x_upp-x_low)/10) 843 x_upp = int(x_upp+(x_upp-x_low)/10)
821 if x_upp >= x_sum.size: 844 if x_upp >= x_sum.size:
822 x_upp = x_sum.size 845 x_upp = x_sum.size
846 tmp = np.copy(self.tbf)
847 tmp_max = tmp.max()
848 tmp[x_low,:] = tmp_max
849 tmp[x_upp-1,:] = tmp_max
850 title = 'Bright field'
851 msnc.quickImshow(tmp, title=title)
852 del tmp
823 msnc.quickPlot((range(x_sum.size), x_sum), 853 msnc.quickPlot((range(x_sum.size), x_sum),
824 ([x_low, x_low], [x_sum.min(), x_sum.max()], 'r-'), 854 ([x_low, x_low], [x_sum_min, x_sum_max], 'r-'),
825 ([x_upp, x_upp], [x_sum.min(), x_sum.max()], 'r-'), 855 ([x_upp, x_upp], [x_sum_min, x_sum_max], 'r-'),
826 title='sum over theta and y') 856 title='sum over theta and y')
827 print(f'lower bound = {x_low} (inclusive)\nupper bound = {x_upp} (exclusive)]') 857 print(f'lower bound = {x_low} (inclusive)\nupper bound = {x_upp} (exclusive)]')
828 use_fit = pyip.inputYesNo('Accept these bounds ([y]/n)?: ', blank=True) 858 use_fit = pyip.inputYesNo('Accept these bounds ([y]/n)?: ', blank=True)
829 if use_fit == 'no': 859 if use_fit == 'no':
830 img_x_bounds = msnc.selectArrayBounds(x_sum, img_x_bounds[0], img_x_bounds[1], 860 img_x_bounds = msnc.selectArrayBounds(x_sum, img_x_bounds[0], img_x_bounds[1],
835 logging.warning('Image bounds and pixel size prevent seamless stacking') 865 logging.warning('Image bounds and pixel size prevent seamless stacking')
836 msnc.quickPlot(range(img_x_bounds[0], img_x_bounds[1]), 866 msnc.quickPlot(range(img_x_bounds[0], img_x_bounds[1]),
837 x_sum[img_x_bounds[0]:img_x_bounds[1]], 867 x_sum[img_x_bounds[0]:img_x_bounds[1]],
838 title='sum over theta and y', path=self.output_folder, 868 title='sum over theta and y', path=self.output_folder,
839 save_fig=self.save_plots, save_only=True) 869 save_fig=self.save_plots, save_only=True)
870 del x_sum
840 logging.debug(f'img_x_bounds: {img_x_bounds}') 871 logging.debug(f'img_x_bounds: {img_x_bounds}')
841 872
842 if self.save_plots_only: 873 if self.save_plots_only:
843 msnc.clearFig('bright field') 874 msnc.clearFig('bright field')
844 msnc.clearFig('sum over theta and y') 875 msnc.clearFig('sum over theta and y')
1067 # Update config and save to file 1098 # Update config and save to file
1068 stack['preprocessed'] = True 1099 stack['preprocessed'] = True
1069 stack.pop('reconstructed', 'reconstructed not found') 1100 stack.pop('reconstructed', 'reconstructed not found')
1070 find_center = self.config.get('find_center') 1101 find_center = self.config.get('find_center')
1071 if find_center: 1102 if find_center:
1103 find_center.pop('completed', 'completed not found')
1072 if self.test_mode: 1104 if self.test_mode:
1073 find_center.pop('completed', 'completed not found')
1074 find_center.pop('lower_center_offset', 'lower_center_offset not found') 1105 find_center.pop('lower_center_offset', 'lower_center_offset not found')
1075 find_center.pop('upper_center_offset', 'upper_center_offset not found') 1106 find_center.pop('upper_center_offset', 'upper_center_offset not found')
1076 else:
1077 if find_center.get('center_stack_index', -1) == index:
1078 self.config.pop('find_center')
1079 self.cf.saveFile(self.config_out) 1107 self.cf.saveFile(self.config_out)
1080 1108
1081 if self.tdf.size: 1109 if self.tdf.size:
1082 del tdf 1110 del tdf
1083 del tbf 1111 del tbf
1084 1112
1085 def _reconstructOnePlane(self, tomo_plane_T, center, thetas_deg, eff_pixel_size, 1113 def _reconstructOnePlane(self, tomo_plane_T, center, thetas_deg, eff_pixel_size,
1086 cross_sectional_dim, plot_sinogram=True): 1114 cross_sectional_dim, plot_sinogram=True, num_core=1):
1087 """Invert the sinogram for a single tomography plane. 1115 """Invert the sinogram for a single tomography plane.
1088 """ 1116 """
1089 # tomo_plane_T index order: column,theta 1117 # tomo_plane_T index order: column,theta
1090 assert(0 <= center < tomo_plane_T.shape[0]) 1118 assert(0 <= center < tomo_plane_T.shape[0])
1091 center_offset = center-tomo_plane_T.shape[0]/2 1119 center_offset = center-tomo_plane_T.shape[0]/2
1115 t0 = time() 1143 t0 = time()
1116 # recon_sinogram = filters.gaussian(recon_sinogram, 3.0) 1144 # recon_sinogram = filters.gaussian(recon_sinogram, 3.0)
1117 recon_sinogram = spi.gaussian_filter(recon_sinogram, 0.5) 1145 recon_sinogram = spi.gaussian_filter(recon_sinogram, 0.5)
1118 recon_clean = np.expand_dims(recon_sinogram, axis=0) 1146 recon_clean = np.expand_dims(recon_sinogram, axis=0)
1119 del recon_sinogram 1147 del recon_sinogram
1120 recon_clean = tomopy.misc.corr.remove_ring(recon_clean, rwidth=17) 1148 recon_clean = tomopy.misc.corr.remove_ring(recon_clean, rwidth=17, ncore=num_core)
1121 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!')
1122 return recon_clean 1150 return recon_clean
1123 1151
1124 def _plotEdgesOnePlane(self, recon_plane, base_name, weight=0.001): 1152 def _plotEdgesOnePlane(self, recon_plane, title, name=None, weight=0.001):
1125 # 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
1126 edges = denoise_tv_chambolle(recon_plane, weight = weight) 1154 edges = denoise_tv_chambolle(recon_plane, weight = weight)
1127 vmax = np.max(edges[0,:,:]) 1155 vmax = np.max(edges[0,:,:])
1128 vmin = -vmax 1156 vmin = -vmax
1129 msnc.quickImshow(edges[0,:,:], f'{base_name} coolwarm', path=self.output_folder, 1157 if self.galaxy_flag:
1130 save_fig=self.save_plots, save_only=self.save_plots_only, cmap='coolwarm') 1158 msnc.quickImshow(edges[0,:,:], title, name=name, save_fig=True, save_only=True,
1131 msnc.quickImshow(edges[0,:,:], f'{base_name} gray', path=self.output_folder, 1159 cmap='gray', vmin=vmin, vmax=vmax)
1132 save_fig=self.save_plots, save_only=self.save_plots_only, cmap='gray', 1160 else:
1133 vmin=vmin, vmax=vmax) 1161 msnc.quickImshow(edges[0,:,:], f'{title} coolwarm', path=self.output_folder,
1162 save_fig=self.save_plots, save_only=self.save_plots_only, cmap='coolwarm')
1163 msnc.quickImshow(edges[0,:,:], f'{title} gray', path=self.output_folder,
1164 save_fig=self.save_plots, save_only=self.save_plots_only, cmap='gray',
1165 vmin=vmin, vmax=vmax)
1134 del edges 1166 del edges
1135 1167
1136 def _findCenterOnePlane(self, sinogram, row, thetas_deg, eff_pixel_size, cross_sectional_dim, 1168 def _findCenterOnePlane(self, sinogram, row, thetas_deg, eff_pixel_size, cross_sectional_dim,
1137 tol=0.1): 1169 tol=0.1, num_core=1, recon_pngname=None):
1138 """Find center for a single tomography plane. 1170 """Find center for a single tomography plane.
1139 """ 1171 """
1140 # sinogram index order: theta,column 1172 # sinogram index order: theta,column
1141 # need index order column,theta for iradon, so take transpose 1173 # need index order column,theta for iradon, so take transpose
1142 sinogram_T = sinogram.T 1174 sinogram_T = sinogram.T
1143 center = sinogram.shape[1]/2 1175 center = sinogram.shape[1]/2
1144 1176
1145 # try automatic center finding routines for initial value 1177 # try automatic center finding routines for initial value
1146 tomo_center = tomopy.find_center_vo(sinogram) 1178 tomo_center = tomopy.find_center_vo(sinogram, ncore=num_core)
1147 center_offset_vo = tomo_center-center 1179 center_offset_vo = tomo_center-center
1148 if not self.test_mode: 1180 if self.test_mode or self.galaxy_flag:
1181 logging.info(f'Center at row {row} using Nghia Vo’s method = {center_offset_vo:.2f}')
1182 if self.test_mode:
1183 del sinogram_T
1184 return float(center_offset_vo)
1185 else:
1149 print(f'Center at row {row} using Nghia Vo’s method = {center_offset_vo:.2f}') 1186 print(f'Center at row {row} using Nghia Vo’s method = {center_offset_vo:.2f}')
1187 if recon_pngname:
1188 logging.warning('Ignoring recon_pngname in _findCenterOnePlane (only for Galaxy)')
1150 recon_plane = self._reconstructOnePlane(sinogram_T, tomo_center, thetas_deg, 1189 recon_plane = self._reconstructOnePlane(sinogram_T, tomo_center, thetas_deg,
1151 eff_pixel_size, cross_sectional_dim, False) 1190 eff_pixel_size, cross_sectional_dim, False, num_core)
1152 if not self.test_mode: 1191 if self.galaxy_flag:
1153 base_name=f'edges row{row} center_offset_vo{center_offset_vo:.2f}' 1192 assert(isinstance(recon_pngname, str))
1154 self._plotEdgesOnePlane(recon_plane, base_name) 1193 title = os.path.basename(recon_pngname)
1155 use_phase_corr = 'no' 1194 self._plotEdgesOnePlane(recon_plane, title, name=recon_pngname)
1156 if not self.test_mode: 1195 del sinogram_T
1157 use_phase_corr = pyip.inputYesNo('Try finding center using phase correlation '+ 1196 del recon_plane
1158 '(y/[n])? ', blank=True) 1197 return float(center_offset_vo)
1159 if use_phase_corr == 'yes': 1198 else:
1199 title = f'edges row{row} center_offset_vo{center_offset_vo:.2f}'
1200 self._plotEdgesOnePlane(recon_plane, title)
1201 if (pyip.inputYesNo('Try finding center using phase correlation '+
1202 '(y/[n])? ', blank=True) == 'yes'):
1160 tomo_center = tomopy.find_center_pc(sinogram, sinogram, tol=0.1, 1203 tomo_center = tomopy.find_center_pc(sinogram, sinogram, tol=0.1,
1161 rotc_guess=tomo_center) 1204 rotc_guess=tomo_center)
1162 error = 1. 1205 error = 1.
1163 while error > tol: 1206 while error > tol:
1164 prev = tomo_center 1207 prev = tomo_center
1166 rotc_guess=tomo_center) 1209 rotc_guess=tomo_center)
1167 error = np.abs(tomo_center-prev) 1210 error = np.abs(tomo_center-prev)
1168 center_offset = tomo_center-center 1211 center_offset = tomo_center-center
1169 print(f'Center at row {row} using phase correlation = {center_offset:.2f}') 1212 print(f'Center at row {row} using phase correlation = {center_offset:.2f}')
1170 recon_plane = self._reconstructOnePlane(sinogram_T, tomo_center, thetas_deg, 1213 recon_plane = self._reconstructOnePlane(sinogram_T, tomo_center, thetas_deg,
1171 eff_pixel_size, cross_sectional_dim, False) 1214 eff_pixel_size, cross_sectional_dim, False, num_core)
1172 base_name=f'edges row{row} center_offset{center_offset:.2f}' 1215 title = f'edges row{row} center_offset{center_offset:.2f}'
1173 self._plotEdgesOnePlane(recon_plane, base_name) 1216 self._plotEdgesOnePlane(recon_plane, title)
1174 accept_center = 'yes' 1217 if (pyip.inputYesNo('Accept a center location ([y]) or continue '+
1175 if not self.test_mode: 1218 'search (n)? ', blank=True) != 'no'):
1176 accept_center = pyip.inputYesNo('Accept a center location ([y]) or continue '+
1177 'search (n)? ', blank=True)
1178 if accept_center != 'no':
1179 del sinogram_T 1219 del sinogram_T
1180 del recon_plane 1220 del recon_plane
1181 if self.test_mode: 1221 center_offset = pyip.inputNum(
1222 f' Enter chosen center offset [{-int(center)}, {int(center)}] '+
1223 f'([{center_offset_vo}])): ', blank=True)
1224 if center_offset == '':
1182 center_offset = center_offset_vo 1225 center_offset = center_offset_vo
1183 else:
1184 center_offset = pyip.inputNum(
1185 f' Enter chosen center offset [{-int(center)}, {int(center)}] '+
1186 f'([{center_offset_vo}])): ', blank=True)
1187 if center_offset == '':
1188 center_offset = center_offset_vo
1189 return float(center_offset) 1226 return float(center_offset)
1190 1227
1191 while True: 1228 while True:
1192 center_offset_low = pyip.inputInt('\nEnter lower bound for center offset '+ 1229 center_offset_low = pyip.inputInt('\nEnter lower bound for center offset '+
1193 f'[{-int(center)}, {int(center)}]: ', min=-int(center), max=int(center)) 1230 f'[{-int(center)}, {int(center)}]: ', min=-int(center), max=int(center))
1202 min=1, max=center_offset_upp-center_offset_low) 1239 min=1, max=center_offset_upp-center_offset_low)
1203 for center_offset in range(center_offset_low, center_offset_upp+center_offset_step, 1240 for center_offset in range(center_offset_low, center_offset_upp+center_offset_step,
1204 center_offset_step): 1241 center_offset_step):
1205 logging.info(f'center_offset = {center_offset}') 1242 logging.info(f'center_offset = {center_offset}')
1206 recon_plane = self._reconstructOnePlane(sinogram_T, center_offset+center, 1243 recon_plane = self._reconstructOnePlane(sinogram_T, center_offset+center,
1207 thetas_deg, eff_pixel_size, cross_sectional_dim, False) 1244 thetas_deg, eff_pixel_size, cross_sectional_dim, False, num_core)
1208 base_name=f'edges row{row} center_offset{center_offset}' 1245 title = f'edges row{row} center_offset{center_offset}'
1209 self._plotEdgesOnePlane(recon_plane, base_name) 1246 self._plotEdgesOnePlane(recon_plane, title)
1210 if pyip.inputInt('\nContinue (0) or end the search (1): ', min=0, max=1): 1247 if pyip.inputInt('\nContinue (0) or end the search (1): ', min=0, max=1):
1211 break 1248 break
1212 1249
1213 del sinogram_T 1250 del sinogram_T
1214 del recon_plane 1251 del recon_plane
1264 'num_iter':secondary_iter} 1301 'num_iter':secondary_iter}
1265 tomo_recon_stack = tomopy.recon(tomo_stack, thetas, centers, 1302 tomo_recon_stack = tomopy.recon(tomo_stack, thetas, centers,
1266 init_recon=tomo_recon_stack, options=options, sinogram_order=True, 1303 init_recon=tomo_recon_stack, options=options, sinogram_order=True,
1267 algorithm=tomopy.astra, ncore=num_core) 1304 algorithm=tomopy.astra, ncore=num_core)
1268 if True: 1305 if True:
1269 tomopy.misc.corr.remove_ring(tomo_recon_stack, rwidth=rwidth, out=tomo_recon_stack) 1306 tomopy.misc.corr.remove_ring(tomo_recon_stack, rwidth=rwidth, out=tomo_recon_stack,
1307 ncore=num_core)
1270 return tomo_recon_stack 1308 return tomo_recon_stack
1271 1309
1272 def findImageFiles(self): 1310 def findImageFiles(self):
1273 """Find all available image files. 1311 """Find all available image files.
1274 """ 1312 """
1356 # Safe updated config 1394 # Safe updated config
1357 if self.is_valid: 1395 if self.is_valid:
1358 self.cf.saveFile(self.config_out) 1396 self.cf.saveFile(self.config_out)
1359 1397
1360 return dark_files, bright_files, tomo_stack_files 1398 return dark_files, bright_files, tomo_stack_files
1399
1400 def loadTomoStacks(self, input_name):
1401 """Load tomography stacks (only for Galaxy).
1402 """
1403 assert(self.galaxy_flag)
1404 t0 = time()
1405 logging.info(f'Loading preprocessed tomography stack from {input_name} ...')
1406 stack_info = self.config['stack_info']
1407 stacks = stack_info.get('stacks')
1408 assert(len(self.tomo_stacks) == stack_info['num'])
1409 with np.load(input_name) as f:
1410 for i,stack in enumerate(stacks):
1411 self.tomo_stacks[i] = f[f'set_{stack["index"]}']
1412 logging.info(f'loaded stack {i}: index = {stack["index"]}, shape = '+
1413 f'{self.tomo_stacks[i].shape}')
1414 logging.info(f'... done in {time()-t0:.2f} seconds!')
1361 1415
1362 def genTomoStacks(self, tdf_files=None, tbf_files=None, tomo_stack_files=None, 1416 def genTomoStacks(self, tdf_files=None, tbf_files=None, tomo_stack_files=None,
1363 dark_field_pngname=None, bright_field_pngname=None, tomo_field_pngname=None, 1417 dark_field_pngname=None, bright_field_pngname=None, tomo_field_pngname=None,
1364 detectorbounds_pngname=None, output_name=None): 1418 detectorbounds_pngname=None, output_name=None):
1365 """Preprocess tomography images. 1419 """Preprocess tomography images.
1457 pixel_size *= img_x_bounds[0] 1511 pixel_size *= img_x_bounds[0]
1458 for stack in stacks: 1512 for stack in stacks:
1459 stack['ref_height'] = stack['ref_height']+pixel_size 1513 stack['ref_height'] = stack['ref_height']+pixel_size
1460 self.cf.saveFile(self.config_out) 1514 self.cf.saveFile(self.config_out)
1461 1515
1462 def findCenters(self): 1516 def findCenters(self, row_bounds=None, center_rows=None, recon_low_pngname=None,
1517 recon_upp_pngname=None, num_core=None):
1463 """Find rotation axis centers for the tomography stacks. 1518 """Find rotation axis centers for the tomography stacks.
1464 """ 1519 """
1520 if num_core is None:
1521 num_core = self.num_core
1465 logging.debug('Find centers for tomography stacks') 1522 logging.debug('Find centers for tomography stacks')
1466 stacks = self.config['stack_info']['stacks'] 1523 stacks = self.config['stack_info']['stacks']
1467 available_stacks = [stack['index'] for stack in stacks if stack.get('preprocessed', False)] 1524 available_stacks = [stack['index'] for stack in stacks if stack.get('preprocessed', False)]
1468 logging.debug('Available stacks: {available_stacks}') 1525 logging.debug('Available stacks: {available_stacks}')
1526 if self.galaxy_flag:
1527 assert(isinstance(row_bounds, list) and len(row_bounds) == 2)
1528 assert(isinstance(center_rows, list) and len(center_rows) == 2)
1529 assert(isinstance(recon_low_pngname, str))
1530 assert(isinstance(recon_upp_pngname, str))
1531 else:
1532 if row_bounds:
1533 logging.warning('Ignoring row_bounds in findCenters (only for Galaxy)')
1534 row_bounds = None
1535 if center_rows:
1536 logging.warning('Ignoring center_rows in findCenters (only for Galaxy)')
1537 center_rows = None
1538 if recon_low_pngname:
1539 logging.warning('Ignoring recon_low_pngname in findCenters (only for Galaxy)')
1540 recon_low_pngname = None
1541 if recon_upp_pngname:
1542 logging.warning('Ignoring recon_upp_pngname in findCenters (only for Galaxy)')
1543 recon_upp_pngname = None
1469 1544
1470 # Check for valid available center stack index 1545 # Check for valid available center stack index
1471 find_center = self.config.get('find_center') 1546 find_center = self.config.get('find_center')
1472 center_stack_index = None 1547 center_stack_index = None
1473 if find_center and 'center_stack_index' in find_center: 1548 if find_center and 'center_stack_index' in find_center:
1479 if self.test_mode: 1554 if self.test_mode:
1480 assert(find_center.get('completed', False) == False) 1555 assert(find_center.get('completed', False) == False)
1481 else: 1556 else:
1482 print('\nFound calibration center offset info for stack '+ 1557 print('\nFound calibration center offset info for stack '+
1483 f'{center_stack_index}') 1558 f'{center_stack_index}')
1484 if (pyip.inputYesNo('Do you want to use this again (y/n)? ') == 'yes' and 1559 if (pyip.inputYesNo('Do you want to use this again ([y]/n)? ',
1485 find_center.get('completed', False) == True): 1560 blank=True) != 'no' and find_center.get('completed', False) == True):
1486 return 1561 return
1487 1562
1488 # Load the required preprocessed stack 1563 # Load the required preprocessed stack
1489 # preprocessed stack order: row,theta,column 1564 # preprocessed stack order: row,theta,column
1490 num_tomo_stacks = self.config['stack_info']['num'] 1565 num_tomo_stacks = self.config['stack_info']['num']
1501 center_stack = self.tomo_stacks[0] 1576 center_stack = self.tomo_stacks[0]
1502 if not center_stack.size: 1577 if not center_stack.size:
1503 stacks[0]['preprocessed'] = False 1578 stacks[0]['preprocessed'] = False
1504 raise OSError('Unable to load the required preprocessed tomography stack') 1579 raise OSError('Unable to load the required preprocessed tomography stack')
1505 assert(stacks[0].get('preprocessed', False) == True) 1580 assert(stacks[0].get('preprocessed', False) == True)
1581 elif self.galaxy_flag:
1582 logging.error('CHECK/FIX FOR GALAXY')
1583 #center_stack_index = stacks[int(num_tomo_stacks/2)]['index']
1506 else: 1584 else:
1507 while True: 1585 while True:
1508 if not center_stack_index: 1586 if not center_stack_index:
1509 center_stack_index = pyip.inputInt('\nEnter tomography stack index to get ' 1587 center_stack_index = pyip.inputInt('\nEnter tomography stack index to get '
1510 f'rotation axis centers {available_stacks}: ') 1588 f'rotation axis centers {available_stacks}: ')
1526 if find_center is None: 1604 if find_center is None:
1527 self.config['find_center'] = {'center_stack_index' : center_stack_index} 1605 self.config['find_center'] = {'center_stack_index' : center_stack_index}
1528 find_center = self.config['find_center'] 1606 find_center = self.config['find_center']
1529 else: 1607 else:
1530 find_center['center_stack_index'] = center_stack_index 1608 find_center['center_stack_index'] = center_stack_index
1609 if not self.galaxy_flag:
1610 row_bounds = find_center.get('row_bounds', None)
1611 center_rows = [find_center.get('lower_row', None),
1612 find_center.get('upper_row', None)]
1531 1613
1532 # Set thetas (in degrees) 1614 # Set thetas (in degrees)
1533 theta_range = self.config['theta_range'] 1615 theta_range = self.config['theta_range']
1534 theta_start = theta_range['start'] 1616 theta_start = theta_range['start']
1535 theta_end = theta_range['end'] 1617 theta_end = theta_range['end']
1543 pixel_size = self.config['detector']['pixel_size'] 1625 pixel_size = self.config['detector']['pixel_size']
1544 if pixel_size is None: 1626 if pixel_size is None:
1545 raise ValueError('Detector pixel size unavailable') 1627 raise ValueError('Detector pixel size unavailable')
1546 eff_pixel_size = 100.*pixel_size/zoom_perc 1628 eff_pixel_size = 100.*pixel_size/zoom_perc
1547 logging.debug(f'eff_pixel_size = {eff_pixel_size}') 1629 logging.debug(f'eff_pixel_size = {eff_pixel_size}')
1548 tomo_ref_heights = [stack['ref_height'] for stack in stacks]
1549 if num_tomo_stacks == 1: 1630 if num_tomo_stacks == 1:
1550 n1 = 0 1631 accept = 'yes'
1551 height = center_stack.shape[0]*eff_pixel_size 1632 if not self.test_mode and not self.galaxy_flag:
1552 if not self.test_mode and pyip.inputYesNo('\nDo you want to reconstruct the full '+ 1633 accept = 'no'
1553 f'sample height ({height:.3f} mm) (y/n)? ') == 'no': 1634 print('\nSelect bounds for image reconstruction')
1554 height = pyip.inputNum('\nEnter the desired reconstructed sample height '+ 1635 if msnc.is_index_range(row_bounds, 0, center_stack.shape[0]):
1555 f'in mm [0, {height:.3f}]: ', min=0, max=height) 1636 a_tmp = np.copy(center_stack[:,0,:])
1556 n1 = int(0.5*(center_stack.shape[0]-height/eff_pixel_size)) 1637 a_tmp_max = a_tmp.max()
1557 else: 1638 a_tmp[row_bounds[0],:] = a_tmp_max
1639 a_tmp[row_bounds[1]-1,:] = a_tmp_max
1640 print(f'lower bound = {row_bounds[0]} (inclusive)')
1641 print(f'upper bound = {row_bounds[1]} (exclusive)')
1642 msnc.quickImshow(a_tmp, title=f'center stack theta={theta_start}',
1643 aspect='auto')
1644 del a_tmp
1645 accept = pyip.inputYesNo('Accept these bounds ([y]/n)?: ', blank=True)
1646 if accept == 'no':
1647 (n1, n2) = msnc.selectImageBounds(center_stack[:,0,:], 0,
1648 title=f'center stack theta={theta_start}')
1649 else:
1650 n1 = row_bounds[0]
1651 n2 = row_bounds[1]
1652 else:
1653 logging.error('CHECK/FIX FOR GALAXY')
1654 tomo_ref_heights = [stack['ref_height'] for stack in stacks]
1558 n1 = int((1.+(tomo_ref_heights[0]+center_stack.shape[0]*eff_pixel_size- 1655 n1 = int((1.+(tomo_ref_heights[0]+center_stack.shape[0]*eff_pixel_size-
1559 tomo_ref_heights[1])/eff_pixel_size)/2) 1656 tomo_ref_heights[1])/eff_pixel_size)/2)
1560 n2 = center_stack.shape[0]-n1 1657 n2 = center_stack.shape[0]-n1
1561 logging.info(f'n1 = {n1}, n2 = {n2} (n2-n1) = {(n2-n1)*eff_pixel_size:.3f} mm') 1658 logging.info(f'n1 = {n1}, n2 = {n2} (n2-n1) = {(n2-n1)*eff_pixel_size:.3f} mm')
1562 if not center_stack.size: 1659 if not center_stack.size:
1563 RuntimeError('Center stack not loaded') 1660 RuntimeError('Center stack not loaded')
1564 if not self.test_mode: 1661 if not self.test_mode and not self.galaxy_flag:
1565 msnc.quickImshow(center_stack[:,0,:], title=f'center stack theta={theta_start}', 1662 tmp = center_stack[:,0,:]
1566 path=self.output_folder, save_fig=self.save_plots, 1663 tmp_max = tmp.max()
1567 save_only=self.save_plots_only) 1664 tmp[n1,:] = tmp_max
1665 tmp[n2-1,:] = tmp_max
1666 if msnc.is_index_range(center_rows, 0, tmp.shape[0]):
1667 tmp[center_rows[0],:] = tmp_max
1668 tmp[center_rows[1]-1,:] = tmp_max
1669 extent = [0, tmp.shape[1], tmp.shape[0], 0]
1670 msnc.quickImshow(tmp, title=f'center stack theta={theta_start}',
1671 path=self.output_folder, extent=extent, save_fig=self.save_plots,
1672 save_only=self.save_plots_only, aspect='auto')
1673 del tmp
1674 #extent = [0, center_stack.shape[2], n2, n1]
1675 #msnc.quickImshow(center_stack[n1:n2,0,:], title=f'center stack theta={theta_start}',
1676 # path=self.output_folder, extent=extent, save_fig=self.save_plots,
1677 # save_only=self.save_plots_only, show_grid=True, aspect='auto')
1568 1678
1569 # Get cross sectional diameter in mm 1679 # Get cross sectional diameter in mm
1570 cross_sectional_dim = center_stack.shape[2]*eff_pixel_size 1680 cross_sectional_dim = center_stack.shape[2]*eff_pixel_size
1571 logging.debug(f'cross_sectional_dim = {cross_sectional_dim}') 1681 logging.debug(f'cross_sectional_dim = {cross_sectional_dim}')
1572 1682
1574 logging.info('Determine center offset at sample row boundaries') 1684 logging.info('Determine center offset at sample row boundaries')
1575 1685
1576 # Lower row center 1686 # Lower row center
1577 use_row = 'no' 1687 use_row = 'no'
1578 use_center = 'no' 1688 use_center = 'no'
1579 row = find_center.get('lower_row') 1689 row = center_rows[0]
1580 if msnc.is_int(row, n1, n2-2): 1690 if msnc.is_int(row, n1, n2-2):
1581 if self.test_mode: 1691 if self.test_mode or self.galaxy_flag:
1582 assert(row is not None) 1692 assert(row is not None)
1583 use_row = 'yes' 1693 use_row = 'yes'
1584 else: 1694 else:
1585 msnc.quickImshow(center_stack[:,0,:], title=f'theta={theta_start}', aspect='auto') 1695 msnc.quickImshow(center_stack[:,0,:], title=f'theta={theta_start}', aspect='auto')
1586 use_row = pyip.inputYesNo('\nCurrent row index for lower center = ' 1696 use_row = pyip.inputYesNo('\nCurrent row index for lower center = '
1603 row = n1 1713 row = n1
1604 if self.save_plots_only: 1714 if self.save_plots_only:
1605 msnc.clearFig(f'theta={theta_start}') 1715 msnc.clearFig(f'theta={theta_start}')
1606 # center_stack order: row,theta,column 1716 # center_stack order: row,theta,column
1607 center_offset = self._findCenterOnePlane(center_stack[row,:,:], row, thetas_deg, 1717 center_offset = self._findCenterOnePlane(center_stack[row,:,:], row, thetas_deg,
1608 eff_pixel_size, cross_sectional_dim) 1718 eff_pixel_size, cross_sectional_dim, num_core=num_core,
1609 logging.info(f'Lower center offset = {center_offset}') 1719 recon_pngname=recon_low_pngname)
1610 1720
1611 # Update config and save to file 1721 # Update config and save to file
1612 find_center['row_bounds'] = [n1, n2] 1722 find_center['row_bounds'] = [n1, n2]
1613 find_center['lower_row'] = row 1723 find_center['lower_row'] = row
1614 find_center['lower_center_offset'] = center_offset 1724 find_center['lower_center_offset'] = center_offset
1616 lower_row = row 1726 lower_row = row
1617 1727
1618 # Upper row center 1728 # Upper row center
1619 use_row = 'no' 1729 use_row = 'no'
1620 use_center = 'no' 1730 use_center = 'no'
1621 row = find_center.get('upper_row') 1731 row = center_rows[1]
1622 if msnc.is_int(row, lower_row+1, n2-1): 1732 if msnc.is_int(row, lower_row+1, n2-1):
1623 if self.test_mode: 1733 if self.test_mode or self.galaxy_flag:
1624 assert(row is not None) 1734 assert(row is not None)
1625 use_row = 'yes' 1735 use_row = 'yes'
1626 else: 1736 else:
1627 msnc.quickImshow(center_stack[:,0,:], title=f'theta={theta_start}', aspect='auto') 1737 msnc.quickImshow(center_stack[:,0,:], title=f'theta={theta_start}', aspect='auto')
1628 use_row = pyip.inputYesNo('\nCurrent row index for upper center = ' 1738 use_row = pyip.inputYesNo('\nCurrent row index for upper center = '
1645 row = n2-1 1755 row = n2-1
1646 if self.save_plots_only: 1756 if self.save_plots_only:
1647 msnc.clearFig(f'theta={theta_start}') 1757 msnc.clearFig(f'theta={theta_start}')
1648 # center_stack order: row,theta,column 1758 # center_stack order: row,theta,column
1649 center_offset = self._findCenterOnePlane(center_stack[row,:,:], row, thetas_deg, 1759 center_offset = self._findCenterOnePlane(center_stack[row,:,:], row, thetas_deg,
1650 eff_pixel_size, cross_sectional_dim) 1760 eff_pixel_size, cross_sectional_dim, num_core=num_core,
1761 recon_pngname=recon_upp_pngname)
1651 logging.info(f'upper_center_offset = {center_offset}') 1762 logging.info(f'upper_center_offset = {center_offset}')
1652 del center_stack 1763 del center_stack
1653 1764
1654 # Update config and save to file 1765 # Update config and save to file
1655 find_center['upper_row'] = row 1766 find_center['upper_row'] = row
1729 del plane1, plane2, plane1_shifted 1840 del plane1, plane2, plane1_shifted
1730 1841
1731 # Update config file 1842 # Update config file
1732 self.config = msnc.update('config.txt', 'check_centers', True, 'find_centers') 1843 self.config = msnc.update('config.txt', 'check_centers', True, 'find_centers')
1733 1844
1734 def reconstructTomoStacks(self): 1845 def reconstructTomoStacks(self, output_name=None, num_core=None):
1735 """Reconstruct tomography stacks. 1846 """Reconstruct tomography stacks.
1736 """ 1847 """
1848 if num_core is None:
1849 num_core = self.num_core
1737 logging.debug('Reconstruct tomography stacks') 1850 logging.debug('Reconstruct tomography stacks')
1738 stacks = self.config['stack_info']['stacks'] 1851 stacks = self.config['stack_info']['stacks']
1739 assert(len(self.tomo_stacks) == self.config['stack_info']['num']) 1852 assert(len(self.tomo_stacks) == self.config['stack_info']['num'])
1740 assert(len(self.tomo_stacks) == len(stacks)) 1853 assert(len(self.tomo_stacks) == len(stacks))
1741 assert(len(self.tomo_recon_stacks) == len(stacks)) 1854 assert(len(self.tomo_recon_stacks) == len(stacks))
1855 if self.galaxy_flag:
1856 assert(isinstance(output_name, str))
1857 else:
1858 if output_name:
1859 logging.warning('Ignoring output_name in reconstructTomoStacks '+
1860 '(only used in Galaxy)')
1742 1861
1743 # Get rotation axis rows and centers 1862 # Get rotation axis rows and centers
1744 find_center = self.config['find_center'] 1863 find_center = self.config['find_center']
1745 lower_row = find_center.get('lower_row') 1864 lower_row = find_center.get('lower_row')
1746 if lower_row is None: 1865 if lower_row is None:
1782 for i,stack in enumerate(stacks): 1901 for i,stack in enumerate(stacks):
1783 # Check if stack can be loaded 1902 # Check if stack can be loaded
1784 # reconstructed stack order for each one in stack : row/z,x,y 1903 # reconstructed stack order for each one in stack : row/z,x,y
1785 # preprocessed stack order for each one in stack: row,theta,column 1904 # preprocessed stack order for each one in stack: row,theta,column
1786 index = stack['index'] 1905 index = stack['index']
1787 available = False 1906 if not self.galaxy_flag:
1788 if stack.get('reconstructed', False): 1907 available = False
1789 self.tomo_recon_stacks[i], available = self._loadTomo('recon stack', index) 1908 if stack.get('reconstructed', False):
1790 if self.tomo_recon_stacks[i].size or available: 1909 self.tomo_recon_stacks[i], available = self._loadTomo('recon stack', index)
1791 if self.tomo_stacks[i].size: 1910 if self.tomo_recon_stacks[i].size or available:
1792 self.tomo_stacks[i] = np.array([]) 1911 if self.tomo_stacks[i].size:
1793 assert(stack.get('preprocessed', False) == True) 1912 self.tomo_stacks[i] = np.array([])
1794 assert(stack.get('reconstructed', False) == True) 1913 assert(stack.get('preprocessed', False) == True)
1914 assert(stack.get('reconstructed', False) == True)
1915 continue
1916 stack['reconstructed'] = False
1917 if not self.tomo_stacks[i].size:
1918 self.tomo_stacks[i], available = self._loadTomo('red stack', index,
1919 required=True)
1920 if not self.tomo_stacks[i].size:
1921 logging.error(f'Unable to load tomography stack {index} for reconstruction')
1922 stack[i]['preprocessed'] = False
1923 load_error = True
1795 continue 1924 continue
1796 else: 1925 assert(0 <= lower_row < upper_row < self.tomo_stacks[i].shape[0])
1797 stack['reconstructed'] = False 1926 center_offsets = [lower_center_offset-lower_row*center_slope,
1798 if not self.tomo_stacks[i].size: 1927 upper_center_offset+(self.tomo_stacks[i].shape[0]-1-upper_row)*center_slope]
1799 self.tomo_stacks[i], available = self._loadTomo('red stack', index, 1928 t0 = time()
1800 required=True) 1929 self.tomo_recon_stacks[i]= self._reconstructOneTomoStack(self.tomo_stacks[i],
1801 if not self.tomo_stacks[i].size: 1930 thetas, center_offsets=center_offsets, sigma=0.1, num_core=num_core,
1802 logging.error(f'Unable to load tomography stack {index} for reconstruction') 1931 algorithm='gridrec', run_secondary_sirt=True, secondary_iter=25)
1803 stack[i]['preprocessed'] = False 1932 logging.info(f'Reconstruction of stack {index} took {time()-t0:.2f} seconds!')
1804 load_error = True 1933 if not self.test_mode and not self.galaxy_flag:
1805 continue 1934 row_slice = int(self.tomo_stacks[i].shape[0]/2)
1806 assert(0 <= lower_row < upper_row < self.tomo_stacks[i].shape[0]) 1935 title = f'{basetitle} {index} slice{row_slice}'
1807 center_offsets = [lower_center_offset-lower_row*center_slope, 1936 msnc.quickImshow(self.tomo_recon_stacks[i][row_slice,:,:], title=title,
1808 upper_center_offset+(self.tomo_stacks[i].shape[0]-1-upper_row)*center_slope] 1937 path=self.output_folder, save_fig=self.save_plots,
1809 t0 = time() 1938 save_only=self.save_plots_only)
1810 self.tomo_recon_stacks[i]= self._reconstructOneTomoStack(self.tomo_stacks[i], 1939 msnc.quickPlot(self.tomo_recon_stacks[i]
1811 thetas, center_offsets=center_offsets, sigma=0.1, num_core=self.num_core, 1940 [row_slice,int(self.tomo_recon_stacks[i].shape[2]/2),:],
1812 algorithm='gridrec', run_secondary_sirt=True, secondary_iter=25) 1941 title=f'{title} cut{int(self.tomo_recon_stacks[i].shape[2]/2)}',
1813 logging.info(f'Reconstruction of stack {index} took {time()-t0:.2f} seconds!') 1942 path=self.output_folder, save_fig=self.save_plots,
1814 if not self.test_mode: 1943 save_only=self.save_plots_only)
1815 row_slice = int(self.tomo_stacks[i].shape[0]/2) 1944 self._saveTomo('recon stack', self.tomo_recon_stacks[i], index)
1816 title = f'{basetitle} {index} slice{row_slice}' 1945 # else:
1817 msnc.quickImshow(self.tomo_recon_stacks[i][row_slice,:,:], title=title, 1946 # np.savetxt(self.output_folder+f'recon_stack_{index}.txt',
1818 path=self.output_folder, save_fig=self.save_plots, 1947 # self.tomo_recon_stacks[i][row_slice,:,:], fmt='%.6e')
1819 save_only=self.save_plots_only) 1948 self.tomo_stacks[i] = np.array([])
1820 msnc.quickPlot(self.tomo_recon_stacks[i] 1949
1821 [row_slice,int(self.tomo_recon_stacks[i].shape[2]/2),:], 1950 # Update config and save to file
1822 title=f'{title} cut{int(self.tomo_recon_stacks[i].shape[2]/2)}', 1951 stack['reconstructed'] = True
1823 path=self.output_folder, save_fig=self.save_plots, 1952 combine_stacks = self.config.get('combine_stacks')
1824 save_only=self.save_plots_only) 1953 if combine_stacks and index in combine_stacks.get('stacks', []):
1825 self._saveTomo('recon stack', self.tomo_recon_stacks[i], index) 1954 combine_stacks['stacks'].remove(index)
1826 # else: 1955 self.cf.saveFile(self.config_out)
1827 # np.savetxt(self.output_folder+f'recon_stack_{index}.txt', 1956
1828 # self.tomo_recon_stacks[i][row_slice,:,:], fmt='%.6e') 1957 # Save reconstructed tomography stack to file
1829 self.tomo_stacks[i] = np.array([]) 1958 if self.galaxy_flag:
1830 1959 t0 = time()
1831 # Update config and save to file 1960 logging.info(f'Saving reconstructed tomography stack to {output_name} ...')
1832 stack['reconstructed'] = True 1961 save_stacks = {f'set_{stack["index"]}':tomo_stack
1833 combine_stacks = self.config.get('combine_stacks') 1962 for stack,tomo_stack in zip(stacks,self.tomo_recon_stacks)}
1834 if combine_stacks and index in combine_stacks.get('stacks', []): 1963 np.savez(output_name, **save_stacks)
1835 combine_stacks['stacks'].remove(index) 1964 logging.info(f'... done in {time()-t0:.2f} seconds!')
1836 self.cf.saveFile(self.config_out)
1837 1965
1838 def combineTomoStacks(self): 1966 def combineTomoStacks(self):
1839 """Combine the reconstructed tomography stacks. 1967 """Combine the reconstructed tomography stacks.
1840 """ 1968 """
1841 # stack order: stack,row(z),x,y 1969 # stack order: stack,row(z),x,y
2093 parser.add_argument('-t', '--test_mode', 2221 parser.add_argument('-t', '--test_mode',
2094 action='store_true', 2222 action='store_true',
2095 default=False, 2223 default=False,
2096 help='Test mode flag') 2224 help='Test mode flag')
2097 parser.add_argument('--num_core', 2225 parser.add_argument('--num_core',
2226 type=int,
2098 default=-1, 2227 default=-1,
2099 help='Number of cores') 2228 help='Number of cores')
2100 args = parser.parse_args() 2229 args = parser.parse_args()
2101 2230
2102 if args.config is None: 2231 if args.config is None: