c - C 中的光线追踪 - 圆柱体、法向 vector 和旋转

Question

我的第一篇文章在这里。 :)我目前正在用 C 语言为我的学校项目编写光线追踪器。我已经可以显示带有一些灯光效果的球体、三角形和平面。现在我想显示圆柱体(然后是圆锥体，但首先是圆柱体!)。我选择有一个平行于 Y 轴的圆柱体。所以我必须解决:x² + z² = r²。所以从技术上讲，我返回相机和交点之间距离的函数是:

double          n_ray_cylinder(t_ray *ray, t_cylinder *cylinder, t_data *d)
{
    double a;
    double b;
    double c;
    double delta;
    double root;

    a = ray->dir->x * ray->dir->x + ray->dir->z * ray->dir->z;
    b = 2 * ray->dir->x * (ray->ori->x - cylinder->base->x) + 2 * ray->dir->z * (ray->ori->z - cylinder->base->z);
    c = (ray->ori->x - cylinder->base->x) * (ray->ori->x - cylinder->base->x) + (ray->ori->z - cylinder->base->z) * (ray->ori->z - cylinder->base->z) - cylinder->radius * cylinder->radius;

    delta = b * b - 4 * a * c;
    if (delta > ACC)
    {
            root = (-1 * b - sqrt(delta)) / 2 * a - ACC;
            if (root <= ACC)
                    root = (-1 * b + sqrt(delta)) / 2 * a - ACC;
            return (root);
    }
    return (-1);
}

其中 x = (射线->ori->x - 圆柱体->基体->x) et z = (射线->ori->z - 圆柱体->基体->z) et r = 圆柱体->半径。
我计算法 vector 的函数是:

t_vect          *cylinder_normal_at(t_cylinder *cylinder, t_vect *intersection)
{
    t_vect *base_tmp;
    t_vect *normal;
    t_vect *intersection_tmp;

    base_tmp = copy_vector(cylinder->base);
    intersection_tmp = copy_vector(intersection);
    base_tmp->y = intersection_tmp->y;
    normal = init_vector(intersection->x - base_tmp->x, intersection->y - base_tmp->y, intersection-> z - base_tmp->z);
    normalize_vector(normal);
    return (normal);
}

使用这些函数，结果是:Image 1

反射“看起来”不错，形状也不错，但光线不好。正常问题吗？我和我的一个 friend 谈论过这个问题，他告诉我像他一样:在我的第一个函数中删除一些数字。所以它变成:

double          n_ray_cylinder(t_ray *ray, t_cylinder *cylinder, t_data *d)
{
    double a;
    double b;
    double c;
    double delta;
    double root;

    a = ray->dir->x * ray->dir->x + ray->dir->z * ray->dir->z;
    b = ray->dir->x * (ray->ori->x - cylinder->base->x) +
            ray->dir->z * (ray->ori->z - cylinder->base->z);
    c = (ray->ori->x - cylinder->base->x) * (ray->ori->x - cylinder->base->x) +
            (ray->ori->z - cylinder->base->z) * (ray->ori->z - cylinder->base->z) -
            cylinder->radius * cylinder->radius;
    delta = b * b - a * c;
    if (delta > ACC)
    {
            root = (-1 * b - sqrt(delta)) / a - ACC;
            if (root <= ACC)
                    root = (-1 * b + sqrt(delta)) / a - ACC;
            return (root);
    }
    return (-1);
}

正常功能保持不变。然后光线就OK了!

第二个功能似乎工作正常。但为什么 ？第一个不是圆柱方程的精确“应用”吗？

这是第二个问题。我想绕 Ox 轴旋转圆柱体。x²+(y.cosθ+z.sinθ)² = r²。

建立方程后，我得到了这个函数(基于我 friend 的第二个函数)

double          nn_ray_cylinder(t_ray *ray, t_cylinder *cylinder, t_data *d)
{
    double a;
    double b;
    double c;
    double delta;
    double root;
    double deg = M_PI * 45 / 180;

    a = ray->dir->x * ray->dir->x + cos(deg) * cos(deg) * ray->dir->z * ray->dir->z +
            2 * ray->dir->z * cos(deg) * ray->dir->y * sin(deg) + ray->dir->y * ray->dir->y *
            sin(deg) * sin(deg);
    b =  (ray->ori->x - cylinder->base->x) * ray->dir->x + cos(deg) * cos(deg) * (ray->ori->z - cylinder->base->z) * ray->dir->z +
            2 * (ray->ori->z - cylinder->base->z) * cos(deg) * ray->dir->y * sin(deg) + 2 * ray->dir->z * cos(deg) *
            (ray->ori->y - cylinder->base->y) * sin(deg) + 2 * (ray->ori->y - cylinder->base->y) * ray->dir->y * sin(deg) * sin(deg);
    c = (ray->ori->x - cylinder->base->x) * (ray->ori->x - cylinder->base->x) + (ray->ori->z - cylinder->base->z) * (ray->ori->z - cylinder->base->z)* cos(deg) * cos(deg) +
            2 * (ray->ori->z - cylinder->base->z) * cos(deg) * (ray->ori->y - cylinder->base->y) * sin(deg) + (ray->ori->y - cylinder->base->y) * (ray->ori->y - cylinder->base->y) *
            sin(deg) * sin(deg) - cylinder->radius * cylinder->radius;
    delta = b * b -  a * c;
    if (delta > ACC)
    {
            root = (-1 * b - sqrt(delta)) / a - ACC;
            if (root <= ACC)
                    root = (-1 * b + sqrt(delta)) / a - ACC;
            return (root);
    }
    return (-1);
}

轮换成功!但现在我必须改变这个交点上的法线 vector 。为此，我对交点应用反向旋转，然后像以前一样计算法线 vector ，然后重新对该点应用旋转以找到旋转圆柱体的法线 vector 。

t_vect          *cylinder_normal_at(t_cylinder *cylinder, t_vect *intersection)
{
    t_vect *base_tmp;
    t_vect *normal;
    t_vect *intersection_tmp;
    t_matrix *rotate = init_rotation_matrix(45, 1, 0, 0);
    t_matrix *rotate_inverted = init_rotation_matrix(-45, 1, 0, 0);

    base_tmp = copy_vector(cylinder->base);
    intersection_tmp = copy_vector(intersection);
    apply_matrix(rotate_inverted, intersection_tmp);
    base_tmp->y = intersection_tmp->y;
    apply_matrix(rotate, intersection_tmp);
    apply_matrix(rotate, base_tmp);
    normal = init_vector(intersection->x - base_tmp->x,
                    intersection->y - base_tmp->y,
                    intersection->z - base_tmp->z);
    normalize_vector(normal);
    return (normal);
}

例如，对于基本圆柱体和旋转 90 度的圆柱体，结果似乎不错。但是45度的反射完全错误......

Cylinder 45 degrees

那么我的错误在哪里呢？是正常功能有问题，还是其他功能有问题？为什么？非常感谢那些能够帮助我的人。我沉浸在数学中......

编辑:

谢谢 chux，二次错误编码现已更正。法 vector 的问题依然存在。

这里我添加了一些旋转功能:

t_matrix        *init_rotation_matrix(double theta, double x, double y, double z)
{
    t_matrix *matrix;
    double rad;

    if ((matrix = malloc(sizeof *matrix)) == NULL)
            return (NULL);
    rad = theta * M_PI / 180;
    matrix->m11 = x * x * (1 - cos(rad)) + cos(rad);
    matrix->m12 = x * y * (1 - cos(rad)) - z * sin(rad);
    matrix->m13 = x * z * (1 - cos(rad)) + y * sin(rad);
    matrix->m14 = 0;
    matrix->m21 = y * x * (1 - cos(rad)) + z * sin(rad);
    matrix->m22 = y * y * (1 - cos(rad)) + cos(rad);
    matrix->m23 = y * z * (1 - cos(rad)) - x * sin(rad);
    matrix->m24 = 0;
    matrix->m31 = x * z * (1 - cos(rad)) - y * sin(rad);
    matrix->m32 = y * z * (1 - cos(rad)) + x * sin(rad);
    matrix->m33 = z * z * (1 - cos(rad)) + cos(rad);
    matrix->m34 = 0;
    matrix->m41 = 0;
    matrix->m42 = 0;
    matrix->m43 = 0;
    matrix->m44 = 1;
    return (matrix);
}


void    apply_matrix(t_matrix *matrix, t_vect *v)
{
    double x;
    double y;
    double z;

    x = v->x;
    y = v->y;
    z = v->z;
    v->x = matrix->m11 * x + matrix->m12 * y + matrix->m13 * z + matrix->m14;
    v->y = matrix->m21 * x + matrix->m22 * y + matrix->m23 * z + matrix->m24;
    v->z = matrix->m31 * x + matrix->m32 * y + matrix->m33 * z + matrix->m34;
}

编辑2:

在计算法线 vector 的函数中，我将 45 替换为 -45，将 45 替换为 -45。现在它可以工作了...我的 z 轴一定有问题。好像正z和负z颠倒了……

Answer 1

OP 似乎错误地编码了二次方程
(-1 * b - sqrt(delta))/2 * a 而不是
(-1 * b - sqrt(delta))/(2 * a)。

建议使用辅助函数，因为该方程在代码中重复使用。

示例快速、未经测试的代码。

#include <assert.h>
#include <math.h>

int quadratic(double a, double b, double c, double x[2]) {
  if (a == 0.0) return 0;
  double d = b*b - 4*a*c;
  if (d < 0.0) return -1;
  d = sqrt(d);
  x[0] = (-b + d) / (2 * a);
  x[1] = (-b - d) / (2 * a);
  return 2;
}