[GRASS-dev] [GRASS GIS] #465: r.proj.seg thins along null areas and raster bounds for bilinear and cubic methods

GRASS GRASS developer
1

#465: r.proj.seg thins along null areas and raster bounds for bilinear and cubic
methods
-----------------------+----------------------------------------------------
Reporter: kyngchaos | Owner: grass-dev@lists.osgeo.org
     Type: defect | Status: new
Priority: major | Milestone:
Component: Raster | Version: svn-develbranch6
Keywords: | Platform: All
      Cpu: All |
-----------------------+----------------------------------------------------
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.

--
Ticket URL: <http://trac.osgeo.org/grass/ticket/465&gt;
GRASS GIS <http://grass.osgeo.org>

2

#465: r.proj.seg thins along null areas and raster bounds for bilinear and cubic
methods
------------------------+---------------------------------------------------
  Reporter: kyngchaos | Owner: grass-dev@lists.osgeo.org
      Type: defect | Status: new
  Priority: major | Milestone:
Component: Raster | Version: svn-develbranch6
Resolution: | Keywords:
  Platform: All | Cpu: All
------------------------+---------------------------------------------------
Comment (by 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.

--
Ticket URL: <http://trac.osgeo.org/grass/ticket/465#comment:1&gt;
GRASS GIS <http://grass.osgeo.org>

3

#465: r.proj.seg thins along null areas and raster bounds for bilinear and cubic
methods
------------------------+---------------------------------------------------
  Reporter: kyngchaos | Owner: grass-dev@lists.osgeo.org
      Type: defect | Status: new
  Priority: major | Milestone:
Component: Raster | Version: svn-develbranch6
Resolution: | Keywords:
  Platform: All | Cpu: All
------------------------+---------------------------------------------------
Comment (by 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)

--
Ticket URL: <http://trac.osgeo.org/grass/ticket/465#comment:2&gt;
GRASS GIS <http://grass.osgeo.org>

4

#465: r.proj.seg thins along null areas and raster bounds for bilinear and cubic
methods
------------------------+---------------------------------------------------
  Reporter: kyngchaos | Owner: grass-dev@lists.osgeo.org
      Type: defect | Status: new
  Priority: major | Milestone:
Component: Raster | Version: svn-develbranch6
Resolution: | Keywords:
  Platform: All | Cpu: All
------------------------+---------------------------------------------------
Comment (by glynn):

Replying to [ticket:465 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.

--
Ticket URL: <http://trac.osgeo.org/grass/ticket/465#comment:3&gt;
GRASS GIS <http://grass.osgeo.org>

5

#465: r.proj.seg thins along null areas and raster bounds for bilinear and cubic
methods
------------------------+---------------------------------------------------
  Reporter: kyngchaos | Owner: grass-dev@lists.osgeo.org
      Type: defect | Status: new
  Priority: major | Milestone:
Component: Raster | Version: svn-develbranch6
Resolution: | Keywords:
  Platform: All | Cpu: All
------------------------+---------------------------------------------------
Comment (by kyngchaos):

Replying to [comment:3 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.

--
Ticket URL: <http://trac.osgeo.org/grass/ticket/465#comment:4&gt;
GRASS GIS <http://grass.osgeo.org>

6

#465: r.proj.seg thins along null areas and raster bounds for bilinear and cubic
methods
--------------------------+-------------------------------------------------
  Reporter: kyngchaos | Owner: grass-dev@lists.osgeo.org
      Type: enhancement | Status: new
  Priority: major | Milestone:
Component: Raster | Version: svn-develbranch6
Resolution: | Keywords:
  Platform: All | Cpu: All
--------------------------+-------------------------------------------------
Changes (by glynn):

  * type: defect => enhancement

Comment:

Replying to [comment:4 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).

--
Ticket URL: <http://trac.osgeo.org/grass/ticket/465#comment:5&gt;
GRASS GIS <http://grass.osgeo.org>

7

#465: r.proj.seg thins along null areas and raster bounds for bilinear and cubic
methods
--------------------------+-------------------------------------------------
  Reporter: kyngchaos | Owner: grass-dev@lists.osgeo.org
      Type: enhancement | Status: new
  Priority: major | Milestone:
Component: Raster | Version: svn-develbranch6
Resolution: | Keywords:
  Platform: All | Cpu: All
--------------------------+-------------------------------------------------
Comment (by kyngchaos):

Replying to [comment:5 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).

--
Ticket URL: <http://trac.osgeo.org/grass/ticket/465#comment:6&gt;
GRASS GIS <http://grass.osgeo.org>

8

#465: r.proj.seg thins along null areas and raster bounds for bilinear and cubic
methods
--------------------------+-------------------------------------------------
  Reporter: kyngchaos | Owner: grass-dev@lists.osgeo.org
      Type: enhancement | Status: new
  Priority: major | Milestone:
Component: Raster | Version: svn-develbranch6
Resolution: | Keywords:
  Platform: All | Cpu: All
--------------------------+-------------------------------------------------
Comment (by glynn):

Replying to [comment:6 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.

--
Ticket URL: <http://trac.osgeo.org/grass/ticket/465#comment:7&gt;
GRASS GIS <http://grass.osgeo.org>

9

#465: r.proj.seg thins along null areas and raster bounds for bilinear and cubic
methods
--------------------------+-------------------------------------------------
  Reporter: kyngchaos | Owner: grass-dev@lists.osgeo.org
      Type: enhancement | Status: new
  Priority: major | Milestone:
Component: Raster | Version: svn-develbranch6
Resolution: | Keywords:
  Platform: All | Cpu: All
--------------------------+-------------------------------------------------
Comment (by kyngchaos):

Replying to [comment:7 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 :frowning:

> 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?

--
Ticket URL: <http://trac.osgeo.org/grass/ticket/465#comment:8&gt;
GRASS GIS <http://grass.osgeo.org>

10

#465: r.proj.seg thins along null areas and raster bounds for bilinear and cubic
methods
--------------------------+-------------------------------------------------
  Reporter: kyngchaos | Owner: grass-dev@lists.osgeo.org
      Type: enhancement | Status: new
  Priority: major | Milestone:
Component: Raster | Version: svn-develbranch6
Resolution: | Keywords:
  Platform: All | Cpu: All
--------------------------+-------------------------------------------------
Comment (by glynn):

Replying to [comment:8 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.

--
Ticket URL: <http://trac.osgeo.org/grass/ticket/465#comment:9&gt;
GRASS GIS <http://grass.osgeo.org>

11

#465: r.proj.seg thins along null areas and raster bounds for bilinear and cubic
methods
--------------------------+-------------------------------------------------
  Reporter: kyngchaos | Owner: grass-dev@lists.osgeo.org
      Type: enhancement | Status: new
  Priority: major | Milestone:
Component: Raster | Version: svn-develbranch6
Resolution: | Keywords:
  Platform: All | Cpu: All
--------------------------+-------------------------------------------------
Comment (by 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.

--
Ticket URL: <http://trac.osgeo.org/grass/ticket/465#comment:10&gt;
GRASS GIS <http://grass.osgeo.org>

12

#465: r.proj.seg thins along null areas and raster bounds for bilinear and cubic
methods
--------------------------+-------------------------------------------------
  Reporter: kyngchaos | Owner: grass-dev@lists.osgeo.org
      Type: enhancement | Status: new
  Priority: major | Milestone:
Component: Raster | Version: svn-develbranch6
Resolution: | Keywords:
  Platform: All | Cpu: All
--------------------------+-------------------------------------------------
Comment (by glynn):

Replying to [comment:10 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.

--
Ticket URL: <http://trac.osgeo.org/grass/ticket/465#comment:11&gt;
GRASS GIS <http://grass.osgeo.org>

13

#465: r.proj.seg thins along null areas and raster bounds for bilinear and cubic
methods
--------------------------+-------------------------------------------------
  Reporter: kyngchaos | Owner: grass-dev@lists.osgeo.org
      Type: enhancement | Status: new
  Priority: major | Milestone:
Component: Raster | Version: svn-develbranch6
Resolution: | Keywords:
  Platform: All | Cpu: All
--------------------------+-------------------------------------------------
Comment (by kyngchaos):

Replying to [comment:11 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.

--
Ticket URL: <http://trac.osgeo.org/grass/ticket/465#comment:12&gt;
GRASS GIS <http://grass.osgeo.org>

14

[Replying here as Trac is not responding at this time.]

GRASS GIS wrote:

> 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.

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.

--
Glynn Clements <glynn@gclements.plus.com>

15

#465: r.proj.seg thins along null areas and raster bounds for bilinear and cubic
methods
--------------------------+-------------------------------------------------
  Reporter: kyngchaos | Owner: grass-dev@lists.osgeo.org
      Type: enhancement | Status: new
  Priority: major | Milestone:
Component: Raster | Version: svn-develbranch6
Resolution: | Keywords:
  Platform: All | Cpu: All
--------------------------+-------------------------------------------------
Comment (by 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.

--
Ticket URL: <http://trac.osgeo.org/grass/ticket/465#comment:13&gt;
GRASS GIS <http://grass.osgeo.org>

16

#465: r.proj.seg thins along null areas and raster bounds for bilinear and cubic
methods
--------------------------+-------------------------------------------------
  Reporter: kyngchaos | Owner: grass-dev@lists.osgeo.org
      Type: enhancement | Status: new
  Priority: major | Milestone:
Component: Raster | Version: svn-develbranch6
Resolution: | Keywords:
  Platform: All | Cpu: All
--------------------------+-------------------------------------------------
Comment (by 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)

--
Ticket URL: <http://trac.osgeo.org/grass/ticket/465#comment:14&gt;
GRASS GIS <http://grass.osgeo.org>

17

#465: r.proj.seg thins along null areas and raster bounds for bilinear and cubic
methods
--------------------------+-------------------------------------------------
  Reporter: kyngchaos | Owner: grass-dev@lists.osgeo.org
      Type: enhancement | Status: new
  Priority: major | Milestone:
Component: Raster | Version: svn-develbranch6
Resolution: | Keywords:
  Platform: All | Cpu: All
--------------------------+-------------------------------------------------
Comment (by 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.

--
Ticket URL: <http://trac.osgeo.org/grass/ticket/465#comment:15&gt;
GRASS GIS <http://grass.osgeo.org>

18

#465: r.proj.seg thins along null areas and raster bounds for bilinear and cubic
methods
--------------------------+-------------------------------------------------
  Reporter: kyngchaos | Owner: grass-dev@lists.osgeo.org
      Type: enhancement | Status: new
  Priority: major | Milestone:
Component: Raster | Version: svn-develbranch6
Resolution: | Keywords:
  Platform: All | Cpu: All
--------------------------+-------------------------------------------------
Comment (by hamish):

Replying to [comment:15 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

--
Ticket URL: <https://trac.osgeo.org/grass/ticket/465#comment:16&gt;
GRASS GIS <http://grass.osgeo.org>

19

#465: r.proj.seg thins along null areas and raster bounds for bilinear and cubic
methods
--------------------------+-------------------------------------------------
  Reporter: kyngchaos | Owner: grass-dev@lists.osgeo.org
      Type: enhancement | Status: new
  Priority: major | Milestone:
Component: Raster | Version: svn-develbranch6
Resolution: | Keywords:
  Platform: All | Cpu: All
--------------------------+-------------------------------------------------
Comment (by 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.

--
Ticket URL: <https://trac.osgeo.org/grass/ticket/465#comment:17&gt;
GRASS GIS <http://grass.osgeo.org>

20

#465: r.proj.seg thins along null areas and raster bounds for bilinear and cubic
methods
--------------------------+-------------------------------------------------
  Reporter: kyngchaos | Owner: grass-dev@lists.osgeo.org
      Type: enhancement | Status: new
  Priority: major | Milestone: 6.5.0
Component: Raster | Version: svn-develbranch6
Resolution: | Keywords: r.proj
  Platform: All | Cpu: All
--------------------------+-------------------------------------------------
Changes (by hamish):

  * keywords: => r.proj
  * milestone: => 6.5.0

--
Ticket URL: <http://trac.osgeo.org/grass/ticket/465#comment:18&gt;
GRASS GIS <http://grass.osgeo.org>