Skip to content

Commit 7cd865f

Browse files
abellgithubabellgithub
committed
Checkpoint.
1 parent 2b8cbaa commit 7cd865f

File tree

3 files changed

+185
-55
lines changed

3 files changed

+185
-55
lines changed

‎alg/viewshed/viewshed_executor.cpp‎

Lines changed: 53 additions & 52 deletions
Original file line numberDiff line numberDiff line change
@@ -16,6 +16,8 @@
1616
#include <cassert>
1717
#include <cmath>
1818
#include <limits>
19+
//ABELL
20+
#include <iostream>
1921

2022
#include "viewshed_executor.h"
2123
#include "progress.h"
@@ -263,11 +265,11 @@ LineLimits ViewshedExecutor::adjustHeight(int nYOffset,
263265
double dfR2 = dfX * dfX + dfY * dfY;
264266
if (dfR2 > m_dfMaxDistance2)
265267
{
268+
if (dfR2 < m_dfMinDistance2)
269+
ll.leftMin++;
266270
ll.left = nXOffset + m_nX + 1;
267271
break;
268272
}
269-
if (dfR2 < m_dfMinDistance2)
270-
ll.leftMin++;
271273

272274
CheckNoData(*pdfHeight);
273275
*pdfHeight -= m_dfHeightAdjFactor * dfR2 + m_dfZObserver;
@@ -283,11 +285,11 @@ LineLimits ViewshedExecutor::adjustHeight(int nYOffset,
283285
double dfR2 = dfX * dfX + dfY * dfY;
284286
if (dfR2 > m_dfMaxDistance2)
285287
{
288+
if (dfR2 < m_dfMinDistance2)
289+
ll.rightMin++;
286290
ll.right = nXOffset + m_nX;
287291
break;
288292
}
289-
if (dfR2 < m_dfMinDistance2)
290-
ll.rightMin++;
291293

292294
CheckNoData(*pdfHeight);
293295
*pdfHeight -= m_dfHeightAdjFactor * dfR2 + m_dfZObserver;
@@ -338,22 +340,18 @@ bool ViewshedExecutor::processFirstLine(std::vector<double> &vLastLineVal)
338340
if (oOpts.outputMode == OutputMode::DEM)
339341
vResult = vThisLineVal;
340342

341-
// iLeft and iRight are the processing limits for the line.
342343
LineLimits ll = adjustHeight(nYOffset, vThisLineVal);
343344

344345
if (!oCurExtent.containsY(m_nY))
345346
processFirstLineTopOrBottom(ll, vResult, vThisLineVal);
346347
else
347348
{
348349
CPLJobQueuePtr pQueue = m_pool.CreateJobQueue();
349-
pQueue->SubmitJob(
350-
[&, left = ll.left]() {
351-
processFirstLineLeft(m_nX - 1, left - 1, vResult, vThisLineVal);
352-
});
350+
pQueue->SubmitJob([&]()
351+
{ processFirstLineLeft(ll, vResult, vThisLineVal); });
353352

354353
pQueue->SubmitJob(
355-
[&, right = ll.right]()
356-
{ processFirstLineRight(m_nX + 1, right, vResult, vThisLineVal); });
354+
[&]() { processFirstLineRight(ll, vResult, vThisLineVal); });
357355
pQueue->WaitCompletion();
358356
}
359357

@@ -392,14 +390,16 @@ void ViewshedExecutor::processFirstLineTopOrBottom(
392390

393391
/// Process the part of the first line to the left of the observer.
394392
///
395-
/// @param iStart X coordinate of the first cell to the left of the observer to be procssed.
396-
/// @param iEnd X coordinate one past the last cell to be processed.
393+
/// @param ll Line limits for masking.
397394
/// @param vResult Vector in which to store the visibility/height results.
398395
/// @param vThisLineVal Height of each cell in the line being processed.
399-
void ViewshedExecutor::processFirstLineLeft(int iStart, int iEnd,
396+
void ViewshedExecutor::processFirstLineLeft(const LineLimits &ll,
400397
std::vector<double> &vResult,
401398
std::vector<double> &vThisLineVal)
402399
{
400+
int iEnd = ll.left - 1;
401+
int iStart = m_nX - 1; // One left of the observer.
402+
403403
// If end is to the right of start, everything is taken care of by right processing.
404404
if (iEnd >= iStart)
405405
return;
@@ -426,20 +426,32 @@ void ViewshedExecutor::processFirstLineLeft(int iStart, int iEnd,
426426
double dfZ = CalcHeightLine(nXOffset, *(pThis + 1));
427427
setOutput(vResult[iPixel], *pThis, dfZ);
428428
}
429-
// For cells outside of the [start, end) range, set the outOfRange value.
430-
std::fill(vResult.begin(), vResult.begin() + iEnd + 1, oOpts.outOfRangeVal);
429+
430+
maskLineLeft(vResult, ll);
431+
}
432+
433+
void ViewshedExecutor::maskLineLeft(std::vector<double> &vResult,
434+
const LineLimits &ll)
435+
{
436+
// Mask cells from the left edge to the left limit.
437+
std::fill(vResult.begin(), vResult.begin() + ll.left, oOpts.outOfRangeVal);
438+
// Mask cells from the left min to the observer.
439+
std::fill(vResult.begin() + ll.leftMin - 1, vResult.begin() + m_nX,
440+
oOpts.outOfRangeVal);
431441
}
432442

433443
/// Process the part of the first line to the right of the observer.
434444
///
435-
/// @param iStart X coordinate of the first cell to the right of the observer to be processed.
436-
/// @param iEnd X coordinate one past the last cell to be processed.
445+
/// @param ll Line limits
437446
/// @param vResult Vector in which to store the visibility/height results.
438447
/// @param vThisLineVal Height of each cell in the line being processed.
439-
void ViewshedExecutor::processFirstLineRight(int iStart, int iEnd,
448+
void ViewshedExecutor::processFirstLineRight(const LineLimits &ll,
440449
std::vector<double> &vResult,
441450
std::vector<double> &vThisLineVal)
442451
{
452+
int iStart = m_nX + 1;
453+
int iEnd = ll.right;
454+
443455
// If start is to the right of end, everything is taken care of by left processing.
444456
if (iStart >= iEnd)
445457
return;
@@ -466,23 +478,38 @@ void ViewshedExecutor::processFirstLineRight(int iStart, int iEnd,
466478
double dfZ = CalcHeightLine(nXOffset, *(pThis - 1));
467479
setOutput(vResult[iPixel], *pThis, dfZ);
468480
}
481+
482+
maskLineRight(vResult, ll);
483+
469484
// For cells outside of the [start, end) range, set the outOfRange value.
470-
std::fill(vResult.begin() + iEnd, vResult.end(), oOpts.outOfRangeVal);
485+
//std::fill(vResult.begin() + iEnd, vResult.end(), oOpts.outOfRangeVal);
486+
}
487+
488+
void ViewshedExecutor::maskLineRight(std::vector<double> &vResult,
489+
const LineLimits &ll)
490+
{
491+
// Mask cells from the right limit to the right edge.
492+
std::fill(vResult.begin() + ll.right, vResult.end(), oOpts.outOfRangeVal);
493+
// Mask cells from the observer to right min.
494+
std::fill(vResult.begin() + m_nX + 1, vResult.begin() + ll.rightMin,
495+
oOpts.outOfRangeVal);
471496
}
472497

473498
/// Process a line to the left of the observer.
474499
///
475500
/// @param nYOffset Offset of the line being processed from the observer
476-
/// @param iStart X coordinate of the first cell to the left of the observer to be processed.
477-
/// @param iEnd X coordinate one past the last cell to be processed.
501+
/// @param ll Line limits
478502
/// @param vResult Vector in which to store the visibility/height results.
479503
/// @param vThisLineVal Height of each cell in the line being processed.
480504
/// @param vLastLineVal Observable height of each cell in the previous line processed.
481-
void ViewshedExecutor::processLineLeft(int nYOffset, int iStart, int iEnd,
505+
void ViewshedExecutor::processLineLeft(int nYOffset, LineLimits &ll,
482506
std::vector<double> &vResult,
483507
std::vector<double> &vThisLineVal,
484508
std::vector<double> &vLastLineVal)
485509
{
510+
int iStart = m_nX - 1;
511+
int iEnd = ll.left - 1;
512+
486513
// If start to the left of end, everything is taken care of by processing right.
487514
if (iStart <= iEnd)
488515
return;
@@ -524,8 +551,7 @@ void ViewshedExecutor::processLineLeft(int nYOffset, int iStart, int iEnd,
524551
setOutput(vResult[iPixel], *pThis, dfZ);
525552
}
526553

527-
// For cells outside of the [start, end) range, set the outOfRange value.
528-
std::fill(vResult.begin(), vResult.begin() + iEnd + 1, oOpts.outOfRangeVal);
554+
maskLineLeft(vResult, ll);
529555
}
530556

531557
/// Process a line to the right of the observer.
@@ -628,10 +654,8 @@ bool ViewshedExecutor::processLine(int nLine, std::vector<double> &vLastLineVal)
628654
// process left half then right half of line
629655
CPLJobQueuePtr pQueue = m_pool.CreateJobQueue();
630656
pQueue->SubmitJob(
631-
[&, left = ll.left]()
632-
{
633-
processLineLeft(nYOffset, m_nX - 1, left - 1, vResult, vThisLineVal,
634-
vLastLineVal);
657+
[&]() {
658+
processLineLeft(nYOffset, ll, vResult, vThisLineVal, vLastLineVal);
635659
});
636660
pQueue->SubmitJob(
637661
[&, right = ll.right]()
@@ -650,29 +674,6 @@ bool ViewshedExecutor::processLine(int nLine, std::vector<double> &vLastLineVal)
650674
return oProgress.lineComplete();
651675
}
652676

653-
/**
654-
void ViewshedExecutor::maskLine(int nLine, std::vector<double> vLine)
655-
{
656-
// If there is no min distance or angle masking or the angle is in the lower half and
657-
// we are processing a line in the upper half (or vice versa), return.
658-
if (minDistance == 0)
659-
{
660-
if (startAngle == endAngle)
661-
return;
662-
if (m_nX - nLine >= 0) // Upper
663-
if (!rayBetween(0, M_PI, oOpts.startAngle) &&
664-
!rayBetween(0, M_PI, oOpts.endAngle))
665-
return;
666-
else // Lower
667-
if (!rayBetween(M_PI, 2 * M_PI, oOpts.startAngle) &&
668-
!rayBetween(M_PI, 2 * M_PI, oOpts.endAngle))
669-
return;
670-
}
671-
672-
for (x
673-
}
674-
**/
675-
676677
/// Run the viewshed computation
677678
/// @return Success as true or false.
678679
bool ViewshedExecutor::run()

‎alg/viewshed/viewshed_executor.h‎

Lines changed: 5 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -74,16 +74,16 @@ class ViewshedExecutor
7474
bool writeLine(int nLine, std::vector<double> &vResult);
7575
bool processLine(int nLine, std::vector<double> &vLastLineVal);
7676
bool processFirstLine(std::vector<double> &vLastLineVal);
77-
void processFirstLineLeft(int iStart, int iEnd,
77+
void processFirstLineLeft(const LineLimits &ll,
7878
std::vector<double> &vResult,
7979
std::vector<double> &vThisLineVal);
80-
void processFirstLineRight(int iStart, int iEnd,
80+
void processFirstLineRight(const LineLimits &ll,
8181
std::vector<double> &vResult,
8282
std::vector<double> &vThisLineVal);
8383
void processFirstLineTopOrBottom(const LineLimits &ll,
8484
std::vector<double> &vResult,
8585
std::vector<double> &vThisLineVal);
86-
void processLineLeft(int nYOffset, int iStart, int iEnd,
86+
void processLineLeft(int nYOffset, LineLimits &ll,
8787
std::vector<double> &vResult,
8888
std::vector<double> &vThisLineVal,
8989
std::vector<double> &vLastLineVal);
@@ -92,6 +92,8 @@ class ViewshedExecutor
9292
std::vector<double> &vThisLineVal,
9393
std::vector<double> &vLastLineVal);
9494
LineLimits adjustHeight(int iLine, std::vector<double> &thisLineVal);
95+
void maskLineLeft(std::vector<double> &vResult, const LineLimits &ll);
96+
void maskLineRight(std::vector<double> &vResult, const LineLimits &ll);
9597
};
9698

9799
} // namespace viewshed

‎autotest/cpp/test_viewshed_internal.cpp‎

Lines changed: 127 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,127 @@
1+
///////////////////////////////////////////////////////////////////////////////
2+
//
3+
// Project: Internal C++ Test Suite for GDAL/OGR
4+
// Purpose: Test viewshed algorithm
5+
// Author: Andrew Bell
6+
//
7+
///////////////////////////////////////////////////////////////////////////////
8+
/*
9+
* SPDX-License-Identifier: MIT
10+
****************************************************************************/
11+
12+
#include <math.h>
13+
#include <utility>
14+
#include <iostream>
15+
16+
//#include "gdal_unit_test.h"
17+
18+
#include "gtest_include.h"
19+
20+
#include "viewshed/viewshed.h"
21+
#include "viewshed/viewshed_types.h"
22+
#include "viewshed/util.h"
23+
24+
namespace gdal
25+
{
26+
namespace viewshed
27+
{
28+
29+
TEST(ViewshedInternal, angle)
30+
{
31+
EXPECT_DOUBLE_EQ(M_PI / 2, normalizeAngle(0));
32+
EXPECT_DOUBLE_EQ(M_PI / 4, normalizeAngle(45));
33+
EXPECT_DOUBLE_EQ(0.0, normalizeAngle(90));
34+
EXPECT_DOUBLE_EQ(7 * M_PI / 4, normalizeAngle(135));
35+
EXPECT_DOUBLE_EQ(3 * M_PI / 2, normalizeAngle(180));
36+
EXPECT_DOUBLE_EQ(M_PI, normalizeAngle(270));
37+
}
38+
39+
TEST(ViewshedInternal, between)
40+
{
41+
EXPECT_TRUE(rayBetween(M_PI, 0, M_PI / 2));
42+
EXPECT_FALSE(rayBetween(M_PI, 0, 3 * M_PI / 2));
43+
EXPECT_TRUE(rayBetween(0, 3 * M_PI / 2, 7 * M_PI / 4));
44+
EXPECT_TRUE(rayBetween(M_PI / 4, 7 * M_PI / 4, 0));
45+
EXPECT_FALSE(rayBetween(M_PI / 4, 7 * M_PI / 4, M_PI));
46+
}
47+
48+
TEST(ViewshedInternal, intersect)
49+
{
50+
// Top side
51+
EXPECT_TRUE(std::isnan(horizontalIntersect(0, 0, 0, -2)));
52+
EXPECT_TRUE(std::isnan(horizontalIntersect(M_PI, 0, 0, -2)));
53+
EXPECT_DOUBLE_EQ(horizontalIntersect(M_PI / 2, 0, 0, -2), 0);
54+
EXPECT_TRUE(std::isnan(horizontalIntersect(3 * M_PI / 2, 0, 0, -2)));
55+
EXPECT_DOUBLE_EQ(horizontalIntersect(M_PI / 4, 0, 0, -2), 2);
56+
EXPECT_DOUBLE_EQ(horizontalIntersect(3 * M_PI / 4, 0, 0, -2), -2);
57+
EXPECT_TRUE(std::isnan(horizontalIntersect(5 * M_PI / 4, 0, 0, -2)));
58+
EXPECT_DOUBLE_EQ(horizontalIntersect(M_PI / 6, 0, 0, -2), 2 * std::sqrt(3));
59+
60+
// Bottom side
61+
EXPECT_TRUE(std::isnan(horizontalIntersect(0, 0, 0, 2)));
62+
EXPECT_TRUE(std::isnan(horizontalIntersect(M_PI, 0, 0, 2)));
63+
EXPECT_TRUE(std::isnan(horizontalIntersect(M_PI / 2, 0, 0, 2)));
64+
EXPECT_DOUBLE_EQ(horizontalIntersect(3 * M_PI / 2, 0, 0, 2), 0);
65+
66+
EXPECT_DOUBLE_EQ(horizontalIntersect(5 * M_PI / 4, 0, 0, 2), -2);
67+
EXPECT_DOUBLE_EQ(horizontalIntersect(7 * M_PI / 4, 0, 0, 2), 2);
68+
EXPECT_TRUE(std::isnan(horizontalIntersect(3 * M_PI / 4, 0, 0, 2)));
69+
EXPECT_NEAR(horizontalIntersect(7 * M_PI / 6, 0, 0, 2), -2 * std::sqrt(3),
70+
1e-10);
71+
72+
// Right side
73+
EXPECT_DOUBLE_EQ(verticalIntersect(0, 0, 0, 2), 0);
74+
EXPECT_TRUE(std::isnan(verticalIntersect(M_PI, 0, 0, 2)));
75+
EXPECT_TRUE(std::isnan(verticalIntersect(M_PI / 2, 0, 0, 2)));
76+
EXPECT_TRUE(std::isnan(verticalIntersect(3 * M_PI / 2, 0, 0, 2)));
77+
EXPECT_TRUE(std::isnan(verticalIntersect(5 * M_PI / 4, 0, 0, 2)));
78+
EXPECT_DOUBLE_EQ(verticalIntersect(M_PI / 4, 0, 0, 2), -2);
79+
EXPECT_DOUBLE_EQ(verticalIntersect(7 * M_PI / 4, 0, 0, 2), 2);
80+
EXPECT_DOUBLE_EQ(verticalIntersect(M_PI / 6, 0, 0, 2), -2 / std::sqrt(3));
81+
82+
// Left side
83+
EXPECT_DOUBLE_EQ(verticalIntersect(M_PI, 0, 0, -2), 0);
84+
EXPECT_TRUE(std::isnan(verticalIntersect(0, 0, 0, -2)));
85+
EXPECT_TRUE(std::isnan(verticalIntersect(M_PI / 2, 0, 0, -2)));
86+
EXPECT_TRUE(std::isnan(verticalIntersect(3 * M_PI / 2, 0, 0, -2)));
87+
EXPECT_TRUE(std::isnan(verticalIntersect(3 * M_PI / 4, 0, 0, 2)));
88+
EXPECT_DOUBLE_EQ(verticalIntersect(3 * M_PI / 4, 0, 0, -2), -2);
89+
EXPECT_DOUBLE_EQ(verticalIntersect(5 * M_PI / 4, 0, 0, -2), 2);
90+
EXPECT_DOUBLE_EQ(verticalIntersect(5 * M_PI / 6, 0, 0, -2),
91+
-2 / std::sqrt(3));
92+
}
93+
94+
void testShrinkWindowForAngles(Window &w, int, int, double, double);
95+
96+
TEST(ViewshedInternal, shrinkbox)
97+
{
98+
auto testExtent = [](double start, double stop, Window expected)
99+
{
100+
Window extent{-3, 3, -2, 2};
101+
testShrinkWindowForAngles(extent, 0, 0, start, stop);
102+
EXPECT_EQ(extent, expected);
103+
};
104+
105+
testExtent(3 * M_PI / 4, M_PI / 4, {-2, 2, -2, 0});
106+
testExtent(M_PI / 4, 3 * M_PI / 4, {-3, 3, -2, 2});
107+
testExtent(0.321750554, 2 * M_PI - 0.321750554,
108+
{0, 3, -1, 1}); // <2, 1>, <2, -1>
109+
testExtent(2 * M_PI - 0.321750554, 0.321750554,
110+
{-3, 3, -2, 2}); // <2, -1>, <2, 1>
111+
testExtent(7 * M_PI / 4, 5 * M_PI / 4, {-2, 2, 0, 2});
112+
testExtent(5 * M_PI / 4, 7 * M_PI / 4, {-3, 3, -2, 2});
113+
testExtent(M_PI + 0.321750554, M_PI - 0.321750554,
114+
{-3, 0, -1, 1}); // <-2, -1>, <-2, 1>
115+
testExtent(M_PI - 0.321750554, M_PI + 0.321750554,
116+
{-3, 3, -2, 2}); // <-2, 1>, <-2, -1>
117+
testExtent(M_PI / 4, 0.321750554, {0, 3, -2, 0}); // <2, 2>, <2, 1>
118+
testExtent(0.321750554, M_PI / 4, {-3, 3, -2, 2}); // <2, 1>, <2, 2>
119+
testExtent(M_PI / 4, 7 * M_PI / 4, {0, 3, -2, 2});
120+
testExtent(M_PI / 4, M_PI + 0.321750554,
121+
{-3, 3, -2, 2}); // <2, 2>, <-2, -1>
122+
testExtent(M_PI + 0.321750554, M_PI / 4,
123+
{-3, 2, -2, 1}); // <-2, -1>, <2, 2>
124+
}
125+
126+
} // namespace viewshed
127+
} // namespace gdal

0 commit comments

Comments
 (0)