问题四十五：怎么画ray tracing图形中的blending and joining surface

当两个曲面在空间相交时，怎么连接相交处？

比如，两个相互垂直的柱面在空间相交如下：

怎么将黑色圈内的相交处画得平滑些？也就是将相交处包起来，添加一个类似水管的四通接口。

45.1 数学推导

45.2 看C++代码实现

----------------------------------------------quartic_blend_cylinder.h ------------------------------------------

quartic_blend_cylinder.h

#ifndef QUARTIC_BLEND_CYLINDER_H
#define QUARTIC_BLEND_CYLINDER_H#include <hitable.h>
#include "material.h"
#include "log.h"class quartic_blend_cylinder : public hitable
{public:quartic_blend_cylinder() {}quartic_blend_cylinder(vec3 cen1, float a1, float b1, float lh1, vec3 cen2, float a2, float b2, float lh2, float a3, float b3, material *m) {center1 = cen1;intercept_x1 = a1;intercept_z1 = b1;length_half_y1 = lh1;//这个参数无实际用途，只是为了切割看其内部结构center2 = cen2;intercept_y2 = a2;intercept_z2 = b2;length_half_x2 = lh2; //这个参数无实际用途，只是为了切割看其内部结构intercept_s1 = a3;intercept_s2 = b3;ma = m;}virtual bool hit(const ray& r, float tmin, float tmax, hit_record& rec) const;vec3 center1;float intercept_x1, intercept_z1, length_half_y1;vec3 center2;float intercept_y2, intercept_z2, length_half_x2, intercept_s1, intercept_s2;material *ma;
};#endif // QUARTIC_BLEND_CYLINDER_H

----------------------------------------------quartic_blend_cylinder.cpp ------------------------------------------

quartic_blend_cylinder.cpp

#include "quartic_blend_cylinder.h"#include <iostream>
using namespace std;bool quartic_blend_cylinder::hit(const ray& r, float t_min, float t_max, hit_record& rec) const {
#if QUARTIC_BLEND_CYLINDER_LOG == 1std::cout << "-------------quartic_blend_cylinder::hit----------------" << endl;
#endif // QUARTIC_BLEND_CYLINDER_LOGfloat a1_square = intercept_x1*intercept_x1;float b1_square = intercept_z1*intercept_z1;float a1b1_square_a3 = a1_square*b1_square+intercept_s1;float a2_square = intercept_y2*intercept_y2;float b2_square = intercept_z2*intercept_z2;float a2b2_square_b3 = a2_square*b2_square+intercept_s2;float a3_square = intercept_s1*intercept_s1;float b3_square = intercept_s2*intercept_s2;float a3b3_sqare = a3_square*b3_square;float xd = r.direction().x();float yd = r.direction().y();float zd = r.direction().z();float xoc1 = r.origin().x() - center1.x();float zoc1 = r.origin().z() - center1.z();float yoc2 = r.origin().y() - center2.y();float zoc2 = r.origin().z() - center2.z();float A1 = b1_square*xd*xd + a1_square*zd*zd;float B1 = 2*b1_square*xd*xoc1 + 2*a1_square*zd*zoc1;float C1 = b1_square*xoc1*xoc1 + a1_square*zoc1*zoc1;float A2 = A1*A1*b3_square;float B2 = 2*A1*B1*b3_square;float C2 = (B1*B1 + 2*A1*C1 - 2*a1b1_square_a3*A1)*b3_square;float D2 = (2*B1*C1 - 2*a1b1_square_a3*B1)*b3_square;float E2 = (C1*C1 - 2*a1b1_square_a3*C1 + a1b1_square_a3*a1b1_square_a3)*b3_square;float A3 = b2_square*yd*yd + a2_square*zd*zd;float B3 = 2*b2_square*yd*yoc2 + 2*a2_square*zd*zoc2;float C3 = b2_square*yoc2*yoc2 + a2_square*zoc2*zoc2;float A4 = A3*A3*a3_square;float B4 = 2*A3*B3*a3_square;float C4 = (B3*B3 + 2*A3*C3 - 2*a2b2_square_b3*A3)*a3_square;float D4 = (2*B3*C3 - 2*a2b2_square_b3*B3)*a3_square;float E4 = (C3*C3 - 2*a2b2_square_b3*C3 + a2b2_square_b3*a2b2_square_b3)*a3_square - a3b3_sqare;float roots[5];roots_quartic_equation2(A2+A4, B2+B4, C2+C4, D2+D4, E2+E4, roots);float temp;if (roots[0] > 0.0001) {for (int i=1; i<int(roots[0]); i++) {for (int j=i+1; j<int(roots[0])+1; j++) {if (roots[i] > roots[j]) {temp = roots[i];roots[i] = roots[j];roots[j] = temp;}}}//对所有实根从小到大排序for (int k=1; k<int(roots[0])+1; k++) {if (roots[k] < t_max && roots[k] > t_min) {rec.t = roots[k];rec.p = r.point_at_parameter(rec.t);vec3 pc1 = rec.p - center1;vec3 pc2 = rec.p - center2;if (((pc1.y()>= -length_half_y1) && (pc1.y() <= length_half_y1)) && ((pc2.x() >= -length_half_x2) && (pc2.x() <= length_half_x2))) {float nx = 2*b3_square*(b1_square*pc1.x()*pc1.x() + a1_square*pc1.z()*pc1.z() - a1b1_square_a3)*2*b1_square*pc1.x();float ny = 2*a3_square*(b2_square*pc2.y()*pc2.y() + a2_square*pc2.z()*pc2.z() - a2b2_square_b3)*2*b2_square*pc2.y();float nz = 2*b3_square*(b1_square*pc1.x()*pc1.x() + a1_square*pc1.z()*pc1.z() - a1b1_square_a3)*2*a1_square*pc1.z()+ 2*a3_square*(b2_square*pc2.y()*pc2.y() + a2_square*pc2.z()*pc2.z() - a2b2_square_b3)*2*a2_square*pc2.z();rec.normal = unit_vector(vec3(nx, ny, nz));if(dot(r.direction(), rec.normal) > 0) {rec.normal = - rec.normal;}rec.mat_ptr = ma;rec.u = -1.0;rec.v = -1.0;return true;}}}}return false;
}

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

main.cpp

        hitable *list[3];list[0] = new quadratic_cylinder_all(vec3(0, 3, 0), 1, 1, 1.00001, 3,new lambertian(vec3(0.0, 1.0, 0.0)), vec3(0, 1, 0), 0);list[1] = new quadratic_cylinder_all(vec3(0, 3, 0), 1, 1, 1.00001, 3,new lambertian(vec3(0.0, 1.0, 0.0)), vec3(1, 0, 0), 0);list[2] = new quartic_blend_cylinder(vec3(0, 3, 0), 1, 1, 3, vec3(0, 3, 0), 1, 1, 3, 0.5, 0.5, new lambertian(vec3(0.0, 1.0, 0.0)));hitable *world = new hitable_list(list,3);vec3 lookfrom(-1.5, 4, 5);vec3 lookat(0, 3, 0);float dist_to_focus = (lookfrom - lookat).length();float aperture = 0.0;camera cam(lookfrom, lookat, vec3(0,1,0), 80, float(nx)/float(ny), aperture, 0.7*dist_to_focus);

输出图片如下：

换个角度看看：

vec3 lookfrom(0.0001, 8, 0);

vec3 lookat(0, 3, 0);

vec3 lookfrom(5, 3, 0);

vec3 lookat(0, 3, 0);

45.3 blending surface的结构图

vec3 lookfrom(0, 3, 5);

vec3 lookat(0, 3, 0);

vec3 lookfrom(0.0001, 8, 0);

vec3 lookat(0, 3, 0);

vec3 lookfrom(5, 3, 0);

vec3 lookat(0, 3, 0);

vec3 lookfrom(0, 6, 3);

vec3 lookat(0, 3, 0);

vec3 lookfrom(3, 6, 3);

vec3 lookat(0, 3, 0);