Vectorization of non-recursive equation made difficult by NaN

I am performing some simple manipulations in a 2-dimensional grid where certain nodes should not be taken into account. These are specified in a similar sized matrix with NaN on to be excluded nodes.

The calculations were originally performed in a double loop. A very simplified version would look like this:

 new = zeros(5000,2000);

 for j=2:size(new,2)-1
   for i=2:size(new,1)-1
     if(~isnan(nanmat(i,j)))
       if(isnan(nanmat(i-1,j))), old(i-1,j)=old(i,j); end  
       if(isnan(nanmat(i+1,j))), old(i+1,j)=old(i,j); end
       if(isnan(nanmat(i,j-1))), old(i,j-1)=old(i,j); end
       if(isnan(nanmat(i,j+1))), old(i,j+1)=old(i,j); end
      
       new(i,j) = old(i,j) + (old(i+1,j) + old(i-1,j) + old(i,j+1) + old(i,j-1))/4;  
     end
   end
 end

Not really elegant...

So I thought it wouldn't be too hard to eliminate these loops, which is true for most of the code...:

 new = zeros(5000,2000);

 is = 2:size(new,1)-1;
 js = 2:size(new,2)-1;
 notNaN = ~isnan(nanmat);

 new(is,js) = old(is,js) + notNaN*(old(is+1,js) + old(is-1,js) + old(is,js+1) + old(is,js-1))/4;

the multiplication with notNaN compensates for the ~isnan check, but I cannot really seem to find a solution to account for the if-statements that check neighboring nodes. Is there an easy way to do so? These calculations are repeated very often and vectorizing this part of my code saves me about 15-20% in calculation time.

Thanks!

Answer 1

[Niels]

I seem to have found a solution that does the job by determining multiplication factors prior to the calculations, i.e.;

 new(is,js) = old(is,js) + notNaN(is,js).*(A0.*old(is,js) + A1.*old(is+1,js) + A2.*old(is-1,js) + A3.*old(is,js+1) + A4.*old(is,js-1))/4;

For instance, if my notNaN matrix is

 0 0 0 0 0 0 0
 0 0 0 1 1 1 0
 0 0 1 1 0 1 0
 0 0 1 1 0 1 0
 0 1 1 1 1 1 0
 0 0 1 1 1 0 0
 0 0 0 0 0 0 0

I would get A1 to be

 0 1 1 1 0
 1 1 0 1 0
 1 1 0 1 0
 1 1 1 1 0
 1 1 1 0 0

Which is equal to notNaN(is+1,js).

Consequently, that would mean

 A1 = notNaN(is+1,js);
 A2 = notNaN(is-1,js);
 A3 = notNaN(is,js+1);
 A4 = notNaN(is,js-1);
 A0 = 4 - (A1+A2+A3+A4);

I could even eliminate the multiplication with notNaN(is,js) by performing that multiplication up front as well.