from rsf.proj import *
from math import pi

# import Bouguer gravity residual map
Fetch('continued.txt','1508_Mapping_and_validating_lineaments',
      server='https://raw.githubusercontent.com',
      top='seg/tutorials-2015/master')

Flow('continued','continued.txt',
     '''
     echo n1=81 n2=81 data_format=ascii_float in=$SOURCE 
     o1=8 d1=-0.1 label1=Y o2=0 d2=0.1 label2=X |
     dd form=native
     ''')

Result('continued','grey title="Bouguer anomaly" screenratio=1 scalebar=y transp=n')

for case in ('xc','yc'):
    txt = case+'.txt'
    Fetch(txt,'1508_Mapping_and_validating_lineaments',
          server='https://raw.githubusercontent.com',
          top='seg/tutorials-2015/master')
    Flow(case,case+'.txt',
         '''
         echo n1=81 n2=81 data_format=ascii_float in=$SOURCE |
         dd form=native
         ''')

Flow('dx','continued','transp | smooth rect1=3 | deriv scale=1 | transp')
Flow('dy','continued','smooth rect1=3 | deriv scale=1')

Flow('left','continued','window n2=1 | spray axis=2 n=81')
Flow('right','continued','window n2=1 f2=-1 | spray axis=2 n=81')
Flow('pad2','left continued right','cat axis=2 ${SOURCES[1:3]}')

Flow('top','pad2','window n1=1 | spray axis=1 n=81')
Flow('bottom','pad2','window n1=1 f1=-1 | spray axis=1 n=81')
Flow('pad','top pad2 bottom','cat axis=1 ${SOURCES[1:3]}')

Flow('dz','pad',
     '''
     fft1 | fft3 axis=2 | 
     sfmath output="input*%g*sqrt(x1*x1+x2*x2)" |
     fft3 axis=2 inv=y | fft1 inv=y |
     smooth rect1=3 rect2=3 |
     window n1=81 n2=81 f1=81 f2=81
     ''' % (2*pi))

for case in ('dx','dy','dz'):
    Result(case,'grey title="%s" screenratio=1 scalebar=y transp=n' % case)

Flow('tdx','dx dy','math dx=${SOURCES[0]} dy=${SOURCES[1]} output="sqrt(dx*dx+dy*dy)" ')

def bplot(title):
    return '''
    grey title="%s" 
    screenratio=1 scalebar=y mean=y color=B transp=n 
    ''' % title

Result('tdx',bplot('Total Horizontal Derivative'))

Flow('tilt','dz tdx','math tdx=${SOURCES[1]} output="input&tdx" ')
Flow('dxt','tilt','transp | smooth rect1=3 | deriv scale=1 | transp')
Flow('dyt','tilt','smooth rect1=3 | deriv scale=1')

Flow('thdr','dxt dyt','math dx=${SOURCES[0]} dy=${SOURCES[1]} output="sqrt(dx*dx+dy*dy)" ')

Result('thdr',bplot('THDR'))

Flow('theta','tdx dx dy dz',
     '''
     math dx=${SOURCES[1]} dy=${SOURCES[2]} dz=${SOURCES[3]} 
     output="input/sqrt(dx*dx + dy*dy +dz*dz)"
     ''')

Result('theta',bplot('Theta'))

# import NSTD, normalized standard deviation of data
Fetch('nstd.txt','1508_Mapping_and_validating_lineaments',
      server='https://raw.githubusercontent.com',
      top='seg/tutorials-2015/master')

Flow('nstd','nstd.txt',
     '''
     echo n1=73 n2=73 data_format=ascii_float in=$SOURCE |
     dd form=native
     ''')

Result('nstd',bplot('NSTD'))

Plot('dzcont','dz','contour screenratio=1 transp=n allpos=n nc=1 c0=0 wanttitle=n wantaxis=n plotcol=7 plotfat=3 scalebar=y')
Result('dzcont','Overlay')

Flow('tdxn','tdx dz','math dz=${SOURCES[1]} output="input&abs(dz)" ')

Result('tdxn',bplot('TDXN'))

# rescale to 0-1 range

Flow('tdxn0','tdxn','add add=-0.000685951 | scale axis=2')

color='cube1_0-1.csv'

Fetch(color,'1508_Mapping_and_validating_lineaments',
      server='https://raw.githubusercontent.com',
      top='seg/tutorials-2015/master')

Plot('color',['continued',color],'grey title="Bouguer anomaly" screenratio=1 scalebar=y transp=n color=${SOURCES[1]} pclip=100')
Result('color','Overlay')

# Mac to Unix text conversion
Flow('color',color,"awk '{ gsub(\"\\r\", \"\\n\"); print $0;}' | csv2rsf")

Flow('byte','continued','byte | sfdd type=float | put n2=1 n1=%d' % (81*81))

# RGB to HSV

Flow('rgb','color byte','transp | inttest1 coord=${SOURCES[1]} interp=lag nw=1')

Flow('r','rgb','window n2=1 f2=0')
Flow('g','rgb','window n2=1 f2=1')
Flow('b','rgb','window n2=1 f2=2')

Flow('cmax','rgb','max axis=2')
Flow('cmin','rgb','min axis=2')
Flow('cdif','cmax cmin','add scale=1,-1 ${SOURCES[1]}')

for case in 'rgb':
    Flow(case+'max',[case,'cmax'],
         'add scale=1,-1 ${SOURCES[1]} | mask min=0 max=0 | dd type=float')

Flow('rh','rmax g b cdif',
     'math g=${SOURCES[1]} b=${SOURCES[2]} d=${SOURCES[3]} output="input*(g-b)/d" ')
Flow('gh','gmax b r cdif',
     'math b=${SOURCES[1]} r=${SOURCES[2]} d=${SOURCES[3]} output="input*(2+(b-r)/d)" ')
Flow('bh','bmax r g cdif',
     'math r=${SOURCES[1]} g=${SOURCES[2]} d=${SOURCES[3]} output="input*(4+(r-g)/d)" ')
Flow('h','rh gh bh','add ${SOURCES[1:3]} | math output="input/6" ')

Flow('s','cdif cmax','div ${SOURCES[1]}')

Flow('v','tdxn0','math output=1-input | put n2=1 n1=%d' % (81*81))

# HSV to RGB

Flow('i','h','math output="6*input-0.5" | dd type=int | dd type=float')
Flow('f','i h','math h=${SOURCES[1]} output="6*h-input" ')
Flow('p','v s','math s=${SOURCES[1]} output="input*(1-s)" ')
Flow('q','v s f','math s=${SOURCES[1]} f=${SOURCES[2]} output="input*(1-s*f)" ')
Flow('t','v s f','math s=${SOURCES[1]} f=${SOURCES[2]} output="input*(1-s*(1-f))" ')

masks = []
mline = 'v=${SOURCES[6]} p=${SOURCES[7]} q=${SOURCES[8]} t=${SOURCES[9]}'
for i in range(6):
    mask = 'mask%d' % i
    Flow(mask,'i','math output=input-%g | mask min=0 max=0 | dd type=float' % i)
    masks.append(mask)
    mline += ' m%d=${SOURCES[%d]}' % (i,i)

Flow('r2',masks+Split('v p q t'),'math %s output="m0*v+m1*q+m2*p+m3*p+m4*t+m5*v" ' % mline)
Flow('g2',masks+Split('v p q t'),'math %s output="m0*t+m1*v+m2*v+m3*q+m4*p+m5*p" ' % mline)
Flow('b2',masks+Split('v p q t'),'math %s output="m0*p+m1*p+m2*t+m3*v+m4*v+m5*q" ' % mline)

Flow('plot.jpg','r2 g2 b2',
     'cat axis=2 ${SOURCES[1:3]} | transp | put n2=81 n3=81 | byte allpos=y clip=1 | byte2jpg')

Result('colorcont','color dzcont','Overlay')

# Part 2 - extract lineaments from enhanced maps

otsus = []
for case in ('theta','tdxn0','nstd'):
    otsu = case + '-otsu'
    Flow(otsu+'.par',case,'histogram n1=101 o1=0 d1=0.01 | otsu | sed s/threshold/min/')
    Flow(otsu,[case,otsu+'.par'],'mask par=${SOURCES[1]}')
    Plot(otsu,'dd type=float | grey allpos=y clip=1 wanttitle=n screenratio=1 crowd=1 wantaxis=n transp=n')
    otsus.append(otsu)

Result('otsu',otsus,'SideBySideIso')

threshold = dict(theta=0.9,tdxn0=0.75,nstd=0.45)

bins = []
for case in ('theta','tdxn0','nstd'):
    bin = case + '-bin'
    Flow(bin,case,'mask min=%g' % threshold[case])
    Plot(bin,'dd type=float | grey allpos=y clip=1 wanttitle=n screenratio=1 crowd=1 wantaxis=n transp=n')
    bins.append(bin)

Result('bin',bins,'SideBySideIso')

#identify and remove all white objects smaller than a specific size

lbls = []
bins = []
clss = []
skts = []
for case in ('theta','tdxn0','nstd'):
    # Label
    lbl = case + '-lbl'
    Flow(lbl,case+'-bin','dd type=float | byte pclip=100 allpos=y clip=1 | label')
    Plot(lbl,'dd type=float | grey allpos=y pclip=100 color=j wanttitle=n screenratio=1 crowd=1 wantaxis=n transp=n')
    lbls.append(lbl)

    # Remove small regions 
    msks = []
    for label in range(60):
        msk = lbl + '-msk%d' % label
        Flow(msk+'in',lbl,'mask min=%d max=%d' % (label,label))
        Flow(msk,msk+'in',
             '''
             dd type=float | 
             stack norm=n axis=2 | stack norm=n axis=1 | 
             spray axis=1 n=%d | spray axis=2 n=%d |
             mask min=%d | mul ${SOURCE}
             ''' % ((81,73)[case=='nstd'],(81,73)[case=='nstd'],75))
        msks.append(msk)
    msk = case + '-msk'
    Flow(msk,msks+[case+'-bin'],'add ${SOURCES[1:60]} | mul ${SOURCES[60]}')
    Plot(msk,'dd type=float | grey allpos=y clip=1 wanttitle=n screenratio=1 crowd=1 wantaxis=n transp=n')
    bins.append(msk)

    # Morphological closing
    cls = case + '-cls'
    Flow(cls,msk,'morph what=closing')
    Plot(cls,'dd type=float | grey allpos=y clip=1 wanttitle=n screenratio=1 crowd=1 wantaxis=n transp=n')
    clss.append(cls)

    # Morphological skeleton
    skt = case + '-skt'
    Flow(skt,cls,'morph what=skeleton')
    Plot(skt,'dd type=float | grey allpos=y clip=1 wanttitle=n screenratio=1 crowd=1 wantaxis=n transp=n')
    skts.append(skt)
Flow('nstd-skt1','nstd-skt','pad beg1=4 end1=4 beg2=4 end2=4')

Result('lbl',lbls,'SideBySideIso')
Result('msk',bins,'SideBySideIso')
Result('cls',clss,'SideBySideIso')
Result('skt',skts,'SideBySideIso')

# process using thresholding

thr = 0.27*255
Flow('thre','thdr','byte pclip=100 allpos=y | dd type=float | mask min=%g' % thr)
Result('thre','dd type=float | grey allpos=y clip=1 wanttitle=n screenratio=1 crowd=1 wantaxis=n transp=n')

# Label
Flow('tlbl','thre','dd type=float | byte pclip=100 allpos=y clip=1 | label')
Result('tlbl',
       'dd type=float | grey allpos=y pclip=100 color=j wanttitle=n screenratio=1 crowd=1 wantaxis=n transp=n')

# Remove small regions 
msks = []
for label in range(60):
    msk = 'tlbl-msk%d' % label
    Flow(msk+'in','tlbl','mask min=%d max=%d' % (label,label))
    Flow(msk,msk+'in',
         '''
         dd type=float | 
         stack norm=n axis=2 | stack norm=n axis=1 | 
         spray axis=1 n=%d | spray axis=2 n=%d |
         mask min=%d | mul ${SOURCE}
         ''' % (81,81,75))
    msks.append(msk)
Flow('tmsk',msks+['thre'],'add ${SOURCES[1:60]} | mul ${SOURCES[60]}')
Result('tmsk','dd type=float | grey allpos=y clip=1 wanttitle=n screenratio=1 crowd=1 wantaxis=n transp=n')

# Morphological closing
Flow('tcls','tmsk','morph what=closing')
Result('tcls','dd type=float | grey allpos=y clip=1 wanttitle=n screenratio=1 crowd=1 wantaxis=n transp=n')

# Morphological skeleton
Flow('tskt','tcls','morph what=skeleton')
Result('tskt','dd type=float | grey allpos=y clip=1 wanttitle=n screenratio=1 crowd=1 wantaxis=n transp=n')

# Part 3 - make a lineament collocation map

# enlarge the skeletons a bit using dilation

fats = []
for skt in ('theta-skt','tdxn0-skt','nstd-skt1'):
    fat = skt + '-fat'
    Flow(fat,skt,'morph what=dilation')
    Plot(fat,'dd type=float | grey allpos=y clip=1 wanttitle=n screenratio=1 crowd=1 wantaxis=n transp=n')
    fats.append(fat)
Result('fat',fats,'SideBySideIso')

# collocation

Flow('col',fats,'add ${SOURCES[1:3]}')

# color map

Result('col',
       '''
       dd type=float | 
       grey allpos=y pclip=100 wanttitle=n screenratio=1  
       scalebar=y color=hot wantaxis=n transp=n
       ''')

# collocation map with a mask applied so as to display lineaments only where at least two methods collocate.

Flow('col2','col','mask min=2 | mul $SOURCE')

Result('col2','col',
       '''
       dd type=float | 
       grey allpos=y pclip=100 wanttitle=n screenratio=1  
       scalebar=y color=hot wantaxis=n transp=n bias=1
       ''')

# display the collocated edges superimposed on the TDXN 

Result('edge','col2 tdxn',
        '''
        dd type=float | add ${SOURCES[1]} |
        grey title="TDXN" color=BC clip=1.58 
        screenratio=1 allpos=y transp=n 
        ''')

End()