# 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