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.