# Sparse-view Cone Beam CT Reconstruction using Data-consistent Supervised and
Adversarial Learning from Scarce Training Data

Anish Lahiri1, Marc L. Klasky2, Jeffrey A. Fessler1 , and Saiprasad
Ravishankar3 Senior Member, IEEE
1A. Lahiri and J. A. Fessler are with the Department of Electrical and
Computer Engineering, University of Michigan, Ann Arbor, MI 48014 USA. Emails:
<EMAIL_ADDRESS><EMAIL_ADDRESS>
3S. Ravishankar is with the Department of Computational Mathematics, Science
and Engineering, and the Department of Biomedical Engineering, Michigan State
University, East Lansing, MI 48824 USA. Email<EMAIL_ADDRESS>
2M. L. Klaksy is with the Theoretical Division, Los Alamos National
Laboratory, Los Alamos, NM 87545 USA. Email<EMAIL_ADDRESS>
###### Abstract

Reconstruction of CT images from a limited set of projections through an
object is important in several applications ranging from medical imaging to
industrial settings. As the number of available projections decreases,
traditional reconstruction techniques such as the FDK algorithm and model-
based iterative reconstruction methods perform poorly. Recently, data-driven
methods such as deep learning-based reconstruction have garnered a lot of
attention in applications because they yield better performance when enough
training data is available. However, even these methods have their limitations
when there is a scarcity of available training data. This work focuses on
image reconstruction in such settings, i.e., when both the number of available
CT projections and the training data is extremely limited. We adopt a
sequential reconstruction approach over several stages using an adversarially
trained shallow network for ‘destreaking’ followed by a data-consistency
update in each stage. To deal with the challenge of limited data, we use image
subvolumes to train our method, and patch aggregation during testing. To deal
with the computational challenge of learning on 3D datasets for 3D
reconstruction, we use a hybrid 3D-to-2D mapping network for the ‘destreaking’
part. Comparisons to other methods over several test examples indicate that
the proposed method has much potential, when both the number of projections
and available training data are highly limited.

###### Index Terms:

Sparse-views, Computed tomography, Machine learning, Deep learning, Image
reconstruction.

##  I Introduction

Computed Tomography (CT) is an important imaging modality across applications
in medicine, industry, science and security. In this work, we develop an
iterative machine learning-based approach for 3D cone beam CT reconstruction
from very limited measurements or projections, and using limited training
data. In the following, we first review some background in limited-view CT
reconstruction before highlighting the contributions of this work.

###  I-A Background

Cone Beam CT (CBCT) is a CT-based technique that allows for three-dimensional
imaging of an object using X-rays diverging from a source. In CBCT, an entire
3D image volume is reconstructed from a set of 2D projections through the
corresponding object. These projections/measurements are obtained at different
angles or ‘views’ around the object, and are collectively dubbed a sinogram.
There are several approaches for the inverse problem of obtaining an image
from these measurements. A classical method for this task is the analytical
Feldkamp-Davis-Kress (FDK) algorithm [1]. More sophisticated methods for 2D or
3D reconstruction involve model-based reconstruction using iterative
algorithms [2, 3, 4], and data-driven algorithms [5, 6].

Model-based image reconstruction (MBIR) or statistical image reconstruction
(SIR) methods exploit sophisticated models for the physics of the imaging
system and models for sensor and noise statistics as well as for the
underlying object. These methods iteratively optimize for the underlying image
based on the system forward model, measurement statistical model, and assumed
prior for the underlying object [7, 8, 9, 10]. In particular, penalized
weighted least squares (PWLS) approaches have been popular for CT image
reconstruction that optimize a combination of a statistically weighted
quadratic data-fidelity term (capturing the forward and noise model) and a
regularizer penalty that captures prior information of the object [11]. MBIR
methods have often used simple regularizers [12] such as edge-preserving
regularization involving nonquadratic functions of differences between
neighboring pixels [13] (implying image gradients may be sparse) or other
improved regularizers [14, 15, 16, 17].

Within the class of data-driven approaches, dictionary learning [12] and deep
learning based methods for reconstruction have gained popularity in recent
years due to their demonstrated effectiveness in removing artifacts from
images in a variety of modalities, great flexibility and the availability of
curated datasets for training [18, 19, 20].

However, in several applications, acquiring many projections or ‘views’
through the object may be undesirable or impossible. This constraint may be to
reduce exposure to radiation in medical imaging applications, or due to only
pre-set limited or sparse views being possible in industrial or security
applications. Moreover, in dynamic imaging applications, where the object is
changing while being imaged, we would also be limited to fewer views per
temporal state, to prevent blurring. While total variation (TV)-based MBIR
methods have been extensively applied to such sparse-view and sparse angle
reconstruction problems [21, 22, 23, 24, 25, 26], other model-based CT
reconstruction algorithms rely upon learned prior-based regularizers [27, 28,
29]. These methods often require many iterations to converge, leading to large
runtimes, and also require careful selection of regularization parameters to
obtain reasonable image quality trade-offs.

Deep learning-based algorithms have also found considerable use in many
problems ranging from artifact correctionto combination with model-based
reconstruction [30, 31, 32, 33]. Deep learning approaches could be supervised
or unsupervised or mixed [34, 35] and include image-domain (denoising)
methods, sensor-domain methods, AUTOMAP, as well as hybrid-domain methods (cf.
reviews in [12, 36]). Hybrid-domain methods are gaining increasing interest
and enforce data consistency (i.e., the reconstruction should be consistent
with the measurement model) during training and reconstruction to improve
stability and performance. Deep learning methods often require large training
data sets and long training times to work well. They may also struggle to
generalize to data with novel features or that are obtained with different
experimental settings.

When reconstructing 3D objects from extremely limited tomographic views or
projections, many of the aforementioned approaches fail. The FDK algorithm
yields reconstructions that are severely ridden with streak artifacts. While
conventional iterative methods perform better, the quality of reconstructions
leave a lot to be desired, and there is often poor bias-variance trade-offs.
Deep learning-based approaches have the potential to perform better in this
scenario, but still perform poorly when there is a scarcity of available data
for training, such as in national security applications where experimental
data is limited and accurate simulations are expensive [37]. While there are
approaches that reconstruct from very limited projections, they either do not
target 3D CBCT imaging [19, 25, 38, 39], or rely upon many paired training
image volumes [40, 18, 20].

###  I-B Contributions

This paper focuses on developing a method that can provide quality
reconstructions from extremely limited projections even with very extremely
limited training data. The proposed reconstruction approach works across
multiple stages, similar to an unrolled-loop algorithm [12], where each stage
consists of a shallow CNN block trained using a combined supervised and
adversarial loss, followed by a data-consistency block. The adversarial
component of our loss yields destreaked images that have more realistic
texture. To mitigate the challenge of reduced training data, we reduce the
scope of our learning to patches or image sub-volumes. This approach allows us
to provide several training examples from even a single training image.
Destreaked patches are aggregated before data-consistency is applied to the
whole volume. Furthermore, we prime our method using an edge-preserving
regularized reconstruction as input.

We compare our methods to a variety of techniques including the FDK algorithm,
edge-preserving regularized reconstruction, and deep CNN-based reconstruction
without data-consistency. Simulation results suggest that the proposed method
provides much better image quality than previous techniques with extremely
limited (four or eight) views of 3D objects.

###  I-C Organization

The rest of this paper is organized as follows. Section II describes the
proposed approach in detail. Section III explains our choices for various
algorithm parameters as well as our experimental setup. Section IV presents
the results of our comparisons to other algorithms as well as other
experiments that offer insights into the process of our reconstruction.
Section V elaborates upon these observations. Finally, Section VI states our
conclusions and offers some avenues for future research.

##  II Algorithm & Problem Setup

Figure 1: Flow diagram depicting the overall pipeline of our algorithm, where
$\bm{x}_{\text{FDK}}$ is the FDK reconstruction, $\bm{x}_{\text{EP}}$ is an
edge-preserving regularized reconstruction, $\bm{x}_{\text{G}}$ is the
generator’s output after slice aggregation, and $\bm{x}_{1}$ is the output of
the first stage.

Our proposed method for CBCT reconstruction focuses on addressing two primary
challenges: a very limited number of available views, and limited number of
available training objects. We address the former through a combination of
three aspects: (1) using an edge-preserving regularized reconstruction [4] to
initialize our iterative-type algorithm; (2) including an adversarial
component to the training loss function for our learned destreaking networks
(similar to generative adversarial networks or GANs); and (3) including data-
consistency blocks that reinforce acquired measurements in the destreaked 3D
reconstruction.

The problem of scarce training data is addressed primarily by two approaches.
First, we split an entire image volume into patches in the form of overlapping
subvolumes. Essentially, this step localizes the scope of CNN-based
destreaking to a comparatively smaller neighborhood, while allowing us to
generate many training examples from a single image volume. Second, we use a
shallow destreaking CNN to avoid overfitting to the training data. To reduce
computation time associated with multiple 3D convolutions and subsequent patch
aggregation, the CNN is designed to map 3D subvolumes to 2D slices [41]. This
approach enables using 3D contextual information for the destreaking task,
while removing the need for patch averaging (of overlapping 3D patches) and
associated artifacts during aggregation.

Many learning-based reconstruction approaches in applications like MRI work
through end-to-end training, whereas such approaches are less practicle in CT
due to the complexity of the system matrix. Thus, our developed algorithm
operates as a multi-stage greedy approach similar to works like [42, 43]. Each
stage is composed of a CNN that maps each 3D subvolume with streaking
artifacts to a clean 2D slice corresponding to the slice at the centre of the
subvolume. Because the objects considered here have finite support, we treat
the slices at the edge of the volume in the direction of aggregation are being
empty, and set them to zero. One could use other boundary conditions for long
objects [44]. Once an entire image volume has been aggregated from individual
clean slices, this volume is passed through a data-consistency update to
reinforce acquired measurements and reduce any ‘hallucinations’ introduced by
the network. This output subvolume is then provided as input to the next
stage. Fig. 1 depicts the process that is akin to algorithm unrolling [45]. As
mentioned earlier, to reduce noise and streaks, the input to the first stage
of our method is an edge-preserving regularized (iteratively obtained)
reconstruction, using the regularizer in [4] and the algorithm in [46], which
in turn is initialized with an FDK reconstruction for faster convergence.

We train the CNN parameters separately for each stage. The training loss for
the destreaking CNN in each stage consists of a weighted combination of a
masked mean squared error term calculated over a region of interest using the
ground truth training image slices, and an adversarial input from another CNN
that acts as a discriminator for the output of the destreaking CNN (also
specific to the stage). Adversarial training often is posed as a min-max
optimization problem [47], but the practical implementation involves
alternating between updating generator (destreaking) network $G$ parameters
$\phi^{k}$ and discriminator network $D$ parameters $\theta^{k}$, with more
frequent updates of the generator parameters. Our approach to updating the
generator weights for the $k$th stage is mathematically expressed as:

$\displaystyle\begin{split}\hat{\theta}^{k}=\arg\min_{\theta}&\
{-\lambda}\,\mathbb{E}\big{[}D_{\phi^{k}}(G_{\theta}(\mathcal{P}_{3}\bm{x}_{k-1}))\big{]}\\\
&+\mathbb{E}\big{[}\|\mathcal{P}_{3}^{2,\text{mid}}\bm{x}_{\text{GT}}-G_{\theta}(\mathcal{P}_{3}\bm{x}_{k-1})\|_{2}^{2}\big{]},\end{split}$
(1)

where $\bm{x}_{k-1}$ is the output of the $(k-1)$th stage of our algorithm,
$\bm{x}_{0}$ is set to be $\bm{x}_{\text{EP}}$ (the edge-preserving
regularized reconstruction), $\bm{x}_{\text{GT}}$ is the ground truth,
$\lambda$ is a regularization parameter that varies as the weights of $G^{k}$
are updated (see Section III-B), $\mathcal{P}_{3}$ is a 3D patch or subvolume
extraction operator, and $\mathcal{P}_{3}^{2,\text{mid}}$ is an operator that
extracts the 2D central slice from an image subvolume, where the position of
the subvolume is determined by $\mathcal{P}_{3}$. We restricted the
reconstructions to a region in the image volume containing the object of
interest. The expectation $\mathbb{E}$ is taken over the set of training
examples.

Our approach to updating the the discriminator network parameters is likewise
given as:

$\displaystyle\begin{split}\hat{\phi}^{k}=\arg\min_{\phi}~{}&\mathbb{E}\big{[}\big{(}D_{\phi}(G_{\theta^{k}}(\mathcal{P}_{3}\bm{x}_{k-1}))-0\big{)}^{2}\big{]}\\\
&+\mathbb{E}\big{[}\big{(}D_{\phi}(\mathcal{P}_{3}^{2,\text{mid}}\bm{x}_{\text{GT}})-1\big{)}^{2}\big{]}.\end{split}$
(2)

The data-consistency update involves seeking an image that is consistent with
the acquired measurements while still being ‘close’ to the slice-aggregated
destreaked image. The optimization problem for this step is framed as:

$\displaystyle\bm{x}_{k}=\underset{\bm{x}}{\textnormal{arg
min}}\|\bm{A}\bm{x}-\bm{y}\|_{2}^{2}+\beta\|\bm{x}-\bm{x}_{k,\text{G}}\|_{2}^{2},$
(3)

where $\bm{A}$ is the CBCT system matrix, implemented with the separable-
footprint projector [48], $\bm{y}$ denotes the projections or acquired
measurements, $\beta>0$ is a regularization parameter, and
$\bm{x}_{k,\text{G}}$ is the output of the generator after slice aggregation
at the $k$th stage. We used an ordinary least-squares (LS) data-fit term
rather than a weighted LS (WLS) term because the focus here is on sparse views
rather than low-dose imaging, but the method generalizes directly to the WLS
case. We used 50 conjugate gradient (CG) iterations to minimize (3).

##  III Methods

###  III-A Dataset and Experimental Settings

To train and test our method, we used the publicly available 3D walnut CT
dataset [49]. To study the ability to learn from very limited data, we used a
single walnut for training our method, and tested our algorithm on 5 different
walnuts. Furthermore, extremely limited data with 8 or 4 views/projections
through the walnuts were used in training and testing the network. Separate
networks were trained for reconstructing image volumes from 4 and 8 views,
respectively. These CBCT views were generated using the MIRT [50] package, and
were equally spaced over 360 degrees. The distance from source to detector was
set to be 20 cm, and the distance from the object to the detector was 4.08 cm,
and each projection view was $150\times 150$ pixels of size $\approx 0.4$ mm
square. Because the CBCT system simulated here has a small cone angle (almost
parallel beam), 8 views over 360∘ probably has only a bit more information
than 4 views over $180^{\circ}$. The image volume for each walnut was
$501\times 501\times 501$. The dimensions of each voxel were approx.
$0.12\times 0.12\times 0.12$ mm3.

###  III-B Hyperparameters and Network Architectures

The subvolume size for our scheme was chosen to be $500\times 500\times 8$.
The parameter $\beta$ was chosen to be 1, and $\lambda$ was changed
dynamically as the network weights were updated according to
$10^{\lfloor\log(r)\rfloor}$, where
$r=\mathbb{E}\big{[}\|\mathcal{P}_{3}^{2,\text{mid}}\bm{x}_{\text{GT}}-G_{\theta}(\mathcal{P}_{3}\bm{x}_{k-1})\|_{2}^{2}\big{]}$.
This was done to maintain balance between the $\ell_{2}$ and adversarial loss
components during generator training. The number of stages in our method was
set to 4, and the networks in each stage were trained for 40 epochs. The
weights for the discriminator were updated once for every 10 times the weights
of the destreaking CNN were updated.

Fig. 1 depicts the generator architecture we used. The kernel size for 3D
convolutions was ($3\times 3\times 3$), and ($3\times 3$) for the 2D
convolutions. The discriminator was akin to a classifier with two
convolutional layers and with 8 filters in each convolutional layer with
($3\times 3$) kernel size with stride 1 followed by fully connected layers
with (1152,8,8) nodes respectively, with a sigmoid activation at the final
output to constrain the output to be between 0 and 1. The batch size during
training the destreaking CNNs was 6. Training the destreaking CNN in each
stage of our algorithm took approx. 10 hours on 3 NVIDIA Quadro RTX 5000 GPUs,
while at test time, each walnut volume required 7 minutes to reconstruct with
a batch size of 3 on two of the same GPUs. The data consistency update
required an additional 3 minutes on a workstation with Intel(R) Xeon(R) Silver
4214 CPU @ 2.20GHz with 48 cores.

###  III-C Compared Methods

To assess the performance of our method, we used 4 stages of our proposed
algorithm to compromise between image quality and runtime. We compared the
output for all 5 test walnuts to the conventional FDK reconstruction, an edge-
preserving regularized (MBIR) reconstruction [4], as well as the slice-
aggregated output from a single stage of our destreaking CNN without data
consistency.

###  III-D Performance Metrics

We primarily used the normalized mean absolute error (NMAE) as a metric for
evaluating the performance of various methods. The error is evaluated over the
voxels within the region-of-interest (ROI) of a three-dimensional mask
obtained by dilating a ground truth segmentation of the walnut being
reconstructed. The masked region includes all voxels within the shell of the
walnut. The NMAE normalization used the mean intensity of the ground truth
voxels within this mask. Essentially,
$E_{\text{NMAE}}(x_{GT},x_{O},\mathscr{M})=\|\mathscr{M}\odot(x_{GT}-x_{O})\|_{1}/\|\mathscr{M}\odot
x_{GT}\|_{1}$, where $x_{GT}$ is the ground truth image volume, $x_{O}$ is the
reconstruction whose quality is being evaluated, and $\mathscr{M}$ is a binary
mask specific to the test walnut, which excludes any pixel not within the a
dilation of the outer boundary of the walnut. These are obtained by a
histogram-based thresholding of the corresponding ground truth volumes for the
test walnuts.

Another metric that is used for comparison in our work is the normalized high-
frequency error norm or NHFEN [51]. We computed the HFEN for every slice of
the reconstructed walnut as the $\ell_{2}$ norm of the difference of masked
edges (obtained through a high-pass filtering) between the input and reference
images. The masking is done similarly as described earlier. A Laplacian of
Gaussian (LoG) filter was used as the edge detector. The kernel size was set
to $15\times 15$, with a standard deviation of 1.5 pixels. The normalization
was performed over the high frequency components of the ground truth image
over the masked ROI. Mathematically, this metric is calculated as
$E_{\text{NHFEN}}(x_{GT},x_{O},\mathscr{M})=\frac{1}{N}\sum_{i}\|\mathscr{H}(\mathscr{M}[:,:,i]\odot
x_{GT}[:,:,i])-(\mathscr{H}(\mathscr{M}[:,:,i]\odot
x_{O}[:,:,i])\|_{2}/\|\mathscr{H}(\mathscr{M}[:,:,i]\odot
x_{GT}[:,:,i])\|_{2}$, where $\mathscr{H}$ denotes the LoG filter described
earlier, $i$ indexes the slices of the image volume in the $z$ direction,
where $N$ is the total number of slices in that direction, and the other
symbols have their usual meaning, as described previously. An advantage of
using such normalized metrics is that it allows for the evaluation of the
reconstruction quality only in areas of interest in the volume, disregarding
the effect of empty spaces around it.

##  IV Results

Table I compares the reconstruction performance of various methods (including
our own) described in the previous sections. The proposed approach
substantially improves the NMAE and NHFEN compared to the reference methods
for both 8 and 4 acquired projections. As expected, the quality of
reconstructions using 4 acquired projections was worse than when 8 projections
were acquired for reconstruction.

8 views
---
Walnut # | FDK recon. | EP recon. | CNN destreaking | Proposed
| NMAE | NHFEN | NMAE | NHFEN | NMAE | NHFEN | NMAE | NHFEN
1 | 0.77 | 0.90 | 0.45 | 0.58 | 0.40 | 0.59 | 0.26 | 0.54
2 | 0.77 | 0.88 | 0.45 | 0.57 | 0.38 | 0.58 | 0.25 | 0.53
3 | 0.79 | 0.90 | 0.49 | 0.61 | 0.42 | 0.62 | 0.30 | 0.58
4 | 0.79 | 0.96 | 0.45 | 0.62 | 0.39 | 0.62 | 0.27 | 0.58
5 | 0.82 | 0.98 | 0.48 | 0.65 | 0.41 | 0.66 | 0.27 | 0.61
4 views
Walnut # | FDK recon. | EP recon. | CNN destreaking | Proposed
| NMAE | NHFEN | NMAE | NHFEN | NMAE | NHFEN | NMAE | NHFEN
1 | 1.08 | 1.11 | 0.62 | 0.61 | 0.61 | 0.61 | 0.44 | 0.60
2 | 1.23 | 1.18 | 0.65 | 0.60 | 0.63 | 0.60 | 0.45 | 0.60
3 | 1.22 | 1.15 | 0.67 | 0.64 | 0.66 | 0.64 | 0.48 | 0.64
4 | 1.22 | 1.22 | 0.68 | 0.64 | 0.67 | 0.65 | 0.46 | 0.65
5 | 1.23 | 1.27 | 0.65 | 0.68 | 0.65 | 0.68 | 0.46 | 0.68
TABLE I: Comparison of the performance of our proposed method against FDK
reconstruction, edge-preserving (EP) regularized reconstruction, a single
stage of CNN-based destreaking without data-consistency (using the same
architecture as in our method) and our proposed multistage reconstruction
algorithm for 8 (top) and 4 (bottom) acquired projections. The metrics used
for performance are the normalized mean absolute error (NMAE) and the
normalized high-frequency error norm (NHFEN).

Fig. 2 gives an example reconstruction for 8 views (walnut 1 in Table I),
showing central slices through the walnut in all three orientations (sagittal,
coronal and transverse). The proposed algorithm provides significantly higher
quality reconstructions than the other methods. This is particularly evident
in the extent to which our algorithm is able to restore the finer features of
the walnuts, and has fewer artifacts.

Ground Truth Test Volume | Ground Truth Training Volume
---|---
|
(a) (NMAE) | (b) (n/a)
FDK | EP Regularized Recon.
|
(c) (0.77) |  (d) (0.45)
Destreaking CNN | Proposed
|
(e) (0.40) | (f) (0.26)
Figure 2: Comparison of the quality of reconstruction of our proposed
algorithm (f) for walnut 1 (8 views) in Table I to various reference methods.
Each subfigure depicts slices through the center of the walnut volume in the
sagittal, coronal and transverse orientations. The normalized mean absolute
errors have also been shown underneath each subfigure. The central slices
corresponding to the ground truth training walnut volume have also been shown
in (b).

Fig. 3 compares the reconstructions from 4 views. The quality of the
reconstruction is poorer compared to that using 8 views, though the proposed
approach still visibly outperforms the other methods.

Ground Truth Test Volume |
---|---
|
(a) (NMAE) |
FDK | EP Regularized Recon.
|
(b) (1.23) |  (c) (0.65)
Destreaking CNN | Proposed
|
(d) (0.63) | (e) (0.45)
Figure 3: Comparison of the quality of reconstruction of our proposed
algorithm (e) for walnut 2 (4 views) in Table I to various reference methods.
Each subfigure depicts slices through the center of the walnut volume in the
sagittal, coronal and transverse orientations. The normalized mean absolute
errors have also been shown underneath each subfigure.

Fig. 4 shows the results of the intermediate steps of the first stage of our
reconstruction for walnut 2 (8 views) in the test dataset. While the patch-
based destreaking compensates for the blurring introduced by the EP-
regularized reconstruction and ‘fills-in’ details in the reconstructed volume,
the data-consistency plays a key role in mitigating hallucinations introduced
by the CNN, and reinforces image features that are consistent with the
acquired measurements.

Input EP Recon. | Post CNN-based Destreaking | Post Data-consistency
---|---|---
|  |
(NMAE:0.45) | (NMAE:0.38) | (NMAE:0.32)
Figure 4: Central slices corresponding to all three orientations of test
walnut 2 (8 views) after various steps during a single stage (Stage 1) of our
proposed algorithm.

Fig. 5 shows walnut 1 from our test dataset being progressively reconstructed
from 8 projections across the stages of our algorithm; as the stages progress,
more features are restored in the reconstructed walnut, until the improvements
become incremental. The residual streaking artifacts outside the walnut are
mitigated in the reconstructions from the third and fourth stages.

Stage 1 | Stage 2
---|---
|
(a) (MAE: 0.32) |  (b) (MAE: 0.29)
Stage 3 | Stage 4
|
(c) (MAE: 0.27) | (d) (MAE: 0.26)
Figure 5: Central slices through the reconstructions of walnut 1 (8 views) in
all three orientations at the end of each of the 4 stages of our algorithm
(post data-consistency). It is evident that the quality of the reconstruction
progressively improves across the stages. The normalized mean absolute error
(NMAE) for each figure is also provided underneath

Fig. 6 (a) and (b) show how the quality of the reconstructions from 8 acquired
views using our method varies when the orientation of the test walnuts is
changed, in comparison to training-time. This is achieved by changing the
position of the acquired (equidistant) projections, as this is akin to
rotating the test walnut. For this purpose, essentially the position of the
first acquired projection is shifted by a specified angle. For the three test
walnuts shown in the figure, the angle of rotation was changed between
+22.5°and -22.5°in intervals of 7.5°. Both the normalized mean absolute error
(NMAE) and the normalized high-frequency error norm (NHFEN) were used as
quality metrics in this experiment. Our method seems to be fairly robust to
rotated test-data.

|
---|---
(a) | (b)
|
(c) | (d)
Figure 6: (a) and (b) show the effects of rotating the test walnuts (by
changing the angle of acquired projections) on the quality of the
reconstructions for 8 acquired views, while (c) and (d) capture the effect of
varying the scale of features of the test walnut compared to the scale used in
training. (a) and (c) depict the normalized MAE metric while (b) and (d)
depict the normalized HFEN error metric for respective experiments.

We also studied the effect of varying the scale of the features of test
walnuts on reconstructions from our method (from 8 projections), which was
still trained on a single scale. For this purpose, we used the three walnuts
which were used in the previous experiment. From Fig. 6 (c) and (d), it is
evident that while there are differences in the NMAE and NHFEN across scales,
our method holds up reasonably well across multi-scale test data.
Particularly, we notice that while higher-frequency features are better
reconstructed when the test walnuts are scaled down, the best NMAE is observed
when the scale of the test walnut is 1.1 times its original scale. We surmise
that this happens because at a scale of 1.1, the scale of features of the test
walnuts matches that of the training walnut very well. (We excluded scale 1.3
for walnut 3 because at that scale, the walnut exceeded the $501\times
501\times 501$ voxel grid.)

##  V Discussion

Our observations in Section IV indicate that the proposed physics-aware
learning-based approach for limited-view CBCT reconstruction is able to
improve upon the quality of reconstructions yielded not only by the FDK
algorithm, but also traditional prior-based iterative reconstruction, which
serves a crucial role in initializing our algorithm, and ensures that the
subsequent CNN-based destreaking is afforded a reasonable mapping to learn
when trained with ground truth images as targets. Furthermore, in Fig. 4, we
also demonstrated the importance of data-consistency when reconstructing from
limited measurements. The data consistency update could easily be the most
crucial step in our algorithm because it corrects for hallucinations
introduced by the CNN-based destreaking step, which are very likely to occur
given the limited availability of training data. This allows our algorithm to
be repeated for several stages, which is key to the improvement provided by
the proposed method over a simple image domain CNN-based denoiser. The
generalizability of our approach is evident from the consistent improvements
yielded by our approach in Table I, which is a consequence of using subvolumes
or patches in training our method as well as the relative shallowness of the
CNNs used therein. These measures effectively reduce the chances of
overfitting when training on extremely limited full-view data. The chosen
patch size effectively allowed the learned destreaking CNN to utilize three-
dimensional context while also reducing: (1) the total number of overlapping
subvolumes that are forward propagated through the destreaking CNN at test
time and (2) the memory requirements for training the network. The hybrid 3D
to 2D mapping also plays an important role by reducing the computational and
time demands of our algorithm, and by removing any chances of artifacts that
may be introduced during (conventional overlapping) patch aggregation.

We noted that the final reconstruction quality for some test walnuts was
better than others, and also that the final reconstruction quality from our
algorithm depended on the quality of the initial EP reconstruction. We think
this may be because the acquired views (equidistant over $360^{\circ}$) may
not be the best choice for higher fidelity reconstructions for some walnuts,
and because these walnuts may be less similar to the training walnut.

Regarding the performance of our reconstructions in terms of the NHFEN metric,
we observe a lack of sensitivity with respect to other compared techniques for
both 8 and 4 views. We surmise that this is because neither of the networks
was trained for improved NHFEN, i.e., a NHFEN term was not part of the
training loss.

While the experiments with rotated projections and varying scale of features
showcased some robustness of our proposed approach to both rotation and
up/downscaling of features– it also brought to light the necessity for data
augmentation for improved performance. The consistent degradation of the high
frequency features (observed through the NHFEN metric) in reconstructions with
increase in scale was interesting to note, and its contrast with the trends in
NMAE points towards a sharpness/fidelity trade-off that may need more
investigation.

However, it is key to realize that the improved performance in the NHFEN
metric at smaller scales using our method is likely due to fewer edges that
need to be reproduced in the image, and does not suggest that reconstructions
of increasingly sharper quality may be obtained by continually shrinking the
scale of walnuts.

##  VI Conclusions and Future Work

This paper developed a method to provide high quality reconstructions from
extremely limited CBCT projections and scarce training data. The key features
of our approach were the multi-stage approach of alternating between learning-
based destreaking and data consistency and the use of subvolume-based learning
and shallower (adversarially trained) CNNs to combat over-fitting. In the
future, we will focus on extending our method to dynamic imaging applications
in CT, as well as being able to jointly segment and reconstruct three-
dimensional objects. Towards this end, we are also interested in finding a
better metric than the mean absolute error to assess the fidelity and quality
of our reconstructions for specific tasks.

The code for reproducing the results in this paper will be available on github
after the paper is accepted.

## References

  * [1] L. A. Feldkamp, L. C. Davis, and J. W. Kress, “Practical cone beam algorithm,” _J. Opt. Soc. Am. A_ , vol. 1, no. 6, pp. 612–9, Jun. 1984\.
  * [2] J.-B. Thibault, K. Sauer, C. Bouman, and J. Hsieh, “A three-dimensional statistical approach to improved image quality for multi-slice helical CT,” _Med. Phys._ , vol. 34, no. 11, pp. 4526–44, Nov. 2007.
  * [3] J. Xu and B. M. W. Tsui, “Interior and sparse-view image reconstruction using a mixed region and voxel-based ML-EM algorithm,” _IEEE Trans. Nuc. Sci._ , vol. 59, no. 5, pp. 1997–2007, Oct. 2012.
  * [4] J. H. Cho and J. A. Fessler, “Regularization designs for uniform spatial resolution and noise properties in statistical image reconstruction for 3D X-ray CT,” _IEEE Trans. Med. Imag._ , vol. 34, no. 2, pp. 678–89, Feb. 2015.
  * [5] G. Wang, “A perspective on deep imaging,” _IEEE Access_ , vol. 4, pp. 8914–24, Nov. 2016.
  * [6] I. Y. Chun, X. Zheng, Y. Long, and J. A. Fessler, “Sparse-view X-ray CT reconstruction using $\ell_{1}$ regularization with learned sparsifying transform,” in _Proc. Intl. Mtg. on Fully 3D Image Recon. in Rad. and Nuc. Med_ , 2017, pp. 115–9.
  * [7] I. A. Elbakri and J. A. Fessler, “Statistical image reconstruction for polyenergetic X-ray computed tomography,” _IEEE Trans. Med. Imag._ , vol. 21, no. 2, pp. 89–99, Feb. 2002.
  * [8] K. Sauer and C. Bouman, “A local update strategy for iterative reconstruction from projections,” _IEEE Trans. Sig. Proc._ , vol. 41, no. 2, pp. 534–48, Feb. 1993.
  * [9] J. A. Fessler, “Statistical image reconstruction methods for transmission tomography,” in _Handbook of Medical Imaging, Volume 2. Medical Image Processing and Analysis_ , M. Sonka and J. M. Fitzpatrick, Eds. Bellingham: SPIE, 2000, pp. 1–70.
  * [10] J.-B. Thibault, C. A. Bouman, K. D. Sauer, and J. Hsieh, “A recursive filter for noise reduction in statistical iterative tomographic imaging,” in _Proc. SPIE 6065 Computational Imaging IV_ , 2006, p. 60650X.
  * [11] A. R. De Pierro, “On the relation between the ISRA and the EM algorithm for positron emission tomography,” _IEEE Trans. Med. Imag._ , vol. 12, no. 2, pp. 328–33, Jun. 1993.
  * [12] S. Ravishankar, J. C. Ye, and J. A. Fessler, “Image reconstruction: from sparsity to data-adaptive methods and machine learning,” _Proc. IEEE_ , vol. 108, no. 1, pp. 86–109, Jan. 2020.
  * [13] S. Ahn, S. G. Ross, E. Asma, J. Miao, X. Jin, L. Cheng, S. D. Wollenweber, and R. M. Manjeshwar, “Quantitative comparison of OSEM and penalized likelihood image reconstruction using relative difference penalties for clinical PET,” _Phys. Med. Biol._ , vol. 60, no. 15, pp. 5733–52, Aug. 2015.
  * [14] W. Yu, C. Wang, X. Nie, M. Huang, and L. Wu, “Image reconstruction for few-view computed tomography based on $\ell_{0}$ sparse regularization,” _Procedia Computer Science_ , vol. 107, pp. 808–813, 2017, advances in Information and Communication Technology: Proceedings of 7th International Congress of Information and Communication Technology (ICICT2017).
  * [15] C. Xu, B. Yang, F. Guo, W. Zheng, and P. Poignet, “Sparse-view CBCT reconstruction via weighted Schatten p-norm minimization,” _Opt. Express_ , vol. 28, no. 24, pp. 35 469–35 482, Nov 2020.
  * [16] X. Zhang, Y. Zhou, W. Zhang, J. Sun, and J. Zhao, “Adaptive Prior Patch Size Based Sparse-View CT Reconstruction Algorithm,” in _2020 IEEE 17th International Symposium on Biomedical Imaging (ISBI)_ , 2020, pp. 1–4.
  * [17] Y. Kim, H. Kudo, K. Chigita, and S. Lian, “Image reconstruction in sparse-view CT using improved nonlocal total variation regularization,” in _Developments in X-Ray Tomography XII_ , vol. 11113, International Society for Optics and Photonics. SPIE, 2019, pp. 329–337.
  * [18] R. Anirudh, H. Kim, J. J. Thiagarajan, K. A. Mohan, K. Champley, and T. Bremer, “Lose the views: limited angle CT reconstruction via implicit sinogram completion,” in _Proc. IEEE Conf. on Comp. Vision and Pattern Recognition_ , 2017, pp. 6343–52.
  * [19] Y. Han, J. Kang, and J. C. Ye, “Deep learning reconstruction for 9-view dual energy CT baggage scanner,” in _Proc. 5th Intl. Mtg. on Image Formation in X-ray CT_ , 2018, pp. 407–10.
  * [20] H. Kim, R. Anirudh, K. A. Mohan, and K. Champley, “Extreme few-view ct reconstruction using deep inference,” _arXiv preprint arXiv:1910.05375_ , 2019.
  * [21] E. Y. Sidky, K.-M. Kao, and X. Pan, “Accurate image reconstruction from few-views and limited-angle data in divergent-beam CT,” _J. X-Ray Sci. and Technol._ , vol. 14, no. 2, pp. 119–39, 2006. [Online]. Available: http://iospress.metapress.com/content/1jduv1cll3f9e2br/(notfree)
  * [22] H. Yu and G. Wang, “Compressed sensing based interior tomography,” _Phys. Med. Biol._ , vol. 54, no. 9, pp. 2791–806, May 2009.
  * [23] J. Bian, J. H. Siewerdsen, X. Han, E. Y. Sidky, J. L. Prince, C. A. Pelizzari, and X. Pan, “Evaluation of sparse-view reconstruction from flat-panel-detector cone-beam CT,” _Phys. Med. Biol._ , vol. 55, no. 22, p. 6575, Nov. 2010.
  * [24] S. Ramani and J. A. Fessler, “A splitting-based iterative algorithm for accelerated statistical X-ray CT reconstruction,” _IEEE Trans. Med. Imag._ , vol. 31, no. 3, pp. 677–88, Mar. 2012.
  * [25] G. T. Herman and R. Davidi, “Image reconstruction from a small number of projections,” _Inverse Prob._ , vol. 24, no. 4, p. 045011, Aug. 2008.
  * [26] G.-H. Chen, J. Tang, and S. Leng, “Prior image constrained compressed sensing (PICCS): A method to accurately reconstruct dynamic CT images from highly undersampled projection data sets,” _Med. Phys._ , vol. 35, no. 2, pp. 660–3, Feb. 2008.
  * [27] L. Pfister and Y. Bresler, “Model-based iterative tomographic reconstruction with adaptive sparsifying transforms,” in _Proc. SPIE 9020 Computational Imaging XII_ , 2014, p. 90200H.
  * [28] C. Zhang, T. Zhang, M. Li, C. Peng, Z. Liu, and J. Zheng, “Low-dose CT reconstruction via L1 dictionary learning regularization using iteratively reweighted least-squares,” _BioMedical Engineering OnLine_ , vol. 15, no. 66, 2016.
  * [29] X. Zheng, S. Ravishankar, Y. Long, and J. A. Fessler, “PWLS-ULTRA: An efficient clustering and learning-based approach for low-dose 3D CT image reconstruction,” _IEEE Trans. Med. Imag._ , vol. 37, no. 6, pp. 1498–510, Jun. 2018.
  * [30] H. Chen, Y. Zhang, Y. Chen, J. Zhang, W. Zhang, H. Sun, Y. Lv, P. Liao, J. Zhou, and G. Wang, “LEARN: Learned experts: assessment-based reconstruction network for sparse-data CT,” _IEEE Trans. Med. Imag._ , vol. 37, no. 6, pp. 1333–47, Jun. 2018.
  * [31] D. Wu, K. Kim, G. E. Fakhri, and Q. Li, “Iterative low-dose CT reconstruction with priors trained by artificial neural network,” _IEEE Trans. Med. Imag._ , vol. 36, no. 12, pp. 2479–86, Dec. 2017.
  * [32] I. Y. Chun, H. Lim, Z. Huang, and J. A. Fessler, “Fast and convergent iterative image recovery using trained convolutional neural networks,” in _Allerton Conf. on Comm., Control, and Computing_ , 2018, pp. 155–9, invited.
  * [33] I. Y. Chun, Z. Huang, H. Lim, and J. A. Fessler, “Momentum-Net: Fast and convergent iterative neural network for inverse problems,” 2019. [Online]. Available: http://arxiv.org/abs/1907.11818
  * [34] S. Ye, Z. Li, M. T. McCann, Y. Long, and S. Ravishankar, “Unified Supervised-Unsupervised (SUPER) Learning for X-Ray CT Image Reconstruction,” _IEEE Transactions on Medical Imaging_ , vol. 40, no. 11, pp. 2986–3001, 2021.
  * [35] A. Lahiri, G. Wang, S. Ravishankar, and J. A. Fessler, “Blind primed supervised (BLIPS) learning for MR image reconstruction,” _IEEE Trans. Med. Imag._ , vol. 40, no. 11, pp. 3113–24, Nov. 2021.
  * [36] G. Ongie, A. Jalal, C. A. Metzler, R. G. Baraniuk, A. G. Dimakis, and R. Willett, “Deep learning techniques for inverse problems in imaging,” _IEEE Journal on Selected Areas in Information Theory_ , vol. 1, no. 1, pp. 39–56, may 2020.
  * [37] M. Hossain, B. T. Nadiga, O. Korobkin, M. L. Klasky, J. L. Schei, J. W. Burby, M. T. McCann, T. Wilcox, S. De, and C. A. Bouman, “High-precision inversion of dynamic radiography using hydrodynamic features,” 2021. [Online]. Available: http://arxiv.org/abs/2112.01627
  * [38] E. T. Quinto, “Exterior and limited-angle tomography in non-destructive evaluation,” _Inverse Problems_ , vol. 14, no. 2, p. 339, 1998.
  * [39] S. Guan, A. A. Khan, S. Sikdar, and P. V. Chitnis, “Limited-view and sparse photoacoustic tomography for neuroimaging with deep learning,” _Scientific reports_ , vol. 10, no. 1, pp. 1–12, 2020.
  * [40] M. Xiang, “Deep learning-based reconstruction of volumetric CT images of vertebrae from a single view X-ray image,” 2020, theses and Dissertations. 2629\. University of Wisconsin Milwaukee. [Online]. Available: https://dc.uwm.edu/etd/2629
  * [41] A. Ziabari, D. H. Ye, S. Srivastava, K. D. Sauer, J. Thibault, and C. A. Bouman, “2.5D deep learning for CT image reconstruction using A multi-GPU implementation,” in _asccs_ , 2018, pp. 2044–9.
  * [42] H. Lim, I. Y. Chun, Y. K. Dewaraja, and J. A. Fessler, “Improved low-count quantitative PET reconstruction with a variational neural network,” _IEEE Trans. Med. Imag._ , vol. 39, no. 11, pp. 3512–22, Nov. 2020.
  * [43] G. Corda-D’Incan, J. A. Schnabel, and A. J. Reader, “Memory-efficient training for fully unrolled deep learned PET image reconstruction with iteration-dependent targets,” _IEEE Trans. Radiation and Plasma Med. Sci._ , 2021.
  * [44] M. Magnusson, P.-E. Danielsson, and J. Sunnegardh, “Handling of long objects in iterative improvement of nonexact reconstruction in helical cone-beam CT,” _IEEE Trans. Med. Imag._ , vol. 25, no. 7, pp. 935–40, Jul. 2006\.
  * [45] V. Monga, Y. Li, and Y. C. Eldar, “Algorithm unrolling: interpretable, efficient deep learning for signal and image processing,” _IEEE Sig. Proc. Mag._ , vol. 38, no. 2, pp. 18–44, Mar. 2021.
  * [46] H. Nien and J. A. Fessler, “Relaxed linearized algorithm for faster X-ray CT image reconstruction,” in _Proc. Intl. Mtg. on Fully 3D Image Recon. in Rad. and Nuc. Med_ , 2015, pp. 260–3. [Online]. Available: http://web.eecs.umich.edu/~fessler/papers/files/proc/15/web/nien-15-rla.pdf
  * [47] I. J. Goodfellow, J. Pouget-Abadie, M. Mirza, B. Xu, D. Warde-Farley, S. Ozair, A. Courville, and Y. Bengio, “Generative adversarial networks,” 2014. [Online]. Available: http://arxiv.org/abs/1406.2661
  * [48] Y. Long, J. A. Fessler, and J. M. Balter, “3D forward and back-projection for X-ray CT using separable footprints,” _IEEE Trans. Med. Imag._ , vol. 29, no. 11, pp. 1839–50, Nov. 2010.
  * [49] H. D. Sarkissian, F. Lucka, M. . Eijnatten, G. Colacicco, S. B. Coban, and K. J. Batenburg, “A cone-beam X-ray computed tomography data collection designed for machine learning,” _Sci. Data_ , vol. 6, p. 215, 2019.
  * [50] J. A. Fessler, “Michigan image reconstruction toolbox (MIRT) for Matlab,” 2016, available from http://web.eecs.umich.edu/ fessler/irt/irt.
  * [51] S. Ravishankar and Y. Bresler, “MR image reconstruction from highly undersampled k-space data by dictionary learning,” _IEEE Trans. Med. Imag._ , vol. 30, no. 5, pp. 1028–41, May 2011.


11institutetext: Utrecht University, the Netherlands, 11email:
<EMAIL_ADDRESS>
# Test Model Coverage Analysis under Uncertainty

I. S. W. B. Prasetya 11 0000−0002−3421−4635 Rick Klomp 11

###### Abstract

In model-based testing (MBT) we may have to deal with a non-deterministic
model, e.g. because abstraction was applied, or because the software under
test itself is non-deterministic. The same test case may then trigger multiple
possible execution paths, depending on some internal decisions made by the
software. Consequently, performing precise test analyses, e.g. to calculate
the test coverage, are not possible. This can be mitigated if developers can
annotate the model with estimated probabilities for taking each transition. A
probabilistic model checking algorithm can subsequently be used to do simple
probabilistic coverage analysis. However, in practice developers often want to
know what the achieved aggregate coverage, which unfortunately cannot be re-
expressed as a standard model checking problem. This paper presents an
extension to allow efficient calculation of probabilistic aggregate coverage,
and moreover also in combination with $k$-wise coverage.

###### Keywords:

probabilistic model based testing, probabilistic test coverage, testing non-
deterministic systems

This is a preprint. The final version of this paper is published in the
Proceedings of the International Conference on Software Engineering and Formal
Methods SEFM 2019, https://doi.org/10.1007/978-3-030-30446-1_12.

##  1 Introduction

Model based testing (MBT) is considered as one of the leading technologies for
systematic testing of software [17, 5, 6]. It has been used to test different
kinds of software, e.g. communication protocols, web applications, and
automotive control systems. In this approach, a model describing the intended
behavior of the system under test (SUT) is first constructed [27], and then
used to guide the tester, or a testing algorithm, to systematically explore
and test the SUT’s states. Various automated MBT tools are available, e.g.
JTorX [4, 26], Phact [11], OSMO [14], APSL [24], and RT-Tester [17].

There are situations where we end up with a non-deterministic model [17, 23,
13], for example when the non-determinism within the system under test, e.g.
due to internal concurrency, interactions with an uncontrollable environment
(e.g. as in cyber physical systems), or use of AI, leads to observable effects
at the model level. Non-determinism can also be introduced as byproduct when
we apply abstraction on an otherwise too large model [20]. Models mined from
executions logs [21, 7, 28] can also be non-deterministic, because log files
only provide very limited information about a system’s states.

MBT with a non-deterministic model is more challenging. The tester cannot
fully control how the SUT would traverse the model, and cannot thus precisely
determine the current state of the SUT. Obviously, this makes the task
deciding which trigger to send next to the SUT harder. Additionally, coverage,
e.g. in terms of which states in the model have been visited by a series of
tests, cannot be determined with 100% certainty either. This paper will focus
on addressing the latter problem —readers interested in test cases generation
from non-deterministic models are referred to e.g. [25, 16, 13]. Rather than
just saying that a test sequence may cover some given state, we propose to
calculate the probability of covering a given coverage goal, given modelers’
estimation on the local probability of each non-deterministic choice in a
model.

Given a probabilistic model of the SUT, e.g. in the form of a Markov Decision
Process (MDP) [3, 22], and a test $\sigma$ in the form of a sequence of
interactions on the SUT, the most elementary type of coverage goal in MBT is
for $\sigma$ is to cover some given state $s$ of interest in the model.
Calculating the probability that this actually happens is an instance of the
probabilistic reachability problem which can be answered using e.g. a
probabilistic model checker [10, 3, 15]. However, in practice coverage goals
are typically formulated in an ’aggregate’ form, e.g. to cover at least 80% of
the states, without being selective on which states to include. Additionally,
we may want to know the aggregate coverage over pairs of states (the
transitions in the LTS), or vectors of states, as in $k$-wise coverage [1], as
different research showed that $k$-wise greatly increases the fault finding
potential of a test suite [18, 9]. Aggregate goals cannot be expressed in LTL
or CTL, which are the typical formalisms in model checking. Furthermore, both
types of goals (aggregate and $k$-wise) may lead to combinatorial explosion.

This paper contributes: (1) a concept and definition of probabilistic test
coverage; as far as we know this has not been covered in the literature
before, and (2) an algorithm to calculate probabilistic coverage, in
particular of aggregate $k$-wise coverage goals.

Paper structure. Section 2 introduces relevant basic concepts. Section 3
introduces the kind of coverage goals we want to be able to express and how
their probabilistic coverage can be calculated. Section 4 presents our
algorithm for efficient coverage calculation. Section 5 shows the results of
our benchmarking. Related work is discussed in Section 6. Section 7 concludes.

##  2 Preliminary: Probabilistic Models and Simple Coverage

As a running example, consider the labelled transition system (LTS) [2] in
Fig. 1 as a model of some SUT. The transitions are labelled with actions, e.g.
$a$ and $b$. A non-$\tau$ action represents an interaction between the SUT and
its environment. In our set up such an action is assumed to occur
synchronously a la CSP [12] (for an action $a$ to take place, both the SUT and
the environment first need to agree on doing $a$; then they will do $a$
together). The action $\tau$ represents an internal action by the SUT, that is
not visible to the environment.

0$\sf
EX_{1}$:123456$(0.5)a$$(0.5)a$$(0.1)b$$(0.9)b$$b$$(0.9)\tau$$(0.1)\tau$$a$$c$
Figure 1:  An example of a probabilistic model of some SUT called $\sf
EX_{1}$.

To test the SUT, the tester controls the SUT by insisting on which action it
wants to synchronize; e.g. if on the state $t$ the SUT is supposed to be able
to either do $a$ or $b$, the tester can insist on doing $a$. If the SUT fails
to go along with this, it is an error. The tester can also test if in this
state the SUT can be coerced to do an action that it is not supposed to
synchronize; if so, the SUT is incorrect. We will assume a black box setup.
That is, the tester cannot actually see the SUT’s state, though he/she can try
to infer this based on information visible to him/her, e.g. the trace of the
external actions done so far. For example after doing $a$ on the SUT $\sf
EX_{1}$ above, the tester cannot tell whether it then goes to the state 1 or
2. However, if the tester manages to do $abc$ he/she would on the hind sight
know that the state after $a$ must have been 1.

When a state $s$ has multiple outgoing transitions with the same label, e.g.
$a$, this implies non-determinism, since the environment cannot control which
$a$ the SUT will take (the environment can only control whether or not it
wants to do $a$). We assume the modeler is able estimate the probability of
taking each of these $a$-transition and annotate this on the transition. E.g.
in Fig. 1 we see that in state 1, two $a$-transitions are possible, leading to
different states, each with the probability of 0.5. Similarly, in state 3
there are two $\tau$-transitions leading to states 4 and 5, with the
probability of 0.9 and 0.1 respectively. A probabilistic model such as in
Figure 1 is also called a Markov Decision Process (MDP) [3].

Let in the sequel $M$ be an MDP model, with finite number of transitions, and
a single initial state. Let $s,t$ be states, and $a$ an action. We write
$s{\in}M$ to mean that $s$ is a state in $M$. The notation
$s{\stackrel{{\scriptstyle a}}{{\longrightarrow}}}t$ denotes a transition that
goes from the state $s$ to $t$ and is labelled with $a$. We write
$s{\stackrel{{\scriptstyle a}}{{\longrightarrow}}}t\in M$ to mean that
$s{\stackrel{{\scriptstyle a}}{{\longrightarrow}}}t$ is a transition in $M$.
$P_{M}(s{\stackrel{{\scriptstyle a}}{{\longrightarrow}}}t)$ denotes the
probability that $M$ will take this particular transition when it synchronizes
over $a$ on the state $s$.

To simplify calculation over non-deterministic actions, we will assume that
$M$ is $\tau$-normalized in the following sense. First, a state cannot have a
mix of $\tau$ and non-$\tau$ outgoing transitions. E.g. a state $s$ with two
transitions
$\\{s{\stackrel{{\scriptstyle\tau}}{{\longrightarrow}}}t,s{\stackrel{{\scriptstyle
a}}{{\longrightarrow}}}u\\}$ should first be re-modelled as
$\\{s{\stackrel{{\scriptstyle\tau}}{{\longrightarrow}}}t,s{\stackrel{{\scriptstyle\tau}}{{\longrightarrow}}}s^{\prime},s^{\prime}{\stackrel{{\scriptstyle
a}}{{\longrightarrow}}}u\\}$ by introducing an intermediate state
$s^{\prime}$, and the modeler should provide estimation on the probability of
taking each of the two $\tau$ transitions. Second, $M$ should have no state
whose all incoming and outgoing transitions are $\tau$ transitions. Such a
state is considered not interesting for our analyses. Third, $M$ should not
contain a cycle that consists of only $\tau$ transitions. In a
$\tau$-normalized model, non-determinism can only be introduced if there is a
state $s$ with multiple outgoing transitions labelled by the same action
(which can be $\tau$).

We define an execution of the SUT as a finite path $\rho$ through the model
starting from its initial state. A trace is a finite sequence of external
actions. The trace of $\rho$, ${\sf tr}(\rho)$, is the sequence external
actions induced by $\rho$. A legal trace is a trace that can be produced by
some execution of the SUT. A test-case is abstractly modeled by a trace. We
will restrict to test-cases that form legal traces, e.g. $ab$, $aba$, and
$ababc$ are possible test cases for $\sf Ex1$ in Fig. 1. A set of test cases
is also called a test suite. Since the model can be non-deterministic, the
same test case may trigger multiple possible executions which are
indistinguishable from their trace. If $\sigma$ is a trace, ${\sf
exec}(\sigma)$ denotes the set of all executions $\rho$ such that ${\sf
tr}(\rho){=}\sigma$, and moreover each such $\rho$ is $\tau$-maximal: it
cannot be extended without breaking the property ${\sf tr}(\rho){=}\sigma$.

###  2.1 Representing a test case: execution model

The probability that a test case $\sigma$ covers some goal $\phi$ (e.g. a
particular state $s$) can in principle be calculated by quantifying over ${\sf
exec}(\sigma)$. However, if $M$ is highly non-deterministic, the size of ${\sf
exec}(\sigma)$ can be exponential with respect to the length of $\sigma$. To
facilitate more efficient coverage calculation we will represent $\sigma$ with
the subgraph of $M$ that $\sigma$ induces, called the execution model of
$\sigma$, denoted by ${\sf E}(\sigma)$. ${\sf E}(\sigma)$ forms a Markov
chain; each branch in ${\sf E}(\sigma)$ is annotated with the probability of
taking the branch, under the premise that $\sigma$ has been observed. Since a
test case is always of finite length and $M$ is assumed to have no
$\tau$-cycle, ${\sf E}(\sigma)$ is always acyclic. Typically the size or ${\sf
E}(\sigma)$ (its number of nodes) is much less than the size of ${\sf
exec}(\sigma)$. For example, the execution model of the test case $aba$ on
$\sf EX_{1}$ is shown in Fig. 2. An artificial state denoted with $\sharp$ is
added so that ${\sf E}(\sigma)$ has a single exit node, which is convenient
for later.

$u_{0}(0)$$u_{1}(1)$$u_{2}(2)$$u_{4}(0)$$u_{3}(3)$$u_{7}(2)$$u_{5}(4)$$u_{6}(1)$$u_{8}(\sharp)$$(0.5)a$$(0.5)a$$b$$(0.1)b$$(0.9)b$$(0.5)a$$(0.5)a$$\tau$$a$
Figure 2: The execution model of the test case $aba$ on $\sf EX_{1}$.

To identify the states in ${\sf E}(\sigma)$ we assign IDs to them
($u_{0}...u_{8}$ in Fig. 2). We write $u.{\sf st}$ to denote $u$’s state
label, which is the ID of a state in $M$ that $u$ represents (so, $u.{\sf
st}\in M$); in Fig. 2 this is denoted by the number between brackets in every
node.

Importantly, notice that the probability of the transitions in ${\sf
E}(\sigma)$ may be different than the original probability in $M$. For
example, the transition
$u_{3}{\stackrel{{\scriptstyle\tau}}{{\longrightarrow}}}u_{5}$ in the above
execution model has probability 1.0, whereas in the original model $\sf
EX_{1}$ this corresponds to the transition
$3{\stackrel{{\scriptstyle\tau}}{{\longrightarrow}}}4$ whose probability is
0.9. This is because the alternative
$3{\stackrel{{\scriptstyle\tau}}{{\longrightarrow}}}5$ could not have taken
place, as it leads to an execution whose trace does not correspond to the test
case $aba$ (which is assumed to have happened).

More precisely, when an execution in the model ${\sf E}(\sigma)$ reaches a
node $u$, the probability of extending this execution with the transition
$u{\stackrel{{\scriptstyle\alpha}}{{\longrightarrow}}}v$ can be calculated by
taking the conditional probability of the corresponding transition in the
model $M$, given that only the outgoing transitions specified by ${\sf
E}(\sigma)$ could happen. So, $P_{{\sf
E}(\sigma)}(u{\stackrel{{\scriptstyle\alpha}}{{\longrightarrow}}}v)$ is
$P_{M}(u.{\sf st}{\stackrel{{\scriptstyle\alpha}}{{\longrightarrow}}}v.{\sf
st})$ divided by the the sum of $P_{M}(u.{\sf
st}{\stackrel{{\scriptstyle\alpha}}{{\longrightarrow}}}w.{\sf st})$ of all $w$
such that $u{\stackrel{{\scriptstyle\alpha}}{{\longrightarrow}}}w\in{\sf
E}(\sigma)$.

Let $E={\sf E}(\sigma)$. Since $E$ is thus acyclic, the probability that SUT
traverses a path/ execution $\rho$ in ${\sf E}(\sigma)$ when it is given
$\sigma$ can be obtained by multiplying the probability of all the transitions
in the path:

$P_{E}(\rho)\ =\
\prod_{s{\stackrel{{\scriptstyle\alpha}}{{\longrightarrow}}}t\in\rho}P_{E}(s{\stackrel{{\scriptstyle\alpha}}{{\longrightarrow}}}t)$
(1)

#### Simple coverage analyses.

As an example of a simple analysis, let’s calculate the probability that a
test case $\sigma$ produces an execution that passes through a given state
$s$, denoted by $P(\langle s\rangle\;|\;\sigma)$. This would then just be the
sum of the probability of all full executions in ${\sf E}(\sigma)$ that
contain $s$. So:

$P(\langle s\rangle\;|\;\sigma)\ =\ \sum_{\rho\;\mbox{s.t.}\;\rho\in{\sf
E}(\sigma)\wedge s\in\rho}P_{{\sf E}(\sigma)}(\rho)$ (2)

For example, on the execution model $\sf EX_{1}$, $P(\langle
1\rangle\;|\;aba)=0.75$, $P(\langle 2\rangle\;|\;aba)=0.725$, $P(\langle
4\rangle\;|\;aba)=0.05$, whereas $P(\langle 5\rangle\;|\;aba)=0$.

##  3 Coverage under Uncertainty

Coverage goals posed in practice are however more complex than goals
exemplified above. Let us first introduce a language for expressing ’goals’;
we will keep it simple, but expressive enough to express what is latter called
’aggregate $k$-wise’ goals. A goal of the form $\langle 0,2,0\rangle$ is
called a word, expressing an intent to cover the subpath $\langle
0,2,0\rangle$ in the MDP model. We will also allow disjunctions of words and
sequences of words (called sentences) to appear as goals. For example:
$(\langle 0,2\rangle\vee\langle 1,0\rangle)\;;\;\langle 1\rangle$ formulates a
goal to first cover the edge $0{\rightarrow}2$ or $1{\rightarrow}0$, and then
(not necessarily immediately) the node 1.

The typical goal people have in practice is to cover at least $p\%$ of the
states. This is called an aggregate goal. We write this a bit differently: a
goal of the form ${}^{1}{\geq}N$ express an intent to cover at least $N$
different states. Covering at least $p\%$ can be expressed as
${}^{1}{\geq}\lfloor p*K/100\rfloor$ where $K$ is the number of states in the
model. To calculate probabilistic coverage in $k$-wise testing [1], the goal
${}^{k}{\geq}N$ expresses an intent to cover at least $N$ different words of
length $k$. Now, more formally:

###### Definition 1

A coverage goal is a formula $\phi$ with this syntax:

$\begin{array}[]{lclr}\phi&::=&S\ |\ A&\mbox{(goal)}\\\ S&::=&C\ |\
C;S&\mbox{(sentence)}\\\ A&::=&{}^{k}{\geq}N&\mbox{(aggregate goal), with
$k{\geq}1$}\\\ C&::=&W\ |\ W{\vee}C&\mbox{(clause)}\\\ W&::=&\langle
s_{0},...,s_{k{-}1}\rangle&\mbox{(word), with $k{\geq}1$}\end{array}$

A sentence is a sequence $C_{0};C_{1};...$. Each $C_{i}$ is called a clause,
which in turns consists of one or more words. A word is denoted by $\langle
s_{0},s_{1},...\rangle$ and specifies one or more connected states in an MDP.

Let $\rho$ be an execution. If $\phi$ is a goal, we write $\rho\vdash\phi$ to
mean that $\rho$ covers $\phi$. Checking this is decidable. For a word $W$,
$\rho\vdash W$ if $W$ is a segment of $\rho$. For a clause $C=W_{0}\vee
W_{2}\vee...$, $\rho\vdash C$ if $\rho\vdash W_{k}$ for some $k$. Roughly, a
sentence $C_{0};C_{1};...$ is covered by $\rho$ if all clauses $C_{i}$ are
covered by $\rho$, and furthermore they are covered in the order as specified
by the sentence. We will however define it more loosely to allow consecutive
clauses to overlap, as follows:

######  Definition 2 (Sentence Coverage)

Let $S$ be a sentence. (1) An empty $\rho$ does not cover $S$. (2) If $S$ is a
just a single clause $C$, then $\rho\vdash S$ iff $\rho\vdash C$. (3) If
$S=C;S^{\prime}$ and a prefix of $\rho$ matches one of the words in $C$, then
$\rho\vdash S$ iff $\rho\vdash S^{\prime}$. If $\rho$ has no such prefix, then
$\rho\vdash S$ iff ${\sf tail}(\rho)\vdash S$.

An aggregate goal of the form ${}^{k}{\geq}N$ is covered by $\rho$ if $\rho$
covers at least $N$ different words of size $k$. While sentences are
expressible in temporal logic, aggregate goals are not. This has an important
consequence discussed later.

Let $\phi$ be a coverage goal and $\sigma$ a test case. Let’s write
$P(\phi\;|\;\sigma)$ to denote the probability that $\phi$ is covered by
$\sigma$, which can be calculated analogous to (2) as follows:

###### Definition 3

$P(\phi\;|\;\sigma)$ is equal to $P(\phi\;|\;E)$ where $E={\sf E}(\sigma)$,
$P(\phi\;|\;E)\ =\ \sum_{\rho\;\mbox{s.t.}\;\rho\in{\sf
exec}(E)\;\wedge\;\rho\vdash\phi}P_{E}(\rho)$, and where $P_{E}(\rho)$ is
calculated as in (1).

For example, consider the test case $aba$ on the SUT ${\sf EX}_{1}$. Fig. 2
shows the execution model of $aba$. $P(\langle 2,0\rangle\;|\;aba)$ is the
probability that $aba$’s execution passes through the transition
$2{\rightarrow}0$; this probability is 0.5. $P((\langle 2\rangle\vee\langle
3\rangle);\langle 1\rangle\;|\;aba)$ is the probability that $aba$ first
visits the state 2 or 3, and sometime later 1; this probability is 0.75.
$P(^{1}{\geq}4\;|\;aba)$ is the probability that the execution of $aba$ visits
at least four different states; this is unfortunately only 0.05.

Due to non-determinism, the size of ${\sf exec}(\sigma)$ could be exponential
with respect to the length of $\sigma$. Simply using using the formula in Def.
3 would then be expensive. Below we present a much better algorithm to do the
calculation.

##  4 Efficient coverage calculation

Coverage goals in the form of sentences are actually expressible in
Computation Tree Logic (CTL) [3]. E.g. $\langle s,t\rangle;\langle u\rangle$
corresponds to ${\sf EF}(s\wedge t\wedge{\sf EF}u)$. It follows, that the
probability of covering a sentence can be calculated through probabilistic CTL
model checking [10, 3]. Unfortunately, aggregate goals are not expressible
CTL. Latter we will discuss a modification of probabilistic model checking to
allow the calculation of aggregate goals. We first start with the calculation
of simple sentences whose words are all of the length one.

Let $S$ be a simple sentence, $\sigma$ a test case, and $E={\sf E}(\sigma)$.
In standard probabilistic model checking, $P(S|\sigma)$ would be calculated
through a series multiplications over a probability matrix [3]. We will
instead do it by performing labelling on the nodes of $E$, resembling more to
non-probabilistic CTL model checking. This approach is more generalizable to
later handle aggregate goals.

Notice that any node $u$ in $E$ induces a unique subgraph, denoted by $E@u$,
rooted in $u$. It represents the remaining execution of $\sigma$, starting at
$u$. When we label $E$ with some coverage goal $\psi$, the labelling will
proceed in such a way that when it terminates every node $u$ in $E$ is
extended with labels of the form $u.{\sf lab}(\psi)$ containing the value of
$P(\psi\;|\;E@u)$. The labelling algorithm is shown in Fig. 3, namely the
procedure ${\sf label}(..)$ —we will explain it below. In any case, after
calling ${\sf label}(E,S)$, the value of $P(S\;|\;\sigma)$ can thus be
obtained simply by inspecting the ${\sf lab}(S)$ of $E$’s root node. This is
done by the procedure $\sf calcSimple$.

1:procedure calcSimple($E,S$)

2: ${\sf label}(E,S)$

3: return ${\sf root}(E).{\sf lab}(S)$

4:end procedure

5:

6:procedure label($E,S$)

7: $u_{0}\leftarrow{\sf root}(E)$

8: ${\bf case}\;S\;{\bf of}$

9: $C\hskip 17.07164pt\rightarrow{\sf label1}(u_{0},C)$

10: $C;S^{\prime}\ \rightarrow{\sf label}(E,S^{\prime})\;;\;{\sf
label1}(u_{0},S)$

11:end procedure

12:

13:procedure checkClause($u,C$)

14: $\rhd$ the clause $C$ is assumed to be of this form, with $k\geq 1$ :

15: let $\langle s_{0}\rangle\vee...\vee\langle s_{k{-}1}\rangle=C$

16: $isCovered\leftarrow u.{\sf st}\in\\{s_{0},...,s_{k{-}1}\\}$

17: return $isCovered$

18:end procedure

19:

20:procedure label1($u,S$)

21: $\rhd$ recurse to $u$’s successors :

22: ${\bf forall}\;v\in u.{\sf next}\ \rightarrow\ {\sf label1}(v,S)$

23: $\rhd$ pre-calculate $u$’s successors’ total probability to cover $S$ :

24: $q^{\prime}\leftarrow\sum_{v{\in}u.{\sf next}}u.{\sf pr}(v)*v.{\sf
lab}(S)$

25: $\rhd$ calc. $u$’s probability to cover $S$ :

26: ${\bf case}\;S\;{\bf of}$

27: $C\hskip 17.07164pt\rightarrow\begin{array}[t]{l}{\bf if}\;{\sf
checkClause}(u,C)\\\ {\bf then}\;q\leftarrow 1\\\ {\bf else}\;q\leftarrow
q^{\prime}\end{array}$

28: $C;S^{\prime}\ \rightarrow\begin{array}[t]{l}{\bf if}\;{\sf
checkClause}(u,C)\\\ {\bf then}\;q\leftarrow u.{\sf lab}(S^{\prime})\\\ {\bf
else}\;q\leftarrow q^{\prime}\end{array}$

29: ${\bf end\;case}$

30: $\rhd$ add the calculated probability as a new label to $u$ :

31: $u.{\sf lab}(S)\leftarrow q$

32:end procedure

Figure 3: The labeling algorithm of to calculate the probability of simple
sentences.

Since $S$ is a sentence, it is a sequence of clauses. The procedure ${\sf
label}(E,S)$ first recursively labels $E$ with with the tail $S^{\prime}$ of
$S$ (line 9), then we proceed with the labelling with $S$ itself, which is
done by the procedure $\sf label1$. In $\sf label1$, the following notations
are used. Let $u$ be a node in $E$. Recall that $u.{\sf st}$ denotes the ID of
the state in $M$ that $u$ represents. We write $u.{\sf next}$ to denote the
set of $u$’s successors in $E$ (and not in $M$!). For such successor $v$,
$u.{\sf pr}(v)$ denotes the probability annotation that $E$ puts on the arrow
$u{\rightarrow}v$. A label is a pair $(\psi,p)$ where $\psi$ is a coverage
goal and $p$ is probability in $[0..1]$. The notation $u.{\sf lab}$ denotes
the labels put so far to the node $u$. The assignment $u.{\sf
lab}(\psi)\leftarrow p$ adds the label $(\psi,p)$ to $u$, and the expression
$u.{\sf lab}(\psi)$ returns now the value of $p$.

The procedure ${\sf label1}(\psi)$ will perform the labelling node by node
recursively in the bottom-up direction over the structure of $E$ (line 19).
Since $E$ is acyclic, only a single pass of this recursion is needed. For
every node $u\in E$, ${\sf label1}(u,S)$ has to add a new label $(S,q)$ to the
node $u$ where $q$ is the probability that the goal $S$ is covered by the part
of executions of $\sigma$ that starts in $u$ (in other words, the value of
$P(S\;|\;E@u)$). The goal $S$ will be in one of these two forms:

  1. 1.

$S$ is just a single clause $C$ (line 24). Because $S$ is a simple sentence,
$C$ is a disjunction of singleton words $\langle
s_{0}\rangle\vee...\vee\langle s_{k{-}1}\rangle$, where each $s_{i}$ is an ID
of a state in $M$. If $u$ represents one of these states, the probability that
$E@u$ covers $C$ would be 1. Else, it is the sum of the probability to cover
$C$ through $u$’s successors (line 20). As an example, Fig. 4 (left) shows how
the labeling of a simple sentence $\langle 1\rangle$ on the execution model in
Fig. 2 proceeds.

  2. 2.

$S$ is a sentence with more than one clauses; so it is of the form
$C;S^{\prime}$ (line 25) where $C$ is a clause and $S^{\prime}$ is the rest of
the sentence, we calculate the coverage probability by $E@u$ by basically
following the third case in Def. 2.

As an example, Fig. 4 (right) shows how the labeling of $S=\langle
0\rangle{;}\langle 1\rangle$ proceeds. At every node $u$ we first check if $u$
covers the first word, namely $\langle 0\rangle$. If this is the case, the
probability that $E@u$ covers $S$ would be the same as the probability that it
covers the rest of $S$, namely $\langle 1\rangle$. The probability of the
later is by now known, calculated by $\sf label$ in its previous recursive
call. The result can be inspected in $u.{\sf lab}(\langle 1\rangle)$.

If $u$ does not cover $S$, the probability that ${\sf E}(u)$ covers $S$ would
be the sum of the probability to cover $S$ through $u$’s successors
(calculated in line 21).

$\begin{array}[]{c}u_{0}(0)\\\ \langle 1\rangle:0.75\end{array}$$\begin{array}[]{c}u_{1}(1)\\\ \langle 1\rangle:1\end{array}$$\begin{array}[]{c}u_{2}(2)\\\ \langle 1\rangle:0.5\end{array}$$\begin{array}[]{c}u_{3}(3)\\\ \langle 1\rangle:1\end{array}$$\begin{array}[]{c}u_{4}(0)\\\ \langle 1\rangle:0.5\end{array}$$\begin{array}[]{c}u_{5}(4)\\\ \langle 1\rangle:1\end{array}$$\begin{array}[]{c}u_{6}(1)\\\ \langle 1\rangle:1\end{array}$$\begin{array}[]{c}u_{7}(2)\\\ \langle 1\rangle:0\end{array}$$\sharp$$0.5$$0.5$$0.1$$0.9$$0.5$$0.5$ | $\begin{array}[]{c}u_{0}(0)\\\ \langle 0\rangle{;}\langle 1\rangle{:}0.75\end{array}$$\begin{array}[]{c}u_{1}(1)\\\ \langle 0\rangle{;}\langle 1\rangle{:}0.45\end{array}$$\begin{array}[]{c}u_{2}(2)\\\ \langle 0\rangle{;}\langle 1\rangle{:}0.5\end{array}$$\begin{array}[]{c}u_{3}(3)\\\ \langle 0\rangle{;}\langle 1\rangle{:}0\end{array}$$\begin{array}[]{c}u_{4}(0)\\\ \langle 0\rangle{;}\langle 1\rangle{:}0.5\end{array}$$\begin{array}[]{c}u_{5}(4)\\\ \langle 0\rangle{;}\langle 1\rangle{:}0\end{array}$$\begin{array}[]{c}u_{6}(1)\\\ \langle 0\rangle{;}\langle 1\rangle{:}0\end{array}$$\begin{array}[]{c}u_{7}(2)\\\ \langle 0\rangle{;}\langle 1\rangle{:}0\end{array}$$\sharp$$0.5$$0.5$$0.1$$0.9$$0.5$$0.5$
---|---
Figure 4: The left graph shows the result of ${\sf label}(\langle 1\rangle)$
on the execution model of $aba$ in Fig.2. For simplicity, the action labels on
the arrows are removed. The probability annotation is kept. In turn, ${\sf
label}()$ calls $\sf label1$, which then performs the labelling recursively
from right to left. The nodes $u_{6}$ and $u_{7}$ (yellow) are base cases. The
probabilities of $\langle 1\rangle$ on them are respectively 1 and 0. This
information is then added as the labels of these nodes. Next, $\sf label1$
proceeds with the labelling of $u_{4}$ and $u_{5}$. E.g. on $u_{4}$ (orange),
because $u_{4}.{\sf st}$ is not 1, for $u_{4}$ to cover $\langle 1\rangle$ we
need an execution that goes through $u_{6}$, with the probability of 0.5. So
the probability of $\langle 1\rangle$ on $u_{4}$ is 0.5. The right graph shows
the result of ${\sf label}(\langle 0\rangle{;}\langle 1\rangle)$ on the same
execution model. This will first call ${\sf label}(\langle 1\rangle)$, thus
producing the labels as shown in the left graph, then proceeds with ${\sf
label1}(\langle 0\rangle{;}\langle 1\rangle)$. Again, $\sf label1$ performs
the labelling recursively from right to left. The base cases $u_{6}$ and
$u_{7}$ do not cover $\langle 0\rangle{;}\langle 1\rangle$, so the
corresponding probability there is 0. Again, this information is added as
labels of the corresponding nodes. Node $u_{4}$ (orange) has $u_{4}.{\sf
st}=0$. So, any execution that starts from there and covers $\langle 1\rangle$
would also cover $\langle 0\rangle{;}\langle 1\rangle$. The probability that
$u_{4}$ covers $\langle 1\rangle$ is already calculated in the left graph,
namely 0.5. So this is also the probability that it covers $\langle
0\rangle{;}\langle 1\rangle$.

Assuming that checking if a node locally covers a clause (the procedure $\sf
checkClause$ in Fig. 3) takes a unit time, the time complexity of $\sf label1$
is ${\mathcal{O}}(|E|)$, where $|E|$ is the size of $E$ in terms of its the
number of edges. The complexity of $\sf label$ is thus
${\mathcal{O}}(|E|*|S|)$, where $|S|$ is the size of the goal $S$ in terms of
the number of clauses it has. The size of $E$ is typically just linear to the
length of the test case: ${\mathcal{O}}(N_{sucs}*|\sigma|)$, where $N_{sucs}$
is the average number of successors that each state in $M$ has. This is
significant improvement compared to the exponential run time that we would get
if we simply use Def. 3.

###  4.1 Non-simple sentences

Coverage goals in $k$-wise testing would require sentences with words of
length $k{>}1$ to express. These are thus non-simple sentences. We will show
that the algorithm in Fig. 3 can be used to handle these sentences as well.

Consider as an example the sentence $\langle 0,2,0\rangle;\langle
4,1,\sharp\rangle$. The words are of length three, so the sentence is non-
simple. Suppose we can treat these words as if they are singletons. E.g. in
$\langle 0,2,0\rangle$ the sequence $0,2,0$ is treated as a single symbol, and
hence the word is a singleton. From this perspective, any non-aggregate goal
is thus a simple sentence, and therefore the algorithm in Fig. 3 can be used
to calculate its coverage probability. We do however need to pre-process the
execution model to align it with this idea.

The only part of the algorithm in Fig. 3 where the size of the words matters
is in the procedure ${\sf checkClause}$. Given a node $u$ in the given
execution model $E$ and a clause $C$, ${\sf checkClause}(u,C)$ checks if the
clause $C$ is covered by $E$’s executions that start at $u$. If the words in
$C$ are all of length one, $C$ can be immediately checked by knowing which
state in $M$ $u$ represents. This information is available in the attribute
$u.{\sf st}$. Clauses with longer words can be checked in a similar way. For
simplicity, assume that the words are all of length $k$. We first restructure
$E$ such that the ${\sf st}$ attribute of every node $u$ in the new $E$
contains a word of length $k$ that would be covered if the execution of $E$
arrives at $u$. We call this restructuring step $k$-word expansion. Given a
base execution model $E$, the produced new execution model will be denoted by
$E^{k}$. As an example, the figure below shows the word expansion with $k{=}3$
of the execution model in Fig. 2 (for every node $v$ we only show its $v.{\sf
st}$ label, which is an execution segment of length 3). Artificial initial and
terminal states are added to the new execution model, labelled with $\sharp$.
When a word of length $k$ cannot be formed, because the corresponding segment
has reached the terminal state $\sharp$ in $E$, we pad the word with
$\sharp$’s on its the end until its length is $k$.

$\sharp$[0,1,3][1,3,4][3,4,1][0,1,0][1,0,1]$[0,1,\sharp]$[1,0,2]$[0,2,\sharp]$[0,2,0][2,0,1]$[0,1,\sharp]$[2,0,2]$[0,2,\sharp]$$[4,1,\sharp]$$[1,\sharp,\sharp]$$[2,\sharp,\sharp]$$\sharp$0.050.450.50.50.50.50.5

###  4.2 Coverage of aggregate goals

We will only discuss the calculation of aggregate goals of the form
${}^{k}{\geq}N$ where $k{=}1$. If $k{>}1$ we can first apply a $k$-word
expansion (Section 4.1) on the given execution model $E$, then we calculate
${}^{1}{\geq}N$ on the expanded execution model.

Efficiently calculating ${}^{1}{\geq}N$ is more challenging. The algorithm
below proceeds along the same idea as how we handled simple sentences, namely
by recursing over $E$. We first need to extend every node $u$ in $E$ with a
new label $u.{\sf A}$. This label is a set containing pairs of the form
$V\mapsto p$ where $V$ is a set of $M$’s states and $p$ is the probability
that $E@u$ would cover all the states mentioned in $V$. Only $V$’s whose
probability is non-zero need to be included in this mapping. After all nodes
in $E$ is labelled like this, the probability ${}^{1}{\geq}N$ can be
calculated from the ${\sf A}$ of the root node $u_{0}$:

$P(^{1}{\geq}N\;|\;\sigma)\ =\ \sum_{V{\mapsto}p\;\in\;u_{0}.A}{\bf
if}\;|V|\geq N\;{\bf then}\;p\;{\bf else}\;0$ (3)

The labelling is done recursively over $E$ as follows:

  1. 1.

The base case is the terminal node $\\#$. The $\sf A$ label of $\\#$ is just
$\emptyset$.

  2. 2.

For every node $u\in E$, we first recurse to all its successors. Then, we
calculate a preliminary mapping for $u$ in the following multi-set
$A^{\prime}$:

$A^{\prime}\ =\ \\{\;V{\cup}\\{u.{\sf st}\\}\mapsto
p{*}P_{E}(u{\rightarrow}v)\ |\ v\in u.{\sf next},\ V{\mapsto}p\in v.{\sf
A}\;\\}$

As a multi-set note that $A^{\prime}$ may contain duplicates, e.g. two
instances of $V\mapsto p_{0}$. Additionally, it may contain different maps
that belong to the same $V$, e.g. $V\mapsto p_{1}$ and $V\mapsto p_{2}$. All
these instances of $V$ need to be merged by summing up their $p$’s, e.g. the
above instances is to be merged to $V\mapsto p_{0}{+}p_{0}{+}p_{1}{+}p_{2}$
The function $\sf merge$ will do this. The label $u.{\sf A}$ is then just:
$u.{\sf A}\ =\ {\sf merge}(A^{\prime})$

The recursion terminates because $E$ is acyclic.

The above algorithm can however perform worse than a direct calculation via
Def. 3. The reason is that $\sf merge$ is an expensive operation if we do it
literally at every node. If we do not merge at all, and make the ${\sf A}$’s
multi-sets instead of sets, we will end up with $u_{0}.{\sf A}$ that contains
as many elements as the number of paths in $E$, so we are not better of
either. Effort to merge is well spent if it delivers large reduction in the
size of the resulting set, otherwise the effort is wasted. Unfortunately it is
hard to predict the amount of reduction we would get for each particular
merge. We use the following merge policy. We only merge on nodes at the
$B{-}1$-th position of ’bridges’ where $B$ is the length of the bridge at
hand. A bridge is a sequence of nodes $v_{0},...,v_{B{-}1}$ such that: (1)
every $v_{i}$ except the last one has only one outgoing edge, leading to
$v_{i{+}1}$, and (2) the last node $v_{B{-}1}$ should have more than one
successor. A bridge forms thus a deterministic section of $E$, that leads to a
non-deterministic section. Merging on a bridge is more likely to be cost
effective. Furthermore, only one merge is needed for an entire bridge. Merging
on a non-deterministic node (a node with multiple successors) is risky. This
policy takes a conservative approach by not merging at all on such nodes. The
next section will discuss the performance of our algorithm.

##  5 Experimental Results

In the following experiment we benchmark the algorithm from Section 4 against
the ’brute force’ way to calculate coverage using Def. 3. We will use a family
of models $M_{m}$ in Fig. 5. Despite its simplicity, $M_{m}$ is highly non-
deterministic and is designed to generate a large number of executions and
words.

We generate a family of execution models $E(i,m)$ by applying a test case
$tc^{i}$ on the model $M_{m}$ where $m\in\\{0,2,8\\}$. The test case is:

$tc^{i}\ =\ ac^{i}ab^{i}ac^{i}a$

The table in Fig. 6 (left) shows the statistics of all execution models used
in this experiment. Additionally we also construct $E(i,m)^{3}$ (applying
3-word expansion). The last column in the table shows the number of nodes in
the corresponding $E(i,m)^{3}$ (the number of executions stays the same, of
course).

0$M_{m}:$1234$t_{0}$…$t_{m{-}1}$5$a$$(0.3)c$$(0.7)c$$a$$(0.7)c$$(0.3)c$$b$$a$$(q)c$$(p)c$$(p)c$$a$$c$$c$
Figure 5:  The model $M_{m}$ used for the benchmarking. If $m{=0}$ then there
is no states $t_{i}$ and $q{=}1$. If $m{>}0$ then we have states
$t_{0}...t_{m{-}1}$; $p{=}0.3/m$ and $q=0.7$.

| $|tc|$ | $\\#nodes$ | $\\#paths$ | $\\#nodes^{3}$
---|---|---|---|---
$E(5,0)$ | 20 | 26 | 16 | 103(4)
$E(6,0)$ | 23 | 30 | 32 | 144(5)
$E(7,0)$ | 26 | 34 | 64 | 223(7)
$E(8,0)$ | 29 | 38 | 128 | 381(10)
$E(9,0)$ | 32 | 42 | 256 | 422(10)
$E(10,0)$ | 35 | 46 | 512 | 501(11)
$E(11,0)$ | 38 | 50 | 1024 | 659(13)
$E(12,0)$ | 41 | 54 | 2048 | 700(13)
$E(5,2)$ | 20 | 34 | 336 | 185(5)
$E(6,2)$ | 23 | 40 | 1376 | 306(8)
$E(7,2)$ | 26 | 46 | 5440 | 435(9)
$E(8,2)$ | 29 | 52 | 21888 | 695(13)
$E(9,2)$ | 32 | 58 | 87296 | 944(16)
$E(10,2)$ | 35 | 64 | 349696 | 1073(17)
$E(11,2)$ | 38 | 70 | 1397760 | 1333(19)
$E(12,2)$ | 41 | 76 | 5593088 | 1582(21)
$E(5,8)$ | 20 | 58 | 3600 | 863(15)
$E(6,8)$ | 23 | 70 | 29984 | 2760(39)
$E(7,8)$ | 26 | 82 | 175168 | 4287(52)
$E(8,8)$ | 29 | 94 | 1309824 | 8261(88)
$E(9,8)$ | 32 | 106 | 8225024 | 23726(224)

| $f_{1}$ | $f_{2}$ | $f_{3}$ | $f_{4}$
---|---|---|---|---
E(5,0) | 0.001 | 0.002 | 0.001 | 0.002
E(6,0) | 0.001 | 0.002 | 0.001 | 0.002
E(7,0) | 0.001 | 0.003 | 0.001 | 0.003
E(8,0) | 0.001 | 0.004 | 0.001 | 0.005
E(9,0) | 0.001 | 0.005 | 0.002 | 0.006
E(10,0) | 0.001 | 0.006 | 0.003 | 0.008
E(11,0) | 0.001 | 0.008 | 0.004 | 0.012
E(12,0) | 0.001 | 0.008 | 0.009 | 0.024
E(5,2) | 0.001 | 0.002 | 0.002 | 0.004
E(6,2) | 0.001 | 0.004 | 0.002 | 0.01
E(7,2) | 0.001 | 0.005 | 0.003 | 0.039
E(8,2) | 0.001 | 0.01 | 0.005 | 0.138
E(9,2) | 0.001 | 0.014 | 0.01 | 0.44
E(10,2) | 0.001 | 0.012 | 0.019 | 1.09
E(11,2) | 0.001 | 0.018 | 0.041 | 3.13
E(12,2) | 0.001 | 0.023 | 0.091 | 10.68
E(5,8) | 0.001 | 0.011 | 0.006 | 0.032
E(6,8) | 0.001 | 0.04 | 0.034 | 0.279
E(7,8) | 0.001 | 0.076 | 0.073 | 1.38
E(8,8) | 0.002 | 0.154 | 0.266 | 12.04
E(9,8) | 0.002 | 0.46 | 0.539 | 219

Figure 6: Left: the execution models used in the benchmark. $\\#nodes$ and
$\\#paths$ are the number of nodes and full paths (executions) in the
corresponding execution model; $\\#nodes^{3}$ is the number of nodes in the
resulting 3-word expansion model. The number between brackets is
$\\#nodes^{3}/\\#nodes$. Right: the run time (seconds) of our coverage
calculation algorithm on different execution models and coverage goals.

The number of possible executions in the execution models correspond to their
degree of non-determinism. The test case $tc^{i}$ has been designed as such
that increasing $i$ exponentially increases the non-determinism of the
corresponding execution model (we can see this in Figure 6 by comparing
$\\#paths$ with the $i$ index of the corresponding $E(i,m)$).

All the models used ($M_{0}$, $M_{2}$, and $M_{8}$) are non-deterministic:
$M_{0}$ is the least non-deterministic one whereas $M_{8}$ is very non-
deterministic. This is reflected in the number of possible executions in their
corresponding execution models, with $E(i,8)$ having far more possible
executions than $E(i,0)$.

The following four coverage goals are used:

$\scriptsize\begin{array}[]{|l|c|c|}\hline\cr\mbox{goal}&\mbox{type}&\mbox{word
expansion}\\\ \hline\cr f_{1}:\langle 2\rangle;\langle
t_{0}\rangle&\mbox{simple sentence}&\mbox{no}\\\ \hline\cr f_{2}:\langle
1,1,1\rangle;\langle 4,4,4\rangle&\mbox{non-simple sentence}&\mbox{3-word}\\\
\hline\cr f_{3}:\;^{1}{\geq}_{8}&\mbox{aggregate}&\mbox{no}\\\ \hline\cr
f_{4}:\;^{3}{\geq}_{8}&\mbox{aggregate}&\mbox{3-word}\\\ \hline\cr\end{array}$

We let our algorithm calculates the coverage of each of the above goals on the
execution models $E(5,0)...E(9,8)$ and measure the time it takes to finish the
calculation. For the merging policy, $n$ is set to 1 when the goal does not
need word expansion, and else it is set to be equal to the expansion
parameter. The experiment is run on a Macbook Pro with 2,7 GHz Intel i5 and 8
GB RAM. Fig. 6 (right) shows the results. For example, we can see that $f_{1}$
can be calculated in just a few milli seconds, even on $E(12,m)$ and $E(i,8)$.
In contrast, brute force calculation using Def. 3 on e.g.
$E(11,2),E(12,2),E(8,8)$, and $E(9,8)$ would be very expensive, because it has
to quantify over more than a million paths in each of these models.

|
---|---
|
Figure 7:  The graphs show our algorithm’s speedup with respect to the brute
force calculation on four different goals: $f_{1}$ (top left), $f_{2}$ (top
right), $f_{3}$ (bottom left), and $f_{4}$ (bottom right). $f_{1}$ and $f_{2}$
are non-aggregate, whereas $f_{3}$ and $f_{4}$ are aggregate goals.
Calculating $f_{1}$ and $f_{3}$ does not use word expansion, whereas $f_{2}$
and $f_{4}$ require 3-word expansion. Each graph shows the speedup with
respect to three families of execution models: $E(i,0)$, $E(i,2)$, and
$E(i,8)$. These models have increasing degree of non-determinism, with models
from $E(i,8)$ being the most non-deterministic ones compared to the models
from other families (with the same $i$). The horizontal axes represent the $i$
parameter, which linearly influences the length of the used test case. The
vertical axes show the speedup in the logarithmic scale.

Fig.7 shows the speedup of our algorithm with respect to the brute force
calculation —note that the graphs are set in logarithmic scale. We can see
that in almost all cases the speedup grows exponentially with respect to the
length of the test case, although the growth rate is different in different
situations. We can notice that the speed up on $E(i,0)$ is much lower (though
we still have speedup, except for $f_{4}$ which we will discuss below). This
is because $E(i,0)$’s are not too non-deterministic. They all induce less than
2100 possible executions. The brute force approach can easily handle such
volume. Despite the low speedup, on all $E(i,0)$’s our algorithm can do the
task in just few milli seconds (1 - 24 ms).

The calculation of $f_{1}$ is very fast (less than 2 ms). This is expected,
because $f_{1}$ is a simple sentence. The calculation of $f_{2}$, on the other
hand, which is a non-simple sentence, must be executed on the corresponding
$3$-word expanded execution model, which can be much larger than the original
execution model. E.g. $E(9,8)^{3}$ is over 200 times larger (in the number of
nodes) than $E(9,8)$. Despite this we see the algorithm performs pretty well
on $f_{2}$.

$f_{3}$ and $f_{4}$ are both aggregate goals. The calculation of $f_{3}$ is
not problematical, however we see that $f_{4}$ becomes expensive on the models
$E(12,2),E(8,8)$, and $E(9,8)$ (see Fig. 6 right). In fact, on $E(9,8)$ the
calculation of $f_{4}$ is even worse than brute force (the dip in the red line
in Fig. 7). Recall that $f_{4}=\;^{3}{\geq}_{8}$; so, calculating its coverage
requires us to sum over different sets of words of size 3 that the different
executions can generate. $E(12,2),E(8,8)$, and $E(9,8)$ are large (over 70
states) and highly non-deterministic. Inevitably, they generate a lot of words
of size 3, and therefore the number of possible sets of these words explodes.
E.g. on $E(8,8)$ and $E(9,8)$ our algorithm ends up with about 1.2M an 6.7M
sets of words to sum over. In contrast, the number of full paths in these
models are about respectively 1.3M and 8.2M. At this ratio, there is not much
to gain with respect to the brute force approach that simply sums over all
full paths, whereas our algorithm also has to deal with the overhead of book
keeping and merging. Hypothetically, if we always merge, the number of final
sets of words can be reduced to respectively about 500K and 2M, so summing
over them would be faster. We should not do this though, because merging is
expensive, but the numbers do suggest that there is room for improvement if
one can figure out how to merge more smartly.

##  6 Related Work

To the best of our knowledge the concept of probabilistic coverage has not
been well addressed in the literature on non-deterministic MBT, or even in the
literature on probabilistic automata. A paper by Zu, Hall, and May [30] that
provides a comprehensive discussion on various coverage criteria does not
mention the concept either. This is a bit surprising since coverage is a
concept that is quite central in software testing. We do find its mentioning
in literature on statistical testing, e.g. [8, 29]. In [29] Whittaker and
Thomason discussed the use of Markov chains to encode probabilistic behavioral
models. The probabilities are used to model the usage pattern of the SUT. This
allows us to generate test sequences whose distribution follows the usage
pattern (so-called ’statistical testing’). Techniques from Markov chain are
then used to predict properties of the test sequences if we are to generate
them in this way, e.g. the probability to obtain a certain level of node or
edge coverage, or conversely the expected number of test runs needed to get
that level of coverage. In contrast, in our work probabilities are used to
model SUT’s non-determinism, rather than its usage pattern. We do not concern
ourselves with how the tester generates the test sequences, and focuses purely
on the calculation of coverage under the SUT’s non-determinism. Our coverage
goal expressions are more general than [29] by allowing words of arbitrary
length (rather than just words of length one or two, which would represent
state and respectively edge coverage), clauses, and sentences to be specified
as coverage goals. Coverage calculation in both [8, 29] basically comes down
to the brute force calculation in Def. 3.

Our algorithm to calculate the coverage of simple sentences has some
similarity with the probabilistic model checking algorithm for Probabilistic
Computation Tree Logic (PCTL) [10, 15]. Although given a formula $f$ a model
checking algorithm tries to decide whether or not $f$ is valid on the given
behavior model, the underlying probabilistic algorithm also labels for every
state in the model with the probability that any execution that starts from
that state would satisfy $f$. Since we only need to calculate over execution
models, which are acyclic, there is no need to do a fix point iteration as in
[15]. From this perspective, our algorithm can be seen as an instance of [15].
However we also add $k$-word expansion. In addition to simplifying the
algorithm when dealing with non-simple sentences, the expansion also serves as
a form of memoisation (we do not have to keep calculating the probability for
a state $u$ to lead to a word $w$). In particular the calculation of aggregate
coverage goals benefits from this memoisation. Though, the biggest difference
between our approach with a model checking algorithm is that the latter does
not deal with aggregate properties (there is no concept of aggregate formulas
in PCTL). Our contribution can also be seen as opening a way to extend a
probabilistic model checking algorithm to calculate such properties. We
believe it is also possible to generalize over the aggregation so that the
same algorithm can be used to aggregate arbitrary state attributes that admit
some aggregation operator (e.g. the cost of staying in various states, which
can be aggregated with the ’+’ operator).

In this paper we have focused on coverage analyses. There are other analyses
that are useful to mention. In this paper we abstract away from the data that
may have been exchanged during the interactions with the SUT. In practice many
systems do exchange data. In this situation we may also want to do data-
related analyses as well. E.g. the work by Prasetya [19] discussed the use of
an extended LTL to query temporal relations between the data exchanged through
the test sequences in a test suite. This is useful e.g. to find test sequences
of a specific property, or to check if a certain temporal scenario has been
covered. The setup is non-probabilistic though (a query can only tell whether
a temporal property holds or not), so an extension would be needed if we are
interested in probabilistic judgement. Another example of analyses is risk
analyses as in the work by Stoelinga and Timmer [23]. When testing a non-
deterministic system, we need to keep in mind that although executing a test
suite may report no error, there might still be lurking errors that were not
triggered due to internal non-determinism. Stoelinga and Timmer propose to
annotate each transition in a model with the estimated probability that it is
incorrectly implemented and the entailed cost if the incorrect behavior
emerges111We gloss over the complication that the transition might be in a
cycle. A test case may thus exercise it multiple times. Each time, exercising
it successfully would arguably decrease the probability that it still hides
some hidden erroneous behavior. This requires a more elaborate treatment, see
[23] for more details. . This then allows us to calculate the probability that
a successful execution of a test suite still hides errors, and the expected
cost (risk) of these hidden errors.

##  7 Conclusion

We have presented a concept of probabilistic coverage that is useful to
express the coverage of a test suite in model-based testing when the used
model is non-deterministic, but has been annotated with estimation on the
probability of each non-deterministic choice. Both aggregate and non-aggregate
coverage goals can be expressed, and we have presented an algorithm to
efficiently calculate the probabilistic coverage of such goals. Quite
sophisticated coverage goals can be expressed, e.g. sequence (words) coverage
and sequence of sequences (sentences) coverage. We have shown that in most
cases the algorithm is very efficient. A challenge still lies on calculating
aggregate $k$-wise test goals on test cases that repeatedly trigger highly
non-deterministic parts of the model. Such a situation bounds to generate
combinatoric explosion on the possible combinations of words that need to be
taken into account. Beyond a certain point, the explosion becomes too much for
the merging policy used in our algorithm to handle. Analyses on the data
obtained from our benchmarking suggests that in theory there is indeed room
for improvement, though it is not yet clear what the best course to proceed.
This is left for future work.

## References

  * [1] Ammann, P., Offutt, J.: Introduction to software testing. Cambridge University Press (2016)
  * [2] Arnold, A.: Finite transition systems. international series in computer science (1994)
  * [3] Baier, C., Katoen, J.P., Larsen, K.G.: Principles of model checking. MIT press (2008)
  * [4] Belinfante, A.: JTorX: Exploring Model-Based Testing. Ph.D. thesis, University of Twente (2014)
  * [5] Bringmann, E., Krämer, A.: Model-based testing of automotive systems. In: 2008 1st international conference on software testing, verification, and validation. pp. 485–493. IEEE (2008)
  * [6] Craggs, I., Sardis, M., Heuillard, T.: Agedis case studies: Model-based testing in industry. In: Proc. 1st Eur. Conf. on Model Driven Software Engineering. pp. 129–132 (2003)
  * [7] Dallmeier, V., Lindig, C., Wasylkowski, A., Zeller, A.: Mining object behavior with adabu. In: Proceedings of the International Workshop on Dynamic Systems Analysis (WODA). pp. 17–24. ACM (2006). https://doi.org/10.1145/1138912.1138918, http://doi.acm.org/10.1145/1138912.1138918
  * [8] Denise, A., Gaudel, M.C., Gouraud, S.D.: A generic method for statistical testing. In: 15th Int. Symp. on Software Reliability Engineering ISSRE. pp. 25–34. IEEE (2004)
  * [9] Grindal, M., Offutt, J., Andler, S.F.: Combination testing strategies: a survey. Software Testing, Verification and Reliability 15(3), 167–199 (2005)
  * [10] Hansson, H., Jonsson, B.: A logic for reasoning about time and reliability. Formal aspects of computing 6(5), 512–535 (1994)
  * [11] Heerink, L., Feenstra, J., Tretmans, J.: Formal test automation: The conference protocol with phact. In: Testing of Communicating Systems, pp. 211–220. Springer (2000)
  * [12] Hoare, C.A.R.: Communicating sequential processes. Prentice Hall (2004)
  * [13] Jard, C., Jéron, T.: Tgv: theory, principles and algorithms. International Journal on Software Tools for Technology Transfer 7(4), 297–315 (2005)
  * [14] Kanstrén, T., Puolitaival, O.P.: Using built-in domain-specific modeling support to guide model-based test generation. Model-Driven Engineering of Information Systems: Principles, Techniques, and Practice pp. 295–319 (2012)
  * [15] Kwiatkowska, M., Norman, G., Parker, D.: Stochastic model checking. In: International School on Formal Methods for the Design of Computer, Communication and Software Systems. pp. 220–270. Springer (2007)
  * [16] Nachmanson, L., Veanes, M., Schulte, W., Tillmann, N., Grieskamp, W.: Optimal strategies for testing nondeterministic systems. In: ACM SIGSOFT Software Engineering Notes. vol. 29, pp. 55–64. ACM (2004)
  * [17] Peleska, J.: Industrial-strength model-based testing - state of the art and current challenges. In: Proceedings 8th Workshop on Model-Based Testing (MBT). pp. 3–28 (2013). https://doi.org/10.4204/EPTCS.111.1, https://doi.org/10.4204/EPTCS.111.1
  * [18] Petke, J., Cohen, M.B., Harman, M., Yoo, S.: Practical combinatorial interaction testing: Empirical findings on efficiency and early fault detection. IEEE Transactions on Software Engineering 41(9), 901–924 (2015)
  * [19] Prasetya, I.: Temporal algebraic query of test sequences. Journal of Systems and Software 136, 223–236 (2018)
  * [20] Pretschner, A., Philipps, J.: 10 methodological issues in model-based testing. In: Model-based testing of reactive systems, pp. 281–291. Springer (2005)
  * [21] Schur, M., Roth, A., Zeller, A.: Mining behavior models from enterprise web applications. In: Proceedings of the 9th Joint Meeting on Foundations of Software Engineering. pp. 422–432. ACM (2013). https://doi.org/10.1145/2491411.2491426, http://doi.acm.org/10.1145/2491411.2491426
  * [22] Stoelinga, M.: An introduction to probabilistic automata. Bulletin of the EATCS 78(2), 176–198 (2002)
  * [23] Stoelinga, M., Timmer, M.: Interpreting a successful testing process: risk and actual coverage. In: 3rd Int. Symp. on Theoretical Aspects of Software Engineering TASE. pp. 251–258. IEEE (2009)
  * [24] Tervoort, T., Prasetya, I.: APSL: A light weight testing tool for protocols with complex messages. In: Haifa Verification Conference. pp. 241–244. Springer, (2017)
  * [25] Tretmans, G.J.: A formal approach to conformance testing. Ph.D. thesis, Twente Univ. (1992)
  * [26] Tretmans, J., Brinksma, E.: TorX: Automated model-based testing. In: 1ST European Conf. on Model-Driven Software Engineering (2003)
  * [27] Utting, M., Pretschner, A., Legeard, B.: A taxonomy of model-based testing approaches. Software Testing, Verification and Reliability 22(5), 297–312 (2012)
  * [28] Vos, T., Tonella, P., Prasetya, W., Kruse, P.M., Bagnato, A., Harman, M., Shehory, O.: Fittest: A new continuous and automated testing process for future internet applications. In: Software Maintenance, Reengineering and Reverse Engineering (CSMR-WCRE), 2014 Software Evolution Week-IEEE Conference on. pp. 407–410. IEEE (2014)
  * [29] Whittaker, J.A., Thomason, M.G.: A markov chain model for statistical software testing. IEEE Transactions on Software Engineering 20(10), 812–824 (1994)
  * [30] Zhu, H., Hall, P.A., May, J.H.: Software unit test coverage and adequacy. ACM Computing Surveys 29(4), 366–427 (1997)