This section illustrates how to use the C++ interface to apply an operator to an IGF90 data object. Operators, like the IGF90 class, are also derived from a class in the HCL library. This ensures the correct interaction between operators and data objects. The implementation of the convolution operator is encapsulated in the IGF90_ConvOp class. The subscript IGF90_ reminds us that this operator is to be applied to IGF90 objects. The IGF90_ConvOp class is derived from the HCL_LinearOpAdj (Linear Operators with an Adjoint) class. The details of how to build the interface of an operator will be presented in a later section. The example below shows the invocation of the forward transform of a convolution operator:

#include "iostream.h" #include "IGF90.h" #include "IGF90_ConvOp.h" ... // Create model vector and a pointer to its vector space IGF90 x( "model" ); // read from disk IGF90Space *domain_space = &( x.Space() ); // pointer to header info IGF90Space *range_space = domain_space; // Create output vector IGF90 y( *range_space ); y.Random(); // Create the convolution operator float eps; if(!fetch("eps","f",&eps)) eps=1. ; IGF90 filter("filt"); // scale filter filter.Mul(eps); IGF90_ConvOp C( domain_space , range_space, filter ) ; // Try out operator C.Image( x, y ); // y = C x y.Write(); // output result

The syntax `IGF90_ConvOp C(domain_space, range_space, filter)`

constructs the convolution operator `C`

, given a pointer to the
domain and range vector spaces and the filter coefficients. The
operator `C`

is an instance (an object) of the class
IGF90_ConvOp. It can only be applied to IGF90 objects that have the
same vector space as domain_space, and the output can only be saved
in an IGF90 object that has the same vector space as
range_space. This ensures that the geometry and the grid information
of the input and output data sets correspond to what the operator
expects. In the special case of the convolution operator, the input
and output data sets have regular geometry.

The syntax `C.Image( x, y )`

applies the forward mapping of the
convolution operator `C`

to the input `x`

and outputs the
result in `y`

. Before applying the mapping, the `Image()`

member function checks that the input and output IGF90 objects being
passed agree with the domain and range spaces being handed over to the
`C`

operator at construction time.

If the adjoint is to be performed, the `Adjoint()`

member
function should be called on the `C`

operator before the
`Image()`

member function is called. Note that the `Image()`

member function expects as first argument the input and as second
argument the output, independent of the operator being the forward or
adjoint mapping. The following lines of code added to the above
example call the forward and adjoint operations of the operator. The
example implements a **dot product test** on the convolution
operator:

IGF90 xp( *domain_space ); IGF90 yp( *range_space ); y.Random(); x.Random(); // Applied the Forward mapping C.Image( x, yp ); // yp = C x // Applied the Adjoint mapping C.Adjoint().Image( y, xp ); // xp = C' y // Dot product for convolution operator cerr << " " << yp.Inner(y)<< " " << xp.Inner(x) << endl;

Trying out another operator consists of changing a few lines of
code. Assuming we keep the IGF90 model vector used in the previous
example, `x`

, we can write a new routine using a linear
interpolation operator:

#include "iostream.h" #include "IGF90.h" #include "IGF90_ConvOp.h" ... IGF90 data( "data" ); IGF90Space *Dspace = &( (IGF90Space&) data.Space() ); // Create linear interpolation operator, // interpolate along header keys "sx" and "sy". IGF90_LinInterpOp L( domain_space, Dspace, "sx", "sy" ); L.Adjoint().Image( data, x ); // x = L' data x.Write();The implementation of the linear interpolation operator is embedded within the IGF90_LinInterpOp class, also derived from HCL_LinearOpAdj. The constructor requires passing the domain and range in which the object will operate. The other two arguments indicate along which headers the interpolation will take place.

We can use the above example to interpolate the Seabeam data set
Claerbout (1994). The data set is a collection of depth
measurements at miscellaneous locations in the (*x*,*y*) plane, which
are represented in a SEPlib90 structure as a set of one-sample-long
traces. We use the IGF90 `data`

object to hold this irregular
data set. The left panel of Figure 1 shows the result of
calling the adjoint interpolation operator on the data set. The
gridding size is determined by the model vector `x`

. The right
panel in Figure 1 shows the result of convolving the
interpolated Seabeam data set with a Laplacian filter. In the next
subsection, we will see how to use the HCL library to chain these two
operators into one.

Figure 1

The dot product test for the linear interpolation operator is simply coded as

IGF90 dp( *Dspace ); IGF90 xp( *domain_space ); data.Random(); x.Random(); // Applied the Forward mapping L.Image( x, dp ); // dp = L x // Applied the Adjoint mapping L.Adjoint().Image( data, xp ); // xp = L' data // Dot product for linear interpolation operator cerr << " " << dp.Inner(data)<< " " << xp.Inner(x) << endl;

11/11/1997