Subtractive NCE‐MRA: Improved background suppression using robust regression‐based weighted subtraction

To correct the intensity difference of static background signal between bright blood images and dark blood images in subtractive non‐contrast–enhanced MR angiography using robust regression, thereby improving static background signal suppression on subtracted angiograms.


| INTRODUCTION
Contrast-enhanced magnetic resonance angiography (CE-MRA) 1 is widely used for the clinical evaluation of vascular anatomy and pathology, 2,3 with the advantage of avoiding ionizing radiation exposure. However, gadolinium-based contrast agents have raised concerns over their association with nephrogenic systemic fibrosis in patients with renal failure, and recently over long-term deposition of gadolinium in the brain even in patients without severe renal dysfunction. 4,5 Such concerns have stimulated the renaissance of non-contrast-enhanced (NCE) MRA techniques, which yield vascular images without injection of exogenous contrast agents. 6,7 Subtractive NCE-MRA is a class of techniques that acquires two image sets (dark blood images [DBIs] and bright blood images [BBIs]) with different vascular signal intensity (SI), which are later subtracted to generate angiograms. Typical subtractive NCE-MRA techniques include flow-sensitive dephasing (FSD) [8][9][10][11][12] and fresh blood imaging (FBI). 13,14 FSD dephases moving blood spins by applying specially designed flow-sensitive preparation schemes. Velocitysensitized (delay alternating with nutation for tailored excitation [DANTE] 15,16 and improved motion-sensitized driven-equilibrium [iMSDE] 8,17 ) or acceleration-sensitized preparations (acceleration dependent vascular anatomy for non-contrast-enhanced MRA [ADVANCE-MRA] 10 ) can be used to suppress either arterial signal or both arterial and venous signals. FBI uses a half-Fourier fast spin-echo pulse sequence to acquire data sets in systole and diastole. The arterial signal is suppressed during systole because of fast flow, which leads to dephasing of moving spins, whereas the signal is preserved in diastole because of slow flow.
In these conventional NCE-MRA techniques, the static background SIs from DBIs and BBIs should ideally be identical, leading to a complete absence of background signal in the subtracted images. However, in practice, because of different imaging settings of the two acquisitions and the varying effective TR between the two acquisitions, some static background tissues show slightly different signal levels in BBIs and DBIs. 9,[18][19][20][21][22][23][24] For example, fat signal, if not sufficiently suppressed, can appear higher on BBIs for both FSD and FBI, leading to a residual fat signal and stripe artifacts in the subtracted images. 9,20,25 Muscle SI is normally higher on BBIs for FSD, 18 but it could be slightly higher on DBIs for FBI. 22,23 These residual background signals can potentially obscure the vascular signal of interest and affect clinical diagnosis. Venous contamination was also reported in a FBI arteriogram, 18,21,24 which can make angiographic interpretation challenging as a result of the close proximity of paired veins alongside peripheral arteries. 6 Residual background signals can potentially be suppressed using a weighted subtraction. Weighted subtraction has been used in some MRA techniques, such as MRA using steadystate free precession (SSFP), 26 hybrid of opposite-contrast MRA (HOP-MRA), 27 and interleaved double-echo MRA and venography, 28 to maximize the blood-to-background contrast. However, the weighting factor was selected manually and empirically. Our initial investigations demonstrated that the weighting factor for the subtraction in FBI and FSD can be obtained adaptively and automatically by performing a linear regression of the SIs. 18,19 However, simple linear regression methods, such as ordinary least squares regression (OLS) and principal component analysis, can be affected by outliers from vascular tissues such as the heart and large arteries. 18,19 Also, some background tissues with high SI on the DBI and BBIs are not given sufficient weighting in the regression procedure(s) to ensure that they are suppressed on the subtracted angiograms.
The purpose of this study is to develop robust regression methods to correct the SI difference of background tissues between the BBIs and DBIs and to improve the background suppression. First, the signal levels of background tissues on BBIs and DBIs are evaluated in several different NCE-MRA techniques. Then, robust regression models, using iteratively reweighted least squares, are proposed to acquire the regression coefficient of the SI of background tissues on BBIs and DBIs-with the weighting function based on either the Euclidean distance or the deviation angle relative to the estimated regression line. Results from these two robust regression models, together with OLS, are compared with reference values subjectively determined by two observers over several different imaging sequences and in different vascular imaging areas. Figure 1A depicts the scatter plot of the SI of each pixel on DBIs versus the SI on BBIs for one single slice ( Figure 1B,C, and 3D, femoral artery FBI). Based on their anatomical location and intensity characteristics, all the pixels were manually categorized into three anatomical types: artery, bladder, and other background tissues. It can be seen that the artery K E Y W O R D S background suppression, MR angiography, non-contrast-enhanced, robust regression signal pixels are mainly located on the left of the scatter map, whereas the background tissue signal pixels distribute along a determined regression line y = x. For conventional direct subtraction, the SI for a pixel P is obtained by subtracting its SI on DBIs from its SI on BBIs (I = I ℎ − I ), which is proportional to the distance between the point P and the green line y = x(line of unity) on the scatter plot (d 1 ). In many cases, the SIs of background tissues appear higher on BBIs than DBIs, which leads to residual background signal on subtracted angiograms because there is still a considerable distance between the background pixel plot and the line of unity (y = x), particularly for background tissues with high SI such as the bladder.

| Characteristics of residual signal
An alternative method that could weight the subtraction to the actual distributions of background SIs would be expected to improve the suppression of background tissues (I = I ℎ − I , no offset was added along the y-axis in all the models used in this study). Using such a method, a regression line (yellow line) for background pixels should be used, and the distance to the regression line would then be d 2 instead of d 1 . Background tissue can thereby be reduced to close to zero, whereas the arterial signal would be less affected, as d 1 and d 2 of arterial pixels are very similar. Figure 2 shows two examples of scatter plots ( Figure  2A,C) with normalized density (shown in the color scale) and the corresponding maximum intensity projections (MIPs) of original raw images, directly subtracted angiograms, and expected angiograms with improved background suppression ( Figure 2B,D). Figure 2A,B are from thoracic 3D FSD (DANTE-balanced steady-state free precession [bSSFP]), 15,25 and Figure 2B,D are from femoral 3D FBI. Reference angiograms with optimal background suppression and vascular visualization were obtained by manual adjustment of the slope of the regression line to match the observed points on the scatter plot. Residual muscle signal (yellow arrowheads) can be seen in the direct subtracted thoracic FSD angiogram, whereas residual venous (yellow arrowheads) and testis signal (blue arrowheads) can be seen in direct subtracted femoral FBI angiograms.
The regression line can be obtained by different linear regression methods, but there are two key challenges to be addressed. The first challenge is the sensitivity to outliers. In some cases, such as thoracic MRA (Figure 2A,B), the presence of the heart leads to many flowing blood pixels, which in the scatter plot lie in the region with high SI on BBIs and low SI on DBIs (the red-dashed region in Figure 2A). These points, which are considered as outliers in the linear regression, have large distances to the model prediction, and thus would be given large weights in models like OLS or principal component analysis. 18 Therefore, the outliers would easily displace the regression line and generate an overlarge slope (the green line in Figure 2A), which can lead to excessive suppression and impair vascular signal.

F I G U R E 1
The scatter plot of the dark blood image (DBI) and bright blood image (BBI) signal intensity (SI) for each pixel in one single slice (A) and the corresponding bright artery (B) and dark artery slice (C). The slice is selected from a femoral 3D fresh blood imaging-magnetic resonance angiography image set. The red points correspond to the artery signal, the purple points correspond to the bladder signal, and blue points correspond to other background tissues such as veins and muscles. Using conventional direct subtraction (green line), the SI of pixel P on the subtracted angiogram is proportional to d 1 , whereas for the optimal subtraction (yellow line), the SI should be proportional to d 2

F I G U R E 2
Two examples of thoracic 3D flow-sensitive dephasing. A and B, Delay alternating with nutation for tailored excitation-balanced steady-state free precession (DANTE-bSSFP), intended to show both arteries and veins) and femoral 3D fresh blood imaging (FBI) (C and D, intended to show arteries only). A and C, Scatter plots of the voxels from the whole 3D data set. The colored scale refers to the pixel number of each point divided by the maximum pixel number among all the pixels, which can be regarded as normalized density or frequency of the points. B and D, Corresponding maximum intensity projections of original raw images, directly subtracted angiograms, and expected angiograms with improved background suppression. The green lines show the results of ordinary least squares regression (OLS), which fails to obtain the regression line of background voxels. The red lines are the expected results determined subjectively by manual selection. SI, signal intensity Second, the model may not be sensitive enough to pixels with high SI. In FBI ( Figure 2C,D), muscle signal appears similar or even slightly lower 22,23 on BBIs than DBIs. Its corresponding pixels lie in the region with low SI on both BBIs and DBIs (the red-dashed region in Figure 2C). In many cases, such as femoral MRA, muscle is a much larger background component than other tissues and thus dominates the regression process, generating a relatively small regression coefficient (the green line in Figure 2C) for OLS. Therefore, the background tissues with high SI on both BBIs and DBIs, such as the bladder and veins, are not fully suppressed.

| Robust regression
For OLS, the linear regression coefficient β can be solved by minimizing the objective function which is given by where n is the number of samples and e is the error of estima- The OLS model assumes a Gaussian aussian distribution of errors e i . However, in the presence of outliers, the long-tail error distribution may lead to a biased estimate of the regression coefficient. Robust regression methods have been proposed to downweight the influence of outliers by modifying the objective function, accommodating more general error distributions, and reducing the sensitivity to the magnitude of the residuals. The most common general method of robust regression is M-estimator (the name comes from the maximum-likelihood estimation), which attempts to minimize the sum of a chosen function ρ(⋅) of the residual errors. 29 The function ρ(⋅) gives the contribution of each residual to the objective function.
Formally defined, M-estimators are given by where τ is a scale parameter determined empirically by the median absolute deviation estimator: The constant 0.6745 makes the estimate unbiased for the normal distribution. 30 The function is minimized by setting the first partial derivatives of ρ(⋅) with respect to β to zero, resulting in a nonlinear equation where (e) = � (e) is called the influence function. The weight function is now defined as Once the format of the weighted function is chosen, the equation can be solved using a numerical method called iteratively reweighted least squares, which iteratively estimates the weighted least squares fit. The steps are as follows: 1. Use the least squares estimate from Equation (1)  where x is the model column vector, with x i as its i-th element, is the current weight matrix. Repeat steps 2 and 3 until the estimated coefficients satisfy the converge criterion: where ε is the floating-point relative accuracy corresponding to the distance from 1.0 to the next largest floating-point number (2 −52 for double precision and 2 −23 for single precision).

| Weight function
The selection of the weight function w(⋅) is the key to achieve robust regression for a specific problem. A typical function in M-estimation, the Welsch redescender 29 is used in this study, which is defined as: The tuning constant c gives coefficient estimates that have 95% asymptotic efficiency with respect to OLS at the Gaussian aussian distribution. 29 Further decreasing its value increases the downweight assigned to outlier points far from the regression line, and vice versa. Figure 3 shows the weight maps of four weight functions investigated in this study. For OLS (or principal component analysis), the weight function is w = 1. The weights are the same everywhere on the scatter map ( Figure 3A), such as point P and Q. In this case, point P has a much larger influence on the regression result because of its longer distance to the current regression line.
For conventional robust regression (cRR) based on the Euclidean distance to the estimated regression line ( Figure  3B), the weight function is given by w (e) = Welsch (e), where e = y − x . The points far from the regression line (point P) are given much smaller weights than the points closer to the line (point Q), which reduces the impact of the outliers corresponding to the arterial signal. However, cRR is not always effective at suppressing tissues with high SI on both BBIs and DBIs. For example, if the points of background tissues distribute along the line y = x ( > 1) (red line), for points R and S, although both are located on y = x, point S has a greater distance to the regression line. Therefore, point S would generate a heavier residual signal on the subtracted image and should ideally be given greater emphasis in the determination of . However, cRR determines the weight based on Euclidean distance and gives point S a smaller weight, making the model less sensitive to the points with large SI values.
F I G U R E 3 Weight maps. A, Ordinary least squares regression (OLS). B, Conventional robust regression using the Euclidean distance (cRR). C, Robust regression using the deviation angle (RRDA). D, Improved RRDA (α = 0.25) when the current regression line is y = x. All of the robust regression models downweight points far from the model prediction. RRDA gives the points with large signal intensity (SI) larger weights compared with cRR. Point P is an example point far from the current regression line, and point Q is a point close to the regression line. Points R and S both locate on the red line and have the same deviation angle, but point S has a larger distance to the current regression line. All of the robust regression models downweight points far from the model prediction (point P). RRDA gives the points with large signal intensity (point S in comparison with point R) larger weights compared with cRR A potential improvement would be to use the polar coordinate system, determining weights from the polar angle instead of cartesian distance. For robust regression based on deviation angle (RRDA), the weight function is given by w ( ) = Welsch ( ), where θ is the deviation angle to the regression line given by = arctan y x − arctan ( ). The contribution of high-SI points can still be too low because of their relatively small number in comparison with low-SI points corresponding to muscle and background air. Therefore, the normalized radial distance r can be used to further increase its sensitivity to points with large values: where ( ) = arctan y x − arctan ( ) and r = √ x 2 + y 2 , α is a parameter controlling the influence of the normalized radial distance. It can be seen that point S has the same weight as R on Figure 3C and outweighs R on Figure 3D. This improved version of RRDA was used for the evaluation in this study. Table 1 summarizes the weight functions and characteristics of the different regression methods mentioned above.

| Study population and imaging protocols
Multiple data sets with differing NCE-MRA acquisition techniques were used in this study, including 36 coronal thoracic FSD-MRA data sets (DANTE-bSSFP) 15,25 from 16 healthy volunteers and 12 patients with central venous obstruction or restricted venous access; 13 coronal iliac FSD-MRA data sets (three DANTE-bSSFP, seven iMSDE-bSSFP, 9,10 three iMSDE-fast spin echo [FSE] 31 ) from six healthy volunteers; and 26 coronal femoral FBI-MRA (electrocardiographically gated FSE) data sets from 17 healthy volunteers. The details of protocols are listed in Table 2.
All the images were acquired using 1.5T MRI systems (Discovery MR450 or Optima MR450w; GE Healthcare, Waukesha, WI). Studies were approved by the local research ethics committee, and all participants gave informed consent.

| Model parameters
In this study, OLS, cRR, and RRDA methods were evaluated for all the data sets. In RRDA, we adopted α = 1 for femoral FBI, which has a large number of muscle voxels with low SIs, and α = 0.5 for thoracic and iliac MRA images.
The 3D data sets in this study have a large number of voxels (8.4 × 10 6 -1.4 × 10 8 ), resulting in a large computational burden for real-time online processing. Therefore, a subset of voxels was randomly sampled for regression to increase computational efficiency. In a pilot study, we evaluated the performance of using different amounts of the randomly sampled data in all the data sets. The number of sampled voxels was increased from 2 10 to 2 23 , multiplying by two in each step, and the error between the results of partial data and full data was calculated. It was found that when the sampled voxels were larger than 2 17 (131 072), the mean error in the regression coefficient caused by this subsampling reduced to <0.01, which is visually undetectable in the intensity corrected angiograms. Therefore, the number of sampled voxels was fixed to 5 × 10 5 (500 000) for all data sets in this study. This corresponded to sampling ratios of 4.88%, 3.70%, and 1.15%, and an estimated mean error of 0.0023, 0.0017, and 0.0017 in the regression line slope for thoracic FSD, iliac FSD, and femoral FBI, respectively. Compared with using full data, the use of partial data reduced the average computation times of thoracic FSD, iliac FSD, and femoral FBI data from 20.2, 30.0, and 55.2 seconds to 0.6, 1.0, and 0.9 seconds, respectively (4-core 3.4 GHz CPU and 16 GB RAM).

| Image assessment
To evaluate the background SI difference between BBIs and DBIs, the SIs of specific background target tissues, such as  for the diastolic acquisition) and flow-spoiled gradients (10% of one-half the area of the readout gradient for both systolic and diastolic acquisitions) 14 were both used to increase flow sensitivity Abbreviations: cRR, conventional robust regression using the Euclidean distance; FBI, fresh blood imaging; FSD, flow-sensitive dephasing; OLS, ordinary least squares regression; RRDA, robust regression using the deviation angle. muscle, liver, bladder, testis, and veins (only in arteriograms) were measured on the individual raw images before subtraction. Matched regions of interest (ROIs) were drawn in representative regions of target tissues. One ROI was drawn on three selected slices for each data set and for each background tissue. The slices were selected randomly, but with the requirement of a large coverage of target tissues. The mean value of each background tissue was calculated for each data set. The calculated regression coefficients were compared with reference values determined by two trained observers with more than 4 years' experience in vascular MRI. Subjective manual determination of regression coefficients was performed based on both MIPs of subtracted images and scatter plots, using an interactive graphical interface developed in MATLAB (R2019a, MathWorks Inc., Natick, MA). Manually assessed optimal regression coefficients were determined by two observers independently, based on the following method: The regression line was first positioned in the center of the distribution of background signal voxels on the scatter plot and then adjusted to achieve optimal background signal suppression without impairing vascular signal on the MIP. The reference regression value was taken as the mean of the two manually determined optimal regression coefficients.
The results obtained by each regression model were compared with the reference values over different subjects using a paired, two-sided Student t test. Pearson product-moment correlation analysis and Bland-Altman analysis were also used to evaluate the agreement between the results of each regression method and the reference values. A P value < .05 was considered statistically significant.
SI ratios of background tissue to vascular signal were also calculated for images obtained by direct subtraction and RRDA. SIs were measured from matched ROIs drawn in representative regions in the target tissues on the MIPs of subtracted images. For FBI-MR arteriography, two vascular ROIs were drawn in the lumen regions of arteries for each MIP; for FSD-MRA visualizing both arteries and veins, one vascular ROI each was drawn from the lumen regions of arteries and veins, respectively, and their mean value used as the vascular SI. Two ROIs were drawn for each type of background tissue, and the mean values were calculated for each data set. Figure 4 shows the signal levels of representative static tissues on DBIs and BBIs in the different study populations. SI differences can be seen in all these tissues, especially for organs containing plenty of water such as bladder and liver. The signal levels of most static tissues are significantly higher on BBIs than DBIs (P < .05), leading to residual background signal on the subtracted angiograms. However, for femoral FBI, there is no significant difference between the muscle SI on BBIs and DBIs, which can bias the regression towards the distribution of other tissues such as the bladder and veins. Figure 5(A-E) shows an example of thoracic FSD-MRA (DANTE-bSSFP) from a patient with central venous obstruction. The performance of RRDA is compared with OLS and direct subtraction. OLS (yellow solid line in Figure 5E) was affected by outliers and overestimated the regression coefficient, leading to signal loss in small vessels (arrowheads in Figure 5A-D). In contrast, both RRDA and cRR generated the correct regression line (red and green solid line in Figure 5E) and thus suppressed residual muscle signal. while preserving F I G U R E 4 Signal levels of background tissues in bright blood images (BBIs) and dark blood images (DBIs). A, Thoracic flow-sensitive dephasing (FSD). B, Iliac FSD. C, Femoral fresh blood imaging (FBI). Asterisk (*) indicates significance at P < .05 level. Signal levels are significantly lower on DBIs in all cases except for the muscle signal in femoral FBI the blood signal ( Figure 5C,D). Similar results can be found in Figure 5F-J showing an example of iliac iMSDE-FSE from a healthy volunteer.

| RESULTS
An example of femoral FBI is shown in Figure 5K-O. The regression in cRR and OLS was dominated by the muscle voxels with a large number of voxels and low SI: Their regression coefficients were insufficient to suppress veins (yellow arrowheads in Figure 5L,M) or background tissues (blue arrowhead in Figure 5L,M). RRDA was less sensitive to the muscle voxels (red solid line in Figure 5J): Background tissues such as testis (the blue arrowhead in Figure 5N) and veins (yellow arrowheads in Figure 5N) were suppressed, whereas the arterial signal is preserved. Note that some small vessels denoted by the green arrowheads in Figure 5M are actually veins, which can be easily confused with arterial branches, but can be suppressed by RRDA.
The statistical analysis results are listed in Table 3. Figure 6 shows the Bland-Altman plots and boxplots comparing the performance of different regression methods for different types of images. The iliac data sets using different sequences were categorized into one category because of their similar image characteristics and regression results.
OLS has high mean absolute errors and low correlation coefficients in all three image types. In particular, based on its sensitivity to outliers, OLS tends to overestimate the regression coefficients in thoracic FSD, leading to a large bias of 0.293 and also large variation. cRR is more resistant to outliers. It has excellent performance in thoracic and iliac FSD. However, cRR has a bias of −0.089 in femoral FBI. The regression coefficients of cRR are smaller than the reference values in all the femoral FBI data sets evaluated in this study. RRDA achieved accurate and robust results in all the three types of images, with small bias, small mean absolute error, and high correlation coefficients. No significant differences were found between the RRDA results and reference values for any of the three methods.
A comparison of the tissue-to-vessel SI ratios between RRDA and direct subtraction is shown in Figure 7. In thoracic F I G U R E 5 Example maximum intensity projections and the corresponding scatter plots of thoracic 3D delay alternating with nutation for tailored excitation-balanced steady-state free precession (A-E), iliac 3D improved motion-sensitized driven-equilibrium-fast spin echo (F-J), and femoral 3D fresh blood imaging (K-O). E, J, and O, Scatter plots of the voxels from the whole 3D data set with normalized density as the color scale. The grey dashed lines on the scatter plots are the 1:1 line for the direct subtraction; the yellow, green, and red solid lines are the regression lines estimated by ordinary least squares regression (OLS), conventional robust regression using the Euclidean distance (cRR), and robust regression using the deviation angle (RRDA), respectively. The red boxes denote the regression results that agree most closely with the reference value. The arrowheads on maximum intensity projections denote the suppressed background tissue signal. The green arrowheads on M denote small venous branches that are easily confused with arterial branches, but can be suppressed by RRDA. RRDA had good regression results in all the three cases. OLS generated an overlarge regression coefficient in thoracic and iliac flow-sensitive dephasing (FSD), which led to signal loss in small vessels. cRR failed to suppress the residual background signal in femoral FBI, but agreed with RRDA in the other cases FSD, the muscle and liver-to-vessel SI ratios in RRDA reduced to 54.4% (from 0.18 ± 0.08 to 0.10 ± 0.05) and 93.0% (from 0.50 ± 0.23 to 0.47 ± 0.25) of their values in direct subtraction, respectively; in iliac FSD, the muscle-, bladderand testis-to-vessel SI ratios reduced to 63.1% (from 0.20 ± 0.10 to 0.13 ± 0.07), 76.5% (from 0.51 ± 0.18 to 0.39 ± 0.18), and 62.2% (from 0.35 ± 0.20 to 0.22 ± 0.15), respectively; in femoral FBI, the muscle-, bladder-, testis-and vein-to-artery SI ratios reduced to 36.0% (from 0.44 ± 0.16 to 0.16 ± 0.07), 35.2% (from 0.21 ± 0.18 to 0.07 ± 0.04), 45.1% (from 0.26 ± 0.07 to 0.12 ± 0.04), and 41.5% (from 0.20 ± 0.08 to 0.08 ± 0.06), respectively.

| DISCUSSION
This study has developed a robust regression method to compensate for the SI differences of static background tissues between BBIs and DBIs in subtractive NCE-MRA. Signallevel differences were found in almost all the static tissues for the subtractive NCE-MRA techniques investigated in this study. These differences lead to residual background signal when using direct subtraction, but this residual signal can be substantially reduced using weighted subtraction, if an appropriate weighting factor is determined, corresponding to the regression coefficient between the background tissue signals.
FSD produces DBIs using flow-sensitive preparation modules to suppress blood with fast flow. Although these modules are velocity-dependent, they can still impair the signal level of static residual background tissues. With MSDE or iMSDE, some static signal loss is inevitable, resulting from inherent T 2 decay, T 1 steady-state decay, diffusion, 16 eddy currents from crusher gradients, and imperfections in the 180° pulses. DANTE pulse trains cause a spoiling effect in flowing spins and suppress their signal. Static spins may form a steady state, but still experience some signal decay, especially when eddy current effects distort the gradient waveforms and thus impair the phase preservation of the static signal. 16

F I G U R E 6
Bland-Altman plots and boxplots comparing the regression results of ordinary least squares regression (OLS), conventional robust regression using the Euclidean distance (cRR), and robust regression using the deviation angle (RRDA). FBI, fresh blood imaging; FSD, flowsensitive dephasing

F I G U R E 7
Comparison of the tissue-to-artery/vein signal intensity ratios between robust regression using the deviation angle (RRDA) and direct subtraction in thoracic flow-sensitive dephasing (FSD) (A), iliac FSD (B), and femoral fresh blood imaging (FBI; C). Statistically significant differences can be seen in all the cases (P < .05) FBI attenuates flowing spins caused by flow dephasing and flow void effect of the FSE sequence. Variable flip angles 32 or spoiling gradients 14 are normally employed to increase its sensitivity to blood flow velocity. Unfortunately, both of these approaches can also slightly reduce the signal from background tissues containing water or blood, such as veins, bladder, and testis. The SI differences between BBIs and DBIs can be more serious when the variable flip angles and flow-spoiled gradients are only employed in systolic acquisitions.
Another factor leading to the background SI difference is the varying TRs of systolic and diastolic acquisitions, 21,22 which may be caused by two effects. First, the FBI sequence in this study uses an interleaved black/bright blood acquisition, which alternates between brightblood and dark blood acquisitions in consecutive heartbeats. Because the diastolic acquisition has a longer trigger delay, that is a longer period of T 1 recovery, background tissues would have slightly higher SI on the diastolic images than the systolic images. This is especially obvious for tissues with a long T 1 and T 2 such as the bladder, as their signals are not fully recovered within a TR of two or three R-R intervals and thus greatly fluctuate between BBIs and DBIs. In contrast, the original FBI technique 6,14 uses a sequential scheme that acquires consecutive slice encoding in bright blood and dark blood series to maintain the same TR and avoid this problem. Also, the FSD sequences used in this study acquire all images with the same cardiac trigger delay and do not have this issue. Another reason for varying TRs is the large R-R variations caused by the irregular heartbeats of subjects. This is more serious for the sequential acquisition used in the original FBI technique and the thoracic FSD sequence in this study. The influence of irregular heartbeats on background artifacts has been evaluated in volunteers with caffeine consumption or with cardiac arrhythmia by Kim et al. 21 The sensitivity to outliers is the main limitation for using OLS. For a regression biased towards background tissues, vascular and heart voxels with different signal levels on BBIs and DBIs can be regarded as outliers and can affect the regression. cRR using iteratively reweighted least squares can minimize the effect of outliers by assigning a weight to each data point according to its distance from the regression line. The results of our study show that cRR achieved accurate and robust results for thoracic and iliac FSD-MRA data. However, in femoral FBI, the muscle signal, which has similar low SI on BBIs and DBIs, has overwhelmingly the largest number of voxels and thus dominates the regression. To solve this problem, we adopted a polar system and deviation angle (in RRDA) instead of the Euclidean distance used in cRR. The reasons are as follows: First, on the scatter plot, the SI difference between BBIs and DBIs is reflected as a deviation angle between the distribution of background points and the line y = x (rotation), rather than a straight distance (translation). Second, compared with low SI tissues such as muscle, high SI tissues have higher residual SI on subtracted images and therefore should have a larger weight to achieve an optimal background suppression, but they are given smaller weights in cRR because of their longer distance to the regression line. A normalized radial distance with an exponential parameter α was added into the weighting function of RRDA to further increase its sensitivity to the points with large values. To achieve the optimal performance, the parameter α may need to be adjusted according to the voxel number of low SI background tissues.
Although background tissues were suppressed using weighted subtraction, the desired vascular signals were only minimally affected. In this study, obvious vascular signal loss was only seen when the weighting factor is overestimated, for example, the OLS results in thoracic FSD-MRA. This is because of the distribution of vascular signal and background signal on the scatter plot as can be seen in Figure 1A. For the arterial pixels in Figure 1A, which are located on the bottom left of the scatter plot, their distance to the line y = x (d 1 ) and the regression line (d 2 ) are very similar, whereas d 1 is much larger than d 2 for background pixels located on the upper right. Therefore, using d 2 instead of d 1 does not change arterial SI much, but can substantially reduce the SI of background tissues.
It should be noted that in conventional subtractive NCE-MRA techniques, background static signal suppression in unsubtracted raw images is still necessary. For example, short-time inversion recovery, magnetization transfer contrast, 33 and longer TE 32 have been used in FBI for background suppression to reduce the SI of background tissues and thereby any residual background SIs on subtracted angiograms. The weighted subtraction method developed in this study gives subtractive NCE-MRA techniques a larger tolerance of background SI difference between BBIs and DBIs. This potentially allows the imaging sequence design to maximize the arterial signal even at the cost of increasing background signals, eg, using minimum TE in FBI and selecting the flip angle maximizing the arterial signal in bSSFP. Another potential sequence design is to generate images with maximum arterial signal, but different background SIs on BBIs and DBIs, for example, using strong flow dephasing for the systolic acquisition and weak flow dephasing for the diastolic acquisition in FBI.
Another potential advantage of reducing the static background signal in angiograms is that it improves image sparsity, which we have recently shown can facilitate compressed sensing reconstruction and potentially improve reconstruction accuracy. 34 Although the regression results and background suppression effect of robust regression methods were analyzed in volunteer and patient images, the diagnostic performance was not evaluated in this study. Future work will assess if clinical diagnosis or diagnostic confidence can be improved by background suppression using robust regression. Further work is also needed to evaluate the methods in other subtractive imaging sequences and anatomical regions.

| CONCLUSIONS
In conclusion, robust regression approaches have been developed to correct the SI difference between background tissues on BBIs and DBIs for several subtractive NCE-MRA techniques, reducing residual background signal on subtracted angiograms. Compared with OLS and the conventional robust regression method, RRDA proved more resistant to outliers across the range of different NCE-MRA methods investigated and more sensitive to background signal with high intensity. In this initial study, it achieved an accurate and robust performance applied to several different subtractive NCE-MRA techniques and different body regions.