Skip to content
GitLab
Commit 5f007fac authored by Pierre Paleo's avatar Pierre Paleo
Browse files

Change title of notebook tutorial

parent 85eca274
Pipeline #24506 passed with stage
in 48 seconds
Hide whitespace changes
Inline Side-by-side
doc/Tutorials/hello_nabu.ipynb
View file @ 5f007fac
%% Cell type:markdown id: tags:
# Nabu reconstruction demo
# Using the Nabu library from Python
This notebook shows how to use the Nabu software for performing a basic reconstruction of a tomography dataset.
The computations are done on a local machine with a GPU and Cuda available.
%% Cell type:markdown id: tags:
## 1 - Load the dataset informations
We must provide `nabu` with the the configuration file (`nabu.conf`), describing the path to the dataset and the processing steps. This is the equivalent of the `.par` file in PyHST2. In this file, no information is given on the detector size, energy, distance, etc: these informations are extracted from the dataset metadata.
%% Cell type:code id: tags:
``` python
from nabu.resources.processconfig import ProcessConfig
```
%% Cell type:code id: tags:
``` python
conf = ProcessConfig("/home/pierre/notebooks/nabu/nabu.conf")
```
%% Cell type:code id: tags:
``` python
# We can easily get information on the processing steps.
nabu_config = conf.nabu_config
print(nabu_config)
# The same can be done with the dataset structure
dataset_infos = conf.dataset_infos
print([getattr(dataset_infos, attr) for attr in ["energy", "distance", "n_angles", "radio_dims"]])
```
%% Output
{'dataset': {'name': 'my_dataset', 'location': '/data/pierre/id19/crayon/5.06_crayon_W150_60_Al2_W0.25_xc1000_', 'file_prefix': 'my_dataset_', 'num_first_radio': 0, 'num_last_radio': 1499, 'flat_file_prefix': 'refHST', 'dark_file_prefix': 'dark', 'binning': 1, 'binning_z': 1, 'projections_subsampling': 1}, 'preproc': {'flatfield_enabled': True, 'ccd_filter_enabled': True, 'ccd_filter_threshold': 0.04, 'clip_value': 10.0, 'take_logarithm': True}, 'phase': {'method': 'paganin', 'save_output': None, 'file_prefix': 'pag_', 'paganin_delta_beta': 100.0, 'paganin_lmicron': None, 'paganin_marge': 50, 'unsharp_coeff': 0.1, 'unsharp_sigma': 1.0, 'paganin_padding_type': 'zeros'}, 'reconstruction': {'angles_file': None, 'rotation_axis_position': None, 'axis_correction_file': None, 'angle_offset': 0.0, 'method': 'fbp', 'fbp_filter_type': 'ramlak', 'padding_type': 'edges', 'enable_halftomo': False, 'start_x': 0, 'end_x': 2047, 'start_y': 0, 'end_y': 2047, 'start_z': 0, 'end_z': 999, 'iterations': 200, 'method_name': 'tv', 'optim_algorithm': 'cp', 'weight_tv': 0.01, 'preconditioning_filter': True, 'positivity_constraint': True}, 'output': {'location': '/tmp/out', 'file_prefix': 'my_dataset', 'file_format': 'hdf5', 'multiframe': False, 'regroup_files': False}, 'resources': {'method': 'slurm', 'nodes': 1, 'gpus_per_node': 2, 'memory_per_node': ((100.0, None), True), 'cores_per_node': ((100.0, None), True), 'walltime': (1, 0, 0)}, 'about': {'nabu_version': '2019.1.0', 'nabu_config_version': '2019.1.0', 'verbosity': 1}}
[63.6, 1.0, 1500, (2048, 1000)]
%% Cell type:markdown id: tags:
## 2 - Chunk processing
Nabu processes data by chunks of radios (see the [documentation](http://www.silx.org/pub/nabu/doc/definitions.html#radios-and-sinograms) for more explanations).
In a first step, we define how to read chunks of radios.
%% Cell type:code id: tags:
``` python
from nabu.io.reader import ChunkReader
```
%% Cell type:markdown id: tags:
What is the largest chunk size we can process ?
The answer is given by inspecting the current GPU memory, and the processing steps.
%% Cell type:code id: tags:
``` python
from nabu.cuda.utils import get_gpu_memory
from nabu.resources.computations import estimate_chunk_size
```
%% Cell type:code id: tags:
``` python
chunk_size = estimate_chunk_size(
get_gpu_memory(0),
nabu_config,
dataset_infos
)
print("Chunk_size = %d" % chunk_size)
```
%% Output
Chunk_size = 300
%% Cell type:code id: tags:
``` python
# Load the first 'chunk_size' lines of all the radios
sub_region = (None, None, None, chunk_size) # start_x, end_x, start_z, end_z
chunk_reader = ChunkReader(dataset_infos.projections, sub_region=sub_region)
```
%% Cell type:code id: tags:
``` python
# Load the current chunk
chunk_reader.load_files() # takes some time
```
%% Cell type:code id: tags:
``` python
print(chunk_reader.files_data.shape)
print(chunk_reader.files_data.dtype)
```
%% Output
(1500, 300, 2048)
uint16
%% Cell type:markdown id: tags:
## 3 - Initialize the GPU
Most of the processing can be done on GPU (or many-core CPU if using OpenCL).
With `pycuda.gpuarray` (or its OpenCL counterpart `pyopencl.array`), we manipulate array objects with memory residing on device. This allows to avoid extraneous host <-> device copies.
%% Cell type:code id: tags:
``` python
import pycuda.gpuarray as garray
from nabu.cuda.utils import get_cuda_context
import numpy as np
```
%% Cell type:code id: tags:
``` python
# Create a Cuda context on device ID 0
# By default, all following GPU processings will be bound on this context
ctx = get_cuda_context(device_id=0)
```
%% Cell type:code id: tags:
``` python
radios = chunk_reader.files_data.astype(np.float32)
n_angles, n_z, n_x = radios.shape
# transfer the chunk on GPU
d_radios = garray.to_gpu(radios)
```
%% Cell type:markdown id: tags:
## 4 - Pre-processing
Pre-processing utilities are available in the `nabu.preproc` module.
Utilities available with the cuda backend are implemented in a module with a `_cuda` suffix.
%% Cell type:markdown id: tags:
### 4.1 - Flat-field
%% Cell type:code id: tags:
``` python
from nabu.preproc.ccd_cuda import CudaFlatField
```
%% Cell type:code id: tags:
``` python
# Configure the `FlatField` processor
cuda_flatfield = CudaFlatField(d_radios, dataset_infos.flats, dataset_infos.darks, sub_region=sub_region)
```
%% Cell type:code id: tags:
``` python
# Perform the normalization on GPU
if nabu_config["preproc"]["flatfield_enabled"]:
print("Doing flat-field")
cuda_flatfield.normalize_radios()
```
%% Output
Doing flat-field
%% Cell type:markdown id: tags:
### 4.2 - Hotspots correction
%% Cell type:code id: tags:
``` python
from nabu.preproc.ccd_cuda import CudaCCDCorrection
```
%% Cell type:code id: tags:
``` python
# This processing is not done in-place, so another array has to be created
cuda_ccd = CudaCCDCorrection(d_radios)
d_radios_corr = garray.zeros_like(d_radios)
```
%% Cell type:code id: tags:
``` python
if nabu_config["preproc"]["ccd_filter_enabled"]:
print("Doing hotspots correction")
cuda_ccd.median_clip_correction(output=d_radios_corr)
```
%% Output
Doing hotspots correction
%% Cell type:markdown id: tags:
### 4.3 - Phase retrieval
%% Cell type:code id: tags:
``` python
from nabu.preproc.phase_cuda import CudaPaganinPhaseRetrieval
```
%% Cell type:code id: tags:
``` python
# Phase retrieval is done on each radio individually, with the sub-region specified above
if nabu_config["phase"]["method"] != None:
print("Doing phase retrieval")
cudapaganin = CudaPaganinPhaseRetrieval(
(n_z, n_x),
distance=dataset_infos.distance * 1e2,
energy=dataset_infos.energy,
delta_beta=nabu_config["phase"]["paganin_delta_beta"],
pixel_size=dataset_infos.pixel_size
)
for i in range(n_angles):
cudapaganin.apply_filter(d_radios_corr[i], output=d_radios[i])
```
%% Output
Doing phase retrieval
%% Cell type:markdown id: tags:
### 4.4 - (Optional) Unsharp mask
%% Cell type:code id: tags:
``` python
from nabu.misc.unsharp_cuda import CudaUnsharpMask
```
%% Cell type:code id: tags:
``` python
if nabu_config["phase"]["unsharp_coeff"] > 0:
print("Doing unsharp")
cuda_unsharp = CudaUnsharpMask((n_z, n_x), 1.0, 5.)
for i in range(n_angles):
cuda_unsharp.unsharp(d_radios[i], d_radios_corr[i])
d_radios.set(d_radios_corr)
```
%% Output
Doing unsharp
%% Cell type:markdown id: tags:
### 4.5 - Logarithm
%% Cell type:code id: tags:
``` python
from nabu.preproc.ccd_cuda import CudaLog
```
%% Cell type:code id: tags:
``` python
if nabu_config["preproc"]["take_logarithm"]:
print("Taking logarithm")
cuda_log = CudaLog(d_radios, clip_min=0.01)
cuda_log.take_logarithm(output=d_radios)
```
%% Output
Taking logarithm
%% Cell type:markdown id: tags:
## 5 - Reconstruction
We use the `Reconstructor` helper class, which is a front-end for `nabu.reconstruction.fbp`.
%% Cell type:code id: tags:
``` python
from nabu.reconstruction.reconstructor_cuda import CudaReconstructor
```
%% Cell type:code id: tags:
``` python
# Reconstruct all n_z=100 horizontal slices (i.e perpendicular to z axis)
#R = CudaReconstructor(
# d_radios.shape, (0, None), axis="z", vol_type="projections",
# extra_options={"padding_mode": "edges"}
#)
#d_recs = garray.zeros((n_z, n_x, n_x), "f")
# Reconstruct 100 vertical slices perpendicular to the "y" axis
R = CudaReconstructor(
d_radios.shape, (500, 600), axis="y", vol_type="projections",
extra_options={"padding_mode": nabu_config["reconstruction"]["padding_type"]}
)
d_recs = garray.zeros((n_z, 100, 2048), "f")
```
%% Cell type:code id: tags:
``` python
R.reconstruct(d_radios, output=d_recs)
# Finally, retrieve the data from GPU
recs = d_recs.get()
```
%% Cell type:code id: tags:
``` python
```
%% Cell type:markdown id: tags:
## 6 - Visualize
%% Cell type:code id: tags:
``` python
%pylab nbagg
```
%% Output
Populating the interactive namespace from numpy and matplotlib
%% Cell type:code id: tags:
``` python
#imshow(recs[0], cmap="gray")
imshow(recs[:, 50, :], cmap="gray")
```
%% Output
<matplotlib.image.AxesImage at 0x7fe0b2878da0>
%% Cell type:code id: tags:
``` python
```
......
Supports Markdown
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