changeset 3:f9c52762c32c draft

"planemo upload for repository https://github.com/rolfverberg/galaxytools commit 7dce44d576e4149f31bdc2ee4dce0bb6962badb6"
author rv43
date Tue, 05 Apr 2022 18:23:54 +0000
parents b8977c98800b
children f221b7ba4829
files __pycache__/msnc_tools.cpython-39.pyc __pycache__/tomo.cpython-39.pyc msnc_tools.py run_tomo_find_center run_tomo_reconstruct run_tomo_setup tomo.py tomo_find_center.py tomo_find_center.xml tomo_reconstruct.py tomo_reconstruct.xml tomo_setup.xml
diffstat 12 files changed, 514 insertions(+), 145 deletions(-) [+]
line wrap: on
line diff
Binary file __pycache__/msnc_tools.cpython-39.pyc has changed
Binary file __pycache__/tomo.cpython-39.pyc has changed
--- a/msnc_tools.py	Thu Mar 31 20:48:17 2022 +0000
+++ b/msnc_tools.py	Tue Apr 05 18:23:54 2022 +0000
@@ -297,7 +297,7 @@
     plt.close(fig=re.sub(r"\s+", '_', title))
 
 def quickImshow(a, title=None, path=None, name=None, save_fig=False, save_only=False,
-            clear=True, **kwargs):
+            clear=True, extent=None, show_grid=False, grid_color='w', grid_linewidth=1, **kwargs):
     if title is not None and not isinstance(title, str):
         illegal_value('title', title, 'quickImshow')
         return
@@ -327,24 +327,32 @@
             path = name
         else:
             path = f'{path}/{name}'
+    if extent is None:
+        extent = (0, a.shape[1], a.shape[0], 0)
     if clear:
         plt.close(fig=title)
     if save_only:
         plt.figure(title)
-        plt.imshow(a, **kwargs)
+        plt.imshow(a, extent=extent, **kwargs)
+        if show_grid:
+            ax = plt.gca()
+            ax.grid(color=grid_color, linewidth=grid_linewidth)
         plt.savefig(path)
         plt.close(fig=title)
         #plt.imsave(f'{title}.png', a, **kwargs)
     else:
         plt.ion()
         plt.figure(title)
-        plt.imshow(a, **kwargs)
+        plt.imshow(a, extent=extent, **kwargs)
+        if show_grid:
+            ax = plt.gca()
+            ax.grid(color=grid_color, linewidth=grid_linewidth)
         if save_fig:
             plt.savefig(path)
         plt.pause(1)
 
 def quickPlot(*args, title=None, path=None, name=None, save_fig=False, save_only=False,
-        clear=True, **kwargs):
+        clear=True, show_grid=False, **kwargs):
     if title is not None and not isinstance(title, str):
         illegal_value('title', title, 'quickPlot')
         return
@@ -383,6 +391,9 @@
                plt.plot(*y, **kwargs)
         else:
             plt.plot(*args, **kwargs)
+        if show_grid:
+            ax = plt.gca()
+            ax.grid(color='k')#, linewidth=1)
         plt.savefig(path)
         plt.close(fig=title)
     else:
@@ -393,6 +404,9 @@
                plt.plot(*y, **kwargs)
         else:
             plt.plot(*args, **kwargs)
+        if show_grid:
+            ax = plt.gca()
+            ax.grid(color='k')#, linewidth=1)
         if save_fig:
             plt.savefig(path)
         plt.pause(1)
@@ -536,15 +550,17 @@
             illegal_value('upp', upp, 'selectImageBounds')
             return None
     bounds = [low, upp]
-    a_tmp = a
+    a_tmp = np.copy(a)
+    a_tmp_max = a.max()
     if axis:
-        a_tmp[:,bounds[0]] = a.max()
-        a_tmp[:,bounds[1]] = a.max()
+        a_tmp[:,bounds[0]] = a_tmp_max
+        a_tmp[:,bounds[1]-1] = a_tmp_max
     else:
-        a_tmp[bounds[0],:] = a.max()
-        a_tmp[bounds[1],:] = a.max()
+        a_tmp[bounds[0],:] = a_tmp_max
+        a_tmp[bounds[1]-1,:] = a_tmp_max
     print(f'lower bound = {low} (inclusive)\nupper bound = {upp} (exclusive)')
     quickImshow(a_tmp, title=title)
+    del a_tmp
     if pyip.inputYesNo('Accept these bounds ([y]/n)?: ', blank=True) == 'no':
         bounds = selectImageBounds(a, axis, low=low_save, upp=upp_save, num_min=num_min_save,
             title=title)
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/run_tomo_find_center	Tue Apr 05 18:23:54 2022 +0000
@@ -0,0 +1,3 @@
+#!/bin/bash
+
+python tomo_find_center.py -i red_stacks.npz -c 'config_l_center.yaml' --row_bounds '620 950' --center_rows '670 890' --output_config 'output.yaml' --recon_center_low 'recon_center_low.png' --recon_center_upp 'recon_center_upp.png'
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/run_tomo_reconstruct	Tue Apr 05 18:23:54 2022 +0000
@@ -0,0 +1,3 @@
+#!/bin/bash
+
+python tomo_reconstruct.py -i red_stacks.npz -c 'config_l_recon.yaml' --output_config 'output.yaml' --output_data 'output.npz'
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/run_tomo_setup	Tue Apr 05 18:23:54 2022 +0000
@@ -0,0 +1,3 @@
+#!/bin/bash
+
+python tomo_setup.py -i inputfiles_l.txt -c 'config_l_setup.yaml' --theta_range '0.0 180.0 354' --output_config 'output.yaml' --output_data 'output.npz' --dark 'dark.png' --bright 'bright.png' --tomo 'tomo.png' --detectorbounds 'detectorbounds.png' 1 20 1 354
--- a/tomo.py	Thu Mar 31 20:48:17 2022 +0000
+++ b/tomo.py	Tue Apr 05 18:23:54 2022 +0000
@@ -34,13 +34,16 @@
 
     def __init__(self, num_core):
         cpu_count = mp.cpu_count()
+        logging.debug(f'start: num_core={num_core} cpu_count={cpu_count}')
         if num_core is None or num_core < 1 or num_core > cpu_count:
             self.num_core = cpu_count
         else:
             self.num_core = num_core
+        logging.debug(f'self.num_core={self.num_core}')
 
     def __enter__(self):
         self.num_core_org = ne.set_num_threads(self.num_core)
+        logging.debug(f'self.num_core={self.num_core}')
 
     def __exit__(self, exc_type, exc_value, traceback):
         ne.set_num_threads(self.num_core_org)
@@ -647,12 +650,12 @@
         logging.debug(f'tdf_cutoff = {tdf_cutoff}')
         logging.debug(f'tdf_mean = {tdf_mean}')
         np.nan_to_num(self.tdf, copy=False, nan=tdf_mean, posinf=tdf_mean, neginf=0.)
-        if not self.test_mode and not self.galaxy_flag:
+        if self.galaxy_flag:
+            msnc.quickImshow(self.tdf, title='dark field', name=dark_field_pngname,
+                    save_fig=True, save_only=True)
+        elif not self.test_mode:
             msnc.quickImshow(self.tdf, title='dark field', path=self.output_folder,
                     save_fig=self.save_plots, save_only=self.save_plots_only)
-        elif self.galaxy_flag:
-            msnc.quickImshow(self.tdf, title='dark field', name=dark_field_pngname,
-                    save_fig=True, save_only=True)
 
     def _genBright(self, tbf_files, bright_field_pngname):
         """Generate bright field.
@@ -681,12 +684,12 @@
             self.tbf -= self.tdf
         else:
             logging.warning('Dark field unavailable')
-        if not self.test_mode and not self.galaxy_flag:
+        if self.galaxy_flag:
+            msnc.quickImshow(self.tbf, title='bright field', name=bright_field_pngname,
+                    save_fig=True, save_only=True)
+        elif not self.test_mode:
             msnc.quickImshow(self.tbf, title='bright field', path=self.output_folder,
                     save_fig=self.save_plots, save_only=self.save_plots_only)
-        elif self.galaxy_flag:
-            msnc.quickImshow(self.tbf, title='bright field', name=bright_field_pngname,
-                    save_fig=True, save_only=True)
 
     def _setDetectorBounds(self, tomo_stack_files, tomo_field_pngname, detectorbounds_pngname):
         """Set vertical detector bounds for image stack.
@@ -743,14 +746,16 @@
             tomo_stack = msnc.loadImageStack(tomo_stack_files[0], self.config['data_filetype'],
                 stacks[0]['img_offset'], 1)
             x_sum = np.sum(tomo_stack[0,:,:], 1)
+            x_sum_min = x_sum.min()
+            x_sum_max = x_sum.max()
             title = f'tomography image at theta={self.config["theta_range"]["start"]}'
             msnc.quickImshow(tomo_stack[0,:,:], title=title, name=tomo_field_pngname,
-                    save_fig=True, save_only=True)
+                    save_fig=True, save_only=True, show_grid=True)
             msnc.quickPlot((range(x_sum.size), x_sum),
-                    ([img_x_bounds[0], img_x_bounds[0]], [x_sum.min(), x_sum.max()], 'r-'),
-                    ([img_x_bounds[1]-1, img_x_bounds[1]-1], [x_sum.min(), x_sum.max()], 'r-'),
+                    ([img_x_bounds[0], img_x_bounds[0]], [x_sum_min, x_sum_max], 'r-'),
+                    ([img_x_bounds[1]-1, img_x_bounds[1]-1], [x_sum_min, x_sum_max], 'r-'),
                     title='sum over theta and y', name=detectorbounds_pngname,
-                    save_fig=True, save_only=True)
+                    save_fig=True, save_only=True, show_grid=True)
             
             # Update config and save to file
             if preprocess is None:
@@ -758,6 +763,7 @@
             else:
                 preprocess['img_x_bounds'] = img_x_bounds
             self.cf.saveFile(self.config_out)
+            del x_sum
             return
 
         # For one tomography stack only: load the first image
@@ -777,9 +783,17 @@
             x_sum = np.sum(tomo_stack[0,:,:], 1)
             use_bounds = 'no'
             if img_x_bounds[0] is not None and img_x_bounds[1] is not None:
+                tmp = np.copy(tomo_stack[0,:,:])
+                tmp_max = tmp.max()
+                tmp[img_x_bounds[0],:] = tmp_max
+                tmp[img_x_bounds[1]-1,:] = tmp_max
+                msnc.quickImshow(tmp, title=title)
+                del tmp
+                x_sum_min = x_sum.min()
+                x_sum_max = x_sum.max()
                 msnc.quickPlot((range(x_sum.size), x_sum),
-                        ([img_x_bounds[0], img_x_bounds[0]], [x_sum.min(), x_sum.max()], 'r-'),
-                        ([img_x_bounds[1]-1, img_x_bounds[1]-1], [x_sum.min(), x_sum.max()], 'r-'),
+                        ([img_x_bounds[0], img_x_bounds[0]], [x_sum_min, x_sum_max], 'r-'),
+                        ([img_x_bounds[1]-1, img_x_bounds[1]-1], [x_sum_min, x_sum_max], 'r-'),
                         title='sum over theta and y')
                 print(f'lower bound = {img_x_bounds[0]} (inclusive)\n'+
                         f'upper bound = {img_x_bounds[1]} (exclusive)]')
@@ -799,11 +813,20 @@
                         save_fig=self.save_plots, save_only=True)
         else:
             x_sum = np.sum(self.tbf, 1)
+            x_sum_min = x_sum.min()
+            x_sum_max = x_sum.max()
             use_bounds = 'no'
             if img_x_bounds[0] is not None and img_x_bounds[1] is not None:
+                tmp = np.copy(self.tbf)
+                tmp_max = tmp.max()
+                tmp[img_x_bounds[0],:] = tmp_max
+                tmp[img_x_bounds[1]-1,:] = tmp_max
+                title = 'Bright field'
+                msnc.quickImshow(tmp, title=title)
+                del tmp
                 msnc.quickPlot((range(x_sum.size), x_sum),
-                        ([img_x_bounds[0], img_x_bounds[0]], [x_sum.min(), x_sum.max()], 'r-'),
-                        ([img_x_bounds[1]-1, img_x_bounds[1]-1], [x_sum.min(), x_sum.max()], 'r-'),
+                        ([img_x_bounds[0], img_x_bounds[0]], [x_sum_min, x_sum_max], 'r-'),
+                        ([img_x_bounds[1]-1, img_x_bounds[1]-1], [x_sum_min, x_sum_max], 'r-'),
                         title='sum over theta and y')
                 print(f'lower bound = {img_x_bounds[0]} (inclusive)\n'+
                         f'upper bound = {img_x_bounds[1]} (exclusive)]')
@@ -820,9 +843,16 @@
                     x_upp = int(x_upp+(x_upp-x_low)/10)
                     if x_upp >= x_sum.size:
                         x_upp = x_sum.size
+                    tmp = np.copy(self.tbf)
+                    tmp_max = tmp.max()
+                    tmp[x_low,:] = tmp_max
+                    tmp[x_upp-1,:] = tmp_max
+                    title = 'Bright field'
+                    msnc.quickImshow(tmp, title=title)
+                    del tmp
                     msnc.quickPlot((range(x_sum.size), x_sum),
-                            ([x_low, x_low], [x_sum.min(), x_sum.max()], 'r-'),
-                            ([x_upp, x_upp], [x_sum.min(), x_sum.max()], 'r-'),
+                            ([x_low, x_low], [x_sum_min, x_sum_max], 'r-'),
+                            ([x_upp, x_upp], [x_sum_min, x_sum_max], 'r-'),
                             title='sum over theta and y')
                     print(f'lower bound = {x_low} (inclusive)\nupper bound = {x_upp} (exclusive)]')
                     use_fit =  pyip.inputYesNo('Accept these bounds ([y]/n)?: ', blank=True)
@@ -837,6 +867,7 @@
                         x_sum[img_x_bounds[0]:img_x_bounds[1]],
                         title='sum over theta and y', path=self.output_folder,
                         save_fig=self.save_plots, save_only=True)
+            del x_sum
         logging.debug(f'img_x_bounds: {img_x_bounds}')
 
         if self.save_plots_only:
@@ -1069,13 +1100,10 @@
             stack.pop('reconstructed', 'reconstructed not found')
             find_center = self.config.get('find_center')
             if find_center:
+                find_center.pop('completed', 'completed not found')
                 if self.test_mode:
-                    find_center.pop('completed', 'completed not found')
                     find_center.pop('lower_center_offset', 'lower_center_offset not found')
                     find_center.pop('upper_center_offset', 'upper_center_offset not found')
-                else:
-                    if find_center.get('center_stack_index', -1) == index:
-                        self.config.pop('find_center')
             self.cf.saveFile(self.config_out)
 
         if self.tdf.size:
@@ -1083,7 +1111,7 @@
         del tbf
 
     def _reconstructOnePlane(self, tomo_plane_T, center, thetas_deg, eff_pixel_size,
-            cross_sectional_dim, plot_sinogram=True):
+            cross_sectional_dim, plot_sinogram=True, num_core=1):
         """Invert the sinogram for a single tomography plane.
         """
         # tomo_plane_T index order: column,theta
@@ -1117,24 +1145,28 @@
         recon_sinogram = spi.gaussian_filter(recon_sinogram, 0.5)
         recon_clean = np.expand_dims(recon_sinogram, axis=0)
         del recon_sinogram
-        recon_clean = tomopy.misc.corr.remove_ring(recon_clean, rwidth=17)
+        recon_clean = tomopy.misc.corr.remove_ring(recon_clean, rwidth=17, ncore=num_core)
         logging.debug(f'filtering and removing ring artifact took {time()-t0:.2f} seconds!')
         return recon_clean
 
-    def _plotEdgesOnePlane(self, recon_plane, base_name, weight=0.001):
+    def _plotEdgesOnePlane(self, recon_plane, title, name=None, weight=0.001):
         # RV parameters for the denoise, gaussian, and ring removal will be different for different feature sizes
         edges = denoise_tv_chambolle(recon_plane, weight = weight)
         vmax = np.max(edges[0,:,:])
         vmin = -vmax
-        msnc.quickImshow(edges[0,:,:], f'{base_name} coolwarm', path=self.output_folder,
-                save_fig=self.save_plots, save_only=self.save_plots_only, cmap='coolwarm')
-        msnc.quickImshow(edges[0,:,:], f'{base_name} gray', path=self.output_folder,
-                save_fig=self.save_plots, save_only=self.save_plots_only, cmap='gray',
-                vmin=vmin, vmax=vmax)
+        if self.galaxy_flag:
+            msnc.quickImshow(edges[0,:,:], title, name=name, save_fig=True, save_only=True,
+                    cmap='gray', vmin=vmin, vmax=vmax)
+        else:
+            msnc.quickImshow(edges[0,:,:], f'{title} coolwarm', path=self.output_folder,
+                    save_fig=self.save_plots, save_only=self.save_plots_only, cmap='coolwarm')
+            msnc.quickImshow(edges[0,:,:], f'{title} gray', path=self.output_folder,
+                    save_fig=self.save_plots, save_only=self.save_plots_only, cmap='gray',
+                    vmin=vmin, vmax=vmax)
         del edges
 
     def _findCenterOnePlane(self, sinogram, row, thetas_deg, eff_pixel_size, cross_sectional_dim,
-            tol=0.1):
+            tol=0.1, num_core=1, recon_pngname=None):
         """Find center for a single tomography plane.
         """
         # sinogram index order: theta,column
@@ -1143,20 +1175,31 @@
         center = sinogram.shape[1]/2
 
         # try automatic center finding routines for initial value
-        tomo_center = tomopy.find_center_vo(sinogram)
+        tomo_center = tomopy.find_center_vo(sinogram, ncore=num_core)
         center_offset_vo = tomo_center-center
-        if not self.test_mode:
+        if self.test_mode or self.galaxy_flag:
+            logging.info(f'Center at row {row} using Nghia Vo’s method = {center_offset_vo:.2f}')
+            if self.test_mode:
+                del sinogram_T
+                return float(center_offset_vo)
+        else:
             print(f'Center at row {row} using Nghia Vo’s method = {center_offset_vo:.2f}')
+            if recon_pngname:
+                logging.warning('Ignoring recon_pngname in _findCenterOnePlane (only for Galaxy)')
         recon_plane = self._reconstructOnePlane(sinogram_T, tomo_center, thetas_deg,
-                eff_pixel_size, cross_sectional_dim, False)
-        if not self.test_mode:
-            base_name=f'edges row{row} center_offset_vo{center_offset_vo:.2f}'
-            self._plotEdgesOnePlane(recon_plane, base_name)
-        use_phase_corr = 'no'
-        if not self.test_mode:
-            use_phase_corr = pyip.inputYesNo('Try finding center using phase correlation '+
-                    '(y/[n])? ', blank=True)
-        if use_phase_corr == 'yes':
+                eff_pixel_size, cross_sectional_dim, False, num_core)
+        if self.galaxy_flag:
+            assert(isinstance(recon_pngname, str))
+            title = os.path.basename(recon_pngname)
+            self._plotEdgesOnePlane(recon_plane, title, name=recon_pngname)
+            del sinogram_T
+            del recon_plane
+            return float(center_offset_vo)
+        else:
+            title = f'edges row{row} center_offset_vo{center_offset_vo:.2f}'
+            self._plotEdgesOnePlane(recon_plane, title)
+        if (pyip.inputYesNo('Try finding center using phase correlation '+
+                '(y/[n])? ', blank=True) == 'yes'):
             tomo_center = tomopy.find_center_pc(sinogram, sinogram, tol=0.1,
                     rotc_guess=tomo_center)
             error = 1.
@@ -1168,24 +1211,18 @@
             center_offset = tomo_center-center
             print(f'Center at row {row} using phase correlation = {center_offset:.2f}')
             recon_plane = self._reconstructOnePlane(sinogram_T, tomo_center, thetas_deg,
-                    eff_pixel_size, cross_sectional_dim, False)
-            base_name=f'edges row{row} center_offset{center_offset:.2f}'
-            self._plotEdgesOnePlane(recon_plane, base_name)
-        accept_center = 'yes'
-        if not self.test_mode:
-            accept_center = pyip.inputYesNo('Accept a center location ([y]) or continue '+
-                    'search (n)? ', blank=True)
-        if accept_center != 'no':
+                    eff_pixel_size, cross_sectional_dim, False, num_core)
+            title = f'edges row{row} center_offset{center_offset:.2f}'
+            self._plotEdgesOnePlane(recon_plane, title)
+        if (pyip.inputYesNo('Accept a center location ([y]) or continue '+
+                'search (n)? ', blank=True) != 'no'):
             del sinogram_T
             del recon_plane
-            if self.test_mode:
+            center_offset = pyip.inputNum(
+                    f'    Enter chosen center offset [{-int(center)}, {int(center)}] '+
+                    f'([{center_offset_vo}])): ', blank=True)
+            if center_offset == '':
                 center_offset = center_offset_vo
-            else:
-                center_offset = pyip.inputNum(
-                        f'    Enter chosen center offset [{-int(center)}, {int(center)}] '+
-                        f'([{center_offset_vo}])): ', blank=True)
-                if center_offset == '':
-                    center_offset = center_offset_vo
             return float(center_offset)
 
         while True:
@@ -1204,9 +1241,9 @@
                         center_offset_step):
                 logging.info(f'center_offset = {center_offset}')
                 recon_plane = self._reconstructOnePlane(sinogram_T, center_offset+center,
-                        thetas_deg, eff_pixel_size, cross_sectional_dim, False)
-                base_name=f'edges row{row} center_offset{center_offset}'
-                self._plotEdgesOnePlane(recon_plane, base_name)
+                        thetas_deg, eff_pixel_size, cross_sectional_dim, False, num_core)
+                title = f'edges row{row} center_offset{center_offset}'
+                self._plotEdgesOnePlane(recon_plane, title)
             if pyip.inputInt('\nContinue (0) or end the search (1): ', min=0, max=1):
                 break
 
@@ -1266,7 +1303,8 @@
                     init_recon=tomo_recon_stack, options=options, sinogram_order=True,
                     algorithm=tomopy.astra, ncore=num_core)
         if True:
-            tomopy.misc.corr.remove_ring(tomo_recon_stack, rwidth=rwidth, out=tomo_recon_stack)
+            tomopy.misc.corr.remove_ring(tomo_recon_stack, rwidth=rwidth, out=tomo_recon_stack,
+                    ncore=num_core)
         return tomo_recon_stack
 
     def findImageFiles(self):
@@ -1359,6 +1397,22 @@
 
         return dark_files, bright_files, tomo_stack_files
 
+    def loadTomoStacks(self, input_name):
+        """Load tomography stacks (only for Galaxy).
+        """
+        assert(self.galaxy_flag)
+        t0 = time()
+        logging.info(f'Loading preprocessed tomography stack from {input_name} ...')
+        stack_info = self.config['stack_info']
+        stacks = stack_info.get('stacks')
+        assert(len(self.tomo_stacks) == stack_info['num'])
+        with np.load(input_name) as f:
+            for i,stack in enumerate(stacks):
+                self.tomo_stacks[i] = f[f'set_{stack["index"]}']
+                logging.info(f'loaded stack {i}: index = {stack["index"]}, shape = '+
+                        f'{self.tomo_stacks[i].shape}')
+        logging.info(f'... done in {time()-t0:.2f} seconds!')
+
     def genTomoStacks(self, tdf_files=None, tbf_files=None, tomo_stack_files=None,
             dark_field_pngname=None, bright_field_pngname=None, tomo_field_pngname=None,
             detectorbounds_pngname=None, output_name=None):
@@ -1459,13 +1513,34 @@
             stack['ref_height'] = stack['ref_height']+pixel_size
         self.cf.saveFile(self.config_out)
 
-    def findCenters(self):
+    def findCenters(self, row_bounds=None, center_rows=None, recon_low_pngname=None,
+            recon_upp_pngname=None, num_core=None):
         """Find rotation axis centers for the tomography stacks.
         """
+        if num_core is None:
+            num_core = self.num_core
         logging.debug('Find centers for tomography stacks')
         stacks = self.config['stack_info']['stacks']
         available_stacks = [stack['index'] for stack in stacks if stack.get('preprocessed', False)]
         logging.debug('Available stacks: {available_stacks}')
+        if self.galaxy_flag:
+            assert(isinstance(row_bounds, list) and len(row_bounds) == 2)
+            assert(isinstance(center_rows, list) and len(center_rows) == 2)
+            assert(isinstance(recon_low_pngname, str))
+            assert(isinstance(recon_upp_pngname, str))
+        else:
+            if row_bounds:
+                logging.warning('Ignoring row_bounds in findCenters (only for Galaxy)')
+                row_bounds = None
+            if center_rows:
+                logging.warning('Ignoring center_rows in findCenters (only for Galaxy)')
+                center_rows = None
+            if recon_low_pngname:
+                logging.warning('Ignoring recon_low_pngname in findCenters (only for Galaxy)')
+                recon_low_pngname = None
+            if recon_upp_pngname:
+                logging.warning('Ignoring recon_upp_pngname in findCenters (only for Galaxy)')
+                recon_upp_pngname = None
 
         # Check for valid available center stack index
         find_center = self.config.get('find_center')
@@ -1481,8 +1556,8 @@
                 else:
                     print('\nFound calibration center offset info for stack '+
                             f'{center_stack_index}')
-                    if (pyip.inputYesNo('Do you want to use this again (y/n)? ') == 'yes' and
-                            find_center.get('completed', False) == True):
+                    if (pyip.inputYesNo('Do you want to use this again ([y]/n)? ',
+                            blank=True) != 'no' and find_center.get('completed', False) == True):
                         return
 
         # Load the required preprocessed stack
@@ -1503,6 +1578,9 @@
                 stacks[0]['preprocessed'] = False
                 raise OSError('Unable to load the required preprocessed tomography stack')
             assert(stacks[0].get('preprocessed', False) == True)
+        elif self.galaxy_flag:
+            logging.error('CHECK/FIX FOR GALAXY')
+            #center_stack_index = stacks[int(num_tomo_stacks/2)]['index']
         else:
             while True:
                 if not center_stack_index:
@@ -1528,6 +1606,10 @@
             find_center = self.config['find_center']
         else:
             find_center['center_stack_index'] = center_stack_index
+        if not self.galaxy_flag:
+            row_bounds = find_center.get('row_bounds', None)
+            center_rows = [find_center.get('lower_row', None),
+                    find_center.get('upper_row', None)]
 
         # Set thetas (in degrees)
         theta_range = self.config['theta_range']
@@ -1545,26 +1627,54 @@
             raise ValueError('Detector pixel size unavailable')
         eff_pixel_size = 100.*pixel_size/zoom_perc
         logging.debug(f'eff_pixel_size = {eff_pixel_size}')
-        tomo_ref_heights = [stack['ref_height'] for stack in stacks]
         if num_tomo_stacks == 1:
-            n1 = 0
-            height = center_stack.shape[0]*eff_pixel_size
-            if not self.test_mode and pyip.inputYesNo('\nDo you want to reconstruct the full '+
-                    f'sample height ({height:.3f} mm) (y/n)? ') == 'no':
-                height = pyip.inputNum('\nEnter the desired reconstructed sample height '+
-                        f'in mm [0, {height:.3f}]: ', min=0, max=height)
-            n1 = int(0.5*(center_stack.shape[0]-height/eff_pixel_size))
+             accept = 'yes'
+             if not self.test_mode and not self.galaxy_flag:
+                 accept = 'no'
+                 print('\nSelect bounds for image reconstruction')
+                 if msnc.is_index_range(row_bounds, 0, center_stack.shape[0]):
+                     a_tmp = np.copy(center_stack[:,0,:])
+                     a_tmp_max = a_tmp.max()
+                     a_tmp[row_bounds[0],:] = a_tmp_max
+                     a_tmp[row_bounds[1]-1,:] = a_tmp_max
+                     print(f'lower bound = {row_bounds[0]} (inclusive)')
+                     print(f'upper bound = {row_bounds[1]} (exclusive)')
+                     msnc.quickImshow(a_tmp, title=f'center stack theta={theta_start}',
+                         aspect='auto')
+                     del a_tmp
+                     accept = pyip.inputYesNo('Accept these bounds ([y]/n)?: ', blank=True)
+             if accept == 'no':
+                 (n1, n2) = msnc.selectImageBounds(center_stack[:,0,:], 0,
+                         title=f'center stack theta={theta_start}')
+             else:
+                 n1 = row_bounds[0]
+                 n2 = row_bounds[1]
         else:
+            logging.error('CHECK/FIX FOR GALAXY')
+            tomo_ref_heights = [stack['ref_height'] for stack in stacks]
             n1 = int((1.+(tomo_ref_heights[0]+center_stack.shape[0]*eff_pixel_size-
                 tomo_ref_heights[1])/eff_pixel_size)/2)
-        n2 = center_stack.shape[0]-n1
+            n2 = center_stack.shape[0]-n1
         logging.info(f'n1 = {n1}, n2 = {n2} (n2-n1) = {(n2-n1)*eff_pixel_size:.3f} mm')
         if not center_stack.size:
             RuntimeError('Center stack not loaded')
-        if not self.test_mode:
-            msnc.quickImshow(center_stack[:,0,:], title=f'center stack theta={theta_start}',
-                    path=self.output_folder, save_fig=self.save_plots,
-                    save_only=self.save_plots_only)
+        if not self.test_mode and not self.galaxy_flag:
+            tmp = center_stack[:,0,:]
+            tmp_max = tmp.max()
+            tmp[n1,:] = tmp_max
+            tmp[n2-1,:] = tmp_max
+            if msnc.is_index_range(center_rows, 0, tmp.shape[0]):
+                tmp[center_rows[0],:] = tmp_max
+                tmp[center_rows[1]-1,:] = tmp_max
+            extent = [0, tmp.shape[1], tmp.shape[0], 0]
+            msnc.quickImshow(tmp, title=f'center stack theta={theta_start}',
+                    path=self.output_folder, extent=extent, save_fig=self.save_plots,
+                    save_only=self.save_plots_only, aspect='auto')
+            del tmp
+            #extent = [0, center_stack.shape[2], n2, n1]
+            #msnc.quickImshow(center_stack[n1:n2,0,:], title=f'center stack theta={theta_start}',
+            #        path=self.output_folder, extent=extent, save_fig=self.save_plots,
+            #        save_only=self.save_plots_only, show_grid=True, aspect='auto')
 
         # Get cross sectional diameter in mm
         cross_sectional_dim = center_stack.shape[2]*eff_pixel_size
@@ -1576,9 +1686,9 @@
         # Lower row center
         use_row = 'no'
         use_center = 'no'
-        row = find_center.get('lower_row')
+        row = center_rows[0]
         if msnc.is_int(row, n1, n2-2):
-            if self.test_mode:
+            if self.test_mode or self.galaxy_flag:
                 assert(row is not None)
                 use_row = 'yes'
             else:
@@ -1605,8 +1715,8 @@
                     msnc.clearFig(f'theta={theta_start}')
             # center_stack order: row,theta,column
             center_offset = self._findCenterOnePlane(center_stack[row,:,:], row, thetas_deg,
-                    eff_pixel_size, cross_sectional_dim)
-        logging.info(f'Lower center offset = {center_offset}')
+                    eff_pixel_size, cross_sectional_dim, num_core=num_core,
+                    recon_pngname=recon_low_pngname)
 
         # Update config and save to file
         find_center['row_bounds'] = [n1, n2]
@@ -1618,9 +1728,9 @@
         # Upper row center
         use_row = 'no'
         use_center = 'no'
-        row = find_center.get('upper_row')
+        row = center_rows[1]
         if msnc.is_int(row, lower_row+1, n2-1):
-            if self.test_mode:
+            if self.test_mode or self.galaxy_flag:
                 assert(row is not None)
                 use_row = 'yes'
             else:
@@ -1647,7 +1757,8 @@
                     msnc.clearFig(f'theta={theta_start}')
             # center_stack order: row,theta,column
             center_offset = self._findCenterOnePlane(center_stack[row,:,:], row, thetas_deg,
-                    eff_pixel_size, cross_sectional_dim)
+                    eff_pixel_size, cross_sectional_dim, num_core=num_core,
+                    recon_pngname=recon_upp_pngname)
         logging.info(f'upper_center_offset = {center_offset}')
         del center_stack
 
@@ -1731,14 +1842,22 @@
         # Update config file
         self.config = msnc.update('config.txt', 'check_centers', True, 'find_centers')
 
-    def reconstructTomoStacks(self):
+    def reconstructTomoStacks(self, output_name=None, num_core=None):
         """Reconstruct tomography stacks.
         """
+        if num_core is None:
+            num_core = self.num_core
         logging.debug('Reconstruct tomography stacks')
         stacks = self.config['stack_info']['stacks']
         assert(len(self.tomo_stacks) == self.config['stack_info']['num'])
         assert(len(self.tomo_stacks) == len(stacks))
         assert(len(self.tomo_recon_stacks) == len(stacks))
+        if self.galaxy_flag:
+            assert(isinstance(output_name, str))
+        else:
+            if output_name:
+                logging.warning('Ignoring output_name in reconstructTomoStacks '+
+                    '(only used in Galaxy)')
 
         # Get rotation axis rows and centers
         find_center = self.config['find_center']
@@ -1784,56 +1903,65 @@
             # reconstructed stack order for each one in stack : row/z,x,y
             # preprocessed stack order for each one in stack: row,theta,column
             index = stack['index']
-            available = False
-            if stack.get('reconstructed', False):
-                self.tomo_recon_stacks[i], available = self._loadTomo('recon stack', index)
-            if self.tomo_recon_stacks[i].size or available:
-                if self.tomo_stacks[i].size:
-                    self.tomo_stacks[i] = np.array([])
-                assert(stack.get('preprocessed', False) == True)
-                assert(stack.get('reconstructed', False) == True)
-                continue
-            else:
-                stack['reconstructed'] = False
-                if not self.tomo_stacks[i].size:
-                    self.tomo_stacks[i], available = self._loadTomo('red stack', index,
-                            required=True)
-                if not self.tomo_stacks[i].size:
-                    logging.error(f'Unable to load tomography stack {index} for reconstruction')
-                    stack[i]['preprocessed'] = False
-                    load_error = True
+            if not self.galaxy_flag:
+                available = False
+                if stack.get('reconstructed', False):
+                    self.tomo_recon_stacks[i], available = self._loadTomo('recon stack', index)
+                if self.tomo_recon_stacks[i].size or available:
+                    if self.tomo_stacks[i].size:
+                        self.tomo_stacks[i] = np.array([])
+                    assert(stack.get('preprocessed', False) == True)
+                    assert(stack.get('reconstructed', False) == True)
                     continue
-                assert(0 <= lower_row < upper_row < self.tomo_stacks[i].shape[0])
-                center_offsets = [lower_center_offset-lower_row*center_slope,
-                        upper_center_offset+(self.tomo_stacks[i].shape[0]-1-upper_row)*center_slope]
-                t0 = time()
-                self.tomo_recon_stacks[i]= self._reconstructOneTomoStack(self.tomo_stacks[i],
-                        thetas, center_offsets=center_offsets, sigma=0.1, num_core=self.num_core,
-                        algorithm='gridrec', run_secondary_sirt=True, secondary_iter=25)
-                logging.info(f'Reconstruction of stack {index} took {time()-t0:.2f} seconds!')
-                if not self.test_mode:
-                    row_slice = int(self.tomo_stacks[i].shape[0]/2) 
-                    title = f'{basetitle} {index} slice{row_slice}'
-                    msnc.quickImshow(self.tomo_recon_stacks[i][row_slice,:,:], title=title,
-                            path=self.output_folder, save_fig=self.save_plots,
-                            save_only=self.save_plots_only)
-                    msnc.quickPlot(self.tomo_recon_stacks[i]
-                            [row_slice,int(self.tomo_recon_stacks[i].shape[2]/2),:],
-                            title=f'{title} cut{int(self.tomo_recon_stacks[i].shape[2]/2)}',
-                            path=self.output_folder, save_fig=self.save_plots,
-                            save_only=self.save_plots_only)
-                    self._saveTomo('recon stack', self.tomo_recon_stacks[i], index)
-#                else:
-#                    np.savetxt(self.output_folder+f'recon_stack_{index}.txt',
-#                            self.tomo_recon_stacks[i][row_slice,:,:], fmt='%.6e')
-                self.tomo_stacks[i] = np.array([])
+            stack['reconstructed'] = False
+            if not self.tomo_stacks[i].size:
+                self.tomo_stacks[i], available = self._loadTomo('red stack', index,
+                        required=True)
+            if not self.tomo_stacks[i].size:
+                logging.error(f'Unable to load tomography stack {index} for reconstruction')
+                stack[i]['preprocessed'] = False
+                load_error = True
+                continue
+            assert(0 <= lower_row < upper_row < self.tomo_stacks[i].shape[0])
+            center_offsets = [lower_center_offset-lower_row*center_slope,
+                    upper_center_offset+(self.tomo_stacks[i].shape[0]-1-upper_row)*center_slope]
+            t0 = time()
+            self.tomo_recon_stacks[i]= self._reconstructOneTomoStack(self.tomo_stacks[i],
+                    thetas, center_offsets=center_offsets, sigma=0.1, num_core=num_core,
+                    algorithm='gridrec', run_secondary_sirt=True, secondary_iter=25)
+            logging.info(f'Reconstruction of stack {index} took {time()-t0:.2f} seconds!')
+            if not self.test_mode and not self.galaxy_flag:
+                row_slice = int(self.tomo_stacks[i].shape[0]/2) 
+                title = f'{basetitle} {index} slice{row_slice}'
+                msnc.quickImshow(self.tomo_recon_stacks[i][row_slice,:,:], title=title,
+                        path=self.output_folder, save_fig=self.save_plots,
+                        save_only=self.save_plots_only)
+                msnc.quickPlot(self.tomo_recon_stacks[i]
+                        [row_slice,int(self.tomo_recon_stacks[i].shape[2]/2),:],
+                        title=f'{title} cut{int(self.tomo_recon_stacks[i].shape[2]/2)}',
+                        path=self.output_folder, save_fig=self.save_plots,
+                        save_only=self.save_plots_only)
+                self._saveTomo('recon stack', self.tomo_recon_stacks[i], index)
+#            else:
+#                np.savetxt(self.output_folder+f'recon_stack_{index}.txt',
+#                        self.tomo_recon_stacks[i][row_slice,:,:], fmt='%.6e')
+            self.tomo_stacks[i] = np.array([])
 
-                # Update config and save to file
-                stack['reconstructed'] = True
-                combine_stacks = self.config.get('combine_stacks')
-                if combine_stacks and index in combine_stacks.get('stacks', []):
-                    combine_stacks['stacks'].remove(index)
-                self.cf.saveFile(self.config_out)
+            # Update config and save to file
+            stack['reconstructed'] = True
+            combine_stacks = self.config.get('combine_stacks')
+            if combine_stacks and index in combine_stacks.get('stacks', []):
+                combine_stacks['stacks'].remove(index)
+            self.cf.saveFile(self.config_out)
+
+        # Save reconstructed tomography stack to file
+        if self.galaxy_flag:
+            t0 = time()
+            logging.info(f'Saving reconstructed tomography stack to {output_name} ...')
+            save_stacks = {f'set_{stack["index"]}':tomo_stack
+                    for stack,tomo_stack in zip(stacks,self.tomo_recon_stacks)}
+            np.savez(output_name, **save_stacks)
+            logging.info(f'... done in {time()-t0:.2f} seconds!')
 
     def combineTomoStacks(self):
         """Combine the reconstructed tomography stacks.
@@ -2095,6 +2223,7 @@
             default=False,
             help='Test mode flag')
     parser.add_argument('--num_core',
+            type=int,
             default=-1,
             help='Number of cores')
     args = parser.parse_args()
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/tomo_find_center.py	Tue Apr 05 18:23:54 2022 +0000
@@ -0,0 +1,74 @@
+#!/usr/bin/env python3
+
+import logging
+
+import sys
+import re
+import argparse
+
+from tomo import Tomo
+
+def __main__():
+
+    # Parse command line arguments
+    parser = argparse.ArgumentParser(
+            description='Find the center axis for a tomography reconstruction')
+    parser.add_argument('-i', '--input_stacks',
+            help='Preprocessed image file stacks')
+    parser.add_argument('-c', '--config',
+            help='Input config')
+    parser.add_argument('--row_bounds',
+            help='Reconstruction row bounds')
+    parser.add_argument('--center_rows',
+            help='Center finding rows')
+    parser.add_argument('--output_config',
+            help='Output config')
+    parser.add_argument('--recon_center_low',
+            help='Lower reconstructed slice center')
+    parser.add_argument('--recon_center_upp',
+            help='Upper reconstructed slice center')
+    parser.add_argument('-l', '--log', 
+            type=argparse.FileType('w'),
+            default=sys.stdout,
+            help='Log file')
+    args = parser.parse_args()
+
+    indexRegex = re.compile(r'\d+')
+    row_bounds = [int(i) for i in indexRegex.findall(args.row_bounds)]
+    center_rows = [int(i) for i in indexRegex.findall(args.center_rows)]
+
+    # Set basic log configuration
+    logging_format = '%(asctime)s : %(levelname)s - %(module)s : %(funcName)s - %(message)s'
+    log_level = 'INFO'
+    level = getattr(logging, log_level.upper(), None)
+    if not isinstance(level, int):
+        raise ValueError(f'Invalid log_level: {log_level}')
+    logging.basicConfig(format=logging_format, level=level, force=True,
+            handlers=[logging.StreamHandler()])
+
+    logging.debug(f'input_stacks = {args.input_stacks}')
+    logging.debug(f'config = {args.config}')
+    logging.debug(f'row_bounds = {args.row_bounds} {row_bounds}')
+    logging.debug(f'center_rows = {args.center_rows} {center_rows}')
+    logging.debug(f'output_config = {args.output_config}')
+    logging.debug(f'recon_center_low = {args.recon_center_low}')
+    logging.debug(f'recon_center_uppm = {args.recon_center_upp}')
+    logging.debug(f'log = {args.log}')
+    logging.debug(f'is log stdout? {args.log is sys.stdout}')
+
+    # Instantiate Tomo object
+    tomo = Tomo(config_file=args.config, config_out=args.output_config, log_level=log_level,
+            log_stream=args.log, galaxy_flag=True)
+    if not tomo.is_valid:
+        raise ValueError('Invalid config file provided.')
+    logging.debug(f'config:\n{tomo.config}')
+
+    # Load preprocessed image files
+    tomo.loadTomoStacks(args.input_stacks)
+
+    # Find centers
+    tomo.findCenters(row_bounds, center_rows, args.recon_center_low, args.recon_center_upp)
+
+if __name__ == "__main__":
+    __main__()
+
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/tomo_find_center.xml	Tue Apr 05 18:23:54 2022 +0000
@@ -0,0 +1,48 @@
+<tool id="tomo_find_center" name="Tomo Find Center Axis" version="0.1.0" python_template_version="3.9">
+    <description>Find the center axis for a tomography reconstruction</description>
+    <requirements>
+        <requirement type="package" version="1.11.0">tomopy</requirement>
+        <requirement type="package" version="3.6.0">h5py</requirement>
+    </requirements>
+    <command detect_errors="exit_code"><![CDATA[
+        $__tool_directory__/tomo_find_center.py
+        -i '$red_stacks'
+        -c '$config'
+        --row_bounds '$row_bound_low $row_bound_upp'
+        --center_rows '$lower_row $upper_row'
+        --output_config 'output_config.yaml'
+        --recon_center_low 'recon_center_low.png'
+        --recon_center_upp 'recon_center_upp.png'
+        -l '$log'
+    ]]></command>
+    <inputs>
+        <param name="red_stacks" type='data' format='npz' optional='false' label="Reduced stacks"/>
+        <param name="config" type='data' format='yaml' optional='false' label="Input config"/>
+        <section name="row_bounds" title="Reconstruction row bounds">
+            <param name="row_bound_low" type="integer" label="Lower bound"/>
+            <param name="row_bound_upp" type="integer" label="Upper bound"/>
+        </section>
+        <section name="recon_centers" title="Reconstruction rows to establish center axis">
+            <param name="lower_row" type="integer" label="Lower row"/>
+            <param name="upper_row" type="integer" label="Upper row"/>
+        </section>
+    </inputs>
+    <outputs>
+        <data name="output_config" format="yaml" label="Output config" from_work_dir="output_config.yaml"/>
+        <data name="recon_center_low" format="png" label="Recontructed slice lower center" from_work_dir="recon_center_low.png"/>
+        <data name="recon_center_upp" format="png" label="Recontructed slice upper center" from_work_dir="recon_center_upp.png"/>
+        <data name="log" format="txt" label="Log"/>
+    </outputs>
+    <help><![CDATA[
+        Preprocess tomography images.
+    ]]></help>
+    <citations>
+        <citation type="bibtex">
+@misc{githubsum_files,
+  author = {Verberg, Rolf},
+  year = {2022},
+  title = {Tomo Find Center Axis},
+}</citation>
+    </citations>
+
+</tool>
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/tomo_reconstruct.py	Tue Apr 05 18:23:54 2022 +0000
@@ -0,0 +1,60 @@
+#!/usr/bin/env python3
+
+import logging
+
+import sys
+import argparse
+
+from tomo import Tomo
+
+def __main__():
+
+    # Parse command line arguments
+    parser = argparse.ArgumentParser(
+            description='Perfrom a tomography reconstruction')
+    parser.add_argument('-i', '--input_stacks',
+            help='Preprocessed image file stacks')
+    parser.add_argument('-c', '--config',
+            help='Input config')
+    parser.add_argument('--output_config',
+            help='Output config')
+    parser.add_argument('--output_data',
+            help='Reconstructed tomography data')
+    parser.add_argument('-l', '--log', 
+            type=argparse.FileType('w'),
+            default=sys.stdout,
+            help='Log file')
+    args = parser.parse_args()
+
+    # Set basic log configuration
+    logging_format = '%(asctime)s : %(levelname)s - %(module)s : %(funcName)s - %(message)s'
+    log_level = 'INFO'
+    level = getattr(logging, log_level.upper(), None)
+    if not isinstance(level, int):
+        raise ValueError(f'Invalid log_level: {log_level}')
+    logging.basicConfig(format=logging_format, level=level, force=True,
+            handlers=[logging.StreamHandler()])
+
+    logging.debug(f'input_stacks = {args.input_stacks}')
+    logging.debug(f'config = {args.config}')
+    logging.debug(f'output_config = {args.output_config}')
+    logging.debug(f'output_data = {args.output_data}')
+    logging.debug(f'log = {args.log}')
+    logging.debug(f'is log stdout? {args.log is sys.stdout}')
+
+    # Instantiate Tomo object
+    tomo = Tomo(config_file=args.config, config_out=args.output_config, log_level=log_level,
+            log_stream=args.log, galaxy_flag=True)
+    if not tomo.is_valid:
+        raise ValueError('Invalid config file provided.')
+    logging.debug(f'config:\n{tomo.config}')
+
+    # Load preprocessed image files
+    tomo.loadTomoStacks(args.input_stacks)
+
+    # Reconstruct tomography stacks
+    tomo.reconstructTomoStacks(args.output_data)
+
+if __name__ == "__main__":
+    __main__()
+
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/tomo_reconstruct.xml	Tue Apr 05 18:23:54 2022 +0000
@@ -0,0 +1,36 @@
+<tool id="tomo_reconstruct" name="Tomo Reconstruction" version="0.1.0" python_template_version="3.9">
+    <description>Perform a tomography reconstruction</description>
+    <requirements>
+        <requirement type="package" version="1.11.0">tomopy</requirement>
+        <requirement type="package" version="3.6.0">h5py</requirement>
+    </requirements>
+    <command detect_errors="exit_code"><![CDATA[
+        $__tool_directory__/tomo_reconstruct.py
+        -i '$red_stacks'
+        -c '$config'
+        --output_data 'output.npz'
+        --output_config 'output_config.yaml'
+        -l '$log'
+    ]]></command>
+    <inputs>
+        <param name="red_stacks" type='data' format='npz' optional='false' label="Reduced stacks"/>
+        <param name="config" type='data' format='yaml' optional='false' label="Input config"/>
+    </inputs>
+    <outputs>
+        <data name="output_data" format="npz" label="Reconstructed tomography stacks" from_work_dir="output_data.npz"/>
+        <data name="output_config" format="yaml" label="Output config" from_work_dir="output_config.yaml"/>
+        <data name="log" format="txt" label="Log"/>
+    </outputs>
+    <help><![CDATA[
+        Reconstruct tomography images.
+    ]]></help>
+    <citations>
+        <citation type="bibtex">
+@misc{githubsum_files,
+  author = {Verberg, Rolf},
+  year = {2022},
+  title = {Tomo Reconstruction},
+}</citation>
+    </citations>
+
+</tool>
--- a/tomo_setup.xml	Thu Mar 31 20:48:17 2022 +0000
+++ b/tomo_setup.xml	Tue Apr 05 18:23:54 2022 +0000
@@ -60,9 +60,6 @@
   author = {Verberg, Rolf},
   year = {2022},
   title = {Tomo Setup},
-  publisher = {GitHub},
-  journal = {GitHub repository},
-  url = {somewhere},
 }</citation>
     </citations>