Next: Linear interpolation Up: EXAMPLES OF LINEAR OPERATORS Previous: Causal and leaky integration

## Nearest-neighbor interpolation

Binning irregularly sampled data to a regular grid is another important operator, often found in practical applications. The adjoint of binning is interpolation, and the simplest case of interpolation is nearest-neighbor interpolation, where each data point is assigned with the nearest grid point value.

Irregular data typically come in ordinate-coordinate pairs. According to the minimal interface design, coordinates as well as the grid information belong to the initialization interface, while the corresponding ordinates go directly to the operator interface. At the initialization step, module bin1 transforms the grid and coordinate information into two data-size arrays: the integer array x with the nearest bin locations and the logical array m, which masks the data points outside of the grid. Additionally, there is a real array w with the weighting function, optionally supplied by the user.

```module bin1

integer,                             private :: nd
integer, dimension (:), allocatable, private :: x
real,    dimension (:), allocatable, private :: w
logical, dimension (:), allocatable, private :: m

contains

subroutine bin1_init (coord, o1, d1, n1, weight)
real, dimension (:), intent (in) :: coord, weight
real,                intent (in) :: o1, d1
integer,             intent (in) :: n1
optional                         :: weight

nd = size (coord)
if (.not. allocated (x)) allocate (x (nd), m (nd), w (nd))

x = 1.5 + (coord - o1) / d1
m = (x < 1 .or. x > n1)

if (present (weight)) then
m = m .or. weight == 0.
where (.not. m) w = weight
else
where (.not. m) w = 1.
end if
end subroutine bin1_init

integer              :: stat
real, dimension (:)  :: mod, ord

integer              :: i

do i = 1, nd ; if (m (i)) cycle
mod (x (i)) = mod (x (i)) + w (i) * ord (   i )
else
ord (   i ) = ord (   i ) + w (i) * mod (x (i))
end if
end do

stat = 0
end function bin1_op

subroutine bin1_close ()
deallocate (x, m, w)
end subroutine bin1_close

end module bin1
```

The output of the interpolation operator is ``smooth'', each data point obtains a unique value from the nearest bin. In a shorthand Fortran 90 notation, we could write the forward operation in one line as where (.not. m) ord = ord + w * mod (x). However, writing the adjoint operation as where (.not. m) mod (x) = mod (x) + w * ord is not correct, because the output of the adjoint process (binning) is ``rough'': each bin may receive zero, one, or more than one input from the data.

Next: Linear interpolation Up: EXAMPLES OF LINEAR OPERATORS Previous: Causal and leaky integration
Stanford Exploration Project
11/11/1997