from rsf.proj import *

g = 'grey crowd1=0.96 crowd2=0.85 wantaxis=n title="%s" '

# Seismic horizon
Fetch('horizon.asc','hall')
Flow('horizon','horizon.asc',
     '''
     echo in=$SOURCE data_format=ascii_float n1=3 n2=57036 | 
     dd form=native | window n1=1 f1=-1 | 
     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
     ''')

# CT-scan slice
Fetch('slice.rsf','ctscan')
Flow('circle','slice','dd type=float')
Flow('square','circle',
     'window window n1=366 n2=366 f1=73 f2=73')

# Satellite data
Fetch('mod.dat','ying')
Flow('sat','mod.dat',
     '''
     echo in=$SOURCE n1=1200 n2=1200 n3=7
     data_format=native_int |
     put d1=1 d2=1 o1=0 o2=0 d3=1 o3=0 |
     window n3=1 f3=3 | dd type=float |
     scale axis=2
     ''', stdin=0)

# Your own data
Fetch('your.dat','ying')
Flow('your','your.dat','dd type=float')

textures = ('horizon','square','sat','your')

for txt in textures:
    # Remove trend
    #-------------
    m = txt+'-one'
    x = txt+'-x'
    y = txt+'-y'
    Flow(m,txt,'math output=1')
    Flow(x,txt,'math output=x1 | scale axis=2')
    Flow(y,txt,'math output=x2 | scale axis=2')
    flt = txt+'-flt'
    trd = txt+'-trd'
    Flow([flt,trd],[m,x,y,txt],
         '''
         cat ${SOURCES[1:3]} | 
         lpf match=${SOURCES[3]} pred=${TARGETS[1]}
         rect1=1000 rect2=1000 
         ''')
    t = txt+'-dtrd'
    Flow(t,[txt,trd],'add scale=1,-1 ${SOURCES[1]}')

    # Original
    #---------
    Plot(t,g % '(a) Training Image')

    # Estimate PEF
    #-------------
    pef = txt+'-pef'
    lag = txt+'-lag'
    Flow([pef,lag],t,
         'hpef niter=50 a=10,10 lag=${TARGETS[1]}')

    # PEF residual
    #-------------
    wht = txt+'-wht'
    Flow(wht,[t,pef],'helicon filt=${SOURCES[1]}')
    Plot(wht,g % '(b) Training Image * PEF')

    # Reconstruct original
    #---------------------
    rec = txt+'-rec'
    Flow(rec,[wht,pef],'helicon filt=${SOURCES[1]} div=y')
    Plot(rec,g % '(c) Training Image * PEF / PEF')

    # Synthesized image
    #------------------
    syn = txt+'-syn'
    Flow(syn,[t,pef],
              '''
              noise rep=y seed=2014 | 
              helicon filt=${SOURCES[1]} div=y
              ''')
    Plot(syn,g % '(d) Random Numbers / PEF')

    Result(txt,[t,wht,rec,syn],'TwoRows')

End()