辐角原理判断点和多边形的关系

\qquad 标题写的多边形，然而在现实情况中多边形只是该问题的一个子集；我们先讲点和连通域的关系，点和多边形的关系自然就包含在其中了。
\qquad 工作中碰到的这个问题的时候第一时间就想到了辐角原理：没有误差（除了边界点，问题直接就是连续的），实现简单（边界条件比其他算法少）。

一、辐角原理

二、适用条件

\qquad 能想到的复杂形状都能判别，多个连通域、连通域嵌套也可以判别。
\qquad 有时间再写，先把代码贴出来😍。

三、算法效果

实验图像：

实验结果：

四、matlab代码实现

\qquad 程序编写会比较简单，但是会用到图像边界的计算的程序；简单写写吧，这个代码对于单连通域是没有问题的，有兴趣的话可以稍微尝试修改，然后就可以得到对于多连域的程序了

%% 辐角原理的在 164 行
clear
close all
clc;
%% 算适用于连通区域为0的区域
%预处理阶段
a = imread('white10.bmp');
a(a<125) = 0;
a(a>=125) = 255;

%处理原图像使之符合要求
a = 255-a;
a = double(a);
b1 = zeros(size(a));
se1 =  [1 1 1;1 1 1;1 1 1];
for ii = 2:size(a,1)-1
    for jj = 2:size(a,2)-1
        val = sum(sum(a(ii - 1:ii+1,jj-1:jj+1).*se1));
        if val == 255 * 9
            b1(ii,jj) = 255;
        end
    end
end
a = b1;

subplot(221)
imshow(a);
title('原图');
a = double(a);
%% 计算边界
b = zeros(size(a));
se =  [1 1 1;1 1 1;1 1 1];
for ii = 2:size(a,1)-1
    for jj = 2:size(a,2)-1
        val = sum(sum(a(ii - 1:ii+1,jj-1:jj+1).*se));
        if val == 255 * 9
            b(ii,jj) = 255;
        end
    end
end
%erode
subplot(222)
imshow(b)
title('腐蚀');
c = a - b;%白色边界
subplot(223)
imshow(c)
title('边界');
subplot(224)
imshow(a)
title('边界连通域');
hold on;
% 顺时针记录连通域边界
x_pos = [];
y_pos = [];
pos_count = 1;
%% 连通域、边界长度
xId = [-1 0 1 1 1 0 -1 -1];
yId = [-1 -1 -1 0 1 1 1 0];
indexMap = zeros(size(c));
markMap  = zeros(size(c));
len = 0;%初始边缘长度
time = 0.0001;%绘制时间
%************搜索模板，搜索起点的值为9,顺时针搜索
%|1   2   3 |
%|8   9   4 |
%|7   6   5 |
%************
for i = 1:size(c,1)*size(c,2)
    i
    if c(i) ~= 0
        x= floor(i/size(c,1)) + 1;
        y= i-(x-1)*size(c,1);
        xStart = x
        yStart = y
        %         x = 68; %设置调试起点,自己设置的
        %         y = 144;
        indexMap(y,x) = 5;%索引起始点
        markMap(y,x)  = 2;%设置起始访问标记值
        plot(x,y,'r*');
        x_pos(pos_count) = x;
        y_pos(pos_count) = y;
        pos_count = pos_count + 1;
        tmp = 1;    %搜索位置起点:1~9
        tmp = tmp +2; %一般设置为1（按道理怎么设置都可以）
        flag = 0;%停止标志
        while 1     %进入该循环之后开始寻找连通域
            if flag
                break;
            end
            if tmp-3 == 0%数组循环,3是回退数量，保证搜索的位置全部包含
                start = tmp;
            else
                start = mod(tmp-3,9);
            end
            for j = start:start + 7 %在当前点位搜索的起始点
                if  j == 9
                    j = 8;
                else
                    jT = mod(j,9);
                end
                xTmp = x + xId(jT);
                yTmp = y + yId(jT);
                %                 c(yTmp,xTmp);%标记点
                %                 if indexMap(yTmp,xTmp) >= 9 %停止条件
                if yTmp == yStart && xTmp ==xStart  %停止条件
                    xEnd = x
                    yEnd = y
                    flag = 1;
                    break;
                end
                if c(yTmp,xTmp) >=125 && markMap(yTmp,xTmp) == 0 %未搜寻且为边界
                    if mod(jT,2) == 1 %在奇数位置进入下一个位置的寻找
                        for z = 1:7
                            if jT + z ==9
                                start2 = 8;
                            else
                                start2 = mod(jT + z,9);
                            end
                            if c(y + yId(start2),x + xId(start2)) == 0
                                break;
                            end
                        end
                    else
                        z = 1;
                    end
                    if jT + z - 1 == 9
                        start1 = 8;
                    else
                        start1 = mod(jT + z-1,9);
                    end
                    xTmp = x + xId(start1);
                    yTmp = y + yId(start1);
                    if xId(start1)*yId(start1)
                        len = len+2^(1/2);
                    else
                        len = len+1;
                    end
                    indexMap(y,x) = indexMap(y,x) + start1;%连通域索引
                    markMap(y,x)  = 1;%标记占用
                    plot(xTmp,yTmp,'r.');   %绘制连通域点
                    pause(time);
                    x = xTmp;
                    y = yTmp;
                    x_pos(pos_count) = x;
                    y_pos(pos_count) = y;
                    pos_count = pos_count + 1;
                    tmp = jT + z -1;
                    break;
                end
            end
        end
        break;
    end
end

figure
plot(x_pos,-y_pos)
% %% 圆形度
% len  = len + sqrt((x-xStart)^2 + (y-yStart)^2);
% s = len^2/4/pi;
% Ratio = area / s
% str = num2str(Ratio);
% handleT = title(['圆形度：',str]);
% handle = get(handleT,'Position');

%% 辐角原理判断点与连通域的关系（包含、不包含、在边界）
margin = 0.1;   %角度差值的裕度
margin2 = 0.0000002;   %判断点是否包含的裕度（除去边界，误差理论上只有计算机的舍入误差），margin2可以比较小，
x_pos = [x_pos,x_pos(1)]; %连通域闭环
y_pos = [y_pos,y_pos(1)];
result = 255 * ones(size(a));
if size(x_pos,2) > 2
    for j = 1:size(a,1)
        for i = 1:size(a,2)
%             i = x_pos(end - 1) % 调试使用
%             j = y_pos(end - 1)
            dx = x_pos(1) - i;
            dy = y_pos(1) - j;
            theta = atan2(dy,dx); % 初始角度
            accumulator = 0; %角度变化累加
            % 该循环计算每个像素点的累加角度积分
            for k = 2:size(x_pos,2) 
                theta_tmp = theta; %记录上一个角度
                dx = x_pos(k) - i;
                dy = y_pos(k) - j;
                flag = 0;
                if dx == 0 && dy == 0
                    flag = 1;
                    break;
                end
                %计算theta
                theta = atan2(dy,dx);
                d_theta = theta - theta_tmp;
                if abs(d_theta)> pi / 2 + margin % 注意连通域上两个连续点之间的角度差值不会大于90°即pi/2（理论值）,除去点在边界的情况最大差值为 pi
                    if d_theta > 0
                        d_theta = d_theta - 2 * pi;
                    else
                        d_theta = d_theta + 2 * pi;
                    end
                end
%                 d_theta
                accumulator = accumulator + d_theta;
            end
            %判断包含关系
            if flag == 0
                if abs(accumulator) < margin2 % 判断为不包含
                    result(j,i) = 128;
                else
                    while abs(accumulator)>pi * 2 - margin2
                        if accumulator > 0
                            accumulator = accumulator - pi * 2;
                        else
                            accumulator = accumulator + pi * 2;
                        end
                    end
                    if abs(accumulator) < margin2 % 判断为包含
                        result(j,i) = 0;
                    else
                        disp('程序编写很大概率出错，并没有达到margin2精度要求！');
                    end
                end
            end
        end
    end
else
    disp('Short of Points!');
end
figure
imshow(uint8(result));title('灰度128判断在连通域外，0判断在边界，255判断在内部','fontsize',16);