Initially, we describe the principle of SWIDeX to obtain a projection representing the second-order derivative of the line integral of δ. Figure 1 is a schematic diagram of SWIDeX measurement system using LAA. In this figure, the xy-coordinate system is set for the subject, and the refractive index distribution \(\delta \left(x,y\right)\) of the subject is placed on the xy-plane. The pq-coordinate system is obtained by rotating the xy-coordinate system by an angle φ around the origin. The incident x-ray beam propagates in the xy-plane parallel to the p-axis. The beam propagating through the subject is impinged on the LAA, which is installed parallel to the q-axis. When the beam enters the LAA at an angle of incidence near the Bragg angle \(_\), it is divided at the back of the LAA into a Forward-Diffraction (FD), which propagates in the same direction as the direction of incidence, and a Diffraction (D), which propagates in the \(2_\) direction. The FD and D X-ray intensities can be calculated theoretically from the respective rocking curves, which are functions of the beam incidence angle into the LAA. Each rocking curve has a steep gradient near the Bragg angle, which can be used to produce refraction-contrast images. When the incident X-ray energy is high and the LAA thickness is thin, the sum of the FD and D intensities is approximately equal to the incident X-ray intensity from the energy conservation law, and the FD and D produce contrasts that are inverted from each other. Each X-ray intensity can be defined using the functions \(_\) and \(_\) for each rocking curve as follows:

$$ I_ \left( q \right) = R_ \left( + \alpha \left( q \right)} \right) \cdot I\left( q \right) , $$

(1)

$$ I_ \left( q \right) = R_ \left( + \alpha \left( q \right)} \right) \cdot I\left( q \right) $$

(2)

where \(_\) is the angle of incidence of the X-ray beam into the LAA when the sample is removed and \(I\left(q\right)\) is the X-ray intensity including X-ray absorption in the sample. Although it is possible to separate the absorption components from \(_\left(q\right)\) and \(_\left(q\right)\) by multiple measurements, for simpler derivation, we assume that \(I\left(q\right)=_\) under the assumption that the sample behaves as a pure phase object. \(\alpha \left(q\right)\) is the angular deviation after the beam has propagated through the sample and represents the first derivative of the line integral of δ on the subject:

$$ \alpha \left( q \right) = \int\limits_}^ \right)}}}dp} $$

(3)

A schematic diagram of SWIDeX measurement system

In conventional XDFI, the distance between the LAA and the detector is sufficiently large so that the FD and D images do not overlap, and only the FD signal is measured when little absorption contrast is observed, such as in biological soft tissue. Given that the refraction angle \(\alpha \left(q\right)\) is very small, the observed signal can be described as follows:

$$ \begin I_} \left( q \right) & \cong I_ \left( } \left( } \right) + \frac} }}}\alpha \left( q \right)} \right) \\ \quad & \equiv C\int\limits_}^ \right)}}}dp + D} \\ \end $$

(4)

where \(C=_\frac_}\) and \(D=__\left(_\right)\) are constants and can be estimated from the rocking curve. Thus, XDFI can measure the first derivative of the line integral of δ as in Eq. (4). SWIDeX acquires an image superimposing \(_\left(q\right)\) and \(_\left(q\right)\) by placing the scintillator immediately after the LAA. Defining \(f\left(\theta \right)=_ _\left(\theta \right)\) and \(g\left(\theta \right)=_ _\left(\theta \right)\) in \(\theta \in _\), where \(_\) is the region used for measurement where the rocking curve is monotonic. Equations (1) and (2) are represented as \(_\left(q\right)=f\left(_+\alpha \left(q\right)\right)\), \(_\left(q\right)=g\left(_+\alpha \left(q\right)\right)\), and \(f\left(\theta \right)+g\left(\theta \right)=_\) from the energy conservation law. If \(_\left(q\right)\) and \(_\left(q\right)\) completely overlap, the refraction contrast disappears, but here we consider the situation where the LAA and the scintillator are separated by a small distance \(\Delta p\). Due to the \(\Delta p\) gap, \(_\left(q\right)\) shifts in the q direction by \(\Delta q=\Delta p\text2_\), so the X-ray intensity observed on the scintillator is

$$ I_ \left( q \right) = I_ \left( q \right) + I_ \left( \right) = f\left( + \alpha \left( q \right)} \right) + g\left( + \alpha \left( \right)} \right) = f\left( + \alpha \left( q \right)} \right) - f\left( + \alpha \left( \right)} \right) + I_ $$

(5)

Substituting \(f\left( + \alpha \left( \right)} \right) \cong f\left( + \alpha \left( q \right)} \right) - \frac}\frac}\Delta q\) to Eq. (5),

$$ I_ \left( q \right) \cong \frac}\frac}\Delta q + I_ = C^ \frac} + I_ $$

(6)

where \(C^ = \Delta q\frac}\left( + \alpha \left( q \right)} \right)\), which can be regarded as a constant, i.e., \(\frac}\left( + \alpha \left( q \right)} \right) \approx \frac}\left( } \right)\). This expression is a first-order approximation; a higher-order approximation may be used to represent the output signal with higher precision.

Thus, Eq. (6) is given as

$$ I_} \left( q \right) = C^ \int\limits_}^ \delta \left( \right)}} }}dp + I_ ,} $$

(7)

This means that the contrast of the projection formed on the scintillator represents the second derivative of the line integral of δ. It is also generally known that the Radon transform of a Laplacian image is a second derivative projection:

$$ \mathop \int \limits_^ \Delta \delta \left( \right)dp = \mathop \int \limits_^ \left( }} }} + \frac }} }}} \right)\delta dp = \mathop \int \limits_^ \left( }} }} + \frac }} }}} \right)\delta dp = \mathop \int \limits_^ \frac \delta }} }}dp + \left[ }} \right]_^ = \mathop \int \limits_^ \frac \delta }} }}dp, $$

(8)

After obtaining the projection of Eq. (7) from various directions of the subject using SWIDeX, the Laplacian image with high sparsity is obtained by reconstruction of these projection views.

In conventional XDFI, the projection representing the first derivative is acquired as shown in Eq. (4), so the Laplacian image reconstruction requires prior derivative processing in the q direction. This process is carried out by software differencing, which leads to noise amplification and spatial resolution degradation. With SWIDeX, Eq. (7) can be obtained directly, eliminating the need for extra preprocessing.

The projection reduction by ART + TV method is effective in reconstructing CT images with high sparsity, and it is expected to have a high projection reduction effect when applied to the second-order differential projection of SWIDeX. We describe a method for reconstructing a Laplacian image \(\Delta \delta (x,y)\) of δ using the ART + TV method. First, a digital image \(D(m,n)\) \((1\le m,n\le K)\) of \(\Delta \delta (x,y)\) is introduced. Also, N-dimensional vector representing \(D(m,n)\) is defined as \(\text=(_,_,\cdots ,_)\), where \(N=K\times K\) and \(_=D(m, n)\). Next, Eq. (2) is digitized as follows:

$$ \beginc} ^ w_ d_ = h_ } & \\ \end $$

(9)

where \(_\) corresponds to the right-hand side of Eq. (8). M is the total number of rays, and the weight matrix \(_\) is equal to the line segment where the i-th ray crosses the j-th pixel and represents the contribution of the i-th line integral to the j-th pixel. Equation (9) can be expressed as \( } = } \) where W is an \(M\times N\) matrix with \(_\) as the (i, j)th element and h is an M-dimensional vector with \(_\) as the i-th component. We formulate the reconstruction for D with TV regularization as follows:

$$ \mathop \limits_ ||} - }||^ \, +\, \lambda TV\left( } \right) $$

(10)

$$ TV\left( } \right) = \sum\limits_} \right) - d\left( \right)} \right|^ + \left| \right) - d\left( \right)} \right|^ } } $$

(11)

where \(\lambda \) controls the relative weights of the data fidelity and regularization terms.

To obtain the reconstructed image, the method proposed in [25] is applied. The algorithm obtains a solution from the minimization problem (10) by alternately repeating the ART step and the TV minimization step. In the ART step, the standard ART algorithm reconstructs the image by iteration. The reconstructed image output by the ART step is then supplied to the TV step as the initial value. In the TV step, the objective function in Eq. (11) is minimized using the standard steepest descent method with derivatives of the TV terms and D is the solution of \(\Delta \delta \) when Eq. (10) falls below a certain threshold.

Finally, \(\delta \left(x,y\right)\) is obtained by solving numerically the Poisson equation \(\left(^/\partial ^+^/\partial ^\right)\delta \left(x,y\right)=\widetilde\left(x,y\right)\), where \(\widetilde\) is a Laplacian image, under the Dirichlet boundary condition that \(\delta \left(x,y\right)=0\) for \(\left(x,y\right)\in \partial\Omega \). To solve the Poisson equation, we use the successive over-relaxation (SOR) method [26]. The three-dimensional δ distribution of an object can be obtained by stacking two-dimensional CT reconstruction images.