Mercurial > repos > rv43 > tomo

comparison tomo.py @ 21:3caba2116858 draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
"planemo upload for repository https://github.com/rolfverberg/galaxytools commit e32e37ee46f47b185f924884685eb1d464761e86"
author rv43
date Mon, 18 Apr 2022 19:43:15 +0000
parents c033f8843dc0
children ac2f726f9054
comparison
equal deleted inserted replaced
20:c033f8843dc0 21:3caba2116858
1881 self.config = msnc.update('config.txt', 'check_centers', True, 'find_centers') 1881 self.config = msnc.update('config.txt', 'check_centers', True, 'find_centers')
1882 1882
1883 def reconstructTomoStacks(self, galaxy_param=None, num_core=None): 1883 def reconstructTomoStacks(self, galaxy_param=None, num_core=None):
1884 """Reconstruct tomography stacks. 1884 """Reconstruct tomography stacks.
1885 """ 1885 """
1886 print('OK1')
1887 if num_core is None: 1886 if num_core is None:
1888 num_core = self.num_core 1887 num_core = self.num_core
1888 if self.galaxy_flag:
1889 assert(galaxy_param)
1890 if not os.path.exists('center_slice_pngs'):
1891 os.mkdir('center_slice_pngs')
1889 logging.debug('Reconstruct tomography stacks') 1892 logging.debug('Reconstruct tomography stacks')
1890 stacks = self.config['stack_info']['stacks'] 1893 stacks = self.config['stack_info']['stacks']
1891 assert(len(self.tomo_stacks) == self.config['stack_info']['num']) 1894 assert(len(self.tomo_stacks) == self.config['stack_info']['num'])
1892 assert(len(self.tomo_stacks) == len(stacks)) 1895 assert(len(self.tomo_stacks) == len(stacks))
1893 assert(len(self.tomo_recon_stacks) == len(stacks)) 1896 assert(len(self.tomo_recon_stacks) == len(stacks))
1894 print('OK2')
1895 if self.galaxy_flag: 1897 if self.galaxy_flag:
1896 assert(isinstance(galaxy_param, dict)) 1898 assert(isinstance(galaxy_param, dict))
1897 # Get rotation axis centers 1899 # Get rotation axis centers
1898 center_offsets = galaxy_param['center_offsets'] 1900 center_offsets = galaxy_param['center_offsets']
1899 assert(isinstance(center_offsets, list) and len(center_offsets) == 2) 1901 assert(isinstance(center_offsets, list) and len(center_offsets) == 2)
1905 if galaxy_param: 1907 if galaxy_param:
1906 logging.warning('Ignoring galaxy_param in reconstructTomoStacks (only for Galaxy)') 1908 logging.warning('Ignoring galaxy_param in reconstructTomoStacks (only for Galaxy)')
1907 galaxy_param = None 1909 galaxy_param = None
1908 lower_center_offset = None 1910 lower_center_offset = None
1909 upper_center_offset = None 1911 upper_center_offset = None
1910 print('OK3')
1911 1912
1912 # Get rotation axis rows and centers 1913 # Get rotation axis rows and centers
1913 find_center = self.config['find_center'] 1914 find_center = self.config['find_center']
1914 lower_row = find_center.get('lower_row') 1915 lower_row = find_center.get('lower_row')
1915 if lower_row is None: 1916 if lower_row is None:
1932 logging.error('Unable to read upper_center_offset from config') 1933 logging.error('Unable to read upper_center_offset from config')
1933 return 1934 return
1934 center_slope = (upper_center_offset-lower_center_offset)/(upper_row-lower_row) 1935 center_slope = (upper_center_offset-lower_center_offset)/(upper_row-lower_row)
1935 1936
1936 # Set thetas (in radians) 1937 # Set thetas (in radians)
1937 print('OK4')
1938 theta_range = self.config['theta_range'] 1938 theta_range = self.config['theta_range']
1939 theta_start = theta_range['start'] 1939 theta_start = theta_range['start']
1940 theta_end = theta_range['end'] 1940 theta_end = theta_range['end']
1941 num_theta = theta_range['num'] 1941 num_theta = theta_range['num']
1942 num_theta_skip = self.config['preprocess'].get('num_theta_skip', 0) 1942 num_theta_skip = self.config['preprocess'].get('num_theta_skip', 0)
1943 thetas = np.radians(np.linspace(theta_start, theta_end, 1943 thetas = np.radians(np.linspace(theta_start, theta_end,
1944 int(num_theta/(num_theta_skip+1)), endpoint=False)) 1944 int(num_theta/(num_theta_skip+1)), endpoint=False))
1945 1945
1946 print('OK5')
1947 # Reconstruct tomo stacks 1946 # Reconstruct tomo stacks
1948 zoom_perc = self.config['preprocess'].get('zoom_perc', 100) 1947 zoom_perc = self.config['preprocess'].get('zoom_perc', 100)
1949 if zoom_perc == 100: 1948 if zoom_perc == 100:
1950 basetitle = 'recon stack full' 1949 basetitle = 'recon stack fullres'
1951 else: 1950 else:
1952 basetitle = f'recon stack {zoom_perc}p' 1951 basetitle = f'recon stack {zoom_perc}p'
1953 load_error = False 1952 load_error = False
1954 stacks = self.config['stack_info']['stacks'] 1953 stacks = self.config['stack_info']['stacks']
1955 for i,stack in enumerate(stacks): 1954 for i,stack in enumerate(stacks):
1969 continue 1968 continue
1970 stack['reconstructed'] = False 1969 stack['reconstructed'] = False
1971 if not self.tomo_stacks[i].size: 1970 if not self.tomo_stacks[i].size:
1972 self.tomo_stacks[i], available = self._loadTomo('red stack', index, 1971 self.tomo_stacks[i], available = self._loadTomo('red stack', index,
1973 required=True) 1972 required=True)
1974 print(f'self.tomo_stacks.shape = {self.tomo_stacks[i].shape}')
1975 if not self.tomo_stacks[i].size: 1973 if not self.tomo_stacks[i].size:
1976 logging.error(f'Unable to load tomography stack {index} for reconstruction') 1974 logging.error(f'Unable to load tomography stack {index} for reconstruction')
1977 stack[i]['preprocessed'] = False 1975 stack[i]['preprocessed'] = False
1978 load_error = True 1976 load_error = True
1979 continue 1977 continue
1983 t0 = time() 1981 t0 = time()
1984 self.tomo_recon_stacks[i]= self._reconstructOneTomoStack(self.tomo_stacks[i], 1982 self.tomo_recon_stacks[i]= self._reconstructOneTomoStack(self.tomo_stacks[i],
1985 thetas, center_offsets=center_offsets, sigma=0.1, num_core=num_core, 1983 thetas, center_offsets=center_offsets, sigma=0.1, num_core=num_core,
1986 algorithm='gridrec', run_secondary_sirt=True, secondary_iter=25) 1984 algorithm='gridrec', run_secondary_sirt=True, secondary_iter=25)
1987 logging.info(f'Reconstruction of stack {index} took {time()-t0:.2f} seconds!') 1985 logging.info(f'Reconstruction of stack {index} took {time()-t0:.2f} seconds!')
1988 if not self.test_mode and not self.galaxy_flag: 1986 if self.galaxy_flag:
1989 row_slice = int(self.tomo_stacks[i].shape[0]/2) 1987 x_slice = int(self.tomo_stacks[i].shape[0]/2)
1990 title = f'{basetitle} {index} slice{row_slice}' 1988 title = f'{basetitle} {index} xslice{x_slice}'
1991 msnc.quickImshow(self.tomo_recon_stacks[i][row_slice,:,:], title=title, 1989 msnc.quickImshow(self.tomo_recon_stacks[i][x_slice,:,:], title=title,
1990 path='center_slice_pngs', save_fig=True, save_only=True)
1991 y_slice = int(self.tomo_stacks[i].shape[0]/2)
1992 title = f'{basetitle} {index} yslice{y_slice}'
1993 msnc.quickImshow(self.tomo_recon_stacks[i][:,y_slice,:], title=title,
1994 path='center_slice_pngs', save_fig=True, save_only=True)
1995 z_slice = int(self.tomo_stacks[i].shape[0]/2)
1996 title = f'{basetitle} {index} zslice{z_slice}'
1997 msnc.quickImshow(self.tomo_recon_stacks[i][:,:,z_slice], title=title,
1998 path='center_slice_pngs', save_fig=True, save_only=True)
1999 elif not self.test_mode:
2000 x_slice = int(self.tomo_stacks[i].shape[0]/2)
2001 title = f'{basetitle} {index} xslice{x_slice}'
2002 msnc.quickImshow(self.tomo_recon_stacks[i][x_slice,:,:], title=title,
1992 path=self.output_folder, save_fig=self.save_plots, 2003 path=self.output_folder, save_fig=self.save_plots,
1993 save_only=self.save_plots_only) 2004 save_only=self.save_plots_only)
1994 msnc.quickPlot(self.tomo_recon_stacks[i] 2005 msnc.quickPlot(self.tomo_recon_stacks[i]
1995 [row_slice,int(self.tomo_recon_stacks[i].shape[2]/2),:], 2006 [x_slice,int(self.tomo_recon_stacks[i].shape[2]/2),:],
1996 title=f'{title} cut{int(self.tomo_recon_stacks[i].shape[2]/2)}', 2007 title=f'{title} cut{int(self.tomo_recon_stacks[i].shape[2]/2)}',
1997 path=self.output_folder, save_fig=self.save_plots, 2008 path=self.output_folder, save_fig=self.save_plots,
1998 save_only=self.save_plots_only) 2009 save_only=self.save_plots_only)
1999 self._saveTomo('recon stack', self.tomo_recon_stacks[i], index) 2010 self._saveTomo('recon stack', self.tomo_recon_stacks[i], index)
2000 # else:
2001 # np.savetxt(self.output_folder+f'recon_stack_{index}.txt',
2002 # self.tomo_recon_stacks[i][row_slice,:,:], fmt='%.6e')
2003 self.tomo_stacks[i] = np.array([]) 2011 self.tomo_stacks[i] = np.array([])
2004 print('OK6')
2005 2012
2006 # Update config and save to file 2013 # Update config and save to file
2007 stack['reconstructed'] = True 2014 stack['reconstructed'] = True
2008 combine_stacks = self.config.get('combine_stacks') 2015 combine_stacks = self.config.get('combine_stacks')
2009 if combine_stacks and index in combine_stacks.get('stacks', []): 2016 if combine_stacks and index in combine_stacks.get('stacks', []):
2010 combine_stacks['stacks'].remove(index) 2017 combine_stacks['stacks'].remove(index)
2011 self.cf.saveFile(self.config_out) 2018 self.cf.saveFile(self.config_out)
2012 2019
2013 print('OK7')
2014 # Save reconstructed tomography stack to file 2020 # Save reconstructed tomography stack to file
2015 if self.galaxy_flag: 2021 if self.galaxy_flag:
2016 t0 = time() 2022 t0 = time()
2017 output_name = galaxy_param['output_name'] 2023 output_name = galaxy_param['output_name']
2018 logging.info(f'Saving reconstructed tomography stack to {output_name} ...') 2024 logging.info(f'Saving reconstructed tomography stack to {output_name} ...')