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..f4a60b0f 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,11 @@ accel inline Size Solver ::trace_ray_imaging(const Model& model, const Vector3D& } } + // 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; } @@ -920,53 +914,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 +973,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; }