commit bdddf1785101f45288aed7555f66f13fbb9103d1
parent cbe54e2e4f21cc541f89c2824180cee2b122c038
Author: Egor Achkasov <eaachkasov@edu.hse.ru>
Date: Wed, 5 Feb 2025 14:16:24 +0100
Impl new scene format to support volumes
Diffstat:
9 files changed, 422 insertions(+), 439 deletions(-)
diff --git a/compute_paths.c b/compute_paths.c
@@ -1,4 +1,5 @@
-#include "compute_paths.h"
+#include "compute_paths.h" /* for compute_paths */
+#include "scene.h" /* for HRT_Scene, HRT_mesh, HRT_Material */
#include <stddef.h> /* for size_t */
#include <stdlib.h> /* for exit, malloc, free */
@@ -113,11 +114,6 @@ const float* BINOMIAL_ALPHA_K[20] = {
}
};
-/* ==== PRECOMPUTED GLOBALS ==== */
-
-/* 4.f * PI * carrier_frequency * 1e9 / SPEED_OF_LIGHT */
-float g_free_space_loss_multiplier;
-
/* ==== STRUCTS ==== */
typedef struct {
@@ -128,35 +124,38 @@ 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;
- /* scattering */
- float s; /* scattering coefficient */
- float r; /* reflection reduction factor (=sqrt(1-s^2)) */
- float alpha; /* directive pattern lobe width parameter */
-} RadioMaterial;
+ /* r = 1.0 - s */
+ float r;
+} Material;
typedef struct {
- Vec3 *vertices;
- size_t num_vertices;
- RadioMaterial *rms;
- size_t num_rms;
- int32_t *indices;
- int32_t *rm_indices;
- size_t num_indices;
- Vec3 *normals;
+ uint32_t num_vertices;
+ Vec3 *vs;
+ uint32_t num_triangles;
+ uint32_t *is;
+ /* Normals. Size [num_triangles] */
+ Vec3 *ns;
+ uint32_t material_index;
} Mesh;
-void free_mesh(Mesh *mesh)
-{
- free(mesh->vertices);
- free(mesh->rms);
- free(mesh->indices);
- free(mesh->rm_indices);
- free(mesh->normals);
-}
+
+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 ==== */
@@ -222,179 +221,94 @@ Vec3 vec3_normalize(const Vec3 *a)
return (Vec3){a->x / norm, a->y / norm, a->z / norm};
}
-/* ==== MESH LOADING ==== */
+/* ==== SCENE LOADING ==== */
-/** Load a mesh from a PLY file.
- *
- * ply
- * format binary_little_endian 1.0
- * element vertex %d\n", num_vertices);
- * property float x
- * property float y
- * property float z
- * element radio_material %d\n", num_materials);
- * property list uint char name
- * property float a
- * property float b
- * property float c
- * property float d
- * property float s
- * property char alpha
- * element face %d\n", num_faces);
- * property list uchar uint vertex_index
- * property uint material_index
- * end_header
+/** Load a mesh from a HRT file.
*
- * All the faces must be triangles.
- *
- * \param mesh_filepath path to the PLY file
+ * \param scene_filepath path to the HRT file
* \param carrier_frequency the carrier frequency in GHz
- * \return the loaded mesh
+ * \return the loaded scene
*/
-Mesh load_mesh_ply(const char *mesh_filepath, float carrier_frequency)
+Scene load_scene(const char *scene_filepath, float carrier_frequency)
{
- FILE *f = fopen(mesh_filepath, "rb");
+ /* Open the HRT file */
+ FILE *f = fopen(scene_filepath, "rb");
if (!f) {
- fprintf(stderr, "Could not open file %s\n", mesh_filepath);
+ fprintf(stderr, "Could not open file %s\n", scene_filepath);
exit(8);
}
- char buff[128];
- size_t num_vertices = 0;
- size_t num_faces = 0;
- size_t num_materials = 0;
-
- /* HEADER */
- /* ply */
- if (fread(buff, 1, 4, f) != 4) exit(8);
- if (strncmp(buff, "ply\n", 4)) exit(8);
- /* format binary_little_endian 1.0 */
- if (fread(buff, 1, 32, f) != 32) exit(8);
- if (strncmp(buff, "format binary_little_endian 1.0\n", 32)) exit(8);
-
- /* element vertex <num_vertices> */
- if (fgets(buff, 128, f) == NULL) exit(8);
- if (strncmp(buff, "element vertex ", 15)) exit(8);
- sscanf(buff, "element vertex %zu", &num_vertices);
- /* property float x */
- /* property float y */
- /* property float z */
- if (fread(buff, 1, 17*3, f) != 17*3) exit(8);
- if (strncmp(buff, "property float x\nproperty float y\nproperty float z\n", 17*3))
- exit(8);
-
- /* element radio_material */
- if (fgets(buff, 128, f) == NULL) exit(8);
- if (strncmp(buff, "element radio_material ", 23)) exit(8);
- sscanf(buff, "element radio_material %zu", &num_materials);
- /* property list uint char name */
- if (fread(buff, 1, 29, f) != 29) exit(8);
- if (strncmp(buff, "property list uint char name\n", 29)) exit(8);
- /* property float a */
- /* property float b */
- /* property float c */
- /* property float d */
- if (fread(buff, 1, 17*4, f) != 17*4) exit(8);
- if (strncmp(buff, "property float a\nproperty float b\nproperty float c\nproperty float d\n", 17*4))
- exit(8);
- /* property float s */
- /* property char alpha */
- if (fread(buff, 1, 17*2+3, f) != 17*2+3) exit(8);
- if (strncmp(buff, "property float s\nproperty char alpha\n", 17*2+3))
- exit(8);
-
- /* element face <num_faces> */
- if (fgets(buff, 128, f) == NULL) exit(8);
- if (strncmp(buff, "element face ", 13)) exit(8);
- sscanf(buff, "element face %zu", &num_faces);
- /* property list uchar int vertex_index */
- if (fread(buff, 1, 38, f) != 38) exit(8);
- if (strncmp(buff, "property list uchar uint vertex_index\n", 38))
- exit(8);
- /* property uint material_index */
- if (fread(buff, 1, 29, f) != 29) exit(8);
- if (strncmp(buff, "property uint material_index\n", 29)) exit(8);
-
- /* end_header */
- if (fread(buff, 1, 11, f) != 11) exit(8);
- if (strncmp(buff, "end_header\n", 11)) exit(8);
-
- /* Init mesh */
- Mesh mesh;
- mesh.num_vertices = num_vertices;
- mesh.num_rms = num_materials;
- mesh.num_indices = num_faces * 3;
- mesh.vertices = (Vec3*)malloc(num_vertices * sizeof(Vec3));
- mesh.rms = (RadioMaterial*)malloc(num_materials * sizeof(RadioMaterial));
- mesh.indices = (int32_t*)malloc(mesh.num_indices * sizeof(int32_t));
- mesh.rm_indices = (int32_t*)malloc(num_faces * sizeof(int32_t));
- mesh.normals = (Vec3*)malloc(num_faces * sizeof(Vec3));
-
- /* VERTICES */
- for (size_t i = 0; i < mesh.num_vertices; ++i)
- if (fread(&mesh.vertices[i], sizeof(Vec3), 1, f) != 1)
- exit(8);
-
- /* RADIO_MATERIALS */
- /* read a, b, c, d, s, alpha
- and calculate relative permitivity */
- unsigned int name_size;
- float abcd[4];
- RadioMaterial *rm;
- for (size_t i = 0; i < num_materials; ++i) {
- rm = &mesh.rms[i];
- /* skip name */
- if (fread(&name_size, sizeof(unsigned int), 1, f) != 1) exit(8);
- fseek(f, name_size, SEEK_CUR);
- /* read a, b, c and d */
- if (fread(abcd, sizeof(float)*4, 1, f) != 1) exit(8);
- /* read s and alpha */
- if (fread(&rm->s, sizeof(float), 1, f) != 1) exit(8);
- if (fread(&rm->alpha, sizeof(uint8_t), 1, f) != 1) exit(8);
- /* calculate r */
- rm->r = sqrtf(1.f - rm->s * rm->s);
- /* calculate eta */
- rm->eta_re = abcd[0] * powf(carrier_frequency, abcd[1]);
- /* eq. 12 */
- rm->eta_im = (abcd[2] * powf(carrier_frequency, abcd[3]))
- / (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);
+ /* Parse the file */
+ /* MAGIC */
+ char magic[3];
+ fread(magic, 1, 3, f);
+ if (strncmp(magic, "HRT", 3)) exit(8);
+ /* SCENE */
+ Scene scene;
+ /* num_meshes */
+ fread(&scene.num_meshes, sizeof(uint32_t), 1, f);
+ /* 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 */
+ fread(&mesh->num_vertices, sizeof(uint32_t), 1, f);
+ /* vertices */
+ mesh->vs = (Vec3*)malloc(mesh->num_vertices * sizeof(Vec3));
+ fread(mesh->vs, sizeof(Vec3), mesh->num_vertices, f);
+ /* num_triangles */
+ fread(&mesh->num_triangles, sizeof(uint32_t), 1, f);
+ /* triangles */
+ mesh->is = (uint32_t*)malloc(mesh->num_triangles * 3 * sizeof(uint32_t));
+ fread(mesh->is, sizeof(uint32_t), mesh->num_triangles * 3, f);
+ /* material_index */
+ fread(&mesh->material_index, sizeof(uint32_t), 1, f);
}
+ fclose(f);
- /* INDICES */
- /* also calculate normals */
- uint8_t n;
- Vec3 *v1, *v2, *v3, u, v;
- for (size_t i = 0; i < mesh.num_indices; i += 3) {
- if (fread(&n, 1, 1, f) != 1) exit(8);
- if (n != 3) {
- fprintf(stderr, "Only trianglar faces are supported\n");
- exit(8);
+ /* 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]);
}
- if (fread(&mesh.indices[i], sizeof(int32_t), 3, f) != 3) exit(8);
- if (fread(&mesh.rm_indices[i / 3], sizeof(int32_t), 1, f) != 1) exit(8);
- /* calculate normal */
- v1 = &mesh.vertices[mesh.indices[i]];
- v2 = &mesh.vertices[mesh.indices[i + 1]];
- v3 = &mesh.vertices[mesh.indices[i + 2]];
- u = vec3_sub(v2, v1);
- v = vec3_sub(v3, v1);
- mesh.normals[i / 3] = vec3_cross(&u, &v);
- mesh.normals[i / 3] = vec3_normalize(&mesh.normals[i / 3]);
}
- fclose(f);
- return mesh;
+ return scene;
}
/* ==== RT ==== */
@@ -402,16 +316,18 @@ Mesh load_mesh_ply(const char *mesh_filepath, float carrier_frequency)
/** Compute Moeleer-Trumbore intersection algorithm.
*
* \param ray ray to cast to the mesh
- * \param mesh triangle 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 ind output index of the hit triangle. If no hit, i 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 Mesh *mesh,
+ IN Scene *scene,
OUT float *t,
- OUT int32_t *ind,
+ OUT uint32_t *mesh_ind,
+ OUT uint32_t *face_ind,
OUT float *theta
)
{
@@ -421,33 +337,37 @@ void moeller_trumbore(
float d, u, v;
float dist = 1e9;
float dist_tmp;
- for (size_t i = 0; i < mesh->num_indices; i += 3) {
- v1 = mesh->vertices[mesh->indices[i]];
- v2 = mesh->vertices[mesh->indices[i + 1]];
- v3 = mesh->vertices[mesh->indices[i + 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;
- *ind = i / 3;
- *theta = acos(vec3_dot(&mesh->normals[i / 3], &ray->d));
- if (*theta > PI / 2.)
- *theta = PI - *theta;
+ 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;
+ }
}
}
}
@@ -456,7 +376,7 @@ void moeller_trumbore(
*
* Implements eqs. (31a)-(31b) from ITU-R P.2040-3.
*
- * \param rm the radio material of the hit triangle
+ * \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}
@@ -464,7 +384,7 @@ void moeller_trumbore(
* \param r_tm_im output imaginary part of R_{eTM}
*/
void refl_coefs(
- IN RadioMaterial *rm,
+ IN uint32_t material_index,
IN float theta1,
OUT float *r_te_re,
OUT float *r_te_im,
@@ -472,6 +392,7 @@ void refl_coefs(
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;
@@ -510,14 +431,11 @@ void refl_coefs(
/** Calculate scattering coefficients.
*
- * Directive model pattern is used.
- *
* Implements eqs. TODO from doi 10.1109/TAP.2006.888422
*
* \param theta_s scattering angle
* \param theta_i angle of incidence
- * \param s scattering coefficient
- * \param alpha directive pattern lobe width parameter
+ * \param material_index the index of the material
* \param a_te_re output real part of a_{eTE}
* \param a_te_im output imaginary part of a_{eTE}
* \param a_tm_re output real part of a_{eTM}
@@ -526,8 +444,7 @@ void refl_coefs(
void scat_coefs(
IN float theta_s,
IN float theta_i,
- IN float s,
- IN uint8_t alpha,
+ IN uint32_t material_index,
OUT float *a_te_re,
OUT float *a_te_im,
OUT float *a_tm_re,
@@ -538,7 +455,9 @@ void scat_coefs(
float theta_i_sin = sinf(theta_i);
float theta_i_cos = cosf(theta_i);
float f, I_k;
- for (uint8_t k = 0; k <= alpha; ++k) {
+ /* Directive model pattern */
+ HRT_Material *hrt_mat = &g_hrt_materials[material_index];
+ for (uint8_t k = 0; k <= hrt_mat->s1_alpha; ++k) {
if (k % 2) {
I_k = 0.f;
for (uint8_t w = 0; w <= (k-1)/2; ++w)
@@ -547,17 +466,18 @@ void scat_coefs(
} else
I_k = 1.f;
I_k *= 2.f * PI / (k + 1);
- F_alpha += BINOMIAL_ALPHA_K[alpha][k] * I_k;
+ F_alpha += BINOMIAL_ALPHA_K[hrt_mat->s1_alpha][k] * I_k;
}
- F_alpha /= powf(2.f, alpha);
- f = powf(.5f + theta_i_cos / 2.f, alpha) / F_alpha;
- *a_te_re = *a_te_im = *a_tm_re = *a_tm_im = s * f;
+ F_alpha /= powf(2.f, hrt_mat->s1_alpha);
+ f = powf(.5f + theta_i_cos / 2.f, hrt_mat->s1_alpha) / F_alpha;
+ *a_te_re = *a_te_im = *a_tm_re = *a_tm_im = hrt_mat->s * f;
+ /* TODO add s2 and s3 */
}
/* ==== MAIN FUNCTION ==== */
void compute_paths(
- IN const char *mesh_filepath, /* path to the mesh file */
+ 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) */
@@ -588,7 +508,7 @@ void compute_paths(
size_t off, off_scat;
/* Load the scene */
- Mesh mesh = load_mesh_ply(mesh_filepath, carrier_frequency);
+ Scene scene = load_scene(scene_filepath, carrier_frequency);
/* Precompute globals */
/* g_free_space_loss_multiplier */
@@ -649,7 +569,7 @@ void compute_paths(
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 */
- int32_t ind;
+ uint32_t mesh_ind, face_ind;
Ray r;
Vec3 *h = (Vec3*)malloc(sizeof(Vec3));
Vec3 n;
@@ -672,7 +592,7 @@ void compute_paths(
for (j = 0; j != num_tx; ++j) {
off = i * num_tx + j;
t = -1.f;
- ind = -1;
+ mesh_ind = face_ind = -1;
/* Create a ray from the tx to the rx */
r.o = tx_pos_v[j];
@@ -692,8 +612,8 @@ void compute_paths(
}
/* Is there any abstacle between the tx and the rx? */
- moeller_trumbore(&r, &mesh, &t, &ind, &theta);
- if (ind != -1 && t <= 1.f) {
+ moeller_trumbore(&r, &scene, &t, &mesh_ind, &face_ind, &theta);
+ if (mesh_ind != -1 && t <= 1.f) {
/* An obstacle between the tx and the rx has been hit */
a_te_re_los[off] = a_tm_re_los[off] = tau_los[off] = 0.f;
continue;
@@ -737,18 +657,16 @@ void compute_paths(
continue;
/* Init */
t = -1.f;
- ind = -1;
+ mesh_ind = face_ind = -1;
/* Find the hit point and trinagle and the angle of incidence */
- moeller_trumbore(&rays[off], &mesh, &t, &ind, &theta);
- if (ind == -1) { /* Ray hit nothing */
+ moeller_trumbore(&rays[off], &scene, &t, &mesh_ind, &face_ind, &theta);
+ if (mesh_ind == -1) { /* Ray hit nothing */
active[off / 8] &= ~(1 << (off % 8));
--num_active;
continue;
}
- /* Get the hit primitive's normal */
- n = mesh.normals[ind];
/* Calculate the reflection coefficients R_{eTE} and R_{eTM} */
- refl_coefs(&mesh.rms[mesh.rm_indices[ind]],
+ refl_coefs(scene.meshes[mesh_ind].material_index,
theta,
&r_te_re, &r_te_im,
&r_tm_re, &r_tm_im);
@@ -778,6 +696,7 @@ void compute_paths(
*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);
@@ -796,9 +715,9 @@ void compute_paths(
off_scat = i * num_rx * num_tx * num_paths + j * num_tx * num_paths + off;
r.d = vec3_sub(&rx_pos_v[j], &r.o);
r.d = vec3_normalize(&r.d);
- ind = -1;
- moeller_trumbore(&r, &mesh, &t, &ind, &theta);
- if (ind != -1 && t <= 1.f) {
+ mesh_ind = face_ind = -1;
+ moeller_trumbore(&r, &scene, &t, &mesh_ind, &face_ind, &theta);
+ if (mesh_ind != -1 && t <= 1.f) {
/* An obstacle between the hit point and the rx has been hit */
a_te_re_scat[off_scat] = a_te_im_scat[off_scat] = a_tm_re_scat[off_scat]
= a_tm_im_scat[off_scat] = tau_scat[off_scat] = 0.f;
@@ -812,11 +731,9 @@ void compute_paths(
/* Calculate the scattering angle */
theta_s = acosf(vec3_dot(&r.d, &n) / sqrtf(vec3_dot(&r.d, &r.d)));
/* Calculate the scattering coefficients */
- scat_coefs(theta_s, theta,
- mesh.rms[mesh.rm_indices[ind]].s,
- mesh.rms[mesh.rm_indices[ind]].alpha,
- &a_te_re_new, &a_te_im_new,
- &a_tm_re_new, &a_tm_im_new);
+ scat_coefs(theta_s, theta, scene.meshes[mesh_ind].material_index,
+ &a_te_re_new, &a_te_im_new,
+ &a_tm_re_new, &a_tm_im_new);
/* Calculate the distance */
t = sqrtf(vec3_dot(&r.d, &r.d));
/* Calculate the delay */
@@ -847,5 +764,6 @@ void compute_paths(
free(h);
free(active);
free(rays);
- free_mesh(&mesh);
+ /* Free the scene */
+ /* TODO */
}
diff --git a/compute_paths.h b/compute_paths.h
@@ -6,11 +6,14 @@
#define IN
#define OUT
-/** Compute gains and delays between in a Mitsuba scene.
+/** Compute gains and delays between tx and rx in a 3D scene.
*
- * Scene must be defined in a specific PLY format. See README for details.
+ * Scene must be defined in a specific format.
+ * The format is defined in scene.h.
+ * The HRT is a binary file cointaining a dereferenced Scene structure.
+ * See README for details.
*
- * \param mesh_filepath path to a mesh .ply file
+ * \param scene_filepath path to a scene .hrt file
* \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)
@@ -40,7 +43,7 @@
* \param tau_scat output array of delays for scatter in seconds, shape (num_bounces, num_rx, num_tx, num_paths)
*/
void compute_paths(
- IN const char *mesh_filepath, /* path to the mesh file */
+ 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) */
diff --git a/scene.h b/scene.h
@@ -0,0 +1,170 @@
+#include <stdint.h> /* for uint32_t */
+
+/*******************************************************************/
+/* Scene */
+/*******************************************************************/
+
+typedef struct {
+ /* Number of vertices */
+ uint32_t num_vertices;
+ /* Vertices coordinates. Size [num_vertices * 3] */
+ float *vs;
+ /* Number of triangles */
+ uint32_t num_triangles;
+ /* Triangles indices. Size [num_triangles * 3] */
+ uint32_t *is;
+ uint32_t material_index;
+} HRT_Mesh;
+
+typedef struct {
+ uint32_t num_meshes;
+ HRT_Mesh *meshes;
+} HRT_Scene;
+
+/*******************************************************************/
+/* Materials */
+/*******************************************************************/
+
+typedef struct {
+ /* Number of characters in the name */
+ uint32_t name_sz;
+ /* Name of the material. Not null-terminated. */
+ char *name;
+ /* Relative permitivity and conductivity properties.
+ * Refer to ITU-R P.2040-3 Table 3.
+ */
+ float a, b, c, d;
+ /* Scattering coefficient.
+ * Ratio between specular and diffuse reflection power.
+ * Must be in [0, 1].
+ */
+ float s;
+ /* Scattering pattern distribution ratios.
+ * Each must be in [0, 1].
+ * The sum of all ratios must be 1.
+ * s1: around specular (Directive Model)
+ * s2: around surface normal (Lambertian Model)
+ * s3: around incident direction (Backward Lobe Model)
+ */
+ float s1, s2, s3;
+ /* Directive model parameters.
+ * s1_alpha: lobe width integer parameter (TODO ref eq).
+ * Must be > 0.
+ */
+ uint8_t s1_alpha;
+ /* Backward lobe model parameters.
+ * s3_alpha: lobe width integer parameter (TODO ref eq).
+ * Must be > 0.
+ */
+ uint8_t s3_alpha;
+} HRT_Material;
+
+/* 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}
+};
+typedef enum {
+ MATERIAL_AIR = 0,
+ MATERIAL_CONCRETE = 1,
+ MATERIAL_BRICK = 2,
+ MATERIAL_PLASTERBOARD = 3,
+ MATERIAL_WOOD = 4,
+ MATERIAL_GLASS1 = 5,
+ MATERIAL_GLASS2 = 6,
+ MATERIAL_CEILING_BOARD1 = 7,
+ MATERIAL_CEILING_BOARD2 = 8,
+ MATERIAL_CHIPBOARD = 9,
+ MATERIAL_PLYWOOD = 10,
+ MATERIAL_MARBLE = 11,
+ MATERIAL_FLOORBOARD = 12,
+ MATERIAL_METAL = 13,
+ MATERIAL_VERY_DRY_GROUND = 14,
+ MATERIAL_MEDIUM_DRY_GROUND = 15,
+ MATERIAL_WET_GROUND = 16
+} MaterialIndex;
+
diff --git a/scene_fromSionna.c b/scene_fromSionna.c
@@ -0,0 +1,75 @@
+/* vim: set tabstop=2:softtabstop=2:shiftwidth=2:noexpandtab */
+
+#include "scene.h" /* for HRT_Scene, HRT_Mesh, HRT_Material */
+
+#include <stdio.h> /* for FILE, fopen, fclose, fseek, fgets, sscanf, printf */
+#include <stdlib.h> /* for malloc */
+#include <stdint.h> /* for uint8_t, uint32_t */
+
+int main(int argc, char *argv[]) {
+ if (argc != 2) {
+ printf("Usage: %s <scene.xml>\n", argv[0]);
+ return 1;
+ }
+
+ /* TODO: parse the XML file */
+ HRT_Mesh mesh_box;
+ mesh_box.num_vertices = 8;
+ float vs_temp[] = {
+ 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,
+ };
+ mesh_box.vs = (float*)vs_temp;
+ mesh_box.num_triangles = 12;
+ uint32_t is_temp[] = {
+ 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_box.is = (uint32_t*)is_temp;
+ mesh_box.material_index = MATERIAL_CONCRETE;
+
+ HRT_Scene scene;
+ scene.num_meshes = 1;
+ scene.meshes = &mesh_box;
+
+ FILE *fp = fopen("scene.hrt", "wb");
+ if (fp == NULL) {
+ printf("Error: cannot open file\n");
+ return 1;
+ }
+
+ /* 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 */
+ fwrite(&scene.meshes[i].num_vertices, sizeof(uint32_t), 1, fp);
+ fwrite(scene.meshes[i].vs, sizeof(float), 3 * scene.meshes[i].num_vertices, fp);
+ fwrite(&scene.meshes[i].num_triangles, sizeof(uint32_t), 1, fp);
+ fwrite(scene.meshes[i].is, sizeof(uint32_t), 3 * scene.meshes[i].num_triangles, fp);
+ fwrite(&scene.meshes[i].material_index, sizeof(uint32_t), 1, fp);
+ }
+
+ fclose(fp);
+
+ return 0;
+}
diff --git a/script/create_material_bin.c b/script/create_material_bin.c
@@ -1,36 +0,0 @@
-#include "material.h" /* for Material */
-
-#include <stdio.h>
-#include <stdint.h> /* for uint8_t */
-
-Material mats[2] = {
- {5, "metal", 1.f, 0.f, 10000000.f, 0.f, 0.f, 0},
- {8, "concrete", 5.24f, 0.f, 0.0462f, 0.7822f, 0.5f, 4}
-};
-
-/* Configure your scene's materials here */
-#define NUM_MATERIALS 1 /* number of materials in the scene */
-#define NUM_INDICES 6 /* number of material indices in the scene (= number of meshes) */
-int material_indices[NUM_INDICES] = {1, 1, 1, 1, 1, 1};
-
-int main() {
- FILE *fp = fopen("materials.bin", "wb");
- int num_materials = NUM_MATERIALS;
- int num_indices = NUM_INDICES;
- fwrite(&num_materials, sizeof(int), 1, fp);
- fwrite(&num_indices, sizeof(int), 1, fp);
- fwrite(material_indices, sizeof(int), NUM_INDICES, fp);
- for (int i = 0; i < NUM_MATERIALS; ++i) {
- fwrite(&mats[i].name_size, sizeof(int), 1, fp);
- fwrite(mats[i].name, sizeof(char), mats[i].name_size, fp);
- fwrite(&mats[i].a, sizeof(float), 1, fp);
- fwrite(&mats[i].b, sizeof(float), 1, fp);
- fwrite(&mats[i].c, sizeof(float), 1, fp);
- fwrite(&mats[i].d, sizeof(float), 1, fp);
- fwrite(&mats[i].s, sizeof(float), 1, fp);
- fwrite(&mats[i].alpha, sizeof(uint8_t), 1, fp);
- }
- fclose(fp);
-
- return 0;
-}
diff --git a/script/create_ply.c b/script/create_ply.c
@@ -1,138 +0,0 @@
-/* vim: set tabstop=2:softtabstop=2:shiftwidth=2:noexpandtab */
-
-#include "material.h" /* for Material */
-
-#include <stdio.h>
-#include <stdlib.h>
-#include <stdint.h>
-
-typedef struct {
- unsigned int num_vertices;
- float *vertices;
- unsigned int num_faces;
- int *faces;
-} Mesh;
-
-int main(int argc, char *argv[]) {
- if (argc < 3) {
- printf("Usage: %s <materials.bin> <path_to_mesh1.ply> [path_to_mesh2.ply] ...\n", argv[0]);
- return 1;
- }
-
- /* materials */
- FILE *fp = fopen(argv[1], "rb");
- if (!fp) {
- printf("Cannot open file %s\n", argv[2]);
- return 1;
- }
- int num_materials;
- int num_indices;
- fread(&num_materials, sizeof(int), 1, fp);
- fread(&num_indices, sizeof(int), 1, fp);
- Material *materials = (Material *)malloc(sizeof(Material) * num_materials);
- int *material_indices = (int *)malloc(sizeof(int) * num_indices);
- fread(material_indices, sizeof(int), num_indices, fp);
- for (int i = 0; i < num_materials; ++i) {
- fread(&materials[i].name_size, sizeof(int), 1, fp);
- materials[i].name = (char *)malloc(sizeof(char) * materials[i].name_size);
- fread(materials[i].name, sizeof(char), materials[i].name_size, fp);
- fread(&materials[i].a, sizeof(float), 1, fp);
- fread(&materials[i].b, sizeof(float), 1, fp);
- fread(&materials[i].c, sizeof(float), 1, fp);
- fread(&materials[i].d, sizeof(float), 1, fp);
- fread(&materials[i].s, sizeof(float), 1, fp);
- fread(&materials[i].alpha, sizeof(uint8_t), 1, fp);
- }
- fclose(fp);
-
- /* meshes */
- char line[256];
- unsigned int num_vertices;
- unsigned int num_faces;
- unsigned int num_meshes = argc - 2;
- Mesh *meshes = (Mesh *)malloc(sizeof(Mesh) * num_meshes);
- for (int i = 2; i < argc; ++i) {
- Mesh *mesh = &meshes[i - 2];
- if ((fp = fopen(argv[i], "rb")) == NULL) {
- printf("Cannot open file %s\n", argv[1]);
- return 1;
- }
- fseek(fp, 36, SEEK_SET);
- fgets(line, 32, fp);
- sscanf(line, "element vertex %d\n", &num_vertices);
- fseek(fp, 85, SEEK_CUR);
- fgets(line, 32, fp);
- sscanf(line, "element face %d\n", &num_faces);
- fseek(fp, 39, SEEK_CUR);
- fgets(line, 12, fp);
- sscanf(line, "end_header\n");
- mesh->vertices = (float *)malloc(12 * num_vertices);
- mesh->faces = (int *)malloc(sizeof(int) * 3 * num_faces);
- mesh->num_vertices = num_vertices;
- mesh->num_faces = num_faces;
- for (int i = 0; i < num_vertices; ++i) {
- fread(&mesh->vertices[3 * i], 4, 3, fp);
- fseek(fp, 8, SEEK_CUR);
- }
- for (int i = 0; i < num_faces; ++i) {
- fseek(fp, 1, SEEK_CUR);
- fread(&mesh->faces[3 * i], 4, 3, fp);
- }
- fclose(fp);
- }
-
- /* unite meshes */
- num_vertices = 0;
- num_faces = 0;
- for (int i = 0; i < num_meshes; ++i) {
- num_vertices += meshes[i].num_vertices;
- num_faces += meshes[i].num_faces;
- }
- unsigned int sum = meshes[0].num_vertices;
- for (int i = 1; i < num_meshes; ++i) {
- for (unsigned int j = 0; j != meshes[i].num_faces; ++j)
- for (int k = 0; k < 3; ++k)
- meshes[i].faces[3 * j + k] += sum;
- sum += meshes[i].num_vertices;
- }
-
- fp = fopen("mesh.ply", "wb");
- fprintf(fp, "ply\n");
- fprintf(fp, "format binary_little_endian 1.0\n");
- fprintf(fp, "element vertex %d\n", num_vertices);
- fprintf(fp, "property float x\n");
- fprintf(fp, "property float y\n");
- fprintf(fp, "property float z\n");
- fprintf(fp, "element radio_material %d\n", num_materials);
- fprintf(fp, "property list uint char name\n");
- fprintf(fp, "property float a\n");
- fprintf(fp, "property float b\n");
- fprintf(fp, "property float c\n");
- fprintf(fp, "property float d\n");
- fprintf(fp, "property float s\n");
- fprintf(fp, "property char alpha\n");
- fprintf(fp, "element face %d\n", num_faces);
- fprintf(fp, "property list uchar uint vertex_index\n");
- fprintf(fp, "property uint material_index\n");
- fprintf(fp, "end_header\n");
- for (unsigned int j = 0; j != num_meshes; ++j)
- fwrite(meshes[j].vertices, 4, 3 * meshes[j].num_vertices, fp);
- for (int i = 0; i != num_materials; ++i) {
- fwrite(&materials[i].name_size, sizeof(int), 1, fp);
- fwrite(materials[i].name, sizeof(char), materials[i].name_size, fp);
- fwrite(&materials[i].a, sizeof(float), 1, fp);
- fwrite(&materials[i].b, sizeof(float), 1, fp);
- fwrite(&materials[i].c, sizeof(float), 1, fp);
- fwrite(&materials[i].d, sizeof(float), 1, fp);
- fwrite(&materials[i].s, sizeof(float), 1, fp);
- fwrite(&materials[i].alpha, sizeof(uint8_t), 1, fp);
- }
- char n = 3;
- for (unsigned int j = 0; j != num_meshes; ++j)
- for (unsigned int k = 0; k != meshes[j].num_faces; ++k) {
- fwrite(&n, 1, 1, fp);
- fwrite(&meshes[j].faces[3 * k], 4, 3, fp);
- fwrite(&material_indices[j], 4, 1, fp);
- }
- fclose(fp);
-}
diff --git a/script/material.h b/script/material.h
@@ -1,9 +0,0 @@
-#include <stdint.h>
-
-typedef struct {
- int name_size;
- char *name;
- float a, b, c, d;
- float s;
- uint8_t alpha;
-} Material;
diff --git a/test.py b/test.py
@@ -3,7 +3,7 @@ import matplotlib.pyplot as plt
from rt import compute_paths
# Define inputs
-mesh_filepath = __file__[:__file__.rfind('/') + 1] + 'scenes/box.ply'
+mesh_filepath = __file__[:__file__.rfind('/') + 1] + 'scenes/box.hrt'
rx_positions = np.array([[0., 0., 2.5]], dtype=np.float64)
tx_positions = np.array([[0., 0., 2.51]], dtype=np.float64)
rx_velocities = np.array([[0., 0., 0.]], dtype=np.float64)
diff --git a/script/viz_ray_bin.py b/viz_ray_bin.py
