Skip to content

Commit c57e878

Browse files
abellgithubabellgithub
committed
Potentially working with SD plane.
1 parent 783aed0 commit c57e878

File tree

6 files changed

+161
-25
lines changed

6 files changed

+161
-25
lines changed

‎alg/viewshed/viewshed.cpp‎

Lines changed: 9 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -310,6 +310,13 @@ bool Viewshed::calcExtents(int nX, int nY,
310310
bool Viewshed::run(GDALRasterBandH band, GDALProgressFunc pfnProgress,
311311
void *pProgressArg)
312312
{
313+
DatasetPtr sdDataset;
314+
if (oOpts.sdFilename.size())
315+
{
316+
GDALDatasetH ds = GDALOpen(oOpts.sdFilename.c_str(), GA_ReadOnly);
317+
sdDataset.reset(GDALDataset::FromHandle(ds));
318+
pSdBand = sdDataset->GetRasterBand(1);
319+
}
313320
pSrcBand = static_cast<GDALRasterBand *>(band);
314321

315322
std::array<double, 6> adfFwdTransform;
@@ -348,8 +355,8 @@ bool Viewshed::run(GDALRasterBandH band, GDALProgressFunc pfnProgress,
348355

349356
// Execute the viewshed algorithm.
350357
GDALRasterBand *pDstBand = poDstDS->GetRasterBand(1);
351-
ViewshedExecutor executor(*pSrcBand, *pDstBand, nX, nY, oOutExtent,
352-
oCurExtent, oOpts, oProgress);
358+
ViewshedExecutor executor(*pSrcBand, *pSdBand, *pDstBand, nX, nY,
359+
oOutExtent, oCurExtent, oOpts, oProgress);
353360
executor.run();
354361
oProgress.emit(1);
355362
return static_cast<bool>(poDstDS);

‎alg/viewshed/viewshed.h‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -87,6 +87,7 @@ class Viewshed
8787
Window oCurExtent;
8888
DatasetPtr poDstDS;
8989
GDALRasterBand *pSrcBand;
90+
GDALRasterBand *pSdBand;
9091

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

‎alg/viewshed/viewshed_executor.cpp‎

Lines changed: 131 additions & 19 deletions
Original file line numberDiff line numberDiff line change
@@ -32,6 +32,9 @@
3232
#include <cassert>
3333
#include <limits>
3434

35+
//ABELL
36+
#include <iostream>
37+
3538
#include "viewshed_executor.h"
3639
#include "progress.h"
3740

@@ -111,6 +114,7 @@ double doMax(int nXOffset, int nYOffset, double dfThisPrev, double dfLast,
111114

112115
/// Constructor -- the viewshed algorithm executor
113116
/// @param srcBand Source raster band
117+
/// @param sdBand Destination raster band
114118
/// @param dstBand Destination raster band
115119
/// @param nX X position of observer
116120
/// @param nY Y position of observer
@@ -119,13 +123,30 @@ double doMax(int nXOffset, int nYOffset, double dfThisPrev, double dfLast,
119123
/// @param opts Configuration options.
120124
/// @param progress Reference to the progress tracker.
121125
ViewshedExecutor::ViewshedExecutor(GDALRasterBand &srcBand,
126+
GDALRasterBand &sdBand,
122127
GDALRasterBand &dstBand, int nX, int nY,
123128
const Window &outExtent,
124129
const Window &curExtent, const Options &opts,
125130
Progress &progress)
126-
: m_pool(4), m_srcBand(srcBand), m_dstBand(dstBand), oOutExtent(outExtent),
127-
oCurExtent(curExtent), m_nX(nX - oOutExtent.xStart), m_nY(nY),
128-
oOpts(opts), oProgress(progress),
131+
: m_pool(4), m_srcBand(srcBand), m_sdBand(sdBand), m_dstBand(dstBand),
132+
oOutExtent(outExtent), oCurExtent(curExtent),
133+
m_nX(nX - oOutExtent.xStart), m_nY(nY), oOpts(opts), oProgress(progress),
134+
m_dfMaxDistance2(opts.maxDistance * opts.maxDistance)
135+
{
136+
if (m_dfMaxDistance2 == 0)
137+
m_dfMaxDistance2 = std::numeric_limits<double>::max();
138+
m_srcBand.GetDataset()->GetGeoTransform(m_adfTransform.data());
139+
}
140+
141+
// No SD Band
142+
ViewshedExecutor::ViewshedExecutor(GDALRasterBand &srcBand,
143+
GDALRasterBand &dstBand, int nX, int nY,
144+
const Window &outExtent,
145+
const Window &curExtent, const Options &opts,
146+
Progress &progress)
147+
: m_pool(4), m_srcBand(srcBand), m_sdBand(m_srcBand), m_dstBand(dstBand),
148+
oOutExtent(outExtent), oCurExtent(curExtent),
149+
m_nX(nX - oOutExtent.xStart), m_nY(nY), oOpts(opts), oProgress(progress),
129150
m_dfMaxDistance2(opts.maxDistance * opts.maxDistance)
130151
{
131152
if (m_dfMaxDistance2 == 0)
@@ -162,17 +183,48 @@ double ViewshedExecutor::calcHeightAdjFactor()
162183
/// dfCellVal Reference to the current cell height. Replace with observable height.
163184
/// dfZ Minimum observable height at cell.
164185
void ViewshedExecutor::setOutput(double &dfResult, double &dfCellVal,
165-
double dfZ)
186+
double &dfSdCellVal, double dfZ, double dfSdZ)
166187
{
167188
if (oOpts.outputMode != OutputMode::Normal)
168189
{
169190
dfResult += (dfZ - dfCellVal);
170191
dfResult = std::max(0.0, dfResult);
192+
dfSdCellVal = std::max(dfSdCellVal, dfSdZ);
171193
}
172194
else
173-
dfResult = (dfCellVal + oOpts.targetHeight < dfZ) ? oOpts.invisibleVal
174-
: oOpts.visibleVal;
175-
dfCellVal = std::max(dfCellVal, dfZ);
195+
{
196+
/**
197+
std::cerr << "cell height/dfz/dfsdz/sd = " <<
198+
(dfCellVal + oOpts.targetHeight) << "/" << dfZ << "/" << dfSdZ << "/" <<
199+
dfSdCellVal;
200+
**/
201+
if (dfZ <= dfCellVal + oOpts.targetHeight)
202+
{
203+
// std::cerr << " - visible!\n";
204+
dfResult = oOpts.visibleVal;
205+
if (dfSdCellVal > 1) // Changing SD to height.
206+
dfSdCellVal = dfSdZ; // See-through case.
207+
else
208+
dfSdCellVal = dfCellVal; // Blocking case
209+
}
210+
else if (dfSdZ <= dfCellVal + oOpts.targetHeight)
211+
{
212+
// std::cerr << " - maybe visible!\n";
213+
dfResult = oOpts.maybeVisibleVal;
214+
if (dfSdCellVal > 1) // Changing SD to height.
215+
dfSdCellVal = dfSdZ; // See-through case.
216+
else
217+
dfSdCellVal = dfCellVal; // Blocking case
218+
dfCellVal = dfZ;
219+
}
220+
else
221+
{
222+
// std::cerr << " - not visible!\n";
223+
dfResult = oOpts.invisibleVal;
224+
dfCellVal = dfZ;
225+
dfSdCellVal = dfSdZ;
226+
}
227+
}
176228
}
177229

178230
/// Read a line of raster data.
@@ -199,9 +251,29 @@ bool ViewshedExecutor::readLine(int nLine, std::vector<Cell> &vThisLine)
199251
oOutExtent.xStart, nLine, oOutExtent.xSize(), 1);
200252
return false;
201253
}
202-
203254
for (size_t i = 0; i < line.size(); ++i)
204255
vThisLine[i].val = line[i];
256+
257+
int sdStatus = 0;
258+
double nodata = std::numeric_limits<double>::max();
259+
{
260+
std::lock_guard g(iMutex);
261+
nodata = m_sdBand.GetNoDataValue();
262+
sdStatus = m_sdBand.RasterIO(
263+
GF_Read, oOutExtent.xStart, nLine, oOutExtent.xSize(), 1,
264+
line.data(), oOutExtent.xSize(), 1, GDT_Float64, 0, 0, nullptr);
265+
}
266+
if (sdStatus)
267+
{
268+
CPLError(CE_Failure, CPLE_AppDefined,
269+
"RasterIO error when reading SD band at position (%d,%d), "
270+
"size (%d,%d)",
271+
oOutExtent.xStart, nLine, oOutExtent.xSize(), 1);
272+
return false;
273+
}
274+
for (size_t i = 0; i < line.size(); ++i)
275+
vThisLine[i].sd = (line[i] == nodata ? 1000.0 : line[i]);
276+
205277
return true;
206278
}
207279

@@ -358,7 +430,6 @@ bool ViewshedExecutor::processFirstLine(std::vector<Cell> &vLastLine)
358430
pQueue->SubmitJob(
359431
[&, left = iLeft]()
360432
{ processFirstLineLeft(m_nX - 1, left - 1, vThisLine); });
361-
362433
pQueue->SubmitJob(
363434
[&, right = iRight]()
364435
{ processFirstLineRight(m_nX + 1, right, vThisLine); });
@@ -390,7 +461,8 @@ void ViewshedExecutor::processFirstLineTopOrBottom(int iLeft, int iRight,
390461
if (oOpts.outputMode == OutputMode::Normal)
391462
c.result = oOpts.visibleVal;
392463
else
393-
setOutput(c.result, c.val, c.val);
464+
//ABELL - Check SD height here. Perhaps 0 if SD < 1?
465+
setOutput(c.result, c.val, c.sd, c.val, c.val);
394466
}
395467
for (int i = 0; i < iLeft; ++i)
396468
vThisLine[i].result = oOpts.outOfRangeVal;
@@ -417,9 +489,12 @@ void ViewshedExecutor::processFirstLineLeft(int iStart, int iEnd,
417489
{
418490
Cell &c = vThisLine[iStart];
419491
if (oOpts.outputMode == OutputMode::Normal)
492+
{
420493
c.result = oOpts.visibleVal;
494+
c.sd = c.val; // Set the SD height to the height.
495+
}
421496
else
422-
setOutput(c.result, c.val, c.val);
497+
setOutput(c.result, c.val, c.sd, c.val, c.val);
423498
iStart--;
424499
}
425500

@@ -429,7 +504,8 @@ void ViewshedExecutor::processFirstLineLeft(int iStart, int iEnd,
429504
Cell &c = vThisLine[iPixel];
430505
int nXOffset = std::abs(iPixel - m_nX);
431506
double dfZ = CalcHeightLine(nXOffset, vThisLine[iPixel + 1].val);
432-
setOutput(c.result, c.val, dfZ);
507+
double dfSdZ = CalcHeightLine(nXOffset, vThisLine[iPixel + 1].sd);
508+
setOutput(c.result, c.val, c.sd, dfZ, dfSdZ);
433509
}
434510
// For cells outside of the [start, end) range, set the outOfRange value.
435511
for (int i = 0; i < iEnd + 1; ++i)
@@ -459,7 +535,8 @@ void ViewshedExecutor::processFirstLineRight(int iStart, int iEnd,
459535
if (oOpts.outputMode == OutputMode::Normal)
460536
c.result = oOpts.visibleVal;
461537
else
462-
setOutput(c.result, c.val, c.val);
538+
//ABELL Perhaps set the last value to zero if SD < 1.
539+
setOutput(c.result, c.val, c.sd, c.val, c.val);
463540
iStart++;
464541
}
465542

@@ -469,7 +546,8 @@ void ViewshedExecutor::processFirstLineRight(int iStart, int iEnd,
469546
Cell &c = vThisLine[iPixel];
470547
int nXOffset = std::abs(iPixel - m_nX);
471548
double dfZ = CalcHeightLine(nXOffset, vThisLine[iPixel - 1].val);
472-
setOutput(c.result, c.val, dfZ);
549+
double dfSdZ = CalcHeightLine(nXOffset, vThisLine[iPixel - 1].sd);
550+
setOutput(c.result, c.val, c.sd, dfZ, dfSdZ);
473551
}
474552
// For cells outside of the [start, end) range, set the outOfRange value.
475553
for (int i = iEnd; i < oOutExtent.xSize(); ++i)
@@ -503,7 +581,8 @@ void ViewshedExecutor::processLineLeft(int nYOffset, int iStart, int iEnd,
503581
if (oOpts.outputMode == OutputMode::Normal)
504582
c.result = oOpts.visibleVal;
505583
else
506-
setOutput(c.result, c.val, c.val);
584+
//ABELL Perhaps set the last value to zero if SD < 1.
585+
setOutput(c.result, c.val, c.sd, c.val, c.val);
507586
iStart--;
508587
}
509588

@@ -512,17 +591,29 @@ void ViewshedExecutor::processLineLeft(int nYOffset, int iStart, int iEnd,
512591
{
513592
int nXOffset = std::abs(iPixel - m_nX);
514593
double dfZ;
594+
double dfSdZ;
515595
if (nXOffset == nYOffset)
516596
{
517597
if (nXOffset == 1)
598+
{
518599
dfZ = vThisLine[iPixel].val;
600+
dfSdZ = vThisLine[iPixel].val; //ABELL - Check this.
601+
}
519602
else
603+
{
520604
dfZ = CalcHeightLine(nXOffset, vLastLine[iPixel + 1].val);
605+
dfSdZ = CalcHeightLine(nXOffset, vLastLine[iPixel + 1].sd);
606+
}
521607
}
522608
else
609+
{
523610
dfZ = oZcalc(nXOffset, nYOffset, vThisLine[iPixel + 1].val,
524611
vLastLine[iPixel].val, vLastLine[iPixel + 1].val);
525-
setOutput(vThisLine[iPixel].result, vThisLine[iPixel].val, dfZ);
612+
dfSdZ = oZcalc(nXOffset, nYOffset, vThisLine[iPixel + 1].sd,
613+
vLastLine[iPixel].sd, vLastLine[iPixel + 1].sd);
614+
}
615+
setOutput(vThisLine[iPixel].result, vThisLine[iPixel].val,
616+
vThisLine[iPixel].sd, dfZ, dfSdZ);
526617
}
527618

528619
// For cells outside of the [start, end) range, set the outOfRange value.
@@ -557,7 +648,8 @@ void ViewshedExecutor::processLineRight(int nYOffset, int iStart, int iEnd,
557648
if (oOpts.outputMode == OutputMode::Normal)
558649
c.result = oOpts.visibleVal;
559650
else
560-
setOutput(c.result, c.val, c.val);
651+
//ABELL - Here too.
652+
setOutput(c.result, c.val, c.sd, c.val, c.val);
561653
iStart++;
562654
}
563655

@@ -566,17 +658,29 @@ void ViewshedExecutor::processLineRight(int nYOffset, int iStart, int iEnd,
566658
{
567659
int nXOffset = std::abs(iPixel - m_nX);
568660
double dfZ;
661+
double dfSdZ;
569662
if (nXOffset == nYOffset)
570663
{
571664
if (nXOffset == 1)
665+
{
572666
dfZ = vThisLine[iPixel].val;
667+
dfSdZ = vThisLine[iPixel].val; //ABELL - Here.
668+
}
573669
else
670+
{
574671
dfZ = CalcHeightLine(nXOffset, vLastLine[iPixel - 1].val);
672+
dfSdZ = CalcHeightLine(nXOffset, vLastLine[iPixel - 1].sd);
673+
}
575674
}
576675
else
676+
{
577677
dfZ = oZcalc(nXOffset, nYOffset, vThisLine[iPixel - 1].val,
578678
vLastLine[iPixel].val, vLastLine[iPixel - 1].val);
579-
setOutput(vThisLine[iPixel].result, vThisLine[iPixel].val, dfZ);
679+
dfSdZ = oZcalc(nXOffset, nYOffset, vThisLine[iPixel - 1].sd,
680+
vLastLine[iPixel].sd, vLastLine[iPixel - 1].sd);
681+
}
682+
setOutput(vThisLine[iPixel].result, vThisLine[iPixel].val,
683+
vThisLine[iPixel].sd, dfZ, dfSdZ);
580684
}
581685

582686
// For cells outside of the [start, end) range, set the outOfRange value.
@@ -613,11 +717,19 @@ bool ViewshedExecutor::processLine(int nLine, std::vector<Cell> &vLastLine)
613717
if (iLeft < iRight)
614718
{
615719
double dfZ;
720+
double dfSdZ;
616721
if (std::abs(nYOffset) == 1)
722+
{
617723
dfZ = vThisLine[m_nX].val;
724+
dfSdZ = vThisLine[m_nX].val;
725+
}
618726
else
727+
{
619728
dfZ = CalcHeightLine(nYOffset, vLastLine[m_nX].val);
620-
setOutput(vThisLine[m_nX].result, vThisLine[m_nX].val, dfZ);
729+
dfSdZ = CalcHeightLine(nYOffset, vLastLine[m_nX].sd);
730+
}
731+
setOutput(vThisLine[m_nX].result, vThisLine[m_nX].val,
732+
vThisLine[m_nX].sd, dfZ, dfSdZ);
621733
}
622734
else
623735
vThisLine[m_nX].result = oOpts.outOfRangeVal;

‎alg/viewshed/viewshed_executor.h‎

Lines changed: 8 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -52,9 +52,14 @@ class ViewshedExecutor
5252
{
5353
double val{};
5454
double result{};
55+
double sd{};
5556
};
5657

5758
public:
59+
ViewshedExecutor(GDALRasterBand &srcBand, GDALRasterBand &sdBand,
60+
GDALRasterBand &dstBand, int nX, int nY,
61+
const Window &oOutExtent, const Window &oCurExtent,
62+
const Options &opts, Progress &oProgress);
5863
ViewshedExecutor(GDALRasterBand &srcBand, GDALRasterBand &dstBand, int nX,
5964
int nY, const Window &oOutExtent, const Window &oCurExtent,
6065
const Options &opts, Progress &oProgress);
@@ -63,6 +68,7 @@ class ViewshedExecutor
6368
private:
6469
CPLWorkerThreadPool m_pool;
6570
GDALRasterBand &m_srcBand;
71+
GDALRasterBand &m_sdBand;
6672
GDALRasterBand &m_dstBand;
6773
const Window oOutExtent;
6874
const Window oCurExtent;
@@ -79,7 +85,8 @@ class ViewshedExecutor
7985
double (*oZcalc)(int, int, double, double, double){};
8086

8187
double calcHeightAdjFactor();
82-
void setOutput(double &dfResult, double &dfCellVal, double dfZ);
88+
void setOutput(double &dfResult, double &dfCellVal, double &sdCellVal,
89+
double dfZ, double dfSdZ);
8390
bool readLine(int nLine, std::vector<Cell> &vThisLine);
8491
bool writeLine(int nLine, std::vector<Cell> &vThisLine);
8592
bool processLine(int nLine, std::vector<Cell> &vLastLine);

‎alg/viewshed/viewshed_types.h‎

Lines changed: 5 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -71,9 +71,10 @@ struct Point
7171
*/
7272
struct Options
7373
{
74-
Point observer{0, 0, 0}; //!< x, y, and z of the observer
75-
double visibleVal{255}; //!< raster output value for visible pixels.
76-
double invisibleVal{0}; //!< raster output value for non-visible pixels.
74+
Point observer{0, 0, 0}; //!< x, y, and z of the observer
75+
double visibleVal{255}; //!< raster output value for visible pixels.
76+
double invisibleVal{0}; //!< raster output value for non-visible pixels.
77+
double maybeVisibleVal{2}; //!< raster output value for non-visible pixels.
7778
double outOfRangeVal{
7879
0}; //!< raster output value for pixels outside of max distance.
7980
double nodataVal{-1}; //!< raster output value for pixels with no data
@@ -83,6 +84,7 @@ struct Options
8384
double curveCoeff{.85714}; //!< coefficient for atmospheric refraction
8485
OutputMode outputMode{OutputMode::Normal}; //!< Output information.
8586
//!< Normal, Height from DEM or Height from ground
87+
std::string sdFilename{}; //!< sd raster filename
8688
std::string outputFormat{}; //!< output raster format
8789
std::string outputFilename{}; //!< output raster filename
8890
CPLStringList creationOpts{}; //!< options for output raster creation

0 commit comments

Comments
 (0)