image = multiscaleVesselEnhancement3dImg (inImg3d,inMVEParams)

Detailed Description

Multiscale vessel enhancement using Frangi's approach on 3d images.

The MultiscaleVesselEnhancement3dImg algorithm is the extension of the Multiscale vessel enhancement 2d algotihm to the 3d case. It is an iterative approach used to enhance tube-like features in an image where each scale focuses on a specific tube width range.

The algorithm steps are applied for each voxel, for this reason, we avoid the notation of the coordinates $(\textbf{x})$ in the following. For a given scale $s$ , the algorithm computes at each voxel the Hessian matrix $H$ . Its eigen values $\lambda_i, i \in \left\{1, 2, 3\right\}$ are then extracted and sorted in order to identify :

$\vert \lambda_1 \vert \leq \vert \lambda_2 \vert \leq \vert \lambda_3 \vert$

For an ideal tubular structure, we have:

$\begin{eqnarray*} \vert \lambda_1 \vert & \approx & 0 \\ \vert \lambda_1 \vert & \ll & \vert \lambda_2 \vert \\ \vert \lambda_2 \vert & \approx & \vert \lambda_3 \vert \end{eqnarray*}$

A blobness coefficient $R_B$ and a structureness coefficient $S$ are then calculated from $\lambda_1$ , $\lambda_2$ and $\lambda_3$ :

$\begin{eqnarray*} R_B & = & \frac{\vert \lambda_1 \vert}{\sqrt{ \vert\lambda_2 \lambda_3 \vert}} \\ S & = & \sqrt{\lambda_1^2 + \lambda_2^2 + \lambda_3^2} \end{eqnarray*}$

A third coefficient $R_A$ , referring to the aspect ratio of the two largest second order derivatives, is also calculated:

$R_A = \frac{\vert \lambda_2 \vert}{\vert \lambda_3 \vert}$

This value allows to distinguish plate-like structures from line-like structures ( $R_A \approx 0$ in the latter case).

The vesselness for a scale $s$ , noted $\nu_s$ , is finally computed as a combination of these coefficients. For bright features on dark background, this measure can be expressed as:

$\nu_s = \begin{cases} 0 & \text{if } \lambda_2 > 0 \text{ or } \lambda_3 > 0\\ \left( 1 - \exp \left(-\frac{R_A^2}{2 \alpha^2} \right) \right) \exp \left(-\frac{R_B^2}{2 \beta^2} \right) \left( 1 - \exp \left(-\frac{S^2}{2c^2} \right) \right) & \text{otherwise} \end{cases}$

Where $\alpha$ , $\beta$ and $c$ are sensitivity thresholds for each coefficients. $\alpha$ and $\beta$ are input parameters and $c$ is calculated as:

$c = \frac{1}{2} \max_{\textbf{x} \in \Omega} \left( \Vert H(\textbf{x}) \Vert \right)$

Where $\Omega$ is the image domain and $H(\textbf{x})$ is the Hessian matrix of the image at voxel $\textbf{x}$ .

For dark features on bright background, the condition $\lambda_2 > 0 \text{ or } \lambda_3 > 0$ becomes $\lambda_2 < 0 \text{ or } \lambda_3 < 0$ .

The final result is an image containing the maximum value of $\nu$ along all the scales on each pixel:

$\nu(\textbf{x}) = \max_{s \in Scales}\left( \nu_s(\textbf{x}) \right)$

Where $Scales$ is the collection of scales given as input parameter.

The algorithm parametrization is done using the MVEParams data item, which contains:

the mandatory BlackWhite field, set to true if the user wants to enhance bright vessels on dark background ; it equals false for dark vessels on bright background,
the mandatory Scales field, a collection of Real64 values containing the scale list,
BlobSensitivityThreshold, an optionnal parameter corresponding to the $\beta$ threshold (by default, $\beta = 0.5$ ),
AreaSensitivityThresold, an optionnal parameter corresponding to the $\alpha$ threshold (by default, $\alpha = 0.5$ ).

Moreover, the algorithm allows the user to provide a class image, with type UInt8. If this image is set, the algorithm will fill it by the scale index corresponding to the maximum $\nu_s$ :

$OutOptClassImg(\textbf{x}) = argmax_{s \in Scales}\left( \nu_s(\textbf{x}) \right)$

See Multiscale vessel enhancement 2d for a 2d illustration of multiscale vessel enhancement.

Example of Python code :

Example imports

import PyIPSDK

import PyIPSDK.IPSDKIPLFiltering as filter

Code Example

    # opening of input images
    inImg = PyIPSDK.loadTiffImageFile(inputImgPath)
    
    # Create the class image
    outVesselnessImg = PyIPSDK.createImage(PyIPSDK.eImageBufferType.eIBT_Real32, inImg.getSizeX(), inImg.getSizeY(), inImg.getSizeZ())
    outClassImg = PyIPSDK.createImage(PyIPSDK.eImageBufferType.eIBT_UInt8, inImg.getSizeX(), inImg.getSizeY(), inImg.getSizeZ())
    
    mveParams = PyIPSDK.createMVEParameter([1, 3, 5, 7], True)
    filter.multiscaleVesselEnhancement3dImg(inImg, mveParams, outVesselnessImg, outClassImg)

Example of C++ code :

Example informations

Header file

#include <IPSDKIPL/IPSDKIPLFiltering/Processor/MultiscaleVesselEnhancement3dImg/MultiscaleVesselEnhancement3dImg.h>

Code Example

        // Definition of the algorithm parameters
        MVEParamsPtr pMVEParams = createMVEParameter(vScales, false);
        // Calculation of the vesselness image only
        // ========================================
        // Sample with a generated output image
        // ------------------------------------
        ImagePtr pAutoOutImg = multiscaleVesselEnhancement3dImg(pInImg, pMVEParams);
        // Sample with a provided output image
        // -----------------------------------
        // create output image
        ImageGeometryPtr pOutputImageGeometry = geometry3d(eImageBufferType::eIBT_Real32, sizeX, sizeY, sizeZ);
        boost::shared_ptr<MemoryImage> pOutImg(boost::make_shared<MemoryImage>());
        pOutImg->init(*pOutputImageGeometry);
        // compute vesselness
        multiscaleVesselEnhancement3dImg(pInImg, pMVEParams, pOutImg);
        // Calculation of the vesselness and class images
        // ==============================================
        // create output image
        ImageGeometryPtr pOutputClassImageGeometry = geometry3d(eImageBufferType::eIBT_UInt8, sizeX, sizeY, sizeZ);
        boost::shared_ptr<MemoryImage> pVesselnessImg(boost::make_shared<MemoryImage>());
        boost::shared_ptr<MemoryImage> pClassImg(boost::make_shared<MemoryImage>());
        pVesselnessImg->init(*pOutputImageGeometry);
        pClassImg->init(*pOutputClassImageGeometry);
        // compute vesselness and the class image
        multiscaleVesselEnhancement3dImg(pInImg, pMVEParams, pVesselnessImg, pClassImg);