commit cbe54e2e4f21cc541f89c2824180cee2b122c038
parent 2f0804aa2c350c0b80e73bdbf3a1701dfff878df
Author: Egor Achkasov <eaachkasov@edu.hse.ru>
Date: Tue, 4 Feb 2025 16:20:09 +0100
rm hit_points out; add directions outs; fix memory dealloc in pybind
Diffstat:
5 files changed, 129 insertions(+), 83 deletions(-)
diff --git a/compute_paths.c b/compute_paths.c
@@ -567,14 +567,15 @@ 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 *hit_points, /* output array of hit points (num_rx, num_tx, num_paths, 3) */
/* LoS */
+ OUT float *directions_los, /* output array of directions (num_rx, num_tx, 3) */
OUT float *a_te_re_los, /* output array real parts of TE gains (num_rx, num_tx) */
OUT float *a_te_im_los, /* output array imaginary parts of TE gains (num_rx, num_tx) */
OUT float *a_tm_re_los, /* output array real parts of TM gains (num_rx, num_tx) */
OUT float *a_tm_im_los, /* output array imaginary parts of TM gains (num_rx, num_tx) */
OUT float *tau_los, /* output array of delays (num_rx, num_tx) */
/* Scatter */
+ OUT float *directions_scat, /* output array of directions (num_rx, num_tx, num_paths, 3) */
OUT float *a_te_re_scat, /* output array real parts of TE gains (num_bounces, num_rx, num_tx, num_paths) */
OUT float *a_te_im_scat, /* output array imaginary parts of TE gains (num_bounces, num_rx, num_tx, num_paths) */
OUT float *a_tm_re_scat, /* output array real parts of TM gains (num_bounces, num_rx, num_tx, num_paths) */
@@ -584,7 +585,7 @@ void compute_paths(
{
/* Init loop indices and offset variables */
size_t i, j;
- size_t off, off_scat, off_hit_points;
+ size_t off, off_scat;
/* Load the scene */
Mesh mesh = load_mesh_ply(mesh_filepath, carrier_frequency);
@@ -658,19 +659,38 @@ void compute_paths(
float free_space_loss_multiplier = 4.f * PI * carrier_frequency * 1e9 / SPEED_OF_LIGHT;
float free_space_loss;
- /* Calculate LoS paths */
+ /***************************************************************************/
+ /* LoS paths */
+ /***************************************************************************/
/* Consider imaginary parts of the LoS gains to be 0 in any way */
for (off = 0; off != num_rx * num_tx; ++off)
a_te_im_los[off] = a_tm_im_los[off] = 0.f;
+
/* Calculate LoS paths directly from tx to rx */
for (i = 0; i != num_rx; ++i)
for (j = 0; j != num_tx; ++j) {
off = i * num_tx + j;
t = -1.f;
ind = -1;
+
+ /* Create a ray from the tx to the rx */
r.o = tx_pos_v[j];
r.d = vec3_sub(&rx_pos_v[i], &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) */
+ directions_los[off * 3] = 1.f;
+ directions_los[off * 3 + 1] = 0.f;
+ directions_los[off * 3 + 2] = 0.f;
+ /* Set gains to 1+0j */
+ a_te_re_los[off] = a_tm_re_los[off] = 1.f;
+ /* Set delay to 0 */
+ tau_los[off] = 0.f;
+ continue;
+ }
+
/* Is there any abstacle between the tx and the rx? */
moeller_trumbore(&r, &mesh, &t, &ind, &theta);
if (ind != -1 && t <= 1.f) {
@@ -678,9 +698,14 @@ void compute_paths(
a_te_re_los[off] = a_tm_re_los[off] = tau_los[off] = 0.f;
continue;
}
- /* No triangle has been hit, the path is LoS */
+
+ /* No obstacle has been hit, the path is LoS */
/* Calculate the distance */
t = sqrtf(vec3_dot(&r.d, &r.d));
+ /* Calculate the direction */
+ directions_los[off * 3] = -r.d.x / t;
+ directions_los[off * 3 + 1] = -r.d.y / t;
+ directions_los[off * 3 + 2] = -r.d.z / t;
/* Calculate the free space loss */
free_space_loss = free_space_loss_multiplier * t;
if (free_space_loss > 1.f) {
@@ -692,35 +717,41 @@ void compute_paths(
tau_los[off] = t / 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 refraction rays with gains as per eqs. (31c)-(31d) 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
*/
for (i = 0; i < num_bounces; ++i) {
for (off = 0; off != num_tx * num_paths; ++off) {
+ /***********************************************************************/
+ /* Reflect the rays */
+ /***********************************************************************/
/* Check if the ray is active */
if (!(active[off / 8] & (1 << (off % 8))))
continue;
/* Init */
t = -1.f;
ind = -1;
- /* Find the hit point and trinagle.
- Calculate an angle of incidence */
+ /* 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 */
active[off / 8] &= ~(1 << (off % 8));
--num_active;
continue;
}
- /* Get the normal */
+ /* Get the hit primitive's normal */
n = mesh.normals[ind];
- /* Calculate the reflection coefficients
- R_{eTE} and R_{eTM} */
+ /* Calculate the reflection coefficients R_{eTE} and R_{eTM} */
refl_coefs(&mesh.rms[mesh.rm_indices[ind]],
- theta,
- &r_te_re, &r_te_im,
- &r_tm_re, &r_tm_im);
+ 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;
if (free_space_loss > 1.f) {
@@ -729,7 +760,7 @@ void compute_paths(
r_tm_re /= free_space_loss;
r_tm_im /= free_space_loss;
}
- /* Update the gains */
+ /* Update the gains as a' = a * refl_coefs */
a_te_re_new = a_te_re_refl[off] * r_te_re - a_te_im_refl[off] * r_te_im;
a_te_im_new = a_te_re_refl[off] * r_te_im + a_te_im_refl[off] * r_te_re;
a_tm_re_new = a_tm_re_refl[off] * r_tm_re - a_tm_im_refl[off] * r_tm_im;
@@ -740,22 +771,31 @@ void compute_paths(
a_tm_im_refl[off] = a_tm_im_new;
/* Update the delay */
tau_t[off] += t / SPEED_OF_LIGHT;
+
+ /* Move and reflect the ray */
+ r = rays[off];
/* Advance the ray to the hit point */
- *h = vec3_scale(&rays[off].d, t);
- rays[off].o = vec3_add(h, &rays[off].o);
- /* Reflect the ray's direction */
- t = vec3_dot(&rays[off].d, &n);
+ *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 */
+ t = vec3_dot(&r.d, &n);
*h = vec3_scale(&n, 2.f * t);
- rays[off].d = vec3_sub(&rays[off].d, h);
- /* Advance the ray a bit to avoid self-intersection */
- *h = vec3_scale(&rays[off].d, 1e-4f);
- rays[off].o = vec3_add(&rays[off].o, h);
+ 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);
+ /* Save the reflected ray */
+ rays[off] = r;
+
+ /***********************************************************************/
+ /* Scatter the rays */
+ /***********************************************************************/
/* Scatter a path from the hit point towards the rx */
- r.o = rays[off].o;
for (j = 0; j != num_rx; ++j) {
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) {
@@ -765,11 +805,10 @@ void compute_paths(
continue;
}
/* No obstacle has been hit */
- /* Save the hit point coordinates */
- off_hit_points = off_scat * 3;
- hit_points[off_hit_points] = r.o.x;
- hit_points[off_hit_points + 1] = r.o.y;
- hit_points[off_hit_points + 2] = r.o.z;
+ /* Calculate the direction */
+ directions_scat[off_scat * 3] = -r.d.x;
+ directions_scat[off_scat * 3 + 1] = -r.d.y;
+ directions_scat[off_scat * 3 + 2] = -r.d.z;
/* Calculate the scattering angle */
theta_s = acosf(vec3_dot(&r.d, &n) / sqrtf(vec3_dot(&r.d, &r.d)));
/* Calculate the scattering coefficients */
@@ -791,6 +830,11 @@ void compute_paths(
a_tm_im_scat[off_scat] /= free_space_loss;
}
}
+
+ /***********************************************************************/
+ /* Refract the rays */
+ /***********************************************************************/
+ /* TODO */
}
}
diff --git a/compute_paths.h b/compute_paths.h
@@ -23,9 +23,8 @@
*
* Outputs:
*
- * \param hit_points global cartesean coordinates of the hit points, shape (num_bounces, num_rx, num_tx, num_paths, 3)
- *
* LoS:
+ * \param directions_los output array of directions of incidence for LoS, shape (num_rx, num_tx, 3)
* \param a_te_re_los output array of real parts of transverse electric gains for LoS, shape (num_rx, num_tx)
* \param a_te_im_los output array of imaginary parts of transverse electric gains for LoS, shape (num_rx, num_tx)
* \param a_tm_re_los output array of real parts of transverse magnetic gains for LoS, shape (num_rx, num_tx)
@@ -33,6 +32,7 @@
* \param tau_los output array of delays for LoS in seconds, shape (num_rx, num_tx)
*
* Scatter:
+ * \param directions_scat output array of directions of incidence for scatter, shape (num_rx, num_tx, num_paths, 3)
* \param a_te_re_scat output array of real parts of transverse electric gains for scatter, shape (num_bounces, num_rx, num_tx, num_paths)
* \param a_te_im_scat output array of imaginary parts of transverse electric gains for scatter, shape (num_bounces, num_rx, num_tx, num_paths)
* \param a_tm_re_scat output array of real parts of transverse magnetic gains for scatter, shape (num_bounces, num_rx, num_tx, num_paths)
@@ -50,14 +50,15 @@ 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 *hit_points, /* bouncing hit points, shape (num_bounces, num_rx, num_tx, num_paths, 3) */
/* LoS */
+ OUT float *directions_los, /* output array of directions of incidence (num_rx, num_tx, 3) */
OUT float *a_te_re_los, /* output array real parts of TE gains (num_rx, num_tx) */
OUT float *a_te_im_los, /* output array imaginary parts of TE gains (num_rx, num_tx) */
OUT float *a_tm_re_los, /* output array real parts of TM gains (num_rx, num_tx) */
OUT float *a_tm_im_los, /* output array imaginary parts of TM gains (num_rx, num_tx) */
OUT float *tau_los, /* output array of delays (num_rx, num_tx) */
/* Scatter */
+ OUT float *directions_scat, /* output array of directions of incidence (num_rx, num_tx, num_paths, 3) */
OUT float *a_te_re_scat, /* output array real parts of TE gains (num_bounces, num_rx, num_tx, num_paths) */
OUT float *a_te_im_scat, /* output array imaginary parts of TE gains (num_bounces, num_rx, num_tx, num_paths) */
OUT float *a_tm_re_scat, /* output array real parts of TM gains (num_bounces, num_rx, num_tx, num_paths) */
diff --git a/compute_paths_pybind11.cpp b/compute_paths_pybind11.cpp
@@ -11,11 +11,20 @@ extern "C" {
namespace py = pybind11;
+// Helper function to create a numpy array from a C array
+py::array_t<float> make_py_array(std::vector<size_t> shape, float *data) {
+ return py::array_t<float>(
+ //shape, data, py::capsule(data, [](float *ptr) { delete[] ptr; })
+ shape, data
+ );
+}
+
std::tuple<
- py::array_t<float>, // hit_points
+ py::array_t<float>, // directions_los
py::array_t<float>, py::array_t<float>, // a_te_re_los, a_te_im_los
py::array_t<float>, py::array_t<float>, // a_tm_re_los, a_tm_im_los
py::array_t<float>, // tau_los
+ py::array_t<float>, // directions_scat
py::array_t<float>, py::array_t<float>, // a_te_re_scat, a_te_im_scat
py::array_t<float>, py::array_t<float>, // a_tm_re_scat, a_tm_im_scat
py::array_t<float> // tau_scat
@@ -39,12 +48,13 @@ compute_paths_wrapper(
py::buffer_info tx_vel_info = tx_velocities.request();
// Output
- float *hit_points = new float[num_bounces * num_rx * num_tx * num_paths * 3];
+ float *directions_los = new float[num_rx * num_tx * 3];
float *a_te_re_los = new float[num_rx * num_tx];
float *a_te_im_los = new float[num_rx * num_tx];
float *a_tm_re_los = new float[num_rx * num_tx];
float *a_tm_im_los = new float[num_rx * num_tx];
float *tau_los = new float[num_rx * num_tx];
+ float *directions_scat = new float[num_rx * num_tx * num_paths * 3];
float *a_te_re_scat = new float[num_bounces * num_rx * num_tx * num_paths];
float *a_te_im_scat = new float[num_bounces * num_rx * num_tx * num_paths];
float *a_tm_re_scat = new float[num_bounces * num_rx * num_tx * num_paths];
@@ -63,48 +73,36 @@ compute_paths_wrapper(
(size_t)num_tx,
(size_t)num_paths,
(size_t)num_bounces,
- hit_points, // Hit points
- a_te_re_los, a_te_im_los, a_tm_re_los, a_tm_im_los, tau_los, // LoS outputs
- a_te_re_scat, a_te_im_scat, a_tm_re_scat, a_tm_im_scat, tau_scat // Scatter outputs
+ // LoS outputs
+ directions_los,
+ a_te_re_los,
+ a_te_im_los,
+ a_tm_re_los,
+ a_tm_im_los,
+ tau_los,
+ // Scatter outputs
+ directions_scat,
+ a_te_re_scat,
+ a_te_im_scat,
+ a_tm_re_scat,
+ a_tm_im_scat,
+ tau_scat
);
- // Convert output arrays into numpy arrays for easy use in Python
- // TODO prevent copying data
- // TODO construct np.complex64 array for a
- py::array_t<float> hit_points_array = py::array_t<float>({num_bounces, num_rx, num_tx, num_paths, 3}, hit_points);
- py::array_t<float> a_te_re_los_array = py::array_t<float>({num_rx, num_tx}, a_te_re_los);
- py::array_t<float> a_te_im_los_array = py::array_t<float>({num_rx, num_tx}, a_te_im_los);
- py::array_t<float> a_tm_re_los_array = py::array_t<float>({num_rx, num_tx}, a_tm_re_los);
- py::array_t<float> a_tm_im_los_array = py::array_t<float>({num_rx, num_tx}, a_tm_im_los);
- py::array_t<float> tau_los_array = py::array_t<float>({num_rx, num_tx}, tau_los);
- py::array_t<float> a_te_re_scat_array = py::array_t<float>({num_bounces, num_rx, num_tx, num_paths}, a_te_re_scat);
- py::array_t<float> a_te_im_scat_array = py::array_t<float>({num_bounces, num_rx, num_tx, num_paths}, a_te_im_scat);
- py::array_t<float> a_tm_re_scat_array = py::array_t<float>({num_bounces, num_rx, num_tx, num_paths}, a_tm_re_scat);
- py::array_t<float> a_tm_im_scat_array = py::array_t<float>({num_bounces, num_rx, num_tx, num_paths}, a_tm_im_scat);
- py::array_t<float> tau_scat_array = py::array_t<float>({num_bounces, num_rx, num_tx, num_paths}, tau_scat);
-
- // Deallocate arrays
- delete[] hit_points;
- delete[] a_te_re_los;
- delete[] a_te_im_los;
- delete[] a_tm_re_los;
- delete[] a_tm_im_los;
- delete[] tau_los;
- delete[] a_te_re_scat;
- delete[] a_te_im_scat;
- delete[] a_tm_re_scat;
- delete[] a_tm_im_scat;
- delete[] tau_scat;
-
- // Return the results as a tuple (gains, delays)
+ // Cast to numpy arrays and return
return std::make_tuple(
- hit_points_array,
- a_te_re_los_array, a_te_im_los_array,
- a_tm_re_los_array, a_tm_im_los_array,
- tau_los_array,
- a_te_re_scat_array, a_te_im_scat_array,
- a_tm_re_scat_array, a_tm_im_scat_array,
- tau_scat_array
+ make_py_array({num_rx, num_tx, 3}, directions_los),
+ make_py_array({num_rx, num_tx}, a_te_re_los),
+ make_py_array({num_rx, num_tx}, a_te_im_los),
+ make_py_array({num_rx, num_tx}, a_tm_re_los),
+ make_py_array({num_rx, num_tx}, a_tm_im_los),
+ make_py_array({num_rx, num_tx}, tau_los),
+ make_py_array({num_rx, num_tx, num_paths, 3}, directions_scat),
+ make_py_array({num_bounces, num_rx, num_tx, num_paths}, a_te_re_scat),
+ make_py_array({num_bounces, num_rx, num_tx, num_paths}, a_te_im_scat),
+ make_py_array({num_bounces, num_rx, num_tx, num_paths}, a_tm_re_scat),
+ make_py_array({num_bounces, num_rx, num_tx, num_paths}, a_tm_im_scat),
+ make_py_array({num_bounces, num_rx, num_tx, num_paths}, tau_scat)
);
}
diff --git a/test.c b/test.c
@@ -19,12 +19,13 @@ int main(int argc, char **argv)
size_t num_tx = 1;
size_t num_paths = 10000;
size_t num_bounces = 3;
- float *hit_points = (float*)malloc(num_bounces * num_rx * num_tx * num_paths * 3 * sizeof(float));
+ 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_scat = (float*)malloc(num_rx * num_tx * num_paths * 3 * sizeof(float));
float *a_te_re_scat = (float*)malloc(num_bounces * num_rx * num_tx * num_paths * sizeof(float));
float *a_te_im_scat = (float*)malloc(num_bounces * num_rx * num_tx * num_paths * sizeof(float));
float *a_tm_re_scat = (float*)malloc(num_bounces * num_rx * num_tx * num_paths * sizeof(float));
@@ -38,9 +39,8 @@ int main(int argc, char **argv)
tx_velocities,
carrier_frequency,
num_rx, num_tx, num_paths, num_bounces,
- hit_points,
- a_te_re_los, a_te_im_los, a_tm_re_los, a_tm_im_los, tau_los,
- a_te_re_scat, a_te_im_scat, a_tm_re_scat, a_tm_im_scat, tau_scat
+ directions_los, a_te_re_los, a_te_im_los, a_tm_re_los, a_tm_im_los, tau_los,
+ directions_scat, a_te_re_scat, a_te_im_scat, a_tm_re_scat, a_tm_im_scat, tau_scat
);
/* Save the results */
@@ -51,12 +51,13 @@ int main(int argc, char **argv)
fwrite(data, sizeof(float), size, f); \
fclose(f);
- WRITE_BIN("hit_points.bin", hit_points, num_bounces * num_rx * num_tx * num_paths * 3);
+ 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_scat.bin", directions_scat, num_rx * num_tx * num_paths * 3);
WRITE_BIN("a_te_re_scat.bin", a_te_re_scat, num_bounces * num_rx * num_tx * num_paths);
WRITE_BIN("a_te_im_scat.bin", a_te_im_scat, num_bounces * num_rx * num_tx * num_paths);
WRITE_BIN("a_tm_re_scat.bin", a_tm_re_scat, num_bounces * num_rx * num_tx * num_paths);
diff --git a/test.py b/test.py
@@ -5,7 +5,7 @@ from rt import compute_paths
# Define inputs
mesh_filepath = __file__[:__file__.rfind('/') + 1] + 'scenes/box.ply'
rx_positions = np.array([[0., 0., 2.5]], dtype=np.float64)
-tx_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)
tx_velocities = np.array([[0., 0., 0.]], dtype=np.float64)
carrier_frequency = 3.0
@@ -27,17 +27,19 @@ output = compute_paths(
num_paths,
num_bounces
)
-hit_points = output[0]
+directions_los = output[0]
a_te_los = output[1] + 1.j*output[2]
a_tm_los = output[3] + 1.j*output[4]
tau_los = output[5]
-a_te_scat = output[6] + 1.j*output[7]
-a_tm_scat = output[8] + 1.j*output[9]
-tau_scat = output[10]
+directions_scat = output[6]
+a_te_scat = output[7] + 1.j*output[8]
+a_tm_scat = output[9] + 1.j*output[10]
+tau_scat = output[11]
print("#####################")
print("LoS:")
print("#####################")
+print(f"shape(directions_los): {directions_los.shape}")
print(f"shape(a_te_los): {a_te_los.shape}")
print(f"Delays: {tau_los}")
print(f"a_te_los.real min, max: {np.min(a_te_los.real)}, {np.max(a_te_los.real)}")
@@ -48,6 +50,7 @@ print(f"a_tm_los.imag min, max: {np.min(a_tm_los.imag)}, {np.max(a_tm_los.imag)}
print("\n#####################")
print("Scatter:")
print("#####################")
+print(f"shape(directions_scat): {directions_scat.shape}")
print(f"shape(a_te_scat): {a_te_scat.shape}")
print(f"Delays: {tau_scat}")
print(f"a_te_scat.real min, max: {np.min(a_te_scat.real)}, {np.max(a_te_scat.real)}")
@@ -56,6 +59,5 @@ print(f"a_tm_scat.real min, max: {np.min(a_tm_scat.real)}, {np.max(a_tm_scat.rea
print(f"a_tm_scat.imag min, max: {np.min(a_tm_scat.imag)}, {np.max(a_tm_scat.imag)}")
# Asssert shapes
-assert hit_points.shape == (num_bounces, num_rx, num_tx, num_paths, 3)
assert a_te_los.shape == a_tm_los.shape == tau_los.shape == (num_rx, num_tx)
assert a_te_scat.shape == a_tm_scat.shape == tau_scat.shape == (num_bounces, num_rx, num_tx, num_paths)
