一、问题描述
给定平面
n
n
n个点，计算最小覆盖圆，使得所有给定点均在圆周以内(包括在圆周上)。博主commonc详细给出了最小覆盖圆算法步骤、算法时间复杂度及其证明、空间复杂度。博主Howe_Young书写的代码简洁、形式优雅，本文在其代码上做了小修改：(1)通过判断距离的平方而非距离，来判断点是否在圆内，从而避免平方根运算，有效提升计算速度 (2)尽量采用乘法代替除法，提升计算速度。下面是改进后的最小覆盖圆算法的MATLAB、LUA、C++完整代码。
2个点示例：
3个点示例：
10个点示例：
100个点示例：
二.MATLAB代码
%{
Function: get_sign
Description: 求实数x的符号
Input: 实数x
Output: 实数x的符号y
Author: Marc Pony(marc_pony@163.com)
%}
function y = get_sign(x)
if abs(x) < 1.0e-8
y = 0;
else
if x < 0.0
y = -1;
else
y = 1;
end
end
end
%{
Function: get_distance_square
Description: 求平面两点之间距离的平方
Input: 平面两点a，b
Output: 平面两点之间距离的平方
Author: Marc Pony(marc_pony@163.com)
%}
function distanceSquare = get_distance_square(a, b)
distanceSquare = (a(1) - b(1))^2 + (a(2) - b(2))^2;
end
%{
Function: get_circle_center
Description: 求三角形外接圆的圆心
Input: 平面三个点a，b，c
Output: 三角形外接圆的圆心center
Author: Marc Pony(marc_pony@163.com)
%}
function center = get_circle_center(a, b, c)
center = zeros(2, 1);
a1 = b(1) - a(1);
b1 = b(2) - a(2);
c1 = 0.5 * (a1 * a1 + b1 * b1);
a2 = c(1) - a(1);
b2 = c(2) - a(2);
c2 = 0.5 * (a2 * a2 + b2 * b2);
d = a1 * b2 - a2 * b1;
center(1) = a(1) + (c1 * b2 - c2 * b1) / d;
center(2) = a(2) + (a1 * c2 - a2 * c1) / d;
end
%{
Function: min_cover_circle
Description: 求平面pointCount个点的最小覆盖圆
Input: 平面pointCount个点的坐标(x,y)，点个数pointCount
Output: 平面pointCount个点的最小覆盖圆圆心center，半径radius
Author: Marc Pony(marc_pony@163.com)
%}
function [center, radius] = min_cover_circle(x, y, pointCount)
p = [x(:)'; y(:)'];
p = p(:, randperm(pointCount)); %随机打乱数据
center = p(:, 1);
radiusSquare = 0.0;
for i = 2 : pointCount
if get_sign(get_distance_square(p(:, i), center) - radiusSquare) > 0
center = p(:, i);
radiusSquare = 0.0;
for j = 1 : i - 1
if get_sign(get_distance_square(p(:, j), center) - radiusSquare) > 0
center = 0.5 * (p(:, i) + p(:, j));
radiusSquare = get_distance_square(p(:, j), center);
for k = 1 : j - 1
if get_sign(get_distance_square(p(:, k), center) - radiusSquare) > 0
center = get_circle_center(p(:, i), p(:, j), p(:, k));
radiusSquare = get_distance_square(p(:, i), center);
end
end
end
end
end
end
radius = sqrt(radiusSquare);
end
%{
Function: draw_circle
Description: 画圆周
Input: 圆心center，圆周半径radius，线型/颜色等设置参数options
Output: 无
Author: Marc Pony(marc_pony@163.com)
%}
function draw_circle(center, radius, options)
theta = 0.0 : 0.001 : 2.0 * pi;
x = center(1) + radius * cos(theta);
y = center(2) + radius * sin(theta);
plot(x, y, options)
end
clc
clear
close all
for pointCount = [1, 2, 3, 10, 100]
x = randi([-100, 100], pointCount, 1);
y = randi([-100, 100], pointCount, 1);
[center, radius] = min_cover_circle(x, y, pointCount);
figure('color', 'w')
draw_circle(center, radius, 'k--')
hold on
plot(x, y, 'r+')
axis equal off
end
三.LUA代码
--[[-----------------------
Function: get_sign
Description: 求实数x的符号
Input: 实数x
Output: 实数x的符号
Author: Marc Pony(marc_pony@163.com)
--]]------------------------
function get_sign(x)
if math.abs(x) < 1.0e-8 then
return 0
else
if x < 0.0 then
return -1
else
return 1
end
end
end
--[[-----------------------
Function: get_distance_square
Description: 求平面两点之间距离的平方
Input: 平面两点a，b
Output: 平面两点之间距离的平方
Author: Marc Pony(marc_pony@163.com)
--]]------------------------
function get_distance_square(a, b)
return (a[1] - b[1])^2 + (a[2] - b[2])^2
end
--[[-----------------------
Function: get_circle_center
Description: 求三角形外接圆的圆心
Input: 平面三个点a，b，c
Output: 三角形外接圆的圆心center
Author: Marc Pony(marc_pony@163.com)
--]]------------------------
function get_circle_center(a, b, c)
local center = {}
local a1 = b[1] - a[1]
local b1 = b[2] - a[2]
local c1 = 0.5 * (a1 * a1 + b1 * b1)
local a2 = c[1] - a[1]
local b2 = c[2] - a[2]
local c2 = 0.5 * (a2 * a2 + b2 * b2)
local d = a1 * b2 - a2 * b1
center[1] = a[1] + (c1 * b2 - c2 * b1) / d
center[2] = a[2] + (a1 * c2 - a2 * c1) / d
return center
end
--[[-----------------------
Function: get_random_array
Description: 将数组随机排列
Input: 数组array，数组起始索引startIndex，结束索引endIndex
Output: 随机排列后数组randomArray
Author: Marc Pony(marc_pony@163.com)
--]]------------------------
function get_random_array(array, startIndex, endIndex)
local i, index, length
local randomArray = {}
math.randomseed(tostring(os.time()) : reverse() : sub(1, 7)) --设置时间种子
length = endIndex - startIndex + 1
for i = 1, endIndex - startIndex + 1 do
index = math.floor(startIndex + (endIndex - startIndex + 1 - i) * math.random() + 0.5) --执行round功能
randomArray[length] = array[index]
array[index] = array[length]
length = length - 1
end
return randomArray
end
--[[-----------------------
Function: min_cover_circle
Description: 求平面pointCount个点的最小覆盖圆
Input: 平面pointCount个点的坐标(x,y)，点个数pointCount
Output: 平面pointCount个点的最小覆盖圆圆心center，半径radius
Author: Marc Pony(marc_pony@163.com)
--]]------------------------
function min_cover_circle(x, y, pointCount)
local i, j, k, radius, radiusSquare
local center = {}
local p_i = {}
local p_j = {}
local p_k = {}
local xx = {}
local yy = {}
local array = {}
local randomArray = {}
--将坐标随机排列
for i = 1, pointCount do
array[i] = i
xx[i] = x[i]
yy[i] = y[i]
end
randomArray = get_random_array(array, 1, pointCount)
for i = 1, pointCount do
x[i] = xx[randomArray[i]]
y[i] = yy[randomArray[i]]
end
center[1] = x[1]
center[2] = y[1]
radiusSquare = 0.0
for i = 2, pointCount do
p_i[1] = x[i]
p_i[2] = y[i]
if get_sign(get_distance_square(p_i, center) - radiusSquare) > 0 then
center[1] = p_i[1]
center[2] = p_i[2]
radiusSquare = 0.0
for j = 1, i - 1 do
p_j[1] = x[j]
p_j[2] = y[j]
if get_sign(get_distance_square(p_j, center) - radiusSquare) > 0 then
center[1] = 0.5 * (p_i[1] + p_j[1])
center[2] = 0.5 * (p_i[2] + p_j[2])
radiusSquare = get_distance_square(p_j, center)
for k = 1, j - 1 do
p_k[1] = x[k]
p_k[2] = y[k]
if get_sign(get_distance_square(p_k, center) - radiusSquare) > 0 then
center = get_circle_center(p_i, p_j, p_k)
radiusSquare = get_distance_square(p_i, center)
end
end
end
end
end
end
radius = math.sqrt(radiusSquare)
return center, radius
end
-- Test random array
local array = {1, 2, 3, 4, 5, 6, 7, 8, 9, 10}
local startIndex = 1
local endIndex = 10
local randomArray = get_random_array(array, startIndex, endIndex)
-- Test min cover circle
local x = {-63, 35, 35, 17, -74, 38, -97, -78, 91, -49}
local y = {91, 2, -12, -51, 74, -22, -82, -34, 56, 66}
local n = 10
local radius
local center = {}
center, radius = min_cover_circle(x,y, n)
print("min cover circle's center : (", center[1], ", ", center[2], ")")
print("min cover circle's radius : ", radius)
四.C++代码
#include <cstdio>
#include <iostream>
#include <cstring>
#include <cmath>
#include <cstdlib>
#include <algorithm>
#define MAX_POINT_COUNT 100
using namespace std;
struct point
{
double x, y;
};
/*************************************************
Function: get_sign
Description: 求实数x的符号
Input: 实数x
Output: 无
Return: 实数x的符号(0, -1, 1)
Author: Marc Pony(marc_pony@163.com)
*************************************************/
int get_sign(double x)
{
if (fabs(x) < 1.0e-8)
{
return 0;
}
else
{
return (x < 0) ? -1 : 1;
}
}
/*************************************************
Function: get_distance_square
Description: 求平面两点之间距离的平方
Input: 平面两点a，b
Output: 无
Return: 平面两点之间距离的平方
Author: Marc Pony(marc_pony@163.com)
*************************************************/
double get_distance_square(const point a, const point b)
{
return (a.x - b.x) * (a.x - b.x) + (a.y - b.y) * (a.y - b.y);
}
/*************************************************
Function: get_circle_center
Description: 求三角形外接圆的圆心
Input: 平面三个点a，b，c
Output: 无
Return: 三角形外接圆的圆心center
Author: Marc Pony(marc_pony@163.com)
*************************************************/
point get_circle_center(const point a, const point b, const point c)
{
point center;
double a1 = b.x - a.x;
double b1 = b.y - a.y;
double c1 = 0.5 * (a1 * a1 + b1 * b1);
double a2 = c.x - a.x;
double b2 = c.y - a.y;
double c2 = 0.5 * (a2 * a2 + b2 * b2);
double d = a1 * b2 - a2 * b1;
center.x = a.x + (c1 * b2 - c2 * b1) / d;
center.y = a.y + (a1 * c2 - a2 * c1) / d;
return center;
}
/*************************************************
Function: min_cover_circle
Description: 求平面pointCount个点的最小覆盖圆
Input: 平面pointCount个点的数组首元素地址p，点个数pointCount
Output: 平面pointCount个点的最小覆盖圆圆心center，半径radius
Return: 无
Author: Marc Pony(marc_pony@163.com)
*************************************************/
void min_cover_circle(point *p, int pointCount, point ¢er, double &radius)//找最小覆盖圆(这里没有用全局变量p[], 因为是为了封装一个函数便于调用)
{
double radiusSquare;
random_shuffle(p, p + pointCount);//随机函数,使用了之后使程序更快点,也可以不用
center = p[0];
radiusSquare = 0.0;
for (int i = 1; i < pointCount; i++)
{
if (get_sign(get_distance_square(p[i], center) - radiusSquare) > 0)//如果p[i]在当前圆的外面, 那么以当前点为圆心开始找
{
center = p[i];//圆心为当前点
radiusSquare = 0.0;//这时候这个圆只包括他自己.所以半径为0
for (int j = 0; j < i; j++)//找它之前的所有点
{
if (get_sign(get_distance_square(p[j], center) - radiusSquare) > 0)//如果之前的点有不满足的, 那么就是以这两点为直径的圆
{
center.x = 0.5 * (p[i].x + p[j].x);
center.y = 0.5 * (p[i].y + p[j].y);
radiusSquare = get_distance_square(p[j], center);
for (int k = 0; k < j; k++)
{
if (get_sign(get_distance_square(p[k], center) - radiusSquare) > 0)//找新作出来的圆之前的点是否还有不满足的, 如果不满足一定就是三个点都在圆上了
{
center = get_circle_center(p[i], p[j], p[k]);
radiusSquare = get_distance_square(p[i], center);
}
}
}
}
}
}
radius = sqrt(radiusSquare);
}
int main()
{
#define POINT_COUNT 10
point p[MAX_POINT_COUNT], center;
double radius;
double x[POINT_COUNT] = { -63, 35, 35, 17, -74, 38, -97, -78, 91, -49 };
double y[POINT_COUNT] = { 91, 2, -12, -51, 74, -22, -82, -34, 56, 66 };
for (int i = 0; i < POINT_COUNT; i++)
{
p[i].x = x[i];
p[i].y = y[i];
}
min_cover_circle(p, POINT_COUNT, center, radius);
printf("min cover circle's center: (%.4f, %.4f)\n", center.x, center.y);
printf("min cover circle's radius: %.4f\n", radius);
getchar();
return 0;
}
2531
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
