Image-guided Interpolation

aigi/imports.py

@@ -0,0 +1,22 @@
+"""
+Common imports for analysis of Teapot Dome data.
+"""
+import sys
+
+from java.awt import *
+from java.io import *
+from java.lang import *
+from java.util import *
+from javax.swing import *
+
+from edu.mines.jtk.awt import *
+from edu.mines.jtk.dsp import *
+from edu.mines.jtk.interp import *
+from edu.mines.jtk.io import *
+from edu.mines.jtk.mosaic import *
+from edu.mines.jtk.ogl.Gl import *
+from edu.mines.jtk.sgl import *
+from edu.mines.jtk.util import *
+from edu.mines.jtk.util.ArrayMath import *
+
+from aigi import *

aigi/notes.txt

@@ -0,0 +1,3 @@
+Image-guided Interpolation (R) and its adjoint (R^T) flow:
+  compute eigentensors: tensors.py -> tpet.dat
+  apply RR^T: rrt.py

aigi/rrt.py

@@ -0,0 +1,96 @@
+"""
+Interpolates scattered data, such as data from well logs.
+"""
+from utils import *
+
+setupForSubset()
+s1,s2,s3 = getSamplings()
+n1,n2,n3 = s1.count,s2.count,s3.count
+
+tpDir = "/Users/yma/work/data/tp/"
+sfile = tpDir+"tpsts" # seismic image (e.g., gradient) to be operated by RR^T
+efile = "./tpet" # eigen-tensors (structure tensors)
+esfile = "./tpets" # eigen-tensors scaled by semblances
+gfile = "./tpg" # known sample locations, null for unknown samples
+tfile = "./times" # times to nearest known samples
+mfile = "./marks" # marks to nearest known samples
+agfile = "./tpag" # result of applying R^T
+pfile = "./tpp" # nearest-neighbor interpolation
+qfile = "./tpq" # image-guided blendend interpolation (result of applying R)
+
+rrtDir = "./"
+
+pnull = -99.9
+
+def main(args):
+#  goTimeMarkers() # computing time and markers 
+#  goAdjointIGI()  # apply R^T
+#  goIGI()         # apply R
+  display(sfile)
+  display(pfile)
+  display(qfile)
+
+def goTimeMarkers():
+  e = getEigenTensors()
+  g = goSamples()
+  writeImage(gfile,g)
+  igi = IGI3(e,rrtDir)
+  igi.saveTimeMarker(pnull,g)
+
+def goSamples():
+  # Generate known sample locations, marked with 0.0 values, and
+  # others(unknown) are marked with pnull.
+  # In this test, I simply use a set of uniform sample locations, and
+  # you can use more advanced methods to generate sample locations. 
+  s = readImage(sfile,n1,n2,n3)
+  g = zerofloat(n1,n2,n3)
+  g = add(g,pnull)
+  k1,k2,k3 = 20,9,2
+  j1,j2,j3 = 24,20,20
+  f1,f2,f3 = 10,10,5
+  for i3 in range(k3):
+    for i2 in range(k2):
+      for i1 in range(k1):
+        g[f3+i3*j3][f2+i2*j2][f1+i1*j1] = 0.0
+  return g  
+
+def goAdjointIGI():
+  e = getEigenTensors()
+  aigi = AdjointIGI3(e,rrtDir)
+  g = readImage(gfile,n1,n2,n3)
+  x = readImage(sfile,n1,n2,n3)
+  ax = zerofloat(n1,n2,n3)
+  aigi.gridBlended(x,ax)
+  ag = aigi.adjointNearest(pnull,g,ax)
+  writeImage(agfile,ag)
+
+def goIGI():
+  e = getEigenTensors()
+  igi = IGI3(e,rrtDir)
+  g = readImage(gfile,n1,n2,n3)
+  ag = readImage(agfile,n1,n2,n3)
+  p = igi.gridNearest(pnull,g,ag)
+  writeImage(pfile,p)
+  q = zerofloat(n1,n2,n3)  
+  igi.gridBlended(p,q)
+  writeImage(qfile,q)  
+
+def getEigenTensors():
+  e = readTensors(esfile)
+  return e
+
+def display(sfile):
+  s = readImage(sfile,n1,n2,n3)
+  world = World()
+  ipg = addImageToWorld(world,s)
+  makeFrame(world)
+
+def display2(sfile,tmfile):
+  s = readImage(sfile,n1,n2,n3)
+  tm = readImage(tmfile,n1,n2,n3)
+  world = World()
+  ipg = addImage2ToWorld(world,s,tm)
+  makeFrame(world)
+
+#############################################################################
+run(main)

aigi/tensors.py

@@ -0,0 +1,100 @@
+"""
+Computes structure eigen-tensors.
+"""
+from utils import *
+
+setupForSubset()
+s1,s2,s3 = getSamplings()
+n1,n2,n3 = s1.count,s2.count,s3.count
+
+tpDir = "/scratch/data/seis/tp/"
+sfile = tpDir+"tpsts" # seismic image, for computing structure tensors
+efile = "./tpet" # eigen-tensors
+esfile = "./tpets" # eigen-tensors scaled by semblance
+s1file = "./tps1" # semblance w,uv
+s2file = "./tps2" # semblance vw,u
+s3file = "./tps3" # semblance uvw,
+
+lsf1 = LocalSemblanceFilter(2,2)
+lsf2 = LocalSemblanceFilter(2,8)
+lsf3 = LocalSemblanceFilter(16,0)
+
+def main(args):
+  makeTensors()
+  semblence()
+  scaleTensors(0.001) # modify eigenvalues with semblences
+#  display()
+
+def makeTensors():
+  sigma = 8.0
+  s = readImage(sfile,n1,n2,n3)
+  lof = LocalOrientFilter(sigma)
+  e = lof.applyForTensors(s)
+  writeTensors(efile,e)
+
+def semblence():
+  semblance1()
+  semblance2()
+  semblance3()
+  #display2(sfile,s1file)
+  #display2(sfile,s2file)
+  #display2(sfile,s3file)
+
+def semblance1():
+  e = readTensors(efile)
+  s = readImage(sfile,n1,n2,n3)
+  s1 = lsf1.semblance(LocalSemblanceFilter.Direction3.W,e,s)
+  writeImage(s1file,s1)
+
+def semblance2():
+  e = readTensors(efile)
+  s = readImage(sfile,n1,n2,n3)
+  s2 = lsf2.semblance(LocalSemblanceFilter.Direction3.VW,e,s)
+  writeImage(s2file,s2)
+
+def semblance3():
+  e = readTensors(efile)
+  s = readImage(sfile,n1,n2,n3)
+  s3 = lsf3.semblance(LocalSemblanceFilter.Direction3.UVW,e,s)
+  writeImage(s3file,s3)
+
+def scaleTensors(eps):
+  e = readTensors(efile)
+  s1 = readImage(s1file,n1,n2,n3); print "s1 min =",min(s1)," max =",max(s1)
+  s2 = readImage(s2file,n1,n2,n3); print "s2 min =",min(s2)," max =",max(s2)
+  s3 = readImage(s3file,n1,n2,n3); print "s3 min =",min(s3)," max =",max(s3)
+  pow(s1,4.0,s1)
+  pow(s2,8.0,s2)
+  pow(s3,4.0,s3)
+  s1 = clip(eps,1.0,s1)
+  s2 = clip(eps,1.0,s2)
+  s3 = clip(eps,1.0,s3)
+  e.setEigenvalues(s3,s2,s1)
+  writeTensors(esfile,e)
+
+def display():
+  s = readImage(sfile,n1,n2,n3)
+  et = readTensors(esfile)
+  world = World()
+  ipg = addImageToWorld(world,s)
+  ipg.setClips(-5.5,5.5)
+  addTensorsInImage(ipg.getImagePanel(Axis.X),et,20)
+  addTensorsInImage(ipg.getImagePanel(Axis.Y),et,20)
+  addTensorsInImage(ipg.getImagePanel(Axis.Z),et,20)
+  frame = makeFrame(world)
+  frame.setSize(1460,980)
+  frame.orbitView.setAzimuth(-65.0)
+  background = Color(254,254,254)
+  frame.viewCanvas.setBackground(background)
+
+def display2(sfile,smfile):
+  s = readImage(sfile,n1,n2,n3)
+  sm = readImage(smfile,n1,n2,n3)
+  print "semblance: min =",min(sm)," max =",max(sm)
+  world = World()
+  ipg = addImage2ToWorld(world,s,sm)
+  ipg.setClips2(0,1)
+  makeFrame(world)
+
+#############################################################################
+run(main)