CN104392422B - A MRI Inhomogeneous Field Estimation Method 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

A MRI Inhomogeneous Field Estimation Method Based on Image Gradient Fitting

technical field

Quickly estimate and correct for inhomogeneous fields in MRI images due to magnetic inhomogeneity or coil sensitivity.

Background technique

Magnetic resonance imaging is widely used in medical diagnosis and scientific research because of its high resolution, no ionizing radiation damage, and imaging at any angle. However, due to the interference of magnetic field inhomogeneity and coil sensitivity, etc., the reconstructed image of the MRI system shows certain inhomogeneity, which will produce certain images for doctor diagnosis and computer-aided analysis such as registration and classification. With the demand for higher-resolution scanning images, the magnetic field strength of the scanner is getting higher and higher, and the magnetic field gradient is getting finer. It's getting serious. The inhomogeneous field is the deviation field caused by the inhomogeneity of the transmitted spatial magnetic field or the inhomogeneity of the receiving coil sensitivity. This deviation field is generally assumed to be a smooth, slowly changing multiplicative deviation field, which will lead to gray images. There is a certain deviation between the degree value and the real value. Some strong non-uniform fields will reduce the contrast of the image and drown out the details of the lesion, which will lead to wrong diagnosis results or registration and segmentation errors. Therefore, correction of inhomogeneous fields is essential for every MRI image.

The non-uniform field correction method generally obtains the non-uniform field estimation by performing surface fitting on the non-uniform characteristics of the image. A non-uniformity estimation based on gradient filtering, which uses the characteristics that the low-frequency gradient of MRI has nothing to do with the tissue model, directly obtains the second-order partial derivative of the image and performs low-pass filtering operation, and re-integrates the filtered gradient image, but Filtering cannot fully extract the gradient high-frequency information caused by the edge of the image, and the remaining high-frequency gradient information will also cause large errors in the reintegration. The second method is a fast non-uniform field estimation based on second-order gradient fitting. It also uses the characteristics that the low-frequency gradient of MRI has nothing to do with the tissue model in the first method, and selects the smaller value in the second-order derivative map of the image as the non-uniform field. The gradient caused by the field is estimated, and the non-uniform field is obtained by direct fitting of a low-order polynomial. Its disadvantage is that the gradient is selected by the threshold method, and it is easy to mistakenly divide the image details into it. If the non-uniform field of the image is small, and the determined fitting data contains more image details, the resulting second-order gradient fitting error will be larger, and the non-uniformity of the corrected image may be more serious. Another non-uniform field correction based on row-by-column fitting assumes that the non-uniform field causes a small change in image grayscale, then the smaller value in the gradient field represents the cause of the non-uniform field, and the larger value is the image detail . By extracting the smaller value of each row and column in the gradient field of the image, as an estimate of the gradient of the non-uniform field, the non-uniform field estimation in the x direction and the y direction is obtained by one-dimensional polynomial fitting, and the image is obtained by multiplying the two non-uniform field. The disadvantage is that the alignment between rows cannot be guaranteed, and for some rows whose high-frequency characteristics account for the majority, it is impossible to extract enough fitting points, resulting in excessive fitting errors. The present invention still utilizes the characteristic that the low-frequency gradient in the image is irrelevant to the tissue model, improves the above-mentioned method, and has the following main features:

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

2) The estimation method does not depend on other equipment and prior knowledge.

3) The calculation speed of the estimation method meets the real-time requirements.

Contents of the invention

The invention proposes an MRI inhomogeneous field estimation method based on image gradient fitting, which is used for fast correction of the gray level inhomogeneity in the nuclear magnetic resonance image caused by the inhomogeneity of the magnetic field or the sensitivity of the coil. This method first calculates the logarithm of the image, converts it to the logarithmic domain, and turns the multiplicative non-uniform field into an additive, thus decoupling the image information and the non-uniform field information. Then select a series of areas with little difference in gray value, extract the gradients in the x-direction and y-direction, and minimize a special objective function through the least square method to directly obtain the fitting estimation of the non-uniform field of the whole image. The specific process is as follows:

Step 1, denoising processing

The original NMR image is obtained by Fourier transform of the collected k-space signal. Noise due to equipment and environment. The present invention is gradient-based non-uniform field estimation, and the gradient is very sensitive to noise. Before the estimation, the original image should be denoised. The preprocessing also includes the contour extraction of the image, removing the background area of the edge and some low signal-to-noise ratio areas.

Step 2. Calculate the gradient field

Estimating an inhomogeneous field requires a fitting operation on the gradient values in the region. Therefore, the gradient field of the image is first calculated, including the x-direction gradient and the y-direction gradient. The image gradient is generally obtained by difference method or sobel operator. In order to reduce noise interference, the present invention proposes a new gradient operator based on Gaussian kernel to calculate the gradient of the image. When calculating the gradient of a certain point, the surrounding points are considered influences. Taking the x-direction as an example, assuming that the original image is v, after the logarithmic operation is v _log , the calculation formula is:

Among them, △ _x is the differential symbol, which means to differentiate the image in the x direction, v _log is the logarithmic domain image of the original image v, (x, y) represents the pixel coordinates of the image, m, n are the size of the gradient operator, ω _i,j are Gaussian coefficients, satisfying

Step 3. Determination of the region of interest

Once the gradient field has been calculated, it can be used to identify regions of interest. Using the region growing method, the region of interest is first initialized by determining some seed points, and whether the neighbors of these seed points are similar to the seeds is used to determine whether to add it to the region of interest. The method of determining the region of interest is as follows:

First, select the seed region

The seed points form the seed area, and the seed points belong to the same tissue area with a greater probability. An indicator M(r) is defined to mark whether a voxel belongs to the seed area. When a voxel r belongs to the seed area, M(r )=1, otherwise, M(r)=0, M(r) is determined by the following rules:

Among them, I ₀ (r) represents the gray value of voxel r in the initial image, the initial image is obtained from the image reconstructed in the last iteration, and the first iteration is obtained by simple average of each coil image, p is the removal of The peak after the background area, σ is the noise variance of the synthetic image, and the final seed area is expressed as This is the initialization region for the region growing algorithm.

Second, regional growth

Once the seed region is determined, the seed point can be used as the initial region of interest, and similar points can be continuously added to expand the region of interest, by comparing the points in each region of interest with its eight neighbor points , if the gradient difference between the two is less than a certain threshold, it is considered that this point is similar to the point in the region of interest, and it will be added to the region of interest. If it encounters a boundary or a critical point with other tissues, due to the gradient If the value is too large, it will not be added as a similar point. Through continuous iteration, the region of interest will continue to grow until it no longer changes. Assuming that the Rth region of interest is determined and iterated to the nth time, then the nth The region of interest is updated as follows:

in, Is the seed area used for initialization, δ is the maximum gradient allowed, its value is determined by experience, can be adjusted according to different images, set to half of the gradient value of all points in the seed area. Neigh(q) represents the point in the eight-neighborhood of point q, and Grad(p) is the gradient value of point p calculated using the above gradient operator Stop growing the region when the number of newly added points no longer changes.

Step 4. Inhomogeneous field estimation

Assuming a non-uniform field b, the real image is u, the original image is v, and the noise is n, then their relationship is v(x)=b(x)u(x)+n(x), after denoising processing can be Noise effects are ignored. Logarithmic operation is performed on the relationship among the three, and v _log (x)=b _log (x)+u _log (x) is obtained. Since the gray values of the pixels contained in the same region of interest are similar to each other, the edges and details of the image are excluded, so the gradient value of the points in the region of interest M can be assumed to be the gradient of the non-uniform field.

Assuming that the non-uniform field can be represented by a 2-dimensional k-order polynomial, it contains K=(k+1)(k+2)/2 polynomial bases, which are expressed as x ^p y ^q , where p+q≤K , p≥0, q≥0, we record it as F ⁱ (x,y), 0<i<K, assuming that the coefficients are w _i , i=1,...,K, then the non uniform field It can be expressed as:

Assume that there are N points in the determined region of interest M, which are respectively denoted as r ₁ , r ₂ ,..., r _l , and _ri represents ( _xi , y _i ), in order to estimate the parameter w _i , i=1, ...,K, consider minimizing the following least squares equation:

in with are the x-direction derivative and y-direction derivative of the polynomial basis F ⁱ respectively, △ _x and △ _y are the differentials in the x-direction and y-direction of the image, v _log is the logarithmic domain image of the original image v, (x, y) represents The pixel coordinates of the image. To solve this least squares problem, we consider the following linear relationship:

ΔB _log = ΔFW,

in The solution to the above formula is as follows:

W＝(△F ^t △F) ^-1 △F ^t △B _log .

where ΔF ^t represents the transpose of matrix ΔF.

Step 5. Calibration

After obtaining the polynomial coefficients, the non-uniform estimation of the full image can be obtained by extrapolation,

The corrected image is:

Step 6. Iteration

Although the non-uniform field of the corrected image has been improved, sometimes the non-uniform field has not been completely eliminated. The corrected image is used for the next round of correction. According to experience, better correction effect can be obtained after two iterations.

Advantages and positive effects of the technical solution of the present invention

1) By estimating the non-uniform field on the gradient of the number domain, the fitting elements are independent of the tissue type, and the error caused by the classification error of the tissue region on the non-uniform field estimation is eliminated.

2) By designing a special objective function to minimize the gradient errors in the x-direction and y-direction at the same time, the error caused by the heavy integration is avoided and the estimation accuracy is improved.

The method is easy to implement, has low complexity, and is extremely time-consuming in practical application.

Description of drawings

Figure 1: Algorithm flow chart of the present invention

Figure 2: Comparison of the correction effects of the brain, spine, shoulders and abdomen, the first row is the non-uniform image, and the second row is the corrected image.

detailed description

Step 1, denoising processing

The original NMR image is obtained by Fourier transform of the collected k-space signal. Noise due to equipment and environment. The present invention is gradient-based non-uniform field estimation, and the gradient is very sensitive to noise. Before the estimation, the original image should be denoised. The preprocessing also includes the contour extraction of the image, removing the background area of the edge and some low signal-to-noise ratio areas.

Step 2. Calculate the gradient field

Estimating an inhomogeneous field requires a fitting operation on the gradient values in the region. Therefore, the gradient field of the image is first calculated, including the x-direction gradient and the y-direction gradient. The image gradient is generally obtained by difference method or sobel operator. In order to reduce noise interference, the present invention proposes a new gradient operator based on Gaussian kernel to calculate the gradient of the image. When calculating the gradient of a certain point, the surrounding points are considered influences. Taking the x-direction as an example, assuming that the original image is v, after the logarithmic operation is v _log , the calculation formula is:

Among them, △ _x is the differential symbol, which means to differentiate the image in the x direction, v _log is the logarithmic domain image of the original image v, (x, y) represents the pixel coordinates of the image, m, n are the size of the gradient operator, ω _i,j are Gaussian coefficients, satisfying

Step 3. Determination of the region of interest

Once the gradient field has been calculated, it can be used to identify regions of interest. Using the region growing method, the region of interest is first initialized by determining some seed points, and whether the neighbors of these seed points are similar to the seeds is used to determine whether to add them to the region of interest. The method of determining the region of interest is as follows:

First, select the seed region

The seed points form the seed area, and the seed points belong to the same tissue area with a greater probability. An indicator M(r) is defined to mark whether a voxel belongs to the seed area. When a voxel r belongs to the seed area, M(r )=1, otherwise, M(r)=0, M(r) is determined by the following rules:

Among them, I ₀ (r) represents the gray value of voxel r in the initial image, the initial image is obtained from the image reconstructed in the last iteration, and the first iteration is obtained by simple average of each coil image, p is the removal of The peak after the background area, σ is the noise variance of the synthetic image, and the final seed area is expressed as This is the initialization region for the region growing algorithm.

Second, regional growth

Once the seed region is determined, the seed point can be used as the initial region of interest, and similar points can be continuously added to expand the region of interest, by comparing the points in each region of interest with its eight neighbor points , if the gradient difference between the two is less than a certain threshold, it is considered that this point is similar to the point in the region of interest, and it will be added to the region of interest. If it encounters a boundary or a critical point with other tissues, due to the gradient If the value is too large, it will not be added as a similar point. Through continuous iteration, the region of interest will continue to grow until it no longer changes. Assuming that the Rth region of interest is determined and iterated to the nth time, then the nth The region of interest is updated as follows:

in, Is the seed area used for initialization, δ is the maximum gradient allowed, its value is determined by experience, can be adjusted according to different images, set to half of the gradient value of all points in the seed area. Neigh(q) represents the point in the eight-neighborhood of point q, and Grad(p) is the gradient value of point p calculated using the above gradient operator Stop growing the region when the number of newly added points no longer changes.

Step 4. Inhomogeneous field estimation

Assuming a non-uniform field b, the real image is u, the original image is v, and the noise is n, then their relationship is v(x)=b(x)u(x)+n(x), after denoising processing can be Noise effects are ignored. Logarithmic operation is performed on the relationship among the three, and v _log (x)=b _log (x)+u _log (x) is obtained. Since the gray values of the pixels contained in the same region of interest are similar to each other, the edges and details of the image are excluded, so the gradient value of the points in the region of interest M can be assumed to be the gradient of the non-uniform field.

Assuming that the non-uniform field can be represented by a 2-dimensional k-order polynomial, it contains K=(k+1)(k+2)/2 polynomial bases, which are expressed as x ^p y ^q , where p+q≤K , p≥0, q≥0, we record it as F ⁱ (x,y), 0<i<K, assuming that the coefficients are w _i , i=1,...,K, then the non uniform field It can be expressed as:

Assuming that there are a total of N points in the determined region of interest M, denoted as r ₁ , r ₂ ,..., r _N , r _i represents ( _xi , y _i ), in order to estimate the parameter w _i , i= 1,...,K, consider minimizing the following least squares equation:

in with are the x-direction derivative and y-direction derivative of the polynomial basis F ⁱ respectively, △ _x and △ _y are the differentials in the x-direction and y-direction of the image, v _log is the logarithmic domain image of the original image v, (x, y) represents The pixel coordinates of the image. To solve this least squares problem, we consider the following linear relationship:

ΔB _log = ΔFW,

in The solution to the above formula is as follows:

W＝(△F ^t △F) ^-1 △F ^t △B _log .

where ΔF ^t represents the transpose of matrix ΔF.

Step 5. Calibration

After obtaining the polynomial coefficients, the non-uniform estimation of the full image can be obtained by extrapolation,

The corrected image is:

Step 6. Iteration

Although the non-uniform field of the corrected image has been improved, sometimes the non-uniform field has not been completely eliminated. The corrected image is used for the next round of correction. According to experience, better correction effect can be obtained after two iterations.

Claims (2)

1. An MRI inhomogeneous field estimation method based on image gradient fitting is used to quickly correct the gray level inhomogeneity caused by magnetic field inhomogeneity or coil sensitivity in the nuclear magnetic resonance image; it is characterized in that comprising the steps: Step 1, denoising processing The original NMR image is obtained by Fourier transform of the collected k-space signal; the original image is denoised before estimation, and the preprocessing also includes image contour extraction, removing edge background areas and some low signal-to-noise ratio areas; Step 2. Calculate the gradient field Estimating the inhomogeneous field requires a fitting operation on the gradient value in the region; therefore, first calculate 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 algorithm based on a Gaussian kernel is used When calculating the gradient of a certain point, consider the influence of its surrounding points; for the x-direction, assuming that the original image is v, after the logarithmic operation is v _log , the calculation formula is: ${\mathrm{<span\; class="notranslate">\; \&Delta;</span>}}_{\mathrm{<span\; class="notranslate">\; x</span>}}{\mathrm{<span\; class="notranslate">\; v</span>}}_{\mathrm{<span\; class="notranslate">\; log</span>}}<span\; class="notranslate">\; =</span>\underset{\mathrm{<span\; class="notranslate">\; j</span>}<span\; class="notranslate">\; =</span><span\; class="notranslate">\; -</span>\mathrm{<span\; class="notranslate">\; no</span>}}{\overset{\mathrm{<span\; class="notranslate">\; no</span>}}{<span\; class="notranslate">\; \&Sigma;</span>}}\underset{\mathrm{<span\; class="notranslate">\; i</span>}<span\; class="notranslate">\; =</span>\mathrm{<span\; class="notranslate">\; 1</span>}}{\overset{\mathrm{<span\; class="notranslate">\; m</span>}}{<span\; class="notranslate">\; \&Sigma;</span>}}{\mathrm{<span\; class="notranslate">\; \&omega;</span>}}_{\mathrm{<span\; class="notranslate">\; i</span>}<span\; class="notranslate">\; ,</span>\mathrm{<span\; class="notranslate">\; j</span>}}<span\; class="notranslate">\; \&lsqb;</span>{\mathrm{<span\; class="notranslate">\; v</span>}}_{\mathrm{<span\; class="notranslate">\; log</span>}}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; x</span>}<span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; i</span>}<span\; class="notranslate">\; ,</span>\mathrm{<span\; class="notranslate">\; the\; y</span>}<span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; j</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; -</span>{\mathrm{<span\; class="notranslate">\; v</span>}}_{\mathrm{<span\; class="notranslate">\; log</span>}}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; x</span>}<span\; class="notranslate">\; -</span>\mathrm{<span\; class="notranslate">\; i</span>}<span\; class="notranslate">\; ,</span>\mathrm{<span\; class="notranslate">\; the\; y</span>}<span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; j</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; \&rsqb;</span>$ Among them, △ _x is the differential symbol, which means to differentiate the image in the x direction, v _log is the logarithmic domain image of the original image v, (x, y) represents the pixel coordinates of the image, m, n are the size of the gradient operator, ω _{i, j} are Gaussian coefficients, satisfying Step 3. Determination of the region of interest After the gradient field is calculated, it is used to determine the region of interest; using the region growing method, first the region of interest is initialized by determining some seed points, and whether to add it to the sensor by judging whether the neighbors of these seed points are similar to the seed area of interest; Step 4. Inhomogeneous field estimation Suppose the non-uniform field is b, the real image is u, the original image is v, and the noise is n, then their relationship is v(x)=b(x)u(x)+n(x), after denoising Ignore the influence of noise; perform logarithmic operation on the relationship between the three, and get v _log (x)=b _log (x)+u _log (x); because the gray values of the pixels contained in the same region of interest are similar to each other, the images are excluded The edges and details of , so the gradient value of the points in the region of interest M is assumed to be the gradient of the non-uniform field; Assuming that the non-uniform field is represented by a 2-dimensional k-order polynomial, it contains a total of K=(k+1)(k+2)/2 polynomial bases, expressed as x ^c y ^d , where c+d≤K, c≥0,d≥0, denoted as F ⁱ (x,y), 0<i≤K, assuming that the coefficients are w _i , i=1,...,K, then the non-uniform field in the logarithmic domain Expressed as: ${\stackrel{<span\; class="notranslate">\; ^</span>}{\mathrm{<span\; class="notranslate">\; b</span>}}}_{\mathrm{<span\; class="notranslate">\; log</span>}}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; x</span>}<span\; class="notranslate">\; ,</span>\mathrm{<span\; class="notranslate">\; the\; y</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; =</span>\underset{\mathrm{<span\; class="notranslate">\; i</span>}<span\; class="notranslate">\; =</span>\mathrm{<span\; class="notranslate">\; 1</span>}}{\overset{\mathrm{<span\; class="notranslate">\; K</span>}}{<span\; class="notranslate">\; \&Sigma;</span>}}{\mathrm{<span\; class="notranslate">\; w</span>}}_{\mathrm{<span\; class="notranslate">\; i</span>}}{\mathrm{<span\; class="notranslate">\; f</span>}}^{\mathrm{<span\; class="notranslate">\; i</span>}}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; x</span>}<span\; class="notranslate">\; ,</span>\mathrm{<span\; class="notranslate">\; the\; y</span>}<span\; class="notranslate">\; )</span>$ Assuming that there are N points in the determined region of interest M, which are recorded as r ₁ , r ₂ ,..., r _N , r _l represents (x _l , y _l ), in order to estimate the parameter w _i , i= 1,...,K, consider minimizing the following least squares equation: in with are the x-direction derivative and y-direction derivative of the polynomial basis F ⁱ respectively, △ _x and △ _y are the differentials of the image in the x-direction and y-direction respectively, v _log is the logarithmic domain image of the original image v, (x, y) Represents the pixel coordinates of the image; in order to solve the least squares problem, consider the following linear relationship: ΔB _log = ΔFW, in The solution to the above formula is as follows: W＝(△F ^t △F) ^-1 △F ^t △B _log Where △F ^t represents the transpose of matrix △F; Step 5. Calibration After obtaining the polynomial coefficients, the non-uniform estimation of the full image is obtained by extrapolation, ${\stackrel{<span\; class="notranslate">\; ^</span>}{\mathrm{<span\; class="notranslate">\; b</span>}}}_{\mathrm{<span\; class="notranslate">\; log</span>}}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; x</span>}<span\; class="notranslate">\; ,</span>\mathrm{<span\; class="notranslate">\; the\; y</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; =</span>\left[\begin{array}{cccc}{\mathrm{<span\; class="notranslate">\; f</span>}}^{\mathrm{<span\; class="notranslate">\; 1</span>}}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; x</span>}<span\; class="notranslate">\; ,</span>\mathrm{<span\; class="notranslate">\; the\; y</span>}<span\; class="notranslate">\; )</span>& {\mathrm{<span\; class="notranslate">\; f</span>}}^{\mathrm{<span\; class="notranslate">\; 2</span>}}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; x</span>}<span\; class="notranslate">\; ,</span>\mathrm{<span\; class="notranslate">\; the\; y</span>}<span\; class="notranslate">\; )</span>& \mathrm{<span\; class="notranslate">\; ...</span>}& {\mathrm{<span\; class="notranslate">\; f</span>}}^{\mathrm{<span\; class="notranslate">\; K</span>}}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; x</span>}<span\; class="notranslate">\; ,</span>\mathrm{<span\; class="notranslate">\; the\; y</span>}<span\; class="notranslate">\; )</span>\end{array}\right]\mathrm{<span\; class="notranslate">\; W</span>}<span\; class="notranslate">\; ,</span>$ The corrected image is: $\stackrel{<span\; class="notranslate">\; ^</span>}{\mathrm{<span\; class="notranslate">\; u</span>}}<span\; class="notranslate">\; =</span>{\mathrm{<span\; class="notranslate">\; e</span>}}^{{\mathrm{<span\; class="notranslate">\; v</span>}}_{\mathrm{<span\; class="notranslate">\; log</span>}}<span\; class="notranslate">\; -</span>{\stackrel{<span\; class="notranslate">\; ^</span>}{\mathrm{<span\; class="notranslate">\; b</span>}}}_{\mathrm{<span\; class="notranslate">\; log</span>}}}$ Step 6. Iteration Although the non-uniform field of the corrected image has improved, there is still a situation where the non-uniform field has not been completely eliminated. At this time, the corrected image is used for the next round of correction, and iterated twice to obtain a better correction effect . 2. a kind of MRI inhomogeneous field estimation method based on image gradient fitting according to claim 1, is characterized in that: the determination method of region of interest is as follows: First, select the seed region The seed points form the seed area, and the seed points belong to the same tissue area with a high probability. An indicator M(r) is defined to mark whether a voxel belongs to the seed area. When a voxel r belongs to the seed area, M(r )=1, otherwise, M(r)=0, M(r) is determined by the following rules: $\mathrm{<span\; class="notranslate">\; m</span>}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; r</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; =</span>\left\{\begin{array}{cc}\mathrm{<span\; class="notranslate">\; 1</span>}<span\; class="notranslate">\; ,</span>& \mathrm{<span\; class="notranslate">\; a</span>}<span\; class="notranslate">\; -</span>\mathrm{<span\; class="notranslate">\; \&sigma;</span>}<span\; class="notranslate">\; <</span>{\mathrm{<span\; class="notranslate">\; I</span>}}_{\mathrm{<span\; class="notranslate">\; 0</span>}}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; r</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; <</span>\mathrm{<span\; class="notranslate">\; a</span>}<span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; \&sigma;</span>}\\ \mathrm{<span\; class="notranslate">\; 0</span>}<span\; class="notranslate">\; ,</span>& \mathrm{<span\; class="notranslate">\; o</span>}\mathrm{<span\; class="notranslate">\; t</span>}\mathrm{<span\; class="notranslate">\; h</span>}\mathrm{<span\; class="notranslate">\; e</span>}\mathrm{<span\; class="notranslate">\; r</span>}\end{array}\right.$ Among them, I ₀ (r) represents the gray value of voxel r in the initial image. The initial image is obtained from the image reconstructed in the last iteration. The first iteration is obtained by simple average of each coil image. a is the histogram of the composite image. The peak after the background area, σ is the noise variance of the synthetic image, and the final seed area is expressed as This is the initialization region for the region growing algorithm; Second, regional growth Once the seed region is determined, the seed point is used as the initial region of interest, and similar points are continuously added to expand the region of interest. By comparing the points in each region of interest with its eight neighbor points, If the gradient difference between the two is less than a certain threshold, it is considered that the eight neighbor points are similar to the points in the region of interest, and it is added to the region of interest. If it encounters a boundary or a critical point with other tissues, then Because the gradient value is too large, it will not be added as a similar point. Through continuous iterations, the region of interest will continue to grow until it no longer changes. Assuming that the Rth region of interest is determined and iterated to the fth time, then The region of interest for the fth time is updated as follows: ${\mathrm{<span\; class="notranslate">\; m</span>}}_{\mathrm{<span\; class="notranslate">\; R</span>}}^{<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; f</span>}<span\; class="notranslate">\; )</span>}<span\; class="notranslate">\; =</span>{\mathrm{<span\; class="notranslate">\; m</span>}}_{\mathrm{<span\; class="notranslate">\; R</span>}}^{<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; f</span>}<span\; class="notranslate">\; -</span>\mathrm{<span\; class="notranslate">\; 1</span>}<span\; class="notranslate">\; )</span>}<span\; class="notranslate">\; \&cup;</span><span\; class="notranslate">\; \{</span>\mathrm{<span\; class="notranslate">\; p</span>}<span\; class="notranslate">\; |</span>\mathrm{<span\; class="notranslate">\; p</span>}<span\; class="notranslate">\; \&Element;</span>\mathrm{<span\; class="notranslate">\; N</span>}\mathrm{<span\; class="notranslate">\; e</span>}\mathrm{<span\; class="notranslate">\; i</span>}\mathrm{<span\; class="notranslate">\; g</span>}\mathrm{<span\; class="notranslate">\; h</span>}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; q</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; ,</span>\mathrm{<span\; class="notranslate">\; G</span>}\mathrm{<span\; class="notranslate">\; r</span>}\mathrm{<span\; class="notranslate">\; a</span>}\mathrm{<span\; class="notranslate">\; d</span>}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; p</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; <</span>\mathrm{<span\; class="notranslate">\; \&delta;</span>}<span\; class="notranslate">\; ,</span>\mathrm{<span\; class="notranslate">\; q</span>}<span\; class="notranslate">\; \&Element;</span>{\mathrm{<span\; class="notranslate">\; m</span>}}_{\mathrm{<span\; class="notranslate">\; R</span>}}^{<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; f</span>}<span\; class="notranslate">\; -</span>\mathrm{<span\; class="notranslate">\; 1</span>}<span\; class="notranslate">\; )</span>}<span\; class="notranslate">\; \}</span>$ in, is the seed area used for initialization, δ is the maximum gradient allowed, its value is determined by the experience value, adjusted according to different images, set to half of the gradient value of all points in the seed area; Neigh(q) represents the point q Points in the eight-neighborhood, Grad(p) is the gradient value of point p calculated using the above gradient operator Stop growing the region when the number of newly added points no longer changes.