Skip to content

Commit c9af480

Browse files
abellgithubabellgithub
committed
Add pitch masking.
1 parent 91d4b45 commit c9af480

File tree

5 files changed

+175
-14
lines changed

5 files changed

+175
-14
lines changed

‎alg/viewshed/viewshed.cpp‎

Lines changed: 6 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -404,22 +404,22 @@ bool Viewshed::run(GDALRasterBandH band, GDALProgressFunc pfnProgress,
404404
"End angle out of range. Must be [0, 360).");
405405
return false;
406406
}
407-
if (oOpts.maxPitch > 90)
407+
if (oOpts.highPitch > 90)
408408
{
409409
CPLError(CE_Failure, CPLE_AppDefined,
410-
"Invalid maxPitch. Cannot be greater than 90.");
410+
"Invalid highPitch. Cannot be greater than 90.");
411411
return false;
412412
}
413-
if (oOpts.minPitch < -90)
413+
if (oOpts.lowPitch < -90)
414414
{
415415
CPLError(CE_Failure, CPLE_AppDefined,
416-
"Invalid maxPitch. Cannot be less than -90.");
416+
"Invalid lowPitch. Cannot be less than -90.");
417417
return false;
418418
}
419-
if (oOpts.maxPitch <= oOpts.minPitch)
419+
if (oOpts.highPitch <= oOpts.lowPitch)
420420
{
421421
CPLError(CE_Failure, CPLE_AppDefined,
422-
"Invalid pitch. maxPitch must be > minPitch");
422+
"Invalid pitch. highPitch must be > lowPitch");
423423
return false;
424424
}
425425

‎alg/viewshed/viewshed_executor.cpp‎

Lines changed: 53 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -130,6 +130,10 @@ ViewshedExecutor::ViewshedExecutor(GDALRasterBand &srcBand,
130130
{
131131
if (m_dfMaxDistance2 == 0)
132132
m_dfMaxDistance2 = std::numeric_limits<double>::max();
133+
if (opts.lowPitch != -90.0)
134+
m_lowTanPitch = std::tan(oOpts.lowPitch * (2 * M_PI / 360.0));
135+
if (opts.highPitch != 90.0)
136+
m_highTanPitch = std::tan(oOpts.highPitch * (2 * M_PI / 360.0));
133137
m_srcBand.GetDataset()->GetGeoTransform(m_adfTransform.data());
134138
int hasNoData = false;
135139
m_noDataValue = m_srcBand.GetNoDataValue(&hasNoData);
@@ -427,7 +431,11 @@ void ViewshedExecutor::processFirstLineLeft(const LineLimits &ll,
427431
if (oOpts.outputMode == OutputMode::Normal)
428432
vResult[iStart] = oOpts.visibleVal;
429433
else
434+
{
435+
maskLowPitch(*pThis, 1, 0);
430436
setOutput(vResult[iStart], *pThis, *pThis);
437+
}
438+
maskHighPitch(vResult[iStart], *pThis, 1, 0);
431439
iStart--;
432440
pThis--;
433441
}
@@ -437,7 +445,9 @@ void ViewshedExecutor::processFirstLineLeft(const LineLimits &ll,
437445
{
438446
int nXOffset = std::abs(iPixel - m_nX);
439447
double dfZ = CalcHeightLine(nXOffset, *(pThis + 1));
448+
maskLowPitch(dfZ, nXOffset, 0);
440449
setOutput(vResult[iPixel], *pThis, dfZ);
450+
maskHighPitch(vResult[iPixel], dfZ, nXOffset, 0);
441451
}
442452

443453
maskLineLeft(vResult, ll, m_nY);
@@ -615,7 +625,11 @@ void ViewshedExecutor::processFirstLineRight(const LineLimits &ll,
615625
if (oOpts.outputMode == OutputMode::Normal)
616626
vResult[iStart] = oOpts.visibleVal;
617627
else
628+
{
629+
maskLowPitch(*pThis, 1, 0);
618630
setOutput(vResult[iStart], *pThis, *pThis);
631+
}
632+
maskHighPitch(vResult[iStart], *pThis, 1, 0);
619633
iStart++;
620634
pThis++;
621635
}
@@ -625,7 +639,9 @@ void ViewshedExecutor::processFirstLineRight(const LineLimits &ll,
625639
{
626640
int nXOffset = std::abs(iPixel - m_nX);
627641
double dfZ = CalcHeightLine(nXOffset, *(pThis - 1));
642+
maskLowPitch(dfZ, nXOffset, 0);
628643
setOutput(vResult[iPixel], *pThis, dfZ);
644+
maskHighPitch(vResult[iPixel], dfZ, nXOffset, 0);
629645
}
630646

631647
maskLineRight(vResult, ll, m_nY);
@@ -664,7 +680,11 @@ void ViewshedExecutor::processLineLeft(int nYOffset, LineLimits &ll,
664680
if (oOpts.outputMode == OutputMode::Normal)
665681
vResult[iStart] = oOpts.visibleVal;
666682
else
683+
{
684+
maskLowPitch(*pThis, m_nX - iStart, nYOffset);
667685
setOutput(vResult[iStart], *pThis, *pThis);
686+
maskHighPitch(vResult[iStart], *pThis, m_nX - iStart, nYOffset);
687+
}
668688
iStart--;
669689
pThis--;
670690
pLast--;
@@ -685,7 +705,9 @@ void ViewshedExecutor::processLineLeft(int nYOffset, LineLimits &ll,
685705
else
686706
dfZ =
687707
oZcalc(nXOffset, nYOffset, *(pThis + 1), *pLast, *(pLast + 1));
708+
maskLowPitch(dfZ, nXOffset, nYOffset);
688709
setOutput(vResult[iPixel], *pThis, dfZ);
710+
maskHighPitch(vResult[iPixel], dfZ, nXOffset, nYOffset);
689711
}
690712

691713
maskLineLeft(vResult, ll, nLine);
@@ -724,7 +746,11 @@ void ViewshedExecutor::processLineRight(int nYOffset, LineLimits &ll,
724746
if (oOpts.outputMode == OutputMode::Normal)
725747
vResult[iStart] = oOpts.visibleVal;
726748
else
749+
{
750+
maskLowPitch(*pThis, m_nX, nYOffset);
727751
setOutput(vResult[0], *pThis, *pThis);
752+
maskHighPitch(vResult[0], *pThis, m_nX, nYOffset);
753+
}
728754
iStart++;
729755
pThis++;
730756
pLast++;
@@ -745,7 +771,9 @@ void ViewshedExecutor::processLineRight(int nYOffset, LineLimits &ll,
745771
else
746772
dfZ =
747773
oZcalc(nXOffset, nYOffset, *(pThis - 1), *pLast, *(pLast - 1));
774+
maskLowPitch(dfZ, nXOffset, nYOffset);
748775
setOutput(vResult[iPixel], *pThis, dfZ);
776+
maskHighPitch(vResult[iPixel], dfZ, nXOffset, nYOffset);
749777
}
750778

751779
maskLineRight(vResult, ll, nLine);
@@ -802,7 +830,9 @@ bool ViewshedExecutor::processLine(int nLine, std::vector<double> &vLastLineVal)
802830
dfZ = vThisLineVal[m_nX];
803831
else
804832
dfZ = CalcHeightLine(nYOffset, vLastLineVal[m_nX]);
833+
maskLowPitch(dfZ, 0, nYOffset);
805834
setOutput(vResult[m_nX], vThisLineVal[m_nX], dfZ);
835+
maskHighPitch(vResult[m_nX], dfZ, 0, nYOffset);
806836
}
807837
else
808838
vResult[m_nX] = oOpts.outOfRangeVal;
@@ -919,5 +949,28 @@ bool ViewshedExecutor::run()
919949
return true;
920950
}
921951

952+
void ViewshedExecutor::maskLowPitch(double &dfZ, int nXOffset, int nYOffset)
953+
{
954+
if (std::isnan(m_lowTanPitch))
955+
return;
956+
957+
double dfDist = std::sqrt(nXOffset * nXOffset + nYOffset * nYOffset);
958+
double dfZmask = dfDist * m_lowTanPitch;
959+
if (dfZmask > dfZ)
960+
dfZ = dfZmask;
961+
}
962+
963+
void ViewshedExecutor::maskHighPitch(double &dfResult, double dfZ, int nXOffset,
964+
int nYOffset)
965+
{
966+
if (std::isnan(m_highTanPitch))
967+
return;
968+
969+
double dfDist = std::sqrt(nXOffset * nXOffset + nYOffset * nYOffset);
970+
double dfZmask = dfDist * m_highTanPitch;
971+
if (dfZmask < dfZ)
972+
dfResult = oOpts.outOfRangeVal;
973+
}
974+
922975
} // namespace viewshed
923976
} // namespace gdal

‎alg/viewshed/viewshed_executor.h‎

Lines changed: 5 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -67,6 +67,8 @@ class ViewshedExecutor
6767
std::mutex oMutex{};
6868
std::array<double, 6> m_adfTransform{0, 1, 0, 0, 0, 1};
6969
std::array<double, 5> m_testAngle{};
70+
double m_lowTanPitch{std::numeric_limits<double>::quiet_NaN()};
71+
double m_highTanPitch{std::numeric_limits<double>::quiet_NaN()};
7072
double (*oZcalc)(int, int, double, double, double){};
7173

7274
double calcHeightAdjFactor();
@@ -100,6 +102,9 @@ class ViewshedExecutor
100102
int nLine);
101103
void maskLineRight(std::vector<double> &vResult, const LineLimits &ll,
102104
int nLine);
105+
void maskLowPitch(double &dfZ, int nXOffset, int nYOffset);
106+
void maskHighPitch(double &dfResult, double dfZ, int nXOffset,
107+
int nYOffset);
103108
void calcTestAngles();
104109
};
105110

‎alg/viewshed/viewshed_types.h‎

Lines changed: 3 additions & 8 deletions
Original file line numberDiff line numberDiff line change
@@ -70,10 +70,10 @@ struct Options
7070
0.0}; //!< minimum distance from observer to compute value.
7171
double startAngle{0.0}; //!< start angle of observable range
7272
double endAngle{0.0}; //!< end angle of observable range
73-
double maxPitch{
74-
90.0}; //!< maximum pitch (vertical angle) of observable points
75-
double minPitch{
73+
double lowPitch{
7674
-90.0}; //!< minimum pitch (vertical angle) of observable points
75+
double highPitch{
76+
90.0}; //!< maximum pitch (vertical angle) of observable points
7777
double curveCoeff{.85714}; //!< coefficient for atmospheric refraction
7878
OutputMode outputMode{OutputMode::Normal}; //!< Output information.
7979
//!< Normal, Height from DEM or Height from ground
@@ -88,11 +88,6 @@ struct Options
8888
{
8989
return startAngle != endAngle;
9090
}
91-
92-
bool pitchMasking() const
93-
{
94-
return maxPitch != 90 || minPitch != -90;
95-
}
9691
};
9792

9893
/**

‎autotest/cpp/test_viewshed.cpp‎

Lines changed: 108 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -118,6 +118,114 @@ TEST(Viewshed, min_max_mask)
118118
}
119119
**/
120120

121+
/**
122+
TEST(Viewshed, high_mask)
123+
{
124+
const int xlen = 30;
125+
const int ylen = 30;
126+
std::array<int8_t, xlen * ylen> in;
127+
in.fill(0);
128+
in[465] = 1;
129+
in[469] = 5;
130+
in[470] = 7;
131+
in[471] = 9;
132+
in[472] = 11;
133+
in[473] = 13;
134+
in[474] = 15;
135+
in[475] = 17;
136+
in[476] = 19;
137+
in[477] = 21;
138+
in[478] = 23;
139+
in[479] = 25;
140+
141+
SCOPED_TRACE("min_max_mask");
142+
Options opts(stdOptions(15, 15));
143+
144+
opts.highPitch = 58;
145+
146+
DatasetPtr output = runViewshed(in.data(), xlen, ylen, opts);
147+
148+
std::array<uint8_t, xlen * ylen> out;
149+
GDALRasterBand *band = output->GetRasterBand(1);
150+
151+
int xOutLen = band->GetXSize();
152+
int yOutLen = band->GetYSize();
153+
std::cerr << "Xsize = " << band->GetXSize() << "!\n";
154+
std::cerr << "Ysize = " << band->GetYSize() << "!\n";
155+
CPLErr err = band->RasterIO(GF_Read, 0, 0, xOutLen, yOutLen, out.data(), xOutLen,
156+
yOutLen, GDT_Int8, 0, 0, nullptr);
157+
EXPECT_EQ(err, CE_None);
158+
159+
uint8_t *p = out.data();
160+
for (int y = 0; y < yOutLen; ++y)
161+
{
162+
for (int x = 0; x < xOutLen; ++x)
163+
{
164+
char c;
165+
if (*p == 0)
166+
c = '*';
167+
else if (*p == 127)
168+
c = '.';
169+
else
170+
c = '?';
171+
std::cerr << c;
172+
p++;
173+
}
174+
std::cerr << "\n";
175+
}
176+
std::cerr << "\n";
177+
}
178+
**/
179+
180+
/**
181+
TEST(Viewshed, low_mask)
182+
{
183+
const int xlen = 30;
184+
const int ylen = 30;
185+
std::array<int8_t, xlen * ylen> in;
186+
in.fill(0);
187+
in[465] = 10;
188+
189+
SCOPED_TRACE("min_max_mask");
190+
Options opts(stdOptions(15, 15));
191+
192+
opts.lowPitch = -45;
193+
opts.outputMode = OutputMode::DEM;
194+
195+
DatasetPtr output = runViewshed(in.data(), xlen, ylen, opts);
196+
197+
std::array<int8_t, xlen * ylen> out;
198+
GDALRasterBand *band = output->GetRasterBand(1);
199+
200+
int xOutLen = band->GetXSize();
201+
int yOutLen = band->GetYSize();
202+
std::cerr << "Xsize = " << band->GetXSize() << "!\n";
203+
std::cerr << "Ysize = " << band->GetYSize() << "!\n";
204+
CPLErr err = band->RasterIO(GF_Read, 0, 0, xOutLen, yOutLen, out.data(), xOutLen,
205+
yOutLen, GDT_Int8, 0, 0, nullptr);
206+
EXPECT_EQ(err, CE_None);
207+
208+
int8_t *p = out.data();
209+
for (int y = 0; y < yOutLen; ++y)
210+
{
211+
for (int x = 0; x < xOutLen; ++x)
212+
{
213+
char c;
214+
if (*p > 0)
215+
c = '+';
216+
else if (*p < -9)
217+
c = '-';
218+
else
219+
c = '0' + -(*p);
220+
std::cerr << c;
221+
p++;
222+
}
223+
std::cerr << "\n";
224+
}
225+
std::cerr << "\n";
226+
}
227+
**/
228+
121229
TEST(Viewshed, all_visible)
122230
{
123231
// clang-format off

0 commit comments

Comments
 (0)