mirror of
https://github.com/openmc-dev/openmc.git
synced 2026-07-26 13:15:39 -04:00
Raytrace plots (#2655)
Co-authored-by: Patrick Shriwise <pshriwise@gmail.com> Co-authored-by: Paul Romano <paul.k.romano@gmail.com>
This commit is contained in:
Gavin Ridley
GitHub
parent
a5b26de041
commit
81b7388624
25 changed files with 1487 additions and 476 deletions
674
src/plot.cpp
674
src/plot.cpp
|
|
@ -204,18 +204,21 @@ void read_plots_xml(pugi::xml_node root)
|
|||
int id = std::stoi(id_string);
|
||||
if (check_for_node(node, "type")) {
|
||||
std::string type_str = get_node_value(node, "type", true);
|
||||
if (type_str == "slice")
|
||||
if (type_str == "slice") {
|
||||
model::plots.emplace_back(
|
||||
std::make_unique<Plot>(node, Plot::PlotType::slice));
|
||||
else if (type_str == "voxel")
|
||||
} else if (type_str == "voxel") {
|
||||
model::plots.emplace_back(
|
||||
std::make_unique<Plot>(node, Plot::PlotType::voxel));
|
||||
else if (type_str == "projection")
|
||||
model::plots.emplace_back(std::make_unique<ProjectionPlot>(node));
|
||||
else
|
||||
} else if (type_str == "wireframe_raytrace") {
|
||||
model::plots.emplace_back(
|
||||
std::make_unique<WireframeRayTracePlot>(node));
|
||||
} else if (type_str == "solid_raytrace") {
|
||||
model::plots.emplace_back(std::make_unique<SolidRayTracePlot>(node));
|
||||
} else {
|
||||
fatal_error(
|
||||
fmt::format("Unsupported plot type '{}' in plot {}", type_str, id));
|
||||
|
||||
}
|
||||
model::plot_map[model::plots.back()->id()] = model::plots.size() - 1;
|
||||
} else {
|
||||
fatal_error(fmt::format("Must specify plot type in plot {}", id));
|
||||
|
|
@ -264,8 +267,8 @@ void Plot::create_image() const
|
|||
}
|
||||
data(x, y) = colors_[model::material_map[id]];
|
||||
} // color_by if-else
|
||||
} // x for loop
|
||||
} // y for loop
|
||||
}
|
||||
}
|
||||
|
||||
// draw mesh lines if present
|
||||
if (index_meshlines_mesh_ >= 0) {
|
||||
|
|
@ -1036,26 +1039,49 @@ RGBColor random_color(void)
|
|||
int(prn(&model::plotter_seed) * 255), int(prn(&model::plotter_seed) * 255)};
|
||||
}
|
||||
|
||||
ProjectionPlot::ProjectionPlot(pugi::xml_node node) : PlottableInterface(node)
|
||||
RayTracePlot::RayTracePlot(pugi::xml_node node) : PlottableInterface(node)
|
||||
{
|
||||
set_output_path(node);
|
||||
set_look_at(node);
|
||||
set_camera_position(node);
|
||||
set_field_of_view(node);
|
||||
set_pixels(node);
|
||||
set_opacities(node);
|
||||
set_orthographic_width(node);
|
||||
set_wireframe_thickness(node);
|
||||
set_wireframe_ids(node);
|
||||
set_wireframe_color(node);
|
||||
set_output_path(node);
|
||||
|
||||
if (check_for_node(node, "orthographic_width") &&
|
||||
check_for_node(node, "field_of_view"))
|
||||
fatal_error("orthographic_width and field_of_view are mutually exclusive "
|
||||
"parameters.");
|
||||
|
||||
// Get centerline vector for camera-to-model. We create vectors around this
|
||||
// that form a pixel array, and then trace rays along that.
|
||||
auto up = up_ / up_.norm();
|
||||
Direction looking_direction = look_at_ - camera_position_;
|
||||
looking_direction /= looking_direction.norm();
|
||||
if (std::abs(std::abs(looking_direction.dot(up)) - 1.0) < 1e-9)
|
||||
fatal_error("Up vector cannot align with vector between camera position "
|
||||
"and look_at!");
|
||||
Direction cam_yaxis = looking_direction.cross(up);
|
||||
cam_yaxis /= cam_yaxis.norm();
|
||||
Direction cam_zaxis = cam_yaxis.cross(looking_direction);
|
||||
cam_zaxis /= cam_zaxis.norm();
|
||||
|
||||
// Cache the camera-to-model matrix
|
||||
camera_to_model_ = {looking_direction.x, cam_yaxis.x, cam_zaxis.x,
|
||||
looking_direction.y, cam_yaxis.y, cam_zaxis.y, looking_direction.z,
|
||||
cam_yaxis.z, cam_zaxis.z};
|
||||
}
|
||||
|
||||
void ProjectionPlot::set_wireframe_color(pugi::xml_node plot_node)
|
||||
WireframeRayTracePlot::WireframeRayTracePlot(pugi::xml_node node)
|
||||
: RayTracePlot(node)
|
||||
{
|
||||
set_opacities(node);
|
||||
set_wireframe_thickness(node);
|
||||
set_wireframe_ids(node);
|
||||
set_wireframe_color(node);
|
||||
}
|
||||
|
||||
void WireframeRayTracePlot::set_wireframe_color(pugi::xml_node plot_node)
|
||||
{
|
||||
// Copy plot background color
|
||||
if (check_for_node(plot_node, "wireframe_color")) {
|
||||
|
|
@ -1068,7 +1094,7 @@ void ProjectionPlot::set_wireframe_color(pugi::xml_node plot_node)
|
|||
}
|
||||
}
|
||||
|
||||
void ProjectionPlot::set_output_path(pugi::xml_node node)
|
||||
void RayTracePlot::set_output_path(pugi::xml_node node)
|
||||
{
|
||||
// Set output file path
|
||||
std::string filename;
|
||||
|
|
@ -1089,33 +1115,7 @@ void ProjectionPlot::set_output_path(pugi::xml_node node)
|
|||
path_plot_ = filename;
|
||||
}
|
||||
|
||||
// Advances to the next boundary from outside the geometry
|
||||
// Returns -1 if no intersection found, and the surface index
|
||||
// if an intersection was found.
|
||||
int ProjectionPlot::advance_to_boundary_from_void(GeometryState& p)
|
||||
{
|
||||
constexpr double scoot = 1e-5;
|
||||
double min_dist = {INFINITY};
|
||||
auto coord = p.coord(0);
|
||||
Universe* uni = model::universes[model::root_universe].get();
|
||||
int intersected_surface = -1;
|
||||
for (auto c_i : uni->cells_) {
|
||||
auto dist = model::cells.at(c_i)->distance(coord.r, coord.u, 0, &p);
|
||||
if (dist.first < min_dist) {
|
||||
min_dist = dist.first;
|
||||
intersected_surface = dist.second;
|
||||
}
|
||||
}
|
||||
if (min_dist > 1e300)
|
||||
return -1;
|
||||
else { // advance the particle
|
||||
for (int j = 0; j < p.n_coord(); ++j)
|
||||
p.coord(j).r += (min_dist + scoot) * p.coord(j).u;
|
||||
return std::abs(intersected_surface);
|
||||
}
|
||||
}
|
||||
|
||||
bool ProjectionPlot::trackstack_equivalent(
|
||||
bool WireframeRayTracePlot::trackstack_equivalent(
|
||||
const std::vector<TrackSegment>& track1,
|
||||
const std::vector<TrackSegment>& track2) const
|
||||
{
|
||||
|
|
@ -1125,7 +1125,7 @@ bool ProjectionPlot::trackstack_equivalent(
|
|||
return false;
|
||||
for (int i = 0; i < track1.size(); ++i) {
|
||||
if (track1[i].id != track2[i].id ||
|
||||
track1[i].surface != track2[i].surface) {
|
||||
track1[i].surface_index != track2[i].surface_index) {
|
||||
return false;
|
||||
}
|
||||
}
|
||||
|
|
@ -1152,7 +1152,7 @@ bool ProjectionPlot::trackstack_equivalent(
|
|||
if (t1_i == track1.size() && t2_i == track2.size())
|
||||
break;
|
||||
// Check if surface different
|
||||
if (track1[t1_i].surface != track2[t2_i].surface)
|
||||
if (track1[t1_i].surface_index != track2[t2_i].surface_index)
|
||||
return false;
|
||||
|
||||
// Pretty sure this should not be used:
|
||||
|
|
@ -1160,7 +1160,7 @@ bool ProjectionPlot::trackstack_equivalent(
|
|||
// t1_i != track1.size() - 1 &&
|
||||
// track1[t1_i+1].id != track2[t2_i+1].id) return false;
|
||||
if (t2_i != 0 && t1_i != 0 &&
|
||||
track1[t1_i - 1].surface != track2[t2_i - 1].surface)
|
||||
track1[t1_i - 1].surface_index != track2[t2_i - 1].surface_index)
|
||||
return false;
|
||||
|
||||
// Check if neighboring cells are different
|
||||
|
|
@ -1176,44 +1176,58 @@ bool ProjectionPlot::trackstack_equivalent(
|
|||
}
|
||||
}
|
||||
|
||||
void ProjectionPlot::create_output() const
|
||||
std::pair<Position, Direction> RayTracePlot::get_pixel_ray(
|
||||
int horiz, int vert) const
|
||||
{
|
||||
// Get centerline vector for camera-to-model. We create vectors around this
|
||||
// that form a pixel array, and then trace rays along that.
|
||||
auto up = up_ / up_.norm();
|
||||
Direction looking_direction = look_at_ - camera_position_;
|
||||
looking_direction /= looking_direction.norm();
|
||||
if (std::abs(std::abs(looking_direction.dot(up)) - 1.0) < 1e-9)
|
||||
fatal_error("Up vector cannot align with vector between camera position "
|
||||
"and look_at!");
|
||||
Direction cam_yaxis = looking_direction.cross(up);
|
||||
cam_yaxis /= cam_yaxis.norm();
|
||||
Direction cam_zaxis = cam_yaxis.cross(looking_direction);
|
||||
cam_zaxis /= cam_zaxis.norm();
|
||||
|
||||
// Transformation matrix for directions
|
||||
std::vector<double> camera_to_model = {looking_direction.x, cam_yaxis.x,
|
||||
cam_zaxis.x, looking_direction.y, cam_yaxis.y, cam_zaxis.y,
|
||||
looking_direction.z, cam_yaxis.z, cam_zaxis.z};
|
||||
|
||||
// Now we convert to the polar coordinate system with the polar angle
|
||||
// measuring the angle from the vector up_. Phi is the rotation about up_. For
|
||||
// now, up_ is hard-coded to be +z.
|
||||
// Compute field of view in radians
|
||||
constexpr double DEGREE_TO_RADIAN = M_PI / 180.0;
|
||||
double horiz_fov_radians = horizontal_field_of_view_ * DEGREE_TO_RADIAN;
|
||||
double p0 = static_cast<double>(pixels_[0]);
|
||||
double p1 = static_cast<double>(pixels_[1]);
|
||||
double vert_fov_radians = horiz_fov_radians * p1 / p0;
|
||||
double dphi = horiz_fov_radians / p0;
|
||||
double dmu = vert_fov_radians / p1;
|
||||
|
||||
// focal_plane_dist can be changed to alter the perspective distortion
|
||||
// effect. This is in units of cm. This seems to look good most of the
|
||||
// time. TODO let this variable be set through XML.
|
||||
constexpr double focal_plane_dist = 10.0;
|
||||
const double dx = 2.0 * focal_plane_dist * std::tan(0.5 * horiz_fov_radians);
|
||||
const double dy = p1 / p0 * dx;
|
||||
|
||||
std::pair<Position, Direction> result;
|
||||
|
||||
// Generate the starting position/direction of the ray
|
||||
if (orthographic_width_ == C_NONE) { // perspective projection
|
||||
Direction camera_local_vec;
|
||||
camera_local_vec.x = focal_plane_dist;
|
||||
camera_local_vec.y = -0.5 * dx + horiz * dx / p0;
|
||||
camera_local_vec.z = 0.5 * dy - vert * dy / p1;
|
||||
camera_local_vec /= camera_local_vec.norm();
|
||||
|
||||
result.first = camera_position_;
|
||||
result.second = camera_local_vec.rotate(camera_to_model_);
|
||||
} else { // orthographic projection
|
||||
|
||||
double x_pix_coord = (static_cast<double>(horiz) - p0 / 2.0) / p0;
|
||||
double y_pix_coord = (static_cast<double>(vert) - p1 / 2.0) / p1;
|
||||
|
||||
result.first = camera_position_ +
|
||||
camera_y_axis() * x_pix_coord * orthographic_width_ +
|
||||
camera_z_axis() * y_pix_coord * orthographic_width_;
|
||||
result.second = camera_x_axis();
|
||||
}
|
||||
|
||||
return result;
|
||||
}
|
||||
|
||||
void WireframeRayTracePlot::create_output() const
|
||||
{
|
||||
size_t width = pixels_[0];
|
||||
size_t height = pixels_[1];
|
||||
ImageData data({width, height}, not_found_);
|
||||
|
||||
// This array marks where the initial wireframe was drawn.
|
||||
// We convolve it with a filter that gets adjusted with the
|
||||
// wireframe thickness in order to thicken the lines.
|
||||
// This array marks where the initial wireframe was drawn. We convolve it with
|
||||
// a filter that gets adjusted with the wireframe thickness in order to
|
||||
// thicken the lines.
|
||||
xt::xtensor<int, 2> wireframe_initial({width, height}, 0);
|
||||
|
||||
/* Holds all of the track segments for the current rendered line of pixels.
|
||||
|
|
@ -1242,16 +1256,13 @@ void ProjectionPlot::create_output() const
|
|||
const int n_threads = num_threads();
|
||||
const int tid = thread_num();
|
||||
|
||||
GeometryState p;
|
||||
p.u() = {1.0, 0.0, 0.0};
|
||||
|
||||
int vert = tid;
|
||||
for (int iter = 0; iter <= pixels_[1] / n_threads; iter++) {
|
||||
|
||||
// Save bottom line of current work chunk to compare against later
|
||||
// I used to have this inside the below if block, but it causes a
|
||||
// spurious line to be drawn at the bottom of the image. Not sure
|
||||
// why, but moving it here fixes things.
|
||||
// Save bottom line of current work chunk to compare against later. This
|
||||
// used to be inside the below if block, but it causes a spurious line to
|
||||
// be drawn at the bottom of the image. Not sure why, but moving it here
|
||||
// fixes things.
|
||||
if (tid == n_threads - 1)
|
||||
old_segments = this_line_segments[n_threads - 1];
|
||||
|
||||
|
|
@ -1259,126 +1270,48 @@ void ProjectionPlot::create_output() const
|
|||
|
||||
for (int horiz = 0; horiz < pixels_[0]; ++horiz) {
|
||||
|
||||
// Projection mode below decides ray starting conditions
|
||||
Position init_r;
|
||||
Direction init_u;
|
||||
|
||||
// Generate the starting position/direction of the ray
|
||||
if (orthographic_width_ == 0.0) { // perspective projection
|
||||
double this_phi =
|
||||
-horiz_fov_radians / 2.0 + dphi * horiz + 0.5 * dphi;
|
||||
double this_mu =
|
||||
-vert_fov_radians / 2.0 + dmu * vert + M_PI / 2.0 + 0.5 * dmu;
|
||||
Direction camera_local_vec;
|
||||
camera_local_vec.x = std::cos(this_phi) * std::sin(this_mu);
|
||||
camera_local_vec.y = std::sin(this_phi) * std::sin(this_mu);
|
||||
camera_local_vec.z = std::cos(this_mu);
|
||||
init_u = camera_local_vec.rotate(camera_to_model);
|
||||
init_r = camera_position_;
|
||||
} else { // orthographic projection
|
||||
init_u = looking_direction;
|
||||
|
||||
double x_pix_coord = (static_cast<double>(horiz) - p0 / 2.0) / p0;
|
||||
double y_pix_coord = (static_cast<double>(vert) - p1 / 2.0) / p0;
|
||||
|
||||
init_r = camera_position_;
|
||||
init_r += cam_yaxis * x_pix_coord * orthographic_width_;
|
||||
init_r += cam_zaxis * y_pix_coord * orthographic_width_;
|
||||
}
|
||||
|
||||
// Resets internal geometry state of particle
|
||||
p.init_from_r_u(init_r, init_u);
|
||||
|
||||
bool hitsomething = false;
|
||||
bool intersection_found = true;
|
||||
int loop_counter = 0;
|
||||
// RayTracePlot implements camera ray generation
|
||||
std::pair<Position, Direction> ru = get_pixel_ray(horiz, vert);
|
||||
|
||||
this_line_segments[tid][horiz].clear();
|
||||
ProjectionRay ray(
|
||||
ru.first, ru.second, *this, this_line_segments[tid][horiz]);
|
||||
|
||||
int first_surface =
|
||||
-1; // surface first passed when entering the model
|
||||
bool first_inside_model = true; // false after entering the model
|
||||
while (intersection_found) {
|
||||
bool inside_cell = false;
|
||||
|
||||
int32_t i_surface = p.surface_index();
|
||||
if (i_surface > 0 &&
|
||||
model::surfaces[i_surface]->geom_type() == GeometryType::DAG) {
|
||||
#ifdef DAGMC
|
||||
int32_t i_cell = next_cell(i_surface,
|
||||
p.cell_last(p.n_coord() - 1), p.lowest_coord().universe);
|
||||
inside_cell = i_cell >= 0;
|
||||
#else
|
||||
fatal_error(
|
||||
"Not compiled for DAGMC, but somehow you have a DAGCell!");
|
||||
#endif
|
||||
} else {
|
||||
inside_cell = exhaustive_find_cell(p);
|
||||
}
|
||||
|
||||
if (inside_cell) {
|
||||
|
||||
// This allows drawing wireframes with surface intersection
|
||||
// edges on the model boundary for the same cell.
|
||||
if (first_inside_model) {
|
||||
this_line_segments[tid][horiz].emplace_back(
|
||||
color_by_ == PlotColorBy::mats ? p.material()
|
||||
: p.lowest_coord().cell,
|
||||
0.0, first_surface);
|
||||
first_inside_model = false;
|
||||
}
|
||||
|
||||
hitsomething = true;
|
||||
intersection_found = true;
|
||||
auto dist = distance_to_boundary(p);
|
||||
this_line_segments[tid][horiz].emplace_back(
|
||||
color_by_ == PlotColorBy::mats ? p.material()
|
||||
: p.lowest_coord().cell,
|
||||
dist.distance, std::abs(dist.surface));
|
||||
|
||||
// Advance particle
|
||||
for (int lev = 0; lev < p.n_coord(); ++lev) {
|
||||
p.coord(lev).r += dist.distance * p.coord(lev).u;
|
||||
}
|
||||
p.surface() = dist.surface;
|
||||
p.n_coord_last() = p.n_coord();
|
||||
p.n_coord() = dist.coord_level;
|
||||
if (dist.lattice_translation[0] != 0 ||
|
||||
dist.lattice_translation[1] != 0 ||
|
||||
dist.lattice_translation[2] != 0) {
|
||||
cross_lattice(p, dist);
|
||||
}
|
||||
|
||||
} else {
|
||||
first_surface = advance_to_boundary_from_void(p);
|
||||
intersection_found =
|
||||
first_surface != -1; // -1 if no surface found
|
||||
}
|
||||
loop_counter++;
|
||||
if (loop_counter > MAX_INTERSECTIONS)
|
||||
fatal_error("Infinite loop in projection plot");
|
||||
}
|
||||
ray.trace();
|
||||
|
||||
// Now color the pixel based on what we have intersected...
|
||||
// Loops backwards over intersections.
|
||||
Position current_color(
|
||||
not_found_.red, not_found_.green, not_found_.blue);
|
||||
const auto& segments = this_line_segments[tid][horiz];
|
||||
for (unsigned i = segments.size(); i-- > 0;) {
|
||||
|
||||
// There must be at least two cell intersections to color, front and
|
||||
// back of the cell. Maybe an infinitely thick cell could be present
|
||||
// with no back, but why would you want to color that? It's easier to
|
||||
// just skip that edge case and not even color it.
|
||||
if (segments.size() <= 1)
|
||||
continue;
|
||||
|
||||
for (int i = segments.size() - 2; i >= 0; --i) {
|
||||
int colormap_idx = segments[i].id;
|
||||
RGBColor seg_color = colors_[colormap_idx];
|
||||
Position seg_color_vec(
|
||||
seg_color.red, seg_color.green, seg_color.blue);
|
||||
double mixing = std::exp(-xs_[colormap_idx] * segments[i].length);
|
||||
double mixing =
|
||||
std::exp(-xs_[colormap_idx] *
|
||||
(segments[i + 1].length - segments[i].length));
|
||||
current_color =
|
||||
current_color * mixing + (1.0 - mixing) * seg_color_vec;
|
||||
RGBColor result;
|
||||
result.red = static_cast<uint8_t>(current_color.x);
|
||||
result.green = static_cast<uint8_t>(current_color.y);
|
||||
result.blue = static_cast<uint8_t>(current_color.z);
|
||||
data(horiz, vert) = result;
|
||||
}
|
||||
|
||||
// save result converting from double-precision color coordinates to
|
||||
// byte-sized
|
||||
RGBColor result;
|
||||
result.red = static_cast<uint8_t>(current_color.x);
|
||||
result.green = static_cast<uint8_t>(current_color.y);
|
||||
result.blue = static_cast<uint8_t>(current_color.z);
|
||||
data(horiz, vert) = result;
|
||||
|
||||
// Check to draw wireframe in horizontal direction. No inter-thread
|
||||
// comm.
|
||||
if (horiz > 0) {
|
||||
|
|
@ -1451,9 +1384,8 @@ void ProjectionPlot::create_output() const
|
|||
#endif
|
||||
}
|
||||
|
||||
void ProjectionPlot::print_info() const
|
||||
void RayTracePlot::print_info() const
|
||||
{
|
||||
fmt::print("Plot Type: Projection\n");
|
||||
fmt::print("Camera position: {} {} {}\n", camera_position_.x,
|
||||
camera_position_.y, camera_position_.z);
|
||||
fmt::print("Look at: {} {} {}\n", look_at_.x, look_at_.y, look_at_.z);
|
||||
|
|
@ -1462,7 +1394,13 @@ void ProjectionPlot::print_info() const
|
|||
fmt::print("Pixels: {} {}\n", pixels_[0], pixels_[1]);
|
||||
}
|
||||
|
||||
void ProjectionPlot::set_opacities(pugi::xml_node node)
|
||||
void WireframeRayTracePlot::print_info() const
|
||||
{
|
||||
fmt::print("Plot Type: Wireframe ray-traced\n");
|
||||
RayTracePlot::print_info();
|
||||
}
|
||||
|
||||
void WireframeRayTracePlot::set_opacities(pugi::xml_node node)
|
||||
{
|
||||
xs_.resize(colors_.size(), 1e6); // set to large value for opaque by default
|
||||
|
||||
|
|
@ -1492,7 +1430,7 @@ void ProjectionPlot::set_opacities(pugi::xml_node node)
|
|||
}
|
||||
}
|
||||
|
||||
void ProjectionPlot::set_orthographic_width(pugi::xml_node node)
|
||||
void RayTracePlot::set_orthographic_width(pugi::xml_node node)
|
||||
{
|
||||
if (check_for_node(node, "orthographic_width")) {
|
||||
double orthographic_width =
|
||||
|
|
@ -1503,7 +1441,7 @@ void ProjectionPlot::set_orthographic_width(pugi::xml_node node)
|
|||
}
|
||||
}
|
||||
|
||||
void ProjectionPlot::set_wireframe_thickness(pugi::xml_node node)
|
||||
void WireframeRayTracePlot::set_wireframe_thickness(pugi::xml_node node)
|
||||
{
|
||||
if (check_for_node(node, "wireframe_thickness")) {
|
||||
int wireframe_thickness =
|
||||
|
|
@ -1514,7 +1452,7 @@ void ProjectionPlot::set_wireframe_thickness(pugi::xml_node node)
|
|||
}
|
||||
}
|
||||
|
||||
void ProjectionPlot::set_wireframe_ids(pugi::xml_node node)
|
||||
void WireframeRayTracePlot::set_wireframe_ids(pugi::xml_node node)
|
||||
{
|
||||
if (check_for_node(node, "wireframe_ids")) {
|
||||
wireframe_ids_ = get_node_array<int>(node, "wireframe_ids");
|
||||
|
|
@ -1529,7 +1467,7 @@ void ProjectionPlot::set_wireframe_ids(pugi::xml_node node)
|
|||
std::sort(wireframe_ids_.begin(), wireframe_ids_.end());
|
||||
}
|
||||
|
||||
void ProjectionPlot::set_pixels(pugi::xml_node node)
|
||||
void RayTracePlot::set_pixels(pugi::xml_node node)
|
||||
{
|
||||
vector<int> pxls = get_node_array<int>(node, "pixels");
|
||||
if (pxls.size() != 2)
|
||||
|
|
@ -1539,19 +1477,19 @@ void ProjectionPlot::set_pixels(pugi::xml_node node)
|
|||
pixels_[1] = pxls[1];
|
||||
}
|
||||
|
||||
void ProjectionPlot::set_camera_position(pugi::xml_node node)
|
||||
void RayTracePlot::set_camera_position(pugi::xml_node node)
|
||||
{
|
||||
vector<double> camera_pos = get_node_array<double>(node, "camera_position");
|
||||
if (camera_pos.size() != 3) {
|
||||
fatal_error(
|
||||
fmt::format("look_at element must have three floating point values"));
|
||||
fatal_error(fmt::format(
|
||||
"camera_position element must have three floating point values"));
|
||||
}
|
||||
camera_position_.x = camera_pos[0];
|
||||
camera_position_.y = camera_pos[1];
|
||||
camera_position_.z = camera_pos[2];
|
||||
}
|
||||
|
||||
void ProjectionPlot::set_look_at(pugi::xml_node node)
|
||||
void RayTracePlot::set_look_at(pugi::xml_node node)
|
||||
{
|
||||
vector<double> look_at = get_node_array<double>(node, "look_at");
|
||||
if (look_at.size() != 3) {
|
||||
|
|
@ -1562,7 +1500,7 @@ void ProjectionPlot::set_look_at(pugi::xml_node node)
|
|||
look_at_.z = look_at[2];
|
||||
}
|
||||
|
||||
void ProjectionPlot::set_field_of_view(pugi::xml_node node)
|
||||
void RayTracePlot::set_field_of_view(pugi::xml_node node)
|
||||
{
|
||||
// Defaults to 70 degree horizontal field of view (see .h file)
|
||||
if (check_for_node(node, "field_of_view")) {
|
||||
|
|
@ -1576,6 +1514,352 @@ void ProjectionPlot::set_field_of_view(pugi::xml_node node)
|
|||
}
|
||||
}
|
||||
|
||||
SolidRayTracePlot::SolidRayTracePlot(pugi::xml_node node) : RayTracePlot(node)
|
||||
{
|
||||
set_opaque_ids(node);
|
||||
set_diffuse_fraction(node);
|
||||
set_light_position(node);
|
||||
}
|
||||
|
||||
void SolidRayTracePlot::print_info() const
|
||||
{
|
||||
fmt::print("Plot Type: Solid ray-traced\n");
|
||||
RayTracePlot::print_info();
|
||||
}
|
||||
|
||||
void SolidRayTracePlot::create_output() const
|
||||
{
|
||||
size_t width = pixels_[0];
|
||||
size_t height = pixels_[1];
|
||||
ImageData data({width, height}, not_found_);
|
||||
|
||||
#pragma omp parallel for schedule(dynamic) collapse(2)
|
||||
for (int horiz = 0; horiz < pixels_[0]; ++horiz) {
|
||||
for (int vert = 0; vert < pixels_[1]; ++vert) {
|
||||
// RayTracePlot implements camera ray generation
|
||||
std::pair<Position, Direction> ru = get_pixel_ray(horiz, vert);
|
||||
PhongRay ray(ru.first, ru.second, *this);
|
||||
ray.trace();
|
||||
data(horiz, vert) = ray.result_color();
|
||||
}
|
||||
}
|
||||
|
||||
#ifdef USE_LIBPNG
|
||||
output_png(path_plot(), data);
|
||||
#else
|
||||
output_ppm(path_plot(), data);
|
||||
#endif
|
||||
}
|
||||
|
||||
void SolidRayTracePlot::set_opaque_ids(pugi::xml_node node)
|
||||
{
|
||||
if (check_for_node(node, "opaque_ids")) {
|
||||
auto opaque_ids_tmp = get_node_array<int>(node, "opaque_ids");
|
||||
|
||||
// It is read in as actual ID values, but we have to convert to indices in
|
||||
// mat/cell array
|
||||
for (auto& x : opaque_ids_tmp)
|
||||
x = color_by_ == PlotColorBy::mats ? model::material_map[x]
|
||||
: model::cell_map[x];
|
||||
|
||||
opaque_ids_.insert(opaque_ids_tmp.begin(), opaque_ids_tmp.end());
|
||||
}
|
||||
}
|
||||
|
||||
void SolidRayTracePlot::set_light_position(pugi::xml_node node)
|
||||
{
|
||||
if (check_for_node(node, "light_position")) {
|
||||
auto light_pos_tmp = get_node_array<double>(node, "light_position");
|
||||
|
||||
if (light_pos_tmp.size() != 3)
|
||||
fatal_error("Light position must be given as 3D coordinates");
|
||||
|
||||
light_location_.x = light_pos_tmp[0];
|
||||
light_location_.y = light_pos_tmp[1];
|
||||
light_location_.z = light_pos_tmp[2];
|
||||
} else {
|
||||
light_location_ = camera_position();
|
||||
}
|
||||
}
|
||||
|
||||
void SolidRayTracePlot::set_diffuse_fraction(pugi::xml_node node)
|
||||
{
|
||||
if (check_for_node(node, "diffuse_fraction")) {
|
||||
diffuse_fraction_ = std::stod(get_node_value(node, "diffuse_fraction"));
|
||||
if (diffuse_fraction_ < 0.0 || diffuse_fraction_ > 1.0) {
|
||||
fatal_error("Must have 0 <= diffuse fraction <= 1");
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
void Ray::compute_distance()
|
||||
{
|
||||
boundary() = distance_to_boundary(*this);
|
||||
}
|
||||
|
||||
void Ray::trace()
|
||||
{
|
||||
// To trace the ray from its origin all the way through the model, we have
|
||||
// to proceed in two phases. In the first, the ray may or may not be found
|
||||
// inside the model. If the ray is already in the model, phase one can be
|
||||
// skipped. Otherwise, the ray has to be advanced to the boundary of the
|
||||
// model where all the cells are defined. Importantly, this is assuming that
|
||||
// the model is convex, which is a very reasonable assumption for any
|
||||
// radiation transport model.
|
||||
//
|
||||
// After phase one is done, we can starting tracing from cell to cell within
|
||||
// the model. This step can use neighbor lists to accelerate the ray tracing.
|
||||
|
||||
// Attempt to initialize the particle. We may have to enter a loop to move
|
||||
// it up to the edge of the model.
|
||||
bool inside_cell = exhaustive_find_cell(*this, settings::verbosity >= 10);
|
||||
|
||||
// Advance to the boundary of the model
|
||||
while (!inside_cell) {
|
||||
advance_to_boundary_from_void();
|
||||
inside_cell = exhaustive_find_cell(*this, settings::verbosity >= 10);
|
||||
|
||||
// If true this means no surface was intersected. See cell.cpp and search
|
||||
// for numeric_limits to see where we return it.
|
||||
if (surface() == std::numeric_limits<int>::max()) {
|
||||
warning(fmt::format("Lost a ray, r = {}, u = {}", r(), u()));
|
||||
return;
|
||||
}
|
||||
|
||||
// Exit this loop and enter into cell-to-cell ray tracing (which uses
|
||||
// neighbor lists)
|
||||
if (inside_cell)
|
||||
break;
|
||||
|
||||
// if there is no intersection with the model, we're done
|
||||
if (boundary().surface == SURFACE_NONE)
|
||||
return;
|
||||
|
||||
event_counter_++;
|
||||
if (event_counter_ > MAX_INTERSECTIONS) {
|
||||
warning("Likely infinite loop in ray traced plot");
|
||||
return;
|
||||
}
|
||||
}
|
||||
|
||||
// Call the specialized logic for this type of ray. This is for the
|
||||
// intersection for the first intersection if we had one.
|
||||
if (boundary().surface != SURFACE_NONE) {
|
||||
// set the geometry state's surface attribute to be used for
|
||||
// surface normal computation
|
||||
surface() = boundary().surface;
|
||||
on_intersection();
|
||||
if (stop_)
|
||||
return;
|
||||
}
|
||||
|
||||
// reset surface attribute to zero after the first intersection so that it
|
||||
// doesn't perturb surface crossing logic from here on out
|
||||
surface() = 0;
|
||||
|
||||
// This is the ray tracing loop within the model. It exits after exiting
|
||||
// the model, which is equivalent to assuming that the model is convex.
|
||||
// It would be nice to factor out the on_intersection at the end of this
|
||||
// loop and then do "while (inside_cell)", but we can't guarantee it's
|
||||
// on a surface in that case. There might be some other way to set it
|
||||
// up that is perhaps a little more elegant, but this is what works just
|
||||
// fine.
|
||||
while (true) {
|
||||
|
||||
compute_distance();
|
||||
|
||||
// There are no more intersections to process
|
||||
// if we hit the edge of the model, so stop
|
||||
// the particle in that case. Also, just exit
|
||||
// if a negative distance was somehow computed.
|
||||
if (boundary().distance == INFTY || boundary().distance == INFINITY ||
|
||||
boundary().distance < 0) {
|
||||
return;
|
||||
}
|
||||
|
||||
// See below comment where call_on_intersection is checked in an
|
||||
// if statement for an explanation of this.
|
||||
bool call_on_intersection {true};
|
||||
if (boundary().distance < 10 * TINY_BIT) {
|
||||
call_on_intersection = false;
|
||||
}
|
||||
|
||||
// DAGMC surfaces expect us to go a little bit further than the advance
|
||||
// distance to properly check cell inclusion.
|
||||
boundary().distance += TINY_BIT;
|
||||
|
||||
// Advance particle, prepare for next intersection
|
||||
for (int lev = 0; lev < n_coord(); ++lev) {
|
||||
coord(lev).r += boundary().distance * coord(lev).u;
|
||||
}
|
||||
surface() = boundary().surface;
|
||||
n_coord_last() = n_coord();
|
||||
n_coord() = boundary().coord_level;
|
||||
if (boundary().lattice_translation[0] != 0 ||
|
||||
boundary().lattice_translation[1] != 0 ||
|
||||
boundary().lattice_translation[2] != 0) {
|
||||
cross_lattice(*this, boundary(), settings::verbosity >= 10);
|
||||
}
|
||||
|
||||
// Record how far the ray has traveled
|
||||
traversal_distance_ += boundary().distance;
|
||||
inside_cell = neighbor_list_find_cell(*this, settings::verbosity >= 10);
|
||||
|
||||
// Call the specialized logic for this type of ray. Note that we do not
|
||||
// call this if the advance distance is very small. Unfortunately, it seems
|
||||
// darn near impossible to get the particle advanced to the model boundary
|
||||
// and through it without sometimes accidentally calling on_intersection
|
||||
// twice. This incorrectly shades the region as occluded when it might not
|
||||
// actually be. By screening out intersection distances smaller than a
|
||||
// threshold 10x larger than the scoot distance used to advance up to the
|
||||
// model boundary, we can avoid that situation.
|
||||
if (call_on_intersection) {
|
||||
on_intersection();
|
||||
if (stop_)
|
||||
return;
|
||||
}
|
||||
|
||||
if (!inside_cell)
|
||||
return;
|
||||
|
||||
event_counter_++;
|
||||
if (event_counter_ > MAX_INTERSECTIONS) {
|
||||
warning("Likely infinite loop in ray traced plot");
|
||||
return;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
void ProjectionRay::on_intersection()
|
||||
{
|
||||
// This records a tuple with the following info
|
||||
//
|
||||
// 1) ID (material or cell depending on color_by_)
|
||||
// 2) Distance traveled by the ray through that ID
|
||||
// 3) Index of the intersected surface (starting from 1)
|
||||
|
||||
line_segments_.emplace_back(
|
||||
plot_.color_by_ == PlottableInterface::PlotColorBy::mats
|
||||
? material()
|
||||
: lowest_coord().cell,
|
||||
traversal_distance_, boundary().surface_index());
|
||||
}
|
||||
|
||||
void PhongRay::on_intersection()
|
||||
{
|
||||
// Check if we hit an opaque material or cell
|
||||
int hit_id = plot_.color_by_ == PlottableInterface::PlotColorBy::mats
|
||||
? material()
|
||||
: lowest_coord().cell;
|
||||
|
||||
// If we are reflected and have advanced beyond the camera,
|
||||
// the ray is done. This is checked here because we should
|
||||
// kill the ray even if the material is not opaque.
|
||||
if (reflected_ && (r() - plot_.camera_position()).dot(u()) >= 0.0) {
|
||||
stop();
|
||||
return;
|
||||
}
|
||||
|
||||
// Anything that's not opaque has zero impact on the plot.
|
||||
if (plot_.opaque_ids_.find(hit_id) == plot_.opaque_ids_.end())
|
||||
return;
|
||||
|
||||
if (!reflected_) {
|
||||
// reflect the particle and set the color to be colored by
|
||||
// the normal or the diffuse lighting contribution
|
||||
reflected_ = true;
|
||||
result_color_ = plot_.colors_[hit_id];
|
||||
Direction to_light = plot_.light_location_ - r();
|
||||
to_light /= to_light.norm();
|
||||
|
||||
// TODO
|
||||
// Not sure what can cause a surface token to be invalid here, although it
|
||||
// sometimes happens for a few pixels. It's very very rare, so proceed by
|
||||
// coloring the pixel with the overlap color. It seems to happen only for a
|
||||
// few pixels on the outer boundary of a hex lattice.
|
||||
//
|
||||
// We cannot detect it in the outer loop, and it only matters here, so
|
||||
// that's why the error handling is a little different than for a lost
|
||||
// ray.
|
||||
if (surface() == 0) {
|
||||
result_color_ = plot_.overlap_color_;
|
||||
stop();
|
||||
return;
|
||||
}
|
||||
|
||||
// Get surface pointer
|
||||
const auto& surf = model::surfaces.at(surface_index());
|
||||
|
||||
Direction normal = surf->normal(r_local());
|
||||
normal /= normal.norm();
|
||||
|
||||
// Need to apply translations to find the normal vector in
|
||||
// the base level universe's coordinate system.
|
||||
for (int lev = n_coord() - 2; lev >= 0; --lev) {
|
||||
if (coord(lev + 1).rotated) {
|
||||
const Cell& c {*model::cells[coord(lev).cell]};
|
||||
normal = normal.inverse_rotate(c.rotation_);
|
||||
}
|
||||
}
|
||||
|
||||
// use the normal opposed to the ray direction
|
||||
if (normal.dot(u()) > 0.0) {
|
||||
normal *= -1.0;
|
||||
}
|
||||
|
||||
// Facing away from the light means no lighting
|
||||
double dotprod = normal.dot(to_light);
|
||||
dotprod = std::max(0.0, dotprod);
|
||||
|
||||
double modulation =
|
||||
plot_.diffuse_fraction_ + (1.0 - plot_.diffuse_fraction_) * dotprod;
|
||||
result_color_ *= modulation;
|
||||
|
||||
// Now point the particle to the camera. We now begin
|
||||
// checking to see if it's occluded by another surface
|
||||
u() = to_light;
|
||||
|
||||
orig_hit_id_ = hit_id;
|
||||
|
||||
// OpenMC native CSG and DAGMC surfaces have some slight differences
|
||||
// in how they interpret particles that are sitting on a surface.
|
||||
// I don't know exactly why, but this makes everything work beautifully.
|
||||
if (surf->geom_type() == GeometryType::DAG) {
|
||||
surface() = 0;
|
||||
} else {
|
||||
surface() = -surface(); // go to other side
|
||||
}
|
||||
|
||||
// Must fully restart coordinate search. Why? Not sure.
|
||||
clear();
|
||||
|
||||
// Note this could likely be faster if we cached the previous
|
||||
// cell we were in before the reflection. This is the easiest
|
||||
// way to fully initialize all the sub-universe coordinates and
|
||||
// directions though.
|
||||
bool found = exhaustive_find_cell(*this);
|
||||
if (!found) {
|
||||
fatal_error("Lost particle after reflection.");
|
||||
}
|
||||
|
||||
// Must recalculate distance to boundary due to the
|
||||
// direction change
|
||||
compute_distance();
|
||||
|
||||
} else {
|
||||
// If it's not facing the light, we color with the diffuse contribution, so
|
||||
// next we check if we're going to occlude the last reflected surface. if
|
||||
// so, color by the diffuse contribution instead
|
||||
|
||||
if (orig_hit_id_ == -1)
|
||||
fatal_error("somehow a ray got reflected but not original ID set?");
|
||||
|
||||
result_color_ = plot_.colors_[orig_hit_id_];
|
||||
result_color_ *= plot_.diffuse_fraction_;
|
||||
stop();
|
||||
}
|
||||
}
|
||||
|
||||
extern "C" int openmc_id_map(const void* plot, int32_t* data_out)
|
||||
{
|
||||
|
||||
|
|
|
|||
Loading…
Add table
Add a link
Reference in a new issue