将线数据插值到网格matlab

Question

我有一个由 5 列组成的矩阵，第一列和第二列用于行的 x_start & y_start，第三列和第四列用于 x_end & y_end 第五个是-这条线的污染物浓度-
例如：

%  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);

Answer 1

我将从一个较小的域开始展示一个示例。这将允许显示中间矩阵，有助于理解代码在做什么。但是，该代码完全可以扩展到更大的域，甚至可以根据需要扩展到曲线而不是直线。

首先我需要定义域限制然后我创建一个网格。

%% // 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 table 中获取线坐标） .

然后我使用 interp1 到 re-interpolate 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)

有趣的部分来了。您会注意到，对于每个用红色圆圈突出显示的点，我应该照亮右侧的网格块 和左侧的网格块 （当然域边界除外）。所以我所要做的就是找到域的哪条网格线段与你的线相交。这可以通过比较网格点的 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) 其他地方。
因此，现在分配您的浓度很容易，只需将 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 表示线下方或线穿过的域块。不幸的是，我只想要 "traversed" 一个，而不是下面的那个，所以我必须只保留 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 块穿过更多 X 个块。如果您的线斜率 > 45 度，最好使用相同的方法但反转轴（重新插入您的线，例如您在 [=21] 处获得 X 值（xi1） =] 网格交叉点，并用 Xq 而不是 Yq 进行比较。如果你的曲线同时具有垂直和水平梯度（例如，如果我放大幅度，就像我上面的正弦波），那么你将不得不使用 both 方法（超过 X 和超过 Y），然后将结果 Z 合并为一个（类似于 Z= Zx | Zy）。

将线数据插值到网格matlab

Interpolate line data to grid matlab

matlab

interpolation

line

解释：

重要提示：