Introduction

Generalized sampling [1, 2] is a framework for estimating representations of functions in different bases and frames in a numerically stable manner. This can for example be relevant in processing MRI images, where hardware often enforces initial frequency measurements, but where wavelets may be better suited for representing the final image. Simply using a discrete inverse Fourier transform often results in the well-known Gibbs phenomenon, with unwanted oscillating behavior in the reconstructed image. This paper documents a toolbox for performing generalized sampling in Julia [3].

The theory of generalized sampling does not restrict the type of bases to consider, but the applications have focused on Fourier bases and multiscale representations like wavelets.

Mathematically, samples of a function f in a Hilbert space H with respect to a sample basis {s_n}_n∈N consists of inner products {〈f, s_n〉}_n∈N. In generalized sampling we want to use these samples to estimate the inner products {〈f, r_n〉}_n∈N, where {r_n}_n∈N is another basis for ℋ. The basis {r_n}_n∈N is used for reconstructing f. In practice we only consider a finite number of sampling and reconstruction functions, i.e., we have access to the samples $w_m^s=\langle f,s_m\rangle$M1 \documentclass[10pt]{article} \usepackage{wasysym} \usepackage[substack]{amsmath} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage[mathscr]{eucal} \usepackage{mathrsfs} \usepackage{pmc} \usepackage[Euler]{upgreek} \pagestyle{empty} \oddsidemargin -1.0in \begin{document} \[ w_m^s = \left\langle {f,{s_m}} \right\rangle \] \end{document} , 1 ≤ m ≤ M . From these samples we estimate the coefficients in the reconstruction basis, ${\tilde{w}}_n^r\approx w_n^r=\langle f,r_n\rangle$M2 \documentclass[10pt]{article} \usepackage{wasysym} \usepackage[substack]{amsmath} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage[mathscr]{eucal} \usepackage{mathrsfs} \usepackage{pmc} \usepackage[Euler]{upgreek} \pagestyle{empty} \oddsidemargin -1.0in \begin{document} \[ \widetilde w_n^r \approx w_n^r = \left\langle {f,{r_n}} \right\rangle \] \end{document} , 1 ≤ n ≤ N, which are used to compute an approximation of f,

An especific example that will be used later is where 〈f, s_m〉 represent frequency samples of f and the reconstruction basis ${\left\{r_n\right\}}_{n=1}^N$M4 \documentclass[10pt]{article} \usepackage{wasysym} \usepackage[substack]{amsmath} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage[mathscr]{eucal} \usepackage{mathrsfs} \usepackage{pmc} \usepackage[Euler]{upgreek} \pagestyle{empty} \oddsidemargin -1.0in \begin{document} \[ \left\{ {{r_n}} \right\}_{n = 1}^N \] \end{document} is given by translates of a Daubechies scaling function φ yielding a system on the form

that has been adapted to the interval $\left[-\frac{1}{2},\frac{1}{2}\right]$M6 \documentclass[10pt]{article} \usepackage{wasysym} \usepackage[substack]{amsmath} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage[mathscr]{eucal} \usepackage{mathrsfs} \usepackage{pmc} \usepackage[Euler]{upgreek} \pagestyle{empty} \oddsidemargin -1.0in \begin{document} \[ \left[ { - {\textstyle{1 \over 2}},{\textstyle{1 \over 2}}} \right] \] \end{document} following [4]). The simplest Daubechies scaling function is the Haar scaling function given by

Other Daubechies scaling functions are at least continuous, but cannot be written in closed form in terms of elementary functions.

The actual computation of the reconstruction coefficients ${\tilde{w}}^r={\left\{{\tilde{w}}_n^r\right\}}_{n=1}^N$M8 \documentclass[10pt]{article} \usepackage{wasysym} \usepackage[substack]{amsmath} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage[mathscr]{eucal} \usepackage{mathrsfs} \usepackage{pmc} \usepackage[Euler]{upgreek} \pagestyle{empty} \oddsidemargin -1.0in \begin{document} \[ {\tilde w^r} = \{ \tilde w_n^r\} _{n = 1}^N \] \end{document} is performed by solving a least squares problem. The infinite change of basis matrix between the sampling and reconstruction subspaces has (i, j)’th entry (r_j, s_i). We consider a finite M × N section of this matrix, denoted by T. The reconstruction coefficients are computed as the least squares solution

For large M and N , e.g. of the size used in 2D examples with images, T is not accessible as a stored matrix. Instead we need to work with T as an operator, i.e., store a representation that allow us to compute products with T and T* and thereby solving (1) with e.g. a conjugate gradient algorithm [6, p. 637].

Examples

We have two goals with the GeneralizedSampling package: It should be fast and easy to use. We have therefore put effort into providing only a few necessary high-level functions and hiding the lower level details. In essence, we need to compute a representation of a change of basis matrix T and solve least squares problems like (1) with this representation.

An example included in the package is the reconstruction of a truncated cosine (with inspiration from [5]). The results are seen in Fig. 1 using a different number of scaling functions to illustrate how this affects the reconstruction. In this example the higher order, continuous Daubechies scaling functions are not well suited to represent the discontinuity.

A reconstruction like this can be computed in only a few lines of code once the frequency measurements are available. For the truncated cosine we have a closed-form expression for the Fourier transform:

In Julia this can be implemented in the following way with a ternary operator (the line break is only for typographical purposes):

julia > ftcos (xi) = ( abs(xi) == 1 ? 0.25 :
im*xi *(1 + exp (-pi*xi*im ))/(2* pi *(1 - xi ^2)) )
julia > @vectorize_1arg Float64 ftcos

With the Fourier transform defined the reconstruction can be computed as follows:

julia > using GeneralizedSampling
julia > xi = grid (128 , 0.5)
julia > fhat = ftcos (xi)
julia > T = Freq2Wave (xi , " haar ", 6)
julia > wcoef = T \ fhat
julia > using IntervalWavelets
julia > x, y = weval ( real ( wcoef ), " haar ", 10)

To compute the reconstruction with the Daubechies 4 scaling functions, one simply has to replace “haar” with “db4”.

As an example of reconstruction of 2D functions/images, consider reconstruction of a simulated brain made with the Matlab [8] software released along with [7] (available at http://bigwww.epfl.ch/algorithms/mriphantom). Here we do not have access to a “ground truth”, but only the frequency measurements. From 512² frequency measurements we reconstruct 256² Daubechies 4 scaling functions and the result is seen in Fig. 2.

Runtime and technical comparison

Julia is a dynamic language with a fast JIT compiler. A consequence is that a function is compiled the first time it is called, causing an overhead in terms of time and memory. All subsequent calls are, however, without this compiling overhead. The runtimes reported in this section are not for the first run.

The examples were carried out on a normal laptop (2.60 GHz Intel Core i7, 8 GB RAM) running GNU/Linux and are summarized in Table 1. The background for the experiments are as follows: We compare the GeneralizedSampling package with the Matlab code released with [5] (available at Error! Hyperlink reference not valid./research/afha/code). To the best of our knowledge this is the only other publicly available software for generalized sampling. In this connection two comparisons are relevant: Computing the representation of the change of matrix (initialization/“init”) and using this representation to compute the least squares solution (solution/“sol”). In both cases the initialization step is fast for Haar scaling functions: All Fourier transforms have simple, closed-form expressions that are easily vectorized. For higher order Daubechies scaling functions all computations rely on iterative procedures.

In both packages the solution step is based on a conjugate gradient-like algorithm, where the computational cost is dominated by multiplication with the change of basis matrix and its adjoint. In Matlab the built-in lsqr function is used and in Julia a custom implementation of the conjugate gradient algorithm is used.

In GeneralizedSampling we have relied on Julia’s ability to write functions that modify their arguments in-place to drastically reduce the memory consumption in an iterative algorithm like conjugate gradients. Julia’s @time macro and the benchmarking mentioned below makes it easy to estimate the memory allocation of a function. In Matlab we have used an undocumented feature for the profiler, namely profile -memory on to measure the allocated memory.

Especially for fast runtimes it is not accurate to rely on timing a single run of a function (using e.g. tic and toc in Matlab). In Matlab the times are obtained with the built-in timeit function and in Julia we use the benchmark package BenchmarkTools available at https://github.com/JuliaCI/BenchmarkTools.jl. For small problems the speed of the Julia and Matlab code are comparable, but for large problems the Julia code is significantly faster. The one exception is for the “Uniform 2D” example: The Julia package is using a more general algorithm for non-equidistant fast Fourier transforms, whereas the Matlab package is considering a special case where the standard fast Fourier transform is applicable.

The total memory allocation of the solver is much smaller in the Julia code than in the Matlab code. In most cases, the memory allocation of the initialization code is smaller the Julia code than in the Matlab code, although this is probably not as important since it is only performed once.

Note that the “sol” times and number of iterations are not directly comparable, since the stopping criteria for the least squares solver may be different and [5] use L²([0, 1]^d) as the reconstruction space, whereas we use $L^2({[-\frac{1}{2},\frac{1}{2}]}^d)$M11 \documentclass[10pt]{article} \usepackage{wasysym} \usepackage[substack]{amsmath} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage[mathscr]{eucal} \usepackage{mathrsfs} \usepackage{pmc} \usepackage[Euler]{upgreek} \pagestyle{empty} \oddsidemargin -1.0in \begin{document} \[ {L^2}({[ - {\textstyle{1 \over 2}},{\textstyle{1 \over 2}}]^d}) \] \end{document} . But the time per iteration (“t/n”) is comparable.

Quality control

An extensive automated test suite has been developed along with the main code and is tested with Travis CI and AppVeyor, allowing tests to be performed on both Unix-like platforms and Windows. All public functions have doc-strings and a detailed documentation is included with the package in reStructuredText that can be exported to a multitude of formats. An HTML version of the documentation is hosted on Read the Docs. All details and links are available on the Generalized-Sampling package’s GitHub page.

Furthermore, all examples have been tested by all authors for bugs and readability.

Operating system

Any system capable of running Julia version 0.4 or higher.

Programming language

The GeneralizedSampling package has been tested with Julia version 0.4 and 0.5.

Additional system requirements

None.

Dependencies

The GeneralizedSampling package depends on the following Julia packages:

• NFFT (version 0.1.4)
• Wavelets (version 0.5.1)
• ArrayViews (version 0.6.4)
• VoronoiCells (version 0.1.3)
• IntervalWavelets (version 0.0.4)

List of contributors

• Jacobsen, Robert Dahl (Aalborg University) – Development.
• Nielsen, Morten (Aalborg University) – Testing and theoretical discussions.
• Rasmussen, Morten Grud (Aalborg University) – Testing and theoretical discussions.

Software location

Code repository

Name: GitHub

Persistent identifier: https://github.com/robertdj/GeneralizedSampling.jl

Licence: MIT “Expat”

Date published: 22/12/2016

Language

English.