计算几何中的编程练习

序

本文参考了Dsp Tian的博客, 主要通过MATLAB程序解决了以下几个小问题：

• 构造简单多边形
• 射线法判断点与多边形关系
• 最小包围矩形/Ritter’s最小包围圆
• 多线段交点

构造简单多边形

输入：$n$边形顶点组成的矩阵($n\times 2$).

输出：排序后的顶点坐标并绘制相应的多边形.

步骤1 计算出所有顶点的中心位置;

步骤2 求每个顶点与中心点的极角;

步骤3 按极角对顶点坐标进行排序;

步骤4 将排序后的顶点顺次连接.

function p=creatPolygon(V)
% 绘制简单多边形
% 输入数据V是顶点组成的矩阵
n=size(V,1); % 顶点个数
center=mean(V);
angle=atan2(V(:,1)-center(1),V(:,2)-center(2)); % 每个点到坐标中心极角

V=[V,angle];
V=sortrows(V,3); % 按极角排序

p=V(:,1:2);
% p(n+1,:)=p(1,:); % 将第一个点添加到最后一个点后面, 方便连线
% plot(p(:,1),p(:,2))
end

n=4;
V=rand(n,2)
p=creatPolygon(V);

V =
    0.8147    0.6324
    0.9058    0.0975
    0.1270    0.2785
    0.9134    0.5469
p =

0.3922    0.0318   -2.7431
0.1712    0.0462   -2.1721
0.7431    0.7060    0.5204
0.6555    0.2769    1.5001
0.3922    0.0318   -2.7431

射线法判断点与多边形的关系

**基本原理：**从待判断点引出一条射线, 射线与多边形相交, 若交点为偶数, 则点不在多边形内; 若交点为奇数, 则点在多边形内.

输入：多边形各顶点组成的$n\times 2$的矩阵及$num$个数据点$(2\times num)$.

输出：多边形内部包含的数据点.

步骤1 遍历所有的数据点, 为每个点生成一个计数器, 从当前点向做竖直的射线;

步骤2 计算多边形各边对应的直线方程;

步骤3 判断射线与多边形各边的相交情况, 每有一边相交, 计数器增加一次;

步骤4 根据计数情况, 将交点个数为偶数的点判定为外部点, 用红点表示; 将交点个数为奇数的点判定为内部点, 用绿色星号表示, 并记录下来.

function [poly,Internal]=Raymethod(points,V)

polyn=length(V); % 顶点个数n
poly=creatPolygon(V);

%将第一个点再添加到末尾, 方便判断
polyn=polyn+1;
poly(polyn,:)=poly(1,:);

% 采样点/矢量个数
num=size(points,2);

figure;
plot(poly(:,1),poly(:,2),'b');
hold on
Internal=[];
for i=1:num
    flag=0;
    for j=2:polyn
        x1=poly(j-1,1);  % 多边形前后两个点
        y1=poly(j-1,2);
        x2=poly(j,1);
        y2=poly(j,2);
        
        % 多边形一条边直线
        k=(y1-y2)/(x1-x2);
        b=y1-k*x1;
        x=points(1,i); 
        y=k*x+b; % 从当前点做竖直向上的射线看与多边形交点
        
        if min([x1 x2])<= x && x<=max([x1 x2]) && ...
                min([y1 y2])<= y && y<=max([y1 y2]) && y>=points(2,i)
            flag=flag+1;
        end
    end
    
    if mod(flag,2)==0 % 偶数则在外部
        plot(points(1,i),points(2,i),'r.');
    else % 奇数在内部
        plot(points(1,i),points(2,i),'g*');
        Internal=[Internal,points(:,i)];
    end
end

end 

V=rand(4,2);
points=rand(500,2)';
[poly,Internal]=Raymethod(points,V)

最小包围矩形

输入：离散数据点.

输出：最小包围矩形的四个顶点坐标.

步骤1 求数据点的凸包;

(注：可调用MATLAB的凸包计算函数$k=covhull(X,Y)$返回凸包上数据点的索引)

步骤2 将凸包中两相邻的点连线作为矩形的一条边;

步骤3 寻找凸包上距离已知边最远的点, 过该点作平行线, 得到矩形的另一条平行边;

步骤4 将凸包上的点向已求得的边投影, 求得相距最远的两个投影点, 过该两点作垂直方向的直线作为矩形另外两条边;

步骤5 遍历凸包所有相邻两点, 重复步骤2-4, 将面积最小的矩形作为最终结果, 输出其四个点的坐标. 通常情况下, 最小包围矩形会过5个数据点.

function pbar=minBar(points)
% 最小包围矩形
% points输入数据点
% 输出最小包围矩形的四个顶点
% 调用凸包计算函数k=convhull(X,Y), 返回凸包上数据点的索引

n=length(points);

k=convhull(points(:,1),points(:,2));
K=length(k);

hull=points(k,:); %凸包上数据点索引

area=inf;

for i=2:K
    p1=hull(i-1,:); % 凸包上两个点
    p2=hull(i,:);
    
    k1=(p1(2)-p2(2))/(p1(1)-p2(1)); % 连接两个点的直线, 作为矩形的一条边
    b1=p1(2)-k1*p1(1);
    d=abs(hull(:,1)*k1-hull(:,2)+b1)/sqrt(k1^2+1); % 所有凸包上的点到k1,b1直线的距离

    [h,ind]=max(d); % 得到距离最大的点的索引ind,距离即为宽

    b2=hull(ind,2)-k1*hull(ind,1); % 与k1,b1直线平行的另一条边k1,b2

    k2=-1/k1; % 已求直线的垂线斜率


    b=hull(:,2)-k2*hull(:,1); % 过凸包所有点构成的k2,b直线系
    x1=-(b1-b)/(k1-k2); % 凸包上所有点在已求的第一条边上的投影
    y1=-(-b*k1+b1*k2)/(k1-k2);

    x2=-(b2-b)/(k1-k2);
    y2=-(-b*k1+b2*k2)/(k1-k2);

    [~,indmax1]=max(x1);  % 投影在第一条边上x方向最大与最小值
    [~,indmin1]=min(x1);


    [~,indmax2]=max(x2);  % 投影在第二条边上x方向最大与最小值
    [~,indmin2]=min(x2);

    w=sqrt((x1(indmax1)-x1(indmin1))^2+(y1(indmax1)-y1(indmin1))^2); % 矩形的长

    if area>=h*w
        area=h*w;
        pbar=[x1(indmax1) y1(indmax1);
              x2(indmax2) y2(indmax2);
              x2(indmin2) y2(indmin2);
              x1(indmin1) y1(indmin1)]; % 矩形四个角点
    end
end
    pbar(5,:)=pbar(1,:);

    plot(points(:,1),points(:,2),'.');
    hold on
    plot(pbar(:,1),pbar(:,2),'r');
    axis equal;
    
end

Ritter’s最小包围圆

输入：离散数据点.

输出：最小包围圆的圆心$(x_0,y_0)$及半径$r$.

步骤1 从离散数据点中随机选出两个点作为直径作初始圆;

步骤2 判断下一个点$p$是否在圆中, 如果在则继续判断, 否则进入步骤3;

步骤3 将点$p$作为新圆的一个边界点, 另一个边界点为距离点$p$最远的旧圆上的点, 使用这两个点作为直径构造新圆;

步骤4 返回步骤2, 直至遍历完所有的数据点.

function [x0,y0,r]=Ritter(p)
n=size(p,1);
p1=p(1,:);
p2=p(2,:);
r=sqrt((p1(1)-p2(1))^2+(p1(2)-p2(2))^2)/2;
cenp=(p1+p2)/2;

for i=3:n

    newp=p(i,:);
    d=sqrt((cenp(1)-newp(1))^2+(cenp(2)-newp(2))^2);
    if d>r
        r=(r+d)/2;
        cenp=cenp+(d-r)/d*(newp-cenp);
    end
end
figure;
plot(p(:,1),p(:,2),'r.');
hold on
x0=cenp(1);
y0=cenp(2);
theta=0:0.01:2*pi;
x=x0+r*cos(theta);
y=y0+r*sin(theta);
plot(x,y,'b-',x0,y0,'g*');
axis equal

多线段交点

基本思路：根据线段的端点求两条直线的交点, 再判断直线的交点是否在两条线段上.

function Mlsip(p)

% p是n*4的矩阵每一行表示线段左右端点的坐标(x1,y1,x2,y2)

n=size(p,1);
figure;
for i=1:n
    pbar=p(i,:);
    pbar=reshape(pbar,[2 2]);
    line(pbar(1,:),pbar(2,:));
end
hold on
for i=1:n-1
    p1=p(i,:);
    k1=(p1(2)-p1(4))/(p1(1)-p1(3));
    b1=p1(2)-k1*p1(1);
    for j=i+1:n
        p2=p(j,:);
	k2=(p2(2)-p2(4))/(p2(1)-p2(3));
	b2=p2(2)-k2*p2(1);
	
	% 求两直线交点
	x=-(b1-b2)/(k1-k2); 
	y=-(-b2*k1+b1*k2)/(k1-k2);
	
	% 判断交点是否在两线段上
	if min(p1(1),p1(3))<=x && x<= max(p1(1),p1(3)) && ...
	   min(p1(2),p1(4))<=y && y<= max(p1(2),p1(4)) && ...
	   min(p2(1),p2(3))<=x && y<= max(p2(1),p2(3)) && ...
	   min(p2(2),p2(4))<=y && y<=max(p2(2),p2(4))
		plot(x,y,'r*')
	end

    end
end