The present invention relates to methods of processing MRI (magnetic resonance imaging) scans, particularly DWI (diffusion-weighted images) MRI scans.
Generally, there are two types of errors [1] in any observation: systematic and random. Systematic errors tend to shift all measurements in a particular direction. Some of the main reasons of such errors are incorrect calibration of an instrument, improper use of the instrument, etc. Large systematic errors can be often be eliminated (e.g. by applying zero correction of the instrument or repeating the experiment), but small systematic errors will always be present since no instrument can ever be calibrated perfectly. This is the reason why several independent confirmations of experimental results should be performed, preferably using different techniques.
If an experiment is performed several times with all experimental conditions constant, the outcome is still different. These fluctuations in the outcome are called random errors (or statistical errors). The value of the outcome is taken as the mean of observations and the standard deviation is taken as the error on the mean. The standard deviation can sometimes be obtained by repeating the experiment, but in some practical situations it is impossible to repeat experiments. In these situations, the knowledge of the distribution of outcome is applied to predict statistical errors. The outcome usually follows certain known distributions depending on the nature of the experiment e.g. the Poisson distribution is a common outcome in experiments which include a count. For the Poisson distribution, standard deviation (σ) is related to mean (μ) as σ=√{square root over (μ)} [1]. Because of the relationship between μ and σ, one is able to predict an error from the outcome of the experiment (where the outcome is a result of counts per unit time); for example, Ref. [2] used a Poisson distribution-based noise removal technique for nuclear medical imaging since such imaging involves a number of decays per unit time.
The process of MRI acquisition is very complicated (http://www.easymeasure.co.uk/principlesmri.aspx, http://www.sunnybrook.ca/research/groups/cardiac_mri/MR_background, accessed Oct. 23, 2007). MRI signal intensity has a complicated dependence on many parameters including the count of magnetically excited protons in the voxel during the image acquisition time. Since intensity is in part also related to the count of magnetically excited protons, the Poisson distribution is used to predict the distribution of intensity of each pixel. In the light of this assumption, an error on the pixel intensity can be predicted. Thus, the reported pixel value can be assumed to have error equal to √{square root over (pμ)} where pμ is the mean pixel intensity of some hypothetical observations.
There are various kinds of acquisition artifacts associated with MRI scans (some are describe at http://www.mritutor.org/mritutor/artifact.htm, accessed Oct. 23, 2007) such as motion artifacts, aliasing artifacts, susceptibility artifacts etc. Some known methods to remove artifacts and reduce noise are as follows. Ref [3] presents a Wavelet-based Rician noise removal for MRI. Ref [4] describes an approach to noise filtering in multi-dimensional data using a partial volume data density model. Ref [5] suggests correcting bulk in-plane motion artifacts in MRI using a point spread function. A more elaborate list of these methods is included in (http://iris.usc.eduNision-Notes/bibliography/medical891.html, accessed Oct. 23, 2007).
Computer aided detection (CAD) plays a significant role in aiding accurate medical image interpretations in different areas [e.g. 6-10]. The present inventors have developed a suite of CAD systems for acute ischemic and hemorrhagic strokes [11-13]. One of key algorithms is segmentation of infarcts. Its accuracy depends on correct discrimination of infarct from artifacts. Accurate and rapid quantification of infarcts from DWI scans is critical in acute ischemic strokes. Acquisition artifacts lead to hyperintense regions in DWI MR scans resulting in false positives. Discriminating infarcts and artifacts helps to reduce infarct segmentation errors.
The present invention relates to post-processing segmented MRI images to increase the accuracy of infarct delineation.
In general terms, the algorithm proposes that an MRI image of a brain, such as a 3D DWI image comprising a plurality of 2D DWI scans, which has been segmented based on the intensity of the pixels in the scan to identify hyperintense regions of a brain which are candidates to correspond to infarct tissue, is processed to eliminate identified regions for which this identification was incorrect. This is done by one or more of: eliminating identified regions which are determined to be similar to the region of the scan which corresponds to the identified region reflected in the mid-sagittal plane (MSP) of the brain; and eliminating regions which are determined not to have corresponding identified regions in one or more of the other scans.
The proposed algorithm may make it possible to discriminate between infarcts and artifacts in DWI scans, and thereby reduce errors in morphological measurements.
The criterion for evaluating the similarity of symmetrically-related hyperintense regions may employ a numerical parameter which is related to the Poisson error in the intensity of each pixel. This is because the expected error in the intensity of each pixel relative to a perfect measurement of the intensity (the “intensity space” of the pixel) is typically given by a normal distribution independently of the nature of the experiment.
Two applications of the present technique are: determination that there is insufficient evidence that a given 2D scan exhibits an infarct (e.g. if, following one or both of the elimination processes proposed above, and in particular the step of eliminating symmetric regions, the amount of the remaining infarct regions does not meet a threshold); and, in an 2D scan which does exhibit an infarct, removing regions which are erroneously identified as infarcts.
The algorithm has the potential to remove artifacts from any infarct processing system. In particular, this approach has application to investigations of thrombolysis using DWT scans, and to quantify morphological properties of a newly discovered infarct. Once the algorithm above has been used to produce a post-processed image, that image may be used to quantify (i) the diffusion perfusion mismatch and (ii) size of infarct to that of MCA ratio.
Note that apart from DWI, other data acquisition techniques such as FLAIR (fluid attenuation inversion recovery), T2, ADC (apparent diffusion coefficient) are used for infarct staging. It is presently considered that the present techniques are of most interest in quantifying newly-identified infarcts, and for these DWI is most valuable, but the technique is applicable wherever the “signal of interest” or “detection of a disease” is brighter than the rest of the image. So irrespective of the type of images this technique is applicable is applicable.
The present algorithm may be implemented by a computer system. If so, it is typically performed automatically (which is here used to mean that, although human interaction may initiate the algorithm, human interaction is not required while the algorithm is carried out). The algorithm might alternatively be performed semi-automatically (in which case there is human interaction with the computer during the processing).
A specific expression of the invention is a method of processing an MRI image of a brain, the MRI image comprising a plurality of 2D MRI scans corresponding to respective planes of the brain, the method including identifying in each scan one or more hyperintense regions which are candidates to correspond to infarct tissue in the brain;
Embodiments of the invention will now be described, for the sake of example only, with reference to the following drawings, in which:
Two processes which are embodiments of the invention will now be described. After this, we will discuss two applications which employ one or both of the processes as part of more complex algorithms, which are also embodiments of the invention. Finally we discuss experimental results of the these two applications.
The input to the first process is a 2D DWI scan (or a plurality of such slices, such as a plurality of axial scans of slices at respective heights in a patient's brain).
The flowchart of the first process (symmetric artifact removal) is presented in
The input 2D DWI scan is labeled in
In a second step 3, the hyperintense regions in the infarct hemisphere are labeled. This can be done by obtaining an intensity histogram of the infarct hemisphere. As is known from the prior art, a typical intensity histogram of an MRI image containing infarct material is as shown in
In step 4, for each hyperintense region in the infarct hemisphere, a corresponding mirror region (at the same distance from the MSP and of the same shape) in the non-infarct hemisphere is examined. The size of region is calculated.
The set of steps indicated as 5 in
First, in step 6, it is determined if the size of segmented region of the infract hemisphere is less than 5% of the total image size (excluding the background).
If the result of the determination of step 6 is “no”, then the situation is as in
Let Dj=pj−pj′ be the difference of intensities of pixels j and j′. From the Law of propagation of errors, [18], the total error on the difference of intensity is obtained (step 8) as:
where
are partial derivatives and δpj(=√{square root over (pj)}), δpj′, (=√{square root over (pj′)}) are errors on intensity of pixels j and j′.
We use this error to estimate the 95% confidence interval around the difference of intensities equal to 0. The pixels (j and j′) are considered to have similar intensities (i.e. there is no evidence that the intensity of one is due to an infarct), if it is determined (step 10) that the difference of their intensities lies in the 95% confidence region around zero i.e. Dj≦1.96Tj (http://mathworld.wolfram.com/ConfidenceInterval.html, accessed Oct. 23, 2007). More generally, the similarity criterion for regarding two pixels as having similar intensities can be written as Dj≦λ1Tj, where λ1 is a similarity parameter. Below we investigate the effects of varying λ1, which is equivalent to exploring other confidence intervals.
Alternatively, if the result of the determination of step 6 is “yes”, then the situation is as in
The error Ek in Rk is derived from the Law of propagation of errors [18] as:
where
is the partial derivative of Rk with respect to pj and δpj is the error on the intensity of jth pixel.
The difference of mean intensities of two regions is calculated (step 13) as:
D
k
=R
k
−R
k′.
The total error in the difference of mean intensities of regions is calculated as:
Here, δRk and δRk, are errors on Rk and Rk′ which are defined as Ek and Ek′. Any two regions k and k′ are considered to have similar intensities if it is determined (step 14) that Dk≦1.96Tk. Similar regions are symmetric regions with similar intensities. More generally, the similarity criterion can be varied, such that it is expressed as Dk≦λ1Tk where λ1 is again the similarity parameter, to explore other confidence intervals.
The symmetric regions and symmetric pixels with similar intensities are considered as artifacts. Specifically, if the determination in steps 9 and 14 is negative, the pixel in the infarct hemisphere is excluded from the set of identified infarct pixels (steps 10 and 15 respectively). Otherwise, it is confirmed that the pixel is indeed an infarct pixel.
The flowchart diagram of different steps of determining 3-D spatial coherence is given in
In step 22 we perform the following set of image processing sub-steps.
First, image dilation [19-20] is performed using a structuring element obtained by taking into account the spatial error around each pixel [http://www.cis.rit.edu/htbooks/mri/, http://www.sunnybrook.ca/research/groups/cardiac_mri/MR_background]. The surrounding region around each pixel can be regarded as the spatial error region.
The structuring element is illustrated in
Second, we determine 3-D connected regions in the volume [21]. Note that in each region, the dilated regions have a slightly different shape. The criterion for deciding that regions in consecutive scans are connected is that at least one pixel must be in common between the hyperintense regions of the consecutive scans.
Third, we calculate the number of slices in which a 3-D connected region occurs continuously, which is called slice frequency ν. For example in
Fourth, we determine the maximum ν which is denoted ν.
In step 23 we determine whether νmax/total infarct slices is a greater than a parameter indicating a significant fraction of the total number of slices. For example, for cases in which the number of interfarct slices is greater than one, we may take the significant fraction as 0.9. If this determination is negative, the process stops (step 24).
Otherwise, in step 25 find any regions with ν equal to 1. Note that a region may have ν equal to 1 even if a similar region appears at the same location (in the 2D space of the scans) after a gap of one of more slices (e.g. regions 2, 4 and 6 in
The first application of the processes above (especially the first process) is for identification of slices for which in fact there is insufficient evidence that infarct material is present. A flowchart to show the application is displayed in
The input to the application is a set of slices which have been identified as likely to contain infarct material, for example by an existing automatic slice identification algorithm [14] which also obtains the hemisphere in which the infarct is likely to be. This existing algorithm can be regarded as a first step 31 of the application, and corresponds to part of step 2 of
In step 32 the hyperintense regions in infarcted hemisphere are obtained by excluding the pixels below a threshold value, say ψ. The value of ψ is obtained as follows. As mentioned above, the second peak in the intensity distribution of a DWI scan (e.g.
We now set ψ=RH+λ2EH where λ2 is a second similarity parameter, and exclude all pixels with a lower intensity (step 34). Steps 33 and 34 correspond to step 3 of
We now perform symmetric region identification (step 35, corresponding to step 4 of
In step 38, we determine the number of infract pixels remaining in the slice after the exclusion, and whether this number of residual pixels is above or below a tolerance parameter. If the number is below the tolerance parameter, the slice is a false positive slice. If the number is above the tolerance parameter, the slice is confirmed as being an infarct slice.
In our experiments, we have taken the tolerance parameter as 0.01% of the total number of pixels in the image after excluding the background.
The second application is illustrated in
A first step 41 of the algorithm of
The next step 42 of the algorithm is to obtain the hemisphere which contains the infarct (this can be done, for example, by the method disclosed in [14]), and exclude all the hyperintense segmented regions in the other (“non-infarct”) hemisphere. That is, any regions of the non-infarct hemisphere which had previously been considered be candidate infarct regions, are relabeled such that they no longer are. This corresponds broadly to steps 2-3 of
The algorithm next (step 43) identifies symmetric artifacts and (step 44) excludes them. This corresponds to steps 4 and 5 of
The algorithm next (step 45) identifies further artifacts based on 3-D spatial coherence (i.e. the process of
The steps of artifact removal are displayed in
We now present experimental results using the processes and applications described above. Fifty one DWI scans were used in this study. This is data the we had used before.
(i) To test automatic slice identification (i.e. the application of
(ii) To test automatic infarct segmentation (i.e. the application of
(iii) 15 additional data set were used to demonstrate application of the proposed algorithm to improve results of a third algorithm [16]. The DWI scans had in-plane resolutions of 1.17 mm×1.17 mm to 2.42 mm×2.42 mm, slice thickness of 6.5-7 mm, and number of slices from 15-20.
Ground truth for all the data sets was marked by an expert.
The automatic Slice and Hemisphere identification algorithm [14] was aimed at automatically identifying the infarct slices and infarct hemisphere. The results of automatic detection of slice identification [14] were: Sensitivity=0.981, specificity=0.514, DSI [22]=0.665. After processing the data with the current technique (that is, the process of
Out of 36 cases, 26 cases showed an improvement in results due to false positive slice removal. The results of remaining 10 cases were unaffected by the processing. If we consider only those cases in which the results changed, the change in results is as follows: initial results for 26 data: (sensitivity, specificity, DSI)=(0.982, 0.474, 0.586). After processing the results for 26 data: (sensitivity, specificity, DSI)=(0.958, 0.664, 0.677). Thus, an increase in specificity and DSI is observed by 19% and 9.1%, respectively, with a decrease in sensitivity by 2.4%. The false negative slices which were removed had an insignificant fraction of area as compared to the maximum area of infarct in a slice. By using the proposed algorithm we are able to remove 31% of the false positive results.
a) shows the effect of changing λ1 and λ2(which, as described above, are employed in the criteria Dk≦λ1Tk and RH+λ2EH, respectively) on the sensitivity of infarct slice identification, and
The overall increase in (specificity, DSI) is (15.2%, 6.9%) with decrease in the sensitivity by only 1.5%.
The results of infarct segmentation algorithm in [15] were as follows: sensitivity=0.81, specificity=0.99 and DSI=0.60. After processing the data with proposed algorithm (i.e. the application of
Out of 13 volumes, all the cases showed an improvement in results due to artifact removal. Only 6 cases had DSI<0.5 out of the total of 13 cases. The average DSI for these 6 cases before processing with current algorithm was 18.3%. After processing with the proposed algorithm, the average increase of DSI is 11.2%. The seven cases which had an average DSI of 74.7% after post processing increased by 4.7%. Thus effect of current processing is more significant in cases where there are a large number of artifacts. The fraction of false positive pixels removed by the current algorithm is 71%.
Since the number of true negative pixels is of the order of total slice pixels and is much larger than the false positive pixels, specificity is always very large. It is hardly affected by any change in the number of false positive pixels. That is why no change in specificity is observed in
In [16] we segmented 15 data (5 each of low, medium and high artifact density). In another test of the application of
One of the goals of stroke CAD is to accurately and automatically identify, segment and measure the stroke region. This is important (a) in context of thrombolysis which requires quantifying the diffusion perfusion mismatch and size of infarct to that of MCA ratio (b) to provide input parameters for studies involving prognostic information like quantifying the impact of infarct location on stroke severity [23], quantification of patterns of DWI lesions [24] etc. While the state-of-the-art algorithms are being developed for achieving the final goal of stroke CAD, the presently proposed algorithms have stand alone applications in related areas of research. The embodiments make it possible to discriminate infarcts and artifacts based on the following two observational properties in DWI scans. They are motivated by two observations:
(i) A first observation is that a normal DWI scan in an axial plane shows both the hemispheres to have approximately similar features in terms of intensity, shape etc [e.g. 25]. Thus, if a DWI scan shows symmetric hyperintense regions in both hemispheres, they are most probably artifacts. An infarct caused by vessel occlusion most likely occurs in a single hemisphere so it will be much more hyperintense than the symmetric region in the opposite hemisphere. The embodiments quantify significant difference of intensity using the Poisson distribution in the intensity space of each pixel.
(ii) The second observation is that the infarct regions in different slices show spatial coherence. The regions that occur at distant locations from the spatially coherent regions are most probably artifacts. The embodiments detect spatial coherence by determining the overlap of dilated regions in different slices.
In many cases even the artifacts show 3-D spatial coherence. For that reason, the embodiment of
This is illustrated in
The embodiments also make use of the observation that the error in the intensity of each pixel is given by a normal distribution which is independent of the nature of the experiment and is generally associated with the randomness of the outcome of the experiment. However, since the mean and standard deviation of the normal distribution are independent, this distribution cannot be used to predict errors in the cases where it is not possible to repeat the outcome of experiment [1]. In fact, the present inventors initially considered using the FWHM of the background peak (as shown in
Identifying pixels which are symmetric about the MSP and have similar intensities to remove symmetric artifacts is very sensitive to errors such as: inherent asymmetry in hemispheres, the inter-hemispheric fissure being curved to a greater extent, etc. For that reason, the present embodiments by preference consider symmetric regions instead of individual pixels for the purpose of comparing the intensities. Even while considering the symmetric regions, due to inherent asymmetry of hemispheres about the MSP, regions lying close to the cortical surface boundary or too close to the ventricles or cerebrospinal fluid (CSF) may contain part of background or parenchyma. This is shown in
Based on the quantification of small, middle and large artifacts in [14], large infarct regions are expected to be those which are above about 5% of the image (excluding the background). Such large regions might have greater errors due to inhomogeneity (e.g. caused by intra slice variation during data observation). For this reason if a region is considered as a large region, the embodiment of
The 3-D spatial coherence is tested only in the cases where the ratio of νmax to the total number of infarct slices is at least equal to a ratio which is considered significant. The choice of high significant ratio as 0.9 is considered to establish the confidence that the infarct occurs at the similar location in the slices (cases with number of infarct slices=1 are excluded for this analysis). This is done to avoid cases in which the infarct occurs at multiple locations and there may be a significant chance that a spatially isolated region is an infarct. In the case that an infarct occurs mainly in one region, the chances that spatially isolated regions represent artifacts are high.
One of the assumptions for this algorithm is that the artifacts are symmetric about the MSP. Infarcts that are caused by vessel occlusion are most likely to be located in one hemisphere. In some rare cases (e.g., caused by a sudden drop in blood pressure in the presence of bilateral stenosis), infarcts may be present in both the hemispheres almost symmetrically. In such cases, only the spatial coherence test of the regions is performed. Symmetry based artifact removal should not be performed then. During data acquisition, the head may be tilted such that there is a significant roll angle. In such cases, the symmetry about the midsagittal plane is violated and 2-D symmetry based comparison of hyperintense regions may be erroneous.
When an artifact is symmetric to an infarct region, there may be a chance of losing an infarct region. In addition, there are various reasons because of which artifacts may be retained by the present embodiments, including:
Nevertheless, it is clear from the experimental results that the embodiments provide a fast and pragmatic approach for artifact removal. They do not utilize anatomy related information whose processing may be challenging and more time consuming. The present approach also does not take into account the location of infarct, which is critical [23-24], and influences the outcome.
The embodiments can be applied as a stand alone post processor or be a part of a stroke CAD system, and can provide a fast post-processing tool to reduce artifacts in infarct processing applications. In fact, the processing time for the present embodiments, when implemented on a matlab platform, was under 1 second.
Filing Document | Filing Date | Country | Kind | 371c Date |
---|---|---|---|---|
PCT/SG09/00009 | 1/6/2009 | WO | 00 | 7/9/2010 |
Number | Date | Country | |
---|---|---|---|
61020187 | Jan 2008 | US |