# Cross variogram with a moving window

I need to generate cross variograms of images using moving windows. For that I use the following equation:

$$\gamma_{jk}(h)=\frac{1}{2n(h)}\sum_{i=1}^{n(h)}\Big\{\big[dn_j(x_i)-dn_j(x_i+h)\big]\cdot\big[dn_k(x_i)-dn_k(x_i+h)\big]\Big\}$$

The first part stands for one band(j) and next part of band k. To illustrate with sample matrices,

j =  1     2     3     4
5     6     7     8
9    10    11    12
13    14    15    16

k = 17    18    19    20
21    22    23    24
25    26    27    28
29    30    31    32


In actual case I am using 7 X 7 windows for large satellite images.

I also had to generate variograms from images for this work. For generation of variograms I had to consider only one band of data. For that case I used nlfilter for moving window and created a function to select and calculate values.

But for cross variograms, I am not able to decide upon what function to use. For this case the calculations will go like this:

(1 - 2)(17 - 18) + (2 - 3)(18 - 19)

and so on.

-
As cross-posting is not encouraged, it would be better to remove your question on SO if you feel you'll be more likely to get a response there. As your question was asked few hours ago, it is fair to wait a little bit before cross-posting, though. Anyway, let's wait and see. If you get an acceptable response here, don't forget to remove the SO duplicate. –  chl Apr 25 '11 at 9:23
@chl definitely I will do that. I will remove one post once I get a favorable response in either of the sites. –  Chethan S. Apr 25 '11 at 9:31
As promised I have deleted this question from SO. Therefore deleted the first part in the question now - "I have also posted this question in StackOverflow...." –  Chethan S. Apr 25 '11 at 17:12
Thanks. –  chl Apr 25 '11 at 17:23

I prefer a slight change of notation due to the many $n$'s appearing in the original. Let $\alpha$ and $\beta$ designate the images. Let $i$ and $j$ each designate pairs of indexes into the image rows and columns. (Indexing goes from $1$ to $m$ for rows and $1$ to $n$ for columns.) Let $h$ designate a relative index pair (so that its two entries are integers, either of which can be negative), also known as an offset. Then, by definition, the value of the experimental cross-variogram of these images at an offset $h$ is

$$\gamma_{\alpha,\beta}(h)=\frac{1}{2n(h)}\sum_{i}\left(\alpha[i+h]-\alpha[i]\right)\left(\beta[i+h] - \beta[i]\right).$$

The sum ranges over all indexes $i$ for which both $i$ and $i+h$ are valid indexes into both images; $n(h)$ is the number of such indexes (easily computed in the same way by taking a similar sum of $1$'s).

By expanding the summand algebraically the calculation is reduced to the problem of obtaining

$$\sum_{i}\alpha[i+h]\beta[i]$$

for various $h$, both positive and negative, ranging from $(1-m,1-n)$ through $(m-1,n-1)$.

Let us say that the reversal of an image negates the indexes; that is, the value of the reversal of $\alpha$ at the pixel $h$ is the value of $\alpha$ at $(m+1,n+1)-h$.

Such a sum can be seen as the reversal of the convolution of the reversal of $\alpha$ with $\beta$. It is best computed using discrete Fourier transforms after first padding each image to the right and down with zeros. The padding must extend to the range of the largest $h$ for which $\gamma$ needs to be computed. Convolutions with Fourier transforms are obtained by taking the inverse Fourier transform (itself a scalar multiple of the FT) of the product of the FTs.

Direct computation of the variogram via its definition for a pair of $m$ by $n$ images requires up to $m n$ products and sums for each value of $h$. Typically $O(m n)$ values of $h$ are needed. The direct algorithm therefore has $O(m^2 n^2)$ computational cost, which is ridiculously large for moderate (megapixel) images. The discrete Fourier transform costs at most $O(2m 2n \log(2m 2n))$ (assuming the maximum range of offsets $h$) and has to be applied only a constant number of times (3). The reversals and paddings cost $O(2m 2n)$. Thus the total cost is still only $O(12 m n \log(4 m n))$, a huge improvement.

As a simple example, take $\alpha$ and $\beta$ to be the matrices

1 2
3 4


and

5 6
7 8


After padding with zeros to the right and down (by two columns and two rows) and reversing $\alpha$, multiplying these two 4 by 4 matrices componentwise, and taking the inverse Fourier transform, we get

 8 0  14 23
0 0   0  0
18 0  20 39
30 0  38 70


Rotating this right by 2 columns and 2 rows and reversing gives

0  0  0  0
0  8 23 14
0 30 70 38
0 18 39 20


If you think of the new row and column indexes ranging $-2, -1, 0, 1$, this new matrix is exactly $\sum_{i}\alpha[i+h]\beta[i]$ (indexed by $h$). For example, the $h = (0,1)$ entry is 38 and indeed

$$\alpha[1,2]\beta[1,1] + \alpha[2,2]\beta[2,1] = 2 \times 5 + 4 \times 7 = 38.$$

-
can you please elaborate on how \sum_{i}\alpha[i+h]\beta[i] is obtained. I am not able to understand how to get it! Also it would greatly help me if you can point me to a resource which explains the above concepts in a simple manner. –  Chethan S. Apr 26 '11 at 10:07
@Chethan The last example shows explicitly how this sum is computed. I learned about the FFT method from the geostatistics literature long ago, but it's a stretch to characterize this as "simple." You can find an explicit description by Denis Marcotte in a 1996 paper in Computers & Geosciences. A 2010 PowerPoint presentation by Jim Jennings summarizes the approach (starting p. 14). –  whuber Apr 26 '11 at 14:33