CN105513058A - Brain active region detection method and device - Google Patents

Abstract

The invention provides a brain active region detection method, and the method comprises the following steps: obtaining fMRI data; constructing a two-dimensional neighborhood characteristic space S of the fMRI data; carrying out the clustering search of the characteristic space S through employing a mean value drift algorithm; and obtaining an active region detection result. The invention also provides a corresponding brain active region detection device. The method and device provided by the invention can detect active voxels with stronger activity, also can sensitively detect a brain active region composed of voxels which are weak in time domain performance and stronger in frequency domain performance, and are good in noise prevention capability and high in sensitivity.

Description

Brain activation region detection method and device

Technical Field

The invention relates to the field of medical signal processing, in particular to a brain activation region detection method and a brain activation region detection device.

Background

Functional magnetic resonance imaging (fMRI) is a technique for obtaining a neuro-activation distribution based on a hemodynamic (BOLD) mechanism by measuring a magnetic change of hemoglobin caused by a change in oxygen concentration in blood. However, since the BOLD-fMRI signal varies very weakly in amplitude almost in accordance with the noise fluctuation, so that it is particularly difficult to separate the signal from the noise, the BOLD-fMRI signal extraction (also called active region detection) method can be generally divided into two categories, namely model-based and data-based.

The model-based algorithm comprises cross-correlation analysis, a general linear model and the like; the basic principle of the cross-correlation analysis is that a reference waveform simulating a cerebral blood flow dynamics response function is predefined according to prior knowledge, then correlation coefficients of a time process of a voxel and the reference waveform are sequentially calculated, whether the voxel is activated or not is determined by thresholding the correlation coefficients, namely a threshold value is determined, if the correlation coefficients of the voxel are larger than the threshold value, the voxel is judged to be in an activated state, and if the correlation coefficients are smaller than the threshold value, the voxel is judged to be in a resting state. In a general linear model, correlation between variables is calculated, after a linear relation between two variables is detected, the model is established, parameters in the model are estimated, then the model is checked by using a t-test method or an f-test method, and a brain activation image corresponding to a threshold value is obtained according to the set threshold value, so that whether the brain area is in an activation state or not is judged.

The data-based algorithm comprises a k-means clustering algorithm, independent component analysis and the like; the basic principle of the k-means clustering algorithm (k-means) is that the number of classes of data is designated as k, k sample points are randomly selected from an fMRI data set to serve as initial clustering centers of each class, then the distance from each sample to the k clustering centers is calculated, the samples are classified into the class where the clustering center closest to the sample is located, and the clustering result is expressed by the value of k. And then, an iterative updating method is adopted, each iterative process is carried out along the direction of reducing the target function based on the given clustering target function until the target function obtains the minimum value, the algorithm is converged, and the detection of the brain activation region is completed. Independent Component Analysis (ICA) is a blind signal separation method, which aims to decompose observed data to extract independent components and find information components implicit in the data. The general approach to processing fMRI data with ICA is: the brain activation region caused by stimulation is detected by performing two experiments under the same condition in the same stimulation mode to obtain two signals of each voxel as a mixed signal, separating signal components related to an event by using ICA, calculating the Z fraction of each voxel, and regarding the voxel with the value larger than a given threshold value as an activation voxel.

Algorithms based on model class, such as cross-correlation analysis, general linear model and the like, are generally based on a priori assumed model, and the quality of a detection result is directly related to the satisfaction degree of data to the model. Secondly, the above methods all belong to a univariate statistical method, whether a certain voxel in the fMRI data is activated or not is determined through analysis of the voxel, and the correlation between adjacent voxels in the fMRI data space, namely the neighborhood information of the voxel, is not considered, so that the method has a certain limitation in the detection of an activation region, and particularly under the condition of low signal-to-noise ratio, the detection sensitivity of the univariate statistical method to the activation region is low.

And (4) data class-based algorithms, such as a k-means clustering method. First, the method requires pre-specifying the value of the number of clusters K. However, the number of data set classes cannot be determined in advance, and therefore, the optimal value of the k value is difficult to accurately select. Secondly, an initial clustering center needs to be determined in the method to perform initial partitioning on the data. The detection result of the clustering is sensitive to the selection of the initial clustering center, and once the initial value is not well selected, the final convergence effect of the clustering can be influenced, and the reliability and the accuracy of the detection result are reduced. The ICA has the defects in the application range, and whether all functional nuclear magnetic resonance imaging data can be processed by the ICA method or not can be realized. The ICA method, while consistent with the results of commonly used activation region detection algorithms, is difficult to process with ICA for the detection of complex brain high-level activities.

Disclosure of Invention

The invention aims to solve the problem of low signal-to-noise ratio of the existing brain activation region detection method.

The purpose of the invention is realized by adopting the following technical scheme.

A brain activation region detection method comprises the following steps:

step S1, acquiring fMRI data;

step S2, constructing a two-dimensional neighborhood feature space S of the fMRI data;

step S3, performing clustering search on the feature space S by adopting a mean shift algorithm;

and step S4, obtaining an activation area detection result.

In a preferred embodiment of the present invention, step S2 further includes the following steps:

in step S21, a set of fMRI volume data V ═ V is given, which includes T time points_t1, 2., T }, where V_tFor fMRI volume data at time T, the time series I (p) of voxels p { I (p, T) | T ═ 1,2,.., T, p ∈ V }, where I (p, T) is the voxel p at V_tSignal strength of (1); in the fMRI imaging process, the external stimulus is a stimulus signal at a specific time point, and is set as a stimulus function e (T) (1, …, T), so that the correlation coefficient r between the time series i (p) of the voxel p and the stimulus function e (T)₁(p) is represented by:

$r_1\left(p\right)=\frac{\underset{t=1}{\overset{T}{\Σ}}(I(p,t)-\stackrel{\‾}{I}\left(p\right))\·(E\left(t\right)-\stackrel{\‾}{E})}{{\[\underset{t=1}{\overset{T}{\Σ}}{\left(I\right(p,t)-\stackrel{\‾}{I}(p\left)\right)}^2\underset{t=1}{\overset{T}{\Σ}}{\left(E\right(t)-\stackrel{\‾}{E})}^2\]}^{1/2}}---\left(1\right)$

wherein,is the mean value of the signal intensities of the voxels p,is the mean value of the stimulation function;

step S22, including neighborhood information around voxel p, defining neighborhood of voxel p as n (p), n (p) ═ pⁱ1,2, …, n is the number of voxels in the neighborhood N (p), and the average value of the correlation coefficient of the time sequence and the stimulation function of all the voxels in each voxel neighborhood N (p) is defined asExpressed as:

$\stackrel{\‾}{r_1}\left(p\right)=\frac{1}{n}\underset{i=1}{\overset{n}{\Σ}}{r}_1\left(p^i\right)---\left(2\right)$

constructing first dimension features of fMRI data

In step S23, a correlation coefficient r between a voxel p and a voxel in the voxel neighborhood is obtained₂(p), expressed as:

$r_2\left(p\right)=\frac{1}{n}\underset{i=1}{\overset{n}{\Σ}}\frac{\underset{t=1}{\overset{T}{\Σ}}(I(p,t)-\stackrel{\‾}{I}(p\left)\right)\·(I(p^i,t)-\stackrel{\‾}{I}\left(p^i\right))}{{\[\underset{t=1}{\overset{T}{\Σ}}{\left(I\right(p,t)-\stackrel{\‾}{I}(p\left)\right)}^2\underset{t=1}{\overset{T}{\Σ}}{\left(I\right(p^i,t)-\stackrel{\‾}{I}(p^i\left)\right)}^2\]}^{1/2}}---\left(3\right)$

wherein,is the mean value of the sequence and is,constructing a second dimension characteristic R of the fMRI data characteristic space for the mean value of the corresponding voxels in the neighborhood₂＝{r₂(p)|p∈V}；

Step S24, constructing the feature space S, where S ═ { R ═ R₁,R₂}。

In a preferred embodiment of the present invention, in the feature space S, the kernel function of the sampling points in the given space is k (x) and the allowable error, and the step S3 further includes the following steps:

step S31, arbitrarily selecting an initial search area circle center O in the feature space S, wherein the radius is the kernel width h;

step S32, calculating the mean shift M in the search area_h(x)：

$M_h\left(x\right)=\frac{\underset{i=1}{\overset{n}{\Σ}}{x}_{i}g\left(\right|\left|\frac{x-x_i}{h}\right|{|}^2)}{\underset{i=1}{\overset{n}{\Σ}}g\left(\right|\left|\frac{x-x_i}{h}\right|{|}^2)}-x---\left(2\right)$

The vector always points in the direction of increasing density;

step S33, if the modulus of the mean shift vector is less than the tolerance error, | M_h(x)||<The iterative algorithm ends(ii) a Otherwise, calculating x by using the following formula (3) to obtain a new circle center O', and returning to execute the step S32;

$x=\frac{\underset{i=1}{\overset{n}{\&Sigma;}}{x}_{i}g\left(\right|\left|\frac{x-x_i}{h}\right|{|}^2)}{\underset{i=1}{\overset{n}{\&Sigma;}}g\left(\right|\left|\frac{x-x_i}{h}\right|{|}^2)}---\left(3\right)$

in a preferred embodiment of the present invention, step S4 further includes the following steps: repeating steps S1-S3 until M_h(x)||<And when the formula (2) converges to a local density maximum point and x drifts to a local maximum point, calculating sample points converged to the same point in the characteristic space S as a class to obtain the detection result of the activation region.

In a preferred embodiment of the present invention, step S1 further includes the following steps: and (3) formulating an experimental stimulation scheme, collecting fMRI experimental data and preprocessing the experimental data.

A brain activation region detection device, comprising:

the data grabbing unit is used for acquiring fMRI data;

the construction unit is used for constructing a two-dimensional neighborhood characteristic space S of the fMRI data;

the searching unit is used for carrying out clustering search on the characteristic space S by adopting a mean shift algorithm; and

and the output unit is used for obtaining the detection result of the activation area.

In a preferred embodiment of the present invention, the building unit further includes the following sub-units:

a data selection subunit for giving a set of fMRI volume data V ═ V comprising T time points_t1, 2., T }, where V_tFor fMRI volume data at time T, the time series I (p) of voxels p { I (p, T) | T ═ 1,2,.., T, p ∈ V }, where I (p, T) is the voxel p at V_tSignal strength of (1); in the fMRI imaging process, the external stimulus is a stimulus signal at a specific time point, and is set as a stimulus function e (T) (1, …, T), so that the correlation coefficient r between the time series i (p) of the voxel p and the stimulus function e (T)₁(p) is represented by:

$r_1\left(p\right)=\frac{\underset{t=1}{\overset{T}{\&Sigma;}}\left(I\right(p,t)-\stackrel{\&OverBar;}{I}\left(p\right))\&CenterDot;\left(E\right(t)-\stackrel{\&OverBar;}{E})}{{\&lsqb;\underset{t=1}{\overset{T}{\&Sigma;}}{\left(I\right(p,t)-\stackrel{\&OverBar;}{I}(p\left)\right)}^2\underset{t=1}{\overset{T}{\&Sigma;}}{\left(E\right(t)-\stackrel{\&OverBar;}{E})}^2\&rsqb;}^{1/2}}---\left(1\right)$

wherein,is the mean value of the signal intensities of the voxels p,is the mean value of the stimulation function;

a first dimension feature construction subunit, configured to introduce neighborhood information around a voxel p, and define a neighborhood of the voxel p as n (p), where n (p) ═ pⁱ1,2, …, n is the number of voxels in the neighborhood N (p), and the average value of the correlation coefficient of the time sequence and the stimulation function of all the voxels in each voxel neighborhood N (p) is defined asExpressed as:

$\stackrel{\&OverBar;}{r_1}\left(p\right)=\frac{1}{n}\underset{i=1}{\overset{n}{\&Sigma;}}{r}_1\left(p^i\right)---\left(2\right)$

constructing first dimension features of fMRI data

A second dimension feature construction subunit for obtaining the correlation coefficient r of the voxel p and the voxels in the voxel neighborhood₂(p), expressed as:

$r_2\left(p\right)=\frac{1}{n}\underset{i=1}{\overset{n}{\&Sigma;}}\frac{\underset{t=1}{\overset{T}{\&Sigma;}}\left(I\right(p,t)-\stackrel{\&OverBar;}{I}\left(p\right))\&CenterDot;\left(I\right(p^i,t)-\stackrel{\&OverBar;}{I}\left(p^i\right))}{{\&lsqb;\underset{t=1}{\overset{T}{\&Sigma;}}{\left(I\right(p,t)-\stackrel{\&OverBar;}{I}(p\left)\right)}^2\underset{t=1}{\overset{T}{\&Sigma;}}{\left(I\right(p^i,t)-\stackrel{\&OverBar;}{I}(p^i\left)\right)}^2\&rsqb;}^{1/2}}---\left(3\right)$

constructing a second dimension feature R of the fMRI data feature space₂＝{r₂(p) | p ∈ V }, wherein,is the mean value of the sequence and is,is the mean value of the corresponding voxels in the neighborhood; and

a feature space construction subunit configured to construct the feature space S, where S ═ { R ═ R₁,R₂}。

In a preferred embodiment of the present invention, in the feature space S, the kernel function of the sampling points in the given space is k (x) and the allowable error, and the searching unit further includes the following sub-units:

the condition setting subunit is used for randomly selecting the circle center O of the initial search area in the feature space S, and the radius is the kernel width h;

a calculation subunit for calculating a mean shift M within the search area_h(x)：

$M_h\left(x\right)=\frac{\underset{i=1}{\overset{n}{\&Sigma;}}{x}_{i}g\left(\right|\left|\frac{x-x_i}{h}\right|{|}^2)}{\underset{i=1}{\overset{n}{\&Sigma;}}g\left(\right|\left|\frac{x-x_i}{h}\right|{|}^2)}-x---\left(2\right)$

The vector always points in the direction of increasing density;

a judging subunit, configured to judge: if the modulus of the mean shift vector is less than the tolerance error, | M_h(x)||<The iterative algorithm ends; otherwise, calculating x by using the following formula (3) to obtain a new circle center O ', and returning the new circle center O' to the calculating subunit for continuous calculation;

$x=\frac{\underset{i=1}{\overset{n}{\&Sigma;}}{x}_{i}g\left(\right|\left|\frac{x-x_i}{h}\right|{|}^2)}{\underset{i=1}{\overset{n}{\&Sigma;}}g\left(\right|\left|\frac{x-x_i}{h}\right|{|}^2)}---\left(3\right)$

in a preferred embodiment of the present invention, the output unit is further configured to calculate a sample point converging to the same point in the feature space S as a class, so as to obtain the detection result of the active region.

In a preferred embodiment of the present invention, the data capturing unit is further configured to pre-process the fMRI experimental data.

Compared with the prior art, the method and the device for detecting the brain activation region provided by the invention can detect the activation voxels with stronger activity, can sensitively detect the brain activation region formed by the voxels with weak activity but highly related spatial domain, and have good anti-noise capability and high sensitivity.

The foregoing description is only an overview of the technical solutions of the present invention, and in order to make the technical means of the present invention more clearly understood, the present invention may be implemented in accordance with the content of the description, and in order to make the above and other objects, features, and advantages of the present invention more clearly understood, the following preferred embodiments are described in detail with reference to the accompanying drawings.

Drawings

Fig. 1 is a schematic flow chart of a method for detecting a brain activation region according to a first embodiment of the present invention.

Fig. 2 is a schematic structural diagram of a brain activation region detection apparatus according to a second embodiment of the present invention.

Fig. 3 is a schematic structural diagram of a construction unit of a brain activation region detection apparatus according to a second embodiment of the present invention.

Fig. 4 is a schematic structural diagram of a search unit of a brain activation region detection apparatus according to a second embodiment of the present invention.

Detailed Description

To facilitate an understanding of the invention, the invention will now be described more fully with reference to the accompanying drawings. Preferred embodiments of the present invention are shown in the drawings. This invention may, however, be embodied in many different forms and should not be construed as limited to the embodiments set forth herein. Rather, these embodiments are provided so that this disclosure will be thorough and complete.

Unless defined otherwise, all technical and scientific terms used herein have the same meaning as commonly understood by one of ordinary skill in the art to which this invention belongs. The terminology used herein in the description of the invention is for the purpose of describing particular embodiments only and is not intended to be limiting of the invention. As used herein, the term "and/or" includes any and all combinations of one or more of the associated listed items.

Referring to fig. 1, fig. 1 is a schematic flow chart of a method for detecting a brain activation region according to a first embodiment of the present invention, including the following steps:

in step S1, fMRI data is acquired.

Specifically, an experimental stimulation scheme is firstly formulated, and fMRI experimental data are collected and preprocessed. The purpose of the preprocessing is to detect and repair artifacts produced by the magnetic resonance scanner or the subject person at the time of data acquisition or to prepare for subsequent data processing.

Step S2, a two-dimensional neighborhood feature space S of the fMRI data is constructed.

Specifically, the method comprises the following steps:

in step S21, a set of fMRI volume data V ═ V is given, which includes T time points_t1, 2., T }, where V_tFor fMRI volume data at time T, the time series I (p) of voxels p { I (p, T) | T ═ 1,2,.., T, p ∈ V }, where I (p, T) is the voxel p at V_tSignal strength of (1); in the fMRI imaging process, the external stimulus is a stimulus signal at a specific time point, and is set as a stimulus function e (T) (1, …, T), so that the correlation coefficient r between the time series i (p) of the voxel p and the stimulus function e (T)₁(p) can be obtained from cross-correlation analysis and is expressed as:

$r_1\left(p\right)=\frac{\underset{t=1}{\overset{T}{\&Sigma;}}\left(I\right(p,t)-\stackrel{\&OverBar;}{I}\left(p\right))\&CenterDot;\left(E\right(t)-\stackrel{\&OverBar;}{E})}{{\&lsqb;\underset{t=1}{\overset{T}{\&Sigma;}}{\left(I\right(p,t)-\stackrel{\&OverBar;}{I}(p\left)\right)}^2\underset{t=1}{\overset{T}{\&Sigma;}}{\left(E\right(t)-\stackrel{\&OverBar;}{E})}^2\&rsqb;}^{1/2}}---\left(1\right)$

wherein,is the mean value of the signal intensities of the voxels p,mean value of stimulation function. Coefficient of correlation r₁The size of (p) reflects the degree of relevance of the voxel to participate in task stimulation.

Step S22, including neighborhood information around voxel p, defining neighborhood of voxel p as n (p), n (p) ═ pⁱ1,2, …, n is the number of voxels in the neighborhood N (p), and the average value of the correlation coefficient of the time sequence and the stimulation function of all the voxels in each voxel neighborhood N (p) is defined asExpressed as:

$\stackrel{\&OverBar;}{r_1}\left(p\right)=\frac{1}{n}\underset{i=1}{\overset{n}{\&Sigma;}}{r}_1\left(p^i\right)---\left(2\right)$

constructing first dimension features of fMRI data

In step S23, a correlation coefficient r between a voxel p and a voxel in the voxel neighborhood is obtained₂(p), correlation coefficient r₂(p) reflects the consistency of the fMRI data neighborhood space, expressed as:

$r_2\left(p\right)=\frac{1}{n}\underset{i=1}{\overset{n}{\&Sigma;}}\frac{\underset{t=1}{\overset{T}{\&Sigma;}}\left(I\right(p,t)-\stackrel{\&OverBar;}{I}\left(p\right))\&CenterDot;\left(I\right(p^i,t)-\stackrel{\&OverBar;}{I}\left(p^i\right))}{{\&lsqb;\underset{t=1}{\overset{T}{\&Sigma;}}{\left(I\right(p,t)-\stackrel{\&OverBar;}{I}(p\left)\right)}^2\underset{t=1}{\overset{T}{\&Sigma;}}{\left(I\right(p^i,t)-\stackrel{\&OverBar;}{I}(p^i\left)\right)}^2\&rsqb;}^{1/2}}---\left(3\right)$

wherein,is the mean value of the sequence and is,constructing a second dimension characteristic R of the fMRI data characteristic space for the mean value of the corresponding voxels in the neighborhood₂＝{r₂(p)|p∈V}；

Step S24, constructing the feature space S, where S ═ { R ═ R₁,R₂}。

In other embodiments, the value range of n may be set according to practical situations, and may be 8 neighborhoods, 4 neighborhoods, or the like.

And step S3, performing clustering search on the feature space S by adopting a mean shift algorithm.

Specifically, in the feature space S, the kernel function of the sampling points in the given space is k (x) and the allowable error, and further includes the steps of:

step S31, arbitrarily selecting an initial search area circle center O in the feature space S, wherein the radius is the kernel width h;

step S32, calculating the mean shift M in the search area_h(x)：

$M_h\left(x\right)=\frac{\underset{i=1}{\overset{n}{\&Sigma;}}{x}_{i}g\left(\right|\left|\frac{x-x_i}{h}\right|{|}^2)}{\underset{i=1}{\overset{n}{\&Sigma;}}g\left(\right|\left|\frac{x-x_i}{h}\right|{|}^2)}-x---\left(2\right)$

The vector always points in the direction of increasing density;

step S33, if the modulus of the mean shift vector is less than the tolerance error, | M_h(x)||<The iterative algorithm ends; otherwise, calculating x by using the following formula (3) to obtain a new circle center O', and returning to execute the step S32;

$x=\frac{\underset{i=1}{\overset{n}{\&Sigma;}}{x}_{i}g\left(\right|\left|\frac{x-x_i}{h}\right|{|}^2)}{\underset{i=1}{\overset{n}{\&Sigma;}}g\left(\right|\left|\frac{x-x_i}{h}\right|{|}^2)}---\left(3\right)$

in step S31, the key parameter is kernel width h, which determines the size of the search area, and the difference in the value of h affects the result of the algorithm.

And step S4, obtaining an activation area detection result.

Specifically, the method comprises the following steps: the steps S1 to S3 (in the present embodiment, the steps S1, S21 to S23, S31 to S33) are repeated until | | M_h(x)||<When the formula (2) converges to the local density maximum point and x drifts to the local maximum point, the formula (2) is used for converting the x into the local density maximum pointSample points converging to the same point in the feature space S are calculated as a class, and an activation area detection result of a VN-MSC (VN-MSC) algorithm is obtained.

The method for detecting the brain activation region provided by the embodiment utilizes the characteristic that voxels in the same region of the brain have similar properties, namely, the brain functional regions are distributed in a block shape and are clusters formed by connected voxels in fMRI data space, but are not distributed in a dot shape formed by single voxels, so that the neighborhood correlation characteristics of fMRI data voxel space are fully utilized to construct a feature space, and the mean shift clustering algorithm is combined to extract the brain activation regions under different stimulation conditions, so that not only the activated voxels with strong activity can be detected, but also the brain activation regions formed by the voxels with weak activity and high spatial domain correlation can be sensitively detected, and the method has good anti-noise capability and high sensitivity.

Referring to fig. 2, fig. 2 is a schematic structural diagram of a brain activation region detection apparatus 100 according to a second embodiment of the present invention.

The brain activation region detection apparatus 100 includes: a data capture unit 10 for acquiring fMRI data; the constructing unit 20 is configured to construct a two-dimensional neighborhood feature space S of fMRI data; the searching unit 30 is used for carrying out clustering search on the characteristic space S by adopting a mean shift algorithm; and an output unit 40 for obtaining an active area detection result.

Referring to fig. 3, fig. 3 is a schematic structural diagram of the building unit 20. The building unit 20 further comprises the following sub-units:

a data selection subunit 201 for giving a set of fMRI volume data V ═ V for T time points_t1, 2., T }, where V_tFor fMRI volume data at time T, the time series I (p) of voxels p { I (p, T) | T ═ 1,2,.., T, p ∈ V }, where I (p, T) is the voxel p at V_tSignal strength of (1); in the fMRI imaging process, the external stimulus is a stimulus signal at a specific time point, and is set as a stimulus function e (T) (1, …, T), so that the correlation coefficient r between the time series i (p) of the voxel p and the stimulus function e (T)₁(p) is represented by:

$r_1\left(p\right)=\frac{\underset{t=1}{\overset{T}{\&Sigma;}}\left(I\right(p,t)-\stackrel{\&OverBar;}{I}\left(p\right))\&CenterDot;\left(E\right(t)-\stackrel{\&OverBar;}{E})}{{\&lsqb;\underset{t=1}{\overset{T}{\&Sigma;}}{\left(I\right(p,t)-\stackrel{\&OverBar;}{I}(p\left)\right)}^2\underset{t=1}{\overset{T}{\&Sigma;}}{\left(E\right(t)-\stackrel{\&OverBar;}{E})}^2\&rsqb;}^{1/2}}---\left(1\right)$

wherein,is the mean value of the signal intensities of the voxels p,is the mean value of the stimulation function;

a first dimension feature constructing subunit 202, configured to introduce neighborhood information around the voxel p, and define a neighborhood of the voxel p as n (p), where n (p) { p { (p) }ⁱ1,2, …, n is the number of voxels in the neighborhood N (p), and the average value of the correlation coefficient of the time sequence and the stimulation function of all the voxels in each voxel neighborhood N (p) is defined asExpressed as:

$\stackrel{\&OverBar;}{r_1}\left(p\right)=\frac{1}{n}\underset{i=1}{\overset{n}{\&Sigma;}}{r}_1\left(p^i\right)---\left(2\right)$

constructing first dimension features of fMRI data

A second dimension feature construction subunit 203, configured to obtain a correlation coefficient r between the voxel p and voxels in the voxel neighborhood₂(p), expressed as:

$r_2\left(p\right)=\frac{1}{n}\underset{i=1}{\overset{n}{\&Sigma;}}\frac{\underset{t=1}{\overset{T}{\&Sigma;}}\left(I\right(p,t)-\stackrel{\&OverBar;}{I}\left(p\right))\&CenterDot;\left(I\right(p^i,t)-\stackrel{\&OverBar;}{I}\left(p^i\right))}{{\&lsqb;\underset{t=1}{\overset{T}{\&Sigma;}}{\left(I\right(p,t)-\stackrel{\&OverBar;}{I}(p\left)\right)}^2\underset{t=1}{\overset{T}{\&Sigma;}}{\left(I\right(p^i,t)-\stackrel{\&OverBar;}{I}(p^i\left)\right)}^2\&rsqb;}^{1/2}}---\left(3\right)$

constructing a second dimension feature R of the fMRI data feature space₂＝{r₂(p) | p ∈ V }, wherein,is the mean value of the sequence and is,is the mean value of the corresponding voxels in the neighborhood; and

a feature space construction subunit 204, configured to construct the feature space S, where S ═ { R ═ R₁,R₂}。

In the feature space S, given that the kernel function of the sampling points in the space is k (x) and the tolerance error, referring to fig. 4, the searching unit 30 further includes the following sub-units:

a condition setting subunit 301, configured to arbitrarily select an initial search area in the feature space S, where a circle center is O and a radius is a kernel width h;

a calculation subunit 302 for calculating a mean shift M in the search area_h(x)：

$M_h\left(x\right)=\frac{\underset{i=1}{\overset{n}{\&Sigma;}}{x}_{i}g\left(\right|\left|\frac{x-x_i}{h}\right|{|}^2)}{\underset{i=1}{\overset{n}{\&Sigma;}}g\left(\right|\left|\frac{x-x_i}{h}\right|{|}^2)}-x---\left(2\right)$

The vector always points in the direction of increasing density;

a judging subunit 303, configured to judge: if the modulus of the mean shift vector is less than the tolerance error, | M_h(x)||<The iterative algorithm ends; otherwise, calculating x by using the following formula (3) to obtain a new circle center O ', and returning the new circle center O' to the calculating subunit 302 for continuous calculation;

$x=\frac{\underset{i=1}{\overset{n}{\&Sigma;}}{x}_{i}g\left(\right|\left|\frac{x-x_i}{h}\right|{|}^2)}{\underset{i=1}{\overset{n}{\&Sigma;}}g\left(\right|\left|\frac{x-x_i}{h}\right|{|}^2)}---\left(3\right)$

the output unit 40 is further configured to calculate a sample point converged to the same point in the feature space S as a class when the formula (2) converges to the local density maximum point, so as to obtain an active region detection result.

The data capture unit 10 is also used for preprocessing the fMRI experimental data.

The brain activation region detection device provided by the embodiment utilizes the characteristic that voxels in the same region of the brain have similar properties, namely, the brain functional regions are distributed in a block shape and are clusters formed by connected voxels in fMRI data space, but are not distributed in a dot shape formed by single voxels, so that the neighborhood correlation characteristics of fMRI data voxel space are fully utilized to construct a feature space, and the brain activation regions under different stimulation conditions are extracted by combining a mean shift clustering algorithm, so that not only can the activated voxels with strong activity be detected, but also the brain activation regions formed by the voxels with weak activity but highly correlated spatial domains can be sensitively detected, and the brain activation region detection device has good anti-noise capability and high sensitivity.

The above-mentioned embodiments only express several embodiments of the present invention, and the description thereof is more specific and detailed, but not construed as limiting the scope of the present invention. It should be noted that, for a person skilled in the art, several variations and modifications can be made without departing from the inventive concept, which falls within the scope of the present invention. Therefore, the protection scope of the present patent shall be subject to the appended claims.

Claims (10)

1. A brain activation region detection method is characterized by comprising the following steps:

step S1, acquiring fMRI data;

step S2, constructing a two-dimensional neighborhood feature space S of the fMRI data;

step S3, performing clustering search on the feature space S by adopting a mean shift algorithm;

and step S4, obtaining an activation area detection result.

2. The detection method according to claim 1, wherein the step S2 includes:

in step S21, a set of fMRI volume data V ═ V is given, which includes T time points_t1, 2., T }, where V_tFor fMRI volume data at time T, the time series I (p) of voxels p { I (p, T) | T ═ 1,2,.., T, p ∈ V }, where I (p, T) is the voxel p at V_tSignal strength of (1); in the fMRI imaging process, the external stimulus is a stimulus signal at a specific time point, and is set as a stimulus function e (T) (1.,. T.), and the time series i (p) of the voxels p and the correlation coefficient r of the stimulus function e (T) are set as₁(p) is represented by:

$r_1\left(p\right)=\frac{\underset{t=1}{\overset{T}{\&Sigma;}}\left(I\right(p,t)-\stackrel{\&OverBar;}{I}(p\left)\right)\&CenterDot;\left(E\right(t)-\stackrel{\&OverBar;}{E})}{{\&lsqb;\underset{t=1}{\overset{T}{\&Sigma;}}{\left(I\right(p,t)-\stackrel{\&OverBar;}{I}(p\left)\right)}^2\underset{t=1}{\overset{T}{\&Sigma;}}{\left(E\right(t)-\stackrel{\&OverBar;}{E})}^2\&rsqb;}^{1/2}}---\left(1\right)$

wherein,is the mean value of the signal intensities of the voxels p,is the mean value of the stimulation function;

step S22, including neighborhood information around voxel p, defining neighborhood of voxel p as n (p), n (p) ═ pⁱ1, 2.. n }, where n is the number of voxels in the neighborhood n (p), and each voxel is takenThe average value of the correlation coefficient of the time series of all the voxels in the neighborhood N (p) and the stimulus function is defined asExpressed as:

$\stackrel{\&OverBar;}{r_1}\left(p\right)=\frac{1}{n}\underset{i=1}{\overset{n}{\&Sigma;}}{r}_1\left(p^i\right)---\left(2\right)$

constructing first dimension features of fMRI data

In step S23, a correlation coefficient r between a voxel p and a voxel in the voxel neighborhood is obtained₂(p), expressed as:

$r_2\left(p\right)=\frac{1}{n}\underset{i=1}{\overset{n}{\&Sigma;}}\frac{\underset{t=1}{\overset{T}{\&Sigma;}}\left(I\right(p,t)-\stackrel{\&OverBar;}{I}(p\left)\right)\&CenterDot;\left(I\right(p^i,t)-\stackrel{\&OverBar;}{I}(p^i\left)\right)}{{\&lsqb;\underset{t=1}{\overset{T}{\&Sigma;}}{\left(I\right(p,t)-\stackrel{\&OverBar;}{I}(p\left)\right)}^2\underset{t=1}{\overset{T}{\&Sigma;}}{\left(I\right(p^i,t)-\stackrel{\&OverBar;}{I}(p^i\left)\right)}^2\&rsqb;}^{1/2}}---\left(3\right)$

wherein,is the mean value of the sequence and is,constructing a second dimension characteristic R of the fMRI data characteristic space for the mean value of the corresponding voxels in the neighborhood₂＝{r₂(p)|p∈V}；

Step S24, constructing the feature space S, where S ═ { R ═ R₁,R₂}。

3. The detection method according to claim 1, wherein in the feature space S, the kernel function of the sampling points in a given space is k (x) and the allowable error, and the step S3 includes:

step S31, selecting an initial search area in the feature space S arbitrarily, wherein the center of the circle O and the radius are the kernel width h;

step S32, calculating the mean shift M in the search area_h(x)：

$M_h\left(x\right)=\frac{\underset{i=1}{\overset{n}{\&Sigma;}}{x}_{i}g\left(\right|\left|\frac{x-x_i}{h}\right|{|}^2)}{\underset{i=1}{\overset{n}{\&Sigma;}}g\left(\right|\left|\frac{x-x_i}{h}\right|{|}^2)}-x---\left(2\right)$

The vector always points in the direction of increasing density;

step S33, if the modulus of the mean shift vector is less than the tolerance error, | M_h(x)||<The iterative algorithm ends; otherwise, calculating x by using the following formula (3) to obtain a new circle center O', and returning to execute the step S32;

$x=\frac{\underset{i=1}{\overset{n}{\&Sigma;}}{x}_{i}g\left(\right|\left|\frac{x-x_i}{h}\right|{|}^2)}{\underset{i=1}{\overset{n}{\&Sigma;}}g\left(\right|\left|\frac{x-x_i}{h}\right|{|}^2)}---\left(3\right)$

4. the detection method according to any one of claims 1 to 3, wherein the step S4 includes: repeating steps S1-S3 until M_h(x)||<And when the formula (2) converges to a local density maximum point and x drifts to a local maximum point, calculating sample points converged to the same point in the characteristic space S as a class to obtain the detection result of the activation region.

5. The detection method according to claim 1, wherein the step S1 includes: and (3) formulating an experimental stimulation scheme, collecting fMRI experimental data and preprocessing the experimental data.

6. A brain activation region detecting device, comprising:

the data grabbing unit is used for acquiring fMRI data;

the construction unit is used for constructing a two-dimensional neighborhood characteristic space S of the fMRI data;

the searching unit is used for carrying out clustering search on the characteristic space S by adopting a mean shift algorithm; and

and the output unit is used for obtaining the detection result of the activation area.

7. The detection apparatus of claim 6, wherein the construction unit comprises:

a data selection subunit for giving a set of fMRI volume data V ═ V comprising T time points_t1, 2., T }, where V_tFor fMRI volume data at time T, the time series I (p) of voxels p { I (p, T) | T ═ 1,2,.., T, p ∈ V }, where I (p, T) is the voxel p at V_tSignal strength of (1); in the fMRI imaging process, the external stimulus is a stimulus signal at a specific time point, and is set as a stimulus function e (T) (1.,. T.), and the time series i (p) of the voxels p and the correlation coefficient r of the stimulus function e (T) are set as₁(p) is represented by:

$r_1\left(p\right)=\frac{\underset{t=1}{\overset{T}{\&Sigma;}}\left(I\right(p,t)-\stackrel{\&OverBar;}{I}(p\left)\right)\&CenterDot;\left(E\right(t)-\stackrel{\&OverBar;}{E})}{{\&lsqb;\underset{t=1}{\overset{T}{\&Sigma;}}{\left(I\right(p,t)-\stackrel{\&OverBar;}{I}(p\left)\right)}^2\underset{t=1}{\overset{T}{\&Sigma;}}{\left(E\right(t)-\stackrel{\&OverBar;}{E})}^2\&rsqb;}^{1/2}}---\left(1\right)$

wherein,is the mean value of the signal intensities of the voxels p,is the mean value of the stimulation function;

a first dimension feature construction subunit, configured to introduce neighborhood information around a voxel p, and define a neighborhood of the voxel p as n (p), where n (p) ═ pⁱ1, 2.. n.n.n is the number of voxels in the neighborhood N (p), and the average value of the correlation coefficients of the time series of all the voxels in each voxel neighborhood N (p) and the stimulation function is defined asExpressed as:

$\stackrel{\&OverBar;}{r_1}\left(p\right)=\frac{1}{n}\underset{i=1}{\overset{n}{\&Sigma;}}{r}_1\left(p^i\right)---\left(2\right)$

constructing first dimension features of fMRI data

A second dimension feature construction subunit for obtaining the correlation coefficient r of the voxel p and the voxels in the voxel neighborhood₂(p), expressed as:

$r_2\left(p\right)=\frac{1}{n}\underset{i=1}{\overset{n}{\&Sigma;}}\frac{\underset{t=1}{\overset{T}{\&Sigma;}}\left(I\right(p,t)-\stackrel{\&OverBar;}{I}(p\left)\right)\&CenterDot;\left(I\right(p^i,t)-\stackrel{\&OverBar;}{I}(p^i\left)\right)}{{\&lsqb;\underset{t=1}{\overset{T}{\&Sigma;}}{\left(I\right(p,t)-\stackrel{\&OverBar;}{I}(p\left)\right)}^2\underset{t=1}{\overset{T}{\&Sigma;}}{\left(I\right(p^i,t)-\stackrel{\&OverBar;}{I}(p^i\left)\right)}^2\&rsqb;}^{1/2}}---\left(3\right)$

constructing a second dimension feature R of the fMRI data feature space₂＝{r₂(p) | p ∈ V }, wherein,is the mean value of the sequence and is,is the mean value of the corresponding voxels in the neighborhood; and

a feature space construction subunit configured to construct the feature space S, where S ═ { R ═ R₁,R₂}。

8. The detecting device according to claim 1, wherein in the feature space S, the kernel function of the sampling points in a given space is k (x) and the tolerance error, the searching unit includes:

the condition setting subunit is used for randomly selecting an initial search area in the feature space S, wherein the circle center is O, and the radius is the kernel width h;

a calculation subunit for calculating a mean shift M within the search area_h(x)：

$M_h\left(x\right)=\frac{\underset{i=1}{\overset{n}{\&Sigma;}}{x}_{i}g\left(\right|\left|\frac{x-x_i}{h}\right|{|}^2)}{\underset{i=1}{\overset{n}{\&Sigma;}}g\left(\right|\left|\frac{x-x_i}{h}\right|{|}^2)}-x---\left(2\right)$

The vector always points in the direction of increasing density;

a judging subunit, configured to judge: if the modulus of the mean shift vector is less than the tolerance error, | M_h(x)||<The iterative algorithm ends; otherwise, calculating x by using the following formula (3) to obtain a new circle center O ', and returning the new circle center O' to the calculating subunit for continuous calculation;

$x=\frac{\underset{i=1}{\overset{n}{\&Sigma;}}{x}_{i}g\left(\right|\left|\frac{x-x_i}{h}\right|{|}^2)}{\underset{i=1}{\overset{n}{\&Sigma;}}g\left(\right|\left|\frac{x-x_i}{h}\right|{|}^2)}---\left(3\right)$

9. the detecting device according to any one of claims 6 to 8, wherein the output unit is further configured to calculate sample points converging to the same point in the feature space S as a class to obtain the active region detecting result.

10. The testing device of claim 6, wherein the data capture unit is further configured to pre-process the fMRI experimental data.