Matlab研究小问题：如何计算一条线段所经过的网格区域和各区域内的长度

问题描述

在一个10*10的网格里，判断一条线段经过的网格位置（判断序号）。并计算经过的每个网格内线段的长度。

出于问题简单化的需要，当线段正好在网格线上时，不计算经过的网格和长度。

分析

这个问题的复杂度在于，需要计算线段与网格线的交点，进而判断线段经过网格的序号。由于线段是以两个端点的形式给出，除了垂直于横轴以外，均可以通过斜截式来表示直线。因此需要对这两种情况分别进行处理，加上不计算网格和长度的情形，一共需要分三种情况进行处理。
设计一个calLength函数来进行相关计算，设计segs结构体来存储结果。
segs.index —— 经过网格的序号；
segs.index_x —— 经过网格的横轴序号；
segs.index_y —— 经过网格的纵轴序号；
segs.length —— 经过网格内的线段长度。
分段结果以结构数组形式存储。

代码

function segs = calLength(P1,P2)

segs.length = [];
segs.index_x = [];
segs.index_y = [];
segs.index = [];

P1_x = P1(1);
P1_y = P1(2);
P2_x = P2(1);
P2_y = P2(2);

xmin = min(P1_x,P2_x);
ymin = min(P1_y,P2_y);
xmax = max(P1_x,P2_x);
ymax = max(P1_y,P2_y);

if (P1_x == P2_x) && (round(P1_x) == P1_x)...
        || (P1_y == P2_y) && (round(P1_y) == P1_y)
    return;
end

pos = reshape(1:100,10,10);

if (P1_x == P2_x) && (round(P1_x) ~= P1_x)
    SP = unique([ymin,ymax,ceil(ymin):floor(ymax)]);
    for t = 1 : size(SP,2)-1
        segs(t).length = SP(t+1)-SP(t);
        segs(t).index_x = ceil(P1_x);
        segs(t).index_y = max(ceil(SP(t+1)),ceil(SP(t)));
        segs(t).index = pos(segs(t).index_x, segs(t).index_y);
    end
end

if (P1_x ~= P2_x)
    K = polyfit([P1_x,P2_x],[P1_y,P2_y],1);
    xpx = []; xpy = [];
    for i = ceil(xmin):floor(xmax)
        xpx(i-ceil(xmin)+1) = i;
        xpy(i-ceil(xmin)+1) = K(1)*i+K(2);
    end

    ypx = []; ypy = [];
    for j = ceil(ymin):floor(ymax)
        ypy(j-ceil(ymin)+1) = j;
        syms x;
        ypx(j-ceil(ymin)+1) = double(solve(K(1)*x + K(2) - j, x));
    end

    SP = unique([P1_x,P2_x,xpx,ypx;P1_y,P2_y,xpy,ypy]','rows');

    L = @(x) sqrt((SP(x+1,1) - SP(x,1)).^2 + (SP(x+1,2) - SP(x,2)).^2);

    for t = 1 : size(SP,1)-1
        segs(t).length = L(t);
        segs(t).index_x = max(ceil(SP(t+1,1)),ceil(SP(t,1)));
        segs(t).index_y = max(ceil(SP(t+1,2)),ceil(SP(t,2)));
        segs(t).index = pos(segs(t).index_x, segs(t).index_y);
    end
end

n = [];
for i = 1 : size(segs,2)
    if segs(i).length < eps
        n = [n i];
    end
end
segs(n) = [];

clear;clc;clf;
x = linspace(0,10,11);
y = linspace(0,10,11);
[X,Y] = meshgrid(x,y);
line(X,Y,'color','b');
line(X',Y','color','b');
axis equal;
axis([0 10 0 10]);
set(gca,'xtick',0:10);

gridindex = reshape(1:100,10,10)';
numposx = 0.5*(X(1:end-1,2:end)+X(1:end-1,1:end-1))-0.1;
numposy = 0.5*(Y(2:end,1:end-1)+Y(1:end-1,1:end-1));
for i = 1 : 10
    for j = 1 : 10
        text(numposx(i,j),numposy(i,j),num2str(gridindex(i,j)));
    end
end

P1 = input('P1=');
P2 = input('P2=');
segs = calLength(P1,P2);
line([P1(1) P2(1)],[P1(2) P2(2)],'color','r');

display('所经过的网格序号\长度分别为:');
for i = 1 : size(segs,2)
    display(['序号: ' num2str(segs(i).index)]);
    display(['长度: ' num2str(segs(i).length)]);
end

示例结果

P1=[8.7,9]
P2=[3.5,2]
所经过的网格序号\长度分别为:
序号: 24
长度: 0.83847
序号: 25
长度: 0.40726
序号: 35
长度: 1.2457
序号: 45
长度: 0.023956
序号: 46
长度: 1.2218
序号: 56
长度: 0.45517
序号: 57
长度: 0.79056
序号: 67
长度: 0.88638
序号: 68
长度: 0.35934
序号: 78
长度: 1.2457
序号: 88
长度: 0.071869
序号: 89
长度: 1.1739

经核对结果是正确的。
注：本例中对特殊情况也进行了测试，详细结果限于篇幅不在此做验证。