Homework 3 |
hole
Figure 8. Seismic depth slice after removing selected parts of the data. | |
---|---|
The goal of the next exercise is to figure out if one can use compactness of the Fourier transform to reconstruct missing data. The missing parts are created artificially by cutting holes in the original data (Figure 8).
fft-data,fft-hole
Figure 9. Fourier transform of the original data (a) and data with holes with holes (b). The absolute value is displayed |
---|
fft-mask
Figure 10. Fourier-domain mask for selecting a convex set. | |
---|---|
Figures 9a and 9b show the digital Fourier transform of the original data and the data with holes. We observe again that the support of the data in the Fourier domain is compact thanks to the data smoothness. Cutting holes in the physical domain creates discontinuities that make the Fourier response spread beyond the original support. Figure 10 shows a Fourier-domain mask designed to contain the support of the original data.
To accomplish the task of missing data interpolation, we will use an iterative method known as POCS (projection onto convex sets). By definition, a convex set is a set of functions such that, for any and from the set, (for ) also belongs to the set. A projection onto a convex set means finding a function in the set that is of the shortest distance to the given function. The POCS theorem states that if one wants to find a function that belongs to the intersection of two convex sets and , the task can be accomplished iteratively by alternating projections onto the two sets (Youla and Webb, 1982).
In our example, is the set of all functions that are equal to the known data outside of the holes. is the set of all functions that have a predifined compact support in the Fourier domain (and therefore are smooth in the physical domain). The algorithm consists of the following steps:
pocs
Figure 11. Missing data interpolated by iterative projection onto convex sets. | |
---|---|
from rsf.proj import * # Download data Fetch('horizon.asc','hall') # Convert format Flow('data','horizon.asc', ''' echo in=$SOURCE data_format=ascii_float n1=3 n2=57036 | dd form=native | window n1=1 f1=-1 | add add=-65 | put n2=291 o2=35.031 d2=0.01 label2=y unit2=km n1=196 o1=33.139 d1=0.01 label1=x unit1=km | costaper nw1=25 nw2=25 ''') # Display def plot(title): return ''' grey color=j title="%s" transp=y yreverse=n clip=14 ''' % title Result('data',plot('Horizon')) # Cut three square holes (!!! CHANGE ME !!!) cut = ''' cut n1=20 n2=20 f1=125 f2=150 | cut n1=20 n2=20 f1=150 f2=50 | cut n1=20 n2=20 f1=75 f2=175 ''' Flow('hole','data',cut) Flow('mask','data', 'math output=1 | %s | math output=1-input' % cut) Plot('hole',plot('Input')) Result('hole','Overlay') # Fourier transform forw = 'rtoc | fft3 axis=1 pad=1 | fft3 axis=2 pad=1' back = 'fft3 axis=2 inv=y | fft3 axis=1 inv=y | real' for data in ('data','hole'): fft = 'fft-'+data Flow(fft,data,forw) Result(fft, ''' math output="abs(input)" | real | grey allpos=y title="Fourier Transform (%s)" screenratio=1 ''' % data.capitalize()) # Create Fourier mask Flow('fft-mask','fft-hole', ''' real | math output="x1*x1+x2*x2" | mask min=50 | dd type=float | math output=1-input | smooth rect1=5 rect2=5 repeat=3 | rtoc ''') Result('fft-mask', ''' real | grey allpos=y title="Fourier Mask" screenratio=1 ''') # POCS iterations niter=5 # !!! CHANGE ME !!! data = 'hole' plots = ['hole'] for iter in range(niter): old = data data = 'data%d' % iter # 1. Forward FFT # 2. Multiply by Fourier mask # 3. Inverse FFT # 4. Multiply by space mask # 5. Add data outside of hole Flow(data,[old,'fft-mask','mask','hole'], ''' %s | mul ${SOURCES[1]} | %s | mul ${SOURCES[2]} | add ${SOURCES[3]} ''' % (forw,back)) Plot(data,plot('Iteration %d' % (iter+1))) plots.append(data) # Put frames in a movie Plot('pocs',plots,'Movie',view=1) # Last frame Result('pocs',data, plot('POCS interpolation (%d iterations)' % niter)) End() |
Your task:
scons viewto reproduce the figures on your screen.
scons pocs.vplto see a movie of different iterations.
Homework 3 |