Skip to content

Commit

Permalink
Merge pull request #275 from Magritte-code/imager_correction
Browse files Browse the repository at this point in the history 
Initial commit for better imaging ray tracing.
  • Loading branch information
@ThomasCeulemans
ThomasCeulemans authored Oct 16, 2024
2 parents c89eaf9 + d5c2e2b commit 329ac79
Show file tree
Hide file tree
Showing 4 changed files with 94 additions and 111 deletions.
46 changes: 23 additions & 23 deletions src/model/geometry/geometry.tpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -58,9 +58,9 @@ accel inline Size Geometry ::get_next_general_geometry<Imagetracer>(
const Size n_nbs = points.n_neighbors[c];
const Size cum_n_nbs = points.cum_n_neighbors[c];

double dmin = std::numeric_limits<Real>::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<Real>::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];
Expand All @@ -76,28 +76,28 @@ accel inline Size Geometry ::get_next_general_geometry<Imagetracer>(
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;
Expand Down
16 changes: 16 additions & 0 deletions src/model/geometry/points/points.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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<Real>::max();
Real max_x = std::numeric_limits<Real>::min();
Real min_y = std::numeric_limits<Real>::max();
Real max_y = std::numeric_limits<Real>::min();
Real min_z = std::numeric_limits<Real>::max();
Real max_z = std::numeric_limits<Real>::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());
Expand Down
2 changes: 2 additions & 0 deletions src/model/geometry/points/points.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -10,6 +10,8 @@ struct Points {
Vector<Vector3D> position; ///< position vectors of each point
Vector<Vector3D> velocity; ///< velocity vectors of each point

Vector3D center; ///< center of the model

Vector<Size> cum_n_neighbors; ///< cumulative number of neighbors
Vector<Size> n_neighbors; ///< number of neighbors
Vector<Size> neighbors; ///< neighbors of each point
Expand Down
141 changes: 53 additions & 88 deletions src/solver/solver.tpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -852,53 +852,47 @@ accel inline Size Solver ::trace_ray_imaging(const Model& model, const Vector3D&

Size nxt = model.geometry.get_next<Imagetracer>(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<frame>(origin, origin_velocity, raydir, crt, Z, false);
model.geometry.get_shift<frame>(origin, origin_velocity, raydir, crt, Z - dZ, false);
double shift_nxt =
model.geometry.get_shift<frame>(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<Imagetracer>(
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<Imagetracer>(origin, raydir, nxt, Z, dZ);

if (!model.geometry.valid_point(nxt)) {
break;
}

shift_nxt =
model.geometry.get_shift<frame>(origin, origin_velocity, raydir, nxt, Z, false);

set_data(model, crt, nxt, shift_crt, shift_nxt, dZ, increment, id1, id2);
}
}

// 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;
}

Expand All @@ -920,53 +914,43 @@ accel inline Size Solver ::trace_ray_imaging_for_line(const Model& model, const

Size nxt = model.geometry.get_next<Imagetracer>(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<frame>(origin, origin_velocity, raydir, crt, Z, false);
double shift_nxt =
model.geometry.get_shift<frame>(origin, origin_velocity, raydir, nxt, Z, false);
model.geometry.get_shift<frame>(origin, origin_velocity, raydir, crt, Z - dZ, false);
double shift_nxt = model.geometry.get_shift<frame>(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<Imagetracer>(
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<Imagetracer>(origin, raydir, nxt, Z, dZ);
shift_nxt =
model.geometry.get_shift<frame>(origin, origin_velocity, raydir, nxt, Z, false);

if (!model.geometry.valid_point(nxt)) {
break;
}

shift_nxt = model.geometry.get_shift<frame>(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;
}

Expand All @@ -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<Imagetracer>(
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<Imagetracer>(origin, raydir, nxt, Z, dZ);
if (!model.geometry.valid_point(nxt)) {
break;
}
l += interp_helper.get_n_interp(model, crt, nxt);
}
}

return l;
}

Expand Down

0 comments on commit 329ac79

Please sign in to comment.