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 學習下相關的知識。