from rsf.proj import *
import random, string

random.seed(2005)

Fetch('lena.img','imgs')
Flow('lena','lena.img',
     '''
     echo n1=512 n2=513 in=$SOURCE data_format=native_uchar |
     dd type=float |
     window f2=1
     ''',stdin=0)

def grey(title,allpos=1):
    return '''
    grey transp=n allpos=%d title="%s"
    screenratio=1 wantaxis=n
    ''' % (allpos,title)

Result('lena',grey('Lena'))

Flow('haar1','lena','transp | dwt type=h | transp')
Result('haar1',grey('1-D Haar Transform'))

Flow('haar2','haar1','dwt type=h')
Result('haar2',grey('2-D Haar Transform'))

Flow('linear1','lena','transp | dwt type=l inv=y unit=y | transp')
Result('linear1',grey('1-D Wavelet Transform'))

Flow('linear2','linear1','dwt type=l inv=y unit=y')
Result('linear2',grey('2-D Wavelet Transform'))

for scale in (2,4,8,16,32,64,128,256):
    lena = 'lena%d' % scale
    Flow(lena,'linear2',
         '''
         cut f2=%d f1=%d |
         dwt type=l inv=y adj=y unit=y |
         transp | dwt type=l inv=y adj=y unit=y | transp
         ''' % (scale,scale))
    Result(lena,grey('Scale=%d' % scale))

Flow('slena','lena',
     '''
     transp |
     put d1=0.004 d2=0.01 o2=0 o1=0 |
     fft1 | fft3 |
     dipfilter v1=-1.5 v2=-1 v3=1 v4=1.5 taper=2 pass=0 |
     fft3 inv=y | fft1 inv=y |
     bandpass flo=15 fhi=45 |
     transp
     ''')
Result('slena',grey('Seismic Lena',0))


nsp = 100
k2 = string.join(map(lambda x: str(random.randint(1,512)), xrange(nsp)),',')
k1 = string.join(map(lambda x: str(random.randint(1,512)), xrange(nsp)),',')

Flow('imps',None,
     '''
     spike nsp=%d k1=%s k2=%s n1=512 n2=512 |
     dwt type=l inv=y adj=y |
     transp | dwt type=l inv=y adj=y | transp
     ''' % (nsp,k1,k2))
Result('imps',grey('Wavelets',0))

k1 = string.join(map(lambda x: str(x+1), xrange(512)),',')

Flow('spikes',None,
     'math n1=128 n2=128 output="abs(x1-x2)" | mask max=0.5 | dd type=float')
for type in ('haar','linear'):
    Type = string.capitalize(type)
    imps = 'i'+type
    Flow(imps,'spikes','dwt type=%s inv=y adj=y' % type)
    Result(imps,
           'grey title="%s Wavelets" label1=Time label2=Scale' % Type)
    fft = 'f' + type
    Flow(fft,imps,'spectra all=n')
    Result(fft,'''
    grey title="%s Wavelets Spectra" label1=Frequency label2=Scale
    allpos=y color=j
    ''' % Type)

title = {'lena': 'Lena',
         'slena': 'Seismic Lena',
         'imps': 'Wavelets'}

for case in title.keys():
    Result('f'+case,case,
       '''
       transp |
       put d1=0.004 d2=0.01 o2=0 o1=0 |
       spectra2 |
       sfgrey color=j title="%s Spectrum"
       label2=Wavenumber label1=Frequency allpos=y
       ''' % title[case])

Flow('sdip','slena','transp | dip rect1=20 rect2=20 order=1')
Result('sdip','grey color=j title="Seismic Slope" scalebar=y')

Flow('spwd','slena sdip','transp | pwd order=1 dip=${SOURCES[1]} | transp')
Result('spwd',grey('Seismic Residual',0) + ' clip=34.4275')

Flow('seis','slena sdip',
     'transp | seislet dip=${SOURCES[1]} eps=0.01 adj=y inv=y unit=y | transp')
Result('seis',grey('Seislet Transform',0))

Flow('sinv','seis sdip',
     'transp | seislet dip=${SOURCES[1]} eps=0.01 inv=y unit=y | transp')
Result('sinv',grey('Inverse Seislet Transform',0))

#for scale in (2,4,8,16,32,64,128,256):
#    slena = 'slena%d' % scale
#    Flow(slena,'seis sdip',
#         '''
#         transp | cut f2=%d | seislet dip=${SOURCES[1]} | transp
#         ''' % scale)
#    Result(slena,grey('Scale=%d' % scale,0))

for c in (1,5):
    rec = 'srec%d' % c
    Flow(rec,'seis sdip',
         '''
         threshold pclip=%d |
         transp |
         seislet dip=${SOURCES[1]} inv=y unit=y | transp
         ''' % c)
    Result(rec,grey('Inverse Seislet Transform (%d%%)' % c,0))

    wrec = 'wrec%d' % c
    Flow(wrec,'linear2',
         '''
         threshold pclip=%d |
         dwt adj=y inv=y unit=y | transp |
         dwt adj=y inv=y unit=y | transp
         ''' % c)
    Result(wrec,grey('Inverse Wavelet Transform (%d%%)' % c))

    srec = 'wsrec%d' % c
    Flow(srec,'slena',
         'dwt | threshold pclip=%d | dwt adj=y inv=y' % c)
    Result(srec,grey('Inverse Wavelet Transform (%d%%)' % c,0))

End()