Skip to content

Commit

Permalink
Merge pull request #10170 from rouault/fix_10167
Browse files Browse the repository at this point in the history 
Contour: include minimum raster value as contour line when it exactly matches the first level
  • Loading branch information
@rouault
rouault committed Jun 18, 2024
2 parents d39e40d + 68b400f commit 0decb44
Show file tree
Hide file tree
Showing 10 changed files with 291 additions and 137 deletions.
63 changes: 41 additions & 22 deletions alg/contour.cpp
View file
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
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

0 comments on commit 0decb44

Please sign in to comment.