CN103632367B - A kind of MRI coil sensitivity estimation method based on the matching of many tissue regions - Google Patents

Abstract

The invention provides a kind of MRI coil sensitivity estimation method based on the matching of many tissue regions, the method key point is: 1) region growth method determination area-of-interest, and the area-of-interest determined by region growth method can ensure abundant point and enough accurate gray scale difference goes to estimate sensitivity function; 2) determination of many tissue regions: iteration is chosen multiple tissue regions and prevented extrapolation error; 3) the polynomial expression accumulation in unicoil estimation, accurately can estimate coil sensitivity functions by accumulation polynomial expression while not amplifying noise.Coil image does not bring next iteration into by carrying out correction to it due to low signal-to-noise ratio, and existing method is all re-use the estimation that lower order polynomial expressions goes as sensitivity in each iteration, but can not simulate the variation tendency of sensitivity so completely.A kind of method that the present invention utilizes level to accumulate, when ensureing to accumulate the polynomial expression of successive ignition when not amplifying coil image noise, more accurately can estimate coil sensitivity.

Description

A kind of MRI coil sensitivity estimation method based on the matching of many tissue regions

Technical field

The present invention relates to the technical field of pattern-recognition, Medical Image Processing, be specifically related to a kind of MRI coil sensitivity estimation method based on the matching of many tissue regions.

Background technology

Application background of the present invention is: Magnetic resonance imaging after it is invented due to its high resolving power, be widely used in medical diagnosis without the advantage such as ionization radiation injury and arbitrarily angled imaging.Current multi-coil imaging technique accelerates due to its imaging and improves the advantage of signal noise ratio (snr) of image and be widely used on imaging apparatus.But owing to needing to use multiple surface coils collection signal in multi-coil imaging technique, and these coils have uneven coil sensitivity, which results between the gray-scale value of coil image and actual value and there is certain deviation.Coil sensitivity causes because the ability of receiving coil Received signal strength is relevant to locus, when the tissue of excitation signal and the distance of coil larger time, the signal that coil receives is just more weak, and the signal of the tissue nearer with coil distance is stronger, the gray scale that the voxel which results in same tissue is corresponding is in the picture inconsistent.The uneven meeting of coil sensitivity causes the gray scale of composograph uneven (see document [1] RoemerPB, EdelsteinWA, HayesCE, SouzaSP, MuellerOM.TheNMRphasedarray.MagnResonMed1990; 16:1992-225.).In addition, coil sensitivity is also the important parameter rebuilding image in some parallel MRI imagings.Some parallel imaging techniques utilize the spatial information of coil array to carry out space encoding, add fast scan speed by reducing phase step number encoder.Only have and accurately calculate coil sensitivity and could restore complete coil image (see document [2] RoemerPB, EdelsteinWA, HayesCE, SouzaSP, MuellerOM.TheNMRphasedarray.MagnResonMed1990; 16:1992-225.).

Related art is:

Scheme one

Scenario Name: based on the coil sensitivity estimation of body coil (see document [2] RoemerPB, EdelsteinWA, HayesCE, SouzaSP, MuellerOM.TheNMRphasedarray.MagnResonMed1990; 16:1992-225.)

Algorithm idea: before formally scanning, uses special sensitivity calibration device first to obtain the full FOV image of the low resolution of body coil; Then the full FOV image of each low resolution in phased array coils is obtained.Using body coil as with reference to image, obtain the original sensitivities of each coil by each coil image divided by body coil, then carry out smooth treatment.

Algorithm shortcomings: carry out prescan to patient, needs to use extra hardware device, increases sweep time.And due to will twice sweep be carried out, if patient there occurs change in the position of formal scanning and prescan, make to there is deviation between coil sensitivity and sweep object, cause motion artifacts.

Scheme two

Scenario Name: based on the coil sensitivity estimation of low-pass filtering (see document [3] BrinkmannBH; ManducaA, RobbRA.OptimizedhomomorphicunsharpmaskingforMRgrayscalei nhomogeneitycorrection.IEEETransMedImaging1998; 17:161-171.)

Algorithm idea: deviation field that coil sensitivity is regarded as a low frequency, that slowly change, carry out low-pass filtering to coil image, the low-frequency information extracting image is the estimation of coil sensitivity

Algorithm shortcomings: the low-frequency information of coil sensitivity and image exists aliasing, low pass filtered wave frequency is bad to be determined.There is boundary effect in the high-contrast area of image, the details of boundary is lost, and needs smoothing process.

Scheme three

Scenario Name: based on Main Tissues extract coil sensitivity estimation (see document [4] VemuriP, KholmovskiE, ParkerD, ChapmanB.CoilsensitivityestimationforoptimalSNRreconstru ctionandintensityinhomogeneitycorrectioninphasedarrayMRI imaging.In:Proceedingsofthe19thInternationalConferenceof IPMI, GlenwoodSprings, CO, USA, 2005.p405-421.)

Algorithm idea: utilize the hypothesis that same in-house voxel gray values is consistent, first unicoil imagery exploitation SoS method is synthesized, the character utilizing Tissue distribution grey level histogram to obey Gaussian mixtures extracts a tissue regions roughly, according to the sensitivity function in the extrapolated whole FOV region of the grey scale change in this tissue regions from the image of synthesis.This Main Tissues region is upgraded by iteration.

Algorithm shortcomings: the distribution of each tissue is destroyed by the heterogeneity that coil sensitivity causes, and no longer obeys Gaussian mixtures, the region difference of searching is large, larger on the impact of extrapolation.And if organizing is not be uniformly distributed in the picture, then the error of extrapolating is very large.Cannot reach accurate to the estimation of coil sensitivity.Iteration only have updated tissue regions, at every turn still again with second order polynomial extrapolation, accurately cannot approach coil sensitivity functions.

Summary of the invention

The object of the invention is: 1), comparatively accurately can estimate coil sensitivity; 2), method of estimation does not rely on other equipment; 3), method of estimation arithmetic speed meets requirement of real time.

Technical solution of the present invention is: a kind of MRI coil sensitivity estimation method based on the matching of many tissue regions, utilize the gradient information determination tissue regions of image, not only ensure that the accuracy in region but also maintain larger gray scale difference, utilize polynomial expression to accumulate simultaneously and increase fitting precision, the idiographic flow of the method is as follows:

1, pre-service

Will carry out denoising to original coil image before estimating, pre-service further comprises the Outside contour extraction of image, rejects the background area at edge;

2, the determination of area-of-interest

Utilize region growth method, first area-of-interest is by determining that some Seed Points carry out initialization, and by judging whether the neighbours of these Seed Points similarly to seed to determine whether being added area-of-interest, concrete steps are as follows:

2.1, the choosing of seed region

Seed Points constitutes seed region, Seed Points belongs to same tissue regions with larger probability, define indicator M (r) to mark a voxel and whether belong to seed region, when a voxel r belongs to seed region, M (r)=1, otherwise M (r)=0, determines M (r) by following rule:

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

Wherein I ₀r () represents the gray-scale value of voxel r in initial pictures, the image that initial pictures is rebuild by last iteration obtains, iteration utilizes the simple average of each coil image to obtain for the first time, p be remove background area in composograph histogram after peak value, σ is the noise variance of composograph, and final seed region is expressed as the initialization area of Here it is region growing algorithm;

2.2, region increases

Once determine seed region, can be just initial area-of-interest with Seed Points, constantly similar point is added to come in expand area-of-interest, by comparing with its eight neighborhood point is capable the point in each area-of-interest, if the gradient difference of the two is less than certain threshold value, then think that this point is similar to the point in area-of-interest, added area-of-interest, if run into border or the critical part with its hetero-organization, then similitude can not be it can be used as to add because Grad is excessive, by continuous iteration, area-of-interest constantly increases, till it no longer changes, suppose determining R area-of-interest, iterate to n-th time, then the region of interest area update of n-th time is as follows:

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

Wherein, be for initialized seed region, δ is the greatest gradient allowed, and its value is determined by empirical value, can regulate according to different images, to be set in seed region the half of Grad a little.Point in the eight neighborhood of Neigh (q) representative point q, Grad (p) is the Grad of the some p utilizing sobel operator to calculate, and when the quantity of the point newly added no longer changes, stop area increases;

The determination of 2.3 multiple semi-cylindrical hills

Estimate to reduce evaluated error to sensitivity by finding multiple semi-cylindrical hills;

3, the estimation of coil sensitivity

Due to the division that area-of-interest is to histological types, the gray-scale value of the point therefore in each region should be consistent, for each coil image, the change of the gray-scale value in area-of-interest can be thought to be caused by coil sensitivity, but the true gray-scale value of these tissue regions is unknowable, suppose that we find out m area-of-interest, be labeled as S _i(i=1 ..., m), in i-th region, the gray-scale value of some x is S _i(x), S _iinterior gray average is set to μ _i, so the coil sensitivity at x point place can be estimated as roughly:

${\hat{g}}_i\left(x\right)=\frac{S_i\left(x\right)}{{\mathrm{\&mu;}}_i}$

Due to the error that region increases, S _icontain to interior possible errors the tissue points of different tissues, so just likely cause the high-contrast between different tissues to be reflected in acute variation on, this does not just meet the hypothesis that coil sensitivity is slowly change, can bring very large error, in order to suppress this impact, to S when matching yet _iin point carrying out a low-pass filtering treatment, obtain a slightly good estimation

${\stackrel{\&OverBar;}{g}}_i\left(x\right)=\mathrm{LPF}\left({\hat{g}}_i(\&CenterDot;)\right)$

Mark the set of all area-of-interests supposing in S that one co-exists in Q point, is r respectively ₁, r ₂..., r _q, so the coil sensitivity initial estimation of these points is (namely corresponding value), suppose that choosing k rank polynomial expression carries out matching, comprises every being respectively wherein l+m≤k, l>=0, m>=0(r _xwith r _ythe coordinate of level and vertical direction), so total K=(k+1) (k+2)/2 such item, is labeled as, F _i(r), i=1,2 ..., K, if F in polynomial expression _ir the coefficient of () is ω _i, then these coefficients meet following relation:

$\stackrel{\&OverBar;}{G}=\mathrm{FW}$

Wherein, $\stackrel{\&OverBar;}{G}=\left[\begin{array}{c}\stackrel{\&OverBar;}{g}\left(r_1\right)\\ \stackrel{\&OverBar;}{g}\left(r_2\right)\\ .\\ .\\ .\\ \stackrel{\&OverBar;}{g}\left(r_Q\right)\end{array}\right],$ $W=\left[\begin{array}{c}{\mathrm{\&omega;}}_1\\ {\mathrm{\&omega;}}_2\\ .\\ .\\ .\\ {\mathrm{\&omega;}}_K\end{array}\right].$

Utilize least square method, obtain:

W＝(F ^tF) ^-1F ^tS

Wherein t and-1 represents transpose of a matrix and inverse respectively, solves W, can derive the coil sensitivity at full images each some r:

g(r)＝[F ₁(r),F ₂(r),…,F _K(r)]W

Be more than the overall process being estimated coil sensitivity g (r) by fitting of a polynomial, by this process markup be:

$g\left(r\right)=\mathrm{poly}\left(\stackrel{\&OverBar;}{g}\left(r\right)\right)$

4, iteration

By iteration repeatedly lower order polynomial expressions approach accurate coil sensitivity, in jth time iteration, estimated the jth order polynomial β of c coil by following formula _c ^(j)(r), that is:

${{\mathrm{\&beta;}}_c}^{\left(j\right)}\left(r\right)=\mathrm{poly}({\stackrel{\&OverBar;}{g}}_c\left(r\right)/{g_c}^{(j-1)}\left(r\right))$

Wherein be the initial sensitivity estimation after the low-pass filtering of c coil, g _c ^(j-1)r () is the accumulation of front j-1 order polynomial:

${g_c}^{(j-1)}\left(r\right)={\mathrm{\&Pi;}}_{i=1}^{j-1}{\mathrm{\&beta;}}_c^i\left(r\right),{g_c}^{\left(0\right)}\left(r\right)=1.$

After obtaining the sensitivity estimation of all coils, synthesize the new reconstruction image of a width by following formula:

$I_{\mathrm{re}}\left(r\right)=\frac{R\left(r\right){\mathrm{\&Psi;}}^{-1}{\stackrel{\&OverBar;}{B}}^H\left(r\right)}{\stackrel{\&OverBar;}{B}\left(r\right){\mathrm{\&Psi;}}^{-1}{\stackrel{\&OverBar;}{B}}^H\left(r\right)}$

Wherein R (r) represents the gray scale vector of the coil image at pixel r place, i.e. R (r)=[R ₁(r), R ₂(r) ..., R _l(r)], for coil sensitivity vector $\stackrel{\&OverBar;}{B}\left(r\right)=[{\stackrel{\&OverBar;}{b}}_1\left(r\right),{\stackrel{\&OverBar;}{b}}_2\left(r\right),\&CenterDot;\&CenterDot;\&CenterDot;,{\stackrel{\&OverBar;}{b}}_L\left(r\right)],$ Ψ is coil coupling matrix, is generally set to unit matrix, and the image after reconstruction goes to upgrade area-of-interest as the initial pictures of next iteration, and stopping criterion for iteration substantially no longer changes when area-of-interest.

The advantage of technical solution of the present invention and good effect are:

(1) when, determining that multiple tissue regions prevents matching, extrapolation error is excessive

Technical solution of the present invention is chosen multiple tissue regions and is carried out fitting of a polynomial simultaneously, effectively prevent single organization's areal distribution inequality and causes sensitivity to estimate the undue problem relying on extrapolation.

(2), by region growth method determination tissue regions

The method utilizing histogram to delimit threshold value determination area-of-interest has very large defect, first coil image is caused deformation by sensitivity function interference, the gray-scale value of different tissues is likely twisted same gray scale rank, and the tissue regions determined with larger scope so likely comprises into different tissues.Although and the probability that less scope comprises into the pixel of same tissue is comparatively large, very few gray scale difference causes to extract sensitivity information completely.And region growth rule effectively avoids this problem, can ensure to choose enough accurate pixel, can ensure that again selected point has enough large gray scale difference to simulate sensitivity function.

(3), coil sensitivity is approached by polynomial expression accumulation

Coil image does not bring next iteration into by carrying out correction to it due to low signal-to-noise ratio, and existing method is all re-use the estimation that lower order polynomial expressions goes as sensitivity in each iteration, but can not simulate the variation tendency of sensitivity so completely.A kind of method that the present invention utilizes level to accumulate, when ensureing to accumulate the polynomial expression of successive ignition when not amplifying coil image noise, more accurately can estimate coil sensitivity.

Accompanying drawing explanation

Fig. 1 is a kind of MRI coil sensitivity estimation method flow schematic diagram based on the matching of many tissue regions.

Fig. 2 uses the image effect comparison diagram estimating coil sensitivity synthesis.Wherein, different human body tissue effect contrasts, and left column is for using SoS synthesis, and centre is classified as and uses single region not use, and the sensitivity of multimodal accumulation estimates that the image of synthesis, the right side are classified as the image of our innovatory algorithm synthesis.

Embodiment

The present invention is further illustrated below in conjunction with accompanying drawing 1 and embodiment.

The present invention finds multiple tissue regions (area-of-interest) by utilizing region growth method, because the grey scale change in each tissue regions is caused by same coil sensitivity, therefore this method voxel intensity of combining in each tissue regions solves same polynomial expression, carrys out matching coil sensitivity functions.In an iterative process, this method not only have updated tissue regions, and adds up to the polynomial expression of iterative fitting before, thus finally obtains comparatively accurate coil sensitivity estimation (detail flowchart is shown in 5 joints).Idiographic flow is as follows:

1, pre-service

Because the present invention has used region growing algorithm, larger to the gradient dependence of image pixel.But gradient information receives the interference of the high-frequency signals such as noise, therefore denoising to be carried out to original coil image before estimating.Pre-service further comprises the Outside contour extraction of image, reject the background area at edge, because only have noise in background area, and sensitivity function can be extrapolated to background area by the extrapolation that sensitivity is used in estimating, thus be exaggerated noise when rebuilding image, so need to extract an accurate picture appearance profile, the simple histogram thresholding of the present invention extracts image outline.

2, the determination of area-of-interest

A tissue in each area-of-interest correspondence image.Grey scale pixel value in identical tissue should be consistent, but the impact of coil sensitivity is subject to when acknowledge(ment) signal, in region there is the difference of slowly change in gray-scale value, therefore by finding area-of-interest to extract the coil sensitivity of local, then can be extrapolated to by matching the coil sensitivity that full images can obtain estimation.Utilize region growth method, first area-of-interest is by determining that some Seed Points carry out initialization, by judging whether the neighbours of these Seed Points similarly to seed to determine whether being added area-of-interest.The voxel of different tissues can be effectively avoided to be marked as the mistake of the same area like this.Concrete steps are as follows:

2.1, the choosing of seed region

Seed Points constitutes seed region.Seed Points belongs to same tissue regions with larger probability.We define indicator M (r) and mark a voxel and whether belong to seed region.When a voxel r belongs to seed region, M (r)=1, otherwise, M (r)=0.We determine M (r) by following rule:

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

Wherein I ₀r () represents the gray-scale value of voxel r in initial pictures, the image that initial pictures is rebuild by last iteration obtains, first time iteration utilize the simple average of each coil image to obtain, p be remove background area in composograph histogram after peak value, σ is the noise variance of composograph.Because the gray-value variation in seed region is no more than noise variance, we believe the very large probability of point in this region all belong to same tissue.But be subject to the impact that gray scale is uneven, still the voxel of its hetero-organization may be included, in order to reduce this impact as far as possible, we carry out a filtering process to seed region again, weed out too sparse point, the point only having those enough dense just can be retained in seed region.Final seed region is expressed as the initialization area of Here it is region growing algorithm.

2.2, region increases

Once determine seed region, can be just initial area-of-interest with Seed Points, constantly similar point be added to come in expand area-of-interest.We are by comparing with its eight neighborhood point is capable the point in each area-of-interest, if the gradient difference of the two is less than certain threshold value, then think that this point is similar to the point in area-of-interest, added area-of-interest, if run into border or the critical part with its hetero-organization, then similitude can not be it can be used as to add because Grad is excessive.By continuous iteration, area-of-interest constantly increases, till it no longer changes.Suppose that we are determining R area-of-interest, iterate to n-th time, then the region of interest area update of n-th time is as follows:

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

Wherein, be for initialized seed region, δ is the greatest gradient allowed, and its value is determined by empirical value, can regulate according to different images, we to be set in seed region the half of Grad a little.Point in the eight neighborhood of Neigh (q) representative point q, Grad (p) is the Grad of the some p utilizing sobel operator to calculate.When the quantity of the point newly added no longer changes, stop area increases.

The determination of 2.3 multiple semi-cylindrical hills

We estimate to reduce evaluated error to sensitivity by finding multiple semi-cylindrical hills.Once determine an area-of-interest, we just it are removed from image and repeat region increases searching second area-of-interest.Tissue regions all in image can be found in theory, but due to region growth method be subject to initial point selection and empirical value to affect application condition large, we arrange and only find two area-of-interests in the application.

3, the estimation of coil sensitivity

Due to the division that area-of-interest is to histological types, the gray-scale value of the point therefore in each region should be consistent.For each coil image, the change of the gray-scale value in area-of-interest can be thought to be caused by coil sensitivity, but the true gray-scale value of these tissue regions is unknowable.We are used as substituting of true gray-scale value by the average of these points in actual applications.Suppose that we find out m area-of-interest, be labeled as S _i(i=1 ..., m), in i-th region, the gray-scale value of some x is S _i(x), S _iinterior gray average is set to μ _i.So the coil sensitivity at x point place can be estimated as roughly:

${\hat{g}}_i\left(x\right)=\frac{S_i\left(x\right)}{{\mathrm{\&mu;}}_i}$

Due to the error that region increases, S _icontain to interior possible errors the tissue points of different tissues, so just likely cause the high-contrast between different tissues to be reflected in acute variation on, this does not just meet the hypothesis that coil sensitivity is slowly change.Also can bring very large error when matching, in order to suppress this impact, we are to S _iin point carrying out a low-pass filtering treatment, obtain a slightly good estimation

${\stackrel{\&OverBar;}{g}}_i\left(x\right)=\mathrm{LPF}\left({\hat{g}}_i(\&CenterDot;)\right)$

Mark the set of all area-of-interests supposing in S that one co-exists in Q point, is r respectively ₁, r ₂..., r _q, so the coil sensitivity initial estimation of these points is (namely corresponding value).Suppose that we choose k rank polynomial expression and carry out matching, comprise every being respectively wherein l+m≤k, l>=0, m>=0(r _xwith r _ythe coordinate of level and vertical direction), so total K=(k+1) (k+2)/2 such item, we are labeled as, F _i(r), i=1,2 ..., K.If F in polynomial expression _ir the coefficient of () is ω _i, then these coefficients meet following relation:

$\stackrel{\&OverBar;}{G}=\mathrm{FW}$

Wherein, $\stackrel{\&OverBar;}{G}=\left[\begin{array}{c}\stackrel{\&OverBar;}{g}\left(r_1\right)\\ \stackrel{\&OverBar;}{g}\left(r_2\right)\\ .\\ .\\ .\\ \stackrel{\&OverBar;}{g}\left(r_Q\right)\end{array}\right],$ $W=\left[\begin{array}{c}{\mathrm{\&omega;}}_1\\ {\mathrm{\&omega;}}_2\\ .\\ .\\ .\\ {\mathrm{\&omega;}}_K\end{array}\right].$

Utilize least square method, obtain:

W＝(F ^tF) ^-1F ^tS

Wherein t and-1 represents transpose of a matrix and inverse respectively.Solve W, the coil sensitivity at full images each some r can be derived:

g(r)＝[F ₁(r),F ₂(r),…,F _K(r)]W

Be more than the overall process being estimated coil sensitivity g (r) by fitting of a polynomial, this process markup is by we:

$g\left(r\right)=\mathrm{poly}\left(\stackrel{\&OverBar;}{g}\left(r\right)\right)$

4, iteration

Although coil sensitivity is a deviation field slowly changed, it can not go to represent with the polynomial expression of a Ge Di class completely.So the estimated result that above-mentioned estimation procedure can not be satisfied with completely.We can by iteration repeatedly lower order polynomial expressions approach accurate coil sensitivity.In jth time iteration, we estimate the jth order polynomial β of c coil by following formula _c ^(j)(r), that is:

${{\mathrm{\&beta;}}_c}^{\left(j\right)}\left(r\right)=\mathrm{poly}({\stackrel{\&OverBar;}{g}}_c\left(r\right)/{g_c}^{(j-1)}\left(r\right))$

Wherein it is the initial sensitivity estimation after the low-pass filtering of c coil.G _c ^(j-1)r () is the accumulation of front j-1 order polynomial:

${g_c}^{(j-1)}\left(r\right)={\mathrm{\&Pi;}}_{i=1}^{j-1}{\mathrm{\&beta;}}_c^i\left(r\right),{g_c}^{\left(0\right)}\left(r\right)=1.$

After obtaining the sensitivity estimation of all coils, we synthesize the new reconstruction image of a width by following formula:

$I_{\mathrm{re}}\left(r\right)=\frac{R\left(r\right){\mathrm{\&Psi;}}^{-1}{\stackrel{\&OverBar;}{B}}^H\left(r\right)}{\stackrel{\&OverBar;}{B}\left(r\right){\mathrm{\&Psi;}}^{-1}{\stackrel{\&OverBar;}{B}}^H\left(r\right)}$

Wherein R (r) represents the gray scale vector of the coil image at pixel r place, i.e. R (r)=[R ₁(r), R ₂(r) ..., R _l(r)], for coil sensitivity vector $\stackrel{\&OverBar;}{B}\left(r\right)=[{\stackrel{\&OverBar;}{b}}_1\left(r\right),{\stackrel{\&OverBar;}{b}}_2\left(r\right),\&CenterDot;\&CenterDot;\&CenterDot;,{\stackrel{\&OverBar;}{b}}_L\left(r\right)],$ Ψ is coil coupling matrix, is generally set to unit matrix.Image after reconstruction goes to upgrade area-of-interest as the initial pictures of next iteration.Stopping criterion for iteration substantially no longer changes (variation range is less than 0.1%) when area-of-interest.

Notice that we are not direct in an iterative process and " correction " (that is, I is carried out to raw coil image ₁(r)/g ₁ ^(j)(r) ..., I _l(r)/g _l ^(j)(r)) reach polynomial expression accumulation object, but when each iteration, the polynomial expression of iteration before is first accumulated to together, use the original gradation of coil divided by accumulation polynomial expression again, matching is carried out to the business of remainder, the main cause done like this has a large amount of low signal-to-noise ratio regions in coil image, if each iteration all corrects coil image, then low signal-to-noise ratio area grayscale is by too high, also be greatly exaggerated noise simultaneously, last composograph has flooded with regard to the full noise be exaggerated of set, therefore the bearing calibration adopting above-mentioned " level accumulation " is only had, the method of the first accumulation polynomial expression matching again when each iteration, just can avoid the problem of amplifying noise.

5, testing process of the present invention

The flow process that the present invention carries out coil sensitivity estimation is as follows:

STEP1): pre-service, the denoising of coil image and the extraction of image outline is comprised;

STEP2): determine multiple tissue regions (area-of-interest), first each tissue regions by determining some Seed Points, and then region increases acquisition.Then reject this region in the picture to enter iteration and determine next tissue regions.

STEP3): carry out low-pass filtering to each area-of-interest, the impact that the error reducing region growth is brought.

STEP4): the local coil sensitivity in associating multiple semi-cylindrical hills is by the coil sensitivity of the extrapolated global image of fitting of a polynomial.

STEP5): the coil sensitivity utilizing current iteration to obtain rebuilds image, for the renewal of area-of-interest in next iteration.

Claims (1)

1. the MRI coil sensitivity estimation method based on the matching of many tissue regions, it is characterized in that, utilize the gradient information determination tissue regions of image, not only ensure that the accuracy in region but also maintain larger gray scale difference, utilize polynomial expression to accumulate simultaneously and increase fitting precision, the idiographic flow of the method is as follows:

1, pre-service

Will carry out denoising to original coil image before estimating, pre-service further comprises the Outside contour extraction of image, rejects the background area at edge;

2, the determination of area-of-interest

Utilize region growth method, first area-of-interest is by determining that some Seed Points carry out initialization, and by judging whether the neighbours of these Seed Points similarly to seed to determine whether being added area-of-interest, concrete steps are as follows:

2.1, the choosing of seed region

Seed Points constitutes seed region, Seed Points belongs to same tissue regions with larger probability, define indicator M (r) to mark a voxel and whether belong to seed region, when a voxel r belongs to seed region, M (r)=1, otherwise M (r)=0, determines M (r) by following rule:

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

Wherein I ₀r () represents the gray-scale value of voxel r in initial pictures, the image that initial pictures is rebuild by last iteration obtains, iteration utilizes the simple average of each coil image to obtain for the first time, p be remove background area in composograph histogram after peak value, σ is the noise variance of composograph, and final seed region is expressed as the initialization area of Here it is region growing algorithm;

2.2, region increases

Once determine seed region, can be just initial area-of-interest with Seed Points, constantly similar point is added to come in expand area-of-interest, by comparing with its eight neighborhood point is capable the point in each area-of-interest, if the gradient difference of the two is less than certain threshold value, then think that this point is similar to the point in area-of-interest, added area-of-interest, if run into border or the critical part with its hetero-organization, then similitude can not be it can be used as to add because Grad is excessive, by continuous iteration, area-of-interest constantly increases, till it no longer changes, suppose determining R area-of-interest, iterate to n-th time, then the region of interest area update of n-th time is as follows:

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

Wherein, for initialized seed region, δ is the greatest gradient allowed, its value is determined by empirical value, can regulate according to different images, to be set in seed region the half of Grad a little, the point in the eight neighborhood of Neigh (q) representative point q, Grad (p) is the Grad of the some p utilizing sobel operator to calculate, when the quantity of the point newly added no longer changes, stop area increases;

The determination of 2.3 multiple semi-cylindrical hills

Estimate to reduce evaluated error to sensitivity by finding multiple semi-cylindrical hills;

3, the estimation of coil sensitivity

Due to the division that area-of-interest is to histological types, the gray-scale value of the point therefore in each region should be consistent, for each coil image, the change of the gray-scale value in area-of-interest can be thought to be caused by coil sensitivity, but the true gray-scale value of these tissue regions is unknowable, suppose that we find out m area-of-interest, be labeled as S _i(i=1 ..., m), in i-th region, the gray-scale value of some x is S _i(x), S _iinterior gray average is set to μ _i, so the coil sensitivity at x point place can be estimated as roughly:

${\hat{g}}_i\left(x\right)=\frac{S_i\left(x\right)}{{\mathrm{\&mu;}}_i}$

Due to the error that region increases, S _icontain to interior possible errors the tissue points of different tissues, so just likely cause the high-contrast between different tissues to be reflected in acute variation on, this does not just meet the hypothesis that coil sensitivity is slowly change, can bring very large error, in order to suppress this impact, to S when matching yet _iin point carrying out a low-pass filtering treatment, obtain a slightly good estimation

${\stackrel{\&OverBar;}{g}}_i\left(x\right)=LPF\left({\hat{g}}_i\right(x\left)\right)$

Mark the set of all area-of-interests supposing in S that one co-exists in Q point, is r respectively ₁, r ₂..., r _q, so the coil sensitivity initial estimation of these points is namely corresponding value, supposes that choosing k rank polynomial expression carries out matching, comprises every being respectively wherein l+m≤k, l>=0, m>=0, r _xwith r _ybe the coordinate of level and vertical direction, so total K=(k+1) (k+2)/2 such item, is labeled as, F _i(r), i=1,2 ..., K, if F in polynomial expression _ir the coefficient of () is ω _i, then these coefficients meet following relation:

$\stackrel{\&OverBar;}{G}=FW$

Wherein, $\stackrel{\&OverBar;}{G}=\left[\begin{array}{c}\stackrel{\&OverBar;}{g}\left(r_1\right)\\ \stackrel{\&OverBar;}{g}\left(r_2\right)\\ \&CenterDot;\\ \&CenterDot;\\ \&CenterDot;\\ \stackrel{\&OverBar;}{g}\left(r_Q\right)\end{array}\right],$ $W=\left[\begin{array}{c}{\mathrm{\&omega;}}_1\\ {\mathrm{\&omega;}}_2\\ \&CenterDot;\\ \&CenterDot;\\ \&CenterDot;\\ {\mathrm{\&omega;}}_K\end{array}\right].$

Utilize least square method, obtain:

W＝(F ^tF) ^-1F ^tS

Wherein t and-1 represents transpose of a matrix and inverse respectively, solves W, can derive the coil sensitivity at full images each some r:

g(r)＝[F ₁(r),F ₂(r),…,F _K(r)]W

Be more than the overall process being estimated coil sensitivity g (r) by fitting of a polynomial, by this process markup be:

$g\left(r\right)=poly\left(\stackrel{\&OverBar;}{g}\right(r\left)\right)$

4, iteration

By iteration repeatedly lower order polynomial expressions approach accurate coil sensitivity, in jth time iteration, estimated the jth order polynomial β of c coil by following formula _c ^(j)(r), that is:

${{\mathrm{\&beta;}}_c}^{\left(j\right)}\left(r\right)=poly\left({\stackrel{\&OverBar;}{g}}_c\right(r)/{g_c}^{(j-1)}(r\left)\right)$

Wherein be the initial sensitivity estimation after the low-pass filtering of c coil, g _c ^(j-1)r () is the accumulation of front j-1 order polynomial:

${g_c}^{(j-1)}\left(r\right)={\mathrm{\&Pi;}}_{i=1}^{j-1}{\mathrm{\&beta;}}_c^i\left(r\right),$ g _c ⁽⁰⁾(r)＝1，

After obtaining the sensitivity estimation of all coils, synthesize the new reconstruction image of a width by following formula:

$I_{re}\left(r\right)=\frac{R\left(r\right){\mathrm{\&Psi;}}^{-1}{\stackrel{\&OverBar;}{B}}^H\left(r\right)}{\stackrel{\&OverBar;}{B}\left(r\right){\mathrm{\&Psi;}}^{-1}{\stackrel{\&OverBar;}{B}}^H\left(r\right)}$

Wherein R (r) represents the gray scale vector of the coil image at pixel r place, i.e. R (r)=[R ₁(r), R ₂(r) ..., R _l(r)], for coil sensitivity vector Ψ is coil coupling matrix, is generally set to unit matrix, and the image after reconstruction goes to upgrade area-of-interest as the initial pictures of next iteration, and stopping criterion for iteration substantially no longer changes when area-of-interest.