任意多边形费马点点群中位中心求解
点群中位中心求解
1问题讨论
从上两式不难看出，空间点（x0,y0）可以取空间中任意一点有无穷种可能，因此难以求得精确。只能使用迭代的方式来求无限逼近的解，当然求用一些奇特的方法求近似解也是可以的。
更加一般的，点群的中位中心问题是p-median问题d特例，即1-median问题。P-中值问题可以描述。而同样的问题放在多边形中，就变成了费马点的求解。
各个GIS工具软件中均有求中位中心的功能，而在matlab中也可以通过最优化函数fminsearch直接求解。本文接下来讨论点群中位中心求解的几个算法，并在matlab环境下给出其实现。
本文主要从底层讨论三种算法并给出其实现，同时指出其他算法可能性。
2算法
2.1 二分互垂近似算法
基本思想为：做一组互垂的直线，两条直线均把把点群分成等数量两部分。取多组这样的直线，每组直线交点的平均中心即为点群中位中心的近似。
算法步骤为：
STEP 1: 将点群副本绕原点顺时针旋转k*90o/n，横纵两个方向，用二分法扫描点群，找到横纵方向能把点群分为两个等数量部分直线，将直线交点逆时针旋转k*90o/n,记录旋转后的点;( k指旋转的次数,初始化为0。n-1指最大旋转次数也代表着求解精度。所谓二分法就是，根据横轴或纵轴区间对半切，判断那边的点多，对多的那一边再对半切，直到切完后直线两边点的数量相等)
STEP 2: 判断k是否不大于n-1，是则k++返回STEP1，否则输出交点的平均中心。
2.2 网格点近似算法
基本思想为：将点群所在的空间划分成小的网格，所有网格点中离点群总距离最近的点就是要求的点群中位中心的近似。
算法步骤为：
STEP 1: 扫描点群，找出其四至；
STEP 2: 根据输入精度n，计算点群四至矩形范围内均匀分布的格网点坐标；
STEP 3: 计算每个网格点与点群的总距离，取总距离最小的网格点作为近似的中位中心。
2.3 四方向变步长贪婪算法
基本原理：点群的中位中心和平均中心其实相隔的并不远，可以点群平均中心为起始点，以一定步长向四周移动，从中选择最优的点。如果当前点比所有移动点更优，说明移动的过分了，这时需要缩小步长继续移动搜索，直到精度达到要求。
算法步骤：
STEP 1：计算点群平均中心，作为起始点pt。计算点群的横纵方向的跨度，取其小者的一半作为起始步长；
STEP 2：以step为步长，向上下左右四个方向移动pt，得到pt1-4；
STEP 3：计算pt，pt1-4与点群的总距离，设总距离最小者为新的pt；
STEP 4: 判断进度是否达到要求，是则输出pt，否则继续；
STEP 5：如果当前pt与上一步pt相同，则令step=step*factor，其中factor是步长变化系数，可以取0.5等。返回第二步。
2.4 其他算法
除了上述三种算法外，此问题还可以使用模拟退火算法或牛顿迭代法进行求解。相关代码思路，在参考文献中已经给出。限于时间精力，这里不再进一步讨论。关于上述三种算法的matlab完整代码在附录中给出。
3算例
以全国300多个主要城市的经纬度坐标作为输入，利用上述三种算法分别计算，得到的结果如下图所示。
各个结果的对比如下表所示，可以看出不同算法差别较小，精度较好：
参考文献
[1]Matlab选址问题[DB/OL].20170312:
[2]空间分析,郭仁忠[M].北京：高等教育出版社.p81-p82
[3]选址问题及其模型与研究[DB/OL].20170312:http://www.docin.com/p-1328378124.html
[4]白话空间统计番外篇：中位数中心算法[DB/OL].20170312:
http://blog.csdn.net/allenlu2008/article/details/47752943
[5]模拟退火求二维费马点[DB/OL].20170312:
http://www.cnblogs.com/hxsyl/archive/2013/08/12/3252773.html
https://wenku.baidu.com/view/ddac681fc281e53a5902ff05.html
[7]求n个点的费马点的花式乱搞[DB/OL].20170312:
http://blog.csdn.net/wjf_wzzc/article/details/24126257
[8]二维费马点[DB/OL].20170312:http://java-mans.iteye.com/blog/1644473
%点群的中位中心计算程序，武汉大学孙一璠于20170312实现
function mediancenter= medianCenter_calcu(pts,method,isdisplay)
%主函数，用于计算点群的中位中心
%pts为输入点群[x y]；
%method 为所采用的算法；
%isdisplay为逻辑值，1则显示计算过程，2不显示计算过程；
if method==1mediancenter=in_linediv_calcu(pts,15,isdisplay);
end
if method==2mediancenter=in_gridpt_calcu(pts,100,isdisplay);
end
if method==3mediancenter=in_stepwise_calcu(pts,0.05,isdisplay);
end
endfunction mediancenter=in_linediv_calcu(pts,n,isdisplay)
%算法1：通过互垂分割先近似求点群中位中心
%pts为输入点群，n为直线旋转次数（精度）
linepts=[0 0];resultpts=[0 0];
for k=0:n-1
pts2=in_rotate(pts,[0 0],k*pi/(2*n));
[up,bottom,left,right]=in_fideborder(pts2);
y=(up+bottom)/2;
[upnumber,downnumber]=in_calcudiv_y(pts2,y);
x=(left+right)/2;
[leftnumber,rightnumber]=in_calcudiv_x(pts2,x);
%二分法迭代求等分线,这个地方要设置一个0.05做为不再二分的条件，防止运算爆炸
while upnumber~=downnumber && abs(up-bottom)>=0.05if upnumber>downnumberbottom=y;y=(up+bottom)/2;[upnumber,downnumber]=in_calcudiv_y(pts2,y); else up=y;y=(up+bottom)/2;[upnumber,downnumber]=in_calcudiv_y(pts2,y); end
end
while leftnumber~=rightnumber && abs(left-right)>=0.05if leftnumber>rightnumberright=x;x=(left+right)/2;[leftnumber,rightnumber]=in_calcudiv_x(pts2,x);endif leftnumber<rightnumberleft=x;x=(left+right)/2;[leftnumber,rightnumber]=in_calcudiv_x(pts2,x);end
end
[up,bottom,left,right]=in_fideborder(pts2);
linepts_temp=in_rotate([left y; right y;x up;x bottom],[0 0],-k*pi/(2*n));
resultpts_temp=in_rotate([x y],[0 0],-k*pi/(2*n));
linepts=[linepts;linepts_temp];
resultpts=[resultpts;resultpts_temp];
end
linepts=linepts(2:end,:);
resultpts=resultpts(2:end,:);
mediancenter=in_meanCenter_calcu(resultpts);%计算过程与计算结果可视化if(isdisplay==true)figure;hold on;axis equal;[m,~]=size(linepts);for i=1:m/4plot(linepts((4*i-3):(4*i-2),1),linepts((4*i-3):(4*i-2),2),'linewidth',0.1,'linestyle',':','color','y');plot(linepts((4*i-1):(4*i),1),linepts((4*i-1):(4*i),2),'linewidth',0.1,'linestyle',':','color','y');endplot(pts(:,1),pts(:,2),'linestyle','none','marker','O','markeredgecolor','k','markersize',3);plot(resultpts(:,1),resultpts(:,2),'linestyle','none','marker','.','markeredgecolor','b');plot(mediancenter(1),mediancenter(2),'linestyle','none','marker','*','markeredgecolor','r');end
endfunction mediancenter=in_gridpt_calcu(pts,n,isdisplay)
%通过网格点来近似和点群的中位中心
%pts是输入点群，n表示网格精度，isdisplay表示是否展示计算过程和结果
[up,bottom,left,right]=in_fideborder(pts);
xspace=(right-left)/n;
yspace=(up-bottom)/n;
gridpts=ones((n+1)^2,2);
for i=1:(n+1)for j=1:(n+1)gridpts((i-1)*(n+1)+j,1)=left+xspace*(i-1);gridpts((i-1)*(n+1)+j,2)=bottom+yspace*(j-1);end
end
distance=in_totaldis_calcu(pts,gridpts);
[~,index]=min(distance);
mediancenter=gridpts(index,:);%绘制计算结果if(isdisplay==true)figure;hold on;axis equal;plot(gridpts(:,1),gridpts(:,2),'linestyle','none','marker','.','markersize',1);plot(pts(:,1),pts(:,2),'linestyle','none','marker','O','markersize',3,'markeredgecolor','k');plot(mediancenter(:,1),mediancenter(:,2),'linestyle','none','marker','*','markeredgecolor','r');end
endfunction mediancenter=in_stepwise_calcu(pts,steplength,isdisplay)
%四方向变步长贪心算法
%其原理是：点群的中位中心和平均中心其实相隔的并不远，可以点群平均中心为起始点，以一定步长向四周移动，从中选择最优的点。如果当前点比所有移动点更优，说明移动的过分了，这时需要缩小步长继续移动搜索，直到精度达到要求。
[up,bottom,left,right]=in_fideborder(pts);
%定义初始步长step和步长变化因子factor,计算点群的平均中心，初始化搜索路径序列
step=0.5*min([up-bottom right-left]);
factor=0.5;
pt=in_meanCenter_calcu(pts);
roadpts=pt;
while step>steplengthpt1=[pt(1) pt(2)+step];pt2=[pt(1)+step pt(2)];pt3=[pt(1) pt(2)-step];pt4=[pt(1)-step pt(2)];searchpts=[pt;pt1;pt2;pt3;pt4];distance=in_totaldis_calcu(pts,searchpts);[~,index]=min(distance);if index==1step=step*factor;elsept=searchpts(index,:);roadpts=[roadpts;pt];end
end
mediancenter=pt;%绘制结果与过程
if(isdisplay==true)figure;hold on;axis equal;plot(roadpts(:,1),roadpts(:,2),'linestyle','-','marker','.','color','y','markeredgecolor','b');plot(pts(:,1),pts(:,2),'linestyle','none','marker','O','markersize',3,'markeredgecolor','k');plot(mediancenter(:,1),mediancenter(:,2),'linestyle','none','marker','*','markeredgecolor','r');
end
endfunction distance=in_totaldis_calcu(pts,inputpts)
%计算输入点到点群的总距离
[m,~]=size(inputpts);
distance=ones(m,1);
for i=1:m
singledis=sqrt((inputpts(i,1)-pts(:,1)).^2+(inputpts(i,2)-pts(:,2)).^2);
distance(i)=sum(singledis);
end
endfunction [up,bottom,left,right]=in_fideborder(pts)
%寻找点群四至值
minvalue=min(pts);
left=minvalue(1);bottom=minvalue(2);
maxvalue=max(pts);
right=maxvalue(1);up=maxvalue(2);
endfunction [upnumber,downnumber]=in_calcudiv_y(pts,y)
%计算水平线上下分别有多少个点
[m,~]=size(pts);
count=0;
for i=1:mif pts(i,2)>=ycount=count+1;end
end
upnumber=count;
downnumber=m-count;
endfunction [leftnumber,rightnumber]=in_calcudiv_x(pts,x)
%计算竖直线左右分别有多少个点
[m,~]=size(pts);
count=0;
for i=1:mif pts(i,1)<=xcount=count+1;end
end
leftnumber=count;
rightnumber=m-count;
endfunction pts2=in_rotate(pts,center,theta)
%将点群pts绕center点逆时针旋转theta度（弧度制）
x=pts(:,1);y=pts(:,2);
x0=center(1);y0=center(2);
x2=(x-x0).*cos(theta)-(y-y0).*sin(theta);
y2=(x-x0).*sin(theta)+(y-y0).*cos(theta);
pts2=[x2 y2];
endfunction meancenter=in_meanCenter_calcu(pts)
%计算点群的平均中心
[m,~]=size(pts);
x=sum(pts(:,1))/m;
y=sum(pts(:,2))/m;
meancenter=[x y];
end
任意多边形费马点点群中位中心求解相关推荐
- 基于matlab 求多边费马点,POJ2420（求多边形费马点） | 学步园
题意:题目的意思就是给你N个点,在平面上寻找一个点,使得这个点到其他点的距离之和最小,问你最小的距离是多 少? 分析:在三角形内部这个点叫做费马点(费马点定义).那么这道题目就是求一个多变形的费马点. ...
- POJ2420（求多边形费马点）
题目:题目链接 题意:题目的意思就是给你N个点,在平面上寻找一个点,使得这个点到其他点的距离之和最小,问你最小的距离是多 少? 分析:在三角形内部这个点叫做费马点(费马点定义).那么这道题目就是求一个 ...
- c语言中余数恒等于1,费马小定理_KANGMANG201102_新浪博客
费马小定理是数论中的一个重要定理,其内容为: 假如p是质数,且(a,p)=1,那么 a^(p-1) ≡1(mod p) 假如p是质数,且a,p互质,那么 a的(p-1)次方除以p的余数恒等于1 费马小 ...
- 神奇的交际圈！这位17世纪的法国神父结交的好朋友，竟然都是一流数学牛人：笛卡尔、费马、加森迪······
全世界只有3.14 % 的人关注了 爆炸吧知识 话说,在近代数学史上,人们惊讶地发现17至18世纪的法国竟然产生了众多一流的数学家. 然而,最早想到要培养一波优秀人才的,成就这段群星璀璨的传奇历史,并 ...
- c语言生成两位随机素数算法,[算法]费马小定理求质数的算法之Miller-Rabin算法，C语言实现 | 李大仁博客...
今天讲点比较高级的算法,目的也很简单,求质数,但是应用一种新的算法Miller-Rabin算法,这是一种利用了概率和费马小定理的算法设计,有点玄乎吧,其实本人也是刚接触这种算法,这是一种纯数学的解法, ...
- python计算多边形的面积并保留两位小数_计算任意多边形面积的Python实现
最近需要实现一个计算非凸多边形面积的功能,需要输入是顺次排序的多边形顶点坐标,假设输入的多边形顶点是V={v0, v1, v2, -, vn-1},则多边形的边为E={, , ,...,, }.要求输 ...
- 费马小定理在ACM中的应用
费马小定理 假如p是素数,且 (a,p)=1 ( a , p ) = 1 (a,p)=1,那么a^(p-1)≡1(mod p) ① 判断素数,对于大素数的判定,Miller-Rabin 素数判定 ②求 ...
- 循环小数与费马小定理
循环小数与费马小定理 17/05/29 22:30:51 | Snakes 背景 题目出自之前亮灯问题.杨辉三角与Sierpinski三角形提及的生日题中的第三.四.五题. 题目 第三题 证明:对于任 ...
- python 费马检测
费马小定理: 如果 n 是一个素数,a 是小于 n 的任意正整数,那么 a 的 n 次方与 a 模 n 同余.(两个数称为是模 n 的同余,如果它们除以 n 的余数相同.数 a 除以 n 的余数称为 ...
最新文章
- linux salt命令 -e,linux 下 Salt 命令的疑难杂症
- python 调用class不指定函数_python调用另一个.py中的类或函数
- 互联网员工桌子上的药
- 初识ABP vNext（1）：开篇计划基础知识
- Z-Stack通过按键中断实现长按功能
- 怎么自学python-你是如何自学 Python 的？
- 高速掌握Lua 5.3 —— 扩展你的程序 (1)
- WINDOWS调用出错后，得到信息字串
- 细数Android原生工程接入EasyAR-SurfaceTracking遇到的坑
- Amesim学习——弹球仿真
- python if实现对话_Python生成微信对话生成器（四）
- 如何学好高中数学 提高高中数学成绩秒杀技巧（这几点很重要）
- 怎么学习CAD？初学CAD如何入门
- Python爬虫练习：爬取软科世界大学学术排名
- 红旗 linux 在哪儿 看 版本,简介红旗Linux不同版本的介绍
- 【调参15】如何配置神经网络的学习率
- 爬取（明星网）明星面部数据
- 利用RAW格式处理大光比照片
- Mac 下Qt 设计师模式下菜单栏不显示问题
- 缓冲区溢出攻击（Buffer Overflows实验笔记）