- 积分
- 0
- 注册时间
- 2009-8-23
- 仿真币
-
- 最后登录
- 1970-1-1
|
本帖最后由 messenger 于 2009-8-25 17:47 编辑
在坛子里没找到答案,只好自己写了,花了我宝贵的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
%%%%%%%%%%%%%%%(完)%%%%%%%%%%%%%%%%%%%%%%%%% |
本帖子中包含更多资源
您需要 登录 才可以下载或查看,没有账号?注册
×
评分
-
1
查看全部评分
-
|