|
1
|
1 #!/usr/bin/env python3
|
|
|
2
|
|
|
3 # -*- coding: utf-8 -*-
|
|
|
4 """
|
|
|
5 Created on Fri Dec 10 09:54:37 2021
|
|
|
6
|
|
|
7 @author: rv43
|
|
|
8 """
|
|
|
9
|
|
|
10 import logging
|
|
|
11
|
|
|
12 import os
|
|
|
13 import sys
|
|
|
14 import getopt
|
|
|
15 import re
|
|
|
16 import io
|
|
|
17 import pyinputplus as pyip
|
|
|
18 import numpy as np
|
|
|
19 import numexpr as ne
|
|
|
20 import multiprocessing as mp
|
|
|
21 import scipy.ndimage as spi
|
|
|
22 import tomopy
|
|
|
23 from time import time
|
|
|
24 from skimage.transform import iradon
|
|
|
25 from skimage.restoration import denoise_tv_chambolle
|
|
|
26
|
|
|
27 import msnc_tools as msnc
|
|
|
28
|
|
|
29 class set_numexpr_threads:
|
|
|
30
|
|
|
31 def __init__(self, nthreads):
|
|
|
32 cpu_count = mp.cpu_count()
|
|
|
33 if nthreads is None or nthreads > cpu_count:
|
|
|
34 self.n = cpu_count
|
|
|
35 else:
|
|
|
36 self.n = nthreads
|
|
|
37
|
|
|
38 def __enter__(self):
|
|
|
39 self.oldn = ne.set_num_threads(self.n)
|
|
|
40
|
|
|
41 def __exit__(self, exc_type, exc_value, traceback):
|
|
|
42 ne.set_num_threads(self.oldn)
|
|
|
43
|
|
|
44 class ConfigTomo(msnc.Config):
|
|
|
45 """Class for processing a config file.
|
|
|
46 """
|
|
|
47
|
|
|
48 def __init__(self, config_file=None, config_dict=None):
|
|
|
49 super().__init__(config_file, config_dict)
|
|
|
50
|
|
|
51 def _validate_txt(self):
|
|
|
52 """Returns False if any required config parameter is illegal or missing.
|
|
|
53 """
|
|
|
54 is_valid = True
|
|
|
55
|
|
|
56 # Check for required first-level keys
|
|
|
57 pars_required = ['tdf_data_path', 'tbf_data_path', 'detector_id']
|
|
|
58 pars_missing = []
|
|
|
59 is_valid = super().validate(pars_required, pars_missing)
|
|
|
60 if len(pars_missing) > 0:
|
|
|
61 logging.error(f'Missing item(s) in config file: {", ".join(pars_missing)}')
|
|
|
62 self.detector_id = self.config.get('detector_id')
|
|
|
63
|
|
|
64 # Find tomography dark field images file/folder
|
|
|
65 self.tdf_data_path = self.config.get('tdf_data_path')
|
|
|
66
|
|
|
67 # Find tomography bright field images file/folder
|
|
|
68 self.tbf_data_path = self.config.get('tbf_data_path')
|
|
|
69
|
|
|
70 # Check number of tomography image stacks
|
|
|
71 self.num_tomo_stacks = self.config.get('num_tomo_stacks', 1)
|
|
|
72 if not msnc.is_int(self.num_tomo_stacks, 1):
|
|
|
73 self.num_tomo_stacks = None
|
|
|
74 msnc.illegal_value('num_tomo_stacks', self.num_tomo_stacks, 'config file')
|
|
|
75 return False
|
|
|
76 logging.info(f'num_tomo_stacks = {self.num_tomo_stacks}')
|
|
|
77
|
|
|
78 # Find tomography images file/folders and stack parameters
|
|
|
79 tomo_data_paths_indices = sorted({key:value for key,value in self.config.items()
|
|
|
80 if 'tomo_data_path' in key}.items())
|
|
|
81 if len(tomo_data_paths_indices) != self.num_tomo_stacks:
|
|
|
82 logging.error(f'Incorrect number of tomography data path names in config file')
|
|
|
83 is_valid = False
|
|
|
84 self.tomo_data_paths = [tomo_data_paths_indices[i][1] for i in range(self.num_tomo_stacks)]
|
|
|
85 self.tomo_data_indices = [msnc.get_trailing_int(tomo_data_paths_indices[i][0])
|
|
|
86 if msnc.get_trailing_int(tomo_data_paths_indices[i][0]) else None
|
|
|
87 for i in range(self.num_tomo_stacks)]
|
|
|
88 tomo_ref_height_indices = sorted({key:value for key,value in self.config.items()
|
|
|
89 if 'z_pos' in key}.items())
|
|
|
90 if self.num_tomo_stacks > 1 and len(tomo_ref_height_indices) != self.num_tomo_stacks:
|
|
|
91 logging.error(f'Incorrect number of tomography reference heights in config file')
|
|
|
92 is_valid = False
|
|
|
93 if len(tomo_ref_height_indices):
|
|
|
94 self.tomo_ref_heights = [
|
|
|
95 tomo_ref_height_indices[i][1] for i in range(self.num_tomo_stacks)]
|
|
|
96 else:
|
|
|
97 self.tomo_ref_heights = [0.0]*self.num_tomo_stacks
|
|
|
98
|
|
|
99 # Check tomo angle (theta) range
|
|
|
100 self.start_theta = self.config.get('start_theta', 0.)
|
|
|
101 if not msnc.is_num(self.start_theta, 0.):
|
|
|
102 msnc.illegal_value('start_theta', self.start_theta, 'config file')
|
|
|
103 is_valid = False
|
|
|
104 logging.debug(f'start_theta = {self.start_theta}')
|
|
|
105 self.end_theta = self.config.get('end_theta', 180.)
|
|
|
106 if not msnc.is_num(self.end_theta, self.start_theta):
|
|
|
107 msnc.illegal_value('end_theta', self.end_theta, 'config file')
|
|
|
108 is_valid = False
|
|
|
109 logging.debug(f'end_theta = {self.end_theta}')
|
|
|
110 self.num_thetas = self.config.get('num_thetas')
|
|
|
111 if not (self.num_thetas is None or msnc.is_int(self.num_thetas, 1)):
|
|
|
112 msnc.illegal_value('num_thetas', self.num_thetas, 'config file')
|
|
|
113 self.num_thetas = None
|
|
|
114 is_valid = False
|
|
|
115 logging.debug(f'num_thetas = {self.num_thetas}')
|
|
|
116
|
|
|
117 return is_valid
|
|
|
118
|
|
|
119 def _validate_yaml(self):
|
|
|
120 """Returns False if any required config parameter is illegal or missing.
|
|
|
121 """
|
|
|
122 is_valid = True
|
|
|
123
|
|
|
124 # Check for required first-level keys
|
|
|
125 pars_required = ['dark_field', 'bright_field', 'stack_info', 'detector']
|
|
|
126 pars_missing = []
|
|
|
127 is_valid = super().validate(pars_required, pars_missing)
|
|
|
128 if len(pars_missing) > 0:
|
|
|
129 logging.error(f'Missing item(s) in config file: {", ".join(pars_missing)}')
|
|
|
130 self.detector_id = self.config['detector'].get('id')
|
|
|
131
|
|
|
132 # Find tomography dark field images file/folder
|
|
|
133 self.tdf_data_path = self.config['dark_field'].get('data_path')
|
|
|
134
|
|
|
135 # Find tomography bright field images file/folder
|
|
|
136 self.tbf_data_path = self.config['bright_field'].get('data_path')
|
|
|
137
|
|
|
138 # Check number of tomography image stacks
|
|
|
139 stack_info = self.config['stack_info']
|
|
|
140 self.num_tomo_stacks = stack_info.get('num', 1)
|
|
|
141 if not msnc.is_int(self.num_tomo_stacks, 1):
|
|
|
142 self.num_tomo_stacks = None
|
|
|
143 msnc.illegal_value('stack_info:num', self.num_tomo_stacks, 'config file')
|
|
|
144 return False
|
|
|
145 logging.info(f'num_tomo_stacks = {self.num_tomo_stacks}')
|
|
|
146
|
|
|
147 # Find tomography images file/folders and stack parameters
|
|
|
148 stacks = stack_info.get('stacks')
|
|
|
149 if stacks is None or len(stacks) is not self.num_tomo_stacks:
|
|
|
150 msnc.illegal_value('stack_info:stacks', stacks, 'config file')
|
|
|
151 return False
|
|
|
152 self.tomo_data_paths = []
|
|
|
153 self.tomo_data_indices = []
|
|
|
154 self.tomo_ref_heights = []
|
|
|
155 for stack in stacks:
|
|
|
156 self.tomo_data_paths.append(stack.get('data_path'))
|
|
|
157 self.tomo_data_indices.append(stack.get('index'))
|
|
|
158 self.tomo_ref_heights.append(stack.get('ref_height'))
|
|
|
159
|
|
|
160 # Check tomo angle (theta) range
|
|
|
161 theta_range = self.config.get('theta_range')
|
|
|
162 if theta_range is None:
|
|
|
163 self.start_theta = 0.
|
|
|
164 self.end_theta = 180.
|
|
|
165 self.num_thetas = None
|
|
|
166 else:
|
|
|
167 self.start_theta = theta_range.get('start', 0.)
|
|
|
168 if not msnc.is_num(self.start_theta, 0.):
|
|
|
169 msnc.illegal_value('theta_range:start', self.start_theta, 'config file')
|
|
|
170 is_valid = False
|
|
|
171 logging.debug(f'start_theta = {self.start_theta}')
|
|
|
172 self.end_theta = theta_range.get('end', 180.)
|
|
|
173 if not msnc.is_num(self.end_theta, self.start_theta):
|
|
|
174 msnc.illegal_value('theta_range:end', self.end_theta, 'config file')
|
|
|
175 is_valid = False
|
|
|
176 logging.debug(f'end_theta = {self.end_theta}')
|
|
|
177 self.num_thetas = theta_range.get('num')
|
|
|
178 if self.num_thetas and not msnc.is_int(self.num_thetas, 1):
|
|
|
179 msnc.illegal_value('theta_range:num', self.num_thetas, 'config file')
|
|
|
180 self.num_thetas = None
|
|
|
181 is_valid = False
|
|
|
182 logging.debug(f'num_thetas = {self.num_thetas}')
|
|
|
183
|
|
|
184 return is_valid
|
|
|
185
|
|
|
186 def validate(self):
|
|
|
187 """Returns False if any required config parameter is illegal or missing.
|
|
|
188 """
|
|
|
189 is_valid = True
|
|
|
190
|
|
|
191 # Check work_folder (shared by both file formats)
|
|
|
192 work_folder = os.path.abspath(self.config.get('work_folder', ''))
|
|
|
193 if not os.path.isdir(work_folder):
|
|
|
194 msnc.illegal_value('work_folder', work_folder, 'config file')
|
|
|
195 is_valid = False
|
|
|
196 logging.info(f'work_folder: {work_folder}')
|
|
|
197
|
|
|
198 # Check data filetype (shared by both file formats)
|
|
|
199 self.data_filetype = self.config.get('data_filetype', 'tif')
|
|
|
200 if not isinstance(self.data_filetype, str) or (self.data_filetype != 'tif' and
|
|
|
201 self.data_filetype != 'h5'):
|
|
|
202 msnc.illegal_value('data_filetype', self.data_filetype, 'config file')
|
|
|
203
|
|
|
204 if self.suffix == '.yml' or self.suffix == '.yaml':
|
|
|
205 is_valid = self._validate_yaml()
|
|
|
206 elif self.suffix == '.txt':
|
|
|
207 is_valid = self._validate_txt()
|
|
|
208 else:
|
|
|
209 logging.error(f'Undefined or illegal config file extension: {self.suffix}')
|
|
|
210
|
|
|
211 # Find tomography bright field images file/folder
|
|
|
212 if self.tdf_data_path:
|
|
|
213 if self.data_filetype == 'h5':
|
|
|
214 if isinstance(self.tdf_data_path, str):
|
|
|
215 if not os.path.isabs(self.tdf_data_path):
|
|
|
216 self.tdf_data_path = os.path.abspath(
|
|
|
217 f'{work_folder}/{self.tdf_data_path}')
|
|
|
218 else:
|
|
|
219 msnc.illegal_value('tdf_data_path', tdf_data_fil, 'config file')
|
|
|
220 is_valid = False
|
|
|
221 else:
|
|
|
222 if isinstance(self.tdf_data_path, int):
|
|
|
223 self.tdf_data_path = os.path.abspath(
|
|
|
224 f'{work_folder}/{self.tdf_data_path}/nf')
|
|
|
225 elif isinstance(self.tdf_data_path, str):
|
|
|
226 if not os.path.isabs(self.tdf_data_path):
|
|
|
227 self.tdf_data_path = os.path.abspath(
|
|
|
228 f'{work_folder}/{self.tdf_data_path}')
|
|
|
229 else:
|
|
|
230 msnc.illegal_value('tdf_data_path', self.tdf_data_path, 'config file')
|
|
|
231 is_valid = False
|
|
|
232 logging.info(f'dark field images path = {self.tdf_data_path}')
|
|
|
233
|
|
|
234 # Find tomography bright field images file/folder
|
|
|
235 if self.tbf_data_path:
|
|
|
236 if self.data_filetype == 'h5':
|
|
|
237 if isinstance(self.tbf_data_path, str):
|
|
|
238 if not os.path.isabs(self.tbf_data_path):
|
|
|
239 self.tbf_data_path = os.path.abspath(
|
|
|
240 f'{work_folder}/{self.tbf_data_path}')
|
|
|
241 else:
|
|
|
242 msnc.illegal_value('tbf_data_path', tbf_data_fil, 'config file')
|
|
|
243 is_valid = False
|
|
|
244 else:
|
|
|
245 if isinstance(self.tbf_data_path, int):
|
|
|
246 self.tbf_data_path = os.path.abspath(
|
|
|
247 f'{work_folder}/{self.tbf_data_path}/nf')
|
|
|
248 elif isinstance(self.tbf_data_path, str):
|
|
|
249 if not os.path.isabs(self.tbf_data_path):
|
|
|
250 self.tbf_data_path = os.path.abspath(
|
|
|
251 f'{work_folder}/{self.tbf_data_path}')
|
|
|
252 else:
|
|
|
253 msnc.illegal_value('tbf_data_path', self.tbf_data_path, 'config file')
|
|
|
254 is_valid = False
|
|
|
255 logging.info(f'bright field images path = {self.tbf_data_path}')
|
|
|
256
|
|
|
257 # Find tomography images file/folders and stack parameters
|
|
|
258 tomo_data_paths = []
|
|
|
259 tomo_data_indices = []
|
|
|
260 tomo_ref_heights = []
|
|
|
261 for data_path, index, ref_height in zip(self.tomo_data_paths, self.tomo_data_indices,
|
|
|
262 self.tomo_ref_heights):
|
|
|
263 if self.data_filetype == 'h5':
|
|
|
264 if isinstance(data_path, str):
|
|
|
265 if not os.path.isabs(data_path):
|
|
|
266 data_path = os.path.abspath(f'{work_folder}/{data_path}')
|
|
|
267 else:
|
|
|
268 msnc.illegal_value(f'stack_info:stacks:data_path', data_path, 'config file')
|
|
|
269 is_valid = False
|
|
|
270 data_path = None
|
|
|
271 else:
|
|
|
272 if isinstance(data_path, int):
|
|
|
273 data_path = os.path.abspath(f'{work_folder}/{data_path}/nf')
|
|
|
274 elif isinstance(data_path, str):
|
|
|
275 if not os.path.isabs(data_path):
|
|
|
276 data_path = os.path.abspath(f'{work_folder}/{data_path}')
|
|
|
277 else:
|
|
|
278 msnc.illegal_value(f'stack_info:stacks:data_path', data_path, 'config file')
|
|
|
279 is_valid = False
|
|
|
280 data_path = None
|
|
|
281 tomo_data_paths.append(data_path)
|
|
|
282 if index is None:
|
|
|
283 if self.num_tomo_stacks > 1:
|
|
|
284 logging.error('Missing stack_info:stacks:index in config file')
|
|
|
285 is_valid = False
|
|
|
286 index = None
|
|
|
287 else:
|
|
|
288 index = 1
|
|
|
289 elif not isinstance(index, int):
|
|
|
290 msnc.illegal_value(f'stack_info:stacks:index', index, 'config file')
|
|
|
291 is_valid = False
|
|
|
292 index = None
|
|
|
293 tomo_data_indices.append(index)
|
|
|
294 if ref_height is None:
|
|
|
295 if self.num_tomo_stacks > 1:
|
|
|
296 logging.error('Missing stack_info:stacks:ref_height in config file')
|
|
|
297 is_valid = False
|
|
|
298 ref_height = None
|
|
|
299 else:
|
|
|
300 ref_height = 0.
|
|
|
301 elif not msnc.is_num(ref_height):
|
|
|
302 msnc.illegal_value(f'stack_info:stacks:ref_height', ref_height, 'config file')
|
|
|
303 is_valid = False
|
|
|
304 ref_height = None
|
|
|
305 # Set reference heights relative to first stack
|
|
|
306 if (len(tomo_ref_heights) and msnc.is_num(ref_height) and
|
|
|
307 msnc.is_num(tomo_ref_heights[0])):
|
|
|
308 ref_height = (round(ref_height-tomo_ref_heights[0], 3))
|
|
|
309 tomo_ref_heights.append(ref_height)
|
|
|
310 tomo_ref_heights[0] = 0.0
|
|
|
311 logging.info('tomography data paths:')
|
|
|
312 for i in range(self.num_tomo_stacks):
|
|
|
313 logging.info(f' {tomo_data_paths[i]}')
|
|
|
314 logging.info(f'tomography data path indices: {tomo_data_indices}')
|
|
|
315 logging.info(f'tomography reference heights: {tomo_ref_heights}')
|
|
|
316
|
|
|
317 # Update config in memory
|
|
|
318 if self.suffix == '.txt':
|
|
|
319 self.config = {}
|
|
|
320 dark_field = self.config.get('dark_field')
|
|
|
321 if dark_field is None:
|
|
|
322 self.config['dark_field'] = {'data_path' : self.tdf_data_path}
|
|
|
323 else:
|
|
|
324 self.config['dark_field']['data_path'] = self.tdf_data_path
|
|
|
325 bright_field = self.config.get('bright_field')
|
|
|
326 if bright_field is None:
|
|
|
327 self.config['bright_field'] = {'data_path' : self.tbf_data_path}
|
|
|
328 else:
|
|
|
329 self.config['bright_field']['data_path'] = self.tbf_data_path
|
|
|
330 detector = self.config.get('detector')
|
|
|
331 if detector is None:
|
|
|
332 self.config['detector'] = {'id' : self.detector_id}
|
|
|
333 else:
|
|
|
334 detector['id'] = self.detector_id
|
|
|
335 self.config['work_folder'] = work_folder
|
|
|
336 self.config['data_filetype'] = self.data_filetype
|
|
|
337 stack_info = self.config.get('stack_info')
|
|
|
338 if stack_info is None:
|
|
|
339 stacks = []
|
|
|
340 for i in range(self.num_tomo_stacks):
|
|
|
341 stacks.append({'data_path' : tomo_data_paths[i], 'index' : tomo_data_indices[i],
|
|
|
342 'ref_height' : tomo_ref_heights[i]})
|
|
|
343 self.config['stack_info'] = {'num' : self.num_tomo_stacks, 'stacks' : stacks}
|
|
|
344 else:
|
|
|
345 stack_info['num'] = self.num_tomo_stacks
|
|
|
346 stacks = stack_info.get('stacks')
|
|
|
347 for i,stack in enumerate(stacks):
|
|
|
348 stack['data_path'] = tomo_data_paths[i]
|
|
|
349 stack['index'] = tomo_data_indices[i]
|
|
|
350 stack['ref_height'] = tomo_ref_heights[i]
|
|
|
351 if self.num_thetas:
|
|
|
352 theta_range = {'start' : self.start_theta, 'end' : self.end_theta,
|
|
|
353 'num' : self.num_thetas}
|
|
|
354 else:
|
|
|
355 theta_range = {'start' : self.start_theta, 'end' : self.end_theta}
|
|
|
356 self.config['theta_range'] = theta_range
|
|
|
357
|
|
|
358 # Cleanup temporary validation variables
|
|
|
359 del self.tdf_data_path
|
|
|
360 del self.tbf_data_path
|
|
|
361 del self.detector_id
|
|
|
362 del self.data_filetype
|
|
|
363 del self.num_tomo_stacks
|
|
|
364 del self.tomo_data_paths
|
|
|
365 del self.tomo_data_indices
|
|
|
366 del self.tomo_ref_heights
|
|
|
367 del self.start_theta
|
|
|
368 del self.end_theta
|
|
|
369 del self.num_thetas
|
|
|
370
|
|
|
371 return is_valid
|
|
|
372
|
|
|
373 class Tomo:
|
|
|
374 """Processing tomography data with misalignment.
|
|
|
375 """
|
|
|
376
|
|
|
377 def __init__(self, config_file=None, config_dict=None, config_out=None, output_folder='.',
|
|
|
378 log_level='INFO', log_stream='tomo.log', galaxy_flag=False, test_mode=False):
|
|
|
379 """Initialize with optional config input file or dictionary
|
|
|
380 """
|
|
|
381 self.ncore = mp.cpu_count()
|
|
|
382 self.config_out = config_out
|
|
|
383 self.output_folder = output_folder
|
|
|
384 self.galaxy_flag = galaxy_flag
|
|
|
385 self.test_mode = test_mode
|
|
|
386 self.save_plots = True # Make input argument?
|
|
|
387 self.save_plots_only = True # Make input argument?
|
|
|
388 self.cf = None
|
|
|
389 self.config = None
|
|
|
390 self.is_valid = True
|
|
|
391 self.tdf = np.array([])
|
|
|
392 self.tbf = np.array([])
|
|
|
393 self.tomo_stacks = []
|
|
|
394 self.tomo_recon_stacks = []
|
|
|
395
|
|
|
396 # Set log configuration
|
|
|
397 logging_format = '%(asctime)s : %(levelname)s - %(module)s : %(funcName)s - %(message)s'
|
|
|
398 if self.test_mode:
|
|
|
399 self.save_plots_only = True
|
|
|
400 if isinstance(log_stream, str):
|
|
|
401 logging.basicConfig(filename=f'{log_stream}', filemode='w',
|
|
|
402 format=logging_format, level=logging.WARNING, force=True)
|
|
|
403 elif isinstance(log_stream, io.TextIOWrapper):
|
|
|
404 logging.basicConfig(filemode='w', format=logging_format, level=logging.WARNING,
|
|
|
405 stream=log_stream, force=True)
|
|
|
406 else:
|
|
|
407 raise ValueError(f'Invalid log_stream: {log_stream}')
|
|
|
408 logging.warning('Ignoring log_level argument in test mode')
|
|
|
409 else:
|
|
|
410 level = getattr(logging, log_level.upper(), None)
|
|
|
411 if not isinstance(level, int):
|
|
|
412 raise ValueError(f'Invalid log_level: {log_level}')
|
|
|
413 if log_stream is sys.stdout:
|
|
|
414 logging.basicConfig(format=logging_format, level=level, force=True,
|
|
|
415 handlers=[logging.StreamHandler()])
|
|
|
416 else:
|
|
|
417 if isinstance(log_stream, str):
|
|
|
418 logging.basicConfig(filename=f'{log_stream}', filemode='w',
|
|
|
419 format=logging_format, level=level, force=True)
|
|
|
420 elif isinstance(log_stream, io.TextIOWrapper):
|
|
|
421 logging.basicConfig(filemode='w', format=logging_format, level=level,
|
|
|
422 stream=log_stream, force=True)
|
|
|
423 else:
|
|
|
424 raise ValueError(f'Invalid log_stream: {log_stream}')
|
|
|
425 stream_handler = logging.StreamHandler()
|
|
|
426 logging.getLogger().addHandler(stream_handler)
|
|
|
427 stream_handler.setLevel(logging.WARNING)
|
|
|
428 stream_handler.setFormatter(logging.Formatter(logging_format))
|
|
|
429
|
|
|
430 # Set output config file name
|
|
|
431 if self.config_out is None:
|
|
|
432 if self.config is None:
|
|
|
433 self.config_out = 'config.yaml'
|
|
|
434 else:
|
|
|
435 self.config_out = config_file
|
|
|
436
|
|
|
437 logging.info(f'ncore = {self.ncore}')
|
|
|
438 logging.debug(f'config_file = {config_file}')
|
|
|
439 logging.debug(f'config_dict = {config_dict}')
|
|
|
440 logging.debug(f'config_out = {self.config_out}')
|
|
|
441 logging.debug(f'output_folder = {self.output_folder}')
|
|
|
442 logging.debug(f'log_stream = {log_stream}')
|
|
|
443 logging.debug(f'log_level = {log_level}')
|
|
|
444 logging.debug(f'galaxy_flag = {self.galaxy_flag}')
|
|
|
445 logging.debug(f'test_mode = {self.test_mode}')
|
|
|
446
|
|
|
447 # Create config object and load config file
|
|
|
448 self.cf = ConfigTomo(config_file, config_dict)
|
|
|
449 if not self.cf.load_flag:
|
|
|
450 self.is_valid = False
|
|
|
451 return
|
|
|
452
|
|
|
453 if self.galaxy_flag:
|
|
|
454 self.ncore = 1 #RV can I set this? mp.cpu_count()
|
|
|
455 assert(self.output_folder == '.')
|
|
|
456 assert(self.test_mode is False)
|
|
|
457 self.save_plots = True
|
|
|
458 self.save_plots_only = True
|
|
|
459 else:
|
|
|
460 # Input validation is already performed during link_data_to_galaxy
|
|
|
461
|
|
|
462 # Check config file parameters
|
|
|
463 self.is_valid = self.cf.validate()
|
|
|
464
|
|
|
465 # Load detector info file
|
|
|
466 df = msnc.Detector(self.cf.config['detector']['id'])
|
|
|
467
|
|
|
468 # Check detector info file parameters
|
|
|
469 if df.validate():
|
|
|
470 pixel_size = df.getPixelSize()
|
|
|
471 num_rows, num_columns = df.getDimensions()
|
|
|
472 if not pixel_size or not num_rows or not num_columns:
|
|
|
473 self.is_valid = False
|
|
|
474 else:
|
|
|
475 pixel_size = None
|
|
|
476 num_rows = None
|
|
|
477 num_columns = None
|
|
|
478 self.is_valid = False
|
|
|
479
|
|
|
480 # Update config
|
|
|
481 self.cf.config['detector']['pixel_size'] = pixel_size
|
|
|
482 self.cf.config['detector']['rows'] = num_rows
|
|
|
483 self.cf.config['detector']['columns'] = num_columns
|
|
|
484 logging.debug(f'pixel_size = self.cf.config["detector"]["pixel_size"]')
|
|
|
485 logging.debug(f'num_rows: {self.cf.config["detector"]["rows"]}')
|
|
|
486 logging.debug(f'num_columns: {self.cf.config["detector"]["columns"]}')
|
|
|
487
|
|
|
488 # Safe config to file
|
|
|
489 if self.is_valid:
|
|
|
490 self.cf.saveFile(self.config_out)
|
|
|
491
|
|
|
492 # Initialize shortcut to config
|
|
|
493 self.config = self.cf.config
|
|
|
494
|
|
|
495 # Initialize tomography stack
|
|
|
496 num_tomo_stacks = self.config['stack_info']['num']
|
|
|
497 if num_tomo_stacks:
|
|
|
498 self.tomo_stacks = [np.array([]) for _ in range(num_tomo_stacks)]
|
|
|
499 self.tomo_recon_stacks = [np.array([]) for _ in range(num_tomo_stacks)]
|
|
|
500
|
|
|
501 logging.debug(f'save_plots = {self.save_plots}')
|
|
|
502 logging.debug(f'save_plots_only = {self.save_plots_only}')
|
|
|
503
|
|
|
504 def findImageFiles(self):
|
|
|
505 """Find all available image files.
|
|
|
506 """
|
|
|
507 self.is_valid = True
|
|
|
508
|
|
|
509 # Find dark field images
|
|
|
510 dark_field = self.config['dark_field']
|
|
|
511 img_start, num_imgs, dark_files = msnc.findImageFiles(
|
|
|
512 dark_field['data_path'], self.config['data_filetype'], 'dark field')
|
|
|
513 if img_start < 0 or num_imgs < 1:
|
|
|
514 logging.error('Unable to find suitable dark field images')
|
|
|
515 if dark_field['data_path']:
|
|
|
516 self.is_valid = False
|
|
|
517 dark_field['num'] = num_imgs
|
|
|
518 dark_field['img_start'] = img_start
|
|
|
519 logging.info(f'Number of dark field images = {dark_field["num"]}')
|
|
|
520 logging.info(f'Dark field image start index = {dark_field["img_start"]}')
|
|
|
521
|
|
|
522 # Find bright field images
|
|
|
523 bright_field = self.config['bright_field']
|
|
|
524 img_start, num_imgs, bright_files = msnc.findImageFiles(
|
|
|
525 bright_field['data_path'], self.config['data_filetype'], 'bright field')
|
|
|
526 if img_start < 0 or num_imgs < 1:
|
|
|
527 logging.error('Unable to find suitable bright field images')
|
|
|
528 self.is_valid = False
|
|
|
529 bright_field['num'] = num_imgs
|
|
|
530 bright_field['img_start'] = img_start
|
|
|
531 logging.info(f'Number of bright field images = {bright_field["num"]}')
|
|
|
532 logging.info(f'Bright field image start index = {bright_field["img_start"]}')
|
|
|
533
|
|
|
534 # Find tomography images
|
|
|
535 tomo_stack_files = []
|
|
|
536 for stack in self.config['stack_info']['stacks']:
|
|
|
537 index = stack['index']
|
|
|
538 img_start, num_imgs, tomo_files = msnc.findImageFiles(
|
|
|
539 stack['data_path'], self.config['data_filetype'], f'tomography set {index}')
|
|
|
540 if img_start < 0 or num_imgs < 1:
|
|
|
541 logging.error('Unable to find suitable tomography images')
|
|
|
542 self.is_valid = False
|
|
|
543 stack['num'] = num_imgs
|
|
|
544 stack['img_start'] = img_start
|
|
|
545 logging.info(f'Number of tomography images for set {index} = {stack["num"]}')
|
|
|
546 logging.info(f'Tomography set {index} image start index = {stack["img_start"]}')
|
|
|
547 tomo_stack_files.append(tomo_files)
|
|
|
548 del tomo_files
|
|
|
549
|
|
|
550 # Safe updated config
|
|
|
551 if self.is_valid:
|
|
|
552 self.cf.saveFile(self.config_out)
|
|
|
553
|
|
|
554 return dark_files, bright_files, tomo_stack_files
|
|
|
555
|
|
|
556 def selectImageRanges(self, available_stacks=None):
|
|
|
557 """Find and check all required image files.
|
|
|
558 """
|
|
|
559 self.is_valid = True
|
|
|
560 stack_info = self.config['stack_info']
|
|
|
561 if available_stacks is None:
|
|
|
562 available_stacks = [False]*stack_info['num']
|
|
|
563 elif len(available_stacks) != stack_info['num']:
|
|
|
564 logging.warning('Illegal dimension of available_stacks in getImageFiles '+
|
|
|
565 f'({len(available_stacks)}');
|
|
|
566 available_stacks = [False]*stack_info['num']
|
|
|
567
|
|
|
568 # Check number of tomography angles/thetas
|
|
|
569 num_thetas = self.config['theta_range'].get('num')
|
|
|
570 if num_thetas is None:
|
|
|
571 num_thetas = pyip.inputInt('\nEnter the number of thetas (>0): ', greaterThan=0)
|
|
|
572 elif not msnc.is_int(num_thetas, 0):
|
|
|
573 msnc.illegal_value('num_thetas', num_thetas, 'config file')
|
|
|
574 self.is_valid = False
|
|
|
575 return
|
|
|
576 self.config['theta_range']['num'] = num_thetas
|
|
|
577 logging.debug(f'num_thetas = {self.config["theta_range"]["num"]}')
|
|
|
578
|
|
|
579 # Find tomography dark field images
|
|
|
580 dark_field = self.config['dark_field']
|
|
|
581 img_start = dark_field.get('img_start', -1)
|
|
|
582 img_offset = dark_field.get('img_offset', -1)
|
|
|
583 num_imgs = dark_field.get('num', 0)
|
|
|
584 if not self.test_mode:
|
|
|
585 img_start, img_offset, num_imgs = msnc.selectImageRange(img_start, img_offset,
|
|
|
586 num_imgs, 'dark field')
|
|
|
587 if img_start < 0 or num_imgs < 1:
|
|
|
588 logging.error('Unable to find suitable dark field images')
|
|
|
589 if dark_field['data_path']:
|
|
|
590 self.is_valid = False
|
|
|
591 dark_field['img_start'] = img_start
|
|
|
592 dark_field['img_offset'] = img_offset
|
|
|
593 dark_field['num'] = num_imgs
|
|
|
594 logging.debug(f'Dark field image start index: {dark_field["img_start"]}')
|
|
|
595 logging.debug(f'Dark field image offset: {dark_field["img_offset"]}')
|
|
|
596 logging.debug(f'Number of dark field images: {dark_field["num"]}')
|
|
|
597
|
|
|
598 # Find tomography bright field images
|
|
|
599 bright_field = self.config['bright_field']
|
|
|
600 img_start = bright_field.get('img_start', -1)
|
|
|
601 img_offset = bright_field.get('img_offset', -1)
|
|
|
602 num_imgs = bright_field.get('num', 0)
|
|
|
603 if not self.test_mode:
|
|
|
604 img_start, img_offset, num_imgs = msnc.selectImageRange(img_start, img_offset,
|
|
|
605 num_imgs, 'bright field')
|
|
|
606 if img_start < 0 or num_imgs < 1:
|
|
|
607 logging.error('Unable to find suitable bright field images')
|
|
|
608 if bright_field['data_path']:
|
|
|
609 self.is_valid = False
|
|
|
610 bright_field['img_start'] = img_start
|
|
|
611 bright_field['img_offset'] = img_offset
|
|
|
612 bright_field['num'] = num_imgs
|
|
|
613 logging.debug(f'Bright field image start index: {bright_field["img_start"]}')
|
|
|
614 logging.debug(f'Bright field image offset: {bright_field["img_offset"]}')
|
|
|
615 logging.debug(f'Number of bright field images: {bright_field["num"]}')
|
|
|
616
|
|
|
617 # Find tomography images
|
|
|
618 for i,stack in enumerate(stack_info['stacks']):
|
|
|
619 # Check if stack is already loaded or available
|
|
|
620 if self.tomo_stacks[i].size or available_stacks[i]:
|
|
|
621 continue
|
|
|
622 index = stack['index']
|
|
|
623 img_start = stack.get('img_start', -1)
|
|
|
624 img_offset = stack.get('img_offset', -1)
|
|
|
625 num_imgs = num_thetas
|
|
|
626 if not self.test_mode:
|
|
|
627 img_start, img_offset, num_imgs = msnc.selectImageRange(img_start, img_offset,
|
|
|
628 num_imgs, f'tomography stack {index}', num_thetas)
|
|
|
629 if img_start < 0 or num_imgs != num_thetas:
|
|
|
630 logging.error('Unable to find suitable tomography images')
|
|
|
631 self.is_valid = False
|
|
|
632 stack['img_start'] = img_start
|
|
|
633 stack['img_offset'] = img_offset
|
|
|
634 stack['num'] = num_imgs
|
|
|
635 logging.debug(f'Tomography stack {index} image start index: {stack["img_start"]}')
|
|
|
636 logging.debug(f'Tomography stack {index} image offset: {stack["img_offset"]}')
|
|
|
637 logging.debug(f'Number of tomography images for stack {index}: {stack["num"]}')
|
|
|
638 available_stacks[i] = True
|
|
|
639
|
|
|
640 # Safe updated config to file
|
|
|
641 if self.is_valid:
|
|
|
642 self.cf.saveFile(self.config_out)
|
|
|
643
|
|
|
644 return
|
|
|
645
|
|
|
646 def _genDark(self, tdf_files, dark_field_pngname):
|
|
|
647 """Generate dark field.
|
|
|
648 """
|
|
|
649 # Load the dark field images
|
|
|
650 logging.debug('Loading dark field...')
|
|
|
651 dark_field = self.config['dark_field']
|
|
|
652 tdf_stack = msnc.loadImageStack(tdf_files, self.config['data_filetype'],
|
|
|
653 dark_field['img_offset'], dark_field['num'])
|
|
|
654
|
|
|
655 # Take median
|
|
|
656 self.tdf = np.median(tdf_stack, axis=0)
|
|
|
657 del tdf_stack
|
|
|
658
|
|
|
659 # RV make input of some kind (not always needed)
|
|
|
660 tdf_cutoff = 21
|
|
|
661 self.tdf[self.tdf > tdf_cutoff] = np.nan
|
|
|
662 tdf_mean = np.nanmean(self.tdf)
|
|
|
663 logging.debug(f'tdf_cutoff = {tdf_cutoff}')
|
|
|
664 logging.debug(f'tdf_mean = {tdf_mean}')
|
|
|
665 np.nan_to_num(self.tdf, copy=False, nan=tdf_mean, posinf=tdf_mean, neginf=0.)
|
|
|
666 if not self.test_mode and not self.galaxy_flag:
|
|
|
667 msnc.quickImshow(self.tdf, title='dark field', path=self.output_folder,
|
|
|
668 save_fig=self.save_plots, save_only=self.save_plots_only)
|
|
|
669 elif self.galaxy_flag:
|
|
|
670 msnc.quickImshow(self.tdf, title='dark field', name=dark_field_pngname,
|
|
|
671 save_fig=True, save_only=True)
|
|
|
672
|
|
|
673 def _genBright(self, tbf_files, bright_field_pngname):
|
|
|
674 """Generate bright field.
|
|
|
675 """
|
|
|
676 # Load the bright field images
|
|
|
677 logging.debug('Loading bright field...')
|
|
|
678 bright_field = self.config['bright_field']
|
|
|
679 tbf_stack = msnc.loadImageStack(tbf_files, self.config['data_filetype'],
|
|
|
680 bright_field['img_offset'], bright_field['num'])
|
|
|
681
|
|
|
682 # Take median
|
|
|
683 """Median or mean: It may be best to try the median because of some image
|
|
|
684 artifacts that arise due to crinkles in the upstream kapton tape windows
|
|
|
685 causing some phase contrast images to appear on the detector.
|
|
|
686 One thing that also may be useful in a future implementation is to do a
|
|
|
687 brightfield adjustment on EACH frame of the tomo based on a ROI in the
|
|
|
688 corner of the frame where there is no sample but there is the direct X-ray
|
|
|
689 beam because there is frame to frame fluctuations from the incoming beam.
|
|
|
690 We don’t typically account for them but potentially could.
|
|
|
691 """
|
|
|
692 self.tbf = np.median(tbf_stack, axis=0)
|
|
|
693 del tbf_stack
|
|
|
694
|
|
|
695 # Subtract dark field
|
|
|
696 if self.tdf.size:
|
|
|
697 self.tbf -= self.tdf
|
|
|
698 else:
|
|
|
699 logging.warning('Dark field unavailable')
|
|
|
700 if not self.test_mode and not self.galaxy_flag:
|
|
|
701 msnc.quickImshow(self.tbf, title='bright field', path=self.output_folder,
|
|
|
702 save_fig=self.save_plots, save_only=self.save_plots_only)
|
|
|
703 elif self.galaxy_flag:
|
|
|
704 msnc.quickImshow(self.tbf, title='bright field', name=bright_field_pngname,
|
|
|
705 save_fig=True, save_only=True)
|
|
|
706
|
|
|
707 def _setDetectorBounds(self, tomo_stack_files, tomo_field_pngname, detectorbounds_pngname):
|
|
|
708 """Set vertical detector bounds for image stack.
|
|
|
709 """
|
|
|
710 preprocess = self.config.get('preprocess')
|
|
|
711 if preprocess is None:
|
|
|
712 img_x_bounds = [0, self.tbf.shape[0]]
|
|
|
713 else:
|
|
|
714 img_x_bounds = preprocess.get('img_x_bounds', [0, self.tbf.shape[0]])
|
|
|
715 if self.test_mode:
|
|
|
716 # Update config and save to file
|
|
|
717 if preprocess is None:
|
|
|
718 self.cf.config['preprocess'] = {'img_x_bounds' : img_x_bounds}
|
|
|
719 else:
|
|
|
720 preprocess['img_x_bounds'] = img_x_bounds
|
|
|
721 self.cf.saveFile(self.config_out)
|
|
|
722 return
|
|
|
723
|
|
|
724 # Check reference heights
|
|
|
725 pixel_size = self.config['detector']['pixel_size']
|
|
|
726 if pixel_size is None:
|
|
|
727 raise ValueError('Detector pixel size unavailable')
|
|
|
728 if not self.tbf.size:
|
|
|
729 raise ValueError('Bright field unavailable')
|
|
|
730 num_x_min = None
|
|
|
731 num_tomo_stacks = self.config['stack_info']['num']
|
|
|
732 stacks = self.config['stack_info']['stacks']
|
|
|
733 if num_tomo_stacks > 1:
|
|
|
734 delta_z = stacks[1]['ref_height']-stacks[0]['ref_height']
|
|
|
735 for i in range(2, num_tomo_stacks):
|
|
|
736 delta_z = min(delta_z, stacks[i]['ref_height']-stacks[i-1]['ref_height'])
|
|
|
737 logging.debug(f'delta_z = {delta_z}')
|
|
|
738 num_x_min = int(delta_z/pixel_size)+1
|
|
|
739 logging.debug(f'num_x_min = {num_x_min}')
|
|
|
740 if num_x_min > self.tbf.shape[0]:
|
|
|
741 logging.warning('Image bounds and pixel size prevent seamless stacking')
|
|
|
742 num_x_min = self.tbf.shape[0]
|
|
|
743
|
|
|
744 # Select image bounds
|
|
|
745 if self.galaxy_flag:
|
|
|
746 if num_x_min is None or num_x_min >= self.tbf.shape[0]:
|
|
|
747 img_x_bounds = [0, self.tbf.shape[0]]
|
|
|
748 else:
|
|
|
749 margin = int(num_x_min/10)
|
|
|
750 offset = min(0, int((self.tbf.shape[0]-num_x_min)/2-margin))
|
|
|
751 img_x_bounds = [offset, num_x_min+offset+2*margin]
|
|
|
752 tomo_stack = msnc.loadImageStack(tomo_stack_files[0], self.config['data_filetype'],
|
|
|
753 stacks[0]['img_offset'], 1)
|
|
|
754 x_sum = np.sum(tomo_stack[0,:,:], 1)
|
|
|
755 title = f'tomography image at theta={self.config["theta_range"]["start"]}'
|
|
|
756 msnc.quickImshow(tomo_stack[0,:,:], title=title, name=tomo_field_pngname,
|
|
|
757 save_fig=True, save_only=True)
|
|
|
758 msnc.quickPlot((range(x_sum.size), x_sum),
|
|
|
759 ([img_x_bounds[0], img_x_bounds[0]], [x_sum.min(), x_sum.max()], 'r-'),
|
|
|
760 ([img_x_bounds[1], img_x_bounds[1]], [x_sum.min(), x_sum.max()], 'r-'),
|
|
|
761 title='sum over theta and y', name=detectorbounds_pngname,
|
|
|
762 save_fig=True, save_only=True)
|
|
|
763
|
|
|
764 # Update config and save to file
|
|
|
765 if preprocess is None:
|
|
|
766 self.cf.config['preprocess'] = {'img_x_bounds' : img_x_bounds}
|
|
|
767 else:
|
|
|
768 preprocess['img_x_bounds'] = img_x_bounds
|
|
|
769 self.cf.saveFile(self.config_out)
|
|
|
770 return
|
|
|
771
|
|
|
772 # For one tomography stack only: load the first image
|
|
|
773 title = None
|
|
|
774 msnc.quickImshow(self.tbf, title='bright field')
|
|
|
775 if num_tomo_stacks == 1:
|
|
|
776 tomo_stack = msnc.loadImageStack(tomo_stack_files[0], self.config['data_filetype'],
|
|
|
777 stacks[0]['img_offset'], 1)
|
|
|
778 title = f'tomography image at theta={self.config["theta_range"]["start"]}'
|
|
|
779 msnc.quickImshow(tomo_stack[0,:,:], title=title)
|
|
|
780 tomo_or_bright = pyip.inputNum('\nSelect image bounds from bright field (0) or '+
|
|
|
781 'from first tomography image (1): ', min=0, max=1)
|
|
|
782 else:
|
|
|
783 print('\nSelect image bounds from bright field')
|
|
|
784 tomo_or_bright = 0
|
|
|
785 if tomo_or_bright:
|
|
|
786 x_sum = np.sum(tomo_stack[0,:,:], 1)
|
|
|
787 use_bounds = 'no'
|
|
|
788 if img_x_bounds[0] is not None and img_x_bounds[1] is not None:
|
|
|
789 if img_x_bounds[0] < 0:
|
|
|
790 msnc.illegal_value('img_x_bounds[0]', img_x_bounds[0], 'config file')
|
|
|
791 img_x_bounds[0] = 0
|
|
|
792 if not img_x_bounds[0] < img_x_bounds[1] <= x_sum.size:
|
|
|
793 msnc.illegal_value('img_x_bounds[1]', img_x_bounds[1], 'config file')
|
|
|
794 img_x_bounds[1] = x_sum.size
|
|
|
795 tomo_tmp = tomo_stack[0,:,:]
|
|
|
796 tomo_tmp[img_x_bounds[0],:] = tomo_stack[0,:,:].max()
|
|
|
797 tomo_tmp[img_x_bounds[1],:] = tomo_stack[0,:,:].max()
|
|
|
798 title = f'tomography image at theta={self.config["theta_range"]["start"]}'
|
|
|
799 msnc.quickImshow(tomo_stack[0,:,:], title=title)
|
|
|
800 msnc.quickPlot((range(x_sum.size), x_sum),
|
|
|
801 ([img_x_bounds[0], img_x_bounds[0]], [x_sum.min(), x_sum.max()], 'r-'),
|
|
|
802 ([img_x_bounds[1], img_x_bounds[1]], [x_sum.min(), x_sum.max()], 'r-'),
|
|
|
803 title='sum over theta and y')
|
|
|
804 print(f'lower bound = {img_x_bounds[0]} (inclusive)\n'+
|
|
|
805 f'upper bound = {img_x_bounds[1]} (exclusive)]')
|
|
|
806 use_bounds = pyip.inputYesNo('Accept these bounds ([y]/n)?: ', blank=True)
|
|
|
807 if use_bounds == 'no':
|
|
|
808 img_x_bounds = msnc.selectImageBounds(tomo_stack[0,:,:], 0,
|
|
|
809 img_x_bounds[0], img_x_bounds[1], num_x_min,
|
|
|
810 f'tomography image at theta={self.config["theta_range"]["start"]}')
|
|
|
811 if num_x_min is not None and img_x_bounds[1]-img_x_bounds[0]+1 < num_x_min:
|
|
|
812 logging.warning('Image bounds and pixel size prevent seamless stacking')
|
|
|
813 tomo_tmp = tomo_stack[0,:,:]
|
|
|
814 tomo_tmp[img_x_bounds[0],:] = tomo_stack[0,:,:].max()
|
|
|
815 tomo_tmp[img_x_bounds[1],:] = tomo_stack[0,:,:].max()
|
|
|
816 title = f'tomography image at theta={self.config["theta_range"]["start"]}'
|
|
|
817 msnc.quickImshow(tomo_stack[0,:,:], title=title, path=self.output_folder,
|
|
|
818 save_fig=self.save_plots, save_only=True)
|
|
|
819 msnc.quickPlot(range(img_x_bounds[0], img_x_bounds[1]),
|
|
|
820 x_sum[img_x_bounds[0]:img_x_bounds[1]],
|
|
|
821 title='sum over theta and y', path=self.output_folder,
|
|
|
822 save_fig=self.save_plots, save_only=True)
|
|
|
823 else:
|
|
|
824 x_sum = np.sum(self.tbf, 1)
|
|
|
825 use_bounds = 'no'
|
|
|
826 if img_x_bounds[0] is not None and img_x_bounds[1] is not None:
|
|
|
827 if img_x_bounds[0] < 0:
|
|
|
828 msnc.illegal_value('img_x_bounds[0]', img_x_bounds[0], 'config file')
|
|
|
829 img_x_bounds[0] = 0
|
|
|
830 if not img_x_bounds[0] < img_x_bounds[1] <= x_sum.size:
|
|
|
831 msnc.illegal_value('img_x_bounds[1]', img_x_bounds[1], 'config file')
|
|
|
832 img_x_bounds[1] = x_sum.size
|
|
|
833 msnc.quickPlot((range(x_sum.size), x_sum),
|
|
|
834 ([img_x_bounds[0], img_x_bounds[0]], [x_sum.min(), x_sum.max()], 'r-'),
|
|
|
835 ([img_x_bounds[1], img_x_bounds[1]], [x_sum.min(), x_sum.max()], 'r-'),
|
|
|
836 title='sum over theta and y')
|
|
|
837 print(f'lower bound = {img_x_bounds[0]} (inclusive)\n'+
|
|
|
838 f'upper bound = {img_x_bounds[1]} (exclusive)]')
|
|
|
839 use_bounds = pyip.inputYesNo('Accept these bounds ([y]/n)?: ', blank=True)
|
|
|
840 if use_bounds == 'no':
|
|
|
841 fit = msnc.fitStep(y=x_sum, model='rectangle', form='atan')
|
|
|
842 x_low = fit.get('center1', None)
|
|
|
843 x_upp = fit.get('center2', None)
|
|
|
844 if x_low is not None and x_low >= 0 and x_upp is not None and x_low < x_upp < x_sum.size:
|
|
|
845 x_low = int(x_low-(x_upp-x_low)/10)
|
|
|
846 if x_low < 0:
|
|
|
847 x_low = 0
|
|
|
848 x_upp = int(x_upp+(x_upp-x_low)/10)
|
|
|
849 if x_upp >= x_sum.size:
|
|
|
850 x_upp = x_sum.size
|
|
|
851 msnc.quickPlot((range(x_sum.size), x_sum),
|
|
|
852 ([x_low, x_low], [x_sum.min(), x_sum.max()], 'r-'),
|
|
|
853 ([x_upp, x_upp], [x_sum.min(), x_sum.max()], 'r-'),
|
|
|
854 title='sum over theta and y')
|
|
|
855 print(f'lower bound = {x_low} (inclusive)\nupper bound = {x_upp} (exclusive)]')
|
|
|
856 use_fit = pyip.inputYesNo('Accept these bounds ([y]/n)?: ', blank=True)
|
|
|
857 if use_fit == 'no':
|
|
|
858 img_x_bounds = msnc.selectArrayBounds(x_sum, img_x_bounds[0], img_x_bounds[1],
|
|
|
859 num_x_min, 'sum over theta and y')
|
|
|
860 else:
|
|
|
861 img_x_bounds = [x_low, x_upp]
|
|
|
862 if num_x_min is not None and img_x_bounds[1]-img_x_bounds[0]+1 < num_x_min:
|
|
|
863 logging.warning('Image bounds and pixel size prevent seamless stacking')
|
|
|
864 msnc.quickPlot(range(img_x_bounds[0], img_x_bounds[1]),
|
|
|
865 x_sum[img_x_bounds[0]:img_x_bounds[1]],
|
|
|
866 title='sum over theta and y', path=self.output_folder,
|
|
|
867 save_fig=self.save_plots, save_only=True)
|
|
|
868 logging.debug(f'img_x_bounds: {img_x_bounds}')
|
|
|
869
|
|
|
870 if self.save_plots_only:
|
|
|
871 msnc.clearFig('bright field')
|
|
|
872 msnc.clearFig('sum over theta and y')
|
|
|
873 if title:
|
|
|
874 msnc.clearFig(title)
|
|
|
875
|
|
|
876 # Update config and save to file
|
|
|
877 if preprocess is None:
|
|
|
878 self.cf.config['preprocess'] = {'img_x_bounds' : img_x_bounds}
|
|
|
879 else:
|
|
|
880 preprocess['img_x_bounds'] = img_x_bounds
|
|
|
881 self.cf.saveFile(self.config_out)
|
|
|
882
|
|
|
883 def _setZoomOrSkip(self):
|
|
|
884 """Set zoom and/or theta skip to reduce memory the requirement for the analysis.
|
|
|
885 """
|
|
|
886 preprocess = self.config.get('preprocess')
|
|
|
887 zoom_perc = 100
|
|
|
888 if not self.galaxy_flag:
|
|
|
889 if preprocess is None or 'zoom_perc' not in preprocess:
|
|
|
890 if pyip.inputYesNo(
|
|
|
891 '\nDo you want to zoom in to reduce memory requirement (y/[n])? ',
|
|
|
892 blank=True) == 'yes':
|
|
|
893 zoom_perc = pyip.inputInt(' Enter zoom percentage [1, 100]: ',
|
|
|
894 min=1, max=100)
|
|
|
895 else:
|
|
|
896 zoom_perc = preprocess['zoom_perc']
|
|
|
897 if msnc.is_num(zoom_perc, 1., 100.):
|
|
|
898 zoom_perc = int(zoom_perc)
|
|
|
899 else:
|
|
|
900 msnc.illegal_value('zoom_perc', zoom_perc, 'config file')
|
|
|
901 zoom_perc = 100
|
|
|
902 num_theta_skip = 0
|
|
|
903 if not self.galaxy_flag:
|
|
|
904 if preprocess is None or 'num_theta_skip' not in preprocess:
|
|
|
905 if pyip.inputYesNo(
|
|
|
906 'Do you want to skip thetas to reduce memory requirement (y/[n])? ',
|
|
|
907 blank=True) == 'yes':
|
|
|
908 num_theta_skip = pyip.inputInt(' Enter the number skip theta interval'+
|
|
|
909 f' [0, {self.num_thetas-1}]: ', min=0, max=self.num_thetas-1)
|
|
|
910 else:
|
|
|
911 num_theta_skip = preprocess['num_theta_skip']
|
|
|
912 if not msnc.is_int(num_theta_skip, 0):
|
|
|
913 msnc.illegal_value('num_theta_skip', num_theta_skip, 'config file')
|
|
|
914 num_theta_skip = 0
|
|
|
915 logging.info(f'zoom_perc = {zoom_perc}')
|
|
|
916 logging.info(f'num_theta_skip = {num_theta_skip}')
|
|
|
917
|
|
|
918 # Update config and save to file
|
|
|
919 if preprocess is None:
|
|
|
920 self.cf.config['preprocess'] = {'zoom_perc' : zoom_perc,
|
|
|
921 'num_theta_skip' : num_theta_skip}
|
|
|
922 else:
|
|
|
923 preprocess['zoom_perc'] = zoom_perc
|
|
|
924 preprocess['num_theta_skip'] = num_theta_skip
|
|
|
925 self.cf.saveFile(self.config_out)
|
|
|
926
|
|
|
927 def _loadTomo(self, base_name, index, required=False):
|
|
|
928 """Load a tomography stack.
|
|
|
929 """
|
|
|
930 # stack order: row,theta,column
|
|
|
931 zoom_perc = None
|
|
|
932 preprocess = self.config.get('preprocess')
|
|
|
933 if preprocess:
|
|
|
934 zoom_perc = preprocess.get('zoom_perc')
|
|
|
935 if zoom_perc is None or zoom_perc == 100:
|
|
|
936 title = f'{base_name} fullres'
|
|
|
937 else:
|
|
|
938 title = f'{base_name} {zoom_perc}p'
|
|
|
939 title += f'_{index}'
|
|
|
940 tomo_file = re.sub(r"\s+", '_', f'{self.output_folder}/{title}.npy')
|
|
|
941 load_flag = 'no'
|
|
|
942 available = False
|
|
|
943 if os.path.isfile(tomo_file):
|
|
|
944 available = True
|
|
|
945 if required:
|
|
|
946 load_flag = 'yes'
|
|
|
947 else:
|
|
|
948 load_flag = pyip.inputYesNo(f'\nDo you want to load {tomo_file} (y/n)? ')
|
|
|
949 stack = np.array([])
|
|
|
950 if load_flag == 'yes':
|
|
|
951 t0 = time()
|
|
|
952 logging.info(f'Loading {tomo_file} ...')
|
|
|
953 try:
|
|
|
954 stack = np.load(tomo_file)
|
|
|
955 except IOError or ValueError:
|
|
|
956 stack = np.array([])
|
|
|
957 logging.error(f'Error loading {tomo_file}')
|
|
|
958 logging.info(f'... done in {time()-t0:.2f} seconds!')
|
|
|
959 if stack.size:
|
|
|
960 msnc.quickImshow(stack[:,0,:], title=title, path=self.output_folder,
|
|
|
961 save_fig=self.save_plots, save_only=self.save_plots_only)
|
|
|
962 return stack, available
|
|
|
963
|
|
|
964 def _saveTomo(self, base_name, stack, index=None):
|
|
|
965 """Save a tomography stack.
|
|
|
966 """
|
|
|
967 zoom_perc = None
|
|
|
968 preprocess = self.config.get('preprocess')
|
|
|
969 if preprocess:
|
|
|
970 zoom_perc = preprocess.get('zoom_perc')
|
|
|
971 if zoom_perc is None or zoom_perc == 100:
|
|
|
972 title = f'{base_name} fullres'
|
|
|
973 else:
|
|
|
974 title = f'{base_name} {zoom_perc}p'
|
|
|
975 if index:
|
|
|
976 title += f'_{index}'
|
|
|
977 tomo_file = re.sub(r"\s+", '_', f'{self.output_folder}/{title}.npy')
|
|
|
978 t0 = time()
|
|
|
979 logging.info(f'Saving {tomo_file} ...')
|
|
|
980 np.save(tomo_file, stack)
|
|
|
981 logging.info(f'... done in {time()-t0:.2f} seconds!')
|
|
|
982
|
|
|
983 def _genTomo(self, tomo_stack_files, available_stacks):
|
|
|
984 """Generate tomography fields.
|
|
|
985 """
|
|
|
986 stacks = self.config['stack_info']['stacks']
|
|
|
987 assert(len(self.tomo_stacks) == self.config['stack_info']['num'])
|
|
|
988 assert(len(self.tomo_stacks) == len(stacks))
|
|
|
989 if len(available_stacks) != len(stacks):
|
|
|
990 logging.warning('Illegal dimension of available_stacks in _genTomo'+
|
|
|
991 f'({len(available_stacks)}');
|
|
|
992 available_stacks = [False]*self.num_tomo_stacks
|
|
|
993
|
|
|
994 preprocess = self.config.get('preprocess')
|
|
|
995 if preprocess is None:
|
|
|
996 img_x_bounds = [0, self.tbf.shape[0]]
|
|
|
997 img_y_bounds = [0, self.tbf.shape[1]]
|
|
|
998 zoom_perc = preprocess['zoom_perc']
|
|
|
999 num_theta_skip = preprocess['num_theta_skip']
|
|
|
1000 else:
|
|
|
1001 img_x_bounds = preprocess.get('img_x_bounds', [0, self.tbf.shape[0]])
|
|
|
1002 img_y_bounds = preprocess.get('img_y_bounds', [0, self.tbf.shape[1]])
|
|
|
1003 zoom_perc = 100
|
|
|
1004 num_theta_skip = 0
|
|
|
1005
|
|
|
1006 if self.tdf.size:
|
|
|
1007 tdf = self.tdf[img_x_bounds[0]:img_x_bounds[1],img_y_bounds[0]:img_y_bounds[1]]
|
|
|
1008 else:
|
|
|
1009 logging.warning('Dark field unavailable')
|
|
|
1010 if not self.tbf.size:
|
|
|
1011 raise ValueError('Bright field unavailable')
|
|
|
1012 tbf = self.tbf[img_x_bounds[0]:img_x_bounds[1],img_y_bounds[0]:img_y_bounds[1]]
|
|
|
1013
|
|
|
1014 for i,stack in enumerate(stacks):
|
|
|
1015 # Check if stack is already loaded or available
|
|
|
1016 if self.tomo_stacks[i].size or available_stacks[i]:
|
|
|
1017 continue
|
|
|
1018
|
|
|
1019 # Load a stack of tomography images
|
|
|
1020 t0 = time()
|
|
|
1021 tomo_stack = msnc.loadImageStack(tomo_stack_files[i], self.config['data_filetype'],
|
|
|
1022 stack['img_offset'], self.config['theta_range']['num'], num_theta_skip,
|
|
|
1023 img_x_bounds, img_y_bounds)
|
|
|
1024 tomo_stack = tomo_stack.astype('float64')
|
|
|
1025 logging.debug(f'loading took {time()-t0:.2f} seconds!')
|
|
|
1026
|
|
|
1027 # Subtract dark field
|
|
|
1028 if self.tdf.size:
|
|
|
1029 t0 = time()
|
|
|
1030 with set_numexpr_threads(self.ncore):
|
|
|
1031 ne.evaluate('tomo_stack-tdf', out=tomo_stack)
|
|
|
1032 logging.debug(f'subtracting dark field took {time()-t0:.2f} seconds!')
|
|
|
1033
|
|
|
1034 # Normalize
|
|
|
1035 t0 = time()
|
|
|
1036 with set_numexpr_threads(self.ncore):
|
|
|
1037 ne.evaluate('tomo_stack/tbf', out=tomo_stack, truediv=True)
|
|
|
1038 logging.debug(f'normalizing took {time()-t0:.2f} seconds!')
|
|
|
1039
|
|
|
1040 # Remove non-positive values and linearize data
|
|
|
1041 t0 = time()
|
|
|
1042 cutoff = 1.e-6
|
|
|
1043 with set_numexpr_threads(self.ncore):
|
|
|
1044 ne.evaluate('where(tomo_stack<cutoff, cutoff, tomo_stack)', out=tomo_stack)
|
|
|
1045 with set_numexpr_threads(self.ncore):
|
|
|
1046 ne.evaluate('-log(tomo_stack)', out=tomo_stack)
|
|
|
1047 logging.debug('removing non-positive values and linearizing data took '+
|
|
|
1048 f'{time()-t0:.2f} seconds!')
|
|
|
1049
|
|
|
1050 # Get rid of nans/infs that may be introduced by normalization
|
|
|
1051 t0 = time()
|
|
|
1052 np.where(np.isfinite(tomo_stack), tomo_stack, 0.)
|
|
|
1053 logging.debug(f'remove nans/infs took {time()-t0:.2f} seconds!')
|
|
|
1054
|
|
|
1055 # Downsize tomography stack to smaller size
|
|
|
1056 tomo_stack = tomo_stack.astype('float32')
|
|
|
1057 if not self.galaxy_flag:
|
|
|
1058 index = stack['index']
|
|
|
1059 title = f'red stack fullres {index}'
|
|
|
1060 if not self.test_mode:
|
|
|
1061 msnc.quickImshow(tomo_stack[0,:,:], title=title, path=self.output_folder,
|
|
|
1062 save_fig=self.save_plots, save_only=self.save_plots_only)
|
|
|
1063 if zoom_perc != 100:
|
|
|
1064 t0 = time()
|
|
|
1065 logging.info(f'Zooming in ...')
|
|
|
1066 tomo_zoom_list = []
|
|
|
1067 for j in range(tomo_stack.shape[0]):
|
|
|
1068 tomo_zoom = spi.zoom(tomo_stack[j,:,:], 0.01*zoom_perc)
|
|
|
1069 tomo_zoom_list.append(tomo_zoom)
|
|
|
1070 tomo_stack = np.stack([tomo_zoom for tomo_zoom in tomo_zoom_list])
|
|
|
1071 logging.info(f'... done in {time()-t0:.2f} seconds!')
|
|
|
1072 del tomo_zoom_list
|
|
|
1073 if not self.galaxy_flag:
|
|
|
1074 title = f'red stack {zoom_perc}p {index}'
|
|
|
1075 if not self.test_mode:
|
|
|
1076 msnc.quickImshow(tomo_stack[0,:,:], title=title, path=self.output_folder,
|
|
|
1077 save_fig=self.save_plots, save_only=self.save_plots_only)
|
|
|
1078
|
|
|
1079 # Convert tomography stack from theta,row,column to row,theta,column
|
|
|
1080 tomo_stack = np.swapaxes(tomo_stack, 0, 1)
|
|
|
1081
|
|
|
1082 # Save tomography stack to file
|
|
|
1083 if not self.galaxy_flag:
|
|
|
1084 if not self.test_mode:
|
|
|
1085 self._saveTomo('red stack', tomo_stack, index)
|
|
|
1086 else:
|
|
|
1087 np.savetxt(self.output_folder+f'red_stack_{index}.txt',
|
|
|
1088 tomo_stack[0,:,:], fmt='%.6e')
|
|
|
1089
|
|
|
1090 # Combine stacks
|
|
|
1091 t0 = time()
|
|
|
1092 self.tomo_stacks[i] = tomo_stack
|
|
|
1093 logging.debug(f'combining nstack took {time()-t0:.2f} seconds!')
|
|
|
1094
|
|
|
1095 # Update config and save to file
|
|
|
1096 stack['preprocessed'] = True
|
|
|
1097 self.cf.saveFile(self.config_out)
|
|
|
1098
|
|
|
1099 if self.tdf.size:
|
|
|
1100 del tdf
|
|
|
1101 del tbf
|
|
|
1102
|
|
|
1103 def _reconstructOnePlane(self, tomo_plane_T, center, thetas_deg, eff_pixel_size,
|
|
|
1104 cross_sectional_dim, plot_sinogram=True):
|
|
|
1105 """Invert the sinogram for a single tomography plane.
|
|
|
1106 """
|
|
|
1107 # tomo_plane_T index order: column,theta
|
|
|
1108 assert(0 <= center < tomo_plane_T.shape[0])
|
|
|
1109 center_offset = center-tomo_plane_T.shape[0]/2
|
|
|
1110 two_offset = 2*int(np.round(center_offset))
|
|
|
1111 two_offset_abs = np.abs(two_offset)
|
|
|
1112 max_rad = int(0.5*(cross_sectional_dim/eff_pixel_size)*1.1) # 10% slack to avoid edge effects
|
|
|
1113 dist_from_edge = max(1, int(np.floor((tomo_plane_T.shape[0]-two_offset_abs)/2.)-max_rad))
|
|
|
1114 if two_offset >= 0:
|
|
|
1115 logging.debug(f'sinogram range = [{two_offset+dist_from_edge}, {-dist_from_edge}]')
|
|
|
1116 sinogram = tomo_plane_T[two_offset+dist_from_edge:-dist_from_edge,:]
|
|
|
1117 else:
|
|
|
1118 logging.debug(f'sinogram range = [{dist_from_edge}, {two_offset-dist_from_edge}]')
|
|
|
1119 sinogram = tomo_plane_T[dist_from_edge:two_offset-dist_from_edge,:]
|
|
|
1120 if plot_sinogram:
|
|
|
1121 msnc.quickImshow(sinogram.T, f'sinogram center offset{center_offset:.2f}',
|
|
|
1122 path=self.output_folder, save_fig=self.save_plots,
|
|
|
1123 save_only=self.save_plots_only, aspect='auto')
|
|
|
1124
|
|
|
1125 # Inverting sinogram
|
|
|
1126 t0 = time()
|
|
|
1127 recon_sinogram = iradon(sinogram, theta=thetas_deg, circle=True)
|
|
|
1128 logging.debug(f'inverting sinogram took {time()-t0:.2f} seconds!')
|
|
|
1129 del sinogram
|
|
|
1130
|
|
|
1131 # Removing ring artifacts
|
|
|
1132 # RV parameters for the denoise, gaussian, and ring removal will be different for different feature sizes
|
|
|
1133 t0 = time()
|
|
|
1134 # recon_sinogram = filters.gaussian(recon_sinogram, 3.0)
|
|
|
1135 recon_sinogram = spi.gaussian_filter(recon_sinogram, 0.5)
|
|
|
1136 recon_clean = np.expand_dims(recon_sinogram, axis=0)
|
|
|
1137 del recon_sinogram
|
|
|
1138 recon_clean = tomopy.misc.corr.remove_ring(recon_clean, rwidth=17)
|
|
|
1139 logging.debug(f'filtering and removing ring artifact took {time()-t0:.2f} seconds!')
|
|
|
1140 return recon_clean
|
|
|
1141
|
|
|
1142 def _plotEdgesOnePlane(self, recon_plane, base_name, weight=0.001):
|
|
|
1143 # RV parameters for the denoise, gaussian, and ring removal will be different for different feature sizes
|
|
|
1144 edges = denoise_tv_chambolle(recon_plane, weight = weight)
|
|
|
1145 vmax = np.max(edges[0,:,:])
|
|
|
1146 vmin = -vmax
|
|
|
1147 msnc.quickImshow(edges[0,:,:], f'{base_name} coolwarm', path=self.output_folder,
|
|
|
1148 save_fig=self.save_plots, save_only=self.save_plots_only, cmap='coolwarm')
|
|
|
1149 msnc.quickImshow(edges[0,:,:], f'{base_name} gray', path=self.output_folder,
|
|
|
1150 save_fig=self.save_plots, save_only=self.save_plots_only, cmap='gray',
|
|
|
1151 vmin=vmin, vmax=vmax)
|
|
|
1152 del edges
|
|
|
1153
|
|
|
1154 def _findCenterOnePlane(self, sinogram, row, thetas_deg, eff_pixel_size, cross_sectional_dim,
|
|
|
1155 tol=0.1):
|
|
|
1156 """Find center for a single tomography plane.
|
|
|
1157 """
|
|
|
1158 # sinogram index order: theta,column
|
|
|
1159 # need index order column,theta for iradon, so take transpose
|
|
|
1160 sinogram_T = sinogram.T
|
|
|
1161 center = sinogram.shape[1]/2
|
|
|
1162
|
|
|
1163 # try automatic center finding routines for initial value
|
|
|
1164 tomo_center = tomopy.find_center_vo(sinogram)
|
|
|
1165 center_offset_vo = tomo_center-center
|
|
|
1166 print(f'Center at row {row} using Nghia Vo’s method = {center_offset_vo:.2f}')
|
|
|
1167 recon_plane = self._reconstructOnePlane(sinogram_T, tomo_center, thetas_deg,
|
|
|
1168 eff_pixel_size, cross_sectional_dim, False)
|
|
|
1169 base_name=f'edges row{row} center_offset_vo{center_offset_vo:.2f}'
|
|
|
1170 self._plotEdgesOnePlane(recon_plane, base_name)
|
|
|
1171 if pyip.inputYesNo('Try finding center using phase correlation (y/[n])? ',
|
|
|
1172 blank=True) == 'yes':
|
|
|
1173 tomo_center = tomopy.find_center_pc(sinogram, sinogram, tol=0.1,
|
|
|
1174 rotc_guess=tomo_center)
|
|
|
1175 error = 1.
|
|
|
1176 while error > tol:
|
|
|
1177 prev = tomo_center
|
|
|
1178 tomo_center = tomopy.find_center_pc(sinogram, sinogram, tol=tol,
|
|
|
1179 rotc_guess=tomo_center)
|
|
|
1180 error = np.abs(tomo_center-prev)
|
|
|
1181 center_offset = tomo_center-center
|
|
|
1182 print(f'Center at row {row} using phase correlation = {center_offset:.2f}')
|
|
|
1183 recon_plane = self._reconstructOnePlane(sinogram_T, tomo_center, thetas_deg,
|
|
|
1184 eff_pixel_size, cross_sectional_dim, False)
|
|
|
1185 base_name=f'edges row{row} center_offset{center_offset:.2f}'
|
|
|
1186 self._plotEdgesOnePlane(recon_plane, base_name)
|
|
|
1187 if pyip.inputYesNo('Accept a center location ([y]) or continue search (n)? ',
|
|
|
1188 blank=True) != 'no':
|
|
|
1189 del sinogram_T
|
|
|
1190 del recon_plane
|
|
|
1191 center_offset = pyip.inputNum(
|
|
|
1192 f' Enter chosen center offset [{-int(center)}, {int(center)}] '+
|
|
|
1193 f'([{center_offset_vo}])): ', blank=True)
|
|
|
1194 if center_offset == '':
|
|
|
1195 center_offset = center_offset_vo
|
|
|
1196 return float(center_offset)
|
|
|
1197
|
|
|
1198 while True:
|
|
|
1199 center_offset_low = pyip.inputInt('\nEnter lower bound for center offset '+
|
|
|
1200 f'[{-int(center)}, {int(center)}]: ', min=-int(center), max=int(center))
|
|
|
1201 center_offset_upp = pyip.inputInt('Enter upper bound for center offset '+
|
|
|
1202 f'[{center_offset_low}, {int(center)}]: ',
|
|
|
1203 min=center_offset_low, max=int(center))
|
|
|
1204 if center_offset_upp == center_offset_low:
|
|
|
1205 center_offset_step = 1
|
|
|
1206 else:
|
|
|
1207 center_offset_step = pyip.inputInt('Enter step size for center offset search '+
|
|
|
1208 f'[1, {center_offset_upp-center_offset_low}]: ',
|
|
|
1209 min=1, max=center_offset_upp-center_offset_low)
|
|
|
1210 for center_offset in range(center_offset_low, center_offset_upp+center_offset_step,
|
|
|
1211 center_offset_step):
|
|
|
1212 logging.info(f'center_offset = {center_offset}')
|
|
|
1213 recon_plane = self._reconstructOnePlane(sinogram_T, center_offset+center,
|
|
|
1214 thetas_deg, eff_pixel_size, cross_sectional_dim, False)
|
|
|
1215 base_name=f'edges row{row} center_offset{center_offset}'
|
|
|
1216 self._plotEdgesOnePlane(recon_plane, base_name)
|
|
|
1217 if pyip.inputInt('\nContinue (0) or end the search (1): ', min=0, max=1):
|
|
|
1218 break
|
|
|
1219
|
|
|
1220 del sinogram_T
|
|
|
1221 del recon_plane
|
|
|
1222 center_offset = pyip.inputNum(f' Enter chosen center offset '+
|
|
|
1223 f'[{-int(center)}, {int(center)}]: ', min=-int(center), max=int(center))
|
|
|
1224 return float(center_offset)
|
|
|
1225
|
|
|
1226 def _reconstructOneTomoStack(self, tomo_stack, thetas, row_bounds=None,
|
|
|
1227 center_offsets=[], sigma=0.1, rwidth=30, ncore=1, algorithm='gridrec',
|
|
|
1228 run_secondary_sirt=False, secondary_iter=100):
|
|
|
1229 """reconstruct a single tomography stack.
|
|
|
1230 """
|
|
|
1231 # stack order: row,theta,column
|
|
|
1232 # thetas must be in radians
|
|
|
1233 # centers_offset: tomography axis shift in pixels relative to column center
|
|
|
1234 # RV should we remove stripes?
|
|
|
1235 # https://tomopy.readthedocs.io/en/latest/api/tomopy.prep.stripe.html
|
|
|
1236 # RV should we remove rings?
|
|
|
1237 # https://tomopy.readthedocs.io/en/latest/api/tomopy.misc.corr.html
|
|
|
1238 # RV: Add an option to do (extra) secondary iterations later or to do some sort of convergence test?
|
|
|
1239 if row_bounds is None:
|
|
|
1240 row_bounds = [0, tomo_stack.shape[0]]
|
|
|
1241 else:
|
|
|
1242 if not (0 <= row_bounds[0] <= row_bounds[1] <= tomo_stack.shape[0]):
|
|
|
1243 raise ValueError('Illegal row bounds in reconstructOneTomoStack')
|
|
|
1244 if thetas.size != tomo_stack.shape[1]:
|
|
|
1245 raise ValueError('theta dimension mismatch in reconstructOneTomoStack')
|
|
|
1246 if not len(center_offsets):
|
|
|
1247 centers = np.zeros((row_bounds[1]-row_bounds[0]))
|
|
|
1248 elif len(center_offsets) == 2:
|
|
|
1249 centers = np.linspace(center_offsets[0], center_offsets[1],
|
|
|
1250 row_bounds[1]-row_bounds[0])
|
|
|
1251 else:
|
|
|
1252 if center_offsets.size != row_bounds[1]-row_bounds[0]:
|
|
|
1253 raise ValueError('center_offsets dimension mismatch in reconstructOneTomoStack')
|
|
|
1254 centers = center_offsets
|
|
|
1255 centers += tomo_stack.shape[2]/2
|
|
|
1256 if True:
|
|
|
1257 tomo_stack = tomopy.prep.stripe.remove_stripe_fw(tomo_stack[row_bounds[0]:row_bounds[1]],
|
|
|
1258 sigma=sigma, ncore=ncore)
|
|
|
1259 else:
|
|
|
1260 tomo_stack = tomo_stack[row_bounds[0]:row_bounds[1]]
|
|
|
1261 tomo_recon_stack = tomopy.recon(tomo_stack, thetas, centers, sinogram_order=True,
|
|
|
1262 algorithm=algorithm, ncore=ncore)
|
|
|
1263 if run_secondary_sirt and secondary_iter > 0:
|
|
|
1264 #options = {'method':'SIRT_CUDA', 'proj_type':'cuda', 'num_iter':secondary_iter}
|
|
|
1265 #RV: doesn't work for me: "Error: CUDA error 803: system has unsupported display driver /
|
|
|
1266 # cuda driver combination."
|
|
|
1267 #options = {'method':'SIRT', 'proj_type':'linear', 'MinConstraint': 0, 'num_iter':secondary_iter}
|
|
|
1268 #SIRT did not finish while running overnight
|
|
|
1269 #options = {'method':'SART', 'proj_type':'linear', 'num_iter':secondary_iter}
|
|
|
1270 options = {'method':'SART', 'proj_type':'linear', 'MinConstraint': 0, 'num_iter':secondary_iter}
|
|
|
1271 tomo_recon_stack = tomopy.recon(tomo_stack, thetas, centers, init_recon=tomo_recon_stack,
|
|
|
1272 options=options, sinogram_order=True, algorithm=tomopy.astra, ncore=ncore)
|
|
|
1273 if True:
|
|
|
1274 tomopy.misc.corr.remove_ring(tomo_recon_stack, rwidth=rwidth, out=tomo_recon_stack)
|
|
|
1275 return tomo_recon_stack
|
|
|
1276
|
|
|
1277 def genTomoStacks(self, tdf_files=None, tbf_files=None, tomo_stack_files=None,
|
|
|
1278 dark_field_pngname=None, bright_field_pngname=None, tomo_field_pngname=None,
|
|
|
1279 detectorbounds_pngname=None, output_name=None):
|
|
|
1280 """Preprocess tomography images.
|
|
|
1281 """
|
|
|
1282 # Try loading any already preprocessed stacks (skip in Galaxy)
|
|
|
1283 # preprocessed stack order for each one in stack: row,theta,column
|
|
|
1284 stack_info = self.config['stack_info']
|
|
|
1285 stacks = stack_info['stacks']
|
|
|
1286 num_tomo_stacks = stack_info['num']
|
|
|
1287 assert(num_tomo_stacks == len(self.tomo_stacks))
|
|
|
1288 available_stacks = [False]*num_tomo_stacks
|
|
|
1289 if self.galaxy_flag:
|
|
|
1290 assert(tdf_files is None or isinstance(tdf_files, list))
|
|
|
1291 assert(isinstance(tbf_files, list))
|
|
|
1292 assert(isinstance(tomo_stack_files, list))
|
|
|
1293 assert(num_tomo_stacks == len(tomo_stack_files))
|
|
|
1294 assert(isinstance(dark_field_pngname, str))
|
|
|
1295 assert(isinstance(bright_field_pngname, str))
|
|
|
1296 assert(isinstance(tomo_field_pngname, str))
|
|
|
1297 assert(isinstance(detectorbounds_pngname, str))
|
|
|
1298 assert(isinstance(output_name, str))
|
|
|
1299 else:
|
|
|
1300 if tdf_files:
|
|
|
1301 logging.warning('Ignoring tdf_files in genTomoStacks (only for Galaxy)')
|
|
|
1302 if tbf_files:
|
|
|
1303 logging.warning('Ignoring tbf_files in genTomoStacks (only for Galaxy)')
|
|
|
1304 if tomo_stack_files:
|
|
|
1305 logging.warning('Ignoring tomo_stack_files in genTomoStacks (only for Galaxy)')
|
|
|
1306 tdf_files, tbf_files, tomo_stack_files = self.findImageFiles()
|
|
|
1307 if not self.is_valid:
|
|
|
1308 return
|
|
|
1309 for i,stack in enumerate(stacks):
|
|
|
1310 if not self.tomo_stacks[i].size and stack.get('preprocessed', False):
|
|
|
1311 self.tomo_stacks[i], available_stacks[i] = \
|
|
|
1312 self._loadTomo('red stack', stack['index'])
|
|
|
1313 if dark_field_pngname:
|
|
|
1314 logging.warning('Ignoring dark_field_pngname in genTomoStacks (only for Galaxy)')
|
|
|
1315 if bright_field_pngname:
|
|
|
1316 logging.warning('Ignoring bright_field_pngname in genTomoStacks (only for Galaxy)')
|
|
|
1317 if tomo_field_pngname:
|
|
|
1318 logging.warning('Ignoring tomo_field_pngname in genTomoStacks (only for Galaxy)')
|
|
|
1319 if detectorbounds_pngname:
|
|
|
1320 logging.warning('Ignoring detectorbounds_pngname in genTomoStacks '+
|
|
|
1321 '(only used in Galaxy)')
|
|
|
1322 if output_name:
|
|
|
1323 logging.warning('Ignoring output_name in genTomoStacks '+
|
|
|
1324 '(only used in Galaxy)')
|
|
|
1325
|
|
|
1326 # Preprocess any unloaded stacks
|
|
|
1327 if False in available_stacks:
|
|
|
1328 logging.debug('Preprocessing tomography images')
|
|
|
1329
|
|
|
1330 # Check required image files (skip in Galaxy)
|
|
|
1331 if not self.galaxy_flag:
|
|
|
1332 self.selectImageRanges(available_stacks)
|
|
|
1333 if not self.is_valid:
|
|
|
1334 return
|
|
|
1335
|
|
|
1336 # Generate dark field
|
|
|
1337 if tdf_files:
|
|
|
1338 self._genDark(tdf_files, dark_field_pngname)
|
|
|
1339
|
|
|
1340 # Generate bright field
|
|
|
1341 self._genBright(tbf_files, bright_field_pngname)
|
|
|
1342
|
|
|
1343 # Set vertical detector bounds for image stack
|
|
|
1344 self._setDetectorBounds(tomo_stack_files, tomo_field_pngname, detectorbounds_pngname)
|
|
|
1345
|
|
|
1346 # Set zoom and/or theta skip to reduce memory the requirement
|
|
|
1347 self._setZoomOrSkip()
|
|
|
1348
|
|
|
1349 # Generate tomography fields
|
|
|
1350 self._genTomo(tomo_stack_files, available_stacks)
|
|
|
1351
|
|
|
1352 # Save tomography stack to file
|
|
|
1353 if self.galaxy_flag:
|
|
|
1354 t0 = time()
|
|
|
1355 logging.info(f'Saving preprocessed tomography stack to file ...')
|
|
|
1356 save_stacks = {f'set_{stack["index"]}':tomo_stack
|
|
|
1357 for stack,tomo_stack in zip(stacks,self.tomo_stacks)}
|
|
|
1358 np.savez(output_name, **save_stacks)
|
|
|
1359 logging.info(f'... done in {time()-t0:.2f} seconds!')
|
|
|
1360
|
|
|
1361 del available_stacks
|
|
|
1362
|
|
|
1363 # Adjust sample reference height, update config and save to file
|
|
|
1364 preprocess = self.config.get('preprocess')
|
|
|
1365 if preprocess is None:
|
|
|
1366 img_x_bounds = [0, self.tbf.shape[0]]
|
|
|
1367 else:
|
|
|
1368 img_x_bounds = preprocess.get('img_x_bounds', [0, self.tbf.shape[0]])
|
|
|
1369 pixel_size = self.config['detector']['pixel_size']
|
|
|
1370 if pixel_size is None:
|
|
|
1371 raise ValueError('Detector pixel size unavailable')
|
|
|
1372 pixel_size *= img_x_bounds[0]
|
|
|
1373 for stack in stacks:
|
|
|
1374 stack['ref_height'] = stack['ref_height']+pixel_size
|
|
|
1375 self.cf.saveFile(self.config_out)
|
|
|
1376
|
|
|
1377 def findCenters(self):
|
|
|
1378 """Find rotation axis centers for the tomography stacks.
|
|
|
1379 """
|
|
|
1380 logging.debug('Find centers for tomography stacks')
|
|
|
1381 stacks = self.config['stack_info']['stacks']
|
|
|
1382 available_stacks = [stack['index'] for stack in stacks if stack.get('preprocessed', False)]
|
|
|
1383 logging.debug('Avaliable stacks: {available_stacks}')
|
|
|
1384
|
|
|
1385 # Check for valid available center stack index
|
|
|
1386 find_center = self.config.get('find_center')
|
|
|
1387 if find_center and 'center_stack_index' in find_center:
|
|
|
1388 center_stack_index = find_center['center_stack_index']
|
|
|
1389 if (not isinstance(center_stack_index, int) or
|
|
|
1390 center_stack_index not in available_stacks):
|
|
|
1391 msnc.illegal_value('center_stack_index', center_stack_index, 'config file')
|
|
|
1392 else:
|
|
|
1393 if self.test_mode:
|
|
|
1394 find_center['completed'] = True
|
|
|
1395 self.cf.saveFile(self.config_out)
|
|
|
1396 return
|
|
|
1397 print('\nFound calibration center offset info for stack '+
|
|
|
1398 f'{center_stack_index}')
|
|
|
1399 if pyip.inputYesNo('Do you want to use this again (y/n)? ') == 'yes':
|
|
|
1400 find_center['completed'] = True
|
|
|
1401 self.cf.saveFile(self.config_out)
|
|
|
1402 return
|
|
|
1403
|
|
|
1404 # Load the required preprocessed stack
|
|
|
1405 # preprocessed stack order: row,theta,column
|
|
|
1406 num_tomo_stacks = self.config['stack_info']['num']
|
|
|
1407 assert(len(stacks) == num_tomo_stacks)
|
|
|
1408 assert(len(self.tomo_stacks) == num_tomo_stacks)
|
|
|
1409 if num_tomo_stacks == 1:
|
|
|
1410 center_stack_index = stacks[0]['index']
|
|
|
1411 if not self.tomo_stacks[0].size:
|
|
|
1412 self.tomo_stacks[0], available = self._loadTomo('red stack', center_stack_index,
|
|
|
1413 required=True)
|
|
|
1414 center_stack = self.tomo_stacks[0]
|
|
|
1415 if not center_stack.size:
|
|
|
1416 logging.error('Unable to load the required preprocessed tomography stack')
|
|
|
1417 return
|
|
|
1418 assert(stacks[0].get('preprocessed', False) == True)
|
|
|
1419 else:
|
|
|
1420 while True:
|
|
|
1421 center_stack_index = pyip.inputInt('\nEnter tomography stack index to get '
|
|
|
1422 f'rotation axis centers {available_stacks}: ')
|
|
|
1423 while center_stack_index not in available_stacks:
|
|
|
1424 center_stack_index = pyip.inputInt('\nEnter tomography stack index to get '
|
|
|
1425 f'rotation axis centers {available_stacks}: ')
|
|
|
1426 tomo_stack_index = available_stacks.index(center_stack_index)
|
|
|
1427 if not self.tomo_stacks[tomo_stack_index].size:
|
|
|
1428 self.tomo_stacks[tomo_stack_index], available = self._loadTomo(
|
|
|
1429 'red stack', center_stack_index, required=True)
|
|
|
1430 center_stack = self.tomo_stacks[tomo_stack_index]
|
|
|
1431 if not center_stack.size:
|
|
|
1432 logging.error(f'Unable to load the {center_stack_index}th '+
|
|
|
1433 'preprocessed tomography stack, pick another one')
|
|
|
1434 else:
|
|
|
1435 break
|
|
|
1436 assert(stacks[tomo_stack_index].get('preprocessed', False) == True)
|
|
|
1437 if find_center is None:
|
|
|
1438 self.config['find_center'] = {'center_stack_index' : center_stack_index}
|
|
|
1439 find_center = self.config['find_center']
|
|
|
1440
|
|
|
1441 # Set thetas (in degrees)
|
|
|
1442 theta_range = self.config['theta_range']
|
|
|
1443 theta_start = theta_range['start']
|
|
|
1444 theta_end = theta_range['end']
|
|
|
1445 num_theta = theta_range['num']
|
|
|
1446 num_theta_skip = self.config['preprocess'].get('num_theta_skip', 0)
|
|
|
1447 thetas_deg = np.linspace(theta_start, theta_end, int(num_theta/(num_theta_skip+1)),
|
|
|
1448 endpoint=False)
|
|
|
1449
|
|
|
1450 # Get non-overlapping sample row boundaries
|
|
|
1451 zoom_perc = self.config['preprocess'].get('zoom_perc', 100)
|
|
|
1452 pixel_size = self.config['detector']['pixel_size']
|
|
|
1453 if pixel_size is None:
|
|
|
1454 raise ValueError('Detector pixel size unavailable')
|
|
|
1455 eff_pixel_size = 100.*pixel_size/zoom_perc
|
|
|
1456 logging.debug(f'eff_pixel_size = {eff_pixel_size}')
|
|
|
1457 tomo_ref_heights = [stack['ref_height'] for stack in stacks]
|
|
|
1458 if num_tomo_stacks == 1:
|
|
|
1459 n1 = 0
|
|
|
1460 height = center_stack.shape[0]*eff_pixel_size
|
|
|
1461 if pyip.inputYesNo('\nDo you want to reconstruct the full samply height '+
|
|
|
1462 f'({height:.3f} mm) (y/n)? ') == 'no':
|
|
|
1463 height = pyip.inputNum('\nEnter the desired reconstructed sample height '+
|
|
|
1464 f'in mm [0, {height:.3f}]: ', min=0, max=height)
|
|
|
1465 n1 = int(0.5*(center_stack.shape[0]-height/eff_pixel_size))
|
|
|
1466 else:
|
|
|
1467 n1 = int((1.+(tomo_ref_heights[0]+center_stack.shape[0]*eff_pixel_size-
|
|
|
1468 tomo_ref_heights[1])/eff_pixel_size)/2)
|
|
|
1469 n2 = center_stack.shape[0]-n1
|
|
|
1470 logging.info(f'n1 = {n1}, n2 = {n2} (n2-n1) = {(n2-n1)*eff_pixel_size:.3f} mm')
|
|
|
1471 if not center_stack.size:
|
|
|
1472 RuntimeError('Center stack not loaded')
|
|
|
1473 msnc.quickImshow(center_stack[:,0,:], title=f'center stack theta={theta_start}',
|
|
|
1474 path=self.output_folder, save_fig=self.save_plots, save_only=self.save_plots_only)
|
|
|
1475
|
|
|
1476 # Get cross sectional diameter in mm
|
|
|
1477 cross_sectional_dim = center_stack.shape[2]*eff_pixel_size
|
|
|
1478 logging.debug(f'cross_sectional_dim = {cross_sectional_dim}')
|
|
|
1479
|
|
|
1480 # Determine center offset at sample row boundaries
|
|
|
1481 logging.info('Determine center offset at sample row boundaries')
|
|
|
1482
|
|
|
1483 # Lower row center
|
|
|
1484 use_row = False
|
|
|
1485 use_center = False
|
|
|
1486 row = find_center.get('lower_row')
|
|
|
1487 if msnc.is_int(row, n1, n2-2):
|
|
|
1488 msnc.quickImshow(center_stack[:,0,:], title=f'theta={theta_start}', aspect='auto')
|
|
|
1489 use_row = pyip.inputYesNo('\nCurrent row index for lower center = '
|
|
|
1490 f'{row}, use this value (y/n)? ')
|
|
|
1491 if self.save_plots_only:
|
|
|
1492 msnc.clearFig(f'theta={theta_start}')
|
|
|
1493 if use_row:
|
|
|
1494 center_offset = find_center.get('lower_center_offset')
|
|
|
1495 if msnc.is_num(center_offset):
|
|
|
1496 use_center = pyip.inputYesNo('Current lower center offset = '+
|
|
|
1497 f'{center_offset}, use this value (y/n)? ')
|
|
|
1498 if not use_center:
|
|
|
1499 if not use_row:
|
|
|
1500 msnc.quickImshow(center_stack[:,0,:], title=f'theta={theta_start}', aspect='auto')
|
|
|
1501 row = pyip.inputInt('\nEnter row index to find lower center '+
|
|
|
1502 f'[[{n1}], {n2-2}]: ', min=n1, max=n2-2, blank=True)
|
|
|
1503 if row == '':
|
|
|
1504 row = n1
|
|
|
1505 if self.save_plots_only:
|
|
|
1506 msnc.clearFig(f'theta={theta_start}')
|
|
|
1507 # center_stack order: row,theta,column
|
|
|
1508 center_offset = self._findCenterOnePlane(center_stack[row,:,:], row, thetas_deg,
|
|
|
1509 eff_pixel_size, cross_sectional_dim)
|
|
|
1510 logging.info(f'Lower center offset = {center_offset}')
|
|
|
1511
|
|
|
1512 # Update config and save to file
|
|
|
1513 find_center['row_bounds'] = [n1, n2]
|
|
|
1514 find_center['lower_row'] = row
|
|
|
1515 find_center['lower_center_offset'] = center_offset
|
|
|
1516 self.cf.saveFile(self.config_out)
|
|
|
1517 lower_row = row
|
|
|
1518
|
|
|
1519 # Upper row center
|
|
|
1520 use_row = False
|
|
|
1521 use_center = False
|
|
|
1522 row = find_center.get('upper_row')
|
|
|
1523 if msnc.is_int(row, lower_row+1, n2-1):
|
|
|
1524 msnc.quickImshow(center_stack[:,0,:], title=f'theta={theta_start}', aspect='auto')
|
|
|
1525 use_row = pyip.inputYesNo('\nCurrent row index for upper center = '
|
|
|
1526 f'{row}, use this value (y/n)? ')
|
|
|
1527 if self.save_plots_only:
|
|
|
1528 msnc.clearFig(f'theta={theta_start}')
|
|
|
1529 if use_row:
|
|
|
1530 center_offset = find_center.get('upper_center_offset')
|
|
|
1531 if msnc.is_num(center_offset):
|
|
|
1532 use_center = pyip.inputYesNo('Current upper center offset = '+
|
|
|
1533 f'{center_offset}, use this value (y/n)? ')
|
|
|
1534 if not use_center:
|
|
|
1535 if not use_row:
|
|
|
1536 msnc.quickImshow(center_stack[:,0,:], title=f'theta={theta_start}', aspect='auto')
|
|
|
1537 row = pyip.inputInt('\nEnter row index to find upper center '+
|
|
|
1538 f'[{lower_row+1}, [{n2-1}]]: ', min=lower_row+1, max=n2-1, blank=True)
|
|
|
1539 if row == '':
|
|
|
1540 row = n2-1
|
|
|
1541 if self.save_plots_only:
|
|
|
1542 msnc.clearFig(f'theta={theta_start}')
|
|
|
1543 # center_stack order: row,theta,column
|
|
|
1544 center_offset = self._findCenterOnePlane(center_stack[row,:,:], row, thetas_deg,
|
|
|
1545 eff_pixel_size, cross_sectional_dim)
|
|
|
1546 logging.info(f'upper_center_offset = {center_offset}')
|
|
|
1547 del center_stack
|
|
|
1548
|
|
|
1549 # Update config and save to file
|
|
|
1550 find_center['upper_row'] = row
|
|
|
1551 find_center['upper_center_offset'] = center_offset
|
|
|
1552 find_center['completed'] = True
|
|
|
1553 self.cf.saveFile(self.config_out)
|
|
|
1554
|
|
|
1555 def checkCenters(self):
|
|
|
1556 """Check centers for the tomography stacks.
|
|
|
1557 """
|
|
|
1558 #RV TODO load all stacks and check at all stack boundaries
|
|
|
1559 return
|
|
|
1560 logging.debug('Check centers for tomography stacks')
|
|
|
1561 center_stack_index = self.config.get('center_stack_index')
|
|
|
1562 if center_stack_index is None:
|
|
|
1563 raise ValueError('Unable to read center_stack_index from config')
|
|
|
1564 center_stack_index = self.tomo_stacks[self.tomo_data_indices.index(center_stack_index)]
|
|
|
1565 lower_row = self.config.get('lower_row')
|
|
|
1566 if lower_row is None:
|
|
|
1567 raise ValueError('Unable to read lower_row from config')
|
|
|
1568 lower_center_offset = self.config.get('lower_center_offset')
|
|
|
1569 if lower_center_offset is None:
|
|
|
1570 raise ValueError('Unable to read lower_center_offset from config')
|
|
|
1571 upper_row = self.config.get('upper_row')
|
|
|
1572 if upper_row is None:
|
|
|
1573 raise ValueError('Unable to read upper_row from config')
|
|
|
1574 upper_center_offset = self.config.get('upper_center_offset')
|
|
|
1575 if upper_center_offset is None:
|
|
|
1576 raise ValueError('Unable to read upper_center_offset from config')
|
|
|
1577 center_slope = (upper_center_offset-lower_center_offset)/(upper_row-lower_row)
|
|
|
1578 shift = upper_center_offset-lower_center_offset
|
|
|
1579 if lower_row == 0:
|
|
|
1580 logging.warning(f'lower_row == 0: one row offset between both planes')
|
|
|
1581 else:
|
|
|
1582 lower_row -= 1
|
|
|
1583 lower_center_offset -= center_slope
|
|
|
1584
|
|
|
1585 # stack order: stack,row,theta,column
|
|
|
1586 if center_stack_index:
|
|
|
1587 stack1 = self.tomo_stacks[center_stack_index-1]
|
|
|
1588 stack2 = self.tomo_stacks[center_stack_index]
|
|
|
1589 if not stack1.size:
|
|
|
1590 logging.error(f'Unable to load required tomography stack {stack1}')
|
|
|
1591 elif not stack2.size:
|
|
|
1592 logging.error(f'Unable to load required tomography stack {stack1}')
|
|
|
1593 else:
|
|
|
1594 assert(0 <= lower_row < stack2.shape[0])
|
|
|
1595 assert(0 <= upper_row < stack1.shape[0])
|
|
|
1596 plane1 = stack1[upper_row,:]
|
|
|
1597 plane2 = stack2[lower_row,:]
|
|
|
1598 for i in range(-2, 3):
|
|
|
1599 shift_i = shift+2*i
|
|
|
1600 plane1_shifted = spi.shift(plane2, [0, shift_i])
|
|
|
1601 msnc.quickPlot((plane1[0,:],), (plane1_shifted[0,:],),
|
|
|
1602 title=f'stacks {stack1} {stack2} shifted {2*i} theta={self.start_theta}',
|
|
|
1603 path=self.output_folder, save_fig=self.save_plots,
|
|
|
1604 save_only=self.save_plots_only)
|
|
|
1605 if center_stack_index < self.num_tomo_stacks-1:
|
|
|
1606 stack1 = self.tomo_stacks[center_stack_index]
|
|
|
1607 stack2 = self.tomo_stacks[center_stack_index+1]
|
|
|
1608 if not stack1.size:
|
|
|
1609 logging.error(f'Unable to load required tomography stack {stack1}')
|
|
|
1610 elif not stack2.size:
|
|
|
1611 logging.error(f'Unable to load required tomography stack {stack1}')
|
|
|
1612 else:
|
|
|
1613 assert(0 <= lower_row < stack2.shape[0])
|
|
|
1614 assert(0 <= upper_row < stack1.shape[0])
|
|
|
1615 plane1 = stack1[upper_row,:]
|
|
|
1616 plane2 = stack2[lower_row,:]
|
|
|
1617 for i in range(-2, 3):
|
|
|
1618 shift_i = -shift+2*i
|
|
|
1619 plane1_shifted = spi.shift(plane2, [0, shift_i])
|
|
|
1620 msnc.quickPlot((plane1[0,:],), (plane1_shifted[0,:],),
|
|
|
1621 title=f'stacks {stack1} {stack2} shifted {2*i} theta={start_theta}',
|
|
|
1622 path=self.output_folder, save_fig=self.save_plots,
|
|
|
1623 save_only=self.save_plots_only)
|
|
|
1624 del plane1, plane2, plane1_shifted
|
|
|
1625
|
|
|
1626 # Update config file
|
|
|
1627 self.config = msnc.update('config.txt', 'check_centers', True, 'find_centers')
|
|
|
1628
|
|
|
1629 def reconstructTomoStacks(self):
|
|
|
1630 """Reconstruct tomography stacks.
|
|
|
1631 """
|
|
|
1632 logging.debug('Reconstruct tomography stacks')
|
|
|
1633
|
|
|
1634 # Get rotation axis rows and centers
|
|
|
1635 find_center = self.config['find_center']
|
|
|
1636 lower_row = find_center.get('lower_row')
|
|
|
1637 if lower_row is None:
|
|
|
1638 logging.error('Unable to read lower_row from config')
|
|
|
1639 return
|
|
|
1640 lower_center_offset = find_center.get('lower_center_offset')
|
|
|
1641 if lower_center_offset is None:
|
|
|
1642 logging.error('Unable to read lower_center_offset from config')
|
|
|
1643 return
|
|
|
1644 upper_row = find_center.get('upper_row')
|
|
|
1645 if upper_row is None:
|
|
|
1646 logging.error('Unable to read upper_row from config')
|
|
|
1647 return
|
|
|
1648 upper_center_offset = find_center.get('upper_center_offset')
|
|
|
1649 if upper_center_offset is None:
|
|
|
1650 logging.error('Unable to read upper_center_offset from config')
|
|
|
1651 return
|
|
|
1652 logging.debug(f'lower_row = {lower_row} upper_row = {upper_row}')
|
|
|
1653 assert(lower_row < upper_row)
|
|
|
1654 center_slope = (upper_center_offset-lower_center_offset)/(upper_row-lower_row)
|
|
|
1655
|
|
|
1656 # Set thetas (in radians)
|
|
|
1657 theta_range = self.config['theta_range']
|
|
|
1658 theta_start = theta_range['start']
|
|
|
1659 theta_end = theta_range['end']
|
|
|
1660 num_theta = theta_range['num']
|
|
|
1661 num_theta_skip = self.config['preprocess'].get('num_theta_skip', 0)
|
|
|
1662 thetas = np.radians(np.linspace(theta_start, theta_end,
|
|
|
1663 int(num_theta/(num_theta_skip+1)), endpoint=False))
|
|
|
1664
|
|
|
1665 # Reconstruct tomo stacks
|
|
|
1666 zoom_perc = self.config['preprocess'].get('zoom_perc', 100)
|
|
|
1667 if zoom_perc == 100:
|
|
|
1668 basetitle = 'recon stack full'
|
|
|
1669 else:
|
|
|
1670 basetitle = f'recon stack {zoom_perc}p'
|
|
|
1671 load_error = False
|
|
|
1672 stacks = self.config['stack_info']['stacks']
|
|
|
1673 for i,stack in enumerate(stacks):
|
|
|
1674 # Check if stack can be loaded
|
|
|
1675 # reconstructed stack order for each one in stack : row/z,x,y
|
|
|
1676 # preprocessed stack order for each one in stack: row,theta,column
|
|
|
1677 index = stack['index']
|
|
|
1678 available = False
|
|
|
1679 if stack.get('reconstructed', False):
|
|
|
1680 self.tomo_recon_stacks[i], available = self._loadTomo('recon stack', index)
|
|
|
1681 if self.tomo_recon_stacks[i].size or available:
|
|
|
1682 if self.tomo_stacks[i].size:
|
|
|
1683 self.tomo_stacks[i] = np.array([])
|
|
|
1684 assert(stack.get('reconstructed', False) == True)
|
|
|
1685 continue
|
|
|
1686 else:
|
|
|
1687 stack['reconstructed'] = False
|
|
|
1688 if not self.tomo_stacks[i].size:
|
|
|
1689 self.tomo_stacks[i], available = self._loadTomo('red stack', index,
|
|
|
1690 required=True)
|
|
|
1691 if not self.tomo_stacks[i].size:
|
|
|
1692 logging.error(f'Unable to load tomography stack {index} for reconstruction')
|
|
|
1693 load_error = True
|
|
|
1694 continue
|
|
|
1695 assert(0 <= lower_row < upper_row < self.tomo_stacks[i].shape[0])
|
|
|
1696 center_offsets = [lower_center_offset-lower_row*center_slope,
|
|
|
1697 upper_center_offset+(self.tomo_stacks[i].shape[0]-1-upper_row)*center_slope]
|
|
|
1698 t0 = time()
|
|
|
1699 self.tomo_recon_stacks[i]= self._reconstructOneTomoStack(self.tomo_stacks[i],
|
|
|
1700 thetas, center_offsets=center_offsets, sigma=0.1, ncore=self.ncore,
|
|
|
1701 algorithm='gridrec', run_secondary_sirt=True, secondary_iter=25)
|
|
|
1702 logging.info(f'Reconstruction of stack {index} took {time()-t0:.2f} seconds!')
|
|
|
1703 if not self.test_mode:
|
|
|
1704 row_slice = int(self.tomo_stacks[i].shape[0]/2)
|
|
|
1705 title = f'{basetitle} {index} slice{row_slice}'
|
|
|
1706 msnc.quickImshow(self.tomo_recon_stacks[i][row_slice,:,:], title=title,
|
|
|
1707 path=self.output_folder, save_fig=self.save_plots,
|
|
|
1708 save_only=self.save_plots_only)
|
|
|
1709 msnc.quickPlot(self.tomo_recon_stacks[i]
|
|
|
1710 [row_slice,int(self.tomo_recon_stacks[i].shape[2]/2),:],
|
|
|
1711 title=f'{title} cut{int(self.tomo_recon_stacks[i].shape[2]/2)}',
|
|
|
1712 path=self.output_folder, save_fig=self.save_plots,
|
|
|
1713 save_only=self.save_plots_only)
|
|
|
1714 self._saveTomo('recon stack', self.tomo_recon_stacks[i], index)
|
|
|
1715 # else:
|
|
|
1716 # np.savetxt(self.output_folder+f'recon_stack_{index}.txt',
|
|
|
1717 # self.tomo_recon_stacks[i][row_slice,:,:], fmt='%.6e')
|
|
|
1718 self.tomo_stacks[i] = np.array([])
|
|
|
1719
|
|
|
1720 # Update config and save to file
|
|
|
1721 stack['reconstructed'] = True
|
|
|
1722 self.cf.saveFile(self.config_out)
|
|
|
1723
|
|
|
1724 def combineTomoStacks(self):
|
|
|
1725 """Combine the reconstructed tomography stacks.
|
|
|
1726 """
|
|
|
1727 # stack order: stack,row(z),x,y
|
|
|
1728 logging.debug('Combine reconstructed tomography stacks')
|
|
|
1729 # Load any unloaded reconstructed stacks
|
|
|
1730 stacks = self.config['stack_info']['stacks']
|
|
|
1731 for i,stack in enumerate(stacks):
|
|
|
1732 if not self.tomo_recon_stacks[i].size and stack.get('reconstructed', False):
|
|
|
1733 self.tomo_recon_stacks[i], available = self._loadTomo('recon stack',
|
|
|
1734 stack['index'], required=True)
|
|
|
1735 if not self.tomo_recon_stacks[i].size or not available:
|
|
|
1736 logging.error(f'Unable to load reconstructed stack {stack["index"]}')
|
|
|
1737 return
|
|
|
1738 if i:
|
|
|
1739 if (self.tomo_recon_stacks[i].shape[1] != self.tomo_recon_stacks[0].shape[1] or
|
|
|
1740 self.tomo_recon_stacks[i].shape[1] != self.tomo_recon_stacks[0].shape[1]):
|
|
|
1741 logging.error('Incompatible reconstructed tomography stack dimensions')
|
|
|
1742 return
|
|
|
1743
|
|
|
1744 # Get center stack boundaries
|
|
|
1745 row_bounds = self.config['find_center']['row_bounds']
|
|
|
1746 if not msnc.is_index_range(row_bounds, 0, self.tomo_recon_stacks[0].shape[0]):
|
|
|
1747 msnc.illegal_value('row_bounds', row_bounds, 'config file')
|
|
|
1748 return
|
|
|
1749
|
|
|
1750 # Selecting xy bounds
|
|
|
1751 tomosum = 0
|
|
|
1752 #RV FIX :=
|
|
|
1753 #[tomosum := tomosum+np.sum(tomo_recon_stack, axis=(0,2)) for tomo_recon_stack in
|
|
|
1754 # self.tomo_recon_stacks]
|
|
|
1755 combine_stacks =self.config.get('combine_stacks')
|
|
|
1756 if combine_stacks and 'x_bounds' in combine_stacks:
|
|
|
1757 x_bounds = combine_stacks['x_bounds']
|
|
|
1758 if not msnc.is_index_range(x_bounds, 0, self.tomo_recon_stacks[0].shape[1]):
|
|
|
1759 msnc.illegal_value('x_bounds', x_bounds, 'config file')
|
|
|
1760 elif not self.test_mode:
|
|
|
1761 msnc.quickPlot(tomosum, title='recon stack sum yz')
|
|
|
1762 if pyip.inputYesNo('\nCurrent image x-bounds: '+
|
|
|
1763 f'[{x_bounds[0]}, {x_bounds[1]}], use these values ([y]/n)? ',
|
|
|
1764 blank=True) == 'no':
|
|
|
1765 if pyip.inputYesNo(
|
|
|
1766 'Do you want to change the image x-bounds ([y]/n)? ',
|
|
|
1767 blank=True) == 'no':
|
|
|
1768 x_bounds = [0, self.tomo_recon_stacks[0].shape[1]]
|
|
|
1769 else:
|
|
|
1770 x_bounds = msnc.selectArrayBounds(tomosum, title='recon stack sum yz')
|
|
|
1771 else:
|
|
|
1772 msnc.quickPlot(tomosum, title='recon stack sum yz')
|
|
|
1773 if pyip.inputYesNo('\nDo you want to change the image x-bounds (y/n)? ') == 'no':
|
|
|
1774 x_bounds = [0, self.tomo_recon_stacks[0].shape[1]]
|
|
|
1775 else:
|
|
|
1776 x_bounds = msnc.selectArrayBounds(tomosum, title='recon stack sum yz')
|
|
|
1777 if False and self.test_mode:
|
|
|
1778 np.savetxt(self.output_folder+'recon_stack_sum_yz.txt',
|
|
|
1779 tomosum[x_bounds[0]:x_bounds[1]], fmt='%.6e')
|
|
|
1780 if self.save_plots_only:
|
|
|
1781 msnc.clearFig('recon stack sum yz')
|
|
|
1782 tomosum = 0
|
|
|
1783 #RV FIX :=
|
|
|
1784 #[tomosum := tomosum+np.sum(tomo_recon_stack, axis=(0,1)) for tomo_recon_stack in
|
|
|
1785 # self.tomo_recon_stacks]
|
|
|
1786 if combine_stacks and 'y_bounds' in combine_stacks:
|
|
|
1787 y_bounds = combine_stacks['y_bounds']
|
|
|
1788 if not msnc.is_index_range(x_bounds, 0, self.tomo_recon_stacks[0].shape[1]):
|
|
|
1789 msnc.illegal_value('y_bounds', y_bounds, 'config file')
|
|
|
1790 elif not self.test_mode:
|
|
|
1791 msnc.quickPlot(tomosum, title='recon stack sum xz')
|
|
|
1792 if pyip.inputYesNo('\nCurrent image y-bounds: '+
|
|
|
1793 f'[{y_bounds[0]}, {y_bounds[1]}], use these values ([y]/n)? ',
|
|
|
1794 blank=True) == 'no':
|
|
|
1795 if pyip.inputYesNo(
|
|
|
1796 'Do you want to change the image y-bounds ([y]/n)? ',
|
|
|
1797 blank=True) == 'no':
|
|
|
1798 y_bounds = [0, self.tomo_recon_stacks[0].shape[1]]
|
|
|
1799 else:
|
|
|
1800 y_bounds = msnc.selectArrayBounds(tomosum, title='recon stack sum yz')
|
|
|
1801 else:
|
|
|
1802 msnc.quickPlot(tomosum, title='recon stack sum xz')
|
|
|
1803 if pyip.inputYesNo('\nDo you want to change the image y-bounds (y/n)? ') == 'no':
|
|
|
1804 y_bounds = [0, self.tomo_recon_stacks[0].shape[1]]
|
|
|
1805 else:
|
|
|
1806 y_bounds = msnc.selectArrayBounds(tomosum, title='recon stack sum xz')
|
|
|
1807 if False and self.test_mode:
|
|
|
1808 np.savetxt(self.output_folder+'recon_stack_sum_xz.txt',
|
|
|
1809 tomosum[y_bounds[0]:y_bounds[1]], fmt='%.6e')
|
|
|
1810 if self.save_plots_only:
|
|
|
1811 msnc.clearFig('recon stack sum xz')
|
|
|
1812
|
|
|
1813 # Combine reconstructed tomography stacks
|
|
|
1814 logging.info(f'Combining reconstructed stacks ...')
|
|
|
1815 t0 = time()
|
|
|
1816 num_tomo_stacks = self.config['stack_info']['num']
|
|
|
1817 if num_tomo_stacks == 1:
|
|
|
1818 low_bound = row_bounds[0]
|
|
|
1819 else:
|
|
|
1820 low_bound = 0
|
|
|
1821 tomo_recon_combined = self.tomo_recon_stacks[0][low_bound:row_bounds[1]:,
|
|
|
1822 x_bounds[0]:x_bounds[1],y_bounds[0]:y_bounds[1]]
|
|
|
1823 if num_tomo_stacks > 2:
|
|
|
1824 tomo_recon_combined = np.concatenate([tomo_recon_combined]+
|
|
|
1825 [self.tomo_recon_stacks[i][row_bounds[0]:row_bounds[1],
|
|
|
1826 x_bounds[0]:x_bounds[1],y_bounds[0]:y_bounds[1]]
|
|
|
1827 for i in range(1, num_tomo_stacks-1)])
|
|
|
1828 if num_tomo_stacks > 1:
|
|
|
1829 tomo_recon_combined = np.concatenate([tomo_recon_combined]+
|
|
|
1830 [self.tomo_recon_stacks[num_tomo_stacks-1][row_bounds[0]:,
|
|
|
1831 x_bounds[0]:x_bounds[1],y_bounds[0]:y_bounds[1]]])
|
|
|
1832 logging.info(f'... done in {time()-t0:.2f} seconds!')
|
|
|
1833 tomosum = np.sum(tomo_recon_combined, axis=(1,2))
|
|
|
1834 if self.test_mode:
|
|
|
1835 zoom_perc = self.config['preprocess'].get('zoom_perc', 100)
|
|
|
1836 filename = 'recon combined sum xy'
|
|
|
1837 if zoom_perc is None or zoom_perc == 100:
|
|
|
1838 filename += ' fullres.dat'
|
|
|
1839 else:
|
|
|
1840 filename += f' {zoom_perc}p.dat'
|
|
|
1841 msnc.quickPlot(tomosum, title='recon combined sum xy',
|
|
|
1842 path=self.output_folder, save_fig=self.save_plots,
|
|
|
1843 save_only=self.save_plots_only)
|
|
|
1844 if False:
|
|
|
1845 np.savetxt(self.output_folder+'recon_combined_sum_xy.txt',
|
|
|
1846 tomosum, fmt='%.6e')
|
|
|
1847 np.savetxt(self.output_folder+'recon_combined.txt',
|
|
|
1848 tomo_recon_combined[int(tomo_recon_combined.shape[0]/2),:,:], fmt='%.6e')
|
|
|
1849 combine_stacks =self.config.get('combine_stacks')
|
|
|
1850
|
|
|
1851 # Update config and save to file
|
|
|
1852 if combine_stacks:
|
|
|
1853 combine_stacks['x_bounds'] = x_bounds
|
|
|
1854 combine_stacks['y_bounds'] = y_bounds
|
|
|
1855 else:
|
|
|
1856 self.config['combine_stacks'] = {'x_bounds' : x_bounds, 'y_bounds' : y_bounds}
|
|
|
1857 self.cf.saveFile(self.config_out)
|
|
|
1858 return
|
|
|
1859 msnc.quickPlot(tomosum, title='recon combined sum xy')
|
|
|
1860 if pyip.inputYesNo(
|
|
|
1861 '\nDo you want to change the image z-bounds (y/[n])? ',
|
|
|
1862 blank=True) != 'yes':
|
|
|
1863 z_bounds = [0, tomo_recon_combined.shape[0]]
|
|
|
1864 else:
|
|
|
1865 z_bounds = msnc.selectArrayBounds(tomosum, title='recon combined sum xy')
|
|
|
1866 if z_bounds[0] != 0 or z_bounds[1] != tomo_recon_combined.shape[0]:
|
|
|
1867 tomo_recon_combined = tomo_recon_combined[z_bounds[0]:z_bounds[1],:,:]
|
|
|
1868 logging.info(f'tomo_recon_combined.shape = {tomo_recon_combined.shape}')
|
|
|
1869 if self.save_plots_only:
|
|
|
1870 msnc.clearFig('recon combined sum xy')
|
|
|
1871
|
|
|
1872 # Plot center slices
|
|
|
1873 msnc.quickImshow(tomo_recon_combined[int(tomo_recon_combined.shape[0]/2),:,:],
|
|
|
1874 title=f'recon combined xslice{int(tomo_recon_combined.shape[0]/2)}',
|
|
|
1875 path=self.output_folder, save_fig=self.save_plots,
|
|
|
1876 save_only=self.save_plots_only)
|
|
|
1877 msnc.quickImshow(tomo_recon_combined[:,int(tomo_recon_combined.shape[1]/2),:],
|
|
|
1878 title=f'recon combined yslice{int(tomo_recon_combined.shape[1]/2)}',
|
|
|
1879 path=self.output_folder, save_fig=self.save_plots,
|
|
|
1880 save_only=self.save_plots_only)
|
|
|
1881 msnc.quickImshow(tomo_recon_combined[:,:,int(tomo_recon_combined.shape[2]/2)],
|
|
|
1882 title=f'recon combined zslice{int(tomo_recon_combined.shape[2]/2)}',
|
|
|
1883 path=self.output_folder, save_fig=self.save_plots,
|
|
|
1884 save_only=self.save_plots_only)
|
|
|
1885
|
|
|
1886 # Save combined reconstructed tomo stacks
|
|
|
1887 base_name = 'recon combined'
|
|
|
1888 combined_stacks = []
|
|
|
1889 for stack in stacks:
|
|
|
1890 base_name += f' {stack["index"]}'
|
|
|
1891 combined_stacks.append(stack['index'])
|
|
|
1892 self._saveTomo(base_name, tomo_recon_combined)
|
|
|
1893
|
|
|
1894 # Update config and save to file
|
|
|
1895 if combine_stacks:
|
|
|
1896 combine_stacks['x_bounds'] = x_bounds
|
|
|
1897 combine_stacks['y_bounds'] = y_bounds
|
|
|
1898 combine_stacks['stacks'] = combined_stacks
|
|
|
1899 else:
|
|
|
1900 self.config['combine_stacks'] = {'x_bounds' : x_bounds, 'y_bounds' : y_bounds,
|
|
|
1901 'stacks' : combined_stacks}
|
|
|
1902 self.cf.saveFile(self.config_out)
|
|
|
1903
|
|
|
1904 def runTomo(config_file=None, config_dict=None, output_folder='.', log_level='INFO',
|
|
|
1905 test_mode=False):
|
|
|
1906 """Run a tomography analysis.
|
|
|
1907 """
|
|
|
1908 # Instantiate Tomo object
|
|
|
1909 tomo = Tomo(config_file=config_file, output_folder=output_folder, log_level=log_level,
|
|
|
1910 test_mode=test_mode)
|
|
|
1911 if not tomo.is_valid:
|
|
|
1912 raise ValueError('Invalid config and/or detector file provided.')
|
|
|
1913
|
|
|
1914 # Preprocess the image files
|
|
|
1915 num_tomo_stacks = tomo.config['stack_info']['num']
|
|
|
1916 assert(num_tomo_stacks == len(tomo.tomo_stacks))
|
|
|
1917 preprocess = tomo.config.get('preprocess', None)
|
|
|
1918 preprocessed_stacks = []
|
|
|
1919 if preprocess:
|
|
|
1920 preprocessed_stacks = [stack['index'] for stack in tomo.config['stack_info']['stacks']
|
|
|
1921 if stack.get('preprocessed', False)]
|
|
|
1922 if len(preprocessed_stacks):
|
|
|
1923 tomo.genTomoStacks()
|
|
|
1924 if not tomo.is_valid:
|
|
|
1925 IOError('Unable to load all required image files.')
|
|
|
1926 find_center = tomo.config.get('find_center')
|
|
|
1927 if find_center and find_center.get('completed', False):
|
|
|
1928 center_stack_index = find_center['center_stack_index']
|
|
|
1929 if not center_stack_index in preprocessed_stacks:
|
|
|
1930 find_center['completed'] = False
|
|
|
1931 #RV FIX
|
|
|
1932 # tomo.config.pop('check_center', 'check_center not found')
|
|
|
1933 # combined_stacks = tomo.config.get('combined_stacks')
|
|
|
1934 # if combined_stacks:
|
|
|
1935 # combined_stacks['completed'] = False
|
|
|
1936 tomo.cf.saveFile(self.config_out)
|
|
|
1937
|
|
|
1938 # Find centers
|
|
|
1939 find_center = tomo.config.get('find_center')
|
|
|
1940 if find_center is None or not find_center.get('completed', False):
|
|
|
1941 tomo.findCenters()
|
|
|
1942
|
|
|
1943 # Check centers
|
|
|
1944 #if num_tomo_stacks > 1 and not tomo.config.get('check_centers', False):
|
|
|
1945 # tomo.checkCenters()
|
|
|
1946
|
|
|
1947 # Reconstruct tomography stacks
|
|
|
1948 if len(tomo.config.get('reconstructed_stacks', [])) != tomo.config['stack_info']['num']:
|
|
|
1949 tomo.reconstructTomoStacks()
|
|
|
1950
|
|
|
1951 # Combine reconstructed tomography stacks
|
|
|
1952 combined_stacks = tomo.config.get('combined_stacks')
|
|
|
1953 if combined_stacks is None or not combined_stacks.get('completed', False):
|
|
|
1954 tomo.combineTomoStacks()
|
|
|
1955
|
|
|
1956 #%%============================================================================
|
|
|
1957 if __name__ == '__main__':
|
|
|
1958 # Parse command line arguments
|
|
|
1959 arguments = sys.argv[1:]
|
|
|
1960 config_file = None
|
|
|
1961 output_folder = '.'
|
|
|
1962 log_level = 'INFO'
|
|
|
1963 test_mode = False
|
|
|
1964 try:
|
|
|
1965 opts, args = getopt.getopt(arguments,"hc:o:l:t")
|
|
|
1966 except getopt.GetoptError:
|
|
|
1967 print('usage: tomo.py -c <config_file> -o <output_folder> -l <log_level> -t')
|
|
|
1968 sys.exit(2)
|
|
|
1969 for opt, arg in opts:
|
|
|
1970 if opt == '-h':
|
|
|
1971 print('usage: tomo.py -c <config_file> -o <output_folder> -l <log_level> -t')
|
|
|
1972 sys.exit()
|
|
|
1973 elif opt in ("-c"):
|
|
|
1974 config_file = arg
|
|
|
1975 elif opt in ("-o"):
|
|
|
1976 output_folder = arg
|
|
|
1977 elif opt in ("-l"):
|
|
|
1978 log_level = arg
|
|
|
1979 elif opt in ("-t"):
|
|
|
1980 test_mode = True
|
|
|
1981 if config_file is None:
|
|
|
1982 if os.path.isfile('config.yaml'):
|
|
|
1983 config_file = 'config.yaml'
|
|
|
1984 else:
|
|
|
1985 config_file = 'config.txt'
|
|
|
1986
|
|
|
1987 # Set basic log configuration
|
|
|
1988 logging_format = '%(asctime)s : %(levelname)s - %(module)s : %(funcName)s - %(message)s'
|
|
|
1989 if not test_mode:
|
|
|
1990 level = getattr(logging, log_level.upper(), None)
|
|
|
1991 if not isinstance(level, int):
|
|
|
1992 raise ValueError(f'Invalid log_level: {log_level}')
|
|
|
1993 logging.basicConfig(format=logging_format, level=level, force=True,
|
|
|
1994 handlers=[logging.StreamHandler()])
|
|
|
1995
|
|
|
1996 logging.debug(f'config_file = {config_file}')
|
|
|
1997 logging.debug(f'output_folder = {output_folder}')
|
|
|
1998 logging.debug(f'log_level = {log_level}')
|
|
|
1999 logging.debug(f'test_mode = {test_mode}')
|
|
|
2000
|
|
|
2001 # Run tomography analysis
|
|
|
2002 runTomo(config_file=config_file, output_folder=output_folder, log_level=log_level,
|
|
|
2003 test_mode=test_mode)
|
|
|
2004
|
|
|
2005 #%%============================================================================
|
|
|
2006 input('Press any key to continue')
|
|
|
2007 #%%============================================================================
|
