问题描述
如图所示，一条任意直线AB（严格来说应该是线段）经过四边形网格，需要求出AB所经过的网格标号以及其在每个网格内的长度，即标号①至⑩的网格中的线段长度。
在达成以上要求的基础上，尽可能提高算法速度。
可能会有疑问的名词解释：
网格线：x = N（N∈Z）及 y = M（M∈Z）所表示的直线
格点：两条横纵方向网格线的交点
问题分析
首先把问题拆解一下，可以发现目标主要有两个：
1. 求直线经过的每个网格的标号
2. 求每个网格内的线段长度
现在将这两个问题分别求解，并存放于两个不同的数组中便于对照。
首先求是每个网格内线段的长度，这个比较简单，只要算出AB与所有网格线的交点坐标（若点A,B不在网格线上，则也包含A,B本身）就能计算了。假设A的坐标为（x,y），B的坐标为（m,n），那么AB经过的x轴和y轴的范围为：
X = ceil(x):floor(m); % m>x
Y = ceil(y):floor(n); % n>y
网格线其实就是函数表达式为x = X 或 y = Y 的直线，计算交点就是简单的方程组求解，假设斜率k存在，则
k = (P2(2)-P1(2))/(P2(1)-P1(1));
grid_Y = ((P2(2)-P1(2))*X+P2(1)*P1(2)-P2(2)*P1(1))/(P2(1)-P1(1)); %交点的y轴坐标
grid_X = ((P2(1)-P1(1))*Y-P2(1)*P1(2)+P2(2)*P1(1))/(P2(2)-P1(2)); %交点的x轴坐标
% 合并交点集并排序，便于之后网格点序号的计算
intersec_X = sort([X,grid_X]);
intersec_Y = sort([Y,grid_Y]);
*以上代码中P1即为A，P1=[x,y]；P2即为B，P2=[m,n]；懒得修改了，并且之后都使用P1,P2。
这样计算会有一个问题，就是当有某个点位于格点时，这个点会被重复计算，因为在沿着x轴和y轴分别求交点的时候这个点都被计算到了。一开始我是使用unique函数去除重复点的，但是后来发现由于浮点数计算精度问题，相等的点会变得“不相等”，比如说交点（3，4），在x = 3时计算得y = 4, 但在y = 4时x的计算结果却可能是2.99999999，这个时候使用unique函数就是无效的，尽管它们是同一个点。所以我现在的思路是找出2.99999999这样与整数十分接近的数值并直接处理（设一个小一点的阈值，就算计算结果真的是准确的2.99999999也影响不大）。
% find computation error
Error_X = diff(intersec_X); %数组差分
error_id_x = find(Error_X<1e-8); %找到十分接近的数值位置
if ~isempty(error_id_x)
intersec_X(error_id_x) = round(intersec_X(error_id_x)); % 将浮点数转为整数
intersec_Y(error_id_x) = round(intersec_Y(error_id_x));
intersec_X(error_id_x+1) = []; % 删去重复数值
intersec_Y(error_id_x+1) = [];
end
至此我们已经得到了直线与网格线的所有交点坐标，其横坐标好和纵坐标分别保存在数组intersec_X和数组intersec_Y中，然后用平面距离公式即可求出每个网格内的线段长度：
delta_X = diff(intersec_X);
delta_Y = diff(intersec_Y);
% 每个网格内的线段长度
distance = sqrt(delta_X.^2 + delta_Y.^2);
接下来解决求网格的标号的问题，在此过程中需要用到刚才算得的交点坐标。
首先可以发现一个规律：直线经过的网格数 = 直线与网格线的交点数（含端点） + 1
既然刚才我们已经算出了直线与网格线的所有交点，那么直线经过的网格数也是已知的了。然后我们沿着直线前进的方向分析一下直线从一个交点延伸到另一个交点的过程中，其所经网格序号会有什么样的变化。
假设直线的斜率k存在且0<k<1，即问题描述中图片所示，取其中任意两个相邻的交点，可以发现，当没有特殊情况时，直线所经过的下一网格的坐标总是按（+1，0）（0，+1）的方式交错增加的，比如假设网格②的坐标为（2,1），那么网格③的坐标就为（2，2），网格④的坐标则为（3，2）。到了网格⑤，⑤的坐标为（4,2），不符合刚才说的规律了，这是因为⑤就是一种“特殊情况”。可以看到，网格④左右的两个交点，其横坐标的差为1，这就导致了直线延伸到下一个网格时，坐标的增长会强制变为（+1，0），但经过这个网格之后，⑥（4，3），⑦（5，3）又回到了刚才所说的（+1，0）（0，+1）规律。还有一种特殊情况是交点刚好也是格点，这时下一网格坐标的增长就会强制变为（+1，+1）。
总的来说，就是在当前网格坐标已知的情况下，对下一网格坐标相对于当前网格坐标的增长情况分了以下三类（0<k<1的情况下）：
- 交点坐标的横坐标差值为1：对应网格坐标增长强制变为（+1,0）
- 交点恰好是格点：对应网格坐标增长强制变为（+1,+1）
- 不属于上述两种情况：网格坐标增长为（0,+1）（+1,0）交替循环
对于k的其余请况，也是分为同样的三类，只不过具体的增长数值会不同，这里就不赘述了。
回到代码部分，首先求出情况二，也就是交点恰好是格点时交点的位置：
% find grid intersection & loccation
gridPoint_X = intersect(X,grid_X); % grid_X是y轴整数坐标带入求得的x轴坐标，求其与x轴整数坐标的交集
f = @(x)find(intersec_X==x);
gridPoint_loc = arrayfun(f,gridPoint_X); % 在交点集中定位格点所在位置
因为直线经过的网格数 = 直线与网格线的交点数（含端点） + 1，在起始坐标已知的情况下，只需要开一个与distance（见上文代码，保存每个网格内的线段长度）长度相同的数组来记录网格坐标增长情况就可以了。
% record route
route = zeros(length(intersec_X),2);
然后判断一下起始坐标的增长，给route赋一个起始值。
% start point
if ~mod(intersec_X(1),1)
route(1,1) = 1;
end
if ~mod(intersec_Y(1),1)
route(1,2) = 1;
end
分类判断：
if k < 1
id = find(delta_X==1); % 情况一
route(id+1,1) = 1;
route(gridPoint_loc,:) = 1;% 情况二
temp = find(all(route==0,2));% 情况三
for i = 1:length(temp)
route(temp(i),1) = route(temp(i)-1,2);
route(temp(i),2) = route(temp(i)-1,1);
end
elseif k > 1
id = find(delta_Y==1);
route(id+1,2) = 1;
route(gridPoint_loc,:) = 1;
temp = find(all(route==0,2));
for i = 1:length(temp)
route(temp(i),1) = route(temp(i)-1,2);
route(temp(i),2) = route(temp(i)-1,1);
end
else % k = 1
if ~isempty(gridPoint_loc)
route(gridPoint_loc,:) = 1;
end
temp = find(all(route==0,2));
for i = 1:length(temp)
route(temp(i),1) = route(temp(i)-1,2);
route(temp(i),2) = route(temp(i)-1,1);
end
end
之后在已知起始坐标的情况下，累加route就可以得到直线所经过的每个网格的坐标，并且可以和之前保存的每个网格内的线段长度distance一一对应。
为了加快算法速度，我尽可能不使用for循环，而是替换成了矩阵运算。
经老师提点，求经过的网格坐标有更简单的方法，只要把交点所在的所有网格都遍历并保存，然后删去其中的重复点就可以了。不过当时我这个已经写完了就没有再去实现，应该是更快捷的方法，而且更容易理解 。
代码
仅限于k>=0的情况，k<0的我另写了一份，没整合，就不放了。
function grid = grid_compute(P1,P2,d,H) % only for k>=0
% both P1 and P2 >=0
% k does not exist;
if P2(1) == P1(1)
Y = ceil(P1(2)+1e-8):ceil(P2(2));
grid = zeros(d,H);
grid(ceil(P2(1)),Y) = 1;
if mod(P1(2),1)
grid(ceil(P2(1)),Y(1)) = Y(1)-P1(2);
end
if mod(P2(2),1)
grid(ceil(P2(1)),Y(end)) = P2(2)-Y(end)+1;
end
return;
end
% k = 0
if P2(2) == P1(2)
X = ceil(P1(1)+1e-8):ceil(P2(1));
grid = zeros(d,H);
grid(X,ceil(P2(2))) = 1;
if mod(P1(1),1)
grid(X(1),ceil(P2(2))) = X(1)-P1(1);
end
if mod(P2(1),1)
grid(X(end),ceil(P2(2))) = P2(1)-X(end)+1;
end
return;
end
X = ceil(P1(1)):floor(P2(1));
Y = ceil(P1(2)):floor(P2(2));
k = (P2(2)-P1(2))/(P2(1)-P1(1));
grid_Y = ((P2(2)-P1(2))*X+P2(1)*P1(2)-P2(2)*P1(1))/(P2(1)-P1(1));
grid_X = ((P2(1)-P1(1))*Y-P2(1)*P1(2)+P2(2)*P1(1))/(P2(2)-P1(2));
% intersection of x-axis and y-axis
intersec_X = sort([X,grid_X]);
intersec_Y = sort([Y,grid_Y]);
% find computation error
Error_X = diff(intersec_X);
error_id_x = find(Error_X<1e-8);
if ~isempty(error_id_x)
intersec_X(error_id_x) = round(intersec_X(error_id_x));
intersec_Y(error_id_x) = round(intersec_Y(error_id_x));
intersec_X(error_id_x+1) = [];
intersec_Y(error_id_x+1) = [];
end
delta_X = diff(intersec_X);
delta_Y = diff(intersec_Y);
% find grid intersection & loccation
gridPoint_X = intersect(X,grid_X);
% gridPoint_Y = intersect(Y,grid_Y);
f = @(x)find(intersec_X==x);
gridPoint_loc = arrayfun(f,gridPoint_X);
% record route, length+1 for easy computation
route = zeros(length(intersec_X),2);
% record grid
distance = sqrt(delta_X.^2 + delta_Y.^2);
% start point
if ~mod(intersec_X(1),1)
route(1,1) = 1;
end
if ~mod(intersec_Y(1),1)
route(1,2) = 1;
end
if k < 1
id = find(delta_X==1); % for delta_x = 1
route(id+1,1) = 1;
route(gridPoint_loc,:) = 1;% for grid point
temp = find(all(route==0,2));% rest of point
for i = 1:length(temp)
route(temp(i),1) = route(temp(i)-1,2);
route(temp(i),2) = route(temp(i)-1,1);
end
elseif k > 1
id = find(delta_Y==1);
route(id+1,2) = 1;
route(gridPoint_loc,:) = 1;
temp = find(all(route==0,2));
for i = 1:length(temp)
route(temp(i),1) = route(temp(i)-1,2);
route(temp(i),2) = route(temp(i)-1,1);
end
else % k = 1
if ~isempty(gridPoint_loc)
route(gridPoint_loc,:) = 1;
end
temp = find(all(route==0,2));
for i = 1:length(temp)
route(temp(i),1) = route(temp(i)-1,2);
route(temp(i),2) = route(temp(i)-1,1);
end
end
% if P2 is not on the grid
if mod(P2(1),1) && mod(P2(2),1)
distance = [distance,sqrt((P2(1)-intersec_X(end))^2 + (P2(2)-intersec_Y(end))^2)];
else
route = route(1:end-1,:); % cut route
end
% draw route
grid = zeros(d,H);
% the first grid
start_x = ceil(intersec_X(1)+1e-8);
start_y = ceil(intersec_Y(1)+1e-8);
for i = 2:length(distance)
grid(start_x,start_y) = distance(i-1);
start_x = start_x+route(i,1);
start_y = start_y+route(i,2);
end
grid(start_x,start_y) = distance(end);
% if P1 is not on the grid
if mod(P1(1),1) && mod(P1(2),1)
grid(ceil(P1(1)),ceil(P1(2))) = sqrt((P1(1)-intersec_X(1))^2 + (P1(2)-intersec_Y(1))^2);
end
1万+
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
