Raytracing in C - Cylinder, normal vector and rotation

Question

My first post here. :) I'm currently coding a raytracer for my school project in C. I already can display spheres, triangles and planes with some light effects. Now I want to display cylinder (then cones, but cylinder first !). I choose to have a cylinder parallel to the Y-Axis. So I have to solve : x² + z² = r². So technically, my function that returns the distance between my camera and the intersection point is :

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);
}

with x = (ray->ori->x - cylinder->base->x) et z = (ray->ori->z - cylinder->base->z) et r = cylinder->radius.
And my function that calculates the normal vector is:

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);
}

With these functions, the result is : Image 1

The reflects "seems" good, the shape too, but the light is not good. Normal problem ? I talked about this to a friend of mine, and he told me to do like him : delete some numbers in my first function. So it becomes :

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);
}

The normal function stays the same. And then the light is ok !

The second function seems to work fine. But why ? Isn't the first one the exact 'application' of the cylinder equation ?

Here is a second problem. I want to rotate my cylinder around the Ox axis. x²+(y.cosθ+z.sinθ)² = r².

After developing the equation, I have this function (based on my second function, of my 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);
}

The rotation is a success ! But now I have to change the normal vector on this intersection point. For this, I apply an inverted rotation to the intersection point, then I calculate the normal vector like before, then I reapply the rotation to the point to find the normal vector of the rotated cylinder.

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);
}

The result seems good for a basic cylinder and a cylinder rotated by 90 degrees for example. But the reflects are totally wrong with 45 degrees ...

Cylinder 45 degrees

So where are my mistakes ? Is the normal function wrong, or the other one ? And why ? Thanks a lot to those who can help me. I drown myself in mathematics ...

Edit :

Thanks chux, the quadratic mis-coding is now corrected. The problem about the normal vector is still here.

Here I add some my rotation functions :

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;
}

Edit 2 :

In the function that calculates the normal vector, I replaced 45 by -45 and 45 by -45. And now it works ... I must have a problem with my z axis. It seems that the positive z and negative z are inverted ...

Perhaps `(-1 * b - sqrt(delta)) / 2 * a` --> `(-1 * b - sqrt(delta)) / (2 * a)`? Looks like a mis-coding of the quadratic equation. Oui ou non? — chux - Reinstate Monica, Sep 23 '16 at 22:59
Oui. :) Thanks !! But reflections are still wrong at 45 degrees... — Seybol, Sep 23 '16 at 23:48
For the rotation part, normal vectors do not transform the same as tangent vectors. Technically it counts as a Pseudovector or axial vector. (https://en.wikipedia.org/wiki/Pseudovector). A safe way to construct the normal is to take two tangent vectors to the surface, apply the transformation to those and then us the cross product to construct the new normal. A better way is to use the [transpose of the inverse of the matrix] (http://stackoverflow.com/questions/13654401/what-is-the-logic-behind-transforming-normals-with-the-transpose-of-the-inverse) — Salix alba, Sep 24 '16 at 09:34
It sounds a bit difficult ... In my head, my idea is ok. It should work. I need an example so I see why it doesn't work... — Seybol, Sep 24 '16 at 10:33
Small idea: `t_matrix *rotate = init_rotation_matrix(45, 1, 0, 0); ... base_tmp->y = intersection_tmp->y;` --> why does rotation use the `x` unit vector, but the later code adjust `y`? — chux - Reinstate Monica, Sep 24 '16 at 13:52
Imagine a cylinder on Y-axis. Now take an intersection point (any point on the cylinder, it's a infinite cylinder). The base is (0, 0, 0). To get the normal, you take the base, then add intersection->y. The base is now at the same "height" as the intersection point. Then you take the vector between the new base and the intersection point and you got the normal. (after normalizing of course). With the rotation, I do the same, but before all these steps, I rotate the intersection point by -45 degrees. I calculate the new base and then rotate the two points to find my new normal. Why am I wrong? — Seybol, Sep 24 '16 at 14:05
I just did some modifications in my code. I have a problem with my z axis ... — Seybol, Sep 24 '16 at 15:47

score 2 · Answer 1 · answered Sep 23 '16 at 23:31

2

OP appears to have mis-coded the quadratic equation with
(-1 * b - sqrt(delta)) / 2 * a rather than
(-1 * b - sqrt(delta)) / (2 * a).

Suggest using a helper function as that equation is used repetitively in code.

Sample quick, untested code.

#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;
}

answered Sep 23 '16 at 23:31

chux - Reinstate Monica

143,097
13
135
256

Yes, thanks a lot. I mis-coded the quadratic. The lights are now fine !! But when I rotate by 45 degrees, the reflections are still wrong ... nearly the same as in my image... 90 degrees is ok, but another one is wrong. Do I calculate the normal the right way ? – Seybol Sep 23 '16 at 23:46
@Seybol Too many unknowns, unposted code. Perhaps `init_rotation_matrix(` uses radians, not degrees, je ne sais pas. – chux - Reinstate Monica Sep 24 '16 at 00:01
I added my rotation functions in main post. – Seybol Sep 24 '16 at 00:23
@Seybol Cannot help much more now. Suggest trying angle 0 degree first, then 5, then 30 degrees. – chux - Reinstate Monica Sep 24 '16 at 00:48

score 0 · Answer 2 · answered Sep 24 '16 at 18:20

I find my answer. The function n_ray_cylinder() is wrong. Actually, it display a cylinder with a negative angle. The formula that I had was wrong. When I put 45 degrees, it calculated -45 degrees. The right function is :

double          n_ray_cylinder(t_ray *ray, t_cylinder *cylinder, t_data *d)
{
    double a;
    double b;
    double c;
    double delta;
    double root;
    double deg;
    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 =  2 * (ray->ori->x - cylinder->base->x) * ray->dir->x + 2 * 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 -  (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);
}

Some "2" become "-2". Now I have to search how to do 2 another rotations on the same cylinder. Thanks to all !

Raytracing in C - Cylinder, normal vector and rotation

2 Answers2