CN104392422B - A kind of MRI non homogen field methods of estimation based on image gradient fitting - Google Patents

Abstract

The present invention proposes a kind of MRI non homogen field methods of estimation based on image gradient fitting, processes image in log-domain first, and actual image signal and non homogen field signal are decoupled.It is fitted by the gradient in the more uniform tissue regions of the gray scale in image, obtains the estimation of non homogen field.Due to carrying out gradient treatment in log-domain so that unrelated with organization type where this pixel in the Grad obtained by a certain pixel, and it is only relevant with non homogen field.Devise special object function in fitting to allow to directly obtain non homogen field estimation by the x directions of image and the first derivative in y directions, without integrating or calculating second order gradient again.

Description

MRI (magnetic resonance imaging) inhomogeneous field estimation method based on image gradient fitting

Technical Field

Non-uniform fields in the magnetic resonance image due to magnetic field non-uniformity or coil sensitivity are quickly estimated and corrected.

Background

Magnetic resonance imaging is widely used in medical diagnosis and scientific research due to its high resolution, no ionizing radiation damage, and imaging at any angle. However, because the magnetic resonance imaging system is interfered by non-uniform magnetic field, coil sensitivity and the like, the reconstructed image shows certain non-uniformity, and certain images can be generated for doctor diagnosis and computer-aided analysis such as registration, classification and the like. With the demand for higher resolution scanning images, the magnetic field strength and the magnetic field gradient of the scanner are higher and finer, and the problem that the magnetic resonance images are subjected to interference of non-uniform fields is more and more serious. An inhomogeneous field is a deviant field caused by the inhomogeneity of the transmitted spatial magnetic field or the inhomogeneous sensitivity of the receiving coils, which is generally assumed to be a smooth, slowly varying multiplicative deviant field, resulting in a certain deviation between the grey values of the image and the true values. Some strong non-uniform fields can reduce the contrast of the image and overwhelm the lesion details, thereby leading to erroneous diagnosis or registration, segmentation errors. Therefore, correction of the inhomogeneous field is essential for each nuclear magnetic resonance image.

The non-uniform field correction method generally obtains a non-uniform field estimate by surface fitting the non-uniform characteristics of the image. The non-uniformity estimation based on gradient filtering utilizes the characteristic that low-frequency gradient of MRI is irrelevant to a tissue model, directly obtains second-order partial derivative of an image, performs low-pass filtering operation, and performs re-integration on the filtered gradient image, but filtering cannot completely provide gradient high-frequency information caused by image edges, and residual high-frequency gradient information can also cause large errors to the re-integration. And secondly, the fast non-uniform field estimation based on second-order gradient fitting is carried out, the characteristic that the low-frequency gradient of MRI in the first method is irrelevant to a tissue model is also utilized, a smaller value in a second derivative diagram of the image is selected as the estimation of the gradient caused by the non-uniform field, and the non-uniform field is obtained by directly fitting a low-order polynomial. The method has the defect that the gradient is selected by adopting a threshold value method, and the image details are easily and wrongly divided. If the non-uniform field of the image is small and the determined fitting data contains more image details, the fitting error of the second-order gradient is large, and the non-uniformity of the corrected image is more serious. Another non-uniform field correction based on row-by-row fitting assumes that the non-uniform field causes a slight change in the image gray scale, and the smaller value in the gradient field represents the image detail caused by the non-uniform field, and the larger value is the image detail. The small value of each line in the gradient field of the image is extracted to be used as the estimation of the gradient of the inhomogeneous field, the inhomogeneous field estimation in the x direction and the inhomogeneous field estimation in the y direction are obtained through one-dimensional polynomial fitting, and the inhomogeneous field of the image is obtained through multiplication of the inhomogeneous field estimation in the x direction and the inhomogeneous field estimation in the y direction. The defects are that the alignment between rows cannot be ensured, and enough fitting points cannot be extracted when the high-frequency characteristics of some rows account for most of the rows, so that the fitting error is overlarge. The invention still utilizes the characteristic that the low-frequency gradient in the image is irrelevant to the tissue model, improves the method and has the main characteristics that:

1) the invention can more accurately estimate the non-uniform field and correct the image.

2) The estimation method does not rely on other equipment and a priori knowledge.

3) The operation speed of the estimation method meets the real-time requirement.

Disclosure of Invention

The invention provides an MRI non-uniform field estimation method based on image gradient fitting, which is used for quickly correcting gray-scale non-uniformity caused by non-uniform magnetic field or coil sensitivity in a nuclear magnetic resonance image. The method firstly carries out logarithm solving on the image, converts the logarithm solving into a logarithm domain, and changes a multiplicative non-uniform field into an additive field, so that the image information and the non-uniform field information are decoupled. Then, a series of regions with slightly different gray values are selected, gradients in the x-direction and the y-direction are extracted, a special objective function is minimized through a least square method, fitting estimation of the non-uniform field of the full image is directly obtained, and the specific flow is as follows:

step one, denoising treatment

The nuclear magnetic resonance raw image is obtained by performing Fourier transform on the collected k-space signals. Subject to noise due to equipment and environmental influences. The invention is a gradient-based non-uniform field estimation, whereas the gradient is very sensitive to noise. The original image is denoised before estimation, and the preprocessing further comprises the steps of extracting the outline of the image, and removing the background area of the edge and some low signal-to-noise ratio areas.

Step two, calculating a gradient field

Estimating the inhomogeneous field requires a fitting operation on the gradient values within the region. The gradient field of the image, including the x-direction gradient and the y-direction gradient, is thus first calculated. The image gradient is generally obtained by adopting a difference method or a sobel operator, in order to reduce noise interference, the invention provides a novel gradient operator based on a Gaussian kernel to calculate the gradient of an image, and the influence of surrounding points is considered when the gradient of a certain point is calculated. Taking the x-direction as an example, assume that the original image is v, and v is obtained after the logarithm operation_logThe calculation formula is as follows:

△ therein_xFor differentiating sign, meaning differentiating the image in the x-direction, v_logIs a logarithmic domain image of an original image v, (x, y) represents the coordinates of pixel points of the image, m, n is the size of a gradient operator, and omega_i,jIs a Gaussian coefficient, satisfies

Step three, determining the region of interest

After the gradient field has been calculated, it can be used to determine a region of interest. By using a region growing method, firstly, the region of interest is initialized by determining a plurality of seed points, and whether the seed points are similar to the seeds is judged by judging the neighbors of the seed points to determine whether the seed points are added into the region of interest, wherein the method for determining the region of interest comprises the following steps:

firstly, selecting a seed region

The seed points constitute a seed region, the seed points belong to the same tissue region with a high probability, an indicator m (r) is defined to mark whether a certain voxel belongs to the seed region, when a certain voxel r belongs to the seed region, m (r) is 1, otherwise, m (r) is 0, and m (r) is determined by the following rule:

wherein I₀(r) represents the gray value of the voxel r in the initial image obtained from the image reconstructed from the previous iteration, the first iteration is obtained by simple averaging of the coil images, p is the peak value of the histogram of the composite image after the background region is removed, σ is the noise variance of the composite image, and the final seed region is represented asThis is the initialization region of the region growing algorithm.

Second, region growing

Once the seed regions are determined, the seed points can be used as initial regions of interest, similar points are continuously added to expand the regions of interest, the points in each region of interest are compared with the eight neighborhood points, if the gradient difference between the two points is less than a certain threshold value, the points are considered to be similar to the points in the region of interest, the points are added to the region of interest, if a boundary or a critical point with other tissues is encountered, the points are not added as similar points due to too large gradient values, the region of interest continuously grows until the point of interest does not change, and if the R-th region of interest is determined, the n-th region of interest is iterated to the nth time, and the n-th region of interest is updated as follows:

wherein,is the seed region for initialization, is the maximum gradient allowed, and its value is determined by experience value, and can be adjusted according to different images, and is set as half of the gradient value of all points in the seed region. Neigh (q) represents points in the eight neighborhood of point q, and Grad (p) is the gradient value of point p calculated by the gradient operatorThe stop zone grows when the number of newly added points no longer changes.

Step four, non-uniform field estimation

Assuming that the inhomogeneous field b is u, the real image is v, and the noise is n, the relationship between v (x) ═ b (x) u (x) + n (x), and the noise effect can be ignored after the denoising process. Carrying out logarithm operation on the relationship of the three to obtain v_log(x)＝b_log(x)+u_log(x) In that respect Because the gray values of the pixel points contained in the same region of interest are similar to each other, the edges and the details of the image are excluded, and therefore the gradient value of the point in the region of interest M can be assumed to be the gradient of the inhomogeneous field.

Assuming that the inhomogeneous field can be expressed by a polynomial of order K in 2 dimensions, the inhomogeneous field contains K ═ K +1) (K +2)/2 polynomial bases, each expressed as x^py^qWhere p + q is equal to or less than K, p is equal to or greater than 0, and q is equal to or greater than 0, let us note Fⁱ(x,y),0<i<K, assuming each coefficient as w_iK, then a non-uniform field in the log domainCan be expressed as:

assuming that there are N points in the determined region of interest M, respectively denoted as r₁,r₂,…,r_l，r_iIs represented by (x)_i,y_i) To estimate the parameter w_i1, K, consider the least squares equation that minimizes:

whereinAndare each a polynomial radical FⁱThe x-direction derivative and the y-direction derivative of △_xAnd △_yIs the differential of the image in the x and y directions, v_logIs a logarithmic domain image of the original image v, and (x, y) represents the pixel point coordinates of the image. To solve this least squares problem, we consider the following linear relationship:

△B_log＝△FW，

whereinThe solution of the above equation is as follows:

W＝(△F^t△F)^-1△F^t△B_log.

wherein △ F^tRepresenting the transpose of matrix △ F.

Step five, correcting

After the polynomial coefficients are obtained, the non-uniformity estimate of the full image can be obtained by extrapolation,

the corrected image is:

step six, iteration

The corrected image may have an improved inhomogeneous field, but the inhomogeneous field is not completely eliminated. And (3) using the corrected image for next round of correction, and iterating for 2 times according to experience to obtain a better correction effect.

The technical scheme of the invention has the advantages and positive effects

1) Non-uniform field estimation is carried out through the logarithm domain gradient, so that the fitting elements are independent of the tissue type, and errors caused by classification errors of tissue regions to the non-uniform field estimation are eliminated.

2) By designing a special objective function and simultaneously minimizing gradient errors in the x-direction and the y-direction, errors caused by the calculus are avoided, and estimation precision is improved.

The method is easy to realize, low in complexity and extremely low in time consumption in practical application.

Drawings

FIG. 1: algorithm flow chart of the invention

FIG. 2: the correction effects of the brain, the spine, the shoulders and the abdomen are compared, the first action is a non-uniform image, and the second action is a corrected image.

Detailed Description

Step one, denoising treatment

The nuclear magnetic resonance raw image is obtained by performing Fourier transform on the collected k-space signals. Subject to noise due to equipment and environmental influences. The invention is a gradient-based non-uniform field estimation, whereas the gradient is very sensitive to noise. The original image is denoised before estimation, and the preprocessing further comprises the steps of extracting the outline of the image, and removing the background area of the edge and some low signal-to-noise ratio areas.

Step two, calculating a gradient field

Estimating the inhomogeneous field requires a fitting operation on the gradient values within the region. The gradient field of the image, including the x-direction gradient and the y-direction gradient, is thus first calculated. The image gradient is generally obtained by adopting a difference method or a sobel operator, in order to reduce noise interference, the invention provides a novel gradient operator based on a Gaussian kernel to calculate the gradient of an image, and the influence of surrounding points is considered when the gradient of a certain point is calculated. Taking the x-direction as an example, assume that the original image is v, and v is obtained after the logarithm operation_logThe calculation formula is as follows:

△ therein_xFor differentiating sign, meaning differentiating the image in the x-direction, v_logIs a logarithmic domain image of an original image v, (x, y) represents the coordinates of pixel points of the image, m, n is the size of a gradient operator, and omega_i,jIs a Gaussian coefficient, satisfies

Step three, determining the region of interest

After the gradient field has been calculated, it can be used to determine a region of interest. By using a region growing method, firstly, the region of interest is initialized by determining a plurality of seed points, and whether the seed points are similar to the seeds is judged by judging the neighbors of the seed points to determine whether the seed points are added into the region of interest, wherein the method for determining the region of interest comprises the following steps:

firstly, selecting a seed region

The seed points constitute a seed region, the seed points belong to the same tissue region with a high probability, an indicator m (r) is defined to mark whether a certain voxel belongs to the seed region, when a certain voxel r belongs to the seed region, m (r) is 1, otherwise, m (r) is 0, and m (r) is determined by the following rule:

wherein I₀(r) represents the gray value of the voxel r in the initial image obtained from the image reconstructed from the previous iteration, the first iteration is obtained by simple averaging of the coil images, p is the peak value of the histogram of the composite image after the background region is removed, σ is the noise variance of the composite image, and the final seed region is represented asThis is the initialization region of the region growing algorithm.

Second, region growing

Once the seed regions are determined, the seed points can be used as initial regions of interest, similar points are continuously added to expand the regions of interest, the points in each region of interest are compared with the eight neighborhood points, if the gradient difference between the two points is less than a certain threshold value, the points are considered to be similar to the points in the region of interest, the points are added to the region of interest, if a boundary or a critical point with other tissues is encountered, the points are not added as similar points due to too large gradient values, the region of interest continuously grows until the point of interest does not change, and if the R-th region of interest is determined, the n-th region of interest is iterated to the nth time, and the n-th region of interest is updated as follows:

wherein,is the seed region for initialization, is the maximum gradient allowed, and its value is determined by experience value, and can be adjusted according to different images, and is set as half of the gradient value of all points in the seed region. Neigh (q) represents points in the eight neighborhood of point q, and Grad (p) is the gradient value of point p calculated by the gradient operatorThe stop zone grows when the number of newly added points no longer changes.

Step four, non-uniform field estimation

Assuming that the inhomogeneous field b is u, the real image is v, and the noise is n, the relationship between v (x) ═ b (x) u (x) + n (x), and the noise effect can be ignored after the denoising process. Carrying out logarithm operation on the relationship of the three to obtain v_log(x)＝b_log(x)+u_log(x) In that respect Because the gray values of the pixel points contained in the same region of interest are similar to each other, the edges and the details of the image are excluded, and therefore the gradient value of the point in the region of interest M can be assumed to be the gradient of the inhomogeneous field.

Assuming that the inhomogeneous field can be expressed by a polynomial of order K in 2 dimensions, the inhomogeneous field contains K ═ K +1) (K +2)/2 polynomial bases, each expressed as x^py^qWhere p + q is equal to or less than K, p is equal to or greater than 0, and q is equal to or greater than 0, let us note Fⁱ(x,y),0<i<K, assuming each coefficient as w_iK, then a non-uniform field in the log domainCan be expressed as:

assuming that there are N points in the determined region of interest M, respectively denoted as r₁,r₂,...,r_N，r_iIs represented by (x)_i,y_i) To estimate the parameter w_i1, K, consider the least squares equation that minimizes:

whereinAndare each a polynomial radical FⁱThe x-direction derivative and the y-direction derivative of △_xAnd △_yIs the differential of the image in the x and y directions, v_logIs a logarithmic domain image of the original image v, and (x, y) represents the pixel point coordinates of the image. To solve this least squares problem, we consider the following linear relationship:

△B_log＝△FW，

whereinThe solution of the above equation is as follows:

W＝(△F^t△F)^-1△F^t△B_log.

wherein △ F^tRepresenting the transpose of matrix △ F.

Step five, correcting

After the polynomial coefficients are obtained, the non-uniformity estimate of the full image can be obtained by extrapolation,

the corrected image is:

step six, iteration

The corrected image may have an improved inhomogeneous field, but the inhomogeneous field is not completely eliminated. And (3) using the corrected image for next round of correction, and iterating for 2 times according to experience to obtain a better correction effect.

Claims (2)

1. An MRI non-uniform field estimation method based on image gradient fitting is used for quickly correcting gray-scale non-uniformity caused by magnetic field non-uniformity or coil sensitivity in a nuclear magnetic resonance image; the method is characterized by comprising the following steps:

step one, denoising treatment

The nuclear magnetic resonance original image is obtained by carrying out Fourier transform on collected k-space signals; carrying out denoising processing on an original image before estimation, wherein the preprocessing also comprises the outline extraction of the image, and the background area of the edge and a plurality of low signal-to-noise ratio areas are removed;

step two, calculating a gradient field

Estimating the inhomogeneous field requires fitting the gradient values in the region; thus, firstly, calculating the gradient field of the image, including the gradient in the x-direction and the gradient in the y-direction; in order to reduce noise interference, a gradient operator based on a Gaussian kernel is adopted to calculate the gradient of the image, and the influence of surrounding points is considered when the gradient of a certain point is calculated; for the x-direction, assume the original image is v, and after the logarithmic operation is v_logThe calculation formula is as follows:

${\mathrm{\&Delta;}}_x{v}_{\log }=\underset{j=-n}{\overset{n}{\&Sigma;}}\underset{i=1}{\overset{m}{\&Sigma;}}{\mathrm{\&omega;}}_{i,j}\&lsqb;v_{\log }(x+i,y+j)-v_{\log }(x-i,y+j)\&rsqb;$

△ therein_xFor differentiating sign, meaning differentiating the image in the x-direction, v_logIs a logarithmic domain image of an original image v, (x, y) represents the coordinates of pixel points of the image, m, n is the size of a gradient operator, and omega_i，jIs a Gaussian coefficient, satisfies

Step three, determining the region of interest

After the gradient field is calculated, the region of interest is determined by using the gradient field; firstly, initializing an interested region by determining a plurality of seed points by using a region growing method, and determining whether the interested region is added into the interested region by judging whether the neighbors of the seed points are similar to the seeds;

step four, non-uniform field estimation

Assuming that the inhomogeneous field is b, the real image is u, the original image is v, and the noise is n, the relationship between v (x) and b (x) u (x) and n (x) is n, and the noise influence is ignored after the denoising process is performed; carrying out logarithm operation on the relationship of the three to obtain v_log(x)＝b_log(x)+u_log(x) (ii) a Because the gray values of the pixel points contained in the same region of interest are similar to each other, the edge and the details of the image are eliminated, and therefore the gradient value of the point in the region of interest M is assumed to be the gradient of the non-uniform field;

assuming that the inhomogeneous field is expressed by a polynomial of order K in 2 dimensions, the inhomogeneous field contains K ═ K +1) (K +2)/2 polynomial bases, each expressed as x^cy^dWhere c + d is equal to or less than K, c is equal to or greater than 0, and d is equal to or greater than 0, denoted as Fⁱ(x,y),0<i is less than or equal to K, and each coefficient is assumed to be w_iK, then a non-uniform field in the log domainExpressed as:

${\hat{b}}_{\log }(x,y)=\underset{i=1}{\overset{K}{\&Sigma;}}{w}_i{F}^i(x,y)$

assuming that there are N points in the determined region of interest M, respectively denoted as r₁,r₂,...,r_N，r_lIs represented by (x)_l,y_l) To estimate the parameter w_i1, K, consider the least squares equation that minimizes:

whereinAndare each a polynomial radical FⁱThe x-direction derivative and the y-direction derivative of △_xAnd △_yDifferential, v, in the x-and y-directions of the image, respectively_logIs a logarithmic domain image of an original image v, and (x, y) represents the pixel point coordinates of the image; to solve the least squares problem, consider the following linear relationship:

△B_log＝△FW，

wherein

The solution of the above equation is as follows:

W＝(△F^t△F)^-1△F^t△B_log

wherein △ F^tRepresents a transpose of matrix △ F;

step five, correcting

After obtaining the polynomial coefficients, the non-uniform estimation of the whole image is obtained through extrapolation,

${\hat{b}}_{\log }(x,y)=\left[\begin{array}{cccc}{F}^1(x,y)& F^2(x,y)& \mathrm{...}& F^K(x,y)\end{array}\right]W,$

the corrected image is:

$\hat{u}=e^{v_{\log }-{\hat{b}}_{\log }}$

step six, iteration

Although the non-uniform field of the corrected image is improved, the non-uniform field is not completely eliminated, and at the moment, the corrected image is used for carrying out the next correction, and the better correction effect is obtained after 2 times of iteration.

2. The MRI inhomogeneous field estimation method based on image gradient fitting according to claim 1, characterized by: the determination method of the region of interest is as follows:

firstly, selecting a seed region

The seed points constitute a seed region, the seed points belong to the same tissue region with a high probability, an indicator m (r) is defined to mark whether a certain voxel belongs to the seed region, when a certain voxel r belongs to the seed region, m (r) is 1, otherwise, m (r) is 0, and m (r) is determined by the following rule:

$M\left(r\right)=\left\{\begin{array}{cc}1,& a-\mathrm{\&sigma;}<I_0\left(r\right)<a+\mathrm{\&sigma;}\\ 0,& other\end{array}\right.$

wherein I₀(r) represents the gray value of the voxel r in the initial image obtained from the image reconstructed from the previous iteration, the first iteration is obtained by simple averaging of the coil images, a is the peak value of the histogram of the composite image after the background region is removed, σ is the noise variance of the composite image, and the final seed region is represented asThis is the initialization region of the region growing algorithm;

second, region growing

Once the seed region is determined, using the seed point as the initial region of interest, continuously adding similar points to expand the region of interest, comparing the point in each region of interest with its eight neighboring points, if the gradient difference between the two is less than a certain threshold, considering the eight neighboring points as being similar to the points in the region of interest, adding them into the region of interest, if a boundary or a critical point with other tissues is encountered, adding them as similar points due to too large gradient value, and continuously iterating to increase the region of interest until it does not change any more, assuming that the R th region of interest is determined, iterating to the f th time, then the f th region of interest is updated as follows:

$M_R^{\left(f\right)}=M_R^{(f-1)}\&cup;\left\{p\right|p\&Element;Neigh\left(q\right),Grad\left(p\right)<\mathrm{\&delta;},q\&Element;M_R^{(f-1)}\}$

wherein,the seed area is used for initialization, the maximum gradient is allowed, the value of the maximum gradient is determined by empirical values, the adjustment is carried out according to different images, and the maximum gradient is set to be half of the gradient value of all points in the seed area; neigh (q) represents points in the eight neighborhood of point q, and Grad (p) is the gradient value of point p calculated by the gradient operatorThe stop zone grows when the number of newly added points no longer changes.