[GRASS-user] integration of sample.c algorithms into r.resample

Tue Aug 15 14:12:01 EDT 2006

On Monday 14 August 2006 18:06, Glynn Clements wrote:
> Dylan Beaudette wrote:
> > After a bit more testing I can see some interesting differences.
> >
> > using the customized r.bilinear code previously posted:
> >
> > #set the region
> > g.region res=10 -ap
> >
> > #cutout a sample
> > r.mapcalc "c = elev_meters"
> >
> > #"resample" 10m DEM data to 1m
> >
> > # 1. via r.bilinear method=nearest
> > r.bilinear in=c out=c_n_1m method=nearest
> > # 2. via r.bilinear method=bilinear
> > r.bilinear in=c out=c_bc_1m method=bicubic
> > # 3. via r.bilinear method=bicubic
> > r.bilinear in=c out=c_b_1m method=bilinear
> >
> > #visual comparison of slope computed at 1m
> > r.slope.aspect elev=c_b_1m slope=c_b_1m.slope
> > r.slope.aspect elev=c_bc_1m slope=c_bc_1m.slope
> > r.slope.aspect elev=c_n_1m slope=c_n_1m.slope
> >
> > #reset the colors for each map
> > r.colors map=c_n_1m.slope color=gyr
> > r.colors map=c_b_1m.slope color=gyr
> > r.colors map=c_bc_1m.slope color=gyr
> >
> > #see images:
> > ##NN
> > http://169.237.35.250/~dylan/temp/raster_sampling/c_n_1m.slope.png
> >
> > ##Bilinear
> > http://169.237.35.250/~dylan/temp/raster_sampling/c_b_1m.slope.png
> >
> > ##Bicubic
> > http://169.237.35.250/~dylan/temp/raster_sampling/c_bc_1m.slope.png
> >
> >
> > note that the slope computed from the bicubic algorithm looks very
> > similar to the NN slope map -- all computed via the new r.bilinear code.
> > I was under the impression that the cubic output would be similar to the
> > bilinear, but "smoother"...
>
> The bicubic interpolation formula was taken directly from
> lib/gis/sample.c; it wouldn't surprise me if it was wrong.
>
> Hmm.
>
> c = v * (v * (v * (c3 - c2 + c1 - c0) + (c2 - c3 - 2 * c1 + 2 * c0)) + (c2
> - c0)) + c1 = v^3 * (c3 - c2 + c1 - c0) + v^2 * (c2 - c3 - 2*c1 + 2*c0) + v
> * (c2 - c0) + c1
>
> at v=0 => c1
> at v=1 => c2
>
> OK.
>
> dc/dv = 3 * v^2 * (c3 - c2 + c1 - c0) + 2 * v * (c2 - c3 - 2*c1 + 2*c0) +
> (c2 - c0) = v^2 * (3*c3 - 3*c2 + 3*c1 - 3*c0) + v * (2*c2 - 2*c3 - 4*c1 +
> 4*c0) + (c2 - c0)
>
> at v=0 => c2 - c0
> at v=1 => c3 - c1
>
> Not OK.
>
> The actual gradients should be half that (the distance between c2 and
> c0 and between c3 and c1 is 2 grid cells). It's doubling the slope of
> the boundary tangents.
>
> We want c = k3 * v^3 + k2 * v^2 + k1 * v + k0 such that:
>
> 	c(0) = c1
> 	c(1) = c2
> 	c'(0) = (c2 - c0)/2
> 	c'(1) = (c3 - c1)/2
>
> The defintion of c gives:
>
> 	c(0) = k0
> 		=> k0 = c1
> 	c(1) = k3 + k2 + k1 + k0
> 		=> k3 + k2 + k1 + k0 = c2
>
> Differentiating w.r.t. v giv
>
> 	c' = 3*k3 * v^2 + 2*k2 * v + k1
>
> 	c'(0) = k1
> 		=> k1 = (c2 - c0)/2
> 	c'(1) = 3*k3 + 2*k2 + k1
> 		=> 3*k3 + 2*k2 + k1 = (c3 - c1)/2
>
> IOW:
>
> 	3*k3 + 2*k2 + k1 = (c3 - c1)/2
> 	k3 + k2 + k1 + k0 = c2
> 	k1 = (c2 - c0)/2
> 	k0 = c1
>
> solving gives:
>
> 	k0 = c1
> 	k1 = c2/2 - c0/2
> 	k2 = -c3/2 + 2*c2 - 5*c1/2 + c0
> 	k3 = c3/2 - 3*c2/2 + 3*c1/2 - c0/2
>
> Try changing the relevant lines to:
>
> 	DCELL c0 = (u * (u * (u * (c30 - 3*c20 + 3*c10 - c00) + (-c30 + 2*c20  -
> 5*c10 + 2*c00)) + (c20 - c00)) + 2*c10)/2; DCELL c1 = (u * (u * (u * (c31 -
> 3*c21 + 3*c11 - c01) + (-c31 + 2*c21  - 5*c11 + 2*c01)) + (c21 - c01)) +
> 2*c11)/2; DCELL c2 = (u * (u * (u * (c32 - 3*c22 + 3*c12 - c02) + (-c32 +
> 2*c22  - 5*c12 + 2*c02)) + (c22 - c02)) + 2*c12)/2; DCELL c3 = (u * (u * (u
> * (c33 - 3*c23 + 3*c13 - c03) + (-c33 + 2*c23  - 5*c13 + 2*c03)) + (c23 -
> c03)) + 2*c13)/2; DCELL c  = (v * (v * (v * (c3  - 3*c2  + 3*c1  - c0 ) +
> (-c3  + 2*c2   - 5*c1  + 2*c0 )) + (c2  - c0 )) + 2*c1 )/2;

Some further refinements:

1. when using your new equations, i get correct results when my input raster 
actually extends beyond the current region, and when I move from small res to 
large res : 10m to 100m. 

2. when I try the same, but with a raster that _does not_ extend beyond the 
current region, i get odd results, from what looks like accessing an array 
with a bad index (?):

(
printed by code [1]
)

[North East]
[4039050.000000 665550.000000]
[R C] [V U]
[0 7] [0 1071644672]
[R0 R1 R2 R3] maprow_f 
[-21 -20 -19 -18] -19.500000
0.000000 0.000000 0.000000 0.000000 
0.000000 0.000000 0.000000 0.000000 
0.000000 0.000000 0.000000 0.000000 
0.000000 
73518217792023544336403454714494196428503373132106328150814719669703774177598362215711000331154672582307817565157691
43262327814427811332838925345097477936939540677015455495010019939540775036113228922258615076583985086816808075264.000000 
0.00
0000 0.000000 
----------------------------------------
c0 c1 c2 c3   ::  c
-0.000000 -459488861200147152102521591965588727678146082075664550942591997935648588609989763848193752069716703639423859782235
571453895488401738208302432834068592371058721292313465968438126246221298439757076807641163442286499067926050504704.000000 -0.
000000 -0.000000   ::  -25846248442508277305766839548064365931895717116756130990520799883880233109311924216460898553921564579
71759211275075089428162122259777421701184691635832087205307269263246072464460134994803723633557042981544362861557257084034088
96.000000

As this is getting beyond my level of know-how I am going to have to put this 
aside for now... I'll take another look later tonight.

Cheers, and thanks for all the pointers

1. code snippet for printing debug info:

printf("[North East]\n");
printf("[%f %f]\n", north, east);

printf("[R C] [V U]\n");
printf("[%i %i] [%i %i]\n", row, col, v, u );

printf("[R0 R1 R2 R3] maprow_f \n");
printf("[%i %i %i %i] %f\n", maprow0, maprow1, maprow2, maprow3, maprow_f );

printf("%f %f %f %f \n", c00, c01, c02, c03);
printf("%f %f %f %f \n", c10, c11, c12, c13);
printf("%f %f %f %f \n", c20, c21, c22, c23);
printf("%f %f %f %f \n", c30, c31, c32, c33);

printf("----------------------------------------\n");
printf("c0 c1 c2 c3   ::  c\n");
printf("%f %f %f %f   ::  %f\n", c0, c1, c2, c3, c);

printf("\n");

-- 
Dylan Beaudette
Soils and Biogeochemistry Graduate Group
University of California at Davis
530.754.7341