diff --git a/include/openmc/plot.h b/include/openmc/plot.h index 7e27679eabb..02ff2848ded 100644 --- a/include/openmc/plot.h +++ b/include/openmc/plot.h @@ -88,14 +88,21 @@ const RGBColor BLACK {0, 0, 0}; * is designed to be implemented by classes that produce plot-relevant data * which can be visualized. */ + +typedef xt::xtensor ImageData; class PlottableInterface { +public: + PlottableInterface() = default; + + void set_default_colors(); + private: void set_id(pugi::xml_node plot_node); int id_; // unique plot ID void set_bg_color(pugi::xml_node plot_node); void set_universe(pugi::xml_node plot_node); - void set_default_colors(pugi::xml_node plot_node); + void set_color_by(pugi::xml_node plot_node); void set_user_colors(pugi::xml_node plot_node); void set_overlap_color(pugi::xml_node plot_node); void set_mask(pugi::xml_node plot_node); @@ -107,9 +114,15 @@ class PlottableInterface { public: enum class PlotColorBy { cells = 0, mats = 1 }; + // Generates image data based on plot parameters and returns it + virtual ImageData create_image() const = 0; + // Creates the output image named path_plot_ virtual void create_output() const = 0; + // Write populated image data to file + void write_image(const ImageData& data) const; + // Print useful info to the terminal virtual void print_info() const = 0; @@ -117,20 +130,19 @@ class PlottableInterface { std::string& path_plot() { return path_plot_; } int id() const { return id_; } int level() const { return level_; } + PlotColorBy color_by() const { return color_by_; } // Public color-related data PlottableInterface(pugi::xml_node plot_node); virtual ~PlottableInterface() = default; - int level_; // Universe level to plot - bool color_overlaps_; // Show overlapping cells? + int level_ {-1}; // Universe level to plot + bool color_overlaps_ {false}; // Show overlapping cells? PlotColorBy color_by_; // Plot coloring (cell/material) RGBColor not_found_ {WHITE}; // Plot background color RGBColor overlap_color_ {RED}; // Plot overlap color vector colors_; // Plot colors }; -typedef xt::xtensor ImageData; - struct IdData { // Constructor IdData(size_t h_res, size_t v_res); @@ -166,6 +178,11 @@ class SlicePlotBase { enum class PlotBasis { xy = 1, xz = 2, yz = 3 }; + // Accessors + + const std::array& pixels() const { return pixels_; } + std::array& pixels() { return pixels_; } + // Members public: Position origin_; //!< Plot origin in geometry @@ -270,11 +287,11 @@ class Plot : public PlottableInterface, public SlicePlotBase { public: // Add mesh lines to ImageData void draw_mesh_lines(ImageData& data) const; - void create_image() const; + ImageData create_image() const override; void create_voxel() const; - virtual void create_output() const; - virtual void print_info() const; + void create_output() const override; + void print_info() const override; PlotType type_; //!< Plot type (Slice/Voxel) int meshlines_width_; //!< Width of lines added to the plot @@ -294,17 +311,32 @@ class Plot : public PlottableInterface, public SlicePlotBase { */ class RayTracePlot : public PlottableInterface { public: + RayTracePlot() = default; RayTracePlot(pugi::xml_node plot); // Standard getters. No setting since it's done from XML. const Position& camera_position() const { return camera_position_; } + Position& camera_position() { return camera_position_; } const Position& look_at() const { return look_at_; } + Position& look_at() { return look_at_; } + const double& horizontal_field_of_view() const { return horizontal_field_of_view_; } + double& horizontal_field_of_view() { return horizontal_field_of_view_; } + + void print_info() const override; + + const std::array& pixels() const { return pixels_; } + std::array& pixels() { return pixels_; } - virtual void print_info() const; + const Direction& up() const { return up_; } + Direction& up() { return up_; } + + //! brief Updates the cached camera-to-model matrix after changes to + //! camera parameters. + void update_view(); protected: Direction camera_x_axis() const @@ -330,8 +362,6 @@ class RayTracePlot : public PlottableInterface { */ std::pair get_pixel_ray(int horiz, int vert) const; - std::array pixels_; // pixel dimension of resulting image - private: void set_look_at(pugi::xml_node node); void set_camera_position(pugi::xml_node node); @@ -341,8 +371,8 @@ class RayTracePlot : public PlottableInterface { double horizontal_field_of_view_ {70.0}; // horiz. f.o.v. in degrees Position camera_position_; // where camera is - Position look_at_; // point camera is centered looking at - + Position look_at_; // point camera is centered looking at + std::array pixels_; // pixel dimension of resulting image Direction up_ {0.0, 0.0, 1.0}; // which way is up /* The horizontal thickness, if using an orthographic projection. @@ -377,8 +407,9 @@ class WireframeRayTracePlot : public RayTracePlot { public: WireframeRayTracePlot(pugi::xml_node plot); - virtual void create_output() const; - virtual void print_info() const; + ImageData create_image() const override; + void create_output() const override; + void print_info() const override; private: void set_opacities(pugi::xml_node node); @@ -434,10 +465,22 @@ class SolidRayTracePlot : public RayTracePlot { friend class PhongRay; public: + SolidRayTracePlot() = default; + SolidRayTracePlot(pugi::xml_node plot); - virtual void create_output() const; - virtual void print_info() const; + ImageData create_image() const override; + void create_output() const override; + void print_info() const override; + + const std::unordered_set& opaque_ids() const { return opaque_ids_; } + std::unordered_set& opaque_ids() { return opaque_ids_; } + + const Position& light_location() const { return light_location_; } + Position& light_location() { return light_location_; } + + const double& diffuse_fraction() const { return diffuse_fraction_; } + double& diffuse_fraction() { return diffuse_fraction_; } private: void set_opaque_ids(pugi::xml_node node); @@ -496,7 +539,7 @@ class ProjectionRay : public Ray { : Ray(r, u), plot_(plot), line_segments_(line_segments) {} - virtual void on_intersection() override; + void on_intersection() override; private: /* Store a reference to the plot object which is running this ray, in order @@ -519,7 +562,7 @@ class PhongRay : public Ray { result_color_ = plot_.not_found_; } - virtual void on_intersection() override; + void on_intersection() override; const RGBColor& result_color() { return result_color_; } diff --git a/src/plot.cpp b/src/plot.cpp index 2cadc48cefe..daeb0bde7c6 100644 --- a/src/plot.cpp +++ b/src/plot.cpp @@ -15,6 +15,7 @@ #include #endif +#include "openmc/cell.h" #include "openmc/constants.h" #include "openmc/container_util.h" #include "openmc/dagmc.h" @@ -123,11 +124,21 @@ extern "C" int openmc_plot_geometry() return 0; } +void PlottableInterface::write_image(const ImageData& data) const +{ +#ifdef USE_LIBPNG + output_png(path_plot(), data); +#else + output_ppm(path_plot(), data); +#endif +} + void Plot::create_output() const { if (PlotType::slice == type_) { // create 2D image - create_image(); + ImageData image = create_image(); + write_image(image); } else if (PlotType::voxel == type_) { // create voxel file for 3D viewing create_voxel(); @@ -170,9 +181,9 @@ void Plot::print_info() const fmt::print("Basis: YZ\n"); break; } - fmt::print("Pixels: {} {}\n", pixels_[0], pixels_[1]); + fmt::print("Pixels: {} {}\n", pixels()[0], pixels()[1]); } else if (PlotType::voxel == type_) { - fmt::print("Voxels: {} {} {}\n", pixels_[0], pixels_[1], pixels_[2]); + fmt::print("Voxels: {} {} {}\n", pixels()[0], pixels()[1], pixels()[2]); } } @@ -234,11 +245,10 @@ void free_memory_plot() // creates an image based on user input from a plots.xml // specification in the PNG/PPM format -void Plot::create_image() const +ImageData Plot::create_image() const { - - size_t width = pixels_[0]; - size_t height = pixels_[1]; + size_t width = pixels()[0]; + size_t height = pixels()[1]; ImageData data({width, height}, not_found_); @@ -275,12 +285,7 @@ void Plot::create_image() const draw_mesh_lines(data); } -// create image file -#ifdef USE_LIBPNG - output_png(path_plot(), data); -#else - output_ppm(path_plot(), data); -#endif + return data; } void PlottableInterface::set_id(pugi::xml_node plot_node) @@ -348,17 +353,17 @@ void Plot::set_output_path(pugi::xml_node plot_node) vector pxls = get_node_array(plot_node, "pixels"); if (PlotType::slice == type_) { if (pxls.size() == 2) { - pixels_[0] = pxls[0]; - pixels_[1] = pxls[1]; + pixels()[0] = pxls[0]; + pixels()[1] = pxls[1]; } else { fatal_error( fmt::format(" must be length 2 in slice plot {}", id())); } } else if (PlotType::voxel == type_) { if (pxls.size() == 3) { - pixels_[0] = pxls[0]; - pixels_[1] = pxls[1]; - pixels_[2] = pxls[2]; + pixels()[0] = pxls[0]; + pixels()[1] = pxls[1]; + pixels()[2] = pxls[2]; } else { fatal_error( fmt::format(" must be length 3 in voxel plot {}", id())); @@ -447,23 +452,31 @@ void PlottableInterface::set_universe(pugi::xml_node plot_node) } } -void PlottableInterface::set_default_colors(pugi::xml_node plot_node) +void PlottableInterface::set_color_by(pugi::xml_node plot_node) { - // Copy plot color type and initialize all colors randomly + // Copy plot color type std::string pl_color_by = "cell"; if (check_for_node(plot_node, "color_by")) { pl_color_by = get_node_value(plot_node, "color_by", true); } if ("cell" == pl_color_by) { color_by_ = PlotColorBy::cells; - colors_.resize(model::cells.size()); } else if ("material" == pl_color_by) { color_by_ = PlotColorBy::mats; - colors_.resize(model::materials.size()); } else { fatal_error(fmt::format( "Unsupported plot color type '{}' in plot {}", pl_color_by, id())); } +} + +void PlottableInterface::set_default_colors() +{ + // Copy plot color type and initialize all colors randomly + if (PlotColorBy::cells == color_by_) { + colors_.resize(model::cells.size()); + } else if (PlotColorBy::mats == color_by_) { + colors_.resize(model::materials.size()); + } for (auto& c : colors_) { c = random_color(); @@ -710,7 +723,8 @@ PlottableInterface::PlottableInterface(pugi::xml_node plot_node) set_id(plot_node); set_bg_color(plot_node); set_universe(plot_node); - set_default_colors(plot_node); + set_color_by(plot_node); + set_default_colors(); set_user_colors(plot_node); set_mask(plot_node); set_overlap_color(plot_node); @@ -857,27 +871,27 @@ void Plot::draw_mesh_lines(ImageData& data) const int ax2_min, ax2_max; if (axis_lines.second.size() > 0) { double frac = (axis_lines.second.back() - ll_plot[ax2]) / width[ax2]; - ax2_min = (1.0 - frac) * pixels_[1]; + ax2_min = (1.0 - frac) * pixels()[1]; if (ax2_min < 0) ax2_min = 0; frac = (axis_lines.second.front() - ll_plot[ax2]) / width[ax2]; - ax2_max = (1.0 - frac) * pixels_[1]; - if (ax2_max > pixels_[1]) - ax2_max = pixels_[1]; + ax2_max = (1.0 - frac) * pixels()[1]; + if (ax2_max > pixels()[1]) + ax2_max = pixels()[1]; } else { ax2_min = 0; - ax2_max = pixels_[1]; + ax2_max = pixels()[1]; } // Iterate across the first axis and draw lines. for (auto ax1_val : axis_lines.first) { double frac = (ax1_val - ll_plot[ax1]) / width[ax1]; - int ax1_ind = frac * pixels_[0]; + int ax1_ind = frac * pixels()[0]; for (int ax2_ind = ax2_min; ax2_ind < ax2_max; ++ax2_ind) { for (int plus = 0; plus <= meshlines_width_; plus++) { - if (ax1_ind + plus >= 0 && ax1_ind + plus < pixels_[0]) + if (ax1_ind + plus >= 0 && ax1_ind + plus < pixels()[0]) data(ax1_ind + plus, ax2_ind) = rgb; - if (ax1_ind - plus >= 0 && ax1_ind - plus < pixels_[0]) + if (ax1_ind - plus >= 0 && ax1_ind - plus < pixels()[0]) data(ax1_ind - plus, ax2_ind) = rgb; } } @@ -887,27 +901,27 @@ void Plot::draw_mesh_lines(ImageData& data) const int ax1_min, ax1_max; if (axis_lines.first.size() > 0) { double frac = (axis_lines.first.front() - ll_plot[ax1]) / width[ax1]; - ax1_min = frac * pixels_[0]; + ax1_min = frac * pixels()[0]; if (ax1_min < 0) ax1_min = 0; frac = (axis_lines.first.back() - ll_plot[ax1]) / width[ax1]; - ax1_max = frac * pixels_[0]; - if (ax1_max > pixels_[0]) - ax1_max = pixels_[0]; + ax1_max = frac * pixels()[0]; + if (ax1_max > pixels()[0]) + ax1_max = pixels()[0]; } else { ax1_min = 0; - ax1_max = pixels_[0]; + ax1_max = pixels()[0]; } // Iterate across the second axis and draw lines. for (auto ax2_val : axis_lines.second) { double frac = (ax2_val - ll_plot[ax2]) / width[ax2]; - int ax2_ind = (1.0 - frac) * pixels_[1]; + int ax2_ind = (1.0 - frac) * pixels()[1]; for (int ax1_ind = ax1_min; ax1_ind < ax1_max; ++ax1_ind) { for (int plus = 0; plus <= meshlines_width_; plus++) { - if (ax2_ind + plus >= 0 && ax2_ind + plus < pixels_[1]) + if (ax2_ind + plus >= 0 && ax2_ind + plus < pixels()[1]) data(ax1_ind, ax2_ind + plus) = rgb; - if (ax2_ind - plus >= 0 && ax2_ind - plus < pixels_[1]) + if (ax2_ind - plus >= 0 && ax2_ind - plus < pixels()[1]) data(ax1_ind, ax2_ind - plus) = rgb; } } @@ -928,9 +942,9 @@ void Plot::create_voxel() const { // compute voxel widths in each direction array vox; - vox[0] = width_[0] / static_cast(pixels_[0]); - vox[1] = width_[1] / static_cast(pixels_[1]); - vox[2] = width_[2] / static_cast(pixels_[2]); + vox[0] = width_[0] / static_cast(pixels()[0]); + vox[1] = width_[1] / static_cast(pixels()[1]); + vox[2] = width_[2] / static_cast(pixels()[2]); // initial particle position Position ll = origin_ - width_ / 2.; @@ -952,18 +966,18 @@ void Plot::create_voxel() const // Write current date and time write_attribute(file_id, "date_and_time", time_stamp().c_str()); - array pixels; - std::copy(pixels_.begin(), pixels_.end(), pixels.begin()); - write_attribute(file_id, "num_voxels", pixels); + array h5_pixels; + std::copy(pixels().begin(), pixels().end(), h5_pixels.begin()); + write_attribute(file_id, "num_voxels", h5_pixels); write_attribute(file_id, "voxel_width", vox); write_attribute(file_id, "lower_left", ll); // Create dataset for voxel data -- note that the dimensions are reversed // since we want the order in the file to be z, y, x hsize_t dims[3]; - dims[0] = pixels_[2]; - dims[1] = pixels_[1]; - dims[2] = pixels_[0]; + dims[0] = pixels()[2]; + dims[1] = pixels()[1]; + dims[2] = pixels()[0]; hid_t dspace, dset, memspace; voxel_init(file_id, &(dims[0]), &dspace, &dset, &memspace); @@ -971,11 +985,11 @@ void Plot::create_voxel() const pltbase.width_ = width_; pltbase.origin_ = origin_; pltbase.basis_ = PlotBasis::xy; - pltbase.pixels_ = pixels_; + pltbase.pixels() = pixels(); pltbase.slice_color_overlaps_ = color_overlaps_; ProgressBar pb; - for (int z = 0; z < pixels_[2]; z++) { + for (int z = 0; z < pixels()[2]; z++) { // update z coordinate pltbase.origin_.z = ll.z + z * vox[2]; @@ -993,7 +1007,7 @@ void Plot::create_voxel() const // update progress bar pb.set_value( - 100. * static_cast(z + 1) / static_cast((pixels_[2]))); + 100. * static_cast(z + 1) / static_cast((pixels()[2]))); } voxel_finalize(dspace, dset, memspace); @@ -1052,7 +1066,10 @@ RayTracePlot::RayTracePlot(pugi::xml_node node) : PlottableInterface(node) check_for_node(node, "field_of_view")) fatal_error("orthographic_width and field_of_view are mutually exclusive " "parameters."); +} +void RayTracePlot::update_view() +{ // 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(); @@ -1079,6 +1096,7 @@ WireframeRayTracePlot::WireframeRayTracePlot(pugi::xml_node node) set_wireframe_thickness(node); set_wireframe_ids(node); set_wireframe_color(node); + update_view(); } void WireframeRayTracePlot::set_wireframe_color(pugi::xml_node plot_node) @@ -1182,8 +1200,8 @@ std::pair RayTracePlot::get_pixel_ray( // 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(pixels_[0]); - double p1 = static_cast(pixels_[1]); + double p0 = static_cast(pixels()[0]); + double p1 = static_cast(pixels()[1]); double vert_fov_radians = horiz_fov_radians * p1 / p0; // focal_plane_dist can be changed to alter the perspective distortion @@ -1219,10 +1237,10 @@ std::pair RayTracePlot::get_pixel_ray( return result; } -void WireframeRayTracePlot::create_output() const +ImageData WireframeRayTracePlot::create_image() const { - size_t width = pixels_[0]; - size_t height = pixels_[1]; + 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 @@ -1245,11 +1263,11 @@ void WireframeRayTracePlot::create_output() const std::vector>> this_line_segments( n_threads); for (int t = 0; t < n_threads; ++t) { - this_line_segments[t].resize(pixels_[0]); + this_line_segments[t].resize(pixels()[0]); } // The last thread writes to this, and the first thread reads from it. - std::vector> old_segments(pixels_[0]); + std::vector> old_segments(pixels()[0]); #pragma omp parallel { @@ -1257,7 +1275,7 @@ void WireframeRayTracePlot::create_output() const const int tid = thread_num(); int vert = tid; - for (int iter = 0; iter <= pixels_[1] / n_threads; iter++) { + for (int iter = 0; iter <= pixels()[1] / n_threads; iter++) { // 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 @@ -1266,9 +1284,9 @@ void WireframeRayTracePlot::create_output() const if (tid == n_threads - 1) old_segments = this_line_segments[n_threads - 1]; - if (vert < pixels_[1]) { + if (vert < pixels()[1]) { - for (int horiz = 0; horiz < pixels_[0]; ++horiz) { + for (int horiz = 0; horiz < pixels()[0]; ++horiz) { // RayTracePlot implements camera ray generation std::pair ru = get_pixel_ray(horiz, vert); @@ -1330,7 +1348,7 @@ void WireframeRayTracePlot::create_output() const // Now that the horizontal line has finished rendering, we can fill in // wireframe entries that require comparison among all the threads. Hence // the omp barrier being used. It has to be OUTSIDE any if blocks! - if (vert < pixels_[1]) { + if (vert < pixels()[1]) { // Loop over horizontal pixels, checking intersection stack of upper // neighbor @@ -1340,7 +1358,7 @@ void WireframeRayTracePlot::create_output() const else top_cmp = &this_line_segments[tid - 1]; - for (int horiz = 0; horiz < pixels_[0]; ++horiz) { + for (int horiz = 0; horiz < pixels()[0]; ++horiz) { if (!trackstack_equivalent( this_line_segments[tid][horiz], (*top_cmp)[horiz])) { wireframe_initial(horiz, vert) = 1; @@ -1357,8 +1375,8 @@ void WireframeRayTracePlot::create_output() const } // end omp parallel // Now thicken the wireframe lines and apply them to our image - for (int vert = 0; vert < pixels_[1]; ++vert) { - for (int horiz = 0; horiz < pixels_[0]; ++horiz) { + for (int vert = 0; vert < pixels()[1]; ++vert) { + for (int horiz = 0; horiz < pixels()[0]; ++horiz) { if (wireframe_initial(horiz, vert)) { if (wireframe_thickness_ == 1) data(horiz, vert) = wireframe_color_; @@ -1369,19 +1387,21 @@ void WireframeRayTracePlot::create_output() const if (i * i + j * j < wireframe_thickness_ * wireframe_thickness_) { // Check if wireframe pixel is out of bounds - int w_i = std::max(std::min(horiz + i, pixels_[0] - 1), 0); - int w_j = std::max(std::min(vert + j, pixels_[1] - 1), 0); + int w_i = std::max(std::min(horiz + i, pixels()[0] - 1), 0); + int w_j = std::max(std::min(vert + j, pixels()[1] - 1), 0); data(w_i, w_j) = wireframe_color_; } } } } -#ifdef USE_LIBPNG - output_png(path_plot(), data); -#else - output_ppm(path_plot(), data); -#endif + return data; +} + +void WireframeRayTracePlot::create_output() const +{ + ImageData data = create_image(); + write_image(data); } void RayTracePlot::print_info() const @@ -1391,7 +1411,7 @@ void RayTracePlot::print_info() const fmt::print("Look at: {} {} {}\n", look_at_.x, look_at_.y, look_at_.z); fmt::print( "Horizontal field of view: {} degrees\n", horizontal_field_of_view_); - fmt::print("Pixels: {} {}\n", pixels_[0], pixels_[1]); + fmt::print("Pixels: {} {}\n", pixels()[0], pixels()[1]); } void WireframeRayTracePlot::print_info() const @@ -1473,8 +1493,8 @@ void RayTracePlot::set_pixels(pugi::xml_node node) if (pxls.size() != 2) fatal_error( fmt::format(" must be length 2 in projection plot {}", id())); - pixels_[0] = pxls[0]; - pixels_[1] = pxls[1]; + pixels()[0] = pxls[0]; + pixels()[1] = pxls[1]; } void RayTracePlot::set_camera_position(pugi::xml_node node) @@ -1521,6 +1541,7 @@ SolidRayTracePlot::SolidRayTracePlot(pugi::xml_node node) : RayTracePlot(node) set_opaque_ids(node); set_diffuse_fraction(node); set_light_position(node); + update_view(); } void SolidRayTracePlot::print_info() const @@ -1529,15 +1550,15 @@ void SolidRayTracePlot::print_info() const RayTracePlot::print_info(); } -void SolidRayTracePlot::create_output() const +ImageData SolidRayTracePlot::create_image() const { - size_t width = pixels_[0]; - size_t height = pixels_[1]; + 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) { + for (int horiz = 0; horiz < pixels()[0]; ++horiz) { + for (int vert = 0; vert < pixels()[1]; ++vert) { // RayTracePlot implements camera ray generation std::pair ru = get_pixel_ray(horiz, vert); PhongRay ray(ru.first, ru.second, *this); @@ -1546,11 +1567,13 @@ void SolidRayTracePlot::create_output() const } } -#ifdef USE_LIBPNG - output_png(path_plot(), data); -#else - output_ppm(path_plot(), data); -#endif + return data; +} + +void SolidRayTracePlot::create_output() const +{ + ImageData data = create_image(); + write_image(data); } void SolidRayTracePlot::set_opaque_ids(pugi::xml_node node)