from rsf.proj import *
from rsf.prog import RSFROOT
import os

Flow('random',None,'math n1=8 n2=8 output=1 | noise rep=y seed=2014 type=n mean=0.5')

Result('spectral','random',
       '''
       grey pclip=100 mean=y color=spectral screenratio=0.9 scalebar=y 
       title=spectral
       ''')

Flow('bar',None,'math n1=16 output=x1')

Result('bar',
       '''
       grey allpos=y pclip=100 color=spectral title=Colors wantaxis=n 
       screenht=2 screenwd=14 transp=n wheretitle=t
       ''')

# Get color pallette

for case in ('spectral','linearlfb'):
       csv = os.path.join(RSFROOT,'include',case+'.csv')
       Flow(case,csv,'csv2rsf')

# Interpolate 16 numbers

Flow('num16','spectral','transp | remap1 o1=0 n1=16 d1=%g | scale dscale=255' % (255.0/15))

# Compute intensity

Flow('intensity','num16','stack scale=0.2989,0.5870,0.1140 norm=n')

Plot('intensity',
     '''
     grey allpos=y pclip=100 title=Intensity wantaxis=n 
     screenht=2 screenwd=14 transp=n wheretitle=t
     ''')

Result('intensity','Overlay')

Plot('intensity2','intensity',
     '''
     graph symbol=o min2=0 max2=255 wanttitle=n label2= unit2= parallel2=n
     screenht=2 screenwd=14 wantaxis1=n labelsz=5 symbolsz=5 plotcol=2
     ''')

Result('intensity2','intensity intensity2','Overlay')

# Sort by intensity

Result('isorted','intensity',
       '''
       transp | headersort head=$SOURCE | window |
       grey allpos=y pclip=100
       title="Sorted color intensity" wantaxis=n 
       screenht=2 screenwd=14 transp=n wheretitle=t
       ''')

Flow('sorted','intensity',
     '''
     math output=x1 |
     transp | headersort head=$SOURCE | window
     ''')

Result('sorted',
       '''
       grey allpos=y pclip=100 color=spectral
       title="Colors sorted by intensity" wantaxis=n 
       screenht=2 screenwd=14 transp=n wheretitle=t
       ''')

# Plotted as Euclidean distance

Flow('euclid','intensity','bin1 x0=0 dx=1 nx=256 head=$SOURCE')

Result('euclid',
       '''
       grey allpos=y pclip=100 
       title="Color intensity, Euclidean distance" wantaxis=n 
       screenht=2 screenwd=14 transp=n wheretitle=t
       ''')

# Great Pyramid of Giza

Result('linear','random',
       '''
       grey pclip=100 mean=y color=linearlfb screenratio=0.9 scalebar=y 
       title=linearlfb
       ''')

Flow('pyramid',None,
     'math n1=101 n2=101 o1=-50 o2=-50 output="100-abs(x1+x2)-abs(x1-x2)"')

Flow('depth',None,'math n1=256 output=x1 | byte allpos=y pclip=100')

for case in ('spectral','linearlfb'):
       Result('pyramid-'+case,'pyramid',
            'grey screenratio=1 scalebar=y allpos=y pclip=100 color=%s wanttitle=n' % case)
       Result('pyramid3d-'+case,'pyramid',
              'plsurf allpos=y pclip=100 color=%s mesh=n' % case)

       # Lightness profile

       Flow(case+'_m',case,'mask min=0.04045 | dd type=float')
       Flow(case+'1',[case,case+'_m'],
            'math mask=${SOURCES[1]} output="mask*((input+0.055)/1.055)^2.4" ')
       Flow(case+'2',[case,case+'_m'],
            'math mask=${SOURCES[1]} output="(1-mask)*input/12.92" ')
       Flow(case+'_y',[case+'1',case+'2'],
            '''
            add ${SOURCES[1]} | stack axis=1 norm=n scale=0.2126,0.7152,0.0722
            ''')
       Flow(case+'y_m',case+'_y','mask min=0.008856 | dd type=float')
       Flow(case+'_y1',[case+'_y',case+'y_m'],
            'math mask=${SOURCES[1]} output="mask*input^(1/3)" ')
       Flow(case+'_y2',[case+'_y',case+'y_m'],
            'math mask=${SOURCES[1]} output="(1-mask)*((7.787*input) + (16/116))" ')

       Flow(case+'-l',[case+'_y1',case+'_y2'],
            'add ${SOURCES[1]} | math output="116*input - 16" ')

       Result(case+'-l',[case+'-l','depth'],
              '''
              graph depth=${SOURCES[1]} color=%s plotfat=25 title=%s
              label2="Lightness L*" label1=Colors min2=0 max2=100 unit1= unit2=
              ''' % (case,case.capitalize()))

# Plot horizon

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 | lapfill grad=y verb=y niter=500 |
     put label=Time unit=ms label1=Inline label2=Crossline
     ''')

for case in ('spectral','linearlfb'):
       Result('horizon-'+case,'horizon',
              'grey color=%s bias=950 clip=133 scalebar=y title=Horizon yreverse=n transp=n' % case)

End()