问题四十三：对ray tracing圆环图形中的细微问题进行修正

第一步：将上一章节最后一张图放大，看局部

有两处问题：

其一：蓝色圈内，出现多余的红色像素点

其二：黄色圈内，出现多余的白色像素点（也就是原本应该出现的红色像素点没有出现）

 

第二步：去掉其它圆环和圆柱面，只留下红色圆环

放大（16倍）看局部：

由于图片像素少，所以没有出现白色的缝（小图中白色不足一次一个像素点，所以，白色的效果就是和红色叠加产生最后的颜色较浅的红色像素）。另外，小图中看不到大图蓝色圈内的多余红色像素点（太少了，小图中的光线根本没有撞击到）。

第三步：“多余”或“缺少”像素点对应一元四次方程的根

当前的问题是“问题四十”的残留问题。
截取“问题四十：第七步：2”如下：

 

将对应代码修改如下：

bool roots_quartic_equation2(float a, float b, float c, float d, float e, float (&roots)[5]) {
    //the first element is the number of the real roots, and other elements are the real roots.
    //Descartes's Method.
    if (a == 0) {
        float *roots3;
        roots3 = roots_cubic_equation(b, c, d, e);
        for (int i=0; i<int(roots3[0])+1; i++) {
            roots[i] = roots3[i];
        }
        delete [] roots3;
    }
    else {
        float a1 = b/a;
        float b1 = c/a;
        float c1 = d/a;
        float d1 = e/a;
        float p = (-3*a1*a1)/8 + b1;
        float q = (a1*a1*a1)/8 - (a1*b1)/2 + c1;
        float r = (-3*a1*a1*a1*a1)/256 + (a1*a1*b1)/16 - (a1*c1)/4 + d1;
        float sa = 2*p;
        float sb = p*p - 4*r;
        float sc = -q*q;
        float k, m, t;
        float *roots_s = roots_cubic_equation(1.0, sa, sb, sc);
        float temp;
        for (int i=1; i<int(roots_s[0]); i++) {
            for (int j=i+1; j<int(roots_s[0])+1; j++) {
                if (roots_s[i] > roots_s[j]) {
                    temp = roots_s[i];
                    roots_s[i] = roots_s[j];
                    roots_s[j] = temp;
                }
            }
        }
        if (roots_s[int(roots_s[0])] > 0) {
            k = sqrt(roots_s[int(roots_s[0])]);
            delete [] roots_s;
        }

        else{

           if (fabs(q) < 0.0000001) {

               k = 0.0;

               delete [] roots_s;

           }

           else {

               roots[0] = 0.0;

               delete [] roots_s;

               return true;

           }

        }

        if (fabs(k) == 0.0) {

           if (fabs(q) < 0.0000001) {

                float *roots2;

                roots2 =roots_quadratic_equation(1, -p, r);

                if (roots2[0]== 0.0) {

                    roots[0] =0.0;

                    delete []roots2;

                    returntrue;

                }

                else {

                    m =roots2[1];

                    t =roots2[2];

                }

                delete []roots2;

           }

           else {

               roots[0] = 0.0;

               return true;

           }

        }
        else {
            m = (k*k*k + k*p + q) / (2*k);
            t = (k*k*k + k*p - q) / (2*k);
        }

        float *roots_y1;
        float *roots_y2;
        roots_y1 = roots_quadratic_equation(1.0, k, t);
        roots_y2 = roots_quadratic_equation(1.0, -k, m);
        if (roots_y1[0] != 0.0) {
            for (int i=1; i<roots_y1[0]+1; i++) {
                roots[i] = roots_y1[i] - a1/4;
            }
        }
        if (roots_y2[0] != 0.0) {
            int roots_y1_number = int(roots_y1[0]);
            for (int j=1; j<roots_y2[0]+1; j++) {
                roots[roots_y1_number+j] =roots_y2[j] - a1/4;
            }
        }
        roots[0] = roots_y1[0] + roots_y2[0];
        delete [] roots_y1;
        delete [] roots_y2;
    }
    return true;
}

 

“多根”的问题搞定了，现在再查查看哪个地方导致了“少根”

 

第四步：分析一元三次方程

拦到这样一组数据：

但是，用网页上的一元四次方程计算器求解是有实根的：

这样看来，是这个一元四次方程求错了。

根据之前的分析，应该是对应一元三次方程的求解除了问题。

 

在main()函数单独测试一个一元四次方程的求解：

 

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

main.cpp

    int main(){

        float roots[5];
        if (roots_quartic_equation2(1.4656595, -1.43020415, 4.95183134, -2.24647474, 0.0000393861628, roots)) {
            for (int i=0; i<(roots[0]+1); i++) {
                std::cout << "roots[" << i << "]=" << roots[i] << endl;
            }
        }

 

另外，在一元四次方程求解函数中添加log：

 

对应的一元三次方程系数上面已经标出，现在在一元三次方程求解函数中设置断点，另截图如下：

 

对应的一元三次方程系数上面已经标出，现在在一元三次方程求解函数中设置断点，另截图如下：

对应的一元三次方程有一个实根，是负的，而且是10^-8数量级，尼玛，和计算过程的0是一个数量级哈。

一个很小很小很小非常非常接近0的负数，但也是个负数啊，所以导致原一元四次方程没有实根。

 

网页计算器求解该一元三次方程结果如下：

尼玛，尼玛，什么鬼，什么鬼？？？

 

发现，网页上的正确答案是一个10^-8数量级的正实根，我们求解出来的一个10^-8数量级的负实根，两个实根相差很小很小（10^-8数量级）。

所以，我们有理由猜测是计算精度的问题。

所以，所以，我们将一元三次方程求解过程的数据全部换成“double”类型，函数修改如下：

float* roots_cubic_equation2(float a, float b, float c, float d) {
    //the first element is the number of the real roots, and other elements are the real roots.
    //Shengjin's formula
    float *roots = new float[4];
    if (a == 0) {
        float *roots2;
        roots2 = roots_quadratic_equation(b, c, d);
        for (int i=0; i<int(roots2[0])+1; i++) {
            roots[i] = roots2[i];
        }
        delete [] roots2;
    }
    else {
        double A = double(b*b - 3*a*c);
        double B = double(b*c - 9*a*d);
        double C = double(c*c - 3*b*d);
        double deita = double(B*B - 4*A*C);

        if ((A == B) && (A == 0)) {
            //the three roots are the same
            if (a != 0) {
                roots[1] = -b/(3*a);
            }
            else {
                if (b != 0) {
                    roots[1] = -c/b;
                }
                else {
                    if (c != 0) {
                        roots[1] = -3*d/c;
                    }
                }
            }
            roots[2] = roots[1];
            roots[3] = roots[1];
            roots[0] = 3;
        }
        else if (deita > 0) {
            //only one real root
            double y1 = double(A*b + (3*a)*(-B + sqrt(deita))/2);
            double y2 = double(A*b + (3*a)*(-B - sqrt(deita))/2);
            double pow_y1, pow_y2;
            if (y1 < 0) {
            //for pow(a,b), when b is not int, a should not be negative.
                pow_y1 = - pow(double(-y1), double(1.0/3.0));
            }
            else {
                pow_y1 = pow(double(y1), double(1.0/3.0));
            }
            if (y2 < 0) {
                pow_y2 = - pow(double(-y2), double(1.0/3.0));
            }
            else {
                pow_y2 = pow(double(y2), double(1.0/3.0));
            }
            roots[1] = float((-b - pow_y1 - pow_y2) / (3*a));
            roots[0] = 1;
        }
        else if (deita == 0) {
            //three real roots and two of them are the same
            double K = B/A;
            roots[1] = float(-b/a + K);
            roots[2] = float(-K/2);
            roots[3] = float(-K/2);
            roots[0] = 3;
        }
        else if (deita < 0) {
            //three different real roots
            double theta = acos((2*A*b-3*a*B) / (2*pow(A, 1.5)));
            roots[1] = float((-b - 2*sqrt(A)*cos(theta/3)) / (3*a));
            roots[2] = float((-b + sqrt(A) * (cos(theta/3) + sqrt(3)*sin(theta/3))) / (3*a));
            roots[3] = float((-b + sqrt(A) * (cos(theta/3) - sqrt(3)*sin(theta/3))) / (3*a));
            roots[0] = 3;
        }
    }
    return roots;
}

修改后，运行结果如下：

第五步：验证修改一元三次方程求解函数后图片效果

输出图片如下：

对比之前的图片：

 

放大截图：

对比之前的截图：

 

组合图：

大图~大图~看大图：

哦也~搞定~