[GRASS-dev] bilinear and bicubic raster interpolation

In i.rectify, I have implemented the interpolation methods available
in r.proj, and while testing discovered an offset of one cell (row +
1, col + 1 from source to target) between the source and target
imagery after rectifying (40+ automatically collected GCPs, order=1,
sub-pixel RMS errors, source image not distorted). I corrected this
offset in i.rectify and suggest that all modules using
Rast_interp_linear/Rast_interp_bilinear/Rast_interp_cubic/Rast_interp_bicubic
should be checked. I am pretty sure that at least r.proj needs this
fix (if a one-cell shift is something to worry about). I think this
kind of correction is also done in r.resamp.interp.

Explanations in
https://trac.osgeo.org/grass/browser/grass/trunk/imagery/i.rectify/nearest.c#L23
https://trac.osgeo.org/grass/browser/grass/trunk/imagery/i.rectify/bilinear.c#L35
https://trac.osgeo.org/grass/browser/grass/trunk/imagery/i.rectify/cubic.c#L36

Markus M

Markus Metz wrote:

In i.rectify, I have implemented the interpolation methods available
in r.proj, and while testing discovered an offset of one cell (row +
1, col + 1 from source to target) between the source and target
imagery after rectifying (40+ automatically collected GCPs, order=1,
sub-pixel RMS errors, source image not distorted). I corrected this
offset in i.rectify and suggest that all modules using
Rast_interp_linear/Rast_interp_bilinear/Rast_interp_cubic/Rast_interp_bicubic
should be checked. I am pretty sure that at least r.proj needs this
fix (if a one-cell shift is something to worry about).

For MODIS satellite data the cells can be 4km tall, for GRIB data 1/2 degree
or full degree. So once reprojected from eg lat/lon to UTM, the cells are
still quite big and a 1 cell offset can be significant when trying to
ground-truth against in situ measurements.

also, i.rectify and gdalwarp share the same code source for order= matrix
calcs, so maybe we should check for the same there if that was touched.

thanks,
Hamish

Hamish wrote:

Markus Metz wrote:

In i.rectify, I have implemented the interpolation methods available
in r.proj, and while testing discovered an offset of one cell (row +
1, col + 1 from source to target) between the source and target
imagery after rectifying (40+ automatically collected GCPs, order=1,
sub-pixel RMS errors, source image not distorted). I corrected this
offset in i.rectify and suggest that all modules using
Rast_interp_linear/Rast_interp_bilinear/Rast_interp_cubic/Rast_interp_bicubic
should be checked. I am pretty sure that at least r.proj needs this
fix (if a one-cell shift is something to worry about).

For MODIS satellite data the cells can be 4km tall, for GRIB data 1/2 degree
or full degree. So once reprojected from eg lat/lon to UTM, the cells are
still quite big and a 1 cell offset can be significant when trying to
ground-truth against in situ measurements.

also, i.rectify and gdalwarp share the same code source for order= matrix
calcs, so maybe we should check for the same there if that was touched.

The transformation matrix code (crs.[c|h]) wasn't touched. This cell
shift really only applies to linear and cubic interpolation. My point
is that a reprojected or transformed or recalculated (in case of
resampling within the same projection) cell center (result is floating
point) of e.g. row 2.1 is in between row centers 1.5 and 2.5 because
1.5 < 2.1 < 2.5, and not in between row centers 2.5 and 3.5, thus (for
linear ip) rows 1 and 2 should be used and not rows 2 and 3, same for
column indices and cubic ip.

Markus M

Markus Metz wrote:

In i.rectify, I have implemented the interpolation methods available
in r.proj, and while testing discovered an offset of one cell (row +
1, col + 1 from source to target) between the source and target
imagery after rectifying (40+ automatically collected GCPs, order=1,
sub-pixel RMS errors, source image not distorted). I corrected this
offset in i.rectify and suggest that all modules using
Rast_interp_linear/Rast_interp_bilinear/Rast_interp_cubic/Rast_interp_bicubic
should be checked. I am pretty sure that at least r.proj needs this
fix (if a one-cell shift is something to worry about).

Actually, I think that r.proj[.seg] have half-cell errors, i.e. they
treat a cell's value as a sample at its north-west corner rather than
its centre. A test case with r.proj (7.0) appears to confirm this.
This applies to both bilinear and bicubic interpolation.

A column coordinate of e.g. 7.3 should interpolate between sample
values at 6.5 and 7.5 with an interpolation parameter of 0.8. This
uses cells 6 and 7 for bilinear interpolation and cells 5,6,7,8 for
bicubic interpolation.

This should be fixed by r43943.

However, there is another issue with r.proj (7.0) and r.proj.seg
(6.x): bilinear and cubic interpolation give a small number of garbage
results due to incorrect use of the caching mechanism (blocks are
being replaced while references remain). I'll fix this shortly.

In the linear and

I think this kind of correction is also done in r.resamp.interp.

r.resamp.interp and Rast_get_sample_* appear to be correct.

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

Glynn Clements wrote:

However, there is another issue with r.proj (7.0) and r.proj.seg
(6.x): bilinear and cubic interpolation give a small number of garbage
results due to incorrect use of the caching mechanism (blocks are
being replaced while references remain). I'll fix this shortly.

This should be fixed by r43944.

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

On Sun, Oct 17, 2010 at 11:40 AM, Glynn Clements
<glynn@gclements.plus.com> wrote:

Glynn Clements wrote:

However, there is another issue with r.proj (7.0) and r.proj.seg
(6.x): bilinear and cubic interpolation give a small number of garbage
results due to incorrect use of the caching mechanism (blocks are
being replaced while references remain). I'll fix this shortly.

This should be fixed by r43944.

The common backport question comes up for this one and r43943...
Indications welcome.

Markus

Markus Neteler wrote:

>> However, there is another issue with r.proj (7.0) and r.proj.seg
>> (6.x): bilinear and cubic interpolation give a small number of garbage
>> results due to incorrect use of the caching mechanism (blocks are
>> being replaced while references remain). I'll fix this shortly.
>
> This should be fixed by r43944.

The common backport question comes up for this one and r43943...
Indications welcome.

Yes to both. Note that 7.0's r.proj is r.proj.seg in 6.x, which may
complicate matters.

Also, I believe that the off-by-half issue (r43943) affects the 6.x
r.proj, which will need a similar fix (but which will need to be done
manually). r43944 is only applicable to the newer version (r.proj.seg
and 7.0), not the 6.x r.proj.

FWIW, I found the errors using:

r.plane --o name=plane dip=30 azimuth=45 easting=599505 northing=4921010 elevation=1453 type=double
r.proj --o --q in=plane location=spearfish57 mapset=glynn output=plane.n method=nearest
r.proj --o --q in=plane location=spearfish57 mapset=glynn output=plane.l method=bilinear
r.proj --o --q in=plane location=spearfish57 mapset=glynn output=plane.c method=cubic
r.mapcalc --o 'diff.n = plane.n - plane'
r.mapcalc --o 'diff.l = plane.l - plane'
r.mapcalc --o 'diff.c = plane.c - plane'
r.info -r diff.n
r.info -r diff.l
r.info -r diff.c

[I had to modify r.proj not to complain about the source and
destination being identical. There's no reason why this should cause
problems, although I can't see any use for it other than testing.]

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

On Sun, Oct 17, 2010 at 6:56 PM, Glynn Clements
<glynn@gclements.plus.com> wrote:

Markus Neteler wrote:

>> However, there is another issue with r.proj (7.0) and r.proj.seg
>> (6.x): bilinear and cubic interpolation give a small number of garbage
>> results due to incorrect use of the caching mechanism (blocks are
>> being replaced while references remain). I'll fix this shortly.
>
> This should be fixed by r43944.

The common backport question comes up for this one and r43943...
Indications welcome.

Yes to both. Note that 7.0's r.proj is r.proj.seg in 6.x, which may
complicate matters.

Slightly... ok I have done both as

G6.4:
http://trac.osgeo.org/grass/changeset/43953
http://trac.osgeo.org/grass/changeset/43956

G6.5:
http://trac.osgeo.org/grass/changeset/43954
http://trac.osgeo.org/grass/changeset/43955

Please verify.

Also, I believe that the off-by-half issue (r43943) affects the 6.x
r.proj, which will need a similar fix (but which will need to be done
manually). r43944 is only applicable to the newer version (r.proj.seg
and 7.0), not the 6.x r.proj.

AFAIK we have r.proj.seg everywhere.

FWIW, I found the errors using:

r.plane --o name=plane dip=30 azimuth=45 easting=599505 northing=4921010 elevation=1453 type=double
r.proj --o --q in=plane location=spearfish57 mapset=glynn output=plane.n method=nearest
r.proj --o --q in=plane location=spearfish57 mapset=glynn output=plane.l method=bilinear
r.proj --o --q in=plane location=spearfish57 mapset=glynn output=plane.c method=cubic
r.mapcalc --o 'diff.n = plane.n - plane'
r.mapcalc --o 'diff.l = plane.l - plane'
r.mapcalc --o 'diff.c = plane.c - plane'
r.info -r diff.n
r.info -r diff.l
r.info -r diff.c

[I had to modify r.proj not to complain about the source and
destination being identical. There's no reason why this should cause
problems, although I can't see any use for it other than testing.]

Markus

Markus Neteler wrote:

> Also, I believe that the off-by-half issue (r43943) affects the 6.x
> r.proj, which will need a similar fix (but which will need to be done
> manually). r43944 is only applicable to the newer version (r.proj.seg
> and 7.0), not the 6.x r.proj.

AFAIK we have r.proj.seg everywhere.

Right; in 6.x, raster/r.proj isn't built, and raster/r.proj.seg is
built as r.proj. In which case, raster/r.proj should probably be
removed if it isn't going to be fixed.

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

On Mon, Oct 18, 2010 at 12:04 PM, Glynn Clements
<glynn@gclements.plus.com> wrote:

Markus Neteler wrote:

> Also, I believe that the off-by-half issue (r43943) affects the 6.x
> r.proj, which will need a similar fix (but which will need to be done
> manually). r43944 is only applicable to the newer version (r.proj.seg
> and 7.0), not the 6.x r.proj.

AFAIK we have r.proj.seg everywhere.

Right; in 6.x, raster/r.proj isn't built, and raster/r.proj.seg is
built as r.proj. In which case, raster/r.proj should probably be
removed if it isn't going to be fixed.

Of there are no objections, I'll do so.

Perhaps we should also remove the other deprecated modules from 6.4.svn:

[neteler@north grass64]$ ls -1 */*/DEP*
display/d.font.freetype/DEPRECATED
display/d.rast.edit/DEPRECATED
display/d.text/DEPRECATED
display/d.text.freetype/DEPRECATED
raster/r.bilinear/DEPRECATED
raster/r.proj/DEPRECATED
raster/r.sun/DEPRECATED
vector/v.buffer/DEPRECATED
vector/v.parallel/DEPRECATED

?
Markus

Glynn Clements wrote:

Glynn Clements wrote:

However, there is another issue with r.proj (7.0) and r.proj.seg
(6.x): bilinear and cubic interpolation give a small number of garbage
results due to incorrect use of the caching mechanism (blocks are
being replaced while references remain). I'll fix this shortly.

This should be fixed by r43944.

and for i.rectify (7.0) in 43973

Markus M

Markus Neteler wrote:

Perhaps we should also remove the other deprecated modules from 6.4.svn:

raster/r.bilinear/DEPRECATED

raster/Makefile still includes r.bilinear in SUBDIRS. Although the
documentation labels it as deprecated, it was included in 6.4.0, so it
may be too late to remove it from 6.4.x.

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

Markus Neteler wrote:
> Perhaps we should also remove the other deprecated
> modules from 6.4.svn:

...

> raster/r.bilinear/DEPRECATED

Glynn:

raster/Makefile still includes r.bilinear in SUBDIRS.
Although the documentation labels it as deprecated, it was
included in 6.4.0, so it may be too late to remove it from
6.4.x.

s/may be/is/

until bugs in the new versions are fixed, I'd vote to leave these
two in place:
vector/v.buffer/DEPRECATED
vector/v.parallel/DEPRECATED

Hamish

Hamish wrote:

until bugs in the new versions are fixed, I'd vote to leave these
two in place:
vector/v.buffer/DEPRECATED
vector/v.parallel/DEPRECATED

That only helps people who build from source. If there are problems
with the new versions, perhaps the old versions should be built as
well?

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

Hamish:

> until bugs in the new versions are fixed, I'd vote to leave
> these two in place:
> vector/v.buffer/DEPRECATED
> vector/v.parallel/DEPRECATED

Glynn:

That only helps people who build from source.

right.

If there are problems with the new versions, perhaps the old
versions should be built as well?

I don't think that's needed, and anyone using it for comparison
& able to do anything about it would be building from source
anyway. ie I mean to use the old version as a delta, not a
fallback. or moreover, I'm not sure if the old version would
help.

semi-offline for a while,
Hamish