离散点凸包边界的搜索函数

注：本文为greation原创，如需转载请注明出处http://hi.baidu.com/greation。

本人最近对离散数据很是着迷，离散数据的边界问题应该是一个基本问题。所谓离散点的凸包边界，就是一个可以包围所有离散点的，凸包边界一般由离散点最外围的点连接而成。一般基于delaunay找出最大张角的向量。

function linepoints_id=chaline(cs_xss,isplot_flag)
% 本过程需要提供两个参数： cs_xss和isplot_flag
% cs_xss 为水深点，格式为矩阵，每个点一行，每行的一般格式为：
% 点序号, x坐标, y坐标, z坐标, 备注信息
% isplot_flag为绘图开关，非零为绘制，0为不绘制
% greation 2008-6-15 17:43 Lianyugang
if nargin==1
isplot_flag=0;
end
x=cs_xss(:,2);
y=cs_xss(:,3);
x=x';
y=y';
tri = delaunay(x,y);
[minx,xi]=min(x);
nextp=xi;
linepoints=[xi];
p=xi;
while(length(p)~=0)
npoint=nextp;
[he,le]=find(tri==npoint); % 得到关注点的相关三角形所在位置，即所在的行he和列le
ta=tri(he,:); % 得到得到关注点的相关三角形，三个点号组成
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% 添加根据三角形边长改正的过程 %%%%%% 待续 %%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
taz=ta'; % 处理前转置，这是一个很巧妙的地方：因为MATLAB是以列为遍历方向的，这样为下面的三角形还原
taz(find(taz==npoint))=[]; % 删除冗余信息，对于以找出的三角形而言，点号npoint就是冗余
t=reshape(taz,2,prod(size(taz))/2)'; % #三角形还原：# 还原为原三角形形式，但是已为两点式，即仅有除npoint外的二点
% 遍历该节点npoint上向量查找计数仅一次的向量，一般返回两个
ip=0; % 用以存储找出的边界点序号
p=[]; % 用以存储找出的边界点
for ik=1:prod(size(t))
if length(find(t==t(ik)))==1
ip=ip+1;
p(ip)=t(ik);
end
end
% 一般p不为空
% 若在linepoints中找到相同的p值则销毁之
if not(isempty(find(linepoints==p(1)))) & not(isempty(find(linepoints==p(2))))
p=[];
elseif not(isempty(find(linepoints==p(1))))
p(1)=[];
elseif not(isempty(find(linepoints==p(2))))
p(2)=[];
end
linepoints=[linepoints,p];
% 对第一个点的特殊照顾
if length(linepoints)==3
linepoints(1)=linepoints(2);
linepoints(2)=xi;
end
nextp=linepoints(end);
end
linepoints_id=linepoints;
if isplot_flag~=0
closelinex=[x(linepoints(end)),x(linepoints)];
closeliney=[y(linepoints(end)),y(linepoints)];
plot(x,y,'.');
hold on;
plot(closelinex,closeliney,'-r','LineWidth',3);
%triplot(tri,x,y);
% 标注点号
for k=1:size(cs_xss,1)
strtxt=num2str(cs_xss(k,1));
text(x(k),y(k),strtxt);
end
end
return
%endfunction

------------------------------------------------------------------------------------------------------------------

测试数据：

>>rss=[1 34.2 62.8 4.6
2 57 74.8 9.6
3 14.62 81.35 3
4 62.86 51.09 9.9
5 43.7 70.42 4.8
6 17.82 66.39 1.4
7 43.7 46.39 7.1
8 75.46 56.98 3.2
9 51.43 60.17 7.8
10 78.15 38.99 6.7
11 67.4 36.98 7.8
12 92.77 29.24 4.5
13 79.66 27.56 6
14 31.43 80.34 6.1
15 35.13 89.92 5
16 1.18 88.74 7.5
17 20.34 88.57 3.9
18 25.71 54.45 2.4
19 25.71 64.2 3.2
20 9.41 71.26 4.5
21 47.9 56.98 5.1
22 35.63 67.9 1.3
23 20.17 74.79 6.4
24 41.01 77.48 4.5
25 59.66 64.03 7.2
26 73.78 51.09 6.9
27 59.66 42.86 9.6
28 52.94 34.96 2.9
29 69.08 26.72 8.9
30 90.08 40.34 1.5
31 85.21 20.34 9.9
32 51.77 51.26 5.7
33 19.16 94.79 4.4
34 4.03 81.51 4.6
35 44.2 82.02 6.2
36 65.38 63.53 5.4
37 58.49 54.96 2.4
38 87.23 33.11 2.5
39 70.76 22.86 2
40 59.5 30.08 4.1
41 37.65 47.06 3
42 37.31 56.64 9.4]

>>chaline(rss) % 仅返回边界点点号

>>chaline(rss,1) % 同时绘制示意图