Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Contour: include minimum raster value as contour line when it exactly matches the first level #10170

Merged
merged 1 commit into from
Jun 18, 2024
Merged

Contour: include minimum raster value as contour line when it exactly matches the first level #10170

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
63 changes: 41 additions & 22 deletions alg/contour.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -43,6 +43,8 @@
#include "ogr_srs_api.h"
#include "ogr_geometry.h"

#include <limits>

static CPLErr OGRPolygonContourWriter(double dfLevelMin, double dfLevelMax,
const OGRMultiPolygon &multipoly,
void *pInfo)
Expand Down Expand Up @@ -670,24 +672,25 @@ CPLErr GDALContourGenerateEx(GDALRasterBandH hBand, void *hLayer,
GDALGetGeoTransform(hSrcDS, oCWI.adfGeoTransform);
oCWI.nNextID = 0;

int bSuccessMin = FALSE;
double dfMinimum = GDALGetRasterMinimum(hBand, &bSuccessMin);
int bSuccessMax = FALSE;
double dfMaximum = GDALGetRasterMaximum(hBand, &bSuccessMax);
if ((!bSuccessMin || !bSuccessMax))
{
double adfMinMax[2];
if (GDALComputeRasterMinMax(hBand, false, adfMinMax) == CE_None)
{
dfMinimum = adfMinMax[0];
dfMaximum = adfMinMax[1];
}
}

bool ok = false;
try
{
if (polygonize)
{
int bSuccessMin = FALSE;
double dfMinimum = GDALGetRasterMinimum(hBand, &bSuccessMin);
int bSuccessMax = FALSE;
double dfMaximum = GDALGetRasterMaximum(hBand, &bSuccessMax);
if ((!bSuccessMin || !bSuccessMax) && !fixedLevels.empty())
{
double adfMinMax[2];
if (GDALComputeRasterMinMax(hBand, false, adfMinMax) == CE_None)
{
dfMinimum = adfMinMax[0];
dfMaximum = adfMinMax[1];
}
}
if (!fixedLevels.empty())
{
// If the minimum raster value is larger than the first requested
Expand All @@ -711,8 +714,12 @@ CPLErr GDALContourGenerateEx(GDALRasterBandH hBand, void *hLayer,
RingAppender appender(w);
if (!fixedLevels.empty())
{
FixedLevelRangeIterator levels(&fixedLevels[0],
fixedLevels.size(), dfMaximum);
// Do not provide the actual minimum value to level iterator
// in polygonal case, otherwise it can result in a polygon
// with a degenerate min=max range.
FixedLevelRangeIterator levels(
&fixedLevels[0], fixedLevels.size(),
-std::numeric_limits<double>::infinity(), dfMaximum);
SegmentMerger<RingAppender, FixedLevelRangeIterator> writer(
appender, levels, /* polygonize */ true);
ContourGeneratorFromRaster<decltype(writer),
Expand All @@ -722,7 +729,11 @@ CPLErr GDALContourGenerateEx(GDALRasterBandH hBand, void *hLayer,
}
else if (expBase > 0.0)
{
ExponentialLevelRangeIterator levels(expBase);
// Do not provide the actual minimum value to level iterator
// in polygonal case, otherwise it can result in a polygon
// with a degenerate min=max range.
ExponentialLevelRangeIterator levels(
expBase, -std::numeric_limits<double>::infinity());
SegmentMerger<RingAppender, ExponentialLevelRangeIterator>
writer(appender, levels, /* polygonize */ true);
ContourGeneratorFromRaster<decltype(writer),
Expand All @@ -732,7 +743,12 @@ CPLErr GDALContourGenerateEx(GDALRasterBandH hBand, void *hLayer,
}
else
{
IntervalLevelRangeIterator levels(contourBase, contourInterval);
// Do not provide the actual minimum value to level iterator
// in polygonal case, otherwise it can result in a polygon
// with a degenerate min=max range.
IntervalLevelRangeIterator levels(
contourBase, contourInterval,
-std::numeric_limits<double>::infinity());
SegmentMerger<RingAppender, IntervalLevelRangeIterator> writer(
appender, levels, /* polygonize */ true);
ContourGeneratorFromRaster<decltype(writer),
Expand All @@ -746,8 +762,8 @@ CPLErr GDALContourGenerateEx(GDALRasterBandH hBand, void *hLayer,
GDALRingAppender appender(OGRContourWriter, &oCWI);
if (!fixedLevels.empty())
{
FixedLevelRangeIterator levels(&fixedLevels[0],
fixedLevels.size());
FixedLevelRangeIterator levels(
&fixedLevels[0], fixedLevels.size(), dfMinimum, dfMaximum);
SegmentMerger<GDALRingAppender, FixedLevelRangeIterator> writer(
appender, levels, /* polygonize */ false);
ContourGeneratorFromRaster<decltype(writer),
Expand All @@ -757,7 +773,7 @@ CPLErr GDALContourGenerateEx(GDALRasterBandH hBand, void *hLayer,
}
else if (expBase > 0.0)
{
ExponentialLevelRangeIterator levels(expBase);
ExponentialLevelRangeIterator levels(expBase, dfMinimum);
SegmentMerger<GDALRingAppender, ExponentialLevelRangeIterator>
writer(appender, levels, /* polygonize */ false);
ContourGeneratorFromRaster<decltype(writer),
Expand All @@ -767,7 +783,8 @@ CPLErr GDALContourGenerateEx(GDALRasterBandH hBand, void *hLayer,
}
else
{
IntervalLevelRangeIterator levels(contourBase, contourInterval);
IntervalLevelRangeIterator levels(contourBase, contourInterval,
dfMinimum);
SegmentMerger<GDALRingAppender, IntervalLevelRangeIterator>
writer(appender, levels, /* polygonize */ false);
ContourGeneratorFromRaster<decltype(writer),
Expand Down Expand Up @@ -804,7 +821,9 @@ struct ContourGeneratorOpaque
double dfNoDataValue, double dfContourInterval,
double dfContourBase, GDALContourWriter pfnWriter,
void *pCBData)
: levels(dfContourBase, dfContourInterval), writer(pfnWriter, pCBData),
: levels(dfContourBase, dfContourInterval,
-std::numeric_limits<double>::infinity()),
writer(pfnWriter, pCBData),
merger(writer, levels, /* polygonize */ false),
contourGenerator(nWidth, nHeight, bNoDataSet != 0, dfNoDataValue,
merger, levels)
Expand Down
49 changes: 35 additions & 14 deletions alg/marching_squares/level_generator.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -96,9 +96,10 @@ class FixedLevelRangeIterator
public:
typedef RangeIterator<FixedLevelRangeIterator> Iterator;

FixedLevelRangeIterator(const double *levels, size_t count,
double maxLevel = Inf)
: levels_(levels), count_(count), maxLevel_(maxLevel)
FixedLevelRangeIterator(const double *levels, size_t count, double minLevel,
double maxLevel)
: levels_(levels), count_(count), minLevel_(minLevel),
maxLevel_(maxLevel)
{
}

Expand All @@ -107,13 +108,15 @@ class FixedLevelRangeIterator
if (min > max)
std::swap(min, max);
size_t b = 0;
for (; b != count_ && levels_[b] < fudge(levels_[b], min); b++)
for (; b != count_ && levels_[b] < fudge(min, minLevel_, levels_[b]);
b++)
;
if (min == max)
return Range<Iterator>(Iterator(*this, int(b)),
Iterator(*this, int(b)));
size_t e = b;
for (; e != count_ && levels_[e] <= fudge(levels_[e], max); e++)
for (; e != count_ && levels_[e] <= fudge(max, minLevel_, levels_[e]);
e++)
;
return Range<Iterator>(Iterator(*this, int(b)),
Iterator(*this, int(e)));
Expand All @@ -126,10 +129,16 @@ class FixedLevelRangeIterator
return levels_[size_t(idx)];
}

double minLevel() const
{
return minLevel_;
}

private:
const double *levels_;
size_t count_;
double maxLevel_;
const double minLevel_;
const double maxLevel_;
};

struct TooManyLevelsException : public std::exception
Expand All @@ -150,8 +159,8 @@ struct IntervalLevelRangeIterator
typedef RangeIterator<IntervalLevelRangeIterator> Iterator;

// Construction by a offset and an interval
IntervalLevelRangeIterator(double offset, double interval)
: offset_(offset), interval_(interval)
IntervalLevelRangeIterator(double offset, double interval, double minLevel)
: offset_(offset), interval_(interval), minLevel_(minLevel)
{
}

Expand All @@ -165,7 +174,7 @@ struct IntervalLevelRangeIterator
if (!(df_i1 >= INT_MIN && df_i1 < INT_MAX))
throw TooManyLevelsException();
int i1 = static_cast<int>(df_i1);
double l1 = fudge(level(i1), min);
double l1 = fudge(min, minLevel_, level(i1));
if (l1 > min)
{
df_i1 = ceil((l1 - offset_) / interval_);
Expand All @@ -183,7 +192,7 @@ struct IntervalLevelRangeIterator
if (!(df_i2 >= INT_MIN && df_i2 < INT_MAX))
throw TooManyLevelsException();
int i2 = static_cast<int>(df_i2);
double l2 = fudge(level(i2), max);
double l2 = fudge(max, minLevel_, level(i2));
if (l2 > max)
{
df_i2 = floor((l2 - offset_) / interval_) + 1;
Expand All @@ -206,18 +215,24 @@ struct IntervalLevelRangeIterator
return idx * interval_ + offset_;
}

double minLevel() const
{
return minLevel_;
}

private:
const double offset_;
const double interval_;
const double minLevel_;
};

class ExponentialLevelRangeIterator
{
public:
typedef RangeIterator<ExponentialLevelRangeIterator> Iterator;

ExponentialLevelRangeIterator(double base)
: base_(base), base_ln_(std::log(base_))
ExponentialLevelRangeIterator(double base, double minLevel)
: base_(base), base_ln_(std::log(base_)), minLevel_(minLevel)
{
}

Expand All @@ -234,7 +249,7 @@ class ExponentialLevelRangeIterator
std::swap(min, max);

int i1 = index1(min);
double l1 = fudge(level(i1), min);
double l1 = fudge(min, minLevel_, level(i1));
if (l1 > min)
i1 = index1(l1);
Iterator b(*this, i1);
Expand All @@ -243,7 +258,7 @@ class ExponentialLevelRangeIterator
return Range<Iterator>(b, b);

int i2 = index2(max);
double l2 = fudge(level(i2), max);
double l2 = fudge(max, minLevel_, level(i2));
if (l2 > max)
i2 = index2(l2);
Iterator e(*this, i2);
Expand All @@ -256,6 +271,11 @@ class ExponentialLevelRangeIterator
return Range<Iterator>(b, e);
}

double minLevel() const
{
return minLevel_;
}

private:
int index1(double plevel) const
{
Expand All @@ -280,6 +300,7 @@ class ExponentialLevelRangeIterator
// exponentiation base
const double base_;
const double base_ln_;
const double minLevel_;
};

} // namespace marching_squares
Expand Down
Loading
Loading