How to find correlation of an image?

Question

There is an image A of fixed size 256*256*3 (RGB). The mathematical formula for covariance between two adjacent pixels values x,y in an image is popularly known to be:

cov(x,y) = 1/n summation from i = 1 to n of [E(x_i-E(x))(y_i-E(y))]

r_xy = cov(x,y) / (sqrt(D(x)*D(y)))

D(x) = 1/n summation  from i = 1 to n of square[(x_i - E(x))]

E(x) = 1/n summation from i = 1 to n of (x_i)

where r_xy is the correlation coefficients between two horizontally, vertically, and diagonally adjacent pixels of these two images.

Q1: How to do the above computation in MATLAB?

Q2: How to randomly select say 5000 pairs of two horizontally adjacent pixels from the images and then plot the distribution of these two horizontally adjacent pixels?

Answer 1

As is typical with image processing for truecolor RGB images, there are a couple of key issues to first address. I mentioned these in my answer to your other question involving a different image processing algorithm, but they bear repeating here:

• Figuring out how to deal with the third dimension: An RGB image is actually a set of three 2-D matrices (one each for the red, green, and blue color components of the pixels) concatenated along a third dimension. When performing pixel-wise operations, you have to decide whether you are going to perform the operations three times (i.e. once for each color plane) or whether you are going to somehow collapse the values along the third dimension (i.e. converting to a grayscale intensity image with functions like RGB2GRAY) to give you a set of 2-D image data to operate on.
• Be mindful of data types: Image data loaded into MATLAB is typically in the form of an unsigned 8-bit integer, but sometimes it can be an unsigned 16-bit integer or a double precision type. When dealing with integer types, conversion to double precision is usually desired before performing certain operations in order to avoid certain aspects of integer arithmetic such as round-off and saturation.

Alright, now that those formalities are out of the way, I see your problem above as comprising two steps. First, you have to select subsets of paired pixels from the image, such as all horizontally-paired pixels. Second, you have to apply the statistical formulae you have above. In the examples below, I'll assume the operations are being performed on the red (i.e. first) color plane of the matrix A:

1. Selecting subsets of paired pixels: Let's start with the set of unique horizontal pairings of pixels. If I select the pixels in the first column of A and place them in the subset x, then the horizontally adjacent pixels will be those in the second column of A, and these will be placed in the subset y. I can also add the pixels in the second column to the subset x, and the horizontally adjacent pixels in the third column would then be placed in the subset y. Repeating this for all columns in A, we can see that the pixels in columns 1 through 255 will be in subset x, and the pixels in columns 2 through 256 will be in the subset y. The matrix indexing would therefore look like this:

   x = A(:,1:end-1,1);  %# All rows and columns 1 through 255 from red plane
   y = A(:,2:end,1);    %# All rows and columns 2 through 256 from red plane
   

   Following similar logic as above, you can construct the entire set of unique vertical pairings of pixels in this fashion:

   x = A(1:end-1,:,1);  %# Rows 1 through 255 and all columns from red plane
   y = A(2:end,:,1);    %# Rows 2 through 256 and all columns from red plane
   

   And likewise for the set of unique diagonal pairings of pixels, where "diagonal" runs from top left to bottom right in the matrix:

   x = A(1:end-1,1:end-1,1);  %# All but the last row and column
   y = A(2:end,2:end,1);      %# All but the first row and column
   

   Or for "anti-diagonals", where "diagonal" runs from bottom left to top right in the matrix:

   x = A(2:end,1:end-1,1);  %# All but the first row and last column
   y = A(1:end-1,2:end,1);  %# All but the last row and first column
   

   Now, you can choose any one of these sets of x and y data to perform the statistical calculations you want for the red color plane. You can repeat the above substituting 2 or 3 for the last index in each line to get the calculation for the green and blue color planes, respectively.

2. Performing the statistical tests: This part is simple. There is already a built-in function CORRCOEF for computing the correlation coefficient in MATLAB. You may have to reshape the subsets of pixel values x and y into column vectors first using single-colon indexing:

   r_xy = corrcoef(x(:),y(:));
   

   Functions also exist for the other formulae as well: MEAN for E(x), VAR for D(x), and COV for cov(x,y).

In regard to your second question, you can first create x and y as I did above for all unique pairs of horizontally adjacent pixels, then create a vector with a random permutation of the integer indices into x and y using the function RANDPERM. Selecting the first 5000 entries of those randomly permuted indices will give you 5000 random indices into x and y:

randIndex = randperm(numel(x));  %# A random permutation of the integers
                                 %#   from 1 to numel(x)
randIndex = randIndex(1:5000);   %# Pick the first 5000 indices
xRand = x(randIndex);            %# 5000 random values from x
yRand = y(randIndex);            %# The corresponding 5000 values from y

This will give you your 5000 pairs of horizontally adjacent pixel values in x and y. However, it is unclear what you mean by "plot the distribution". I'm guessing you will either end up using the function HIST or perhaps the function SCATTER for this purpose.