diff --git a/coretrace/simulation_data/aperture.cpp b/coretrace/simulation_data/aperture.cpp index 5982dd51..a107f1c5 100644 --- a/coretrace/simulation_data/aperture.cpp +++ b/coretrace/simulation_data/aperture.cpp @@ -670,6 +670,7 @@ namespace SolTrace::Data x3(x3), y3(y3), x4(x4), y4(y4) { + ensure_valid_diagonal(); } IrregularQuadrilateral::IrregularQuadrilateral(const nlohmann::ordered_json &jnode) @@ -683,6 +684,7 @@ namespace SolTrace::Data this->y3 = jnode.at("y3"); this->x4 = jnode.at("x4"); this->y4 = jnode.at("y4"); + ensure_valid_diagonal(); } double IrregularQuadrilateral::aperture_area() const @@ -750,6 +752,54 @@ namespace SolTrace::Data return inquad(x1, y1, x2, y2, x3, y3, x4, y4, x, y); } + void IrregularQuadrilateral::ensure_valid_diagonal() + { + // Returns true when x2 and x4 are on opposite sides of the x1-x3 line, + // i.e. the x1-x3 diagonal is interior to the polygon. + const auto diagonal_ok = [this]() + { + const double ex = x3 - x1, ey = y3 - y1; + const double d2 = ex * (y2 - y1) - ey * (x2 - x1); + const double d4 = ex * (y4 - y1) - ey * (x4 - x1); + return d2 * d4 < 0.0; + }; + + if (!diagonal_ok()) + { + // Vertices are not in simple-polygon order (self-intersecting / + // bowtie input). Re-sort by angle from the centroid to recover a + // valid CCW traversal. + const double cx = 0.25 * (x1 + x2 + x3 + x4); + const double cy = 0.25 * (y1 + y2 + y3 + y4); + using P2 = std::array; + std::array pts = {P2{x1, y1}, P2{x2, y2}, P2{x3, y3}, P2{x4, y4}}; + std::sort(pts.begin(), pts.end(), + [cx, cy](const P2 &a, const P2 &b) + { + return std::atan2(a[1] - cy, a[0] - cx) < + std::atan2(b[1] - cy, b[0] - cx); + }); + x1 = pts[0][0]; y1 = pts[0][1]; + x2 = pts[1][0]; y2 = pts[1][1]; + x3 = pts[2][0]; y3 = pts[2][1]; + x4 = pts[3][0]; y4 = pts[3][1]; + + // The sort produces a valid CCW polygon but the starting vertex is + // arbitrary. If the reflex vertex landed at position x1 or x3, the + // x1-x3 diagonal is still exterior. One cyclic shift swaps which + // pair of opposite vertices forms the diagonal (x1-x3 becomes the + // old x2-x4), which is guaranteed to be interior for a simple polygon. + if (!diagonal_ok()) + { + const double tx = x1, ty = y1; + x1 = x2; y1 = y2; + x2 = x3; y2 = y3; + x3 = x4; y3 = y4; + x4 = tx; y4 = ty; + } + } + } + aperture_ptr IrregularQuadrilateral::make_copy() const { // Invokes the implicit copy constructor diff --git a/coretrace/simulation_data/aperture.hpp b/coretrace/simulation_data/aperture.hpp index c460234f..bef66f4c 100644 --- a/coretrace/simulation_data/aperture.hpp +++ b/coretrace/simulation_data/aperture.hpp @@ -709,7 +709,14 @@ namespace SolTrace::Data struct IrregularQuadrilateral : public Aperture { - // Locations of the 4 vertices + // Locations of the 4 vertices. + // + // The quadrilateral must be a SIMPLE (non-self-intersecting) polygon. + // The vertices may be supplied in any order; the constructors call + // ensure_valid_diagonal() which sorts them by angle from the centroid + // into a consistent CCW traversal and ensures the x1-x3 diagonal is + // interior. All triangle-decomposition logic (inquad, the flat and + // parabolic intersection kernels) depends on this diagonal being interior. double x1; double y1; double x2; @@ -786,6 +793,16 @@ namespace SolTrace::Data */ virtual std::tuple, std::vector> triangulation() const override; + + /** + * @brief Ensure the x1-x3 diagonal is interior to the quadrilateral. + * All runners decompose the quad into triangles (x1,x2,x3) and + * (x1,x3,x4) along this diagonal. If x2 and x4 lie on the same side + * of the x1-x3 line the diagonal is exterior and the decomposition + * is incorrect; swapping x2<->x4 moves the reflex vertex onto the + * shared diagonal endpoints and makes x1-x3 interior. + */ + void ensure_valid_diagonal(); }; /** diff --git a/coretrace/simulation_runner/optix_runner/OptixCSP/src/core/CspElement.cpp b/coretrace/simulation_runner/optix_runner/OptixCSP/src/core/CspElement.cpp index 7ade3585..8c8be291 100644 --- a/coretrace/simulation_runner/optix_runner/OptixCSP/src/core/CspElement.cpp +++ b/coretrace/simulation_runner/optix_runner/OptixCSP/src/core/CspElement.cpp @@ -224,10 +224,11 @@ GeometryDataST CspElement::toDeviceGeometryData() const if (surface_type == SurfaceType::PARABOLIC) { - Vec3d edge_x = v1 * (float)(-width); - Vec3d edge_y = v2 * (float)height; + Vec3d edge_x = v1 * (float)(width); + Vec3d edge_y = v2 * (float)(height); - Vec3d local_anchor(x_coord + width, y_coord, 0.0); + // Lower left corner + Vec3d local_anchor(x_coord, y_coord, 0.0); // float3 anchor = OptixCSP::toFloat3(m_origin - v1 * 0.5 - v2 * 0.5); Vec3d global_anchor = rotation_matrix * local_anchor + m_origin; @@ -264,16 +265,40 @@ GeometryDataST CspElement::toDeviceGeometryData() const v2 = tri.get_v1(); v3 = tri.get_v2(); - // given the origin and rotation, compute global coordinates of the triangle vertices - Matrix33d rotation_matrix = get_rotation_matrix(); // L2G rotation matrix - Vec3d v1_global = rotation_matrix * v1 + m_origin; - Vec3d v2_global = rotation_matrix * v2 + m_origin; - Vec3d v3_global = rotation_matrix * v3 + m_origin; - - GeometryDataST::Triangle_Flat heliostat(OptixCSP::toFloat3(v1_global), - OptixCSP::toFloat3(v2_global), - OptixCSP::toFloat3(v3_global)); - geometry_data.setTriangle_Flat(heliostat); + if (surface_type == SurfaceType::FLAT) + { + // given the origin and rotation, compute global coordinates of the triangle vertices + Matrix33d rotation_matrix = get_rotation_matrix(); // L2G rotation matrix + Vec3d v1_global = rotation_matrix * v1 + m_origin; + Vec3d v2_global = rotation_matrix * v2 + m_origin; + Vec3d v3_global = rotation_matrix * v3 + m_origin; + + GeometryDataST::Triangle_Flat heliostat(OptixCSP::toFloat3(v1_global), + OptixCSP::toFloat3(v2_global), + OptixCSP::toFloat3(v3_global)); + geometry_data.setTriangle_Flat(heliostat); + } + + if (surface_type == SurfaceType::PARABOLIC) + { + Matrix33d rotation_matrix = get_rotation_matrix(); // L2G rotation matrix + float3 vx = OptixCSP::toFloat3(rotation_matrix.get_x_basis()); + float3 vy = OptixCSP::toFloat3(rotation_matrix.get_y_basis()); + + float cx = (float)(m_surface->get_curvature_1()); + float cy = (float)(m_surface->get_curvature_2()); + float3 o = OptixCSP::toFloat3(m_origin); + + // Triangle vertices are in local element XY frame (z=0). + // The 2D local x,y coords for the aperture test equal the + // Vec3d x and y components directly (since rotation is orthonormal). + float2 lv0 = make_float2((float)v1[0], (float)v1[1]); + float2 lv1 = make_float2((float)v2[0], (float)v2[1]); + float2 lv2 = make_float2((float)v3[0], (float)v3[1]); + + GeometryDataST::Triangle_Parabolic trip(o, vx, vy, cx, cy, lv0, lv1, lv2); + geometry_data.setTriangle_Parabolic(trip); + } } if (aperture_type == ApertureType::QUADRILATERAL) @@ -287,18 +312,38 @@ GeometryDataST CspElement::toDeviceGeometryData() const p3 = quad.get_p2(); p4 = quad.get_p3(); - // given the origin and rotation, compute global coordinates of the triangle vertices - Matrix33d rotation_matrix = get_rotation_matrix(); // L2G rotation matrix - Vec3d p1_global = rotation_matrix * p1 + m_origin; - Vec3d p2_global = rotation_matrix * p2 + m_origin; - Vec3d p3_global = rotation_matrix * p3 + m_origin; - Vec3d p4_global = rotation_matrix * p4 + m_origin; - - GeometryDataST::Quadrilateral_Flat heliostat(OptixCSP::toFloat3(p1_global), - OptixCSP::toFloat3(p2_global), - OptixCSP::toFloat3(p3_global), - OptixCSP::toFloat3(p4_global)); - geometry_data.setQuadrilateral_Flat(heliostat); + if (surface_type == SurfaceType::FLAT) + { + // given the origin and rotation, compute global coordinates of the triangle vertices + Matrix33d rotation_matrix = get_rotation_matrix(); // L2G rotation matrix + Vec3d p1_global = rotation_matrix * p1 + m_origin; + Vec3d p2_global = rotation_matrix * p2 + m_origin; + Vec3d p3_global = rotation_matrix * p3 + m_origin; + Vec3d p4_global = rotation_matrix * p4 + m_origin; + + GeometryDataST::Quadrilateral_Flat heliostat(OptixCSP::toFloat3(p1_global), + OptixCSP::toFloat3(p2_global), + OptixCSP::toFloat3(p3_global), + OptixCSP::toFloat3(p4_global)); + geometry_data.setQuadrilateral_Flat(heliostat); + } + + if (surface_type == SurfaceType::PARABOLIC) + { + Matrix33d rotation_matrix = get_rotation_matrix(); // L2G rotation matrix + float3 vx = OptixCSP::toFloat3(rotation_matrix.get_x_basis()); + float3 vy = OptixCSP::toFloat3(rotation_matrix.get_y_basis()); + + float cx = (float)(m_surface->get_curvature_1()); + float cy = (float)(m_surface->get_curvature_2()); + float3 o = OptixCSP::toFloat3(m_origin); + + GeometryDataST::Quadrilateral_Parabolic heliostat( + o, vx, vy, cx, cy, + make_float2(p1[0], p1[1]), make_float2(p2[0], p2[1]), + make_float2(p3[0], p3[1]), make_float2(p4[0], p4[1])); + geometry_data.setQuadrilateral_Parabolic(heliostat); + } } if (aperture_type == ApertureType::CIRCLE) @@ -307,39 +352,79 @@ GeometryDataST CspElement::toDeviceGeometryData() const float r = circ.get_radius(); float3 o = OptixCSP::toFloat3(m_origin); float3 n = normalize(OptixCSP::toFloat3(m_aim_point - m_origin)); + if (surface_type == SurfaceType::FLAT) { GeometryDataST::Circle_Flat heliostat(o, n, r); geometry_data.setCircle_Flat(heliostat); } + + if (surface_type == SurfaceType::PARABOLIC) + { + Matrix33d rotation_matrix = get_rotation_matrix(); // L2G rotation matrix + float3 v1 = OptixCSP::toFloat3(rotation_matrix.get_x_basis()); + float3 v2 = OptixCSP::toFloat3(rotation_matrix.get_y_basis()); + float cx = (float)(m_surface->get_curvature_1()); + float cy = (float)(m_surface->get_curvature_2()); + GeometryDataST::Circle_Parabolic heliostat(o, v1, v2, cx, cy, r); + geometry_data.setCircle_Parabolic(heliostat); + } } if (aperture_type == ApertureType::HEXAGON) { + Matrix33d rotation_matrix = get_rotation_matrix(); // L2G rotation matrix + float3 vx = OptixCSP::toFloat3(rotation_matrix.get_x_basis()); + float3 vy = OptixCSP::toFloat3(rotation_matrix.get_y_basis()); + ApertureHexagon hex = static_cast(*m_aperture); float s = hex.get_side_length(); float3 o = OptixCSP::toFloat3(m_origin); float3 n = normalize(OptixCSP::toFloat3(m_aim_point - m_origin)); + if (surface_type == SurfaceType::FLAT) { - GeometryDataST::Hexagon_Flat hex(o, n, s); + GeometryDataST::Hexagon_Flat hex(o, n, vx, vy, s); geometry_data.setHexagon_Flat(hex); } + + if (surface_type == SurfaceType::PARABOLIC) + { + float cx = (float)(m_surface->get_curvature_1()); + float cy = (float)(m_surface->get_curvature_2()); + GeometryDataST::Hexagon_Parabolic hex(o, vx, vy, cx, cy, s); + geometry_data.setHexagon_Parabolic(hex); + } } if (aperture_type == ApertureType::ANNULUS) { + Matrix33d rotation_matrix = get_rotation_matrix(); // L2G rotation matrix + float3 vx = OptixCSP::toFloat3(rotation_matrix.get_x_basis()); + float3 vy = OptixCSP::toFloat3(rotation_matrix.get_y_basis()); + ApertureAnnulus anf = static_cast(*m_aperture); float radius_in = anf.get_radius_inner(); float radius_out = anf.get_radius_outer(); float arc = anf.get_arc(); float3 o = OptixCSP::toFloat3(m_origin); float3 n = normalize(OptixCSP::toFloat3(m_aim_point - m_origin)); + if (surface_type == SurfaceType::FLAT) { - GeometryDataST::Annulus_Flat anf(o, n, radius_in, radius_out, arc); + GeometryDataST::Annulus_Flat anf(o, n, vx, vy, + radius_in, radius_out, arc); geometry_data.setAnnulus_Flat(anf); } + + if (surface_type == SurfaceType::PARABOLIC) + { + float cx = (float)(m_surface->get_curvature_1()); + float cy = (float)(m_surface->get_curvature_2()); + GeometryDataST::Annulus_Parabolic anp(o, vx, vy, cx, cy, + radius_in, radius_out, arc); + geometry_data.setAnnulus_Parabolic(anp); + } } geometry_data.id = this->m_id; diff --git a/coretrace/simulation_runner/optix_runner/OptixCSP/src/core/geometry_manager.cpp b/coretrace/simulation_runner/optix_runner/OptixCSP/src/core/geometry_manager.cpp index 118f69b9..3d49bf15 100644 --- a/coretrace/simulation_runner/optix_runner/OptixCSP/src/core/geometry_manager.cpp +++ b/coretrace/simulation_runner/optix_runner/OptixCSP/src/core/geometry_manager.cpp @@ -46,6 +46,10 @@ void GeometryManager::collect_geometry_info(const std::vector(OpticalEntityType::ANNULUS_FLAT); } + else if (element->get_surface_type() == SurfaceType::PARABOLIC) + { + sbt_offset = static_cast(OpticalEntityType::ANNULUS_PARABOLIC); + } else { std::stringstream ss; @@ -63,6 +67,10 @@ void GeometryManager::collect_geometry_info(const std::vector(OpticalEntityType::CIRCLE_FLAT); } + else if (element->get_surface_type() == SurfaceType::PARABOLIC) + { + sbt_offset = static_cast(OpticalEntityType::CIRCLE_PARABOLIC); + } else { std::stringstream ss; @@ -80,6 +88,10 @@ void GeometryManager::collect_geometry_info(const std::vector(OpticalEntityType::HEXAGON_FLAT); } + else if (element->get_surface_type() == SurfaceType::PARABOLIC) + { + sbt_offset = static_cast(OpticalEntityType::HEXAGON_PARABOLIC); + } else { std::stringstream ss; @@ -122,6 +134,10 @@ void GeometryManager::collect_geometry_info(const std::vector(OpticalEntityType::TRIANGLE_FLAT); } + else if (element->get_surface_type() == SurfaceType::PARABOLIC) + { + sbt_offset = static_cast(OpticalEntityType::TRIANGLE_PARABOLIC); + } else { std::stringstream ss; @@ -139,6 +155,10 @@ void GeometryManager::collect_geometry_info(const std::vector(OpticalEntityType::QUADRILATERAL_FLAT); } + else if (element->get_surface_type() == SurfaceType::PARABOLIC) + { + sbt_offset = static_cast(OpticalEntityType::QUADRILATERAL_PARABOLIC); + } else { std::stringstream ss; diff --git a/coretrace/simulation_runner/optix_runner/OptixCSP/src/core/pipeline_manager.cpp b/coretrace/simulation_runner/optix_runner/OptixCSP/src/core/pipeline_manager.cpp index ade2345d..7f6b40ca 100644 --- a/coretrace/simulation_runner/optix_runner/OptixCSP/src/core/pipeline_manager.cpp +++ b/coretrace/simulation_runner/optix_runner/OptixCSP/src/core/pipeline_manager.cpp @@ -45,7 +45,12 @@ const std::map IntersectionKernelMap = { {OpticalEntityType::QUADRILATERAL_FLAT, "__intersection__quadrilateral_flat"}, {OpticalEntityType::CIRCLE_FLAT, "__intersection__circle_flat"}, {OpticalEntityType::HEXAGON_FLAT, "__intersection__hexagon_flat"}, - {OpticalEntityType::ANNULUS_FLAT, "__intersection__annulus_flat"}}; + {OpticalEntityType::ANNULUS_FLAT, "__intersection__annulus_flat"}, + {OpticalEntityType::CIRCLE_PARABOLIC, "__intersection__circle_parabolic"}, + {OpticalEntityType::HEXAGON_PARABOLIC, "__intersection__hexagon_parabolic"}, + {OpticalEntityType::TRIANGLE_PARABOLIC, "__intersection__triangle_parabolic"}, + {OpticalEntityType::ANNULUS_PARABOLIC, "__intersection__annulus_parabolic"}, + {OpticalEntityType::QUADRILATERAL_PARABOLIC, "__intersection__quadrilateral_parabolic"}}; pipelineManager::pipelineManager(SoltraceState &state) : m_state(state) {} diff --git a/coretrace/simulation_runner/optix_runner/OptixCSP/src/core/soltrace_type.h b/coretrace/simulation_runner/optix_runner/OptixCSP/src/core/soltrace_type.h index d98569a8..ff8703fa 100644 --- a/coretrace/simulation_runner/optix_runner/OptixCSP/src/core/soltrace_type.h +++ b/coretrace/simulation_runner/optix_runner/OptixCSP/src/core/soltrace_type.h @@ -31,12 +31,12 @@ namespace OptixCSP SurfaceType surfaceType; ApertureType apertureType; - SurfaceApertureMap() : surfaceType(SurfaceType::FLAT), - apertureType(ApertureType::RECTANGLE) - { - } + SurfaceApertureMap() : surfaceType(SurfaceType::FLAT), + apertureType(ApertureType::RECTANGLE) + { + } - SurfaceApertureMap(SurfaceType surf, ApertureType ap) : surfaceType(surf), + SurfaceApertureMap(SurfaceType surf, ApertureType ap) : surfaceType(surf), apertureType(ap) { } @@ -48,19 +48,19 @@ namespace OptixCSP return (surfaceType == map.surfaceType) && (apertureType == map.apertureType); } - bool operator<(const SurfaceApertureMap &b) const - { - return surfaceType < b.surfaceType || - (surfaceType == b.surfaceType && apertureType < b.apertureType); - } + bool operator<(const SurfaceApertureMap &b) const + { + return surfaceType < b.surfaceType || + (surfaceType == b.surfaceType && apertureType < b.apertureType); + } - friend std::ostream& operator<<(std::ostream &os, const SurfaceApertureMap &sam); + friend std::ostream &operator<<(std::ostream &os, const SurfaceApertureMap &sam); }; - inline std::ostream& operator<<(std::ostream &os, const SurfaceApertureMap &sam) + inline std::ostream &operator<<(std::ostream &os, const SurfaceApertureMap &sam) { - os << "SurfApMap -- Surface: " << static_cast(sam.surfaceType) - << " Aperture: " << static_cast(sam.apertureType); - return os; + os << "SurfApMap -- Surface: " << static_cast(sam.surfaceType) + << " Aperture: " << static_cast(sam.apertureType); + return os; } } diff --git a/coretrace/simulation_runner/optix_runner/OptixCSP/src/shaders/GeometryDataST.h b/coretrace/simulation_runner/optix_runner/OptixCSP/src/shaders/GeometryDataST.h index 217e64dd..97608211 100644 --- a/coretrace/simulation_runner/optix_runner/OptixCSP/src/shaders/GeometryDataST.h +++ b/coretrace/simulation_runner/optix_runner/OptixCSP/src/shaders/GeometryDataST.h @@ -25,7 +25,12 @@ namespace OptixCSP QUADRILATERAL_FLAT = 6, CIRCLE_FLAT = 7, HEXAGON_FLAT = 8, - ANNULUS_FLAT = 9 + ANNULUS_FLAT = 9, + CIRCLE_PARABOLIC = 10, + HEXAGON_PARABOLIC = 11, + TRIANGLE_PARABOLIC = 12, + ANNULUS_PARABOLIC = 13, + QUADRILATERAL_PARABOLIC = 14 }; struct Parallelogram @@ -40,10 +45,10 @@ namespace OptixCSP this->v2 *= 1.0f / dot(v2, v2); plane = make_float4(normal, d); } - float4 plane; - float3 v1; - float3 v2; - float3 anchor; + float4 plane; // plane equation: (normal, dot(anchor, normal)) + float3 v1; // edge vector 1, stored as v1/dot(v1,v1) + float3 v2; // edge vector 2, stored as v2/dot(v2,v2) + float3 anchor; // corner point of the parallelogram }; // same as parallelogram, however defined with different attributes @@ -58,12 +63,12 @@ namespace OptixCSP plane = make_float4(normal, d); } - float4 plane; - float3 center; - float3 x; - float3 y; - float width; - float height; + float4 plane; // normal unit vector, dot(center, normal) + float3 center; // center in global coordinates + float3 x; // local x axis unit vector + float3 y; // local y axis unit vector + float width; // full width along x + float height; // full height along y }; struct Cylinder_Y @@ -75,11 +80,11 @@ namespace OptixCSP assert(dot(base_x, base_z) < 1e-3f); } - float3 center; - float radius; - float half_height; - float3 base_x; // x axis of the cylinder - float3 base_z; // z axis of the cylinder + float3 center; // center of the cylinder in global coordinates + float radius; // radius + float half_height; // half the height along the local y axis + float3 base_x; // x axis of the cylinder + float3 base_z; // z axis of the cylinder }; struct Rectangle_Parabolic @@ -96,13 +101,13 @@ namespace OptixCSP plane = make_float4(normal, d); } - float4 plane; - float3 v1; - float3 v2; - float3 anchor; + float4 plane; // plane equation of the base rectangle: (normal, dot(anchor, normal)) + float3 v1; // edge vector 1, stored as v1/dot(v1,v1) + float3 v2; // edge vector 2, stored as v2/dot(v2,v2) + float3 anchor; // corner point of the base rectangle // float3 focus; - float curv_x; - float curv_y; + float curv_x; // curvature along local x axis + float curv_y; // curvature along local y axis }; struct Triangle_Flat @@ -116,8 +121,8 @@ namespace OptixCSP } float3 v0; // base vertex float3 e1, e2; // edges - float3 normal; - float d; // plane distance + float3 normal; // normal unit vector + float d; // plane distance }; struct Quadrilateral_Flat @@ -142,40 +147,169 @@ namespace OptixCSP Circle_Flat(const float3 &origin, const float3 &normal, const float &radius) : r(radius), center(origin) { - plane = make_float4(normalize(normal), dot(center, normal)); + const float3 n = normalize(normal); + plane = make_float4(n, dot(center, n)); } - float4 plane; - float3 center; - float r; + float4 plane; // normal unit vector, dot(center, normal) + float3 center; // local origin in global coordinates + float r; // radius }; struct Hexagon_Flat { Hexagon_Flat() = default; - Hexagon_Flat(const float3 &origin, const float3 &normal, const float &side_length) - : s(side_length), center(origin) + Hexagon_Flat(const float3 &origin, const float3 &normal, + const float3 &x_ax, const float3 &y_ax, + const float &side_length) + : center(origin), x_axis(x_ax), y_axis(y_ax), s(side_length) { - plane = make_float4(normalize(normal), dot(center, normal)); + const float3 n = normalize(normal); + plane = make_float4(n, dot(center, n)); } - float4 plane; - float3 center; - float s; + float4 plane; // normal unit vector, dot(center, normal) + float3 center; // local origin in global coordinates + float3 x_axis; // unit vector + float3 y_axis; // unit vector + float s; // side length }; struct Annulus_Flat { Annulus_Flat() = default; Annulus_Flat(const float3 &origin, const float3 &normal, - const float &r_inner, const float &r_outer, const float &arc) - : center(origin), ri(r_inner), ro(r_outer), arc(arc) + const float3 &x_ax, const float3 &y_ax, + const float &r_inner, const float &r_outer, + const float &arc) + : center(origin), x_axis(x_ax), y_axis(y_ax), + ri(r_inner), ro(r_outer), arc(arc) + { + const float3 n = normalize(normal); + plane = make_float4(n, dot(center, n)); + } + float4 plane; // normal unit vector, dot(center, normal) + float3 center; // local origin in global coordinates + float3 x_axis; // local x axis unit vector (arc is centered about this axis) + float3 y_axis; // local y axis unit vector + float ri; // inner radius + float ro; // outer radius + float arc; // total arc angle in radians, centered on x_axis + }; + + struct Circle_Parabolic + { + Circle_Parabolic() = default; + Circle_Parabolic(const float3 &origin, const float3 &x_ax, const float3 &y_ax, + const float &curv_x, const float &curv_y, const float &r) + : center(origin), x_axis(x_ax), y_axis(y_ax), + cx(curv_x), cy(curv_y), radius(r) + { + } + float3 center; + float3 x_axis; + float3 y_axis; + float cx; + float cy; + float radius; + }; + + struct Hexagon_Parabolic + { + Hexagon_Parabolic() = default; + Hexagon_Parabolic(const float3 &origin, const float3 &x_ax, const float3 &y_ax, + const float &curv_x, const float &curv_y, const float &side_len) + : center(origin), x_axis(x_ax), y_axis(y_ax), + cx(curv_x), cy(curv_y), s(side_len) { - plane = make_float4(normalize(normal), dot(center, normal)); } - float4 plane; float3 center; + float3 x_axis; + float3 y_axis; + float cx; + float cy; + float s; + }; + + struct Triangle_Parabolic + { + Triangle_Parabolic() = default; + // Vertices v0, v1, v2 are in the local XY aperture frame. + // The constructor precomputes the barycentric inverse transform so + // the aperture test reduces to two dot products with no per-ray division. + Triangle_Parabolic(const float3 &origin, const float3 &x_ax, const float3 &y_ax, + const float &curv_x, const float &curv_y, + const float2 &v0, const float2 &v1, const float2 &v2) + : center(origin), x_axis(x_ax), y_axis(y_ax), + cx(curv_x), cy(curv_y) + { + const float2 e1 = make_float2(v1.x - v0.x, v1.y - v0.y); + const float2 e2 = make_float2(v2.x - v0.x, v2.y - v0.y); + const float inv_det = 1.0f / (e1.x * e2.y - e1.y * e2.x); + // u = dot(utest, float3(px, py, 1.0f)) + utest = make_float3(e2.y, -e2.x, v0.y * e2.x - v0.x * e2.y) * inv_det; + // v = dot(vtest, float3(px, py, 1.0f)) + vtest = make_float3(-e1.y, e1.x, v0.x * e1.y - v0.y * e1.x) * inv_det; + } + float3 center; // element origin in global coordinates + float3 x_axis; // local x axis unit vector + float3 y_axis; // local y axis unit vector + float cx; // curvature along local x axis + float cy; // curvature along local y axis + float3 utest; // precomputed row: u = dot(utest, float3(px, py, 1.0f)) + float3 vtest; // precomputed row: v = dot(vtest, float3(px, py, 1.0f)) + }; + + struct Annulus_Parabolic + { + Annulus_Parabolic() = default; + Annulus_Parabolic(const float3 &origin, const float3 &x_ax, const float3 &y_ax, + const float &curv_x, const float &curv_y, + const float &r_inner, const float &r_outer, const float &arc) + : center(origin), x_axis(x_ax), y_axis(y_ax), + cx(curv_x), cy(curv_y), + ri(r_inner), ro(r_outer), arc(arc) + { + } + float3 center; + float3 x_axis; + float3 y_axis; + float cx; + float cy; float ri; float ro; - float arc; // Arc angle in radians with x-axis in the middle + float arc; // in radians + }; + + struct Quadrilateral_Parabolic + { + Quadrilateral_Parabolic() = default; + Quadrilateral_Parabolic(const float3 &origin, const float3 &x_ax, const float3 &y_ax, + const float &curv_x, const float &curv_y, + const float2 &v0, const float2 &v1, + const float2 &v2, const float2 &v3) + : center(origin), x_axis(x_ax), y_axis(y_ax), cx(curv_x), cy(curv_y) + { + const float2 e1 = make_float2(v1.x - v0.x, v1.y - v0.y); + const float2 e2 = make_float2(v2.x - v0.x, v2.y - v0.y); + const float inv_det = 1.0f / (e1.x * e2.y - e1.y * e2.x); + // u = dot(utest, float3(px, py, 1.0f)) + u1test = make_float3(e2.y, -e2.x, v0.y * e2.x - v0.x * e2.y) * inv_det; + // v = dot(vtest, float3(px, py, 1.0f)) + v1test = make_float3(-e1.y, e1.x, v0.x * e1.y - v0.y * e1.x) * inv_det; + + const float2 e3 = make_float2(v3.x - v0.x, v3.y - v0.y); + const float inv_det2 = 1.0f / (e2.x * e3.y - e2.y * e3.x); + u2test = make_float3(e3.y, -e3.x, v0.y * e3.x - v0.x * e3.y) * inv_det2; + v2test = make_float3(-e2.y, e2.x, v0.x * e2.y - v0.y * e2.x) * inv_det2; + } + float3 center; // element origin in global coordinates + float3 x_axis; // local x axis unit vector + float3 y_axis; // local y axis unit vector + float cx; // curvature along local x axis + float cy; // curvature along local y axis + float3 u1test; // precomputed row: u = dot(utest, float3(px, py, 1.0f)) + float3 v1test; // precomputed row: v = dot(vtest, float3(px, py, 1.0f)) + float3 u2test; + float3 v2test; }; GeometryDataST() = default; @@ -297,6 +431,71 @@ namespace OptixCSP return annulus_flat; } + void setCircle_Parabolic(const Circle_Parabolic &circp) + { + assert(type == UNKNOWN_TYPE); + type = CIRCLE_PARABOLIC; + circle_parabolic = circp; + } + + __host__ __device__ const Circle_Parabolic &getCircle_Parabolic() const + { + assert(type == CIRCLE_PARABOLIC); + return circle_parabolic; + } + + void setHexagon_Parabolic(const Hexagon_Parabolic &hexp) + { + assert(type == UNKNOWN_TYPE); + type = HEXAGON_PARABOLIC; + hexagon_parabolic = hexp; + } + + __host__ __device__ const Hexagon_Parabolic &getHexagon_Parabolic() const + { + assert(type == HEXAGON_PARABOLIC); + return hexagon_parabolic; + } + + void setTriangle_Parabolic(const Triangle_Parabolic &tp) + { + assert(type == UNKNOWN_TYPE); + type = TRIANGLE_PARABOLIC; + triangle_parabolic = tp; + } + + __host__ __device__ const Triangle_Parabolic &getTriangle_Parabolic() const + { + assert(type == TRIANGLE_PARABOLIC); + return triangle_parabolic; + } + + void setAnnulus_Parabolic(const Annulus_Parabolic &ap) + { + assert(type == UNKNOWN_TYPE); + type = ANNULUS_PARABOLIC; + annulus_parabolic = ap; + } + + __host__ __device__ const Annulus_Parabolic &getAnnulus_Parabolic() const + { + assert(type == ANNULUS_PARABOLIC); + return annulus_parabolic; + } + + void setQuadrilateral_Parabolic(const Quadrilateral_Parabolic &qp) + { + assert(type == UNKNOWN_TYPE); + type = QUADRILATERAL_PARABOLIC; + quadrilateral_parabolic = qp; + } + + __host__ __device__ const Quadrilateral_Parabolic &getQuadrilateral_Parabolic() const + { + assert(type == QUADRILATERAL_PARABOLIC); + return quadrilateral_parabolic; + } + Type type = UNKNOWN_TYPE; int32_t id = OptixCSP::kElementIdUnassigned; @@ -313,6 +512,11 @@ namespace OptixCSP Circle_Flat circle_flat; Hexagon_Flat hexagon_flat; Annulus_Flat annulus_flat; + Circle_Parabolic circle_parabolic; + Hexagon_Parabolic hexagon_parabolic; + Triangle_Parabolic triangle_parabolic; + Annulus_Parabolic annulus_parabolic; + Quadrilateral_Parabolic quadrilateral_parabolic; }; }; } diff --git a/coretrace/simulation_runner/optix_runner/OptixCSP/src/shaders/Soltrace.h b/coretrace/simulation_runner/optix_runner/OptixCSP/src/shaders/Soltrace.h index 1c4c87ab..5426b0e6 100644 --- a/coretrace/simulation_runner/optix_runner/OptixCSP/src/shaders/Soltrace.h +++ b/coretrace/simulation_runner/optix_runner/OptixCSP/src/shaders/Soltrace.h @@ -36,14 +36,19 @@ namespace OptixCSP{ }; enum OpticalEntityType : unsigned int { - RECTANGLE_FLAT = 0, - RECTANGLE_PARABOLIC = 1, - CYLINDRICAL = 2, - TRIANGLE_FLAT = 3, - QUADRILATERAL_FLAT = 4, - CIRCLE_FLAT = 5, - HEXAGON_FLAT = 6, - ANNULUS_FLAT = 7, + RECTANGLE_FLAT = 0, + RECTANGLE_PARABOLIC = 1, + CYLINDRICAL = 2, + TRIANGLE_FLAT = 3, + QUADRILATERAL_FLAT = 4, + CIRCLE_FLAT = 5, + HEXAGON_FLAT = 6, + ANNULUS_FLAT = 7, + CIRCLE_PARABOLIC = 8, + HEXAGON_PARABOLIC = 9, + TRIANGLE_PARABOLIC = 10, + ANNULUS_PARABOLIC = 11, + QUADRILATERAL_PARABOLIC = 12, NUM_OPTICAL_ENTITY_TYPES }; diff --git a/coretrace/simulation_runner/optix_runner/OptixCSP/src/shaders/intersection.cu b/coretrace/simulation_runner/optix_runner/OptixCSP/src/shaders/intersection.cu index b9688d39..092619e7 100644 --- a/coretrace/simulation_runner/optix_runner/OptixCSP/src/shaders/intersection.cu +++ b/coretrace/simulation_runner/optix_runner/OptixCSP/src/shaders/intersection.cu @@ -9,12 +9,235 @@ extern "C" __constant__ OptixCSP::LaunchParams params; } +/**************** Surface Helper Functions ****************/ + +// ----------------------------------------------------------------------- +// Shared helpers for a planar (flat) surface +// +// All planar surfaces have the equation +// = 0 +// where n is the normal vector, p0 = (x0, y0, z0) is a point in the plane, +// generally the origin of the aperture. +// ----------------------------------------------------------------------- + +// Return the ray parameter t at which the ray (ro + t*rd) intersects the plane. +// The plane is encoded as float4(nx, ny, nz, d) where nx/ny/nz is the unit normal +// and d = dot(n, p0) for any point p0 on the plane. extern "C" __device__ __inline__ float ray_distance_to_plane(float3 ro, float3 rd, float4 plane) { const float3 n = make_float3(plane); return (plane.w - dot(n, ro)) / dot(rd, n); } +// ----------------------------------------------------------------------- +// Shared helpers for parabolic surface intersections. +// +// All parabolic surfaces share the same quadric equation: +// z = (cx/2)*x^2 + (cy/2)*y^2 +// in a local frame (center, x_ax, y_ax, n=cross(x_ax,y_ax)). +// The three helpers below factor out the ray transform, quadratic solve, +// and normal computation. Each kernel only supplies the aperture test. +// ----------------------------------------------------------------------- + +// Transform a world-space ray into the local parabolic frame. +// Outputs the frame normal n = normalize(cross(x_ax, y_ax)) and +// the local ray origin (ox,oy,oz) and direction (dx,dy,dz). +extern "C" __device__ __inline__ void parabolic_ray_to_local( + const float3 &ray_orig, const float3 &ray_dir, + const float3 ¢er, + const float3 &x_ax, const float3 &y_ax, + float3 &n, + float &ox, float &oy, float &oz, + float &dx, float &dy, float &dz) +{ + n = normalize(cross(x_ax, y_ax)); + const float3 d = ray_orig - center; + ox = dot(d, x_ax); + oy = dot(d, y_ax); + oz = dot(d, n); + dx = dot(ray_dir, x_ax); + dy = dot(ray_dir, y_ax); + dz = dot(ray_dir, n); +} + +// Solve A*t^2 + B*t + C = 0 for the paraboloid-ray intersection and return +// up to two hits within [ray_tmin, ray_tmax], ordered by ascending t. +// Returns the number of valid hits (0, 1, or 2). +// t_out[i], lx_out[i], ly_out[i] give the ray parameter and local (x,y) of each hit. +extern "C" __device__ __inline__ int parabolic_solve( + float ox, float oy, float oz, + float dx, float dy, float dz, + float cx, float cy, + float ray_tmin, float ray_tmax, + float t_out[2], float lx_out[2], float ly_out[2]) +{ + const float A = 0.5f * cx * dx * dx + 0.5f * cy * dy * dy; + const float B = cx * ox * dx + cy * oy * dy - dz; + const float C = 0.5f * cx * ox * ox + 0.5f * cy * oy * oy - oz; + + const float eps = 1e-12f; + int count = 0; + + if (fabsf(A) < eps) + { + if (fabsf(B) > eps) + { + const float t = -C / B; + if (t >= ray_tmin && t <= ray_tmax) + { + t_out[0] = t; + lx_out[0] = ox + t * dx; + ly_out[0] = oy + t * dy; + count = 1; + } + } + } + else + { + const float discr = B * B - 4.0f * A * C; + if (discr >= 0.0f) + { + const float sq = sqrtf(discr); + // A > 0 (physical curvature), so ta <= tb is guaranteed. + const float ta = -0.5f * (B + sq) / A; + const float tb = -0.5f * (B - sq) / A; + if (ta >= ray_tmin && ta <= ray_tmax) + { + t_out[count] = ta; + lx_out[count] = ox + ta * dx; + ly_out[count] = oy + ta * dy; + ++count; + } + if (tb >= ray_tmin && tb <= ray_tmax) + { + t_out[count] = tb; + lx_out[count] = ox + tb * dx; + ly_out[count] = oy + tb * dy; + ++count; + } + } + } + return count; +} + +// Compute the world-space unit normal at a parabolic surface hit. +// x_hit, y_hit : local coordinates of the hit point +// cx, cy : curvature parameters +// x_ax, y_ax : local frame unit vectors +// n : normalize(cross(x_ax, y_ax)) +extern "C" __device__ __inline__ float3 parabolic_world_normal( + float x_hit, float y_hit, + float cx, float cy, + const float3 &x_ax, const float3 &y_ax, const float3 &n) +{ + const float3 N_local = make_float3(-cx * x_hit, -cy * y_hit, 1.0f); + return N_local.x * x_ax + N_local.y * y_ax + N_local.z * n; +} + +// ----------------------------------------------------------------------- + +/**************** Aperture Helper Functions ****************/ + +// ----------------------------------------------------------------------- +// Shared helpers for hexagon apertures +// +// A regular hexagon with circumradius s (vertex-to-center distance) is +// decomposed into three vertical strips in the local (x, y) frame: +// Left cap : x in [-s, -s/2) — bounded by the two left edges +// Center : x in [-s/2, s/2] — full-height rectangular band +// Right cap : x in ( s/2, s] — bounded by the two right edges +// The flat (top/bottom) edges are horizontal at y = ±(sqrt(3)/2)*s and +// the diagonal edges satisfy |y| = sqrt(3)*(|x| - s) on each cap. +// ----------------------------------------------------------------------- + +// Return true if the point (px, py) lies within a flat-top regular hexagon +// centered at the origin with circumradius s (vertex-to-center distance). +// px and py must be expressed in the hexagon's local x/y frame. +extern "C" __device__ __inline__ bool hexagon_contains(float px, float py, float s) +{ + bool is_in = false; + const float xl = 0.5f * s; + const float yl = 0.5f * sqrtf(3.0f) * s; + if (-xl <= px && px <= xl && -yl <= py && py <= yl) + { + // Center + is_in = true; + } + else if (-s <= px && px < -xl) + { + // Left side + float y1 = sqrtf(3.0f) * (px + s); + float y2 = -y1; + if (y2 <= py && py <= y1) + { + is_in = true; + } + } + else if (xl < px && px <= s) + { + // Right side + float y1 = sqrtf(3.0f) * (px - s); + float y2 = -y1; + if (y1 <= py && py <= y2) + { + is_in = true; + } + } + return is_in; +} + +// ----------------------------------------------------------------------- +// Shared helpers for annulus apertures +// +// An annular sector aperture is defined in the local (x, y) frame by: +// Inner radius ri : points closer than ri to the origin are excluded +// Outer radius ro : points farther than ro from the origin are excluded +// Arc angle arc : full sweep angle (radians) centered on the +x axis; +// a point at angle theta is included when |theta| <= arc/2 +// A full annulus (no arc clipping) has arc = 2*pi. +// ----------------------------------------------------------------------- + +// Return true if the point (px, py) lies within an annular sector centered at the +// origin. The sector is defined by inner radius ri, outer radius ro, and a full +// arc angle arc (in radians) symmetric about the local x-axis (i.e., theta = 0). +// px and py must be expressed in the aperture's local x/y frame. +extern "C" __device__ __inline__ bool annulus_contains(float px, float py, float ri, float ro, float arc) +{ + const float rsq = px * px + py * py; + if (rsq < ri * ri || rsq > ro * ro) + return false; + const float theta = atan2f(py, px); + return fabsf(theta) <= 0.5f * arc; +} + +// ----------------------------------------------------------------------- +// Shared helpers for triangle apertures +// +// The triangle aperture is encoded as two precomputed row vectors utest and +// vtest such that the barycentric coordinates of a point (px, py) are: +// u = dot(utest, float3(px, py, 1.0f)) +// v = dot(vtest, float3(px, py, 1.0f)) +// The point is inside when u >= 0, v >= 0, and u + v <= 1. +// These rows are computed once at setup time from the three triangle vertices. +// ----------------------------------------------------------------------- + +// Return true if the 2-D point (px, py) lies inside the triangle whose +// barycentric inverse transform is encoded in utest and vtest. +extern "C" __device__ __inline__ bool triangle_contains(float px, float py, + float3 utest, float3 vtest) +{ + const float3 p = make_float3(px, py, 1.0f); + const float u = dot(utest, p); + if (u < 0.0f || u > 1.0f) + return false; + const float v = dot(vtest, p); + return v >= 0.0f && (u + v) <= 1.0f; +} + +// ----------------------------------------------------------------------- + +/**************** Optix Intersection Functions ****************/ + extern "C" __global__ void __intersection__parallelogram() { int i = optixGetPrimitiveIndex(); @@ -303,171 +526,6 @@ extern "C" __global__ void __intersection__cylinder_y_capped() __float_as_uint(world_normal.z)); } -// For a parabolic surface rectangle aperture where -// the base (normal projection) is defined by the center and its two unit edge vectors -// In a local coordinate system (with origin at the anchor) the flat rectangle covers: -// x in [0, L1] and y in [0, L2], -// where L1 and L2 are the lengths of the original edge vectors. -// The parabolic surface is given by: -// z = (curv_x/2)*x^2 + (curv_y/2)*y^2 -// and the ray (in local coordinates) is: -// (ox,oy,oz) + t*(dx,dy,dz) -// We solve for t such that: -// oz + t*dz = (curv_x/2)*(ox+t*dx)^2 + (curv_y/2)*(oy+t*dy)^2 -// which expands into a quadratic: A*t^2 + B*t + C = 0. -// After finding the valid t, we compute the local hit (x,y) and then check that -// 0 <= x <= L1 and 0 <= y <= L2. -// Finally, we compute the surface normal from the paraboloid derivative -// f_x = curv_x * x and f_y = curv_y * y, -// so that the (unnormalized) local normal is (-f_x, -f_y, 1). -// -// The local hit point is then transformed back to world space for reporting. -extern "C" __global__ void __intersection__rectangle_parabolic() -{ - const OptixCSP::GeometryDataST::Rectangle_Parabolic &rect = params.geometry_data_array[optixGetPrimitiveIndex()].getRectangleParabolic(); - // Get ray information. - const float3 ray_orig = optixGetWorldRayOrigin(); - const float3 ray_dir = optixGetWorldRayDirection(); - const float ray_tmin = optixGetRayTmin(); - const float ray_tmax = optixGetRayTmax(); - - // - // Build the local coordinate system. - // - // We assume that the rectangle was defined with an anchor at its corner and - // two edge vectors. The stored rect.v1 and rect.v2 are the reciprocals: - // stored_v1 = original_v1 / dot(original_v1, original_v1) - // Thus, the original edge lengths are: - // Note rect.v1 had the size of 1/original_v1_length - float L1 = 1.0f / length(rect.v1); - float L2 = 1.0f / length(rect.v2); - // And the unit edge directions are: - float3 e1 = rect.v1 * L1; // recovers the original direction of edge 1, unit vector - float3 e2 = rect.v2 * L2; // recovers the original direction of edge 2, unit ve - // The flat (undeformed) rectangle's normal is: - float3 n = normalize(cross(e2, e1)); - - // - // Transform ray into local coordinates. - // The local coordinates (x,y,z) are defined such that: - // - The origin is at rect.anchor. - // - The x-axis is e1. - // - The y-axis is e2. - // - The z-axis is n. - // - // Compute the rectangle center (shifting from the lower-right corner) - float3 rect_center = rect.anchor + (L1 / 2.0f) * e1 + (L2 / 2.0f) * e2; - - float3 d = ray_orig - rect_center; - float ox = dot(d, e1); - float oy = dot(d, e2); - float oz = dot(d, n); - - float dx = dot(ray_dir, e1); - float dy = dot(ray_dir, e2); - float dz = dot(ray_dir, n); - - // Retrieve curvature parameters. - const float curv_x = rect.curv_x; - const float curv_y = rect.curv_y; - - float A = (curv_x * 0.5f) * (dx * dx) + (curv_y * 0.5f) * (dy * dy); - float B = curv_x * (ox * dx) + curv_y * (oy * dy) - dz; - float C = (curv_x * 0.5f) * (ox * ox) + (curv_y * 0.5f) * (oy * oy) - oz; - - float t = 0.0f; - float t1 = 0.0f; - float t2 = 0.0f; - float x_hit = 0.0f; - float y_hit = 0.0f; - const float eps = 1e-12f; - bool valid = false; - - if (fabsf(A) < eps) { - if (fabsf(B) > eps) { - t1 = -C / B; - const float p1x = ox + t1 * dx; - const float p1y = oy + t1 * dy; - const float a1 = p1x / (L1 / 2.0f); - const float a2 = p1y / (L2 / 2.0f); - - if (t1 > 0.0f && t1 >= ray_tmin && t1 <= ray_tmax - && !(a1 < -1.0f || a1 > 1.0f || a2 < -1.0f || a2 > 1.0f)) { - t = t1; - x_hit = p1x; - y_hit = p1y; - valid = true; - } - } - } - else - { - float discr = B * B - 4.0f * A * C; - if (discr >= 0.0f) - { - float sqrt_discr = sqrtf(discr); - t1 = -0.5f * (B + sqrt_discr) / A; - t2 = -0.5f * (B - sqrt_discr) / A; - - const float p1x = ox + t1 * dx; - const float p1y = oy + t1 * dy; - const float p2x = ox + t2 * dx; - const float p2y = oy + t2 * dy; - const float a1_1 = p1x / (L1 / 2.0f); - const float a2_1 = p1y / (L2 / 2.0f); - const float a1_2 = p2x / (L1 / 2.0f); - const float a2_2 = p2y / (L2 / 2.0f); - - if (t1 > 0.0f && t1 >= ray_tmin && t1 <= ray_tmax - && !(a1_1 < -1.0f || a1_1 > 1.0f || a2_1 < -1.0f || a2_1 > 1.0f)) { - t = t1; - x_hit = p1x; - y_hit = p1y; - valid = true; - } - else if (t2 > 0.0f && t2 >= ray_tmin && t2 <= ray_tmax - && !(a1_2 < -1.0f || a1_2 > 1.0f || a2_2 < -1.0f || a2_2 > 1.0f)) { - t = t2; - x_hit = p2x; - y_hit = p2y; - valid = true; - } - } - } - - if (!valid) { - return; - } - - // - // Compute the surface normal at the hit on the paraboloid. - // The height function is: - // f(x,y) = (curv_x/2)*x^2 + (curv_y/2)*y^2 - // so its partial derivatives are: - // f_x = curv_x * x and f_y = curv_y * y. - // Then the (unnormalized) local normal is: - // N_local = (-f_x, -f_y, 1) = ( -curv_x*x_hit, -curv_y*y_hit, 1 ). - // - float3 N_local = normalize(make_float3(-curv_x * x_hit, - -curv_y * y_hit, - 1.0f)); - // Transform the normal back to world coordinates. - float3 world_normal = normalize(N_local.x * e1 + - N_local.y * e2 + - N_local.z * n); - - // Compute the hit point in world space. - float3 world_hit = ray_orig + t * ray_dir; - - // Report the intersection. - // Here, the two reported extra attributes are the parametric coordinates (a1, a2), - // encoded as unsigned integers. - optixReportIntersection(t, 0, - __float_as_uint(world_normal.x), - __float_as_uint(world_normal.y), - __float_as_uint(world_normal.z)); -} - // intersection algorithm for a flat triangle based on "Fast, Minimum Storage Ray/Triangle Intersection" by M�ller and Trumbore (1997) // code from here: https://en.wikipedia.org/wiki/M%C3%B6ller%E2%80%93Trumbore_intersection_algorithm extern "C" __device__ __inline__ float _triangle_intersect( @@ -505,51 +563,6 @@ extern "C" __device__ __inline__ float _triangle_intersect( return t; } -// // intersection algorithm for a flat triangle based on "Fast, Minimum Storage Ray/Triangle Intersection" by M�ller and Trumbore (1997) -// // code from here: https://en.wikipedia.org/wiki/M%C3%B6ller%E2%80%93Trumbore_intersection_algorithm -// extern "C" __global__ void __intersection__triangle_flat() -// { -// const OptixCSP::GeometryDataST::Triangle_Flat& tri = params.geometry_data_array[optixGetPrimitiveIndex()].getTriangle_Flat(); - -// const float3 ro = optixGetObjectRayOrigin(); -// const float3 rd = optixGetObjectRayDirection(); - -// //printf("Ray origin: (%f,%f,%f), direction: (%f,%f,%f)\n", ro.x, ro.y, ro.z, rd.x, rd.y, rd.z); - -// const float3 edge1 = tri.e1; -// const float3 edge2 = tri.e2; - -// const float3 pvec = cross(rd, edge2); -// const float det = dot(edge1, pvec); - -// // Backface culling + parallel rejection -// // (det must be strictly positive and not tiny) -// const float eps = 1e-8f; -// if (det <= eps) return; - -// const float inv_det = 1.0f / det; - -// const float3 tvec = ro - tri.v0; -// const float u = dot(tvec, pvec) * inv_det; -// if (u < 0.0f || u > 1.0f) return; - -// const float3 qvec = cross(tvec, edge1); -// const float v = dot(rd, qvec) * inv_det; -// if (v < 0.0f || (u + v) > 1.0f) -// return; - -// const float t = dot(edge2, qvec) * inv_det; -// if (t < optixGetRayTmin() || t > optixGetRayTmax()) return; - -// float3 world_normal = tri.normal; - -// optixReportIntersection(t, 0, -// __float_as_uint(world_normal.x), -// __float_as_uint(world_normal.y), -// __float_as_uint(world_normal.z)); - -// } - // intersection algorithm for a flat triangle based on "Fast, Minimum Storage Ray/Triangle Intersection" by M�ller and Trumbore (1997) // code from here: https://en.wikipedia.org/wiki/M%C3%B6ller%E2%80%93Trumbore_intersection_algorithm extern "C" __global__ void __intersection__triangle_flat() @@ -579,18 +592,17 @@ extern "C" __global__ void __intersection__quadrilateral_flat() const float3 ro = optixGetObjectRayOrigin(); const float3 rd = optixGetObjectRayDirection(); + // Decompose as (p0,p1,p2) ∪ (p0,p2,p3) — same p0-p2 diagonal as the parabolic kernel. const float3 p0 = quad.p0; - const float3 p2 = quad.p2; + const float3 e02 = quad.p2 - p0; float3 e1 = quad.p1 - p0; - float3 e2 = quad.p3 - p0; - float t = _triangle_intersect(p0, e1, e2, ro, rd); + float t = _triangle_intersect(p0, e1, e02, ro, rd); if (t < optixGetRayTmin() || t > optixGetRayTmax()) { - e1 = quad.p1 - p2; - e2 = quad.p3 - p2; - t = _triangle_intersect(p2, e1, e2, ro, rd); + const float3 e3 = quad.p3 - p0; + t = _triangle_intersect(p0, e02, e3, ro, rd); } if (t < optixGetRayTmin() || t > optixGetRayTmax()) @@ -624,6 +636,9 @@ extern "C" __global__ void __intersection__circle_flat() // Verify intersection distance and Report ray intersection point if (t > ray_tmin && t < ray_tmax) { + // Since everything is in global coordinates (e.g., the ray intersection coordinates), + // and the circle is rotationally symmetric, we don't need to worry about + // any rotation of the circle in local coordinates float3 p = ray_orig + ray_dir * t; float d = length(p - circ.center); if (d <= circ.r) @@ -652,40 +667,13 @@ extern "C" __global__ void __intersection__hexagon_flat() // Verify intersection distance and Report ray intersection point if (t > ray_tmin && t < ray_tmax) { - // TODO: Need to adjust for possible rotation... - bool is_in = false; float3 p = ray_orig + ray_dir * t - hex.center; - // float d = length(p - circ.center); - float s = hex.s; - float xl = 0.5f * s; - float yl = 0.5f * sqrtf(3.0f) * s; - if (-xl <= p.x && p.x <= xl && -yl <= p.y && p.y <= yl) - { - // Center - is_in = true; - } - else if (-s <= p.x && p.x < xl) - { - // Left side - float y1 = sqrtf(3.0f) * (p.x + s); - float y2 = -y1; - if (y2 <= p.y && p.y <= y1) - { - is_in = true; - } - } - else if (xl < p.x && p.x <= s) - { - // Right side - float y1 = sqrtf(3.0f) * (p.x - s); - float y2 = -y1; - if (y1 <= p.y && p.y <= y2) - { - is_in = true; - } - } + // Project onto the local x and y axes which are unit vectors + const float px = dot(p, hex.x_axis); + const float py = dot(p, hex.y_axis); + const float s = hex.s; - if (is_in) + if (hexagon_contains(px, py, s)) { optixReportIntersection(t, 0, @@ -711,20 +699,270 @@ extern "C" __global__ void __intersection__annulus_flat() // Verify intersection distance and Report ray intersection point if (t > ray_tmin && t < ray_tmax) { - // TODO: Need to adjust for possible rotation... float3 p = ray_orig + ray_dir * t - anf.center; - float d = length(p); - if (anf.ri <= d && d <= anf.ro) + const float px = dot(p, anf.x_axis); + const float py = dot(p, anf.y_axis); + if (annulus_contains(px, py, anf.ri, anf.ro, anf.arc)) { - float theta = atan2f(p.y, p.x); - if (fabsf(theta) <= 0.5f * anf.arc) - { - optixReportIntersection(t, - 0, - __float_as_uint(n.x), - __float_as_uint(n.y), - __float_as_uint(n.z)); - } + optixReportIntersection(t, + 0, + __float_as_uint(n.x), + __float_as_uint(n.y), + __float_as_uint(n.z)); + } + } +} + +// Parabolic surface, rectangle aperture. +// z = (curv_x/2)*x^2 + (curv_y/2)*y^2, aperture: |x| <= L1/2, |y| <= L2/2. +extern "C" __global__ void __intersection__rectangle_parabolic() +{ + const OptixCSP::GeometryDataST::Rectangle_Parabolic &rect = + params.geometry_data_array[optixGetPrimitiveIndex()].getRectangleParabolic(); + + const float3 ray_orig = optixGetWorldRayOrigin(); + const float3 ray_dir = optixGetWorldRayDirection(); + const float ray_tmin = optixGetRayTmin(); + const float ray_tmax = optixGetRayTmax(); + + // Recover unit edge vectors and half-lengths from the stored reciprocal vectors. + // rect.v1 = original_v1 / dot(original_v1, original_v1), so |rect.v1| = 1/L1. + const float L1 = 1.0f / length(rect.v1); + const float L2 = 1.0f / length(rect.v2); + const float3 e1 = rect.v1 * L1; + const float3 e2 = rect.v2 * L2; + const float3 center = rect.anchor + (0.5f * L1) * e1 + (0.5f * L2) * e2; + + float3 n; + float ox, oy, oz, dx, dy, dz; + parabolic_ray_to_local(ray_orig, ray_dir, center, e1, e2, + n, ox, oy, oz, dx, dy, dz); + + float ts[2], lxs[2], lys[2]; + const int nc = parabolic_solve(ox, oy, oz, dx, dy, dz, + rect.curv_x, rect.curv_y, + ray_tmin, ray_tmax, ts, lxs, lys); + + const float half_L1 = 0.5f * L1; + const float half_L2 = 0.5f * L2; + for (int i = 0; i < nc; ++i) + { + if (lxs[i] >= -half_L1 && lxs[i] <= half_L1 && + lys[i] >= -half_L2 && lys[i] <= half_L2) + { + const float3 wn = parabolic_world_normal(lxs[i], lys[i], + rect.curv_x, rect.curv_y, + e1, e2, n); + optixReportIntersection(ts[i], 0, + __float_as_uint(wn.x), + __float_as_uint(wn.y), + __float_as_uint(wn.z)); + return; + } + } +} + +// Parabolic surface, circle aperture. +// z = (cx/2)*x^2 + (cy/2)*y^2, aperture: x^2 + y^2 <= radius^2. +extern "C" __global__ void __intersection__circle_parabolic() +{ + const OptixCSP::GeometryDataST::Circle_Parabolic &circp = + params.geometry_data_array[optixGetPrimitiveIndex()].getCircle_Parabolic(); + + const float3 ray_orig = optixGetWorldRayOrigin(); + const float3 ray_dir = optixGetWorldRayDirection(); + const float ray_tmin = optixGetRayTmin(); + const float ray_tmax = optixGetRayTmax(); + + float3 n; + float ox, oy, oz, dx, dy, dz; + parabolic_ray_to_local(ray_orig, ray_dir, + circp.center, circp.x_axis, circp.y_axis, + n, ox, oy, oz, dx, dy, dz); + + float ts[2], lxs[2], lys[2]; + const int nc = parabolic_solve(ox, oy, oz, dx, dy, dz, + circp.cx, circp.cy, + ray_tmin, ray_tmax, + ts, lxs, lys); + + const float r2 = circp.radius * circp.radius; + for (int i = 0; i < nc; ++i) + { + if (lxs[i] * lxs[i] + lys[i] * lys[i] <= r2) + { + const float3 wn = parabolic_world_normal(lxs[i], lys[i], + circp.cx, circp.cy, + circp.x_axis, circp.y_axis, n); + optixReportIntersection(ts[i], 0, + __float_as_uint(wn.x), + __float_as_uint(wn.y), + __float_as_uint(wn.z)); + return; + } + } +} + +extern "C" __global__ void __intersection__hexagon_parabolic() +{ + const OptixCSP::GeometryDataST::Hexagon_Parabolic &hexp = + params.geometry_data_array[optixGetPrimitiveIndex()].getHexagon_Parabolic(); + + const float3 ray_orig = optixGetWorldRayOrigin(); + const float3 ray_dir = optixGetWorldRayDirection(); + const float ray_tmin = optixGetRayTmin(); + const float ray_tmax = optixGetRayTmax(); + + float3 n; + float ox, oy, oz, dx, dy, dz; + parabolic_ray_to_local(ray_orig, ray_dir, + hexp.center, hexp.x_axis, hexp.y_axis, + n, ox, oy, oz, dx, dy, dz); + + float ts[2], lxs[2], lys[2]; + const int nc = parabolic_solve(ox, oy, oz, dx, dy, dz, + hexp.cx, hexp.cy, + ray_tmin, ray_tmax, + ts, lxs, lys); + + for (int i = 0; i < nc; ++i) + { + float3 p = ray_orig + ray_dir * ts[i] - hexp.center; + // Project onto the local x and y axes which are unit vectors + const float px = dot(p, hexp.x_axis); + const float py = dot(p, hexp.y_axis); + const float s = hexp.s; + + if (hexagon_contains(px, py, s)) + { + const float3 wn = parabolic_world_normal(lxs[i], lys[i], + hexp.cx, hexp.cy, + hexp.x_axis, hexp.y_axis, n); + + optixReportIntersection(ts[i], + 0, + __float_as_uint(wn.x), + __float_as_uint(wn.y), + __float_as_uint(wn.z)); + + return; + } + } +} + +extern "C" __global__ void __intersection__triangle_parabolic() +{ + const OptixCSP::GeometryDataST::Triangle_Parabolic &trip = + params.geometry_data_array[optixGetPrimitiveIndex()].getTriangle_Parabolic(); + + const float3 ray_orig = optixGetWorldRayOrigin(); + const float3 ray_dir = optixGetWorldRayDirection(); + const float ray_tmin = optixGetRayTmin(); + const float ray_tmax = optixGetRayTmax(); + + float3 n; + float ox, oy, oz, dx, dy, dz; + parabolic_ray_to_local(ray_orig, ray_dir, + trip.center, trip.x_axis, trip.y_axis, + n, ox, oy, oz, dx, dy, dz); + + float ts[2], lxs[2], lys[2]; + const int nc = parabolic_solve(ox, oy, oz, dx, dy, dz, + trip.cx, trip.cy, + ray_tmin, ray_tmax, + ts, lxs, lys); + + for (int i = 0; i < nc; ++i) + { + if (triangle_contains(lxs[i], lys[i], trip.utest, trip.vtest)) + { + const float3 wn = parabolic_world_normal(lxs[i], lys[i], + trip.cx, trip.cy, + trip.x_axis, trip.y_axis, + n); + optixReportIntersection(ts[i], 0, + __float_as_uint(wn.x), + __float_as_uint(wn.y), + __float_as_uint(wn.z)); + return; + } + } +} + +extern "C" __global__ void __intersection__annulus_parabolic() +{ + const OptixCSP::GeometryDataST::Annulus_Parabolic &anap = + params.geometry_data_array[optixGetPrimitiveIndex()].getAnnulus_Parabolic(); + + const float3 ray_orig = optixGetWorldRayOrigin(); + const float3 ray_dir = optixGetWorldRayDirection(); + const float ray_tmin = optixGetRayTmin(); + const float ray_tmax = optixGetRayTmax(); + + float3 n; + float ox, oy, oz, dx, dy, dz; + parabolic_ray_to_local(ray_orig, ray_dir, + anap.center, anap.x_axis, anap.y_axis, + n, ox, oy, oz, dx, dy, dz); + + float ts[2], lxs[2], lys[2]; + const int nc = parabolic_solve(ox, oy, oz, dx, dy, dz, + anap.cx, anap.cy, + ray_tmin, ray_tmax, + ts, lxs, lys); + + for (int i = 0; i < nc; ++i) + { + if (annulus_contains(lxs[i], lys[i], anap.ri, anap.ro, anap.arc)) + { + const float3 wn = parabolic_world_normal(lxs[i], lys[i], + anap.cx, anap.cy, + anap.x_axis, anap.y_axis, + n); + optixReportIntersection(ts[i], 0, + __float_as_uint(wn.x), + __float_as_uint(wn.y), + __float_as_uint(wn.z)); + return; + } + } +} + +extern "C" __global__ void __intersection__quadrilateral_parabolic() +{ + const OptixCSP::GeometryDataST::Quadrilateral_Parabolic &quap = + params.geometry_data_array[optixGetPrimitiveIndex()].getQuadrilateral_Parabolic(); + + const float3 ray_orig = optixGetWorldRayOrigin(); + const float3 ray_dir = optixGetWorldRayDirection(); + const float ray_tmin = optixGetRayTmin(); + const float ray_tmax = optixGetRayTmax(); + + float3 n; + float ox, oy, oz, dx, dy, dz; + parabolic_ray_to_local(ray_orig, ray_dir, + quap.center, quap.x_axis, quap.y_axis, + n, ox, oy, oz, dx, dy, dz); + + float ts[2], lxs[2], lys[2]; + const int nc = parabolic_solve(ox, oy, oz, dx, dy, dz, + quap.cx, quap.cy, + ray_tmin, ray_tmax, + ts, lxs, lys); + + for (int i = 0; i < nc; ++i) + { + if (triangle_contains(lxs[i], lys[i], quap.u1test, quap.v1test) || + triangle_contains(lxs[i], lys[i], quap.u2test, quap.v2test)) + { + const float3 wn = parabolic_world_normal(lxs[i], lys[i], + quap.cx, quap.cy, + quap.x_axis, quap.y_axis, n); + optixReportIntersection(ts[i], 0, + __float_as_uint(wn.x), + __float_as_uint(wn.y), + __float_as_uint(wn.z)); + return; } } } diff --git a/coretrace/simulation_runner/optix_runner/OptixCSP/src/shaders/materials.cu b/coretrace/simulation_runner/optix_runner/OptixCSP/src/shaders/materials.cu index 004dbf73..2b73f562 100644 --- a/coretrace/simulation_runner/optix_runner/OptixCSP/src/shaders/materials.cu +++ b/coretrace/simulation_runner/optix_runner/OptixCSP/src/shaders/materials.cu @@ -58,7 +58,7 @@ extern "C" __device__ __inline__ float3 orthonormal_vector(float3 v) // w = make_float3(0.0f, 1.0f, 0.0f); // u = cross(v, w); u = make_float3(-v.z, 0.0f, v.x); - } + } return normalize(u); } @@ -104,10 +104,13 @@ extern "C" __global__ void __closesthit__element() const bool optical_errors = params.optical_errors; // Fetch the normal vector from the hit attributes passed by OptiX - float3 object_normal = make_float3(__uint_as_float(optixGetAttribute_0()), __uint_as_float(optixGetAttribute_1()), + float3 object_normal = make_float3(__uint_as_float(optixGetAttribute_0()), + __uint_as_float(optixGetAttribute_1()), __uint_as_float(optixGetAttribute_2())); - // Transform the object-space normal to world space using OptiX built-in function - float3 world_normal = normalize(optixTransformNormalFromObjectToWorldSpace(object_normal)); + // // Transform the object-space normal to world space using OptiX built-in function + // float3 world_normal = normalize(optixTransformNormalFromObjectToWorldSpace(object_normal)); + // All reports from the intersection functions are in world coordinates + float3 world_normal = normalize(object_normal); // Compute the facing normal, which handles the direction of the normal based on the incoming ray direction const float3 ray_dir = optixGetWorldRayDirection(); @@ -140,8 +143,8 @@ extern "C" __global__ void __closesthit__element() const float transmissivity = material.transmissivity; const bool use_transmissivity = material.use_refraction; const float reflectivity = material.reflectivity; - const float normal_sigma = 1e-3f * material.slope_error; // Convert mrad to rad - const float spec_sigma = 1e-3f * material.specularity_error; // Convert mrad to rad + const float normal_sigma = 1e-3f * material.slope_error; // Convert mrad to rad + const float spec_sigma = 1e-3f * material.specularity_error; // Convert mrad to rad const uint8_t error_type = material.optical_dist; // Surface normal (macro-surface) errors @@ -288,7 +291,7 @@ extern "C" __global__ void __closesthit__element() else if (!absorbed) { // Ray hit an element but max depth was reached; count it (absorption at this depth does not count). - atomicAdd(reinterpret_cast(params.d_depth_exceeded_count), 1ULL); + atomicAdd(reinterpret_cast(params.d_depth_exceeded_count), 1ULL); } setPayload(prd); diff --git a/coretrace/simulation_runner/optix_runner/optix_runner.cpp b/coretrace/simulation_runner/optix_runner/optix_runner.cpp index 61136def..a72fe809 100644 --- a/coretrace/simulation_runner/optix_runner/optix_runner.cpp +++ b/coretrace/simulation_runner/optix_runner/optix_runner.cpp @@ -345,10 +345,19 @@ RunnerStatus OptixRunner::setup_elements(const SimulationData *data) auto el_aperture = std::dynamic_pointer_cast(el->get_aperture()); assert(el_aperture != nullptr); + OptixCSP::Vec3d p0(el_aperture->x1, el_aperture->y1, 0.0); OptixCSP::Vec3d p1(el_aperture->x2, el_aperture->y2, 0.0); OptixCSP::Vec3d p2(el_aperture->x3, el_aperture->y3, 0.0); + // Ensure CCW winding (right-hand rule) required by ApertureTriangle + // Aperture always lies in x-y-plane with positive z-axis corresponding to + // the front of the geometric element so CCW winding always gives the front + // of the element. + const double signed_area = (p1[0] - p0[0]) * (p2[1] - p0[1]) - (p1[1] - p0[1]) * (p2[0] - p0[0]); + if (signed_area < 0.0) + std::swap(p1, p2); + auto aperture = std::make_shared(p0, p1, p2); optix_el->set_aperture(aperture); @@ -364,6 +373,16 @@ RunnerStatus OptixRunner::setup_elements(const SimulationData *data) OptixCSP::Vec3d p2(el_aperture->x3, el_aperture->y3, 0.0); OptixCSP::Vec3d p3(el_aperture->x4, el_aperture->y4, 0.0); + // Ensure CCW winding using the shoelace signed-area formula. + // For a simple (non-self-intersecting) quad the sign of 2*A tells + // the winding without decomposing into triangles. + const double signed_area2 = (p0[0] * p1[1] - p1[0] * p0[1]) + + (p1[0] * p2[1] - p2[0] * p1[1]) + + (p2[0] * p3[1] - p3[0] * p2[1]) + + (p3[0] * p0[1] - p0[0] * p3[1]); + if (signed_area2 < 0.0) + std::swap(p1, p3); // reverse winding, keeping p0 and p2 fixed + auto aperture = std::make_shared(p0, p1, p2, p3); optix_el->set_aperture(aperture); diff --git a/google-tests/unit-tests/simulation_data/aperture_test.cpp b/google-tests/unit-tests/simulation_data/aperture_test.cpp index ddf0f642..50944ef1 100644 --- a/google-tests/unit-tests/simulation_data/aperture_test.cpp +++ b/google-tests/unit-tests/simulation_data/aperture_test.cpp @@ -600,6 +600,93 @@ TEST(Aperture, IrregularQuadrilateral) std::cout << "Expected fraction of hits: " << frac << std::endl; } +TEST(Aperture, IrregularQuadrilateral_VertexOrder) +{ + // Verify that the same parallelogram described with four different vertex + // orderings — CCW, CW (reverse), a cyclic rotation, and a self-intersecting + // "bowtie" ordering — all produce identical area and is_in results after + // ensure_valid_diagonal() normalises the representation. + // + // Canonical CCW vertices: (0,0) -> (3,0) -> (4,2) -> (1,2) + const double TOL = 1e-12; + const double AREA = 6.0; // parallelogram: base 3 * height 2 + // One interior and one exterior point, chosen away from edges. + const double IX = 2.5, IY = 1.0; // clearly inside + const double OX = 4.5, OY = 1.0; // clearly outside + + // Helper: extract the four normalised vertices as an array of (x,y) pairs. + using V4 = std::array, 4>; + auto get_verts = [](const aperture_ptr &ap) -> V4 { + const auto *q = dynamic_cast(ap.get()); + return {{ {q->x1,q->y1}, {q->x2,q->y2}, {q->x3,q->y3}, {q->x4,q->y4} }}; + }; + // Returns true when both arrays contain the same four (x,y) pairs (any order). + auto same_vertex_set = [&](const V4 &a, const V4 &b) -> bool { + for (const auto &va : a) { + bool found = false; + for (const auto &vb : b) + if (std::abs(va.first - vb.first) < TOL && + std::abs(va.second - vb.second) < TOL) + { found = true; break; } + if (!found) return false; + } + return true; + }; + // Returns true when the x1-x3 diagonal of q is interior to the quad, + // i.e. x2 and x4 are on strictly opposite sides of the x1-x3 line. + auto diagonal_interior = [&](const aperture_ptr &ap) -> bool { + const auto *q = dynamic_cast(ap.get()); + const double ex = q->x3 - q->x1, ey = q->y3 - q->y1; + const double d2 = ex * (q->y2 - q->y1) - ey * (q->x2 - q->x1); + const double d4 = ex * (q->y4 - q->y1) - ey * (q->x4 - q->x1); + return d2 * d4 < 0.0; + }; + + // 1) CCW (canonical) + auto q_ccw = make_aperture( + 0.0, 0.0, 3.0, 0.0, 4.0, 2.0, 1.0, 2.0); + EXPECT_NEAR(q_ccw->aperture_area(), AREA, TOL) << "CCW area"; + EXPECT_TRUE( q_ccw->is_in(IX, IY)) << "CCW inside"; + EXPECT_FALSE(q_ccw->is_in(OX, OY)) << "CCW outside"; + EXPECT_TRUE(diagonal_interior(q_ccw)) << "CCW diagonal interior"; + const V4 canonical = get_verts(q_ccw); + + // 2) CW (full reversal): (0,0) -> (1,2) -> (4,2) -> (3,0) + // The x1-x3 diagonal (0,0)-(4,2) is already interior for this ordering, + // so ensure_valid_diagonal() is a no-op and the vertices stay CW. + auto q_cw = make_aperture( + 0.0, 0.0, 1.0, 2.0, 4.0, 2.0, 3.0, 0.0); + EXPECT_NEAR(q_cw->aperture_area(), AREA, TOL) << "CW area"; + EXPECT_TRUE( q_cw->is_in(IX, IY)) << "CW inside"; + EXPECT_FALSE(q_cw->is_in(OX, OY)) << "CW outside"; + EXPECT_TRUE(diagonal_interior(q_cw)) << "CW diagonal interior"; + EXPECT_TRUE(same_vertex_set(canonical, get_verts(q_cw))) << "CW vertex set"; + + // 3) Cyclic rotation of CCW: start from (3,0) + // (3,0) -> (4,2) -> (1,2) -> (0,0) + // The x1-x3 diagonal (3,0)-(1,2) is interior, so ensure_valid_diagonal() + // is again a no-op. + auto q_rot = make_aperture( + 3.0, 0.0, 4.0, 2.0, 1.0, 2.0, 0.0, 0.0); + EXPECT_NEAR(q_rot->aperture_area(), AREA, TOL) << "rotated area"; + EXPECT_TRUE( q_rot->is_in(IX, IY)) << "rotated inside"; + EXPECT_FALSE(q_rot->is_in(OX, OY)) << "rotated outside"; + EXPECT_TRUE(diagonal_interior(q_rot)) << "rotated diagonal interior"; + EXPECT_TRUE(same_vertex_set(canonical, get_verts(q_rot))) << "rotated vertex set"; + + // 4) Self-intersecting "bowtie" order: (0,0) -> (4,2) -> (3,0) -> (1,2) + // The x1-x3 diagonal test fails (same-sign cross products), so + // ensure_valid_diagonal() angle-sorts and cyclic-shifts to recover a + // valid simple polygon with an interior x1-x3 diagonal. + auto q_bowtie = make_aperture( + 0.0, 0.0, 4.0, 2.0, 3.0, 0.0, 1.0, 2.0); + EXPECT_NEAR(q_bowtie->aperture_area(), AREA, TOL) << "bowtie area"; + EXPECT_TRUE( q_bowtie->is_in(IX, IY)) << "bowtie inside"; + EXPECT_FALSE(q_bowtie->is_in(OX, OY)) << "bowtie outside"; + EXPECT_TRUE(diagonal_interior(q_bowtie)) << "bowtie diagonal interior"; + EXPECT_TRUE(same_vertex_set(canonical, get_verts(q_bowtie))) << "bowtie vertex set"; +} + TEST(Aperture, MakeApertureFromType) { const double TOL = 1e-12; diff --git a/google-tests/unit-tests/simulation_runner/optix_runner/geometry_intersection_test.cpp b/google-tests/unit-tests/simulation_runner/optix_runner/geometry_intersection_test.cpp index b3f065d8..80795be1 100644 --- a/google-tests/unit-tests/simulation_runner/optix_runner/geometry_intersection_test.cpp +++ b/google-tests/unit-tests/simulation_runner/optix_runner/geometry_intersection_test.cpp @@ -1,10 +1,12 @@ #include +#include #include #include #include #include +using SolTrace::Data::D2R; using SolTrace::Runner::RunnerStatus; const double Z_ELEM = 50.0; @@ -12,9 +14,49 @@ const double Z_BACKSTOP = Z_ELEM - 0.5 * Z_ELEM; const double TOL = 1e-6; const uint_fast64_t NRAYS = 10000; +// Adapted from CspElement::set_bounding_box_local. +// Computes the global axis-aligned bounding box for a local bounding box +// transformed by a z-axis rotation (in degrees) and origin offset. +static void compute_global_bounding_box( + const glm::dvec3 &lower_local, + const glm::dvec3 &upper_local, + double zrot_deg, // z-axis rotation in degrees (local-to-global) + const glm::dvec3 &origin, + glm::dvec3 &lower_global, + glm::dvec3 &upper_global) +{ + const double cr = cos(zrot_deg * D2R); + const double sr = sin(zrot_deg * D2R); + // glm is column-major: columns are (cr, sr, 0), (-sr, cr, 0), (0, 0, 1) + const glm::dmat3 rotation(cr, sr, 0.0, + -sr, cr, 0.0, + 0.0, 0.0, 1.0); + const glm::dvec3 c0 = lower_local; + const glm::dvec3 c7 = upper_local; + const glm::dvec3 corners[8] = { + rotation * glm::dvec3(c0[0], c0[1], c0[2]) + origin, + rotation * glm::dvec3(c0[0], c0[1], c7[2]) + origin, + rotation * glm::dvec3(c0[0], c7[1], c0[2]) + origin, + rotation * glm::dvec3(c7[0], c0[1], c0[2]) + origin, + rotation * glm::dvec3(c0[0], c7[1], c7[2]) + origin, + rotation * glm::dvec3(c7[0], c0[1], c7[2]) + origin, + rotation * glm::dvec3(c7[0], c7[1], c0[2]) + origin, + rotation * glm::dvec3(c7[0], c7[1], c7[2]) + origin, + }; + + lower_global = glm::dvec3(std::numeric_limits::max()); + upper_global = glm::dvec3(std::numeric_limits::lowest()); + for (const auto &c : corners) + { + lower_global = glm::min(lower_global, c); + upper_global = glm::max(upper_global, c); + } +} + element_id set_default_sd(SimulationData &sd, surface_ptr surf, - aperture_ptr ap) + aperture_ptr ap, + double rotation) { sd.clear(); @@ -36,6 +78,7 @@ element_id set_default_sd(SimulationData &sd, el->set_optical_property_set(opt_ref); el->set_name("el"); + el->set_zrot(rotation); // Add element to stage element_id id = sd.add_element(el); @@ -43,14 +86,24 @@ element_id set_default_sd(SimulationData &sd, // Back stop element that is bigger than the created element so that the // testing element casts a shadow on this big thing. element_ptr stop = make_element(); - double xlb, xub, ylb, yub; + double xlb, xub, ylb, yub, zlb, zub; ap->bounding_box(xlb, xub, ylb, yub); - const double sx = std::max(fabs(xlb), fabs(xub)) + 1.0; - const double sy = std::max(fabs(ylb), fabs(yub)) + 1.0; + surf->bounding_box(xlb, xub, ylb, yub, zlb, zub); + + // Transform local AABB into global AABB accounting for zrot + glm::dvec3 lower_global, upper_global; + compute_global_bounding_box(glm::dvec3(xlb, ylb, zlb), + glm::dvec3(xub, yub, zub), + rotation, + glm::dvec3(0.0, 0.0, Z_ELEM), + lower_global, upper_global); + + const double sx = std::max(fabs(lower_global[0]), fabs(upper_global[0])) + 2.0; + const double sy = std::max(fabs(lower_global[1]), fabs(upper_global[1])) + 2.0; stop->set_origin(0, 0, Z_BACKSTOP); stop->set_aim_vector(0, 0, 100); stop->set_surface(make_surface()); - stop->set_aperture(make_aperture(sx, sy)); + stop->set_aperture(make_aperture(2.0 * sx, 2.0 * sy)); stop->set_optical_property_set(opt_ref); sd.add_element(stop); @@ -69,11 +122,14 @@ TEST(OptixRunner, FlatRectangle) { const double XL = 10.0; const double YL = 5.0; + const double ROT_DEG = -10.0; auto surf = make_surface(); auto aper = make_aperture(XL, YL); + uint_fast64_t fpos = 0, fneg = 0, hits = 0, misses = 0; + SimulationData sd; - element_id test_elid = set_default_sd(sd, surf, aper); + element_id test_elid = set_default_sd(sd, surf, aper, ROT_DEG); SimulationResult result; OptixRunner runner; @@ -88,6 +144,8 @@ TEST(OptixRunner, FlatRectangle) ASSERT_EQ(result.get_number_of_records(), sd.get_simulation_parameters().number_of_rays); + const double cos_rot = cos(ROT_DEG * D2R); + const double sin_rot = sin(ROT_DEG * D2R); for (int i = 0; i < (int)result.get_number_of_records(); ++i) { auto rr = result[i]; @@ -98,32 +156,48 @@ TEST(OptixRunner, FlatRectangle) auto id = rr->get_element(1); EXPECT_NEAR(p0[0], p1[0], TOL) << "ray " << i; EXPECT_NEAR(p0[1], p1[1], TOL) << "ray " << i; + const double lx = p1[0] * cos_rot - p1[1] * sin_rot; + const double ly = p1[0] * sin_rot + p1[1] * cos_rot; if (id == test_elid) { // We hit the test element. Check that the height is as expected EXPECT_NEAR(p1[2], Z_ELEM, TOL * Z_ELEM) << "ray " << i; // And that we are in the aperture - EXPECT_TRUE(aper->is_in(p1[0], p1[1])); + EXPECT_TRUE(aper->is_in(lx, ly)); + ++hits; + if (!aper->is_in(lx, ly)) + ++fpos; } else { // We hit the back stop element. Check that the height is as expected EXPECT_NEAR(p1[2], Z_BACKSTOP, TOL * Z_ELEM); // And that we are not in the aperture. - EXPECT_FALSE(aper->is_in(p1[0], p1[1])); + EXPECT_FALSE(aper->is_in(lx, ly)); + ++misses; + if (aper->is_in(lx, ly)) + ++fneg; } } + EXPECT_GT(hits, 0u); + EXPECT_GT(misses, 0u); + std::cout << "hits: " << hits << ", misses: " << misses + << ", false positives: " << fpos + << ", false negatives: " << fneg << std::endl; } TEST(OptixRunner, FlatEquilateralTriangle) { const double d = 4.0; + const double ROT_DEG = 90.0; auto surf = make_surface(); auto aper = make_aperture(d); + uint_fast64_t fpos = 0, fneg = 0, hits = 0, misses = 0; + SimulationData sd; - element_id test_elid = set_default_sd(sd, surf, aper); + element_id test_elid = set_default_sd(sd, surf, aper, ROT_DEG); SimulationResult result; OptixRunner runner; @@ -138,6 +212,8 @@ TEST(OptixRunner, FlatEquilateralTriangle) ASSERT_EQ(result.get_number_of_records(), sd.get_simulation_parameters().number_of_rays); + const double cos_rot = cos(ROT_DEG * D2R); + const double sin_rot = sin(ROT_DEG * D2R); for (int i = 0; i < (int)result.get_number_of_records(); ++i) { auto rr = result[i]; @@ -148,29 +224,117 @@ TEST(OptixRunner, FlatEquilateralTriangle) auto id = rr->get_element(1); EXPECT_NEAR(p0[0], p1[0], TOL) << "ray " << i; EXPECT_NEAR(p0[1], p1[1], TOL) << "ray " << i; + const double lx = p1[0] * cos_rot - p1[1] * sin_rot; + const double ly = p1[0] * sin_rot + p1[1] * cos_rot; if (id == test_elid) { EXPECT_NEAR(p1[2], Z_ELEM, TOL * Z_ELEM) << "ray " << i; - EXPECT_TRUE(aper->is_in(p1[0], p1[1])); + EXPECT_TRUE(aper->is_in(lx, ly)); + ++hits; + if (!aper->is_in(lx, ly)) + ++fpos; } else { EXPECT_NEAR(p1[2], Z_BACKSTOP, TOL * Z_ELEM); - EXPECT_FALSE(aper->is_in(p1[0], p1[1])); + EXPECT_FALSE(aper->is_in(lx, ly)); + ++misses; + if (aper->is_in(lx, ly)) + ++fneg; } } + EXPECT_GT(hits, 0u); + EXPECT_GT(misses, 0u); + std::cout << "hits: " << hits << ", misses: " << misses + << ", false positives: " << fpos + << ", false negatives: " << fneg << std::endl; } TEST(OptixRunner, FlatTriangle) { + // Same triangle as FlatTriangle but vertices supplied in CCW order: + // (0,0)->(2,0)->(1,2). The optix_runner's cross-product sign check sees + // a positive area and does NOT swap, so the CCW path is taken directly. + // Results should be identical to FlatTriangle. + const double x1 = 0.0, x2 = 2.0, x3 = 1.0; + const double y1 = 0.0, y2 = 0.0, y3 = 2.0; + const double ROT_DEG = 110.0; + auto surf = make_surface(); + auto aper = make_aperture(x1, y1, x2, y2, x3, y3); + + uint_fast64_t fpos = 0, fneg = 0, hits = 0, misses = 0; + + SimulationData sd; + element_id test_elid = set_default_sd(sd, surf, aper, ROT_DEG); + SimulationResult result; + + OptixRunner runner; + RunnerStatus sts = runner.initialize(); + ASSERT_EQ(sts, RunnerStatus::SUCCESS); + sts = runner.setup_simulation(&sd); + ASSERT_EQ(sts, RunnerStatus::SUCCESS); + sts = runner.run_simulation(); + ASSERT_EQ(sts, RunnerStatus::SUCCESS); + sts = runner.report_simulation(&result, 0); + ASSERT_EQ(sts, RunnerStatus::SUCCESS); + + ASSERT_EQ(result.get_number_of_records(), + sd.get_simulation_parameters().number_of_rays); + const double cos_rot = cos(ROT_DEG * D2R); + const double sin_rot = sin(ROT_DEG * D2R); + for (int i = 0; i < (int)result.get_number_of_records(); ++i) + { + auto rr = result[i]; + ASSERT_GE(rr->get_number_of_interactions(), 2); + glm::dvec3 p0, p1; + rr->get_position(0, p0); + rr->get_position(1, p1); + auto id = rr->get_element(1); + EXPECT_NEAR(p0[0], p1[0], TOL) << "ray " << i; + EXPECT_NEAR(p0[1], p1[1], TOL) << "ray " << i; + const double lx = p1[0] * cos_rot - p1[1] * sin_rot; + const double ly = p1[0] * sin_rot + p1[1] * cos_rot; + + if (id == test_elid) + { + EXPECT_NEAR(p1[2], Z_ELEM, TOL * Z_ELEM) << "ray " << i; + EXPECT_TRUE(aper->is_in(lx, ly)); + ++hits; + if (!aper->is_in(lx, ly)) + ++fpos; + } + else + { + EXPECT_NEAR(p1[2], Z_BACKSTOP, TOL * Z_ELEM); + EXPECT_FALSE(aper->is_in(lx, ly)); + ++misses; + if (aper->is_in(lx, ly)) + ++fneg; + } + } + EXPECT_GT(hits, 0u); + EXPECT_GT(misses, 0u); + std::cout << "hits: " << hits << ", misses: " << misses + << ", false positives: " << fpos + << ", false negatives: " << fneg << std::endl; +} + +TEST(OptixRunner, FlatTriangle_CW) +{ + // Vertices are in CW order (signed area = -2). The optix_runner's + // cross-product sign check detects this and swaps p1<->p2 to produce + // CCW winding before building the geometry. const double x1 = 0.0, x2 = 1.0, x3 = 2.0 * x2; const double y1 = 0.0, y2 = 2.0, y3 = y1; + const double ROT_DEG = 110.0; auto surf = make_surface(); auto aper = make_aperture(x1, y1, x2, y2, x3, y3); + uint_fast64_t fpos = 0, fneg = 0, hits = 0, misses = 0; + SimulationData sd; - element_id test_elid = set_default_sd(sd, surf, aper); + element_id test_elid = set_default_sd(sd, surf, aper, ROT_DEG); SimulationResult result; OptixRunner runner; @@ -185,6 +349,8 @@ TEST(OptixRunner, FlatTriangle) ASSERT_EQ(result.get_number_of_records(), sd.get_simulation_parameters().number_of_rays); + const double cos_rot = cos(ROT_DEG * D2R); + const double sin_rot = sin(ROT_DEG * D2R); for (int i = 0; i < (int)result.get_number_of_records(); ++i) { auto rr = result[i]; @@ -195,18 +361,31 @@ TEST(OptixRunner, FlatTriangle) auto id = rr->get_element(1); EXPECT_NEAR(p0[0], p1[0], TOL) << "ray " << i; EXPECT_NEAR(p0[1], p1[1], TOL) << "ray " << i; + const double lx = p1[0] * cos_rot - p1[1] * sin_rot; + const double ly = p1[0] * sin_rot + p1[1] * cos_rot; if (id == test_elid) { EXPECT_NEAR(p1[2], Z_ELEM, TOL * Z_ELEM) << "ray " << i; - EXPECT_TRUE(aper->is_in(p1[0], p1[1])); + EXPECT_TRUE(aper->is_in(lx, ly)); + ++hits; + if (!aper->is_in(lx, ly)) + ++fpos; } else { EXPECT_NEAR(p1[2], Z_BACKSTOP, TOL * Z_ELEM); - EXPECT_FALSE(aper->is_in(p1[0], p1[1])); + EXPECT_FALSE(aper->is_in(lx, ly)); + ++misses; + if (aper->is_in(lx, ly)) + ++fneg; } } + EXPECT_GT(hits, 0u); + EXPECT_GT(misses, 0u); + std::cout << "hits: " << hits << ", misses: " << misses + << ", false positives: " << fpos + << ", false negatives: " << fneg << std::endl; } TEST(OptixRunner, FlatQuadrilateral) @@ -214,12 +393,15 @@ TEST(OptixRunner, FlatQuadrilateral) // Parallelogram const double x1 = 0.0, x2 = 3.0, x3 = (x2 - x1) + 1.0, x4 = x3 - x2 + x1; const double y1 = 0.0, y2 = y1, y3 = 2.0, y4 = y3; + const double ROT_DEG = -45.0; auto surf = make_surface(); auto aper = make_aperture( x1, y1, x2, y2, x3, y3, x4, y4); + uint_fast64_t fpos = 0, fneg = 0, hits = 0, misses = 0; + SimulationData sd; - element_id test_elid = set_default_sd(sd, surf, aper); + element_id test_elid = set_default_sd(sd, surf, aper, ROT_DEG); SimulationResult result; OptixRunner runner; @@ -234,6 +416,8 @@ TEST(OptixRunner, FlatQuadrilateral) ASSERT_EQ(result.get_number_of_records(), sd.get_simulation_parameters().number_of_rays); + const double cos_rot = cos(ROT_DEG * D2R); + const double sin_rot = sin(ROT_DEG * D2R); for (int i = 0; i < (int)result.get_number_of_records(); ++i) { auto rr = result[i]; @@ -244,18 +428,104 @@ TEST(OptixRunner, FlatQuadrilateral) auto id = rr->get_element(1); EXPECT_NEAR(p0[0], p1[0], TOL) << "ray " << i; EXPECT_NEAR(p0[1], p1[1], TOL) << "ray " << i; + const double lx = p1[0] * cos_rot - p1[1] * sin_rot; + const double ly = p1[0] * sin_rot + p1[1] * cos_rot; if (id == test_elid) { EXPECT_NEAR(p1[2], Z_ELEM, TOL * Z_ELEM) << "ray " << i; - EXPECT_TRUE(aper->is_in(p1[0], p1[1])); + EXPECT_TRUE(aper->is_in(lx, ly)); + ++hits; + if (!aper->is_in(lx, ly)) + ++fpos; } else { EXPECT_NEAR(p1[2], Z_BACKSTOP, TOL * Z_ELEM); - EXPECT_FALSE(aper->is_in(p1[0], p1[1])); + EXPECT_FALSE(aper->is_in(lx, ly)); + ++misses; + if (aper->is_in(lx, ly)) + ++fneg; + } + } + EXPECT_GT(hits, 0u); + EXPECT_GT(misses, 0u); + std::cout << "hits: " << hits << ", misses: " << misses + << ", false positives: " << fpos + << ", false negatives: " << fneg << std::endl; +} + +TEST(OptixRunner, FlatQuadrilateral_CW) +{ + // Same parallelogram as FlatQuadrilateral but vertices supplied in + // clockwise winding order: (0,0)->(1,2)->(4,2)->(3,0). + // The x1-x3 diagonal is already interior for this ordering, so + // ensure_valid_diagonal() leaves the vertices unchanged (CW). + // The optix_runner's shoelace sign check detects the negative area and + // swaps p1<->p3 to restore CCW winding before building the geometry. + // Results should be identical to FlatQuadrilateral. + const double x1 = 0.0, x2 = 1.0, x3 = 4.0, x4 = 3.0; + const double y1 = 0.0, y2 = 2.0, y3 = 2.0, y4 = 0.0; + const double ROT_DEG = -45.0; + auto surf = make_surface(); + auto aper = make_aperture( + x1, y1, x2, y2, x3, y3, x4, y4); + + uint_fast64_t fpos = 0, fneg = 0, hits = 0, misses = 0; + + SimulationData sd; + element_id test_elid = set_default_sd(sd, surf, aper, ROT_DEG); + SimulationResult result; + + OptixRunner runner; + RunnerStatus sts = runner.initialize(); + ASSERT_EQ(sts, RunnerStatus::SUCCESS); + sts = runner.setup_simulation(&sd); + ASSERT_EQ(sts, RunnerStatus::SUCCESS); + sts = runner.run_simulation(); + ASSERT_EQ(sts, RunnerStatus::SUCCESS); + sts = runner.report_simulation(&result, 0); + ASSERT_EQ(sts, RunnerStatus::SUCCESS); + + ASSERT_EQ(result.get_number_of_records(), + sd.get_simulation_parameters().number_of_rays); + const double cos_rot = cos(ROT_DEG * D2R); + const double sin_rot = sin(ROT_DEG * D2R); + for (int i = 0; i < (int)result.get_number_of_records(); ++i) + { + auto rr = result[i]; + ASSERT_GE(rr->get_number_of_interactions(), 2); + glm::dvec3 p0, p1; + rr->get_position(0, p0); + rr->get_position(1, p1); + auto id = rr->get_element(1); + EXPECT_NEAR(p0[0], p1[0], TOL) << "ray " << i; + EXPECT_NEAR(p0[1], p1[1], TOL) << "ray " << i; + const double lx = p1[0] * cos_rot - p1[1] * sin_rot; + const double ly = p1[0] * sin_rot + p1[1] * cos_rot; + + if (id == test_elid) + { + EXPECT_NEAR(p1[2], Z_ELEM, TOL * Z_ELEM) << "ray " << i; + EXPECT_TRUE(aper->is_in(lx, ly)); + ++hits; + if (!aper->is_in(lx, ly)) + ++fpos; + } + else + { + EXPECT_NEAR(p1[2], Z_BACKSTOP, TOL * Z_ELEM); + EXPECT_FALSE(aper->is_in(lx, ly)); + ++misses; + if (aper->is_in(lx, ly)) + ++fneg; } } + EXPECT_GT(hits, 0u); + EXPECT_GT(misses, 0u); + std::cout << "hits: " << hits << ", misses: " << misses + << ", false positives: " << fpos + << ", false negatives: " << fneg << std::endl; } TEST(OptixRunner, ParabolaRectangle) @@ -265,11 +535,14 @@ TEST(OptixRunner, ParabolaRectangle) constexpr double FX = 0.5 / CX; constexpr double FY = 0.5 / CY; const double XL = 10.0, YL = 5.0; + const double ROT_DEG = -135.0; auto surf = make_surface(FX, FY); auto aper = make_aperture(XL, YL); + uint_fast64_t fpos = 0, fneg = 0, hits = 0, misses = 0; + SimulationData sd; - element_id test_elid = set_default_sd(sd, surf, aper); + element_id test_elid = set_default_sd(sd, surf, aper, ROT_DEG); SimulationResult result; OptixRunner runner; @@ -284,6 +557,8 @@ TEST(OptixRunner, ParabolaRectangle) ASSERT_EQ(result.get_number_of_records(), sd.get_simulation_parameters().number_of_rays); + const double cos_rot = cos(ROT_DEG * D2R); + const double sin_rot = sin(ROT_DEG * D2R); for (int i = 0; i < (int)result.get_number_of_records(); ++i) { auto rr = result[i]; @@ -294,30 +569,46 @@ TEST(OptixRunner, ParabolaRectangle) auto id = rr->get_element(1); EXPECT_NEAR(p0[0], p1[0], TOL) << "ray " << i; EXPECT_NEAR(p0[1], p1[1], TOL) << "ray " << i; + const double lx = p1[0] * cos_rot - p1[1] * sin_rot; + const double ly = p1[0] * sin_rot + p1[1] * cos_rot; if (id == test_elid) { - const double z1 = Z_ELEM + 0.5 * CX * p1[0] * p1[0] + 0.5 * CY * p1[1] * p1[1]; + const double z1 = Z_ELEM + 0.5 * CX * lx * lx + 0.5 * CY * ly * ly; EXPECT_NEAR(p1[2], z1, TOL * Z_ELEM) << "ray " << i; - EXPECT_TRUE(aper->is_in(p1[0], p1[1])); + EXPECT_TRUE(aper->is_in(lx, ly)); + ++hits; + if (!aper->is_in(lx, ly)) + ++fpos; } else { EXPECT_NEAR(p1[2], Z_BACKSTOP, TOL * Z_ELEM); - EXPECT_FALSE(aper->is_in(p1[0], p1[1])); + EXPECT_FALSE(aper->is_in(lx, ly)); + ++misses; + if (aper->is_in(lx, ly)) + ++fneg; } } + EXPECT_GT(hits, 0u); + EXPECT_GT(misses, 0u); + std::cout << "hits: " << hits << ", misses: " << misses + << ", false positives: " << fpos + << ", false negatives: " << fneg << std::endl; } TEST(OptixRunner, Cylinder) { const double R = 5.0; const double YL = 3.0; // Total cylinder length + const double ROT_DEG = 25.0; auto surf = make_surface(R); auto aper = make_aperture(2 * R, YL); + uint_fast64_t fpos = 0, fneg = 0, hits = 0, misses = 0; + SimulationData sd; - element_id test_elid = set_default_sd(sd, surf, aper); + element_id test_elid = set_default_sd(sd, surf, aper, ROT_DEG); SimulationResult result; OptixRunner runner; @@ -332,6 +623,8 @@ TEST(OptixRunner, Cylinder) ASSERT_EQ(result.get_number_of_records(), sd.get_simulation_parameters().number_of_rays); + const double cos_rot = cos(ROT_DEG * D2R); + const double sin_rot = sin(ROT_DEG * D2R); for (int i = 0; i < (int)result.get_number_of_records(); ++i) { auto rr = result[i]; @@ -342,29 +635,45 @@ TEST(OptixRunner, Cylinder) auto id = rr->get_element(1); EXPECT_NEAR(p0[0], p1[0], TOL) << "ray " << i; EXPECT_NEAR(p0[1], p1[1], TOL) << "ray " << i; + const double lx = p1[0] * cos_rot - p1[1] * sin_rot; + const double ly = p1[0] * sin_rot + p1[1] * cos_rot; if (id == test_elid) { - const double z1 = Z_ELEM + sqrt(R * R - p1[0] * p1[0]); + const double z1 = Z_ELEM + sqrt(R * R - lx * lx); EXPECT_NEAR(p1[2], z1, TOL * Z_ELEM) << "ray " << i; - EXPECT_TRUE(aper->is_in(p1[0], p1[1])); + EXPECT_TRUE(aper->is_in(lx, ly)); + ++hits; + if (!aper->is_in(lx, ly)) + ++fpos; } else { EXPECT_NEAR(p1[2], Z_BACKSTOP, TOL * Z_ELEM); - EXPECT_FALSE(aper->is_in(p1[0], p1[1])); + EXPECT_FALSE(aper->is_in(lx, ly)); + ++misses; + if (aper->is_in(lx, ly)) + ++fneg; } } + EXPECT_GT(hits, 0u); + EXPECT_GT(misses, 0u); + std::cout << "hits: " << hits << ", misses: " << misses + << ", false positives: " << fpos + << ", false negatives: " << fneg << std::endl; } TEST(OptixRunner, FlatCircle) { const double R = 5.0; + const double ROT_DEG = 10.0; // Should make no difference auto surf = make_surface(); auto aper = make_aperture(2 * R); + uint_fast64_t fpos = 0, fneg = 0, hits = 0, misses = 0; + SimulationData sd; - element_id test_elid = set_default_sd(sd, surf, aper); + element_id test_elid = set_default_sd(sd, surf, aper, ROT_DEG); SimulationResult result; OptixRunner runner; @@ -394,23 +703,37 @@ TEST(OptixRunner, FlatCircle) { EXPECT_NEAR(p1[2], Z_ELEM, TOL * Z_ELEM) << "ray " << i; EXPECT_TRUE(aper->is_in(p1[0], p1[1])); + ++hits; + if (!aper->is_in(p1[0], p1[1])) + ++fpos; } else { EXPECT_NEAR(p1[2], Z_BACKSTOP, TOL * Z_ELEM); EXPECT_FALSE(aper->is_in(p1[0], p1[1])); + ++misses; + if (aper->is_in(p1[0], p1[1])) + ++fneg; } } + EXPECT_GT(hits, 0u); + EXPECT_GT(misses, 0u); + std::cout << "hits: " << hits << ", misses: " << misses + << ", false positives: " << fpos + << ", false negatives: " << fneg << std::endl; } TEST(OptixRunner, FlatHexagon) { const double S = 5.0; + const double ROT_DEG = 30.0; auto surf = make_surface(); auto aper = make_aperture(2 * S); + uint_fast64_t fpos = 0, fneg = 0, hits = 0, misses = 0; + SimulationData sd; - element_id test_elid = set_default_sd(sd, surf, aper); + element_id test_elid = set_default_sd(sd, surf, aper, ROT_DEG); SimulationResult result; OptixRunner runner; @@ -425,6 +748,8 @@ TEST(OptixRunner, FlatHexagon) ASSERT_EQ(result.get_number_of_records(), sd.get_simulation_parameters().number_of_rays); + const double cos_rot = cos(ROT_DEG * D2R); + const double sin_rot = sin(ROT_DEG * D2R); for (int i = 0; i < (int)result.get_number_of_records(); ++i) { auto rr = result[i]; @@ -435,30 +760,46 @@ TEST(OptixRunner, FlatHexagon) auto id = rr->get_element(1); EXPECT_NEAR(p0[0], p1[0], TOL) << "ray " << i; EXPECT_NEAR(p0[1], p1[1], TOL) << "ray " << i; + const double lx = p1[0] * cos_rot - p1[1] * sin_rot; + const double ly = p1[0] * sin_rot + p1[1] * cos_rot; if (id == test_elid) { EXPECT_NEAR(p1[2], Z_ELEM, TOL * Z_ELEM) << "ray " << i; - EXPECT_TRUE(aper->is_in(p1[0], p1[1])); + EXPECT_TRUE(aper->is_in(lx, ly)); + ++hits; + if (!aper->is_in(lx, ly)) + ++fpos; } else { EXPECT_NEAR(p1[2], Z_BACKSTOP, TOL * Z_ELEM); - EXPECT_FALSE(aper->is_in(p1[0], p1[1])); + EXPECT_FALSE(aper->is_in(lx, ly)); + ++misses; + if (aper->is_in(lx, ly)) + ++fneg; } } + EXPECT_GT(hits, 0u); + EXPECT_GT(misses, 0u); + std::cout << "hits: " << hits << ", misses: " << misses + << ", false positives: " << fpos + << ", false negatives: " << fneg << std::endl; } -TEST(OptixRunner, FlatAnnulus) +TEST(OptixRunner, FlatAnnulus_FullArc) { const double R0 = 5.0; - const double R1 = 10.0; - const double ARC = 2 * PI; + const double R1 = 180.0; + const double ARC = 360.0; + const double ROT_DEG = -15.0; // Should make no difference auto surf = make_surface(); auto aper = make_aperture(R0, R1, ARC); + uint_fast64_t fpos = 0, fneg = 0, hits = 0, misses = 0; + SimulationData sd; - element_id test_elid = set_default_sd(sd, surf, aper); + element_id test_elid = set_default_sd(sd, surf, aper, ROT_DEG); SimulationResult result; OptixRunner runner; @@ -488,11 +829,449 @@ TEST(OptixRunner, FlatAnnulus) { EXPECT_NEAR(p1[2], Z_ELEM, TOL * Z_ELEM) << "ray " << i; EXPECT_TRUE(aper->is_in(p1[0], p1[1])); + ++hits; + if (!aper->is_in(p1[0], p1[1])) + ++fpos; } else { EXPECT_NEAR(p1[2], Z_BACKSTOP, TOL * Z_ELEM); EXPECT_FALSE(aper->is_in(p1[0], p1[1])); + ++misses; + if (aper->is_in(p1[0], p1[1])) + ++fneg; + } + } + EXPECT_GT(hits, 0u); + EXPECT_GT(misses, 0u); + std::cout << "hits: " << hits << ", misses: " << misses + << ", false positives: " << fpos + << ", false negatives: " << fneg << std::endl; +} + +TEST(OptixRunner, FlatAnnulus_PartialArc) +{ + const double R0 = 5.0; + const double R1 = 180.0; + const double ARC = 90.0; + const double ROT_DEG = -15.0; + auto surf = make_surface(); + auto aper = make_aperture(R0, R1, ARC); + + uint_fast64_t fpos = 0, fneg = 0, hits = 0, misses = 0; + + SimulationData sd; + element_id test_elid = set_default_sd(sd, surf, aper, ROT_DEG); + SimulationResult result; + + OptixRunner runner; + RunnerStatus sts = runner.initialize(); + ASSERT_EQ(sts, RunnerStatus::SUCCESS); + sts = runner.setup_simulation(&sd); + ASSERT_EQ(sts, RunnerStatus::SUCCESS); + sts = runner.run_simulation(); + ASSERT_EQ(sts, RunnerStatus::SUCCESS); + sts = runner.report_simulation(&result, 0); + ASSERT_EQ(sts, RunnerStatus::SUCCESS); + + ASSERT_EQ(result.get_number_of_records(), + sd.get_simulation_parameters().number_of_rays); + const double cos_rot = cos(ROT_DEG * D2R); + const double sin_rot = sin(ROT_DEG * D2R); + for (int i = 0; i < (int)result.get_number_of_records(); ++i) + { + auto rr = result[i]; + ASSERT_GE(rr->get_number_of_interactions(), 2); + glm::dvec3 p0, p1; + rr->get_position(0, p0); + rr->get_position(1, p1); + auto id = rr->get_element(1); + EXPECT_NEAR(p0[0], p1[0], TOL) << "ray " << i; + EXPECT_NEAR(p0[1], p1[1], TOL) << "ray " << i; + const double lx = p1[0] * cos_rot - p1[1] * sin_rot; + const double ly = p1[0] * sin_rot + p1[1] * cos_rot; + + if (id == test_elid) + { + EXPECT_NEAR(p1[2], Z_ELEM, TOL * Z_ELEM) << "ray " << i; + EXPECT_TRUE(aper->is_in(lx, ly)); + ++hits; + if (!aper->is_in(lx, ly)) + ++fpos; + } + else + { + EXPECT_NEAR(p1[2], Z_BACKSTOP, TOL * Z_ELEM); + EXPECT_FALSE(aper->is_in(lx, ly)); + ++misses; + if (aper->is_in(lx, ly)) + ++fneg; + } + } + EXPECT_GT(hits, 0u); + EXPECT_GT(misses, 0u); + std::cout << "hits: " << hits << ", misses: " << misses + << ", false positives: " << fpos + << ", false negatives: " << fneg << std::endl; +} + +TEST(OptixRunner, ParabolicHexagon) +{ + const double S = 5.0; + const double FOCAL_X = 10.0; + const double FOCAL_Y = 30.0; + const double ROT_DEG = 30.0; + auto surf = make_surface(FOCAL_X, FOCAL_Y); + auto aper = make_aperture(2 * S); + + uint_fast64_t fpos = 0, fneg = 0, hits = 0, misses = 0; + + SimulationData sd; + element_id test_elid = set_default_sd(sd, surf, aper, ROT_DEG); + SimulationResult result; + + OptixRunner runner; + RunnerStatus sts = runner.initialize(); + ASSERT_EQ(sts, RunnerStatus::SUCCESS); + sts = runner.setup_simulation(&sd); + ASSERT_EQ(sts, RunnerStatus::SUCCESS); + sts = runner.run_simulation(); + ASSERT_EQ(sts, RunnerStatus::SUCCESS); + sts = runner.report_simulation(&result, 0); + ASSERT_EQ(sts, RunnerStatus::SUCCESS); + + ASSERT_EQ(result.get_number_of_records(), + sd.get_simulation_parameters().number_of_rays); + + const double cos_rot = cos(ROT_DEG * D2R); + const double sin_rot = sin(ROT_DEG * D2R); + + for (int i = 0; i < (int)result.get_number_of_records(); ++i) + { + auto rr = result[i]; + ASSERT_GE(rr->get_number_of_interactions(), 2); + glm::dvec3 p0, p1; + rr->get_position(0, p0); + rr->get_position(1, p1); + auto id = rr->get_element(1); + EXPECT_NEAR(p0[0], p1[0], TOL) << "ray " << i; + EXPECT_NEAR(p0[1], p1[1], TOL) << "ray " << i; + + const double lx = p1[0] * cos_rot - p1[1] * sin_rot; + const double ly = p1[0] * sin_rot + p1[1] * cos_rot; + + if (id == test_elid) + { + double zsol = lx * lx / (4.0 * FOCAL_X) + ly * ly / (4.0 * FOCAL_Y) + Z_ELEM; + EXPECT_NEAR(p1[2], zsol, TOL * Z_ELEM) << "ray " << i; + EXPECT_TRUE(aper->is_in(lx, ly)); + ++hits; + if (!aper->is_in(lx, ly)) + ++fpos; + } + else + { + EXPECT_NEAR(p1[2], Z_BACKSTOP, TOL * Z_ELEM); + EXPECT_FALSE(aper->is_in(lx, ly)); + ++misses; + if (aper->is_in(lx, ly)) + ++fneg; + } + } + EXPECT_GT(hits, 0u); + EXPECT_GT(misses, 0u); + std::cout << "hits: " << hits << ", misses: " << misses + << ", false positives: " << fpos + << ", false negatives: " << fneg << std::endl; +} + +TEST(OptixRunner, ParabolicTriangle) +{ + const double d = 4.0; + const double FOCAL_X = 10.0; + const double FOCAL_Y = 30.0; + const double ROT_DEG = 90.0; + auto surf = make_surface(FOCAL_X, FOCAL_Y); + auto aper = make_aperture(d); + + uint_fast64_t fpos = 0, fneg = 0, hits = 0, misses = 0; + + SimulationData sd; + element_id test_elid = set_default_sd(sd, surf, aper, ROT_DEG); + SimulationResult result; + + OptixRunner runner; + RunnerStatus sts = runner.initialize(); + ASSERT_EQ(sts, RunnerStatus::SUCCESS); + sts = runner.setup_simulation(&sd); + ASSERT_EQ(sts, RunnerStatus::SUCCESS); + sts = runner.run_simulation(); + ASSERT_EQ(sts, RunnerStatus::SUCCESS); + sts = runner.report_simulation(&result, 0); + ASSERT_EQ(sts, RunnerStatus::SUCCESS); + + ASSERT_EQ(result.get_number_of_records(), + sd.get_simulation_parameters().number_of_rays); + + const double cos_rot = cos(ROT_DEG * D2R); + const double sin_rot = sin(ROT_DEG * D2R); + + for (int i = 0; i < (int)result.get_number_of_records(); ++i) + { + auto rr = result[i]; + ASSERT_GE(rr->get_number_of_interactions(), 2); + glm::dvec3 p0, p1; + rr->get_position(0, p0); + rr->get_position(1, p1); + auto id = rr->get_element(1); + EXPECT_NEAR(p0[0], p1[0], TOL) << "ray " << i; + EXPECT_NEAR(p0[1], p1[1], TOL) << "ray " << i; + + const double lx = p1[0] * cos_rot - p1[1] * sin_rot; + const double ly = p1[0] * sin_rot + p1[1] * cos_rot; + + if (id == test_elid) + { + double zsol = lx * lx / (4.0 * FOCAL_X) + ly * ly / (4.0 * FOCAL_Y) + Z_ELEM; + EXPECT_NEAR(p1[2], zsol, TOL * Z_ELEM) << "ray " << i; + EXPECT_TRUE(aper->is_in(lx, ly)); + ++hits; + if (!aper->is_in(lx, ly)) + ++fpos; + } + else + { + EXPECT_NEAR(p1[2], Z_BACKSTOP, TOL * Z_ELEM); + EXPECT_FALSE(aper->is_in(lx, ly)); + ++misses; + if (aper->is_in(lx, ly)) + ++fneg; + } + } + EXPECT_GT(hits, 0u); + EXPECT_GT(misses, 0u); + std::cout << "hits: " << hits << ", misses: " << misses + << ", false positives: " << fpos + << ", false negatives: " << fneg << std::endl; +} + +TEST(OptixRunner, ParabolicAnnulus) +{ + const double R0 = 2.0; + const double R1 = 5.0; + const double ARC = 180.0; + const double FOCAL_X = 10.0; + const double FOCAL_Y = 30.0; + const double ROT_DEG = -15.0; + auto surf = make_surface(FOCAL_X, FOCAL_Y); + auto aper = make_aperture(R0, R1, ARC); + + uint_fast64_t fpos = 0, fneg = 0, hits = 0, misses = 0; + + SimulationData sd; + element_id test_elid = set_default_sd(sd, surf, aper, ROT_DEG); + SimulationResult result; + + OptixRunner runner; + RunnerStatus sts = runner.initialize(); + ASSERT_EQ(sts, RunnerStatus::SUCCESS); + sts = runner.setup_simulation(&sd); + ASSERT_EQ(sts, RunnerStatus::SUCCESS); + sts = runner.run_simulation(); + ASSERT_EQ(sts, RunnerStatus::SUCCESS); + sts = runner.report_simulation(&result, 0); + ASSERT_EQ(sts, RunnerStatus::SUCCESS); + + ASSERT_EQ(result.get_number_of_records(), + sd.get_simulation_parameters().number_of_rays); + + const double cos_rot = cos(ROT_DEG * D2R); + const double sin_rot = sin(ROT_DEG * D2R); + + for (int i = 0; i < (int)result.get_number_of_records(); ++i) + { + auto rr = result[i]; + ASSERT_GE(rr->get_number_of_interactions(), 2); + glm::dvec3 p0, p1; + rr->get_position(0, p0); + rr->get_position(1, p1); + auto id = rr->get_element(1); + EXPECT_NEAR(p0[0], p1[0], TOL) << "ray " << i; + EXPECT_NEAR(p0[1], p1[1], TOL) << "ray " << i; + + const double lx = p1[0] * cos_rot - p1[1] * sin_rot; + const double ly = p1[0] * sin_rot + p1[1] * cos_rot; + + if (id == test_elid) + { + double zsol = lx * lx / (4.0 * FOCAL_X) + ly * ly / (4.0 * FOCAL_Y) + Z_ELEM; + EXPECT_NEAR(p1[2], zsol, TOL * Z_ELEM) << "ray " << i; + EXPECT_TRUE(aper->is_in(lx, ly)); + ++hits; + if (!aper->is_in(lx, ly)) + ++fpos; + } + else + { + EXPECT_NEAR(p1[2], Z_BACKSTOP, TOL * Z_ELEM); + EXPECT_FALSE(aper->is_in(lx, ly)); + ++misses; + if (aper->is_in(lx, ly)) + ++fneg; + } + } + EXPECT_GT(hits, 0u); + EXPECT_GT(misses, 0u); + std::cout << "hits: " << hits << ", misses: " << misses + << ", false positives: " << fpos + << ", false negatives: " << fneg << std::endl; +} + +TEST(OptixRunner, ParabolicQuadrilateral) +{ + // Parallelogram + const double x1 = 0.0, x2 = 3.0, x3 = (x2 - x1) + 1.0, x4 = x3 - x2 + x1; + const double y1 = 0.0, y2 = y1, y3 = 2.0, y4 = y3; + const double FOCAL_X = 10.0; + const double FOCAL_Y = 30.0; + const double ROT_DEG = -45.0; + auto surf = make_surface(FOCAL_X, FOCAL_Y); + auto aper = make_aperture(x1, y1, x2, y2, x3, y3, x4, y4); + + uint_fast64_t fpos = 0, fneg = 0, hits = 0, misses = 0; + + SimulationData sd; + element_id test_elid = set_default_sd(sd, surf, aper, ROT_DEG); + SimulationResult result; + + OptixRunner runner; + RunnerStatus sts = runner.initialize(); + ASSERT_EQ(sts, RunnerStatus::SUCCESS); + sts = runner.setup_simulation(&sd); + ASSERT_EQ(sts, RunnerStatus::SUCCESS); + sts = runner.run_simulation(); + ASSERT_EQ(sts, RunnerStatus::SUCCESS); + sts = runner.report_simulation(&result, 0); + ASSERT_EQ(sts, RunnerStatus::SUCCESS); + + ASSERT_EQ(result.get_number_of_records(), + sd.get_simulation_parameters().number_of_rays); + + const double cos_rot = cos(ROT_DEG * D2R); + const double sin_rot = sin(ROT_DEG * D2R); + + for (int i = 0; i < (int)result.get_number_of_records(); ++i) + { + auto rr = result[i]; + ASSERT_GE(rr->get_number_of_interactions(), 2); + glm::dvec3 p0, p1; + rr->get_position(0, p0); + rr->get_position(1, p1); + auto id = rr->get_element(1); + EXPECT_NEAR(p0[0], p1[0], TOL) << "ray " << i; + EXPECT_NEAR(p0[1], p1[1], TOL) << "ray " << i; + + const double lx = p1[0] * cos_rot - p1[1] * sin_rot; + const double ly = p1[0] * sin_rot + p1[1] * cos_rot; + + if (id == test_elid) + { + double zsol = lx * lx / (4.0 * FOCAL_X) + ly * ly / (4.0 * FOCAL_Y) + Z_ELEM; + EXPECT_NEAR(p1[2], zsol, TOL * Z_ELEM) << "ray " << i; + EXPECT_TRUE(aper->is_in(lx, ly)); + ++hits; + if (!aper->is_in(lx, ly)) + ++fpos; + } + else + { + EXPECT_NEAR(p1[2], Z_BACKSTOP, TOL * Z_ELEM); + EXPECT_FALSE(aper->is_in(lx, ly)); + ++misses; + if (aper->is_in(lx, ly)) + ++fneg; + } + } + EXPECT_GT(hits, 0u); + EXPECT_GT(misses, 0u); + std::cout << "hits: " << hits << ", misses: " << misses + << ", false positives: " << fpos + << ", false negatives: " << fneg << std::endl; +} + +TEST(OptixRunner, ParabolicCircle) +{ + const double R = 5.0; + const double FOCAL_X = 10.0; + const double FOCAL_Y = 30.0; + const double ROT_DEG = -45.0; + auto surf = make_surface(FOCAL_X, FOCAL_Y); + auto aper = make_aperture(2 * R); + + // z(x,y) = x^2/(2*FOCAL_X) + y^2/(2*FOCAL_Y). Maximum over the circle x^2+y^2<=R^2 + // is R^2 / (2 * min(FOCAL_X, FOCAL_Y)), achieved along the axis with lower focal length. + // const double Z_MAX_OFFSET = R * R / (2.0 * std::min(FOCAL_X, FOCAL_Y)); + + uint_fast64_t fpos = 0, fneg = 0, hits = 0, misses = 0; + + SimulationData sd; + element_id test_elid = set_default_sd(sd, surf, aper, ROT_DEG); + SimulationResult result; + + OptixRunner runner; + RunnerStatus sts = runner.initialize(); + ASSERT_EQ(sts, RunnerStatus::SUCCESS); + sts = runner.setup_simulation(&sd); + ASSERT_EQ(sts, RunnerStatus::SUCCESS); + sts = runner.run_simulation(); + ASSERT_EQ(sts, RunnerStatus::SUCCESS); + sts = runner.report_simulation(&result, 0); + ASSERT_EQ(sts, RunnerStatus::SUCCESS); + + ASSERT_EQ(result.get_number_of_records(), + sd.get_simulation_parameters().number_of_rays); + + const double cos_rot = cos(ROT_DEG * D2R); + const double sin_rot = sin(ROT_DEG * D2R); + + for (int i = 0; i < (int)result.get_number_of_records(); ++i) + { + auto rr = result[i]; + ASSERT_GE(rr->get_number_of_interactions(), 2); + glm::dvec3 p0, p1; + rr->get_position(0, p0); + rr->get_position(1, p1); + auto id = rr->get_element(1); + EXPECT_NEAR(p0[0], p1[0], TOL) << "ray " << i; + EXPECT_NEAR(p0[1], p1[1], TOL) << "ray " << i; + + const double lx = p1[0] * cos_rot - p1[1] * sin_rot; + const double ly = p1[0] * sin_rot + p1[1] * cos_rot; + + if (id == test_elid) + { + // z is curved: Z_ELEM <= z <= Z_ELEM + Z_MAX_OFFSET + // EXPECT_GE(p1[2], Z_ELEM - TOL * Z_ELEM) << "ray " << i; + // EXPECT_LE(p1[2], Z_ELEM + Z_MAX_OFFSET + TOL * Z_ELEM) << "ray " << i; + double zsol = lx * lx / (4.0 * FOCAL_X) + ly * ly / (4.0 * FOCAL_Y) + Z_ELEM; + EXPECT_NEAR(p1[2], zsol, TOL * Z_ELEM); + EXPECT_TRUE(aper->is_in(lx, ly)); + ++hits; + if (!aper->is_in(lx, ly)) + ++fpos; + } + else + { + EXPECT_NEAR(p1[2], Z_BACKSTOP, TOL * Z_ELEM); + EXPECT_FALSE(aper->is_in(lx, ly)); + ++misses; + if (aper->is_in(lx, ly)) + ++fneg; } } + EXPECT_GT(hits, 0u); + EXPECT_GT(misses, 0u); + std::cout << "hits: " << hits << ", misses: " << misses + << ", false positives: " << fpos + << ", false negatives: " << fneg << std::endl; }