up [pdf]
from rsf.proj import *
from rsf.prog import RSFROOT

def Grey(data,other):
        Result(data,'grey labelsz=10 titlesz=10 clip=8501 color=g label1=Time unit1=s %s'%other)

Flow('gath','elfgath','dd form=native | put unit1=s unit2=m label2=Offset d2=25 o2=100')
Plot('gath','grey labelsz=12 titlesz=12.5 title="(a) Input" ')

Flow('mask','gath','math output=1 | pad n1=1000')
Flow('gpad','gath','pad n1=1000')

Flow('dip','gpad mask','dip rect1=10 rect2=10 p0=0 pmin=0 mask=${SOURCES[1]} order=2')
Plot('dip','window n1=800 | grey labelsz=12 titlesz=12.5 color=j title="(b) Slope" allpos=y')

Result('dip','window n1=800 | grey labelsz=10 titlesz=10 color=j label1=Time unit1=s title="Local slope" scalebar=y screenratio=1.2')

Flow('trace','dip','window n2=1 | math output=x1')
Flow('pick','dip trace','pwpaint seed=${SOURCES[1]} eps=0.05 order=2')
Plot('pick0','pick',
     '''
     window n1=800 |
     grey labelsz=12 titlesz=12.5 color=j allpos=y
     title="(d) Relative Age" clip=3.2
     ''')
Plot('pick1','pick',
     '''
     window n1=800 |
     contour c0=0 dc=0.05 nc=64
     transp=y yreverse=y plotcol=7 plotfat=5
     wanttitle=n wantaxis=n
     ''')
Plot('pick','pick0 pick1','Overlay')

Flow('flat','gpad pick','iwarp warp=${SOURCES[1]} n1=800')
Flow('gath-rec','flat pick','iwarp warp=${SOURCES[1]} inv=n | window n1=800')

Plot('flat','grey labelsz=12 titlesz=12.5 title="(d) Flattened" ')
Plot('gath-rec','grey labelsz=12 titlesz=12.5 title="EMD filtered" ')

Flow('cont','gpad',
     'window n2=1 | envelope | max1 | window n1=20 | real')

Flow('k1.p','cont',
     '''
     math output="1.5+input/0.004" | 
     dd type=int form=ascii format="%d," line=100 --out=$TARGET
     ''',stdout=0)
Command('k1.par','k1.p',
         'printf "k1=" > $TARGET && cat $SOURCE >> $TARGET')

Flow('spikes','k1.par',
     '''
     spike n1=1000 mag=6000 
     nsp=20 par=$SOURCE |
     smooth rect1=2 repeat=2 
     ''',stdin=0)

Flow('paint','dip spikes','pwpaint seed=${SOURCES[1]} eps=0.05 order=2')
Plot('paint','paint gath',
     '''
     window n1=800 | mutter v0=1200 half=n |
     add ${SOURCES[1]} scale=12,1 | 
     grey labelsz=12 titlesz=12.5 color=G title="(c) Painted" allpos=y
     ''')

Result('pgath','gath dip paint flat gath-rec','SideBySideAniso')

########################################################################
# INITIALIZATION
########################################################################
matlab         = WhereIs('matlab')
matROOT = '../matfun'
matfun = 'FX_EMD'
matlabpath = os.environ.get('MATLABPATH',os.path.join(RSFROOT,'lib'))

if not matlab:
    sys.stderr.write('\nCannot find Matlab.\n')
    sys.exit(1)
n1=800
n2=128
dt=0.004
lf=5
hf=120
N=1
verb=0
############################################################
## with parameter
############################################################
Flow('flat-emd0',[os.path.join(matROOT,matfun+'.m'),'flat'],
     '''MATLABPATH=%(matlabpath)s %(matlab)s 
     -nosplash -nojvm -r "addpath %(matROOT)s;%(matfun)s('${SOURCES[1]}','${TARGETS[0]}',%(n1)d,%(n2)d,%(dt)g,%(lf)g,%(hf)g,%(N)d,%(verb)d);quit"
     '''%vars(),stdin=0,stdout=-1)
Flow('flat-emd','flat-emd0','put d2=25 d1=0.004 o2=100 o1=0 label2=Offset unit2=m')
Flow('flat-mf','flat','transp | mean nfw=10 | transp')

Flow('flat-emd-rec0','flat-emd pick','iwarp warp=${SOURCES[1]} inv=n | window n1=800')
Flow('flat-mf-rec0','flat-mf pick','iwarp warp=${SOURCES[1]} inv=n | window n1=800')

Flow('dif0','gath flat-emd-rec0','add scale=1,-1 ${SOURCES[1]}')
Flow('dif0-mf','gath flat-mf-rec0','add scale=1,-1 ${SOURCES[1]}')

Flow('dif-t1','dif0','window n1=287')
Flow('dif-t2','dif0','window f1=322')
Flow('dif-t3-0','dif0','window f1=287 n1=35')
Flow('dif-t3','dif-t3-0','threshold1 ifperc=1 thr=20 | add scale=-1,1 ${SOURCES[0]}')
Flow('dif','dif-t1 dif-t3 dif-t2','cat axis=1 ${SOURCES[1:3]}')
Flow('flat-emd-rec','gath dif','add scale=1,-1 ${SOURCES[1]}')
Flow('flat-mf-rec','gath dif0-mf','add scale=1,-1 ${SOURCES[1]}')

Grey('gath','title="Original" screenratio=1.4')
Grey('flat','title="Flattened" screenratio=1.4')
Grey('flat-emd','title="EMD" screenratio=1.4')
Grey('flat-emd-rec','title="EMD Filtered" screenratio=1.4')
Grey('flat-mf','title="MF" screenratio=1.4')
Grey('flat-mf-rec','title="MF Filtered" screenratio=1.4')
Grey('dif','title="Diff (EMD)" screenratio=1.4')
Grey('dif0-mf','title="Diff (MF)" screenratio=1.4')

Flow('zooma-1','gath','window f1=250 n1=200 f2=5 n2=25')
Flow('zooma-2','flat-emd-rec','window f1=250 n1=200 f2=5 n2=25')
Flow('zooma-3','flat-mf-rec','window f1=250 n1=200 f2=5 n2=25')

Flow('zoomb-1','gath','window f1=550 n1=200 f2=30 n2=25')
Flow('zoomb-2','flat-emd-rec','window f1=550 n1=200 f2=30 n2=25')
Flow('zoomb-3','flat-mf-rec','window f1=550 n1=200 f2=30 n2=25')

Grey('zooma-1','title="Zoomed A (original)" screenratio=1.4')
Grey('zooma-2','title="Zoomed A (Filtered)" screenratio=1.4')
Grey('zooma-3','title="Zoomed A (MF)" screenratio=1.4')
Grey('zoomb-1','title="Zoomed B (original)" screenratio=1.4')
Grey('zoomb-2','title="Zoomed B (Filtered)" screenratio=1.4')
Grey('zoomb-3','title="Zoomed B (MF)" screenratio=1.4')

## Creating framebox1
x=50
y=0.64
w=80
w1=0.56

Flow('frame1.asc',None,'echo %s n1=10 data_format=ascii_float in=$TARGET'% \
        ' '.join(map(str,(x,y,x+w,y,x+w,y+w1,x,y+w1,x,y))))
Plot('frame1','frame1.asc',
        '''
        dd type=complex form=native |
        graph min1=0 max1=471 min2=0 max2=2.048 pad=n plotfat=15 plotcol=4 
        wantaxis=n wanttitle=n yreverse=y screenratio=1.4 
        ''')

## Creating framebox2
x=140
y=1.40
w=80
w1=0.56

Flow('frame2.asc',None,'echo %s n1=10 data_format=ascii_float in=$TARGET'% \
        ' '.join(map(str,(x,y,x+w,y,x+w,y+w1,x,y+w1,x,y))))
Plot('frame2','frame2.asc',
        '''
        dd type=complex form=native |
        graph min1=0 max1=471 min2=0 max2=2.048 pad=n plotfat=15 plotcol=2 
        wantaxis=n wanttitle=n yreverse=y screenratio=1.4 
        ''')

## Create label A
Plot('labela',None,
        '''
        box x0=2.8 y0=5.8 label="A" xt=0.5 yt=-0.5 length=0.75 
        ''')

## Create label B
Plot('labelb',None,
        '''
        box x0=3.5 y0=3.8 label="B" xt=-0.5 yt=0.5 length=0.75
        ''')

Result('gath0','Fig/gath.vpl frame1 frame2 labela labelb','Overlay')
Result('flat-emd-rec0','Fig/flat-emd-rec.vpl frame1 frame2 labela labelb','Overlay')
Result('flat-mf-rec0','Fig/flat-mf-rec.vpl frame1 frame2 labela labelb','Overlay')


End()

sfdd
sfput
sfgrey
sfmath
sfpad
sfdip
sfwindow
sfpwpaint
sfcontour
sfiwarp
sfenvelope
sfmax1
sfreal
sfspike
sfsmooth
sfmutter
sfadd
sftransp
sfmean
sfthreshold1
sfcat
sfgraph
sfbox