From 5ec250047f2021374e043395c72fb5a546178e04 Mon Sep 17 00:00:00 2001 From: Thomas Ceulemans Date: Fri, 11 Oct 2024 15:21:00 +0200 Subject: [PATCH 1/2] Initial commit for better imaging ray tracing. Now only stops when encountering boundary points on the opposite side of the domain. Will now also not compute anything if only boundary points are encountered. --- src/model/geometry/geometry.tpp | 46 ++++----- src/model/geometry/points/points.cpp | 16 +++ src/model/geometry/points/points.hpp | 2 + src/solver/solver.tpp | 139 ++++++++++----------------- 4 files changed, 92 insertions(+), 111 deletions(-) diff --git a/src/model/geometry/geometry.tpp b/src/model/geometry/geometry.tpp index 5d10d435..fdff3578 100644 --- a/src/model/geometry/geometry.tpp +++ b/src/model/geometry/geometry.tpp @@ -58,9 +58,9 @@ accel inline Size Geometry ::get_next_general_geometry( const Size n_nbs = points.n_neighbors[c]; const Size cum_n_nbs = points.cum_n_neighbors[c]; - double dmin = std::numeric_limits::max(); // Initialize to "infinity" - Size next = parameters->npoints(); // return npoints when there is no next - double maxdist_neighbors2 = 0.0; + double dmin = std::numeric_limits::max(); // Initialize to "infinity" + Size next = parameters->npoints(); // return npoints when there is no next + // double maxdist_neighbors2 = 0.0; for (Size i = 0; i < n_nbs; i++) { const Size n = points.neighbors[cum_n_nbs + i]; @@ -76,28 +76,28 @@ accel inline Size Geometry ::get_next_general_geometry( dZ = Z_new - Z; // such that dZ > 0.0 } } - const Vector3D diff_nxt_crt = points.position[n] - points.position[c]; - const double dist_neighbor2 = diff_nxt_crt.dot(diff_nxt_crt); - if (dist_neighbor2 > maxdist_neighbors2) { - maxdist_neighbors2 = dist_neighbor2; - } + // const Vector3D diff_nxt_crt = points.position[n] - points.position[c]; + // const double dist_neighbor2 = diff_nxt_crt.dot(diff_nxt_crt); + // if (dist_neighbor2 > maxdist_neighbors2) { + // maxdist_neighbors2 = dist_neighbor2; + // } } - const Vector3D R_c = points.position[c] - origin; - const double Z_c = R_c.dot(raydir); - const double distance_curr_from_ray2 = R_c.dot(R_c) - Z_c * Z_c; - - // Precaution against tracing stuff along the boundary, stopping when we move - // farther from the ray than the distance traveled using a safety factor of 2; - // TO DO: actually implement some manner of 'surface unit vector' to determine - // when to stop this ray; then this entire condition might be replaced by this - // (maybe also merge with 'regular' get_next_general_geometry) we also check - // whether the distance from the ray increases (as irregularly placed boundary - // points might result in stopping too early) - if ((!not_on_boundary(c)) && (2.0 * maxdist_neighbors2 < dmin) - && (distance_curr_from_ray2 < dmin)) { - return parameters->npoints(); - } + // const Vector3D R_c = points.position[c] - origin; + // const double Z_c = R_c.dot(raydir); + // const double distance_curr_from_ray2 = R_c.dot(R_c) - Z_c * Z_c; + + // // Precaution against tracing stuff along the boundary, stopping when we move + // // farther from the ray than the distance traveled using a safety factor of 2; + // // TO DO: actually implement some manner of 'surface unit vector' to determine + // // when to stop this ray; then this entire condition might be replaced by this + // // (maybe also merge with 'regular' get_next_general_geometry) we also check + // // whether the distance from the ray increases (as irregularly placed boundary + // // points might result in stopping too early) + // if ((!not_on_boundary(c)) && (2.0 * maxdist_neighbors2 < dmin) + // && (distance_curr_from_ray2 < dmin)) { + // return parameters->npoints(); + // } // Update distance along ray Z += dZ; diff --git a/src/model/geometry/points/points.cpp b/src/model/geometry/points/points.cpp index 6c77adea..f7950f49 100644 --- a/src/model/geometry/points/points.cpp +++ b/src/model/geometry/points/points.cpp @@ -19,12 +19,28 @@ void Points ::read(const Io& io) { io.read_array(prefix + "position", position_buffer); io.read_array(prefix + "velocity", velocity_buffer); + Real min_x = std::numeric_limits::max(); + Real max_x = std::numeric_limits::min(); + Real min_y = std::numeric_limits::max(); + Real max_y = std::numeric_limits::min(); + Real min_z = std::numeric_limits::max(); + Real max_z = std::numeric_limits::min(); + for (Size p = 0; p < parameters->npoints(); p++) { position[p] = Vector3D(position_buffer[p][0], position_buffer[p][1], position_buffer[p][2]); velocity[p] = Vector3D(velocity_buffer[p][0], velocity_buffer[p][1], velocity_buffer[p][2]); + + if (position[p].x() < min_x) min_x = position[p].x(); + if (position[p].x() > max_x) max_x = position[p].x(); + if (position[p].y() < min_y) min_y = position[p].y(); + if (position[p].y() > max_y) max_y = position[p].y(); + if (position[p].z() < min_z) min_z = position[p].z(); + if (position[p].z() > max_z) max_z = position[p].z(); } + center = Vector3D((min_x + max_x) / 2, (min_y + max_y) / 2, (min_z + max_z) / 2); + const Size totnnbs = io.get_length(prefix + "neighbors"); cum_n_neighbors.resize(parameters->npoints()); diff --git a/src/model/geometry/points/points.hpp b/src/model/geometry/points/points.hpp index ad64c9b1..edaef3c7 100644 --- a/src/model/geometry/points/points.hpp +++ b/src/model/geometry/points/points.hpp @@ -10,6 +10,8 @@ struct Points { Vector position; ///< position vectors of each point Vector velocity; ///< velocity vectors of each point + Vector3D center; ///< center of the model + Vector cum_n_neighbors; ///< cumulative number of neighbors Vector n_neighbors; ///< number of neighbors Vector neighbors; ///< neighbors of each point diff --git a/src/solver/solver.tpp b/src/solver/solver.tpp index 40e1f45a..b1c0d644 100644 --- a/src/solver/solver.tpp +++ b/src/solver/solver.tpp @@ -852,46 +852,35 @@ accel inline Size Solver ::trace_ray_imaging(const Model& model, const Vector3D& Size nxt = model.geometry.get_next(origin, raydir, start_bdy, Z, dZ); + bool accessed_non_boundary = false; + Size original_id1 = + id1; // initial ray index to return in case of only boundary points on the ray + if (model.geometry.valid_point(nxt)) { + // Size crt = o; double shift_crt = - model.geometry.get_shift(origin, origin_velocity, raydir, crt, Z, false); + model.geometry.get_shift(origin, origin_velocity, raydir, crt, Z - dZ, false); double shift_nxt = model.geometry.get_shift(origin, origin_velocity, raydir, nxt, Z, false); set_data(model, crt, nxt, shift_crt, shift_nxt, dZ, increment, id1, id2); - // Due to boundary points begin slightly annoying to - // start a ray from, we must make sure the ray only end - // when we are sure that we stay on the boundary - while (true) { - if (!model.geometry.not_on_boundary(nxt)) { - Size curr_cons_bdy = 1; - Size temp_nxt = nxt; - double temp_Z = Z; - double temp_dZ = dZ; - - while (true) { - temp_nxt = model.geometry.get_next( - origin, raydir, temp_nxt, temp_Z, temp_dZ); - if ((!model.geometry.valid_point(temp_nxt)) - || (curr_cons_bdy == MAX_CONSECUTIVE_BDY)) { - return id1; // the ray ends if we cannot find - // any points anymore, or we have - // too many boundary points after - // eachother - } - if (model.geometry.not_on_boundary(temp_nxt)) { - break; // the ray continues, as we find once - // again non-boundary points - } - curr_cons_bdy += 1; - } - } - + while (model.geometry.not_on_boundary(nxt) + || (model.geometry.points.position[nxt] - model.geometry.points.center).dot(raydir) + < 0) { crt = nxt; shift_crt = shift_nxt; + if (model.geometry.not_on_boundary(nxt)) { + accessed_non_boundary = true; + } + nxt = model.geometry.get_next(origin, raydir, nxt, Z, dZ); + + if (!model.geometry.valid_point(nxt)) { + break; + } + shift_nxt = model.geometry.get_shift(origin, origin_velocity, raydir, nxt, Z, false); @@ -899,6 +888,9 @@ accel inline Size Solver ::trace_ray_imaging(const Model& model, const Vector3D& } } + if (!accessed_non_boundary) { + return original_id1; + } return id1; } @@ -920,53 +912,43 @@ accel inline Size Solver ::trace_ray_imaging_for_line(const Model& model, const Size nxt = model.geometry.get_next(origin, raydir, start_bdy, Z, dZ); + bool accessed_non_boundary = false; + Size original_id1 = + id1; // initial ray index to return in case of only boundary points on the rays + if (model.geometry.valid_point(nxt)) { double shift_crt = - model.geometry.get_shift(origin, origin_velocity, raydir, crt, Z, false); - double shift_nxt = - model.geometry.get_shift(origin, origin_velocity, raydir, nxt, Z, false); + model.geometry.get_shift(origin, origin_velocity, raydir, crt, Z - dZ, false); + double shift_nxt = model.geometry.get_shift(origin, origin_velocity, raydir, nxt, Z); set_data_for_line(model, l, crt, nxt, shift_crt, shift_nxt, dZ, increment, id1, id2); - // Due to boundary points begin slightly annoying to - // start a ray from, we must make sure the ray only end - // when we are sure that we stay on the boundary - while (true) { - if (!model.geometry.not_on_boundary(nxt)) { - Size curr_cons_bdy = 1; - Size temp_nxt = nxt; - double temp_Z = Z; - double temp_dZ = dZ; - - while (true) { - temp_nxt = model.geometry.get_next( - origin, raydir, temp_nxt, temp_Z, temp_dZ); - if ((!model.geometry.valid_point(temp_nxt)) - || (curr_cons_bdy == MAX_CONSECUTIVE_BDY)) { - return id1; // the ray ends if we cannot find - // any points anymore, or we have - // too many boundary points after - // eachother - } - if (model.geometry.not_on_boundary(temp_nxt)) { - break; // the ray continues, as we find once - // again non-boundary points - } - curr_cons_bdy += 1; - } - } - + while (model.geometry.not_on_boundary(nxt) + || (model.geometry.points.position[nxt] - model.geometry.points.center).dot(raydir) + < 0) { crt = nxt; shift_crt = shift_nxt; + if (model.geometry.not_on_boundary(nxt)) { + accessed_non_boundary = true; + } + nxt = model.geometry.get_next(origin, raydir, nxt, Z, dZ); - shift_nxt = - model.geometry.get_shift(origin, origin_velocity, raydir, nxt, Z, false); + + if (!model.geometry.valid_point(nxt)) { + break; + } + + shift_nxt = model.geometry.get_shift(origin, origin_velocity, raydir, nxt, Z); set_data_for_line(model, l, crt, nxt, shift_crt, shift_nxt, dZ, increment, id1, id2); } } + if (!accessed_non_boundary) { + return original_id1; + } + return id1; } @@ -989,37 +971,18 @@ accel inline Size Solver ::get_ray_length_new_imager( if (model.geometry.valid_point(nxt)) { l += interp_helper.get_n_interp(model, crt, nxt); - // Due to boundary points begin slightly annoying to - // start a ray from, we must make sure the ray only end - // when we are sure that we stay on the boundary - while (true) { - if (!model.geometry.not_on_boundary(nxt)) { - Size curr_cons_bdy = 1; - Size temp_nxt = nxt; - double temp_Z = Z; - double temp_dZ = dZ; - while (true) { - temp_nxt = model.geometry.get_next( - origin, raydir, temp_nxt, temp_Z, temp_dZ); - if ((!model.geometry.valid_point(temp_nxt)) - || (curr_cons_bdy == MAX_CONSECUTIVE_BDY)) { - return l; // the ray ends if we cannot find any - // points anymore, or we have too many - // boundary points after eachother - } - if (model.geometry.not_on_boundary(temp_nxt)) { - break; // the ray continues, as we find once - // again non-boundary points - } - curr_cons_bdy += interp_helper.get_n_interp(model, temp_nxt, temp_nxt); - } - } - + while (model.geometry.not_on_boundary(nxt) + || (model.geometry.points.position[nxt] - model.geometry.points.center).dot(raydir) + < 0) { crt = nxt; nxt = model.geometry.get_next(origin, raydir, nxt, Z, dZ); + if (!model.geometry.valid_point(nxt)) { + break; + } l += interp_helper.get_n_interp(model, crt, nxt); } } + return l; } From d5c2e2b6e5b890f5548e4a79e60c38ecb44fce14 Mon Sep 17 00:00:00 2001 From: Thomas Ceulemans Date: Wed, 16 Oct 2024 15:58:06 +0200 Subject: [PATCH 2/2] Added boundary treatment fix for 1D imager. --- src/solver/solver.tpp | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/src/solver/solver.tpp b/src/solver/solver.tpp index b1c0d644..f4a60b0f 100644 --- a/src/solver/solver.tpp +++ b/src/solver/solver.tpp @@ -888,7 +888,9 @@ accel inline Size Solver ::trace_ray_imaging(const Model& model, const Vector3D& } } - if (!accessed_non_boundary) { + // Spherically symmetric models have a correct boundary treatment already, thus no need to + // reduce their accuracy in a region between the outer boundary and the nearest point. + if (!accessed_non_boundary && !model.parameters->spherical_symmetry()) { return original_id1; } return id1;