Geophysical Computing Notebook
 
Joseph Jennings

Table Of Contents

This Page

Structure-oriented smoothing of 2D & 3D images on the GPU

March 20, 2015

This page describes the code I wrote for my undergraduate thesis while in my final year at Colorado School of Mines. The code is found in two different git repositories on my GitHub page. The final working code and demo is found within my fork of the Mines Java Tool Kit created by Professor Dave Hale (my advisor for this project). This is the jtk repository and my test and other code is found within the Senior_Design repository.

Installation

The instructions for installing the working code and demo can be found in the README file on this web page.

Introduction

Geoscientists use seismic images to better understand the subsurface geology of a given region. The proper interpretation of these images is key in creating an accurate geologic model to make predictions on the location of hydrocarbon reservoirs. This interpretation can be time-consuming and a major bottleneck in the process of creating geologic models [Fehmers-Hocker]. Applying image processing algorithms can aid in speeding up this interpretation as they can remove unwanted noise and help delineate subtle yet important geologic features. One of these algorithms is known as structure-oriented smoothing [Fehmers-Hocker]. The application of this filter to a seismic image will smooth along the geologic structure in the image rather than across it in order to enhance the structural features. By clicking Figure 1 below, you can see the results of structure-oriented smoothing of a 2D seismic image.

Figure 1: Click the figure above to toggle between the input and smoothed image.

Structure-oriented smoothing

The underlying algorithm that performs this adaptive smoothing is known as anisotropic diffusion filtering and has been used extensively in image processing applications [Weickert]. By solving the anisotropic diffusion equation using the input seismic image as an initial condition, we simulate the anistropic diffusion process which diffuses the seismic amplitudes of the input image [Fehmers-Hocker]. This diffusion is guided along the structure present in the image by computed structure tensors that serve as input parameters in the anisotropic diffusion process. [Hale], describes an implementation of this algorithm by implicitly solving a modified version of the diffusion equation via conjugate gradient iterations.

The software that implements this conjugate gradient method for solving the anisotropic diffsuion equation can be found within the solve method in the LocalSmoothingFilter class of the Mines JTK. This method, written by Hale, is shown in the listing below:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
  private void solve(Operator2 a, float[][] b, float[][] x) {
    int n1 = b[0].length;
    int n2 = b.length;
    float[][] d = new float[n2][n1];
    float[][] q = new float[n2][n1];
    float[][] r = new float[n2][n1];
    scopy(b,r);
    a.apply(x,q);
    saxpy(-1.0f,q,r); // r = b-Ax
    scopy(r,d); // d = r
    float delta = sdot(r,r); // delta = r'r
    float bnorm = sqrt(sdot(b,b));
    float rnorm = sqrt(delta);
    float rnormBegin = rnorm;
    float rnormSmall = bnorm*_small;
    int iter;
    log.fine("solve: bnorm="+bnorm+" rnorm="+rnorm);
    for (iter=0; iter<_niter && rnorm>rnormSmall; ++iter) {
      log.finer("  iter="+iter+" rnorm="+rnorm+" ratio="+rnorm/rnormBegin);
      a.apply(d,q); // q = Ad
      float dq = sdot(d,q); // d'q = d'Ad
      float alpha = delta/dq; // alpha = r'r/d'Ad
      saxpy( alpha,d,x); // x = x+alpha*d
      saxpy(-alpha,q,r); // r = r-alpha*q
      float deltaOld = delta;
      delta = sdot(r,r); // delta = r'r
      float beta = delta/deltaOld;
      sxpay(beta,r,d); // d = r+beta*d
      rnorm = sqrt(delta);
    }
    log.fine("  iter="+iter+" rnorm="+rnorm+" ratio="+rnorm/rnormBegin);
  }

GPU Implementation