# A Neural Template Matching Method to Detect Knee Joint Areas

Juha Tiirola Email<EMAIL_ADDRESS>Research Unit of Mathematical
Sciences, University of Oulu, Finland

###### Abstract

In this paper, new methods are considered to detect knee joint areas in
bilateral PA fixed flexion knee X-ray images. The methods are of template
matching type where the distance criterion is based on the negative normalized
cross-correlation. The manual annotations are made on only one side of a
single bilateral image when the templates are selected. The best matching
patch search is formulated as an unconstrained continuous domain minimization
problem. For the minimization problem different optimization methods are
considered. The main method of the paper is a trainable optimizer where the
method is taught to take zoomed and possibly rotated patches from its input
images which look like the template. In the experiments, we compare the
minimum values found by different optimization methods. We also look at some
test images to examine the correspondence between the minimum value and how
well the knee area is localized. It seems that making annotations only to a
single image enables to detect knee joint areas quite precisely.

Keywords: Automatic Detection, Convolutional Neural Networks, Spatial
Transformer

##  1 Introduction

Knee osteoarthritis (OA) refers to a disease where the knee articular
cartilage is degreded. The major features to assess the disease are joint
space narrowing, formation of osteophytes and sclerosis. The classical X-ray
imaging based grading system to knee osteoarthritis is the KL grading system,
proposed by Kellgren and Lawrence in [11] which gives an integer score from
$0$ to $4$ to each knee image. In this system, the higher the score, the more
severe the disease. The KL grading system is based on the features found in
the knee joint area. Thus if one wants to exclude the influence of the
background on the KL classification, a natural first step is to localize the
knee joint areas from the bilateral image. In [4], an overview is given about
the disease and the problems of localizing the joint area and the automatic
assesment of OA severity. In [15], template matching is used to automatically
detect knee joint areas. They use the Euclidean distance as the closeness
criterion and $20$ template images. In [4], the template matching method of
[15] is also discussed and they state that the detections by the method may
depend much on the choice of the template images. Other approaches also exist,
in [13] the whole bilateral X-ray image was used as an input for learning to
measure the severity of osteoarthritis. In [12], the fast computation of
normalized cross-correlation is considered. If template matching is formulated
by comparing the template to each possible patch in the input image in a
sliding window way, it is not so easy to incorporate other prior information,
for instance that the left and right knee joint areas are nearly of equal
size, and still get fast computation methods. In contrast, if the template
matching is formulated as a continuous domain minimization problem, it is easy
to add regularizing terms to the template matching energy and one can use the
methods of gradient based optimization instead of brute-force based approaches
to find the optimal patch.

In [5], Sobel horizontal image gradients are used as the input features for
detecting knee joints. In training, the knee joint areas are manually selected
in each training image. This leads to positive and negative samples depending
on whether the center of the labeled Region of Interest (ROI) is inside the
sample. Using a linear SVM classifier they learn to classify negative and
positive samples. In the test phase, the maximal score is taken as the center
of the detected ROI and a fixed size neighbourhood is taken as the detected
ROI. In [3], binary images with masks specifying the manually annotated knee
joint areas are used as the training labels. Then a fully convolutional
network is learnt to map the input image to a soft pixel-based classification
map. In the test phase, the estimated ROIs are obtained from the soft output
by contour detection. In [16], a two-step method is proposed to the knee joint
area detection where first knee joint area proposals are generated based on
the limb anatomy. Then these proposals are scored based on their Histogram of
oriented gradients features. Training of the scoring step is based on the SVM
classifier with positive and negative samples depending on the amount of
overlap with the manually selected bounding boxes. In [6] and [17], YOLOv2
method of [14] is finetuned in a manually annotated data set to find knee
joint areas. In [8], Yolo v3-tiny of [1] is trained to detect knee joint
areas.

In [10], the spatial transformer module was introduced to warp feature maps in
a differentiable way. The module can be inserted into neural network
architectures without a need for extra labeling. In the method, given a
parametric family of transformations, the localization net analyses its input
feature map and outputs a parameter vector. The transformation corresponding
to the parameter vector is then used to warp the input feature map. If the set
of affine transformations is used as the parametric family, the spatial
transformer module enables to take zoomed and rotated patches from its input
feature map in a differentiable way.

Gradient descent is the classical optimization approach also in the training
of neural networks. For instance in [2], also the design of optimization
algorithms is cast as a learning problem. They proposed to learn a neural
network optimizer when different but related minimization problems are
available as the training data, and they show that the learnt optimizer can
generalize to unseen, but nevertheless training-related minimization problems.
One example problem considered in [2] is the style transfer problem [7]. In
the problem, given a style patch, an input image is transformed in a structure
preserving way into a same style image as the given style patch. The same
style image is obtained as the solution of a minimization problem. If the
input image is changed, a different but related minimization problem is
obtained.

In this paper, we consider a neural network based method to the knee joint
area localization problem. As the training data, we use $500$ bilateral knee
X-ray images and manual annotations are made only on a single X-ray image. The
annotation includes selecting a bounding box which forms the template image.
In addition, we select two smaller boxes inside the template image at the left
and right sides of the knee joint borders. We show the manual annotation done
in Fig. 2. Our method uses the sampling mechanism of the spatial transformer
to obtain a patch from the input image in a differentiable way. The spatial
transformer module warps the input image such that the warped output looks
like the template patch. The corresponding transformation parameter vector is
found by solving a minimization problem where the energy expression to be
minimized is based on the negative, normalized cross-correlation between the
extracted patches and the template image plus some regularization. The
minimization problem has $8$ parameters. Since for each training image there
is a minimization problem, similar to [2], we learn the optimizer. The
architecture of the optimizer is a Resnet18 backbone followed by two
additional 1-1 convolutional layers followed by a regression network which
consists of a three layer fully connected neural network.

For comparison, we also consider to solve each minimization problem
separately. The minimization problem is highly non-convex so the outputs of
the gradient based optimization algorithms depend heavily on the
initialization. To minimize the energy for a given bilateral image, we run a
gradient based optimization for many different initializations. The
initializations are based on a grid. First, the grid has a few scales. For a
given scale, many translations are tested for the initial center point of the
candidate optimal patch. When the gradient based optimization is run, we
record all loss values. Finally, we select the parameter corresponding to the
smallest loss encountered. This method is time consuming since for a given
bilateral X-ray image, depending on the grid there may be thousands of
optimization problems to be solved.

##  2 Pre-processing

The input images are bilateral X-ray knee images as $u$ in Fig. 1. The goal is
to detect the knee joint areas from the input image.

Figure 1: The deterministic splitting and random augmentation steps. The image
$u$ is split first into somewhat overlapping halves $u_{l}$ and $u_{r}$. Then
zeros are added to make the ratio of the height over width of each half a
constant. This deterministic step is followed by random translation, the
horizontal flip of the left half and resizing both halves to $800\times 500$
pixels.

After the preprocessing we make the problem simplifies somewhat in the sense
that there is only one knee joint area from each half to be found. An overview
of the split process and augmentation is shown in Fig. 1.

We use a down-sampled version of $u$ to roughly divide $u$ into vertically
symmetric halves. If $h_{0}$ is the horizontal coordinate of the center point
of the downsampled image $\hat{u}$, then for each
$a\in\\{h_{0}-0.125width(\hat{u}),h_{0}+0.125width(\hat{u})\\}$ we do the
following: let $w:=\min\\{a,\text{width}(\hat{u})-a\\}$. Compute the average
Euclidean distance between $\hat{u}[:,a-w:a]$ and
$\text{horizontal}\\_\text{flip}(\hat{u}[:,a:a+w])$. Select the $a$ giving the
smallest average distance as the location where the split is done. Then widen
the non-overlapping halves by the factor $1.1$ to get the slightly overlapping
halves $u_{l}$ and $u_{r}$. Let $w_{l}$ and $w_{r}$ be the width of $u_{l}$
and $u_{r}$, respectively. If $w_{l}>w_{r}$, then we widen $u_{r}$ by adding
zeros to its both sides by equal amounts so that the augmented $u_{r}$ has
width $w_{l}$. Otherwise, we widen $u_{l}$. As a result both $u_{l}$ and
$u_{r}$ are of the same size $h\times w$. If $\frac{h}{w}>1.6$, we widen both
$u_{l}$ and $u_{r}$ by adding zeros to both sides such that $\frac{h}{w}$
becomes $1.6$. If $\frac{h}{w}<1.6$, we add zeros both on top and bottom
parts. The whole process has been so far deterministic. Next we add a small
random translation in the training step. After this we resize both halves to
size $800\times 500$ pixels and we flip the left part horizontally. We get
images which we call $u^{left}$ and $u^{right}$.

##  3 Parametrization

The extraction of the knee joint areas from $u^{left}$ and $u^{right}$ is
based on the use of the sampling mechanism of the spatial transformer. Given a
parametric family of transformations the spatial transformer has a
localization net which takes a feature map as an input and outputs
transformation parameters. The transformation corresponding to the parameters
is used to produce sampling points in the input feature map domain. Then a
sampling mechanism is used to produce a warped feature map by interpolating
the input feature map at the sampling points. Here we would like to place a
possibly rotated rectangle inside each input image such that the rectangle
defines the joint area. A parameter vector specifies the scale, location and
rotation of the rectangle. Since the area must be inside the input image, the
parametrization is constrained. The parameter vector $\boldsymbol{\theta}^{s}$
determines the corresponding transformation $A(\boldsymbol{\theta}^{s})$ which
is used to warp a regular grid to obtain the sampling points. Each half of
Fig. 1 has its own transformation parameter vector. The regular grid of the
input $u^{s}$ lies in $[-1,1]^{2}$ and the sampling points lie inside the
warped square $\\{\,A(\boldsymbol{\theta}^{s})\begin{pmatrix}\boldsymbol{x}\\\
1\end{pmatrix}\,|\,\boldsymbol{x}\in[-1,1]^{2}\,\\}$. We want the warped
square to be inside $[-1,1]^{2}$. We denote by
$u^{s}(\boldsymbol{\theta}^{s})$ the corresponding warped image obtained by
applying the sampler of the spatial transformer to interpolate $u^{s}$ at the
sampling points. The template matching tries to find a suitable
$\boldsymbol{\theta}^{s}$ such that the distance between
$u^{s}(\boldsymbol{\theta}^{s})$ and the template is as small as possible.

Let $f:=\text{height}(\text{template})/\text{width}(\text{template})$. We
assume that the template is such that $f\geq 1$. A warped image has its own
regular normalized coordinates $[-1,1]^{2}$. This square is multiplied
coordinatewise by the scale vector
$(s_{1},s_{2}):=(\boldsymbol{\theta}^{s}_{1},\boldsymbol{\theta}^{s}_{1}/f)$,
possibly rotated slightly by $\boldsymbol{\theta}^{s}_{4}$ radians, translated
by a vector $(\boldsymbol{\theta}^{s}_{2},\boldsymbol{\theta}^{s}_{3})$, and
the resulting transformed coordinate square must lie inside the $[-1,1]^{2}$
coordinate square of the input half image. Thus, $\boldsymbol{\theta}^{s}$ is
constrained and required to be meaningfully interpretable.

Instead of searching the transformation parameters directly as the solution of
a constrained minimization problem, we consider an unconstrained optimization
problem where an arbitrary, unconstrained parameter vector
$\boldsymbol{v}^{s}\in\mathbb{R}^{4}$ is mapped to an interpretable
transformation vector $\boldsymbol{\theta}^{s}\in\mathbb{R}^{4}$.

Let $\boldsymbol{v}^{s}\in\mathbb{R}^{4}$ be arbitrary. First we consider how
the scale$\in(0,1)$ is computed from $\boldsymbol{v}^{s}_{1}$. In Fig. 1 we
see a typical input halves $u^{left}$ and $u^{right}$ from which we try to
extract patches which look like the template patch $\boldsymbol{T}$ of Fig. 2.
The sizes of the extracted patches cannot be arbitrarily small, instead the
width of the patch is typically like a third of the width of the halves. Thus
we map the domain $\mathbb{R}$ of $\boldsymbol{v}^{s}_{1}$ into a meaningful
interval subset of $(0,1)$. We set $\alpha_{0}=0.15$, $\beta_{0}=0.8$, and

$\boldsymbol{\theta}^{s}_{1}:=\alpha_{0}+\beta_{0}\left(\frac{1}{2}\left(1+\tanh(\boldsymbol{v}^{s}_{1})\right)\right).$
(1)

We see that
$\boldsymbol{\theta}^{s}_{1}\in[\alpha_{0},\alpha_{0}+\beta]\subset[0,1]$.

Next we come to the choice of the translation parameters. For each scale, the
possible locations for the center of the zoomed patch are constrained such
that the extracted patch lies fully inside the input image. If
$\boldsymbol{\theta}^{s}_{2}$ gives the horizontal center location of the
rectangle of width $2s_{1}$, then we get constraints
$-1\leq\boldsymbol{\theta}^{s}_{2}-s_{1}$ and
$\boldsymbol{\theta}^{s}_{2}+s_{1}\leq 1$, and thus
$-1+s_{1}\leq\boldsymbol{\theta}^{s}_{2}\leq 1-s_{1}$. If we define

$a_{13}:=\frac{1}{2}(1+\tanh(\boldsymbol{v}^{s}_{2}))\text{ and
}a_{23}:=\frac{1}{2}(1+\tanh(\boldsymbol{v}^{s}_{3}))$

and then set

$\begin{split}\boldsymbol{\theta}^{s}_{2}&:=-1+s_{1}+a_{13}\cdot
2(1-s_{1})=(1-s_{1})(-1+2a_{13}),\\\
\boldsymbol{\theta}^{s}_{3}&:=-1+s_{2}+a_{23}\cdot
2(1-s_{2})=(1-s_{2})(-1+2a_{23})\end{split}$ (2)

we get the above constraints satisfied.

In addition, we add a possibility for small rotation. We set

$\boldsymbol{\theta}^{s}_{4}:=0.13\tanh(\boldsymbol{v}^{s}_{4}).$

The constraint for the location was for a non-rotated rectangle so in theory a
part of a rotated rectangle might be outside $[-1,1]^{2}$ but this was
negligible in the experiments.

We have described how to transform an unconstrained
$\boldsymbol{v}^{s}\in\mathbb{R}^{4}$ to an interpretable transformation
vector $\boldsymbol{\theta}^{s}\in\mathbb{R}^{4}$. We denote by

$D:\mathbb{R}^{4}\to\mathbb{R}^{4}$ (3)

the above map $\boldsymbol{v}^{s}\mapsto\boldsymbol{\theta}^{s}$.

If $\boldsymbol{\theta}^{s}=D(\boldsymbol{v}^{s})$, the transformation matrix
corresponding to $\boldsymbol{v}^{s}$ is given by

$A(\boldsymbol{\theta}^{s}):=\begin{bmatrix}\boldsymbol{\theta}^{s}_{1}\cos(\boldsymbol{\theta}^{s}_{4})&-\boldsymbol{\theta}^{s}_{1}\sin(\boldsymbol{\theta}^{s}_{4})&\boldsymbol{\theta}^{s}_{2}\\\
(\boldsymbol{\theta}^{s}_{1}/f)\sin(\boldsymbol{\theta}^{s}_{4})&(\boldsymbol{\theta}^{s}_{1}/f)\cos(\boldsymbol{\theta}^{s}_{4})&\boldsymbol{\theta}^{s}_{3}\end{bmatrix}.$
(4)

When the above matrix is applied to a vector $\boldsymbol{x}\in[-1,1]^{2}$, it
is understood that $\boldsymbol{x}$ is augmented first by adding $1$ as the
third dimension, so the result of applying the matrix to the vector is
$A(\boldsymbol{\theta}^{s})\begin{pmatrix}\boldsymbol{x}\\\
1\end{pmatrix}\in\mathbb{R}^{2}$, see [10].

##  4 Loss function

We manually select a template $\boldsymbol{T}$ and its two sub-templates from
a training image, see Fig. 2. If we have a loss function $L$, which measures
the distance between warped half images and the template, then a way to find
the warping transformations is to solve the unconstrained minimization problem

$\min_{(\boldsymbol{v}^{left},\boldsymbol{v}^{right})\in\mathbb{R}^{8}}L(u^{left}(D(\boldsymbol{v}^{left})),\boldsymbol{T})+L(u^{right}(D(\boldsymbol{v}^{right})),\boldsymbol{T}).$

Next we consider how to choose the loss function.

Figure 2: A template $\boldsymbol{T}$ and its two sub-templates shown by red
and green boxes. The images $u^{s}(\boldsymbol{\theta}^{s})$ at the regular
grid points $\boldsymbol{x}$ are obtained by interpolating $u^{s}$ at the
points $A(\boldsymbol{\theta}^{s})[\boldsymbol{x},1]^{T}$.

Let $\mathbf{u}$ and $\mathbf{v}$ be two images. Let $\tilde{\mathbf{u}}$ and
$\tilde{\mathbf{v}}$ be the same images with their means removed. Define a
global registration cost using the negative normalized cross-correlation as

$c^{gl}(u,v):=1-\frac{\tilde{\mathbf{u}}\cdot\tilde{\mathbf{v}}}{||\tilde{\mathbf{u}}||\,||\tilde{\mathbf{v}}||}.$

In addition to the global window indicated by the blue rectangle at the
boundary in Fig. 2, we also consider costs over two sub-windows whose fixed
locations are indicated by the red and green windows in Fig. 2. We define sub-
costs

$c^{r}(u,v):=1-\frac{\widetilde{\mathbf{u}|_{rw}}\cdot\widetilde{\mathbf{v}|_{rw}}}{||\widetilde{\mathbf{u}|_{rw}}||\,||\widetilde{\mathbf{v}|_{rw}}||},~{}c^{g}(u,v):=1-\frac{\widetilde{\mathbf{u}|_{gw}}\cdot\widetilde{\mathbf{v}|_{gw}}}{||\widetilde{\mathbf{u}|_{gw}}||\,||\widetilde{\mathbf{v}|_{gw}}||}$

where for instance $u|_{rw}$ denotes $u$ restricted to the interior of the red
bounding box. We add the global cost and sub-costs and define

$\tilde{c}(u,v):=\frac{1}{2}\left(c^{gl}(u,v)+\max\\{c^{r}(u,v),c^{g}(u,v)\\}\right).$

Since there are also image negatives in the training data, we take this into
account and define the loss function by

$L(u,v):=\min\\{\tilde{c}(u,v),\tilde{c}(-u,v)\\}.$

As visible also in Fig. 2, the knee joint areas between halves are quite equal
sized and in addition the vertical coordinates of the center of the knee joint
areas are near each other. Since $\boldsymbol{\theta}^{s}_{1}$ gives the size
and $\boldsymbol{\theta}^{s}_{3}$ the vertical coordinate of the center of the
extracted patch, we would like to have
$\boldsymbol{\theta}_{1}^{left}\approx\boldsymbol{\theta}_{1}^{right}$ and
$\boldsymbol{\theta}_{3}^{left}\approx\boldsymbol{\theta}_{3}^{right}.$ If we
have a bilateral image $u$ and
$\boldsymbol{v}:=(\boldsymbol{v}^{left},\boldsymbol{v}^{right})\in\mathbb{R}^{8}$,
we define a regularized cost as

$\begin{split}f^{\boldsymbol{T}}(u,\boldsymbol{v}):=&\mathcal{L}^{left}+\mathcal{L}^{right}+\mathcal{L}^{reg}\\\
:=&L(u^{left}(D(\boldsymbol{v}^{left})),\boldsymbol{T})+L(u^{right}(D(\boldsymbol{v}^{right})),\boldsymbol{T})\\\
+&||D(\boldsymbol{v}^{left})_{1}-D(\boldsymbol{v}^{right})_{1}||_{2}^{2}+||D(\boldsymbol{v}^{left})_{3}-D(\boldsymbol{v}^{right})_{3}||_{2}^{2}\\\
=:&L(u^{left}(D(\boldsymbol{v}^{left})),\boldsymbol{T})+L(u^{right}(D(\boldsymbol{v}^{right})),\boldsymbol{T})\\\
+&L^{reg}(D(\boldsymbol{v}^{left}),D(\boldsymbol{v}^{right}))\\\
=:&L(u^{left},u^{right},\boldsymbol{v},\boldsymbol{T}).\end{split}$ (5)

If for a fixed $u$, standard gradient based optimization is used to find a
candidate for the minimizer $\boldsymbol{v}$ of
$f^{\boldsymbol{T}}(u,\boldsymbol{v})$, the vector found depends on the
initial point. One can try to minimize $\boldsymbol{v}\mapsto
f^{\boldsymbol{T}}(u,\boldsymbol{v})$ by doing optimization with several
different initial points, where the initial points are constructed from a
grid, and finally taking the argument which gives the smallest recorded cost.
This is highly time consuming since depending on the density of the grid there
may be thousands of optimization problems for each $u$. We consider a possible
choice of the grid in Appendix.

##  5 Learning the minimizer

As in [2] we consider how to learn the minimizer of (5). The learnt minimizer
$G(\cdot;\boldsymbol{\phi})$ takes $u^{s}$ as input and outputs
$\boldsymbol{v}^{s}:=G(u^{s};\boldsymbol{\phi})$, $s\in\\{right,left\\}.$ The
architecture of $G$ shown in Fig. 3 is of Siamese type so the weights of the
left and right branches are shared. The fully connected part has the last
layer $x\mapsto 6(-0.5+\sigma(x))$ so that when the $\tanh$ nonlinearity is
applied to the output of the fully connected part later in $D$, the output of
$\tanh$ has its range near $(-1,1)$. In Fig. 3, $D\circ
G(\cdot;\boldsymbol{\phi})$ is called the Localisation net in [10].

Figure 3: In the siamese network, the backbone consists of the convolutional
part of the Resnet18 network [9] followed by two additional $1\times 1$
convolutional layers with a ReLU nonlinearity between the additional layers.
If $N$ is the batch size, the shape of the Resnet part output changes as
$(N,512,7,7)\overset{\text{additional
layer}}{\longrightarrow}(N,128,7,7)\overset{\text{ReLU}}{\longrightarrow}(N,128,7,7)\overset{\text{additional
layer}}{\longrightarrow}(N,32,7,7)\overset{\text{ReLU}}{\longrightarrow}(N,32,7,7)\overset{\text{Normalize,
$dim=1$}}{\longrightarrow}(N,32,7,7)$. The last tensor is reshaped to
$(N,32\cdot 7\cdot 7)$ and goes through a feed-forward network, which has
three fully connected layers:
$(N,1568)\overset{\text{linear}}{\longrightarrow}(N,512)\overset{\text{ReLU}}{\longrightarrow}(N,512)\overset{\text{linear}}{\longrightarrow}(N,32)\overset{\text{ReLU}}{\longrightarrow}(N,32)\overset{\text{linear}}{\longrightarrow}(N,4)\overset{6(-\frac{1}{2}+\sigma(\cdot))}{\longrightarrow}(N,4).$
The last output corresponding to input $u^{s}$ is denoted by
$\boldsymbol{v}^{s}$, so $\boldsymbol{v}^{s}=G(u^{s};\boldsymbol{\phi})$.
Also, $\boldsymbol{\theta}^{s}:=D(\boldsymbol{v}^{s}).$ Figure 4: Outline of
the two-phase algorithm. Template $\boldsymbol{T}$ and its sub-template
$\tilde{\boldsymbol{T}}$ are extracted from the same X-ray image.

We denote
$G(u;\boldsymbol{\phi}):=(G(u^{left};\boldsymbol{\phi}),G(u^{right};\boldsymbol{\phi}))$.
We learn the parameters $\boldsymbol{\phi}$ by minimizing the function

$\mathcal{L}(\boldsymbol{\phi}):=E_{u}[c_{u}f^{\boldsymbol{T}}(u,G(u;\boldsymbol{\phi}))]$
(6)

where $E$ denotes the expectation and $u$ is sampled from the space of
bilateral X-ray knee images. We add a scaling factor $c_{u}$, $c_{u}>0$, in
the expectation since the value of
$\min_{\boldsymbol{v}}f^{\boldsymbol{T}}(u,\boldsymbol{v})$ is not constant
with respect to $u$ and as a result, in training, if $B$ is a minibatch, then
in $1/|B|\sum_{i\in B}f^{\boldsymbol{T}}(u_{i};G(u_{i};\boldsymbol{\phi}))$ a
possible much larger than others summand can dominate and training tries to
make this large value as small as possible. To treat minimization problems
more equally, we add the scaling $c_{u}$. A possible choice of $c_{u}$ could
be $c_{u}=1/\min_{\boldsymbol{v}}f^{\boldsymbol{T}}(u,\boldsymbol{v})$. But,
as was argued earlier, this is slow to compute, so we use an approximation to
it.

We do the search of the optimal patch in two phases. We outline this in Fig. 4
and Algorithm 2. First, we learn to search the template $\boldsymbol{T}$ from
the input images $u^{left}$, $u^{right}$. This phase produces a network
$G(\cdot;\boldsymbol{\phi}_{coarse})$. The input images for the second phase
are slightly enlargened $u^{left}(D(G(u^{left};\boldsymbol{\phi}_{coarse})))$
and $u^{right}(D(G(u^{right};\boldsymbol{\phi}_{coarse})))$. We denote the
enlargened images by $\widetilde{u^{left}}$ and $\widetilde{u^{right}}.$ From
these images the second network $G(\cdot;\boldsymbol{\phi})$ is taught to
extract patches that look like the template $\tilde{\boldsymbol{T}}$ of Fig.
4.

Algorithm 1 A Phase Training

0: $\text{Training data }\\{u_{i}\\}_{i=1}^{N},\text{ number of outer
iterations }$

$T,\text{ number of epochs }M,\text{ template }\boldsymbol{T}$

$t\leftarrow 0$

while $t\neq T$ do

if $t==0:$ then

Set $c_{u_{i}}=1$ for all $u_{i}$

else

For each $u_{i}:=(u^{left}_{i},u^{right}_{i})$ do sharpening: use
$\boldsymbol{v}_{0}:=(G(u^{left}_{i};\boldsymbol{\phi}^{t}),G(u^{right}_{i};\boldsymbol{\phi}^{t}))$
as an initial point for the minimization of
$f^{\boldsymbol{T}}(u_{i},\boldsymbol{v}).$ Use the Adam optimizer with $300$
iterations. Record all values of
$f^{\boldsymbol{T}}(u_{i},\boldsymbol{v}_{j})$ for $j=0,1,\dotsc,300$, and set
$c_{u_{i}}=1/\min_{j}f^{\boldsymbol{T}}(u_{i},\boldsymbol{v}_{j})$.

end if

$t\leftarrow t+1$

$\text{Update }\boldsymbol{\phi}\text{ in }\eqref{eq:L_phi}\text{ by running
}M\text{ epochs with minibatch}$

$\text{size }8\text{ over the training data.}$ We get
$G(\cdot;\boldsymbol{\phi}^{t}).$

end while

return $G(\cdot;\boldsymbol{\phi}^{T})$

Algorithm 2 Two Phases Training

Phase 1: Using $\\{u_{i}:=(u^{left}_{i}$, $u^{right}_{i})\\}_{i=1}^{500}$ and
the template $\boldsymbol{T}$ from Fig. 4, apply Algorithm 1 with $T=3$,
$M=15$, to obtain $G(\cdot;\boldsymbol{\phi}^{coarse})$.

Phase 2: Enlarge the images
$u_{i}^{left}(D(G(u_{i}^{left};\boldsymbol{\phi}^{coarse})))$ and
$u_{i}^{right}(D(G(u_{i}^{right};\boldsymbol{\phi}^{coarse})))$ by $1.2$ for
all $i=1,\dotsc,500,$ and rename the enlarged pairs as $u_{i}$, select the
template $\widetilde{\boldsymbol{T}}$ from Fig. 4 as $\boldsymbol{T}$, $T=3$,
$M=15$, and apply Algorithm 1 to obtain $G(\cdot;\boldsymbol{\phi})$.

We train the neural model using $500$ bilateral Osteoarthritis Initiative
(OAI) images. The images are from the baseline i.e. from the first visit. We
minimize the energy using the Adam optimizer with batch size $8$, the learning
rate $1e-5$ for the Resnet18 backbone parameters and $1e-4$ for the rest of
the parameters. We use PyTorch and we initialize the Resnet18 parameters using
its pre-trained version.

##  6 Experimental results

We compare four methods in the numerical experiments. In the experiments, we
first compare minimums found by different methods by examining the magnitudes
of the non-regularized loss values
$\mathcal{L}^{left}(P^{left},\tilde{\boldsymbol{T}})$ and
$\mathcal{L}^{right}(P^{right},\tilde{\boldsymbol{T}})$ in (5), and also the
sum values
$\mathcal{L}^{left}(P^{left},\tilde{\boldsymbol{T}})+\mathcal{L}^{right}(P^{right},\tilde{\boldsymbol{T}})$
where $P^{left}$ and $P^{right}$ are the detected areas. We also visually
compare whether there is a correspondence between how small the loss value is
and how good the obtained knee joint area localization is. We do not consider
for instance the classical intersection over union values, since the
continuous domain methods we consider also allow for rotations, and it is
quite difficult and laborious to annotate possibly rotated bounding boxes, the
step we wanted to avoid in the first place.

The baseline template matching method we consider is based on Python’s scikit-
image’s match$\\_$template function which is based on [12]. We use the method
to compare in a sliding window manner the template $\tilde{\boldsymbol{T}}$ in
Fig. 4 to five differently scaled versions of the input image and we extract
the patches corresponding to the smallest left plus right loss value. The
second method we consider is the neural network method in Algorithm 2. The
third method is the second method followed by the sharpening step described
inside Algorithm 1, so the trained neural network method is used to get a good
initialization point which is followed by $300$ iterations with the Adam
optimizer and the patch corresponding to the smallest value recorded is taken
as the output patch of the third method. The fourth method is described at the
end of Section 4 and in Appendix.

In the first experiment we consider how small the minimums found by different
methods are on average between different methods. We consider images from the
OAI dataset and the Multicenter Osteoarthritis Study (MOST) dataset. The OAI
images are from the first visit but different from the training data. In
Tables 1 and 2 we consider the averages of the values
$\mathcal{L}^{left}(P^{left},\tilde{\boldsymbol{T}})+\mathcal{L}^{right}(P^{right},\tilde{\boldsymbol{T}})$
by three different methods. It seems that in the OAI dataset the baseline
template matching method does not have dense enough scale grid or since the
method is non-regularized and does not allow for rotations, the method does
not find as small minimums as the neural method or its sharpening step. In the
MOST dataset the neural method does not directly find as small minimums as in
the OAI dataset, but the method produces such good initial points that the
sharpening step again finds much smaller minimums.

In Fig. 5 we show some typical test detections by different methods and the
corresponding one-sided values
$\mathcal{L}^{left}(P^{left},\tilde{\boldsymbol{T}})$ and
$\mathcal{L}^{right}(P^{right},\tilde{\boldsymbol{T}})$. From the experiments
we made, it seems that in the baseline template matching method, there is more
variability in the patches found and there is more vertical translations
between the sides compared to other methods. In contrast, the neural method
and its sharpening step produce detections which are more uniform. This is
again due to regularization incorporated into the continuous domain methods.
It seems that the very time consuming fourth method based on the grid
initialized search does not necessarily find smaller loss values than the
neural network method followed by the sharpening step. Also, in the neural
network method, the parameters describing the knee joint areas are based on
the last convolutional layer features of the Resnet18 architecture which are
not necessarily so sensitive to small variations in the input image. This
might explain why the fourth method does not produce as uniform detections as
the neural method.

Table 1: The sides combined average values over $400$ test images in the OAI dataset. Method |  Left$+$right loss
---|---
Template matching |  0.221
Neural method |  0.159
Neural method + Sharpening step |  0.141
Table 2: The sides combined average values over $400$ test images in the MOST dataset. Method |  Left$+$right loss
---|---
Template matching |  0.249
Neural method |  0.224
Neural method + Sharpening step |  0.157

|
---|---
|
|
|
|
|
|
|

$0.1212$

$0.0604$

$0.0563$

$0.056$

$0.0717$

$0.0686$

$0.0622$

$0.0626$

$0.0892$

$0.078$

$0.0658$

$0.0654$

$0.1083$

$0.0766$

$0.0715$

$0.0714$

$0.099$

$0.0663$

$0.062$

$0.0628$

$0.1064$

$0.0949$

$0.0822$

$0.0628$

$0.2258$

$0.106$

$0.0873$

$0.0803$

$0.1474$

$0.1047$

$0.0905$

$0.0803$

$0.0872$

$0.0917$

$0.0362$

$0.0366$

$0.0985$

$0.1649$

$0.0667$

$0.065$

$0.0941$

$0.0587$

$0.0528$

$0.0525$

$0.0753$

$0.0723$

$0.0586$

$0.0573$

$0.1164$

$0.1393$

$0.1142$

$0.0966$

$0.3015$

$0.1196$

$0.0939$

$0.0964$

$0.08$

$0.232$

$0.0878$

$0.0727$

$0.0863$

$0.1122$

$0.0729$

$0.0669$

Figure 5: The detections and the corresponding one-sided loss values obtained
by four different methods on eight test images. The vertical bar separates the
sides. The first and fifth column: the template matching method. The second
and sixth column: the neural method. The third and seventh column: the neural
method followed by the sharpening step. The fourth and eighth column: the grid
initialized search algorithm.

##  7 Conclusion

In this paper, we considered detecting knee joint areas from bilateral knee
X-ray images. For the annotations, we used just one side of a single bilateral
image where we manually selected the templates. Otherwise, the training data
in the teaching of the neural method was bilateral X-ray knee images without
any labels. In the experiments it seemed that the neural network based method
can give detections which are more uniform than detections from a classical
version of the template matching. Overall, since the knee joint area varies
quite rigidly between different knee images and there are no big deformations
due to the knee anatomy across the images, it seems that in the knee joint
area detection, it is not necessary to have many manually annotated images,
one template image may be enough, provided the task-dependent energy function
is properly chosen.

## Acknowledgments

This work was done in a project funded by the Research Unit of Mathematical
Sciences and the Research Unit of Medical Imaging, Physics and Technology at
the University of Oulu, Finland.

## Appendix: Grid initialized search algorithm

In this section, we consider a grid initialized search algorithm to minimize
$f^{\boldsymbol{T}}(u,\boldsymbol{v})$ in (5). We consider the template
$\boldsymbol{T}$ and the input halves as in Fig. 2. Let
$s\in\\{left,right\\}$. If $\boldsymbol{\theta}^{s}=D(\boldsymbol{v}^{s}),$
then $\boldsymbol{v}^{s}$ is an unconstrained vector and
$\boldsymbol{\theta}^{s}$ describes the scale, translation and rotation. We
form a grid for $\boldsymbol{\theta}^{s}$ and use the corresponding
$\boldsymbol{v}^{s}=D^{-1}(\boldsymbol{\theta}^{s})$ as an initial point for
the unconstrained optimization problem. Here $\boldsymbol{\theta}^{s}_{1}$
determines the scale of the patch we try to match to the template,
$(\boldsymbol{\theta}_{2}^{s},\boldsymbol{\theta}_{3}^{s})$ determines the
center of the patch. We set $\boldsymbol{\theta}^{s}_{4}=0$ so there is no
rotation in the initialization. We consider a few choices for the scale and
for each scale we consider a dense enough grid for the translations such that
adjacent rectangles are close enough. If the scale is large, i.e.
$\boldsymbol{\theta}^{s}_{1}\approx 1$, then there is no need for so many
translation grid points compared to when the scale is small, i.e.
$\boldsymbol{\theta}^{s}_{1}\approx 0$.

First, we consider scales. We split the interval
$[\alpha_{0},\alpha_{0}+\beta_{0}]$, which is the range of (1), into $N_{s}$
number of equal distant points. Then we get
$s_{1}=\alpha_{0},\dots,s_{N_{s}}=\alpha_{0}+\beta_{0}$. How to initialize
$d_{j}:=(\boldsymbol{v}^{s}_{j})_{1}$ such that we get $s_{j}$? Now

$\displaystyle s_{j}$
$\displaystyle=\alpha_{0}+\beta_{0}\left(\frac{1}{2}\left(1+\tanh(d_{j})\right)\right)$
$\displaystyle\Leftrightarrow d_{j}$
$\displaystyle=\tanh^{-1}\left(\frac{2(s_{j}-\alpha_{0})}{\beta_{0}}-1\right).$

For translations, for a given scale $s_{i}$, we must have that
$-s_{i}+\boldsymbol{\theta}^{s}_{2},s_{i}+\boldsymbol{\theta}^{s}_{2}\in[-1,1]$
and
$-s_{i}/f+\boldsymbol{\theta}^{s}_{3},s_{i}/f+\boldsymbol{\theta}^{s}_{3}\in[-1,1]$
in order for the extracted patch to lie inside the input image.

We get that we must have $s_{i}-1\leq\boldsymbol{\theta}^{s}_{2}\leq 1-s_{i}$.
We split the interval $[s_{i}-1,1-s_{i}]$ into $N_{\theta}$ equally spaced
points $\theta_{13}^{i,j}$, $j=1,\dots,N_{\theta}$, such that
$\theta_{13}^{i,1}=s_{i}-1,\dotsc,\theta_{13}^{i,N_{\theta}}=1-s_{i}$ and
$\Delta^{1}$ denotes the distance between two consecutive grid points. How to
choose the number of points $N_{\theta}$? We have
$\theta_{13}^{i,j}=(s_{i}-1)+(j-1)\cdot\Delta^{1}$. Two $2s_{i}$ sized
adjacent line segments have an overlap of size $2s_{i}-\Delta^{1}$, when
$\Delta^{1}<2s_{i}$. If two adjacent line segments are required to have an
overlap whose length is at least $2s_{i}-r2s_{i}$, where $0<r<1$, we get that
$\Delta^{1}\leq r2s_{i}$. We set $\Delta^{1}=r2s_{i}$ with $r=0.25$ which
leads to the choice $N_{\theta}=1+ceil(\frac{2-2s_{i}}{r2s_{i}})$. We note
that $N_{\theta}$ increases as $s_{i}$ decreases which is what was desired.

How to initialize $(\boldsymbol{v}^{s}_{i,j})_{2}$ such that
$\theta_{13}=(\boldsymbol{\theta}_{i,j}^{s})_{2}=\theta_{13}^{i,j}$? From (2),

$\displaystyle\theta_{13}^{i,j}$
$\displaystyle=2(1-s_{i})a_{13}^{i,j}-1+s_{i}$ $\displaystyle\Leftrightarrow
a_{13}^{i,j}$ $\displaystyle=\frac{\theta_{13}^{i,j}+1-s_{i}}{2(1-s_{i})}$
$\displaystyle\Leftrightarrow(\boldsymbol{v}_{i,j}^{s})_{2}$
$\displaystyle=\tanh^{-1}\left(\frac{\theta_{13}^{i,j}+1-s_{i}}{1-s_{i}}-1\right)$

and analogously for $\boldsymbol{v}_{3}$.

We finally consider the combination of the sides initializations. If
$[-s_{i}^{left}+\theta^{left,i,j}_{1,3},s^{left}_{1}+\theta^{left,i,j}_{1,3}]\subset[-1,1]$
is the projection of a left initialization rectangle into the horizontal axis,
we initialize the patch in the right side by using the same vertical
initialization, but horizontally we test every initialization rectangle whose
horizontal center is within the interval
$[\theta^{left,i,j}_{1,3}-\frac{1}{3},\theta^{left,i,j}_{1,3}+\frac{1}{3}]$.

## References

  * [1] P. Adarsh, P. Rathi, and M. Kumar.  Yolo v3-tiny: Object detection and recognition using one stage improved model.  In 2020 6th International Conference on Advanced Computing and Communication Systems (ICACCS), pages 687–694, 2020.
  * [2] M. Andrychowicz, M. Denil, S. Gómez, M. W. Hoffman, D. Pfau, T. Schaul, B. Shillingford, and N. de Freitas.  Learning to learn by gradient descent by gradient descent.  In D. Lee, M. Sugiyama, U. Luxburg, I. Guyon, and R. Garnett, editors, Advances in Neural Information Processing Systems, volume 29. Curran Associates, Inc., 2016.
  * [3] J. Antony, K. McGuinness, K. Moran, and N. E. O’Connor.  Automatic detection of knee joints and quantification of knee osteoarthritis severity using convolutional neural networks.  In MLDM, 2017.
  * [4] J. Antony, K. McGuinness, K. Moran, and N. E. O’Connor.  Feature Learning to Automatically Assess Radiographic Knee Osteoarthritis Severity, pages 9–93.  Springer International Publishing, Cham, 2020.
  * [5] J. Antony, K. McGuinness, N. E. O’Connor, and K. A. Moran.  Quantifying radiographic knee osteoarthritis severity using deep convolutional neural networks.  2016 23rd International Conference on Pattern Recognition (ICPR), pages 1195–1200, 2016.
  * [6] P. Chen, L. Gao, X. Shi, K. Allen, and L. Yang.  Fully automatic knee osteoarthritis severity grading using deep neural networks with a novel ordinal loss.  Computerized Medical Imaging and Graphics, 75:84–92, 2019.
  * [7] L. A. Gatys, A. S. Ecker, and M. Bethge.  Image style transfer using convolutional neural networks.  In 2016 IEEE Conference on Computer Vision and Pattern Recognition (CVPR), pages 2414–2423, 2016.
  * [8] H. Gu, K. Li, R. J. Colglazier, J. Yang, M. Lebhar, J. O’Donnell, W. A. Jiranek, R. C. Mather, R. J. French, N. Said, J. Zhang, C. Park, and M. A. Mazurowski.  Automated grading of radiographic knee osteoarthritis severity combined with joint space narrowing.  CoRR, abs/2203.08914, 2022.
  * [9] K. He, X. Zhang, S. Ren, and J. Sun.  Deep residual learning for image recognition.  In 2016 IEEE Conference on Computer Vision and Pattern Recognition (CVPR), pages 770–778, 2016.
  * [10] M. Jaderberg, K. Simonyan, A. Zisserman, and k. kavukcuoglu.  Spatial transformer networks.  In C. Cortes, N. Lawrence, D. Lee, M. Sugiyama, and R. Garnett, editors, Advances in Neural Information Processing Systems, volume 28. Curran Associates, Inc., 2015.
  * [11] J. H. Kellgren and J. S. Lawrence.  Radiological assessment of osteo-arthrosis.  Ann Rheum Dis., 16(4):494–502, 1957.
  * [12] J. Lewis.  Fast normalized cross-correlation.  Ind. Light Magic, 10, 2001.
  * [13] E. Pierson, D. Cutler, J. Leskovec, S. Mullainathan, and Z. Obermeyer.  An algorithmic approach to reducing unexplained pain disparities in underserved populations.  Nature Medicine, 27:136–140, 01 2021.
  * [14] J. Redmon and A. Farhadi.  Yolo9000: Better, faster, stronger.  In 2017 IEEE Conference on Computer Vision and Pattern Recognition (CVPR), pages 6517–6525, 2017.
  * [15] L. Shamir, S. Ling, W. Scott, M. Hochberg, L. Ferrucci, and I. Goldberg.  Early detection of radiographic knee osteoarthritis using computer-aided analysis.  Osteoarthritis and Cartilage, 17(10):1307–1312, 2009.
  * [16] A. Tiulpin, J. Thevenot, E. Rahtu, and S. Saarakkala.  A novel method for automatic localization of joint area on knee plain radiographs.  In Scandinavian Conference on Image Analysis, pages 290–301. Springer, 2017.
  * [17] Y. Wang, Z. Bi, Y. Xie, T. Wu, X. Zeng, S. Chen, and D. Zhou.  Learning from highly confident samples for automatic knee osteoarthritis severity assessment: Data from the osteoarthritis initiative.  IEEE Journal of Biomedical and Health Informatics, pages 1–1, 2021\.