WEIGHTING AND RECONSTRUCTING

The conjugate of extracting all the patches is adding them back. Because of overlaps, the conjugate is not the inverse. In many applications, inverse parcelling is required, i.e. patching things back together seamlessly. This can be done with weighting functions. You can have any weighting function you wish and I will provide you a reconstruction like

$$\tilde {\bf d} \quad = \quad [ \bold W_{\rm wall} \bold P' \bold W_{\rm wind} \bold P ] \bold d$$ (1)

where $\bold d$ is your initial data, $\tilde{\bf d}$ is the reconstructed data, $\bold P$ is the parcelling operator, $\bold P'$ is conjugate parcelling (adding the patches). $\bold W_{\rm wind}$ is your chosen weighting function in the window, and $\bold W_{\rm wall}$ is the weighting function for the whole wall. You specify any $\bold W_{\rm wind}$ you like, and subroutine mkwallwt() below builds the weighting function $\bold W_{\rm wall}$that you need to apply to your wall of reconstructed data to undo the nasty effects of the overlap of windows and the shape of your window-weighting function. You do not need to change your window weighting function when you increase or decrease the amount of overlap between windows because $\bold W_{\rm wall}$takes care of it. The method is to add the weights of each window onto the wall using conjugate patch() , and finally to invert the sum wherever it is non-zero. (You lose data wherever the sum is zero).

# make wall weights from window weights.
#
subroutine mkwallwt( k1,k2, windwt, w1,w2,  wallwt, n1,n2)
integer i1,i2,       k1,k2,         w1,w2,          n1,n2
real                        windwt( w1,w2), wallwt( n1,n2)
call null(                                  wallwt, n1*n2)
do i1= 1, k1 {
do i2= 1, k2 {
	call patch( 1, 1, i1,i2, k1,k2, wallwt, n1,n2, windwt,w1,w2 )
	}}
do i2= 1, n2 {
do i1= 1, n1 {
	if(  wallwt(i1,i2) != 0.)
	     wallwt(i1,i2)  = 1. / wallwt(i1,i2)
	}}
return; end

No matrices are needed to show this method succeeds because data values are never mixed with one another. An equation for any reconstructed data value $\tilde d$as a function of the original value d and the weights w_i that hit d is $\tilde d = (\sum_i w_i d) / \sum_i w_i = d$.

To demonstrate the program, I begin by making a random weighting function for each window with subroutine mkrandwt(). The only requirement is that the numbers be positive. Zero valued weights are OK too, but you would lose data there unless it was picked up by an overlapping window. The general strategy allows different weights in different windows. That flexibility adds clutter, however, so here we work with the same weighting function in each window.

# Put positive random numbers in a window.
#
subroutine mkrandwt( windwt, w1,w2)
integer i1,i2,	             w1,w2,	iseed
real		     windwt( w1,w2),	rand01
					iseed = 1992
do i1= 1, w1 {
do i2= 1, w2 {
		windwt(i1,i2) =  rand01( iseed)
		}}
return; end

Typically, we apply some linear operator, say $\bold F$ in each patch. Mathematically, this is:

$$\tilde {\bf d} \quad = \quad [ \bold W_{\rm wall} \bold P' ... ... \bold P ]\ \bold d \quad = \quad [ {\tt Idempatch}] \ \bold d$$ (2)

The first code of this nature we will examine (subroutine idempatch()) uses the trivial identity matrix as the linear operator $\bold F$.(Later, after you understand the clutter, we go on to install 2-D filtering as $\bold F$.) To perform the weighting operation we use an adaptation of subroutine diag() from PVI. To show you where to install data processing with any linear operator, I included some trivial processing, multiplying by 1.0. This is done by subroutine ident() , also adapted from PVI. Since the composite operator is also a linear operator I also include code for the conjugate. Including the conjugate nearly doubles the length of the code which means that you hardly need to think about the second half which is a mirror image.

# idempatch -- cut up section and paste back together, idempotently.
#
subroutine idempatch( conj, add, w1,w2, k1,k2, data, n1,n2, fout)
integer i1,i2,	      conj, add, w1,w2, k1,k2,       n1,n2
real					       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 mkrandwt(        windwt, w1,w2)
call mkwallwt( k1,k2, windwt, w1,w2,  wallwt, n1,n2)
if( conj == 0 ) {
	do i2= 1, k2 {
	do i1= 1, k1 {
		call patch( 0, 0,  i1,i2, k1,k2, data,n1,n2, pdat,w1,w2 )
		call ident( 0, 0,   1.0,  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 i2= 1, k2 {
	do i1= 1, k1 {
		call patch( 0, 0,  i1,i2, k1,k2, copy,n1,n2, pfou,w1,w2 )
		call diag(  1, 0, windwt, w1*w2, pfou,       pfou       )
		call ident( 1, 0,   1.0,  w1*w2, pdat,       pfou       )
		call patch( 1, 1,  i1,i2, k1,k2, data,n1,n2, pdat,w1,w2 )
		}}
	}
return; end

Figure 2 shows the result when a plane of identical values is passed through the idempatch() subroutine. The result is constant on the 2-axis which confirms that despite the random function as a weight, idempatch() performs as advertised and that there is adequate sampling on the 2-axis. On the 1-axis the output is constant except for being zero in gaps because the windows do not overlap on the 1-axis. Curiously, the conjugate operation produces the same result, i.e. we can do the operations in reverse order.

 

idempatch Figure 2 A plane of identical values passed through subroutine idempatch(). Results are shown for the same parameters as Figure 1.