From 251ea560e504b2bf3f1b79dbbb9f73bae3aa9e7c Mon Sep 17 00:00:00 2001 From: Julien Marrec Date: Mon, 12 Aug 2024 12:27:56 +0200 Subject: [PATCH] Fix #10656 - Rewrite the CLIPLINE algorithm Use the Liang-Barsky algorithm * Loosely adapted from https://en.wikipedia.org/wiki/Liang%E2%80%93Barsky_algorithm * Rewrote it until performance is optimal I used microbenchmarking: https://github.com/jmarrec/CppBenchmarks/blob/94f0e1594ff17f9833a0547f0725d68c561311a1/bench_clipline.cpp BM_CLIPLINE is the current E+ algorithm. BM_liang_barsky_clipper3 is my final version that I used here, and which is 15% faster than the original while being actually correct. in the edge cases ``` ------------------------------------------------------------------- Benchmark Time CPU Iterations ------------------------------------------------------------------- BM_CLIPLINE 567 ns 566 ns 1231826 BM_Clip2 703 ns 703 ns 991672 BM_liang_barsky_clipper 588 ns 588 ns 1199750 BM_liang_barsky_clipper2 490 ns 490 ns 1430600 BM_liang_barsky_clipper3 482 ns 482 ns 145236 ``` --- src/EnergyPlus/SolarShading.cc | 233 ++++++++--------------- src/EnergyPlus/SolarShading.hh | 2 +- tst/EnergyPlus/unit/SolarShading.unit.cc | 10 +- 3 files changed, 84 insertions(+), 161 deletions(-) diff --git a/src/EnergyPlus/SolarShading.cc b/src/EnergyPlus/SolarShading.cc index 4f01b2b0583..df6c21bf79b 100644 --- a/src/EnergyPlus/SolarShading.cc +++ b/src/EnergyPlus/SolarShading.cc @@ -46,6 +46,7 @@ // POSSIBILITY OF SUCH DAMAGE. // C++ Headers +#include #include #include #include @@ -3868,162 +3869,99 @@ inline bool d_eq(Real64 a, Real64 b) return std::abs(a - b) < 2.0; } -void CLIPLINE(Real64 &x1, Real64 &x2, Real64 &y1, Real64 &y2, Real64 maxX, Real64 minX, Real64 maxY, Real64 minY, bool &visible, bool &rev) +void CLIPLINE(Real64 &x0, Real64 &x1, Real64 &y0, Real64 &y1, Real64 maxX, Real64 minX, Real64 maxY, Real64 minY, bool &visible) { + // Line segment clipping // Reference: - // Slater, M., Barsky, B.A. + // Liang, Y.D., Barsky, B.A., Slater, M. // 2D line and polygon clipping based on space subdivision. // The Visual Computer 10, 407–422 (1994). - Real64 dx, dy, e, xinc, yinc, tempVar; - bool needX = true, needY = true; - int c1, c2; - - if (x1 > x2) { // reverse for efficiency - tempVar = x1; - x1 = x2; - x2 = tempVar; - tempVar = y1; - y1 = y2; - y2 = tempVar; + + // Tweaked via microbenchmarking to improve efficiency + + bool rev = false; + if (x0 > x1) { // reverse for efficiency + std::swap(x0, x1); + std::swap(y0, y1); rev = true; } - if (x1 > maxX || x2 < minX) return; // x is positive - if (x1 < minX) { - if (y1 < minY) { - if (y2 < minY) return; - c1 = 0; - dx = x2 - x1; - dy = y2 - y1; - e = dy * (minX - x1) + dx * (y1 - minY); - } else if (y1 > maxY) { - if (y2 > maxY) return; - c1 = 6; - dx = x2 - x1; - dy = y2 - y1; - e = dy * (minX - x1) + dx * (y1 - maxY); - } else { - c1 = 3; - dx = x2 - x1; - dy = y2 - y1; - if (dy > 0) { - e = dy * (minX - x1) + dx * (y1 - maxY); - } else { - e = dy * (minX - x1) + dx * (y1 - minY); - } + + if (x0 > maxX || x1 < minX) { + // Both points are outside the clip window, so they can't cross it + return; + } + + // defining variables + Real64 const dx = x1 - x0; // >= 0 + Real64 const dy = y1 - y0; + + Real64 const q1 = x0 - minX; + Real64 const q2 = maxX - x0; + Real64 const q3 = y0 - minY; + Real64 const q4 = maxY - y0; + + Real64 u1 = 0; + Real64 u2 = 1; + + if ((dx == 0 && (q1 < 0 || q2 < 0)) || (dy == 0 && (q3 < 0 || q4 < 0))) { + // Line is parallel to clipping window + return; + } + if (dx != 0) { + Real64 const r1 = q1 / -dx; + if (r1 > u1) { + u1 = r1; } - } else { - if (y1 < minY) { - if (y2 < minY) return; - c1 = 1; - dx = x2 - x1; - dy = y2 - y1; - e = dy * (maxX - x1) + dx * (y1 - minY); - } else if (y1 > maxY) { - if (y2 > maxY) return; - c1 = 7; - dx = x2 - x1; - dy = y2 - y1; - e = dy * (maxX - x1) + dx * (y1 - maxY); - } else { - visible = true; - if (x2 <= maxX && (y2 >= minY && y2 <= maxY)) return; - c1 = 4; - dx = x2 - x1; - dy = y2 - y1; - if (dy > 0) { - e = dy * (maxX - x1) + dx * (y1 - maxY); - } else { - e = dy * (maxX - x1) + dx * (y1 - minY); - } + Real64 const r2 = q2 / dx; + if (r2 < u2) { + u2 = r2; } } - c2 = c1; - if (dy > 0) { - while (true) { - if (e < 0.0) { - if (c2 == 1) - return; - else if (c2 == 3) { - visible = true; - x1 = minX; - y1 = maxY + e / dx; - if (x2 <= maxX && y2 <= maxY) return; - } else if (c2 == 4) { - x2 = maxX; - y2 = maxY + e / dx; - return; - } - if (needX) { - xinc = dy * (maxX - minX); - needX = false; - } - e += xinc; - c2 += 1; - } else { - if (c2 == 3) - return; - else if (c2 == 1) { - visible = true; - x1 = maxX - e / dy; - y1 = minY; - if (x2 <= maxX && y2 <= maxY) return; - } else if (c2 == 4) { - x2 = maxX - e / dy; - y2 = maxY; - return; - } - if (needY) { - yinc = dx * (maxY - minY); - needY = false; - } - e -= yinc; - c2 += 3; + if (dy != 0) { + Real64 const r3 = q3 / -dy; + Real64 const r4 = q4 / dy; + if (dy > 0) { + if (r3 > u1) { + u1 = r3; } - } - } else { - while (true) { - if (e >= 0.0) { - if (c2 == 7) - return; - else if (c2 == 3) { - visible = true; - x1 = minX; - y1 = minY + e / dx; - if (x2 <= maxX && y2 >= minY) return; - } else if (c2 == 4) { - x2 = maxX; - y2 = minY + e / dx; - return; - } - if (needX) { - xinc = dy * (maxX - minX); - needX = false; - } - e += xinc; - c2 += 1; - } else { - if (c2 == 3) - return; - else if (c2 == 7) { - visible = true; - x1 = maxX - e / dy; - y1 = maxY; - if (x2 <= maxX && y2 >= minY) return; - } else if (c2 == 4) { - x2 = maxX - e / dy; - y2 = minY; - return; - } - if (needY) { - yinc = dx * (maxY - minY); - needY = false; - } - e += yinc; - c2 -= 3; + if (r4 < u2) { + u2 = r4; + } + } else { + if (r4 > u1) { + u1 = r4; + } + if (r3 < u2) { + u2 = r3; } } } + + if (u1 > u2) { // reject + // Line is outside the clipping window + return; + } + + visible = true; + + Real64 const xn0 = x0 + dx * u1; + Real64 const yn0 = y0 + dy * u1; + + Real64 const xn1 = x0 + dx * u2; + Real64 const yn1 = y0 + dy * u2; + + if (rev) { + x0 = xn1; + y0 = yn1; + x1 = xn0; + y1 = yn0; + } else { + x0 = xn0; + y0 = yn0; + x1 = xn1; + y1 = yn1; + } } void CLIPRECT(EnergyPlusData &state, int const NS2, int const NV1, int &NV3) @@ -4063,20 +4001,11 @@ void CLIPRECT(EnergyPlusData &state, int const NS2, int const NV1, int &NV3) Real64 x1 = x_1, x2 = x_2, y1 = y_1, y2 = y_2; bool visible = false; - bool rev = false; - CLIPLINE(x_1, x_2, y_1, y_2, maxX, minX, maxY, minY, visible, rev); + CLIPLINE(x_1, x_2, y_1, y_2, maxX, minX, maxY, minY, visible); if (visible) { if ((x_1 != x1 || y_1 != y1) || (x_2 != x2 || y_2 != y2)) { INTFLAG = true; } - if (rev) { // undo reverse - Real64 tempVar = x_1; - x_1 = x_2; - x_2 = tempVar; - tempVar = y_1; - y_1 = y_2; - y_2 = tempVar; - } // if line on edge, or inside, add both points if (arrc == 0 || ((neq(arrx[arrc - 1], x_1) || neq(arry[arrc - 1], y_1)) && (neq(arrx[0], x_1) || neq(arry[0], y_1)))) { arrx[arrc] = x_1; diff --git a/src/EnergyPlus/SolarShading.hh b/src/EnergyPlus/SolarShading.hh index f91a571901b..31878be083d 100644 --- a/src/EnergyPlus/SolarShading.hh +++ b/src/EnergyPlus/SolarShading.hh @@ -136,7 +136,7 @@ namespace SolarShading { void CLIP(EnergyPlusData &state, int const NVT, Array1D &XVT, Array1D &YVT, Array1D &ZVT); - void CLIPLINE(Real64 &x0, Real64 &x1, Real64 &y0, Real64 &y1, Real64 maxX, Real64 minX, Real64 maxY, Real64 minY, bool &visible, bool &rev); + void CLIPLINE(Real64 &x0, Real64 &x1, Real64 &y0, Real64 &y1, Real64 maxX, Real64 minX, Real64 maxY, Real64 minY, bool &visible); void CTRANS(EnergyPlusData &state, int const NS, // Surface number whose vertex coordinates are being transformed diff --git a/tst/EnergyPlus/unit/SolarShading.unit.cc b/tst/EnergyPlus/unit/SolarShading.unit.cc index 302a684be5b..68d481be86b 100644 --- a/tst/EnergyPlus/unit/SolarShading.unit.cc +++ b/tst/EnergyPlus/unit/SolarShading.unit.cc @@ -6217,9 +6217,8 @@ TEST_F(EnergyPlusFixture, CLIPLINE_Throw) Real64 y0 = 4.5; Real64 y1 = 1.0; bool visible = false; - bool rev = false; - EXPECT_NO_THROW(CLIPLINE(x0, x1, y0, y1, maxX, minX, maxY, minY, visible, rev)); + EXPECT_NO_THROW(CLIPLINE(x0, x1, y0, y1, maxX, minX, maxY, minY, visible)); EXPECT_DOUBLE_EQ(maxX, x0); EXPECT_DOUBLE_EQ(4.5, y0); @@ -6271,12 +6270,7 @@ TEST_F(EnergyPlusFixture, CLIPLINE_Full) std::string const msg = fmt::format("From ({}, {}) to ({}, {})", t.line_ori.p0.x, t.line_ori.p0.y, t.line_ori.p1.x, t.line_ori.p1.y); bool visible = false; - bool rev = false; - CLIPLINE(x0, x1, y0, y1, maxX, minX, maxY, minY, visible, rev); - if (rev) { - std::swap(x0, x1); - std::swap(y0, y1); - } + CLIPLINE(x0, x1, y0, y1, maxX, minX, maxY, minY, visible); if (t.visible) { EXPECT_TRUE(visible) << msg; EXPECT_DOUBLE_EQ(t.line_new.p0.x, x0) << msg;