Nonstationarity: patching

PATCHING TECHNOLOGY

A plane of information, either data or an image, say wall(nwall1, nwall2), will be divided up into an array of overlapping windows each window of size (nwind1,nwind2). To choose the number of windows, you specify (npatch1,npatch2). Overlap on the 2-axis is measured by the fraction (nwind2*npatch2)/nwall2. We turn to the language of F90 which allows us to discuss $N$ -dimensional hypercubes almost as easily as two-dimensional spaces. We define an $N$ -dimensional volume (like the wall) with the vector nwall= (nwall1, nwall2, ...). We define subvolume size (like a 2-D window) with the vector nwind=(nwind1, nwind2, ...). The number of subvolumes on each axis is npatch=(npatch1, npatch2, ...). The operator patch simply grabs one patch from the wall, or when used in adjoint form, it puts the patch back on the wall. The number of patches on the wall is product(npatch). Getting and putting all the patches is shown later in module patching .

The $i$ -th patch is denoted by the scalar counter ipatch. Typical patch extraction begins by taking ipatch, a C linear index, and converting it to a multidimensional subscript jj each component of which is less than npatch. The patches cover all edges and corners of the given data plane (actually the hypervolume) even where nwall/npatch is not an integer, even for axes whose length is not an integer number of the patch length. Where there are noninteger ratios, the spacing of patches is slightly uneven, but we'll see later that it is easy to reassemble seamlessly the full plane from the patches, so the unevenness does not matter. You might wish to review the utilities line2cart and cart2line which convert between multidimensional array subscripts and the linear memory subscript before looking at the patch extraction-putback code:

user/gee/patch.c

void patch_lop (bool adj, bool add, int nx, int ny, 
		float* wall, float* wind)
/*< patch operator >*/
{
    int i, j, shift;
 
    sf_adjnull (adj, add, nx, ny, wall, wind);
    sf_line2cart(dim, npatch, ipatch, jj);  
    for(i = 0; i < dim; i++) {
	if(npatch[i] == 1) {
	    jj[i] = 0;
	} else if (jj[i] == npatch[i]-1) {
	    jj[i] = nwall[i] - nwind[i];
	} else {	    
	    jj[i] = jj[i]*(nwall[i] - nwind[i])/(npatch[i] - 1.0);
	}
    }

    /* shift till the patch start */
    shift = sf_cart2line(dim, nwall, jj); 
    for(i = 0; i < ny; i++) {
	sf_line2cart(dim, nwind, i, ii); 
	j = sf_cart2line(dim, nwall, ii) + shift;   
	if (adj) wall[j] += wind[i];
	else     wind[i] += wall[j];
    }
}

The cartesian vector jj points to the beginning of a patch, where on the wall the (1,1,..) coordinate of the patch lies. Obviously this begins at the beginning edge of the wall. Then we pick jj so that the last patch on any axis has its last point exactly abutting the end of the axis. The formula for doing this would divide by zero for a wall with only one patch on it. This case arises legitimately where an axis has length one. Thus we handle the case npatch=1 by abutting the patch to the beginning of the wall and forgetting about its end. As in any code mixing integers with floats, to guard against having a floating-point number, say 99.9999, rounding down to 99 instead of up to 100, the rule is to always add .5 to a floating point number the moment before converting it to an integer. Now we are ready to sweep a window to or from the wall. The number of points in a window is size(wind) or equivalently product(nwind).

Figure 2 shows an example with five nonoverlapping patches on the 1-axis and many overlapping patches on the 2-axis.

parcel Figure 2. A plane of identical values after patches have been cut and then added back. Results are shown for nwall=(100,30), nwind=(17,6), npatch=(5,11). For these parameters, there is gapping on the horizontal axis and overlap on the depth axis.

Nonstationarity: patching