我如何使用MATLAB的'isosurface'函数创建三角形球体 [英] How can I create a triangulated sphere using 'isosurface' function of MATLAB
本文介绍了我如何使用MATLAB的'isosurface'函数创建三角形球体的处理方法，对大家解决问题具有一定的参考价值，需要的朋友们下面随着小编来一起学习吧！
问题描述
我想要像这样的东西，
http://imageshack.us/a/img198/5041/ 71183923.png
我搜索了一下，看到我可以使用MATLAB的函数 isosurface ，但三角形没有相同的面积和形式，它们基本上是用方形对角线分成两半的正方形。
这是我的代码：
>> [X，Y，Z] =网格（-10：10，-10：10，-10：10）;
>> V = sqrt（X. ^ 2 + Y. ^ 2 + Z. ^ 2）;
>> fv =等值面（X，Y，Z，V）
>> p = patch（fv）;
>> set（p，'EdgeColor'，'g'）
，结果如下：
http://imageshack.us/a /img818/7125/72103297.png
解决方案使用此功能从文件交换
函数[p，t] = TriSphere（N，R）
％TRISPHERE：使用
返回球体的三角化模型％icosaedron细分方法。
％
％INPUT：
％N（整数）表示细分的数量，
％可以使用0-inf之间的值。 N越大越好，表面看上去
％，但表面结构花费的时间会更多，而三角形的
％将会放在输出模型中。
％
％OUTPUT：
％在p（nx3）和t（mx3）中，我们可以找到点和三角形索引
％的模型。该球体应该是单位半径，并且以
％（0,0,0）为中心。要获得以不同位置为中心的球体，或以
％不同半径的球体，只需转换和缩放
％转化。这些操作并没有被这个代码执行，因为它是
％非常方便，以时间性能的顺序，在函数中执行这个操作
％，避免每次调用构造步骤。
％
％注意：
％这个函数比matlab命令球体更加高效，以
％的模型尺寸/重建精度为准。这是由于
％的良好traingulated模型，需要少量的补丁
％相同的几何recostruction准确性。在时间执行和模型细分的灵活性方面可能的改进是
％。
％
％例如：
％
％N = 5;
％
％[p，t] = TriSphere（N）; （t，p（：，1），p（：，2），p（：，3））;
％
％ axis vis3d
％view（3）
％作者：Giaccari Luigi创建日期：2009年25月04日％
％信息/错误/问题/建议：giaccariluigi @ msn。 （速度/可读性）
％原始名称：BUILDSPHERE
％
％由Rody Oldenhuis调整（速度/可读性）
误差陷阱
错误（nargchk（1， 1，nargin））;
错误（nargoutchk（1,2，nargout））;
if〜isscalar（N）
错误（'Buildsphere：N_mustbe_scalar'，...
'输入N必须是标量'。
end
如果round（N）〜= N
错误（'Buildsphere：N_mustbe_scalar'，...
'输入N必须是一个整数值'。
end
％从Jon Leech的文件中获取坐标
％单位球面上的十二面体顶点
tau = 0.8506508083520400; ％t =（1 + sqrt（5））/ 2，tau = t / sqrt（1 + t ^ 2）
one = 0.5257311121191336; ％one = 1 / sqrt（1 + t ^ 2）（unit sphere）
p = [
+ tau + one +0％ZA
-tau + one +0％ZB
-tau -one + 0％ZC
+ tau -one + 0％ZD
+ one +0 + tau％YA
+ one +0 -tau％YB
- 1 + 0 -tau％YC
-one +0 + tau％YD
+0 + tau + one％XA
+0 -tau + one％XB
+0 - tau -one％XC
+0 + tau -one]; ％XD
％单位二十面体的结构
t = [
5 8 9
5 10 8
6 12 7
6 7 11
1 4 5
1 6 4
3 2 8
3 7 2
9 12 1
9 2 12
10 4 11
10 11 3
9 1 5
12 6 1
5 4 10
6 11 4
8 2 9
7 12 2
8 10 3
7 3 11];
％可能的快速退出
如果N == 0，返回，结束
％加载预生成的trispheres（现在最多8个...）
如果N <= 8
S = load（['TriSphere'，num2str（N），'.mat']，'pts'，'idx'）;
p = S.pts; t = S.idx;
如果nargin == 2，p = p * R;如果还有更多的请求（为什么你会在地球上？），请确保从最大预加载三球$ b $开始
％返回
else
％ b S = load（'TriSphere8.mat'，'pts'，'idx'）;
p = S.pts; t = S.idx;清除S; N0 = 10;
end
％我们有多少个三角形/顶点？
nt = size（t，1）; np = size（p，1）; totp = np;
％计算ii = N0时的最终分数
：N，totp = 4 * totp - 6;结束
％初始化点阵列
p = [p;零（totp-12，3）];
％确定三角测量指标的适当类别
numbits = 2 ^ ceil（log（log（totp + 1）/ log（2））/ log（2））;
castToInt = ['uint'，num2str（numbits）];
％需要时发出警告
如果numbits> 64
警告（'TriSphere：too_many_notes'，...
['给定的迭代次数需要％s类准确'，...
'表示三角测量指标。双反而;期待'，...
'奇怪的结果！']）;
castToInt = @double;
else
castToInt = str2func（castToInt）;
结束
％细化二十面体N次
用于ii = N0：N
％初始化内循环
tell = t;
t =零（nt * 4，3）;
％使用稀疏。是的，它的循环速度较慢，但对于N = 6，尺寸为
％的尺寸已经达到了10,000x10,000，增长了4倍，每增加
％增加N;它的内存密集程度太高以致于无法使用零（）。
peMap = sparse（np，np）;
ct = 1;
％循环遍历所有旧三角形
for j = 1：nt
％一些辅助变量
p1 = tell（j，1）;
p2 = tell（j，2）;
p3 =被告知（j，3）;
x1 = p（p1,1）; x2 = p（p2,1）; x3 = p（p3,1）;
y1 = p（p1,2）; y2 = p（p2,2）; y3 = p（p3,2）;
z1 = p（p1,3）; z2 = p（p2,3）; z3 = p（p3,3）;
％first edge
％ - = - = - = - = = = - = - = - = - = - = - = - = - = - =
％保留三角形方向
if p1 < p2，p1m = p1; p2m = p2;否则p2m = p1; p1m = p2;结束
％如果该点不存在，则计算新点
p4 = peMap（p1m，p2m）;
if p4 == 0
np = np + 1;
p4 = np;
peMap（p1m，p2m）= np;％＃ok
p（np，1）=（x1 + x2）/ 2;
p（np，2）=（y1 + y2）/ 2;
p（np，3）=（z1 + z2）/ 2;
end
％second edge
％ - = - = - = - = - = - = - = - = - = - = - = - = - - = - = - = - = - = - =
％保留三角形方向
if p2 < P3; p2m = p2; p3m = p3;否则p2m = p3; p3m = p2;结束
％如果该点不存在，则计算新点
p5 = peMap（p2m，p3m）;
如果p5 == 0
np = np + 1;
p5 = np;
peMap（p2m，p3m）= np;％＃ok
p（np，1）=（x2 + x3）/ 2;
p（np，2）=（y2 + y3）/ 2;
p（np，3）=（z2 + z3）/ 2;
end
％第三边缘
％ - = - = - = - = - = - = - = - = - - = - = - = - = - = - = =
％保留三角形方向
if p1 < P3; p1m = p1; p3m = p3;否则p3m = p1; p1m = p3;结束
％如果该点不存在，则计算新点
p6 = peMap（p1m，p3m）;
如果p6 == 0
np = np + 1;
p6 = np;
peMap（p1m，p3m）= np;％＃ok
p（np，1）=（x1 + x3）/ 2;
p（np，2）=（y1 + y3）/ 2;
p（np，3）=（z1 + z3）/ 2;
结束
％分配新三角形
％细化索引
％p1
％/ \
％/ t1\
％p6 / ____ \p4
％/ \ / \
％/ t4\t2 / t3\
％/ ____ \ / ____ \
％p3 p5 p2
t（ct，1）= p1; t（ct，2）= p4; t（ct，3）= p6; ct = ct + 1;
t（ct，1）= p4; t（ct，2）= p5; t（ct，3）= p6; ct = ct + 1;
t（ct，1）= p4; t（ct，2）= p2; t（ct，3）= p5; ct = ct + 1;
t（ct，1）= p6; t（ct，2）= p5; t（ct，3）= p3; ct = ct + 1;
end％end subloop
％更新三角形数量
nt = ct-1;
end％end main loop
$ b $％将所有点归一化为1（或R）
p = bsxfun（@rdivide，p，sqrt（sum（p。^ 2,2） ））;
if（nargin == 2），p = p * R;结束
％将t转换为适当的整数类
t = castToInt（t）;
end％funciton TriSphere
How can I create a Triangulated sphere with faces of triangles with the same area each one.
I want something like this,
http://imageshack.us/a/img198/5041/71183923.png
and I searched and saw I could use the function isosurfaceof MATLAB, but the triangles don't have equal area and form and they are essentially squares divided in two with the square diagonal.
here's my code:
>> [X,Y,Z] = meshgrid(-10:10,-10:10,-10:10);
>> V = sqrt(X.^2+Y.^2+Z.^2);
>> fv = isosurface(X,Y,Z,V)
>> p = patch(fv);
>> set(p,'EdgeColor','g')
and the result is below:
http://imageshack.us/a/img818/7125/72103297.png
解决方案
Use this function from file exchange
function [p, t] = TriSphere(N, R)
% TRISPHERE: Returns the triangulated model of a sphere using the
% icosaedron subdivision method.
%
% INPUT:
% N (integer number) indicates the number of subdivisions,
% it can assumes values between 0-inf. The greater N the better will look
% the surface but the more time will be spent in surface costruction and
% more triangles will be put in the output model.
%
% OUTPUT:
% In p (nx3) and t(mx3) we can find points and triangles indexes
% of the model. The sphere is supposed to be of unit radius and centered in
% (0,0,0). To obtain spheres centered in different location, or with
% different radius, is just necessary a traslation and a scaling
% trasformation. These operation are not perfomed by this code beacuse it is
% extrimely convinient, in order of time perfomances, to do this operation
% out of the function avoiding to call the costruction step each time.
%
% NOTE:
% This function is more efficient than the matlab command sphere in
% terms of dimension fo the model/ accuracy of recostruction. This due to
% well traingulated model that requires a minor number of patches for the
% same geometrical recostruction accuracy. Possible improvement are possible
% in time execution and model subdivision flexibilty.
%
% EXAMPLE:
%
% N=5;
%
% [p,t] = TriSphere(N);
%
% figure(1) axis equal hold on trisurf(t,p(:,1),p(:,2),p(:,3)); axis vis3d
% view(3)
% Author: Giaccari Luigi Created:25/04/2009%
% For info/bugs/questions/suggestions : giaccariluigi@msn.com
% ORIGINAL NAME: BUILDSPHERE
%
% Adjusted by Rody Oldenhuis (speed/readability)
% error traps
error(nargchk(1,1,nargin));
error(nargoutchk(1,2,nargout));
if ~isscalar(N)
error('Buildsphere:N_mustbe_scalar',...
'Input N must be a scalar.');
end
if round(N) ~= N
error('Buildsphere:N_mustbe_scalar',...
'Input N must be an integer value.');
end
% Coordinates have been taken from Jon Leech' files
% Twelve vertices of icosahedron on unit sphere
tau = 0.8506508083520400; % t = (1+sqrt(5))/2, tau = t/sqrt(1+t^2)
one = 0.5257311121191336; % one = 1/sqrt(1+t^2) (unit sphere)
p = [
+tau +one +0 % ZA
-tau +one +0 % ZB
-tau -one +0 % ZC
+tau -one +0 % ZD
+one +0 +tau % YA
+one +0 -tau % YB
-one +0 -tau % YC
-one +0 +tau % YD
+0 +tau +one % XA
+0 -tau +one % XB
+0 -tau -one % XC
+0 +tau -one]; % XD
% Structure for unit icosahedron
t = [
5 8 9
5 10 8
6 12 7
6 7 11
1 4 5
1 6 4
3 2 8
3 7 2
9 12 1
9 2 12
10 4 11
10 11 3
9 1 5
12 6 1
5 4 10
6 11 4
8 2 9
7 12 2
8 10 3
7 3 11 ];
% possible quick exit
if N == 0, return, end
% load pre-generated trispheres (up to 8 now...)
if N <= 8
S = load(['TriSphere', num2str(N), '.mat'],'pts','idx');
p = S.pts; t = S.idx;
if nargin == 2, p = p*R; end
return
else
% if even more is requested (why on Earth would you?!), make sure to START
% from the maximum pre-loadable trisphere
S = load('TriSphere8.mat','pts','idx');
p = S.pts; t = S.idx; clear S; N0 = 10;
end
% how many triangles/vertices do we have?
nt = size(t,1); np = size(p,1); totp = np;
% calculate the final number of points
for ii=N0:N, totp = 4*totp - 6; end
% initialize points array
p = [p; zeros(totp-12, 3)];
% determine the appropriate class for the triangulation indices
numbits = 2^ceil(log(log(totp+1)/log(2))/log(2));
castToInt = ['uint',num2str(numbits)];
% issue warning when required
if numbits > 64
warning('TriSphere:too_many_notes',...
['Given number of iterations would require a %s class to accurately ',...
'represent the triangulation indices. Using double instead; Expect ',...
'strange results!']);
castToInt = @double;
else
castToInt = str2func(castToInt);
end
% refine icosahedron N times
for ii = N0:N
% initialize inner loop
told = t;
t = zeros(nt*4, 3);
% Use sparse. Yes, its slower in a loop, but for N = 6 the size is
% already ~10,000x10,000, growing by a factor of 4 with every
% increasing N; its simply too memory intensive to use zeros().
peMap = sparse(np,np);
ct = 1;
% loop trough all old triangles
for j = 1:nt
% some helper variables
p1 = told(j,1);
p2 = told(j,2);
p3 = told(j,3);
x1 = p(p1,1); x2 = p(p2,1); x3 = p(p3,1);
y1 = p(p1,2); y2 = p(p2,2); y3 = p(p3,2);
z1 = p(p1,3); z2 = p(p2,3); z3 = p(p3,3);
% first edge
% -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
% preserve triangle orientation
if p1 < p2, p1m = p1; p2m = p2; else p2m = p1; p1m = p2; end
% If the point does not exist yet, calculate the new point
p4 = peMap(p1m,p2m);
if p4 == 0
np = np+1;
p4 = np;
peMap(p1m,p2m) = np;%#ok
p(np,1) = (x1+x2)/2;
p(np,2) = (y1+y2)/2;
p(np,3) = (z1+z2)/2;
end
% second edge
% -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
% preserve triangle orientation
if p2 < p3; p2m = p2; p3m = p3; else p2m = p3; p3m = p2; end
% If the point does not exist yet, calculate the new point
p5 = peMap(p2m,p3m);
if p5 == 0
np = np+1;
p5 = np;
peMap(p2m,p3m) = np;%#ok
p(np,1) = (x2+x3)/2;
p(np,2) = (y2+y3)/2;
p(np,3) = (z2+z3)/2;
end
% third edge
% -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
% preserve triangle orientation
if p1 < p3; p1m = p1; p3m = p3; else p3m = p1; p1m = p3; end
% If the point does not exist yet, calculate the new point
p6 = peMap(p1m,p3m);
if p6 == 0
np = np+1;
p6 = np;
peMap(p1m,p3m) = np;%#ok
p(np,1) = (x1+x3)/2;
p(np,2) = (y1+y3)/2;
p(np,3) = (z1+z3)/2;
end
% allocate new triangles
% refine indexing
% p1
% /\
% /t1\
% p6/____\p4
% /\ /\
% /t4\t2/t3\
% /____\/____\
% p3 p5 p2
t(ct,1) = p1; t(ct,2) = p4; t(ct,3) = p6; ct = ct+1;
t(ct,1) = p4; t(ct,2) = p5; t(ct,3) = p6; ct = ct+1;
t(ct,1) = p4; t(ct,2) = p2; t(ct,3) = p5; ct = ct+1;
t(ct,1) = p6; t(ct,2) = p5; t(ct,3) = p3; ct = ct+1;
end % end subloop
% update number of triangles
nt = ct-1;
end % end main loop
% normalize all points to 1 (or R)
p = bsxfun(@rdivide, p, sqrt(sum(p.^2,2)));
if (nargin == 2), p = p*R; end
% convert t to proper integer class
t = castToInt(t);
end % funciton TriSphere
这篇关于我如何使用MATLAB的'isosurface'函数创建三角形球体的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！
查看全文
![http://imageshack.us/a/img198/5041/ 71183923.png](\"http://imageshack.us/a/img198/5041/71183923.png\"){kind=link}
![http://imageshack.us/a /img818/7125/72103297.png](\"http://imageshack.us/a/img818/7125/72103297.png\"){kind=link}
