Skip to content

Commit

Permalink
debug ray-grid-traversal
Browse files Browse the repository at this point in the history
  • Loading branch information
@relleums
relleums committed May 27, 2024
1 parent 16acc12 commit 45f9e97
Show file tree
Hide file tree
Showing 4 changed files with 122 additions and 54 deletions.
18 changes: 11 additions & 7 deletions libs/mli/src/mli_ray_grid_traversal.c
View file
Original file line number Diff line number Diff line change
@@ -1,5 +1,7 @@
/* Copyright 2018-2024 Sebastian Achim Mueller */
#include "mli_ray_grid_traversal.h"
#include "mli_math.h"
#include <assert.h>

/* Inspired by:
* A Fast Voxel Traversal Algorithm for Ray Tracing
Expand Down Expand Up @@ -117,17 +119,17 @@ struct mliVec mliAxisAlignedGridTraversal_first_plane(

double calc_t_for_x_plane(const double x_plane, const struct mliRay* ray)
{
return (ray->support.x + x_plane) / (ray->direction.x);
return (ray->support.x + x_plane) / (ray->direction.x);
}

double calc_t_for_y_plane(const double y_plane, const struct mliRay* ray)
{
return (ray->support.y + y_plane) / (ray->direction.y);
return (ray->support.y + y_plane) / (ray->direction.y);
}

double calc_t_for_z_plane(const double z_plane, const struct mliRay* ray)
{
return (ray->support.z + z_plane) / (ray->direction.z);
return (ray->support.z + z_plane) / (ray->direction.z);
}

struct mliAxisAlignedGridTraversal mliAxisAlignedGridTraversal_start(
Expand Down Expand Up @@ -155,9 +157,9 @@ struct mliAxisAlignedGridTraversal mliAxisAlignedGridTraversal_start(
traversal.tMax.y = calc_t_for_y_plane(first_plane.y, ray);
traversal.tMax.z = calc_t_for_z_plane(first_plane.z, ray);

traversal.tDelta.x = grid->bin_width.x / ray->direction.x;
traversal.tDelta.y = grid->bin_width.y / ray->direction.y;
traversal.tDelta.z = grid->bin_width.z / ray->direction.z;
traversal.tDelta.x = fabs(grid->bin_width.x / ray->direction.x);
traversal.tDelta.y = fabs(grid->bin_width.y / ray->direction.y);
traversal.tDelta.z = fabs(grid->bin_width.z / ray->direction.z);
}
return traversal;
}
Expand Down Expand Up @@ -208,8 +210,10 @@ void mliAxisAlignedGridTraversal_fprint(FILE* f, struct mliAxisAlignedGridTraver
{
struct mliAxisAlignedGridTraversal* t = traversal;
fprintf(f, " grid.bounds.upper: [%f, %f, %f]\n", t->grid->bounds.upper.x, t->grid->bounds.upper.y, t->grid->bounds.upper.z);
fprintf(f, " grid.bounds.lower: [%f, %f, %f]\n", t->grid->bounds.lower.x, t->grid->bounds.lower.y, t->grid->bounds.upper.z);
fprintf(f, " grid.bounds.lower: [%f, %f, %f]\n", t->grid->bounds.lower.x, t->grid->bounds.lower.y, t->grid->bounds.lower.z);
fprintf(f, " grid.num_bins: [%ld, %ld, %ld]\n", t->grid->num_bins.x, t->grid->num_bins.y, t->grid->num_bins.z);

fprintf(f, " valid: %d\n", t->valid);
fprintf(f, " voxel: [%ld, %ld, %ld]\n", t->voxel.x, t->voxel.y, t->voxel.z);
fprintf(f, " step: [%f, %f, %f]\n", t->step.x, t->step.y, t->step.z);
fprintf(f, " tMax: [%f, %f, %f]\n", t->tMax.x, t->tMax.y, t->tMax.z);
Expand Down
3 changes: 2 additions & 1 deletion libs/mli/src/mli_ray_grid_traversal.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -4,7 +4,8 @@

#include <stdint.h>
#include "mliVec.h"
#include "mliRay_AABB.h"
#include "mliAABB.h"
#include "mliRay.h"

struct mliIdx3 {
int64_t x;
Expand Down
32 changes: 32 additions & 0 deletions libs/mli/src/mli_ray_grid_traversal.test.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -216,3 +216,35 @@ CASE("ray and grid traversal simple")
CHECK(rc == traversal.valid);
CHECK(!traversal.valid);
}


CASE("Actual example from simulated shower")
{
int rc = 10;
struct mliAxisAlignedGrid grid;
struct mliAxisAlignedGridTraversal traversal;
struct mliRay ray;

grid = mliAxisAlignedGrid_set(
mliAABB_set(
mliVec_init(-5119268.215672, -5119974.912158, -5000.000000),
mliVec_init(5120731.784328, 5120025.087842, 5000.000000)
),
mliIdx3_set(1024, 1024, 1)
);

ray.support = mliVec_init(4.171968e+03, 4.857704e+03, 0.000000e+00);
ray.direction = mliVec_init(-3.894984e-02, -6.049007e-01, 7.953478e-01);

traversal = mliAxisAlignedGridTraversal_start(
&grid,
&ray
);
mliAxisAlignedGridTraversal_fprint(stderr, &traversal);

while (traversal.valid && rc > 0) {
rc -= 1;
mliAxisAlignedGridTraversal_next(&traversal);
mliAxisAlignedGridTraversal_fprint(stderr, &traversal);
}
}
123 changes: 77 additions & 46 deletions libs/mli_corsika/apps/ground_grid.main.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -9,6 +9,14 @@
#include "../../mli/src/mliIo.h"
#include "../../mli/src/mliStr_numbers.h"
#include "../../mli/src/mliDynArray.h"
#include "../../mli/src/mli_ray_grid_traversal.h"
#include <time.h>


double clock2second(const clock_t t)
{
return ((double)t) / CLOCKS_PER_SEC;
}


int read_config(
Expand Down Expand Up @@ -75,40 +83,51 @@ int main(int argc, char *argv[])
float runh[273] = {0.0};
float evth[273] = {0.0};
float raw[8] = {0.0};
double xystart = -500 * 50e2;
double xystop = 500 * 50e2;

unsigned int num_bins_each_axis = 1000;
uint64_t numover = (uint64_t)(1.5 * (double)num_bins_each_axis);

struct mliCorsikaHistogram2d hist = mliCorsikaHistogram2d_init();

struct mliCorsikaPhotonBunch bunch;
struct mliAxisAlignedGrid grid;
struct mliAxisAlignedGridTraversal traversal;

struct mliDynDouble x_bin_edges = mliDynDouble_init();
struct mliDynDouble y_bin_edges = mliDynDouble_init();
struct mliDynDouble z_bin_edges = mliDynDouble_init();

struct mliDynUint32 x_idxs = mliDynUint32_init();
struct mliDynUint32 y_idxs = mliDynUint32_init();
struct mliDynUint32 z_idxs = mliDynUint32_init();
struct mliDynDouble overlaps = mliDynDouble_init();

mli_linspace(xystart, xystop, x_bin_edges.array, x_bin_edges.size);
mli_linspace(xystart, xystop, y_bin_edges.array, y_bin_edges.size);
mli_linspace(-50e2, 50e2, z_bin_edges.array, z_bin_edges.size);

mliDynUint32_malloc(&x_idxs, numover);
mliDynUint32_malloc(&y_idxs, numover);
mliDynUint32_malloc(&z_idxs, numover);
mliDynDouble_malloc(&overlaps, numover);

chk_msg(argc == 4,
"Usage: ground_grid event_tape_path out_path config_path");

chk_msg(read_config(argv[3], &x_bin_edges, &y_bin_edges, &z_bin_edges),
"Can not read config with bin edges.");

chk_msg(z_bin_edges.size == 2, "Expected z_bin_edges.size == 2.");
chk_msg(z_bin_edges.array[0] < 0, "Expected z_bin_edges[0] < 0");
chk_msg(z_bin_edges.array[1] > 0, "Expected z_bin_edges[1] > 0");

grid = mliAxisAlignedGrid_set(
mliAABB_set(
mliVec_init(
x_bin_edges.array[0],
y_bin_edges.array[0],
z_bin_edges.array[0]
),
mliVec_init(
x_bin_edges.array[x_bin_edges.size - 1],
y_bin_edges.array[y_bin_edges.size - 1],
z_bin_edges.array[z_bin_edges.size - 1]
)
),
mliIdx3_set(
x_bin_edges.size - 1,
y_bin_edges.size - 1,
z_bin_edges.size - 1
)
);

mliDynDouble_free(&x_bin_edges);
mliDynDouble_free(&y_bin_edges);
mliDynDouble_free(&z_bin_edges);

istream = fopen(argv[1], "rb");
ostream = fopen(argv[2], "wt");
chk(mliEventTapeReader_begin(&arc, istream));
Expand All @@ -118,7 +137,8 @@ int main(int argc, char *argv[])
chk(mliIo_malloc_capacity(&buff, 10 * 1000));

while (mliEventTapeReader_read_evth(&arc, evth)) {
uint64_t bunch_index = 0;
clock_t t_ray_voxel = 0;
clock_t t_avl_histogram = 0;
uint64_t xx, yy, ii, jj = 0;

if (hist.dict.capacity > 10 * 1000 * 1000) {
Expand All @@ -135,31 +155,47 @@ int main(int argc, char *argv[])
}

while (mliEventTapeReader_read_cherenkov_bunch(&arc, raw)) {
unsigned int oo;
clock_t t1, t2, t3, t4, t5;
int num_overlaps = 0;
struct mliRay ray;
mliCorsikaPhotonBunch_set_from_raw(&bunch, raw);
mli_corsika_overlap_of_ray_with_voxels(
&bunch,
&x_bin_edges,
&y_bin_edges,
&z_bin_edges,
&x_idxs,
&y_idxs,
&z_idxs,
&overlaps);

for (oo = 0; oo < overlaps.size; oo ++)
{
ray.support.x = bunch.x_cm;
ray.support.y = bunch.y_cm;
ray.support.z = 0.0;
ray.direction.x = mli_corsika_ux_to_cx(bunch.ux);
ray.direction.y = mli_corsika_vy_to_cy(bunch.vy);
ray.direction.z = mli_corsika_restore_direction_z_component(
ray.direction.x,
ray.direction.y
);

t1 = clock();
traversal = mliAxisAlignedGridTraversal_start(
&grid,
&ray
);
/*mliAxisAlignedGridTraversal_fprint(stderr, &traversal);*/

t2 = clock();
t_ray_voxel += (t2 - t1);

while (traversal.valid) {
num_overlaps += 1;
t3 = clock();
chk(
mliCorsikaHistogram2d_assign(
&hist,
x_idxs.array[oo],
y_idxs.array[oo],
traversal.voxel.x,
traversal.voxel.y,
bunch.weight_photons
)
);
t4 = clock();
t_avl_histogram += (t4 - t3);
mliAxisAlignedGridTraversal_next(&traversal);
t5 = clock();
t_ray_voxel += (t5 - t4);
}

bunch_index += 1;
}
chk(mliCorsikaHistogram2d_dumps(&hist, &buff));
buff.pos = 0;
Expand All @@ -169,22 +205,17 @@ int main(int argc, char *argv[])
buff.size,
ostream
);

fprintf(stdout, "t_ray_voxel: %es\n", clock2second(t_ray_voxel));
fprintf(stdout, "t_avl_histogram: %es\n", clock2second(t_avl_histogram));

}
mliCorsikaHistogram2d_free(&hist);

chk(mliEventTapeReader_finalize(&arc));
fclose(istream);
fclose(ostream);

mliDynDouble_free(&x_bin_edges);
mliDynDouble_free(&y_bin_edges);
mliDynDouble_free(&z_bin_edges);

mliDynUint32_free(&x_idxs);
mliDynUint32_free(&y_idxs);
mliDynUint32_free(&z_idxs);
mliDynDouble_free(&overlaps);

return EXIT_SUCCESS;
chk_error:
return EXIT_FAILURE;
Expand Down

0 comments on commit 45f9e97

Please sign in to comment.