Bone-tissue decomposition of a single X-ray image via solving a Laplace equation

X-ray images

Dense structures like bones block more X-rays, making them appear darker on the sensor compared to other areas. On the other hand, less dense and thinner soft tissue allows more X-rays to pass through, resulting in brighter areas on the sensor.

To improve the visibility of bone regions, modern X-ray images often use a subtraction process. This involves subtracting a constant maximum value, caused by the X-ray dose, from the original image. This operation makes bone regions appear brighter and soft tissue regions appear darker, making the bone region more distinguishable.

In clinical applications, the window technique is commonly used to limit the intensity of bone within a specific range. This effectively removes soft tissue regions that do not overlap with bones. However, this method cannot remove soft tissue that overlaps with bone regions.

To achieve better visualization of bones, other methods aim to enhance image contrast through techniques like histogram equalization, Contrast Limited Adaptive Histogram Equalization (CLAHE), and more (Gong et al., 2019; Aldoury et al., 2023; Yan et al., 2023). However, these methods often enhance both bone and soft tissue simultaneously, resulting in visually improved images but introducing a complex or unknown relationship between image intensity and the actual X-ray dose. Such a complex or unknown relationship causes artifacts and makes the results more difficult to interpret.

To address these challenges, this article proposes a method to separate an X-ray image into a soft tissue image and a bone image. This decomposition technique maintains the linear relationship between image intensity and the physical properties of the objects being imaged. Figure 1 showcases two examples of our method, demonstrating the guaranteed bone enhancement. Further details about our method will be explained in subsequent sections.

Scattered light in physics

As shown in the left column of Fig. 1, bones are usually surrounded by soft tissue, which is similar to many natural scenes. One example is foggy weather, as depicted in Fig. 2A, where the fog can be considered as “soft tissue” (low density) and the buildings as “bone” (high density).

The phenomenon behind this is called light scattering (Van de Hulst, 1958), which occurs in different situations, including X-ray images in clinics, foggy weather in natural scenes, and fluorescence images in biological imaging (Gong & Sbalzarini, 2016; Gong, 2015). Scattered light can affect the quality of an image, such as when soft tissue scatters X-rays, making bone details less clear.

The study of light scattering dates back to 1871 when Rayleigh examined this phenomenon for light wavelengths larger than the particles’ radius in the medium (Van de Hulst, 1958). Mie scattering, which involves spherical particles, was later studied by Mie and has been named after him. These fundamental principles of physics form the basis for modern image dehazing techniques.

Scattered light in natural images

A typical example of scattering light in nature is the fog. As a result, the image from foggy natural scene is not clear. Image processing algorithms that virtually remove the fog are called dehazing (as shown in the top row of Fig. 2).

For natural images, the dehazing mathematical model is Fattal (2008), He, Sun & Tang (2011)

(1) $$f(x,y)=J(x,y)t(x,y)+A(1-t(x,y)),$$where $f$ is the observed image, J is the unknown clear image to be estimated, $t$ is the transmission map to be estimated, and A is the global atmospheric light to be estimated.

In recent years, there has been significant progress in improving the visibility of hazy images. These methods can be divided into three categories: simple contrast enhancement methods (Tan, 2008; Fattal, 2008), dark channel based methods (He, Sun & Tang, 2011) and deep learning methods (Cai et al., 2016).

Initially, the focus was on enhancing the contrast of foggy images, assuming that they lacked sufficient contrast (Tan, 2008; Fattal, 2008). However, these methods often require extensive computation and may result in noticeable artifacts.

A significant breakthrough came with the introduction of the dark channel prior in dehazing methods (He, Sun & Tang, 2011). The dark channel prior assumes the existence of a region with low-intensity values in the local neighborhood. This prior can be efficiently solved using the guided image filter, which has popularized the use of the dark channel prior in dehazing.

Deep learning is another approach for reducing scattering light in images (Cai et al., 2016; Engin, Genc & Ekenel, 2018). It assumes the availability of paired clear and foggy images, and a neural network can be trained to learn the mapping from foggy images to their corresponding clear images. However, in practical applications, acquiring paired images can often be challenging.

Bone suppression or enhancement

Instead of interested by bones, some applications focus on the soft tissue such as pneumonia. For these applications, they try to reduce the visualization of bones. Such task is called bone suppression (Suzuki et al., 2006; Chen & Suzuki, 2014; Li et al., 2020).

In such bone suppression, previous approaches assume that the observed image $f$ is linearly composed by soft tissue $f_{\mathrm{t}\mathrm{i}\mathrm{s}\mathrm{s}\mathrm{u}\mathrm{e}}$ and $f_{\mathrm{b}\mathrm{o}\mathrm{n}\mathrm{e}}$ as Von Berg et al. (2016)

(2) $$f=f_{\mathrm{t}\mathrm{i}\mathrm{s}\mathrm{s}\mathrm{u}\mathrm{e}}+f_{\mathrm{b}\mathrm{o}\mathrm{n}\mathrm{e}}.$$

Surely, such linear composition is fundamentally different from the nonlinear model in Eq. (1).

Moreover, such methods require strong geometric prior information about the imaging objects, such as the rib shapes for chest X-ray images. And they usually require to exactly find the bone boundaries (bone segmentation). Due to such geometric prior and accurate segmenation requirements, these methods are difficult to be extended from one imaging object to another. For example, the methods developed for ribs can not easily be used for feet or knees.

Even with the accurate bone segmentation, the resulting soft tissue images may have obvious artifacts because their linear composition assumption, Eq. (2), is not always valid. In contrast, the nonlinear model Eq. (1) has been shown effective in removing the haze in natural images.

These limitations motivate us to develop a new and generic mathematical model. Instead of suppressing bones or soft tissue, our model decompose one X-ray image into one soft tissue image and one bone image. These two images have exactly the same imaging domain. Therefore, our task is fundamentally different from bone enhancement task and bone suppression task.

Motivation and contributions

The soft tissue in human body usually scatters the X-ray, severely reducing the quality of bone details in the resulting images. This fact motivates us to construct a novel mathematical model that can decompose bones and the soft tissue in X-ray images. The decomposed soft tissue image can be used for its related study such as pneumonia. The bone image can be adopted for its related research such as bone fracture and bone age estimation.

The scattering light in X-ray images by the soft tissue shares the same physical law as the fog in natural images (Gong & Sbalzarini, 2016; Gong, 2015). Thus, the dehazing methods that have been developed for natural foggy images must be also valid on X-ray images (Gong et al., 2019). We transform the dehazing model into a new form for X-ray images with proper assumptions.

Different from the bone enhancement or suppression, we propose to decompose the input X-ray image into one bone image and one soft tissue image. Such task is named as bone and tissue decomposition (BTD). We construct a new mathematical model that can effectively decompose a single X-ray image into one soft tissues image and one bone image. Be aware the difference between our model and the bone segmentation task. Bone segmentation separates the imaging domain into bone region and background region (without overlap). However, our background and bone images share the same imaging domain (exactly overlapped with the same imaging domain). Such difference is illustrated in Fig. 3.

Our contributions are in the following folds:

• We propose a new image processing task named BTD.

• We propose a new mathematical model for BTD. This model is based on the well known image dehazing model, but with proper improvements for X-ray images.

• With proper assumptions, the BTD model leads to a standard Laplace equation, which can be efficiently solved.

• The resulting bone image is theoretically guaranteed to have larger image contrast than the original image.

Portions of this text were previously published as part of a preprint (https://arxiv.org/abs/2007.14510).

Mathematical equation

In this article, we propose a new mathematical model for X-ray image decomposition. Our model is mathematically rooted in the scatted light physics and the dehazing model for natural images in Eq. (1). More specifically, we propose following model for X-ray images:

(3) $$f(x,y)=\frac{1}{\alpha }B(x,y)(1-T(x,y))+T(x,y),$$where $f(x,y)\in [0,1]$ is the observed image, B is the unknown bone image, $T(x,y)$ is the soft tissue image, $\alpha \ge 0$ is a scalar parameter.

When $f(x,y)=T(x,y)$, it would force $B(x,y)=0$. It means that the observation only comes from the soft tissue and there is no bones in the image. When $f(x,y)=\frac{1}{\alpha }B(x,y)$, it would force $T(x,y)=0$, which indicates that the observation only comes from bones. Otherwise, the observation is composed by the tissue and bones images as the similar way in scattered light.

Parameter $\text{α}$ analysis

We use $\alpha$ as global constant variable, instead of spatially varying $\alpha (x,y)$. Although $\alpha (x,y)$ might achieve better visual result, it could introduce artifacts and it would lose the relationship between actual dose and image intensity in X-ray image. But when we use spatially constant $\alpha$, such linear scaling will keep such relationship between the actual physics and the intensity in X-ray images.

In later section, we will prove that $\alpha \ge 1$ (Eq. (11)), which theoretically guarantees to increase the image contrast. This property becomes clear when we set the background $T(x,y)=0$. That is $f(x,y)=\frac{1}{\alpha }B(x,y)$. It means $\mathrm{\nabla }B(x,y)=\alpha \mathrm{\nabla }f(x,y)$, where $\mathrm{\nabla }$ is the standard gradient operator. Since $\alpha \ge 1$, the contrast in bone image $B(x,y)$ is theoretically larger than the contrast in the input image $f(x,y)$. This theoretical property is numerically confirmed by all our experiments.

Relationship with the dehazing model

Our model is different from Eq. (1) with two important changes. First, we define the soft tissue image (also called background image in this article) as

(4) $$T(x,y)=A(1-t(x,y)),\mathrm{w}\mathrm{h}\mathrm{e}\mathrm{r}\mathrm{e}\mathrm{A}=1.$$

Here, we assume $A=1$. The reason is that only X-ray can reach the sensors (there is no other light resource). Second, we define the unknown bone image $B(x,y)$ as a linear scaling of $J(x,y)$

(5) $$B(x,y)=\frac{1}{\alpha }J(x,y),$$where $\alpha \ge 0$ is a scalar parameter. Such linear scaling keeps the physical meaning of the image intensity in X-ray images.

Comparison with the linear model

Our model is nonlinear and inspired by scattered light models in physics and dehazing models for natural images. Once input X-ray image and its bone image are given, the soft tissue image can be computed as

(6) $$T(x,y)=\frac{f-\frac{B}{\alpha }}{1-\frac{B}{\alpha }},$$which is fundamentally different from the linear model, Eq. (2), in bone suppression (Von Berg et al., 2016)

(7) $$f_{\mathrm{t}\mathrm{i}\mathrm{s}\mathrm{s}\mathrm{u}\mathrm{e}}=f-f_{\mathrm{b}\mathrm{o}\mathrm{n}\mathrm{e}}.$$

To quantitatively demonstrate this differential performance, we implement both the nonlinear formulation (Eq. (6)) and its linear counterpart (Eq. (7)) using identical input radiographs $f$ and standardized bone images ( $f_{\mathrm{b}\mathrm{o}\mathrm{n}\mathrm{e}}=\frac{B}{\alpha }$). The experimental validation employs a publicly accessible benchmark dataset (https://www.kaggle.com/raddar/digitally-reconstructed-radiographs-drr-bones) comprising 193 paired X-ray projections and their corresponding osseous segmentation maps. Subsequent soft tissue reconstructions derived from Eqs. (6) and (7) exhibit marked qualitative differences, as exemplified in Fig. 4. Critical evaluation of the resultant images, as annotated by directional indicators (red/green arrows), reveals that the nonlinear paradigm achieves significantly superior bone suppression efficacy relative to the linear approximation model.

Since the linear model and our model take the same input X-ray image and the same bone image, the difference in the estimated soft tissue images can only come from the models themselves. Therefore, the results in Fig. 4 numerically confirm that our model is better than the linear model Eq. (2).

Assumptions

Since our model, Eq. (3), is ill-posed, we have to make some assumptions to solve this model. We make the following four assumptions to simplify the solving process.

• First of all, we assume $T(x,y)\le f(x,y)$. This assumption makes sure that $U(x,y)\ge 0$.

• Second, we assume $0\le T(x,y)<1$, which avoids the denominator to be zero. As a result, $\frac{1}{1\text{-}T(x,y)}>1$, which helps in improving the bone image contrast as shown in later sections (Eq. (11)).

• Third, we assume that $T(x,y)$ is second order differentiable (Gong & Goksel, 2019; Gong & Sbalzarini, 2017). Such smoothness assumption is reasonable because the physical configuration of soft tissue is smooth.

• Fourth, we assume that the maximum value in $B(x,y)$ is one. Such assumption is used to determine the value of $\alpha$. In later section, based on this assumption, we can prove that $\alpha \ge 1$. Such theoretical result guarantees the bone image contrast enhancement as shown in later section.

A flexible bone mask

Now, we have enough assumptions to find $T(x,y)$. We estimate the $T(x,y)$ by a two-step strategy. First, we roughly estimate a mask $M(x,y)$ that covers bones. Be aware that the mask $M(x,y)$ only needs to cover the bones. It does not necessarily be exactly aligned with bone boundaries as the image segmentation task.

Therefore, there are several ways to obtain such mask. First, it can be easily obtained by a simple threshold method followed by morphology operations. Second, it can also be estimated by active contour methods. Third, it can even be given interactively by users’ input. In short, the way of obtaining this mask is flexible.

Moreover, the mask $M(x,y)$ itself is also flexible. We use three different masks for the same input image, as shown in each row of Fig. 5. Although the masks are varying, the resulting soft tissue images and bone images are similar. Be aware that our mask only needs to cover the bones. Its boundary is not necessarily aligned with the bone boundary (the image segmentation task).

Soft tissue image

With a bone mask $M(x,y)$, we can find the soft tissue image $T(x,y)$ by solving a Laplace equation. Let $M(x,y)$ denote our mask. Now, we need to estimate the soft tissue intensity in this mask. This problem can be modeled as following minimization task

(8) $$min\int {\int }_M||\mathrm{\nabla }T|{|}^2\mathrm{d}x\mathrm{d}y,\mathrm{s}.\mathrm{t}.T_{\partial M}=f_{\partial M},$$where $\mathrm{\partial }$ denotes the boundary. The optimal solution of this energy is the standard Laplace equation

(9) $$\mathrm{\Delta }{T}_M=0,\mathrm{s}.\mathrm{t}.T_{\mathrm{\partial }M}=f_{\mathrm{\partial }M}.$$

There are several efficient Poisson solvers available for this equation. We summarize these solvers in Table 1, where each algorithm’s computational complexity is also shown.

We use the convolution pyramid method (Farbman, Fattal & Lischinski, 2011) to solve the this equation for two reasons. First, it is a direct method and so we do not have to iterate. Second, the pyramid method has linear computational complexity.

Thanks to the efficiency of Pyramid method, the solution of Eq. (9) can be easily obtained. The estimated soft tissue image is shown in Fig. 5C. The running time is 0.1 s in MATLAB on a ThinkPad P1 laptop with Intel Xeon E2176 CPU with 2.70 GHz. The image resolution is $1\text{,}022\times 757$. Such performance (7.7 Mpixel/second) is fast enough for clinical applications in practice.

Bone image

After estimating $T(x,y)$, we need to estimate $\alpha$ for bone image $B(x,y)$ estimation. As mentioned, we assume the maximum value in $B(x,y)$ is one. Therefore, we define

(10) $$\alpha \equiv \frac{1}{\underset{(x,y)}{max}\left\{\frac{f(x,y)-T(x,y)}{1-T(x,y)}\right\}}.$$

Since $0\le f(x,y)\le 1$, we have $\frac{f(x,y)-T(x,y)}{1\text{-}T(x,y)}\le 1$. As a result, we can prove

(11) $$\alpha \ge 1.$$

This parameter linearly increases the contrast in the bone image. As mentioned, such linearity can keep the physical meaning of intensity in X-ray images.

Finally, the bone image can be computed as (shown in Fig. 5D)

(12) $$B(x,y)=\alpha \frac{f(x,y)-T(x,y)}{1-T(x,y)}.$$

One example is shown in the first row of Fig. 5. And the middle line intensity profiles of original and our results are shown in Fig. 6. In this case, $\alpha =1.44$ and the image contrast is enhanced. As shown in Eq. (11), $\alpha \ge 1$ and the enhancement is theoretically guaranteed. All our experiments also confirm this property.

Guaranteed contrast enhancement

With our model, we can prove the following theorem

Proof. To show the image contrast enhancement, we need to study the gradient of the resulting bone image. From Eq. (12), we can get the gradient of B

(13) $$\mathrm{\nabla }B=\alpha [\frac{\mathrm{\nabla }f}{1-T}-\frac{1-f}{{(1-T)}^2}\mathrm{\nabla }T].$$

Since we assume $T(x,y)$ is smooth, we know that $\mathrm{\nabla }T(x,y)\approx 0$ for most locations (Gong & Sbalzarini, 2016) (also see the gradient statistics in Fig. 5 from Gong & Goksel (2019)). Therefore, we have

(14) $$\mathrm{\nabla }B\approx \alpha \frac{\mathrm{\nabla }f}{1-T}\ge \alpha \mathrm{\nabla }f\ge \mathrm{\nabla }f.$$

The last inequality comes from the $\alpha \ge 1$ in Eq. (11). This result indicates that the bone image has better image contrast than the input image for most of pixels.

Such theoretical guarantee is important for practical applications, where the robustness of the method is concerned. Our method makes sure that the bone image is clearer than the original input. Moreover, the larger $\alpha$, the better bone image contrast. The $\alpha$ depends on the X-ray source and the imaging objects.

The complete algorithm

In summary, our model in Eq. (3) can be efficiently solved by Algorithm 1.

Different imaging objects

Our model is not restricted by any imaging object. It can work on various bone X-ray images. Several results from our method are shown in Fig. 7, including knees, arm, hand, etc. The left column is the original input image. The right two columns are the resulting soft tissue and bone image from our method, respectively. It can be told that the soft tissue image is smooth as we assumed. Meanwhile, the bone image has better image contrast as desired.

Moreover, our method can reach real-time performance on these X-ray images. The running time of our method on these images is reported in Table 2. Our method is fast enough for these low resolution images. For high resolution images, the running time will be shown in later section. The experiments are performed in MATLAB 2022b (The MathWorks, Natick, MA, USA) on a laptop with Intel Xeon E2176 CPU with 2.70 GHz.

Comparison with dehazing and enhancement

We further compare our method with a classical image enhancement method and a dehazing method for natural images (He, Sun & Tang, 2011), which uses dark channel prior. We tested on ten images and two of them are shown in Fig. 8. The classical image enhancement method (histogram equalization) enhances both the soft tissue and bones. Such nonlinear process losses the relationship between image intensity and physical X-ray dose.

Our model is also different from the dehazing model. The dehazing method for natural images can not completely remove the soft tissue in X-ray image, as shown by the red arrows in Fig. 8. In contrast, our method does not have this issue. This is because we estimate a better soft tissue image. Moreover, our bone image has better image contrast, which is theoretically guaranteed as described.

Comparison with bone suppression and bone enhancement

We compare our soft image with a bone suppression method (Gozes & Greenspan, 2020) and compare our bone image with a bone enhancement method (Dougherty & Lotufo, 2003). We tested them on ten images and three results are shown in Fig. 9, where the left column is the original input image.

For the soft tissue image, our method obtains a result which does not have obvious bones on various X-ray images, including knees, humerus, and hands. In contrast, the method from Gozes & Greenspan (2020) does not suppress the bone much.

For the bone image, our method gives an image that only contains bones. In contrast, the bone enhancement method (Dougherty & Lotufo, 2003) enhances the bones but still keeps some soft tissue in the result.

High resolution images

Finally, we applied our method on a hand X-ray image data set (RSNA), which contains more than 10,000 hand X-ray images. The image has high resolution (usually larger than $1\text{,}514\times 2\text{,}044$). These images are collected from clinical applications. Therefore, we can test the performance of our method on practical images, showing its efficiency and effectiveness.

In each panel of Fig. 10, the input image (left) is decomposed into soft tissue (middle) and bone image (right) by our method. Although we only show several images from the data set, the results for the rest images are similar.

The bone images have better image contrast since the parameter $\alpha \ge 1$ is theoretically guaranteed. Since details on bones become clear, such enhancement can benefit the bone diagnosis in practice.

Moreover, the running time of our method on such high resolution images is less than half second in the MATLAB language on a laptop. Therefore, it can be easily deployed in real applications. If higher performance is required, our model can be solved by the parallel Laplace equation solver on a modern graphic process unit (GPU), which usually has thousands of cores.

We believe that such bone and soft tissue decomposition model is important for X-ray images, bone study, soft tissue diagnosis, etc. Despite the nice mathematical properties of the model, it can be very efficiently solved by solving a standard Laplace equation.