问题五十六：怎么用ray tracing画参数方程表示的曲面（3）—

接“问题五十四”，已经简单学习了bezier曲线曲面，知道了双三次bezier曲面的矩阵表示形式，同时也以此画出了曲面图形。

这一章节主要以对比bezier曲线曲面和b-spline曲线曲面的方式来简单学习b-spline曲线曲面。

56.1 b-spline曲线

56.2 b-spline曲面

56.3 b-spline的C++代码实现

----------------------------------------------parametric_surface.h ------------------------------------------

parametric_surface.h

#ifndef PARAMETRIC_SURFACE_H
#define PARAMETRIC_SURFACE_H#include <hitable.h>
#include "material.h"
#include "log.h"#define TYPE 3
/*
TYPE=1: sphere;
TYPE=2: horn;
TYPE=3: bezier3;
*/
#define CURVE 1
/*
CURVE=1: bezier
CURVE=2: b-spline
*/class parametric_surface : public hitable
{public:parametric_surface() {}
#if ((TYPE == 1) || (TYPE ==2))parametric_surface(vec3 cen, float a, float b, float c, float hy, material *m, float r, float t, float n) : center(cen), intercept_x(a), intercept_y(b), intercept_z(c), height_half_y(hy), ma(m), rho(r), theta(t), num(n) {}
#endif // TYPE
#if TYPE == 3parametric_surface(material *m, float r, float t, float n, vec3 *cp) : ma(m), rho(r), theta(t), num(n) {
/*相比bicubic Bezier surface，bicubic b-spline surface只需要改动如下黄底部分代码*/vec3 ctrl_points[4][4];for (int i=0; i<16; i++) {ctrl_points[(i/4)][(i%4)] = cp[i];}
#if CURVE == 1//bezierfloat matrix_t[4][4] = {{ 1,  0,  0, 0},{-3,  3,  0, 0},{ 3, -6,  3, 0},{-1,  3, -3, 1}};float matrix[4][4] = {{1, -3,  3, -1},{0,  3, -6,  3},{0,  0,  3, -3},{0,  0,  0,  1}};
#endif // CURVE

#if CURVE == 2//b-spline

float matrix_t_6[4][4] = {{ 1, 4, 1,0},

{-3, 0, 3,0},

{ 3,-6, 3, 0},

{-1, 3, -3, 1}};

float matrix_t[4][4], matrix[4][4];

for (int i=0; i<4; i++) {

for (int j=0; j<4; j++) {

matrix_t[i][j] =matrix_t_6[i][j] / 6.0;

matrix[j][i] = matrix_t[i][j];

}

#endif // CURVE

            float points_x[4][4], points_y[4][4], points_z[4][4], points_x_t[4][4], points_y_t[4][4], points_z_t[4][4];get_bezier_matrix_xyz_4_4(ctrl_points, points_x, points_y, points_z);matrix_4_4_multiply_4_4(matrix_t, points_x, points_x_t);matrix_4_4_multiply_4_4(points_x_t, matrix, matrix_c_x);matrix_4_4_multiply_4_4(matrix_t, points_y, points_y_t);matrix_4_4_multiply_4_4(points_y_t, matrix, matrix_c_y);matrix_4_4_multiply_4_4(matrix_t, points_z, points_z_t);matrix_4_4_multiply_4_4(points_z_t, matrix, matrix_c_z);float min_x = points_x[0][0];float max_x = points_x[0][0];float min_y = points_y[0][0];float max_y = points_y[0][0];float min_z = points_z[0][0];float max_z = points_z[0][0];for (int i=0; i<4; i++) {for (int j=0; j<4; j++) {if (min_x > points_x[i][j]) {min_x = points_x[i][j];}if (max_x < points_x[i][j]) {max_x = points_x[i][j];}if (min_y > points_y[i][j]) {min_y = points_y[i][j];}if (max_y < points_y[i][j]) {max_y = points_y[i][j];}if (min_z > points_z[i][j]) {min_z = points_z[i][j];}if (max_z < points_z[i][j]) {max_z = points_z[i][j];}}}bezier3_bl = vec3(min_x, min_y, min_z);bezier3_bh = vec3(max_x, max_y, max_z);}
#endif // TYPEvirtual bool hit(const ray& r, float tmin, float tmax, hit_record& rec) const;
#if ((TYPE == 1) || (TYPE ==2))vec3 center;float intercept_x;float intercept_y;float intercept_z;float height_half_y;
#endif // TYPEmaterial *ma;float rho, theta, num;
#if TYPE == 3float matrix_c_x[4][4], matrix_c_y[4][4], matrix_c_z[4][4];vec3 bezier3_bl, bezier3_bh;
#endif // TYPE
};#endif // PARAMETRIC_SURFACE_H

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

main.cpp

/*这个文件bicubic Bezier surface，bicubic b-spline surface是完全一样的*/

//1vec3 ctrl_points[16] = {vec3(-3, 2,  0), vec3(-2, 3,  0), vec3(2, 3,  0), vec3(3, 2,  0),vec3(-3, 5, -2), vec3(-2, 5, -2), vec3(2, 5, -1), vec3(3, 5, -2),vec3(-3, 5, -4), vec3(-2, 5, -4), vec3(2, 5, -3), vec3(3, 5, -4),vec3(-3, 2, -6), vec3(-2, 3, -6), vec3(2, 3, -4), vec3(3, 2, -6)};hitable *list[1];list[0] = new parametric_surface(new lambertian(vec3(1.0, 0.0, 0.0)), 0.05, 20, 200);hitable *world = new hitable_list(list,1);vec3 lookfrom(10, 10, 10);vec3 lookat(0, 3, -2);float dist_to_focus = (lookfrom - lookat).length();float aperture = 0.0;camera cam(lookfrom, lookat, vec3(0,1,0), 20, float(nx)/float(ny), aperture, 0.7*dist_to_focus);

输出图片bicubic Bezier surface（前），bicubic b-spline surface（后）对比：

#if TYPE == 3

#define UN 16

#define VN 16

#endif // TYPE

//1

vec3 ctrl_points[4][4] = {vec3(-3, 2, 0), vec3(-2, 3, 0), vec3(2, 3, 0), vec3(3, 2, 0)},

{vec3(-3,5, -2), vec3(-2, 5, -2), vec3(2, 5, -1), vec3(3, 5, -2)},

{vec3(-3,5, -4), vec3(-2, 5, -4), vec3(2, 5, -3), vec3(3, 5, -4)},

{vec3(-3,2, -6), vec3(-2, 3, -6), vec3(2, 3, -4), vec3(3, 2, -6)};

//2

vec3 ctrl_points[16] = {vec3(-3, 2, 0), vec3(-2, 3, 0), vec3(2, 3, 0), vec3(3, 2, 0)},

{vec3(-3,5, -2), vec3(-2, 5, -2), vec3(2, 4, -1), vec3(3,3, -2)},

{vec3(-3,5, -4), vec3(-2, 5, -4), vec3(2, 4, -3), vec3(3,3, -4)},

{vec3(-3,2, -6), vec3(-2, 3, -6), vec3(2, 3, -4), vec3(3, 2, -6)};

//3

vec3 ctrl_points[16] = {vec3(-3,2, 0), vec3(-2, 3, 0), vec3(2, 3, 0), vec3(3, 2, 0)},

{vec3(-3,5, -2), vec3(-2, 5, -2), vec3(2, 2, -1), vec3(3,1, -2)},

{vec3(-3,5, -4), vec3(-2, 5, -4), vec3(2, 2, -3), vec3(3,1, -4)},

{vec3(-3,2, -6), vec3(-2, 3, -6), vec3(2, 3, -4), vec3(3, 2, -6)};

//4

vec3ctrl_points[16] = {vec3(-3, 2, 0),vec3(-2, 3, 0), vec3(2, 5, 0), vec3(3, 5, 0)},

{vec3(-3,5, -2), vec3(-2, 5, -2), vec3(2, 2, -1), vec3(3, 1, -2)},

{vec3(-3,5, -4), vec3(-2, 5, -4), vec3(2, 2, -3), vec3(3, 1, -4)},

{vec3(-3,2, -6), vec3(-2, 3, -6), vec3(2, 5, -4), vec3(3,5, -6)};

//5

vec3 ctrl_points[16] = {vec3(-3, 2, 0), vec3(-2, 3, 0), vec3(2, 5, 0), vec3(3, 5, 0)},

{vec3(-3,5, -2), vec3(-2, 5, -2), vec3(2, 0, -1), vec3(3, -1, -2)},

{vec3(-3,5, -4), vec3(-2, 5, -4), vec3(2, 0, -3), vec3(3,-1, -4)},

{vec3(-3,2, -6), vec3(-2, 3, -6), vec3(2, 5, -4), vec3(3, 5, -6)};

//11

vec3 ctrl_points[16] = {vec3(-3,2, 0), vec3(-2, 3, 0), vec3(2, 5, 0), vec3(-6, 3, 0),

vec3(-3, -1,-2), vec3(-2, 5, -2), vec3(2, 0, -1), vec3(-6, 3, -2),

vec3(-3, -1,-4), vec3(-2, 5, -4), vec3(2, 0, -3), vec3(-6, 3, -4),

vec3(-3, 2, -6), vec3(-2, 3, -6), vec3(2,5, -4), vec3(-6, 3, -6)};