We used the two-stage data infill method Claerbout (1996). Our
results are similar to Claerbout's; the interpolated data shows a loss
of power near the edge of the rectangle (see Figure 1).
Claerbout's Ratfor implementation uses two routines, `pef` and
`mis`, which solve for the filter and the data, respectively. The
abstraction of C++ and the HCL framework allows a single routine to
be used for both the data and the filter:

// Solve the missing data problem LMx = -Lk, that is L(Mx + k) = 0 // k contains the known data (and is zero where unknown) // M is a mask operator which zeroes the known components // MissSolve can be used for finding missing data or an unknown filtervoid MissSolve( HCL_LinearOpAdj &L, const HCL_Vector &k, // known data with zeroes where unknown HCL_Vector &x, // initial guess (should be pre-masked!) // on exit x contains u, the unknown data HCL_LinearOpAdj &M, // mask which zeroes the known data HCL_LinearSolver &S ) // solver used with the normal equations { HCL_Vector *b = L.Range().Member(); // Allocate space for b L.Image( k, *b ); // Set b = Lk (*b).Neg(); // Set b = -Lk HCL_CompLinearOpAdj A( &M, 0, &L, 0 ); // Set A = LM // Solve LMx = -Lk // (using the normal equations) NormalSolve( A, *b, x, S ); delete b; }

The first step of two-stage data infill is to create two RGFs which
consist of zeroes and ones. The first has the shape of the data
region, but with ones where the data is known. The second has the
shape of the prediction-error filter, but with ones where the filter
values are modifiable. Internal convolution is then performed with
these two RGFs to give an RGF which we call the influence mask. Then
we use `MissSolve` to find the prediction-error filter using a
linear operator which is the composition of two operators.
Specifically, it consists of `icaf` with the known data, followed
by masking with the influence mask. (Notice that `icaf` is an
operator which takes a filter as the argument).

After finding the filter, we fill in the missing data by again using
`MissSolve`. This time, the operator is `tcai` with the
newly-calculated filter. We employed HCL's conjugate gradient solver
to find both the filter and the missing data.

11/11/1997