Field_II囊肿模型进行颈动脉血流速度剖面超声仿真--MATLAB

• 实验题目

· 根据颈动脉超声仿真静态模型，建立颈动脉血流速度剖面超声仿真动态模型：

        血流速度要求：中心最大速度设置为1.2 m/s。

        自上管壁至下管壁，血流速度剖面为抛物线型（如图1）：

图1 血流速度图

• · 在该仿真动态模型中，自上管壁至下管壁，将血管腔平均分为30份，每份一个血流速度值，速度计算方法如下： 

  

图2 血流建模示意图

 · 扫描帧频设置为200Hz，生成5张B超图像。

• 方法步骤

· 建立血管的静态模型，这里使用的是静态血管中的建模方法。

http://t.csdn.cn/ZbUly

r1=7/2/1000;%血流
r2=9/2/1000;%血管
yc=0/1000;
zc=9.5/1000+z_start;
inside = ( ((y-yc).^2 + (z-zc).^2) <r1^2);
amp = amp .* (1-inside)+ 0.1*amp .* inside;

       设置血液的超声反射幅度为0.1，让血流的反射强度为周围组织的1/10。

yc=0/1000;
zc=9.5/1000+z_start;
inside = ( ((y-yc).^2 + (z-zc).^2) < r2^2&((y-yc).^2 + (z-zc).^2) > r1^2) ;
amp = amp .* (1-inside) + 5*amp .* inside;

       设置血管壁的超声反射幅度为5，让血管壁的反射强度为周围组织的5倍，血流的50倍。

· 通过题中所给出的公式求解出不同血流层的血流流速。

i=1;
for r=0:7.5/30:3.5
v(i)=1.2*(1-(r/1000)^2/((3.5/1000)^2));
i=i+1;
end
r=0:7.5/30:3.5;

        由于我采用的分层方式是圆环分层，即在血管中心为圆心，不断的增加半径来改变层数的方法，且血流流速呈抛物线分布按血管中心线对称，所以程序只需要循环15次即可得到距离血管中心不同距离的血流流速。

· 添加血流每层的位移。  

for j=1:5
    x1=x;
    for i=1:15
        if i==1
            inside = (((y-yc).^2 + (z-zc).^2) ==r(1)^2);
            s=inside.*v(i)*t*j;
            x1=x1+s;%平移数组
            inside=x1>19/1000;
            x1=x1-38/1000.*inside;
        else
            inside = (((y-yc).^2 + (z-zc).^2) < r(i)^2&((y-yc).^2 + (z-zc).^2) > r(i-1)^2) ;
            s=inside.*v(i)*t*j;
            x1=x1+s;%平移数组
            inside=x1>19/1000;
            x1=x1-38/1000.*inside;
        end
    end
    xresult(:,j)=x1;
end

        在上面的代码片段中，我通过内层for循环把血管分成了中心线对称的30层。然后通过s=v*t来计算出每一层血流的位移，再判断扫描到的血流部分是否是该血流流速对应的部分，如果是则对这部分原有的横坐标x1进行：x1=x1+s的位移处理，如果不是则没有位移。按照实验要求：血流在血管内能循环流动，所以需要对x1的值进行一个移位拼接处理。如果x1在移动后的大小超出了扫描到的血管长度范围[-19mm,19mm]，我们需要把血管的右端超出的部分拼接到血管的左端。内层for循环处理的是血流的每一层，而外层的for循环处理的是我们需要生成的5幅图像的x数据。由于在内层循环中有对每次的位移做时间变化处理，所以在处理完一幅图像后要把x1还原到初始位置，所以在外层循环中有x1=x的赋值操作。最后把每幅图像处理得到的x数据拼接保存到xresult中。

· 调用建模函数建立血流模型。

[phantom_positions, phantom_amplitudes] = cyst_pht(10000);%100000
ph1=phantom_positions(:,[1 6 7]);
ph2=phantom_positions(:,[2 6 7]);
ph3=phantom_positions(:,[3 6 7]);
ph4=phantom_positions(:,[4 6 7]);
ph5=phantom_positions(:,[5 6 7]);

        这个部分是调用建模函数进行血流建模，然后接收建模函数回传的参数。设置散射体的数量为10000，分解建模函数回传的模型位置参数phantom_positions；phantom_positions参数回传的是10000×7的矩阵，前5列是不同时间段（1/200，2/200，…，5/200s）的模型x坐标，后两列是模型的y，z坐标。由于在后续的模型扫描过程中一次只能扫描一组[x,y,z]，所以使用临时变量ph1，…，ph5分解出五组[x,y,x]。

phantom_positions=ph1;
% phantom_positions=ph2;
% phantom_positions=ph3;
% phantom_positions=ph4;
% phantom_positions=ph5;
save pht_data.mat phantom_positions phantom_amplitudes

        这部分代码需要单独运行，每运行一行赋值命令就要运行一次保存命令，然后再运行模型扫描程序扫描出图像数据保存到rf_data文件夹中。更改已有的文件夹名称，新建一个rf_data文件夹命令行执行一条赋值语句和保存语句。重复过程5次生成5幅图像数据。

· 运行成像程序生成图像并保存。

        依次更改建模生成的文件夹名称为rf_data，然后运行成像程序生成血流模型并保存。也可以不更改，但是每生成一幅图像后需要清空这个rf_data文件夹的内容。

注意：一定要保存生成的图像。

• 实验结果

图3 建模生成的坐标

        上图中的数据是调用建模程序后生成的模型坐标数据。前5列是五组x的坐标，可以看出前五列中的x在z=[15mm,24mm]范围内依次递增。但是可以发现在z=15和z=24附近的x没有依次递增，这是因为在靠近血管壁的位置血流速度最慢，按题目中的速度方程可以解出这两个位置附近的血流流速为0，所以x没有依次递增。

以下给出在不添加血流循环流动时的成像：

        从上面的结果可以看到血液在向右流动的过程中在血管的中心线上的速度最快，靠近血管壁的部分流速最慢，血流整体的流动呈抛物线型。所以可以得出结论模型的建立正确。

以下给出在添加血流循环流动时的成像：

 模型建立源码：--cyst_pht.m

%  Create a computer model of a cyst phantom. The phantom contains
%  fiven point targets and 6, 5, 4, 3, 2 mm diameter waterfilled cysts, 
%  and 6, 5, 4, 3, 2 mm diameter high scattering regions. All scatterers 
%  are situated in a box of (x,y,z)=(50,10,60) mm and the box starts 
%  30 mm from the transducer surface.
%
%  Calling: [positions, amp] = cyst_phantom (N);
%
%  Parameters:  N - Number of scatterers in the phantom
%
%  Output:      positions  - Positions of the scatterers.
%               amp        - amplitude of the scatterers.
%
%  Version 2.2, April 2, 1998 by Joergen Arendt Jensen

function [positions, amp] = cyst_phantom (N)

x_size = 38/1000;   %  Width of phantom [mm]
y_size = 0/1000;   %  Transverse width of phantom [mm]
z_size = 20/1000;   %  Height of phantom [mm]
z_start = 10/1000;  %  Start of phantom surface [mm];

%  Create the general scatterers

x = (rand (N,1)-0.5)*x_size;
y = (rand (N,1)-0.5)*y_size;
z = rand (N,1)*z_size + z_start;

%  Generate the amplitudes with a Gaussian distribution

amp=randn(N,1);

%  Make the cyst and set the amplitudes to zero inside

%  6 mm cyst
r1=7/2/1000;%血流
r2=9/2/1000;%血管

yc=0/1000;
zc=9.5/1000+z_start;
inside = ( ((y-yc).^2 + (z-zc).^2) <r1^2);
amp = amp .* (1-inside)+ 0.1*amp .* inside; 

yc=0/1000;
zc=9.5/1000+z_start;  
inside = ( ((y-yc).^2 + (z-zc).^2) < r2^2&((y-yc).^2 + (z-zc).^2) > r1^2) ;
amp = amp .* (1-inside) + 5*amp .* inside;

%求V(r)
i=1;
for r=0:7.5/30:3.5
    v(i)=1.2*(1-(r/1000)^2/((3.5/1000)^2));
    i=i+1;
end
r=0:7.5/30:3.5;
r=r/1000;
t=1/200;

for j=1:5
    x1=x;
    for i=1:15
        if i==1
            inside = (((y-yc).^2 + (z-zc).^2) ==r(1)^2);
            s=inside.*v(i)*t*j;
            x1=x1+s;%平移数组
            inside=x1>19/1000;
            x1=x1-38/1000.*inside;
        else
            inside = (((y-yc).^2 + (z-zc).^2) < r(i)^2&((y-yc).^2 + (z-zc).^2) > r(i-1)^2) ;
            s=inside.*v(i)*t*j;
            x1=x1+s;%平移数组
            inside=x1>19/1000;
            x1=x1-38/1000.*inside;
        end
    end
    xresult(:,j)=x1;
end

positions=[xresult y z];

 调用模型建立部分--mk_pht.m

注意：这个程序（mk_pht）只运行一次，这是因为建模部分已经一次生成了五幅图像的数据不需要运行，还有一个主要原因是建模部分有rand函数（它每次调用得到的随机数是不一样的）会造成结果呈现周围组织、血管壁、血流都在移动的效果。

所以，本程序只运行一次，此后每生成一幅图像并保存后再在命令行窗口运行一条phantom_positions赋值命令和一条保存命令。

%  Make the scatteres for a simulation and store
%  it in a file for later simulation use

%   Joergen Arendt Jensen, Feb. 26, 1997

[phantom_positions, phantom_amplitudes] = cyst_pht(10000);%100000
ph1=phantom_positions(:,[1 6 7]);
ph2=phantom_positions(:,[2 6 7]);
ph3=phantom_positions(:,[3 6 7]);
ph4=phantom_positions(:,[4 6 7]);
ph5=phantom_positions(:,[5 6 7]);

phantom_positions=ph1;
% phantom_positions=ph2;
% phantom_positions=ph3;
% phantom_positions=ph4;
% phantom_positions=ph5;
save pht_data.mat phantom_positions phantom_amplitudes