# Narrow artificial intelligence with machine learning for real time
estimation of a mobile agent’s location using Hidden Markov Models

Cédric Beaulac and Fabrice Larribe
Département de mathématiques, Université du Québec à Montréal
201, avenue du Président-Kennedy, Montréal, Québec, Canada

###### Abstract

We propose to use a supervised machine learning technique to track the
location of a mobile agent in real time. Hidden Markov Models are used to
build artificial intelligence that estimates the unknown position of a mobile
target moving in a defined environment. This narrow artificial intelligence
performs two distinct tasks. First, it provides real-time estimation of the
mobile agent’s position using the forward algorithm. Second, it uses the Baum-
Welch algorithm as a statistical learning tool to gain knowledge of the mobile
target. Finally, an experimental environment is proposed, namely a video game
that we use to test our artificial intelligence. We present statistical and
graphical results to illustrate the efficiency of our method.

###### Index Terms:

Hidden Markov Model, User tracking, State estimation, Statistical learning,
Non player characters

##  I Motivation

In this paper we address the problem of locating or tracking down a mobile
target; specifically, we are interested in estimating the unknown position of
a mobile agent. To do so, we rely on two main sources of information: our
knowledge of the agent’s usual behaviour and the information obtained from the
environment (e.g. through vision, hearing). We build a proper solution to this
problem by creating a mathematical model, an associated algorithm, and finally
programming this algorithm; the result is autonomous artificial intelligence
(AI).

A common problem in video game environments when multiple games are played in
succession is that a human player will improve with each game played against
the same opponent because he or she learns the opponent’s strategies, while
the typical video game’s AI would not adapt to its opponent. Yannakakis &
Togelius [1] identify core research areas within the field of AI in games.
According to the list of areas they identify, our problem is to both create a
believable agent and to provide it with behavior learning abilities.

We were inspired by Hladky and Bulitko [2] to use a Hidden Markov Model (HMM)
to achieve this goal. In that context, the hidden state is the unknown
position of the mobile agent, and our problem is one of state estimation. Our
main idea is to use the Baum-Welch algorithm[3] as a machine learning tool
that uses the AI’s observations to build up knowledge on a mobile agent over
time. This knowledge is represented by the transition matrix that is used by
our AI to estimate the hidden state at any time. This matrix contains the
various probabilities that the mobile agent moves from one position to
another, i.e. from one state to another. By construction, this matrix is
central in state estimation; thus, the better our AI estimates the transition
matrix, the more accurately it can predict the mobile target’s movements. Our
goal is to demonstrate that this learning algorithm can be efficient in that
precise context. By using this methodology our AI should adapt to and learn
from various situations, which makes it suitable for various competitive video
games.

##  II Related work

Particle models were first used in the context of state estimation problems in
video game environments by Bererton [4]. His article describes the actual
state of AI in video games, defines precisely the major problems and
demonstrates that particle filters are useful in order to solve these
problems. This paper inspired a number of other papers that resort to particle
filters for addressing similar problems. For example, Weber & al [5] used a
particle model with multiple layers in order to represent the state of a
complex game (e.g. an RTS game) and proved the feasibility of the estimation
of this state. Southey & al [6] were among the first to use hidden Markov
models in a video game AI context. They proved that these models could be
useful in tracking the movement of various units in an RTS. Finally, Hladky
and Bulitko [2] decided to compare the efficiency of particle models with
hidden Markov models in order to estimate the position of a moving agent and
noticed that Markov models might produce more precise estimates in certain
situations. A limitation of Hladky’s model was that it had no behavior
learning abilities.

More recently, Stanescu & Certicky[7] used Answer Set Programming in order to
predict an opponent’s production in a RTS video game. Their problem shares
many similarities with ours; for example, they only obtain incomplete
knowledge on the current situation and they are also trying to predict some
elements of the current game state. In order to predict unit production, they
generate every valid unit combination and pick the most probable one as their
estimate. As we will see in section IV-B, the recursive structure of the
forward algorithm will naturally do all of these steps more efficiently in
order to predict the opponent’s location.

Reinforcement learning (RL) is also a very popular approach to build strong
video game AI. Lately, Wang & Tan [8] and Emigh & al. [9] approached the wider
problem of building artificial intelligence that can win a video game base
solely on RL. Even though both of these articles and many others demonstrate
the incredible potential of AI using RL, these models also have their fair
share of weaknesses. When it comes to building AI for a rather complex video
game, both the extremely large state space and high dimensionality will cause
the learning process to be very slow. The uses of RL will also result in AI
that takes actions that maximize the expected reward without really
understanding why. In a certain way, this AI never truly understands the goal
of the game nor the multiple distinct components that can cause a player to
win or lose a game. Wang & Tang [8] tackles briefly these problems by creating
two distinct learning processes: learning the appropriate behavior and
focusing on optimal weapon selection.

Our approach is to further segment AI learning process in order to create a
believable agent and to reduce the potential dimensionality problems. For
instance, for a FPS video game, we would like our AI to understand that it
must aim well, establish a good strategy, estimate its opponent’s location as
accurately as possible and that these subtasks are independent of one another.
Because of the results of [2], we decided to approach the estimation of an
opponent’s location with hidden Markov models. The Markov assumption is very
natural when it comes to modelling movements and, therefore, utilizing a
Markov model would be more appropriate for solving this precise problem as
opposed to other articles that are trying to solve larger problems. To further
improve the actual models, we will use an Expectation-Maximisation algorithm
in order to give the AI the ability to build up it’s own knowledge based on
its experience. This article contains the key results that are explained in
full detail in the associated thesis [10].

##  III Hidden Markov Model techniques and notations

In the following section, we define a Hidden Markov Model, the notations we
use and the algorithms related to that model. This model consists of a Markov
process and an observation function. First, let us define $X_{t}$ as the
Markov process with state space $S=\\{s_{1},s_{2},...,s_{n}\\}$. Because we
will work with discrete-time Markov processes, it is natural to define a
transition function $a_{ij}=P(x_{t}=s_{j}|x_{t-1}=s_{i})$, and an associated
transition matrix A. Finally, the Markov process also consists of an initial
distribution function $\mu_{i}=P(x_{1}=s_{i})$ which is frequently represented
as a vector ${\bf\mu}=\\{\mu_{1},\mu_{2},...,\mu_{n}\\}$.

As a hidden Markov process, the realizations of the Markov process $x_{t}$ are
hidden, i.e. they are not observed. What we obtain instead, the observation,
is a random variable $Y_{t}$ dependent on the hidden state. To complete the
definition of an HMM, we need an observation function:
$b_{i}(y_{t})=P(Y_{t}=y_{t}|x_{t}=s_{i})$.

The challenge is to use a sample of observations $\\{y_{1},y_{2},...,y_{T}\\}$
to make a statistical inference on the hidden Markov process $X_{t}$. For
example, it would be useful to be able to evaluate $P(x_{t}=s_{i}\mid
y_{1},y_{2},...,y_{t})=P(y_{1},y_{2},...,y_{t}\cap
x_{t}=s_{i})/P(y_{1},y_{2},...,y_{t})$. To proceed, we use forward and
backward values as defined by Rabiner [11]. These values will be used in
inference and real-time estimation of the state. We adopt the classic
definition of the forward values, i.e.
$\alpha_{t}(i)=P(y_{1},y_{2},...,y_{t}\cap x_{t}=s_{i})$. The forward values
can be calculated recursively:

$\begin{split}\alpha_{t+1}(j)&=\sum_{i}^{n}\alpha_{t}(i)\ a_{ij}\
b_{j}(x_{t+1}).\\\ \end{split}$ (1)

The initial values are $\alpha_{1}(i)=\mu_{i}b_{i}(x_{1})$ $\forall i$.
Backward values are denoted as usual,
$\beta_{t}(i)=P(y_{t+1},y_{t+2},...,y_{T}|x_{t}=s_{i})$, and they can also be
calculated recursively in a similar manner:

$\begin{split}\beta_{t}(i)=\sum_{j=1}^{n}a_{ij}b_{j}(x_{t+1})\beta_{t+1}(j)\\\
\end{split}$ (2)

using the initial values $\beta_{T}(i)=1$ $\forall i$.

Finally, the Baum-Welch algorithm is needed for parameter estimation. This EM
(Expectation-Maximization) algorithm, adapted for Hidden Markov Model
estimation, is the best known and most widely used algorithm related to Hidden
Markov Model inference. With this tool, every parameter of our model can be
estimated. Note that in our context the initial distribution $\mu$ and the
observation function $b$ will be known and, thus, need not be estimated.

Forward and backward values are calculated as described earlier using the
initial values submitted to the Baum-Welch algorithm. We can use these values
to compute various expected values, two of which are of particular interest.
First:

$\begin{split}\gamma_{t}(i)&=P[x_{t}=S_{i}\mid y_{1},y_{2},...,y_{T},]\\\
&=\alpha_{t}(i)\ \beta_{t}(i)/L_{T},\\\ \end{split}$ (3)

where $L_{T}$ is the likelihood of the observations. Using the forward and
backward values we can also compute the following expected value:

$\begin{split}\xi_{t}(i,j)&=P[x_{t-1}=s_{i},x_{t}=s_{j}\mid
y_{1},y_{2},...,y_{T},]\\\
&=\alpha_{t-1}(i)a_{i,j}b_{j}(y_{t})\beta_{t}(j)/L_{T}.\\\ \end{split}$ (4)

Finally, after we’ve computed $\gamma_{t}(i)$ and $\xi_{t}(i,j)$ for all $t,i$
and $j$, we can perform the maximization step. As mentioned, only the
transition function has to be estimated. We have to find the transition matrix
that maximizes the likelihood of our observations. We thus obtain an estimator
very similar to the classic Markov model maximum-likelihood estimator:

$\begin{split}\hat{a}_{i,j}&=\frac{\sum_{t=2}^{T}\xi_{t}(i,j)}{\sum_{t=2}^{T}\gamma_{t}(i)}=\frac{\sum_{t=2}^{T}\xi_{t}(i,j)}{\sum_{t=2}^{T}\sum_{j}^{n}\xi_{t}(i,j)}\\\
&=\frac{\sum_{t=2}^{T}P[x_{t-1}=s_{i},x_{t}=s_{j}\mid
y_{1},y_{2},...,y_{T}]}{\sum_{t=2}^{T}P[x_{t}=s_{i}\mid
y_{1},y_{2},...,y_{T}]}.\end{split}$ (5)

To conclude this section, we present the details of how the algorithm works.
The first step is to submit the initial values of the Hidden Markov Model’s
parameters to the algorithm. These consist of the initial distribution
function, the transition function and the observation function. The second
step is to calculate forward and backward values using the parameters we
currently have. The third step involves computing the expected values as we
did in equation (3) and (4) using the $\alpha$’s and $\beta$’s we computed in
the previous step. Finally, the maximization step consists of using both of
these expected values to estimate the parameters of the Hidden Markov Models.
We then repeat this process starting from the second step until the
parameter’s estimates are stable.

To summarize, here are the steps of the algorithm.

1) Submit the initial values of the parameters A, $\mu$ and $b$.

2) Compute $\alpha_{t}(i)$ and $\beta_{t}(i)$, $\forall t$ and $\forall i$
using the forward and backward algorithms.

3) Compute the expected values $\gamma_{t}(i)$ and $\xi_{t}(i,j)$, $\forall t$
and $\forall i,j$ using the forward and backward values.

4) Estimate the parameters A, $\mu$ and $b$ using the expected values as in
(5).

5) Return to step 2) and repeat until a desired level of convergence is
attained.

##  IV Main idea of our approach

###  IV-A Mathematical modelling

First, we define what the Markov process represents. Recall that our main goal
is to estimate the unknown position of a mobile agent in a restricted area; in
our experimental environment, which is a competitive video game, the mobile
agent we are tracking is our opponent. Because its position is unknown, we
suppose that its movements follow a Markov chain, i.e. the location where the
agent will be at time $t+1$ depends only on its location at time $t$. More
formally, we define our Markov process, $X_{t}$, as the position of the mobile
agent, at time $t$. Because this position is unknown we never observe the
realizations of this process, but we obtain observations that are dependent on
these hidden realizations. These observations are the realizations of a Hidden
Markov Model.

Next, we define the numerous parameters we need to work with an HMM. First, we
need to define the state space of the Markov process. To this end we must
introduce our defined environment, since we are trying to estimate the unknown
position of a mobile agent in a known restricted area. In the context of a
video game, we call this the map, which is the equivalent of a sports field;
it is where the action takes places. For example, in the context of smart
vehicle technologies the map could be a certain city. Our goal is to create a
methodology that efficiently estimates a mobile agent’s location on this map.
To begin, we grid this map, creating possible positions, and the set of all
these positions forms our finite state space $S=\\{s_{1},s_{2},...,s_{n}\\}$.

Now we describe how we define our initial distribution, $\mu$. In our video
games application, we will assume that the initial distribution is directly
implemented in the game, and known by all players. As most competitive video
games have fixed starting locations for each player, a comparison can be made
to a game of chess where each pawn starts at a precise location. Some video
games have a set of possible starting locations, one of which is chosen
randomly. Generally, these possible locations are near one another and are
known by every experienced player. Thus, the initial distribution is
considered to be known by our AI in our experimental set up.

As mentioned earlier, Hidden Markov Models involve an unobservable realization
of the underlying Markov process $x_{t}$, and an observation which is a random
variable $Y_{t}$ dependent on the hidden realization. In other words, the
observations represent the stream of information received from the
environment. In our context, the observations are the set of positions we can
eliminate at each unit of time and, thus, are random and dependent on the
hidden state (the actual position of the mobile agent). The problem consists
of using these observations, via a well-constructed observation function, in a
manner that will help us track down the agent. In our experimental
environment, our methodology is used by an AI battling against a human player.
In that context, the observation is the set of positions observed by the AI’s
avatar inside the video game at each time frame.

Figure 1: A representation of a vision radius, common in video games. The dark
brown square represents our avatar, the red square the mobile agent, and the
green tiles the positions we observe.

Figure 1 depicts two situations faced when tracking down a mobile agent. We
either directly observe the agent or not. If we see the mobile target, as in
the right image in Figure 1, the estimate of the current hidden state is quite
simple: it is not an estimate, it is agent’s actual position. The model gets
interesting when we actually have to estimate our opponent’s position, as in
the left image in Figure 1.

Suppose we are in a situation where estimation is required, and let us define
the observation function $b$. To proceed, we introduce $W$, which is the set
of all visible positions, namely the green tiles in Figure 1. In other
contexts, these observations can be obtained with human eyes, camera,
satellite or even other senses such as hearing. At this precise time, these
are the locations where we know the agent is not found. They form a set of
states we can eliminate as the possible hidden state of the Markov process.
What we observe is that the agent is not at these positions. We develop the
following observation function:

$b_{i}(y_{t})=P(Y_{t}=y_{t}\mid X_{t}=s_{i})=1-\bf{1}[s_{i}\in W].$ (6)

Note that we are in a situation where we observe a set of positions $W$ but we
do not observe the mobile agent we are tracking. This simple but logical
definition of $b$ implies that the probability of obtaining our observation,
which is not seeing the agent in any elements of $W$, conditional on the fact
that he is there is 0. It also implies that the probability of not seeing the
agent in the set of locations we observed knowing he is not in this set is 1.
This binary observation function not only cancels the set of visible locations
but also eliminates every path that would have passed by these positions at
this particular time, due to the recursive nature of the forward algorithm.
We’ll explain this in detail in the next subsection. Note that if the hidden
state is $s_{i}$ at time $t$ and that at this very time, $s_{i}\in W$, we
simply define $b_{i}(y_{t})=1$ and $b_{j}(y_{t})=0$ $\forall j\neq i$.

Notice that this observation function could be modified easily in order to
include uncertainties about certain position by defining a function that take
values somewhere between 0 and 1 at some point. This could prove useful if we
were interested in using other types of observations such as the sounds
produced by the mobile agent, which are useful observations in multiple video
games and other applications.

To build the desired AI we need to use a Hidden Markov Model for two distinct
tasks, namely: real-time estimation of the mobile agent’s position, and
machine learning for improved future estimation. It should be noted that
Hladky and Bultiko [2] and Southey & al. [6] used a hidden semi-Markov Model
instead of a Hidden Markov Model in their research. Even though we considered
using this model, we found that the gain in realism was not worth the
programming cost because it directly affected the efficiency of the learning
component of our AI. The Baum-Welch algorithm for a Hidden semi-Markov Model
was indeed slower to compute. Using such an algorithm in the context of a
video game might create enormous loading time between games, and the gain of
realism might not help the AI that much. For the benefit of the machine
learning part of the AI we are designing, we decided to use a standard Hidden
Markov Model.

###  IV-B Real-time state estimation

Recall that if we observe the mobile agent, estimation is useless. Therefore,
let us suppose we are in a situation where estimation is required, i.e. we de
not observe the mobile agent. We evaluate the probability that the Markov
process is in state $s_{i}$ at time $t$ considering current and past
observations for all states at all discrete moments. This probability is
denoted $P_{t}(s_{i})$ and is calculated using the forward algorithm:

$\begin{split}P_{t}(s_{i})&=P(x_{t}=s_{i}\mid
Y_{1}=y_{1},Y_{2}=y_{2},...,Y_{t}=y_{t})\\\ &=\frac{P(x_{t}=s_{i}\cap
Y_{1}=y_{1},Y_{2}=y_{2},...,Y_{t}=y_{t})}{P(Y_{1}=y_{1},Y_{2}=y_{2},...,Y_{t}=y_{t})}\\\
&=\frac{\alpha_{t}(i)}{\sum_{i=1}^{n}\alpha_{t}(i)}.\\\ \end{split}$ (7)

Let us suppose that we observe the position $s_{i}$ at time $t$. Because of
the way the observation function is set up, $b_{i}(y_{t})=0$ and thus
$P_{t}(s_{i})=0$, which is desired and which contributes to eliminating a set
of possible positions at time $t$. Specifically, this methodology eliminates
every state $s_{i}\in W$ if the agent is not observed.

That being said, the important part is the fact that knowing the hidden state
was not $s_{i}$ at time $t$ implies not only that $P_{t}(s_{i})=0$ but also
that $\alpha_{t}(i)=0$. Because the $\alpha_{t+1}$’s are computed recursively
using the $\alpha_{t}$’s, it simplifies further calculations by eliminating
elements of the summations. In fact, it eliminates every state at time $t+1$
that requires the agent to be at state $s_{i}$ at time $t$. Additionally, from
our tracking point of view, it also eliminates every path that passes through
any positions $s_{i}\in W$ for every turn. This is the main driver for this
methodology’s efficiency.

###  IV-C Machine learning

Recall that a common problem in video game environments when multiple games
happen consecutively on the same map is that a human player will improve with
each game played against the same opponent because he or she learns the
opponent’s strategies, while the typical video game’s AI would not adapt to
its opponent. We refer to this ability of human players as adaptive memory. In
any context where multiple experiments take place in the same environment, the
machine learning process we are using could aid the artificial intelligence’s
learning from each experiment. Given that observations only eliminate slowly
some possible positions throughout the experiment, real-time estimation relies
heavily on the transition matrix A. It represents all the information about
the mobile agent the AI possesses. Our goal is to update this information
after each set of observations is obtained. We choose to use the Baum-Welch
[3] algorithm using multiple independent sets of observations to create the
transition matrix that represents the agent’s behaviour. We built a slightly
modified version from the one presented in equation 5. To build that equation,
we assume each set of observations to be independent and we were inspired by
the maximum likelihood estimator of a regular Markov model. At each iteration
the transition matrix is re-estimated using this formula:

${\hat{a}}_{i,j}=\frac{\sum^{m}\sum_{t=2}^{T}\xi_{t}(i,j)}{\sum^{m}\sum_{t=2}^{T}\sum_{j}^{n}\xi_{t}(i,j)},$
(8)

where $m$ is the number of experiments realized.

Because we want our AI to be flexible and to detect changes in strategies, we
will not use every game in our machine learning phase. Let us say the AI has
played 100 games against a certain player. The Baum-Welch algorithm will use
all of these 100 observation sequences to estimate the transition matrix A.
Therefore, if the human player changes his strategy drastically, it will take
multiple games until this new strategy has enough weight in the estimation to
change the estimate of the transition matrix. For this reason, the matrix used
by the AI for real-time state estimation will be a mix between a general long-
term knowledge of the agent and short-term knowledge represented by a
transition matrix built upon only the last few games. This matrix will be more
detailed in the next section.

##  V Experimental set up and results

After developing the algorithm and implementing it in $C$++ we were interested
in testing this new AI with adaptive memory. For this we had to design our own
experimental environment, a video game, which was also implemented in $C$++.
We designed a turn-based game in which the two players play one at a time and
make one move per turn. The human player can see up to two tiles away, while
the AI only see the positions directly surrounding it, as seen in Figure 1.

In the game we created for testing purposes, the goal of the human player is
to get to a particular position and stay there for a known number of turns. In
contrast, the goal of the AI-controlled character is to prevent that from
happening by touching its opponent. We chose this specific objective because
it was directly affected by the quality of the estimations and because it did
not require that we code a complicated strategy surrounding the AI’s movement.
We invented a pirate concept to intrigue the players who ran tests on our
game. In this concept, the player is a pirate trying to steal a treasure from
the captain on a remote island.

One game mechanic we added was parrots, to go with our pirate concept. They
are scattered on various tiles on the map, and the human player knows where
they are. They serve as “cameras” for the captain, as they constantly give the
AI observations to inform it whether the mobile agent is on one of these
tiles. They not only help the AI by giving it more observations, but they also
motivate the player to build complex strategies that take the cameras into
account.

Figure 2: An image of the video game we designed as an experimental
environment. The player controls the humanoid on top, while the AI controls
the humanoid that starts at the bottom. Goals are marked by the red Xs on the
map.

Figure 2 represents the map we designed. In this figure, players are currently
at their respective initial positions.

The AI obtains observations from its environment as defined by equation 6. At
each turn, a vector of observations is created. All elements of this vector
are 1’s except the elements representing observed positions, which are 0’s.
The set of observed positions is the set of tiles in the vision radius of the
AI, as illustrated in Figure 1, and the tiles observed by the cameras,
corresponding to the parrots in our concept. These observations are used by
the AI along with the transition matrix, which represents information on the
opponent’s behavior, and the initial distribution to estimate the unknown
position of its opponent as described in equation 7.

For this experiment, we also had to code a path-finding algorithm. Even though
our primary goal was to estimate the position of a mobile agent, we needed the
AI to move autonomously in order to obtain data and conduct tests. Using
equation 7, the AI can find the most probable hidden state and define it as
its estimate of its opponent’s position at each turn. This position is
generally defined as the goal in the path-finding literature. Considering each
state as a node in a graph, we implemented the Dijkstra’s algorithm [12] to
let the AI determine the shortest path to that estimate and start moving
towards it. At each turn, the AI determines a new goal, calculates the
distance between the goal and each of its surrounding positions, and chooses a
position to move from among these surrounding positions.

To prevent convergence problems, we used various logarithmic transformations,
that were inspired by Rabiner [11] and that are well explained by both Durbin
et al. [13] and Beaulac [10].

Finally, note that the forward and backward algorithms were coded using a
matrix notation as presented by Zucchini [14] to simplify computations. Given
that these values are a sum of a product, matrix multiplication can be useful.
The forward vector $\alpha_{t}$, can be obtained by calculating the following:
$\alpha_{t}=\mu\mbox{{B}}(y_{1})\mbox{{A}}\mbox{{B}}(y_{2})\mbox{{A}}\cdot\cdot\cdot\mbox{{A}}\mbox{{B}}(y_{t})$
and the backward values:
$\beta_{t}=\mbox{{A}}\mbox{{B}}(y_{t+1})\mbox{{A}}\mbox{{B}}(y_{t+2})\cdot\cdot\cdot\mbox{{A}}\mbox{{B}}(y_{T})$
where :

$\textbf{B($y_{t}$)}=\begin{bmatrix}b_{1}(y_{t})&0&0&\cdots&0\\\
0&b_{2}(y_{t})&0&\cdots&0\\\ \vdots&\vdots&\vdots&\ddots&\vdots\\\
0&0&0&\cdots&b_{n}(y_{t})\end{bmatrix}.$

We tested our AI under different scenarios to determine whether or not the AI
could learn from its experiences to estimate the position of the mobile agent
more accurately. To this end, we used the mean distance through out a game
between the real position and the estimate as our statistic. We wanted to
compare the performance of our AI with adaptive memory with an AI similar to
Hladky’s [2], which is not capable of machine learning, using the same
statistic. We couldn’t compare the precision of our estimates with estimates
produced by a more recent AI using reinforcement learning since those AI’s
aren’t actually trying to predict their opponent’s location. We played
multiple games using different strategies against various AI’s to collect
data. For testing purposes we created multiple variations of the experiment.

First, we developed multiple strategies for the player. By a strategy, we mean
a particular path from the starting location to one of the treasures that
initially leads to victory. For example, a strategy could be to start by
touching the top left parrots in order to gets the pirate to believe the
player is going that way, then directly head to the right treasure by using
the top lane.

The first test set up was to use the exact same path for ten games in row.
We’ve run this experiment with two different winning strategy. Note that
intelligence with no learning abilities would lose those ten games. We did
this to evaluate the capacity of our AI with machine learning to adapt, to
estimate the position of its opponent more precisely and finally to win. The
second experiment consists of using a path for eight games then changing
drastically to another winning strategy. We needed 4 winning precise path to
run this experiment twice. This test was intended to verify the ability of our
intelligence to recognize a change in strategy and adapt to the new tactic.
Finally, we tried alternating between two winning strategies. We could thus
examine how efficiently our AI could learn two paths simultaneously. We have
produced various graphics to demonstrate the evolution of our statistic
throughout the experiment.

First, we present our results when we used the exact same tactic ten times in
a row. We repeated this process with two different winning strategies.

Figure 3: These graphics represent the 12th, 16th, 24th and 30th turn of two
games. The left column is the first game against a new player, while the
graphics in the right column are from the tenth game against the same player.
The blue square represents the human player, the intelligence’s avatar is the
dark brown square, the green triangles represent the parrots and objectives
are represented by the golden squares. The player used the exact same path in
all ten games. Figure 4: Average distance for each game. Each figure
represents a precise path. The green line represents an AI with no machine
learning and which used a uniformly distributed transition matrix. The black
line represents the same AI but this time one that used a transition matrix
built using game experience, as Hladky did. Finally, the red line was obtained
by playing against our AI with adaptive memory.

We represent the real-time estimations and learning process elegantly via heat
maps of probabilities. These maps are a representation of our experimental
video game map in which we indicate the estimates performed by the artificial
intelligence for each position at each turn. Rather than writing the exact
probabilities we represented them using colors: black tiles are walls and do
not count as possible positions; dark red tiles represent a high probability;
white tiles represent the fact that the AI’s estimate for this position is 0.
The player is identified as the blue square, and the intelligence’s avatar is
the dark brown square. The green triangles represent the parrots. Finally,
objectives are represented by the golden squares.

The Figure 3 contains heat maps from two different games at the exact same
turns. The left column represents the first game, and the right column
represents the tenth, after the learning process. Not only can the graphics in
Figure 3 help us comprehend how this AI works, that is how it sees the map,
but they can also help us visualize the learning process. We can see the
evolution of the state estimation obtained from the AI with adaptive memory.
In the left column, the player is frequently in the pale red tile, meaning the
AI does not consider this position to be highly probable. In contrast, the
right column clearly shows that the player’s real position is often painted in
dark red. This illustrates that after multiple games, the AI learned the path
used by its opponent and now considers the agent’s real position highly
probable. We can also notice the difference in the path used by the
intelligence to pursue its opponent. The AI’s avatar is much closer to the
agent in the right column throughout the game.

We compare the performance of our AI with machine learning to that of the AI
with no learning component using the statistics discussed earlier. Figure 4
represents the evolution of the mean distance between the real position of the
mobile agent and the estimation of the position from three different AI’s. We
observe that for both of these strategies, the AI with no learning component
performs in exactly the same way in every single game. However, the mean
distance between the mobile agent’s position and the estimate from the AI we
designed decreases constantly. This indicates that our AI’s estimates are more
precise as it learns from past games.

We finally tested that the statistic of interest was significantly smaller for
the AI we designed compared to Hladky’s AI. Our AI’s estimates were more
precise in both of these situations ( t = -7.3049 with p-value = $6.225\times
10^{-09}$ and t-value = -3.2673 with p-value = 0.00619).

Figure 5: Average distance for each game. Each figure represents a different
run of the experiment. The green line represents an AI with no machine
learning and which used a uniformly distributed transition matrix. The black
line represents the same AI but this time one that used a transition matrix
built using game experience, as Hladky did. Finally, the red line was obtained
by playing against our AI with adaptive memory.

It was important to verify how our AI performed in other scenarios. Because a
human player can adapt to multiple strategies, we had to be sure that our
intelligence would adapt to an abrupt change in tactics. That is why we played
eight games using the same strategy, then completely changed our strategy to a
new one for the eight following games. Figure 5 contains the evolution of the
test statistics we designed. It illustrates that the AI with adaptive memory
produces superior estimates after the learning process. Even though the change
in strategy surprises the AI at the ninth confrontation causing poor
performance, we can clearly see that our AI adapts quickly to this new tactic.

Note that the AI we developed had more precise estimates after the learning
process. The difference was statistically significant for the eighth game
(t-value = -6.3822 with p-value = $9.97\times 10^{-08}$ and t-value = -7.2901
with p-value = $4.362\times 10^{-08}$) and for the 16th game (t-value =
-5.0932 with p-value = $2.241\times 10^{-05}$ and t-value = -8.0063 with
p-value = $1.387\times 10^{-09}$). Finally, we did not want the learning
process to become a handicap. It was important for us that to ensure that
after learning a certain strategy, our AI would not produce poor estimates for
every other strategy. That is why we also compared our test statistic right
after the change of strategy. We can confirm that even after learning a
tactic, our AI’s estimates are as precise as an AI with no machine learning
(t-value = -1.1967 with p-value = 0.2386 and t-value = -0.2683 with p-value =
0.7894).

Our last experiment was set up to check how our AI would perform against an
opponent that alternated between two initially winning strategies. We wanted
to verify our AI’s ability to learn more than one strategy simultaneously. To
demonstrate the learning process we present another set of heat maps, shown in
Figure 6.

Once again, we can notice a substantial difference in the shade of red of the
player’s tile, as shown in Figure 6. We can also easily observe a difference
in the path that the AI’s avatar used. This time it was also much closer to
the human player in the 11th game than in the first game.

Figure 6: These graphics represent the 6th, 10th, 12th and 14th turn of two
games. The left column is the first game against a new player, while the right
column represents the eleventh game against the same player. The blue square
represents the human player, the intelligence’s avatar is the dark brown
square, the green triangles represent the parrots and objectives are
represented by the golden squares. The player alternated between two
strategies. Figure 7: Average distance for each game. Each figure represents a
different run of the experiment. The green line represents an AI with no
machine learning and which used a uniformly distributed transition matrix. The
black line represents the same AI but this time one that used a transition
matrix built using game experience, as Hladky did. Finally, the red line was
obtained by playing against our AI with adaptive memory.

Lastly, we graphically illustrate the test statistic for this experiment.
Figure 7 indicated that our AI slowly improved its estimates. The process is a
bit more unstable and slower than the process with only one strategy, but it
is still highly effective. We can see that our AI is able to learn at least
two strategies simultaneously.

We now verify that the intelligence with adaptive memory produced more precise
estimations after the learning process took place. Since it is constantly
switching between two strategies, we examined the 11th and the 12th matches.
The mean distance between the estimate and the real position of our bot was
smaller both for the 11th game ( t-value = -5.3526 with p-value = $6.451\times
10^{-06}$ and t-value = -6.8703 with p-value = $1.983\times 10^{-08}$) and the
12th game ( t-value = -7.8525 with p-value = $1.783\times 10^{-09}$ and
t-value = -6.2715 with p-value = $2.553\times 10^{-07}$).

##  VI Conclusion

We have developed a way to represent the problem of tracking a mobile agent
mathematically, specifically, through the use of a Hidden Markov Model. After
finding a way to efficiently program the HMM’s algorithms, we built artificial
intelligence that could use this model to accurately estimate the unknown
position of a mobile agent in a restricted area.

We programmed our own video game for testing purposes, and constructed test
statistics and multiple graphical tools to visualize the effect of our work.
We found that AI demonstrated very strong results; the real-time estimation
procedure combined with the learning process made it difficult for human
players to win after a number of rounds. Furthermore, the processing time
involved in the real-time estimation component was efficient to the point of
not being noticeable, which affords it the potential to be implemented in
numerous applications. Finally, the adaptive memory also produced very strong
results. Although it requires more time to adjust the various settings
precisely, our AI can be a powerful tool in many applications, especially
video games applications.

Much work remains before statistics and probability models can be used
efficiently in the construction of artificial intelligence, but we believe
that our work described above is evidence that this kind of interdisciplinary
research has great potential. We do hope that statistics become part of the
culture of development of new artificial intelligence, and that our work has
contributed to that trend.

## Acknowledgments

The authors personally thank Louis-Alexandre Vallières-Lavoie and Guillaume
Racicot for their impressive work on the program. The authors also gratefully
acknowledge the financial support from the NSERC and the FRQNT. Finally, the
authors would like to thank Sorana Froda for her important contribution and
Ali Al-Aradi for his most helpful suggestions.

## References

  * [1] G. N. Yannakakis and J. Togelius, “A panorama of artificial and computational intelligence in games,” _IEEE Transactions on Computational Intelligence and AI in Games_ , vol. 7, no. 4, pp. 317–335, Dec 2015.
  * [2] S. Hladky and V. Bulitko, “An evaluation of models for predicting opponent positions in first-person shooter video games,” _Computational Intelligence and Games, 2008. CIG ’08. IEEE Symposium On_ , pp. 39–46, 2008.
  * [3] L. E. Baum, T. Petrie, G. Soules, and N. Weiss, “A maximization technique occurring in the statistical analysis of probabilistic functions of markov chains,” _The Annals of Mathematical Statistics_ , vol. 41(1), pp. 164–171, 1970.
  * [4] C. Bererton, “State estimation for game ai using particle filters,” _Proceedings of the AAAI 2004 Workshop on Challenges in Game AI_ , pp. 36–40, 2005.
  * [5] B. Weber, M. Mateas, and A. Jhala, “A particle model for state estimation in real-time strategy games,” _Proceedings of the Seventh AAAI Conference on Artificial Intelligence and Interactive Digital Entertainment_ , pp. 103–108, 2011.
  * [6] F. Southey, W. Loh, and D. Wilkinson, “Inferring complex agent motions from partial trajectory observations,” _Proceedings of the Eighth International Joint Conference on Artificial Intelligence_ , pp. 2631–2637, 2007\.
  * [7] M. Stanescu and M. Čertický, “Predicting opponent’s production in real-time strategy games with answer set programming,” _IEEE Transactions on Computational Intelligence and AI in Games_ , vol. 8, no. 1, pp. 89–94, March 2016\.
  * [8] D. Wang and A. H. Tan, “Creating autonomous adaptive agents in a real-time first-person shooter computer game,” _IEEE Transactions on Computational Intelligence and AI in Games_ , vol. 7, no. 2, pp. 123–138, June 2015.
  * [9] M. S. Emigh, E. G. Kriminger, A. J. Brockmeier, J. C. Príncipe, and P. M. Pardalos, “Reinforcement learning in video games using nearest neighbor interpolation and metric learning,” _IEEE Transactions on Computational Intelligence and AI in Games_ , vol. 8, no. 1, pp. 56–66, March 2016.
  * [10] C. Beaulac, “Intelligence artificielle avec apprentissage automatique pour l’estimation de la position d’un agent mobile en utilisant les modèles de markov cachés,” Master’s thesis, Université du Québec à Montréal, 201, avenue du Président-Kennedy, 2015.
  * [11] L. R. Rabiner, “Numerically stable hidden markov model implementation,” _Proceedings of the IEEE_ , vol. 77(2), pp. 257–286, 1989.
  * [12] E. W. Dijkstra, “A note on two problems in connexion with graphs,” _Numerische Mathematik_ , vol. 1, pp. 269–271, 1959.
  * [13] R. Durbin, S. Edyy, A. Krogh, and G. Mitchison, _Biological Sequence Analysis : Probabilistic Models of Proteins and Nucleic Acids_. Cambridge University Press, 1998.
  * [14] W. Zucchini and I. MacDonald, _Hidden Markov Models for Time Series: An Introduction Using R_ , ser. Chapman & Hall/CRC Monographs on Statistics & Applied Probability. CRC Press, 2009. [Online]. Available: http://books.google.ca/books?id=LDDzvCsdVs8C


# Inner product regularized multi-energy X-ray tomography for material
decomposition

Salla-Maaria Latva-Äijö  Department of Mathematics and Statistics, University
of Helsinki, Helsinki, Finland  Filippo Zanetti  School of Mathematics,
University of Edinburgh, Scotland, UK  Ari-Pekka Honkanen  Department of
Physics, University of Helsinki, Helsinki, Finland  Simo Huotari  Department
of Physics, University of Helsinki, Helsinki, Finland  Jacek Gondzio  School
of Mathematics, University of Edinburgh, Scotland, UK  Matti Lassas
Department of Mathematics and Statistics, University of Helsinki, Helsinki,
Finland  Samuli Siltanen  Department of Mathematics and Statistics, University
of Helsinki, Helsinki, Finland

(June 2023)

## Abstract

Multi-energy X-ray tomography is studied for decomposing three materials using
three X-ray energies and a classical energy-integrating detector. A novel
regularization term comprises inner products between the material distribution
functions, penalizing any overlap of different materials. The method is tested
on real data measured of a phantom embedded with Na2SeO3, Na2SeO4, and
elemental selenium. It is found that the two-dimensional distributions of
selenium in different oxidation states can be mapped and distinguished from
each other with the new algorithm. The results have applications in material
science, chemistry, biology and medicine.

###### Contents

  1. 1 Introduction
  2. 2 Measurement model
  3. 3 Preconditioned interior point method
  4. 4 Experimental data
    1. 4.1 The phantom
    2. 4.2 X-ray measurement
    3. 4.3 K-edges in X-ray attenuation
  5. 5 Results
    1. 5.1 Error measures of the simulations
    2. 5.2 Effect of the new regularization term
    3. 5.3 Results with simulated data
    4. 5.4 Results with experimental X-ray data
  6. 6 Conclusions

##  1 Introduction

We study multi-energy X-ray tomography for material decomposition,
generalizing the reconstruction method we introduced in [5]. There we
demonstrated with simulated data the feasibility of one inner product
regularizer term in the context of dual-energy X-ray tomography for a two-
material object. Here we work with three different X-ray energies and
materials, applying our generalized reconstruction method to measured
monochromatic X-ray data. We note that our methods assume classical energy-
integrating sensors and do not need photon counting detectors.

The new regularizer takes the form

$\langle\mathbf{g}^{(1)},\mathbf{g}^{(2)}\rangle+\langle\mathbf{g}^{(1)},\mathbf{g}^{(3)}\rangle+\langle\mathbf{g}^{(2)},\mathbf{g}^{(3)}\rangle,$
(1)

where the non-negative vector $\mathbf{g}^{(j)}$ models the distribution of
material $j$ inside the target. Now the penalty (1) promotes the a priori
information that the target consists of just one material at any point, since
(1) becomes zero if at least one function in each inner product vanishes at a
pixel.

Clinical computed tomography (CT) devices are usually built with conventional
X-ray tubes, which produce polychromatic broad-bandwidth beams. Such beams
lead to beam-hardening artifacts and insensitivity to the chemical composition
of the imaged object. These limitations can be overcome with monochromatic
X-ray beams. They also allow taking advantage of sharp changes in attenuation
as function of energy (K-edge variations) of certain materials.

Monochromatic data was produced with a new method, where a polychromatic X-ray
beam is directed through a spherically bent crystal analyser [8]. It
monochromatises and refocuses the beam at the Rowland circle which acts as a
secondary source. This makes it possible to adjust the energy of the
diffracted photons and record the spatially resolved changes in the
attenuation coefficients. All of this can be done easily in a traditional
X-ray laboratory (as opposed to a synchrotron facility). See [8].

Previous literature on the topic is mostly based on the discrete tomography
approach [1, 2, 3, 6, 7]. Our algorithm is different from them as it is based
on the regularizing penalty based on inner products. We remark that using an
energy-sensitive photon counting detector opens up further possibilities, see
e.g. [9, 15, 12]. We do not discuss such avenues of investigation further in
this paper.

Our reconstruction results suggest that the inner product regularizer is
promising for material decomposition. There are potential applications in
chemistry, biology and medicine.

##  2 Measurement model

We restrict here to the intersection of a physical body with a two-dimensional
square $\Omega\subset\mathbb{R}^{2}$. The measured X-rays are thus assumed to
travel in the plane determined by $\Omega$. This restriction is only for the
simplicity of the exposition and computation; the method generalizes also to
higher dimensions.

We discretize $\Omega$ into $N\times N$ square-shaped pixels. There are three
unknowns: non-negative $N\times N$ matrices $G^{(1)}$, $G^{(2)}$ and $G^{(3)}$
modelling the distributions of material 1, material 2 and material 3,
respectively. The number $G^{(\ell)}_{i,j}\geq 0$ represents the concentration
of material $\ell$ in pixel $(i,j)$, where $i$ is the row index and $j$ is the
column index. In numerical computations, we represent the elements of the
material matrices $(G^{(1)},G^{(2)},G^{(3)})\in(\mathbb{R}^{N\times N})^{3}$,
as a vertical vector

$\mathbf{g}=\left[\\!\\!\begin{array}[]{l}\mathbf{g}^{(1)}\\\
\mathbf{g}^{(2)}\\\
\mathbf{g}^{(3)}\end{array}\\!\\!\right]\in\mathbb{R}^{3N^{2}},$

where the vertical vector $\mathbf{g}^{(j)}\in\mathbb{R}^{N^{2}}$ is formed by
stacking the columns of the matrix $G^{(j)}\in\mathbb{R}^{N\times N}$.

We record X-ray transmission data with three different energies, calling them
low, middle and high energy. Measurements are resulting in three
$M$-dimensional data vectors called $\mathbf{m}^{L}$, $\mathbf{m}^{M}$ and
$\mathbf{m}^{H}$. For example, the low-energy measurement is given by

$\mathbf{m}^{L}=c_{11}A^{L}\mathbf{g}^{(1)}+c_{12}A^{L}\mathbf{g}^{(2)}+c_{13}A^{L}\mathbf{g}^{(2)},$
(2)

and similar form applies to all the different energies.

The materials attenuate the X-rays with individual strengths described by the
energy dependent attenuation constants $c_{11}>0$, $c_{12}>0$, $c_{13}>0$, …
and so on. The empirical values of $c_{ij}$ can be determined by measuring
pure samples of each of three known materials. The $D{\times}N^{2}$ matrix
$A^{L}$ encodes the geometry of the tomographic measurement using the standard
pencil-beam model. Here $D$ is the number of sensors (pixels) in the linear
X-ray detector. The geometric system matrices $A^{M}$, $A^{H}$ and $A^{L}$ can
be different from each other. We have kept the different letters for all the
system matrices in the theory part, but in the actual calculations, we used in
the same measurement geometry and system matrix for all the measurements:
$A^{H}=A^{M}=A^{L}$.

We can combine all the measurements in a unified linear system:

$\mathbf{m}=\left[\\!\\!\begin{array}[]{l}\mathbf{m}^{L}\\\ \mathbf{m}^{M}\\\
\mathbf{m}^{H}\end{array}\\!\\!\right]=\begin{bmatrix}c_{11}A^{L}&c_{12}A^{L}&c_{13}A^{L}\\\
c_{21}A^{M}&c_{22}A^{M}&c_{23}A^{M}\\\ c_{31}A^{H}&c_{32}A^{H}&c_{33}A^{H}\\\
\end{bmatrix}\left[\\!\\!\begin{array}[]{l}\mathbf{g}^{(1)}\\\
\mathbf{g}^{(2)}\\\
\mathbf{g}^{(3)}\end{array}\\!\\!\right]=\mathcal{A}\mathbf{g}.$ (3)

The core idea in this triple-energy X-ray tomography for material
decomposition is to choose the three energies so that three materials respond
to them differently. The solution of (3) is rather analogous to solving a
system of three linear equations for three variables.

We use our novel variational regularization approach in the space
$\mathbf{g}\in\mathbb{R}^{3N^{2}}$, including a non-negativity constraint:

$\widetilde{\mathbf{g}}_{\alpha,\beta}=\mathop{\mbox{arg\,min}}_{\mathbf{g}^{(j)}\geq
0}\left\\{\|\mathbf{m}-\mathcal{A}\mathbf{g}\|_{2}^{2}+\alpha\mathcal{R}(\mathbf{g})+\beta\mathcal{S}(\mathbf{g})\right\\},$
(4)

where $\alpha,\beta>0$ are regularization parameters, $\mathbf{g}^{(j)}\geq 0$
means that the elements of the vector are non-negative numbers and the
regularizer $\mathcal{R}$ can be any of the standard choices such as the
Tikhonov penalty

$\mathcal{R}(\mathbf{g})=\|\mathbf{g}\|_{2}^{2}.$ (5)

The novelty arises from the term that penalises all Inner Products (IP)
between vectors $\mathbf{g}^{(1)}$, $\mathbf{g}^{(2)}$ and
$\mathbf{g}^{(3)}\in\mathbb{R}^{N^{2}}$. The regularization term
$\mathcal{S}(g)$ is defined as double the sum of those:

$\begin{split}\mathcal{S}(\mathbf{g})=\mathcal{S}(\left[\\!\\!\begin{array}[]{l}\mathbf{g}^{(1)}\\\
\mathbf{g}^{(2)}\\\
\mathbf{g}^{(3)}\end{array}\\!\\!\right]):=2(\langle\mathbf{g}^{(1)},\mathbf{g}^{(2)}\rangle+\langle\mathbf{g}^{(1)},\mathbf{g}^{(3)}\rangle+\langle\mathbf{g}^{(2)},\mathbf{g}^{(3)}\rangle)\\\
=2\sum_{i=1}^{N^{2}}\mathbf{g}^{(1)}_{i}\mathbf{g}^{(2)}_{i}+2\sum_{i=1}^{N^{2}}\mathbf{g}^{(1)}_{i}\mathbf{g}^{(3)}_{i}+2\sum_{i=1}^{N^{2}}\mathbf{g}^{(2)}_{i}\mathbf{g}^{(3)}_{i}.\end{split}$
(6)

Together with the non-negativity constraint, $\mathcal{S}$ promotes the point-
wise separation between three materials: at each pixel, at least two of the
images, $G^{(1)}$, $G^{(2)}$ or $G^{(3)}$, needs to have a zero value to make
$\mathcal{S}$ minimal.

The quadratic program resulting from the application of the novel variational
regularization is solved using an Interior Point Method (IPM) [14, 4]; we
develop an efficient preconditioner for the normal equations which guarantees
the spectrum of the preconditioned matrix to remain independent of the IPM
iteration.

##  3 Preconditioned interior point method

In this section we discuss optimization with preconditioned interior point
method (IPM).

By combining the use of Tikhonov regularizer (5) and the Inner Product
regularizer (6), which promotes the point-wise separation of three materials,
we arrive at the constrained quadratic programming task

$\mathop{\mbox{arg\,min}}_{\mathbf{g}^{(j)}\geq
0}\left\\{\|\mathbf{m}-\mathcal{A}\mathbf{g}\|_{2}^{2}+\alpha\|\mathbf{g}\|_{2}^{2}+\beta\,\mathbf{g}^{T}L\mathbf{g}\right\\},$
(7)

where

$L=\begin{bmatrix}0&I&I\\\ I&0&I\\\ I&I&0\\\ \end{bmatrix},$

with nine blocks of size $N^{2}{\times}N^{2}$ each. In general, using the
regularizer with $k$ materials is equivalent to adding a regularization term
$\mathbf{g}^{T}L\mathbf{g}$, with

$L_{k}=E_{k}\otimes I_{N^{2}}-I_{kN^{2}}$

where $E_{k}$ is the matrix of all ones of size $k$,
$E_{k}=[\mathbf{e}\,\mathbf{e}\,\dots\mathbf{e}]$

The problem may be written as an explicit quadratic program with inequality
(non-negativity) constraints

$\mathop{\mbox{arg\,min}}_{\mathbf{g}^{(j)}\geq
0}-\mathbf{m}^{T}\mathcal{A}\mathbf{g}+\frac{1}{2}\mathbf{g}^{T}(Q_{1}+Q_{2})\mathbf{g}$
(8)

where

$\displaystyle Q_{1}$
$\displaystyle=\begin{bmatrix}c_{11}^{2}&c_{11}c_{12}&c_{11}c_{13}\\\
c_{11}c_{12}&c^{2}_{12}&c_{11}c_{13}\\\
c_{11}c_{13}&c_{12}c_{13}&c^{2}_{13}\end{bmatrix}\otimes(A^{L})^{T}A^{L}+$ (9)
$\displaystyle+\begin{bmatrix}c_{21}^{2}&c_{21}c_{22}&c_{21}c_{23}\\\
c_{21}c_{22}&c^{2}_{22}&c_{22}c_{23}\\\
c_{21}c_{23}&c_{22}c_{23}&c^{2}_{23}\end{bmatrix}\otimes(A^{M})^{T}A^{M}+$
$\displaystyle+\begin{bmatrix}c_{31}^{2}&c_{31}c_{32}&c_{31}c_{33}\\\
c_{31}c_{32}&c^{2}_{32}&c_{32}c_{33}\\\
c_{31}c_{33}&c_{32}c_{33}&c^{2}_{33}\end{bmatrix}\otimes(A^{H})^{T}A^{H},$
$Q_{2}=\begin{bmatrix}\alpha&\beta&\beta\\\ \beta&\alpha&\beta\\\
\beta&\beta&\alpha\end{bmatrix}\otimes I_{N^{2}},$ (10)

where $\otimes$ represents the Kronecker product.

For an arbitrary number of materials $k$, $Q_{1}$ is still equal to
$\mathcal{A}^{T}\mathcal{A}$ and $Q_{2}$ is equal to

$Q_{2}=\beta L_{k}+\alpha
I_{kN^{2}}=\underbrace{\begin{bmatrix}\alpha&\beta&\beta&\dots&\beta\\\
\beta&\alpha&\beta&\dots&\beta\\\ \beta&\beta&\alpha&\ddots&\vdots\\\
\vdots&\vdots&\ddots&\ddots&\beta\\\
\beta&\beta&\dots&\beta&\alpha\end{bmatrix}}_{Q^{\prime}_{2}\in\mathbb{R}^{k\times
k}}\otimes I_{N^{2}}$ (11)

Recall this important property about the Kronecker product:

######  Lemma 3.1.

Given two square matrices $T$ and $Z$, the eigenvalues of the Kronecker
product $T\otimes Z$ are given by $t\cdot z$, where $t$ is an eigenvalue of
$T$ and $z$ is an eigenvalue of $Z$.

To study the convexity of problem (8), let us first analyze the eigenvalues of
matrix $Q_{2}^{\prime}$.

######  Lemma 3.2.

When considering $k$ materials, matrix $Q_{2}^{\prime}$ has eigenvalues
$\alpha-\beta$ with multiplicity $k-1$ and $\alpha+(k-1)\beta$ with
multiplicity $1$.

###### Proof.

The result follows from noticing that

$Q_{2}^{\prime}-(\alpha-\beta)I_{k}=\begin{bmatrix}\beta&\dots&\beta\\\
\vdots&\ddots&\vdots\\\ \beta&\dots&\beta\end{bmatrix}$

has rank $1$ and

$Q_{2}^{\prime}-(\alpha+(k-1)\beta)I_{k}=\begin{bmatrix}-(k-1)\beta&\beta&\dots&\beta\\\
\beta&-(k-1)\beta&\beta&\vdots\\\ \vdots&&\ddots&\\\
\beta&\dots&\beta&-(k-1)\beta\end{bmatrix}$

has rank $k-1$.

∎

We can now give a sufficient condition for the convexity of problem (8).

######  Lemma 3.3.

If $\alpha\geq\beta$, problem (8) is convex, for any number of materials $k$.

###### Proof.

If $\alpha\geq\beta$, then matrix $Q_{2}^{\prime}$ has only positive
eigenvalues, due to Lemma 3.2. Therefore, matrix $Q_{2}$ is positive semi-
definite, due to Lemma 3.1. Since $\mathcal{A}^{T}\mathcal{A}$ is always
positive semi-definite, $Q=\mathcal{A}^{T}\mathcal{A}+Q_{2}$ is positive semi-
definite if $Q_{2}$ is. Therefore, the problem is a convex quadratic program.

∎

Therefore, in the following we assume that $\alpha\geq\beta$.

We now generalize the eigenvalue bounds for the preconditioned matrix found in
[5] to the case of $k$ materials. We recall that the matrix of the Newton
system that arises at each interior point iterations can be written as
$M=Q_{1}+Q_{2}+\Theta$, where $Q_{1}$ is given in (9), $Q_{2}$ is given in
(10) and $\Theta$ is a diagonal matrix with strictly positive entries. We
restrict to the case of imaging using the same angles for all energies (i.e.
$A^{L}=A^{M}=A^{H}=:A$); then, in the case of $k$ materials, $Q_{1}$ can be
written as

$Q_{1}=(\mathcal{C}^{T}\mathcal{C})\otimes(A^{T}A),\qquad\mathcal{C}=\begin{bmatrix}c_{11}&c_{12}&\dots&c_{1k}\\\
c_{21}&c_{22}&\dots&c_{2k}\\\ \vdots&&\ddots&\\\
c_{k1}&c_{k2}&\dots&c_{kk}\end{bmatrix}.$

We suppose to precondition matrix $M$ with
$P=(\mathcal{C}^{T}\mathcal{C})\otimes(\rho I_{N^{2}})+Q_{2}+\Theta$, where
$\rho\in\mathbb{R}$ is an approximation of the mean diagonal element of
$A^{T}A$. Notice that this matrix has a $k\times k$ block structure with
diagonal blocks and is thus easy to invert and apply to a vector.

Then, the following eigenvalue bounds hold:

######  Lemma 3.4.

Given $M=(\mathcal{C}^{T}\mathcal{C})\otimes(A^{T}A)+Q2^{\prime}\otimes
I_{N^{2}}+\Theta$ and $P=(\mathcal{C}^{T}\mathcal{C})\otimes(\rho
I_{N^{2}})+Q2^{\prime}\otimes I_{N^{2}}+\Theta$, the eigenvalues of $P^{-1}M$
are real and lie in the interval

$\Big{[}\frac{\alpha-\beta}{\rho\sigma_{\text{max}}^{2}(\mathcal{C})+\alpha+(k-1)\beta},\frac{\sigma_{\text{max}}^{2}(A)\sigma_{\text{max}}^{2}(\mathcal{C})+\alpha+(k-1)\beta}{\rho\sigma_{\text{min}}^{2}(\mathcal{C})+\alpha-\beta}\Big{]}.$

###### Proof.

For a certain eigenvector $\mathbf{v}$ of unit norm, the eigenvalues of
$P^{-1}M$ can be written as

$\lambda=\frac{\mathbf{v}^{T}M\mathbf{v}}{\mathbf{v}^{T}P\mathbf{v}}=\frac{\mathbf{v}^{T}\big{(}(\mathcal{C}^{T}\mathcal{C})\otimes(A^{T}A)\big{)}\mathbf{v}+\mathbf{v}^{T}\big{(}Q_{2}^{\prime}\otimes
I_{N^{2}}\big{)}\mathbf{v}+\mathbf{v}^{T}\Theta\mathbf{v}}{\mathbf{v}^{T}\big{(}(\mathcal{C}^{T}\mathcal{C})\otimes(\rho
I_{N^{2}})\big{)}\mathbf{v}+\mathbf{v}^{T}\big{(}Q_{2}^{\prime}\otimes
I_{N^{2}}\big{)}\mathbf{v}+\mathbf{v}^{T}\Theta\mathbf{v}}.$

Using basic properties of the Kronecker product (see [5] for more details),
notice that

$\mathbf{v}^{T}\big{(}(\mathcal{C}^{T}\mathcal{C})\otimes(A^{T}A)\big{)}\mathbf{v}\,\in\,\big{[}0,\sigma_{\text{max}}^{2}(A)\sigma_{\text{max}}^{2}(\mathcal{C})\big{]},$
$\mathbf{v}^{T}\big{(}(\mathcal{C}^{T}\mathcal{C})\otimes(\rho
I_{N^{2}})\big{)}\mathbf{v}\,\in\,\big{[}\rho\sigma_{\text{min}}^{2}(\mathcal{C}),\rho\sigma_{\text{max}}^{2}(\mathcal{C})\big{]},$
$\mathbf{v}^{T}\big{(}Q_{2}^{\prime}\otimes
I_{N^{2}}\big{)}\mathbf{v}\,\in\,\big{[}\alpha-\beta,\alpha+(k-1)\beta\big{]}.$

With similar arguments as in [5] the thesis can be easily derived. ∎

Similar bounds can be found in the case $A^{L}\neq A^{M}\neq A^{H}$. Notice
that when $k$ becomes larger, the effect of the inner product regularization
coefficient $\beta$ on the eigenvalues of $P^{-1}M$ becomes more significant
with respect to $\alpha$. In particular, if $k$ grows, for the same values of
$\alpha$ and $\beta$, the spectrum of $P^{-1}M$ widens, worsening the efficacy
of the preconditioner. A possible solution to this issue is to increase the
gap between $\alpha$ and $\beta$ when the number of materials grows.

##  4 Experimental data

We use monochromatic X-ray data described in [8]. There, three-dimensional
distributions of selenium in different oxidation states is mapped and
distinguished from background and each other with absorption edge contrast
tomography.

###  4.1 The phantom

In [8], three holes were drilled into a cuboid acrylic glass (PMMA) phantom.
The holes where filled with selenium compounds: elemental Se, Na2SeO3 and
Na2SeO4 and starch. See Figure 1 for a schematic drawing of the phantom, where
we show the materials with different colors. Figure 2 shows a photograph of
the sample. Note how there seems to be no visible difference between the
measured materials: they all have same color and structure.

Figure 1: The figure shows a schematic drawing of the phantom. The cuboid PMMA
(grey) was drilled with three holes each of which was filled with a mixture of
a selenium compound [elemental Se (red), Na2SeO3 (green) and Na2SeO4 (blue)]
and starch. The holes were capped with tissue paper (very light grey). Figure
2: Here is a photograph of the tomography setup. The detector in the picture
is a different TimePIX-based model than the one used in this work. There is no
visible difference between the materials.

###  4.2 X-ray measurement

The phantom described in Section 4.1 was measured with a XAS-CT setup, see
Figure 3. This kind of arrangement monochromatises the polychromatic X-rays
from the X-ray tube when they travel through a spherically bent crystal
analyser.

Figure 3: The schematic drawing of the XAS-CT setup. The polychromatic X-rays
produced by the X-ray tube are monochromatised with the spherically bent
crystal analyser. The sample to be imaged is illuminated by the
monochromatised beam by moving it away from the Rowland circle so that the
defocused beam completely covers it. The beam transmitted through the sample
is recorded with a position-sensitive detector.[8]

Figure 4: Two random example projection images of the sample. These projection
images were taken with energies 12.658 keV and 12.662 keV. Figure 5: Example
sinogram of the data. The row 80 (81 in Matlab syntax) was picked from all the
projection images to the sinograms. There was all in all three different
sinograms, corresponding the energies 12.658 keV, 12.662 keV and 12.685 keV.
These energies correspond the low, middle and high energies in our measurement
model. This sinogram is produced from 12.662 keV projection images.

###  4.3 K-edges in X-ray attenuation

Our IP method works on its best when we investigate materials with X-ray
energies which are near the K-edges of the materials. In X-ray absorption
spectroscopy, the K-edge is a sudden increase in X-ray absorption, occurring
when the energy of the X-rays is just above the binding energy of the atoms’
electron shell. The term is based on the International Union of Pure and
Applied Chemistry (IUPAC) X-ray notation, where the atomic orbitals are
labelled with letters K, L, M, N… and so on. The innermost electron shell is
known as the K-shell.

Physically, the sudden increase in attenuation is caused by the photoelectric
absorption of photons. For this interaction to occur, the photons must have
more energy than the binding energy of the K-shell electrons (K-edge). A
photon having an energy just above the binding energy of the electron is more
likely to be absorbed than a photon having an energy just below the binding
energy or significantly above it.

Figure 6 illustrates this. In the figure, we have the different attenuation
values $\frac{\mu}{\rho}$ as a function of photon energies in keV:s. The plots
show how the attenuation varies with three different materials, which are pure
selenium Se (blue), Na2SeO3 (red) and PMMA (green). We can see the sudden
increase in attenuation of the selenium containing samples at the measured
X-ray energy range (black box). We can compare with the PMMAs attenuation
values, which are very smooth at that energy range. Clear k-edges at the
measured energy range make the differentiation of materials much easier for
the used material decomposition. The attenuation values in figure 6 are based
on the attenuation value tables on the NIST (National Institute of Standards
and Technology) X-ray attenuation online simulator [11].

Figure 6: The mass attenuation coefficient, $\frac{\mu}{\rho}$, as a function
of photon energy, for elemental Selenium Se, Na2O3 and PMMA. In the plot, the
K-edges of the Se and Na2O3 are clearly visible, which makes the decomposing
of the materials easier with the IP method. In the energy window that we use
(12,54-12,80 keV), which is marked with a black box in the figure, the
attenuation values of PMMA are smooth, which means that it appears similar in
every energy images. So we can subtract it out from all the different energy
sinograms, by subtracting low energy sinogram from all the others.

##  5 Results

In this section we show the results with the simulations and finally with real
data. The section also includes a brief demonstration of the effect of the new
regularization term. We have measured the quality of the simulations with some
error measures, which will be explained below.

###  5.1 Error measures of the simulations

We were avoiding the inverse crime in our simulations by adding some noise and
modelling error into our data. The noise level in the simulations was selected
to be one percent, and we added modelling error by rotating the phantom 45
degrees and interpolating the difference.

We calculated quality measures for the simulations by comparing the original
phantoms with the resulting reconstructions. The classical errors were
calculated correspondingly than in [5]. We defined the $L_{2}$-error

$L_{2}\text{-error}=\frac{\text{norm}(\text{phantom}(:)-\text{reconstruction}(:))}{\text{norm(phantom(:))}}.$
(12)

In addition to that we calculated the structural similarity index (SSIM) [13]
and Haar wavelet-based perceptual similarity index (HaarPSI) [10]. We have
written down the mean values of the errors for all the materials, by summing
up the numerical errors of the different materials and dividing by three. The
same holds for the mean pixel error, which we have then listed as a percentage
of correctly classified pixels (Corr.pix %) in table 2.

###  5.2 Effect of the new regularization term

In this section we show how regularization effects to our reconstruction
results with IP method, when using simulations. We use digital phantoms,
because the point is just to demonstrate the effect of the inner product
regularization parameter.

The coefficients which we used for the simulated materials with different
energies are listed in the table 1.

Energy | M1 | M2 | M3
---|---|---|---
E1 | 22.73 | 8.56 | 3.51
E2 | 5.95 | 12.32 | 10.88
E3 | 7.81 | 3.51 | 27.77
Table 1: Attenuation coefficients $c_{11}$, $c_{12}$, … of the demonstration
materials, with increasing photon energies E1, E2 and E3.

We chose resolution $N=128$ for the reconstructions. The number of used
projection images was sparse in these simulations: 15 equispaced angles.

After choosing the most optimal $\alpha$ value, (here $\alpha$ = 20000) we
start to increase the $\beta$ parameter from zero, to see how it affects the
result. See figure 8. The errors of the reconstructions, compared with the
ground truth are listed in the table 2. We can see from the results that
increasing the $\beta$ term reduces the error and improves the quality of the
reconstructions. The range of $\beta$ is limited to be smaller than the
$\alpha$ parameter to ensure that the problem stays convex.

Mean errors: | L2 | SSIM | haarPSI | Corr.pix % | $\alpha$ | $\beta$
---|---|---|---|---|---|---
First row: | 0.60 | 0.25 | 0.14 | 96 | 20000 | 0
Second row: | 0.56 | 0.26 | 0.17 | 96.4 | 20000 | 8000
Third row: | 0.55 | 0.39 | 0.19 | 96.3 | 20000 | 16000
Fourth row: | 0.59 | 0.41 | 0.21 | 95.6 | 20000 | 19500
Ground truth: | 0 | 1 | 1 | 100 |  |
Table 2: Error table for different $\beta$ values. We can see that the
$L_{2}$-errors12 are getting smaller, structural similarity is increasing and
also the HaarPSI index is improving, when we increase the $\beta$
regularization. We can not increase $\beta$ more than this, because the
problem might turn to unstable. Unfortunately the $\beta$ regularization does
not reduce the random noise spots around the actual target, which is natural
consequence if one thinks how the term works. It penalises only the
overlapping materials in different images and random errors in different
images rarely overlap.

Figure 7: Beta demonstration. First row: $\beta$=0, second row: $\beta$=8000,
third row: $\beta$=16000, fourth row: $\beta$ = 19500. We may notice from the
figures that when we increase the $\beta$ regularization, the materials
sharpen their shape in their own images, even though the effect is not very
strong.

Figure 8: Beta demonstration. First row: $\beta$=0, second row: $\beta$=8000,
third row: Ground truth. We may notice from the figures that when we increase
the $\beta$ regularization, the materials sharpen their shape in their own
images, even though the effect is not very strong.

###  5.3 Results with simulated data

We have simulated the real data case virtually by using the measured
attenuation values of the selenium compounds 4. We have the error measures
separately for each material in table 3.

Figure 9: Results with simulated data. The first row is the simulated result for the different selenium samples. Here we have totally ignored the plastic container, which was present in the practical real world measurements. The second row shows the perfectly separated materials, which we have used as ground truths in the error calculations. Table 3: In this table, we have error calculations for the simulated phantom. We are aiming for as small as possible $L_{2}$-error12, but for SSIM and HPSI, the closer the value gets to one, the better, because it means that the images are more structurally similar or visually alike, which is the goal, when comparing the simulated results with the actual ground truth. Material | Method |  mmL2mm | SSIM | HPSI
---|---|---|---|---
Se | IP | 0.27 | 0.87 | 0.42
SeO3 | IP | 0.25 | 0.81 | 0.47
SeO4 | IP | 0.20 | 0.89 | 0.54

###  5.4 Results with experimental X-ray data

In this section, we show the reconstruction results with the standard filtered
back projection (FBP) method and with the IP-method, see figure 10.

The number of used projection images is 200 with the real data. We used the
X-ray attenuation coefficients determined in [8] in our system matrix, in
place of the $c_{ij}$ coefficients. The X-ray attenuation coefficients used
with different X-ray energies are listed in Table 4.

Energy [keV] | Se | SeO3 | SeO4
---|---|---|---
12.658 | 8.4734 | 2.4380 | 3.5091
12.662 | 7.9902 | 11.859 | 10.885
12.685 | 7.4175 | 8.0033 | 27.777
Table 4: The X-ray attenuation coefficients for the different materials at
different energies. The coefficients were determined in the article of
Honkanen & al.[8].

Our simulation studies reported in Section 5.3 suggest that $\alpha=20000$ and
$\beta=19500$ are suitable choices for the regularization parameters.

The reconstruction images in figure 10 have been black point and white point
corrected for better visual impression. In these reconstructions, we have used
all the projection images taken, so we have used data from 200 different
projection angles. The first row of the figure shows FBP reconstructions made
from the sinograms, which were measured with energies 12,658 keV, 12,662 keV
and 12,685 keV. The end result is quite noisy and there are only small
differences between the reconstructions. We have the PMMA phantom as a fourth
material on the first row (and there was also some potassium as fulfilling
material). The attenuation coefficients of these "extra" materials are quite
same with every measured energy (they have no K-edges, see figure 6). This
for, we subtracted one low energy sinogram 12,645 keV (which was below all the
K-edges) from the other sinograms before reconstruction. The result can be
seen on the second and third row in the figure, where we have got rid of the
phantom matrix. We have the FBP reconstructions on the second row and the IP
reconstructions of the different materials on the third row.

FBP with energy 1Material 1FBP with energy 2Material 2FBP with energy
3Material 3(a)(b)(c) Figure 10: Reconstruction results. Row (a) represents the
FBP reconstructions, with the phantom matrix. Row (b) shows the result with
the FBP method, when we have subtracted one measurement (taken below all the
K-edges) from all sinograms before reconstruction. This subtraction removes
the PMMA phantom matrix from the reconstructions. Row (b): IP reconstructions
(also with subtracted phantom matrix).

##  6 Conclusions

The Inner Product reconstruction approach seems to be capable for material
decomposition in the case of real, monochromatic X-ray data. The
reconstruction results show clearly where the different materials are located
in the sample, even though there is some noticeable dotted noise around the
samples. The theory of this method has been shown to work for multiple
materials and X-ray energies and in this article we demonstrated the working
in the case of three energies and three materials. As a drawback, we can
notice that there must be clear difference in the attenuation coefficients of
the materials with different energies for the proper working of the method.
This is achieved using monochromatic X-ray beam and by selecting the used
X-ray energies wisely just below and over the K-edge peak of the materials to
be decomposed. When we keep these limitations in mind, the method works as a
reliable and simple tool for material separation and identification.

## References

  * [1] Kees Joost Batenburg, Sara Bals, J Sijbers, C Kübel, PA Midgley, JC Hernandez, U Kaiser, ER Encina, EA Coronado, and G Van Tendeloo.  3d imaging of nanomaterials by discrete tomography.  Ultramicroscopy, 109(6):730–740, 2009.
  * [2] Kees Joost Batenburg and Jan Sijbers.  Dart: a practical reconstruction algorithm for discrete tomography.  IEEE Transactions on Image Processing, 20(9):2542–2553, 2011.
  * [3] Joachim Baumann, Zoltán Kiss, Sven Krimmel, Attila Kuba, Antal Nagy, Lajos Rodek, Burkhard Schillinger, and Jürgen Stephan.  Discrete tomography methods for nondestructive testing.  In Advances in discrete tomography and its applications, pages 303–331. Springer, 2007.
  * [4] J. Gondzio.  Interior point methods 25 years later.  European Journal of Operational Research, 218:587–601, 2012.
  * [5] Jacek Gondzio, Matti Lassas, Salla-Maaria Latva-Äijö, Samuli Siltanen, and Filippo Zanetti.  Material-separating regularizer for multi-energy x-ray tomography.  Inverse Problems, 38(2):025013, 2022.
  * [6] Gabor T Herman and Attila Kuba.  Advances in discrete tomography and its applications.  Springer Science & Business Media, 2008.
  * [7] Gabor T Herman and Attila Kuba.  Discrete tomography: Foundations, algorithms, and applications.  Springer Science & Business Media, 2012.
  * [8] Ari-Pekka Honkanen and Simo Huotari.  Monochromatic computed tomography using laboratory-scale setup.  Scientific Reports, 13(1):363, 2023.
  * [9] Yong Long and Jeffrey A Fessler.  Multi-material decomposition using statistical image reconstruction for spectral ct.  IEEE transactions on medical imaging, 33(8):1614–1626, 2014.
  * [10] Rafael Reisenhofer, Sebastian Bosse, Gitta Kutyniok, and Thomas Wiegand.  A haar wavelet-based perceptual similarity index for image quality assessment.  Signal Processing: Image Communication, 61:33–43, 2018.
  * [11] EB Saloman, JH Hubbell, and JH Scofield.  X-ray attenuation cross sections for energies 100 ev to 100 kev and elements z= 1 to z= 92.  Atomic Data and Nuclear Data Tables, 38(1):1–196, 1988.
  * [12] Rasmus Solem, Till Dreier, Isabel Gonçalves, and Martin Bech.  Material decomposition in low-energy micro-ct using a dual-threshold photon counting x-ray detector.  Frontiers in Physics, 9:673843, 2021.
  * [13] Zhou Wang, Alan C Bovik, Hamid R Sheikh, and Eero P Simoncelli.  Image quality assessment: from error visibility to structural similarity.  IEEE transactions on image processing, 13(4):600–612, 2004.
  * [14] S. J. Wright.  Primal-Dual Interior-Point Methods.  SIAM, 1997.
  * [15] Weiwen Wu, Haijun Yu, Peijun Chen, Fulin Luo, Fenglin Liu, Qian Wang, Yining Zhu, Yanbo Zhang, Jian Feng, and Hengyong Yu.  Dictionary learning based image-domain material decomposition for spectral ct.  Physics in Medicine & Biology, 65(24):245006, 2020.