Construction of arterial geometry from angiography video

Angiograms are routinely used by vascular surgeons to assess artery shape and flow. This research addresses the problem of quantitatively determining the geometry of moderately complex arterial systems from a time sequence of contrast agent concentrations. The contrast agent is assumed to be a passive scalar transported by the blood flow. Measurements of the contrast agent concentration q are available on a number of projection planes { π j , j = 1 , … , n } . The contrast agent satisfies an advection-diffusion transport equation q t + u ⁡ ( x , y , t ) ⁢ q x + v ⁡ ( x , y , t ) ⁢ q y = α ⁢ ∇ 2 q . The frames from an angiography video furnish the concentration field q ⁡ ( x , y , t ) for a number of time levels t ∈ { t 1 , t 2 , … , t n } and for a number of spatial cells ( x i , y j ) that correspond to the recorded digital pixel image. The diffusion coefficient α is assumed to known. The goal is to reconstruct the flow field ( u , v ) and the arterial system geometry from this data. A typical sequence of snapshots is shown below.

Failure of edge detection algorithms

One approach to extracting geometric information from angiograms is to use standard edge detection algorithms. Consider two successive frames as shown below:

The difference between the images

The difference image shows considerable high-frequency noise which can be filtered out (Fourier or z-transform)

The main artery is brought into sharper contrast. Nonetheless, standard edge detection algorithms lead to poor results in identifying artery geometry

Flow based geometry extraction

The problem is that other anatomical features are contributing to the image. These cannot be realistically eliminated in practical situations. However we have not yet used the fact that the artery system should form a closed flow network.

Least-square velocity field, geometry identification

Let q be the dye concentration, integrated along the line of sight. Dye is transported by an advection-diffusion process ∂ q ∂ t + ∂ ( u ⁢ q ) ∂ x + ∂ ( v ⁢ q ) ∂ y = [ ∂ ∂ x ( α ⁢ ∂ q ∂ x ) + ∂ ∂ y ( α ⁢ ∂ q ∂ y ) ]

Discretize the equation. Taylor series expansion in time, central space derivatives:

1. First order q ⁡ ( x , y , t + Δ t ) = q + Δ t 1 ⁢ ∂ q ∂ t + Δ t 2 2 ⁢ ∂ 2 q ∂ t 2 + O ⁢ ( Δ t 3 ) ∂ q ∂ t = − ∂ ( u ⁢ q ) ∂ x − ∂ ( v ⁢ q ) ∂ y + ∂ ∂ x ( α ⁢ ∂ q ∂ x ) + ∂ ∂ y ( α ⁢ ∂ q ∂ y ) Q i ⁣ j n + 1 = Q i ⁣ j n − σ ⁢ [ U i + 1 , j n ⁢ Q i + 1 , j n − U i − 1 , j n ⁢ Q i − 1 , j n + V i , j + 1 n ⁢ Q i , j + 1 n − V i , j − 1 n ⁢ Q i , j − 1 n ] + β ⁢ ( Q i + 1 , j n + Q i − 1 , j n + Q i , j + 1 n + Q i , j − 1 n − 4 ⁢ Q i ⁣ j n ) σ = Δ t 2 ⁢ Δ x , β = α ⁢ Δ t Δ x 2

2. Second order ∂ 2 q ∂ t 2 = − ∂ u ∂ t ⁢ ∂ q ∂ x − ∂ v ∂ t ⁢ ∂ q ∂ y − u ⁢ ∂ ∂ x [ − u ⁢ ∂ q ∂ x − v ⁢ ∂ q ∂ y + α ⁢ ( ∂ 2 q ∂ x 2 + ∂ 2 q ∂ y 2 ) ] − v ⁢ ∂ ∂ y [ − u ⁢ ∂ q ∂ x − v ⁢ ∂ q ∂ y + α ⁢ ( ∂ 2 q ∂ x 2 + ∂ 2 q ∂ y 2 ) ] + α ⁢ ( ∂ 2 ∂ x 2 + ∂ 2 ∂ y 2 ) ⁢ [ − u ⁢ ∂ q ∂ x − v ⁢ ∂ q ∂ y + α ⁢ ( ∂ 2 q ∂ x 2 + ∂ 2 q ∂ y 2 ) ]

First order least squares identification

Let Q i ⁣ j n + 1 be the measured dye concentration, Q ˜ i ⁣ j n + 1 the one predicted by the advection diffusion equation Q ˜ i ⁣ j n + 1 = Q i ⁣ j n − σ ⁢ [ U i + 1 , j n ⁢ Q i + 1 , j n − U i − 1 , j n ⁢ Q i − 1 , j n + V i , j + 1 n ⁢ Q i , j + 1 n − V i , j − 1 n ⁢ Q i , j − 1 n ] + β ⁢ ( Q i + 1 , j n + Q i − 1 , j n + Q i , j + 1 n + Q i , j − 1 n − 4 ⁢ Q i ⁣ j n )

S ⁡ ( U , V , β ) = ∑ i , j ( Q ˜ i ⁣ j n + 1 − Q i ⁣ j n + 1 ) 2 ∂ S ∂ U l ⁣ k = 2 ⁢ ∑ i , j ( Q ˜ i ⁣ j n + 1 − Q i ⁣ j n + 1 ) ⁢ ∂ Q ˜ i ⁣ j n + 1 ∂ U l ⁣ k = 0 ∂ S ∂ V l ⁣ k = 2 ⁢ ∑ i , j ( Q ˜ i ⁣ j n + 1 − Q i ⁣ j n + 1 ) ⁢ ∂ Q ˜ i ⁣ j n + 1 ∂ V l ⁣ k = 0 ∂ S ∂ β = 2 ⁢ ∑ i , j ( Q ˜ i ⁣ j n + 1 − Q i ⁣ j n + 1 ) ⁢ ( Q i + 1 , j n + Q i − 1 , j n + Q i , j + 1 n + Q i , j − 1 n − 4 ⁢ Q i ⁣ j n ) = 0

Solving the above equations gives grid values for ( U l ⁣ k , V l ⁣ k ) . We can now trace out streamlines ⅆ x ⅆ t = u , ⅆ y ⅆ t = v

Flow modeling

Consider each artery formed of many small fibers or flow tubes. Blood flow is assumed to be described by the incompressible Navier-Stokes equations ∇⋅ u ⃗ = 0 u ⃗ t + ( u ⃗ ⋅ ∇ ) ⁢ u ⃗ = − ∇ p + 1 Re ⁢ ∇ 2 u ⃗ + g ⃗ Consider the Frenet triad defined by a flow tube ( e ⃗ T , e ⃗ N , e ⃗ B ) .

In this local coordinate system the velocity is given by u ⃗ ⁡ ( T ) = u ⁡ ( T ) ⁢ e ⃗ T ⁢ ( T ) . The differential nabla operator is ∇ = e ⃗ T L T ⁢ ∂ ∂ T + e ⃗ N L N ⁢ ∂ ∂ N + e ⃗ B L B ∂ ∂ B with ( L T , L N , L B ) the metric coefficients (Lame parameters) of the transformation induced by the flow tube geometry. The continuity equation reduces to ∂ ∂ T ( u T ⁢ L N ⁢ L B ) = 0 The curvilinear parameter T is taken to be the arclength (hence L T = 1 ). We have the Frenet relations d ⁢ e ⃗ T d ⁢ T = κ ⁢ e ⃗ N , d ⁢ e ⃗ N d ⁢ T = − κ ⁢ e ⃗ T + τ ⁢ e ⃗ B , d ⁢ e ⃗ B d ⁢ T = − τ ⁢ e ⃗ N The arc length along the flow tube can be expressed in terms of the flow velocity through T 1 − T 0 = ∫ t 0 t 1 u ⁡ ( T ⁡ ( t ) ) ⁢ ⅆ t or ⅆ T ⅆ t = u . The Frenet relations imply ⅆ e ⃗ T ⅆ t = κ ⁢ u ⁢ e ⃗ N , ⅆ e ⃗ N ⅆ t = ( − κ ⁢ e ⃗ T + τ ⁢ e ⃗ B ) ⁢ u , d ⁢ e ⃗ B d ⁢ T = − τ ⁢ u ⁢ e ⃗ N

Insert the above geometrical relations into the Navier-Stokes equations ( ∂ u ∂ t + u ⁢ ∂ u ∂ T ) ⁢ e ⃗ T + 2 ⁢ u 2 ⁢ κ ⁢ e ⃗ N = ∇ p + ∇ 2 u ⃗ + g ⃗

Introduce a splitting of the pressure p ⁡ ( T , N , B ) = P ⁡ ( T ) + p 1 ⁡ ( T , N , B ) , p 1 ≪ P , p 1 ⁡ ( T , 0 , 0 ) = 0. The streamwise component of the Navier-Stokes equation is ∂ u ∂ t + u ⁢ ∂ u ∂ T = ∂ P ∂ T + 1 Re ⁢ [ ∂ ∂ T ( L N ⁢ L B ⁢ ∂ u ∂ T ) − κ 2 ⁢ L N ⁢ L B ⁢ u ] + g ⃗ ⋅ e ⃗ T