找回密码
 注册
Simdroid-非首页
查看: 297|回复: 15

[编程进阶] 一道代码测试题: N 点最小内接圆问题

[复制链接]
发表于 2010-4-21 00:14:26 | 显示全部楼层 |阅读模式 来自 北京
本帖最后由 FreddyMusic 于 2010-4-23 18:54 编辑

有二十个随机点,存在很多圆可以使得所有的点都在圆内(包括圆周上).求满足条件的最小的圆的圆心和半径,当然满足条件的圆可能不止一个,大家看看如何找到其中一个或多个最小的圆.
  1. (*为了简单期间,下面用来测试的坐标均为整数值*)
  2. weizhi={{99, 2}, {43, 64}, {44, 28}, {80, 0}, {28, 94}, {8, 15}, {12,
  3.   89}, {1, 30}, {81, 41}, {10, 3}, {68, 11}, {78, 43}, {7, 17}, {64,
  4.   99}, {67, 11}, {65, 19}, {61, 41}, {27, 68}, {76, 73}, {45, 79}}
复制代码

评分

1

查看全部评分

发表于 2010-4-21 01:55:11 | 显示全部楼层 来自 黑龙江哈尔滨
Simdroid开发平台
你想测试什么?

这里有个类似的问题,用mmtc的mathoptpropackage做的。

http://library.wolfram.com/infocenter/TechNotes/6202/
回复 不支持

使用道具 举报

发表于 2010-4-21 02:13:51 | 显示全部楼层 来自 黑龙江哈尔滨
你想测试什么?

这里有个类似的问题,用mmtc的mathoptpropackage做的。

http://library.wolfram.com/infocenter/TechNotes/6202/
TBE_Legend 发表于 2010-4-21 01:55

本帖子中包含更多资源

您需要 登录 才可以下载或查看,没有账号?注册

×
回复 不支持

使用道具 举报

发表于 2010-4-21 08:28:47 | 显示全部楼层 来自 北京海淀
本帖最后由 shamohu 于 2010-4-22 15:20 编辑

用1stOpt试了下,似乎很简单:

  1. Constant x=[99,43,44,80,28,8,12,1,81,10,68,78,7,64,67,65,61,27,76,45],
  2.          y=[2,64,28,0,94,15,89,30,41,3,11,43,17,99,11,19,41,68,73,79];
  3. MinFunction R;                                            //目标函数:半径最小
  4. For(j=1:20)(sqrt((x0-x[j])^2+(y0-y[j])^2)<=R);//约束:圆心到各点的距离不大于半径
复制代码
结果:
目标函数值(最小): 61.522731470598
R: 61.522731470598
x0: 54.9772727272726
y0: 44.9772727272726
回复 不支持

使用道具 举报

发表于 2010-4-21 11:30:28 | 显示全部楼层 来自 江苏苏州
g 的题目出的不错,继续出题,由浅入深,引导会员。继续多出题。

T 的动作够快,反应敏捷,找的方向还是对的,但是我不知这个packing是Wolfram 开发的还是外部人员开发的。能否看见所有的源码?我记得Demo 上也有类似 packing 的演示。

s 有解不错,下一步能否可视化,或者一般化?
回复 不支持

使用道具 举报

 楼主| 发表于 2010-4-21 12:08:23 | 显示全部楼层 来自 北京
3# TBE_Legend
"代码测试"
(1)是测试你对该问题的解题思路;
(2)是测试按照你的思路写出的代码的运算时间的.
我比较倾向于用用最基本的函数实现想要实现的功能.

当然所谓的最基本也是因人而异的,
至少我不希望动不动就使用Reduce函数,尽管有时候觉得它的功能的确强大,运算时间也短,但是我无法理解它的计算原理,所以我并不将它看作是最基本的函数.

计算机如何计算两数的乘积也是我无法理解的,但是我还是将它看作是基本的函数.

当然处理实际问题时,我还是很乐意用所有可以用到的函数的.
毕竟我是业余爱好者,和你们的想法未必相同.
回复 不支持

使用道具 举报

 楼主| 发表于 2010-4-21 12:13:45 | 显示全部楼层 来自 北京
本帖最后由 ggggwhw 于 2010-4-21 16:16 编辑

4# shamohu
shamohu同志的代码倒是很短啊,结果也很精确,但是很可惜,我无法验证结果,因为我没有使用过1stop,也不知道速度如何,据说1stop的运算速度是很快的.
我更想知道,这个圆心点是如何找到的,能用语言描述一下它的寻找过程吗?
肯定不是穷举平面内所有点比较得到的哦,呵呵

我在8楼给出的Mathematica 代码给出了精确值,
计算结果是:
圆心{2419/44, 1979/44}≈{54.97727272727273`, 44.97727272727273`}
半径(5 Sqrt[146557/2])/22≈61.52273147059796`

可见1stop还真是厉害啊.
回复 不支持

使用道具 举报

 楼主| 发表于 2010-4-21 18:01:00 | 显示全部楼层 来自 北京
本帖最后由 ggggwhw 于 2010-4-21 18:14 编辑

虽然代码还没有完成, 但是现在已经出来精确解了, 下面是基本思路, 先传上来大家分享一下, 我需要一些时间完善一下, 现在需要改善的是再添几个判断, 如果还不行的话, 就要加上循环语句了.
  1. weizhi = {{99, 2}, {43, 64}, {44, 28}, {80, 0}, {28,
  2.    94}, {8, 15}, {12, 89}, {1, 30}, {81, 41}, {10, 3}, {68, 11}, {78,
  3.    43}, {7, 17}, {64, 99}, {67, 11}, {65, 19}, {61, 41}, {27,
  4.    68}, {76, 73}, {45, 79}}
  5. (*预定义三角形外心公式*)
  6. waixin[{x1_, y1_}, {x2_, y2_}, {x3_,
  7.     y3_}] := {(x2^2 y1 - x3^2 y1 - x1^2 y2 + x3^2 y2 - y1^2 y2 +
  8.       y1 y2^2 + x1^2 y3 - x2^2 y3 + y1^2 y3 - y2^2 y3 - y1 y3^2 +
  9.       y2 y3^2)/(2 (x2 y1 - x3 y1 - x1 y2 + x3 y2 + x1 y3 -
  10.         x2 y3)), -(-x1^2 x2 + x1 x2^2 + x1^2 x3 - x2^2 x3 - x1 x3^2 +
  11.        x2 x3^2 - x2 y1^2 + x3 y1^2 + x1 y2^2 - x3 y2^2 - x1 y3^2 +
  12.        x2 y3^2)/(2 (x2 y1 - x3 y1 - x1 y2 + x3 y2 + x1 y3 - x2 y3))};
  13. (*先将这二十个点的重心作为圆心的近似解*)
  14. yuanxin = Mean[weizhi];
  15. weizhi = Sort[weizhi, Norm[#1 - yuanxin] >= Norm[#2 - yuanxin] &];
  16. (*将离圆心最远的一个点作为最小圆上的一个点doc[1]*)
  17. doc[1] = weizhi[[1]];
  18. banjing = Norm[doc[1] - yuanxin];
  19. (*作图,第一个近似最小圆*)
  20. tu[0] = ListPlot[{yuanxin}, PlotStyle -> Directive[Red]];
  21. tu[1] = ListPlot[weizhi, PlotStyle -> Directive[Blue]];
  22. tu[2] = Graphics[Circle[yuanxin, banjing]];
  23. Show[Table[tu[j], {j, 0, 2}],
  24. AxesOrigin -> {0, 0},
  25. AspectRatio -> Automatic,
  26. PlotRange ->
  27.   Transpose@{yuanxin - banjing - 10, yuanxin + banjing + 10}]
  28. "第一个近似最小圆的圆心:"
  29. yuanxin
  30. yuanxin // N
  31. "第一个近似最小圆的半径:"
  32. banjing
  33. banjing // N
  34. (*将离doc[1]最远的一个点作为最小圆上的另一个点doc[2]*)
  35. weizhi = Sort[weizhi, Norm[#1 - doc[1]] >= Norm[#2 - doc[1]] &];
  36. doc[2] = weizhi[[1]];
  37. (*将线段 {doc[1],doc[2]}近似看作最小圆的一条直径*)
  38. banjing = Norm[doc[1] - doc[2]]/2;
  39. yuanxin = (doc[1] + doc[2])/2;
  40. (*作图,第二个近似最小圆*)
  41. tu[0] = ListPlot[{yuanxin}, PlotStyle -> Directive[Red]];
  42. tu[1] = ListPlot[weizhi, PlotStyle -> Directive[Blue]];
  43. tu[2] = Graphics[Circle[yuanxin, banjing]];
  44. Show[Table[tu[j], {j, 0, 2}],
  45. AxesOrigin -> {0, 0},
  46. AspectRatio -> Automatic,
  47. PlotRange ->
  48.   Transpose@{yuanxin - banjing - 10, yuanxin + banjing + 10}]
  49. "第二个近似最小圆的圆心:"
  50. yuanxin
  51. yuanxin // N
  52. "第二个近似最小圆的半径:"
  53. banjing
  54. banjing // N
  55. (*将离圆心最远的点作为最小圆上的第三个点doc[3]*)
  56. weizhi = Sort[weizhi, Norm[#1 - yuanxin] >= Norm[#2 - yuanxin] &];
  57. doc[3] = weizhi[[1]];
  58. If[! MemberQ[{doc[1], doc[2]}, doc[3]],
  59.   (*将三角形 {doc[1],doc[2],doc[3]}的外接圆作为最小圆*)
  60.   yuanxin = waixin[doc[1], doc[2], doc[3]];
  61.   banjing = Norm[doc[1] - yuanxin]
  62.   ];
  63. tu[0] = ListPlot[{yuanxin}, PlotStyle -> Directive[Red]];
  64. tu[1] = ListPlot[weizhi, PlotStyle -> Directive[Blue]];
  65. tu[2] = Graphics[Circle[yuanxin, banjing]];
  66. Show[Table[tu[j], {j, 0, 2}],
  67. AxesOrigin -> {0, 0},
  68. AspectRatio -> Automatic,
  69. PlotRange ->
  70.   Transpose@{yuanxin - banjing - 10, yuanxin + banjing + 10}]
  71. "满足条件的最小圆的圆心:"
  72. yuanxin
  73. yuanxin // N
  74. "满足条件的最小圆的半径:"
  75. banjing
  76. banjing // N
复制代码
回复 不支持

使用道具 举报

发表于 2010-4-21 21:52:05 | 显示全部楼层 来自 上海宝山区
g,

s 的代码和方法还是不错的,以最小化半径/面积为目标,
1opt 的代码仔细耐心的看你可以看懂。
回复 不支持

使用道具 举报

 楼主| 发表于 2010-4-22 20:09:04 | 显示全部楼层 来自 北京
本帖最后由 ggggwhw 于 2010-4-22 20:24 编辑

下面的代码是将8楼的代码的正确率提高后的结果,
该代码的优点是计算时间较短,如果结果正确,那么结果将是精确值,而不是近似值.
缺点是正确率小于100%(可以确定正确率大于50%,具体是多少我没有计算过,也没有实验测定过),于是后面给出了结果是否正确的说明.
  1. (*由随机函数得到20个随机点坐标*)"下面是20个随机点:"
  2. weizhi = Round[
  3. Table[{100*Cos[2 Pi*Random[]], 100*Sin[2 Pi*Random[]]}, {j, 1, 20}],
  4. 1/1000]
  5. lenw = Length[weizhi];
  6. (*预定义三角形外心公式*)
  7. waixin[{x1_, y1_}, {x2_, y2_}, {x3_,
  8. y3_}] := {(x2^2 y1 - x3^2 y1 - x1^2 y2 + x3^2 y2 - y1^2 y2 +
  9. y1 y2^2 + x1^2 y3 - x2^2 y3 + y1^2 y3 - y2^2 y3 - y1 y3^2 +
  10. y2 y3^2)/(2 (x2 y1 - x3 y1 - x1 y2 + x3 y2 + x1 y3 -
  11. x2 y3)), -(-x1^2 x2 + x1 x2^2 + x1^2 x3 - x2^2 x3 - x1 x3^2 +
  12. x2 x3^2 - x2 y1^2 + x3 y1^2 + x1 y2^2 - x3 y2^2 - x1 y3^2 +
  13. x2 y3^2)/(2 (x2 y1 - x3 y1 - x1 y2 + x3 y2 + x1 y3 - x2 y3))};
  14. jiaodu[{x_, y_}] := ArcTan[x[[1]] - y[[1]], x[[2]] - y[[2]]];
  15. (*先将这二十个点的重心作为圆心的近似解*)
  16. yuanxin = Mean[weizhi];
  17. weizhi = Sort[weizhi, Norm[#1 - yuanxin] >= Norm[#2 - yuanxin] &];
  18. (*将离圆心最远的一个点作为最小圆上的一个点doc1*)
  19. doc1 = weizhi[[1]];
  20. banjing = Norm[doc1 - yuanxin];
  21. (*将离doc1最远的一个点作为最小圆上的另一个点doc2*)
  22. weizhi = Sort[weizhi, Norm[#1 - doc1] >= Norm[#2 - doc1] &];
  23. doc2 = weizhi[[1]];
  24. (*将线段 {doc1,doc2}近似看作最小圆的一条直径*)
  25. banjing = Norm[doc1 - doc2]/2;
  26. yuanxin = (doc1 + doc2)/2;
  27. (*将离圆心最远的点作为最小圆上的第三个点doc3*)
  28. weizhi = Sort[weizhi, Norm[#1 - yuanxin] >= Norm[#2 - yuanxin] &];
  29. doc3 = weizhi[[1]];
  30. If[! MemberQ[{doc1, doc2}, doc3],(*将三角形 {doc1,doc2,doc3}的外接圆作为最小圆*)
  31. yuanxin = waixin[doc1, doc2, doc3];
  32. banjing = Norm[doc1 - yuanxin];];
  33. (*判断一下上面的外接圆是否满足条件,如果前面的外接圆不是最小的外接圆,寻找更合适的外接圆*)
  34. weizhi = Sort[weizhi, Norm[#1 - yuanxin] >= Norm[#2 - yuanxin] &];
  35. doc4 = weizhi[[1]];
  36. For[j = 1, j <= lenw, j++,
  37. weizhi = Sort[weizhi, Norm[#1 - yuanxin] >= Norm[#2 - yuanxin] &];
  38. doc4 = weizhi[[1]];
  39. If[! MemberQ[{doc1, doc2, doc3}, doc4],
  40. yuanxin = waixin[doc1, doc3, doc4];
  41. banjing = Norm[doc1 - yuanxin];
  42. If[Norm[doc2 - yuanxin] > banjing,
  43. yuanxin = waixin[doc2, doc3, doc4];
  44. banjing = Norm[doc2 - yuanxin];
  45. doc1 = doc4, doc2 = doc4;], Break[];];];
  46. (*判断一下最后的结果是否满足条件*)
  47. yuanhu = Select[weizhi, Norm[# - yuanxin] == banjing &];
  48. jiaodu1 =
  49. Table[ArcTan[yuanhu[[j, 1]] - yuanxin[[1]],
  50. yuanhu[[j, 2]] - yuanxin[[2]]], {j, 1, Length[yuanhu]}];
  51. jiaodu1 = Sort[jiaodu1, #1 >= #2 &];
  52. jiaodu2 = RotateRight[jiaodu1, 1];
  53. jiaodu2[[1]] = jiaodu2[[1]] + 2 Pi;
  54. jiaodu3 = jiaodu2 - jiaodu1;
  55. jiaodu3 = Sort[jiaodu3, #1 >= #2 &];
  56. If[j > lenw, Print["下面的结果是可疑的:"],
  57. If[jiaodu3[[1]] > 180*Degree, Print["下面的结果是错误的:"],
  58. Print["下面的结果是正确的:"];];];
  59. "满足条件的最小圆的圆心:"
  60. yuanxin
  61. yuanxin // N
  62. "满足条件的最小圆的半径:"
  63. banjing
  64. banjing // N
  65. tu[0] = ListPlot[{yuanxin}, PlotStyle -> Directive[Blue]];
  66. tu[1] = ListPlot[weizhi, PlotStyle -> Directive[Red]];
  67. tu[2] = Graphics[Circle[yuanxin, banjing]];
  68. Show[Table[tu[j], {j, 0, 2}], AxesOrigin -> {0, 0},
  69. AspectRatio -> Automatic,
  70. PlotRange ->
  71. Transpose@{yuanxin - banjing - 10, yuanxin + banjing + 10}]
复制代码
回复 不支持

使用道具 举报

发表于 2010-4-23 14:57:51 | 显示全部楼层 来自 上海宝山区
这题还是挺复杂的,有时候两点确定最小圆,有时候三点确定最小圆,优化方法同s。
回复 不支持

使用道具 举报

发表于 2010-4-23 18:54:22 | 显示全部楼层 来自 上海宝山区
本帖最后由 FreddyMusic 于 2010-4-23 18:58 编辑

Check this.
http://mathworld.wolfram.com/MinimalEnclosingCircle.html

http://www.cs.mcgill.ca/~cs507/projects/1998/jacob/problem.html

My style, at the moment.

  1. Manipulate[
  2. SeedRandom[seed];pts=RandomReal[{-1,1},{size,2}];
  3. (* solve *)
  4. equations=Table[(x-pts[][[1]])^2+(y-pts[][[2]])^2<= R^2,{i,Length@pts}];
  5. expr=Prepend[Prepend[equations,R>=0],R];
  6. {R0,solu}=Quiet@NMinimize[expr,{x,y,R},AccuracyGoal->4,PrecisionGoal->4,MaxIterations->100];
  7. {x0,y0,r}={x,y,R}/.solu;
  8. (* decide *)
  9. ptsToCenter=Sort[pts,((x0-#1[[1]])^2+(y0-#1[[2]])^2>=(x0-#2[[1]])^2+(y0-#2[[2]])^2)&];
  10. If[Length@pts==2,
  11. ptsOnCircle=Take[ptsToCenter,2],
  12. candidates=Take[ptsToCenter,3];
  13. distance=EuclideanDistance[{x0,y0},#]&/@candidates;
  14. error=10^(-5);
  15. ptsOnCircle=Take[candidates,If[(distance[[1]]-distance[[2]]<error)&&(distance[[2]]-distance[[3]]<error),3,2]]
  16. ];
  17. (* visulize *)
  18. Graphics[{
  19. If[show,{Cyan,Thickness[.01],Circle[{x0,y0},r],
  20. Pink,Thickness[.005],Dashed,Line[{{x0,y0},#}&/@ptsOnCircle],
  21. Red,PointSize[.02],Tooltip[Point[{x0,y0}],"circle center"]},{}],
  22. Gray,PointSize[.015],Point[pts],Tooltip[Point[ptsOnCircle],"points on circle"]
  23. },Axes->True,AxesLabel->{"x","y"},AxesOrigin->{0,0},PlotRange->{{-1.5,1.5},{-1.5,1.5}}, ImageSize-> 400],
  24. {{seed,12,"new random seed"},1,100,1,Appearance->"Labeled"},
  25. {{size,2,"number of points"},2,100,1,Appearance->"Labeled"},
  26. {{show,True,"show circle"},{False,True},ControlType->Checkbox},
  27. ControlPlacement->Top,TrackedSymbols->Manipulate,SaveDefinitions->True]
复制代码

本帖子中包含更多资源

您需要 登录 才可以下载或查看,没有账号?注册

×
回复 不支持

使用道具 举报

 楼主| 发表于 2010-4-25 19:00:35 | 显示全部楼层 来自 北京
F贴出来的长代码几乎每一次都无法运行,或许是他认为原样贴出会被人学去,或许……,
总之我很讨厌它.
回复 不支持

使用道具 举报

发表于 2010-4-25 20:50:31 | 显示全部楼层 来自 上海宝山区
“每次”啊?我有这么糟糕嘛?

看附件吧,代码老了。

本帖子中包含更多资源

您需要 登录 才可以下载或查看,没有账号?注册

×
回复 不支持

使用道具 举报

发表于 2010-4-25 21:26:26 | 显示全部楼层 来自 上海宝山区
这个问题是个经典的计算几何的问题。
推而广之,可以三维空间球体的包容,见我的另一个演示。


但是,我的算法和s的算法都是属于,数值计算的全局优化,算法的效率有一点,但都是比较低的。
换言之,我们的算法属于一种通用算法,还不是为这个问题量身定做的高效算法,或是基于严格的数学证明。

这个问题的也有较好的算法,网上有些论文可以找到,
但我也没仔细研究,究竟哪个算法才是经典高效。

你要是对此问题有新的认识和发现,经管发帖探讨。

本帖子中包含更多资源

您需要 登录 才可以下载或查看,没有账号?注册

×
回复 不支持

使用道具 举报

发表于 2010-4-28 15:45:48 | 显示全部楼层 来自 广西北海
去掉园内多余的点,形成凸多边形后再计算,速度要快不少。
回复 不支持

使用道具 举报

您需要登录后才可以回帖 登录 | 注册

本版积分规则

Archiver|小黑屋|联系我们|仿真互动网 ( 京ICP备15048925号-7 )

GMT+8, 2024-3-29 13:31 , Processed in 0.044988 second(s), 15 queries , Gzip On, MemCache On.

Powered by Discuz! X3.5 Licensed

© 2001-2024 Discuz! Team.

快速回复 返回顶部 返回列表