Euclidean Group Invariant Computation of Stochastic Completion Fields Using Shiftable-Twistable Functions.
We describe a method for computing the likelihood that a completion joining two contour fragments passes through any given position and orientation in the image plane. Like computations in primary visual cortex (and unlike all previous models of contour completion), the output of our computation is invariant under rotations and translations of the input pattern. This is achieved by representing the input, output, and intermediate states of the computation in a basis of shiftable-twistable functions.
Introduction
Any computational model of human visual information processing must reconcile two apparently
contradictory observations. First, computations in primary visual cortex are largely
Euclidean invariant|an arbitrary rotation and translation of the input pattern of light falling
on the retina produces an identical rotation and translation of the output of the computa-
tion. Second, simple calculations based on the size of primary visual cortex (60 mm  80
mm) and the observed density of cortical hypercolumns (4/mm 2 ) suggest that the discrete
spatial sampling of the visual eld is exceedingly sparse [29]. The apparent contradiction becomes
clear when we ask the following questions: How is this remarkable invariance achieved
in computations performed by populations of cortical neurons with broadly tuned receptive
elds centered at so few locations? Why doesn't our perception of the world change
dramatically when we tilt our head by 5 degrees? 1
asks a related question in a recent Nature paper [8]:
\On average, a region of just 1 mm 2 on the surface of the cortex will contain all possible orientation
preferences, and, accordingly, can analyze orientation for one small area of the visual
eld. This topographical arrangement allows closely spaced objects with dierent orientations
to interact. But it also means that a continuous line across the whole visual eld would be
cortically depicted in a patchy, discontinuous fashion. How can the spatially separated elements
be bound together functionally?"
(a) (b)

Figure

1: (a) Ehrenstein Figure. (b) Kanizsa Triangle.
One of the main goals of our research is to show how the sparse and seemingly haphazard
nature of the sampling of the visual eld can be reconciled with the Euclidean invariance of
visual computations. To realize this goal, we introduce the notion of a shiftable-twistable
basis of functions on the space, R 2  S 1 , of positions and directions. This notion is a
generalization of the notion of a shiftable-steerable basis of functions on the plane, R 2 ,
introduced by Freeman, Adelson, Simoncelli, and Heeger in two seminal papers [9, 26].
Freeman and Adelson [9] clearly appreciated the importance of the issues raised above when
they devised the notion of a steerable basis to implement rotationally invariant computations.
In fact, for computations in the plane the contradictions discussed above were largely resolved
with the introduction by Simoncelli et al. [26] of the shiftable-steerable pyramid transform,
which was specically designed to perform Euclidean invariant computations on R 2 . The
basis functions in the shiftable-steerable pyramid are very similar to simple cell receptive
elds in primary visual cortex. However, many computations in V1 and V2 likely operate
on functions of the space of positions and directions, R 2  S 1 , rather than on functions
of the plane, R 2 (e.g., [11, 12, 18, 22, 24, 30, 31, 34]). Consequently, we propose that
shiftability-twistability (in addition to shiftability-steerability) is the property which binds
sparsely distributed receptive elds together functionally to perform Euclidean invariant
computations in visual cortex.
In this article, we describe a new algorithm for completing the boundaries of partially
occluded objects. This algorithm is based on a computational theory of contour completion
in primary and secondary visual cortex developed in recent years by Williams and
colleagues [27, 28, 30, 31]. Like computations in V1 and V2, and unlike previous models
of illusory contour formation, our computation is Euclidean invariant. This invariance is
achieved by representing the input, output, and intermediate states of the computation in a
basis of shiftable-twistable functions.
Mumford [21] proposed that the probability distribution of natural shapes can be modeled
by particles traveling with constant speed in directions given by Brownian motions.
More recently, Williams and Jacobs [30] dened the stochastic completion eld to be the
distribution of particle trajectories joining pairs of position and direction constraints, and
showed how it could be computed in a neural network.
The neural network described in [31] is based on Mumford's observation that the evolution
in time of the probability density function (p.d.f.) representing the position, (x; y), and
direction, , of the particle can be modeled as a set of independent advection equations
acting in the (x; y) dimension coupled in the  dimension by the diusion equation [21].
Unfortunately, solutions of this Fokker-Planck equation computed by numerical integration
on a rectangular grid do not exhibit the robust invariance under rotations and translations
which characterizes the output of computations performed in primary visual cortex. Nor
does any other existing model of contour completion, sharpening, or saliency (e.g., [11, 12,
Our new algorithm computes stochastic completion elds in a Euclidean invariant manner

Figure

2 (left) is a picture of the stochastic completion eld due to the Kanizsa Triangle
stimulus in Figure 1(b). Figure 2 (right) shows the stochastic completion eld due to a
rotation and translation of the (input) Kanizsa Triangle. The Euclidean invariance of our
algorithm can be seen by observing that the (output) stochastic completion eld on the right
in

Figure

2 is itself a rotation and translation of the stochastic completion eld on the left,
by the same amount.
An important fact which has become clear to us is that Euclidean invariance is a property
only of computations dened in the continuum, and can be achieved in discrete computations
only when the discrete computations implement Euclidean invariant continuous computa-
tions. This statement might seem contradictory, but we will see that it is sometimes the case
that a Euclidean invariant computation in the continuum can be implemented using a nite
number of operations on a nite number of basis functions. Done properly, the resulting
discrete computation preserves the Euclidean invariance of the computation in the contin-
uum. 2 The general approach we adopt in this paper involves explicitly using the shiftability
and steerability of the basis functions to transform a Euclidean invariant computation in the
continuum into a discrete computation on a lattice. The discrete computation operates on
the coecients of a wavelet-like transform of the function to be computed. We will see that,
as is the case with the solution of the 2D advection equation discussed in this paper, a Euclidean
invariant discrete implementation of a continuous computation can be signicantly
dierent than the discrete computation which results from straightforward subsampling of
functions and substituting sums for integrals and dierences for derivatives.
neuroscience
Our new Euclidean invariant algorithm was motivated, in part, by the following experimental
ndings. To begin with, the receptive elds of simple cells, which have been traditionally
described as edge (or bar) detectors, can be accurately modeled using two-dimensional Gabor
functions [6, 19], which are the product of a Gaussian (localized in position) and a harmonic
grating (localized in orientation and spatial frequency). Gabor functions are well suited to
the purpose of encoding visual information, since, by the Heisenberg Uncertainty Principle,
they are the unique functions which are maximally localized in both space and frequency.
The sampling of the visual eld in V1 is quite sparse|there are about about 100  100
hypercolumns, with receptive elds of about 5 scales and 16 orientations in each hypercol-
umn. Neglecting size (and phase), a simple cell receptive eld can be parameterized by
2 It is important to note that simply representing the states of a discrete computation using shiftable-
steerable functions will not, in general, yield a Euclidean invariant computation. Discrete implementations
of continuous computations must be carefully designed to preserve Euclidean invariance.

Figure

2: Stochastic completion elds. Left: Of Kanizsa Triangle. Right: After the initial conditions
have been rotated and translated.
its position and orientation. The spatial distribution of these two parameters, known as
orientation preference structure, is an attempt (on the part of evolution) to smoothly map
the three-dimensional parameter space, R 2  S 1 , of edge positions and orientations onto
the two-dimensional surface, of the visual cortex, R 2 . Due to the dierences in dimension-
ality, orientation preference structure is punctuated by so-called pinwheels, which are the
singularities in this mapping [1].
As a rst approximation, a neuron's response to an arbitrary grey-level image can be
modeled as the L 2 -inner product of the image with the neuron's receptive eld. These
experimental observations suggested to Daugman [7] (and others, e.g. [23, 17]) that an
ensemble of simple cell receptive elds can be regarded as performing a wavelet transform
of the image, in which the responses of the neurons correspond to the transform coecients
and the receptive elds correspond to the basis functions.
Recent experiments have demonstrated that the response of simple cells in V1 can be
modulated by stimuli outside the classical receptive elds. Apart from underscoring the
limitations of the classical (linear) model, they suggest a function for the long-range connections
which have been observed between simple cells. For example, in a recent experiment,
Gilbert [10] has demonstrated that a short horizontal bar stimulus can modulate the response
of simple cells whose receptive elds are located at a signicant horizontal distances from
the bar, and which have a similar orientation preference to the bar. Non-linear long-range
eects have also been observed in secondary visual cortex. For example, von der Heydt et
al.[13] reported that the ring rate of certain neurons in V2 increases when their \receptive
elds" are crossed by illusory contours (of specic orientations) which are induced by pairs
of bars
anking the receptive eld. Signicantly, these neurons do not respond to these same
bars presented in isolation|they only respond to pairs. 3
3 These experiments suggests that the source and sink elds, which are intermediate representations in
Williams and Jacobs model of illusory contour formation, could be represented by populations of simple cells
Although our new contour completion algorithm does not provide a model for illusory
contour formation in the brain which is realistic in every respect, it does have several features
which are biologically plausible, none of which are found in previous algorithms, e.g.,
[11, 12, 18, 22, 24, 30, 31, 34]. These features are that (1) all states of the computation be represented
in a wavelet-like basis of functions which are localized in both space and frequency
(spatial localization allows the computation to be performed in parallel); (2) the computation
operates on the coecients in the wavelet-like transform and can be implemented in a
neural network; (3) the computation is Euclidean invariant; and (4) it is accomplished using
basis functions with centers lying on a (relatively) sparse grid in the image plane.
Shiftable-twistable bases
Many visual and image processing tasks are most naturally formulated in the continuum
and are invariant under a group of symmetries of the continuum. The Euclidean group,
of rotations and translations, is one example of a continuous symmetry group. However,
because discrete lattices are not preserved by the action of continuous symmetry groups,
the natural invariance of a computation can be easily lost when it is performed in a discrete
network. In this section we will introduce the notion of a shiftable-twistable basis and show
how it can be used to implement discrete computations on the continuous space of positions
and directions in a way which preserves their natural invariance.
In image processing, the input and output are functions on R 2 , and the appropriate notion
of the invariance of computations is Euclidean invariance|any rotation and translation of
the input should produce an identical rotation and translation of the output. Simoncelli
et al. [26, 9] introduced the notion of a shiftable-steerable basis of functions on R 2 , and
showed how it can be used to achieve Euclidean invariance in discrete computations for
image enhancement, stereo disparity measurement, and scale-space analysis.
Given the nature of simple cell receptive elds, the input and output of computations in
primary visual cortex are more naturally thought of as functions dened on the continuous
space, R 2  S 1 , of positions, in the plane, R 2 , and directions, , in the circle,
. For such computations the appropriate notion of invariance is determined by those
symmetries, T ~x
, of R 2  S 1 , which perform a shift in R 2 by ~x 0 , followed by a twist in
through an angle,  0 . A twist through an angle,  0 , consists of two parts: (1) a
rotation, R 0 , of R 2 and (2) a translation in S 1 , both by  0 . The symmetry, T ~x
, which is
called a shift-twist transformation 4 , is given by the formula,
(R  0
A visual computation on R 2  S 1 is called shift-twist invariant if, for all (~x
a shift-twist of the input by (~x produces an identical shift-twist of the output.
Correspondingly, we dene a shiftable-twistable basis 5 of functions on R 2  S 1 to be a
set of functions on R 2  S 1 with the property that whenever a function, P (~x; ), is in their
in V1, and that the stochastic completion eld, which is the product of the source and sink elds, could be
represented in V2.
4 The relationship between shift-twist transformations and computations in V1 was described by Williams
and Jacobs in [30] and more recently by Kalitzin et al. [16] and Cowan [3].
5 We use this terminology even though the basis functions need not be linearly independent.
span, then so is P (T ~x
(~x; )), for every choice of (~x As such, the notion of
a shiftable-twistable basis on R 2  S 1 generalizes that of a shiftable-steerable basis on R 2 .
Shiftable-twistable bases can be constructed as follows. First we recall Simoncelli's concept
of the shiftability of a function, which is closely related to the Shannon-Whittaker
Sampling Theorem. A periodic function, (x), of period X, is shiftable if there is an inte-
ger, K, such that the shift of by an arbitrary amount, x 0 , can be expressed as a linear
combination of K basic shifts of , i.e., if there exist interpolation functions, b k
that
is the basic shift amount. The simplest shiftable function in one dimension
is a pure harmonic signal, e i!x , in which case K = 1. More generally, Simoncelli et al. [26]
proved that any band-limited function is shiftable. In fact, if the set of non-zero Fourier series
frequencies of is (a subset of) then can be shifted using
the K interpolation is the complex conjugate of
the perfect bandpass lter constructed from the set of K frequencies, B. In particular, note
that the interpolation functions only depend on the set of non-zero frequencies of , and not
on itself.
Strictly speaking, since they are not band-limited, functions such as Gabors are not
shiftable. Nevertheless, for all intents and purposes, they can be shifted by choosing the
set, B, to consist of all Fourier series frequencies, !, of , such that the Fourier amplitude,
(!)j, exceeds some small threshold value. Such functions will be called eectively shiftable.
Let (~x; ) be a function on R 2  S 1 which is periodic (with period X) in both spatial
variables, ~x. In analogy with the denition of a shiftable-steerable function on R 2 , we say
that is shiftable-twistable on R 2  S 1 if there are integers, K and M , and interpolation
the shift-twist of by (~x
is a linear combination of a nite number of basic shift-twists of by amounts ( ~ k; m  ),
i.e., if
~ k;m
Here is the basic shift amount and  is the basic twist amount. The sum
in equation (2.3) is taken over all pairs of integers, ~
and all integers, m, in the range, 0  m < M . As we will show, for many shiftable-
twistable bases, the interpolation functions, b ~ k;m (~x are dened in terms of Simoncelli's
one-dimensional interpolation functions, b k
The simplest shiftable-twistable functions are those which can be twisted with
basic twists: (R 0 (~x);   0 functions will be called self-twistable.
The following Proposition shows how to construct a shiftable-twistable basis from a single
shiftable-twistable function.
Proposition 2.1 Let be a shiftable-twistable function with interpolation functions, b ~ k;m (~x
Then the collection of functions, ~ k;m , dened by
where  and   are the basic shift and twist amounts, form a shiftable-twistable basis. More
precisely, if
~ ';n
c ~ ';n ~ ';n (~x; ); (2.5)
then
(R
~ k;m
where
~ l;n
The proof of the Proposition is a straightforward application of the following composition
rule for shift-twist transformations:
Equation (2.8) can be explained as follows. The left hand side is a formula for the action of
the shift-twist, T ~x
, on T ~x
(R  0
and then rotates the result by  1 to obtain
(R  1
(R  0
This is the same as the single shift-twist transformation on the right hand side of Equation
(2.8), which shifts ~x by ~x 0 +R 0 (~x 1 ) and then rotates the result by  . Details of
the proof of Proposition 2.1 can be found in the Appendix.
This paper is concerned with visual computations whose input and output are functions
on the continuum, R 2  S 1 , and which are shift-twist invariant. We propose the following
general framework for performing such computations in a shift-twist invariant manner in a
discrete network. First, the various states of the computation, P (~x; ), are to be expressed
in a shiftable-twistable basis, ~ k;m , as
~ k;m
c ~ k;m ~ k;m (~x; ); (2.10)
is the state or coecient vector. The input state vector is transformed to
the output state vector by a feedforward or recurrant neural network. This is an appropriate
general framework for biologically plausible visual computation.
The biological plausibility of a given computation also depends on the specic choice
of shiftable-twistable basis. We conclude this section by presenting several examples of
shiftable-twistable bases, in order of increasing biological plausibility and increasing com-
plexity. In each case, the basis functions are separable, i.e., they are the product of a
periodic shiftable-steerable function, (~x), on R 2 and a shiftable function, f(), on
Example A. The Gaussian-Fourier basis
The Gaussian-Fourier basis, G ~ k;! , is the product of a shiftable-steerable basis of Gaussians
in ~x and a Fourier series basis in . Let
be a radial Gaussian of standard
deviation, . We regard g as a periodic function of period, X, which is chosen to be much
larger than , so that g( X
its derivatives are essentially zero. For each frequency, !,
Also, given a choice of a shift amount, , so that
is an integer, we dene the Gaussian-Fourier basis functions, G ~ k;! , by
G ~ k;! (~x;
To show that G ! is shiftable-twistable, we argue as follows. First, G ! is self-twistable
since the pure harmonic, e i! , is self-shiftable, e i!( and the Gaussian, being
isotropic, is self-steerable, g(R  0
separable functions (such as Gaussians)
are shiftable if and only if each of their factors, shiftable. Furthermore,
the interpolation function of a separable function is the product of the interpolation functions
of its factors, b ~ k (~x 0
(y 0 ). These observations immediately imply the following
result:
Proposition 2.2 The periodic function, G ! , (of period X) is eectively shiftable-twistable.
More precisely, let let K be the number of essentially non-zero Fourier series
coecients of the factor, g X
, of g(~x). Then,
where the interpolation functions are given by
Here b ~ k (~x 0
~ k), where b(~x 0 ) is the complex conjugate of the perfect bandpass lter
constructed from the set of K 2 essentially non-zero Fourier series coecients, ~, of g(~x).
For the Gaussian-Fourier basis, G ~ k;! , because of the isotropy of the Gaussian, Proposition
2.1 takes the following special form.
Corollary 2.3 The Gaussian-Fourier basis, G ~ k;! (~x; ), dened by equation (2.11), is shiftable-
twistable. More precisely, if
~ ';!
then
(R  0
~ k;!
where
In particular, if
each frequency, !, the transformed coecient vector, c ! (~x is the circular
convolution of the original coecient vector, c g, with the interpolation function
vector,
Although it is not as biologically plausible as the other bases we describe, the Gaussian-
Fourier basis can be a very ecient basis to use in computer vision because of equation (2.17)
and because in certain cases the input can be easily represented in the basis. For example,
suppose that the input is modeled as a linear combination of ne scale three-dimensional
Gaussians, centered at arbitrary points, (~x Since the input is the product
of a Gaussian in ~x and a Gaussian in  it can be represented in a single scale Gaussian-
Fourier basis as follows. First, the Gaussian in  is represented in the Fourier basis using the
standard analysis and synthesis formulae for Fourier series. Second, if the two-dimensional
input Gaussians in ~x are chosen to be shifts of the basis function, g(~x), then we can use
Proposition 2.2 to represent the input Gaussians in ~x in the Gaussian basis.
Example B. The complex directional derivative of Gaussian (CDDG){Fourier
basis
This example is very similar to the previous one, except that the Gaussian, g(~x), is replaced
by its complex directional derivative in the direction of the complex valued vector, [1; i] T .
We use the notation, (~x), to represent this directional derivative. Whereas g(~x) is not a
wavelet, (~x) is. The formula for (~x) is
@x
@g
@y
As in Example A we regard as a periodic function of period, X. Pictures of (~x) are
shown in Figure 3. The complex directional derivative of Gaussian (CDDG){Fourier basis,
~ k;! , is dened by
Freeman and Adelson [9] showed that the complex directional derivative of a Gaussian
is self-steerable: (R  0
Similarly the function, (~x;
twistable: (R 0 (~x);   0 As such, is the simplest non-isotropic
shiftable-twistable function which is localized in the spatial variables, ~x.
The analogies of Proposition 2.2 and Corollary 2.3 hold for the basis, ~ k;! . The only
dierence is that in this case, the interpolation functions in equations (2.13) and (2.16) are
given by

Figure

3: The real part of the complex directional derivative of a Gaussian (CDDG), steered
clockwise by 45  . Left: Real part with standard deviation, (with 40:0  40:0 display
region). The imaginary part is the rotation of the real part, counter clockwise by 90  . Right: Real
part at the ner scale, 0:25. The completion elds in Figure 2 were computed by solving the
Fokker-Planck equation in a CDDG-Fourier basis of this scale.
One property of the CDDG-Fourier basis is that the function, (~x), (which is the -
marginal of 0 ) looks much more like the receptive eld of a simple cell in V1 than a
Gaussian does. Unlike the Gaussian-Fourier basis, the CDDG-Fourier basis has the property
that its spatial factor, (~x), is a wavelet, which means that arbitrary input functions can be
represented in a multi-scale shiftable-twistable basis constructed from the functions, ! .
The multi-scale method of representing the input function exploits the fact that the
function, (~x), is a mother wavelet which generates an (approximately) self-inverting, over-complete
wavelet basis [26] (i.e., an approximately tight frame [4]). Each of the 2D wavelet
basis functions, p; ~ k (~x), is dened to be the translation by 2 p=2  ~ k of the scaling of by 2 p=2
(i.e., two voices per octave),

where the (scale dependent) period, X, of the wavelet, p; ~ k , is chosen so that
and  is chosen so that each scale of the wavelet-Fourier basis, p; ~ k;! (~x; ), forms a shiftable-
twistable basis. Even though the basis functions are not linearly independent, a large class
of functions, P (~x; ), can be well approximated in the basis using the analysis and synthesis
formulae,
c
ZZZ
c
Example C. The complex directional derivative of Gaussian (CDDG){Gaussian
coupled basis
The complex directional derivative of a Gaussian, (~x), is a simple edge detector. Since the
real part of has the preferred orientation, , an edge centered at ~x 0 with orientation,  0 ,
will elicit a large response from the real part of the shift-twist of by (~x
problems in vision require that the direction or orientation of edges in an image, I(~x), be
explicitly encoded in a third variable, . This can be done using a coupled basis consisting of
complex directional derivatives of Gaussians in ~x and Gaussians in , where the center of the
Gaussian in  is coupled to the preferred orientation of the complex directional derivative of
Gaussian.
be a Gaussian whose standard deviation, , is chosen so that g
can be regarded as a periodic function of period, 2. Let M be the number of essentially
non-zero frequencies of g, and let
M be the basic twist amount. Let (~x) be the complex
directional derivative of Gaussian dened in equation (2.18). The complex directional
derivative of Gaussian (CDDG)-Gaussian coupled basis functions are dened by
(~x ~ k).
Given an image, I(~x), the function
~ k;m
explicitly encodes the position and direction of edges in I, where
c ~
ZZ
the interpolation functions required to shift g() and (~x)
respectively. Then the functions, ~ k;m , form a shiftable-twistable basis with interpolation
functions 6
Also observe that the -marginal of ~ k;m (~x; ) is ~ k (Rm
(~x)), which is similar in shape
to the prole of a simple cell receptive eld in V1.
Example D. The Gabor-Gaussian basis
The two-dimensional Gabor function is often used to model two-dimensional simple cell receptive
elds in V1. Its formula is is a radial
Gaussian and e i 0 (x+y) is a harmonic grating of frequency,  0 . Like the complex directional
6 The factor e i(0 m  ) steers (R 0 (~x)) to (Rm
derivative of Gaussian, the Gabor possesses a preferred direction, and so one might attempt
to construct a coupled Gabor-Gaussian shiftable-twistable basis. Unfortunately, such
a coupled basis in not shiftable-twistable since, unlike the complex directional derivative
of Gaussian, the Gabor function is not self-steerable. 7 Nevertheless the uncoupled basis,
shiftable-twistable since the Gabor functions, ~ k;m (~x), form
a shiftable-steerable (multiscale wavelet) basis and the Gaussians, g n (), form a shiftable
basis.
3 Stochastic completion elds
In their computational theory of illusory contour formation, Williams and Jacobs [30] argued
that, given a prior probability distribution of possible completion shapes, the visual system
computes the local image plane statistics of the distribution of all possible completions, rather
than simply the most probable completion. This view is in accord with human experience|
some illusory contours are more salient than others, and some appear sharper than others.
They dened the notion of a stochastic completion eld to model illusory contours in a
probabilistic manner. The stochastic completion eld is a probability density function (p.d.f.)
on the space, R 2  S 1 , of positions, in the plane, R 2 , and directions, , in the
circle, S 1 . It is dened in terms of a set of position and direction constraints representing
the beginning and ending points of a set of contour fragments (called sources and sinks), and
a prior probability distribution of completion shapes, which is modeled as the set of paths
followed by particles traveling with constant speed in directions described by Brownian
motions [21]. The magnitude of the stochastic completion eld, C(~x; ), is the probability
that a completion from the prior probability distribution will pass through (~x; ) on a path
joining two of the contour fragments. Williams and Jacobs [30] showed that the stochastic
completion eld could be factored into a source eld and a sink eld. The source eld,
represents the probability that a contour beginning at a source will pass through
(~x; ) and the sink eld, Q 0 (~x; ), represents the probability that a contour beginning at (~x; )
will reach a sink. The completion eld is
The source (or sink) eld itself is obtained by integrating a probability density function,
positive times, t, where P (~x;  ; t) represents the probability that a particle
beginning at a source reaches (~x; ) at time t,
Mumford [21] observed that P evolves according to a Fokker-Planck equation of the form,
@t
@x
sin
@y
7 We speculate that a coupled basis whose basis functions are similar in shape to Gabors in ~x and Gaussians
in  exists, but that the basis functions are not separable, i.e., they are not of the form, (~x)g().
where the initial probability distribution of sources (or sinks) is described by P (~x;  ; 0). This
partial dierential equation can be viewed as a set of independent advection equations in
(the rst and second terms) coupled in the  dimension by the diusion equation
(the third term). The advection equations translate probability mass in direction  with
unit speed, while the diusion term models the Brownian motion in direction, with diusion
parameter, . The combined eect of these three terms is that particles tend to travel in
straight lines, but over time they drift to the left or right by an amount proportional to  2 .
Finally, the eect of the fourth term is that particles decay over time, with a half life given by
the decay constant,  . This represents our prior expectation on the length of gaps|most are
quite short. In [31] stochastic completion elds were computed by solving the Fokker-Planck
equation using a standard nite dierencing scheme on a regular grid.
In [30, 31] the initial sources and sinks were extracted automatically from the input
image using steerable lters [25] to identify corners and measure orientations. However, this
method should only be regarded as an interim solution to the problem of specifying initial
conditions given a brightness image, since we believe that in the long term, the solution lies
in the generalization of the notion of a completion eld described by Thornber and Williams
in [28]. For this reason and for convenience's sake, in all experiments in this paper, the
initial sources and sinks were specied by hand.
4 Description of algorithm
One of the main goals of this paper is to derive a discrete numerical algorithm to compute
stochastic completion elds in a shift-twist invariant manner. This invariance is achieved
by rst evolving the Fokker-Planck equation in a shiftable-twistable basis of R 2  S 1 to
obtain representations of the source and sink elds in the basis, and then multiplying these
representations in a shift-twist invariant manner to obtain a representation of the completion
eld in a shiftable-twistable basis.
We observe that a discrete Dirac basis, consisting of functions, ~ k;m (~x;
~ m) is a triple of integers, is not shiftable-twistable. This is
because a Dirac function located o the grid of Dirac basis functions is not in their span.
A major shortcoming of all previous contour completion algorithms [11, 12, 18, 22, 24, 30,
31, 34] is that they perform computations in this basis. As a consequence, initial conditions
which do not lie directly on the grid cannot be accurately represented. This problem is
often skirted by researchers in this area by choosing input patterns which match their choice
of sampling rate and phase. For example, Li [18] used only six orientations (including 0  )
and Heitger and von der Heydt [12], only twelve (including 0  , 60  and 120  ). Li's rst
test pattern was a line of orientation, 0  , while Heitger and von der Heydt used a Kanizsa
Triangle with sides of 0  , There is no reason to believe that the
experimental results they show would be the same if their input patterns were rotated by as
little as 5  . 8
8 Nor are we blameless in this respect. Williams and Jacobs [30, 31] used 36 directions (including 0  ,
and 120  ) and demonstrated their computation with a Kanizsa Triangle with sides of 0  , 60  and 120
orientation.
In addition to the problem of representing the input, the computation itself must be
Euclidean invariant. Stochastic completion elds computed using the nite dierencing
scheme of [31] exhibit marked anisotropic spatial smoothing due to the manner in which 2D
advection is performed on a grid (see Figures 7,8 and 9). Although probability mass advects
perfectly in either of the two principal coordinate directions, mass which is moving at an
angle to the grid gradually disperses, since, at each time step, bilinear interpolation is used to
place the mass on the grid. One way to restore the isotropy of the advection transformation
is to carefully add extra spatial smoothing [33]. Higher order nite dierencing schemes have
also proved extremely eective [2].
Stochastic completion elds can be computed using any of the shiftable-twistable bases
in Section 2. For reasons of simplicity, in this paper, we chose to perform the computation in
the Gaussian-Fourier basis. The initial conditions for the Fokker-Planck initial value problem
are modeled by ne scale, three-dimensional Gaussians, whose centers are determined by the
locations and directions of the edge fragments to be completed. We use the single scale
method discussed in Example A of Section 2 to represent the initial conditions in the basis.
To solve the Fokker-Planck equation, we express its solution in terms of the basis func-
tions, G ~ k;! (~x; ), as
~ k;!
c ~ k;! (t) G ~ k;! (~x; ); (4.1)
where the coecients, c ~ k;! (t), depend on time. Then, we derive a linear transformation,
to evolve the coecient vector in time. This transformation is
the composition of an advection transformation, A, which has the eect of transporting
probability mass in directions , and a diusion-decay transformation, D, which implements
both the diusion of mass in , and the decay of mass over time. Representations of source
or sink elds in the basis are obtained by integrating the coecient vector, c(t), over time,
where the initial coecient vector represents the initial sources or sinks.
The shiftability-twistability of the basis functions is used in two distinct ways to obtain
shift-twist invariant source and sink elds. First, it enables any two initial conditions, which
are related by an arbitrary transformation, T ~x
, to be represented equally well in the ba-
sis. Second, it is used to derive a shift-twist invariant advection transformation, A, thereby
eliminating the grid orientation artifacts described above. In summary, given a desired resolution
at which to represent the initial conditions, our new algorithm produces source and
sink elds, at the given resolution, which transform appropriately under arbitrary Euclidean
transformations of the input image. In contrast, in all previous contour completion algo-
rithms, the degree of failure of Euclidean invariance is highly dependent on the resolution of
the grid, and can be quite large relative to the grid resolution.
The nal step in our shift-twist invariant algorithm is to compute the completion eld
(the product of the source and sink elds) in a shiftable-twistable basis. The particular basis
used to represent completion elds is the same as the one used to represent the source and
sink elds, except that the variance of the Gaussian basis functions in R 2 needs to be halved.
The need to use a slightly dierent basis to represent completion elds is not biologically
implausible, since the experimental evidence described in Section 1 suggests that the neural
locus of the source and sink elds could be V1, while completion elds are more likely located
5 The solution of the Fokker-Planck equation
In this section we derive a shift-twist invariant linear transformation,
of the coecient vector which evolves the Fokker-Planck equation in a shiftable-twistable
basis. The derivation holds for any shiftable-twistable basis constructed from shiftable-
twistable functions of the form, ! (~x; (~x). Since the
transformation, A  D, will only involve interactions between functions, (~x), at dierent
positions ~ k, and not at dierent scales or orientations, the basis functions and coecients
will be denoted by ~ k;! and c ~ k;! (t) respectively.
To derive an expression,
~ k;!
for the advection transformation, A, in the basis, ~ k;! , we exploit the fact that spatial advection
can be done perfectly using shiftable basis functions, ~ k (~x), in R 2 , and the continuous
variable,  2 S 1 . Suppose that P is given in the form,
~ k;!
c ~ k;! (t) ~ k (~x) e
c ~ k; (t) ~ k (~x); (5.2)
where  c(t) is related to c(t) by the standard synthesis formula for Fourier series,  c ~
which we denote by  the translation of P in direction, , at unit
speed, for time, t, is given by
where the second equation follows from equation (5.2). The shiftability of then implies
that
A ~ ';; ~ k; (t)  c ~ k; (t); (5.5)
where
A ~ ';; ~ k;
Finally, the advection transformation, A, in the basis, ~ k;! , is given by the similarity
denotes the standard analysis formula for Fourier
series,
R 2f()e i! d. Since we have the following result.
Theorem 5.1 In the basis, ~ k;! , the advection transformation, A, is given by
~ k;!
where
In particular, the transformation, A, is shift-twist invariant and is a convolution operator
on the vector space of coecients, c ~ k;! .
The proof of the shift-twist invariance of A is given in the Appendix.
The degree of accuracy with which the discrete advection transformation, A, models the
continuous advection process is determined by the number of basis functions in the shiftable-
twistable basis. In order to prevent aliasing in the representation of P in the basis, the set
of Fourier series frequencies, !, must be large enough to capture the -frequency content of
the p.d.f., P (~x; ; t). Suppose that the initial p.d.f. is the Gaussian-Fourier basis function,
. Then the advection, P
of P can be described as follows. For each xed angle, , the function, g(~x) is translated in
R 2 by t units in the direction, , and weighted by the factor, e i! . Consequently, at time, t,
the p.d.f. is supported on a neighborhood of a helix 9 of radius, t, oriented about the -axis
in R 2  S 1 .
Consider the case that the initial p.d.f., P (~x; ; 0), is the Gaussian-Fourier basis function,
e 12i (of period where the basis consists of
shifts in the x and y variables and 176 frequencies, !, in the  variable.

Figure

4 (left) shows the -integral of the advection of P (~x; ; 0) at time,
R 2P (~x; ; 14:0) d. In particular, note that
R 2P (~x; ; 14:0) d is supported on a circle of
radius 14.0, which is the projection from R 2 S 1 to R 2 of the helical support of P (~x; ; 14:0).
Coupling the -diusion process to the advection process decreases the -frequency content
of the p.d.f., P (~x; ; t), and so fewer Fourier series frequencies are required in the basis.
For example, if the diusion parameter is then the number of -frequencies can be
reduced from 176 to 92. Figure 4 (right) shows the -integral of the source eld obtained by
evolving the initial p.d.f., G ~ 0;12
according to the Fokker-Planck equation, with diusion
parameter, decay constant,
The -frequency content of the helical p.d.f. and the accuracy of the advection process
are analyzed in the appendix. This analysis implies that if the spatial resolution of the
Gaussian initial conditions is increased by a factor of two in all three variables, (~x; ), then
the number of basis functions must be multiplied by 16 to ensure the same degree of accuracy
in the advection of the initial conditions.
The diusion-decay transformation, c(t which we use, is given in the
following Proposition, the proof of which is in the Appendix.
Proposition 5.2 Let N be the number of Fourier series frequencies, !, used in the shiftable-
twistable basis, and let . The diusion-decay transformation, D, is given by
c ~ k;! (t
. Furthermore, D is shift-twist invariant. 10
9 The parameterization of the helix is 2.
If the transformations, A and D, are shift-twist invariant, then so is A  D.

Figure

4: Left: The -integral of the advection of the Gaussian-Fourier basis function, G ~ 0;12 (~x; ),
at time Right: The -integral of the source eld, P
R 14:0
is the solution of the Fokker-Planck equation, with initial conditions given by the basis
The equation (5.9) is the transfer function for the
diusion equation, and the factor, e t= implements the decay. If c ~ k;! (t) evolves according
to equation (5.9), then P (~x;  ; t) evolves according to the standard, explicit, 3-point stencil
nite dierencing scheme for diusion in  and exponential decay in time. In particular, if
t is chosen so that   0:5, then this nite dierencing scheme is stable.
Theorem 5.1 and Proposition 5.2 imply that the computation of source and sink elds
can be performed in a recurrent neural network using a xed set of units as described in [31].
Furthermore, the resulting source and sink eld coecient vectors only depend on the initial
data and on the set of (essentially) non-zero Fourier series frequencies of the basis function,
(~x).
Since the advection transformation, A, is a convolution operator on the space of coe-
cients, for eciency's sake we implemented both A and D in the 3D Fourier domain of the
coecient vector. In this domain, A is given by a diagonal matrix and D by a circulant
tridiagonal matrix.
6 Completion elds in the basis
The representation of completion elds in the Gaussian-Fourier basis uses basis functions
on R 2 which have half the variance of those used to represent source and sink elds. 11 The
be used to refer to the ner scale basis. Set e
2, e
2,
e
2K and e
. Then the corresponding
11 The reason we need to use ner scale basis functions is that the product of the Gaussian, e x 2 =2 , with
itself is the ner scale Gaussian, e x 2
Gaussian-Fourier basis, e
G ~ j; , is shiftable-twistable, with interpolation functions, e b ~ j (~x
Recall that the completion eld, C(~x; ), is the product of the source and sink elds,
the expressions for the source and sink elds in the Gaussian-
Fourier basis be denoted by
~ k;!
and
where the source and sink coecient vectors, p 0
~ k;!
and q 0
, are computed from the initial
source and sink coecient vectors, p ~ k;! (0) and q ~ '; (0), by integrating the solution coecient
vectors, p ~ k;! (t) and q ~ '; (t), of the Fokker-Planck initial value problem over time.
To represent the completion eld in the Gaussian-Fourier basis, e
G ~ j; , it suces to express
the product, G ~ k;! G ~ '; , of two basis functions in the basis, e
G ~ j; . In the Appendix we prove
the product formula 12
G ~ k;! G ~
2e
e
Combining this product formula with equations (6.1) and (6.2) gives the following Theorem.
Theorem 6.1 The expression for the completion eld in the Gaussian-Fourier basis is given
by
e
where the completion eld coecient vector, C ~ j; , is given by
2e
~ k; ~ '
~ k;
e
~ k;!
and
are the source and sink eld coecient vectors.
In particular, the calculation of C ~ j; from the initial source and sink coecient vectors,
can be performed in a recurrent neural network using a xed set of units
as described in [31].
12 The interpolation functions e b ~ j are evaluated at e
since G ~ k;! G ~ '; is a Gaussian
centered at ( ~ k
It is helpful to observe two features of equation (6.5). First, the presence of the Gaussian
involving k ~ k ~ 'k means that in practice, the sum can be taken only over those indices
~ k and ~ ' for which k ~ k ~ 'k is quite small, i.e., only spatially proximal pairs of source and
sink eld coecients need to interact to compute the completion eld coecient vector.
Second, the circular convolution of p 0
~ k;
and q 0
in the frequency variable, , corresponds to
multiplication in the  variable, which occurs when the completion eld is constructed from
the source and sink elds.
7 Completion elds in the CDDG{Fourier basis
The algorithm we have presented for computing stochastic completion elds in the Gaussian-
Fourier basis can be interpreted as a computation in the more biologically plausible complex
directional derivative of Gaussian (CDDG)-Fourier basis. As in the Gaussian-Fourier basis, a
single scale basis, ~ k;! (~x; ), is used to represent the initial conditions and solve the Fokker-Planck
equation.
Rather than solving the Fokker-Planck equation with Gaussian initial conditions, we
take a directional derivative in the spatial variables, ~x, of the initial conditions, (as in
equation (2.18)) to obtain new initial conditions which are sums of functions of the form,
shift of the basis function, (~x), and g is a Gaussian.
The shiftability of the basis functions, ~ k;! , can be used to accurately represent the initial
conditions in the basis. Using the directional derivatives of the initial sources 13 as initial
conditions, a source eld coecient vector, p 0
~ k;!
, is obtained by solving the Fokker-Planck
equation in the basis. In this manner we obtain a representation of a derivative of the source
eld in the CDDG-Fourier basis:
~ k;!
~ k;!
where D is the operator,
@y , in equation (2.18). A representation of the source
eld itself can be obtained by integrating equation (7.1) as follows. Let G ~ k;! (~x; ) be the
Gaussian-Fourier basis function dened in Equation (2.11). By equation (2.18), ~
Therefore, by equation (7.1), the source eld can be represented in the Gaussian-
Fourier basis as 14
~ k;!
~ k;! G ~ k;! (~x;
13 The idea of representing the initial conditions by their derivatives is not biologically implausible. Indeed
wavelet basis functions are more suited to representing the derivatives of the initial conditions than the
initial conditions themselves. This is due to the fact that wavelets have average value zero. Consequently,
the initial probability density functions, which are inherently positive, cannot be accurately represented in a
wavelet basis as positive functions, whereas their derivatives, like the wavelets used to represent them, have
average value zero.
14 The rationale for this method of calculating the source eld is that, for each xed , the Fokker-Planck
equation is a constant coecient linear equation.
Note that the two source eld coecient vectors, p 0
~ k;!
, obtained by solving the Fokker-Planck
equation in the two bases, ~ k;! and G ~ k;! are identical, since the frequency content of G ! and
are essentially the same and since the source eld coecient vector, p 0
~ k;!
, only depends
on the interpolation functions, b ~ k (~x 0 ;  0 ), and on the locations and directions of the source
initial conditions.
Finally, the expression for the directional derivative, D C, of the completion eld in the
CDDG-Fourier basis is given by
e
where C ~ j; is given by equation (6.5).
8 Experimental results
We present four experiments demonstrating the Euclidean invariance of our algorithm. In
each experiment, the Gaussian-Fourier basis consisted of in each spatial
variable of a Gaussian (of period signals of
in the angular variable, for a total of 2:355 functions. Pictures of completion
elds were obtained by analytically integrating over  and rendering the completion eld on
a 256  256 grid.
In the rst experiment, we computed source elds using the new algorithm. The diusion
parameter was 0:12, the decay constant, and the time increment,
0:1. The left column of Figure 5 shows the p.d.f. and source eld due to a single three-dimensional
Gaussian initial condition centered at (~x 0 while in the
right column the initial condition has been shifted and twisted to be centered at (~x 0
The p.d.f.'s are shown at time, and the source elds were
obtained by integrating the p.d.f. up to time, 30:0. The source elds were clipped above
at
In the remaining experiments, we compare the new algorithm with the nite dierencing
scheme of [31]. For the nite dierencing scheme, the 40:0  40:0  2 space was discretized
using a 256256 spatial grid with 36 discrete orientations, for a total of 2:35910 6 Dirac basis
functions. The intent was to use approximately the same number of basis functions for both
algorithms. The initial conditions were represented on the grid using tri-linear interpolation
and pictures of the completion elds were obtained by summing over the discrete angles.
The same parameters were used for both algorithms. The decay constant was
the time increment, 0:1. The diusion parameter was for the second and
third experiments and for the last. 15 In Figures 7, 8 and 9 the completion elds
constructed using the algorithm of [31] are in the left column, while those constructed using
the new algorithm are in the right column.
In the second experiment, we computed straight line completion elds joining two diametrically
opposed points on a circle of radius, 16:0, with initial directions normal to
15 The diusion parameter, , was required to be larger in the last experiment because of the high curvature
circles in the Kanizsa triangle.
the circle. That is, given an angle, , the initial stimulus consisted of the two points,
(left). The completion elds are shown in Figure 7,
with those in the left column, computed using the method of [31], clipped above at 210 6 .
To compare the degree of Euclidean invariance of the two algorithms, we extracted a
section of each completion eld along the diameter of the circle normal to the direction of
the completion eld. In Figure 6 (right), we plot the mean of each section as a function of
the angle . The dashed line indicates the means computed using the new algorithm, and
the solid line shows the means computed using the algorithm of [31]. 16 The fact that the
dashed line graph is constant provides solid evidence for the Euclidean invariance of the new
algorithm. The solid line graph demonstrates the two major sources of the lack of Euclidean
invariance in the method of [31]. First, the rapid oscillation of period 10  is due to the initial
conditions coming in and out of phase with the angular grid. This 10  periodicity can be
seen in the periodicity of the general shape of the completion elds in the left column of

Figure

7. Second, the large spikes at 90  intervals are due to the anisotropic manner in which
the advection transformation was solved on the spatial grid. These large spikes correspond
to the very bright horizontal line artifacts in the rst two completion elds in the left column
of

Figure

7.
In the third experiment, we computed completion elds due to rotations of the Ehrenstein
initial stimulus in Figure 1(a). Pictures of the completion elds are shown in Figure 8. 17 The
top row shows the completion elds due to the Ehrenstein stimulus in Figure 1(a), while in
subsequent rows, the initial conditions have been rotated clockwise through angles,
. The completion elds computed using the method of [31] were clipped above
at . For our nal experiment, we compute completion elds due to rotations and
translations of the Kanizsa Triangle stimulus in Figure 1(b). Completion elds are shown in

Figure

2, which was discussed in the Introduction, and in Figure 9. The top row of Figure 9
shows completion elds due to the Kanizsa Triangle in Figure 1(b). In the third row the
initial conditions have been rotated clockwise by 5  . The second and fourth rows show the
regions inside the boxes in the rst and third rows, magnied times. The completion
elds computed using the method of [31] were clipped above at 9  10 5 .
The completion elds in the right columns of Figures 8 and 9, and in Figure 2, demonstrate
the Euclidean invariance of our new algorithm. This is in marked contrast with the
obvious lack of Euclidean invariance in the completion elds in the left columns of Figures 8
and 9. The visible straight line artifacts in these completion elds, which are oriented along
the coordinate axes, are due to the anisotropic nature of the advection process in the algorithm
of [31], and (to a lesser extent), to the way in which the initial conditions were
represented on the grid.
The angles, , were taken in 5  increments from 0  to 45  . For illustration purposes the -axis was
extended to 360  so as to re
ect the symmetry of the grid. Both graphs were normalized to have average
value one.
17 Because of the periodicity in the spatial variables, ~x, to avoid wrap around in this experiment, for the
new algorithm the computation was performed on a 80:0  80:0  2 space with

Figure

5: Probability density functions and source elds due to a single 3D Gaussian initial
condition centered at a point, (~x computed using the new algorithm. The p.d.f.'s are shown
at time, and the source elds are integrated out to time,
condition centered at (~x 0
R 30:0
eld with (~x
9 Future Directions
In a recent paper, Thornber and Williams [28] have generalized the denition of stochastic
completion eld in several important ways. Instead of computing the distribution of contours
joining a pair of position and direction constraints, they now compute the scale-invariant
distribution of closed contours satisfying a set of position constraints. The sources and
sinks for the stochastic completion eld, instead of serving as input to the computation,
are themselves the solutions of an eigenvector-eigenvalue problem. The eigenvector is the
xed-point of a continuous state, discrete time recurrent neural network. The kernel of the
weight matrix for the neural network is a solution of the Fokker-Planck equation.
One of our original motivations in developing shiftable-twistable functions was the need
to simulate the continuous state, discrete time neural network using a xed number of
parallel operations on a xed set of basis functions. One property of the solution in the
continuum, and of eigenvector-eigenvalue problems in general, is that relatively small matrix
values, through iteration, can eventually become the dominant factors in determining the
components of the eigenvector. We were aware of this fact and also aware of the fact that
the anisotropic error introduced by the straightforward method of solving the Fokker-Planck
section
source
sink

Figure

Left: The geometry of the straight line completion eld experiment. Right: Graph
of the mean along a section normal to the straight line completion eld as a function of the
direction . Dashed line: Our new algorithm. Solid line: The algorithm of [31].
equation [31] on a regular grid was on the order of 100 percent (see Figure 6). It was obvious
that without eliminating this huge relative error, that it would be impossible to accurately
simulate the continuous state, discrete time neural network. Now that we have succeeded in
our goal of developing a Euclidean invariant method for solving the Fokker-Planck equation
on a discrete grid, we are in a position to develop a local parallel network for computing the
generalization of the stochastic completion eld described in Thornber and Williams [28].
An important initial stage in the analysis of a scene requires completion of the boundaries
of partially occluded objects. Williams and Jacobs introduced the notion of the stochastic
completion eld which measures the probability distribution of completed boundary shapes
in a given scene. In this article we have described a new, parallel, algorithm for computing
stochastic completion elds. As is required of any computational model of human visual
information processing, our algorithm attempts to reconcile the apparent contradiction between
the Euclidean invariance of human early visual computations on the one hand, and
the observed sparseness of the discrete spatial sampling of the visual eld by primary and
secondary visual cortex on the other hand. The new algorithm reconciles these two contradictions
by performing the computation in a basis of separable functions with spatial
components similar to the receptive elds of simple cells in primary visual cortex. In partic-
ular, the Euclidean invariance of the computation is achieved by exploiting the shiftability
and twistability of the basis functions.

Figure

7: Straight line completion elds due to an initial stimulus consisting of two points
on a circle with direction  normal to the circle, with Using the

Figure

8: Completion elds due to rotations of the Ehrenstein initial stimulus in Figure 1(a).
From top to bottom, the initial conditions are rotated clockwise through angles,
15  and 45  . Left: Using the algorithm of [31]. Right: Using the new algorithm.

Figure

9: Top row: Completion elds due to the Kanizsa triangle initial stimulus in Figure
1(b). Left: Using the algorithm of [31]. Right: Using the new algorithm. Second row:
The regions inside the boxes, magnied 4. Third row: Initial conditions rotated clockwise
by 5  . Bottom row: The regions inside the boxes, magnied 4.
In this paper, we have described three basic results. First, we have generalized Simoncelli
et al.'s notion of shiftability and steerability in R 2 to a more general notion of shiftability
and twistability in R 2  S 1 . The notion of shiftability and twistability mirrors the coupling
between the advection and diusion terms in the Fokker-Planck equation, and at a deeper
level, basic symmetries in the underlying random process characterizing the distribution of
completion shapes. Second, we described a new method for numerical solution of the Fokker-Planck
equation in a shiftable-twistable basis. Finally, we used this solution to compute
stochastic completion elds, and demonstrated, both theoretically and experimentally, the
invariance of our computation under translations and rotations of the input pattern.

Appendix


Proof of Proposition 2.1
First observe that, by equation (2.5),
~ ';n
c ~ ';n ~ ';n (T ~x 0
and that the composition rule (2.8) for shift-twist transformations implies that
The result now follows by applying the dening equation (2.3) for shiftability-twistability to
express the right hand side of equation (11.2) in terms of the basis functions, ~ k;m , and then
substituting the resulting equation into equation (11.1).
Analysis of the accuracy and shift-twist invariance of the advection process
In the continuum, the advection process is shift-twist invariant. Consequently, the greater
the accuracy with which the continuous advection process (for time, t) is modeled by
the discrete advection transformation, A, the greater the degree to which A is shift-twist
invariant.
If the coecient vector, c ~ k;! , evolves according to equation (5.7), then the the p.d.f.,
advect in the continuum according to equation (5.4) provided that the basis, ~ k;! ,
is perfectly shiftable-twistable and includes all (essentially) non-zero -frequencies of the
interpolation of the continuous solution, P
of the advection process. 19
A numerical study showed that in the case that the initial p.d.f is a three-dimensional
Gaussian,
, the number, N , of essentially non-zero -Fourier
Since the basis functions are only eectively shiftable-twistable, there is inevitably some error in the
advection process. This error is not Euclidean invariant.
19 For small t  0:1, the number of nonzero -frequencies of the interpolation functions is much less than
that of P .
series coecients of P (~x; ; t) is given by the formula
at
~x
for some constants a and b. In particular if the spatial resolution of the initial conditions is
increased by a factor of two in all three variables, then the number of frequencies, !, must
be multiplied by four to prevent aliasing in the advection process, and the number of spatial
basis functions, ~ k (~x), must also be multiplied by four to ensure that the basis is shiftable.
Proof of Proposition 5.2
The standard explicit nite dierence scheme for the -diusion of P (~x;  ; t) on the grid of
points,
Substituting equation
~ k;! c ~ k;! (t) ~ k (~x)e i! and equating coecients
of ~ k (~x) yields the formula
c ~ k;! (t +t)e
where we have used the fact that  Equation (5.9) now follows by equating
coecients of e i!n .
The following argument shows that the diusion transformation, D, is shift-twist invari-
ant. First, a similar argument to the one just presented shows that equation (5.9) implies
that
for all  2 S 1 . Equation (11.6) implies that D of the shift-twist by (~x of P is given by
(R  0
(R  0
(R  0
which is equal to the shift-twist by (~x of D of P , i.e., D is shift-twist invariant.
Proof of Equation (6.3)
Equation (6.3) follows immediately from the two equations
2e
e
and the shiftability formula
To verify equation (11.8), observe that
and that, by completing the square, k~x ~ kk 2 +k~x ~ 'k

Acknowledgments

J.W.Z. was supported (in part) by the Albuquerque High Performance Computing Center.



--R

Putative Strategies of Scene Segmentation in Monkey Visual Cortex
Minimization of Grid Orientation E
Neurodynamics and Brain Mechanisms
Ten Lectures on Wavelets

Uncertainty Relation for Resolution in Space
Complete Discrete 2-D Gabor Transforms by Neural Networks for Image Analysis and Compression
Turning a Corner in Vision Research
The Design and Use of Steerable Filters

Neural Dynamics of Form Perception: Boundary Completion
A Computational Model of Neural Contour Pro- cessing

The curve of least energy
Toward Discrete Geometric Models for Early Vision

Image Representation Using 2D Gabor Wavelets
A Neural Model of Contour Integration in Primary Visual Cortex

A Computational Investigation into the Human Representation and Processing of Visual Information
Elastica and Computer Vision
Curvature Consistency and Curve De- tection
The Generalized Gabor Scheme of Image Representation in Biological and Machine Vision
The Detection of Globally Salient Structures Using a Locally Connected Network
Steerable Wedge Filters for Local Orientation Analysis
Shiftable Multiscale Trans- forms
Analytic Solution of

Foundations of Vision
A Neural Model of Illusory Contour Shape and Salience
Local Parallel Computation of
A Comparison of Measures for Detecting Natural Shapes in Cluttered Backgrounds
Computing
Salient Contour Extraction by Temporal Binding in a Cortically- Based Network
--TR
visual cortex;fokker-planck equation;boundary completion;euclidean invariant computation;shiftable basis