Homework 2

Derivative filters

In this part of the assignment, we return to the digital elevation map of the of the of the Mount St. Helens area, shown in Figure 4.

data Figure 4. Digital elevation map of Mount St. Helens area.

Figure 5 shows a directional derivative, a digital approximation to

$$\cos{\alpha}\,\frac{\partial}{\partial x_1} + \sin{\alpha}\,\frac{\partial}{\partial x_2}\;,$$ (3)

applied to the data. A directional derivative highlights the structure of the mountain as if illuminating it with a light source.

der Figure 5. Directional derivative of elevation.

Figure 6 shows an application of helical derivative, a filter designed by spectral factorization of the Laplacian filter

$$L(Z_1,Z_2) = 4 - Z_1 - 1/Z_1 - Z_2 - 1/Z_2\;.$$ (4)

To invert the Laplacian filter, we put on a helix, where it takes the form

$$L_H(Z) = 4 - Z - Z^{-1} - Z^{N_1} - Z^{-N_1}\;,$$ (5)

and factor it into two minimum-phase parts $L_H(Z) = D(Z)\,D(1/Z)$ using the Wilson-Burg algorithm. The helical derivative $D(Z)$ enhances the image but does not have a preferential direction.

helder Figure 6. Helix derivative of elevation.

1. Change directory to hw2/helix.
2. Run

   scons view
   

   to reproduce the figures on your screen.
3. Edit the SConstruct file. Find the parameter that corresponds to $\alpha$ in equation (3) and try to modify it until you create the most interesting image. After changing the parameter, you can view the result by running

   scons der.view
   

4. EXTRA CREDIT for suggesting and implementing a method for finding optimal $\alpha$ automatically.
5. A more accurate version of the Laplacian filter is
   

   $$\hat{L}_2(Z_1,Z_2) = 20 - 4\,Z_1 - 4\,Z_1^{-1} - 4\,Z_2 - 4\,Z_2^{-1}$$

   $$- Z_1\,Z_2 - Z_1\,Z_2^{-1} - Z_2\,Z_1^{-1} - Z_1^{-1}\,Z_2^{-1}\;.$$ (6)

   
   
   Modify the SConstruct file to use filter (6) instead of (4).
6. EXTRA CREDIT An even more accurate version of the Laplacian filter involves polynomial division in addition to polynomial multiplication. Find the coefficient for the polynomial division and modify the SConstruct file to implement the corresponding helical derivative.

from rsf.proj import *
import math

# Download data
txt = 'st-helens_after.txt' 
Fetch(txt,'data',
      server='https://raw.githubusercontent.com',
      top='agile-geoscience/notebooks/master')
Flow('data.asc',txt,'/usr/bin/tail -n +6')

# Convert to RSF format
Flow('data','data.asc',
     '''
     echo in=$SOURCE data_format=ascii_float
     label=Elevation unit=m
     n1=979  o1=557.805  d1=0.010030675 label1=X
     n2=1400 o2=5107.965 d2=0.010035740 label2=Y |
     dd form=native | 
     clip2 lower=0 | lapfill grad=y niter=200 
     ''')

Result('data','grey title="Digital Elevation Map" allpos=y')

# Vertical and horizontal derivatives
Flow('der1','data','igrad')
Flow('der2','data','transp | igrad | transp')

# !!! CHANGE BELOW !!!
alpha=30

# Directional derivative
Flow('der','der1 der2',
     '''
     add ${SOURCES[1]} scale=%g,%g
     ''' % (math.cos(alpha*math.pi/180),
            math.sin(alpha*math.pi/180)))

Result('der',
       '''
       grey title="Directional Derivative at %gô_" 
       ''' % alpha)

# Laplacian filter on a helix

Flow('slag0.asc',None,
     '''echo 1 1000 n=1000,1000 n1=2 in=$TARGET 
     data_format=ascii_int
     ''')
Flow('slag','slag0.asc','dd form=native')

Flow('ss0.asc','slag',
     '''echo -1 -1 a0=2 n1=2
     lag=$SOURCE in=$TARGET data_format=ascii_float''')
Flow('ss','ss0.asc','dd form=native')

# Wilson-Burg factorization

na=50 # filter length

lags = range(1,na+1) + range(1001-na,1001)

Flow('alag0.asc',None,
     '''echo %s n=1000,1000 n1=%d in=$TARGET 
     data_format=ascii_int
     ''' % (' '.join(map(str,lags)),2*na))
Flow('alag','alag0.asc','dd form=native')

Flow('hflt hlag','ss alag',
     'wilson lagin=${SOURCES[1]} lagout=${TARGETS[1]}')

# Helical derivative

Flow('helder','data hflt hlag','helicon filt=${SOURCES[1]}')
Result('helder','grey title="Helical Derivative" ')

End()

Homework 2