As illustrated in Figure.1, serial images derived from multiple spin-echo scans are utilized as input. The processing can be considered in two ways. As diagrammed on the top line, the first image is subjected to operations including smoothing and segmentation. On the bottom line, a weighted least square fitting method is used to produce the T2 parametric image, which is converted into a confidence image based on the Gaussian kernel derived from the phantom study. Finally, the confidence image is cast into the different regions in the segmented image, and groups of fat pixels are distinguished according to their regional confidence scores.
The key techniques used here will be described in detail according to the processing steps including (A) the fitting of parametric image, (B) segmentation of first image, (C) evaluation of phantom and (D) fat extraction.
A: Fitting of parametric image
To measure the transverse relaxation time for each pixel of the parametric image, a curve is fitted to the decay of intensity with increasing TEs [9]. The non-linear least squares is most commonly adopted [10, 11]. For the fitting of fat data, we used the weighted least square method with baseline subtraction and a least point constraint.
a1.Fitting model
After testing both mono and bi-exponential models, the transverse relaxation time measurement of fat appears to satisfy a mono-exponential physical model known as: [12]
Si(S0,T2)=S0*exp(-Tei/T2)
with Tei=i*te and S0 is the pseudo-proton density, which is relative to the true proton density, T1 value and receiver coil response.
a2.Weighted least square
To fit a set of experimental data to the parametric model, a cost function and optimization method was selected.
For the cost function, a least square method (LS) is the common way of fitting the curve, which consists in minimizing the quadratic distance Φ2 between the fitted curve to the curve represented by the raw data.
Where Ii is the scanned intensity of intensity images in ith echo.
Weighted least square method (WLS) takes other factors into account using the weight with merit function:
Where wi, the weight of ith point, should represent the confidence of the signal. As the low intensity signals are more likely to be noise and according have less certainty, here wi was simply set as proportion to the intensity (wi= Ii).
Second, for the optimization, a Marquardt-Levenberg algorithm[13], which is specially designed for multi-dimensional non-linear least squares fitting was selected.
B: Image segmentation of first image
Considerable effort has been devoted to the intesnity image with the aim to segmenting the MRI reconstructed image into tissues. Unfortunately, inhomogeneity, bias, edge blurring and other ill-posed problems require complicated mathematical models. However, unlike the previous methods, only regions within one tissue is our goal for the segmentation (i.e. inhomogeneous fat tissue may be divided into several regions). These regions served as the casting board for the confidence image and fat i ultimately determined by the confidence scores.
To improve the precision our morphological segmentation, bilateral filtering and mean shift methods were used and described below. Here, the image obtained in the first echo is chosen due to its higher SNR relative to the other echoes.
b1.MRI filtering
MRI filtering was implemented as a preliminary step to decrease the noise, which replaces the signal of a pixel according to the neighbouring pixels. Filtered image If can be regarded as a convolution between original image I0 and kernel K [14]: If = K ⊗ I0. The filter kernels, typically a matrix of size M*M, represents the number of pixels nearby taken into account. Each element ky (i,j∈[1,M]) at different positions represents the weight at a given point.
Linear and nonlinear are two typical filters used in image processing. In linear convolution filters, the weighting coefficients ki,j only takes into consideration the relative position in kernel K, and remains constant throughout the whole image filtering. Nonlinear filters are relative to the target pixel and the coefficients are calculated as a function of local variations of the signal [15]. In the linear filter class, average and Gaussian filers are often used.
Among the nonlinear filters, the median filter is popular. As well, a selective blurring filter [16] is used [11, 17], which emphases the pixels with similar intensity to the target pixel. A bilateral filter [18] is an edge preserving technique proposed recently and has been widely used in image processing. In comparison to the selective blurring filter, not only is the intensity similarity taken into account, but also the spatial similarity after processing with a uniform or Gaussian kernel. A comparison of the effects of different filters on the image data is demonstrated later.
b2.Mean shift segments
After smoothing, the first image is segmented into different regions by a mean shift algorithm. Mean shift [19, 20] is a nonparametric estimator of density, which can cluster all the pixels according to both their distribution in space and their intensity thus creating a unique feature space. Given a data point x in the first image, after the first shift, the new vector x(1) is:
Where xi are data points around x and the function g() is related to a kernel expression, which defines different influences of xi according to their distance from x. The parameter h is the bandwidth used to control the kernel size.
The shift is iteratively performed until x converges to a stable mode point. Because the similar data points will converge to the same or nearby mode points, all the data points in the first image can be segmented by merging the mode points based on their distances. More details of the mean shift segmentation can be found in [20]
Compared to conventional MR segmentation methods, mean shift is better for analyzing fat images. Mean shift is a local method making it insensitive to non-uniform intensities. Second, unlike the K-mean or EM methods [21], the number of clusters does not need to be predefined. In addition, the bandwidth parameters can be adjusted for the different sizes of kernel functions, which provide a flexible way to define the scale according to research objectives. After this segmentation process is complete, the first image will be divided into regions.
C:Phantom evaluation
The fitted T2 values of fat pixels in a given region are not exactly the same but rather present a distribution. To evaluate this distribution, a phantom study was done.
c1. SNR effect on T2 distribution
From previous literature [22], it is known that the signal to noise ratio (SNR) can influence this distribution, whether or not variances are caused by the animal or acquisition technique. We excluded variance of individual animals, only the factors related to the acquisition technique were taken into account, which included the effect of the imaging protocol and the instrument.
From the aspect of the MRI protocol, SNR is decided by many parameters, including ratio of echo time, spatial resolution, thickness of slices, and receiver band width. From the instrumentation aspect, it is more complicated to understand the physical factors, which include coil setup, magnetic field bias, RF in homogeneities, and phase deviation. Therefore, to evaluate the relationship between SNR and T2 distribution is a complicated task. A practical way is to determine SNR effects is to use a homogeneous sample with the similar acquisition parameters in the same instrument [22].
c2. Phantom simulation
A tube filled with uniform lard was used as the phantom (Figure 2) under the identical scan protocol as used in the animals. To simulate the noise effects, white Gaussian noise with increasing normalized variation was added to the original phantom images. With the increasing step size of 0.0001, 100 trails were performed and SNR decreased from 36.48 to 2.46 according to the definition:
SNR=avg(S1(S0,T2))/σ.
Where σ denotes the root mean square (RMS) noise in the background and avg() is the average signal value in lard area.
The noise is assumed to follow a Gaussian distribution when SNR>2 because the actual non-zero average Rician noise will become a quasi-Gaussian distribution in both real and imaginary components [12].
In the simulated data, weighted least square fitting was utilized for pixels in the center area of tube as defined by a manually selected ROI. The calculated T2 values were recorded as shown in Figure 3 (All values exceeding 200 were trimmed). The results indicated that when SNR<5, the images were corrupted by noise and T2 values were obviously erroneous. Thus, the standard deviations of T2 are only calculated on valid results in Figure 4. An exponential function was fitted and served as the calibration curve for the confidence image.
D: Fat extraction with confidence image
d1. Confidence image
From the phantom evaluation, the measured T2 value in the fat region varies relative to the SNR. Assuming a Gaussian distribution of T2 values, a SNR based Gaussian kernel is utilized here to draw the confidence image. With this kernel, mean and standard deviation need to be determined.
In the T2 parametric image, the T2 histogram shows the different T2 distributions, where different peak mean the different tissues with different biological characters. The T2 value in a peak point means the highest distribution of one region. It is possible to approach the real value. Here, the mode is found by performing a 1-D mean shift method.
With the Gaussian kernel, a weighted filter was used with the calculated T2 value pixel by pixel. For each pixel P, the cast weight wp can be calculated as:
wp=exp(-(T2,p-T2,m)2/(2σSNR))
with T2,p is the T2 value of P, T2,m is the mode value P belongs to, and σSNR means the standard deviation value corresponding to current SNR determined in the phantom study.
From Figure 5, it can be found that the higher SNR, the narrower the kernel will be, which means the T2 value distribution will be closer to the peak value. Also, after the weight, the pixel with the T2 value near the peak value will approach 1, and those with a difference more than four standard deviations will make no sense in the confidence image.
Therefore, the confidence images are defined by the convolution of T2 and SNR based Gaussian kernel: Ic=Iw ⊗ w.
d2. Fat extraction
Finally, confidence images are cast into the segmented first images pixel by pixel. And confidence scores of each region are summed by the confidence associated with pixels located in that region.
The probability of each region belonging to adipose tissue is determined by the confidence scores. Fat can be extracted by setting a predefined or real-time threshold on these confidence scores.
