# Subpixel motion estimation for video-based target-free vibration monitoring under complex environmental conditions

Zhilei Luo, Boualem Merainani, Vincent Baltazart, Qinghua Zhang, Michael Döhler

<figure><img src="image_1.png" /></figure>

## ► To cite this version:

Zhilei Luo, Boualem Merainani, Vincent Baltazart, Qinghua Zhang, Michael Döhler. Subpixel motion estimation for video-based target-free vibration monitoring under complex environmental conditions. Mechanical Systems and Signal Processing, 2025, 226, pp.112342. 10.1016/j.ymssp.2025.112342. hal-04894400

HAL Id: hal-04894400

https://inria.hal.science/hal-04894400v1

Submitted on 17 Jan 2025

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Subpixel Motion Estimation for Video-based Target-free Vibration Monitoring under Complex Environmental Conditions

Zhilei Luo<sup>a</sup>, Boualem Merainani<sup>b</sup>, Vincent Baltazart<sup>b</sup>, Qinghua Zhang<sup>a</sup>, Michael Döhler<sup>a,*</sup>

<sup>a</sup>Univ. Gustave Eiffel, Inria, COSYS-SII, 14S, 35042 Rennes, France

<sup>b</sup>Univ. Gustave Eiffel, Inria, COSYS-SII, 14S, 44344 Bouguenais, France

## Abstract

An emerging technique to measure structural vibrations is based on motion signals extracted from video images. The use of video cameras offers multiple advantages over traditional mechanical sensors: contactless measurement, large coverage, easy installation and maintenance. This paper proposes a new method for real-time motion signal extraction from video images with subpixel accuracy. It aims to address the challenges posed by complex operating conditions, namely illumination variations and background interference. Illumination robustness is achieved by efficiently combining image intensity interpolation with an affine brightness and contrast tuning transformation. Background robustness is obtained by automatically selecting active pixels in the processed images. Moreover, the on-line numerical computations are fast enough for real-time applications. Results of theoretical analysis ensure that the considered optimization criteria are well-posed and that, in the main step of subpixel motion estimation, the involved nonlinear optimization problem is efficiently solved in closed-form through a linear least squares problem. The performance of the proposed method is evaluated both on simulated images and on laboratory experiments with a target-free cantilever beam in comparison with existing methods. The reported results demonstrate the robustness and computational efficiency of the proposed method under complex environmental conditions, allowing real-time computation on a standard laptop for vibration monitoring with more than 100 virtual sensors at 600 frames per second.

**Keywords:** video images, motion signal extraction, subpixel accuracy, real-time, illumination and background robustness, vibration monitoring

## 1. Introduction

### 1.1. Background

Digital image-based motion extraction has been developed as a powerful photogrammetric tool for displacement measurement on real world structures, with the advantages of its flexibility and efficiency in

<footnote>*Corresponding author; E-mail address: michael.doehler@inria.fr</footnote>
operation [1]. It has been applied for various practical purposes, such as motion tracking and compensation [2], remote sensing [3], satellite imagery [4]. Recently, this methodology has particularly regained attention of researchers in the field of vision-based vibration analysis [5–9], e.g., for the purpose of civil structure health monitoring. Using cameras in such applications has multiple advantages, like contact-less measurement, large coverage and easy maintenance. However, it is sensitive to fluctuations of ambient light conditions, which imply image brightness and contrast variations. Complex background of the scene captured by the camera also interferes with vision-based motion estimation. Robustness to real world environmental conditions is of particular importance for reliable long term and long range monitoring of civil engineering structures. Moreover, algorithms designed for vibration analysis must be fast enough for real time applications.

Digital image correlation (DIC) techniques have been widely developed for displacement extraction and surface deformation measurement [10–16]. With such techniques, subpixel accuracy is obtained by extending image cross-correlation to non integer image shifts in pixel unit, either by upsampling somehow the processed images, or by interpolating the 2D integer grid of image cross-correlations [17–19]. In recent years, DIC techniques have become popular for applications to vibration measurement and analysis [20–24]. However, these typical DIC schemes involve time-consuming iterative nonlinear optimization algorithms, which are incompatible with real time vibration monitoring.

Recently, phase-based methods have been developed for vibration monitoring [25, 26]. Among them, the phase only correlation method [27, 28] considers the global phase information based on the Fourier transform of the processed image, discarding completely the amplitude information. Besides, an alternative approach is based on local phase in images extracted by two-dimensional band-pass filters, such as Gabor filter and steerable filter [29–33]. Although the phase based methods are impressive for their detectability of subtle vibration [34, 35] and insensitivity to image intensity fluctuations [36–38], the performance of these techniques highly depends on image phase extraction, of which robustness and efficiency need further investigations [39, 40].

Besides DIC and phase based methods, optical flow methods are also well spread for vision based motion extraction [41, 42]. The fundamental principle of these methods relies on the assumption of stationary illumination conditions so that the image brightness level and contrast corresponding to the same point on the observed target remain constant within the processed image sequence. Particularly, two main techniques are prevalently adopted for vibration estimation: image intensity gradient oriented optical flow estimation [43–45], and Kanade-Lucas-Tomasi (KLT) method [46–48]. Nevertheless, these optical flow methods still need to be fully studied and demonstrated for their robustness and efficiency in vibration estimation applications.

*The purpose of this paper is to propose a new motion estimation method robust to light variations and interference background patterns, with a numerical efficiency compatible with real-time processing of video images for vibration monitoring of mechanical or civil structures as detailed below. This proposed method*
<figure><img src="image_3.png" /><figcaption>Figure 1: Motion estimation from digital images. Active pixels are automatically selected within pre-selected ROIs. One of the ROIs is illustrated with a yellow square. The proposed ALISSE method estimates displacement through image correlation analysis based on selected active pixels.</figcaption></figure>

is shortly outlined in the next subsection, before being presented in detail in Section 2.

## 1.2. Problem statement and summary of the proposed subpixel motion extraction scheme

Consider a vibrating structure observed by a camera, as illustrated in Figure 1. Within the captured image sequence, pre-selected regions of interests (ROI) are analyzed for motion extraction. Traditionally, displacement, velocity or acceleration measurements were obtained by mechanical sensors installed on the monitored structure. The information collected by each mechanical sensor would correspond to the motion extracted from an ROI by processing the video sequence. The vibrating structure is locally rigid enough so that its shape deformation observed within each ROI is insignificant. In this way, it is possible to replace mechanical sensors by a camera with specially developed image processing algorithms. It is assumed that rotational motion is negligible, as well as the motion along the camera optical axis direction.

This paper proposes a new digital image correlation method for vision based vibration monitoring, with a particular attention focused on robustness to complex operational conditions, including illumination variations and interfering background patterns. It is fast enough for real time processing on an ordinary laptop computer, and can be easily parallelized for large scale applications. For example, more than 100 ROIs could be processed in real time when processing video frames recorded at 600 Hz.

This new method, called Adaptive Light intensity Interpolation for Subpixel Shift Estimation (ALISSE), is intended to work in the presence of image brightness and contrast variations, as well as interfering backgrounds, while accommodating both large and small vibration amplitudes. At every sampling instant, displacements are estimated between the current frame and a reference frame (typically the first frame of the processed video sequence). For each considered ROI, the pixel level displacement is first estimated by maximizing the cross-correlation between the current frame and the reference frame. The current frame is then shifted so that the remaining displacement between the shifted frame and the reference frame is smaller than one pixel. At the subpixel level, the inter-pixel gray level of the reference image frame is
interpolated in order to extend the cross-correlation to subpixel shifts. Instead of maximizing this extended cross-correlation, the brightness and contrast of the interpolated reference image are adjusted by an affine transformation, while matching the transformed interpolated image to the current image frame through a subpixel shift in the least squares sense. It is the simultaneous optimization of the subpixel shift and the brightness-contrast adjustment that leads to a subpixel displacement estimation robust to brightness and contrast variations. Robustness to interfering background patterns is gained by computing the least-squares fit of two images with automatically selected active pixels in each ROI.

The focus of this paper is on algorithm design and performance evaluation for motion estimation. For reliable real-time applications, the main steps of the designed algorithm are theoretically analyzed, to ensure that the considered optimization criteria are well-posed and that they do allow solving the formulated motion estimation problem. Since the final goal of the proposed ALISSE method is vision-based vibration monitoring, its accuracy is evaluated under complex operation conditions involving light variations and interfering background, in comparison with existing algorithms developed for the same purpose, based on both simulated data and laboratory experiments. Numerical computation costs are also taken into account in these comparisons. The estimated displacement data will serve as input to a multitude of methods for structural vibration analysis, see e.g. [49–51].

This paper is organized as follows. The proposed new motion estimation algorithm is presented in detail in Section 2. Important properties of the main steps of the method are theoretically analyzed in Section 3. Then, in Section 4, the performance of the algorithm is evaluated through numerically simulated video images. Section 5 reports laboratory experiments conducted on a target-free cantilever beam. Finally, concluding remarks are drawn in Section 6.

## 2. The proposed motion estimation method: ALISSE

For the purpose of vibration analysis, motion signals of the monitored structure are extracted from a video sequence at pre-selected positions, where mechanical sensors would be installed or embedded in traditional vibration monitoring technologies. To extract the displacement signal at every pre-selected position, the real-time processing of each new image frame mainly consists of two steps: the pixel level displacement is first estimated by cross-correlation maximization, then the result is refined at the subpixel level based on image inter-pixel interpolation, by efficiently combining subpixel displacement estimation with brightness and contrast adjustment. Both steps rely on the processing of active pixels automatically pre-selected before the real-time processing of the video sequence. This section presents in detail ALISSE, along with a summarizing flow chart diagram in Figure 2.
<figure><img src="image_4.png" /><figcaption>Figure 2: Flow chart diagram of ALISSE, the proposed motion estimation method.</figcaption></figure>

## 2.1. Active pixel selection

Within the scope of Structural Health Monitoring (SHM) applications, motion estimation is typically conducted on selected regions of interest (ROI) along the mechanical structure being monitored. ROIs generally include image textures or spatial features to facilitate motion estimation. In certain scenarios, motion estimation is focused on a set of effective pixels corresponding to spatial features, such as corner points, edges, or contours on the moving structure [52–54]. These effective pixels can be identified from the magnitude of the image gradient or through steerable filters [30].

In contrast, the proposed ALISSE method leverages a preprocessing technique described in [55] to select active pixels within the ROI. This technique employs an FFT-based spectral analysis of the gray-level time variations for each pixel within the ROI. The energy spectrum generated from this analysis identifies the most vibrating pixels, thereby enabling more efficient motion estimation in terms of both accuracy and computational load. The method in [55] further enhances the edge pixels of the moving structure, effectively reducing interference from complex static backgrounds. Moreover, the spatial distribution of these selected active pixels can roughly infer the motion direction and amplitude range.

Each selected active pixel being referenced by its coordinates $(i, j)$ in pixel unit, in each ROI, the set of $N$ active pixels will be referred to with

$$A = \{(i_1, j_1), (i_2, j_2), \dots, (i_N, j_N)\}. \qquad (1)$$
## 2.2. Pixel level estimation based on cross-correlation maximization

When a video sequence is processed for motion extraction, the displacement to be estimated is defined between the current frame and a reference frame, typically the first one of the sequence. For ease of presentation, the time instants of the reference and current frames are denoted by $t_1$ and $t_2$, respectively. While $t_1$ is a constant for the processed video, $t_2$ generically represents the time instant of any current frame in the image sequence. It is assumed that each of the considered ROIs contains (part of) a moving object. The displacement is estimated based on active pixels selected within each ROI.

The pixel level displacement, i.e., the displacement rounded to an integer in pixel unit, is estimated by maximizing the cross-correlation between the corresponding sets of active pixels within the considered ROI in the current image frame and in the reference image frame. The zero-mean normalized cross-correlation (ZNCC) [56, 57] is adopted for better accuracy. In each ROI, let $I(i, j, t_1)$ denote the gray level of the image pixel indexed by the integer pair $(i, j)$ at instant $t_1$, where $(i, j) \in A$, with $A$ representing the set of coordinates of active pixels of the ROI, as in (1), and similarly $I(i, j, t_2)$ for the image at instant $t_2$. Then ZNCC is defined, for integer horizontal and vertical shifts $m$ and $n$ (in pixel unit), as:

$$Z(m, n) = \frac{\sum_{(i,j) \in A} [I(i, j, t_1) - \bar{I}_{t_1}] [I(i + m, j + n, t_2) - \bar{I}_{t_2}]}{\sqrt{\sum_{(i,j) \in A} [I(i, j, t_1) - \bar{I}_{t_1}]^2 \sum_{(i,j) \in A} [I(i + m, j + n, t_2) - \bar{I}_{t_2}]^2}}, \quad (2)$$

where $\bar{I}_{t_1}$ and $\bar{I}_{t_2}$ are respectively the empirical mean gray levels of the active pixels of the ROI at instants $t_1$ and $t_2$. The integer pair $(m, n)$ maximizing $Z(m, n)$ corresponds to the estimated pixel level displacement, namely:

$$(m^*, n^*) = \arg \max_{(m,n) \in S} Z(m,n), \quad (3)$$

where $S$ is a chosen set of integer pairs defining the search area for $(m, n)$ maximizing $Z(m, n)$.

## 2.3. Local linear interpolation for subpixel accuracy

After pixel level displacement estimation, the ROI in the current image frame (at instant $t_2$) is shifted accordingly, so that the remaining shift between the current image and the reference image is smaller than one pixel. This shift is illustrated in Figure 3, as commented in detail below.

On the left, it is an ROI in the reference frame (at $t_1$), containing part of the monitored moving structure, and part of the noisy background; on the right, it is the corresponding ROI in the frame captured at $t_2$, both the original one and the shifted one. The intersecting points of thin vertical and horizontal lines represent pixels, while the square areas between these lines are the inter-pixel regions, where the light intensity should be interpolated. The highlighted orange rectangles represent the active pixels in the considered ROI. For the frame at $t_2$, in order to compensate the pixel level motion (horizontal in this example), the original
<figure><img src="image_5.png" /><figcaption>Figure 3: On the left, an ROI in the reference frame (at $t_1$), containing part of the monitored moving structure, and part of the noisy background; on the right, the corresponding ROI in the frame captured at $t_2$, both the original one and the shifted one. The intersecting points of thin vertical and horizontal lines represent pixels, while the square areas between these lines are the inter-pixel regions, where the image intensity should be interpolated. The highlighted orange rectangles represent the active pixels in the considered ROI. For the frame at $t_2$, in order to compensate the pixel level motion (horizontal in this example), the original ROI (dashed yellow box) is shifted according to the estimated pixel level displacement. The shifted ROI is represented by the dark blue box.</figcaption></figure>

ROI (dashed yellow box) is shifted according to the estimated pixel level displacement. The shifted ROI is represented by the dark blue box. The dashed vertical magenta line indicates an edge of the structure, while the red thin vertical line corresponds to the position of the edge at $t_1$. Of course, at $t_1$ these two lines are at the same position. At $t_2$, their distance has a big part of several pixel units, and the remaining part smaller than one pixel unit (subpixel displacement). This subpixel displacement is represented, in the right side picture, by the distance between the dashed vertical magenta line and the red thin vertical line. From $t_1$ to $t_2$, the active pixels maintain their relative positions.

After the ROI shift, the shifted image intensity at $t_2$ will be denoted by:

$$I_s(i, j, t_2) \triangleq I(i + m^*, j + n^*, t_2), \quad (i, j) \in \Omega. \quad (4)$$

It corresponds to a change of pixel coordinates. $\Omega$ represents the area of the considered ROI at $t_1$. Its position has a shift ($m^*, n^*$) at $t_2$. The new notation $I_s(i, j, t_2)$ will be used instead of $I(i + m^*, j + n^*, t_2)$ in subpixel shift estimation, in order to simplify expressions involving the integer values $m^*$ and $n^*$.

The image intensity $I(i, j, t_1)$ is only available on an integer grid of pixels indexed with $(i, j)$, so is
I$_s(i, j, t_2)$. For the purpose of subpixel shift estimation, the image intensity within each ROI in the reference image $I(i, j, t_1)$ will be interpolated to inter-pixel intensity values. In the proposed ALISSE method, there is no need to interpolate every image frame in the processed video sequence. For numerical efficiency, only the reference image frame will be interpolated.

Let $\hat{I}(i+x, j+y, t_1)$ denote the intensity of the reference image interpolated around an active pixel indexed by $(i, j)$, with $(x, y)$ corresponding to real number coordinates in a neighborhood of the origin $(0, 0)$, i.e.,

$$ (x, y) \in \mathbb{H} \triangleq \{(x, y) \in \mathbb{R}^2 : |x| < 1, |y| < 1\}. \quad (5) $$

The way in which the inter-pixel intensity $\hat{I}(i+x, j+y, t_1)$ is computed from available values of $I(i, j, t_1)$ is essential in this approach, and its choice strongly impacts the performance of the resulting subpixel displacement estimation method. In order to efficiently estimate subpixel displacements, in ALISSE the computation of $\hat{I}(i+x, j+y, t_1)$ is based on a local linear interpolation in the neighborhood of $(i, j)$, i.e., for $(x, y) \in \mathbb{H}$,

$$ \hat{I}(i+x, j+y, t_1) = a_{i,j}x + b_{i,j}y + c_{i,j} \quad (6) $$

where $a_{i,j}$, $b_{i,j}$ and $c_{i,j}$ are local linear interpolation coefficients determined by fitting this linear model to the available image intensity values in the neighborhood of $(i, j)$. Because this interpolation model is fitted and applied locally around $(i, j)$, the coefficients $a_{i,j}$, $b_{i,j}$ and $c_{i,j}$ are indexed by the integer coordinates $i$ and $j$. The advantage in terms of numerical efficiency due to this choice of linear interpolation will become clear later in the following subsections.

The linear interpolation (6) is similar to the first order Taylor expansion as used in [18], in the sense that they both use a plane to locally approximate the image intensity surface $I(i+x, j+y, t_1)$. The Taylor expansion corresponds to a tangent plane of $I(i+x, j+y, t_1)$, usually estimated by finite difference involving the 5 pixels in the neighborhood of $(i, j)$. In the proposed ALISSE method, the linear interpolation (6) will be based on the plane fitted in the least squares sense to the same 5 pixels, namely the pixels indexed by $(i, j)$, $(i+1, j)$, $(i-1, j)$, $(i, j+1)$ and $(i, j-1)$. More specifically, given an active pixel indexed by $(i, j) \in \mathbb{A}$, the local linear model (6) is fitted by minimizing the following least squares criterion:

$$ \min_{a_{i,j}, b_{i,j}, c_{i,j}} \sum_{(p,q) \in \mathcal{V}} [I(i+p, j+q, t_1) - (a_{i,j}p + b_{i,j}q + c_{i,j})]^2 \quad (7) $$

where

$$ \mathcal{V} = \{(0,0), (-1,0), (1,0), (0,-1), (0,1)\}. \quad (8) $$

This is a standard linear least squares problem, whose closed-form solution is well known. It will be explicitly written later in equation (40), when the well-posedness of this least squares problem is analyzed.
The resulting linear interpolation (6) is valid locally for $(x, y) \in \mathbb{H}$, with $\mathbb{H}$ defined in (5). Its purpose is the interpolation around $(i, j)$, i.e., for computing $\hat{I}(i+x, j+y, t_1)$ following (6), with $(x, y) \in \mathbb{H}$. To interpolate the image around all the active pixels, the linear interpolation (6) has to be fitted and implemented for every active pixel indexed by $(i, j) \in \mathbb{A}$.

## 2.4. Subpixel shift estimation with robustness to brightness and contrast variations

To deal with image intensity fluctuation in practical applications, a simple solution is image intensity normalization of each image frame [58]. However, this pre-processing only normalizes each image at the pixel level. To gain robustness at the subpixel level, the proposed new method ALISSE will involve an affine transformation adjusting image brightness and contrast. For the interpolated reference image intensity around active pixels $\hat{I}(i+x, j+y, t_1)$, the affine transformation consists of:

$$ \hat{I}(i+x, j+y, t_1) \rightarrow \alpha \hat{I}(i+x, j+y, t_1) + \beta, \qquad (9) $$

in which $\alpha$ and $\beta$ are scalar real numbers adjusting the contrast and the brightness of the image, respectively. See, e.g., [59], for a similar brightness and contrast parametrization.

In order to efficiently estimate the subpixel shift while adjusting $\alpha$ and $\beta$, the following least squares problem is considered:

$$ \min_{x,y,\alpha,\beta} \sum_{(i,j) \in \mathbb{A}} [I_s(i,j,t_2) - \alpha(a_{i,j}x + b_{i,j}y + c_{i,j}) - \beta]^2, \qquad (10) $$

in which the linear interpolation $\hat{I}(i+x, j+y, t_1) = a_{i,j}x + b_{i,j}y + c_{i,j}$ has been inserted, with the coefficients $a_{i,j}$, $b_{i,j}$ and $c_{i,j}$ already computed by solving the least squares problem (7) for each active pixel $(i, j) \in \mathbb{A}$. Through the optimization problem (10), the interpolated reference image $\hat{I}(i+x, j+y, t_1)$ is fitted to the current image $I_s(i, j, t_2)$ by simultaneously adjusting the subpixel shift $(x, y)$ and the contrast-brightness parameters $\alpha, \beta$. Remark that the two optimization problems (7) and (10), though similar, have been formulated for different purposes. The first one computing the interpolator coefficients $a_{i,j}$, $b_{i,j}$, $c_{i,j}$ fits the interpolator to the neighborhood of each active pixel, whereas the second one estimates the subpixel displacement $(x, y)$ while adjusting the contrast and brightness parameters $\alpha$ and $\beta$ for all active pixels in the ROI.

It appears that Equation (10) refers to a nonlinear minimization problem, because of the products $\alpha x$ and $\alpha y$ involving the unknowns $\alpha, x$ and $y$. Instead of directly solving such a nonlinear problem, let us define two auxiliary variables:

$$ x^a \triangleq \alpha x, \quad y^a \triangleq \alpha y, \qquad (11) $$
and consider the minimization problem:

$$ (\hat{x}^a, \hat{y}^a, \hat{\alpha}, \hat{\beta}) = \arg \min_{x^a, y^a, \alpha, \beta} \sum_{(i,j) \in A} [I_s(i, j, t_2) - (a_{i,j}x^a + b_{i,j}y^a + c_{i,j}\alpha) - \beta]^2, \quad (12) $$

in which the previous products of unknowns $\alpha x$ and $\alpha y$ have been replaced by $x^a$ and $y^a$, so that the nonlinear terms $\alpha a_{i,j}x + \alpha b_{i,j}y$ (involving products of unknowns $\alpha x$ and $\alpha y$) become the linear terms $a_{i,j}x^a + b_{i,j}y^a$. It results in a standard quadratic criterion in terms of the unknowns $x^a, y^a, \alpha, \beta$ expressed in (12). Its minimization is easily solved in closed-form by the well known least squares solution. Define the matrix containing the interpolation coefficients already compute by solving (7):

$$ X = \begin{bmatrix} a_{i_1,j_1} & b_{i_1,j_1} & c_{i_1,j_1} & 1 \\ a_{i_2,j_2} & b_{i_2,j_2} & c_{i_2,j_2} & 1 \\ \vdots & \vdots & \vdots & \vdots \\ a_{i_N,j_N} & b_{i_N,j_N} & c_{i_N,j_N} & 1 \end{bmatrix}, \quad (13) $$

and the vector

$$ \Upsilon(t_2) = \begin{bmatrix} I_s(i_1, j_1, t_2) \\ I_s(i_2, j_2, t_2) \\ \vdots \\ I_s(i_N, j_N, t_2) \end{bmatrix}, \quad (14) $$

in which the index pairs $(i_1, j_1), \dots, (i_N, j_N)$ belong to the set of active pixels $A$ as defined in (1). Then, denote the pseudo-inverse of matrix $X$ by

$$ X^\dagger = (X^T X)^{-1} X^T, \quad (15) $$

the solution $(\hat{x}^a, \hat{y}^a, \hat{\alpha}, \hat{\beta})$ minimizing (12) is given by

$$ [\hat{x}^a \quad \hat{y}^a \quad \hat{\alpha} \quad \hat{\beta}]^T = X^\dagger \Upsilon(t_2). \quad (16) $$

Since only the reference image $I(i, j, t_1)$ is interpolated in ALISSE, the interpolation coefficients in $X$ computed by solving (7) depend only on the reference image frame. In contrast, the vector $\Upsilon(t_2)$ involves the current image frame at $t_2$. For this reason, the matrix $X$ and its pseudo inverse $X^\dagger$ are computed only once before the online processing of current frames, and the online computation for processing each image frame is mainly the last matrix product in (16).

If $\hat{\alpha} \neq 0$, following (11), let:

$$ \hat{x} \triangleq \frac{\hat{x}^a}{\hat{\alpha}}, \quad \hat{y} \triangleq \frac{\hat{y}^a}{\hat{\alpha}}. \quad (17) $$

Then the resulting $\hat{x}, \hat{y}, \hat{\alpha}, \hat{\beta}$ solve the nonlinear minimization problem (10), provided $\hat{\alpha} \neq 0$. The case $\hat{\alpha} = 0$ is unlikely in practice, since it would imply a zero contrast fitted image. If ever $\hat{\alpha} = 0$ happens, then the
criterion in (10) is minimized by $\alpha^* = 0$ and $\beta^*$ equal to the average value of $I_s(i, j, t_2)$ for $(i, j) \in A$. See Proposition 3 of the next subsection for explanations. When the fitted image has zero contrast, subpixel displacement cannot be estimated, then by default the subpixel displacement $(x^*, y^*) = (0, 0)$. As already stated above, in practice, it is unlikely that the optimally fitted image has a zero contrast.

The transformed optimization problem (12) involves the local linear interpolator $\hat{I}(i + x, j + y, t_1)$. Actually, according to (11),

$$a_{i,j}x^a + b_{i,j}y^a + c_{i,j}\alpha = (a_{i,j}x + b_{i,j}y + c_{i,j})\alpha \quad (18)$$

$$= \hat{I}(i + x, j + y, t_1)\alpha. \quad (19)$$

The linearity of the interpolator, i.e., the fact that $\hat{I}(i + x, j + y, t_1) = a_{i,j}x + b_{i,j}y + c_{i,j}$ depends linearly on $(x, y)$, is important so that the minimized criterion in (12) is quadratic in its unknowns $x^a, y^a, \alpha, \beta$. It is this quadratic form of the criterion that leads to an efficient analytic solution of the optimization problem with ALISSE, i.e., the standard closed-form least squares solution (16), in the online computation. If a nonlinear interpolator was used, then less efficient iterative nonlinear optimization algorithms would have to be applied for minimizing the resulting nonlinear criterion instead of (10) or (12).

## 2.5. When the motion direction is known

In the previous subsections, there was no need to know the motion direction to estimate the horizontal and vertical displacement components. In some applications, the motion direction is known in advance, then it can be used to simplify computations. It is assumed in this subsection that the angle between the motion direction and the x-axis is $\theta \in (-\pi/2, \pi/2]$, with positive values of $\theta$ corresponding to the first quadrant of the Cartesian coordinates system. Displacements along this direction can be positive or negative, hence there is no need to consider larger angles.

## 2.6. Pixel level estimation

For the pixel level displacement estimation by ZNCC maximization (Equation (2)), when angle $\theta$ of the motion direction is known, the search of the maximum is reduced to 1D instead of the 2D search in the general case. First assume that $|\theta| < \pi/4$. The case $|\theta| \ge \pi/4$ will be considered later.

Let $\bar{m}$ be the chosen integer value so that the ZNCC maximum will be searched for integer values of $m$ within the set:

$$M \triangleq \{-\bar{m}, -\bar{m} + 1, -\bar{m} + 2, \dots, \bar{m}\}. \quad (20)$$

The corresponding integer values of $n$ so that the integer pairs $(m, n)$ stay around the straight-line in the direction of $\theta$ are $\lfloor m \tan(\theta) \rfloor$ and $\lfloor m \tan(\theta) \rfloor$, which denote respectively the largest integer not larger than
$m \tan(\theta)$ and the smallest integer not smaller than $m \tan(\theta)$. Then the ZNCC maximum will be searched as:

$$ (m^*, n^*) = \arg \max_{m \in \mathcal{S}_\theta} Z(m, n) \qquad (21) $$

with

$$ \mathcal{S}_{\theta} \triangleq \{(m, n) : m \in \mathbb{U}, n = \lfloor m \tan(\theta) \rfloor\} \cup \{(m, n) : m \in \mathbb{U}, n = \lceil m \tan(\theta) \rceil\}. \qquad (22) $$

In the case with $|\theta| \ge \pi/4$, the definition of $\mathcal{S}_\theta$ is replaced by:

$$ \mathcal{S}_{\theta} \triangleq \{(m, n) : n \in \mathbb{N}, m = \lfloor n \cot(\theta) \rfloor\} \cup \{(m, n) : n \in \mathbb{N}, m = \lceil n \cot(\theta) \rceil\}. \qquad (23) $$

with

$$ \mathbb{N} \triangleq \{-\bar{n}, -\bar{n} + 1, -\bar{n} + 2, \dots, \bar{n}\} \qquad (24) $$

and some chosen integer $\bar{n}$.

In the particular case of horizontal motion with $\theta = 0$, there is no vertical displacement, hence $n^* = 0$ and ZNCC maximization is simplified by searching only in the horizontal direction. Then Equation (3) is replaced by:

$$ m^* = \arg \max_{m \in \mathbb{M}} Z(m, 0) \qquad (25) $$

where $\mathbb{M}$ denotes the integers corresponding to the horizontal search range, and $Z(m, 0)$ is as defined in (2).

Similarly, in the case of vertical motion with $\theta = \pi/2$, $m^* = 0$ and:

$$ n^* = \arg \max_{n \in \mathbb{N}} Z(0, n). \qquad (26) $$

## Subpixel level estimation

At the subpixel level, the local linear interpolation presented in Section 2.3 remains unchanged. For the next step, the least squares problem (12) is replaced by:

$$ (\hat{z}^a, \hat{\alpha}, \hat{\beta}) = \arg \min_{z^a, \alpha, \beta} \sum_{(i,j) \in \mathcal{A}} [I_s(i, j, t_2) - (a_{i,j} z^a \cos(\theta) + b_{i,j} z^a \sin(\theta) + c_{i,j} \alpha) - \beta]^2. \qquad (27) $$

Then the subpixel displacement is estimated by:

$$ \hat{x} = \frac{\hat{z}^a \cos(\theta)}{\hat{\alpha}}, \quad \hat{y} = \frac{\hat{z}^a \sin(\theta)}{\hat{\alpha}}, \qquad (28) $$

provided $\hat{\alpha} \ne 0$.

In the particular case with $\theta = 0$, the subpixel estimation is simplified as:

$$ (\hat{z}^a, \hat{\alpha}, \hat{\beta}) = \arg \min_{z^a, \alpha, \beta} \sum_{(i,j) \in \mathcal{A}} [I_s(i, j, t_2) - (a_{i,j} z^a + c_{i,j} \alpha) - \beta]^2. \qquad (29) $$
$$ \hat{x} = \frac{\hat{z}^a}{\hat{\alpha}}, \quad \hat{y} = 0. \qquad (30) $$

Similarly, when $\theta = \pi/2$,

$$ (\hat{z}^a, \hat{\alpha}, \hat{\beta}) = \arg \min_{z^a, \alpha, \beta} \sum_{(i,j) \in A} [I_s(i, j, t_2) - (b_{i,j}z^a + c_{i,j}\alpha) - \beta]^2, \qquad (31) $$

$$ \hat{x} = 0, \quad \hat{y} = \frac{\hat{z}^a}{\hat{\alpha}}. \qquad (32) $$

### 3. Theoretical analysis of least squares estimation in the proposed method

The essential numerical computations in the proposed subpixel displacement estimation method are the solutions of the two least squares problems formulated in (7) and (12). It is important to ensure that these least squares problems are well-posed, in the sense that there exists a solution and the solution is unique, for each of the two considered problems. Moreover, it should be formally proved that the nonlinear optimization problem (10) is indeed solved through the standard least squares problem (12). Such results are presented below through 3 propositions and their formal proofs.

**Proposition 1.** For any shifted image pixel intensity value $I(i, j, t_1)$ and any intensity values at the 4 surrounding pixels, or in other words, for any $I(i + p, j + q, t_1)$ with $(p, q) \in \mathbb{V}$ as defined in (8), the least squares problem formulated in (7) has a unique solution.

Before proving this result, let us first introduce some notations. The elements in the set $\mathbb{V}$ defined in (8) are referenced by the double index $(p, q)$. In order to use a single index in some analysis, it is redefined as

$$ \mathbb{V} = \{v_1 = (0,0), v_2 = (-1,0), v_3 = (1,0), v_4 = (0,-1), v_5 = (0,1)\}. \qquad (33) $$

Then the single index $k = 1, 2, \dots, 5$ can be used for referring to $v_k$. Accordingly, vectors $\gamma_k$ and a matrix $\Gamma$ are defined such that

$$ \Gamma = \begin{bmatrix} \gamma_1^T \\ \gamma_2^T \\ \gamma_3^T \\ \gamma_4^T \\ \gamma_5^T \end{bmatrix} = \begin{bmatrix} 0 & 0 & 1 \\ -1 & 0 & 1 \\ 1 & 0 & 1 \\ 0 & -1 & 1 \\ 0 & 1 & 1 \end{bmatrix}. \qquad (34) $$

#### Proof of Proposition 1

To apply Lemma 1 (see Appendix A) to the least squares problem formulated in (7), the linear interpolation model involved in (7), which is related to the linear regression formulated in (A.1), is now written
in vector form as

$$a_{i,j}p + b_{i,j}q + c_{i,j} = \begin{bmatrix} p & q & 1 \end{bmatrix} \begin{bmatrix} a_{i,j} \\ b_{i,j} \\ c_{i,j} \end{bmatrix}. \qquad (35)$$

Here the pair $(i, j)$ is fixed. In fact, the linear interpolation is made around each active pixel $(i, j) \in \mathbb{A}$, but let us focus on one of them and the associated least squares problem (7). However, the pair $(p, q)$ varies within $\mathbb{V}$, which was redefined in (33). To explicitly show the varying $(p, q)$, let

$$(p(k), q(k)) = v_k, \qquad (36)$$

with $v_k$ as defined in (33), then (35) becomes

$$a_{i,j}p(k) + b_{i,j}q(k) + c_{i,j} = \begin{bmatrix} p(k) & q(k) & 1 \end{bmatrix} \begin{bmatrix} a_{i,j} \\ b_{i,j} \\ c_{i,j} \end{bmatrix} = \gamma_k^T \begin{bmatrix} a_{i,j} \\ b_{i,j} \\ c_{i,j} \end{bmatrix} \qquad (37)$$

with $\gamma_k$ as introduced in (34).

290

With respect to the linear regression (A.1) in Lemma 1, the regressor vector $\psi_k$ corresponds to $\gamma_k$, and the regression coefficient vector $\eta$ corresponds to $[a_{i,j}, b_{i,j}, c_{i,j}]^T$. Then the matrix sum expressed in (A.3) is

$$\sum_{k=1}^{5} \psi_k \psi_k^T = \sum_{k=1}^{5} \gamma_k \gamma_k^T = \Gamma^T \Gamma = \begin{bmatrix} 2 & 0 & 0 \\ 0 & 2 & 0 \\ 0 & 0 & 5 \end{bmatrix} \qquad (38)$$

where $\Gamma$ is as in (34). This result is clearly an invertible matrix. Therefore, the considered least squares problem (7) has a unique solution, according to Lemma 1. $\square$

By the way, let $\zeta(i, j, t_1)$ denote the vector composed of the light intensity $I(i, j, t_1)$ and those of the 4 surrounding pixels, i.e.,

$$\zeta(i, j, t_1) \triangleq \begin{bmatrix} I(i, j, t_1) \\ I(i-1, j, t_1) \\ I(i+1, j, t_1) \\ I(i, j-1, t_1) \\ I(i, j+1, t_1) \end{bmatrix}, \qquad (39)$$

then the solution of the least squares problem formulated in (7) is given by

$$\begin{bmatrix} a_{i,j} \\ b_{i,j} \\ c_{i,j} \end{bmatrix} = (\Gamma^T \Gamma)^{-1} \Gamma^T \zeta(i, j, t_1). \qquad (40)$$

This result will be useful later.
**Proposition 2.** Let $\zeta(i, j, t_1)$ denote the vector composed of the image intensity $I(i, j, t_1)$ and those of the 4 surrounding pixels, as defined in (39). Consider the set of such vectors corresponding to $(i, j) \in A$. Assume that the empirical covariance matrix of these vectors $\zeta(i, j, t_1)$ is invertible, then the least squares problem formulated in (12) has a unique solution.

While in the previous proposition the well-posedness of the least squares problem (7) is always true, the well-posedness of the least squares problem (12) requires some condition. Intuitively, the processed images should contain motion information to enable displacement estimation. A trivial counterexample is an image whose pixels all have the same light intensity. The condition on the empirical covariance matrix of the vectors $\zeta(i, j, t_1)$ means that the light intensities of the active pixels and their surrounding pixels should vary sufficiently. More specifically, let $\zeta_1, \zeta_2, \dots, \zeta_N$ denote the $N$ vectors $\zeta(i, j, t_1)$ corresponding to $(i, j) \in A$, then their empirical covariance matrix is

$$ \frac{1}{N-1} \sum_{k=1}^{N} (\zeta_k - \bar{\zeta})(\zeta_k - \bar{\zeta})^T \qquad (41) $$

with

$$ \bar{\zeta} \triangleq \frac{1}{N} \sum_{k=1}^{N} \zeta_k. \qquad (42) $$

In (41), each term $(\zeta_k - \bar{\zeta})(\zeta_k - \bar{\zeta})^T$ either increases the rank of the matrix sum by one or leave the rank unchanged. Adding such a term can never decrease the matrix rank. When the rank reaches 5, the $5 \times 5$ empirical covariance matrix is invertible. Therefore, for sufficiently varying $\zeta_k$ and sufficiently large $N$, the empirical covariance matrix (41) is invertible in most practical situations.

The proof of Proposition 2, technically more involved than the previous one, is presented in Appendix B.

**Proposition 3.** Under the same assumption as in Proposition 2, which ensures the existence and the uniqueness of the solution to the standard least squares problem (12), if $\hat{\alpha} \neq 0$ in the solution of (12), then the nonlinear optimization problem (10) has a unique solution $(\hat{x}, \hat{y}, \hat{\alpha}, \hat{\beta})$. This solution is related to the unique solution of the linear least squares problem (12) through (17). If $\hat{\alpha} = 0$ in the solution of (12), then the solution of (10) is not unique. Each of the solutions in this case, say $(x^*, y^*, \alpha^*, \beta^*)$, contains $\alpha^* = 0$ and

$$ \beta^* = \arg \min_{\beta} \sum_{(i,j) \in A} [I_s(i, j, t_2) - \beta]^2. \qquad (43) $$

The condition that $\hat{\alpha} \neq 0$ needs some explanation. According to Proposition 3, $\hat{\alpha} = 0$ in the solution of the linear optimization problem (12) implies that, in the solution of the nonlinear optimization problem (10), $\alpha^* = 0$, hence the optimally fitted image has a zero contrast. Such a situation is unlikely in practice.
The proof of Proposition 3 is presented in Appendix B.

The above propositions and their formal proofs confirm that the two least squares problems (7) and (12) are well posed, and that the nonlinear optimization problem (10) is indeed solved through the standard least squares problem (12). After this theoretic consolidation of the proposed new method, its performance will be evaluated in the next section, both through simulated data and through laboratory experiments.

# 4. Performance assessment on simulated data sets

In this section, the subpixel accuracy of the proposed Adaptive Light intensity Interpolation for Subpixel Shift Estimation (ALISSE) method will be evaluated on simulated image datasets, with a particular attention to its robustness regarding image intensity variations and interfering background patterns. The numerical computations are carried out in MATLAB using a laptop equipped with Intel processor i9-11950H.

## 4.1. Simulation and performance criteria

In each of the test cases presented in this section, grayscale images with dynamic range from 0 to 255 (8 digital bits) are numerically simulated. The image size is of 101×101 pixels, and a rectangular window of 59×59 pixels is chosen as region of interest (ROI), as shown in Figure 4(a). The simulated moving object is a rectangular bar of size 51×11 pixels at gray level 205, while the background is at gray level 60. Its intensity profile along the middle horizontal axis is presented in Figure 4(b). An image formation model proposed by Chojnacki and Szpak [60] is employed to generate the images with Gaussian noises, with respect to the decibel signal-to-noise (SNR) levels defined as follows [61]:

$$ \text{SNR} = 20 \log_{10} \frac{\max(I) - \min(I)}{\sigma_N} \qquad (44) $$

where $I$ is the processed image, and $\max(I)$ and $\min(I)$ are respectively its maximum and minimum values. $\sigma_N$ is the Gaussian noise level, and three noise levels are defined in the tests: $\sigma_N = 2, 4, 6$ and $8$, corresponding to SNR = 25.17 dB, 27.66 dB, 31.19 dB and 37.21 dB.

The objective motion is simulated as an oscillation. The motion direction angle $\theta = 2\pi/9$ (40°) is used in the simulation for data generation, but the value of $\theta$ is unknown to the evaluated motion estimation algorithm. The horizontal and vertical components of the simulated displacement, denoted respectively by $u(t)$ and $v(t)$, follow the equations:

$$ u(t) = \sum_{n}^{N_f} \frac{A_m}{N_f} \sin(2\pi f_n t) \cos(\theta), \quad v(t) = \sum_{n}^{N_f} \frac{A_m}{N_f} \sin(2\pi f_n t) \sin(\theta), \qquad (45) $$

where $A_m$ is the oscillation amplitude, $f_n$ are angular frequencies, and $N_f$ is the total number of frequency components. To simulate a subpixel vibration behavior, the following values are used to generate data: $A_m = 0.5$ pixel, $N_f=8$, $f_n=1.05$ Hz, 6.58 Hz, 18.44 Hz, 36.12 Hz, 59.71 Hz, 89.19 Hz, 124.55 Hz and 165.79
Hz (these frequencies are selected following the finite element simulation of the cantilever beam used in the subsequent experimental application [26]). The oscillation equations (45) and the involved numerical values of $\theta$, $A_m$, $f_n$, $N_f$ are unknown to the motion estimation method. The sampling rate of simulated data sets is set to 600 Hz, and data is simulated for a duration of 10 seconds, generating a 6,000-frames video (image sequence). For numerical efficiency, the computation of active pixels in the proposed algorithm only takes 600 frames out of the entire sequence, recording around 1 cycle of beam vibration.

The performances of displacement estimation are evaluated with the root mean square error (RMSE) in pixel unit:

$$ \text{RMSE} = \sqrt{\frac{\sum_{t=1}^{N_s} [\hat{u}(t) - u(t)]^2 + [\hat{v}(t) - v(t)]^2}{N_s}} \quad (46) $$

where $t$ is the image frame index in the simulated image sequence, $\hat{u}(t)$ and $\hat{v}(t)$ are the estimates of horizontal and vertical displacements $u(t)$ and $v(t)$, and $N_s$ is the total number of image frames.

In some simulations, the image intensity and background of frames in the simulated sequence are disturbed to test the robustness of the algorithm to light variations and to non-uniform background. Global contrast and brightness of each image frame are modulated as in (9) with

$$ \alpha(t) = A_{\alpha} \sin(2\pi f_{\alpha} t + \phi_{\alpha}) + Q_{\alpha}, \quad \beta(t) = A_{\beta} \sin(2\pi f_{\beta} t + \phi_{\beta}) + Q_{\beta}, \quad (47) $$

where the variation amplitudes and frequencies are set to $A_{\alpha} = 0.040$, $A_{\beta} = 10$, $f_{\alpha} = f_{\beta} = 1$ Hz. The parameter values $\phi_{\alpha}=\phi_{\beta} = 0.305$, $Q_{\alpha} = 0.988$ and $Q_{\beta} = -3$ are chosen such that $\alpha(0) = 1$ and $\beta(0) = 0$. The resulting image contrast and brightness evolution are shown in Figure 4(d). Simulations are first made with a uniform background. Then, in another simulated scenario, as exhibited in Figure 4(c), triangular patterns at gray level 145 are added into the background in each simulated image frame, to investigate the performance of the proposed ALISSE method to interfering non-uniform background. Some of these extra background patterns are very close to or even touch the objective bar region in the images, simulating a difficult situation for motion estimation.

## 4.2. Results under different simulation conditions

The results in 4 different cases will be presented below, starting from the basic case, then with one more difficulty added to each next case: light variation, non-uniform background, and random noise.

### The basic case: no light variation, uniform background, noise-free

The proposed ALISSE method is first tested on the simplest data set: the image sequence with uniform background, without light variation and without noise corruption. In this test, motion estimation is carried out twice, first with the selected 359 active pixels (shown in Figure 5(a)), then with the entire ROI, for the purpose of comparison. The estimated horizontal and vertical motion curves (from frames 3,000–3100 in the dataset) are presented in the first row of Figure 6. It is shown that, in both directions, the estimated
<figure><img src="image_6.png" /><figcaption>Figure 4: (a) Simulated image (101×101 pixels) and the region of interest (ROI, 59×59 pixels). The bright rectangular bar in the ROI is the image of the moving object. Its motion direction, shown by the yellow arrow, is along $2\pi/9$ (40°) with respect to the horizontal axis. (b) Image intensity profile simulated with a uniform background, along the horizontal line cutting the middle of the ROI. (c) Several static triangular background patterns are added to simulate an interfering background. (d) Evolution of the parameters $\alpha$ and $\beta$ of the simulated image sequences, as introduced in equation (9).</figcaption></figure>

motions based on active pixels and ROI all accurately match the true value. The corresponding RMSE and RE results presented in the first line of Table 1. The estimation errors with active pixel selection yield the RMSE of merely 0.006 pixel, whereas the errors with the whole ROI are about 3 times as large. These first results show the advantage of using active pixels over the whole ROI.

Regarding the computation time, 0.13 seconds are consumed for the active pixels selection, then 6 milliseconds are spent for other procedures of the off-line computation. The on-line motion estimation based on the active pixels takes 0.19 seconds in total, corresponding to 0.03 milliseconds per frame on average. In this sense, if dealing with one megapixel (1,000,000 pixels), the total on-line time by the proposed algorithm would be 85.9 milliseconds. For the motion estimation using the entire ROI, there is no active pixel selection. 7 milliseconds are used for off-line computation, and a total time of 0.49 seconds is consumed for on-line processing the video sequence, corresponding to 0.08 milliseconds per frame.

390

<figure><img src="image_7.png" /><figcaption>Figure 5: Active pixels selected within ROI in cases of (a): no light variation, uniform background, and no noise corruption, (b): with light variation, uniform background, and no noise corruption, (c): with light variation, non-uniform background, and no noise corruption, (d): with light variation, non-uniform background, and noise corruption (signal-to-noise ratio is 25.17 dB)</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>Figure 6: Estimated horizontal (left) and vertical (right) displacements by the proposed method from frames 3,000–3,100 in the simulated datasets under different conditions.</figcaption></figure>

<table><caption>Table 1: RMSE in pixel of motion estimation by ALISSE without noise corruption.</caption><thead><tr><td>Test conditions</td><td>Active pixels</td><td>ROI</td></tr></thead><tbody><tr><td>No light variation<br>Uniform background</td><td>0.006</td><td>0.015</td></tr><tr><td>Light variation<br>Uniform background</td><td>0.005</td><td>0.014</td></tr><tr><td>Light variation<br>Non-uniform background</td><td>0.010</td><td>0.061</td></tr></tbody></table>

### *Light variation*

The robustness of the proposed method to image brightness and contrast variations is assessed. The simulated periodic evolution of image intensity variations is described in equation (47). The other simulation parameters are the same as the previous test. In this case, 348 active pixels are selected as illustrated in Figure 5(b), which has no significant difference from the previous case shown in 5(a). In this test, the results are also compared between using active pixels and the whole ROI. The estimated motion curves are exhibited in the second row of Figure 6, which shows no significant difference from the previous case. As shown in the second line of Table 1, the RMSE is slightly smaller than in the previous test. Again the results using active pixels are more accurate than using the whole ROI. These results demonstrate the robustness
of the ALISSE method to brightness and contrast variations.

### *Light variation + non-uniform background*

In addition to illumination variations, in the next case, the robustness to non-uniform background of the proposed method is tested. As shown in Figure 5(c), in this test without noise corruption, 342 active pixels are selected within ROI, which are mostly similar to those of the previous cases, with some difference close to the static background patterns touching the moving object. The estimated motion curves under this condition are shown in the third row of Figure 6. In this case, the active pixel based estimated motion still matches with the true value accurately, whereas the ROI based result shows significant bias, especially in the vertical direction. This result is also indicated quantitatively by the estimation errors presented in the last row of Table 1. Due to the interfering non-uniform background, the error levels of estimation on active pixels and ROI are larger than in previous results. The result based on active pixels is much more accurate than those using the whole ROI, with a factor of about 6, at around 0.01 pixel.

### *Light variation + non-uniform background + random noise*

In the final test case, Gaussian noises are added to the simulated images at different SNR levels as defined in equation (44). The selected active pixels under the highest chosen noise level (SNR=25.17 dB) are presented in Figure 5(d). It shows no significant difference to the case without noise corruption, which illustrates the noise tolerance of the active pixel selection. The estimated motion curves under the severest noise condition (SNR=25.17 dB) are shown in the last row of Figure 6. The active pixel based estimation starts to show slight bias, and the ROI based result is greatly corrupted, which is similar to the previous test. According to this result, only the active pixel based estimation errors with respect to different noise levels are considered in Table 2. It can be observed that the RMSE increases from 0.010 pixel to 0.034 pixel with different noise levels, which seems reasonable given also the simultaneous perturbations by light variation and non-uniform background. These results further confirm the robustness of the proposed motion estimation method.

<table><caption>Table 2: RMSE in pixel unit of simulated motion estimation by the proposed ALISSE method on active pixels under different signal-to-noise ratios, with light variation and non-uniform background.</caption><thead><tr><th>No noise</th><th>SNR=37.21 dB</th><th>SNR=31.19 dB</th><th>SNR=27.66 dB</th><th>SNR=25.17 dB</th></tr></thead><tbody><tr><td>0.010</td><td>0.013</td><td>0.019</td><td>0.025</td><td>0.034</td></tr></tbody></table>
# 5. Experimental validation

In order to further evaluate the performance of the proposed motion estimation method, experiments have been conducted in the Structures and Integrated Instrumentation Laboratory at Gustave Eiffel University in Nantes, France. Two experimental conditions have been tested, first with no light variation and a uniform background, then with light variation and a non-uniform background.

## 5.1. *Experimental setup and evaluation criteria*

The experiments are carried out on a steel cantilever beam structure of size 30 mm × 2,000 mm × 5mm without any targets or markers to promote the detection of beam vibration. Material properties of the beam are presented in Table 3. A high speed 21MPx Mikroton EoGens camera is installed on a tripod at a distance of 10 m from the beam, as illustrated in Figure 7(a). The beam is installed on and driven by a shake table imposing a horizontal random excitation perpendicular to the camera axis, which is band-limited between 1 Hz and 200 Hz. The vibration of the beam is simultaneously measured by a camera and three laser displacement sensors (LDS), so that the displacements estimated from the video images can be compared with the laser measurements. The resolution of the recorded video frames is 1,632×128 pixels, at a sampling rate of 600 Hz. The duration of each measurement is 2 minutes, corresponding to 72,000 video frames. A scaling factor of 1.14 mm/pixel is measured by calculating the ratio between the image size and the physical size of the captured scene. Two Keyence IL-030 [62] and one Wenglor YP11MGVL80 [63] laser sensors are used to provide simultaneous measurement. Their characteristics are presented in Table 4, and their positions alongside the beam are shown in Figure 8(a). Three ROIs of size 20×10 pixels are selected in the video, referred to as “ROI 1”, “ROI 2” and “ROI 3” corresponding to the laser measurement positions, exhibited in Figure 8(b).

<table><caption>Table 3: Material properties of the tested cantilever beam.</caption><thead><tr><td>Properties</td><td>Value</td></tr></thead><tbody><tr><td>Young's modules (Pa)</td><td>2.1 × 10<sup>11</sup></td></tr><tr><td>Mass density (kg/m<sup>3</sup>)</td><td>7850</td></tr><tr><td>Poisson's ratio</td><td>0.21</td></tr></tbody></table>

The different experimental conditions have been implemented as follows. A blank curtain is placed behind the beam for the experiment with a uniform background (see Figure 7(b)), it is then replaced with a steel honeycomb grid structure for the experiment with a complex non-uniform background (see Figure 8(b)). This complex background contains mesh-shaped features in a laboratory environment, which is challenging for motion estimation algorithms. Three lamps of 1,500 W are utilized to provide light, and light variations have been implemented by modulating the illumination intensity of the lamps between 50%
<table><caption>Table 4: Parameters of laser sensors.</caption><thead><tr><td>Parameters</td><td>Keyence IL-030</td><td>Wenglor YP11MGVL80</td></tr></thead><tbody><tr><td>Reference distance</td><td>30 mm</td><td>50 mm</td></tr><tr><td>Measurement range</td><td>20 to 45 mm</td><td>50 to 100 mm</td></tr><tr><td>Wavelength of light source</td><td>655 nm</td><td>655 nm</td></tr><tr><td>Repeatability</td><td>1 µm</td><td>6 µm</td></tr><tr><td>Resolution</td><td>1 µm</td><td>25 µm</td></tr><tr><td>Light spot diameter</td><td>200×750 µm</td><td>1 mm</td></tr></tbody></table>

<figure><img src="image_9.png" /><figcaption>Figure 7: (a) Experimental setup illustration from the view along the camera axis, as well as a view of the beam from the camera. (b) Camera is placed 10 m away from the beam. (c) Background intensity variation in the experimental video image.</figcaption></figure>

and 100% periodically at 1 Hz, which is close to the first natural frequency of the beam. The temporal evolution of light illumination in the video, estimated from the image intensities in a small area on the static background, is shown in Figure 7(c). Experiments have been carried out with shake table settings yielding small vibration amplitudes (up to 2 pixels in ROI 1, and 1 pixel in ROI 2 and ROI 3) and large vibration amplitudes (up to 3.5 pixels in ROI 1, and 2 pixels in ROI 2 and ROI 3).

Our proposed method ALISSE is applied using the active pixel selection, and compared to two other relevant methods for motion estimation. Since the robustness of phase-based motion estimation has been investigated in previous works [35, 37], our method is compared to phase-based estimation using steerable filters (ST) [26] and the phase only correlation method (POC) [64]. The details of these two algorithms
<figure><img src="image_10.png" /><figcaption>Figure 8: (a) Cantilever beam and laser sensors, and (b) ROIs (20×10 pixels) corresponding to three laser sensor positions within video frames (1632×128 pixels) of uniform background (left) and non-uniform background (right).</figcaption></figure>

are described in Appendix C. Both ALISSE and ST apply active pixel selection, for which 6,000 frames of the video sequence (corresponding to around 10 cycles of beam vibration) have been used for efficiency, whereas no active pixel selection is available for POC. As a steerable filter-based method, ST assumes that the motion direction is known. For fair comparisons, the vibration direction is assumed to be known for all the evaluated algorithms, with $\theta = 0$ (see Section 2.5), i.e., only the horizontal displacement of the beam is estimated. The computations of RMSE are accordingly simplified as $RMSE_u$

$$RMSE_u = \sqrt{\frac{\sum_{t=1}^{N_s} [\hat{u}(t) - u_L(t)]^2}{N_s}}, \qquad (48)$$

<footnote n="470">470 where $\hat{u}(t)$ is the estimated horizontal displacement by each tested method, $u_L(t)$ is the LDS measurement of horizontal beam displacement.</footnote>
<table><caption>Table 5: RMSE<sub>u</sub> in pixel for each ROI evaluated under small vibration and large vibrations, both with constant lighting and uniform background, and the averaged time (in milliseconds per frame) for on-line computation on each video frame.</caption><thead><tr><th rowspan="2">Methods</th><th colspan="3">Small vibration</th><th colspan="3">Large vibration</th><th rowspan="2">On-line time ms/frame</th></tr><tr><th>ROI 1</th><th>ROI 2</th><th>ROI 3</th><th>ROI 1</th><th>ROI 2</th><th>ROI 3</th></tr></thead><tbody><tr><td>ST</td><td>0.028</td><td>0.022</td><td>0.023</td><td>0.108</td><td>0.065</td><td>0.030</td><td>6.41</td></tr><tr><td>POC</td><td>0.029</td><td>0.020</td><td>0.023</td><td>0.105</td><td>0.058</td><td>0.032</td><td>0.24</td></tr><tr><td>ALISSE</td><td>0.029</td><td>0.019</td><td>0.023</td><td>0.105</td><td>0.058</td><td>0.034</td><td>0.06</td></tr></tbody></table>

## 5.2. No light variation and uniform background

The beam vibration experiment is conducted without light variation and with a uniform background, first with the small vibration amplitude and then with the large vibration amplitude. For both cases, the RMSE<sub>u</sub> is reported in Table 5. In both cases of small and large vibration, the tested methods have very close performance, and the ST method behaves slightly worse in ROIs 1 and 2 but better in ROI 3 than the other two methods.

Besides the accuracy, the computational efficiency is also shown in Table 5, where the average time for on-line processing of each single video frame is reported. ALISSE requires 0.06 milliseconds per frame, corresponding to 192.2 milliseconds per megapixel, which is almost one fourth of the cost of POC and 1/100 cost of ST. In addition to the on-line processing of each video frame, each of the methods performs also some off-line computations for the whole video sequence. However, these computations are performed only once and the consumed time is usually negligible. For the proposed method ALISSE, the main cost of the preprocessing is the active pixel selection, requiring merely around 0.30 seconds. In total, less than 5 seconds are required for processing the selected three ROIs of size 20×10 pixels in the whole 72,000-frame video with ALISSE, which represents a 2 minutes vibration record. These results validate that the real-time implementation of vibration estimation can be realized by using the proposed method.

The proposed method ALISSE is designed with robustness to light variations and non-uniform background. While experimental results under such more realistic conditions are reported in the next subsection, the results presented above show that ALISSE also has a good performance under more ideal working conditions.

## 5.3. Light variation and non-uniform background corruptions

The robustness of displacement estimation to light variations and non-uniform background is investigated in this section, where performance evaluation is also carried out under both small and large vibrations. First, the displacement signals estimated by each of the tested methods are shown alongside the measured displacements by the laser sensor in ROI 1 for 200 frames in the upper part of Figure 9. The estimates
from the proposed method ALISSE seem to match the measured displacements the best, whereas there are larger discrepancies with ST and POC. This is also visible in the errors shown in the lower part of Figure 9, where ALISSE exhibits the smallest errors. The quantitative results of RMSE<sub>u</sub> over the entire video sequence are shown in Table 6. In both cases of small and large vibrations, the proposed method ALISSE clearly has the best performance, with errors being one third to a half lower than the ones of ST and POC. These results confirm the robustness of the proposed method to complex environmental conditions that it has been designed for. It should be noted that the contrast between the moving body and the background was particularly low in ROIs 1 and 3 (see Figure 8(b)), but errors with ALISSE were still low in these ROIs with 0.3 and 0.1 pixels, respectively. Because the left side of the background in ROI 1 is darker than the right side, the dark beam has a lower contrast with the background when it moves to the left side. For this reason the curves in Figure 9 show larger errors when the displacement is negative.

As in the previous case, the computation times are also shown Table 6, confirming the computational performance from the previous section. The average on-line computation time of 0.06 milliseconds per frame (corresponding to 252.1 milliseconds per megapixel) for three ROIs is fast enough for real time processing of the video sequence, which is sampled at 600 frames per second.

<table><caption>Table 6: RMSE<sub>u</sub> in pixel for each ROI evaluated under small and large vibrations, both with light variation and non-uniform background, and the averaged time (in milliseconds per frame) for on-line computation on each video frame.</caption><thead><tr><th rowspan="2">Methods</th><th colspan="3">Small vibration</th><th colspan="3">Large vibration</th><th rowspan="2">On-line time ms/frame</th></tr><tr><th>ROI 1</th><th>ROI 2</th><th>ROI 3</th><th>ROI 1</th><th>ROI 2</th><th>ROI 3</th></tr></thead><tbody><tr><td><strong>ST</strong></td><td>0.269</td><td>0.087</td><td>0.163</td><td>0.447</td><td>0.119</td><td>0.217</td><td>6.05</td></tr><tr><td><strong>POC</strong></td><td>0.267</td><td>0.046</td><td>0.138</td><td>0.425</td><td>0.074</td><td>0.382</td><td>0.24</td></tr><tr><td><strong>ALISSE</strong></td><td>0.205</td><td>0.030</td><td>0.051</td><td>0.296</td><td>0.063</td><td>0.104</td><td>0.06</td></tr></tbody></table>

## 6. Conclusions

A new real-time subpixel motion estimation method has been presented in this paper for video image-based vibration monitoring. The proposed method actively compensates for image brightness and contrast variations, while preventing interference from non-uniform background features.

The performance of the proposed method has been evaluated both with simulated video sequences and through a laboratory experiment on a cantilever beam. The following results were obtained:

* The simulation results show that the estimation error remain below 0.025 pixels, despite light variation, a non-uniform background, and high noise conditions (signal-to-noise ratio is 27.66 dB). In the
<figure><img src="image_11.png" /><figcaption>Figure 9: (a) Measured displacements with laser sensors (LDS) and estimated displacements from video frames 55900–56100 in ROI 1 with the three tested methods, for small (left) and large (right) vibration experiments with light variations and non-uniform background. (b) The corresponding errors with respect to laser measurements.</figcaption></figure>

laboratory experiments under complex environmental conditions, combining a 1 Hz ambient lighting variation and a non-uniform background, the errors with respect to laser sensor measurements range from 0.051 pixels to 0.296 pixels, depending on the contrast between the moving body and the non-uniform background at the measurement positions.

* Thanks to the efficient closed-form solutions of the optimization problems involved, the online computation for video sequence processing is fast enough for real-time applications. In the simulation example, the average computation time for processing each image frame containing a single ROI is 0.03 milliseconds. In the reported laboratory experiments, the on-line processing per frame containing 3 ROIs takes 0.06 milliseconds, while the video sampling rate of 600 Hz corresponds to 1.67 milliseconds per frame. These numerical computations have been conducted in Matlab on a laptop computer. There is no doubt that it is possible to process thousands of ROIs in real-time for large scale vibration
monitoring when parallelizing and optimizing the computations on adequate equipment.

The outcomes of this study provide substantial evidence that the proposed motion signal extraction method is promising for highly efficient and robust vibration measurement estimation in real-world applications.

# Appendix A. Some lemmas

**Lemma 1 (Least squares problem well-definedness).** Consider the linear regression

$$z_k = \psi_k^T \eta + e_k \quad (A.1)$$

with the regressor vector $\psi_k \in \mathbb{R}^n$, the regression coefficient vector $\eta \in \mathbb{R}^n$, the dependent variable $z_k \in \mathbb{R}$, the regression error $e_k \in \mathbb{R}$, and the observation index $k = 1, 2, \dots, N$. The least squares problem

$$\min_{\eta} (z_k - \psi_k^T \eta)^2 \quad (A.2)$$

has a unique solution, if and only if the matrix

$$\sum_{k=1}^{N} \psi_k \psi_k^T \quad (A.3)$$

is invertible.

This well known result is due to the fact that the sum expressed in (A.3) is the Hessian matrix (second-order derivative of the minimized criterion in $\eta$). This matrix is always symmetric positive semi-definite, and is positive definite when it is invertible.

**Lemma 2 (Well-definedness of the least squares problem with offset).** Consider the linear regression

$$z_k = \phi_k^T \eta + \beta + e_k \quad (A.4)$$

with the regressor vector $\phi_k \in \mathbb{R}^n$, the regression coefficient vector $\eta \in \mathbb{R}^n$, the dependent variable $z_k \in \mathbb{R}$, the offset $\beta \in \mathbb{R}$, the regression error $e_k \in \mathbb{R}$, and the observation index $k = 1, 2, \dots, N$. The least squares problem

$$\min_{\eta, \beta} (z_k - \phi_k^T \eta - \beta)^2 \quad (A.5)$$

has a unique solution, if and only if empirical covariance matrix of the regressors $\phi_1, \phi_2, \dots, \phi_N$ is invertible. This empirical covariance matrix is computed as

$$\frac{1}{N-1} \sum_{k=1}^{N} (\phi_k - \bar{\phi})(\phi_k - \bar{\phi})^T \quad (A.6)$$

with

$$\bar{\phi} \triangleq \frac{1}{N} \sum_{k=1}^{N} \phi_k. \quad (A.7)$$
## Proof.

First introduce some matrix and vector notations:

$$ \Phi \triangleq \left[ \begin{matrix} \phi_1^T \\ \phi_2^T \\ \vdots \\ \phi_N^T \end{matrix} \right], \quad \vec{1} = \left[ \begin{matrix} 1 \\ 1 \\ \vdots \\ 1 \end{matrix} \right] \right\} \text{N rows.} \qquad (A.8) $$

To rewrite the linear regression with offset (A.4) in the form of (A.1) so that Lemma 1 can be applied, let

$$ \psi_k = \begin{bmatrix} \phi_k \\ 1 \end{bmatrix}, \quad \Psi = \begin{bmatrix} \psi_1^T \\ \psi_2^T \\ \vdots \\ \psi_N^T \end{bmatrix} = \begin{bmatrix} \phi_1^T & 1 \\ \phi_2^T & 1 \\ \vdots & \vdots \\ \phi_N^T & 1 \end{bmatrix} = [\Phi \quad \vec{1}], \qquad (A.9) $$

then the linear regression with offset (A.4) is rewritten in the form of (A.1), and

$$ \sum_{k=1}^{N} \psi_k \psi_k^T = \Psi^T \Psi. \qquad (A.10) $$

According to Lemma 1, the considered least squares problem has a unique solution if and only if the matrix $\Psi^T \Psi$ is invertible, as in (A.10).

With $\bar{\phi}$ as in (A.7), define

$$ \tilde{\Phi} \triangleq \begin{bmatrix} \phi_1^T - \bar{\phi}^T \\ \phi_2^T - \bar{\phi}^T \\ \vdots \\ \phi_N^T - \bar{\phi}^T \end{bmatrix} = \Phi - \vec{1} \bar{\phi}^T. \qquad (A.11) $$

For $\Psi$ built in (A.9), $\Psi = [\Phi \quad \vec{1}]$, then some simple matrix computations yield

$$ \Psi \begin{bmatrix} I & 0 \\ -\bar{\phi}^T & 1 \end{bmatrix} = [\Phi \quad \vec{1}] \begin{bmatrix} I & 0 \\ -\bar{\phi}^T & 1 \end{bmatrix} = [\Phi - \vec{1} \bar{\phi}^T \quad \vec{1}] = [\tilde{\Phi} \quad \vec{1}]. \qquad (A.12) $$

In these last equalities, the block-triangular matrix is in fact triangular, with the diagonal filled with 1's. This square matrix is thus invertible, then the matrix $\Psi$ has the same column rank as $[\tilde{\Phi} \quad \vec{1}]$.

560

The invertibility of $\Psi^T \Psi$ is equivalent to the full column rankness of $\Psi$, which is in turn equivalent to the full column rankness of $[\tilde{\Phi} \quad \vec{1}]$.

To summarize, with the shorthand FCR=“Full Column Rank”, it is already proved above that

$$ \text{FCR} \left[ \tilde{\Phi} \quad \vec{1} \right] \Leftrightarrow \text{FCR} \Psi \Leftrightarrow \text{Invertible} \ \Psi^T \Psi = \sum_{k=1}^{N} \psi_k \psi_k^T \qquad (A.13) $$
On the other hand, in the condition of Lemma 2, the matrix sum (A.6)

$$ \frac{1}{N-1} \sum_{k=1}^{N} (\phi_k - \bar{\phi})(\phi_k - \bar{\phi})^T = \frac{1}{N-1} \tilde{\Phi}^T \tilde{\Phi}. \quad (A.14) $$

To apply Lemma 1 requiring an invertible matrix in (A.3), which is reproduced at the end of (A.13), it remains to show that

$$ \text{Invertible } \tilde{\Phi}^T \tilde{\Phi} \Leftrightarrow \text{FCR } [\tilde{\Phi} \quad \vec{1}]. \quad (A.15) $$

Let us first prove the sufficiency, i.e. “⇒” in (A.15).

With $\bar{\phi}$ as in (A.7), $\tilde{\Phi}$ as in (A.11),

$$ \vec{1}^T \tilde{\Phi} = \sum_{k=1}^{N} \phi_k^T - N \bar{\phi} = 0. \quad (A.16) $$

Suppose that a vector

$$ \begin{bmatrix} w \\ v \end{bmatrix} \in \mathbb{R}^{n+1}, \quad (A.17) $$

570

with $W \in \mathbb{R}^n$ and $v \in \mathbb{R}$, satisfies

$$ [\tilde{\Phi} \quad \vec{1}] \begin{bmatrix} w \\ v \end{bmatrix} = 0. \quad (A.18) $$

Then

$$ \vec{1}^T [\tilde{\Phi} \quad \vec{1}] \begin{bmatrix} w \\ v \end{bmatrix} = 0, \quad (A.19) $$

hence

$$ 0 = \vec{1}^T [\tilde{\Phi} \quad \vec{1}] \begin{bmatrix} w \\ v \end{bmatrix} = [0 \quad N] \begin{bmatrix} w \\ v \end{bmatrix} = Nv, \quad (A.20) $$

therefore, $v = 0$, and it follows from (A.18) that

$$ 0 = [\tilde{\Phi} \quad \vec{1}] \begin{bmatrix} w \\ v \end{bmatrix} = \tilde{\Phi}w. \quad (A.21) $$

The zero vector $w = 0$ is the unique solution of $\tilde{\Phi}w = 0$ if and only if $\tilde{\Phi}$ has a full column rank.

Because (A.18) implies $v = 0$, it follows from the equality in (A.21) that $[w^T \ v^T] = 0$ is the unique solution (A.18) if and only if $w = 0$ is the unique solution of $\tilde{\Phi}w = 0$. It is then proved that

$$ \text{FCR } [\tilde{\Phi} \quad \vec{1}] \Leftrightarrow \text{FCR } \tilde{\Phi} \Leftrightarrow \text{Invertible } \tilde{\Phi}^T \tilde{\Phi}. \quad (A.22) $$

Hence (A.15) is confirmed. Then the proof of Lemma 2 is completed. □
**Lemma 3.** Given two sets of vectors $\phi_k \in \mathbb{R}^n$ and $\omega_k \in \mathbb{R}^m$, both indexed by $k = 1, 2, \dots, N$, related by

$$ \phi_k = M\omega_k \qquad (A.23) $$

for all $k = 1, 2, \dots, N$ and for some full row rank matrix $M \in \mathbb{R}^{n \times m}$. If $\omega_k$ has an invertible empirical covariance, then so does $\phi_k$.

**Proof.** The empirical covariance of $\phi_k$ is computed through (A.6) and (A.7). With the matrix notations defined in (A.8) and (A.11),

$$ \tilde{\Phi} = \Phi - \tilde{1} \bar{\phi}^T, \qquad (A.24) $$

then the empirical covariance matrix of $\phi_k$ is given by

$$ \frac{1}{N-1} \tilde{\Phi}^T \tilde{\Phi} \qquad (A.25) $$

Similarly, $\Omega, \tilde{\Omega}$ and $\bar{\omega}$ are defined from $\omega_k$, and

$$ \tilde{\Omega} = \Omega - \tilde{1} \bar{\omega}^T. \qquad (A.26) $$

The equality (A.23) for all $k = 1, 2, \dots, N$ implies

$$ \Phi^T = M\Omega^T \qquad (A.27) $$

$$ \bar{\phi} = M\bar{\omega}. \qquad (A.28) $$

Therefore,

$$ \tilde{\Phi}^T = \Phi^T - \bar{\phi} \tilde{1}^T \qquad (A.29) $$

$$ = M\Omega^T - M\bar{\omega} \tilde{1}^T \qquad (A.30) $$

$$ = M\tilde{\Omega}^T, \qquad (A.31) $$

then

$$ \frac{1}{N-1} \tilde{\Phi}^T \tilde{\Phi} = M \left( \frac{1}{N-1} \tilde{\Omega}^T \tilde{\Omega} \right) M^T, \qquad (A.32) $$

which means

$$ \mathrm{Cov}(\phi_k) = M \mathrm{Cov}(\omega_k) M^T \qquad (A.33) $$

with “Cov” denoting the empirical covariance.

Because $\mathrm{Cov}(\omega_k)$ is invertible (assumed condition of this lemma), it is also symmetric positive definite.
Moreover, $M$ has full row rank, then $\mathrm{Cov}(\phi_k)$ is also positive definite, thus invertible. $\square$
# Appendix B: Technical proofs

## Proof of Proposition 2.

To prove this result with Lemma 2 (see Appendix A), the considered least squares problem (12) is associated to a linear regression as in (A.4), i.e.,

$$z_k = \phi_k^T \eta + \beta + e_k \qquad (B.1)$$

in which $\phi_k^T \eta + \beta$ corresponds to the terms $a_{i,j}x^a + b_{i,j}y^a + c_{i,j}\alpha + \beta$ in (12). When the linear interpolation (6) was built, the linear coefficients $a_{i,j}$, $b_{i,j}$, $c_{i,j}$ were fitted to the image intensities of the pixel $(i, j)$ and its 4 surrounding pixels by solving (7). In fact, it was made for each active pixel indexed by $(i, j) \in A$. To convert the double index $(i, j)$ to a single index $k$, in agreement with (B.1), choose an arbitrary order of the elements of $A$ so that they are indexed by $(i_k, j_k)$ as

$$A = \{(i_1, j_1), (i_2, j_2), \dots, (i_N, j_N)\} \qquad (B.2)$$

with $N$ the number of active pixels. Then in the considered least squares problem (12), the linear interpolation part $a_{i,j}x^a + b_{i,j}y^a + c_{i,j}\alpha$ becomes

$$a_{i_k, j_k} x^a + b_{i_k, j_k} y^a + c_{i_k, j_k} \alpha = \begin{bmatrix} a_{i_k, j_k} & b_{i_k, j_k} & c_{i_k, j_k} \end{bmatrix} \begin{bmatrix} x^a \\ y^a \\ \alpha \end{bmatrix} \qquad (B.3)$$

$$= \phi^T \eta \qquad (B.4)$$

with

$$\phi_k = \begin{bmatrix} a_{i_k, j_k} \\ b_{i_k, j_k} \\ c_{i_k, j_k} \end{bmatrix}, \quad \eta = \begin{bmatrix} x^a \\ y^a \\ \alpha \end{bmatrix}. \qquad (B.5)$$

corresponding to the notations in (B.1).

To apply Lemma 2, we need to show that its assumption, i.e., the invertibility of the covariance matrix expressed in (A.6) and (A.7) with $\phi_k$ specified in (B.5), is ensured by the assumption made in Proposition 2, i.e., the empirical covariance of the vectors $\zeta(i, j, t_1)$ defined in (39) is invertible.

The coefficients $a_{i_k, j_k}$, $b_{i_k, j_k}$, $c_{i_k, j_k}$ were computed in (40). Let us rewrite it with the single index $k$:

$$\phi_k = \begin{bmatrix} a_{i_k, j} \\ b_{i_k, j} \\ c_{i_k, j} \end{bmatrix} = (\Gamma^T \Gamma)^{-1} \Gamma^T \zeta(i_k, j_k, t_1) \qquad (B.6)$$
With $\Gamma$ as in (34),

$$ (\Gamma^T \Gamma)^{-1} \Gamma^T = \begin{bmatrix} 0 & -1/2 & 1/2 & 0 & 0 \\ 0 & 0 & 0 & -1/2 & 1/2 \\ 1/5 & 1/5 & 1/5 & 1/5 & 1/5 \end{bmatrix}. \qquad (B.7) $$

This matrix has a full row rank (rank=3). Based on equality (B.6), and according to Lemma 3, the invertibility of the empirical covariance matrix of $\zeta(i_k, j_k, t_1)$ (assumed in Proposition 2) implies the invertibility of the empirical covariance matrix of $\phi_k$. It then follows from Lemma 2 that the considered least squares problem (12) has a unique solution. □

### Proof of Proposition 3

Consider the criteria

$$ \rho(x, y, \alpha, \beta) = \sum_{(i,j) \in \mathcal{A}} [I_s(i, j, t_2) - \alpha(a_{i,j}x + b_{i,j}y + c_{i,j}) - \beta]^2, \qquad (B.8) $$

and

$$ \rho^a(x^a, y^a, \alpha, \beta) = \sum_{(i,j) \in \mathcal{A}} [I_s(i, j, t_2) - (a_{i,j}x^a + b_{i,j}y^a + c_{i,j}\alpha) - \beta]^2, \qquad (B.9) $$

which are respectively minimized in (10) and (12).

For any $(x, y, \alpha, \beta) \in \mathbb{R}^4$, not only those minimizing (10), the auxiliary variables $x^a$ and $y^a$ defined in (11) ensure that

$$ \rho(x, y, \alpha, \beta) = \rho^a(x^a, y^a, \alpha, \beta). \qquad (B.10) $$

Therefore, any value reached by $\rho(x, y, \alpha, \beta)$ is also reachable by $\rho^a(x^a, y^a, \alpha, \beta)$, then, for any $(x, y, \alpha, \beta) \in \mathbb{R}^4$,

$$ \rho(x, y, \alpha, \beta) \geq \underline{\rho} \triangleq \min_{(x^a, y^a, \alpha, \beta) \in \mathbb{R}^4} \rho^a(x^a, y^a, \alpha, \beta). \qquad (B.11) $$

If a quadruple $(x^*, y^*, \alpha^*, \beta^*)$ leads to

$$ \rho(x^*, y^*, \alpha^*, \beta^*) = \underline{\rho}, \qquad (B.12) $$

no value of $\rho(x, y, \alpha, \beta)$ can be smaller than $\underline{\rho}$, then $(x^*, y^*, \alpha^*, \beta^*)$ minimizes $\rho(x, y, \alpha, \beta)$.

### The case $\hat{\alpha} \neq 0$ in the solution of (12)

Let $\hat{x}^a, \hat{y}^a, \hat{\alpha}, \hat{\beta}$ denote the solution minimizing (12), and let $\hat{x}$ and $\hat{y}$ be computed through (17). Then, as a particular case of (B.10),

$$ \rho(\hat{x}, \hat{y}, \hat{\alpha}, \hat{\beta}) = \rho^a(\hat{x}^a, \hat{y}^a, \hat{\alpha}, \hat{\beta}) = \min \rho^a(x^a, y^a, \alpha, \beta) = \underline{\rho}, \qquad (B.13) $$

hence $(\hat{x}, \hat{y}, \hat{\alpha}, \hat{\beta})$ minimizes $\rho(x, y, \alpha, \beta)$, according to the remark immediately following (B.11).
Still in the case with $\hat{\alpha} \neq 0$, to show the uniqueness of $(\hat{x}, \hat{y}, \hat{\alpha}, \hat{\beta})$ minimizing $\rho(x, y, \alpha, \beta)$, suppose that another quadruple $(\hat{x}', \hat{y}', \hat{\alpha}', \hat{\beta}')$ minimizes also $\rho(x, y, \alpha, \beta)$. Then

$$ \rho(\hat{x}', \hat{y}', \hat{\alpha}', \hat{\beta}') = \underline{\rho} = \min \rho^a(x^a, y^a, \alpha, \beta). \qquad (B.14) $$

Let

$$ \hat{x}'^a = \hat{\alpha}' \hat{x}', \quad \hat{y}'^a = \hat{\alpha}' \hat{y}'. \qquad (B.15) $$

Because (B.10) holds for any $(x, y, \alpha, \beta) \in \mathbb{R}^4$ with the corresponding auxiliary variables, hence

$$ \rho(\hat{x}', \hat{y}', \hat{\alpha}', \hat{\beta}') = \rho^a(\hat{x}'^a, \hat{y}'^a, \hat{\alpha}', \hat{\beta}') = \min \rho^a(x^a, y^a, \alpha, \beta). \qquad (B.16) $$

Due to Proposition 2, the minimum of $\rho^a(x^a, y^a, \alpha, \beta)$ is unique, hence

$$ (\hat{x}'^a, \hat{y}'^a, \hat{\alpha}', \hat{\beta}') = (\hat{x}^a, \hat{y}^a, \hat{\alpha}, \hat{\beta}) \qquad (B.17) $$

and

$$ \hat{\alpha}' = \hat{\alpha} \neq 0. \qquad (B.18) $$

These results then imply, through (11) and (B.15), that

$$ (\hat{x}', \hat{y}', \hat{\alpha}', \hat{\beta}') = (\hat{x}, \hat{y}, \hat{\alpha}, \hat{\beta}). \qquad (B.19) $$

It is then proved that $(\hat{x}, \hat{y}, \hat{\alpha}, \hat{\beta})$ minimizing $\rho(x, y, \alpha, \beta)$ is unique when $\hat{\alpha} \neq 0$.

### The case $\hat{\alpha} = 0$ in the solution of (12)

When $\hat{\alpha} = 0$, it is not possible to divide $\hat{x}^a$ and $\hat{y}^a$ and by $\hat{\alpha}$, then the solution of (12) cannot be used to solve the nonlinear optimization problem (10).

Let $(x^*, y^*, \alpha^*, \beta^*)$ denote a solution solving the nonlinear optimization problem (10), i.e.,

$$ \rho(x^*, y^*, \alpha^*, \beta^*) = \min \rho(x, y, \alpha, \beta). \qquad (B.20) $$

In order to prove that $\alpha^* = 0$ holds certainly in this case (when $\hat{\alpha} = 0$), by means of contradiction, suppose that $\alpha^* \neq 0$.

Let

$$ x^{*a} = \alpha^{*} x^{*} \qquad (B.21) $$

$$ y^{*a} = \alpha^{*} y^{*}. \qquad (B.22) $$

Because $\hat{\alpha} = 0$ in the present case, $\rho^a(x^a, y^a, \alpha, \beta)$ is minimized by $(\hat{x}^a, \hat{y}^a, 0, \hat{\beta})$, which is the unique minimizing solution according to Proposition 2. Moreover, $\alpha^* \neq 0$ has been assumed (for proof by contradiction), then $(x^{*a}, y^{*a}, \alpha^*, \beta^*) \neq (\hat{x}^a, \hat{y}^a, 0, \hat{\beta})$. Then the uniqueness of the minimizing solution implies

$$ \rho^a(x^{*a}, y^{*a}, \alpha^*, \beta^*) > \rho^a(\hat{x}^a, \hat{y}^a, 0, \hat{\beta}). \qquad (B.23) $$
The standard least squares criterion $\rho^a(x^a, y^a, \alpha, \beta)$ is strictly convex under the assumption made in the present proposition. Define

$$ \bar{x}^a = \frac{1}{2} x^{*a} + \frac{1}{2} \hat{x}^a \qquad (B.24a) $$

$$ \bar{y}^a = \frac{1}{2} y^{*a} + \frac{1}{2} \hat{y}^a \qquad (B.24b) $$

$$ \bar{\alpha}^a = \frac{1}{2} \alpha^* + \frac{1}{2} \hat{\alpha} = \frac{1}{2} \alpha^* \neq 0 \qquad (B.24c) $$

$$ \bar{\beta}^a = \frac{1}{2} \beta^* + \frac{1}{2} \hat{\beta}, \qquad (B.24d) $$

then the strict convexity implies

$$ \rho^a(\bar{x}^a, \bar{y}^a, \bar{\alpha}, \bar{\beta}) < \frac{1}{2} \rho^a(x^{*a}, y^{*a}, \alpha^*, \beta^*) + \frac{1}{2} \rho^a(\hat{x}^a, \hat{y}^a, 0, \hat{\beta}). \qquad (B.25) $$

This result, together with (B.23), then leads to

$$ \rho^a(\bar{x}^a, \bar{y}^a, \bar{\alpha}, \bar{\beta}) < \rho^a(x^{*a}, y^{*a}, \alpha^*, \beta^*). \qquad (B.26) $$

650

With $\bar{\alpha} \neq 0$ (according to (B.24c)), let

$$ \bar{x} = \frac{\bar{x}^a}{\bar{\alpha}} \qquad (B.27) $$

$$ \bar{y} = \frac{\bar{y}^a}{\bar{\alpha}}. \qquad (B.28) $$

Because (B.10) holds for any $(x, y, \alpha, \beta) \in \mathbb{R}^4$ with the corresponding auxiliary variables,

$$ \rho(\bar{x}, \bar{y}, \bar{\alpha}, \bar{\beta}) = \rho^a(\bar{x}^a, \bar{y}^a, \bar{\alpha}, \bar{\beta}) \qquad (B.29) $$

$$ \rho(x^*, y^*, \alpha^*, \beta^*) = \rho^a(x^{*a}, y^{*a}, \alpha^*, \beta^*). \qquad (B.30) $$

It then follows from (B.26) that

$$ \rho(\bar{x}, \bar{y}, \bar{\alpha}, \bar{\beta}) < \rho(x^*, y^*, \alpha^*, \beta^*). \qquad (B.31) $$

This inequality is in contradiction with the assumed $(x^*, y^*, \alpha^*, \beta^*)$ minimizing $\rho(x, y, \alpha, \beta)$ with $\alpha^* \neq 0$. Therefore, $\alpha^* = 0$ must hold when $\hat{\alpha} = 0$.

When $\alpha = 0$, (B.8) becomes

$$ \rho(x, y, 0, \beta) = \sum_{(i,j) \in \mathcal{A}} [I_s(i, j, t_2) - \beta]^2. \qquad (B.32) $$

Therefore, $x$ and $y$ have no effect on $\rho(x, y, 0, \beta)$, and the minimizing $\beta$ is given by (43), which is simply equal to the average value of $I_s(i, j, t_2)$ for $(i, j) \in \mathcal{A}$.
Appendix C. Tested phased based methods on the experimental data

This appendix briefly presents the two alternative motion estimation methods which have been used as a reference to compare with the proposed ALISSE method. The two methods are phase-based and thus are known to be robust against lighting variations at some extent [26, 27].

C1. Steerable Filter (ST): The steerable filter has been introduced for different image processing purposes in [30]. It has been applied to motion estimation in [31]. The original subpixel motion range limitation of the ST method was mitigated in [25] by a multiscale pyramidal approach, and in [26] by a step-wise cumulative approach.

The steerable filter allows to express the local spatial grey level variations in terms of magnitude $A^{\Theta}(i, j, t)$ and phase $\varphi^{\Theta}(i, j, t)$ as follows:

$$ \mathcal{I}(i, j, t) = \mathcal{S}^{\Theta}(i, j, t) \otimes I(i, j, t) = A^{\Theta}(i, j, t) e^{i\varphi^{\Theta}(i, j, t)}, \quad (C.1) $$

where $\mathcal{S}$ denotes the steerable filter, and $\Theta$ represents the direction angle of the steerable filter with respect to the horizontal axis on image, $\mathcal{I}$ is the filtered image, $A$ and $\varphi$ are the amplitude and phase of the filtered image, $(i, j)$ is the spatial coordinate, and $t$ denotes time instant.

For motion estimation, the ST is applied on each image frame. The local phase is then related to the motion according to the assumption of constant contour of local phase [29], the velocity ($\varpi$ in horizontal $v$ in vertical directions) of selected pixels can be calculated with the phase difference as follows:

$$ \varpi(i, j, t) = - \left[ \frac{\partial \varphi^0(i, j, t)}{\partial x} \right]^{-1} \frac{\partial \varphi^0(i, j, t)}{\partial t}, \quad v(i, j, t) = - \left[ \frac{\partial \varphi^{\pi/2}(i, j, t)}{\partial y} \right]^{-1} \frac{\partial \varphi^{\pi/2}(i, j, t)}{\partial t}, \quad (C.2) $$

and the subpixel displacement ($d_x$ in horizontal and $d_y$ in vertical direction) between frames $t_1$ and $t_2$ is then determined by integration over time $t$:

$$ \begin{aligned} x(i, j, t_2) &= \left[ \frac{\partial \varphi^0(i, j, t_2)}{\partial x} \right]^{-1} (\varphi^0(i, j, t_2) - \varphi^0(i, j, t_1)), \\ y(i, j, t_2) &= \left[ \frac{\partial \varphi^{\pi/2}(i, j, t_2)}{\partial y} \right]^{-1} (\varphi^{\pi/2}(i, j, t_2) - \varphi^{\pi/2}(i, j, t_1)). \end{aligned} \quad (C.3) $$

In [26], the basic ST approach is extended to motion amplitude larger than one pixel at a single image scale. To this aim, the set of active pixels is shifted between each frame according to the following condition:

$$ |\Delta x| + |\Delta y| \ge 1, \quad (C.4) $$

where $\Delta x$ (resp. $\Delta y$) are the horizontal (resp. vertical) pixel shift computed from the zero-mean normalized cross correlation, namely, ZNCC, between the ROI on the reference frame ($t_1$) and the one on the sequent frame at $t_2$. When this condition is true, the selected active pixels are translated with the updated pixel level shift ($\Delta x, \Delta y$), and the subpixel displacement is then estimated using Equation (C.3). Otherwise, only
the subpixel displacement in Equation (C.3) is considered. In the first case, the overall estimated motion composes of pixel level and subpixel displacements is:

$$ 
\begin{aligned}
d_x(i, j, t) &= \Delta x(t) + x(i, j, t), \\
d_y(i, j, t) &= \Delta y(t) + y(i, j, t).
\end{aligned} 
\qquad (C.5) $$

Finally, the displacement is averaged over the selected active pixels.

**C2. Phase Only Correlation (POC)**: Phase correlation methods are based on the Fourier shift theorem. Assuming the current image at time $t_2$ is the translated copy of the reference image frame $I(i, j, t_1)$ as follows:

$$ I(i, j, t_2) = I(i + d_x, j + d_y, t_1), \qquad (C.6) $$

the Fourier transforms of both images $\mathcal{F}(f, g, t)$ are related to each other according to:

$$ \mathcal{F}(f, g, t_2) = \mathcal{F}(f, g, t_1)e^{-i \cdot 2\pi (fd_x+gd_y)}, \qquad (C.7) $$

where $f$ and $g$ are spectral coordinates after transform, while $dx$ and $dy$ are the spatial translation. To directly extract the translation $d_x$ and $d_y$, the following normalized Fourier product between frames is computed:

$$ \mathcal{N}(f, g) = \frac{\mathcal{F}(f, g, t_2)\mathcal{F}(f, g, t_1)^*}{|\mathcal{F}(f, g, t_2)||\mathcal{F}(f, g, t_1)|} = e^{-i \cdot 2\pi (fd_x+gd_y)}, \qquad (C.8) $$

Then, the translation between both images is the solution of the following least squares plane fitting estimation:

$$ \min_{d_x, d_y} \sum_{f,g \in \text{ROI}} [\arg\{\mathcal{N}(f,g)\} + 2\pi (fd_x + gd_y)]^2. \qquad (C.9) $$

In practice, the estimation above is only valid for subpixel translation due to the phase wrapping [65], and the normalized cross power spectrum estimated on digital image is not ideally equivalent to Equation (C.8) owing to aliasing problem [27]. Thus, prior operations are further developed in POC algorithms, including estimating pixel level displacement and mitigating aliasing effect by choosing only low frequency components on images. In addition, some spectral magnitude weighting strategies are also performed in the POC algorithms to mitigate the influence of noise on small magnitudes spectral components [27, 64]. Performance comparison between the most well-known existing spectral-based and spatial-based POC techniques is available in [3].

## References

[1] X. Wang, M. Jaiswal, A. M. Corpuz, S. Paudel, A. Balu, A. Krishnamurthy, J. Yan, M.-C. Hsu, Photogrammetry-based computational fluid dynamics, Computer Methods in Applied Mechanics and Engineering 417 (2023) 116311.

[2] A. Alba, J. F. Vigueras-Gomez, E. R. Arce-Santana, R. M. Aguilar-Ponce, Phase correlation with sub-pixel accuracy: a comparative study in 1D and 2D, Computer Vision and Image Understanding 137 (2015) 76-87.
710

[3] X. Tong, Z. Ye, Y. Xu, S. Gao, H. Xie, Q. Du, S. Liu, X. Xu, S. Liu, K. Luan, U. Stilla, Image registration with Fourier-based image correlation: A comprehensive review of developments and applications, IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing 12 (10) (2019) 4062-4081.

[4] S. A. Mohamed, A. K. Helmi, M. A. Fkirin, S. M. Badwai, Subpixel accuracy analysis of phase correlation shift measurement methods applied to satellite imagery, International Journal of Advanced Computer Science and Applications 3 (12) (2012).

[5] J. Baqersad, P. Poozesh, C. Niezrecki, P. Avitabile, Photogrammetry and optical methods in structural dynamics – a review, Mechanical Systems and Signal Processing 86 (2017) 17–34.

[6] D. Feng, M. Q. Feng, Computer vision for SHM of civil infrastructure: From dynamic response measurement to damage detection – a review, Engineering Structures 156 (2018) 105–117.

[7] B. F. Spencer, V. Hoskere, Y. Narazaki, Advances in computer vision-based civil infrastructure inspection and monitoring, Engineering 5 (2) (2019) 199–222.

[8] Y. Fu, Y. Shang, W. Hu, B. Li, Q. Yu, Non-contact optical dynamic measurements at different ranges: a review, Acta Mechanica Sinica 37 (2021) 537–553.

720

[9] C. Ferraris, G. Amprimo, G. Pettiti, Computer vision and image processing in structural health monitoring: overview of recent applications, Signals 4 (3) (2023) 539–574.

[10] G. Besnard, F. Hild, S. Roux, “Finite-element” displacement fields analysis from digital images: application to Portevin–Le Châtelier bands, Experimental Mechanics 46 (2006) 789–803.

[11] B. Pan, Reliability-guided digital image correlation for image deformation measurement, Applied Optics 48 (8) (2009) 1535–1542.

[12] B. Pan, D. Wu, Y. Xia, An active imaging digital image correlation method for deformation measurement insensitive to ambient light, Optics & Laser Technology 44 (1) (2012) 204–209.

[13] V. Belloni, R. Ravanelli, A. Nascetti, M. Di Rita, D. Mattei, M. Crespi, py2DIC: a new free and open source software for displacement and strain measurements in the field of experimental mechanics, Sensors 19 (18) (2019) 3832.

730

[14] Y. Huang, J. Ji, K.-M. Lee, Model-based digital image correlation for noncontact deformation measurement of strain field and mechanical property, IEEE Transactions on Industrial Informatics 15 (9) (2019) 5109–5118.

[15] V. F. Sciuti, R. B. Canto, J. Neggers, F. Hild, On the benefits of correcting brightness and contrast in global digital image correlation: monitoring cracks during curing and drying of a refractory castable, Optics and Lasers in Engineering 136 (2021) 106316.

[16] Y. Zhou, Q. Zuo, L. Zhou, B. Yang, Z. Liu, Y. Liu, L. Tang, S. Dong, Z. Jiang, Image feature based quality assessment of speckle patterns for digital image correlation measurement, Measurement 222 (2023) 113590.

[17] M. Guizar-Sicairos, S. T. Thurman, J. R. Fienup, Efficient subpixel image registration algorithms, Optics Letters 33 (2) (2008) 156–158.

740

[18] D. Zhang, J. Guo, X. Lei, C. Zhu, A high-speed vision-based sensor for dynamic vibration analysis using fast motion extraction algorithms, Sensors 16 (4) (2016) 572.

[19] B. Xiong, Q. Zhang, V. Baltazart, On quadratic interpolation of image cross-correlation for subpixel motion extraction, Sensors 22 (3) (2022) 1274.

[20] X. W. Ye, T.-H. Yi, C. Z. Dong, T. Liu, Vision-based structural displacement measurement: system performance evaluation and influence factor analysis, Measurement 88 (2016) 372–384.

[21] R. Janeliukstis, X. Chen, Review of digital image correlation application to large-scale composite structure testing, Composite Structures 271 (2021) 114143.

[22] P. Frankovský, I. Delyová, P. Sivák, J. Bocko, J. Živčák, M. Kicko, Modal analysis using digital image correlation technique, Materials 15 (16) (2022) 5658.
750

[23] E. Cai, Y. Zhang, X. Lu, P. Li, T. Zhao, G. Lin, W. Guo, A target-free video structural motion estimation method based on multi-path optimization, *Mechanical Systems and Signal Processing* 198 (2023) 110452.

[24] Y. Wang, Z. Huang, P. Zhu, R. Zhu, T. Hu, D. Zhang, D. Jiang, Effects of compressed speckle image on digital image correlation for vibration measurement, *Measurement* 217 (2023) 113041.

[25] N. Wadhwa, M. Rubinstein, F. Durand, W. T. Freeman, Phase-based video motion processing, *ACM Transactions on Graphics* (2013).

[26] B. Merainani, B. Xiong, V. Baltazart, M. Döhler, J. Dumoulin, Q. Zhang, Subspace-based modal identification and uncertainty quantification from video image flows, *Journal of Sound and Vibration* 569 (2024) 117957.

[27] H. Stone, M. Orchard, E.-C. Chang, S. Martucci, A fast direct Fourier-based algorithm for subpixel registration of images, *IEEE Transactions on Geoscience and Remote Sensing* 39 (10) (2001) 2235-2243.

760

[28] K. Takita, T. Aoki, Y. Sasaki, T. Higuchi, K. Kobayashi, High-accuracy subpixel image registration based on phase-only correlation, *IEICE Transactions On Fundamentals Of Electronics, Communications and Computer Sciences* 86-A (2003) 1925-1934.

[29] D. J. Fleet, A. D. Jepson, Computation of component image velocity from local phase information, *International Journal of Computer Vision* 5 (1) (1990) 77-104.

[30] W. Freeman, E. Adelson, The design and use of steerable filters, *IEEE Transactions on Pattern Analysis and Machine Intelligence* 13 (9) (1991) 891-906.

[31] J. G. Chen, A. Davis, N. Wadhwa, F. Durand, W. T. Freeman, O. Büyüköztürk, Video camera-based vibration measurement for civil infrastructure applications, *Journal of Infrastructure Systems* 23 (3) (2017) B4016013.

770

[32] Y. Miao, Y. Kong, H. Nam, S. Lee, G. Park, Phase-based vibration imaging for structural dynamics applications: Marker-free full-field displacement measurements with confidence measures, *Mechanical Systems and Signal Processing* 198 (2023) 110418.

[33] M. Li, G. Liu, Z. Mao, Phased-based motion estimation through short-distance Hilbert transform, *Mechanical Systems and Signal Processing* 211 (2024) 111219.

[34] K. Luo, X. Kong, J. Li, J. Hu, L. Deng, Motion magnification for video-based vibration measurement of civil structures: A review, *Mechanical Systems and Signal Processing* 220 (2024) 111681.

[35] N. A. Valente, Z. Mao, C. Niezrecki, Holistically Nested Edge Detection and particle filtering for subtle vibration extraction, *Mechanical Systems and Signal Processing* 204 (2023) 110753.

[36] A. Oppenheim, J. Lim, The importance of phase in signals, in: *Proceedings of the IEEE*, Vol. 69, 1981, pp. 529-541.

[37] Y. Wang, W. Hu, J. Teng, Y. Xia, Phase-based motion estimation in complex environments using the illumination-invariant log-Gabor filter, *Mechanical Systems and Signal Processing* 186 (2023).

780

[38] X. Bai, Q. Zhu, X. Wang, Q. Zhang, Y. Du, Modal-Weighted Super-Sensitive phase optical flow method for structural Micro-Vibration modal identification, *Mechanical Systems and Signal Processing* 224 (2025) 112095.

[39] S. Liu, L. Yu, W. Niu, J. Wang, Z. Zhong, J. Huang, M. Shan, Fast and accurate visual vibration measurement via derivative-enhanced phase-based optical flow, *Mechanical Systems and Signal Processing* 209 (2024) 111089.

[40] Y. Wang, W. Hu, J. Teng, Y. Xia, Motion estimation from noisy data with unknown distributions using multi-frame phase-preserving denoising, *Mechanical Systems and Signal Processing* 206 (2024) 110924.

[41] B. K. P. Horn, B. G. Schunck, Determining optical flow, *Artificial Intelligence* 17 (1) (1981) 185-203.

[42] Z. Luo, X. Hu, L. Wang, Y. Lu, Measurement of large, discontinuous displacement from digital images, *Strain* 56 (4) (2020) e12349.

790

[43] J. Javh, J. Slavič, M. Boltežar, The subpixel resolution of optical-flow-based modal analysis, *Mechanical Systems and Signal Processing* 88 (2017) 89-99.

[44] S. Li, F. A. Azad, Y. Yang, Super-sensitivity full-field measurement of structural vibration with an adaptive incoherent
optical method, Mechanical Systems and Signal Processing 202 (2023) 110666.

[45] K. Čufar, J. Slavič, M. Boltežar, Mode-shape magnification in high-speed camera measurements, Mechanical Systems and Signal Processing 213 (2024) 111336.

[46] C. Xin, C. Wang, Z. Xu, M. Qin, M. He, Marker-free vision-based method for vibration measurements of RC structure under seismic vibration, Earthquake Engineering & Structural Dynamics 51 (8) (2022) 1773–1793.

[47] C. Xiu, Y. Weng, W. Shi, Vision and vibration data fusion-based structural dynamic displacement measurement with test validation, Sensors (Basel, Switzerland) 23 (9) (2023).

[48] X. Pan, T. Y. Yang, Y. Xiao, H. Yao, H. Adeli, Vision-based real-time structural vibration measurement through deep-learning-based detection and tracking methods, Engineering Structures 281 (2023).

[49] S. Das, P. Saha, S. K. Patro, Vibration-based damage detection techniques used for health monitoring of structures: a review, Journal of Civil Structural Health Monitoring 6 (3) (2016) 477–507.

[50] S. Pereira, E. Reynders, F. Magalhães, A. Cunha, J. P. Gomes, The role of modal parameters uncertainty estimation in automated modal identification, modal tracking and data normalization, Engineering Structures 224 (2020) 111208.

[51] O. Avci, O. Abdeljaber, S. Kiranyaz, M. Hussein, M. Gabbouj, D. J. Inman, A review of vibration-based damage detection in civil structures: from traditional methods to machine learning and deep learning applications, Mechanical Systems and Signal Processing 147 (2021) 107077.

[52] J. Shi, C. Tomasi, Good features to track, in: Proceedings of IEEE Conference on Computer Vision and Pattern Recognition, 1994, pp. 593–600.

[53] D. G. Lowe, Distinctive image features from scale-invariant keypoints, International Journal of Computer Vision 60 (2) (2004) 91–110.

[54] R. Yazdi, H. Khotanlou, E. Alighardash, M. Zolfaghari, Edge detection method based on the differences in intensities of rotating kernel borders, in: 6th International Conference on Pattern Recognition and Image Analysis (IPRIA), 2023.

[55] B. Merainani, B. Xiong, V. Baltazart, J. Dumoulin, M. Döhler, Q. Zhang, Experimental investigation of structural modal identification using pixels intensity and motion signals from video-based imaging devices: performance, comparison and analysis, in: Multimodal Sensing and Artificial Intelligence: Technologies and Applications II, SPIE, 2021.

[56] J. P. Lewis, Fast normalized cross-correlation, Vision interface 10 (1995) 120–123.

[57] R. M. Haralick, L. G. Shapiro, Computer and robot vision, Vol. 1, Addison-wesley Reading, MA, 1992.

[58] B. Pan, K. Li, W. Tong, Fast, robust and accurate digital image correlation calculation without redundant computations, Experimental Mechanics 53 (7) (2013) 1277–1289.

[59] M. A. Sutton, J.-J. Orteu, H. W. Schreier, Digital image correlation (DIC), in: Image Correlation for Shape, Motion and Deformation Measurements: Basic Concepts, Theory and Applications, Springer US, Boston, MA, 2009, pp. 81–117.

[60] W. Chojnacki, Z. L. Szpak, Determining ellipses from low-resolution images with a comprehensive image formation model, Journal of the Optical Society of America A 36 (2) (2019) 212–233.

[61] X. Bian, Video-based vibration analysis for structural health monitoring in civil engineering, PhD thesis, Université de Rennes 1 (Dec. 2021).

[62] KEYENCE (UK) Ltd, Sensor heads - IL-030, https://www.keyence.co.uk/products/sensor/positioning/il/models/il-030/.

[63] Wenglor, YP11MGVL80 Laser Distance Sensor Triangulation, https://www.wenglor.com/en/Sensors/Photoelectronic-Sensors/Laser-Sensors-for-distance-measurement/Laser-Distance-Sensors-Triangulation/Laser-Distance-Sensor-Triangulation/p/YP11MGVL80.

[64] D. Yatsenko, measureShift, https://www.mathworks.com/matlabcentral/fileexchange/57578-measureshift (2016).

[65] H. An, Y. Cao, H. Li, H. Zhang, Temporal phase unwrapping based on unequal phase-shifting code, IEEE Transactions on Image Processing 32 (2023) 1432–1441.