Horizon holes filling¶

Author: M. Ravasi, KAUST

Welcome to the Matrix-free inverse problems with PyLops tutorial!

The aim of this tutorial is to fill holes in seismic horizons (also known in the computer vision community as impainting). As a by-product you will learn to:

Familiarize with the pylops.signalprocessing submodule, and more specifically with the 2D-FFT (FFT2D) and Wavelet transform (DWT);
Learn how to use sparse solvers within the pylops.optimization.sparsity submodule, and more specifically the FISTA solver.

From a mathematical point of view we write the impainting problem as follows:

J = arg min_{\mathbf{p}} ||\mathbf{d} - \mathbf{M} \mathbf{P} \mathbf{p}||_2 + \epsilon ||\mathbf{p}||_1

(1)

where $\mathbf{M}$ is a masking operator, $\mathbf{P}$ is a sparsyfing transform, $\mathbf{d}$ is the horizon with holes, and $\mathbf{m}=\mathbf{P}\mathbf{p}$ is the filled horizon we wish to obtain.

Let's first import the libraries we need in this tutorial

# Run this when using Colab (will install the missing libraries)
# !pip install pylops scooby

%matplotlib inline

import numpy as np
import matplotlib.pyplot as plt
import scooby

from pylops.basicoperators import *
from pylops.optimization.sparsity import *
from pylops.signalprocessing import FFT2D, DWT, DWT2D
from pylops.utils.tapers import taper2d

Data loading¶

Let's start by loading a seismic horizon

# You can download this horizon using Madagascar:
# Fetch('horizon.asc','hall')
f = np.loadtxt('data/horizon.asc')
ig = f[:, -1]-65
ig = ig.reshape(291, 196).T
nyorig, nxorig = ig.shape
y = 33.139 + np.arange(nyorig) * 0.01
x = 35.031 + np.arange(nxorig) * 0.01

# We apply some tapering on the edges
#tap = taper2d(nxorig, nyorig, (25, 25))
#ig = ig * tap

# Pad to closest multiple of 2 (for DWT)
ig = np.pad(ig, ((0, 256-nyorig), (0, 512-nxorig)))
ny, nx = ig.shape

plt.figure(figsize=(15, 12))
plt.imshow(ig[:nyorig, :nxorig], cmap='jet', vmin=-14, vmax=14,
           extent=(x[0], x[-1], y[0], y[1]), origin='lower')
plt.xlabel('y (km)')
plt.ylabel('x (km)')
plt.title('Horizon')
plt.axis('tight');

Let's now create some holes in the horizon

mask = np.ones_like(ig)
mask[125:145, 150:170] = 0
mask[150:170, 50:70] = 0
mask[75:95, 175:195] = 0
igholes = ig * mask

plt.figure(figsize=(15, 12))
plt.imshow(igholes[:nyorig, :nxorig], cmap='jet', vmin=-14, vmax=14,
           extent=(x[0], x[-1], y[0], y[1]), origin='lower')
plt.xlabel('y (km)')
plt.ylabel('x (km)')
plt.title('Horizon with holes')
plt.axis('tight');

Reconstruction with DWT¶

First of all we need to create the sparsyfing operator

wavkind = 'sym9'

# Cascaded DWT
Wopy, Wopx = DWT((ny,nx), 1, wavelet=wavkind, level=5), DWT((ny,nx), 0, wavelet=wavkind, level=5)
Wop = Wopy*Wopx
dimswav = Wopy.dimsd

# DWT2D
#Wop = DWT2D((ny, nx), wavelet=wavkind, level=5)
#dimswav = Wop.dimsd

ig_wav = Wop * ig.ravel()
iginv = Wop.H * ig_wav
iginv = iginv.reshape(ny, nx)

fig, axs = plt.subplots(1, 3, figsize=(16, 4))
axs[0].imshow(ig[:nyorig, :nxorig], cmap='jet', vmin=-14, vmax=14,
              extent=(x[0], x[-1], y[0], y[1]), origin='lower')
axs[0].set_title('Image')
axs[0].axis('tight')
axs[1].imshow(ig_wav.reshape(dimswav), cmap='gray', vmin=0, vmax=5e-1)
axs[1].set_title('DWT2 coefficients')
axs[1].axis('tight')
axs[2].imshow(iginv[:nyorig, :nxorig], cmap='jet', vmin=-14, vmax=14,
              extent=(x[0], x[-1], y[0], y[1]), origin='lower')
axs[2].set_title('Reconstructed image')
axs[2].axis('tight');

/Users/ravasim/opt/anaconda3/envs/pylops/lib/python3.8/site-packages/pywt/_multilevel.py:43: UserWarning: Level value of 5 is too high: all coefficients will experience boundary effects.
  warnings.warn(

Let's now solve the inverse problem

Op = Diagonal(mask.ravel())

igdata = Op * ig.ravel()
igfilled_wav, niter, cost = FISTA(Op * Wop.H, igdata, niter=100, 
                                  #threshkind='soft', eps=1e1, 
                                  threshkind='soft-percentile', perc=10, 
                                  returninfo=True, show=True)
igfilled = np.real((Wop.H * igfilled_wav).reshape(ny, nx))

plt.figure(figsize=(12, 3))
plt.plot(cost, 'k', lw=2)
plt.title('Cost function');

FISTA optimization (soft-percentile thresholding)
-----------------------------------------------------------
The Operator Op has 131072 rows and 131072 cols
eps = 1.000000e-01	tol = 1.000000e-10	niter = 100
alpha = 1.000000e+00	perc = 10.0
-----------------------------------------------------------

   Itn       x[0]        r2norm     r12norm     xupdate
     1  -2.13550e+02   3.678e+02   2.794e+03   1.180e+03
     2  -2.13552e+02   3.622e+02   2.783e+03   3.994e+00
     3  -2.13554e+02   3.575e+02   2.773e+03   4.574e+00
     4  -2.13555e+02   3.551e+02   2.763e+03   5.078e+00
     5  -2.13556e+02   3.527e+02   2.753e+03   5.554e+00
     6  -2.13556e+02   3.509e+02   2.743e+03   5.982e+00
     7  -2.13557e+02   3.492e+02   2.733e+03   6.428e+00
     8  -2.13558e+02   3.468e+02   2.723e+03   6.798e+00
     9  -2.13559e+02   3.443e+02   2.713e+03   7.095e+00
    10  -2.13560e+02   3.417e+02   2.703e+03   7.279e+00
    11  -2.13560e+02   3.400e+02   2.695e+03   7.347e+00
    21  -2.13565e+02   3.312e+02   2.628e+03   7.819e+00
    31  -2.13567e+02   3.279e+02   2.598e+03   6.431e+00
    41  -2.13567e+02   3.272e+02   2.585e+03   5.294e+00
    51  -2.13567e+02   3.274e+02   2.584e+03   2.515e+00
    61  -2.13567e+02   3.273e+02   2.584e+03   1.386e+00
    71  -2.13567e+02   3.269e+02   2.583e+03   9.770e-01
    81  -2.13567e+02   3.265e+02   2.583e+03   9.608e-01
    91  -2.13567e+02   3.271e+02   2.583e+03   6.146e-01
    92  -2.13567e+02   3.271e+02   2.583e+03   5.876e-01
    93  -2.13567e+02   3.271e+02   2.583e+03   5.628e-01
    94  -2.13567e+02   3.272e+02   2.583e+03   5.353e-01
    95  -2.13567e+02   3.272e+02   2.583e+03   5.003e-01
    96  -2.13567e+02   3.272e+02   2.583e+03   4.660e-01
    97  -2.13567e+02   3.272e+02   2.583e+03   4.299e-01
    98  -2.13567e+02   3.272e+02   2.583e+03   3.862e-01
    99  -2.13567e+02   3.271e+02   2.583e+03   3.394e-01
   100  -2.13567e+02   3.272e+02   2.583e+03   2.938e-01

Iterations = 100        Total time (s) = 3.21
---------------------------------------------------------

fig, axs = plt.subplots(1, 3, figsize=(16, 4))
axs[0].imshow(ig[:nyorig, :nxorig], cmap='jet', vmin=-14, vmax=14,
              extent=(x[0], x[-1], y[0], y[1]), origin='lower')
axs[0].set_title('Image')
axs[0].axis('tight')
axs[1].imshow(igfilled[:nyorig, :nxorig], cmap='jet', vmin=-14, vmax=14,
              extent=(x[0], x[-1], y[0], y[1]), origin='lower')
axs[1].set_title('Reconstructed image')
axs[1].axis('tight');
axs[2].imshow(ig[:nyorig, :nxorig]-igfilled[:nyorig, :nxorig], cmap='jet', vmin=-14, vmax=14,
              extent=(x[0], x[-1], y[0], y[1]), origin='lower')
axs[2].set_title('Error image')
axs[2].axis('tight');

plt.figure(figsize=(15, 12))
plt.imshow(igfilled[:nyorig, :nxorig], cmap='jet', vmin=-14, vmax=14,
           extent=(x[0], x[-1], y[0], y[1]), origin='lower')
plt.xlabel('y (km)')
plt.ylabel('x (km)')
plt.title('Filled Horizon')
plt.axis('tight');

At this point you may wonder if the DWT is the best sparsifying transform for this data.

I encourage you to experiment with more transforms, for example:

FFT2D (as shown below);
Curvelet transform (see curvelops for a PyLops wrapper);
Other transforms such as the Shearlet transform (you try to wrap pyshearlab into PyLops).

Wop = FFT2D((ny, nx))
ig_wav = Wop * ig.ravel()
igholes_wav = Wop * igholes.ravel()
iginv = Wop.H * ig_wav
iginv = np.real(iginv).reshape(ny, nx)

fig, axs = plt.subplots(1, 3, figsize=(16, 4))
axs[0].imshow(ig[:nyorig, :nxorig], cmap='jet', vmin=-14, vmax=14,
              extent=(x[0], x[-1], y[0], y[1]), origin='lower')
axs[0].set_title('Image')
axs[0].axis('tight')
axs[1].imshow(np.fft.fftshift(np.abs(np.reshape(ig_wav, (ny,nx)))), 
              cmap='gray', vmin=0, vmax=5)
axs[1].set_title('DWT2 coefficients')
axs[1].axis('tight')
axs[2].imshow(iginv[:nyorig, :nxorig], cmap='jet', vmin=-14, vmax=14,
              extent=(x[0], x[-1], y[0], y[1]), origin='lower')
axs[2].set_title('Reconstructed image')
axs[2].axis('tight');

Op = Diagonal(mask.ravel())

igdata = Op * ig.ravel()
igfilled_wav, niter, cost = FISTA(Op * Wop.H, igdata.astype(np.complex), niter=100, 
                                  #threshkind='soft', eps=1e1, 
                                  threshkind='soft-percentile', perc=10, 
                                  returninfo=True, show=True)
igfilled = np.real((Wop.H * igfilled_wav).reshape(ny, nx))

plt.figure(figsize=(12, 3))
plt.plot(cost, 'k', lw=2)
plt.title('Cost function');

FISTA optimization (soft-percentile thresholding)
-----------------------------------------------------------
The Operator Op has 131072 rows and 131072 cols
eps = 1.000000e-01	tol = 1.000000e-10	niter = 100
alpha = 1.000000e+00	perc = 10.0
-----------------------------------------------------------

   Itn       x[0]        r2norm     r12norm     xupdate
     1   9.37598e+00   3.312e+03   6.465e+03   1.165e+03
     2   9.72815e+00   2.777e+03   5.926e+03   2.491e+01
     3   1.00406e+01   2.461e+03   5.607e+03   2.177e+01
     4   1.03056e+01   2.307e+03   5.442e+03   1.914e+01
     5   1.05361e+01   2.224e+03   5.346e+03   1.724e+01
     6   1.07374e+01   2.179e+03   5.289e+03   1.577e+01
     7   1.09111e+01   2.161e+03   5.259e+03   1.457e+01
     8   1.10651e+01   2.148e+03   5.235e+03   1.364e+01
     9   1.12026e+01   2.132e+03   5.213e+03   1.286e+01

/Users/ravasim/Desktop/KAUST/OpenSource/pylops/pylops/optimization/sparsity.py:1099: ComplexWarning: Casting complex values to real discards the imaginary part
  msg = '%6g  %12.5e  %10.3e   %9.3e  %10.3e' % \

    10   1.13216e+01   2.129e+03   5.202e+03   1.217e+01
    11   1.14275e+01   2.124e+03   5.192e+03   1.159e+01
    21   1.21523e+01   2.117e+03   5.165e+03   4.671e+00
    31   1.25355e+01   2.119e+03   5.166e+03   1.236e+00
    41   1.23711e+01   2.120e+03   5.166e+03   1.134e+00
    51   1.21924e+01   2.119e+03   5.166e+03   5.629e-01
    61   1.22937e+01   2.118e+03   5.165e+03   3.963e-01
    71   1.23834e+01   2.118e+03   5.165e+03   2.951e-01
    81   1.23040e+01   2.119e+03   5.166e+03   1.610e-01
    91   1.22670e+01   2.119e+03   5.165e+03   1.657e-01
    92   1.22706e+01   2.119e+03   5.165e+03   1.698e-01
    93   1.22752e+01   2.119e+03   5.165e+03   1.698e-01
    94   1.22807e+01   2.119e+03   5.165e+03   1.660e-01
    95   1.22868e+01   2.119e+03   5.165e+03   1.585e-01
    96   1.22933e+01   2.118e+03   5.165e+03   1.479e-01
    97   1.23000e+01   2.118e+03   5.165e+03   1.346e-01
    98   1.23067e+01   2.118e+03   5.165e+03   1.193e-01
    99   1.23132e+01   2.118e+03   5.165e+03   1.029e-01
   100   1.23193e+01   2.118e+03   5.165e+03   8.630e-02

Iterations = 100        Total time (s) = 2.42
---------------------------------------------------------

fig, axs = plt.subplots(1, 3, figsize=(16, 4))
axs[0].imshow(ig[:nyorig, :nxorig], cmap='jet', vmin=-14, vmax=14,
              extent=(x[0], x[-1], y[0], y[1]), origin='lower')
axs[0].set_title('Image')
axs[0].axis('tight')
axs[1].imshow(igfilled[:nyorig, :nxorig], cmap='jet', vmin=-14, vmax=14,
              extent=(x[0], x[-1], y[0], y[1]), origin='lower')
axs[1].set_title('Reconstructed image')
axs[1].axis('tight');
axs[2].imshow(ig[:nyorig, :nxorig]-igfilled[:nyorig, :nxorig], cmap='jet', vmin=-14, vmax=14,
              extent=(x[0], x[-1], y[0], y[1]), origin='lower')
axs[2].set_title('Error image')
axs[2].axis('tight');

scooby.Report(core='pylops')

Horizons

Horizon holes filling¶

Data loading¶

Reconstruction with DWT¶