问题四十三：对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 formulafloat *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 sameif (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 rootdouble 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 samedouble 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 rootsdouble 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;
}

修改后，运行结果如下：

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

输出图片如下：

对比之前的图片：

放大截图：

对比之前的截图：

组合图：

大图~大图~看大图：

哦也~搞定~