Imaging model for the scintillator and its application to digital radiography image enhancement

Qian Wang; Yining Zhu; Hongwei Li

doi:10.1364/OE.23.033753

1. Introduction

Microscopic computed tomography (Micro-CT) systems have been widely used in many fields such as biomedicine, materials and electronic packaging. Two kinds of detectors [1, 2], flat panel detector (FPD) and optical coupling detector (OCD) are typically used in Micro-CT systems. FPD, consisting of an x-ray scintillator and amorphous-silicon photodiodes, shown in Fig. 1(a), adopts the essential idea of geometrical magnification for high resolution imaging [3]. The OCD is based on optical magnifying principle, and consists of an x-ray scintillator, a set of optical lenses and a Charge-coupled Device (CCD) camera (or CMOS/sCMOS) [4], seen Fig. 1(b). Fig. 1(c) illustrates a real OCD-based Micro-CT system called nanoVoxel-2700 and produced by Sanying Precision Instrucments Ltd. Generally speaking, although the efficiency of OCD is much lower than FPD, the imaging resolution is much higher [5]. The OCD-based image formation can be divided into three stages: the conversion of x-ray beam to visible light photons by scintillator, the magnification of visible light by optical lenses and the capture of visible light by CCD-camera [6]. The combination of magnifying optical lenses and high resolution CCD camera can achieve the resolution of micron or submicron.

Fig. 1 Imaging systems of Micro-CT: (a) FPD-based system schematic; (b) OCD-based system schematic; (c) Real OCD-based system (model: nanoVoxel-2700; manufacture: Sanying Precision Instrucments Ltd.).

Download Full Size | PDF

Scintillator and optical devices are widely employed in indirectly detecting high-energy radiations in decades. Daugherty et al. [7] utilized scintillator in magnetic spectrometer to count signals efficiently. Then, in the study of monoenergetic synchrotron radiation based Micro-CT, Bonse et al. [8, 9, 10] employed fluorescent screen in x-ray-to-light conversion, which is 10 to 50 times more efficient than single crystal scintillator. Moreover, two-stage optical amplifi-cation mechanism, x-ray-optical magnification by magnifying groove crystal and light-optical magnification by lens system, is used to achieve higher spatial resolution. To achieve submicron resolution x-ray imaging, Koch et al. [11] combined transparent luminescent screens with optical microscope system to magnify the x-ray image onto a CCD camera. Parallel x-ray beam impinges perpendicularly onto the scintillator, then visible-light image is generated. A microscope objective magnifies the object plane, and a photo eyepiece corrects for field curvature and magnifies the object plane further. Between the lens elements, a real inverted image is formed. In order to improve the potential accessibility of high-resolution detector in existing facilities with low cost, Williams et al. [12] designed a system with long-working distance lenses and applied it in neutron imaging. The proposed infinity corrected optical design produces parallel light between the objective and imaging lenses, which is different from Kochs work. In addition, all the micro-CT systems above employ a mirror to remove visible light from the x-ray path before entering optical lenses. Our study is based on the detecting system shown in Fig. 2, which is somewhat different from the apparatuses above. Visible light is obtained from the conversion of cone-beam x-ray by scintillator, then magnified and parallelized by objectives. After reflected by a mirror, a tube lens makes the parallel light focus on the CCD-camera.

Fig. 2 The optical magnifying process of OCD-based system.

Download Full Size | PDF

The DR images obtained by OCD-based systems may suffer from low contrast, especially when the low-Z materials, e.g. core samples and bio-samples, are tested. It is because the voltage of x-ray source in laboratorial Micro-CT is usually higher than 20 kV, which it too high to obtain high-contrast DR images. In this case, some kinds of enhancement methods are needed to improve the image quality. However, conventional contrast enhancement methods were proposed for natural images and may no longer suit for DR images. For example, piecewise-linear intensity level transformation [13] is one of the most popular method for contrast-stretching. By selective line stretch, it expands the range of centralized intensity levels and enhances the contrast of interest. This method maps the intensity levels in the whole image and may do a poor job in local patches. Contrast limited adaptive histogram equalization algorithm (CLAHE) [14, 15] can obtain an image that shows a great deal of intensity levels in each local patch, based on local histogram equalization. Meanwhile, by limiting the height of the local histogram, possible noise amplification can be avoided. However, the CLAHE method pays more attention to the local histogram equalization rather than the global intensity varieties, so that the intensity of the processed image may not be loyal to real absorption values.

Meanwhile, because of the blurring effect, spatial resolution of DR images is decreased. This attribute can be described by the point spread function (PSF) mathematically. Generally speaking, according to whether prior knowledge, e.g. blurring kernel, is provided, the approaches for resolution improvement can be classified into two classes, non-blind deblurring [16] and blind deblurring [17, 18, 19]. A typical strategy in blind deblurring is to estimate the blurring kernel first. Thus, the problem is converted into a non-blind deblurring which is a inverse problem in mathematics. Due to the ill-posedness, some conventional approaches, such as inverse filtering, Wiener filtering [20] and constrained least squares filtering [13], cannot perform well. Thus, some regularization techniques, such as Tikhonov [21] and total variantion (TV) [22, 23], are employed to ensure the stability of the deblurring process. Particularly, TV based methods can effectively suppress the noise while maintain sharp edges.

To overcome the limitations of conventional contrast enhancement methods for DR images, we studied the image formation in the scintillator of OCD-based systems. Based on classical radiative transfer equation (RTE), the transferring process of visible light in the scintillator is analyzed. An image model is then established, which simulates the image formation under real lightening conditions, i.e., multiple light sources along a single path while the light spectra varies in different paths. It turns out that the resulting imaging model is similar to the haze imaging model. For haze imaging model, there is an assumption that the specific intensity of the incident light should be uniformly distributed. This is not the case for the Micro-CT system with cone-beam x-rays and flat-type detectors. To meet this assumption, a pre-correction technique is proposed and applied to DR images before utilizing the contrast enhancement methods designed for the haze imaging model.

In [24], dark channel prior was introduced by He et al. to solve the haze imaging model based inverse problem effectively. This idea was further elaborated in [25]. The rationale behind is that most local patches in haze-free outdoor images contain some pixels with low intensities near zero in at least one color channel [24]. In Micro-CT systems, the specific intensity as well as the energy of the x-rays are rather low compared to industrial CT systems, e.g., micro-CT power range is 10–40 W and industrial CT power range is generally 1–2 kW. This leads to low transmission, and a small fraction of high density materials could lead to very low intensity values in DR images. Usually, the density distribution of the objects under DR scanning, such as core samples and bio-samples, could be thought of being homogeneous. So low intensity values in DR images should also be homogeneously distributed, which means most local patches in DR images contain some pixels with low intensities near zero. This exactly corresponds to the Dark Channel Prior. So the methods proposed in [24] and [25] are extended to solve the inverse problem based on the proposed imaging model.

Further, in our circumstances, the PSF in scintillator is unknown, so effective blind deblurring method is needed. According to our previous study [26], a effective blincd deblurring method based on ENR maximization principle [27, 28, 29] and modified ROF deblurring model [26] is employed to improve the resolution of DR images.

The remainder of this paper is organized as follows. In section 2, We introduce classical RTE and give some explanations. In section 3, the scintillator imaging model is developed. To solve the correlative inverse problem, in section 4, pre-correction technique is proposed to reduce the intensity inhomogeneity of the x-rays. Then, dark channel prior based haze removal method shall be employed to improve the contrast. Further, to reduce the blurring effect, blind deblurring method based on ENR maximization principle and modified ROF deblurring model is utilized. In section 5, experiments are performed to compare the proposed method with other methods, which validate the effectiveness of our approach. In section 6, we conclude this paper.

2. Primary knowledge of RTE

The propagation of visible light in scintillator can be recognized as a radiative transfer process based on the assumption that the refractive index of the media is constant [30, 31]. In order to establish the imaging model in scintillator, we start by explaining the mathematical model of radiative transfer process of visible light, i.e., radiative transfer equation:

\begin{array}{l} \frac{1}{c} \cdot \frac{\partial I_{λ} (\vec{r}, \vec{v}, t)}{\partial t} & = - \vec{v} \cdot \vec{\nabla} I_{λ} (\vec{r}, \vec{v}, t) - b_{abs} (λ, \vec{r}) \cdot I_{λ} (\vec{r}, \vec{v}, t) - b_{sca} (λ, \vec{r}) \cdot I_{λ} (\vec{r}, \vec{v}, t) \\ + b_{sca} (λ, \vec{r}) \int \frac{P (λ, \vec{r}, \vec{u}, \vec{v})}{4 π} \cdot I_{λ} (\vec{r}, \vec{u}, t) d \vec{u} + S_{λ} (\vec{r}, \vec{v}, t), \end{array}

where the symbols are defined in Table 1.

Table 1. Symbol definition for Eq. (1)

View Table | View all tables in this article

The radiative transfer Eq. (1) mathematically describes the physical phenomenon of energy transfer in the form of electromagnetic radiation. One of its derivation methods, i.e., Eulerian method [32], considers the change of the photon number with respect to time in a fixed five-dimensional “cube”, of which the volume is

Δ V = Δ x Δ y Δ z Δ μ Δ ϕ,

where μ ≡ cosθ and θ stands for polar angle; ϕ is azimuthal angle; x, y and z represent Cartesian coordinates. Moreover,

\begin{array}{l} d \vec{r} = d x d y d z, \\ d \vec{v} = sin θ d θ d ϕ = d μ d ϕ . \end{array}

According to the conservation law, one can obtain the following equation

\begin{array}{l} \frac{\partial ξ_{λ} (\vec{r}, \vec{v}, t)}{\partial t} Δ V & = - \vec{v} \cdot \vec{\nabla} ξ_{λ} (\vec{r}, \vec{v}, t) Δ V - b_{abs} (λ, \vec{r}) \cdot ξ_{λ} (\vec{r}, \vec{v}, t) Δ V - b_{sca} (λ, \vec{r}) \cdot ξ_{λ} (\vec{r}, \vec{v}, t) Δ V \\ + b_{sca} (λ, \vec{r}) \int \frac{P (λ, \vec{r}, \vec{u}, \vec{v})}{4 π} \cdot ξ_{λ} (\vec{r}, \vec{u}, t) d \vec{u} Δ V + s_{λ} (\vec{r}, \vec{v}, t) Δ V, \end{array}

where ξ_λ(r⃗, x⃗, t) is the distribution function and s_λ(r⃗, v⃗, t) denotes the averaged emission density. Each component of the right-hand side of Eq. (2) corresponds to a physical process involved in the propagation of radiation through a medium, such as absorption, emission and scattering, which are specified in Table 2.

Table 2. the physical process in a volume element

View Table | View all tables in this article

The relationship between the intensity I_λ(r⃗, v⃗, t) and the distribution function ξ_λ(r⃗, x⃗, t) reads

ξ_{λ} (\vec{r}, \vec{v}, t) = I_{λ} (\vec{r}, \vec{v}, t) \cdot \frac{λ}{h \cdot c^{2}},

where h = 6.6262 × 10⁻³⁴Js is the Planck’s constant and c represents the speed of light in vacuum. By substituting Eq. (3) into Eq. (2) and defining the function of the source photons as

S_{λ} (\vec{r}, \vec{v}, t) = h \cdot \frac{c}{λ} \cdot s_{λ} (\vec{r}, \vec{v}, t),

Eq. (1) can be derived.

3. Scintillator imaging model

In scintillator, visible light photons are generated by the radiative de-excitation of the excited atoms when penetrated by x-rays [33, 34, 35, 36]. Because of the thin thickness of the scintillator, the attenuation effect on x-rays is reasonably negligible. Moreover, x-ray source is heterogenetic in practice, and the spectrum distribution of object-attenuated x-rays are different along dissimilar path. Thus, It is reasonable to assume that the intensity of the visible light is the same along a single x-ray path while could vary along different paths. Based on this observation, Eq. (1) can be simplified. To do so, we put Eq. (1) into a spatial Euclidean coordinate system and analyze each physical process described in Table 2. As shown in Fig. 3, the center of the incident surface of the scintillator is taken to be the origin, the x-axis specifies the source-detector direction and the z-axis is perpendicular to the plane of the trajectory. Due to the assumption of homogeneous medium, the visible light sources can be thought of being approximated uniformly distributed. This allows us to segment the scintillator into cubic elements. As shown in Fig. 3, each of the cubic element contains only one visible light source. For a single element, along the travel direction v⃗, the physical processes specified in Table 2 could be analyzed as follows:

streaming The photon number of in-streaming is slightly larger than that of out-streaming, which means the photon number of net-streaming from a pointolite is always positive. This is because the solid angle of the incidence surface is larger than that of the emergence surface, which is demonstrated on the left part of Fig. 4. Moreover, the photon number of net-streaming decreases as the pointolite-element distance increases. Hence, we just consider the visible radiation from the pointolites located within a specified distance and neglect the others outside of this distance. Depending on the assumption of uniformly distributed visible light sources, for an interior element, i.e. case A in Fig. 4, the photon number of net-streaming can be reasonably approximated as zero. For any considered pointolite (located in the orange circle), its centrosymmetric pointolite about the element is still located in the circle. Thus, the photon number of net-streaming caused by this pair of pointolites is approximately zero. For the elements near the boundaries, i.e. case B in Fig. 4, the net-streaming cannot be assumed zero. However, the number of boundary elements is very few compared to the number of total elements. So their contributions could be neglected when studying the propagation of visible light in the scintillator.
absorption and out-scattering The photons propagating along v⃗ decrease due to absorption and out-scattering (primary scattering) [37]. Since the absorption effect of the scintillator to visible light is negligible compared with the scattering effect, the absorption-caused attenuation can be ignored. Thus, we approximate the scattering coefficient as the extinction coefficient. According to the Lambert-Beer Law, this attenuation depends on the wavelength and the propagation distance.
in-scattering Each emitted photon is subjected to multiple scattering. Due to the randomness of the scattering direction, one could assume that all the elements along the direct v⃗ are exposed to an uniformly distributed field of scattered photons.
emission The photons in an element are mainly emitted by the central light source. Thus, the distribution function for the volume element could be approximated by the averaged density function determined by the central light source.

Based on the above analyses, Eq. (2) can be approximated as

\begin{array}{l} \frac{\partial ξ_{λ} (\vec{r}, \vec{v}, t)}{\partial t} Δ V & = - b_{sca} (λ, \vec{r}) \cdot s_{λ} (\vec{r}, \vec{v}, t) Δ V + b_{sca} (λ, \vec{r}) \cdot a_{\bar{λ}} (t) Δ V + s_{λ} (\vec{r}, \vec{v}, t) Δ V \\ = (1 - b_{sca} (λ, \vec{r})) \cdot s_{λ} (\vec{r}, \vec{v}, t) Δ V + b_{sca} (λ, \vec{r}) \cdot a_{\bar{λ}} (t) Δ V, \end{array}

where, a_λ̄ is the average multiple scattering density function. By substituting Eq. (3) and Eq. (4) into Eq. (5), Eq. (1) can be simplified as

\frac{1}{c} . \frac{\partial I_{λ} (\vec{r}, \vec{v}, t)}{\partial t} = (1 - b_{sca} (λ, \vec{r})) \cdot S_{λ} (\vec{r}, \vec{v}, t) + b_{sca} (λ, \vec{r}) \cdot A_{\bar{λ}} (t),

where

A_{\bar{λ}} (t) = h \cdot \frac{c}{λ} \cdot a_{\bar{λ}} (t),

is the function of multiple scattering.

Fig. 3 Illustration of the multiple light sources model. A local box of the scintillator is magnified and shown in the yellow wireframe.

Download Full Size | PDF

Fig. 4 Illustration of streaming analyses. The difference between in-streaming and out-streaming is shown on the left, which decreases with the increase of the pointolite-element distance. As is shown on the right, the elements of the scintillator can be seperated into two cases: (A) interior element; (B) boundary element.

Download Full Size | PDF

For any element, its energy transfer is approximately isotropous. However, due to the flatness of the CCD array, the detected intensity for each detection unit varies. Especially, the highest specific intensity is located at the same ray-penetrated detection unit (only x coordinate of the element and the detection unit is different), and the intensities for other detection units are relatively lower, due to the following reasons: a. the solid angle of other detection unit is smaller than the same ray-penetrated detection; b. the light path to other detection unit is longer than the same ray-penetrated detection, which leads to more attenuation; c. the light path to other detection unit is not perpendicular to the emergent surface of the scintillator, so that reflection and refraction will happen; d. since the scintillator is optically denser medium comparing to air, when the incident angle is greater than the critical angle, total reflection will happen.

Admittedly, for any element, excluding its intensity mainly obtained by the detector unit in the same ray path, some radiation influences the neighbor detection units. Inspired by our previous study on optical blurring in [26], this effect attributing to the reduction of resolution can be described by PSF mathematically. Then, the function of PSF as a linear convolution operation, makes the complicated image formation in scintillator separable: contrast decrease subprocess and resolution reduction subprocess.

3.1. Contrast decrease subprocess

In the contrast decrease subprocess, we consider the radiative transfer process along a single x-ray path. Thus, for any element, we only focus on its intensity for the same ray-penetrated detection unit, i.e., the case that the travel direction is x⃗ = (1, 0, 0). The intensity for each detector unit can be obtained by integration along the corresponding x-ray path in the scintillator. According to Eq. (3), the photon distribution is proportional to intensity, so that the photon number and the integral of intensity maintain the same relation. Especially, the photon number can be calculated by summing the photon number of all the same ray-penetrated elements. The photon number of each element is the product of the corresponding average photon distribution function and the volume of the cube. In an exposure time, the photon numbers collected by the detector units correspond to the DR data, i.e.,

I (\vec{r}) = \int_{t}^{t + Δ t} \sum_{i = 0}^{N - 1} \frac{\partial ξ_{λ} (\vec{r} + i Δ d \vec{x}, \vec{x}, \hat{t})}{\partial t} Δ V d \hat{t},

where Δt is the exposure time, r⃗ + iΔdx⃗ is the location of the i-th element (i = 0, 1,···, N − 1) in the x-ray path noted as r⃗ + lx⃗, l ∈ (−∞, ∞). As is shown in Fig. 5, the thickness of scintillator is defined as d, and the center-to-center distance between adjacent elements in the same path is noted as Δd. N = d/Δd is the number of elements in a single x-ray path. The travel distance of the visible light generated in the i-th element can be expressed as (d − iΔd). Particularly, along a single x-ray path, the average emitted photon distribution of each element is approximately identical, which can be noted as s(r⃗). By substituting Eq. (5) and s(r⃗), we rewrite Eq. (7) as

\begin{array}{l} I (\vec{r}) & = \int_{t}^{t + Δ t} \sum_{i = 0}^{N - 1} [(1 - b_{sca} (λ, \vec{r} + i Δ d \vec{x})) \cdot s_{λ} (\vec{r} + i Δ d \vec{x}, \vec{x}, \hat{t}) Δ V + b_{sca} (λ, \vec{r} + i Δ d \vec{x}) \cdot a_{\bar{λ}} (\hat{t}) Δ V] d \hat{t} \\ = \int_{t}^{t + Δ t} [s_{λ} (\vec{r}, \vec{x}, \hat{t}) Δ V \sum_{i = 0}^{N - 1} (1 - b_{sca} (λ, \vec{r} + i Δ d \vec{x})) + a_{\bar{λ}} (\hat{t}) Δ V \sum_{i = 0}^{N - 1} b_{sca} (λ, \vec{r} + i Δ d \vec{x})] d \hat{t} \\ = \int_{t}^{t + Δ t} [s_{λ} (\vec{r}, \vec{x}, \hat{t}) N Δ V + (a_{\bar{λ}} (\hat{t}) - s_{λ} (\vec{r}, \vec{x}, \hat{t})) Δ V \sum_{i = 0}^{N - 1} b_{sca} (λ, \vec{r} + i Δ d \vec{x})] d \hat{t} \\ = \int_{t}^{t + Δ t} [s_{λ} (\vec{r}, \vec{x}, \hat{t}) N Δ V + (a_{\bar{λ}} (\hat{t}) - s_{λ} (\vec{r}, \vec{x}, \hat{t})) N Δ V \frac{1}{N} \sum_{i = 0}^{N - 1} μ {(\vec{r})}^{N - i}] d \hat{t}, \end{array}

where

μ {(\vec{r})}^{N - i} = exp (- β (\vec{r}) (N - i) Δ d) = exp (- β (\vec{r}) (d - i Δ d)) \approx b_{sca} (λ, \vec{r} + i Δ d \vec{x}),

is an approximately exponential model corresponding to the Lambert-Beer Law, β(r⃗) is the “average” attenuation coefficient of the visible light with average wavelength along the x-ray path defined as r⃗ + lx⃗, l ∈ (−∞, ∞). Summing out Eq. (9), we have

\begin{array}{l} \frac{1}{N} \sum_{i = 0}^{N - 1} μ {(\vec{r})}^{N - i} & = \frac{μ (\vec{r}) (1 - μ {(\vec{r})}^{N})}{N (1 - μ (\vec{r}))} \\ = \frac{exp (- β (\vec{r}) Δ d) (1 - exp (- β (\vec{r}) N Δ d))}{N (1 - exp (- β (\vec{r}) Δ d))} \\ = \frac{exp (- β (\vec{r}) \frac{d}{N}) (1 - exp (- β (\vec{r}) d))}{N (1 - exp (- β (\vec{r}) \frac{d}{N}))} . \end{array}

Fig. 5 Illustration of the imaging model in a single x-ray path. Depending on the multiple light sources model shown in Fig. 3, the elements are uniformly distributed. Along a single x-ray path, the average emitted distribution function of each element is equal.

Download Full Size | PDF

Due to the thin thickness of the scintillator and the large number of elements along a x-ray path, d/N ≈ 0. Thus, Eq. (10) can be approximated as

\begin{array}{l} \frac{1}{N} \sum_{i = 0}^{N - 1} μ {(\vec{r})}^{N - i} & \approx \frac{exp (- β (\vec{r}) \frac{d}{N}) (1 - exp (- β (\vec{r}) d))}{N [1 - (1 - β (\vec{r}) \frac{d}{N})]} \\ = \frac{1 - exp (- β (\vec{r}) d)}{β (\vec{r}) d} . \end{array}

By replacing exp(−β(r⃗)d) with its second order approximation, Eq. (11) can be further approximated as

\begin{array}{l} \frac{1}{N} \sum_{i = 0}^{N - 1} μ {(\vec{r})}^{N - i} & \approx \frac{1 - [1 - β (\vec{r}) d + \frac{1}{2} {(β (\vec{r}) d)}^{2}]}{β (\vec{r}) d} \\ = 1 - \frac{1}{2} β (\vec{r}) d . \end{array}

Substituting Eq. (12) into Eq. (8), we get

\begin{array}{l} I (\vec{r}) & = \int_{t}^{t + Δ t} [s_{λ} (\vec{r}, \vec{x}, \hat{t}) N Δ V + (a_{\bar{λ}} (\hat{t}) - s_{λ} (\vec{r}, \vec{x}, \hat{t})) N Δ V (1 - \frac{1}{2} β (\vec{r}) d)] d \hat{t} \\ = \int_{t}^{t + Δ t} [\frac{1}{2} β (\vec{r}) d s_{λ} (\vec{r}, \vec{x}, \hat{t}) N Δ V + (1 - \frac{1}{2} β (\vec{r}) d) a_{\bar{λ}} (\hat{t}) N Δ V] d \hat{t} \\ = B (\vec{r}) \int_{t}^{t + Δ t} s_{λ} (\vec{r}, \vec{x}, \hat{t}) N Δ V d \hat{t} + (1 - B (\vec{r})) \int_{t}^{t + Δ t} a_{\bar{λ}} (\hat{t}) N Δ V d \hat{t} \\ = B (\vec{r}) J (\vec{r}) + (1 - B (\vec{r})) K, \end{array}

where

B (\vec{r}) = \frac{1}{2} β (\vec{r}) d, J (\vec{r}) = \int_{t}^{t + Δ t} s_{λ} (\vec{r}, \vec{x}, \hat{t}) N Δ V d \hat{t}, K = \int_{t}^{t + Δ t} a_{\bar{λ}} (\hat{t}) N Δ V d \hat{t} .

If the scattering process is nonexistent, b_sca(λ, r⃗ + iΔdx⃗) in Eq. (8) can be treated as 0 for i = 0,···, N − 1. Based on this assumption, from Eq. (8), we can deduce that

I {(\vec{r})}^{*} = \int_{t}^{t + Δ t} s_{λ} (\vec{r}, \vec{x}, \hat{t}) N Δ V d \hat{t} = J (\vec{r}),

which indicates that J(r⃗) is the ideal DR data without scattering-reduced contrast. In addition, B(r⃗) in Eq. (13) represents the effective transmission of visible light in the scintillator, and K is the photon number corresponding to the multiple scattering process.

In the outdoor image formation, the atmospheric absorption and scattering caused by turbid medium (e.g., dust, water-drops, particles) lead to hazed images with low contrast. As shown in Fig. 6, the light received by camera consists of the reflected light from the scene point and the ambient light [38]. Particularly, the reflected light is attenuated along the light path. This imaging model can be expressed as [39, 40, 24, 41, 42]:

I (t) = t (t) J (r) + (1 - t (r)) A,

which is called the haze imaging model and widely used in computer vision and computer graphics. In Eq. (14), r represents the position of a pixel, I is the observed intensity, J is the scene radiance, A is the ambient light, and t is a scalar in [0, 1], which presents the transmission of scene radiance. Although the mathematical models of the haze imaging model and the model in Eq. (13) are similar, the analyses and derivations are quite different.

Fig. 6 Illustration of haze imaging model. The light received by camera consists of the reflected light from the scene point and the ambient light. Particularly, the reflected light is attenuated along the light path.

Download Full Size | PDF

3.2. Resolution reduction subprocess

In the resolution reduction subprocess, we focus on the radiation influence from the other different ray-penetrated elements, and combine it into model Eq. (13). Conventionally, the PSF in image formation can be approximately characterized as Gaussian function. By denoting the Gaussian function as G_σ₀, of which the standard deviation is σ₀, we can get the blurred version of Eq. (13) as follows,

\begin{array}{l} \dot{I} (\vec{r}) & = \int_{\vec{x} \in Ω (\vec{r})} G_{σ_{0}} (\vec{x} - \vec{r}) I (\vec{x}) d \vec{x} = G_{σ_{0}} (\vec{r}) \otimes I (\vec{r}), \\ = G_{σ_{0}} (\vec{r}) \otimes (B (\vec{r}) J (\vec{r}) + (1 - B (\vec{r})) K), \end{array}

where Ω(r⃗) is a neighborhood of position r⃗; the ⊗ symbol stands for convolution operation; İ(r⃗) is the DR image with reduced resolution. We call Eq. (15) the scintillator imaging model.

In real CT system, DR images are usually polluted by Poisson noise [43]. Thus, Eq. (15) can be further transformed into the following formula,

\overset{̆}{I} (\vec{r}) ~ Poission (\dot{I} (\vec{r})) = Poission (G_{σ_{0}} (\vec{r}) \otimes (B (\vec{r}) J (\vec{r}) + (1 - B (\vec{r})) K)),

where Ĭ(r⃗) is the degraded DR images with low contrast, reduced resolution, and noise influence. The goal of enhancement is to recover J from Ĭ, which is an inverse problem in mathematics.

4. Algorithms

In this section, some algorithms are developed and extended to solve the scintillator imaging model based inverse problem. First, an innovative pre-correction technique is proposed to reduce the intensity inhomogeneity of the x-rays of Micro-CT systems. Then, the solution of the scintillator imaging model is divided into two processes: to improve the contrast, the haze removal method based on dark channel prior [24, 25] is employed; to reduce the blurring effect, the blind deblurring method [26] based on ENR maximization principle and modified ROF deblurring model is utilized.

4.1. Pre-correction for intensity of x-ray

In the haze imaging model, the reflected light comes from the sun at infinity, whose intensity follows uniform distribution. In Micro-CT systems, however, the energy of cone-beam x-rays captured by the flat array detectors is nonuniform. To eliminate the bright field effect, additional DR image of original specific intensity is needed conventionally, which doubles both the scan number and the experimental time. For the sake of reducing experimental cost and maintaining same performance, a pre-correction technique, based on quintic polynomial approximation of the intensity distribution, is applied to the DR images. It is generally accepted that the original intensity of x-rays obeys an approximate Gaussian distribution [44]. However, experiments show that Gaussian model can not describe the variation of the intensity sufficiently.

In the following experiments, we use the DR image of original intensity collected by OCD-based system with 4× lens, shown in Fig. 7(a). The fitting results of these two models are illustrated in Fig. 7(b) and Fig. 7(c). To evaluate the fitness, we compare their ratio images and difference images separately. In the analysis of ratio images, Fig. 8(a) has some obvious variations. Especially, the middle area circled in red and its corresponding the bulged area circled in blue shown in Fig. 7(b), shows a weak fit of Gaussian function. On the contrary, quintic polynomial has a better description, seen from Fig. 8(b). When comparing the difference images in Fig. 9, we can find that the difference of Gaussian function is higher than quintic polynomial. As a result, quintic polynomial is more reasonable in pre-correction.

Fig. 7 The DR image of original specific intensity and its pre-correction fitting surfaces: (a) the image of original specific intensity; fitting surface by using (b) Gaussian function and (c) quintic polynomial.

Download Full Size | PDF

Fig. 8 Ratio images between Fig. 7(a) and its fitting surfaces by using (a) Gaussian function and (b) quintic polynomial.

Download Full Size | PDF

Fig. 9 Difference images between Fig. 7(a) and its fitting surfaces by using (a) Gaussian function and (b) quintic polynomial.

Download Full Size | PDF

In order to compare the performance of pre-correction with ideal result in standard flat field, a spider DR image obtained in the same magnification is employed, shown in Fig. 10(a). The ideal result in standard flat field and quintic polynomial based pre-correction result are shown in Fig. 10(b) and Fig. 10(c) respectively. From the perspective of visional evaluation, both of them have uniform background. Moreover, the correlative histograms in Fig. 10(e) and Fig. 10(f) suggest analogous performances as well.

Fig. 10 Compared illustrations of pre-correction: (a) original spider DR image; (b) ideal result in standard flat field; (c) pre-correction result by using quintic polynomial. Correlative histograms are shown in (d)–(f).

Download Full Size | PDF

4.2. DR image enhancement

To solve scintillator imaging model based inverse problem effectively, we analyze the feature of the DR images which is caused by the Micro-CT systems. We find that, in Micro-CT systems, the intensity as well as the energy of the x-rays are rather low compared to industrial CT systems, which leads to low transmission. And a small fraction of high density materials could lead to very low intensity values in DR images. Moreover, the density distribution of the objects under DR scanning, such as core samples and bio-samples, could be usually thought of homogeneous. Fig. 11 shows typical samples, shale and bamboo stick, whose local patches in the DR image have similar structures and homogeneous density distribution. So the low intensity values in DR images should be homogeneously distributed, which means most local patches in DR images contain some pixels with low intensities near zero. This exactly corresponds to the Dark Channel Prior, an effective prior method to solve the haze imaging model based inverse problem. Dark channel prior, a kind of statistics of haze-free outdoor images, is based on a key observation - most local patches in haze-free outdoor images contain some pixels with low intensities near zero in at least one color channel [24]. Due to the similar prior, we extended and employed the methods proposed in [24] and [25] to solve our problem.

Fig. 11 Typical samples: shale (on the left) and bamboo stick (on the right).

Download Full Size | PDF

Formally, for any DR image J, its dark channel is defined as follows,

J^{dark} (\vec{r}) = min_{\vec{m} \in Ω (\vec{r})} J (\vec{m}),

where Ω(r⃗) is a local patch centered at r⃗. Without degradation, the dark channel of any pixel r⃗ in DR image tends to zero, so

J^{dark} (\vec{r}) \approx 0 .

According to the definition in Eq. (16), we compute the dark channel on both sides of the normalized Eq. (15) with the assumption that the transmission of a local patch Ω(r⃗) is a constant B̈(r⃗), i.e.,

\begin{array}{l} min_{\vec{m} \in Ω (\vec{r})} (\frac{\dot{I} (\vec{m})}{K}) & = min_{\vec{m} \in Ω (\vec{r})} (\frac{G_{σ_{0}} (\vec{m}) \otimes (\ddot{B} (\vec{m}) J (\vec{m}) + (1 - \ddot{B} (\vec{m})) K)}{K}) \\ = min_{\vec{m} \in Ω (\vec{r})} (\ddot{B} (\vec{m}) \frac{G_{σ_{0}} (\vec{m}) \otimes J (\vec{m})}{K} + 1 - \ddot{B} (\vec{m})) \\ = \ddot{B} (\vec{r}) min_{\vec{m} \in Ω (\vec{r})} \frac{G_{σ_{0}} (\vec{m}) \otimes J (\vec{m})}{K} + 1 - \ddot{B} (\vec{r}) . \end{array}

By concerning the local smoothness of DR images and the fact that K as the multiple scattering photon number along a x-ray path in the scintillator is always positive, according to Eq. (17), we have

min_{\vec{m} \in Ω (\vec{r})} \frac{G_{σ_{0}} (\vec{m}) \otimes J (\vec{m})}{K} \approx 0 .

If K is given, by evaluating Eq. (18) with Eq. (19), we can estimate the transmission B̈(r⃗) simply by

\ddot{B} (\vec{r}) = 1 - min_{\vec{m} \in Ω (\vec{r})} \frac{\dot{I} (\vec{m})}{K} .

Moreover, the dark channel of DR image can be expressed as

{\dot{I}}^{dark} (\vec{r}) = \ddot{B} (\vec{r}) {(G_{σ_{0}} \otimes J)}^{dark} (\vec{r}) + K (1 - \ddot{B} (\vec{r})) .

Based on Eq. (19), we obtain

{\dot{I}}^{dark} (\vec{r}) \approx K (1 - \ddot{B} (\vec{r})) .

As B̈(r⃗) is a scalar in [0, 1], according to Eq. (20), the lowest transmission B̈(r⃗) corresponds to the best approximation of the multiple scattering photon number along a x-ray path, i.e., K. Therefore, in applications, we pick the brightest pixel in the dark channel of DR images as the approximation of K. With K and B̈(r⃗), we can get

G_{σ_{0}} (\vec{r}) \otimes J (\vec{r}) = \frac{\dot{I} (\vec{r}) - K}{\ddot{B} (\vec{r})} + K .

Due to the locally-uniform-transmission assumption, transmission B̈(r⃗) in Eq. (21) can not conform to the actual situation. As a result, halo artifacts appear near edges in the recovered image J. In this case, the guided filter [25] based on the linear ridge regression model [45, 46], is applied to refine B̈(r⃗). It has also been proven useful in image matting [47] and image super-resolution [48]. With the guidance I and the filtering input B̈(r⃗), we can compute the filtering output B by

B_{i} = {\bar{a}}_{i} {\dot{I}}_{i} + {\bar{b}}_{i},

where

{\bar{a}}_{i} = \frac{1}{| ω |} \sum_{k \in ω_{i}} \frac{\frac{1}{| ω_{i} |} {\dot{I}}_{i} {\ddot{B}}_{i} - μ_{k} - \bar{{\ddot{B}}_{k}}}{σ_{k}^{2} + ε}, {\bar{b}}_{i} = \frac{1}{| ω |} \sum_{k \in ω_{i}} (\bar{{\ddot{B}}_{k}} - a_{k} μ_{k}),

μ_k and

σ_{k}^{2}

are the mean and variance of İ in ω_k, |ω| is the number of pixels in ω_k, and

\bar{{\ddot{B}}_{k}}

is the mean of B̈ in ω_k. In Appendix, we utilize the experiment to show the performance of guided filter.

In addition, when B is close to zero, the noise of recovered G_σ₀ ⊗ J is enhanced. So that, we improve Eq. (21) by

G_{σ_{0}} (\vec{r}) \otimes J (\vec{r}) = \frac{\dot{I} (\vec{r}) - K}{max (B (\vec{r}), B_{0})} + K,

where B₀ is a lower bound. Denoting

\dot{J} (\vec{r}) = \frac{\dot{I} (\vec{r}) - K}{max (B (\vec{r}), B_{0})} + K,

we have

G_{σ_{0}} \otimes J = \dot{J} .

According to Eq. (22), we can recovery J̇ from the achieved DR image Ĭ. In order to get undegraded image J from Eq. (23), a effective deblurring method is needed. In practice, the standard deviation of Gaussian function is unknown, and the noise in real CT system is subject to Poisson distribution. Thus, the blind deblurring method proposed in our previous work [26] can be employed, which solves this problem in 2 steps: the estimation of blurring kernel and the non-blind deblurring process.

In the step of estimating of blurring kernel, ENR maximization principle proposed in [28, 29, 27] is utilized, which is defined as

arg max_{σ} ENR (σ) = \frac{E (σ)}{N (σ)},

where

E (σ) = I (R_{σ} (J), G_{σ} \otimes R_{σ} (J)), N (σ) = I (J, G_{σ} \otimes R_{σ} (J)),

R_σ(.) is the deblurring operator, and I is the following discrepancy measure based on I-divergence (generalized Kullback distance) [49]

I (u, v) = \sum_{x} u (x) log \frac{u (x)}{v (x)} - \sum_{x} [u (x) - v (x)],

where u(x) and v(x) are two nonnegative distributions. ENR maximization principle reduces the parameter estimation problem into an 1-D optimization problem, which can be solved by dichotomy, golden section search and so on.

In the step of non-blind deblurring, for any estimated σ, a modified ROF deblurring model proposed in [26] and extended in [50] is employed, which reads

min {\frac{1}{2} {‖ \frac{1}{\sqrt{G_{σ} \otimes J}} (G_{σ} \otimes J - \dot{J}) ‖}_{L^{2}}^{2} + μ {‖ \nabla J ‖}_{L_{1}}},

where μ > 0, is a regularization parameter to balance regularity and fidelity; ∇ stands for gradient operator. This model inspired by the ROF deblurring model proposed by Rudin and Osher [23] and the Poisson-ROF model proposed by Zhu et al. [43] is designed in light of Gaussian blurring and Poisson noise. Based on nonlinearity approximation [26] and split Bregman iteration [51, 52], model (24) can be solved by Algorithm 1.

Algorithm 1. The Algorithm of Solving Model (24)

View Table | View all tables in this article

In Algorithm 1, Δ is the Laplace operator; the tilde symbol stands for adjoint operator; the hat symbol stands for fast Fourier transform; IFFT means inverse fast Fourier transform;

shrink (x, a) = {\begin{array}{l} x - a, & x \in (a, + \infty), \\ 0, & x \in [- a, a], \\ x + a, & x \in (- \infty, - a) . \end{array}

5. Experimental results

In order to prove the validity of our imaging model and the superiority of our enhancement method, we compare our method with piecewise-linear intensity level transformation and CLAHE, where blind deblurring is combined into piecewise-linear intensity level transformation method to improve its performance in spatial resolution. In our experiments, we utilize two DR images of typical biological samples, the spider DR image and the beetle DR image. Both of them are collected by nanoVoxel-2700 system manufactured by Sanying Precision In-strucments Ltd. The optical lens is 4× magnification and the voltage of x-ray source is 40 kV. Piecewise-linear intensity level transformation can be expressed as:

f (x) = {\begin{array}{l} a x, & 0 \leq x < d, \\ b x - c, & d \leq x \geq 255 . \end{array}

And in CLAHE algorithm, the limited ratio is 0.01, and the number of local patches is less than 8 × 8 with 256 bins for the histogram.

In the first experiment, we enhance the spider DR image shown in Fig. 10(a) by using blind deblurring, the combination of blind deblurring and piecewise-linear intensity level transformation, CLAHE and the proposed approach respectively. Correlative results, zoom-in patches and histograms are shown in Fig. 12. For piecewise-linear intensity level transformation, the parameters are set as a = 0.8, b = 1.33, c = 80, and d = 160. In our approach, patchSize = 15 × 15, μ = 800, and λ = 100. Deblurred result in Fig. 12(a) suggests a noticeable improvement in spatial resolution and a low contrast vision. By using piecewise-linear intensity level transformation, as is shown in Fig. 12(b), global contrast is stretched effectively. However, local patches still concentrated in small grey range and some grey levels are lost. According to the last two columns of Fig. 12, both CLAHE and the proposed approach have dramatic performances on global contrast-stretching and detail distinguish (inside the blue circle). Admittedly, the enhanced result by CLAHE shows obvious local contrast, significant physical information (inside the red oval-ring) is changed meanwhile. On the contrary, our method can keep the property of relative intensity and obtain a satisfying enhancement. Same profiles (row: 230; column: 610–620) are extracted from different results (marked in green), shown in Fig. 13(a). For better illustration, gain calibration is employed on each of them. It is obvious that the proposed method can achieve greater varieties in both spatial dimension and grey dimension than other compared methods.

Fig. 12 Enhanced results of Fig. 10(a) by using (a) blind deblurring; (b) the combination of blind deblurring and piecewise-linear intensity level transformation; (c) CLAHE; (d) the proposed method. The zoom-in patches of (a)–(d) (inside the red box) are illustrated in (e)–(h). The correlative histograms of (a)–(d) are illustrated in (i)–(l).

Download Full Size | PDF

Fig. 13 Profiles of enhanced results: (a) is spider DR image (marked by green); (b) is beetle DR image (marked by green).

Download Full Size | PDF

The second experiment is based on the beetle DR image shown in Fig. 14(a). The results of piecewise-linear intensity level transformation, combination of blind deblurring and piecewise-linear intensity level transformation, CLAHE, the proposed approach (deblurring free) and the proposed approach are illustrated in Fig. 14(b)–(f). The correlative histograms and zoom-in patches are shown in Fig. 14(g)–(l) and Fig. 15. For the proposed approach, the parameters are set as patchSize = 15 × 15, μ = 800, and λ = 100. Comparing Fig. 14(b) and Fig. 14(c), blind deblurring method improves the spatial resolution effectively, but both of them suffer low local contrast and lose some gray levels in histograms. These drawbacks reflect the limitation of piecewise-linear intensity level transformation. Coherently, CLAHE and the proposed method perform well in light of contrast-stretching. In magnified details, red arrows indicate the enhanced effect of the proposed method is more obvious than of CLAHE and the boundary in blur circle is more reasonable as well, which verifies a better performance of the proposed approach. The shaper edges in Fig. 14(f) than in Fig. 14(e) demonstrate the function of blind de-blurring, i.e., the effectively removed blurred vision. Gain calibration processed profiles (row: 1085:1105; column: 250) of different results (marked in green) are shown in Fig. 13(a), which validates the conspicuous superiority of the propose method coherently.

Fig. 14 Real DR image and its enhanced results: (a) original beetle DR image; enhanced results by using (b) piecewise-linear intensity level transformation; (c) the combination of blind deblurring and piecewise-linear intensity level transformation; (d) CLAHE; (e) the proposed method without blind deblurring process; (f) the proposed method. The correlative histograms of (a)–(f) are illustrated in (g)–(l).

Download Full Size | PDF

Fig. 15 The zoom-in patches of Fig. 14(a)–(f) (inside the red box): (a)–(f) is the first group; (g)–(l) is the second group.

Download Full Size | PDF

6. Conclusion

In this paper, an imaging model is established to describe the image formation in scintillator from the view of the radiative transfer process of visible photons. When solving the corresponding inverse problem, an innovative pre-correction technique for specific intensity of x-rays is developed. In this way, the proposed imaging model coincides with the haze imaging model. By analyzing the feature of DR images caused by OCD-based Micro-CT system, we find the potential dark channel prior which is first proposed in haze removal of outdoor images. Then we extend and apply the haze removal technique, which is based on dark channel prior and guided filter, to deal with DR images. Further, by considering the blurring effect, an effective blind deblurring method is employed to improve the spatial resolution. With the comparison against the conventional contrast enhancement algorithms, our approach performs well in both visual enhancements and physical description. Consequently, our method can improve the quality of DR images effectively.

Acknowledgments

This work was supported by the National Natural Science Foundation of China under the grants 61501310 and 61127003. This work was also supported by the Beijing Center for Mathematics and Information Interdisciplinary Sciences (BCMIIS).

References and links

1. A. D. Maidment and M. J. Yaffe, “Analysis of the spatial-frequency-dependent DQE of optically coupled digital mammography detectors,” Medical physics 21, 721–729 (1994). [CrossRef]

2. M. J. Yaffe and J. Rowlands, “X-ray detectors for digital radiography,” Physics in Medicine and Biology 42, 1 (1997). [CrossRef] [PubMed]

3. E. Kotter and M. Langer, “Digital radiography with large-area flat-panel detectors,” European radiology 12, 2562–2570 (2002). [CrossRef] [PubMed]

4. G. F. Knoll, Radiation Detection and Measurement (John Wiley & Sons, 2010).

5. Y. Zhu, M. Zhao, H. Li, and P. Zhang, “Micro-CT artifacts reduction based on detector random shifting and fast data inpainting,” Medical physics 40, 031114 (2013). [CrossRef] [PubMed]

6. M. Nikl, “Scintillation detectors for X-rays,” Measurement Science and Technology 17, R37 (2006). [CrossRef]

7. J. K. Daugherty, R. C. Hartman, and P. J. Schmidt, “A measurement of cosmic-ray positron and negatron spectra between 50 and 800 MV,” The Astrophysical Journal 198, 493–505 (1975). [CrossRef]

8. U. Bonse, R. Nusshardt, F. Busch, R. Pahl, J. H. Kinney, Q. C. Johnson, R. A. Saroyan, and M. C. Nichols, “X-ray tomographic microscopy of fibre-reinforced materials,” Journal of materials science 26, 4076–4085 (1991). [CrossRef]

9. U. Bonse, F. Busch, O. Günnewig, F. Beckmann, R. Pahl, G. Delling, M. Hahn, and W. Graeff, “3D computed X-ray tomography of human cancellous bone at 8 μm spatial and 10⁻⁴ energy resolution,” Bone and mineral 25, 25–38 (1994). [CrossRef]

10. U. Bonse and F. Busch, “X-ray computed microtomography (μ CT) using synchrotron radiation (SR),” Progress in biophysics and molecular biology 65, 133–169 (1996). [CrossRef]

11. A. Koch, C. Raven, P. Spanne, and A. Snigirev, “X-ray imaging with submicrometer resolution employing transparent luminescent screens,” JOSA A 15, 1940–1951 (1998). [CrossRef]

12. S. H. Williams, A. Hilger, N. Kardjilov, I. Manke, M. Strobl, P. A. Douissard, T. Martin, H. Riesemeier, and J. Banhart, “Detection system for microimaging with neutrons,” Journal of Instrumentation 7, P02014 (2012). [CrossRef]

13. R. C. Gonzalez and R. E. Woods, Digital Imaging Processing (Addison-Wesley, Massachusetts1992).

14. S. M. Pizer, E. P. Amburn, J. D. Austin, R. Cromartie, A. Geselowitz, T. Greer, B. ter Haar Romeny, J. B. Zimmerman, and K. Zuiderveld, “Adaptive histogram equalization and its variations,” Computer vision, graphics, and image processing 39, 355–368 (1987). [CrossRef]

15. K. Zuiderveld, “Contrast limited adaptive histogram equalization,” in “Graphics gems IV,” (Academic Press Professional, Inc., 1994), pp. 474–485. [CrossRef]

16. P. C. Hansen, J. G. Nagy, and D. P. O’leary, Deblurring Images: Matrices, Spectra, and Filtering, vol. 3 (SIAM, 2006). [CrossRef]

17. D. Kundur and D. Hatzinakos, “Blind image deconvolution revisited,” IEEE Signal Processing Magazine , 13, 61–63 (1996). [CrossRef]

18. D. Kundur and D. Hatzinakos, “Blind image deconvolution,” IEEE Signal Processing Magazine , 13, 43–64 (1996). [CrossRef]

19. P. Campisi and K. Egiazarian, Blind Image Deconvolution: Theory and Applications (CRC press, 2007). [CrossRef]

20. N. Wiener, Extrapolation, Interpolation, and Smoothing of Stationary Time Series, vol. 2 (MIT press, 1949).

21. A. N. Tikhonov and V. Y. Arsenin, Solutions of Ill-posed Problems (V.H.Winston and Sons, 1977).

22. L. Rudin, S. Osher, and E. Fatemi, “Nonlinear total variation based noise removal algorithms,” Physica D 60, 259–268 (1992). [CrossRef]

23. L. Rudin and S. Osher, “Total variation based image restoration with free local constraints,” in “Image Processing, 1994. Proceedings. ICIP-94., IEEE International Conference,” (IEEE, 1994), pp. 31–35.

24. K. He, J. Sun, and X. Tang, “Single image haze removal using dark channel prior,” in “Computer Vision and Pattern Recognition, 2009. CVPR 2009. IEEE Conference on,” (Miami, FL, 2009), pp. 1956–1963.

25. K. He, J. Sun, and X. Tang, “Guided Image Filtering,” Pattern Analysis and Machine Intelligence, IEEE Transactions on Pattern Analysis and Machine Intelligence, IEEE Transactions on 35, 1397–1409 (2013). [CrossRef] [PubMed]

26. Q. Wang, Y. Zhu, H. Li, M. Li, P. Yu, and Y. Zhao, “Resolution Improvement of Digital Radiography Image Based on a Modified ROF Deblurring Model,” The 13th International Meeting on Fully Three-Dimensional Image Reconstruction in Radiology and Nuclear Medicine pp. 671–674 (2015).

27. M. Jiang, G. Wang, M. W. Skinner, J. T. Rubinstein, and M. W. Vannier, “Blind deblurring of spiral CT imagesl-comparative studies on edge-to-noise ratios,” Medical physics 29, 821–829 (2002). [CrossRef]

28. M. Jiang and G. Wang, “Development of blind image deconvolution and its applications,” Journal of X-ray Science and Technology 11, 13–19 (2003). [PubMed]

29. M. Jiang, G. Wang, M. W. Skinner, J. T. Rubinstein, and M. W. Vannier, “Blind deblurring of spiral CT images,” Medical Imaging, IEEE Transactions on 22, 837–845 (2003). [CrossRef]

30. I. W. Busbridge, The Mathematics of Radiative Transfer, 50 (University Press, 1960).

31. S. Chandrasekhar, Radiative Transfer (Courier Corporation, 2013).

32. G. C. Pomraning, The Equations of Radiation Hydrodynamics (Courier Corporation, 2005).

33. C. W. Van Eijk, “Inorganic scintillators in medical imaging,” Physics in medicine and biology 47, R85–R106 (2002). [CrossRef]

34. P. Lecoq, A. Annenkov, A. Gektin, M. Korzhik, and C. Pedrini, Inorganic Scintillators for Detector Systems: Physical Principles and Crystal Engineering (Springer, 2006).

35. P. A. Rodnyi, Physical Processes in Inorganic Scintillators, vol. 14(CRC press, 1997).

36. E. Samei and M. J. Flynn, “An experimental comparison of detector performance for direct and indirect digital radiography systems,” Medical physics 30, 608–622 (2003). [CrossRef]

37. A. Ishimaru, Electromagnetic Wave Propagation, Radiation, and Scattering (Prentice-Hall, 1991).

38. H. Koschmieder, Theorie der horizontalen Sichtweite: Kontrast und Sichtweite (Keim & Nemnich, 1925).

39. S. G. Narasimhan and S. K. Nayar, “Chromatic framework for vision in bad weather,” in “Computer Vision and Pattern Recognition, 2000. Proceedings. IEEE Conference on,” vol. 1 (IEEE, 2000), vol. 1, pp. 598–605.

40. R. Fattal, “Single image dehazing,” in “ACM Transactions on Graphics (TOG),” (ACM, 2008), p. 72.

41. R. T. Tan, “Visibility in bad weather from a single image,” in “Computer Vision and Pattern Recognition, 2008. CVPR 2008. IEEE Conference on,” (IEEE, 2008), pp. 1–8.

42. S. G. Narasimhan and S. K. Nayar, “Vision and the atmosphere,” International Journal of Computer Vision 48, 233–254 (2002). [CrossRef]

43. Y. Zhu, M. Zhao, Y. Zhao, H. Li, and P. Zhang, “Noise reduction with low dose CT data based on a modified ROF model,” Opt. Express 20, 17987–18004 (2012). [CrossRef]

44. E. R. Howells, D. C. Phillips, and D. Rogers, “The probability distribution of X-ray intensities. ii. Experimental investigation and the X-ray detection of centres of symmetry,” Acta Crystallographica 3, 210–214 (1950). [CrossRef]

45. N. R. Draper and H. Smith, Applied Regression Analysis2nded. (John Wiley, 1981).

46. T. Hastie, R. Tibshirani, J. Friedman, T. Hastie, J. Friedman, and R. Tibshirani, The Elements of Statistical Learning, vol. 2 (Springer, 2009). 1. [CrossRef]

47. A. Levin, D. Lischinski, and Y. Weiss, “A closed-form solution to natural image matting,” Pattern Analysis and Machine Intelligence, IEEE Transactions on 30, 228–242 (2008). [CrossRef]

48. A. Zomet and S. Peleg, “Multi-sensor super-resolution,” in “Applications of Computer Vision, 2002.(WACV 2002). Proceedings. Sixth IEEE Workshop on,” (IEEE, 2002), pp. 27–31.

49. I. Csiszar, “Why least squares and maximum entropy? An axiomatic approach to inference for linear inverse problems,” The annals of statistics 192032–2066 (1991). [CrossRef]

50. D. Chen, H. Li, Q. Wang, P. Zhang, and Y. Zhu, “Computed tomography for high-speed rotation object,” Optics Express 23, 13423–13442 (2015). [CrossRef] [PubMed]

51. S. Osher, M. Burger, D. Goldfarb, J. Xu, and W. Yin, “An iterative regularization method for total variation-based image restoration,” Multiscale Modeling & Simulation 4, 460–489 (2005). [CrossRef]

52. T. Goldstein and S. Osher, “The split Bregman method for L1-regularized problems,” SIAM Journal on Imaging Sciences 2, 323–343 (2009). [CrossRef]

Symbol	Definition	Symbol	Definition
r⃗	position vector	b_abs(λ, r⃗)	absorption coefficient
v⃗	travel direction	b_sca(λ, r⃗)	scattering coefficient
t	current time	P(λ, r⃗, u⃗, v⃗)	the probability of direction change
λ	wavelength of visible light	S_λ(r⃗, x⃗, t)	photon production rate
I_λ (r⃗, v⃗, t)	photon intensity

Physical Process	Mathematical Model
streaming	−v⃗ ·∇⃗ξ_λ(r⃗, v⃗, t)ΔV
absorption	−b_abs(λ, r⃗) · ξ_λ(r⃗, v⃗, t)ΔV
out-scattering	−b_sca(λ, r⃗) · ξ_λ(r⃗, v⃗, t)ΔV
in-scattering	$b_{sca} (λ, \vec{r}) \int \frac{P (λ, \vec{r}, \vec{u}, \vec{v})}{4 π} \cdot ξ_{λ} (\vec{r}, \vec{u}, \vec{t}) d \vec{u} Δ V$
emission	s_λ(r⃗, v⃗, t)ΔV

Symbol	Definition	Symbol	Definition
r⃗	position vector	b_abs(λ, r⃗)	absorption coefficient
v⃗	travel direction	b_sca(λ, r⃗)	scattering coefficient
t	current time	P(λ, r⃗, u⃗, v⃗)	the probability of direction change
λ	wavelength of visible light	S_λ(r⃗, x⃗, t)	photon production rate
I_λ (r⃗, v⃗, t)	photon intensity

Physical Process	Mathematical Model
streaming	−v⃗ ·∇⃗ξ_λ(r⃗, v⃗, t)ΔV
absorption	−b_abs(λ, r⃗) · ξ_λ(r⃗, v⃗, t)ΔV
out-scattering	−b_sca(λ, r⃗) · ξ_λ(r⃗, v⃗, t)ΔV
in-scattering	$b_{sca} (λ, \vec{r}) \int \frac{P (λ, \vec{r}, \vec{u}, \vec{v})}{4 π} \cdot ξ_{λ} (\vec{r}, \vec{u}, \vec{t}) d \vec{u} Δ V$
emission	s_λ(r⃗, v⃗, t)ΔV

Imaging model for the scintillator and its application to digital radiography image enhancement

Abstract

1. Introduction

2. Primary knowledge of RTE

3. Scintillator imaging model

3.1. Contrast decrease subprocess

3.2. Resolution reduction subprocess

4. Algorithms

4.1. Pre-correction for intensity of x-ray

4.2. DR image enhancement

5. Experimental results

6. Conclusion

Acknowledgments

References and links

Cited By

Figures (15)

Tables (3)

Equations (40)

Optics Express