cozis/ray_tracing
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

first commit

Browse Source
This commit is contained in:
cozis
2024-09-16 13:13:56 +02:00
commit ddd25a798e
245 changed files with 46862 additions and 0 deletions
Download Patch File Download Diff File Expand all files Collapse all files
src/camera.c
+93
View File
@@ -0,0 +1,93 @@
#include <math.h>
#include "camera.h"
static bool first_mouse = true;
static float yaw = -90.0f;
static float pitch = 0.0f;
static float last_x = 800.0f / 2.0;
static float last_y = 600.0f / 2.0;
static float fov = 30.0f;
static float delta_time = 0.0f;
static float last_frame = 0.0f;
static Vector3 camera_pos = {5, 5, 5};
static Vector3 camera_front = {-1, -1, -1};
static Vector3 camera_up = {0, 1, 0};
Vector3 get_camera_pos(void)
{
return camera_pos;
}
void rotate_camera(double mouse_x, double mouse_y)
{
float x = mouse_x;
float y = mouse_y;
if (first_mouse) {
last_x = x;
last_y = y;
first_mouse = false;
}
float dx = x - last_x;
float dy = last_y - y;
last_x = x;
last_y = y;
float sensitivity = 0.1f;
dx *= sensitivity;
dy *= sensitivity;
yaw += dx;
pitch += dy;
if (pitch > 89.0f) pitch = 89.0f;
if (pitch < -89.0f) pitch = -89.0f;
float yaw_rad = deg2rad(yaw);
float pitch_rad = deg2rad(pitch);
Vector3 front;
front.x = cos(yaw_rad) * cos(pitch_rad);
front.y = sin(pitch_rad);
front.z = sin(yaw_rad) * cos(pitch_rad);
front = normalize(front);
camera_front = front;
}
void move_camera(Direction dir, float speed)
{
switch (dir) {
case UP : camera_pos = combine(camera_pos, camera_front, 1, +speed); break;
case DOWN : camera_pos = combine(camera_pos, camera_front, 1, -speed); break;
case LEFT : camera_pos = combine(camera_pos, normalize(cross(camera_front, camera_up)), 1, -speed); break;
case RIGHT: camera_pos = combine(camera_pos, normalize(cross(camera_front, camera_up)), 1, +speed); break;
}
}
Matrix4 camera_pov(void)
{
return lookat_matrix(camera_pos, combine(camera_pos, camera_front, 1, 1), camera_up);
}
Ray ray_through_screen_at(float px, float py, float aspect_ratio)
{
Vector3 w = normalize(scale(camera_front, -1));
Vector3 u = normalize(cross(camera_up, w));
Vector3 v = cross(w, u);
float screen_h = 2 * tan(fov / 2);
float screen_w = aspect_ratio * screen_h;
Vector3 horizontal = scale(u, screen_w);
Vector3 vertical = scale(v, screen_h);
Vector3 lower_left_corner = combine4(camera_pos, horizontal, vertical, w, 1, -0.5, -0.5, -1);
Vector3 dir = combine4(lower_left_corner, horizontal, vertical, camera_pos, 1, px, py, -1);
return (Ray) {.origin=camera_pos, .direction=dir};
}
src/camera.h
+11
View File
@@ -0,0 +1,11 @@
#include "vector.h"
typedef enum {
UP, DOWN, LEFT, RIGHT,
} Direction;
Matrix4 camera_pov(void);
void move_camera(Direction dir, float speed);
void rotate_camera(double mouse_x, double mouse_y);
Vector3 get_camera_pos(void);
Ray ray_through_screen_at(float u, float v, float aspect_ratio);
src/main.c
+573
View File
@@ -0,0 +1,573 @@
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <assert.h>
#include <float.h> // FLT_MAX
#include <glad/glad.h>
//#define GLFW_INCLUDE_NONE
#include <GLFW/glfw3.h>
#include "utils.h"
#include "camera.h"
#include "vector.h"
#include "mesh.h"
int screen_w;
int screen_h;
float maxf(float x, float y) { return x > y ? x : y; }
float minf(float x, float y) { return x < y ? x : y; }
static unsigned int
compile_shader(const char *vertex_file,
const char *fragment_file)
{
int success;
char infolog[512];
char *vertex_str = load_file(vertex_file, NULL);
if (vertex_str == NULL) {
fprintf(stderr, "Couldn't load file '%s'\n", vertex_file);
return 0;
}
char *fragment_str = load_file(fragment_file, NULL);
if (fragment_str == NULL) {
fprintf(stderr, "Couldn't load file '%s'\n", fragment_file);
free(vertex_str);
return 0;
}
unsigned int vertex_shader = glCreateShader(GL_VERTEX_SHADER);
glShaderSource(vertex_shader, 1, &vertex_str, NULL);
glCompileShader(vertex_shader);
glGetShaderiv(vertex_shader, GL_COMPILE_STATUS, &success);
if(!success) {
glGetShaderInfoLog(vertex_shader, sizeof(infolog), NULL, infolog);
fprintf(stderr, "Couldn't compile vertex shader '%s' (%s)\n", vertex_file, infolog);
free(vertex_str);
free(fragment_str);
return 0;
}
unsigned int fragment_shader = glCreateShader(GL_FRAGMENT_SHADER);
glShaderSource(fragment_shader, 1, &fragment_str, NULL);
glCompileShader(fragment_shader);
glGetShaderiv(fragment_shader, GL_COMPILE_STATUS, &success);
if(!success) {
glGetShaderInfoLog(fragment_shader, sizeof(infolog), NULL, infolog);
fprintf(stderr, "Couldn't compile fragment shader '%s' (%s)\n", fragment_file, infolog);
free(vertex_str);
free(fragment_str);
return 0;
}
unsigned int shader_program = glCreateProgram();
glAttachShader(shader_program, vertex_shader);
glAttachShader(shader_program, fragment_shader);
glLinkProgram(shader_program);
glGetProgramiv(shader_program, GL_LINK_STATUS, &success);
if(!success) {
glGetProgramInfoLog(shader_program, sizeof(infolog), NULL, infolog);
fprintf(stderr, "Couldn't link shader program (%s)\n", infolog);
free(vertex_str);
free(fragment_str);
return 0;
}
glDeleteShader(vertex_shader);
glDeleteShader(fragment_shader);
free(vertex_str);
free(fragment_str);
return shader_program;
}
static void set_uniform_m4(unsigned int program, const char *name, Matrix4 value)
{
int location = glGetUniformLocation(program, name);
if (location < 0) {
printf("Can't set uniform '%s'\n", name);
abort();
}
glUniformMatrix4fv(location, 1, GL_FALSE, (float*) &value);
}
static void set_uniform_v3(unsigned int program, const char *name, Vector3 value)
{
int location = glGetUniformLocation(program, name);
if (location < 0) {
printf("Can't set uniform '%s' (program %d, location %d)\n", name, program, location);
abort();
}
glUniform3f(location, value.x, value.y, value.z);
}
static void set_uniform_i(unsigned int program, const char *name, int value)
{
int location = glGetUniformLocation(program, name);
if (location < 0) {
printf("Can't set uniform '%s'\n", name);
abort();
}
glUniform1i(location, value);
}
static void set_uniform_f(unsigned int program, const char *name, float value)
{
int location = glGetUniformLocation(program, name);
if (location < 0) {
printf("Can't set uniform '%s'\n", name);
abort();
}
glUniform1f(location, value);
}
static void error_callback(int error, const char* description)
{
fprintf(stderr, "Error: %s\n", description);
}
static void key_callback(GLFWwindow* window, int key, int scancode, int action, int mods)
{
if (key == GLFW_KEY_ESCAPE && action == GLFW_PRESS)
glfwSetWindowShouldClose(window, GLFW_TRUE);
}
void framebuffer_size_callback(GLFWwindow* window, int width, int height)
{
glViewport(0, 0, width, height);
}
void cursor_pos_callback(GLFWwindow *window, double x, double y)
{
rotate_camera(x, y);
}
typedef struct {
Vector3 origin;
Vector3 size;
} Cube;
bool intersect_cube(Ray r, Cube c, float *tnear, float *tfar, Vector3 *normal)
{
float txmin, txmax;
float tymin, tymax;
float tzmin, tzmax;
float tn;
float tf;
Vector3 a = c.origin;
Vector3 b = combine(c.origin, c.size, 1, 1);
int hit_axis = 0; // 0=x, 1=y, 2=z
if (r.direction.x >= 0) {
txmin = (a.x - r.origin.x) / r.direction.x;
txmax = (b.x - r.origin.x) / r.direction.x;
} else {
txmax = (a.x - r.origin.x) / r.direction.x;
txmin = (b.x - r.origin.x) / r.direction.x;
}
if (r.direction.y >= 0) {
tymin = (a.y - r.origin.y) / r.direction.y;
tymax = (b.y - r.origin.y) / r.direction.y;
} else {
tymax = (a.y - r.origin.y) / r.direction.y;
tymin = (b.y - r.origin.y) / r.direction.y;
}
if (txmin > tymax || tymin > txmax)
return false;
if (tymin > txmin) { txmin = tymin; hit_axis = 1; }
if (tymax < txmax) txmax = tymax;
if (r.direction.z >= 0) {
tzmin = (a.z - r.origin.z) / r.direction.z;
tzmax = (b.z - r.origin.z) / r.direction.z;
} else {
tzmax = (a.z - r.origin.z) / r.direction.z;
tzmin = (b.z - r.origin.z) / r.direction.z;
}
if (txmin > tzmax || tzmin > txmax)
return false;
if (tzmin > txmin) { txmin = tzmin; hit_axis = 2; };
if (tzmax < txmax) txmax = tzmax;
if (tnear) *tnear = txmin;
if (tfar) *tfar = txmax;
if (normal) {
switch (hit_axis) {
case 0: *normal = r.direction.x > 0 ? (Vector3) {-1, 0, 0} : (Vector3) {1, 0, 0}; break;
case 1: *normal = r.direction.y > 0 ? (Vector3) {0, -1, 0} : (Vector3) {0, 1, 0}; break;
case 2: *normal = r.direction.z > 0 ? (Vector3) {0, 0, -1} : (Vector3) {0, 0, 1}; break;
}
}
return true;
}
bool intersect_sphere(Ray r, Sphere s, float *t)
{
/*
* Any point of the ray can be written as
*
* P(t) = O + t * D
*
* with O origin and D direction.
*
* All points P=(x,y,z) of a sphere can be described as
* those (and only those) that satisfy the equation
*
* x^2 + y^2 + z^2 = R^2
* P^2 - R^2 = 0
*
* with R radius of the sphere. The sphere here is centered
* at the origin.
*
* Intersection points of the ray with the sphere must satisfy
* both:
*
* P(t) = O + t * D
* P^2 - R^2 = 0
*
* => (O + tD)^2 - R^2 = 0
* => t^2 * D^2 + t * 2OD + O^2 - R^2 = 0
*
* we can use the quadratic formula here, and more specifically
* the discriminant to check if solutions exist and how many
*/
Vector3 oc = combine(s.center, r.origin, 1, -1);
float a = dotv(r.direction, r.direction);
float b = -2 * dotv(oc, r.direction);
float c = dotv(oc, oc) - s.radius * s.radius;
float discr = b*b - 4*a*c;
if (discr > 0) {
float s0 = (- b + sqrt(discr)) / (2 * a);
float s1 = (- b - sqrt(discr)) / (2 * a);
if (s0 > s1) {
float tmp = s0;
s0 = s1;
s1 = tmp;
}
if (s0 < 0) {
s0 = s1;
if (s0 < 0) return false;
}
if (t) *t = s0;
return true;
}
// Zero solutions
return false;
}
typedef enum {
OBJECT_CUBE,
OBJECT_SPHERE,
} ObjectType;
typedef struct {
ObjectType type;
union {
Sphere sphere;
Cube cube;
};
Vector3 color;
} Object;
Object cube(Vector3 color, Vector3 origin, Vector3 size) { return (Object) {.color=color, .type=OBJECT_CUBE, .cube=(Cube) {.origin=origin, .size=size}}; }
Object sphere(Vector3 color, Vector3 origin, float radius) { return (Object) {.color=color, .type=OBJECT_SPHERE, .sphere=(Sphere) {.center=origin, .radius=radius}}; }
bool intersect_object(Ray r, Object o, float *t, Vector3 *normal)
{
switch (o.type) {
case OBJECT_CUBE:
return intersect_cube(r, o.cube, t, NULL, normal);
case OBJECT_SPHERE:
if (intersect_sphere(r, o.sphere, t)) {
if (normal) {
Vector3 hit_point = combine(r.origin, r.direction, 1, *t);
*normal = normalize(combine(hit_point, o.sphere.center, 1, -1));
}
return true;
}
return false;
}
return false;
}
float random_float(void)
{
return (float) rand() / RAND_MAX;
}
Vector3 random_vector(void)
{
return (Vector3) {
.x = random_float(),
.y = random_float(),
.z = random_float(),
};
}
Vector3 random_direction(void)
{
return normalize(random_vector());
}
Vector3 reflect(Vector3 dir, Vector3 normal)
{
float f = -2 * dotv(normal, dir);
return combine(dir, normal, 1, f);
}
#define MAX_OBJECTS 32
Object objects[MAX_OBJECTS];
int num_objects = 0;
void add_object(Object o)
{
if (num_objects < MAX_OBJECTS)
objects[num_objects++] = o;
}
typedef struct {
float distance;
Vector3 point;
Vector3 normal;
int object;
} HitInfo;
HitInfo trace_ray(Ray ray)
{
ray.direction = normalize(ray.direction);
float nearest_t = FLT_MAX;
int nearest_object = -1;
Vector3 nearest_normal;
for (int i = 0; i < num_objects; i++) {
float t;
Vector3 n;
if (!intersect_object(ray, objects[i], &t, &n))
continue;
if (t >= 0 && t < nearest_t) {
nearest_t = t;
nearest_object = i;
nearest_normal = n;
}
}
if (nearest_object == -1) {
HitInfo result;
result.distance = -1;
result.normal = (Vector3) {0, 0, 0};
result.point = (Vector3) {0, 0, 0};
result.object = -1;
return result;
} else {
HitInfo result;
result.distance = nearest_t;
result.normal = nearest_normal;
result.point = combine(ray.origin, ray.direction, 1, nearest_t);
result.object = nearest_object;
return result;
}
}
Vector3 pixel(float x, float y)
{
Ray ray = ray_through_screen_at(x, y, (float) screen_w/screen_h);
float multiplier = 1;
Vector3 color = {0, 0, 0};
int bounces = 2;
for (int i = 0; i < bounces; i++) {
HitInfo hit = trace_ray(ray);
if (hit.object == -1) {
Vector3 sky_color = {0, 0, 0};
color = combine(color, sky_color, 1, multiplier);
break;
}
Vector3 light_dir = normalize((Vector3) {-1, -1, -1});
#if 1
float light_intensity = 0;
if (trace_ray((Ray) {combine(hit.point, light_dir, 1, -0.001), scale(light_dir, -1)}).object == -1)
light_intensity = maxf(dotv(hit.normal, scale(light_dir, -1)), 0);
#elif 0
float light_intensity = maxf(dotv(hit.normal, scale(light_dir, -1)), 0);
#else
float light_intensity = 1;
#endif
color = combine(color, objects[hit.object].color, 1, light_intensity * multiplier);
multiplier *= 0.7;
Vector3 new_dir = reflect(ray.direction, scale(hit.normal, -1));
ray = (Ray) { combine(hit.point, new_dir, 1, 0.001), new_dir };
}
return color;
}
Vector3 *frame = NULL;
int frame_w = 0;
int frame_h = 0;
unsigned int frame_texture;
void update_frame_texture(float scale)
{
if (frame_w != scale * screen_w || frame_h != scale * screen_h) {
frame_w = scale * screen_w;
frame_h = scale * screen_h;
if (frame) free(frame);
frame = malloc(sizeof(Vector3) * frame_w * frame_h);
if (!frame) { printf("OUT OF MEMORY\n"); abort(); }
}
for (int j = 0; j < frame_h; j++)
for (int i = 0; i < frame_w; i++) {
float u = (float) i / (frame_w - 1);
float v = (float) j / (frame_h - 1);
u = 1 - u;
v = 1 - v;
frame[j * frame_w + i] = pixel(u, v);
}
glBindTexture(GL_TEXTURE_2D, frame_texture);
glTexImage2D(GL_TEXTURE_2D, 0, GL_RGB, frame_w, frame_h, 0, GL_RGB, GL_FLOAT, frame);
glBindTexture(GL_TEXTURE_2D, 0);
}
int main(void)
{
add_object(cube((Vector3) {0.6, 0.6, 0.8}, (Vector3) {0, 0, 0}, (Vector3) {10, 0.1, 10})),
add_object(cube((Vector3) {0.3, 0, 0}, (Vector3) {0, 0, 0}, (Vector3) {10, 0.1, 0.1})),
add_object(cube((Vector3) {0, 0.3, 0}, (Vector3) {0, 0, 0}, (Vector3) {0.1, 10, 0.1})),
add_object(cube((Vector3) {0, 0, 0.3}, (Vector3) {0, 0, 0}, (Vector3) {0.1, 0.1, 10})),
add_object(cube((Vector3) {0.3, 0, 0}, (Vector3) {7, 0, 8}, (Vector3) {1, 1, 1})),
add_object(cube((Vector3) {0.3, 0, 0.3}, (Vector3) {6, 0, 7}, (Vector3) {1, 1, 1})),
add_object(sphere((Vector3) {0.3, 0, 0}, (Vector3) {3, 1, 3}, 1)),
add_object(sphere((Vector3) {0, 0.3, 0}, (Vector3) {5, 2, 5}, 1)),
glfwSetErrorCallback(error_callback);
if (!glfwInit())
return -1;
glfwWindowHint(GLFW_CONTEXT_VERSION_MAJOR, 3);
glfwWindowHint(GLFW_CONTEXT_VERSION_MINOR, 3);
glfwWindowHint(GLFW_OPENGL_PROFILE, GLFW_OPENGL_CORE_PROFILE);
GLFWwindow *window = glfwCreateWindow(640, 480, "Path Trace", NULL, NULL);
if (!window) {
glfwTerminate();
return -1;
}
glfwGetWindowSize(window, &screen_w, &screen_h);
glfwSetKeyCallback(window, key_callback);
glfwSetFramebufferSizeCallback(window, framebuffer_size_callback);
glfwSetCursorPosCallback(window, cursor_pos_callback);
glfwSetInputMode(window, GLFW_CURSOR, GLFW_CURSOR_DISABLED);
glfwMakeContextCurrent(window);
if (!gladLoadGLLoader((GLADloadproc)glfwGetProcAddress)) {
printf("Failed to initialize GLAD\n");
return -1;
}
glfwSwapInterval(1);
unsigned int screen_program = compile_shader("assets/screen.vs", "assets/screen.fs");
if (!screen_program) { printf("Couldn't compile program\n"); return -1; }
set_uniform_i(screen_program, "screenTexture", 0);
unsigned int vao, vbo;
{
float vertices[] = {
// positions // texCoords
-1.0f, 1.0f, 0.0f, 1.0f,
-1.0f, -1.0f, 0.0f, 0.0f,
1.0f, -1.0f, 1.0f, 0.0f,
-1.0f, 1.0f, 0.0f, 1.0f,
1.0f, -1.0f, 1.0f, 0.0f,
1.0f, 1.0f, 1.0f, 1.0f
};
glGenVertexArrays(1, &vao);
glGenBuffers(1, &vbo);
glBindVertexArray(vao);
glBindBuffer(GL_ARRAY_BUFFER, vbo);
glBufferData(GL_ARRAY_BUFFER, sizeof(vertices), &vertices, GL_STATIC_DRAW);
glEnableVertexAttribArray(0);
glVertexAttribPointer(0, 2, GL_FLOAT, GL_FALSE, 4 * sizeof(float), (void*)0);
glEnableVertexAttribArray(1);
glVertexAttribPointer(1, 2, GL_FLOAT, GL_FALSE, 4 * sizeof(float), (void*)(2 * sizeof(float)));
}
glGenTextures(1, &frame_texture);
glBindTexture(GL_TEXTURE_2D, frame_texture);
glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_S, GL_CLAMP_TO_EDGE);
glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_T, GL_CLAMP_TO_EDGE);
glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER, GL_LINEAR);
glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MAG_FILTER, GL_LINEAR);
while (!glfwWindowShouldClose(window)) {
glfwGetWindowSize(window, &screen_w, &screen_h);
float speed = 0.5;
if (glfwGetKey(window, GLFW_KEY_W) == GLFW_PRESS) move_camera(UP, speed);
if (glfwGetKey(window, GLFW_KEY_S) == GLFW_PRESS) move_camera(DOWN, speed);
if (glfwGetKey(window, GLFW_KEY_A) == GLFW_PRESS) move_camera(LEFT, speed);
if (glfwGetKey(window, GLFW_KEY_D) == GLFW_PRESS) move_camera(RIGHT, speed);
Vector3 clear_color = {1, 1, 1};
update_frame_texture(0.4);
glViewport(0, 0, screen_w, screen_h);
glClearColor(clear_color.x, clear_color.y, clear_color.z, 1.0f);
glClearStencil(0);
glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT | GL_STENCIL_BUFFER_BIT);
glUseProgram(screen_program);
glActiveTexture(GL_TEXTURE0);
glBindTexture(GL_TEXTURE_2D, frame_texture);
glBindVertexArray(vao);
glDrawArrays(GL_TRIANGLES, 0, 6);
glBindVertexArray(0);
glfwSwapBuffers(window);
glfwPollEvents();
}
glfwDestroyWindow(window);
glfwTerminate();
return 0;
}
src/mesh.c
+345
View File
@@ -0,0 +1,345 @@
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include "utils.h"
#include "mesh.h"
#define TINYOBJ_LOADER_C_IMPLEMENTATION
#include "tinyobj_loader_c.h"
void append_vertex(VertexArray *array, Vertex v)
{
if (array->size == array->capacity) {
if (array->capacity == 0) {
array->data = malloc(8 * sizeof(Vertex));
array->capacity = 8;
} else {
array->data = realloc(array->data, 2 * array->capacity * sizeof(Vertex));
array->capacity *= 2;
}
if (!array->data) {
printf("OUT OF MEMORY\n");
abort();
}
}
array->data[array->size++] = v;
}
static Vector3 get_sphere_point(float angle_x, float angle_y, float radius)
{
Vector3 p;
p.x = radius * sin(angle_y) * cos(2 * angle_x);
p.y = radius * cos(angle_y);
p.z = radius * sin(angle_y) * sin(2 * angle_x);
return p;
}
typedef struct {
Vertex a;
Vertex b;
Vertex c;
} Triangle;
static void calculate_and_set_normals(Triangle *T)
{
#define X1 (T->b.x - T->a.x)
#define Y1 (T->b.y - T->a.y)
#define Z1 (T->b.z - T->a.z)
#define X2 (T->c.x - T->a.x)
#define Y2 (T->c.y - T->a.y)
#define Z2 (T->c.z - T->a.z)
/*
* xnormal = y1*z2 - z1*y2
* ynormal = z1*x2 - x1*z2
* znormal = x1*y2 - y1*x2
*/
Vector3 n;
n.x = Y1 * Z2 - Z1 * Y2;
n.y = Z1 * X2 - X1 * Z2;
n.z = X1 * Y2 - Y1 * X2;
T->a.nx = n.x;
T->a.ny = n.y;
T->a.nz = n.z;
T->b.nx = n.x;
T->b.ny = n.y;
T->b.nz = n.z;
T->c.nx = n.x;
T->c.ny = n.y;
T->c.nz = n.z;
#undef X1
#undef Y1
#undef Z1
#undef X2
#undef Y2
#undef Z2
}
static Triangle make_triangle(Vector3 a, Vector3 b, Vector3 c)
{
Triangle T;
T.a = (Vertex) {a.x, a.y, a.z};
T.b = (Vertex) {b.x, b.y, b.z};
T.c = (Vertex) {c.x, c.y, c.z};
calculate_and_set_normals(&T);
return T;
}
VertexArray make_sphere_mesh_2(float radius, int num_segms, bool fake_normals)
{
VertexArray vertices = {0, 0, 0};
int x_num_segms = num_segms;
int y_num_segms = num_segms;
for (int i = 0; i < x_num_segms; i++)
for (int j = 0; j < y_num_segms; j++) {
int g = j;
if (g == y_num_segms-1)
g = 0;
Vector3 p1 = get_sphere_point((i + 0) * 3.14 / x_num_segms, (g + 0) * 3.14 / y_num_segms, radius);
Vector3 p2 = get_sphere_point((i + 1) * 3.14 / x_num_segms, (g + 0) * 3.14 / y_num_segms, radius);
Vector3 p3 = get_sphere_point((i + 1) * 3.14 / x_num_segms, (g + 1) * 3.14 / y_num_segms, radius);
Vector3 p4 = get_sphere_point((i + 0) * 3.14 / x_num_segms, (g + 1) * 3.14 / y_num_segms, radius);
Triangle t1 = make_triangle(p1, p2, p3);
if (fake_normals) {
t1.a.nx = p1.x;
t1.a.ny = p1.y;
t1.a.nz = p1.z;
t1.b.nx = p2.x;
t1.b.ny = p2.y;
t1.b.nz = p2.z;
t1.c.nx = p3.x;
t1.c.ny = p3.y;
t1.c.nz = p3.z;
}
append_vertex(&vertices, t1.a);
append_vertex(&vertices, t1.b);
append_vertex(&vertices, t1.c);
Triangle t2 = make_triangle(p4, p1, p3);
if (fake_normals) {
t2.a.nx = p4.x;
t2.a.ny = p4.y;
t2.a.nz = p4.z;
t2.b.nx = p1.x;
t2.b.ny = p1.y;
t2.b.nz = p1.z;
t2.c.nx = p3.x;
t2.c.ny = p3.y;
t2.c.nz = p3.z;
}
append_vertex(&vertices, t2.a);
append_vertex(&vertices, t2.b);
append_vertex(&vertices, t2.c);
}
return vertices;
}
VertexArray make_sphere_mesh(float radius)
{
return make_sphere_mesh_2(radius, 32, true);
}
VertexArray make_cube_mesh(void)
{
float vertices[] = {
1.0, 1.0, 0.0, 0.0, 0.0, -1.0,
1.0, 0.0, 0.0, 0.0, 0.0, -1.0,
0.0, 0.0, 0.0, 0.0, 0.0, -1.0,
0.0, 1.0, 0.0, 0.0, 0.0, -1.0,
1.0, 1.0, 0.0, 0.0, 0.0, -1.0,
0.0, 0.0, 0.0, 0.0, 0.0, -1.0,
0.0, 0.0, 1.0, 0.0, 0.0, 1.0,
1.0, 0.0, 1.0, 0.0, 0.0, 1.0,
1.0, 1.0, 1.0, 0.0, 0.0, 1.0,
0.0, 0.0, 1.0, 0.0, 0.0, 1.0,
1.0, 1.0, 1.0, 0.0, 0.0, 1.0,
0.0, 1.0, 1.0, 0.0, 0.0, 1.0,
0.0, 0.0, 0.0, 0.0, -1.0, 0.0,
1.0, 0.0, 0.0, 0.0, -1.0, 0.0,
1.0, 0.0, 1.0, 0.0, -1.0, 0.0,
0.0, 0.0, 0.0, 0.0, -1.0, 0.0,
1.0, 0.0, 1.0, 0.0, -1.0, 0.0,
0.0, 0.0, 1.0, 0.0, -1.0, 0.0,
1.0, 1.0, 1.0, 0.0, 1.0, 0.0,
1.0, 1.0, 0.0, 0.0, 1.0, 0.0,
0.0, 1.0, 0.0, 0.0, 1.0, 0.0,
0.0, 1.0, 1.0, 0.0, 1.0, 0.0,
1.0, 1.0, 1.0, 0.0, 1.0, 0.0,
0.0, 1.0, 0.0, 0.0, 1.0, 0.0,
0.0, 1.0, 1.0, -1.0, 0.0, 0.0,
0.0, 1.0, 0.0, -1.0, 0.0, 0.0,
0.0, 0.0, 0.0, -1.0, 0.0, 0.0,
0.0, 0.0, 1.0, -1.0, 0.0, 0.0,
0.0, 1.0, 1.0, -1.0, 0.0, 0.0,
0.0, 0.0, 0.0, -1.0, 0.0, 0.0,
1.0, 0.0, 0.0, 1.0, 0.0, 0.0,
1.0, 1.0, 0.0, 1.0, 0.0, 0.0,
1.0, 1.0, 1.0, 1.0, 0.0, 0.0,
1.0, 0.0, 0.0, 1.0, 0.0, 0.0,
1.0, 1.0, 1.0, 1.0, 0.0, 0.0,
1.0, 0.0, 1.0, 1.0, 0.0, 0.0,
};
VertexArray result = {0, 0, 0};
for (int i = 0; i < (int) (sizeof(vertices)/sizeof(vertices[0])); i += 6) {
Vertex v;
v.x = vertices[i + 0];
v.y = vertices[i + 1];
v.z = vertices[i + 2];
v.nx = vertices[i + 3];
v.ny = vertices[i + 4];
v.nz = vertices[i + 5];
v.tx = 0;
v.ty = 0;
append_vertex(&result, v);
}
return result;
}
static void
get_file_data_callback(void *context, const char *filename,
const int is_mtl, const char *obj_filename, char **data, size_t *len)
{
(void) context;
if (filename == NULL) {
fprintf(stderr, "null filename\n");
(*data) = NULL;
(*len) = 0;
return;
}
size_t data_len = 0;
printf("Reading file '%s'\n", filename);
*data = load_file(filename, &data_len);
(*len) = data_len;
char **free_me = context;
if (*free_me == NULL) *free_me = *data;
}
bool load_mesh_from_file(const char *file, VertexArray *result)
{
tinyobj_attrib_t attrib;
tinyobj_shape_t* shapes = NULL;
size_t num_shapes;
tinyobj_material_t* materials = NULL;
size_t num_materials;
char *free_me = NULL;
unsigned int flags = TINYOBJ_FLAG_TRIANGULATE;
int ret = tinyobj_parse_obj(&attrib, &shapes, &num_shapes, &materials, &num_materials, file, get_file_data_callback, &free_me, flags);
if (ret != TINYOBJ_SUCCESS) {
if (free_me) free(free_me);
printf("Failed loading '%s'\n", file);
return false;
}
*result = (VertexArray) {0, 0, 0};
size_t face_offset = 0;
for (int i = 0; i < (int) attrib.num_face_num_verts; i++) {
assert(attrib.face_num_verts[i] % 3 == 0); /* assume all triangle faces. */
for (size_t f = 0; f < (size_t)attrib.face_num_verts[i] / 3; f++) {
tinyobj_vertex_index_t idx0 = attrib.faces[face_offset + f * 3 + 0];
tinyobj_vertex_index_t idx1 = attrib.faces[face_offset + f * 3 + 1];
tinyobj_vertex_index_t idx2 = attrib.faces[face_offset + f * 3 + 2];
Vertex v0;
Vertex v1;
Vertex v2;
/*
* Positions
*/
v0.x = attrib.vertices[idx0.v_idx * 3 + 0];
v0.y = attrib.vertices[idx0.v_idx * 3 + 1];
v0.z = attrib.vertices[idx0.v_idx * 3 + 2];
v1.x = attrib.vertices[idx1.v_idx * 3 + 0];
v1.y = attrib.vertices[idx1.v_idx * 3 + 1];
v1.z = attrib.vertices[idx1.v_idx * 3 + 2];
v2.x = attrib.vertices[idx2.v_idx * 3 + 0];
v2.y = attrib.vertices[idx2.v_idx * 3 + 1];
v2.z = attrib.vertices[idx2.v_idx * 3 + 2];
/*
* Normals
*/
v0.nx = attrib.normals[idx0.vn_idx * 3 + 0];
v0.ny = attrib.normals[idx0.vn_idx * 3 + 1];
v0.nz = attrib.normals[idx0.vn_idx * 3 + 2];
v1.nx = attrib.normals[idx1.vn_idx * 3 + 0];
v1.ny = attrib.normals[idx1.vn_idx * 3 + 1];
v1.nz = attrib.normals[idx1.vn_idx * 3 + 2];
v2.nx = attrib.normals[idx2.vn_idx * 3 + 0];
v2.ny = attrib.normals[idx2.vn_idx * 3 + 1];
v2.nz = attrib.normals[idx2.vn_idx * 3 + 2];
/*
* Texture coordinates
*/
v0.tx = attrib.normals[idx0.vt_idx * 2 + 0];
v0.ty = attrib.normals[idx0.vt_idx * 2 + 1];
v1.tx = attrib.normals[idx1.vt_idx * 2 + 0];
v1.ty = attrib.normals[idx1.vt_idx * 2 + 1];
v2.tx = attrib.normals[idx2.vt_idx * 2 + 0];
v2.ty = attrib.normals[idx2.vt_idx * 2 + 1];
append_vertex(result, v0);
append_vertex(result, v1);
append_vertex(result, v2);
}
face_offset += (size_t)attrib.face_num_verts[i];
}
if (free_me)
free(free_me);
tinyobj_attrib_free(&attrib);
tinyobj_shapes_free(shapes, num_shapes);
tinyobj_materials_free(materials, num_materials);
return true;
}
src/mesh.h
+21
View File
@@ -0,0 +1,21 @@
#include "vector.h"
typedef struct {
float x, y, z;
float nx, ny, nz;
float tx, ty;
} Vertex;
typedef struct {
Vertex *data;
int size;
int capacity;
} VertexArray;
void append_vertex(VertexArray *array, Vertex v);
VertexArray make_sphere_mesh(float radius);
VertexArray make_sphere_mesh_2(float radius, int num_segms, bool fake_normals);
VertexArray make_cube_mesh(void);
bool load_mesh_from_file(const char *file, VertexArray *result);
src/tinyobj_loader_c.h
+1735
View File
File diff suppressed because it is too large Load Diff
src/utils.c
+31
View File
@@ -0,0 +1,31 @@
#include <stdio.h>
#include <stdlib.h>
#include "utils.h"
char *load_file(const char *file, size_t *size)
{
FILE *stream = fopen(file, "rb");
if (stream == NULL) return NULL;
fseek(stream, 0, SEEK_END);
long size2 = ftell(stream);
fseek(stream, 0, SEEK_SET);
char *dst = (char*) malloc(size2+1);
if (dst == NULL) {
fclose(stream);
return NULL;
}
fread(dst, 1, size2, stream);
if (ferror(stream)) {
free(dst);
fclose(stream);
return NULL;
}
dst[size2] = '\0';
fclose(stream);
if (size) *size = size2;
return dst;
}
src/utils.h
+5
View File
@@ -0,0 +1,5 @@
#include <stddef.h>
#define COUNTOF(X) (int) (sizeof(X) / sizeof((X)[0]))
char *load_file(const char *file, size_t *size);
src/vector.c
+785
View File
@@ -0,0 +1,785 @@
#include <math.h>
#include <assert.h>
#include <stdio.h>
#include <string.h>
#include <stdlib.h>
#include "vector.h"
#define EPSILON 0.00001
float deg2rad(float deg)
{
return 3.14159265358979323846 * deg / 180;
}
float norm2_of(Vector3 v)
{
return v.x * v.x + v.y * v.y + v.z * v.z;
}
float norm_of(Vector3 v)
{
return sqrt(norm2_of(v));
}
Vector3 normalize(Vector3 v)
{
float norm = norm_of(v);
if (norm < EPSILON && norm > -EPSILON)
return v;
v.x /= norm;
v.y /= norm;
v.z /= norm;
return v;
}
Vector3 scale(Vector3 v, float f)
{
v.x *= f;
v.y *= f;
v.z *= f;
return v;
}
Vector3 combine(Vector3 u, Vector3 v, float a, float b)
{
Vector3 r;
r.x = u.x * a + v.x * b;
r.y = u.y * a + v.y * b;
r.z = u.z * a + v.z * b;
return r;
}
Vector3 combine4(Vector3 u, Vector3 v, Vector3 g, Vector3 t, float a, float b, float c, float d)
{
Vector3 r;
r.x = u.x * a + v.x * b + g.x * c + t.x * d;
r.y = u.y * a + v.y * b + g.y * c + t.y * d;
r.z = u.z * a + v.z * b + g.z * c + t.z * d;
return r;
}
Vector3 cross(Vector3 u, Vector3 v)
{
Vector3 r;
r.x = u.y * v.z - u.z * v.y;
r.y = u.z * v.x - u.x * v.z;
r.z = u.x * v.y - u.y * v.x;
return r;
}
Matrix4 translate_matrix(Vector3 v, float f)
{
Matrix4 m;
m.data[0][0] = 1;
m.data[0][1] = 0;
m.data[0][2] = 0;
m.data[0][3] = 0;
m.data[1][0] = 0;
m.data[1][1] = 1;
m.data[1][2] = 0;
m.data[1][3] = 0;
m.data[2][0] = 0;
m.data[2][1] = 0;
m.data[2][2] = 1;
m.data[2][3] = 0;
m.data[3][0] = f * v.x;
m.data[3][1] = f * v.y;
m.data[3][2] = f * v.z;
m.data[3][3] = 1;
return m;
}
Matrix4 scale_matrix(Vector3 v)
{
Matrix4 m;
m.data[0][0] = v.x;
m.data[0][1] = 0;
m.data[0][2] = 0;
m.data[0][3] = 0;
m.data[1][0] = 0;
m.data[1][1] = v.y;
m.data[1][2] = 0;
m.data[1][3] = 0;
m.data[2][0] = 0;
m.data[2][1] = 0;
m.data[2][2] = v.z;
m.data[2][3] = 0;
m.data[3][0] = 0;
m.data[3][1] = 0;
m.data[3][2] = 0;
m.data[3][3] = 1;
return m;
}
Matrix4 rotate_matrix_x(float angle)
{
Matrix4 m;
m.data[0][0] = 1;
m.data[1][0] = 0;
m.data[2][0] = 0;
m.data[3][0] = 0;
m.data[0][1] = 0;
m.data[1][1] = cos(angle);
m.data[2][1] = -sin(angle);
m.data[3][1] = 0;
m.data[0][2] = 0;
m.data[1][2] = sin(angle);
m.data[2][2] = cos(angle);
m.data[3][2] = 0;
m.data[0][3] = 0;
m.data[1][3] = 0;
m.data[2][3] = 0;
m.data[3][3] = 1;
return m;
}
Matrix4 rotate_matrix_y(float angle)
{
Matrix4 m;
m.data[0][0] = cos(angle);
m.data[1][0] = 0;
m.data[2][0] = sin(angle);
m.data[3][0] = 0;
m.data[0][1] = 0;
m.data[1][1] = 1;
m.data[2][1] = 0;
m.data[3][1] = 0;
m.data[0][2] = -sin(angle);
m.data[1][2] = 0;
m.data[2][2] = cos(angle);
m.data[3][2] = 0;
m.data[0][3] = 0;
m.data[1][3] = 0;
m.data[2][3] = 0;
m.data[3][3] = 1;
return m;
}
Matrix4 rotate_matrix_z(float angle)
{
Matrix4 m;
m.data[0][0] = cos(angle);
m.data[1][0] = -sin(angle);
m.data[2][0] = 0;
m.data[3][0] = 0;
m.data[0][1] = sin(angle);
m.data[1][1] = cos(angle);
m.data[2][1] = 0;
m.data[3][1] = 0;
m.data[0][2] = 0;
m.data[1][2] = 0;
m.data[2][2] = 1;
m.data[3][2] = 0;
m.data[0][3] = 0;
m.data[1][3] = 0;
m.data[2][3] = 0;
m.data[3][3] = 1;
return m;
}
Matrix4 transpose(Matrix4 m)
{
Matrix4 r;
r.data[0][0] = m.data[0][0];
r.data[0][1] = m.data[1][0];
r.data[0][2] = m.data[2][0];
r.data[0][3] = m.data[3][0];
r.data[1][0] = m.data[0][1];
r.data[1][1] = m.data[1][1];
r.data[1][2] = m.data[2][1];
r.data[1][3] = m.data[3][1];
r.data[2][0] = m.data[0][2];
r.data[2][1] = m.data[1][2];
r.data[2][2] = m.data[2][2];
r.data[2][3] = m.data[3][2];
r.data[3][0] = m.data[0][3];
r.data[3][1] = m.data[1][3];
r.data[3][2] = m.data[2][3];
r.data[3][3] = m.data[3][3];
return r;
}
Matrix4 identity_matrix(void)
{
Matrix4 m;
m.data[0][0] = 1;
m.data[0][1] = 0;
m.data[0][2] = 0;
m.data[0][3] = 0;
m.data[1][0] = 0;
m.data[1][1] = 1;
m.data[1][2] = 0;
m.data[1][3] = 0;
m.data[2][0] = 0;
m.data[2][1] = 0;
m.data[2][2] = 1;
m.data[2][3] = 0;
m.data[3][0] = 0;
m.data[3][1] = 0;
m.data[3][2] = 0;
m.data[3][3] = 1;
return m;
}
Vector4 ldotv(Vector4 v, Matrix4 m)
{
Vector4 r;
r.x = m.data[0][0] * v.x + m.data[0][1] * v.y + m.data[0][2] * v.z + m.data[0][3] * v.w;
r.y = m.data[1][0] * v.x + m.data[1][1] * v.y + m.data[1][2] * v.z + m.data[1][3] * v.w;
r.z = m.data[2][0] * v.x + m.data[2][1] * v.y + m.data[2][2] * v.z + m.data[2][3] * v.w;
r.w = m.data[3][0] * v.x + m.data[3][1] * v.y + m.data[3][2] * v.z + m.data[3][3] * v.w;
return r;
}
Vector4 rdotv(Matrix4 m, Vector4 v)
{
Vector4 r;
r.x = m.data[0][0] * v.x + m.data[1][0] * v.y + m.data[2][0] * v.z + m.data[3][0] * v.w;
r.y = m.data[0][1] * v.x + m.data[1][1] * v.y + m.data[2][1] * v.z + m.data[3][1] * v.w;
r.z = m.data[0][2] * v.x + m.data[1][2] * v.y + m.data[2][2] * v.z + m.data[3][2] * v.w;
r.w = m.data[0][3] * v.x + m.data[1][3] * v.y + m.data[2][3] * v.z + m.data[3][3] * v.w;
return r;
}
Matrix4 dotm(Matrix4 a, Matrix4 b)
{
Matrix4 r;
for (int i = 0; i < 4; i++)
for (int j = 0; j < 4; j++) {
r.data[i][j] = 0;
for (int k = 0; k < 4; k++)
r.data[i][j] += a.data[k][j] * b.data[i][k];
}
return r;
}
float dotv(Vector3 u, Vector3 v)
{
return u.x * v.x + u.y * v.y + u.z * v.z;
}
Vector3 mulv(Vector3 a, Vector3 b)
{
return (Vector3) {
.x = a.x * b.x,
.y = a.y * b.y,
.z = a.z * b.z,
};
}
Matrix4 lookat_matrix(Vector3 eye, Vector3 center, Vector3 up)
{
Vector3 forward = combine(center, eye, 1, -1);
forward = normalize(forward);
Vector3 right = cross(forward, up);
right = normalize(right);
up = cross(right, forward);
up = normalize(up);
Matrix4 m;
m.data[0][0] = right.x;
m.data[0][1] = up.x;
m.data[0][2] = -forward.x;
m.data[0][3] = 0;
m.data[1][0] = right.y;
m.data[1][1] = up.y;
m.data[1][2] = -forward.y;
m.data[1][3] = 0;
m.data[2][0] = right.z;
m.data[2][1] = up.z;
m.data[2][2] = -forward.z;
m.data[2][3] = 0;
m.data[3][0] = -dotv(right, eye);
m.data[3][1] = -dotv(up, eye);
m.data[3][2] = dotv(forward, eye);
m.data[3][3] = 1;
return m;
}
Matrix4 ortho_matrix(float left, float right, float bottom, float top, float near, float far)
{
Matrix4 m;
for (int i = 0; i < 4; i++)
for (int j = 0; j < 4; j++)
m.data[i][j] = 0;
m.data[0][0] = (float) 2 / (right - left);
m.data[1][1] = (float) 2 / (top - bottom);
m.data[2][2] = (float) -2 / (far - near);
m.data[3][3] = 1;
m.data[3][0] = -(right + left) / (right - left);
m.data[3][1] = -(top + bottom) / (top - bottom);
m.data[3][2] = -(far + near) / (far - near);
return m;
}
Matrix4 perspective_matrix(float fov, float aspect, float near, float far)
{
Matrix4 m;
for (int i = 0; i < 4; i++)
for (int j = 0; j < 4; j++)
m.data[i][j] = 0;
float tan_half_fov = tanf(fov / 2.0);
m.data[0][0] = (float) 1 / (aspect * tan_half_fov);
m.data[1][1] = (float) 1 / tan_half_fov;
m.data[2][2] = -(far + near) / (far - near);
m.data[2][3] = -1;
m.data[3][2] = -(2 * far * near) / (far - near);
return m;
}
void print_matrix(Matrix4 m)
{
printf(
"| %f %f %f %f |\n"
"| %f %f %f %f |\n"
"| %f %f %f %f |\n"
"| %f %f %f %f |\n\n",
m.data[0][0],
m.data[0][1],
m.data[0][2],
m.data[0][3],
m.data[1][0],
m.data[1][1],
m.data[1][2],
m.data[1][3],
m.data[2][0],
m.data[2][1],
m.data[2][2],
m.data[2][3],
m.data[3][0],
m.data[3][1],
m.data[3][2],
m.data[3][3]);
}
void print_matrix3(Matrix3 m)
{
printf(
"| %f %f %f |\n"
"| %f %f %f |\n"
"| %f %f %f |\n\n",
m.data[0][0],
m.data[0][1],
m.data[0][2],
m.data[1][0],
m.data[1][1],
m.data[1][2],
m.data[2][0],
m.data[2][1],
m.data[2][2]);
}
/*
void lina_transpose(float *A, float *B, int m, int n)
{
assert(m > 0 && n > 0);
assert(A != NULL && B != NULL);
if(m == 1 || n == 1) {
// For a matrix with height or width of 1
// row-major and column-major order coincide,
// so the stransposition doesn't change the
// the memory representation. A simple copy
// does the job.
if(A != B) // Does the copy or the branch cost more?
memcpy(B, A, sizeof(A[0]) * m * n);
} else if(m == n) {
// Iterate over the upper triangular portion of
// the matrix and switch each element with the
// corresponding one in the lower triangular portion.
// NOTE: We're assuming A,B might be the same matrix.
// If A,B are the same matrix, then the diagonal
// is copied onto itself. By removing the +1 in
// the inner loop, the copying of the diagonal
// is avoided.
for(int i = 0; i < n; i += 1)
for(int j = 0; j < i+1; j += 1) {
float temp = A[i*n + j];
B[i*n + j] = A[j*n + i];
B[j*n + i] = temp;
}
} else {
// Not only the matrix needs to be transposed
// assuming the destination matrix is the same
// as the source matrix, but the memory representation
// of the matrix needs to switch from row-major
// to col-major, so it's not as simple as switching
// value's positions.
// This algorithm starts from the A[0][1] value and
// moves it where it needs to go, then gets the value
// that was at that position and puts that in it's
// new position. This process is iterated until the
// starting point A[0][1] is overwritten with the
// new value. In this process the first and last
// value of the matrix never move.
B[0] = A[0];
B[m*n - 1] = A[m*n - 1];
float item = A[1];
int next = m;
while(next != 1) {
float temp = A[next];
B[next] = item;
item = temp;
next = (next % n) * m + (next / n);
}
B[1] = item;
}
}
int lina_decompLUP(float *A, float *L,
float *U, int *P,
int n)
{
assert(n > 0);
assert(A != L && A != U && L != U);
for (int i = 0; i < n; i++)
P[i] = i;
int swaps = 0;
for (int i = 0; i < n; i++) {
int v = P[i];
float max_v = A[v * n + i];
int max_i = i;
for (int j = i+1; j < n; j++) {
int u = P[j];
float abs = fabs(A[u * n + j]);
if (abs > max_v) {
max_v = abs;
max_i = j;
}
}
if (max_i != i) {
// Swap rows
int temp = P[i];
P[i] = P[max_i];
P[max_i] = temp;
swaps++;
}
}
for (int i = 0; i < n; i++)
for (int j = 0; j < n; j++)
U[i * n + j] = A[P[i] * n + j];
memset(L, 0, sizeof(float) * n * n);
for (int i = 0; i < n; i++)
L[i * n + i] = 1;
for (int i = 0; i < n; i++)
for (int j = i+1; j < n; j++) {
float u = U[i * n + i];
L[j * n + i] = U[j * n + i] / u;
for (int k = 0; k < n; k++)
U[j * n + k] -= L[j * n + i] * U[i * n + k];
}
return swaps;
}
// Function: lina_det
//
// Calculates the determinant of the n by n matrix A
// and returns it throught the output parameter [det].
//
// If not enough memory is available, false is returned,
// else true is returned.
//
// Notes:
// - The output parameter [det] is optional. (you can
// ignore the result by passing NULL).
//
bool lina_det(float *A, int n, float *det)
{
// Allocate the space for the L,U matrices.
// I can't think of a version of this algorithm
// where a temporary buffer isn't necessary.
float *T = (float*) malloc(sizeof(float) * n * n * 2 + sizeof(int) * n);
if (T == NULL)
return false;
// Do the decomposition
float *L = T;
float *U = L + (n * n);
int *P = (int*) (U + (n * n));
int swaps = lina_decompLUP(A, L, U, P, n);
if (swaps < 0) {
free(T);
return false;
}
// Knowing that
//
// A = LU
//
// then
//
// det(A) = det(LU) = det(L)det(U)
//
// Since L and U are triangular, their
// determinant is the product of their
// diagonals, so the product of the
// determinants is the product of both
// the diagonals.
float prod = 1;
for (int i = 0; i < n; i++) {
float l = L[i * n + i];
float u = U[i * n + i];
prod *= l * u;
}
if (swaps & 1)
prod = -prod;
if (det)
*det = prod;
free(T);
return true;
}
// Create the n-1 by n-1 matrix D obtained by
// removing the [del_col] column and [del_row]
// frow the n by n matrix M.
static void
copyMatrixWithoutRowAndCol(float *M, float *D, int n,
int del_col, int del_row)
{
// Copy the upper-left portion of matrix M
// that comes before the deleted column and
// row.
for (int i = 0; i < del_row; i++)
for (int j = 0; j < del_col; j++)
D[i * (n-1) + j] = M[i * n + j];
// Copy the lower left portion that comes
// after both the deleted column and row.
for (int i = del_row+1; i < n; i++)
for (int j = del_col+1; j < n; j++)
D[(i-1) * (n-1) + (j-1)] = M[i * n + j];
// Copy the bottom portion that comes after
// the deleted row but before the deleted column.
for (int i = del_row+1; i < n; i++)
for (int j = 0; j < del_col; j++)
D[(i-1) * (n-1) + j] = M[i * n + j];
// Copy the right portion that comes after
// the deleted column but before the deleted row.
for (int i = 0; i < del_row; i++)
for (int j = del_col+1; j < n; j++)
D[i * (n-1) + (j-1)] = M[i * n + j];
}
bool lina_inverse(Matrix4 M, Matrix4 *D)
{
float det;
if (!lina_det((float*) &M, 4, &det))
return false;
printf("det=%f\n", det);
if (det == 0)
return false; // The matrix can't be inverted
Matrix3 T;
Matrix4 M_t = transpose(M);
for (int i = 0; i < 4; i++)
for (int j = 0; j < 4; j++) {
copyMatrixWithoutRowAndCol((float*) &M_t, (float*)&T, 4, j, i);
float det2;
if (!lina_det((float*) &T, 3, &det2))
return false;
printf("-----------------\n");
print_matrix3(T);
printf("det2=%f\n", det2);
// If the determinant of M isn't zero,
// neither is this!
assert(det2 != 0);
bool i_is_odd = i & 1;
bool j_is_odd = j & 1;
int sign = (i_is_odd == j_is_odd) ? 1 : -1;
D->data[i][j] = sign * det2 / det;
}
return true;
}
Matrix4 invert(Matrix4 m)
{
Matrix4 r;
if (!lina_inverse(m, &r)) {
printf("Couldn't invert!\n");
abort();
}
return r;
}
*/
static bool gluInvertMatrix(const float m[16], float invOut[16])
{
float inv[16], det;
int i;
inv[0] =
m[5] * m[10] * m[15] -
m[5] * m[11] * m[14] -
m[9] * m[6] * m[15] +
m[9] * m[7] * m[14] +
m[13] * m[6] * m[11] -
m[13] * m[7] * m[10];
inv[4] = -m[4] * m[10] * m[15] +
m[4] * m[11] * m[14] +
m[8] * m[6] * m[15] -
m[8] * m[7] * m[14] -
m[12] * m[6] * m[11] +
m[12] * m[7] * m[10];
inv[8] = m[4] * m[9] * m[15] -
m[4] * m[11] * m[13] -
m[8] * m[5] * m[15] +
m[8] * m[7] * m[13] +
m[12] * m[5] * m[11] -
m[12] * m[7] * m[9];
inv[12] = -m[4] * m[9] * m[14] +
m[4] * m[10] * m[13] +
m[8] * m[5] * m[14] -
m[8] * m[6] * m[13] -
m[12] * m[5] * m[10] +
m[12] * m[6] * m[9];
inv[1] = -m[1] * m[10] * m[15] +
m[1] * m[11] * m[14] +
m[9] * m[2] * m[15] -
m[9] * m[3] * m[14] -
m[13] * m[2] * m[11] +
m[13] * m[3] * m[10];
inv[5] = m[0] * m[10] * m[15] -
m[0] * m[11] * m[14] -
m[8] * m[2] * m[15] +
m[8] * m[3] * m[14] +
m[12] * m[2] * m[11] -
m[12] * m[3] * m[10];
inv[9] = -m[0] * m[9] * m[15] +
m[0] * m[11] * m[13] +
m[8] * m[1] * m[15] -
m[8] * m[3] * m[13] -
m[12] * m[1] * m[11] +
m[12] * m[3] * m[9];
inv[13] = m[0] * m[9] * m[14] -
m[0] * m[10] * m[13] -
m[8] * m[1] * m[14] +
m[8] * m[2] * m[13] +
m[12] * m[1] * m[10] -
m[12] * m[2] * m[9];
inv[2] = m[1] * m[6] * m[15] -
m[1] * m[7] * m[14] -
m[5] * m[2] * m[15] +
m[5] * m[3] * m[14] +
m[13] * m[2] * m[7] -
m[13] * m[3] * m[6];
inv[6] = -m[0] * m[6] * m[15] +
m[0] * m[7] * m[14] +
m[4] * m[2] * m[15] -
m[4] * m[3] * m[14] -
m[12] * m[2] * m[7] +
m[12] * m[3] * m[6];
inv[10] = m[0] * m[5] * m[15] -
m[0] * m[7] * m[13] -
m[4] * m[1] * m[15] +
m[4] * m[3] * m[13] +
m[12] * m[1] * m[7] -
m[12] * m[3] * m[5];
inv[14] = -m[0] * m[5] * m[14] +
m[0] * m[6] * m[13] +
m[4] * m[1] * m[14] -
m[4] * m[2] * m[13] -
m[12] * m[1] * m[6] +
m[12] * m[2] * m[5];
inv[3] = -m[1] * m[6] * m[11] +
m[1] * m[7] * m[10] +
m[5] * m[2] * m[11] -
m[5] * m[3] * m[10] -
m[9] * m[2] * m[7] +
m[9] * m[3] * m[6];
inv[7] = m[0] * m[6] * m[11] -
m[0] * m[7] * m[10] -
m[4] * m[2] * m[11] +
m[4] * m[3] * m[10] +
m[8] * m[2] * m[7] -
m[8] * m[3] * m[6];
inv[11] = -m[0] * m[5] * m[11] +
m[0] * m[7] * m[9] +
m[4] * m[1] * m[11] -
m[4] * m[3] * m[9] -
m[8] * m[1] * m[7] +
m[8] * m[3] * m[5];
inv[15] = m[0] * m[5] * m[10] -
m[0] * m[6] * m[9] -
m[4] * m[1] * m[10] +
m[4] * m[2] * m[9] +
m[8] * m[1] * m[6] -
m[8] * m[2] * m[5];
det = m[0] * inv[0] + m[1] * inv[4] + m[2] * inv[8] + m[3] * inv[12];
if (det == 0)
return false;
det = 1.0 / det;
for (i = 0; i < 16; i++)
invOut[i] = inv[i] * det;
return true;
}
bool invert(Matrix4 a, Matrix4 *b)
{
return gluInvertMatrix((float*) &a, (float*) b);
}
src/vector.h
+68
View File
@@ -0,0 +1,68 @@
#ifndef VECTOR_INCLUDED
#define VECTOR_INCLUDED
#include <stdbool.h>
typedef struct {
float x;
float y;
float z;
} Vector3;
typedef struct {
float x;
float y;
float z;
float w;
} Vector4;
typedef struct {
float data[4][4];
} Matrix4;
typedef struct {
float data[3][3];
} Matrix3;
typedef struct {
Vector3 origin;
Vector3 direction;
} Ray;
typedef struct {
Vector3 center;
float radius;
} Sphere;
float deg2rad(float deg);
Matrix4 translate_matrix(Vector3 v, float f);
Matrix4 identity_matrix(void);
Matrix4 scale_matrix(Vector3 v);
Matrix4 rotate_matrix_x(float angle);
Matrix4 rotate_matrix_y(float angle);
Matrix4 rotate_matrix_z(float angle);
Matrix4 lookat_matrix(Vector3 pos, Vector3 front, Vector3 up);
Matrix4 ortho_matrix(float left, float right, float bottom, float top, float near, float far);
Matrix4 perspective_matrix(float fov, float aspect, float near, float far);
void print_matrix(Matrix4 m);
float norm2_of(Vector3 v);
float norm_of(Vector3 v);
Vector3 normalize(Vector3 v);
Vector3 scale(Vector3 v, float f);
Vector3 combine(Vector3 u, Vector3 v, float a, float b);
Vector3 combine4(Vector3 u, Vector3 v, Vector3 g, Vector3 t, float a, float b, float c, float d);
Vector3 cross(Vector3 u, Vector3 v);
Vector3 mulv(Vector3 a, Vector3 b);
Vector4 ldotv(Vector4 v, Matrix4 m);
Vector4 rdotv(Matrix4 m, Vector4 v);
float dotv(Vector3 u, Vector3 v);
Matrix4 dotm(Matrix4 a, Matrix4 b);
Matrix4 transpose(Matrix4 m);
bool invert(Matrix4 a, Matrix4 *inv);
#endif