# Image Compression using only Attention based Neural Networks

Natacha Luka 1<EMAIL_ADDRESS>Thanks to Agence Innovation Défense (AID)
and École des Ponts ParisTech for funding. This work was granted access to the
HPC resources of IDRIS under the allocation 20XX-AD011011290 made by GENCI
Romain Negrel1<EMAIL_ADDRESS>David Picard1<EMAIL_ADDRESS>1LIGM, École des Ponts, Université Gustave Eiffel, ESIEE Paris  CNRS, Marne-
La-Vallée, France

###### Abstract

In recent research, Learned Image Compression has gained prominence for its
capacity to outperform traditional handcrafted pipelines, especially at low
bit-rates. While existing methods incorporate convolutional priors with
occasional attention blocks to address long-range dependencies, recent
advances in computer vision advocate for a transformative shift towards fully
transformer-based architectures grounded in the attention mechanism. This
paper investigates the feasibility of image compression exclusively using
attention layers within our novel model, _QPressFormer_. We introduce the
concept of learned _image queries_ to aggregate patch information via cross-
attention, followed by quantization and coding techniques. Through extensive
evaluations, our work demonstrates competitive performance achieved by
convolution-free architectures across the popular Kodak, DIV2K, and CLIC
datasets.

Keywords — Learned Image Compression, Vision Transformers, Visual Attention

Our:
---
bpp=0.31
PSNR=23.42
lpips=0.2402
Conv+MSE:
---
bpp=0.31
PSNR=27.48
lpips=0.4023
Original BPG:
---
bpp=0.32
PSNR=29.31
lpips=0.3716
Our:
---
bpp=0.31
PSNR=23.48
lpips=0.2391
Conv+MSE:
---
bpp=0.31
PSNR=27.48
lpips=0.4023
Original BPG:
---
bpp=0.32
PSNR=29.31
lpips=0.3716
Figure 1: Image reconstruction from DIV2K. Zoom comparisons with original,
convolutional equivalent model from [1] and BPG are shown on the sides. Even
at lower PSNR, our method significantly improves on the details of the
penguin’s fur, as captured by the lpips metrics. Examples on Kodak and CLIC in
appendix A.

##  1 Introduction

Image compression stands as a pivotal research domain, driven by the ubiquity
of images in today’s digital landscape. The persistent pursuit of higher image
quality, while simultaneously optimizing bit size for storage and
transmission, underscores its paramount importance. For many years, image
compression codecs were meticulously handcrafted, yielding iconic solutions
like the widely-used JPEG algorithm. In the post-JPEG era, a lineage of
handcrafted codecs emerged, including JPEG2000 [2], BPG format [3] (leveraging
the HEVC video compression standard [4]), and Google’s WebP [5]. Nevertheless,
recent years have witnessed a transformative shift towards learned image
compression, propelled by deep learning approaches. In the contemporary
landscape, these learned codecs have not only achieved parity with their
handcrafted counterparts but have also demonstrated the potential to surpass
them.

In both hand-designed and learned image compression paradigms, common
components are shared within the encoding and decoding processes. During
encoding, an analysis transform maps the image from its pixel space
representation to a more compression-friendly domain, followed by
quantization, a pivotal step in lossy codecs. Subsequently, a lossless coding
scheme translates the image data into a compact bitstream. In the decoding
phase, inverse operations unfold: the received bits are decoded, followed by
synthesis through a suitable transform, ultimately restoring the image to its
pixel space. A distinctive advantage inherent in learned compression methods
arises from their holistic optimization approach. In this paradigm, all
constituent blocks are collectively learned and end-to-end optimized. This
optimization strategy diligently works to minimize distortion while
concurrently minimizing bit-rate consumption. Conversely, traditional hand-
designed methods require the separate optimization of individual blocks, which
may not inherently synergize when used together. Furthermore, hand-designed
approaches involve human-centric decisions in selecting the transform,
introducing inherent bias. In stark contrast, learned methods offer greater
flexibility, as the network autonomously adapts both the analysis and
synthesis transforms, determining the optimal compromise between rate and
quality within specified constraints.

In learned image compression, numerous approaches have emerged for the
analysis stage. Remarkably, transformer-based architectures, introduced by
Vaswani et al. [6], have remained relatively unexplored, despite their
pervasive adoption in Computer Vision where they have demonstrated their
supremacy over CNNs in many computer vision tasks [7]. It is interesting to
note that while attention or transformer blocks find occasional integration,
they are almost invariably utilized in conjunction with other usual
convolution blocks, in both the analysis and synthesis processes [8, 9].
Furthermore, a critical prerequisite for effective image compression is the
reduction of spatial redundancy, and the transformer’s attention mechanism,
which aggregates related elements within a sequence, appears especially well-
suited for this purpose. With this context in mind, we propose in this paper a
comprehensive investigation of the potential benefits of relying exclusively
on the transformer’s attention mechanism for image compression, thereby
eliminating convolutional priors from the equation.

The remaining of this paper is structured as follows: First, we present recent
achievements in learned image compression, as well as basics on attention
mechanisms. Then, we present and explain our transformer-based
analysis/synthesis model based on our _image queries_ concept. These _image
queries_ are learned prototypes which aggregate information from image step by
step thanks to transformer cross-attention blocks. Finally, we show and
analyze results of the proposed architecture on several datatsets to give some
insights about the inner working of our proposed model _QPressFormer_ , before
we conclude.

##  2 Related Work

###  2.1 Learned Image Compression

In the realm of image compression, conventional lossless and lossy codecs,
typified by ubiquitous standards like JPEG and its successors such as JPEG2000
[2], WebP [5], and BPG [3] , have long prevailed as the industry norm. These
well-established codecs continue to be indispensable in contemporary
applications, exemplified by the ongoing development of the VVC codec [10] for
video, with adaptions tailored for image compression. Nonetheless, the recent
years have witnessed a burgeoning interest in deep learning-driven approaches
for both lossless and near-lossless image compression [11], marking a
noteworthy departure from conventional methodologies.

Todericci et al. introduced an innovative approach using recurrent neural
networks (RNNs) [12, 13] for variable-rate image compression. Their method
involves multiple passes of the original image through the network, generating
an $m$-bit representation at each step to control the final bit-rate. While
they explored various scenarios, they did not incorporate attention
mechanisms.

Another strategy entails training multiple networks for distinct bit-rates.
Ballé et al. introduced an innovative approach [14, 1] centered on learning
the probability distribution of images within a latent space, with a
concurrent focus on minimizing both the distortion and the entropy of this
learned distribution. As per Shannon’s theory, entropy serves as a lower bound
for code length, rendering it synonymous with bit-rate minimization. The
delicate equilibrium between these two objectives delineates the
distortion/bit-rate trade-off. Their methodology incorporates convolutional
transforms, paired with normalization functions [15] for image compression. To
enhance the estimation of the latent probability distribution, Ballé and
colleagues introduced side information through latent variable models, where
the variables obtained from the analysis transform are treated as observed
variables. A prior, referred to as the hyper-prior, is learned on these
variables, effectively serving as a ’prior of the prior’ to better account for
spatial dependencies within the latent space. The ensuing hyper-latent
variables are subsequently learned, contributing to the distribution in the
hyper-latent space whose entropy is minimized to minimize bit-rate of this
side-information. These hyper-variables allowing the variances inference of
the original latent variables modeled as zero-mean Gaussian variables.

In recent works [16], significant advancements were made in optimizing the
image compression pipeline, with two notable steps. Causal consideration of
context began to be incorporated to enhance the inference of parameters within
the distribution of subsequent latent variables. Building upon this
foundation, attention mechanisms were integrated to encompass contextual
information, wherein latent variables were modeled as Gaussian mixture models
(GMMs), with parameter inferences facilitated by hyper latent variables [17].
Concurrently, several works embraced attention mechanisms to more effectively
account for contextual nuances, particularly concerning bit allocation. These
mechanisms dynamically allocate more bits to regions of higher attention,
signifying their increased complexity and importance [18, 19]. Notably, while
transformers have occasionally served as attention modules within these
architectures [9, 8], they have been employed alongside convolutional
components, distinguishing our approach.

Finally, generative models have also been investigated [20, 21, 22]. This
results in better perception quality image especially in low bit-rate setup
even if it is sometimes not reflected by distortion metrics like PSNR or MS-
SSIM due to the Perception-Distortion trade-off [23]. Similarly to other
approaches, even when attention is used it is in addition to convolutions.

###  2.2 Transformers

Transformers, initially introduced by Vaswani et al. in the domain of Natural
Language Processing (NLP) for translation tasks [6], have since permeated
diverse NLP applications, from BERT [24] to LlaMA [25] models. Within Computer
Vision (CV), Dosovitskiy et al. inaugurated the integration of transformers
[26], which are now prevalent in various vision tasks encompassing object
detection [27], image synthesis [28], or cross-modal endeavors such as Text-
Image Retrieval and Visual Question Answering [29]. Notably, transformers have
supplanted traditional CNNs [7].

The transformer architecture relies on an attention mechanism, facilitating
interactions between two distinct sequences of tokens. Its fundamental
objective is to iteratively aggregate tokens from a sequence denoted as
_values_ , organized in a matrix $V$, into another sequence known as _queries_
, stored in a matrix $Q$. This aggregation is guided by the similarity between
elements in $Q$ and $V$. To achieve this, the _values_ are subject to $h$
projections, distributed across different sub-units termed _heads_ , yielding
$h$ sets of _keys_ , amalgamated within a matrix $K$. For each set of _keys_ ,
a similarity score $s_{i,j}$ is computed using the dot product between $q_{i}$
and $k_{j}$, where $k_{j}$ corresponds to specific features of $v_{j}$. These
similarity scores are gathered within a matrix, normalized through the
$\mathrm{softmax}$ function, and designated as the attention matrix $A$. Thus,
at each iteration $n$, $h$ attention matrices are computed using the following
equation:

$\displaystyle A^{n}$
$\displaystyle=\sigma\left(Q^{n}\left(K^{n}\right)^{T}\right)$ (1)

Once the attention matrices are obtained, they are employed to compute a
weighted sum of the _values_ , subsequently added to the _queries_.
Consequently, the _queries_ accumulate the most akin _values_ corresponding to
$h$ specific features, as determined through the projection of _values_ into
$h$ sets of _keys_. This process is concisely captured in the following
equation, where $[\dots]_{h}$ denotes the aggregation of the $h$ sub-units:

$\displaystyle Q^{n+1}$
$\displaystyle=Q^{n}+\left[A^{n}V^{n}\right]_{h}=Q^{n}+\left[\left(\sum_{k}a^{n}_{ik}v^{n}_{kj}\right)_{i,j}\right]_{h}$
(2)

The equation 2 correspond to the main module of transformer called Multi-Head-
Attention denoted $\operatorname{MHA}$ in the following. In order to preserve
the dimension $d$ of the whole transformer, all elements in the $h$ sub-units
are projected in a space of dimension $d_{h}=\frac{d}{h}$, which implies
$d\equiv 0\left(mod\ h\right)$.

There are 2 different $\operatorname{MHA}$ modules. The Self-Attention module,
$\operatorname{SA}\left[Q\right]$, correspond to the case where each _query_
aggregates itself with its most similar _queries_ according to particular
features. In this case _queries_ and _values_ are identical. The other case
correspond to the Cross-Attention , $\operatorname{CA}\left[Q,V\right]$. Here,
the set of _queries_ and _values_ are different. The sequence of _queries_
aggregates information from another sequence. More details are available in
[6].

##  3 Method

In this section, we introduce our convolution-free architecture, denoted
QPressFormer, which exclusively relies on attention blocks. For the
compression component, we adopt a factorized prior model similar to [1]. To
ensure comparability, our aim is to develop a model with similar target
bitrate, while replacing convolutional analysis and synthesis transforms with
attention-based counterparts. The whole pipeline of our QPressFormer model is
illustrated in Figure 2.

Analysis
---
$\varphi$
---

$I$

$E$
---
$N$ Learnable Images Queries
---

$\varphi\left(I\right)$

$V_{E}$
---
$K_{E}$
---
$Q_{E}^{in}$
---
$Q^{out}_{E}$
---
Synthesis
---

$\hat{I}$

$D$
---
Learnable Patches
---
$Q_{D}^{in}$
---
$K_{D}$
---
$V_{D}$
---
$\psi$
---
$Q_{D}^{out}$
---
Quantization
---
Coding
---
Decoding
---

1010…1110

Coding/Decoding Module
---
Bitstream
---
Original PatchesLearned parametersOutpout queries: Latent image
representationRecons Patches Figure 2: QPressFormer : Our Proposed
Architecture.

An image $I$ is first divided into patches, which are unrolled and augmented
with positional encodings corresponding to their patch index using the
transformation denoted as $\varphi$. The resulting set of vectors serves as
the _values_ /_keys_ ($V_{E}$/$K_{E}$) for the encoder $E$, which employs
attention blocks exclusively. Additionally, a set of $N$ learned vectors,
referred to as ’Image queries’ ($Q_{E}^{in}$), serves as _queries_. The set
$Q_{E}$ aggregates information from the image at each step $n$ through the
attention process outlined in Section 2.2. In-depth, the blocks of $E$ consist
of transformer decoder blocks [6], governed by the following equations for
each layer $n$:

$\displaystyle Q^{n+1/2}$
$\displaystyle=\operatorname{CA}\left[\operatorname{SA}\left[Q^{n}\right],\varphi(I)\right]$
(3) $\displaystyle Q^{n+1}$
$\displaystyle=Q^{n+1/2}+\mathrm{FFW}^{n}(Q^{n+1/2})$

At the final stage, the output _queries_ ($Q_{E}^{out}$) containing the
compressed image. In this study, we work with 256x256 RGB images divided into
16x16 patches. This results in 256 patches, flattened as 768-dimensional
vectors, which sets the embedded dimension of our architecture. We opt for
$N=64$ _Image Queries_. $Q_{E}^{out}$ corresponds to the latent variables so
the probability distribution of $Q_{E}^{out}$ is estimated using a network and
its entropy minimized to reduce the bit-rate, as discussed in Section 2.1.
Importantly, our reconstruction perception metric is not the Mean-Square-Error
(MSE) but lpips [30], a perceptual metric that aligns better with human
perception.

The decoding process mirrors the encoding procedure. The vectors
$\left\\{Q_{E}^{out}\right\\}$, containing the compressed image, serve as
_values_ ($V_{D}$)/_keys_ ($K_{D}$) for the image decoder $D$. The _queries_
in $D$ are learned vectors corresponding to patch prototypes. They aggregate
the compressed image progressively, and in the final stage, are reshaped and
reassembled to yield the decompressed image $\hat{I}$.

The architecture is trained on the ImageNet using Adam with a learning rate of
1e-4. Beyond 100,000 steps, optimizing lpips becomes more challenging. We
observed that upscaling images to 512x512 before using lpips enhances results
without requiring other changes. Both $E$ and $D$ have 12 heads and a depth of
12.

##  4 Results

###  4.1 Quantitative Results

In Table 1, we present a comparative analysis between our model, BPG and
equivalent convolutional factorized prior models from [1] with pretrained
weights obtained from CompressAI [31] on three datasets (Kodak [32], CLIC [33]
and DIV2K [34]). For the MSE and MS-SSIM optimized models, we select weights
corresponding to the closest bit per pixel (bpp) to our model. Our evaluation
employs a range of metrics, including PSNR, MS-SSIM [35], lpips [30], and FID
[36], which assesses the realism of generated images. Given that our
architecture was trained on 256x256 images, we decompose all images into
256x256 segments and input each 256x256 part to the network before
reassembling them. In cases where the image dimensions are not divisible by
256x256, we employ center cropping to obtain the largest segment divisible by
256x256. To ensure a fair comparison with convolutional counterparts or hand-
designed BPG format, we apply the same image decomposition for reconstruction.

As shown in Table 1, the convolutional model optimized on MSE (respectively on
MS-SSIM) performs better in terms of PSNR (respectively on MS-SSIM). In the
same manner, our model performs significantly better on lpips. Moreover, our
model performs better on FID even if it was not trained specifically for it.
Besides, even if the PSNR and MS-SSIM of our model is lower than
convolutionnal counterpart, it captures better perceptual details – as
expected with the lpips optimization – as shown in figure 1. Convolution bias
is thus not a prerequisite for Image compression and a transformer
architecture seems at least equivalent or better – in terms of perceptual
metrics FID and lpips – than convolutional networks and BPG.

Method | bpp | PSNR$\uparrow$ | MS-SSIM$\uparrow$ | lpips$\downarrow$ | FID$\downarrow$
---|---|---|---|---|---
|  | Kodak
QpressFormer | 0.296 | 27.19 | 0.9203 | 0.2126 | 68.60
Conv+mse [1] | 0.258 | 31.59 | 0.9608 | 0.3386 | 154.06
Conv+mssim [1] | 0.320 | 30.66 | 0.9745 | 0.2913 | 117.94
BPG [3] | 0.301 | 32.78 | 0.9626 | 0.2872 | 110.90
|  | DIV2K (val-set)
QpressFormer | 0.300 | 25.50 | 0.9159 | 0.2267 | 19.19
Conv+mse [1] | 0.290 | 30.10 | 0.9642 | 0.3429 | 63.12
Conv+mssim [1] | 0.322 | 28.94 | 0.9716 | 0.3237 | 49.78
BPG [3] | 0.307 | 30.84 | 0.9604 | 0.3155 | 54.73
|  | CLIC (test-set)
QpressFormer | 0.286 | 27.96 | 0.9339 | 0.2244 | 9.63
Conv+mse [1] | 0.230 | 32.33 | 0.9657 | 0.3820 | 51.85
Conv+mssim [1] | 0.277 | 31.32 | 0.9758 | 0.3536 | 36.83
BPG [3] | 0.297 | 34.31 | 0.9712 | 0.3216 | 40.78
Table 1: Quantitative results on Kodak, CLIC and DIV2K datasets.

###  4.2 Visual analysis

To understand how the attention blocks process the image during encoding, we
perform several visualizations. During the attention process, the attention
weights are computed for a given pair of _query_ and _key_. As each _key_
corresponds to a patch, we can see which part of the image one _query_ pays
attention to by plotting its attention matrix. By combining the attention
matrices of all queries, layers and heads (with maximum), we obtain a heatmap
showing which parts of the image have been processed (Figure 3). As expected
for good reconstruction, the whole image is attended, almost uniformly. In
other words, at least one query pays attention only once in each image area.

Low attn(0.015) High attn(1.000)

Figure 3: Attention map (max attn scores for all queries/layers/heads) for a
Kodak image.

To better understand what information is captured by a _query_ , we visualize
the layers and heads mean attention map for each _query_ individually. An
example is given in Figure 9. It is worth noting the maximum of attention maps
and the mean of attention maps able us to make two distinct observations. The
mean correspond to the area which is the most watched by the query whereas the
maximum allows us to ensure that all areas are watched by queries at least one
time. An original 256x256 crop from Kodak and the decoded patch are given in
Figure 9 and 9. The figure 9 shows the mean attention map – across heads and
layers – for query 11\. The attention is more focused on a part of the roof.
To verify that this specific query encodes this particular region, we
reconstruct the image without it. We can then see on the figure 9 that the
area where the removed query paid attention to is not reconstructed properly.
Interestingly, the area is inferred and in-painted with the roof texture,
indicating that information from others patches able the encoder to infer
information in that area.

[t]0.19  [t]0.19  [t]0.19  Low(0.005) High(0.251) [t]0.19  [t]0.19

Figure 4: 2562 patch from Kodak Figure 5: Reconstruction Figure 6: Mean Attn
of $\left\\{Q_{E}\right\\}_{11}$ Figure 7: Reconstruction† Figure 8: Mean
recons. error† Figure 9: 256x256 Kodak Image Crop Reconstruction with all
$Q_{E}$ or with $\left\\{Q_{E}\right\\}_{11}$ removed (denoted †).

In Figure 9, the reconstruction error without this particular query and
averaged over the Kodak dataset shows that the spatial information encode by
the _query_ is the same for all images. In fact, after looking at the 64
averaged reconstruction error, nearly all _queries_ encode spatial
information. It is worth noting that this spatial encoding property in the
queries is an emergent behavior of the model, as it was not designed in the
architecture or in the optimization.

However, spatial information is not the whole information encoded by queries
$Q_{E}^{out}$. First, a look at maximum attention maps of queries allow us to
discover that the same query pays attention to areas in different locations
for different 256x256 crops of the same kodak image. Figure 27 in appendix C
gives more details. Second, if each query was specialized spacially, in-
painting wouldn’t be possible.

To better visualize the elements encoded by queries, we perform a PCA on the
whole queries for a Kodak image to reduce the dimensions of queries from 768
to 3. In the following we refers to these reduce dimensional queries as meta-
queries. We can then visualize a RGB image corresponding to the main, or at
least the most important, information decoded at each layer of the decoder.
Remember that the output queries of the encoder $Q_{E}^{out}$ serve as
keys/values for the decoder wherin the queries are learned patchs prototypes.
After reducing the dimensionality of $Q_{E}^{out}$, we can project on it the
mean across heads of cross-attention matrices (instead of the decoder values
obtained from the full encoder queries) for a given layer of the decoder with
the following equation:

$\displaystyle I_{visu}^{n}$ $\displaystyle=\frac{1}{nb_{heads}}\sum_{h\in
heads}\left(A_{h}^{n}\right)Q_{E;\left[:3\right]}^{out}$ (4)

In equation 4 $Q_{E;\left[:3\right]}^{out}$ denotes the meta-queries – i.e.
output encoder queries obtained after PCA for reducing dimension to 3. The
figure 14 shows clearly that queries contain spatial information extracted at
some specific layers: figure 14 shows the same colors across different
positions for different 256x256 blocks of an kodak image in addition to smooth
transitioning between colors. In contrast figure 14 show a nearly segmentation
of the original kodak image demonstrating some semantic information is
extracted at others layers in the decoder, and so present in the queries.
Thus, to perform a full reconstruction, queries encodes at the same time
semantic and spatial information. More visualizations are given in appendix B

[t]0.49  [t]0.49  [t]0.49  [t]0.49

Figure 10: Kodak original Image Figure 11: First layer of decoder cross-
attention matrix projected on meta queries Figure 12: Decoded Kodak Image with
QPressFormer Figure 13: Tenth layer of decoder cross-attention matrix
projected on meta queries Figure 14: Visualization of information encoded by
output encoder queries and extracted at layer 1 and 11 of the decoder. 14
shows the layer 1 extract more spatial information from the queries whereas 14
shows more semantics extraction with a nearly segmentation of image.

##  5 Conclusion

In this study, we explore learned image compression exclusively through
attention blocks, eliminating convolutional biases. Our model, _QPressFormer_
, introduces ”image queries” for cross-attention-based information
aggregation. Our results show competitive results with convolutional
counterparts with improved perceptual metrics. Future research should further
investigate attention-based image compression. While we focused on analysis
and synthesis components with a basic prior model, incorporating enhancements
on prior model from convolutional models into attention-based architectures
can assess transformer-based compression against state-of-the-art
convolutional counterparts.

## References

  * [1] Johannes Ballé, David Minnen, Saurabh Singh, Sung Jin Hwang, and Nick Johnston,  “Variational image compression with a scale hyperprior,”  in ICLR, 2018.
  * [2] Michael W Marcellin and David S Taubman,  “Jpeg2000: image compression fundamentals, standards, and practice,”  International Series in Engineering and Computer Science, Secs, 2002.
  * [3] Fabrice Ballard,  “Bpg image format,” https://bellard.org/bpg/.
  * [4] Gary J. Sullivan, Jens-Rainer Ohm, Woo-Jin Han, and Thomas Wiegand,  “Overview of the high efficiency video coding (hevc) standard,”  IEEE Transactions on Circuits and Systems for Video Technology, 2012.
  * [5] Google LLC,  “An image format for the web —webp,” https://developers.google.com/speed/webp/.
  * [6] Ashish Vaswani, Noam Shazeer, Niki Parmar, Jakob Uszkoreit, et al.,  “Attention is all you need,”  in NIPS, 2017.
  * [7] Ze Liu, Yutong Lin, Yue Cao, Han Hu, Yixuan Wei, Zheng Zhang, Stephen Lin, and Baining Guo,  “Swin transformer: Hierarchical vision transformer using shifted windows,”  in ICCV, 2021.
  * [8] Binglin Li, Jie Liang, and Jingning Han,  “Variable-rate deep image compression with vision transformers,”  IEEE Access, 2022.
  * [9] Ming Lu, Peiyao Guo, Huiqing Shi, Chuntong Cao, and Zhan Ma,  “Transformer-based image compression,”  in DCC, 2022.
  * [10] Gary Sullivan,  “Versatile video coding (vvc) arrives,”  in 2020 IEEE International Conference on Visual Communications and Image Processing (VCIP), 2020.
  * [11] Fabian Mentzer, Eirikur Agustsson, Michael Tschannen, Radu Timofte, and Luc Van Gool,  “Practical full resolution learned lossless image compression,”  in CVPR, 2019.
  * [12] George Toderici, Sean M. O’Malley, Sung Jin Hwang, Damien Vincent, David Minnen, Shumeet Baluja, Michele Covell, and Rahul Sukthankar,  “Variable rate image compression with recurrent neural networks,”  in ICLR, 2016.
  * [13] George Toderici, Damien Vincent, Nick Johnston, Sung Jin Hwang, et al.,  “Full resolution image compression with recurrent neural networks,”  in CVPR, 2017.
  * [14] Johannes Ballé, Valero Laparra, and Eero P. Simoncelli,  “End-to-end optimized image compression,”  in ICLR, 2017.
  * [15] Johannes Ballé, Valero Laparra, and Eero P. Simoncelli,  “Density modeling of images using a generalized normalization transformation,”  in ICLR, 2016.
  * [16] David Minnen, Johannes Ballé, and George D Toderici,  “Joint autoregressive and hierarchical priors for learned image compression,”  in Advances in Neural Information Processing Systems, 2018.
  * [17] Zhengxue Cheng, Heming Sun, Masaru Takeuchi, and Jiro Katto,  “Learned image compression with discretized gaussian mixture likelihoods and attention modules,”  in CVPR, 2020.
  * [18] Tong Chen, Haojie Liu, Zhan Ma, Qiu Shen, Xun Cao, and Yao Wang,  “End-to-end learnt image compression via non-local attention optimization and improved context modeling,”  IEEE Transactions on Image Processing, 2021.
  * [19] Mu Li, Wangmeng Zuo, Shuhang Gu, Debin Zhao, and David Zhang,  “Learning convolutional networks for content-weighted image compression,”  in CVPR, 2018.
  * [20] Eirikur Agustsson, Michael Tschannen, Fabian Mentzer, Radu Timofte, and Luc Van Gool,  “Generative adversarial networks for extreme learned image compression,”  in ICCV, 2019.
  * [21] Fabian Mentzer, George D Toderici, Michael Tschannen, and Eirikur Agustsson,  “High-fidelity generative image compression,”  in Advances in Neural Information Processing Systems. 2020, Curran Associates, Inc.
  * [22] Eirikur Agustsson, David Minnen, George Toderici, and Fabian Mentzer,  “Multi-realism image compression with a conditional generator,”  in CVPR, 2023.
  * [23] Yochai Blau and Tomer Michaeli,  “The perception-distortion tradeoff,”  in CVPR, 2018.
  * [24] Jacob Devlin, Ming-Wei Chang, Kenton Lee, and Kristina Toutanova,  “BERT: pre-training of deep bidirectional transformers for language understanding,”  in Association for Computational Linguistics, 2019.
  * [25] Hugo Touvron, Thibaut Lavril, Gautier Izacard, Xavier Martinet, Lachaux, et al.,  “Llama: Open and efficient foundation language models,”  arXiv preprint arXiv:2302.13971, 2023.
  * [26] Alexey Dosovitskiy, Lucas Beyer, Alexander Kolesnikov, Dirk Weissenborn, et al.,  “An image is worth 16x16 words: Transformers for image recognition at scale,”  ICLR, 2021.
  * [27] Nicolas Carion, Francisco Massa, Gabriel Synnaeve, Nicolas Usunier, Alexander Kirillov, and Sergey Zagoruyko,  “End-to-end object detection with transformers,”  in ECCV. Springer, 2020.
  * [28] Patrick Esser, Robin Rombach, and Bjorn Ommer,  “Taming transformers for high-resolution image synthesis,”  in CVPR, 2021.
  * [29] Yen-Chun Chen, Linjie Li, Licheng Yu, Ahmed El Kholy, Faisal Ahmed, Zhe Gan, et al.,  “Uniter: Universal image-text representation learning,”  in ECCV. Springer, 2020.
  * [30] Richard Zhang, Phillip Isola, Alexei A. Efros, Eli Shechtman, and Oliver Wang,  “The unreasonable effectiveness of deep features as a perceptual metric,”  in CVPR, June 2018.
  * [31] Jean Bégaint, Fabien Racapé, Simon Feltman, and Akshay Pushparaja,  “Compressai: a pytorch library and evaluation platform for end-to-end compression research,”  arXiv preprint arXiv:2011.03029, 2020.
  * [32] The Eastman Kodak Company,  “The kodak lossless true color image suite,” https://r0k.us/graphics/kodak/.
  * [33] George Toderici, Wenzhe Shi, Radu Timofte, Lucas Theis, Johannes Balle, Eirikur Agustsson, Nick Johnston, and Fabian Mentzer,  “Workshop and challenge on learned image compression (clic2020),” 2020.
  * [34] Eirikur Agustsson and Radu Timofte,  “Ntire 2017 challenge on single image super-resolution: Dataset and study,”  in CVPR Workshops, 2017.
  * [35] Z. Wang, E.P. Simoncelli, and A.C. Bovik,  “Multiscale structural similarity for image quality assessment,”  in The Thrity-Seventh Asilomar Conference on Signals, Systems & Computers, 2003.
  * [36] Martin Heusel, Hubert Ramsauer, Thomas Unterthiner, Bernhard Nessler, and Sepp Hochreiter,  “Gans trained by a two time-scale update rule converge to a local nash equilibrium,”  in Advances in Neural Information Processing Systems, 2017.

Appendix

##  Appendix A Reconstruction example on CLIC and Kodak

An example reconstruction on CLIC is given in figure 15.

An example reconstruction on Kodak is given in figure 16.

Our:
---
bpp=0.26
PSNR=33.07
lpips=0.1230
Original Conv+MSE:
---
bpp=0.29
PSNR=38.46
lpips=0.4682
Our:
---
bpp=0.26
PSNR=33.07
lpips=0.1230
Original Conv+MSE:
---
bpp=0.29
PSNR=38.46
lpips=0.4682
Figure 15: Image reconstruction from CLIC. Zoom comparisons with original and
convolutional equivalent model from [1] are shown on the sides. Even at lower
PSNR, our method improves on the details of the leaves in dark areas. Our:
---
bpp=0.30
PSNR=25.16
lpips=0.2217
Conv+MSE:
---
bpp=0.3461
PSNR=28.37
lpips=0.3161
Original Conv+MS-SSIM:
---
bpp=0.28
PSNR=25.92
lpips=0.3473
Our:
---
bpp=0.30
PSNR=25.16
lpips=0.2217
Conv+MSE:
---
bpp=0.31
PSNR=27.48
lpips=0.4023
Original Conv+MS-SSIM:
---
bpp=0.28
PSNR=25.92
lpips=0.3473
Figure 16: Image reconstruction from Kodak. Zoom comparisons with original and
convolutional equivalent models (optimised with MSE and MS-SSIM) from [1] are
shown on the sides.

##  Appendix B More PCA visualization

[t]0.32  [t]0.32  [t]0.32  [t]0.32  [t]0.32  [t]0.32  [t]0.32  [t]0.32
[t]0.32  [t]0.32  [t]0.32  [t]0.32  [t]0.32  [t]0.32  [t]0.32

Figure 17: Column: Kodak original Image Figure 18: Column: $1^{st}$ layer of
decoder cross-attention matrix projected on meta queries Figure 19: Column:
$10^{th}$ layer of decoder cross-attention matrix projected on meta queries
Figure 20: Visualization of information encoded by output encoder queries and
extracted at layer 1 and 11 of the decoder. Column 20 shows the layer 1
extract more spatial information from the queries whereas column 20 shows more
semantics extraction with a nearly segmentation of image.

Figure 20 shows that the first layer of the decoder extracts spatial
information from the queries. Indeed by visualizing the projection of
attention matrix on a YCbCr space111YCbCr follows by conversion in RGB gives
better visualization colors (then easily switched into RGB space) compute
thanks to 3 first principal directions in the output encoder queries space, we
can see transitions between colors are smooth, with really similar or even
same colors between all 256x256 crops of a same Kodak image (i.e. with the
same meta-queries and so the same 3 directions used for constructing the RGB
space). In an other layer of the decoder, semantics is extracted from queries
which able us to nearly segment the original image.

##  Appendix C Maximum attention map vs mean attention map

We show here the difference between mean attention maps – i.e. attention map
obtained with the average of attention maps across heads and layers – and the
maximum attention maps – i.e. attention map obtained with the maximum of
attention maps across heads and layers. In a case of mean attention maps, as
many areas have nearly 0 attention somewhere in the encoding process, the mean
below all attention values keeping only the ones which are high during the
whole processing. We claim with this visualization we can see spatial encoding
of queries: attention weights of areas where a specific query looks always
during the encoding will not be reduced by the mean compared to other areas.
In the opposite, the maximum attention maps allow us to see all areas where
the query looked at, but information encoded are more difficult to interpret,
as shown on figure 27.

[t]0.49  [t]0.49  [t]0.49  [t]0.49  [t]0.49  [t]0.49

Figure 21: Kodak original Image Figure 22: Maximum attention map of
$\left\\{Q_{E}\right\\}_{36}$ Figure 23: Decoded Kodak Image with QPressFormer
Figure 24: Mean attention map of $\left\\{Q_{E}\right\\}_{36}$ Figure 25:
Image decoded without $\left\\{Q_{E}\right\\}_{36}$ Figure 26: Mean Error of
reconstruction without $\left\\{Q_{E}\right\\}_{36}$ Figure 27: Difference
between mean attention map and maximum attention map. Whereas mean attention
map 27 shows clearly the spatial encoding of the query 36, the maximum
attention map 27 is more difficult to interpret. But we can see that spatial
information is not the only information encoded into queries since different
areas have high attention scores into different 265x256 blocks of the same
image. 27 and 27 corroborate the spatial encoding of queries.


# Estimation of the optimal parameters for K-edge subtraction imaging with
PixiRad-2/PixieIII photon counting detector on a conventional laboratory X-ray
micro-tomograph.

R. Granger  L. Salvo  S. Rolland Du Roscoat  P. Lhuissier

###### Abstract

Photon Counting Detectors (PCDs) open new opportunities in X-ray imaging.
Pixie III is a PCD using simultaneously two energy thresholds. This enables to
acquire images from two distinct energy bins in a single exposure. This is
particularly suited to perform K-Edge Subtraction (KES) imaging with
laboratory sources. In that context, one has however to deal with an energy
bin optimization: narrow energy bins leads to high KES signal, at the expense
of higher noise, while wider energy bins leads to poor KES signal but better
statistics. This work presents a model that aims at finding the optimal
thresholds and source voltage in order to retrieve the best Contrast to Noise
Ratio (CNR) for a given sample. The model also optimizes the configurations
for conventional absorption modality and compares both modalities. We noticed
that the input flux and the energy difference between the thresholds influence
the noise on image. We included this in the model using phenomenological laws.
The model is then compared to empirical optimization by experimental screening
of the parameters on model materials composed of barium, iodine and water.
Finally, a study of the predicted CNRs has function of the sample composition
is presented as an example of usage of the model.

##  1 Introduction

X-ray imaging and computed tomography (CT) are sets of characterization
techniques daily used in a wide variety of fields. Material science, medical
and biological science, industries are examples (see e.g. [1]). The
conventional and most employed technique is based on absorption contrast,
discriminating material by their difference of absorption at a given energy
(monochromatic beam) or by their difference of integrated absorption over an
energy range (with a polychromatic beam). However, it may happen that this
contrast is not large enough to differentiate the materials, a possibility to
enlarge contrast in that case is to exploit phase contrast relative techniques
[2, 3, 4] however these techniques are not easy to handle on laboratory
tomograph.

A variation of the absorption contrast approach is to build contrast from
variation of absorption coefficient measured at two different energies (or
over two distinct energy ranges). The interest of doing so, is to increase the
contrast between two materials when one of them see its K-edge, that
corresponds to a strong variation of absorption coefficient, lying between the
two spectral measurements points. Thus, it is named K-edge Subtraction (KES)
imaging [5]. The classical example of application is the use of contrast agent
as iodine or barium in biological sample, in medical or biological context.
Historically, the technique was developed under somehow monochromatic
condition with laboratory source, using secondary source fluorescence peaks
[6, 7] or additional monochromator [8]. This permits to perform measurements
right before and right after the K-edge. However, these approaches were
restricted to specific chemical elements and limited to low flux. Practically,
this technique was mainly used at synchrotron [9, 10] as the beam offer high
monochromaticity and high flux at the expense of the accessibility of the
technique. However, recent development of new technologies offers mean to
overcome that limitation. Recently, a Compact Light Source has been employed
to perform KES imaging for medical application[11]. While this kind of sources
definitely offers promising perspectives, it is currently not widespread and
correspond to scale of size and cost intermediates between laboratory sources
and synchrotron, but still significantly larger than laboratory sources.
Another strategy to develop KES imaging relies on the recent development of
photon counting detectors, offering detectors with energy-discriminating
capabilities [12]. One of them is the PixieIII photon counting detector [13] :
it has the capacity to count incoming photons over two distinct energy ranges
simultaneously, allowing to acquire two images in a single exposure time when
illuminated with a polychromatic beam. This technology allows taking benefits
of the polychromatic aspects of laboratory X-ray sources, and their associated
advantages.

However, as in that case the absorption coefficient measurements correspond to
mean over energy ranges, and not at specific energies, one has to deal with an
optimization problem. This problem is depicted in figure 1.

Figure 1:  Scheme of the optimization problem of KES imaging with lab source.
Green and purple solid curves are absorbance spectra of the two materials. a)
Taking narrow energy bins: good signal, low statistic, b) large bins: good
statistics, low signal, c) optimal energy bins.

Consider two absorbance spectra of two distinct materials, one with K-edge at
energy $E_{K}$ and another one without K-edge represented by the green and
purple solid curves respectively. The detector thresholds $E_{L}$ and $E_{H}$
and the source voltage U defines the energy range $\Omega_{1}$ and
$\Omega_{2}$ over which the two spectral images are formed (represented by the
blue and red shaded area respectively). In order to get the strongest KES
signal, one may be tempted to use narrow energy ranges on each side of the
K-edge energy as sketched by figure 1.a). However, since laboratory sources
offer limited photon flux, this would result in low statistic images.
Inversely, in order to get lower noise, one may choose large energy bins as in
figure 1.c). In this case, the KES signal is low as the averaging over large
energy bin somehow smooths the K-edge. So the optimal tuning of the energy
ranges over which to perform the K-Edge imaging is a trade-off between high
KES signal and good statistic in the images (fig. 1.b))

To address this problem, He et al. [14, 15] presented a theoretical approach
to predict the optimal width of the energy range. Their model is based on
optimization of contrast to noise ratio of the projection images, the noise
being Poisson noise from the beam, and the contrast signal obtained from Beer-
Lambert law through the sample. Their results assumed an equal bin width of
the two energy bins. Later, Meng et al. [16] used the same approach, the
optimization being based on the Contrast to Noise Ration (CNR) of the
reconstructed images, and showed that bin width are similar when optimizing
CNR on projections or reconstructed images. While this approach answer to the
optimization problem for the general case, detector-specific features like
detector noise or spectral resolution are not taken into account. For
PixieIII, Brun et al. [17] presented results of empirical optimization of one
of the two thresholds on a test phantom made of iodine, water and barium. The
objective of this work is to present an optimization model that includes both
theoretical modelling of general behaviour of detector and empirical detector-
specific characteristics (noise, spectral resolution, counter depth) in order
to refine predictions when using PixieIII detector. By the way, the model aims
at being as less constrained as possible, and optimizes simultaneously: source
voltage and current, both energy thresholds, and predicts if it is better to
use KES imaging or conventional absorption imaging for the considered sample.
This paper is structured as follows. The first section presents the detector.
The second section is dedicated to the model, and is divided in three sub-
model : (i) a source model, (ii) a contrast model and (iii) a noise model. The
third section presents results obtained with the model and compares them with
experimental optimizations. Finally, the fourth section shows an application
example where optimal CNRs for both absorption and KES modality are compared
over a range of chemical composition of the sample.

##  2 Materials

PixiRad2-PixieIII is a photon counting detector with energy-discriminating
capabilities. The detector is made of 2 PixieIII module tiled side-by-side and
made from a 750 µm CdTe crystal bonded to a processing ASIC. This offers a
grid of 1024x402 squared pixels of 62 µm. The detector offers two thresholds –
that will be named low threshold, at energy $E_{L}$ and high threshold at
energy $E_{H}$ – associated with two 15 bit counters, this enabled to acquire
2 images, on two distinct energy bins simultaneously. The so called low-energy
counter counts photons with energy in range $\Omega_{1}=[E_{L},E_{H}]$ while
the high energy counters count photons with energy in range
$\Omega_{2}=[E_{H},\infty[$. The detector has the specificity to offer 3
levels of correction of charge sharing at the hardware level. Throughout this
article, the detector was operated in NPISUM mode, which corresponds to the
finest level of correction of charge sharing, and corrects both the number of
counts, and the energy of the interaction events. ([13, 18, 19]).

For all acquisitions, we used a microfocus reflection source (HAMAMATSU
L12161-07) with tungsten anode and operated with an acceleration voltage of 50
kV.

To develop the absorption model and to evaluate its performance, we used a
phantom sample made with three capillary glass of 1 mm opening diameter,
filled with a solution of 0.3 g/cm3 of KI, barium sulfate a solution of 0.46
g/cm3 of BaSO4 and distilled water for reference (figure 2).

Figure 2: An example of test phantom made with three capillary glass tubes of
1 mm of inner diameter, filled with KI solution, BaSO4 solution and water.

##  3 Model

The objective of the model is to help the operator to choose configuration
parameters based on the sample he wants to image and the objectives of the
acquisition. Given the materials $A$ and $B$ to be contrasted in the sample,
with respective densities $\rho_{a}$ and $\rho_{b}$, the model predicts the
source voltage $U$, current $i$, the two detector energy thresholds ($E_{L}$
and $E_{H}$) and the modality to use (absorption or KES imaging) in order to
achieve the highest Contrast-to-Noise Ratio (CNR) on the obtained radiographs.
Actually, current $i$ and exposure time have equivalent influence on the
variable of interest, namely it acts on the number of counted photon during
exposure. Here the reasoning is done talking with current $i$ only for
simplicity, thinking exposure time as fixed. But the current variable could
then be thought as the product of current times exposure time or exposure
while keeping current fixed without modifying the conclusions. The objective
of the model is represented by the function $f$ in the equation 3.1.

$(U,E_{L},E_{H},i,mode)=f(A,B,\rho_{A},\rho_{B})$ (3.1)

For the context of computed Tomography, [16] showed that optimizing CNR for
radiographs is close to optimizing on reconstructed images. Thus, we focused
in predicting CNR on radiographs, however, as shown by relation 3.1, the
objective is rather to focus on wider optimization by optimizing several
parameters as the choice of these parameters are linked, and particularly, the
two thresholds are optimized.

###  3.1 Analytical expression of CNR

The metric to evaluate CNR for a given modality MOD is defined as the
difference of mean signal between material $A$ and $B$ divided by a standard
deviation representing noise computed on the image with that modality:

$CNR_{MOD}=\frac{|\langle MOD_{A}\rangle-\langle
MOD_{B}\rangle|}{\sigma_{MOD}}$ (3.2)

Where $MOD_{m}$ represent the image obtained with the modality MOD for
material m, $\langle\cdot\rangle$ stand for mean of pixel values over a
representative region of Interest (ROI) and $\sigma_{MOD}$ is the mean of the
standard deviations measured for material A and B.

As it will be developed in that section, the absorption image and KES image
are combinations of raw images. In that section we derive the expression of
absorption contrast $C_{ABS}$ and noise standard deviation $\sigma_{ABS}$ as
well as KES contrast $C_{KES}$ and standard deviation $\sigma_{KES}$ from the
mean and the standard deviation of count levels of acquired images.

For both modality, raw data consists of an image $I_{0}^{\Omega}$ taken
without the sample, and a sample image containing a set $I_{A}^{\Omega}$ of
pixels filled with material $A$ and a set $I_{B}^{\Omega}$ of pixels filled
with material $B$. The input flux on the detector at a given energy $E$ will
be denoted $I_{0}(E)$. Assuming no origins of dispersion other than random
noise, the image without sample is flat with a mean counts given by

$\overline{I_{0}^{\Omega}}=\int_{\Omega}I_{0}(E)dE$ (3.3)

for the image with the sample, the mean count is

$\overline{I_{m}^{\Omega}}=\int_{\Omega}I_{0}(E)e^{-A_{m}(E)}dE$ (3.4)

for pixels of a sample image filled with material $m\in\\{A,B\\}$. Here
$A_{m}(E)$ refer to the absorbance of material $m$ and is given by :

$A_{m}(E)=\frac{4\pi El_{m}}{hc}\Im(n_{m}(E))$ (3.5)

Where $\Im(n_{m}(E))$ is the imaginary part of the refractive index $n_{m}(E)$
at energy $E$ for material $m$ with density $\rho_{m}$ and thickness $l_{m}$
$h$ is the Planck constant and $c$ is the speed of light. The refractive index
$n_{m}(E)$ is taken from xraylib database [20], which gives data with 1 keV
spectral resolution.

As in practice the flat image is an image acquired with high statistics in
order to correct for systematic dispersion of the gray levels, dispersion of
values on flat image will be assumed negligible in comparison with the one of
sample image, which standard deviation will be noted $\sigma_{m}^{\Omega}$,
evaluated in $I_{m}^{\Omega}$.

Let start by establishing the expression of the contrast and noise for
absorption modality. The absorption image is typically obtained by dividing
the sample image with the flat image. The absorption, gathering counts
obtained with both registers of Pixie detector, counting over $\Omega_{1}$ and
$\Omega_{2}$ respectively, will be expressed by

$ABS_{m}=\frac{I_{m}^{\Omega_{1}}+I_{m}^{\Omega_{2}}}{I_{0}^{\Omega_{1}}+I_{0}^{\Omega_{2}}}$
(3.6)

Note that as pixieIII is a photon counting detector, with energy
discriminating thresholds, images without X-Rays are perfectly darks, so there
is no need to correct from any dark levels as it may usually be done with
charge integrating detectors. The mean may be obtained by Taylor expansion
[21] of a bivariate function for the ratio function, to give :

$\langle
ABS_{m}\rangle=\frac{\overline{I_{m}^{\Omega_{1}}}+\overline{I_{m}^{\Omega_{2}}}}{\overline{I_{0}^{\Omega_{1}}}+\overline{I_{0}^{\Omega_{2}}}}$
(3.7)

as $\frac{\operatorname*{Var}(I_{0}^{\Omega_{1}}+I_{0}^{\Omega_{2}})}{\langle
I_{0}^{\Omega_{1}}+I_{0}^{\Omega_{2}}\rangle^{2}}\ll 1$ by assumption, where
$\operatorname*{Var}(X)$ denotes the variance of the statistical variable $X$.
Similarly, assuming
$\frac{\operatorname*{Var}(I_{0}^{\Omega_{1}}+I_{0}^{\Omega_{2}})}{\langle
I_{0}^{\Omega_{1}}+I_{0}^{\Omega_{2}}\rangle}\ll\frac{\operatorname*{Var}(I_{m}^{\Omega_{1}}+I_{m}^{\Omega_{2}})}{\langle
I_{m}^{\Omega_{1}}+I_{m}^{\Omega_{2}}\rangle}$ (low noise in flats compared to
sample images), allows expressing the associated variance as:

$\operatorname*{Var}(ABS_{m})=\frac{(\sigma_{m}^{\Omega_{1}})^{2}+(\sigma_{m}^{\Omega_{2}})^{2}}{(\overline{I_{0}^{\Omega_{1}}}+\overline{I_{0}^{\Omega_{2}}})^{2}}$
(3.8)

Then, the CNR for absorption modality is computed as:

$CNR_{ABS}=\frac{|\langle ABS_{A}\rangle-\langle
ABS_{B}\rangle|}{\sigma_{ABS}}$ (3.9)

where
$\sigma_{ABS}=0.5\times(\sqrt{\operatorname*{Var}(ABS_{A})}+\sqrt{\operatorname*{Var}(ABS_{B})})$
is the mean of standard deviation of absorption image computed between pixels
filled with materials $A$ and $B$ respectively.

KES image may be defined as

$KES_{m}=\log\frac{I_{m}^{\Omega_{2}}}{I_{0}^{\Omega_{2}}}-\log\frac{I_{m}^{\Omega_{1}}}{I_{0}^{\Omega_{1}}}$
(3.10)

Following similar assumption than the ones used to obtain equations 3.7 and
3.8, one obtains:

$\left\langle\frac{I_{m}^{\Omega_{i}}}{I_{0}^{\Omega_{i}}}\right\rangle=\frac{\overline{I_{m}^{\Omega_{i}}}}{\overline{I_{0}^{\Omega_{i}}}}$
(3.11)

and

$\operatorname*{Var}\left(\frac{I_{m}^{\Omega_{i}}}{I_{0}^{\Omega_{i}}}\right)=\frac{(\sigma_{m}^{\Omega_{i}})^{2}}{\overline{I_{0}^{\Omega_{i}}}^{2}}$
(3.12)

where i=1 or 2. Then, using Taylor expansion at first order of the $\log$
function [21], one may estimate:

$\left\langle\log\frac{I_{m}^{\Omega_{i}}}{I_{0}^{\Omega_{i}}}\right\rangle=\log\frac{\overline{I_{m}^{\Omega_{i}}}}{\overline{I_{0}^{\Omega_{i}}}}-\frac{\operatorname*{Var}\left(\frac{I_{m}^{\Omega_{i}}}{I_{0}^{\Omega_{i}}}\right)}{2\cdot\left(\frac{\overline{I_{m}^{\Omega_{i}}}}{\overline{I_{0}^{\Omega_{i}}}}\right)^{2}}$
(3.13)

and

$\operatorname*{Var}\left(\log\frac{I_{m}^{\Omega_{i}}}{I_{0}^{\Omega_{i}}}\right)=\frac{\operatorname*{Var}\left(\frac{I_{m}^{\Omega_{i}}}{I_{0}^{\Omega_{i}}}\right)}{\left(\frac{\overline{I_{m}^{\Omega_{i}}}}{\overline{I_{0}^{\Omega_{i}}}}\right)^{2}}-\frac{\operatorname*{Var}\left(\frac{I_{m}^{\Omega_{i}}}{I_{0}^{\Omega_{i}}}\right)^{2}}{4\cdot\left(\frac{\overline{I_{m}^{\Omega_{i}}}}{\overline{I_{0}^{\Omega_{i}}}}\right)^{4}}$
(3.14)

this lead to

$\langle
KES_{m}\rangle=\left\langle\log\frac{I_{m}^{\Omega_{2}}}{I_{0}^{\Omega_{2}}}\right\rangle-\left\langle\log\frac{I_{m}^{\Omega_{1}}}{I_{0}^{\Omega_{1}}}\right\rangle$
(3.15)

and

$\operatorname*{Var}(KES_{m})=\operatorname*{Var}\left(\log\frac{I_{m}^{\Omega_{2}}}{I_{0}^{\Omega_{2}}}\right)+\operatorname*{Var}\left(\log\frac{I_{m}^{\Omega_{1}}}{I_{0}^{\Omega_{1}}}\right)$
(3.16)

Thus, the CNR for KES modality may be expressed as:

$CNR_{KES}=\frac{|\langle KES_{A}\rangle-\langle
KES_{B}\rangle|}{\sigma_{KES}}$ (3.17)

where
$\sigma_{KES}=0.5\cdot(\sqrt{\operatorname*{Var}(KES_{A})}+\sqrt{\operatorname*{Var}(KES_{B})})$.

Equation 3.9 and 3.17 give expression of CNR computed from input spectrum
$I_{0}(E)$, gray level $\overline{I_{m}^{\Omega_{i}}}$ for both material and
noise $\sigma_{m}^{\Omega_{i}}$ on raw data. The remaining of that section is
dedicated to describe models for these quantities.

###  3.2 Source model

In the above expression, one needs to evaluate $I_{0}(E)$ over the considered
energy ranges. There are basically two possibilities for that: using simulated
spectrum or an empiric one. We used SpekPy([22, 23] to simulate the input
spectrum that illuminates the sample. This model simulates tungsten anode
reflection sources, at desired voltage and current and any source-detector
distances.

Alternatively, we can also obtain an input spectrum by doing a threshold scan
with Pixie-III (i.e., acquiring flat images with increasing threshold value
and then differentiating successive images to obtain a spectrum).

###  3.3 Absorption model

In order to evaluate the terms $\overline{I_{m}^{\Omega_{i}}}$, one needs
evaluations of the absorbance spectra $A_{m}(E)$ for any $E$. As described in
the last section, this is basically done using xraylib. However, due to finite
resolution of the detector, the measured K-edges are smoother than the
tabulated values. In order to take into account that effect, we measure
absorbance spectra with the phantom sample described in section 2. Images with
and without the phantom were acquired with the detector using the low register
only with $E_{L}$ varying from 15 to 50 keV by steps of 1 keV. The mean number
of count per pixel for images without the sample ranges from 20 000 for
$E_{L}$=15keV to 500 for $E_{L}=40$ keV. Then, successive images were
differentiated to obtain the number of counts per energy bins of 1 keV width.
For each energy bins, the image with the sample was divided by the image
without, and the opposite of the logarithm of the result was computed to
obtain spectral image series of the absorbance. For each of the three tube, a
mean spectrum was computed as the mean value over an ROI consisting of a
column of pixel centred in the centre of the tube so that the thickness of
solution correspond to the tubes’ diameter. Finally, the spectrum obtained
from the water tube was subtracted from those obtained from the Iodine and
Barium tubes in order to remove contributions from the glass tube and water.
This results in the absorbance spectra $A_{KI}(E)$ and $A_{BaSO_{4}}(E)$ of
iodine and barium sulfate respectively.

Figure 3: Absorbance spectra for Iodine ans BaSO4 solution measured with
PixiRad-2/PixieIII (solid blue curve) compared against computed ones from
tabulated value of refractive index (orange dashed line), and computed one
with Gaussian kernel of 1.8 keV (red dashed line). a) For iodine and b) for
barium solution

The obtained spectra are the light blue solid curves in figure 3. The
tabulated spectra obtained using xraylib are displayed in orange dashed line.
As one can see, the position of the K-edge measure by the detector is correct,
as well as the global evolution of the absorbance. However, one need to
account for the spectral resolution of the detector which is necessarily
finite and lead to a smoothing of the edges, which is then not measured as a
strict vertical edge. This is done by applying a Gaussian smoothing with a 1.8
keV kernel. This value was chosen in order to reproduce the experimental slope
of the edge. The result is represented by the red dashed line in figure 3. As
one can see, the resulting spectra are really close to the measured one. In
the model, smoothed value of $A_{m}(E)$ (eq. 3.5) are used to compute
$\overline{I_{m}^{\Omega_{i}}}$ terms.

###  3.4 Noise model

The last part of the model is the prediction of the standard deviation
$\sigma_{m}^{\Omega_{i}}$ that one can expect on the acquired images.

As mentioned in the introduction, previous papers ([14, 15, 16]) focused on
modelling Poisson noise from the source. Here we aimed at studying noise
arising from acquisition procedure. Particularly, we focus on the dependence
of the noise to thresholds values. To that purpose, we acquired a set of
images without sample, taking all possible combinations of thresholds, between
25 keV and 46 keV by step of 1 keV and using the low energy register, i.e.
counting photon in range $[E_{L},E_{H}]$. For each combination of threshold,
50 images of 1s of exposure time were taken.

Then we processed the data as follows: for each combination of thresholds, a
flat field corrected image was obtained by dividing the first image of the
series of 50 images by the median of these 50 images, and the resulting image
is then multiplied by the mean value of the median image. This correction
permits to compensate the dispersion of the gray values arising from the
inhomogeneity of the beam, and pixel to pixel systematic variations, letting
noise as the remaining source of dispersion of the gray values, sot that the
image fit the assumptions of the model described in section 3.1

Pixel with a zero value in the median image were labelled as defective. Then,
for each combination of thresholds, mean count level $I$ and standard
deviation $\sigma$ were computed on a ROI of 1850 pixels (370x5) of the flat-
field corrected image, ignoring defective pixels.

Figure 4: Evolution of standard deviation with mean counts. Colour refers to
the mean of the two energy thresholds and size of dots to the difference of
the two energy thresholds.

Figure 4 shows the variation of the relation between mean count level and
standard deviation when varying mean energy $E_{mean}=0.5*(E_{L}+E_{H})$ and
energy bin width $E_{w}=E_{H}-E_{L}$. As one may see, the relation between
noise and input flux is dependent on the values of the thresholds. In order to
investigate that dependence, we modelled the measured noise as being the beam
intrinsic $\sigma_{intr}$ (proportional to the square root of the input flux)
noise affected by a prefactor $\Gamma$ that may depend on the thresholds
values and possibly flux:

$\sigma=\sigma_{intr}\cdot\Gamma(E_{mean},E_{gap},I)$ (3.18)

Representing the intrinsic noise by $\sigma_{intr}=\sqrt{I}$ (thus letting the
proportionality factor in the prefactor) permits to compute the prefactor
with:

$\Gamma=\frac{\sigma}{\sqrt{I}}$ (3.19)

This quantity is represented in figure 5, as a function of $I$, for three
value of $E_{w}$, and colour corresponding to $E_{m}$.

Figure 5: Evolution of the prefactor computed with 3.19 with mean counts.
Colour refers to the mean of the two energy thresholds and size of dots to the
difference of the two energy thresholds.

One can see a first dependence of $\sigma$ with $I$ and/or $E_{mean}:$
$\sigma$ increases with $I$ while $E_{mean}$ decreases. However, as the
experiment is performed without filter on the source, it is expected that the
mean count decreases with the bin mean energy. Additionally, a clear
dependence can be observed on $E_{w}$: points of different size appear to
align along different curves. We suppose that this effect might be due to
finite spectral resolution of the detector: Ideally, when a photon interact
with the detector, it should be counted if and only if its energy lies in the
energy bin defined by threshold values, i.e. the evolution of the probability
for a photon to be counted presents a jump from 0 to 1 at the thresholds
values. However, in practice, that probability evolves smoothly over a finite
energy range centred on the thresholds values. As a result, some photons with
an energy lying in that transition range may be badly counted. This results in
an additional noise which increases with decreasing the bin width, as the
proportion of badly classified photons increases (since the mean difference
between photons energies and threshold value decreases).

Figure 6: a) evolution of prefactor with mean count level for different bin
with $I$ (dots), and fitted model from eq. 3.20. b) and c) evolution of the
fitting parameter amplitude $B$ and characteristic intensity $I_{c}$ of this
model with $E_{w}$ (dots), and fitted parametrization given in eq 3.21 (orange
solid line).

In order to parametrize that dependence, the left panel of figure 6.a) present
the evolution of the prefactor against mean count $I$ for different fixed
$E_{w}$ (dots). From the evolution, to describe the relation between prefactor
$\Gamma$, $E_{w}$ and $I$ we propose a function of the form:

$\Gamma(E_{w},I)=B(E_{w})(1-e^{-\frac{I}{Ic(E_{w})}})$ (3.20)

where $B$ and $I_{c}$ are fitting parameters (amplitude and characteristic
Intensity) that depends on the bin width $E_{w}$. This shape has been fitted
on the data point using least square optimization (solid curves on left panel)
and the evolution of the fitting parameters $B$, and $I_{c}$ against $E_{w}$
are displayed on panels b) and c) of the figure. Evolution for $B$ is
monotonically decreasing. $I_{c}$ is approximately increasing up to a given
value then does not evolve. In fact, from left panel, one see that for largest
bin width, the prefactor depends on I only for really low count rate. Thus, we
suggest parametrizing those evolutions with a decreasing exponential shape for
$B$ and a piecewise constant function for $I_{c}$. Least square fitting of
these shapes gives:

$\begin{split}B(I)&=36.9+18.8e^{-E_{w}/2.28}\\\
I_{c}(I)&=\begin{cases}544&\text{ if }E_{w}\leq 10\text{keV}\\\
147&\text{else}\end{cases}\\\ \end{split}$ (3.21)

for $E_{w}$ in keV.

As a summary, equations 3.18, 3.20 and 3.21 forms a model to estimates the
noise $\sigma$ from the mean count on detector $I$, and the detectors
thresholds, or more precisely, the energy bin width $E_{w}$.

We check the goodness of fit of the parametrization by plotting the estimated
standard deviation from $I$ and $E_{w}$ against the measured one in figure 7.
On this figure, one can see that the point align along the dashed line which
the $y=x$ line. It means that the developed parametrization succeed in
explaining the variability of the dataset.

Figure 7: Validation of the model: estimated standard deviation compared
against measured standard deviation from dataset.

##  4 Results

###  4.1 Resolution

Now we get an expression for estimating contrast and noise from experimental
configuration parameters, we need to specify the procedure to find the optimum
and the associated constrains.

The procedure is detailed by algorithm 1. As a first implementation, the
principle is exhaustive, i.e. it is to evaluate CNR at all points of discrete
search ranges $\Omega_{U}$ (for voltage), $\Omega_{E_{L}}$ (for $E_{L}$) and
$\Omega_{E_{H}}$ (for $E_{H}$). Optionally, during the computation of CNR, the
current is set so that a prescribed count level $CL$ is obtained on the
detector when the sample is out of the beam for each set of voltage and
thresholds.

Algorithm 1 Algorithm to optimize $(U,E_{L},E_{H},\text{mode},i)$ for maximum
$CNR$

for all $(U,E_{L},E_{H})\in\Omega_{U}\times\Omega_{E_{L}}\times\Omega_{E_{H}}$
do

Optionally, compute source spectrum $I_{0}(E)\forall E\in[0,U]$ for a seed
current $i_{s}$

Compute current to reach desired count level $CL$ on registers:
$i(U,E_{L},E_{H})=\frac{CL}{max(I_{0}^{\Omega_{1}},I_{0}^{\Omega_{2}})}\cdot
i_{s}$

Estimate mean gray level on raw images $\overline{I_{m}^{\Omega_{1}}}$,
$\overline{I_{m}^{\Omega_{2}}}$, $\overline{I_{0}^{\Omega_{1}}}$,
$\overline{I_{0}^{\Omega_{2}}}$

Estimate noise on raw images $\sigma_{m}^{\Omega_{1}}$ and
$\sigma_{m}^{\Omega_{2}}$ using equations 3.18 3.21

Compute absorption CNR using eq 3.9

Compute KES CNR using equation 3.17

end for

Find optimum $(U_{opt},E_{L,opt},E_{H,opt})$, mode and current
$i(U_{opt},E_{L,opt},E_{H,opt})$ for highest CNR.

Then, raw images mean counts $\overline{I_{m}^{\Omega_{i}}}$ and
$\overline{I_{0}^{\Omega_{i}}}$, and noise and raw images
$\sigma_{m}^{\Omega_{i}}$ are estimated. From these variables, CNR for both
absorption and KES modality are computed. Once these computations for all
combinations of source voltage, and detector thresholds are performed, the
best configuration is extracted by searching for the maximum CNR obtained.

###  4.2 Validation

In order to evaluate the model, we compared predictions in terms of voltage
and thresholds to empirical optimization performed on a test sample by
evaluating K-Edge CNR for all combination of thresholds in given ranges. The
dataset was obtained using the phantom sample used described in section 2.
This sample allows testing model on three contrast : and I/BaSO4, water/BaSO4,
water/KI. For that acquisition, the source was operated with 50 kV
acceleration voltage, a target current of 14µA, exposing for 1s at an SDD =
281.12 mm, and SOD = 25.08 mm, thresholds varying from 25 to 45 keV by step of
1keV. For each couple of low and high thresholds within that range and for
each energy bin, 50 beam images were acquired and 2 images of the sample were
taken. The images were averaged in order to get one beam image and one sample
image per couple of thresholds from which KES images were computed.

Then, CNR was computed from mean count level and standard deviation measured
over ROIs focusing on each material. The ROIs were made of 3600 pixels
(360x10) centred on their corresponding tubes. While this implies a
heterogeneous thickness of tube traversed by the X-rays, it is necessary to
get sufficient statistics. In parallel, the simulation of that problem was
performed using the model presented in the previous section. As the model is
made from different sub-models, let us first compare some intermediate
variables to evaluate these different sub-models. The top raw of figure 8
compares $\overline{I_{KI}^{\Omega_{1}}}$, the mean gray level on the ROI
focusing on KI tube measured on the low energy register; bottom:
$\sigma_{KI}^{\Omega_{1}}$, the standard deviation measured on that same ROI
and low energy register.

Figure 8: Comparison of intermediates variable as function of the two energy
thresholds. Top row: $\overline{I_{KI}^{\Omega_{1}}}$, bottom row:
$\sigma_{KI}^{\Omega_{1}}$. Left column: simulated quantities, middle column:
experimental quantities, right column: ratio between simulated and
experimental.

We can see that for both variables, simulated values are close to the
experimentally measured ones (last 2 figures of figure 8). The shapes of the
maps of these variables are respected, and quantitatively, the value are
globally well respected, particularly where $E_{w}$ and the number of counts
are highest. For $E_{w}$ getting close to 1, the error increases and the ratio
of simulated values over experimental one is globally about 5-6.

Figure 9: Comparison of simulated (left column) and experimental (right
column) CNR maps for K-edge modality obtained with BaSO4 / KI / water test
phantom: BaSO4/ KI contrast (first row), water/ BaSO4 contrast (second row)
and water/KI (third row)

Then we compare KES CNR maps with respect to the two threshold values (see
figure 9). Optimal configuration were then defined as the maximum of theses
CNR maps. The empirical optimum is represented by the red dot and the
simulated one by the red cross. Additionally, we represent the value of the
corresponding K-edges by grey dashed lines.

Let’s begin the comparison with the top row, i.e. for BaSO4/KI contrast. The
general shape of the modelled map is close to the experimental one. From top
right to bottom left: CNR is low for any configuration with both threshold
above both K-edges, a principal maximum is reached for low thresholds between
both k-edges, and high threshold above Ba K-edge. A secondary maximum at the
bottom left part of the map where both thresholds are below both K-edges and a
narrow “valley” of minima separate both maximum. Keeping in memory that KES
signal measures the variation in absorption between the two energy bins, the
principal maximum corresponds to the Ba K-edge as this implies strong
variation for BaSO4 while moderate ones for KI. The secondary maximum can be
understood by seeing the KES CNR as the difference of absorption CNR between
the two energy bin. This maximum corresponds to a configuration where the high
energy bin includes both K-edges so that the absorption CNR between both
material on that energy bin is low. It remains that the KES CNR in that
configuration is made only from the absorption CNR on the lower energy bin
which focus on the monotonically decreasing part of the spectra. The minimum
between these maxima corresponds to a configuration where the lower energy bin
is centred on the lower K-edge and the higher energy bin on the highest
K-edges. In that configuration, contrast vanishes as for each bin absorption
is stronger for one material on the first part of the bin, and stronger for
the second material in the second part of the bin so that absorption CNR is
low on both energy bin. Quantitatively, the model predicts an optimum at (33
keV, 40 keV) with a CNR of 0.4 while the empirical optimum is (34 keV, 40 keV)
with a CNR of 3.1. The secondary maximum for model and experiment is
respectively about 0.3 and 0.4.

Maps for BaSO4/water or KI/water CNR exhibit simpler shapes as the K-edge
split the map in 3 area: Both maxima are reached with low thresholds clearly
below the K-edge and high threshold at the K-edge value and CNR get low values
and even vanishes if both thresholds are above (upper right part) or below
(bottom left part) the K-edge. Quantitatively, for BaSO4/water contrast, the
model predicts an optimum at (27 keV, 38 keV) with a CNR of 0.6 where it is
(34 keV, 39 keV) with a CNR of 3.0 for the experiment. For KI/water contrast,
the model predicts an optimum at (25 keV, 33 keV) with a CNR of 0.5 where it
is (29 keV, 34 keV) with a CNR of 0.9 for the experiment. The quantitative
differences between model and experiment described in that section will be
commented in section 5.

For a qualitative appreciation of these results, figure 10 shows images
obtained a) with conventional approach, using the images from the high energy
register with $E_{H}=26$ keV (1 phantom images and 50 flat images for flat
field correction) b) KES image obtained by procedure described above at the
predicted optimum for BaSO4/KI contrast, i.e. $E_{L}$=33 keV, $E_{H}$=40 keV,
c) KES image at empirical optimum for BaSO4/KI contrast, i.e. $E_{L}$=34 keV,
$E_{H}$=40 keV, d) KES image at the predicted optimum for BaSO4/water
contrast, i.e. $E_{L}$=27 keV, $E_{H}$=38 keV and e) KES image at the
empirical optimum for BaSO4/water contrast, i.e. $E_{L}$=34 keV, $E_{H}$=39
keV Two observations can be made: i) KES imaging permits to considerably
increase the contrast between KI and BaSO4, even removing water and glass
capillary tubes from the image, and ii) KES images obtained at predicted and
empirical optima are visually very close.

Figure 10: Comparison of the phantom images with a) simple absorption (high
energy bin, $E_{H}=$26 keV, b) KES images, modelled optima or BaSO4/KI
contrast ($E_{L}$=33 keV, $E_{H}$=40 keV) and c) KES images at empirical
optima for BaSO4/KI contrast ($E_{L}$=34 keV, $E_{H}$=40 keV) c) KES images at
modelled optima for BaSO4/water contrast ($E_{L}$=27 keV, $E_{H}$=38 keV) d)
KES images at empirical optima for BaSO4/water contrast ($E_{L}$=34 keV,
$E_{H}$=39 keV). Dark vertical lines are columns of defective pixels, ignored
in the measurements of quantities related to gray level.

###  4.3 Application example

In order to illustrate the application of the model, we show here optimization
results with BaSO4 and KI, fixing density of BaSO4 and varying density of KI.

The materials whose contrast has to be optimized are 1 mm of BaSO4 at
0.46g/cm3, and 1 mm of KI at density in range [0.1g/cm3, 0.5g/cm3]. For each
density of KI, algorithm 1 is applied with $\Omega_{U}$ = [50,80,100],
$\Omega_{E_{L}}=\Omega_{E_{H}}=\llbracket 25\text{ keV},45\text{
keV}\rrbracket$ and $CL=10000$.

Figure 11: Evolution of optimal absorption and KES CNR between BaSO4 and KI as
function of KI density for BaSO4 density fixed at 0.46 g/cm3. Letters
correspond to the different panels of figure 12, plotting absorbance spectra
at difference densities represented by the gray dashed vertical lines.

Figure 11 shows the evolution of the absorption and KES CNR computed from
equations 3.9 and 3.17 respectively. As one may see, the modality to obtain
the highest CNR is not necessarily KES depending on the concentration of the
materials. KES is recommended for density of iodine in range [0.24 g/cm3, 0.37
g/cm3] and absorption is recommended elsewhere. The curve for absorption can
be split in two parts: a first decreasing and a second increasing with iodine
density. The plots a) to d) of figure 12 display the optimal choice of
thresholds for both modality (represented by vertical dashed lines: blue for
absorption and red for KES) in relation with the absorbance spectra of the two
material (solid curves: orange for barium, green for iodine) as well as the
optimal voltage (black solid vertical line) which was found to be 50kV for all
configurations.

Figure 12: Plot of absorbance spectra of both material at different
characteristic densities of iodine, corresponding to vertical gray dashed
lines on top figure: a) 0.1g./cm3, b) 0.3g/cm3 c) 0.4 g/cm3 d) 0.5 g/cm3. Blue
and red vertical dashed lines represent the optimal thresholds for absorption
and KES modality respectively at the corresponding density. The vertical dark
solid line represents the highest photon energy in the source spectrum with
voltage at 50kV, the optimal voltage.

Let us focus on the absorption first: blue solid curve on figure 11, and blue
dashed line on figure 12. For low KI density, BaSO4 has the highest contrast,
and the difference is at its highest above Ba K-edge, the strategy for point
a) and b) is then to focus on that part of the spectrum (E${}_{L}\simeq
40keV$) and choosing the high threshold so both detector count as much as they
can (i.e. $E_{H}$ is set so that CL is reached on both detectors). As the
iodine density increases, the spectral contrast above the Ba K-edge reduces
and gets inverted (iodine more absorbent than barium) between the two K-edges.
This explains why the absorption CNR decreases with iodine density up to 0.3
$g/cm^{3}$. At this point, the best strategy for absorption modality becomes
to focus on the part of the energy domain where KI is the most absorbent.
Then, the thresholds are set so that the low energy bin focus on the part
where iodine is the most absorbent. From here, since the CNR is built from the
fact that KI is the most absorbent, the absorption CNR keeps increasing with
KI density. Now focus on the KES modality: red solid curve on top panel, and
red dashed line on a)-d) panels.

The strategy for KES is to focus the low energy bin on the part between both
K-edges, and the high energy bin above the Ba K-edge. One may notice that the
curve for KES CNR is much flatter than for absorption, this may be an
interesting fact to ensure an acceptable CNR for example sample containing
areas with variable density of KI or unknown composition.

Thus, this example illustrates how the present model permits to take into
consideration different strategies and put them into balance.

##  5 Discussion

First, it has to be precised that the version of the PixiRad2-PixieIII we own
exhibits a salt and pepper noise when used in NPISUM (see figure 10) mode.
This artefact is not a normal behaviour of the detector, but is known to be
linked to errors showing up randomly in the 15bit counters of the detector. A
first consequence of this, is that the current noise level is particularly
high and globally the CNR values described above somewhat underestimate the
performances the technology has to offer. In the future, this extra noise is
brought to disappear with the resolution of this issue, thus the noise model
will need to be updated. However, the current dependencies described by that
model (on the mean count level and $E_{w}$ for the low energy register) are
explained on the basis of other phenomenons intrinsic to the detector
technology. Thus, we may expect that the current shape of the noise model is
not affected by this artefact. Without that noise artefact, one may need to
re-adjust the parameters of the parametrization by reproducing the exposed
experiment, seeing it as a calibration procedure of the current model.

The presented model allows reproducing qualitatively well the variations of
CNR with energy thresholds of the detector. Particularly, the shapes of the
maps on figure 9 are correct, and the predicted optimum is quantitatively
correct for the case where KES imaging is the most relevant. However,
noticeable errors remain. Both flux and noise standard deviation on detector
after going through the materials shows increased error for small energy bins
(i.e. for $E_{H}$ and $E_{L}$ close). As the mean intensity is used by the
noise model to estimate the noise standard deviation, it is expected that the
error seen in $\overline{I_{m}^{\Omega_{i}}}$ propagates in
$\sigma_{m}^{\Omega_{i}}$ and at the end in the computed CNR. The errors seen
on $\overline{I_{m}^{\Omega_{i}}}$ originates from the input flux model
$I_{0}^{\Omega_{i}}(E)$ and sample absorbance model $A_{m}(E)$. Additionally,
the CNR computed in the 4.2 appeared underestimated up to a factor $\simeq 10$
for KI/BaSO4 contrast. This quantitative errors may have several explanations:

  * •

The exposed results demonstrate the performances when used without prior
empirical knowledge of the other elements from the imagery chain, i.e. this
demonstrates the minimal performance of the model. The results may be improved
if one dispose of i) empirical source spectra $I_{0}^{\Omega_{i}}(E)$ ii)
empirical absorbance spectrum of the samples $A_{m}(E)$. Both these spectra
may be measured using Pixie-III by performing a threshold scan.

  * •

Secondly, as it has been shown in section 3.4, the noise depend on both input
flux and thresholds, i.e. it shows a spectral sensibility. However, although
the developed noise model fits well to the data for images without the sample,
the terms $\sigma_{m}^{\Omega_{i}}$ are computed taking into account the
effect of the sample on the input flux only, but not on the spectrum shape.
Thus, it may happen that the parametrization does not give accurate
estimations for those terms, particularly when $E_{w}$ is small. This would be
a point to be studied in finer details.

  * •

The ratio $\frac{\overline{I_{m}^{\Omega_{i}}}}{\sigma_{m}^{\Omega_{i}}}$
ranges from 1 to 5 in the validation example, and is about 1.2 for the
coordinates of the optimum. Additionally, the specific salt and pepper noise
described at the beginning of that section is made of values far from the
mean. For both reasons, the presented Taylor expansions may not hold. For
example, experimentally, in eq. 3.14, the first order term (second term in the
right-hand side) represents 0.17 times the zero order term, which is not
completely negligible yet.

From figure 9, particularly for BaSO4 contrast, one can see that the model
predicts a wider low energy bin than found experimentally (lower $E_{L}$).
Taking into account efficiency of the detector may improve that aspect.
Indeed, at the moment, the input flux model does not take into account
detector specificities : typically efficiency and linearity. It is anticipated
that efficiency decreases significantly for energy above the Cd and Te K-edges
(26.7 keV and 31.8 keV respectively) since part of the incoming photons leads
to fluorescence. The bin width results from a trade-off between KES signal and
statistics. Thus, in practice, energies just below the K-edges are under-
weighted in the measurement of the absorbance on the low energy bin in
comparison to the lowest energies of that bin. The absorbance measured on that
bin is then larger than expected (as spectral absorbance is monotonically
decreasing below the K-edge) and thus the KES signal is reduced. As a
consequence, this pushes to take narrower low energy bin in order to
compensate that loss of KES signal. The optimum is then moved towards higher
$E_{L}$ value. This effect should increase with the K-edge energy of the
sample as the efficiency decreases, this is actually the case when comparing
KI/water contrast to BaSO4 contrast. Additionally, this effect is not
observable for the KI/BaSO4 as in that case the trade-off is to not mix KES
signal from both K-edges as they would annihilate each other. In that case,
both model and experiment results in taking $E_{L}$ just above the lower
K-edge.

Finally, the optimization procedure described by algorithm 1 is quite simple.
One may consider developing it into a refined procedure by adding some
constrains on the product exposure time times source current intensity.
Indeed, when asking a count level CL in the image, the algorithm adapt that
product to reach that level. In practice both parameters are limited. The
source is physically limited by the design of the source and the exposure time
is in practice limited to get a reasonable acquisition time or by the need in
time resolution for in-situ experiment for example. Additionally, one may
include additional modalities in the comparison. For example, it would be
straightforward to consider the absorption image obtained uniquely on the low
energy bin. This would consider an absorption modality where some low and high
energies are ignored. This may be interesting for a case like the one depicted
in figure 12.b where the contrast is mainly due to the absorption difference
in the low energy bin, but the high energy image is contrasted in the opposite
manner and finally reduces the global contrast if summing the counts from both
energy bins.

##  6 Conclusions and perspectives

We designed a model that estimates CNR of both absorption and KES modalities
for Pixie-III photon counting detector. The model optimizes the setup
parameters : voltage U, the product exposure time times intensity, thresholds
$E_{L}$ and $E_{H}$ and acquisition modality (conventional absorption or KES).
The model can be broken down into sub-models that can be used and improved
independently. Particularly, we have presented:

  * •

A noise model, that evaluate the standard deviation of count levels on a
representative ROI of homogeneous material. Besides accounting for intrinsic
Poisson noise which evolves as the square root of the mean intensity, the
noise model includes an empirical parametrization of Pixie-III that reveals
and takes into account a specific behaviour of the detector for small bin
width $E_{w}$. The resulting parametrization reproduces quantitatively the
data.

  * •

A noise propagation model that evaluate CNRs for absorption and KES modalities
on the projection image at the end of the processing chain.

Computationally speaking, the model reduces to simple operations on spectral
integrals which lead to instantaneous computations. The optimization consists
in an exhaustive exploration of the space of acquisition parameters given a
desired count level on the images. The output typically consists of maps of
CNR in the space of acquisition parameters. Confrontation with experimental
exploration of the parameters shows that the model reproduces well the shapes
of the maps. The prediction in terms of thresholds compares within few keV
with experimental. Globally, the model tends to predict lower threshold values
than the ones suggested by experiment. Quantitatively, the error tends to
increase with small energy bins. Additionally, it has to be kept in mind that
the detector currently in our possession suffer from a random noise artefact.
That should not impact the predicted optimal parameters, but may require
recalibration of the noise model for quantitative predictions. Finally, it has
been illustrated how the model enables to compare different strategies to
enhance the contrast and how it allows to understand the underlying trade-offs
the operator has to make.

By these developments, this work constitutes a base from which additional
improvements may be considered. This includes:

  * •

to push forward the noise model by taking into account the effect of sample on
the spectrum of the X-ray beam.

  * •

to account for linearity and efficiency of the detector. Particularly, this
would improve the predictions for sample with K-edge above CdTe K-edges.

  * •

to design more advanced resolution algorithms, this includes : additions of
derivate modalities in the comparison, constrains on the space of set up
parameters, computational improvement of the optimization scheme.

## Supplementary Materials

The model as been implemented in Python and is available on gitlab :
https://gricad-gitlab.univ-grenoble-alpes.fr/TomoX_SIMaP/px3opt

## Acknowledgments

This project was funded by the French National Research Agency
ANR-18-CE42-0005. 3SR is part of LabEx Tec 21 - ANR-11-LABX-0030 and of
Institut Carnot PolyNat (ANR16-CARN-0025). SIMAP is part of LabEx CEMAM
(ANR-10-LABX-44-01).

## References

  * [1] E. Maire and P. J. Withers, _Quantitative X-ray tomography_ , _International Materials Reviews_ 59 (Jan., 2014) 1–43.
  * [2] K. A. Nugent, T. E. Gureyev, D. F. Cookson, D. Paganin and Z. Barnea, _Quantitative phase imaging using hard x rays_ , _Physical review letters_ 77 (1996) 2961.
  * [3] D. Paganin, _Coherent X-ray Optics_.  No. 6. Oxford University Press on Demand, 2006.
  * [4] G. R. Myers, S. C. Mayo, T. E. Gureyev, D. M. Paganin and S. W. Wilkins, _Polychromatic cone-beam phase-contrast tomography_ , _Physical Review A_ 76 (Oct., 2007) 045804.
  * [5] W. Thomlinson, _K-edge subtraction synchrotron X-ray imaging in bio-medical research_ , _Physica Medica_ (2018) 19.
  * [6] B. Jacobson, _Dichromatic Absorption Radiography. Dichromography_ , _Acta Radiologica_ 39 (June, 1953) 437–452.
  * [7] B. K. Rutt, I. A. Cunningham and A. Fenster, _Selective iodine imaging using lanthanum K fluorescence_ , _Medical Physics_ 10 (1983) 801–808.
  * [8] Z. Zhong, D. Chapman, R. Menk, J. Richardson, S. Theophanis and W. Thomlinson, _Monochromatic energy-subtraction radiography using a rotating anode source and a bent Laue monochromator_ , _Physics in Medicine and Biology_ 42 (Sept., 1997) 1751–1762.
  * [9] E. Rubenstein, _Medical imaging with synchrotron radiation_ , _Nuclear Instruments and Methods in Physics Research_ 222 (May, 1984) 302–307.
  * [10] H. Elleaume, A. M. Charvet, S. Corde, F. Est ve and J. F. L. Bas, _Performance of computed tomography for contrast agent concentration measurements with monochromatic x-ray beams: Comparison of K-edge versus temporal subtraction_ , _Physics in Medicine and Biology_ 47 (Sept., 2002) 3369–3385.
  * [11] S. Kulpe, M. Dierolf, E. Braig, B. Günther, K. Achterhold, B. Gleich et al., _K-edge subtraction imaging for coronary angiography with a compact synchrotron X-ray source_ , _PLOS ONE_ 13 (Dec., 2018) e0208446.
  * [12] R. Ballabriga, J. Alozy, M. Campbell, E. Frojdh, E. H. M. Heijne, T. Koenig et al., _Review of hybrid pixel detector readout ASICs for spectroscopic X-ray imaging_ , _Journal of Instrumentation_ 11 (2016) .
  * [13] R. Bellazzini, A. Brez, G. Spandre, M. Minuti, M. Pinchera, P. Delogu et al., _PIXIE III: A very large area photon-counting CMOS pixel ASIC for sharp X-ray spectral imaging_ , _Journal of Instrumentation_ 10 (Jan., 2015) C01032–C01032.
  * [14] P. He, B. Wei, W. Cong and G. Wang, _Optimization of K-edge imaging with spectral CT_ , _Medical Physics_ 39 (Nov., 2012) 6572–6579.
  * [15] P. He, B. Wei, P. Feng, M. Chen and D. Mi, _Material Discrimination Based on K-edge Characteristics_ , _Computational and Mathematical Methods in Medicine_ 2013 (2013) 308520.
  * [16] B. Meng, W. Cong, Y. Xi, B. De Man and G. Wang, _Energy Window Optimization for X-ray K-edge Tomographic Imaging_ , _IEEE transactions on bio-medical engineering_ 63 (Aug., 2016) 1623–1630.
  * [17] F. Brun, V. D. Trapani, J. Albers, P. Sacco, D. Dreossi, L. Brombal et al., _Single-shot K-edge subtraction x-ray discrete computed tomography with a polychromatic source and the Pixie-III detector_ , _Physics in Medicine & Biology_ 65 (Mar., 2020) 055016.
  * [18] V. Di Trapani, L. Brombal, S. Donato, B. Golosio, R. Longo, P. Oliva et al., _36\. Optimization of the acquisition threshold of Photon Counting Detectors (PCDs) used in X-ray medical imaging_ , _Physica Medica_ 56 (Dec., 2018) 84.
  * [19] V. Di Trapani, A. Bravin, F. Brun, D. Dreossi, R. Longo, A. Mittone et al., _Characterization of the acquisition modes implemented in Pixirad-1/Pixie-III X-ray Detector: Effects of charge sharing correction on spectral resolution and image quality_ , _Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment_ 955 (Mar., 2020) 163220.
  * [20] T. Schoonjans, A. Brunetti, B. Golosio, M. S. del Rio, V. A. Solé, C. Ferrero et al., _The xraylib library for X-ray–matter interactions. Recent developments_ , _Spectrochimica Acta Part B: Atomic Spectroscopy_ 66 (2011) 776–784.
  * [21] H. Benaroya, S. M. Han and M. Nagurka, _Probability Models in Engineering and Science_.  CRC Press, June, 2005.
  * [22] G. Poludniowski, A. Omar, R. Bujila and P. Andreo, _Technical Note: SpekPy v2.0—a software toolkit for modeling x-ray tube spectra_ , _Medical Physics_ 48 (2021) 3630–3637.
  * [23] R. Bujila, A. Omar and G. Poludniowski, _A validation of SpekPy: A software toolkit for modelling X-ray tube spectra_ , _Physica Medica_ 75 (July, 2020) 44–54.