An axially-variant kernel imaging model applied to ultrasound image reconstruction

1 An axially-v ariant kernel imaging model applied to ultrasound image reconstruction Mihai I. Florea, Adrian Basarab, Denis K ouam ´ e, and Sergiy A. V orobyov Abstract —Existing ultrasound decon volution approaches un- realistically assume, primarily f or computational reasons, that the con volution model relies on a spatially in variant kernel and circulant boundary conditions. W e discard both r estrictions and introduce an image f ormation model applicable to ultrasound imaging and decon volution based on an axially varying kernel, that accounts for arbitrary boundary conditions. Our model has the same computational complexity as the one employing spatially in variant con volution and has negligible memory requirements. T o accommodate state-of-the-art decon volution approaches when applied to a variety of in verse problem f ormulations, we also pro vide an equally efficient adjoint expression of our model. Simulation results confirm the tractability of our model f or the deconv olution of large images. Moreover , the quality of reconstruction using our model is superior to that obtained using spatially in variant con volution. Index T erms —Axially varying, deconv olution, f orward model, kernel, matrix-free, point-spr ead function, ultrasound I . I N T RO D U C T I O N Ultrasound imaging is a medical imaging modality widely adopted due to its ef ficiency , lo w cost, and safety . These adv an- tages come at the expense of image quality . Consequently , the accurate estimation of the tissue reflecti vity function (TRF) from ultrasound images is the subject of activ e research. Generally , existing approaches assume that the formation of ultrasound images follows a 2D conv olution model between the TRF and the system kernel. The con volution model is further constrained for computational reasons to hav e spa- tially in variant kernel and circulant boundary conditions (see e.g. [1]–[6]). Pulse-echo emission of focused wav es still remains the most widely used acquisition scheme in ultrasound imaging. It con- sists of sequentially transmitting narro w focused beams. For each transmission centered at a lateral position, the raw data is used to beamform one RF signal. Gi ven the repeatability of the process in the lateral direction, the k ernels do not v ary laterally . Howe ver , despite dynamic focusing in reception and time gain compensation, the kernels become wider as we move away from the focal depth, thus degrading the spatial resolution and motiv ating the proposed kernel variation model. Previous works account for this variation by assuming local regions of kernel in variance and performing decon volution M. I. Florea and S. A. V orobyo v are with Aalto Univ ersity , Dept. Signal Processing and Acoustics, FI-00076, AAL TO, Finland (e-mail: mihai.florea@aalto.fi; sergiy .vorobyov@aalto.fi). A. Basarab and D. Kouam ´ e are with IRIT , CNRS UMR 5505, Univ ersity of T oulouse (e-mail: basarab@irit.fr; kouame@irit.fr). Part of this work has been supported by the Academy of Finland, under grant no. 299243, and CIMI Labex, T oulouse, France, under grant ANR-11-LABX-0040-CIMI within the program ANR-11-IDEX-0002-02. block-wise (e.g. [7]). V ery recently , ultrasound imaging con- volution models with continuously varying kernels were pro- posed in [8] and [9]. Howe ver , [8] makes the ov erly restrictive assumption that the spatially varying kernel is obtained from a constant reference kernel modulated by the exponential of a fixed discrete generator scaled by the varying kernel center image coordinates. Therefore, it does not take into account the depth-dependent spatial resolution degradation explained previously . On the other hand, the deconv olution proposed in [9] has an iteration complexity proportional to the cube of the number of pixels in the image, limiting its applicability to very small images. The contributions of this work are as follows. (i) W e propose a nov el axially-variant kernel ultrasound image for- mation model (Section III). (ii) Our model is linear and may be implemented as a matrix. Howe ver , the matrix form does not scale because its complexity is proportional to the square of the number of pixels in the image. Therefore we provide an efficient matrix-free implementation of axially varying conv olution that entails the same computational cost as spatially in variant con volution (Subsection III-B). (iii) The decon volution problem is ill-posed and many deconv olution models can only be solved approximately using proximal splitting methods (see [10], [11], and references therein) that compute at e very iteration the gradient of a data fidelity term. The data fidelity gradient e xpression includes calls to both the model operator and its adjoint. W e express this adjoint operator in a form of equal complexity to that of the forward model operator (Section IV). (iv) W e confirm using simulation results that decon volution with our model is tractable e ven for large images and produces results superior to those obtained using the spatially in v ariant model (Section V). I I . N OTA T I O N In this work, images (ultrasound images and TRFs) are vectorized in column-major order but referenced in 2D form. For instance, image v ∈ R m v n v corresponds to a m v × n v 2D image and has the pixel v alue at coordinates ( i, j ) given by v m v ( j − 1)+ i . Howe ver , for clarity of exposition, we denote it as 2D object v ∈ R m v × n v , with the pixel value at location ( i, j ) giv en by v i,j . Bold marks this artificial indexing. Similarly , linear operators are matrices b ut referred to as 4D tensors, e.g., O : R m v × n v → R m w × n w denotes O ∈ R m v n v × m w n w . In the sequel, we define sev eral classes of linear operators that constitute the mathematical building blocks of our pro- posed model and analysis. Note that these are more general than normal linear operators because their dimensions not 2 k ∗ 1 a k ∗ 2 a k a (a) a b = W a ( i 1 , i 2 ) a c = Z a ( i 1 , i 2 ) b m a n a i 1 i 2 m a n a (b) Fig. 1. (a) Con v olving test image a with a Gaussian k ernel k . The inner rectangle represents v alid con v olution whereas the outer marks full con v olution; (b) Applying the full-width wi ndo w operator , follo wed by a full-width zero-padding operator on a test image a . Here, black and white correspond to v alues of 1 and 0 , respecti v ely . K ernel k is displayed after min-max normalization. only depend on those of their parameters b ut also of their ar guments. A. Con volution oper ator s F or all m k , n k ≥ 1 , all k ernels k ∈ R m k × n k , and all m a ≥ m k , n a ≥ n k , we define the linear operators C 1 ( k ) and C 2 ( k ) as C 1 ( k ) a def = k ∗ 1 a , C 2 ( k ) a def = k ∗ 2 a , for all a ∈ R m a × n a , where operations ∗ 1 and ∗ 2 denote (discrete) v alid con v olution and full con v olution, respecti v ely , defined as ( k ∗ 1 a ) i,j def = m k X p =1 n k X q =1 k p,q a i − p + m k ,j − q + n k , i ∈ { 1 , ..., m a − m k + 1 } , j ∈ { 1 , ..., n a − n k + 1 } , ( k ∗ 2 a ) i,j def = ¯ p i X p = p i ¯ q j X q = q j k p,q a i − p +1 ,j − q +1 , i ∈ { 1 , ..., m a + n k − 1 } , j ∈ { 1 , ..., n a + n k − 1 } , p i = max { 1 , i − m a + 1 } , ¯ p i = min { i, m k } , q j = max { 1 , j − n a + 1 } , ¯ q j = min { j , n k } . The dif ference between the tw o forms of con v olution is e x em- plified in Fig. 1(a). V alid con v olution is thereby the subset of full con v olution where e v ery output pix el is e xpressed using the entire k ernel k . B. A uxiliary oper ator s F or conciseness, we also introduce the follo wing auxiliary operators. None in v olv e an y computation in practice. Let the rotation operator R ( k ) be gi v en by ( R ( k )) i,j def = k m k − i +1 ,n k − j +1 , i ∈ { 1 , ..., m k } , j ∈ { 1 , ..., n k } . T o further simplify notation, for all indices 1 ≤ a ≤ b ≤ c , we denote the e xception inde x set I ( a, b, c ) as I ( a, b, c ) def = { 1 , ..., c } \ { a, ..., b } . The full-width windo w and zero padding operators are defined as ( W s ( i 1 , i 2 ) a ) i,j def = a i + i 1 ,j , i ∈ { 0 , ..., i 2 − i 1 } , ( Z s ( i 1 , i 2 ) a ) i,j def =  a i − i 1 ,j , i ∈ { i 1 , ..., i 2 } , 0 , i ∈ I ( i 1 , i 2 , m s ) , where j ∈ { 1 , ..., n s } and inde x s ∈ { t, p } stands for image size quantities m t , m p , n t , and n p . Their ef fect on a test image is sho wn in Fig. 1(b). I I I . A X I A L L Y - V A R I A N T K E R N E L B A S E D U L T R A S O U N D I M A G I N G M O D E L W e propose the follo wing image formation model y = H P x + n , (1) where x , y , n ∈ R m t × n t denote the TRF to be reco v ered, the observ ed radio-frequenc y (RF) image, and independent identically distrib uted (i.i.d.) additi v e white Gaussian noise (A WGN), respecti v ely . A. P adding Operator P : R m t × n t → R m p × n p pads the TRF with a boundary of width n r and height m r , yielding an image of size m p = m t + 2 m r times n p = n t + 2 n r . P adding in our ultrasound imaging model allo ws us to reconstruct a TRF of the same size as the observ ed RF image. T o this end, we must simulate the ef fect the surrounding tissues ha v e on the imaged tissues. P adding is an estimation of the surrounding tissues using information from the imaged TRF . This estimation only af fects the border of the reconstructed TRF . If this border information is not required, the reconstructed TRF can simply be cropped accordingly . The addition of padding to our model brings the adv antage of accommodating both options. F or computational reasons, P is assumed linear and separa- ble along the dimensions of the image. Separability translates to P = P m P n . Here, P m pads e v ery column of the image independently by applying the 1D padding (linear) operator P ( m t , m r ) . Consequently , when n t = 1 and n r = 0 , oper - ators P and P ( m t , m r ) are equi v alent. The ro w component P n treats e v ery ro w as a column v ector , applies P ( n t , n r ) to it, and turns the result back into a ro w . P adding, either in 1D or 2D, can be performed without e xplicitly deri ving an operator matrix. Ho we v er , the matrix form f acilitates the formulation of the corresponding adjoint operator . Common matrix forms of operator P ( m t , n t ) are sho wn in Fig. 2 for m t = 1 0 and m r = 3 . These e xamples demonstrate that the matrix form of P ( m t , m r ) can be easily generated programatically and, due to its sparsity , can be stored in memory e v en for v ery lar ge v alues of m t and m r . These properties e xtend to the matrix form of the 2D padding operator P by virtue of the follo wing result. Theor em 1. P adding oper ator P can be obtained pr o gr amat- ically in the form of a spar se matrix as P = P ( n t , n r ) ⊗ P ( m t , m r ) , wher e ⊗ denotes the Kr onec k er pr oduct. Pr oof: See Appendix A. 3 Zero Circular Replicate Symmetric Fig. 2. Common matrix forms of 1D padding operators P (10 , 3) . Black denotes a value of 1 and white denotes 0 . B. Axially varying con volution Linear operator H : R m p × n p → R m t × n t performs the axially-v ariant kernel conv olution. W e define axially- variant conv olution as the linear operation whereby each r ow i h ∈ { 1 , ..., m t } of the output image is obtained by the valid con volution between the kernel pertaining to that row k ( i h ) ∈ R m k × n k , where m k = 2 m r + 1 and n k = 2 n r + 1 , and the corresponding patch in the input (padded) TRF . The auxiliary operators defined in Section II enable us to write H as a sum of linear operators based on the observation that the concatenation of output rows has the same effect as the summation of the rows appropriately padded with zeros. Analytically , this translates to H = m t X i h =1 Z t ( i h , i h ) C 1 ( k ( i h )) W p ( i h , i h + 2 m r ) . (2) In matrix form, operator H would need to store m p n p m t n t coefficients and its inv ocation would entail an equal number of multiplications. Its complexity would thus be greater than the square of the number of pixels in the image, limiting its applicability to medium sized images. Using the matrix-free expression in (2), operator H performs m k n k m t n t multipli- cations and has negligible memory requirements. Therefore, in ultrasound imaging, the matrix-free representation is not only vastly superior to its matrix counterpart (because the kernel is much smaller than the image), but has the same computational complexity as spatially in v ariant con volution (excluding the unrealistic circulant boundary case). Unlike the forward model which, by utilizing operators H and P , can be computed exactly with great ef ficiency , many deconv olution models can only be solved approximately using proximal splitting methods that optimize an objectiv e containing a data fidelity term φ ( H P x − y ) . These methods employ at ev ery iteration the gradient of the data fidelity term, giv en by ∇ ( φ ( H P x − y )) = P T H T ( ∇ φ )( H P x − y ) . (3) Note that, under our A WGN assumption, φ is the square of the ` 2 -norm b ut the results in this work may be applied to other additiv e noise models. The gradient expression in (3) depends on H and P as well as their adjoints. In the follo wing, we deri ve computationally efficient e xpressions for adjoint operators H T and P T . I V . A D J O I N T O F M O D E L O P E R A T O R By taking the adjoint in (2), we get H T = m t X i h =1 ( W p ( i h , i h + 2 m r )) T ( C 1 ( k ( i h ))) T ( Z t ( i h , i h )) T . T o obtain a matrix-free representation of H T , we need the corresponding matrix-free expression for the adjoints of the con volution and auxiliary operators. First, it trivially holds that the windo w operator and corresponding zero padding operator are mutually adjoint, expressed as ( W s ( i 1 , i 2 )) T = Z s ( i 1 , i 2 ) . (4) The adjoint of valid con volution can be linked to full con vo- lution as follows. Theorem 2. The adjoint of valid convolution is full corr elation (con volution with the r otated kernel), namely ( C 1 ( k )) T = C 2 ( R ( k )) . Pr oof: See Appendix B. Theorem 2 and (4) yield a matrix-free expression for H T in the form of H T = m t X i h =1 Z p ( i h , i h + 2 m r ) C 2 ( R ( k ( i h ))) W t ( i h , i h ) . (5) Therefore, operators H and H T hav e equal computational complexity . Moreov er , they exhibit two levels of parallelism. The con v olution operators themselves are fully parallel and the computations pertaining to each ro w i h can be performed concurrently . Thus, in matrix-free from, both operators benefit from parallelization in the same way as their matrix counter- parts. The adjoint of the padding operator, P T , can be obtained either directly through sparse matrix transposition or by ap- plying transposition in Theorem 1 as P T = ( P ) T = ( P ( n t , n r )) T ⊗ ( P ( m t , m r )) T . (6) Finally , note that whereas the column-major order assump- tion can be made without loss of generality for operator H , it is not the case for the padding operator P . In particular , row-major order rev erses the terms in the Kronecker product. V . E X P E R I M E N TA L R E S U LT S W e ha ve tested our model on a simulated ultrasound image decon volution problem. Three TRFs (TRF 1, TRF 2, and TRF 3) were generated following the procedure used in com- puting the kidney phantom within the Field II simulator [12], [13]. Each TRF is an interpolation of Gaussian distrib uted random scatterers with standard deviations determined by a pixel intensity map. The map for TRF 1 is a patch from a human kidney tissue optical scan included in the Field II distribution. The maps for TRF 2 and TRF 3 are patches from a single slice (visual identifier 3272 ) of the female dataset provided by the V isible Human Project [14]. Decon volution experiments were performed for each of the 3 TRFs, with the padding size, matching the kernel radii, giv en by m r × n r . The 4 Fig. 3. (a) Ground truth (in B-mode) of the tissue reflectivity function (TRF); (b) Demodulated kernels k ( i h ) for twenty depths at regularly spaced intervals of 2 mm; (c) Observed B-mode image simulated following the proposed axially-variant convolution model; (d) Axially-in variant decon volution result AI (in B-mode) obtained with a fix ed kernel equal to k ( m t / 2) (the center kernel of the axially-variant model); (e) Axially-v ariant deconv olution result A V (in B-mode) using our model. All images are displayed using a dynamic range of 40 dB. Row (a1)-(e1) pertains to TRF 1, ro w (a2)-(e2) to TRF 2, and row (a3)-(e3) to TRF 3. T ABLE I D E CO N VO L U TI O N E X P ER I M EN T PA R AM E T E RS F O R E AC H T R F TRF m t n t Axial size Lateral size m r n r TRF 1 2480 480 94 mm 95 mm 9 50 TRF 2 2598 480 100 mm 95 mm 7 35 TRF 3 2598 480 100 mm 95 mm 5 25 TRFs are shown in Fig. 3(a1), (a2), and (a3), respectively , and their parameters are listed in T able I. For every row i h ∈ { 1 , ..., m t } , we have defined the k ernel k ( i h ) in (2) as k ( i h ) i,j = ρ µ z ,σ z ( i ) ρ µ x ,σ x ( i h ) ( j ) cos(2 π f 0 /f s ( i − µ z )) , where ρ µ,σ ( x ) is a normalized Gaussian window , given by ρ µ,σ ( x ) = 1 √ 2 π σ exp  − ( x − µ ) 2 2 σ 2  , and parameters µ z and µ x are the center coordinates of the kernel. Axial stan- dard deviation (SD) w as set to σ z = σ 1 and lateral SD to σ x ( i h ) = p ((2 i h ) /m t − 1) 2 ( σ 2 2 − σ 2 1 ) + σ 2 1 , with σ 1 = m r / 3 and σ 2 = n r / 3 . Here, f 0 = 3 MHz and f s = 20 MHz are the ultrasound central and sampling frequencies, respectiv ely . For each TRF decon volution experiment, the depth-dependent width variation of the kernel simulates the lateral spatial resolution degradation when moving away from the focus point, located in the center of the image ( 47 mm from the probe for TRF 1, 50 mm for TRF 2 and 3). The en velopes of these kernels at regular intervals across the image are sho wn in Fig. 3(b1), (b2), and (b3). Whereas the TRFs are of approximately the same size, the intensity of the blur differs in each experiment in order to display both the high and low frequency reconstruction capabilities of our method. W e chose symmetric padding, as illustrated in Fig. 2, because it is more realistic than circular and zero padding and, by using a larger number of pixels from the TRF , more robust to noise than replicate padding. A small amount of noise was added to each TRF such that the signal to noise ratio is 40 dB. The ultrasound 5 images produced from each TRF using our forw ard model in (1) are shown in Fig. 3(c1), (c2), and (c3). T o estimate the TRFs, we have considered an elastic net [15] regularized least-squares (based on the A WGN assumption) decon volution model min x 1 2 k H P x − y k 2 2 + λ 1 k x k 1 + λ 2 2 k x k 2 2 . with manually tuned parameters λ 1 = 2e − 3 and λ 2 = 1e − 4 . For deconv olution, we have employed the Accelerated Com- posite Gradient Method (A CGM) [16], [17] on account of its low resource usage, applicability , adaptability , and near- optimal linear conv ergence rate on elastic net regularized optimization problems. As the term proximal splitting suggests, A CGM splits the objectiv e into two functions f ( x ) = 1 2 k H P x − y k 2 2 , Ψ( x ) = λ 1 k x k 1 + λ 2 2 k x k 2 2 . Every iteration of A CGM is dominated by the computationally intensiv e data fidelity gradient function ∇ f ( x ) = P T H T ( H P x − y ) . All other calculations performed by A CGM are either negligible when compared to ∇ f ( x ) or can be reduced to subexpressions of ∇ f ( x ) . Due to the efficient matrix-free expressions of H in (2) and H T in (5) as well as the sparse matrix implementation of P and P T (easily precomputed using Theorem 1 and (6), respectiv ely), deconv olution with our model entails the same computational cost as with a fixed kernel model. The results of axially-in variant decon volution (AI) are shown in Fig. 3(d1), (d2), and (d3) and using our axially- variant model (A V) in Fig. 3(e1), (e2), and (e3), all after 150 iterations. Our approach achie ves almost perfect lo w frequency reconstruction for each TRF . For higher frequencies, blurring destroys visual information. Therefore, reconstruction quality at this lev el depends on the size of the kernels. F or TRF 3, axial standard de viation v alues are the smallest and the blurring process has not suppressed all minor details. Here, our method manages to reconstruct ev en small, high contrast features. For all TRFs, the gain in reconstruction quality is e vident especially in the upper and lower extremities, as can be discerned from Fig. 3. Interestingly , ev en though the two models dif fer only slightly in the center of the image, our model performs better in that region as well. V I . C O N C L U S I O N S In this work, we have proposed an axially varying con- volution forward model for ultrasound imaging. The physics of ultrasound image formation as well as our decon v olution simulation results sho w the superiority of our model over the traditional fixed kernel model. Our matrix-free formulas for the adjoints of the con volution and auxiliary operators, necessary for the implementation of decon v olution using proximal splitting techniques, also constitute a solid theoretical foundation for decon volution methodologies using more sophisticated models, particularly those where the kernel also varies along the lateral direction. Furthermore, our theoretical results and methodology are not restricted to ultrasound imaging and may be extrapolated to other imaging modalities and applications. A P P E N D I X A P RO O F O F T H E O R E M 1 Since image columns are staked one on top of the other, P m will ha ve a block diagonal structure, with each block gi ven by P ( m t , m r ) , namely P m = diag { P ( m t , m r ) , P ( m t , m r ) , ..., P ( m t , m r ) } = I n t ⊗ P ( m t , m r ) . (7) Matrix P n lacks this block structure. It instead has the ele- ments of P ( n t , n r ) strided horizontally by m t and replicated along the diagonal. This can be expressed succinctly as P m =    P ( n t , n r ) 1 , 1 I m t · · · P ( n t , n r ) 1 ,n t I m t . . . . . . . . . P ( n t , n r ) m p , 1 I m t · · · P ( n t , n r ) m p ,n t I m t    = P ( n t , n r ) ⊗ I m t . (8) By combining the results along the dimensions in (7) and (8), we can compute P in closed form P = P m P n = ( I n t ⊗ P ( m t , m r ))( P ( n t , n r ) ⊗ I m t ) = P ( n t , n r ) ⊗ P ( m t , m r ) . (9) A P P E N D I X B T H E A D J O I N T O F V A L I D C O N VO L U T I O N I S F U L L C O R R E L A T I O N A. Additional notation W ithin the context of this section, we define more general window and zero padding operators, respecti vely , as W L,H a def = a i + m L ,j + n L , i ∈ { 0 , ..., m H − m L } , j ∈ { 0 , ..., n H − n L } , Z L,H a def =  a i − m L ,j − n L , i ∈ { m L , ..., m H } , j ∈ { n L , ..., n H } , 0 , i ∈ I ( m L , m H , m N ) , j ∈ I ( n L , n H , n N ) , where indices L, H ∈ { 1 , k , a, M , N } stand for quantities m 1 , m k , m a , etc., with m 1 = n 1 = 1 and m M = m a − m k + 1 , n M = n a − n k + 1 , m N = m a + m k − 1 , n N = n a + n k − 1 . Central to our proof is circular (periodic) con volution [18]. It only applies to arguments of equal size. When k and a are not the same size, they can only be circularly con volved by appropriately padding them with zeros. W e consider ar guments ¯ k and ¯ a of equal size m by n . For brevity , we introduce the circular sum ⊕ and dif ference  , respectiv ely gi ven by a ⊕ c b def = (( a + b − 2) mod c ) + 1 , a  c b def = (( a − b ) mod c ) + 1 , for all c ≥ 1 and a, b ∈ { 1 , ..., c } . 6 Circular con volution ~ and its corresponding linear operator C ( ¯ k ) are defined as ( ¯ k ~ ¯ a ) i,j def = m X p =1 n X q =1 ¯ k p,q ¯ a i  m p,j  n q , i ∈ { 1 , ..., m } , j ∈ { 1 , ..., n } , C ( ¯ k ) ¯ a def = ¯ k ~ ¯ a . W e further define the circular shift operator S ( p, q ) as S ( p, q ) ¯ a def = e ( p, q ) ~ ¯ a , where e ( p, q ) is the standard basis vector image e ( p, q ) i,j =  1 , i = p, j = q , 0 , otherwise. It follows that S (1 , 1) is the identity operator . B. Supporting r esults Lemma 1. Shift operator s cumulate and can be taken out of con volutions as C ( S ( p 1 , q 1 ) ¯ k ) S ( p 2 , q 2 ) = S ( p 1 ⊕ m p 2 , q 1 ⊕ n q 2 ) C ( ¯ k ) . Pr oof: The result follows tri vially from the commutativity and associativity properties of circular con volution C ( S ( p 1 , q 1 ) ¯ k ) S ( p 2 , q 2 ) ¯ a = ( e ( p 1 , q 1 ) ~ ¯ k ) ~ ( e ( p 2 , q 2 ) ~ ¯ a ) = ( e ( p 1 , q 1 ) ~ e ( p 2 , q 2 )) ~ ¯ k ~ ¯ a = e ( p 1 ⊕ m p 2 , q 1 ⊕ n q 2 ) ~ ( ¯ k ~ ¯ a ) = S ( p 1 ⊕ m p 2 , q 1 ⊕ n q 2 ) C ( ¯ k ) ¯ a . Lemma 2. The adjoint of circular con volution is cir cular cr oss-correlation cir cularly shifted forwar d by one position ( C ( ¯ k )) T = S (2 , 2) C ( R ( ¯ k )) . Pr oof: The conv olution theorem [19] states that ¯ k ~ ¯ a = F H (( F ¯ k )  ( F ¯ a )) = F H diag( F ¯ k ) F ¯ a , (10) where F is the Discrete Fourier Transform (DFT), ( . ) H de- notes the Hermitian adjoint,  represents Hadamard (element- wise) product and diag( x ) produces a diagonal matrix (linear operator) with the entries of x . Therefore, the circular con v o- lution operator is diagonalized in Fourier domain as C ( ¯ k ) = F H diag ( F ¯ k ) F . (11) T aking the adjoint in (11), we obtain that ( C ( ¯ k )) T = ( C ( ¯ k )) H = F H diag (( F ¯ k ) ∗ ) F = F H diag ( F ∗ ¯ k ) F , (12) where ( . ) ∗ denotes the complex element-wise conjugate. The DFT matrix has conjugate symmetry [18], namely F ∗ ¯ k = F S (2 , 2) R ( ¯ k ) . (13) Substituting (13) in (12) we have ( C ( ¯ k )) T = F H diag ( F S (2 , 2) R ( ¯ k )) F , which can be expressed using con volution theorem (10) as ( C ( ¯ k )) T = C ( S (2 , 2) R ( ¯ k )) . (14) The application of Lemma 1 with p 1 = q 1 = 2 , p 2 = q 2 = 1 , in (14) completes the proof. C. Pr oof of Theorem 2 All necessary theoretical results pertaining to circular con- volution ha ve already been deriv ed. W e can now re vert to kernel k and image a . Their dimensions are as specified in Section II. Using the more general window and zero padding operators, we can express the valid and full conv olution linear operators, respectiv ely , as C 1 ( k ) = W k,a C ( Z 1 ,k k ) Z 1 ,a , (15) C 2 ( k ) = W 1 ,a C ( Z 1 ,k k ) Z 1 ,M . (16) Note that if C 1 ( k ) takes as arguments images of size m a × n a , its adjoint expression in volv es an operator C 2 ( k ) with argu- ments of size m M × n M . W e apply the adjoint in (15) and obtain ( C 1 ( k )) T = ( Z 1 ,a ) T ( C ( Z 1 ,k k )) T ( W k,a ) T . (17) Like their full-width counterparts, the window operator and corresponding zero padding operators introduced in this sec- tion are also mutually adjoint, namely ( W L,H ) T = Z L,H . (18) Applying (18) in (17) we obtain ( C 1 ( k )) T = W 1 ,a ( C ( Z 1 ,k k )) T Z k,a . (19) Lemma 2 giv es ( C ( Z 1 ,k k )) T = S (2 , 2) C ( R ( Z 1 ,k k )) = S (2 , 2) C ( Z a,N R ( k )) = S (2 , 2) C ( S ( m a , n a ) Z a,N R ( k )) = S ( m a + 1 , n a + 1) C ( Z 1 ,k R ( k )) . (20) Using (20) in (19) and applying Lemma 1 we obtain that ( C 1 ( k )) T = W 1 ,a S ( m a + 1 , n a + 1) C ( Z 1 ,k R ( k )) S ( m k , n k ) Z 1 ,M = W 1 ,a S (1 , 1) C ( Z 1 ,k R ( k )) Z 1 ,M = W 1 ,a C ( Z 1 ,k R ( k )) Z 1 ,M . (21) The desired result follo ws by observing that the right-hand sides of (21) and (16) are identical. 7 R E F E R E N C E S [1] J. Ng, R. Prager, N. Kingsbury , G. Treece, and A. Gee, “W avelet restora- tion of medical pulse-echo ultrasound images in an EM framework, ” IEEE T rans. Ultrason. F err oelectr . Fr eq. Contr ol , vol. 54, no. 3, pp. 550–568, 2007. [2] R. Rangarajan, C. V . Krishnamurthy , and K. Balasubramaniam, “Ultra- sonic imaging using a computed point spread function, ” IEEE T rans. Ultrason. F err oelectr . F req. Contr ol , vol. 55, no. 2, pp. 451–464, Feb. 2008. [3] H.-C. Shin, R. Prager, J. Ng, H. Gomersall, N. Kingsbury , G. T reece, and A. Gee, “Sensiti vity to point-spread function parameters in medical ultrasound image deconv olution, ” Ultrasonics , vol. 49, no. 3, pp. 344 – 357, 2009. [4] M. Alessandrini, S. Maggio, J. Poree, L. D. Marchi, N. Speciale, E. Franceschini, O. Bernard, and O. Basset, “ A restoration framew ork for ultrasonic tissue characterization, ” IEEE T rans. Ultrason. F err oelectr . F req. Contr ol , vol. 58, no. 11, pp. 2344–2360, 2011. [5] C. Dalitz, R. Pohle-Frohlich, and T . Michalk, “Point spread functions and deconvolution of ultrasonic images, ” IEEE T rans. Ultrason. F erro- electr . Fr eq. Contr ol , vol. 62, no. 3, pp. 531–544, Mar . 2015. [6] N. Zhao, A. Basarab, D. K ouam ´ e, and J.-Y . T ourneret, “Joint seg- mentation and deconv olution of ultrasound images using a hierarchical Bayesian model based on generalized Gaussian priors, ” IEEE T rans. Image Pr ocess. , vol. 25, no. 8, pp. 3736–3750, 2016. [7] J. G. Nagy and D. P . O’Leary , “Restoring images degraded by spatially variant blur, ” SIAM J. Sci. Comput. , v ol. 19, no. 4, pp. 1063–1082, Jul. 1998. [8] O. V . Michailovich, “Non-stationary blind decon volution of medical ultrasound scans, ” in Proc. SPIE , v ol. 101391C, Mar . 2017. [9] L. Roquette, M. M. J.-A. Simeoni, P . Hurley , and A. G. J. Besson, “On an analytical, spatially-varying, point-spread-function, ” in 2017 IEEE International Ultr asound Symposium (IUS) , Sep. 2017, W ashington D.C., USA. [10] P . L. Combettes and J.-C. Pesquet, “Proximal splitting methods in signal processing, ” in F ixed-point algorithms for inverse pr oblems in science and engineering . Springer, 2011, pp. 185–212. [11] N. Parikh, S. P . Boyd et al. , “Proximal algorithms, ” F ound. T r ends Optim. , vol. 1, no. 3, pp. 127–239, 2014. [12] J. A. Jensen, “Field: A program for simulating ultrasound systems, ” Med. Biol. Eng. Comput. , v ol. 34, pp. 351–353, 1996. [13] J. A. Jensen and N. B. Svendsen, “Calculation of pressure fields from arbitrarily shaped, apodized, and excited ultrasound transducers, ” IEEE T rans. Ultrason. F err oelectr . F r eq. Control , vol. 39, no. 2, pp. 262–267, Mar . 1992. [14] M. J. Ackerman, “The visible human project, ” Proceedings of the IEEE , vol. 86, no. 3, pp. 504–511, 1998. [15] H. Zou and T . Hastie, “Regularization and variable selection via the elastic net, ” J . R. Stat. Soc. Ser . B. Methodol. , vol. 67, no. 2, pp. 301– 320, 2005. [16] M. I. Florea and S. A. V orobyov , “ An accelerated composite gradient method for large-scale composite objective problems, ” arXiv pr eprint arXiv:1612.02352 , 2016. [17] ——, “ A generalized accelerated composite gradient method: Unit- ing Nesterov’ s fast gradient method and FIST A, ” arXiv pr eprint arXiv:1705.10266 , 2017. [18] B. P . Lathi, Signal Pr ocessing and Linear Systems . Oxford University Press, New Y ork, 1998. [19] R. C. Gonzalez and R. E. W oods, Digital Image Pr ocessing , 3rd ed. Pearson Prentice Hall, Upper Saddle River , NJ, USA, 2008.