Commit c0a2fdb4 authored by Pierre Paleo's avatar Pierre Paleo

Merge branch 'release_2020.5.0' into 'master'

Prepare release 2020.5.0

Closes #204, #206, #207, and #205

See merge request !98
parents 57425323 fde07d4a
Pipeline #37577 passed with stages
in 5 minutes and 17 seconds
# Change Log
## 2020.5.0
Major version (November 2020).
### Application
- Add cuda backend for histogram, also for the `nabu-histogram` command.
- Enable histogram computation for other output formats than HDF5 (ex. tiff)
- The files "nabu_processes.h5" and "nabu.log" are now named as a function of the dataset prefix. This prevents conflicts when several datasets are reconstruction simultaneously from the same folder.
- Add nabu-generate-info CLI tool, which generates a ".info" for pluggin ESRF legacy pipeline post-processing tools.
- Add a utility to exclude projections from being processed: `exclude_projections` in `[dataset]`
- Add new methods for finding the center of rotation
- Mitigate issue of nabu hanging forever to allocate memory on big data chunks
- Fix unsharp mask when used without phase retrieval
### Library
- Fix half-tomography sinograms when center of rotation is on the left
- Add misc.histogram_cuda: Cuda backend for partial volume histograms
- Add `preproc.alignment.CenterOfRotationAdaptiveSearch`, `preproc.alignment.CenterOfRotationSlidingWindow` and `preproc.alignment.CenterOfRotationGrowingWindow` for automatic CoR estimation.
## 2020.4.2
This version brings new features and fixes.
......@@ -88,4 +111,4 @@ This is a "preview release", preparing an important release foreseen for Septemb
## 2020.1.0: 2020/04/29
This is the first official release of Nabu. Please see the [features overview](www.silx.org/pub/nabu/doc/features.html)
\ No newline at end of file
This is the first official release of Nabu. Please see the [features overview](www.silx.org/pub/nabu/doc/features.html)
nabu.resources.cli.validate module
==================================
nabu.misc.histogram\_cuda module
================================
.. automodule:: nabu.resources.cli.validate
.. automodule:: nabu.misc.histogram_cuda
:members:
:undoc-members:
:show-inheritance:
......@@ -10,6 +10,7 @@ Submodules
nabu.misc.binning
nabu.misc.fourier_filters
nabu.misc.histogram
nabu.misc.histogram_cuda
nabu.misc.unsharp
nabu.misc.unsharp_cuda
nabu.misc.unsharp_opencl
......
nabu.resources.cli.generate\_header module
==========================================
.. automodule:: nabu.resources.cli.generate_header
:members:
:undoc-members:
:show-inheritance:
......@@ -9,11 +9,11 @@ Submodules
nabu.resources.cli.bootstrap
nabu.resources.cli.cli_configs
nabu.resources.cli.generate_header
nabu.resources.cli.histogram
nabu.resources.cli.nx_z_splitter
nabu.resources.cli.reconstruct
nabu.resources.cli.utils
nabu.resources.cli.validate
Module contents
---------------
......
......@@ -162,12 +162,13 @@ Mind that jpeg2000 support currently has the following limitations:
Nabu provides a method to find the half-shift between two images. The center of axial vertical rotation is obtained when the fist image is a radiography at the rotation angle 0 and the second image is given by the radiography at the rotation angle 180 after flipping the image horizontally. The rotation axis position is the center of the image plus the found shift.
Configuration file: section `[reconstruction]`, key `rotation_axis_position`. Values can be:
- A number (known CoR)
- Empty: the CoR is set to the middle of the detector
- `centered` (or `auto`): this searches for a CoR around the center, but does not work for half-acquisition
- `global`: new method, should return a CoR in any setting.
- Empty (default): the CoR is set to the middle of the detector: `(detector_width - 1)/2.0`
- A number (known CoR)
- `centered`: a fast and simple auto-CoR method. It only works when the CoR is not far from the middle of the detector. It does not work for half-tomography.
- `global`: a slow but robust auto-CoR.
- `sliding-window`: semi-automatically find the CoR with a sliding window. You have to specify on which side the CoR is (left, center, right).
- `growing-window` : automatically find the CoR with a sliding-and-growing window. You can tune the option with the parameter 'cor_options'.
API: [CenterOfRotation](apidoc/nabu.preproc.alignment.rst#nabu.preproc.alignment.CenterOfRotation)
......
......@@ -16,6 +16,7 @@ Quick links:
Latest news:
* November 2020: `Version 2020.5.0 released <v2020_5_0.md>`_
* October 2020: `Version 2020.4.0 released <v2020_4_0.md>`_
* August 2020: `Version 2020.3.0 released <v2020_3_0.md>`_
* June 2020: `Version 2020.2.0 released <v2020_2_0.md>`_
......@@ -83,6 +84,8 @@ Release history
.. toctree::
:maxdepth: 3
v2020_5_0.md
v2020_4_0.md
v2020_3_0.md
v2020_2_0.md
......
# Version 2020.5.0
Version 2020.5.0 is a major version aiming at stabilizing fixes and features added over weeks.
``` note:: The changelog is available at https://gitlab.esrf.fr/tomotools/nabu/-/blob/master/CHANGELOG.md
```
## Highlights
This section highlights some of the available [features](features.md).
### New methods for estimating the center of rotation
In the configuration, nabu now provides four methods for estimating the Center Of Rotation (CoR):
- `centered`: a fast and simple auto-CoR method. It only works when the CoR is not far from the middle of the detector. It does not work for half-tomography.
- `global`: a slow but robust auto-CoR.
- `sliding-window`: semi-automatically find the CoR with a sliding window. You have to specify on which side the CoR is (left, center, right).
- `growing-window` : automatically find the CoR with a sliding-and-growing window. You can tune the option with the parameter 'cor_options'.
Advanced options for each of these methods can be tuned with the `cor_options` parameter.
### Histogram improvements
- When the slice/volume histogram is computed, the processing is now done on GPU, dramatically speeding up this step.
- The histogram can now be computed for output file formats other than HDF5.
### Excluding projections
Sometimes it is useful to exclude certain projection images from being processes. It can be done with two different means:
- The CLI tool `nxtomomill patch-nx` (available since `nxtomomill 0.4.0`) can patch a HDF5-Nexus dataset inplace.
- In the nabu configuration, you can provide `exclude_projections = excluded.txt` were the text file contains one radio index to ignore per line.
## Changes
### Unique files `nabu_processes.h5` and `nabu.log`
The files `nabu_processes.h5` (containing various processing steps dumped for further reuse) and `nabu.log` (default file name for logs) are now named `dataset_prefix_nabu_processes.h5` and `dataset_prefix_nabu.log`.
The reason is, sometimes several reconstruction are launched simultaneously on dataset files placed in the same folder, giving conflicts in both files.
### Use of "magnified pixel size"
The HDF5-NX file used as an input contains two pixel sizes: detector pixel size, and "magnified pixel size" accounting for optics effects. Until version `2020.4.1`, the former was taken for scaling the reconstruction and computing the Paganin filter. From `2020.4.1` and on, the magnified pixel size will be used.
This has mainly two effects:
- For phase retrieval, this change impacts the "blurring effect" but not the gray values scale.
When using a smaller pixel size, the delta/beta value will have to be much smaller to obtain a result similar to the output of previous versions. More precisely, if `magnified_pixel_size = pixel_size/10`, then a result obtained with the parameters `(pixel_size, delta_beta)` using `nabu < 2020.4.1` can be obtained by using parameters `(magnified_pixel_size, delta_beta/100)` in `nabu >= 2020.4.1`.
- For reconstruction, this change impacts the final gray values range, as the reconstructions values are normalized with (horizontal) pixel size.
\ No newline at end of file
__version__ = "2020.4.2-dev"
__version__ = "2020.5.0-beta2"
__nabu_modules__ = [
"app",
"cuda",
......
......@@ -141,7 +141,7 @@ class FullFieldPipeline:
we need to load extra data around the edges of each image. Otherwise,
a default padding type is applied.
"""
if "phase" not in self.processing_steps:
if not(self.use_radio_processing_margin):
self._phase_margin = None
return
up_margin = self._phase_margin_up
......@@ -151,8 +151,13 @@ class FullFieldPipeline:
self._phase_margin = ((up_margin, down_margin), (left_margin, right_margin))
@property
def use_radio_processing_margin(self):
return ("phase" in self.processing_steps) or ("unsharp_mask" in self.processing_steps)
def _get_phase_margin(self):
if "phase" not in self.processing_steps:
if not(self.use_radio_processing_margin):
return ((0, 0), (0, 0))
return self._phase_margin
......@@ -263,7 +268,7 @@ class FullFieldPipeline:
return shape
def _get_phase_output_shape(self):
if "phase" not in self.processing_steps:
if not(self.use_radio_processing_margin):
self._radios_cropped_shape = self.radios.shape
return
((up_margin, down_margin), (left_margin, right_margin)) = self._phase_margin
......@@ -281,7 +286,7 @@ class FullFieldPipeline:
def _allocate_recs(self, ny, nx):
self.n_slices = self.radios.shape[1] # TODO modify with vertical shifts
if "phase" in self.processing_steps:
if self.use_radio_processing_margin:
self.n_slices -= sum(self.phase_margin[0])
self.recs = self._allocate_array((self.n_slices, ny, nx), "f", name="recs")
......@@ -454,10 +459,10 @@ class FullFieldPipeline:
@use_options("reconstruction", "reconstruction")
def _prepare_reconstruction(self):
options = self.processing_options["reconstruction"]
x_s, x_e = options["start_x"], options["end_x"]+1
y_s, y_e = options["start_y"], options["end_y"]+1
self._rec_roi = (x_s, x_e, y_s, y_e)
self._allocate_recs(y_e - y_s, x_e - x_s)
x_s, x_e = options["start_x"], options["end_x"]
y_s, y_e = options["start_y"], options["end_y"]
self._rec_roi = (x_s, x_e + 1, y_s, y_e + 1)
self._allocate_recs(y_e - y_s + 1, x_e - x_s + 1)
@use_options("reconstruction", "reconstruction")
......
......@@ -235,7 +235,7 @@ class CudaFullFieldPipelineLimitedMemory(CudaFullFieldPipeline):
return shape
def _get_phase_output_shape(self):
if "phase" not in self.processing_steps:
if not(self.use_radio_processing_margin):
self._radios_cropped_shape = self.radios_group_shape
return
((up_margin, down_margin), (left_margin, right_margin)) = self._phase_margin
......
......@@ -106,8 +106,9 @@ class LocalReconstruction:
def _compute_phase_margin(self):
unsharp_margin = self._compute_unsharp_margin()
if "phase" not in self.process_config.processing_steps:
self._phase_margin = (0, 0)
self._phase_margin = unsharp_margin
self._margin_v = self._phase_margin[0]
return
radio_shape = self.process_config.dataset_infos.radio_dims[::-1]
......@@ -124,11 +125,24 @@ class LocalReconstruction:
pixel_size=opts["pixel_size_microns"],
padding=opts["padding_type"]
)
margin_v = max(margin_v, unsharp_margin[0])
self._phase_margin = (margin_v, margin_h)
self._margin_v = self._phase_margin[0]
self.logger.info("Estimated phase margin: %d pixels" % self._margin_v)
def _compute_unsharp_margin(self):
if "unsharp_mask" not in self.process_config.processing_steps:
return (0, 0)
opts = self.process_config.processing_options["unsharp_mask"]
sigma = opts["unsharp_sigma"]
# nabu uses cutoff = 4
cutoff = 4
gaussian_kernel_size = int(ceil(2 * cutoff * sigma + 1))
self.logger.debug("Unsharp mask margin: %d pixels" % gaussian_kernel_size)
return (gaussian_kernel_size, gaussian_kernel_size)
def _get_pipeline_class(self):
self._limited_mem = False
# Actually less in some cases (margin_far_up + margin_far_down instead of 2*margin_v).
......@@ -313,8 +327,14 @@ class LocalReconstruction:
of nabu config will be taken.
"""
out_cfg = self.process_config.nabu_config["output"]
out_dir = out_cfg["location"]
# prevent issue when out_dir is empty, which happens only if dataset/location is a relative path.
# this should be prevented earlier
if out_dir is None or len(out_dir.strip()) == 0:
out_dir = dirname(dirname(self.results[list(self.results.keys())[0]]))
#
if output_file is None:
output_file = join(out_cfg["location"], out_cfg["file_prefix"]) + ".hdf5"
output_file = join(out_dir, out_cfg["file_prefix"]) + ".hdf5"
if isfile(output_file):
msg = str("File %s already exists" % output_file)
if out_cfg["overwrite_results"]:
......@@ -346,7 +366,7 @@ class LocalReconstruction:
"nabu_config": self.process_config.nabu_config,
"processing_options": self.process_config.processing_options,
},
base_dir=out_cfg["location"],
base_dir=out_dir,
overwrite=out_cfg["overwrite_results"]
)
self.merge_histograms(output_file=output_file)
......
......@@ -161,6 +161,12 @@ class CudaKernel(object):
def call(self, *args, **kwargs):
block, grid = self.get_block_grid(*args, **kwargs)
# pycuda crashes when any element of block/grid is not a python int (ex. numpy.int64).
# A weird behavior once observed is "data.shape" returning (np.int64, int, int) (!).
# Ensure that everything is a python integer.
block = tuple(int(x) for x in block)
grid = tuple(int(x) for x in grid)
#
args = self.follow_gpu_arr(args)
if self.prepared:
......
......@@ -117,3 +117,14 @@ def get_h5_value(fname, h5_path, default_ret=None):
except KeyError:
val_ptr = default_ret
return val_ptr
def get_h5_str_value(dataset_ptr):
"""
Get a HDF5 field which can be bytes or str (depending on h5py version !).
"""
data = dataset_ptr[()]
if isinstance(data, str):
return data
else:
return bytes.decode(data)
......@@ -7,7 +7,7 @@ computations.py: determine computational needs, chunking method to be used, etc.
from silx.image.tomography import get_next_power
def estimate_required_memory(process_config, chunk_size=None, radios_and_slices=True):
def estimate_required_memory(process_config, chunk_size=None, radios_and_slices=True, warn_from_GB=None):
"""
Estimate the memory (RAM) needed for a reconstruction.
......@@ -15,6 +15,17 @@ def estimate_required_memory(process_config, chunk_size=None, radios_and_slices=
-----------
radios_and_slices: bool
Whether radios and slices will co-exist in memory (meaning more memory usage)
warn_from_GB: float, optional
Amount of memory in GB from where a warning flag will be raised.
Default is None
If set to a number, the result will be in the form (estimated_memory_GB, warning)
where 'warning' is a boolean indicating whether memory allocation/transfer might be problematic.
Notes
-----
It seems that Cuda does not allow allocating and/or transferring more than 16384 MiB (17.18 GB).
If warn_from_GB is not None, then the result is in the form (estimated_memory_GB, warning)
where warning is a boolean indicating wheher memory allocation/transfer might be problematic.
"""
dataset = process_config.dataset_infos
nabu_config = process_config.nabu_config
......@@ -25,6 +36,7 @@ def estimate_required_memory(process_config, chunk_size=None, radios_and_slices=
Ny = chunk_size
total_memory_needed = 0
memory_warning = False
# Read data
# ----------
......@@ -33,6 +45,8 @@ def estimate_required_memory(process_config, chunk_size=None, radios_and_slices=
data_volume_size = Nx * Ny * Na * 4
data_image_size = Nx * Ny * 4
total_memory_needed += data_volume_size
if (warn_from_GB is not None) and data_volume_size/1e9 > warn_from_GB:
memory_warning = True
# CCD processing
# ---------------
......@@ -82,7 +96,10 @@ def estimate_required_memory(process_config, chunk_size=None, radios_and_slices=
reconstructed_volume_size = Nx_rec * Ny_rec * Nz_rec * 4 # float32
total_memory_needed += reconstructed_volume_size
return total_memory_needed
if warn_from_GB is None:
return total_memory_needed
else:
return (total_memory_needed, memory_warning)
def estimate_chunk_size(available_memory_GB, process_config, chunk_step=50):
"""
......@@ -109,11 +126,14 @@ def estimate_chunk_size(available_memory_GB, process_config, chunk_step=50):
chunk_size = chunk_step
last_good_chunk_size = chunk_size
while True:
required_mem = estimate_required_memory(
process_config, chunk_size=chunk_size, radios_and_slices=radios_and_slices
(required_mem, mem_warn) = estimate_required_memory(
process_config,
chunk_size=chunk_size,
radios_and_slices=radios_and_slices,
warn_from_GB=17 # 2**32 elements - see estimate_required_memory docstring note
)
required_mem_GB = required_mem / 1e9
if required_mem_GB > available_memory_GB or chunk_size > max_dz:
if required_mem_GB > available_memory_GB or chunk_size > max_dz or mem_warn:
break
last_good_chunk_size = chunk_size
chunk_size += chunk_step
......
......@@ -50,12 +50,13 @@ class NabuValidator(object):
return res
def _get_nx_ny(self):
def _get_nx_ny(self, binning=1):
if self.nabu_config["reconstruction"]["enable_halftomo"]:
cor = int(round(self.dataset_infos.axis_position))
nx = max(2*cor, 2 * (self.dataset_infos.radio_dims[0] - 1 - cor))
cor = int(round(self.dataset_infos.axis_position / binning))
nx = self.dataset_infos.radio_dims[0] // binning
nx = max(2*cor, 2 * (nx - 1 - cor))
else:
nx, nz = self.dataset_infos.radio_dims
nx = self.dataset_infos.radio_dims[0] // binning
ny = nx
return nx, ny
......@@ -270,22 +271,22 @@ class NabuValidator(object):
self.dataset_infos.n_angles //= subsampling_factor
if self.binning != (1, 1):
bin_x, bin_z = self.binning
bin_y = bin_x # square slices
# Update end_x, end_y, end_z
rec_cfg = self.nabu_config["reconstruction"]
end_x, end_y = self.get_end_xy() # Not so trivial. See function documentation
rec_cfg["end_x"] = end_x
rec_cfg["end_y"] = end_y
rec_cfg["end_z"] = (rec_cfg["end_z"] + 1) // bin_z - 1
rec_cfg["start_x"] = (rec_cfg["start_x"] // bin_x)
rec_cfg["start_y"] = (rec_cfg["start_y"] // bin_y)
rec_cfg["start_z"] = (rec_cfg["start_z"] // bin_z)
# Update radio_dims
nx, nz = self.dataset_infos.radio_dims
nx //= bin_x
ny = nx # square slices
nz //= bin_z
self.dataset_infos.radio_dims = (nx, nz)
rec_cfg = self.nabu_config["reconstruction"]
end_x = rec_cfg["end_x"]
end_y = rec_cfg["end_y"]
end_z = rec_cfg["end_z"]
bin_y = bin_x # square slices
rec_cfg["end_x"] = (end_x // bin_x) - ((nx % bin_x) != 0)
rec_cfg["end_y"] = (end_y // bin_y) - ((ny % bin_y) != 0)
rec_cfg["end_z"] = (end_z // bin_z) - ((nz % bin_z) != 0)
rec_cfg["start_x"] = (rec_cfg["start_x"] // bin_x)
rec_cfg["start_y"] = (rec_cfg["start_y"] // bin_y)
rec_cfg["start_z"] = (rec_cfg["start_z"] // bin_z)
# Update the Rotation center
# TODO axis_corrections
rot_c = self.dataset_infos.axis_position
......@@ -316,3 +317,74 @@ class NabuValidator(object):
# TODO handle other modes
def get_end_xy(self):
"""
Return the "end_x" value for reconstruction, accounting for binning.
There are three situations:
1. Normal setting: Nx_rec = Nx
2. Half acquisition, CoR on the right side: Nx_rec = 2 * |c|
3. Half acquisition, CoR on the left side: Nx_rec = 2 * (Nx - 1 - |c|)
where |c| = int(round(c)).
**Without binnig**
Let e0 denote the default value for "end_x", without user modification.
By default, all the slice is reconstructed, so
e0 = Nx_rec - 1
Let e denote the user value for "end_x". By default e = e0, but the user might
want to tune it.
Let d denote the distance between e and e0: d = e0 - e. By default d = 0.
**With binning**
Let b denote the binning value in x.
1. Nx' = Nx // b
2. Nx' = 2 * |c/b|
3. Nx' = 2 * (Nx//b - 1 - |c/b|)
With the same notations, with a prime suffix, we have:
* e0' = Nx' - 1 is the default value of "end_x" accounting for binning
* e' is the user value for "end_x" accounting for binning
* d' = e0' - e' is the distance between e0' and e'
In the standard setting (no half tomography), computing e' (i.e end_x) is straightforward:
e' = (e+1)//b - 1
With half tomography, because of the presence of |c| = int(floor(c)), the things
are a bit more difficult.
We enforce the following invariant in all settings:
(I1) dist(e0', e') = dist(e0, e) // b
Another possible invariant is
(I2) delta_x' = delta_x // b
Which is equivalent to (I1) only in setting (1) (i.e no half-tomography).
In the half-tomography setting, (I1) and (I2) are not equivalent anymore, so we have
to choose between the two.
We believe it is more consistent to preserve "end_x", so that a user not modying "end_x"
ends up with all the range [0, Nx - 1] reconstructed.
Therefore:
e' = e0' - d//b
"""
b, _ = self.binning
end_xy = []
for i in range(2):
what = ["end_x", "end_y"][i]
e0 = self._get_nx_ny()[i] - 1
e0p = self._get_nx_ny(binning=b)[i] - 1
d = e0 - self.nabu_config["reconstruction"][what]
ep = e0p - d//b
end_xy.append(ep)
return tuple(end_xy)
......@@ -5,7 +5,7 @@ from tomoscan.io import HDF5File
from tomoscan.esrf.hdf5scan import ImageKey
from ..utils import check_supported
from ..io.writer import NXProcessWriter
from ..io.utils import get_first_hdf5_entry, hdf5_entry_exists
from ..io.utils import get_first_hdf5_entry, hdf5_entry_exists, get_h5_str_value
from .logger import LoggerOrPrint
val_to_nxkey = {
......@@ -30,8 +30,6 @@ def replace_h5_entry(data_url, new_entry):
)
class NXFlatField:
"""
A helper class to load flats and darks, or to compute the final ones.
......@@ -155,11 +153,12 @@ class NXFlatField:
res = True
try:
for key in ["flats_reduction_method", "darks_reduction_method"]:
res &= (my_config[key] == bytes.decode(cfg[key][()]))
res &= (my_config[key] == get_h5_str_value(cfg[key]))
res &= np.allclose(my_config["data_shape"], cfg["data_shape"][()])
res &= np.allclose(my_config["image_key"], cfg["image_key"][()])
du1 = DataUrl(path=my_config["input_file"])
du2 = DataUrl(path=bytes.decode(cfg["input_file"][()]))
du2 = DataUrl(path=get_h5_str_value(cfg["input_file"]))
if not(ignore_filenames):
res &= (
os.path.basename(du1.file_path()) == os.path.basename(du2.file_path())
......
......@@ -69,6 +69,8 @@ class ProcessConfig:
extra_options={
"force_flatfield": self.nabu_config["preproc"]["flatfield_enabled"] == "forced",
"exclude_projections": self.nabu_config["dataset"]["exclude_projections"],
"output_dir": self.nabu_config["output"]["location"],
},
logger=self.logger
)
......
......@@ -31,7 +31,7 @@ def get_values_from_file(fname, n_values=None, shape=None, sep=None, any_size=Fa
"""
if not(any_size) and not((n_values is not None) ^ (shape is not None)):
raise ValueError("Please provide either n_values or shape")
arr = np.loadtxt(fname)
arr = np.loadtxt(fname, ndmin=1)
if (n_values is not None) and (arr.shape[0] != n_values):
raise ValueError("Expected %d values, but could get %d values" % (n_values, arr.shape[0]))
if (shape is not None) and (arr.shape != shape):
......
......@@ -60,9 +60,8 @@ def setup_package():
'silx >= 0.12.0',
'distributed', # >= 2.10.0',
'dask_jobqueue',
# ~ 'pycuda',
# ~ 'pyopencl',
'tomoscan >= 0.3.2',
'tomoscan >= 0.4.0',
'h5py < 3.0',
],
long_description = """
Nabu - Tomography software
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment