一种误差较小的轮廓面积计算算法

1.背景

        基于微分思想的轮廓面积计算方法之一是将多边形轮廓边与X轴会Y轴进行围合，形成一个个梯形，每个梯形的面积有符号，累计求和即得到多边形轮廓的面积。详见博主之前的文章，

记录导致计算轮廓面积出错的一个坑点-CSDN博客文章浏览阅读377次，点赞4次，收藏9次。计算轮廓面积是常见的几何算法话题，获取轮廓面积、计算轮廓法线等场景会涉及到。计算轮廓面积的方法有很多，一种常用的是微积分思路的分段求和办法，即组成轮廓的每条线段与X轴或Y轴进行有向投影，轮廓边线与X轴或Y轴的投影之和即为轮廓的有向面积。其中第0、1、6的投影为正，第2、3、4、5的投影为负，累积即为轮廓面积（绿色填充部分）。如上所示，如果关心面积的正负时，一般将轮廓移到第一或第四象限再求面积。满足上述条件时，逆时针轮廓的面积为正，顺时针轮廓的面积为负。https://blog.csdn.net/baidu_38621657/article/details/140853728

        如果出现轮廓跨Y轴的情况，则需要将轮廓平移到Y轴右侧，此时面积符号为逆正顺负（当然平移到Y轴左侧也可以，如果那样所计算面积的符号为“顺正逆负”）。

        但是这种方法可能会出现精度误差较大，从而导致面积值和符号计算错误的情况，究其原因是因为对轮廓做了偏置处理，而这种处理会带来误差，偏置的距离越大误差可能越大，详细原因见博主此前的文章：《关于float浮点值二进制存储和运算精度损失的话题》。

关于float浮点值二进制存储和运算精度损失的话题_float的二进制存储和输出-CSDN博客文章浏览阅读1k次，点赞8次，收藏18次。浮点值的存储、运算都可能会带来精度损失，了解精度损失背后的机制原因方便我们更好的了解什么情况下会发生精度损失、什么情况下精度损失较大，以及思考怎么避免或减少精度损失。_float的二进制存储和输出https://blog.csdn.net/baidu_38621657/article/details/141027014        需要尽量减少和简化计算过程，保持算法的简便整洁，注意这里的简化不是接受更大的误差的意思，指的是简化逻辑，使用精简准确的公式或算法过程。

对于此种情况（轮廓跨越Y轴）该怎么处理，避免不必要的操作过程，从而减少精度损失呢？

2. 公式推导

博主比较懒， 推导过程就不详细讲述了，有兴趣可以和博主联系或关注博主公众号。

3. 代码

/// <summary>
/// 求有向面积，顺时针为负，逆时针为正
/// </summary>
/// <param name="polygon"></param>
/// <returns></returns>
double GetArea2(const list<Line>& polygon, const Transform& trsW2L)
{
    //  目前只支持边为线段的轮廓
    double dArea = 0;
    for (auto curve : polygon)
    {
        Vector3f pt0, pt1;
        Transform::MultPoint(trsW2L, curve.pt0, pt0);
        Transform::MultPoint(trsW2L, curve.pt1, pt1);

        double areaThis = 0.0;
        if (pt0.X >= 0.f)
        {
            if (pt1.X >= 0.f)
                areaThis = (pt0.X + pt1.X) * (pt1.Y - pt0.Y) * 0.5;
            else
            {
                double s1 = (pt1.Y - pt0.Y) * pt0.X * pt0.X * 0.5 / (pt0.X - pt1.X);
                double s2 = -(pt1.Y - pt0.Y) * pt1.X * pt1.X * 0.5 / (pt0.X - pt1.X);
                areaThis = s1 + s2;
            }
        }
        else
        {
            if (pt1.X >= 0.f)
            {
                double s1 = -(pt1.Y - pt0.Y) * pt0.X * pt0.X * 0.5 / (pt1.X - pt0.X);
                double s2 = (pt1.Y - pt0.Y) * pt1.X * pt1.X * 0.5 / (pt1.X - pt0.X);
                areaThis = s1 + s2;
            }
            else
            {
                areaThis = (pt0.X + pt1.X) * (pt1.Y - pt0.Y) * 0.5;
            }
        }
         dArea += areaThis;
    }

    return dArea;
}

各种象限和跨象限情况的单元测试均通过。

欢迎交流，相互学习，公众号：geometrylib