Snakemake workflow: Reconstructing raw tomography data using tomopy.

public public 1yr ago Version: Version 1 0 bookmarks

Snakemake workflow: Reconstructing raw tomography data

A Snakemake worfklow for tomographically reconstructing raw data using tomopy .

Installation

First download this repo and navigate to it

git clone https://codebase.helmholtz.cloud/gernha62/reconstructing-raw-tomography-data.git cd /path/to/repo

(Optional) Download the example folder with:

wget -m -np https://doi2.psi.ch/datasets/das/work/p15/p15869/compression/MI04_02/tif

Create a virtual environment and install all necessary packages (requires conda):

conda env create --name reconstr_env --file workflow/envs/reconstr.yml

Activate the new virtual environment:

conda activate reconstr_env

Configuration

To configure the workflow, adapt the config file found at config/config.yaml . The config looks as follows:

number_of_darks: 50 number_of_flats: 100 number_of_projections: 501 rotation_center: 508.77 raw_data: MI04_02: doi2.psi.ch/datasets/das/work/p15/p15869/compression/MI04_02/tif

In the config, adjust number_of_darks , number_of_flats , number_of_projections and rotation_center to the number of darks, flats, projections and the rotation center of your dataset. The necessary information can usually be found in the .log file of the folder that contains the raw data.

MI04_02: doi2.psi.ch/datasets/das/work/p15/p15869/compression/MI04_02/tif denotes the path to the example folder used for reconstruction and the keyword MI04_02 will be used to name the output (e.g. in this case the output folder will be named recon_dir_MI04_02 ). Replace the examle path with the path to the dataset you want to reconstruct. Additionally, if you want the name of the output folder to have a different suffix, replace the keyword MI04_02 with a name you prefer.

Run the workflow

If the .tif files contain a numerical prefix that is not separated from the actual image index, it is best to first rename the files. The files will be renamed to 00001.tif , 00002.tif and so on. If the renaming is needed, run:

snakemake --cores 1 'logs/renamefile_MI04_02.log'

If you replaced the keyword MI04_02 in the config file then adjust the command accordingly (e.g. if you replaced the keyword with Tomo_dataset then the command should be snakemake --cores 1 'logs/renamefile_Tomo_dataset.log' ).

Before trying to compute the reconstructions, make sure you have enough memory available (ideally more than 60 GB). To compute the reconstructions using one core, use the command:

snakemake --cores 1

If you want to use all available cores instead, use:

snakemake --cores all

This creates a folder in results with the reconstructed data.

Code Snippets

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
import tomopy
import dxchange
import numpy as np
from glob import glob

def reconstruction(path_to_output, tifFile, numDarks, numFlats, numProjs, rot_center_init):
    #Load Raw Data into Memory
    tifRef = tifFile
    firstIndex = 1
    #Read Darks
    darkStart = firstIndex
    darkStop = darkStart + numDarks 
    dark = dxchange.reader.read_tiff_stack(tifRef, [i for i in range(darkStart, darkStop)])
    #Read Flats
    flatStart = darkStop
    flatStop = flatStart + numFlats
    flat = dxchange.reader.read_tiff_stack(tifRef, [i for i in range(flatStart, flatStop)])
    #Read Projections
    projStart = flatStop
    projStop = projStart + numProjs
    proj = dxchange.reader.read_tiff_stack(tifRef, [i for i in range(projStart, projStop)])
    print ("Shape of projection cube", proj.shape)
    #Define equiangular angles (angles here not part of raw data)
    theta = tomopy.angles(proj.shape[0])
    #Start Reconstruction
    #Correct for dark and flat
    threads = 12
    proj = tomopy.normalize(proj, flat, dark, ncore=threads)
    proj = tomopy.minus_log(proj, ncore=threads)
    rot_center = rot_center_init
    #Inverse Radon
    recon = tomopy.recon(proj, theta, center=rot_center, algorithm='gridrec', sinogram_order=False, ncore=threads)
    #Cosmetics
    recon = tomopy.circ_mask(recon, axis=0, ratio=0.95, ncore=threads)
    print("Reconstruction finished")
    dxchange.write_tiff_stack(recon, fname=path_to_output + '/recon')
    print ("Finished")

def main():
    path_to_directory = snakemake.input[0]
    fnames = glob(path_to_directory + '/' + '*.tif')
    tiff_file_example = fnames[0]
    path_to_output = snakemake.output[0]
    numDarks_magma = snakemake.params.numDarks
    numFlats_magma = snakemake.params.numFlats
    numProjs_magma = snakemake.params.numProjs
    rot_center_magma = snakemake.params.rotCenter
    reconstruction(path_to_output, tiff_file_example, numDarks_magma, numFlats_magma, numProjs_magma, rot_center_magma)

if __name__ == '__main__':
    main()
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
import os
from glob import glob
import sys

#numerical prefix that is not separated from the actual image index confuses the dxchange filename parser, 
# so it is better to first rename the files. This function renames the files into 00000.tif, 00001.tif, ...
def rename_tifs(path_to_directory):
    fnames = glob(path_to_directory + '/' + '*.tif')
    fnames.sort()
    for f in range(len(fnames)):
        new_name = path_to_directory + '/' + '{:0>5}'.format(str(f+1)) + '.tif'
        os.rename(fnames[f], new_name)
    with open(snakemake.log[0], "w") as f:
        sys.stderr = sys.stdout = f
        f.write('success')

def main():
    path_to_directory = snakemake.input[0]
    rename_tifs(path_to_directory)

if __name__ == '__main__':
    main()
20
21
script:
    'scripts/rename_tifs.py'
33
34
script:
    'scripts/reconstructs_tomo_dataset.py'
ShowHide 3 more snippets with no or duplicated tags.

Login to post a comment if you would like to share your experience with this workflow.

Do you know this workflow well? If so, you can request seller status , and start supporting this workflow.

Free

Created: 1yr ago
Updated: 1yr ago
Maitainers: public
URL: https://codebase.helmholtz.cloud/gernha62/reconstructing-raw-tomography-data.git
Name: reconstructing-raw-tomography-data
Version: Version 1
Badge:
workflow icon

Insert copied code into your website to add a link to this workflow.

Downloaded: 0
Copyright: Public Domain
License: MIT License
  • Future updates

Related Workflows

cellranger-snakemake-gke
snakemake workflow to run cellranger on a given bucket using gke.
A Snakemake workflow for running cellranger on a given bucket using Google Kubernetes Engine. The usage of this workflow ...