next up previous

Evaluation of velocity measurements with a
video camera

Stefan Keller
Institute of Low Temperature Science
Hokkaido University, Sapporo, Japan


Date: 8, August 1996


Contents

Introduction

Gravity currents with balls of uniform sizes are carried out in chutes and open, unchanneled slopes. As balls, very light pingpongballs as well as heavy golfballs are used. This text deals with the measurement of the vertical velocity ditribution with a video camera, which is positioned perpendicular above the flow. The theoretical background of such flows will not be discussed here.

Measurement setup

The video camera is set perpendicular to the ground. The measurement principle is as follows: The balls close to the lens of the camera appear bigger than those which are further away. Because all balls in one flow have an uniform diameter, the distance from the camera to a ball can be calculated by measuring the size of a ball in the video picture. In the same way, by measuring the distance the ball moves within a certain time, the velocity of the ball can be calculated. Not only the two dimensional velocity vector, parallel to the ground, can be calculated, but also the component in the z-direction, by taking into account the different ball heights before and after this time interval. This time can either be the time interval of two successive frames, ti (``two picture measurement''), or the time of the shutter speed, ts (``one picture measurement''). In the first case, the ball must be seen in both frames, in the second one, the entire path of the ball, looking as a long, rather faint line of the width of one ball diameter must be recognized and measured.

The height hcam of the camera (resp. lens) over the ground is dependent on the opening angle of the camera lens, on the maximum ball velocity vmax, the maximum flow height hmax and the principle of measurement (one or two picture measurement). In order to get a reasonable resolution of the ball height, the camera should be very close to the flow (see figure 1).


  
Figure 1: Size of a pingpongball in pixels in relation to the distance to the lens. The maximum wide angle position of the lens, without and with additional wide angle adapter is presented. The closer the ball to the lens, the steeper the curve and the better is the height resolution. The size in pixel is referred to a picture size of 512 pixel (long side).
\includegraphics[width=120mm]{BILDER/resolution.eps}

Notation


  
Figure 2: Schematic view of the measurement setup. The notations are explained in the text.
\includegraphics[width=120mm]{BILDER/video.eps}

A drawing of the setup with some basic notations can be seen in figure 2.


  db diameter of ball, 
zb height over ground of ball,
$\vec{v}(z)$ velocity of ball at height z over ground,
u velocity component in flow direction (x),
v velocity component in y - direction,
w vertical velocity component (z - direction),
umax maximum velocity of ball, in flow direction,
hmax maximum flow height,
hlen distance ground - lens,
dlb = hlen - hb, distance lens - ball,
dmin = hlen - hmax, minimum distance lens - maximum flow height,
m, b geometrical characteristics of lenses,
l(dlb) $= p_{max} \cdot m (d_{lb} + b)$, length of representation at distance dlb,
lmin minimum required length of representation at maximum flow height hmax
(dependent of vmax and principle of measurement),
ti frame interval of camera,
ts shutter speed of camera,
pmax number of pixels in the length side of the video frame,
p(dlb) frame fraction (from length side) of ball at distance dlb from lens.

These notations are related to a camera, which is fixed perpendicular to the ground, with the long frame side in the flow direction. However, there might be the situation, that it is not possible to fix it in that way, either that the camera axis is not vertical to the ground or the orientation of the camera frame is not parallel to the flow direction. This has to be taken into account by the evaluation.

Coordinate systems

The coordinate system of the camera frame in pixels (let's say px-py - plane or p-plane with p(px,py) as a point on it) shall be defined with zero in the centre of the frame with the two axis px and py parallel to the edges, with px somewhat in the in main flow direction, i.e., along the long side of the frame (figure 3).


 

 
Figure 3: Coordinate system of the video frame, defined as the px-py - plane (p-plane). It can be regarded as the two-dimensional pixel frame. The origin, 0, is situated in the centre of the frame. The t-plane, where the correction of the lens-distortion is taken into account, is very similar.
\includegraphics[width=50mm]{BILDER/pxpyplane.eps}  

Basically equivalent to the p-plane is the system with the correction of the distortion, the tx-ty - plane or, shorter, the t-plane, with t(tx,ty) as a point on it. It has the same origin and orientation, but the axes have other scalings. The unit shall be called newpixel. A second coordinate system of the camera can be defined as a three-dimensional space (say xc-yc-zc - space or c-space, with c(xc,yc, zc) as a point in the space) with the same orientation as the pixel-coordinates. The units of this system shall be the usual length units, as at the ground coordinate system. The origin 0 shall be in the centre of the video frame, i.e., in the focus of the lens system. In spite of the same orientation as the px-py - plane, the zc - direction to the ground is negatif (figure 4). This is in order to keep the same main orientation as the coordinate system of the measurement setup, with usual units and x-y-z - notation, with the x - axis in the main flow direction and w upwards. The origin of the x-y-z - space (x-space, with g(x,y,z) as a point in it) shall be at the base of the vertical line from the ground plane (x-y) to the focus of the video camera.
 
 
Figure 4: xc-yc-zc - coordinate system (c-space) of the video camera. The origin is defined in the focus of the lens. The orientation is the same is in the pixelframe (figure 3) and in the ground x-y-z - coordinate system, therefor the zc - direction towards the ground is negatif. The units are the common length units.
\includegraphics[width=50mm]{BILDER/xcyczcspace.eps}  

Calibration of video geometry

The calibration includes two parts: The properties of the lens system at the used focal distance with the distortion towards the corners of the frame, and secondly the setup of the camera within the system of the experiment. The first part can be made in the laboratory, and is valid as long as the lens system is unchanged and a constant focal distance is used. The second calibration has to be made at each new experimentel setup, in order to get the orientation of the camera.

Properties of the lens system

For the calibration of the lens system (at a fixed setting) two mappings are necessary:

1.
Distortion fdt: $p \longmapsto t$, $f_{dt}(p(p_x,p_y)) = t(t_x,t_y)\,.$
2.
Focus ffoc: $t \longmapsto c$, $f_{foc}(t(t_x,t_y)) = c(x_c,y_c,z_c)\,.$

Distortion

The calibration of the distortion, fdt, can be made with a grid, which is fixed some distance away from the lens. This grid must be positioned vertical to the axis of the lens. This can be easily made by using a mirror, parallel on the plate with the grid. The camera is now fixed in a position, where the mirror image of the lens is exact in the middle of the whole frame. The orientation and the origin of the grid coordinates must be the same as on the p- and t-plane The pixelcoordinates of different points of the grid can be determined (p-plane, figure 5). Either one can find a closed solution for the mapping fdt, or a grid mapping and interpolation can be used.

The latter may look like the follows: In the p-plane, the coordinates of the points of intersections of the grid lines are measured (such as point p in figure 5). These points can be easily mapped to the t-plane by making the lines parallel and vertical (in the figure, point p is mapped to t). To get the mapping fdt(m(p)) = m(t)of any measurement point m(p) in the p-plane, the grid field, where m(p) is located is searched. The coordinates of the corners of the field are known in the t-plane, therefore, m(t) can be determined by keeping the same location within the the grid field, now in the t-plane.


 

 
Figure 5: Schematical presentation of the distortion. The grid with the parallel and vertical lines is seen on the p-plane as a distorted grid (thick lines). The mapping fdt brings the point p(px,py) in video frame (in pixel-units) to t(tx,ty) (in newpixel-units). In the t-plane (thin lines), all in the original parallel lines must be parallel again. The p- and the t-planes have the same origin and the same orientations.
\includegraphics[width=60mm]{BILDER/distortion.eps}  


 

 
Figure 6: Interpolation from the p-frame to the t-frame. Here, in the p-frame, straight grid lines are assumed. Then, the location of M(p) within the grid is given with certain ratios, which can be used for the calculation of M in the t-frame.
\includegraphics[width=60mm]{BILDER/lineardist.eps}  

The location within the grid can be calculated like the following, assuming straight grid lines (see figure 6; actually, they are bent): The line $\overline{Ma\,Mc}$ goes through M(p) in a way, that $\overline{A\,Ma}/\overline{A\,B} =
\overline{D\,Mc}/\overline{D\,C}$. For the ratio $\overline{Md\,M}/\overline{Md\,Mb}$ on the analogous line between $\overline{A\,B}$ and $\overline{D\,C}$, the same ratio is valid: $\overline{Md\,M}/\overline{Md\,Mb} =
\overline{A\,Ma}/\overline{A\,B}$. Thus, calculating this ratio, it can be used in the t-frame to calculate M(t).

A more accurate method uses a two-dimensional spline interpolation, which is carried out from the rectangular t-frame to the distorted p-frame. With an iteration, f-1dt(m(t)) comes as close to m(p), the starting value, as wished.

Focal length


   
Figure 7: Sketch of the calibration of the coordinate system of the camera (c-space). Big letters refer to calibration values, b denotes the distance from the edge of the lens to the focus. The aim of the calibration is the mapping of a point (such as a) in the t-plane with a given distance d from the lens to a point in the c-space.
\includegraphics[width=105mm]{BILDER/camerasyst.eps}  

The calibration of the focus distance, i.e., of the mapping ffoc: $t \longmapsto c$, can be made with a length scale, fixed at different distances from the lens, vertical to the lens axis in the middle of the frame, in the direction of the tx-axis. Here, the same grid as above can be used, but only the scale in one direction is of importance. At a distance D from the lens, the newpixels T of the length L on the scale are measured. At another distance D2, the newpixels Tm of the length Lm are measured (see figure 7). To get same value T in newpixels at both distances D and D2, it is

\begin{displaymath}L_2 = L_m \,\frac{T}{T_m}\,.
\end{displaymath} (1)

From the proportionality follows

\begin{displaymath}\frac{L}{D + b} = \frac{L_2}{D_2 + b}
\end{displaymath} (2)

and for b

 \begin{displaymath}b = \frac{L \cdot D_2 - L_2 \cdot D}{L_2 - L}\,.
\end{displaymath} (3)

For a point a with a coordinate tx in the t-plane and a given distance d from the lens (i.e., d is the distance to the plane, where the point a is located on, and which is vertical to the axis of the lens), the new x-coordinate in the c-space, xc, is

\begin{displaymath}x_c = \frac{t_x}{T}\,l_x\,,
\end{displaymath} (4)

with lx the equivalent to L and L2, the length proportional to T newpixels. For lx it is, as above,

\begin{displaymath}l_x = \frac{L}{D + b}\cdot (d + b)
\end{displaymath} (5)

and therefore

\begin{displaymath}x_c = \frac{t_x}{T}\,\frac{L}{D + b}\cdot (d + b)\,.
\end{displaymath} (6)

All in all, with

 \begin{displaymath}m = \frac{L}{T\,(D+b)}
\end{displaymath} (7)

and b from equation 3, the lens system with a fixed focus is defined. Because of the radial symmetrie of the lens, m, b and the above calculation for the x-coordinate is valid for the y-coordinate, too. The direct dependances of m and b on the calibration values L, D and T and Lm, Dm and Tm look like the following:
  
b = $\displaystyle \frac{L \, D_m \, T_m - L_m \, D \, T}{L_m \, T - L \, T_m},$ (8)
m = $\displaystyle \frac{L_m \, T - L \, T_m}{T\, T_m \, (D_m - D)}\,.$ (9)

Therefore, the mapping ffoc: $t \longmapsto c$ of a point a(tx, ty) in the t-plane with a distance d from the lens is given with

   
xc = $\displaystyle t_x \cdot m \,(d + b)\,,$ (10)
yc = $\displaystyle t_y \cdot m \,(d + b)\,,$ (11)
zc = $\displaystyle - (d + b)\,,$ (12)

with b and m from equations 9 and 9, respectively.

Two remarks can be added:

1.
The constant m includes both, the lens characterization and the maximum number of pixel pmax in the frame.
2.
With the maximum number of pixels pmax in the length of the frame, the ``view''-length l(d) at the distance d from the lens is given: $l(d) = p_{max} \cdot m (d + b)$.

Setup of the camera in the system of the experiment

This calibration has to be made at each experiment, where the camera is new installed or their position is changed. The procedure is as follows: A scale is fixed at the ground, the direction is the main flow direction and defines the x-axis of the x-space. The scale must go through the point, where the vertical line to the focus of the camera meets the ground plane (V). Again, this can be easily fixed with a mirror, as stated above. Beside of the point V (see figures 8 and 9), two other points on the scale, towards both sides of the frame, shall be fixed while looking through the camera. Their names shall be A and B and the measured distance $\overline{AB}$ is $\Delta D$ ($\Delta D$ remains constant in the x- and c-space). In the pixelframe, the p-plane, the corresponding coordinates pV, pA and pB can be determined.

The aim of this calibration is now to define the mapping from the c-space to the x-space. This includes the determination of the vertical distance hf from the focus to the ground plane and the different rotation angles. In the case that the axis of the lens of the camera is fixed perpendicular to the ground and the x-axis of the p-, resp. t-plane and the c-space is in the direction of the main flow, no rotation to the x-space occurs and only the vertical distance hf has to be determined.

The procedure looks like the following. The points pV, pA and pB are measured in the p-plane. With the mapping fdt: $p \longmapsto t$ we get tV, tA and tB, all points should be on one line yet. The following mapping ffoc: $t \longmapsto c$ for these points can not be calculated, as long as the distances hV, hA and hB are not known (see figure 9). Let's imagine a rotation of the x-space around the y-axis of the c-space, the rotation angle $\epsilon_y$ is given in way, that the new x-coordinate of the point V is 0. This rotation is only related the the x-coordinates and leaves the y-values unchanged.


   
Figure 8: Ground sketch of the calibration of the camera in the system of the experiment. In the p-plane the coordinates of V, A and B, are measured, all of them points in the ground plane, whereas V is the point of intersection of the straight line, which is perpendicular to the ground plane and which goes to the focus of the lens. The direction from A to B determines the main flow direction and the x-axis of the x-space. Since pV is not in the origin of the p-plane, the x-y-plane in the x-space is not parallel to the xc-yc-plane of the c-space. As a consequence, different rotations are needed.
\includegraphics[width=90mm]{BILDER/camposground.eps}  


   
Figure 9: Side view of the calibration of the camera in the system of the experiment. Beside of the different rotation angles, the distance hf from V to the focus of the lens must be determined.
\includegraphics[width=105mm]{BILDER/camposside.eps}  


   
Figure 10: Sketch of the rotations of V to Vy and to Vn. V is the point of intersection between the ground plane and the line from the focus perpendicular to the ground. After the rotations, Vn is situated on the z-axis of the c-space.
\includegraphics[width=83mm]{BILDER/vrotation.eps}  


   
Figure 11: Sketch of the rotations of A to Ay and to An. A lies on the ground plane and has after the rotations the same z-coordinate as Vn, hf.
\includegraphics[width=105mm]{BILDER/protation.eps}  

For the point V, following the mapping from t to c (equation 10), it is

 \begin{displaymath}Vx_c = t_{Vx} \cdot m \cdot h_V\,,
\end{displaymath} (13)

hV, unknown yet, is the distance from V to the plane, which goes through the focus of the camera and which is perpendicular to the axis of the lens (z-axis of the c-space). The rotation angle $\epsilon_y$ is given with (see figure 10)

 \begin{displaymath}\tan \epsilon_y = \frac{Vx_c}{h_{V}} = t_{Vx} \cdot m\,.
\end{displaymath} (14)

After the rotation the x-coordinate of V is zero. For any point $A \neq V$ in the ground plane (see figure 11) it is

\begin{displaymath}Ax_c = t_{Ax} \cdot m \cdot h_{A}\,.
\end{displaymath} (15)

The angle $\alpha_y$ is given as above with

\begin{displaymath}\tan \alpha_y = \frac{Ax_c}{h_{A}} = t_{Ax} \cdot m\,.
\end{displaymath} (16)

After the rotation with the angle $\epsilon_y$ it is

\begin{displaymath}\tan(\alpha_y - \epsilon_y) = \frac{Ax_c^y}{h_{A}^y}\,,
\end{displaymath} (17)

hAy is the new distance to the plane, which goes through the focus and which is perpendicular to the axis of the lens. If dAy is the distance from A and Ay to the y-axis, it is

\begin{displaymath}{d_{Ay}}^2 = {h_A}^2 + (Ax_c)^2 = (h_A^y)^2 + (Ax_c^y)^2\,.
\end{displaymath} (18)

The new x-coordinates in the c-space and the t-plane, respectively, are given with

 \begin{displaymath}Ax_c^y = \tan(\alpha_y - \epsilon_y) \cdot h_A^y
\end{displaymath} (19)

and

\begin{displaymath}t_{Ax}^y = \frac{Ax_c^y}{m \cdot h_p^y} =
\frac{\tan(\alpha_y - \epsilon_y)}{m}\,.
\end{displaymath} (20)

Since $\epsilon_y$ is given with equation 14, tAxy is determined, but not Axcy because of the unknown hAy (equation 19 is not determined).

The next step is the analogous rotation around the x-axis of the c-space. This rotation angle $\epsilon_x$ is given with the condition, that the new y-coordinate of V is 0. This rotation affects only y-values. Again it is

 \begin{displaymath}Vy_c = t_{Vy} \cdot m \cdot h_V\,,
\end{displaymath} (21)

the rotation is not related anymore to the distance hV but

 \begin{displaymath}\tan \epsilon_x = \frac{Vy_c}{d_{Vy}} =
\frac{t_{Vy}\, m\, h_V}{d_{Vy}}\,,
\end{displaymath} (22)

dVy is the distance from V and Vy to the y-axis (analogous to dAy), and it is

\begin{displaymath}{d_{Vy}}^2 = {h_V}^2 + {Vx_c}^2\,.
\end{displaymath} (23)

The relation between hf, the wanted distance from V to the focus, and hV is given with

 \begin{displaymath}{h_f}^2 = {h_V}^2 + {Vx_c}^2 + {Vy_c}^2\,.
\end{displaymath} (24)

Together with the equations 13 and 14, it follows for hV and dVy
hV2 = $\displaystyle {h_f}^2\frac{1}{1 + {t_{Vx}}^2 m^2 + {t_{Vy}}^2 m^2}\,,$ (25)
dVy2 = $\displaystyle {h_f}^2
\frac{1 + {t_{Vx}}^2 m^2}{1 + {t_{Vx}}^2 m^2 + {t_{Vy}}^2 m^2}\,.$ (26)

For the rotation angle $\epsilon_x$ (equation 22) we get therefore

 \begin{displaymath}\tan \epsilon_x =
\frac{t_{Vy} \cdot m}{(1 + {t_{Vx}}^2 m^2)^{\frac{1}{2}}}\,,
\end{displaymath} (27)

i.e., the rotation angle $\epsilon_x$ is independent on the distance hf. For the point A it is

\begin{displaymath}Ay_c = t_{Ay} \cdot m \cdot h_{A}\,,
\end{displaymath} (28)

the angle between Ay = (Axcy,Ayc), the x-axis and the x-z-plane is

\begin{displaymath}\tan \alpha_x = \frac{Ay_c}{h_{A}^y} = t_{Ay} \cdot m \frac{h_A}{h_A^y}\,.
\end{displaymath} (29)

For ha and hAy it is

\begin{displaymath}\frac{h_A}{d_{Ay}} = \sin(\frac{\pi}{2}- \alpha_y) = \cos \alpha_y\,,
\end{displaymath} (30)


\begin{displaymath}\frac{h_{Ay}}{d_{Ay}} = \sin(\frac{\pi}{2}- (\alpha_y - \epsilon_y)
= \cos(\alpha_y - \epsilon_y)\,.
\end{displaymath} (31)

Therefore, the angle $\alpha_x$ is given with

\begin{displaymath}\tan \alpha_x =
t_{Ay} \cdot m \,\frac{\cos \alpha_y}{\cos(\alpha_y - \epsilon_y)}\,.
\end{displaymath} (32)

After the rotation it is
 
$\displaystyle \tan(\alpha_x - \epsilon_x)$ = $\displaystyle \frac{Ay_c^n}{h_{f}}\,,$ (33)
Aycn = $\displaystyle \tan(\alpha_x - \epsilon_x) \cdot h_{f}\,,$ (34)
tAyn = $\displaystyle \frac{Ay_c^n}{m \cdot h_{f}} =
\frac{\tan(\alpha_x - \epsilon_x)}{m}\,.$ (35)

Here, we have hf, since An lies on the ground plane and this is, after these two rotations, parallel to the x-y-plane and goes through Vn, which lies on the z-axis. The relation between hf and hAy, which determines Axcy in equation19 is given with the distance dAyx from Ay to the x-axis:

 \begin{displaymath}\frac{h_f}{\cos (\alpha_x - \epsilon_x)} = d_{A^{y}x} =
\frac{h_A^y}{\cos \alpha_x}\,.
\end{displaymath} (36)

With these two rotations around the y- and the x-axis of the c-space and with the angles $\epsilon_y$ and $\epsilon_x$, given in equations 14 and 27, the ground plane with V is now parallel to the x-y-plane of the c-space. But the distance hf from the ground plane to the focus is not known yet.

hf is related to the measured length $\Delta D = \overline{AB}$. In newpixel units, the length is $\overline{t_A^n t_B^n}$, with equations 10 and 11 it is with tAxn = tAxy and tBxn = tBxy:

 \begin{displaymath}\Delta D = \overline{t_A^n t_B^n} \cdot m \cdot h_f\,,
\end{displaymath} (37)

and with equation 35

\begin{displaymath}\begin{split}
{\overline{t_A^n t_B^n}}^2 & =
(t_{Bx}^n - t...
...\tan (\alpha _{x} - \epsilon_x)}{m}
\biggr) ^2\,.
\end{split}\end{displaymath} (38)

Here, all angles are defined by the newpixel units of V, A and B, and, therefore, the whole system is determined.

The final step is a rotation around the z-axis of the c-space in order to get the angle $\epsilon_z$ of the flow direction (given with $\overrightarrow{AB} $) to the x-axis of the c-space. Because Vxc = Vyc = 0, $\epsilon_z$ is only dependend on A (or B):

\begin{displaymath}\tan \epsilon_z = \frac{Ay_c^n}{Ax_c^n} =
\frac{\tan (\alph...
..._x) \cdot h_f}
{\tan (\alpha_y - \epsilon_y) \cdot h_A^y}\,.
\end{displaymath} (39)

The new coordinates of A are dependend on hAy and hf, with the equation 36 it is finally

\begin{displaymath}\tan \epsilon_z =
\frac{\sin (\alpha_x - \epsilon_x)}
{\tan (\alpha_y - \epsilon_y) \cos \alpha_x}\,.
\end{displaymath} (40)

Therefore, the last rotation around the z-axis is independent on hf, too.

The same calculation can be made for the situation, where V is not situated on the straight line AB. This is even more conveniant for the procedure of the calibration.

Summary of the main equations

Rotation around the y-axis

\begin{displaymath}\tan \epsilon_y = t_{Vx} \cdot m\,.
\end{displaymath}

Rotation around the x-axis

\begin{displaymath}\tan \epsilon_x =
\frac{t_{Vy} \cdot m}{(1 + {t_{Vx}}^2 m^2)^{\frac{1}{2}}}\,.
\end{displaymath}

Relevant angles for a point A on the ground plane

\begin{displaymath}\tan \alpha_y = t_{Ax} \cdot m\,,
\end{displaymath}


\begin{displaymath}\tan \alpha_x =
t_{Ay} \cdot m \,\frac{\cos \alpha_y}{\cos(\alpha_y - \epsilon_y)}\,.
\end{displaymath}

Rotation around the z-axis

\begin{displaymath}\tan \epsilon_z =
\frac{\sin (\alpha_x - \epsilon_x)}
{\tan (\alpha_y - \epsilon_y) \cos \alpha_x}\,.
\end{displaymath}

Distance hf from V to the focus

\begin{displaymath}h_f = \frac{\Delta D}{\,{\overline{t_A^n t_B^n}\,}^*\,}\,,
\end{displaymath}


\begin{displaymath}{({\overline{t_A^n t_B^n}\,}^*)}^2 = \bigl(
\tan (\beta _{y...
...- \epsilon_x) - \tan (\alpha _{x} - \epsilon_x)
\bigr) ^2\,.
\end{displaymath}

Mapping from the system of the camera to the system of the experiment

This mapping consists of the three rotations around the axes of the c-space, the system of the camera, and a shifting along the x-axis of the x-space. Special attention has to be made at the definitions of the angles and the directions of the rotations. It follows, that the necessary rotations around the x-axis and the z-axis are opposite to the definitions of $\epsilon_x$ and $\epsilon_z$. In the following calculation, the correct directions of this angles are already assumed, i.e., the rotations are always in the positive direction.

The point in the c-space shall be $\overrightarrow{c} = (x_c, y_c, z_c)$. The first step is a rotation around the y-axis, with the angle $\epsilon_y$, the value of yc remains untouched. The rotation is given with

\begin{displaymath}\overrightarrow{c} ^y
=
\begin{pmatrix}
\cos \epsilon_y &...
... 0 & \cos \epsilon_y
\end{pmatrix} \,
\overrightarrow{c} \,.
\end{displaymath} (41)

The second rotation is around the x-axis:

\begin{displaymath}\overrightarrow{c} ^x
=
\begin{pmatrix}
1 & 0 & 0 \\
0 ...
... \cos \epsilon_x
\end{pmatrix} \,
\overrightarrow{c} ^y
\,.
\end{displaymath} (42)

The third rotation is around the z-axis, leaving the value of zcx unchanged:

\begin{displaymath}\overrightarrow{c} ^z
=
\begin{pmatrix}
\cos \epsilon_z &...
...0 \\
0 & 0 & 1
\end{pmatrix} \,
\overrightarrow{c} ^x
\,.
\end{displaymath} (43)

Finally, there is the shifting of +hf along the new z-axis

\begin{displaymath}\overrightarrow{x} =
\begin{pmatrix}
0 \\
0 \\
h_f
\end{pmatrix} +
\overrightarrow{c} ^z
\,.
\end{displaymath} (44)

Together, this gives for the new coordinates in the x-space

 \begin{displaymath}\begin{split}
x & = x_c \,(\cos \epsilon_y \, \cos \epsilon_...
...silon_x \, \cos \epsilon_y \, \sin \epsilon_z)\, ,
\end{split}\end{displaymath} (45)


 \begin{displaymath}\begin{split}
y & = x_c \, (\cos \epsilon_y \, \sin \epsilon...
...silon_x \, \cos \epsilon_y \, \cos \epsilon_z) \,,
\end{split}\end{displaymath} (46)


 \begin{displaymath}\begin{split}
z & = - x_c \, \cos \epsilon_x \, \sin \epsilo...
...s \epsilon_x \cos \epsilon_y \\
& \quad + h_f\,.
\end{split}\end{displaymath} (47)

Determination of the location of a ball

In the p-frame, the diameter of the ball is measured using two end points of a diameter in any orientation, D1 and D2. The basic procedure is as the following, but additionally, it needs a correction for the real ball diameter and a correction for the real visibility of the ball.

\begin{displaymath}\begin{split}
D1_p & = (p_{D1x}, p_{D1y})\,,\\
D2_p & = (p_{D2x}, p_{D2y})\,.
\end{split}\end{displaymath} (48)

If the ball is rather a long line because of a high shutter speed (which is the case at the ``one picture measurement''), the shortest diameter at a representative place, i.e., at the beginning or at the end of the line, has to be choosen. These two points will be transformed to the t-frame and we get tD1 and tD2. Now, their middle point tS and the length dbt of $\overline{t_{D1} t_{D2}}$ have to be calculated:

\begin{displaymath}\begin{split}
t_{Sx} & = \frac{t_{D1x} + t_{D2x}}{2}\,,\\
t_{Sy} & = \frac{t_{D1y} + t_{D2y}}{2}\,,
\end{split}\end{displaymath} (49)


\begin{displaymath}d_b^t = \bigl(
(t_{D2x} - t_{D1x})^2 + (t_{D2y} - t_{D1y})^2
\bigr)^{\frac{1}{2}}\,.
\end{displaymath} (50)

Since the real diameter db of the ball is known, the transformation to the c-space follows from equations 10 and 11. The distance from the focus the centre of the ball, hS, is given with

 \begin{displaymath}d_b = d_b^t \cdot m \cdot h_S\,,
\end{displaymath} (51)

and the coordinates of tS in the c-space are

 \begin{displaymath}\begin{split}
Sx_c & = t_{Sx} \, \frac{d_b}{d_b^t}\,,\\
Sy...
...\\
Sz_c & = -h_S = -\frac{d_b}{d_b^t \cdot m}\,.
\end{split}\end{displaymath} (52)

With the equations 45 to 47 we have the mapping from the coordinate system of the camera to the reference system on the ground and therefore, the location of the ball is determined.

Correction of the real ball diameter

Actually, the above mentioned calculation is only valid for a ball diameter, which is measured vertical to the line connecting the centre of the ball with the origin of the p- and t-frame. In all the other situations, the ball diameter appears too big and finally it results a too close distance to the camera. The cause of this error can be understand with figure 12: What should be measured is the real diameter of the ball, $\overline{D1r\,D2r}$. Actually, a projection $\overline{D1\,D2}$ to a ground plane G, which is parallel to the t-frame and which shall be defined to go through the centre of the ball, is measured. In the situation, where the diameter points D1r and D2r have the same distance from the camera, i.e., the diameter is parallel to G (and is on this ground plane, actually), D1 and D1r are congruent.

In the figure (see as well figure 13), beside the mentioned ground plane G, two more planes are used: The D-plane, connecting the origin with the two diameter points D1 and D2, and the V-plane, which goes through the centre of the ball, C, and which is vertical to co, the line from the origin to C. dg is the intersection line between the D and G planes, dv is the intersection line between D and V, and dl is the intersection line between the G and V planes (in figure 12, dl is congruent with the point C).


   
Figure: Representation in the D-plane, which goes through the origin and the two measured diameter pints D1 and D2. The length $\overline{D1\,D2}$ does not represent the real diameter $\overline{D1r\,D2r}$. The necessary relation is given with the angle $\delta $, which is defined with the lines dg and dv. This lines are the intersection lines of two planes with the D-plane: G is the ground plane, which is vertical to the axis of the lens and goes through the centre of the ball, with dg as the intersection line, and V is the plane vertical to co, the line from the origin to the centre of the ball, the intersection line between D and V is dv.
\includegraphics[width=60mm]{BILDER/realdiam1.eps}  


   
Figure: Calculation of the angle $\delta $ in the c-space: For this purpose, the point VP is needed and the angles $\gamma $, $\alpha $ and $\beta $.
\includegraphics[width=100mm]{BILDER/realdiam2s.eps}  

For $\overline{D1r\,D2r}$, it is:

\begin{displaymath}\overline{D1r\,D2r} = \overline{D1\,D2}\cdot \cos(\delta)\,,
\end{displaymath} (53)

where $\delta $ is the angle between dv and dg. $\delta $ can be obtained with the help of co, the line from the centre of the ball to the origin. The angle between dv and co is always $\frac{\pi}{2}$, therefore, only $\gamma $, the angle between co and dg has to be calculated. For this, first the intersection line dg has to be determined. dg is given with D1 and D2 on the G plane, for any (x, y) on dg it is

 \begin{displaymath}y = a \cdot x + c\,,
\end{displaymath} (54)

with
a = $\displaystyle \frac{D2_y - D1_y}{D2_x - D1_x}\,,$ (55)
c = $\displaystyle \frac{D1_y \cdot D2_x - D1_x \cdot D2_y}{D2_x - D1_x}\,.$ (56)

On dg, we need the point VP = (VPx, VPy), where $\overrightarrow{VP\,O} $ is vertical to dg. It is

\begin{displaymath}\overrightarrow{dg} =
\begin{pmatrix}
D2_x - D1_x\\
D2_y - D1_y
\end{pmatrix} \,,
\end{displaymath}


\begin{displaymath}\overrightarrow{VP\,O} =
\begin{pmatrix}
O_x - VP_x\\
O_y - VP_y
\end{pmatrix} \,,
\end{displaymath}

with Ox = Oy = 0 as the origin. The scalar product must be zero:

\begin{displaymath}\langle\overrightarrow{dg} , \overrightarrow{VP\,O}\rangle = 0\,,
\end{displaymath} (57)


\begin{displaymath}(D2_x - D1_x) ( -VP_x) + (D2_y - D1_y) ( -VP_y) = 0\,.
\end{displaymath} (58)

Together with equation 54, this gives for VPx and VPy of the point VP
  
VPx = $\displaystyle \frac{c ( D2_y - D1_y)}{D1_x - D2_x + a(D1_y - D2_y)}$ (59)
VPy = $\displaystyle a\cdot VP_x + c\,.$ (60)

For the angle $\gamma $ it is

\begin{displaymath}\sin \gamma = \frac{\,\overline{VP\,O}\,}{\overline{C\,O}}\,,
\end{displaymath} (61)

with (for the angles, see as well figure 13)
$\displaystyle \overline{C\,O}$ = $\displaystyle \frac{hc}{\cos\alpha}\,,$ (62)
$\displaystyle \overline{VP\,O}$ = $\displaystyle \frac{hc}{\cos\beta}\,.$ (63)

This gives for $\gamma $ and $\delta = \frac{\pi}{2} - \gamma$

\begin{displaymath}\sin\gamma = \frac{\cos\alpha}{\cos\beta} = \cos\delta\,.
\end{displaymath} (64)

I.e., we need to know the angles between the z-axis and CO and $VP\,O$, $\alpha $ and $\beta $, respectively.


   
Figure: Determination of the angle $\phi $ between the z-axis and the line from the origin to a point in the c-space.
\includegraphics[width=60mm]{BILDER/centreangle.eps}  

They can be obtained with the mapping from the t-frame to the c-space. In figure 14 it is for a point P = (Px, Py, Pz): $P_x = t_{Px} \cdot m \cdot hp$ and $P_y = t_{Py} \cdot m \cdot hp$ (m follows from the lens geometry, equation 7).

For $dp = \overline{P\,O}$ it is

 \begin{displaymath}dp^2 = (t_{Px}\cdot m\cdot hp)^2 + (t_{Py}\cdot m\cdot hp)^2 + hp^2\,,
\end{displaymath} (65)

For $\phi $ it follows

\begin{displaymath}\cos\phi = \frac{hp}{dp}
= \frac{1}{\sqrt{m^2({t_{Px}}^2 + {t_{Py}}^2) + 1}\,}\,.
\end{displaymath} (66)

Finally, the angle $\delta $ is determined: With VPx and VPy from equation 59 and 60 and with equation 7, tVx and tVy can be calculated and the angle $\beta $ is obtained. C is given with the diameter points D1 and D2, which gives the necessary angle $\alpha $.

Effective location of the centre of the ball


   
Figure 15: Effective location of the centre of the ball. In case of the diameter is not measured vertical to the line connecting the origin and the ball centre (i.e., where the D- and the V-planes are not congruent), M, the middle of the line from D1 to D2 is not equal to the real centre of the ball.
\includegraphics[width=60mm]{BILDER/effectcentre.eps}  

In the situation, where the diameter of the ball is not measured in the direction, which is vertical to the line co, which goes to the origin, the centre of the ball C = (Cx, Cy) is not simply in the middle of the visible diameter, M = (Mx, My) (figure 15):

\begin{displaymath}M_i = D1_i + \frac{D2_i - D1_i}{2} \neq C_i
\end{displaymath} (67)

(i = x, y). With $x = \overline{C\,D2}$ and r as the ball diameter, it is
$\displaystyle \cos(\delta - \varepsilon)$ = $\displaystyle \frac{r}{x}\,,$ (68)
$\displaystyle \cos(\delta + \varepsilon)$ = $\displaystyle \frac{r}{\,\overline{D1\,D2}-x}\,.$ (69)

For the ratio $x / \overline{D1\,D2}$ it follows

\begin{displaymath}\frac{x}{\,\overline{D1\,D2}\,} =
\frac{\cos(\delta + \vare...
...}
{\cos(\delta - \varepsilon)- \cos(\delta + \varepsilon)}\,.
\end{displaymath} (70)

The angle $\delta $ is the same as in the previous section, and $\varepsilon$ can be obtained with (r = db/2, equation 51, the derivation of dc is like in equation 65)

\begin{displaymath}\begin{split}
\sin \varepsilon &= \frac{r}{dc}\\
&= \frac{...
... \bigl((t_{Cx}\,m)^2 + (t_{Cy}\,m)^2 + 1\bigr)}\,.
\end{split}\end{displaymath} (71)

Here, for tdball, the values from the real ball diameter from the previous section, in newpixel units, have to be used, and since the real C = (Cx, Cy) is not known yet, M can be used. This error can be neglected. Altogether, we get the ratio $x / \overline{D1\,D2}$, which leads to the effective ball centre between D1 and D2.

Effective visibility of the ball


   
Figure 16: Real visibility of the ball diameter: Not the real diameter 2r is seen but a shorter length 2l. The left picture shows a front view of the situation, the right one is a side view.
\includegraphics[width=105mm]{BILDER/visibility.eps}  

From a ball, we do not see the real ball diameter 2r, but the shorter length 2l (figure 16). Instead of

\begin{displaymath}2r = t_{2r} \cdot m \cdot hr
\end{displaymath} (72)

we measure

 \begin{displaymath}2l = t_{2l} \cdot m \cdot hl\,.
\end{displaymath} (73)

The angle $\varepsilon$ is the same as in the previous section:

\begin{displaymath}\sin \varepsilon = \frac{t_{dball} \cdot m}
{2\, \bigl((t_{Cx}\,m)^2 + (t_{Cy}\,m)^2 + 1\bigr)}\,.
\end{displaymath} (74)

Related to $\varepsilon$, there follow different equations:
  
$\displaystyle \tan \varepsilon$ = $\displaystyle \frac{D}{l}\,,$ (75)
$\displaystyle \sin \varepsilon$ = $\displaystyle \frac{r}{dl + D}\,.$ (76)

Together with equation 73, this gives

\begin{displaymath}\frac{t_{2l}\cdot m \cdot hl}{2} =
\frac{r - dl\,\sin\varepsilon}{\sin\varepsilon \, \tan \varepsilon}\,.
\end{displaymath} (77)

The relation from dl to hl has already been used above (equation 65):

\begin{displaymath}hl = \frac{dl}{\sqrt{(t_{Lx}\,m)^2 + (t_{Ly}\,m)^2 + 1}\,}\,
\end{displaymath} (78)

from now on, the root expression shall be abreviated with $\sqrt{\cdots}$. Finally, we get for dl

\begin{displaymath}dl = \frac{2\,r\,\sqrt{\cdots}}
{t_{2l}\cdot m \sin\varepsil...
...tan \varepsilon
+ 2\, \sin\varepsilon \, \sqrt{\cdots}\,}\,.
\end{displaymath} (79)

Together with equation 76, D is obtained and therefore, with dr = dl + D the real (direct, not vertical) distance of the ball is determined. To get hr = hl + D', we obtain for D'

\begin{displaymath}D' = D\,\frac{hl}{dl} = \frac{D}{\sqrt{\cdots}\,}
\end{displaymath} (80)

and there is, finally,

\begin{displaymath}hr = \frac{dl + D}{\sqrt{\cdots}\,}
\end{displaymath} (81)

the real vertical distance from the ball to the lens. The centre of the ball has to be shifted as well, it will be slightly further away from the origin (with i = x, y):

\begin{displaymath}T_{Ci} = T_{Li}\,\frac{dl + D}{dl}\,.
\end{displaymath} (82)

One picture measurement


   
Figure: Representation of the one picture measurement. The shutter speed is set to a low value (such as 1/30s) and the balls appear as a streak. The balls very close to the camera appear big (balls 1 and 2), the longer the streak, the higher the velocity (2 and 3). At each end of the faible streak the diameter, marked with two arrows, has to be measured. To get more accuracy related to the length of the streak, the end points should be measured too, as it is marked at the ball 1.
\includegraphics[width=80mm]{BILDER/onepicture.eps}  

At the one picture measurement, the velocity of a ball is measured within one frame of the video shot. Using a low shutter speed ts, i.e., values of 1/200 to 1/30s, the shape of one ball is not distinct but stretched to a long streak, which is related to the velocity of the ball (figure 17). The longer the streak, the higher the velocity of the ball, the larger the streak, the closer is the ball to the camera. By determining the diameters at the beginning and at the end point of the ball (B1, B2 and E1, E2 at the ball 3 in the figure), the locations SB and SE can be calculated, the three-dimensional velocity $\overrightarrow{v_s} $ of the ball is given with

\begin{displaymath}\overrightarrow{v_s} = \frac{\overrightarrow{S_B S_E} }{t_s}\,.
\end{displaymath} (83)

It is easy to see: to get the velocity, the whole length of the streak has to be recognized. It follows that in the case of a dense flow, only the velocities of the balls in the upper part can be measured.

Since the ball is stretched to a long line, it looses a lot of contrast; it might be very difficult to determine the end points, especially in dense flows, where all the balls have the same color. Using different colors might be very helpful.

Characterization of the setup

In order to get a good resolution of the ball distance, the ball should be very close to the camera (see figure 1). Assuming a maximum flow height over ground hmax and a maximum detectable velocity vmax, it follows for the length of one streak ls

\begin{displaymath}l_s = v_{max} \cdot t_s + d_b\,.
\end{displaymath} (84)

Since the ball should be recognized in one or the following frame, $v_{max} \cdot t_i$ (ti frame interval of the camera) has to be added to ls. The minimum length lmin, which has to be recognized at the flow height of hmax is therefore

\begin{displaymath}l_{min} = v_{max} \,(t_s + t_i) + d_b\,.
\end{displaymath} (85)

The length of representation $l(d_{min}) = p_{max}\cdot m \cdot d_{min}$, dependant on the distance dmin to the focus of the camera and the maximum number of pixels pmax in the video frame in the flow direction, must be equal to lmin. It follows for the minimum required distance from the maximum flow height to the focus of the camera

\begin{displaymath}d_{min} = \frac{v_{max} \,(t_s + t_i) + d_b}{p_{max} \cdot m}\,.
\end{displaymath} (86)

The camera must be fixed at least at the distance hmax + dmin over the ground. Here it is assumed the axis of the camera is perpendicular to the ground.

Two picture measurement


   
Figure: Velocity measurement with two video frames. The time interval from one frame to the other one is ti, the shutter speed is set to a high value (1/10000 - 1/2000s) to get a clear shape of the balls.
\includegraphics[width=53mm]{BILDER/twopicture1.eps} \includegraphics[width=53mm]{BILDER/twopicture2.eps}  

Here the velocity of a ball is measured with the time difference ti from one video shot to the next one (figure 18). The shutter speed must be set very high in order to get a distinct shape of the ball ( $t_s \approx 1/10000 \, \dots \, 1/2000$s). Beside of this, the procedure is analogous to the one picture measurement: The location of the ball is measured in the first frame, with the diameter-points B1 and B2, and in the next frame, E1 and E2.

\begin{displaymath}\overrightarrow{v_i} = \frac{\overrightarrow{S_B S_E} }{t_i}\,.
\end{displaymath} (87)

The balls, or at least, some balls, have to be distinguished clearly from one frame to the other one. This goes better when the balls have different colors or have specific marks. Related to the ``view depth'', the same remarks can be added as above: the denser the flow the less deep the video camera can see.

In addition to the measurement of the velocities, single video frames made with a high shutter speed can be used for the estimation of the ball densities at different heights. Since the balls must not be tracked from one frame to the other one, it is easier to see into the flow and to determine the flow heights of the balls.

Characterization of the setup

The same equations as above are valid, only that the shutter speed ts has to be replaced with ti, since the balls have to be tracked during a longer time. It is:

\begin{displaymath}l_{min} = 2\cdot v_{max} \cdot t_i + d_b\,,
\end{displaymath} (88)


\begin{displaymath}d_{min} = \frac{2 \cdot v_{max} \cdot t_i + d_b}{p_{max} \cdot m}\,.
\end{displaymath} (89)

It follows that the camera must be fixed a bit higher above the flow as at the one picture measurement.

Accuracy of the measurement

The accuracy is mostly dependent on the pixel resolution of the video frame pmax and on the lens characterization. Both informations are included in the constant m. Let's assume an error of $\Delta t$ in pixel or newpixel units (for this calculation the difference can be neglected) at the measurement of any point in the p-plane. The error of the middle point of a ball, given from the end points of the diameter, is $\Delta t/\sqrt{2}$. Related to the diameter, the error is $(2 \, \Delta t)/\sqrt{2}$. From equation 51 follows the error $\Delta h_S$ in the calculation of the distance hS from the focus of the camera:

\begin{displaymath}\Delta h_S = \left \vert \frac{d_b}{(d_b^t \pm (2\,\Delta t/\sqrt{2})) \, m}
- h_S \right \vert \,.
\end{displaymath} (90)

In dependance on the distance of the camera hS it is

\begin{displaymath}\Delta h_S = \left \vert \frac{d_b \cdot h_S}
{d_b \pm (2\,\Delta t/\sqrt{2}) \cdot m \cdot h_S} - h_S \right \vert \,,
\end{displaymath} (91)

the result can be seen in figure 19. The error $\Delta S_c$ in the x- or y-coordinate is (equation 52)

\begin{displaymath}\Delta S_c =
\left \vert
(t_{S} \pm \Delta t/\sqrt{2}) \, ...
...d_b}{d_b^t \pm (2\,\Delta t/\sqrt{2})} - S_c
\right \vert \,,
\end{displaymath} (92)

which leads to

\begin{displaymath}\Delta S_c =
\left \vert
\frac{S_c \cdot d_b \pm (\Delta t...
...\,\Delta t/\sqrt{2}) \cdot m \cdot h_S} - S_c
\right \vert\,.
\end{displaymath} (93)

The error $\Delta S_c$ grows from a minimum at 0 to higher values with increasing Sc (figure 20).


   
Figure: Error $\Delta h_S$ in the calculation of the distance of the ball from the focus. The curves for the errors of $\Delta t = \pm 1, 2, 3$ in the measurement of each pixel coordinate are shown. The parameters are set to common values: db = 38mm and m = 0.00225 (wide angle lens with additional adapter).
\includegraphics[width=72mm]{BILDER/hserror.eps}  


   
Figure: Error $\Delta S_c$ (in mm) in the horizontal component of the coordinate, in dependance on the horizontal coordinate and the vertical distance from the focus. The part to the upper right of the diagonal can not be seen from the camera. It is assumed, that there occurs an error of $\Delta t = \pm 1$ in the measurement of each pixel coordinate, while the other parameters are set to common values: db = 38mm and m = 0.00225.
\includegraphics[width=72mm]{BILDER/hsscerror.eps}  

The error $\Delta v_h$ of the horizontal velocity component is related to the beginning and end point, Bc and Ec, of the vector:

\begin{displaymath}\Delta v_h = \frac{1}{\sqrt{2}}\,
\left \vert \frac{\Delta B_c + \Delta E_c}{t_d}\right \vert\,,
\end{displaymath} (94)

with td equals to ts or ti, related to the measurement method. It can be seen: The higher the time interval td the lower the error. For the error $\Delta v_z$ of the vertical velocity component it is similar

\begin{displaymath}\Delta v_z = \frac{1}{\sqrt{2}}\,
\left\vert\frac{\Delta h_B + \Delta h_E}{t_d}\right\vert\,.
\end{displaymath} (95)

The estimation of the errors in measuring the coordinates shows very clearly the strong increasing of the inaccuracy with increasing distance from the camera and with increasing distance from the axis of the lens. As a conclusion for the measurement it follows that the balls should be as close to the camera as possible, and balls only situated near the axis of the camera should be used for the evaluation of the position and the velocity.

Acknowledgements

This work was carried out in the Institute of Low Temperature Science at Hokkaido University in Sapporo/Japan. It took part on the avalanche dynamics project, which includes large and small scale studies on snow flows and ping pong ball flows in the chute of the Shinjo branch of Snow and Ice studies, NIED, Japan, and on the Miyanomori 70m ski jump field in Sapporo. This project includes as well systematic observations of real avalanches in Kurobe Canyon (Japan). I want to express my thanks to Dr. K. Nishimura, who gratefully supported my stay in this institute. This work was enabled by the Japanese Society for the Promotion of Science (JSPS).

About this document ...

Evaluation of velocity measurements with a
video camera

This document was generated using the LaTeX2HTML translator Version 98.1p1 release (March 2nd, 1998)

Copyright © 1993, 1994, 1995, 1996, 1997, Nikos Drakos, Computer Based Learning Unit, University of Leeds.

The command line arguments were:
latex2html -split 0 veval.tex.

The translation was initiated by Jim McElwaine on 1999-10-27


next up previous
Jim McElwaine
1999-10-27