up [pdf]
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()

sfdd
sfgrey
sftransp
sfsmooth
sfderiv
sfwindow
sfspray
sfcat
sffft1
sffft3
sfmath
sfcontour
sfadd
sfscale
sfcsv2rsf
sfbyte
sfput
sfinttest1
sfmax
sfmin
sfmask
sfdiv
sfbyte2jpg
sfhistogram
sfotsu
sflabel
sfstack
sfmul
sfmorph
sfpad

https://raw.githubusercontent.com/seg/tutorials-2015/master/1508_Mapping_and_validating_lineaments/continued.txt
https://raw.githubusercontent.com/seg/tutorials-2015/master/1508_Mapping_and_validating_lineaments/xc.txt
https://raw.githubusercontent.com/seg/tutorials-2015/master/1508_Mapping_and_validating_lineaments/yc.txt
https://raw.githubusercontent.com/seg/tutorials-2015/master/1508_Mapping_and_validating_lineaments/nstd.txt
https://raw.githubusercontent.com/seg/tutorials-2015/master/1508_Mapping_and_validating_lineaments/cube1_0-1.csv