## 3 B-Mode Ultrasound Velocimetry Based on Optical Flow
### Implementation of the Proposed Algorithm
#### 3.4.3 Optical Flow Calculation Iteration  
The relationship between the matching error \(\epsilon^{L}\) and the residual value of optical flow \(d^{L}=[d^{L}_{x},d^{L}_{y}]^{T}\) in the \(L\) layer can be expressed as follows:  
\[\epsilon^{L}(d^{L}_{x},d^{L}_{y})=\sum_{x,y}(I(x,y)-J(x+g^{L}_{x}+d^{L}_{x},y+g^ {L}_{y}+d^{L}_{y}))^{2} \tag{12}\]  
\[d^{L}=[d^{L}_{x},d^{L}_{y}]^{T}=arg\,min\{\epsilon^{L}\} \tag{13}\]  
Let \(A(x,y)=i^{l}(x,y)\), \(B(x,y)=j^{l}(x,y)\), then Equation (12) can be expressed as:  
\[\epsilon(d^{L}_{x},d^{L}_{y})=\sum_{x,y}(A(x,y)-B(x+d^{L}_{x},y+d^{L}_{y}))^{2} \tag{14}\]  
Solving the above equation using the L-K algorithm, the first-order partial derivative of \(\epsilon\) with respect to \((d^{L}_{x},d^{L}_{y})\) is zero when \(\epsilon\) is the minimum, and we have:  
\[\frac{\partial\epsilon}{\partial d^{L}}\bigg{|}_{d^{L}=d^{L}_{opt}}=[0,0] \tag{15}\]  
Substituting Equation (14) into Equation (15) yields:  
\[\frac{\partial\epsilon}{\partial d^{L}}=-2\sum_{x,y}(A(x,y)-B(x+d^{L}_{x},y+d^ {L}_{y}))\bigg{[}\frac{\partial B}{\partial x}\quad\frac{\partial B}{\partial y }\bigg{]} \tag{16}\]  
Expanding Equation (16) with the first-order Taylor formula yields:  
\[B(x+d^{L}_{x},y+d^{L}_{y})=B(x,y)+\bigg{[}\frac{\partial B}{\partial x}\quad \frac{\partial B}{\partial y}\bigg{]}\cdot\Big{[}d^{L}_{x},d^{L}_{y}\Big{]}^{ T}+o(x,y) \tag{17}\]  
Substituting Equation (17) into Equation (16) yields:  
\[\frac{\partial\epsilon}{\partial d^{L}}\approx-2\sum_{x,y}(A(x,y)-B(x,y)- \bigg{[}\frac{\partial B}{\partial x}\quad\frac{\partial B}{\partial y}\bigg{]} \cdot d^{L})\bigg{[}\frac{\partial B}{\partial x}\quad\frac{\partial B}{ \partial y}\bigg{]} \tag{18}\]  
The value of \(A(x,y)-B(x,y)\) can be considered as the derivative of the ultrasound image at the point \((0,0)\) with respect to time, denoted as:  
\[\delta I(x,y)=A(x,y)-B(x,y) \tag{19}\]  
The gradient vector of the ultrasound image can be defined as:  
\[\
abla I=\bigg{[} \begin{array}{c}I_{x}\\ I_{y}\end{array} \bigg{]}=\bigg{[}\frac{\partial B}{\partial x}\quad\frac{ \partial B}{\partial y}\bigg{]}^{T} \tag{20}\]  
The gradient vector \([I_{x},I_{y}]\) of the image \(A(x,y)\) can be calculated from the neighborhood of the image at the pixel \(M\). The gradient of the image can be expressed in terms of the central difference as:  
\[ \begin{split} I_{x}(x,y)&=\frac{\partial A(x,y)}{\partial x }=\frac{A(x+1,y)-A(x-1,y)}{2}\\ I_{y}(x,y)&=\frac{\partial A(x,y)}{\partial y}=\frac{A(x,y+1)-A(x,y-1 )}{2}\end{split}  \tag{21}\]Substituting Equations (20) and (21) into Equation (18), we get:  
\[ \begin{array}{rcl}\frac{1}{2}\frac{\partial\epsilon}{\partial d^{L}}&\approx&( \
abla I^{T}\cdot d^{L}-\delta I)\
abla I^{T}\\ &=&\sum\limits_{x,y}\left(\left[\begin{array}{cc}I_{x}^{2}&I_{x}I_{y}\\ I_{x}I_{y}&I_{y}^{2}\end{array}\right]-\left[\begin{array}{cc}\delta I\cdot I _{x}\\ \delta I\cdot I_{y}\end{array}\right]\right)\\ &=&G\cdot d^{L}-b\end{array}  \tag{22}\]  
where,  
\[G=\sum\limits_{x,y}\left[ \begin{array}{cc}I_{x}^{2}&I_{x}I_{y}\\ I_{x}I_{y}&I_{y}^{2}\end{array} \right],\quad b=\left[ \begin{array}{cc}\delta I \cdot I_{x}\\ \delta I\cdot I_{y}\end{array} \right] \tag{23}\]  
When the error \(\epsilon\) achieves its minimum value, the optimal optical flow residual displacement vector \(d_{opt}^{L}\) of the \(L\) layer can be derived as follows:  
\[d_{opt}^{L}=G^{-1}b \tag{24}\]  
To further solve the optical flow residual displacement vector \(d^{L}\) in the \(L\) layer accurately, the Newton-Raphson method can be used to improve the exact solution of the optical flow. Let the optical flow value after \(k-1\) iterations at the \(L\) layer be \(\zeta^{k-1}=[\zeta_{x}^{k-1},\zeta_{y}^{k-1}]^{T}\).