Image Reconstruction from Undersampled Fourier Data Using the Polynomial Annihilation Transform

18 Pages·2015·2.56 MB·English

Noname manuscript No.
(will be inserted by the editor)
Image Reconstruction from Undersampled Fourier Data Using the Polynomial
Annihilation Transform
Rick Archibald · Anne Gelb · Rodrigo B. Platte
thedateofreceiptandacceptanceshouldbeinsertedlater
Abstract Fourier samples are collected in a variety of applications including magnetic resonance imaging
(MRI) and synthetic aperture radar (SAR). The data are typically under-sampled and noisy. In recent years,
l1 regularizationhasreceivedconsiderableattentionindesigningimagereconstructionalgorithmsfromunder-
sampled and noisy Fourier data. The underlying image is assumed to have some sparsity features, that is,
some measurable features of the image have sparse representation. The reconstruction algorithm is typically
designed to solve a convex optimization problem, which consists of a fidelity term penalized by one or more l1
regularization terms. The Split Bregman Algorithm provides a fast explicit solution for the case when TV is
used for the l1 regularization terms. Due to its numerical efficiency, it has been widely adopted for a variety
of applications.
AwellknowndrawbackinusingTV asanl1regularizationtermisthatthereconstructedimagewilltendto
defaulttoapiecewiseconstantimage.Thisissuehasbeenaddressedinseveralways.Recently,thepolynomial
annihilationedgedetectionmethodwasusedtogenerateahigherordersparsifyingtransform,andwascoined
the “polynomial annihilation (PA) transform.” This paper adapts the Split Bregman Algorithm for the case
when the PA transform is used as the l1 regularization term. In so doing, we achieve a more accurate image
reconstruction method from under-sampled and noisy Fourier data. Our new method compares favorably to
the TV Split Bregman Algorithm, as well as to the popular TGV combined with shearlet approach.
Keywords Fourier Data · l1 regularization · Split Bregman · Edge Detection · Polynomial Annihilation
1 Introduction
Data are acquired as partial Fourier samples in several applications, including magnetic resonance imaging
(MRI) and synthetic aperture radar (SAR). In an idealized situation, recovering images from partial Fourier
data may be done simply and efficiently by using the inverse fast Fourier transform (FFT). In practice, the
dataacquisitionsystemisusuallyunder-prescribedandnoisy.Moreover,theFourierdomainisnotwellsuited
for recovering the underlying image, which is generally only piecewise smooth.
In recent years l1 regularization has received considerable attention in designing image reconstruction
algorithms from under-sampled and noisy data for images that have some sparsity properties, that is, some
measurable features of the image have sparse representation.1 Also, l1 regularization provides a formulation
thatiscompatiblewithcompressedsensing(CS)applications,specifically,whenanimagecanbereconstructed
ComputerScienceandMathematicsDivision,OakRidgeNationalLaboratory,OakRidge,TN37831([email protected])
SchoolofMathematicalandStatisticalSciences,ArizonaStateUniversity,Tempe,AZ,85287([email protected])
SchoolofMathematicalandStatisticalSciences,ArizonaStateUniversity,Tempe,AZ,85287([email protected])
Address(es)ofauthor(s)shouldbegiven
1 Technically,thel0 normofanexpressionisabettermeasureofsparsity.However,thel0 normdoesnotmeettheconvexity
requirementsandisveryslowtocompute.Additionaldetailonusingthel1 norminplaceofthel0 inordertomeasuresparsity
canbefoundin[7].
2 RickArchibaldetal.
fromaverysmallnumberofmeasurements[4–6].Inparticular,thegoalforreconstructinganimagefromSAR
or MRI data is to solve
minJ(f) such that ||Ff −ˆf|| =0, (1)
2
f
where ˆf consists of samples of the Fourier transform of the unknown image, f. F contains a subset of rows of
a Fourier matrix, and J is an appropriate l1 regularization term, [10,20–23]. Typically for measured data the
related (TV) “denoising” problem,
minJ(f) such that ||Ff −ˆf|| <σ, (2)
2
f
is solved. It is in general still difficult to develop efficient and robust techniques for solving (2). The Split
Bregman Algorithm, [12], is a numerically efficient and stable algorithm that has successfully solved (2) for
a variety of applications. In this paper we use the Split Bregman Algorithm as a launching point to develop
a new technique for solving (2) based on the polynomial annihilation l1 regularization introduced in [27]. We
will demonstrate that our method yields improved accuracy in regions away from discontinuities, especially
in the case of under-sampled data. We adopt the standardizations and terminology from [12] to describe our
algorithm.
A well known drawback in using TV as an l1 regularization term is that the reconstructed image defaults
to a piecewise constant approximation. While suitable for some applications, in others it is desirable to see
more details. This has been addressed in several ways. For example, total generalized variation (TGV), which
generates a piecewise (typically quadratic) polynomial approximation in smooth regions, was developed in
[2]. Multi-wavelets have also been used to formulate sparsifying transforms, [24]. The polynomial annihilation
(PA) transform, which exploits the sparsity of the underlying image in the jump discontinuity domain, was
introduced in [27]. It was demonstrated there that generating a sparsifying transform based on the sparsity of
edges in the underlying image (as expressed in (11)), yields improved accuracy and convergence properties for
both image reconstruction and edge identification. In particular, high order accuracy is possible in regions of
smoothness, which has two important consequences, also true in multi-dimensions. First, it is possible to see
more variation in the underlying image, and second, fewer data points are needed to reconstruct an image.
In [27], the MATLAB CVX package [15,14] was used to implement (2) for the PA transform l1 regular-
ization. Although suitable for one-dimensional problems, CVX is not efficient enough for higher dimensional
problems. Because of this limitation, the technique introduced in [27] was not pursued for other applications,
suchascompressionorreducingthedimensionalityofthedataforefficientprocessing.Thispaperthusseeksto
expand the recent results in [27] in two ways. First, we will improve the efficiency of the numerical algorithms
in multi-dimensions. In this regard, we will adapt the Split Bregman Algorithm, [12], to the PA transform.
Once this is accomplished, we will demonstrate how the PA transform is an effective tool for reconstructing
piecewise smooth images from under-sampled Fourier data. The rest of the paper is organized as follows. In
Section 2 we describe the reconstruction problem in one dimension and discuss how the polynomial annihila-
tionedgedetectionmethodisusedtoconstructanl1 regularizationtermforsolving(2).InSection3wereview
the Split Bregman Algorithm for sparse Fourier data using the TV operator and demonstrate how it can be
adapted for the PA transform operator without increasing computational cost. In Section 4 we provide some
numerical examples and demonstrate that our method is robust to noise and undersampling. We compare our
results to those using TV and TGV (combined with shearlet regularization), based on the code given in [17].2
Concluding remarks are provided in Section 5.
2 Preliminaries
We begin by describing the one-dimensional problem. Let f ∈ R be a piecewise smooth function on [−1,1],
with f(−1)=f(1). Suppose we are given its first 2N +1 (normalized) Fourier coefficients,
1(cid:90) 1
fˆ(k)= f(x)e−ikπxdx, fork =−N,...,N. (3)
2
−1
2 In [27] our method compared favorably to multi-wavelet constructed regularization terms, [24]. We do not repeat those
experimentshere.
ImageReconstructionfromUndersampledFourierDataUsingthePolynomialAnnihilationTransform 3
We wish to recover f from (3) at a finite set of grid points, x = −1+j∆x, j = 0,··· ,2N, with ∆x = 1.3
j N
The standard Fourier partial sum reconstruction used to approximate periodic replicas of f is defined as
N
S f(x)= (cid:88) fˆeikπx. (4)
N k
k=−N
At the gridpoints x=x , we can write (4) as the linear system
j
Ff =ˆf, (5)
where F denotes the discrete Fourier transform matrix,
F =e−ikπxj, 0≤j ≤2N, −N ≤k ≤N. (6)
k,j
Hereˆf is the 2N +1 vector of Fourier coefficients given in (3), and f = {f }2N is the calculated solution to
j j=0
f(x ),j =0,...,2N.TheFFTcanbeusedtoefficientlysolve(5),[13].Sincetheunderlyingfunctionf isonly
j
piecewise smooth, the approximation given by (4) will yield the Gibbs phenomenon. Filtering is often used to
alleviate the Gibbs phenomenon and to reduce the effects of high frequency noise. The filtered approximation
is given by
N
Sσf(x)= (cid:88) σ fˆeikπx, (7)
N k k
k=−N
where σ is an admissible filter, [13]. The corresponding linear system at the gridpoints x=x is then
k j
Fσf =ˆf, (8)
where Fσ denotes the filtered discrete Fourier transform matrix
Fσ =σ e−ikπxj, 0≤j ≤2N, −N ≤k ≤N, (9)
k,j k
andonceagaintheFFTcanbeemployed.AswillbeseeninSection4,filteringmaycausetoomuchsmoothing
over the jump discontinuities, especially if the Gibbs oscillations are to be completely removed. It is also
evident that filtering can not address the poor reconstruction quality due to undersampling. Thus we seek
other methods of regularization. Since noise is inherent to all sampling systems, we will demonstrate that our
technique is effective for noisy input data,˜f =ˆf +(cid:15).
2.1 Sparse Representation in the Jump Function Domain
As mentioned in the introduction, we seek to regularize (2) by enforcing the sparsity of edges in the image
domain. To do this, we must first define the jump function of a piecewise smooth function. Note that the
words “jump” and “edge” are used interchangeably throughout the remainder of this paper.
Definition 1 Let f : [−1,1] :→ R. For all x ∈ (−1,1), let f(x−) and f(x+) denote its left and right hand
limits. The jump function of f is defined at each x as
[f](x)=f(x+)−f(x−). (10)
3 For simplicity we choose 2N +1 equally spaced grid points to match the number of Fourier coefficients. The techniques
describedinthispaperareeasilyextendedtodifferentgriddingschemesintheimagedomain.
4 RickArchibaldetal.
From the above definition, we see that [f] is zero everywhere that f is continuous, while it takes on the jump
value at each jump discontinuity location. We make the assumption that there is at most one jump within a
cell I =[x ,x ). Thus, if [f](x ) is the value of the jump that occurs within the cell I , we can write
j j j+1 j j
2N−1
(cid:88)
[f](x)= [f](x )χ (x), (11)
j j
j=0
where χ (x) is defined as
j
(cid:26)
1 if x∈I
χ (x)= j
j 0for all other x.
For simplicity, the numerical algorithms used in this investigation all place the jump discontinuity at the left
boundary of its corresponding cell. (This also means that the solution to (2) can be written as the expansion
coefficients of the standard basis.)
Since [f](x) = 0 for almost all grid point values x = x , it is apparent that (11) has only a few nonzero
j
coefficients. Therefore we say that [f](x) has sparse representation, or equivalently, that the jump function
domain of f is sparse. Hence we seek to regularize (5) in the form of (2) using (11). This will require an
approximation to (11) that is suitable for convex optimization problems.
2.2 Convex Optimization Using Sparsity of Edges
Whilethereareavarietyofwaystoapproximate(11),forourpurposewewillusethepolynomial annihilation
edge detection method,[1].Theadvantageinusingthepolynomialannihilationmethodisthatitishighorder,
meaning that in regions of smoothness the coefficients of the approximation of (11) will indeed be sparse.
Moreover, it is simple to generate a transform matrix for l1 regularization. This was accomplished in [27],
where the polynomial annihilation (PA) transform matrix was introduced.
The polynomial annihilation edge detection method, [1], is defined as
1 (cid:88)
Lmf(x)= c (x)f(x ), (12)
qm(x) j j
xj∈Sx
whereS isthelocalsetofm+1gridpointsfromthesetofgivengridpointsaboutx,c (x)arethepolynomial
x j
annihilation edge detection coefficients, (13), and qm(x) is the normalization factor, (14). Each parameter of
the method can be further described as:
– S : For any particular cell I =[x ,x ), there are m possible stencils, S of size m+1, that contain the
x j j j+1 x
interval I . For simplicity, we assume that the stencils are centered around the interval of interest, I , and
j j
are given by
SIj ={xj−m2 ,··· ,xj+m2 }, SIj ={xj−m2+1,··· ,xj+m2−1}
for m even and odd respectively. For non-periodic solutions the stencils are adapted to be more one sided
as the boundaries of the interval are approached, [1]. To avoid cumbersome notation, we write S as the
x
generic stencil unless further clarification is needed.
– c (x): The polynomial annihilation edge detection coefficients, c (x), j =1,...,m+1, are constructed to
j j
annihilate polynomials up to degree m. They are obtained by solving the system
(cid:88) c (x)p (x )=p(m)(x), j =1,...,m+1, (13)
j (cid:96) j (cid:96)
xj∈Sx
where p , (cid:96)=0,...,m, is a basis of for the space of polynomials of degree ≤m.
(cid:96)
– qm(x): The normalization factor, qm(x), normalizes the approximation to assure the proper convergence
of Lmf to the jump value at each discontinuity. It is computed as
(cid:88)
qm(x)= c (x), (14)
j
xj∈Sx+
where S+ is the set of points x ∈S such that x ≥x.
x j x j
ImageReconstructionfromUndersampledFourierDataUsingthePolynomialAnnihilationTransform 5
It was shown in [1] that the the polynomial annihilation edge detection method has mth order convergence
away from the jump discontinuities. More precisely, this accuracy is seen in the region outside the stencil that
contains the jump. Oscillations will naturally develop in the region of the jump, which will increase with m.
The consequences of these oscillations on our method will be discussed more in Section 4.
Ifthesolutionvectorf isonuniformpointsin[−1,1],thenwesolve(12)onthesetofpointsx =−1+j∆x,
j
j = 0,...,2N, with ∆x = 1. In this case there is an explicit formula for the polynomial annihilation edge
N
detection coefficients, independent of location x, computed as ([1])
m!
c = , j =1,...,m+1. (15)
j (cid:81)m+1 (j−k)∆x
k=1,k(cid:54)=j
We can now define the polynomial annihilation (PA) transform matrix, Lm, as
c(j,l)
Lm = , 0≤l≤2N, 0≤j <2N, (16)
j,l qm(x )
l
where
(cid:26)c 0<j−l−(cid:98)m(cid:99)+s(j,l)≤m+1
c(j,l)= j−l−(cid:98)m2(cid:99) 2
0 otherwise
and
l−(cid:98)m(cid:99) l≤(cid:98)m(cid:99)
 2 2
s(j,l)= l+m−(cid:98)m(cid:99)−2N l+m−(cid:98)m(cid:99)>2N
2 2
0 otherwise.
For example, not assuming periodicity, the banded matrix Lm for m=4 can be written as
 
3 −12 18 −12 3
−1 4 −6 4 −1 
 
 1 −4 6 −4 1 
 
 1 −4 6 −4 1 
1 . 
Lm =  . .
3 . 
 
 1 −4 6 −4 1 
 
 1 −4 6 −4 1 
 
 1 −4 6 −4 1 
−3 12 −18 12 −3
The PA transform produces a vector with small, nonzero values in the smooth regions of f and large
values at the jump locations. Therefore, by minimizing ||Lmf|| , we encourage a solution f that has sparse
1
representation in the jump function domain, as given in (11). As was demonstrated in [27], using the PA
transformin(2)reducestheGibbsoscillationswithoutsmearingoverthediscontinuitiesorcausingastaircasing
effect.
Remark 1 When m=1, applying the PA transform is equivalent to TV regularization (up to a normalization
constant).AsisevidentinFigure2,usingTV asaminimizingconstraintcausesa“staircasing”effect,sinceit
encourages a solution that minimizes the differences between the solution vector components. In this regard,
the PA transform can be seen as high order, meaning that the accuracy away from the jump discontinuities is
O(∆xm).
As mentioned in the introduction, in [27], the MATLAB CVX package was used to implement (2) with
(16).Whilesuitableforonedimensionalproblems,itisnotefficientformultipledimensions.Belowwedescribe
how the split Bregman algorithm, [12], can be adapted to incorporate the PA transform, (16), into (2), and
thus improve the overall efficiency.
6 RickArchibaldetal.
3 The Split Bregman Algorithm for l1 regularization of Sparse Fourier Data
Bregmaniteration,firstdevelopedtofindextremaofconvexfunctionals,[3],hasbeenusedinmanyapplications
wheretheoptimizationproblemisoftheform(1)or(2),[18,22,28].TheSplitBregmanAlgorithm,developed
in[12],wasshowntobeequivalenttoBregmaniteration.Itspopularityisowedtothefactthatitisveryfast,
typicallyhavingthecomputationalcostontheorderofanFFT,andalsothatthenonlinearstepsinvolveonly
soft thresholding, while all other aspects of the algorithm involve solving invertible linear systems.
In recent years the Split Bregman Algorithm has been used for solving a broad class of l1 regularized
optimization problems. In particular, the Split Bregman Algorithm was used to solve (2) for the two dimen-
sional Fourier transform and termed the “sparse MRI data reconstruction problem” in [18,20,25].4 In this
case we assume that f : R2 → R is a periodic piecewise smooth function on [−1,1]2 for which we are given
the (normalized) Fourier coefficients,
1(cid:90) 1 (cid:90) 1
fˆ(k,l)= f(x,y)e−iπ(kx+ly)dydx, fork,l=−N,...,N. (17)
4
−1 −1
We wish to recover f from (17) at a finite set of uniform grid points, (x ,y ), for i,j = 0,··· ,2N. The
i j
fidelitytermin(2)isdirectlyextendedwithF beingtheanalogoustwo-dimensionalFouriertransformmatrix
of (6). Although the polynomial annihilation edge detection is an inherently multi-dimensional method, [1],
our results indicate that using a dimension by dimension construction, J and J for the two dimensional
x y
regularization term J, is more efficient for uniformly spaced data.
We now seek f ={f(x ,y ):0≤i,j ≤2N} that solves the convex optimization problem
i j
min(J f +J f) such that ||Ff −ˆf|| ≤σ (18)
x y 2
f
In the case where only partial information is given, the corresponding convex optimization problem we solve
is
min(J f +J f) such that ||MFf −ˆf|| ≤σ. (19)
x y 2
f
Here the matrix M represents a ‘row selector’ matrix, which comprises a subset of the rows of an identity
matrix, corresponding to the known set of Fourier samples. We note that the one dimensional version posed
in (2) can also be similarly adapted. Fast algorithms were developed for (19) as a means for approximating
(18) in [28] and [12]. It was first demonstrated in [28] that (19) could be solved using a Bregman iteration of
the sequence of two unconstrained problems of the form
µ
fk+1 =min||J f|| +||J f|| + ||MFf −ˆfk|| , (20a)
f x 1 y 1 2 2
ˆfk+1 =ˆfk+ˆf −MFfk+1, (20b)
where µ > 0 is an optimization parameter. The computational challenge comes in solving (20a). Using the
SplitBregmanmethodology,afastexplicitsolutionwasdevelopedin[12]forthecasewhereTV isusedforthe
l1 regularization terms, that is, when J =∇ and J =∇ . We outline the basic steps of the Split Bregman
x x y y
Algorithm for the TV denoising model below and then demonstrate in Section 3.2 how it can be extended to
the case where J and J are the polynomial annihilation transforms, resulting in a more accurate solution
x y
without additional computational cost. Since we do not explicitly motivate the algorithm, we refer interested
readers to [12] for more general details.
4 AlthoughdesignedforMRIdata,itisapplicablewheneverdataaresampledintheFourierdomain,especiallywhendimension
reductionisdesirable.
ImageReconstructionfromUndersampledFourierDataUsingthePolynomialAnnihilationTransform 7
3.1 The Split Bregman Algorithm for TV Denoising
Following the framework in [12], for J = ∇ and J = ∇ , we make the replacements d ← ∇ f and
x x y y x x
d ←∇ f. Using
y y
(cid:88)(cid:113)
||∇v|| = |∇ v |2+|∇ v |2,
1 x i,j y i,j
i,j
we obtain as the augmented optimization of (20a):
(cid:88)(cid:113) λ
min |d |2+|d |2 + ||d −∇ f −b ||
f,dx,dy i,j x,i,j y,i,j 2 x x x 2
λ µ
+ ||d −∇ f −b || + ||MFf −ˆf|| . (21)
2 y y y 2 2 2
The variables b and b arise from derivation of the Split Bregman Algorithm, [12], and their calculation is
x y
given in Algorithm 1. Solving (21) is accomplished in steps. When f is held fixed, the exact optimization of
d and d can be calculated using the shrink operator, [26], as
x y
∇ f +b ∇ f +b
dopt =max(s−1/λ,0) x x and dopt =max(s−1/λ,0) y y, (22)
x s y s
where
(cid:113)
s= |∇ f +b |2+|∇ f +b |2. (23)
x x y y
When d and d are held fixed, optimizing (21) reduces to the l2 minimization problem
x y
λ λ µ
min ||d −∇ f −b || + ||d −∇ f −b || + ||MFf −ˆf|| . (24)
f 2 x x x 2 2 y y y 2 2 2
Since the subproblem in (24) is differentiable, the optimal solution f can be found by differentiating with
respect to f and setting the result equal to zero, arriving at
(µFTMTMF +λ∇T∇ +λ∇T∇ )f =rhs, (25)
x x y y
where
rhs=µFTMTˆf +λ∇T(d −b )+λ∇T(d −b ). (26)
x x x y y y
Using the identities ∇T∇=−(cid:52) and FT =F−1 produces the inverse problem
F−1KFf =rhs, (27)
whereK =(µMTM−F(cid:52)F−1).NotethatK isadiagonaloperator.Therefore,theoptimalsolutionf to(24)
can be calculated at the cost of two Fourier transforms as
fopt =F−1K−1F(rhs). (28)
Incorporatingtheabovecalculationsinto(19)leadstoafastalgorithmfortheTV denoisingproblem,[12]:
Algorithm 1 The Split Bregman Algorithm for TV Denoising
Define J =∇ and J =∇ in (19).
x x y y
Initialize k =0, f0 =F−1MTˆf, and b0 =b0 =d0 =d0 =0
x y x y
while ||MFfk−ˆf|| >σ
2
fk+1 =F−1K−1Frhsk
dk+1 =max(sk−1/λ,0)∇xfk+bkx
x sk
dk+1 =max(sk−1/λ,0)∇yfk+bky
y sk
bk+1 =bk+1+(∇ fk+1−dk+1)
x x x x
bk+1 =bk+1+(∇ fk+1−dk+1)
y y y y
8 RickArchibaldetal.
ˆfk+1 =ˆfk+ˆf −MFfk+1
k =k+1
end
The iteration terms sk and rhsk are provided in (23) and (26), where all variables are superscripted with k.
The first five steps of Algorithm 1 can alternatively be iterated in a short inner loop before updating ˆfk+1.
We note that Algorithm 1 works for the one dimensional optimization problem given in (2) by using the one
dimensional Fourier transform with J as the TV operator and eliminating the y components, b and d .
y y
3.2 The Split Bregman Algorithm for PA l1 Regularization
Asnotedintheintroduction,themainlimitationinusingtheTV operatorasanl1 regularizationtermisthat
it will generate a piecewise constant solution for the underlying signal or image. In [27] it was demonstrated
that the TV operator is equivalent to the PA transform operator, (12), with m=1, and that better accuracy
may be achieved when m > 1, especially if the underlying signal or image has significant variation between
discontinuities. We now demonstrate that solving (2) with the PA transform operator as the l1 regularization
term can be made efficient via an extension of the Split Bregman Algorithm for TV denoising.
We begin by defining J := Lm and J := Lm in (19), where Lm and Lm are the respective directional
x x y y x y
PA transform operators defined by in one dimension by (16). Bregman iteration of (19) is then given by the
sequence of two unconstrained optimization problems of the form
µ
fk+1 =min||Lmf|| +||Lmf|| + ||MFf −ˆfk|| , (29a)
f x 1 y 1 2 2
ˆfk+1 =ˆfk+ˆf −MFfk+1, (29b)
where µ > 0 is again an optimization parameter. Following the technique described in Section 3.1, we write
the augmented problem corresponding to (29a) and (29b) as
(cid:88)(cid:113) λ
min |d |2+|d |2 + ||d −Lmf −b ||
f,dx,dy i,j x,i,j y,i,j 2 x x x 2
λ µ
+ ||d −Lmf −b || + ||MFf −ˆf|| . (30)
2 y y x 2 2 2
When f is held fixed, the exact optimization of d and d can be calculated using the shrink operator, [26],
x y
similar to the derivation of (22), as
dopt =max(s−1/λ,0)Lmx f +bx and dopt =max(s−1/λ,0)Lmy f +by, (31)
x s y s
where s is given by
(cid:113)
s= |Lmf +b |2+|Lmf +b |2. (32)
x x y y
When d and d are held fixed, (30) reduces to the l2 minimization problem
x y
λ λ µ
min ||d −Lmf −b || + ||d −Lmf −b || + ||MFf −ˆf|| . (33)
f 2 x x x 2 2 y y y 2 2 2
Onceagain,sincethesubproblemin(33)isdifferentiable,theoptimalsolutionf canbefoundbydifferentiating
with respect to f and setting the result equal to zero, arriving at
(cid:16) (cid:17)
µFTMTMF +λ(Lm)TLm+λ(Lm)TLm f =rhs, (34)
x x y y
where
rhs=µFTMTˆf +λ(Lm)T(d −b )+λ(Lm)T(d −b ). (35)
x x x y y y
ImageReconstructionfromUndersampledFourierDataUsingthePolynomialAnnihilationTransform 9
Of significant importance here is that
(cid:16) ∂2m ∂2m (cid:17)
(Lm)TLm =(Lm)TLm+(Lm)TLm =−1m +
x x y y ∂x2m ∂y2m
is a diagonal operator in the Fourier domain. Hence the inverse problem is solved by
(cid:0)F−1K F(cid:1)f =rhs. (36)
LM
Note that K = (µMTM +F(Lm)TLmF−1) is a diagonal operator. Therefore, the optimal solution f to
LM
(24) can be calculated at the cost of two Fourier transforms as
fopt =F−1K−1 Frhs. (37)
LM
Thus, just as in the case for the TV denoising model, solving (19) with PA l1 regularization, (12), can be
made numerically efficient via the Split Bregman Algorithm:
Algorithm 2 The Split Bregman Algorithm for PA l1 Regularization
Define J =Lm and J =Lm in (19).
x x y y
Initialize k =0, f0 =F−1MTˆf, and b0 =b0 =d0 =d0 =0
x y x y
while ||MFfk−ˆf|| >σ
2
fk+1 =F−1K−1 Frhsk
LM
dk+1 =max(sk−1/λ,0)Lmxfk+bkx
x sk
dk+1 =max(sk−1/λ,0)Lmyfk+bky
y sk
bk+1 =bk+1+(Lmfk+1−dk+1)
x x x x
bk+1 =bk+1+(Lmfk+1−dk+1)
y y y y
ˆfk+1 =ˆfk+ˆf −MFfk+1
k =k+1
end
The iteration terms sk and rhsk are provided in (32) and (35), where all variables are superscripted with k.
The first five steps of Algorithm 2 can alternatively be iterated in a short inner loop before updating ˆfk+1.
As before, Algorithm 2 works for the one dimensional optimization problem (2) by using the one dimensional
FouriertransformwithJ =Lm andeliminatingthey componentsb andd .Finally,wenotethatAlgorithm
y y
2 is equivalent to Algorithm 1 when m=1.
4 Numerical Results
We are now ready to demonstrate the Split Bregman Algorithm for PA l1 regularization, given by Algorithm
2, for solving the one and two dimensional optimization problems. We will assume that we are given a finite
number of Fourier coefficients for the underlying signal or image. The wave numbers of the chosen coefficients
were then drawn from a normal distribution with standard deviation 2N/6 (rounded to the nearest integer),
where 2N +1 is the number of recovered function values.
For our one dimensional experiments, we consider the two test functions defined on [−1,1]:
Example 1
(cid:26)−1−xif −1≤x<0 cosπ2x if 1≤x<−12
f (x)= ; f (x)= cos3πx if −1 ≤x< 1
a 1−x otherwise b 2 2 2
cos7πx if 1 ≤x≤1
2 2
10 RickArchibaldetal.

Test Function Fourier Reconstruction Filtered Fourier Reconstruction
1 Uniform Grid 1 Test Function 1 Test Function
0.8 0.8 0.8
0.6 0.6 0.6
0.4 0.4 0.4
0.2 0.2 0.2
0 0 0
−0.2 −0.2 −0.2
−0.4 −0.4 −0.4
−0.6 −0.6 −0.6
−0.8 −0.8 −0.8
−1 −1 −1
−1 −0.5 0 0.5 1 −1 −0.5 0 0.5 1 −1 −0.5 0 0.5 1
(a) fa(x) (b) SNfa(x) (c) SNσfa(x)

1 1 1
0.8 0.8 0.8
0.6 0.6 0.6
0.4 0.4 0.4
0.2 0.2 0.2
0 0 0
−0.2 −0.2 −0.2
−0.4 −0.4 −0.4
−0.6 −0.6 −0.6
−0.8 −0.8 −0.8
−1 Test Function −1 Fourier Reconstruction −1 Filtered Fourier Reconstruction
Uniform Grid Test Function Test Function
−1 −0.5 0 0.5 1 −1 −0.5 0 0.5 1 −1 −0.5 0 0.5 1
(d) fb(x) (e) SNfb(x) (f) SNσfb(x)
Fig. 1 (a)fa and(d)fb plottedon2N+1=129uniformgridpoints;(b)and(e)Fourierpartialsumapproximation,(4),forfa
and f respectively; (c) and (f) Filtered Fourier partial sum approximation, (7), with filter σ =exp(−α(k/N)2p), α=32 and
b k
2p=4,forfa andfb respectively.

PA m=1 PA m=2 PA m=3
1 Test Function 1 Test Function 1 Test Function
0.8 0.8 0.8
0.6 0.6 0.6
0.4 0.4 0.4
0.2 0.2 0.2
0 0 0
−0.2 −0.2 −0.2
−0.4 −0.4 −0.4
−0.6 −0.6 −0.6
−0.8 −0.8 −0.8
−1 −1 −1
−1 −0.5 0 0.5 1 −1 −0.5 0 0.5 1 −1 −0.5 0 0.5 1
(a) PA(m=1)forfa (b) PA(m=2)forfa (c) PA(m=3)forfa

1 1 1
0.8 0.8 0.8
0.6 0.6 0.6
0.4 0.4 0.4
0.2 0.2 0.2
0 0 0
−0.2 −0.2 −0.2
−0.4 −0.4 −0.4
−0.6 −0.6 −0.6
−0.8 −0.8 −0.8
−1 PA m=1 −1 PA m=2 −1 PA m=3
Test Function Test Function Test Function
−1 −0.5 0 0.5 1 −1 −0.5 0 0.5 1 −1 −0.5 0 0.5 1
(d) PA(m=1)forf (e) PA(m=2)forf (f) PA(m=3)forf
b b b
Fig. 2 ResultsusingAlgorithm2with2N =128,µ=1(29a),λ=.015andm=1,2and3.
Figures 1(a) and (d) show f and f plotted on a uniform grid of 2N +1 = 129 points in [−1,1]. The
a b
corresponding Fourier and filtered Fourier partial sum approximations, (4) and (7), are illustrated in Figure
1(b) and (e), and in Figure 1(c) and (f) respectively. In this reconstruction we used 50% of the original
2N + 1 = 129 Fourier coefficients sampled. It is not surprising that having a sparse sampling of Fourier
data severely limits the reconstruction quality when directly applying either the standard or filtered Fourier
approximation. The convex optimization framework is clearly better suited in this case.
As illustrated in Figure 2, the approximation is readily improved by employing the PA transform, (16),
for the l1 regularization term in (2). Algorithm 2 was used in all cases with the number of iterations ranging
from 30-40.