如何获取线矩形相交的线段？

Question

我希望找到代数重建方法的权重矩阵。为此，我必须找到与网格的交点。我可以找到直线与直线的交点，但我必须按网格编号存储相交的线段。因此，假设第一个方格不与网格相交，则在权重矩阵的第一个元素上放置零。

这是我尝试进行线交叉的代码：

ak = 3:6
aka = 3:6
x = zeros(size(aka))
y = zeros(size(ak))
for k = 1:length(ak)
  line([ak(1) ak(end)], [aka(k) aka(k)],'color','r')
end

% Vertical grid
for k = 1:length(aka)
  line([ak(k) ak(k)], [aka(1) aka(end)],'color','r')
end
hold on;
 X =[0 15.5]
 Y = [2.5 8.5] 
 m = (Y(2)-Y(1))/(X(2)-X(1)) ;
 c = 2.5 ; 
 plot(X,Y)
axis([0 10 0 10])
axis square
% plotting y intercept
for i = 1:4
    y(i) = m * ak(i) + c
    if y(i)<2 || y(i)>6
        y(i) = 0
    end
end
% plotting x intercept
for i = 1:4
   x(i) = (y(i) - c)/m 
    if x(i)<2 || x(i)>6
        x(i) = 0
    end
end  
z = [x' y']

我有一条由参数m，h定义的线，其中y = m * x + h。这条线横跨网格（即像素）。

对于网格的每个方块(a, b)（即方块[a, a+1]x[b, b+1]），我想确定给定的线是否横跨此方块，并且如果是，则该段在方块中的长度是多少，以便我可以构建权重矩阵，这对于代数重建方法至关重要。

- Parth Patel

---

我不知道你需要什么帮助。你的代码画了一条线和一个网格，网格大小为1，从3到6。并返回3、4、5和6处的y值。你想要什么样的最终结果，确切的问题是什么？ - The Minion

@TheMinion 我有一条线，由参数m、h定义，其中y = m*x + h。这条线穿过一个网格（即像素）。对于网格（即[a，a+1] x [b，b+1]的正方形）的每个正方形（a，b），我想确定给定的线是否穿过该正方形，如果是，则该线段在正方形中的长度是多少。这样我就可以构建权重矩阵，这对于代数重建方法至关重要。 - Parth Patel

@ParthPatel 考虑将问题的标题更改为“如何在matlab中获取线矩形相交线段”。我已经回答了你的问题，并给出了如何做的示例。 - DontCareBear

@DontCareBear 谢谢 :) 我用了不同的方法 :) 但是你的方法看起来更有效 :) - Parth Patel

在Matlab中找到不使用循环解决问题的方法总是很好的。你可以将我给你的代码应用于任何四边形或三角形网格。玩得开心 :) - DontCareBear

Answer 1

这是一种不错的方法，可以将一条线与矩形网格相交，并获取每个交点段的长度：我使用了第三个答案中伪代码的线线交点link。

% create some line form the equation y=mx+h
m = 0.5; h = 0.2;
x = -2:0.01:2;
y = m*x+h;
% create a grid on the range [-1,1]
[X,Y] = meshgrid(linspace(-1,1,10),linspace(-1,1,10));
% create a quad mesh on this range
fvc = surf2patch(X,Y,zeros(size(X)));
% extract topology
v = fvc.vertices(:,[1,2]);
f = fvc.faces;
% plot the grid and the line
patch(fvc,'EdgeColor','g','FaceColor','w'); hold on;
plot(x,y);
% use line line intersection from the link
DC = [f(:,[1,2]);f(:,[2,3]);f(:,[3,4]);f(:,[4,1])];
D = v(DC(:,1),:);
C = v(DC(:,2),:);
A = repmat([x(1),y(1)],size(DC,1),1);
B = repmat([x(end),y(end)],size(DC,1),1);
E = A-B;
F = D-C;
P = [-E(:,2),E(:,1)];
h = dot(A-C,P,2)./dot(F,P,2);
% calc intersections
idx = (0<=h & h<=1);
intersections = C(idx,:)+F(idx,:).*repmat(h(idx),1,2);
intersections = uniquetol(intersections,1e-8,'ByRows',true);
% sort by x axis values
[~,ii] = sort(intersections(:,1));
intersections = intersections(ii,:);
scatter(intersections(:,1),intersections(:,2));
% get segments lengths
directions = diff(intersections);
lengths = sqrt(sum(directions.^2,2));
directions = directions./repmat(sqrt(sum(directions.^2,2)),1,2);
directions = directions.*repmat(lengths,1,2);
quiver(intersections(1:end-1,1),intersections(1:end-1,2),directions(:,1),directions(:,2),'AutoScale','off','Color','k');

这是结果（图像中箭头的长度是线段长度）

- DontCareBear