Interpolate line data to grid matlab

Question

i have a matrix consist of 5 columns the first and second columns are for x_start & y_start of the line, the third and fourth are for x_end & y_end the fifth is -concentration of contaminant in this line-
for example:

%  x_start  y_start  x_end   y_end concentration
frac= [0         0      1       1     0.3
       1         1      3       3     0.6
       3         3      10      2     1.2
       3         3      10      8     0.5];

if i assume that i have a domain of interest 10mx10m and this domain is divided by finite difference cell size 1mx1m (i.e domain is 10 cells by 10 cells)
and i want to interpolate the above mentioned line values so that each cell intersects with a certain line can take line's concentration value

Here's an example:

_{(source: 0zz0.com)}

for schematic preview of the discretized domain
the line has a concentration value and the intersecting cells should take same value

i knew about such functions in matlab can do this but unfortunately for scattered points such as

   interp2 & griddata

i tried a small code for this problem but it still gives error and the interp2 function is designed for scatter interpolation (points not lines)

   [X,Y] = meshgrid(0:10);
   conc=interp2(frac(:,1),frac(:,2),frac(:,5),X,Y); 

Answer 1

I'll show an example with a smaller domain to start with. This will allow showing the intermediate matrices, useful to understand what the code is doing. However, the code is completely scalable to bigger domains and even to curves instead of lines if necessary.

First I need to define the domain limits Then I create a grid.

%% // Domain Parameters
domainProps.Xlim = [0 5] ; 
domainProps.Ylim = [0 5] ;
domainProps.nCellX = 5 ;       %// number of cells into the domain (X axis)
domainProps.nCellY = 5 ;       %// number of cells nito the domain (Y axis)

xd = linspace( domainProps.Xlim(1) , domainProps.Xlim(2) , domainProps.nCellX+1 ) ;
yd = linspace( domainProps.Ylim(1) , domainProps.Ylim(2) , domainProps.nCellY+1 ) ;
[Xq,Yq,Zq] = meshgrid( xd, yd, 0 ) ;

Nothing too exotic so far. Then I'll define an arbitrary line (you will have to take the line coordinates from your frac table).

I then use interp1 to re-interpolate the Y values of the line on the X grid to have the points intersecting the Y grid lines.

%% // example for first line
xl1 = [0 5] ;      %// X for line 1 - take these data from your table "frac"
yl1 = [1 3.8] ;    %// Y for line 1 - take these data from your table "frac"

%// reinterp the Y values over the X-Grid defining the domain
yi1 = interp1( xl1 , yl1 , xd ) ;

So so far we have:

%% // display what we've got so far
figure ; grid on ; hold on
hp = plot(xl1,yl1) ;
hpi = plot(xd , yi1,'or')
set(gca,'Xlim',domainProps.Xlim,'YLim',domainProps.Ylim,'XTick',xd,'YTick',yd)

Now comes the fun part. You'll notice that for each point highlighted by a red circle, I should illuminate the grid block on the right and the grid block on the left (except on the domain border of course). So all I have to do is to locate which grid line segments of the domain are intersected by your line. This can be easily done by comparing the Y coordinates of the grid points (we have that in the matrix Yq defined above) with the Y coordinates of the line intersecting them (this is yi1). So on we go:

%// find the elements of the grid lower or equal to the corresponding value of yi
Z = bsxfun(@le,Yq,yi1) ;
%// keep only the "boundary" values (between "0" and "1" domains)
Z(1:end-1,:) = xor( Z(1:end-1,:) , Z(2:end,:) ) ;
%// now replicate the pattern to the left (because each point at
%// intersection must illuminate both block on the left and on the right
Z(:,1:end-1) = Z(:,1:end-1) | Z(:,2:end) ;

*details of what is happening here will be given at the end of the answer.

Now you have a logical matrix Z, the size of your domain, which contains 1 (true) for every block traversed by the line, and 0 (false) everywhere else.
So to assign your concentration is now easy, simply multiply Z by the desired concentration value and you are done:

%% // Assign value to block traversed by the line
conc = Z .* 0.8 ;           %// use the concentration value from your table
figure
hpc = pcolor(Xq,Yq,conc) ;  %// plot the domain with the illuminated blocks

%% // cosmetic changes
set(hpc,'EdgeColor',[0.5 0.5 0.5] )
caxis([0 1])
colorbar
hold on
hpp = plot(xl1,yl1,':k','LineWidth',2) ;
hpi = plot(xd , yi1,'or')

et voila:

---

As I said earlier, this is completely scalable. It will work for larger domain size ... as for not so straight lines:

---

Explanations:

(with a domain size of 5 blocks)

>> Z = bsxfun(@le,Yq,yi1)
Z =
     1     1     1     1     1     1
     1     1     1     1     1     1
     0     0     1     1     1     1
     0     0     0     0     1     1
     0     0     0     0     0     0
     0     0     0     0     0     0

But matlab coordinate system is down to top, while the matrix indexing is top to bottom, so to represent the matrix values as you would see the picture, I flip the matrix vertically (just for display though, do not do that in the real calculations). So:

>> flipud( Z )
ans =
     0     0     0     0     0     0
     0     0     0     0     0     0
     0     0     0     0     1     1
     0     0     1     1     1     1
     1     1     1     1     1     1
     1     1     1     1     1     1

All these ones represent the blocks of the domain which are below or traversed by the line. Unfortunately, I only want the "traversed" one, not the one below, so I have to keep only the interface between the 1s and the 0s. This is done by shifting the matrix by one row and apply a xor filter (again, do not actually use the flipud):

>> Z(1:end-1,:) = xor( Z(1:end-1,:) , Z(2:end,:) ) ;
>> flipud(Z)
ans =
     0     0     0     0     0     0
     0     0     0     0     0     0
     0     0     0     0     1     1
     0     0     1     1     0     0
     1     1     0     0     0     0
     0     0     0     0     0     0

And since we said above that both left and right blocks of each intersection have to be highlighted, we also replicate the pattern one column to the left:

>> Z(:,1:end-1) = Z(:,1:end-1) | Z(:,2:end) ;
>> flipud(Z)
ans =
     0     0     0     0     0     0
     0     0     0     0     0     0
     0     0     0     1     1     1
     0     1     1     1     0     0
     1     1     0     0     0     0
     0     0     0     0     0     0

we're done :-). Multiply that by the value you want, only the highlighted cells will have a value, the rest will stay at 0.

---

Important note:

So far I only use the Y intersections. It worked very well because the Y gradient is smooth (the line having an angle < 45 degrees, it always crosses more X blocks than Y blocks. If your line slopes are > 45 degrees, it will be better to use the very same method but invert the axis (reinterpolate your line such as you obtain the X values (xi1) at the Y grid intersections, and do the comparison with Xq instead of Yq. If you have a curve which have both vertical and horizontal gradients (like my sine wave above if I amplify the magnitude for example), then you will have to use both methods (over X and overY), then combine the resulting Z into one (something like Z= Zx | Zy).