commit aee0e840afa15dd73aaf45ee355ca8bf0097ab56
parent 12b9729956fe5f071988397df27cb013202c8f76
Author: Egor Achkasov <eaachkasov@edu.hse.ru>
Date: Mon, 23 Dec 2024 17:52:22 +0100
Comment shapes fix
Diffstat:
2 files changed, 237 insertions(+), 237 deletions(-)
diff --git a/compute_paths.c b/compute_paths.c
@@ -12,53 +12,53 @@
#define EPS 1e-6
typedef struct {
- float x, y, z;
+ float x, y, z;
} Vec3;
typedef struct {
- Vec3 o, d;
+ Vec3 o, d;
} Ray;
typedef struct {
- float a, b, c, d;
+ float a, b, c, d;
} RadioMaterial;
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;
+ 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;
} Mesh;
/* ==== 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};
+ 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};
+ 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
- };
+ 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;
+ 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};
+ return (Vec3){a->x * s, a->y * s, a->z * s};
}
/* ==== MESH LOADING ==== */
@@ -89,115 +89,115 @@ Vec3 vec3_scale(const Vec3 *a, float s)
*/
Mesh load_mesh_ply(const char *mesh_filepath)
{
- FILE *f = fopen(mesh_filepath, "rb");
- if (!f) {
- fprintf(stderr, "Could not open file %s\n", mesh_filepath);
- exit(8);
- }
+ FILE *f = fopen(mesh_filepath, "rb");
+ if (!f) {
+ fprintf(stderr, "Could not open file %s\n", mesh_filepath);
+ exit(8);
+ }
- char buff[128];
- size_t num_vertices = 0;
- size_t num_faces = 0;
- size_t num_materials = 0;
+ 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);
+ /* 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 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);
+ /* 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);
- /* 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);
+ /* 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);
+ /* 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));
+ /* 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);
+ /* 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 */
- unsigned int name_size;
- for (size_t i = 0; i < num_materials; ++i) {
- /* skip name */
- if (fread(&name_size, sizeof(unsigned int), 1, f) != 1) exit(8);
- fseek(f, name_size, SEEK_CUR);
- if (fread(&mesh.rms[i], sizeof(RadioMaterial), 1, f) != 1) exit(8);
- }
+ /* RADIO_MATERIALS */
+ unsigned int name_size;
+ for (size_t i = 0; i < num_materials; ++i) {
+ /* skip name */
+ if (fread(&name_size, sizeof(unsigned int), 1, f) != 1) exit(8);
+ fseek(f, name_size, SEEK_CUR);
+ if (fread(&mesh.rms[i], sizeof(RadioMaterial), 1, f) != 1) exit(8);
+ }
- /* 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);
- }
- 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);
+ /* 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);
}
+ 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);
+ }
- fclose(f);
- return mesh;
+ fclose(f);
+ return mesh;
}
/* ==== RT ==== */
@@ -211,132 +211,132 @@ Mesh load_mesh_ply(const char *mesh_filepath)
*/
void moeller_trumbore(Ray *ray, Mesh *mesh, float *t, int32_t *ind)
{
- /* 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 (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 > -EPS && d < EPS) continue;
- s = vec3_sub(&ray->o, &v1);
- u = vec3_dot(&s, &re2_cross) / d;
- if ((u < 0. && fabs(u) > EPS)
- || (u > 1. && fabs(u - 1.) > EPS))
- continue;
- se1_cross = vec3_cross(&s, &e1);
- v = vec3_dot(&ray->d, &se1_cross) / d;
- if ((v < 0. && fabs(v) > EPS)
- || (u + v > 1. && fabs(u + v - 1.) > EPS))
- continue;
- dist_tmp = vec3_dot(&e2, &se1_cross) / d;
- if (dist_tmp > EPS && dist_tmp < dist) {
- dist = dist_tmp;
- *t = dist;
- *ind = i;
- }
+ /* 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 (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 > -EPS && d < EPS) continue;
+ s = vec3_sub(&ray->o, &v1);
+ u = vec3_dot(&s, &re2_cross) / d;
+ if ((u < 0. && fabs(u) > EPS)
+ || (u > 1. && fabs(u - 1.) > EPS))
+ continue;
+ se1_cross = vec3_cross(&s, &e1);
+ v = vec3_dot(&ray->d, &se1_cross) / d;
+ if ((v < 0. && fabs(v) > EPS)
+ || (u + v > 1. && fabs(u + v - 1.) > EPS))
+ continue;
+ dist_tmp = vec3_dot(&e2, &se1_cross) / d;
+ if (dist_tmp > EPS && dist_tmp < dist) {
+ dist = dist_tmp;
+ *t = dist;
+ *ind = i;
}
+ }
}
/* ==== MAIN FUNCTION ==== */
void compute_paths(
- IN const char *mesh_filepath, /* path to the mesh file */
- IN const float *rx_positions, /* shape (num_rx, 3) */
- IN const float *tx_positions, /* shape (num_tx, 3) */
- IN const float *rx_velocities, /* shape (num_rx, 3) */
- IN const float *tx_velocities, /* shape (num_tx, 3) */
- IN float carrier_frequency, /* > 0.0 */
- 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 float *a_re, /* output array real parts of gains (num_tx, num_paths) */
- OUT float *a_im, /* output array imaginary parts of gains (num_tx, num_paths) */
- OUT float *tau /* output array of delays (num_tx, num_paths) */
+ IN const char *mesh_filepath, /* path to the mesh file */
+ IN const float *rx_positions, /* shape (num_rx, 3) */
+ IN const float *tx_positions, /* shape (num_tx, 3) */
+ IN const float *rx_velocities, /* shape (num_rx, 3) */
+ IN const float *tx_velocities, /* shape (num_tx, 3) */
+ IN float carrier_frequency, /* > 0.0 */
+ 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 float *a_re, /* output array real parts of gains (num_rx, num_tx, num_paths) */
+ OUT float *a_im, /* output array imaginary parts of gains (num_rx, num_tx, num_paths) */
+ OUT float *tau /* output array of delays (num_rx, num_tx, num_paths) */
)
{
- /* Load the scene */
- Mesh mesh = load_mesh_ply(mesh_filepath);
+ /* Load the scene */
+ Mesh mesh = load_mesh_ply(mesh_filepath);
- /* Calculate fibonacci sphere */
- Vec3 *ray_directions = (Vec3*)malloc(num_paths * sizeof(Vec3));
- float k, phi, theta;
- for (size_t i = 0; i < num_paths; ++i) {
- k = (float)i + .5;
- phi = acos(1. - 2. * k / num_paths);
- theta = PI * (1. + sqrt(5.)) * k;
- ray_directions[i] = (Vec3){
- cos(theta) * sin(phi),
- sin(theta) * sin(phi),
- cos(phi)
- };
- }
+ /* Calculate fibonacci sphere */
+ Vec3 *ray_directions = (Vec3*)malloc(num_paths * sizeof(Vec3));
+ float k, phi, theta;
+ for (size_t i = 0; i < num_paths; ++i) {
+ k = (float)i + .5;
+ phi = acos(1. - 2. * k / num_paths);
+ theta = PI * (1. + sqrt(5.)) * k;
+ ray_directions[i] = (Vec3){
+ cos(theta) * sin(phi),
+ sin(theta) * sin(phi),
+ cos(phi)
+ };
+ }
- /* Create num_path rays for each tx */
- /* Shape (num_tx, num_paths) */
- Ray *rays = (Ray*)malloc(num_tx * num_paths * sizeof(Ray));
- Vec3 tx_position;
- for (size_t i = 0; i < num_tx; ++i) {
- tx_position = (Vec3){tx_positions[i * 3], tx_positions[i * 3 + 1], tx_positions[i * 3 + 2]};
- for (size_t j = 0; j < num_paths; ++j) {
- rays[i * num_paths + j].o = tx_position;
- rays[i * num_paths + j].d = ray_directions[j];
- }
+ /* Create num_path rays for each tx */
+ /* Shape (num_tx, num_paths) */
+ Ray *rays = (Ray*)malloc(num_tx * num_paths * sizeof(Ray));
+ Vec3 tx_position;
+ for (size_t i = 0; i < num_tx; ++i) {
+ tx_position = (Vec3){tx_positions[i * 3], tx_positions[i * 3 + 1], tx_positions[i * 3 + 2]};
+ for (size_t j = 0; j < num_paths; ++j) {
+ rays[i * num_paths + j].o = tx_position;
+ rays[i * num_paths + j].d = ray_directions[j];
}
- free(ray_directions);
+ }
+ free(ray_directions);
- /* Calculate the paths:
- - sum distances for each path into tau
- - a???
- */
- /* TODO calculate a and tau in moeller_trumbore */
- /* shape (num_tx, num_paths) */
- Vec3 *hits = (Vec3*)malloc(num_tx * num_paths * sizeof(Vec3));
- int32_t *hit_indices = (int32_t*)malloc(num_tx * num_paths * sizeof(int32_t));
- tau = (float*)memset(tau, 0, num_tx * num_paths * sizeof(float));
- float t;
- int32_t ind;
- Ray *r;
- Vec3 *h;
- size_t off;
- for (size_t i = 0; i < num_bounces; ++i) {
- off = 0;
- for (size_t j = 0; j < num_tx; ++j)
- for (size_t k = 0; k < num_paths; ++k) {
- t = -1.;
- ind = -1;
- r = &rays[off];
- moeller_trumbore(r, &mesh, &t, &ind);
- tau[off] += t;
- hit_indices[off] = ind;
- h = &hits[off];
- *h = vec3_scale(&r->d, t);
- *h = vec3_add(h, &r->o);
- r->o = *h;
- off += 1;
- }
- }
+ /* Calculate the paths:
+ - sum distances for each path into tau
+ - a???
+ */
+ /* TODO calculate a and tau in moeller_trumbore */
+ /* shape (num_tx, num_paths) */
+ Vec3 *hits = (Vec3*)malloc(num_tx * num_paths * sizeof(Vec3));
+ int32_t *hit_indices = (int32_t*)malloc(num_tx * num_paths * sizeof(int32_t));
+ tau = (float*)memset(tau, 0, num_tx * num_paths * sizeof(float));
+ float t;
+ int32_t ind;
+ Ray *r;
+ Vec3 *h;
+ size_t off;
+ for (size_t i = 0; i < num_bounces; ++i) {
+ off = 0;
+ for (size_t j = 0; j < num_tx; ++j)
+ for (size_t k = 0; k < num_paths; ++k) {
+ t = -1.;
+ ind = -1;
+ r = &rays[off];
+ moeller_trumbore(r, &mesh, &t, &ind);
+ tau[off] += t;
+ hit_indices[off] = ind;
+ h = &hits[off];
+ *h = vec3_scale(&r->d, t);
+ *h = vec3_add(h, &r->o);
+ r->o = *h;
+ off += 1;
+ }
+ }
- /* Calculate tau */
- for (size_t i = 0; i < num_paths; ++i)
- tau[i] /= SPEED_OF_LIGHT;
+ /* Calculate tau */
+ for (size_t i = 0; i < num_paths; ++i)
+ tau[i] /= SPEED_OF_LIGHT;
- /* Free */
- free(hits);
- free(hit_indices);
- free(rays);
- /* Free the mesh */
- free(mesh.vertices);
- free(mesh.rms);
- free(mesh.indices);
- free(mesh.rm_indices);
- free(mesh.normals);
+ /* Free */
+ free(hits);
+ free(hit_indices);
+ free(rays);
+ /* Free the mesh */
+ free(mesh.vertices);
+ free(mesh.rms);
+ free(mesh.indices);
+ free(mesh.rm_indices);
+ free(mesh.normals);
}
diff --git a/compute_paths.h b/compute_paths.h
@@ -22,9 +22,9 @@
* \param num_rx number of receivers. Must be > 0
* \param num_tx number of transmitters. Must be > 0
* \param num_paths number of paths to compute. Must be > 0
- * \param a_re output array of real parts of gains, shape (num_tx, num_paths)
- * \param a_im output array of imaginary parts of gains, shape (num_tx, num_paths)
- * \param tau output array of delays in seconds, shape (num_tx, num_paths)
+ * \param a_re output array of real parts of gains, shape (num_rx, num_tx, num_paths)
+ * \param a_im output array of imaginary parts of gains, shape (num_rx, num_tx, num_paths)
+ * \param tau output array of delays in seconds, shape (num_rx, num_tx, num_paths)
*/
void compute_paths(
IN const char *mesh_filepath, /* path to the mesh file */
@@ -37,9 +37,9 @@ void compute_paths(
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 float *a_re, /* output array real parts of gains (num_tx, num_paths) */
- OUT float *a_im, /* output array imaginary parts of gains (num_tx, num_paths) */
- OUT float *tau /* output array of delays (num_tx, num_paths) */
+ OUT float *a_re, /* output array real parts of gains (num_rx, num_tx, num_paths) */
+ OUT float *a_im, /* output array imaginary parts of gains (num_rx, num_tx, num_paths) */
+ OUT float *tau /* output array of delays (num_rx, num_tx, num_paths) */
);
#endif /* COMPUTE_PATHS_H */
