next up previous [pdf]

Next: 2-D filtering in patches Up: PATCHING TECHNOLOGY Previous: PATCHING TECHNOLOGY

Weighting and reconstructing

The adjoint of extracting all the patches is adding them back. Because of the overlaps, the adjoint is not the inverse. In many applications, inverse patching 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 the patching reconstruction operator $ \tilde{\bold I}_p$ in

$\displaystyle \tilde {\bold d} \quad = \quad [ \bold W_{\rm wall} \bold P\T \bold W_{\rm wind} \bold P ] \bold d \quad = \quad \tilde{\bold I}_p \, \bold d$ (1)

where $ \bold d$ is your initial data, $ \tilde{\bold d}$ is the reconstructed data, $ \bold P$ is the patching operator, $ \bold P\T$ is adjoint patching (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 module mkwallwt below builds the weighting function $ \bold W_{\rm wall}$ that you need to apply to your wall of reconstructed data, so it will 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 use adjoint patch [*] to add the weights of each window onto the wall and finally to invert the sum wherever it is non-zero. (You lose data wherever the sum is zero).

user/gee/mkwallwt.c
void mkwallwt(int dim     /* number of dimensions */, 
	      int* npatch /* number of patches [dim] */, 
	      int* nwall  /* data size [dim] */, 
	      int* nwind  /* patch size [dim] */, 
	      float* windwt /* window weighting (input) */, 
	      float* wallwt /* wall weighting (output) */)
/*< make wall weight >*/
{
    int i, j, ip, np, n, nw;

    np = 1; 
    n = 1;
    nw = 1;

    for (j=0; j < dim; j++) {
	np *= npatch[j];
	n *= nwall[j];
	nw *= nwind[j];
    }

    for (i = 0; i < n; i++) {
	wallwt[i] = 0.;
    }

    patch_init(dim, npatch, nwall, nwind);
  
    for (ip=0; ip < np; ip++) {
	patch_lop(true, true, n, nw, wallwt, windwt);
	patch_close ();
    }

    for (i = 0; i < n; i++) {
	if ( wallwt[i] != 0.) wallwt[i] = 1. / wallwt[i];
    }
}

No matrices are needed to show that 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$ . Thus, our process is simply a ``partition of unity.''

To demonstrate the program, I made a random weighting function to use in each window with positive random numbers. The general strategy allows us to use different weights in different windows. That flexibility adds clutter, however, so here we simply use the same weighting function in each window.

The operator $ \tilde{\bold I}_p$ is called ``idempotent.'' The word ``idempotent'' means ``self-power,'' because for any $ N$ , $ 0^N=0$ and $ 1^N=1$ , thus the numbers 0 and 1 share the property that raised to any power they remain themselves. Likewise, the patching reconstruction operator multiplies every data value by either one or zero. Figure 3 shows the result obtained when a plane of identical constant values $ \bold d$ is passed into the patching reconstruction operator $ \tilde{\bold I}_p$ . The result is constant on the 2-axis, which confirms that there is adequate sampling on the 2-axis, and although the weighting function is made of random numbers, all trace of random numbers has disappeared from the output. On the 1-axis the output is constant, except for being zero in gaps, because the windows do not overlap on the 1-axis.

idempatch
Figure 3.
A plane of identical values passed through the idempotent patching reconstruction operator. Results are shown for the same parameters as Figure 2.
idempatch
[pdf] [png] [scons]

Module patching assists in reusing the patching technique. It takes a linear operator $ \bold F$ . as its argument and applies it in patches. Mathematically, this is $ [\bold W_{\rm wall} \bold P\T \bold W_{\rm wind} \bold F \bold P ] \bold d$ . It is assumed that the input and output sizes for the operator oper are equal.

user/gee/patching.c
    for (i=0; i < n; i++) data[i] = 0.;

    patch_init(dim, npatch, nwall, nwind);
    for (ip = 0; ip < np; ip++) {
	/* modl -> winmodl */
	patch_lop(false, false, n, nw, modl, winmodl);
	/* winmodl -> windata */
	oper(false, false, nw, nw, winmodl, windata);
	/* apply window weighting */
	for (iw=0; iw < nw; iw++) windata[iw] *= windwt[iw];
	/* data <- windata */
	patch_lop(true, true, n, nw, data, windata);
	patch_close();
    }

    /* windwt -> wallwt */
    mkwallwt(dim, npatch, nwall, nwind, windwt, wallwt);

    /* apply wall weighting */
    for (i=0; i < n; i++) data[i] *= wallwt[i];


next up previous [pdf]

Next: 2-D filtering in patches Up: PATCHING TECHNOLOGY Previous: PATCHING TECHNOLOGY

2013-07-26