Stats

Popular Posts

Followers

Mathematica 教學:Convex Hull

戴忠淵 於 2020年8月25日星期二 下午2:09 發表

splitdata[datatemp___]:=Block[{eqn,temp,x,y,group,vertex,ans},
temp=SortBy[datatemp,Left][[{1,-1}]];
eqn[x_,y_]=y-temp[[1,2]]-
(temp[[1,2]]-temp[[2,2]])/(temp[[1,1]]-temp[[2,1]])(x-temp[[1,1]]);
group={Select[datatemp,eqn[#[[1]],#[[2]]]> 0&],
Select[datatemp,eqn[#[[1]],#[[2]]]< 0&]};

vertex={SortBy[group[[1]],Abs@eqn[#[[1]],#[[2]]]&][[-1]],
SortBy[group[[2]],Abs@eqn[#[[1]],#[[2]]]&][[-1]]};

ans={{{temp[[1]],vertex[[1]]},
Select[group[[1]],RegionMember[
Polygon[Insert[temp,vertex[[1]],-1]],#]==False&&#[[1]]< vertex[[1,1]]&]},
{{temp[[2]],vertex[[1]]},
Select[group[[1]],RegionMember[
Polygon[Insert[temp,vertex[[1]],-1]],#]==False&&#[[1]]> vertex[[1,1]]&]},
{{temp[[1]],vertex[[2]]},
Select[group[[2]],RegionMember[
Polygon[Insert[temp,vertex[[2]],-1]],#]==False&&#[[1]]< vertex[[2,1]]&]},
{{temp[[2]],vertex[[2]]},
Select[group[[2]],RegionMember[
Polygon[Insert[temp,vertex[[2]],-1]],#]==False&&#[[1]]> vertex[[2,1]]&]}
};

If[Length[#[[2]]]==0,{#[[1]],{#[[1,1]]}},{#[[1]],#[[2]]}]&/@ans
];

convexhull[{p1__,p2__},datatemp___]:=Block[{eqn,temp,x,y,group,vertex,ans},
temp=SortBy[{p1,p2},Left];
If[Length[datatemp]==0,{{temp,{}},{temp,{}}},
If[Length[datatemp]==1,{{{temp[[1]],datatemp[[1]]},{}},
{{temp[[2]],datatemp[[1]]},{}}},

eqn[x_,y_]=y-temp[[1,2]]-(temp[[1,2]]-temp[[2,2]])/
(temp[[1,1]]-temp[[2,1]])(x-temp[[1,1]]);

vertex=SortBy[datatemp,Abs@eqn[#[[1]],#[[2]]]&][[-1]];

group={Select[datatemp,#[[1]]< vertex[[1]]&],
Select[datatemp,#[[1]] > vertex[[1]]&]};

ans={{{temp[[1]],vertex},
Select[group[[1]],RegionMember[Polygon[Insert[temp,vertex,-1]],#]==
False&&#[[1]]< vertex[[1]]&]},
{{vertex,temp[[2]]},
Select[group[[2]],RegionMember[Polygon[Insert[temp,vertex,-1]],#]==
False&&#[[1]]> vertex[[1]]&]}}
]]
];

data=RandomReal[NormalDistribution[0,2],{250,2}];
Needs["ComputationalGeometry`"]

test=splitdata[data];
xx=Nest[Flatten[convexhull@@@#,1]&,test,3];
myline=Line/@Tally[xx[[All,1]]][[All,1]];
mypts=Point@Tally[xx[[All,1]]][[All,1,1]];

ListPlot[data,AspectRatio->1,
Epilog->{Red,Thickness[0.015],myline,Blue,PointSize[0.05],mypts}]
 
Tags: , ,

讀者回應 ( 0 意見 )

發佈留言

Please leave your name and tell me what you thought about this site. Comments, suggestions and views are welcomed.

如果這篇文章對你有幫助,那請留個訊息給我~