commit 5b5adef14a58975360f9a790aaa78ae0090a2fab
parent b7f19dffbb1e09485c00ed0dc3f985e220d17604
Author: Egor Achkasov <eaachkasov@edu.hse.ru>
Date: Wed, 4 Dec 2024 10:43:08 +0100
Change signature; fread error handling; continous arrays...
Change ints to size_t in compute_paths signature
Call exit(8) on incorrect fread and fgets return
Make hits and hit_indices continuous
Diffstat:
2 files changed, 68 insertions(+), 60 deletions(-)
diff --git a/compute_paths.c b/compute_paths.c
@@ -6,6 +6,7 @@
#include <string.h> /* for strlen, sprintf */
#include <stdio.h> /* for fopen, FILE, fclose */
#include <math.h> /* for sin, cos, sqrt */
+#include <stdint.h>
#define PI 3.14159265358979323846
#define EPS 1e-6
@@ -101,52 +102,52 @@ Mesh load_mesh_ply(const char *mesh_filepath)
/* HEADER */
/* ply */
- fread(buff, 1, 4, f);
+ if (fread(buff, 1, 4, f) != 4) exit(8);
if (strncmp(buff, "ply\n", 4)) exit(8);
/* format binary_little_endian 1.0 */
- fread(buff, 1, 32, f);
+ 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> */
- fgets(buff, 128, f);
+ 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 */
- fread(buff, 1, 17*3, f);
+ 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 */
- fgets(buff, 128, f);
+ 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 */
- fread(buff, 1, 29, f);
+ 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 */
- fread(buff, 1, 17*4, f);
+ 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> */
- fgets(buff, 128, f);
+ 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 */
- fread(buff, 1, 38, f);
+ 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 */
- fread(buff, 1, 29, f);
+ if (fread(buff, 1, 29, f) != 29) exit(8);
if (strncmp(buff, "property uint material_index\n", 29)) exit(8);
/* end_header */
- fread(buff, 1, 11, f);
+ if (fread(buff, 1, 11, f) != 11) exit(8);
if (strncmp(buff, "end_header\n", 11)) exit(8);
/* Init mesh */
@@ -162,15 +163,16 @@ Mesh load_mesh_ply(const char *mesh_filepath)
/* VERTICES */
for (size_t i = 0; i < mesh.num_vertices; ++i)
- fread(&mesh.vertices[i], sizeof(Vec3), 1, f);
+ 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 */
- fread(&name_size, sizeof(unsigned int), 1, f);
+ if (fread(&name_size, sizeof(unsigned int), 1, f) != 1) exit(8);
fseek(f, name_size, SEEK_CUR);
- fread(&mesh.rms[i], sizeof(RadioMaterial), 1, f);
+ if (fread(&mesh.rms[i], sizeof(RadioMaterial), 1, f) != 1) exit(8);
}
/* INDICES */
@@ -178,13 +180,13 @@ Mesh load_mesh_ply(const char *mesh_filepath)
uint8_t n;
Vec3 *v1, *v2, *v3, u, v;
for (size_t i = 0; i < mesh.num_indices; i += 3) {
- fread(&n, 1, 1, f);
+ if (fread(&n, 1, 1, f) != 1) exit(8);
if (n != 3) {
fprintf(stderr, "Only trianglar faces are supported\n");
exit(8);
}
- fread(&mesh.indices[i], sizeof(int32_t), 3, f);
- fread(&mesh.rm_indices[i / 3], sizeof(int32_t), 1, f);
+ 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]];
@@ -207,7 +209,7 @@ Mesh load_mesh_ply(const char *mesh_filepath)
* \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.
*/
-void moeller_trumbore(Ray *ray, Mesh *mesh, float *t, size_t *ind)
+void moeller_trumbore(Ray *ray, Mesh *mesh, float *t, int32_t *ind)
{
/* TODO BVH */
/* for each triangle */
@@ -252,10 +254,10 @@ void compute_paths(
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 int num_rx, /* number of receivers */
- IN int num_tx, /* number of transmitters */
- IN int num_paths, /* number of paths */
- IN int num_bounces, /* number of bounces */
+ 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, /* output array of gains (num_paths,) */
OUT float *tau /* output array of delays (num_paths,) */
)
@@ -268,7 +270,7 @@ void compute_paths(
float k, phi, theta;
for (size_t i = 0; i < num_paths; ++i) {
k = (float)i + .5;
- phi = 1. - 2. * k / num_paths;
+ phi = acos(1. - 2. * k / num_paths);
theta = PI * (1. + sqrt(5.)) * k;
ray_directions[i] = (Vec3){
cos(theta) * sin(phi),
@@ -279,14 +281,13 @@ void compute_paths(
/* Create num_path rays for each tx */
/* Shape (num_tx, num_paths) */
- Ray **rays = (Ray**)malloc(num_tx * sizeof(Ray*));
+ Ray *rays = (Ray*)malloc(num_tx * num_paths * sizeof(Ray));
Vec3 tx_position;
for (size_t i = 0; i < num_tx; ++i) {
- rays[i] = (Ray*)malloc(num_paths * sizeof(Ray));
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][j].o = tx_position;
- rays[i][j].d = ray_directions[j];
+ rays[i * num_paths + j].o = tx_position;
+ rays[i * num_paths + j].d = ray_directions[j];
}
}
free(ray_directions);
@@ -294,31 +295,38 @@ void compute_paths(
/* Calculate the paths */
/* TODO calculate a and tau in moeller_trumbore */
/* shape (num_bounces, num_tx, num_paths) */
- Vec3 ***hits = (Vec3***)malloc(num_bounces * sizeof(Vec3**));
- size_t ***hit_indices = (size_t***)malloc(num_bounces * sizeof(size_t**));
+ Vec3 *hits = (Vec3*)malloc(num_bounces * num_tx * num_paths * sizeof(Vec3));
+ int32_t *hit_indices = (int32_t*)malloc(num_bounces * num_tx * num_paths * sizeof(int32_t));
float t;
- size_t ind;
- for (size_t i = 0; i < num_bounces; ++i) {
- hits[i] = (Vec3**)malloc(num_tx * sizeof(Vec3*));
- hit_indices[i] = (size_t**)malloc(num_tx * sizeof(size_t*));
- for (size_t j = 0; j < num_tx; ++j) {
- hits[i][j] = (Vec3*)malloc(num_paths * sizeof(Vec3));
- hit_indices[i][j] = (size_t*)malloc(num_paths * sizeof(size_t));
+ int32_t ind;
+ Ray *r;
+ Vec3 *h;
+ for (size_t i = 0; i < num_bounces; ++i)
+ for (size_t j = 0; j < num_tx; ++j)
for (size_t k = 0; k < num_paths; ++k) {
t = -1.;
ind = -1;
- moeller_trumbore(&rays[j][k], &mesh, &t, &ind);
- hits[i][j][k] = vec3_scale(&rays[j][k].d, t);
- hits[i][j][k] = vec3_add(&hits[i][j][k], &rays[j][k].o);
- rays[j][k].o = hits[i][j][k];
- hit_indices[i][j][k] = ind;
+ r = &rays[j * num_paths + k];
+ moeller_trumbore(r, &mesh, &t, &ind);
+ hit_indices[i * num_tx * num_paths + j * num_paths + k] = ind;
+ h = &hits[i * num_tx * num_paths + j * num_paths + k];
+ *h = vec3_scale(&r->d, t);
+ *h = vec3_add(h, &r->o);
+ r->o = *h;
}
- }
- }
+ FILE *f = fopen("hits.bin", "wb");
+ fwrite(hits, sizeof(Vec3), num_bounces * num_tx * num_paths, f);
+ fclose(f);
/* TODO Remove */
- a = hits;
- tau = hit_indices;
+ for (size_t i = 0; i < num_bounces; ++i)
+ for (size_t j = 0; j < num_rx; ++j)
+ for (size_t k = 0; k < num_paths; ++k) {
+ h = &hits[i * num_tx * num_paths + j * num_paths + k];
+ a[i * num_rx * num_paths * 3 + j * num_paths * 3 + k * 3] = h->x;
+ a[i * num_rx * num_paths * 3 + j * num_paths * 3 + k * 3 + 1] = h->y;
+ a[i * num_rx * num_paths * 3 + j * num_paths * 3 + k * 3 + 2] = h->z;
+ }
/* Free */
/* Free the hits */
@@ -335,13 +343,13 @@ void compute_paths(
free(hit_indices);
*/
/* Free the rays */
- for (size_t i = 0; i < num_tx; ++i)
- free(rays[i]);
- free(rays);
- /* Free the mesh */
- free(mesh.vertices);
- free(mesh.rms);
- free(mesh.indices);
- free(mesh.rm_indices);
- free(mesh.normals);
+ // for (size_t i = 0; i < num_tx; ++i)
+ // free(rays[i]);
+ // 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
@@ -1,14 +1,14 @@
#ifndef COMPUTE_PATHS_H
#define COMPUTE_PATHS_H
-#include <stdint.h> /* for int32_t */
+#include <stddef.h> /* for size_t */
#define IN
#define OUT
/** Compute gains and delays between in a Mitsuba scene.
*
- * Scene must be defined in a Mitsuba .xml JSON format.
+ * Scene must be defined in a specific PLY format. See README for details.
* The meshes must be present in the same directory as the .xml file,
* under meshes/ directory in .ply files in PLY format.
* Meshes filenames must be the same as the name of the object in the .xml file.
@@ -32,10 +32,10 @@ void compute_paths(
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 int num_rx, /* number of receivers */
- IN int num_tx, /* number of transmitters */
- IN int num_paths, /* number of paths */
- IN int num_bounces, /* number of bounces */
+ 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, /* output array of gains (num_paths,) */
OUT float *tau /* output array of delays (num_paths,) */
);
