2-D FILTERING

A way to do time and space variable filtering is to do invariant filtering within each patch. A convenient subroutine for two-dimensional filtering is cinloi() from PVI page 192 (or SEP-71). For this no-end-effect convolution routine I built a triangular weighting routine cinloiwt() that tapers from the center of the patch of the filter's outputs towards the edges. This is complicated by (1) the constraint that the filter must not move off the edge of the inputs and (2) the alignment of the input and the output. You need not study the details of cinloiwt() because we have seen that the shape of the weighting function is not crucial to fitting things back together. (If you insist on looking at the details, first notice in cinloi() that when the filter index b1 equals lag1 that the output 1-axis pointer y1 matches the input pointer x1 and likewise for the 2-axis.)

# cinloiwt -- get triangle tent weights for 2-D convolution (cinloi()) output  
#
subroutine cinloiwt( a1,a2, lag1,lag2, windwt, w1,w2)
integer 	     a1,a2, lag1,lag2,         w1,w2,     i1,i2, s1,s2, e1,e2
real 				       windwt( w1,w2)
real mid1,mid2, wide1,wide2, x, y
			    call null( windwt, w1*w2)
s1= 1+a1-lag1;	e1= w1-lag1+1;  mid1=(e1+s1)/2.;  wide1=(e1-s1+1.)/2.
s2= 1+a2-lag2;	e2= w2-lag2+1;  mid2=(e2+s2)/2.;  wide2=(e2-s2+1.)/2.

do i1= s1, e1 {	    x = abs((i1-mid1)/wide1)
do i2= s2, e2 {     y = abs((i2-mid2)/wide2)
#	windwt(i1,i2) = amax1( 0., 2. - abs(x+y) - abs(x-y) )		# Cheop
	windwt(i1,i2) = amax1( 0., 1. - abs(x))  * 
			amax1( 0., 1. - abs(y)) 
	}}
return; end

In applications where triangle weights are needed on the inputs (or where we can work on a patch without having interference with edges) we can get ``triangle tent'' weights from cinloiwt() by setting filter dimensions and lags to unity such as shown in Figure 3.

 

windwt Figure 3 Window weights from cinloiwt() with w1=61, w2=19 a1=1, a2=1, lag1=1, lag2=1 .

Triangle weighting functions can sum to a constant if the spacing is such that the midpoint of one triangle is at the beginning of the next. I imagined in 2-D that something similar would happen with shapes like Egyptian pyramids, 2-|x-y|+|x+y|. Instead, this property is possessed by the tent-like shapes in Figure 3 with the equation (1-|x|)(1-|y|) as is shown in Figure 4 where, to add interest, I separated the windows by a little more than the precisely matching distance. In practice we may chose window shapes and overlaps for other reasons than the constancy of the sum of weights because mkwallwt() accounts for that.

 

wallwt Figure 4 (Inverse) wall weights with n1=100, w1=61, k1=2, n2=30, w2=19, k2=2

Inserting the calls to internal convolution cinloi() and cinloiwt() into idempatch() gives the subroutine cinloip() that does the time and space variable filtering by invariant filtering in patches.

The basic rule of computer program management is that we should have no duplicated code. Otherwise improvements might be made to one copy but not the other, and inconsistency could grow. It is an aggravation to me, and a defect in the Fortran language that I must present to you two copies of nearly identical code, idempatch() and cinloip() , and when you or I prepare the next application, we cannot link this code from a library, but we will need to create another copy changing the call to ident() to a call to our next application subroutine. I plan to study the possibility of using interlude subroutines that pick up the required variables through common and the use of other languages.

Finally is an example of filtering a plane of uniform constants with an impulse function. The impulse function is surrounded by zeros so the filter output patches are smaller than the input patches in Figure 2. In Figure 5, both axes need more window density.

 

cinloip Figure 5 Output of cinloi() with the same parameters as Figures 1 and 2. The filter parameters are a1=11 a2=5 lag1=6 lag2=1 . Thus, windows are centered on the 1-axis and pushed back out the 2-axis.

# cinloip -- apply time and space variable 2-D filter, conjugate to input.
#
subroutine cinloip( conj,add, w1,w2,lag1,lag2, aa,a1,a2,k1,k2, data,n1,n2, fout)
integer i1,i2,      conj,add, w1,w2,lag1,lag2,    a1,a2,k1,k2,      n1,n2
real 				     aa( a1,a2, k1,k2), data(n1,n2), fout(n1,n2)
temporary real pdat(w1,w2),   copy(n1,n2)
temporary real pfou(w1,w2), wallwt(n1,n2), windwt(w1,w2)
call conjnull(      conj,add,                           data,n1*n2,  fout,n1*n2)
call cinloiwt( a1,a2, lag1,lag2,             windwt,w1,w2)
call mkwallwt( k1,k2, windwt, w1,w2, wallwt,n1,n2)
if( conj == 0 ) {
	do i1= 1, k1 {
	do i2= 1, k2 {
		call patch(  0, 0, i1,i2, k1,k2, data,n1,n2, pdat,w1,w2 )
		call cinloi( 0, 0, lag1,lag2, a1,a2, aa(1,1,i1,i2),
					  w1,w2, pdat,       pfou       )
		call diag(   0, 0,windwt, w1*w2, pfou,       pfou       )
		call patch(  1, 1, i1,i2, k1,k2, copy,n1,n2, pfou,w1,w2 )
		}}
	call diag(           0, 1,wallwt, n1*n2, copy,       fout       )
	}
else {	call diag(           1, 0,wallwt, n1*n2, copy,       fout       )
	do i1= 1, k1 {
	do i2= 1, k2 {
		call patch(  0, 0, i1,i2, k1,k2, copy,n1,n2, pfou,w1,w2 )
		call diag(   1, 0,windwt, w1*w2, pfou,       pfou       )
		call cinloi( 1, 0, lag1,lag2, a1,a2, aa(1,1,i1,i2),
				          w1,w2, pdat,       pfou       )
		call patch(  1, 1, i1,i2, k1,k2, data,n1,n2, pdat,w1,w2 )
		}}
	}
return; end