This project was done by:
Jhonas Leal - https://github.com/Jburlama
Uatilla Almeida - https://github.com/Uatilla
My first RayTracer with miniLibX the last C project from 42 Common Core Program. The main goal of this project is about implementing a raytracing algorithm. We must implement scene aspects as ambient lightining, camera, light, and objects (sphere, plane, cylinder, and cone (bonus only)). Besides these elements, we must apply moviments and rotations into the each element of the scene. For more detailed information, look at the subject of this project.
One of our learnings after building this project is that everyone can apply mathematics or physics without being a mathematician or physician.
If you want to run the project yourself:
- First clone this git repository to your machine.
- Enter the repo and run ``` make ```
- The first time will download and set up the minilibx
- Run ./minirt <scene>.rt
-
As will be explained on Matrices Transformations, every element on the screen has two structures one called transformation and another called inverse these two matrices has the function to deal with modifications on how the scene element is displayed on the screen. Every time the program catches a movement or rotation event, it updates these elements to express the behavior desired by the user:
-
Here you can see the possible movements/rotation for each element:
-
Element type selection:
-
Remember once you hit 'o' all objects were selected, then you can use tab to change between the objects on the screen or click on home to remove the element type selection:
The effects of modidifications on scene elements are:
Center of the image:
Move to the right (using 'c' to select the camera, using right key to move the camera):
Rotating the camera to the left on its own position (not necessary to hit 'c' again, once the camera was already selected, using 'q' key to rotate counter clockwork).
Bringing only the yellow sphere to the camera (using 'o' to select the objects (all objects were selected), hitting 'tab' to select only the sphere and '+' to zoom the sphere).
Moving the light to the left side (using 'l' to select the light and then left arrow to move it to the left).
Norminette:
The Norminette isΒ a tool provided by 42 Network to check whether the source code complies with the school's norm. The rules include the prohibition of βforβ loops, a limit of 25 lines, and no more than 4 parameters in each function, among other specifications.
Here is a link to more details about what is evaluated by the Norminette tool.
MiniLibX:
-
This project uses the minilibx-linux, a API written in c by 42 students, to interact with the X11 window system.
-
A tuple is just an ordered list of things, like numbers.
-
For this project, we used two types of tuples:
- points and vector
-
the way we can tell if a tuple (x, y, z) is a point or a vector, we introduce a new component w;
- (x, y, z, w);
- if w == 0, the tuple is a vector;
- if w == 1, the tuple is a point;
typedef struct s_tuple
{
float x;
float y;
float z;
float w;
} t_tuple;t_tuple creating_tuple(float x, float y, float z, float w)
{
t_tuple ret;
ret.x = x;
ret.y = y;
ret.z = z;
ret.w = w;
return (ret);
}Now that we have tuples, we need to be able to make operations with them.
- Let's say you start from a point (3, -2, 5, 1) and you want to know where you will be if you follow a vector (-2, 3, 1, 0) from that point, you can add those tuples together, and the result will be a new point.
t_tuple sum_tuples(t_tuple *a, t_tuple *b)
{
t_tuple res;
res.x = a->x + b->x;
res.y = a->y + b->y;
res.z = a->z + b->z;
res.w = a->w + b->w;
return (res);
}-
You can find the vector that goes from one point to another by subtracting them.
-
for example:
-
point(3, 2, 1, 1) - point(5, 6, 7, 1) = vector(-2, -4, -6, 0).
-
notice that the w components cancel each other out, resulting in a vector;
-
Similarly, you can subtract a vector (w of 0) from a point (w of 1) and get another tuple with a w of 1, a point. Conceptually, this is just moving backward by the given vector.
-
point(3, 2, 1, 1) - vector(5, 6, 7, 0) = vector(-2, -4, -6, 0)
-
Lastly, subtracting two vectors gives us another vector, representing the change in direction between the two.
-
vector(3, 2, 1, 0) - vector(5, 6, 7, 0) = vector(-2, -4, -6, 0).
t_tuple subtrac_tuples(t_tuple *a, t_tuple *b)
{
t_tuple res;
res.x = a->x - b->x;
res.y = a->y - b->y;
res.z = a->z - b->z;
res.w = a->z - b->w;
return (res);
}- Sometimes youβll want to know what the opposite of some vector is.
- Mathematically, you can get it by subtracting the vector from the tuple (0, 0, 0, 0).
- (0, 0, 0, 0) - (1, -2, 3, 0) = (-1, 2, -3, 0).
t_tuple negating_tuple(t_tuple *a)
{
return (subtrac_tuples(&(t_tuple){0, 0, 0, 0}, a));
}-
Now letβs say you have some vector and you want to know what point lies 3.5 times farther in that direction.
-
Multiplying the vector by 3.5 does just what you need. The 3.5 here is a scalar value because multiplying by it scales the vector (changes its length uniformly).
-
To do it, you multiply each tuple component by the scalar.
-
tuple(1, -2, 3, -4) * 3.5 = (3.5, -7, 10.5, -14).
-
The same works for division
-
tuple(1, -2, 3, -4) / 2 = (0.5, -1, 1.5, -2).
-
Another way of doing division is to divide by a fraction
-
tuple(1, -2, 3, -4) * 0.5 = (0.5, -1, 1.5, -2).
t_tuple mult_tuple_scalar(t_tuple *a, float sc)
{
t_tuple res;
res.x = a->x * sc;
res.y = a->y * sc;
res.z = a->z * sc;
res.w = a->w * sc;
return (res);
}t_tuple div_tuple_scalar(t_tuple *a, float sc)
{
t_tuple res;
res.x = a->x / sc;
res.y = a->y / sc;
res.z = a->z / sc;
res.w = a->w / sc;
return (res);
}- The distance represented by a vector is called its magnitude, or length.
- We can use the Pythagorean Theorem to find the magnitude.
float magnitude(t_tuple *a)
{
t_tuple res;
res.x = a->x * a->x;
res.y = a->y * a->y;
res.z = a->z * a->z;
res.w = a->w * a->w;
return (sqrtf(res.x + res.y + res.z + res.w));
}- Vectors with magnitudes of 1 are called a unit vectors, and these will be super handy.
- Normalization is the process of taking an arbitrary vector and converting it into a unit vector.
- It will keep the calculations anchored relative to a common scale (the unit vector), which is pretty important. If we were to skip normalizing the ray vectors or the surface normals, the calculations would be scaled differently for every ray we cast, and the scenes would look terrible (if they rendered at all).
- You normalize a tuple by dividing each of its components by its magnitude.
t_tuple normalize(t_tuple *a)
{
return (creating_tuple(a->x / magnitude(a), a->y / magnitude(a),
a->z / magnitude(a), a->w / magnitude(a)));
}- The dot product takes two vectors and returns a scalar value.
- Given those two vectors, the dot product is computed as the sum of the products of the corresponding components of each vector.
float dot_product(t_tuple *a, t_tuple *b)
{
t_tuple res;
res.x = a->x * b->x;
res.y = a->y * b->y;
res.z = a->z * b->z;
res.w = a->w * b->w;
return (res.x + res.y + res.z + res.w);
}- the smaller the dot product, the larger the angle between the vectors. For example, given two unit vectors, a dot product of 1 means the vectors are identical, and a dot product of -1 means they point in opposite directions.
- The cross product is another vector operation, but unlike the dot product, it returns another vector instead of a scalar
- For the cross product the order of operands matters because if you change the order you change the direction of the vector
t_tuple cross_product(t_tuple *a, t_tuple *b)
{
return (creating_tuple((a->y * b->z) - (a->z * b->y), \
(a->z * b->x) - (a->x * b->z), \
(a->x * b->y) - (a->y * b->x), 0));
}- This returns a new vector that is perpendicular to both of the original vectors.
- A matrix is a grid of numbers that you can manipulate as a single unit.
- For this project, we used 4Γ4 matrices.
t_matrix *mtx_create(t_minirt *data, int rows, int cols)
{
t_matrix *mtx_struct;
int curr_row;
curr_row = 0;
mtx_struct = (t_matrix *)ft_calloc(sizeof(t_matrix), 1);
if (!mtx_struct)
clean_matrix(data, mtx_struct, errno);
mtx_struct->rows = rows;
mtx_struct->cols = cols;
mtx_struct->mtx = (float **)ft_calloc(sizeof(float *), rows);
if (!mtx_struct->mtx)
clean_matrix(data, mtx_struct, errno);
while (curr_row < rows)
{
mtx_struct->mtx[curr_row] = (float *)ft_calloc(sizeof(float), cols);
if (!mtx_struct->mtx[curr_row])
clean_matrix(data, mtx_struct, errno);
curr_row++;
}
return (mtx_struct);
}- Is with matrices that we'll be able to do transformations in our objects.
- The transformations we'll use are:
- *Translation.
- *Scaling.
- *Rotation.
- The way we multiply matrices is:
- for each element of the resultant matrix, we add the products of the multiplication of each pair from the first matrix row, and the second matrix collum.
t_matrix *mtx_multiply(t_minirt *mrt, t_matrix *mtx_a, t_matrix *mtx_b)
{
t_matrix *mtx_res;
int row;
int col;
if (!mtx_size_compare(mtx_a, mtx_b))
return (NULL);
mtx_res = mtx_create(mrt, mtx_a->rows, mtx_a->cols);
row = -1;
while (++row < mtx_res->rows)
{
col = -1;
while (++col < mtx_res->cols)
mtx_res->mtx[row][col] = mult_mtx_row_col(mtx_a, mtx_b, row, col);
}
clean_matrix(mrt, mtx_b, 0);
clean_matrix(mrt, mtx_a, 0);
return (mtx_res);
}- We are also able to multiply matrices with tuples to produce another tuple.
- We treat the tuple as a 1 column and 4 row matrix.
t_tuple mtx_mult_tuple(t_matrix *mtx_a, t_tuple *tup)
{
t_tuple tup_res;
float *ptr_res;
int i;
ptr_res = (float *)&tup_res;
i = -1;
while (++i < mtx_a->rows)
ptr_res[i] = mult_mtx_row_tuple(mtx_a, tup, i);
return (tup_res);
}- Every shape will have a transformation matrix, and the default without any transformation is the identity matrix.
- If you multiply any matrix by the identity matrix, you'll get the matrix itself.
void fill_idnty_mtx(t_matrix *mtx)
{
int i;
int j;
i = -1;
while (++i < mtx->rows)
{
j = -1;
while (++j < mtx->cols)
{
if (i == j)
mtx->mtx[i][j] = 1;
}
}
}-
The matrix transpose will be used when normal vectors, between object space and world space.
-
The transpose of the identity matrix always gives you the identity matrix.
t_matrix *mtx_transpose(t_minirt *mrt, t_matrix *mtx)
{
t_matrix *mtx_res;
int i;
int j;
int k;
mtx_res = mtx_create(mrt, mtx->rows, mtx->cols);
if (!mtx_res)
clean_matrix(mrt, mtx_res, errno);
i = -1;
while (++i < mtx->rows)
{
j = -1;
while (++j < mtx->cols)
mtx_res->mtx[i][j] = mtx->mtx[j][i];
}
return (mtx_res);
}- If you multiply some matrix A by another matrix B, producing C, you can multiply C by the inverse of B to get A again.
- Well employ a method known as cofactor expansion.
- The determinant is a number that is derived from the elements of a matrix.
- If the determinant is zero, then the corresponding system of equations has no solution.
- A submatrix is what is left when you delete a single row and column from a matrix.
- The submatrix of a 4x4 matrix is 3x3, and the submatrix of a 3x3 matrix is 2x2, and we know how to find the determinant of 2x2 matrices.
t_matrix *submatrix(t_minirt *mrt, t_matrix *mtx, int excl_row, int excl_col)
{
t_matrix *sub_mtx_res;
if (excl_row >= mtx->rows || excl_col >= mtx->cols)
return (NULL);
sub_mtx_res = mtx_create(mrt, mtx->rows - 1, mtx->cols - 1);
if (!sub_mtx_res)
clean_matrix(mrt, sub_mtx_res, errno);
copy_val_to_submtx(mtx, sub_mtx_res, excl_row, excl_col);
return (sub_mtx_res);
}- The minor is the determinant of the submatrix at (i, j).
float minor(t_minirt *mrt, t_matrix *mtx, int row, int col)
{
t_matrix *submtx;
float determ;
submtx = submatrix(mrt, mtx, row, col);
determ = determinant(mrt, submtx);
if (submtx)
clean_matrix(mrt, submtx, 0);
return (determ);
}- Cofactors minors that have (possibly) had their sign changed.
- first well compute the minor at the given row and column. Then well consider that row and column to determine whether or not to negate the result.
- If the row and column identify a spot with a +, then the minorβs sign doesnβt change. If the row and column identify a spot with a -, then you negate the minor.
float cofactor(t_minirt *mrt, t_matrix *mtx, int row, int col)
{
if ((row + col) % 2 == 0)
return (minor(mrt, mtx, row, col));
else
return (minor(mrt, mtx, row, col) * -1);
}- Finding the determinant of matrices larger than 2x2 works recursively.
- To find the determinant, look at any one of the rows or columns. It doesnβt matter which, so letβs just choose the first row.
- Then, for each of those elements, youβll multiply the element by its cofactor, and add the products together.
float determinant(t_minirt *mrt, t_matrix *mtx)
{
float determ;
int i;
int col;
(void)mrt;
determ = 0;
i = -1;
if (mtx->rows == 2 && mtx->cols == 2)
return ((mtx->mtx[0][0] * mtx->mtx[1][1]) -
(mtx->mtx[0][1] * mtx->mtx[1][0]));
else
{
col = -1;
while (++col < mtx->cols)
determ += mtx->mtx[0][col] * cofactor(mrt, mtx, 0, col);
}
return (determ);
}-
Not every matrix is invertible. If the determinant is ever 0, the matrix is not invertible. Anything else is okay.
- We start the algorithm with the construction of a matrix of cofactors.
- Then we transpose the cofactor matrix.
- Finally, divide each of the resulting elements by the determinant of the original matrix.
t_matrix *mtx_inverse(t_minirt *mrt, t_matrix *mtx)
{
t_matrix *mtx_res;
t_matrix *mtx_trans;
float determ;
int row;
int col;
determ = determinant(mrt, mtx);
if (compare_float(determ, 0))
return (NULL);
mtx_res = mtx_create(mrt, mtx->rows, mtx->cols);
if (!mtx_res)
clean_matrix(mrt, mtx_res, errno);
row = -1;
while (++row < mtx->rows)
{
col = -1;
while (++col < mtx->cols)
mtx_res->mtx[row][col] = cofactor(mrt, mtx, row, col) / determ;
}
mtx_trans = mtx_transpose(mrt, mtx_res);
clean_matrix(mrt, mtx_res, 0);
return (mtx_trans);
}- Calculating the inverse of a matrix is very computationally intense, is more efficient to do it before the render.
- With transformation matrices will be able to move and deform objects.
- The below functions take the identity matrix as a parameter.
- Translation is a transformation that moves a point by adding to it.
void mtx_translation(t_matrix *mtx, t_tuple *tup_transl)
{
int row;
float *tup_arr;
tup_arr = (float *)tup_transl;
row = -1;
while (++row < 3)
mtx->mtx[row][3] = tup_arr[row];
}void mtx_scaling(t_matrix *mtx, t_tuple *tup_scale)
{
int row;
int col;
float *tup_arr;
tup_arr = (float *)tup_scale;
row = -1;
while (++row < 3)
{
col = -1;
while (++col < 3)
{
if (row == col)
mtx->mtx[row][col] = tup_arr[row];
}
}
}- Multiplying a tuple by a rotation matrix will rotate that tuple around an axis.
- The rotation function will take as a parameter the angle in radians.
void mtx_rotation_x(t_matrix *mtx, float rot_deg)
{
mtx->mtx[1][1] = cosf(rot_deg);
mtx->mtx[1][2] = sinf(rot_deg) * -1;
mtx->mtx[2][1] = sinf(rot_deg);
mtx->mtx[2][2] = cosf(rot_deg);
}void mtx_rotation_y(t_matrix *mtx, float rot_deg)
{
mtx->mtx[0][0] = cosf(rot_deg);
mtx->mtx[0][2] = sinf(rot_deg);
mtx->mtx[2][0] = sinf(rot_deg) * -1;
mtx->mtx[2][2] = cosf(rot_deg);
}void mtx_rotation_z(t_matrix *mtx, float rot_deg)
{
mtx->mtx[0][0] = cosf(rot_deg);
mtx->mtx[0][1] = sinf(rot_deg) * -1;
mtx->mtx[1][0] = sinf(rot_deg);
mtx->mtx[1][1] = cosf(rot_deg);
}- We can perform multiple transformations by multiplying the above matrices.
- Note that the order of the multiplications is important! Matrix multiplication is associative, but not commutative.
- if you want a single matrix that rotates, and then scales, and then translates, you can multiply the translation matrix by the scaling matrix, and then by the rotation matrix. That is to say, you must concatenate the transformations in reverse order to have them applied in the order you want!
t_matrix *mtx_multiply(t_minirt *mrt, t_matrix *mtx_a, t_matrix *mtx_b)
{
t_matrix *mtx_res;
int row;
int col;
if (!mtx_size_compare(mtx_a, mtx_b))
return (NULL);
mtx_res = mtx_create(mrt, mtx_a->rows, mtx_a->cols);
row = -1;
while (++row < mtx_res->rows)
{
col = -1;
while (++col < mtx_res->cols)
mtx_res->mtx[row][col] = mult_mtx_row_col(mtx_a, mtx_b, row, col);
}
clean_matrix(mrt, mtx_b, 0);
clean_matrix(mrt, mtx_a, 0);
return (mtx_res);
}- Ray casting is the process of creating a ray, or line, and finding the intersections of that ray with the objects in a scene
- Each ray will have a starting point called the origin, and a vector called the direction which says where it points.
typedef struct s_ray
{
t_point origin;
t_vector direction;
} t_ray;- With the origin (O) and direction (D), we can find any point (P) that lies long the ray, t units.
- P = O + Dt
t_point position(t_ray *ray, float t)
{
t_point point;
point = mult_tuple_scalar(&ray->direction, t);
point = sum_tuples(&point, &ray->origin);
return (point);
}-
when casting a ray, there are three possible different scenarios:
- the ray misses the sphere and has 0 intersections.
- the ray hits a tangent and only has one intersection.
- the hits two different points in the sphere.
-
Assuming the sphere is centered at the world origin (0, 0, 0) with radii of 1.
-
For any point (P), if the point minus the sphere center (C) is equal to the radius of the sphere (1), then the point belongs to the sphere.
$P - C = 1$
-
To find if a ray intersects the sphere, we check if the ray has any possible t, being P a point in the ray:
$(O + Dt) - C = 1$
-
$(O + Dt)^2 - C^2 = 1^2 =>$ -
$O^2 + 2ODt + (Dt)^2 - C^2 = 1 =>$ -
$O^2 + 2ODt + D^2 * t^2 - C^2 - 1 = 0 =>$ -
$2ODt + D^2 * t^2 + O^2 - C^2 - 1 = 0$ -
$O - C$ is a vector from the ray origin and the sphere center, let's call W. -
$2ODt + D^2 * t^2 + W^2 - 1 = 0$ => -
$D^2 * t^2 + 2ODt + W^2 - 1 = 0$ -
What we have left is a quadratic equation, and we can solve it in function of t:
-
$a = D^2$ -
$b = 2OD$ -
$c = W^2 - 1$ -
a = the dot product of the ray direction with the ray direction
-
b = 2 * the dot product of the ray direction with the ray origin
-
c = the dot product of w with w minus 1.
-
First, we need to solve what's inside the square root, to see if the equation has a solution.
- if we get a negative value, the ray misses the sphere.
- if we get zero than we will only have one solution, and the ray hits a tangent.
- if we get a value bigger than zero, the ray intercepts the sphere in two points.
int8_t ray_sphere_intersect(t_ray *ray, float *t)
{
float a;
float b;
float c;
t_tuple sphere_to_ray;
float discriminant;
if (!ray || !t)
return (0);
sphere_to_ray = subtrac_tuples(&ray->origin, &(t_point){0, 0, 0, 1});
a = dot_product(&ray->direction, &ray->direction);
b = 2 * dot_product(&ray->direction, &ray->origin);
c = dot_product(&sphere_to_ray, &sphere_to_ray) - 1;
discriminant = (b * b) - 4 * a * c;
if (discriminant < 0)
return (0);
t[0] = (-b - sqrtf(discriminant)) / (2 * a);
t[1] = (-b + sqrtf(discriminant)) / (2 * a);
if (t[0] == t[1])
return (1);
return (2);
}-
A plane is a perfectly flat surface that extends infinitely in two dimensions, in our case, the plane will extend in both x and z dimensions, passing through the origin. Using the transformation matrix will be able to rotate and translate the plane in any orientation we like.
-
The normal vector is the same in every point of the plane, (0,1,0).
-
If the ray is parallel to the plane there will be no intersections.
-
Also no intersection if the ray is coplanar with the plane, which is to say that the rayβs origin is on the plane, and the rayβs direction is parallel to the plane.
int8_t ray_plane_intersect(t_ray *ray, float *t)
{
if (fabs(ray->direction.y) < EPSILON)
return (0);
t[0] = -ray->origin.y / ray->direction.y;
t[1] = -ray->origin.y / ray->direction.y;
return (1);
}-
At first, we assume the cylinder goes extends infinitely in the y direction, and have a radius of 1.
-
We can check if the ray intersects the cylinder by joining the Ray formula
$(O + Dt)$ with the circle formula$(x^2 + z^2 = 1)$ ; -
$(Ox + Dxt)^2 + (Oz + Dzt)^2 = 1$ -
$Ox^2 + 2OxDxt + Dxt^2 + Oz^2 + 2OzDzt + Dzt^2 = 1$ -
$Dxt^2 + Dzt^2 + 2OxDxt + 2OzDzt + Ox^2 + Oz^2 = 1$ -
$t^2(Dx^2 + Dz^2) + t(OxDxt + OzDzt) + Ox^2 + Oz^2 - 1 = 0$ -
Again just like the circle we can use the quadratic equation to solve in function of t:
$a = Dx^2 + Dz^2$ $b = 2OxDxt + 2OzDzt$ $c = Ox^2 + Oz^2 - 1$
int8_t ray_cylinder_intersect(t_ray *ray, float *t, t_shape *obj)
{
float a;
float b;
float c;
float discriminant;
bool count[2];
a = (ray->direction.x * ray->direction.x) + (ray->direction.z
* ray->direction.z);
if (fabs(a) < EPSILON)
return (0);
b = (2 * ray->origin.x * ray->direction.x) + (2 * ray->origin.z
* ray->direction.z);
c = (ray->origin.x * ray->origin.x) + (ray->origin.z
* ray->origin.z) - 1;
discriminant = (b * b) - 4 * a * c;
if (discriminant < 0)
return (0);
t[0] = (-b - sqrt(discriminant)) / (2 * a);
t[1] = (-b + sqrt(discriminant)) / (2 * a);
if (t[0] > t[1])
swap(t);
count[0] = check_cy_range(ray, t[0], obj);
count[1] = check_cy_range(ray, t[1], obj);
return (intercections_count(count, t));
}- before checking the intersections we cap the cylinder in a range, set in the object materials
- Then in the check_cy_range we check if the intersection point is between that range
bool check_cy_range(t_ray *ray, float t, t_shape *obj)
{
float y;
float cap_t;
y = ray->origin.y + (ray->direction.y * t);
if (obj->material.min < y && y < obj->material.max)
return (true);
return (false);
}- The intersections_count, will return how many points the ray hits.
int8_t intercections_count(bool *count, float *t)
{
if (count[0] && count[1])
{
if (t[0] == t[1])
return (1);
else
return (2);
}
else if (!count[0] && !count[1])
return (0);
else if (count[0])
{
t[1] = t[0];
return (1);
}
else if (count[1])
{
t[0] = t[1];
return (1);
}
return (0);
}- For the top and bottom of the cylinder we treated them like another object, planes. All we need to check is if it's within the y set in the t_material, and if it's within the radius of the circle.
int8_t ray_cap_inter(t_ray *ray, float *t, t_shape *obj)
{
bool count[2];
if (obj->material.closed == false || fabs(ray->direction.y) < EPSILON)
return (0);
t[0] = (obj->material.min - ray->origin.y) / ray->direction.y;
t[1] = (obj->material.max - ray->origin.y) / ray->direction.y;
count[0] = check_cap(ray, t[0], obj, 1);
count[1] = check_cap(ray, t[1], obj, 2);
return (intercections_count(count, t));
}
bool check_cap(t_ray *ray, float t, t_shape *obj, int8_t order)
{
float x;
float z;
x = ray->origin.x + (t * ray->direction.x);
z = ray->origin.z + (t * ray->direction.z);
return (((x * x) + (z * z)) <= 1);
}- The cone works just like the cylinder but the a, b, and c are computed differently:
$a = d2x β d2y + d2z$ $b = 2ox dx β 2oy dy + 2oz dz$ $c = o2x β o2y + o2z$
int8_t ray_cone_intersect(t_ray *ray, float *t, t_shape *obj)
{
float a;
float b;
float c;
float discriminant;
bool count[2];
set_cone_intersect(&a, &b, &c, ray);
if (a == 0 && b == 0)
return (0);
if (a == 0)
{
t[0] = -c / 2 * b;
t[1] = -c / 2 * b;
return (1);
}
discriminant = (b * b) - 4 * a * c;
if (discriminant < 0)
return (0);
t[0] = (-b - sqrt(discriminant)) / (2 * a);
t[1] = (-b + sqrt(discriminant)) / (2 * a);
if (t[0] > t[1])
swap(t);
count[0] = check_cy_range(ray, t[0], obj);
count[1] = check_cy_range(ray, t[1], obj);
return (intercections_count(count, t));
}
void set_cone_intersect(float *a, float *b, float *c, t_ray *ray)
{
*a = (ray->direction.x * ray->direction.x) - (ray->direction.y
* ray->direction.y) + (ray->direction.z * ray->direction.z);
*b = (2 * ray->origin.x * ray->direction.x) - (2 * ray->origin.y
* ray->direction.y) + (2 * ray->origin.z * ray->direction.z);
*c = (ray->origin.x * ray->origin.x) - (ray->origin.y
* ray->origin.y) + (ray->origin.z * ray->origin.z);
}- This algorithm will check for intersections on a double-napped cone.
- If both a and b are 0, then the ray misses the cone.
- When a is 0 and b isn't the ray is parallel to one of the cone halves.
-
This still means that the ray will intersect the other half, so we use the following formula to find the single point of intersection:
$t = -c / 2 * b$
-
To check the cap for the cone we use the same logic as the cylinder, but the radius of the circle is not the same everywhere, and will change along the y coordinates.
bool check_cap(t_ray *ray, float t, t_shape *obj, int8_t order)
{
float x;
float z;
x = ray->origin.x + (t * ray->direction.x);
z = ray->origin.z + (t * ray->direction.z);
if (obj->type == CY)
return (((x * x) + (z * z)) <= 1);
else if (obj->type == CONE)
{
if (order == 1)
return (((x * x) + (z * z)) <= fabs(obj->material.min \
* obj->material.min));
else if (order == 2)
return (((x * x) + (z * z)) <= fabs(obj->material.max \
* obj->material.max));
}
return (false);
}- All shapes will be represented by the t_shape struct, which will be put into a linked list.
typedef struct s_shape
{
t_material material;
t_ray trans_ray;
t_matrix *mtx_trans;
t_matrix *mtx_inver;
t_point center;
t_angle angle;
enum e_id type;
float amb_ratio;
int id;
void *next;
} t_shape;
typedef t_shape t_sphere;
typedef t_shape t_plane;
typedef t_shape t_cyl;
typedef t_shape t_cone;- They will all have in common:
- a material:
typedef struct s_material
{
t_pattern pattern;
t_color color;
t_color color_sec;
float ambient;
float diffuse;
float specular;
float shininess;
float reflective;
float min;
float max;
bool closed;
bool is_bump;
} t_material;- The min, max, and closed values will be only used by the cone and cylinder.
- a transformation matrix, which will be set as the identity matrix, if the object doesn't have a transformation, and the inverse of that matrix.
- The transformed ray. The ray that multiplied by the inverse of the object transformation matrix.
- Every time a ray hits an object will add then to a list of intersections
typedef struct s_intersections
{
t_point point;
float t[2];
float hit;
t_shape *obj;
struct s_intersections *next;
} t_intersections;-
after check if the ray intersects with all the objects, we keep a reference to the intersection with the smallest positive hit value, being that the first point the ray hits in a object.
-
The intersection's node keeps track of:
- the t values the ray hits, if the ray only hits one point, t[0] will be equal to t[1];
- the t value the ray hit first in hit;
- the point of intersection;
- a pointer to the object the ray intersects;
- You can notice that the calculations for the ray intersections will always assume that the object is in the center of the world, point(0,0,0), and with radii of 1 in the case of the sphere.
- The way we move objects is not by moving the object itself but the ray.
- before we check the intersections, we first transform the ray by multiplying it by the object's inverse matrix.
- This will give us the correct location of the points of the objects.
t_ray ray_trasform(t_ray *ray, t_matrix *mtx)
{
t_ray new_ray;
new_ray.origin = mtx_mult_tuple(mtx, &ray->origin);
new_ray.direction = mtx_mult_tuple(mtx, &ray->direction);
return (new_ray);
}void intersections(t_minirt *data, t_ray *ray)
{
t_shape *obj;
obj = data->world.objs;
while (obj)
{
obj->trans_ray = ray_trasform(ray, obj->mtx_inver);
ray_intersections(data, obj, &obj->trans_ray, ray);
obj = obj->next;
}
first_hit(data);
}We use the rgb representation to determine the color of the pixel.
- Each value will range from [0-1] being 0, completely absent of that color, and 1 the full value of that color.
- We will represent color as a tuple, and we will also need to be able to do tuple operations with them.
typedef struct s_tuple
{
union
{
struct
{
float x;
float y;
float z;
float w;
};
struct
{
float r;
float g;
float b;
};
};
} t_tuple;
typedef t_tuple t_point;
typedef t_tuple t_vector;
typedef t_tuple t_color;- We used union to make the code more readable, and more distinguishable when we are working with r,g,b colors or x,y,z coordinates.
- Every call to point light will add a new light source to the world and will put it on a list of lights.
void point_light(t_point *pos, t_color *intensity, t_world *world)
{
t_light *light;
if (world->light == NULL)
{
world->light = ft_calloc(sizeof(*world->light), 1);
if (world->light == NULL)
exit (errno);
world->light->position = *pos;
world->light->intensity = *intensity;
world->light->next = NULL;
return ;
}
light = ft_calloc(sizeof(*light), 1);
if (light == NULL)
exit(errno);
light->position = *pos;
light->intensity = *intensity;
light->next = world->light;
world->light = light;
}-
To simulate the reflection of light on a surface, we'll implement the Phong algorithm.
-
For the algorithm, we need 4 vectors.
-
The eye vector. P (Points from the point of intersection), to the origin of the ray.
- To find E, we negate the rayβs direction vector, turning it around to point back at its origin.
-
The light vector. Pointing from P to the origin of the light source.
- To find L, we subtract P from the position of the light source, giving us the vector pointing toward the light.
-
The surface normal, a vector that is perpendicular to the surface at P.
-
The reflection vector, pointing in the direction that incoming light would bounce, or reflect.
-
- First, we need to undo the transformations of the point, making it easier to calculate the normal.
- We achieve that by multiplying the point by the object inverse matrix.
- The reason we use the inverse matrix, instead of the object transformation matrix, is because the transformation matrix wont preserve the correct y component, changing the size of the vector.
- So we multiply the normal by the inverse transpose matrix instead.
t_vector normal_at(t_shape *obj, t_point *point, t_minirt *data)
{
t_point local_point;
t_vector local_normal;
t_vector world_normal;
t_matrix *transpose;
local_point = mtx_mult_tuple(obj->mtx_inver, point);
local_normal = local_normal_at(obj, &local_point);
transpose = mtx_transpose(data, obj->mtx_inver);
world_normal = mtx_mult_tuple(transpose, &local_normal);
world_normal.w = 0;
clean_matrix(data, transpose, 0);
return (normalize(&world_normal));- after undoing the transformations of the point, we will calculate the normal for each object.
-
We subtract the point by the center, giving us the normalized normal of the tangent plane of that point.
if (obj->type == SP)
local_normal = subtrac_tuples(local_point, &(t_point){0, 0, 0, 1});- the local normal at a plane is always perpendicular to the y axis.
else if (obj->type == PL)
local_normal = (t_vector){0, 1, 0, 0};t_vector normal_at_cone(t_point *point, t_shape *obj)
{
float dist;
float y;
t_vector normal;
dist = (point->x * point->x) + (point->z * point->z);
if (dist < obj->material.max * obj->material.max \
&& point->y >= obj->material.max - EPSILON)
return ((t_vector){0, 1, 0, 0});
else if (dist < obj->material.min * obj->material.min \
&& point->y <= obj->material.min + EPSILON)
return ((t_vector){0, -1, 0, 0});
y = sqrtf(dist);
if (y > EPSILON)
y = -y;
normal = (t_vector){point->x, y, point->z, 0};
return (normal);
}- This is the reflected vector around the normal.
- the vector coming in will have the same angle as the vector going out.
t_vector reflect(t_vector *in, t_vector *normal)
{
t_tuple vect;
vect = mult_tuple_scalar(normal, 2);
vect = mult_tuple_scalar(&vect, dot_product(in, normal));
vect = subtrac_tuples(in, &vect);
return (vect);
}-
To simulate the reflection of light we implement the Phong reflection model.
-
This algorithm simulates the interaction between 3 types of lights.
- Ambient: Is the background lightning, this value will be a constant.
- Diffuse: Is the light reflected from the material surface, this only depends on the angle between the light source and the surface normal.
- Specular: Is the reflection of the light source itself, the bright spot on a curved surface, It depends only on the angle between the reflection vector and the eye vector and is controlled by a parameter that we call shininess. The higher the shininess, the smaller and tighter the specular highlight.
typedef struct s_phong
{
t_color ambient;
t_color diffuse;
t_color spec;
} t_phong;t_color lighting(t_comps *comps, t_light *light, t_world *world)
{
t_color color;
t_phong phong;
float light_normal_dot;
color = color_multiply(&comps->obj->material.color, &world->ambient_light);
comps->lightv = subtrac_tuples(&light->position, &comps->point);
comps->lightv = normalize(&comps->lightv);
phong.ambient = mult_tuple_scalar(&color, comps->obj->material.ambient);
light_normal_dot = dot_product(&comps->lightv, &comps->normalv);
if (light_normal_dot < 0 || comps->is_shadown)
light_is_behind_obj(&phong.diffuse, &phong.spec);
else
{
color = color_multiply(&color, &light->intensity);
phong.diffuse = mult_tuple_scalar(&color,
comps->obj->material.diffuse * light_normal_dot);
bump(&phong, comps->obj);
comps->reflectv = negating_tuple(&comps->lightv);
comps->reflectv = reflect(&comps->reflectv, &comps->normalv);
get_specular(comps, light, &phong);
}
return (add_color3(&phong.ambient, &phong.diffuse, &phong.spec));
}
void light_is_behind_obj(t_color *diffuse, t_color *specular)
{
*diffuse = (t_color){0, 0, 0, 999999};
*specular = (t_color){0, 0, 0, 999999};
}
void get_specular(t_comps *comps, t_light *light, t_phong *phong)
{
float ref_dot_eye;
ref_dot_eye = dot_product(&comps->reflectv, &comps->eyev);
if (ref_dot_eye <= 0)
phong->spec = (t_color){0, 0, 0, 999999};
else
phong->spec = specular(&comps->obj->material, light, ref_dot_eye);
}- If the point is in shadow, in the lightning function, we will ignore the specular and diffuse values and use only the ambient to color the point.
- We will cast a shadow ray, going through the point of intersection to the light source, and if that ray intersects an object, then that point is on shadow.
bool is_shadowed(t_world *w, t_light *light, t_point *p)
{
t_minirt data;
t_vector v;
float distance;
t_ray ray;
ft_memset(&data, 0, sizeof(data));
data.world = *w;
v = subtrac_tuples(&light->position, p);
distance = magnitude(&v);
ray.direction = normalize(&v);
ray.origin = *p;
intersections(&data, &ray);
if (data.first_hit && distance > data.first_hit->hit)
{
clear_ray_inter(&data);
return (true);
}
clear_ray_inter(&data);
return (false);
}-
We Measure the distance from point to the light source by subtracting point from the light position and taking the magnitude of the resulting vector. Call this distance.
-
We Create a ray from point toward the light source by normalizing the vector from step 1.
-
We intersect the world with that ray.
-
We check to see if there was a hit, and if so, whether t is less than distance. If so, the hit lies between the point and the light source, and the point is in shadow.
- The camera will represent our eye.
- The camera will have a transformation matrix, that will pretty much orientate the world relative to our eye.
- The view_transformation function will return the camera transformation matrix and take 3 parameters.
- From: The point of origin of the camera.
- To: To were the camera is pointing to.
- Up: The vector indicating which direction is up.
typedef struct s_view
{
t_vector forward;
t_vector upn;
t_vector left;
t_vector true_uper;
} t_view;t_matrix *view_transformation(t_point *from, t_point *to, t_vector *up)
{
t_view view;
t_matrix *orientation;
t_matrix *trans;
t_matrix *ret;
view.forward = subtrac_tuples(to, from);
view.forward = normalize(&view.forward);
view.upn = normalize(up);
view.left = cross_product(&view.forward, &view.upn);
view.true_uper = cross_product(&view.left, &view.forward);
orientation = view_orientation(&view.left, &view.true_uper, &view.forward);
trans = mtx_create(NULL, 4, 4);
fill_idnty_mtx(trans);
mtx_translation(trans, &(t_point){-from->x, -from->y, -from->z, 1});
ret = mtx_multiply(NULL, orientation, trans);
return (ret);
}t_matrix *view_orientation(t_vector *left, t_vector *up, t_vector *forward)
{
t_matrix *orientation;
orientation = mtx_create(NULL, 4, 4);
orientation->mtx[0][0] = left->x;
orientation->mtx[0][1] = left->y;
orientation->mtx[0][2] = left->z;
orientation->mtx[1][0] = up->x;
orientation->mtx[1][1] = up->y;
orientation->mtx[1][2] = up->z;
orientation->mtx[2][0] = -forward->x;
orientation->mtx[2][1] = -forward->y;
orientation->mtx[2][2] = -forward->z;
orientation->mtx[3][3] = 1;
return (orientation);
}-
We compute the forward vector by subtracting from, from to, and normalizing the result.
-
We compute the left vector by taking the cross product of forward and the normalized up vector.
-
We compute the true_up vector by taking the cross product of left and forward. This allows the original up vector to be only approximately up, which makes framing scenes a lot easier.
-
With these left, true_up, and forward vectors, we can now construct a matrix that represents the orientation transformation:
![[Pasted image 20240926174653.png]]
-
The last step is to append a translation to that transformation to move the scene into place before orienting it. Multiplying orientation by translation(-from.x, -from.y, -from.z).
typedef struct s_camera
{
t_matrix *trans;
t_matrix *inver;
t_point center;
t_point direct_center;
t_vector up;
int hsize;
int vsize;
float half_width;
float half_height;
float pixel_size;
float fov;
} t_camera;- The camera is defined by the following four attributes:
- hsize: is the horizontal size (in pixels) of the canvas that the picture will be rendered to.
- vsize: is the canvasβs vertical size (in pixels).
- field_of_view (fov): is an angle that describes how much the camera can see.
- transform (trans): is a matrix describing how the world should be oriented relative to the camera. We achieve this with view transformation.
t_camera camera_construct(size_t hsize, size_t vsize, float fov)
{
t_camera camera;
float half_view;
float aspect;
camera.hsize = hsize;
camera.vsize = vsize;
camera.fov = fov;
half_view = tan(fov / 2);
aspect = (float)hsize / (float)vsize;
if (aspect >= 1)
{
camera.half_width = half_view;
camera.half_height = half_view / aspect;
}
else
{
camera.half_width = half_view * aspect;
camera.half_height = half_view;
}
camera.pixel_size = (camera.half_width * 2) / hsize;
return (camera);- The ray_for_pixel() function computes the world coordinates at the center of the given pixel, and then constructs a ray that passes through that point.
- It takes as a parameter the camera, px (the x position of the pixel), and py (the y position of the pixel).
t_ray ray_for_pixel(t_camera *camera, size_t px, size_t py)
{
t_ray ray;
t_point pixel;
t_rfp rpf;
rpf.xoffset = (px + 0.5) * camera->pixel_size;
rpf.yoffset = (py + 0.5) * camera->pixel_size;
rpf.w_x = camera->half_width - rpf.xoffset;
rpf.w_y = camera->half_height - rpf.yoffset;
pixel = mtx_mult_tuple(camera->inver, &(t_point){rpf.w_x, rpf.w_y, -1, 1});
ray.origin = mtx_mult_tuple(camera->inver, &(t_point){0, 0, 0, 1});
ray.direction = subtrac_tuples(&pixel, &ray.origin);
ray.direction = normalize(&ray.direction);
return (ray);
}- The render function will iterate, through the pixels of the scream.
- With the help of the functions ray_for_pixel(), and color_at(), will create the ray, and then compute the color of that pixel.
- Then write the pixel in the window using the minilibx API.
void render(t_minirt *data)
{
t_color color;
t_ray ray;
int x;
int y;
y = -1;
while (++y < data->camera.vsize - 1)
{
x = -1;
while (++x < data->camera.hsize - 1)
{
ray = ray_for_pixel(&data->camera, x, y);
color = color_at(data, &ray, 5);
write_pixel(&data->canvas, x, y, &color);
}
}
printf("RENDER\t\t[OK]\n");
manage_interface(data);
}- The reflection works by spawning an additional ray at the point of intersection and recursively following it to determine the color at that point.
t_color reflected_color(t_comps *comps, t_minirt *data, int8_t remaining)
{
t_color color;
t_ray reflected_ray;
if (comps->obj->material.reflective == 0 || remaining == 0)
return ((t_color){0, 0, 0, 999999});
reflected_ray.origin = comps->over_point;
reflected_ray.direction = comps->reflectv;
color = color_at(data, &reflected_ray, remaining - 1);
color = mult_tuple_scalar(&color, comps->obj->material.reflective);
return (color);
}- The remaining parameter, is do that we don't have an infinity recursion if we have two reflective objects pointing at each other.
- the reflective attribute in the object material, will determine how reflective is the object, being 0 non-reflective and 1 a perfect mirror.
- We'll add the color returning from this function to the color of the object intersected, to have the final color.
-
If the object has a pattern, will change the default color of the object according to the point, and the pattern the object has.
-
The possible patterns the objects can have are:
- STR, stripe pattern;
- PC, point color;
- GR, gradient;
- RNG, ring;
- CHK, checkers.
enum e_p
{
STR = 0,
PC = 1,
GR = 2,
RNG = 3,
CHK = 4,
};void set_pattern(t_intersections *inter, t_point *point)
{
if (inter->obj->material.pattern.has)
{
inter->obj->material.color = pattern_at(&inter->obj->material.pattern, \
point, inter->obj, inter->obj->material.pattern.type);
}
}- As the x coordinate changes, the pattern alternates between the two colors.
t_color stripe_at(t_pattern *patterns, t_point *point)
{
t_color c;
if ((int)point->x % 2 == 0)
c = patterns->a;
else
c = patterns->b;
return (c);
}- A gradient pattern is like stripes, but instead of discrete steps from one color to the next, the function returns a blend of the two colors, linearly interpolating from one to the other as the x coordinate changes.
t_color gradient(t_pattern *pattern, t_point *point)
{
t_color distance;
float fraction;
t_color ret;
distance = subtrac_tuples(&pattern->a, &pattern->b);
fraction = point->x - floor(point->x);
ret = mult_tuple_scalar(&distance, fraction);
return (sum_tuples(&pattern->a, &ret));
}- A ring pattern depends on two dimensions, x, and z, to decide which color to return. It works similarly to stripes, but instead of testing the distance of the point in just x, it tests the distance of the point in both x and z.
t_color ring_patt(t_pattern *pattern, t_point *point)
{
t_color c;
if ((int)sqrtf((point->x * point->x) + (point->z * point->z)) % 2 == 0)
c = pattern->a;
else
c = pattern->b;
return (c);
}- The function of this pattern is very much like that of stripes, but instead of relying on a single dimension, it relies on the sum of all three dimensions, x, y, and z.
t_color checker_patt(t_pattern *pattern, t_point *point)
{
t_color c;
if ((int)(floor(point->x) + floor(point->y) + floor(point->z)) % 2 == 0)
c = pattern->a;
else
c = pattern->b;
return (c);
}This project would not be possible without the guidelines of the book:
The Ray Tracer Challenger - by Jamis Buck
Jamis Buck - The Ray Tracer Challenge | PDF
We strongly recommend to anyone who wants to build this project from scratch, that deep dive on the reading of this book to understand the concepts, there is a lot of valuable information there.
Contributions to the project are welcome! If you have any ideas, improvements, or bug fixes, please submit them as issues or pull requests to this repository.
- π I'm
Uatilla Viana Almeida. - π± I'm currently studying Common Core at 42 Porto after successfully completing the Piscine.
- π I'm interested in
Bitcoin,Macroeconomy,AdventureandSportsof all kinds. - π Additionally, I have an interest in
Blockchain,Python,Data AnalysisandMachine Learning. - π« You can reach me on LinkedIn.
- π€ Feeling amazing about how the blockchain technology will change our lives.









































