我将展示一个以较小域开头的示例。这将允许显示中间矩阵，有助于理解代码在做什么。但是，代码完全可以扩展到更大的域，如果需要，甚至可以扩展到曲线而不是直线。
首先我需要定义域限制然后我创建一个网格。
%% // 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)。
因此，现在分配您的注意力很容易，只需将 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 代表了在该行下方或经过该行的域的块。不幸的是，我只想要“遍历”的那个，而不是下面的那个，所以我只需要保留1s 和0s 之间的接口。这是通过将矩阵移动一行并应用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 渐变是平滑的（角度X 块而不是Y 块。
如果您的线斜率 > 45 度，最好使用相同的方法但反转轴（重新插入您的线，例如在Y 网格交叉点获得X 值（xi1），以及与Xq 而不是Yq 进行比较。
如果您有一条同时具有垂直 和 水平梯度的曲线（例如，如果我放大幅度，就像我上面的正弦波），那么您将不得不使用 both 方法（ X 和 Y)，然后将生成的 Z 合并为一个（类似于 Z= Zx | Zy）。