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
  use adj_mod

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

function bin1_op (adj, add, mod, ord) result (stat)
    integer              :: stat
    logical, intent (in) :: adj, add
    real, dimension (:)  :: mod, ord

integer              :: i

call adjnull (adj, add, mod, ord)

do i = 1, nd ; if (m (i)) cycle
       if (adj) then
          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.