将线数据插值到网格Matlab [英] Interpolate line data to grid matlab
问题描述
例如:
% 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];
如果我假设我有一个感兴趣的域10mx10m，并且此域除以有限差分像元大小1mx1m(即，域是10个像元乘以10个像元)
我想对上述线值进行插值，以便每个单元格与某条线相交可以获取线的浓度值
这是一个例子:
(来源: 0zz0.com )
用于离散域的示意图预览
该行具有浓度值，并且相交的单元格应具有相同的值
我知道matlab中的此类功能可以做到这一点，但不幸的是，对于散布点，例如
interp2 & griddata
我为这个问题尝试了一个小代码，但是它仍然会给出错误，并且interp2函数被设计用于散点插值(点而不是线)
[X,Y] = meshgrid(0:10);
conc=interp2(frac(:,1),frac(:,2),frac(:,5),X,Y);
我将显示一个以较小域开头的示例.这将允许显示中间矩阵，这对于理解代码在做什么很有用.但是，该代码可以完全扩展到更大的域，甚至在必要时甚至可以扩展到曲线而不是直线.
首先，我需要定义域限制，然后创建一个网格.
%% // 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 ) ;
到目前为止，没有什么太过异国情调了.然后，我将定义一条任意线(您将不得不从您的 frac 表中获取该线的坐标).
我使用interp1重新插入X网格上线的Y值，以使点与Y网格线相交.
%% // 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 ) ;
到目前为止，我们有:
%% // 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)
现在来了有趣的部分.您会注意到，对于用红色圆圈突出显示的每个点，我应该分别照亮右侧和左侧的网格块(当然，在域边界上除外).因此，我要做的就是找到您的线与该域相交的网格线 segments .通过比较网格点的Y坐标(在上面定义的矩阵Yq中具有该坐标)和与它们相交的线的Y坐标(即yi1)，可以轻松地做到这一点.所以我们继续:
%// 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) ;
*答案末尾将给出此处发生的详细情况.
现在您有了一个逻辑矩阵Z，即域的大小，该逻辑矩阵包含该行遍历的每个块的1(true)，以及其他所有地方的0(false). br>
因此，现在可以轻松地分配浓度，只需将Z乘以所需的浓度值就可以了:
%% // 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')
等一下:
正如我之前所说，这是完全可扩展的.它将适用于更大的域大小...因为不是那么直线:
说明:
(域大小为5个块)
>> 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
但是，matlab坐标系是从上到下的，而矩阵索引是从上到下的，因此，为了代表矩阵值，就像您看到的图片一样，我将矩阵垂直翻转(尽管只是为了显示，但不要这样做).实际计算). 所以:
>> 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
所有这些ones代表该行以下或遍历的域块.不幸的是，我只想要遍历"的一个，而不是下面的一个，因此我只需要保留1和0之间的接口.这是通过将矩阵移动一行并应用xor过滤器来完成的(再次，实际上不使用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
并且由于我们在上面说过，每个交叉点的左和右图块都必须突出显示，因此我们也将模式复制到左侧的一列:
>> 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
我们完成了:-).将该值乘以所需的值，只有突出显示的单元格才会有一个值，其余的将保持为0.
重要提示:
到目前为止，我仅使用Y交集.由于Y渐变是平滑的(角度为<45度的线，它与Y块的交点总是比Y块多)，因此效果很好.
如果线斜率> 45度，则最好使用相同的方法，但将轴反转(重新插入线，例如在Y网格交点处获得X值(xi1)，并且用Xq而不是Yq进行比较.
如果您的曲线同时具有垂直的和梯度(例如，如果放大幅度，例如上面的正弦波)，那么您将必须同时使用两者方法( X和over Y)，然后将生成的Z合并为一个(类似于Z= Zx | Zy).
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);
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).
这篇关于将线数据插值到网格Matlab的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！
