问题三十二：怎么用ray tracing画多边形（polygon, triangle）

画多边形主要分为两步：

1，光线和多边形所在的平面相交，求得交点；

2，判断交点是否在多边形内；

 

32.1 光线和多边形所在的平面相交

分别定义光线、多边形和多边形所在平面。光线和平面的方程如下：（注意：用到法向量的地方一般都将其标准化。尽管有时候没有必要，但是在某些公式或算法中只能使用标准法向量（而不是法向量）。所以，为了避免特别记忆，统一都进行标准化。）

 

光线方程：

多边形：

平面方程：

如上Pn=[A, B, C]为平面的单位法向量。D是原点[0, 0, 0]到平面的距离。

可以先通过多边形任意两条不平行的边对应的向量的叉乘求得平面的法向量；然后代入任意顶点的坐标可以求出常数D。

 

将（C1）展开，然后代入（C2），如下可以求得t.

 

用向量表示是这样的：

通常，我们要求的平面法向量是面向光线的那一条（而不是远离光线的那一条），所以，我们将法向量的方向作如下调整：

32.2 判断交点是否在多边形内

在“32.1”中求得：

多边形的法向量Pn:

判断交点是否在多边形内的步骤如下：

 

第一步：将所有点（顶点和交点）投影到二维UV平面。

（投影是不会改变图的拓扑结构的）

 

根据法向量Pn=[A, B, C]的各个坐标值的绝对值的大小，选出绝对值最大的作为主坐标。比如：Pn=[-5, 3, 1]则X坐标-5的绝对值最大，则X为主坐标；

Pn=[-5, 8, 1]则Y坐标8的绝对值最大，则Y为主坐标；

Pn=[-5, 3, 10]则Z坐标10的绝对值最大，则Z为主坐标；

 

扔掉所有点的三维坐标中的主坐标，剩余的两个坐标即为二维坐标中的U和V。比如：

若，Pn=[-5, 3, 1]则X坐标-5的绝对值最大，则X为主坐标

所有点的X坐标都扔掉，剩余Y、Z坐标：

交点：Ri=[Xi, Yi, Zi]--------------------------[Yi, Zi]即为[Ui, Vi]

顶点：Gn=[Xn, Yn, Zn]----------------------[Yn, Zn]即为[Un, Vn]

 

第二步：将UV平面上交点[Ui, Vi]平移到UV的原点[0, 0]。

则：

交点：Ri=[Xi, Yi, Zi]--------------------------[Yi, Zi]即为[Ui,Vi]-------------------[0, 0]

顶点：Gn=[Xn, Yn, Zn]----------------------[Yn, Zn]即为[Un,Vn]----------------[Un-Ui, Vn-Vi]

 

第三步：计算从原点出发、沿+U轴的射线与多边形交点的个数。

 

交点个数为奇数，则Ri在原多边形内；

交点个数为偶数，则Ri在原多边形外。

先说明一点，像F点这样的顶点在坐标轴上，我们认为是在坐标轴上方的（即EF和坐标轴相交，FA不和坐标轴相交）

具体怎么数多边形和+U轴的交点个数呢？

假设多边形的顶点顺序为ABCDEF（即六条边为：AB、BC、CD、DE、EF、FA。注意：第一个顶点A会被使用两次）。一条边一条边依次判断是否和+U轴相交，每条边的判断依据如下：

 

1，边要与+U轴相交，前提是必须与U轴（含-U和+U）相交：也就是前后两个顶点的V坐标必须是异号的。

2，前后两顶点异号的边与U轴相交有如下8种情况：

32.3 看C++代码实现

----------------------------------------------polygon.h------------------------------------------

polygon.h

#ifndef POLYGON_H
#define POLYGON_H

#include <hitable.h>


class polygon : public hitable
{
    public:
        polygon() {}
        polygon(vec3 *v, int n, material *m) : vertexes(v), number(n), ma(m) {}
        virtual bool hit(const ray& r, float tmin, float tmax, hit_record& rec) const;
        vec3 *vertexes;/*多边形顶点数组/列表首地址*/
        int number;/*多边形顶点的个数*/
        material *ma;/*多边形的材质*/
};

#endif // POLYGON_H

----------------------------------------------polygon.cpp------------------------------------------

polygon.cpp

#include "polygon.h"
#include "vec2.h"
#include "log.h"

#include <iostream>
#include <fstream>

using namespace std;

bool in_polygon_test(vec2 *vertexes_uv, int number) {
/*根据多边形在UV平面的投影与+U轴交点的个数来判断交点是否在多边形内*/
        int sh, nsh;
        int nc = 0;
        if(vertexes_uv[0].v() < 0) { sh = -1;}
        else { sh = 1;}
        for(int j=0; j<number; j++) {
            if(vertexes_uv[j+1].v() < 0) { nsh = -1;}
            else { nsh = 1;}
            if(sh != nsh) {//前后两个顶点的纵坐标V异号
                if((vertexes_uv[j].u() > 0) && (vertexes_uv[j+1].u() >0)) { nc++;}
/*前后两个顶点的横坐标U都大于零，这条边必定过+U轴*/
                else
                    if((vertexes_uv[j].u() > 0) || (vertexes_uv[j+1].u() >0)) {
/*前后两个顶点的横坐标异号，这条边过+U轴的条件是：Uj - Vj /k > 0*/
                        if(vertexes_uv[j].u() - (vertexes_uv[j].v())*(vertexes_uv[j+1].u()-vertexes_uv[j].u())/(vertexes_uv[j+1].v()-vertexes_uv[j].v()) > 0) { nc++;}
                    }
            }
            sh = nsh;
        }
        if((nc)%(2)) {return true;}//交点数为奇数，原交点在多边形内，返回true
        else {return false;}
}

bool polygon::hit(const ray& r, float t_min, float t_max, hit_record& rec) const {
        vec3 poly_n;
        for(int i=0; i<number-2; i++) {
/*以连续三个顶点构成的两个平行的向量的叉乘得到平面的法向量，并对法向量进行标准化*/
            poly_n = unit_vector(cross((vertexes[i]-vertexes[i+1]), (vertexes[i+1]-vertexes[i+2])));
//determine the normal of the plane
            if(!vector_equ(poly_n, vec3(0,0,0))) {//法向量不为0，即所选的两个向量不平行
                break;
            }
        }
        float poly_d = -(dot(poly_n, vertexes[0]));
/*determine the distance from the origin to the plane。随便找一个顶点（这里找的是第一个顶点），代入方程求得参数D（即为原点到平面的距离）*/
        float vd = dot(poly_n, r.direction());
        float v0 = -(dot(poly_n, r.origin()) + poly_d);

        if(vd == 0) {//the ray is parallel to the polygon plane
            return false;
        }
        else {
            float t = v0/vd;//determine t and intersection pi
            vec3 pi = r.point_at_parameter(t);//get the coordinates of pi

            float temp = abs(poly_n.x());
            int i = 1;
            if(temp <= abs(poly_n.y())) {
                temp = abs(poly_n.y());
                i++;
            }
            if(temp <= abs(poly_n.z())) {
                i++;
            }
            /*find the dominant coordinate, X, Y, or Z?
            i=1: means that X is the dominant coordinate;
            i=2: means that Y is the dominant coordinate;
            i=3: means that Z is the dominant coordinate;
            */

/*the following “switch…case…” throw the dorminant coordinate of 3-d vector, then we get 2-d vector in uv-plane*/
            vec2 vertexes_uv[number+1];
//声明二维向量数组。（需要创建一个vec2的类，类似与vec3，不赘述）
            switch (i) {
            case 1://throw X coordinates
                for(int i=0; i<number; i++) {
                    vertexes_uv[i] = vec2(vertexes[i].y(),vertexes[i].z());
                }
                vertexes_uv[number] = vec2(pi.y(),pi.z());
                break;
            case 2: //throw Y coordinates
                for(int i=0; i<number; i++) {
                    vertexes_uv[i] = vec2(vertexes[i].x(),vertexes[i].z());
                }
                vertexes_uv[number] = vec2(pi.x(),pi.z());
                break;
           case 3: //throw Z coordinates
                for(int i=0; i<number; i++) {
                    vertexes_uv[i] = vec2(vertexes[i].x(),vertexes[i].y());
                }
                vertexes_uv[number] = vec2(pi.x(),pi.y());
                break;
            }

            /*this “for” loop moves intersection uv-coordinate to origin.
            so all the vertexes substract intersection uv-coordinate.*/
            for(int i=0; i<number; i++) {
                vertexes_uv[i] = vertexes_uv[i] - vertexes_uv[number];
            }
            vertexes_uv[number] = vertexes_uv[0];
/*the first vertex will be used twice, so set the first vertex to the last position of the array, so that we get the whole vertexes loop*/
            if(in_polygon_test(vertexes_uv,number)) {
/*check if the intersection locates inside the polygon or not. If the intersection is inside the polygon, the current ray hits the polygon, otherwise, misses the polygon. Hits: store related information about the polygon and return true, note that the poly_n should be the unit vector. Misses: return false.*/
                if (t < t_max && t > t_min) {
                    rec.t = t;
                    rec.p = r.point_at_parameter(rec.t);
                    rec.normal = poly_n;
                    rec.mat_ptr = ma;
                    return true;
                }
                 return false;
            }
            return false;
        }
}

----------------------------------------------main.cpp------------------------------------------

main.cpp

main()函数中

//triangle1, the light pink metal one
        vec3 vertexes3_1[3];
        vertexes3_1[0] = vec3(3,-0.5,-2);
        vertexes3_1[1] = vec3(-1,-0.5,-4);
        vertexes3_1[2] = vec3(0,5,-3);
//triangle2, the green lambertian one
        vec3 vertexes3_2[3];
        vertexes3_2[0] = vec3(-1,1.5,-1);
        vertexes3_2[1] = vec3(-3,1.5,-1);
        vertexes3_2[2] = vec3(-2,3.5,-1);
//pentagon
        vec3 vertexes5[5];
        vertexes5[0] = vec3(3.3000,1.3000,0.0000);
        vertexes5[1] = vec3(2.3489,0.6090,0.0000);
        vertexes5[2] = vec3(2.7122,-0.5090,0.0000);
        vertexes2[3] = vec3(3.8878,-0.5090,0.0000);
        vertexes5[4] = vec3(4.2511,0.6090,0.0000);

        hitable *list[5];
        list[0] = new polygon(vertexes3_2, 3, new lambertian(vec3(0.3, 0.8, 0.0)));
        list[1] = new sphere(vec3(0,-100.5,-1), 100, new lambertian(vec3(0.8, 0.8, 0.0)));
        list[2] = new sphere(vec3(-2,0.5,-1), 1, new lambertian(vec3(0.5, 0.7, 0.6)));
        list[3] = new polygon(vertexes3_1, 3, new metal(vec3(0.8, 0.6, 0.5), 0.0));
        list[4] = new polygon(vertexes5, 5, new lambertian(vec3(1.0, 0.0, 0.0)));
        hitable *world = new hitable_list(list,5);

        vec3 lookfrom(0,0,10);
        vec3 lookat(0,1,-1);
        float dist_to_focus = (lookfrom - lookat).length();
        float aperture = 0.0;
        camera cam(lookfrom, lookat, vec3(0,1,0), 30, float(nx)/float(ny), aperture, 0.7*dist_to_focus);

两个漫射球，一个大的（半径为100），一个小的（半径为1）；

两个三角形，一个漫射材料，一个镜面材料；

一个五边形，漫射材料。

 

输出图片如下：

意外发现：多边形顶点的定义顺序决定着最终的多边形形状。

比如：

当五边形的顶点定义顺序为ABCDE时（即之前的图片），赋值为：

//pentagon

        vec3 vertexes5[5];

        vertexes5[0] =vec3(3.3000,1.3000,0.0000);

        vertexes5[1] =vec3(2.3489,0.6090,0.0000);

        vertexes5[2] =vec3(2.7122,-0.5090,0.0000);

        vertexes2[3] =vec3(3.8878,-0.5090,0.0000);

        vertexes5[4] =vec3(4.2511,0.6090,0.0000);

若将五边形的顶点定义顺序改为ACEBD，会是什么图片呢？？赋值为：

//pentacle

        vec3 vertexes5[5];

        vertexes5[0] = vec3(3.3000,1.3000,0.0000);

        vertexes5[1] =vec3(2.7122,-0.5090,0.0000);

        vertexes5[2] =vec3(4.2511,0.6090,0.0000);

        vertexes5[3] =vec3(2.3489,0.6090,0.0000);

        vertexes5[4] =vec3(3.8878,-0.5090,0.0000);

输出图片是这样的：

咦~咦，五角星怎么是空心的呢？

根据我们的算法，从原点沿着+U轴的射线与多边形的交点个数为奇数时，表示原交点在多边形内，否则在多边形外。

如上图，从O沿着+U轴交点（蓝点）是两个，不是奇数，所以说O点在多边形外。上图中心的蓝色五边形都是这种情况，所以，这个蓝色五边形不属于多边形ACEBD。

这个结果和实际情况有异，我们希望这个蓝色五边形属于多边形ACEBD。

怎么修改我们的算法呢？？？

以多边形ACEBD为例，该多边形有五条边：AC、CE、EB、BD、DA。与+U轴相交的是CE和DA。

对于CE，是从-V轴到+V轴，所以NC=NC-1；

对于DA，是从-V轴到+V轴，所以NC=NC-1；

所以NC的最终结果是-2，不是0，所以O点在多边形ACEBD内。

 

顺便提一下（与当前问题无关）：什么情况下NC=NC+1呢？将上面两条边反过来EC、AD，则是从+V轴到-V轴，所以NC=NC+1。

 

怎么修改代码呢？？？

 

----------------------------------------------polygon.cpp------------------------------------------

polygon.cpp

bool in_polygon_test2(vec2 *vertexes_uv, int number) {
/*根据winding number是否为0来判断原交点是否在多边形内*/
        int sh, nsh;
        int nc = 0;
        if(vertexes_uv[0].v() < 0) { sh = -1;}
        else { sh = 1;}
        for(int j=0; j<number; j++) {
            if(vertexes_uv[j+1].v() < 0) { nsh = -1;}
            else { nsh = 1;}
            if(sh != nsh) {
/*前后两个顶点的纵坐标V异号。所以，大的V值为正，小的V值为负。前顶点的比后顶点的大，则是从+V到-V，则NC=NC+1；反之，从-V到+V，NC=NC-1*/
                if((vertexes_uv[j].u() > 0) && (vertexes_uv[j+1].u() >0)) {
                    if(vertexes_uv[j].v() > vertexes_uv[j+1].v()) { nc = nc + 1;}
                    else { nc = nc - 1;}
                }
                else
                    if((vertexes_uv[j].u() > 0) || (vertexes_uv[j+1].u() >0)) {
                        if(vertexes_uv[j].u() - (vertexes_uv[j].v())*(vertexes_uv[j+1].u()-vertexes_uv[j].u())/(vertexes_uv[j+1].v()-vertexes_uv[j].v()) > 0) {
                            if(vertexes_uv[j].v() > vertexes_uv[j+1].v()) { nc = nc + 1;}
                            else { nc = nc - 1;}
                        }
                    }
            }
            sh = nsh;
        }
        if(nc != 0) {return true;}
        else {return false;}
}

采用新的判定算法后，输出图片为：

接下来再在图片中添加一个介质三角形（因为漫射材料、镜面材料的多边形多有了，只差介质材料了）

//triangle1, the light pink metal one
        vec3 vertexes3_1[3];
        vertexes3_1[0] = vec3(3,-0.5,-2);
        vertexes3_1[1] = vec3(-1,-0.5,-4);
        vertexes3_1[2] = vec3(0,5,-3);
//triangle2, the green lambertian one
        vec3 vertexes3_2[3];
        vertexes3_2[0] = vec3(-1,1.5,-1);
        vertexes3_2[1] = vec3(-3,1.5,-1);
        vertexes3_2[2] = vec3(-2,3.5,-1);
//triangle3, the bigger dielectric one
        vec3 vertexes3_4[3];
        vertexes3_4[0] = vec3(4,-0.5,4);
        vertexes3_4[1] = vec3(-4,-0.5,4);
        vertexes3_4[2] = vec3(0,4,4);
//pentacle
        vec3 vertexes5[5];
        vertexes5[0] = vec3(3.3000,1.3000,0.0000);
        vertexes5[1] = vec3(2.7122,-0.5090,0.0000);
        vertexes5[2] = vec3(4.2511,0.6090,0.0000);
        vertexes5[3] = vec3(2.3489,0.6090,0.0000);
        vertexes5[4] = vec3(3.8878,-0.5090,0.0000);

        hitable *list[6];
        list[0] = new polygon(vertexes3_2, 3, new lambertian(vec3(0.3, 0.8, 0.0)));
        list[1] = new sphere(vec3(0,-100.5,-1), 100, new lambertian(vec3(0.8, 0.8, 0.0)));
        list[2] = new sphere(vec3(-2,0.5,-1), 1, new lambertian(vec3(0.5, 0.7, 0.6)));
        list[3] = new polygon(vertexes3_1, 3, new metal(vec3(0.8, 0.6, 0.5), 0.0));
        list[4] = new polygon(vertexes5, 5, new lambertian(vec3(1.0, 0.0, 0.0)));
        list[5] = new polygon(vertexes3_4, 3, new dielectric(1.5));
        hitable *world = new hitable_list(list,6);

        vec3 lookfrom(0,0,10);
        vec3 lookat(0,1,-1);
        float dist_to_focus = (lookfrom - lookat).length();
        float aperture = 0.0;
        camera cam(lookfrom, lookat, vec3(0,1,0), 30, float(nx)/float(ny), aperture, 0.7*dist_to_focus);

加介质三角形之前：

折射系数为1.5：

折射系数为3：

折射系数为4.5：

折射系数为9：

上面加的介质三角形比较大，下面将介质三角形改小一点。

//triangle3, the smaller dielectric one

        vec3vertexes3_3[3];

       vertexes3_3[0] = vec3(1.5,-0.5,4);

       vertexes3_3[1] = vec3(-1.5,-0.5,4);

       vertexes3_3[2] = vec3(0,2,4);

加介质三角形之前（还是这张图）：

折射系数为1.5：

折射系数为3：

折射系数为4.5：

折射系数为9：