[GSIF] what API to call for a simple reboxing test?

summary: I'm trying to test reboxing, or 3D interpolation, using
package=GSIF. The idea is to

* assume a uniform input field

* create a uniform input grid with box/voxel dimensions=2x2x2 with
  each box having value=1

* "rebox" to a uniform output grid with dimensions=1x1x8: if output
  box values approximate input box values, declare test to pass.

But I'm not seeing what GSIF API to call to calculate the output box
values, and my geostatistical background (or lack thereof) does not
prepare me for

http://gsif.r-forge.r-project.org/tutorial_eberg.php

details:

# 1. Create 3D input grid and data such that each of its regular boxes
#    has dims=2x2x2 and value=1

# horizontal grid sequences to produce input 2D grid unit area=4
hor.from.in <- 0
hor.to.in   <- 8
hor.by.in   <- 2
# grid corners
corners.hor.in <- seq(from=hor.from.in, to=hor.to.in, by=hor.by.in)
#corners.hor.in
# grid midpoints
midpts.hor.in <-
  seq(from=(hor.from.in+(hor.by.in/2)),
      to=(hor.to.in-(hor.by.in/2)),
      by=hor.by.in)
#midpts.hor.in

# vertical grid sequence to produce input 3D grid unit volume=8
ver.from.in <- 0
ver.to.in   <- 8
ver.by.in   <- 2
corners.ver.in <- seq(from=ver.from.in, to=ver.to.in, by=ver.by.in)
#corners.ver.in
midpts.ver.in <-
  seq(from=(ver.from.in+(ver.by.in/2)),
      to=(ver.to.in-(ver.by.in/2)),
      by=ver.by.in)
#midpts.ver.in

grid3D.in <-
  expand.grid(x=midpts.hor.in, y=midpts.hor.in, z=midpts.ver.in)
# examined below

val.in <- 1.0 # data value for each input box of volume=8
data3D.in <-
  data.frame(v=rep(val.in, times=length(row.names(grid3D.in))))
coordinates(data3D.in) <- grid3D.in
summary(data3D.in)
#> Object of class SpatialPointsDataFrame
#> Coordinates:
#>   min max
#> x   1   7
#> y   1   7
#> z   1   7
#> Is projected: NA 
#> proj4string : [NA]
#> Number of points: 64
#> Data attributes:
#>   Min. 1st Qu.  Median  Mean 3rd Qu.  Max. 
#>      1       1       1     1       1     1 

gridded(grid3D.in)=~x+y+z
summary(grid3D.in)
#> Object of class SpatialPixels
#> Coordinates:
#>   min max
#> x   0   8
#> y   0   8
#> z   0   8
#> Is projected: NA 
#> proj4string : [NA]
#> Number of points: 64
#> Grid attributes:
#>   cellcentre.offset cellsize cells.dim
#> x                 1        2         4
#> y                 1        2         4
#> z                 1        2         4

# So there are 64 box/voxels each of volume=8 in the input grid
# and each box has value=1.

# 2. Create 3D output grid such that each of its regular boxes
#    has dims=1x1x8

# Previous grid setup method fails, so instead do

grid3D.out <-
  SpatialGrid(
    GridTopology(
      c(0.5,0.5,4.0), # centers
      c(1.0,1.0,8.0), # sizes
      c(8,  8,  1)),  # dim/number
    proj4string = CRS(as.character(NA)))
summary(grid3D.out)
# gets
#> Object of class SpatialGrid
#> Coordinates:
#>      min max
#> [1,]   0   8
#> [2,]   0   8
#> [3,]   0   8
#> Is projected: NA 
#> proj4string : [NA]
#> Grid attributes:
#>   cellcentre.offset cellsize cells.dim
#> 1               0.5        1         8
#> 2               0.5        1         8
#> 3               4.0        8         1

# So there are 64 box/voxels each of volume=8 in the output grid.

Since the input and output spaces are uniform, and each contains total
values=64, I'm assuming I should be able to calculate (aka predict) the
values in each output box, and furthermore that the output box values
should approximate the input box values. Am I missing something? If not:

How to call GSIF API to compute the values of the output boxes from the
input grid and its values?

Your assistance is appreciated, Tom Roche <Tom_Roche at pobox.com>