当前位置：首页 > news >正文

Matlab 求不规则图形的内切圆和外接圆函数

news 来源：原创 2024/5/2 6:08:57

1、全代码：

%getZhongxin
function [zhongxin1,zhongxin2,smallR,bigR]=getZhongxin(varargin)
%% 此函数用于计算已知边界的不规则图形的最大内切圆和最小外接圆
% 输出：
%     zhongxin1 最大内切圆圆心
%     zhongxin2 最小外接圆圆心
%     smallR    最大内切圆半径
%     bigR      最小外接圆半径
% 输入：
%     bianjie   不规则图形的边界，2 x n
%     cuX       包裹不规则边界的矩形 x 坐标
%     cuY       包裹不规则边界的矩形 y 坐标
% 调用：
%     [zhongxin1,zhongxin2,smallR,bigR]=getZhongxin(bianjie,X2,Y2);
%     [zhongxin1,zhongxin2,smallR,bigR]=getZhongxin(bianjie);
%-------------------------------------------------------------------
    %%%%    Authors:   Bill O'Hanlon
    %%%%    EMAIL:     ohanlon@qq.com
    %%%%    DATE:      24-08-2020
%% 计算出 cuX,cuY，有输入则跳过
bianjie=varargin{1};
XArray=bianjie(1,:);
YArray=bianjie(2,:);
Xmax1=max(XArray);%未取整
Xmin1=min(XArray);
Ymax1=max(YArray);
Ymin1=min(YArray);

Xmax=(floor(Xmax1*2)+1)*0.5;
Ymax=(floor(Ymax1*2)+1)*0.5;
Xmin=(floor(Xmin1*2))*0.5;
Ymin=(floor(Ymin1*2))*0.5;%这些已经取整了。
Xmax=ceil(Xmax);
Ymax=ceil(Ymax);
Xmin=floor(Xmin);
Ymin=floor(Ymin);
if nargin==3
    cuX=varargin{2}; %注意，这里是{, (会是元胞
    cuY=varargin{3};
elseif nargin==1
    Y=Ymin:0.0083333333:Ymax; %分辨率0.0083333333°  1°=120
    X=Xmin:0.0083333333:Xmax; %求逻辑矩阵用到
    a=size(X,2);b=size(Y,2);
    Y2=repmat(Y',1,a);
    cuX=repmat(X,b,1);
    cuY=flipud(Y2);                %这个纬度得上下翻转一下。
else
    disp('参数不合要求！');
    return;
end
XArray=bianjie(1,:);
YArray=bianjie(2,:);
[a,b]=size(cuX);
zhongdian=[mean(XArray);mean(YArray)];
R=(max(bianjie(1,:))-min(bianjie(1,:))+max(bianjie(2,:))-min(bianjie(2,:)))/6;%当内切圆不对，可以改这个数
zoom=max(a/20,b/20);%最后只保留900左右的点
smallX=imblizoom(cuX,1/zoom);
smallY=imblizoom(cuY,1/zoom);
sX=smallX(:);
sY=smallY(:);
px=[sX';sY'];
n=size(px,2);
sR=zeros(n,1);bR=zeros(n,1);
m=size(bianjie,2);
dis=zeros(m,1);
for i=1:n
    xy=px(:,i);
    det=bianjie-xy;
    for j=1:m
        dis(j)=norm(det(:,j));
    end
    sR(i)=min(dis);
    bR(i)=max(dis);
end
smallR=max(sR);%最大内切圆
bigR=min(bR);  %最小外接圆
zhongxin1=[sX(sR==smallR);sY(sR==smallR)];
zhongxin2=[sX(bR==bigR);sY(bR==bigR)];
for i=1:n
    smallR=max(sR);
    zhongxin1=[sX(sR==smallR);sY(sR==smallR)];
    if norm(zhongxin1-zhongdian)<=R
        break;
    else
        sR(sR==smallR)=0;
        continue;
    end
end
%% 画图
left1=zhongxin1-smallR;  %画内切圆左下角矩形坐标
left2=zhongxin2-bigR;    %画外接圆左下角矩形坐标
figure;
bigcircle=rectangle('Position',[left2',bigR*2,bigR*2],'Curvature',[1,1]);%画圆
hold on;
smallcircle=rectangle('Position',[left1',smallR*2,smallR*2],'Curvature',[1,1]);%画圆
hold on;
plot(bianjie(1,:),bianjie(2,:),zhongxin1(1),zhongxin1(2),'r*',zhongxin2(1),zhongxin2(2),'b*');
bigcircle.FaceColor='r';smallcircle.FaceColor='w';
hold on
detx=Xmax-Xmin; dety=Ymax-Ymin;
rectangle('Position',[Xmin,Ymin,detx,dety]);
axis equal;
axis([Xmin Xmax Ymin Ymax]);
end

依赖函数：

https://blog.csdn.net/Gou_Hailong/article/details/108206521

2、调用

俩效果一样。

[zhongxin1,zhongxin2,smallR,bigR]=getZhongxin(bianjie,X2,Y2);
[zhongxin1,zhongxin2,smallR,bigR]=getZhongxin(bianjie);

在这里插入图片描述
注：这是北京市的shp 搞出的不规则图形。

补于2021-3-3

有的朋友对这个代码的原理感兴趣，有的朋友调试过程遇到了bug…bug的话我没办法解决，毕竟这个代码是根据我当时独特的需求写的，像排除bug自己跑起来的话只好自己逐步调试了。关于原理，我回顾了下，补充如下：

算法原理：
先将边界稀疏化（分辨率降低些，提高计算效率），再计算出来不规则边界的最小外包矩形。逐个取外包矩形中的点和边界上的每一个点计算一下距离，这些距离中最大的距离作为备选最小外接圆半径存到数组bR中，这些距离中最小的距离作为备选最大内切圆半径存到数组sR中；将外包矩形中所有的点计算完毕之后，取bR数组中最小的值作为最小外接圆半径，取sR数组中最大的值作为最大内切圆半径。

注：由于时间较久，不知描述是否符合之前的本意，如有错误欢迎批评指正！