commit 5dfc8fb72fbb7a1af122abf85a0b146edd02c956
parent 2e7acb2ebe5382af48399387efecfb88b1d42a3a
Author: Egor Achkasov <eaachkasov@edu.hse.ru>
Date: Wed, 26 Feb 2025 20:46:35 +0100
Impl rays visualizer; Fix active mask bug; Add rays and active log
Diffstat:
| M | compute_paths.c | | | 32 | ++++++++++++++++++++++++++------ |
| M | test.c | | | 63 | ++++++++++++++++++++++++++++++--------------------------------- |
| A | test.h | | | 11 | +++++++++++ |
| D | viz_ray_bin.py | | | 17 | ----------------- |
| A | vizrays.c | | | 317 | +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ |
5 files changed, 384 insertions(+), 56 deletions(-)
diff --git a/compute_paths.c b/compute_paths.c
@@ -673,18 +673,34 @@ void compute_paths(
* - 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) {
- size_t active_byte_index = 0, active_bit_pos = 0; /* active[active_byte_index] & (1 << active_bit_pos) */
+ /* 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) {
+ for (size_t path = 0; path < num_paths; ++path, ++off_tx_path, active_bit <<= 1) {
/* Check if the ray is active */
- if (active_bit_pos == 8) {
- active_bit_pos = 0;
+ if (!active_bit) {
+ active_bit = 1;
++active_byte_index;
}
- if (!(active[active_byte_index] & (1 << active_bit_pos)))
+ if (!(active[active_byte_index] & active_bit))
continue;
/***********************************************************************/
@@ -697,7 +713,7 @@ void compute_paths(
/* 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] &= ~(1 << active_bit_pos);
+ active[active_byte_index] &= ~active_bit;
--num_active;
continue;
}
@@ -800,7 +816,11 @@ void compute_paths(
/***********************************************************************/
/* 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);
diff --git a/test.c b/test.c
@@ -1,4 +1,5 @@
#include "compute_paths.h"
+#include "test.h"
#include <stdio.h> /* for fprintf */
#include <stdlib.h> /* for malloc */
@@ -10,28 +11,24 @@ int main(int argc, char **argv)
return 1;
}
- float rx_positions[3] = {0.0, 0.0, 2.5};
- float tx_positions[3] = {0.0, 0.0, 2.5};
+ 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 */
- size_t num_rx = 1;
- size_t num_tx = 1;
- size_t num_paths = 10000;
- size_t num_bounces = 3;
- float *directions_los = (float*)malloc(num_rx * num_tx * 3 * sizeof(float));
- float *a_te_re_los = (float*)malloc(num_rx * num_tx * sizeof(float));
- float *a_te_im_los = (float*)malloc(num_rx * num_tx * sizeof(float));
- float *a_tm_re_los = (float*)malloc(num_rx * num_tx * sizeof(float));
- float *a_tm_im_los = (float*)malloc(num_rx * num_tx * sizeof(float));
- float *tau_los = (float*)malloc(num_rx * num_tx * sizeof(float));
- float *directions_rx_scat = (float*)malloc(num_rx * num_tx * num_bounces * num_paths * 3 * sizeof(float));
- float *directions_tx_scat = (float*)malloc(num_paths * 3 * sizeof(float));
- float *a_te_re_scat = (float*)malloc(num_rx * num_tx * num_bounces * num_paths * sizeof(float));
- float *a_te_im_scat = (float*)malloc(num_rx * num_tx * num_bounces * num_paths * sizeof(float));
- float *a_tm_re_scat = (float*)malloc(num_rx * num_tx * num_bounces * num_paths * sizeof(float));
- float *a_tm_im_scat = (float*)malloc(num_rx * num_tx * num_bounces * num_paths * sizeof(float));
- float *tau_scat = (float*)malloc(num_rx * num_tx * num_bounces * num_paths * sizeof(float));
+ float *directions_los = (float*)malloc(g_numRx * g_numTx * 3 * sizeof(float));
+ float *a_te_re_los = (float*)malloc(g_numRx * g_numTx * sizeof(float));
+ float *a_te_im_los = (float*)malloc(g_numRx * g_numTx * sizeof(float));
+ float *a_tm_re_los = (float*)malloc(g_numRx * g_numTx * sizeof(float));
+ float *a_tm_im_los = (float*)malloc(g_numRx * g_numTx * sizeof(float));
+ float *tau_los = (float*)malloc(g_numRx * g_numTx * sizeof(float));
+ float *directions_rx_scat = (float*)malloc(g_numRx * g_numTx * g_numBounces * g_numPaths * 3 * sizeof(float));
+ float *directions_tx_scat = (float*)malloc(g_numPaths * 3 * sizeof(float));
+ float *a_te_re_scat = (float*)malloc(g_numRx * g_numTx * g_numBounces * g_numPaths * sizeof(float));
+ float *a_te_im_scat = (float*)malloc(g_numRx * g_numTx * g_numBounces * g_numPaths * sizeof(float));
+ float *a_tm_re_scat = (float*)malloc(g_numRx * g_numTx * g_numBounces * g_numPaths * sizeof(float));
+ float *a_tm_im_scat = (float*)malloc(g_numRx * g_numTx * g_numBounces * g_numPaths * sizeof(float));
+ float *tau_scat = (float*)malloc(g_numRx * g_numTx * g_numBounces * g_numPaths * sizeof(float));
compute_paths(
argv[1],
rx_positions,
@@ -39,7 +36,7 @@ int main(int argc, char **argv)
rx_velocities,
tx_velocities,
carrier_frequency,
- num_rx, num_tx, num_paths, num_bounces,
+ g_numRx, g_numTx, g_numPaths, g_numBounces,
directions_los, a_te_re_los, a_te_im_los, a_tm_re_los, a_tm_im_los, tau_los,
directions_rx_scat, directions_rx_scat,
a_te_re_scat, a_te_im_scat, a_tm_re_scat, a_tm_im_scat, tau_scat
@@ -53,17 +50,17 @@ int main(int argc, char **argv)
fwrite(data, sizeof(float), size, f); \
fclose(f);
- WRITE_BIN("directions_los.bin", directions_los, num_rx * num_tx * 3);
- WRITE_BIN("a_te_re_los.bin", a_te_re_los, num_rx * num_tx);
- WRITE_BIN("a_te_im_los.bin", a_te_im_los, num_rx * num_tx);
- WRITE_BIN("a_tm_re_los.bin", a_tm_re_los, num_rx * num_tx);
- WRITE_BIN("a_tm_im_los.bin", a_tm_im_los, num_rx * num_tx);
- WRITE_BIN("tau_los.bin", tau_los, num_rx * num_tx);
- WRITE_BIN("directions_rx_scat.bin", directions_rx_scat, num_rx * num_tx * num_bounces * num_paths * 3);
- WRITE_BIN("directions_tx_scat.bin", directions_tx_scat, num_paths * 3);
- WRITE_BIN("a_te_re_scat.bin", a_te_re_scat, num_rx * num_tx * num_bounces * num_paths);
- WRITE_BIN("a_te_im_scat.bin", a_te_im_scat, num_rx * num_tx * num_bounces * num_paths);
- WRITE_BIN("a_tm_re_scat.bin", a_tm_re_scat, num_rx * num_tx * num_bounces * num_paths);
- WRITE_BIN("a_tm_im_scat.bin", a_tm_im_scat, num_rx * num_tx * num_bounces * num_paths);
- WRITE_BIN("tau_scat.bin", tau_scat, num_rx * num_tx * num_bounces * num_paths);
+ WRITE_BIN("directions_los.bin", directions_los, g_numRx * g_numTx * 3);
+ WRITE_BIN("a_te_re_los.bin", a_te_re_los, g_numRx * g_numTx);
+ WRITE_BIN("a_te_im_los.bin", a_te_im_los, g_numRx * g_numTx);
+ WRITE_BIN("a_tm_re_los.bin", a_tm_re_los, g_numRx * g_numTx);
+ WRITE_BIN("a_tm_im_los.bin", a_tm_im_los, g_numRx * g_numTx);
+ WRITE_BIN("tau_los.bin", tau_los, g_numRx * g_numTx);
+ WRITE_BIN("directions_rx_scat.bin", directions_rx_scat, g_numRx * g_numTx * g_numBounces * g_numPaths * 3);
+ WRITE_BIN("directions_tx_scat.bin", directions_tx_scat, g_numPaths * 3);
+ WRITE_BIN("a_te_re_scat.bin", a_te_re_scat, g_numRx * g_numTx * g_numBounces * g_numPaths);
+ WRITE_BIN("a_te_im_scat.bin", a_te_im_scat, g_numRx * g_numTx * g_numBounces * g_numPaths);
+ WRITE_BIN("a_tm_re_scat.bin", a_tm_re_scat, g_numRx * g_numTx * g_numBounces * g_numPaths);
+ WRITE_BIN("a_tm_im_scat.bin", a_tm_im_scat, g_numRx * g_numTx * g_numBounces * g_numPaths);
+ WRITE_BIN("tau_scat.bin", tau_scat, g_numRx * g_numTx * g_numBounces * g_numPaths);
}
diff --git a/test.h b/test.h
@@ -0,0 +1,11 @@
+#ifndef TEST_H
+#define TEST_H
+
+#include <stdint.h>
+
+uint32_t g_numPaths = 10000;
+uint32_t g_numBounces = 3;
+uint32_t g_numTx = 1;
+uint32_t g_numRx = 1;
+
+#endif
diff --git a/viz_ray_bin.py b/viz_ray_bin.py
@@ -1,17 +0,0 @@
-import open3d as o3d
-import numpy as np
-from sys import argv
-
-pcs = []
-for arg in argv[1:]:
- f = open(arg, 'rb')
- #a = np.fromfile(f, dtype=np.float32).reshape((-1, 3))
- #a = o3d.geometry.PointCloud(o3d.utility.Vector3dVector(a))
- #pcs.append(a)
- a = np.fromfile(f, dtype=np.float32).reshape((-1, 10000, 3))
- f.close()
- for ai in a:
- ai = o3d.geometry.PointCloud(o3d.utility.Vector3dVector(ai))
- pcs.append(ai)
-pcs.append(o3d.io.read_triangle_mesh('scenes/box.ply'))
-o3d.visualization.draw(pcs, show_ui=True)
diff --git a/vizrays.c b/vizrays.c
@@ -0,0 +1,317 @@
+#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>
+
+/**
+ * Structs
+ */
+
+typedef struct {
+ float x, y, z;
+} Vec3;
+
+typedef struct {
+ Vec3 o;
+ 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, 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));
+ fread(g_rays, sizeof(float), fileSize, file);
+ 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));
+ fread(g_active, sizeof(uint8_t), activeSize, file);
+ 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];
+ fread(magic, 1, 3, file);
+ if (magic[0] != 'H' || magic[1] != 'R' || magic[2] != 'T') {
+ printf("Invalid file format\n");
+ exit(8);
+ }
+
+ /* Scene */
+ fread(&g_scene.num_meshes, sizeof(uint32_t), 1, file);
+ 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)
+ {
+ fread(&mesh->num_vertices, sizeof(uint32_t), 1, file);
+ mesh->vs = (Vec3*)malloc(mesh->num_vertices * sizeof(Vec3));
+ fread(mesh->vs, sizeof(Vec3), mesh->num_vertices, file);
+ fread(&mesh->num_triangles, sizeof(uint32_t), 1, file);
+ mesh->is = (uint32_t*)malloc(mesh->num_triangles * 3 * sizeof(uint32_t));
+ fread(mesh->is, sizeof(uint32_t), mesh->num_triangles * 3, file);
+ }
+ 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];
+ 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("Skipped %d paths\n", num_skipped);
+
+ /* 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 display() {
+ glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
+ glLoadIdentity();
+
+ /* Apply camera transformation */
+ glRotatef(g_cam.pitch, 1.0f, 0.0f, 0.0f);
+ glRotatef(g_cam.yaw, 0.0f, 1.0f, 0.0f);
+ glTranslatef(-g_cam.o.x, -g_cam.o.y, -g_cam.o.z);
+
+ drawScene();
+ drawRays();
+
+ glutSwapBuffers();
+}
+
+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 keyboard(unsigned char key, int x, int y) {
+ float speed = 0.5f;
+ switch (key) {
+ case 'w': // Forward
+ g_cam.o.x += sin(g_cam.yaw * M_PI/180.0f) * speed;
+ g_cam.o.z -= cos(g_cam.yaw * M_PI/180.0f) * speed;
+ break;
+ case 's': // Backward
+ g_cam.o.x -= sin(g_cam.yaw * M_PI/180.0f) * speed;
+ g_cam.o.z += cos(g_cam.yaw * M_PI/180.0f) * speed;
+ break;
+ case 'a': // Left
+ g_cam.o.x -= cos(g_cam.yaw * M_PI/180.0f) * speed;
+ g_cam.o.z -= sin(g_cam.yaw * M_PI/180.0f) * speed;
+ break;
+ case 'd': // Right
+ g_cam.o.x += cos(g_cam.yaw * M_PI/180.0f) * speed;
+ g_cam.o.z += sin(g_cam.yaw * M_PI/180.0f) * speed;
+ break;
+ case 'q': // Up
+ g_cam.o.y += speed;
+ break;
+ case 'e': // Down
+ g_cam.o.y -= speed;
+ break;
+ case 'x': // Next bounce
+ if (g_bounce_cur != g_numBounces)
+ g_bounce_cur++;
+ break;
+ case 'z': // Previous bounce
+ if (g_bounce_cur != 0)
+ g_bounce_cur--;
+ break;
+ case 27: // Escape key
+ free(g_rays);
+ for (uint32_t i = 0; i < g_scene.num_meshes; i++) {
+ free(g_scene.meshes[i].vs);
+ free(g_scene.meshes[i].is);
+ }
+ free(g_scene.meshes);
+ exit(0);
+ }
+ glutPostRedisplay();
+}
+
+void motion(int x, int y) {
+ if (!g_mouseDown) return;
+ g_cam.yaw += (x - g_cam.lastX) * 0.2f;
+ g_cam.pitch += (y - g_cam.lastY) * 0.2f;
+ g_cam.pitch = fmax(-89.0f, fmin(89.0f, g_cam.pitch));
+ g_cam.lastX = x;
+ g_cam.lastY = y;
+ glutPostRedisplay();
+}
+
+void mouse(int button, int state, int x, int y) {
+ if (button == GLUT_LEFT_BUTTON) {
+ if (state == GLUT_DOWN) {
+ g_mouseDown = 1;
+ g_cam.lastX = x;
+ g_cam.lastY = y;
+ } else {
+ g_mouseDown = 0;
+ }
+ }
+}
+
+/**
+ * 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(motion);
+
+ glutMainLoop();
+
+ return 0;
+}
