mask a matrix based on values in two arrays

I am looking for a more effective means of masking a matrix based on values held in two arrays. Ultimately I'm looking to mask off chunks of a surface to alter their values. I have been able to achieve this based on the logic in the following loop, but for my actual model runs, which typically have meshgrid sizes in the millions of elements, and array lengths in the 1000's I feel like there could be a less memory intensive (and read/write intensive) solution.

[x y]=meshgrid(0:1:100,0:1:100);
u=floor(rand(10,1)*100);
v=floor(rand(10,1)*100);

padu=5; padv=5;

for i=1:size(u,1)
    if i==1
        mask=1*(x>u(i)-padu & x<u(i)+padu & y>v(i)-padv & y<v(i)+padv);
    else
        maskh=(x>u(i)-padu & x<u(i)+padu & y>v(i)-padv & y<v(i)+padv);
        mask(maskh==1)=1;
    end
end

surf(x,y,mask)

is there a more effective solution? I've searched but perhaps using the wrong terms, as I don't have an answer.

Thank you.

[Jan]

It would be useful if you write down in words, how the mask should be created. Of course Matlab can parse this easily from the code, but even experiences programmers can use some help by reading the idea of the author.

[oran]

Thanks Jan. So the basis of my question is this:

I am working to modify topography for large construction sites. The x & y coordinates are in the state plane coordinate system - think 10e7 sized numbers in x&y that do not match - i.e., x may be in the 3e7 range and y might be in 5e8 range, etc. Throughout the construction site I have 100's to 1000's of individual points that I need to mask off (including the buffer 'padu, padv') so that the resulting mask allows me to overwrite the unmodified topography with the modified topography that I've constructed - contained in a scatteredInterpolant.

I can do this with relative success based on the format of the code that I included, but I am taking a performance hit (I believe) as I am writing the mask and maskh variable, which often contain 10e8 or so values once per loop, on a loop that may have 10e5 elements.

Let me know if this is more clear.

I appreciate your help.

Answer 1

[Jan]

n = 1e4;
u = floor(rand(10,1)*100);
v = floor(rand(10,1)*100);

tic  % Original version:
[x y] = meshgrid(linspace(0, 100, n), linspace(0, 100, n));
padu  = 5;
padv  = 5;
for i=1:size(u,1)
    if i==1
        mask=1*(x>u(i)-padu & x<u(i)+padu & y>v(i)-padv & y<v(i)+padv);
    else
        maskh=(x>u(i)-padu & x<u(i)+padu & y>v(i)-padv & y<v(i)+padv);
        mask(maskh==1)=1;
    end
end
toc

Elapsed time is 3.186382 seconds.

tic;  % Version 1:
x = linspace(0, 100, n);   % Avoid MESHGRID, because if creates huge
y = linspace(0, 100, n).'; % redundant matrices
padu  = 5;
padv  = 5;
% Use a logical array as mask instead of a double array:
mask2 = (x>u(1)-padu & x<u(1)+padu) & (y>v(1)-padv & y<v(1)+padv);
on    = true;  % Slightly faster to use a constant instead of a function
for i = 2:size(u,1)
    % Avoid check of i==1 in each iteration my moving this exception
    % out of the loop.
    maskh = (x>u(i)-padu & x<u(i)+padu) & (y>v(i)-padv & y<v(i)+padv);
    mask2(maskh) = on;
end
toc

Elapsed time is 1.137238 seconds.

assert(isequal(mask, mask2), 'Different results!')

It is faster to call

mask = or(mask, newmask)

then

newmask = xy;

mask(newmask) = true;

See:

tic;  % Version 2
x = linspace(0, 100, n).';
y = linspace(0, 100, n).';
padu = 5;
padv = 5;
xm   = (x > u.'-padu & x < u.' + padu);  % Or: abs(x - u.') < padu;
ym   = (y > v.'-padv & y < v.' + padv);  % Or: abs(y - v.') < padv;
on   = true;
mask3 = false(n, n);
for i = 1:size(u, 1)
    mask3 = mask3 | (xm(:, i).' & ym(:, i));
end
toc

Elapsed time is 0.652382 seconds.

assert(isequal(mask, mask3), 'Different results!')

In my measurements this has almost the same speed:

tic;  % Version 3
x = linspace(0, 100, n);   % Avoid MESHGRID, because if creates huge
y = linspace(0, 100, n).'; % redundant matrices
padu  = 5;
padv  = 5;
mask4 = false;
for i = 1:size(u,1)
    mask4 = mask4 | (abs(x-u(i)) < padu) & (abs(y-v(i)) < padv);
    % Or:
    % mask4 = mask4 | (x>u(i)-padu & x<u(i)+padu) & (y>v(i)-padv & y<v(i)+padv);
end
toc

Elapsed time is 0.592103 seconds.

assert(isequal(mask, mask4), 'Different results!')

Check the speed with the original data on your local machine. The timings from the forum seem to be strange sometimes.

[oran]

Thank you so much - this gives me a lot of stuff to try. I appreciate it! I'll report back with what works best.

Answer 2

[Catalytic]

A=accumarray([u,v],1,[100,100]);

k=ones(2*padu,2*padv);

mask=conv2(A,k,'same');