Control concepts for image-based structure tracking with ultrafast electron beam X-ray tomography

Abstract

In this paper, a novel approach for tracking moving structures in multiphase flows over larger axial ranges is presented, which at the same time allows imaging the tracked structures and their environment. For this purpose, ultrafast electron beam X-ray computed tomography (UFXCT) is being extended by an image-based position control. Application is scanning and tracking of, for example, bubbles, particles, waves and other features of multiphase flows within vessels and pipes. Therefore, the scanner has to be automatically traversed with the moving structure basing on real-time scanning, image reconstruction and image data processing. In this paper, requirements and different strategies for reliable object tracking in dual image plane imaging mode are discussed. Promising tracking strategies have been numerically implemented and evaluated.

Keywords

Ultrafast X-ray computed tomography structure tracking multiphase flow model predictive control

Introduction

Fast imaging techniques are essential experimental tools in fluid dynamics research, particularly for the investigation of multiphase flows. Such flows are to be found in, for example, chemical reactors, thermal power plants, oil and gas processing, refrigeration systems, fluid transport systems and others. Multiphase flows are inherently complex due to the co-existence of multiple phases with different physical properties and highly deformable interfaces; for example, in gas-liquid and liquid-liquid flows. The modelling and numerical simulation of multiphase flows with computational fluid dynamics (Brennen, 2005: 20 ff.) require validation data from experiments. For such measurement techniques with high temporal and spatial resolution are needed. Particularly, non-invasive imaging is of great attractiveness (Reinecke et al., 1998). However, today there is no multiphase flow imaging technique available, which gives a full three-dimensional (3D) view of opaque multiphase flows at high speed and resolution. As an example, we address a simple bubbly flow in a pipe or vessel, that is, a coexistent flow of continuous liquid and disperse bubble in a flow channel. As such, it would be ideal to “see” all bubbles in a given volume, that is, pipe or column section, at a spatial and temporal resolution that discloses the dynamics of the bubbles themselves, the dynamics of their interfaces and their geometric relations to other bubbles. Of course, it would also be of interest to know the local continuous liquid velocities. As this is all not possible today, an alternative would be to have a fast cross-sectional, that is, two-dimensional (2D) imaging, with a capability to track a moving structure, like a bubble. This would at least allow to study a single bubble and its environment dynamically. In the following, we will describe our approach to solve this task by tracking structures in a flow with an ultrafast X-ray tomography device. Before that, we will briefly review dynamic imaging techniques for multiphase flows, to provide a good start into the subject.

Today, there are a few fast flow measurement techniques available for scientists and engineers. Most prominent are high-speed cameras. However, they are seldom applicable for multiphase flows as such are inherently opaque. The same holds for ultrasound-based imaging techniques. Techniques, which are in principle suitable for multiphase flows, are electrical tomography, radiation-based imaging techniques and magnetic resonance imaging. Hence, we will first briefly review these techniques with a focus on their current capability to deliver a full picture of, for example, a bubble flow in a pipe or column.

Electrical tomography, with its variants electrical capacitance tomography (ECT), electrical resistance tomography (ERT) and electrical impedance tomography (EIT), is able to recover structures with different electrical properties in multiphase flows. The techniques can be made quite fast and have a certain 3D imaging capability (Xie et al., 1995). However, their spatial resolution is very limited, which hampers tracking of dedicated structures.

Transmission imaging techniques based on radiation can be classified into radiographic and tomographic imaging. Fast radiography is known for X-rays (Boden et al., 2014) and neutrons (Kaestner et al., 2011a, 2011b). X-ray radiography may be used and has been implemented in the form of X-ray stereoscopic particle tracking (see below). Conventional X-ray tomography, that is being used in medicine and for non-destructive testing, is generally too slow as it requires mechanical movement of system components (Kak and Slaney, 1988). An exception is ultrafast X-ray tomography, which we address in this study (Fischer et al., 2008).

Magnetic resonance imaging (MRI) provides high spatial resolution and 3D imaging capacity. However, the fastest MRI scanners achieve up to 50 fps (frames per second) (Müller et al., 2007; Tayler et al., 2012b), which is still quite slow compared with the time scales of multiphase flows. Recent applications of ultrashort echo time MRI have attempted to track moving structures in fluidized beds (Fabich et al., 2017) and bubbly flows (Tayler et al., 2012a) at 40 fps. A main advantage of MRI is the simultaneous measurement of liquid-phase velocities that cannot be achieved by radiation-based tomography methods. However, the relatively small volume of investigation limits the application to snapshots of bubble swarms rather than visualizing their paths and dynamics over a larger range.

Besides imaging techniques, there exist also pure structure tracking techniques for multiphase flow. Such are CARPT (Degaleesan et al., 2001), PEPT (Cole et al., 2010), and X-ray stereoscopy (Nadeem and Heindel, 2018). In PEPT and CARPT, radioactively marked tracer particles are tracked over a longer period. As they are pure tracking techniques the surroundings of the tracked structures (particles) are not disclosed. X-ray stereoscopy obtains the position of high-contrast particles from triangulation. To some degree surrounding structures are disclosed in the associated X-ray images. However, as these images provide no good depth information, the 3D imaging capacity is rather low (Nadeem and Heindel, 2018).

Ultrafast X-ray tomography is a relatively new fast imaging technique that comes very close to the given requirements. It bases on X-rays, provides high frame rate of up to 8000 fps, 1 mm spatial resolution and dual plane imaging (Bieberle et al., 2012). As this technique yet has only 2D imaging capability, though in two planes simultaneously, structure tracking out of the imaging planes is not possible. To qualify this technique for structure tracking, a fast positioning system for the scanner has to be implemented and a fast data analysis strategy for controlling the system has to be developed. The latter development has been started already in our group with the transfer of the data processing onto graphic processing units (GPU) (Bieberle et al., 2017; Frust et al., 2017). This now gives short enough latency times for real-time axial scanner position control.

In the scope of the control concept the two imaging planes of the X-ray scanner act as a pair of binary structure sensors. Many control schemes for binary sensors or binary sensor networks (BSN) are known from the literature. However, their application is limited to either unidirectional continuous processes, for example, liquid level control (Wang et al., 2003), or multi-dimensional positioning problems for zero-dimensional objects, for example, tracking of persons across various rooms, (Bai et al., 2015). Thus, Liu et al. (2013) proposed a strategy for extraction of velocity and size of a circular object moving at constant speed using a set of transmitter-receiver style binary sensors. However, the extracted information is used neither to control the process nor to track it actively by repositioning the sensors.

There are also various strategies for image-based control, for example, based on high-speed camera footage (Hutchinson et al., 1996). They are, however, only applicable for controlling processes within the imaging region. Combining the available concepts, we propose a control strategy for image-based control of processes occurring normal to the imaging plane(s). Here, instead of a BSN, only a pair of controllably traversable binary sensors is used.

To track a moving structure, an effective position control of the ultrafast X-ray CT scanner must be provided, employing a low-latency trajectory generation. For that a so-called model predictive control approach is most suited in which the target’s position changes are suitably modelled. Furthermore, constraints regarding the maximal applicable velocity and acceleration need to be included. This requires iterative computation techniques that are critical in terms of real-time operation (Kim et al., 2007; Neunert et al., 2016). Current analytic approaches do not include such a movement model (Haschke et al., 2008; Ruppel et al., 2011). However, for tracking of a temporarily non-visible moving target, such as a bubble moving between imaging planes, this is indispensable. Thus, a time-optimal, analytic trajectory generation is proposed that complies with maximum velocity and acceleration constraints. A double-setpoint controller is employed based on the periodic modelling of the structure’s movement with a parabolic position profile. Similar to the receding horizon approach, the model parameters are updated in every time-step and only the optimal control output for the next time step is used. Different from classic receding horizon control, the analytic solution allows modelling the movement up to the current control target instead of exploring only a fixed prediction horizon.

Materials

The ultrafast electron beam X-ray CT scanner

As shown in Figure 1a, the ultrafast electron beam X-ray CT scanners use an electron beam that is focussed and deflected onto a circular tungsten target to generate a rapidly rotating X-ray focal spot along two staggered paths without any mechanical movement of components. Data acquisition is performed by a radiation detector that comprises two distinct rings of CdTe detector pixels arranged concentric to the X-ray target. The CT scanner under consideration offers a circular imaging area with a diameter of $190 mm$ and two distinct imaging planes with a geometric distance $Δ h_{c} = 13 mm$ (Figure 1b). The alternating scanning in two planes with a small axial offset allows the determination of structure velocities by time-of-flight methods (Barthel et al., 2015). Currently, each detector ring comprises 432 seamlessly arranged pixels, whose analogue signals are digitized with a sampling frequency of $f_{samp} = 10^{6}$ samples/s. Cross-sectional images are acquired alternately in both imaging planes (see Figure 1c) with a maximal imaging rate of up to $f_{im} = 8000 fps$ .

Figure 1.

(a) Principle sketch of the ultrafast electron beam X-ray CT, (b) vertical view of the object area. Horizontal lines indicate the imaging planes and (c) principal order of alternating cross-sectional images per plane over time. Latency is sketched for the real-time image reconstruction.

Image reconstruction is currently performed offline. That is, the digitized signals of the radiation detectors are temporarily stored in the random access memory of the detector electronics before they are being transferred to the host computer after completion of a scan (Bieberle et al., 2017). There, image reconstruction and post-processing is done at a later time.

For structure tracking, real-time image reconstruction as well as control strategies to traverse the CT scanner are needed. Data processing pipeline systems (Frust et al., 2017; Kopmann et al., 2016) seem to be proper tools for that. For that we evaluated three different traverse control concepts, as explained in the following.

Tracking strategy

A moving structure, in this case represented by a rising gas bubble, is fully described by the following set of one-dimensional parameters (Figure 2a): the current position $s_{s}$ of the structure “s”, its velocity $v_{s}$ , its acceleration $a_{s}$ and its axial length $l_{s}$ . Furthermore, the CT scanner “c” is also fully described by its current position of the upper imaging plane (called scanning position) $s_{c}$ , velocity $v_{c}$ and acceleration $a_{c}$ . As constraints we have a maximal absolute velocity $v_{c}^{\max}$ , a maximal absolute acceleration $a_{c}^{\max}$ and position limits $s_{c}^{\min}$ and $s_{c}^{\max}$ . Moreover, we have the constant geometric distance $Δ h_{c}$ between both imaging planes $h_{0}$ and $h_{1}$ and the time interval $Δ t_{cycle} = 2 / f_{im}$ for the acquisition of an image-pair. The latter serves as a master clock for parameter extraction and storage. With each newly acquired pair of images $n \in [1, . ., N_{\max}]$ the system clock is updated to $t_{n} = t_{n - 1} + Δ t_{cycle}$ .

Figure 2.

(a) Description of important system parameters. (b) Qualitative scanner movement for the different control strategies.

To track the moving structure, the scanning position $s_{c}$ has to be adapted by controlling the scanner’s movement based on these parameters. To determine this target scanning position $s_{c}^{target}$ , three strategies have been considered. For this, an additional offset parameter $p < 1$ is introduced which is used to calculate $s_{c}^{target}$ based on the current structure length $l_{s}$ using $s_{c}^{target} = s_{s} + s_{i}^{off}$ for $i \in {1, 2}$ with $s_{1}^{off} = - p \cdot l_{s}$ and $s_{2}^{off} = - (Δ h_{c} + (1 - p) \cdot l_{s})$ as described below. A value $0 < p < 1$ leads to a target scanning position within the tracked structure at portion $p$ of its length $l_{s}$ , that is, the structure remains visible, whereas $p < 0$ leads to deliberate overshooting by $| p | \cdot l_{s}$ , that is, temporary non-visibility of the structure. The default values $p$ for each of the considered tracking strategies were determined empirically (see Table 1) and may be adapted within the respective given ranges (see below) for each of the following strategies.

Table 1.

Parameter space for the numeric simulations of the behaviour of the complete system, that is, CT scanner and control unit, as depicted in Figure 9.

Parameter	Symbol	Values
Input delay	$Δ t_{in}$	$[0, 5, 10] ms$
Output delay	$Δ t_{out}$	$[0, 1, 2] ms$
CT scanner position uncertainty	$Δ s_{c}^{meas}$	$[0, 0.1, 0.5] mm / sample$
Control strategy	−	$SIPT (p = 0.5)$ $DIPT (p = 0.25)$ $DIRPT (p = - 0.25)$
Movement profiles	−	(1) $v_{s} (t_{0}) = 50 mm \cdot s^{- 1}, a_{s} = 0$ (2) $v_{s} (t_{0}) = 300 mm \cdot s^{- 1}, a_{s} = 0$ (3) $v_{s} (t_{0}) = 50 mm \cdot s^{- 1}, a_{s} = 2 mm \cdot s^{- 2}$ (4) $v_{s} (t_{0}) = 300 mm \cdot s^{- 1}, a_{s} = 10 mm \cdot s^{- 2}$ (5) $v_{s} (t_{0}) = 50 mm \cdot s^{- 1}, a_{s} = - 0.2 mm \cdot s^{- 2}$ (6) $v_{s} (t_{0}) = 300 mm \cdot s^{- 1}, a_{s} = - 5 mm \cdot s^{- 2}$ (7) $up / down with \bar{\| v_{s} \|} = 200 mm \cdot s^{- 1}$
Length profiles	−	(1) $l_{s} = 8 mm$ (2) $l_{s} = 18 mm$ (3) $l_{s} = 8 mm + 4 mm \cdot s_{s} \cdot (2500 mm)^{- 1}$ (4) $l_{s} = 18 mm + 9 mm \cdot s_{s} \cdot (2500 mm)^{- 1}$ (5) $l_{s} = 8 mm - 4 mm \cdot s_{s} \cdot (2500 mm)^{- 1}$ (6) $l_{s} = 18 mm - 9 mm \cdot s_{s} \cdot (2500 mm)^{- 1}$ (7) $l_{s} = 8 mm + 2 mm \cdot \sin (t \cdot 2 s^{- 1} \cdot 2 π)$ (8) $l_{s} = 18 mm + 4.5 mm \cdot \sin (t \cdot 2 s^{- 1} \cdot 2 π)$

Single image plane tracking (SIPT, $0 < p < 1$ ): In the first strategy, the CT scanner’s position is controlled such that the structure remains visible in the lower image plane $h_{0}$ . Therefore, the parameters $s_{s}, l_{s}$ and $v_{s}$ are initially determined while the structure passes both imaging planes first. Afterwards, the CT scanner drives to its target scanning position $s_{c}^{target}$ at a defined portion of the structure’s length $l_{s}$ within the image plane $h_{0}$ with a final velocity of $v_{c} = v_{s}$ , as shown in Figure 2b (SIPT). The target position is therefore defined to be $s_{c}^{target} = s_{s} - p \cdot l_{s}$ . In case the velocity of the structure $v_{s}$ does not change over the section to be investigated, the scanner has not to be accelerated anymore. Anyway, in the more likely case, where $v_{s} \neq const .$ , the velocity of the scanner $v_{c}$ has to be iteratively adapted correspondingly as described in chapter 2.3. With the structure visible most of the time in one imaging plane, this strategy is expected to provide detailed trajectory information. However, size and shape changes of the structure cannot be detected.

Dual image plane tracking (DIPT, $0 < p < 1$ ): As in SIPT, the structure is initially completely scanned in both imaging planes to determine its current size and velocity. Then the CT scanner moves upward, thereby overtaking the structure with the upper imaging plane and partially with the lower imaging plane until the lower image plane is a distance $p \cdot l_{s}$ above the bottom of the structure (Figure 2b-DIPT). There, the CT scanner accelerates with $- a_{c}^{\max}$ such that the upper image plane is positioned at distance $p \cdot l_{s}$ below the top of the rising structure. For very large or slowly rising structures, this leads to intentional downwards movement of the scanner. Thereby, the number of sampling points $t_{k}$ is maximized (see section 2.3). This procedure is then repeated. With DIPT the structure is always in the axial region between the upper and lower image plane and visible for almost all the time. By letting the structure frequently pass through both imaging planes, its shape and size can be continuously measured. However, this strategy is still somewhat sensitive to fast changes of the structure’s geometry that may result from, for example, flow disturbances.

Dual image plane re-tracking (DIPRT, $p < 0$ ): In the third strategy, the structure is repeatedly scanned completely by the CT scanner in both image planes. As described for the DIPT strategy, the target scanning position again alternates between the two target offsets $s_{1}^{off}$ and $s_{2}^{off}$ but this time with the offset parameter $p < 0$ . Thereby, the scanner drives above the structure and scans again the structure’s parameters. Afterwards, the scanner accelerates with $- a_{c}^{\max}$ to let the structure pass completely through both imaging planes once more and repeats the cycle. Employing this strategy, the structure’s shape and size changes can be characterised. However, information about its trajectory is limited to the intermediate scans. Further, the structure leaves the scanning planes multiple times and needs to be rediscovered.

Parameter extraction of the moving structure

No matter which tracking strategy will be applied, parameters of the investigated object have to be extracted initially from the image pair data. However, as recognizable from Figure 1c, in a single image plane, a small and slowly rising structure leads to the same image sequence as a tall and fast rising structure. Thus, for the extraction of the structure velocity, the images from the second imaging plane $h_{1}$ need to be employed to track the boundaries of the structure whilst passing through both imaging planes (see Figure 3a).

Figure 3.

Determination of structure velocity and axial dimension using dual-image mode of the CT scanner.

The time interval $Δ t_{1} = t_{2} - t_{1}$ (see Figure 3a) is initially determined as the time between the structure’s front entering the lower plane $h_{0}$ and upper plane $h_{1}$ . Afterwards, the velocity of the structure’s boundary can be calculated by

v_{s} = \frac{Δ h_{c} + s_{c} (t_{2}) - s_{c} (t_{1})}{Δ t_{1}}

(1)

Secondly, the time interval $Δ t_{2} = t_{3} - t_{1}$ between the entering of the structure’s front at the lower image plane $h_{0}$ and the structure’s back at the lower image plane $h_{1}$ is determined. As can be seen in Figure 3b, the current axial length of the structure $l_{s}$ can then finally be calculated by

l_{s} = v_{s} \cdot Δ t_{2} - [s_{c} (t_{3}) - s_{c} (t_{1})]

(2)

using the latest previously estimated velocity $v_{s}$ . Similarly, the structure’s lower boundary can also be used to determine these parameters.

By following an axially moving structure, its boundaries are tracked every time one crosses one of both image planes. Between these crossing events at times $t_{k}$ , the structure’s length is assumed to be constant. At each such event, an updated structure velocity $v_{s}$ value can be calculated and, thus, the structure’s acceleration $a_{s} = [v_{s} (t_{k}) - v_{s} (t_{k - 1})] / (t_{k} - t_{k - 1})$ . Consequently, continuous updates of the estimated structure parameters $s_{s}, v_{s}, a_{s}$ and $l_{s}$ are realized by evaluating the crossing events together with both their corresponding time-stamp $t_{k}$ and their respective CT scanner position $s_{c}$ (see Figure 7).

To distinguish between the structure’s upper and lower boundary, six different CT scanner states $S$ have been defined, as shown in Figure 4.

Figure 4.

Defined CT scanner states with corresponding cross-sectional images to distinguish between upper and lower structure boundary.

States $1, 2^{b}$ and $3$ can be classified solely from visibility information in the reconstructed cross-sections and are therefore defined as certain states. To differentiate between states $0, 2^{a}$ and $4$ , the state history of the CT scanner must be considered. These states are therefore classified as uncertain states. Possible transitions between states are shown in Figure 5.

Figure 5.

State transition graph with valid state transitions and their respective transition conditions.

Each transition is used to identify the tracked structure’s upper or lower boundary, respectively. For example, transition $0 \to 1$ indicates the structures upper boundary and $3 \to 4$ the structures lower boundary at the time of scanning. However, because uncertain states can only be determined based on the latest estimated structure properties, like $v_{s}$ and $l_{s}$ , false state estimations may occur due to disturbances like turbulence or vibrations. This may lead to invalid transitions as shown in Figure 6 (dotted line).

Figure 6.

Exemplary retroactive correction for invalid state transitions.

In case such an invalid transition is detected, the system falls back to the latest certain state at time $t_{m}$ and retroactively corrects all state estimates and structure property extractions up to the current time-stamp $t_{n}$ . Thus, the CT scanner states must additionally be acquired and saved at each state transition with their respective time-stamp. Each of these events is considered as a sample point $t_{k}$ for the parameter extraction. Based on the continuously updated structure parameters the CT scanner’s acceleration profile can be properly adapted (see Figure 7).

Figure 7.

Concept of evaluating structure’s front and back using image-pairs and the history of the CT scanner states using the example of the initial investigation of structure parameters required at the beginning of each of the aforementioned structure tracking strategies.

Motion planning for the CT scanner

Depending on the above given tracking strategy and the latest parameter estimation of the structure, a new target position $s_{c}^{target}$ is calculated for each clock cycle (see Figure 2). To reach the target height in the minimum positioning time interval $Δ t_{trav}$ , acceleration and deceleration phases have to be performed with the systems maximal applicable acceleration value $a_{c}^{\max}$ . Thus, the generated trajectory for the CT scanner’s movement needs to comply with these boundary conditions

s_{c} (t_{0} + Δ t_{trav}) = s_{c}^{target} (t_{0} + Δ t_{trav})

(3)

v_{c} (t_{0} + Δ t_{trav}) = v_{s} (t_{0} + Δ t_{trav})

(4)

| v_{c} (t) | \leq v_{c}^{\max}

(5)

a_{c} (t) \in {0, a_{c}^{\max}, - a_{c}^{\max}}

(6)

The structure’s movement profile is described by

s_{s} (t) = \frac{a_{s}}{2} \cdot t^{2} + v_{s} (t_{0}) \cdot t + s_{s} (t_{0}),

(7)

v_{s} (t) = a_{s} \cdot t + v_{s} (t_{0}),

(8)

a_{s} (t) = const .

(9)

For better readability and without loss of generality, terms in the form $s_{s} (t_{0})$ are shortened to $s_{s}$ and the initial time $t_{0}$ is set to $0$ in the following. Based on the structure’s movement profile, a similar profile for the target position $s_{c}^{target}$ of the CT scanner can be calculated using $s_{c}^{target} (t) = s_{s} (t) + s_{i}^{off}$ with $i \in {1, 2}$ as explained in section 2.2. Because the three discrete possible acceleration values $a_{c}$ for the CT scanner are predefined (see equation (6)), only the time intervals for acceleration and deceleration phase need to be calculated to provide the time-optimal trajectory. This optimal acceleration time $Δ t_{acc}$ is calculated explicitly by solving the given parabolic position and linear velocity profiles (see Figure 8)

Δ t_{acc} = \frac{dir \cdot (v_{s} - v_{c}) + \sqrt{(\frac{dir \cdot a_{s}}{2 \cdot a_{c}^{\max}} + \frac{1}{2}) \cdot [{(v_{c} - v_{s})}^{2} + 2 \cdot (a_{s} - dir \cdot a_{c}^{\max}) \cdot (s_{c} - s_{c}^{target})]}}{a_{c}^{\max} - dir \cdot a_{s}} .

(10)

Therein, the directional term $dir \in {- 1, 1}$ is used to switch the order of acceleration and deceleration depending on the current situation at each master clock event. Using the wrong $dir$ -value leads to a transition time $Δ t_{acc} < 0$ . If this is the case, the calculation for $Δ t_{acc}$ is repeated with $dir = - dir$ . If the calculated acceleration interval $Δ t_{acc}$ does not lead to a violation of the predefined maximum velocity $v_{c}^{\max}$ , the corresponding deceleration time

Δ t_{dec} = \frac{2 \cdot dir \cdot a_{c}^{\max} \cdot Δ t_{acc} + v_{c} - v_{s}}{a_{s} + dir \cdot a_{s}^{\max}} - Δ t_{acc}

(11)

can be calculated. Otherwise, an intermediate phase segment with constant and maximum velocity of duration $Δ t_{lin}$ must be added between acceleration and deceleration interval (Figure 8b). Therefore, the acceleration time interval $Δ t_{acc}$ must accordingly be updated at first by

Δ t_{acc} = \frac{v_{c}^{\max}}{a_{c}^{\max}} - dir \cdot \frac{v_{c}}{a_{c}^{\max}}

(12)

considering the current velocity $v_{c}$ of the CT scanner. Secondly, the corresponding deceleration interval is calculated by

Δ t_{dec} = \sqrt{\frac{{[v_{c} (t_{acc}) - v_{s} (t_{acc})]}^{2} + 2 \cdot a_{s} \cdot [s_{c} (t_{acc}) - s_{c}^{target} (t_{acc})]}{dir \cdot a_{s} \cdot a_{c}^{\max} + {(a_{c}^{\max})}^{2}}} .

(13)

Lastly, the intermediate time interval $Δ t_{lin}$ is calculated from the given movement profiles by

Δ t_{lin} = {\begin{matrix} \frac{v_{c} (t_{acc}) - v_{s} (t_{acc}) - Δ t_{dec} \cdot (a_{s} + dir \cdot a_{c}^{\max})}{a_{s}} for a_{s} \neq 0 \\ \frac{s_{c}^{target} (t_{acc}) - s_{c} (t_{acc})}{v_{c} (t_{acc}) - v_{s} (t_{acc})} + \frac{v_{s} (t_{acc}) - v_{c} (t_{acc})}{2 \cdot dir \cdot a_{c}^{\max}} for a_{s} = 0 \end{matrix}

(14)

Finally, equation (15) compiles the minimum traverse time $t^{(trav)}$

Δ t_{trav} = {\begin{matrix} Δ t_{acc} + Δ t_{dec} & for & | v_{c} (t_{acc}) | < v_{c}^{\max} \\ Δ t_{acc} + Δ t_{lin} + Δ t_{dec} & for & | v_{c} (t_{acc}) | = v_{c}^{\max} \end{matrix}

(15)

on which the scanner will hit its target scanning position. The trajectory is updated every master clock cycle and passed as input to the scanner’s positioning unit.

Figure 8.

Velocity profile for the CT scanner (a) without and (b) with a constant velocity segment to be able to follow a moving structure.

After presenting different tracking strategies and the parameter extraction procedure for the observation of dynamic structures, the entire control concept (Figure 9) can be concluded. As input parameters the continuously reconstructed image-pair stream as well as the current position of its imaging planes $s_{c}$ are required.

Figure 9.

Basic tracking sequence to determine a structure’s rising velocity and length using a dual-plane CT scanner.

Each time the CT scanner state $S$ is changed, the structure’s parameters are updated. Based on these extracted parameters as well as the limits of the positioning system, the CT scanner’s motion is planned using one of the three introduced tracking strategies.

Results and discussion

To investigate the feasibility of the proposed control strategies the behaviour of the entire system, that is, CT scanner and control unit as depicted in Figure 9, has been simulated for various flow scenarios. For this, an artificial 1D logical phantom vector which resembles predefined structure motion and length profiles (see Figure 10) was used.

Figure 10.

Simulated structure movements. a) Movement profiles with their respective indices (linear (1,2), accelerated (3,4), decelerated (5,6), up and down (7)) and b) length profiles (constant (1,2), linear increase (3,4), linear decrease (5,6), fluctuating (7,8)) along the test section. Dotted lines represent $Δ h_{c}$ . Quantitative definitions are given in Table 1.

To cover a wide range of technical applications different movement and structure length profiles have been simulated (see Figure 10). This set of parameters includes linear, accelerated, decelerated and up/down movement as well as constant, linearly increasing / decreasing and fluctuating structure lengths. The used values for the study are compiled in Table 1 with their respective indices. Combining all parameter variations, a total number of 4536 different parameter sets have been simulated.

During the simulation, the path and the velocity of the structure phantom is calculated for each current time step $t_{n}$ considering also the input delay $Δ t_{in}$ , which accounts for data transfer and data processing delays. According to Frust (2016), the image reconstruction latency is expected to be $4.3 ms$ at maximum for the default imaging rate of $f_{im} = 2000 fps$ . Data transfer times are estimated to be $0.2 ms$ based on the available bandwidth. For structural feature extraction an additional latency of $0.5 ms$ is assumed based on the usage of a threshold-based binarization data processing stage. The total input latency for the default image-pair rate is, therefore, estimated to be $Δ t_{in} = 5.0 ms$ (see Figure 9). A larger input delay $Δ t_{in}$ is also considered because the recognition and rediscovery of a structure may require multiple pairs of images. Furthermore, the output delay $Δ t_{out}$ , that is, the time interval between generating the acceleration signal and realizing the corresponding movement, is estimated to be $1 ms$ based on commercially available motor controllers (SEW-EURODRIVE, 2010).

In the real positioning system, the position $s_{c}$ is determined using a wire potentiometer with a measurement uncertainty of $Δ s_{c}^{meas} = \pm 0.1 mm$ . This offset was added as Gaussian noise to the current system position $s_{c}$ at each clock cycle in the simulations. An additional larger position uncertainty $Δ s_{c}^{meas} = \pm 0.5 mm$ was simulated to account for mechanical offsets, for example, vibrations or sagging of the wire, in the real setup. Height information time delay is negligible. A time step discretization, that is, cycle time of $Δ t_{cycle} = 1 ms$ is used to fit the typical acquisition rate of $f_{im} = 2000 fps$ . The CT scanner’s position limits are

s_{c}^{\min} = 500 mm and s_{c}^{\max} = 2500 mm .

The CT scanner starts at standstill, that is, $v_{c} (t_{0}) = 0$ , at an initial height $s_{c} (t_{0}) = 500 mm$ and comes to a complete stop at its maximum height $s_{c}^{\max}$ . Its maximum velocity and acceleration were set to

$v_{c}^{\max} = 500 mm \cdot s^{- 1} and a_{c}^{\max} = 500 mm \cdot s^{- 2}$ , respectively. The structure’s initial height is $s_{s} (t_{0}) = 0 mm$ . All simulations were carried out in Matlab. All parameter sets with $Δ s_{c}^{meas} > 0$ were simulated 100 times to get statistically significant information about the influence of position noise $Δ s_{c}^{meas}$ whilst retaining manageable simulation times (about 5 seconds / simulation).

Structure parameter extraction

To evaluate the tracking strategies regarding velocity information extraction, the momentary velocity estimations $v_{s} (t_{k})$ at each sample point $t_{k}$ (see Figure 5) were compared with the momentary nominal velocities $v_{s}^{nominal}$ , that is, the predefined velocities of the phantom vector, for each set of parameters being simulated. The median velocity estimation ${\tilde{v}}_{s} (t_{k})$ of each simulation was selected and averaged with the estimate of all 100 simulations of the same parameter set. The median value was selected instead of the mean value because single estimations may become implausibly large or small due to short durations between sample points or numeric limitations, for example, division by length scales close to floating point precision. Results are shown in a parity plot in Figure 11.

Figure 11.

Median momentary velocity estimates for different movement profiles (standard deviation < 5%).

Results at all sample points show an agreement between nominal velocity and median measured velocity ${\tilde{v}}_{s} (t_{k})$ with deviations below $2 %$ , which is acceptable for this application. Furthermore, velocity estimates have been identified to be insensitive to structure length $l_{s}$ , system delay $Δ t_{in / out}$ and position noise $Δ s_{c}^{meas}$ in the tested ranges.

To evaluate the length estimation quality the momentary length estimations $l_{s} (t_{k})$ for each simulation were compared with the momentary nominal structure length $l_{s}^{nominal}$ as explained previously for the velocity estimates. Median length estimation results ${\tilde{l}}_{s} (t_{k})$ show good agreement with typical deviations better than $10 %$ (see Figure 12). Length estimation deviations tend to increase for fluctuating structure lengths which results from the limited number of sample points for the length determination, that is, relatively low sample frequency compared with the fluctuation frequency. Length estimates are insensitive to system delays $Δ t_{in}$ and $Δ t_{out}$ as well as position uncertainty $Δ s_{c}^{meas}$ in the tested ranges.

Figure 12.

Median momentary length estimates for different movement profiles (standard deviations < 10%).

Comparison of control strategies

As depicted in Figure 13a, all control strategies deliver similar velocity estimation results of about ${\tilde{v}}_{s} (t_{k}) / v_{s}^{nominal} = 0.999 \pm 0.019$ averaged over all parameter combinations. However, length estimates reveal larger relative errors for the single image plane tracking (SIPT) strategy with an average result of ${\tilde{l}}_{s} (t_{k}) / l_{s}^{nominal} = 0.985 \pm 0.084$ (Figure 13b). This is due to false state estimates which cannot be corrected with the limited data provided by SIPT, especially for structure lengths shorter than the distance between the imaging planes $Δ h_{c}$ because a differentiation between states $0$ and $2^{a}$ is impossible (see Figure 4). Dual image plane tracking (DIPT) and dual image plane re-tracking (DIRPT) show more narrow distributions of length estimates with respective results of ${\tilde{l}}_{s} (t_{k}) / l_{s}^{nominal} = 0.992 \pm 0.060$ and ${\tilde{l}}_{s} (t_{k}) / l_{s}^{nominal} = 0.992 \pm 0.063$ .

Figure 13.

Comparison of feature extraction quality for different control strategies "SIPT", "DIPT" and "DIPRT". (a) velocity determination, (b) length determination.

Moreover, DIPRT’s length estimates show a bimodal distribution (see Figure 13b). Further investigations of length estimations for DIPRT in different structure length profiles reveal that this bimodal distribution occurs only for fluctuating structure length profiles (see Figure 14a).

Figure 14.

Length estimation distribution for (a) different structure length profiles using DIPRT strategy and (b) fluctuating length profiles using different control strategies.

This bimodal distribution occurs for all simulated movement profiles, delays $Δ t_{in}$ and $Δ t_{out}$ as well as position noise levels $Δ s_{c}^{meas}$ in combination with fluctuating structure length (see Figure 14b). It is especially apparent for initial structure velocities $v_{s}^{nominal} (t_{0}) = 50 mm \cdot s^{- 1}$ . The reason for this distribution is a systematic underestimation of fluctuating structures with mean length $l_{s} < Δ h_{c}$ combined with an overestimation of fluctuating structures with mean length $l_{s} \geq Δ h_{c}$ . This tendency is most apparent in DIPRT due to its narrower distribution of length estimates. In SIPT and DIPT, the wider distributions do overlap and thereby conceal this tendency when averaging all test cases.

Increasing input latency $Δ t_{in}$ leads to an overall increase in the velocity estimate’s standard deviation in each control strategy from $σ ({\tilde{v}}_{s} (t_{k}) / v_{s}^{nominal}) = 0.015$ for $Δ t_{in} \in {0, 0.005 s}$ to $σ ({\tilde{v}}_{s} (t_{k}) / v_{s}^{nominal}) = 0.035$ for $Δ t_{in} = 0.010 s$ . Variation of output delay $Δ t_{out}$ and position noise $Δ s_{c}$ do not show significant deviations ( $σ ({\tilde{v}}_{s} (t_{k}) / v_{s}^{nominal}) \leq 0.03$ ) in the simulated range of values. An acceleration switching frequency of $10 Hz$ was sufficient for all simulated parameter sets, that is, the electro-mechanical positioning system requires a corresponding cut-off frequency. The influence of different time discretization $Δ t_{cycle}$ was only briefly investigated and did not show noteworthy differences.

Conclusion

Different tracking strategies for the application of ultrafast X-ray tomography for studying moving structures in multiphase flows have been evaluated. A control scheme was proposed that provides a defined CT scanner positioning movement based on currently acquired image-pair data. Therefore, three different tracking strategies, namely Single-Image Plane Tracking (SIPT), Dual-Image Plane Tracking (DIPT) and Dual-Image Plane Re-Tracking (DIPRT), were discussed and evaluated by numerical simulation. All strategies promise very good velocity tracking results. However, SIPT is preferred for cases with constant structure shapes and lengths, for example, tracking of a tracer particle. To extract length and shape information for time-variable structure length and shapes, DIPT is more suitable for low disturbances. The DIPRT approach is preferred for detailed shape determination at high turbulence cases. Simulations have shown that DIPT and DIPRT were more robust than SIPT concerning the influence of time delays and position uncertainty.

Future work will focus on reliable structure recognition from the image pair data. Moreover, further post-processing steps to increase state estimation quality, and therefore feature extraction quality, using computationally complex techniques like particle filtering are of high interest.

Footnotes

Declaration of conflicting interests

The author(s) declared no potential conflict of interests with respect to the research, authorship and/or publication of this article.

Funding

The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: Financial support by Deutsche Forschungsgemeinschaft (DFG) is gratefully acknowledged (BI 1770/2-1).

ORCID iDs

Dominic Windisch

Martina Bieberle

References

Bai

Baidoo-Williams

Mudumbai

, et al. (2015) Robust tracking of piecewise linear trajectories with binary sensor networks. Automatica 61: 134–145. https://doi.org/10.1016/j.automatica.2015.07.012

Barthel

Bieberle

Hoppe

, et al. (2015) Velocity measurement for two-phase flows based on ultrafast X-ray tomography. Flow Measurement and Instrumentation 46: 196–203.

Bieberle

Frust

Wagner

, et al. (2017) Data processing performance analysis for ultrafast electron beam X-ray CT using parallel processing hardware architectures. Flow Measurement and Instrumentation 53: 180–188. https://doi.org/10.1016/j.flowmeasinst.2016.04.004

Bieberle

Barthel

Hoppe

, et al. (2012) Ultrafast electron beam X-ray computed tomography for 2D and 3D two-phase flow imaging. In: 2012 IEEE International Conference on Imaging Systems and Techniques Proceedings pp. 605–610. https://doi.org/10.1109/IST.2012.6295548

Boden

Dos Santos Rolo

Baumbach

, et al. (2014) Synchrotron radiation microtomography of Taylor bubbles in capillary two-phase flow. Experiments in Fluids 55(7): 1768.

Brennen

(2005) Fundamentals of Multiphase Flow. Fundamentals of Multiphase Flow. Cambridge: Cambridge University Press.

Cole

Waters

Parker

, et al. (2010) PEPT combined with high speed digital imaging for particle tracking in dynamic foams. Chemical Engineering Science 65(5): 1887–1890.

Degaleesan

Dudukovic

Pan

(2001) Experimental study of gas-induced liquid-flow structures in bubble columns. AIChE Journal 47(9): 1913–1931.

Fabich

Sederman

Holland

(2017) Study of bubble dynamics in gas-solid fluidized beds using ultrashort echo time (UTE) magnetic resonance imaging (MRI). Chemical Engineering Science 172: 476–486.

10.

Fischer

Hoppe

Schleicher

, et al. (2008) An ultra fast electron beam x-ray tomography scanner. Measurement Science and Technology 19(9): 094002. https://doi.org/10.1088/0957-0233/19/9/094002.

11.

Frust

(2016) Entwicklung und Umsetzung eines echtzeitfähigen Datenverarbeitungs- und Rekonstruktionsalgorithmus für die Ultraschnelle Elektronenstrahl-Röntgen-CT.

12.

Frust

Wagner

Stephan

, et al. (2017) Rapid Data processing for ultrafast electron beam x-ray computed tomography using scalable and modular CUDA based pipelines. Computer Physics Communications 219: 353–360. https://doi.org/10.1016/j.cpc.2017.05.025

13.

Haschke

Weitnauer

Ritter

(2008) On-Line planning of time-optimal, jerk-limited trajectories. In: 2008 IEEE/RSJ International Conference on Intelligent Robots and Systems, 22 September 2008, pp. 3248–3253. https://doi.org/10.1109/IROS.2008.4650924

14.

Hutchinson

Hager

Corke

, et al. (1996) A tutorial on visual servo control. IEEE Transactions on Robotics and Automation 12(5): 651–670.

15.

Kaestner

Munch

Trtik

(2011a) Spatiotemporal computed tomography of dynamic processes. Optical Engineering 50(12): 123201.

16.

Kaestner

Hartmann

Kühne

, et al. (2011b) The ICON beamline A facility for cold neutron imaging at SINQ. Nuclear Instruments and Methods in Physics Research, Section A: Accelerators, Spectrometers, Detectors and Associated Equipment 659(1): 387–393.

17.

Kak

Slaney

Wang

(1988) Principles of Computerized Tomographic Imaging. IEEE Press, 1988.

18.

Kim

Chen

(2007) Planning algorithms for S-curve trajectories. In: IEEE/ASME International Conference on Advanced Intelligent Mechatronics, AIM, Zurich, Switzerland, 4–7 September 2007, pp. 1–6. IEEE. https://doi.org/10.1109/AIM.2007.4412440

19.

Kopmann

Chilingaryan

Vogelgesang

, et al. (2016) UFO – a scalable platform for high-speed synchrotron X-ray imaging. In: 2016 IEEE Nuclear Science Symposium, Medical Imaging Conference and Room-Temperature Semiconductor Detector Workshop (NSS/MIC/RTSD), Strasbourg, France, October 2016, pp. 1–5. IEEE. https://doi.org/10.1109/NSSMIC.2016.8069895

20.

Liu

(2013) Target Tracking with Binary Sensor Networks. Master’s Thesis. Electrical and Computer Engineering. University of Kentucky, USA. http://uknowledge.uky.edu/ece_etds/14

21.

Müller

Davidson

Dennis

, et al. (2007) Oscillations in gas-fluidized beds: Ultra-fast magnetic resonance imaging and pressure sensor measurements. Powder Technology 177(2): 87–98.

22.

Nadeem

Heindel

(2018) Review of noninvasive methods to characterize granular mixing. Powder Technology 332: 331–350.

23.

Neunert

Crousaz

Furrer

, et al. (2016) Fast Nonlinear Model Predictive Control for Unified Trajectory Optimization and Tracking. In: 2016 IEEE International Conference on Robotics and Automation (ICRA), Stockholm, Sweden, June 2016, pp. 1398–1404. https://doi.org/10.1109/ICRA.2016.7487274.

24.

Reinecke

Petritsch

Schmitz

, et al. (1998) Tomographic measurement techniques–visualization of multiphase flows. Chemical Engineering & Technology 21(1): 7–18.

25.

Ruppel

Knierim

Sawodny

(2011) Analytical Multi-Point Trajectory Generation for Differentially Flat Systems with Output Constraints. IFAC Proceedings Volumes 44(1): 950–955.

26.

SEW-EURODRIVE (2010) System Manual - MOVIDRIVE MDX60B/61B. Notes.

27.

Tayler

Holland

Sederman

, et al. (2012a) Applications of ultra-fast MRI to high voidage bubbly flow: Measurement of bubble size distributions, interfacial area and hydrodynamics. Chemical Engineering Science 71: 468–483.

28.

Tayler

Holland

Sederman

, et al. (2012b) Exploring the origins of turbulence in multiphase flow using compressed sensing MRI. Physical Review Letters 108(26): 1–5.

29.

Wang

Member

Zhang

, et al. (2003) System identification using binary sensors. IEEE Transactions on Automatic Control 48(11): 1892–1907. https://doi.org/10.1109/TAC.2003.819073

30.

Xie

Reinecke

Beck

, et al. (1995) Electrical tomography techniques for process engineering applications. The Chemical Engineering Journal and The Biochemical Engineering Journal 56(3): 127–133. https://doi.org/10.1016/0923-0467(94)02907-5