C++图形学之线性插值

原文链接

// 图形学线性插值：已知三角形三个顶点的坐标以及需要插入的点P,首先判断P是否在三角形内，然后求解权重。

#include <iostream>
#include <math.h>

const double eps = 1e-6;
using namespace std;

struct  node
{
    double x, y;
    double w;
}p[5];

double dis(node A,node B)
{
    return sqrt(pow((A.x - B.x), 2) + pow((A.y-B.y),2));
}

bool equ(double A, double B)
{
    if (fabs(A - B) < eps) return true;
    return false;
}

//由边长求面积
double area(node A,node B, node C)
{
    double len1 = dis(A, B), len2 = dis(B, C),len3 = dis(A,C);
    double p = (len1 + len2 + len3) / 2;
    return sqrt(p*(p-len1)*(p-len2)*(p-len3));
}

int main(int argc, char *argv[])
{
    /*
    前三行每行包括三个浮点数，分别是x,y,w,表示（x,y）这个坐标的三角顶点权重为w
    第四行包括两个浮点数，分别是x,y,表示需要计算的插值点坐标。  0 <= x,y,w <= 100
    输入：0 0 0
          0 1 1
          1 0 2
          0.25 0.25
    输出：0.75

          0 0 0
          0 1 1
          1 0 2
          2 2
    输出：-1
    */
    while (cin >> p[0].x >> p[0].y >> p[0].w)
    {
        for (int i = 1; i < 3; i++)
        {
            cin >> p[i].x >> p[i].y >> p[i].w;
        }
        cin >> p[3].x >> p[3].y; //第三个点

        for (int i = 3; i >= 0; i--)
        {
            p[i].x -= p[0].x;
            p[i].y -= p[0].y;
        }

        //判断插值点是否在三角形内部
        if (equ(area(p[0],p[1],p[2]),area(p[0],p[1],p[3])+area(p[0],p[2],p[3])+area(p[1],p[2],p[3])) == false)
        {
            cout << "-1" << endl;
            continue;
        }

#if 0
        double B = (p[1].x*p[3].y - p[3].x*p[1].y) / (p[1].x*p[2].y - p[2].x*p[1].y);
        double A;
        if (equ(p[1].x, 0) == false) A = (p[3].x - B*p[2].x) / p[1].x;
        else A = (p[3].y - B*p[2].y) / p[1].y;
        printf("%.2f\n", A*p[1].w + B*p[2].w + (1.0 - A - B)*p[0].w);
#else
        /*
        p[3].x = (1-B-A)*p[0].x + B*p[1].x + A*p[2].x
        p[3].y = (1-B-A)*p[0].y + B*p[1].y + A*p[2].y

        求解
        */

        double B = (p[3].x*p[2].y - p[2].x*p[3].y) / (p[1].x*p[2].y - p[2].x*p[1].y);
        double A;
        if (equ(p[2].x, 0) == false) A = (p[3].x - B*p[1].x) / p[2].x;
        else A = (p[3].y - B*p[1].y) / p[2].y;
        printf("%.2f\n", A*p[2].w + B*p[1].w + (1.0 - A - B)*p[0].w);
#endif

    }

    return 0;
}