Skip to content

Commit

Permalink
Fix #10656 - Rewrite the CLIPLINE algorithm
Browse files Browse the repository at this point in the history 
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
```
  • Loading branch information
@jmarrec
jmarrec committed Aug 12, 2024
1 parent 55edc57 commit 251ea56
Show file tree
Hide file tree
Showing 3 changed files with 84 additions and 161 deletions.
233 changes: 81 additions & 152 deletions src/EnergyPlus/SolarShading.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -46,6 +46,7 @@
// POSSIBILITY OF SUCH DAMAGE.

// C++ Headers
#include <algorithm>
#include <cassert>
#include <cmath>
#include <memory>
Expand Down Expand Up @@ -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)
Expand Down Expand Up @@ -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;
Expand Down
2 changes: 1 addition & 1 deletion src/EnergyPlus/SolarShading.hh
View file
Original file line number Diff line number Diff line change
Expand Up @@ -136,7 +136,7 @@ namespace SolarShading {

void CLIP(EnergyPlusData &state, int const NVT, Array1D<Real64> &XVT, Array1D<Real64> &YVT, Array1D<Real64> &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
Expand Down
10 changes: 2 additions & 8 deletions tst/EnergyPlus/unit/SolarShading.unit.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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);
Expand Down Expand Up @@ -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;
Expand Down

4 comments on commit 251ea56

@nrel-bot-2c
Copy link

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

10656_CLIPLINE (jmarrec) - x86_64-Linux-Ubuntu-22.04-gcc-11.4: OK (3696 of 3697 tests passed, 0 test warnings)

Failures:\n

EnergyPlusFixture Test Summary

  • Passed: 1578
  • Failed: 1

Build Badge Test Badge

@nrel-bot-3
Copy link

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

10656_CLIPLINE (jmarrec) - x86_64-MacOS-10.18-clang-15.0.0: OK (3655 of 3656 tests passed, 0 test warnings)

Failures:\n

EnergyPlusFixture Test Summary

  • Passed: 1578
  • Failed: 1

Build Badge Test Badge

@nrel-bot-2b
Copy link

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

10656_CLIPLINE (jmarrec) - x86_64-Linux-Ubuntu-22.04-gcc-11.4-UnitTestsCoverage-Debug: OK (2071 of 2072 tests passed, 0 test warnings)

Failures:\n

EnergyPlusFixture Test Summary

  • Passed: 1578
  • Failed: 1

Build Badge Test Badge Coverage Badge

@nrel-bot-2
Copy link

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

10656_CLIPLINE (jmarrec) - x86_64-Linux-Ubuntu-22.04-gcc-11.4-IntegrationCoverage-Debug: OK (795 of 795 tests passed, 0 test warnings)

Build Badge Test Badge Coverage Badge

Please sign in to comment.