Mercurial > repos > rv43 > tomo
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: |