WO2004104617A2

WO2004104617A2 - Method of providing noise cancellation in fmri data

Info

Publication number: WO2004104617A2
Application number: PCT/US2004/016312
Authority: WO
Inventors: Anders M. Dale; Douglas Nowlin Greve
Original assignee: The General Hospital Corporation
Priority date: 2003-05-21
Filing date: 2004-05-21
Publication date: 2004-12-02
Also published as: WO2004104617A3

Abstract

Methods of performing noise cancellation are presented. One method includes acquiring a plurality of data sets, analyzing a subset of the plurality of data sets, generating a filter from the subset of data sets, and applying the filter to a data set which results in a noise cancelled data set. Alternately a regressor is generated instead of a filter, the regressor is applied to the model, and then the data set is applied to the model. Another method includes applying an excitation pulse to a subject, acquiring a baseline data set before a data set induced by the excitation pulse has developed, waiting a predetermined period of time, acquiring the data set induced by the excitation pulse, generating a regressor from the baseline data set, applying the regressor to a model, and applying a dataset to the model to obtain a noise cancelled data set.

Description

TITLE OF THE INVENTION

Method of Providing Noise Cancellation In fMRI Data

FIELD OF THE INVENTION The present invention relates generally to functional Magnet Resonance Imaging

(fMRI) analysis and more particularly to a method for applying filtering to a preprocessing stage of an fMRI analysis stream thereby resulting in scans that are more reliable.

BACKGROUND OF THE INVENTION Functional MRI (fMRI) is known to those of ordinary skill in the art. fMRI is a technique which can be used to image brain activity related to a specific task or sensory process. fMRI is based upon increase in blood flow to local vasculature which accompanies neural activity in the brain. This increase in blood flow results in a corresponding local reduction in deoxyhemoglobin since the increase in blood flow occurs without an increase of similar magnitude in oxygen extraction. Since deoxyhemoglobin is paramagnetic, it alters the magnetic resonance image signal. Using an appropriate imaging sequence, human cortical functions can be observed on a clinical scanner without the use of exogenous contrast enhancing agents.

One advantage of fMRI is that it is not necessary to inject radioactive isotopes as is required by certain other imaging techniques. Additionally, another advantage of fMRI imaging is that the total scan time required can be relatively short.

An fMRI scan produces an image comprising a plurality of voxels. A voxel is a three-dimensional element which subdivides a volume of space and can be represented using a three-dimensional coordinate system (e.g., x, y, z coordinates) and an intensity value. Spatially correlated noise in fMRI is indicated by two or more voxels whose non- task-related temporal waveforms have high cross-correlation. This phenomenon has been observed in the study of functional connectivity and resting-state networks.

Lack of measurement reliability or repeatability is in fMRI is attributable, at least in part, to variance induced by noise. While there are many sources of noise in fMRI, they can generally be divided into two categories: measurement noise and physiological noise. The presence of measurement noise is in a large part due to thermal noise generated in the electronics used to detect a very weak electrical signal (it is also known as Johnson noise). The presence of physiological noise results from an interaction between the physiology of the subject and the way that the fMRI signal is measured. Specifically, the fMRI blood oxygenation level dependent (BOLD) contrast can be contaminated by changes in heart rate, respiration, motion, and other effects unrelated to changes in the neural firing of the brain itself.

Measurement noise tends to be spatially and temporally uncorrelated (white), whereas physiological noise tends to be spatially and temporally correlated (colored). The effect of physiological noise relative to measurement noise increases as the fMRI magnetic field strength increases or as the sensitivity of the receiver coil improves, both of which are currently strong trends in fMRI. Measurement noise, however, plays a role in the inability of fMRI to reliably detect activation in susceptibility regions. Susceptibility regions result from a strong inhomogeneity within a voxel, generally due to the presence of a nearby air cavity (e.g., sinuses). This inhomogeneity causes the signal in those regions to decay very rapidly such that the signals are swamped by the measurement noise. As a result, fMRI is essentially blind in those regions. Along the same lines, fMRI data directly from the scanner is usually converted into an image using an inverse Fourier transform. The assumptions that make this operation work often do not hold very well, resulting in intensity and metric distortion. Intensity distortion results in the image being too bright or not bright enough. Metric distortion results in the image being stretched or warped, for example. Changes in these intensity and metric distortions over time can result in noise.

Several techniques have been developed to make fMRI measurements more reliable. A first technique is known as spatial smoothing. In the spatial smoothing technique, values of neighboring voxels are averaged together. The averaging reduces spatially uncorrelated (e.g., thermal) noise but does little to reduced spatially correlated noise and can actually enhance such noise. The spatial smoothing technique also has the effect of reducing the spatial resolution. Another technique used to make fMRI measurements more reliable is known as temporal smoothing. In the temporal smoothing technique, neighboring time points are averaged together. Again, this reduces the variability due to temporally uncorrelated noise but does little to reduced temporally correlated noise. This also has the effect of reducing temporal resolution.

Still another technique is known as temporal whitening. In the temporal whitening technique, the temporal correlations are modeled and filtered out thus reducing the overall variance. While this can be effective, the technique models the temporal noise as a stationary process, which is usually not the case in fMRI. Also, the noise covariance structure must be estimated from a very limited amount of data.

It would, therefore, be desirable to provide a method for filtering fMRI data to provide more reliable images.

SUMMARY OF THE INVENTION

A method for providing noise cancellation in measurements is presented. The method includes acquiring a plurality of data sets, analyzing a subset of the plurality of data sets, generating a filter from the subset of data sets, and applying the filter to a data set to provide a noise cancelled data set.

Another embodiment includes acquiring a plurality of data sets, analyzing the data sets, generating a regressor from the data sets, applying the regressor to a model, and applying a dataset to the model which results in a noise cancelled data set.

Another method for providing noise cancellation in measurements includes applying an excitation pulse to a subject, acquiring a baseline data set before a data set induced by the excitation pulse has developed, waiting a predetermined period of time, acquiring a data set induced by the excitation pulse, generating a regressor from the baseline data set, applying the regressor to a model, and applying a dataset to the model to obtain a noise cancelled data set. A further embodiment includes applying an excitation pulse to a subject, acquiring a baseline data set before a data set induced by the excitation pulse has developed, waiting a predetermined period of time, acquiring a data set induced by applying an excitation pulse, and applying the baseline dataset and the dataset to a model which results in a noise cancelled data set. Still another embodiment for providing noise cancellation in measurements involves acquiring a plurality of data sets from a single excitation and applying the data to a model. The reconstruction of the image is constrained to account for intensity distortion and metric distortion to allow for recovery of the signal from susceptibility regions as well as a general reduction in measurement noise.

A system for providing noise cancellation in fMRI measurements includes an MRI scanner coupled to an MRI control and storage system. An image processing system is coupled to the MRI control and storage system, as is a user interface. The system further includes a noise cancellation system coupled to the image processing system and the MRI control and storage system. The noise cancellation system includes either a filter processor, a regression processor or both.

BRIEF DESCRIPTION OF THE DRAWINGS

The foregoing features of this invention, as well as the invention itself, may be more fully understood from the following description of the drawings in which:

Figure 1 is a block diagram of the system used in performing spatial noise cancellation; Figure 2 is a flow diagram of a first embodiment of spatial noise cancellation;

Figure 3 is a state diagram showing the interconnection of various states of the process of the present invention;

Figure 4 is a flow diagram of a second embodiment of performing spatial noise cancellation; Figure 5 is a graph of variance in rest runs of four different subjects;

Figure 6 is a graph showing the effect of spatial noise cancellation on a temporal power spectrum averaged over in-brain voxels;

Figure 7 is a brain scan showing a reduction in residual error variance due to spatial noise cancellation; Figure 8 is a brain scan showing a percent increase in signal to noise ratio due to spatial noise cancellation; and

Figure 9 is a flow diagram of a method for performing short-post excitation delay regression analysis.

DETAILED DESCRIPTION OF THE INVENTION

Referring now to Figure 1, a magnetic resonance imaging (MRI) system 10 utilized to provide noise cancellation of functional MRI (fMRI) data is shown. The MRI system 10 includes an MRI scanner 12 which includes a magnet having gradient coils and radio frequency (RF) coils disposed thereabout in a particular manner to provide a magnet system. In response to control signals provided from a controller processor, a transmitter provides a transmit signal to the RF coil through an RF power amplifier. A gradient amplifier provides a signal to the gradient coils also in response to signals provided by a control system 14.

The magnet system is driven by the transmitter and amplifiers. The transmitter generates a steady magnetic field and the gradient amplifier provides a magnetic field gradient which may have an arbitrary direction. For generating a uniform, steady magnetic field required for MRI, the magnet system may be provided having a resistance or superconducting coils and which are driven by a generator. The magnetic fields are generated in an examination or scanning space or region in which the object to be examined »is disposed. For example, if the object is a person or patient to be examined, the person or portion of the person to be examined is disposed in the scanning region. The transmitter / amplifier drive the coil. After activation of the transmitter coil, spin resonance signals are generated in the object situated in the examination space, which signals are detected and are applied to a receiver. Depending upon the measuring technique to be executed, the same coil can be used for the transmitter coil and the receiver coil or use can be made of separate coils for transmission and reception. The detected resonance signals are sampled, digitized in a digitizer. Digitizer converts the analog signals to a stream of digital bits which represent the measured data and provides the bit stream to the control system 14 for storage and/or further processing.

The control system 14 processes the resonance signals measured so as to obtain an image of the excited part of the object. An image processing system 16 is coupled to the control system 14 and can be used to display the reconstructed image. The display may be provided for example as a monitor or a terminal, such as a CRT or flat panel display. A user provides scan and display operation commands and parameters to the control system 14 through an interface 18 which provides means for a user to interface with and control the operating parameters of the MRI system 10 in a manner well known to those of ordinary skill in the art. The control system 14 is used, among other things, to specify magnetic field gradient directions for each scan and to generate images for each magnetic field gradient direction. The control system 14 also has coupled thereto a noise cancellation system 19 which includes a filter processor and a regression reconstruction processor. Each of the components depicted in FIG. 1, except for the noise cancellation system 19, are standard equipment in commercially available MRI systems. It should be appreciated that the MRI system must be capable of acquiring the data which can be used by the noise cancellation system 19 in the manner to be described hereinbelow. fMRI waveforms from spatially disparate brain regions can have very high temporal correlations even when the subject is not performing any task. While the waveforms themselves change over time, the spatial pattern of voxels that are temporally correlated tends to remain constant. This phenomenon has been observed in the study of functional connectivity and resting-state networks.

Those techniques seek to explain brain function through the study of spatial correlations; the present invention exploits these spatial correlations to improve the estimation and detection of fMRI activation, thereby providing more reliable scans. If the noise at one voxel is significantly correlated with the noise at a second voxel, and if this correlation is reliable over time, then the waveform at the second voxel can be used to predict the noise at the first voxel (and vice versa). The prediction can then be subtracted from the first voxel to result in a noise waveform with a lower variance. This subtraction process is referred to as noise cancellation. Lower noise variance makes the estimation and detection of activation more robust. The overall effect of the noise cancellation is to spatially whiten the noise, thereby providing more reliable scans.

The present invention makes blood oxygenation level dependent (BOLD) functional magnetic resonance imaging (fMRI) measurements of the brain more reliable and repeatable by reducing spatially and temporally correlated noise. The present invention further serves to reduce the intensity and metric distortion inherent in fMRI, especially from brain regions typically lost due to susceptibility artifacts.

The present invention achieves the reduction in spatial and temporal correlated noise by use of at least one of the following techniques. A first technique is referred to as Spatial Noise Cancellation (SNC), a second technique is referred to as Short Post- Excitation Delay (SPED) regression analysis, and a third technique is referred to as Physics-Constrained Reconstruction (PCR) or Time Domain Reconstruction (TDR).

The SNC technique will now be described in detail. In Spatial Noise Cancellation, the actual noise waveforms are estimated directly from the data by exploiting the fact that noise tends to form repeatable spatial patterns. Once estimated, the noise waveforms are removed (cancelled) from the data. While SNC is specifically designed to remove spatially correlated noise, this technique also has the effect of removing temporally correlated noise because the two tend to be inseparable. An advantage associated with the SNC technique is that it will work even if the noise is temporally non-stationary. SNC does not require any changes in the way the fMRI data are acquired. SNC removes long- range spatial correlations from fMRI data. This spatial whitening also reduces temporally non-stationary fMRI noise. SNC is a type of spatial filtering and is applied in the preprocessing stage of the fMRI analysis stream.

There are several hurdles that must be overcome in order to implement the SNC method. First, it would be desirable to use all the voxels to predict the waveform at each voxel. On the surface, the complexity of this operation appears to be square in the number of voxels, which would make it too computationally intense for practical implementation using current state of the art processing devices and techniques. Also, because the number of voxels is much greater than the number of time points, the computation of the filter is grossly underdetermined. Second, in order to avoid circularities, the spatial pattern of weights that are used to compute the predicted noise waveform are preferably derived from a data set that is independent of the data set to which cancellation will be applied. This requires that the spatial noise patterns be repeatable.

Third, in runs where subjects are performing a task, some voxels will have task- related signal which will be spatially correlated, and this signal spatial correlation should not be represented in the cancellation filter (i.e., the filter should be derived from data without signal).

Fourth, even if the filter is derived from signal-free data, the signal and noise may still be spatially correlated (i.e., the filter derived from signal-free data may look very similar to a filter derived from data with signal). The implication of this is that one would be subtracting one voxel with signal from another voxel with signal which will result in a loss of signal. If this signal loss is greater than the noise loss, then there will be a net loss in the signal-to-noise ratio, making it more difficult to detect activation. Additionally, it would be desirable to provide a methodology in which there is no need to collect special data (e.g., rest runs) or to throw data away (e.g., to avoid circularity).

To derive the filter, a singular value decomposition (SVD) of an fMRI data set is used. The SVD may be represented mathematically as shown in Equation (1) below:

D = U-S-V¹ Equation (1)

where D is a data matrix, U is a matrix whose columns are temporal eigenvectors of D, S is a diagonal matrix of singular values, and V is a matrix whose columns are the spatial eigenvectors of D.

The spatial pattern of noise correlations is embedded in V. Use of the SVD is equivalent to principal component analysis (PC A). The data matrix D has a number of rows equal to the number of time points (Nt_p) and a number of columns equal to the number of voxels (N_v). Under these conditions, the U matrix is N_tp rows by N_tp columns (which may be expressed as Ntp x N_tp), S is Ntp x N_v, and V is t x N_v. Simply computing the SVD would be computationally intractable because V is so large. However, only the components of V that actually account for variance in D need to be considered. The amount of variance spanned by an eigenvector is equal to the singular value that corresponds to that eigenvector. In this case, only Np of the singular values are non-zero, so the relevant dimensionality of V is only N_v x tp which is computable. Further, as described below, only the first 20 or 30 spatial eigenvectors are of interest.

The prediction of the noise N_j is computed by projecting the raw data into and back out of the spatial eigenvector space. This may be accomplished using Equation (2)

N_J = (D_J - V_k) - Vl Equation (2) where D_j is a set of raw data values taken from a data set j, and V_k are the spatial eigenvectors derived from an independent data set k. The parentheses assure that no intermediate computations will create a matrix having N_v rows and N_v columns (i.e. an N_v x N_v matrix). This noise estimate is subtracted (cancelled) from the original data set as shown in Equation (3): Dj = D_J -N =D_J -(D_J -V_k)-V_k ^τ Equation (3)

It should be noted here that the spatial noise cancellation works by orthogonalizing the raw data with respect to the noise spatial eigenvectors; i.e., Dj is in the null space of V . This operation can also be interpreted as a spatial filtering stage as expressed by Equation (4):

Dj = D_}F_k,F_k = (I-V_kVl) Equation (4)

where F is a spatial filter having Ν_v columns and N_v rows. It is not practical for computational reasons to compute the spatial filter F explicitly. This part of the derivation is only meant to demonstrate that SNC is formally equivalent to a spatial filtering operation (indeed, it is applied at the same point in the functional analysis stream that one would apply spatial smoothing).

The SNC filter has been derived directly from data set D, implicitly assuming that D contained no signal. In what follows, the use of data sets with task-related signal are described, as well as the computation of independent filters without sacrificing any of the data.

Flow charts of the presently disclosed methods are depicted in Figures 1, 3, 8 and 9. The rectangular elements are herein denoted "processing blocks" and represent computer software instructions or groups of instructions which may be executed by a computing device (e.g. a personal computer, a general purpose computer or a processor). The diamond shaped elements, are herein denoted "decision blocks," and represent computer software instructions, or groups of instructions which affect the execution of the computer software instructions represented by the processing blocks. Alternatively, the processing and decision blocks represent processing performed by functionally equivalent circuits such as a digital signal processor circuit or an application specific integrated circuit (ASIC). It should be appreciated that the flow diagrams do not depict the syntax of any particular programming language. Rather, the flow diagrams illustrate the functional information one of ordinary skill in the art requires to fabricate circuits or to generate computer software to perform the desired processing. It should be noted that many routine program elements, such as initialization of loops and variables and the use of temporary variables are not shown.

Referring now to Figures 2 and 3, a flow diagram 20 (Figure 2) and a state diagram 32 (Figure 3) relating to a first process for performing SNC are shown. By way of example only, four hypothetical data sets 34a, 34b, 34c, 34d (also referred to herein as task runs) are described, though it should be understood that any number of task runs could be used. The process of performing spatial noise cancellation begins at processing block 22 in which data sets (task runs) are acquired. The data sets comprise run A (labeled 34a in Figure 3), run B (labeled 34b in Figure 3), run C (labeled 34c in Figure 3) and run D (labeled 34d in Figure 3).

Once the data sets are acquired, processing proceeds to processing block 24 in which a subset of the task runs 34a-34d are analyzed. In the example of Figure 3, task runs 34b-34d (i.e. runs B, C, and D) are segregated from run A and are analyzed separately. In the process of analysis (labeled 36 in Figure 3) with a General Linear Model (GLM), the raw waveforms are separated into signal waveforms and residual waveforms. A General Linear Model (GLM) is a model of the observed signal changes over time. In particular, the observed signal is hypothesized to be a weighted (linear) sum of task and nuisance regression vectors. In one embodiment, the weights are determined emperically by fitting the model to the data using a least-means-square (LMS) algorithm. Other techniques (including other empirical techniques as well as analytical techniques), however, can also be used to determine the weights. The regression vectors define the signal subspace. Everything not fit by the GLM is the residual error. Residual waveforms 38b, 38c, 38d are provided for each data run (also labeled in Figure 3 as E_b , E_c , and E_d)

Additionally, as shown in Figure 3, an omnibus activation map 40 (i.e. a map of voxels activated by any component of the task) is computed and is later used to remove voxels with signal from the computation of the estimation of the spatial noise. This map is also referred to as a "Projection Mask", and provides a manner for compensating for spatial correlations in the signal and noise. It should be noted that all voxels receive cancellation regardless of the projection mask because the projection mask only effects which voxels are used as the source of cancellation.

Processing then proceeds to processing block 26 in which a filter is generated from the subsets. In the exemplary state diagram of Figure 3, the filter is labeled as F_bcd associated with reference number 42 and is computed from the set of residuals 38b-38d, which reduces the influence of the signal on the filter 42.

Processing then proceeds to processing block 28 in which the filter is applied to a data set which was not used to generate the filter. In the exemplary state diagram of Figure 3, the filter is applied to run A (reference numeral 34a) and results in a spatial noise cancelled data set A 44. This process is repeated successively leaving out one run in each repetition (this technique is sometimes referred to as "jackknifmg"). That is, filter F_acd is derived from data sets A, C and D and is applied to run B to get spatial noise cancelled data set B ; filter F_abd is derived from data sets A, B and D and is applied to run C to get spatial noise cancelled data set C ; and filter F_abC is derived from data sets A, B and C and is applied to run D to get spatial noise cancelled data set D . In such a manner, an independent filter is applied to each run and all runs are used.

In processing block 30, the filtered data set is analyzed. In the exemplary embodiment of Figure 3, the analysis is represented by reference numeral 46. In this

example, the four spatial noise cancelled data sets (A,B,C,andD) are processed processed using any technique for the detection of activation, such as a GLM.

An alternative method of performing spatial noise cancellation comes from the fact that the columns of a matrix Rj (which may be computed as Dj ^* V_k) are a set of temporal vectors which will span the temporal space of N, . These vectors are used as regressors in

the general linear model (GLM) instead of subtracting N_} from Dj. There is some degree of circularity here because the regressors used to model a data set are derived from the data set.

Referring now to Figures 3 and 4 an alternate embodiment for performing spatial noise cancellation (SΝC) is shown. A flow diagram 60 of the processing to perform spatial noise cancellation begins at processing block 62 in which data sets are acquired.

The data sets comprise comprise run A (labeled 34a in Figure 3), run B (labeled 34b in

Figure 3), run C (labeled 34c in Figure 3) and run D (labeled 34d in Figure 3).

In processing block 64 the task runs are analyzed. In the process of analysis with the GLM, the raw waveforms are separated into signal waveforms and residual waveforms. The residual waveforms are provided for each data run and are labeled E_b (38b in Figure 3), E_c (38c in Figure 3), and E_d(38d in Figure 3). Not shown is the residual EA for data run A 34a.

In processing block 66 regressors are generated from the residual waveforms 38b, 38c, 38d. The regressors will be applied to the GLM.

In processing block 68 the regressors are applied to the GLM. In this example, the regressors derived from data sets A, B, C, and D are applied to the GLM.

In processing block 70 a dataset (data set A in the example shown in Figure 3) is applied to the model including the regressors generated from datasets A, B, C, and D. This results in filtered data set A (denoted by reference number 44 in Figure 3).

In processing block 72, the filtered data set is analyzed. In this example, the spatial noise cancelled data set A is processed in the normal way.

This process is repeated successively for each data run. That is, run B is applied to the GLM having the regressors to get spatial noise cancelled data set B ; run C is applied

to the GLM having the regressors to get spatial noise cancelled data set C ; and run D is applied to the GLM having the regressors to get spatial noise cancelled data set D .

To test the SNC methodology, several fMRI data sets were acquired using several scanners. The main data set (subject DG) used for testing the method was acquired using a Siemens 3.OT Allegra scanner (TR = 2 sec, TE = 30 ms, 3.125 mm in-plane, 6 mm slice thickness, including 1 mm skip, with 23 transaxial slices). Twelve runs (180 time points per run) were acquired, with the six odd runs being task runs and the six even runs being rest runs. During the task runs, the subject performed an event-related visual semantic association task. Briefly, the subject was shown either 2 or 4 probe words followed by a single target word. The subject then had to respond (with a key press) as to which probe word was most related to the target word. The entire trial lasted only 3 seconds. The probe-target pairs were either loosely or highly related for a total of 4 event types. Seventy- two such trials were presented with random order and stimulus onset asynchrony for each task run. During rest runs the subject as instructed to lay still with his/her eyes closed. Several other rest-only data sets were collected. These were subjects JL and SA (GE 3.OT Signa scanner, voxel size 3.125 X 3.125 X 8 mm³, Nψ= 128), subject NH (Siemens 3.OT Allegra, voxel size 3.125 X 3.125 X 5 mm³, N_φ = 116), and subject AD (Siemens 3.OT Allegra, voxel size 2.7 X 2.7 X 3 mm³, N_tp= 128). For all data sets, the TR was 2.0 sec and TE was 30 ins. All data sets were motion corrected by applying correction software (e.g., Analysis of Functional Neuro Images software) to the first image of the first run.

The method was also applied to rest data in order to demonstrate the repeatability of the spatial noise patterns and to document its temporal effects. The first step in the evaluation of this method is to verify that the spatial patterns are repeatable in rest data, otherwise the filter computed from one data set will be inappropriate for an independent data set. To do this the amount of whole-brain variance in run j spanned by each of the spatial eigenvectors computed from all runs excluding j (this is a jack-knife cross- validation scheme) was computed. This was repeated for all runs, and then the variance spanned was then averaged across all runs and rescaled to a percentage of whole-brain variance to give a measure of repeatability. More formally, the percent variance spanned (PVS) by the n spatial eigenvector is given by

PVS(n) = Equation (5)

2 where _(jD is the variance of run j measured across all time points and brain voxels and

_j is the same after applying SNC as described by Equation 3 but using only the n^to column of V_k. Any low-frequency fluctuations in the data were removed using a 5 order polynomial prior to computation or application of the SNC filter.

The cumulative percent variance spanned (CPVS) is the amount of variance spanned by the first through the n spatial eigenvectors. The amount of variance spanned is also a good measure because it equals the reduction expected in the residual error variance when SNC is applied to task data, and so directly relates to the expected improvement in statistical power since the SNC operation does not remove any signal.

Referring now to Figure 5, the CPVS curves 76-82 are shown for subjects JL (labeled 76), SA (labeled 78), NH (labeled 80), and AD (labeled 82). The results show that the dimensionality of the noise is very low with as few as 10 spatial components removing as much as 25-30% of whole-brain noise. The effect is quite consistent across subject and scanner manufacturer. The variance explained for subject AD 82 is lower because the voxel size is smaller making the relative contribution of spatially white (instrument) noise much greater. The implication of this is that the relative contribution of spatially correlated noise can be expected to grow with higher field strength and better coils. The CPVS appears to asymptote after about 20-30 eigenvectors. Actually, the CPVS curve will slowly rise and reach 100% only after the number of eigenvectors reaches the number of brain voxels.

This exercise helps resolve another question about how to compute the SNC filter, namely how many spatial eigenvectors (i.e., columns of Vk) should be including in the computation of the filter. The answer is to set the number of spatial eigenvectors at the number of EVs corresponding to the knee of CPVS curve. Referring now to Figure 6, the effect of SNC on the power spectrum is shown for subject DG. The curves were computed by averaging the power spectrums at each voxel before ("Brain" 84) or after ("Brain-SNC-20" 86) the SNC filter computed from the first 20 eigenvectors was applied. The curve labeled "Air" 88 was computed in the same fashion across all voxels outside of the head; this gives a measure of the noise floor. The power spectrum of the unfiltered in-brain voxels has the I/f shape typically found in fMRI waveforms; it also has large fluctuations across frequency. The power spectrum after SNC is different than that of the raw data in three important respects. First, the noise power is lower across all frequencies; this indicates that the total noise power is reduced as would be expected from the CPVS results above. Second, it has a much smaller low-frequency component. Third, the overall spectrum is much smoother. These last two points indicate that purely spatial filtering of SNC has the effect of temporally whitening the noise.

The SNC methodology was also applied to task data to demonstrate its efficacy. SNC was applied to a data set in which the subject was performing an event-related semantic association task as described above. The data were analyzed with and without SNC using a GLM in which the shape of the hemodynamic response was assumed to be a gamma function; low frequency drift was removed by including 5 order polynomial regressors in the GLM (separate regressors for each run). The SNC filter had an order of 20. All data were spatially smoothed (FWHIM=5mm). It took less than 1 hour of CPU time (1GHz Pentium III) to compute and apply the SNC filter for the entire data set.

The effect of SNC on residual error variance was evaluated by computing the percentage by which it was reduced relative to the no-SNC analysis. The result for a single slice is shown in Figure 7. The reduction is mostly in and around the cortical gray matter, where one would, for physiological reasons, expect spatially correlated noise to exist. This is also where it is most important to reduce the noise. The reduction can exceed 50% in many voxels.

The effect of SNC on the contrast-to-noise ratio (CNR) was evaluated by examining the omnibus test F-ratio of voxels that were significant at a p<.0001 level under either analysis method. This was done to account for voxels which had no signal present. For the voxels which had a signal present, the percent increase in the CNR (i.e., omnibus F-ratio) was computed. The result (Figure 8) is an increase in CNR by as much as 50% in the primary visual areas. An increase of 50% is quite dramatic as it is equivalent to presenting 50% more stimuli (i.e., scanning for 50% longer). The CNR dropped in some areas, mainly in white matter and around the ventricles, though other areas fMRI noise is highly spatially correlated. As few as 10 spatial principal components can account for as much as 30% of rest noise. A spatial noise cancellation techniques was introduced to estimate and remove the spatial correlations across all brain voxels. The SNC technique is computationally tractable and can be applied to task data without collecting additional scans. The SNC technique does not make any assumptions about the shape of the spatial correlation function. SNC reduces residual error variance by 50% in cortical gray matter. Spatial Noise is probably in legitimate cortical gray matter. Cancellation can increase contrast-to-noise ratio by 50% in primary visual areas. Another method used to render fMRI scans more reliable is referred to as the SPED method. In the SPED method the data acquisition is changed. Typically, MRI data is collected by applying an RF excitation pulse and then measuring the signal induced by the excitation. In fMRI, the signal is typically sampled over a time range that brackets the time at which one would expect the maximum BOLD contrast to develop. Prior to that time, the scanner is (for the most part) idle. SPED utilizes information collected immediately after the RF excitation pulse to characterize noise. Data from this time period is especially good for characterizing the noise because the delay at which the data are collected is so short that no BOLD-related contrast has had a chance to develop, and so any variation during the SPED time is mainly related to noise. The SPED data does not need to be an entire image.

The SPED data is decomposed using Principal Component Analysis, or PCA (also known as singular value decomposition), in order to extract the most relevant temporal principal components. These components are then used as regressors in the standard general linear model based analysis of fMRI.

Referring now to Figure 9, a process 90 for performing SPED to provide noise cancellation is shown. In processing block 92 an excitation pulse is applied to a subject. fMRI systems typically utilize an RF pulse as the excitation pulse.

In processing block 94 a baseline data set is acquired immediately after the excitation pulse is supplied. This data set is acquired before any data that is induced by the excitation pulse has had a chance to develop.

In processing block 96 data collection is performed after waiting a predetermined period of time. This predetermined period of time is required to allow the signals induced by the excitation pulse to develop.

In processing block 98 a data set is acquired. This data set is a result of the reaction to the excitation pulse.

In processing block 100 a regressor is generated from the baseline data set. This regressor is used to cancel the effects of noise in the data.

In processing block 102 the regressor is applied to the model. The baseline data set is comprised primarily of noise, and thus is applied to the model as a regressor in order to cancel the effects of the noise from the data acquired in response to the excitation signal.

In processing block 104 the dataset acquired in response to the excitation pulse is applied to the model and regressor. This results in the provision a noise-cancelled data set..

Yet another method to render fMRI scans more reliable is referred to as the Physics Constrained Reconstruction (PCR) technique. PCR removes the metric and intensity distortion from fMRI signals while additionally attenuating the white noise and recovering the signal from the susceptibility regions. Similar to the SPED technique, the PCR technique utilizes information collected immediately after the RF excitation pulse to characterize noise. PCR reduces intensity and metric distortion by constraining the reconstruction of the image in such a way so as to assure that the underlying signal evolved according to physics-based models of the fMRI environment. This allows recovery of the fMRI signal from susceptibility regions as well as a general reduction in measurement noise. In the PCR method, physical models of how the signal will evolve over the post- excitation period are utilized in order to reconstruct the image with less intensity and metric distortion. Currently, the most prominent method of reconstruction assumes that the collected data is a perfect Fourier transform of the image and so applies an inverse Fourier transform in order to reconstruct the image. Some of the implicit assumptions in the inverse Fourier reconstruction are that the signal at each voxel does not change over the course of acquisition, that k-space is sampled uniformly, and that the phase developed at a point in k-space is only due to contrast in the image. In fMRI, these assumptions are violated which causes intensity and metric distortions, however there are models of how these violations will affect the post-excitation signal. In PCR, data collected acquired across the entire range of post-excitation delays and is used to fit this more complicated model. Note that the data acquired for PCR can also be used for SPED.

Image reconstruction in MRI is usually implemented by arranging the observed data into a k-space image and then applying an inverse FFT. There are two assumptions underlying this method: 1) the intensity of the underlying image does not change during readout; and 2) each value during readout is the result of a convolution of the intensity image with a spatial kernel that is consistent with the FFT. If these conditions are not met, then the reconstructed image will be distorted, both in space and intensity, with respect to the true underlying intensity image.

In conventional MRI, the first assumption is usually met because the readout time is short (i.e., a few milliseconds) relative to the decay of the underlying intensity image. The second condition may not be met due to nonlinearities in the gradients. In fast MR imaging (e.g., EPI and spiral), the readout duration is much longer (e.g., 30-50 ms). During this time, the intensity of the image can drop substantially (by 70% or more). This can cause blurring in the phase encode direction when a simple inverse FFT is used to reconstruct the image. The spatial convolution kernel can also deviate substantially from that of an FFT due to long readout in the presence of B0 (the static magnetic field plus gradients) field inhomogeneity (deviation of the magnetic field from the average value of the field). B0 inhomogeneity causes the resonant frequency to be slightly off, which, over the course of a long readout, can result in substantial deviation of the spin from its ideal phase. This results in substantial warping of the EPI images. BO gradients within a functional voxel (due, e.g., to local susceptibility effects) can produce a rapid de-phasing of the spins within the voxel, leading to signal "drop-outs". A technique for addressing these problems is presented. The technique uses the full time-domain model of the readout signal, which is given by

where y(τ) is the (complex) readout signal, τ is the time since the last RF excitation (i.e., the post-excitation delay, or PED), y(τ,v) is the contribution from voxel v which spans volume element R_v over which the integral is computed. p(r) represents the transverse magnetization immediately after the RF-pulse at location r , g(τ) is the vector gradient moment, R^(r) is the relaxation rate, and Δω₀(r) is the deviation of the resonance frequency from the ideal due to inhomogeneities in the BO field. Without the relaxation or BO effects, the above equation reduces to a simple FFT of p(r) . This is approximately the case when the readout is short (i.e., the interval of τ is small).

In the case when the interval of τ is not small, the relaxation and BO effects are taken into account in order to reconstruct the image without metric and intensity distortion. In order to do this, the above equation is simplified by assuming that R , p , and g are constant within the voxel element. This results in the following equation

y(τ) = ∑ p(v)e"^8~ (e^{' Mτ} (7)

All the effects of relaxation and B0 are encapsulated in the part of the equation in parentheses. This is referred to as the decay function: d(τ,v) - e-^R*^{v)τ ^e-^ ^r)Tdr (8)

Again, without the decay portion, equation 7 is just the FFT. If the decay function is known, then the transverse magnetization image p(r) can be reconstructed without distortion. The decay function is measured explicitly using a multi-echo FLASH sequence in which the same line of k-space is read out repeatedly (e.g., 64 echoes) for each excitation. When all the lines corresponding to the nth echo are assembled into a k-space image, the inverse FFT gives an estimate of the free induction decay (FID) at nth echo time for each voxel. This data should be acquired with the same spatial resolution and slice prescription as the functional scans. Once this decay function has been mapped for each voxel, the estimation of the image ( p(r) ) reduces to a straightforward linear one, using the appropriate linear basis functions.

The method for reconstructing the p(r) image without the distorting influences of relaxation and BO inhomogeneities has been described. However, in functional BOLD imaging, maps of changes in relaxation rates R^(r) as a function of time should be obtained. Accordingly, the above method is extended to reconstruct the Δ/^ image. This is shown in the equation below. _e-RX^,»r _{= e}-(ϊξ(v⁾-AR ^(v,_l))τ _ _e-^^(») Δ^(^». ⁾r „ _e M* (\ - £jζ(V,t)τ) (9) where t is the shot-to-shot time (not to be confused with τ which is the time since the last shot), and R^ has been separated into a component that does not vary with time t (R^(v) , the mean R^) and a part that varies across time t (Δ ζ(v,t)). Finally, if ζ(v,t) is small, then e^{~ 2 v}'^{, r} i_s approximately I- AR^(v,t)τ . If equation 9 is substituted into equation 7, the readout equation can be separated into a time-invariant component associated with the average relaxation time and a time-varying component due to fluctuations in relaxation rate:

y(τ,t) = y_t (τ)-y^(τ,t) (10) where

and y_j, (τ) is the readout signal created by the average relaxation rate; this is what is

Ri measured by the multi-echo FLASH calibration scan. Thus, when the same reconstruction algorithm is applied that was used for equation 7, an undistorted image of the product of the magnetization transfer and the change in relaxation rate are provided. Note that the decay function has changed: ^d,Aτ) = (™ ^■Rι( )τ -iAω₀(r)τ flf^« dr) (12)

This is the same as equation 8 multiplied by τ and is dependent upon the average relaxation rate, which is what will be measured by the multi-echo FLASH calibration scan described above.

Note that this approach applies to any arbitrary k-space traversal g(τ) , including traversals that redundantly sample k-space (e.g., multi-echo EPI acquisitions). In most fMRI studies using EPI acquisition, only a single echo is acquired, using an echo-time that represents a compromise between BOLD contrast-to-noise (which is optimal at TE=T₂* of blood, approximately 50ms at 1.5T) and drop-outs due to susceptibility gradients (reduced by lowering the TE). By collecting multiple echoes, starting immediately after the RF excitation, and using the time-domain model above, estimates can be obtained for R^(v) which optimally combine the measurements at all time delays ( τ ) for every voxel, thus greatly improving the contrast-to-noise of the resulting functional activation maps.

To show this, equation 7 is rewritten in matrix format: y = (F » D)p = Xp, X = F * D (13) where y is the vector of readout values where each row represents a different readout time (i.e., τ ) which also corresponds to a different point in k-space (i.e., g(τ)), p is the vector of image intensities, F is a matrix in which each row is the spatial FFT kernel for the point in k-space indicated by that row, D is a matrix where each column is the decay function (equation 8) for the corresponding voxel in p , and • indicates an element by element multiplication. This linear equation can be solved for the image intensities:

p = (X'Xy^lX'y (14)

Note that when the k-space traversal is redundant (i.e., there are more observations (rows in y) than voxels (rows in p )) and the noise in y is white, then this will yield the optimal unbiased estimate of the image intensity. Equation 11 can also be recast into matrix format:

where AR^ is the vector of changes in relaxation rate, and D . is the matrix where each column is the decay function (equation 12) for the corresponding voxel in p . The optimal estimate ofz is then given by:

z = (X' X X-^lX' y . (16)

Note that this estimate is optimal under the white noise assumption and takes into account the fact that the BOLD contrast evolves as τe^~Rιτ .

This formulation can also be used to determine an optimal redundant k-space traversal (i.e., a traversal that will result in the lowest variance of z ). The mean-square error (MSE) in z is given by MSE - tracedX^X^γ¹ ^ . = aσ) . (17)

— —Vi ^JάR2 <U!2 where a is the variance reduction factor; this is the amount that the noise inherent in the measurement will be reduced by redundant measurements. Thus, the best estimate of z arises when a is smallest. Note that or only depends upon F and the decay function. F is fixed by the FFT, and the k-space traversal is indicated by the order of the rows of F. If one can make a reasonable estimate of the decay function for the object being scanned, then a can be computed for many k-space traversals, and the one that yields the smallest a can be used.

Methods of performing noise cancellation have been described. The method includes acquiring a plurality of data sets, analyzing a subset of the plurality of data sets, generating a filter from the subset of data sets, and applying the filter to a data set which results in a noise cancelled data set. Alternately a regressor is generated instead of a filter, the regressor is applied to the model, and then the data set is applied to the model. Another method for providing noise cancellation in measurements includes applying an excitation pulse to a subject, acquiring a baseline data set before a data set induced by the excitation pulse has developed, waiting a predetermined period of time, acquiring a data set induced by the excitation pulse, generating a regressor from the baseline data set, applying the regressor to a model, and applying a dataset to the model to obtain a noise cancelled data set. In an alternate embodiment, instead of generating a regressor both datasets are applied the model. Having described preferred embodiments of the invention it will now become apparent to those of ordinary skill in the art that other embodiments incorporating these concepts may be used. Additionally, the software included as part of the invention may be embodied in a computer program product that includes a computer useable medium. For example, such a computer usable medium can include a readable memory device, such as a hard drive device, a CD-ROM, a DVD-ROM, or a computer diskette, having computer readable program code segments stored thereon. The computer readable medium can also include a communications link, either optical, wired, or wireless, having program code segments carried thereon as digital or analog signals. Accordingly, it is submitted that that the invention should not be limited to the described embodiments but rather should be limited only by the spirit and scope of the appended claims.

Claims

CLAIMSWhat is claimed is:

1. A method for providing noise cancellation in measurements comprising: acquiring a plurality of data sets; analyzing a subset of said plurality of data sets; generating a filter from said subset of data sets; and applying said filter to a data set not included in said subset of data sets, said applying resulting in a noise cancelled data set.

2. The method of claim 1 wherein said plurality of data sets includes data from one or more performances of a task.

3. The method of claim 1 wherein said plurality of data sets includes at least one of a rest data set.

4. The method of claim 1 wherein said analyzing comprises: determining an average of said subset of data sets; subtracting said average from each data set of said subset of data sets to obtain a plurality of residual data sets.

5. The method of claim 4 wherein said generating a filter comprises generating a filter from said plurality of residual data sets. ^'

6. The method of claim 1 wherein said analyzing further comprises generating a compensation mask, said compensation mask compensating for spatial correlations in the signal and noise.

7. The method of claim 6 wherein said generating a filter from said subset of data sets comprises generating a filter from said subset of data sets and said compensation mask.

8. The method of claim 1 wherein said acquiring a plurality of data sets comprises acquiring said plurality of data sets with an fMRI system.

9. A method for providing noise cancellation in measurements comprising: acquiring a plurality of data sets; analyzing said plurality of data sets; generating a regressor from said plurality of data sets; applying said regressor to a model; and applying a dataset to said model, said applying a dataset resulting in a noise cancelled data set.

10. The method of claim 9 wherein said plurality of data sets includes data from one or more performances of a task.

11. The method of claim 9 wherein said plurality of data sets includes at least one of a rest data set.

12. The method of claim 1 wherein said analyzing comprises: determining an average of said plurality of data sets; subtracting said average from each data set of said plurality of data sets to obtain a plurality of residual data sets.

13. The method of claim 12 wherein said generating a regressor comprises generating a regressor from said plurality of residual data sets.

14. The method of claim 9 wherein said acquiring a plurality of data sets comprises acquiring said plurality of data sets with an fMRI system.

15. A method for providing noise cancellation in measurements comprising: applying an excitation pulse to a subject; acquiring a baseline data set before a data set induced by said applying an excitation pulse has developed; waiting a predetermined period of time; acquiring a data set induced by said applying an excitation pulse; generating a regressor from said baseline data set; applying said regressor to a model; and applying a dataset to said model, said applying a dataset resulting in a noise cancelled data set.

16. The method of claim 15 wherein said data set comprises data from one or more performances of a task.

17. The method of claim 15 wherein said data set comprises a rest data set.

18. The method of claim 15 wherein said acquiring a data set comprises acquiring said of data set with an fMRI system.

19. A method for providing noise cancellation in measurements comprising: applying an excitation pulse to a subject; acquiring a baseline data set before a data set induced by said applying an excitation pulse has developed; waiting a predetermined period of time; acquiring a data set induced by said applying an excitation pulse applying said baseline dataset and said dataset to a model, said applying said baseline data set and said dataset resulting in a noise cancelled data set.

20. The method of claim 19 wherein said data set comprises data from one or more performances of a task.

21. The method of claim 19 wherein said data set comprises a rest data set.

22. The method of claim 19 wherein said acquiring a data set comprises acquiring said data set with an fMRI system.

23. The method of claim 19 wherein said model comprises a General Linear Model (GLM).

24. A system for performing noise cancellation of fMRI data comprising: an MRI scanner; an MRI control and storage system in communication with said MRI scanner; an image processing system in communication with said MRI control and storage system; a user interface in communication with said MRI control and storage system; and a noise cancellation system in communication with said image processing system and said MRI control and storage system.

25. The system of claim 24 wherein said noise cancellation system includes a filter processor.

26. The system of claim 24 wherein said noise cancellation system includes a regression processor.