A Novel Similarity Measure for Image Sequences

TL;DR

This paper introduces a new global similarity measure for image sequences based on Schatten-q-norms of gradient matrices, enabling unbiased registration of dynamic and serial images with promising initial results.

Contribution

It proposes a novel similarity measure that considers all images simultaneously, improving registration accuracy for sequences and serial sections.

Findings

Effective registration of DCE-MRI sequences demonstrated.

Successful serial section alignment shown in experiments.

Preservation of global data structure in registration results.

Abstract

Quantification of image similarity is a common problem in image processing. For pairs of two images, a variety of options is available and well-understood. However, some applications such as dynamic imaging or serial sectioning involve the analysis of image sequences and thus require a simultaneous and unbiased comparison of many images. This paper proposes a new similarity measure, that takes a global perspective and involves all images at the same time. The key idea is to look at Schatten-q-norms of a matrix assembled from normalized gradient fields of the image sequence. In particular, for q = 0, the measure is minimized if the gradient information from the image sequence has a low rank. This global perspective of the novel SqN-measure does not only allow to register sequences from dynamic imaging, e.g. DCE-MRI, but is also a new opportunity to simultaneously register serial sections, e.g. in histology. In this way, an accumulation of small, local registration errors may be avoided. First numerical experiments show very promising results for a DCE-MRI sequence of a human kidney as well as for a set of serial sections. The global structure of the data used for registration with SqN is preserved in all cases.

Equations

$$A=U\mathrm{diag}(S)V^{\top},\mbox{ with }U^{\top}U=E_{n},\ V^{\top}V=E_{T}.\$$

$$\Vert A{\Vert }_{S,q}^q:=\Vert S{\Vert }_q^q={\sum }_i{\sigma }_i^q\text{mbox}forq\ge 0.$$

$$\mathrm{S}q\mathrm{N}(I_1,\dots ,I_T)$$

$$A(x)$$

$${\eta }_t$$

$$A^{\top}A=\Big{(}\begin{array}[]{rr}1&\alpha\\ \alpha&1\end{array}\Big{)}\mbox{ with singular values }\sigma_{1}^{2}=1+\alpha\mbox{ and }\sigma_{2}^{2}=1-\alpha.$$

$$\mathcal{J}(y_1,\dots ,y_T)=\mathrm{S}q\mathrm{N}(I_1(y_1),\dots ,I_T(y_T))+\alpha {\sum }_{t=1}^T\mathcal{S}(y_t).$$

$$\mathcal{J}^{\mathrm{NGF}}(y_{1},\dots,y_{T})=\textstyle\sum_{t=1}^{T-1}\Big{\{}\ \mathrm{NGF}(I_{t}(y_{t}),I_{t+1}(y_{t+1}))+\alpha\ \mathcal{S}(y_{t})\Big{\}}$$

$$y_t^{k+1}={\mathrm{argmin}}_{y_t}{\mathcal{J}}^{\mathrm{NGF}}(y_1^{k+1},\dots ,y_{t-1}^{k+1},y_t,y_{t+1}^k,\dots ,y_T^k),t=1,\dots ,T.$$

Full text

This paper proposes a new similarity measure, that takes a global perspective and involves all images at the same time. The key idea is to look at Schatten-$q$-norms of a matrix assembled from normalized gradient fields of the image sequence. In particular, for $q=0$, the measure is minimized if the gradient information from the image sequence has a low rank.

This global perspective of the novel $\mathrm{S}q\mathrm{N}$-measure does not only allow to register sequences from dynamic imaging, e.g. DCE-MRI, but is also a new opportunity to simultaneously register serial sections, e.g. in histology. In this way, an accumulation of small, local registration errors may be avoided.

First numerical experiments show very promising results for a DCE-MRI sequence of a human kidney as well as for a set of serial sections. The global structure of the data used for registration with $\mathrm{S}q\mathrm{N}$ is preserved in all cases.

1 Introduction

Quantification of image similarity is a common problem in image processing. For pairs of two images, a variety of options such as sum of squared differences (SSD), normalized gradient fields (NGF), or mutual information (MI) exists and these measures are well-understood; see e.g. [12, 13, 17]. However, some applications such as dynamic imaging or serial sectioning involve the analysis of image sequences and thus require a simultaneous and unbiased comparison of many images.

In some of these applications, image intensity may changes over time. For example, the glomerular filtration rate (GFR) is an important parameter for kidney malfunction [18]. The GFR might be determined on the basis of a time series of dynamic contrast-enhanced magnetic resonance images (DCE-MRI). This sequence then needs to be registered in order to correct for motion artifacts; see, e.g., [7, 9]. Another example is the analysis of a histological serial sectioning. Here, the staining of sections might be different and/or can express severe variations.

For applications like these, a proper image similarity measure is crucial. A standard approach is to perform a sequential comparison of pairs of two images from the sequence. However, a sequential registration is restricted to local rather than global information. Moreover, there might be issues choosing a suited starting image and determining the order of the sequence. Results may depend on these choices. There also exists a variety of statistical approaches for global image registration; see e.g. [5, 10, 16, 21]. However, we focus on a new global measure which is based on deterministic image features.

In this paper, we propose the new $\mathrm{S}q\mathrm{N}$ similarity measure, that is designed to compare a complete image sequence simultaneously and thus automatically distributing information in a global way. Our key idea is to look at Schatten-$q$-norms (more precisely: Schatten-$q$-quasinorms) of a matrix that is assembled from normalized gradient fields of the image sequence. Particularly for $q=0$, the measure is minimized for sequences of images where the gradient matrix has low rank and the approach is thus connected to principal component analysis and sparsity.

Our idea is motivated from color image denoising; see Möllenhoff et al. [14]. In that context, similar concepts are used as a regularization for TV denoising, and the gradient matrix is formed directly from gradients of the three color channels. In our paper, we interpret the individual images from a sequence as individual channels and use a Schatten-$q$-norm of the matrix of normalized gradients as a data fitting term rather than as a regularizer. As we will show, this can be viewed as a natural extension of NGF [6] and can thus deal with multi-modal frames. We remark that the concept also relates to ideas in video compression.

Our paper is organized as follows: At first we describe the novel distance measure and its relation to NGF as well as its characteristics. Moreover, we show numerical results for DCE-MRI time series and a H&E stained histological serial section of a mouse brain. We compare the performance of $\mathrm{S}q\mathrm{N}$ with NGF (DCE-MRI) and the well-known SSD (serial sectioning). Our examples show that $\mathrm{S}q\mathrm{N}$ results at least comparable registration results but is about six times faster as the competitive approaches. To conclude our paper, we discuss the numerical results and give a brief outlook on what our next steps are.

2 The novel similarity measure $\mathrm{S}q\mathrm{N}$

Our new distance measure is motivated by a regularizer for color image denoising; see [14]. The underlying idea is that in natural images, gradients of the different color channels are linearly dependent; see also Fig. 1. Therefore, an appropriate measure of dependency such as Schatten-$q$-norms [24] are excellent regularizers in color image denoising. This idea can be generalized to more than three channels and is therefore useful in applications such as parameter estimation in DCE-MRI [8] or the registration of multiple images as they appear for example in serial sectioning or time series [11].

We motivate our extension starting with a conceptual simpler but computational infeasible approach. The main point is to motivate the use of Schatten-$q$-norms. We then present a computational tractable version that is based on local image gradients. In contrast to [14] where a similar measure is used as regularizer, we also propose to use our new functional as a distance measure.

We recall that any matrix $A\in\mathbb{R}^{n,T}$ has a singular value decomposition (SVD) [4],

[TABLE]

Here, $E_{d}$ denotes the $d$-by-$d$ identity matrix, $S=(\sigma_{1},\ldots,\sigma_{\min\{n,T\}})$ is a vector with ordered entries $\sigma_{j}\geq\sigma_{j+1}\geq 0$ and $\mathrm{diag}(S)\in\mathbb{R}^{n,T}$ denotes a diagonal matrix. Using the SVD, the Schatten-$q$-(quasi)-norm of $A$ is then defined as

[TABLE]

An optimal choice of $q$ is obviously application dependent. In particular for registration problems, it is topic of current research. Note that for $q=0$, the measure counts non-zero entries and is not a norm but a so-called quasinorm. With $q=0$, we thus promote sparsity of $S$ and hence low rank of $A$. However, optimization of the [math]-quasinorm is non-trivial and it is typically replaced by the minimization of the $1$-norm; see [2]. Following [14], we use $q=1/2$ in this paper.

We are now ready to describe our novel similarity measure. To this end, we assume that $T\in\mathbb{N}$ discrete images $I_{t}$ are given, where each $I_{t}$ is of dimension $m_{1}\times\cdots\times m_{d}\in\mathbb{N}^{d}$ and $d\in\mathbb{N}$ denotes the spatial dimension. Let $n:=m_{1}\cdots m_{d}$ and $I_{t}$ be reshaped as $n\times 1$ array. The first idea is to look at the rank of $A=[I_{1},\ldots,I_{T}]\in\mathbb{R}^{n,T}$ as an indicator for the linear dependency of the images. Note that this approach is similar to a principal component analysis of the data. However, this approach is not suitable for images with varying intensities such as DCE-MRI or serial sections. Although this approach is conceptually appealing, it is computationally challenging as the complexity of the SVD is $\mathcal{O}\left(\min\{nT^{2},Tn^{2}\}\right)$ [4].

We escape the complexity problem by applying the measure to local structures. More precisely, we define our new distance measure $\mathrm{S}q\mathrm{N}$ for a sequence of $T$ images by

[TABLE]

Here, $\eta>0$ is a parameter discriminating signal from noise; see also [6]. If the noise level is unknown, we pick a small value, e.g. $\eta=10^{-5}$ in our experiments where every entry of $I_{t}$ lies in $[0,256)$.

We now outline the connection to NGF [6]. For two images $I_{1}$ and $I_{2}$, fixed $x$, $\eta=0$ and with $\alpha:=|\cos\angle(\nabla I_{1},\nabla I_{2})|$, we have

[TABLE]

Particularly for $q=0$, the distance is minimal if the image gradients are collinear, i.e. $\alpha=1$. Therefore, our new measure might be interpreted as a generalization of the normalized gradient field based distance measure [6]. Analogous arguments hold for $0\leq q<2$. Remarkably, for $q=2$, the energy function is constant and the measure thus meaningless. For $q>2$, the function is minimal for perpendicular gradients.

3 Numerical Results for Dynamic Imaging and Serial Sections

We now present results for the registration of DCE-MRI sequences of a human kidney and a histological serial sectioning of a mouse brain.

We start with registrations of DCE-MRI sequences of a human kidney; data courtesy of Jarle Rørvik, Haukeland University Hospital Bergen, Norway. Here, 3D images are taken at 49 time points. The objective is to register the images while maintaining the dynamics. More precisely, we use 146-by-82 coronal slices of a 146-by-82-by-52-by-49 volume for z-Slice 25; see Fig. 2.

Our registration scheme is based on the variational image registration framework FAIR [13]. More precisely, we minimize

[TABLE]

For ease of presentation, we use the curvature regularizer $\mathcal{S}(y)=\int(\Delta y)^{2}\ dx$ with regularization parameter $\alpha=0.1$ [1]. Optimization is performed in a standard way using a Gauß-Newton algorithm with Armijo linesearch [15] within a multi-level framework [13]. All computations are performed using MATLAB.

Fig. 2 (top row) shows three representative coronal slices of the original dataset. The slices correspond to different times during acquisition. Different phases of contrast agent uptake are visible, particularly within the kidney.

Fig. 3 displays sagittal and axial slices of the same volume, two each. The top row shows non-registered slices where motion artefacts are clearly visible. The middle row shows corresponding $\mathrm{S}q\mathrm{N}$-registration results for $q=0.5$ and the bottom row shows results for sequential NGF, respectively. More precisely, we optimized

[TABLE]

using alternating optimization, i.e.

[TABLE]

As to be expected, the $\mathrm{S}q\mathrm{N}$ and NGF results are very similar. However, within our non-optimized MATLAB framework, the $\mathrm{S}q\mathrm{N}$-registration is about six times faster than the NGF registration — even if the alternating optimization approach is limited to a forward-backward sweep.

Fig. 2 also illustrates the intensity variations in the original and registered data. It is apparent that intensity variations due to motion have been reduced tremendously. Note that the variations, in particular in kidney and bladder (partly visible), are still visible and therefore the schemes allow for a subsequent dynamical analysis.

The $\mathrm{S}q\mathrm{N}$ measure is also capable of serial section registration. Fig. 5 displays results for a H&E stained histological serial section of a mouse brain; 189 sections of 512-by-512 pixel; data courtesy of O. Schmitt, University of Rostock, Germany; see [19] for experimental details. As this data has been normalized, we compare an $\mathrm{S}q\mathrm{N}$ registration with a sequential SSD based approach. The figure displays the original data (top row), $\mathrm{S}q\mathrm{N}$ results (middle row), and SSD results (bottom row). It is apparent that registration can reconstruct the local brain structure. Here, the SSD results appear to be slightly more blurred.

Note that a sequential approach involves an alternating optimization framework and is based on a fixed initial and final slice to avoid the so-called banana-effect [3, 19, 20, 23], i.e. a global drift of structures due to sequential registraton. The sequential registration process may accumulate small errors that can cause a major drift of the overall structure; see [23] for examples.

In contrast, the $\mathrm{S}q\mathrm{N}$ approach enables a global optimization which in addition can be performed in just one pass. Moreover, our experiments indicate that $\mathrm{S}q\mathrm{N}$ does not introduce global drifts.

Fig. 4 shows the results of a translation experiment for the distance measures SSD, MI, NGF as wells as for $\mathrm{S}q\mathrm{N}$ for three different configurations of the parameter $q$. It is apparent that $q$ has an impact on the global minimum of the energy in this experiment. All measures have the same minimizer which is achieved for the template image in its origin when it is not shifted as shown in the difference image of Fig. 4.

4 Discussion and Conclusions

The novel $\mathrm{S}q\mathrm{N}$ image similarity measure has been proposed. The new measure is motivated from color image denoising [14]. The main idea is to quantify structural image information expressed by normalized intensity gradients with Schatten-q-(quasi)-norms. For $q=0$, the measure quantifies sparsity, i.e. redundancy of information in an image sequence. Moreover, using the normalized gradient fields as structural information, the focus is on image structures and not on intensity. Therefore, the new measure is effective in a multi-modal setup and it might be interpreted as an extension of NGF [6].

The novel measure considers the image sequence as a whole and therefore has no bias towards a particular ordering, which is common in a sequential setup. Therefore, the new measure provides also a global information transport, which might be beneficial for particular applications. In particular, it omits an unwanted drift of structure as it is very common in sequential approaches; so-called banana-effect [20].

We show the potential of the new measure in a registration setup. More precisely, we show how easily the new measure can be integrated in an existing registration framework such as FAIR [13]. Both, the measure and its analytic derivative are computed and thus, efficient optimization techniques can be used.

Exemplarily, we demonstrate the power of the measure in a registration problem for a dynamic contrast enhanced MRI sequence and histological serial sectioning. As it turns out, the new scheme produces results of at least comparable quality much faster (about six times faster in our experiments).

Particularly for serial section registrations, we escape an additional iterative process on the sequences of sections and make use of global information transport.

Our next steps include to come up with a more efficient implementation which is suitable in a 3D setup. We remark that the main computational cost is not the computation of the SVD of the gradient matrix. This involves the computation of eigenvalues of a $T$-by-$T$ matrix, where $T$ is the number of images in the sequence. In particular, $T$ does not dependent on the spatial dimension of the data or the spatial resolution. However, the computation of the gradient matrix $A$ does.

Moreover, we will investigate the impact of the distance measure in the analysis of the dynamic information within the DCE-MRI setting. With this knowledge, we would be able to quantify an optimal choice for the parameter $q$. In our current experiments, it appears that the impact of $q$ is rather small unless we pick values close to $q=2$. We also will compare the measure to global statistical measures (e.g. [5, 10, 16, 21]) already used in this field of registration and image processing.

An implementation of $\mathrm{S}q\mathrm{N}$ will soon be available within the FAIR-Toolbox; see https://github.com/C4IR.

Acknowledgement

The authors acknowledge the financial support by the Federal Ministry of Education and Research of Germany in the framework of MED4D (project number 05M16FLA).