diff --git a/src/slic3r/GUI/CAD/DesignSketchTool.cpp b/src/slic3r/GUI/CAD/DesignSketchTool.cpp index d5f0399d6e..525a396bea 100644 --- a/src/slic3r/GUI/CAD/DesignSketchTool.cpp +++ b/src/slic3r/GUI/CAD/DesignSketchTool.cpp @@ -53,6 +53,7 @@ #include #include #include +#include #include "libslic3r/AppConfig.hpp" #include "libslic3r/Color.hpp" #include "slic3r/GUI/GLSelectionRectangle.hpp" @@ -5384,6 +5385,144 @@ void DesignSketchTool::drag_rib_handle(GLCanvas3D& canvas, const wxMouseEvent& e } // ---- Reference/base planes (Onshape-style default planes) ----------------------------- +namespace { +// One level of a BSP tree whose splitters are the planes themselves, in order: the painter's +// algorithm made exact for polygons that cross. Pieces of plane k (and of any plane lying in it) +// are in the splitter; every other piece is wholly on one side or is cut in two. Far side, then the +// splitter's own pieces, then the near side is back to front. A piece is only cut when it has +// corners strictly on both sides, so neither half can be degenerate. +void paint_back_to_front(std::vector&& in, size_t k, const std::vector& planes, double eps, + const Vec3d& eye, const Vec3d& forward, bool perspective, std::vector& out) +{ + if (in.empty()) + return; + if (k == planes.size()) { // not reached: every piece is in the splitter at its own plane's level + std::move(in.begin(), in.end(), std::back_inserter(out)); + return; + } + // The square's own normal: SketchPlane::normal is reversed on XZ, and the side tests only need + // one consistent choice. + const Vec3d n = planes[k].x_axis.cross(planes[k].y_axis).normalized(); + const double d = n.dot(planes[k].origin); + std::vector front, back, on; + for (PlanePiece& piece : in) { + if (piece.plane == int(k)) { + on.push_back(std::move(piece)); + continue; + } + std::vector s; + int sides = 0; + for (const Vec3d& c : piece.corners) { + s.push_back(n.dot(c) - d); + sides |= s.back() > eps ? 1 : s.back() < -eps ? 2 : 0; + } + if (sides == 0) + on.push_back(std::move(piece)); + else if (sides == 1) + front.push_back(std::move(piece)); + else if (sides == 2) + back.push_back(std::move(piece)); + else { + PlanePiece f{ piece.plane, {} }, b{ piece.plane, {} }; + const size_t m = piece.corners.size(); + for (size_t i = 0; i < m; ++i) { + const size_t j = (i + 1) % m; + const Vec3d& a = piece.corners[i]; + if (s[i] >= -eps) f.corners.push_back(a); + if (s[i] <= eps) b.corners.push_back(a); + if ((s[i] > eps && s[j] < -eps) || (s[i] < -eps && s[j] > eps)) { + const Vec3d x = a + (piece.corners[j] - a) * (s[i] / (s[i] - s[j])); + f.corners.push_back(x); + b.corners.push_back(x); + } + } + front.push_back(std::move(f)); + back.push_back(std::move(b)); + } + } + // An orthographic eye is at infinity behind the view direction. Camera::get_position() is a + // finite point there and can sit on the wrong side of a plane, so only the direction counts. + const bool eye_in_front = perspective ? n.dot(eye) - d > 0. : n.dot(forward) < 0.; + paint_back_to_front(std::move(eye_in_front ? back : front), k + 1, planes, eps, eye, forward, perspective, out); + std::move(on.begin(), on.end(), std::back_inserter(out)); + paint_back_to_front(std::move(eye_in_front ? front : back), k + 1, planes, eps, eye, forward, perspective, out); +} + +Vec2d in_frame(const SketchPlane& p, const Vec3d& x) { return Vec2d((x - p.origin).dot(p.x_axis), (x - p.origin).dot(p.y_axis)); } + +double plane_distance(const SketchPlane& p, const Vec3d& x) { return p.x_axis.cross(p.y_axis).normalized().dot(x - p.origin); } + +// Shrink the segment ab, which lies in p, to its part inside p's square (Liang-Barsky). False when +// nothing is left. +bool clip_to_square(Vec3d& a, Vec3d& b, const SketchPlane& p, double half, double eps) +{ + const Vec2d s = in_frame(p, a), d = in_frame(p, b) - s; + double t0 = 0., t1 = 1.; + for (int axis = 0; axis < 2; ++axis) + for (double sign : { -1., 1. }) { // keep sign * (s + t * d)[axis] <= half + const double room = half + eps - sign * s[axis], rate = sign * d[axis]; + if (rate > 0.) + t1 = std::min(t1, room / rate); + else if (rate < 0.) + t0 = std::max(t0, room / rate); + else if (room < 0.) + return false; + } + if (t1 - t0 < 1e-9) + return false; + const Vec3d ab = b - a; + b = a + ab * t1; + a = a + ab * t0; + return true; +} +} // namespace + +std::vector planes_back_to_front(const std::vector& planes, double half, + const Vec3d& eye, const Vec3d& forward, bool perspective) +{ + std::vector squares; + for (int i = 0; i < int(planes.size()); ++i) { + const SketchPlane& p = planes[i]; + squares.push_back({ i, { p.to_world(Vec2d(-half, -half)), p.to_world(Vec2d(half, -half)), + p.to_world(Vec2d(half, half)), p.to_world(Vec2d(-half, half)) } }); + } + const double eps = 1e-6 * std::max(half, 1.); + std::vector out; + paint_back_to_front(std::move(squares), 0, planes, eps, eye, forward, perspective, out); + + // What to outline: the piece's share of its square's border, and of where it crosses another + // square. Every other edge is a cut lying in the plane that made it, but the cut ran along that + // whole infinite plane, so it is clipped to that plane's square: a datum that never reaches a base + // plane gets no line across it. A piece is convex, so an edge with both ends on one side line of + // its square lies along that side. + for (PlanePiece& piece : out) { + const SketchPlane& own = planes[piece.plane]; + for (size_t i = 0; i < piece.corners.size(); ++i) { + const Vec3d& a = piece.corners[i]; + const Vec3d& b = piece.corners[(i + 1) % piece.corners.size()]; + const Vec2d fa = in_frame(own, a), fb = in_frame(own, b); + bool border = false; + for (int axis = 0; axis < 2; ++axis) + for (double side : { -half, half }) + border = border || (std::abs(fa[axis] - side) <= eps && std::abs(fb[axis] - side) <= eps); + if (border) { + piece.lines.emplace_back(a, b); + continue; + } + for (int k = 0; k < int(planes.size()); ++k) { + auto in_k = [&](const Vec3d& x) { return std::abs(plane_distance(planes[k], x)) <= 2. * eps; }; + // Not a plane the piece itself lies in: clipped to its own square, a cut is kept whole. + if (k == piece.plane || !in_k(a) || !in_k(b) || std::all_of(piece.corners.begin(), piece.corners.end(), in_k)) + continue; + Vec3d ca = a, cb = b; + if (clip_to_square(ca, cb, planes[k], half, eps)) + piece.lines.emplace_back(ca, cb); + } + } + } + return out; +} + void DesignSketchTool::set_base_pick(std::vector planes, std::vector bases, std::vector labels) { @@ -5419,49 +5558,72 @@ double DesignSketchTool::dbp_half_extent() const return half; } -// Draw the reference planes as large translucent labelled squares; the hovered one brightens. +// Draw the reference planes as labelled translucent squares outlined in their own hue; the hovered +// one brightens. Depth testing is off (they overlay the bed and any bodies), so draw order is the +// blend order — and the planes cross, so they go down piece by piece, back to front. void DesignSketchTool::render_base_pick() { if (!m_dbp_active || m_dbp_planes.empty()) return; using EPT = GLModel::Geometry::EPrimitiveType; using EVL = GLModel::Geometry::EVertexLayout; const double H = dbp_half_extent(); - // Onshape-ish per-plane tints: XY blue, XZ green, YZ red (keyed by base index 0/1/2; datums grey). - auto tint = [](int base, bool hot) -> ColorRGBA { - float a = hot ? 0.10f : 0.047f; // base planes kept faint (reduced ~2/3 from 0.30/0.14) - if (base == 0) return ColorRGBA(0.30f, 0.55f, 0.95f, a); - if (base == 1) return ColorRGBA(0.35f, 0.80f, 0.45f, a); - if (base == 2) return ColorRGBA(0.92f, 0.42f, 0.42f, a); - return ColorRGBA(0.70f, 0.72f, 0.78f, a); + // Two layers of alpha a blended in either order differ by only a^2 of their colour difference, + // so a faint fill hides which plane is in front however well the pieces are sorted: at 0.16 that + // is under 3%, and the crossing planes read as one grey smear. At 0.35 it is ~12%, and the bed + // grid still reads through all three. + constexpr float kFillAlpha = 0.35f, kFillAlphaHot = 0.50f; + constexpr float kLineAlpha = 0.90f, kLineAlphaHot = 1.00f; + // Onshape-ish per-plane hues: XY blue, XZ green, YZ red (keyed by base index 0/1/2; datums grey). + auto hue = [](int base) -> ColorRGBA { + if (base == 0) return ColorRGBA(0.30f, 0.55f, 0.95f, 1.0f); + if (base == 1) return ColorRGBA(0.35f, 0.80f, 0.45f, 1.0f); + if (base == 2) return ColorRGBA(0.92f, 0.42f, 0.42f, 1.0f); + return ColorRGBA(0.70f, 0.72f, 0.78f, 1.0f); }; + const Camera& cam = wxGetApp().plater()->get_camera(); + const Vec3d vd = cam.get_dir_forward(); + const double hw = 1.5 / std::max(cam.get_zoom(), 1e-6); // outline ribbon, as on datum planes glsafe(::glDisable(GL_DEPTH_TEST)); glsafe(::glDisable(GL_CULL_FACE)); glsafe(::glEnable(GL_BLEND)); // alpha is ignored without this glsafe(::glBlendFunc(GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA)); - const SketchPlane saved_plane = m_plane; - for (size_t i = 0; i < m_dbp_planes.size(); ++i) { - const SketchPlane& p = m_dbp_planes[i]; - const Vec3d q0 = p.to_world(Vec2d(-H, -H)), q1 = p.to_world(Vec2d(H, -H)), - q2 = p.to_world(Vec2d(H, H)), q3 = p.to_world(Vec2d(-H, H)); - GLModel::Geometry quad; quad.format = { EPT::Triangles, EVL::P3 }; - quad.add_vertex((Vec3f)q0.cast()); quad.add_vertex((Vec3f)q1.cast()); - quad.add_vertex((Vec3f)q2.cast()); quad.add_vertex((Vec3f)q3.cast()); - quad.add_triangle(0, 1, 2); quad.add_triangle(0, 2, 3); - GLModel m; m.init_from(std::move(quad)); - const bool hot = (int(i) == m_dbp_hover); - const int base = (i < m_dbp_base.size()) ? m_dbp_base[i] : -1; - m.set_color(tint(base, hot)); - m.render(); + for (const PlanePiece& piece : planes_back_to_front(m_dbp_planes, H, cam.get_position(), vd, + cam.get_type() == Camera::EType::Perspective)) { + const bool hot = piece.plane == m_dbp_hover; + ColorRGBA col = hue(piece.plane < int(m_dbp_base.size()) ? m_dbp_base[piece.plane] : -1); + GLModel::Geometry fill; fill.format = { EPT::Triangles, EVL::P3 }; + for (const Vec3d& q : piece.corners) fill.add_vertex((Vec3f)q.cast()); + for (unsigned int i = 1; i + 1 < piece.corners.size(); ++i) fill.add_triangle(0, i, i + 1); // convex: a fan + GLModel fm; fm.init_from(std::move(fill)); + col.a(hot ? kFillAlphaHot : kFillAlpha); + fm.set_color(col); + fm.render(); - // Label near the top-left corner, drawn in the plane (draw_text lifts through m_plane). - if (i < m_dbp_labels.size() && !m_dbp_labels[i].empty()) { - m_plane = p; - const double th = H * 0.10; - const ColorRGBA lc = tint(base, true); ColorRGBA lcs(lc.r(), lc.g(), lc.b(), 1.0f); - draw_text(m_line_model, m_dbp_labels[i], Vec2d(-H + th * 2.0, H - th * 1.6), th, lcs); + // Its outline and crossing lines go down with it, so a line behind another plane is tinted + // by it exactly like the plane it lies on. + std::vector> strokes; + for (const auto& [a, b] : piece.lines) strokes.push_back({ a, b }); + GLModel::Geometry outline; + append_ribbons(outline, -1, strokes, vd, Vec3d::Zero(), hw); // body -1: already world coordinates + if (!outline.is_empty()) { + GLModel om; om.init_from(std::move(outline)); + col.a(hot ? kLineAlphaHot : kLineAlpha); + om.set_color(col); + om.render(); } } - m_plane = saved_plane; // draw_text renders each label immediately (draw_strokes self-renders) + + // Label near the top-left corner, drawn in the plane (draw_text lifts through m_plane). Labels + // are ImGui chips, on top of the planes whatever the order here. + const SketchPlane saved_plane = m_plane; + const double th = H * 0.10; + for (size_t i = 0; i < m_dbp_planes.size() && i < m_dbp_labels.size(); ++i) { + if (m_dbp_labels[i].empty()) continue; + m_plane = m_dbp_planes[i]; + draw_text(m_line_model, m_dbp_labels[i], Vec2d(-H + th * 2.0, H - th * 1.6), th, + hue(i < m_dbp_base.size() ? m_dbp_base[i] : -1)); + } + m_plane = saved_plane; glsafe(::glDisable(GL_BLEND)); } @@ -5474,7 +5636,7 @@ int DesignSketchTool::hit_test_base_pick(GLCanvas3D& canvas, const wxMouseEvent& // THE LABEL WINS, and it has to. Each plane's name is a screen-space chip centred on its // own in-plane anchor, and it is the one part of a base plane a user aims at deliberately — - // the quads are near-transparent and overlap everywhere. Ray-casting the quads alone made + // the quads are translucent and overlap everywhere. Ray-casting the quads alone made // the labels pure decoration: on a fresh document at 1920x1060, clicking "XY" reported // "XZ plane selected", because the XZ quad happens to sit in front at that pixel. Nothing // about the click was ambiguous to the user; they clicked the word XY. diff --git a/src/slic3r/GUI/CAD/DesignSketchTool.hpp b/src/slic3r/GUI/CAD/DesignSketchTool.hpp index 91adb207a0..8d575241e9 100644 --- a/src/slic3r/GUI/CAD/DesignSketchTool.hpp +++ b/src/slic3r/GUI/CAD/DesignSketchTool.hpp @@ -54,6 +54,22 @@ inline ColorRGBA design_idle_face_color() { return ColorRGBA(0.72f, 0.76f, 0.80f, 0.14f); } + +// A convex piece of the square drawn on planes[plane], and the segments to outline with it: its share +// of the square's border and of the lines where it crosses the other squares. +struct PlanePiece +{ + int plane; + std::vector corners; + std::vector> lines; +}; +// The squares of half-extent `half` on `planes`, cut where they cross one another and ordered back +// to front for an eye at `eye` (perspective) or looking along `forward` (orthographic). Translucent +// planes that cross cannot be drawn in any per-plane order: each is partly in front of and partly +// behind the others. Drawn piece by piece in this order, each one tints only what is behind it. +std::vector planes_back_to_front(const std::vector& planes, double half, + const Vec3d& eye, const Vec3d& forward, bool perspective); + class DesignSketchTool { public: enum class Mode { Select, Dimension, Polyline, Line, CornerRect, CenterRect, ObliqueRect, diff --git a/tests/slic3rutils/test_design_sketch_tool.cpp b/tests/slic3rutils/test_design_sketch_tool.cpp index 7b99cde012..3446d813d2 100644 --- a/tests/slic3rutils/test_design_sketch_tool.cpp +++ b/tests/slic3rutils/test_design_sketch_tool.cpp @@ -13,8 +13,19 @@ #endif #include +#include +#include #include +#include +#include +#include +#include +#include +#include +#include +#include + #include "libslic3r/BoundingBox.hpp" #include "libslic3r/CAD/SketchEngine.hpp" #include "libslic3r/Point.hpp" @@ -24,6 +35,7 @@ using namespace Slic3r; using namespace Slic3r::GUI; using Catch::Matchers::WithinAbs; +using Catch::Matchers::WithinRel; namespace { @@ -43,6 +55,112 @@ void show_two_sketches(DesignSketchTool& tool) { { circle({ -60., 0. }, 5.) }, SketchPlane::XY(), 2 } }); } +// The Design tab's base planes, all through one origin, as DesignPanel shows them. +std::vector base_planes() { return { SketchPlane::XY(), SketchPlane::XZ(), SketchPlane::YZ() }; } + +struct View +{ + Vec3d eye; + Vec3d forward; + bool perspective; +}; + +// Eyes in four octants, off every plane, plus two orthographic directions. +const View kViews[] = { + { Vec3d(300., -400., 250.), Vec3d(-300., 400., -250.).normalized(), true }, + { Vec3d(-350., -200., 300.), Vec3d(350., 200., -300.).normalized(), true }, + { Vec3d(250., 300., -200.), Vec3d(-250., -300., 200.).normalized(), true }, + { Vec3d(-300., 350., -250.), Vec3d(300., -350., 250.).normalized(), true }, + { Vec3d::Zero(), Vec3d(-0.5, 0.7, -0.5).normalized(), false }, + { Vec3d::Zero(), Vec3d(0.3, 0.4, 0.85).normalized(), false }, +}; + +double area(const PlanePiece& piece, const SketchPlane& plane) +{ + Vec3d sum = Vec3d::Zero(); + for (size_t i = 0; i < piece.corners.size(); ++i) + sum += piece.corners[i].cross(piece.corners[(i + 1) % piece.corners.size()]); + return 0.5 * std::abs(sum.dot(plane.x_axis.cross(plane.y_axis))); +} + +// How far along the ray (from, unit dir) it crosses `piece`, or nothing if it misses. +std::optional hit_distance(const PlanePiece& piece, const SketchPlane& plane, const Vec3d& from, const Vec3d& dir) +{ + const Vec3d n = plane.x_axis.cross(plane.y_axis); + const double dn = n.dot(dir); + if (std::abs(dn) < 1e-9) + return std::nullopt; + const double t = n.dot(piece.corners.front() - from) / dn; + if (t <= 0.) + return std::nullopt; + const Vec3d x = from + dir * t; + double lo = 0., hi = 0.; + for (size_t i = 0; i < piece.corners.size(); ++i) { + const Vec3d& a = piece.corners[i]; + const Vec3d& b = piece.corners[(i + 1) % piece.corners.size()]; + const double s = (b - a).cross(x - a).dot(n); + lo = std::min(lo, s); + hi = std::max(hi, s); + } + if (lo < 0. && hi > 0.) + return std::nullopt; // outside one of the edges + return t; +} + +// Whether x lies on one of the segments the pieces are outlined with. +bool outlined(const std::vector& pieces, const Vec3d& x) +{ + for (const PlanePiece& piece : pieces) + for (const auto& [a, b] : piece.lines) { + const Vec3d ab = b - a; + const double t = std::clamp((x - a).dot(ab) / ab.squaredNorm(), 0., 1.); + if ((a + ab * t - x).norm() < 1e-6) + return true; + } + return false; +} + +std::vector pieces_of(const std::vector& planes) +{ + const View& view = kViews[0]; + return planes_back_to_front(planes, 75., view.eye, view.forward, view.perspective); +} + +struct SightLines +{ + int overlapping = 0; // sight lines through two or more pieces, where draw order matters + int out_of_order = 0; // ...of which meet a nearer piece before a farther one +}; + +// Translucent pieces blend correctly only if, along every line of sight, each piece is drawn after +// every piece behind it. +SightLines sight_lines(const std::vector& planes, double half, const View& view) +{ + const std::vector pieces = planes_back_to_front(planes, half, view.eye, view.forward, view.perspective); + std::mt19937 rng(7); + std::uniform_real_distribution coord(-0.95 * half, 0.95 * half); + SightLines seen; + for (int r = 0; r < 500; ++r) { + const Vec3d target(coord(rng), coord(rng), coord(rng)); + const Vec3d from = view.perspective ? view.eye : Vec3d(target - view.forward * (10. * half)); + const Vec3d dir = (target - from).normalized(); + double last = std::numeric_limits::max(); + int hits = 0; + bool ok = true; + for (const PlanePiece& piece : pieces) + if (const std::optional t = hit_distance(piece, planes[piece.plane], from, dir)) { + ok = ok && *t <= last + 1e-6 * half; + last = *t; + ++hits; + } + if (hits >= 2) { + ++seen.overlapping; + seen.out_of_order += ok ? 0 : 1; + } + } + return seen; +} + } // namespace TEST_CASE("Fit frames the picked sketch region, not the other sketches", "[DesignSketchTool]") @@ -108,3 +226,82 @@ TEST_CASE("Fit frames nothing when the Design tab shows nothing", "[DesignSketch GLVolumeCollection no_bodies; CHECK_FALSE(tool.fit_box(no_bodies).defined); } + +TEST_CASE("Crossing base planes are drawn back to front from any viewpoint", "[DesignSketchTool]") +{ + const View& view = kViews[GENERATE(range(0, int(std::size(kViews))))]; + const SightLines seen = sight_lines(base_planes(), 75., view); + CHECK(seen.overlapping >= 200); // of the 500: most lines of sight into the planes cross two + CHECK(seen.out_of_order == 0); +} + +TEST_CASE("Datum planes are drawn back to front among the base planes", "[DesignSketchTool]") +{ + // A datum parallel to XY 30 mm up, and one tilted 30 degrees about X through (0, 0, 10). + std::vector planes = base_planes(); + SketchPlane raised = SketchPlane::XY(); + raised.origin = Vec3d(0., 0., 30.); + SketchPlane tilted; + tilted.origin = Vec3d(0., 0., 10.); + tilted.y_axis = Vec3d(0., std::cos(M_PI / 6.), std::sin(M_PI / 6.)); + tilted.normal = tilted.x_axis.cross(tilted.y_axis); + planes.push_back(raised); + planes.push_back(tilted); + + const View& view = kViews[GENERATE(range(0, int(std::size(kViews))))]; + const SightLines seen = sight_lines(planes, 75., view); + CHECK(seen.overlapping >= 200); + CHECK(seen.out_of_order == 0); +} + +TEST_CASE("Base planes are outlined along every line where they cross", "[DesignSketchTool]") +{ + // XY, XZ and YZ cross along the three axes, all through the middle of each 75 mm half-square. + const std::vector pieces = pieces_of(base_planes()); + int missed = 0; + for (double t = -70.; t <= 70.; t += 10.) + for (const Vec3d& axis : { Vec3d(1., 0., 0.), Vec3d(0., 1., 0.), Vec3d(0., 0., 1.) }) + missed += outlined(pieces, axis * t) ? 0 : 1; + CHECK(missed == 0); +} + +TEST_CASE("A datum clear of the base planes gets no lines across it", "[DesignSketchTool]") +{ + // Parallel to YZ at x = 100, past the 75 mm half-squares of XY and XZ: the infinite XY and XZ + // planes still cut it, along z = 0 and y = 0, but the squares never meet. + std::vector planes = base_planes(); + SketchPlane beyond = SketchPlane::YZ(); + beyond.origin = Vec3d(100., 0., 0.); + planes.push_back(beyond); + + const std::vector pieces = pieces_of(planes); + CHECK_FALSE(outlined(pieces, Vec3d(100., 30., 0.))); + CHECK_FALSE(outlined(pieces, Vec3d(100., 0., 30.))); +} + +TEST_CASE("A line where two squares cross stops where either square ends", "[DesignSketchTool]") +{ + // Parallel to XY 30 mm up and moved 60 mm along X, so it spans x = -15..135. It meets XZ along + // y = 0, z = 30, but XZ's square only reaches x = 75. + std::vector planes = base_planes(); + SketchPlane raised = SketchPlane::XY(); + raised.origin = Vec3d(60., 0., 30.); + planes.push_back(raised); + + const std::vector pieces = pieces_of(planes); + CHECK(outlined(pieces, Vec3d(0., 0., 30.))); + CHECK(outlined(pieces, Vec3d(70., 0., 30.))); + CHECK_FALSE(outlined(pieces, Vec3d(100., 0., 30.))); +} + +TEST_CASE("Cutting the base planes along each other keeps every plane whole", "[DesignSketchTool]") +{ + const std::vector planes = base_planes(); + const double half = 75.; + const View& view = kViews[0]; + std::vector covered(planes.size(), 0.); + for (const PlanePiece& piece : planes_back_to_front(planes, half, view.eye, view.forward, view.perspective)) + covered[piece.plane] += area(piece, planes[piece.plane]); + for (double a : covered) + CHECK_THAT(a, WithinRel(4. * half * half, 1e-9)); +}