Matlab限定Voronoi输出的泰森多边形范围

Matlab提供的voronoi函数可以用来计算给定点集的泰森多边形

[vx,vy]=voronoi(...)

但是这一函数返回的泰森多边形线集存在一个明显的问题，当某一数据点在点集的边缘区域时，其泰森多边形可能不是一个闭合区域，甚至有一些边的长度比点集的范围尺度要大很多。

为了解决这个问题，许多网友的选择是用xlim或者axis此类函数去限定绘图区域，这种做法适用于只需要泰森多边形结果图片的人，像笔者这样需要对voronoi函数返回值进行进一步操作的人来说，没有真正改变线集矩阵的操作是不能满足要求的。

此处笔者通过一些简单的数学方法来将voronoi返回的线集进行修改，实现将泰森多变形束缚在指定的范围内，并添加一些边来使多边形闭合。

[voronoiX,voronouY]=voronoi(...)

以下是相关matlab代码，其中设置的边界范围为x∈(−300,900),y∈(−100,500)

xBound1=-300;xBound2=900;yBound1=-100;yBound2=500;[a b]=size(voronoiX);for indexI=1:b       %循环找出越界的边，并根据直线方程将其缩小到边界处    node_x1=voronoiX(1,indexI);    node_x2=voronoiX(2,indexI);    node_y1=voronoiY(1,indexI);    node_y2=voronoiY(2,indexI);    if node_x1<xBound1        indexY=(node_y2-node_y1)/(node_x2-node_x1)*(xBound1-node_x1)+node_y1;        if indexY<=yBound2&&indexY>=yBound1            voronoiX(1,indexI)=xBound1;            voronoiY(1,indexI)=indexY;        end    end    if node_x2<xBound1        indexY=(node_y1-node_y2)/(node_x1-node_x2)*(xBound1-node_x2)+node_y2;        if indexY<=yBound2&&indexY>=yBound1            voronoiX(2,indexI)=xBound1;            voronoiY(2,indexI)=indexY;        end    end    %%%%%%%%%%%%%    %%%%%%%%%%%%%    if node_x1>xBound2        indexY=(node_y2-node_y1)/(node_x2-node_x1)*(xBound2-node_x1)+node_y1;        if indexY<=yBound2&&indexY>=yBound1            voronoiX(1,indexI)=xBound2;            voronoiY(1,indexI)=indexY;        end    end    if node_x2>xBound2        indexY=(node_y1-node_y2)/(node_x1-node_x2)*(xBound2-node_x2)+node_y2;        if indexY<=yBound2&&indexY>=yBound1            voronoiX(2,indexI)=xBound2;            voronoiY(2,indexI)=indexY;        end    end    %%%%%%%%%%%%%    %%%%%%%%%%%%%    if node_y1<yBound1        indexX=(node_x2-node_x1)/(node_y2-node_y1)*(yBound1-node_y1)+node_x1;        if indexX<=xBound2&&indexX>=xBound1            voronoiX(1,indexI)=indexX;            voronoiY(1,indexI)=yBound1;        end    end    if node_y2<yBound1        indexX=(node_x1-node_x2)/(node_y1-node_y2)*(yBound1-node_y2)+node_x2;        if indexX<=xBound2&&indexX>=xBound1            voronoiX(2,indexI)=indexX;            voronoiY(2,indexI)=yBound1;        end    end    %%%%%%%%%%%%%    %%%%%%%%%%%%%    if node_y1>yBound2        indexX=(node_x2-node_x1)/(node_y2-node_y1)*(yBound2-node_y1)+node_x1;        if indexX<=xBound2&&indexX>=xBound1            voronoiX(1,indexI)=indexX;            voronoiY(1,indexI)=yBound2;        end    end    if node_y2>yBound2        indexX=(node_x1-node_x2)/(node_y1-node_y2)*(yBound2-node_y2)+node_x2;        if indexX<=xBound2&&indexX>=xBound1            voronoiX(2,indexI)=indexX;            voronoiY(2,indexI)=yBound2;        end    endend%%%去除由于越界边均在界外计算出的孤立点index=abs(voronoiX(1,:)-voronoiX(2,:))<0.00001&abs(voronoiY(1,:)-voronoiY(2,:))<0.00001;sum(index)voronoiX(:,index)=[];voronoiY(:,index)=[];%% %循环将边界点连接起来，作为新边index=voronoiX(1,:)==xBound1|voronoiX(1,:)==xBound2;mUnionX=[voronoiX(1,index)];mUnionY=[voronoiY(1,index)];index=voronoiX(2,:)==xBound1|voronoiX(2,:)==xBound2;mUnionX=[mUnionX voronoiX(2,index)];mUnionY=[mUnionY voronoiY(2,index)];index=voronoiY(1,:)==yBound1|voronoiY(1,:)==yBound2;mUnionX=[mUnionX voronoiX(1,index)];mUnionY=[mUnionY voronoiY(1,index)];index=voronoiY(2,:)==yBound1|voronoiY(2,:)==yBound2;mUnionX=[mUnionX voronoiX(2,index)];mUnionY=[mUnionY voronoiY(2,index)];%在找到的边界点序列中加入限定区域的四个顶点mUnionX=[mUnionX xBound1 xBound2 xBound2 xBound1];mUnionY=[mUnionY yBound2 yBound2 yBound1 yBound1];%对边界点进行排序,每条边分开,顺时针方向排序index=mUnionY(:)==yBound2;indexX=mUnionX(index);indexY=mUnionY(index);indexP=sortrows([indexX' indexY'],1);index=mUnionX(:)==xBound2;indexX=mUnionX(index);indexY=mUnionY(index);indexP=[indexP;sortrows([indexX' indexY'],-2)];index=mUnionY(:)==yBound1;indexX=mUnionX(index);indexY=mUnionY(index);indexP=[indexP;sortrows([indexX' indexY'],-1)];index=mUnionX(:)==xBound1;indexX=mUnionX(index);indexY=mUnionY(index);indexP=[indexP;sortrows([indexX' indexY'],2)];%循环去除重复点[a,b]=size(indexP);pointLeft=a;k=0;for indexI=1:a    if indexI>=pointLeft        break;    end    k=0;    for indexI2=indexI+1:pointLeft        if indexP(indexI,1)==indexP(indexI2-k,1) && indexP(indexI,2)==indexP(indexI2-k,2)            indexP(indexI2-k,:)=[];            k=k+1;            pointLeft=pointLeft-1;        end    endend%%%将新增边加入到泰森多边形矩阵中[a b]=size(indexP);for indexI=1:a-1    voronoiX=[voronoiX indexP(indexI:indexI+1,1)];    voronoiY=[voronoiY indexP(indexI:indexI+1,2)];endvoronoiX=[voronoiX [indexP(a,1);indexP(1,1)]];voronoiY=[voronoiY [indexP(a,2);indexP(1,2)]];plot(voronoiX,voronoiY,'b');

笔者使用的方法比较无脑简单，下图为实验结果，如果有更好的操作请分享出来帮助笔者进步。