# Enhancing sparse representation of color images by cross channel
transformation

Laura Rebollo-Neira
Mathematics Department
Aston University
B4 7ET, Birmingham, UK
Aurelien Inacio
ENSIIE, 1 Rue de la Résistance,
91000 Évry-Courcouronnes, France

###### Abstract

Transformations for enhancing sparsity in the approximation of color images by
2D atomic decomposition are discussed. The sparsity is firstly considered with
respect to the most significant coefficients in the wavelet decomposition of
the color image. The discrete cosine transform is singled out as an effective
transformation for this purpose. The enhanced feature is further exploited by
approximating the transformed arrays using an effective greedy strategy with a
separable highly redundant dictionary. The relevance of the achieved sparsity
is illustrated by a simple encoding procedure. On a set of typical test images
the compression at high quality recovery is shown to significantly improve
upon JPEG and WebP formats. The results are competitive with those produced by
the JPEG2000 standard.

##  2 Introduction

In the signal processing field sparse representation usually refers to the
approximation of a signal as superposition of elementary components, called
atoms, which are members of a large redundant set, called a dictionary [1].
The superposition, termed atomic decomposition, aims at approximating the
signal involving as few atoms as possible [1, 3, 2, 4]. Sparsity is also
relevant to data collection. Within the emerging theory of sampling known as
compressive sensing (CS) this property is key for reconstructing signals from
a reduced number of measures [5, 6, 7]. In particular, distributed compressive
sensing (DCS) algorithms exploit inter signal correlation structure for
multiple signal recovery [8, 9]

Sparse image representation using redundant dictionaries has been considered
in numerous works e.g. [11, 12, 13] and in the context of applications such as
image restoration [14, 15, 16], feature extraction [17, 18, 19, 20] and super
resolution [21, 22, 23, 24].

The sparsity property of some types of 3D images benefits from 3D processing
[25, 26, 27, 28]. In particular most RGB (Red-Green-Blue) images admit a
sparser atomic decomposition if approximated by 3D atoms [29]. Within the 3D
framework the gain in sparsity comes at expenses of computational cost though.

The purpose of this work is to show that the application of a transformation
across the direction of the colors improve sparsity in the wavelet domain
representation of the 2D color channels. The relevance of this feature is
demonstrated by a simple encoding procedure rendering good compression results
in comparison to the state of the art.

The contributions of the paper are listed below.

  * •

Highlighting of a distinctive feature concerning sparse representation of RGB
images. Namely, that transformations in the direction of the color channels
significantly improve the sparsity of the 2D atomic decomposition of the
image.

  * •

Demonstration that the discrete cosine transform (dct) can improve sparsity
with respect to the optimal decorrelation transform, i.e. that given by the
eigenvectors of the correlation matrix, which is thereby image dependent.

  * •

Showcasing of the relevance of sparsity to image compression with high quality
recovery. On a set of 15 typical test images the achieved compression is shown
to significantly improve upon JPEG and WebP formats. The results are
competitive with those produced by JPEG2000.

The work is organised as follows: Sec. 3 introduces the mathematical notation.
Sec. 4 compares several cross color transformations enhancing sparsity. A
numerical example on a large data set is used to illustrate the suitability of
the dct for that purpose. Sec. 5 demonstrates the gain in sparsity obtained by
the atomic decomposition of color images when the dct is applied across the
RGB channels. Sec. 6 illustrates the relevance of the approach to image
compression with high quality recovery. The conclusions are summarised in Sec.
7.

##  3 Mathematical Notation

Throughout the paper we use the following notational convention. $\mathbb{R}$
represents the set of real numbers. Boldface letters indicate Euclidean
vectors, 2D and 3D arrays. Standard mathematical fonts indicate components,
e.g., $\mathbf{d}\in\mathbb{R}^{N}$ is a vector of components
$d(i)\in\mathbb{R},\,i=1,\ldots,N$. The elements of a 2D array
$\mathbf{I}\in\mathbb{R}^{N_{x}\times N_{y}}$ are indicated as
$I(i,j),\,i=1,\ldots,N_{x},\,j=1,\ldots,N_{y}$ and the color channels of
$\mathbf{I}\in\mathbb{R}^{N_{x}\times N_{y}\times 3}$ as $I(:,:,z),\,z=1,2,3$.
The transpose of a matrix, $\mathbf{G}$ say, is indicated as
$\mathbf{G}^{\top}$. A set of, say $M$, color images is represented by the
arrays $\mathbf{I}\\{m\\}\in\mathbb{R}^{N_{x}\times N_{y}\times
3},m=1,\ldots,M.$

The inner product between arrays in $\mathbb{R}^{N_{x}\times N_{y}}$ is given
by the Frobenius inner product
$\left\langle\cdot,\cdot\right\rangle_{\mathrm{F}}$ as

$\left\langle\mathbf{G},\mathbf{I}\right\rangle_{\mathrm{F}}=\sum_{i=1}^{N_{x}}\sum_{j=1}^{N_{y}}G(i,j)I(i,j).$

Consequently, the Frobenius norm $\|\cdot\|_{\mathrm{F}}$ is calculated as

$\|\mathbf{G}\|_{\mathrm{F}}=\sqrt{\sum_{i=1}^{N_{x}}\sum_{j=1}^{N_{y}}G(i,j)^{2}}.$

The inner produce between arrays in $\mathbb{R}^{N}$ is given by the Euclidean
inner product $\left\langle\cdot,\cdot\right\rangle$ as

$\left\langle\mathbf{d},\mathbf{g}\right\rangle=\sum_{i=1}^{N}d(i)g(i).$

##  4 Cross color transformations

Given a color image $I(i,j,k),i=1,\dots,N_{x},j=1,\dots,N_{y},\,k=1,2,3$ the
processing of the 3 channels can be realised either in the pixel/intensity or
in the wavelet domain. Since the representation of most images is sparser in
the wavelet domain [30, 31, 32, 33, 29] we approximate in that domain and
reconstruct the approximated image by the inverse wavelet transform. Thus,
using a $3\times 3$ matrix $\mathbf{T}$, we construct the transformed arrays
$\mathbf{U}$ and $\mathbf{W}$ as follows

$\displaystyle U(:,:,z)$ $\displaystyle\\!\\!=$
$\displaystyle\\!\\!\sum_{l=1}^{3}I(:,:,l)T(l,z)\quad z=1,2,3.$ (1)
$\displaystyle W(:,:,z)$ $\displaystyle\\!\\!=\\!\\!$
$\displaystyle\mathbf{dwt}(U(:,:,z))\,\,\,\,\,\,\,\quad z=1,2,3,$ (2)

where $\mathbf{dwt}$ indicates the 2D wavelet transform. For the
transformation $\mathbf{T}$ we consider the following cases

  * (i)

The dct.

  * (ii)

The reversible YCbCr color space transform.

  * iii)

The principal components (PC) transform.

  * iv)

A transformation learned from an independent set of images.

The dct is given by the matrix

$\begin{pmatrix}\frac{1}{\sqrt{3}}&\frac{\sqrt{2}}{\sqrt{3}}\cos(\frac{\pi}{6})&\frac{\sqrt{2}}{\sqrt{3}}\cos(\frac{\pi}{3})\\\
\frac{1}{\sqrt{3}}&0&\frac{\sqrt{2}}{\sqrt{3}}\cos(\pi)\\\
\frac{1}{\sqrt{3}}&\frac{\sqrt{2}}{\sqrt{3}}\cos(\frac{5\pi}{6})&\frac{\sqrt{2}}{\sqrt{3}}\cos(\frac{5\pi}{6})\end{pmatrix}$

The YCbCr color space transform is given by the matrix below [35]

$\begin{pmatrix}0.299&-0.169&0.5\\\ 0.587&-0.331&-0.419\\\
0.114&0.5&-0.0813\\\ \end{pmatrix}$

The columns of the principal components transform are the normalised
eigenvectors of covariance matrix of the RGB pixels. Thus, the transformation
is image dependent and optimal with respect to decorrelating the color
channels.

As a first test, the approximation of the transformed channels is realised by
keeping a fixed number of the largest absolute value entries, and setting the
others equal to zero. In relation to this, for an image of size $L_{x}\times
L_{y}\times 3$ we define the Sparsity Ratio (SR) as follows

$\text{SR}=\frac{\text{$L_{x}\cdot L_{y}\cdot 3$}}{\text{Number of nonzero
entries in the three channels}.}$ (3)

The quality of a reconstructed image $\tilde{\mathbf{I}}$, with respect to the
original 8-bit image $\mathbf{I}$, is compared using the Peak Signal-to-Noise
Ratio ($\mathrm{PSNR}$),

$\mathrm{PSNR}=10\log_{10}\left(\frac{255^{2}}{\mathrm{MSE}}\right),$ (4)

where

$\mathrm{MSE}=\frac{1}{L_{x}\cdot L_{y}\cdot
3}\sum_{i,j,z=1}^{L_{x},L_{y},3}(I(i,j,z)-\tilde{I}(i,j,z))^{2}.$

For the numerical examples below the transformation corresponding to case (iv)
is learned from a set of images ${\mathbf{I}}\\{m\\},\,m=1,\ldots,M$ all of
the same size.Starting from an invertible $3\times 3$ matrix $\mathbf{T}^{k}$,
with $k=1$, the learning algorithm proceeds through the following
instructions.

  * 1)

Use $\mathbf{T}^{k}$ to transform the arrays
${\mathbf{I}}\\{m\\}\to\mathbf{U}^{k}\\{m\\}\to\mathbf{W}^{k}\\{m\\}$ as in
(1) and (2).

  * 2)

Approximate each transformed array $\mathbf{W}^{k}\\{m\\}$ to obtain
$\tilde{\mathbf{W}}^{k}\\{m\\}$ by keeping the largest $K$ absolute value
entries.

  * 3)

Apply the inverse 2D wavelet transform $\mathbf{idwt}$ to reconstruct the
approximated arrays $\tilde{\mathbf{U}}^{k}\\{m\\},\,m=1,\ldots,M$ as

$\displaystyle\tilde{U}^{k}\\{m\\}(:,:,z)$ $\displaystyle\\!\\!\\!=\\!\\!\\!$
$\displaystyle\mathbf{idwt}(\tilde{W}^{k}\\{m\\}(:,:,z)),\,z=1,2,3.$

  * 4)

Use the original images ${\mathbf{I}}\\{m\\},\,m=1,\ldots,M$ to find
$\mathbf{G}=\mathbf{T}^{-1}$ by least squares fitting, i.e

$\mathbf{G}^{\ast}=\operatorname*{arg\,max}_{\begin{subarray}{c}G(z,l)\\\
z,l=1,2,3\end{subarray}}{\mathcal{E}}(\mathbf{G}),$

where

${\mathcal{E}}(\mathbf{G})=\sum_{m=1}^{M}\sum_{\begin{subarray}{c}i,j,l=1\end{subarray}}^{L_{x},L_{y},3}(I\\{m\\}(i,j,l)-\tilde{I}\\{m\\}(i,j,l))^{2}$

with

$\tilde{I}\\{m\\}(i,j,l))=\sum_{z=1}^{3}\tilde{U}^{k}\\{m\\}(i,j,z)G(z,l).$

  * 5)

While $\cal{E}$ decreases, or the maximum number of allowed iterations has not
been reached, set $k\to k+1$, $\mathbf{T}^{k}=(\mathbf{G}^{\ast})^{-1}$ and
repeat steps 1) – 5).

Given the arrays $\tilde{\mathbf{U}}^{k}\\{m\\},\,m=1,\ldots,M$ the least
squares problem for determining the transformation $\mathbf{T}^{k}$ has a
unique solution. However, the joint optimisation with respect to the arrays
$\tilde{\mathbf{U}}^{k}\\{m\\},m=1,\ldots,M$ and the transformation
$\mathbf{T}^{k}$ is not convex. Hence, the algorithm’s outcome depends on the
initial value.

The transformation (iv) used in the numerical examples of Secs. 4.1 and 5.1
has been learned from $M=100$ images, all of size $384\times 512\times 3$,
from the UCID data set [34], which contains images of buildings, places and
cars. The learning curves for two random orthonormal initialisations is shown
in Fig. 1.

Figure 1: Learning curves for the transformation (iv) corresponding to two
different random orthogonal transforms initialising the process. The mean
value PSNR with respect to the 100 images in the training set corresponds to
SR=15 for all the images.

It is worth mentioning that, as shown in Fig. 1, learning is richer when
starting comparatively distant from a local minimum (left graph in Fig. 1).
However, since the convergence is to a local minimum other random
initialisations, even if generating less learning, may give better results
(right graph in Fig. 1).

The aim of the next numerical test is to demonstrate the effect on the SR
(c.f. (3)) produced by the above (i) - (iv) transformations across the color
channels.

###  4.1 Numerical Test I

Using the whole Berkeley data set [36], consisting of 300 images all of size
$321\times 481\times 3$ we proceed with each image as in (1) and (2). The
$\mathbf{dwt}$ corresponds to the Cohen-Daubechies-Feauveau 9$/$7 (CDF97)
wavelet family. Fig. 3 shows the transformed channels of the image in Fig.2,
including the dct transformation across channels (right graph) and without
$\mathbf{T}$ transformation (left graph).

Figure 2: One of the RGB images in the Berkeley’s data set Figure 3: Magnitude
of the entries in the array $\mathbf{W}$ constructed as in (2) from the image
of Fig. 2, with $\mathbf{T}$ the dct (right graph) and without $\mathbf{T}$
(left graph).

For the numerical test the approximations are realised fixing SR=20 and SR=10.
The reconstructed images are obtained for the approximated arrays $\tilde{W}$
as

$\displaystyle\tilde{U}(:,:,z)$ $\displaystyle=$
$\displaystyle\mathbf{idwt}(\tilde{W}(:,:,z)),\quad\;\;\;\;\;\;\;z=1,2,3$ (5)
$\displaystyle\tilde{I}(:,:,z)$ $\displaystyle=$
$\displaystyle\\!\\!\\!\sum_{k=1}^{3}\tilde{U}(:,:,k)(T^{-1}(z,k)),\;z=1,2,3$
(6)

where $\mathbf{T}^{-1}$ is the inverse of $\mathbf{T}$. When no transformation
$\mathbf{T}$ is considered the image is reconstructed directly from (5).

The 2nd and 4th columns of Table 1 show the mean value $\mathrm{PSNR}$
($\overline{\mathrm{PSNR}}$), with respect to the whole data set, for SR=20
and SR=10 respectively, corresponding to the transformations $\mathbf{T}$
listed in the 1st column of the table. The 3rd and 5th columns give the
standard deviations (std).

Transf. | $\overline{\mathrm{PSNR}}$ | std | $\overline{\mathrm{PSNR}}$ | std
---|---|---|---|---
(i) dct | 33.9 | 5.0 | 38.9 | 5.1
(ii) YCbCr | 33.7 | 5.0 | 38.7 | 5.1
(iii) PC | 33.7 | 4.9 | 38.4 | 5.0
(iv) Learned | 33.7 | 4.9 | 38.7 | 5.0
(v) No transf. | 29.7 | 4.8 | 32.9 | 5.2
Table 1: Mean value PSNR, with respect to the 300 images in the Berkeley data
set. The approximation of each image is realised by setting the least
significant entries in the arrays $\mathbf{W}\\{m\\}=1,\ldots,300$ equal zero,
in order to obtain SR=20 for all the images (2nd column) and SR=10 for all the
images (4th column).

While Table 1 shows that all (i) - (iv) transformations render equivalent
superior results, with respect to not applying a transformation (case (v)),
the dct slightly exceeds the others. Case (iv) refers to the best results
achieved when initialising the learning algorithm with 500 different random
orthonormal transformations. When initialising the algorithm with
transformations (i) - (iii) it appears that each of these transformations is
close to a local minimiser of the algorithm. This stems from the fact that
such initialisations do not generate significant learning.

The common feature of most of the 300 images in the data set used in this
numerical example is the correlation property of the three color channels.
This property is assessed by the correlation coefficients

$\displaystyle r(z)$ $\displaystyle=$
$\displaystyle\frac{\sum_{i=1}^{L_{x}}\sum_{j=1}^{L_{y}}\Gamma(i,j,z)\Gamma(i,j,z+1)}{\sigma(z)\sigma(z+1)},\quad
z=1,2,$ $\displaystyle r(3)$ $\displaystyle=$
$\displaystyle\frac{\Gamma(1)\Gamma(3)}{\sigma(1)\sigma(3)},$

where

$\displaystyle\Gamma(i,j,z)$ $\displaystyle\\!\\!=\\!\\!$
$\displaystyle(I(i,j,z)-\overline{I(:,:,z)}),$ (7) $\displaystyle\sigma(z)$
$\displaystyle\\!\\!=$
$\displaystyle\\!\\!\sqrt{\sum_{i=1}^{L_{x}}\sum_{j=1}^{L_{y}}(I(i,j,z)-\overline{I(:,:,z)})^{2}},$
(8)

and $\overline{I(:,:,z)}$ indicates the mean value of channel $z$.

As seen in the 3 histograms of Fig. 4 for most images in the data set the
correlation of the color channels is high. It is surprising then than the PC
transformation (iii), which completely decorrelates the channels, does not
overperform the other transformations, on the contrary. This feature has also
been noticed in the context of bit allocation for subband color image
compression [37].

Figure 4: Histograms of the correlation coefficients between the RGB channels.

##  5 Approximations by atomic decomposition

We have seen (c.f. Table 1) that by transformation of channels it is possible
to gain quality when reducing nonzero entries in the channels. Now we discuss
how to improve quality further by approximating the 2D arrays (2) by an atomic
decomposition, other than just by neglecting their less significant entries.
For the approximation to be successful it is important to use an appropriate
dictionary. To this end, one possibility could be to learn the dictionary from
training data [38, 39, 40, 41, 42, 43, 44, 45]. However as demonstrated in
previous works [32, 33, 29] a separable dictionary, which is easy to
construct, is well suited for the purposes of achieving sparsity and delivers
a fast implementation of the approach. Since we use that dictionary in the
numerical examples, below we describe the method for constructing the atomic
decomposition of the array $\mathbf{W}$ considering specifically a separable
dictionary.

Firstly we concatenate the 3 planes $W(:,:,z),z=1,2,3$ into an extended 2D
array $\mathbf{W}^{\prime}\in\mathbb{R}^{3L_{x}\times L_{y}}$ and divide this
array in small blocks $\mathbf{W}^{\prime}_{q},\,q=1,\ldots,Q$. Without loss
of generality the blocks are assumed to be square of size $N_{b}\times N_{b}$
say, and are approximated using separable dictionaries
$\mathcal{D}^{x}=\\{\mathbf{d}^{x}_{n}\in\mathbb{R}^{N_{b}},\,\|\mathbf{d}^{x}_{n}\|_{2}=1\\}_{n=1}^{M_{x}}$
and
$\mathcal{D}^{y}=\\{\mathbf{d}^{y}_{m}\in\mathbb{R}^{N_{b}},\,\|\mathbf{d}^{y}_{n}\|_{2}=1\\}_{m=1}^{M_{y}}$,
with $M_{x}M_{y}=M$.

For $q=1,\ldots,Q$ every element ${\mathbf{W}^{\prime}}_{q}$ is approximated
by an atomic decomposition as below:

${\mathbf{W}^{\prime}}_{q}^{k_{q}}=\sum_{n=1}^{k_{q}}c^{k_{q},q}(n)\mathbf{d}^{x}_{\ell^{x,q}_{n}}(\mathbf{d}^{y}_{\ell^{y_{,}q}_{n}})^{T},$
(9)

where ${\ell^{y_{,}q}_{n}}$ is the index in the set $\\{1,2,\ldots,M_{y}\\}$
corresponding to the label of the atom in the dictionary $\mathcal{D}^{y}$
contributing to the $n$-th term in the approximation of the $q$-th block. The
index ${\ell^{x_{,}q}_{n}}$ has the equivalent description. The assembling of
the approximated blocks gives rise to the approximated array
$\mathbf{W}^{\prime K}$, i.e.,
${\mathbf{W}^{\prime}}^{K}=\hat{\operatorname{J}}_{q=1}^{Q}{\mathbf{W}^{\prime}}_{q}^{k_{q}}$,
where $K=\sum_{q=1}^{Q}k_{q}$ and $\hat{\operatorname{J}}$ represents the
assembling operation.

The approximation of the partition ${\mathbf{W}^{\prime}}_{q}=1,\ldots,Q$ is
carried out iteratively as a two step process which selects i) the atoms in
the atomic decomposition (9) and ii) the sequence in which the blocks in the
partition are approximated. The procedure is called Hierarchized Block Wise
(HBW) implementation of greedy pursuit strategies [31, 46]. For the selection
of the atoms we apply the Orthogonal Matching Pursuit (OMP) approach [47]
dedicated to 2D with separable dictionaries (OMP2D) [48]. Thus, the whole
algorithm is termed HBW-OMP2D [31]. The method iterates as described below.

On setting $k_{q}=0$ and $\mathbf{R}_{q}^{0}={\mathbf{W}^{\prime}}_{q}$ at
iteration $k_{q}+1$ the algorithm selects the indices $\ell^{x,q}_{k_{q}+1}$
and $\ell^{y,q}_{k_{q}+1}$, as follows:

$\ell^{x,q}_{k_{q}+1},\ell^{y,q}_{k_{q}+1}=\operatorname*{arg\,max}_{\begin{subarray}{c}n=1,\ldots,M_{x}\\\
m=1,\ldots,M_{y}\end{subarray}}\left|\left\langle\mathbf{d}^{x}_{n},\mathbf{R}_{q}^{k_{q}}\mathbf{d}^{y}_{m}\right\rangle_{F}\right|,\,\,\,q=1,\ldots,Q,$
(10)

where
$\mathbf{R}_{q}^{k_{q}}={\mathbf{W}^{\prime}}_{q}-{\mathbf{W}^{\prime}}_{q}^{k_{q}}$.
The calculation of ${\mathbf{W}^{\prime}}_{q}^{k_{q}}$ is realised in order to
minimise $\|\mathbf{R}_{q}^{k_{q}}\|_{F}$, which is equivalent to finding the
orthogonal projection onto the subspace spanned by the selected
atoms$\\{\mathbf{A}_{n}=\mathbf{d}^{x}_{\ell^{x,q}_{n}}(\mathbf{d}^{y}_{\ell^{y,q}_{n}})^{T}\\}_{n=1}^{k_{q}}$.
In our implementation, the calculation of the coefficients
$c^{q}(n),\,n=1,\ldots,k_{q}$ in (9) is realised as

$c^{k_{q},q}(n)=\left\langle\mathbf{B}_{n}^{k},{\mathbf{W}^{\prime}}_{q}^{k_{q}}\right\rangle_{\mathrm{F}},\quad
n=1,\ldots,k_{q},$ (11)

where the set $\\{\mathbf{B}_{n}^{k_{q}}\\}_{n=1}^{k_{q}}$ is biorthogonal to
the set $\\{\mathbf{A}_{n}\\}_{n=1}^{k_{q}}$ and needs to be upgraded and
updated to account for each newly selected atom. Starting from
$\mathbf{B}_{1}^{1}=\mathbf{Q}_{1}=\mathbf{A}_{1}=\mathbf{d}^{x}_{\ell^{x}_{1}}(\mathbf{d}^{y}_{\ell^{y}_{1}})^{\top}$
the updating and upgrading is realised through the recursive equations [48,
49]:

$\begin{split}\mathbf{B}_{n}^{k_{q}+1}&=\mathbf{B}_{n}^{k_{q}}-\mathbf{B}_{k_{q}+1}^{k_{q}+1}\left\langle\mathbf{A}_{k_{q}+1},\mathbf{B}_{n}^{k}\right\rangle_{\mathrm{F}},\quad
n=1,\ldots,k_{q}\\\
\mathbf{B}_{k_{q}+1}^{k_{q}+1}&=\mathbf{Q}_{k_{q}+1}/\|\mathbf{Q}_{k_{q}+1}\|_{\mathrm{F}}^{2},\,\,\text{where}\\\
\mathbf{Q}_{k_{q}+1}&=\mathbf{A}_{k_{q}+1}-\sum_{n=1}^{k_{q}}\frac{\mathbf{Q}_{n}}{\|\mathbf{Q}_{n}\|_{\mathrm{F}}^{2}}\left\langle\mathbf{Q}_{n},\mathbf{A}_{k_{q}+1}\right\rangle_{\mathrm{F}},\end{split}$
(12)

with the additional re-orthogonalisation step

$\mathbf{Q}_{k_{q}+1}\leftarrow\mathbf{Q}_{k_{q}+1}-\sum_{n=1}^{k_{q}}\frac{\mathbf{Q}_{n}}{\|\mathbf{Q}_{n}\|_{\mathrm{F}}^{2}}.\left\langle\mathbf{Q}_{n},\mathbf{Q}_{k_{q}+1}\right\rangle_{\mathrm{F}}.$
(13)

As discussed in [31, 46], for $\ell^{x,q}_{k_{q}+1}$ and
$\ell^{y,q}_{k_{q}+1},\,q=1,\ldots,Q$ the indices resulting from (10), the
block to be approximated in the next iteration corresponds to the value
$q^{\star}$ such that

$q^{\star}=\operatorname*{arg\,max}_{q=1,\ldots,Q}\left|\left\langle\mathbf{d}^{x}_{\ell^{x,q}_{k_{q}+1}},\mathbf{R}_{q}^{k_{q}}\,\mathbf{d}^{y}_{\ell^{y,q}_{k_{q}+1}}\right\rangle\right|.$

The algorithm stops when the required total number of $K$ atoms has been
selected. This number can be fixed using the SR, which is now calculated as

$\text{SR}=\frac{\text{Number of elements in the arrays}}{K}.$ (14)

######  Remark 1.

The above described implementation of HBW-OMP2D is very effective in terms of
speed, but demanding in terms of memory (the partial outputs corresponding to
all the blocks in the partition need to be stored at every iteration). An
alternative implementation, termed HBW Self Projected Matching Pursuit (HBW-
SPMP) [32, 50], would enable the application of the identical strategy to much
larger images than the ones considered in this work.

###  5.1 Numerical example II

For this and the next numerical example, we use a mixed dictionary consisting
of two classes of sub-dictionaries of different nature:

  * I)

The trigonometric dictionaries $\mathcal{D}_{C}^{x}$ and
$\mathcal{D}_{S}^{x}$, defined below,

$\displaystyle\mathcal{D}_{C}^{x}$
$\displaystyle\\!\\!\\!\\!=\\!\\!\\!\\!\\!\\!$
$\displaystyle\\{w_{c}(n)\cos{\frac{{\pi(2i-1)(n-1)}}{2M}},i=1,\ldots,N_{b}\\}_{n=1}^{M}$
$\displaystyle\mathcal{D}_{S}^{x}$
$\displaystyle\\!\\!\\!\\!=\\!\\!\\!\\!\\!\\!$
$\displaystyle\\{w_{s}(n)\sin{\frac{{\pi(2i-1)(n)}}{2M}},i=1,\ldots,N_{b}\\}_{n=1}^{M},$

where $w_{c}(n)$ and $w_{s}(n),\,n=1,\ldots,M$ are normalisation factors.

  * II)

The dictionary $\mathcal{D}_{L}^{x}$, which is constructed by translation of
the prototype atoms in Fig. 5.

The mixed dictionary $\mathcal{D}^{x}$ is built as
$\mathcal{D}^{x}=\mathcal{D}_{C}^{x}\cup\mathcal{D}_{S}^{x}\cup\mathcal{D}_{L}^{x}$
and $\mathcal{D}^{y}=\mathcal{D}^{x}$.

Figure 5: Prototypes (each in different color) that generate the dictionaries
$\mathcal{D}_{L}^{x}$ by sequential translations of one point.

Table 2 shows the improvement in $\overline{\mathrm{PSNR}}$ achieved by atomic
decompositions using the mixed dictionary for SR= 20 and SR= 10.

Transf. | $\overline{\mathrm{PSNR}}$ | std | $\overline{\mathrm{PSNR}}$ | std
---|---|---|---|---
(i) dct | 40.5 | 5.0 | 48.1 | 4.4
(ii) YCbCr | 40.3 | 5.0 | 47.8 | 4.4
(iii) PC | 39.6 | 4.8 | 46.3 | 4.5
(iv) Learned | 40.3 | 4.9 | 47.4 | 4.3
(v) No transf. | 34.0 | 5.2 | 39.1 | 5.6
Table 2: Mean value PSNR, with respect to 300 images in the Berkeley data set,
produced by 2D atomic decomposition of the arrays
$\mathbf{W}\\{m\\}=1,\ldots,300$ in order to obtain SR=20 (2nd column) and
SR=10 (4th column).

Notice that while case (v), which does not include any $\mathbf{T}$
transformation, gives superior results than by disregarding entries (c.f.
Table 1) when applying any of the transformations (i) – (iv) the results
improve further. The PC transform, however, appears significantly less
effective than the others. In addition to rendering the best results, the dct
brings along the additional advantage of being orthonormal. Consequently, it
does not magnify errors at the inversion step. Because of this, we single out
the dct as the most convenient cross color transformation out of the four
considered here.

##  6 Application to image compression

In order achieve compression by filing an atomic decomposition we need to
address two issues. Namely, the quantisation of the coefficients
$c_{q}(n)\,n=1,\ldots,k_{q},\,q=1,\ldots,Q$ in (9) and the storage of the
indices $(\ell^{x,q}_{n},\ell^{y,q}_{n}),n=1,\ldots,k_{q},\,q=1,\ldots,Q$. We
tackle the matters by simple but effective procedures [33].

For $q=1,\ldots,Q$ the absolute value coefficients
$|c_{q}(n)|,\,n=1,\ldots,k_{q}$ are converted to integers through uniform
quantisation as follows

$c^{\Delta}_{q}(n)=\begin{cases}\lceil\frac{|c_{q}(n)|-\theta}{\Delta}\rceil,&\mbox{if}\quad|c_{q}(n)|\geq\theta\\\
&0\quad\mbox{otherwise}.\end{cases}$ (15)

The signs of the coefficient are encoded separately as a vector,
$\mathbf{s}_{q}$, using a binary alphabet. Each pair of indices
$(\ell^{x,q}_{n},\ell^{y,q}_{n})$ corresponding to the atoms in the
decompositions of the block ${\mathbf{W}}^{\prime}_{q}$ is mapped into a
single index $o_{q}(n)$. The set $o_{q}(1),\ldots,o_{q}(k_{q})$ is sorted in
ascending order $o_{q}(n)\rightarrow\tilde{o}_{q}(n),\,n=1,\ldots,k_{q}$ to
take the differences
$\delta_{q}(n)=\tilde{o}_{q}(n)-\tilde{o}_{q}(n-1),\,n=2,\ldots,k_{q}$ and
construct the string of non-negative numbers
${\tilde{o}_{q}(1),\delta_{q}(2),\ldots,\delta_{q}(k_{q})}$. The order of the
set $\tilde{o}_{q}(n),\,n=1,\ldots,k_{q}$ induces order in the unsigned
coefficients, $c^{\Delta}_{q}(n)\rightarrow\tilde{c}^{\Delta}_{q}(n)$, and in
the corresponding signs $s_{q}(n)\rightarrow\tilde{s}_{q}(n)$.

For each $q$ the number 0 is added at the end of the indices
$\tilde{o}_{q}(n),\,n=1,\ldots,k_{q}$ before concatenation, to be able to
separate strings corresponding to different blocks. Each sequence of strings
corresponding to $q=1,\ldots,Q$ is concatenated and encoded using the off-the-
shelf MATLAB function Huff06 [51], which implements Huffman coding.

The compression rate is given in bits-per-pixel (bpp) which is defined as

$\mathrm{bpp}=\frac{\text{Size of the file in bits}}{\text{Number of intensity
pixels in a single channel}}.$

At the reconstruction stage the indices
$(\tilde{\ell}^{x,q}_{n},\tilde{\ell}^{y,q}_{n}),\,n=1\,\ldots k_{q}$ are
recovered from the string of differences $\delta_{q}(n)),\,n=2,\ldots,k_{q}$.
The signs of the coefficients are read from the binary string. The quantised
unsigned coefficients are read and transformed into real numbers as:

$|\tilde{c}_{q}^{\mathrm{r}}(n)|=\Delta\cdot{\tilde{c}}_{q}^{\Delta}(n)+(\theta-\Delta/2)\,\quad
n=1\,\ldots k_{q}.$

The codec for reproducing the examples in the next sections has been made
available on [52].

###  6.1 Numerical Example III

The relevance to image compression of the achieved sparsity by dct cross color
transformation is illustrated in this section by comparison with results
yielded by the compression standards JPEG, WebP, and JPEG2000, on the 15
images in Table 3. These are typical images, used for compression tests,
available in ppm or png format. The first 9 images are classic test images
taken from [53]. The last six images are portions of $1024\times 1024\times 3$
pixels from very large high resolution images available on [54].

All the results have been obtained in the MATLAB environment (version R2019a),
using a machine with CPU Intel(R) Core(TM) i7-3520M RAM 8GB CPU @ 2.90GHz. For
the image approximation the HBW-OMP2D method was implemented with a C++ MEX
file. All the channels and all the images were partitioned into blocks of size
$16\times 16$. The approximation times to produce the results in Table 4 are
displayed in the last column of Table 3.

Figure 6: Illustration of the $1024\times 1024\times 3$ panels from the six
high resolution test images (No 10–15) listed in Table 3.

For realising the comparison we proceed as follows: we set the required value
of PSNR as that produced by JPEG at quality=95 and tune compression with the
other methods to produce the same PSNR. In our codec the tuning is realised by
approximating the image up to $\mathrm{PSNR}_{o}=1.025\cdot\mathrm{PSNR}$
(where $\mathrm{PSNR}$ is the targeted quality) and setting the quantisation
parameter $\Delta$ so as to reproduce the targeted quality. For compression
with JPEG and JPEG2000 we use the MATLAB imwrite command. The comparison with
WebP was realised using the software for Ubuntu distributed on [55].

No |  Image | Size | time (secs)
---|---|---|---
1 | Lenna | $512\times 512\times 3$ | 1.8
2 | Goldhill | $576\times 720\times 3$ | 3.8
3 | Barbara | $576\times 720\times 3$ | 3.2
4 | Baboon | $512\times 512\times 3$ | 2.7
5 | Zelda | $576\times 784\times 3$ | 2.5
6 | Sailboat | $512\times 512\times 3$ | 1.8
7 | Boy | $512\times 768\times 3$ | 3.8
8 | Jupiter | $1072\times 1376\times 3$ | 2.9
9 | Saturn | $1200\times 1488\times 3$ | 2.6
10 | Building | $1024\times 1024\times 3$ | 9.5
11 | Cathedral | $1024\times 1024\times 3$ | 7.6
12 | Flower | $1024\times 1024\times 3$ | 3.6
13 | Spider-web | $1024\times 1024\times 3$ | 3.8
14 | Bridge | $1024\times 1024\times 3$ | 8.5
15 | Deer | $1024\times 1024\times 3$ | 5.8
Table 3: Test Images. The last column gives the approximation times to produce the results in Table 4 I | dB | $\mathrm{bjp}$ | $\mathrm{bwb}$ | $\mathrm{bjp2}$ | $\mathrm{bsr}$
---|---|---|---|---|---
1 | 35.9 | 3.28 | 2.21 | 1.40 | 1.37
2 | 36.6 | 3.63 | 2.47 | 2.11 | 2.01
3 | 37.2 | 3.67 | 2.80 | 1.91 | 1.77
4 | 28.8 | 5.80 | 4.66 | 2.74 | 2.48
5 | 39.3 | 2.61 | 1.81 | 1.08 | 1.07
6 | 30.9 | 4.49 | 3.21 | 1.62 | 1.51
7 | 32.6 | 4.34 | 3.07 | 2.36 | 2.22
8 | 48.2 | 0.60 | 0.51 | 0.24 | 0.20
9 | 49.0 | 0.34 | 0.36 | 0.15 | 0.12
10 | 37.4 | 3.41 | 2.35 | 1.84 | 1.75
11 | 38.5 | 2.84 | 1.83 | 1.26 | 1.20
12 | 41.5 | 1.78 | 1.08 | 0.56 | 0.53
13 | 34.5 | 3.57 | 2.48 | 1.61 | 1.48
14 | 45.0 | 1.55 | 1.10 | 0.58 | 0.57
15 | 30.9 | 3.76 | 2.59 | 1.04 | 0.90
Table 4: Compression rate ($\mathrm{bpp}$) corresponding to JPEG
($\mathrm{bjp}$), WebP ($\mathrm{bwb}$), JPEG2000 ($\mathrm{bjp2}$) and the
proposed sparse representation ($\mathrm{bsr}$), for the values of PSNR given
in the $2\mathrm{nd}$ column.

##  7 Conclusions

The application of a cross color transformation for enhancing sparsity in the
atomic decomposition of RGB images has been proposed. It was demonstrated that
the effect of the transformation is to re-distribute the most significant
values in the dwt of the 2D channels. As a result, when approximating the
arrays by disregarding the less significant entries, the quality of the
reconstructed image improves with respect to disabling the cross color
transformation. Four transformations have been considered: (i) a 3 point dct,
(ii) the reversible YCbCr color space transform, (iii) the PC transform, (iv)
a transformation learned from an independent set of images.

The quality of the image approximation was improved further by approximating
the transformed arrays by an atomic decomposition using a separable dictionary
and the greedy pursuit strategy HBW-OMP2D. The dct was singled out as the most
convenient cross color transformation for approximating RGB color images in
the wavelet domain.

The approximation approach has been shown to be relevant for image
compression. By means of a simple coding strategy the achieved compression
considerably improves upon the compression standards JPEG and WebP. On a set
of 15 typical test images the results are at least as good as those produced
by JPEG2000. We feel confident that the findings communicated here will
motivate the consideration of dct as a cross color transformations for other
image processing applications which also benefit from the sparsity of a
representation.

## References

  * [1] S. Mallat and Z. Zhang, “Matching pursuit with time-frequency dictionaries,” IEEE Trans. Signal Process., Vol (41,12) 3397–3415 (1993).
  * [2] S. Chen, D. Donoho and M. Saunders “Atomic Decomposition by Basis Pursuit”, SIAM Review, 129–159 (2001).
  * [3] S. Mallat, A Wavelet Tour of Signal Processing: The Sparse Way, Academic Press (2009).
  * [4] Y. Eldar, P. Kuppinger and H. Biölcskei, “Block-Sparse Signals: Uncertainty Relations and Efficient Recovery”, IEEE Trans. Signal Process., 58, 3042–3054 (2010).
  * [5] E. Candès and M. Wakin,  “An introduction to compressive sampling”,  IEEE Signal Processing Magazine 25, 21 – 30 (2008).
  * [6] J. Romberg,  “Imaging via compressive sampling”,  IEEE Signal Processing Magazine 25, 14 – 20 (2008).
  * [7] R. Baraniuk,  “More Is less: Signal processing and the data deluge”,  Science 331, 717 – 719 (2011).
  * [8] M.F. Duarte, S. Sarvotham, D. Baron, M. B. Wakin, and R.G. Baraniuk, “Distributed compressed sensing of jointly sparse signals”, Proc. Asilomar Conf. Signals, Syst. Comput. 1537–1541 (2005).
  * [9] R. Torkamani, H. Zayyani, and R. Sadeghzadeh,  ‘Model-based decentralized Bayesian algorithm for distributed compressed sensing’, 95, (2021), 116212.
  * [10] Jua Xiao, Rui Zhao, and Kim-Man Lam, “Bayesian sparse hierarchical model for image denoising”, Signal Processing: Image Communication (2021), 116299
  * [11] J. Wright, Yi Ma, J. Mairal, G. Sapiro, T.S. Huang, and S. Yan,  “Sparse Representation for Computer Vision and Pattern Recognition”,  Proc. of the IEEE 98, 1031 – 1044 (2010).
  * [12] M. Elad, Sparse and Redundant Representations: From Theory to Applications in Signal and Image Processing, Springer (2010).
  * [13] Z. Zhang, Y. Xu, J. Yang, X. Li, and D. Zhang, “A survey of sparse representation: algorithms and applications”, IEEE access (2015).
  * [14] J. Mairal, M. Eldar, and G. Sapiro,  “Sparse Representation for Color Image Restoration”,  IEEE Trans. Image Proces., 17, 53 – 69 (2008).
  * [15] W. Dong, L. Zhang, G. Shi, and X. Li, “Nonlocally Centralized Sparse Representation for Image Restoration”,  IEEE Trans. Image Proces., 22, 1620 – 1630 (2013).
  * [16] Zhenming Su, Simiao Zhu, Xin Lv, and Yi Wan, “Image restoration using structured sparse representation with a novel parametric data-adaptive transformation matrix”,  Signal Processing: Image Communication, 52, 151–172 (2017).
  * [17] J. Wright, A. Y. Yang, and A. Ganesh, “Robust Face Recognition via Sparse Representation”,  IEEE Trans. Pattern Analysis and Machine Intelligence, 31, 210 – 227 (2009).
  * [18] XT. Yuan, X. Liu, and S. Yan, “Visual Classification With Multitask Joint Sparse Representation”,  IEEE Trans. Image Proces., 21, 4349 – 4360 (2012).
  * [19] D. Heinsohn, E. Villalobos, L. Prieto, and D. Mery “Face recognition in low-quality images using adaptive sparse representations”, Image and Vision Computing, 85, 46–58 (2019).
  * [20] K. Raja, R. Ramachandra, and C. Busch “Collaborative representation of blur invariant deep sparse features for periocular recognition from smartphones”. Image and Vision Computing, 101, 46–58 (2020).
  * [21] J. Yang, J. Wright, and T. Huang, “Image Super-Resolution via Sparse Representation,” IEEE Trans. Image Proces., 19, 2861 – 2873 (2010).
  * [22] Y. Zhang, J. Liu, W. Yang, and Z. Guo, “Image Super-Resolution Based on Structure-Modulated Sparse Representation”, IEEE Trans. Image Proces., 9, 2797 – 2810 (2015).
  * [23] Chaopeng Zhang, Weirong Liu, Jie Lui, Chaorong Liu, and Changhong Shi, “Sparse representation and adaptive mixed samples regression for single image super-resolution”, Signal Processing: Image Communication, 67, 79–89 (2018).
  * [24] Xuesong Li, Guo Cao, Youqiang Zhang, Ayesha Shafique, and Peng Fu, “Combining synthesis sparse with analysis sparse for single image super-resolution”, Signal Processing: Image Communication, 83 115805 (2020).
  * [25] C. F. Caiafa and A. Cichocki, “Computing sparse representations of multidimensional signals Using Kronecker Bases”, Neural computation, 25, 186 – 220 (2013).
  * [26] A. Cichocki, D. Mandic, L. De Lathauweri, G. Zhou, Q. Zhao, C. Caiafai, and H. A. Phan, “Tensor decompositions for signal processing applications: From two-way to multiway component analysis”, IEEE Signal Processing Magazine, 32, 145–163 (2015).
  * [27] Q. Dai, S. Yoo, and A. Kappeler “Sparse Representation-Based Multiple Frame Video Super-Resolution”, IEEE Trans. Image Proces., 26, 2080 – 2095 (2017).
  * [28] H. S. Mousavi and V. Monga, “Sparsity-based Color Image Super Resolution via Exploiting Cross Channel Constraints”, IEEE Trans. Image Proces., 26, 5094 – 5106 (2017).
  * [29] L. Rebollo-Neira and D. Whitehouse, “Sparse representation of 3D images for piecewise dimensionality reduction with high quality reconstruction”, Array, 1, doi = 10.1016/j.array.2019.100001 (2019).
  * [30] U. Grenander and A. Srivastava, “Probability models for clutter in natural images”, IEEE Trans. Pattern Anal. Mach. Intell., 23, 424–429 (2001) 10.1109/34.917579
  * [31] L. Rebollo-Neira, R. Matiol, and S. Bibi, “Hierarchized block wise image approximation by greedy pursuit strategies,” IEEE Signal Process. Letters, 20, 1175–1178 (2013).
  * [32] L. Rebollo-Neira, “ Effective sparse representation of X-Ray medical images”, International Journal for Numerical Methods in Biomedical Engineering, in press (2017).
  * [33] L. Rebollo-Neira, “A competitive scheme for storing sparse representation of X-Ray medical images”, PLoS ONE, 13(8): e0201455. https://doi.org/10.1371/journal.pone.0201455 (2018)
  * [34] G. Schaefer and M. Stich, “UCID: an uncompressed color image database”, Proceedings Volume 5307, Storage and Retrieval Methods and Applications for Multimedia 2004; (2003) https://doi.org/10.1117/12.525375.
  * [35] D. S. Taubman and M. W. Marcellin, JPEG2000 Image Compression, Fundamentals, Standards and Practice, Kluwer Academic Publishers (2002).
  * [36] https:www2.eecs.berkeley.edu/Research/Projects/CS/vision/bsds/  (Last access May 2021).
  * [37] E. Gershikov and M. Porat, “On color transforms and bit allocation for optimal subband image compression”, Signal Processing: Image Communication, 22, 1–18 (2007).
  * [38] K. Kreutz-Delgado, J. F. Murray, B. D. Rao, K. Engan, Te-Won Lee, and T. J. Sejnowski, “Dictionary Learning Algorithms for Sparse Representation”, Neurocomputing, 15, 349–396 (2003).
  * [39] M. Aharon, M. Elad, and A. Bruckstein, “K-SVD: An Algorithm for Designing Overcomplete Dictionaries for Sparse Representation”, IEEE Trans. Signal Process. 54, 4311–4322 (2006).
  * [40] R. Rubinstein, M. Zibulevsky, and M. Elad, “Double Sparsity: Learning Sparse Dictionaries for Sparse Signal Approximation”, IEEE Trans. Signal Process. 58, 1553–1564 (2010).
  * [41] I. Tos̆ić and P. Frossard, “Dictionary Learning: What is the right representation for my signal?”, IEEE Signal Processesing Magazine, 28, 27–38 (2011).
  * [42] J. Zepeda, C. Guillemot, and E. Kijak, “Image Compression Using Sparse Representations and the Iteration-Tuned and Aligned Dictionary”, IEEE Journal of Selected Topics in Signal Processing, 5,1061–1073 (2011).
  * [43] S. Hawe, M. Seibert, and M. Kleinsteuber, “Separable dictionary learning”, Proceedings of the 2013 IEEE Conference on Computer Vision and Pattern Recognition, 438–445 (2013).
  * [44] M. Srinivas, R. R. Naidu, C.S. Sastry, and C. Krishna Mohana, “Content based medical image retrieval using dictionary learning”, Neurocomputing, 168, 880–895 (2015).
  * [45] B. Wen, S. Ravishankar, and Y. Bresler, “‘Structured Overcomplete Sparsifying Transform Learning with Convergence Guarantees and Applications”, International Journal of Computer Vision, 114, 137–167 (2015).
  * [46] L. Rebollo-Neira, “Cooperative greedy pursuit strategies for sparse signal representation by partitioning”, Signal Processing, 125, 365–375 (2016).
  * [47] Y.C. Pati, R. Rezaiifar, and P.S. Krishnaprasad,  “Orthogonal matching pursuit: recursive function approximation with applications to wavelet decomposition,”  Proc. of the 27th ACSSC,1,40–44 (1993).
  * [48] L. Rebollo-Neira, J. Bowley, A. Constantinides and A. Plastino,  “Self contained encrypted image folding”,  Physica A, 391, 5858–5870 (2012).
  * [49] L. Rebollo-Neira and D. Lowe, “Optimized orthogonal matching pursuit approach”, IEEE Signal Process. Letters, 9, 137–140 (2002).
  * [50] L. Rebollo-Neira, M. Rozložník, and P. Sasmal, “Analysis of the Self Projected Matching Pursuit Algorithm”,  Journal of The Franklin Institute 8980–8994 (2020).
  * [51] K. Skretting, https://www.mathworks.com/matlabcentral/fileexchange/
2818-huffman-coding-and-arithmetic-codin (Last access May 2021).

  * [52] http://www.nonlinear-approx.info/examples/node015.html (Last access May 2020).
  * [53] https://www.hlevkin.com/06testimages.htm (Last access May 2021).
  * [54] https://imagecompression.info/test_images (Last access May 2021).
  * [55] https://developers.google.com/speed/webp (Last access May 2021).


J

[a]Phani<EMAIL_ADDRESS>if different
from Xiong Nguyen Collins Magazzeni Wilkinson

[a]Department of Materials, University of Oxford, Parks Road, Oxford – OX1
3PH, United Kingdom [b]Safran SA, Safran Tech, Department of Materials and
Process, Rue des Jeunes Bois – Châteaufort, 78772 Many-les-Hameaux, France
[c]School of Metallurgy and Materials, University of Birmingham, Edgbaston,
Birmingham – B15 2TT, United Kingdom

# Weighted ellipse fitting routine for spotty or incomplete Debye-Scherrer
rings on a 2D detector

Yi  Chi-Toan  David M  Christopher M  Angus J

###### Abstract

We introduce a weighted ellipse fitting routine to measure Debye Scherrer
rings acquired on 2D area detectors and demonstrate its use in strain
determination. The method is relatively robust against incomplete rings due to
low number of grains in the diffraction volume (spotty rings), or strong
texture (intensity depletion in some azimuths). The method works by applying
an annular mask around each diffraction ring and fitting an ellipse, using all
pixel positions and their diffracted intensity as weights in the minimisation.
We compare this method to the more popular cake integration method, and show
that the weighted ellipse method works when the cake integration method fails
or works poorly. The lattice strain sensitivity from spotty diffraction rings
is in the order or $2\times 10^{-5}$ or better. The algorithm has been made
available for public use and works with 2D diffraction patterns acquired in a
laboratory scale XRD equipment, TEM or a synchrotron.

An algorithm for a weighted ellipse fitting routine to measure D-S diffraction
ring diameters on 2D detectors in electron or X-ray diffraction.

##  1 Introduction

Microscale deformation is a governing factor for microstructural phenomena
like strain accumulation, dislocation motion, crack initiation amongst others
that determine how a material behaves under load[Fitzpatrick_2003].
Undemanding laboratory based X-ray diffraction methods for residual
strain/stress measurement lack the ability to provide information about the
spatial distribution/concentration or have quite low sufficient spatial
resolution [WithersBhadeshia]. The most popular techniques at a laboratory
scale that provide quite useful information about strain/stress hotspots seem
to be Digital Image correlation [MCCORMICK201052] or diffraction based
electron microscopy methods [WILKINSON2006307, WILKINSON2012366, BECHE201310].
However, these techniques are limited mostly to surface measurement, and DIC
captures strain changes under mechanical/thermal loading rather than initial
residual stress levels. High energy synchrotron sources and neutron sources
allow these measurements for samples with a much greater thickness, and permit
acquisition in a transmission geometry. The nearly parallel beam in a
synchrotron also removes the need for a focussing geometry setup and reduced
low angle peak broadening effects [Willmott2011].

The framework for stress and strain measurement using synchrotron based X-ray
techniques has been quite popular for several decades now [Willmott2011,
CullityB.D2014Eoxd]. Incoming X-ray beam entering a material diffracts when
the Bragg condition is satisfied. A 2D or area detector placed orthogonal to
the incoming X-ray beam can be used to record the conic projections (Debye
Scherrer rings) as whole or portions of rings, depending on the detector size,
distance to sample and detector lateral position. The analysis of X-ray
intensity of these rings and positions can then be used to determine the
lattice plane spacing within the material. With a knowledge of the reference
or initial lattice spacing, these spacing measurements can be converted to
lattice strain (and through elastic constants, Stress)[ERINOSHO20161,
WARWICK20126720]. A geometrically similar situation occurs in the transmission
electron microscope (TEM) in the collection of selected area electron
diffraction (SAED) patterns.

As Erinosho _et al._ [ERINOSHO20161] describe in their paper, good quality and
representative data can be obtained only when a statistically large number of
grains are sampled in the diffracting volume. A statistically large number of
grains is satisfied when the individual diffraction spots cannot be
distinguished from respective grains and the lattice strain measurement from
many grains is representative of the bulk material response. Lattice planes
that do not satisfy the Bragg condition remain undetected, leading to a loss
of information. Textured materials have an even fewer diffracting grains in
some orientations, which often makes the lower multiplicity reflections
unsuitable for use [ERINOSHO20161, WARWICK20126720]. Thus, spotty diffraction
rings arise when the size distribution of diffracting grains is a large
fraction of the incident beam diameter and leads to an orientation
distribution lean/skewed. The most common approach to analysing 2D XRD data
seems to be adapting a powerful ‘cake remapping’ approach that was introduced
by researchers at European Synchrotron (ESRF), built into ‘FIT2D’
[Hammersley1996], and then re-used in many other software packages that are
commonly used to handle and analyse data acquired at a synchrotron source like
DAWN [Filik:vg5068], DIOPTAS [Prescher2015] and pyFAI [Kieffer:gy5006]. This
method involves integrating each diffraction pattern along azimutal sectors of
defined angular spacing so as to convert them into 1D line profiles of
Intensity vs ring radius (or diffraction angle, 2$\theta$) over that azimuthal
sector. This integration process requires that the centre of the diffraction
ring has already been established with good accuracy. An example 2D
diffraction pattern and the resulting 1D intensity profiles generated in two
sectors is shown in Figure 1 (a) and (b).

This approach seems often works reasonably well and the post processing
involves measuring shifts in peak positions for all the reflections, which can
then be converted to lattice strain with respect to a reference [Willmott2011,
CullityB.D2014Eoxd]. The peak position is usually established by identifying
an interval for each lattice plane and fitting a curve (usually a pseudo-Voigt
function) to determine the peak [ERINOSHO20161]. However, for samples with a
strong texture or relatively large grain sizes, some of these sectors may not
have any projected intensity of observable reflections or too much scatter in
the spots that do not yield a workable 1D profile. This effect of texture is
shown in Figure 1 (c) and (d). This peak corresponds to the 0002 reflection or
the 2nd ring (as measured from the centre) in Figure 1 (a). The presence of
texture in the 0002 planes does not produce reasonable intensity along the
Y-direction and hence results in a poor peak fitting routine, but a good fit
along the X-direction. This is also reflected in the intensity of the peak
being much lower along the Y direction.

As _Bronfenbrenner_ [Bronfenbrenner] shows, any three-dimensional strain
manifests as an ellipsoid in the reciprocal space and the projection of this
‘strain ellipsoid’ on a orthogonal 2D area detector takes the shape of an
ellipse. Determining the deviation of shape from a circle to an ellipse can be
used to extract the change in lattice spacing from a strain-free lattice. This
change in the shape of this distorted ellipse with respect to an initial
reference can be used to calculate the variation in lattice spacing and
lattice strain. In this paper, we explore the possibility, sensitivity and
accuracy of fitting an ellipse to these diffraction rings in situations where
the ‘cake remapping’ approach may not be appropriate, without a compromise on
the data obtained.

There has been some interest in fitting Debye-Scherrer rings directly in the
2D data as ellipses rather than manipulating the data to 1D intensity
profiles. Most of these cases have focussed on application to detector
calibration rather than strain determination. Hart _et al._ [Hart2013] and
Borbely _et al._ [Borbely2014] have developed hybrid 1D-2D procedures in which
many 1D profiles spaced around a diffraction ring are used to fit peak
positions which are then used as input data for a 2D analysis to fit an
ellipse. Such procedures require assumptions to be made about the geometry so
that the 1D intensity profiles are extracted on approximately radial paths.
Hart and Drakopoulos [hart2013weighted] suggest that $\sim$700 profiles be
used for continuous rings, but this has to be increased to $\sim$3000 for
spotty rings so that sufficient spots are intercepted. Weighting the
importance of the peak positions using their intensity then emphasises the
stronger parts of the ring in the 2D ellipse fitting. We believe the
approaches in [Hart2013, hart2013weighted] have been incorporated in the
ellipse fitting part of the DAWN [Filik:vg5068] analysis software package.

A more direct 2D approach was outlined by Hanan _et al._ [hanan2004new], in
which positions of all pixels associated with a diffraction ring were used
directly in a least squares fit to an ellipse. The pixels associated with a
particular ring were established by (i) a mask of user defined width around
the estimated ring (to avoid overlap with other nearby rings), and (ii) an
intensity threshold (rejecting low intensely pixels with the masked region).
This was performed for a series of different intensity thresholds and the
intensity weighted averaged used as the best fit for that ring.

The method introduced here is closer to that by Hanan _et al._ [hanan2004new],
in that it uses a mask to define pixels associated with individual ellipses
but then uses a single fit to an ellipse with pixels weighted by their
normalised diffraction intensity using an algorithm like that used by Hart and
Drakopoulos [hart2013weighted], in turn a modification of Gal’s implementation
of a least squares fit to an ellipse [OhadGal]. This method may well be
applied for other techniques that involve ellipse fitting to diffraction rings
in lab-scale XRD, $\mu$-XRD or TEM based measurement of lattice parameters.
The superiority of this method in comparison to the 1D integration routine is
shown here for experimentally obtained patterns with incomplete and spotty
rings.

A Matlab routine to perform this analysis has been made available for public
use on GitHub [PhaniKaramchedGitHub:Online]. A similar approach in python has
also been made available for open source users in the directory. No distortion
correction has been applied to the images in this algorithm. The script takes
as input the raw diffraction images as #.tif and needs an approximate
knowledge of the pattern centre (to apply a mask around individual rings). As
suggested by other researchers, a self-convolution of the image can be applied
to locate this approximate pattern centre [hart2013weighted,
Cervellino:aj5105]. The algorithm estimates a more precise centre of the
individual rings, along with the major and minor axes lengths and rotation of
the major axis of the weighted ellipse fit. The precise lattice dimensions
along the X and Y directions can then be calculated.

##  2 Experimental method

Diffraction patterns were obtained as a part of a larger experiment at the I12
beamline at Diamond Light source [Drakopoulos:ie5138]. We first demonstrate
the effectiveness of the weighted ellipse fitting routine using a series of
ten powder diffraction patterns from a calibration standard CeO2 (NIST,
SRM674b) sample. During the time of recording these patterns, no experimental
variation in the sample or geometry is applied. We generate spotty and noisy
patterns from this sample and apply the ellipse fitting routine. Noise and
spotty pattern have been generated in three ways, (a) by adding white gaussian
noise, (b) removing a fixed percentage of points in the pattern at random and
(c) removing data along the loading and transverse directions. We compare the
cake integration method with a definition of $\pm 5^{\circ}$ sectors and
fitting a pseudo-Voigt function to identify the peak position as this seems to
be the most popular current approach in the literature [ERINOSHO20161,
COLLINS201546] against the weighted ellipse fit method we introduce in this
paper. Upon removing data in sectors close to the loading or transverse axis
(to simulate the effect of textured samples), the weighted ellipse fit is
still able to measure the lattice parameter along the X and Y directions, but
the cake integration method fails to work due to absence of intensity data
along that direction.

As a part of the larger experiment, various Ti samples under were deformed
under tensile load in-situ using an Instron Electro-Thermal-Mechanical Tester
(ETMT) system (described elsewhere in [XIONG2020561, XIONG2021116937]). A
schematic of the setup is shown in Figure 2(a). A typical spotty diffraction
pattern is also shown in Figure 2(b). Such rings as already discussed are a
result of a smaller number of grains/orientations being sampled in the
diffraction volume, either due to a larger grain size or texture. An ETMT was
placed orthogonal to the incident X-ray beam to perform uniaxial tensile
experiments. Tensile load was applied as shown in Figure 2(a) along the
Y-direction and the Debye-Scherrer rings recorded on a Thales Pixium RF4343 2D
area detector, which has a $2880\times 2881$ X-ray sensitive array of pixels,
of size 148 $\mu$m [Drakopoulos:ie5138]. The image acquisition time per
pattern was fixed at 1 s and the incident beam size was set to $1\times 1\
mm^{2}$ on a sample of thickness 1 mm. This weighted ellipse approach has
already been used to evaluate the lattice strains for several other
reflections in four different cp-titanium grades, the results of which are
available elsewhere [XIONG2020561, XIONG2021116937].

Two sufficiently large grained titanium samples are also compared in this
paper. The first dataset is from a cp-titanium (grade-2) sample that has a
grainsize of $\sim 40\mu m$ and spotty rings (ring-2 and ring-4). The ring-2
and ring-4 correspond to the 0002 and 01-12 reflections for the cp-titanium.
The second dataset discussed is from a tensile test of the same sample of cp-
titanium (grade-2) and the strain measurement from the relaxation part of the
test is demonstrated, which normally is a result of cold-creep in these
alloys.

The proposed method to fit a weighted ellipse to the diffraction rings uses an
approximate ring centre to apply an adjustable annular mask around each ring,
and then uses the intensity values as weights to fit an ellipse. This method
is almost directly based on the Fitzgibbon approach [Fitzgibbon] and is
mathematically simple, easy to apply and computationally very efficient.
Analysis of 10 rings in 16-bit depth images on a regular desktop Intel i5
($8^{th}$ gen) processor and 8gb RAM tends to take of the order of 1 s per
image, without parallel computing. The weighted ellipse fitting routine
reports an accurate ellipse centre (independent of the input ring centre for
the mask), major and minor axis lengths and the semi major axis angle.

###  2.1 Ellipse fitting routine

In the cartesian system, the equation for a general conic section follows the
form:

$F(p,x)=p\cdot x={ax^{2}+bxy+cy^{2}+dx+ey+f=0}$ (1)

where $p=[a\ b\ c\ d\ e\ f]^{T}$ and $x=[x^{2}\ xy\ y^{2}\ x\ y\ 1]^{T}$. The
fitting of a general conic can be approached by minimizing the sum of squared
algebraic distances of the curve to the points. This methods can be better
used by applying appropriate constraints to the conic section and fit an
ellipse [Fitzgibbon]. This implementation has already been presented as a
remarkedly compact MATLAB implementation by Fitzgibbon _et al._ [Fitzgibbon],
and we modify these equations to make a weighted fit.

The set of points associated with the ring $(x_{i},y_{i})$ are represented as
a design matrix defined by $D$ and converted to a weighted matrix $D1$ as
shown in equations 2.

$\begin{gathered}D=[x.^{\wedge}\ 2,\ x.*y,\ y.^{\wedge}\ 2,\ x,\ y,\ 1]\\\
w=repmat(intensities,[1,6])\\\ D1=D.*w\end{gathered}$ (2)

The original set of points are assigned to a design matrix $D$ and converted
to a weighted design matrix $D1$, after multiplying with normalized weights as
a vector $w$, obtained from the intensities at the positions around the
diffraction ring. Mathematical solution follows in line with the Fitzgibbon
_et al._ approach [Fitzgibbon] as follows:

$\begin{gathered}\%Construction\ of\ a\ scatter\ matrix\ S\\\
S=D1^{\prime}*D1\\\ \%Build\ a\ 6\times 6\ constraint\ matrix\ C\\\ C\ =\
zeros(6);\ C(1,3)\ =\ 2;\ C(3,1)\ =\ 2;\ C(2,2)\ =\ -1;\\\ \%Solve\ the\
Eigen\ system\\\ [E,V]\ =\ eig(S\backslash C);\\\ \%Find\ the\ positive\
eigenvalue\\\ [PosR,PosC]\ =\ max(max(abs(V)));\\\ \%Extract\ eigenvector\
corresponding\ to\ positive\ eigenvalue\\\ a=E(:,PosC);\\\ \end{gathered}$ (3)

The original MATLAB based implementation has been preserved in this paper to
highlight the difference in the design matrix used from Fitzgibbon’s approach
[Fitzgibbon]. The logical procedure is described as a comment for each
individual command in equation 3 as a comment (indicated by a % symbol at the
beginning of each line). The method uses a generalised eigen-system and
introducing constraints to specifically fit an ellipse. Amongst the six
eigenvalue - eigenvector pairs for this constrained solution, exactly one
positive eigenvalue exists and this unique solution defines the desired
elliptical fit. More mathematically trained users are directed to [Fitzgibbon]
for further reading.

##  3 Results and Discussion

Figure 3 demonstrates that the weighted ellipse fitting routines works even in
the presence of noisy patterns from a powdered calibration standard sample of
CeO2. For raw patterns obtained, the standard deviation in ring-1
(corresponding to the 111 reflection) semi-length measurement is 0.0013 (in
pixel units) via the weighted ellipse fit method and 0.0062 via the cake
integration approach. In the presence of white noise, the standard deviation
in the peak-1 semi-length measured (in pixels) is 0.3274 using the cake
integration method and a pseudo-Voigt fit. The standard deviation measured by
the weighted ellipse fitting routine for the corresponding ring-1 is 0.1784.
We compare a series of datasets, where randomly a percentage of pixels in the
patterns have been brought down to zero. Both the techniques compare extremely
well as can be seen in Figure 4 even when 50% of the points have been zeroed.
This demonstrates that the weighted fitting routine is not only self-
consistent, but also consistent with the cake-integration approach. When the
intensity along the viscinity of the horizontal and vertical axes are set of
zero, the cake integration fails to work, due to the absence of any intensity
usable for a pseudo-Voigt fit, but the weighted ellipse fit routine still
manages to index as shown in Figure 3(d). The standard deviation measured in
the ring-1 semi-length by the ellipse fit is 0.0017 (in pixel units) and
closely compares to measurement from the raw patterns.

The two methods are next compared on diffraction patterns from polycrystal
(non-powdered) cp-titanium samples (grade-2). Figure 2(b) shows the quality of
the weighted ellipse fit for spotty diffraction rings and with poor definition
of high intensity spots. As with the CeO2 sample, we measure the deviation on
the cp-titanium (grade-2) sample from a recorded set of 30 patterns. During
this time sequence of patterns obtained, no experimental parameters in the
sample or geometry were applied. We choose to make this measurement for poorly
defined set of intensity on the reciprocal lattice spacing of the 0002 plane,
or ring-2 with sparsely populated high intensity spots. This particular
orientation in titanium alloys plays a significant role in fatigue and creep
deformation, especially when oriented along the loading direction
[Dunne2008_facet, XIONG2020561]. The ellipse fitting routine used to estimate
the variation in the ellipse semi-lengths, shows the sensitivity of this
measurement from the spotty ring (ring-2) from the 0002 plane in Figure 5.

The semi-length for the elliptical ring measurement sensitivity of the
weighted ellipse technique seems to be better than 0.01 of a pixel. The
standard deviation measured in the movement of the centre of the ellipse
translated to better than 0.005 in X and 0.007 of a pixel in Y. The standard
deviation of the measured lengths along the X and Y axes of the spotty ring-2
of the cp-titanium sample is $0.0065$ and $0.0075$ (in pixel units)
respectively. This sensitivity in the ellipse determination translates to a
strain sensitivity ($\Delta d/d_{0}$) better than $2\times 10^{-5}$ in these
cp-titanium samples, along the direction (Y) with the poorest diffraction
information.

Upon using cake-integration, the 1D intensity profiles generated from the same
pattern into 36 10° wide sectors have clearly defined peaks for only 6 of the
36 profiles. Most data points do not have sufficient information to define a
peak and follow peak shift measurements. When the data is divided into smaller
sectors, the peak definition along the X and Y directions is completely
missing or has large errors for lattice strain estimation. This missing
definition for the 0002 reflection (or ring-2) is clearly visible in Figure
2(b). With a well-defined ellipse, the strains along any direction orthogonal
to the beam direction can be estimated, but without an appropriate peak
location, the cake integration method is unusable. The issues with the peak-
fitting routine from 1-D profiles for a similar sample (having a better
intensity definition) have already been shown in Figure 1(c) and (d).

For a second comparison, we use data from a sample of cp-titanium (grade-2),
which has been deformed in an ETMT in the strain relaxed regime. In this
comparison we also remove the intensity points along the loading axis of this
sample and show that the weighted fit algorithm tracks the ellipse
appropriately where the cake integration method is known to fail.

Figure 6(a) shows the loading, relaxation and unloading regime for the cp-
titanium sample that was subjected to in the ETMT. After loading slightly
beyond the yield point, the stroke arm was held constantly for 5 min. This
load (and hence strain) drop during the relaxation period is commonly referred
to as cold-creep in these alloys and our primary interest. During this test,
the 2D diffraction Debye-Scherrer rings were recorded and lattice spacing
measured to estimate lattice strain. We now compare the cake integration vs
ellipse fitting routines for ring-2 and also show that despite masking the
intensity points along the loading direction, the weighted ellipse fitting
routine is able to measure the spacing (and changes) along the Y-direction
(loading direction).

It is expected that the lattice spacing and strains along the X and
Y-directions also follow a trend (either tensile or compressive depending on
the Poisson’s compression and the influence of other crystal planes) as the
relaxation shown in Figure 6(a). This relaxation is tracked after a time
period of 110 seconds after the start of the test until complete unloading. A
comparison is shown in Figure 7.

As discussed already, we make this comparison for the measurements of the
reciprocal lattice spacing of the 0002 plane. It is quite clear that in the
X-direction, where this particular plane has a poor definition, the
sensitivity of the weighted ellipse method is superior to the cake-integration
approach. This also holds true when the intensity points are masked along the
Y-direction and here the cake-integration approach fails completely.
Furthermore, the cake integration method is highly dependent on the
calibration used to perform this analysis, that identifies the Beam centre.
The weighted ellipse fit is independent of this and a moderately accurate beam
centre definition is required only to apply a mask. The centre is then found
more precisely from the ellipse fitting. This value of the measured peak
position or ring semi-lengths (in pixels) can be converted to lattice spacing
(and strain) upon further knowledge of geometry of the setup and is described
in quite a bit of detail in literature [Fitzpatrick_2003, CullityB.D2014Eoxd,
COLLINS201546].

##  4 Conclusions

We have thus demonstrated that a weighted ellipse fitting of diffraction rings
is a superior, computationally efficient and a more elegant method for
identifying/measuring lattice spacing and measuring strain from 2-D
diffraction patterns, when compared to the popular cake-integration method.
The approach presented here is a combination of several other already existing
methods for fitting ellipses and works on raw patterns without a critical need
for calibration routines. In cases where the diffraction rings are spotty (due
of a smaller number of diffracting crystals in the volume or due to the use of
textured samples), the cake integration method may not always be suitable for
processing by this method. The weighted ellipse fitting routine presented here
and made available for use [PhaniKaramchedGitHub:Online] is more accurate and
usually faster.

The weighted ellipse method to measure the ring semi-lengths does not need a
calibrated pattern center as is required in the case of the cake integration
method. The method only relies on the knowledge of an approximate pattern
center to apply an annular ring mask. An accurate pattern center is measured
by the fitting routine. Thus, this method can also be applied on a sample
regularly used for calibration or for precisely measuring lattice spacing in
other methods such as laboratory based XRD or TEM based diffraction. The
weighted-ellipse fitting routine also does not assimilate nearby data in
comparison to the cake-fitting method that needs a defined azimuthal
averaging. The sensitivity in ring semi-lengths measured by the weighted
ellipse method is better than $0.01$ (in pixels units) for spotty rings on a
cp-titanium sample and the strain sensitivity $2\times 10^{-5}$ or better. For
a standard powder sample of CeO2, the standard deviation in the measured value
of ring semi-length corresponds to a value of 0.0015 (in pixel units) and in
terms of measured d-spacing, this value translates to better than $1.5\times
10^{-5}$.

##  5 Acknowledgements

The authors acknowledge funding from the EPSRC through the HexMat programme
grant (EP/K034332/1) and the Diamond Light Source for beam time under
experiment EE17222. We would like to thank Dr.Thomas Connolley, Dr.Robert
Atwood and Dr.Stefan Michalik for their friendly and patient help at the
beamline I12 and Dr.Nicolo Grilli for the helpful discussion.

##  6 Author Contributions

This work formed part of a larger testing campaign proposed by AJW and carried
out at Diamond Light Source by YX, PK, CTN, DMC, CMM and AJW. PK oversaw
preparation of test pieces and preliminary characterisation undertaken by YX,
CMM and PK. CTN and PK setup and guided use of the ETMT system. PK, YX, DMC
and AJW developed data analysis procedures. PK wrote the script for pseudo-
voigt fitting and ellipse fitting routines. PK drafted the initial manuscript,
with all authors contributing to and agreeing the final version.

##  7 List of Figures

Figure 1: (a) 2D diffraction pattern acquired on an area detector (plotted as
a logarithmic intensity), depicting sectors $\pm 5^{\circ}$ around an
azimuthal angle $\psi=0^{\circ}$(i.e. parallel to X in blue) and
$\psi=90^{\circ}$(i.e. parallel to Y in red) from a cp-titanium sample
(grade-4), (b) azimuthally integrated in the 2 sectors to show 1D profile of
Intensity vs ring radius along X and Y. (c) & (d) show Pseudo-Voigt peak
fitting for 0002 plane along $\psi=0^{\circ}$(i.e. parallel to X) and
$\psi=90^{\circ}$(i.e. parallel to X) Figure 2: (a) Experimental setup at I12
Diamond Light Source, showing a diffraction pattern from a cp-titanium
(grade-4) sample on the Pixium detector (b) Typical spotty diffraction rings
(zoomed-in) from a sample of cp-titanium (grade-2) and the working weighted
ellipse fit routine. Figure 3:  5 rings indexed in (a) 2D diffraction pattern
from a CeO2 (NIST, SRM674b) sample (b) the presence of white gaussian noise
(c) 50% of the pixels reduced to zero intensity (d) points removed along
horizontal and vertical directions. Figure 4: Ellipse major axis semi-
length(in blue) or peak positions along Y-direction (in red), measured in
pixels via the weighted ellipse fitted method and the cake integration
approach for a fixed percentage of randomly zeroed pixels in the image for the
111 reflection in the CeO2 sample . The error bars depict the standard
deviation in 10 patterns for the first ring/peak. The values on the horizontal
axis have been deliberately offset for clarity. Figure 5: Measured ellipse
semi-lengths (ring-2) from 30 patterns in the 0002 plane of cp-titanium
(grade-2). Figure 6: (a) Load (N) vs Time (s) recorded from the ETMT. (b) Mask
applied to remove data from the loading direction and weighted ellipse fit.
Figure 7: Demonstrating the comparison of the two methods (cake integration vs
weighted ellipse fit method) of the reciprocal lattice spacing (in pixel
units) of the 0002 plane along (a) X-direction (b) Y-direction or the loading
direction.

[mybibfile]