-
Notifications
You must be signed in to change notification settings - Fork 825
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
base: master
Are you sure you want to change the base?
WIP: gdal: Remove special logic for finding overviews. #3975
Conversation
//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); |
There was a problem hiding this comment.
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:
mapnik/plugins/input/gdal/gdal_featureset.cpp
Lines 200 to 202 in 333ef9f
//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); |
Have you been looking on visual test results? Many tests are broken.
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. |
There are many tests with minor visual differences due to resampling and missing
Of course, enjoy your vacation. |
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
So I have to fix that in my local setup. |
Great, I will try to take a look soon.
If you used |
Visual tests look good now. I had to add forgotten
Visual test mapnik/plugins/input/gdal/gdal_datasource.cpp Lines 88 to 93 in 2929c4a
|
I had thought that 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 |
@flippmoke See It really make some sense for output quality, compare results from gdal-filter-factor-1-bilinear.xml and gdal-filter-factor-2-bilinear.xml I know that at least in mapnik/python-mapnik#186 a user was complaining about blurry results until
|
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 I would also be interested to see how gdal-filter-factor-1-bilinear.xml would look like if resampling is done in gdal. |
I am against the way that |
It's indeed strange logic. It seems to me as |
I just ran the tests with the resampling method hardcoded to bilinear in gdal.
This is the output of and this is the output of 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 |
I also see it as the way to go. |
So I guess the way forward for this MR is the following:
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 |
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.
Simple
I would do it the same way as |
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? |
@cpaulik Yes, both visual and unit tests can run single test case. Visual testsThis 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
$ test/visual/run ~/some/path/test.xml For more information about visual test runner options run Unit testsIt works similarly with unit tests: $ test/unit/run "gdal"
===============================================================================
All tests passed (6 assertions in 1 test case) |
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. |
Looks good. I'm thinking what to do with 7aa06d9. @flippmoke is right that We can then remove these two visual tests:
These tests need a bit of attention, especially the first one:
|
61ff94f
to
aa8fc60
Compare
I've started working on this again. Sorry for letting this sit so long. The two tests
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? |
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. |
ecef7c5
to
21b2666
Compare
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. |
21b2666
to
ab56586
Compare
@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. |
@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 see you followed the example of Since this brought attention to
|
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? |
ab56586
to
a47c749
Compare
@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. |
@cpaulik - looks good 👍 Lets see if I can manage to update test images :) My plan is:
|
referencegdal-do-not-use-special-mapnik-logic-for-overviewsI have 140 visual tests failing, most look OK-ish but I see some artifacts like second letter /cc @talaj @lightmare |
If the test uses nearest-neighbor resampling and the overview does not fit exactly then I think this is what should be expected. |
The |
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 ( |
Thanks @lightmare I'll try to figure out what is happening. |
I get the following log output:
It seems that a 514x514 image is requested from gdal which should be 512x512. I tried hardcoding this to So the question is where do the margin and offset calculations come from and can I change them? |
a47c749
to
90815d3
Compare
90815d3
to
9d1d095
Compare
I've rebased this to include the change for |
@cpaulik - yes, I think |
I've taken a look at the issue with margin. First I tried setting it to zero just to see what happens --> Let's annotate some local variables: There's an issue with intermediate image (that is the image we'll get from GDAL) dimensions Then we can mitigate some rounding errors, by requesting floating-point bbox from GDAL: bFloatingPointWindowValidity
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 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
{ |
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):
|
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 I'm also not sure if this would have any impact on issues like |
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.