Add minimal demo-scripts

examples/common_data_processing_workflow/data_reduction_cli.py

@@ -108,6 +108,7 @@
         output_format = "folder"
         if os.path.isdir(output_path):
             num = 1
+            output_path = os.path.normpath(output_path)
             parent_path = os.path.dirname(output_path)
             last_folder = os.path.basename(output_path)
             while True:
@@ -182,7 +183,7 @@
             dsp_folder = folder_path + "_dsp_tmp/"
             dsp_path = dsp_folder + file_name
         else:
-            dsp_folder = dsp_path = output_path + "_dsp_tmp/"
+            dsp_folder = dsp_path = os.path.normpath(output_path) + "_dsp_tmp/"
         try:
             if rescale == 32:
                 post.downsample_dataset(input_path, dsp_path, downsample,

examples/common_data_processing_workflow/find_center_manual_cli.py

@@ -121,10 +121,10 @@
 
 if view == "sino":
     util.find_center_visual_sinograms(sinogram, output_base, start_center,
-                                      stop_center, step=1, zoom=1.0)
+                                      stop_center, step=step_center, zoom=1.0)
 else:
     util.find_center_visual_slices(sinogram, output_base, start_center,
-                                   stop_center, 1.0, zoom=1.0, method=method,
+                                   stop_center, step_center, zoom=1.0, method=method,
                                    gpu=False, angles=angles, ratio=1.0,
                                    filter_name=None)
 

examples/common_data_processing_workflow/minimal_scripts/data_reduction.py

@@ -0,0 +1,126 @@
+import os
+import shutil
+import time
+import timeit
+import numpy as np
+import h5py
+
+import algotom.io.loadersaver as losa
+import algotom.post.postprocessing as post
+
+
+"""
+This script is used for data reduction of a reconstructed volume: 
+rescaling, downsampling, cropping, reslicing.
+"""
+
+input_path = "/tomography/scan_00008/full_reconstruction.hdf"
+# input_path = "/tomography/scan_00008/full_reconstruction/" # For tifs
+
+output_path = "/tomography/tmp/scan_00008/data_reduction"
+
+crop = (0, 0, 0, 0, 0, 0)  #  (top, bottom, left, right, front, back)
+rescale = 8  # 8-bit
+downsample = 2  # Downsample
+reslice = 1  # Reslice along axis-1
+rotate_angle = 0.0  # Rotate the volume is reslice != 0
+hdf_key = "entry/data/data"
+
+print("\n====================================================================")
+print("          Run the script for data reduction")
+print("          Time: {}".format(time.ctime(time.time())))
+print("====================================================================\n")
+print("Crop-(top, bottom, left, right, front, back) = {}".format(crop))
+print("Rescale: {}-bit".format(rescale))
+print("Downsample: {}".format(downsample))
+print("Reslice: axis-{}".format(reslice))
+print("Rotate: {}-degree".format(rotate_angle))
+print("...")
+print("...")
+print("...")
+
+# Check output pathw
+if output_path.lower().endswith(('.hdf', '.nxs', '.h5')):
+    output_format = "hdf"
+    if os.path.isfile(output_path):
+        output_path = losa.make_file_name(output_path)
+else:
+    _, file_extension = os.path.splitext(output_path)
+    if file_extension != "":
+        raise ValueError("Output-path must be a folder or a hdf/nxs/h5 file!!!")
+    else:
+        output_format = "folder"
+        if os.path.isdir(output_path):
+            num = 0
+            output_path = os.path.normpath(output_path)
+            parent_path = os.path.dirname(output_path)
+            last_folder = os.path.basename(output_path)
+            while True:
+                name = ("0000" + str(num))[-5:]
+                new_path = parent_path + "/" + last_folder + "_" + name
+                if os.path.isdir(new_path):
+                    num = num + 1
+                else:
+                    break
+            output_path = new_path
+
+t_start = timeit.default_timer()
+downsample = int(np.clip(downsample, 1, 30))
+if reslice == 0:
+    if downsample > 1:
+        if rescale == 32:
+            post.downsample_dataset(input_path, output_path, downsample, key_path=hdf_key,
+                                    method='mean', rescaling=False,
+                                    skip=None, crop=crop)
+        else:
+            post.downsample_dataset(input_path, output_path, downsample, key_path=hdf_key,
+                                    method='mean', rescaling=True, nbit=rescale,
+                                    skip=None,crop=crop)
+    else:
+        post.rescale_dataset(input_path, output_path, nbit=rescale, crop=crop,
+                                key_path=hdf_key)
+else:
+    if downsample == 1:
+        if rescale == 32:
+            post.reslice_dataset(input_path, output_path, rescaling=False, key_path=hdf_key,
+                                 rotate=rotate_angle, axis=reslice, crop=crop,
+                                 chunk=40, show_progress=True, ncore=None)
+        else:
+            post.reslice_dataset(input_path, output_path, rescaling=True, key_path=hdf_key,
+                                 rotate=rotate_angle, nbit=rescale, axis=reslice, crop=crop,
+                                 chunk=40, show_progress=True, ncore=None)
+    else:
+        if output_format == "hdf":
+            file_name = os.path.basename(output_path)
+            folder_path = os.path.splitext(output_path)[0]
+            dsp_folder = folder_path + "_dsp_tmp/"
+            dsp_path = dsp_folder + file_name
+        else:
+            dsp_folder = dsp_path = os.path.normpath(output_path) + "_dsp_tmp/"
+        try:
+            if rescale == 32:
+                post.downsample_dataset(input_path, dsp_path, downsample, key_path=hdf_key,
+                                        method='mean', rescaling=False,
+                                        skip=None, crop=crop)
+            else:
+                post.downsample_dataset(input_path, dsp_path, downsample, key_path=hdf_key,
+                                        method='mean', rescaling=True, nbit=rescale,
+                                        skip=None,crop=crop)
+            if rescale == 32:
+                post.reslice_dataset(dsp_path, output_path, rescaling=False, key_path="entry/data",
+                                    rotate=rotate_angle, axis=reslice,
+                                    chunk=40, show_progress=True, ncore=None)
+            else:
+                post.reslice_dataset(dsp_path, output_path, rescaling=True, key_path="entry/data",
+                                    rotate=rotate_angle, nbit=rescale, axis=reslice,
+                                    chunk=40, show_progress=True, ncore=None)
+            shutil.rmtree(dsp_folder)
+        except Exception as e:
+            shutil.rmtree(dsp_folder)
+            raise ValueError(e)
+
+t_stop = timeit.default_timer()
+print("\n====================================================================")
+print("Output: {}".format(output_path))
+print("!!! All Done. Total time cost {} !!!".format(t_stop - t_start))
+print("\n====================================================================")

examples/common_data_processing_workflow/minimal_scripts/find_center_manual.py

@@ -0,0 +1,87 @@
+import time
+import timeit
+import numpy as np
+import algotom.io.loadersaver as losa
+import algotom.prep.correction as corr
+import algotom.prep.removal as remo
+import algotom.prep.filtering as filt
+import algotom.util.utility as util
+
+"""
+This script is used to find the center of rotation manually:
+https://algotom.readthedocs.io/en/latest/toc/section4/section4_5.html#finding-the-center-of-rotation
+"""
+
+proj_file = "/tomography/raw_data/scan_00008/scan_00008.nxs"
+flat_file = "/tomography/raw_data/scan_00009/flat_00000.hdf"
+dark_file = "/tomography/raw_data/scan_00009/dark_00000.hdf"
+
+output_base = "/tomography/tmp/scan_00008/find_center/"
+
+slice_index = 1000
+
+start_center = 1550
+stop_center = 1650
+step_center = 2
+crop_left = 0
+crop_right = 0
+
+t_start = timeit.default_timer()
+print("\n====================================================================")
+print("            Run the script for finding center manually")
+print("            Time: {}".format(time.ctime(time.time())))
+print("====================================================================\n")
+
+# Keys to hdf/nxs/h5 datasets
+proj_key = "entry/data/data"
+flat_key = "entry/data/data"
+dark_key = "entry/data/data"
+angle_key = "entry/data/rotation_angle"
+
+# Load data, average flat and dark images
+proj_obj = losa.load_hdf(proj_file, proj_key)  # hdf object
+(depth, height, width) = proj_obj.shape
+left = crop_left
+right = width - crop_right
+width1 = right - left
+angles = np.deg2rad(losa.load_hdf(proj_file, angle_key)[:])
+flat_field = np.mean(np.asarray(losa.load_hdf(flat_file, flat_key)), axis=0)
+dark_field = np.mean(np.asarray(losa.load_hdf(dark_file, dark_key)), axis=0)
+
+if slice_index < 0 or slice_index > height - 1:
+    raise ValueError(f"Index is out of the range [0, {height - 1}]")
+if start_center < 0 or start_center > width1 - 1:
+    raise ValueError(f"Incorrect starting value, given image-width: {width1}")
+if stop_center < 1 or stop_center > width1:
+    raise ValueError(f"Incorrect stopping value, given image-width: {width1}")
+
+sinogram = corr.flat_field_correction(proj_obj[:, slice_index, left:right],
+                                      flat_field[slice_index, left:right],
+                                      dark_field[slice_index, left:right])
+# Apply zinger removal
+# sinogram = remo.remove_zinger(sinogram, 0.08)
+
+# Apply ring removal
+sinogram = remo.remove_stripe_based_normalization(sinogram, 15)
+# sinogram = remo.remove_stripe_based_sorting(sinogram, 21)
+# sinogram = remo.remove_all_stripe(sinogram, 2.0, 51, 21)
+
+# Apply contrast enhancement
+sinogram = filt.fresnel_filter(sinogram, 100)
+
+
+## Visual using sinogram
+# util.find_center_visual_sinograms(sinogram, output_base, start_center,
+#                                     stop_center, step=step_center, zoom=1.0)
+
+## Visual using reconstructed image
+util.find_center_visual_slices(sinogram, output_base, start_center,
+                               stop_center, step_center, zoom=1.0,
+                               method="gridrec", gpu=False, angles=angles,
+                               ratio=1.0, filter_name=None)
+
+t_stop = timeit.default_timer()
+print("====================================================================\n")
+print("All done! Time cost {}".format(t_stop - t_start))
+print("Output: {}".format(output_base))
+print("====================================================================\n")