【原创】外圈voronoi多边形面积计算的代码

wanyou9

在坛子里没找到答案，只好自己写了，花了我宝贵的1天时间完成。由于刚接触matlab一个月，代码中不规范的地方，还请各位多指教！附件包含四个文件：两个坐标文本文件，一个m代码文件，一个fig成图文件。
%%————————————————————————————————
% 代码名称：外圈voronoi多边形面积计算的代码
% 作者：蜗牛wan
% 单位：中科院地理所
% 时间：20090825
% 基本思路：
% 1、遍历voronoi多边形顶点集合c，若集合包含编号1的顶点（暂定以为外圈voronoi多边形），
% 则根据1前后的顶点从vx、vy数组中找到与之对应的无穷远端点
% （判断标准是无穷远端点不会有其他连接点，只在vx或vy数组中出现一次），
% 注意区分1前面front和1后面behind两类无穷远端点；
% 2、得到的无穷远端点可能会有多个，但每个voronoi多边形中只会保留两个
% （判断标准是与原voronoi多边形顶点进行组合，后查看多边形内点个数是否为1，
% 且点的编号i是否与多边形顶点集合c编号一致），
% 结果可能是2个front点，可能是两个behind点，更通常的情况是1个fron+1个behind点；
% 3、处理多边形的端点顺序（非必须）；
% 4、用叠加操作，计算每个外圈voronoi多边形与边界形状叠加后的面积。
%%——————————————————————————————
clear;
clf;
%%%%%%%%%%%边界线：可定以为任意封闭曲线%%%%%%%%%
boundary=[-60,-40
    20,-40
    20, 60
    -60,60
    -60,-40];
maxX=max(boundary(:,1));
minX=min(boundary(:,1));
maxY=max(boundary(:,2));
minY=min(boundary(:,2));
%%%%%%%%%%%%%%%%%载入坐标数据%%%%%%%%%%%
fid=fopen('abc.txt','r');
geo=[];
while 1
    tline=fgetl(fid);
    if ~ischar(tline),break;end
    tline=str2num(tline);
    geo=[geo;tline];
end
fclose(fid);   
%%%%%%%%%%%%%%%%%%%%%%%%%分离case和control数据%%%%%%%%%%%
fid=fopen('abc.cas','r');
casOrCtl=[];
while 1
    tline=fgetl(fid);
    if ~ischar(tline),break;end
    tline=str2num(tline);
    casOrCtl=[casOrCtl;tline];
end
fclose(fid);   
cas=[];bg=[];
for i=1:length(casOrCtl)
    if casOrCtl(i,2)==1
        cas=[cas;geo(i,:)];
    else
        bg=[bg;geo(i,:)];
    end
end
%bg和cas数组2、3列是坐标数据
plot(bg(:,2),bg(:,3),'g+');
hold on
plot(cas(:,2),cas(:,3),'bo');
cas2=[];
for i=1:caseNum
    cas2=[cas2;cas(i,2),cas(i,3)];
end
%用两种方法创建case数据的voronoi多边形
%v存放了有界voronoi多边形的顶点(带编号)，c是每个voronoi多边形顶点编号的集合，
%c集合中编号1代表此顶点是无界的
%vx,vy中分别存放voronoi多边形的线段，对于无界的顶点也给出一个远端值
[v,c]=voronoin(cas2);
[vx,vy] = voronoi(cas(:,2),cas(:,3));
vxOrig=vx; vyOrig=vy;
%plot(vx,vy,'k-');%可用来验证多边形的最终形状
%%%%%%%%%%%%%%%%%%%%%%%%核心部分代码%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for i=1:caseNum
     %1、将多边形顶点中的NOF点替换为[vx,vy]中对应线段,注意每个多边形的nof点都包含以前以后两个点
    if any(c{i}==1)      
        nof=find(c{i}==1);
        front=nof-1;        %NOF点的前一个点
        behind=nof+1;       %NOF点的后一个点
        if front<1
            front=length(c{i});
        end
        if behind>length(c{i})  
            behind=1;
        end
        frontPnt=c{i}(front);   %frontPnt是与nof点组成直线的另一端点的编号
        behindPnt=c{i}(behind); %behindPnt也是与nof点组成直线的另一端点的编号
        vLen=length(v);       %在加入无穷远端点之前的v长度，用来最后做筛查用
        tempVF=[];tempVB=[];    %存储新添加的顶点编号，分别对应frontPnt和behindPnt
        for ii=1:length(vx)     
            %刚接触matlab不长时间，不知道数字截断函数是什么，只好自己造了一个
            if round(10000*(vx(1,ii)-v(frontPnt,1)))/10000==0 && round(10000*(vy(1,ii)-v(frontPnt,2)))/10000==0
            
%找出与frontPnt相连的nof点，加入frontPnt的后面
                if vx(2,ii)>maxX || vy(2,ii)>maxY ||vx(2,ii)<minX|| vy(2,ii)<minY
                    %v2只出现一次才能证明其是个无穷远点
                    if length(find(vxOrig==vx(2,ii)))==1 && length(find(vyOrig==vy(2,ii)))==1
                        v(length(v)+1,1)=vx(2,ii);
                        v(length(v),2)=vy(2,ii);
                        if nof~=1
                           m=[c{i}(1:front),length(v),c{i}(nof:end)];
                        else
                            m=[c{i}(1:length(c{i})),length(v)];
                        end
                        c{i}=m;
                        tempVF=[tempVF;length(v)];
                    end               
                end
            elseif round(10000*(vx(1,ii)-v(behindPnt,1)))/10000==0 && round(10000*(vy(1,ii)-v(behindPnt,2)))/10000==0
                %找出与behindPnt相连的nof点,加入behindPnt的前面
                if vx(2,ii)>maxX || vy(2,ii)>maxY ||vx(2,ii)<minX|| vy(2,ii)<minY
                    %v2只出现一次才能证明其是个无穷远点
                    if length(find(vxOrig==vx(2,ii)))==1 && length(find(vyOrig==vy(2,ii)))==1   
                       v(length(v)+1,1)=vx(2,ii);
                       v(length(v),2)=vy(2,ii);
                       if behind~=1
                           m=[c{i}(1:nof),length(v),c{i}(behind:end)];
                       else
                           m=[c{i}(1:length(c{i})),length(v)];
                       end
                        c{i}=m;                        
                        tempVB=[tempVB;length(v)];
                   end               
                end
            end
        end
        nof=find(c{i}==1); c{i}(nof)=[];
        %获得的无穷远点可能有多余的，需要根据点在多边形内的个数排除掉不正确的无穷远点      
        %构造临时的多边形，使内部只含一个cas点且编号为i的那个无穷远端点是voronoi多边形的顶点
        %目前还没发处理多个frontPnt对应无穷远端点+多个behindPnt对应无穷远端点的情况，仅考虑1对多的情况
        m=c{i}; found=0;
        if length(tempVF)>1 && length(tempVB)==1
           for j=1:length(tempVF)
             if found~=1
               nof=find(m>vLen & m~=tempVF(j) & m~=tempVB(1));
             for jj=length(nof):-1:1
                 m(nof(jj))=[];
             end         
               m=[m,m(1)];
               in=inpolygon(cas(:,2),cas(:,3),v(m,1),v(m,2));
               m(end)=[];
               if sum(in)==1
                   nof=find(in==1);
                   if nof==i
                       c{i}=m; found=1;
                   end
               else
                   m=c{i};
               end
             end
           end
        elseif length(tempVF)==1 && length(tempVB)>1
           for j=1:length(tempVB)
             if found~=1
               nof=find(m>vLen & m~=tempVF(j) & m~=tempVB(1));
               for jj=length(nof):-1:1
                   m(nof(jj))=[];
               end
               m=[m,m(1)];
               in=inpolygon(cas(:,2),cas(:,3),v(m,1),v(m,2));
               m(end)=[];
               if sum(in)==1
                   nof=find(in==1);
                   if nof==i
                       c{i}=m; found=1;
                   end
               else
                   m=c{i};
               end
             end
           end
        end
%         %处理多边形的端点顺序，需要来回比对，太繁琐，没有成功
%         if ~ispolycw(v(c{i},1),v(c{i},2))
%             [x2,y2] = poly2cw(v(c{i},1),v(c{i},2));
%             c{i}=[x2(:,1),y2(:,1)];
%         end
    end
    %voronoi多边形与区域边界相交的面积
    m=[c{i},c{i}(1)];    %将多边形点首尾相连，以便做图和求面积
    plot(v(m,1),v(m,2),'r-');   %画出带外圈的voronoi多边形图
    [ai,bi]=polybool('intersection',boundary(:,1),boundary(:,2),v(m,1),v(m,2));
    c{i,3}=polyarea(ai,bi);      
end
%%%%%%%%%%%%%%%（完）%%%%%%%%%%%%%%%%%%%%%%%%%

shunfly

1. 
   
2. for i=1:caseNum
   
3. 
   
4.     cas2=[cas2;cas(i,2),cas(i,3)];
   
5. 
   
6. end
   

复制代码

Undefined  variable caseNum

在坛子里没找到答案，只好自己写了，花了我宝贵的1天时间完成。由于刚接触matlab一个月，代码中不规范的地方，还请各位多指教！附件包含四个文件：两个坐标文本文件，一个m代码文件，一个fig成图文件。
%%————— ...
wanyou9 发表于 2009-8-25 15:59

stellazhao

高手来的，很好很强大，看下合不合用，先谢过

hongfeng101

学习了，支持一下，不愧是高手

hehevswalle

找了很久， 谢谢你 1# wanyou9

wcjhit110

谢谢楼主的代码！

xiao20081016

新手受教，多谢！