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.

Answer 1

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)