r.proj.seg thins along null areas and raster bounds for bilinear and cubic methods

[kyngchaos]

Due to the nature of the bilinear and cubic interpolation methods, when any of the surrounding matrix of cells is null, the interpolated cell is set to null, instead of making do with what is available. With the bilinear, thinning occurs to the right and bottom of the interpolated cell, with cubic it will occur on all sides.

Since the bounds of a raster have effectively null cells outside the bounds, it also occurs at raster boundaries. This is especially noticable when projecting latlong rasters across the +-180 meridian -- there will be a line of cells missing where 180 E joins with 180 W.

This is a change from the previous r.proj, which set any surrounding nulls to the "nearest" cell temporarily.

Patches attached to handle both issues. Some notes:

The wraparound case I decided to limit only to latlong input projections. This does not need any option to enable/disable it, but each r.proj interpolation method needs an extra flag in its function call (see r.proj.h). It also only needs to wraparound when the interpolated cell is within 2 cells of +-180 long or +-90 lat, (oops, I just realized that +-90 lat doesn't wrap the same as +-180 long). The +-180 and +-90 values are hardcoded in, maybe there is a better way to handle maximum latlong extents (like the case where it's 0-360 long)?

The null filling case could be optional, with a module flag. This would have to be passed to each r.proj interpolation function with another flag alongside the wraparound flag.

For now, the surrounding cells are filled with the "nearest" input cell to the interpolated cell, like the original r.proj did. This works for the bilinear interp, but could cause odd interpolation in the cubic method, so the cubic null filling needs some work.

[kyngchaos]

I forgot to mention for the wraparound - if there is no wraparound (either latlong bounds less than max latlong, or not latlong input), at the edges of the bounds the interp matrix will "collapse". That is, the border row or column will be duplicated. This feature could be considered null filling, and dependent on a fill null flag, if that is made optional.

[kyngchaos]

OK, N/S wraparound got too complicated. For now I made those collapse.

Problem: the interp function would need to know something about the coordinate extents of the input, but all it knows is cell numbers. (that's why I have the wrap test in main.c and pass a flag)

[glynn]

Replying to kyngchaos:

For now, the surrounding cells are filled with the "nearest" input cell to the interpolated cell, like the original r.proj did. This works for the bilinear interp, but could cause odd interpolation in the cubic method, so the cubic null filling needs some work.

The existing behaviour regarding actual nulls is correct, and will be preserved. If you want nulls filled, use e.g. r.fillnulls. Don't go stuffing ad-hoc fudges into r.proj.

The issue with lat/lon wrapping is valid, although personally I think that this is a fundamental defect with GRASS. GRASS only supports lat/lon to the extent that every single module deals with it by itself.

[kyngchaos]

Replying to glynn:

The existing behaviour regarding actual nulls is correct, and will be preserved. If you want nulls filled, use e.g. r.fillnulls. Don't go stuffing ad-hoc fudges into r.proj.

The problem is not when the target input cell is null - that I don't have a problem with. It's when input cell is non-null, but some or all of the cells in the surrounding interpolation matrix are null. The basic interpolation algorithms need values in all cells, and I'm suggesting an option to fake those surrounding values (in the worst case, the interpolation drops down to a nearest neighbor interp).

It comes down to - if there is a cell in the input, I want a cell there in the output. Hmmm, I suppose this would work (could be done from the main.c loop), and it would not mess with the algorithms:

• if cubic option, interpolate cell cubic
• if no value yet, or bilinear option, interpolate cell bilinear
• if no value yet, or nearest option, interpolate cell nearest

This would attempt lower interpolations if the higher one did not work because of nulls in the matrix.

[glynn]

Replying to kyngchaos:

The existing behaviour regarding actual nulls is correct, and will be preserved. If you want nulls filled, use e.g. r.fillnulls. Don't go stuffing ad-hoc fudges into r.proj.

The problem is not when the target input cell is null - that I don't have a problem with. It's when input cell is non-null, but some or all of the cells in the surrounding interpolation matrix are null. The basic interpolation algorithms need values in all cells,

Yes. Yes, they do need values in all surrounding cells. This is not a bug.

and I'm suggesting an option to fake those surrounding values (in the worst case, the interpolation drops down to a nearest neighbor interp).

It comes down to - if there is a cell in the input, I want a cell there in the output.

This assumes that the output cells correspond to input cells. They don't (at least, not for bilinear and bicubic). They correspond to the input surface, and the surface is undefined within 1 or 2 cells of the centre of any null cell.

Hmmm, I suppose this would work (could be done from the main.c loop), and it would not mess with the algorithms:

• if cubic option, interpolate cell cubic
• if no value yet, or bilinear option, interpolate cell bilinear
• if no value yet, or nearest option, interpolate cell nearest

IOW, the equivalent of running r.proj with each of the three methods, then r.patch-ing the results.

This would attempt lower interpolations if the higher one did not work because of nulls in the matrix.

The requested interpolation always works. Sometimes the (correct) result is null; that isn't the same thing as not working.

In any case, if you want to implement this, the least invasive approach would be to extend the "menu" array with:

    {p_bilinear_fudge, "bilinear_fudge", "bilinear with fallback"},
    {p_cubic_fudge, "cubic_fudge", "cubic with fallback"},

and add the corresponding functions.

This would ensure that the existing code paths are completely unaffected.

I'm wary about doing anything which might possibly affect existing usage, as bogus results typically won't be obvious unless you explicitly check with e.g. r.slope.aspect.

The previous calculation was used for at least a decade without anyone noticing that it was nonsense (see r22190).

[kyngchaos]

Replying to glynn:

This assumes that the output cells correspond to input cells. They don't (at least, not for bilinear and bicubic). They correspond to the input surface, and the surface is undefined within 1 or 2 cells of the centre of any null cell.

Ah, this clarifies the logic of the algorithms.

The requested interpolation always works. Sometimes the (correct) result is null; that isn't the same thing as not working.

...

and add the corresponding functions.

This would ensure that the existing code paths are completely unaffected.

I'm wary about doing anything which might possibly affect existing usage, as bogus results typically won't be obvious unless you explicitly check with e.g. r.slope.aspect.

Understood, though I was thinking of a flag, but that would have gotten real messy.

OK, I'll work on separate methods. Is it OK for an interp function here to call another interp function? ie, instead of duplicating the code in cubic, bilinear and nearest, I could call each, as I outlined to do in main.c. It would be a little (lot?) slower, but easier to maintain.

Any comments on the wraparound stuff? I noticed Markus Metz's comment on the list about extending latlong locations beyond +-180 deg. This would make my wraparound patch more difficult. Though your idea of Euclidifying regions might solve it (but that sounds like something for GRASS 7).

[glynn]

Replying to kyngchaos:

OK, I'll work on separate methods. Is it OK for an interp function here to call another interp function? ie, instead of duplicating the code in cubic, bilinear and nearest, I could call each, as I outlined to do in main.c.

Yes.

It would be a little (lot?) slower, but easier to maintain.

Any slowdown probably won't be significant.

One option is to inline p_nearest(), as that's trivial, and if it produces a null result, p_bilinear() and p_cubic() will also produce a null result. If p_cubic() succeeds, the overhead of p_cubic() (16 cells and 5 cubic evaluations) will dwarf the (unnecessary) p_nearest() calculation, and if p_cubic() fails, you've saved a fair amount of computation.

Any comments on the wraparound stuff? I noticed Markus Metz's comment on the list about extending latlong locations beyond +-180 deg. This would make my wraparound patch more difficult. Though your idea of Euclidifying regions might solve it (but that sounds like something for GRASS 7).

Adding special-case code for lat/lon wraparound to individual modules doesn't make sense, IMHO. I gave up on trying to get this right in r.resamp.interp (the same issue will also affect r.resamp.stats, and probably similar modules). Pushing this into the libraries should be simple enough, at least for rasters. This will have to be reserved for 7.x, though.

[kyngchaos]

Replying to glynn:

One option is to inline p_nearest(), as that's trivial, and if it produces a null result, p_bilinear() and p_cubic() will also produce a null result. If p_cubic() succeeds, the overhead of p_cubic() (16 cells and 5 cubic evaluations) will dwarf the (unnecessary) p_nearest() calculation, and if p_cubic() fails, you've saved a fair amount of computation.

Sounds good. I know little to nothing about inlining :(

Adding special-case code for lat/lon wraparound to individual modules doesn't make sense, IMHO.

It works well enough for my frequent projection need to wrap at 180. I can leave it at a personal customization, but it would be nice if some form could be worked into dev6, maybe as an option, so it doesn't change current behavior.

Pushing this into the libraries should be simple enough, at least for rasters. This will have to be reserved for 7.x, though.

Should this ticket be split off to separate the null fallback and wraparound ideas?

[glynn]

Replying to kyngchaos:

Should this ticket be split off to separate the null fallback and wraparound ideas?

Wraparound affects far more than just r.proj. I'm not sure that a ticket is the correct approach for fundamental design issues.

[kyngchaos]

OK, something like this? Ignoring the wraparound for now, new set of patches and new source for fallback interp functions.

I made all interp functions return a status value:

• 0 on successfull non-null interpolation OR if the nearest input cell is null -> no further interpolation required
• 1 if the interpolated output is null -> fallback to the next interpolation level

I don't know about inlining, so feel free to suggest how to do that, if you think it will speed it up significantly.

[glynn]

Replying to kyngchaos:

OK, something like this? Ignoring the wraparound for now, new set of patches and new source for fallback interp functions.

I made all interp functions return a status value:

• 0 on successfull non-null interpolation OR if the nearest input cell is null -> no further interpolation required
• 1 if the interpolated output is null -> fallback to the next interpolation level

There is no need to modify the existing functions; just check whether the calculated value is null.

Also, it would have been better to attach the changes as a single patch, i.e. "svn add" the new files then use "svn diff" to create a patch containing all of the changes.

I don't know about inlining, so feel free to suggest how to do that, if you think it will speed it up significantly.

Actually, I doubt that it matters.

[kyngchaos]

Replying to glynn:

There is no need to modify the existing functions; just check whether the calculated value is null.

Checking if the output is null doesn't tell us if it's null because the nearest input cell is null or because it interpolated to null from nulls in the matrix. If the nearest input cell is null, there is no point to try simpler interpolations (this mainly affects trying bilinear after cubic, for saving any computaion time).

I can try some speed tests Monday to see if this is worthwhile.

Also, it would have been better to attach the changes as a single patch, i.e. "svn add" the new files then use "svn diff" to create a patch containing all of the changes.

I could apply them to svn if they're OK.

[kyngchaos]

Replying to glynn off-trac:

Try nearest first; if that returns null, then you know that bilinear and bicubic will also return null, so you can just abort at that point. Otherwise, try bicubic then bilinear.

... ah, I was thinking it would be repeating a call to nearest, to apply that at the end if it actually needed it, but I suppose I could save that first value and reuse it when needed.

[kyngchaos]

Now I remember - to get the value of the cell in my added interp functions, after calling p_nearest to set the cell value, I'd be duplicating a most of p_nearest(): floor the row/col for the cell index, checking if it's out of bounds, getting the value with that CPTR() macro, and testing if it's null.

Hmmm, I suppose duplicate p_nearest completely and skip calling it? Is this what you meant by inlining? (I thought it would be some complex macro)

[kyngchaos]

bilinear and cubic fallback methods added to dev6 (r35734) and trunk (r35735).

Do you think it's OK to add them to release 6.4?

P.S. After doing some timing tests (7500x5000 cells), I realized that disk IO is the limiting factor here. duh. So there is not the 4x and 10x differences from bilinear and cubic to nearest, more like 1.01x. So my worries about speed were silly :S Though larger regions may show more significant differences in speed.

[hamish]

Replying to kyngchaos:

bilinear and cubic fallback methods added to dev6 (r35734) and trunk (r35735).

Do you think it's OK to add them to release 6.4?

can you summarize the new method for those of us who haven't followed this that closely?

is it a critical bug fix? if not then lets play conservative please, as we are already up to RC3.

FWIW I do a lot of reprojecting things like 8-color (ie CELL category 1-8) nautical charts where using anything but nearest cat is wrong, and automatic fall-over to another method will break the result.

thanks, Hamish

[kyngchaos]

I updated the description.html to describe it. As Glynn pointed out, it's not really a bug fix, the original change in behavior that removed this "feature" was itself a fix so that it strictly followed the algorithms.

I can live with leaving it out of 6.4, for my own use I'll just backport it locally then.

It doesn't go up the scale of interpolation complexity. Nearest will always be nearest. Also, it's a separate choice - if you want a strict bilinear or cubic, that's what you'll get. No change in current behavior.