up [pdf]
from rsf.proj import *
from rsf.recipes.beg import server

sgy='BOONS.SGY'

Fetch(sgy,'boonsegy',server)

Flow('boons tboons boons.asc boons.bin',
     sgy,
     'segyread tfile=${TARGETS[1]} hfile=${TARGETS[2]} bfile=${TARGETS[3]}')

#C 3 Line number in bytes 9  - 12      105 - 201 (fldr)
#C 4 Trace number in bytes 21 - 24      74 - 206 (cdp)                       

# 110 foot stacking bins

Flow('boon3','boons','intbin xk=cdp yk=fldr | put label3=Line label2=Trace')

Result('boon3','byte gainpanel=all | grey3 frame1=1000 frame2=76 frame3=45 title="Boonsville Data" flat=n point1=0.7 point2=0.7')

Flow('boont','boon3','window j1=2')

Flow('dip3','boont','fdip rect1=20 rect2=5 rect3=5 verb=y')

Flow('dip1','dip3','window n4=1')
Flow('dip2','dip3','window f4=1')

Result('dip1','byte gainpanel=all bar=bar1.rsf | grey3 frame1=500 frame2=76 frame3=45 title="Inline Dip" color=j scalebar=y flat=n point1=0.7 point2=0.7')
Result('dip2','byte gainpanel=all bar=bar2.rsf | grey3 frame1=500 frame2=76 frame3=45 title="Crossline Dip" color=j scalebar=y flat=n point1=0.7 point2=0.7')

Flow('seed','dip3','window n2=1 n3=1 n4=1 | math output=x1')

Flow('shift1','boont','window f2=1')
Flow('shift2','boont','window f3=1')

Flow('last1','boont','window f2=-1 squeeze=n')
Flow('last2','boont','window f3=-1 squeeze=n')

Flow('ref1','shift1 last1','cat axis=2 ${SOURCES[1]}')
Flow('ref2','shift2 last2','cat axis=3 ${SOURCES[1]}')

Flow('ref1s','ref1','add mode=p $SOURCE | stack axis=1 norm=n')
Flow('ref2s','ref2','add mode=p $SOURCE | stack axis=1 norm=n')

Flow('corr1','ref1 boont','add mode=p ${SOURCES[1]} | stack axis=1 norm=n')
Flow('corr2','ref2 boont','add mode=p ${SOURCES[1]} | stack axis=1 norm=n')

Flow('boon2','boont','add mode=p $SOURCE | stack axis=1 norm=n')

Flow('wcos1','corr1 boon2 ref1s',
     '''
     math s1=${SOURCES[1]} s2=${SOURCES[2]} output="(s1*s2)/(input*input)"
     ''')
Flow('wcos2','corr2 boon2 ref2s',
     '''
     math s1=${SOURCES[1]} s2=${SOURCES[2]} output="(s1*s2)/(input*input)"
     ''')
Flow('wcos','wcos1 wcos2',
     '''
     cat axis=3 ${SOURCES[1]} |
     smooth rect1=40 rect2=40 
     ''')

Flow('time','wcos',
     '''
     mul $SOURCE | stack axis=3 norm=n |
     put o1=0 d1=1 o2=0 d2=1 o3=0 d3=1 |
     eikonal vel=n zshot=76 yshot=45
     ''')

Flow('pick','dip3 seed time',
     'pwpaint2 seed=${SOURCES[1]} cost=${SOURCES[2]} order=1')
Result('pick','byte gainpanel=all allpos=y | grey3 frame1=500 frame2=76 frame3=45 title="Relative Age" color=j flat=n point1=0.7 point2=0.7')

Flow('flat','boont pick','iwarp warp=${SOURCES[1]}')
Result('boont','byte gainpanel=all | grey3 frame1=500 frame2=76 frame3=45 wanttitle=n flat=y')
Result('flat','byte gainpanel=all | grey3 frame1=500 frame2=76 frame3=45 wanttitle=n flat=y')

Flow('wboon','boont','window min1=0.75 max1=1.25')
Flow('wflat','flat', 'window min1=0.75 max1=1.25')

Result('wboon','byte gainpanel=all | grey3 frame1=119 frame2=76 frame3=45 wanttitle=n flat=y color=G')
Result('wflat','byte gainpanel=all | grey3 frame1=119 frame2=76 frame3=45 wanttitle=n flat=y color=G')

# Contour3

Flow('cont','boont',
     '''
     window n2=1 n3=1 f2=76 f3=45 | envelope | 
     max1 | window n1=10 | real
     ''')
Plot('wboon','byte gainpanel=all | grey3 frame1=119 frame2=76 frame3=45 wanttitle=n flat=y')
Plot('cont','pick cont',
     '''
     window min1=0.75 max1=1.25 |
     contour3 frame1=119 frame2=76 frame3=45 flat=y 
     wanttitle=n wantaxis=n plotfat=5 cfile=${SOURCES[1]}
     ''')
Result('cont','wboon cont','Overlay')

# Strongest horizon

Flow('cont0','boont',
     '''
     window n2=1 n3=1 f2=76 f3=45 | envelope | 
     max1 | window n1=1 | real
     ''')
Plot('cont0','pick cont0',
     '''
     window min1=0.75 max1=1.25 |
     contour3 frame1=119 frame2=76 frame3=45 flat=y 
     wanttitle=n wantaxis=n plotfat=5 cfile=${SOURCES[1]}
     ''')
Result('cont0','wboon cont0','Overlay')

# cont1 and cont2 -> mask

t0=1.129
dt=0.1

t1=t0-dt
t2=t0+dt

Flow('mask','pick','mask min=%g max=%g | dd type=float | smooth rect1=10' % (t1,t2))
Flow('boonm','boont mask','mul ${SOURCES[1]}')

Plot('wboonm','boonm',
     '''
     window min1=1.05 max1=1.25 |
     byte gainpanel=all | 
     grey3 frame1=50 frame2=76 frame3=45 wanttitle=n flat=y point1=0.25
     ''')
Plot('contm','pick cont0',
     '''
     window min1=1.05 max1=1.25 |
     contour3 frame1=50 frame2=76 frame3=45 flat=y point1=0.25
     wanttitle=n wantaxis=n plotfat=5 cfile=${SOURCES[1]}
     ''')

Plot('boonm',
     '''
     window min1=0.8 max1=1.4 |
     byte gainpanel=all |
     grey3 frame1=%d frame2=76 frame3=45 wanttitle=n 
     flat=y point1=0.7 point2=0.7
     ''' % int((t0-0.8)/0.002))

Result('boonm','wboonm contm','Overlay')

# previous method - similarity with envelope
Flow('rot','boonm','simenv verb=n rect=1000',split=[3,97],reduce='cat axis=4')

# new method - local skewness
Flow('sqr','boonm','localskew verb=n rect=1000',split=[3,97],reduce='cat axis=4')
Flow('kur','boonm','localskew verb=n rect=1000 const=y inv=n',split=[3,97],reduce='cat axis=4')
Flow('skew','sqr kur','mul ${SOURCES[1]}')

Flow('mrot','sqr','stack axis=1 | transp',split=[4,97],reduce='cat axis=3')

Flow('mpik','mrot','pick rect1=10 rect2=10 | add add=90')

Result('mpik',
       '''
       grey color=j scalebar=y title="Estimated Phase Rotation"
       unit1= unit2= transp=n yreverse=n 
       ''')

Flow('ampl','flat','window n1=1 min1=1.129')
Result('ampl','grey title="Amplitude Before Rotation" unit1= unit2= transp=n yreverse=n color=G')

Flow('npik','mpik','spray axis=1 n=1001 d=0.002')

Flow('boonr','boonm npik','phaserot phase=${SOURCES[1]}')

Result('boonr',
       '''
       window min1=0.8 max1=1.4 |
       byte gainpanel=all |
       grey3 frame1=%d frame2=76 frame3=45 wanttitle=n flat=n point1=0.7 point2=0.7
       ''' % int((t0-0.8)/0.002))

Flow('flatr','boonr pick','iwarp warp=${SOURCES[1]}')

Flow('ampr','flatr','window n1=1 min1=1.129')
Result('ampr','grey title="Amplitude After Rotation" unit1= unit2= transp=n yreverse=n color=G')

Flow('hist','ampl','histogram o1=-80000 n1=251 d1=800')

Result('hist','dd type=float | graph title="Amplitude Distribution Before Rotation" ')

Flow('histr','ampr','histogram o1=-80000 n1=251 d1=800')

Result('histr','dd type=float | graph title="Amplitude Distribution After Rotation" ')

End()

sfsegyread
sfintbin
sfput
sfbyte
sfgrey3
sfwindow
sffdip
sfmath
sfcat
sfadd
sfstack
sfsmooth
sfmul
sfeikonal
sfpwpaint2
sfiwarp
sfenvelope
sfmax1
sfreal
sfcontour3
sfmask
sfdd
sfsimenv
sflocalskew
sftransp
sfpick
sfgrey
sfspray
sfphaserot
sfhistogram
sfgraph