{'fraction_non_alphanumeric': 0.036201117318435755, 'fraction_numerical': 0.010949720670391062, 'mean_word_length': 4.640966386554622, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Vision-based Control of a Soft Robot for Maskless Head and Neck Cancer Radiotherapy

This work presents an on-going investigation of the control of a pneumatic soft-robot actuator addressing accurate patient positioning systems in maskless head and neck cancer radiotherapy. We employ two RGB-D sensors in a sensor fusion scheme to better estimate a patient's head pitch motion. A system identification prediction error model is used to obtain a linear time invariant state space model. We then use the model to design a linear quadratic Gaussian feedback controller to manipulate the patient head position based on sensed head pitch motion. Experiments demonstrate the success of our approach.


I. INTRODUCTION
This paper presents a continuation of our investigation of an image-guided soft robot patient positioning system for use in head and neck (H&N) cancer radiotherapy (RT). In 2014, over 1.6 million patients developed pharynx and oral cavity cancers in the United States, which led to over 580,000 deaths [1]. Typical H&N cancer treatment involves intensitymodulated radiotherapy (IMRT), which delivers high potent dose to tumors while simultaneously minimizing dose to adjacent critical organs such as spinal cord, parotids glands, and optical nerves. Typically, a patient lies on a 6-DOF movable treatment couch, and laser or image-guidance systems are used to ensure the patient is in the proper position.
IMRT requires accurate patient positioning. An examination of patient displacement and beam angle misalignment during IMRT showed errors as small as 3-mm caused 38% decrease in minimum target dose or 41% increase in the maximum spinal cord dose [2]. Image-guided radiotherapy (IGRT) has improved IMRT accuracy while reducing set-up times [3], [4], [5]. However, current IGRT practices focus on using images acquired before treatment to confirm beam placement [6]. The discomfort caused by head masks in prolonged IMRT treatment can increase patients voluntary and involuntary motion. Studies show that translational errors caused by patient motion can be larger than 6mm, and rotational errors can be as high as 2 • [7]. Current motion-tracking systems, such as Cyberknife and Novalis are not compatible with conventional linear particle accelerators used at the majority of cancer centers. Moreover, these two systems are limited to assuming the patient's body is rigid during motion tracking and compensation. Recently, a robotic realtime surface image-guided positioning system was studied for feasibility in frameless and maskless cranial stereotactic radiosurgery [8]. While it achieved similar accuracy as the existing clinical methods, the system may not be suitable to IMRT due to the presence of mechanical and electrical parts in the path of the radiation beam.
Soft robot systems are deformable polymer enclosures with fluid-filled chambers that enable manipulation and locomotion tasks by a proportional control of the amount of fluid in the chamber [9], [10]. Their customizable, deformable nature and compliance make them suitable to biomedical applications as opposed to rigid and stiff mechanical robot components. They can also be made radiotransparent, which is necessary in IMRT.
The long term goal of our work is to address the nonrigid motion compensation during H&N RT. As we continue our initial investigation, we control one degree of freedom, raising or lowering of a generic patient's head, lying in a supine position, to a desired height above a table. The current system consists of a single inflatable air bladder (IAB), a mannequin head and a neck/torso motion simulator, two different Kinect RGB-D cameras to measure patient position, two current-controlled pneumatic valve actuators, and a National Instruments myRIO microcontroller. In this work, we extended and improve our previous work [11]. This paper contributes better vision tracking and localization methods via filtering and fusion of the two RGB-D estimates. We improve on the system identification of the soft-robot system and now incorporate an optimal control network. The result is a much improved motion control.
Section II of this paper briefly presents the design of the soft robot system. Section III discusses the computer vision algorithms to detect the patient's face and fusion of measurements from the RB-D images. Section IV presents results of system identification for the soft-robot system. Section V presents design of the linear quadratic Gaussian (LQG) controller, and Section VI presents several experiments to demonstrate the system.

II. SOFT ROBOT DESIGN OVERVIEW
The soft robot actuation mechanism combines an IAB (19" x 12") made of lightweight, durable and deformable polyester and PVC, two current-controlled proportional solenoid valves, and a pair of silicone rubber tubes (attached to a T-port connector at the orifice of the IAB) in order to convey air in/out of the IAB. A 1HP air compressor supplied regulated air at 30 psi to the inlet actuating valve, while an We use a Kinect Xbox 360, and a Kinect for Windows v2 sensor to estimate head position and velocity. The two sensors use different electronic perception technologies to determine distance of an object from the camera origin. They therefore have different lateral and range resolutions as well as different noise characteristics. Image processing for both cameras is executed on a 22GB RAM mobile workstation with Intel Core i7-4800MQ processor running 64-bit Ubuntu Trusty on a Linux 4.04 kernel.

III. IMAGE-BASED PATIENT POSITION ESTIMATION
We perform recursive filter estimations of the RGB-D measurements and improve position estimates by using an additional sensor to better localize tracked features. We add the Kinect v2 sensor (henceforth called the v2 sensor), based on the time-of-flight (ToF) electronic perception principle. In ToF, light pulses illuminate a scene, and depth is calculated by determining the phase shift of the returned light signals. The active infra-red reduces the dependence on ambient lighting [12], and this sensor has a higher spatial depth resolution of 512×424 pixels at 30Hz interactive rate, compared to the Xbox's 320 × 240 pixels [13]. To minimize the noise due to the limited sensor resolution, the v2 has in-built noise improvement capabilities [14].
The v2 provides a higher depth-map accuracy and lower noise floor compared to Kinect for Xbox, as can be seen from Fig. 1, where the v2 exhibits a noise auto-covariance of 11.4707mm 2 compared with 22.7057mm 2 for the Xbox. Despite the improved performance of the v2, noise remains an issue, as is the case for every electronic perception system. To alleviate this, we employ a multisensor data fusion of both Kinect sensors' observations. We achieved this by local Kalman Filter estimates of each sensor's observations, and we fuse the estimates via a variance-weighted multisensor Kalman filter fusion scheme described later in this section.

A. Face Detection and eye-feature tracking
We approached face detection using Haar Cascade Classifiers (HCC) [15]. HCC's are based on integral image representation, which allow for features evaluation while maintaining high detection rates. The features resemble Haar basis functions. A classifier is formed by choosing a small number of crucial features with AdaBoost, and a weighted sum of individual classifiers is used to construct a strong object detection classifier in a cascade manner. This increases the detector's speed by concentrating on areas within an image with high probability of features of interest.
A drawback of HCC's is the memory consumed on computing devices when searching through image pixels for specific regions of interest. Searching through a 640 × 480 pixels gray-scale image for specific features caused a 90% reduction in the frame rates of either sensor, when the algorithm is run on a CPU. To overcome this, both sensor's images were spatially down-sampled via linear interpolation before HCCs were applied. Face detection was performed on a single NVIDIA Quadro K1100M GPU. We retrieve each detected face from the GPU, and then detect eyes within detected faces using the same procedure.
To achieve robust detection, the minimum number of neighbors in each candidate rectangle feature was determined based on our experience. The search area within an image was chosen to be within the range of (5 × 5) pixels and (20 × 20) pixels. This gave us more than 90% face detection rate for both sensors. A similar approach was used for the eye classifier. The final implementation achieved a frame rate of 15Hz for each sensor running independently on the Linux host computer. Further improvement in frame rates is an avenue for future work.

B. Local Kalman Filters
From Fig. 1, we see that both RGB-D sensors suffer from notable associated noise, which is not suitable for our control requirement. To refine the observation, local Kalman Filter (KF) estimates for each sensor were computed to determine state estimatesx(i) that minimizes the meansquared error to the true state x(i), given a measurement sequence z(1), · · · , z(j), that iŝ where the obtained estimate is the expected value of the state at time i given observations up to time j. The covariance of the estimation error is given by Assuming the model of the state is common to both sensors, and denoting the distance from the v2 to the head as d(k), we define x(k) = [d(k),ḋ(k)] T ∈ R 2 as the state vector of interest, and let ∆T be the time between steps k−1 and k. The model state update equations are given by where F(k) ∈ R 2×2 is the state transition matrix given by is the control input matrix that maps inputs to system states, G(k) ∈ R 2×2 process noise matrix, and w(k) ∈ R 2 is a random variable that models the state uncertainty. In the absence of inputs B k u k = 0, and the model becomes where w k is the effect of an unknown input and G k applies that effect to the state vector, x k . The process noise is assumed unknown and is modeled as uncontrolled forces causing an acceleration a k in the head position (a k is thus a scalar random variable with normal distribution, zero mean and standard deviation σ a ). We model this into (3) by setting G k to identity and set w( . Therefore, we find that Denoting the head displacement at time k as measured by the Xbox and v2 as z 1 (k) and z 2 (k) respectively, the sensors' measurements were mapped to the v2 reference frame and modeled as where H s (k) = 1 0 T maps the system's state space into the observed space, and v s (k) ∈ R is a random variable  Fig. 3: KF results for the Xbox observation that models the sensor error. We define v s (k) as a normally distributed random variable with zero mean and variance σ 2 rs . We assume the random sequences v 1 (k),v 2 (k), w(k) are independent and uncorrelated in time.
At each time step, k, each local KF's priori and posteriori estimates are computed through the typical prediction and update phases Prediction Phase: wherex k|k−1 and P k|k−1 are the state prediction vector and the prediction covariance matrix respectively. Update Phase: where K k ,x k|k , and P k|k are respectively the KF gain, posteriori state estimate and its state covariance matrix. In implementing the KF of (8) and (9), the variance of the process noise/signal noise of each local KF was informed by our knowledge of the physics of both sensors (electronic perception methods, range resolutions and examining each sensor's depth map to understand the data available to the filter), engineering judgment, and kinematics of the process model. We found these values sufficiently modeled the underlying process dynamics

C. Data Fusion
Each local KF estimate was combined at a central fusion site to obtain a track-to-track fused global estimate. To communicate each estimate and associated covariance matrix, we create Unix FIFO special files (i.e. named pipes) on the kernel file system, write the estimates and covariance matrices to the pipes at each local site and retrieve the values at the central site. Named pipes are low-level file I/O systems that can be shared by processes with different ancestry. During data exchange through a FIFO, the kernel forwards all data internally without having to write it to the file system. Since they exist within the kernel and the file system is just an entry serving as a reference point for the processes to access the pipe with a file system name, there is practically no delay in data communication.
Local tracks are generated at each sensor site according to (8), resulting in two local state predictions from the Kalman filters (3). At the central fusion site, we assume a state model common to both sensors given by (9) and adopt a varianceweighted average of each local track in the global track fusion algorithm [16] x F (k|k) = P F (k|k) Fig . 5 shows the output of the fusion scheme compared against the single Kalman filters during a head-raising motion. The fusion of the local tracks produces better estimates, with improved signal to noise ratio. The fused estimate assigns more weight to the less noisy signal from Kinect v2. Through the implementation of the local tracks and a global track KF estimator, we improved the accuracy of the effective signal to be used in our control algorithm to no more than a standard deviation of 0.75mm from the true position of an object. The noise spikes in the fused tracks when the process state estimates are yet to converge as noticeable in Fig. 5 can be attributed to the noisy initialization of pixels in the sensors before they attain their steady state values. On average, it takes approximately 30 seconds for the pixel values in the Kinect sensor to reach their final steady state values [17]. This can be avoided by running the fusion algorithm for at least 2 minutes before the fused signal is used for any control purposes. The code for the multisensor fusion experiment is available on the git repos [18] and [19].

IV. SOFT ROBOT SYSTEM IDENTIFICATION
We approach the modeling procedure with an identification prediction error (PEM) approach, where we estimate a mathematical model, G(t), based on the minimization of the sum of squared errors between estimates of the head height, y(t), and true head height, y(t), from the fusion i.e.
Z N = {u(1) · · · u(N ) y(1) · · · y(N )} is the vector of past input and output (fused estimates) measurements over a bounded interval [1, N ] and θ is the greedy vector of parameters that approximate the model we seek to build. (11) is a special case of the least squares criterion.

A. Model Structure
Following Ljung's formulation in [20, §4.5], we pose the identification problem as determining the "best model" from a set of candidate model sets via an iterative approach that parametrizes the noncountable model sets smoothly over an area with the assumption that the underlying system is linear time-invariant. Here, our model structure is a differentiable mapping from a connected, compact subset D M of R d to a model set M * , such that the gradients of the predictor functions are stable. This procedure is included in the MATLAB system identification toolbox, and since the method is welldocumented in [21] we omit details.
External disturbances and stochastic variables are modeled as additive white noise sequence, e(k), based on lagged inputs and outputs, and our objective is to estimate a stochastic state space model structure of the form where the noise terms w(k) and v(k) compensate for the effect of disturbances beyond frequencies of interest to system dynamics and make the model robust to model uncertainties.
Since u and y alone are measurable in our setup, the states x(k) are estimated and (12) becomes a linear regression problem, where all the unknown matrix entries are linear combinations of the measured inputs and output variables. This can be written as where We assume the noise term is white in order to assure an unbiased model. The parameter estimation problem is then to estimate the A, B, C, and D matrices by the linear least squares regression of (13) assuming no physical insight into the system (i.e. a black box model). E(w(k)) and E(v(k)) are estimated as a sampled sum of squared errors of the residuals.

B. Parameter Estimation
The input, u(k), and output signals, y(k), can be characterized by a linear difference equation of the form where e(k) describes the equation error as a moving average of white noise, and we assume e(k) has a bias-variance term λ. We can rearrange (14) using the vectors The adjustable parameters of (15) are elements of θ. In our prediction model, it is convenient to write (14) as a one-stepahead predictor of the form with G(q, θ) = B(q) which is a complete autoregressive moving average with exogenous input (ARMAX) model. G(q, θ) represents the transfer function from input to output predictions, and H(q, θ) denotes the transfer function of prediction errors to the output model,ŷ(k); q is the z-transform, z −1 , while A(q), B(q), and C(q) are polynomials defined as A(q) = 1 + a 1 q −1 + · · · + a na q −na , C(q) = 1 + c 1 q −1 + · · · + c nc q −nc (18) [22]. The predictor turns out to be a linear filter of the form and where H(q, θ) is the noise model andŷ(k) above can be regarded as the one-step ahead predictor. After rearranging (19),we find that such that the residual errors from (19) become We can consider (21) as passing the prediction errors through a linear filter that allows extra freedom in dealing with nonmomentary properties of the prediction errors. Since the model is that of a linear system, (21) satisfies our objective by approximating the prefilter with the choice of the noise model in (13). The estimation problem is to predict the estimates,ŷ(k|θ) so that the errors, ε(t, θ) = y(t) −ŷ(t|θ) p are minimized by the choice of an appropriate p-norm criterion function, such as the mean squared error proposed in (11).
1) Input Signal Design: The input signal choice for a system identification experiment will determine a system's operating point and model accuracy. Therefore, the input should be rich enough to excite a system and force it to show properties needed for the model's purpose. For the model to be informative across all the desired frequency range, a periodic, persistently exciting uniform Gaussian White noise (UGWN) signal with clipped amplitudes corresponding to the bandwidth of the valves was designed offline, and its frequency spectrum analyzed to ensure it had as small a crest factor as possible (since the asymptotic properties of the model will be mostly influenced by the spectrum rather than the waveform's time-series shape). Gaussian White Noise signals (GWN) and Pseudo-Random Binary Signals (PRBS) are well-known to achieve virtually any signal spectrum without very narrow pass bands. Therefore, pseudo-random uniform white noise sequences were generated using the very-long cycle random number generator algorithm. Given that the probability density function, f (x), of the uniformly distributed uniform white noise is where A is the amplitude. The expected mean, µ, and the expected standard deviation, σ of the sequence are [23] µ = E(x) = 0, The spectrum of the resulting signal in Fig. 6 gives good signal power, which nicely relates to the bandwidth of the pneumatic valves and achieves virtually all signal spectrum with little narrow pass bands.
We therefore use the signal of (22) to model the desired asymptotic estimates,ŷ(t), of (17). We sampled y(t), the fused measurement described in (III-C) well-above the system's Nyquist frequencyand acquired enough samples to make Z N asymptotically approachθ N as N → ∞. The data collection procedure closely follows that described in our previous paper and we refer readers to [11,§IV.A] for a more detailed treatment.
The collected data was separated in a 60:40% ratio for training and testing purposes, respectively, to assure a training model that generalizes well.
2) State Space Realization: If we definê it follows that 1) as N → ∞, there are n-th order minimal state space descriptions of the system if and only if the rank of the matrix of prediction vectors,Ŷ , is equal to n for all r ≥ n; and 2) the state vector of any minimal realization in innovations can be chosen as linear combinations ofŶ r that form a row basis forŶ , i.e., with L being an n × pr matrix (p is the dimension of y(k)) [20, §7.3]. The true prediction is given by (17) with innovations e(j) written as a linear combination of past input-output data. The predictor can thus be expressed as a linear function of u(i), y(i), i ≤ k − 1. In practice, the predictor is approximated so that it depends on a finite amount past data such as s 1 past outputs and s 2 past inputs of the form y(k|k − 1) = α 1 y(k − 1) + · · · + α s1 y(k − s 1 ) (24) + β 1 u(k − 1) + · · · + β s2 u(k − s 2 ).
Piping the identification data through the MATLAB function 'ssest' and testing various model orders based on the ranking of singular values of the Hankel matrix of inputoutput measurements [20], we obtained the results listed in In the frequency-domain, this is the equivalent to having pole-zero cancellations for higher-order models. We therefore conclude there is no useful properties a higherorder model could predict beyond an order of 8. The secondorder model sufficiently approximates the system and is not significantly outperformed by the higher order models -which would contribute higher complexity to the control design. We therefore pick the 2nd order state space model (12) as where T s is the sampling period, e(k) is the modeled zeromean Gaussian white noise with non-zero variance,

V. LQG CONTROL
We employ a LQG controller and estimator to minimize the following cost function subject to the state equation (26) where K is the terminal sampling instant, Q is a symmetric, positive semi-definite matrix that weights the n-states of the A matrix, N specifies a matrix of appropriate dimensions that penalizes the cross-product between the input and state vectors, while R is a symmetric, positive definite weighting 1 where ∆u is the future control sequence and the first element in the sequence is used in the control law at every time instant. We model additive white noise disturbances into the discrete estimator's states; therefore the optimization problem becomes a stochastic optimization problem that must be solved.
The separation theorem ensures that we can construct a state estimator which asymptotically tracks the internal states from observed outputs, y(k), using the algebraic Riccati equation given as where P is an unknown n × n symmetric matrix and A, B, Q, and R are known coefficient matrices as in (26) and (27). We find an optimal control law by solving the minimization of the LQ problem, (27) which we then feed into the states.
In practice, it is a good idea to start with an identity matrix, Q, a zero penalty matrix, N , and tune R till one obtains convergence by the state estimator. The following optimal values were used after a heuristic search We construct a full online estimator for the identified plant as in Fig. 7, whereby the noise processes are assumed to be independent, white, Gaussian, of zero mean and known covariances. The optimal controller gains, K opt , are determined from the equation [25] where P is the solution to the algebraic Riccati equation (29) and E[w(k)w (τ )] = R(k)δ(k − τ ). Therefore, the online optimal estimate,x(k + 1) of x(k) iŝ wherex(k 0 ) = E [x(k 0 )] The observer is equivalent to a discrete stochastic Kalman filter that estimates the optimal statex(k|k) as shown in Fig. 7. The estimator equations are similar to equations (8) and (9) and the online, unbiased estimate iŝ Through heuristics, we found the following variances of the online estimator to be useful:

VI. EXPERIMENTAL RESULTS AND DISCUSSION
The control network was implemented on an NI-myRIO running LabVIEW 2015. We initialized the Kinect sensors to allow for all pixels within the depth cameras to reach steady state under ambient light. We performed multiple experiments to evaluate the developed state space model of IV-B and LQG controller of V. 1 . The input variable is the current that excites the valve, which in turn actuates the bladder; the head moves in response to bladder actuation. The fused estimate of the Kinect sensors are used to estimate the real-time head pitch motion as described in III-C; this is in turn used in a feedback to the LQG controller. Fig. 8 shows the results from a constant reference trajectory, which the head is meant to track. We notice a settling time of approximately 24 seconds before we reach steady state. The delay arises from our design requirements and is not a drawback in clinical trajectory tracking where we must ensure smooth head motion to desired target. It is also seen that the controller exhibits relatively smooth tracking within a 1.5 mm standard deviation over time after a relative overshoot of 5mm in bottom graph of Fig. 8. The overshoot can be explained by the estimator's search for a steady state region based on the time it takes for the pixel values of the sensors to reach steady state. The controller tracks the reference to within ±2mm.
However, we noticed an inconsistency at certain operating ranges in the current LTI model. The applied current based on fusion feedback occasionally reaches a steady state error, as can be seen from Fig. 9. We conjecture this is due to an unmodeled nonlinearity at the inlet valve that maps input currents to system states. To better approximate the nonlinearity from input to output, we will investigate using a Hammerstein block-structured model that better approximates the nonlinearity from inputs to states and states to output of the system.
The fusion algorithm proved useful to cancel jitter in the depth measurement of the sensor over our previous results,  Fig. 9: LQG Controller on Manikin Head but it falls short of the 1mm accuracy in head and neck cancer RT. Having established proof of concept in this investigation, we will begin investigation of better head localization by using sophisticated motion capture systems or laser scanners, such as those actively employed in IGRT.

VII. SUMMARY
We have presented a continuation of or initial investigation into a pneumatically-driven soft robot system for head and neck radiotherapy. Measurements from two Kinect RGB-D cameras are fused to refine the accuracy of observations. System identification on the soft robot was combined with LQG controller design to provide an optimal controller. Experiments showed we could actuate the patient head within 2.5mm accuracy. Future work will investigate accurate multiaxis positioning using better 3D sensor, control of multiple IABs, and nonlinear modeling or model free-methods to overcome limitations of the LTI model.