% copyright (c) 1997 Jon Claerbout\title{Nonstationarity: patching}\author{Jon Claerbout}\maketitle\sx{nonstationarity}\sx{patching}\label{paper:pch}\inputdir{XFig}There are many reasonsfor cutting data planes or image planes into overlapping pieces (patches),operating on the pieces, and then putting them back together again,as depicted in Figure~\ref{fig:antoine}.The earth's dip varies with lateral location and depth. The dip spectrum and spatial spectrum thus also varies.The dip itself is the essence of almost all earth mapping,and its spectrum plays an important rolein the estimation any earth properties.In statistical estimation theory,the word to describe changing statistical properties is ``\bx{nonstationary}''.\plot{antoine}{width=5.00in,height=3.5in}{  Decomposing a wall of information into windows (also called patches).  Left is an example of a 2-D space  input to module \texttt{patch}.  Right shows a close-up of the output (top left corner).}%\activesideplot{rayab2D}{width=3.00in,height=1.5in}{NR}{%        Left is space of inputs and outputs.%        Right is during analysis.%        }\parWe begin this chapter with basic patching conceptsalong with supporting utility code.%As in Chapters \ref{iin/paper:iin} and \ref{gem/paper:gem},%I exhibit two-dimensional subroutines only.%Three-dimensional code is an easy extension,%which is in the library (CD-ROM and web)%but not displayed (to reduce clutter).The language of this chapter,{\it patch,}{\it overlap,}{\it window,}{\it wall,}is two-dimensional,but it may as well be three-dimensional,{\it cube,}{\it subcube,}{\it brick,}or one-dimensional,{\it line,}{\it interval.}We sometimes use the language of windows on a wall.But since we usually want to have overlapping windows,better imagery would be to say we assemble a quilt from patches.\parThe codes are designed to work in any number of dimensions.After developing the infrastructure,we examine some two-dimensional, time- and space-variable applications:adaptive steep-dip rejection,noise reduction by prediction,and segregation of signals and noises.\section{PATCHING TECHNOLOGY}A plane of information, either data or an image,\sx{wall}say {\tt wall(nwall1, nwall2)}, will be divided up intoan array of overlapping \bx{window}seach window of size {\tt (nwind1,nwind2)}.To choose the number of windows, you specify {\tt (npatch1,npatch2)}.Overlap on the 2-axis is measured by the fraction{\tt (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\texttt{nwall= (nwall1, nwall2, ...)}.We define subvolume size (like a 2-D window) with the vector\texttt{nwind=(nwind1, nwind2, ...)}.The number of subvolumes on each axis is\texttt{npatch=(npatch1, npatch2, ...)}.The operator\texttt{patch} \vpageref{lst: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\texttt{product(npatch)}.Getting and putting all the patches is shown later in module\texttt{patching} \vpageref{lst:patching}.\parThe $i$-th patch is denoted by the scalar counter \texttt{ipatch}.Typical patch extraction begins by taking\texttt{ipatch}, a C linear index,and converting it to a multidimensional subscript \texttt{jj}each component of which is less than \texttt{npatch}.The patches cover all edges and corners of the given data plane(actually the hypervolume)even where\texttt{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 thatit is easy to reassemble seamlessly the full plane from the patches,so the unevenness does not matter.You might wish to review the utilities\texttt{line2cart} and\texttt{cart2line} \vpageref{lst:cartesian}which convert between multidimensional array subscriptsand the linear memory subscriptbefore looking at the patch extraction-putback code:\opdex{patch}{extract patches}{41}{67}{user/gee}The cartesian vector \texttt{jj}points to the beginning of a patch, where on the wallthe (1,1,..) coordinate of the patch lies.Obviously this begins at the beginning edge of the wall.Then we pick \texttt{jj} so that the last patch on any axishas its last point exactly abutting the end of the axis.The formula for doing this would divide by zerofor a wall with only one patch on it.This case arises legitimately where an axis has length one.Thus we handle the case \texttt{npatch=1} by abutting the patch to thebeginning 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 numberthe 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 \texttt{size(wind)} or equivalently\texttt{product(nwind)}.\inputdir{patch}Figure~\ref{fig:parcel} shows an example with fivenonoverlapping patches on the 1-axis and many overlapping patcheson the 2-axis.\sideplot{parcel}{width=3.00in}{  A plane of identical values  after patches have been cut  and then added back.  Results are shown for  \texttt{nwall=(100,30)},  \texttt{nwind=(17,6)},  \texttt{npatch=(5,11)}.  For these parameters,  there is gapping on the horizontal axis  and overlap on the depth axis.  % n1=100 w1=17 k1=5 n2=30 w2=6 k2=11}\subsection{Weighting and reconstructing}\sx{weighting patches}\parThe adjoint of extracting all the patches is adding them back.Because of the overlaps, the adjoint is not the inverse.In many applications, \bx{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 wishand I will provide youthe patching reconstruction operator $\tilde{\bold I}_p$ in\begin{equation}\tilde {\bf d}\quad = \quad [ \bold W_{\rm wall} \bold P' \bold W_{\rm wind} \bold P ] \bold d\quad = \quad \tilde{\bold I}_p \, \bold d\label{eqn:idemeqn}\end{equation}where $\bold d$ is your initial data,$\tilde{\bf d}$ is the reconstructed data,$\bold P$ is the patching operator,$\bold P'$ 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 functionfor the whole wall.You specify any $\bold W_{\rm wind}$ you like,and module \texttt{mkwallwt} belowbuilds the weighting function $\bold W_{\rm wall}$that you need to applyto your wall of reconstructed data,so it will undo the nasty effects of the overlap of windowsand the shape of your window-weighting function.You do not need to change your window weighting functionwhen you increase or decrease the amount of overlapbetween windows because$\bold W_{\rm wall}$takes care of it.The method is touse adjoint \texttt{patch} \vpageref{lst:patch}to add the weights of each window onto the walland finally to invert the sum wherever it is non-zero.(You lose data wherever the sum is zero).\moddex{mkwallwt}{make wall weight}{25}{59}{user/gee}\parNo 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.''\parTo demonstrate the program,I made a random weighting functionto 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.%\progdex{mkrandwt}{make random wt.}\parThe operator$\tilde{\bold I}_p$is called ``\bx{idempotent}.''The word ``idempotent'' means ``self-power,'' becausefor any $N$,  $0^N=0$ and $1^N=1$,thus the numbers 0 and 1 share the property that raisedto any power they remain themselves.Likewise, the patching reconstruction operatormultiplies every data value by either one or zero.Figure~\ref{fig:idempatch} shows the resultobtained 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 confirmsthatthere 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.\sideplot{idempatch}{width=3.00in}{  A plane of identical values passed through  the idempotent patching reconstruction operator.  Results are shown for the same parameters  as Figure~\protect\ref{fig:parcel}.}\parModule \texttt{patching} assists in reusing the patching technique. Ittakes a linear operator $\bold F$.as its argument and applies it in patches.Mathematically, this is$ [\bold W_{\rm wall} \bold P' \bold W_{\rm wind} \bold F \bold P ] \bold d$.It is assumed that the input and output sizes for the operator\texttt{oper} are equal.\moddex{patching}{generic patching}{51}{70}{user/gee}\par%The first code of this nature we will examine%(subroutine {\tt idempatch()})%uses the trivial identity matrix as the linear%operator $\bold F$.%(Later, after you understand the clutter,%we proceed to install 2-D filtering as $\bold F$.)%To perform the weighting operation we use subroutine%\GPROG{diag}.%To show you where to install data processing with any linear operator,%I included some trivial processing, multiplying by {\tt 1.0}.%This is done by subroutine \GPROG{ident}.%Because the composite operator is also a linear operator%I also include code for the adjoint.%Including the adjoint nearly doubles the length of the code%which means that you hardly need to think about%the second half which is a mirror image.%\progdex{idempatch}{patch inversion}%%\par%Figure~\FIG{idempatch} shows the result%when a plane of identical values is passed through the%{\tt idempatch()} subroutine.\subsection{2-D filtering in patches}A way to do time- and space-variable filtering\sx{filter ! time and space variable}is to do invariant filtering within each patch.Typically, we apply a filter, say $\bold F_p$, in each patch.The composite operator, filtering in patches,$\tilde{\bold F}$,is given by\begin{equation}\tilde {\bf d}\quad = \quad [ \bold W_{\rm wall} \bold P' \bold W_{\rm wind} \bold F_p \bold P ]\  \bold d\quad = \quad \tilde{\bold F}\ \bold d\label{eqn:patchfilt}\end{equation}%A convenient subroutine for two-dimensional filtering is%\texttt{icaf2()} \vpageref{lst:icaf2}.%For this no-end-effect convolution routine,\inputdir{XFig}I built a triangular weighting routine\texttt{tentn()}that tapers from the center of the patch of the filter's {\it outputs}towards the edges.Accomplishing this weighting is complicated by(1) the constraintthat the filter must not move off the edge of the input patchand(2) the alignment of the input and the output.The layout for prediction-error filters is shownin Figure \ref{fig:rabdomain}.\sideplot{rabdomain}{width=1.50in}{  Domain of inputs and outputs of a two-dimensional  prediction-error filter.}\sideplot{rabtent}{width=1.50in}{  Placement of tent-like weighting  function in the space of filter inputs and outputs.}We need a weighting function that vanishes where the filterhas no outputs.The amplitude of the weighting function is not very importantbecause we have learned how to put signals back togetherproperly for arbitrary weighting functions.We can use any pyramidal or tent-like shapethat drops to zero outside the domain of the filter output.The job is done by subroutine {\tt tentn()}.A new parameter needed by \texttt{tentn}is \texttt{a}, the coordinate of the beginning of the tent.%When you are studying diagrams in%Figure \ref{fig:rabdomain}, %to internal filtering%%with subroutine \texttt{icaf2()} \vpageref{lst:icaf2},%it helps to remember that when the filter index {\tt b1}%equals {\tt lag1},%then the output 1-axis pointer {\tt y1}%matches the input pointer {\tt x1} (and likewise for the 2-axis).\moddex{tent}{tent weights}{45}{58}{user/gee}\inputdir{patch}\parIn applications where triangle weights are needed on the {\it inputs}(or where we can work on a patchwithout having interference with edges),we can get ``triangle tent'' weightsfrom {\tt tentn()} if we set filter dimensions and lags to unity,as shown in Figure~\ref{fig:wind0wt}.\sideplot{wind1wt}{width=3.00in}{  Window weights from {\tt tentn()}  with  %	{\tt         w1=61, w2=19 a1=1, a2=1, lag1=1, lag2=1 }.  \texttt{ nwind=(61,19), center=(31,1), a=(1,1) }.}\par    Triangle weighting functions\sx{triangle weighting functions}