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

Brain active region detection method and device Download PDF

Info

Publication number
CN105513058A
CN105513058A CN201510862020.4A CN201510862020A CN105513058A CN 105513058 A CN105513058 A CN 105513058A CN 201510862020 A CN201510862020 A CN 201510862020A CN 105513058 A CN105513058 A CN 105513058A
Authority
CN
China
Prior art keywords
sigma
overbar
voxel
fmri
voxels
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Granted
Application number
CN201510862020.4A
Other languages
Chinese (zh)
Other versions
CN105513058B (en
Inventor
张睿
高欣
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Suzhou Institute of Biomedical Engineering and Technology of CAS
Original Assignee
Suzhou Institute of Biomedical Engineering and Technology of CAS
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Suzhou Institute of Biomedical Engineering and Technology of CAS filed Critical Suzhou Institute of Biomedical Engineering and Technology of CAS
Priority to CN201510862020.4A priority Critical patent/CN105513058B/en
Publication of CN105513058A publication Critical patent/CN105513058A/en
Application granted granted Critical
Publication of CN105513058B publication Critical patent/CN105513058B/en
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00—Image analysis
    • G06T7/0002—Inspection of images, e.g. flaw detection
    • G06T7/0012—Biomedical image inspection
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00—Indexing scheme for image analysis or image enhancement
    • G06T2207/10—Image acquisition modality
    • G06T2207/10072—Tomographic images
    • G06T2207/10088—Magnetic resonance imaging [MRI]
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00—Indexing scheme for image analysis or image enhancement
    • G06T2207/20—Special algorithmic details
    • G06T2207/20081—Training; Learning
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00—Indexing scheme for image analysis or image enhancement
    • G06T2207/30—Subject of image; Context of image processing
    • G06T2207/30004—Biomedical image processing
    • G06T2207/30016—Brain

Landscapes

  • Engineering & Computer Science (AREA)
  • Health & Medical Sciences (AREA)
  • General Health & Medical Sciences (AREA)
  • Medical Informatics (AREA)
  • Nuclear Medicine, Radiotherapy & Molecular Imaging (AREA)
  • Radiology & Medical Imaging (AREA)
  • Quality & Reliability (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Magnetic Resonance Imaging Apparatus (AREA)

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

一种脑激活区检测方法和装置Method and device for detecting brain activation area

技术领域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 technique

功能磁共振成像(functionalMagneticResonanceImaging,fMRI)是一种通过测量血液中氧浓度变化引起的血红蛋白的磁性改变，得到基于血液动力学(BloodOxygenLevelDependence,BOLD)机制的神经激活分布的技术。然而，由于BOLD-fMRI信号变化幅度非常微弱，与噪声波动几乎一致，使得从噪声中分离信号变得尤为困难目前，BOLD-fMRI信号提取(也称激活区检测)方法一般可分为两类，分别是基于模型类和基于数据类。Functional Magnetic Resonance Imaging (fMRI) is a technology that measures the magnetic changes of hemoglobin caused by changes in oxygen concentration in blood, and obtains the distribution of neural activation based on the hemodynamic (Blood Oxygen Level Dependence, BOLD) mechanism. However, because the variation of BOLD-fMRI signal is very weak, which is almost consistent with the noise fluctuation, it is very difficult to separate the signal from the noise. At present, BOLD-fMRI signal extraction (also called activation area detection) methods can generally be divided into two categories, They are model-based and data-based, respectively.

基于模型的算法有互相关分析、一般线性模型等；其中互相关分析的基本原理是根据先验知识，预先定义一个模拟脑血流动力学响应函数的参考波形，然后依次计算体素的时间过程与该参考波形的相关系数，通过阈值化相关系数来确定体素激活与否，即确定一个阈值，若体素的相关系数大于该阈值，则判定体素处于激活状态，若相关系数小于该阈值，则体素处于静息状态。一般线性模型是通过计算变量之间的相关性，在检测到两个变量之间的线性关系后，建立模型，对其中的参数进行估计，然后用t检验或f检验方法对模型进行检验，根据所设置的阈值，得到与阈值相对应的大脑激活图像，从而判定脑区是否处于激活状态。Model-based algorithms include cross-correlation analysis, general linear model, etc. The basic principle of cross-correlation analysis is to pre-define a reference waveform that simulates the cerebral hemodynamic response function based on prior knowledge, and then calculate the time course of voxels in turn The correlation coefficient with the reference waveform is determined by thresholding the correlation coefficient to determine whether the voxel is activated or not, that is, a threshold is determined. If the correlation coefficient of the voxel is greater than the threshold, it is determined that the voxel is in an active state. If the correlation coefficient is less than the threshold , the voxel is in a resting state. The general linear model is to calculate the correlation between variables, after detecting the linear relationship between two variables, build the model, estimate the parameters, and then use the t test or f test method to test the model, according to The set threshold value is used to obtain the brain activation image corresponding to the threshold value, so as to determine whether the brain region is in an activated state.

基于数据的算法有k均值聚类算法，独立成分分析等；其中k均值聚类算法(k-means)的基本原理是首先指定数据的类别数为k，随机从fMRI数据集中选取k个样本点作为每个类的初始聚类中心，然后计算各个样本到k个聚类中心的距离，把样本归到离它最近的那个聚类中心所在的类，聚类结果由k值表达。之后采用迭代更新的方法，基于给定的聚类目标函数使每一次迭代过程都是沿着目标函数减小的方向进行，直到目标函数取得最小值，算法收敛，完成对脑激活区的检测。独立成分分析(independentcomponentanalysis,ICA)是一种盲信号分离的方法，其目的是将观察到的数据进行分解提取独立成分，发现数据中隐含的信息成分。用ICA处理fMRI数据一般做法是：用同样的刺激方式在同样的情况下作两次实验得到每个体素的两个信号作混合信号，用ICA分离出与事件相关的信号成分，然后计算每个体素的Z分数，把其值大于给定的阈值的体素认为是激活体素，由此检测由刺激引起的脑激活区域。Data-based algorithms include k-means clustering algorithm, independent component analysis, etc.; the basic principle of k-means clustering algorithm is to first specify the number of data categories as k, and randomly select k sample points from the fMRI data set As the initial clustering center of each class, then calculate the distance from each sample to k clustering centers, and classify the samples into the class of the nearest clustering center, and the clustering result is expressed by the k value. Afterwards, the iterative update method is adopted. Based on the given clustering objective function, each iterative process is carried out in the direction of decreasing the objective function until the objective function achieves the minimum value, the algorithm converges, and the detection of the brain activation area is completed. Independent component analysis (ICA) is a method of blind signal separation, the purpose of which is to decompose the observed data to extract independent components and discover the hidden information components in the data. The general method of processing fMRI data with ICA is: use the same stimulation method to conduct two experiments under the same situation to obtain two signals of each voxel as a mixed signal, use ICA to separate the signal components related to the event, and then calculate the signal components of each voxel. The z-score of a voxel, the voxel whose value is greater than a given threshold is considered to be an activated voxel, thereby detecting the brain activation area caused by the stimulus.

基于模型类的算法如互相关分析、一般线性模型等这类算法一般是基于先验假设模型的，而检测结果的好坏直接和数据对模型的满足程度相关。其次，它们都属于一元统计方法，通过对fMRI数据中某个体素的分析确定其是否被激活，并没有考虑fMRI数据空间中相邻体素间的相互关系，即体素的邻域信息，因而在激活区域检测方面具有一定的局限性，特别是在低信噪比的条件下，这类一元统计方法对激活区的检测灵敏度较低。Algorithms based on models, such as cross-correlation analysis, general linear model, etc., are generally based on a priori assumption models, and the quality of the detection results is directly related to the satisfaction degree of the data to the model. Secondly, they all belong to unary statistical methods, which determine whether a voxel is activated by analyzing the fMRI data, and do not consider the relationship between adjacent voxels in the fMRI data space, that is, the neighborhood information of the voxel. There are certain limitations in the detection of activation regions, especially under the condition of low signal-to-noise ratio, and the detection sensitivity of such unary statistical methods to activation regions is low.

基于数据类的算法，如k均值聚类方法。首先该方法需预先指定聚类数目K值。然而一般情况下无法预先判定数据集类别数目，因此k值的最优值是很难准确选择的。其次该方法中需确定一个初始聚类中心对数据进行初始划分。而聚类的检测结果对于初始聚类中心的选择较为敏感，一旦初始值选择的不好，将会影响聚类最终的收敛效果，降低检测结果的可靠性和准确性。而ICA的不足体现在其适用范围，是否所有功能核磁共振成像数据都可以采用ICA方法处理。ICA方法虽然可以与一般常用的激活区检测算法的结果保持一致性，但对复杂的大脑的高级活动的检测来说，用ICA处理较为困难。Algorithms based on data classes, such as the k-means clustering method. First of all, this method needs to pre-specify the K value of the number of clusters. However, in general, the number of categories in a data set cannot be determined in advance, so it is difficult to accurately select the optimal value of k. Secondly, in this method, an initial cluster center needs to be determined to initially divide the data. However, the detection results of clustering are more sensitive to the selection of the initial cluster center. Once the initial value is not well selected, it will affect the final convergence effect of the clustering and reduce the reliability and accuracy of the detection results. The deficiency of ICA is reflected in its scope of application, whether all fMRI data can be processed by ICA method. Although the ICA method can maintain consistency with the results of commonly used activation region detection algorithms, it is difficult to use ICA to detect complex high-level activities of the brain.

发明内容Contents of the invention

本发明的目的在于，解决现有的脑激活区检测方法信噪比低的问题。The purpose of the present invention is to solve the problem of low signal-to-noise ratio in the existing brain activation area detection method.

本发明的目的是采用以下技术方案来实现的。The purpose of the present invention is achieved by adopting the following technical solutions.

一种脑激活区检测方法，包括以下步骤：A method for detecting brain activation areas, comprising the following steps:

步骤S1，获取fMRI数据；Step S1, obtaining fMRI data;

步骤S2，构建所述fMRI数据的二维邻域特征空间S；Step S2, constructing a two-dimensional neighborhood feature space S of the fMRI data;

步骤S3，采用均值漂移算法对所述特征空间S进行聚类搜索；Step S3, performing a cluster search on the feature space S by using the mean shift algorithm;

步骤S4，得到激活区检测结果。Step S4, obtaining the detection result of the active area.

本发明一较佳实施方式中，步骤S2进一步包括以下步骤：In a preferred embodiment of the present invention, step S2 further includes the following steps:

步骤S21，给定一组包含T个时间点的fMRI体数据V＝{Vt|t＝1,2,...,T}，其中Vt为时间点t的fMRI体数据，体素p的时间序列I(p)＝{I(p,t)|t＝1,2,...,T,p∈V}，其中I(p,t)为体素p在Vt中的信号强度；fMRI成像过程中，外部刺激为特定时间点的刺激信号，设定为刺激函数E(t)(t＝1,…,T)，则体素p的时间序列I(p)与刺激函数E(t)的相关系数r1(p)表示为：Step S21, given a set of fMRI volume data V={V t |t=1,2,...,T} containing T time points, where V t is the fMRI volume data at time point t, voxel p The time series I(p)={I(p,t)|t=1,2,...,T,p∈V}, where I(p,t) is the signal of voxel p in V t Intensity; during the fMRI imaging process, the external stimulus is the stimulus signal at a specific time point, which is set as the stimulus function E(t) (t=1,...,T), then the time series I(p) of voxel p and the stimulus function The correlation coefficient r 1 (p) of E(t) is expressed as:

rr 11 (( pp )) == ΣΣ tt == 11 TT (( II (( pp ,, tt )) -- II ‾‾ (( pp )) )) ·· (( EE. (( tt )) -- EE. ‾‾ )) [[ ΣΣ tt == 11 TT (( II (( pp ,, tt )) -- II ‾‾ (( pp )) )) 22 ΣΣ tt == 11 TT (( EE. (( tt )) -- EE. ‾‾ )) 22 ]] 11 // 22 -- -- -- (( 11 ))

其中，为体素p的信号强度均值，为刺激函数均值；in, is the mean signal intensity of voxel p, is the mean value of the stimulus function;

步骤S22，引入体素p周围邻域信息，定义体素p的邻域为N(p)，N(p)＝{pi|i＝1,2,…,n}，n为邻域N(p)内体素的个数，取每个体素邻域N(p)内所有体素的时间序列与刺激函数的相关系数的平均值，定义为表示为：Step S22, introduce the neighborhood information around voxel p, define the neighborhood of voxel p as N(p), N(p)={p i |i=1,2,...,n}, n is the neighborhood N The number of voxels in (p), take the average value of the correlation coefficient between the time series of all voxels in each voxel neighborhood N(p) and the stimulus function, defined as Expressed as:

rr 11 ‾‾ (( pp )) == 11 nno ΣΣ ii == 11 nno rr 11 (( pp ii )) -- -- -- (( 22 ))

构建fMRI数据的第一维特征 Constructing first-dimensional features of fMRI data

步骤S23，获取体素p与该体素邻域中的体素的相关系数r2(p)，表示为：Step S23, obtaining the correlation coefficient r 2 (p) between the voxel p and the voxels in the neighborhood of the voxel, expressed as:

rr 22 (( pp )) == 11 nno ΣΣ ii == 11 nno ΣΣ tt == 11 TT (( II (( pp ,, tt )) -- II ‾‾ (( pp )) )) ·· (( II (( pp ii ,, tt )) -- II ‾‾ (( pp ii )) )) [[ ΣΣ tt == 11 TT (( II (( pp ,, tt )) -- II ‾‾ (( pp )) )) 22 ΣΣ tt == 11 TT (( II (( pp ii ,, tt )) -- II ‾‾ (( pp ii )) )) 22 ]] 11 // 22 -- -- -- (( 33 ))

其中，为序列均值，为邻域中对应体素的均值，构建fMRI数据特征空间的第二维特征R2＝{r2(p)|p∈V}；in, is the sequence mean, Construct the second dimension feature R 2 ={r 2 (p)|p∈V} of the fMRI data feature space as the mean value of the corresponding voxels in the neighborhood;

步骤S24，构建所述特征空间S，其中，S＝{R1,R2}。Step S24, constructing the feature space S, where S={R 1 , R 2 }.

本发明一较佳实施方式中，在所述特征空间S中，给定空间中的采样点的核函数为k(x)和容许误差ε，所述步骤S3进一步包括以下步骤：In a preferred embodiment of the present invention, 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 further includes the following steps:

步骤S31，在所述特征空间S中任意选择初始搜索区域圆心O，半径为核宽h；Step S31, randomly select the center O of the initial search area in the feature space S, and the radius is the kernel width h;

步骤S32，在搜索区域内计算均值漂移Mh(x)：Step S32, calculate the mean shift M h (x) in the search area:

Mm hh (( xx )) == ΣΣ ii == 11 nno xx ii gg (( || || xx -- xx ii hh || || 22 )) ΣΣ ii == 11 nno gg (( || || xx -- xx ii hh || || 22 )) -- xx -- -- -- (( 22 ))

该向量总是指向密度增加的方向；This vector always points in the direction of increasing density;

步骤S33，如果均值漂移向量的模小于容许误差ε，||Mh(x)||<ε，迭代算法结束；否则利用下式(3)计算x得到新的圆心O′，返回执行步骤S32；Step S33, if the modulus of the mean shift vector is less than the allowable error ε, ||M h (x)||<ε, the iterative algorithm ends; otherwise, use the following formula (3) to calculate x to obtain a new center O', and return to step S32 ;

xx == &Sigma;&Sigma; ii == 11 nno xx ii gg (( || || xx -- xx ii hh || || 22 )) &Sigma;&Sigma; ii == 11 nno gg (( || || xx -- xx ii hh || || 22 )) -- -- -- (( 33 ))

本发明一较佳实施方式中，步骤S4进一步包括以下步骤：重复步骤S1至S3，直到||Mh(x)||<ε、公式(2)收敛到局部密度极大值点，x漂移到局部极大点时，将所述特征空间S中收敛到同一个点的样本点算作一类，得到所述激活区检测结果。In a preferred embodiment of the present invention, step S4 further includes the following steps: repeat steps S1 to S3 until ||M h (x)||<ε, formula (2) converges to the local density maximum point, and x drifts When the local maximum point is reached, the sample points in the feature space S that converge to the same point are counted as one class to obtain the activation region detection result.

本发明一较佳实施方式中，步骤S1进一步包括以下步骤：制定实验刺激方案，采集fMRI实验数据并对所述实验数据进行预处理。In a preferred embodiment of the present invention, step S1 further includes the following steps: formulating an experimental stimulation scheme, collecting fMRI experimental data and preprocessing the experimental data.

一种脑激活区检测装置，其包括：A device for detecting brain activation areas, comprising:

数据抓取单元，用于获取fMRI数据；A data capture unit for acquiring fMRI data;

构建单元，用于构建所述fMRI数据的二维邻域特征空间S；A construction unit for constructing the two-dimensional neighborhood feature space S of the fMRI data;

搜索单元，采用均值漂移算法对所述特征空间S进行聚类搜索；以及A search unit, which uses a mean shift algorithm to perform a cluster search on the feature space S; and

输出单元，用于得到激活区检测结果。The output unit is used to obtain the detection result of the activation area.

本发明一较佳实施方式中，所述构建单元进一步包括以下子单元：In a preferred embodiment of the present invention, the construction unit further includes the following subunits:

数据选择子单元，用于给定一组包含T个时间点的fMRI体数据V＝{Vt|t＝1,2,...,T}，其中Vt为时间点t的fMRI体数据，体素p的时间序列I(p)＝{I(p,t)|t＝1,2,...,T,p∈V}，其中I(p,t)为体素p在Vt中的信号强度；fMRI成像过程中，外部刺激为特定时间点的刺激信号，设定为刺激函数E(t)(t＝1,…,T)，则体素p的时间序列I(p)与刺激函数E(t)的相关系数r1(p)表示为：The data selection subunit is used to give a set of fMRI volume data V={V t |t=1,2,...,T} containing T time points, where V t is the fMRI volume data at time point t , the time series of voxel p I(p)={I(p,t)|t=1,2,...,T,p∈V}, where I(p,t) is the voxel p in V The signal intensity in t ; during the fMRI imaging process, the external stimulus is the stimulus signal at a specific time point, which is set as the stimulus function E(t) (t=1,...,T), then the time series I(p ) and the correlation coefficient r 1 (p) of the stimulus function E(t) is expressed as:

rr 11 (( pp )) == &Sigma;&Sigma; tt == 11 TT (( II (( pp ,, tt )) -- II &OverBar;&OverBar; (( pp )) )) &CenterDot;&CenterDot; (( EE. (( tt )) -- EE. &OverBar;&OverBar; )) &lsqb;&lsqb; &Sigma;&Sigma; tt == 11 TT (( II (( pp ,, tt )) -- II &OverBar;&OverBar; (( pp )) )) 22 &Sigma;&Sigma; tt == 11 TT (( EE. (( tt )) -- EE. &OverBar;&OverBar; )) 22 &rsqb;&rsqb; 11 // 22 -- -- -- (( 11 ))

其中，为体素p的信号强度均值，为刺激函数均值；in, is the mean signal intensity of voxel p, is the mean value of the stimulus function;

第一维特征构建子单元，用于引入体素p周围邻域信息，定义体素p的邻域为N(p)，N(p)＝{pi|i＝1,2,…,n}，n为邻域N(p)内体素的个数，取每个体素邻域N(p)内所有体素的时间序列与刺激函数的相关系数的平均值，定义为表示为：The first dimension feature construction subunit is used to introduce the neighborhood information around voxel p, and define the neighborhood of voxel p as N(p), N(p)={p i |i=1,2,...,n }, n is the number of voxels in the neighborhood N(p), taking the average value of the correlation coefficient between the time series and the stimulus function of all voxels in the neighborhood N(p) of each voxel, defined as Expressed as:

rr 11 &OverBar;&OverBar; (( pp )) == 11 nno &Sigma;&Sigma; ii == 11 nno rr 11 (( pp ii )) -- -- -- (( 22 ))

构建fMRI数据的第一维特征 Constructing first-dimensional features of fMRI data

第二维特征构建子单元，用于获取体素p与该体素邻域中的体素的相关系数r2(p)，表示为：The second-dimensional feature construction subunit is used to obtain the correlation coefficient r 2 (p) between the voxel p and the voxels in the neighborhood of this voxel, expressed as:

rr 22 (( pp )) == 11 nno &Sigma;&Sigma; ii == 11 nno &Sigma;&Sigma; tt == 11 TT (( II (( pp ,, tt )) -- II &OverBar;&OverBar; (( pp )) )) &CenterDot;&CenterDot; (( II (( pp ii ,, tt )) -- II &OverBar;&OverBar; (( pp ii )) )) &lsqb;&lsqb; &Sigma;&Sigma; tt == 11 TT (( II (( pp ,, tt )) -- II &OverBar;&OverBar; (( pp )) )) 22 &Sigma;&Sigma; tt == 11 TT (( II (( pp ii ,, tt )) -- II &OverBar;&OverBar; (( pp ii )) )) 22 &rsqb;&rsqb; 11 // 22 -- -- -- (( 33 ))

构建fMRI数据特征空间的第二维特征R2＝{r2(p)|p∈V}，其中，为序列均值，为邻域中对应体素的均值；以及Construct the second dimension feature R 2 ={r 2 (p)|p∈V} of fMRI data feature space, where, is the sequence mean, is the mean value of the corresponding voxels in the neighborhood; and

特征空间构建子单元，用于构建所述特征空间S，其中，S＝{R1,R2}。The feature space construction subunit is used to construct the feature space S, where S={R 1 , R 2 }.

本发明一较佳实施方式中，在所述特征空间S中，给定空间中的采样点的核函数为k(x)和容许误差ε，所述搜索单元进一步包括以下子单元：In a preferred embodiment of the present invention, 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 search unit further includes the following subunits:

条件设置子单元，用于在所述特征空间S中任意选择初始搜索区域圆心O，半径为核宽h；The condition setting subunit is used to arbitrarily select the center O of the initial search area in the feature space S, and the radius is the kernel width h;

计算子单元，用于在搜索区域内计算均值漂移Mh(x)：Calculation subunit for computing the mean shift M h (x) within the search area:

Mm hh (( xx )) == &Sigma;&Sigma; ii == 11 nno xx ii gg (( || || xx -- xx ii hh || || 22 )) &Sigma;&Sigma; ii == 11 nno gg (( || || xx -- xx ii hh || || 22 )) -- xx -- -- -- (( 22 ))

该向量总是指向密度增加的方向；This vector always points in the direction of increasing density;

判断子单元，用于判断：如果均值漂移向量的模小于容许误差ε，||Mh(x)||<ε，迭代算法结束；否则利用下式(3)计算x得到新的圆心O′，返回给所述计算子单元继续计算；The judging sub-unit is used for judging: if the modulus of the mean shift vector is less than the allowable error ε, ||M h (x)||<ε, the iterative algorithm ends; otherwise, use the following formula (3) to calculate x to obtain a new center O' , returning to the calculation subunit to continue calculation;

xx == &Sigma;&Sigma; ii == 11 nno xx ii gg (( || || xx -- xx ii hh || || 22 )) &Sigma;&Sigma; ii == 11 nno gg (( || || xx -- xx ii hh || || 22 )) -- -- -- (( 33 ))

本发明一较佳实施方式中，所述输出单元还用于将所述特征空间S中收敛到同一个点的样本点算作一类，得到所述激活区检测结果。In a preferred embodiment of the present invention, the output unit is further configured to count the sample points converging to the same point in the feature space S as one class to obtain the activation region detection result.

本发明一较佳实施方式中，所述数据抓取单元还用于对所述fMRI实验数据进行预处理。In a preferred embodiment of the present invention, the data capture unit is also used to preprocess the fMRI experiment data.

相较于现有技术，本发明提供的脑激活区检测方法和装置不但可以检测到活动较强的激活体素，还可以敏感地检测到那些由活动微弱但空间域高度相关的体素构成的脑激活区，具有良好的抗噪能力和高灵敏度。Compared with the prior art, the method and device for detecting brain activation regions provided by the present invention can not only detect activation voxels with strong activity, but also sensitively detect those voxels composed of weak activity but highly correlated spatial domains. Brain activation area, with good noise immunity and high sensitivity.

上述说明仅是本发明技术方案的概述，为了能够更清楚了解本发明的技术手段，而可依照说明书的内容予以实施，并且为了让本发明的上述和其他目的、特征和优点能够更明显易懂，以下特举较佳实施例，并配合附图，详细说明如下。The above description is only an overview of the technical solution of the present invention. In order to better understand the technical means of the present invention, it can be implemented according to the contents of the description, and in order to make the above and other purposes, features and advantages of the present invention more obvious and understandable , the following preferred embodiments are specifically cited below, and are described in detail as follows in conjunction with the accompanying drawings.

附图说明Description of drawings

图1是本发明第一实施例提供的脑激活区检测方法的流程示意图。FIG. 1 is a schematic flowchart of a method for detecting brain activation regions provided by the first embodiment of the present invention.

图2是本发明第二实施例提供的脑激活区检测装置的结构示意图。Fig. 2 is a schematic structural diagram of a brain activation area detection device provided by a second embodiment of the present invention.

图3是本发明第二实施例提供的脑激活区检测装置的构建单元的结构示意图。Fig. 3 is a schematic structural diagram of the construction units of the brain activation area detection device provided by the second embodiment of the present invention.

图4是本发明第二实施例提供的脑激活区检测装置的搜索单元的结构示意图。Fig. 4 is a schematic structural diagram of the search unit of the brain activation area detection device provided by the second embodiment of the present invention.

具体实施方式detailed description

为了便于理解本发明，下面将参照相关附图对本发明进行更全面的描述。附图中给出了本发明的较佳实施方式。但是，本发明可以以许多不同的形式来实现，并不限于本文所描述的实施方式。相反地，提供这些实施方式的目的是使对本发明的公开内容理解的更加透彻全面。In order to facilitate the understanding of the present invention, the present invention will be described more fully below with reference to the associated drawings. Preferred embodiments of the invention are shown in the accompanying drawings. However, the present invention can be embodied in many different forms and is not limited to the embodiments described herein. On the contrary, the purpose of providing these embodiments is to make the disclosure of the present invention more thorough and comprehensive.

除非另有定义，本文所使用的所有的技术和科学术语与属于本发明的技术领域的技术人员通常理解的含义相同。本文中在本发明的说明书中所使用的术语只是为了描述具体的实施方式的目的，不是旨在于限制本发明。本文所使用的术语“及/或”包括一个或多个相关的所列项目的任意的和所有的组合。Unless otherwise defined, all technical and scientific terms used herein have the same meaning as commonly understood by one of ordinary skill in the technical field of the invention. The terminology used herein in the description of the present invention is only for the purpose of describing specific embodiments, and is not intended to limit the present invention. As used herein, the term "and/or" includes any and all combinations of one or more of the associated listed items.

请参阅图1，图1是本发明第一实施例提供的脑激活区检测方法的流程示意图，包括以下步骤：Please refer to FIG. 1. FIG. 1 is a schematic flowchart of a method for detecting brain activation regions provided by the first embodiment of the present invention, including the following steps:

步骤S1，获取fMRI数据。Step S1, obtaining fMRI data.

具体地，先制定实验刺激方案，采集fMRI实验数据并对所述实验数据进行预处理。预处理目的是检测和修复由磁共振扫描仪或者被试个人在进行数据采集时产生的伪迹或者为接下来的数据处理做准备。Specifically, an experimental stimulation scheme is formulated first, fMRI experimental data is collected, and the experimental data is preprocessed. The purpose of preprocessing is to detect and repair the artifacts produced by the MRI scanner or individual subjects during data acquisition or to prepare for the next data processing.

步骤S2，构建fMRI数据的二维邻域特征空间S。Step S2, constructing a two-dimensional neighborhood feature space S of fMRI data.

具体地，包括以下步骤：Specifically, the following steps are included:

步骤S21，给定一组包含T个时间点的fMRI体数据V＝{Vt|t＝1,2,...,T}，其中Vt为时间点t的fMRI体数据，体素p的时间序列I(p)＝{I(p,t)|t＝1,2,...,T,p∈V}，其中I(p,t)为体素p在Vt中的信号强度；fMRI成像过程中，外部刺激为特定时间点的刺激信号，设定为刺激函数E(t)(t＝1,…,T)，则体素p的时间序列I(p)与刺激函数E(t)的相关系数r1(p)可由互相关分析获得，表示为：Step S21, given a set of fMRI volume data V={V t |t=1,2,...,T} containing T time points, where V t is the fMRI volume data at time point t, voxel p The time series I(p)={I(p,t)|t=1,2,...,T,p∈V}, where I(p,t) is the signal of voxel p in V t Intensity; during the fMRI imaging process, the external stimulus is the stimulus signal at a specific time point, which is set as the stimulus function E(t) (t=1,...,T), then the time series I(p) of voxel p and the stimulus function The correlation coefficient r 1 (p) of E(t) can be obtained by cross-correlation analysis, expressed as:

rr 11 (( pp )) == &Sigma;&Sigma; tt == 11 TT (( II (( pp ,, tt )) -- II &OverBar;&OverBar; (( pp )) )) &CenterDot;&Center Dot; (( EE. (( tt )) -- EE. &OverBar;&OverBar; )) &lsqb;&lsqb; &Sigma;&Sigma; tt == 11 TT (( II (( pp ,, tt )) -- II &OverBar;&OverBar; (( pp )) )) 22 &Sigma;&Sigma; tt == 11 TT (( EE. (( tt )) -- EE. &OverBar;&OverBar; )) 22 &rsqb;&rsqb; 11 // 22 -- -- -- (( 11 ))

其中，为体素p的信号强度均值，为刺激函数均值。相关系数r1(p)的大小反应了该体素参与任务刺激的相关程度。in, is the mean signal intensity of voxel p, is the mean value of the stimulus function. The magnitude of the correlation coefficient r 1 (p) reflects the correlation degree of the voxel participating in the task stimulus.

步骤S22，引入体素p周围邻域信息，定义体素p的邻域为N(p)，N(p)＝{pi|i＝1,2,…,n}，n为邻域N(p)内体素的个数，取每个体素邻域N(p)内所有体素的时间序列与刺激函数的相关系数的平均值，定义为表示为：Step S22, introduce the neighborhood information around voxel p, define the neighborhood of voxel p as N(p), N(p)={p i |i=1,2,...,n}, n is the neighborhood N The number of voxels in (p), take the average value of the correlation coefficient between the time series of all voxels in each voxel neighborhood N(p) and the stimulus function, defined as Expressed as:

rr 11 &OverBar;&OverBar; (( pp )) == 11 nno &Sigma;&Sigma; ii == 11 nno rr 11 (( pp ii )) -- -- -- (( 22 ))

构建fMRI数据的第一维特征 Constructing first-dimensional features of fMRI data

步骤S23，获取体素p与该体素邻域中的体素的相关系数r2(p)，相关系数r2(p)反映了fMRI数据邻域空间的一致性，表示为：Step S23, obtaining the correlation coefficient r 2 (p) between the voxel p and the voxels in the voxel neighborhood, the correlation coefficient r 2 ( p) reflects the consistency of the fMRI data neighborhood space, expressed as:

rr 22 (( pp )) == 11 nno &Sigma;&Sigma; ii == 11 nno &Sigma;&Sigma; tt == 11 TT (( II (( pp ,, tt )) -- II &OverBar;&OverBar; (( pp )) )) &CenterDot;&Center Dot; (( II (( pp ii ,, tt )) -- II &OverBar;&OverBar; (( pp ii )) )) &lsqb;&lsqb; &Sigma;&Sigma; tt == 11 TT (( II (( pp ,, tt )) -- II &OverBar;&OverBar; (( pp )) )) 22 &Sigma;&Sigma; tt == 11 TT (( II (( pp ii ,, tt )) -- II &OverBar;&OverBar; (( pp ii )) )) 22 &rsqb;&rsqb; 11 // 22 -- -- -- (( 33 ))

其中，为序列均值，为邻域中对应体素的均值，构建fMRI数据特征空间的第二维特征R2＝{r2(p)|p∈V}；in, is the sequence mean, Construct the second dimension feature R 2 ={r 2 (p)|p∈V} of the fMRI data feature space as the mean value of the corresponding voxels in the neighborhood;

步骤S24，构建所述特征空间S，其中，S＝{R1,R2}。Step S24, constructing the feature space S, where S={R 1 , R 2 }.

在其他实施例中，n的取值范围，可以根据实际情况设置，可以是8邻域或者4邻域等等。In other embodiments, the value range of n can be set according to the actual situation, and can be 8 neighborhoods or 4 neighborhoods and so on.

步骤S3，采用均值漂移算法对特征空间S进行聚类搜索。Step S3, cluster search is performed on the feature space S by using the mean shift algorithm.

具体地，在所述特征空间S中，给定空间中的采样点的核函数为k(x)和容许误差ε，且进一步包括以下步骤：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 following steps:

步骤S31，在所述特征空间S中任意选择初始搜索区域圆心O，半径为核宽h；Step S31, randomly select the center O of the initial search area in the feature space S, and the radius is the kernel width h;

步骤S32，在搜索区域内计算均值漂移Mh(x)：Step S32, calculate the mean shift M h (x) in the search area:

Mm hh (( xx )) == &Sigma;&Sigma; ii == 11 nno xx ii gg (( || || xx -- xx ii hh || || 22 )) &Sigma;&Sigma; ii == 11 nno gg (( || || xx -- xx ii hh || || 22 )) -- xx -- -- -- (( 22 ))

该向量总是指向密度增加的方向；This vector always points in the direction of increasing density;

步骤S33，如果均值漂移向量的模小于容许误差ε，||Mh(x)||<ε，迭代算法结束；否则利用下式(3)计算x得到新的圆心O′，返回执行步骤S32；Step S33, if the modulus of the mean shift vector is less than the allowable error ε, ||M h (x)||<ε, the iterative algorithm ends; otherwise, use the following formula (3) to calculate x to obtain a new center O', and return to step S32 ;

xx == &Sigma;&Sigma; ii == 11 nno xx ii gg (( || || xx -- xx ii hh || || 22 )) &Sigma;&Sigma; ii == 11 nno gg (( || || xx -- xx ii hh || || 22 )) -- -- -- (( 33 ))

在步骤S31中，关键参数为核宽h，其决定了搜索区域的大小，h取值的不同会影响算法的结果，在其他实施例中，可根据构成的特征空间S的不同，核宽h的取值也会不同。In step S31, the key parameter is the kernel width h, which determines the size of the search area, and the different values of h will affect the results of the algorithm. In other embodiments, the kernel width h can be determined according to the different feature spaces S formed. The value of will also be different.

步骤S4，得到激活区检测结果。Step S4, obtaining the detection result of the active area.

具体地，包括：重复步骤S1至S3(在本实施例中是指重复步骤S1、S21至S23、S31至S33)，直到||Mh(x)||<ε、公式(2)收敛到局部密度极大值点，x漂移到局部极大点时，将所述特征空间S中收敛到同一个点的样本点算作一类，得到VN-MSC算法(Voxel’sNeighborhood-basedMean-ShiftClustering,VN-MSC)的激活区检测结果。Specifically, it includes: repeating steps S1 to S3 (referring to repeating steps S1, S21 to S23, and S31 to S33 in this embodiment), until ||M h (x)||<ε, formula (2) converges to Local density maximum point, when x drifts to the local maximum point, the sample points that converge to the same point in the feature space S are counted as a class, and the VN-MSC algorithm (Voxel'sNeighborhood-basedMean-ShiftClustering, Activation region detection results of VN-MSC).

本实施例提供的脑激活区检测方法利用了大脑同一区域的体素具有相似性质的特点，即脑功能区都是呈块状分布，是由一些fMRI数据空间相连的体素所组成的团簇，而非单个体素形成的点状分布，因此充分利用了fMRI数据体素空间的邻域相关性特征来构建特征空间，并结合均值漂移聚类算法提取不同刺激条件下的脑激活区，不但可以检测到活动较强的激活体素，还可以敏感地检测到那些由活动微弱但空间域高度相关的体素构成的脑激活区，具有良好的抗噪能力和高灵敏度。The brain activation region detection method provided in this example takes advantage of the fact that the voxels in the same region of the brain have similar properties, that is, the brain functional regions are all distributed in blocks, and are clusters composed of some fMRI data spatially connected voxels , rather than a point-like distribution formed by a single voxel, so the neighborhood correlation feature of the fMRI data voxel space is fully utilized to construct the feature space, and the mean shift clustering algorithm is combined to extract the brain activation area under different stimulation conditions, not only It can detect activated voxels with strong activity, and can also sensitively detect those brain activation areas composed of weakly active but highly correlated spatial domain voxels, with good anti-noise ability and high sensitivity.

请参阅图2，图2是本发明第二实施例提供的脑激活区检测装置100的结构示意图。Please refer to FIG. 2 . FIG. 2 is a schematic structural diagram of a brain activation region detection device 100 according to a second embodiment of the present invention.

脑激活区检测装置100包括：数据抓取单元10，用于获取fMRI数据；构建单元20，用于构建fMRI数据的二维邻域特征空间S；搜索单元30，采用均值漂移算法对所述特征空间S进行聚类搜索；以及输出单元40，用于得到激活区检测结果。The brain activation region detection device 100 includes: a data capture unit 10, used to obtain fMRI data; a construction unit 20, used to construct a two-dimensional neighborhood feature space S of fMRI data; a search unit 30, which uses a mean shift algorithm to search the feature The space S is clustered and searched; and an output unit 40 is used to obtain the detection result of the activation area.

请参阅图3，图3是构建单元20的结构示意图。构建单元20进一步包括以下子单元：Please refer to FIG. 3 , which is a schematic structural diagram of the construction unit 20 . Construction unit 20 further includes the following subunits:

数据选择子单元201，用于给定一组包含T个时间点的fMRI体数据V＝{Vt|t＝1,2,...,T}，其中Vt为时间点t的fMRI体数据，体素p的时间序列I(p)＝{I(p,t)|t＝1,2,...,T,p∈V}，其中I(p,t)为体素p在Vt中的信号强度；fMRI成像过程中，外部刺激为特定时间点的刺激信号，设定为刺激函数E(t)(t＝1,…,T)，则体素p的时间序列I(p)与刺激函数E(t)的相关系数r1(p)表示为：The data selection subunit 201 is used to give a set of fMRI volume data V={V t |t=1,2,...,T} containing T time points, where V t is the fMRI volume at time point t Data, the time series of voxel p I(p)={I(p,t)|t=1,2,...,T,p∈V}, where I(p,t) is the voxel p in The signal intensity in V t ; in the fMRI imaging process, the external stimulus is the stimulus signal at a specific time point, which is set as the stimulus function E(t) (t=1,...,T), then the time series I of voxel p( The correlation coefficient r 1 (p) between p) and the stimulus function E(t) is expressed as:

rr 11 (( pp )) == &Sigma;&Sigma; tt == 11 TT (( II (( pp ,, tt )) -- II &OverBar;&OverBar; (( pp )) )) &CenterDot;&Center Dot; (( EE. (( tt )) -- EE. &OverBar;&OverBar; )) &lsqb;&lsqb; &Sigma;&Sigma; tt == 11 TT (( II (( pp ,, tt )) -- II &OverBar;&OverBar; (( pp )) )) 22 &Sigma;&Sigma; tt == 11 TT (( EE. (( tt )) -- EE. &OverBar;&OverBar; )) 22 &rsqb;&rsqb; 11 // 22 -- -- -- (( 11 ))

其中，为体素p的信号强度均值，为刺激函数均值；in, is the mean signal intensity of voxel p, is the mean value of the stimulus function;

第一维特征构建子单元202，用于引入体素p周围邻域信息，定义体素p的邻域为N(p)，N(p)＝{pi|i＝1,2,…,n}，n为邻域N(p)内体素的个数，取每个体素邻域N(p)内所有体素的时间序列与刺激函数的相关系数的平均值，定义为表示为：The first dimension feature construction subunit 202 is used to introduce the neighborhood information around voxel p, and define the neighborhood of voxel p as N(p), N(p)={p i |i=1,2,..., n}, n is the number of voxels in the neighborhood N(p), take the average value of the correlation coefficient between the time series and the stimulus function of all voxels in the neighborhood N(p) of each voxel, defined as Expressed as:

rr 11 &OverBar;&OverBar; (( pp )) == 11 nno &Sigma;&Sigma; ii == 11 nno rr 11 (( pp ii )) -- -- -- (( 22 ))

构建fMRI数据的第一维特征 Constructing first-dimensional features of fMRI data

第二维特征构建子单元203，用于获取体素p与该体素邻域中的体素的相关系数r2(p)，表示为：The second-dimensional feature construction subunit 203 is used to obtain the correlation coefficient r 2 (p) between the voxel p and the voxels in the neighborhood of the voxel, expressed as:

rr 22 (( pp )) == 11 nno &Sigma;&Sigma; ii == 11 nno &Sigma;&Sigma; tt == 11 TT (( II (( pp ,, tt )) -- II &OverBar;&OverBar; (( pp )) )) &CenterDot;&CenterDot; (( II (( pp ii ,, tt )) -- II &OverBar;&OverBar; (( pp ii )) )) &lsqb;&lsqb; &Sigma;&Sigma; tt == 11 TT (( II (( pp ,, tt )) -- II &OverBar;&OverBar; (( pp )) )) 22 &Sigma;&Sigma; tt == 11 TT (( II (( pp ii ,, tt )) -- II &OverBar;&OverBar; (( pp ii )) )) 22 &rsqb;&rsqb; 11 // 22 -- -- -- (( 33 ))

构建fMRI数据特征空间的第二维特征R2＝{r2(p)|p∈V}，其中，为序列均值，为邻域中对应体素的均值；以及Construct the second dimension feature R 2 ={r 2 (p)|p∈V} of fMRI data feature space, where, is the sequence mean, is the mean value of the corresponding voxels in the neighborhood; and

特征空间构建子单元204，用于构建所述特征空间S，其中，S＝{R1,R2}。The feature space construction subunit 204 is configured to construct the feature space S, where S={R 1 , R 2 }.

在特征空间S中，给定空间中的采样点的核函数为k(x)和容许误差ε，请参阅图4，搜索单元30进一步包括以下子单元：In the feature space S, the kernel function of the sampling point in the given space is k(x) and the allowable error ε, referring to Fig. 4, the search unit 30 further includes the following subunits:

条件设置子单元301，用于在所述特征空间S中任意选择初始搜索区域，圆心为O，半径为核宽h；The condition setting subunit 301 is used to arbitrarily select an initial search area in the feature space S, the center of the circle is O, and the radius is the kernel width h;

计算子单元302，用于在搜索区域内计算均值漂移Mh(x)：Calculation subunit 302, for calculating the mean shift M h (x) in the search area:

Mm hh (( xx )) == &Sigma;&Sigma; ii == 11 nno xx ii gg (( || || xx -- xx ii hh || || 22 )) &Sigma;&Sigma; ii == 11 nno gg (( || || xx -- xx ii hh || || 22 )) -- xx -- -- -- (( 22 ))

该向量总是指向密度增加的方向；This vector always points in the direction of increasing density;

判断子单元303，用于判断：如果均值漂移向量的模小于容许误差ε，||Mh(x)||<ε，迭代算法结束；否则利用下式(3)计算x得到新的圆心O′，返回给所述计算子单元302继续计算；Judgment subunit 303, used to judge: if the modulus of the mean shift vector is less than the allowable error ε, ||M h (x)||<ε, the iterative algorithm ends; otherwise, use the following formula (3) to calculate x to obtain a new circle center O ', return to the calculation subunit 302 to continue calculation;

xx == &Sigma;&Sigma; ii == 11 nno xx ii gg (( || || xx -- xx ii hh || || 22 )) &Sigma;&Sigma; ii == 11 nno gg (( || || xx -- xx ii hh || || 22 )) -- -- -- (( 33 ))

输出单元40还用于将公式(2)收敛到局部密度极大值点时，特征空间S中收敛到同一个点的样本点算作一类，得到激活区检测结果。The output unit 40 is also used to calculate the sample points that converge to the same point in the feature space S when the formula (2) converges to the local density maximum point as one class, and obtain the activation region detection result.

数据抓取单元10还用于对fMRI实验数据进行预处理。The data capture unit 10 is also used for preprocessing the fMRI experiment data.

本实施例提供的脑激活区检测装置利用了大脑同一区域的体素具有相似性质的特点，即脑功能区都是呈块状分布，是由一些fMRI数据空间相连的体素所组成的团簇，而非单个体素形成的点状分布，因此充分利用了fMRI数据体素空间的邻域相关性特征来构建特征空间，并结合均值漂移聚类算法提取不同刺激条件下的脑激活区，不但可以检测到活动较强的激活体素，还可以敏感地检测到那些由活动微弱但空间域高度相关的体素构成的脑激活区，具有良好的抗噪能力和高灵敏度。The brain activation region detection device provided in this embodiment takes advantage of the fact that the voxels in the same region of the brain have similar properties, that is, the brain functional regions are all distributed in blocks, and are clusters composed of voxels connected in some fMRI data spaces , rather than the point-like distribution formed by a single voxel, so the neighborhood correlation feature of the fMRI data voxel space is fully utilized to construct the feature space, and the mean shift clustering algorithm is combined to extract the brain activation area under different stimulation conditions, not only It can detect activated voxels with strong activity, and can also sensitively detect those brain activation areas composed of weakly active but highly correlated spatial domain voxels, with good anti-noise ability and high sensitivity.

以上所述实施例仅表达了本发明的几种实施方式，其描述较为具体和详细，但并不能因此而理解为对本发明专利范围的限制。应当指出的是，对于本领域的普通技术人员来说，在不脱离本发明构思的前提下，还可以做出若干变形和改进，这些都属于本发明的保护范围。因此，本发明专利的保护范围应以所附权利要求为准。The above-mentioned embodiments only express several implementation modes of the present invention, and the description thereof is relatively specific and detailed, but should not be construed as limiting the patent scope of the present invention. It should be pointed out that those skilled in the art can make several modifications and improvements without departing from the concept of the present invention, and these all belong to the protection scope of the present invention. Therefore, the protection scope of the patent for the present invention should be based on 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 pointst1, 2., T }, where VtFor 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 VtSignal 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 as1(p) is represented by:
r 1 ( p ) = &Sigma; t = 1 T ( I ( p , t ) - I &OverBar; ( p ) ) &CenterDot; ( E ( t ) - E &OverBar; ) &lsqb; &Sigma; t = 1 T ( I ( p , t ) - I &OverBar; ( p ) ) 2 &Sigma; t = 1 T ( E ( t ) - E &OverBar; ) 2 &rsqb; 1 / 2 - - - ( 1 )
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) ═ pi1, 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:
r 1 &OverBar; ( p ) = 1 n &Sigma; i = 1 n r 1 ( p i ) - - - ( 2 )
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 obtained2(p), expressed as:
r 2 ( p ) = 1 n &Sigma; i = 1 n &Sigma; t = 1 T ( I ( p , t ) - I &OverBar; ( p ) ) &CenterDot; ( I ( p i , t ) - I &OverBar; ( p i ) ) &lsqb; &Sigma; t = 1 T ( I ( p , t ) - I &OverBar; ( p ) ) 2 &Sigma; t = 1 T ( I ( p i , t ) - I &OverBar; ( p i ) ) 2 &rsqb; 1 / 2 - - - ( 3 )
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 neighborhood2＝{r2(p)|p∈V}；
Step S24, constructing the feature space S, where S ═ { R ═ R1,R2}。
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 areah(x)：
M h ( x ) = &Sigma; i = 1 n x i g ( | | x - x i h | | 2 ) &Sigma; i = 1 n g ( | | x - x i h | | 2 ) - x - - - ( 2 )
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, | Mh(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 = &Sigma; i = 1 n x i g ( | | x - x i h | | 2 ) &Sigma; i = 1 n g ( | | x - x i h | | 2 ) - - - ( 3 )
4. the detection method according to any one of claims 1 to 3, wherein the step S4 includes: repeating steps S1-S3 until Mh(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 pointst1, 2., T }, where VtFor 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 VtSignal 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 as1(p) is represented by:
r 1 ( p ) = &Sigma; t = 1 T ( I ( p , t ) - I &OverBar; ( p ) ) &CenterDot; ( E ( t ) - E &OverBar; ) &lsqb; &Sigma; t = 1 T ( I ( p , t ) - I &OverBar; ( p ) ) 2 &Sigma; t = 1 T ( E ( t ) - E &OverBar; ) 2 &rsqb; 1 / 2 - - - ( 1 )
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) ═ pi1, 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:
r 1 &OverBar; ( p ) = 1 n &Sigma; i = 1 n r 1 ( p i ) - - - ( 2 )
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 neighborhood2(p), expressed as:
r 2 ( p ) = 1 n &Sigma; i = 1 n &Sigma; t = 1 T ( I ( p , t ) - I &OverBar; ( p ) ) &CenterDot; ( I ( p i , t ) - I &OverBar; ( p i ) ) &lsqb; &Sigma; t = 1 T ( I ( p , t ) - I &OverBar; ( p ) ) 2 &Sigma; t = 1 T ( I ( p i , t ) - I &OverBar; ( p i ) ) 2 &rsqb; 1 / 2 - - - ( 3 )
constructing a second dimension feature R of the fMRI data feature space2＝{r2(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 ═ R1,R2}。
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 areah(x)：
M h ( x ) = &Sigma; i = 1 n x i g ( | | x - x i h | | 2 ) &Sigma; i = 1 n g ( | | x - x i h | | 2 ) - x - - - ( 2 )
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, | Mh(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 = &Sigma; i = 1 n x i g ( | | x - x i h | | 2 ) &Sigma; i = 1 n g ( | | x - x i h | | 2 ) - - - ( 3 )
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.
CN201510862020.4A 2015-12-01 2015-12-01 A kind of brain activation area detection method and device Active CN105513058B (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201510862020.4A CN105513058B (en) 2015-12-01 2015-12-01 A kind of brain activation area detection method and device

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201510862020.4A CN105513058B (en) 2015-12-01 2015-12-01 A kind of brain activation area detection method and device

Publications (2)

Publication Number Publication Date
CN105513058A true CN105513058A (en) 2016-04-20
CN105513058B CN105513058B (en) 2019-04-23

Family

ID=55721015

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201510862020.4A Active CN105513058B (en) 2015-12-01 2015-12-01 A kind of brain activation area detection method and device

Country Status (1)

Country Link
CN (1) CN105513058B (en)

Cited By (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN106127726A (en) * 2016-05-18 2016-11-16 中国石油大学(华东) A kind of non-feature extraction and non-parametric 3D rendering registration new algorithm
CN107194934A (en) * 2017-05-08 2017-09-22 西安交通大学 A kind of brain active region detection method based on correlation analysis
CN116583221A (en) * 2021-09-17 2023-08-11 皇家飞利浦有限公司 Determination of subject-specific hemodynamic response functions

Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN101292871A (en) * 2007-04-25 2008-10-29 中国科学院自动化研究所 A Method for Extracting Brain Activation Regions in Magnetic Resonance Imaging Based on Pattern Recognition Classification
CN102508184A (en) * 2011-10-26 2012-06-20 中国科学院自动化研究所 Brain function active region detection method based on moving average time series models
CN103961103A (en) * 2014-05-07 2014-08-06 大连理工大学 Method for performing phase correction on ICA estimation components of plural fMRI data

Patent Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN101292871A (en) * 2007-04-25 2008-10-29 中国科学院自动化研究所 A Method for Extracting Brain Activation Regions in Magnetic Resonance Imaging Based on Pattern Recognition Classification
CN102508184A (en) * 2011-10-26 2012-06-20 中国科学院自动化研究所 Brain function active region detection method based on moving average time series models
CN103961103A (en) * 2014-05-07 2014-08-06 大连理工大学 Method for performing phase correction on ICA estimation components of plural fMRI data

Non-Patent Citations (3)

* Cited by examiner, † Cited by third party
Title
JAMES S.HYDE ET AL.: "Cross-Correlation:An fMRI Signal-Processing Strategy", 《NEUROIMAGE》 *
LEO AI ET AL.: "Application of mean-shift clustering to Blood oxygen level dependent functional MRI activation detection", 《BMC MEDICAL IMAGING》 *
张波 等: "一种改进ICA算法在脑功能区提取中的应用", 《计算机仿真》 *

Cited By (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN106127726A (en) * 2016-05-18 2016-11-16 中国石油大学(华东) A kind of non-feature extraction and non-parametric 3D rendering registration new algorithm
CN106127726B (en) * 2016-05-18 2019-09-06 中国石油大学(华东) A Non-Feature Extraction and Non-parametric 3D Image Registration Method
CN107194934A (en) * 2017-05-08 2017-09-22 西安交通大学 A kind of brain active region detection method based on correlation analysis
CN116583221A (en) * 2021-09-17 2023-08-11 皇家飞利浦有限公司 Determination of subject-specific hemodynamic response functions

Also Published As

Publication number Publication date
CN105513058B (en) 2019-04-23

Similar Documents

Publication Publication Date Title
Edupuganti et al. Uncertainty quantification in deep MRI reconstruction
Marquand et al. Bayesian multi-task learning for decoding multi-subject neuroimaging data
US11776679B2 (en) Methods for risk map prediction in AI-based MRI reconstruction
Dimitriadis et al. Improving the reliability of network metrics in structural brain networks by integrating different network weighting strategies into a single graph
CN103957784B (en) To the method that brain function MR data processes
CN104715241B (en) Tensor decomposition-based fMRI feature extraction and identification method
Anitha et al. Automated detection of white matter lesions in MRI brain images using spatio-fuzzy and spatio-possibilistic clustering models
CN115005798B (en) Brain image feature extraction method based on continuous edge functional connection
Wismüller et al. Large-scale nonlinear Granger causality: A data-driven, multivariate approach to recovering directed networks from short time-series data
Bhatele et al. Glioma Segmentation and Classification System Based on Proposed Texture Features Extraction Method and Hybrid Ensemble Learning.
Afshin‐Pour et al. A mutual information‐based metric for evaluation of fMRI data‐processing approaches
Bendel et al. A regularized conditional GAN for posterior sampling in image recovery problems
CN105513058B (en) A kind of brain activation area detection method and device
Onal et al. A new representation of fMRI signal by a set of local meshes for brain decoding
CN105894493A (en) FMRI data feature selection method based on stability selection
Song et al. Unsupervised spatiotemporal fMRI data analysis using support vector machines
Bush et al. Improving the precision of fMRI BOLD signal deconvolution with implications for connectivity analysis
Kolahkaj et al. A connectome-based deep learning approach for Early MCI and MCI detection using structural brain networks
Tantisatirapong Texture analysis of multimodal magnetic resonance images in support of diagnostic classification of childhood brain tumours
CN105342615A (en) Brain active region detection method and device
Sundararaj et al. An expert system based on texture features and decision tree classifier for diagnosis of tumor in brain MR images
Li et al. Improving medical/biological data classification performance by wavelet preprocessing
Khalid et al. Adaptive 2DCCA based approach for improving spatial specificity of activation detection in functional MRI
KR102516868B1 (en) 3d convolutional neural network for detection of parkinson&#39;s disease
Mahmoud et al. A hybrid deep contractive autoencoder and restricted boltzmann machine approach to differentiate representation of female brain disorder

Legal Events

Date Code Title Description
C06 Publication
PB01 Publication
C10 Entry into substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant