Merge branch 'release_2021.1.0' into 'master'

Release 2021.1.0

Closes #249

See merge request !131

CHANGELOG.md

@@ -21,6 +21,7 @@
 
   - `PaganinPhaseRetrieval` length parameters are now expressed in meters.
   - Add `nabu.estimation` module
+  - The `preproc.FlatField` class now takes arrays as input instead of `DataUrl`. The class `FlatFieldDataUrls` can be used as a replacement.
 
 ## 2020.5.0
 

README.md

@@ -9,16 +9,11 @@ To install the development version:
 ```bash
 pip install [--user] git+https://gitlab.esrf.fr/tomotools/nabu.git
 ```
-The above command should automatically install the following nabu Python dependencies:
+To install the stable version:
 
-- numpy, distributed, pytest
-- For computations acceleration: pycuda, pyopencl
-- silx ([http://www.silx.org/](http://www.silx.org/))
-- tomoscan ([https://gitlab.esrf.fr/tomotools/tomoscan](https://gitlab.esrf.fr/tomotools/tomoscan))
-
-All of them can be installed using `pip`.
-
-Please note that Nabu supports Python >= 3.5.
+```bash
+pip install [--user] nabu
+```
 
 ## Usage
 
@@ -29,7 +24,7 @@ Nabu can be used in several ways:
 
 To get quickly started, launch:
 ```bash
-nabu-config --bootstrap
+nabu-config 
 ```
 Edit the generated configuration file (`nabu.conf`). Then:
 
@@ -41,7 +36,7 @@ will reconstruct the slices 500 to 600, with processing steps depending on `nabu
 
 ## Documentation
 
-The documentation can be found on the silx.org page ([http://www.silx.org/pub/nabu/doc](http://www.silx.org/pub/nabu/doc)).
+The documentation can be found on the silx.org page ([https://www.silx.org/pub/nabu/doc](http://www.silx.org/pub/nabu/doc)).
 
 The latest documentation built by continuous integration can be found here: [https://tomotools.gitlab-pages.esrf.fr/nabu/](https://tomotools.gitlab-pages.esrf.fr/nabu/)
 

doc/features.md

@@ -60,7 +60,7 @@ API: [CCDProcessing.Log](apidoc/nabu.preproc.ccd) and [CudaLog](apidoc/nabu.prep
 
 ### Projections rotation and detector tilt auto-detection
 
-Each rotation can be rotated by a user-defined angle. This can be useful to correct the effects of a detector tilt, or rotation axis tilt in some cases.  
+Each projection can be rotated by a user-defined angle. This can be useful to correct the effects of a detector tilt, or rotation axis tilt in some cases.  
 A command-line interface tool `nabu-rotate` is also available.
 
 Configuration file: `[preproc]`: `rotate_projections` (value in degrees).

nabu/__init__.py

-__version__ = "2020.5.19-dev"
+__version__ = "2021.1.0-rc3"
 __nabu_modules__ = [
     "app",
     "cuda",

nabu/app/fullfield.py

@@ -753,7 +753,8 @@ class FullFieldPipeline:
         el = time() - t0
 
         shp = self.chunk_reader.data.shape
-        GB = np.prod(shp) * self.chunk_reader.data.dtype.itemsize / 1e9
+        itemsize = self.chunk_reader.dtype.itemsize if hasattr(self.chunk_reader, "dtype") else 4
+        GB = np.prod(shp) * itemsize / 1e9
         self.logger.info(
             "Read subvolume %s (%.2f GB) in %.2f s: %.3f GB/s"
             % (str(shp), GB, el, GB/el)

nabu/app/local_reconstruction.py

@@ -486,7 +486,7 @@ class FullFieldReconstructor:
         volume_shape = (
             rec_options["end_z"] - rec_options["start_z"] + 1,
             rec_options["end_y"] - rec_options["start_y"] + 1,
-            rec_options["end_x"] - rec_options["start_x"],
+            rec_options["end_x"] - rec_options["start_x"] + 1,
         )
         writer = NXProcessWriter(
             output_file, entry=output_entry, filemode="a", overwrite=True

nabu/cuda/src/halftomo.cu

@@ -49,49 +49,3 @@ __global__ void halftomo_kernel(
     }
 
 }
-
-
-// Same as previous kernel, but now the CoR is outside the image support
-// i.e rotation_axis_position >= n_x
-// Weigting is different !
-__global__ void halftomo_kernel_cor_outside(
-    float* sinogram,
-    float* output,
-    float* weights,
-    int n_angles,
-    int n_x,
-    int rotation_axis_position
-) {
-    int x = blockDim.x * blockIdx.x + threadIdx.x;
-    int y = blockDim.y * blockIdx.y + threadIdx.y;
-
-    int n_a2 = n_angles / 2;
-    int d = rotation_axis_position - n_x;
-    int n_x2  = 2 * rotation_axis_position;
-
-    if ((x >= n_x2) || (y >= n_a2)) return;
-
-    // output[:, :nx] = sino[:n_a2, :nx]
-    if (x < n_x) {
-        output[y * n_x2 + x] = sinogram[y * n_x + x];
-    }
-
-    // output[:, nx : nx + d] = (1 - weights) * sino[:n_a2, -d :][:, ::-1]
-    else if (x < n_x + d) {
-        float w = weights[x - n_x];
-        // output[y * n_x2 + x] = (1.0f - w) * sinogram[y*n_x + n_x - 1 - (n_x - x)];
-        output[y * n_x2 + x] = (1.0f - w) * sinogram[y*n_x + 2 * n_x - 1 - x];
-    }
-
-    // output[:, nx + d : nx + 2*d] = weights * sino[n_a2:, ::-1][:, :d] = sino[n_a2:, -d-1:-n_x-1:-1]
-    else if (x < n_x + 2 * d) {
-        float w = weights[x - (n_x + d)];
-        output[y * n_x2 + x] = w * sinogram[(n_a2 + y)*n_x + x - 2*d];
-    }
-
-    // output[:, nx+2*d:] = sino[n_a2:, ::-1][:, 2 * d :] = sino[n_a2:, -2*d-1:-n_x-1:-1]
-    else {
-        output[y * n_x2 + x] = sinogram[(n_a2 + y)*n_x + (n_x2 - 1 - x)];
-    }
-
-}
\ No newline at end of file

nabu/estimation/cor_sino.py

@@ -114,7 +114,7 @@ class SinoCor(object):
             self.cor_abs = (imax + window_width - 1.0) / 2
             self.cor_rel = self.cor_abs - self.sx / 2.0 - 1
         if imax < 1:
-            self.logger.warning("sliding width ", window_width, " seems too large!")
+            self.logger.warning(str("sliding width ", window_width, " seems too large!"))
         self.rcor_abs = round(self.cor_abs)
         return self.rcor_abs
 

nabu/preproc/distortion.py

@@ -74,7 +74,6 @@ class DistortionCorrection:
         check_supported(correction_method, self.correction_methods.keys(), "correction method")
         self.corrector = self.correction_methods[correction_method]
         self._corrector_kwargs = correction_kwargs or {}
-        self._corrector_kwargs["fill_value"] = np.nan
 
 
     def estimate_distortion(self, image, reference_image):
@@ -84,8 +83,10 @@ class DistortionCorrection:
 
     def correct_distortion(self, image, coords):
         image_corrected =  self.corrector(image, coords, **self._corrector_kwargs)
-        mask = np.isnan(image_corrected)
-        image_corrected[mask] = image[mask]
+        fill_value = self._corrector_kwargs.get("fill_value", None)
+        if fill_value is not None and np.isnan(fill_value):
+            mask = np.isnan(image_corrected)
+            image_corrected[mask] = image[mask]
         return image_corrected
 
     correct = correct_distortion

nabu/preproc/phase_cuda.py

@@ -91,7 +91,7 @@ class CudaPaganinPhaseRetrieval(PaganinPhaseRetrieval, CudaProcessing):
         ]
         # TODO configurable constant values
         if self.padding == "constant":
-            self.kern_args.extend([ # FIXME self.kernel_args ?
+            self.padding_kernel_args.extend([
                 0, 0, 0, 0
             ])
 

nabu/preproc/sinogram.py

@@ -173,7 +173,7 @@ def convert_halftomo(sino, rotation_axis_position):
     assert (sino.shape[0] % 2) == 0
 
     na, nx = sino.shape
-    if rotation_axis_position >= nx:
+    if rotation_axis_position > nx:
         return _convert_halftomo_right(sino, rotation_axis_position)
 
     na2 = na // 2
@@ -199,24 +199,12 @@ def _convert_halftomo_right(sino, rotation_axis_position):
     setting, with a CoR outside the image support.
     """
     assert sino.ndim == 2
-    assert (sino.shape[0] % 2) == 0
-    assert rotation_axis_position > sino.shape[1]
     na, nx = sino.shape
-    na2 = na // 2
-    r = rotation_axis_position
-    d = r - nx
-    res = np.zeros((na2, 2 * r), dtype="f")
+    assert (na % 2) == 0
+    assert rotation_axis_position > nx
 
-    sino1 = sino[:na2, :]
-    sino2 = sino[na2:, ::-1]
-
-    w = np.linspace(0, 1, d, endpoint=True) # has to be adapted when the sample actually goes beyond this point
-    res[:, :nx] = sino1[:, :]
-    res[:, nx:nx+d] = sino1[:, -d:][:, ::-1] * (1 - w)
-    res[:, nx+d:nx+2*d] = sino2[:, :d][:, ::-1] * w
-    res[:, nx+2*d:] = sino2[:, :]
-
-    return res
+    sino2 = np.pad(sino, ((0, 0), (0, rotation_axis_position - nx)), mode="reflect")
+    return convert_halftomo(sino2, rotation_axis_position)
 
 
 

nabu/preproc/sinogram_cuda.py

 import numpy as np
 import pycuda.gpuarray as garray
 from .sinogram import SinoProcessing, SinoNormalization
+from .sinogram import _convert_halftomo_right # FIXME Temporary patch
 from ..cuda.processing import CudaProcessing
 from ..cuda.kernel import CudaKernel
 from ..utils import get_cuda_srcfile, updiv
@@ -26,10 +27,7 @@ class CudaSinoProcessing(SinoProcessing, CudaProcessing):
     def _init_cuda_halftomo(self):
         if not(self.halftomo):
             return
-        if self._rot_center_int < self.n_x:
-            kernel_name = "halftomo_kernel"
-        else:
-            kernel_name = "halftomo_kernel_cor_outside"
+        kernel_name = "halftomo_kernel"
         self.halftomo_kernel = CudaKernel(
             kernel_name,
             get_cuda_srcfile("halftomo.cu"),
@@ -75,14 +73,8 @@ class CudaSinoProcessing(SinoProcessing, CudaProcessing):
         return output
 
 
-
     # Overwrite parent method
     def _radios_to_sinos_halftomo(self, radios, sinos):
-        # TODO
-        if not(np.isscalar(self.rot_center)):
-            raise NotImplementedError("Half tomo with varying rotation axis position is not implemented yet")
-        #
-
         n_a, n_z, n_x = radios.shape
         n_a2 = n_a // 2
         rc = self._rot_center_int
@@ -98,6 +90,11 @@ class CudaSinoProcessing(SinoProcessing, CudaProcessing):
         else:
             sinos = garray.zeros(self.sinos_halftomo_shape, dtype=np.float32)
 
+        # FIXME: TEMPORARY PATCH, waiting for cuda implementation
+        if self._rot_center_int > self.n_x:
+            return self._radios_to_sinos_halftomo_external_cor(radios, sinos)
+        #
+
         # need to use a contiguous 2D, array otherwise kernel does not work
         d_radio = radios[:, 0, :].copy()
         for i in range(n_z):
@@ -123,6 +120,36 @@ class CudaSinoProcessing(SinoProcessing, CudaProcessing):
         return sinos
 
 
+    def _radios_to_sinos_halftomo_external_cor(self, radios, sinos):
+        """
+        TEMPORARY PATCH waiting to have a cuda implementation
+        Processing is done by getting reach radio on host, which is suboptimal
+        """
+        n_a, n_z, n_x = radios.shape
+        n_a2 = n_a // 2
+        rc = self._rot_center_int
+        out_dwidth = 2 * rc
+
+        # need to use a contiguous 2D, array otherwise kernel does not work
+        sino = radios[:, 0, :].get()
+        for i in range(n_z):
+            radios[:, i, :].get(sino)
+            if self._halftomo_flip:
+                sino = sino[:, ::-1]
+            sino_half = _convert_halftomo_right(sino, self._rot_center_int)
+            sinos[i, :, :] = sino_half[:]
+            if self._halftomo_flip:
+                self.xflip_kernel(
+                    sinos[i], 2*rc, n_a2,
+                    grid=self._xflip_gridsize_2, block=self._xflip_blksize
+                )
+        return sinos
+
+
+
+
+
+
 class CudaSinoNormalization(SinoNormalization, CudaProcessing):
     def __init__(self, kind="chebyshev", sinos_shape=None, radios_shape=None, rot_center=None, halftomo=False, cuda_options=None):
         SinoNormalization.__init__(

nabu/preproc/tests/test_ctf.py

@@ -88,7 +88,7 @@ class TestCtf:
         )
 
         my_flat = scipy.interpolate.interpn(
-            (np.arange(my_flat.shape[1]), np.arange(my_flat.shape[0])),
+            (np.arange(my_flat.shape[0]), np.arange(my_flat.shape[1])),
             my_flat,
             new_coordinates,
             bounds_error=False,
@@ -125,7 +125,7 @@ class TestCtf:
                 "correction_spike_threshold": 3.0
             },
             correction_method="interpn",
-            correction_kwargs=None
+            correction_kwargs={"fill_value": None}
         )
         flats = FlatFieldArrays(
             [1200] + list(self.img_shape_vh),

nabu/resources/cli/bootstrap.py

@@ -2,7 +2,6 @@ from os import path
 from .utils import parse_params_values
 from .cli_configs import BootstrapConfig
 from ...io.config import generate_nabu_configfile
-from ..validators import convert_to_bool
 
 def parse_sections(sections):
     sections = sections.lower()
@@ -38,11 +37,18 @@ def bootstrap():
     prefilled_values = {}
     if args["dataset"] != "":
         prefilled_values["dataset"] = {}
-        prefilled_values["dataset"]["location"] = args["dataset"]
+        user_dataset = args["dataset"]
+        if not path.isabs(user_dataset):
+            user_dataset = path.abspath(user_dataset)
+            print("Warning: using absolute dataset path %s" % user_dataset)
+        if not path.exists(user_dataset):
+            print("Error: cannot find the file or directory %s" % user_dataset)
+            exit(1)
+        prefilled_values["dataset"]["location"] = user_dataset
+
     generate_nabu_configfile(
         args["output"],
         comments=not(no_comments),
         options_level=args["level"],
         prefilled_values=prefilled_values,
     )
-

nabu/resources/cli/reconstruct.py

@@ -5,6 +5,7 @@ from ...io.config import NabuConfigParser, validate_nabu_config
 from .utils import parse_params_values
 from ..utils import is_hdf5_extension
 from .cli_configs import ReconstructConfig
+from ..validators import convert_to_int
 
 
 def update_reconstruction_start_end(conf_dict, user_indices):
@@ -12,24 +13,30 @@ def update_reconstruction_start_end(conf_dict, user_indices):
         return
     rec_cfg = conf_dict["reconstruction"]
     err = None
-    if user_indices in ["first", "middle", "last"]:
+
+    val_int, conv_err = convert_to_int(user_indices)
+    if conv_err is None:
         start_z = user_indices
         end_z = user_indices
-    elif user_indices == "all":
-        start_z = 0
-        end_z = -1
-    elif "-" in user_indices:
-        try:
-            start_z, end_z = user_indices.split("-")
-            start_z = int(start_z)
-            end_z = int(end_z)
-        except Exception as exc:
-            err = "Could not interpret slice indices '%s': %s" % (user_indices, str(exc))
     else:
-        err = "Could not interpret slice indices: %s" % user_indices
-    if err is not None:
-        print(err)
-        exit(1)
+        if user_indices in ["first", "middle", "last"]:
+            start_z = user_indices
+            end_z = user_indices
+        elif user_indices == "all":
+            start_z = 0
+            end_z = -1
+        elif "-" in user_indices:
+            try:
+                start_z, end_z = user_indices.split("-")
+                start_z = int(start_z)
+                end_z = int(end_z)
+            except Exception as exc:
+                err = "Could not interpret slice indices '%s': %s" % (user_indices, str(exc))
+        else:
+            err = "Could not interpret slice indices: %s" % user_indices
+        if err is not None:
+            print(err)
+            exit(1)
     rec_cfg["start_z"] = start_z
     rec_cfg["end_z"] = end_z
 

nabu/resources/dataset_analyzer.py

@@ -287,13 +287,15 @@ class HDF5DatasetAnalyzer(DatasetAnalyzer):
                                                  "tomwer_processes.h5")
         if tomwer_processes_file is None or not (os.path.isfile(tomwer_processes_file)):
                 return False
-        self.logger.info("Loading darks and refs from %s" % tomwer_processes_file)
         new_flats = get_flats_frm_process_file(
             tomwer_processes_file, self.dataset_scanner.entry
         )
         new_darks = get_darks_frm_process_file(
             tomwer_processes_file, self.dataset_scanner.entry
         )
+        if not (len(new_flats) > 0 and len(new_darks) > 0):
+            return False
+        self.logger.info("Loading darks and refs from %s" % tomwer_processes_file)
         self.flats = new_flats
         self.darks = new_darks
         return True

nabu/resources/estimators.py

@@ -260,6 +260,8 @@ class SinoCORFinder:
 
 
     def _check_360(self):
+        if self.dataset_info.dataset_scanner.scan_range == 360:
+            return
         if not is_fullturn_scan(self.dataset_info.rotation_angles):
             raise ValueError("Sinogram-based Center of Rotation estimation can only be used for 360 degrees scans")
 

nabu/resources/nabu_config.py

@@ -212,7 +212,7 @@ nabu_config = {
             "default": "",
             "help": "A file where each line has two values: horizontal and vertical displacements. There should be as many values as there are projection images.",
             "validator": optional_file_location_validator,
-            "type": "optional",
+            "type": "advanced",
         },
         "ctf_advanced_params": {
             "default": "length_scale=1e-5; lim1=1e-5; lim2=0.2; normalize_by_mean=True",

nabu/resources/processconfig.py

 import os
 import posixpath
-from ..utils import copy_dict_items
-from ..io.config import NabuConfigParser, validate_nabu_config
+from ..utils import copy_dict_items, compare_dicts
+from ..io.config import NabuConfigParser, validate_nabu_config, import_h5_to_dict
 from ..io.utils import hdf5_entry_exists, get_h5_value
 from .dataset_analyzer import analyze_dataset, DatasetAnalyzer
 from .dataset_validator import NabuValidator
@@ -552,14 +552,15 @@ class ProcessConfig:
         # TODO add more tests to compare configurations
         if self.resume_from_step is None:
             return None
-        # Get file dataset shape/start_z/end_z
-        cfg_h5_path = posixpath.join(
+
+        # Check dataset shape/start_z/end_z
+        rec_cfg_h5_path = posixpath.join(
             self.dataset_infos.hdf5_entry or "entry",
             self.resume_from_step,
             "configuration/nabu_config/reconstruction"
         )
-        dump_start_z = get_h5_value(process_file, posixpath.join(cfg_h5_path, "start_z"))
-        dump_end_z = get_h5_value(process_file, posixpath.join(cfg_h5_path, "end_z"))
+        dump_start_z = get_h5_value(process_file, posixpath.join(rec_cfg_h5_path, "start_z"))
+        dump_end_z = get_h5_value(process_file, posixpath.join(rec_cfg_h5_path, "end_z"))
         start_z, end_z = self.nabu_config["reconstruction"]["start_z"], self.nabu_config["reconstruction"]["end_z"]
         if not (dump_start_z <= start_z and end_z <= dump_end_z):
             msg = "File %s was built with start_z=%d, end_z=%d but current configuration asks for start_z=%d, end_z=%d" % (process_file, dump_start_z, dump_end_z, start_z, end_z)
@@ -568,6 +569,32 @@ class ProcessConfig:
                 return None
             self.logger.fatal(msg)
             raise ValueError(msg)
+
+        # Check parameters other than reconstruction
+        filedump_nabu_config = import_h5_to_dict(
+            process_file,
+            posixpath.join(
+                self.dataset_infos.hdf5_entry or "entry",
+                self.resume_from_step,
+                "configuration/nabu_config"
+            )
+        )
+        sections_to_ignore = ["reconstruction", "output"]
+        for section in sections_to_ignore:
+            filedump_nabu_config[section] = self.nabu_config[section]
+        # special case of the "save_steps process"
+        filedump_nabu_config["pipeline"]["save_steps"] = self.nabu_config["pipeline"]["save_steps"]
+
+        diff = compare_dicts(filedump_nabu_config, self.nabu_config)
+        if diff is not None:
+            msg = "Nabu configuration in file %s differ from the current one: %s" % (process_file, diff)
+            if not raise_on_error:
+                self.logger.error(msg)
+                return None
+            self.logger.fatal(msg)
+            raise ValueError(msg)
+        #
+
         return (dump_start_z, dump_end_z)
 
 

nabu/resources/utils.py

@@ -94,6 +94,8 @@ def extract_parameters(params_str, sep=";"):
     >>> extract_parameters("window_width=None; median_filt_shape=(3,3); padding_mode='wrap'")
     ... {'window_width': None, 'median_filt_shape': (3, 3), 'padding_mode': 'wrap'}
     """
+    if params_str in ("", None):
+        return {}
     params_list = params_str.strip(sep).split(sep)
     res = {}
     for param_str in params_list: