Diffusion equation quantification: selective enhancement algorithm for bone metastasis lesions in CT images (2025)

Bone metastases (BM) are pathological conditions that induce fracture, paralysis, or pain, and which are highly related with the deterioration of patients' quality of life and overall survival (Ulas et al 2016, Hara et al 2021). Early diagnosis and treatment for BM are therefore crucially important. Lesions of BM (LBM) are mostly observed using x-ray photography, computed tomography (CT), magnetic resonance imaging (MRI), single photon emission CT, positron emission tomography (PET), and bone scintigraphy. Of these, CT is able to scan the whole-body patient state conveniently. Consequently, CT imaging can reveal LBMs of early-stage cancer patients, although the information is not necessarily easy for diagnosticians to identify and use (Rybak and Rosenthal 2001, Yang et al 2011, Łukaszewski et al 2017). Typically, LBMs present complex images with mixed osteoblastic and osteoclastic tissues that entail some difficulty of distinction.

To discern LBMs, diagnosticians must adjust the window level or window width for numerous images, which requires empirical skill or knowledge. Sometimes, no readily apparent signals or findings suggests that LBM exist in the first place. In contrast, machine learning (ML) methods, while incorporating consideration of stochastic characteristics (Shalev-Shwartz and Tewari 2009, Byrd et al 2011, Duchi et al 2011, Ouyang et al 2013, Hariri-Ardebili et al 2022), texture patterns (Hegde et al 2018, Patel et al 2019, Zhang et al 2020, Naeem et al 2021, Jena et al 2022), flows (Heishman and Duric 2007, Fukami et al 2019, Lagemann et al 2021, Peper et al 2022), or mathematical transformations (Livani and Evrenosoglu 2013, Amin et al 2015, Von Lilienfeld et al 2015, Spellings and Glotzer 2018, Martínez-Más et al 2019, Dutordoir et al 2020, Tat Dat et al 2020, Varanasi et al 2020, Anetai et al 2021), can extract more complex feature of images, which are useful for detecting very small lesions. Deep learning (DL) is a powerful method to detect and analyze these hard-to-find and incomprehensible features. To date, several DL methods to find LBMs using CT images have been reported (Roth et al 2015, Chmelik et al 2018, Noguchi et al 2020, 2022, Afnouch et al 2023, Huo et al 2023, Koike et al 2023). Nevertheless, DL methods depend on a dataset that might typically include manually annotated teaching data (Shrestha and Mahmood 2019, Indolia et al 2018, Mikołajczyk and Grochowski 2018). Moreover, the immense amounts of data require enormous effort for annotation (Noguchi et al 2020), or require reinforcement learning to deal adaptively with unknown data which entails computational time costs (Shrestha and Mahmood 2019, Sutton 2018, Szepesvári 2022) despite using graphical processing unit. In addition, inscrutability of the complex networks used for DL presents difficulty when interpreting whether the good results are attributable to a specific network method or dataset (Weber et al 2023). Marcus et al presented ten concerns related to DL (Marcus 2018). Traditional ML methods using mathematical or technical approaches are still very useful for the computer vision (O'Mahony et al 2020). Current DL methods are unable to handle large clinical variations adaptively under rapidly changing conditions (Nisha and Meeral 2021, Ahmad et al 2023).

To provide a technical solution based on principles and general facts, one effective mathematical approach is diffusion-equation (DE) based image processing. Mainly, DE-based image processing executes smoothing (Barash 2002, Fang et al 2008, Li et al 2018), noise reduction (Shih et al 2009, Rafsanjani et al 2016, Liao and Feng 2021), or object segmentation (Geiger and Yuille 1991, Gilboa and Osher 2007, Kim et al 2011) for an image. Generative AI imaging based on diffusion is a state-of-the-art method that is able to fuse various image features naturally (Ho et al 2020, Cao et al 2024, Smolders et al 2024).

Perona and Malik developed an anisotropic diffusion (PMD) model (Perona and Malik 1990). This method has been applied to a wide diversity of images (Lee et al 2006, Maalouf et al 2008, Chen et al 2010, Mirebeau et al 2015) including diagnostic MRI images (Lee et al 2006, Ahmed and Mohamad 2008, Hossain et al 2013, Perkunder and Ivanova 2013). In fact, PMD and its latest customized techniques of optimal diffusion coefficient (Gupta and Lamba 2021), nonlinear kernel method (Maiseli 2023), and combined isotropic diffusion (Wang et al 2018) have shown good performance for complex object images. Nevertheless, although these methods present benefits of natural smoothing over an image retaining the edge of some object, their results exhibit a balancing between denoising and edge enhancing for general images similar to the traditional PMD method. The LBM regions require enhancement of the small image gradient of complex tissues. To overcome these shortcomings, we originally developed a DE quantification (DEQ) model aimed at diffusing a focused region selectively and at enhancing LBM identification (figure 1). Herein, we demonstrate DEQ as an effective technique for improving visibility against LBM. This technique is expected to provide effective assistance of diagnoses and enhancement of computer visualization in terms of contrast while preserving the shape of the image object.

2.1.Image processing based on DE

2.2.Numerical calculation for DE

2.3.Calculation algorithm for DEQ

We numerically calculated DEQ according to equation (10). The calculation algorithm is the following:

Algorithm 1. DEQ calculation.
Data: original CT image matrix for the positive directions of row and column the size of is
Data: 0–1 normalized image , where the time developing is to
Initialization: are set as 1.0
diffusion coefficient
coefficient calculation
coefficient calculation
implicit solving loop is set as 10
for T iterations do:
for 10 iterations do:
for N, M iterations do:
Crank–Nicolson parameter
Crank–Nicolson parameter
end
normalization (0–1)
end
normalization (0–1)
end
normalization (0–1)
return

2.4.Determining the filter function of DEQ for LBM and its evaluation

2.5.Dataset and calculation conditions

We used one CT image dataset (16 bit 512 × 512 pixels, 0.937 mm pixel size, 2 mm slice thickness) of a BM patient who received radiotherapy (Case 00), where the LBM region was diagnosed based on bone scintigraphy. For verification, we also calculated 16 PET CT images (16 bit 341 × 512 pixels, 1.17 mm pixel size, 3.26 mm slice thickness) (Case 01–16) including typical BM of vertebrae, pelvic bone, limbs, ribs, and femoral heads. These CT images followed the Digital Imaging and Communications in Medicine format. The patient datasets were approved by the Institutional Review Board and Independent Ethics Committee of Kansai Medical University (approval number, 2022065). The calculations of DEQ and the PMD models were conducted using the entire image. DEQ, TV, MPM models were calculated using calculation code developed in house using Python 3.9.4 (with numpy 1.23.4, matplotlib 3.6.2, pydicom 2.3.0). For this study, PMD was calculated based on the code of medpy 0.4.0. Computations were conducted by a central processing unit of 11th Gen Core i9-11900H 2.50 GHz (Intel Corp.).

3.1.Comparisons between PMD and DEQ, and determination of optimal DEQ conditions

The PMD model demonstrated nice smoothing with increased iteration. Figure 4 presents a representative CT image that displays an understandable and visible LBM. The PMD performed smoothing and enhanced the boundary between vertebral bone and tissues with its iterating process, whereas the detailed regions, such as between tissues and tissues, were also averaged and became less visible. Especially affected lesions in the bone became less than one-third of the original. This smoothing of similar levels and contrasting of different levels is a well-known effect of PMD or DE imaging, whereas this processing radically exerted the opposite effect against LBM; strong smoothing deteriorates . In the peripheral region, despite nearly 1.3 times improvement of visibility (), boundary expansion occurred because high promoted strong diffusion (a3 of figure 4), which also became the opposite effect of microlesion extraction.

In contrast, DEQ performed various effects with respect to different levels of diffusion. The enhancement of soft tissue diffusion, such as and was radically problematic to the LBM. In contrast, and with respect to , especially in performed the LBM enhancement outstandingly (figure 5). The tendency was more significant in the profile of the LBM region (figure 6). the monotonically increasing types ( and ) performed better in the LBM region compared to the monotonically decreasing types ( and ). Regarding , representing osteoclastic soft tissue region, improved by around 125%–135%. , representing the border regions between bone and soft tissues, were greatly improved to approximately 160% while retaining the original peak location when . An increase of and iteration simply enhanced the diffusion effect. From equation (10), the term can be a negative value when , which accelerated strong smoothing effect of ambient objects related to similar levels of , as shown in figures 4 and 5. Consequently, assisted strong diffusion with an increase in iterations . Of the DEQ calculations, there were calculation cases where deteriorated when high with additional iterations as shown in tables 1–3. However, rather improved . Although the best performer slightly degraded the value, the profile followed the features of the original (figure 6c). The visibility of the soft tissues in the osteoblastic region was apparently rather improved (figure 5(a)). Therefore, for this study, we determined as the optimal condition for the LBM.

Table 1.Visibility criterion for the calculation area of the entire image.

	PMD		Iteration
			Original ()
			0.015 74	0.005 10	0.005 14	0.004 68	0.003 66
	DEQ		Diffusion propagation factor
	Filter function
			0.004 56	0.004 13	0.005 40	0.005 53	0.005 53
			0.004 34	0.005 17	0.005 37	0.005 37	0.005 33
			0.004 24	0.005 13	0.005 17	0.005 00	0.004 90
			0.004 50	0.004 87	0.004 70	0.004 23	0.003 90
			0.005 54	0.004 32	0.004 40	0.003 83	0.003 17
			0.002 76	0.003 37	0.003 18	0.002 94	0.002 72
			0.003 07	0.003 46	0.003 44	0.003 25	0.003 04
			0.003 26	0.003 55	0.003 68	0.003 56	0.003 38
			0.003 13	0.003 71	0.003 83	0.003 67	0.005 14
			0.002 54	0.003 67	0.005 13	0.004 64	0.002 72
			0.004 15	0.005 82	0.006 15	0.006 57	0.006 89
			0.004 20	0.005 60	0.005 75	0.006 07	0.006 41
			0.004 24	0.005 13	0.005 15	0.004 98	0.004 90
			0.005 62	0.004 85	0.004 56	0.004 17	0.003 75
			0.005 83	0.006 35	0.006 55	0.006 50	0.006 39
			0.005 09	0.005 09	0.005 08	0.005 07	0.005 07
			0.004 71	0.004 14	0.003 65	0.003 19	0.002 78
			0.003 26	0.003 55	0.003 68	0.003 56	0.003 38
			0.002 94	0.003 28	0.002 93	0.002 64	0.002 39
			0.002 78	0.003 05	0.002 65	0.002 38	0.002 15

Table 2.Visibility criterion for the calculation area of the entire image.

	PMD		Iteration
			Original ()
			−0.002 10	−0.002 01	−0.001 01	−0.001 19	−0.001 25
	DEQ		Diffusion propagation factor
	Filter function
			−0.001 34	−0.002 36	−0.002 88	−0.003 19	−0.003 91
			−0.001 50	−0.002 55	−0.003 03	−0.003 76	−0.004 00
			−0.001 55	−0.002 54	−0.002 95	−0.003 60	−0.003 74
			−0.001 40	−0.002 24	−0.002 61	−0.002 73	−0.003 06
			−0.002 82	−0.001 57	−0.001 95	−0.002 19	−0.002 31
			−0.001 37	−0.001 84	−0.002 21	−0.002 20	−0.002 08
			−0.001 33	−0.001 87	−0.002 30	−0.002 36	−0.002 31
			−0.001 34	−0.001 94	−0.002 44	−0.002 54	−0.002 53
			−0.001 44	−0.002 05	−0.002 59	−0.002 67	−0.003 06
			−0.001 55	−0.002 30	−0.002 44	−0.002 75	−0.002 56
			−0.002 31	−0.003 99	−0.004 68	−0.005 74	−0.005 66
			−0.00 216	−0.003 24	−0.004 33	−0.004 56	−0.005 33
			−0.001 54	−0.002 54	−0.002 95	−0.003 60	−0.003 74
			−0.002 57	−0.001 40	−0.001 73	−0.002 01	−0.002 23
			−0.002 59	−0.002 88	−0.003 03	−0.003 08	−0.003 11
			−0.002 24	−0.002 24	−0.002 23	−0.002 22	−0.002 21
			−0.001 94	−0.001 19	−0.001 50	−0.001 74	−0.001 93
			−0.001 34	−0.001 94	−0.002 44	−0.002 54	−0.002 53
			−0.001 52	−0.002 16	−0.002 18	−0.002 36	−0.002 11
			−0.001 54	−0.002 08	−0.002 01	−0.002 09	−0.001 79

Table 3.Visibility criterion for the calculation area of the entire image.

	PMD		Iteration
			Original ()
			−0.048 46	−0.066 53	−0.066 88	−0.067 17	−0.067 44
	DEQ		Diffusion propagation factor
	Filter function
			−0.054 73	−0.033 93	−0.022 57	−0.016 70	−0.013 23
			−0.051 97	−0.028 53	−0.018 30	−0.013 33	−0.010 53
			−0.050 80	−0.026 07	−0.016 30	−0.011 67	−0.006 85
			−0.053 23	−0.027 97	−0.017 30	−0.011 93	−0.006 60
			−0.064 00	−0.042 70	−0.029 97	−0.021 93	−0.016 43
			−0.034 11	−0.013 96	−0.007 92	−0.005 37	−0.003 11
			−0.037 89	−0.017 01	−0.009 99	−0.006 85	−0.005 17
			−0.040 56	−0.019 45	−0.011 78	−0.008 20	−0.006 24
			−0.039 56	−0.018 48	−0.011 14	−0.007 78	−0.005 95
			−0.032 92	−0.013 04	−0.007 44	−0.005 12	−0.003 07
			−0.036 38	−0.015 19	−0.007 59	−0.006 57	−0.006 23
			−0.038 09	−0.015 94	−0.010 22	−0.006 35	−0.005 81
			−0.050 82	−0.026 06	−0.016 30	−0.011 66	−0.006 86
			−0.068 28	−0.056 84	−0.045 08	−0.036 63	−0.030 38
			−0.074 50	−0.078 51	−0.078 85	−0.076 61	−0.074 22
			−0.066 42	−0.066 38	−0.066 34	−0.066 30	−0.066 26
			−0.061 91	−0.053 80	−0.047 03	−0.041 41	−0.036 75
			−0.040 56	−0.019 45	−0.011 78	−0.008 20	−0.006 24
			−0.027 92	−0.009 90	−0.005 48	−0.002 91	−0.002 36
			−0.024 63	−0.008 09	−0.003 31	−0.002 42	−0.001 97

3.2.Comparisons between DEQ and PMD models in various clinical cases using PETCT

The DL method is a powerful tool to assist diagnoses of LBMs (Roth et al 2015, Chmelik et al 2018, Noguchi et al 2020, 2022, Afnouch et al 2023, Huo et al 2023, Koike et al 2023). However, DL entails high costs for examining datasets with time- and effort-dependencies for the learning processes. This study provides an effective method in terms of data amounts and time costs to process images to assess LBMs visually using DE. The PMD and its improved models have demonstrated excellent performance to denoise or smooth images using DE. However, the diffusion characteristics of image processing have become obstacles to contrast organ objects in images, particularly between soft tissues as shown in figures 8(a)–(d). In contrast, our method (DEQ) can provide a soft-tissue contrasted image for difficult lesions (figures 5– 8) by selectively diffusing the high pixel value region (figure 9(b)) while maintaining image quality (figure 8(e)). Consequently, DEQ also works in the various thoracic and abdominal parts of CT image to enhance the LBM region over the CT image (appendix figure s2(a)). Surprisingly, this technique is applicable to general 8-bit RGB images. Moreover, this technique outperformed boundary extract of tiny objects in the images, as shown in appendix figure s2(b).

4.1.Filter function

4.2.Calculation strategy of iteration

We adopted step-by-step normalization process for the implicit computation process in the algorithm. Initial condition starts with the , then spatial row and column are calculated over the image according to equation (10). We empirically set 10 iterations using Crank–Nicolson method for implicit computations. However, the normalization process rather makes the descent optimization toward convergence progress slow with respect to cost . Therefore, the calculation result is actually similar to an explicit solution despite using Crank–Nicolson method. You et al demonstrated an infinite number of minima for calculating an image with anisotropic diffusion of PMD (1996). DEQ continuously provides a slightly different solution from iteration to for . Our purpose is not to find the accurate diffusion phenomena but to assess the effects of image processing. Therefore, the time factor, i.e., an iteration, is determined qualitatively. DEQ applied the effect of iterations to the evaluation.

4.3.Diffusion control parameter

4.4.Impacts on diagnosis and clinical usability

4.5.Limitations and prospects for future study

This study demonstrated DEQ as a novel and effective method to enhance the lesions for diagnostic purpose without regard to the dataset size or time costs. Moreover, applying DE to CT images to detect an LBM region is the state-of-the-art for imaging technology. Nevertheless, several limitations are noteworthy. First, this study demonstrated a proof-of-concept in several clinical cases. Although the calculation concept must be effective and prominent, and although the sufficient performance of DEQ was confirmed by our coauthor's five radiation oncologists, further studies including more diverse clinical cases should be conducted beyond a single institution. Second, although the detected enhancements were similar to those of PET fusion, the clinical evaluations for the area should be discussed further. If we deal with CT images with functional information similarly PET fusion, then diagnoses would be more meaningful. Third, some regions showed improved visibility whereas the diffusion characteristic of DEQ processing blurred object-boundary of tissues. Therefore, this result indicates the method was helpful for determining visible lesions and locations, but not for quantification such as the size. Finally, improvement of visibility (-value) or image quality (SSIM) alone would not be sufficient for diagnosis or to prevent human oversight. To be used more in computer vision techniques, some combinations of ML or DL techniques can be explored for automatic diagnosis. Studies conducted by our group show that analyses of radiological reports using natural language processing can facilitate determination of whether an enhanced region is a certain lesion (Doi et al 2023). Moreover, LDIA, which we used to depict the effect of diffusion, is useful to detect object-peripheral enhancement clearly (Anetai et al 2023). Combined with these cutting-edge techniques are expected to bring greater improvement in future studies.

Derived data supporting the findings of this study are available from the corresponding author. Y.A., on request. Calculation source cords regarding DEQ and LDIA are provided in our GitHub repositories, anetaiys/DEQ and anetaiys/LDIA.

The data cannot be made publicly available upon publication because they contain sensitive personal information. The data that support the findings of this study are available upon reasonable request from the authors.

Improved PM diffusion models have been studied extensively. Their validity has been reported for denoising. In this section, three popular calculation models are discussed and compared, i.e., total variation (TV) model (Rudin et al 1992), modified Perona–Malik (MPM) (Wang et al 2013), and coherence-enhancing (CE) diffusion model (Weickert 1999, Xiao et al 2013) (we modified the original model). These methods have similar DE approaches to image processing for denoising. The DE in directional Laplacian form using the diffusion coefficient can be rewritten as

where represents an image, is an arbitrary unit vector, is a Laplacian operator, and represents an image filter function that originally denoted the diffusion coefficient with respect to the image gradient in PM model as follows,

where is a constant that is usually fixed as .

Now, we assume as a normal vector as shown in the following.

where denotes partial derivative operator with respect to Then, Equation (C1) is rewritten as shown below, which is equivalent to the TV model (Wang et al 2013),

The MPM model slightly modifies Equation (C4) with a smoothing effect as

where is Laplacian, and where is a coefficient for the magnitude of the smoothing effect. For this study, was used.

One of another famous approach is coherence of the image (Weickert 1999). Here, the general image structure tensor is given as,

where is convolution operator, and is a convolution filter; usually a Gaussian filter is used. However, strong noise reduction deteriorates mixed tissue regions. Therefore, we used in this study. Then, the eigenvalues and are obtained as

The image coherence is obtained as . The key idea for the diffusion calculation is based on the dependence of the diffusion rate in the image. However, weak coherence areas such as flat areas tend to retain noise (Xiao et al 2013). We applied this coherence logic to the Laplacian term of Equation (C5), which denotes MPM2 model in this study. Therefore, the final form of the diffusion is obtained as

where is the same coefficient for the magnitude of the smoothing effect. Then, we used .

Those models are calculated numerically using explicit method. Therefore, , where is the image processing denoted by the right-hand term of equations (C4), (C5), and (C9). To match the DEQ calculation condition, the 0–1 normalized process was applied to the iteration.

This study used 13-nearest-neighbor points for SSIM calculation, because SSIM for the entire image has insufficient sensitivity to detect image quality. Now, we have th processed image according to a calculation algorithm, then the 13-nearest-neighbor points are obtained as

which we can specifically obtain where and of the -row, -column image array. We ignored two rows and columns at the end of the image array because of 0. Therefore, the 13-nearest-neighbor mean value and standard deviation are obtainable for each pixel. The luminance, contrast, and structure for the images A and B can be represented as

where represents a covariance with respect to the corresponding th and th component of each image. The constant factors , , and are used for this study for the 0–1 normalized image. Then, SSIM is calculated as:

The values of and are obtained as the sum of the components divided by , defining the mean of , , , and , respectively.

LDIA is one of a very useful method to detect the discrepancy using image flow from the original image (Anetai et al 2023). The flow can be defined as a gradient for the image . Let images A (original) and B (calculated) be and , and let their flows be and , where and are tensor suffixes, namely . LDIA calculates the flow deviation from the original flow as, (figure 9(a)), which is expanded as,

Then, we can calculate following values with respect to the deviation field ,

where represents the intensity of deviation caused by a change in flow. denotes vorticity of the deviation field, which is associated with the robustness of deviation vector against ambient change. represents the flow amount, thereby assisting in finding the deviation flow concentration. We calculated the -value map in this study. Specifically, we can calculate the components of the deviation field as

Therein, the -value mean is obtained as the sum of the components divided by .