matlab凸包算法

⑴ MATLAB 一堆散点如何求包络线

你可以用convexHull来找出凸包。

%%x,y弄成列向量

dt=DelaunayTri(x,y)
k=convexHull(dt)
plot(x,y,'.','markersize',10);holdon;
plot(x(k),y(k),'r');holdoff;

⑵ matlab 里面如何用红色十字标记凸包节点

如果没记错你可以是这么写试试：plot(x(s),y(s),'r+')

⑶ 已知点集的matlab 三维凸包 用公式表达出来

用 qhull 计算三维点集的凸包
在计算几何领域，qhull 是个很强大的程序，它可以计算 2 维、3 维，以及4 维以上维度点集的凸包、Delaunay 网格、Voronoi 图，并且 Matlab 和 Octave 都基于它来提供计算几何功能，Mathematica 使用它实现 Delaunay 网格构造。不过，也正是因为它过于强大，所以我在它的源代码中逡巡了好久，也没有看懂。现在我能做到的，就是找个梯子先爬上这个黑箱子。

因为 qhull 是一个复杂的命令行工具箱，用户可以通过各种命令选项去调用适当的程序。比如，要对点集进行 Delaunay 网格化，可以直接使用 "qdelaunay" 命令来实现，也可以通过 "qhull d Qbb" 命令来实现。
在 qhull 工具箱中，要构建点集的凸包，可以调用 "qconvex" 命令来实现，而且 "qhull" 如果在未设定命令选项时，默认调用的程序就是 qconvex。对于我要解决的问题，即使用 qhull 构造三维点集的凸包而言，基本命令格式如下：
$ qconvex [选项] < input_file > output_file
qconvex 程序的行为由一组功能选项来控制，常用的如下：
Qt - 三角化输出，也就是输出由三角面片组合而成的凸包数据
QJ - 对于近似于平面的数据进行一些简化，譬如对于三维点集，
- 可以保证不会出现 4 点共面的情况
Tv - 从结构、凸性以及数据包含等方面对凸包构建结果进行校验
- - 输出 qconvex 所有选项信息
对于凸包构建结果的输出，qconvex 提供了一组输出控制选项，常用的如下：
s - 输出凸包构建结果概要 (default)
i - 输出凸包上每个面片的顶点
n - 输出凸包上每个面的方程系数
p - 输出要进行凸包求解的点集的坐标
Fx - 输出极点（凸包顶点）
FA - 输出凸包的面积和体积
o - 采用 OFF 格式输出凸包构建结果（维度, 顶点数, 面数, 边数）
G - 采用 Geomview 格式输出凸包构建结果（只支持 2 维至 4 维）
m - 采用 Mathematica 格式输出凸包构建结果（只支持 2 维与 3 维）
TO 文件名 - 将凸包构建结果输出到文件
我们来玩真格的。首先准备好一份存放三维数据点信息的文本文件，文件的第一行是点数，其余每行都是一个数据点的 x, y, z 坐标信息。对于下图所示的 venus 点云，其数据文件 venus.asc 格式为：

26218
3.554705 199.173300 8.394049
3.584999 199.553600 10.168500
3.648500 200.610500 9.662390
3.667395 198.429900 10.576820
3.740701 198.771200 12.616080
3.796498 200.566700 7.518420
3.798301 197.899800 9.092709
3.814104 197.907700 12.370720
3.839600 200.038700 12.814610
... ... ... ... ... ...
现在使用 qconvex 对 venus.asc 文件所包含的点集构建凸包，采用 OFF (Object File Format) 格式输出：
$ qconvex o < venus.asc TO result.off
qconvex 输出的 off 格式文件自上而下由三部分构成：
文件头信息，即文件的第一行是数据的维度，第二行的内容是凸包顶点数、面片数和边数；
点表，存放凸包顶点的坐标信息；
面表，按逆时针次序记录面片顶点在点表中的序号（从 0 开始）。
在 off 文件的面表区域，第一列数字用来表示每个面片所含的顶点的个数。
在 linux 下，可以使用 geomview 来显示 OFF 格式文件，但前提是将 qconvex 输出的 off 文件的第一行内容替换为 "OFF"。下面是一份 geomview 所能接受的 OFF 文件格式，显示结果如下图所示。

# 文件头 （本行文本是注释，实际中应当去掉）
OFF
26218 4870 7305
# 点表 （本行文本是注释，实际中应当去掉）
3.584999 199.5536 10.1685
3.740701 198.7712 12.61608
3.796498 200.5667 7.51842
... ... ... ... ... ...
# 面表 （本行文本是注释，实际中应当去掉）
3 9864 9563 9674
3 9563 9864 9833
3 6318 1422 452

⑷ 点在球面上均匀分布 在matlab中求解

今天回去想了一下，有另一种办法生成球面上均匀分布的点

试了一下，速度比之前遗传算法的快，算500个点用了两分钟

也有可能是之间的算法中计算凸包比较耗时间

新的办法是基于一个物理模型，让球面上的N个点互相之间有排斥力

排斥力随点间的距离增加而衰减，例如可以是距离平方衰减的力

可以将问题比拟为球面上N个带相同等量电荷的小点的运动问题

当运动最后平衡，小球停止运动的时候，由于斥力的存在，小球的分布就会是均匀的

这时候问题转变为一个多变量微分方程的数值问题，最终的平衡状态就是我们想要的

用最简单的欧拉差分公式，为每个点设定x，y，z，vx，vy，vz六个状态量，描述运动状态

每一步计算跟新所有点的运动状态，而初始状态随机分配，速度都为0

function[]=main()

N=12;%点数量

a=rand(N,1)*2*pi;%根据随机求面均匀分布，先生成一个初始状态

b=asin(rand(N,1)*2-1);

r0=[cos(a).*cos(b),sin(a).*cos(b),sin(b)];

v0=zeros(size(r0));

G=1e-2;%斥力常数，试验这个值比较不错

forii=1:200%模拟200步，一般已经收敛，其实可以在之下退出

[rn,vn]=countnext(r0,v0,G);%更新状态

r0=rn;v0=vn;

end

plot3(rn(:,1),rn(:,2),rn(:,3),'.');holdon;%画结果

[xx,yy,zz]=sphere(50);

h2=surf(xx,yy,zz);%画一个单位球做参考

set(h2,'edgecolor','none','facecolor','r','facealpha',0.7);

axisequal;

axis([-11-11-11]);

holdoff;

end

function[rnvn]=countnext(r,v,G)%更新状态的函数

%r存放每点的x，y，z数据，v存放每点的速度数据

num=size(r,1);

dd=zeros(3,num,num);%各点间的矢量差

form=1:num-1

forn=m+1:num

dd(:,m,n)=(r(m,:)-r(n,:))';

dd(:,n,m)=-dd(:,m,n);

end

end

L=sqrt(sum(dd.^2,1));%各点间的距离

L(L<1e-2)=1e-2;%距离过小限定

F=sum(dd./repmat(L.^3,[311]),3)';%计算合力

Fr=r.*repmat(dot(F,r,2),[13]);%计算合力径向分量

Fv=F-Fr;%切向分量

rn=r+v;%更新坐标

rn=rn./repmat(sqrt(sum(rn.^2,2)),[13]);

vn=v+G*Fv;%跟新速度

end

这种办法还有个好处，可以看到每步的结果，

以下是用算法某次求12点均匀分布的变化过程

最后得到的20面体很接近正20面体了

⑸ 求二维凸包的增量算法，最好能详细解释一下

#include <iostream>
#include <algorithm>
#include <cmath>
using namespace std;
typedef double Type; // 注意下面的fabs().
const int maxn = 1005;
const double EPS = 1e-8;
const double Pi = acos(-1.0);
struct Point
{
Type x, y;
Point () {}
Point (Type & xx, Type & yy) : x(xx), y(yy) {}
};
// 判断正负.
int dblcmp(Type d)
{
if (fabs(d) < EPS) return 0; // 注意数据类型不同，abs()也不同.
return d > 0 ? 1 : -1;
}
// 叉乘.
// cross proct of (c->a) and (c->b).
Type Cross(Point & c, Point a, Point b)
{
return (a.x - c.x) * (b.y - c.y) - (b.x - c.x) * (a.y - c.y);
}
Type Distance(Point & u, Point & v)
{
return sqrt( 0.0 + (u.x - v.x) * (u.x - v.x) + (u.y - v.y) * (u.y - v.y) );
}
int n;
Point point[maxn];
int Stack[maxn];
int top;
double ans;
bool cmp(const Point & a, const Point & b)
{
if (a.y == b.y) return a.x < b.x;
return a.y < b.y;
}
void graham_scan(void)
{
int i;
int temp_top;
if (n <= 1)
{
top = 0;
return ;
}
sort(point, point + n, cmp); // point[0]即为起点.
// 做右链.
top = -1;
Stack[++top] = 0; Stack[++top] = 1;
for (i = 2; i < n; i++)
{
while (top >= 1 && dblcmp(Cross(point[Stack[top - 1]], point[i], point[Stack[top]])) >= 0) top--; // 如果不能左转，则退栈. 如果只要求极点，则共线的点也是不要的(即要加等于).
Stack[++top] = i;
}
temp_top = top; // 此时的栈顶元素一定是第n个点.
// 做左链.
Stack[++top] = n - 2;
for (i = n - 3; i >= 0; i--)
{
while (top >= temp_top + 1 && dblcmp(Cross(point[Stack[top - 1]], point[i], point[Stack[top]])) >= 0) top--; // 如果不能左转，则退栈. 如果只要求极点，则共线的点也是不要的(即要加等于).
Stack[++top] = i;
}
// 此时的栈顶元素是第1个点.(如果凸包是一条直线，则左右链倒置相同.)
// 凸包的顶点为point[Stack[0]] 到 point[Stack[top - 1]].
}
int main(void)
{
int i;
scanf("%d", &n);
for (i = 0; i < n; i++)
{
scanf("%lf %lf", &point[i].x, &point[i].y);
}
graham_scan();
ans = 0;
for (i = 0; i < top; i++)
{
printf("%lf %lf\n", point[Stack[i]].x, point[Stack[i]].y);
ans += Distance(point[Stack[i]], point[Stack[i + 1]]); // point[Stack[top]] = point[Stack[0]].
}
printf("%lf\n", ans);
return 0;
}
/******************************************************************************************************************************************************/
/******************************************************************************************************************************************************/
/******************************************************************************************************************************************************/
/******************************************************************************************************************************************************/
/******************************************************************************************************************************************************/
/******************************************************************************************************************************************************/
/******************************************************************************************************************************************************/
/******************************************************************************************************************************************************/
/* 按照lrj的黑书来写的.
适用条件：简单多边形(点按顺时针或逆时针给出).
复杂度：O(n).
*/
#include <iostream>
#include <algorithm>
#include <cmath>
using namespace std;
typedef double Type; // 注意下面的fabs().
const int maxn = 1005;
const double EPS = 1e-8;
const double Pi = acos(-1.0);
struct Point
{
Type x, y;
Point () {}
Point (Type & xx, Type & yy) : x(xx), y(yy) {}
};
// 判断正负.
int dblcmp(Type d)
{
if (fabs(d) < EPS) return 0; // 注意数据类型不同，abs()也不同.
return d > 0 ? 1 : -1;
}
// 叉乘.
// cross proct of (c->a) and (c->b).
Type Cross(Point & c, Point a, Point b)
{
return (a.x - c.x) * (b.y - c.y) - (b.x - c.x) * (a.y - c.y);
}
Type Distance(Point & u, Point & v)
{
return sqrt( 0.0 + (u.x - v.x) * (u.x - v.x) + (u.y - v.y) * (u.y - v.y) );
}
int n;
Point point[maxn];
int Stack[2 * maxn]; // 两头栈.
int bot, top; // 栈底，栈顶.
double ans;
void Melkman(void)
{
int i;
int temp;
Stack[n] = 0;
// 注意：前三个点不能是共线的.
for (i = 1; i < n; i++)
{
Stack[n + 1] = i; // 当三点平行时要的是后一个点.
if (dblcmp(Cross(point[Stack[n]], point[Stack[n + 1]], point[i + 1]))) break; // 前三个点不共线.
}
bot = n, top = n + 1;
Stack[++top] = Stack[--bot] = i + 1;
// 保证开始的三个点成右手系，否则交换Stack[n]和Stack[n + 1] .
if (dblcmp(Cross(point[Stack[n]], point[Stack[n + 1]], point[Stack[n + 2]])) < 0)
{
temp = Stack[n]; Stack[n] = Stack[n + 1]; Stack[n + 1] = temp;
}
// 维护栈里的点为右手系(即栈中任意连续三点组成的路径是左旋的，或成直线).
for (i = i + 2; i < n; i++)
{
if (dblcmp(Cross(point[Stack[top - 1]], point[Stack[top]], point[i])) > 0 &&
dblcmp(Cross(point[Stack[bot]], point[Stack[bot + 1]], point[i])) > 0)
{ // 如果该点对于栈顶左旋且对于栈底右旋，则说明该点在凸包内.
continue;
}
while (dblcmp(Cross(point[Stack[top - 1]], point[Stack[top]], point[i])) <= 0) top--;
Stack[++top] = i;
while (dblcmp(Cross(point[Stack[bot]], point[Stack[bot + 1]], point[i])) <= 0) bot++;
Stack[--bot] = i;
}
}
int main(void)
{
int i;
scanf("%d", &n);
for (i = 0; i < n; i++)
{
scanf("%lf %lf", &point[i].x, &point[i].y);
}
Melkman(); // 得到的凸包点序列也是按极角序的.
cout << top - bot << endl;
for (i = bot; i < top; i++)
{
printf("%lf %lf\n", point[Stack[i]].x, point[Stack[i]].y);
ans += Distance(point[Stack[i]], point[Stack[i + 1]]); // Stack[top]为第1个点.
}
printf("%lf\n", ans);
return 0;
}

⑹ matlab如何求大量数据中分布最多的范围及中心点。

楼主是不是想求出一个最小半径的圆，圆内包含所有的点？这个问题很有趣。

寻找这个圆的时候注意一下几点：
1.这个圆必然穿过图中某些靠外围的点，这样才是最小半径的圆。
2.几何中我们知道，三个点可以确定一个圆, 我们就是需要找出这三个点来.

算法如下：1.先求这些点对应的凸包，已经有现成的算法。
2.生成凸包后，在看凸包上哪三点确定的圆可以包含凸包。

当然如果楼主讨论的不是以上所述，而是模式分类的话，建议看看数据分类方法。可以搜索关键字：Gaussian mixtrual model, expectation-maximization algorithm 和 k-mean algorithm 学习下相关的知识。