Acoustic Wave Field Reconstruction from Compressed Measurements with Application in Photoacoustic Tomography
Marta M. Betcke, Ben T. Cox, Nam Huynh, Edward Z. Zhang, Paul C. Beard, Simon R. Arridge
JJOURNAL OF L A TEX CLASS FILES, VOL. 11, NO. 4, DECEMBER 2012 1
Acoustic Wave Field Reconstruction fromCompressed Measurements with Application inPhotoacoustic Tomography
Marta M. Betcke, Ben T. Cox, Nam Huynh, Edward Z. Zhang, Paul C. Beard and Simon R. Arridge
Abstract —We present a method for the recovery of compres-sively sensed acoustic fields using patterned, instead of point-by-point, detection. From a limited number of such compressedmeasurements, we propose to reconstruct the field on the sensorplane in each time step independently assuming its sparsityin a Curvelet frame. A modification of the Curvelet frame isproposed to account for the smoothing effects of data acquisitionand motivated by a frequency domain model for photoacoustictomography. An ADMM type algorithm, SALSA, is used torecover the pointwise data in each individual time step fromthe patterned measurements. For photoacoustic applications, thephotoacoustic image of the initial pressure is reconstructed usingtime reversal in k -Wave Toolbox. Index Terms —ADMM methods, compressed sensing, curveletframe, L1 minimization, photoacoustic tomography.
I. I
NTRODUCTION C OMPRESSED Sensing (CS) is a new measurementparadigm, which allows for the reconstruction of sparsesignals sampled at sub-Nyquist rates. Nowadays, it is commonunderstanding that many digital signals and images admit anadequate representation with far fewer coefficients than theiractual length. This phenomena is known as compressibilityand it has been a driving force in many image processingapplications, most notably the image compression algorithmsJPEG and its successor JPEG 2000. CS emerged from therealization that the signals could in fact be acquired directlyin their compressed form instead of sampling the signal atNyquist rate and then compressing it and while doing sodiscarding most of the laboriously obtained coefficients.Since the seminal works by Donoho [1] and C`andes,Romberg and Tao [2], [3], there has been an explosion ofresults in the field and many applications has been suggestedstarting with the prototype single pixel camera [4].Here we propose a new application of CS in acoustic wavefield sensing. In [5] the authors proved that the acoustic field isalmost optimally sparse in Curvelet frame. As the Curveletsessentially describe the wave front sets their propagation iswell approximated through geometrical optics (high frequencyasymptotic solution to the wave equation). The wave frontsets, and hence the Curvelets, propagate along the geometricalrays which are projections on R d of phase space solutions M. Betcke and S. Arridge are with the Department of Computer Sci-ence, University College London, UK, WC1E 6BT London, UK, e-mail:[email protected]. Cox and N. Huynh and E. Zhang and P. Beard are with Department ofMedical Physics, University College London, UK, WC1E 6BT London, UK. of the corresponding Hamilton-Jacobi equation resulting fromthe high frequency asymptotic. Motivated by this result weinvestigate the Curvelet representation of the cross-section ofthe wave field by the planar ultrasound sensor. While the argu-ments in [5] do not directly apply to this situation, the planarcross-section through the acoustic wave front constitutes asingularity along a smooth curve for which Curvelets havebeen demonstrated to be a nearly optimal representation [6].
A. Contribution
In this paper we focus on the reconstruction problem viadata recovery for a novel way of interrogating the highresolution ultrasound sensor using patterns instead of themore conventional sequential point-by-point interrogation. Anexample of such a system using a Single-Pixel Optical Camera(SPOC) was presented in [7]. The theory of CS predicts thatsubstantially fewer such measurements, of the order k log n ,need to be taken in order to capture a signal of length n andsparsity k leading to a substantial reduction of the acquisitiontime. We propose from such measurements to recover thepressure at the detector at each time step independently usingCS recovery algorithms. We discuss the motivation for usinga Curvelet tranform as the sparsifying transformation for theacoustic field at the detector and we propose its modification:a low-frequency Curvelet tranform tailored to the frequencyrange of the acoustic field on the detector. We investigatethe appropriate choice of interrogation patterns and recoveryalgorithm for this problem. The proposed techniques can beused in various applications ranging from ultrasonic fieldmapping to photoacoustic tomography (PAT). In this work wefocus on the latter to illustrate our method.For PAT applications, the immediate benefit of recoveringdata independently at each time step is the decoupling of theCS reconstruction from the acoustic inversion, allowing forthe recovered time series to be input into any reconstructionalgorithm for PAT for which software is readily available. B. Related work
For PAT, one and two step approaches have been investi-gated in the literature. In one step approaches the initial pres-sure is directly recovered from the compressed measurements.In [8] the authors propose such an approach for a 2D problemwith measurements limited in angle and frequency. It wasnumerically investigated using an analytic approximate inver-sion formula and a number of sparse representations including a r X i v : . [ m a t h . NA ] S e p OURNAL OF L A TEX CLASS FILES, VOL. 11, NO. 4, DECEMBER 2012 2
Wavelets and Curvelets. An approach for 3D imaging using atranslatable circular detection array and Wavelets as a sparserepresentation was presented in [9]. A generic variationalapproach for 2D and 3D PAT reconstruction based on analgebraic adjoint was first discussed in [10]. Two variationalapproaches based on the analytical adjoint were recentlyproposed, one in a FEM-BEM setting [11] and the other in a k -space setting with an efficient k -Wave implementation [12],[13].In contrast, in the context of PAT, the present paper isabout a two step method. We propose to first recover thephotoacoustic data in every time step independently frompattern measurements using sparsity of the data in Curveletbasis, and subsequently to reconstruct it using standard PATreconstruction methods. Recently, another two stage approachhas been suggested in [14], where the authors explore thetemporal sparsity of the data by means of a custom madetransform in time. From the sparse pressure data subsequentlythe initial pressure is reconstructed with the Universal BackProjection formula. C. Outline
The remainder of this paper is organized as follows. InSection II we introduce the forward and inverse problemfor PAT and methods for their solution. Section III brieflyrecapitulates the theory of compressed sensing. An appro-priate multiscale representation of the time series PAT datais considered in Section IV, where we derive the frequencymodel of sensor data and propose a modified version ofCurvelet transform tailored to the range of frequencies ofthe acoustic field. In Section V we discuss specific issuesarising when compressively sensing photoacoustic data usingpatterned interrogation of optical ultrasound detector. Webriefly describe the Single-Pixel Optical Camera based PATscanner. We consider the challenges for Curvelet transformsfor approximation of sensor data over the time series. Wediscuss choice of the interrogation patterns and the algorithmfor recovery of the sensor data. In Section VI we presentrecovery results for the optical sensor data and the final PATimage reconstruction from both simulated and real data.II. P
HOTOACOUSTIC T OMOGRAPHY
Photoacoustic tomography (PAT) is an example of a widerrange of hybrid imaging techniques, in which contrast inducedby one type of wave is read out by another wave. In this way,both high contrast and high resolution can be simultaneouslyachieved, which is often difficult with conventional imagingtechniques that usually provide either one or the other, butnot both. PAT is an emerging biomedical imaging modalitywith both pre-clinical and clinical applications that can providecomplementary information to established imaging techniques[15]–[19]Many PAT applications require a high resolution, threedimensional image e.g. an image of capillaries of a few tens ofmicrons diameter in a cm sized imaging region. Such highlyresolved imaging requires an ultrasound sensor array of tensof thousands of pixels. In one such PAT system [20], the sensor is a Fabry Perot (FP) interferometer interrogated bya laser whose focus is moved to form a raster scan of thedesired resolution. For sequential sampling, such as this, theultimate limit to the data acquisition rate (the rate at whichtime series are collected) is the propagation time for soundto cross the specimen, e.g. it would take µs for the signalto reach the detector from 15mm depth, resulting in kHz acquisition rate (not to be confused with the sampling rate,which might be as high as 100 MHz). No sequential scanner isclose to approaching this limit, making a sequential acquisitiona major practical limitation for high resolution 3D PAT. For in vivo applications the required acquisition times at currentlyachievable rates are typically a few minutes, not only resultingin motion artifacts, but limiting studies to phenomena on suchlong timescales.The principle involved in PAT is to send a short (ns)pulse of near-infrared or visible light into tissue, whereuponabsorption of the photons e.g. by haemoglobin molecules,generates a small local increase in pressure which propagatesto the surface as a broadband, ultrasonic pulse. If the amplitudeof this signal is recorded over an array of sensors at thetissue surface, an image reconstruction algorithm can be usedto estimate the original 3D pressure increase due to opticalabsorption; this is the photoacoustic image, p .Mathematically, under an assumption of free space propa-gation the photoacoustic forward problem is modelled as aninitial value problem for the wave equation [21] c ( x ) ∂ p ( x , t ) ∂t = ρ ( x ) ∇ · (cid:18) ρ ( x ) ∇ (cid:19) p ( x , t ) , x ∈ R d , t ∈ (0 , T ) , (1a) p ( x ,
0) = p ( x ) , (1b) ∂∂t p ( x ,
0) = 0 , (1c)where p ( x , t ) denotes the time dependent acoustic pressure in R d × (0 , T ) , d = 2 , , p ( x ) its initial value and c ( x ) and ρ ( x ) are the ambient speed of sound and density, respectively.The photoacoustic inverse problem is to recover this initialpressure p ( x ) , x ∈ Ω compactly supported in the regionof interest Ω from a time series measurement g ( x , t ) = p ( x , t ) , x ∈ S , t ∈ (0 , T ) on the surface S (e.g. boundaryof Ω ) and it amounts to a solution of the following initialboundary value problem [22] c ( x ) ∂ p ( x , t ) ∂t = ρ ( x ) ∇ · (cid:18) ρ ( x ) ∇ (cid:19) p ( x , t ) , x ∈ Ω , t ∈ (0 , T ) , (2a) p ( x ,
0) = , (2b) ∂∂t p ( x ,
0) = 0 , (2c) p ( x , t ) = g ( x , T − t ) , x ∈ S , t ∈ (0 , T ) , (2d)with PAT data fed backwards in time as boundary values (alsoreferred to as time reversal ). The formulation (2) holds exactlyin 3D for non-trapping smooth sound speed c ( x ) if T hasbeen chosen large enough so that g ( x S , t ) = 0 , t ≥ T andthe wave has left the domain Ω . Furthermore, assuming thatthe measurement surface S surrounds the region of interest OURNAL OF L A TEX CLASS FILES, VOL. 11, NO. 4, DECEMBER 2012 3 Ω containing the support of the initial pressure p , problem(2) has a unique solution. The condition on S to surround p can be relaxed under additional assumptions on S and Ω and smoothess of initial pressure p ∈ H (Ω) ; see [23]and the citations within. The stability of reconstruction canbe obtained from microlocal analysis which provides insightsinto which singularities are visible and which are not [24],[25], which directly translates to stability of reconstruction ofthose singularities [23].In a real experiment, we can only acquire a discrete (bothin time and space) subset of time series measurements, andfrequently it is not possible to acquire measurements on asurface S surrounding the object. An example of a popularsensor violating this assumption is a planar sensor. In practice,such sensor will have finite size, resulting in invisibility ofsome interfaces and in turn in artefacts in the reconstructedimage.We should mention that other approaches to image recon-struction exist. An overview of methods for the case when S is a surface surrounding Ω can be found e.g. in [23].In this paper we directly solve the time-reversal approach(2) using the pseudospectral method implemented in the k -Wave Toolbox [26], which is an efficient numerical schemefor solving the wave equation in domains with heterogeneousacoustic properties, and is exact in the case of homogeneousmedia. Furthermore, the methodology proposed here is tailoredto high resolution detectors which are planar and so we assume S to be a finite rectangular section of the xy -plane S = { ( x, y, z ) : | x | ≤ x d / , | y | ≤ y d / , z = 0 } . and solve the resulting initial partial boundary value problem(2). The rectangular planar detector shape allows us to use anysparsifying transform derived for natural images.III. C OMPRESSED SENSING
Let Ψ be a sparsifying transform Ψ : R n → R N , resulting in possibly an overdetermined representation N ≥ n and Ψ − its inverse (in the frame sense). With f we denotethe transformation of the original signal g ∈ R n f = Ψ g. (3)In compressed sensing the signal g ∈ R n is projected on aseries of sensing vectors φ i ∈ R n , i ∈ , . . . , m , where m (cid:28) n , yielding a vector b ∈ R m of compressive measurements b = Φ g + e with Φ = [ φ , φ , . . . , φ m ] T , (4)where e is the measurement noise, (cid:107) e (cid:107) ≤ ε . The followingrecovery result guarantees that the original signal g , can berobustly recovered from compressed measurements (4) viasolution of the minimization problem for f [3] min f ∈ R n (cid:107) f (cid:107) , s.t. (cid:107) ΦΨ − f − b (cid:107) ≤ ε. (5) Theorem RR (Robust recovery [27], [28]) . Let δ k be theisometry constant of ΦΨ − defined as a smallest positivenumber such that (1 − δ k ) (cid:107) f (cid:107) ≤ (cid:107) ΦΨ − f (cid:107) ≤ (1 + δ k ) (cid:107) f (cid:107) (6) holds for all k -sparse vectors f .If δ k < √ − (relaxed to δ k < √ / in [28]) thenthe error of solution to (5) f ∗ is bounded as follows (cid:107) f − f ∗ (cid:107) ≤ C ε + C (cid:107) f − f k (cid:107) √ k , (7) where f k denotes best k -term approximation, obtained from f selecting its k largest in magnitude coefficients, and C , C are constants dependent only on δ k . The robust recovery Theorem RR holds for any vector f ,however the error norm (cid:107) f − f k (cid:107) of the k -term approxi-mation is only small for k -sparse or compressible vector f (i.e. with fast enough decaying magnitude of the coefficients).Furthermore, for k -sparse vectors we have ( (cid:107) f − f k (cid:107) ) / √ k ≤(cid:107) f − f k (cid:107) and hence the the bound can be expressed in L norm (cid:107) f − f ∗ (cid:107) ≤ C ε + C (cid:107) f − f k (cid:107) . IV. M
ULTISCALE REPRESENTATION OF TIME SERIES DATA
A. Frequency model of sensor data
For constant sound speed and density c ( x ) = c , ρ ( x ) = ρ ,the wave equation (1) becomes homogeneous and it admits ananalytic solution in Fourier domain [29] p ( k , t ) = cos( c | k | t ) p ( k ) , (8)where k is the frequency domain wave vector. Equation (8)formally allows us to calculate the PAT data time series p ( x S , t ) = F − ( p ( k , t )) (cid:12)(cid:12) x = x S (9) = F − (cos( c | k | t ) p ( k )) (cid:12)(cid:12) x = x S = 2 (cid:90) ∞ F − (cid:107) (cos( c | k | t ) p ( k )) d k ⊥ , where we used the tensor decomposition F = F (cid:107) F ⊥ of the3D Fourier transform in, and orthogonal to, the sensor plane,and the fact that for x S we have z = 0 and hence e i k ⊥ z = 1 .Since p ( x ) is a real valued function, p ( k ) and conse-quently p ( k , t ) are symplectic in k i.e. p ( k ) = p ( − k ) ∗ .Furthermore, p ( k , t ) as given by (8) is even in t . As ourdomain Ω is positioned in z ≥ to restore uniqueness weassume p ( x ) to be even w.r.t. z = 0 plane. Consequently, theFourier transform in the z as well as the t direction is equal tothe cosine transform and can be calculated integrating from 0to ∞ ; thus in what follows the two are used interchangeably.In experiments however, p ( x S , t ) is available to us onlyindirectly, through measurements. To account for limitationsof physical equipment, we introduce a degradation operator, D . Note, that here D is not a measurement operator, but itmodels the effect of the finite size and temporal response ofthe measurement system which reads the acoustic time series.In other words, g ( x S , t ) = D p ( x S , t ) describes the filtering of OURNAL OF L A TEX CLASS FILES, VOL. 11, NO. 4, DECEMBER 2012 4 the acoustic pressure p ( x S , t ) by the physical presence of thesensor before the sensing vector, φ i is applied to it to collectthe compressive measurement, b i .In this work we assume that D is a band limiting spatialand temporal filtering operator acting on the time dependentpressure on optical sensor, which for simplicity we describein the frequency domain D [ p ]( k S , ω ) = w t ( ω ) w (cid:107) ( k S ) p ( k S , ω ) , (10)where w (cid:107) ( k S ) and w t ( ω ) are some frequency window func-tions on the sensor and in time.Taking Fourier transform of (9) in both variables we obtain p ( k S , ω ) = F t (cid:0) F (cid:107) ( p ( x S , t )) (cid:1) = F t (cid:16) (cid:82) ∞ F (cid:107) (cid:16) F − (cid:107) (cos( c | k | t ) p ( k )) (cid:17) d k ⊥ (cid:17) = F t (cid:18) (cid:82) ∞ ω/c cos( ωt ) √ ( ω/c ) −| k S | p ( k S , (cid:112) ( ω/c ) − | k S | ) dω (cid:19) = ω/c √ ( ω/c ) −| k S | p ( k S , (cid:112) ( ω/c ) − | k S | ) , (11)where changing the integration variable to ω , ω/c > | k S | , inthe third line allowed us to interpret the integral as inversecosine transform in ω . Equation (11) connects the Fouriertransform of the pressure time series on the detector withthe Fourier transform of initial pressure and is the basis ofthe reconstruction formula derived in [30], [31]. From (11) itis obvious that application of D [ p ]( k S , ω ) corresponds to theapplication of a filter window w t ( c | k | ) w (cid:107) ( k S ) to p ( k ) D [ p ]( k S , ω ) = w t ( ω ) w (cid:107) ( k S ) ω/c √ ( ω/c ) −| k S | p ( k S , (cid:112) ( ω/c ) − | k S | )= w t ( c | k | ) w (cid:107) ( k S ) | k || k ⊥ | c p ( k ) =: D Ω [ p ]( k ) for scalar k ⊥ > . Thus only a smoothed initial pres-sure can be recovered from the data. Conversely, whensimulating PAT data, the initial pressure can be smoothedwith w t ( c | k | ) w (cid:107) ( k S ) before forward propagation instead ofsmoothing the sensor data with D , which can be useful toeliminate Gibbs phenomena in k -space methods; see also [32],[33] where a Blackman window was applied to p . B. Curvelets
In this work we suggest using Curvelets to represent eachof the measured time series data g ( x S , t ) . For the planarsensor S this corresponds to a planar cross-section of thewave field p ( x , t ) at a given time t (which we observe over afinite section of the plane, S ). Due to PAT forward problembeing an initial value problem, the corresponding wave field p ( x , t ) is essentially smooth away from the p -shaped (withcorners smoothed) wave front and the same holds for theirplanar cross-sections. As the time evolves p ( x , t ) developsoverlapping wave fronts. Those however can be treated assuperposition and hence we only need to be able to representan individual wave front.In what follows we assume g ( x S , t ) to be sampled on an n × n grid on the sensor, resulting in an image g t [ i , i ] , i =1 , . . . , n i = 1 , . . . , n , where ( · ) denotes function evalua-tion and [ · ] array indexing (following the original notation in[34]). Curvelets [34], [35] are a multiscale pyramid with manydirections and positions at each scale. Curvelets obey parabolicscaling, meaning that at scale − j Curvelet has an envelopewhich aligns along a ridge of length − j/ and width − j ,resulting in a location, direction and scale dependent frequencyplain tiling. Curvelets have been shown to be nearly optimalrepresentation of objects smooth away from (piecewise) C singularities [36], [37] and hence the error of the k -termcurvelet approximation (corresponding to taking the k largestmagnitude coefficients) can be bounded as [6] (cid:107) g − g k (cid:107) ≤ C · (log k ) · k − . (12)Computation of Curvelets at the finest scale, − J , is notstraight forward because in order to capture the direction ofthe wave it is necessary that the Curvelet is sampled morefinely than the maximal frequency in the image. As a resultthe frequency domain support of a Curvelet at the finest scale J with the direction θ (cid:96) , ˜ U J,(cid:96) , does not fit into the fundamentalcell [ − n / , n / − × [ − n / , n / − . One solutiongiven in [34] is to wrap the Fourier transform back onto thefundamental cell which effectively periodizes it ˜ U J,(cid:96) [ (cid:0) i + n (cid:1) mod n − n , (cid:0) i + n (cid:1) mod n − n ]= ˜ U J,(cid:96) (2 πi , πi ) , where the indices i , are chosen in the support of ˜ U J,(cid:96) ( · , · ) .The periodization in Fourier domain corresponds to undersam-pling in space, which in turn causes aliasing. In [34] this effectwas shown to account for less than 10% of the squared normof the coefficients.After periodization, the Fourier transform is multiplied by a C ∞ partition of unity window i.e. a window which weights theoriginal frequency and its periodic extension such that the sumof their squares is equal to 1. This acts to preserve the normof the signal throughout periodization thereby guaranteeingthat the computed transform is a numerical isometry. The C ∞ window used in [34] is a tensor product of the following onedimensional C ∞ windows w [ i ] = a ( x i ) h ∞ (1 − x i ) , i = 1 , . . . , m , , i = m + 1 , . . . , m + 1 ,a ( x i ) h ∞ ( x i − m − ) , i = 3 m + 2 , . . . , m + 1 , (13)where x i = ( i − / ( m − , i = 1 , . . . m , m = (cid:98) n / (cid:99) , h ∞ is a C ∞ monotonically decreasing function and a ( x ) = ( h ∞ (1 − x ) + h ∞ ( x ) ) − / (14)is a normalizing factor. Expression (13) assumes n mod 3 =0 , and hence m = n / resulting in a window of length m + 1 , while corresponding windows can be derived in thecases n mod 3 = 1 , . C. Low-frequency Curvelets
The sampled time series data, g t , is band limited with thehigh frequencies roll off (see Section IV-A) in contrast tostandard images which are “sharp” i.e. high frequencies arenot damped. This raises the question, if we should representthose rolling off frequencies in the same way as the undamped OURNAL OF L A TEX CLASS FILES, VOL. 11, NO. 4, DECEMBER 2012 5 frequencies, as the standard Curvelets would do at a cost whichscales with n = n n as n log n .As a compromise, we propose here to use a low-frequencyCurvelet representation for the high resolution image g t ,i.e. to compute Curvelets corresponding to a subset of thefrequencies excluding the Cartesian annulus of the high-est frequencies, [ − n LF / , n LF / − × [ − n LF / , n LF / − , n LFl < n l , l = 1 , . The high frequencies in the annulus [ − n / , n LF / − ∪ [ n LF / , n / − × [ − n / , n LF / − ∪ [ n LF / , n / − present in the high resolution image g t can then be used to fill up the frequency range needed forcomputation of the low-frequency Curvelets at the finest scale,instead of having to periodize. In this way, we avoid undersam-pling of the low-frequency Curvelets at the finest scale whichnow however corresponds to frequencies lower than when theCurvelets are computed for the full set of frequencies. Weare going to concentrate on the implementation via wrappingbecause in this case the original Curvelet transform for thefull set of frequencies results in a numerical isometry.In what follows, we are going to denote the standardCurvelet transform for an n × n image with J scales as C n ,n J and the low-frequency Curvelet transform correspond-ing to frequencies up to n LF / × n LF / (we will use suchsimplified notation to denote both the negative and positivefrequencies), n LFl < n l , l = 1 , for the same image with C n ,n J,n LF ,n LF .Let ˜ U LFJ,(cid:96) = [ − N LF / , N LF / − × [ − N LF / , N LF / − with N LFl = 2 (cid:98) m LFl (cid:99) +1 , m LFl = n LFl / , l = 1 , denotethe frequency domain support of the finest scale low-frequencyCurvelet. In general we have two cases:i) N LFl / ≤ n l / , l = 1 , : the frequency support of thefinest scale low-frequency Curvelet ˜ U LFJ,(cid:96) is contained inthe fundamental cell [ − n / , n / − × [ − n / , n / − ;ii) N LFl / > n l / , l = 1 , : the finest scale low-frequencyCurvelet support ˜ U LFJ,(cid:96) extends to frequencies outside ofthe fundamental cell.In the first case, to compute the standard Curvelets corre-sponding to image with maximal frequency n LF / × n LF / , C n LF ,n LF J , one would compute the Fourier transform of g t , re-strict it to frequencies up to n LF / × n LF / and subsequentlyperiodize it while weighting with the C ∞ window (13).Instead, to compute the low-frequency Curvelets, C n ,n J,n LF ,n LF ,we apply to the full Fourier transform of g t (with frequenciesup to n / × n / ) a rectangular low-pass window with a cut-off frequency N LFl / , l = 1 , being the highest frequency inthe support ˜ U LFJ,(cid:96) . Then we proceed as for standard Curvelets C n LF ,n LF J to compute the finest scale coefficients.In the second case, to compute the low-frequency Curvelets, C n ,n J,n LF ,n LF , we still need to use priodization to fill up the sup-port ˜ U LFJ,(cid:96) . However, we only periodically extend the Fouriertransform of g t beyond the range of available frequencies,here n / × n / , rather then n LF / × n LF / as would bedone using standard Curvelets, C n LF ,n LF J . This results in a C ∞ window with shorter (and steeper) flanks as we only need tofill in the frequencies in the range [ − N LF / , − n / − ∪ ψ (64,43)80 100 120 140 160 18080100120140160180 wave front diection (a) ψ ∈ C , ψ (48,33)80 100 120 140 160 18080100120140160180 wave front direction (b) ψ ∈ C , , , (c) log | ψ | , ψ ∈ C , (d) log | ψ | , ψ ∈ C , , , Fig. 1. Zoomed in image of (a) standard Curvelet C , , (b) low-frequency Curvelet C , , , , at the finest level 3. (c), (d) Surface plot oflog amplitude of (a), (b). [ n / , N LF / − × [ − N LF / , − n / − ∪ [ n / , N LF / − . The corresponding window is a tensor product of thefollowing one dimensional C ∞ windows w LF [ i ] = a ( x i ) h ∞ (1 − x i ) , i = 1 , . . . , N LF − n , , i = N LF − n + 1 , . . . , n ,a ( x i ) h ∞ ( x i − n ) , i = n + 1 , . . . , N LF , (15)where x i = ( i − / ( N LF − n − , i = 1 , . . . , N LF − n ,and a and h ∞ are as before.Subsequently, all the lower scales of low-frequencyCurvelets can be computed exactly as for the standardCurvelets C n LF ,n LF J . The construction of low-frequencyCurvelets above results in a numerical isometry on the re-striction to frequency range up to min { N LF / , n / } × min { N LF / , n / } . Thus low-frequency Curvelet transformis a different transformation to the original Curvelet transformwhich however for a choice of n LFl / such that N LFl / ≥ n l / , l = 1 , , is also an isometry on [ − n / , n / − × [ − n / , n / − as the standard Curvelets C n ,n J . Figures1(a),(b) show the finest scale Curvelet and low-frequencyCurvelet, respectively. The plots of the logarithm of the am-plitude in Figures 1(c),(d) illustrate that the decay of the low-frequency Curvelets along the wavefront direction is slowerthan that of the standard Curvelets, which is due to thefrequencies higher then n LFl / , l = 1 , used in computationof Curvelets at the finest scale.Recapitulating, there are two major benefits of such low-frequency Curvelet transform. First, is the super linear reduc-tion in computation cost of the transform, which is particularlybeneficial for solution of the CS recovery problem (5) whichinvolves repeated application of ΦΨ T and its adjoint ΨΦ T .Second, is the ability to effectively represent realistic PAT data,in which due to measurement process the high frequenciesare damped as in the model derived in Section IV-A. This OURNAL OF L A TEX CLASS FILES, VOL. 11, NO. 4, DECEMBER 2012 6 results in amplitudes of low-frequency Curvelet coefficientsbeing higher and exhibiting a quicker decay than those of thestandard Curvelet coefficients (see Figure 2 and the accompa-nying discussion in Section VI-A), and consequently in higherrobustness to noise and imperfect compressibility.V. C
OMPRESSED SENSING OF OPTICAL ULTRASOUNDDETECTOR
A. Single-pixel optical camera
In a series of publications [7], [38], [39] we introduced asingle-pixel optical camera (SPOC) for ultrasonic and pho-toacoustic imaging. With the SPOC, instead of recording thepressure on the detector point-by-point, the entire active areaof the optical ultrasound sensor is illuminated, and using adigital micro-mirror device (DMD) a pattern φ j is applied tothe wide field light reflected from the sensor resulting in acompressed measurement b j = φ T j g ( x S , t i ) , t i ∈ (0 , T ) . (16)Figure 3 shows a sketch of the operational principle of SPOCwhile for technical details of the system we refer to [7]. B. Sparse representation of the sensor data
In PAT the entire time series for one point or pattern isacquired with one excitation, φ T j g ( x , t i ) , x ∈ S , t i ∈ (0 , T ) .This has the consequence that we acquire the same numberof compressed/point measurements of the wave field at eachtime step t i . Furthermore, due to the rate at which the DMDcan change patterns, at least at present, we are limited to useof only one pattern throughout the acoustic propagation.However, as the wave propagates, the complexity of thesensor data varies with time. After a wave front reaches thedetector, its cross-section with the detector plane expands andmore wave fronts, corresponding to features farther away,reach the detector. As a result, in general the complexity ofthe wave field at the sensor over time first increases and atsome point it starts to decrease again corresponding to the tailof the wave field. If the complexity is reflected by the sparsity,the error of the best k -term approximation for the sensor datavaries throughout the time series and with it the recoveryerror bound in the robust recovery Theorem RR. Furthermore,due the sound intensity, which is ∝ p ( x , t ) , obeying theinverse square law, p ( x , t ) decays as inverse distance of x to the source. This means that as with increasing t , p ( x S , t ) encodes information about p further away from the sensor, therecorded pressure amplitudes decrease linearly for the sameinitial pressure amplitude value, ultimately resulting in a lowerpoint wise signal to noise ratio for longer propagation times. C. Sensing patterns
In our experiments we used scrambled Hadamard patterns H sj = P r H j P c , (17)where H j is the j × j Hadamard matrix and P c , P r ∈{ , } j × j denote permutation matrices for columns androws, respectively. For compressed sensing we only select first m (cid:28) n = 2 j rows of H sj . The application of H sj to a vector H sj v = P r H j P c v amounts to application of thepermutation matrix P c to v , performing Hadamard transformon the permuted vector P c v , and subsequently permutingthe rows. Thus the scrambled Hadamard transform can becomputed at essentially the same cost as the fast Hadamardtransform, while scrambled Hadamard matrices have recoveryproperties similar to those of random Bernoulli matrices seee.g. [40].In practice using the DMD, it was only possible to apply bi-nary patterns. In order to make use of properties of Hadamardtransform such as othogonality and self inversion, the experi-mental Hadamard matrix H (0 , j needs to be transformed intothe Hadamard matrix H j using the simple relation H j = (cid:16) H (0 , j − T (cid:17) / √ j , (18)and correspondingly the measured data w = H (0 , j g into H j g = 1 √ j (cid:16) H (0 , j − T (cid:17) g = 1 √ j (2 w − w (1)) . (19)Here, we used that w (1) = T g corresponds to the measure-ment acquired with ‘all-1s’ pattern, which is the first row of H (0 , . It is immediately clear that the same transformation(18), (19) holds for scrambled Hadamard matrices (17) (with1 in w (1) replaced by the ‘all-1s’ row number after rowpermutation).The light reflected from the DMD is integrated by a photo-diode. In order to best utilize the dynamic range of the photodiode, it is necessary to keep the optical power incident onthe photodiode in the same range for each pattern. All but the‘all-1s’ pattern are composed of an equal number of 0 and1s. Therefore, the ‘all-1s’ pattern was replaced with a vectorwhich first half entries are 0 and the second half 1. As thenegative of this vector (0 becomes 1 and vice versa) is exactlythe j / row of H (0 , j , the data corresponding to the‘all-1s’ pattern can be constructed by adding the data fromthe modified first row pattern and (2 j / th row pattern.Again, this construction is not affected by scrambling. D. Recovery of the sensor data
The sampled PAT data at each time step, g ( x S , t ) , canbe recovered by solving the optimization problem (5). Thesensing matrix Φ in (5) is set to be the first m rows of thescrambled Hadamard matrix H s log ( n ) , n = n n , Φ = SH s log ( n ) = SP r H log ( n ) P c , (20)where S ∈ { , } m × n is a binary subsamplig matrix suchthat S T S is a n × n diagonal matrix with ones at positionscorresponding to the chosen and zeros to the ommitted rows,respectively, while SS T = I m × m is an m × m identity matrix.Consequently, we have ΦΦ T = I m × m .The vector of measurements, b , is computed from expe-rimantal measurements using equation (19). The sparsifyingtransform Ψ is chosen as an orthonormal low-pass or standardCurvelet transform on the n × n optical sensor S and f t isthe (sparse) vector of sought for coefficients, f t = Ψ g t . OURNAL OF L A TEX CLASS FILES, VOL. 11, NO. 4, DECEMBER 2012 7
To take advantage of the structure of the problem wesolve (5) using the Split Augmented Lagrangian ShrinkageAlgorithm (SALSA) proposed in [41]. SALSA is an ADMMscheme which solves the unconstrained problem min f ζ ( f ) := 12 (cid:107) ΦΨ T f − b (cid:107) + τ (cid:107) f (cid:107) . (21)A version used in this work is summarized in Algorithm 1. Algorithm 1
Split Augmented Lagrangian Shrinkage Algo-rithm (SALSA), [41]. Choose µ > , v and d i := 0 repeat f i +1 = arg min f (cid:107) ΦΨ T f − b (cid:107) + µ (cid:107) f − v i − d i (cid:107) v i +1 = arg min v τ (cid:107) v (cid:107) + µ/ (cid:107) f i +1 − v − d i (cid:107) d i +1 = d i − ( f i +1 − v i +1 ) i = i + 1 until | ζ ( f i +1 ) − ζ ( f i ) | /ζ ( f i ) < ε The quadratic minimization problem in line 4 leads to thelinear system f i +1 = (ΨΦ T ΦΨ T + µI ) − (ΨΦ T b + µ ( v i + d i )) , (22)where using Sherman Morrison Woodbury formula the inversecan be expressed in terms of application of ΦΨ T and its adjoint (ΨΦ T ΦΨ T + µI ) − = 1 µ (cid:18) I − µ + 1 ΨΦ T ΦΨ T (cid:19) . (23)The proximal operator in line 5 amounts to point wise softthresholding T τ/µ ( x ) = sign ( x )( | x | − τ /µ ) + . (24)VI. PAT RECONSTRUCTION WITH THE RECOVEREDSENSOR DATA
We present reconstruction results for three photoacousticimage reconstruction problems, reconstruction from: simulatedpattern data, synthesized pattern data from a point-by-pointscanner data, and data acquired with the SPOC (Section V-A).After the PAT data has been recovered using the proposedmethod for acoustic field reconstruction, the PAT images arereconstructed using time reversal via first order method in k -Wave Toolbox. A. Simulated data: clock phantom
In our first example we simulate the PAT data for initialpressure distribution p depicted in Figure 4. The “clockphantom” is a collection of balls which give rise to sphericalwave forms expanding uniformly in all directions, which is adifficult case for directional basis like Curvelets.We consider a volume of × × voxels of size . mm and the sensor of matching resolution, × ,placed at z = 0 . We assume homogeneous ambient speed ofsound and density of m/s and kg/m , respectively.The pressure at the sensor is sampled every ns. Both, theforward and inverse PAT problems are solved with first order method from k -Wave Toolbox. Before the PAT time series issimulated, p is smoothed with the 3D Blackman window.We attempt to recover the PAT data from 18% of noise-less compressed measurements obtained with binary scram-bled Hadamard patterns. We solve the recovery problemin each time step using SALSA with the data dependentregularization parameter τ = 0 .
01 max( | ΨΦ T b t | ) and µ =5 max( | ΨΦ T b t | ) / (cid:107) b t (cid:107) . We stop the algorithm if the relativechange in the objective function ζ ( f ) (21) drops below · − or after iterations.We start by examining the utility of the Curvelet transformas a sparsifying transform for PAT data recovery problem.In particular we compare the standard Curvelet transform C , with the low-frequency Curvelet transform C , , , ,the lowest resolution transform which is still an isometry onthe full frequency range (up to / × / ).First, we consider the approximation properties of thestandard Curvelet and low-frequency Curvelet transforms.Figure 2(a-b) shows the decay of the amplitudes of Curveletcoefficients for standard Curvelets C , and low-frequencyCurvelets C , , , of the photoacoustic field at time steps t i = 100 , corresponding to the fields depicted in Figures6(a), 8(a). For the largest amplitude coefficients the low-frequency Curvelets consistently exhibit a quicker decay andcorrespondingly lower approximation error shown in Figure2(c-d). While in the early time steps t i = 100 the decay ratedifference is most pronounced for later time steps the dif-ference gets smaller but is distributed over more coefficients.In all cases eventually the approximation error of standardCurvelet transform falls below that of the low-frequencyCurvelet transform but this is only after all the significantcoefficients has been captured. −3 −2 −1 log k Ψ g ( x S ,t ) k ψ , ψ ∈ C ψ , ψ ∈ C ψ , ψ ∈ C ψ , ψ ∈ C ψ , ψ ∈ C ψ , ψ ∈ C (a) log | Ψ g ( x S , t ) | −3 −2 −1 log k Ψ g ( x S ,t ) k ψ , ψ ∈ C ψ , ψ ∈ C ψ , ψ ∈ C ψ , ψ ∈ C ψ , ψ ∈ C ψ , ψ ∈ C (b) log | Ψ g ( x S , t ) | −1 log kk Ψ g ( x S ,t ) − Ψ ˆ g ( x S ,t ) kk Ψ = C Ψ = C (c) Error of ˆ g ( x S , t ) −1 log kk Ψ g ( x S ,t ) − Ψ ˆ g ( x S ,t ) kk Ψ = C Ψ = C (d) Error of ˆ g ( x S , t ) Fig. 2. Clock phantom. The decay of log amplitude of Curvelet coefficients of g ( x S , t i ) at time steps (a) t , (b) t ; the colors correspond to coefficientsat different scales. (c), (d) The corresponding compression error in log scale. Next, we compare the compression with of coefficientto recovery from of measurements which correspondsto times the assumed sparsity, an empirically chosen factorfor binary scrambled Hadamard patterns of size . Figure5 shows the mean square error (MSE) of the compressed OURNAL OF L A TEX CLASS FILES, VOL. 11, NO. 4, DECEMBER 2012 8 versus recovered PAT data for both transforms. While theMSEs of the reconstructed data are almost identical, forthe data compressed with the low-frequency Curvelets, aftersome initial time steps the MSE becomes lower than for thestandard Curvelets, and most importantly it matches closerthe MSE of the reconstructed data. This demonstrates thatthe low-frequency Curvelets are an adequate (while cheaper)representation of the PAT data. The series of Figures 6,7,8,9shows the PAT data g ( x S , t ) , its compression ˆ g ( x S , t ) andreconstruction ˜ g ( x S , t ) at different time steps. Consistently,we observe that the higher scale coefficients are eliminated bycompression while they partially reappear in the reconstructionhinting that the factor 6 maybe somewhat pessimistic. This isalso evident in the higher frequency appearance of the errorof the reconstruction in comparison to the compression.The PAT image reconstruction from the recovered PAT datafor both transforms is depicted in Figure 10, which are visuallyvery similar with MSE of . · − for the standardCurvelets and . · − for the low-frequency Curvelets,where the reconstruction from full data has been used as theground truth in MSE calculations. Fig. 3. Single-pixel optical camera.
Fig. 4. Clock phantom. −10 −5 t MSE(ˆ g ( x S ,t ))MSE(˜ g ( x S ,t )) (a) Ψ = C , −10 −5 t MSE(ˆ g ( x S ,t ))MSE(˜ g ( x S ,t )) (b) Ψ = C , , , Fig. 5. MSE of the compressed PAT data ˆ g ( x S , t ) versus the reconstructedPAT data ˜ g ( x S , t ) over time, for (a) standard Curvelet transform, C , ;(b) low-frequency Curvelet transform, C , , , . B. Synthesized pattern data: knotted tubes filled with ink
Next, we present reconstruction from compressed measure-ments synthesized from the point-by-point FP sensor mea-surements. The purpose of such data is to demonstrate thereconstruction with realistic noise but good signal to noiseratio, which at present is a limiting factor for SPOC.Two polythene tubes were filled with 10% and 100%ink and tied into a knot, see Figure 11(a). The tubes wereimmersed in a % intralipid solution. The wavelength ofthe excitation laser was nm delivering energy of ap-proximately mJ. A full scan data consists of × locations corresponding to spatial resolution of µm × µm, sampled at time points corresponding to timeresolution of ns. g ( x S ,t )
50 100 150 200 25050100150200250 00.050.10.150.20.250.30.35 (a) g ( x S , t ) ˆ g ( x S ,t )
50 100 150 200 25050100150200250 (b) ˆ g ( x S , t ) ˜ g ( x S ,t )
50 100 150 200 25050100150200250 (c) ˜ g ( x S , t ) Ψ ( g ( x S ,t )) (d) Ψ g ( x S , t ) Ψ (ˆ g ( x S ,t )) (e) Ψˆ g ( x S , t ) Ψ (˜ g ( x S ,t )) (f) Ψ˜ g ( x S , t ) ˆ g ( x S ,t ) − g ( x S ,t )
50 100 150 200 25050100150200250 −0.01−0.00500.0050.01 (g) Error of ˆ g ( x S , t ) ˜ g ( x S ,t ) − g ( x S ,t )
50 100 150 200 25050100150200250 −15−10−505x 10 −3 (h) Error of ˜ g ( x S , t ) Fig. 6. Clock phantom. PAT data at time step t , (a) simulated fulldata g ( x S , t ) , (b) compressed ˆ g ( x S , t ) and (c) reconstructed ˜ g ( x S , t ) datawith standard Curvelet transform, C , . The corresponding Curveletcoefficients (d-f) and (g) compression error ˆ g ( x S , t ) − g ( x S , t ) , (h) recoveryerror ˜ g ( x S , t ) − g ( x S , t ) . As a gold standard we take the reconstruction from afull point-by-point data set shown in Figure 13(a). The lin-ear reconstruction from 25% of patterns obtained by set-ting the missing Hadamard pattern measurements to 0 isshown in Figure 13(b) while the nonlinear reconstructionwith SALSA (with parameters τ = 0 .
01 max( | ΨΦ T b t | ) , µ =5 max( | ΨΦ T b t | ) / (cid:107) b t (cid:107) and stopping tolerance · − or after iterations) in Figure 13(c). The nonlinear reconstructioneffectively restores the contrast lost in the linear 0-paddedreconstruction. The MSE with respect to the full data recon-struction is . · − for the linear 0-padded reconstructionand . · − for the nonlinear reconstruction. C. Experimental pattern data: hair knot
Our last example is a synthetic hair knot phantom ofdiameter ∼
150 µm, immersed in 1% intralipid solution andpositioned approximately mm above the sensor and mmdeep below the intralipid surface. The photo of the phantomis shown in Figure 11(b). The area of DMD corresponding to × micromirrors grouped in × was used to form × pixels. Due to the angle of the optical path, eachsuch pixel corresponds to an area of . × µm on the FPsensor. The measurement was averaged over four excitations.The speed of sound used for time reversal was 1490 m/s.Figure 14(a) shows the PAT image reconstruction obtainedfrom the full set of scrambled binary Hadamard pat-terns (the full point data was computed inverting scrambledHadamard transform). The ∼ . amplified linear reconstruc-tion from 0-padded 18% measurements is shown in Figure OURNAL OF L A TEX CLASS FILES, VOL. 11, NO. 4, DECEMBER 2012 9 g ( x S ,t )
50 100 150 200 25050100150200250 00.050.10.150.20.250.30.35 (a) g ( x S , t ) ˆ g ( x S ,t )
50 100 150 200 25050100150200250 (b) ˆ g ( x S , t ) ˜ g ( x S ,t )
50 100 150 200 25050100150200250 (c) ˜ g ( x S , t ) Ψ ( g ( x S ,t )) (d) Ψ g ( x S , t ) Ψ (ˆ g ( x S ,t )) (e) Ψˆ g ( x S , t ) Ψ (˜ g ( x S ,t )) (f) Ψ˜ g ( x S , t ) ˆ g ( x S ,t ) − g ( x S ,t )
50 100 150 200 25050100150200250 −0.01−0.00500.0050.01 (g) Error of ˆ g ( x S , t ) ˜ g ( x S ,t ) − g ( x S ,t )
50 100 150 200 25050100150200250 −0.015−0.01−0.00500.0050.01 (h) Error of ˜ g ( x S , t ) Fig. 7. Clock phantom. PAT data at time step t , (a) simulated full data g ( x S , t ) , (b) compressed ˆ g ( x S , t ) and (c) reconstructed ˜ g ( x S , t ) data withlow-frequency Curvelet transform, C , , , . The corresponding Curveletcoefficients (d-f) and (g) compression error ˆ g ( x S , t ) − g ( x S , t ) , (h) recoveryerror ˜ g ( x S , t ) − g ( x S , t ) . g ( x S ,t )
50 100 150 200 25050100150200250 −0.3−0.2−0.100.10.20.3 (a) g ( x S , t ) ˆ g ( x S ,t )
50 100 150 200 25050100150200250 (b) ˆ g ( x S , t ) ˜ g ( x S ,t )
50 100 150 200 25050100150200250 (c) ˜ g ( x S , t ) ˆ g ( x S ,t ) − g ( x S ,t )
50 100 150 200 25050100150200250 −0.1−0.0500.050.1 (d) Error of ˆ g ( x S , t ) ˜ g ( x S ,t ) − g ( x S ,t )
50 100 150 200 25050100150200250 −0.06−0.04−0.0200.020.04 (e) Error of ˜ g ( x S , t ) Fig. 8. Clock phantom. PAT data at time step t , (a) simulated full data g ( x S , t ) , (b) compressed ˆ g ( x S , t ) and (c) reconstructed ˜ g ( x S , t ) data withstandard Curvelet transform, C , . (d) Compression error ˆ g ( x S , t ) − g ( x S , t ) , (e) recovery error ˜ g ( x S , t ) − g ( x S , t ) . τ = 0 .
05 max(ΨΦ T b t ) , µ = 0 .
75 max(ΨΦ T b t ) / (cid:107) b t (cid:107) andthe stopping tolerance · − or maximum iterations) isillustrated in Figure 14(c).VII. C ONCLUSIONS AND DISCUSSION
We presented a method for acoustic field reconstructionfrom a limited number of patterned measurements from anoptical ultrasound detector. When compressively sensing pho-toacoustic signals, our method recovered the pressure field on g ( x S ,t )
50 100 150 200 25050100150200250 −0.3−0.2−0.100.10.20.3 (a) g ( x S , t ) ˆ g ( x S ,t )
50 100 150 200 25050100150200250 (b) ˆ g ( x S , t ) ˜ g ( x S ,t )
50 100 150 200 25050100150200250 (c) ˜ g ( x S , t ) ˆ g ( x S ,t ) − g ( x S ,t )
50 100 150 200 25050100150200250 −0.1−0.08−0.06−0.04−0.0200.020.040.060.08 (d) Error of ˆ g ( x S , t ) ˜ g ( x S ,t ) − g ( x S ,t )
50 100 150 200 25050100150200250 −0.06−0.04−0.0200.020.040.06 (e) Error of ˜ g ( x S , t ) Fig. 9. Clock phantom. PAT data at time step t , (a) simulated fulldata g ( x S , t ) , (b) compressed ˆ g ( x S , t ) and (c) reconstructed ˜ g ( x S , t ) datawith low-frequency Curvelet transform, C , , , . (d) Compression error ˆ g ( x S , t ) − g ( x S , t ) , (e) recovery error ˜ g ( x S , t ) − g ( x S , t ) . ˜ p : full data
50 100 150 200 25050100150200250 00.10.20.30.40.50.60.70.80.9 (a)
T R ( g ( x S , t )) ˜ p : 18% data, C ,
50 100 150 200 25050100150200250 (b)
T R (˜ g ( x S , t )) ˜ p : 18% data, C , , ,
50 100 150 200 25050100150200250 (c)
T R (˜ g LF ( x S , t )) Error ˜ p : 18% data, C ,
50 100 150 200 25050100150200250 −0.2−0.15−0.1−0.0500.050.1 (d) Error (b)-(a)
Error ˜ p : 18% data, C , , ,
50 100 150 200 25050100150200250 (e) Error (c)-(a)Fig. 10. Clock phantom. Central slice through reconstructed PAT image ˜ p from (a) full data g ( x S , t ) , (b) data reconstructed with standard Curvelettransform C , , (c) data reconstructed with low-frequency Curvelettransform C , , , . The corresponding error (d),(e). the sensor (the PAT data) at each time step independently usingthe sparsity of the data in a low-frequency Curvelet frame,which is a modification of the standard Curvelets tailored toaccount for the smoothing of the wave front during opticalacquisition. The major advantage of such a scheme is that theseries of problems to solve are standard (2D) CS recoveryproblems and as they are independent they can be solvedin parallel. Furthermore, decoupling the CS and acousticreconstruction affords more flexibility in PAT modeling, forinstance including absorption or nonlinear effects, and allowsfor the use of highly optimized readily available software forphotoacoustic image reconstruction.One of the major challenges for compressed sensing ofphotoacoustic signals is that the same number of interrogationpatterns is applied in each time step. As the sparsity of thewave field on the detector varies, in the proposed scheme thisinevitably results in different quality reconstruction in different OURNAL OF L A TEX CLASS FILES, VOL. 11, NO. 4, DECEMBER 2012 10 (a) (b)Fig. 11. (a) Knotted ink tubes phantom on the FP sensor. (b) Photo of theartificial hair phantom. −4 −3 −2 t MSE(ˆ g ( x S ,t ))MSE(˜ g ( x S ,t )) Fig. 12. Knotted tubes. MSE of the compressed PAT data ˆ g ( x S , t ) versusthe reconstructed PAT data ˜ g ( x S , t ) over time for low-frequency Curvelettransform C , , , . time steps. In particular, we first lose the “sharpness” of thewave front, because this is reflected in the coefficients at thehighest scale which magnitude is generally smaller. This effectis partially counterweight by the proposed low-frequencyCurvelet representation which is tailored to frequency rangeof photoacoustic data and hence boosts those coefficients.Furthermore, in our experiments we observed that the datarecovery errors were partially alleviated during the acousticinversion. As the initial pressure p is mapped to the entiretime series, the entire time series carries the information aboutthe wave front and the acoustic inversion acts to average outdata errors. In [12], [13] this problem was tackled utilizing thesparsity directly in the photoacoustic image p , rather than inthe data. There we solve one large (3D) CS recovery problemwhere the CS sensing operator is a composition of the acousticpropagation and pattern measurements. Such forward operatoris significantly more expensive to apply and its incoherenceproperties have to be analyzed. Finally, nonlinear effects suchas acoustic absorption are out of scope of the standard linearCS framework.In future work, we intend to extend the here proposed acous-tic field reconstruction method to efficiently handle dynamicproblems. A CKNOWLEDGMENT
This research was supported by Engineering and PhysicalSciences Research Council, UK (EP/K009745/1).R
EFERENCES[1] D. L. Donoho, “Compressed sensing,”
IEEE Trans. Inf. Theory , vol. 52,no. 4, pp. 1289–1306, Apr. 2006.[2] E. J. Cand´es, J. Romberg, and T. Tao, “Robust uncertainty principles:exact signal reconstruction from highly incomplete frequency informa-tion,”
IEEE Trans. Inf. Theory , vol. 52, no. 2, pp. 489–509, Feb. 2006. ˜ p : full data
20 40 60 80 100 120 14020406080100120140 0246810 (a)
T R ( g ( x S , t )) ˜ p lin ∗ .
77: 25% data
20 40 60 80 100 120 14020406080100120140 (b)
T R (˜ g lin ( x S , t )) ˜ p : 25% data
20 40 60 80 100 120 14020406080100120140 (c)
T R (˜ g ( x S , t )) Di ff erence
20 40 60 80 100 120 14020406080100120140 0.050.10.150.20.250.30.350.4 (d) Error (b) - (a) Di ff erence
20 40 60 80 100 120 14020406080100120140 0.10.20.30.40.50.60.70.80.9 (e) Error (c) - (a)Fig. 13. Knotted tubes. Maximum intensity projection of reconstructedPAT image ˜ p from (a) full data g ( x S , t ) , (b) linearly reconstructed data ˜ g lin ( x S , t ) (here plotted scaled by a factor . ), (c) nonlinearly recon-structed data ˜ g ( x S , t ) with low-frequency Curvelet transform C , , , , and(d),(e) the corresponding error. ˜ p : full data
20 40 60 80 100 12020406080100120 12345678910x 10 −5 (a) T R ( g ( x S , t )) ˜ p lin ∗ .
12: 18% data
20 40 60 80 100 12020406080100120 (b)
T R (˜ g lin ( x S , t )) ˜ p : 18% data
20 40 60 80 100 12020406080100120 (c)
T R (˜ g ( x S , t )) Fig. 14. Hair knot. Maximum intensity projection of reconstructed PAT image ˜ p from (a) full data g ( x S , t ) , (b) linearly reconstructed data ˜ g lin ( x S , t ) (here plotted scaled by a factor . ), (c) nonlinearly reconstructed data ˜ g ( x S , t ) with low-frequency Curvelet transform C , , , .[3] E. J. Cand´es, J. K. Romberg, and T. Tao, “Stable signal recovery fromincomplete and inaccurate measurements,” Commun. Pure Appl. Math. ,vol. 59, no. 8, pp. 1207–1223, 2006.[4] M. F. Duarte, M. A. Davenport, D. Takbar, J. N. Laska, T. Sun, K. F.Kelly, and R. G. Baraniuk, “Single-Pixel Imaging via CompressiveSampling,”
IEEE Signal Processing Magazine , vol. 25, no. 2, pp. 83–91,Mar. 2008.[5] E. J. Candes and L. Demanet, “The curvelet representation of wavepropagators is optimally sparse,”
Commun Pur Appl Math , vol. 58,no. 11, pp. 1472–1528, 2005.[6] E. J. Cand´es and D. L. Donoho, “New tight frames of curvelets andoptimal representations of objects with piecewise C singularities,” Commun. Pure Appl. Math , vol. 57, no. 2, pp. 219–266, Feb. 2004.[7] N. Huynh, E. Zhang, M. Betcke, S. Arridge, P. Beard, and B. Cox,“Single-pixel optical camera for video rate ultrasonic imaging,”
Optica ,vol. 3, no. 1, pp. 26–29, Jan. 2016.[8] J. Provost and F. Lesage, “The application of compressed sensing forphoto-acoustic tomography,”
IEEE Transactions on Medical Imaging ,vol. 28, no. 4, pp. 585–594, Apr. 2009.[9] Z. Guo, C. Li, L. Song, and L. V. Wang, “Compressed sensing inphotoacoustic tomography in vivo,”
J. Biomed. Opt. , vol. 15, no. 2, p.021311, 2010.[10] C. Huang, K. Wang, L. Nie, L. V. Wang, and M. A. Anastasio, “Full-Wave Iterative Image Reconstruction in Photoacoustic Tomography WithAcoustically Inhomogeneous Media,”
IEEE Transactions on MedicalImaging , vol. 32, no. 6, pp. 1097–1110, Jun 2013.[11] Z. Belhachmi, T. Glatz, and O. Scherzer, “A direct method for photoa-coustic tomography with inhomogeneous sound speed,”
Inverse Prob-lems , vol. 32, no. 4, p. 045005.[12] S. Arridge, M. Betcke, B. Cox, F. Lucka, and B. Treeby, “On the AdjointOperator in Photoacoustic Tomography,”
ArXiv e-prints , Feb. 2016.
OURNAL OF L A TEX CLASS FILES, VOL. 11, NO. 4, DECEMBER 2012 11 [13] S. Arridge, P. Beard, M. Betcke, B. Cox, N. Huynh, F. Lucka, O. Ogun-lade, and E. Zhang, “Accelerated High-Resolution Photoacoustic To-mography via Compressed Sensing,”
ArXiv e-prints , Apr. 2016.[14] M. Haltmeier, T. Berer, S. Moon, and P. Burgholzer, “Compressedsensing and sparsity in photoacoustic tomography,”
ArXiv e-prints , May2016.[15] L. V. Wang,
Photoacoustic Imaging and Spectroscopy . CRC Press,Mar. 2009, vol. 144.[16] P. Beard, “Biomedical photoacoustic imaging.”
Interface Focus , vol. 1,no. 4, pp. 602–631, Aug. 2011.[17] K. S. Valluru, K. E. Wilson, and J. K. Willmann, “Photoacoustic imagingin oncology: Translational preclinical and early clinical experience,”
Radiology , vol. 280, no. 2, pp. 332–349, 2016.[18] Y. Zhou, J. Yao, and L. V. Wang, “Tutorial on photoacoustic tomogra-phy,”
Journal of Biomedical Optics , vol. 21, no. 6, p. 061007, 2016.[19] J. Xia and L. V. Wang, “Small-animal whole-body photoacoustic to-mography: A review,”
IEEE Transactions on Biomedical Engineering ,vol. 61, no. 5, pp. 1380–1389, May 2014.[20] E. Zhang, J. Laufer, and P. Beard, “Backward-mode multiwavelengthphotoacoustic scanner using a planar Fabry-Perot polymer film ultra-sound sensor for high-resolution three-dimensional imaging of biologicaltissues,”
Appl. Opt. , vol. 47, no. 4, pp. 561–577, Feb 2008.[21] B. E. Treeby, J. Jaros, A. P. Rendell, and B. T. Cox, “Modelingnonlinear ultrasound propagation in heterogeneous media with powerlaw absorption using a k-space pseudospectral method.”
The Journal ofthe Acoustical Society of America , vol. 131, no. 6, pp. 4324–36, 2012.[22] D. Finch and R. Sarah K. Patch, “Determining a function from itsmean values over a family of spheres,”
SIAM Journal on MathematicalAnalysis , vol. 35, no. 5, pp. 1213–1240, 2004.[23] P. Kuchment and L. Kunyansky,
Mathematics of Photoacoustic andThermoacoustic Tomography . New York, NY: Springer New York,2011, pp. 817–865.[24] P. Stefanov and G. Uhlmann, “Thermoacoustic tomography with variablesound speed,”
Inverse Problems , vol. 25, no. 7, p. 075011.[25] J. Frikel and E. T. Quinto, “Artifacts in incomplete data tomography withapplications to photoacoustic tomography and sonar,”
SIAM Journal onApplied Mathematics , vol. 75, no. 2, pp. 703–725, 2015.[26] B. E. Treeby and B. T. Cox, “k-Wave: Matlab toolbox for the simulationand reconstruction of photoacoustic wave fields,”
J Biomed Opt , vol. 15,no. 2, p. 021314, 2010.[27] E. J. Cand´es, “The restricted isometry property and its implications forcompressed sensing,”
Comptes Rendus Mathematique , vol. 346, pp. 589–592, 2008.[28] M. Fornasier and H. Rauhut, “Compressive Sensing,” in
Handbook ofMathematical Methods in Imaging . New York, NY: Springer New York,2011, pp. 187–228.[29] B. T. Cox and P. C. Beard, “Fast calculation of pulsed photoacousticfields in fluids using k-space methods,”
The Journal of the AcousticalSociety of America , vol. 117, no. 6, pp. 3616–3627, 2005.[30] K. P. K¨ostli, M. Frenz, H. Bebie, and H. P. Weber, “Temporal backwardprojection of optoacoustic pressure transients using Fourier transformmethods,”
Phys Med Biol , vol. 46, no. 7, pp. 1863–1872, Jul. 2001.[31] Y. Xu, D. Feng, and L. V. Wang, “Exact frequency-domain recon-struction for thermoacoustic tomography. I. Planar geometry,”
IEEETransactions on Medical Imaging , vol. 21, no. 7, pp. 823–828, July2002.[32] M. Tabei, T. D. Mast, and R. C. Waag, “A k-space method for coupledfirst-order acoustic propagation equations,”
J Acoust Soc Am , vol. 111,no. 1, pp. 53–63, Jan. 2002.[33] B. E. Treeby and B. T. Cox, “A k-space Green’s function solution foracoustic initial value problems in homogeneous media with power lawabsorption,”
Journal of the Acoustical Society of America , 2011.[34] E. Cand´es, L. Demanet, D. Donoho, and L. Ying, “Fast discrete curvelettransforms,”
Multiscale Model. Simul. , vol. 5, no. 3, pp. 861–899, 2006.[35] E. J. Cand`es and D. L. Donoho, “New tight frames of curvelets andoptimal representations of objects with piecewise C-2 singularities,”
Commun Pur Appl Math , vol. 57, no. 2, pp. 219–266, Feb. 2004.[36] D. L. Donoho, “Sparse components of images and optimal atomicdecompositions,”
Constr. Approx. , vol. 17, no. 3, pp. 353–382, 2001.[37] B. S. Kashin, “Approximation properties of complete orthonormalsystems,”
Trudy Mat. Inst. Steklov. , vol. 172, pp. 187–191, 353, 1985.[38] N. Huynh, E. Zhang, M. Betcke, S. Arridge, P. Beard, and B. Cox,“Patterned interrogation scheme for compressed sensing photoacousticimaging using a Fabry Perot planar sensor,” in
SPIE BiOS , A. A.Oraevsky and L. V. Wang, Eds. SPIE, Mar. 2014, p. 894327. [39] ——, “A real-time ultrasonic field mapping system using a Fabry Perotsingle pixel camera for 3D photoacoustic imaging,” vol. 9323, 2015, pp.93 231O–93 231O–7.[40] S. Foucart and H. Rauhut,
A Mathematical Introduction to CompressedSensing . Springer New York Heidelberg Dordrecht London, 2010.[41] M. V. Afonso, J. M. Bioucas-Dias, and M. A. T. Figueiredo, “Fast imagerecovery using variable splitting and constrained optimization,”