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.

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

WIP: gdal: Remove special logic for finding overviews. #3975

Open
wants to merge 10 commits into
base: master
Choose a base branch
from

Conversation

cpaulik
Copy link
Contributor

@cpaulik cpaulik commented Aug 23, 2018

Trying to fix #3966

Current status:

We now send the correct image extent to gdal when reading so that resampling and overview selection is done by gdal.

This solves my original problem but does not implement #3966 (comment)

So if #3966 (comment) should also be implemented here I will need some guidance of where to get information about the resampling method and the warnings we should emit.

//calculate actual box2d of returned raster
view_transform t2(current_width, current_height, raster_extent_, 0, 0);
view_transform t2(im_width, im_height, raster_extent_, 0, 0);
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think that another transform view_transform t2 is not needed in our case. I would take inspiration from former implementation, where overviews were working the way we want now:

//calculate actual box2d of returned raster
box2d<double> feature_raster_extent(x_off, y_off, x_off + width, y_off + height);
intersect = t.backward(feature_raster_extent);

@talaj
Copy link
Member

talaj commented Aug 28, 2018

Have you been looking on visual test results? Many tests are broken.

So if #3966 (comment) should also be implemented here I will need some guidance of where to get information about the resampling method and the warnings we should emit.

I would implement it separately. Let's focus on basic functionality first.

@cpaulik
Copy link
Contributor Author

cpaulik commented Aug 28, 2018

I would implement it separately. Let's focus on basic functionality first.

I looked at the visual tests and to me it seemed that a lot of them show differences because of resampling methods. So I would like to avoid a situation where we fix a lot of visual tests now just to revert them in a second step.
I will look at it in detail again next week. On vacation this week and not at the computer a lot.

@talaj
Copy link
Member

talaj commented Aug 28, 2018

I looked at the visual tests and to me it seemed that a lot of them show differences because of resampling methods.

There are many tests with minor visual differences due to resampling and missing filter_factor, but there are also major differences with raster completely missing or in incorrect dimensions. We must ensure that all visual tests look visually good before we can merge this pull request and continue with subsequent steps.

I will look at it in detail again next week. On vacation this week and not at the computer a lot.

Of course, enjoy your vacation.

@cpaulik
Copy link
Contributor Author

cpaulik commented Sep 3, 2018

I think I fixed the incorrect dimensions that were returned in the last commit.

I also noticed that I am not running all visual tests locally because I get errors like

test/data-visual/styles/raster-color-to-alpha4b.xml ERROR (failed to initialize projection with: '+init=epsg:32630' in Map of 'test/data-visual/styles/raster-color-to-alpha4b.xml')

So I have to fix that in my local setup.

@talaj
Copy link
Member

talaj commented Sep 3, 2018

I think I fixed the incorrect dimensions that were returned in the last commit.

Great, I will try to take a look soon.

I also noticed that I am not running all visual tests locally because I get errors

If you used source bootstrap.sh previously, just call source mapnik-settings.env to restore some environment variables useful for testing.

@talaj
Copy link
Member

talaj commented Sep 4, 2018

Visual tests look good now.

I had to add forgotten filter_factor . Can you, please, add commit mapycz@b4a8ef0 to this pull request?

filter_factor is an option of RasterSymbolizer to allow to fetch raster in higher resolution for better quality output. Its default value is 1.0 for scaling="near" and 2.0 for other resampling filters.

Visual test gdal-max-image-area.xml doesn't make sense now as max_image_area option has been removed in this pull request. We should remove this test and remaining code from gdal_datasource.cpp:

// Maximum memory limitation for image will be simply based on the maximum
// area we allow for an image. The true memory footprint therefore will vary based
// on the type of imagery that exists. This is not the maximum size of an image
// on disk but rather the maximum size we will load into mapnik from GDAL.
// max_im_area based on 50 mb limit for RGBA
max_image_area_ = *params.get<mapnik::value_integer>("max_image_area", (50*1024*1024) / 4);

@flippmoke
Copy link
Member

filter_factor is an option of RasterSymbolizer to allow to fetch raster in higher resolution for better quality output.

I had thought that filter_factor did almost nothing now @talaj. I believe it only comes into use for select image rescaling algorithms.

https://github.com/mapnik/mapnik/blob/master/include/mapnik/image_scaling_traits.hpp#L200-L205

I know that gdal plugin has used it for overviews, though I can't recall why this is done.

https://github.com/mapnik/mapnik/blob/master/plugins/input/gdal/gdal_featureset.cpp#L244-L249

I think we should try to keep filter_factor as only a method to change resample algorithms.

@talaj
Copy link
Member

talaj commented Sep 4, 2018

@flippmoke See attribute_collector.hpp#L178-L197. Filter factor is being set to 2.0 for scaling != near.

It really make some sense for output quality, compare results from gdal-filter-factor-1-bilinear.xml and gdal-filter-factor-2-bilinear.xml

gdal-filter-factor-1-bilinear-600-400-1 0-agg-no-ff
gdal-filter-factor-1-bilinear-600-400-1 0-agg-ff

I know that at least in mapnik/python-mapnik#186 a user was complaining about blurry results until filter_factor was put back.

filter_factor originally comes from #563.

@cpaulik
Copy link
Contributor Author

cpaulik commented Sep 4, 2018

I've removed the max image area parameter and updated the changelog.

Also opened mapnik/test-data-visual#72 for the test data.

About the filter factor: If we want to do resampling in gdal then we should not need the filter factor?

Reading #563 it depends what a user would expect. If I provide mapnik with a gdal datasource with nearest neighbor overviews and I then instruct mapnik to use bilinear resampling then I would argue that this is at least a strange request to make, maybe user error.

I would also be interested to see how gdal-filter-factor-1-bilinear.xml would look like if resampling is done in gdal.

@flippmoke
Copy link
Member

I am against the way that filter_factor works currently. There needs to be two different variables for this process, rather then one. One changes the way that resampling algorithms (a few of them) operate, another changes the resolution of you expect GDAL to provide prior to any resampling that is done within the mapnik resampling (as far as I understand).

@talaj
Copy link
Member

talaj commented Sep 4, 2018

It's indeed strange logic. It seems to me as filter_factor was a hack to force resampling in Mapnik. Perhaps we can replace filter_factor in mapnik::query by resampling filter and pass it to RasterIO.

@cpaulik
Copy link
Contributor Author

cpaulik commented Sep 4, 2018

I just ran the tests with the resampling method hardcoded to bilinear in gdal.

GDALRasterIOExtraArg psExtraArg;
INIT_RASTERIO_EXTRA_ARG(psExtraArg);
psExtraArg.eResampleAlg = GRIORA_Bilinear;

This is the output of gdal-filter-factor-1-bilinear-600-400-1 0-agg
gdal-filter-factor-1-bilinear-600-400-1 0-agg

and this is the output of gdal-filter-factor-1-nearest-600-400-1 0-agg
gdal-filter-factor-1-nearest-600-400-1 0-agg

Those look very close to the filter factor 2 billinear when the resampling is done in mapnik instead of gdal.

We will need a lookup table between the mapnik resampling methods and https://www.gdal.org/gdal_8h.html#a640ada511cbddeefac67c548e009d5ac so we can pass it to RasterIO. For that we need the resampling method as a parameter in mapnik::query I guess.

@talaj
Copy link
Member

talaj commented Sep 5, 2018

We will need a lookup table between the mapnik resampling methods and https://www.gdal.org/gdal_8h.html#a640ada511cbddeefac67c548e009d5ac so we can pass it to RasterIO. For that we need the resampling method as a parameter in mapnik::query I guess.

I also see it as the way to go.

@cpaulik
Copy link
Contributor Author

cpaulik commented Sep 5, 2018

So I guess the way forward for this MR is the following:

  • add resampling method as a parameter to mapnik::query
  • Make lookup table between gdal and mapnik resampling methods and display warning if no exact match is found
  • Tell mapnik that no further resampling inside of mapnik is necessary because gdal has already done it?

Do you think we need the last point?

I also could use some pointers on point 1.

What is a good way in mapnik/C++ to implement a lookup table? I can make a series of if else statements but don't know if there isn't a more elegant way?

@talaj
Copy link
Member

talaj commented Sep 5, 2018

@cpaulik

Do you think we need the last point?

Mapnik will not do any further resampling if returned raster is of correct size. So I think there is no need to pass such information from gdal datasource. We can make a unit test that ensures such behavior.

What is a good way in mapnik/C++ to implement a lookup table? I can make a series of if else statements but don't know if there isn't a more elegant way?

Simple switch statement should be appropriate.

I also could use some pointers on point 1.

I would do it the same way as filter_factor is set to mapnik::query. The actual setting is done here. There is a device called attribute_collector which collects all attribute names (such as [name]) used in symbolizers and rule filters and also collects value of filter_factor.

@cpaulik
Copy link
Contributor Author

cpaulik commented Oct 19, 2018

I'm starting to work on this again.

Is there a way for me to run a single test so that I can debug through the code for understanding it better? How do you normally do that when developing in mapnik?

@talaj
Copy link
Member

talaj commented Oct 21, 2018

@cpaulik Yes, both visual and unit tests can run single test case.

Visual tests

This is an example how to run single visual test:

$ test/visual/run --agg -s 1 gdal-overview
.
Visual rendering: 0 failed / 1 passed / 0 overwritten / 0 errors

gdal-overview will run gdal-overview.xml by looking to the path where visual tests are stored. You can also run visual test by path to xml file:

$ test/visual/run ~/some/path/test.xml

For more information about visual test runner options run test/visual/run -h.

Unit tests

It works similarly with unit tests:

$ test/unit/run "gdal"
===============================================================================
All tests passed (6 assertions in 1 test case)

@cpaulik
Copy link
Contributor Author

cpaulik commented Oct 22, 2018

Thanks. I've added two more commits that add the resampling method to the query object and implement a lookup between gdal and mapnik resampling methods.

The first lookup table is based on a very quick reading of what the mapnik resampling methods do so it is probably not yet the best.

If that looks good, I can rebase onto the latest master and we can start talking about which tests to add/change/remove.

@talaj
Copy link
Member

talaj commented Oct 28, 2018

Looks good. I'm thinking what to do with 7aa06d9. @flippmoke is right that filter-factor does two different things at once. I would just remove that commit with no replacement.

We can then remove these two visual tests:

  • gdal-max-image-area.xml
  • gdal-overview-bilinear-scaling-filter-factor.xml

These tests need a bit of attention, especially the first one:

  • tiff-nodata-edge-rgba.xml
  • tiff-resampling.xml
  • tiff-nodata-tolerance.xml

@cpaulik cpaulik force-pushed the gdal-do-not-use-special-mapnik-logic-for-overviews branch from 61ff94f to aa8fc60 Compare July 8, 2019 14:16
@cpaulik
Copy link
Contributor Author

cpaulik commented Jul 8, 2019

I've started working on this again. Sorry for letting this sit so long.

The two tests

  • gdal-max-image-area.xml
  • gdal-overview-bilinear-scaling-filter-factor.xml

have been removed.

I will take a closer look at the three other tests you highlighted. But there are 220 visual tests that fail in total. Did you already look at those and select only the three for closer inspection or should I look at all of them?

@cpaulik
Copy link
Contributor Author

cpaulik commented Jul 9, 2019

With the last commit I fixed all the tests other than the ones related to overviews. I will double check but I think it is expected that they fail now.

@cpaulik cpaulik force-pushed the gdal-do-not-use-special-mapnik-logic-for-overviews branch from ecef7c5 to 21b2666 Compare July 9, 2019 16:20
@cpaulik
Copy link
Contributor Author

cpaulik commented Jul 9, 2019

I think It should be ok now. 128 tests fail but the results look reasonably close and are related to resampling differences as far as I could see.

@cpaulik cpaulik force-pushed the gdal-do-not-use-special-mapnik-logic-for-overviews branch from 21b2666 to ab56586 Compare November 1, 2019 09:01
@cpaulik
Copy link
Contributor Author

cpaulik commented Nov 1, 2019

@lightmare Sorry for tagging you here but you seem to be the most active committer in recent months.

Are you guys still interested in this fix?

I think it is ready but I would like to get the OK from a maintainer before I replace all the test images with slight differences.

@lightmare
Copy link
Contributor

@cpaulik I might not be qualified to weigh in on this. That comment of mine you mentioned above shows I wasn't aware of the intricacies of what I was talking about. Apparently I thought it was simply a matter of mapping Mapnik's scaling enum to GDAL's scaling enum, without realizing that Mapnik's scaling is buried in RasterSymbolizer, whereas GDAL scaling happens before RasterSymbolizer gets into play. But since you asked, I dived into the code to educate myself.

TL;DR:
I'm in favour of removing intricate/flawed logic from gdal plugin, passing as much as possible to libgdal.
I'm not a fan of pumping parameters up from symbolizers to datasource and down to features.

I see you followed the example of filter_factor handling, which I don't like, but understand it was the right way to do this. I was able to come up with two alternative solutions, one clearly inferior from user perspective (adding GDAL scaling parameter to datasource), the other more complicated (lazy loading of images, from RasterSymbolizer). So in the context of existing code, I think this PR is okay.


Since this brought attention to attribute_collector being abused to pump symbolizer attributes into mapnik::query, I'd like to see this hack removed some day.

filter_factor is a parameter whose meaning is defined by the scaling method, but its use in mapnik datasources is dubious. For whatever reason mapnik::raster objects store filter_factor but not scaling_method. At some point filter_factor was also used to expand the source extent, because scaling usually requires some extra pixels around the area to create seamless tiles. But to do this properly you need to know which scaling method you're going to use, or more precisely, how exactly that method's implementation interprets filter_factor (see e.g. in agg lanczos it's always >= 2, whereas in gdal lanczos it's always == 3 (but in this case gdal handles it, so mapnik shouldn't care)).

@cpaulik
Copy link
Contributor Author

cpaulik commented Jan 24, 2020

Thanks. I would like to get this merged.

I realize that I should update the visual test data. But somehow I've not yet managed to run the visual tests locally so that all of them pass. Could one of you maybe run the tests and generate the replacement images?

@cpaulik cpaulik force-pushed the gdal-do-not-use-special-mapnik-logic-for-overviews branch from ab56586 to a47c749 Compare January 24, 2020 14:09
@cpaulik
Copy link
Contributor Author

cpaulik commented Feb 11, 2020

@artemp Would you mind giving your opinion on this one also? It is ready from my side but I do struggle with re-creating all the test images in a consistent way in my development environment.

@artemp
Copy link
Member

artemp commented Feb 11, 2020

@cpaulik - looks good 👍 Lets see if I can manage to update test images :) My plan is:

  • get a fresh mapnik clone
  • ./bootstrap.sh
  • ./configure && make JOBS=4
  • make test JOBS=4 <-- all tests are green

@artemp
Copy link
Member

artemp commented Feb 11, 2020

reference

gdal-overview-gray8-256-256-1 0-agg-reference

gdal-do-not-use-special-mapnik-logic-for-overviews

gdal-overview-gray8-256-256-1 0-agg

I have 140 visual tests failing, most look OK-ish but I see some artifacts like second letter V above ^ Looks like regression or perhaps inevitable result of gdal implementation ?? @cpaulik

/cc @talaj @lightmare

@cpaulik
Copy link
Contributor Author

cpaulik commented Feb 11, 2020

If the test uses nearest-neighbor resampling and the overview does not fit exactly then I think this is what should be expected.

@cpaulik
Copy link
Contributor Author

cpaulik commented Feb 11, 2020

The psExtraArg is now missing in the GDT_Int32 case. Otherwise it looks good.

@lightmare
Copy link
Contributor

I'm getting the same broken "OVERVIEW" image as @artemp above, and I don't think nearest-neighbour resampling explains that satisfactorily. Try flipping between the two images, you'll notice the broken one is shifted up. Next, check the letter 'E' - the middle horizontal line got thicker, while the gap above it shrank. When you resample with nearest-neighbour, you either duplicate some pixels (larger output), or drop some pixels (smaller output), but you never both duplicate and drop pixels in one dimension. Or so I thought.

More importantly, though, resampling shouldn't be happening here. The test (overview-test-gray8) asks for a 2x2 down-scaled image, and there's an overview matching that scale exactly.

@cpaulik
Copy link
Contributor Author

cpaulik commented Feb 12, 2020

Thanks @lightmare I'll try to figure out what is happening.

@cpaulik
Copy link
Contributor Author

cpaulik commented Feb 12, 2020

I get the following log output:

✘Mapnik LOG> 2020-02-12 10:37:24: cairo_renderer: Scale=2
Mapnik LOG> 2020-02-12 10:37:24: cairo_renderer: Start map processing bbox=box2d(0.0000000000000000,0.0000000000000000,512.0000000000000000,512.0000000000000000)
Mapnik LOG> 2020-02-12 10:37:24: cairo_renderer: Start processing layer=dataraster
Mapnik LOG> 2020-02-12 10:37:24: cairo_renderer: -- datasource=0x5570c1447460
Mapnik LOG> 2020-02-12 10:37:24: cairo_renderer: -- query_extent=box2d(0.0000000000000000,0.0000000000000000,512.0000000000000000,512.0000000000000000)
Mapnik LOG> 2020-02-12 10:37:24: cairo_renderer:start style processing
Mapnik LOG> 2020-02-12 10:37:24: gdal_featureset: Next feature in Dataset=0x5570c144afd0
Mapnik LOG> 2020-02-12 10:37:24: gdal_featureset: margin_x=2
Mapnik LOG> 2020-02-12 10:37:24: gdal_featureset: margin_y=2
Mapnik LOG> 2020-02-12 10:37:24: gdal_featureset: x_off=-2
Mapnik LOG> 2020-02-12 10:37:24: gdal_featureset: y_off=3582
Mapnik LOG> 2020-02-12 10:37:24: gdal_featureset: end_x=514
Mapnik LOG> 2020-02-12 10:37:24: gdal_featureset: end_y=4098
Mapnik LOG> 2020-02-12 10:37:24: gdal_featureset: Feature Raster extent=box2d(0.0000000000000000,3582.0000000000000000,514.0000000000000000,4096.0000000000000000)
Mapnik LOG> 2020-02-12 10:37:24: gdal_featureset: Raster extent=box2d(0.0000000000000000,0.0000000000000000,4096.0000000000000000,4096.0000000000000000)
Mapnik LOG> 2020-02-12 10:37:24: gdal_featureset: Feature Raster extent=box2d(0.0000000000000000,0.0000000000000000,514.0000000000000000,514.0000000000000000)
Mapnik LOG> 2020-02-12 10:37:24: gdal_featureset: View extent=box2d(0.0000000000000000,0.0000000000000000,512.0000000000000000,512.0000000000000000)
Mapnik LOG> 2020-02-12 10:37:24: gdal_featureset: Query resolution=0.5,0.5
Mapnik LOG> 2020-02-12 10:37:24: gdal_featureset: StartX=0 StartY=3582 Width=514 Height=514
Mapnik LOG> 2020-02-12 10:37:24: gdal_featureset: IM StartX=0 StartY=3582 Width=256 Height=256
Mapnik LOG> 2020-02-12 10:37:24: gdal_featureset: Requested Size from gdal=(514,514)
Mapnik LOG> 2020-02-12 10:37:24: gdal_featureset: Internal Image Size=(256,256)
Mapnik LOG> 2020-02-12 10:37:24: gdal_featureset: Reading band=1
Mapnik LOG> 2020-02-12 10:37:24: gdal_featureset: requested resampling method= near
Mapnik LOG> 2020-02-12 10:37:24: gdal_featureset: Using GDAL Nearest Neighbor resampling
Mapnik LOG> 2020-02-12 10:37:24: cairo_renderer:end style processing
Mapnik LOG> 2020-02-12 10:37:24: cairo_renderer: End layer processing
Mapnik LOG> 2020-02-12 10:37:24: cairo_renderer: End map processing
Mapnik LOG> 2020-02-12 10:37:24: stroker: Destroy stroker=0x5570c14958b0
Mapnik LOG> 2020-02-12 10:37:24: png_reader: bit_depth=8,color_type=6

It seems that a 514x514 image is requested from gdal which should be 512x512.

I tried hardcoding this to 512x512 and the test passes.

So the question is where do the margin and offset calculations come from and can I change them?

@cpaulik cpaulik force-pushed the gdal-do-not-use-special-mapnik-logic-for-overviews branch from a47c749 to 90815d3 Compare February 12, 2020 11:26
@cpaulik cpaulik force-pushed the gdal-do-not-use-special-mapnik-logic-for-overviews branch from 90815d3 to 9d1d095 Compare February 12, 2020 11:35
@cpaulik
Copy link
Contributor Author

cpaulik commented Feb 14, 2020

I've rebased this to include the change for Int32 images. Do you want to change the margin calculations? If so in this MR or in a separate one?

@artemp
Copy link
Member

artemp commented Feb 14, 2020

@cpaulik - yes, I think margin calculations fix should be in this PR. I didn't have a chance to look into this issue but lets track progress here.

@lightmare
Copy link
Contributor

I've taken a look at the issue with margin. First I tried setting it to zero just to see what happens --> tiff-reprojection-* tests broke, so margin is needed. Next I dived into transforms and rounding.

Let's annotate some local variables:
intersect = query bbox in source srs coordinates
box = query bbox in source pixel coordinates
feature_raster_extent = box + margin in source srs coordinates

There's an issue with intermediate image (that is the image we'll get from GDAL) dimensions im_width and im_height. These are calculated from intersect bbox, but the intermediate image will hold pixels corresponding to feature_raster_extent, i.e. including margin.

Then we can mitigate some rounding errors, by requesting floating-point bbox from GDAL: bFloatingPointWindowValidity

Indicate if dfXOff, dfYOff, dfXSize and dfYSize are set. Mostly reserved from the VRT driver to communicate a more precise source window. Must be such that dfXOff - nXOff < 1.0 and dfYOff - nYOff < 1.0 and nXSize - dfXSize < 1.0 and nYSize - dfYSize < 1.0

They say mostly reserved for VRT driver, but GTiff driver supports it as well. I'm not quite sure about the constraints, though. Assuming they also imply dfXOff >= nXoff && dfYOff >= nYOff, offset constraints mean that the sub-pixel offset must be within the top-most pixel of the integer bbox -- trivial. But the size constraint confuses me. I think that the integer bbox should cover the whole request, so that the driver can tell how many pixels it'll need to load just from looking at integer dimensions, but then the constraint would need to be nXSize - dfXSize < 2.0. I looked into GDAL and didn't find these constraints enforced, so I'm using this relaxed interpretation.

Here's my patch enabling sub-pixel GDAL RasterIO:

commit cbce9ba5e3a8b0c82f3e8e5e143ef9c2745c4dfb
Author: Mickey Rose <mickey.rose@seznam.cz>
Date:   Wed Mar 4 21:28:24 2020 +0100

    gdal: request floating-point bbox

diff --git a/plugins/input/gdal/gdal_featureset.cpp b/plugins/input/gdal/gdal_featureset.cpp
index 04b5966ed..c29de4534 100644
--- a/plugins/input/gdal/gdal_featureset.cpp
+++ b/plugins/input/gdal/gdal_featureset.cpp
@@ -166,6 +166,9 @@ feature_ptr gdal_featureset::get_feature(mapnik::query const& q)
     GDALRasterBand * grey = 0;
     CPLErr raster_io_error = CE_None;
 
+    GDALRasterIOExtraArg psExtraArg;
+    INIT_RASTERIO_EXTRA_ARG(psExtraArg);
+
     /*
 #ifdef MAPNIK_LOG
       double tr[6];
@@ -182,47 +185,40 @@ feature_ptr gdal_featureset::get_feature(mapnik::query const& q)
     box2d<double> box = t.forward(intersect);
     double filter_factor = q.get_filter_factor();
 
+    // resolution is 1 / requested pixel size
+    const double width_res = std::get<0>(q.resolution());
+    const double height_res = std::get<1>(q.resolution());
+
     //size of resized output pixel in source image domain
-    double margin_x = 1.0 / (std::fabs(dx_) * std::get<0>(q.resolution()));
-    double margin_y = 1.0 / (std::fabs(dy_) * std::get<1>(q.resolution()));
+    double inv_margin_x = std::fabs(dx_) * width_res;
+    double inv_margin_y = std::fabs(dy_) * height_res;
+    double margin_x = inv_margin_x >= 1 ? 1 : 1 / inv_margin_x;
+    double margin_y = inv_margin_y >= 1 ? 1 : 1 / inv_margin_y;
     MAPNIK_LOG_DEBUG(gdal) << "gdal_featureset: margin_x=" << margin_x;
     MAPNIK_LOG_DEBUG(gdal) << "gdal_featureset: margin_y=" << margin_y;
-    if (margin_x < 1)
-    {
-        margin_x = 1.0;
-    }
-    if (margin_y < 1)
-    {
-        margin_y = 1.0;
-    }
+
+    //clip to available data
+    double dfXEnd = std::min(box.maxx() + margin_x, static_cast<double>(raster_width_));
+    double dfYEnd = std::min(box.maxy() + margin_y, static_cast<double>(raster_height_));
+    psExtraArg.dfXOff = std::max(box.minx() - margin_x, 0.0);
+    psExtraArg.dfYOff = std::max(box.miny() - margin_y, 0.0);
+    psExtraArg.dfXSize = dfXEnd - psExtraArg.dfXOff;
+    psExtraArg.dfYSize = dfYEnd - psExtraArg.dfYOff;
+    psExtraArg.bFloatingPointWindowValidity = TRUE;
+    box.init(psExtraArg.dfXOff, psExtraArg.dfYOff, dfXEnd, dfYEnd);
 
     //select minimum raster containing whole box
-    int x_off = static_cast<int>(std::floor((box.minx() - margin_x) + .5));
-    int y_off = static_cast<int>(std::floor((box.miny() - margin_y) + .5));
-    int end_x = static_cast<int>(std::floor((box.maxx() + margin_x) + .5));
-    int end_y = static_cast<int>(std::floor((box.maxy() + margin_y) + .5));
+    int x_off = static_cast<int>(std::floor(psExtraArg.dfXOff));
+    int y_off = static_cast<int>(std::floor(psExtraArg.dfYOff));
+    int end_x = static_cast<int>(std::ceil(dfXEnd));
+    int end_y = static_cast<int>(std::ceil(dfYEnd));
     MAPNIK_LOG_DEBUG(gdal) << "gdal_featureset: x_off=" << x_off;
     MAPNIK_LOG_DEBUG(gdal) << "gdal_featureset: y_off=" << y_off;
     MAPNIK_LOG_DEBUG(gdal) << "gdal_featureset: end_x=" << end_x;
     MAPNIK_LOG_DEBUG(gdal) << "gdal_featureset: end_y=" << end_y;
 
-    //clip to available data
-    if (x_off < 0)
-    {
-        x_off = 0;
-    }
-    if (y_off < 0)
-    {
-        y_off = 0;
-    }
-    if (end_x > (int)raster_width_)
-    {
-        end_x = raster_width_;
-    }
-    if (end_y > (int)raster_height_)
-    {
-        end_y = raster_height_;
-    }
+    //calculate actual box2d of returned raster
+    box2d<double> feature_raster_extent = t.backward(box);
 
     // width and height of the portion of the source image we are requesting
     int width = end_x - x_off;
@@ -230,12 +226,8 @@ feature_ptr gdal_featureset::get_feature(mapnik::query const& q)
 
     double im_offset_x = x_off;
     double im_offset_y = y_off;
-
-    // resolution is 1 / requested pixel size
-    const double width_res = std::get<0>(q.resolution());
-    const double height_res = std::get<1>(q.resolution());
-    int im_width = int(width_res * intersect.width() + 0.5);
-    int im_height = int(height_res * intersect.height() + 0.5);
+    int im_width = static_cast<int>(std::round(width_res * feature_raster_extent.width()));
+    int im_height = static_cast<int>(std::round(height_res * feature_raster_extent.height()));
 
     // no upsampling in gdal
     if (im_width > width)
@@ -243,12 +235,6 @@ feature_ptr gdal_featureset::get_feature(mapnik::query const& q)
     if (im_height > height)
       {im_height=height;}
 
-    //calculate actual box2d of returned raster
-    box2d<double> feature_raster_extent(x_off, y_off, x_off + width, y_off + height);
-
-    MAPNIK_LOG_DEBUG(gdal) << "gdal_featureset: Feature Raster extent=" << feature_raster_extent;
-    feature_raster_extent = t.backward(feature_raster_extent);
-
     MAPNIK_LOG_DEBUG(gdal) << "gdal_featureset: Raster extent=" << raster_extent_;
     MAPNIK_LOG_DEBUG(gdal) << "gdal_featureset: Feature Raster extent=" << feature_raster_extent;
     MAPNIK_LOG_DEBUG(gdal) << "gdal_featureset: View extent=" << intersect;
@@ -263,8 +249,6 @@ feature_ptr gdal_featureset::get_feature(mapnik::query const& q)
         MAPNIK_LOG_DEBUG(gdal) << "gdal_featureset: Internal Image Size=(" << im_width << "," << im_height << ")";
         MAPNIK_LOG_DEBUG(gdal) << "gdal_featureset: Reading band=" << band_;
         MAPNIK_LOG_DEBUG(gdal) << "gdal_featureset: requested resampling method=" << scaling_method_to_string(q.get_scaling_method());
-        GDALRasterIOExtraArg psExtraArg;
-        INIT_RASTERIO_EXTRA_ARG(psExtraArg);
         psExtraArg.eResampleAlg = get_gdal_resample_alg(q.get_scaling_method());
         if (band_ > 0) // we are querying a single band
         {

@lightmare
Copy link
Contributor

Results of the above patch are promising, but far from perfect. I'm still getting lots of visual test failures, categorized by severity here (checked boxes mean I consider that not-a-failure):

  • NO APPARENT DIFFERENCES:

    • gdal-filter-factor-1-bilinear FAILED (112895 different pixels)
    • gdal-filter-factor-1-lanczos FAILED (118673 different pixels)
    • gdal-filter-factor-2-bilinear FAILED (112895 different pixels)
    • gdal-filter-factor-2-lanczos FAILED (118673 different pixels)
    • gdal-filter-factor-3-lanczos FAILED (113950 different pixels)
    • gdal-filter-factor-default-bilinear FAILED (112895 different pixels)
    • gdal-filter-factor-default-lanczos FAILED (118673 different pixels)
    • layer-compositing-offset FAILED (13209 different pixels)
    • raster-color-to-alpha2 FAILED (13794 different pixels)
    • tiff-non-alpha-mask-band FAILED (2821 different pixels)
    • tiff-opaque-edge-gdal FAILED (3678 different pixels)
    • tiff-opaque-edge-gdal2 FAILED (8141 different pixels)
    • tiff-reprojection-1 FAILED (19140 different pixels)
  • FEW PIXEL DIFFERENCES:

    • gdal-filter-factor-1-nearest FAILED (291 different pixels)
    • gdal-filter-factor-2-nearest FAILED (291 different pixels)
    • gdal-filter-factor-default-nearest FAILED (291 different pixels)
    • raster-color-to-alpha4 FAILED (44 different pixels)
    • raster-color-to-alpha5 FAILED (53 different pixels)
    • raster_colorizer FAILED (53 different pixels)
    • raster_symbolizer FAILED (239 different pixels)
    • tiff_colortable FAILED (2429 different pixels)
    • tiff_colortable_custom_nodata FAILED (2429 different pixels)
    • vrt_colortable FAILED (293 different pixels)
  • APPARENT OUTLINE SHIFT:

    • raster-colorizer-bilinear-scaling-gdal FAILED (992 different pixels)
    • raster-scaling-nodata FAILED (5500 different pixels)
    • tiff-alpha-gdal-tiled FAILED (2004 different pixels)
    • tiff-edge-alignment-gdal1 FAILED (44195 different pixels)
    • tiff-edge-alignment-gdal1 FAILED (48848 different pixels)
    • tiff-edge-alignment-gdal2 FAILED (44360 different pixels)
    • tiff-edge-alignment-gdal2 FAILED (41495 different pixels)
  • APPARENT DISTORTION:

    • tiff-reprojection-2-tiled FAILED (288020 different pixels)
    • tiff-reprojection-2 FAILED (9029 different pixels)
    • tiff-reprojection-3 FAILED (4906 different pixels)
  • TOTAL WRECK:

    • tiff-nodata-tolerance FAILED (22781 different pixels)

@cpaulik
Copy link
Contributor Author

cpaulik commented May 4, 2021

What can we do to get this merged?

I don't know enough of the mapnik internals to make any decisions here or really help with the margin fixes.

I'm also not sure if this would have any impact on issues like

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

Mapnik does not use overviews in the same way gdal does
5 participants