Method for Obtaining an Image Biomarker That Quantifies the Quality of the Trabecular Structure of Bones

Abstract

A method is disclosed for obtaining an image biomarker that quantifies the quality of the trabecular structure of bones. The method includes: retrieving high-resolution CT (Computed Tomography) and/or MRI (Magnetic Resonance Imaging) trabecular images from an image database; pre-processing and post-processing the high-resolution CT and/or MRI trabecular images and obtaining the unique image biomarker “QTS”. Pre-processing may include: calculating the region of interest “ROI”; calculating the bone fraction map; eliminating the partial volume effect; and, binarizing. Post-processing may include: skeletonisation and extraction of morphological and structural characteristics. Lastly, the unique image biomarker “QTS” is defined as: QTS=0.7137*Comp1+0.2863*Comp2.

Claims

1. A method for obtaining an image biomarker that quantifies the quality of the trabecular structure of bones, the method comprising: retrieving high-resolution trabecular images generated by a technique selected from among Computed Tomography (CT), Magnetic Resonance Imaging (MRI) and a combination of CT and MRI, wherein the high-resolution trabecular images are retrieved from a medical image database with a large quantity of content from trabecular regions; pre-processing the high-resolution trabecular images, wherein pre-processing the high-resolution trabecular images comprises: obtaining a region of interest (ROI); calculating a bone fraction map; removing a partial volume effect; and binarizing; post-processing the high-resolution trabecular images, wherein post-processing the high-resolution trabecular images comprises: skeletonisation; and extracting morphological and structural characteristics; and obtaining a unique image biomarker (QTS) based on the following equation:
QTS=0.7137*Comp1+0.2863*Comp2, where: Comp1=BV/TV.sub.1*BV/TV+TbTh.sub.1*TbTh+TbSp.sub.1*TbSp+TbN.sub.1*TbN+D2D.sub.1*D2D+D3D.sub.1*D3D; Comp2=BV/TV.sub.2*BV/TV+TbTh.sub.2*TbTh+TbSp.sub.2*TbSp+TbN.sub.2*TbN+D2D.sub.2*D2D+D3D.sub.2*D3D, where: BV/TV.sub.1=0.255; TbTh.sub.1=−0.023; TbSp.sub.1=−0.277; TbN.sub.1=0.280; D2D.sub.1=0.246; D3D.sub.1=0.089; BV/TV.sub.2=0.331; TbTh.sub.2=0.670; TbSp.sub.2=0.066; TbN.sub.2=0.123; D2D.sub.2=−0.239; D3D.sub.2=−0.292, and, where: BV/TV=(BV/TV−mean(BV/TV)/std.dev(BV/TV); TbTh=(TbTh−mean (TbTh)/std.dev(TbTh); TbSp=(TbSp−mean(TbSp)/std.dev(TbSp); TbN=(TbN−mean(TbN)/std.dev(TbN); D2D=(D2D−mean(D2D)/std.dev (D2D); D3D=(D3D−mean (D3D)/std.dev(D3D), wherein BV/TV is associated with a trabecular volume; TbTh is associated with a mean trabecular thickness; TbSp.sub.1 is associated with a mean trabecular separation; TbN.sub.1 is associated with a trabecular number; D2D is associated with a 2D fractal dimension; D3D is associated with a 3D fractal dimension; and VARIABLE_NAME.sub.n, where n=1, 2, is the respective VARIABLE_NAME value to calculate “Comp1” and “Comp2” respectively.

2. The method for obtaining an image biomarker that quantifies the quality of the trabecular structure of bones, according to claim 1, further comprising: retrieving high-resolution CT trabecular images with the following parameters and corresponding values: kVp=120, Tube Current (mAs)<150, Voxel size in millimetres (mm), <150, and Slice Thickness in millimetres (mm)<0.03.

3. The method for obtaining an image biomarker that quantifies the quality of the trabecular structure of bones, according to claim 1, further comprising: retrieving high-resolution MRI trabecular images by magnetic resonance equipment with a static field strength of at least 3T and multi-element coils with parallel acquisition capability, wherein the magnetic resonance equipment comply with the following parameters and corresponding values: Acquisition Mode=3D, TR (in milliseconds) (1.5/3.0T)=16, TE (in milliseconds) (1.5T/3.0T)=5, Flip angel (°)=25, Voxel size (in millimetres)<(0.3×0.3), Slice thickness (in millimetres)<0.3, Phase Direction=A-P, Reconstruction Matrix=512×512, Bandwidth (Hz/pix)=220, Parallel image factor=2, and Number of averages (NSA, NEX)=3.

4. The method for obtaining an image biomarker that quantifies the quality of the trabecular structure of bones, according to claim 1, wherein the medical images with a large quantity of content from trabecular regions represent a wrist, a femoral head or a vertebra of an individual.

5. The method for obtaining an image biomarker that quantifies the quality of the trabecular structure of bones, according to claim 2, wherein obtaining the ROI comprises: delimiting the ROI, which is obtained by detecting the trabecular regions.

6. The method for obtaining an image biomarker that quantifies the quality of the trabecular structure of bones, according to claim 5, wherein calculating the bone fraction map comprises: eliminating heterogeneities in a scale of intensities of the high-resolution MRI trabecular image that contains the ROI that was delimited, applying a statistically-based local threshold algorithm, which determines an intensity value of marrow in a vicinity of each voxel based on the statistics of the nearest neighbour, adjusting the threshold of the high-resolution MRI trabecular image and enabling the scaling of the intensity of the voxels partially occupied by trabecula, and producing a bone volume fraction map.

7. The method for obtaining an image biomarker that quantifies the quality of the trabecular structure of bones, according to claim 5, wherein removing the partial volume effect comprises: dividing each voxel into sub-voxels, which are assigned an intensity value according to the voxel itself and its neighbouring voxels with the condition that a quantity of intensity must be maintained.

8. The method for obtaining an image biomarker that quantifies the quality of the trabecular structure of bones, according to claim 7, wherein each voxel is divided into equal sub-voxels.

9. The method for obtaining an image biomarker that quantifies the quality of the trabecular structure of bones, according to claim 7, wherein binarization comprises: discretising the high-resolution trabecular images obtained from removing the partial volume effect into binary code in such a way that the voxels representing the trabeculae are represented by “1's” and the voxels representing bone marrow are represented by “0's”; applying a bimodal threshold algorithm on a histogram of a volume; minimising a variance of intra-class intensities; and obtaining a 3D binarized volume made up only of “1's” and “0's”.

10. The method for obtaining an image biomarker that quantifies the quality of the trabecular structure of bones, according to claim 1, wherein skeletonisation comprises: iteratively eroding an object so that voxels of a contour of the object are eliminated without breaking the object.

11. The method for obtaining an image biomarker that quantifies the quality of the trabecular structure of bones, according to claim 9, wherein extracting morphological and structural characteristics comprises: calculating, from the binarized volume, the BV/TV defined as the fraction or percentage between a number of voxels tagged as “1's” and a total number of voxels comprising the binarized volume; calculating the TbTh defined as the average of the thickness of all the trabeculae present in the binarized volume based on: detecting contours of the binarized volume; performing a 2D skeletonization each of the slices; using a distance-transforming method applied to each point of the skeleton based on performing the 2D skeletonization for each of the slices, the minimum distance from each point to a contour being obtained and multiplied by two, in order to represent the thickness of the trabeculae at each point of the skeleton; and adding all of the distances stored at each point of the skeleton and dividing by the number of points contained in the skeleton to obtain the TbTh, calculating TbSp by repeating the previous two steps, but starting from an inverted binarized volume; and calculating the TbN using the following equation: $\begin{array}{cc}\mathrm{TbN}=\frac{BV/TV}{TbTh}& \left(\mathrm{equation}1\right)\end{array}$ calculating the D2D, wherein calculating the D2D comprises: obtaining the contour of the surface to be characterised; applying an initial mesh to the surface to be characterised, where each plot contains contours; reducing size of each side of the plots by half; repeating reduction of a matrix size iteratively until a pixel size is reached; and solving equation:
log(N)=−D.sup.2D.Math.log(λ)+k where the first equation represents a ratio between a number of plots “N” corresponding to a size “λ”, the D2D and a proportionality constant; and calculating the D3D, wherein calculating the D3D comprises: obtaining the 3D volume to be characterised; applying an initial mesh of cubes on the surface to be characterised; reducing the size of the initial mesh of cubes iteratively until the size of the voxel is reached; and solving a second equation:
log(N)=−D.sup.3D.Math.log(λ)+k

Description

BRIEF DESCRIPTION OF THE DRAWINGS

[0052] The patent or application file contains at least one drawing executed in color. Copies of this patent or patent application publication with color drawing(s) will be provided by the Office upon request and payment of the necessary fee.

[0053] FIG. 1 shows a flow chart of the steps of the method for obtaining an image biomarker that quantifies the quality of the trabecular structure of bones of the present disclosure.

[0054] FIG. 2(a) shows, within the step of pre-processing the retrieved high-resolution CT and/or MRI images, on the region of interest “ROI”, which exclusively contains trabeculae (dark) and bone marrow (bright).

[0055] FIG. 2(b) shows, within the step of pre-processing the retrieved high-resolution CT and/or MRI images, the calculation of the bone fraction map after the correction of heterogeneities, where the intensity of the voxel is proportional to the quantity of bone contained in the voxel.

[0056] FIG. 2(c) shows, within the step of pre-processing the retrieved high-resolution CT and/or MRI images, an interpolated image after increasing the resolution to eliminate partial volume effects.

[0057] FIG. 2(d) shows, within the step of pre-processing the retrieved high-resolution CT and/or MRI images, the binary discretisation, where the voxels have been classified into trabeculae (white) and marrow (black).

[0058] FIG. 3 shows, within the step of post-processing the retrieved high-resolution CT and/or MRI trabecular images, the result of “skeletonisation”, the thinner line that preserves the connectivity and shape of an object without breaking it superimposed on the contour of the trabecula.

[0059] FIG. 4 shows the diagram of obtaining the trabecular thickness from a binarized volume.

[0060] FIG. 5 shows the process described for calculating the 2D fractal dimension, where: (a) shows the contours of the trabeculae on the surface to be characterised; (b) shows an initial meshing of the surface; (c) after reducing the size of each plot side by half, six plots that do not contain any contours are found; and, (d) reducing the side of each plot by half again, multiple plots that do not contain contours appear.

DETAILED DESCRIPTION

[0061] The embodiments of the disclosure in relation to the figures are described below.

[0062] FIG. 1 shows the method of the present disclosure, where the different steps required to obtain the QTS (Quality of Trabecular Structure) image biomarker are shown.

[0063] Initially, it involves retrieving images from a medical image database, preferably in DICOM format. Currently, the images captured by medical devices are archived and transmitted through the “PACS” (Picture Archiving and Communication System). PACS is a computerised system for digital archiving of medical images (nuclear medicine, computed tomography, ultrasound, mammography, etc.) and for the transmission of these images to dedicated display stations or between these stations via a computer network. Therefore, the first step of the method 10 is the retrieval of high-resolution CT (Computed Tomography) and MRI (Magnetic Resonance Imaging) images from the medical image database.

[0064] The high-resolution CT images must fulfil several minimum parameters which are shown in the following table (Table 1).

TABLE-US-00003 TABLE 1 kVp 120 Tube Current (mAs) <150 Voxel size [mm] <0.3 × 0.3 Slice thickness [mm] <0.3

[0065] The high-resolution MRI images must fulfil several minimum parameters which are shown in the following table (Table 2). The high-resolution MRI images must have been generated by clinical magnetic resonance equipment with a static field strength greater than or equal to 3T and multi-element coils with parallel acquisition capacity specific to the anatomical region to be covered.

TABLE-US-00004 TABLE 2 Acquisition Mode 3D Phase Direction A-P TR [ms] (1.5T/3.0T) 16 Reconstructed matrix 512x512 TE [ms] (1.5T/3.0T) 5 Bandwidth [Hz/pix] 220 Flip angle (°) 25 Parallel image factor 2 Voxel size [mm] <0.3 × Number of averages 3 0.3 (NSA, NEX) Slice thickness [mm] <0.3 Temporal Resolution —

[0066] Preferably, high-resolution CT and MRI images must contain the largest quantity of trabeculae. The largest quantity of trabeculae is concentrated in the regions indicated in Table 3.

TABLE-US-00005 TABLE 3 Wrist Level of distal metaphysis of the radius Femoral Metaphysis, femoral head and neck Head Vertebra Vertebra body

[0067] Once the high-resolution CT and MRI images have been retrieved 10, pre-processing 20 and post-processing 30, which are explained below, are applied.

[0068] Therefore, the second step of the method of the present disclosure is to pre-process 20 the retrieved high-resolution CT and/or MRI images. This second step of the method in turn comprises four sub-steps which are listed below: [0069] Obtaining 21 the region of interest “ROI” (FIG. 2a); [0070] Calculating 22 the bone fraction map (FIG. 2b); [0071] Removing 23 the partial volume effect (FIG. 2c); [0072] Binarizing 24 (FIG. 2d).

[0073] The first sub-step 21 in which the region of interest “ROI” (FIG. 2a) of the second step, pre-processing 20, is obtained consists of the delimitation of the “region of interest” “ROI”, which must contain only trabecular volume. To calculate the trabecular volume, the contour of the trabeculae must be calculated. To calculate the contour of the trabecula, the internal contour of the cortical bone is traced. The internal contour of the cortical bone is obtained by detecting the areas in which a large variation in intensity is observed for those same voxels in adjacent slices.

[0074] The second sub-step 22 in which the bone fraction map (FIG. 2b) of the second step, pre-processing 20, is calculated consists of the elimination of the heterogeneities in the scale of intensities of the MRI image containing the delimited ROI. Due to the high resolution of the MRI images, slight intensity modulations are produced by small local field inhomogeneities in MRI. This means that voxels belonging to the same class located in different regions of the volume may be classified into different classes because they present different signal intensities, resulting in erroneous detection of the bone trabeculae in the subsequent binarization sub-step. To correct it, a statistically-based local threshold algorithm is applied, which determines the intensity value of the marrow in the vicinity of each voxel based on the statistics of the nearest neighbour. This adjusts the threshold of the image and enables the scaling of the intensity of the voxels partially occupied by trabecula, finally producing a bone volume fraction map. This step is not necessary in the case of CT images.

[0075] The third sub-step 23 in which the partial volume effect (FIG. 2c) is eliminated consists of dividing each voxel into sub-voxels, which are assigned an intensity value according to the voxel itself and its neighbours, under the assumption that the quantity of intensity must be maintained. The partial volume effect is shown as a blurring on the images which is produced as a result of the spatial variation of the sample, which is comparable in magnitude to the resolution used during high-resolution CT and MRI imaging. The use of a resolution on the same order or higher than the structure to be acquired together with a low signal-to-noise ratio due to the high-resolution causes the intensity of each voxel to represent the bone fraction contained in the bone itself, which means that the third sub-step described increases the reconstructed spatial resolution. Preferably, each voxel is divided into 8 sub-voxels, which are assigned an intensity value according to the voxel itself and its neighbours, under the assumption that the quantity of intensity must be maintained. The increased resolution after this sub-step results in a filtered volume and isometric voxels (if the original voxel was already isometric; otherwise the voxels would not be isometric), which indicate the bone fraction that occupies them with a level of intensity. It is important to note that if a voxel is divided into 8 sub-voxels, each of these will be half the size of the original in each of its three spatial dimensions.

[0076] The fourth sub-step 24 of binarization (FIG. 2d) consists of discretising the high-resolution images obtained from the third sub-step into binary code in such a way that the voxels representing the trabeculae are represented by “1's” and the voxels representing the bone marrow are represented by “0's”, That is, it applies a bimodal thresholding method based on minimisation of intra-class variance, to select only those voxels corresponding to bone trabecula. The method is implemented on the 3D volume histogram, minimising the variance of intra-class intensities, which maximises the variance between classes. The result is a logical 3D volume made up of only “1's” and “0's”, representing the trabeculae and bone marrow, respectively. The 3D logical volume composed of only “1's” and “0's” is called a “binarized volume”.

[0077] After pre-processing 20, the post-processing 30 is applied. Therefore, the third step of the method of the present disclosure is to post-process 30 the retrieved high-resolution CT and/or MRI images.

[0078] This third step 30 of the method in turn comprises four sub-steps which are listed below: [0079] Skeletonisation 31; and, [0080] Extraction of morphological and structural characteristics 32.

[0081] The purpose of the first sub-step 31 of “skeletonisation” is to produce the connected internal structure with the smallest thickness of an object, preserving its topology and maintaining the shape of the object as much as possible. Skeletonisation, see FIG. 3, consists of iteratively eroding the object so that the voxels of its contour are eliminated without breaking the object. FIG. 3 shows the skeleton 311, the thinner line that preserves the connectivity and shape of an object without breaking it. FIG. 3 also shows the contour 312 that surrounds the skeleton.

[0082] The second sub-step 32 of “extraction of morphological and structural characteristics” serves to structurally characterise the bone trabecula. The structural characterisation of the bone trabeculae is largely determined by the extraction of the trabecular volume or percentage of bone volume in relation to the total volume under analysis (BV/TV, “Bone Volume to Total Volume”) and other morphological parameters, such as the mean trabecular thickness (TbTh), the mean trabecular separation (TbSp) or trabecular pore size and trabecular number (TbN) or trabecular index. The trabecular volume BV/TV is obtained as the fraction or percentage between the number of voxels tagged as “1's”, in the logical volume obtained in the previous step, and the total number of voxels that make up the logical volume. The mean trabecular thickness “TbTh” represents the average of the thickness of all the trabeculae present in the volume. The process of obtaining the mean trabecular thickness is shown in FIG. 4. To do this, starting with the binarized volume 321, the contours of the binarized volume are detected 322 and a 2D skeletonisation is performed 323 for each of the slices. Subsequently, the distance-transforming method is applied 324 to the entire slice, and this transformation is multiplied by the skeleton, thereby obtaining the minimum distance from each point to the contour 325. Multiplying this minimum distance by two, the trabecular thickness for each point of the skeleton is obtained. Finally, adding up the trabecular thickness obtained for each point of the skeleton and dividing it by the number of voxels that make up the skeleton, the mean trabecular thickness is obtained. The mean trabecular separation TbSp is determined by the pore size, or average thickness of the bone marrow present between trabeculae, which means that the process of obtaining the mean trabecular separation is complementary to the process of obtaining the mean trabecular thickness. In order to obtain it, the binarized volume is inverted and the previous process is applied. The trabecular number TbN gives an idea of the number of trabeculae present in the volume, so, in order to obtain it, the ratio between the volume of trabeculae present in the volume and the mean trabecular thickness is established, by applying equation 1:

[00002]$\begin{array}{cc}\mathrm{TbN}=\frac{\mathrm{BV}/\mathrm{TV}}{\mathrm{TbTh}}& \left(\mathrm{equation}1\right)\end{array}$

[0083] The fractal characterisation of structures makes it possible to quantify the spatial irregularities of the formation of complex structures, presenting a significant correlation with its fracture resistance and the direction of crack propagation, which means that obtaining it is important to determine trabecular bone quality. In the present method, a 2D and a 3D focus have been implemented for the in vivo characterisation of the trabecular structure. For the 2D approach, a conventional square or plot counting algorithm is used which consists firstly, of obtaining the contour of the structure, slice by slice, with each slice then divided into smaller plots in each iteration, collecting the number of plots containing a contour for each iteration in a vector, as shown in FIG. 5. FIG. 5 shows the process described for the calculation of the 2D fractal dimension, first, (FIG. 5a) the contour 3251 of the surface to be characterised is obtained. Subsequently, an initial meshing 3252 of the surface is carried out (FIG. 5b), where each of the plots contains contours. After reducing 3253 the size of each side of the plot by half (FIG. 5c), six plots that do not contain any contours are found. Reducing 3254 the side of each plot by half again (FIG. 5d) multiple plots that do not contain contours appear.

[0084] Lastly, using a least squares adjustment, the value of the 2D fractal dimension 325 is obtained by following equation 2:

log(N)=−D.sup.2D.Math.log(λ)+k  (equation 2)

[0085] Equation 2 represents the ratio between the number of plots (N) corresponding to a size (λ), 2D fractal dimension and a proportionality constant.

[0086] The same analysis methodology for calculating the 2D fractal dimension has been extrapolated to the 3D fractal dimension by means of cube counting in this case. The number of cubes that contain trabeculae are stored in a vector and finally, by means of the least squares adjustment method, the 3D fractal dimension is cleared.

[0087] With the steps described above, the values BV/TV, TbTh, TbSp, TbN, D2D (2D fractal dimension) and D3D (3D fractal dimension) of a CT and/or MRI study of the regions indicated in Table 3 are obtained. In order to obtain an image biomarker that quantifies the quality of the trabecular structure “QTS”, a principal components multivariate analysis is applied to the set of parameters extracted for a set of N studies analysed from the same region (“ROI”). The purpose of which is to reduce the dimensionality of the set of extracted parameters (BV/TV, TbTh, TbSp, TbN, D2D and D3D) to a single parameter that quantifies the quality of the trabecular structure referred to in the present disclosure as “QTS”.

[0088] Therefore, the last step of the method is performed after extracting the previous parameters (BV/TV, TbTh, TbSp, TbN, D2D and D3D) from a sample of N studies.

[0089] To do this, first, the previously extracted parameters are normalised and centred, by subtracting the mean value (expressed mathematically as “mean(X)”) and dividing by the standard deviation (expressed mathematically as “std.dev(X)”) of each parameter.

BV/TV=(BV/TV−mean(BV/TV)/std.dev(BV/TV);

TbTh=(TbTh−mean(TbTh)/std.dev(TbTh);

TbSp=(TbSp−mean(TbSp)/std.dev(TbSp);

TbN=(TbN−mean(TbN)/std.dev(TbN);

D2D=(D2D−mean(D2D)/std.dev (D2D);

D3D=(D3D−mean(D3D)/std.dev(D3D);

[0090] After normalising the parameters, by studying the correlation matrix of the aforementioned parameters, two main components or principal factors are obtained, which have been named “Comp1” and “Comp2”, such that “Comp1” represents the largest proportion of original variability and “Comp2” the maximum possible variability not represented by “Comp1”. Each of these principal factors is expressed as a linear combination of the original parameters with weight VARIABLE_NAME.sub.n, n=1, 2 being for “Comp1” and “Comp2”, respectively, as shown in the following equation.

Comp1=BV/TV.sub.1*BV/TV+TbTh.sub.1*TbTh+TbSp.sub.1*TbSp+TbN.sub.1*TbN+D2D.sub.1*D2D++D3D.sub.1*D3D;

Comp2=BV/TV.sub.2*BV/TV+TbTh.sub.2*TbTh+TbSp.sub.2*TbSp+TbN.sub.2*TbN+D2D.sub.2*D2D++D3D.sub.2*D3D;

with [0091] BV/TV.sub.1=0.255; [0092] TbTh.sub.1=−0.023; [0093] TbSp.sub.1=−0.277; [0094] TbN.sub.1=0.280; [0095] D2D.sub.1=0.246; [0096] D3D.sub.1=0.089; [0097] BV/TV.sub.2=0.331; [0098] TbTh.sub.2=0.670; [0099] TbSp.sub.2=0.066; [0100] TbN.sub.2=0.123; [0101] D2D.sub.2=−0.239; [0102] D3D.sub.2=−0.292;

[0103] The values for each of the coefficients indicated above may vary in significantly different populations, and they may even undergo slight variations as the study sample is enlarged.

[0104] Lastly, the unique image biomarker that expresses the quality of the trabecular bone “QTS” is obtained by the final equation:

QTS=0.7137*Comp1+0.2863*Comp2;