How to get line rectangle intersection segment?

Question

I want to find weight matrix for Algebraic reconstruction method. For this I have to find the line intersection with grid. I can find direct line intersection with line but I have to store the intersected line segment grid number wise. So suppose if in grid first square don't intersect with grid then put zero on first element of weight matrix.

Here code which I tried for line intersection:

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']

I have a line, defined by the parameters m, h, where y = m*x + h This line goes across a grid (i.e. pixels).

For each square (a, b) of the grid (i.e. the square [a, a+1]x[b, b+1]), I want to determine if the given line crosses this square or not, and if so, what is the length of the segment in the square so that I can construct the weight matrix which is essential for algebraic reconstruction method.

I do not know what it is that you want help with. Your code is drawing a line and a grid with gridsize 1 from 3 to 6. And returns the y-values at 3,4,5 and 6. What is it that you want as a final result and what is the exact problem? — The Minion, Sep 25 '14 at 08:42
@TheMinion I have a line, defined by the parameters m, h, where y = m*x + h This line goes across a grid (i.e. pixels). For each square (a, b) of the grid (ie the square [a, a+1] x [b, b+1]), I want to determine if the given line crosses this square or not, and if so, what is the length of the segment in the square. So that I can construct the weight matrix which is essential for algebraic reconstruction method. — Parth Patel, Sep 25 '14 at 13:00
@ParthPatel Consider changing the title of the question to "How to get line rectangle intersection segment in matlab". I answered your question below and gave an example of how to do it. — DontCareBear, Sep 12 '15 at 01:02
@DontCareBear THank you :) I had done different way :) But your way seems efficient:) — Parth Patel, Sep 13 '15 at 05:07
It's always nice to find ways to solve those problems without loops in matlab. You can use the code I gave you to any quad or triangle mesh. Have fun :) — DontCareBear, Sep 13 '15 at 05:31

score 1 · Accepted Answer · edited May 23 '17 at 12:03

Here's a nice way to intersect a line with grid of rectangles and getting the lengths of each of the intersection segments: I used the line line intersection from the pseudo code in the third answer from this 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');

This is the result (the lengths of the arrows in the image are the segment lengths)

How to get line rectangle intersection segment?

1 Answers1