-
Notifications
You must be signed in to change notification settings - Fork 36
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
Unexpected polygon rasterization results #376
Comments
Thanks for the issue and full MWE! Making issues is never reckless, please report anything you think is even slightly wrong, slow, or confusing! But also it's best to separate your two concerns into two separate issues, or one of them may be forgotten. There are a lot of changes to the rasterization algorithms on the way, I'm just struggling to find time to write tests for them: #336 But I will be sure to add tests for the integer case. As for if the line should contribute, GDAL has the same problem, and a more confusing solution which is to go positive on x and negative on y (or the opposite, I can't remember). We can do better than GDAL as we have all the flexibility of writing it in Julia. We just don't have the manpower to workshop these things much. So if you have time, please make another issue detailing the exact logic (even with plot mockups like above) of how you think rasterization should work in regards to lines and their interactions with the intervals of a pixel. If you write a good, logical spec I would implement it exactly. Working out what should happen is the hardest part and you seem to have more ideas than I do. Edit: We have the option of adding IntervalSets.jl intervals as the values in the lookup index, which should let us specify closed/open intervals for the pixels to make it work however we want. Even per pixel! It's also just something I haven't had time to implement |
Thanks for pointing out that the two issues should be separated. I just thought that they were of the same topic and highly related. Well, to make the discussions consecutive, I choose to give out my answer here instead of opening another issue. Actually I am the author of Mikrubi.jl, a package you probably still remember. (Thank you so much for the suggestions! I haven't make the repo public here just because I am not sure whether I would break any rules if I do so.) In that package, I wrote the rasterization routine carefully. The logic there is:
The problem we discuss here lies in the last step. To illustrate how I handle that, please allow me to paste my core implementation verbatim: function strokeedge(x1, y1, x2, y2, xb, yb)
pixels = Tuple{Int,Int}[]
x1, y1, x2, y2 = trimedge(x1, y1, x2, y2, xb, yb)
isnan(x1) && return pixels
xr = floor(Int, min(x1, x2)) + 1 : ceil(Int, max(x1, x2))
yr = floor(Int, y1) + 1 : ceil(Int, y2)
(length(xr) == 0 || length(yr) == 0) && return pixels
(length(xr) == 1 || length(yr) == 1) && return I.(xr, yr)
y2x = interpolate(y1, y2, x1, x2)
if x1 < x2
xleft = min(xb, max(0, x1))
for y = yr[begin:end-1]
xright = min(xb, max(0, y2x(y)))
append!(pixels, I.(floor(Int, xleft) + 1 : ceil(Int, xright), y))
xleft = xright
end
xright = min(xb, max(0, x2))
return append!(pixels,
I.(floor(Int, xleft) + 1 : ceil(Int, xright), ceil(Int, y2)))
else
xright = min(xb, max(0, x1))
for y = yr[begin:end-1]
xleft = min(xb, max(0, y2x(y)))
append!(pixels, I.(floor(Int, xleft) + 1 : ceil(Int, xright), y))
xright = xleft
end
xleft = min(xb, max(0, x2))
return append!(pixels,
I.(floor(Int, xleft) + 1 : ceil(Int, xright), ceil(Int, y2)))
end
end In this function, Not sure if there are better or smarter ways to do that, but at least I think the codes work well. I hope this can help a little bit. |
Oh right, I have actually read your method! In practice we need a much faster method, thus using the Bresenham style algorithm, for use cases of the scale we use in #336 But we do need to get it right here. If you can, please check out the branch in the pull request (instructions in the description at the top) and see if it fixes your problem. Run this in the REPL to use: (But I will also look at incorporating your algorithm instead, it does seem simpler than mine) |
The internal method is now called: Rasters._burn_line!(raster, (start=(x=x1, y=y1), stop=(x=x2, y=y2)), true) And no it still has this bug. Ok, I will try to fix before next week. Its strange how it only occurs with exact integer values and is otherwise fine. |
Good to know the progress in the morning! It is always reasonable to care about time. I admit that my overall rasterization routine is not as fast as yours, but in this very step (your And yes, I think the Let me explain my idea. If Without loss of generality, let us assume I would be happy if my codes and explanation can be of a little use to you and to the package. But I am still not very familiar with the function interfaces in Rasters.jl, so I am afraid that I would break your design somewhere if I randomly modify your codes by myself. |
The amount of computation can be further reduced. The trick is to check whether the slope of the edge has an absolute value greater than one. When |
Ok, I think I see what you want. If I understand correctly, a vertical line on the line between two pixels does not actually touch any pixel?
I don't think any of these methods is right or wrong, they just serve different requirements. So I think we should just handle all of them? We need to correctly define the intervals we are working with and what "touches" means in relation to the intervals of the pixels. I think you want pixels to cover the area (x1, x2), (y1, y2). We currently choose [x1, x2), [y1, y2), and gdal chooses something like [x1, x2), (y1, y2]. [x1, x2], [y1, y2] could also be useful for some things (edit: and maybe the ArcGIS default). Probably using https://github.com/JuliaMath/IntervalSets.jl is the right way, and it's already a dependency here. We can allow the lookup array to contain intervals, and also add an |
Oh and with the plot, there is a pixel missing for the 0.99999 case because the line is shorter with a step of 0.99999 and does not actually touch that pixel at all. The end is just rounded so it looks like it crosses. |
Yes, I fully agree with you. They are just different needs for different situations. "Touching inside" is only my personal preference for a special application. I didn't mean the current criterion you apply here is incorrect. Anyway, it will be the best if users can freely determine how the package behaves for cases involving pixel boundaries. And I guess this can probably be made by simply changing the directions of rounding ( Sorry I misunderstood the plot titles --- I thought they were values of |
Ok, I'm very happy that we agree on this and have a way forward. I will look into the implementation over the weekend. Making this change also solves why I have taken so long to merge the Now we should be able to get exactly the same results as GDAL and test against them by specifying the same interval rules. In terms of your use case, it will likely take me a month or more to have this properly finished, tested and merged - I hope that is not too long. The fix to the first problem (the actual bug) will be merged soon. |
I am really happy that I learned a lot of new things from your codes and our discussion here. Thank you for considering my request! I am looking forward to the new version. Please take your time --- the current behavior is already sufficient for most practical cases. |
Great. Good to hear there is no rush 😅 If you would like to collaborate on anything here in future (even just by making more issues with your needs) we would benefit from the amount of thought and care you put questions like this. |
No problem! I will come here again if I have any other ideas worth discussing. |
Dear Mr. Schouten and other developers,
I am new to this package. Please forgive me for my recklessness if the following issues are due to my misuse or misunderstanding.
I guess that there may be a bug lying in the rasterization routine for (multi)polygons when
boundary
is set to:touches
, or more specifically, in the function_fill_line!
. Here is an minimal working example (running in Julia v1.8.2):In the given raster grid, the pixels are bounded by adjacent integer coordinates. Thus we may denote a pixel by its right-above coordinates, e.g., (3,5) for the pixel ranging from 2 to 3 on x-axis and from 4 to 5 on y-axis. In this example, the straight line connecting (1.00, 4.50) and (4.75, 0.75) should touches exactly eight pixels including
But the current routine returns a different set of pixels (Figure), i.e.,
The unexpected behavior is observed when
x1
is an integer. Whenx1
is something like0.99
or1.01
, the package returns the correct rasterization results. Therefore I suspect that the problem is caused by rounding of floating numbers, although I am not yet able to locate it.Another issue, which I don't think is a bug like the one above, is about symmetry. Please consider an edge totally overlapping with some line of the grid, e.g., from (1.2, 4.0) to (3.4, 4.0). Now the rasterization functions return the pixels on the positive side of the line, i.e., (2,5), (3,5), and (4,5). But positivity may sometimes be arbitrary --- thus the preference may not be fair. From my point of view, it may be better to have these edges contribute no pixels in the rasterization process.
The text was updated successfully, but these errors were encountered: