matlab求椭圆的弧长,用MATLAB实现求椭球上任意两点的最短弧长
基于法向矢量导向的求椭球上两点的最短弧长
问题分析
求椭球上任意两点间的最短弧长用数学来推算解析解的话十分复杂,因此考虑通过使用计算机来近似求解。问题的难点在于怎样让每一步都是处在最优的状态,以及怎样使每一步的方向都尽量处在该点能够选择的最优方向上。
求椭球上两点的最短路,很容易想到用传统的最短路算法如迪杰斯特拉算法或弗洛伊德算法求解,但是在椭球将步长离散化后,构造以及运用邻接矩阵十分复杂,在保证一定精度的条件下需要极大的运算量,所以最短路算法暂时没有太好的应用条件。
另一种考虑过的思路为采用经典优化算法,如粒子群算法或者蚁群算法等进行迭代求解。这种方法理论上是可行的,但是由于构造出椭球面的离散化极坐标坐标点相当多,以及这类算法本身有较大的运算量以及结果的随机性,暂时不予以考虑。
还曾考虑过用梯度下降算法来求解,就是将椭球旋转至使目标点的法向量为沿z轴负方向时构造出此时的椭球的函数表达式,即可使用梯度下降算法。这种方法是最为理想的方法,但是椭球的旋转涉及到仿射变换,构造变换矩阵非常复杂,所以暂时不予考虑。
根据梯度下降算法,可以很自然的想到既然椭球不转那我们能不能直接利用目标点的外法向量作为一个基准方向来进行导向。由次我们就想出了基于法向矢量导向的椭球上最短路求解。
模型假设
(1) 假设椭球上的网格足够小,运动点在椭球上的运动方向近似为椭球表面的切线方向
(2) 运动点到达距离终点的欧式距离小于一定范围时就视为到达终点
(3) 椭球上的运动点每次移动的欧式距离近似等于在表面移动时的弧长
(4) 椭球上的运动点的每次移动都会选择最优的方向
模型建立与算法设计
步骤
操作
1
给定起点A和终点B的椭球极坐标,计算其直角坐标,并计算B的外法向量记为BW
2
计算A点在二维极坐标中以A为中心的一定宽度方格内的点Mi的极坐标
3
计算所有Mi的直角坐标,并且计算出向量AMi,并将其全部化为单位向量
4
计算BW*AMi的值,并找出使该值最小的Mi点作为下一个A点,重复步骤1
变量与符号说明
A……………………………………………….代指椭球表面的运动点的当前位置
B………………………………………………………………………………….终点
AMi _{i}i……………………….椭球运动点在某次运动中可选择的单位方向向量集合
BW…………………………………………………椭球终点的垂直于表面的法向量
di _{i}i ………………………………………………………每次运动点移动的欧式距离
ϵ \epsilonϵ ………………………………………………………………………….终点判断值
D…………………………………………………………………………总路径的长度
模型建立与算法设计
当A要向B靠近时,求使AMi的单位方向向量与BW乘积最大的AMi作为从当前运动点指向下一个运动点的向量。对于共线的向量AMi的情况 ,取AMi长度较小的一方指向的点作为下一个A点,这样可以通过减小步长来适当降低误差。
注意注意,这里的A点周围的单位向量其实并不在一个平面内,因为我们找的是椭球上A相邻极坐标的点,所以只是椭球上从一个点出发到周围多个点的连线,我们只是用这些连线形成的向量来帮我们判断下一步的方向以及下一步所到的点。
1 给定起点和终点B的椭球极坐标,,A初始值设为起点坐标,计算其直角坐标,并计算B的外法向量记为BW
2 计算A点在二维极坐标中以A为中心的一定宽度方格内的点Mi的极坐标
3 计算所有Mi的直角坐标,并且计算出向量AMi,并将其全部化为单位向量
4 计算BW*AMi的值,并找出使该值最小的Mi点作为下一个A点
5 计算当前A点与B的欧式距离,若该距离小于ϵ \epsilonϵ ,则停止循环,输出D
6 计算下一个点A于当前点A之间的欧式距离di _{i}i,D=D+di _{i}i
代码实现
椭球坐标转化为直角坐标函数:
function [x,y,z]=elip2cart(theta,alpha)
a=6000;
b=5000;
x=a.*sin(theta).*sin(alpha);
y=a.*cos(theta).*sin(alpha);
z=b.*cos(alpha);
end
主程序:
clear all
a=6000;
b=5000;
%x=[-2200 2900];
%y=[-3600 3300];
%z=b*sqrt(1-(x.*x+y.*y)/(a*a));
%comval=zeros(9,33);
pcs=600;%网格精度,不能超过3000,否则电脑会无法存储过多的网格坐标数据而崩溃
cho=8;%AMi的选择个数控制,选择个数为(2*cho+1)^2;设置的越高,曲线就越平滑,但是可能出现跨度比较大的点,必须为正整数,一般设为6以上
g=linspace(0,pi,pcs);h=linspace(0,2*pi,pcs*2);
thetaS=g(pcs*0.3);alphaS=h(pcs*1.3);
thetaP=g(pcs*1.0);alphaP=h(pcs*1.9);%给出S,T的极坐标
sphS=[thetaS,alphaS];%路径点集
tL=g(2)-g(1);aL=tL;
[theta,alpha]=meshgrid(g,h);
[x,y,z]=elip2cart(theta,alpha);%计算椭圆上所有点的直角坐标
[xs, ys ,zs]=elip2cart(thetaS,alphaS);
[xt ,yt ,zt]=elip2cart(thetaP,alphaP); %计算S, T的直角坐标
pw=[xt*2/a^2,yt*2/a^2,zt*2/b^2];%计算T点的外法向量
pw=pw/norm(pw);
num=2;
lasts1=0;lasts2=0;
k=0;%距离终点的欧式距离
while num<10000
%----------------------构造所有的AMi的选择向量的指向点-------------------------------
m=thetaS-cho*tL:tL:thetaS+cho*tL;
n=alphaS-cho*aL:aL:alphaS+cho*aL;%生成网格节点
[moveSt,moveSa]=meshgrid(m,n);
moveSt=reshape(moveSt,[1,(cho*2+1)^2]);
moveSa=reshape(moveSa,[1,(cho*2+1)^2]);
[xmove,ymove,zmove]=elip2cart(moveSt,moveSa);%通过自己写的elip2cart函数将极坐标转为直角坐标
%---------------------进行核心算法判断,来选择椭圆表面移动点的下一个位置
sm=zeros((cho*2+1)^2,3);
cval=zeros(1,(cho*2+1)^2);
for i=1:(cho*2+1)^2
sm(i,:)=[xmove(i)-xs,ymove(i)-ys,zmove(i)-zs]/norm([xmove(i)-xs,ymove(i)-ys,zmove(i)-zs]);
cval(i)=sm(i,:)*pw';%计算每个AMi与PW的乘积
end
nextp=find(cval==max(cval)); %按照设计好的算法来选出最大乘积的AMi
thetaS=moveSt(nextp(1));alphaS=moveSa(nextp(1));%移动,将选出的M点坐标作为下一个起点
[xs, ys ,zs]=elip2cart(thetaS,alphaS);
sphS=[sphS;thetaS,alphaS];
num=num+1;
k=sqrt((xs-xt)^2+(ys-yt)^2+(zs-zt)^2);%计算当前移动后的点与终点的距离,若小于20(也可以自己设)则判读其到达了目的地
if k<20
break
end
end
thetaS=g(pcs*0.3);alphaS=h(pcs*1.3);%这个地方必须要跟前面设置的起点坐标保持一致
[xs, ys ,zs]=elip2cart(thetaS,alphaS);
[xr,yr,zr]=elip2cart(sphS(:,1),sphS(:,2));
td=0;%总距离
for i=2:length(xr)
td=td+sqrt((xr(i)-xr(i-1))^2+((yr(i)-yr(i-1))^2+((zr(i)-zr(i-1))^2)));
end
%comval(cho-1,(pcs-300)/100)=td;
td
realval=17314.47799
mesh(x,y,z);
hold on
plot3(xs, ys ,zs,'r.','markersize',24,'color','blue')
plot3(xt ,yt ,zt,'r.','markersize',24,'color','blue')
plot3(xr,yr,zr,'r.','markersize',5)
%,1:length(comval),realval.*ones(1,length(comval))
运行结果及结果分析
可以从图中看出,基本上已经实现了椭球上最短距离的寻找。接下来选取两个特殊点进行结果的验证和分析。
选取椭球的最高点到最低点的作为起点和终点,根据现有的理论近似公式计算出理论近似值,将其与通过该方法计算出的值进行比较。
理论计算的公式为计算椭圆的周长公式,名为Ramanujan近似公式。
计算得到的理论值为17314.47799,如上图所示,我们通过该方法计算得到的值为1.7283*104,与真实值比较接近。误差的来源可能为每次计算的欧式距离总是小于实际上的弧长,也可能是理论计算公式的带来的近似误差。
当改变网格的精度时,会使结果与理论计算值之间的差异发生变化,计算得到的值近似收敛于某个定值。推测收敛到这个值可能即为最短距离的真实值。
当改变AMi的可选择方向数目时对结果的影响不大。
优缺点及改进方向
本次实验较好地实现了通过计算机求解椭球面上两点的最短距离,得出的结果与理论近似公式得出的结果非常接近,对于椭球上的任意两点能够比较灵敏地自动寻找轨迹。所运用的方法是根据现有的梯度下降法的思想,方法创新性较强。
对于两个影响指标进行了分析,得到了这两个指标对结果的影响,并且给出了可能的解释。对于其中一个影响指标的趋势给出了可能的猜测。
本实验较大的不足在于,主要的算法没有进行严格的数学证明,是仅凭主观意识得出的结论,虽然结果比较理想但是没有严格证明结果的可靠性和正确性。
核心算法有漏洞,可能会导致特殊情况下路径选择的偏离。具体为,当向量AMi的选择范围较大时,比如大于100个方向向量,其选择的方向可能导致前后两个运动点间的跨度过大,从而导致路径偏离和较大的误差。如图所示。
matlab求椭圆的弧长,用MATLAB实现求椭球上任意两点的最短弧长相关推荐
- C语言入门实战(2):求平面上任意两点之间的距离
这是<C语言入门实战>系列的第2篇. 上一篇:C语言入门实战(1):准备开发环境.快速上手main()函数 下一篇:C语言入门实战(3):秒数转换为时:分:秒 文章目录 题目 提示 参考代 ...
- matlab 最小的两个,测定两个超椭球表面之间的最小距离
Determination of the minimum distance between two SuperEllipsoids surfaces. Using Optimization What´ ...
- JOJ2737:狼与羊的故事(求图上任意两点间的桥边)
2737: 狼与羊的故事 Result TIME Limit MEMORY Limit Run Times AC Times JUDGE 3s 65536K 53 10 Standard 村长要召开 ...
- 已知圆上任意两点求圆心和半径_已知圆上三点坐标求圆心和半径
R半径 PCenter圆点坐标 public void GetCircular(PointF P1,PointF P2,PointF P3,ref float R,ref PointF PCenter ...
- 任意两点最短路floyd算法matlab,多源最短路——Floyd算法
Floyd算法 问题的提出:已知一个有向网(或者无向网),对每一对定点vi!=vj,要求求出vi与vj之间的最短路径和最短路径的长度. 解决该问题有以下两种方法: (1)轮流以每一个定点为源点,重复执 ...
- matlab 椭圆方程拟合,matlab中如何插值拟合求椭圆方程
[g_fitting.rar] 使用正交多项式完成数据拟合.程序对读入的gps采样点完成曲线拟合. (2007-08-01, matlab, 1KB, 26次) [曲面拟合.rar] 这是利用matl ...
- matlab 变成圆形坐标,求圆和椭圆上任意角度的点的坐标
圆上任意角度的点的坐标 如上图,给定圆心(Cx,Cy),半径为R, 求θ \thetaθ对应的点的坐标? 此处θ \thetaθ是相对于水平轴的角度. 显然我们可以使用极坐标转换来求: { p x = ...
- MATLAB蒙特卡洛(Monte Carlo)方法求椭圆面积
MATLAB蒙特卡洛方法求椭圆面积 代码 代码 在某个规定的范围内随机打点,找到满足条件的点,并数一下这些点的数量与总的随机点数量的比,就OK了.关键是设置条件. 代码 clear;clc; n=10 ...
- matlab画椭圆 长轴 短轴,跟踪目标的快速椭圆拟合方法
摘 要: 提出一种基于最小外包矩形的快速椭圆拟合方法,该方法利用最小二乘法获得目标的最小外包矩形框,再求取外包矩形框的内切椭圆,该椭圆能有效反映目标的大部分运动信息.本文对该方法进行了目标拟合的有效 ...
最新文章
- 薛其坤、向涛两位院士,担任这家研究院联合院长
- dede使用方法----如何调用最新文章,最热文章,友情链接
- 嵌入式小知识(累积更新)
- 关于工大瑞普Dynamips模拟器
- Codeforces Round #521 (Div.3)题解
- 短视频时代不可忽视的幕后功臣竟然是它!
- linux libusb应用实例,在Linux中使用libusb-1.0作为非root用户访问USB设备
- Error running ‘Tomcat‘: Unable to open debugger port (127.0.0.1:2148): java.net.SocketExceptio
- Julia: arrow,一种革命性的数据格式
- GJB150-2009军用装备实验室环境试验方法新版标准
- bartender连接oracle,设置BarTender打印用的数据库连接
- 无线扫码枪 服务器查询异常,无线扫描枪常见问题及解决方法
- 对话机器人(一)——对话机器人基础知识
- c++学习:多态案例之魔法门英雄无敌
- java钝化_黑马day14 监听器之javaBean对象的活化和钝化
- ROS PGM格式文件详解 | 九七的ROS
- 北斗卫星短信通信与定位详解
- 第一次开发EOS区块链的经验总结
- 【408数据结构】备考常见必会算法图鉴
- 社交消费时代到来,趣享付占领社交营销一线