关于超声和光声中的延迟求和算法的详细说明 Detail description of DAS algorithm
- 背景 background
- 存在问题（个人思考）
- 算法说明及代码示例
背景 background
In LED-base Photoacoustic imaging,DAS algorithm was always used for the image reconstruction. This note is about illustrate the details of DAS algorithm.
LED-based Photoacoustic imaging technology often use a linear probe, generally they have 128 elements in this probe,when the exciting light incident to the samples, if there is an absorb in the sample then it will generate the ultrasound signals and all the 128 elements will detected this signal, but due to the different distance between this absorb to every element, it will have a time-delay in each channel, the small distance will get the signal early and the large distance will get it later.
Therefore if you wanna build the image, you need to know the signal amplitude in every pixel, and each pixel was detected by 128-channels, so you need to sum all the 128-channel’s value together to get the right amplitude of every-pixel.
在LED光声中最基本的图像重建算法就是延迟求和算法DAS，通常我们都会使用128陈元的线性探头。当激励光源照射进成像样品之后，如果样品中存在一个吸收体，那么这个吸收体产生的光声信号会被这128个阵元都检测到，但由于每个阵元到这个吸收体的距离不一样，因此导致每个通道监测到的信号有不同的时间延迟，离吸收体近的阵元将会早一点的接收到超声信号，离得远点的就会晚一点接收到这个信号，但是这个信号就是这个吸收体的，只是在接收的时间上有一定的差异，因此为了正确得到这个吸收体产生的光声信号的幅值应该将这128个阵元监测的信号加起来。
存在问题（个人思考）
There is a condition on this algorithm，because if you just based on the distance difference to determine where is the absorb it will make some confuse problems, we all know that the incident light is scattered as a ball in the tissue, consequently，the generated signal can not existed on a plane, it should positioned as a ball, but in generally we assume all the signal was on a plane, then use the DAS algorithm to build the image.
有一个问题就是入射光进入组织以后会散射向各个方向，因此我们接收到的信号不应该只是来自一个平面上的，而应该是一个立体范围内的，但是为了方便计算成像，我们假设用线阵探头进行光声成像的时候接收到的光声信号来自与一个平面即我们通常所说的B-scan,B-plane.而这种线阵探头之所以成像速度快就是其可以一次就完成一个B-scan,不像Photoacoustic microscopy 需要先扫描A-line 然后组成B-plane.
假设所有的吸收体都来自这个B-plane 以后我们就可以根据这个B-plane每个像素距离每个阵元的距离来判断信号的延迟时间，然后根据每个通道接收到信号的时间差异提取出每个通道处信号的具体位置进行叠加，即可以重建出图像。
算法说明及代码示例
DAS延迟和算法，不难从字面理解，计算延迟+数据求和。
（1）延迟时间：通过计算从目标点P到传感器接收单元的距离，再根据传输介质的速度，以及超声传感器阵列的采集频率，求光声信号从目标点到信号接收所需要的时间，由此定位出所收集的数据点在传感器接收数据集中的位置，找到目标点在某个传感器单元所对应的时间的接收信号。
（2）数据求和：在确定了传感器阵列中每个单元所对应的接收时间及对应的信号强度，将数据全部提前进行求和，即为该目标点的信号强度。
根据这个方法，可以求出图像中每个像素所对应的所有阵列信号之和，由此得到整个图像区域的强度信息。
参考：参考来源
这个参考中的图像表述的不够清楚：手画图方便大家理解
- 首先计算延迟时间
要计算延迟时间就得知道每个像素到每个阵元的距离，因为可以由距离得到信号的飞行时间。那么根据两点间的距离公式，第（i，j）个像素到第（N,0）个阵元之间的距离可以表示为（以图象上边的中心和线性探头的中心为原点）N为阵元数量：
我们现在设成像大小为M x M，每个阵元之间的距离为Pitch_size, 每个像素大小为dx,dy
d = s q r t ( （（ i − N / 2 ） ∗ d x − ( C h a n n e l i − N / 2 ) ∗ P i t c h s i z e ） 2 + ( j ∗ d y − 0 ) 2 ) d = sqrt (（（i-N/2）*dx-(Channel_i - N/2)*Pitch_size）^2+(j*dy - 0)^2) d=sqrt(（（i−N/2）∗dx−(Channeli−N/2)∗Pitchsize）2+(j∗dy−0)2)
飞行时间就是用计算出来的距离处理声音传播的速度。
为了确定在每个通道中采集到的到底是一维信号的哪一个点还需要乘以采样频率，采样频率就代表每秒钟产生了多少个数据点，因此有 Vs是声速，fs是探头的采样频率。
t i m e = ( d / V s ) ∗ f s time = (d/Vs)*fs time=(d/Vs)∗fs
%% 参数定义
%%% 定义图像信息
X = 550; Y = 550; % 定义图像大小
dx = 0.1e-3; % 定义图像像素尺寸，单位：m
dy = 0.1e-3; % 定义图像像素尺寸，单位：m
IM = zeros(X,Y); % 定义一个背景图像% 定义换能器信息
N = 128; % 超声换能器阵列单元数量
fs = 80e6; % 数据的采集频率，单位：Hz
dt = 1/fs; % 数据的采集时间分辨率，单位：s
pitch = 0.3e-3; % 阵列单元间隔，单位：m% 传播介质信息
vs = 1.43e3; % 介质声速，单位：m/s%% DAS图像重构
%%for i=1:Xfor j=1:Yfor channel_i = 1:N% 计算目标点到传感器的距离，基于两点距离公式，探头与图像中间对齐（为保证中间较强信号保留）d = sqrt(((i-X/2)*dx-(channel_i-N/2)*pitch)^2+(j*dy-0)^2); % 计算传输需要的时间（延迟时间）time = floor(d/vs*fs)+1; % 在像素位置对所有阵元的数据叠加处理IM(i,j) = IM(i,j) + sensor_data(channel_i,time);endend
end% 根据实际需要进行图像变换
IM=IM'; % 图像反转
相关解释：
（1）有关图像参数设置，可以以现有数据为参考；
（2）探头参数的定义，以实际情况为准；
（3）注意sensor data中通道和时间的行列分布；
（4）可能存在数据叠加时，报错超出数组范围的情况。起初，所采用的解决办法是调整传感器采集频率。传感器实际采集频率为80M Hz，由于超出数据，程序中使用的是75M Hz。后来，改了图像的大小，从600改到了500，感觉比较合适。整体只需要把图像的大小和像素大小控制在比较合适的范围即可。最后，在叠加语句前加了判定语句，如果时间超过时间检测大小，就不执行叠加；时间在采集范围内则叠加。
修改后的叠加程序：
if time<=size(sensor_data,2)% 在像素位置对所有阵元的数据叠加处理IM(i,j) = IM(i,j) + sensor_data(channel_i,time);end
重点参考：https://zhuanlan.zhihu.com/p/606874661
