LOMOPLAN FORM IN 3-D

Three-dimensional LOMOPLAN is somewhat like multiple passes of two-dimensional LOMOPLAN, i.e. we first LOMOPLAN the (t,x)-plane for each y, and then we LOMOPLAN the (t,y)-plane for each x. Actually, 3-D LOMOPLAN is a little more complicated than this. Each LOMOPLAN filter is designed on all the data in a small (t,x,y) volume.

I began from my earlier two-dimensional code shown in Claerbout (1992c) and made the obvious extensions to three dimensions. For example, the obvious extension to the 2-D code for convolution with a 3-D filter not being allowed to flow outside the faces of the input volume is in subroutine cinlof3() .

# CINLOF -- Convolution INternal with Lags.  Output is conjugate with FILTER.
#
subroutine cinlof3( conj,add, lg1,lg2,lg3,  xx,n1,n2,n3,  bb,nb1,nb2,nb3,  yy)
integer             conj,add, lg1,lg2,lg3,     n1,n2,n3,     nb1,nb2,nb3
real                                        xx(n1,n2,n3), bb(nb1,nb2,nb3)
real                                        yy(n1,n2,n3)
integer y1,y2,y3,  x1,x2,x3,   b1, b2, b3
call conjnull( conj, add,  bb,nb1*nb2*nb3,  yy,n1*n2*n3)

do b3=1,nb3 {  do y3= 1+nb3-lg3, n3-lg3+1 {  x3= y3 - b3 + lg3
do b2=1,nb2 {  do y2= 1+nb2-lg2, n2-lg2+1 {  x2= y2 - b2 + lg2
do b1=1,nb1 {  do y1= 1+nb1-lg1, n1-lg1+1 {  x1= y1 - b1 + lg1
	if( conj == 0 )
		yy( y1,y2,y3) = yy( y1,y2,y3) + bb( b1,b2,b3) * xx( x1,x2,x3)
	else
		bb( b1,b2,b3) = bb( b1,b2,b3) + yy( y1,y2,y3) * xx( x1,x2,x3)
	}}  }}  }}
return; end

Likewise, instead of cutting a window of data into patches, we cut a volume of data into subcubes with patch3() .

# patch3 -- copy the j[123]-th of k[123] subcubes from a volume.
#
subroutine patch3( conj,add, j1,j2,j3, k1,k2,k3, wall, n1,n2,n3, wind, w1,w2,w3)
integer            conj,add, j1,j2,j3, k1,k2,k3,       n1,n2,n3,       w1,w2,w3
integer i1,i2,i3, s1,s2,s3, d1,d2,d3
real                                             wall( n1,n2,n3),wind( w1,w2,n3)
call conjnull(     conj,add,                     wall, n1*n2*n3, wind, w1*w2*w3)

if( k3 != 1) {  s3 = 1.5 + (n3 - w3) * (j3-1.)/(k3-1.)}    else { s3= 1}
if( k2 != 1) {  s2 = 1.5 + (n2 - w2) * (j2-1.)/(k2-1.)}    else { s2= 1}
if( k1 != 1) {  s1 = 1.5 + (n1 - w1) * (j1-1.)/(k1-1.)}    else { s1= 1}

do i3= 1, w3 {	d3= i3 + s3 - 1
do i2= 1, w2 {	d2= i2 + s2 - 1
do i1= 1, w1 {	d1= i1 + s1 - 1
	if( conj == 0 )
		wind( i1,i2,i3) = wind( i1,i2,i3) + wall( d1,d2,d3)
	else
		wall( d1,d2,d3) = wall( d1,d2,d3) + wind( i1,i2,i3)
	}}}
return;	end

Likewise weighting function code and conjugate-gradient code extend in the obvious way. When you finish the obvious extensions, you are left with a little interesting decision making on the shape of the three-dimensional least-squares filters that you want to design. (If you don't like LOMOPLAN, you have prepared yourself an alternative to do any other locally variable three-dimensional least-squares filtering. That is how I came to do steep-dip deconvolution and the time-domain (f,x) filtering described elsewhere in this report.)

For LOMOPLAN, the final step is to specialize to a two-dimensional filter ((t,x) or/and (t,y)-space) in a three-dimensional volume. Luckily, the codes look quite simple and elegant (thanks to the ratfor and sat preprocessors to fortran) and the debugging went quickly.