reproducible version

book/tccs/interpfault/SConstruct

@@ -0,0 +1,3 @@
+from rsf.tex import *
+
+End(options='reproduce',use='hyperref,listings', color='ALL')

book/tccs/interpfault/interpfault-sigmoid/SConstruct

@@ -0,0 +1,155 @@
+from rsf.proj import *
+import math
+
+###### Draw the Box ######
+min1,max1=0.3,0.7
+min2,max2=0.85,1.1
+Flow('box.asc',None,'''echo %s n1=2 n2=5 data_format=ascii_float
+    in=$TARGET'''%' '.join(map(str,(min1,min2,max1,min2,
+                                    max1,max2,min1,max2,
+                                    min1,min2))))
+Plot('box','box.asc','''dd form=native type=complex |
+    window | graph transp=y yreverse=y min1=0 max1=1.2 min2=0 max2=1.28
+    wanttitle=n plotfat=5 plotcol=1 wantaxis=n''')
+
+#####################################
+###### Extract Fault Attribute ######
+#####################################
+Flow('sigmoid',None,'''sigmoid d1=.0032 n1=400 d2=.0032 n2=400 taper=n |
+    smooth rect1=5 diff1=1 | smooth rect1=5''')
+Plot('sigmoid','''window max1=1.2 |
+    grey title="2D Synthetic" label1=Time label2=Distance''')
+Result('sigmoid','sigmoid box','Overlay')
+
+###### Dip Estimation ######
+Flow('dip','sigmoid','dip2 rect1=10 rect2=10 order=4')
+Result('syndip','dip','''window max1=1.2 |
+    grey color=j label2="Distance" title=Dip scalebar=y''')
+
+###### Plane Wave Destruction ######
+Flow('pwd','sigmoid dip','pwd2 dip=${SOURCES[1]} order=3')
+Flow('pws','sigmoid dip','pwsmooth dip=${SOURCES[1]} ns=3 order=3')
+
+###### Sobel Filter ######
+Flow('sobel1','pwd pws','''math pwd=${SOURCES[0]} pws=${SOURCES[1]}
+    output="pwd*pwd+pws*pws"''')
+Result('sobel1','grey title="Old Sobel" allpos=y')
+
+Flow('dip12','dip','transp plane=12')
+Flow('sobel2','pwd dip12',
+     '''transp plane=12 |
+     pwsmooth dip=${SOURCES[1]} ns=3 order=3 |
+     transp plane=12''')
+
+Flow('sobel','sobel1 sobel2','''math s1=${SOURCES[0]} s2=${SOURCES[1]}
+    output="s1*s1+s2*s2"''')
+Result('synsobel','sobel','''window max1=1.2 |
+    grey allpos=y label2="Distance" title="Sobel"''')
+
+Flow('dip2','sobel','odip rect1=10 rect2=10 | transp')
+Result('dip2','transp | grey color=j')
+Flow('smooth','sobel dip2','''transp |
+    pwsmooth dip=${SOURCES[1]} ns=5 order=3 | transp |
+    smooth rect1=5 rect2=5''')
+Result('smooth','''window max1=1.2 |
+    grey allpos=y label2="Distance" title="Smoothed Sobel"''')
+
+
+##########################################
+###### Well Log Predictive Painting ######
+##########################################
+amp=20
+wellLoc = [30,200,370]
+scale=10.0
+r0=100.0
+
+for well in range(3):
+    log = 'log%d'%well
+    Flow(log,'sigmoid','window n2=1 f2=%d'%wellLoc[well])
+    Plot(log,'''scale axis=1 | math output="input*%g+%g" | window max1=1.2 |
+        graph transp=y yreverse=y min2=0 max2=400 min1=0 max1=1.2
+        wherexlabel=top wheretitle=bottom wanttitle=n wantaxis=n
+        plotfat=3'''%(amp,wellLoc[well]))
+
+    paint = 'paint%d'%well
+    Flow(paint,['dip',log],'''pwpaint order=3 seed=${SOURCES[1]} eps=0.1
+        i0=%d'''%wellLoc[well])
+    Plot(paint,'''window max1=1.2 |
+        grey title="Predictive Painting from Log"''')
+    Result(paint,[paint,log,'box'],'Overlay')
+
+    rbfold = 'rbfold%d'%well
+    Flow(rbfold,'smooth','''faultrbf1d useinput=n
+        i0=%d scale=%g r0=%g'''%(wellLoc[well],scale,r0))
+
+    old = 'old%d'%well
+    Flow(old,[paint,rbfold],'add mode=m ${SOURCES[1:2]}')
+
+    rbf = 'rbf%d'%well
+    Flow(rbf,'smooth','''faultrbf1d useinput=y
+        i0=%d scale=%g r0=%g'''%(wellLoc[well],scale,r0))
+    Result(rbf,'''window max1=1.2 | grey color=j scalebar=y bias=0.5
+        title="RBF of log %d"'''%(well+1))
+
+    weighted = 'weighted%d'%well
+    Flow(weighted,[paint,rbf],'add mode=m ${SOURCES[1:2]}')
+
+    seed = 'seed%d'%well
+    Flow(seed,log,'math output="input*0.0"')
+    dist = 'dist%d'%well
+    Flow(dist,['dip',seed,'smooth'],'''distpaint order=3 eps=0.1 i0=%d
+        seed=${SOURCES[1]} flt=${SOURCES[2]} faultscale=5 |
+        clip2 upper=100'''%wellLoc[well])
+    rbfnew = 'rbfnew%d'%well
+    Flow(rbfnew,dist,'''math output="input*300" |
+        math output="1./(1.+(input/%g)*(input/%g))"'''%(r0,r0))
+    Result(rbfnew,'''window max1=1.2 | grey color=j scalebar=y bias=0.5
+        title="New RBF of log %d"'''%(well+1))
+
+    rbfweird = 'rbfweird%d'%well
+    Flow(rbfweird,['dip',seed,'smooth'],'''distpaint order=3 eps=0.1 i0=%d
+        seed=${SOURCES[1]} flt=${SOURCES[2]} faultscale=1 |
+        math output="input*300" |
+        math output="1./(1.+(input/%g)*(input/%g))"
+        '''%(wellLoc[well],r0,r0))
+    Result(rbfweird,'''window max1=1.2 | grey color=j scalebar=y bias=0.5
+        title="RBF of log %d"'''%(well+1))
+
+    weightednew = 'weightednew%d'%well
+    Flow(weightednew,[paint,rbfnew],'add mode=m ${SOURCES[1:2]}')
+
+Flow('rbfoldsum','rbfold0 rbfold1 rbfold2','add ${SOURCES[1:3]}')
+Flow('oldsum','old0 old1 old2','add ${SOURCES[1:3]}')
+Flow('interpold','oldsum rbfoldsum','add mode=d ${SOURCES[1:2]}')
+Plot('interpold','''window max1=1.2 |
+    grey title="RBF interpolation without fault"''')
+Result('interpold','interpold log0 log1 log2 box','Overlay')
+
+Flow('errold','sigmoid interpold',
+    'add scale=1,-1 ${SOURCES[1:2]} | math output="abs(input)"')
+Result('errold','''window max1=1.2 |
+    grey scalebar=y bias=0.0015 clip=0.001 maxval=0.003
+    title="Error of interpolation without fault"''')
+
+Flow('rbfsum','rbf0 rbf1 rbf2','add ${SOURCES[1:3]}')
+Flow('weightedsum','weighted0 weighted1 weighted2','add ${SOURCES[1:3]}')
+Flow('interp','weightedsum rbfsum','add mode=d ${SOURCES[1:2]}')
+Plot('interp','''window max1=1.2 |
+    grey title="RBF interpolation with fault"''')
+Result('interp','interp log0 log1 log2 box','Overlay')
+
+Flow('err','sigmoid interp',
+    'add scale=1,-1 ${SOURCES[1:2]} | math output="abs(input)"')
+Result('err','''window max1=1.2 |
+    grey scalebar=y bias=0.0015 clip=0.001 maxval=0.003
+    title="Error of interpolation with Sobel"''')
+
+Flow('rbfnewsum','rbfnew0 rbfnew1 rbfnew2','add ${SOURCES[1:3]}')
+Flow('weightednewsum','weightednew0 weightednew1 weightednew2',
+    'add ${SOURCES[1:3]}')
+Flow('interpnew','weightednewsum rbfnewsum','add mode=d ${SOURCES[1:2]}')
+Plot('interpnew','''window max1=1.2 |
+    grey title="New RBF interpolation with fault"''')
+Result('interpnew','interpnew log0 log1 log2 box','Overlay')
+
+End()

book/tccs/interpfault/interpfault-teapot/SConstruct

@@ -0,0 +1,129 @@
+from rsf.proj import *
+import math
+
+## ###### Draw the Box ######
+## min1,max1=0.3,0.7
+## min2,max2=0.85,1.1
+## Flow('box.asc',None,'''echo %s n1=2 n2=5 data_format=ascii_float
+##     in=$TARGET'''%' '.join(map(str,(min1,min2,max1,min2,
+##                                     max1,max2,min1,max2,
+##                                     min1,min2))))
+## Plot('box','box.asc','''dd form=native type=complex |
+##     window | graph transp=y yreverse=y min1=0 max1=1.2 min2=0 max2=1.28
+##     wanttitle=n plotfat=5 plotcol=1 wantaxis=n''')
+
+#################################
+###### Import seismic data ######
+#################################
+Flow('teapot','seis','window max1=1.8')
+Result('teapot','''window max1=1.8 |
+    grey title="Teapot Dome"''')
+
+###### Dip Estimation ######
+Flow('dip','teapot','dip2 rect1=10 rect2=10 order=4')
+Result('syndip','dip','''window max1=1.8 |
+    grey color=j label2="Distance" title="Dip Estimation" scalebar=y''')
+
+Flow('fault','flt','window max1=1.8')
+Result('fault','''window max1=1.8 |
+    grey allpos=y pclip=99.9 label2="Distance" title="Fault Attribute"
+    scalebar=y barlabel="Fault Likelihood"''')
+
+
+##########################################
+###### Well Log Predictive Painting ######
+##########################################
+amp=40
+wellLoc = [ 52,203,353, 15, 77,105,133,181,247,246,264,263,285]
+wellMin = [  0, 59,  0,  0,  0,  0,  0, 68, 99,225,103,217,  0]
+wellMax = [313,305,222,111,197,193,120,169,225,290,217,263,113]
+wellLen = [313,246,222,111,197,193,120,101,126, 65,114, 46,113]
+# wellLoc = [ 52,203,353,133,247]
+# wellMin = [  0, 59,  0,  0, 99]
+# wellMax = [313,305,222,120,225]
+# wellLen = [313,246,222,120,126]
+r0=4
+
+logPlotCollect=[]
+log2PlotCollect=[]
+rbfCollect=[]
+weightedCollect=[]
+
+for well in range(13):
+    log = 'log%d'%well
+    Flow(log,'den','window n2=1 f2=%d'%wellLoc[well])
+
+    seed = 'seed%d'%well
+    Flow(seed,log,'mask max=1.0 | dd type=float | math output="input*50"')
+
+    Plot(log,'''window f1=%d n1=%d |
+        scale axis=1 | math output="(input-0.90)*%g+%g" |
+        graph transp=y yreverse=y min2=0 max2=357 min1=0.6 max1=1.8
+        wherexlabel=top wheretitle=bottom wanttitle=n wantaxis=n
+        scalebar=y bartype=v plotcol=7
+        plotfat=3'''%(wellMin[well],wellLen[well],amp,wellLoc[well]))
+    logPlotCollect.append(log)
+
+    log2 = 'logyellow%d'%well
+    Plot(log2,log,'''window f1=%d n1=%d |
+        scale axis=1 | math output="(input-0.90)*%g+%g" |
+        graph transp=y yreverse=y min2=0 max2=357 min1=0.6 max1=1.8
+        wherexlabel=top wheretitle=bottom wanttitle=n wantaxis=n
+        scalebar=n bartype=v plotcol=1
+        plotfat=3'''%(wellMin[well],wellLen[well],amp,wellLoc[well]))
+
+    seis = 'seis%d'%well
+    Result(seis,['teapot',log2],'Overlay')
+
+    log2PlotCollect.append(log2)
+
+    paint = 'paint%d'%well
+    Flow(paint,['dip',log],'''pwpaint order=3 seed=${SOURCES[1]} eps=0.1
+        i0=%d'''%wellLoc[well])
+    Plot(paint,'''window max1=1.8 |
+        grey title="Predictive Painting from Log %d"
+        color=j scalebar=y bartype=v bias=2.3 clip=0.5
+        minval=1.5 maxval=3.0 barlabel="Density (g/cc)"'''%(well+1))
+    Result(paint,[paint,log],'Overlay')
+
+    # thr=1.0
+    # mask = 'mask%d'%well
+    # Flow(mask,paint,'mask min=%g | dd type=float'%thr)
+
+    dist = 'dist%d'%well
+    Flow(dist,['dip',seed,'fault'],'''distpaint order=3 eps=0.1 i0=%d
+        seed=${SOURCES[1]} flt=${SOURCES[2]} faultscale=50 |
+        clip2 upper=500'''%wellLoc[well])
+    Plot(dist,'''window max1=1.8 |
+        grey color=j title="Distance Painting from Log %d"
+        scalebar=y bartype=v barlabel="Geologic Distance"
+        minval=0 maxval=6 bias=2 clip=3'''%(well+1))
+    Result(dist,[dist,log],'Overlay')
+
+    rbf = 'rbf%d'%well
+    Flow(rbf,dist,'math output="1./(1.+(input/%g)*(input/%g))"'%(r0,r0))
+    Plot(rbf,'''window max1=1.8 |
+        grey color=j title="RBF (Center is Log %d)"
+        scalebar=y bartype=v barlabel="RBF Interp Weight"
+        bias=0.8 clip=0.4'''%(well+1))
+
+    Result(rbf,[rbf,log],'Overlay')
+
+    rbfCollect.append(rbf)
+
+    weighted = 'weighted%d'%well
+    Flow(weighted,[paint,rbf],'add mode=m ${SOURCES[1:2]}')
+    weightedCollect.append(weighted)
+
+Flow('rbfsum',rbfCollect,'add ${SOURCES[1:13]}')
+Flow('weightedsum',weightedCollect,'add ${SOURCES[1:13]}')
+Flow('interp','weightedsum rbfsum','add mode=d ${SOURCES[1:2]}')
+Plot('interp','''window max1=1.8 |
+    grey title="Geologic Distance Guided Interpolation" color=j
+    scalebar=y bartype=v minval=1.5 maxval=3.0 barlabel="Density (g/cc)"
+    bias=2.37 clip=0.35 titlesz=9 labelsz=7''')
+Result('interp-tpt','interp '+' '.join(logPlotCollect),'Overlay')
+
+Result('seisall','teapot '+' '.join(log2PlotCollect),'Overlay')
+
+End()

book/tccs/interpfault/interpfault-teapot/den.rsf

@@ -0,0 +1,16 @@
+
+	data_format="native_float"
+	esize=4
+	in="./den.rsf@"
+
+    n1=401
+    d1=0.004
+    o1=0.6
+    label1="Depth"
+    unit1=km
+
+    n2=357
+    d2=0.025
+    o2=0.0
+    label2="Lateral"
+    unit2=km

book/tccs/interpfault/interpfault-teapot/flt.rsf

@@ -0,0 +1,16 @@
+
+	data_format="native_float"
+	esize=4
+	in="./flt.rsf@"
+
+    n1=401
+    d1=0.004
+    o1=0.6
+    label1="Depth"
+    unit1=km
+
+    n2=357
+    d2=0.025
+    o2=0.0
+    label2="Lateral"
+    unit2=km

book/tccs/interpfault/interpfault-teapot/seis.rsf

@@ -0,0 +1,16 @@
+
+	data_format="native_float"
+	esize=4
+	in="./seis.rsf@"
+
+    n1=401
+    d1=0.004
+    o1=0.6
+    label1="Depth"
+    unit1=km
+
+    n2=357
+    d2=0.025
+    o2=0.0
+    label2="Lateral"
+    unit2=km