Make index of elements based on random values in another matrix

I have the following two matrices, one of random numbers, and the other of a random integer 1-8:

 X = rand(10,10);
 r = round(0.5+rand(size(X))*8);

I want the numbers in matrix r to refer to adjacent elements in the matrix X as:

 1    2      3
 4    X(i,j) 5
 6    7      8

and return the indices of those elements. I will avoid repeating elements by only working on every 3rd row and column (1:3:end,1:3:end). I want to use periodic boundary conditions.

So, how can I return the indices of the randomly selected neighbor of each element X(1:3:end,1:3:end)? It should return in a matrix of size(X(1:3:end,1:3:end))

The reason I want to do this is because I want to flip the elements in X between X(i,j) and a neighboring X with a row/column offset determined by the value in matrix r. So if r(i,j)=7, then I can use the index to make a new matrix X_new where X_new(i,j)=X(i+1,j). This needs to be fast, so I am trying to vectorize.

Thank you!!!

Answer 1

[Guillaume]

What you want to do is transform your r values into linear offsets into X. These linear offsets would be

-size(X, 1)-1  -1   size(X,1)-1
-size(X, 1)    N/A  size(X,1)
-size(X, 1)+1   1   size(X,1)+1

So:

offsets = [-size(X,1)-1 -1 size(X,1)-1 -size(X,1) size(X,1) -size(X,1)+1 1 size(X,1)+1];
roffsets = offsets(r);

You then just have to add these offsets to the linear indices of X:

xindices = reshape(1:numel(X), size(X));
swapindices = xindices + roffsets;

swapindices defines the index of each element to replace the original x element, for the whole matrice. You could just do:

newX = X(swapindices);

to get the new matrice. There will be problems at the edges of the matrix since there's no restriction to swap with columns (rows) 0 and numcol (numrow) + 1. You haven't specified what you want to do at the edges.

Note: you can generate r simply with:

r = randi(8, size(X));

[Christopher]

Thanks much! I want the edges to be treated periodically--So left of the first column is the last column, above the first row is the last row, top left of the first element is the last element, etc.

Is there an easy way to incorporate that into your proposed linear index offset, or do I have to pad the matrix and only work on the interior? Thanks again.

[Guillaume]

There are probably many ways to deal with the boundary. You could indeed pad the matrix. Another way would be to deal with in the creation of the linear indices. For the left and right edge, it's easy it's just a modulo operation. For the top and bottom edge, I would extend the indices of r so [9 10 11] means swap with bottom edge and [12 13 14] swap with top edge:

r = randi(8, size(X));
rt = r(1, :); rt(rt <= 3) = rt(rt <= 3) + 8; r(1, :) = rt;
rb = r(size(r, 1), :); rb(rb >= 6) = rb(rb >= 6) + 6; r(size(r, 1), :) = rb;

You then extend the offsets matrix to:

offsets = [-size(X,1)-1 -1 size(X,1)-1 -size(X,1) size(X,1) -size(X,1)+1 1 size(X,1)+1 -1 size(X,1)-1 2*size(X,1)-1 -2*size(X,1)+1 -size(X,1)+1 1];

The roffsets and swapindices calculation is the same:

roffsets = offsets(r);
xindices = reshape(1:numel(X), size(X));
swapindices = xindices + roffsets;

But before you use the indices you pass them through modulo:

swapindices = mod(swapindices-1, numel(X)) + 1;
newX = X(swapindices);

Answer 2

[Image Analyst]

Not sure where the 7 comes into play in your last paragraph. And what are you doing? I didn't understand your algorithm. What do you do with the 3x3, or 7x7, window once you've extracted it? What I'm wondering is if you can use conv2(), imfilter(), or blockproc() to do what you want to do instead of manually scanning the array with the window yourself.

[Christopher]

7 indicates a reference to the element below (linear index +1). So if matrix X are my original values, I will switch the value of the randomly selected neighboring element in the matrix X. The commenter Guillaume has me on the right track. Now I just need to deal with boundaries.

I don't think I can use conv2 because there are 8 different kernels, not one. I'm not sure how the other functions work.

[Image Analyst]

OK that's fine. Glad he got it all working for you and you accepted his answer. I just didn't see why, when you said " So if r(i,j)=7, then I can use the index to make a new matrix X_new where X_new(i,j)=X(i+1,j)" that 7 did not show up in the indexes for X_new or X, and how they would change if the r value was 13 or 42 or 123 or whatever, since i and j didn't have anything at all to do with 7.

I still don't see that you've said what you're going to do with X_new. And since you will have a different X_new at every location, I don't see why you'd need to call conv2() 8 times. For example if you were just going to sum the values in X_new, you could simply do

out = conv2(in, ones(3), 'same');

and do that just one time to process all possible X_news in the entire image.