commit 067b8b25158e91bbbcf8650c8b61411ddfb86e87
parent ac2e2edc450ebfa639408b5a889365359eaac1fc
Author: Egor Achkasov <eaachkasov@edu.hse.ru>
Date: Wed, 26 Mar 2025 20:29:33 +0100
Major refactor
Diffstat:
22 files changed, 2051 insertions(+), 1997 deletions(-)
diff --git a/.gitignore b/.gitignore
@@ -1,4 +1,3 @@
-test
scenes
.vscode
build
diff --git a/Makefile b/Makefile
@@ -0,0 +1,36 @@
+CC := gcc
+CFLAGS := -Wall -Wextra -Wno-builtin-declaration-mismatch
+LDFLAGS := -lm
+
+SRC := src/materials.c src/scene.c src/compute_paths.c test/test.c
+VIZRAYS_SRC := src/scene.c src/vizrays.c
+
+TARGET := test.elf
+VIZRAYS_TARGET := vizrays.elf
+
+.PHONY: all dbg clean vizrays vizraysdbg
+
+all: CFLAGS += -O3
+all: $(TARGET)
+
+dbg: CFLAGS += -g -O0
+dbg: $(TARGET)
+
+vizrays: CFLAGS += -O3
+vizrays: LDFLAGS += -lGL -lGLU -lglut
+vizrays: $(VIZRAYS_TARGET)
+
+vizraysdbg: CFLAGS += -g -O0
+vizraysdbg: LDFLAGS += -lGL -lGLU -lglut
+vizraysdbg: $(VIZRAYS_TARGET)
+
+$(TARGET): $(SRC)
+ $(CC) $(CFLAGS) $^ $(LDFLAGS) -o $@
+
+$(VIZRAYS_TARGET): $(VIZRAYS_SRC)
+ $(CC) $(CFLAGS) $^ $(LDFLAGS) -o $@
+
+clean:
+ rm -f $(TARGET) $(VIZRAYS_TARGET) rt rt.*.so *.bin a.out vizrays
+ rm -rf __pycache__
+
diff --git a/compute_paths.c b/compute_paths.c
@@ -1,884 +0,0 @@
-#include "inc/compute_paths.h" /* for compute_paths */
-#include "inc/scene.h" /* for HRT_Scene, HRT_mesh, HRT_Material */
-#include "inc/materials.h" /* for g_hrt_materials */
-
-#include <stddef.h> /* for size_t */
-#include <stdlib.h> /* for exit, malloc, free */
-#include <string.h> /* for strlen, sprintf */
-#include <stdio.h> /* for fopen, FILE, fclose */
-#include <math.h> /* for sin, cos, sqrt */
-#include <stdint.h> /* for int32_t, uint8_t */
-
-/* ==== CONSTANTS ==== */
-
-#define PI 3.14159265358979323846f /* pi */
-#define SPEED_OF_LIGHT 299792458.0f /* m/s */
-
-/* (2w, w) binomial coefficients for w = 0, 1, ..., 19 */
-const float BINOMIAL_2W_W[20] = {
- 1.f, 2.f, 6.f, 20.f, 70.f,
- 252.f, 924.f, 3432.f, 12870.f, 48620.f,
- 184756.f, 705432.f, 2704156.f, 10400600.f, 40116600.f,
- 155117520.f, 601080390.f, 2333606220.f, 9075135300.f, 35345263800.f
-};
-/* (aplha, k) binomial coefficients for alpha = 0, 1, ..., 19 with k = 0, 1, ..., alpha */
-const float* BINOMIAL_ALPHA_K[20] = {
- (float[]){
- 1.f
- },
- (float[]){
- 1.f, 1.f
- },
- (float[]){
- 1.f, 2.f, 1.f
- },
- (float[]){
- 1.f, 3.f, 3.f, 1.f
- },
- (float[]){
- 1.f, 4.f, 6.f, 4.f, 1.f
- },
- (float[]){
- 1.f, 5.f, 10.f, 10.f, 5.f,
- 1.f
- },
- (float[]){
- 1.f, 6.f, 15.f, 20.f, 15.f,
- 6.f, 1.f
- },
- (float[]){
- 1.f, 7.f, 21.f, 35.f, 35.f,
- 21.f, 7.f, 1.f
- },
- (float[]){
- 1.f, 8.f, 28.f, 56.f, 70.f,
- 56.f, 28.f, 8.f, 1.f
- },
- (float[]){
- 1.f, 9.f, 36.f, 84.f, 126.f,
- 126.f, 84.f, 36.f, 9.f, 1.f
- },
- (float[]){
- 1.f, 10.f, 45.f, 120.f, 210.f,
- 252.f, 210.f, 120.f, 45.f, 10.f,
- 1.f
- },
- (float[]){
- 1.f, 11.f, 55.f, 165.f, 330.f,
- 462.f, 462.f, 330.f, 165.f, 55.f,
- 11.f, 1.f
- },
- (float[]){
- 1.f, 12.f, 66.f, 220.f, 495.f,
- 792.f, 924.f, 792.f, 495.f, 220.f,
- 66.f, 12.f, 1.f
- },
- (float[]){
- 1.f, 13.f, 78.f, 286.f,
- 715.f, 1287.f, 1716.f, 1716.f, 1287.f,
- 715.f, 286.f, 78.f, 13.f, 1.f
- },
- (float[]){
- 1.f, 14.f, 91.f, 364.f, 1001.f,
- 2002.f, 3003.f, 3432.f, 3003.f, 2002.f,
- 1001.f, 364.f, 91.f, 14.f, 1.f
- },
- (float[]){
- 1.f, 15.f, 105.f, 455.f, 1365.f,
- 3003.f, 5005.f, 6435.f, 6435.f, 5005.f,
- 3003.f, 1365.f, 455.f, 105.f, 15.f,
- 1.f
- },
- (float[]){
- 1.f, 16.f, 120.f, 560.f, 1820.f,
- 4368.f, 8008.f, 11440.f, 12870.f, 11440.f,
- 8008.f, 4368.f, 1820.f, 560.f, 120.f,
- 16.f, 1.f
- },
- (float[]){
- 1.f, 17.f, 136.f, 680.f, 2380.f,
- 6188.f, 12376.f, 19448.f, 24310.f, 24310.f,
- 19448.f, 12376.f, 6188.f, 2380.f, 680.f,
- 136.f, 17.f, 1.f
- },
- (float[]){
- 1.f, 18.f, 153.f, 816.f, 3060.f,
- 8568.f, 18564.f, 31824.f, 43758.f, 48620.f,
- 43758.f, 31824.f, 18564.f, 8568.f, 3060.f,
- 816.f, 153.f, 18.f, 1.f
- },
- (float[]){
- 1.f, 19.f, 171.f, 969.f, 3876.f,
- 11628.f, 27132.f, 50388.f, 75582.f, 92378.f,
- 92378.f, 75582.f, 50388.f, 27132.f, 11628.f,
- 3876.f, 969.f, 171.f, 19.f, 1.f
- }
-};
-
-/* ==== STRUCTS ==== */
-
-typedef struct {
- float x, y, z;
-} Vec3;
-
-typedef struct {
- Vec3 o, d;
-} Ray;
-
-/* Radio material eta and additional parameters */
-typedef struct {
- /* eta */
- float eta_re, eta_sqrt_re, eta_inv_re, eta_inv_sqrt_re;
- float eta_im, eta_sqrt_im, eta_inv_im, eta_inv_sqrt_im;
- float eta_abs, eta_abs_pow2, eta_abs_inv_sqrt;
- /* r = 1.0 - s */
- float r;
-} Material;
-
-typedef struct {
- uint32_t num_vertices;
- Vec3 *vs;
- uint32_t num_triangles;
- uint32_t *is;
- /* Normals. Size [num_triangles] */
- Vec3 *ns;
- uint32_t material_index;
- float velocity[3];
-} Mesh;
-
-typedef struct {
- uint32_t num_meshes;
- Mesh *meshes;
-} Scene;
-
-/* ==== PRECOMPUTED GLOBALS ==== */
-
-/* 4.f * PI * carrier_frequency * 1e9 / SPEED_OF_LIGHT */
-float g_free_space_loss_multiplier;
-
-/* Materials etas, calucalted for the given carrier_frequency */
-Material g_materials[NUM_G_MATERIALS];
-
-/* ==== COMPLEX OPERATIONS ==== */
-
-void csqrtf(
- IN float re,
- IN float im,
- IN float abs,
- OUT float *sqrt_re,
- OUT float *sqrt_im
-)
-{
- *sqrt_re = sqrtf((re + abs) / 2.f);
- if (fabsf(im) < __FLT_EPSILON__ && re >= -__FLT_EPSILON__)
- *sqrt_im = 0;
- else {
- *sqrt_im = sqrtf((abs - re) / 2.f);
- if (im < 0.f) *sqrt_im = -*sqrt_im;
- }
-}
-void cdivf(
- IN float a_re,
- IN float a_im,
- IN float b_re,
- IN float b_im,
- OUT float *c_re,
- OUT float *c_im
-)
-{
- float d = b_re * b_re + b_im * b_im;
- *c_re = (a_re * b_re + a_im * b_im) / d;
- *c_im = (a_im * b_re - a_re * b_im) / d;
-}
-
-/* ==== VECTOR OPERATIONS ==== */
-
-Vec3 vec3_sub(const Vec3 *a, const Vec3 *b)
-{
- return (Vec3){a->x - b->x, a->y - b->y, a->z - b->z};
-}
-Vec3 vec3_add(const Vec3 *a, const Vec3 *b)
-{
- return (Vec3){a->x + b->x, a->y + b->y, a->z + b->z};
-}
-Vec3 vec3_cross(const Vec3 *a, const Vec3 *b)
-{
- return (Vec3){
- a->y * b->z - a->z * b->y,
- a->z * b->x - a->x * b->z,
- a->x * b->y - a->y * b->x
- };
-}
-float vec3_dot(const Vec3 *a, const Vec3 *b)
-{
- return a->x * b->x + a->y * b->y + a->z * b->z;
-}
-Vec3 vec3_scale(const Vec3 *a, float s)
-{
- return (Vec3){a->x * s, a->y * s, a->z * s};
-}
-Vec3 vec3_normalize(const Vec3 *a)
-{
- float norm = sqrtf(a->x * a->x + a->y * a->y + a->z * a->z);
- return (Vec3){a->x / norm, a->y / norm, a->z / norm};
-}
-
-/* ==== SCENE LOADING ==== */
-
-/** Load a mesh from a HRT file.
- *
- * \param scene_filepath path to the HRT file
- * \param carrier_frequency the carrier frequency in GHz
- * \return the loaded scene
- */
-Scene load_scene(const char *scene_filepath, float carrier_frequency)
-{
- /* Open the HRT file */
- FILE *f = fopen(scene_filepath, "rb");
- if (!f) {
- fprintf(stderr, "Could not open file %s\n", scene_filepath);
- exit(8);
- }
-
- /* Parse the file */
- /* MAGIC */
- char magic[3];
- if (fread(magic, 1, 3, f) != 3) exit(8);
- if (strncmp(magic, "HRT", 3)) exit(8);
- /* SCENE */
- Scene scene;
- /* num_meshes */
- if (fread(&scene.num_meshes, sizeof(uint32_t), 1, f) != 1) exit(8);
- /* meshes */
- scene.meshes = (Mesh*)malloc(scene.num_meshes * sizeof(Mesh));
- for (uint32_t i = 0; i != scene.num_meshes; ++i) {
- /* MESH */
- Mesh *mesh = &scene.meshes[i];
- /* num_vertices */
- if(fread(&mesh->num_vertices, sizeof(uint32_t), 1, f) != 1) exit(8);
- /* vertices */
- mesh->vs = (Vec3*)malloc(mesh->num_vertices * sizeof(Vec3));
- if (fread(mesh->vs, sizeof(Vec3), mesh->num_vertices, f) != mesh->num_vertices)
- exit(8);
- /* num_triangles */
- if (fread(&mesh->num_triangles, sizeof(uint32_t), 1, f) != 1) exit(8);
- /* triangles */
- mesh->is = (uint32_t*)malloc(mesh->num_triangles * 3 * sizeof(uint32_t));
- if (fread(mesh->is, sizeof(uint32_t), mesh->num_triangles * 3, f) != mesh->num_triangles * 3)
- exit(8);
- /* material_index */
- if (fread(&mesh->material_index, sizeof(uint32_t), 1, f) != 1) exit(8);
- /* velocity */
- if (fread(&mesh->velocity, sizeof(float), 3, f) != 3) exit(8);
- }
- fclose(f);
-
- /* Precompute material permitivity */
- uint8_t material_eta_isdone[NUM_G_MATERIALS] = {0};
- for (uint32_t i = 0; i != scene.num_meshes; ++i)
- if (!material_eta_isdone[scene.meshes[i].material_index]) {
- HRT_Material *hrt_mat = &g_hrt_materials[scene.meshes[i].material_index];
- Material *rm = &g_materials[scene.meshes[i].material_index];
- /* calculate eta */
- rm->eta_re = hrt_mat->a * powf(carrier_frequency, hrt_mat->b);
- /* eq. 12 */
- rm->eta_im = (hrt_mat->c * powf(carrier_frequency, hrt_mat->d))
- / (0.0556325027352135f * carrier_frequency);
- rm->eta_abs_pow2 = rm->eta_re * rm->eta_re + rm->eta_im * rm->eta_im;
- rm->eta_abs = sqrtf(rm->eta_abs_pow2);
- rm->eta_abs_inv_sqrt = 1.f / sqrtf(rm->eta_abs);
- /* calculate sqrt(eta) */
- csqrtf(rm->eta_re, rm->eta_im, rm->eta_abs, &rm->eta_sqrt_re, &rm->eta_sqrt_im);
- /* calculate 1 / eta */
- rm->eta_inv_re = rm->eta_re / rm->eta_abs_pow2;
- rm->eta_inv_im = -rm->eta_im / rm->eta_abs_pow2;
- /* calculate 1 / sqrt(eta) */
- csqrtf(rm->eta_inv_re, rm->eta_inv_im, 1.f / rm->eta_abs,
- &rm->eta_inv_sqrt_re, &rm->eta_inv_sqrt_im);
- /* calculate r */
- rm->r = 1.f - hrt_mat->s;
- }
-
- /* Calculate normals */
- for (uint32_t i = 0; i != scene.num_meshes; ++i) {
- Mesh *mesh = &scene.meshes[i];
- mesh->ns = (Vec3*)malloc(mesh->num_triangles * sizeof(Vec3));
- for (uint32_t j = 0; j != mesh->num_triangles; ++j) {
- Vec3 v1 = mesh->vs[mesh->is[j * 3]];
- Vec3 v2 = mesh->vs[mesh->is[j * 3 + 1]];
- Vec3 v3 = mesh->vs[mesh->is[j * 3 + 2]];
- Vec3 e1 = vec3_sub(&v2, &v1);
- Vec3 e2 = vec3_sub(&v3, &v1);
- mesh->ns[j] = vec3_cross(&e1, &e2);
- mesh->ns[j] = vec3_normalize(&mesh->ns[j]);
- }
- }
-
- return scene;
-}
-
-/* ==== RT ==== */
-
-/** Compute Moeleer-Trumbore intersection algorithm.
- *
- * \param ray ray to cast to the mesh
- * \param scene the scene to cast the ray to
- * \param t output distance to the hit point. If no hit, t is not modified.
- * \param mesh_ind output index of the hit mesh. If no hit, i is not modified.
- * \param face_ind output index of the hit triangle. If no hit, i is not modified.
- * \param theta output angle of incidence
- */
-void moeller_trumbore(
- IN Ray *ray,
- IN Scene *scene,
- OUT float *t,
- OUT uint32_t *mesh_ind,
- OUT uint32_t *face_ind,
- OUT float *theta
-)
-{
- /* TODO BVH */
- /* for each triangle */
- Vec3 v1, v2, v3, e1, e2, re2_cross, s, se1_cross;
- float d, u, v;
- float dist = 1e9;
- float dist_tmp;
- for (uint32_t i = 0; i != scene->num_meshes; ++i) {
- Mesh *mesh = &scene->meshes[i];
- for (uint32_t j = 0; j != mesh->num_triangles; ++j) {
- v1 = mesh->vs[mesh->is[3 * j]];
- v2 = mesh->vs[mesh->is[3 * j + 1]];
- v3 = mesh->vs[mesh->is[3 * j + 2]];
- e1 = vec3_sub(&v2, &v1);
- e2 = vec3_sub(&v3, &v1);
- re2_cross = vec3_cross(&ray->d, &e2);
- d = vec3_dot(&e1, &re2_cross);
- if (d > -__FLT_EPSILON__ && d < __FLT_EPSILON__) continue;
- s = vec3_sub(&ray->o, &v1);
- u = vec3_dot(&s, &re2_cross) / d;
- if ((u < 0. && fabs(u) > __FLT_EPSILON__)
- || (u > 1. && fabs(u - 1.) > __FLT_EPSILON__))
- continue;
- se1_cross = vec3_cross(&s, &e1);
- v = vec3_dot(&ray->d, &se1_cross) / d;
- if ((v < 0. && fabs(v) > __FLT_EPSILON__)
- || (u + v > 1. && fabs(u + v - 1.) > __FLT_EPSILON__))
- continue;
- dist_tmp = vec3_dot(&e2, &se1_cross) / d;
- if (dist_tmp > __FLT_EPSILON__ && dist_tmp < dist) {
- dist = dist_tmp;
- *t = dist;
- *mesh_ind = i;
- *face_ind = j;
- *theta = acos(vec3_dot(&mesh->ns[j], &ray->d));
- if (*theta > PI / 2.)
- *theta = PI - *theta;
- }
- }
- }
-}
-
-/** Compute reflection R_{eTE} and R_{eTM} coefficients.
- *
- * Implements eqs. (31a)-(31b) from ITU-R P.2040-3.
- *
- * \param material_index the index of the material
- * \param theta1 the angle of incidence
- * \param r_te_re output real part of R_{eTE}
- * \param r_te_im output imaginary part of R_{eTE}
- * \param r_tm_re output real part of R_{eTM}
- * \param r_tm_im output imaginary part of R_{eTM}
- */
-void refl_coefs(
- IN uint32_t material_index,
- IN float theta1,
- OUT float *r_te_re,
- OUT float *r_te_im,
- OUT float *r_tm_re,
- OUT float *r_tm_im
-)
-{
- Material *rm = &g_materials[material_index];
- float sin_theta1 = sinf(theta1);
- if (rm->eta_abs_inv_sqrt * sin_theta1 > 1.f - __FLT_EPSILON__) {
- *r_te_re = *r_tm_re = 1.f;
- *r_te_im = *r_tm_im = 0.f;
- return;
- }
-
- /* eq. 33 */
- float sin_theta1_pow2 = sin_theta1 * sin_theta1;
- float cos_theta2_re = sqrtf(1.f + rm->eta_inv_re / rm->eta_abs_pow2 * sin_theta1_pow2);
- float cos_theta2_im = sqrtf(1.f - rm->eta_inv_im / rm->eta_abs_pow2 * sin_theta1_pow2);
-
- /* R_{eTE}, eq. 31a */
- float sqrt_eta_cos_theta2_re = rm->eta_sqrt_re * cos_theta2_re - rm->eta_sqrt_im * cos_theta2_im;
- float sqrt_eta_cos_theta2_im = rm->eta_sqrt_re * cos_theta2_im + rm->eta_sqrt_im * cos_theta2_re;
- float cos_theta1 = cosf(theta1);
- cdivf(cos_theta1 - sqrt_eta_cos_theta2_re, -sqrt_eta_cos_theta2_im,
- cos_theta1 + sqrt_eta_cos_theta2_re, sqrt_eta_cos_theta2_im,
- r_te_re, r_te_im);
-
- /* R_{eTM}, eq. 31b */
- float sqrt_eta_cos_theta1_re = rm->eta_sqrt_re * cos_theta1;
- float sqrt_eta_cos_theta1_im = rm->eta_sqrt_im * cos_theta1;
- cdivf(sqrt_eta_cos_theta1_re - cos_theta2_re,
- sqrt_eta_cos_theta1_im - cos_theta2_im,
- sqrt_eta_cos_theta1_re + cos_theta2_re,
- sqrt_eta_cos_theta1_im + cos_theta2_im,
- r_tm_re, r_tm_im);
-
- /* Apply reflection reduction factor */
- *r_te_re *= rm->r;
- *r_te_im *= rm->r;
- *r_tm_re *= rm->r;
- *r_tm_im *= rm->r;
-}
-
-/** Calculate scattering coefficients.
- *
- * Implements a directive scattering model inspired by Blaunstein et al. (DOI 10.1109/TAP.2006.888422).
- * Models scattering from rough surfaces (e.g., vegetation) with polarization dependence.
- *
- * \param theta_s scattering angle (radians, between scattered direction and surface normal)
- * \param theta_i angle of incidence (radians, between incident ray and surface normal)
- * \param material_index index into g_hrt_materials array
- * \param a_te_re output real part of TE scattering coefficient
- * \param a_te_im output imaginary part of TE scattering coefficient
- * \param a_tm_re output real part of TM scattering coefficient
- * \param a_tm_im output imaginary part of TM scattering coefficient
- */
-void scat_coefs(
- IN float theta_s,
- IN float theta_i,
- IN uint32_t material_index,
- OUT float *a_te_re,
- OUT float *a_te_im,
- OUT float *a_tm_re,
- OUT float *a_tm_im
-)
-{
- HRT_Material *mat = &g_hrt_materials[material_index];
-
- /* Precompute trigonometric values */
- float cos_theta_s = cosf(theta_s);
- float cos_theta_i = cosf(theta_i);
- float sin_theta_i = sinf(theta_i);
-
- /* Scattering directivity pattern: exponential decay with angle difference */
- /* f = s * exp(-alpha * |theta_s - theta_i|) models directive scattering */
- float theta_diff = fabsf(theta_s - theta_i);
- float f = mat->s * expf(-mat->s1_alpha * theta_diff);
-
- /* Roughness factor: reduces specular component as alpha increases */
- float roughness = 1.0f / (1.0f + mat->s1_alpha); // Simplified, 0 to 1
- float specular = roughness * cos_theta_s; // Specular-like term
- float diffuse = (1.0f - roughness) * cos_theta_s; // Diffuse term
-
- /* TE and TM coefficients (real parts) */
- /* TE: perpendicular to plane of incidence, less affected by angle */
- float te_real = f * (specular + diffuse);
- /* TM: parallel to plane of incidence, stronger angular dependence */
- float tm_real = f * (specular * cos_theta_i + diffuse);
-
- /* Phase shift (imaginary parts) due to surface roughness or path difference */
- /* Simplified: assume small phase shift proportional to roughness and angle */
- float phase_shift = mat->s1_alpha * sin_theta_i * 0.1f; // Arbitrary scaling
- float te_imag = te_real * sinf(phase_shift);
- float tm_imag = tm_real * sinf(phase_shift);
-
- /* Normalize to ensure energy conservation (optional, adjust as needed) */
- float norm = sqrtf(te_real * te_real + te_imag * te_imag +
- tm_real * tm_real + tm_imag * tm_imag);
- if (norm > 1e-6f) { // Avoid division by zero
- te_real /= norm;
- te_imag /= norm;
- tm_real /= norm;
- tm_imag /= norm;
- }
-
- /* Output coefficients */
- *a_te_re = te_real;
- *a_te_im = te_imag;
- *a_tm_re = tm_real;
- *a_tm_im = tm_imag;
-
- /* TODO: Incorporate s2 and s3 for additional roughness or frequency effects */
-}
-
-/* ==== MAIN FUNCTION ==== */
-
-void compute_paths(
- IN const char *scene_filepath, /* path to the scene file */
- IN const float *rx_pos, /* shape (num_rx, 3) */
- IN const float *tx_pos, /* shape (num_tx, 3) */
- IN const float *rx_vel, /* shape (num_rx, 3) */
- IN const float *tx_vel, /* shape (num_tx, 3) */
- IN float carrier_frequency, /* > 0.0 (IN GHz!) */
- IN size_t num_rx, /* number of receivers */
- IN size_t num_tx, /* number of transmitters */
- IN size_t num_paths, /* number of paths */
- IN size_t num_bounces, /* number of bounces */
- OUT PathsInfo *los, /* output LoS information. los->num_paths = 1 */
- OUT PathsInfo *scatter /* output scatter information. scatter->num_paths = num_bounces*num_paths */
-)
-{
- /* Load the scene */
- Scene scene = load_scene(scene_filepath, carrier_frequency);
-
- /* Precompute globals */
- /* g_free_space_loss_multiplier */
- g_free_space_loss_multiplier = 4.f * PI * carrier_frequency * 1e9 / SPEED_OF_LIGHT;
-
- /* Calculate a fibonacci sphere */
- Vec3 *ray_directions = (Vec3*)malloc(num_paths * sizeof(Vec3));
- float k, phi, theta;
- for (size_t path = 0; path < num_paths; ++path) {
- k = (float)path + .5f;
- phi = acos(1.f - 2.f * k / num_paths);
- theta = PI * (1.f + sqrtf(5.f)) * k;
- ray_directions[path] = (Vec3){
- cos(theta) * sin(phi),
- sin(theta) * sin(phi),
- cos(phi)
- };
- }
-
- /* Record the tx directions */
- /* This is only valid for the fibbonaci sphere */
- scatter->directions_tx = (float*)memcpy(
- scatter->directions_tx,
- ray_directions,
- num_paths * sizeof(Vec3)
- );
- if (scatter->directions_tx == NULL) {
- fprintf(stderr, "Failed to copy memory to scatter->directions_tx\n");
- exit(70);
- }
-
- /* Cast the positions and velocities to Vec3 */
- Vec3 *rx_pos_v = (Vec3*)rx_pos;
- Vec3 *tx_pos_v = (Vec3*)tx_pos;
- Vec3 *rx_vel_v = (Vec3*)rx_vel;
- Vec3 *tx_vel_v = (Vec3*)tx_vel;
-
- /* Create num_path rays for each tx */
- /* Shape (num_tx, num_paths) */
- Ray *rays = (Ray*)malloc(num_tx * num_paths * sizeof(Ray));
- for (size_t tx = 0, off = 0; tx < num_tx; ++tx) {
- for (size_t path = 0; path < num_paths; ++path, ++off) {
- rays[off].o = tx_pos_v[tx];
- rays[off].d = ray_directions[path];
- }
- }
- free(ray_directions);
-
- /* Initialize variables needed for the RT */
-
- /* Bouncing (reflecting) rays gains and delays */
- /* shape (num_tx, num_paths) */
- float *a_te_re_refl = (float*)malloc(num_tx * num_paths * sizeof(float));
- float *a_te_im_refl = (float*)calloc(num_tx * num_paths, sizeof(float));
- float *a_tm_re_refl = (float*)malloc(num_tx * num_paths * sizeof(float));
- float *a_tm_im_refl = (float*)calloc(num_tx * num_paths, sizeof(float));
- for (size_t i = 0; i < num_tx * num_paths; ++i)
- a_te_re_refl[i] = a_tm_re_refl[i] = 1.f;
- float *tau_t = (float*)calloc(num_tx * num_paths, sizeof(float));
-
- /* Active rays bitmask. 1 if the ray is still traced, 0 if it has left the scene */
- uint8_t *active = (uint8_t*)malloc((num_tx * num_paths / 8 + 1) * sizeof(uint8_t));
- for (size_t i = 0; i < num_tx * num_paths / 8 + 1; ++i)
- active[i] = 0xff;
- size_t num_active = num_tx * num_paths;
-
- /* Temporary variables */
- float t, r_te_re, r_te_im, r_tm_re, r_tm_im;
- float a_te_re_new, a_te_im_new, a_tm_re_new, a_tm_im_new;
- float theta_s; /* scattering angle */
- uint32_t mesh_ind, face_ind, material_index;
- Ray r;
- Vec3 *h = (Vec3*)malloc(sizeof(Vec3));
- Vec3 n;
-
- /* Friis free space loss multiplier.
- Multiply this by a distance and take the second power to get a free space loss. */
- float free_space_loss_multiplier = 4.f * PI * carrier_frequency * 1e9 / SPEED_OF_LIGHT;
- float free_space_loss;
-
- /* Doppler frequency shift carrier frequency multiplier */
- float doppler_shift_multiplier = carrier_frequency * 1e9 / SPEED_OF_LIGHT;
-
- /* Initialize the doppler frequency shift using the tx velocities */
- /* The scatter->freq_shift has the shape of (num_rx, num_tx, num_bounces*num_paths) */
- /* The first num_tx*num_paths elements are the same for each rx and bounce */
- /* so we can initialize them only once and memcpy them to the rest */
- for (size_t tx = 0; tx != num_tx; ++tx)
- for (size_t path = 0; path != num_paths; ++path) {
- size_t off_rays = tx * num_paths + path;
- size_t off_scat = tx * num_paths * num_bounces + path;
- scatter->freq_shift[off_scat] = vec3_dot(&tx_vel_v[tx], &rays[off_rays].d);
- scatter->freq_shift[off_scat] *= doppler_shift_multiplier;
- }
- for (size_t bounce = 1; bounce < num_bounces; ++bounce)
- if (!memcpy(scatter->freq_shift + num_tx * num_paths * bounce,
- scatter->freq_shift, num_tx * num_paths * sizeof(float)))
- exit(70);
- for (size_t rx = 1; rx < num_rx; ++rx)
- if (!memcpy(scatter->freq_shift + num_tx * num_paths * num_bounces * rx,
- scatter->freq_shift, num_tx * num_paths * num_bounces * sizeof(float)))
- exit(70);
-
- /***************************************************************************/
- /* LoS paths */
- /***************************************************************************/
-
- /* Consider imaginary parts of the LoS gains to be 0 in any way */
- for (size_t off = 0; off != num_rx * num_tx; ++off)
- los->a_te_im[off] = los->a_tm_im[off] = 0.f;
-
- /* Calculate LoS paths directly from tx to rx */
- for (size_t rx = 0, off = 0; rx != num_rx; ++rx)
- for (size_t tx = 0; tx != num_tx; ++tx, ++off) {
- t = -1.f;
- mesh_ind = face_ind = -1;
-
- /* Create a ray from the tx to the rx */
- r.o = tx_pos_v[tx];
- r.d = vec3_sub(&rx_pos_v[rx], &r.o);
-
- /* In case the tx and the rx are at the same position */
- if (vec3_dot(&r.d, &r.d) < __FLT_EPSILON__) {
- /* Assume the direction to be (1, 0, 0) */
- los->directions_rx[off * 3] = 1.f;
- los->directions_tx[off * 3] = -1.f;
- los->directions_rx[off * 3 + 1] = los->directions_tx[off * 3 + 1] = 0.f;
- los->directions_rx[off * 3 + 2] = los->directions_tx[off * 3 + 2] = 0.f;
- /* Set gains to 1+0j */
- los->a_te_re[off] = los->a_tm_re[off] = 1.f;
- /* Set delay to 0 */
- los->tau[off] = 0.f;
- /* Set doppler shift to 0. */
- los->freq_shift[off] = 0.f;
- continue;
- }
-
- /* Is there any abstacle between the tx and the rx? */
- moeller_trumbore(&r, &scene, &t, &mesh_ind, &face_ind, &theta);
- if (mesh_ind != (uint32_t)-1 && t <= 1.f) {
- /* An obstacle between the tx and the rx has been hit */
- los->a_te_re[off] = los->a_tm_re[off] = los->tau[off] = 0.f;
- continue;
- }
-
- /* No obstacle has been hit, the path is LoS */
- /* Calculate the distance */
- t = sqrtf(vec3_dot(&r.d, &r.d));
- /* Calculate the direction */
- Vec3 d = {r.d.x / t, r.d.y / t, r.d.z / t};
- los->directions_tx[off * 3] = d.x;
- los->directions_tx[off * 3 + 1] = d.y;
- los->directions_tx[off * 3 + 2] = d.z;
- los->directions_rx[off * 3] = -d.x;
- los->directions_rx[off * 3 + 1] = -d.y;
- los->directions_rx[off * 3 + 2] = -d.z;
- /* Calculate the free space loss */
- free_space_loss = free_space_loss_multiplier * t;
- if (free_space_loss > 1.f) {
- los->a_te_re[off] = 1.f / free_space_loss;
- los->a_tm_re[off] = 1.f / free_space_loss;
- } else
- los->a_te_re[off] =los-> a_tm_re[off] = 1.f;
- /* Calculate the delay */
- los->tau[off] = t / SPEED_OF_LIGHT;
- /* Calculate the doppler shift */
- los->freq_shift[off] = vec3_dot(tx_vel_v, &d) - vec3_dot(rx_vel_v, &d);
- los->freq_shift[off] *= carrier_frequency * 1e9 / SPEED_OF_LIGHT;
- }
-
- /***************************************************************************/
- /* Scatter paths */
- /***************************************************************************/
-
- /* Bounce the rays. On each bounce:
- * - Add path length / SPEED_OF_LIGHT to tau
- * - Update gains using eqs. (31a)-(31b) from ITU-R P.2040-3
- * - Spawn scattering rays towards the rx
- * - TODO Spawn refraction rays with gains as per eqs. (31c)-(31d) from ITU-R P.2040-3
- */
- FILE* ray_file = fopen("rays.bin", "wb");
- if (!ray_file) {
- fprintf(stderr, "Could not open file rays.bin\n");
- exit(8);
- }
- fwrite(rays, sizeof(Ray), num_tx * num_paths, ray_file);
-
- FILE* active_file = fopen("active.bin", "wb");
- if (!active_file) {
- fprintf(stderr, "Could not open file active.bin\n");
- exit(8);
- }
- fwrite(active, sizeof(uint8_t), num_tx * num_paths / 8 + 1, active_file);
-
- for (size_t bounce = 0; bounce < num_bounces; ++bounce) {
- /* active[active_byte_index] & (1 << active_bit_pos) */
- size_t active_byte_index = 0;
- uint8_t active_bit = 1;
- size_t off_tx_path = 0; /* tx * num_paths + path */
- for (size_t tx = 0; tx < num_tx; ++tx)
- for (size_t path = 0; path < num_paths; ++path, ++off_tx_path, active_bit <<= 1) {
-
- /* Check if the ray is active */
- if (!active_bit) {
- active_bit = 1;
- ++active_byte_index;
- }
- if (!(active[active_byte_index] & active_bit))
- continue;
-
- /***********************************************************************/
- /* Reflect the rays */
- /***********************************************************************/
-
- /* Init */
- t = -1.f;
- mesh_ind = face_ind = -1;
- /* Find the hit point and trinagle and the angle of incidence */
- moeller_trumbore(&rays[off_tx_path], &scene, &t, &mesh_ind, &face_ind, &theta);
- if (mesh_ind == (uint32_t)-1) { /* Ray hit nothing */
- active[active_byte_index] &= ~active_bit;
- --num_active;
- continue;
- }
- /* Calculate the reflection coefficients R_{eTE} and R_{eTM} */
- material_index = scene.meshes[mesh_ind].material_index;
- refl_coefs(material_index, theta,
- &r_te_re, &r_te_im,
- &r_tm_re, &r_tm_im);
- /* Calculate the free space loss */
- free_space_loss = free_space_loss_multiplier * t;
- free_space_loss *= free_space_loss;
- if (free_space_loss > 1.f) {
- r_te_re /= free_space_loss;
- r_te_im /= free_space_loss;
- r_tm_re /= free_space_loss;
- r_tm_im /= free_space_loss;
- }
- /* Update the gains as a' = a * refl_coefs */
- a_te_re_new = a_te_re_refl[off_tx_path] * r_te_re - a_te_im_refl[off_tx_path] * r_te_im;
- a_te_im_new = a_te_re_refl[off_tx_path] * r_te_im + a_te_im_refl[off_tx_path] * r_te_re;
- a_tm_re_new = a_tm_re_refl[off_tx_path] * r_tm_re - a_tm_im_refl[off_tx_path] * r_tm_im;
- a_tm_im_new = a_tm_re_refl[off_tx_path] * r_tm_im + a_tm_im_refl[off_tx_path] * r_tm_re;
- a_te_re_refl[off_tx_path] = a_te_re_new;
- a_te_im_refl[off_tx_path] = a_te_im_new;
- a_tm_re_refl[off_tx_path] = a_tm_re_new;
- a_tm_im_refl[off_tx_path] = a_tm_im_new;
- /* Update the delay */
- tau_t[off_tx_path] += t / SPEED_OF_LIGHT;
-
- /* Move and reflect the ray */
- r = rays[off_tx_path];
- /* Advance the ray to the hit point */
- *h = vec3_scale(&r.d, t);
- r.o = vec3_add(h, &r.o);
- /* Reflect the ray's direction as d' = d - 2*(d \cdot n)*n */
- n = scene.meshes[mesh_ind].ns[face_ind];
- t = vec3_dot(&r.d, &n);
- *h = vec3_scale(&n, 2.f * t);
- r.d = vec3_sub(&r.d, h);
- /* Advance the ray origin a bit to avoid intersection with the same triangle */
- *h = vec3_scale(&r.d, 1e-4f);
- r.o = vec3_add(&r.o, h);
-
- /* Calculate the doppler shift caused by the reflection */
- Vec3* mesh_vel = (Vec3*)&scene.meshes[mesh_ind].velocity;
- Vec3 d_rd = vec3_sub(&r.d, &rays[off_tx_path].d);
- scatter->freq_shift[off_tx_path] += vec3_dot(&d_rd, mesh_vel) * doppler_shift_multiplier;
-
- /* Save the reflected ray */
- rays[off_tx_path] = r;
-
- /***********************************************************************/
- /* Scatter the rays */
- /***********************************************************************/
-
- /* Scatter a path from the hit point towards the rx */
- Ray r_scat = { .o = r.o };
- for (size_t rx = 0; rx < num_rx; ++rx) {
- /* [rx, tx, bounce, path] */
- size_t off_scat = ((rx * num_tx + tx) * num_bounces + bounce) * num_paths + path;
- /* Create a ray from the hit point to the rx */
- r_scat.d = vec3_sub(&rx_pos_v[rx], &r_scat.o);
- float d2rx = sqrtf(vec3_dot(&r_scat.d, &r_scat.d));
- r_scat.d = vec3_normalize(&r_scat.d);
- /* Check if the ray has a LoS of the rx */
- t = -1.f;
- mesh_ind = face_ind = -1;
- moeller_trumbore(&r_scat, &scene, &t, &mesh_ind, &face_ind, &theta);
- if (mesh_ind != (uint32_t)-1 && t <= 1.f) {
- /* An obstacle between the hit point and the rx has been hit */
- scatter->a_te_re[off_scat] = scatter->a_te_im[off_scat]
- = scatter->a_tm_re[off_scat]
- = scatter->a_tm_im[off_scat]
- = scatter->tau[off_scat]
- = 0.f;
- continue;
- }
- /* No obstacle has been hit */
- /* Calculate the scattering angle */
- theta_s = acosf(vec3_dot(&r_scat.d, &n));
- /* Calculate the scattering coefficients */
- float a_te_re_new, a_te_im_new, a_tm_re_new, a_tm_im_new;
- scat_coefs(theta_s, theta, material_index,
- &a_te_re_new, &a_te_im_new, &a_tm_re_new, &a_tm_im_new);
- scatter->a_te_re[off_scat] = a_te_re_refl[off_tx_path] * a_te_re_new
- - a_te_im_refl[off_tx_path] * a_te_im_new;
- scatter->a_te_im[off_scat] = a_te_re_refl[off_tx_path] * a_te_im_new
- + a_te_im_refl[off_tx_path] * a_te_re_new;
- scatter->a_tm_re[off_scat] = a_tm_re_refl[off_tx_path] * a_tm_re_new
- - a_tm_im_refl[off_tx_path] * a_tm_im_new;
- scatter->a_tm_im[off_scat] = a_tm_re_refl[off_tx_path] * a_tm_im_new
- + a_tm_im_refl[off_tx_path] * a_tm_re_new;
- /* Calculate the direction */
- scatter->directions_rx[off_scat * 3] = -r_scat.d.x;
- scatter->directions_rx[off_scat * 3 + 1] = -r_scat.d.y;
- scatter->directions_rx[off_scat * 3 + 2] = -r_scat.d.z;
- /* Calculate the delay */
- scatter->tau[off_scat] = tau_t[off_tx_path] + d2rx / SPEED_OF_LIGHT;
- /* Calculate the free space loss */
- free_space_loss = free_space_loss_multiplier * d2rx;
- free_space_loss *= free_space_loss;
- if (free_space_loss > 1.f) {
- scatter->a_te_re[off_scat] /= free_space_loss;
- scatter->a_te_im[off_scat] /= free_space_loss;
- scatter->a_tm_re[off_scat] /= free_space_loss;
- scatter->a_tm_im[off_scat] /= free_space_loss;
- }
- /* Calculate the doppler shift */
- Vec3 d_rd = vec3_sub(&r_scat.d, &rays[off_tx_path].d);
- float freq_shift = vec3_dot(&d_rd, mesh_vel) * doppler_shift_multiplier;
- scatter->freq_shift[off_scat] -= freq_shift;
- }
-
- /***********************************************************************/
- /* Refract the rays */
- /***********************************************************************/
- /* TODO */
- }
- fwrite(rays, sizeof(Ray), num_tx * num_paths, ray_file);
- fwrite(active, sizeof(uint8_t), num_tx * num_paths / 8 + 1, active_file);
- }
- fclose(ray_file);
- fclose(active_file);
-
- /* Free */
- free(a_te_re_refl);
- free(a_te_im_refl);
- free(a_tm_re_refl);
- free(a_tm_im_refl);
- free(tau_t);
- free(h);
- free(active);
- free(rays);
- /* Free the scene */
- /* TODO */
-}
diff --git a/inc/common.h b/inc/common.h
@@ -1,16 +1,25 @@
+#ifndef COMMON_H
+#define COMMON_H
+
#define IN
#define OUT
+#include <stdlib.h> /* for exit */
+
#define FREE_POINTERS(...) \
- do { \
- void *ptrs[] = {__VA_ARGS__};\
- for (size_t i = 0; i < sizeof(ptrs) / sizeof(ptrs[0]); i++) \
- free(ptrs[i]); \
- } while (0)
+ do { \
+ void *ptrs[] = {__VA_ARGS__};\
+ size_t cnt = sizeof(ptrs) / sizeof(ptrs[0]); \
+ for (size_t i = 0; i < cnt; i++) \
+ free(ptrs[i]); \
+ } while (0)
+/* Variadic arguments after rc are the pointers to free */
#define PERROR_CLEANUP_EXIT(msg, rc, ...) \
- do { \
- perror(msg); \
- FREE_POINTERS(__VA_ARGS__); \
- exit(rc); \
- } while (0)
+ do { \
+ perror(msg); \
+ FREE_POINTERS(__VA_ARGS__); \
+ exit(rc); \
+ } while (0)
+
+#endif
diff --git a/inc/compute_paths.h b/inc/compute_paths.h
@@ -1,6 +1,8 @@
#ifndef COMPUTE_PATHS_H
#define COMPUTE_PATHS_H
+#include "scene.h" /* for Scene */
+
#include <stddef.h> /* for size_t */
#include <stdint.h> /* for uint32_t */
@@ -31,7 +33,7 @@ typedef struct {
*
* The output parameters are allocated by the caller (including the fields).
*
- * \param scene_filepath path to a scene .hrt file
+ * \param scene pointer to a loaded scene in HRT format (see scene.h)
* \param rx_pos receiver positions, shape (num_rx, 3)
* \param tx_pos transmitter positions, shape (num_tx, 3)
* \param rx_vel receiver velocities, shape (num_rx, 3)
@@ -46,18 +48,18 @@ typedef struct {
* \param scatter output scatter results. scatter->num_paths = num_bounces * num_paths
*/
void compute_paths(
- IN const char *scene_filepath, /* path to the scene file */
- IN const float *rx_pos, /* shape (num_rx, 3) */
- IN const float *tx_pos, /* shape (num_tx, 3) */
- IN const float *rx_vel, /* shape (num_rx, 3) */
- IN const float *tx_vel, /* shape (num_tx, 3) */
- IN float carrier_frequency, /* > 0.0 (IN GHz!) */
- IN size_t num_rx, /* number of receivers */
- IN size_t num_tx, /* number of transmitters */
- IN size_t num_paths, /* number of paths */
- IN size_t num_bounces, /* number of bounces */
- OUT PathsInfo *los, /* output LoS information */
- OUT PathsInfo *scatter /* output scatter information */
+ IN Scene *scene, /* Pointer to a loaded scene */
+ IN const float *rx_pos, /* shape (num_rx, 3) */
+ IN const float *tx_pos, /* shape (num_tx, 3) */
+ IN const float *rx_vel, /* shape (num_rx, 3) */
+ IN const float *tx_vel, /* shape (num_tx, 3) */
+ IN float carrier_frequency, /* > 0.0 (IN GHz!) */
+ IN size_t num_rx, /* number of receivers */
+ IN size_t num_tx, /* number of transmitters */
+ IN size_t num_paths, /* number of paths */
+ IN size_t num_bounces, /* number of bounces */
+ OUT PathsInfo *los, /* output LoS information */
+ OUT PathsInfo *scatter /* output scatter information */
);
#endif /* COMPUTE_PATHS_H */
diff --git a/inc/materials.h b/inc/materials.h
@@ -1,100 +1,14 @@
#ifndef MATERIALS_H
#define MATERIALS_H
-#include "scene.h" /* for HRT_Material */
+#include "scene.h" /* for Material */
#include <stdint.h> /* for uint32_t */
#include <string.h> /* for strncmp */
/* number of defined materials */
#define NUM_G_MATERIALS 17
-HRT_Material g_hrt_materials[NUM_G_MATERIALS] = {
- /* 0 */
- {3, "air",
- 1.f, 0.f, 0.f, 0.001f,
- 0.1f, 0.5f, 0.3f, 0.2f,
- 2, 2},
- /* 1 */
- {8, "concrete",
- 5.24f, 0.f, 0.0462f, 0.7822f,
- 0.5f, 0.33f, 0.34f, 0.33f,
- 4, 4},
- /* 2 */
- {5, "brick",
- 3.91f, 0.f, 0.0238f, 0.16f,
- 0.4f, 0.4f, 0.3f, 0.3f,
- 3, 3},
- /* 3 */
- {12, "plasterboard",
- 2.73f, 0.f, 0.0085f, 0.9395f,
- 0.3f, 0.4f, 0.4f, 0.2f,
- 3, 3},
- /* 4 */
- {4, "wood",
- 1.99f, 0.f, 0.0047f, 1.0718f,
- 0.2f, 0.5f, 0.3f, 0.2f,
- 2, 2},
- /* 5 */
- {5, "glass",
- 6.31f, 0.f, 0.0036f, 1.3394f,
- 0.3f, 0.4f, 0.4f, 0.2f,
- 3, 3},
- /* 6 */
- {5, "glass",
- 5.79f, 0.f, 0.0004f, 1.658f,
- 0.3f, 0.4f, 0.4f, 0.2f,
- 3, 3},
- /* 7 */
- {13, "ceiling board",
- 1.48f, 0.f, 0.0011f, 1.0750f,
- 0.2f, 0.5f, 0.3f, 0.2f,
- 2, 2},
- /* 8 */
- {13, "ceiling board",
- 1.52f, 0.f, 0.0029f, 1.029f,
- 0.2f, 0.5f, 0.3f, 0.2f,
- 2, 2},
- /* 9 */
- {9, "chipboard",
- 2.58f, 0.f, 0.0217f, 0.7800f,
- 0.4f, 0.4f, 0.3f, 0.3f,
- 3, 3},
- /* 10 */
- {7, "plywood",
- 2.71f, 0.f, 0.33f, 0.f,
- 0.3f, 0.5f, 0.3f, 0.2f,
- 3, 3},
- /* 11 */
- {6, "marble",
- 7.074f, 0.f, 0.0055f,
- 0.9262f, 0.3f, 0.4f, 0.4f, 0.2f,
- 3, 3},
- /* 12 */
- {10, "floorboard",
- 3.66f, 0.f, 0.0044f, 1.3515f,
- 0.3f, 0.4f, 0.4f, 0.2f,
- 3, 3},
- /* 13 */
- {5, "metal",
- 1.f, 0.f, 10000000.f, 0.f,
- 0.f, 0.f, 1.f, 0.f,
- 1, 1},
- /* 14 */
- {15, "very dry ground",
- 3.f, 0.f, 0.00015f, 2.52f,
- 0.4f, 0.3f, 0.4f, 0.3f,
- 4, 4},
- /* 15 */
- {17, "medium dry ground",
- 15.f, -0.1f, 0.035f, 1.63f,
- 0.5f, 0.33f, 0.34f, 0.33f,
- 4, 4},
- /* 16 */
- {10, "wet ground",
- 30.f, -0.4f, 0.15f, 1.30f,
- 0.5f, 0.33f, 0.34f, 0.33f,
- 4, 4}
-};
+extern Material g_materials[NUM_G_MATERIALS];
typedef enum {
MATERIAL_AIR = 0,
MATERIAL_CONCRETE = 1,
@@ -115,38 +29,7 @@ typedef enum {
MATERIAL_WET_GROUND = 16
} MaterialIndex;
-/* Map material name strings to material indices */
-typedef struct {
- /* Material name. Null-terminated string. Lower case. */
- const char *name;
- MaterialIndex index;
-} MaterialMap;
-
-MaterialMap material_map[] = {
- {"air", MATERIAL_AIR},
- {"concrete", MATERIAL_CONCRETE},
- {"brick", MATERIAL_BRICK},
- {"plasterboard", MATERIAL_PLASTERBOARD},
- {"wood", MATERIAL_WOOD},
- {"glass1", MATERIAL_GLASS1},
- {"glass2", MATERIAL_GLASS2},
- {"ceiling_board1", MATERIAL_CEILING_BOARD1},
- {"ceiling_board2", MATERIAL_CEILING_BOARD2},
- {"chipboard", MATERIAL_CHIPBOARD},
- {"plywood", MATERIAL_PLYWOOD},
- {"marble", MATERIAL_MARBLE},
- {"floorboard", MATERIAL_FLOORBOARD},
- {"metal", MATERIAL_METAL},
- {"very_dry_ground", MATERIAL_VERY_DRY_GROUND},
- {"medium_dry_ground", MATERIAL_MEDIUM_DRY_GROUND},
- {"wet_ground", MATERIAL_WET_GROUND}
-};
-
-MaterialIndex get_material_index(const char *name) {
- for (uint32_t i = 0; i < NUM_G_MATERIALS; ++i)
- if (strcmp(material_map[i].name, name) == 0)
- return material_map[i].index;
- return MATERIAL_AIR; /* Invalid material. Default to air */
-}
+/** Get the material index from the material name */
+extern MaterialIndex get_material_index(const char *name);
#endif /* MATERIALS_H */
diff --git a/inc/scene.h b/inc/scene.h
@@ -1,27 +1,34 @@
#ifndef SCENE_H
#define SCENE_H
+#include "common.h" /* for IN */
+#include "vec3.h" /* for Vec3 */
+
#include <stdint.h> /* for uint32_t */
typedef struct {
/* Number of vertices */
uint32_t num_vertices;
- /* Vertices coordinates. Size [num_vertices * 3] */
- float *vs;
+ /* Vertices coordinates. Size [num_vertices] */
+ Vec3 *vs;
/* Number of triangles */
uint32_t num_triangles;
/* Triangles indices. Size [num_triangles * 3] */
uint32_t *is;
/* Index of the material in g_hrt_materials array defined in scene.h */
uint32_t material_index;
- /* Global cartesian velocity vector. Size [3] */
- float velocity[3];
-} HRT_Mesh;
+ /* Global cartesian velocity vector */
+ Vec3 velocity;
+ /* Normals of the triangles. Size [num_triangles] */
+ /* NOTE: This field is not stored in the HRT file.
+ The loader does not load, allocate or calculate it. */
+ Vec3 *ns;
+} Mesh;
typedef struct {
uint32_t num_meshes;
- HRT_Mesh *meshes;
-} HRT_Scene;
+ Mesh *meshes;
+} Scene;
typedef struct {
/* Number of characters in the name */
@@ -55,6 +62,25 @@ typedef struct {
* Must be > 0.
*/
uint8_t s3_alpha;
-} HRT_Material;
+} Material;
+
+/** Save a scene to a HRT file. (See README for details)
+ *
+ * NOTE: This function does not save the normals of the triangles (Mesh.ns field).
+ *
+ * \param scene pointer to the scene to save
+ * \param filepath path to the HRT file. Will be overwritten if exists.
+ */
+void scene_save(IN Scene* scene, IN const char* filepath);
+
+/** Load a mesh from a HRT file.
+ *
+ * NOTE: This function does not load, allocate or calculate
+ * the normals of the triangles (Mesh.ns field).
+ *
+ * \param filepath path to the HRT file
+ * \return the loaded scene
+ */
+Scene scene_load(IN const char *filepath);
#endif /* SCENE_H */
diff --git a/inc/vec3.h b/inc/vec3.h
@@ -0,0 +1,40 @@
+#ifndef VEC3_H
+#define VEC3_H
+
+#include <math.h> /* for sqrtf */
+
+typedef struct {
+ float x, y, z;
+} Vec3;
+
+static inline Vec3 vec3_sub(const Vec3 *a, const Vec3 *b)
+{
+ return (Vec3){a->x - b->x, a->y - b->y, a->z - b->z};
+}
+static inline Vec3 vec3_add(const Vec3 *a, const Vec3 *b)
+{
+ return (Vec3){a->x + b->x, a->y + b->y, a->z + b->z};
+}
+static inline Vec3 vec3_cross(const Vec3 *a, const Vec3 *b)
+{
+ return (Vec3){
+ a->y * b->z - a->z * b->y,
+ a->z * b->x - a->x * b->z,
+ a->x * b->y - a->y * b->x
+ };
+}
+static inline float vec3_dot(const Vec3 *a, const Vec3 *b)
+{
+ return a->x * b->x + a->y * b->y + a->z * b->z;
+}
+static inline Vec3 vec3_scale(const Vec3 *a, float s)
+{
+ return (Vec3){a->x * s, a->y * s, a->z * s};
+}
+static inline Vec3 vec3_normalize(const Vec3 *a)
+{
+ float norm = sqrtf(a->x * a->x + a->y * a->y + a->z * a->z);
+ return (Vec3){a->x / norm, a->y / norm, a->z / norm};
+}
+
+#endif /* VEC3_H */
+\ No newline at end of file
diff --git a/scene.hrt b/scene.hrt
Binary files differ.
diff --git a/scene_fromSionna.c b/scene_fromSionna.c
@@ -1,503 +0,0 @@
-/* vim: set tabstop=2:softtabstop=2:shiftwidth=2:noexpandtab */
-
-#include "inc/common.h" /* for IN, OUT, PERROR_CLEANUP_EXIT */
-#include "inc/scene.h" /* for HRT_Scene, HRT_Mesh, HRT_Material */
-#include "inc/materials.h" /* for g_hrt_materials, MATERIAL_CONCRETE */
-
-#include <stdio.h> /* for FILE, fopen, fclose, fseek, fgets, sscanf, printf */
-#include <stdlib.h> /* for malloc */
-#include <stdint.h> /* for uint8_t, uint32_t */
-#include <string.h> /* for strncmp, strrchr, strcmp */
-
-/**
- * Hardcoded scenes
- */
-
-/* Box */
-float mesh_box_vs[] = {
- 5.f, 5.f, 0.f,
- -5.f, 5.f, 0.f,
- -5.f, -5.f, 0.f,
- 5.f, -5.f, 0.f,
- 5.f, 5.f, 5.f,
- -5.f, 5.f, 5.f,
- -5.f, -5.f, 5.f,
- 5.f, -5.f, 5.f,
-};
-uint32_t mesh_box_is[] = {
- 0, 1, 2,
- 0, 2, 3,
- 0, 4, 5,
- 0, 5, 1,
- 1, 5, 6,
- 1, 6, 2,
- 2, 6, 7,
- 2, 7, 3,
- 3, 7, 4,
- 3, 4, 0,
- 4, 7, 6,
- 4, 6, 5,
-};
-HRT_Mesh mesh_box = {
- .num_vertices = 8,
- .vs = (float*)mesh_box_vs,
- .num_triangles = 12,
- .is = (uint32_t*)mesh_box_is,
- .material_index = MATERIAL_CONCRETE,
- .velocity = {0.f, 0.f, 0.f},
-};
-HRT_Scene scene_box = {
- .num_meshes = 1,
- .meshes = &mesh_box,
-};
-const char* scene_box_filename = "box.xml";
-
-/* Simple reflector */
-float mesh_simpleReflector_vs[] = {
- -0.5f, -0.5f, 0.f,
- 0.5f, -0.5f, 0.f,
- 0.5f, 0.5f, 0.f,
- -0.5f, 0.5f, 0.f,
-};
-uint32_t mesh_simpleReflector_is[] = {
- 0, 1, 2,
- 0, 2, 3,
-};
-HRT_Mesh mesh_simpleReflector = {
- .num_vertices = 4,
- .vs = (float*)mesh_simpleReflector_vs,
- .num_triangles = 2,
- .is = (uint32_t*)mesh_simpleReflector_is,
- .material_index = MATERIAL_CONCRETE,
- .velocity = {0.f, 0.f, 0.f},
-};
-HRT_Scene scene_simpleReflector = {
- .num_meshes = 1,
- .meshes = &mesh_simpleReflector,
-};
-const char* scene_simpleReflector_filename = "simple_reflector.xml";
-
-/**
- * Read a Sionna scene
- */
-
-/* Read a binary PLY file
- *
- * The PLY file have the following header:
- * ply
- * format binary_little_endian 1.0
- * element vertex <num_vertices>
- * property float x
- * property float y
- * property float z
- * property float s
- * property float t
- * element face <num_triangles>
- * property list uchar int vertex_index
- * end_header
- *
- * \param filepath Path to the PLY file
- * \return The mesh
- */
-HRT_Mesh readPly(const char* filepath) {
- FILE* f = fopen(filepath, "rb");
- if (f == NULL)
- PERROR_CLEANUP_EXIT("Error: cannot open file", 8);
-
- char buf[256];
- HRT_Mesh mesh = {
- .num_vertices = 0,
- .vs = NULL,
- .num_triangles = 0,
- .is = NULL,
- .material_index = 0,
- .velocity = {0.f, 0.f, 0.f},
- };
-
- /* read header */
- while(fgets(buf, 256, f)) {
- if (!strncmp(buf, "end_header", 10))
- break;
- if (!strncmp(buf, "element vertex ", 15)) {
- if (sscanf(buf, "element vertex %u", &mesh.num_vertices) != 1)
- PERROR_CLEANUP_EXIT("Error: cannot read number of vertices", 8);
- }
- else if (!strncmp(buf, "element face ", 13))
- if (sscanf(buf, "element face %u", &mesh.num_triangles) != 1)
- PERROR_CLEANUP_EXIT("Error: cannot read number of triangles", 8);
- }
-
- /* check if vertex and face elements are found */
- if (mesh.num_vertices == 0 || mesh.num_triangles == 0)
- PERROR_CLEANUP_EXIT("Error: PLY element vertex or element face not found", 8);
- /* check if the numbers are not too big */
- else if (mesh.num_vertices > 1000000 || mesh.num_triangles > 1000000)
- PERROR_CLEANUP_EXIT("Error: PLY element vertex or element face too big", 8);
-
- /* allocate memory for vertices and faces */
- mesh.vs = (float*)malloc(3 * mesh.num_vertices * sizeof(float));
- mesh.is = (uint32_t*)malloc(3 * mesh.num_triangles * sizeof(uint32_t));
-
- /* read vertices */
- for (uint32_t i = 0; i != mesh.num_vertices * 3; i += 3) {
- if (fread(&mesh.vs[i], sizeof(float), 3, f) != 3)
- PERROR_CLEANUP_EXIT("Error: cannot read vertex", 8, mesh.vs, mesh.is);
- /* skip s and t */
- if (fseek(f, 8, SEEK_CUR) != 0)
- PERROR_CLEANUP_EXIT("Error: cannot skip s and t", 8, mesh.vs, mesh.is);
- }
-
- /* read faces */
- for (uint32_t i = 0; i != mesh.num_triangles * 3; i += 3) {
- uint8_t n;
- if (fread(&n, sizeof(uint8_t), 1, f) != 1)
- PERROR_CLEANUP_EXIT("Error: cannot read number of vertices in face", 8, mesh.vs, mesh.is);
- if (n != 3)
- PERROR_CLEANUP_EXIT("Error: face is not a triangle", 8, mesh.vs, mesh.is);
- if (fread(&mesh.is[i], sizeof(uint32_t), 3, f) != 3)
- PERROR_CLEANUP_EXIT("Error: cannot read face", 8, mesh.vs, mesh.is);
- }
-
- fclose(f);
- return mesh;
-}
-
-/** Read a scene config CSV file.
- *
- * The CSV file must be a text file with the following format:
- * name,material_index,velocity_x,velocity_y,velocity_z
- *
- * The file can have any number of lines except the first one, which is the header.
- * If a mesh is not in the CSV file,
- * it will have the default (0) material index and velocity.
- *
- * \param filepath Path to the CSV file
- * \param num_meshes Output number of meshes in the CSV file
- * \param mesh_filenames Output array of mesh filenames
- * \param material_indicies Output array of material indicies
- * \param velocity Output array of velocities. Size [num_meshes * 3]
- */
-void readCsv(
- IN const char* filepath,
- OUT uint32_t* csv_num_meshes,
- OUT char** csv_mesh_filenames,
- OUT uint32_t* csv_material_indicies,
- OUT float* csv_velocities
-)
-{
- FILE* f = fopen(filepath, "r");
- if (!f)
- PERROR_CLEANUP_EXIT("Error: cannot open file", 8);
-
- /* Check header */
- char buf[256];
- if (!fgets(buf, 256, f))
- PERROR_CLEANUP_EXIT("Error: cannot read header", 8);
- if (strncmp(buf, "name,material_index,velocity_x,velocity_y,velocity_z\n", 48))
- PERROR_CLEANUP_EXIT("Error: invalid header", 8);
-
- /* Count number of lines */
- *csv_num_meshes = 0;
- while (fgets(buf, 256, f))
- ++(*csv_num_meshes);
- if (*csv_num_meshes == 0) {
- fclose(f);
- return;
- }
- /* Return to the beginning of the file */
- rewind(f);
- /* Skip header */
- if (!fgets(buf, 256, f))
- PERROR_CLEANUP_EXIT("Error: cannot read header", 8);
-
- /* Allocate memory for the arrays */
- csv_mesh_filenames = (char**)malloc(*csv_num_meshes * sizeof(char*));
- csv_material_indicies = (uint32_t*)malloc(*csv_num_meshes * sizeof(uint32_t));
- csv_velocities = (float*)malloc(*csv_num_meshes * 3 * sizeof(float));
-
- /* Read lines */
- for (uint32_t i = 0; i != *csv_num_meshes; ++i) {
- /* Read line */
- if (!fgets(buf, 256, f))
- PERROR_CLEANUP_EXIT("Error: cannot read line", 8);
- /* Parse line */
- char name[50];
- int rc = sscanf(
- buf,
- "%49[^,],%u,%f,%f,%f\n",
- name,
- &csv_material_indicies[i],
- &csv_velocities[i * 3],
- &csv_velocities[i * 3 + 1],
- &csv_velocities[i * 3 + 2]
- );
- if (rc != 4)
- PERROR_CLEANUP_EXIT("Error: cannot parse line", 8);
- /* Copy name */
- csv_mesh_filenames[i] = (char*)malloc(strlen(name) + 1);
- if (!csv_mesh_filenames[i])
- PERROR_CLEANUP_EXIT("Error: cannot allocate memory for mesh filename", 8);
- strcpy(csv_mesh_filenames[i], name);
- }
-}
-
-/** Read a Sionna scene XML file.
- *
- * Extracts each shape's name, file path and material.
- *
- * \param filepath Path to the XML file
- * \param xml_num_meshes Output number of meshes in the XML file
- * \param xml_mesh_filepaths Output array of mesh file paths (null-terminated strings)
- * \param xml_mesh_names Output array of mesh names (null-terminated strings)
- * \param xml_mesh_materials Output array of mesh materials (null-terminated strings)
- */
-void readXml(
- IN const char* filepath,
- OUT uint32_t* xml_num_meshes,
- OUT char*** xml_mesh_filepaths,
- OUT char*** xml_mesh_names,
- OUT char*** xml_mesh_materials
-)
-{
- FILE* f = fopen(filepath, "r");
- if (f == NULL)
- PERROR_CLEANUP_EXIT("Error: cannot open the xml file", 8);
-
- fseek(f, 0, SEEK_END);
- long fsize = ftell(f);
- fseek(f, 0, SEEK_SET);
-
- char* buf = (char*)malloc(fsize + 1);
- if (!buf)
- PERROR_CLEANUP_EXIT("Error: cannot allocate memory for the xml file", 8);
- if (fread(buf, 1, fsize, f) != fsize)
- PERROR_CLEANUP_EXIT("Error: cannot read the xml file", 8, buf);
- buf[fsize] = '\0';
-
- /* Count the amount of "<shape" entries */
- char** shape_positions = NULL;
- *xml_num_meshes = 0;
- char* pos = buf, *end;
- while ((pos = strstr(pos, "<shape")) != NULL) {
- ++*xml_num_meshes;
- shape_positions = (char**)realloc(shape_positions, *xml_num_meshes * sizeof(char*));
- if (!shape_positions)
- PERROR_CLEANUP_EXIT("Error: cannot reallocate memory for shape positions", 8, buf);
- shape_positions[*xml_num_meshes - 1] = pos;
- pos += 6;
- }
- if (*xml_num_meshes == 0)
- PERROR_CLEANUP_EXIT("Error: no shapes found in the xml file", 8, buf);
- if (!shape_positions)
- PERROR_CLEANUP_EXIT("Error: cannot allocate memory for shape positions", 8, buf);
-
- /* Allocate memory for the mesh data */
- *xml_mesh_filepaths = (char**)malloc(*xml_num_meshes * sizeof(char*));
- *xml_mesh_names = (char**)malloc(*xml_num_meshes * sizeof(char*));
- *xml_mesh_materials = (char**)malloc(*xml_num_meshes * sizeof(char*));
- if (!*xml_mesh_filepaths || !*xml_mesh_names || !*xml_mesh_materials)
- PERROR_CLEANUP_EXIT("Error: cannot allocate memory for mesh data", 8, buf);
-
- /* For each "<shape" entry, extract the name, file path and material */
- char *e_name, *e_filepath, *e_material;
- for (uint32_t i = 0; i != *xml_num_meshes; ++i) {
- pos = shape_positions[i] + 6;
-
- /* Find mesh name */
- if ((pos = strstr(pos, "name=\"")) == NULL)
- PERROR_CLEANUP_EXIT("Error: cannot find mesh name", 8,
- buf, shape_positions, *xml_mesh_filepaths, *xml_mesh_names, *xml_mesh_materials);
- pos += 6;
- end = strchr(pos, '\"');
- if (!end)
- PERROR_CLEANUP_EXIT("Error: cannot find mesh name end", 8,
- buf, shape_positions, *xml_mesh_filepaths, *xml_mesh_names, *xml_mesh_materials);
- e_name = (char*)malloc(end - pos + 1);
- if (!e_name)
- PERROR_CLEANUP_EXIT("Error: cannot allocate memory for mesh name", 8,
- buf, shape_positions, *xml_mesh_filepaths, *xml_mesh_names, *xml_mesh_materials);
- strncpy(e_name, pos, end - pos);
- e_name[end - pos] = '\0';
-
- /* Find mesh file path */
- if (!(pos = strstr(pos, "<string name=\"filename\"")))
- PERROR_CLEANUP_EXIT("Error: cannot find mesh file path", 8, buf, e_name,
- shape_positions, *xml_mesh_filepaths, *xml_mesh_names, *xml_mesh_materials);
- if (!(pos = strstr(pos, "value=\"")))
- PERROR_CLEANUP_EXIT("Error: cannot find mesh file path value", 8, buf, e_name,
- shape_positions, *xml_mesh_filepaths, *xml_mesh_names, *xml_mesh_materials);
- pos += 7;
- end = strchr(pos, '\"');
- if (!end)
- PERROR_CLEANUP_EXIT("Error: cannot find mesh file path end", 8, buf, e_name,
- shape_positions, *xml_mesh_filepaths, *xml_mesh_names, *xml_mesh_materials);
- e_filepath = (char*)malloc(end - pos + 1);
- if (!e_filepath)
- PERROR_CLEANUP_EXIT("Error: cannot allocate memory for mesh file path", 8, buf, e_name,
- shape_positions, *xml_mesh_filepaths, *xml_mesh_names, *xml_mesh_materials);
- strncpy(e_filepath, pos, end - pos);
- e_filepath[end - pos] = '\0';
-
- /* Find mesh material */
- if (!(pos = strstr(pos, "id=\"mat-itu_")))
- PERROR_CLEANUP_EXIT("Error: cannot find mesh material", 8, buf, e_name, e_filepath,
- shape_positions, *xml_mesh_filepaths, *xml_mesh_names, *xml_mesh_materials);
- pos += 12;
- end = strchr(pos, '\"');
- if (!end)
- PERROR_CLEANUP_EXIT("Error: cannot find mesh material end", 8, buf, e_name, e_filepath,
- shape_positions, *xml_mesh_filepaths, *xml_mesh_names, *xml_mesh_materials);
- e_material = (char*)malloc(end - pos + 1);
- if (!e_material)
- PERROR_CLEANUP_EXIT("Error: cannot allocate memory for mesh material", 8, buf, e_name, e_filepath,
- shape_positions, *xml_mesh_filepaths, *xml_mesh_names, *xml_mesh_materials);
- strncpy(e_material, pos, end - pos);
- e_material[end - pos] = '\0';
-
- /* Save mesh data */
- (*xml_mesh_filepaths)[i] = e_filepath;
- (*xml_mesh_names)[i] = e_name;
- (*xml_mesh_materials)[i] = e_material;
- }
-
- free(shape_positions);
- free(buf);
-}
-
-/* Read a Sionna .xml scene. Assumes meshes are in a "meshes" directory.
- *
- * \param filepath Path to the .xml scene file
- * \param scene The scene to fill
- */
-void readScene(const char* filepath, HRT_Scene* scene) {
- /* Check if filepath ends with ".xml" */
- size_t filepath_len = strlen(filepath);
- if (filepath_len > 4 && strcmp(filepath + filepath_len - 4, ".xml") != 0)
- PERROR_CLEANUP_EXIT("Error: scene file must end with .xml", 8);
-
- /* Read the xml file */
- uint32_t xml_num_meshes;
- char** xml_mesh_filepaths;
- char** xml_mesh_names;
- char** xml_mesh_materials;
- readXml(
- filepath,
- &xml_num_meshes,
- &xml_mesh_filepaths,
- &xml_mesh_names,
- &xml_mesh_materials
- );
-
- /* Read the scene config CSV file */
- char* csv_path = (char*)malloc(filepath_len + 1);
- if (!csv_path)
- PERROR_CLEANUP_EXIT("Error: cannot allocate memory for csv_path", 8);
- strcpy(csv_path, filepath);
- strcpy(csv_path + filepath_len - 4, ".csv");
- uint32_t csv_num_meshes;
- char** csv_mesh_filenames;
- uint32_t* csv_material_indicies;
- float* csv_velocity;
- readCsv(
- csv_path,
- &csv_num_meshes,
- csv_mesh_filenames,
- csv_material_indicies,
- csv_velocity
- );
- free(csv_path);
-
- /* Read the meshes */
-
- /* Get the scene directory path */
- size_t scene_dir_len = strrchr(filepath, '/') - filepath + 1;
- /* if filepath is "/path/to/scene.xml", scene_dir is "/path/to/" */
- char* scene_dir = (char*)malloc(scene_dir_len + 7);
- if (!scene_dir)
- PERROR_CLEANUP_EXIT("Error: failed to allocate scene_dir", 8);
- strncpy(scene_dir, filepath, scene_dir_len);
-
- /* For each ply file from the xml file */
- scene->num_meshes = xml_num_meshes;
- scene->meshes = (HRT_Mesh*)malloc(xml_num_meshes * sizeof(HRT_Mesh));
- if (!scene->meshes && xml_num_meshes > 0)
- PERROR_CLEANUP_EXIT("Error: failed to allocate meshes array", 8);
- for (uint32_t i = 0; i != xml_num_meshes; ++i) {
- /* Get the mesh file path */
- size_t mesh_filepath_len = scene_dir_len + strlen(xml_mesh_filepaths[i]);
- char* mesh_filepath = (char*)malloc(mesh_filepath_len);
- mesh_filepath[0] = '\0';
- strncat(mesh_filepath, scene_dir, scene_dir_len);
- strncat(mesh_filepath + scene_dir_len, xml_mesh_filepaths[i], mesh_filepath_len - scene_dir_len);
- /* Read the mesh */
- scene->meshes[i] = readPly(mesh_filepath);
- free(mesh_filepath);
- /* Set the material index from the material name from xml */
- scene->meshes[i].material_index = get_material_index(xml_mesh_materials[i]);
- /* Set the material index and velocity from CSV if present */
- for (uint32_t j = 0; j != csv_num_meshes; ++j)
- if (strcmp(xml_mesh_names[i], csv_mesh_filenames[j]) == 0) {
- scene->meshes[i].material_index = csv_material_indicies[j];
- scene->meshes[i].velocity[0] = csv_velocity[j * 3];
- scene->meshes[i].velocity[1] = csv_velocity[j * 3 + 1];
- scene->meshes[i].velocity[2] = csv_velocity[j * 3 + 2];
- break;
- }
- }
-
- free(scene_dir);
-}
-
-
-/**
- * Main
- */
-
-int main(int argc, char *argv[]) {
- if (argc != 2) {
- printf("Usage: %s <scene.xml>\n", argv[0]);
- return 1;
- }
- const char* scene_filepath = argv[1];
- const char* scene_filename = strrchr(scene_filepath, '/');
- if (scene_filename == NULL)
- exit(8);
- else
- ++scene_filename;
-
- HRT_Scene* scene_ptr;
-
- /* read scene */
- /* check if the scene is a hardcoded scene */
- if (!strcmp(scene_filename, scene_box_filename))
- scene_ptr = &scene_box;
- else if (!strcmp(scene_filename, scene_simpleReflector_filename))
- scene_ptr = &scene_simpleReflector;
- else {
- scene_ptr = (HRT_Scene*)malloc(sizeof(HRT_Scene));
- readScene(scene_filepath, scene_ptr);
- }
-
- FILE *fp = fopen("scene.hrt", "wb");
- if (!fp)
- PERROR_CLEANUP_EXIT("Error: cannot open file", 8);
-
- /* MAGIC */
- fwrite("HRT", 1, 3, fp);
- /* SCENE */
- /* num_meshes */
- fwrite(&scene_ptr->num_meshes, sizeof(uint32_t), 1, fp);
- /* meshes */
- for (uint32_t i = 0; i != scene_ptr->num_meshes; ++i) {
- /* MESH */
- fwrite(&scene_ptr->meshes[i].num_vertices, sizeof(uint32_t), 1, fp);
- fwrite(scene_ptr->meshes[i].vs, sizeof(float), 3 * scene_ptr->meshes[i].num_vertices, fp);
- fwrite(&scene_ptr->meshes[i].num_triangles, sizeof(uint32_t), 1, fp);
- fwrite(scene_ptr->meshes[i].is, sizeof(uint32_t), 3 * scene_ptr->meshes[i].num_triangles, fp);
- fwrite(&scene_ptr->meshes[i].material_index, sizeof(uint32_t), 1, fp);
- fwrite(&scene_ptr->meshes[i].velocity, sizeof(float), 3, fp);
- }
-
- fclose(fp);
- return 0;
-}
diff --git a/scene_fromSionna.elf b/scene_fromSionna.elf
Binary files differ.
diff --git a/src/compute_paths.c b/src/compute_paths.c
@@ -0,0 +1,799 @@
+#include "../inc/compute_paths.h" /* for compute_paths */
+#include "../inc/scene.h" /* for Scene, Mesh, Material */
+#include "../inc/materials.h" /* for g_materials */
+
+#include <stddef.h> /* for size_t */
+#include <stdlib.h> /* for exit, malloc, free */
+#include <string.h> /* for strlen, sprintf */
+#include <stdio.h> /* for fopen, FILE, fclose */
+#include <math.h> /* for sin, cos, sqrt */
+#include <stdint.h> /* for int32_t, uint8_t */
+
+/* ==== CONSTANTS ==== */
+
+#define PI 3.14159265358979323846f /* pi */
+#define SPEED_OF_LIGHT 299792458.0f /* m/s */
+
+/* (2w, w) binomial coefficients for w = 0, 1, ..., 19 */
+const float BINOMIAL_2W_W[20] = {
+ 1.f, 2.f, 6.f, 20.f, 70.f,
+ 252.f, 924.f, 3432.f, 12870.f, 48620.f,
+ 184756.f, 705432.f, 2704156.f, 10400600.f, 40116600.f,
+ 155117520.f, 601080390.f, 2333606220.f, 9075135300.f, 35345263800.f
+};
+/* (aplha, k) binomial coefficients for alpha = 0, 1, ..., 19 with k = 0, 1, ..., alpha */
+const float* BINOMIAL_ALPHA_K[20] = {
+ (float[]){
+ 1.f
+ },
+ (float[]){
+ 1.f, 1.f
+ },
+ (float[]){
+ 1.f, 2.f, 1.f
+ },
+ (float[]){
+ 1.f, 3.f, 3.f, 1.f
+ },
+ (float[]){
+ 1.f, 4.f, 6.f, 4.f, 1.f
+ },
+ (float[]){
+ 1.f, 5.f, 10.f, 10.f, 5.f,
+ 1.f
+ },
+ (float[]){
+ 1.f, 6.f, 15.f, 20.f, 15.f,
+ 6.f, 1.f
+ },
+ (float[]){
+ 1.f, 7.f, 21.f, 35.f, 35.f,
+ 21.f, 7.f, 1.f
+ },
+ (float[]){
+ 1.f, 8.f, 28.f, 56.f, 70.f,
+ 56.f, 28.f, 8.f, 1.f
+ },
+ (float[]){
+ 1.f, 9.f, 36.f, 84.f, 126.f,
+ 126.f, 84.f, 36.f, 9.f, 1.f
+ },
+ (float[]){
+ 1.f, 10.f, 45.f, 120.f, 210.f,
+ 252.f, 210.f, 120.f, 45.f, 10.f,
+ 1.f
+ },
+ (float[]){
+ 1.f, 11.f, 55.f, 165.f, 330.f,
+ 462.f, 462.f, 330.f, 165.f, 55.f,
+ 11.f, 1.f
+ },
+ (float[]){
+ 1.f, 12.f, 66.f, 220.f, 495.f,
+ 792.f, 924.f, 792.f, 495.f, 220.f,
+ 66.f, 12.f, 1.f
+ },
+ (float[]){
+ 1.f, 13.f, 78.f, 286.f,
+ 715.f, 1287.f, 1716.f, 1716.f, 1287.f,
+ 715.f, 286.f, 78.f, 13.f, 1.f
+ },
+ (float[]){
+ 1.f, 14.f, 91.f, 364.f, 1001.f,
+ 2002.f, 3003.f, 3432.f, 3003.f, 2002.f,
+ 1001.f, 364.f, 91.f, 14.f, 1.f
+ },
+ (float[]){
+ 1.f, 15.f, 105.f, 455.f, 1365.f,
+ 3003.f, 5005.f, 6435.f, 6435.f, 5005.f,
+ 3003.f, 1365.f, 455.f, 105.f, 15.f,
+ 1.f
+ },
+ (float[]){
+ 1.f, 16.f, 120.f, 560.f, 1820.f,
+ 4368.f, 8008.f, 11440.f, 12870.f, 11440.f,
+ 8008.f, 4368.f, 1820.f, 560.f, 120.f,
+ 16.f, 1.f
+ },
+ (float[]){
+ 1.f, 17.f, 136.f, 680.f, 2380.f,
+ 6188.f, 12376.f, 19448.f, 24310.f, 24310.f,
+ 19448.f, 12376.f, 6188.f, 2380.f, 680.f,
+ 136.f, 17.f, 1.f
+ },
+ (float[]){
+ 1.f, 18.f, 153.f, 816.f, 3060.f,
+ 8568.f, 18564.f, 31824.f, 43758.f, 48620.f,
+ 43758.f, 31824.f, 18564.f, 8568.f, 3060.f,
+ 816.f, 153.f, 18.f, 1.f
+ },
+ (float[]){
+ 1.f, 19.f, 171.f, 969.f, 3876.f,
+ 11628.f, 27132.f, 50388.f, 75582.f, 92378.f,
+ 92378.f, 75582.f, 50388.f, 27132.f, 11628.f,
+ 3876.f, 969.f, 171.f, 19.f, 1.f
+ }
+};
+
+/* ==== STRUCTS ==== */
+
+typedef struct {
+ Vec3 o, d;
+} Ray;
+
+/* Radio material eta and additional parameters */
+typedef struct {
+ /* eta */
+ float eta_re, eta_sqrt_re, eta_inv_re, eta_inv_sqrt_re;
+ float eta_im, eta_sqrt_im, eta_inv_im, eta_inv_sqrt_im;
+ float eta_abs, eta_abs_pow2, eta_abs_inv_sqrt;
+ /* r = 1.0 - s */
+ float r;
+} MaterialPrecomputed;
+
+/* ==== COMPLEX OPERATIONS ==== */
+
+void csqrtf(
+ IN float re,
+ IN float im,
+ IN float abs,
+ OUT float *sqrt_re,
+ OUT float *sqrt_im
+)
+{
+ *sqrt_re = sqrtf((re + abs) / 2.f);
+ if (fabsf(im) < __FLT_EPSILON__ && re >= -__FLT_EPSILON__)
+ *sqrt_im = 0;
+ else {
+ *sqrt_im = sqrtf((abs - re) / 2.f);
+ if (im < 0.f) *sqrt_im = -*sqrt_im;
+ }
+}
+void cdivf(
+ IN float a_re,
+ IN float a_im,
+ IN float b_re,
+ IN float b_im,
+ OUT float *c_re,
+ OUT float *c_im
+)
+{
+ float d = b_re * b_re + b_im * b_im;
+ *c_re = (a_re * b_re + a_im * b_im) / d;
+ *c_im = (a_im * b_re - a_re * b_im) / d;
+}
+
+/* ==== PRECOMPUTED GLOBALS ==== */
+
+/* 4.f * PI * carrier_frequency * 1e9 / SPEED_OF_LIGHT */
+float g_free_space_loss_multiplier;
+/* Materials etas, calucalted for the given carrier_frequency */
+MaterialPrecomputed g_materials_precomputed[NUM_G_MATERIALS];
+
+void precompute_free_space_loss_multiplier(
+ IN float carrier_frequency
+)
+{
+ g_free_space_loss_multiplier = 4.f * PI * carrier_frequency * 1e9 / SPEED_OF_LIGHT;
+}
+
+void precompute_materials(
+ IN Scene* scene,
+ IN float carrier_frequency
+)
+{
+ uint8_t material_eta_isdone[NUM_G_MATERIALS] = {0};
+ for (uint32_t i = 0; i != scene->num_meshes; ++i) {
+ uint32_t mat_ind = scene->meshes[i].material_index;
+ if (material_eta_isdone[mat_ind])
+ continue;
+ Material *mat = &g_materials[mat_ind];
+ MaterialPrecomputed *mat_precomp = &g_materials_precomputed[mat_ind];
+ /* calculate eta */
+ mat_precomp->eta_re = mat->a * powf(carrier_frequency, mat->b);
+ /* eq. 12 */
+ mat_precomp->eta_im = (mat->c * powf(carrier_frequency, mat->d))
+ / (0.0556325027352135f * carrier_frequency);
+ mat_precomp->eta_abs_pow2 = mat_precomp->eta_re * mat_precomp->eta_re
+ + mat_precomp->eta_im * mat_precomp->eta_im;
+ mat_precomp->eta_abs = sqrtf(mat_precomp->eta_abs_pow2);
+ mat_precomp->eta_abs_inv_sqrt = 1.f / sqrtf(mat_precomp->eta_abs);
+ /* calculate sqrt(eta) */
+ csqrtf(mat_precomp->eta_re, mat_precomp->eta_im,
+ mat_precomp->eta_abs, &mat_precomp->eta_sqrt_re,
+ &mat_precomp->eta_sqrt_im);
+ /* calculate 1 / eta */
+ mat_precomp->eta_inv_re = mat_precomp->eta_re / mat_precomp->eta_abs_pow2;
+ mat_precomp->eta_inv_im = -mat_precomp->eta_im / mat_precomp->eta_abs_pow2;
+ /* calculate 1 / sqrt(eta) */
+ csqrtf(mat_precomp->eta_inv_re, mat_precomp->eta_inv_im,
+ 1.f / mat_precomp->eta_abs,
+ &mat_precomp->eta_inv_sqrt_re, &mat_precomp->eta_inv_sqrt_im);
+ /* calculate r */
+ mat_precomp->r = 1.f - mat->s;
+ }
+}
+
+void precompute_normals(Scene* scene) {
+ /* Calculate normals */
+ for (uint32_t i = 0; i != scene->num_meshes; ++i) {
+ Mesh *mesh = &scene->meshes[i];
+ mesh->ns = (Vec3*)malloc(mesh->num_triangles * sizeof(Vec3));
+ for (uint32_t j = 0; j != mesh->num_triangles; ++j) {
+ Vec3* v1 = (Vec3*)&mesh->vs[mesh->is[j * 3]];
+ Vec3* v2 = (Vec3*)&mesh->vs[mesh->is[j * 3 + 1]];
+ Vec3* v3 = (Vec3*)&mesh->vs[mesh->is[j * 3 + 2]];
+ Vec3 e1 = vec3_sub(v2, v1);
+ Vec3 e2 = vec3_sub(v3, v1);
+ Vec3 n = vec3_cross(&e1, &e2);
+ n = vec3_normalize(&n);
+ memcpy(&mesh->ns[j], &n, sizeof(Vec3));
+ }
+ }
+}
+
+/* ==== RT ==== */
+
+/** Compute Moeleer-Trumbore intersection algorithm.
+ *
+ * \param ray ray to cast to the mesh
+ * \param scene the scene to cast the ray to
+ * \param t output distance to the hit point. If no hit, t is not modified.
+ * \param mesh_ind output index of the hit mesh. If no hit, i is not modified.
+ * \param face_ind output index of the hit triangle. If no hit, i is not modified.
+ * \param theta output angle of incidence
+ */
+void moeller_trumbore(
+ IN Ray *ray,
+ IN Scene *scene,
+ OUT float *t,
+ OUT uint32_t *mesh_ind,
+ OUT uint32_t *face_ind,
+ OUT float *theta
+)
+{
+ /* TODO BVH */
+ /* for each triangle */
+ Vec3 *v1, *v2, *v3, *n;
+ Vec3 e1, e2, re2_cross, s, se1_cross;
+ float d, u, v;
+ float dist = 1e9;
+ float dist_tmp;
+ for (uint32_t i = 0; i != scene->num_meshes; ++i) {
+ Mesh *mesh = &scene->meshes[i];
+ for (uint32_t j = 0; j != mesh->num_triangles; ++j) {
+ v1 = (Vec3*)&mesh->vs[mesh->is[3 * j]];
+ v2 = (Vec3*)&mesh->vs[mesh->is[3 * j + 1]];
+ v3 = (Vec3*)&mesh->vs[mesh->is[3 * j + 2]];
+ e1 = vec3_sub(v2, v1);
+ e2 = vec3_sub(v3, v1);
+ re2_cross = vec3_cross(&ray->d, &e2);
+ d = vec3_dot(&e1, &re2_cross);
+ if (d > -__FLT_EPSILON__ && d < __FLT_EPSILON__) continue;
+ s = vec3_sub(&ray->o, v1);
+ u = vec3_dot(&s, &re2_cross) / d;
+ if ((u < 0. && fabs(u) > __FLT_EPSILON__)
+ || (u > 1. && fabs(u - 1.) > __FLT_EPSILON__))
+ continue;
+ se1_cross = vec3_cross(&s, &e1);
+ v = vec3_dot(&ray->d, &se1_cross) / d;
+ if ((v < 0. && fabs(v) > __FLT_EPSILON__)
+ || (u + v > 1. && fabs(u + v - 1.) > __FLT_EPSILON__))
+ continue;
+ dist_tmp = vec3_dot(&e2, &se1_cross) / d;
+ if (dist_tmp > __FLT_EPSILON__ && dist_tmp < dist) {
+ dist = dist_tmp;
+ *t = dist;
+ *mesh_ind = i;
+ *face_ind = j;
+ n = (Vec3*)&mesh->ns[j];
+ *theta = acos(vec3_dot(n, &ray->d));
+ if (*theta > PI / 2.)
+ *theta = PI - *theta;
+ }
+ }
+ }
+}
+
+/** Compute reflection R_{eTE} and R_{eTM} coefficients.
+ *
+ * Implements eqs. (31a)-(31b) from ITU-R P.2040-3.
+ *
+ * \param material_index the index of the material
+ * \param theta1 the angle of incidence
+ * \param r_te_re output real part of R_{eTE}
+ * \param r_te_im output imaginary part of R_{eTE}
+ * \param r_tm_re output real part of R_{eTM}
+ * \param r_tm_im output imaginary part of R_{eTM}
+ */
+void refl_coefs(
+ IN uint32_t material_index,
+ IN float theta1,
+ OUT float *r_te_re,
+ OUT float *r_te_im,
+ OUT float *r_tm_re,
+ OUT float *r_tm_im
+)
+{
+ MaterialPrecomputed *rm = &g_materials_precomputed[material_index];
+ float sin_theta1 = sinf(theta1);
+ if (rm->eta_abs_inv_sqrt * sin_theta1 > 1.f - __FLT_EPSILON__) {
+ *r_te_re = *r_tm_re = 1.f;
+ *r_te_im = *r_tm_im = 0.f;
+ return;
+ }
+
+ /* eq. 33 */
+ float sin_theta1_pow2 = sin_theta1 * sin_theta1;
+ float cos_theta2_re = sqrtf(1.f + rm->eta_inv_re / rm->eta_abs_pow2 * sin_theta1_pow2);
+ float cos_theta2_im = sqrtf(1.f - rm->eta_inv_im / rm->eta_abs_pow2 * sin_theta1_pow2);
+
+ /* R_{eTE}, eq. 31a */
+ float sqrt_eta_cos_theta2_re = rm->eta_sqrt_re * cos_theta2_re - rm->eta_sqrt_im * cos_theta2_im;
+ float sqrt_eta_cos_theta2_im = rm->eta_sqrt_re * cos_theta2_im + rm->eta_sqrt_im * cos_theta2_re;
+ float cos_theta1 = cosf(theta1);
+ cdivf(cos_theta1 - sqrt_eta_cos_theta2_re, -sqrt_eta_cos_theta2_im,
+ cos_theta1 + sqrt_eta_cos_theta2_re, sqrt_eta_cos_theta2_im,
+ r_te_re, r_te_im);
+
+ /* R_{eTM}, eq. 31b */
+ float sqrt_eta_cos_theta1_re = rm->eta_sqrt_re * cos_theta1;
+ float sqrt_eta_cos_theta1_im = rm->eta_sqrt_im * cos_theta1;
+ cdivf(sqrt_eta_cos_theta1_re - cos_theta2_re,
+ sqrt_eta_cos_theta1_im - cos_theta2_im,
+ sqrt_eta_cos_theta1_re + cos_theta2_re,
+ sqrt_eta_cos_theta1_im + cos_theta2_im,
+ r_tm_re, r_tm_im);
+
+ /* Apply reflection reduction factor */
+ *r_te_re *= rm->r;
+ *r_te_im *= rm->r;
+ *r_tm_re *= rm->r;
+ *r_tm_im *= rm->r;
+}
+
+/** Calculate scattering coefficients.
+ *
+ * Implements a directive scattering model inspired by Blaunstein et al. (DOI 10.1109/TAP.2006.888422).
+ * Models scattering from rough surfaces (e.g., vegetation) with polarization dependence.
+ *
+ * \param theta_s scattering angle (radians, between scattered direction and surface normal)
+ * \param theta_i angle of incidence (radians, between incident ray and surface normal)
+ * \param material_index index into g_materials array
+ * \param a_te_re output real part of TE scattering coefficient
+ * \param a_te_im output imaginary part of TE scattering coefficient
+ * \param a_tm_re output real part of TM scattering coefficient
+ * \param a_tm_im output imaginary part of TM scattering coefficient
+ */
+void scat_coefs(
+ IN float theta_s,
+ IN float theta_i,
+ IN uint32_t material_index,
+ OUT float *a_te_re,
+ OUT float *a_te_im,
+ OUT float *a_tm_re,
+ OUT float *a_tm_im
+)
+{
+ Material *mat = &g_materials[material_index];
+
+ /* Precompute trigonometric values */
+ float cos_theta_s = cosf(theta_s);
+ float cos_theta_i = cosf(theta_i);
+ float sin_theta_i = sinf(theta_i);
+
+ /* Scattering directivity pattern: exponential decay with angle difference */
+ /* f = s * exp(-alpha * |theta_s - theta_i|) models directive scattering */
+ float theta_diff = fabsf(theta_s - theta_i);
+ float f = mat->s * expf(-mat->s1_alpha * theta_diff);
+
+ /* Roughness factor: reduces specular component as alpha increases */
+ float roughness = 1.0f / (1.0f + mat->s1_alpha); // Simplified, 0 to 1
+ float specular = roughness * cos_theta_s; // Specular-like term
+ float diffuse = (1.0f - roughness) * cos_theta_s; // Diffuse term
+
+ /* TE and TM coefficients (real parts) */
+ /* TE: perpendicular to plane of incidence, less affected by angle */
+ float te_real = f * (specular + diffuse);
+ /* TM: parallel to plane of incidence, stronger angular dependence */
+ float tm_real = f * (specular * cos_theta_i + diffuse);
+
+ /* Phase shift (imaginary parts) due to surface roughness or path difference */
+ /* Simplified: assume small phase shift proportional to roughness and angle */
+ float phase_shift = mat->s1_alpha * sin_theta_i * 0.1f; // Arbitrary scaling
+ float te_imag = te_real * sinf(phase_shift);
+ float tm_imag = tm_real * sinf(phase_shift);
+
+ /* Normalize to ensure energy conservation (optional, adjust as needed) */
+ float norm = sqrtf(te_real * te_real + te_imag * te_imag +
+ tm_real * tm_real + tm_imag * tm_imag);
+ if (norm > 1e-6f) { // Avoid division by zero
+ te_real /= norm;
+ te_imag /= norm;
+ tm_real /= norm;
+ tm_imag /= norm;
+ }
+
+ /* Output coefficients */
+ *a_te_re = te_real;
+ *a_te_im = te_imag;
+ *a_tm_re = tm_real;
+ *a_tm_im = tm_imag;
+
+ /* TODO: Incorporate s2 and s3 for additional roughness or frequency effects */
+}
+
+/* ==== MAIN FUNCTION ==== */
+
+void compute_paths(
+ IN Scene *scene, /* Pointer to a loaded scene */
+ IN const float *rx_pos, /* shape (num_rx, 3) */
+ IN const float *tx_pos, /* shape (num_tx, 3) */
+ IN const float *rx_vel, /* shape (num_rx, 3) */
+ IN const float *tx_vel, /* shape (num_tx, 3) */
+ IN float carrier_frequency, /* > 0.0 (IN GHz!) */
+ IN size_t num_rx, /* number of receivers */
+ IN size_t num_tx, /* number of transmitters */
+ IN size_t num_paths, /* number of paths */
+ IN size_t num_bounces, /* number of bounces */
+ OUT PathsInfo *los, /* output LoS information. los->num_paths = 1 */
+ OUT PathsInfo *scatter /* output scatter information. scatter->num_paths = num_bounces*num_paths */
+)
+{
+ /* Precompute globals */
+ precompute_free_space_loss_multiplier(carrier_frequency);
+ precompute_materials(scene, carrier_frequency);
+ precompute_normals(scene);
+
+ /* Calculate a fibonacci sphere */
+ Vec3 *ray_directions = (Vec3*)malloc(num_paths * sizeof(Vec3));
+ float k, phi, theta;
+ for (size_t path = 0; path < num_paths; ++path) {
+ k = (float)path + .5f;
+ phi = acos(1.f - 2.f * k / num_paths);
+ theta = PI * (1.f + sqrtf(5.f)) * k;
+ ray_directions[path] = (Vec3){
+ cos(theta) * sin(phi),
+ sin(theta) * sin(phi),
+ cos(phi)
+ };
+ }
+
+ /* Record the tx directions */
+ /* This is only valid for the fibbonaci sphere */
+ scatter->directions_tx = (float*)memcpy(
+ scatter->directions_tx,
+ ray_directions,
+ num_paths * sizeof(Vec3)
+ );
+ if (scatter->directions_tx == NULL) {
+ fprintf(stderr, "Failed to copy memory to scatter->directions_tx\n");
+ exit(70);
+ }
+
+ /* Cast the positions and velocities to Vec3 */
+ Vec3 *rx_pos_v = (Vec3*)rx_pos;
+ Vec3 *tx_pos_v = (Vec3*)tx_pos;
+ Vec3 *rx_vel_v = (Vec3*)rx_vel;
+ Vec3 *tx_vel_v = (Vec3*)tx_vel;
+
+ /* Create num_path rays for each tx */
+ /* Shape (num_tx, num_paths) */
+ Ray *rays = (Ray*)malloc(num_tx * num_paths * sizeof(Ray));
+ for (size_t tx = 0, off = 0; tx < num_tx; ++tx) {
+ for (size_t path = 0; path < num_paths; ++path, ++off) {
+ rays[off].o = tx_pos_v[tx];
+ rays[off].d = ray_directions[path];
+ }
+ }
+ free(ray_directions);
+
+ /* Initialize variables needed for the RT */
+
+ /* Bouncing (reflecting) rays gains and delays */
+ /* shape (num_tx, num_paths) */
+ float *a_te_re_refl = (float*)malloc(num_tx * num_paths * sizeof(float));
+ float *a_te_im_refl = (float*)calloc(num_tx * num_paths, sizeof(float));
+ float *a_tm_re_refl = (float*)malloc(num_tx * num_paths * sizeof(float));
+ float *a_tm_im_refl = (float*)calloc(num_tx * num_paths, sizeof(float));
+ for (size_t i = 0; i < num_tx * num_paths; ++i)
+ a_te_re_refl[i] = a_tm_re_refl[i] = 1.f;
+ float *tau_t = (float*)calloc(num_tx * num_paths, sizeof(float));
+
+ /* Active rays bitmask. 1 if the ray is still traced, 0 if it has left the scene */
+ uint8_t *active = (uint8_t*)malloc((num_tx * num_paths / 8 + 1) * sizeof(uint8_t));
+ for (size_t i = 0; i < num_tx * num_paths / 8 + 1; ++i)
+ active[i] = 0xff;
+ size_t num_active = num_tx * num_paths;
+
+ /* Temporary variables */
+ float t, r_te_re, r_te_im, r_tm_re, r_tm_im;
+ float a_te_re_new, a_te_im_new, a_tm_re_new, a_tm_im_new;
+ float theta_s; /* scattering angle */
+ uint32_t mesh_ind, face_ind, material_index;
+ Ray r;
+ Vec3 *h = (Vec3*)malloc(sizeof(Vec3));
+ Vec3 n;
+
+ /* Friis free space loss multiplier.
+ Multiply this by a distance and take the second power to get a free space loss. */
+ float free_space_loss_multiplier = 4.f * PI * carrier_frequency * 1e9 / SPEED_OF_LIGHT;
+ float free_space_loss;
+
+ /* Doppler frequency shift carrier frequency multiplier */
+ float doppler_shift_multiplier = carrier_frequency * 1e9 / SPEED_OF_LIGHT;
+
+ /* Initialize the doppler frequency shift using the tx velocities */
+ /* The scatter->freq_shift has the shape of (num_rx, num_tx, num_bounces*num_paths) */
+ /* The first num_tx*num_paths elements are the same for each rx and bounce */
+ /* so we can initialize them only once and memcpy them to the rest */
+ for (size_t tx = 0; tx != num_tx; ++tx)
+ for (size_t path = 0; path != num_paths; ++path) {
+ size_t off_rays = tx * num_paths + path;
+ size_t off_scat = tx * num_paths * num_bounces + path;
+ scatter->freq_shift[off_scat] = vec3_dot(&tx_vel_v[tx], &rays[off_rays].d);
+ scatter->freq_shift[off_scat] *= doppler_shift_multiplier;
+ }
+ for (size_t bounce = 1; bounce < num_bounces; ++bounce)
+ if (!memcpy(scatter->freq_shift + num_tx * num_paths * bounce,
+ scatter->freq_shift, num_tx * num_paths * sizeof(float)))
+ exit(70);
+ for (size_t rx = 1; rx < num_rx; ++rx)
+ if (!memcpy(scatter->freq_shift + num_tx * num_paths * num_bounces * rx,
+ scatter->freq_shift, num_tx * num_paths * num_bounces * sizeof(float)))
+ exit(70);
+
+ /***************************************************************************/
+ /* LoS paths */
+ /***************************************************************************/
+
+ /* Consider imaginary parts of the LoS gains to be 0 in any way */
+ for (size_t off = 0; off != num_rx * num_tx; ++off)
+ los->a_te_im[off] = los->a_tm_im[off] = 0.f;
+
+ /* Calculate LoS paths directly from tx to rx */
+ for (size_t rx = 0, off = 0; rx != num_rx; ++rx)
+ for (size_t tx = 0; tx != num_tx; ++tx, ++off) {
+ t = -1.f;
+ mesh_ind = face_ind = -1;
+
+ /* Create a ray from the tx to the rx */
+ r.o = tx_pos_v[tx];
+ r.d = vec3_sub(&rx_pos_v[rx], &r.o);
+
+ /* In case the tx and the rx are at the same position */
+ if (vec3_dot(&r.d, &r.d) < __FLT_EPSILON__) {
+ /* Assume the direction to be (1, 0, 0) */
+ los->directions_rx[off * 3] = 1.f;
+ los->directions_tx[off * 3] = -1.f;
+ los->directions_rx[off * 3 + 1] = los->directions_tx[off * 3 + 1] = 0.f;
+ los->directions_rx[off * 3 + 2] = los->directions_tx[off * 3 + 2] = 0.f;
+ /* Set gains to 1+0j */
+ los->a_te_re[off] = los->a_tm_re[off] = 1.f;
+ /* Set delay to 0 */
+ los->tau[off] = 0.f;
+ /* Set doppler shift to 0. */
+ los->freq_shift[off] = 0.f;
+ continue;
+ }
+
+ /* Is there any abstacle between the tx and the rx? */
+ moeller_trumbore(&r, scene, &t, &mesh_ind, &face_ind, &theta);
+ if (mesh_ind != (uint32_t)-1 && t <= 1.f) {
+ /* An obstacle between the tx and the rx has been hit */
+ los->a_te_re[off] = los->a_tm_re[off] = los->tau[off] = 0.f;
+ continue;
+ }
+
+ /* No obstacle has been hit, the path is LoS */
+ /* Calculate the distance */
+ t = sqrtf(vec3_dot(&r.d, &r.d));
+ /* Calculate the direction */
+ Vec3 d = {r.d.x / t, r.d.y / t, r.d.z / t};
+ los->directions_tx[off * 3] = d.x;
+ los->directions_tx[off * 3 + 1] = d.y;
+ los->directions_tx[off * 3 + 2] = d.z;
+ los->directions_rx[off * 3] = -d.x;
+ los->directions_rx[off * 3 + 1] = -d.y;
+ los->directions_rx[off * 3 + 2] = -d.z;
+ /* Calculate the free space loss */
+ free_space_loss = free_space_loss_multiplier * t;
+ if (free_space_loss > 1.f) {
+ los->a_te_re[off] = 1.f / free_space_loss;
+ los->a_tm_re[off] = 1.f / free_space_loss;
+ } else
+ los->a_te_re[off] =los-> a_tm_re[off] = 1.f;
+ /* Calculate the delay */
+ los->tau[off] = t / SPEED_OF_LIGHT;
+ /* Calculate the doppler shift */
+ los->freq_shift[off] = vec3_dot(tx_vel_v, &d) - vec3_dot(rx_vel_v, &d);
+ los->freq_shift[off] *= carrier_frequency * 1e9 / SPEED_OF_LIGHT;
+ }
+
+ /***************************************************************************/
+ /* Scatter paths */
+ /***************************************************************************/
+
+ /* Bounce the rays. On each bounce:
+ * - Add path length / SPEED_OF_LIGHT to tau
+ * - Update gains using eqs. (31a)-(31b) from ITU-R P.2040-3
+ * - Spawn scattering rays towards the rx
+ * - TODO Spawn refraction rays with gains as per eqs. (31c)-(31d) from ITU-R P.2040-3
+ */
+ FILE* ray_file = fopen("rays.bin", "wb");
+ if (!ray_file) {
+ fprintf(stderr, "Could not open file rays.bin\n");
+ exit(8);
+ }
+ fwrite(rays, sizeof(Ray), num_tx * num_paths, ray_file);
+
+ FILE* active_file = fopen("active.bin", "wb");
+ if (!active_file) {
+ fprintf(stderr, "Could not open file active.bin\n");
+ exit(8);
+ }
+ fwrite(active, sizeof(uint8_t), num_tx * num_paths / 8 + 1, active_file);
+
+ for (size_t bounce = 0; bounce < num_bounces; ++bounce) {
+ /* active[active_byte_index] & (1 << active_bit_pos) */
+ size_t active_byte_index = 0;
+ uint8_t active_bit = 1;
+ size_t off_tx_path = 0; /* tx * num_paths + path */
+ for (size_t tx = 0; tx < num_tx; ++tx)
+ for (size_t path = 0; path < num_paths; ++path, ++off_tx_path, active_bit <<= 1) {
+
+ /* Check if the ray is active */
+ if (!active_bit) {
+ active_bit = 1;
+ ++active_byte_index;
+ }
+ if (!(active[active_byte_index] & active_bit))
+ continue;
+
+ /***********************************************************************/
+ /* Reflect the rays */
+ /***********************************************************************/
+
+ /* Init */
+ t = -1.f;
+ mesh_ind = face_ind = -1;
+ /* Find the hit point and trinagle and the angle of incidence */
+ moeller_trumbore(&rays[off_tx_path], scene, &t, &mesh_ind, &face_ind, &theta);
+ if (mesh_ind == (uint32_t)-1) { /* Ray hit nothing */
+ active[active_byte_index] &= ~active_bit;
+ --num_active;
+ continue;
+ }
+ /* Calculate the reflection coefficients R_{eTE} and R_{eTM} */
+ material_index = scene->meshes[mesh_ind].material_index;
+ refl_coefs(material_index, theta,
+ &r_te_re, &r_te_im,
+ &r_tm_re, &r_tm_im);
+ /* Calculate the free space loss */
+ free_space_loss = free_space_loss_multiplier * t;
+ free_space_loss *= free_space_loss;
+ if (free_space_loss > 1.f) {
+ r_te_re /= free_space_loss;
+ r_te_im /= free_space_loss;
+ r_tm_re /= free_space_loss;
+ r_tm_im /= free_space_loss;
+ }
+ /* Update the gains as a' = a * refl_coefs */
+ a_te_re_new = a_te_re_refl[off_tx_path] * r_te_re - a_te_im_refl[off_tx_path] * r_te_im;
+ a_te_im_new = a_te_re_refl[off_tx_path] * r_te_im + a_te_im_refl[off_tx_path] * r_te_re;
+ a_tm_re_new = a_tm_re_refl[off_tx_path] * r_tm_re - a_tm_im_refl[off_tx_path] * r_tm_im;
+ a_tm_im_new = a_tm_re_refl[off_tx_path] * r_tm_im + a_tm_im_refl[off_tx_path] * r_tm_re;
+ a_te_re_refl[off_tx_path] = a_te_re_new;
+ a_te_im_refl[off_tx_path] = a_te_im_new;
+ a_tm_re_refl[off_tx_path] = a_tm_re_new;
+ a_tm_im_refl[off_tx_path] = a_tm_im_new;
+ /* Update the delay */
+ tau_t[off_tx_path] += t / SPEED_OF_LIGHT;
+
+ /* Move and reflect the ray */
+ r = rays[off_tx_path];
+ /* Advance the ray to the hit point */
+ *h = vec3_scale(&r.d, t);
+ r.o = vec3_add(h, &r.o);
+ /* Reflect the ray's direction as d' = d - 2*(d \cdot n)*n */
+ n = scene->meshes[mesh_ind].ns[face_ind];
+ t = vec3_dot(&r.d, &n);
+ *h = vec3_scale(&n, 2.f * t);
+ r.d = vec3_sub(&r.d, h);
+ /* Advance the ray origin a bit to avoid intersection with the same triangle */
+ *h = vec3_scale(&r.d, 1e-4f);
+ r.o = vec3_add(&r.o, h);
+
+ /* Calculate the doppler shift caused by the reflection */
+ Vec3* mesh_vel = &scene->meshes[mesh_ind].velocity;
+ Vec3 d_rd = vec3_sub(&r.d, &rays[off_tx_path].d);
+ scatter->freq_shift[off_tx_path] += vec3_dot(&d_rd, mesh_vel) * doppler_shift_multiplier;
+
+ /* Save the reflected ray */
+ rays[off_tx_path] = r;
+
+ /***********************************************************************/
+ /* Scatter the rays */
+ /***********************************************************************/
+
+ /* Scatter a path from the hit point towards the rx */
+ Ray r_scat = { .o = r.o };
+ for (size_t rx = 0; rx < num_rx; ++rx) {
+ /* [rx, tx, bounce, path] */
+ size_t off_scat = ((rx * num_tx + tx) * num_bounces + bounce) * num_paths + path;
+ /* Create a ray from the hit point to the rx */
+ r_scat.d = vec3_sub(&rx_pos_v[rx], &r_scat.o);
+ float d2rx = sqrtf(vec3_dot(&r_scat.d, &r_scat.d));
+ r_scat.d = vec3_normalize(&r_scat.d);
+ /* Check if the ray has a LoS of the rx */
+ t = -1.f;
+ mesh_ind = face_ind = -1;
+ moeller_trumbore(&r_scat, scene, &t, &mesh_ind, &face_ind, &theta);
+ if (mesh_ind != (uint32_t)-1 && t <= 1.f) {
+ /* An obstacle between the hit point and the rx has been hit */
+ scatter->a_te_re[off_scat] = scatter->a_te_im[off_scat]
+ = scatter->a_tm_re[off_scat]
+ = scatter->a_tm_im[off_scat]
+ = scatter->tau[off_scat]
+ = 0.f;
+ continue;
+ }
+ /* No obstacle has been hit */
+ /* Calculate the scattering angle */
+ theta_s = acosf(vec3_dot(&r_scat.d, &n));
+ /* Calculate the scattering coefficients */
+ float a_te_re_new, a_te_im_new, a_tm_re_new, a_tm_im_new;
+ scat_coefs(theta_s, theta, material_index,
+ &a_te_re_new, &a_te_im_new, &a_tm_re_new, &a_tm_im_new);
+ scatter->a_te_re[off_scat] = a_te_re_refl[off_tx_path] * a_te_re_new
+ - a_te_im_refl[off_tx_path] * a_te_im_new;
+ scatter->a_te_im[off_scat] = a_te_re_refl[off_tx_path] * a_te_im_new
+ + a_te_im_refl[off_tx_path] * a_te_re_new;
+ scatter->a_tm_re[off_scat] = a_tm_re_refl[off_tx_path] * a_tm_re_new
+ - a_tm_im_refl[off_tx_path] * a_tm_im_new;
+ scatter->a_tm_im[off_scat] = a_tm_re_refl[off_tx_path] * a_tm_im_new
+ + a_tm_im_refl[off_tx_path] * a_tm_re_new;
+ /* Calculate the direction */
+ scatter->directions_rx[off_scat * 3] = -r_scat.d.x;
+ scatter->directions_rx[off_scat * 3 + 1] = -r_scat.d.y;
+ scatter->directions_rx[off_scat * 3 + 2] = -r_scat.d.z;
+ /* Calculate the delay */
+ scatter->tau[off_scat] = tau_t[off_tx_path] + d2rx / SPEED_OF_LIGHT;
+ /* Calculate the free space loss */
+ free_space_loss = free_space_loss_multiplier * d2rx;
+ free_space_loss *= free_space_loss;
+ if (free_space_loss > 1.f) {
+ scatter->a_te_re[off_scat] /= free_space_loss;
+ scatter->a_te_im[off_scat] /= free_space_loss;
+ scatter->a_tm_re[off_scat] /= free_space_loss;
+ scatter->a_tm_im[off_scat] /= free_space_loss;
+ }
+ /* Calculate the doppler shift */
+ Vec3 d_rd = vec3_sub(&r_scat.d, &rays[off_tx_path].d);
+ float freq_shift = vec3_dot(&d_rd, mesh_vel) * doppler_shift_multiplier;
+ scatter->freq_shift[off_scat] -= freq_shift;
+ }
+
+ /***********************************************************************/
+ /* Refract the rays */
+ /***********************************************************************/
+ /* TODO */
+ }
+ fwrite(rays, sizeof(Ray), num_tx * num_paths, ray_file);
+ fwrite(active, sizeof(uint8_t), num_tx * num_paths / 8 + 1, active_file);
+ }
+ fclose(ray_file);
+ fclose(active_file);
+
+ /* Free */
+ free(a_te_re_refl);
+ free(a_te_im_refl);
+ free(a_tm_re_refl);
+ free(a_tm_im_refl);
+ free(tau_t);
+ free(h);
+ free(active);
+ free(rays);
+ /* Free the scene */
+ /* TODO */
+}
diff --git a/src/materials.c b/src/materials.c
@@ -0,0 +1,122 @@
+#include "../inc/materials.h"
+
+Material g_materials[NUM_G_MATERIALS] = {
+ /* 0 */
+ {3, "air",
+ 1.f, 0.f, 0.f, 0.001f,
+ 0.1f, 0.5f, 0.3f, 0.2f,
+ 2, 2},
+ /* 1 */
+ {8, "concrete",
+ 5.24f, 0.f, 0.0462f, 0.7822f,
+ 0.5f, 0.33f, 0.34f, 0.33f,
+ 4, 4},
+ /* 2 */
+ {5, "brick",
+ 3.91f, 0.f, 0.0238f, 0.16f,
+ 0.4f, 0.4f, 0.3f, 0.3f,
+ 3, 3},
+ /* 3 */
+ {12, "plasterboard",
+ 2.73f, 0.f, 0.0085f, 0.9395f,
+ 0.3f, 0.4f, 0.4f, 0.2f,
+ 3, 3},
+ /* 4 */
+ {4, "wood",
+ 1.99f, 0.f, 0.0047f, 1.0718f,
+ 0.2f, 0.5f, 0.3f, 0.2f,
+ 2, 2},
+ /* 5 */
+ {5, "glass",
+ 6.31f, 0.f, 0.0036f, 1.3394f,
+ 0.3f, 0.4f, 0.4f, 0.2f,
+ 3, 3},
+ /* 6 */
+ {5, "glass",
+ 5.79f, 0.f, 0.0004f, 1.658f,
+ 0.3f, 0.4f, 0.4f, 0.2f,
+ 3, 3},
+ /* 7 */
+ {13, "ceiling board",
+ 1.48f, 0.f, 0.0011f, 1.0750f,
+ 0.2f, 0.5f, 0.3f, 0.2f,
+ 2, 2},
+ /* 8 */
+ {13, "ceiling board",
+ 1.52f, 0.f, 0.0029f, 1.029f,
+ 0.2f, 0.5f, 0.3f, 0.2f,
+ 2, 2},
+ /* 9 */
+ {9, "chipboard",
+ 2.58f, 0.f, 0.0217f, 0.7800f,
+ 0.4f, 0.4f, 0.3f, 0.3f,
+ 3, 3},
+ /* 10 */
+ {7, "plywood",
+ 2.71f, 0.f, 0.33f, 0.f,
+ 0.3f, 0.5f, 0.3f, 0.2f,
+ 3, 3},
+ /* 11 */
+ {6, "marble",
+ 7.074f, 0.f, 0.0055f,
+ 0.9262f, 0.3f, 0.4f, 0.4f, 0.2f,
+ 3, 3},
+ /* 12 */
+ {10, "floorboard",
+ 3.66f, 0.f, 0.0044f, 1.3515f,
+ 0.3f, 0.4f, 0.4f, 0.2f,
+ 3, 3},
+ /* 13 */
+ {5, "metal",
+ 1.f, 0.f, 10000000.f, 0.f,
+ 0.f, 0.f, 1.f, 0.f,
+ 1, 1},
+ /* 14 */
+ {15, "very dry ground",
+ 3.f, 0.f, 0.00015f, 2.52f,
+ 0.4f, 0.3f, 0.4f, 0.3f,
+ 4, 4},
+ /* 15 */
+ {17, "medium dry ground",
+ 15.f, -0.1f, 0.035f, 1.63f,
+ 0.5f, 0.33f, 0.34f, 0.33f,
+ 4, 4},
+ /* 16 */
+ {10, "wet ground",
+ 30.f, -0.4f, 0.15f, 1.30f,
+ 0.5f, 0.33f, 0.34f, 0.33f,
+ 4, 4}
+};
+
+
+/* Map material name strings to material indices */
+typedef struct {
+ /* Material name. Null-terminated string. Lower case. */
+ const char *name;
+ MaterialIndex index;
+} MaterialMap;
+MaterialMap g_str2material_map[] = {
+ {"air", MATERIAL_AIR},
+ {"concrete", MATERIAL_CONCRETE},
+ {"brick", MATERIAL_BRICK},
+ {"plasterboard", MATERIAL_PLASTERBOARD},
+ {"wood", MATERIAL_WOOD},
+ {"glass1", MATERIAL_GLASS1},
+ {"glass2", MATERIAL_GLASS2},
+ {"ceiling_board1", MATERIAL_CEILING_BOARD1},
+ {"ceiling_board2", MATERIAL_CEILING_BOARD2},
+ {"chipboard", MATERIAL_CHIPBOARD},
+ {"plywood", MATERIAL_PLYWOOD},
+ {"marble", MATERIAL_MARBLE},
+ {"floorboard", MATERIAL_FLOORBOARD},
+ {"metal", MATERIAL_METAL},
+ {"very_dry_ground", MATERIAL_VERY_DRY_GROUND},
+ {"medium_dry_ground", MATERIAL_MEDIUM_DRY_GROUND},
+ {"wet_ground", MATERIAL_WET_GROUND}
+};
+MaterialIndex get_material_index(const char *name) {
+ for (uint32_t i = 0; i < NUM_G_MATERIALS; ++i)
+ if (strcmp(g_str2material_map[i].name, name) == 0)
+ return g_str2material_map[i].index;
+ return MATERIAL_AIR; /* Invalid material. Default to air */
+}
diff --git a/src/scene.c b/src/scene.c
@@ -0,0 +1,84 @@
+#include "../inc/scene.h" /* for Scene, Mesh, Material */
+#include "../inc/materials.h" /* for g_hrt_materials */
+#include "../inc/common.h" /* for PERROT_CLEANUP_EXIT */
+
+#include <stdio.h> /* for fopen, FILE, fclose */
+
+void scene_save(
+ IN Scene* scene,
+ IN const char* filepath
+)
+{
+ FILE *fp = fopen(filepath, "wb");
+ if (!fp)
+ PERROR_CLEANUP_EXIT("Error: cannot open file", 8);
+
+ /* MAGIC */
+ fwrite("HRT", 1, 3, fp);
+ /* SCENE */
+ /* num_meshes */
+ fwrite(&scene->num_meshes, sizeof(uint32_t), 1, fp);
+ /* meshes */
+ for (uint32_t i = 0; i != scene->num_meshes; ++i) {
+ /* MESH */
+ Mesh *mesh = &scene->meshes[i];
+ fwrite(&mesh->num_vertices, sizeof(uint32_t), 1, fp);
+ fwrite(mesh->vs, sizeof(Vec3), mesh->num_vertices, fp);
+ fwrite(&mesh->num_triangles, sizeof(uint32_t), 1, fp);
+ fwrite(mesh->is, sizeof(uint32_t), 3 * mesh->num_triangles, fp);
+ fwrite(&mesh->material_index, sizeof(uint32_t), 1, fp);
+ fwrite(&mesh->velocity, sizeof(Vec3), 1, fp);
+ }
+
+ fclose(fp);
+}
+
+Scene scene_load(IN const char *filepath)
+{
+ /* Open the HRT file */
+ FILE *f = fopen(filepath, "rb");
+ if (!f)
+ PERROR_CLEANUP_EXIT("Could not open scene file", 8);
+
+ /* Parse the file */
+ /* MAGIC */
+ char magic[3];
+ if (fread(magic, 1, 3, f) != 3) exit(8);
+ if (strncmp(magic, "HRT", 3)) exit(8);
+ /* SCENE */
+ Scene scene;
+ /* num_meshes */
+ if (fread(&scene.num_meshes, sizeof(uint32_t), 1, f) != 1) exit(8);
+ if (scene.num_meshes == 0)
+ PERROR_CLEANUP_EXIT("Scene has no meshes", 8);
+ if (scene.num_meshes > 1000)
+ PERROR_CLEANUP_EXIT("Scene has too many meshes", 8);
+ /* meshes */
+ scene.meshes = (Mesh*)malloc(scene.num_meshes * sizeof(Mesh));
+ for (uint32_t i = 0; i != scene.num_meshes; ++i) {
+ /* MESH */
+ Mesh *mesh = &scene.meshes[i];
+ /* num_vertices */
+ if(fread(&mesh->num_vertices, sizeof(uint32_t), 1, f) != 1) exit(8);
+ /* vertices */
+ mesh->vs = (Vec3*)malloc(mesh->num_vertices * sizeof(Vec3));
+ if (fread(mesh->vs, sizeof(float) * 3, mesh->num_vertices, f) != mesh->num_vertices)
+ exit(8);
+ /* num_triangles */
+ if (fread(&mesh->num_triangles, sizeof(uint32_t), 1, f) != 1) exit(8);
+ /* triangles */
+ mesh->is = (uint32_t*)malloc(mesh->num_triangles * 3 * sizeof(uint32_t));
+ if (fread(mesh->is, sizeof(uint32_t), mesh->num_triangles * 3, f) != mesh->num_triangles * 3)
+ exit(8);
+ /* material_index */
+ if (fread(&mesh->material_index, sizeof(uint32_t), 1, f) != 1) exit(8);
+ /* velocity */
+ if (fread(&mesh->velocity, sizeof(Vec3), 1, f) != 1) exit(8);
+ /* normals */
+ mesh->ns = NULL;
+ }
+
+ fclose(f);
+ return scene;
+}
+
diff --git a/src/scene_fromSionna.c b/src/scene_fromSionna.c
@@ -0,0 +1,482 @@
+/* vim: set tabstop=2:softtabstop=2:shiftwidth=2:noexpandtab */
+
+#include "../inc/vec3.h" /* for Vec3 */
+#include "../inc/common.h" /* for IN, OUT, PERROR_CLEANUP_EXIT */
+#include "../inc/scene.h" /* for Scene, Mesh, Material, scene_save */
+#include "../inc/materials.h" /* for g_hrt_materials, MATERIAL_CONCRETE */
+
+#include <stdio.h> /* for FILE, fopen, fclose, fseek, fgets, sscanf, printf */
+#include <stdlib.h> /* for malloc */
+#include <stdint.h> /* for uint8_t, uint32_t */
+#include <string.h> /* for strncmp, strrchr, strcmp */
+
+/**
+ * Hardcoded scenes
+ */
+
+/* Box */
+float mesh_box_vs[] = {
+ 5.f, 5.f, 0.f,
+ -5.f, 5.f, 0.f,
+ -5.f, -5.f, 0.f,
+ 5.f, -5.f, 0.f,
+ 5.f, 5.f, 5.f,
+ -5.f, 5.f, 5.f,
+ -5.f, -5.f, 5.f,
+ 5.f, -5.f, 5.f,
+};
+uint32_t mesh_box_is[] = {
+ 0, 1, 2,
+ 0, 2, 3,
+ 0, 4, 5,
+ 0, 5, 1,
+ 1, 5, 6,
+ 1, 6, 2,
+ 2, 6, 7,
+ 2, 7, 3,
+ 3, 7, 4,
+ 3, 4, 0,
+ 4, 7, 6,
+ 4, 6, 5,
+};
+Mesh mesh_box = {
+ .num_vertices = 8,
+ .vs = (Vec3*)mesh_box_vs,
+ .num_triangles = 12,
+ .is = (uint32_t*)mesh_box_is,
+ .material_index = MATERIAL_CONCRETE,
+ .velocity = {0.f, 0.f, 0.f},
+};
+Scene scene_box = {
+ .num_meshes = 1,
+ .meshes = &mesh_box,
+};
+const char* scene_box_filename = "box.xml";
+
+/* Simple reflector */
+float mesh_simpleReflector_vs[] = {
+ -0.5f, -0.5f, 0.f,
+ 0.5f, -0.5f, 0.f,
+ 0.5f, 0.5f, 0.f,
+ -0.5f, 0.5f, 0.f,
+};
+uint32_t mesh_simpleReflector_is[] = {
+ 0, 1, 2,
+ 0, 2, 3,
+};
+Mesh mesh_simpleReflector = {
+ .num_vertices = 4,
+ .vs = (Vec3*)mesh_simpleReflector_vs,
+ .num_triangles = 2,
+ .is = (uint32_t*)mesh_simpleReflector_is,
+ .material_index = MATERIAL_CONCRETE,
+ .velocity = {0.f, 0.f, 0.f},
+};
+Scene scene_simpleReflector = {
+ .num_meshes = 1,
+ .meshes = &mesh_simpleReflector,
+};
+const char* scene_simpleReflector_filename = "simple_reflector.xml";
+
+/**
+ * Read a Sionna scene
+ */
+
+/* Read a binary PLY file
+ *
+ * The PLY file have the following header:
+ * ply
+ * format binary_little_endian 1.0
+ * element vertex <num_vertices>
+ * property float x
+ * property float y
+ * property float z
+ * property float s
+ * property float t
+ * element face <num_triangles>
+ * property list uchar int vertex_index
+ * end_header
+ *
+ * \param filepath Path to the PLY file
+ * \return The mesh
+ */
+Mesh readPly(const char* filepath) {
+ FILE* f = fopen(filepath, "rb");
+ if (f == NULL)
+ PERROR_CLEANUP_EXIT("Error: cannot open file", 8);
+
+ char buf[256];
+ Mesh mesh = {
+ .num_vertices = 0,
+ .vs = NULL,
+ .num_triangles = 0,
+ .is = NULL,
+ .material_index = 0,
+ .velocity = {0.f, 0.f, 0.f},
+ };
+
+ /* read header */
+ while(fgets(buf, 256, f)) {
+ if (!strncmp(buf, "end_header", 10))
+ break;
+ if (!strncmp(buf, "element vertex ", 15)) {
+ if (sscanf(buf, "element vertex %u", &mesh.num_vertices) != 1)
+ PERROR_CLEANUP_EXIT("Error: cannot read number of vertices", 8);
+ }
+ else if (!strncmp(buf, "element face ", 13))
+ if (sscanf(buf, "element face %u", &mesh.num_triangles) != 1)
+ PERROR_CLEANUP_EXIT("Error: cannot read number of triangles", 8);
+ }
+
+ /* check if vertex and face elements are found */
+ if (mesh.num_vertices == 0 || mesh.num_triangles == 0)
+ PERROR_CLEANUP_EXIT("Error: PLY element vertex or element face not found", 8);
+ /* check if the numbers are not too big */
+ else if (mesh.num_vertices > 1000000 || mesh.num_triangles > 1000000)
+ PERROR_CLEANUP_EXIT("Error: PLY element vertex or element face too big", 8);
+
+ /* allocate memory for vertices and faces */
+ mesh.vs = (Vec3*)malloc(mesh.num_vertices * sizeof(Vec3));
+ mesh.is = (uint32_t*)malloc(3 * mesh.num_triangles * sizeof(uint32_t));
+
+ /* read vertices */
+ for (uint32_t i = 0; i != mesh.num_vertices; ++i) {
+ if (fread(&mesh.vs[i], sizeof(Vec3), 1, f) != 1)
+ PERROR_CLEANUP_EXIT("Error: cannot read vertex", 8, mesh.vs, mesh.is);
+ /* skip s and t */
+ if (fseek(f, 8, SEEK_CUR) != 0)
+ PERROR_CLEANUP_EXIT("Error: cannot skip s and t", 8, mesh.vs, mesh.is);
+ }
+
+ /* read faces */
+ for (uint32_t i = 0; i != mesh.num_triangles * 3; i += 3) {
+ uint8_t n;
+ if (fread(&n, sizeof(uint8_t), 1, f) != 1)
+ PERROR_CLEANUP_EXIT("Error: cannot read number of vertices in face", 8, mesh.vs, mesh.is);
+ if (n != 3)
+ PERROR_CLEANUP_EXIT("Error: face is not a triangle", 8, mesh.vs, mesh.is);
+ if (fread(&mesh.is[i], sizeof(uint32_t), 3, f) != 3)
+ PERROR_CLEANUP_EXIT("Error: cannot read face", 8, mesh.vs, mesh.is);
+ }
+
+ fclose(f);
+ return mesh;
+}
+
+/** Read a scene config CSV file.
+ *
+ * The CSV file must be a text file with the following format:
+ * name,material_index,velocity_x,velocity_y,velocity_z
+ *
+ * The file can have any number of lines except the first one, which is the header.
+ * If a mesh is not in the CSV file,
+ * it will have the default (0) material index and velocity.
+ *
+ * \param filepath Path to the CSV file
+ * \param num_meshes Output number of meshes in the CSV file
+ * \param mesh_filenames Output array of mesh filenames
+ * \param material_indicies Output array of material indicies
+ * \param velocity Output array of velocity vectors. Size [num_meshes]
+ */
+void readCsv(
+ IN const char* filepath,
+ OUT uint32_t* csv_num_meshes,
+ OUT char*** csv_mesh_filenames,
+ OUT uint32_t** csv_material_indicies,
+ OUT Vec3** csv_velocities
+)
+{
+ FILE* f = fopen(filepath, "r");
+ if (!f)
+ PERROR_CLEANUP_EXIT("Error: cannot open file", 8);
+
+ /* Check header */
+ char buf[256];
+ if (!fgets(buf, 256, f))
+ PERROR_CLEANUP_EXIT("Error: cannot read header", 8);
+ if (strncmp(buf, "name,material_index,velocity_x,velocity_y,velocity_z\n", 48))
+ PERROR_CLEANUP_EXIT("Error: invalid header", 8);
+
+ /* Count number of lines */
+ *csv_num_meshes = 0;
+ while (fgets(buf, 256, f))
+ ++(*csv_num_meshes);
+ if (*csv_num_meshes == 0) {
+ fclose(f);
+ return;
+ }
+ /* Return to the beginning of the file */
+ rewind(f);
+ /* Skip header */
+ if (!fgets(buf, 256, f))
+ PERROR_CLEANUP_EXIT("Error: cannot read header", 8);
+
+ /* Allocate memory for the arrays */
+ *csv_mesh_filenames = (char**)malloc(*csv_num_meshes * sizeof(char*));
+ *csv_material_indicies = (uint32_t*)malloc(*csv_num_meshes * sizeof(uint32_t));
+ *csv_velocities = (Vec3*)malloc(*csv_num_meshes * sizeof(Vec3));
+
+ /* Read lines */
+ for (uint32_t i = 0; i != *csv_num_meshes; ++i) {
+ /* Read line */
+ if (!fgets(buf, 256, f))
+ PERROR_CLEANUP_EXIT("Error: cannot read line", 8);
+ /* Parse line */
+ char name[50];
+ int rc = sscanf(
+ buf,
+ "%49[^,],%u,%f,%f,%f\n",
+ name,
+ &(*csv_material_indicies)[i],
+ &(*csv_velocities)[i].x,
+ &(*csv_velocities)[i].y,
+ &(*csv_velocities)[i].z
+ );
+ if (rc != 4)
+ PERROR_CLEANUP_EXIT("Error: cannot parse line", 8);
+ /* Copy name */
+ (*csv_mesh_filenames)[i] = (char*)malloc(strlen(name) + 1);
+ if (!csv_mesh_filenames[i])
+ PERROR_CLEANUP_EXIT("Error: cannot allocate memory for mesh filename", 8);
+ strcpy((*csv_mesh_filenames)[i], name);
+ }
+}
+
+/** Read a Sionna scene XML file.
+ *
+ * Extracts each shape's name, file path and material.
+ *
+ * \param filepath Path to the XML file
+ * \param xml_num_meshes Output number of meshes in the XML file
+ * \param xml_mesh_filepaths Output array of mesh file paths (null-terminated strings)
+ * \param xml_mesh_names Output array of mesh names (null-terminated strings)
+ * \param xml_mesh_materials Output array of mesh materials (null-terminated strings)
+ */
+void readXml(
+ IN const char* filepath,
+ OUT uint32_t* xml_num_meshes,
+ OUT char*** xml_mesh_filepaths,
+ OUT char*** xml_mesh_names,
+ OUT char*** xml_mesh_materials
+)
+{
+ FILE* f = fopen(filepath, "r");
+ if (f == NULL)
+ PERROR_CLEANUP_EXIT("Error: cannot open the xml file", 8);
+
+ fseek(f, 0, SEEK_END);
+ long fsize = ftell(f);
+ fseek(f, 0, SEEK_SET);
+
+ char* buf = (char*)malloc(fsize + 1);
+ if (!buf)
+ PERROR_CLEANUP_EXIT("Error: cannot allocate memory for the xml file", 8);
+ if (fread(buf, 1, fsize, f) != fsize)
+ PERROR_CLEANUP_EXIT("Error: cannot read the xml file", 8, buf);
+ buf[fsize] = '\0';
+
+ /* Count the amount of "<shape" entries */
+ char** shape_positions = NULL;
+ *xml_num_meshes = 0;
+ char* pos = buf, *end;
+ while ((pos = strstr(pos, "<shape")) != NULL) {
+ ++*xml_num_meshes;
+ shape_positions = (char**)realloc(shape_positions, *xml_num_meshes * sizeof(char*));
+ if (!shape_positions)
+ PERROR_CLEANUP_EXIT("Error: cannot reallocate memory for shape positions", 8, buf);
+ shape_positions[*xml_num_meshes - 1] = pos;
+ pos += 6;
+ }
+ if (*xml_num_meshes == 0)
+ PERROR_CLEANUP_EXIT("Error: no shapes found in the xml file", 8, buf);
+ if (!shape_positions)
+ PERROR_CLEANUP_EXIT("Error: cannot allocate memory for shape positions", 8, buf);
+
+ /* Allocate memory for the mesh data */
+ *xml_mesh_filepaths = (char**)malloc(*xml_num_meshes * sizeof(char*));
+ *xml_mesh_names = (char**)malloc(*xml_num_meshes * sizeof(char*));
+ *xml_mesh_materials = (char**)malloc(*xml_num_meshes * sizeof(char*));
+ if (!*xml_mesh_filepaths || !*xml_mesh_names || !*xml_mesh_materials)
+ PERROR_CLEANUP_EXIT("Error: cannot allocate memory for mesh data", 8, buf);
+
+ /* For each "<shape" entry, extract the name, file path and material */
+ char *e_name, *e_filepath, *e_material;
+ for (uint32_t i = 0; i != *xml_num_meshes; ++i) {
+ pos = shape_positions[i] + 6;
+
+ /* Find mesh name */
+ if ((pos = strstr(pos, "name=\"")) == NULL)
+ PERROR_CLEANUP_EXIT("Error: cannot find mesh name", 8,
+ buf, shape_positions, *xml_mesh_filepaths, *xml_mesh_names, *xml_mesh_materials);
+ pos += 6;
+ end = strchr(pos, '\"');
+ if (!end)
+ PERROR_CLEANUP_EXIT("Error: cannot find mesh name end", 8,
+ buf, shape_positions, *xml_mesh_filepaths, *xml_mesh_names, *xml_mesh_materials);
+ e_name = (char*)malloc(end - pos + 1);
+ if (!e_name)
+ PERROR_CLEANUP_EXIT("Error: cannot allocate memory for mesh name", 8,
+ buf, shape_positions, *xml_mesh_filepaths, *xml_mesh_names, *xml_mesh_materials);
+ strncpy(e_name, pos, end - pos);
+ e_name[end - pos] = '\0';
+
+ /* Find mesh file path */
+ if (!(pos = strstr(pos, "<string name=\"filename\"")))
+ PERROR_CLEANUP_EXIT("Error: cannot find mesh file path", 8, buf, e_name,
+ shape_positions, *xml_mesh_filepaths, *xml_mesh_names, *xml_mesh_materials);
+ if (!(pos = strstr(pos, "value=\"")))
+ PERROR_CLEANUP_EXIT("Error: cannot find mesh file path value", 8, buf, e_name,
+ shape_positions, *xml_mesh_filepaths, *xml_mesh_names, *xml_mesh_materials);
+ pos += 7;
+ end = strchr(pos, '\"');
+ if (!end)
+ PERROR_CLEANUP_EXIT("Error: cannot find mesh file path end", 8, buf, e_name,
+ shape_positions, *xml_mesh_filepaths, *xml_mesh_names, *xml_mesh_materials);
+ e_filepath = (char*)malloc(end - pos + 1);
+ if (!e_filepath)
+ PERROR_CLEANUP_EXIT("Error: cannot allocate memory for mesh file path", 8, buf, e_name,
+ shape_positions, *xml_mesh_filepaths, *xml_mesh_names, *xml_mesh_materials);
+ strncpy(e_filepath, pos, end - pos);
+ e_filepath[end - pos] = '\0';
+
+ /* Find mesh material */
+ if (!(pos = strstr(pos, "id=\"mat-itu_")))
+ PERROR_CLEANUP_EXIT("Error: cannot find mesh material", 8, buf, e_name, e_filepath,
+ shape_positions, *xml_mesh_filepaths, *xml_mesh_names, *xml_mesh_materials);
+ pos += 12;
+ end = strchr(pos, '\"');
+ if (!end)
+ PERROR_CLEANUP_EXIT("Error: cannot find mesh material end", 8, buf, e_name, e_filepath,
+ shape_positions, *xml_mesh_filepaths, *xml_mesh_names, *xml_mesh_materials);
+ e_material = (char*)malloc(end - pos + 1);
+ if (!e_material)
+ PERROR_CLEANUP_EXIT("Error: cannot allocate memory for mesh material", 8, buf, e_name, e_filepath,
+ shape_positions, *xml_mesh_filepaths, *xml_mesh_names, *xml_mesh_materials);
+ strncpy(e_material, pos, end - pos);
+ e_material[end - pos] = '\0';
+
+ /* Save mesh data */
+ (*xml_mesh_filepaths)[i] = e_filepath;
+ (*xml_mesh_names)[i] = e_name;
+ (*xml_mesh_materials)[i] = e_material;
+ }
+
+ free(shape_positions);
+ free(buf);
+}
+
+/* Read a Sionna .xml scene. Assumes meshes are in a "meshes" directory.
+ *
+ * \param filepath Path to the .xml scene file
+ * \param scene The scene to fill
+ */
+void readScene(const char* filepath, Scene* scene) {
+ /* Check if filepath ends with ".xml" */
+ size_t filepath_len = strlen(filepath);
+ if (filepath_len > 4 && strcmp(filepath + filepath_len - 4, ".xml") != 0)
+ PERROR_CLEANUP_EXIT("Error: scene file must end with .xml", 8);
+
+ /* Read the xml file */
+ uint32_t xml_num_meshes;
+ char** xml_mesh_filepaths;
+ char** xml_mesh_names;
+ char** xml_mesh_materials;
+ readXml(
+ filepath,
+ &xml_num_meshes,
+ &xml_mesh_filepaths,
+ &xml_mesh_names,
+ &xml_mesh_materials
+ );
+
+ /* Read the scene config CSV file */
+ char* csv_path = (char*)malloc(filepath_len + 1);
+ if (!csv_path)
+ PERROR_CLEANUP_EXIT("Error: cannot allocate memory for csv_path", 8);
+ strcpy(csv_path, filepath);
+ strcpy(csv_path + filepath_len - 4, ".csv");
+ uint32_t csv_num_meshes;
+ char** csv_mesh_filenames;
+ uint32_t* csv_material_indicies;
+ Vec3* csv_velocities;
+ readCsv(
+ csv_path,
+ &csv_num_meshes,
+ &csv_mesh_filenames,
+ &csv_material_indicies,
+ &csv_velocities
+ );
+ free(csv_path);
+
+ /* Read the meshes */
+
+ /* Get the scene directory path */
+ size_t scene_dir_len = strrchr(filepath, '/') - filepath + 1;
+ /* if filepath is "/path/to/scene.xml", scene_dir is "/path/to/" */
+ char* scene_dir = (char*)malloc(scene_dir_len + 7);
+ if (!scene_dir)
+ PERROR_CLEANUP_EXIT("Error: failed to allocate scene_dir", 8);
+ strncpy(scene_dir, filepath, scene_dir_len);
+
+ /* For each ply file from the xml file */
+ scene->num_meshes = xml_num_meshes;
+ scene->meshes = (Mesh*)malloc(xml_num_meshes * sizeof(Mesh));
+ if (!scene->meshes && xml_num_meshes > 0)
+ PERROR_CLEANUP_EXIT("Error: failed to allocate meshes array", 8);
+ for (uint32_t i = 0; i != xml_num_meshes; ++i) {
+ /* Get the mesh file path */
+ size_t mesh_filepath_len = scene_dir_len + strlen(xml_mesh_filepaths[i]);
+ char* mesh_filepath = (char*)malloc(mesh_filepath_len);
+ mesh_filepath[0] = '\0';
+ strncat(mesh_filepath, scene_dir, scene_dir_len);
+ strncat(mesh_filepath + scene_dir_len, xml_mesh_filepaths[i], mesh_filepath_len - scene_dir_len);
+ /* Read the mesh */
+ scene->meshes[i] = readPly(mesh_filepath);
+ free(mesh_filepath);
+ /* Set the material index from the material name from xml */
+ scene->meshes[i].material_index = get_material_index(xml_mesh_materials[i]);
+ /* Set the material index and velocity from CSV if present */
+ for (uint32_t j = 0; j != csv_num_meshes; ++j)
+ if (strcmp(xml_mesh_names[i], csv_mesh_filenames[j]) == 0) {
+ scene->meshes[i].material_index = csv_material_indicies[j];
+ scene->meshes[i].velocity = csv_velocities[j];
+ break;
+ }
+ }
+
+ free(scene_dir);
+}
+
+
+/**
+ * Main
+ */
+
+int main(int argc, char *argv[]) {
+ if (argc != 2) {
+ printf("Usage: %s <scene.xml>\n", argv[0]);
+ return 1;
+ }
+ const char* scene_filepath = argv[1];
+ const char* scene_filename = strrchr(scene_filepath, '/');
+ if (scene_filename == NULL) exit(8);
+ else ++scene_filename;
+
+ Scene* scene_ptr;
+
+ /* read scene */
+ /* check if the scene is a hardcoded scene */
+ if (!strcmp(scene_filename, scene_box_filename))
+ scene_ptr = &scene_box;
+ else if (!strcmp(scene_filename, scene_simpleReflector_filename))
+ scene_ptr = &scene_simpleReflector;
+ else {
+ scene_ptr = (Scene*)malloc(sizeof(Scene));
+ readScene(scene_filepath, scene_ptr);
+ }
+
+ /* Save the scene as HRT */
+ scene_save(scene_ptr, "scene.hrt");
+
+ return 0;
+}
diff --git a/src/vizrays.c b/src/vizrays.c
@@ -0,0 +1,336 @@
+#include "../inc/vec3.h" /* for Vec3 */
+#include "../inc/scene.h" /* for Scene, Mesh, scene_load */
+
+#include "../test/test.h" /* for g_numPaths, g_numBounces, g_numTx, g_numRx */
+
+#include <GL/glut.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+#include <stdint.h>
+#include <ctype.h> /* for tolower */
+
+/**
+ * Structs
+ */
+
+typedef struct {
+ Vec3 o, t;
+ float d;
+ float lastX, lastY;
+ float yaw, pitch, roll;
+} Camera;
+
+/**
+ * Globals
+ */
+
+Camera g_cam = {
+ {0.f, 0.f, 5.f},
+ {0.f, 0.f, 0.f},
+ 5.f,
+ 0.f, 0.f,
+ 0.f, 0.f, 0.f
+};
+uint8_t g_mouseDown = 0;
+
+float* g_rays = NULL;
+uint8_t *g_active = NULL;
+
+Scene g_scene = {0};
+
+uint32_t g_bounce_cur = 0;
+
+/**
+ * Loading functions
+ */
+
+void loadRays(const char* filename) {
+ FILE* file = fopen(filename, "rb");
+ if (!file) {
+ printf("Failed to open file\n");
+ exit(8);
+ }
+
+ uint32_t fileSize = (g_numBounces + 1) * g_numTx * g_numPaths * 2 * 3;
+ g_rays = (float*)malloc(fileSize * sizeof(float));
+ if (fread(g_rays, sizeof(float), fileSize, file) != fileSize) {
+ printf("Failed to read rays\n");
+ exit(8);
+ }
+ fclose(file);
+
+ /* active.bin */
+ file = fopen("active.bin", "rb");
+ if (!file) {
+ printf("Failed to open file\n");
+ exit(8);
+ }
+ uint32_t activeSize = (g_numBounces + 1) * (g_numTx * g_numPaths / 8 + 1);
+ g_active = (uint8_t*)malloc(activeSize * sizeof(uint8_t));
+ if (fread(g_active, sizeof(uint8_t), activeSize, file) != activeSize) {
+ printf("Failed to read active\n");
+ exit(8);
+ }
+ fclose(file);
+}
+
+/**
+ * Draw functions
+ */
+
+void drawScene() {
+ glBegin(GL_TRIANGLES);
+ for (uint32_t i = 0; i < g_scene.num_meshes; i++) {
+ Mesh* mesh = &g_scene.meshes[i];
+ float mesh_color = (float)i / g_scene.num_meshes;
+ glColor3f(mesh_color, 1.f - mesh_color, 0.f);
+
+ for (uint32_t j = 0; j < mesh->num_triangles; j++) {
+ uint32_t idx = mesh->is[j * 3];
+ Vec3 v1 = mesh->vs[idx];
+ idx = mesh->is[j * 3 + 1];
+ Vec3 v2 = mesh->vs[idx];
+ idx = mesh->is[j * 3 + 2];
+ Vec3 v3 = mesh->vs[idx];
+ glVertex3f(v1.x, v1.y, v1.z);
+ glVertex3f(v2.x, v2.y, v2.z);
+ glVertex3f(v3.x, v3.y, v3.z);
+ }
+ }
+ glEnd();
+}
+
+void drawRays() {
+ float color = (float)g_bounce_cur / g_numBounces;
+ glColor3f(color, 0.0f, 1.0f - color);
+
+ glBegin(GL_LINES);
+ uint32_t idx_base = g_bounce_cur * g_numTx * g_numPaths * 2 * 3;
+ uint32_t active_byte = g_bounce_cur * (g_numTx * g_numPaths / 8 + 1);
+ uint8_t active_bit = 1;
+ uint32_t num_skipped = 0;
+ for (uint32_t path = 0; path < g_numPaths; ++path, active_bit <<= 1) {
+ uint32_t idx = idx_base + path * 2 * 3;
+
+ if (!active_bit) {
+ active_byte++;
+ active_bit = 1;
+ }
+ if (!(g_active[active_byte] & active_bit))
+ {
+ num_skipped++;
+ continue;
+ }
+
+ /* Origin */
+ float x1 = g_rays[idx];
+ float y1 = g_rays[idx + 1];
+ float z1 = g_rays[idx + 2];
+
+ /* Destination */
+ float x2 = x1 + g_rays[idx + 3];
+ float y2 = y1 + g_rays[idx + 4];
+ float z2 = z1 + g_rays[idx + 5];
+
+ glVertex3f(x1, y1, z1);
+ glVertex3f(x2, y2, z2);
+ }
+ glEnd();
+ printf("\rSkipped %d paths", num_skipped);
+ fflush(stdout);
+
+ /* Draw points */
+ glPointSize(5.0f);
+ glBegin(GL_POINTS);
+ active_byte = g_bounce_cur * (g_numTx * g_numPaths / 8 + 1);
+ active_bit = 1;
+ for (uint32_t path = 0; path < g_numPaths; ++path) {
+ uint32_t idx = idx_base + path * 2 * 3;
+ if (!active_bit) {
+ active_byte++;
+ active_bit = 1;
+ }
+ if (!(g_active[active_byte] & active_bit))
+ continue;
+ glVertex3f(g_rays[idx], g_rays[idx + 1], g_rays[idx + 2]);
+ }
+ glEnd();
+}
+
+/**
+ * Callbacks
+ */
+
+void reshape(int w, int h) {
+ glViewport(0, 0, w, h);
+ glMatrixMode(GL_PROJECTION);
+ glLoadIdentity();
+ gluPerspective(60.0f, (float)w/h, 0.1f, 1000.0f);
+ glMatrixMode(GL_MODELVIEW);
+}
+
+void display() {
+ glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
+ glLoadIdentity();
+
+ /* Look at target from camera position */
+ gluLookAt(
+ g_cam.o.x, g_cam.o.y, g_cam.o.z,
+ g_cam.t.x, g_cam.t.y, g_cam.t.z,
+ 0.0f, 1.0f, 0.0f
+ );
+
+ drawScene();
+ drawRays();
+
+ glutSwapBuffers();
+}
+
+void updateCamera() {
+ float x = g_cam.d * cosf(g_cam.pitch) * cosf(g_cam.yaw);
+ float y = g_cam.d * sinf(g_cam.pitch);
+ float z = g_cam.d * cosf(g_cam.pitch) * sinf(g_cam.yaw);
+
+ g_cam.o.x = g_cam.t.x + x;
+ g_cam.o.y = g_cam.t.y + y;
+ g_cam.o.z = g_cam.t.z + z;
+
+ gluLookAt(g_cam.o.x, g_cam.o.y, g_cam.o.z,
+ g_cam.t.x, g_cam.t.y, g_cam.t.z,
+ 0.f, 1.f, 0.f);
+}
+
+void keyboard(unsigned char key, int x, int y) {
+ float moveSpeed = 0.1f;
+ float tiltSpeed = 0.15f;
+ /* Speed up if SHIFT is down */
+ if (glutGetModifiers() & GLUT_ACTIVE_SHIFT) {
+ moveSpeed *= 5;
+ key = tolower(key);
+ }
+
+ /* Calculate forward direction */
+ Vec3 dir = {
+ g_cam.t.x - g_cam.o.x,
+ g_cam.t.y - g_cam.o.y,
+ g_cam.t.z - g_cam.o.z
+ };
+ float norm = sqrtf(dir.x * dir.x + dir.y * dir.y + dir.z * dir.z);
+ Vec3 forward = {dir.x / norm, dir.y / norm, dir.z / norm};
+
+ /* Calculate right vector (cross product of forward and up (0,0,1)) */
+ norm = sqrtf(forward.z * forward.z + forward.x * forward.x);
+ Vec3 right = {forward.z / norm, 0.0f, -forward.x / norm};
+
+ switch (key) {
+ case 'w': /* Move target forward */
+ g_cam.t.x += forward.x * moveSpeed;
+ g_cam.t.z += forward.z * moveSpeed;
+ break;
+ case 's': /* Move target backward */
+ g_cam.t.x -= forward.x * moveSpeed;
+ g_cam.t.z -= forward.z * moveSpeed;
+ break;
+ case 'a': /* Move target left */
+ g_cam.t.x += right.x * moveSpeed;
+ g_cam.t.z += right.z * moveSpeed;
+ break;
+ case 'd': /* Move target right */
+ g_cam.t.x -= right.x * moveSpeed;
+ g_cam.t.z -= right.z * moveSpeed;
+ break;
+
+ case 'q': /* Roll left */
+ g_cam.roll += tiltSpeed;
+ break;
+ case 'e': /* Roll right */
+ g_cam.roll -= tiltSpeed;
+ break;
+
+ case 'x': /* Increase g_curBounce */
+ if (g_bounce_cur < g_numBounces) g_bounce_cur++;
+ break;
+ case 'z': /* Decrease g_curBounce */
+ if (g_bounce_cur > 0) g_bounce_cur--;
+ break;
+ }
+
+ updateCamera();
+ glutPostRedisplay();
+}
+
+void mouse(int button, int state, int x, int y) {
+ if (button == GLUT_LEFT_BUTTON) {
+ g_mouseDown = (state == GLUT_DOWN);
+ g_cam.lastX = x;
+ g_cam.lastY = y;
+ }
+ /* Handle scroll wheel */
+ else if (button == 3) { /* Scroll up */
+ g_cam.d = fmaxf(0.1f, g_cam.d - 0.5f); /* Prevent going too close */
+ updateCamera();
+ glutPostRedisplay();
+ }
+ else if (button == 4) { /* Scroll down */
+ g_cam.d += 0.5f;
+ updateCamera();
+ glutPostRedisplay();
+ }
+}
+
+void mouseMotion(int x, int y) {
+ if (!g_mouseDown) return;
+
+ float sensitivity = 0.005f;
+ float dx = (float)(x - g_cam.lastX);
+ float dy = (float)(y - g_cam.lastY);
+
+ g_cam.yaw -= dx * sensitivity;
+ g_cam.pitch -= dy * sensitivity;
+
+ /* Clamp pitch to prevent flipping */
+ if (g_cam.pitch > 1.5f) g_cam.pitch = 1.5f;
+ if (g_cam.pitch < -1.5f) g_cam.pitch = -1.5f;
+
+ g_cam.lastX = x;
+ g_cam.lastY = y;
+
+ updateCamera();
+ glutPostRedisplay();
+}
+
+/**
+ * Main
+ */
+
+int main(int argc, char** argv) {
+ const char* rays_filename = "rays.bin";
+
+ /* Load ray data */
+ loadRays(rays_filename);
+ if (argc > 1)
+ g_scene = scene_load(argv[1]);
+
+ /* Initialize GLUT */
+ glutInit(&argc, argv);
+ glutInitDisplayMode(GLUT_DOUBLE | GLUT_RGB | GLUT_DEPTH);
+ glutInitWindowSize(800, 600);
+ glutCreateWindow("Ray Visualization");
+
+ /* Set up OpenGL */
+ glEnable(GL_DEPTH_TEST);
+ glClearColor(0.f, 0.f, 0.f, 1.0f);
+
+ /* Register callbacks */
+ glutDisplayFunc(display);
+ glutReshapeFunc(reshape);
+ glutKeyboardFunc(keyboard);
+ glutMouseFunc(mouse);
+ glutMotionFunc(mouseMotion);
+
+ updateCamera();
+ glutMainLoop();
+
+ return 0;
+}
diff --git a/test.c b/test.c
@@ -1,63 +0,0 @@
-#include "compute_paths.h"
-#include "test.h"
-
-#include <stdio.h> /* for fprintf */
-#include <stdlib.h> /* for malloc */
-
-int main(int argc, char **argv)
-{
- if (argc < 2) {
- fprintf(stderr, "Usage: %s <path_to_mesh.ply>\n", argv[0]);
- return 1;
- }
-
- float rx_positions[3] = {0.0, 0.0, .5};
- float tx_positions[3] = {0.0, 0.0, .5};
- float rx_velocities[3] = {0.0, 0.0, 0.0};
- float tx_velocities[3] = {0.0, 0.0, 0.0};
- float carrier_frequency = 3.0; /* 3 GHz */
- PathsInfo los = {
- .num_paths = 1,
- .directions_rx = (float*)malloc(g_numRx * g_numTx * 3 * sizeof(float)),
- .directions_tx = (float*)malloc(g_numRx * g_numTx * 3 * sizeof(float)),
- .a_te_re = (float*)malloc(g_numRx * g_numTx * sizeof(float)),
- .a_te_im = (float*)malloc(g_numRx * g_numTx * sizeof(float)),
- .a_tm_re = (float*)malloc(g_numRx * g_numTx * sizeof(float)),
- .a_tm_im = (float*)malloc(g_numRx * g_numTx * sizeof(float)),
- .tau = (float*)malloc(g_numRx * g_numTx * sizeof(float)),
- .freq_shift = (float*)malloc(g_numRx * g_numTx * sizeof(float))
- };
- PathsInfo scatter = {
- .num_paths = g_numBounces * g_numPaths,
- .directions_rx = (float*)malloc(g_numRx * g_numTx * g_numBounces * g_numPaths * 3 * sizeof(float)),
- .directions_tx = (float*)malloc(g_numPaths * 3 * sizeof(float)),
- .a_te_re = (float*)malloc(g_numRx * g_numTx * g_numBounces * g_numPaths * sizeof(float)),
- .a_te_im = (float*)malloc(g_numRx * g_numTx * g_numBounces * g_numPaths * sizeof(float)),
- .a_tm_re = (float*)malloc(g_numRx * g_numTx * g_numBounces * g_numPaths * sizeof(float)),
- .a_tm_im = (float*)malloc(g_numRx * g_numTx * g_numBounces * g_numPaths * sizeof(float)),
- .tau = (float*)malloc(g_numRx * g_numTx * g_numBounces * g_numPaths * sizeof(float)),
- .freq_shift = (float*)malloc(g_numRx * g_numTx * g_numBounces * g_numPaths * sizeof(float))
- };
- compute_paths(
- argv[1],
- rx_positions,
- tx_positions,
- rx_velocities,
- tx_velocities,
- carrier_frequency,
- g_numRx, g_numTx, g_numPaths, g_numBounces,
- &los, &scatter
- );
-
- /* Save the results */
-
- FILE *f;
- #define WRITE_BIN(name, data, size) \
- f = fopen(name, "wb"); \
- fwrite(data, sizeof(float), size, f); \
- fclose(f);
-
- /* Use this macro to write what you would like to inspect to a binary file */
- /* Example: */
- /* WRITE_BIN("los_directions_rx.bin", los.directions_rx, g_numRx * g_numTx * 3); */
-}
diff --git a/test/2cars.c b/test/2cars.c
@@ -0,0 +1,13 @@
+#include "../inc/test.h"
+#include "../inc/common.h"
+#include "../inc/compute_paths.h"
+
+int main(int argc, char *argv[]) {
+ char scene_filepath[] = "scenes/2cars.hrt";
+ float rx_pos[3] = {0.0, 0.0, 0.0};
+ float tx_pos[3] = {0.0, 0.0, 0.0};
+ float rx_vel[3] = {0.0, 0.0, 0.0};
+ float tx_vel[3] = {0.0, 0.0, 0.0};
+ float carrier_frequency = 70.f;
+
+}
diff --git a/test/test.c b/test/test.c
@@ -0,0 +1,66 @@
+#include "../inc/compute_paths.h" /* for compute_paths */
+#include "../inc/scene.h" /* for Scene, scene_load */
+
+#include "test.h" /* for g_numRx, g_numTx, g_numPaths, g_numBounces */
+
+#include <stdio.h> /* for fprintf */
+#include <stdlib.h> /* for malloc */
+
+int main(int argc, char **argv)
+{
+ if (argc < 2) {
+ fprintf(stderr, "Usage: %s <path_to_mesh.ply>\n", argv[0]);
+ return 1;
+ }
+
+ float rx_positions[3] = {0.0, 0.0, .5};
+ float tx_positions[3] = {0.0, 0.0, .5};
+ float rx_velocities[3] = {0.0, 0.0, 0.0};
+ float tx_velocities[3] = {0.0, 0.0, 0.0};
+ float carrier_frequency = 3.0; /* 3 GHz */
+ PathsInfo los = {
+ .num_paths = 1,
+ .directions_rx = (float*)malloc(g_numRx * g_numTx * 3 * sizeof(float)),
+ .directions_tx = (float*)malloc(g_numRx * g_numTx * 3 * sizeof(float)),
+ .a_te_re = (float*)malloc(g_numRx * g_numTx * sizeof(float)),
+ .a_te_im = (float*)malloc(g_numRx * g_numTx * sizeof(float)),
+ .a_tm_re = (float*)malloc(g_numRx * g_numTx * sizeof(float)),
+ .a_tm_im = (float*)malloc(g_numRx * g_numTx * sizeof(float)),
+ .tau = (float*)malloc(g_numRx * g_numTx * sizeof(float)),
+ .freq_shift = (float*)malloc(g_numRx * g_numTx * sizeof(float))
+ };
+ PathsInfo scatter = {
+ .num_paths = g_numBounces * g_numPaths,
+ .directions_rx = (float*)malloc(g_numRx * g_numTx * g_numBounces * g_numPaths * 3 * sizeof(float)),
+ .directions_tx = (float*)malloc(g_numPaths * 3 * sizeof(float)),
+ .a_te_re = (float*)malloc(g_numRx * g_numTx * g_numBounces * g_numPaths * sizeof(float)),
+ .a_te_im = (float*)malloc(g_numRx * g_numTx * g_numBounces * g_numPaths * sizeof(float)),
+ .a_tm_re = (float*)malloc(g_numRx * g_numTx * g_numBounces * g_numPaths * sizeof(float)),
+ .a_tm_im = (float*)malloc(g_numRx * g_numTx * g_numBounces * g_numPaths * sizeof(float)),
+ .tau = (float*)malloc(g_numRx * g_numTx * g_numBounces * g_numPaths * sizeof(float)),
+ .freq_shift = (float*)malloc(g_numRx * g_numTx * g_numBounces * g_numPaths * sizeof(float))
+ };
+ Scene scene = scene_load(argv[1]);
+ compute_paths(
+ &scene,
+ rx_positions,
+ tx_positions,
+ rx_velocities,
+ tx_velocities,
+ carrier_frequency,
+ g_numRx, g_numTx, g_numPaths, g_numBounces,
+ &los, &scatter
+ );
+
+ /* Save the results */
+
+ FILE *f;
+ #define WRITE_BIN(name, data, size) \
+ f = fopen(name, "wb"); \
+ fwrite(data, sizeof(float), size, f); \
+ fclose(f);
+
+ /* Use this macro to write what you would like to inspect to a binary file */
+ /* Example: */
+ /* WRITE_BIN("los_directions_rx.bin", los.directions_rx, g_numRx * g_numTx * 3); */
+}
diff --git a/inc/test.h b/test/test.h
diff --git a/test.py b/test/test.py
diff --git a/vizrays.c b/vizrays.c
@@ -1,394 +0,0 @@
-#include "test.h" /* for g_numPaths, g_numBounces, g_numTx, g_numRx */
-
-#include <GL/glut.h>
-#include <stdio.h>
-#include <stdlib.h>
-#include <math.h>
-#include <stdint.h>
-#include <ctype.h> /* for tolower */
-
-/**
- * Structs
- */
-
-typedef struct {
- float x, y, z;
-} Vec3;
-
-typedef struct {
- Vec3 o, t;
- float d;
- float lastX, lastY;
- float yaw, pitch;
-} Camera;
-
-typedef struct {
- uint32_t num_vertices;
- Vec3* vs;
- uint32_t num_triangles;
- uint32_t* is;
-} Mesh;
-
-typedef struct {
- uint32_t num_meshes;
- Mesh* meshes;
-} Scene;
-
-/**
- * Globals
- */
-
-Camera g_cam = {
- {0.f, 0.f, 5.f},
- {0.f, 0.f, 0.f},
- 5.f,
- 0.f, 0.f,
- 0.f, 0.f
-};
-uint8_t g_mouseDown = 0;
-
-float* g_rays = NULL;
-uint8_t *g_active = NULL;
-
-Scene g_scene = {0};
-
-int g_bounce_cur = 0;
-
-/**
- * Loading functions
- */
-
-void loadRays(const char* filename) {
- FILE* file = fopen(filename, "rb");
- if (!file) {
- printf("Failed to open file\n");
- exit(8);
- }
-
- uint32_t fileSize = (g_numBounces + 1) * g_numTx * g_numPaths * 2 * 3;
- g_rays = (float*)malloc(fileSize * sizeof(float));
- if (fread(g_rays, sizeof(float), fileSize, file) != fileSize) {
- printf("Failed to read rays\n");
- exit(8);
- }
- fclose(file);
-
- /* active.bin */
- file = fopen("active.bin", "rb");
- if (!file) {
- printf("Failed to open file\n");
- exit(8);
- }
- uint32_t activeSize = (g_numBounces + 1) * (g_numTx * g_numPaths / 8 + 1);
- g_active = (uint8_t*)malloc(activeSize * sizeof(uint8_t));
- if (fread(g_active, sizeof(uint8_t), activeSize, file) != activeSize) {
- printf("Failed to read active\n");
- exit(8);
- }
- fclose(file);
-}
-
-void loadScene(const char* filename) {
- FILE* file = fopen(filename, "rb");
- if (!file) {
- printf("Failed to open file\n");
- exit(8);
- }
-
- /* Magic */
- char magic[3];
- if (fread(magic, 1, 3, file) != 3) {
- printf("Failed to read magic\n");
- exit(8);
- }
- if (magic[0] != 'H' || magic[1] != 'R' || magic[2] != 'T') {
- printf("Invalid file format\n");
- exit(8);
- }
-
- /* Scene */
- if (fread(&g_scene.num_meshes, sizeof(uint32_t), 1, file) != 1) {
- printf("Failed to read num_meshes\n");
- exit(8);
- }
- g_scene.meshes = (Mesh*)malloc(g_scene.num_meshes * sizeof(Mesh));
-
- /* Meshes */
- for (Mesh* mesh = g_scene.meshes;
- mesh != g_scene.meshes + g_scene.num_meshes;
- ++mesh)
- {
- if (fread(&mesh->num_vertices, sizeof(uint32_t), 1, file) != 1) {
- printf("Failed to read num_vertices\n");
- exit(8);
- }
- mesh->vs = (Vec3*)malloc(mesh->num_vertices * sizeof(Vec3));
- if (fread(mesh->vs, sizeof(Vec3), mesh->num_vertices, file) != mesh->num_vertices) {
- printf("Failed to read vertices\n");
- exit(8);
- }
- if (fread(&mesh->num_triangles, sizeof(uint32_t), 1, file) != 1) {
- printf("Failed to read num_triangles\n");
- exit(8);
- }
- mesh->is = (uint32_t*)malloc(mesh->num_triangles * 3 * sizeof(uint32_t));
- if (fread(mesh->is, sizeof(uint32_t), mesh->num_triangles * 3, file) != mesh->num_triangles * 3) {
- printf("Failed to read indices\n");
- exit(8);
- }
- fseek(file, sizeof(uint32_t), SEEK_CUR); /* Skip material index */
- }
- fclose(file);
-}
-
-/**
- * Draw functions
- */
-
-void drawScene() {
- glColor3f(1.0f, 1.0f, 1.0f);
- glBegin(GL_TRIANGLES);
- for (uint32_t i = 0; i < g_scene.num_meshes; i++) {
- Mesh* mesh = &g_scene.meshes[i];
- float mesh_color = (float)i / g_scene.num_meshes;
- glColor3f(0.f, mesh_color, 1.f - mesh_color);
-
- for (uint32_t j = 0; j < mesh->num_triangles; j++) {
- uint32_t idx = mesh->is[j * 3];
- Vec3 v1 = mesh->vs[idx];
- idx = mesh->is[j * 3 + 1];
- Vec3 v2 = mesh->vs[idx];
- idx = mesh->is[j * 3 + 2];
- Vec3 v3 = mesh->vs[idx];
- glVertex3f(v1.x, v1.y, v1.z);
- glVertex3f(v2.x, v2.y, v2.z);
- glVertex3f(v3.x, v3.y, v3.z);
- }
- }
- glEnd();
-}
-
-void drawRays() {
- float color = (float)g_bounce_cur / g_numBounces;
- glColor3f(color, 0.0f, 1.0f - color);
-
- glBegin(GL_LINES);
- int idx_base = g_bounce_cur * g_numTx * g_numPaths * 2 * 3;
- uint32_t active_byte = g_bounce_cur * (g_numTx * g_numPaths / 8 + 1);
- uint8_t active_bit = 1;
- uint32_t num_skipped = 0;
- for (int path = 0; path < g_numPaths; ++path, active_bit <<= 1) {
- int idx = idx_base + path * 2 * 3;
-
- if (!active_bit) {
- active_byte++;
- active_bit = 1;
- }
- if (!(g_active[active_byte] & active_bit))
- {
- num_skipped++;
- continue;
- }
-
- /* Origin */
- float x1 = g_rays[idx];
- float y1 = g_rays[idx + 1];
- float z1 = g_rays[idx + 2];
-
- /* Destination */
- float x2 = x1 + g_rays[idx + 3];
- float y2 = y1 + g_rays[idx + 4];
- float z2 = z1 + g_rays[idx + 5];
-
- glVertex3f(x1, y1, z1);
- glVertex3f(x2, y2, z2);
- }
- glEnd();
- printf("\rSkipped %d paths", num_skipped);
- fflush(stdout);
-
- /* Draw points */
- glPointSize(5.0f);
- glBegin(GL_POINTS);
- active_byte = g_bounce_cur * (g_numTx * g_numPaths / 8 + 1);
- active_bit = 1;
- for (int path = 0; path < g_numPaths; ++path) {
- int idx = idx_base + path * 2 * 3;
- if (!active_bit) {
- active_byte++;
- active_bit = 1;
- }
- if (!(g_active[active_byte] & active_bit))
- continue;
- glVertex3f(g_rays[idx], g_rays[idx + 1], g_rays[idx + 2]);
- }
- glEnd();
-}
-
-/**
- * Callbacks
- */
-
-void reshape(int w, int h) {
- glViewport(0, 0, w, h);
- glMatrixMode(GL_PROJECTION);
- glLoadIdentity();
- gluPerspective(60.0f, (float)w/h, 0.1f, 1000.0f);
- glMatrixMode(GL_MODELVIEW);
-}
-
-void display() {
- glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
- glLoadIdentity();
-
- /* Look at target from camera position */
- gluLookAt(
- g_cam.o.x, g_cam.o.y, g_cam.o.z,
- g_cam.t.x, g_cam.t.y, g_cam.t.z,
- 0.0f, 1.0f, 0.0f
- );
-
- drawScene();
- drawRays();
-
- glutSwapBuffers();
-}
-
-void updateCamera() {
- float x = g_cam.d * cosf(g_cam.pitch) * cosf(g_cam.yaw);
- float y = g_cam.d * sinf(g_cam.pitch);
- float z = g_cam.d * cosf(g_cam.pitch) * sinf(g_cam.yaw);
-
- g_cam.o.x = g_cam.t.x + x;
- g_cam.o.y = g_cam.t.y + y;
- g_cam.o.z = g_cam.t.z + z;
-
- gluLookAt(g_cam.o.x, g_cam.o.y, g_cam.o.z,
- g_cam.t.x, g_cam.t.y, g_cam.t.z,
- 0.0f, 1.0f, 0.0f);
-}
-
-void keyboard(unsigned char key, int x, int y) {
- float speed = 0.1f;
- /* Speed up if SHIFT is down */
- if (glutGetModifiers() & GLUT_ACTIVE_SHIFT) {
- speed *= 5;
- key = tolower(key);
- }
-
- /* Calculate forward direction */
- Vec3 dir = {
- g_cam.t.x - g_cam.o.x,
- g_cam.t.y - g_cam.o.y,
- g_cam.t.z - g_cam.o.z
- };
- float norm = sqrtf(dir.x * dir.x + dir.y * dir.y + dir.z * dir.z);
- Vec3 forward = {dir.x / norm, dir.y / norm, dir.z / norm};
-
- /* Calculate right vector (cross product of forward and up (0,0,1)) */
- norm = sqrtf(forward.z * forward.z + forward.x * forward.x);
- Vec3 right = {forward.z / norm, 0.0f, -forward.x / norm};
-
- switch (key) {
- case 'w': /* Move target forward */
- g_cam.t.x += forward.x * speed;
- g_cam.t.z += forward.z * speed;
- break;
- case 's': /* Move target backward */
- g_cam.t.x -= forward.x * speed;
- g_cam.t.z -= forward.z * speed;
- break;
- case 'a': /* Move target left */
- g_cam.t.x += right.x * speed;
- g_cam.t.z += right.z * speed;
- break;
- case 'd': /* Move target right */
- g_cam.t.x -= right.x * speed;
- g_cam.t.z -= right.z * speed;
- break;
- case 'x': /* Increase g_curBounce */
- if (g_bounce_cur < g_numBounces) g_bounce_cur++;
- break;
- case 'z': /* Decrease g_curBounce */
- if (g_bounce_cur > 0) g_bounce_cur--;
- break;
- }
-
- updateCamera();
- glutPostRedisplay();
-}
-
-void mouse(int button, int state, int x, int y) {
- if (button == GLUT_LEFT_BUTTON) {
- g_mouseDown = (state == GLUT_DOWN);
- g_cam.lastX = x;
- g_cam.lastY = y;
- }
- /* Handle scroll wheel */
- else if (button == 3) { /* Scroll up */
- g_cam.d = fmaxf(0.1f, g_cam.d - 0.5f); /* Prevent going too close */
- updateCamera();
- glutPostRedisplay();
- }
- else if (button == 4) { /* Scroll down */
- g_cam.d += 0.5f;
- updateCamera();
- glutPostRedisplay();
- }
-}
-
-void mouseMotion(int x, int y) {
- if (!g_mouseDown) return;
-
- float sensitivity = 0.005f;
- float dx = (float)(x - g_cam.lastX);
- float dy = (float)(y - g_cam.lastY);
-
- g_cam.yaw -= dx * sensitivity;
- g_cam.pitch -= dy * sensitivity;
-
- /* Clamp pitch to prevent flipping */
- if (g_cam.pitch > 1.5f) g_cam.pitch = 1.5f;
- if (g_cam.pitch < -1.5f) g_cam.pitch = -1.5f;
-
- g_cam.lastX = x;
- g_cam.lastY = y;
-
- updateCamera();
- glutPostRedisplay();
-}
-
-/**
- * Main
- */
-
-int main(int argc, char** argv) {
- const char* rays_filename = "rays.bin";
-
- /* Load ray data */
- loadRays(rays_filename);
- if (argc > 1)
- loadScene(argv[1]);
-
- /* Initialize GLUT */
- glutInit(&argc, argv);
- glutInitDisplayMode(GLUT_DOUBLE | GLUT_RGB | GLUT_DEPTH);
- glutInitWindowSize(800, 600);
- glutCreateWindow("Ray Visualization");
-
- /* Set up OpenGL */
- glEnable(GL_DEPTH_TEST);
- glClearColor(0.f, 0.f, 0.f, 1.0f);
-
- /* Register callbacks */
- glutDisplayFunc(display);
- glutReshapeFunc(reshape);
- glutKeyboardFunc(keyboard);
- glutMouseFunc(mouse);
- glutMotionFunc(mouseMotion);
-
- updateCamera();
- glutMainLoop();
-
- return 0;
-}
