MatLab: Matrix with one peak and rest decreasing

Question

I'm trying to create a matrix such that if I define a random number between 0 and 1 and a random location in the matrix, I want all the values around that to "diffuse" out. Here's sort of an example:

0.214  0.432  0.531  0.631  0.593  0.642
0.389  0.467  0.587  0.723  0.654  0.689
0.421  0.523  0.743  0.812  0.765  0.754
0.543  0.612  0.732  0.843  0.889  0.743
0.322  0.543  0.661  0.732  0.643  0.694
0.221  0.321  0.492  0.643  0.521  0.598

if you notice, there's a peak at (4,5) = 0.889 and all the other numbers decrease as they move away from that peak.

I can't figure out a nice way to generate a code that does this. Any thoughts? I need to be able to generate this type of matrix with random peaks and a random rate of decrease...

Answer 1

I ended up creating a code that wraps the way I want based on a dimension, which I provide. Here's the code:

dims = 100;
A = zeros(dims);

b = floor(1+dims*rand(1));
c = floor(1+dims*rand(1));

d = rand(1);

x1 = c;
y1 = b;

A(x1,y1) = d;

for i = 1:dims
    for j = i
        k = 1-j;

        while k <= j
            if x1-j>0 && y1+k>0 && y1+k <= dims
                if A(x1-j,y1+k) == 0
                    A(x1-j,y1+k) = eqn(d,x1-j,y1+k,x1,y1);
                end
            end

            k = k+1;
        end

    end

    for k = i
        j = 1-k;
        while j<=k
            if x1+j>0 && y1+k>0 && y1+k <= dims && x1+j <= dims
                if A(x1+j,y1+k)==0
                    A(x1+j, y1+k) = eqn(d,x1+j,y1+k,x1,y1);
                end
            end

            j = j+1;
        end
    end

    for j = i
        k = 1-j;
        while k<=j
            if x1+j>0 && y1-k>0 && x1+j <= dims && y1-k<= dims
                if A(x1+j,y1-k) == 0
                    A(x1+j,y1-k) = eqn(d,x1+j,y1-k,x1,y1);
                end
            end

            k=k+1;
        end
    end

    for k = i
        j = 1-k;
        while j<=k
            if x1-j>0 && y1-k>0 && x1-j <= dims && y1-k<= dims
                if A(x1-j,y1-k)==0
                    A(x1-j,y1-k) = eqn(d,x1-j,y1-k,x1,y1);
                end
            end

            j = j+1;
        end
    end
end

colormap('hot');
imagesc(A);
colorbar;

If you notice, the code calls a function (I called it eqn), which provided the information for how to changes the values in each cell. The function that I settled on is d/distance (distance being computed using the standard distance formula).

It seems to work pretty well. I'm now just trying to develop a good way to have multiple peaks in the same square without one peak completely overwriting the other.

Answer 2

Without knowing what other constraints you want to implement:

1. Come up with a function z = f(x,y) whose peak value is at (x0,y0) == (0,0) and whose values range between [0,1]. As an example, the PDF for the Normal distribution with mu = 0 and sigma = 1/sqrt(2*pi) has a peak at x == 0 of 1.0, and whose lower bound is zero. Similarly, a bivariate normal PDF with mu = {0,0} and determinate(sigma) == [1/(2*pi)]^2 will have similar characteristics.
2. Any mathematical function may have its domain shifted: f(x-x0, y-y0)

Your code will look something like this:

someFunction = @(x,y) theFunctionYouPicked(x,y);

[x0,y0,peak] = %{ you supply these values %};
myFunction = @(x,y) peak * someFunction(x - x0, y - y0);

[dimX,dimY] = %{ you supply these values %};
mymatrix = bsxfun( myFunction, 0:dimX, (0:dimY)' );

You can read more about bsxfun here; however, here's an example of how it works:

bsxfun( blah, [a b c], [d e f]' )

That should give the following matrix (or its transpose ... I don't have matlab in front of me):

[blah(a,d) blah(a,e) blah(a,f);
 blah(b,d) blah(b,e) blah(b,f);
 blah(c,d) blah(c,e) blah(c,f)]

---

Get a toy example working, then you can tinker with it to be more flexible. If the function dictating how it decreases is random (with the constraint that points closer to (x0,y0) are larger than more distant points), it won't be an issue to make a procedural function instead of using strictly mathematical ones.

---

In response to your answer:

Your equation could be thought of as a model for gravity where an object instantaneously induces a force on another mass, then stops exerting force. Following that logic, it could be modified to a naive vector formulation like this:

% v1 & v2 are vectors that point from the two peak points to the point [ii,jj]
theMatrix(ii,jj) = norm(  (r1 / norm( v1 )) * v1 / norm( v1 ) ...
                        + (r2 / norm( v2 )) * v2 / norm( v2 ) ...
                       );

The most extreme type of corner case you'll run into is one where v1 & v2 point in the same direction as in the following row:

[ . . A X1 X2 . . ]

... where you want a value for A w/respect to X1 & X2. Using the above expression it'll boil down to A = X1 / norm(v1) + X2 / norm(v2), which will definitely exceed the peak value at X1 because norm(v1) == 1. You could certainly do some dirty stuff to Band-Aid it, but personally I'd start looking for a different function.

Along those lines, if you used Newton's Law of Universal Gravitation with a few modifications:

• You wouldn't need an analogue for G, so you could just assume G == 1
• Treat each of the points in the matrix as having mass m2 == 1, so the equation reduces to: F_12 == -1 * (m1 / r^2) * RHAT_12
• Sum the "force" vectors and calculate the norm to get each value

... you'll still run into the same problem. The corner case I laid out above would boil down to A = X1/norm(v1)^2 + X2/norm(v2)^2 == X1 + X2/4. Since it's inversely proportional to the square of the distances, it'd be easier to Band-Aid than the linear one, but I wouldn't recommend it.

Similarly, if you use polynomials it won't scale well; you can design one that won't ever exceed your chosen peaks, but there wouldn't be a lower bound.

You could use the logistic function to help with this:

1 / (1 + E^(-c*x))

Here's an example of using the logistic function on a degree 4 polynomial with peaks at points 2 & 4; you'll note I gave the polynomial a scaling factor to pull the polynomial down to relatively small values so calculated values aren't so close together.