Applying a linear operator to IGF90 objects

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.

 

anis
Figure 1 The left panel shows the Seabeam dataset after being interpolated into a $100 \times 100$ mesh. In the right panel, the interpolation result is convolved with a Laplacian filter. The Seabeam dataset has been decimated to 1/100th the original dataset.

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;