A novel fusion method for X-ray phase contrast imaging based on fast adaptive bidimensional empirical mode decomposition

Abstract

BACKGROUNDS:

X-ray phase contrast imaging (XPCI) can separate the attenuation, refraction, and scattering signals of the object. The application of image fusion enables the concentration of distinctive information into a single image. Some methods have been applied in XPCI field, but wavelet-based decomposition approaches often result in loss of original data.

OBJECTIVE:

To explore the application value of a novel image fusion method for XPCI system and computed tomography (CT) system.

METHODS:

The means of fast adaptive bidimensional empirical mode decomposition (FABEMD) is considered for image decomposition to avoid unnecessary information loss. A parameter δ is proposed to guide the fusion of bidimensional intrinsic mode functions which contain high-frequency information, using a pulse coupled neural network with morphological gradients (MGPCNN). The residual images are fused by the energy attribute fusion strategy. Image preprocessing and enhancement are performed on the result to ensure its quality. The effectiveness of other image fusion methods has been compared, such as discrete wavelet transforms and anisotropic diffusion fusion.

RESULTS:

The δ-guided FABEMD-MGPCNN method achieved either the first or second position in objective evaluation metrics with biological samples, as compared to other image fusion methods. Moreover, comparisons are made with other fusion methods used for XPCI. Finally, the proposed method applied in CT show expected results to retain the feature information.

CONCLUSIONS:

The proposed δ-guided FABEMD-MGPCNN method shows potential feasibility and superiority over traditional and recent image fusion methods for X-ray differential phase contrast imaging and computed tomography systems.

Keywords

X-ray phase contrast imaging X-ray grating interferometry fast adaptive bidimensional empirical mode decomposition pulse coupled neural network image fusion

1 Introduction

Traditional absorption-based X-ray imaging technology has been extensively used in numerous fields in the past few decades such as industrial flaw detection and medical imaging. The generation of the Talbot grating interferometer based on the phase stepping method and the proposal of the Talbot-Lau theory made the X-ray source free from high coherence and long exposure time, which led to the rapid development of X-ray phase contrast imaging (XPCI) and X-ray dark field imaging (XDFI) [1 –5]. In recent years, there have been some developments in the simulation and optimization of Moiré fringes in Talbot grating interferometers [6, 7]. According to the theory of Angular signal radiography, XPCI and XDFI reflect the attenuation, refraction and scattering signals of the object under detection, which can be expressed as attenuation contrast image (AC), differential phase contrast image (DPC), and dark field contrast image (DFC), respectively [8]. Compared with the traditional attenuated X-ray imaging technology which only provide absorption imaging, various characteristics of the object under measurement could be fully indicated due to the complementary detailed information provided by DPC and DFC images. During the recent years, the XPCI has been widely reported as a necessary tool in more fields on account of the advantages proposed above. For instance, in the medical field, XPCI is often introduced to diagnose breast cancer tissue specimens because of its prominent benefits of discriminating trend for cancerous areas with microcalcifications or milder symptoms [9 –11]. In addition, differential phase and dark-field signals of the lungs could be used to identify a variety of types of lung diseases [12 –14]. In the field of industrial detection, new examination possibility in fiber composite materials was offered since that DPC and DFC images are sensitivity to porosity or micro morphology [15 –19].

For the reason that three contrast images represent different physical characteristics, it is not convenient to directly observe in practical applications. Image fusion technology always merges images of different features acquired by a same sample into a single enhanced image. The fusion result always has more detailed information and contains enough characteristics of the images before fusion, making visualization and analysis easier for diagnosis. Image fusion technology has been regularly used in fields which require complementary information to consist a single image, such as the fusion of infrared images and visible light images, the synthesis of multi-focus images, and the combination of images from CT, MRI, SPECT or PET [20 –28].

In the Talbot-Lau grating interferometer, the simultaneous acquisition of signals of different contrasts is particularly convenient for image fusion, since the reconstructed AC, DPC, and DFC images are directly registered in space. Up to now, there has been already lots of means for merging AC and DPC images [29, 30]. Furthermore, Hass et al introduced the concept of centroid from computed graphics subject into the XPCI images fusion and proposed a linearly weighted fusion approach [31]. Zhentian et al proposed a fusion scheme based on principal component analysis (PCA) [32]. Stampanoni et al constructed an image fusion scheme based on a multi-resolution (MR) framework, whose core idea is to convert the AC, DPC, and DFC images into multiple frequency bands [33]. Scholkmann et al introduced a two-step fusion framework based on shift-invariant discrete wavelet transform [34]. Eduardo Coello et al used Fourier domain fusion algorithm, which defined a weighting function to select frequency bands [35]. Gusenbauer et al combined three contrast images after high-pass and low-pass filter with appropriate image masks to obtain fusion result [36, 37].

In the present work, an image fusion framework is constructed to decompose AC, DPC, and DFC images, which based on fast adaptive bidimensional empirical mode decomposition (FABEMD) method. Images will be decomposed into a corresponding set of bidimensional intrinsic mode functions (BIMFs) and the residual images. Subsequently, image preprocessing will be carried out on both the BIMFs and residual images. The BIMFs that contain high-frequency information will be subjected to adaptive Wiener filtering while the lowest-order BIMF component could be discarded optionally to mitigate the influence of noise in the source image. Residual images mainly reflect the distribution of low-frequency energy and can be preprocessed to reduce uneven light intensity distribution in the original image. During the image fusion process, the BIMFs components are fused by using a pulse-coupled neural network with multi-level morphological gradients (MGPCNN), while the residual images are fused using a maximum energy fusion strategy. A new parameter δ is proposed based on the physical meaning of the image obtained by XPCI, which definition is the standard deviation (SD) and average grad (AG) values of the low-order BIMFs signal of the AC and DPC images. The parameter δ reflects the amount of information of the image involved in the fusion and guide the MGPCNN fusion process. Finally, the fused image is acquired by inverse bidimensional empirical mode decomposition (IBEMD) after series of image enhancement means. In order to validate the practicability of this method, various of experimental samples were used to get the results of image fusion. The method proposed in this article is compared with traditional and recent image fusion methods, and further compared with some techniques previously applied for XPCI image fusion. Both subjective image quality judgment and objective image quality evaluation indicators of FABEMD-MGPCNN method show better results. A specialized discussion was also conducted on the limitations of this method. In addition, the proposed method is applied to medical computed tomography and also has an expected effect. This novel image fusion method proposed in this paper is named δ-guided FABEMD-MGPCNN method.

2 Methods

2.1 Principle of X-ray grating interferometer

In this study, AC, DPC, and DFC images of the sample under test were obtained using the phase stepping method with a Talbot-Lau grating interferometer. The instrument makes use of Talbot self-imaging effect and Lau superposition effect [38, 39]. With stepped phase acquisition of the sample and background separately, different structures within the object cause beam diffraction or scattering, resulting in changes in the intensity and shape of the moiré fringes recorded on the detector [40]. Moreover, the contrast changes of the moiré fringes reflect the X-ray absorption characteristics, while the deformation of the moiré fringes is caused by the refraction and small-angle scattering of X-rays passing through the object. According to the theory of Angular signal radiography, the attenuation, refraction and scattering signals of the object under detection can be expressed as attenuation contrast image (AC), differential phase contrast image (DPC), and dark field contrast image (DFC), respectively [8].

The structure of the Talbot-Lau interferometer is illustrated in Fig. 1. The laboratory’s large-focus X-ray tube light source is divided into several sub-light sources via the source grating G0, with each sub-light source capable of generating a Talbot self-imaging pattern behind the G1 grating. By designing specific geometric parameters, the self-imaging patterns of each sub-light source can be separated by one or several periods, as the sub-light sources are incoherent with one another. Additionally, the intensity can be superimposed at the G2 grating to form self-imaging fringes that meet the experimental requirements.

Fig. 1

(a): The structure of the Talbot-Lau X-ray grating interferometer. S denotes the X-ray source, G0 is the source grating, G1 is the phase grating, and G2 is the analyzer grating. O represents the experimental sample to be imaged, while D is the detector. (b): One image with a moiré fringe pattern captured by the detector. (c–e): AC, DPC, DFC images calculated by the moiré fringe images obtained through a five-step phase stepping process.

The formula of the intensity distribution of the moiré fringe patterns for both the sample and background can be expressed as follows: $I_{k}^{r} (x, y) = A_{0}^{r} (x, y) + A_{1}^{r} (x, y) cos [2 π \frac{k}{5} + φ^{r} (x, y)], k = 1, 2, \dots 5,$ (1) $I_{k}^{s} (x, y) = A_{0}^{s} (x, y) + A_{1}^{s} (x, y) cos [2 π \frac{k}{5} + φ^{s} (x, y)], k = 1, 2, \dots 5,$ (2)

Where, $I_{k}^{r}$ and $I_{k}^{s}$ represent the intensity of sample and background image, $A_{0}^{r}$ and $A_{1}^{r}$ represent the zeroth and first-order coefficients after decomposition of $I_{k}^{r}$ , (x, y) represent the horizontal and vertical coordinates, k represents the number of step increments.

The absorption contrast image, dark-field contrast image, and differential phase contrast image are represented as follows:

$Abs (x, y) = \frac{I_{s} (x, y)}{I_{r} (x, y)} = \frac{A_{0}^{s} (x, y)}{A_{0}^{r} (x, y)}$ (3)

$\begin{matrix} Dfc (x, y) \\ = \frac{V_{s} (x, y)}{V_{r} (x, y)} = \frac{A_{1}^{s} (x, y)}{A_{0}^{s} (x, y)} \frac{A_{0}^{r} (x, y)}{A_{1}^{r} (x, y)} = \frac{\sum_{k = 1}^{5} I_{k}^{r} (x, y)}{\sum_{k = 1}^{5} I_{k}^{s} (x, y)} \times \frac{rem \sum_{k = 1}^{5} [I_{k}^{s} (x, y) \times \exp (- 2 π i \frac{k}{5})]}{rem \sum_{k = 1}^{5} [I_{k}^{r} (x, y) \times \exp (- 2 π i \frac{k}{5})]} \end{matrix}$ (4) $\begin{matrix} Dpc (x, y) = \frac{p_{2}}{2 π d} [φ^{s} (x, y) - φ^{r} (x, y)] = \frac{p_{2}}{2 π d} {\arg \sum_{k = 1}^{5} [I_{k}^{s} (x, y) \times exp (- 2 π i \frac{k}{5})] \\ - \arg \sum_{k = 1}^{5} [I_{k}^{s} (x, y) \times exp (- 2 π i \frac{k}{5})]} \end{matrix}$ (5)

Where, Abs (x, y) represents AC image, Dfc (x, y) represents the DFC image, Dpc (x, y) represents the DPC image, φ^s and φ^r represent the image phase of sample and background, p₂ represents the period of the analysis grating, d represents the distance between the phase grating and the analysis grating, arg represents the angle of a complex number and rem represents the modulus of a complex number.

The grating interferometer used to capture the experimental images in the study comprises the following devices: (i) a rotating-anode molybdenum x-ray tube source with a focal spot size of 100μm (Canon, E7290AX) which operated at 39 kV with a mean energy of 17.45 keV, (ii) a CMOS flat-panel detector (Shad-o-Box 6K HS, Teledyne DALSA) with a native element dimension of 49.5×49.5μm², and (iii) three gratings were fabricated by MicroWorks GmbH, Karlsruhe, Germany. The source grating G0, π-phase grating G1 and analyser grating G2 have a period of 6μm with a duty cycle of 0.5, and their total active area is a circular area with a diameter of 70 mm. Moreover, the grating G0 is coated with a 96um±8um gold layer on a silicon substrate and the grating G2 is coated with an 86um±6um gold layer. For more specific information on the equipment used, please refer to the following article [41].

In experimental setup, the Talbot-Lau interferometry is established in a symmetric geometry as Fig. 1. The source grating G0 is located near the beryllium window of the X-ray source. The distance between G0 and G1 is equal to the distance between G1 and G2, which is the third-order Talbot distance (i.e., 38.1 cm). The G2 grating is slightly tilted to generate a moiré pattern with a visibility of 20% along the vertical image direction, and the detector is situated close to the G2 grating. The sample and background are imaged separately to remove the effects of the system’s inherent intensity distribution and noise. For phase stepping method, the complete data is acquired by stepping the G2 grating with five steps, with one second exposure time for each step. The step distance perpendicular to the grating direction is 1.2μm. By processing the acquired moiré fringe images, the AC, DPC, and DFC images of the sample can be obtained.

2.2 Overview of the proposed fusion framework

The image fusion framework proposed in this article is composed of three main parts, as depicted in Fig. 2. (i) Fast adaptive bidimensional empirical mode decomposition is performed for each modal image (AC, DPC, and DFC) of the sample, resulting in BIMFs and a residual image, respectively. A wiener filter is applied to the low-order BIMF component for noise reduction, and the residual image containing vignetting components is corrected for the influence of light distribution. (ii) A two-step fusion method is employed, where the AC image and the DPC image are merged, followed by the fusion of the new image with the DFC image. Specifically, each BIMF of different order is fused by the PCNN with morphological gradients, while the energy attribute (EA) fusion strategy is utilized in residual image fusion. Before the high-frequency information fusion of BIMF components, the SD and PSNR of its low-order components are calculated to obtain parameter δ that reflects the richness of detail information and guides MGPCNN fusion. (iii) Once the fused image’s BIMFs and residual image are obtained, each BIMF is subjected to image enhancement, while the residual image is processed via contrast enhancement and histogram equalization. During the inverse bidimensional empirical mode decomposition process to restore the fused image, the BIMFs and residual image are appropriately linearly weighted to sharpen the image

Fig. 2

The framework proposed in this article for X-ray multimodal image fusion consists of three main parts: decomposition (i), fusion (ii), and image enhancement (iii).

2.3 Fast adaptive bidimensional empirical mode decomposition

Bidimensional empirical mode decomposition (BEMD) is a signal processing method developed by Huang, which is used to extract local characteristics and time-frequency distribution of data, enabling the analysis of nonlinear and non-stationary data [42 –45]. In recent years, new algorithms have been proposed to improve the decomposition speed of each intrinsic mode function (IMF) and the decomposition stopping conditions [46 –51].

The bidimensional empirical mode decomposition (BEMD) method extracts residual image and IMFs in an empirical, adaptive, and data-driven manner [52]. When applied correctly, the decomposition process preserves the original information and avoids the introduction of unnecessary artifacts. BEMD has also been successfully employed in image fusion applications [53, 54]. Bhuiyan et al developed a fast adaptive bidimensional empirical mode decomposition (FABEMD) method, which makes the decomposition of the original data have a faster speed, stronger adaptability, and a better quality of BIMFs.

For a given image I (x, y), FABEMD produces a representation of the sum of BIMFs D_i (x, y) and residual images R_N (x, y): $I (x, y) = \sum_{i = 1}^{N - 1} D_{i} (x, y) + R_{N} (x, y)$ (6)

Let I denotes the original image, D denotes a BIMF, and R denotes the residue. The i-th component of the BIMFs in the decomposition process, denoted as D_i, is obtained from its source image R_i. R_i is a residue image calculated as R_i = R_i-1 - D_i-1, where R₁ = I. The BIMFs are acquired through the following decomposition steps:

Set i = 1. Take the original image I and set R_i = I.

Set j = 1. Set D_{T
_j} = R_i.

Obtain the local maxima map (LMMAX) and the local minima map (LMMIN) of D_{T
_j} denoted as P_j and Q_j, respectively, using the sliding window method with a 3×3 size.

Apply the MAX filter to P_j and the MIN filter to Q_j to create an upper envelope (UE_j) surface and a lower envelope (LE_j) surface.

Smooth UE_j and LE_j using an averaging filter.

Find the mean average envelope (ME) as ME_j = (UE_j + LE_j)/2.

Calculate D_{T
_j+1} as D_{T
_j+1} = D_{T
_j} - ME_j.

Check whether D_{T
_j+1} satisfies the loop stop condition.

$I (x, y) = \sum_{x = 1}^{M} \sum_{y = 1}^{N} \frac{| D_{T_{j + 1}} (x, y) - D_{T_{j}} (x, y) |^{2}}{{| D_{T_{j}} (x, y) |}^{2}}$ (7)

if D_{T
_j+1} meets the criteria given in step (viii), then take D_i = D_{T
_j+1}. Set i = i + 1 first and then calculate R_i = R_i-1 - D_i-1. Go to step (ix). If D_{T
_j+1} does not meet the stopping criteria, then set j = j + 1, go to step (iii) and continue up to step (viii) as before until the criteria are fulfilled [55].

2.4 Preprocessing before image fusion

By adjusting relevant parameters and stopping conditions, different image information can be acquired in BIIMFs. It’s practicable to note that the first-level BIMF reflects most of the noise and contains less detailed image information, while the residual image reflects the light intensity distribution information in the original image. Thus, the quality of the raw images can be improved by such methods. Each first-level BIMF can be either discarded, or alternatively applied by an adaptive Wiener filter to increase the signal-to-noise ratio (SNR). The methods aim to remove the noise as much as possible while preserve edges and high-frequency information in the images. The wiener filter can be expressed as follows, $I_{w} (x, y) = μ + \frac{σ^{2} - n^{2}}{σ^{2}} (I (x, y) - μ)$ (8) where I_w (x, y) means the denoised image after wiener filter, I (x, y) means the input image, μ represents the local mean, σ² represents the local variance and n² represents the noise variance. $μ = \frac{1}{NM} \sum_{x = 1}^{M} \sum_{y = 1}^{N} I (x, y)$ (9) $σ^{2} = \frac{1}{NM} \sum_{x = 1}^{M} \sum_{y = 1}^{N} (I^{2} (x, y) - μ^{2})$ (10)

The parameter M is the total number of rows, and N is the total number of columns in the image.

When the sample is irradiated with an X-ray tube under laboratory conditions, the phenomenon of uneven light intensity distribution often occurs in the resulting image because the light spot can be approximately regarded as a Gaussian distribution. When calculating the AC, DPC, and DFC images, the method of capturing the background significantly reduces the impact of light intensity distribution and inherent noise in the system. During the image capture process, it is unavoidable that the light source exhibits some degree of drift, which can still introduce certain effects on the light intensity distribution. Additionally, random noise may be present in both the sample and background images. In order to further minimize the influence of uneven light intensity distribution when restoring the image, the residual image, which contains few high-frequency information can be discarded.

Figure 3 illustrates the results of using FABEMD to correct image noise and uneven light distribution in a simulation. The image used in this simulation is part of the ‘concordorthophoto.png’ from the MATLAB library.

Fig. 3

The process and results of regulating the images contaminated by noise and vignetting with FABEMD. (a) the original image, (b) the distribution of light, (c) the noisy image, (d) the polluted image by (b) and (c), (e)-(l) the BIMFs and residual component of (d), (m) the improved image, (n) the image optimized by histogram equalization, (o) the histogram distribution of the raw data (a) and the processed data (n).

Figure 3a shows the original image, while Fig. 3b–3d illustrate the distribution of light intensity, noise, and polluted image, respectively. In addition, Fig. 3e–3l show the BIMFs and residual image of Fig. 3d, with seven-level decomposition used for FABEMD. Figure 3m displays the reconstructed image obtained after discarding the residual image and filtering the first-level BIMF, while Fig. 3n shows the result after gamma transformation for the reconstructed image. Figure 3o represents the comparison of the grayscale distribution between Fig. 3n and the original image.

By utilizing the FABEMD method, it is possible to effectively reduce the influence of noise and uneven illumination on the image. After removing the residual component, techniques such as histogram equalization can be employed to enhance the image contrast.

2.5 Principle for image fusion

After decomposing the AC, DPC, and DFC images using FABEMD, the BIMFs of AC and DPC are first fused by PCNN in gradient domain. Then, the EA fusion strategy is utilized in their residual images’ fusion. The same operation is repeated for the fusion of BIMFs and residual images from DFC and the results obtained in the previous step.

2.5.1 Principle for BIMFs fusion

BIMFs mainly contain the edge and detail information of an image. Multi-level morphological gradient (MLMG) is a valid method for conveniently and effectively extracting gradient information from an image. The specific details are expressed as follows [56].

(i) the multi-level structuring element is defined as $Φ_{t} = Φ_{11} \oplus Φ_{12} \oplus \dots \oplus Φ_{1 t}, t \in {1, 2, 3 \dots, N}$ (11) where Φ_1i represents the basic structural unit of the image. The symbol ⊕ represents morphological dilation opretaor, and t represents scale factor.

(ii) the gradient feature G_t can be represented by the operators shows in step (i) from the image I. $G_{t} (x, y) = I (x, y) \oplus Φ_{t} - I (x, y) ⊖ Φ_{t}$ (12) where ⊕ and ⊖ denote the morphological dilation and erosion operators, respectively.

(iii) MLMG can be acquired by weighted sum. $M (x, y) = \sum_{t = 1}^{N} w_{t} \cdot G_{t} (x, y), w_{t} = \frac{1}{t + 1}$ (13)

PCNN offers several advantages, including: (i) No complex training is needed for PCNN [57 –59]. (ii) the spatial feature information of the image is well preserved after iterations due to the high level of connectivity between each neuron and the spatial information. The PCNN model typically comprises three components: the receptive field, the modulation field, and the pulse generator.

The traditional PCNN structure is shown in Fig. 4, while the mathematical structure of PCNN is described by Equation (14) through Equation (21).

Fig. 4

Traditional dual-channel PCNN structure model.

$F_{ij}^{1} (k) = S_{ij}^{1} (k), F_{ij}^{2} (k) = S_{ij}^{2} (k)$ (14)

$V_{O} (k) = {\begin{matrix} 1, & if \sum_{i \in in} Y_{ij} (k - 1) > 0 \\ 0, & otherwise \end{matrix}$ (15) $L_{ij} (k) = e^{- α L} L_{ij} (k - 1) + V_{O} (k)$ (16) $U_{ij}^{1, 2} (k) = F_{ij}^{1, 2} (k) (1 + β \cdot L_{ij} (k) + δ \cdot L_{ij} (k))$ (17) $U_{ij} (k) = Max (U_{ij}^{1} (k), U_{ij}^{2} (k))$ (18) $θ_{ij} (k) = e^{- α T} θ_{ij} (k - 1) + V_{θ} Y_{ij} (k - 1)$ (19) $T_{ij} (k) = T_{ij} (k - 1) + Y_{ij} (k)$ (20) $Y_{ij} (k) = {\begin{matrix} k, & if U_{ij} (k) ⩾ T_{ij} (k + 1) \\ T_{ij} (k - 1), & otherwise \end{matrix}$ (21)

As shown in Fig. 4, $S_{ij}^{1}$ , $S_{ij}^{2}$ in the neural network represent the pixel value of input pre-fusion images after morphological gradients (MG) process at point I₁ (i, j), I₂ (i, j); $F_{ij}^{1}$ , $F_{ij}^{2}$ are the feedback input and L_ij is the link input; β_ij represents link strength; $U_{ij}^{1, 2}$ is the output of modulation field; T_ij is the threshold of step function, θ_ij is the dynamic threshold, and V_θ determines the dynamic threshold amplification factor; α_L and α_F are the time decay constants of the two inputs, respectively. In the practical application process, the high-frequency BIMFs components of the AC image and DPC image after MG process are initially used as inputs I₁ (i, j) and I₂ (i, j). Then, the resulting BIMF components from the fusion process are combined with the BIMFs of the DFC image after MG process as the new I₁ (i, j) and I₂ (i, j) inputs, followed by the second fusion. As shown in the input part of the PCNN model image, $S_{ij}^{1, 2}$ in formula (14) represents the pixel value in I_1,2 (i, j). Formula (15) and (16) are used to calculate L_ij based on preset values. Formula (17) is used to calculate U_ij in the modulation field. Formula (18) gets the max value of U_ij. Formulas (19–21) calculate Y_ij, and determine whether to proceed to the next iteration or output the result based on the preset number of iterations.

In this article, a new parameter δ is proposed to aid in the calculation of the MGPCNN network. The parameter is introduced to calculate physical significance of AC, DPC and DFC images. The parameter δ is defined to guide the iterative calculation of the pulse coupling neural network, and it represents the brightness and detail information of the absorption image (AC), the edge features and gradient information of the phase image (DPC) and dark field image (DFC). $δ = \frac{SD (AC)}{SD (DFC)} \cdot AG (DPC) \cdot AG (DFC)$ (22) where, SD represents the standard deviation of the image, which can be used to characterize the brightness and detail information of the image. AG is the average gradient to acquire the gradient features and edge information of the image.

After the calculation is completed, the δ is substituted into the Y_in in PCNN.

2.5.2 Principle for residual images fusion

The residual images typically contain information about areas with less severe changes and relatively flat gray levels. To effectively use this information, the energy attribute (EA) fusion strategy is employed for each residual image. In the residual image, the low-frequency component contains the primary contour information of the signal, which is equivalent to an approximation of the original signal at a certain scale. The EA fusion strategy is especially suitable for features with slow changes of low frequency. It combines the mean and median values of the image to obtain the relative feature distribution of each point in the image. These distributions are then used as weights in linear weighted fusion. This EA fusion strategy is divided into three steps:

(i) The characteristic attributes of the residual images are computed as $Γ_{A} = (μ_{A}^{2} + M_{A}^{2})^{1 / 2}, Γ_{B} = (μ_{B}^{2} + M_{B}^{2})^{1 / 2}$ (23) where μ and M represent the mean value and the median value of input images I_A and I_B.

(ii) Then, the characteristic attributes are further computed using the following equation: $W_{i} (x, y) = exp (\frac{| I_{i} (x, y) - Γ_{i} |}{ζ}) i = A, B$ (24) where ζ is the modulation parameter.

(iii) Finally, the fused residual image is obtained by taking a weighted mean as follows: $I_{F} (x, y) = \frac{W_{A} (x, y) \times I_{A} (x, y) + W_{B} (x, y) \times I_{B} (x, y)}{W_{A} (x, y) + W_{B} (x, y)}$ (25)

2.6 Enhanced after image fusion

Before performing inverse bidimensional empirical mode decomposition (IBEMD) on the BIMFs and residual image obtained through the previous steps, several operations will be carried out to enhance the quality of the fused image.

To improve the sharpness and increase the detail in the fusion results, a local standard deviation-based method is used on the BIMFS owning to the high-frequency information. Each BIMF component will be divided into contextual regions of size (2m+1) * (2m+1). The pixels in the BIMFs are then enhanced directly by the standard deviation of the context region, which is given by ${BIMF}_{k} (i, j) = μ (i, j) + ξ \times [{BIMF}_{k} (i, j) - μ (i, j)]$ (26) where μ (i, j) is the mean value of contextual region, k is the number of BIMFs, and ξ is a parameter used for enhancement. The value of ξ is determined by the standard deviation of contextual region. $ξ = {\begin{matrix} ξ_{1}, & if 0 < SD (i, j) < {SD}_{h} \\ ξ_{2}, & if {SD}_{h} < SD (i, j) < {SD}_{o} \\ ξ_{3}, & otherwise \end{matrix}$ (27) where SD (i, j) represents the standard deviation of the contextual region. SD_o and SD_h represent one and half of standard deviation of the original image, respectively.

To further sharpen the boundary in the fused image, the Unsharp Masking model described as below is utilized. The model is described as follows: $I_{sharp} = I_{O} + α \times HFilter (I_{O})$ (28) where I_O is the original image, α is the gain parameter, HFilter is a high-pass filter, and I_sharp is the sharpened image.

The residual image mainly contains the background and brightness information of the original image, whose contrast can be enhanced by using histogram modification based function (HMF) and gamma transform. The probability density function (PDF) and the cumulative distribution function (CDF) of the residue can be initially computed as $P (n) = \frac{n_{k}}{N}, n = 1, 2, 3, \dots, L - 1$ (29) $P_{T} (n) = {\begin{matrix} \min {P (n), P (T)}, & if P (n) < T \\ T, & if P (n) ⩾ T \end{matrix}$ (30) $F_{T} (k) = \sum_{n = 1}^{k} P_{T} (n)$ (31) where n represents the intensity level, N represents the total count of pixels, n_k represents the number of pixels having intensity level n, T is a threshold to acquire the P_T (n), while F_T (k) is the CDF.

When reconstructing the fused image using IBEMD, a weighted parameter approach is employed to increase the proportion of high-frequency components, which further enhances the clarity of the image. $I_{fused} = \sum_{k = 1}^{N} ω_{k} \cdot {BIMF}_{k} + residual$ (32) where ω_k can be a scalar (ω_k = 1 when k > 3 and ω_k > 1 when k = 1,2,3.) or a function.

2.7 Fusion image quality evaluation metrics

Except for intuitive visual differences, numerous indicators are used to measure the detail information and the quality of the fusion image in the related papers regarding image fusion [60 –62]. Furthermore, to validate the difference of image fusion techniques for X-ray phase contrast imaging, average grad (AG), standard deviation (SD), and peak signal to noise ratio (PSNR) are considered. SD usually represents the intensity changes before and after image fusion and reflects the amount of information contained in the image, which is defined as [63, 64] $SD = \sqrt{\frac{1}{M \times N} \sum_{i = 1}^{M - 1} \sum_{j = 1}^{N - 1} {(I (i, j) - μ)}^{2}}$ (33)

AG mainly reflects the small and local information, reflecting the clarity of the image. Meanwhile, AG can be utilized to evaluate the contrast and the resolution of fusion image [65]. Its mathematical expressions are as follows, $AG = \frac{1}{(M - 1) (N - 1)} \sum_{i = 1}^{M - 1} \sum_{j = 1}^{N - 1} \sqrt{\frac{{(\frac{\partial I (x, y)}{\partial i})}^{2} + {(\frac{\partial I (x, y)}{\partial j})}^{2}}{2}}$ (34)

PSNR is used to evaluate the spatial reconstruction quality of each band [66, 67]. Higher value of PSNR expresses the better quality of the fusion image.

It could be shown as, $MSE (type) = \frac{1}{M \times N} \sum_{i = 1}^{M} \sum_{j = 1}^{N} (I_{fusion} (i, j) - I_{type} (i, j)) type = AC, DPC or DFC$ (35) $PSNR = 10 \times \log_{10} (\frac{{(max (I_{fusion}))}^{2}}{(\sum MSE) / 3})$ (36)

3 Results

In order to assess the feasibility and effectiveness of δ-guided FABEMD-MGPCNN, this paper provides a comparison of three biological samples (black skirt tetra, tiger barb and Xiphophorus maculatus) with other image fusion methods. The samples were imaged using an X-ray Talbot-Lau grating interferometer system under laboratory condition to obtain their AC, DPC and DFC images. Both subjective and objective evaluations of image quality are demonstrated to compare the results gained by the δ-guided FABEMD-MGPCNN fusion method proposed in this paper with eight traditional or recently introduced image fusion methods. During the FABEMD process, the source image was decomposed into BIMF components 7 times after extensive testing. This separation of information from BIMFs prevents modal aliasing and the influence of ripples, ultimately improving the efficiency of the fusion process. The initial values for the parameters in the morphological gradients PCNN are set as follows, αL = 0.025, αT = 2.5, L₀ = 7, β = 2.8, θ₀ = 1, V_θ = 15, iteration count in MG is 50 and iteration count in PCNN is 150. The specific parameters should be adjusted based on the image quality obtained during the fusion process and the desired retention of detailed features.

The fusion methods used to compare are anisotropic diffusion fusion (ADF), cross bilateral filter fusion (CBF), two-scale image fusion (TIF), stationary wavelet transform fusion (SWT), fourth-order partial differential equation fusion (FPDE), gradient transfer fusion (GTF), multi-scale guided image and video fusion (MGFF) and visual saliency map and weighted least fusion (VSMWLS). The core difference between the above methods are highlighted in Fig. 5.

Fig. 5

The fused results of three samples obtained from various methods. The core difference between the above methods are highlighted. The scale of images is [0 1].

Fish 1, Fish 2, and Fish 3 in Fig. 5. correspond to black skirt tetra, tiger barb, and Xiphophorus maculatus, respectively. The key differences between the fusion results lie in their ability to incorporate the characteristics of the level between the fish bubbles and the bones (Fish1), the structure around the fish eye (Fish2), and the structure around the preliminaries (Fish3). This paper utilizes a combination of objective evaluation metrics and subjective judgment to assess the quality of image fusion results. Several quality evaluation factors were used, namely SD, AG, PSNR and operation time. SD represents the total amount of information while AG reflects the details of the image. PSNR reflects the degree of preservation of image features before and after fusion. The specific calculation results are presented in Table 1. The data highlighted in yellow show higher values, but according to subjective judgment, the corresponding images have significant errors. Therefore, these data should be excluded from the evaluation scope. The red and blue data represent the first and second highest value among the same class. This is consistent with the results of subjective judgment, indicating that the MGFF method, VSMWLS method, and δ-guided FABEMD-MGPCNN method produced better image fusion results. Their fusion images can clearly depict various details present in the original X-ray images while maintaining better contrast and signal-to-noise ratio.

Table 1

Objective evaluation index of different methods

					Fish1
	ADF	CBF	TIF	SWT	FPDE	GTF	MGFF	VSMWLS	Our method
AG	0.018	0.022	0.016	0.021	0.026	0.020	0.016	0.027	0.029
SD	34.9	41.5	42.3	39.1	93.2	49.8	39.7	46.4	55.2
PSNR	13.43	13.89	16.43	15.66	4.575	16.68	16.43	15.88	17.42
time	- -	++	- -	- -	-	+	- -	+	+
					Fish2
	ADF	CBF	TIF	SWT	FPDE	GTF	MGFF	VSMWLS	Our method
AG	0.016	0.021	0.019	0.023	0.022	0.021	0.027	0.030	0.031
SD	38.6	40.5	46.7	41.2	69.7	43.2	56.0	52.1	55.6
PSNR	16.56	14.79	16.68	16.45	20.20	16.27	17.78	17.24	18.47
time	- -	++	- -	- -	-	+	- -	+	+
					Fish3
	ADF	CBF	TIF	SWT	FPDE	GTF	MGFF	VSMWLS	Our method
AG	0.016	0.022	0.016	0.022	0.018	0.015	0.025	0.027	0.026
SD	42.8	40.5	43.9	43.6	91.8	40.7	55.5	53.3	60.8
PSNR	17.14	15.69	17.41	17.45	20.21	16.86	18.33	17.19	18.75
time	- -	++	- -	- -	-	+	- -	+	+

* From - - to ++ in time term means that the time taken by each algorithm is from short to long.

The calculation results demonstrate the strengths of the MGFF method, VSMWLS method and the δ-guided FABEMD-MGPCNN method proposed in this article. In order to further reflect the advantages of the δ-guided FABEMD-MGPCNN method in the fusion process of AC, DPC, and DFC images obtained by XPCI technology, the three methods are implemented on bee and steleophaga plancyi samples.

Figures 6 and 7 display the AC, DPC and DFC images of the bee and steleophaga plancyi samples obtained through XPCI, along with the results of three fusion methods used. In the marked areas of the bee sample, the FABEMD-MGPCNN method produces a result with better contrast and clarity. The PSNR values for the three methods in Fig. 6 are 14.0432 (MGFF, Fig. 6(d)), 14.4235 (VSMWLS, Fig. 6(e)) and 15.0046 (δ-guided FABEMD-MGPCNN, Fig. 6(f)). More detail information can be seen on the legs of the bee sample, such as fluff. In the back and tail area of the bee sample, the MGFF method mainly reflects the brightness value of the AC image, while the detail information and clarity of the VSMWLS method’s result are not as clear as the method proposed in this article.

Fig. 6

(a) AC image, (b) DPC image, (c) DFC image, (d) MGFF method, (e)VSMWLS method, (f) δ-guided FABEMD-MGPCNN method. The scale of image is [0 1].

Fig. 7

(a) AC image, (b) DPC image, (c) DFC image, (d) MGFF method, (e)VSMWLS method, (f) δ-guided FABEMD-MGPCNN method. The scale of image is [0 1].

In the marked area of the steleophaga plancyi sample, the fusion result obtained by FABEMD-MGPCNN exhibits a better signal-to-noise ratio, more detailed information, high contrast and better structural characteristics. The PSNR values for the three methods in Fig. 7 are 12.4913 (MGFF, Fig. 7(d)), 13.0682 (VSMWLS, Fig. 7(e)) and 14.1062 (δ-guided FABEMD-MGPCNN, Fig. 7(f)). In the head and tail area of the sample, the regional characteristics signal of the results acquired by the MGFF method and the VSMWLS method is weak, and the local gradient is not obvious, making it difficult to observe these structures easily.

To further test the applicability of the algorithm proposed in X-ray phase contrast imaging, three contrast images of the bee sample were fused by using three methods cited in previous article. These methods include the linear weighted fusion method proposed by Hass [30], the fusion method based on the shift invariant wavelet discrete transform (SIDWT) framework proposed by Scholkmann [34], and the image fusion method in the Fourier domain proposed by Coello [35].

As shown in Fig. 8, the image fusion result obtained by the δ-guided FABEMD-MGPCNN method preserves more detailed features from differential phase contrast and dark field contrast. The enlarged area of the image with an arrow marking reveals that the proposed method provides clearer texture information and more detailed features.

Fig. 8

The results obtained by four fusion methods used in XPCI, while main difference has been marked in the figure. The enlarged areas in the right image are Fourier-method, Linear weighting, SIDWT and FABEMD-MGPCNN from top to bottom.

In summary, various image fusion methods were compared with the method proposed in this article using three fish samples. The fusion results obtained by different methods were evaluated using subjective and objective evaluation indicators. Among them, the MGFF method, the VSMWLS method and the δ-guided FABEMD-MGPCNN method achieved better results. Further experiments were conducted using bee and steleophaga plancyi samples, which demonstrated that the image fusion results obtained by the FABEMD-MGPCNN method have richer detail information and structural characteristics, while also exhibiting better contrast and signal-to-noise ratio.

4 Discussion

In Section 3, the δ-guided FABEMD-MGPCNN method was found to be effective on the fusion of AC, DPC, and DFC images in 2D imaging XPCI technology. Compared to recent image fusion methods applied in the XPCI field, the proposed method in this paper has also demonstrated advantages in terms of metrics and visual quality. Furthermore, the image fusion method proposed in this paper also yields good results for x-ray 3D reconstruction. Currently, traditional X-ray absorption imaging has been applied in computed tomography (CT) in clinical medicine. However, CT imaging based on phase contrast and dark field contrast has not yet been widely employed. Although the application of XPCI technology in CT can be realized under laboratory conditions, it has not fully entered clinical surgery. To explore the potential application of this method, fish samples were utilized to perform the preliminary attempts of X-ray absorption, phase and dark field CT, and a three-dimensional view of the reconstruction results was obtained, as shown in Fig. 9. Unlike the application of this method in differential phase contrast images discussed earlier, the phase contrast in the 3D reconstruction results obtained by CT is the integral phase. The results demonstrate that the proposed method is also applicable for integrating phase contrast and can achieve good fusion images in XPCI-CT. The different contrast images reflect various characteristics of fish sample marked in the figure. For example, AC images reflect the smooth tissue and the bones of the fish. PC images reflect the teeth and muscle information with high contrast, while DFC images reflect the edge feature information inside the fish. The δ-guided FABEMD-MGPCNN method is used to characterize the feature fusion of the CT reconstruction results. Figure 10 presents that the fusion images involve the characteristic information of the original images such as teeth, bones and muscles as marked in the chart. The expected fusion of CT images highlights the potential application prospects of XPCI technology and the δ-guided FABEMD-MGPCNN image fusion approach in the field of medicine. These methods can aid doctors in acquiring additional pathological characteristics of patients while focusing on these features in a fusion image, thus providing convenience and possibilities for auxiliary diagnosis.

Fig. 9

The three-dimensional view of AC-CT, PC-CT, DFC-CT reconstructions and FABEMD-MGPCNN fusion results of fish sample, while main characteristics have been marked in the figure.

Fig. 10

The impact of changing the number of decomposition layers, the number of iterations in each layer, and the selection of window size on FABEMD method. (i) 1st layer BIMF, (ii) final layer BIMF, (iii) residual image. (a) layers = 5, iterations = 3, type = low size, (b) layers = 3, iterations = 3, type = low size, (c) layers = 7, iterations = 3, type = low size, (d) layers = 5, iterations = 1, type = low size, (e) layers = 5, iterations = 10, type = low size, (f) layers = 5, iterations = 3, type = high size.

Nevertheless, there are some limitations when the δ-guided FABEMD-MGPCNN method used in practical application. For different fusion targets, it is necessary to choose the number of decomposition layers and corresponding iteration times during the mode decomposition process. Figure 10 shows the impact of selecting different parameters on the image decomposition process. In Fig. 10(i–iii) represent the first layer of BIMFs, the last layer of BIMFs, and the residual image. The difference between (a) and (c) lies in the number of decomposition layers, which are 5, 3, and 7 layers, respectively. A higher number of decomposition layers results in closer spatial frequency components of different images in the fusion process. Meanwhile, when there is uneven illumination in the original image, increasing the number of decomposition layers can separate the low-frequency illumination distribution into residual image. However, increasing the number of decomposition layers also greatly raises the computational time. The number of iterations for each layer of (a), (d), and (e) is 3, 1, and 10, with the remaining indicators consistent. In the process of BEMD method, the objective function should be set in the decomposition of each layer in principle. However, the time required is too long in the actual operation process and the number of iterations for each layer is directly set. When the number of iterations is set too high, the decomposed image will produce ripples shown in (e). The window size of the upper and lower pole planes used for reconstructing the image in (a) is the minimum, while (f) is the maximum. Choosing the maximum value as the window size will retain more detailed information in the residual image, which is not suitable for the process of image fusion.

5 Conclusion

This paper proposes a novel image fusion method based on FABEMD to address the requirements of X-ray phase contrast imaging. FABEMD, a data-driven function, is used to extract image features from decomposed BIMFs components. Furthermore, it effectively reduces noise, mitigates the impact of uneven illumination, and enhances image detail quality. The AC, DPC and DFC images are decomposed using FABEMD and image features are extracted from BIMFs and residual components. High-frequency and low-frequency information are fused using the MGPCNN method and the EA strategy, respectively. Comparative evaluations against other fusion methods demonstrate that the proposed approach captures more image features and preserves finer details during the fusion process. In conclusion, the δ-guided FABEMD-MGPCNN exhibits significant improvements over traditional fusion methods in both subjective and objective evaluation criteria. However, there are still areas where further enhancements are needed, such as manually controlling the number of decompositions for different images to achieve optimal quality and requiring a longer period of time in decomposition process.

Funding

This work was funded by the National Natural Science Foundation of China (62205287, 62125504, 61827825 and 61735017), Major Program of the Natural Science Foundation of Zhejiang Province (LD21F050002); Key Research and Development Program of Zhejiang Province (2020C01116), Zhejiang Lab (2020MC0AE01), Zhejiang Provincial Ten Thousand Plan for Young Top Talents (2020R52001).

References

Koyama

, Yoshikawa

, Momose

Phase-contrast X-ray imaging with a triple-bragg-case interferometer, Jpn J Appl Phys Part 2 42(1A-B) (2003), 80–82. doi: 10.1143/JJAP.42.L80.

Momose

, Yashiro

, Takeda

et al. Phase tomography by X-ray Talbot interferometry for biological imaging, Jpn J Appl Phys Part 1 45(1) (2006), 5254–5262. doi: 10.1143/JJAP.45.5254.

Pfeiffer

, Weitkamp

, Bunk

et al. Phase retrieval and differential phase-contrast imaging with low-brilliance X-ray sources, Nat Phys 2(4) (2006), 258–261. doi: 10.1143/JJAP.45.5254.

Pfeiffer

, Bech

, Bunk

et al. Hard-X-ray dark-field imaging using a grating interferometer, Nat Ter 7(2) (2008), 134–137. doi: 10.1038/nmat2096.

Aoki

, Hoshino

, Watanabe

Application of a Zernike-type phase contrast hard X-ray microscope to microtomography, J X-Ray Sci Technol 15(2) (2007), 65–72. doi: 10.1111/j.1365-2818.2006.01592.x.

Tao

S.W.

, Xu

Y.S.

, Bai

et al. Moiré artifacts reduction in Talbot-Lau X-ray phase contrast imaging using a three-step iterative approach, Opt Express 30(20) (2022), 35096–35111. doi: 10.1364/OE.466277.

Yan

A.M.

, Wu

X.Z.

, Liu

Predicting fringe visibility in dual-phase grating interferometry with polychromatic X-ray sources, J X-Ray Sci Technol 28(6) (2020), 1055–1067. doi: 10.3233/XST-200726.

P.Y.

, Zhang

, Bao

et al. Angular signal radiography, Opt Express 24(6) (2016), 5829–5845. doi: 10.1364/OE.24.005829.

Gureyev

T.E.

, Paganin

D.M.

, Arhatari

et al. Dark-field signal extraction in propagation-based phase-contrast imaging, Phys Med Biol 65(21) (2020), doi: 10.1088/1361-6560/abac9d.

10.

Emons

, Fasching

P.A.

, Wunderle

et al. Assessment of the additional clinical potential of X-ray dark-field imaging for breast cancer in a preclinical setup, Ther Adv Med 12 (2020), doi: 10.1177/1758835920957932.

11.

Rauch

, Rieger

, Pelzer

et al. Discrimination analysis of breast calcifications using x-ray dark-field radiography, Med Phys 47(4) (2020), 1813–1826. doi: 10.1002/mp.14043.

12.

Ludwig

, Seifert

, Hauke

et al. Exploration of different x-ray Talbot-Lau setups for dark-field lung imaging examined in a porcine lung, Phys Med Biol 64(6) (2020), doi: 10.1088/1361-6560/ab051c.

13.

Umkehrer

, Morrone

, Dinkel

et al. A proof-of principal study using phase-contrast imaging for the detection of large airway pathologies after lung transplantation, Sci Rep 10(1) (2020), doi: 10.1038/s41598-020-75185-4.

14.

Zhang

, Fowler

A.M.

, Wilke

L.G.

et al. Fast acquisition with seamless stage translation (FASST) for a trimodal x-ray breast imaging system, Med Phys 47 (2020), 4356–4362. doi: 10.1002/mp.14297.

15.

Matsunaga

, Yano

, Endrizzi

et al. Detection of individual sub-pixel features in edge-illumination x-ray phase contrast imaging by means of the dark-field channel, J Phys D: Appl Phys 53(9) (2020), doi: 10.1088/1361-6463/ab5aa0.

16.

Wang

Z.L.

, Shi

X.M.

, Ren

et al. Transmission, refraction and dark-field retrieval in hard X-ray grating interferometry, J Synchrotron Radiat 27 (2020), 494–502. doi: 10.1107/S1600577519017223.

17.

Glinz

, Sleichrt

, Kytyr

et al. Phase-contrast and dark-field imaging for the inspection of resin-rich areas and fiber orientation in non-crimp vacuum infusion carbon-fiber-reinforced polymers, (16), J Mater Sci 56 (2021), 9712–9727. doi: 10.1007/s10853-021-05907-0.

18.

Morimoto

, Kimura

, Shirai

et al. Talbot-Lau interferometry-based x-ray imaging system with retractable and rotatable gratings for nondestructive testing, Rev Sci Instrum 91(2) (2020), doi: 10.1063/1.5131306.

19.

Ludwig

, Seifert

, Niepold

et al. Non-destructive testing of archaeological findings by grating-based X-ray phase-contrast and dark-field, Imaging J 4(4) (2018), doi: 10.3390/jimaging4040058.

20.

Kaur

, Koundal

, Kadyan

Image fusion techniques: A survey arch, Comput Methods Eng 28 (2021), 4425–4447. doi: 10.1007/s11831-021-09540-7.

21.

Chen

Y.L.

, Cheng

L.L.

, Wu

et al. Infrared and visible image fusion based on iterative differential thermal information filter, Opt Lasers Eng 148 (2022), doi: 10.1016/j.optlaseng.2021.106776.

22.

, Kang

X.D.

, Duan

P.H.

et al. Attribute filter based infrared and visible image fusion, Inf Fusion 75 (2021), 41–54. doi: 10.1016/j.inffus.2021.04.005.

23.

Tan

, Zhang

T.P.

, Zhao

L.C.

et al. Multi-focus image fusion with Geometrical Sparse Representation, Signal Processing Image Commun 92 (2021), 116–130. doi: 10.1016/j.image.2020.116130.

24.

Duan

, Zhang

T.P.

, Luo

X.L.

et al. DCKN: Multi-focus image fusion via dynamic convolutional kernel network, Signal Processing 189 (2021), doi: 10.1016/j.sigpro.2021.108282.

25.

Verschure

D.O.

, Verberne

H.J.

Gated SPECT MPI and CT venography fusion: A new approach for appropriate CRT-pacemaker lead placement, J Nucl Cardiol 28 (2021), 1446–1448. doi: 10.1007/s12350-019-01882-x.

26.

Faragallah

O.S.

, El-Hoseny

, El-Shafai

et al. A comprehensive survey analysis for present solutions of medical image fusion and future directions, IEEE Access 9 (2021), 11358–11371. doi: 10.1109/ACCESS.2020.3048315.

27.

Yang

, Park

D.S.

, Huang

S.Y.

et al. Semantic-guided polarization image fusion method based on a dual-discriminator GAN, 18(2) (2010), 157–170. doi: 10.1364/OE.472214.

28.

Shen

, Wen

Y.M.

, Cai

Y.F.

Efficient x-ray image enhancement algorithm using image fusion, J X-Ray Sci Technol 17(3) (2009), 207–220. doi: 10.3233/XST-2009-0223.

29.

Roessl

, Koehler

, Stevendaal

et al. Image fusion algorithm for differential phase contrast imaging, Proceedings of SPIE 8313(54) (2012), doi: 10.1117/12.911255.

30.

Maack

, Roessl

, Koehler

et al. Frequency dependent combination of X-ray images of different modalities, Patient (2016).

31.

Haas

, Polyanskaya

, Bayer

et al. Image fusion in x-ray differential phase-contrast imaging,, Proceedings of Spie (2012), 208–220. doi: 10.1117/12.911223.

32.

Jiang

X.L.

, Zhang

, Wang

Z.T.

et al. Multi-scale image fusion for X-ray grating-based mammography, Proceedings of SPIE 8506 (2012), doi: 10.1117/12.929625.

33.

Wang

, Clavijo

C.A.

, Stampanoni

et al. Image fusion scheme for differential phase contrast mammography, JINST 8 (2013), doi: 10.1088/1748-0221/8/07/C07011.

34.

Scholkmann

, Revol

, Kaufmann

et al. A new method for fusion, denoising and enhancement of x-ray images retrieved from Talbot-Lau grating interferometry, Phys Med Biol 59(6) (2014), 1425–1440. doi: 10.1088/0031-9155/59/6/1425.

35.

Coello

, Sperl

J.I.

, Beque

et al. Fourier domain image fusion for differential X-ray phase-contrast breast imaging, Eur J Radiol 89 (2017), 27–32. doi: 10.1016/j.ejrad.2017.01.019.

36.

Gusenbauer

, Reiter

, Kastner

et al. Detection of non-metallic inclusions in quenched and tempered steel bars by XCT and after fatigue life testing steel res, Int 87 (2016), 386–393. doi: 10.1002/srin.201500083.

37.

Gusenbauer

, Reiter

, Plank

et al. Porosity determination of carbon and glass fibre reinforced polymers using phase-contrast imaging, J Nondestr Eval 38(1) (2019), doi: 10.1007/s10921-018-0529-6.

38.

David

, Nohammer

, Solak

H.H.

et al. Differential x-ray phase contrast imaging using a shearing interferometer, Appl Phys Lett 81 (2002), 3287–3289. doi: 10.1063/1.1516611.

39.

Weitkamp

, Diaz

, David

et al. X-ray phase imaging with a grating interferometer, Opt Express 13(16) (2005), 6296–6304. doi: 10.1364/OPEX.13.006296.

40.

Momose

Recent advances in x-ray phase imaging, (9A), Japan J App Phys 44 (2005), 6355–6367. doi: 10.1143/JJAP.44.6355.

41.

Tao

S.W.

, Tian

Z.H.

, Ba

et al. Grating-based X-ray phase contrast imaging with a polycapillary semilens, IEEE Photon Technol Lett 34(20) (2022), 1104–1107. doi: 10.1109/LPT.2022.3203012.

42.

Huang

N.E.

, Shen

, Long

S.R.

A new view of nonlinear water waves: The Hilbert spectrum, Annu Rev Fluid Mech 31 (1999), 417–457. doi: 10.1146/annurev.fluid.31.1.417.

43.

Huang

N.E.

, Shen

, Long

S.R.

et al. The empirical mode decomposition and the Hilbert spectrum for nonlinear and non-stationary time series analysis, P Roysoc A-Math Phy 454 (1998), 903–995. doi: 10.1098/rspa.1998.0193.

44.

Nunes

J.C.

, Bouaoune

, Delechelle

et al. Image analysis by bidimensional empirical mode decomposition, Image Vision Comput 21(12) (2003), 1019–1026. doi: 10.1016/S0262-8856(03)00094-5.

45.

Damerval

, Meignen

, Perrier

A fast algorithm for bidimensional EMD, IEEE Signal Process Lett 1(10) (2005), 701–704. doi: 10.1109/LSP.2005.855548.

46.

Zhu

, Liu

, Zhou

X.L.

Infrared polarization and intensity image fusion based on bivariate BEMD and sparse representation, Multimed Tools Appl 80(3) (2021), 4455–4471. doi: 10.1007/s11042-020-09860-z.

47.

H.Y.

, Wang

C.J.

, Zhao

An improved EMD and its applications to find the basis functions of EMI signals,, Math Probl Eng (2015), doi: 10.1155/2015/150127.

48.

Long

S.R.

Applications of HHT in image analysis, Interdisciplinary Mathematical Sciences 5 (2005), 289–305. doi: 10.1142/9789812703347_0013.

49.

Bhuiyan

S.M.A.

, Adhamii

R.R.

, Khan

J.F.

Fast and adaptive bidimensional empirical mode decomposition using order-statistics filter based envelope estimation, EURASIP J Adv Sig Pr 728356 (2016), doi: 10.1155/2008/728356.

50.

Nunes

J.C.

, Niang

, Bouaoune

, Delechelle

et al. Bidimensional empirical mode decomposition modified for texture analysis, Lect Notes Artif Intell, 2749 (2003), 171–177. doi: 10.1007/3-540-45103-X_24.

51.

G.T.

, Sang

E.F.

, Liu

Z.F.

et al. BEMD criterion based on the characteristic points’ number and positions, CISP IEEE 462 (2008), doi: 10.1109/CISP.2008.509.

52.

Yan

J.H.

, Lu

Improved Hilbert-Huang transform based weak signal detection methodology and its application on incipient fault diagnosis and ECG signal analysis, Signal Process 98 (2014), 74–87. doi: 10.1016/j.sigpro.2013.11.012.

53.

Cai

H.Y.

, Zhuo

L.R.

, Chen

X.D.

et al. Infrared and visible image fusion based on BEMSD and improved fuzzy set, Infrared Phys Technol 98 (2019), 201–211. doi: 10.1016/j.infrared.2019.03.013.

54.

Zhu

, Liu

, Zhou

X.L.

Infrared polarization and intensity image fusion based on bivariate BEMD and sparse representation, Multimed 80 (2021), 1–17. doi: 10.1007/s11042-020-09860-z.

55.

Y.Z.

, Lu

R.F.

Fast Bi-dimensional empirical mode decomposition as an image enhancement technique for fruit defect detection, Comput Electron Agric 152 (2018), 314–323. doi: 10.1016/j.compag.2018.07.025.

56.

Zhang

, Wang

C.S.

, Chen

G.T.

et al. Infrared and visible image fusion via NSCT and gradient domain PCNN, Proceedings of SPIE 12065 (2021), doi: 10.1117/12.2606088.

57.

Monica

, Kumar

Pulse coupled neural networks and its applications, Expert Systems With Applications 41(8) (2014), 3965–3974. doi: 10.1016/j.eswa.2013.12.02710.1016/j.eswa.2013.12.027.

58.

Ying

, Chen

S.H.

, Pan

G.J.

, He

Z.M.

Artificial intelligence pulse coupled neural network algorithm in the diagnosis and treatment of severe sepsis complicated with acute kidney injury under ultrasound image, Journal of Healthcare Engineering 2021 (6761364), doi: 10.1155/2021/6761364.

59.

Liu

, Liu

, Li

D.F.

Recent advances in pulse-coupled neural networks with applications in image processing, Electronics 11(20) (2022), doi: 10.3390/electronics11203264.

60.

Guo

, Li

, Shao

Z.Q.

et al. Objective image fusion evaluation method for target recognition based on target quality factor, Multimedia Syst 28(2) (2022), 495–510. doi: 10.1007/s00530-021-00850-1.

61.

Liu

S.Q.

, Lu

Y.C.

, Wang

et al. A new focus evaluation operator based on max–min filter and its application in high quality multi-focus image fusion, Multidim Syst Sign P 31(2) (2020), 569–590. doi: 10.1007/s11045-019-00675-2.

62.

Azam

M.A.

, Khan

K.B.

, Salahuddin

et al. A review on multimodal medical image fusion: Compendious analysis of medical modalities, multimodal databases, fusion techniques and quality metrics, 53), Comput Biol Med 144 (2022), doi: 10.1016/j.compbiomed.2022.105253.

63.

Sunitha

T.O.

, Rajalakshmi

Multi-modal image fusion technique for enhancing image quality with multi-scale decomposition algorithm, Comput Method Biomec 9(2) (2020), 1–13. doi: 10.1080/21681163.2020.1830437.

64.

You

C.S.

, Yang

S.Y.

A simple and effective multi-focus image fusion method based on local standard deviations enhanced by the guided filter, Displays 72 (2022), doi: 10.1016/j.displa.2021.102146.

65.

X.X.

, Jin

H.Y.

, Sun

et al. Multispectral image fusion based on laplacian pyramid decomposition with immune co-evolutionary optimization, LNAI 9242 (2015), 583–591. doi: 10.1007/978-3-319-23989-7_59.

66.

Long

, Peng

Y.X.

Blind fusion of hyperspectral multispectral images based on matrix factorization, Remote Sens 13(4219) (2021), doi: 10.3390/rs13214219.

67.

Liu

Y.Q.

, Sheng

Z.H.

, Shen

H.L.

Guided image deblurring by deep multi-modal image fusion,, IEEE Access 10 (2022), 130708–130718. doi: 10.1109/ACCESS.2022.3229056.