ART and CS
We use the algebraic reconstruction technique (ART) and CS as a platform for the sparseview image reconstruction. ART is a minimum mean square error (MMSE) solver to find the image f that best matches the measured projection data g:
\mathbf{A}\mathbf{f}=\mathbf{g}
(1)
where A is the system matrix describing the forward projection in the CT scan [21]. In ART, the above equation is solved in an iterative way that the difference between the projection data measured in the real scan and the projection data calculated from the estimated image is backprojected on to the image estimated at the previous iteration step. ART is known to have better performance than FBP in suppressing streak artifacts in sparseview imaging. Many variants of ART with different iteration schemes have been proposed to improve the image quality and to reduce the computation time [22–25]. In this study, we use the orderedsubset simultaneous algebraic reconstruction technique (OSSART) [24] for an ART solver.
The CSbased image reconstruction methods solve the following constrained optimization problem which has constraints of data fidelity and pixel positivity [26]:
\underset{\mathbf{f}}{\text{arg min}}{\left\Psi \mathbf{f}\right}_{{l}_{1}}\phantom{\rule{0.48em}{0ex}}s.t.\phantom{\rule{0.48em}{0ex}}\mathbf{A}\mathbf{f}=\mathbf{g}\text{,}\phantom{\rule{0.48em}{0ex}}\mathbf{f}\ge 0
(2)
where Ψ is a sparsifying transform operator, and {\left\mathbf{z}\right}_{{l}_{1}}={\sum}_{i=1}^{N}\left{z}_{i}\right is the l
_{1} norm of an Ndimensional vector z. In this study, we use the discrete gradient transform for the sparsifying transform which has been widely used in the CSbased image reconstruction [14]:
\Psi \mathbf{f}=\sqrt{{\left[f\left(i+1,j\right)f\left(i,j\right)\right]}^{2}+{\left[f\left(i,j+1\right)f\left(i,j\right)\right]}^{2}}
(3)
where i and j are the pixel indices in the x and ydirections, respectively, and f is the 2D matrix form of f. The discrete gradient transform is often denoted as TV. We implement the CSbased image reconstruction algorithm using the OSSART to enforce the data fidelity and the steepest descent method to minimize the TV in an alternating manner in the iteration. We summarize the algorithm for the CSbased image reconstruction in the following pseudocode [19, 27].
function\phantom{\rule{0.12em}{0ex}}{\text{CS}}_{\text{TV}}\left(\mathbf{g}\text{,}\beta \text{,}{\beta}_{\text{red}}\text{K}\text{,}{\mathbf{f}}^{\text{init}}\right)
(4)

1.
f^{0} : = f^{init};

2.
for k = 1: 1: K(main loop)

3.
update f^{k} by OSSART from the projection data g;

4.
for l = 1: 1: 10 (TV minimization loop)

5.
compute the steepest decent direction d of TV;

6.
\text{\rho}=max\left(\left{\text{f}}^{k}\right\right)\xf7max\left(\left\text{d}\right\right);

7.
{\mathbf{f}}^{k}={\mathbf{f}}^{k}\text{\beta}\times \text{\rho}\times \mathbf{d};

8.
end

9.
\text{\beta}=\text{\beta}\times {\text{\beta}}_{\text{red}};

10.
end

11.
return f^{k}
The TVminimization step has two control parameters, the maximum step size β in the steepest descent search, the reduction factor β_{red} of the maximum step size after each iteration of the main loop. It is commonly known that the large step size of the steepest descent makes the image look smooth, and the small one makes the image look sharp [13, 14]. In this study, we empirically choose β and β_{red} considering that too large β makes the image weakcontrasted whilst too small β makes the image very similar to the one reconstructed by ART [14].
Streakartifactsuppressed CS image reconstruction (SASCS)
To reduce streak artifacts stemming from highintensity structures like bones or metal implants, we combine CSbased image reconstruction approaches with the conventional FBP. Figure 1 shows the basic idea of the proposed method.
Step 1 : {\mathbf{f}}^{\text{FBP}}=\text{FBP}\left({\mathbf{g}}^{\text{acq}}\right)
To identify the highintensity region that makes strong streak artifacts, we first reconstruct an image using FBP, f^{FBP}, from the acquired sinogram g^{acq}. The resulting image may have streak artifacts around the highintensity structures.
Step 2 : Extracting f^{bone} from f^{FBP}
From f^{FBP}
, we extract the highintensity region, denoted as f^{bone}, by applying a thresholding technique. We manually choose the global threshold T^{bone} by visual inspection of the image histogram.
Step 3 : Computing g^{bone} by forward projecting f^{bone}
From f^{bone}, we compute forward projection of the highintensity region to make the sinogram data g^{bone} that accounts for the highintensity region only.
Step 4 : {\mathbf{g}}^{\text{soft}}={\mathbf{g}}^{\text{acq}}{\mathbf{g}}^{\text{bone}}
We subtract g^{bone} from the measured sinogram, g^{acq}, to exclude the components stemming from the highintensity region.
Step 5 : {\mathbf{f}}^{\text{soft}}={\text{CS}}_{\text{TV}}\left({\mathbf{g}}^{\text{soft}},\text{\beta}=0.0060,{\text{\beta}}_{\text{red}}=0.98,\text{K}=30,{\mathbf{f}}^{\text{init}}=0\right)]
We use the subtracted sinogram g^{soft}, which account for the soft tissues only, for reconstruction of soft tissue images using the CSbased method. In this step of CSbased image reconstruction, we use a uniform image of zeroes as an initial guess of the CSbased image reconstruction.
Step 6 : {\mathbf{f}}^{\text{sum}}={\mathbf{f}}^{\text{bone}}+{\mathbf{f}}^{\text{soft}}
After reconstructing the soft tissue image f^{soft} via CS, we add the highintensity region image f^{bone}, which has been reconstructed by FBP, to the soft tissue image to get the composite image f^{sum}.
Step 7 : {\mathbf{f}}^{\text{final}}={\text{CS}}_{\text{TV}}\left({\mathbf{g}}^{\text{acq}},\text{\beta}=0.0033,{\text{\beta}}_{\text{red}}=0.98,\text{K}=30,{\mathbf{f}}^{\text{init}}={\mathbf{f}}^{\text{sum}}\right).
To further refine the CT image, we perform the CSbased iterations again on the original sinogram with the initial guess of the CT image set to f^{sum}obtained at the last step.
The CSbased image reconstruction in step 5 and 7 solves the constrained minimization problem defined in Eq. (2). Step 5 needs two inputs and three control parameters. The two inputs are the soft tissue sinogram, g^{soft}, and the initial guess of the reconstructed image which is all zeroed. The three control parameters are β and β_{red} defined in the previous section, and K the maximum number of iterations of the main loop. In step 7, we perform the CSbased image reconstruction again using the same procedure as in step 5, but with the data inputs of g^{acq} and f^{sum} which has been obtained in step 6.
Data acquisition
We have performed all the CT scans using the labbuilt micro CT system described in our previous work [28]. The micro CT system consists of a microfocus xray source, a rotating object holder, a CMOS flatpanel detector. The microfocus xray source (L812101, Hamamatsu, Japan) has a fixed tungsten anode having an angle of 25° against the electron beam and a 200 μmthick beryllium exit window. The emitted xray beam has a span angle of 43°. The source has a variable focal spot size from 5 μm to 50 μm depending on the applied tube power. We have operated the microfocus xray source in a continuous mode with a 1 mmthick Al filter. We have used a commercially available flatpanel detector (C7942, Hamamatsu, Japan) as a 2D digital xray imager in the microCT system. The flatpanel detector consists of a 2240 × 2240 active matrix of transistors and photodiodes with a pixel pitch of 50 μm, and a CsI:Tl scintillator.
To validate the proposed method, we have performed CT scans of a contrast phantom and a sacrificed adult rat using the microCT. The contrast phantom consists of seven inserts six of which have physical densities similar to that of water and the rest of which has the boneequivalent physical density. Figure 2 shows a schematic diagram of the contrast phantom along with the physical densities of the inserts. The seven inserts of 5 mm diameter were in a water bath made of an acryl cylinder of 40 mm diameter. We have made the inserts using the commercial electron density phantoms (Model 76430, Nuclear Associates, NY, USA). We applied tube voltage and current of 40 kVp and 0.5 mA for the contrast phantom imaging, and 65 kVp and 0.34 mA for the rat imaging, respectively. To get reference images, we performed fullview scans with the number of views of 900 over 360 degrees. To get sparseview projection data, we decimated the fullview projection data in the view direction.
Image quality evaluation
In order to evaluate the final image quality particularly in terms of streak artifact formation, we use two metrics, one for the total variations stemming from streak artifacts and the other for the relative errors of the final image with respect to the reference image. We define the normalized streak indicator (SI) using the total variation of the difference image [12] :
SI=\frac{\text{TV}\left(\mathbf{f}{\mathbf{f}}^{\text{ref}}\right)}{\text{TV}\left({\mathbf{f}}^{\text{FBP}}{\mathbf{f}}^{\text{ref}}\right)}
(5)
where f is the sparseview image reconstructed by the proposed method, f^{ref} the reference fullview image reconstructed by FBP, and f^{FBP} the sparseview image reconstructed by FBP, respectively. We calculate TV using Eq. (3). For the reference images, we use the 900view images reconstucted by FBP which have little streak artifacts.
To evaluate reconstruction errors as compared to the reference image, we use the relative root mean square error (RRME) defined by [12]:
RRME=\sqrt{\frac{\sum _{i,j}{\left(f\left(i,j\right){f}^{\mathit{ref}}\left(i,j\right)\right)}^{2}}{\sum _{i,j}{f}^{\text{ref}}{\left(i,j\right)}^{2}}}\text{,}
(6)
where f is the matrix form of the image vector f.