SPATIAL AND SPECTRAL ANISOTROPY IN HARP IMAGES: AN AUTOMATED APPROACH

Strain and related tensors play a major role in cardiac function characterization, so correct estimation of the local phase in tagged images is crucial for quantitative myocardial motion studies. We propose an Harmonic Phase related procedure that is adaptive in the spatial and the spectral domains: as for the former, we use an angled-steered analysis window prior to the Fourier Transform; as for the latter, the bandpass ﬁlter is also angle-adaptive. Both of them are narrow in the modulation direction and wide in the orthogonal direction. Moreover, no parameters are manually set since their values are partially based on the information available at the DICOM headers and additional information is estimated from data. The procedure is tested in terms of accuracy (on synthetic data) and reproducibility (on real data) of the deformation gradient tensor, measured by means of the distribution of the Frobenius norm differences between two tensor datasets.


INTRODUCTION
Magnetic Resonance-Tagging (MR-T) is based on the generation of a set of saturated magnetization planes on the imaged volume which may be subsequently tracked throughout the cardiac cycle [1].This modality is of special relevance in the analysis of myocardial motion to provide local functional indicators, such as the strain or the strain rate tensor [2].
Regarding the analysis of MR-T, an important family of methods are based on the extraction of the complex image phase obtained by bandpass filtering the spectral peaks introduced by the applied modulation; these methods are based on the fact that the HARmonic Phase (HARP) is linearly related to a directional component of the true motion [3]; these methods are capable of reconstructing dense displacement This work was partially supported by the Spanish Ministerio de Ciencia e Innovacion under Research Grant TEC2013-44194-P, the Spanish Ministerio de Ciencia e Innovacion and the European Regional Development Fund (ERDF-FEDER) under Research Grant TEC2014-57428-R, the Spanish Junta de Castilla y Leon under Grant VA136U13 and by the Universidad de Valladolid-Banco de Santander grant program.fields accurately grounded on the assumption of constant local phase, which turns out to be more reliable than a constant pixel brightness assumption.
Subsequent works have focused on the nonhomogeneous deformation of the stripe pattern and have made use of the Windowed Fourier Transform (WFT) [4] (hereafter referred to as WHARP) to balance the spatial and spectral localization of the image to be filtered by means of a truncated Gaussian window (as an appropriate trade-off between spatial and frequency resolutions); the result of which is a smooth local phase estimation.However, in this approach, Gaussian windows used prior to the WFT are isotropic, which constitutes a limitation on the capability to adapt to the locally changing stripe pattern.This Gaussian window has been made adaptive in [5], hereafter referred to as AWHARP, to accommodate tag local properties but, with such a design, performance is dependent of the absolute tag orientation of the original stripe pattern, which is a suboptimal design choice as well.
Once the local spectrum is calculated, a HARP bandpass filtering stage should be performed to extract the phase.A circumferential spectral filter centered at the location of the maximum of the spectrum is proposed in [6] following the assumption of an isotropic estimation procedure.In [3] a smooth ellipse is chosen for the bandpass region of the filter because it has a simple geometry appropriate for the gross shape of the spectral peaks.This approach implicitly makes use of the concept of anisotropy, but key parameters, such as the bandwidth or the eccentricity of the ellipse, are manually set.
Thresholding methods based on variance [7,8] or entropy [9] could provide automatic estimates of the appropriate filter parameters.These methods, however, should be complemented with means of producing smooth transfer functions that are capable of dealing with the nonhomogeneities caused by pattern deformations.To the best of our knowledge, the filtering schemes proposed in the HARP literature are not capable of automatically adjusting to these spectral peaks, specially in the case of the aforementioned isotropic windowing procedures where the spectral overlapping caused by the anatomical image could result in an artifacted reconstruction.
In this paper we will show that the assumption of anisotropy along the processing pipeline (both at windowing prior to the FT and at the filtering stage) is beneficial in order to avoid spectral interferences without losing deformation information.In addition,we also estimate the filter parameters direclty from data, so our procedure leads to a fully automatic motion estimation tecnique.The method here proposed will be hereafter referred to as automatic anisotropic WHARP, i.e., AA-WHARP.

Materials
We have acquired a medial slice of a MR SPAMM SENSitive Encoding (SENSE) Turbo Field Echo (TFE) sequence on a Philips Achieva 3T scanner.The images have a spatial resolution of (1.333x1.333)mm 2 and a slice thickness of 8 mm.The acquisition parameters are T E = 3.634 ms, T R = 6.0182 ms and α = 10 • .Tag spacing has been set so that k = 1/λ with λ=7 mm; multiple tag orientations have been obtained that fully span the orientation plane, uniformly from 0 to π with a separation of π/18 radians.We have also acquired a SENSE balanced TFE cine sequence with spatial resolution of (1.25x1.25)mm 2 , a slice thickness of 8 mm, T E = 1.663 ms, T R = 3.325 ms and α = 45 • .The cine sequence is acquired at the same slice location as the tagging sequence.The myocardium is segmented in the end-diastole (ED) and endsystole (ES) phases of the cine sequence and in the ED phase of all the acquired orientations of the tagging sequence.The cine segmentation at ED phases is used to align the tagging orientations to a common reference system to correct for patient breathing motion.The ES segmentation is used to define a region of interest (ROI) on which to compute meaningful measures of the strain tensor.
In order to have a deformation ground truth, heart motion has been estimated from the cine sequence by means of a free-form deformation (FFD) [10].Then, an undeformed Tagging modulation has been added at ED phase to the original cine sequence and the result has been deformed by the synthetic transformation to give rise to a simulated Tagging sequence with motion parameters known beforehand.As it is well-known, FFDs are obtained by manipulating an underlying mesh of control points.Their spacing as well as the B-spline [11] order act as parameters of the FFD [10]; in our case, we have used a dense control point mesh as well as a high order B-spline (control points spacing is set to 5 mm and a 3 rd order B-spline has been selected) so that our model can accommodate highly local nonrigid deformations.Additive noise with zero mean and tag fading have also been included in the simulation.

Methods
As stated in Section 1, the original HARP makes use of a complete FT although the tag pattern in MR-T images suffers local variations, i.e., nonhomogeneities.To alleviate this limitation, WFT methods have been proposed both with fixed size [4] and adaptive spatial windows [5].
The anisotropic windowing technique here proposed is based on the extraction of the local phase of the stripe pattern by applying a WFT along the stripe direction combined with a full FT on the direction perpendicular to the stripe over the image at ES phase.Specifically, it is well known [12] that the equation of the straight line that is orthogonal to the unity vector n = (cos(θ), sin(θ)) and for which the distance to a point p = (p x , p y ) is s, turns out to be Therefore, if the analysis window is designed as: where n is the tag unity vector read from the DICOM headers and p = (p x , p y ) the point at which the local FT is calculated, then the window will have a rapid tapering, controlled by parameter σ, in the modulation direction but will suffer no attenuation in the orthogonal direction.For the window design, σ has been set equal to the nominal tag spacing (as read from the DICOM headers).Once the local spectrum is calculated, HARP bandpass filtering techniques can be directly applied on the spatially localized spectrum of the image to extract the phase.The assumption of anisotropy is a key fact, considering that the spectrum will suffer a significant stretching orthogonal to the modulation direction as compared with the isotropic counterpart, caused by the anisotropic window (see Fig. 1).In addition, spectral interferences seem mitigated in the anisotropic version.To adequately retrieve the shape of the spectral peak, we have resorted to an anisotropic filter composed by a Gaussian band-pass filter along the modulation direction and an allpass filter along the orthogonal direction.This is done by means of a spectral mask as defined in Eq. 2 where, as before, n is the tag unity vector read from the DICOM headers and p is the location of the maximum of the spectrum inside a region in the surroundings of the reference spatial frequency of the tags k i , defined by k : ki 2 ≤ k ≤ 2k i ∧ |θ − k| ≤ π 6 .For the filtering stage, σ has been estimated by means of the Otsu's method [7], i.e., by comparing with an spec-tral threshold that is calculated minimizing the bimodal intraclass variance.To this end, we define a rectangular region in the surroundings of the spectral peak p according to k : |s| < ki 2 ; once the threshold is calculated from the information within that area, a foreground region is obtained.However, as it is only necessary to estimate the width of the Gaussian filter in the modulation direction, the foreground points are projected over that direction, the sample standard deviation is calculated and σ is set as four times that measure.This way we assure that approximately 99% of these points fall within an amplitude 1 √ 2 times the maximum of the filter.Finally, the local phase can be extracted in the spatial domain from the inverse anisotropic WFT of the aforementioned filtered spectrum.Once the local phase images are reconstructed, standard procedures are applied to estimate the material deformation gradient tensor F, at ES phase [3,13].

Validation
Validation is carried out in terms of the variability of the material deformation gradient tensor estimate F at the ES phase; as for accuracy, this variability is calculated between the ground truth and the tensor estimates obtained from each pair of orthogonal stripes used.As for reproducibility, variability is calculated out of the tensor estimates obtained from two different sets of orthogonal stripes pairs.Ideally, the tensors should be equal for all datasets and estimation procedures.Therefore, a natural measure for our purposes is the Frobenius norm difference (FND) of the tensors estimates, which is defined as: where F j m,n (x) stands for the (m, n) component of the tensor F j at pixel x, calculated from method (or dataset) j (j = {1, 2}).
Recall from section 2.1 that 18 stripe directions have been imaged; this gives rise to 9 pairs of orthogonal stripes as well as 9  2 = 36 pairwise comparisons.Therefore, as for the reproducibility experiment on real data, we have stacked together all the FND values for all pixels and all pairwise comparisons.For the experiment on synthetic data, the stack consist of 9 comparisons to the ground-truth tensor in every pixel.

RESULTS
Figure 2 shows the cumulative distribution of F N D, as defined in Eq. 3; the plots are indexed by the method used for the tensor estimation.For the non-automatic methods under test, the filter bandwidth is normalised with respect to the wave number of the applied modulation using the parameter µ = r/k.This parameter has been set to 0.35 according to [6] for both experiments.The figures show that both in terms of accuracy and of reproducibility, the proposed method outperforms previous proposals, i.e., the F N D distributions are closer to the ideal distributions.Notice that no dependence on the orientation is observed since data from different orientations are stacked together.Figure 4 does show such a dependence.Specifically, for each orientation indicated in the horizontal axis (together with its orthogonal accompanying direction) we have calculated the F N D for each of the other 8 orientation pairs; data samples have been stacked together and then the ν(F N D) is calculated.Clearly, the method here proposed shows less orientation sensitiveness than previous proposals.

DISCUSSION
From the results in Section 3 we can see that the methodology here presented is effective in improving both robustness and accuracy when obtaining direct strain measurements from MR-T sequences; as indicated, the distribution of F N D is closer to the ideal pattern and no crossing between curves and our proposal are observed.Interferences, however, so still remain, so further research in this direction should be considered.In any case, the notion of anisotropy both in spectral and  spatial domains seems beneficial with respect to other adaptive and/or isotropic approaches.
Figure 3 shows a considerable improvement in terms of ν(F N D) of reproducibility and precision for the AA-WHARP method with respect to reported methods for a wide range of µ values, specially with higher bandwidths.
The proposed anisotropic design leads to an automated bandwidth estimation method resulting in a smooth Gaussian filter that outperforms other HARP based methods when reconstructing the tensor with arbitrary orthogonal pairs of orientations.
It is worth noting that, as indicated in Figure 4, AWHARP is sensitive to the pattern orientation despite its adaptive condition; the fluctuations in ν(F N D) have considerable higher peaks than those from our AA-WHARP; it should be highlighted that AWHARP shows very similar behaviour to WHARP in the vicinity of 45 • since at this (nominal) orientation the analysis window in AWHARP is close to isotropic.This fact is reflected in both figures.The proposed method, however, remains fairly orientation independent, but for residual interference and interpolation effects.

CONCLUSIONS
A robust and reliable automatic motion estimation technique is achieved thanks to an anisotropic design all along the processing pipeline.The windowing procedure proposed consists in an anisotropic Gaussian window whose shape remains unaltered but it rotates according to the nominal pattern orientation.Since the spectral peaks undergo a significant stretching, we resort to an anisotropic filtering procedure to diminish spectral interference.The proposed spectral design, in addition, lends itself to a full automatic bandwidth estimation method.Results indicate an improving performance with respect to reported methods for a wide range of filter bandwidths.

Fig. 1 :
Fig. 1: Effects of the WFT design over the image spectrum.

Fig. 2 :
Fig. 2: CDF of F N D for the sets considered in Section 2.1.