from rsf.proj import *
from rsf.recipes.helderiv import Helderiv

# Download data

horizon = 'Penobscot_HorB.txt'

Fetch(horizon,'data',
      server='https://raw.githubusercontent.com',
      top='agile-geoscience/notebooks/master')

# Convert from xyz triples to horizon

Flow('xyz',horizon,
     '''
     echo n1=3 n2=254674 data_format=ascii_float in=$SOURCE |
     dd form=native
     ''')

Flow('xy','xyz','window n1=2 | dd type=int')
Flow('horizon','xyz xy',
     '''
     window f1=2 squeeze=n | 
     intbin head=${SOURCES[1]} xkey=0 ykey=1 |
     window | put label=Time unit=ms label1=Inline label2=Crossline
     ''')

Result('horizon','grey color=j bias=750 scalebar=y title=Horizon yreverse=n transp=n')

# Interpolate missing values

Flow('filled','horizon','lapfill grad=y verb=y niter=500')

# Display with the cube helix colormap and the correct aspect ratio

Plot('cubeh','filled',
     '''
     grey color=x bias=950 clip=133 scalebar=y title=Horizon 
     yreverse=n transp=n screenwd=11.9 screenht=9.62
     ''')

Result('cubeh','Overlay')

# Take a window

Flow('window','horizon','window f1=200 n1=125 f2=200 n2=100')

Result('window',
       'grey color=x bias=950 clip=133 scalebar=y title=Horizon yreverse=n transp=n screenratio=0.8')

Flow('wind.asc',None,
     '''
     echo 
     200 200 
     200 300
     325 300
     325 200
     200 200
     n1=2 n2=5 data_format=ascii_int in=$TARGET
     ''')
Plot('wind','wind.asc',
     '''
     dd type=complex form=native | window |
     graph wanttitle=n wantaxis=n plotcol=7 plotfat=10 
     min1=1 max1=595 min2=1 max2=463 scalebar=y screenwd=11.9 screenht=9.62
     ''')

Result('wind','cubeh wind','Overlay')

# Add noise
Flow('noisy','filled','noise type=n seed=2014 range=15')

# Regular smoothing
Flow('smoothed','noisy','smooth rect1=4 rect2=4')

# Anisotropic diffusion
Flow('diffused','noisy','impl2 rect1=10 rect2=10 tau=1')

# Median filter
Flow('median','noisy','despike2 wide1=10 wide2=10')

# Bilateral filter
Flow('bilateral','noisy','bilat2 a1=0.1 a3=0.001 r1=10 r2=10')

for case in Split('filled noisy smoothed diffused median bilateral'):
    Result(case,
           '''
           grey color=j bias=950 clip=133 scalebar=y 
           title="%s Horizon" yreverse=n transp=n
           minval=800 maxval=1100
           ''' % case.capitalize())
    if case != 'filled':
        Result(case+'-slice',['filled',case],
               '''
               cat axis=3 ${SOURCES[1]} |
               window n1=1 f1=300 n2=450 |
               graph plotfat=7,1 wanttitle=n 
               yreverse=y label2=Time unit2=ms
               ''')

# Shaded relief
Flow('dx','filled','deriv')
Flow('dy','filled','transp | deriv | transp')
Flow('shaded','dx dy','add ${SOURCES[1]}')

# Helix derivative
Helderiv('hel',na=10)

Flow('helix','filled hel','helicon filt=${SOURCES[1]}')

for case in Split('shaded helix'):
      Result(case,
           '''
           grey title="%s Horizon" yreverse=n transp=n mean=y
           ''' % case.capitalize())

End()