Skip to content

Commit 76cc722

Browse files
gridleygridley
committed
add option to do orthographic projection
1 parent d51fd87 commit 76cc722

3 files changed

Lines changed: 65 additions & 12 deletions

File tree

‎include/openmc/plot.h‎

Lines changed: 10 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -274,9 +274,6 @@ class Plot : public PlottableInterface, public SlicePlotBase {
274274

275275
class ProjectionPlot : public PlottableInterface {
276276

277-
// TODO add optional choice between perspective or orthographic
278-
// TODO add optionl to turn off wireframe
279-
280277
public:
281278
ProjectionPlot(pugi::xml_node plot);
282279

@@ -290,6 +287,7 @@ class ProjectionPlot : public PlottableInterface {
290287
void set_field_of_view(pugi::xml_node node);
291288
void set_pixels(pugi::xml_node node);
292289
void set_opacities(pugi::xml_node node);
290+
void set_orthographic_width(pugi::xml_node node);
293291

294292
/* Used for drawing wireframe and colors. We record the list of
295293
* surface/cell/material intersections and the corresponding lengths as a ray
@@ -316,6 +314,12 @@ class ProjectionPlot : public PlottableInterface {
316314
Direction up_ {0.0, 0.0, 1.0}; // which way is up
317315
bool wireframe_; // draw wireframe around ID boundaries (material or cell
318316
// based on color_by)
317+
318+
/* The horizontal thickness, if using an orthographic projection.
319+
* If set to zero, we assume using a perspective projection.
320+
*/
321+
double orthographic_width_ {0.0};
322+
319323
RGBColor wireframe_color_ {BLACK}; // wireframe color
320324
std::vector<double>
321325
xs_; // macro cross section values for cell volume rendering
@@ -332,6 +336,9 @@ class ProjectionPlot : public PlottableInterface {
332336
*/
333337
static bool trackstack_equivalent(const std::vector<TrackSegment>& track1,
334338
const std::vector<TrackSegment>& track2);
339+
340+
// Closed form 3x3 matrix inversion
341+
static std::vector<double> invert_matrix(const std::vector<double>& input);
335342
};
336343

337344
//===============================================================================

‎include/openmc/position.h‎

Lines changed: 5 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -83,6 +83,11 @@ struct Position {
8383
return x * other.x + y * other.y + z * other.z;
8484
}
8585
inline double norm() const { return std::sqrt(x * x + y * y + z * z); }
86+
inline Position cross(Position other) const
87+
{
88+
return {y * other.z - z * other.y, z * other.x - x * other.z,
89+
x * other.y - y * other.x};
90+
}
8691

8792
//! Reflect a direction across a normal vector
8893
//! \param[in] other Vector to reflect across

‎src/plot.cpp‎

Lines changed: 50 additions & 9 deletions
Original file line numberDiff line numberDiff line change
@@ -1033,6 +1033,12 @@ ProjectionPlot::ProjectionPlot(pugi::xml_node node) : PlottableInterface(node)
10331033
set_field_of_view(node);
10341034
set_pixels(node);
10351035
set_opacities(node);
1036+
set_orthographic_width(node);
1037+
1038+
if (check_for_node(node, "orthographic_width") &&
1039+
check_for_node(node, "field_of_view"))
1040+
fatal_error("orthographic_width and field_of_view are mutually exclusive "
1041+
"parameters.");
10361042
}
10371043

10381044
void ProjectionPlot::set_output_path(pugi::xml_node node)
@@ -1101,15 +1107,26 @@ void ProjectionPlot::create_output() const
11011107
{
11021108
// Get centerline vector for camera-to-model. We create vectors around this
11031109
// that form a pixel array, and then trace rays along that.
1110+
auto up = up_ / up_.norm();
11041111
Direction looking_direction = look_at_ - camera_position_;
11051112
looking_direction /= looking_direction.norm();
1113+
if (std::abs(std::abs(looking_direction.dot(up)) - 1.0) < 1e-9)
1114+
fatal_error("Up vector cannot align with vector between camera position "
1115+
"and look_at!");
1116+
Direction cam_yaxis = looking_direction.cross(up);
1117+
cam_yaxis /= cam_yaxis.norm();
1118+
Direction cam_zaxis = cam_yaxis.cross(looking_direction);
1119+
cam_zaxis /= cam_zaxis.norm();
1120+
1121+
// Transformation matrix for directions
1122+
std::vector<double> camera_to_model = {looking_direction.x, cam_yaxis.x,
1123+
cam_zaxis.x, looking_direction.y, cam_yaxis.y, cam_zaxis.y,
1124+
looking_direction.z, cam_yaxis.z, cam_zaxis.z};
11061125

11071126
// Now we convert to the polar coordinate system with the polar angle
11081127
// measuring the angle from the vector up_. Phi is the rotation about up_. For
11091128
// now, up_ is hard-coded to be +z.
11101129
constexpr double DEGREE_TO_RADIAN = M_PI / 180.0;
1111-
double phi0 = std::atan2(looking_direction.y, looking_direction.x);
1112-
double mu0 = std::acos(looking_direction.z);
11131130
double horiz_fov_radians = horizontal_field_of_view_ * DEGREE_TO_RADIAN;
11141131
double p0 = static_cast<double>(pixels_[0]);
11151132
double p1 = static_cast<double>(pixels_[1]);
@@ -1154,11 +1171,25 @@ void ProjectionPlot::create_output() const
11541171
for (int vert = 0; vert < pixels_[1]; ++vert) {
11551172
old_segments = this_line_segments;
11561173
for (int horiz = 0; horiz < pixels_[0]; ++horiz) {
1157-
double this_phi = phi0 - horiz_fov_radians / 2.0 + dphi * horiz;
1158-
double this_mu = mu0 - vert_fov_radians / 2.0 + dmu * vert;
1159-
s.u.x = std::cos(this_phi) * std::sin(this_mu);
1160-
s.u.y = std::sin(this_phi) * std::sin(this_mu);
1161-
s.u.z = std::cos(this_mu);
1174+
1175+
// Generate the starting position/direction of the ray
1176+
if (orthographic_width_ == 0.0) { // perspective projection
1177+
double this_phi = -horiz_fov_radians / 2.0 + dphi * horiz;
1178+
double this_mu = -vert_fov_radians / 2.0 + dmu * vert + M_PI / 2.0;
1179+
Direction camera_local_vec;
1180+
camera_local_vec.x = std::cos(this_phi) * std::sin(this_mu);
1181+
camera_local_vec.y = std::sin(this_phi) * std::sin(this_mu);
1182+
camera_local_vec.z = std::cos(this_mu);
1183+
s.u = camera_local_vec.rotate(camera_to_model);
1184+
} else { // orthographic projection
1185+
s.u = looking_direction;
1186+
1187+
double x_pix_coord = (static_cast<double>(horiz) - p0 / 2.0) / p0;
1188+
double y_pix_coord = (static_cast<double>(vert) - p1 / 2.0) / p0;
1189+
s.r = camera_position_ + cam_yaxis * x_pix_coord * orthographic_width_ +
1190+
cam_zaxis * y_pix_coord * orthographic_width_;
1191+
}
1192+
11621193
p.from_source(&s); // put particle at camera
11631194
bool hitsomething = false;
11641195
bool intersection_found = true;
@@ -1238,7 +1269,7 @@ void ProjectionPlot::create_output() const
12381269
draw_wireframe || !trackstack_equivalent(
12391270
this_line_segments[horiz], old_segments[horiz]);
12401271
}
1241-
if (draw_wireframe) {
1272+
if (draw_wireframe && wireframe_) {
12421273
data(horiz, vert) = wireframe_color_;
12431274
}
12441275
}
@@ -1263,7 +1294,6 @@ void ProjectionPlot::print_info() const
12631294

12641295
void ProjectionPlot::set_opacities(pugi::xml_node node)
12651296
{
1266-
12671297
xs_.resize(colors_.size(), 1e6); // set to large value for opaque by default
12681298

12691299
for (auto cn : node.children("color")) {
@@ -1292,6 +1322,17 @@ void ProjectionPlot::set_opacities(pugi::xml_node node)
12921322
}
12931323
}
12941324

1325+
void ProjectionPlot::set_orthographic_width(pugi::xml_node node)
1326+
{
1327+
if (check_for_node(node, "orthographic_width")) {
1328+
double orthographic_width =
1329+
std::stod(get_node_value(node, "orthographic_width", true));
1330+
if (orthographic_width < 0.0)
1331+
fatal_error("Requires positive orthographic_width");
1332+
orthographic_width_ = orthographic_width;
1333+
}
1334+
}
1335+
12951336
void ProjectionPlot::set_pixels(pugi::xml_node node)
12961337
{
12971338
vector<int> pxls = get_node_array<int>(node, "pixels");

0 commit comments

Comments
 (0)