diff --git a/CHANGELOG.md b/CHANGELOG.md index 1231c63e9354ad49a646bf8b08cb0978d1e4a706..8873cfda02b1ccfd33418fae5d9400aa932c7bd2 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -1,5 +1,27 @@ # Change Log +## 2020.4.0 + +This is a version adds a number of features. + +### Application + + - Automatic center of rotation estimation, working for both half-acquisition and regular scans + - Command line tool for splitting a NXTomo file by "z" series + - Volume histogram. Command line tool for merging histogram of multiple volumes. + - Enable to perform flat-field normalization with darks/flats from another dataset + - Sinogram normalization (baseline subtraction) + - Nabu does not need `tomwer_processes.h5` to get the "final" darks/refs anymore. + - Fix `fbp_filter_type = none`: now it will actually disable any filtering + - Fix auto-CoR not working when flatfield is disabled + +### Library + + - Add `misc.histogram` for computing partial histograms and merging them + - Add `preproc.sinogram.SinoNormalization` + - Fix double flat-field when `dff_sigma > 0`  which was giving nonsense results. + - Half-tomography: add support for center of rotation on the left side + ## 2020.3.0 This is a release for ESRF User Service Mode restart. diff --git a/doc/apidoc/nabu.misc.histogram.rst b/doc/apidoc/nabu.misc.histogram.rst new file mode 100644 index 0000000000000000000000000000000000000000..cc3b26fb88280eea48fbfdef60483e080315e10d --- /dev/null +++ b/doc/apidoc/nabu.misc.histogram.rst @@ -0,0 +1,7 @@ +nabu.misc.histogram module +========================== + +.. automodule:: nabu.misc.histogram + :members: + :undoc-members: + :show-inheritance: diff --git a/doc/apidoc/nabu.misc.rst b/doc/apidoc/nabu.misc.rst index 6d99cef25dfca6a2f2bb3b33cc2e63a40af1041b..9279016eb39b27ad9bbeb5e881714e34a3eb5679 100644 --- a/doc/apidoc/nabu.misc.rst +++ b/doc/apidoc/nabu.misc.rst @@ -9,6 +9,7 @@ Submodules nabu.misc.binning nabu.misc.fourier_filters + nabu.misc.histogram nabu.misc.unsharp nabu.misc.unsharp_cuda nabu.misc.unsharp_opencl diff --git a/doc/apidoc/nabu.resources.cli.histogram.rst b/doc/apidoc/nabu.resources.cli.histogram.rst new file mode 100644 index 0000000000000000000000000000000000000000..99358abf326e265d3f447e68ec57fe4f70f37313 --- /dev/null +++ b/doc/apidoc/nabu.resources.cli.histogram.rst @@ -0,0 +1,7 @@ +nabu.resources.cli.histogram module +=================================== + +.. automodule:: nabu.resources.cli.histogram + :members: + :undoc-members: + :show-inheritance: diff --git a/doc/apidoc/nabu.resources.cli.nx_z_splitter.rst b/doc/apidoc/nabu.resources.cli.nx_z_splitter.rst new file mode 100644 index 0000000000000000000000000000000000000000..b56beebe67aa36e19d3897b73c21446e2abb957b --- /dev/null +++ b/doc/apidoc/nabu.resources.cli.nx_z_splitter.rst @@ -0,0 +1,7 @@ +nabu.resources.cli.nx\_z\_splitter module +========================================= + +.. automodule:: nabu.resources.cli.nx_z_splitter + :members: + :undoc-members: + :show-inheritance: diff --git a/doc/apidoc/nabu.resources.cli.rst b/doc/apidoc/nabu.resources.cli.rst index 2842697784673cb56950b333a1ba6aea5c5004bc..c4425da16f237ed815a6b20b2fa4a5d81a75715c 100644 --- a/doc/apidoc/nabu.resources.cli.rst +++ b/doc/apidoc/nabu.resources.cli.rst @@ -9,6 +9,8 @@ Submodules nabu.resources.cli.bootstrap nabu.resources.cli.cli_configs + nabu.resources.cli.histogram + nabu.resources.cli.nx_z_splitter nabu.resources.cli.reconstruct nabu.resources.cli.utils nabu.resources.cli.validate diff --git a/doc/apidoc/nabu.resources.nxflatfield.rst b/doc/apidoc/nabu.resources.nxflatfield.rst new file mode 100644 index 0000000000000000000000000000000000000000..f329f67eeb67d5980a3f15ef4efd1a656c8506a6 --- /dev/null +++ b/doc/apidoc/nabu.resources.nxflatfield.rst @@ -0,0 +1,7 @@ +nabu.resources.nxflatfield module +================================= + +.. automodule:: nabu.resources.nxflatfield + :members: + :undoc-members: + :show-inheritance: diff --git a/doc/apidoc/nabu.resources.rst b/doc/apidoc/nabu.resources.rst index f0f5e6593fa1186c119855ab3351212499ef1f84..a2af0e6e1985aedf11ad027705455e79536a73d4 100644 --- a/doc/apidoc/nabu.resources.rst +++ b/doc/apidoc/nabu.resources.rst @@ -23,6 +23,7 @@ Submodules nabu.resources.logger nabu.resources.machinesdb nabu.resources.nabu_config + nabu.resources.nxflatfield nabu.resources.params nabu.resources.processconfig nabu.resources.utils diff --git a/doc/conf.py b/doc/conf.py index a59960283ad30ce284094093cf3aed331a301f9a..77db0c42f4524aa4d3faffea0f56c25450fc6fc2 100644 --- a/doc/conf.py +++ b/doc/conf.py @@ -22,7 +22,7 @@ copyright = '2020, ESRF' author = 'Pierre Paleo' # The full version, including alpha/beta/rc tags -release = '2020.3.0' +#release = '2020.4.0' # -- General configuration --------------------------------------------------- diff --git a/doc/features.md b/doc/features.md index d96b1a03c55e82b27606842ee17cf5c18f11b499..4ecd543724cc5914bd0ba803b2e9076837435600 100644 --- a/doc/features.md +++ b/doc/features.md @@ -56,6 +56,18 @@ Configuration file: section `[preproc]`: `take_logarithm = 1`, `log_min_clip = 1 API: [CCDProcessing.Log](apidoc/nabu.preproc.ccd.rst#nabu.preproc.ccd.Log) and [CudaLog](apidoc/nabu.preproc.ccd_cuda.rst#nabu.preproc.ccd_cuda.CudaLog) +### Sinogram normalization + +Sinograms can sometimes be "messy" for various reasons. For example, a synchrotron beam refill can take place during the acquisition, and not be compensated properly by flats. +In this case, you can "normalize" the sinogram to get rid of defects. Currently, only a baseline removal is implemented. + +Mind that you probably lose quantativeness by using this additional normalization ! + +Configuration file: `[preproc]` : `sino_normalization = chebyshev`. + +API: [SinoNormalization](apidoc/nabu.preproc.sinogram.rst#nabu.preproc.sinogram.SinoNormalization) + + ## Phase retrieval Phase retrieval is still part of the "pre-processing", although it has a dedicated section in the [configuration file](nabu_config_file). @@ -85,6 +97,7 @@ Configuration file: section `[phase]`: `unsharp_coeff = 1.` and `unsharp_sigma Setting `coeff` to zero (default) disables unsharp masking. + ## Reconstruction Tomographic reconstruction is the process of reconstructing the volume from projections/sinograms. @@ -107,6 +120,16 @@ The purpose of this class is to quickly reconstruct slices over the three main a The Reconstructor enables to reconstruct slices/region of interest without reconstructing the whole volume. +## Post-processing + +### Histogram + +Nabu can compute the histogram of the reconstucted (sub-) volume. As the volume usually does not fit in memory, the histogram is computed by parts, and the final histogram is obtained by merging partial histograms. + +Configuration file: section `[postproc]`: `output_histogram = 1`, `histogram_bins = 1000000`. + +API : [PartialHistogram](apidoc/nabu.misc.histogram.rst#nabu.misc.histogram.PartialHistogram) and [VolumeHistogram](apidoc/nabu.misc.histogram.rst#nabu.misc.histogram.VolumeHistogram) + ## File formats ### HDF5 @@ -118,7 +141,7 @@ When a [multi-stage reconstruction](nabu_cli.md) is performed, the volume is rec ### Tiff Reconstruction can be saved as tiff files by specifying `file_format = tiff` in the configuration `[output]` section. -In the current version (2020.3), tiff support has the following limitations: +Mind that tiff support currently has the following limitations: - One file per slice - Data is saved as `float32` data type, no normalization - No metadata (configuration used to obtain the reconstruction, date, version,...) @@ -126,7 +149,7 @@ In the current version (2020.3), tiff support has the following limitations: ### Jpeg2000 Reconstruction can be saved as jpeg2000 files by specifying `file_format = jpeg2000` in the configuration `[output]` section. -In the current version (2020.3), jpeg2000 support has the following limitations: +Mind that jpeg2000 support currently has the following limitations: - When exporting to `uint16` data type (mandatory for jpeg2000), the normalization from `float32` to `uint16` is done slice-wise instead of volume-wise. This is slightly less accurate. - Only lossless compression is supported. In the future, compression will be tunable through Nabu configuration. - No metadata (configuration used to obtain the reconstruction, date, version,...) @@ -138,6 +161,14 @@ In the current version (2020.3), jpeg2000 support 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. + + + API: [CenterOfRotation](apidoc/nabu.preproc.alignment.rst#nabu.preproc.alignment.CenterOfRotation) ### Detector Translation Along the Beam diff --git a/doc/index.rst b/doc/index.rst index 8ec0aec1083eb278fa98751e162a034b588c377e..d09945d51b36b48929400855ee4937881844e7f4 100644 --- a/doc/index.rst +++ b/doc/index.rst @@ -16,6 +16,7 @@ Quick links: Latest news: + * October 2020: `Version 2020.4.0 released `_ * August 2020: `Version 2020.3.0 released `_ * June 2020: `Version 2020.2.0 released `_ @@ -82,6 +83,7 @@ Release history .. toctree:: :maxdepth: 3 + v2020_4_0.md v2020_3_0.md v2020_2_0.md diff --git a/doc/install.md b/doc/install.md index 6dbd8918d321154f544b0ed1cc1c79862de48f58..48bccf49cb15fb8f631e262037b8140d1d845245 100644 --- a/doc/install.md +++ b/doc/install.md @@ -14,18 +14,8 @@ pip install git+https://gitlab.esrf.fr/tomotools/nabu.git The above `pip install` command should automatically install the nabu dependencies. -For the record, the mandatory dependencies are the following: -- numpy -- pytest -- psutil -- [silx](http://www.silx.org/) -- [tomoscan](https://gitlab.esrf.fr/tomotools/tomoscan) - -All of them can be installed with `pip`. - There are optional dependencides enabling various features: - Computations acceleration (GPU/manycore CPU): `pycuda`, `pyopencl` - Computations distribution: `distributed`, `dask_jobqueue` - Please note that Nabu supports Python >= 3.5. diff --git a/doc/nabu_config_file.md b/doc/nabu_config_file.md index 170f81fa2b53766d7306660b0c5c829d11b43aba..a2fb1244f9e9f3373897dd31ed7187714806e8a8 100644 --- a/doc/nabu_config_file.md +++ b/doc/nabu_config_file.md @@ -52,6 +52,7 @@ Each section describe a usual processing step. In the current version, the avail - `preproc`: phase retrieval, CCD corrections, ... - `phase`: phase retrieval - `reconstruction`: tomography reconstruction +- `postproc`: post-processing: histogram - `output`: output data description - `resources`: computing resources description - `about`: extra information about the current configuration file @@ -71,6 +72,6 @@ Yes. Nabu is foremost a library, meaning that all its component can be accessed ### Compatibility policy -At this early development stage, it is not entirely clear which keys should be in the configuration file. In the future, some keys might be removed, added, or have their name changed. Also, with new features coming, new keys/sections might appear in the configuration file. +During the development of Nabu, some features will be added, leading to new keys in the configuration file. Besides, some keys might be renamed or even deleted if deemed necessary. -Thus, the configuration file might change according to the version of nabu. For this reason, there is a section `[about]` in the configuration file, which defines a `nabu_version` and a `nabu_config_version`. A given version of the nabu software (`nabu_version`) will be compatible with a certain range of nabu configuration files (`nabu_config_version`). +In any case, **a configuration file from an "old" version of nabu will be supported by newer versions** for some time, with a deprecation warning when using obsolete keys. \ No newline at end of file diff --git a/doc/v2020_4_0.md b/doc/v2020_4_0.md new file mode 100644 index 0000000000000000000000000000000000000000..3c482142cd5dad414ad5636230c80833bd9bfa97 --- /dev/null +++ b/doc/v2020_4_0.md @@ -0,0 +1,52 @@ +# Version 2020.4.0 + + +Version 2020.4.0 adds a number of features, notably command line tools for manipulating datasets. + +``` 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). + +### Fully automatic center of rotation estimation + +Various methods exist for estimating automatically the center of rotation (CoR). Nabu had one, but it did not work for half tomography acquisitions. Another method is now available and should work with any kind of acquisition. + +The paramater is `rotation_axis_position` (in section`[reconstruction]`), the value can be: + - A number (known CoR) + - Empty: the CoR is set to the middle of the detector + - `centered` (was `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. + + +### Histograms + +Nabu can now compute a histogram of the reconstructed volume. This histogram can serve for further volume processing (conversion to uint16 for example). The computation is done while the data is still in memory to avoid extraneous disk reads/writes. +The option is available in a new section `[postproc]` with the parameter `output_histogram = 1`. + +A command line tool is also available to merge the histogram of multiple volumes: `nabu-histogram ... `. + +### Use darks/flats from another dataset + +Sometimes, flats cannot be acquired properly during an acquisition. In order to still have the ability to reconstruct the data, a common workaround is to use flats from another dataset. + +Suppose you want to reconstruct a dataset which does not have flats/darks. In this case, nabu will detect that these images are missing and won't do flat-field correction. However you can force the flat-field normalization with `flatfield_enabled = forced` in the `[preproc]` section, and by providing `processes_file = /path/to/nabu_processes.h5`. This `nabu_processes.h5` file (or `tomwer_processes.h5`) can be generated either by tomwer or nabu. + +### Sinogram normalization + +A sinogram normalization is now available in the section `[preproc]` : `sino_normalization = chebyshev`. If enabled, each line of the sinogram undergoes a baseline removal. This can be useful to correct a "refill" taking place during the scan. See [this discussion](https://gitlab.esrf.fr/tomotools/nabu/-/issues/118#note_74284) for more details. + +### Split a NXTomo file by "z series" + +Z-series is the name given to a multi-stage scan, i.e scanning a volume in several stages by moving the sample vertically after each stage (this is different from helical tomography). + +For the time being, the CLI tool `nxtomomill` merges everything in a single file. +Projections have a varying "z" identified by the `entry/sample/z_translation` key (in conformity to the Nexus-Tomo standard). + +The command line tool `nabu-zsplit` enables to split a file into distinct "z" to reconstruct each volume individually. + +``` warning:: This is a temporary solution. This will be natively be handled by nxtomomill in the future. +``` \ No newline at end of file diff --git a/nabu/__init__.py b/nabu/__init__.py index a2e97eadc84aef5a996d5c59f4f1c9c05a2de632..69b4f7f47c18257f3ecc5c9b5f90949f01b093d9 100644 --- a/nabu/__init__.py +++ b/nabu/__init__.py @@ -1,4 +1,4 @@ -__version__ = "2020.3.1" +__version__ = "2020.4.0" __nabu_modules__ = [ "app", "cuda", diff --git a/nabu/app/fullfield.py b/nabu/app/fullfield.py index a3c15ae1ee82eaa715c4b4bc65a9f5ced6e39900..fa365c51f93b665661d073258e76e696db6fbe12 100644 --- a/nabu/app/fullfield.py +++ b/nabu/app/fullfield.py @@ -358,12 +358,13 @@ class FullFieldPipeline: @use_options("double_flatfield", "double_flatfield") def _init_double_flatfield(self): options = self.processing_options["double_flatfield"] + avg_is_on_log = (options["sigma"] is not None) self.double_flatfield = self.DoubleFlatFieldClass( self._get_shape("double_flatfield"), result_url=None, input_is_mlog=False, output_is_mlog=False, - average_is_on_log=False, + average_is_on_log=avg_is_on_log, sigma_filter=options["sigma"] ) @@ -395,7 +396,7 @@ class FullFieldPipeline: margin=margin, ) if "unsharp_mask" not in self.processing_steps: - self.register_callback("phase", self._reshape_radios_after_phase) + self.register_callback("phase", FullFieldPipeline._reshape_radios_after_phase) @use_options("unsharp_mask", "unsharp_mask") def _init_unsharp(self): @@ -405,7 +406,7 @@ class FullFieldPipeline: options["unsharp_sigma"], options["unsharp_coeff"], mode="reflect", method="gaussian" ) - self.register_callback("unsharp_mask", self._reshape_radios_after_phase) + self.register_callback("unsharp_mask", FullFieldPipeline._reshape_radios_after_phase) @use_options("take_log", "mlog") def _init_mlog(self): @@ -458,6 +459,7 @@ class FullFieldPipeline: self._rec_roi = (x_s, x_e, y_s, y_e) self._allocate_recs(y_e - y_s, x_e - x_s) + @use_options("reconstruction", "reconstruction") def _init_reconstruction(self): options = self.processing_options["reconstruction"] @@ -471,6 +473,8 @@ class FullFieldPipeline: rot_center = options["rotation_axis_position"] if options.get("cor_estimated_auto", False): self.logger.info("Estimated center of rotation: %.2f" % rot_center) + if self.sino_builder._halftomo_flip: + rot_center = self.sino_builder.rot_center self.reconstruction = self.FBPClass( self._get_shape("reconstruction"), @@ -484,7 +488,7 @@ class FullFieldPipeline: "axis_correction": options["axis_correction"], } ) - if options["fbp_filter_type"] == "none": + if options["fbp_filter_type"] is None: self.reconstruction.fbp = self.reconstruction.backproj @use_options("histogram", "histogram") diff --git a/nabu/app/fullfield_cuda.py b/nabu/app/fullfield_cuda.py index d800bd8ea669a95f02a8a9efbe508d5f75b44354..29bd36a0f570f955b88d9f35b379e20905c535fd 100644 --- a/nabu/app/fullfield_cuda.py +++ b/nabu/app/fullfield_cuda.py @@ -91,11 +91,11 @@ class CudaFullFieldPipeline(FullFieldPipeline): def _register_callbacks(self): - self.register_callback("read_chunk", self._read_data_callback) + self.register_callback("read_chunk", CudaFullFieldPipeline._read_data_callback) if self.reconstruction is not None: - self.register_callback("reconstruction", self._rec_callback) + self.register_callback("reconstruction", CudaFullFieldPipeline._rec_callback) if self.writer is not None: - self.register_callback("save", self._saving_callback) + self.register_callback("save", CudaFullFieldPipeline._saving_callback) diff --git a/nabu/app/utils.py b/nabu/app/utils.py index 52b5e4dfa5348a2ab72d23ccca892f6185cebea3..d810f45cbcb1a509fd224a6346f19e671837a969 100644 --- a/nabu/app/utils.py +++ b/nabu/app/utils.py @@ -33,7 +33,7 @@ def pipeline_step(step_attr, step_desc): self.logger.debug("End " + step_desc) callback = self._callbacks.get(self._steps_component2name[step_attr], None) if callback is not None: - callback() + callback(self) return res return wrapper return decorator diff --git a/nabu/cuda/medfilt.py b/nabu/cuda/medfilt.py index e3db4ab20387d5d3c4941cc37a0b4947e1ef4a1f..f0432d048d1ebc65c45bf805d1b9d9e376c4b674 100644 --- a/nabu/cuda/medfilt.py +++ b/nabu/cuda/medfilt.py @@ -32,18 +32,18 @@ class MedianFilter(CudaProcessing): Size of the median filter, in the format (y, x). mode: str Boundary handling mode. Available modes are: - "reflect": cba|abcd|dcb - "nearest": aaa|abcd|ddd - "wrap": bcd|abcd|abc - "constant": 000|abcd|000 + - "reflect": cba|abcd|dcb + - "nearest": aaa|abcd|ddd + - "wrap": bcd|abcd|abc + - "constant": 000|abcd|000 Default is "reflect". threshold: float, optional Threshold for the "thresholded median filter". A thresholded median filter only replaces a pixel value by the median if this pixel value is greater or equal than median + threshold. - Cuda Parameters - --------------- + Notes + ------ Please refer to the documentation of the CudaProcessing class for the other parameters. """ diff --git a/nabu/cuda/src/ElementOp.cu b/nabu/cuda/src/ElementOp.cu index 22bd2da195ef1ba4bfc28618d6ac3a89cac76e0f..909a09df8cf16bcc4dec67a92d6ed61df1737356 100644 --- a/nabu/cuda/src/ElementOp.cu +++ b/nabu/cuda/src/ElementOp.cu @@ -86,3 +86,20 @@ __global__ void nlog(float* array, int Nx, int Ny, int Nz, float clip_min, float #endif array[pos] = -logf(val); } + + + +// Reverse elements of a 2D array along "x", i.e: +// arr = arr[:, ::-1] +// launched with grid (Nx/2, Ny) +__global__ void reverse2D_x(float* array, int Nx, int Ny) { + uint x = blockDim.x * blockIdx.x + threadIdx.x; + uint y = blockDim.y * blockIdx.y + threadIdx.y; + if ((x >= Nx/2) || (y >= Ny)) return; + uint pos = y*Nx + x; + uint pos2 = y*Nx + (Nx - 1 - x); + float tmp = array[pos]; + array[pos] = array[pos2]; + array[pos2] = tmp; +} + diff --git a/nabu/io/utils.py b/nabu/io/utils.py index 3f2f01afe7591dfe3c8e27493172e32cd8c56407..60163eacac939003bda426452f561a227a4eaa7d 100644 --- a/nabu/io/utils.py +++ b/nabu/io/utils.py @@ -105,6 +105,11 @@ def get_first_hdf5_entry(fname): return entry +def hdf5_entry_exists(fname, entry): + with HDF5File(fname, "r") as fid: + res = fid.get(entry, None) is not None + return res + def get_h5_value(fname, h5_path, default_ret=None): with HDF5File(fname, "r") as fid: try: diff --git a/nabu/preproc/double_flatfield_cuda.py b/nabu/preproc/double_flatfield_cuda.py index d8b06062048fdbcc74c2bd8d5eb71a6878b47ab8..4d08eca4037b66ec48e7367ceb8b481a6c358eee 100644 --- a/nabu/preproc/double_flatfield_cuda.py +++ b/nabu/preproc/double_flatfield_cuda.py @@ -63,8 +63,12 @@ class CudaDoubleFlatField(DoubleFlatField, CudaProcessing): return o @staticmethod - def _proc_mlog(x, o): - cumath.log(x, out=o) + def _proc_mlog(x, o, min_clip=None): + if min_clip is not None: + garray.maximum(x, min_clip, output=o) + cumath.log(o, out=o) + else: + cumath.log(x, out=o) o *= -1 return o diff --git a/nabu/preproc/sinogram.py b/nabu/preproc/sinogram.py index d1ff14ac3b64269edae5c9562b76da9ad3dfbae5..a9a6242abd86d09b39e9be55726d4a1447d5817e 100644 --- a/nabu/preproc/sinogram.py +++ b/nabu/preproc/sinogram.py @@ -71,7 +71,15 @@ class SinoProcessing(object): ) self.rot_center = rot_center self._rot_center_int = int(round(self.rot_center)) - + # If CoR is on the left: "flip" the logic + self._halftomo_flip = False + rc = self._rot_center_int + if rc < (self.n_x - 1)/2: + rc = self.n_x - 1 - rc + self._rot_center_int = rc + self.rot_center = self.n_x - self.rot_center + self._halftomo_flip = True + # def _set_halftomo(self, halftomo): self.halftomo = halftomo @@ -123,7 +131,12 @@ class SinoProcessing(object): else: sinos = np.zeros(self.sinos_halftomo_shape, dtype=np.float32) for i in range(n_z): - sinos[i][:] = convert_halftomo(radios[:, i, :], self._rot_center_int) + radio = radios[:, i, :] + if self._halftomo_flip: + radio = radio[:, ::-1] + sinos[i][:] = convert_halftomo(radio, self._rot_center_int) + if self._halftomo_flip: + sinos[i][:] = sinos[i][:, ::-1] return sinos diff --git a/nabu/preproc/sinogram_cuda.py b/nabu/preproc/sinogram_cuda.py index aca5caed78e7a58a99bcd6312b8f767c24c08590..309a0c0557d0cb7787d970fbf4e0cc684448514f 100644 --- a/nabu/preproc/sinogram_cuda.py +++ b/nabu/preproc/sinogram_cuda.py @@ -42,6 +42,20 @@ class CudaSinoProcessing(SinoProcessing, CudaProcessing): d = self.n_x - rc # will have to be adapted for varying axis pos self.halftomo_weights = np.linspace(0, 1, d, endpoint=True, dtype="f") self.d_halftomo_weights = garray.to_gpu(self.halftomo_weights) + if self._halftomo_flip: + self.xflip_kernel = CudaKernel( + "reverse2D_x", + get_cuda_srcfile("ElementOp.cu"), + signature="Pii" + ) + blk = (32, 32, 1) + self._xflip_blksize = blk + self._xflip_gridsize_1 = ( + updiv(self.n_x, blk[0]), + updiv(self.n_angles, blk[1]), + 1 + ) + self._xflip_gridsize_2 = self._halftomo_gridsize # Overwrite parent method @@ -80,15 +94,28 @@ class CudaSinoProcessing(SinoProcessing, CudaProcessing): else: sinos = garray.zeros(self.sinos_halftomo_shape, dtype=np.float32) + # need to use a contiguous 2D, array otherwise kernel does not work + d_radio = radios[:, 0, :].copy() for i in range(n_z): + d_radio[:] = radios[:, i, :].copy() + if self._halftomo_flip: + self.xflip_kernel( + d_radio, n_x, n_a, + grid=self._xflip_gridsize_1, block=self._xflip_blksize + ) self.halftomo_kernel( - radios[:, i, :].copy(), # need to copy this 2D, array otherwise kernel does not work + d_radio, sinos[i], self.d_halftomo_weights, n_a, n_x, rc, grid=self._halftomo_gridsize, block=self._halftomo_blksize ) + if self._halftomo_flip: + self.xflip_kernel( + sinos[i], 2*rc, n_a2, + grid=self._xflip_gridsize_2, block=self._xflip_blksize + ) return sinos @@ -125,7 +152,3 @@ class CudaSinoNormalization(SinoNormalization, CudaProcessing): sino, *self._cuda_kernel_args, **self._cuda_kernel_kwargs ) return sino - - - - diff --git a/nabu/preproc/tests/test_halftomo.py b/nabu/preproc/tests/test_halftomo.py index 627706fa0d9a8f0f462e7ffdcab6496d37dd6acf..2a739ed522e89a7851db81d392120a83565f3f11 100644 --- a/nabu/preproc/tests/test_halftomo.py +++ b/nabu/preproc/tests/test_halftomo.py @@ -5,12 +5,13 @@ import h5py from nabu.testutils import compare_arrays, utilstest from nabu.preproc.sinogram import SinoProcessing, convert_halftomo from nabu.cuda.utils import __has_pycuda__ + if __has_pycuda__: import pycuda.gpuarray as garray from nabu.preproc.sinogram_cuda import CudaSinoProcessing -@pytest.fixture(scope='class') +@pytest.fixture(scope="class") def bootstrap(request): cls = request.cls radio, sino_ref, cor = get_data_h5("halftomo.h5") @@ -26,37 +27,87 @@ def bootstrap(request): def get_data_h5(*dataset_path): dataset_relpath = os.path.join(*dataset_path) dataset_path = utilstest.getfile(dataset_relpath) - with h5py.File(dataset_path, 'r') as hf: + with h5py.File(dataset_path, "r") as hf: radio = hf["entry/radio/results/data"][()] sino = hf["entry/sino/results/data"][()] cor = hf["entry/sino/configuration/configuration/rotation_axis_position"][()] return radio, sino, cor - -@pytest.mark.usefixtures('bootstrap') +@pytest.mark.usefixtures("bootstrap") class TestHalftomo: - def test_halftomo(self): sino_processing = SinoProcessing( - radios_shape=self.radios.shape, - rot_center=self.rot_center, - halftomo=True + radios_shape=self.radios.shape, rot_center=self.rot_center, halftomo=True ) sinos_halftomo = sino_processing.radios_to_sinos(self.radios) - _, err = compare_arrays(sinos_halftomo[0], self.sino_ref, self.tol, return_residual=True) - assert err < self.tol, "Something wrong with SinoProcessing.radios_to_sino, halftomo=True" + _, err = compare_arrays( + sinos_halftomo[0], self.sino_ref, self.tol, return_residual=True + ) + assert ( + err < self.tol + ), "Something wrong with SinoProcessing.radios_to_sino, halftomo=True" - @pytest.mark.skipif(not(__has_pycuda__), reason="Need pycuda for this test") + @pytest.mark.skipif(not (__has_pycuda__), reason="Need pycuda for this test") def test_cuda_halftomo(self): sino_processing = CudaSinoProcessing( - radios_shape=self.radios.shape, - rot_center=self.rot_center, - halftomo=True + radios_shape=self.radios.shape, rot_center=self.rot_center, halftomo=True ) d_radios = garray.to_gpu(self.radios) d_sinos = garray.zeros(sino_processing.sinos_halftomo_shape, "f") sino_processing.radios_to_sinos(d_radios, output=d_sinos) sino_halftomo = d_sinos.get()[0] - _, err = compare_arrays(sino_halftomo, self.sino_ref, self.tol, return_residual=True) - assert err < self.tol, "Something wrong with SinoProcessing.radios_to_sino, halftomo=True" + _, err = compare_arrays( + sino_halftomo, self.sino_ref, self.tol, return_residual=True + ) + assert ( + err < self.tol + ), "Something wrong with SinoProcessing.radios_to_sino, halftomo=True" + + @staticmethod + def _flip_array(arr): + if arr.ndim == 2: + return np.fliplr(arr) + res = np.zeros_like(arr) + for i in range(arr.shape[0]): + res[i] = np.fliplr(arr[i]) + return res + + def test_halftomo_left(self): + na, nz, nx = self.radios.shape + left_cor = nx - 1 - self.rot_center + radios = self._flip_array(self.radios) + sino_processing = SinoProcessing( + radios_shape=radios.shape, rot_center=left_cor, halftomo=True + ) + sinos_halftomo = sino_processing.radios_to_sinos(radios) + _, err = compare_arrays( + sinos_halftomo[0], + self._flip_array(self.sino_ref), + self.tol, + return_residual=True, + ) + assert ( + err < self.tol + ), "Something wrong with SinoProcessing.radios_to_sino, halftomo=True" + + @pytest.mark.skipif(not (__has_pycuda__), reason="Need pycuda for this test") + def test_cuda_halftomo_left(self): + na, nz, nx = self.radios.shape + left_cor = nx - 1 - self.rot_center + radios = self._flip_array(self.radios) + sino_processing = CudaSinoProcessing( + radios_shape=radios.shape, rot_center=left_cor, halftomo=True + ) + d_radios = garray.to_gpu(radios) + d_sinos = garray.zeros(sino_processing.sinos_halftomo_shape, "f") + sino_processing.radios_to_sinos(d_radios, output=d_sinos) + sino_halftomo = d_sinos.get()[0] + _, err = compare_arrays( + sino_halftomo, self._flip_array(self.sino_ref), + self.tol, return_residual=True + ) + assert ( + err < self.tol + ), "Something wrong with SinoProcessing.radios_to_sino, halftomo=True" + diff --git a/nabu/resources/computations.py b/nabu/resources/computations.py index ec7b600ba4e7a2f0a7975c55d8336fcfee65e9cb..79a7ff68ed57db65de527f46c48f6e597506c37b 100644 --- a/nabu/resources/computations.py +++ b/nabu/resources/computations.py @@ -69,7 +69,9 @@ def estimate_required_memory(process_config, chunk_size=None, radios_and_slices= if "rotation_axis_position_halftomo" in rec_config: # Slice has a different shape in half acquisition. # Not sure if start_x, ..., end_y are supported - Nx_rec = 2 * rec_config["rotation_axis_position_halftomo"] + rc = rec_config["rotation_axis_position_halftomo"] + Nx_rec = 2 * rc + Nx_rec = max(2 * rc, 2 * (Nx - rc)) Ny_rec = Nx_rec else: Nx_rec = (rec_config["end_x"] - rec_config["start_x"] + 1) diff --git a/nabu/resources/cor.py b/nabu/resources/cor.py index 76bb3463245a1363aa1644df3f3d0a700ac731f4..238bd77403a213fbeb50538e04e37ffe33f43e15 100644 --- a/nabu/resources/cor.py +++ b/nabu/resources/cor.py @@ -7,7 +7,7 @@ class CORFinder: """ An application-type class for finding the Center Of Rotation (COR). """ - def __init__(self, dataset_info, angles=None, halftomo=False): + def __init__(self, dataset_info, angles=None, halftomo=False, do_flatfield=True): """ Initialize a CORFinder object. @@ -23,6 +23,7 @@ class CORFinder: """ self.halftomo = halftomo self.dataset_info = dataset_info + self.do_flatfield = do_flatfield self.shape = dataset_info._radio_dims_notbinned[::-1] self._get_angles(angles) self._init_radios() @@ -78,6 +79,8 @@ class CORFinder: def _init_flatfield(self): + if not(self.do_flatfield): + return self.flatfield = FlatField( self.radios.shape, flats=self.dataset_info.flats, @@ -89,6 +92,8 @@ class CORFinder: def _apply_flatfield(self): + if not(self.do_flatfield): + return self.flatfield.normalize_radios(self.radios) diff --git a/nabu/resources/dataset_analyzer.py b/nabu/resources/dataset_analyzer.py index 224fcede491443be01dd0241aa0ce1a1b91ffb67..3c449ee03624cef560272c4ebb66659364e820b4 100644 --- a/nabu/resources/dataset_analyzer.py +++ b/nabu/resources/dataset_analyzer.py @@ -206,6 +206,11 @@ class HDF5DatasetAnalyzer(DatasetAnalyzer): processes_file = self.processes_file if processes_file is None: processes_file = os.path.join(self.dataset_scanner.path, "nabu_processes.h5") + # By default, write in processes_file. Not sure it is a good idea + results_file = processes_file + # Don't write in processes_file if flatfield = forced + if self.extra_options["force_flatfield"]: + results_file = None lookup_files = [ DataUrl( file_path=processes_file, @@ -216,7 +221,8 @@ class HDF5DatasetAnalyzer(DatasetAnalyzer): self._get_dataset_hdf5_url(), self.dataset_scanner.image_key, lookup_files=lookup_files, - results_file=processes_file, # overwrite processes file + results_file=results_file, + force_load_existing_results=self.extra_options["force_flatfield"], flats_reduction="median", darks_reduction="mean" ) diff --git a/nabu/resources/dataset_validator.py b/nabu/resources/dataset_validator.py index 22dd5d31d542dfd28d61da35a3737b200dfcae2f..4cc2bf9b7abf8f78963fe644f16b49a712b7ed3c 100644 --- a/nabu/resources/dataset_validator.py +++ b/nabu/resources/dataset_validator.py @@ -50,6 +50,16 @@ class NabuValidator(object): return res + def _get_nx_ny(self): + 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)) + else: + nx, nz = self.dataset_infos.radio_dims + ny = nx + return nx, ny + + def convert_negative_indices(self): """ Convert any negative index to the corresponding positive index. @@ -59,8 +69,7 @@ class NabuValidator(object): if self.nabu_config["reconstruction"]["enable_halftomo"]: if self.dataset_infos.axis_position is None: raise ValueError("Cannot use rotation axis position in the middle of the detector when half tomo is enabled") - cor = int(round(self.dataset_infos.axis_position)) - ny = nx = 2*cor + nx, ny = self._get_nx_ny() what = ( ("reconstruction", "start_x", nx), ("reconstruction", "end_x", nx), @@ -199,8 +208,7 @@ class NabuValidator(object): nx, nz = self.dataset_infos.radio_dims rec_params = self.nabu_config["reconstruction"] if rec_params["enable_halftomo"]: - cor = int(round(self.dataset_infos.axis_position)) - ny = nx = 2*cor + ny, nx = self._get_nx_ny() what = ( ("start_x", "end_x", nx), ("start_y", "end_y", nx), diff --git a/nabu/resources/nxflatfield.py b/nabu/resources/nxflatfield.py index 7a3ba19cd288ca3248c9d75ff8f37fd56244e0d8..c2ab9ecc67a297fb3f387d3705f2e782dea46450 100644 --- a/nabu/resources/nxflatfield.py +++ b/nabu/resources/nxflatfield.py @@ -6,6 +6,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 .logger import LoggerOrPrint val_to_nxkey = { @@ -14,6 +15,23 @@ val_to_nxkey = { "projections": ImageKey.PROJECTION.value, } + +def _extract_entry(data_path): + return data_path.lstrip("/").split("/")[0] + +def replace_h5_entry(data_url, new_entry): + split_path = data_url.data_path().lstrip("/").split("/") + split_path[0] = new_entry + new_data_path = "/".join(split_path) + return DataUrl( + file_path=data_url.file_path(), + data_path=new_data_path, + data_slice=data_url.data_slice(), + scheme=data_url.scheme + ) + + + class NXFlatField: """ A helper class to load flats and darks, or to compute the final ones. @@ -166,6 +184,18 @@ class NXFlatField: def _get_existing_results_url(self, data_url): res = {"flats": {}, "darks": {}} + # If entry is incorrect, take first entry + # TODO find a more elegant mechanism + fname = data_url.file_path() + entry = _extract_entry(data_url.data_path()) + if not(hdf5_entry_exists(fname, entry)): + new_entry = get_first_hdf5_entry(fname) + data_url = replace_h5_entry(data_url, new_entry) + self.logger.warning( + "Entry %s in file %s does not exist. Using the first entry %s" + % (entry, fname, new_entry) + ) + # results_path = os.path.join(data_url.data_path(), "results") # NXProcess with HDF5File(data_url.file_path(), "r") as fid: for what in res.keys(): diff --git a/nabu/resources/processconfig.py b/nabu/resources/processconfig.py index 40309d1619843e3c4c5f7c54eb1a62f320469151..8022f7dbb8213ff17c9ee76528bde595fe1ae7a9 100644 --- a/nabu/resources/processconfig.py +++ b/nabu/resources/processconfig.py @@ -82,6 +82,7 @@ class ProcessConfig: self.corfinder = CORFinder( self.dataset_infos, halftomo=self.nabu_config["reconstruction"]["enable_halftomo"], + do_flatfield=self.nabu_config["preproc"]["flatfield_enabled"] ) cor = self.corfinder.find_cor(search_method=cor) self.dataset_infos.axis_position = cor @@ -220,7 +221,7 @@ class ProcessConfig: if rec_options["enable_halftomo"]: rec_options["angles"] = rec_options["angles"][:rec_options["angles"].size//2] cor_i = int(round(rec_options["rotation_axis_position"])) - # New key + # New keys rec_options["rotation_axis_position_halftomo"] = (2*cor_i-1)/2. # New key rec_options["cor_estimated_auto"] = isinstance(nabu_config["reconstruction"]["rotation_axis_position"], str)