[PATCH] Two Bugs in Contour Algorithm

[GregersP]

I have found two bugs; one that was always in the contour algorithm, and one which was introduced by my patch from ticket #3129. A further explanation is given below.

Bug 1 - erroneous directions

The bug is in the contour algorithms way to ensure correct direction (my earlier patch from ticket #3129). In my eager to simplify the if statement, I seem to have oversimplified the situation, such that the checks for "left is high side" were not done diagonally - this meant that contour segments which were exactly horizontal or vertical would have arbitrary direction.

This can be solved by expanding the check for "left is high side" to be done using diagonally opposite points.

Basis: When we check the left side for intersection, we know that the top left point will never be the resulting intersection (due to the nature of the comparisons made in the "Intersect" function). The same is true for the starting point of the other sections due to symmetry - i.e. the bottom left point for the bottom section, the bottom right point for the right section, and the top right point for the right section.

So: In the special case of the contour following an edge exactly, we can not use the former checks of two point on the same side, e.g. the top left corner and bottom left corner both being 10.0 meters does not tell us if "left is high side". But if we check diagonally instead (as I did before I oversimplified my patch) the check of bottom left point vs. top right point will give us the correct answer, when check the property for a "left to bottom" contour.

When it comes to a contour from left to right (possibly horizontal), we can use that the top left point will never be on said line (due to symmetry neither will the bottom right point). So we in this case we can still check using two diagonal points.

Bug 2 - duplicate contour segments

During my work on the problem above, I stumbled upon the other bug: When creating contour segments, some are correctly left out due to symmetry conditions. Though, when the contour exactly follows the border of the rectangle (pixel) currently being processed, it is produced for all/any edges. The result is that one will see extra (duplicate) contour segments on top of the "real" contour. This is due to the fact, that only one of the two duals (the "up segment" or the "down segment") will be used when connecting the segments to a single contour.

In the patch attached, I have arbitrarily chosen to exclude those following the left side and the top side.

[GregersP]

My patch

[rouault]

Assigning to Chaitanya as he's applied the first patch from #3129

[chaitanya]

I will create a test along with the patch.

Gregers, can you give me a sample raster that demonstrates these issues? Please send it to chaitanya(a)osgeo.in if you want to keep it confidential.

[GregersP]

Replying to chaitanya:

I will create a test along with the patch. Gregers, can you give me a sample raster that demonstrates these issues? Please send it to chaitanya(a)osgeo.in if you want to keep it confidential.

Hi Chaitanya,

You should by now have received my example by email.

[chaitanya]

Added the patch and a test in r19878 in trunk. Added the patch in r19879.

[rouault]

Chaitanya, you forgot to commit autotest/utilities/data/contour_orientation.tif :

TEST: test_gdal_contour_5 ... ERROR 4: `data/contour_orientation.tif' does not exist in the file system,
and is not recognised as a supported dataset name.

[chaitanya]

oops! Added the data file in r19880.

[rouault]

Chaitanya,

I observe new failures on 3 slavebots - szekerest-vc71-full, szekerest-vc80x64-full and szekerest-vc90x64-full - that seem to be directly related to the changes done in that ticket. See  http://buildbot.osgeo.org:8500/waterfall

 ------------ Failures ------------
Script: utilities/test_gdal_contour.py
  TEST: test_gdal_contour_2 ... fail
  TEST: test_gdal_contour_3 ... fail
Script: utilities/test_gdal_rasterize.py
  TEST: test_gdal_rasterize_3 ... fail
    did not get expected nodata value
  TEST: test_gdal_rasterize_4 ... fail (blowup)
 ----------------------------------

In fact, 'valgrind gdal_contour ../gdrivers/data/n43.dt0 tmp/n43dt0.shp -i 10 -3d' reveals that :

==22571== Conditional jump or move depends on uninitialised value(s)
==22571==    at 0x54E5B21: GDALContourGenerator::ProcessRect(double, double, double, double, double, double, double, double, double, double, double, double) (contour.cpp:604)
==22571==    by 0x54E6417: GDALContourGenerator::ProcessPixel(int) (contour.cpp:422)
==22571==    by 0x54E6783: GDALContourGenerator::FeedLine(double*) (contour.cpp:791)
==22571==    by 0x54E6BF5: GDALContourGenerate (contour.cpp:1600)
==22571==    by 0x401E7F: main (gdal_contour.cpp:266)
...10==22571== 
==22571== Conditional jump or move depends on uninitialised value(s)
==22571==    at 0x54E5B21: GDALContourGenerator::ProcessRect(double, double, double, double, double, double, double, double, double, double, double, double) (contour.cpp:604)
==22571==    by 0x54E6300: GDALContourGenerator::ProcessPixel(int) (contour.cpp:414)
==22571==    by 0x54E6783: GDALContourGenerator::FeedLine(double*) (contour.cpp:791)
==22571==    by 0x54E6BF5: GDALContourGenerate (contour.cpp:1600)
==22571==    by 0x401E7F: main (gdal_contour.cpp:266)
...20...30...40...50...60...70...80...90...==22571== 
==22571== Conditional jump or move depends on uninitialised value(s)
==22571==    at 0x54E5B21: GDALContourGenerator::ProcessRect(double, double, double, double, double, double, double, double, double, double, double, double) (contour.cpp:604)
==22571==    by 0x54E6527: GDALContourGenerator::ProcessPixel(int) (contour.cpp:430)
==22571==    by 0x54E6783: GDALContourGenerator::FeedLine(double*) (contour.cpp:791)
==22571==    by 0x54E67F4: GDALContourGenerator::FeedLine(double*) (contour.cpp:804)
==22571==    by 0x54E6BF5: GDALContourGenerate (contour.cpp:1600)
==22571==    by 0x401E7F: main (gdal_contour.cpp:266)

The reason is that the CPLErr eErr declared at line 531 is not initialized, and not necessarily set afterwards, which lead to random exits when reaching line 604. I guess it should be initialized to CE_None, right ? I let you check and fix as appropriate.

[chaitanya]

Replying to rouault:

Applied the suggested fix in r19881 in trunk and r19882 in 1.7 branch.

I'm waiting for the buildbots showing the error to come online to close this ticket.