Fast CSF MRI for brain segmentation; Cross-validation by comparison with 3D T1-based brain segmentation methods

Objective In previous work we have developed a fast sequence that focusses on cerebrospinal fluid (CSF) based on the long T2 of CSF. By processing the data obtained with this CSF MRI sequence, brain parenchymal volume (BPV) and intracranial volume (ICV) can be automatically obtained. The aim of this study was to assess the precision of the BPV and ICV measurements of the CSF MRI sequence and to validate the CSF MRI sequence by comparison with 3D T1-based brain segmentation methods. Materials and methods Ten healthy volunteers (2 females; median age 28 years) were scanned (3T MRI) twice with repositioning in between. The scan protocol consisted of a low resolution (LR) CSF sequence (0:57min), a high resolution (HR) CSF sequence (3:21min) and a 3D T1-weighted sequence (6:47min). Data of the HR 3D-T1-weighted images were downsampled to obtain LR T1-weighted images (reconstructed imaging time: 1:59 min). Data of the CSF MRI sequences was automatically segmented using in-house software. The 3D T1-weighted images were segmented using FSL (5.0), SPM12 and FreeSurfer (5.3.0). Results The mean absolute differences for BPV and ICV between the first and second scan for CSF LR (BPV/ICV: 12±9/7±4cc) and CSF HR (5±5/4±2cc) were comparable to FSL HR (9±11/19±23cc), FSL LR (7±4, 6±5cc), FreeSurfer HR (5±3/14±8cc), FreeSurfer LR (9±8, 12±10cc), and SPM HR (5±3/4±7cc), and SPM LR (5±4, 5±3cc). The correlation between the measured volumes of the CSF sequences and that measured by FSL, FreeSurfer and SPM HR and LR was very good (all Pearson’s correlation coefficients >0.83, R2 .67–.97). The results from the downsampled data and the high-resolution data were similar. Conclusion Both CSF MRI sequences have a precision comparable to, and a very good correlation with established 3D T1-based automated segmentations methods for the segmentation of BPV and ICV. However, the short imaging time of the fast CSF MRI sequence is superior to the 3D T1 sequence on which segmentation with established methods is performed.


Introduction
Active gas bearings are an interesting alternative to more commonly used bearings such as ball and magnetic bearings.Compared to ball bearings, active gas bearings have extremely low friction.At the same time, active gas bearings are open loop stable for appropriately low rotational speeds, unlike magnetic bearings, which are always open loop unstable.The active gas film itself delivers low damping.Because of the low damping, the rotor system may become unstable due to self-excited rotor whirling.Control is therefore still required although it is less limited than for unstable systems.Designing models of gas bearings has, however, been shown to be less straightforward, compared to the frictionless alternative of magnetic bearings.3][4] Recently, a low order model able to describe the dynamics of the active gas bearing was presented in Theisen et al. 5 Such low order models greatly simplify the control design phase.
Identification of gas bearings has been a subject of some interest recently. 6,7It has been shown possible to identify appropriately low order models that mimic the behaviour of the gas bearing in a satisfying manner when no feedback loop is active.Usually the gas bearing will operate under conditions at which the damping is so small that feedback control is needed for safe operation.It is therefore of extreme relevance to be able to identify an appropriate model for the gas bearing under conditions where it is only possible to obtain identification data while the gas bearing is part of a closed loop system.
There are plenty of different methods for identifying closed loop systems.One of such methods transforms the closed loop identification problem of the plant into an open loop identification problem by using the Youla parametrisation.Well-known subspace identification methods are in this article used for identifying the Youla deviation system.It has become increasingly popular to use subspace identification methods due to their natural connection to multiple-input and multipleoutput state space models. 8Some of the most wellknown subspace identification methods are the Numerical algorithm for Subspace IDentification (N4SID) proposed in Van Overschee and De Moor 9 and the Multivariable Output-Error State-sPace (MOESP) proposed in Verhaegen and Dewilde. 10Both methods, as with most subspace identification methods, assume the noise and input to be uncorrelated; hence, the methods are based on assumptions that are only valid for open loop identification.Some work has been done on extending the subspace identification methods to identification of closed loop systems.Several different methods for coping with feedback connection have been proposed. 8However, often such methods show the same weakness of depending on being able to identify the noise signal imposed on the system.
One possible indirect identification method is subspace identification of the Youla deviation system.The Youla deviation system is used to reformulate the closed loop to open loop identification problem using a coprime factorisation of the nominal plant model and the controller.This identification procedure is known as the Hansen scheme and was first introduced in Hansen et al. 11 The original version of the Hansen scheme used external excitation signals to indirectly excite the Youla deviation system.A modified version of the Hansen scheme, first introduced in Sekunda et al., 12 is used instead in this article, which makes it possible to directly impose an excitation signal onto the Youla deviation system.A gas bearing easily becomes an unstable system when the rotational speed is increased.However, the Youla deviation system identified by the modified Hansen scheme will always be stable.The open loop identification technique is thus applicable where direct open loop identification of the gas bearing is not.The modified Hansen scheme, however, has the drawback of loss of physical understanding of the system.
The key contribution of this article is to offer insight into when it can be advantageous to identify an active gas bearing using the modified Hansen scheme instead of the more traditional Prediction Error Method (PEM) identification.This article thus focuses on comparing known identification techniques for identification of a gas bearing when being part of a closed loop system.The results are gathered by a combination of simulations and experiments conducted on a laboratory installation.
The article is structured as follows: in the following section, the model of the gas bearing is introduced; in the third section, some preliminary theory is introduced in order to let the article stand alone; in the fourth section, the three methods used for identification of the gas bearing are discussed; in the fifth section, the identification procedure is given and the quality of the nominal model is discussed; in the sixth section, the experimental methods are presented, and the resulting model identified using each of the three different identification methods is shown.Finally, in the last section, a discussion of the results and future possible improvements is given.

Gas bearing system
6][7] These methods, however, have all been both conducted and verified under open loop conditions.All identification in this article are based on the gas bearing being subject to feedback control.The nominal model is based on the model structure introduced in Theisen et al. 5 which showed it possible to model the vertical and horizontal positions of a disc attached to a flexible shaft.The shaft is held in position by the gas bearing, using gas injected into the bearing in both the horizontal and vertical directions using a sixth-order model.A picture of the test rig is presented in Figure 1(a).
An illustration of the gas bearing actuators is shown in Figure 1(b).The air is going into the bearing through the small pipe; a piezo electric actuator is attached, which controls how much air is able to get into the bearing, and thus controls the preasure inside the bearing.
The model of the gas bearing is based on the massspringer-damper model in equation (1) for the nonrotating case Here, l is a vector with the position in the horizontal and vertical directions, and m is the state of the actuator dynamics.K is the specific stiffness matrix given in ½N=kg mm, D is the specific damping matrix in ½Ns=kg mm and B is the actuator gain matrix in ½N=kg V. Using equation (1), the state vector is defined in equation ( 2) The structure of the state space model is given in equations ( 3)- (6), where each element is a 2 3 2 matrix The matrix P is a diagonal matrix with each element p j , defined by equation (7) h j is the first order low pass filter through which the system is actuated.The state space model is thus given in equation ( 8) In order to conduct closed loop identification of the system, the parameters of the model have been identified with the disc not rotating as part of an open loop scheme.The identified matrices are given in equations ( 9)-( 12 Theory The first part of this section is dedicated to presenting some definitions to the reader in order to make the rest of the article stand alone.A nominal plant is given by G(0) and a stabilising controller is given by K.The coprime factorisation of G(0) and K is given by Here, N and M denote the right coprime factorisation of the nominal plant, and Ñ and M denote the left coprime factorisation of the nominal plant.Equivalently, the right coprime factorisation of the controller is given by V and U, and the left coprime factorisation is given by Ṽ and Ũ such that the Bezout identity given in equation ( 15) is satisfied It will be assumed in this article that K is a stabilising controller for both the nominal plant and the real plant G(S).This is the case if equation ( 15) is true.Here, S is the Youla deviation system describing the divergence between the nominal model and the real plant.With the controller and plant model factorisations defined in equations ( 13) and ( 14), the real plant is given by equation ( 16) using the right coprime factorisation and equation ( 17) using the left corpime factorisation With a formulation of the true plant as in equation ( 17), only the Youla deviation system (S) is unknown and needs to be identified.G(S) from equation ( 16) can be formulated using a linear fractional transformation (LFT) given by where Here, the input and output signals of J G are presented in equation ( 19) The closed loop system using the system description in equation ( 18) is shown in Figure 2, where v 1 and v 2 are two possible excitation signals and n is the noise.
Using equation ( 18) and the setup in Figure 2, the vector introduced, e, is given in equation ( 20) Equation (20) shows that the identification of S is straightforward using the signals h and e.The two signals e and h in equation (20) are not directly accessible.However, it is shown that using the available signal vectors v 1 , v 2 , y and u, it is possible to generate e and h.It is therefore not possible to inject h directly for the open loop identification of S. 11 Instead, including the Youla parametrisation of all controllers in the closed loop, it will be possible to have a direct access to e and h as shown further on.Let the Youla parametrisation of all stabilising controllers for the nominal system G(0) be given by the following LFT description where The input and output signals of J K are in equation ( 22) Again it is possible to introduce two new signals in order to determine The two new signals b and a are part of the controller and can thus be directly measured and inserted.The relationship between the signals in the plant and controller can therefore give a significant advantage for direct identification of S.
The relationship between the signals h, e, b and a can be calculated as a Redheffer star product as defined in Zhou and Doyle. 15The resulting relationship is shown in equation ( 24) Here, H denotes the Redheffer star product.The Youla deviation system is therefore possible to identify using the signals found in the controller, which is used for the closed loop identification.Based on the system description given in Figure 2, the Youla deviation system can thus be identified using equation (25 This ability to directly impose the excitation signal on the Youla deviation system S is the advantage of the modified Hansen scheme compared to the original version.The modification is believed to make the identification more intuitive and was first introduced in Sekunda et al. 12

Identification methods
In all tests conducted, the plant is subject to a feedback loop with an a priori designed controller.In order to conduct a proper comparison between the methods, all methods use data from the same experiment.Furthermore, the noise (n) has been found to be approximately white Gaussian from an open loop data set without any excitation.

Identification using the modified Hansen scheme
It was argued in Sekunda et al. 12 that the signals should be obtained directly from the controller by modifying it to the form, shown in Figure 3.
Using the signals a and b, it is possible to cast an identification problem as given in equation ( 25) by taking advantage of the relationship found in equation (24).Since a does not depend on the noise n as shown  in Anderson, 16 it is clear from equation (25) that the noise and excitation signal are not correlated.The system S is identified using the subspace method N4SID for identification.The choice of using a subspace method for identification is twofold.One of the great advantages of using subspace identification is the direct identification of Multiple Input Multiple Output (MIMO) state space models.The other important advantages of choosing a subspace identification method are that the methods are based on zero a priori knowledge.Unlike the true system G(S), it is rather difficult to determine a proper initial structure of the Youla deviation system S.This comes from the fact that the nominal system is believed to be the best initial guess possible.All prior knowledge therefore points towards the Youla devitation system being 0.

Direct grey box identification
This identification method, as the name suggests, does not use any information regarding the feedback loop and the controller to determine the model.The method was presented in Ljung 17 and is based on obtaining input and output data of the plant and determining the model.In this article, it is implemented as a PEM.The method uses the input-output data together with a nominal model as the information for parameter identification.

Preliminary open loop identification
Identification is conducted in two separate steps.First, a nominal model is generated to use as a basis for identification of the plant when part of a closed loop scheme.Such a model is obtained using data acquired from an open loop experiment conducted when the active gas bearing is not rotating.The identification results for the nominal model are based on the method described in Theisen et al. 7 The identified nominal model is able to predict the displacement to an acceptable degree.The parameters identified through the experiment were presented in the first section.In order to give a quantitative measure of the quality of the different identification procedures used during the experiments, a goodness of fit is calculated using equation (26), which is the normalised root mean square error Here, y is the measured output, y is the mean of the measured output and ŷ is the estimated output using the model.Based on the nominal model, a controller is designed and a coprime factorisation is conducted.In the experiment, only the left factorised form is used.In order to show the effect of the closed loop identification, the nominal model has been degraded in several different ways to examine the ability of the different identification procedures when the plant is part of a closed loop scheme.Finally, a controller has been designed using the Linear-Quadratic Regulator (LQR) function in MATLAB and implemented using the controller design presented in Figure 3.

Simulation results
It is in this section investigated which benefits there are from identifying the Youla deviation system.Identifying the Youla deviation system is compared to other more well-known methods.It is decided to compare the method with grey box PEM identification and subspace identification.Normal subspace identification is independent of the nominal model.Grey box PEM identification on the other hand uses information from the nominal model.

Inferior a priori knowledge of the system dynamics
A definition of the inferior model from which the system has to be identified is needed to compare the identification methods.For the first test, inferior knowledge about the system dynamics has been investigated.The system matrix of the initial model is defined in equation ( 27), where A model has inferior knowledge relative to the real system matrix A Here, A is the system matrix of the system to be identified, A model is the system matrix of the nominal model and u is a uniformly distributed random variable in the interval 0-1.The nominal model thus contains a system matrix, where all elements are smaller than for the real plant.A simulation has been conducted, where 500 different initial models were constructed using equation ( 27).For each model, a controller was designed that would stabilise the initial model and the real plant.Noise has been added as output noise and the system has been excited using a series of square waves.For each initial model, an input-output set has been gathered and an identification has been conducted.In each case, a second verification data set has been obtained, where noise was omitted from the system.The identification results are shown in Figure 4 for identification using PEM and subspace identification of the Youla deviation system.
It can be seen from Figure 4 that for small deviations between the real system and the initial model, both methods produce similar results.However, as u is increased, the PEM identification starts to produce a lower R 2 fit.
With the methods compared, it is concluded that identification of the Youla deviation system using subspace identification is resulting in higher R 2 fit when sufficient a priori knowledge of the system dynamics is lacking.Identifying the Youla deviation system instead of the plant directly increases the complexity.The effect has therefore been examined when compared to direct subspace identification of the plant.Again 500 initial models were constructed using equation (27).Identification was conducted using subspace identification to identify the plant directly and to identify the Youla deviation system.A box plot of the resulting fit of the identified models is shown in Figure 5.
It is seen from Figure 5 that when the plant is directly identified, using a subspace identification method, the variance of the R 2 fit increases.Furthermore, the mean of the R 2 fit is higher when identifying the Youla deviation system.Identifying the plant indirectly by identification of the Youla deviation system thus produces better results.Since both methods suffer from the loss of physical understanding of the models produced, it is decided not to look further into direct subspace identification.

Experiments and results
With the controller designed and nominal model identified using the method described in Theisen et al., 5 experiments have been conducted under three different conditions.The results have been summarised in Table 1 for all three experiments.In all experiments, the same two identification methods have been applied: PEM using a fixed a priori model structure and identification of the Youla deviation system.The identification signals are found as shown in equation (28) for direct identification of S. The system has been excited using the input signal denoted a in Figure 3 The Youla devitation system is identified based on equation (25).The first experiment is conducted, where the disc is nonrotating, and the nominal model has  been degraded so that it is not able to mimic the plant dynamic in a satisfying manner.The degradation of the nominal model was achieved by multiplying all elements of the system matrix with 0:5.The result of the identification is shown in Figure 6, where the nominal model and the two identified models have been compared with a verification data set.
It is shown in Figure 6 that with a poor nominal model, identification of the Youla deviation system is able to recover the system dynamics.It is also noticed that the PEM is unable to identify a model of the same quality and obtains a lower R 2 fit when the nominal model is of such low quality.The same trend was observed in section when u is equal to 0.5.The result is summarised in column 2 in Table 1.
To further examine the effect of the identification methods for low-quality nominal models, an experiment is conducted, where the disc is rotating with 2500 r/min.Here, it has instead been chosen to multiply all elements of the input matrix with 0:7 so that the input gain is higher than expected.The verification data together with the predicted output for each of the models are shown in Figure 7.
The nominal model is better than in the previous example at predicting the position of the disc; however, the nominal model is still having a fit lower than what is expected possible to obtain.This is seen from the relatively large increase in the fit for each of the identification methods.The result is summarised in column 1 in Table 1.Finally, an experiment is conducted to  Comparison of a verification signal (blue) with the predicted output given the nominal model (red) and each of the three identified models.The nominal model G was found to have a fit of 3.188% in the horizontal direction and 2.999% in the vertical direction.The PEM (yellow) was found to have a fit of 48.91% in the horizontal direction and 39.47% in the vertical direction.Identification of the Youla deviation system (purple) gave a fit of 76.76% in the horizontal direction and 83.58% in the vertical direction.
examine the identification methods' ability to recover small deviations.Again the nominal model has been tampered with to degrade its performance.The elements of the system matrix have been multiplied by 0:97 so that the dynamics of the real plant is slightly faster than for the nominal model.Identification of the Youla deviation system using equation (28) has thus been conducted, and a plot of the signal b measured and predicted using the identified system is shown in Figure 8.The nominal model and plant are fairly similar because noise is expected to dominate the b signal.This is seen in Figure 8, where the measured signal b is presented.It is clear that the signal contains a dominant noise part in the horizontal direction.It is clear from Figure 8 that the identification should improve the model in the vertical direction mainly.Again, the identified model has been compared with a verification signal, which is shown in Figure 9.
It can be seen from Figure 9 that the nominal model does indeed produce a better fit.Still, both methods are able to improve the model.As clearly seen from Figure 8, the improvement of the identified model was mostly in the vertical direction.A zoom on the verification experiment is shown in Figure 10.
It is clear from Figure 10 that the identified models have improved the prediction accuracy in the vertical direction.Both models show a clear improvement compared to the nominal case as shown in the third column of Table 1.The results of the three experiments for each of the identification methods using equation (26) to determine the quality of the identified model are shown in Table 1 together with the fit of the nominal model.
As seen in Table 1, using PEM or subspace identification of the Youla deviation system based on the modified Hansen scheme (G(S)) gives similar results for small deviations.However, in the case where the nominal model lacks significant insight to the plant (low fit for nominal model), the PEM results in a lower fit.

Model reduction of the identified models
The models identified using the modified Hansen scheme are all of a inconveniently high order.This is a well-known result of the identification method, and it is therefore found useful to investigate the impact on the quality of the identification when model reduction is conducted.Model reduction has, in this article, been conducted by requiring all models to be of the same order as the nominal model (sixth order).The same  systems as identified previously have been compared in Table 1, where those systems identified using the Hansen scheme have been reduced to be of sixth order.
The identified plant after model reduction is given in equation ( 29) Here, modred() denotes the model reduction function used, and G(S) is the originally identifed plant using the modified Hansen scheme.The results of the models for which model reduction has been conducted are shown in Table 1.The model reduction technique used is described in Varga, 18 and the resulting sixthorder approximations are shown to have a similar fit to the full order identified systems as shown in Table 1.It is worth noting that the initially identified models had as much as 40 states, which makes control design problematic.

Conclusion
Identification of a gas bearing as part of a closed loop system has been conducted.It was shown possible to identify the gas bearing, as part of a closed loop system, by reformulating the problem into identification of the Youla deviation system which is a standard open loop identification problem.Furthermore, using an open loop PEM method for identification gives similar results when the nominal model predicts the system dynamics with high accuracy.A simulation example was given which pointed towards the identification of the Youla deviation system to produce better results on average for poor knowledge of the plant dynamics.
Experiments showed that for a poor nominal model, the modified Hansen scheme method did indeed produce better results, with a fit 28% in the horizontal direction and 44% higher in the vertical direction for poor a priori knowledge of the plant dynamics.The PEM method furthermore used considerably longer time to identify a model than the subspace identification method used to determine the Youla deviation system.This result might be due to how the PEM method is restricted to the model structure determined while the modified Hansen scheme is not.The identified systems using the Youla deviation system, however, also have some disadvantages.Any physical understanding related to model parameters of the gas bearing is lost when using the modified Hansen scheme.Furthermore, the methods produced models of higher order than the nominal model and the plant identified using the PEM for identification.This tendency of order increase comes naturally from the construction of the plant when the Youla deviation system has been identified.Unnecessary high-order model solutions are thus wellknown consequences when using the modified Hansen scheme.Model reduction has therefore been conducted which was able to produce acceptable results.It has thus been shown that the method has its clear advantage over PEM identification when model expert knowledge is lacking.

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

Figure 1 .
Figure 1.(a) Picture of the experimental test rig used for conducting active fault detection.The different parts of the test rig are as follows: a, vertical piezo actuator; b, horizontal piezo actuators; c, flexible shaft; d, vertical displacement sensor; e, disc and f, horizontal displacement sensor.(b) Cross-section illustration of the actuator.The further out the pin is moved the more air is able to flow in which increase the pressure.The position of the shaft is thus controlled by changing the flow of air through the inlet valve.

Figure 2 .
Figure 2. Closed loop system with the plant represented as a nominal part ( J G ) and the unknown Youla deviation system S.

Figure 3 .
Figure 3. Representation of a YJKB parametrised controller for generation of signals for identification.

Figure 4 .
Figure 4. Plot of the R 2 fit in the horizontal direction (top plot) and vertical direction (bottom plot) using the PEM for identification and subspace identification of the Youla system.

Figure 5 .
Figure 5. Boxplot of the R 2 fit given direct identification of the plant and identification of the Youla deviation system.The left plot is for the horizontal direction, whereas the right plot is for the vertical direction.

Figure 6 .
Figure 6.Comparison of a verification signal (blue) with the predicted output given the nominal model (red) and each of the three identified models.The nominal model G was found to have a fit of 3.188% in the horizontal direction and 2.999% in the vertical direction.The PEM (yellow) was found to have a fit of 48.91% in the horizontal direction and 39.47% in the vertical direction.Identification of the Youla deviation system (purple) gave a fit of 76.76% in the horizontal direction and 83.58% in the vertical direction.

Figure 7 .Figure 8 .
Figure 7.Comparison of a verification signal (blue) with the predicted output given the nominal model (red) and each of the three identified models.The nominal model G was found to have a fit of 51.36% in the horizontal direction and 54.13% in the vertical direction.The PEM (yellow) was found to have a fit of 62.87% in the horizontal direction and 75.78% in the vertical direction.Identification of the Youla deviation system (purple) gave a fit of 61.04% in the horizontal direction and 69.19% in the vertical direction.

Figure 9 .
Figure 9.Comparison of a verification signal (blue) with the predicted output given the nominal model (red) and each of the three identified models.The nominal model G was found to have a fit of 81.93% in the horizontal direction and 76.12% in the vertical direction.The PEM (yellow) was found to have a fit of 83.9% in the horizontal direction and 84.63% in the vertical direction.Identification of the Youla deviation system (purple) gave a fit of 83.53% in the horizontal direction and 81.18% in the vertical direction.

Figure 10 .
Figure 10.Zoom-in on the comparison of a verification signal (blue) with the predicted output given the nominal model (red) and each of the two identified models, PEM (yellow) and Youla deviation system (purple).

Table 1 .
Model fit using each of the two identification methods.