任意两个多边形（非凸）相交的面积，目前最全的原理讲解及代码

六级0分

已于 2023-07-19 14:31:33 修改

阅读量1.4k

点赞数 4

分类专栏：数字图像处理文章标签：算法数据结构

于 2023-07-14 17:35:19 首次发布

本文链接：https://blog.csdn.net/qq_43021748/article/details/131717230

版权

数字图像处理专栏收录该内容

3 篇文章 0 订阅

订阅专栏

任意两个多边形（非凸）相交的面积，目前最全的原理讲解及代码

前言
一、算法设计
二、代码实现
三、代码的适用范围以及其他情况的分析
四、算法的性能
五、 Bug Fixed!
五、ToDo

前言

项目中遇到一个求两个任意多边形重叠面积的问题，乍一想感觉问题太简单了；然后就出现了下图的尴尬情况：
在这里插入图片描述

这玩意不好求啊，网上搜索找到一个方案
https://blog.csdn.net/lingyunxianhe/article/details/104948684/
采用了Sutherland-Hodgeman-Polygon-Clipping Algorihtm（自行了解），主要解决了凸多变形重叠区域的contour的寻找，算法的思想就是沿着一个凸多边形的边裁切，直到将所有的凸多边形的边全部裁切完成。但是面对非凸多边形而且有多个闭合的相交区域的情况的时候，一条边的裁切可能会将另一条边裁掉，因此算法不适用。

在这里插入图片描述
又到了造轮子的时候了，人麻了

一、算法设计

整体的算法流程比较简单：

轮廓contour的交点求取
一个轮廓的点与另外一个轮廓的位置关系（包含、不包含、在边界）
轮廓点的正确选取（divide and conquer），包括交点和各自包含在对方轮廓内部的点

1. 轮廓contour的交点求取

轮廓的交点就意味着两个线段的交点，废话少说上图
在这里插入图片描述

上图中可以看到两个线段向量由四个端点来表示，交点的情况大致分为三类：

交点都在两条向量的内部
有交点在向量的顶端
没有交点（没画出来）

如果你知道叉乘的话很容易就能想到使用叉乘就能很方便的解决一些向量上的问题，下面我们以线段vec1作为我们的基准向量，使用叉乘来进行计算：
在这里插入图片描述

点A到交点E的向量长度为L1，向量vec1的长度为L2，求出L1即可确定E的位置。
图中定义了两个新的向量chrod（弦长），s1 = chord1×vec1为2倍的三角形ABC的面积，s2 = vec1×chord2为2倍的三角形ABD的面积，s3 = vec(AE)×vec2为二倍三角形ACD的面积。

观察图像可以知道s3 / (s1 + s2) = L1 / L2 = ratio，我们同时求出L1和ratio，现在就可以得到E的坐标和ratio了，ratio在后面会起到排序的作用，因为一条边上可能会有多个交点。我们只取这个方法求得的ratio = 0~1的点，代表着有实交点，其他情况为虚交点。
另外还有两个线段平行的的情况在，自行判断。
见函数lineIntersect

2. 一个轮廓的点与另外一个轮廓的位置关系

这个在网络上有很多种方式，我在这里使用的是辐角原理（argument theorem），这是复变函数上沿着固定虚数进行封闭曲线进行线积分，只关心其实数部分而推导出来的结果，经典的自动控制领域常用来判断传递函数的稳定性。

建议先去了解一下辐角原理再来看这一节。简单介绍一下原理，下面我绘制了三种情况：
在这里插入图片描述
由左到右分别是点p在轮廓上，点p在轮廓内，点p早轮廓外，θn为p到轮廓上的点的向量角度的变化值

p在轮廓上（仅限于上图在平坦的contour位置）：θ1+θ2+θ3 = 180 * N+（N为正整数）
p在轮廓内：θ1+θ2+θ3 = 360 * N+
p在轮廓外：θ1+θ2+θ3 + θ4 + θ5= 0

比较难处理的点是在轮廓上的点，上面例子中的180只是点p在平坦的contour上，如果点p正好在contour的角点上，那么根据辐角原理，整个变化的角度就是这个角点张成角度的大小，所以可以是0~360°中的任意角度（假设整个contour只环绕了一圈，而不是argument theorem允许的任意圈）。

但是————既然我说出来了，那这个问题肯定是被我克服了，哈哈哈哈或或，代码中使用了一个(bool)vertises_on 来决定是否关心角点上的环绕问题，细节请查看函数argument

3. 轮廓点的正确选取（divide and conquer）

这段的实现稍微复杂，一开始也没有想到使用分治的思想，是在写程序的时候想到的，所以程序还是要自己手写的，所不定就激发灵感了。

（1）思考围成两轮廓相交部分的边的构成方式

想到两个歪七扭八的复杂轮廓相交，交出来亲妈都不认识的图形来就很麻烦，因此我们需要找到一种可以研究这个问题本质的方法论——微分。
在这里插入图片描述
a. 以上面右图微分的方式，将所有微分单元（红线部分）的顺时针的线积分相加会得到原始轮廓（黑线部分）
$\oint_c d(c) = \sum_{i = 0}^n \oint_{c_i}d(c_i)$ 这意味这再复杂的相交部分的情况都是由多个微分单元组成的，了解这个情况可以使我们确定研究简单的相交情况是和而复杂的相交的情况是一致的。

b. 在a.的基础上我们接着再向前跨出一步，以微分（对应着连续） 的视角去分析一块相交的部分，此时我们将两个轮廓分别定义为contour和base。
在这里插入图片描述
我将上图中围成相交轮廓的部分成了四部分

A, … , B：contour 穿入base————> 穿出base（省略号代表存在省略的点但是非交点）
B, … , C：base 穿入contour————> 穿出contour
C, … , D：contour 穿入base————> 穿出base
D, … , A：base 穿入contour————> 穿出contour
其中第1、3和第2、4是同一种情况
其实现在就已经可以发现第1、2种情况存在对偶，这正是使用分治策略的端倪

（2）如何选取起始点以及为什么使用分治的策略

现在我们已经了解了相交部分的边的构成情况，我们只要求出第1中情况的边的点那么第3种边也是适用的，又由于1、2情况对偶，那么使contour作为base，而base作为contour也是OK的，那么我们使用分治思想先处理contour上的情况1、3，再处理base上的情况1、3，即可将所有的情况都包含。

我们在选择起点时我们就选择当前边上的第一个交点点位（前面ratio就在这个时候排上了用场，用于确定点位的顺序，使得我们在循环中正确选取第i个位置，详见代码），并添加进入vector。

这里进入分支，

如果满足contour 穿入base：

a.如果该交点并非该直线上的最后一个交点那么我们在结果vector中继续压入下一个交点，并添加标记点（0，0）。

b.如果为最后一个交点，检查下一个contour在base内部的点，并将其添加进入结果vector中，并添加标记点（0，0）。

如果满足contour 穿出base：

do nothing，因为分治的思想将其转换成base穿出contour

（3）如何计算面积

问：你回忆一下我们的任务是什么
答：算面积啊！
问：那我们做了个什么
答：把相交的点位和在内部的点位全放到vector中了，中间还“乱”放了一些（0，0）… …
问怒：八嘎八嘎，这有个毛用?! 良心大大的坏了。

答：别急，确实没什么毛用，但是能算面积用。
想想之前的contour也是存在vector中的，轻易就能算出面积，但是现在给出来的是什么？
：？
自答：是大风车！
瓜摊老板：你是故意找茬是吧，你能不能算吧。
手机死活对不上焦了，看来iphone15要出了手机有bug了。

瓜摊老板：你这是画了个西瓜，你故意找… …
华强打断施法：我这保熟，你那全是生瓜蛋子【看吸铁石】【小刀劈瓜】【ntm皮我瓜是吧】【萨日朗】【小摩托】【嘟噜嘟噜】

瓜摊老板不会算躺下了，不知道你行不行[滑稽]

每个三角加起来不就是罗，但是这些个三角最外面的边连起来的contour一定要是顺时针的，现在（0，0）想一想是不是只是为了构造三角形的一个点。

瓜摊老板坐起来：你这是原点（0， 0）在contour内部的情况，那要是在外面怎么办呢
华强：建议你回原位躺着

提一嘴green公式，有没有想到点什么， OK解决最后一个问题：
$\int_c Pdx + Qdy = \int\int(\frac{ \partial Q }{ \partial x } - \frac{ \partial P }{ \partial y })dxdy$
我们令P = -y, Q = x，得到
$\frac{1}{2}\int_c -ydx + xdy = \int\int dxdy= S$
这是什么意思，面积等于绕着轮廓的一个线积分，观察这个部分：-ydx + xdy，看出来这是个什么玩意了吗，是个叉乘(x, y)×(dx, dy)，这种情况对应的是上图“大风车”中的虚线边和一条实线边(dx, dy)。
都是三角嘛我们选取两条边都是虚线边可以嘛，当然可以，是等效的，但是积分的路线就不是contour了。我们就得到了离散的算面积的方式：叉乘
$\sum_{i=1}^{n} (x_i, y_i)×(x_{i +1}, y_{i+1})，其中x_{n+1}= x_0，y_{n+1}= y_0$
既然等效了 green公式和原点的位置无关系，那我们的也没啥关系， contour随意在原点的什么位置。

恭喜读完了

二、代码实现

文章写的时候只是将原理部分进行了推导，但是代码涉及到了一些边缘性的问题，并没有做充分的测试，如果有bug的话还请下面评论。
代码中使用到了OpenCV中的cv::Point2f的定义，如果不想安装opencv 的话我提供了一个自己定义的point2f，如下：

typedef struct point2f_
{
    float x;
    float y;
    point2f_ operator+(point2f_ ob)
    {
        point2f_ p;
        p.x = this->x + ob.x;
        p.y = this->y + ob.y;
        return p;
    };
    point2f_ operator-(point2f_ ob)
    {
        point2f_ p;

        p.x = this->x - ob.x;
        p.y = this->y - ob.y;
        return p;
    };      
    point2f_ operator*(float multi)
    {
        point2f_ p;

        p.x = this->x * multi;
        p.y = this->y * multi;
        return p;
    };   
    point2f_(float x, float y)
    {
        // std::cout << x << std::endl;
        this->x = x;
        this->y = y;
        // return * this;
    };
    point2f_(){this->x = 0.0, this->y = 0.0;};    
    ~point2f_(){};
}point2f;

头文件

/// class.h
#include <iostream>
#include <opencv2/opencv.hpp>

#define R2A 180 / 3.1415927
//point and the base contour border start point index
typedef struct _IntsecP_
{
    cv::Point2f p;
    int idx;
    float ratio;
}IntsecP;
/**
 * @param surround test if point are surrounded by contours(return even) or not(return 0), 
 *                  additionally return odd number if the point is on the border of the counter
 * @param angle     angle changes proccess
 * @param angle_d   angle changes step
 */
typedef struct _Argument
{
    int surround; 
    std::vector<float> angle;
    std::vector<float> angle_d;
}Argument;
/// @brief find out the intersaction point of the two vector
/// @author Cillian Young
/// @param p00 
/// @param p01 
/// @param p10 
/// @param p11 
/// @param intsc 
/// @param parallel I- whether detect parallel line
/// @return if true then there is intersacted point
int lineIntersect(cv::Point2f p00, cv::Point2f p01,
                  cv::Point2f p10, cv::Point2f p11, 
                  cv::Point2f& intsc, float *ratio, bool parallel = 1);

/// @brief Test if point are surrounded by contours(returns even surround) or not(returns 0 surround).
///        Additionally returns odd surround number if the point is on the border of the counter
/// @author Cillian Young
/// @param contours I- input contour 
/// @param testP I- point to test 
/// @param jump_overlap I- jump overlap points
/// @param vertices_on I- used to show srrounddingstimes that are non-integer-multiply of 180 
/// @return contour half-surround times; if vertises_on == true, mutiply by 1000 and last three decimal number should be zero, unless there is(are) point(s) on the vertises!
/// @note the situation that points on the vertises or boders have not been handled perfectly but satisfy your expectation
///        (you'd better deal with it outside for boundary result, if you understand the argument theorem),
///       they are transition point so that the angle accumulate is not the intege multiple of 180(degree system)
Argument argument(std::vector<cv::Point2f> contours, cv::Point2f testP, bool jump_overlap = 1, bool vertices_on = 0);

/// @brief 
/// @param contour 
/// @param p 
/// @return in(1),on(0),out(-1)
int inContour(std::vector<cv::Point2f> contour, cv::Point2f p);

/// @brief a class to process contours
class contourProcess
{
private:
    /* data */
    std::vector<cv::Point2f> ct;
public:
    contourProcess(/* args */);
    contourProcess(std::vector<cv::Point2f> contours);
    std::vector<cv::Point2f>& getContour();

    Argument argumentTest(cv::Point2f testP, 
                        bool jump_overlap = 1, 
                        bool vertices_on = 0);    
    /**
    * @brief find each point of the test contour whether in the base contour 
    * @param contour 
    * @param jump_overlap 
    * @param vertices_on 
    * @return 
    */
    std::vector<int> containRelation(std::vector<cv::Point2f> contour, 
                        bool inv = 0,
                        bool jump_overlap = 1, 
                        bool vertices_on = 0);

    /**               
    * @brief every contour edge intersect point with the base contour
    * @author Cillian Young
    * @param contour I- input contour
    * @return return value will be the same size with contour
    */
    std::vector<std::vector<IntsecP>> edgeIntersection(std::vector<cv::Point2f> contour, bool inv = 0);

    /**
    * @brief find a kind of special contour to describe two arbitrary contour's intersection part(s)
    * @author Cillian Young
    * @param contour 
    * @return a kind of special contour that looks like radiation, 
    *         but in terms of area calculation it's the same with ordinary contour
    */
    std::vector<cv::Point2f> findContourIoU(std::vector<cv::Point2f> contour);

    /**
     * calculate contour's area
     */
    float contourArea();
    // ok for less complicate situtation 
    // std::vector<cv::Point2f> findContourIoU2(std::vector<cv::Point2f> contour);
    void print(std::string text);
    ~contourProcess();
};

实现

// class.cpp
#include "class.h"
static float epsilon = 10e-5;

contourProcess::contourProcess(/* args */)
{
}


contourProcess::contourProcess(std::vector<cv::Point2f> contour)
{
    ct = contour;
}

std::vector<cv::Point2f>& contourProcess::getContour()
{
    return this->ct;
}
void contourProcess::print(std::string text)
{
    std::cout << text << std::endl;
    for(int i = 0; i < ct.size(); i++)
    {
        std::cout << ct[i] << std::endl;
    }
}



contourProcess::~contourProcess()
{
}






/// @brief Test if point are surrounded by contours(returns even surround) or not(returns 0 surround).
///        Additionally returns odd surround number if the point is on the border of the counter
/// @author Cillian Young
/// @param contours I- input contour 
/// @param testP I- point to test 
/// @param jump_overlap I- jump overlap points
/// @param vertices_on I- used to show srrounddingstimes that are non-integer-multiply of 180 
/// @return contour half-surround times; if vertises_on == true, mutiply by 1000 and last three decimal number should be zero, unless there is(are) point(s) on the vertises!
/// @note the situation that points on the vertises or boders have not been handled perfectly but satisfy your expectation
///        (you'd better deal with it outside for boundery result, if you understand the argument theorem),
///       they are transition point so that the angle accumulate is not the intege multiple of 180(degree system)
Argument argument(std::vector<cv::Point2f> contours, cv::Point2f testP, bool jump_overlap, bool vertices_on)
{// ret Argument should not been set as reference or it will be released automatically
    Argument ag;
    bool corner_flag = 0;
    // float epsilon = 10e-5;
    int size = contours.size();
    float accumulator = 0.0;
    int rewind_f = 0;
    for(int i = 0; i < size; i++)
    {
        cv::Point2f diff1 = contours[(i) % size] - testP;
        int depth = 0; // 
        while(1 && jump_overlap)
        {
            if(contours[(i + size- depth) % size] == contours[(i + 2 + depth) % size])
            {
                depth++;
                rewind_f = 1;
                // std::cout << "-------" << std::endl;
            }
            else
                break;
        }
        i += depth;     
        cv::Point2f diff2 = contours[(i + 1) % size] - testP;
        if(fabs(diff2.x) < epsilon && fabs(diff2.y) < epsilon)
        {
            // printf("jump %d\n", i + 1);
            i++;    // jump the vertises that overlap
            diff2 = diff1;; // set the angle change as 0
        }
        float tmp = (atan2(diff2.y, diff2.x) - atan2(diff1.y, diff1.x)) * R2A;
        if(fabs(fabs(tmp) - 180.0) < epsilon)
            tmp = 0;    // transition point in mathmatic from 180 -> -180,
        // printf("[atant1, atan2, sub] = %f, %f, %f\n", atan2(diff1.y, diff1.x) * R2A, atan2(diff2.y, diff2.x) * R2A, tmp);
        // if(fabs(tmp) > 135.0 + 10)        
        if(fabs(tmp) > 180) // yjy::need to test
        {
            if(tmp > 0)
                tmp = tmp - 360.0;
            else
                tmp = tmp + 360.0;
        }
        // std::cout << "tmp : " << tmp << std::endl;
        if(rewind_f)
        {
            if(tmp > 0) // 正常的变化不会超过180，为了适应opencv的格式，假设三个相邻点重复时按照角度变化不超过360进行累加
            {
                tmp -= 180.0;
                accumulator += 180.0;
                
                ag.angle.push_back(accumulator);
                ag.angle_d.push_back(180.0);
            }
            else
            {
                tmp += 180.0;
                accumulator -= 180.0;
                ag.angle.push_back(accumulator);
                ag.angle_d.push_back(-180.0);
            }
            rewind_f = 0;
        }
        accumulator += tmp;
        ag.angle.push_back(accumulator);
        ag.angle_d.push_back(tmp);
    }
    // printf("accumulator: %f\n", accumulator);
    int angle_ac = 0;
    if(accumulator > 0)
    {
        angle_ac = accumulator - static_cast<int>(accumulator) < 0.5 ? static_cast<int>(accumulator) : static_cast<int>(accumulator) + 1;
    }
    else{
        angle_ac = static_cast<int>(accumulator) - accumulator < 0.5 ? static_cast<int>(accumulator) : static_cast<int>(accumulator) - 1;
    }
    if(vertices_on)
    {
        ag.surround = angle_ac / 180 * 1000 + angle_ac % 180; // contour half-surround times mutiply by 1000
    }
    else{
        ag.surround = angle_ac / 180;
    }
    return ag; 
}


Argument contourProcess::argumentTest(cv::Point2f testP, bool jump_overlap, bool vertices_on)
{
    return argument(ct, testP, jump_overlap, vertices_on);
}
/// @brief find each point of the test contour whether in the base contour 
/// @param contour 
/// @param jump_overlap 
/// @param vertices_on 
/// @return 
std::vector<int> contourProcess::containRelation(std::vector<cv::Point2f> contour, bool inv, bool jump_overlap, bool vertices_on)
{
    int size;
    if(!inv)
        size = contour.size();
    else
        size = ct.size();      
    std::vector<int> containR(size);
    if(!inv)
        for(int idx = 0; idx < size; idx++){
            containR[idx] = this->argumentTest(contour[idx]).surround;
        }
    else
        for(int idx = 0; idx < size; idx++){
            containR[idx] = argument(contour, ct[idx]).surround;
        }
    return containR;  
}

float crossProduct(cv::Point2f p1, cv::Point2f p2)
{
    return p1.x * p2.y - p1.y * p2.x;
}
float dotProduct(cv::Point2f p1, cv::Point2f p2)
{
    return p1.x * p2.x + p1.y * p2.y;

}
/// @brief every contour edge intersect point with the base contour
/// @author Cillian Young
/// @param contour I- input contour
/// @return return value will be the same size with contour
std::vector<std::vector<IntsecP>> contourProcess::edgeIntersection(std::vector<cv::Point2f> contour, bool inv)
{
    std::vector<std::vector<IntsecP>> intsc_on_border;
    int size_base = ct.size();
    int size_contour = contour.size();
    float ratio;
    if(!inv)
    {
        for(int j = 0; j < size_contour; j++)
        {
            std::vector<IntsecP> intsc_v;   
            int cnt = 0;
            for(int i = 0; i < size_base; i++)
            {
                // std::cout << j * size_base + i << std::endl;
                cv::Point2f intsc;
                int if_intsc = lineIntersect(contour[j], contour[(j + 1) % size_contour], ct[i], ct[(i + 1) % size_base], intsc, &ratio);
                if(if_intsc > 0)
                {
                    // std::cout << "ratio:       " << ratio << std::endl;
                    for(auto iter = intsc_v.begin(); iter != intsc_v.end(); iter++)
                    {
                        if(iter->ratio > ratio)
                        {
                            intsc_v.insert(iter--, {intsc, i, ratio});
                            break;
                        }
                    }
                    if(0 == cnt++)
                        intsc_v.push_back({intsc, i, ratio});                
                }
            }

            intsc_on_border.push_back(intsc_v);
        }        
    }
    else
    {
        for(int i = 0; i < size_base; i++)
        {
            std::vector<IntsecP> intsc_v;   
            for(int j = 0; j < size_contour; j++)
            {
                // std::cout << j * size_base + i << std::endl;
                cv::Point2f intsc;
                int if_intsc = lineIntersect(ct[i], ct[(i + 1) % size_base], contour[j], contour[(j + 1) % size_contour], intsc, &ratio);
                if(if_intsc > 0)
                {
                    // std::cout << "ratio:       " << ratio << std::endl;
                    intsc_v.push_back({intsc, i, ratio});
                }
            }
            intsc_on_border.push_back(intsc_v);
        }  
    }
    return intsc_on_border;
}

float contourProcess::contourArea()
{
	float area0 = 0.f;
	for (int i = 0 ; i < ct.size() ; i++ ) 
    {
		int j = (i+1)%ct.size();
		//this is also a more generical green theorem with dynamic d_x and d_y, integral
		area0 += ct[i].x * ct[j].y;
		area0 -= ct[i].y * ct[j].x;
	}
	return 0.5f * fabs(area0);
}


/// @brief 
/// @param contour 
/// @param p 
/// @return in(1),on(0),out(-1)
int inContour(std::vector<cv::Point2f> contour, cv::Point2f p)
{
    int srnd = argument(contour, p).surround;
    if(srnd > 0 && srnd % 2 == 0)   // make sure all correct base contour points has been pushed
    {
        return 1;
    }
    else if(srnd > 0 && srnd %2 != 0)
    {
        return 0;
    }
    else if(srnd == 0)
    {
        return -1;
    }
    std::cout << "inContour error ..." << std::endl;
    abort();
}
///
void findEdgeOnSelf(std::vector<cv::Point2f>& contour, 
                            std::vector<cv::Point2f>& base_contour, 
                            std::vector<std::vector<IntsecP>>& intsc_v,
                            std::vector<cv::Point2f>& IoU)
{
    float lambda = 0.05;
    int con_size = contour.size();
    for(int i = 0; i < con_size; i++)
    {
        int intsc_num = intsc_v[i].size();
        if(intsc_num)
        {
            for(int j = 0; j < intsc_num; j++)
            {
                // IoU.push_back(intsc_v[i][j].p);
                cv::Point2f extend_p = intsc_v[i][j].p + contour[(i + 1) % con_size] * lambda - contour[i] * lambda;
                if(intsc_v[i][j].ratio < 1.0 - epsilon && intsc_v[i][j].ratio > epsilon) // shrink the range
                {// intsc on border of contour

                    int in = inContour(base_contour, extend_p);
                    std::cout << "in: " << in << std::endl;
                    if(in >= 0)
                    {// contour extend in
                        IoU.push_back(intsc_v[i][j].p);
                        if(j < intsc_num - 1) // find ct points
                        {
                            IoU.push_back(intsc_v[i][j + 1].p);
                            IoU.push_back({0.0, 0.0}); // end signal
                        }
                        else
                        {// shift to next edge
                            int contour_idx = (i + 1) % con_size;
                            if(inContour(base_contour, contour[contour_idx]) >= 0)
                            {
                                IoU.push_back(contour[contour_idx]);
                            }
                        }
                    }
                    else
                    {
                        if(0 == j) // edge's last point on itself
                        {
                            IoU.push_back(intsc_v[i][j].p);
                            IoU.push_back({0.0, 0.0}); // end signal
                        }
                    }
                }
                else
                {// treat as vertises, count on other code to collect, or processed later
                }
            }
        }
        else
        {
            if(inContour(base_contour, contour[(i + 1) % con_size]) >= 0)
            {
                IoU.push_back(contour[(i + 1) % con_size]);
            }
        }
    }
    IoU.push_back(IoU[0]); // make the special iou contour addable
    IoU.push_back({0.0, 0.0});
}
std::vector<cv::Point2f> contourProcess::findContourIoU(std::vector<cv::Point2f> contour)
{    
    std::vector<cv::Point2f> IoU;
    std::vector<std::vector<IntsecP>> intsc_v = edgeIntersection(contour);
    std::vector<std::vector<IntsecP>> intsc_inv = edgeIntersection(contour, true);
    findEdgeOnSelf(contour, ct, intsc_v, IoU);
    findEdgeOnSelf(ct, contour, intsc_inv, IoU);
    return IoU;
}
/// @brief find out the intersaction point of the two vector
/// @author Cillian Young
/// @param p00 
/// @param p01 
/// @param p10 
/// @param p11 
/// @param intsc 
/// @param ratio O- location on the vector <p01 - p00>
/// @param parallel I- whether detect parallel line
/// @return if true then there is intersected point
int lineIntersect(cv::Point2f p00, cv::Point2f p01, // ct , contour
                  cv::Point2f p10, cv::Point2f p11, 
                  cv::Point2f& intsc, float* ratio, bool parallel)
{
    /*chord1
    * 
    * ^(p10)    ^(p01)
    * |        /
    * |       /
    * |      /
    * |     / border
    * |    /
    * |   /
    * |  /
    * (p00)———————>(p11) chord2
    * above vector <c1_1 - c1_0> and <c0_1 - c0_0> will intersect with each other, 
    * code below here will figure out diferent situations and find cross point based on the vectors chord1, boder, chord2*/
    cv::Point2f chord1 = p10 - p00;
    cv::Point2f border = p01 - p00;
    cv::Point2f chord2 = p11 - p00;
    int cross1 = crossProduct(chord1, border);
    int cross2 = crossProduct(border, chord2);
    float epsilon_ = -1.0e-5;
    float epsilon_judge = 1.0e-4;
    if(fabs(cross1) < epsilon_judge && fabs(cross2) < epsilon_judge)
    {
        if(fabs(cross1 + cross2) < epsilon_judge && parallel)  //on the border
        {
            int flag = 0;
            if(dotProduct((p10 - p00), (p10 - p01)) < -epsilon_judge) // tight
            {
                intsc = p10;
                *ratio = ((p10 - p00).x + (p10 - p00).y) / ((p01 - p00).x + (p01 - p00).y);
                flag++;
            }
            if(dotProduct((p11 - p00), (p11 - p01)) < -epsilon_judge) //tight
            {
                intsc = p11;
                *ratio = ((p11 - p00).x + (p11 - p00).y) / ((p01 - p00).x + (p01 - p00).y);
                flag++;
            }  
            assert(*ratio < 1.0);

            return flag;
        }
        else if(fabs(cross1) < epsilon_judge && cross2 > epsilon_judge) // on the vertises
        {
            *ratio = ((p10 - p00).x + (p10 - p00).y) / ((p01 - p00).x + (p01 - p00).y);
            if(*ratio < 1.0 + epsilon_ && *ratio > -epsilon_)             {
                intsc = p10;
                return 1;
            }
        }
        else if(fabs(cross2) < epsilon_judge && cross1 > epsilon_judge) // on the vertises
        {
            *ratio = ((p11 - p00).x + (p11 - p00).y) / ((p01 - p00).x + (p01 - p00).y);
            if(*ratio < 1.0 + epsilon_ && *ratio < -epsilon_)             {
                intsc = p11;
                return 1;
            }
        }  
    }
    else if(cross1 * cross2 > epsilon_judge)
    {
        *ratio = (float)crossProduct(chord1, chord2) / (cross1 + cross2); //dividend is 2*area of ct0, contour0, contour1
        // std::cout << "ratio: " << *ratio << std::endl;
        // extend the range a little bit for endpoint
        if(*ratio > 1.0 + epsilon_ || *ratio < -epsilon_) // outside cross, perhaps set '+' as '-' will reduce the same point to output
        {
            // std::cout << "return ratio: " << *ratio << std::endl;
            return 0;
        }
        intsc = cv::Point2f(*ratio * border.x + p00.x, *ratio * border.y + p00.y);
        return 1;
    }

    return 0;
}

调用

// main.cpp
#include <iostream>
#include <string.h>
#include "class.h"
#include <opencv2/opencv.hpp>
void call_IoU(contourProcess cp, contourProcess cp2)
{
    contourProcess iou = cp.findContourIoU(cp2.getContour()); // implicit conversion
    std::cout << "area : " << iou.contourArea() << std::endl;
    iou.print("iou contour: ");
}
int main()
{      
    // auto start = std::chrono::system_clock::now();
    // std::cout << "start: " << start << std::endl;
    contourProcess cp, cp2;
    cp.getContour().push_back(cv::Point2f(0,0));
    cp.getContour().push_back(cv::Point2f(10, 0));
    cp.getContour().push_back(cv::Point2f(10,9));
    cp.getContour().push_back(cv::Point2f(4,4));
    cp.getContour().push_back(cv::Point2f(9,10));
    cp.getContour().push_back(cv::Point2f(0, 10));

    cp2.getContour().push_back(cv::Point2f(5, 5));    
    cp2.getContour().push_back(cv::Point2f(11, 5));
    // cp2.getContour().push_back(cv::Point2f(3, 9));
    // cp2.getContour().push_back(cv::Point2f(0, 5));
    // cp2.getContour().push_back(cv::Point2f(3, 1));
    cp2.getContour().push_back(cv::Point2f(11, 11));
    // cp2.getContour().push_back(cv::Point2f(5, 11));
    // cp2.getContour().push_back(cv::Point2f(9, 4));
    // cp2.getContour().push_back(cv::Point2f(5,8));
    // cp2.getContour().push_back(cv::Point2f(5,6));


    call_IoU(cp, cp2); 
    // call_argument(cp, cp2);
    // call_intersection(cp, cp2);
    return 0;
}

三、代码的适用范围以及其他情况的分析

代码适用的轮廓为简易表示的轮廓就像opencv中的cv::CHAIN_APPROX_SIMPLE，组成轮廓的点不是紧挨着，用于紧挨着的轮廓效率不高，紧密排布的点可以根据该程序更改。

我这里主要想求的是面积，至于获得的轮廓我并未进行排序，但是他是可排序的，去除（0，0）、重复点、被（0，0）包围的孤立点，然后根据辐角原理进行排序即可。

四、算法的性能

算法的复杂度为O(mn)， m、n分别为两个轮廓的点的个数但是要比SHPC算法多几个常数倍，但他不挑轮廓；
算法复杂度与图像尺寸无相关性，这点很重要，求相交的轮廓有人可能想到用十分简单的方法：轮廓的内部都设置值为1，轮廓一和轮廓二相交的部分就是他们相加之后值为2的部分。小图无伤大雅，大图等着哭把，内存和耗时疯涨。

五、 Bug Fixed!

一共有两个bug，这两个bug仅在此处做一次更改上面的代码还是有bug 的状态，使用的话请自己对照下面代码改正。

一个是因为图省事，直接复制的相同的代码片段，结果更改的时候忘了要同时更改，还是不要图省事了，拷贝代码确实不太优雅，关键是有隐藏的bug（当然下面的代码还是重复的复制，太晚了就没有改）。

// fixed version
/// @brief every contour edge intersect point with the base contour
/// @author Cillian Young
/// @param contour I- input contour
/// @return return value will be the same size with contour
std::vector<std::vector<IntsecP> > contourProcess::edgeIntersection(std::vector<point2f> contour, bool inv)
{// drawback: just simply make another copy of the code for "inv" 
    std::vector<std::vector<IntsecP> > intsc_on_border;
    int size_base = ct.size();
    int size_contour = contour.size();
    float ratio;
    if(!inv)
    {
        for(int j = 0; j < size_contour; j++)
        {
            std::vector<IntsecP> intsc_v;   
            int cnt = 0;
            for(int i = 0; i < size_base; i++)
            {
                // std::cout << j * size_base + i << std::endl;
                point2f intsc;
                int if_intsc = lineIntersect(contour[j], contour[(j + 1) % size_contour], ct[i], ct[(i + 1) % size_base], intsc, &ratio);
                if(if_intsc > 0)
                {
                    // std::cout << "ratio:       " << ratio << std::endl;
                    for(auto iter = intsc_v.begin(); iter != intsc_v.end(); iter++)
                    {
                        if(iter->ratio > ratio)
                        {
                            intsc_v.insert(iter--,(IntsecP){intsc, i, ratio});
                            break;
                        }
                    }
                    if(0 == cnt++)
                        intsc_v.push_back((IntsecP){intsc, i, ratio});                
                }
            }

            intsc_on_border.push_back(intsc_v);
        }        
    }
    else
    {
        for(int i = 0; i < size_base; i++)
        {
            std::vector<IntsecP> intsc_v;   
            int cnt = 0;
            for(int j = 0; j < size_contour; j++)
            {
                // std::cout << j * size_base + i << std::endl;
                point2f intsc;
                int if_intsc = lineIntersect(ct[i], ct[(i + 1) % size_base], contour[j], contour[(j + 1) % size_contour], intsc, &ratio);
                if(if_intsc > 0)
                {
                    // std::cout << "ratio:       " << ratio << std::endl;
                    for(auto iter = intsc_v.begin(); iter != intsc_v.end(); iter++)
                    {
                        if(iter->ratio > ratio)
                        {
                            intsc_v.insert(iter--,(IntsecP){intsc, i, ratio});
                            break;
                        }
                    }
                    if(0 == cnt++)
                        intsc_v.push_back((IntsecP){intsc, i, ratio});       
                    // // std::cout << "ratio:       " << ratio << std::endl;
                    // intsc_v.push_back((IntsecP){intsc, i, ratio});
                }
            }
            intsc_on_border.push_back(intsc_v);
        }  
    }
    return intsc_on_border;
}

另一个是编程过程中的遗漏，少添加了一个处点（0,0）作为标记点，关键这个位置还比较隐蔽。人脑还是有局限的，不能考虑到这么多层的东西，或者我的编程技术还有待提高；当然这也从反面体现出来了技巧的重要性，如果技巧到位了debug都方便了，代码会比较清晰

// fixed version
void findEdgeOnSelf(std::vector<point2f>& contour, 
                            std::vector<point2f>& base_contour, 
                            std::vector<std::vector<IntsecP> >& intsc_v,
                            std::vector<point2f>& IoU)
{
    float lambda = 0.05;
    int con_size = contour.size();
    bool flag_point_btwn_intsc = 0; // bug fix for intersected points
    for(int i = 0; i < con_size; i++)
    {
        if(flag_point_btwn_intsc)
        {
            IoU.push_back(point2f(0.0, 0.0)); // end signal
            flag_point_btwn_intsc = 0;
        }
        int intsc_num = intsc_v[i].size();
        if(intsc_num)
        {
            for(int j = 0; j < intsc_num; j++)
            {
                // IoU.push_back(intsc_v[i][j].p);
                point2f extend_p = intsc_v[i][j].p + contour[(i + 1) % con_size] * lambda - contour[i] * lambda;
                
                if(intsc_v[i][j].ratio < 1.0 - epsilon && intsc_v[i][j].ratio > epsilon) // shrink the range
                {// intsc on border of contour

                    int in = inContour(base_contour, extend_p);
                    // std::cout << "in: " << in << std::endl;
                    if(in >= 0)
                    {// contour extend in
                        IoU.push_back(intsc_v[i][j].p);
                        if(j < intsc_num - 1) // find ct points
                        {
                            IoU.push_back(intsc_v[i][j + 1].p);
                            IoU.push_back(point2f(0.0, 0.0)); // end signal
                        }
                        else
                        {// shift to next edge
                            int contour_idx = (i + 1) % con_size;
                            if(inContour(base_contour, contour[contour_idx]) >= 0)
                            {
                                IoU.push_back(contour[contour_idx]); // points that between intersected points on the contour
                                flag_point_btwn_intsc = 1;
                            }
                        }
                    }
                    else
                    {
                        if(0 == j) // edge's last point on itself
                        {
                            IoU.push_back(intsc_v[i][j].p);
                            IoU.push_back(point2f(0.0, 0.0)); // end signal
                        }
                    }
                }
                else
                {// treat as vertises, count on other code to collect, or processed later
                }
            }
        }
        else
        {

            if(inContour(base_contour, contour[(i + 1) % con_size]) >= 0)
            {
                IoU.push_back(contour[(i + 1) % con_size]);
            }
        }
    }
    IoU.push_back(IoU[0]); // make the special iou contour addable
    IoU.push_back(point2f(0.0, 0.0));
}

五、ToDo

本程序生成的放射状的轮廓点中包含一些孤立点：被两个（0， 0）夹住的非（0， 0）点，这些点对轮廓的生成和面积的计算均无影响，可以通过条件判断将其删除，自行尝试更改函数 findEdgeOnSelf 即可。
生成的轮廓点并没有按照顺序排列（项目中并不需要），想进行排序的话调用函数 Argument ag = argument(IoU, cv::Point2f(0, 0)); 得到ag.angle将其限制在范围-180~180按照极角进行排序即可 。

六级0分

关注

4
点赞
踩
12

收藏

觉得还不错? 一键收藏
1
评论
任意两个多边形（非凸）相交的面积，目前最全的原理讲解及代码

项目中遇到一个求两个任意多边形重叠面积的问题，乍一想感觉问题太简单了；然后就出现了下图的尴尬情况：这玩意不好求啊，网上搜索找到一个方案采用了Sutherland-Hodgeman-Polygon-Clipping Algorihtm（自行了解），主要解决了凸多变形重叠区域的contour的寻找，算法的思想就是沿着一个凸多边形的边裁切，直到将所有的凸多边形的边全部裁切完成。但是面对非凸多边形而且有多个闭合的相交区域的情况的时候，一条边的裁切可能会将另一条边裁掉，因此算法不适用。
复制链接

扫一扫