3d Multi Scale Optic Flow Algorithm

Published Date: 02 Nov 2017

Vimal Chandran

1 Introduction

Cardiac diseases are the leading cause of death and disability in

developed countries and therefore a major public health

challenge. If all forms of cardiac disease were eliminated, the

average life expectancy would increase by 10 years [1]. The

cardiac illness is mainly influenced by the deformation of the

cardiac walls. Assesing the variation in movement of the cardiac

wall may provide a quantitative indication of the health of

cardiac muscle. Montioring and quantification of abnormal cardiac

wall motion is a strong clinical predictor of sudden, arrhythmic

and cardiac death [2].

Estimating the motion field between two consecutive images has

been a heavily investigated field of research for decades. Optic

flow is one of the well known techniques used for motion field

estimation. In this context, the motion field is characterized by

the field of vectors, that shows the displacement of points in

the optic field relative to the observer. Optic flow shows the

velocity field of pixels in the image. In literature [3], there

are several approaches have been proposed for optic flow

estimation and in most of these techniques the assumption is

that, the brightness doesnot change by small displacement and the

motion is estimated by solving the Optic Flow Constrained

Equation.

F_{x}u+F_{y}v+F_{t}=0

where F(x,y,t)\::\mathbb{R}^{2}\times\mathbb{R}^{+}\rightarrow\mathbb{R}

is an image sequence with x, y and t representing the spatial

and temporal coordiantes respectively, F_{x},F_{y},F_{t}

are the

spatiotemporal derivatives, u(x,y,t)\: v(x,y,t),\::\mathbb{R}^{2}\times\mathbb{R}^{+}\rightarrow\mathbb{R}

are the two unkown velocity vectors. It is to be noted that

there are two unknowns with one equation, so a unique solution

cannot be found and this is referred as the "aperture problem" of

the Optic Flow estimation. The unique solution can be obtained if

there are equal number of equation for the number of unknowns.

Inorder to find a probable solution, several methods has been

proposed and they can be grouped into three main categories as

region based, frequency domain based and gradient based

approaches. In the region based approach [3,4,5], velocity is

defined as the shift between regions of subsequent images that

minimizes a sum of squared distance measure for finding the best

match. In frequency domain based approach [6], estimation is

carried out in the frequency domain, where velocity is defined as

energy function and minimization the motion energies is done. In

the gradient based approach [7,8,9], velocity is computed from

the spatiotemporal derivative of the image intensity or filtered

version of the image.

The application of optic flow method to tagged MR images was

introduced by Florack et al [10,11], where a robust differential

technique was developed for scalar and density images. This

generalized framework was adapted to cardiac tagged MR images by

Niessen [12, 13] et al, Suibesuaputra [14]. Van Assen et al [15]

and Florack and Van Asses [16] developed a method based on

multiple independent MR tagging acquition, removing altogether

the aperture problem. Researches in measurement of cardiac motion

are generally carried out in 2D, but for modelling true heart

motion, the 3D model is essential. In 2D optic flow algorithm

captures only the expanison whereas the contraction and rotation

of cardiac tissues are not considered, but with 3D optic flow

algorithm it can be estimated. The latest increase in computation

power made it possible to compute 3D optic flow field from tagged

MRI data.

The 3- dimensional version of the optic flow constrained equation

is given as

F_{x}u+F_{y}v+F_{z}w+F_{t}=0

where F(x,y,z,t)\::\mathbb{R}^{3}\times\mathbb{R}^{+}\rightarrow\mathbb{R}

is an image sequence with x, y, z and t representing the spatial

and temporal coordiantes respectively, F_{x},F_{y},F_{z},F_{t}

are the spatiotemporal derivatives, u(x,y,z,t)\: v(x,y,z,t),\: w(x,y,z,t)\::\mathbb{R}^{3}\times\mathbb{R}^{+}\rightarrow\mathbb{R}

are the three unkown velocity vectors. Similar to 2D, there are

three unknown in this case but only one equation, so the aperture

problem exits. To overcome this problem, Barron [17] extended the

concept developed by Horn and Schrunk, Lucas and Kanade to three

dimension. This method, imposes a constant intensity assumption

which in tagged MR images does not hold due to the tag fading by

T1 relaxation. Inorder to solve the aperature problem and to

model true motion of the heart, we have derived the mutiscale

optic flow technique for 3D case from generalized mutiscale

framework proposed by Florack [11]. Multiple independent MR

tagging acquition are used to solve the aperture problem. In the

experiment, the proposed technique is tested on an artificial 3D

image sequnce, the vectorfield components are calculated and

analysed. The method is further extended to the tagged 3D cardiac

data. In Section 2, generated phantom datasets and real datasets

are discussed. The derived the mutiscale optic flow technique for

3D case is presented in Section 3. Finally in Section 4 and 5 we

describe the experiments, the results and discuss future

directions.

2 Material

2.1 Phantom Data

Image phantom are used to analyse the performance of an

algorithm. In this case, three 3D phantom data is constructed

with a sine function varying in x, y and z direction

respectively. The sequence is constructed by extending in the

temporal direction and with sine function incremented by a

constant t

. The generated phantom sequences are given as:

Fx\left(\sin\left(x+t\right),y,z,t\right),\, Fy\left(x,\sin\left(y+t\right),z,t\right),\, Fz\left(x,y,\sin\left(z+t\right),t\right)\::\mathbb{R}^{3}\times\mathbb{R}^{+}\rightarrow\mathbb{R}

where Fx,\, Fy,\, Fz

represents the 3 dimensional phantom in x,

y, z direction respectively. The generated phantom images has a

resolution of 30\times30

pixels, and an image volume of 30\times30\times30

voxels. There are twenty temporal sequences and the overall

dimension of the generated phanton is 30\times30\times30\times20

representing x\times y\times z\times t

dimensions respectively.

The generated sine phantoms are shown as a varying sine function

in a 3D plot (in Figure 1) and as a 2D sine image (in Figure 2).

[float Figure:

<Graphics file: C:/Users/310079996/Documents/Intern_University/Report_18_2/Phantom Plot.jpg>

Figure 1: 3D Plot of the phantom Fx,\: Fy,\: Fz

at the first

sequence, in x, y, z direction respectively

<Graphics file: C:/Users/310079996/Documents/Intern_University/Report_18_2/Phantom image.jpg>

Figure 2: 2D Image of the phantom Fx,\: Fy,\: Fz

at the first

sequence, in x, y, z direction respectively

]

2.2 Phantom Data with Gaussian Noise

Sensitivity of the algorithm to noise is tested by adding noise

to the phantom sequence, carrying out the same experiments of

noiseless phantom sequence and comparing the results. The

gaussian white noise with zero mean and finite variance is added

to the generated sine phantom sequence. The overall dimension of

the phantom data with noise is 30\times30\times30\times20

representing x\times y\times z\times t

dimensions respectively.

The generated sine phantoms with noise are shown as a varying

sine function with gaussian noise in a 3D plot (in Figure 3) and

as a 2D image (in Figure 4).

[float Figure:

<Graphics file: C:/Users/310079996/Documents/Intern_University/Report_18_2/Noise/Motion/3DPlot.jpg>

Figure 3: 3D Plot of the phantom with gaussian noise at the first

sequence, in x, y, z direction respectively

<Graphics file: C:/Users/310079996/Documents/Intern_University/Report_18_2/Noise/Motion/2Dimage.jpg>

Figure 4: 2D Image of the phantom with gaussian noise at the

first sequence, in x, y, z direction respectively

]

2.3 Real Data

Tagged MR imaging is technique that can be used for quantitative

assessment of myocardial contractile function. The tagging

patterns are inherited in the tissue. They move along the tissue

and enables local motion analysis.The 3 dimensional tagged MR

image volume sequence of a patient heart were acquired using a

HARP technique[2] developed at ETH Zurich, Switzerland. It

consists of 23 sequences with a temporal resolution of 30ms.

There are three different views, one short axis and two long axis

views, which are perpendicular with respect to each other.The

images represents a resolution of 112\times112

pixels, and an

image volume of 112\times112\times112

voxels. The overall

dimension of the generated phanton is 112\times112\times112\times23

representing x\times y\times z\times t

dimensions respectively.

The 3D tagged MR image sequence of a patient heart taken one

short axis and two long axis views is shown (in Figure 5 and 6).

[float Figure:

<Graphics file: C:/Users/310079996/Documents/Intern_University/Report_18_2/CardiacViews.jpg>

Figure 5: From left to right: short axis view containing

horizontal tags, two long axis views containing both the vertical

and horizantal tags and the combination of the image plane.

<Graphics file: C:/Users/310079996/Documents/Intern_University/Report_18_2/Cardiacimage.jpg>

Figure 6: From left to right: 2D Image of the short axis view and

two long axis views of the heart.

]

3 Proposed Method

3.1 3D Multi-Scale Optic Flow Algorithm

The 3D image sequence is a real function F(x,y,z,t)\epsilon\mathbb{R}^{3}\times\mathbb{R}^{+}

with x, y, z and t as spatial coordinates and temporal

coordinates respectively. F(x(t),y(t),z(t))

is the trajectory of

the material point projected to the volume plane. The scale space

representation I(x,y,z,t,\sigma,\tau)\epsilon\mathbb{R}^{3}\times\mathbb{R}^{+}

of a 3D image sequence F(x,y,z,t)

is given as the convolution

between the 3D image sequence and the spatiotemporal gaussian

kernel\phi(x,y,z,t,\sigma,\tau)\epsilon\mathbb{R}^{3}\times\mathbb{R}^{+}

.

I\left(x,y,z,t,\sigma,\tau\right)=\left(F\otimes\phi\right)\left(x,y,z,t,\sigma,\tau\right)=\intop_{-\infty}^{\infty}\intop_{-\infty}^{\infty}\intop_{-\infty}^{\infty}\intop_{-\infty}^{\infty}F(x,y,z,t)\phi(x,y,z,t,\sigma,\tau)dxdydzdt

where \phi(x,y,z,t,\sigma,\tau)=\frac{1}{\left(\sqrt{2\pi\sigma}\right)^{3}\left(\sqrt{2\pi\tau^{2}}\right)}\exp\left(-\frac{\left(x^{2}+y^{2}+z^{2}\right)}{2s^{2}}-\frac{t^{2}}{2\tau^{2}}\right)

and x, y and z are the spatial coordinates, whereas \sigma,\tau\epsilon\mathbb{R}^{+}

denotes the spatial and temporal scales of gaussian kernel

respectively The equation gives the blurred version of the image

and strength depends on the choice of the scale. The brightness

constancy assumption states that the pixel intensities F(x(t),y(t),z(t))

remains constant in time. Inorder to have constant intensities

over time, the Lie-derivative of the observation I with respect

to the vectorfield of the motion is zero. The Lie-derivative

captures the variation of space-time quantities along the

integral flow of some vectorfield. The first order Lie-derivative

of the observation L is considered to have a linear model of

optic flow. The Lie-derivative L_{\vec{v}}

of a function I(\phi)

with respect to the vector field\vec{V}

is defined as L_{\vec{V}}I(\phi)

. The Optic Flow Constraint Equation (OFCE) states that the

luminance does not change when we take the derivative along the

vectorfied of motion.

L_{\vec{V}}I(\phi)\equiv0

For scalar images the Lie-derivative of the observed

spatiotemporal image I is defined as

L_{\vec{V}}I(\phi)\equiv I\left(L_{\vec{V}}^{T}\phi\right)

where L_{\vec{V}}^{T}\phi=-\vec{\nabla}.\left(\phi\vec{V}\right)

[11], so we get

-\intop_{-\infty}^{\infty}\intop_{-\infty}^{\infty}\intop_{-\infty}^{\infty}\intop_{-\infty}^{\infty}F(x',y',z',t')\vec{\nabla}\phi(x-x',y-y',z-z',t-t')\vec{V}(x-x',y-y',z-z',t-t')dx'dy'dz'dt'=0

In the equation 7, the optic flow vectorfield \vec{V}

is

unknown, which is computed by deriving a set of equation. The

optic flow vectorfield with \vec{u}(x,y,z,t)

in x direction, \vec{v}(x,y,z,t)

in y direction and \vec{w}(x,y,z,t)

in z direction is

approximated for the first order which gives

I_{x}u+I_{xt}u_{t}+I_{xx}u_{x}+I_{xy}u_{y}+I_{xz}u_{z}+I_{y}v+I_{yt}v_{t}+I_{xy}v_{x}+I_{yy}v_{y}+I_{yz}v_{z}+I_{z}w+I_{zt}w_{t}+I_{xz}w_{x}+I_{yz}w_{y}+I_{zz}w_{z}=-I_{t}

The first order vectorfield gives an equation with 15 unknown

components of vector field

\left\{ u,u_{t},u_{x},u_{y},u_{z},v,v_{t},v_{x},v_{y},v_{z},w,w_{t},w_{x},w_{y},w_{z}\right\}

. This has to solved but with one equation unique solution is not

possible. However, because the Lie-derivative of the image

vanishes identically so do all their partial derivatives. So the

equations may be added for the vanishing Lie-derivatives with

respect to x, y, z and t since it represents the first order

vectorfield. This gives four extra equations. In total there are

five linearly independent equations

\begin{array}{c}

I_{x}u+I_{xt}u_{t}+I_{xx}u_{x}+I_{xy}u_{y}+I_{xz}u_{z}+I_{y}v+I_{yt}v_{t}+I_{xy}v_{x}+I_{yy}v_{y}+I_{yz}v_{z}+I_{z}w+I_{zt}w_{t}+I_{xz}w_{x}+I_{yz}w_{y}+I_{zz}w_{z}=-I_{t}\\

I_{xx}u+I_{xxt}u_{t}+\left(I_{x}+I_{xxx}\right)u_{x}+I_{xxy}u_{y}+I_{xxz}u_{z}+I_{xy}v+I_{xyt}v_{t}+I_{xxy}v_{x}+I_{xyy}v_{y}+I_{xyz}v_{z}+I_{xz}w+I_{xzt}w_{t}+I_{xxz}w_{x}+I_{xyz}w_{y}+I_{xzz}w_{z}=-I_{xt}\\

I_{xy}u+I_{xyt}u_{t}+I_{xyx}u_{x}+I_{xyy}u_{y}+I_{xyz}u_{z}+I_{yy}v+I_{yyt}v_{t}+I_{xyy}v_{x}+\left(I_{y}+I_{yyy}\right)v_{y}+I_{yyz}v_{z}+I_{yz}w+I_{yzt}w_{t}+I_{xyz}w_{x}+I_{yyz}w_{y}+I_{yzz}w_{z}=-I_{yt}\\

I_{xz}u+I_{xzt}u_{t}+I_{xxz}u_{x}+I_{xyz}u_{y}+I_{xzz}u_{z}+I_{yz}v+I_{yzt}v_{t}+I_{xyz}v_{x}+I_{yyz}v_{y}+I_{yzz}v_{z}+I_{zz}w+I_{zzt}w_{t}+I_{xzz}w_{x}+I_{yzz}w_{y}+\left(I_{z}+I_{zzz}\right)w_{z}=-I_{zt}\\

I_{xt}u+\left(I_{x}+I_{xtt}\right)u_{t}+I_{xxt}u_{x}+I_{xyt}u_{y}+I_{xzt}u_{z}+I_{yt}v+\left(I_{y}+I_{ytt}\right)v_{t}+I_{xyt}v_{x}+I_{yyt}v_{y}+I_{yzt}v_{z}+I_{zt}w+\left(I_{z}+I_{ztt}\right)w_{t}+I_{xzt}w_{x}+I_{yzt}w_{y}+I_{zzt}w_{z}=-I_{tt}

\end{array}

There are five linearly independent equation and fifteen

unknowns, so the unique solution is not possible. This is

refereed as the aperture problem, where the number of vectorfield

to be estimated is larger than the number of linearly independent

equations. In other words, there is not enough information to

recover the optic flow at one point by looking at first order

derivative of image intensity.

Inorder to solve this problem, multiple independent MR tagging

acquition are used, in this case, the same 3D image sequence is

taken in all the three directions x, y and z as F^{X}(x,y,z,t)\epsilon\mathbb{R}^{3}\times\mathbb{R}^{+}

, F^{Y}(x,y,z,t)\epsilon\mathbb{R}^{3}\times\mathbb{R}^{+}

and F^{Z}(x,y,z,t)\epsilon\mathbb{R}^{3}\times\mathbb{R}^{+}

respectively. When their equations are solved with the scale

space representation in directions x, y and z as I^{X}(x,y,z,t,\sigma,\tau)\epsilon\mathbb{R}^{3}\times\mathbb{R}^{+}

,I^{Y}(x,y,z,t,\sigma,\tau)\epsilon\mathbb{R}^{3}\times\mathbb{R}^{+}

and I^{Z}(x,y,z,t,\sigma,\tau)\epsilon\mathbb{R}^{3}\times\mathbb{R}^{+}

respectively. Each direction provide five equation, when added

they provide fifteen equation in total. The first order

vectorfield provides fifteen unknowns and there are fifteen

equation, so unique solution is obtained. This is given in the

form of matrix notation as

AV=-a

A =\left(\begin{array}{ccccccccccccccc}

I_{x}^{X} & I_{y}^{X} & I_{z}^{X} & I_{xt}^{X}\tau^{2} & I_{xx}^{X}\sigma^{2} & I_{xy}^{X}\sigma^{2} & I_{xz}^{X}\sigma^{2} & I_{yt}^{X}\tau^{2} & I_{xy}^{X}\sigma^{2} & I_{yy}^{X}\sigma^{2} & I_{yz}^{X}\sigma^{2} & I_{zt}^{X}\tau^{2} & I_{xz}^{X}\sigma^{2} & I_{yz}^{X}\sigma^{2} & I_{zz}^{X}\sigma^{2}\\

I_{xx}^{X} & I_{xy}^{X} & I_{xz}^{X} & I_{xxt}^{X}\tau^{2} & I_{x}^{X}+I_{xxx}^{X}\sigma^{2} & I_{xxy}^{X}\sigma^{2} & I_{xxz}^{X}\sigma^{2} & I_{xyt}^{X}\tau^{2} & I_{y}^{X}+I_{xxy}^{X}\sigma^{2} & I_{xyy}^{X}\sigma^{2} & I_{xyz}^{X}\sigma^{2} & I_{xzt}^{X}\tau^{2} & I_{z}^{X}+I_{xxz}^{X}\sigma^{2} & I_{xyz}^{X}\sigma^{2} & I_{xzz}^{X}\sigma^{2}\\

I_{xy}^{X} & I_{yy}^{X} & I_{yz}^{X} & I_{xyt}^{X}\tau^{2} & I_{xxy}^{X}\sigma^{2} & I_{x}^{X}+I_{xyy}^{X}\sigma^{2} & I_{xyz}^{X}\sigma^{2} & I_{yyt}^{X}\tau^{2} & I_{xyy}^{X}\sigma^{2} & I_{y}^{X}+I_{yyy}^{X}\sigma^{2} & I_{yyz}^{X}\sigma^{2} & I_{yzt}^{X}\tau^{2} & I_{xyz}^{X}\sigma^{2} & I_{z}^{X}+I_{yyz}^{X}\sigma^{2} & I_{yzz}^{X}\sigma^{2}\\

I_{xz}^{X} & I_{yz}^{X} & I_{zz}^{X} & I_{xzt}^{X}\tau^{2} & I_{xxz}^{X}\sigma^{2} & I_{xyz}^{X}\sigma^{2} & I_{x}^{X}+I_{xzz}^{X}\sigma^{2} & I_{yzt}^{X}\tau^{2} & I_{xyz}^{X}\sigma^{2} & I_{yyz}^{X}\sigma^{2} & I_{y}^{X}+I_{yzz}^{X}\sigma^{2} & I_{zzt}^{X}\tau^{2} & I_{xzz}^{X}\sigma^{2} & I_{yzz}^{X}\sigma^{2} & I_{z}^{X}+I_{zzz}^{X}\sigma^{2}\\

I_{xt}^{X} & I_{yt}^{X} & I_{zt}^{X} & I_{x}^{X}+I_{xtt}^{X}\tau^{2} & I_{xxt}^{X}\sigma^{2} & I_{xyt}^{X}\sigma^{2} & I_{xzt}^{X}\sigma^{2} & I_{y}^{X}+I_{ytt}^{X}\tau^{2} & I_{xyt}^{X}\sigma^{2} & I_{yyt}^{X}\sigma^{2} & I_{yzt}^{X}\sigma^{2} & I_{z}^{X}+I_{ztt}^{X}\tau^{2} & I_{xzt}^{X}\sigma^{2} & I_{yzt}^{X}\sigma^{2} & I_{zzt}^{X}\sigma^{2}\\

I_{x}^{Y} & I_{y}^{Y} & I_{z}^{Y} & I_{xt}^{Y}\tau^{2} & I_{xx}^{Y}\sigma^{2} & I_{xy}^{Y}\sigma^{2} & I_{xz}^{Y}\sigma^{2} & I_{yt}^{Y}\tau^{2} & I_{xy}^{Y}\sigma^{2} & I_{yy}^{Y}\sigma^{2} & I_{yz}^{Y}\sigma^{2} & I_{zt}^{Y}\tau^{2} & I_{xz}^{Y}\sigma^{2} & I_{yz}^{Y}\sigma^{2} & I_{zz}^{Y}\sigma^{2}\\

I_{xx}^{Y} & I_{xy}^{Y} & I_{xz}^{Y} & I_{xxt}^{Y}\tau^{2} & I_{x}^{Y}+I_{xxx}^{Y}\sigma^{2} & I_{xxy}^{Y}\sigma^{2} & I_{xxz}^{Y}\sigma^{2} & I_{xyt}^{Y}\tau^{2} & I_{y}^{Y}+I_{xxy}^{Y}\sigma^{2} & I_{xyy}^{Y}\sigma^{2} & I_{xyz}^{Y}\sigma^{2} & I_{xzt}^{Y}\tau^{2} & I_{z}^{Y}+I_{xxz}^{Y}\sigma^{2} & I_{xyz}^{Y}\sigma^{2} & I_{xzz}^{Y}\sigma^{2}\\

I_{xy}^{Y} & I_{yy}^{Y} & I_{yz}^{Y} & I_{xyt}^{Y}\tau^{2} & I_{xxy}^{Y}\sigma^{2} & I_{x}^{Y}+I_{xyy}^{Y}\sigma^{2} & I_{xyz}^{Y}\sigma^{2} & I_{yyt}^{Y}\tau^{2} & I_{xyy}^{Y}\sigma^{2} & I_{y}^{Y}+I_{yyy}^{Y}\sigma^{2} & I_{yyz}^{Y}\sigma^{2} & I_{yzt}^{Y}\tau^{2} & I_{xyz}^{Y}\sigma^{2} & I_{z}^{Y}+I_{yyz}^{Y}\sigma^{2} & I_{yzz}^{Y}\sigma^{2}\\

I_{xz}^{Y} & I_{yz}^{Y} & I_{zz}^{Y} & I_{xzt}^{Y}\tau^{2} & I_{xxz}^{Y}\sigma^{2} & I_{xyz}^{Y}\sigma^{2} & I_{x}^{Y}+I_{xzz}^{Y}\sigma^{2} & I_{yzt}^{Y}\tau^{2} & I_{xyz}^{Y}\sigma^{2} & I_{yyz}^{Y}\sigma^{2} & I_{y}^{Y}+I_{yzz}^{Y}\sigma^{2} & I_{zzt}^{Y}\tau^{2} & I_{xzz}^{Y}\sigma^{2} & I_{yzz}^{Y}\sigma^{2} & I_{z}^{Y}+I_{zzz}^{Y}\sigma^{2}\\

I_{xt}^{Y} & I_{yt}^{Y} & I_{zt}^{Y} & I_{x}^{Y}+I_{xtt}^{Y}\tau^{2} & I_{xxt}^{Y}\sigma^{2} & I_{xyt}^{Y}\sigma^{2} & I_{xzt}^{Y}\sigma^{2} & I_{y}^{Y}+I_{ytt}^{Y}\tau^{2} & I_{xyt}^{Y}\sigma^{2} & I_{yyt}^{Y}\sigma^{2} & I_{yzt}^{Y}\sigma^{2} & I_{z}^{Y}+I_{ztt}^{Y}\tau^{2} & I_{xzt}^{Y}\sigma^{2} & I_{yzt}^{Y}\sigma^{2} & I_{zzt}^{Y}\sigma^{2}\\

I_{x}^{Z} & I_{y}^{Z} & I_{z}^{Z} & I_{xt}^{Z}\tau^{2} & I_{xx}^{Z}\sigma^{2} & I_{xy}^{Z}\sigma^{2} & I_{xz}^{Z}\sigma^{2} & I_{yt}^{Z}\tau^{2} & I_{xy}^{Z}\sigma^{2} & I_{yy}^{Z}\sigma^{2} & I_{yz}^{Z}\sigma^{2} & I_{zt}^{Z}\tau^{2} & I_{xz}^{Z}\sigma^{2} & I_{yz}^{Z}\sigma^{2} & I_{zz}^{Z}\sigma^{2}\\

I_{xx}^{Z} & I_{xy}^{Z} & I_{xz}^{Z} & I_{xxt}^{Z}\tau^{2} & I_{x}^{Z}+I_{xxx}^{Z}\sigma^{2} & I_{xxy}^{Z}\sigma^{2} & I_{xxz}^{Z}\sigma^{2} & I_{xyt}^{Z}\tau^{2} & I_{y}^{Z}+I_{xxy}^{Z}\sigma^{2} & I_{xyy}^{Z}\sigma^{2} & I_{xyz}^{Z}\sigma^{2} & I_{xzt}^{Z}\tau^{2} & I_{z}^{Z}+I_{xxz}^{Z}\sigma^{2} & I_{xyz}^{Z}\sigma^{2} & I_{xzz}^{Z}\sigma^{2}\\

I_{xy}^{Z} & I_{yy}^{Z} & I_{yz}^{Z} & I_{xyt}^{Z}\tau^{2} & I_{xxy}^{Z}\sigma^{2} & I_{x}^{Z}+I_{xyy}^{Z}\sigma^{2} & I_{xyz}^{Z}\sigma^{2} & I_{yyt}^{Z}\tau^{2} & I_{xyy}^{Z}\sigma^{2} & I_{y}^{Z}+I_{yyy}^{Z}\sigma^{2} & I_{yyz}^{Z}\sigma^{2} & I_{yzt}^{Z}\tau^{2} & I_{xyz}^{Z}\sigma^{2} & I_{z}^{Z}+I_{yyz}^{Z}\sigma^{2} & I_{yzz}^{Z}\sigma^{2}\\

I_{xz}^{Z} & I_{yz}^{Z} & I_{zz}^{Z} & I_{xzt}^{Z}\tau^{2} & I_{xxz}^{Z}\sigma^{2} & I_{xyz}^{Z}\sigma^{2} & I_{x}^{Z}+I_{xzz}^{Z}\sigma^{2} & I_{yzt}^{Z}\tau^{2} & I_{xyz}^{Z}\sigma^{2} & I_{yyz}^{Z}\sigma^{2} & I_{y}^{Z}+I_{yzz}^{Z}\sigma^{2} & I_{zzt}^{Z}\tau^{2} & I_{xzz}^{Z}\sigma^{2} & I_{yzz}^{Z}\sigma^{2} & I_{z}^{Z}+I_{zzz}^{Z}\sigma^{2}\\

I_{xt}^{Z} & I_{yt}^{Z} & I_{zt}^{Z} & I_{x}^{Z}+I_{xtt}^{Z}\tau^{2} & I_{xxt}^{Z}\sigma^{2} & I_{xyt}^{Z}\sigma^{2} & I_{xzt}^{Z}\sigma^{2} & I_{y}^{Z}+I_{ytt}^{Z}\tau^{2} & I_{xyt}^{Z}\sigma^{2} & I_{yyt}^{Z}\sigma^{2} & I_{yzt}^{Z}\sigma^{2} & I_{z}^{Z}+I_{ztt}^{Z}\tau^{2} & I_{xzt}^{Z}\sigma^{2} & I_{yzt}^{Z}\sigma^{2} & I_{zzt}^{Z}\sigma^{2}

\end{array}\right)

V=(u,v,w,u_{t},v_{t},w_{t},u_{x},v_{x},w_{x},u_{y},v_{y},w_{y},u_{z},v_{z},w_{z})^{T}

a=(I_{t}^{X},I_{xt}^{X},I_{yt}^{X},I_{zt}^{X},I_{tt}^{X},I_{t}^{Y},I_{xt}^{Y},I_{yt}^{Y},I_{zt}^{Y},I_{tt}^{Y},I_{t}^{Z},I_{xt}^{Z},I_{yt}^{Z},I_{zt}^{Z},I_{tt}^{Z})^{T}

4 Results

4.1 Phantom Data

4.1.1 Experiment 1

Inorder to assess the accuracy of the proposed 3D mutiscale optic

flow algorithm, an experiment was performed on a more realistic

sine phase grid phantom with non regid motion such as contraction

and expansion. The generated sine phantom has the dimension of 30\times30\times30\times20

representing x\times y\times z\times t

respectively. Isotropic

spatial scale of 2 and temporal scale of 1.5 is used for testing

the algorithm. The fifteen vectorfield components are computed

and their vectorfields taken for the fifth time sequence in the

middle of the image volume are shown in Figure 7.

[float Figure:

<Graphics file: C:/Users/310079996/Documents/Intern_University/Report_18_2/Without/Velocity_components/U.jpg>

\;

<Graphics file: C:/Users/310079996/Documents/Intern_University/Report_18_2/Without/Velocity_components/Ut.jpg>

<Graphics file: C:/Users/310079996/Documents/Intern_University/Report_18_2/Without/Velocity_components/Ux.jpg>

\;

<Graphics file: C:/Users/310079996/Documents/Intern_University/Report_18_2/Without/Velocity_components/Vx.jpg>

<Graphics file: C:/Users/310079996/Documents/Intern_University/Report_18_2/Without/Velocity_components/Wx.jpg>

Figure 7: Fifteen vectorfield components computed by the proposed

algorithm on the generated phantom image.

]

The values of the all vectorfield components at specific location

such as I(1,1,1), I(10,10,10), I(20,20,20) and I(30,30,30) are

analysed. Figure 8 shows the results obtained for different

velocity components at different locations in the 3D image. The x

axis represents the image sequence (time), and yaxis represents

the optic flow velocity.

[float Figure:

<Graphics file: C:/Users/310079996/Documents/Intern_University/Report_18_2/Without/RegionSelection/UVW.jpg>

<Graphics file: C:/Users/310079996/Documents/Intern_University/Report_18_2/Without/RegionSelection/UtVtWt.jpg>

<Graphics file: C:/Users/310079996/Documents/Intern_University/Report_18_2/Without/RegionSelection/UxUyUz.jpg>

<Graphics file: C:/Users/310079996/Documents/Intern_University/Report_18_2/Without/RegionSelection/VxVyVz.jpg>

<Graphics file: C:/Users/310079996/Documents/Intern_University/Report_18_2/Without/RegionSelection/WxWyWz.jpg>

Figure 8: Fifteen vectorfield components at different locations

in the 3D image and in different sequences. (x-axis : image

sequences, y-axis : optic flow velocity)

]

[float Figure:

<Graphics file: C:/Users/310079996/Documents/Intern_University/Report_18_2/Without/Isotropic/UVW.jpg>

<Graphics file: C:/Users/310079996/Documents/Intern_University/Report_18_2/Without/Isotropic/UtVtWt.jpg>

<Graphics file: C:/Users/310079996/Documents/Intern_University/Report_18_2/Without/Isotropic/UxUyUz.jpg>

<Graphics file: C:/Users/310079996/Documents/Intern_University/Report_18_2/Without/Isotropic/VxVyVz.jpg>

<Graphics file: C:/Users/310079996/Documents/Intern_University/Report_18_2/Without/Isotropic/WxWyWz.jpg>

Figure 9: Fifteen vectorfield components at I(20,20,20) in the 3D

image with different isotropical scale setting. (x-axis : image

sequences, y-axis : optic flow velocity)

]

From the Figure 8, it can seen that the optic flow velocities are

continuously changing for different time sequences. The major

variation are found for the extreme voxels located at I(1,1,1)

and I(30,30,30), this is because of the image boundary

conditions. The minor variation are found for the voxels located

at I(10,10,10) and I(20,20,20), which is of interest. The

generated sine phantoms is incremented with a constant (1, in

this case) in the temporal direction for every iteration, so the

expected optic flow velocity is 1 in {u, v, w} vectorfield and

zero for the rest, because their derivatives vanishes. As

expected, the vectorfields {u, v, w}located at I(10,10,10) and

I(20,20,20) gives satisfactory results, except some variations in

the starting and ending part of the sequence due to temporal

boundary condition and the values remains constant from 3-17

sequence.

4.1.2 Experiment 2

The scale selection plays a key role in multiscale optic flow

estimation because the strength of blurring depends on it. In

order to evaluate the proposed method for different scale

setting, an experiments is performed and the results are

analysed. The experiment involves the usage of same sine phantom

sequence as described before, with isotropic scales {1,2,3}and

non isotropic scale{1,1.5,2} setting in both spatial and temporal

coordintates and the results are compared. Figure 9 and 10 shows

the results obtained for different velocity components at the

location I(20,20,20) in the 3D image. The x axis represents the

image sequence (time), and yaxis represents the optic flow

velocity. s(x,y,z,t) represents the spatial scales in x, y, z

direction and the temporal scale, that are used for testing the

algorithm.

[float Figure:

<Graphics file: C:/Users/310079996/Documents/Intern_University/Report_18_2/Without/Different/UVW.jpg>

<Graphics file: C:/Users/310079996/Documents/Intern_University/Report_18_2/Without/Different/UtVtWt.jpg>

<Graphics file: C:/Users/310079996/Documents/Intern_University/Report_18_2/Without/Different/UxUyUz.jpg>

<Graphics file: C:/Users/310079996/Documents/Intern_University/Report_18_2/Without/Different/VxVyVz.jpg>

<Graphics file: C:/Users/310079996/Documents/Intern_University/Report_18_2/Without/Different/WxWyWz.jpg>

Figure 10: Fifteen vectorfield components at I(20,20,20) in the

3D image with different non-isotropical scale setting. (x-axis :

image sequences, y-axis : optic flow velocity)

]

From the Figure 9 and 10, the optic flow vectorfields {u, v, w}

gives the values 1 and rest 0, as expected. There are some

variations in the vectorfields in the begining and the end of the

sequences, which is due to the temporal boundary condition before

a constant value is achieved,. The faster the constant value is

achieved the better the scale would be. Isotropic scale results

looks convincing than the non isotrpic results, which is as

expected because of uniform deformation in all the axis.

4.2 Phantom Data with noise

4.2.1 Experiment 3

The sensitivity of the proposed 3D mutiscale optic flow algorithm

is tested with an experiment where the gaussian white noise with

zero mean and finite variance, is added to the generated image

phantom sequences. The gaussian white noise is uncorrelated

between pixels and it resembles tagged cardiac MR, because

different regions show different deformations.

An experiment to test the sensitivity of the proposed algorithm

is performed on generated image phantom sequences with gaussian

noise with different scale settings as before. The optic flow

vectorfields are computed and the results are shown figure 11 and

12, for a particular location I(20,20,20) in the 3D image. The x

axis represents the image sequence (time), and yaxis represents

the optic flow velocity. s(x,y,z,t) represents the spatial scales

in x, y, z direction and the temporal scale.

[float Figure:

<Graphics file: C:/Users/310079996/Documents/Intern_University/Report_18_2/Noise/Isotropic/uvw.jpg>

<Graphics file: C:/Users/310079996/Documents/Intern_University/Report_18_2/Noise/Isotropic/utvtwt.jpg>

<Graphics file: C:/Users/310079996/Documents/Intern_University/Report_18_2/Noise/Isotropic/UxUyUz.jpg>

<Graphics file: C:/Users/310079996/Documents/Intern_University/Report_18_2/Noise/Isotropic/VxVyVz.jpg>

<Graphics file: C:/Users/310079996/Documents/Intern_University/Report_18_2/Noise/Isotropic/WxWyWz.jpg>

Figure 11: Fifteen vectorfield components at I(20,20,20) in the

3D noisy image, with different isotropical scale setting .

(x-axis : image sequences, y-axis : optic flow velocity)

]

[float Figure:

<Graphics file: C:/Users/310079996/Documents/Intern_University/Report_18_2/Noise/Different/UVW.jpg>

<Graphics file: C:/Users/310079996/Documents/Intern_University/Report_18_2/Noise/Different/UtVtWt.jpg>

<Graphics file: C:/Users/310079996/Documents/Intern_University/Report_18_2/Noise/Different/UxUyUz.jpg>

<Graphics file: C:/Users/310079996/Documents/Intern_University/Report_18_2/Noise/Different/VxVyVz.jpg>

<Graphics file: C:/Users/310079996/Documents/Intern_University/Report_18_2/Noise/Different/WxWyWz.jpg>

Figure 12: Fifteen vectorfield components at I(20,20,20) in the

3D noisy image with different non-isotropical scale setting.

(x-axis : image sequences, y-axis : optic flow velocity)

]

From the Figure 11 and 12, the optic flow vectorfields {u, v, w}

should have the values 1 and rest 0, as before. All the

vectorfield components shows major variation, this is because the

algorithm depends entirely on the voxel brightness and no

critical point features are selected. So even small change in the

voxel intensity value can bring about major change in optic flow

velocity values. So, proper preprocessing has to be carried out.

The blurring process by increasing the scale can give better

prediction, which is as expected. The higher the scale the better

the prominent information is captured and better the optic flow

velocity values. In the experiment, It is evident that the

blurring has reduced the gaussian noise and the constant value

with minor variations is achieved faster for higher scale. The

better results are obtained for the isotrpic scale (3,3,3,3) when

compared to the rest. Figure 13, shows the amount of change

between the optic flow velocities obtained for without noise and

ground truth (1 in this case), for the isotropic scale (3,3,3,3).

[float Figure:

<Graphics file: C:/Users/310079996/Documents/Intern_University/Report_18_2/Noise/Noise.jpg>

Figure 13: Variation in Optic flow velocity calculated for

isotropic scale of (3,3,3,3) at I(20,20,20) of Phantom data with

noise

]

The variations in the optic flow velocities seems to be in the

order of 0.2 for the region of interest excluding the pixels near

the boundary. So the proper scale selection is also depend on the

noise. Similar to the above results the isotropic scaling gives

better results that non isotropic because of uniform deformtion

in all the direction. However in real time appropriate scale at

different locations in the cardiac muscle may be different, since

the heart exhibits different deformations in different regions

such as streching and compression.

4.3 Tagged Cardiac MR Data (Not yet Completed this part)

The MR tags are chessboard like patterns constructed from stripes

move along with the moving tissue. The tag fading is an MR

property and occurs due to finite relaxation time T_{1}

and this

property however does not affect the vanishing image gradient.

The proposed 3D mutiscale optic flow algorithm, is tested on the

tagged cardiac MR data from multiple acquisitions. The dimension

of the data used for testing the algorithm is 30\times30\times30\times20

representing x\times y\times z\times t

respectively. The

isotropic spatial scale of 2 and temporal scale of 1.5 is used.

The results are shown in figure.

The appropriate scale at different locations of the heart is

essential since different regions exhibit different deformations.

The optic flow velocities are retrieved at a certain scale

without considering the size of the basis gaussian function. The

choice of scales greater than zero may provide good results with

respect to noise due to the smoothing related to the increase in

scale. Based on the locatation also, the scale has to be chosen.

5 Discussion

We have proposed a new method for estimating the cardiac motion,

by means of 3D mutiscale optic flow. Initially there are five

independent equation computed for the 3D image sequence in one

direction, and the first order OFCE has fifteen unknown vector

field components. To get a unique solution the same 3D image

sequence is obtained from other two directions, which gave

another ten independent equation. So there is a total of fifteen

independent equations and fifteen unknowns, so the aperatre

problem is removed and unique solution is obtained. In real time,

the multiple independent MR tagging acquitions are used for

solving the aperature problem.

The performance of the algorithm is analysed by performing the

experiment 1. The sine phantom is generated with displacement of

1 unit over frame. The optic flow method, finds the instant

velocity in pixels per frame, which is approximately the same as

the displacement over the frame. The vectorfield gave the correct

results for the region of interest excluding the pixels near the

boundary. The vectorfield derivatives vanishes as expected.

The sensitivity of the algorithm with respect to scale is assesed

in experiment 2. The same phantom setting is used, but with

different isotropic and non-isotropic scale. The scale induces

the blurring effect into the image and strength of blurring

depends on scale selection. The same results in above experiment

is expected, and computed vectorfields with isotropic scale gave

the good results compared to non-isotropic scale. Since the image

has uniform deformation in all the axis, the same the scale in

all directions better the results.

The sensitivity of the algorithm with respect to scale is

examined in experiment 3. The sine phantom with gaussian noise

with displacement of 1 unit over frame is used with different

isotropic and non-isotropic scale. The vectorfield components

computed in this method depends entirely on the voxel intensity

values. Small changes in voxel intensity can alter the

vectorfield. The vectorfield obtained from the experiment showes

large variation when compared to previous results. The scale

selection played a key role in removing the noise information.

Higher scale gave better results than the lower scales.

(*The algorithm is tested on the tagged cardiac MR data sequence.

The tagged sequence is obtained from three simultaneous

acquisition is used. *, Not yet Completed this part)

6 Future Work

The Optic Flow Constrained Equation (OFCE) is approximated to the

first order in the proposed algorithm, but the approximation is

feasible for infinite order. The infinite order of approximation

gives the true instant velocity. The first order approximation

assumes the local pattern as linear, and correct solution is

obtained if the linearity condition is satisfied. If the local

pattern is non-linear, the solution depends on the scale to

linearize the pattern, eventhough the scale change doesn't

influence the derivatives. The zeroth order OFCE gives the

velocity components independent of scale and only the derivatives

depends them. So the optic flow velocity error will be higher

than first order and it decreases with increase in order.

The order of OFCE has to altered based on the motion pattern,

because for linear pattern the first order approximation gives

better results, and for non-linear pattern higher OFCE has to

chosen. The higher order derivative have propagation of noise, to

avoid that better scale selection has to be done. Hence in case

of tagged cardiac MR sequences, different regions exhibit

different motion, so in future the algorithm has to been tuned to

different order and different scale depending on the motion

pattern in a particular region.