Regularised iterative reconstruction#

This tutorial covers the main contribution of ToMoBAR - the advanced model-based iterative reconstruction with regularisation. It is the plug-and-play functionality of model-based methods makes them versatile for different types of data with various noise models and artefacts. On the other hand, it can be difficult to find the optimal solution (even when it formally exists) by tediously sweeping through multiple parameters. We suggest here to check various Demos provided, as well as this more compact tutorial.

We start by defining a 3D projection data Numpy array of unsigned integer 16-bit data type (optional) and with axes labels given as ["detY", "angles", "detX"]. We also provide the corresponding flats and darks fields (also 3D arrays of the same axes order).

The first step is to normalise dataRaw using the normaliser function from tomobar.supp.suppTools.normaliser.

from tomobar.supp.suppTools import normaliser

data_norm = normaliser(dataRaw, flats, darks, log=True, method='mean', axis=1)

Instantiate the iterative reconstructor tomobar.methodsIR:

from tomobar.methodsIR import RecToolsIR


detectorVert, angles_number, detectorHoriz = np.shape(data_norm)

Rectools = RecToolsIR(
    DetectorsDimH=detectorHoriz,  # Horizontal detector dimension
    DetectorsDimV=detectorVert,  # Vertical detector dimension
    CenterRotOffset=0.0,  # Center of Rotation (needs to be found)
    AnglesVec=angles_rad,  # A vector of projection angles in radians
    ObjSize=detectorHoriz,  # The reconstructed object dimensions
    datafidelity="LS",  # Data fidelity, choose from LS, KL, PWLS and SWLS
    device_projector="gpu", # Device to perform reconstruction on
    )

Now we have an access to all methods of this particular reconstructor. There are two advanced model-based reconstruction methods available: FISTA and ADMM. The former one is significantly more flexible than the latter, so we recommend it. FISTA has been also accelerated with Ordered-Subsets and requires less memory than ADMM to solve the sub-problem for the data-fidelity term.

Please note that the dictionaries needed for all iterative methods with exact keyword arguments defined in tomobar.supp.dicts.

_data_ = {"projection_norm_data": data_norm,
          "data_axes_labels_order": ["detY", "angles", "detX"],
}  # data dictionary

lc = Rectools.powermethod(_data_) # calculate Lipschitz constant

_algorithm_ = {"iterations": 300, "lipschitz_const": lc} # algorithm dict

_regularisation_ = {
    "method": "PD_TV",  # Selected regularisation method
    "regul_param": 0.000002,  # Regularisation parameter
    "iterations": 150,  # The number of inner regularisation iterations
    "device_regulariser": "gpu",
}

FISTA_Rec = Rectools.FISTA(_data_, _algorithm_, _regularisation_)

Let us consider a faster and slightly more advanced modification of the FISTA algorithm - Penalised Weighted Least Squares (PWLS) Ordered-Subsets FISTA with Total Variation regularisation and WAVELETS thresholding (note that pypwt package needed for wavelets, see Dependencies).

from tomobar.methodsIR import RecToolsIR

detectorVert, angles_number, detectorHoriz = np.shape(data_norm)

Rectools = RecToolsIR(
    DetectorsDimH=detectorHoriz,  # Horizontal detector dimension
    DetectorsDimV=detectorVert,  # Vertical detector dimension
    CenterRotOffset=0.0,  # Center of Rotation (needs to be found)
    AnglesVec=angles_rad,  # A vector of projection angles in radians
    ObjSize=detectorHoriz,  # The reconstructed object dimensions
    datafidelity="PWLS",  # Data fidelity, choose from LS, KL, PWLS and SWLS
    device_projector="gpu", # Device to perform reconstruction on
    )

_data_ = {
    "projection_norm_data": data_norm,  # Normalised projection data
    "projection_raw_data": dataRaw,  # Raw projection data
    "OS_number": 6,  # The number of ordered-subsets
    "data_axes_labels_order": ["detY", "angles", "detX"],
}

lc = Rectools.powermethod(_data_)

_algorithm_ = {"iterations": 25, "lipschitz_const": lc}

_regularisation_ = {
    "method": "PD_TV_WAVELETS",  # Selected regularisation method
    "regul_param": 0.000002,  # Regularisation parameter for PD-TV
    "regul_param2": 0.000002,  # Regularisation parameter for wavelets
    "iterations": 30,  # The number of regularisation iterations
    "device_regulariser": "gpu",
}

FISTA_OS_PWLS_Rec = Rectools.FISTA(_data_, _algorithm_, _regularisation_)

There are hundreds of different data fidelities and regularisation combinations possible in ToMoBAR. Please note, however, that before using a certain combination of data and prior terms, its worth knowing approximately what could be the problem with your data. For instance, you might want to know what is your reconstructed object characteristics (geometry, texture etc.) and if your data contains noise, zingers, stripes, or/and other data inaccuracies? See an example in Synthethic data reconstruction.

One can also operate purely on CuPy arrays if Dependencies are satisfied for the CuPy package. For that one needs to use tomobar.methodsIR_CuPy class instead of tomobar.methodsIR. Note that the array of angles for the CuPy modules should be provided as a Numpy array.