A Spatial Concordance Correlation Coefficient with an Application to Image Analysis

TL;DR

This paper introduces a spatial concordance coefficient for comparing spatial sequences along a fixed lag, explores its properties for specific covariance functions, and applies it to analyze digital images of forest canopies.

Contribution

It defines a new spatial concordance coefficient, establishes its asymptotic normality, and develops a local estimation approach for large digital images.

Findings

The coefficient was explored for Matérn and Wendland covariance functions.

Asymptotic normality was proven for Wendland covariance in increasing domain sampling.

The local approach revealed spatial decay of concordance in forest canopy images.

Abstract

In this work we define a spatial concordance coefficient for second-order stationary processes. This problem has been widely addressed in a non-spatial context, but here we consider a coefficient that for a fixed spatial lag allows one to compare two spatial sequences along a 45-degree line. The proposed coefficient was explored for the bivariate Mat\'ern and Wendland covariance functions. The asymptotic normality of a sample version of the spatial concordance coefficient for an increasing domain sampling framework was established for the Wendland covariance function. To work with large digital images, we developed a local approach for estimating the concordance that uses local spatial models on non-overlapping windows. Monte Carlo simulations were used to gain additional insights into the asymptotic properties for finite sample sizes. As an illustrative example, we applied this methodology to two similar images of a deciduous forest canopy. The images were recorded with different cameras but similar fields-of-view and within minutes of each other. Our analysis showed that the local approach helped to explain a percentage of the non-spatial concordance and to provided additional information about its decay as a function of the spatial lag.

Figures

Equations

$$\text{MSD}={ϵ}^2=\mathbb{E}[D^2]=\mathbb{E}[(Y-X{)}^2].$$

$${\rho }_c=1-\frac{{ϵ}^2}{{ϵ}^2|\rho =0}=\frac{2{\sigma }_{YX}}{{\sigma }_Y^2+{\sigma }_X^2+({\mu }_Y^2-{\mu }_X^2{)}^2}.$$

$${\hat{\rho }}_c=\frac{2s_{YX}}{s_Y^2+s_X^2+(\overline{y}-\overline{x}{)}^2}.$$

$$Z=\frac{1}{2}\left(\frac{1+{\hat{\rho }}_c}{1-{\hat{\rho }}_c}\right)\stackrel{\mathcal{D}}{⟶}\mathcal{N}(\psi ,{\sigma }_Z^2),\text{as}n\to \infty ,$$

$$\psi ={\tanh }^{-1}({\rho }_c)=\frac{1}{2}\left(\frac{1+{\rho }_c}{1-{\rho }_c}\right),$$

$${\sigma }_Z^2=\frac{1}{n-2}\left[\frac{(1-{\rho }^2){\rho }_c^2}{(1-{\rho }_c^2){\rho }^2}+\frac{2v^2(1-{\rho }_c){\rho }_c^3}{(1-{\rho }_c^2{)}^2\rho }+\frac{v^4{\rho }_c^4}{2(1-{\rho }_c^2{)}^2{\rho }^2}\right],$$

$$v^2=\frac{({\mu }_Y-{\mu }_X{)}^2}{{\sigma }_Y{\sigma }_X}.$$

$${\text{H}}_0:{\rho }_c={\rho }_0\text{versus}{\text{H}}_1:{\rho }_c\ne {\rho }_0$$

$${\hat{\rho }}_c\pm z_{\alpha /2}\sqrt{{\sigma }_Z^2}$$

$$C(\mathbf{h})=\left(\begin{array}{cc}{C}_X(\mathbf{h})& C_{XY}(\mathbf{h})\\ C_{YX}(\mathbf{h})& C_Y(\mathbf{h})\end{array}\right).$$

$${\rho }^c(\mathbf{h})$$

$$=\frac{2C_{YX}(\mathbf{h})}{C_X(0)+C_Y(0)+({\mu }_1-{\mu }_2{)}^2}.$$

$$C_X(\mathbf{h})$$

$$C_Y(\mathbf{h})$$

$${\mu }_1$$

$$C_{YX}(\mathbf{h},{\nu }_{12},a_{12})$$

$${\rho }^c(\mathbf{h})$$

$$=\frac{2{\rho }_{12}{\sigma }_1{\sigma }_{2}M(\mathbf{h}|{\nu }_{12},a_{12})}{{\sigma }_1^{2}M(0|{\nu }_1,a_1)+{\sigma }_2^{2}M(0|{\nu }_2,a_2)}$$

$$=\frac{2{\rho }_{12}{\sigma }_1{\sigma }_{2}M(\mathbf{h}|{\nu }_{12},a_{12})}{{\sigma }_1^2+{\sigma }_2^2}$$

$$=\frac{2{\sigma }_1{\sigma }_{2}M(\mathbf{h},{\nu }_{12},a_{12})}{{\sigma }_1^2+{\sigma }_2^2}$$

$$=\eta \cdot {\rho }_{12},$$

$$M(\mathbf{h}|{\nu }_{12},a_{12})=M(\mathbf{h}|n+1/2,a_{12})=\exp (-a_{12}\Vert \mathbf{h}\Vert )\sum _{k=0}^n\frac{(n+k)!}{(2n)!}\left(\begin{array}{c}n\\ k\end{array}\right)(2a_{12}\Vert \mathbf{h}\Vert {)}^{n-k},$$

$${\rho }^c(\mathbf{h})=\frac{2{\sigma }_{12}}{{\sigma }_1^2+{\sigma }_2^2}\exp (-a_{12}\Vert \mathbf{h}\Vert )\sum _{k=0}^n\frac{(n+k)!}{(2n)!}\left(\begin{array}{c}n\\ k\end{array}\right)(2a_{12}\Vert \mathbf{h}\Vert {)}^{n-k}.$$

$${\rho }^c(\mathbf{h})=\frac{2{\sigma }_{12}}{{\sigma }_1^2+{\sigma }_2^2}\exp (-a_{12}\Vert \mathbf{h}\Vert ).$$

$$\mathbf{C}(\mathbf{h})={\left[{\rho }_{ij}{\sigma }_{ii}{\sigma }_{jj}{R}_{ij}(\mathbf{h})\right]}_{i,j=1}^2,$$

$${\tilde{\psi }}_{v,k}(t)={\int }_t^1\frac{u(u^2-t^2{)}^{k-1}(1-u{)}_{+}^v}{B(2k,v+1)}du,0\le t\le 1$$

$${\rho }^c(\mathbf{h})=\frac{2{\rho }_{12}{\sigma }_1{\sigma }_{2}R(\mathbf{h},{\psi }_{12})}{{\sigma }_1^2+{\sigma }_2^2+({\mu }_1-{\mu }_2{)}^2},\mathbf{h}\in {\mathbb{R}}^2.$$

$${\rho }^c(\mathbf{h})=\frac{2{\rho }_{12}{\sigma }_1{\sigma }_2\left(1+l\Vert \mathbf{h}\Vert /b_{12}\right){\left(1-\Vert \mathbf{h}\Vert b_{12}\right)}_{+}^l}{{\sigma }_1^2+{\sigma }_2^2+({\mu }_1-{\mu }_2)}.$$

$${\gamma }_{XY}(\mathbf{h})=\mathbb{E}[(X(\mathbf{s}+\mathbf{h}))-X(\mathbf{s})(Y(\mathbf{s}+\mathbf{h})-(Y(\mathbf{s}))],$$

$${\rho }^c(\mathbf{h})=1-\frac{2{\gamma }_{xy}(\mathbf{h})}{2{\gamma }_{xy}(0)+({\mu }_1-{\mu }_2{)}^2+2C_{xy}(\mathbf{h})}.$$

$${\hat{\rho }}^c(\mathbf{h})={\hat{\rho }}_{12}(\mathbf{h}){\hat{C}}_{ab},$$

$$\mathbf{C}(\mathbf{h})={\left[{\rho }_{ij}{\sigma }_{ii}{\sigma }_{jj}\left(1+(\nu +1)\frac{\Vert \mathbf{h}\Vert }{b_{12}}\right){\left(1-\frac{\Vert \mathbf{h}\Vert }{b_{12}}\right)}_{+}^{\nu +1}\right]}_{i,j=1}^2,$$

$${\left(\nabla g(\mathbf{\theta }{)}^{\top }{\mathbf{F}}_n(\mathbf{\theta }{)}^{-1}\nabla g(\mathbf{\theta })\right)}^{-1/2}(g({\hat{\mathbf{\theta }}}_n)-g(\mathbf{\theta }))\stackrel{D}{\to }\mathcal{N}(0,1),\text{as}n\to \infty ,$$

$$g(\mathbf{\theta })=\frac{2{\rho }_{12}{\sigma }_1{\sigma }_2\left(1+(\nu +1)\frac{\Vert \mathbf{h}\Vert }{b_{12}}\right){\left(1-\frac{\Vert \mathbf{h}\Vert }{b_{12}}\right)}_{+}^{\nu +1}}{{\sigma }_1^2+{\sigma }_2^2},$$

Full text

††journal: Spatial Statistics

1 Introduction

In recent decades, concordance correlation coefficients have been developed in a variety of different contexts. For instance, in assay or instrument validation processes, the reproducibility of the measurements among trials or laboratories is of interest. When a new instrument is developed, it may be relevant to evaluate whether its performance is concordant with other, existing ones, or its results accord with a “gold standard.” There are also situations in which one is interested in comparing two methods without a designated gold standard or target values [Lin et al., 2002]. In the literature, this latter type of concordance has been tackled from different perspectives [Barnhart et al., 2007]. Cohen [1968] discussed this problem in the context of categorical data. Schall and Williams [1996] and Lin [2000] performed similar studies in the context of bioequivalence.

One way to approach the concordance problem for continuous measurements is to construct a scaled summary index that can take on values between $-1$ and 1, analogous to a correlation coefficient. Using this approach, Lin [1989] suggested a concordance correlation coefficient (CCC) that evaluates the agreement between two continuous variables by measuring their joint deviation from a 45∘ line through the origin. There have been some extensions of this CCC that use several measuring instruments and techniques to evaluate the agreement between two instruments; these efforts have led to interesting graphical tools [Hiriote and Chinchilli, 2011, Stevens et al., 2017]. In the context of goodness of fit, Vonesh et al. [1996] proposed a modified Lin’s CCC for choosing models that have a better agreement between observed and the predicted values. Recently Stevens et al. [2017] and Chodhary and Nagaraja [2017] developed the probability of agreement, and Leal et al. [2019] studied the local influence of the CCC and the probability of agreement considering both first- and second-order measures under the case-weight perturbation scheme. Atkinson and Nevill [1997] critiqued the CCC because any correlation coefficient is highly dependent on the measurement range. In general, therefore, CCC is used only when measuring ranges are comparable or when methods are on the same scale.

In this paper, we suggest an approach to assessing the agreement between two continuous responses when the observations of both variables have been georeferenced in space. We define a spatial CCC (SCCC) as a generalization of Lin’s (1989) coefficient that measures the agreement between two spatial variables. For a fixed lag, our SCCC shares the same properties as the original CCC. For an increasing sampling scheme, we establish the asymptotic normality of the sample SCCC for a bivariate Gaussian process with a Wendland covariance function. To improve the behavior of the coefficient, we developed a local approach for estimating it that uses local spatial models on non-overlapping windows. This approach constitutes a new way of thinking about concordance that has not been considered previously, especially for large digital images. Our approach also captures the decay of the SCCC as a function of the norm of the spatial lag. Monte Carlo simulations and numerical experiments with real datasets accompany the exposition of the methodological aspects. An image-analysis example is worked in detail to illustrate the fitting of a local SCCCs. We conclude with a summary of the main findings and an outline of problems to be tackled in future research.

2 Preliminaries and Notation

Assume that $X$ and $Y$ are two continuous random variables such that the joint distribution of $X$ and $Y$ has finite second moments with means $\mu_{X}$ and $\mu_{Y}$, variances $\sigma^{2}_{X}$ and $\sigma_{Y}^{2}$, and covariance $\sigma_{YX}$. The mean squared deviation of $D=Y-X$ is

[TABLE]

It is straightforward to see that $\epsilon^{2}=(\mu_{X}-\mu_{Y})^{2}+\sigma_{Y}^{2}+\sigma_{X}^{2}-2\sigma_{YX}$ and the sample counterpart satisfies $e^{2}=(\overline{y}-\overline{x})^{2}+s_{Y}^{2}+s_{X}^{2}-2s_{XY}.$ Using this framework, Lin [1989] defined a CCC as:

[TABLE]

The CCC satisfies the following properties:

$\rho_{c}=\alpha\cdot\rho,$ where $\alpha=\frac{2}{w+1/w+v^{2}}$ and $w=\frac{\sigma_{Y}}{\sigma_{X}}.$

• 2.

$|\rho_{c}|\leq 1.$

• 3.

$\rho_{c}=0$ if and only if $\rho=0.$

• 4.

$\rho_{c}=\rho$ if and only if $\sigma_{Y}=\sigma_{X}$ and $\mu_{Y}=\mu_{X}$.

The sample estimate of $\rho_{c}$ is given as

[TABLE]

The inference for this coefficient was addressed via Fisher’s transformation. Lin [1989] proved that

[TABLE]

where

[TABLE]

[TABLE]

and

[TABLE]

As a consequence of the asymptotic normality of the sample CCC, an approximate hypothesis testing problem of the form

[TABLE]

for a fixed $\rho_{0}$ can be constructed. Alternatively, an approximate confidence interval of the form

[TABLE]

can be used, where $z_{\alpha/2}$ is the upper quantile of order $\alpha/2$ of the standard normal distribution. Applications and extensions of Lin’s coefficient can be found in Lin et al. [2012], among others.

3 A Spatial Concordance Coefficient and its Properties

We start by extending Lin’s CCC for bivariate second-order spatial processes for a fixed lag in space.

Definition 1**.**

Let $(X(\bm{s}),Y(\bm{s}))^{\top}$ be a bivariate second-order stationary random field with $\bm{s}\in\mathbb{R}^{2}$, mean $(\mu_{1},\mu_{2})^{\top}$, and covariance function

[TABLE]

Then the SCCC is defined as

[TABLE]

Some straightforward properties of this SCCC are:

$\rho^{c}(\bm{h})=\eta\cdot\rho_{YX}(\bm{h}),$ where $\eta=\frac{2\sqrt{C_{X}(\bm{0})C_{Y}(\bm{0})}}{C_{X}(\bm{0})+C_{Y}(\bm{0})+(\mu_{1}-\mu_{2})^{2}}.$

• 2.

$|\rho^{c}(\bm{h})|\leq 1.$

• 3.

$\rho^{c}(\bm{h})=0$ iff $\rho_{YX}(\bm{h})=0.$

• 4.

$\rho^{c}(\bm{h})=\rho_{YX}(\bm{h})$ iff $\mu_{1}=\mu_{2}$ and $C_{X}(\bm{0})=C_{Y}(\bm{0}).$

• 5.

For a bivarite Matérn covariance function defined as [Gneiting et al., 2010]

[TABLE]

where $M(\bm{h},\nu,a)=(a||\bm{h}||)^{\nu}K_{\nu}(a||\bm{h}||)$, $K_{\nu}(\cdot)$ is a modified Bessel function of the second kind, and $\rho_{12}=\text{cor}[X(\bm{s}_{i}),Y(\bm{s}_{j})]$, it follows that

[TABLE]

where $\eta=\frac{2\sigma_{1}\sigma_{2}M(\bm{h},\nu_{12},a_{12})}{\sigma_{1}^{2}+\sigma_{2}^{2}}.$

A special case of the Matérn covariance function is when $\nu_{12}=n+1/2$. Then

[TABLE]

and the SCCC is

[TABLE]

By choosing $n=0$ and $\nu_{12}=1/2$, $M(\bm{h}|1/2,a_{12})=\exp(-a_{12}\|\bm{h}\|)$. This gives the SCCC in its simplest form:

[TABLE]

For illustrative purposes, consider $\sigma_{1}=1$, $\sigma_{2}=2$, $\sigma_{12}=1.8$, $a_{12}=1/2$ and $\nu_{12}=\{\frac{1}{2},\frac{3}{2},\frac{5}{2}\}$. Then for $\bm{h}\in\{0,1,\ldots,15\}$, the patterns of the SCCC for the Matérn covariance function are illustrated in Figure 1. The range of the SCCC increases as $\nu$ increases.

For a bivariate Wendland-Gneiting covariance function [Daley et al., 2015] of the form

[TABLE]

where $R(\bm{h},\psi_{12})=c_{ij}b_{ij}^{\nu+2k+1}B(\nu+2k+1,\gamma_{ij}+1)\tilde{\psi}_{\nu+\gamma_{ij}+1,k}\left(\dfrac{\|\bm{h}\|}{b_{ij}}\right),~{}B(\cdot,\cdot)$ is the beta function, and $\tilde{\psi}_{\nu,k}$ is defined, for $k\geq 1$ [Gneiting, 2002], as

[TABLE]

the SCCC is

[TABLE]

In particular, considering $R_{ij}(\bm{h})=p_{k}(\|\bm{h}\|)(1-\|\bm{h}\|/b_{ij})_{+}^{l}$, where $k=1$, $l=v+\gamma+1,$ $\gamma=0$ and $b_{ij}>0$,

[TABLE]

Using similar arguments as in properties 1–6, the SCCC could be derived for other parametric bivariate correlation functions.

For a bivariate intrinsically stationary random field $(X(\bm{s}),Y(\bm{s}))^{\top},~{}\bm{s}\in\mathbb{R}^{2}$, with cross-variogram given by

[TABLE]

there is another characterizations of the SCCC defined in equation (1). The coefficient can be written as

[TABLE]

Because $\rho^{c}(\cdot)$ in equation (8) depends on $C_{XY}(\cdot)$ and $\gamma_{XY}(\cdot)$, we prefer the representation in equation (1).

4 Inference

In the previous section we proved that for several covariance structures, the spatial concordance correlation coefficient defined in equation (1) can be written as a product of the correlation coefficient and a constant. Thus, we can consider plug-in estimators for the correlation coefficient and the constant.

Let $(X(\bm{s}),Y(\bm{s}))^{\top},~{}s\in D\subset\mathbb{R}^{2}$ be a Gaussian process with mean $\bm{\mu}=(\mu_{1},\mu_{2})^{\top}$ and covariance function $\bm{C}(\bm{h})$, $\bm{s},\bm{h}\in\mathbb{R}^{2}$. Then a sample estimate of the SCCC index (1) is

[TABLE]

where $\widehat{C}_{ab}=((\widehat{a}+1/\widehat{a}+\widehat{b}^{2})/2)^{-1},\quad\widehat{a}=\left(\dfrac{\widehat{C}_{11}(\bm{0})}{\widehat{C}_{22}(\bm{0})}\right)^{1/2},\quad\widehat{b}=\dfrac{\widehat{\mu_{1}}-\widehat{\mu}_{2}}{(\widehat{C}_{11}(\bm{0})\widehat{C}_{22}(\bm{0}))^{1/4}},$ and $\widehat{\mu}_{1}$, $\widehat{\mu}_{2},$ $\widehat{C}_{11}(\bm{0})$ and $\widehat{C}_{22}(\bm{0})$ are the maximum likelihood (ML) estimates of $\mu_{1}$, $\mu_{2},$ $C_{11}(\bm{0})$, and $C_{22}(\bm{0})$, respectively.

The asymptotic properties of an estimator as in equation (9) have been studied in the literature for specific cases. Bevilaqua et al. [2015] studied the asymptotic properties of the ML estimator for a separable Matérn covariance model. They used a result provided by Mardia and Marshall [1984] in an increasing domain sampling framework. Using this theorem and the delta method, we can establish the following result for the Wendland-Gneiting model:

Theorem 1**.**

Let $\bm{Z}(\bm{s})=(X(\bm{s}),Y(\bm{s}))^{\top}$, $s\in D\subset\mathbb{R}^{2}$ be a bivariate Gaussian spatial process with mean $\bm{0}$ and covariance function given by

[TABLE]

for $\nu>0$ fixed. Define $\bm{\theta}=(\sigma_{1}^{2},\sigma_{2}^{2},\rho_{12},b_{12})^{\top}$ and denote $\widehat{\bm{\theta}}_{n}$ the ML estimator of $\bm{\theta}.$ Then

[TABLE]

in an increasing domain sense, where

[TABLE]

$\bm{F}_{n}(\bm{\theta})^{-1}$ is the covariance matrix of $\widehat{\bm{\theta}}_{n},$

[TABLE]

and $f(b_{12})=\left(-\dfrac{(\nu+1)\|\bm{h}\|}{b_{12}^{2}}\right)\left(1-\dfrac{\|\bm{h}\|}{b_{12}}\right)^{\nu+1}_{+}+\left(1+\dfrac{(\nu+1)\|\bm{h}\|}{b_{12}}\right)\left(1-\dfrac{\|\bm{h}\|}{b_{12}}\right)^{\nu}_{+}\dfrac{(\nu+1)\|\bm{h}\|}{b_{12}^{2}}.$

Proof.

See the Appendix. ∎

5 A Local Approach

When the size of the images is large, it is difficult to find a single model fitting reasonably well to an entire image. This has been investigated in the literature for autoregressive processes defined on the plane in the context of image restoration and segmentation. For examples, see Bustos et al. [2009] and Ojeda et al. [2010].

Here we describe a local approach for a bivariate process of the form $\bm{Z}(\bm{s})=(\bm{Z}_{1}(\bm{s}),\bm{Z}_{2}(\bm{s}))^{\top}$, $\bm{s}\in D\subset\mathbb{R}^{2}$, where the observations are located over a rectangular grid of size $n\times m$. The extension to an $l\in\mathbb{N}$-variate process is natural when $l>2$. In this framework, we assume that the whole domain $D$ can be divided into $p$ sub-windows $D_{i}$, such that $\cup_{i=1}^{n}D_{i}=D$, for $i=1,...p$. Then we define $p$ processes of the form $\bm{Z}_{i}(\bm{s})=(\bm{Z}_{i1}(\bm{s}),\bm{Z}_{i2}(\bm{s}))^{\top},~{}\bm{s}\in D_{i}$, where each process has a covariance function given by

[TABLE]

Then for each local process $\bm{Z}_{i}(\cdot)$ we define the local SCCC $\rho_{i}^{c}(\cdot)$ using the theory developed in Section 3:

[TABLE]

Based on the local coefficients $\rho_{i}^{c}(\cdot)$, we suggest two global SCCCs. The first one is the average of the $p$ local coefficients, given by

[TABLE]

The second one considers the average of each parameter in the correlation function such that the global coefficient is

[TABLE]

where $\overline{\sigma}_{1}=\frac{1}{p}\sum_{i=1}^{p}\sigma_{1i}^{2}$, and similarly for $\overline{\sigma}_{2},~{}\overline{\rho}_{12},$ and $\overline{\bm{\psi}}$. As a result we have two global coefficients of spatial concordance depending on averages, the first one is the average of the local coefficients and the second one is a plug in of the parameter averages.

When process $\bm{Z}(\bm{s})$ have been observed in the sites $\bm{s}_{1},\ldots,\bm{s}_{n}$ and all the local coefficients have been computed, the sample versions of $\rho_{1}(\cdot)$ and $\rho_{2}(\cdot)$ are

[TABLE]

where $\widehat{\overline{\sigma}}_{1}$, $\widehat{\overline{\sigma}}_{2}$, $\widehat{\overline{\rho}}_{12}$, and $\widehat{\overline{\bm{\psi}}}$ are the ML estimations of the parameters defined in equation (12).

Considering and increasing domain sampling scheme, the asymptotic normality of $\widehat{\rho}_{1}(\bm{h})$ is straightforward. Indeed, let $\bm{Z}_{i}(\bm{s})=(Z_{i1}(\bm{s}),Z_{i2}(\bm{s}))^{\top},~{}s\in D_{i}$, be a bivariate process with correlation structure given by $\bm{C}_{i}(\bm{h})=\left[\rho^{i}_{jk}\sigma^{i}_{jj}\sigma^{i}_{kk}R_{i}(\bm{h},\psi_{i})\right]_{j,k=1}^{2},~{}i=1,...p.$ Define the parameter vector $\bm{\theta}^{i}=(\rho_{12}^{i},\sigma_{1}^{i},\sigma_{2}^{i},\sigma_{12}^{i},\psi_{i})^{\top}$ associated with $\bm{Z}_{i}(\bm{s})$. If the covariance satisfies the Mardia and Marshall [1984] conditions, then

[TABLE]

where $\bm{F}^{i}_{n}(\bm{\theta}^{i})$ is the covariance matrix of $\widehat{\bm{\theta}}^{i}$. Then for $g(\bm{\theta})=\rho_{i}^{c}(\bm{h})$, we have that

[TABLE]

Now assuming that $\bm{Z}_{i}(\bm{s})$ and $\bm{Z}_{j}(\bm{s})$ are independent for all $i\neq j$, we get

[TABLE]

6 Monte Carlo Simulations

We used Monte Carlo simulation to explore the properties of the SCCC, $\rho^{c}(\cdot)$, for finite samples sizes. The performance of the ML estimations were then analyzed with respect to the true values of the coefficient. We generated 500 replicates from a Gaussian random field sampled on a regular lattice of size $20\times 20$ inside the region $[-\frac{3}{2},\frac{3}{2}]^{2}.$ Each replicate was generated from a bivariate Gaussian random field with mean zero and Wendland-Gneiting covariance function given in equation (6). In each case, we estimated the parameters of the covariance function using ML and used them to compute the SCCC given in equation (6). Three set of parameters were considered, one set each for $\sigma_{1}=\sigma_{2}=1$, $\nu=4$ and $k=1$:

Case 1: $\rho_{12}=-0.15$, $b_{1}=0.5$, $b_{2}=0.4$, and $b_{12}=0.35$. 2. 2.

Case 2: $\rho_{12}=0.25$, $b_{1}=1.2$, $b_{2}=0.9$, and $b_{12}=1$. 3. 3.

Case 3: $\rho_{12}=0.3$, $b_{1}=1.8$, $b_{2}=1.4$, and $b_{12}=1.5$.

In Figure 2 we show a realization of the random field for each case.

The ML estimates of the parameters of the Wendland-Gneiting covariance function had low bias and standard errors, and agreed with previously published results (e.g., Bevilaqua et al., 2019). Using these estimates, we computed the SCCC in each case for $0<||\bm{h}||<2$. The mean square errors of the estimates were bounded by $3.9404\cdot 10^{-5}$, $9.7958\cdot 10^{-5}$, and $0.0002$, respectively, for cases 1–3. $\rho^{c}(\bm{h})$ versus $||h||$ and $\widehat{\rho}^{c}(\bm{h})$ versus $||h||$ are plotted in Figure (3); the true coefficient is drawn with a continuous line. The estimates of the SCCC were reasonably well-behaved but worsened when $||h||$ was close to zero, as is typical of lag-dependent spatial functions computed over a rectangular grid where the minimum distance between coordinates is fixed. The general Monte Carlo simulation study also involved the bivariate Matérn covariance function and the results were similar. The estimate of $\rho_{12}$ was better for the Matérn case in terms of the mean square error. This is important because in both cases, the estimate of $\rho_{12}$ affected the estimate of the SCCC.

Finally, for the same region used in the previous Monte Carlo simulation, we computed the asymptotic variance of $\widehat{\rho}^{c}(\cdot)$. For $0<||\bm{h}||<2$, all variances are bounded by 0.006, and the largest discrepancies between cases 1–3 were seen near the origin.

To gain more insight into the computational time required for computing $\widehat{\rho}^{c}(\cdot)$ for the covariance functions used in this work, we ran similar simulations with different window sizes. We ran 100 simulations, and in each, $\widehat{\rho}^{c}(\cdot)$ was computed for the Matérn and Wendland-Gneiting covariance functions for window sizes = $8\times 8,$ $12\times 12$, $16\times 16$, and $20\times 20$. All computations were done using an HP ProLiant DL380G9 server, equipped with a 2x Intel Xeon E5-2630 v3 2.40 GHz processor, 128 GB DDR4 2.133 Ghz RAM, and 512 GB SSD storage.

Time to run each simulation increased exponentially with window size (Figure 4). Although the time required to compute the Wendland-Gneiting covariance function was always smaller than the time to compute the Matérn covariance function, for real images it is not feasible to compute $\widehat{\rho}^{c}(\cdot)$, at least using an interpreted language like R as we did here. This result further supports the use of the local approach we presented in Section 5, but we will continue to explore ways to optimize and accelerate the computation of $\widehat{\rho}^{c}(\cdot)$.

7 An Application

7.1 Motivation

Our application derives from ecology. In order to track the seasonality (“phenology”) of vegetation activity in different ecosystems, digital cameras have been deployed to record high-frequency images of the canopy at hundreds of research sites around the world [Richardson, 2018]. From each image, color-channel information (e.g., RGB [red-green-blue] values of each pixel) are extracted and converted to a suite of “vegetation indices” derived from linear or nonlinear transformations of the RGB or other color spaces [Sonnentag et al., 2012, Mizunuma et al., 2014, Toomey et al., 2015, Nguy-Robertson et al., 2016]. These indices have been used to identify the timing of seasonal phenomena such as leaf-out, senescence, and abscission, and to monitor how these phenomena are changing in response to ongoing climatic change [Sonnentag et al., 2012]. However, different cameras may render the same scene differently because of the specifics of the imaging sensor being used (e.g., CCD, CMOS) and researchers have used a wide range of different cameras because of considerations including trade-offs between cost and image quality. Additionally, changes in scene illumination (e.g., caused by time-of-day or cloud cover) also may impact the resulting image. Although previous research has shown that diurnal, seasonal, and weather-related changes in illumination can have large effects on estimates of average color (or color index) for the whole image or a region of interest [Sonnentag et al., 2012], spatial information has not been incorporated previously in these estimates.

7.2 Imagery

We focus here on comparing two jpeg images taken of the same scene on 20 October 2010 by two different cameras (Figures 5, 5). These images were taken with, respectively, an outdoor StarDot NetCam XL 3MP camera with a $2048\times 1636$-pixel CMOS sensor (Figure 5) and an outdoor Axis 223M camera with a $1600\times 1200$-pixel CCD sensor (Figure 5). These images were selected from the image archive associated with an experiment, analyzed and reported on previously by Sonnentag et al. [2012], in which images, color time series, and phenological transition dates from eleven different cameras were compared. Although the two images we use here are of the same scene and were taken at the same time, they are not identical. For example, both cameras were pointing due north with an $\approx 20\degree$ tilt angle, but image displacement occurred because the cameras were mounted at different positions on a fixed platform. The resolution and overall field-of-view also differed because of different sensor sizes and lens characteristics. Sonnentag et al. [2012] compared color information averaged across a small “region of interest” in the images. Here, we work with the entire images after correction for differences of field-of-view and displacement.

To account for differences in field-of-view and displacement, the two images were first manually cropped using tools in IrfanView (version 4.38; Skiljan 2014) to equivalent areas and aspect ratios. The resulting images had $2023\times 1444$ pixels for the higher-resolution one taken with the StarDot camera and $1297\times 922$ pixels for the lower-resolution one taken with the Axis camera. The higher-resolution image was then resized and down-sampled in IrfanView so that it had the same number of pixels as the lower-resolution image (Figures 6, 6). These two images were loaded into the R software package (version 3.51; R Core Team, 2018) using the load.image function in the imager package [Urbanek, 2014] and transformed either to gray-scale using the grayscale function in the same package (Figures 7, 7) or to green chromatic coordinates (gcc), which normalizes for brightness ($g_{cc}=\frac{G}{R+G+B}$; Gillespie et al., 1987) (Figures 8, 8). For both the gray-scale and gcc images, the lower-resolution image (Figures 7 and 8, respectively) was then coordinate-registered to the higher-resolution image (Figure 7 and 8, respectively) using the R package RNiftyReg and a linear (affine) transformation with 12 degrees of freedom [Clayton et al., 2018]. Spatial concordance was assessed between the resampled higher-resolution images (Figure 7 or 8) and the coordinate-registered lower-resolution images (Figure 7 or 8).

7.3 Estimating concordance

For each pair of images, we first calculated Lin’s (1989) CCC. We then calculated the SCCC as described in Section 5. We calculated the local concordance coefficient $\rho_{i}(\cdot)$ in small ($12\times 12$-pixel) non-overlapping windows. To fit the local model to each small window, we used a Gaussian process $\bm{Z}(\bm{s})=(Z_{1}(\bm{s}),Z_{2}(\bm{s}))^{\top},\ \bm{s}\in\mathbb{R}^{2}$, with mean $(\mu_{1},\mu_{2})^{\top}$ and the covariance functions described in equations (3)–(6). We used the function GeoFit in the R package GeoModels [Bevilacqua and Morales-Oñate, 2018] to compute the ML estimators of the parameters involved in the models. For computational efficiency, the Matérn and Wendland-Gneiting covariances were estimated for a randomly-selected set of ten $20\times 20$-pixel subimages; the one to be used was selected based on the Akaike and Bayesian Information Critera (AIC and BIC, respectively). In general, the AIC and BIC coefficients were smaller for estimates of the Matérn covariance than for the Wendland-Gneiting covariance, and so we used the Matérn model even though it took somewhat more time to use it to compute the local estimators. Finally, the global SCCCs for each pair of images were estimated using equations (11) and (12).

7.4 Estimates of concordance

Lin’s coefficient was $\rho_{c}=0.1334$ for the grayscale images (Figures 7 and 7) and $\rho_{c}=0.2450$ for the gcc-indexed images (Figures 8 and 8). In Figure 9 we plotted Lin’s coefficient and the two global coefficients as a function of the spatial norm. We observed a rapid decay of $\widehat{\rho}_{2}(\cdot)$ and a slower decay of $\widehat{\rho}_{1}(\cdot)$. The decay was related to the way in which the estimates were computed for each window: $\widehat{\rho}_{1}(\cdot)$ is a coefficient obtained by plugging in the average of the parameters in the concordance function, but $\widehat{\rho}_{2}(\cdot)$ is the average of the concordance using all possible windows.

For $\bm{h}=\bm{0}$, we observed that the SCCC was approximately one-third ($0.08/0.245\times 100\approx 0.33$) of Lin’s CCC. This suggests that Lin’s CCC overestimates the spatial concordance between these two images, and implies that would be inappropriate to use it for modeling spatial data.

The images and all the code used in this paper are available from https://harvardforest.fas.harvard.edu/harvard-forest-data-archive, dataset HF322.

8 Discussion

With the work presented herein, we have extended the standard methodology for estimating concordance into the spatial domain. Our approach consisted in defining a new coefficient that preserves the interpretation of Lin’s (1989) concordance correlation coefficient (CCC) for two spatial variables and for a fixed distance lag. Our new spatial concordance correlation coefficient (SCCC) compares the correlation between two spatial variables with respect to their fit to a 45°line that passes through the origin. The properties of Lin’s (1989) CCC are inherited by our SCCC. The ML estimator of our SCCC for the Wendland-Gneiting covariance function is asymptotically normal for an increasing domain sampling scheme. We defined a local SCCC and established its asymptotic normality for the sample version. From the local SCCC, we derived two estimates for the overall SCCC, one based on the average of the p local coefficients and the other based on the average of the parameters in the correlation function. Deriving the global SCCC from local coefficients estimated in small non-overlapping windows is computationally more efficient and permits the estimation of spatial concordance for large images that are used commonly in a wide range of applications.

The Monte Carlo simulation study presented in Section 6 revealed that for the Matérn and Wendlang-Gneiting covariance functions, the sample version of the SCCC produced accurate estimates of the SCCC that decreased with distance (spatial lag). However, the time required to compute SCCC grows exponentially with window size, implying that for a large image size it is unfeasible to compute $\widehat{\rho}(\cdot)$ using an interpreted language like R. Although we are exploring ways to improve computational efficiency, the local approach introduced here (Section 5) appears to be a straightforward way to estimate SCCC for large images.

The camera comparison experiment conducted by Sonnentag et al. [2012] found that images recorded with a variety of different camera makes and models, all mounted on the top of the same canopy access tower and with a similar field of view, varied in visual appearance, including color balance, saturation, contrast, and brightness. These differences can be attributed to internal differences in sensor design and image processing, and external factors such as lighting. However, Sonnentag et al. [2012] also found that when simple normalized indices were calculated from the image data, and the emphasis was placed on the seasonality—rather than absolute magnitude—of those indices, the phenological information derived from the imagery was extremely similar across all cameras. Notably, their analysis focused on information about the average color across a large “region of interest” drawn across the canopy [Sonnentag et al., 2012]. Although this approach is widely used [Richardson, 2018] and it has the advantage of enabling integration across multiple individuals or species that may comprise a typical forest canopy, it lacks spatial information.

The spatial concordance correlation coefficient we developed and presented here summarizes and accounts for the spatial information in the images, permitting more rigorous characterization of agreement between high-resolution digital images recorded by different sensors. Other applications include using images from different satellite remote-sensing platforms as part of ongoing efforts to harmonize , for example, imagery with different spatial resolution, spectral sensitivity, and angular characteristics [e.g., Landsat-Sentinel efforts Claverie et al., 2018]. Calculation of concordance statistics before and after sensor harmonization could provide critical and objective information about the success of different harmonization methods. There also could be potential applications in the fusion of remotely-sensed data obtained at different spatiotemporal resolutions, such as MODIS, with its 500-m spatial resolution, daily temporal resolution, and Landsat, with its 30-m spatial resolution, 16-day temporal resolution [Gao et al., 2015].

9 Future Work

Several related theoretical and applied problems arise from the methodology suggested in this article that would be fruitful directions for future research. First, SCCC could be applied to images taken at to points in time by the same camera. The decay of the SCCC as a function of the norm would be expected to be similar to that seen in Figure 9 for each sequential pair of images. Another approach for dealing with the same problem would be to consider a sequence of n images taken with the same camera to be a spatiotemporal process. Then, the SCCC and its estimation properties could be studied in that context. This generalization of the SCCC would have applications in, for example, spatiotemporal analysis of satellite images taken weeks, months, or years apart as a way of characterizing patterns of landscape change.

Acknowledgments

This work has been partially supported by the AC3E, UTFSM, under grant FB-0008, and by grants to A.D.R. from the US National Science Foundation (EF-1065029, EF-1702697, DEB-1237491) and the United States Geological Survey (G10AP00129). This is a contribution from the Harvard Forest Long-term Ecological Research (LTER) site, supported most recently by the US National Science Foundation (DEB 18-32210).

Appendix

Mardia and Marshall Theorem

Let $\left\{Y(\bm{s}):\bm{s}\in D\subset{\small\textsf{R}}^{d}\right\}$ be a Gaussian random field such that $Y(\cdot)$ is observed on $D_{n}\subset D$. It is assumed that $D_{n}$ is a non-random set satisfying $\|\bm{s}-\bm{t}\|\geq\gamma>0$ for all $\bm{s},\bm{t}\in D_{n}$. This ensures that the sampling set is increasing as $n$ increases. Denote \bm{Y}=(Y(\bm{s}_{1}),\dots,Y(\mbox{\boldmath{s}}_{n}))^{\top} and assume that \mathbb{E}[\bm{Y}]=\bm{X}\mbox{\boldmath{\beta}}, $\operatorname{cov}(Y(\bm{t}),Y(\bm{s}))=\sigma(\bm{t},\bm{s};\bm{\theta})$, $\bm{X}$ is $n\times p$ with \text{rank}(\bm{X})=p,~{}\mbox{\boldmath{\beta}}\in\mathbb{R}^{p}, and $\bm{\theta}\in\Theta$, where $\Theta$ is an open set of $\in\mathbb{R}^{q}$. Let \mbox{\boldmath{\Sigma}}=\bm{\Sigma}(\mbox{\boldmath{\theta}}) be the covariance matrix of $\bm{Y}$ such that the $ij$-th element of $\bm{\Sigma}$ is $\sigma_{ij}=\sigma(\bm{s}_{i},\bm{s}_{j};\bm{\theta})$. We can estimate $\bm{\theta}$ and $\bm{\beta}$ using ML, by maximizing

[TABLE]

where k is a constant.

Let $\bm{L}_{n}^{(1)}=\nabla L=(\bm{L}_{\bm{\beta}}^{\top},\bm{L}_{\bm{\theta}}^{\top})^{\top}$ and

[TABLE]

be the gradient vector and Hessian matrix, respectively, obtained from equation (A.1). Let $\bm{F}_{n}=-\mathbb{E}[\bm{L}_{n}^{(2)}]$ be the Fisher information matrix with respect to $\bm{\beta}$ and $\bm{\theta}.$ Then, $\bm{F}_{n}=\operatorname{diag}(\bm{F}_{n}{(\bm{\beta})},\bm{F}_{n}{(\bm{\theta}})),$ where $\bm{F}_{n}{(\bm{\beta})}=-\mathbb{E}[\bm{L}_{\bm{\beta}\bm{\beta}}]$ and $\bm{F}_{n}{(\bm{\theta})}=-\mathbb{E}[\bm{L}_{\bm{\theta}\bm{\theta}}]$.

For a twice differentiable covariance function $\sigma(\cdot,\cdot;\bm{\theta})$ on $\Theta$ with continuous second derivatives, Mardia and Marshall [1984] provided sufficient conditions on $\bm{\Sigma}$ and $\bm{X}$ such that the limiting distribution of $(\widehat{\bm{\beta}}_{n}^{\top},\widehat{\bm{\theta}}_{n}^{\top})^{\top}$ is normal, per the following:

Theorem. Let $\lambda_{1}\leq\cdots\leq\lambda_{n}$ be the eigenvalues of $\bm{\Sigma}$, and let those of \bm{\Sigma}_{i}=\dfrac{\partial\mbox{\boldmath{\Sigma}}}{\partial\theta_{i}} and \bm{\Sigma}_{ij}=\dfrac{\partial^{2}\mbox{\boldmath{\Sigma}}}{\partial\theta_{i}\partial\theta_{j}} be $\lambda_{k}^{i}$ and $\lambda_{k}^{ij}$, $k=1,\dots,n,$ such that $|\lambda_{1}^{i}|\leq\cdots\leq|\lambda_{n}^{i}|$ and $|\lambda_{1}^{ij}|\leq\cdots\leq|\lambda_{n}^{ij}|$ for $i,j=1,\cdots,q$. Suppose that as $n\rightarrow\infty$

1. (i)

$\lim\lambda_{n}=C<\infty$, $\lim|\lambda_{n}^{i}|=C_{i}<\infty$ and $\lim|\lambda_{n}^{ij}|=C_{ij}<\infty$ for all $i,j=1,\dots,q$. 2. (ii)

$\|\bm{\Sigma}_{i}\|^{-2}=\mathcal{O}(n^{-\frac{1}{2}-\delta})$ for some $\delta>0$, for $i=1,\dots,q$. 3. (iii)

For all $i,j=1,\dots,q$, $a_{ij}=\lim\left[t_{ij}/(t_{ii}t_{jj})^{\frac{1}{2}}\right]$ exists, where t_{ij}=\operatorname{tr}\left(\mbox{\boldmath{\Sigma}}^{-1}\bm{\Sigma}_{i}\bm{\Sigma}^{-1}\mbox{\boldmath{\Sigma}}_{j}\right) and $\bm{A}=(a_{ij})$ is nonsingular. 4. (iv)

$\lim(\bm{X}^{\top}\bm{X})^{-1}=\bm{0}.$

Then, $(\widehat{\bm{\beta}}_{n}^{\top},\widehat{\bm{\theta}}_{n}^{\top})^{\top}\xrightarrow[]{\mathcal{L}}\mathcal{N}\left((\bm{\beta}^{\top},\bm{\theta}^{\top})^{\top},\bm{F}_{n}^{-1}\right)$ as $n\rightarrow\infty$, in an increasing domain sense.

Proof of Theorem 1

The proof consists of verifying the Mardia and Marshall [1984] conditions. In Theorem 1, $\mathbb{E}[{\bm{Z}(\bm{s})}]=\bm{0}$; thus the fourth condition in Mardia and Marshall’s 1984 theorem (above), $\lim(\bm{X}^{\top}\bm{X})^{-1}=\bm{0}$, is trivially satisfied. Satisfying the first three conditions is somewhat more complex.

For the first two conditions, we start by considering $\nu$ to be fixed. Then

[TABLE]

Let us consider an increasing domain scenario for process $\bm{Z}(\bm{s})$, with points $\bm{s}_{1},...,\bm{s}_{n}$ located in a rectangle $D_{n}\subset\Delta\mathbb{Z}^{d}$, for $0<\Delta<L$, and $D_{n}\subset D_{n+1}$, for all $n$.

Define the distance matrix $\bm{H}_{n}=\left[H_{lq}\right]_{l=q=1}^{n}$, where $H_{lq}=\|\bm{s}_{l}-\bm{s}_{q}\|$, and $\|\cdot\|$ denotes the Euclidean norm. Then the covariance matrix of $(\bm{Z}(\bm{s}_{1})^{\top},\ldots,\bm{Z}(\bm{s}_{n})^{\top})^{\top}$ can be written as

[TABLE]

where $\bm{\Gamma}_{n}=\left[\left(1+\dfrac{(\nu+1)H_{lq}}{b_{12}}\right)\left(1-\dfrac{H_{lq}}{b_{12}}\right)_{+}^{\nu+1}\right]_{l=q=1}^{n}$  and $\bm{\theta}=(\sigma_{1}^{2},\sigma_{2}^{2},\rho_{12},b_{12})^{\top}$. Taking derivatives, we obtain

[TABLE]

where $\bm{S}_{n}$ is given by

[TABLE]

For any matrix norm, the spectral radius $\lambda_{\text{max}}\{\bm{A}\}$ of an $n\times n$ matrix $\bm{A}$ satisfies $\lambda_{\text{max}}\{\bm{A}\}\leq\|\bm{A}\|$. Then, considering the norm $\|\cdot\|_{\infty}$, we have

[TABLE]

One can check that

[TABLE]

Thus $\sup_{n}\lambda_{\text{max}}\{\bm{\Gamma}_{n}\}<\infty$, which implies that $\sup_{n}\lambda_{\text{max}}\{\bm{\Sigma}_{n}(\bm{\theta})\}<\infty$. Because $\bm{\Gamma}_{n}$ is positive definite, $\lambda_{i}\{\bm{\Gamma}_{n}\}>0$, $i=1,...,n$. In particular, $\lambda_{\text{min}}\{\bm{\Gamma}_{n}\}>0$, so $\inf_{n}\lambda_{\text{min}}\{\bm{\Gamma}_{n}\}>0$ and $\inf_{n}\lambda_{\text{min}}\{\bm{\Sigma}_{n}\}>0$. Further,

[TABLE]

Similarly,

[TABLE]

Moreover, $\lambda_{\text{max}}\{\bm{S}_{n}\}\leq\|\bm{S}_{n}\|_{\infty}<\infty$ because of the form of the polynomial in $\bm{s}\in\mathbb{R}^{d}$ and the compact support in $b_{12}$. Then, for $\sigma_{1}^{2},\sigma_{2}^{2},\sigma_{1}\sigma_{2}\rho_{12}<\infty$,

[TABLE]

This implies that,

[TABLE]

The second derivatives are:

[TABLE]

where $\bm{SS}_{n}=\dfrac{\partial}{\partial b_{12}}\bm{S}_{n}$.

Because $\sup_{n}\lambda_{\text{max}}\{\bm{0}\}<\infty$, the compact support of $\bm{SS}_{n}$ in $b_{12}$, and the previous results, $\lambda_{\text{max}}\{\bm{SS}_{n}\}\leq\|\bm{SS}_{n}\|_{\infty}<\infty$. Then, for $\sigma_{1}^{2},\sigma_{2}^{2},\sigma_{1}\sigma_{2}<\infty$,

[TABLE]

In addition,

[TABLE]

This satisfies the first two conditions of Mardia and Marshall’s theorem.

For the third condition, we consider $\bm{A}=[a_{ij}]_{i=j=1}^{p}$, with $a_{ij}=\left\{\dfrac{t_{ij}}{(t_{mm}t_{nn})^{1/2}}\right\},$ and $t_{ij}=\text{tr}\left\{\bm{\Sigma}_{n}(\bm{\theta})^{-1}\dfrac{\partial\bm{\Sigma}_{n}(\bm{\theta})}{\partial\theta_{i}}\bm{\Sigma}_{n}(\bm{\theta})^{-1}\dfrac{\partial\bm{\Sigma}_{n}(\bm{\theta})}{\partial\theta_{j}}\right\}\ \text{for all}\ i,j=1,...,p$; we prove that $\bm{A}$ is non singular.

Notice that

[TABLE]

with $\bm{A}_{n}=\{\bm{\Gamma}_{n}^{-1}\circ\bm{S}_{n}\}$ where the operator $\circ$ denotes the matrix Hadamard product.

Then,

[TABLE]

For matrix $\bm{A}$ in equation (A.2),we have extended the result established by Bevilaqua et al. [2015]. Thus $\bm{A}$ is positive definite. By Mardia and Marshall’s Theorem the ML estimator of $\bm{\theta}=(\sigma_{1}^{2},\sigma_{2}^{2},\rho_{12},b_{12})^{\top}$ is asymptotically normal with variance $\bm{F}_{n}(\bm{\theta})^{-1}$.

Equation (7) implies that

[TABLE]

Fixing $\nu>0$, noting that $g(\cdot)$ is a continuously differentiable function for $\sigma_{1}\neq$ 0 and $\sigma_{2}\neq 0$, and using the multivariate delta method for $g(\cdot)$ we obtain

[TABLE]

where

[TABLE]

and $f(b_{12})=\left(-\dfrac{(\nu+1)\|\bm{h}\|}{b_{12}^{2}}\right)\left(1-\dfrac{\|\bm{h}\|}{b_{12}}\right)^{\nu+1}_{+}+\left(1+\dfrac{(\nu+1)\|\bm{h}\|}{b_{12}}\right)\left(1-\dfrac{\|\bm{h}\|}{b_{12}}\right)^{\nu}_{+}\dfrac{(\nu+1)\|\bm{h}\|}{b_{12}^{2}}.$ $\square$