Skip to content

Commit a94a315

Browse files
abellgithubabellgithub
committed
Working viewshed.
1 parent 589feaf commit a94a315

File tree

8 files changed

+566
-60
lines changed

8 files changed

+566
-60
lines changed

‎alg/viewshed/viewshed.cpp‎

Lines changed: 29 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -373,6 +373,21 @@ bool Viewshed::calcExtents(int nX, int nY, const GDALGeoTransform &invGT)
373373
bool Viewshed::run(GDALRasterBandH band, GDALProgressFunc pfnProgress,
374374
void *pProgressArg)
375375
{
376+
return run(band, nullptr, pfnProgress, pProgressArg);
377+
}
378+
379+
/// Compute the viewshed of a raster band with .
380+
///
381+
/// @param band Pointer to the raster band to be processed.
382+
/// @param sdBand Pointer to the SD raster band to be processed.
383+
/// @param pfnProgress Pointer to the progress function. Can be null.
384+
/// @param pProgressArg Argument passed to the progress function
385+
/// @return True on success, false otherwise.
386+
bool Viewshed::run(GDALRasterBandH band, GDALRasterBandH sdBand,
387+
GDALProgressFunc pfnProgress, void *pProgressArg)
388+
{
389+
if (sdBand)
390+
pSdBand = static_cast<GDALRasterBand *>(sdBand);
376391
pSrcBand = static_cast<GDALRasterBand *>(band);
377392

378393
GDALGeoTransform fwdTransform, invTransform;
@@ -443,10 +458,20 @@ bool Viewshed::run(GDALRasterBandH band, GDALProgressFunc pfnProgress,
443458

444459
// Execute the viewshed algorithm.
445460
GDALRasterBand *pDstBand = poDstDS->GetRasterBand(1);
446-
ViewshedExecutor executor(*pSrcBand, *pDstBand, nX, nY, oOutExtent,
447-
oCurExtent, oOpts, oProgress,
448-
/* emitWarningIfNoData = */ true);
449-
executor.run();
461+
if (pSdBand)
462+
{
463+
ViewshedExecutor executor(*pSrcBand, *pSdBand, *pDstBand, nX, nY,
464+
oOutExtent, oCurExtent, oOpts, oProgress,
465+
/* emitWarningIfNoData = */ true);
466+
executor.run();
467+
}
468+
else
469+
{
470+
ViewshedExecutor executor(*pSrcBand, *pDstBand, nX, nY, oOutExtent,
471+
oCurExtent, oOpts, oProgress,
472+
/* emitWarningIfNoData = */ true);
473+
executor.run();
474+
}
450475
oProgress.emit(1);
451476
return static_cast<bool>(poDstDS);
452477
}

‎alg/viewshed/viewshed.h‎

Lines changed: 5 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -53,6 +53,10 @@ class Viewshed
5353
GDALProgressFunc pfnProgress = GDALDummyProgress,
5454
void *pProgressArg = nullptr);
5555

56+
CPL_DLL bool run(GDALRasterBandH hBand, GDALRasterBandH hSdBand,
57+
GDALProgressFunc pfnProgress = GDALDummyProgress,
58+
void *pProgressArg = nullptr);
59+
5660
/**
5761
* Fetch a pointer to the created raster band.
5862
*
@@ -69,6 +73,7 @@ class Viewshed
6973
Window oCurExtent{};
7074
DatasetPtr poDstDS{};
7175
GDALRasterBand *pSrcBand = nullptr;
76+
GDALRasterBand *pSdBand = nullptr;
7277

7378
DatasetPtr execute(int nX, int nY, const std::string &outFilename);
7479
void setOutput(double &dfResult, double &dfCellVal, double dfZ);

‎alg/viewshed/viewshed_executor.cpp‎

Lines changed: 92 additions & 29 deletions
Original file line numberDiff line numberDiff line change
@@ -267,6 +267,7 @@ void ViewshedExecutor::setOutputSd(Lines &lines, int pos, double dfZ)
267267
if (cellHeight > dfZ)
268268
result = oOpts.maybeVisibleVal;
269269
}
270+
270271
if (sd <= 1)
271272
cur = std::max(dfZ, cur);
272273
else
@@ -505,9 +506,8 @@ bool ViewshedExecutor::processFirstLine(Lines &lines)
505506
lines.result[m_nX] = oOpts.visibleVal;
506507
}
507508

508-
auto process = [this, &ll, &lines](bool sdCalc)
509+
auto process = [this, &ll, &lines]()
509510
{
510-
m_sdCalc = sdCalc;
511511
if (!oCurExtent.containsY(m_nY))
512512
processFirstLineTopOrBottom(ll, lines);
513513
else
@@ -517,14 +517,17 @@ bool ViewshedExecutor::processFirstLine(Lines &lines)
517517
pQueue->SubmitJob([&]() { processFirstLineRight(ll, lines); });
518518
pQueue->WaitCompletion();
519519
}
520-
sdCalc = false;
521520
};
522521

523-
process(false);
522+
process();
523+
lines.prev = lines.cur;
524524
if (sdMode())
525525
{
526+
m_sdCalc = true;
526527
lines.cur = std::move(lines.input);
527-
process(true);
528+
process();
529+
lines.prevTmp = lines.cur;
530+
m_sdCalc = false;
528531
}
529532

530533
if (oOpts.pitchMasking())
@@ -598,7 +601,12 @@ void ViewshedExecutor::processFirstLineLeft(const LineLimits &ll, Lines &lines)
598601
{
599602
double dfZ = lines.cur[iStart];
600603
if (oOpts.outputMode == OutputMode::Normal)
604+
{
601605
lines.result[iStart] = oOpts.visibleVal;
606+
if (m_sdCalc)
607+
if (lines.sd[iStart] > 1)
608+
lines.cur[iStart] = m_dfZObserver;
609+
}
602610
else
603611
setOutput(lines, iStart, dfZ);
604612
iStart--;
@@ -810,7 +818,12 @@ void ViewshedExecutor::processFirstLineRight(const LineLimits &ll, Lines &lines)
810818
{
811819
double dfZ = lines.cur[iStart];
812820
if (oOpts.outputMode == OutputMode::Normal)
821+
{
813822
lines.result[iStart] = oOpts.visibleVal;
823+
if (m_sdCalc)
824+
if (lines.sd[iStart] > 1)
825+
lines.cur[iStart] = m_dfZObserver;
826+
}
814827
else
815828
setOutput(lines, iStart, dfZ);
816829
iStart++;
@@ -916,8 +929,17 @@ void ViewshedExecutor::processLineRight(int nYOffset, LineLimits &ll,
916929
int nXOffset = std::abs(iPixel - m_nX);
917930
double dfZ;
918931
if (nXOffset == nYOffset)
932+
{
933+
if (m_sdCalc && nXOffset == 1)
934+
{
935+
lines.result[iPixel] = oOpts.visibleVal;
936+
if (lines.sd[iPixel] > 1)
937+
lines.cur[iPixel] = m_dfZObserver;
938+
continue;
939+
}
919940
dfZ = CalcHeightLine(nYOffset, lines.cur[iPixel],
920941
lines.prev[iPixel - 1]);
942+
}
921943
else
922944
dfZ = oZcalc(nXOffset, nYOffset, lines.cur[iPixel - 1],
923945
lines.prev[iPixel], lines.prev[iPixel - 1]);
@@ -927,24 +949,41 @@ void ViewshedExecutor::processLineRight(int nYOffset, LineLimits &ll,
927949
maskLineRight(lines.result, ll, nLine);
928950
}
929951

930-
/// Apply angular mask to the initial X position. Assumes m_nX is in the raster.
952+
/// Apply angular/distance mask to the initial X position. Assumes m_nX is in the raster.
931953
/// @param vResult Raster line on which to apply mask.
954+
/// @param ll Line limits.
932955
/// @param nLine Line number.
933-
void ViewshedExecutor::maskInitial(std::vector<double> &vResult, int nLine)
956+
/// @return True if the initial X position was masked.
957+
bool ViewshedExecutor::maskInitial(std::vector<double> &vResult,
958+
const LineLimits &ll, int nLine)
934959
{
960+
// Mask min/max.
961+
if (ll.left >= ll.right || ll.leftMin != ll.rightMin)
962+
{
963+
vResult[m_nX] = oOpts.outOfRangeVal;
964+
return true;
965+
}
966+
935967
if (!oOpts.angleMasking())
936-
return;
968+
return false;
937969

938970
if (nLine < m_nY)
939971
{
940972
if (!rayBetween(oOpts.startAngle, oOpts.endAngle, M_PI / 2))
973+
{
941974
vResult[m_nX] = oOpts.outOfRangeVal;
975+
return true;
976+
}
942977
}
943978
else if (nLine > m_nY)
944979
{
945980
if (!rayBetween(oOpts.startAngle, oOpts.endAngle, 3 * M_PI / 2))
981+
{
946982
vResult[m_nX] = oOpts.outOfRangeVal;
983+
return true;
984+
}
947985
}
986+
return false;
948987
}
949988

950989
/// Process a line above or below the observer.
@@ -962,37 +1001,61 @@ bool ViewshedExecutor::processLine(int nLine, Lines &lines)
9621001
// Adjust height of the read line.
9631002
LineLimits ll = adjustHeight(nYOffset, lines);
9641003

965-
// Handle the initial position on the line.
1004+
auto process = [this, nYOffset, &ll, &lines]()
1005+
{
1006+
CPLJobQueuePtr pQueue = m_pool.CreateJobQueue();
1007+
pQueue->SubmitJob([&]() { processLineLeft(nYOffset, ll, lines); });
1008+
pQueue->SubmitJob([&]() { processLineRight(nYOffset, ll, lines); });
1009+
pQueue->WaitCompletion();
1010+
};
1011+
1012+
bool masked = false;
1013+
// Handle initial position on the line.
9661014
if (oCurExtent.containsX(m_nX))
9671015
{
968-
if (ll.left < ll.right && ll.leftMin == ll.rightMin)
1016+
masked = maskInitial(lines.result, ll, nLine);
1017+
if (!masked)
9691018
{
9701019
double dfZ = CalcHeightLine(std::abs(nYOffset), lines.cur[m_nX],
9711020
lines.prev[m_nX]);
9721021
setOutput(lines, m_nX, dfZ);
9731022
}
974-
else
975-
lines.result[m_nX] = oOpts.outOfRangeVal;
976-
977-
maskInitial(lines.result, nLine);
9781023
}
9791024

980-
auto process = [this, nYOffset, &ll, &lines](bool sdCalc)
981-
{
982-
m_sdCalc = sdCalc;
983-
CPLJobQueuePtr pQueue = m_pool.CreateJobQueue();
984-
pQueue->SubmitJob([&]() { processLineLeft(nYOffset, ll, lines); });
985-
pQueue->SubmitJob([&]() { processLineRight(nYOffset, ll, lines); });
986-
pQueue->WaitCompletion();
987-
m_sdCalc = false;
988-
};
1025+
process();
9891026

990-
process(false);
1027+
// Process SD mode
9911028
if (sdMode())
9921029
{
1030+
m_sdCalc = true;
1031+
1032+
lines.prev = std::move(lines.prevTmp);
1033+
lines.prevTmp = std::move(lines.cur);
9931034
lines.cur = std::move(lines.input);
994-
process(true);
1035+
// Handle initial position on the line.
1036+
if (!masked && oCurExtent.containsX(m_nX))
1037+
{
1038+
if (std::abs(nYOffset) == 1)
1039+
{
1040+
lines.result[m_nX] = oOpts.visibleVal;
1041+
if (lines.sd[m_nX] > 1)
1042+
lines.cur[m_nX] = m_dfZObserver;
1043+
}
1044+
else
1045+
{
1046+
1047+
double dfZ = CalcHeightLine(std::abs(nYOffset), lines.cur[m_nX],
1048+
lines.prev[m_nX]);
1049+
setOutput(lines, m_nX, dfZ);
1050+
}
1051+
}
1052+
process();
1053+
lines.prev = std::move(lines.prevTmp);
1054+
lines.prevTmp = lines.cur;
1055+
m_sdCalc = false;
9951056
}
1057+
else
1058+
lines.prev = lines.cur;
9961059

9971060
if (oOpts.pitchMasking())
9981061
applyPitchMask(lines.result, lines.pitchMask);
@@ -1077,7 +1140,8 @@ bool ViewshedExecutor::run()
10771140
[&]()
10781141
{
10791142
Lines lines(oCurExtent.xSize());
1080-
lines.prev = firstLine.cur;
1143+
lines.prev = firstLine.prev;
1144+
lines.prevTmp = firstLine.prevTmp;
10811145
if (oOpts.pitchMasking())
10821146
lines.pitchMask.resize(
10831147
oOutExtent.xSize(),
@@ -1090,7 +1154,6 @@ bool ViewshedExecutor::run()
10901154
{
10911155
if (!processLine(nLine, lines))
10921156
err = true;
1093-
lines.prev = lines.cur;
10941157
if (oOpts.pitchMasking())
10951158
std::fill(lines.pitchMask.begin(), lines.pitchMask.end(),
10961159
std::numeric_limits<double>::quiet_NaN());
@@ -1102,7 +1165,8 @@ bool ViewshedExecutor::run()
11021165
[&]()
11031166
{
11041167
Lines lines(oCurExtent.xSize());
1105-
lines.prev = firstLine.cur;
1168+
lines.prev = firstLine.prev;
1169+
lines.prevTmp = firstLine.prevTmp;
11061170
if (oOpts.pitchMasking())
11071171
lines.pitchMask.resize(
11081172
oOutExtent.xSize(),
@@ -1115,7 +1179,6 @@ bool ViewshedExecutor::run()
11151179
{
11161180
if (!processLine(nLine, lines))
11171181
err = true;
1118-
lines.prev = lines.cur;
11191182
if (oOpts.pitchMasking())
11201183
std::fill(lines.pitchMask.begin(), lines.pitchMask.end(),
11211184
std::numeric_limits<double>::quiet_NaN());

‎alg/viewshed/viewshed_executor.h‎

Lines changed: 4 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -38,7 +38,8 @@ struct Lines
3838
std::vector<double>
3939
pitchMask; //!< Height/indicator values for pitch masking.
4040
std::vector<double> input; //!< Copy of input data when in SD mode.
41-
std::vector<double> sd; //!< SD mask.
41+
std::vector<double> prevTmp; //!< Saved prev values when in SD mode.
42+
std::vector<double> sd; //!< SD mask.
4243

4344
/// Constructor
4445
Lines() : cur(), result(), prev(), pitchMask(), input(), sd()
@@ -129,7 +130,8 @@ class ViewshedExecutor
129130
void processLineLeft(int nYOffset, LineLimits &ll, Lines &lines);
130131
void processLineRight(int nYOffset, LineLimits &ll, Lines &lines);
131132
LineLimits adjustHeight(int iLine, Lines &lines);
132-
void maskInitial(std::vector<double> &vResult, int nLine);
133+
bool maskInitial(std::vector<double> &vResult, const LineLimits &ll,
134+
int nLine);
133135
bool maskAngleLeft(std::vector<double> &vResult, int nLine);
134136
bool maskAngleRight(std::vector<double> &vResult, int nLine);
135137
void maskLineLeft(std::vector<double> &vResult, const LineLimits &ll,

‎apps/gdal_viewshed.cpp‎

Lines changed: 28 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -28,7 +28,9 @@ struct Options
2828
{
2929
viewshed::Options opts;
3030
std::string osSrcFilename;
31+
std::string osSdFilename;
3132
int nBandIn{1};
33+
// Currently no option for SD band number -- always 1.
3234
bool bQuiet;
3335
};
3436

@@ -160,6 +162,13 @@ Options parseArgs(GDALArgumentParser &argParser, const CPLStringList &aosArgv)
160162
.nargs(1)
161163
.help(_("Spacing between observer cells when using cumulative mode."));
162164

165+
argParser.add_argument("-sd")
166+
.store_into(localOpts.osSdFilename)
167+
.metavar("<value>")
168+
.nargs(1)
169+
.help(_("Name of file whose first band represents the standard "
170+
"deviation of the raster."));
171+
163172
argParser.add_quiet_argument(&localOpts.bQuiet);
164173

165174
argParser.add_argument("src_filename")
@@ -221,7 +230,7 @@ void validateArgs(Options &localOpts, const GDALArgumentParser &argParser)
221230

222231
if (opts.outputMode == viewshed::OutputMode::Cumulative)
223232
{
224-
for (const char *opt : {"-ox", "-oy", "-vv", "-iv", "-md"})
233+
for (const char *opt : {"-ox", "-oy", "-vv", "-iv", "-md", "-sd"})
225234
if (argParser.is_used(opt))
226235
{
227236
std::string err = "Option " + std::string(opt) +
@@ -299,6 +308,21 @@ MAIN_START(argc, argv)
299308
exit(2);
300309
}
301310

311+
/* -------------------------------------------------------------------- */
312+
/* Open source SD raster file. */
313+
/* -------------------------------------------------------------------- */
314+
GDALDatasetH hSdDS = GDALOpen(localOpts.osSdFilename.c_str(), GA_ReadOnly);
315+
if (hSrcDS == nullptr)
316+
exit(2);
317+
318+
GDALRasterBandH hSdBand = GDALGetRasterBand(hSdDS, 1);
319+
if (hSdBand == nullptr)
320+
{
321+
CPLError(CE_Failure, CPLE_AppDefined,
322+
"SD band (1) does not exist on SD dataset.");
323+
exit(2);
324+
}
325+
302326
if (!argParser.is_used("-cc"))
303327
opts.curveCoeff =
304328
gdal::viewshed::adjustCurveCoeff(opts.curveCoeff, hSrcDS);
@@ -320,8 +344,9 @@ MAIN_START(argc, argv)
320344
else
321345
{
322346
viewshed::Viewshed oViewshed(opts);
323-
bSuccess = oViewshed.run(hBand, localOpts.bQuiet ? GDALDummyProgress
324-
: GDALTermProgress);
347+
bSuccess = oViewshed.run(hBand, hSdBand,
348+
localOpts.bQuiet ? GDALDummyProgress
349+
: GDALTermProgress);
325350
hDstDS = GDALDataset::FromHandle(oViewshed.output().release());
326351
GDALClose(hSrcDS);
327352
if (GDALClose(hDstDS) != CE_None)

0 commit comments

Comments
 (0)