hermespy-rt

Minimalistic signal processing ray-tracer in C
git clone https://git.ea.contact/hermespy-rt
Log | Files | Refs

commit 379a405fb71653adf06384ff3aee8f337651b5af
parent dfbae03a54332ee62fa09d1c577559312234e476
Author: Egor Achkasov <eaachkasov@edu.hse.ru>
Date:   Wed, 12 Mar 2025 00:31:57 +0100

Refactor return args; Add mesh velocity; Mv *.h to inc/; Impl scene config

Diffstat:
Mcompute_paths.c | 101+++++++++++++++++++++++++++++++++++++++++--------------------------------------
Dcompute_paths.h | 74--------------------------------------------------------------------------
Mcompute_paths_pybind11.cpp | 156+++++++++++++++++++++++++++++++++++++++++++++++---------------------------------
Ainc/compute_paths.h | 60++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
Ainc/scene.h | 173+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
Rtest.h -> inc/test.h | 0
Dscene.h | 170-------------------------------------------------------------------------------
Mscene_fromSionna.c | 132++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++---
Mtest.c | 53++++++++++++++++++++++++-----------------------------
Mtest.py | 63++++++++++++++++++++++++++++++---------------------------------
10 files changed, 560 insertions(+), 422 deletions(-)

diff --git a/compute_paths.c b/compute_paths.c @@ -142,6 +142,7 @@ typedef struct { /* Normals. Size [num_triangles] */ Vec3 *ns; uint32_t material_index; + float velocity[3]; } Mesh; typedef struct { @@ -266,6 +267,8 @@ Scene load_scene(const char *scene_filepath, float carrier_frequency) exit(8); /* material_index */ if (fread(&mesh->material_index, sizeof(uint32_t), 1, f) != 1) exit(8); + /* velocity */ + if (fread(&mesh->velocity, sizeof(float), 3, f) != 3) exit(8); } fclose(f); @@ -515,21 +518,8 @@ 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 */ - /* 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_rx_scat, /* output array of directions for rx (num_rx, num_tx, num_bounces * num_paths, 3) */ - OUT float *directions_tx_scat, /* output array of directions for tx (num_paths, 3) */ - OUT float *a_te_re_scat, /* output array real parts of TE gains (num_rx, num_tx, num_bounces * num_paths) */ - OUT float *a_te_im_scat, /* output array imaginary parts of TE gains (num_rx, num_tx, num_bounces * num_paths) */ - OUT float *a_tm_re_scat, /* output array real parts of TM gains (num_rx, num_tx, num_bounces * num_paths) */ - OUT float *a_tm_im_scat, /* output array imaginary parts of TM gains (num_rx, num_tx, num_bounces * num_paths) */ - OUT float *tau_scat /* output array of delays (num_rx, num_tx, num_bounces * num_paths) */ + OUT PathsInfo *los, /* output LoS information */ + OUT PathsInfo *scatter /* output scatter information */ ) { /* Load the scene */ @@ -555,7 +545,15 @@ void compute_paths( /* Record the tx directions */ /* This is only valid for the fibbonaci sphere */ - directions_tx_scat = (float*)ray_directions; + scatter->directions_tx = (float*)memcpy( + scatter->directions_tx, + ray_directions, + num_paths * sizeof(Vec3) + ); + if (scatter->directions_tx == NULL) { + fprintf(stderr, "Failed to copy memory to scatter->directions_tx\n"); + exit(70); + } /* Cast the positions and velocities to Vec3 */ Vec3 *rx_pos_v = (Vec3*)rx_pos; @@ -612,7 +610,7 @@ void compute_paths( /* Consider imaginary parts of the LoS gains to be 0 in any way */ for (size_t off = 0; off != num_rx * num_tx; ++off) - a_te_im_los[off] = a_tm_im_los[off] = 0.f; + los->a_te_im[off] = los->a_tm_im[off] = 0.f; /* Calculate LoS paths directly from tx to rx */ for (size_t rx = 0, off = 0; rx != num_rx; ++rx) @@ -627,13 +625,14 @@ void compute_paths( /* 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; + los->directions_rx[off * 3] = 1.f; + los->directions_tx[off * 3] = -1.f; + los->directions_rx[off * 3 + 1] = los->directions_tx[off * 3 + 1] = 0.f; + los->directions_rx[off * 3 + 2] = los->directions_tx[off * 3 + 2] = 0.f; /* Set gains to 1+0j */ - a_te_re_los[off] = a_tm_re_los[off] = 1.f; + los->a_te_re[off] = los->a_tm_re[off] = 1.f; /* Set delay to 0 */ - tau_los[off] = 0.f; + los->tau[off] = 0.f; continue; } @@ -641,7 +640,7 @@ void compute_paths( moeller_trumbore(&r, &scene, &t, &mesh_ind, &face_ind, &theta); if (mesh_ind != (uint32_t)-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; + los->a_te_re[off] = los->a_tm_re[off] = los->tau[off] = 0.f; continue; } @@ -649,18 +648,21 @@ void compute_paths( /* 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; + los->directions_tx[off * 3] = r.d.x / t; + los->directions_tx[off * 3 + 1] = r.d.y / t; + los->directions_tx[off * 3 + 2] = r.d.z / t; + los->directions_rx[off * 3] = -los->directions_tx[off * 3]; + los->directions_rx[off * 3 + 1] = -los->directions_tx[off * 3 + 1]; + los->directions_rx[off * 3 + 2] = -los->directions_tx[off * 3 + 2]; /* Calculate the free space loss */ free_space_loss = free_space_loss_multiplier * t; if (free_space_loss > 1.f) { - a_te_re_los[off] = 1.f / free_space_loss; - a_tm_re_los[off] = 1.f / free_space_loss; + los->a_te_re[off] = 1.f / free_space_loss; + los->a_tm_re[off] = 1.f / free_space_loss; } else - a_te_re_los[off] = a_tm_re_los[off] = 1.f; + los->a_te_re[off] =los-> a_tm_re[off] = 1.f; /* Calculate the delay */ - tau_los[off] = t / SPEED_OF_LIGHT; + los->tau[off] = t / SPEED_OF_LIGHT; } /***************************************************************************/ @@ -775,8 +777,11 @@ void compute_paths( moeller_trumbore(&r_scat, &scene, &t, &mesh_ind, &face_ind, &theta); if (mesh_ind != (uint32_t)-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; + scatter->a_te_re[off_scat] = scatter->a_te_im[off_scat] + = scatter->a_tm_re[off_scat] + = scatter->a_tm_im[off_scat] + = scatter->tau[off_scat] + = 0.f; continue; } /* No obstacle has been hit */ @@ -786,28 +791,28 @@ void compute_paths( float a_te_re_new, a_te_im_new, a_tm_re_new, a_tm_im_new; scat_coefs(theta_s, theta, material_index, &a_te_re_new, &a_te_im_new, &a_tm_re_new, &a_tm_im_new); - a_te_re_scat[off_scat] = a_te_re_refl[off_tx_path] * a_te_re_new - - a_te_im_refl[off_tx_path] * a_te_im_new; - a_te_im_scat[off_scat] = a_te_re_refl[off_tx_path] * a_te_im_new - + a_te_im_refl[off_tx_path] * a_te_re_new; - a_tm_re_scat[off_scat] = a_tm_re_refl[off_tx_path] * a_tm_re_new - - a_tm_im_refl[off_tx_path] * a_tm_im_new; - a_tm_im_scat[off_scat] = a_tm_re_refl[off_tx_path] * a_tm_im_new - + a_tm_im_refl[off_tx_path] * a_tm_re_new; + scatter->a_te_re[off_scat] = a_te_re_refl[off_tx_path] * a_te_re_new + - a_te_im_refl[off_tx_path] * a_te_im_new; + scatter->a_te_im[off_scat] = a_te_re_refl[off_tx_path] * a_te_im_new + + a_te_im_refl[off_tx_path] * a_te_re_new; + scatter->a_tm_re[off_scat] = a_tm_re_refl[off_tx_path] * a_tm_re_new + - a_tm_im_refl[off_tx_path] * a_tm_im_new; + scatter->a_tm_im[off_scat] = a_tm_re_refl[off_tx_path] * a_tm_im_new + + a_tm_im_refl[off_tx_path] * a_tm_re_new; /* Calculate the direction */ - directions_rx_scat[off_scat * 3] = -r_scat.d.x; - directions_rx_scat[off_scat * 3 + 1] = -r_scat.d.y; - directions_rx_scat[off_scat * 3 + 2] = -r_scat.d.z; + scatter->directions_rx[off_scat * 3] = -r_scat.d.x; + scatter->directions_rx[off_scat * 3 + 1] = -r_scat.d.y; + scatter->directions_rx[off_scat * 3 + 2] = -r_scat.d.z; /* Calculate the delay */ - tau_scat[off_scat] = tau_t[off_tx_path] + dist2rx / SPEED_OF_LIGHT; + scatter->tau[off_scat] = tau_t[off_tx_path] + dist2rx / SPEED_OF_LIGHT; /* Calculate the free space loss */ free_space_loss = free_space_loss_multiplier * dist2rx; free_space_loss *= free_space_loss; if (free_space_loss > 1.f) { - a_te_re_scat[off_scat] /= free_space_loss; - a_te_im_scat[off_scat] /= free_space_loss; - a_tm_re_scat[off_scat] /= free_space_loss; - a_tm_im_scat[off_scat] /= free_space_loss; + scatter->a_te_re[off_scat] /= free_space_loss; + scatter->a_te_im[off_scat] /= free_space_loss; + scatter->a_tm_re[off_scat] /= free_space_loss; + scatter->a_tm_im[off_scat] /= free_space_loss; } } diff --git a/compute_paths.h b/compute_paths.h @@ -1,74 +0,0 @@ -#ifndef COMPUTE_PATHS_H -#define COMPUTE_PATHS_H - -#include <stddef.h> /* for size_t */ - -#define IN -#define OUT - -/** Compute gains and delays between tx and rx in a 3D scene. - * - * 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 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) - * \param tx_vel transmitter velocities, shape (num_tx, 3) - * \param carrier_frequency carrier frequency in GHz. Must be > 0.0 - * \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 num_bounces number of bounces to compute. Must be > 0 - * - * Outputs: - * - * 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) - * \param a_tm_im_los output array of imaginary parts of transverse magnetic gains for LoS, shape (num_rx, num_tx) - * \param tau_los output array of delays for LoS in seconds, shape (num_rx, num_tx) - * - * Scatter: - * \param directions_rx_scat output array of directions of incidence for rx for scatter, shape (num_rx, num_tx, num_bounces * num_paths, 3) - * \param directions_tx_scat output array of directions of incidence for tx for scatter, shape (num_paths, 3) - * \param a_te_re_scat output array of real parts of transverse electric gains for scatter, shape (num_rx, num_tx, num_bounces * num_paths) - * \param a_te_im_scat output array of imaginary parts of transverse electric gains for scatter, shape (num_rx, num_tx, num_bounces * num_paths) - * \param a_tm_re_scat output array of real parts of transverse magnetic gains for scatter, shape (num_rx, num_tx, num_bounces * num_paths) - * \param a_tm_im_scat output array of imaginary parts of transverse magnetic gains for scatter, shape (num_rx, num_tx, num_bounces * num_paths) - * \param tau_scat output array of delays for scatter in seconds, shape (num_rx, num_tx, num_bounces * num_paths) -*/ -void compute_paths( - 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) */ - IN const float *tx_vel, /* shape (num_tx, 3) */ - IN float carrier_frequency, /* > 0.0 (IN GHz!) */ - 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 */ - /* LoS */ - OUT float *directions_los, /* output array of directions of incidence for rx (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_rx_scat, /* output array of directions of incidence for rx (num_rx, num_tx, num_bounces * num_paths, 3) */ - OUT float *directions_tx_scat, /* output array of directions of incidence for tx (num_paths, 3) */ - OUT float *a_te_re_scat, /* output array real parts of TE gains (num_rx, num_tx, num_ounces * num_paths) */ - OUT float *a_te_im_scat, /* output array imaginary parts of TE gains (num_rx, num_tx, num_ounces * num_paths) */ - OUT float *a_tm_re_scat, /* output array real parts of TM gains (num_rx, num_tx, num_ounces * num_paths) */ - OUT float *a_tm_im_scat, /* output array imaginary parts of TM gains (num_rx, num_tx, num_ounces * num_paths) */ - OUT float *tau_scat /* output array of delays (num_rx, num_tx, num_ounces * num_paths) */ -); - -#endif /* COMPUTE_PATHS_H */ diff --git a/compute_paths_pybind11.cpp b/compute_paths_pybind11.cpp @@ -1,6 +1,10 @@ #include <pybind11/pybind11.h> #include <pybind11/numpy.h> +#include <complex> +#include <vector> +#include <iostream> + #ifdef _WIN32 extern "C" { #include "compute_paths.h" @@ -11,23 +15,62 @@ 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); +// Helper function to create a complex numpy array from two C arrays +py::array_t<std::complex<float>> make_complex_array( + const float* real, + const float* imag, + size_t num_rx, + size_t num_tx, + size_t num_paths +) { + size_t total_size = num_rx * num_tx * num_paths; + std::complex<float>* complex_data = new std::complex<float>[total_size]; + + for (size_t i = 0; i < total_size; ++i) + complex_data[i] = std::complex<float>(real[i], imag[i]); + + auto capsule = py::capsule(complex_data, [](void* ptr) { + delete[] static_cast<std::complex<float>*>(ptr); + }); + + return py::array_t<std::complex<float>>( + {num_rx, num_tx, num_paths}, complex_data, capsule + ); } -std::tuple< - 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_rx_scat - py::array_t<float>, // directions_tx_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 -> +class PathsInfoPython { +public: + size_t num_paths; + py::array_t<float> directions_rx; + py::array_t<float> directions_tx; + py::array_t<std::complex<float>> a_te; + py::array_t<std::complex<float>> a_tm; + py::array_t<float> tau; + + PathsInfoPython(const PathsInfo& paths, size_t num_rx, size_t num_tx) { + num_paths = paths.num_paths; + + auto capsule_deleter = [](void* ptr) { delete[] static_cast<float*>(ptr); }; + + directions_rx = py::array_t<float>( + std::vector<size_t>{num_rx, num_tx, (size_t)paths.num_paths, 3}, + paths.directions_rx, + py::capsule(paths.directions_rx, capsule_deleter) + ); + directions_tx = py::array_t<float>( + std::vector<size_t>{num_rx, num_tx, (size_t)paths.num_paths, 3}, + paths.directions_tx, + py::capsule(paths.directions_tx, capsule_deleter) + ); + + a_te = make_complex_array(paths.a_te_re, paths.a_te_im, num_rx, num_tx, paths.num_paths); + a_tm = make_complex_array(paths.a_tm_re, paths.a_tm_im, num_rx, num_tx, paths.num_paths); + + tau = py::array_t<float>({num_rx, num_tx, (size_t)paths.num_paths}, paths.tau, py::capsule(paths.tau, capsule_deleter)); + } +}; + +std::tuple<PathsInfoPython, PathsInfoPython> compute_paths_wrapper( const std::string &mesh_filepath, py::array_t<float> rx_positions, @@ -47,19 +90,26 @@ compute_paths_wrapper( py::buffer_info tx_vel_info = tx_velocities.request(); // Output - 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_rx_scat = new float[num_rx * num_tx * num_bounces * num_paths * 3]; - float *directions_tx_scat = new float[num_paths * 3]; - float *a_te_re_scat = new float[num_rx * num_tx * num_bounces * num_paths]; - float *a_te_im_scat = new float[num_rx * num_tx * num_bounces * num_paths]; - float *a_tm_re_scat = new float[num_rx * num_tx * num_bounces * num_paths]; - float *a_tm_im_scat = new float[num_rx * num_tx * num_bounces * num_paths]; - float *tau_scat = new float[num_rx * num_tx * num_bounces * num_paths]; + PathsInfo los = { + .num_paths = 1, + .directions_rx = new float[num_rx * num_tx * 3], + .directions_tx = new float[num_rx * num_tx * 3], + .a_te_re = new float[num_rx * num_tx], + .a_te_im = new float[num_rx * num_tx], + .a_tm_re = new float[num_rx * num_tx], + .a_tm_im = new float[num_rx * num_tx], + .tau = new float[num_rx * num_tx] + }; + PathsInfo scatter = { + .num_paths = (uint32_t)(num_bounces * num_paths), + .directions_rx = new float[num_rx * num_tx * num_bounces * num_paths * 3], + .directions_tx = new float[num_rx * num_tx * num_bounces * num_paths * 3], + .a_te_re = new float[num_rx * num_tx * num_bounces * num_paths], + .a_te_im = new float[num_rx * num_tx * num_bounces * num_paths], + .a_tm_re = new float[num_rx * num_tx * num_bounces * num_paths], + .a_tm_im = new float[num_rx * num_tx * num_bounces * num_paths], + .tau = new float[num_rx * num_tx * num_bounces * num_paths] + }; // Call the C function compute_paths( @@ -73,44 +123,28 @@ compute_paths_wrapper( (size_t)num_tx, (size_t)num_paths, (size_t)num_bounces, - // 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_rx_scat, - directions_tx_scat, - a_te_re_scat, - a_te_im_scat, - a_tm_re_scat, - a_tm_im_scat, - tau_scat + &los, + &scatter ); - // Cast to numpy arrays and return + // Wrap the results into Python objects return std::make_tuple( - 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_bounces * num_paths, 3}, directions_rx_scat), - make_py_array({num_paths, 3}, directions_tx_scat), - make_py_array({num_rx, num_tx, num_bounces * num_paths}, a_te_re_scat), - make_py_array({num_rx, num_tx, num_bounces * num_paths}, a_te_im_scat), - make_py_array({num_rx, num_tx, num_bounces * num_paths}, a_tm_re_scat), - make_py_array({num_rx, num_tx, num_bounces * num_paths}, a_tm_im_scat), - make_py_array({num_rx, num_tx, num_bounces * num_paths}, tau_scat) + PathsInfoPython(los, num_rx, num_tx), + PathsInfoPython(scatter, num_rx, num_tx) ); } PYBIND11_MODULE(rt, m) { + py::class_<PathsInfoPython>(m, "PathsInfo") + .def_readonly("num_paths", &PathsInfoPython::num_paths) + .def_readonly("directions_rx", &PathsInfoPython::directions_rx) + .def_readonly("directions_tx", &PathsInfoPython::directions_tx) + .def_readonly("a_te", &PathsInfoPython::a_te) + .def_readonly("a_tm", &PathsInfoPython::a_tm) + .def_readonly("tau", &PathsInfoPython::tau); + // TODO write a proper docstring - m.def("compute_paths", &compute_paths_wrapper, "Compute gains and delays in PLY scene", + m.def("compute_paths", &compute_paths_wrapper, "Compute gains and delays", py::arg("mesh_filepath"), py::arg("rx_positions"), py::arg("tx_positions"), @@ -121,10 +155,4 @@ PYBIND11_MODULE(rt, m) { py::arg("num_tx"), py::arg("num_paths"), py::arg("num_bounces")); - - // Scene filepaths - // TODO - m.def("get_scene_fp_box", - []() { return __FILE__ + std::string("/scenes/box.ply"); } - ); } diff --git a/inc/compute_paths.h b/inc/compute_paths.h @@ -0,0 +1,60 @@ +#ifndef COMPUTE_PATHS_H +#define COMPUTE_PATHS_H + +#include <stddef.h> /* for size_t */ +#include <stdint.h> /* for uint32_t */ + +#define IN +#define OUT + +/** Raytracing information returned by compute_paths for each path type (LoS, scatter) */ +typedef struct { + uint32_t num_paths; + float *directions_rx; /* shape (num_rx, num_tx, num_paths, 3) */ + float *directions_tx; /* shape (num_rx, num_tx, num_paths, 3) */ + float *a_te_re; /* shape (num_rx, num_tx, num_paths) */ + float *a_te_im; /* shape (num_rx, num_tx, num_paths) */ + float *a_tm_re; /* shape (num_rx, num_tx, num_paths) */ + float *a_tm_im; /* shape (num_rx, num_tx, num_paths) */ + float *tau; /* shape (num_rx, num_tx, num_paths) */ +} PathsInfo; + +/** Compute gains and delays between tx and rx in a 3D scene. + * + * 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. + * + * The output PathInfo structure is defined in compute_paths.h + * + * \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) + * \param tx_vel transmitter velocities, shape (num_tx, 3) + * \param carrier_frequency carrier frequency in GHz. Must be > 0.0 + * \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 num_bounces number of bounces to compute. Must be > 0 + * + * \param los output Line-of-Sight results. los->num_paths = 1 + * \param scatter output scatter results. scatter->num_paths = num_bounces * num_paths +*/ +void compute_paths( + 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) */ + IN const float *tx_vel, /* shape (num_tx, 3) */ + IN float carrier_frequency, /* > 0.0 (IN GHz!) */ + 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 PathsInfo *los, /* output LoS information */ + OUT PathsInfo *scatter /* output scatter information */ +); + +#endif /* COMPUTE_PATHS_H */ diff --git a/inc/scene.h b/inc/scene.h @@ -0,0 +1,173 @@ +#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; + /* Index of the material in g_hrt_materials array defined in scene.h */ + uint32_t material_index; + /* Global cartesian velocity vector. Size [3] */ + float velocity[3]; +} 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. */ + const 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/test.h b/inc/test.h diff --git a/scene.h b/scene.h @@ -1,170 +0,0 @@ -#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. */ - const 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 @@ -44,6 +44,7 @@ HRT_Mesh mesh_box = { .num_triangles = 12, .is = (uint32_t*)mesh_box_is, .material_index = MATERIAL_CONCRETE, + .velocity = {0.f, 0.f, 0.f}, }; HRT_Scene scene_box = { .num_meshes = 1, @@ -68,6 +69,7 @@ HRT_Mesh mesh_simpleReflector = { .num_triangles = 2, .is = (uint32_t*)mesh_simpleReflector_is, .material_index = MATERIAL_CONCRETE, + .velocity = {0.f, 0.f, 0.f}, }; HRT_Scene scene_simpleReflector = { .num_meshes = 1, @@ -110,7 +112,8 @@ HRT_Mesh readPly(const char* filepath) { .vs = NULL, .num_triangles = 0, .is = NULL, - .material_index = 0 + .material_index = 0, + .velocity = {0.f, 0.f, 0.f}, }; /* read header */ @@ -184,12 +187,102 @@ free_and_error: exit(8); } +/** Read a scene config CSV file. + * + * The CSV file must be a text file with the following format: + * material_index,velocity_x,velocity_y,velocity_z + * + * The file can have any number of lines except the first one, which is the header. + * If a mesh is not in the CSV file, + * it will have the default (0) material index and velocity. + * + * \param filepath Path to the CSV file + * \param num_meshes Output number of meshes in the CSV file + * \param mesh_filenames Output array of mesh filenames + * \param material_indicies Output array of material indicies + * \param velocity Output array of velocities. Size [num_meshes * 3] + */ +void readCsv( + const char* filepath, + uint32_t* num_meshes, + char** mesh_filenames, + uint32_t* material_indicies, + float* velocity +) +{ + FILE* f = fopen(filepath, "r"); + if (f == NULL) { + perror("Error: cannot open file"); + exit(8); + } + + /* Check header */ + char buf[256]; + if (fgets(buf, 256, f) == NULL) { + perror("Error: cannot read header"); + exit(8); + } + if (strncmp(buf, "material_index,velocity_x,velocity_y,velocity_z\n", 48)) { + perror("Error: invalid header"); + exit(8); + } + + /* Count number of lines */ + *num_meshes = 0; + while (fgets(buf, 256, f)) + ++(*num_meshes); + /* Return to the beginning of the file */ + rewind(f); + /* Skip header */ + if (fgets(buf, 256, f) == NULL) { + perror("Error: cannot read header"); + exit(8); + } + + /* Allocate memory for the arrays */ + mesh_filenames = (char**)malloc(*num_meshes * sizeof(char*)); + material_indicies = (uint32_t*)malloc(*num_meshes * sizeof(uint32_t)); + velocity = (float*)malloc(*num_meshes * 3 * sizeof(float)); + + /* Read lines */ + for (uint32_t i = 0; i != *num_meshes; ++i) { + /* Read line */ + if (fgets(buf, 256, f) == NULL) { + perror("Error: cannot read line"); + exit(8); + } + /* Parse line */ + char name[50]; + int rc = sscanf( + buf, + "%49[^,],%u,%f,%f,%f\n", + name, + &material_indicies[i], + &velocity[i * 3], + &velocity[i * 3 + 1], + &velocity[i * 3 + 2] + ); + if (rc != 4) { + perror("Error: cannot parse line"); + exit(8); + } + /* Copy name */ + mesh_filenames[i] = (char*)malloc(strlen(name) + 1); + if (mesh_filenames[i] == NULL) { + perror("Error: cannot allocate memory for mesh filename"); + exit(8); + } + strcpy(mesh_filenames[i], name); + } +} + + /* Read a Sionna .xml scene. Assumes meshes are in a "meshes" directory. * * \param filepath Path to the .xml scene file * \param scene The scene to fill */ -void read_scene(const char* filepath, HRT_Scene* scene) { +void readScene(const char* filepath, HRT_Scene* scene) { /* Check if filepath ends with ".xml" */ size_t filepath_len = strlen(filepath); if (filepath_len > 4 && strcmp(filepath + filepath_len - 4, ".xml") != 0) { @@ -197,6 +290,26 @@ void read_scene(const char* filepath, HRT_Scene* scene) { exit(8); } + /* Read the scene config CSV file */ + char* csv_path = (char*)malloc(filepath_len + 1); + if (!csv_path) { + perror("Error: failed to allocate csv_path"); + exit(8); + } + strcpy(csv_path, filepath); + strcpy(csv_path + filepath_len - 4, ".csv"); + uint32_t csv_num_meshes; + char** csv_mesh_filenames; + uint32_t* csv_material_indicies; + float* csv_velocity; + readCsv( + csv_path, + &csv_num_meshes, + csv_mesh_filenames, + csv_material_indicies, + csv_velocity + ); + /* Get the meshes directory path */ const char* last_slash = strrchr(filepath, '/'); size_t base_len = (last_slash ? (last_slash - filepath + 1) : 0); @@ -264,8 +377,18 @@ void read_scene(const char* filepath, HRT_Scene* scene) { /* Read mesh */ scene->meshes[mesh_idx] = readPly(ply_path); - ++mesh_idx; + /* Get material index and velocity from the CSV file */ + for (uint32_t i = 0; i != csv_num_meshes; ++i) + if (strcmp(entry->d_name, csv_mesh_filenames[i]) == 0) { + scene->meshes[mesh_idx].material_index = csv_material_indicies[i]; + scene->meshes[mesh_idx].velocity[0] = csv_velocity[i * 3]; + scene->meshes[mesh_idx].velocity[1] = csv_velocity[i * 3 + 1]; + scene->meshes[mesh_idx].velocity[2] = csv_velocity[i * 3 + 2]; + break; + } + + ++mesh_idx; free(ply_path); } @@ -301,7 +424,7 @@ int main(int argc, char *argv[]) { scene_ptr = &scene_simpleReflector; else { scene_ptr = (HRT_Scene*)malloc(sizeof(HRT_Scene)); - read_scene(scene_filepath, scene_ptr); + readScene(scene_filepath, scene_ptr); } FILE *fp = fopen("scene.hrt", "wb"); @@ -323,6 +446,7 @@ int main(int argc, char *argv[]) { fwrite(&scene_ptr->meshes[i].num_triangles, sizeof(uint32_t), 1, fp); fwrite(scene_ptr->meshes[i].is, sizeof(uint32_t), 3 * scene_ptr->meshes[i].num_triangles, fp); fwrite(&scene_ptr->meshes[i].material_index, sizeof(uint32_t), 1, fp); + fwrite(&scene_ptr->meshes[i].velocity, sizeof(float), 3, fp); } fclose(fp); diff --git a/test.c b/test.c @@ -16,19 +16,26 @@ int main(int argc, char **argv) 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 */ - 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)); + PathsInfo los = { + .num_paths = 1, + .directions_rx = (float*)malloc(g_numRx * g_numTx * 3 * sizeof(float)), + .directions_tx = (float*)malloc(g_numRx * g_numTx * 3 * sizeof(float)), + .a_te_re = (float*)malloc(g_numRx * g_numTx * sizeof(float)), + .a_te_im = (float*)malloc(g_numRx * g_numTx * sizeof(float)), + .a_tm_re = (float*)malloc(g_numRx * g_numTx * sizeof(float)), + .a_tm_im = (float*)malloc(g_numRx * g_numTx * sizeof(float)), + .tau = (float*)malloc(g_numRx * g_numTx * sizeof(float)) + }; + PathsInfo scatter = { + .num_paths = g_numBounces * g_numPaths, + .directions_rx = (float*)malloc(g_numRx * g_numTx * g_numBounces * g_numPaths * 3 * sizeof(float)), + .directions_tx = (float*)malloc(g_numPaths * 3 * sizeof(float)), + .a_te_re = (float*)malloc(g_numRx * g_numTx * g_numBounces * g_numPaths * sizeof(float)), + .a_te_im = (float*)malloc(g_numRx * g_numTx * g_numBounces * g_numPaths * sizeof(float)), + .a_tm_re = (float*)malloc(g_numRx * g_numTx * g_numBounces * g_numPaths * sizeof(float)), + .a_tm_im = (float*)malloc(g_numRx * g_numTx * g_numBounces * g_numPaths * sizeof(float)), + .tau = (float*)malloc(g_numRx * g_numTx * g_numBounces * g_numPaths * sizeof(float)) + }; compute_paths( argv[1], rx_positions, @@ -37,9 +44,7 @@ int main(int argc, char **argv) tx_velocities, carrier_frequency, 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 + &los, &scatter ); /* Save the results */ @@ -50,17 +55,7 @@ int main(int argc, char **argv) fwrite(data, sizeof(float), size, f); \ fclose(f); - 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); + /* Use this macro to write what you would like to inspect to a binary file */ + /* Example: */ + /* WRITE_BIN("los_directions_rx.bin", los.directions_rx, g_numRx * g_numTx * 3); */ } diff --git a/test.py b/test.py @@ -1,11 +1,10 @@ import numpy as np -import matplotlib.pyplot as plt from rt import compute_paths # Define inputs -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) +mesh_filepath = __file__[:__file__.rfind('/') + 1] + 'scenes/simple_reflector.hrt' +rx_positions = np.array([[0., 0., .15]], dtype=np.float64) +tx_positions = np.array([[0., 0., .151]], 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 @@ -15,7 +14,7 @@ num_paths = 10000 num_bounces = 3 # Call compute_paths -output = compute_paths( +los, scatter = compute_paths( mesh_filepath, rx_positions, tx_positions, @@ -27,41 +26,39 @@ output = compute_paths( num_paths, num_bounces ) -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] -directions_rx_scat = output[6] -directions_tx_scat = output[7] -a_te_scat = output[8] + 1.j*output[9] -a_tm_scat = output[10] + 1.j*output[11] -tau_scat = output[12] 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)}") -print(f"a_te_los.imag min, max: {np.min(a_te_los.imag)}, {np.max(a_te_los.imag)}") -print(f"a_tm_los.real min, max: {np.min(a_tm_los.real)}, {np.max(a_tm_los.real)}") -print(f"a_tm_los.imag min, max: {np.min(a_tm_los.imag)}, {np.max(a_tm_los.imag)}") +print(f"shape(los.directions_rx): {los.directions_rx.shape}") +print(f"shape(los.directions_tx): {los.directions_tx.shape}") +print(f"shape(los.a_te): {los.a_te.shape}") +print(f"shape(los.a_tm): {los.a_tm.shape}") +print(f"Delays: {los.tau}") +print(f"los.a_te.real min, max: {np.min(los.a_te.real)}, {np.max(los.a_te.real)}") +print(f"los.a_te.imag min, max: {np.min(los.a_te.imag)}, {np.max(los.a_te.imag)}") +print(f"los.a_tm.real min, max: {np.min(los.a_tm.real)}, {np.max(los.a_tm.real)}") +print(f"los.a_tm.imag min, max: {np.min(los.a_tm.imag)}, {np.max(los.a_tm.imag)}") print("\n#####################") print("Scatter:") print("#####################") -print(f"shape(directions_rx_scat): {directions_rx_scat.shape}") -print(f"shape(directions_tx_scat): {directions_tx_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)}") -print(f"a_te_scat.imag min, max: {np.min(a_te_scat.imag)}, {np.max(a_te_scat.imag)}") -print(f"a_tm_scat.real min, max: {np.min(a_tm_scat.real)}, {np.max(a_tm_scat.real)}") -print(f"a_tm_scat.imag min, max: {np.min(a_tm_scat.imag)}, {np.max(a_tm_scat.imag)}") +print(f"shape(scatter.directions_rx): {scatter.directions_rx.shape}") +print(f"shape(scatter.directions_tx): {scatter.directions_tx.shape}") +print(f"shape(scatter.a_te): {scatter.a_te.shape}") +print(f"shape(scatter.a_tm): {scatter.a_tm.shape}") +print(f"Delays: {scatter.tau}") +print(f"scatter.a_te.real min, max: {np.min(scatter.a_te.real)}, {np.max(scatter.a_te.real)}") +print(f"scatter.a_te.imag min, max: {np.min(scatter.a_te.imag)}, {np.max(scatter.a_te.imag)}") +print(f"scatter.a_tm.real min, max: {np.min(scatter.a_tm.real)}, {np.max(scatter.a_tm.real)}") +print(f"scatter.a_tm.imag min, max: {np.min(scatter.a_tm.imag)}, {np.max(scatter.a_tm.imag)}") # Asssert shapes -assert directions_los.shape == (num_rx, num_tx, 3) -assert directions_rx_scat.shape == (num_rx, num_tx, num_bounces * 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_rx, num_tx, num_bounces * num_paths) +assert los.num_paths == 1 +assert scatter.num_paths == num_bounces * num_paths +assert los.directions_rx.shape == (num_rx, num_tx, 1, 3) +assert los.directions_tx.shape == (num_rx, num_tx, 1, 3) +assert scatter.directions_rx.shape == (num_rx, num_tx, scatter.num_paths, 3) +assert scatter.directions_tx.shape == (num_rx, num_tx, scatter.num_paths, 3) +assert los.a_te.shape == los.a_tm.shape == los.tau.shape == (num_rx, num_tx, 1) +assert scatter.a_te.shape == scatter.a_tm.shape == scatter.tau.shape == (num_rx, num_tx, scatter.num_paths)