**Contents**



## **About the Special Issue Editors**

**Nicola Pio Belfiore**, Full Professor, IEEE Member, teaches Applied Mechanics and Functional Design at the University of Roma Tre, Italy. After being awarded his Ph.D. degree, completed at "Sapienza" in cooperation with the University of Maryland of College Park, he won three international awards in recognition of his research. He has served as Honorary Professor of Obuda University, Hungary, since 2008. The author of three textbooks, three patents, and around one hundred scientific papers, Belfiore has also served as coordinator of numerous scientific projects, both national and European. In 2013, he was Director of the 2nd Level Vocational Master in Energy Conversion Efficiency and Renewable Energy. His current research interests include dynamics, functional design, MEMS, robotics, tribology, and kinematics. In October 2017, he moved from Sapienza to Roma Tre University, where, in December 2020, has was appointed Head of the Degree Programs of Mechanical Engineering, which includes the BSc and MSc in Mechanical Engineering, MSc in Aeronautical Engineering, and BSc and MSc in Marine and Ocean Engineering.

**Pietro Ursi** is General Surgery University Fixed-Term Researcher at the Department of General and Specialized Surgery "Paride Stefanini" of Sapienza University in Rome, where he teaches courses concerning general surgery at MS level. His activities have been funded through research grants for surgical specialties and organ transplantation. After his achievement of the Italian Laurea Degree in Medicine and Surgery in 2000, Pietro Ursi obtained his Ph.D. in "New Technologies in Surgery" from the Consortium of Marche Polytechnic University and Sapienza University in 2009. He has been working as a specialist in digestive system surgery and digestive endoscopy at Sapienza University and has also been a member of the Italian Medical Association since 2001. He carries out outpatient activities, first aid ward assistance, and surgical operating room activities in small, medium, and large surgeries in addition to serving as a first aid doctor on mobile units.

**Andrea Scorza**, Researcher, teaches Industrial Measurements at the University of Roma Tre, Italy. He received his Ph.D. degree in Mechanical Measurements for Engineering from the University of Padova, Padua, Italy, in 2004. He has served as Assistant Professor of Measurements and Clinical Engineering at the Department of Engineering, University of Roma Tre, Roma, Italy, since his appointment in 2012. His current research interests include mechanical and thermal measurement systems and instrumentation, design and testing of biomedical instrumentation, and experimental mechanics applied in biomedical fields.

## **Preface to "Engineering for Surgery"**

The interaction between engineering and surgery serves as a source of progress for both sectors. The high standards of surgical operations, in terms demand of accuracy, reliability, and miniaturization, present a challenge for engineers, while some new achievements in engineering, in turn, have inspired numerous surgeons to improve the efficacy and success rate of their operations.

This is an intrinsically vast and interdisciplinary subject and, therefore, the present Special Issue offers only a little sample of the immense variety of applications, some of which have been successfully applied or are still under development, while we have been offered a preview of others thanks to the fantasy of creative science fiction writers.

The purpose of this Special Issue is, therefore, to stimulate the interest of engineers and surgeons, who will benefit from mutual advantages gained from their cooperation. We have been pleased to receive a number of contributions, and we sincerely appreciate all contributions, although the acceptance rate for this Special Issue was around 50%.

Finally, we gratefully acknowledge the entire staff of MDPI for supporting and trusting our work. Our service for Applied Sciences has made us more aware of the activities related to the communication of scientific results. A huge thank you goes to Luca Shao, who encouraged and supported all of us to develop the current Special Issue from the beginning of last year until today.

> **Nicola Pio Belfiore, Pietro Ursi, Andrea Scorza** *Special Issue Editors*

### *Editorial* **Engineering-Aided Inventive Surgery**

#### **Nicola Pio Belfiore 1,\*, Andrea Scorza <sup>1</sup> and Pietro Ursi <sup>2</sup>**


Received: 1 June 2020; Accepted: 3 June 2020; Published: 7 June 2020

**Abstract:** This Editorial presents a new Special Issue dedicated to some old and new interdisciplinary areas of cooperation between engineering and surgery. The first two sections offer some food for thought, in terms of a brief introductory and general review of the past, present, future and visionary perspectives of the synergy between engineering and surgery. The last section presents a very short and reasoned review of the contributions that have been included in the present Special Issue. Given the vastness of the topic that this Special Issue deals with, we hope that our effort may have offered a stimulus, albeit small, to the development of cooperation between engineering and surgery.

**Keywords:** engineering; surgery; interdisciplinarity

#### **1. Introduction**

For centuries, humankind has been dreaming about how to save lives and pursue immortality, and for this reason medicine and surgery have been always fundamental topics. Fiction and science fiction illustrate clearly the vast variety of expectations that people have from progress in these two important fields. An example of how human vision pushes forward the most secret ambitions is described in *Frankenstein* Mary Shelley's 1818 novel [1], wherein the doctor main character puts together pieces of dead bodies to build a new body and makes him alive with electricity. In the American epic space-opera *Star Wars* (Lucasfilm, 1977) there are also many examples of how medicine is expected to be in a far future: the prosthetic hand that replaces Luke's lost one, in a perfectly equivalent manner; the hibernation techniques; the robotic obstetricians; the fully automated orthopedics apparatus that allows a total replacement of the lower limbs; and Darth Vader's portable automatic respirator. Another interesting example of futuristic surgery has been suggested in *Fantastic Voyage* (20th Century Fox, 1966), which nowadays receives an exaggerated number of mentions at conferences by authors presenting their work in micro or nanosurgery. According to the plot, after an incredible miniaturization process, a submarine about the size of a microbe flows in a patient's ducts to remove a blood clot in his brain. Additionally, the task of producing a correct and fast diagnosis is every doctor's secret dream. In the science fictional series *Star trek* (Desilu Productions, 1966), chief medical officer Dr. McCoy obtains an immediate and detailed diagnosis simply moving a small sensor back and forth over the patient.

Unfortunately, we are far enough away from these goals, and therefore hard, maybe impossible, work remains for us to do along the road ahead. One way to start our endeavor consists of straightening the cooperation between medicine end engineering, because any progress in any technical apparatus that gives enhancements in surgery is based on both the fields of application. Furthermore, a great amount of creativity and interdisciplinary approach is needed to enhance the developments of new tools, which we could refer to as *inventive engineering for surgery* or *engineering-aided inventive surgery*.

The need for cooperation is intrinsically related to the fact that doctors *know what to do*, while engineers *know how to help to make it real* or may even suggest new facilities that open new scenarios for unheard-of operations. At first, progress may have been developed thanks to the creativity of surgeons who had a bit of technical know-how or who dared to experiment with new technical stuff. Again, making reference to another American television drama series, in *The Knick* (AMBEG Screen Products, 2014), Dr. Thackery, chief surgeon at the Knickerbocker Hospital, in case of emergencies, enhances surgical procedures by using the technical equipment in very creative arrangements.

After all these fictional examples, it is a pleasure to mention the exciting and enlightening survey of real cases in surgery recently written by van de Laar [2], where 28 *historical* operations have been described in terms so clear that even an engineer can understand. Among the conclusive remarks, it is very agreeable to assess that, for the moment, *there is as yet no question of computers completely taking over the tasks of human doctors*. However, the described operations show how the cooperation between engineering and medicine has been or could have been important to complete the task with success.

#### **2. From Early Tools to the State of the Art**

One way to explain how the cooperation between Technology and Surgery works consists in interpreting it as a customer-provider relationship where engineering offers new technological developments to the surgery' demand. In order to appreciate how strong the surgery–engineering relationship is, let us consider the following classical and fundamental topics in engineering and some of their representative applications:


Many other branches of engineering are relevant too. All of these capabilities are also necessary to develop most of the new frontiers of surgery, such as


They can be applied to general, lung, gynecologic, head and neck, heart, neuro-spine, vascular, and urological surgery.

#### *Classical Fields of Application*

The above-mentioned relationship between engineering and surgery cannot be described in terms of a simple offer-demand interchange. In fact, this cooperation is quite complex and multifaceted a liaison: any possible development in surgery could be supported by a proper collaboration with the engineering counterparts, while almost any new development in engineering could be successfully applied to improve surgical operations. Both alliances need a strong and very integrated partnership and an enduring team work. For all these reasons, the original call for papers from this Special Issue has been open to the following general topics:


While more specific topics have been also solicited:


The next section offers a very short and reasoned review of the contributions that have been included in the present Special Issue. Given the vastness of the topic that this Special Issue deals with, we hope that our effort may have offered a stimulus, albeit small, to the development of cooperation between engineering and surgery.

#### **3. About the Present Issue**

Nowadays the science of engineering may support surgery in different ways and through synergies that were hardly conceivable until a few years ago, thanks to the recent advances in applied sciences and technology. In particular, engineering contributions may range from pre-operative assessment to post operative care, and from computer aided-surgery to hardware development for performance improvement of consolidated treatments or novel surgical approaches, such as management and characterization of surgical tools and instrumentation.

Therefore, in this special issue some stimulating contributions are proposed for their valuable applications into the pre-operative field, focusing on modern simulation methods and 3D imaging tools for surgical planning, prediction methodologies [22,23] and data acquisition by means of novel wearable devices [24].

Anyway, the operating field in the context of synergies between engineering and surgery provides most of the advanced and promising solutions. In this regard, further interesting applications are described in the issue: from the control of MRI-compatible robots [25] and the guidance of surgical needles [26], to the use of very complex image analysis methods for surgical tool characterization [27] and the development of novel devices with high functional performances [28] and better ergonomic design for laparoscopic applications [29].

Moreover, solutions and innovations become very forward-thinking when engineering science is challenged with the requirements of minimally invasive surgery, as reported, for example, in the paper [30] where cataract surgery is combined with high frequency deep sclerotomy (HFDS). One more article [31] concerns the novel laser assisted in situ Keratomileusis (LASIK) applications for vision correction in myopia and presbyopia diseases.

Finally, a significant part of modern surgery relies on pioneering efforts to bring about advances in nano and microengineering to the lab and surgical activities. In this issue, an example of extreme miniaturization of the tools used in surgery has been provided by two papers, one concerning the development of a new piezo MEMS tweezer for soft materials characterization [32] and one describing some reasonable progress of MEMS for operations in a surgical scenario [33].

**Funding:** This research received no external funding.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**

1. Wollstonecraft Shelley, M. *Frankenstein, or, the Modern Prometheus: The 1818 Text*; Oxford University Press: Oxford, MI, USA; New York, NY, USA, 1998.


c 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## *Review* **3D Soft-Tissue Prediction Methodologies for Orthognathic Surgery—A Literature Review**

#### **Elena Carlotta Olivetti 1, Sara Nicotera 1, Federica Marcolin 1,\*, Enrico Vezzetti 1, Jacqueline P. A. Sotong 2, Emanuele Zavattero <sup>2</sup> and Guglielmo Ramieri <sup>2</sup>**


Received: 2 July 2019; Accepted: 22 October 2019; Published: 26 October 2019

**Abstract:** Three-dimensional technologies have had a wide diffusion in several fields of application throughout the last decades; medicine is no exception and the interest in their introduction in clinical applications has grown with the refinement of such technologies. We focus on the application of 3D methodologies in maxillofacial surgery, where they can give concrete support in surgical planning and in the prediction of involuntary facial soft-tissue changes after planned bony repositioning. The purpose of this literature review is to offer a panorama of the existing prediction methods and software with a comparison of their reliability and to propose a series of still pending issues. Various software are available for surgical planning and for the prediction of tissue displacements, but their reliability is still an unknown variable in respect of the accuracy needed by surgeons. Maxilim, Dolphin and other common planning software provide a realistic result, but with some inaccuracies in specific areas of the face; it also is not totally clear how the prediction is obtained by the software and what is the theoretical model they are based on.

**Keywords:** orthognathic surgery; 3D face analysis; surgical planning; soft tissue prediction; prediction methods

#### **1. Introduction**

Patients presenting dentofacial deformities are commonly subject to combined orthodontic and surgical treatment such as Le Fort I osteotomy (LFI), bilateral sagittal split osteotomy (BSSO), intraoral ramus vertical osteotomy (IVRO), sagittal split ramus osteotomy (SSRO), bimaxillary surgery and genioplasty. These interventions have been commonly planned with two-dimensional methodologies; today, the last challenge is the three-dimensional (3D) surgery planning. The refining of 3D graphics and imaging tools gives the chance to explore surgical planning and prediction of the effects of different clinical approaches; these techniques are based on images acquired with computed tomography (CT), cone bean computed tomography (CBCT) and multi-slice computed tomography (MSCT), which provide volumetric images of facial anatomical structure. Additionally, the surface of the face can be scanned and mapped to underline the effect of changes in facial appearance using 3D laser technology [1] that gives a great contribution to surgeons to decide on the type of surgeries as well as on the magnitude and direction of surgical movements to correct facial dysmorphology. Moreover, a deep interest in the prediction of soft-tissue response to hard tissue movements has been growing and two-dimensional conventional methodologies seem to be insufficient for this aim as they do not take into account the third dimension.

Indeed, the possibility to know the soft-tissue response to surgical operations helps surgeons to plan surgical movements and it gives surgeons more information about the need of orthodontic decompensation. Furthermore, the purpose of these interventions is not only to correct the facial dysmorphology from a functional point of view, but also to obtain an aesthetic enhancement of patients' facial aspect. Therefore, an accurate treatment planning is very important to obtain a good aesthetic and occlusal result [2]. For the aforementioned reasons, having a preview of the soft-tissue arrangement is extremely important. In future, this opportunity could bring to have a surgical planning methodology able to best match patients' aesthetic expectations, additionally to the functional corrective aspect.

Several methods have been considered to forecast soft tissue responses; the most common are the mass-spring model (MSM), finite element model (FEM) and mass tensor model (MTM). Most of the software packages currently adopted in clinical practice are based on these models. Generally, those softwares seem to reach an acceptable overall accuracy, but with inaccuracies for specific areas of the face, for example around the lips.

Studies have been conducted to develop 2D prevision models based on the ratio between facial and bony movements, particularly focusing on face sub regions [3]; by comparing predictive models constructed on this ratio and patients' post-operative conditions, it is stated that traditional approaches show limits in forecasting the outcome of large and complex movements. Moreover, a large number of variables must be selected [4].

Although there are several studies on the soft tissue changes after maxillary osteotomies, few of them report a systematic analysis.

The target of this systematic review is to gather information on the existing prediction methods and software of soft tissue displacements after dysmorphism corrective surgery, and to draw conclusions on the accuracy and reliability of these software in the preview of surgical outcome. Additionally, some problems related to the level of accuracy needed by surgeons, to the predictive imprecision reported in the studies and to the magnitude of the acceptable error are presented.

This work is structured as follows. The section Material and Methods describes the methods that have been adopted for the literature review. The section Results includes the clinical details of patients involved and the details of the articles considered; the articles have been divided in subsections according to the software used by the authors. Finally, the section Discussion and Conclusions discusses the current prediction methods and concludes the work.

#### **2. Materials and Methods**

The research of this systematic review is based on the Population Intervention Control Outcome Study design (PICOS) format (Table 1). PubMed and Scopus are the databases adopted for our research. The considered keywords are: orthognathic surgery, facial dysmorphism, surgical planning, 3D, 3D face analysis, soft tissue, BSSO, bilateral sagittal split osteotomy, Le Fort I, prediction methods.

The research has been set on different combinations of keywords; at first, we focused generically on the orthognathic surgery, then the type of interventions has been specified: BSSO, IVRO, SSRO and LFI. Finally, the research has been limited to soft tissue simulation and prediction.

All the articles found were assessed by three authors and classified in: prospective study (PS), retrospective study (RS), case series study (CS).

Single cases reports have been excluded from our analysis because they have been considered not clinically significant.

Articles published before 2000 have not been considered.

The following data have been recorded for each eligible study: first author and year of publication, journal, study design, sample size, gender, mean age, diagnosis, imaging technology, typologies of surgery, software used for the prediction, soft and hard tissue landmarks considered, time interval from surgery to post-surgical imaging, results and conclusions of the authors.


**Table 1.** PICOS (Population Intervention Control Outcome Study) criteria for the systematic review.

#### **3. Results**

Different responses have been obtained varying the insertion order of the keywords. After having combined the outcomes and removed the duplicates, the remaining articles have been evaluated on the basis of their relevance to the topic. In the end, 24 articles to be deepened and a number of interesting articles to be referenced for a better comprehension of the topic have been selected.

#### *3.1. Clinincal Details*

A total of 24 articles have been compared in this work; 12 are retrospective studies, two are case studies and five are prospective studies. Remaining articles have not declared categorization. The sample of patients involved in each study vary from seven to 100 subjects. The soft tissue prediction has been assessed for Le Fort I osteotomy, BSSO, BSSRO and genioplasty in the correction of different types of facial dysmorphism.

Tables 2 and 3 report a brief summary of the referenced papers we focused on. In particular, the demographic details of the patients are summarized in Table 2. A brief overview of methodologies and results of each article is reported in Table 3.

CBCT has been used for the assessment of soft tissues changes in twelve of the twenty-four studies, with the addition of 3D photographs for one article, while cephalometric radiographs have been used in the remaining eight works, with the addition of CT and MSCT for two of them.

The timing of post-surgical imaging has been stated in all articles.

#### *3.2. Prediction Methodologies*

Several approaches have been considered to make a mathematical three-dimensional prediction of soft tissue changes after orthognathic surgery. MSM, FEM [5,6] and MTM are the most common. These have been developed into software packages which are currently used in clinical practice. The functioning of these software is generally acceptable if we consider the creation of a plausible facial outcome, but this prediction does not necessarily match with the real outcome. Moreover, studies show that the prediction accuracy decreases for specific facial areas, especially around the lip, and with the increasing complexity of the surgery (larger bony repositioning often results in a higher inaccuracy of the prediction). The studies presented in this literature review evaluate soft-tissue predictions obtained with different software packages.



**Table3.**Imagingtechnology,softwareusedforprediction,numberoflandmarksforsofttissueandhardtissue,timeintervalfromsurgerytopost-operative


**Table 3.** *Cont.*


**Table 3.** *Cont.*

#### 3.2.1. Works Based on Traditional (2D) Methodologies

Since a certain number of 3D software currently available for clinical practice are based on algorithm developed to work on 2D and then adapted to 3D representation of data, a section dedicated to traditional soft-tissue prediction is introduced. Moreover, since many surgeons currently rely more on traditional techniques than on three-dimensional methods, it is reasonable to focus on the limitations of this technology; 2D analysis based on lateral cephalograms is subject to errors caused by the impossibility to accurately locate points on both hard and soft tissue relying only on 2D anatomic structures [7]. Analysis performed on cephalograms reported that reasonable prediction can be made on soft-tissue landmarks positions, but with low accuracy concerning the lips and particularly in the vertical prediction of displacements [8]. Moreover, the inaccuracy increases with a higher level of complexity of intervention, which leads to errors both in hard- and soft-tissue simulation [9].

In this section, we particularly focused on study assessing the accuracy of Dolphin Imaging, but we also reported studies that took into consideration other software.

Peterman et al. worked on the quantification of the accuracy of Dolphin VTO (visual treatment objective) in predicting soft tissue changes in patients with a class III deformity and on the validation of its efficacy. In their study, authors considered fourteen patients receiving comprehensive orthodontic treatment and orthognathic surgery including both maxillary advancement and/or mandibular setback. Cephalometric tracing and analysis were performed with Dolphin Imaging software; pre- and post-operative and traced cephalometrics were superimposed using the cranial base as reference. The maxillary movements were recorded at the anterior nasal spine (ANS) and A point, while the mandibular changes were recorded at B point and pogonion (Pg) in both x and y axis. Subsequently, the pre-treatment profile pictures were superimposed to digitally-traced soft-tissue landmarks of pre-treatment cephalometric radiograph to initiate software VTO simulation (Figure 1). Finally, a prediction profile photograph generated by the software was compared to actual post-operative patient's profile photograph; coordinates of nine soft-tissue landmarks were used to quantify on each axis the differences between the predicted and the actual of patients. The results reported in this study are consistent with the prevision error calculated in previous works; Dolphin Imaging software was found to have a much larger standard deviation in the Y direction than in the X, with various degrees of accuracy in both directions for all soft tissue landmarks. The accuracy was 79% and 61% for X-axis and Y-axis respectively, with an acceptable error set at 2.00 mm. Furthermore, for these surgical cases, the least accurate prediction was the one concerning the lower lip region. In conclusion, authors affirmed that VTO program can be useful for rough movement previsions during surgical planning, but surgeons cannot rely on it for precise surgical movements (measurement range lower than 1.00 mm) [10].

**Figure 1.** Super-imposition of pre- and post- and traced cephalometric tracing showing advancement of the maxilla and setback of the mandible [10].

In the study of Pektas et al., the treatment planning for each patient was constituted on the base of clinical and cephalometric evaluations, and pre-operative study models. All patients were treated with pre- and post-surgical fixed orthodontic appliances. Scanned images of cephalometric radiographs were processed using Dolphin Imaging software; the Ricketts and Steiner analysis was selected to start the digitization. To compare actual outcome and software simulation, firstly the superimposition of pre- and post-operative cephalograms was performed on the sella-nasion plane registered at sella. Landmarks located on structures not involved in the surgery were directly transposed from pre-operative to post-operative cephalogram to refine the superimposition process. The objective of an accurate superimposition of pre- and post-operative cephalograms was to construct a spreadsheet for landmarks movements to quantify changes of each point after operation; the position of points was defined in form of x and y. In accordance with the spreadsheet, a simulation of the treatment was produced, and the prediction tracings were obtained from actual post-operative changes amount. Actual post treatment tracings and prediction tracings were superimposed and differences have been measured again according to an x, y coordinate system.

The mean differences between the predicted and the actual result were less than 1 mm in four of seven soft tissue measurements; predictions turn out to be more accurate for the sagittal plane than for the vertical one. The largest errors were reported for the lower lip region, in accordance with other studies considered by authors and investigating other software [11].

Ana de Lourdes sá de Líra et al. studied the accuracy of digital prediction with respect to the actual post-operative outcome with Dolphin program. Patients with class II malocclusion were involved. Eighty Caucasian subjects treated with combined surgical and orthodontic treatment were considered. Cephalometric radiographs of pre- and post-operative were digitized; a time of at least one year from surgical to post-operative was required to rule out any effects of post-operative swelling. Digital tracings were performed with Dolphin imaging software and digital lateral cephalometric radiographs were traced at the same time to minimize the error variance; before tracing, Dolphin program was used to carry out the compensation for the effect of radiographic magnification. An x y coordinates system was constructed with the x axes corresponding to the horizontal Frankfort plane and the y-axes passing through Nasion and perpendicular to x-axes. On x-axes, vertical changes have been measured, while on y-axes authors evaluated horizontal changes. Dalberg formula was used in evaluating the reproducibility of measurements; with regard to the method error, it did not exceed 0.37◦ for angular measurements and 0.29 mm for linear measurements. Results of the statistical analysis were based on mean values found in each studied group. The analysis of the authors underlined that there were no significant differences between soft-tissue values at pre-operative (t1), post-operative (t2) and predicted (t3) stages. In both groups, surgeries have been more extensive than planned; facial convexity and the distance between the lips and the cranial base presented similar values between t2 and t3. The conclusion of the authors was that statistical differences in the considered measurements did not invalidate surgical prediction performed with the software, because forecasted changes were sufficiently close to of the surgical outcomes [12].

Lira et al. followed the same approach of the previous article on two groups of patients; the former underwent Le Fort I without segmental surgery, while on the latter, in addition to LFI surgery, mandibular setback was performed with an SSRO. What came out from this study was that 12 months after the operation surgical reductions in mandibular length and angle were substantially greater than indicated by predictive cephalometric tracings. The analysis of dental and skeletal mandibular changes in anteroposterior and vertical directions revealed that the effect on the profile was greater than those of the maxilla; this is probably due to the fact that the mandible is a moving structure and the upper lip leans upon the lower lip. Maxilla advancements were performed for both groups according to the planning provided by Dolphin imaging software; in group two, mandibular setback was statistically overestimated by the software planning. Nevertheless, authors concluded that statistically the differences in the considered measurements do not invalidate the surgical prediction generated with the software, because the predicted result was close to the surgical outcome. Moreover, they noticed that concerning the mandible, dentoskeletal measurements in horizontal and vertical directions showed a greater correlation to the profile than the maxillary measurements [13].

Magro-Filho et al. worked on a subjective comparison of the soft-tissue surgical simulation of two software: Dentofacial Planner Plus (Dentofacial Software, Toronto, Ontario, Canada) and Dolphin Imaging (version 9.0). Surgeries considered in the study involved retro-positioning of the mandible and advancement of the maxilla, in which linear movements were at least 4 mm for one bone segment or in the sum of both mandibular and maxillary movements, without, or least possible, postsurgical orthodontic movements. Profile predictions were made six months after the surgery; real movements, obtained from the superimposition of the cephalometries of pre- and post-operative of each subject, were used as values of the prediction. Soft-tissue images of pre-operative, real predictive images of Dentofacial Planner Plus (DFP) and Dolphin Imaging (DI) and post-operative images were compared (Figure 2). With Photoshop program (version 6.0, Adobe, San Jose, Calif), pre- and post-operative profile images were digitized and standardized; then, images were imported in DFP and DI. Movements of maxilla and mandible were made using the real values extracted from the comparison between pre- and post-operative lateral cephalograms of the surgical displacements along the vertical and horizontal directions for molars and incisors of the maxilla and for incisors of the mandible. The tip of the nose, the nasolabial angle, upper and lower lips, menton region, the base of the mandible and the complete profile were used to compare simulated and real images. A group of orthodontists, maxillofacial surgeons and general dentists evaluated the similarity among the programs simulated images and the real post-operative images with a scale ranging from "very similar" to "different". For the comparison of soft-tissue cephalometric points, to compare the two programs and to judge the criteria of similarity, chi-square test was used. From their analysis, the authors reported that for nasal tip, chin area and mandibular base DI gave better results, while for nasolabial angle, upper lip and lower lip DFP produced a more accurate prediction; authors also reported that there was no difference in the evaluation of the complete profile. Moreover, in terms of time, working with DI was longer than working with DFP. In conclusion, it was stated that in treating subjects affected by dentofacial deformities, it would be better if orthodontists and oral maxillofacial surgeons were to base on their clinical experience and to use software just as a co adjuvant. Considering the methodology of the study, there were not many differences in obtaining a two-dimensional prediction of patient's profile [14].

**Figure 2.** Profile images of a patient with Class III malocclusion, treated with orthognathic surgery: (**A**) pre-treatment; (**B**) simulation from Dentofacial Planner Plus (DFP); (**C**) simulation from Dolphin Imaging (DI); (**D**) actual post-operative photograph [14].

Other software are available for the prediction of soft tissue changes. OrthoForecast is a data-based prediction software; its database contains information of 400 patients (100 pts with asymmetry, 100 pts with skeletal class II jaw relationship and 200 pts with skeletal class III jaw relationship). All patients underwent SSRO or IVRO, with or without LFI osteotomy (respectively, 352 and 48 cases). It was evaluated by Abe et al.; the study included 15 patients with facial asymmetry, 15 with skeletal class II jaw relationship and 15 with skeletal class III jaw relationship. Patients underwent LFI osteotomy, SSRO and IVRO. Twenty-four landmarks were chosen and digitized. Some angles were also compared between the following lines: line 1, the line connecting the bilateral mandibular angles; line 3, the line connecting the bilateral eyespots; line 2, the line connecting the bilateral oral angles. In the asymmetry group, they measured angle 1, identified by line 1 and line 3 and angle 2, identified by line 2 and line 3 on the frontal facial photograph to assess the asymmetry of the face. In the class II and class III groups, the value of angle 3, identified by the angle formed by lines connecting eyespot, upper lip, and soft tissue B-point on the lateral facial photograph was measured. Angle 3 was the reference for evaluating the convexity of the face, which was originally defined as soft-tissue convexity (Figure 3). The distances measured between landmarks on the actual and on the predicted images were analyzed, and the similarity between real and software-generated images was assessed by 39 evaluators. The mean difference measured between the positions of landmarks is less than 3.4 mm and less than 1.0◦ in the evaluated angles. Moreover, more than half of the evaluators stated that, in all groups, the predicted images were very 'similar' or 'similar' to the actual situation of the patients in the post-operative condition. Less than 6% of the evaluators ranked the predicted images as 'different' from the actual outcome. In conclusion, OrthoForecast is considered by the authors to be a software with high levels of accuracy, reliability and usefulness [15].

**Figure 3.** Construction of feature points and lines (**A**) 1, Angle between line 1 (connecting the bilateral mandibular angles) and line 3 (connecting the bilateral eyespots); 2, angle between line 2 (connecting the bilateral oral angles) and line 3; (**B**) 3, angle formed with lines connecting the eyespot point, upper lip point, and soft tissue B-point [15].

Donatsky et al. evaluated another prediction software, TIOPS. This is a computerized, cephalometric, orthognathic surgical planning system previously used in studies on hard tissue stability and accuracy. As other soft tissue prediction systems, TIOPS is based on predefined ratios of hard to soft tissue movements; its present algorithms partly rely on cephalometric observation of the post-operative situation. The study included 52 patients. Clinical photographs, study models mounted on an articulator (SAM) and standardized lateral cephalometric radiographs of the pre-operative situation were performed. Standardized lateral cephalometric radiographs were performed 5–6 weeks after surgery. The mean accuracy of the planned and predicted results of both hard and soft tissue varied from 0.0 mm to 0.5 mm from one cephalometric reference point to another. In the locus of cephalometric reference points, where it was shown that there were significant differences between planned/predicted hard and soft tissue changing in terms of position, these significant inaccuracies were moderately small and varied from 0.2 mm to 1.1 mm, with the exception of the horizontal position of the lower lip. However, the variability of the previewed hard and soft tissue individual result was evaluated to be relatively high. The study demonstrated moderately high predictability of the immediate post-surgical hard and soft tissue outcome. However, due to the significant variability of

simulated individual results, it would be better to be cautious in presenting the planned and predicted hard and soft-tissue results to the individual patient in the pre-operative [16].

Finally, Chai Hui Koh et al. evaluated the accuracy of soft tissue predictions generated by the CASSOS (Computer-Assisted Simulation System for Orthognathic Surgery) software in Chinese skeletal class III patients underwent bimaxillary surgery. The digitization of pre-surgical and post-surgical lateral cephalograms of 35 patients was performed with the CASSOS program. A simulation of the surgery was performed on the pre-surgical tracing. An analysis was performed on thirty-two linear measurements on the cephalograms superimposition to assess the differences in the soft tissue profile between the post treatment results and the predicted result. It showed differences on 16 linear measurements with the most prediction errors on the upper and lower lip having a mean difference relatively small, with a greatest mean difference of 2 mm in the vertical position of stomion inferius. The authors concluded that CASSOS 2001 produced a clinically meaningful forecast of soft-tissue profile changes following bimaxillary surgery [17].

The same CASSOS software was evaluated by Jones et al. for 33 patients affected by class III skeletal deformities and underwent maxillary advancement (17 patients) or bimaxillary surgery (16 patients) [18]. The post-operative cephalograms were used to determine surgical bony movements to produce the CASSOS profile simulation. Linear differences between the predicted profile and the real outcome were measured at 12 soft tissue landmarks. The authors' conclusion was that the profile predictions obtained with CASSOS could be considered useful for both type of surgeries, even if substantial variations were found. As previously [17], the most inaccuracies were found in the lip region.

#### 3.2.2. Works Based on Maxilim Software (3D)

In this section a collection of studies is presented considering patients' soft tissue prediction obtained with Maxilim software (Medicim—Medical Image Computing, Mechelen, Belgium). This software, based on the mass tensor model algorithm (MTM), allows surgeons to determine bony movements and to see the effect of the procedure [19]. MTM has been introduced by Cotin et al.; the geometry of the anatomical structure is discretized into a tetrahedral mesh in which the displacement vector at a generic internal point is defined as a function of vertices displacement vectors. The elastic force is written as a function of these same four vectors. The displacement for a set of mass points is set fixed and other model points will move as consequence of elastic forces due to the displacements of fixed points; finally, the new rest position of the free points is computed by integrating the Newtonian motion equation [20].

Shafi et al. [21] investigated the accuracy of Maxilim using cone-beam computed tomography (CBCT) scans in pre- and post-surgery phases (6–12 months after surgery) for 13 patients subjected to Le Fort I surgery. A 3D mesh was generated from skeletal movements (predicted model). Subsequently, the mesh generated from the post-operative of the patient and the predicted model were compared to evaluate the accuracy of the prediction. The soft tissue was divided in different areas: nose, right and left nares, right and left paranasal regions upper and lower lip and chin. The absolute distance was calculated between the meshes for each region. For each facial region, the absolute distance was calculated between the meshes.

In almost all the facial regions, the accuracy was significantly less than 3.00 mm (3.00 mm is clinically acceptable). The upper lip area was the exception; in fact, for this region, the accuracy was greater than 3.00 mm. In all cases, Maxilim produced an over prediction of the new position of the upper lip. A possible reason could be a non-linear response of the upper lip to some hard-tissue movements and, consequently, the modeling algorithm used (MTM) hardly previewed this response.

The conclusion of the authors was that the 3D soft tissue prediction for Le Fort I advancements produced by Maxilim was in general clinically satisfactory, but it was associated with marked errors around the region of the upper lip.

Liebregts et al. evaluated the accuracy of Maxilim for 3D simulation of soft tissue changes after bimaxillary osteotomy. The 3D rendered pre- and post-operative scans were matched. Segmented maxilla and mandible were aligned to the post-operative position. In order to calculate the error between the simulation and the actual post-operative condition, authors used 3D distance maps and cephalometric analysis. Concerning the facial profile, the mean absolute error between the 3D simulation and the actual post-operative facial profile was 0.81 ± 0.22 mm for the face as a whole. In this study, the accuracies of the simulation (average absolute error ≤2 mm) for the whole face and for the upper lip, lower lip and chin sub regions were 100%, 93%, 90% and 95%, respectively.

Authors affirmed that the MTM based soft tissue simulation is an accurate model for the prediction of soft tissue changes following bimaxillary surgery. The magnitude of skeletal movements influences the accuracy of the prediction; moreover, the age of the patient and the use of V–Y closure affect the precision of the predicted model. Low predictability on the upper and lower lip regions is registered once again [22].

The same authors evaluated Maxilim performances for soft-tissue simulation in 100 patients underwent BSSO for mandibular advancement [23]. As in the previous case, the accuracy of the simulation was assessed with two methods, a 3D cephalometric analysis and a 3D distance map for the entire face and for specific regions of interest. Their analysis showed that for the entire face the mean absolute 90th percentile error was less than 2 mm (clinically acceptable); the least accuracy resulted in the region of the lower lip, while the most accurate prediction involved the sub nasal region. A possible explanation suggested by authors for labial inaccuracies was the difficulty to replicate labial position during different image acquisition.

The conclusion drawn by authors on the basis of this analysis was that the soft-tissue prediction produced by Maxilim software was clinically acceptable both for the whole face and for restricted regions, but the limitation was the fact that patients considered in the study underwent a single type of surgery instead of combined operations.

In the analysis by Mundluru et al., 13 non-syndromatic adults with a midline deviation of the chin point not less than 2.0 mm underwent Le Fort I or bimaxillary osteotomies (BSSO) to correct facial asymmetry. The accuracy of an innovative concept for the soft-tissue changes prediction was evaluated. The segmentation of the 3D model was performed to identify the following regions: upper lip, lower lip, chin area, right and left paranasal regions, nose and right and left cheeks. To evaluate the prediction reliability, the surface distances between predicted and actual post-operative positions were measured. Particularly, mean (signed) and mean absolute distances were calculated on 3D meshes for each region. Through a directional analysis, the accuracy of the prediction of soft-tissue changes was evaluated. The results showed that the distances between the predicted and the actual post-operative soft tissue models were less than 2.0 mm in all regions. In details, a general tendency to the under prediction was found in the area of the cheek and of the chin (medio-laterally). An over-prediction of changes was found at the inferior border of the mandible (bilaterally) [24].

De Riu, et al. realized a retrospective study with the purpose of giving an objective quantification of accuracy of the Maxilim virtual planning method, through the comparison between planned and actual movements in jaw osteotomy. They started from the assumption that a simple superimposition of simulation and cephalometry results did not allow to define a correlation between positional error and surgical movements. Indeed, a slight positional error could be insignificant talking about large movements, but could be unacceptable dealing with small displacements. A 3D planning was performed, and virtual frontal and lateral cephalometries of pre-operative and simulated surgery were extracted; those were compared to define the predicted movements of the jaw in reference to cranial bones. All surgeries have been planned by the same surgeons using Maxilim software, with digital intermediate splints to guide osteotomies. Cephalometric analysis has been performed with Dolphin. To evaluate the accuracy, the mean linear difference between planned and actual movements was considered. Measurement of the differences between planned and actual movements were deemed accurate. As reported by authors, most of the clinical studies validating the virtual surgical planning

set to 2 mm the criterion for the prediction success, with a success rate slightly lower than 100%. In this study, some significant differences with other works were found; first, differences between planned and achieved anterior facial heights were found (*p* = 0.033). This discrepancy was relevant (*p* = 0.042) in patients that did not undergo genioplasty, while it was not significant in patients who underwent this surgery (0.235). This error is probably due to the approximation of the soft tissue model, which does not allow the management of the vertical dimension. Presumably, genioplasty compensate the error on facial eight. The second error reported by authors was in the differences between actual and planned measurements for SNA (angle between Sella, Nasion and A point) and for SNB (angle between Sella, Nasion and B point), with *p* = 0.008 and *p* = 0.006 respectively. In the opinion of the authors, virtual planning inaccuracies were principally caused by the difficulty of simulate the soft tissue changes. Moreover, they concluded that virtual planning could not relieve surgeons of the necessity of monitoring jaws movements intraoperatively and of a real time supervision of planned and actual outcome, despite a high level of accuracy for most of the analyzed parameters [25].

The double blind prospective study performed by Van Hemelen et al. had the aim of providing a comparison between a traditional planning method and a 3D planning performed with Maxilim, both for hard and soft tissue [26]. In their analysis, authors considered 66 patients affected by class II or class III angle malocclusion underwent bimaxillary osteotomy (46 patients), BSSO (17) patients and LeFort I osteotomy (3 patients). Genioplasty was performed on 28 patients. For the 2D planning, clinical facial examination, lateral and frontal cephalograms were taken; 17 and 16 cephalometric landmarks were considered for 3D and for 2D respectively, to perform the validation of both 2D and 3D planning. The traditional planning was performed with Onyx Ceph Version 3.1.111 and the analyses of accuracy was performed by measuring projections on x and y axis of distances between the 16 selected points on the planned and on the actual post-operative. In the case of the 3D planning, the landmarks were defined in three dimensions and the planned and post-operative data were aligned. Distances between the 17 cephalometric points defined on planned and actual models were evaluated in depth, height and distance in the sagittal plane.

As a result of the analysis, authors reported that, for hard tissue planning, the mean differences occurring between actual and planned outcome evaluated at cephalometric points were 1.71 mm and 1.42 mm in the horizontal direction and 1.69 mm and 1.44 mm in the vertical direction, for 2D and 3D respectively. Concerning soft tissue displacements, the mean values computed in the horizontal direction were 2.29 mm for traditional planning and 1.48 for Maxilim prediction; in the vertical direction, authors noted mean values of 2.07 mm and 1.46 mm for 2D and 3D techniques respectively. These evaluations led authors to conclude that there was a statistically significant difference in the soft-tissue prediction between the two considered methods (chi-squared test and independent t-test were used to determine statistical significance, with *p* < 0.005); particularly, the 3D approach seemed to be more predictive than the traditional one. No statistically significant differences were reported for hard tissue planning.

The conclusion reached by authors was that, concerning the soft tissue prediction, the 3D planning approach here performed with Maxilim software led to a more accurate planning than the traditional method, even with disadvantages such as the cost of the software and of the CBCT scans. Moreover, a higher learning time must be considered when planning with the 3D method. In the authors' opinion, 3D planning will not replace the traditional prediction if those negative aspects are not addressed.

The study by Nadjmi et al. focused on the accuracy of soft-tissue profile changes simulated by Maxilim with respect to the accuracy obtained with Dolphin. Particularly, the goal of the comparison was to assess if the 3D prevision was more accurate than the 2D one. A "natural head distance" was defined as distance between suprasternal notch and soft tissue pogonion to replicate patients' head position. On lateral cephalogram, soft tissue and hard tissue landmarks (15 and 25 points respectively) were located; soft tissue predictions were generated with Dolphin and Maxilim. The two predictive results were superimposed with patient's post-operative profile photograph and differences in terms of linear measurements of landmarks in x and y directions were measured. All comparison steps

were performed by the same operator. The movements of bony structures in the simulation must be the same as performed by the surgeon, to correctly compare the two predictions. Since this study made a comparison between 2D and 3D simulations, even if Maxilim gave in output 3D models, only profile (2D) images were considered. From the analysis of the results, the authors found that both Dolphin and Maxilim predictions were accurate, but the simulation obtained with Maxilim made possible a quantification of volumetric changes (lips aspect) and a prevision of changes in the transverse plane. Even if the study focused on the comparison of accuracies (and for this reason only the lateral prediction was considered) between the two software, authors stated that 3D simulation would be helpful in complicated surgeries thanks to the third-dimensional information [27].

#### 3.2.3. Works Based on Other 3D Software

The work of Resnick et al. investigates the accuracy of Dolphin 3D soft tissue prediction using seven patients who had a single-segment Le Fort I osteotomy. Dolphin 3D Imaging (Dolphin Imaging and Management Solutions, Chatsworth, CA, USA) is based on a landmark photographic morphing algorithm developed for 2D prediction and adapted to 3D. The user must locate 79 landmarks, 47 on the hard tissue and 37 on the soft tissue, on the CT volume; the system generates adaptable curves between the landmarks, similarly to the tracing of a lateral cephalometry.

On the pre-operative segments LF I osteotomy was simulated by Dolphin, then hard and soft tissue landmarks were assigned, morphing curves were adjusted and finally prediction image (Tp) and post-operative image (T1) were aligned and registered for the measurement of the differences. For each patient, the error was calculated as the difference occurring between Tp and T1 at 14 points, six on the midline and eight laterally and at the nasolabial angle (Figure 4). Twelve of these points are standard anthropometric or cephalometric landmarks. Two new points were introduced for the study: lateral ala (LA), obtained from the intersection of the line tangent to endocanthion (EN) and the one tangent to subalare (SBAL), and maxillary buttress (MB), traced at the intersection of lines tangent to exocanthion (EX) and SBAL (Figure 4A).

**Figure 4.** Points for measurement by Resnick et al. (**A**) Frontal view. The two points derived for this study are shown: lateral ala point (LA), defined as the intersection of lines tangent to endocanthion (EN) and subnasale (SBAL), and maxillary buttress (MB) point, and defined as the intersection of lines tangent to exocanthion (EX) and SBAL. (**B**) Lateral view showing all points [28].

In closing, their analysis showed that the capability to predict the soft-tissue changes in the three dimensions after LF I surgery using Dolphin software had some limitations. For linear changes, its accuracy was acceptable, while it was not acceptable for lateral points of the face. Concerning the midline, changes at nasal base were more subject to prediction errors [28].

Bianchi et al. [29] and Marchetti et al. [30] worked on SurgiCase\_CMF using CBCT and MSCT respectively of 10 patients with facial deformities. In the study presented by Bianchi et al., patients underwent LFI, BSSO and genioplasty. The pre- and post-operative CT images were aligned on the top of each other. Firstly, a registration point was made, with the location of two corresponding landmarks on pre- and post-operative conditions. After the manual registration, automatic surface registration was carried out using an iterative closest point algorithm. A crucial point in surfaces registration was to indicate those that did not change after the surgery; authors used the eyes area and the forehead as fixed regions. Finally, the post-operative soft tissue surface was compared with the software virtual simulation. A comparison algorithm was used to measure the distance of every triangle corner of the post-operative surface against the preoperative one. The percentage of error was <2.0 mm for the 86.80% of patients, even if there were important errors in the top lip and chin regions. The results in the group of patients studied with CBCT (which had a reliability of 86.8%) stated that the use of SurgiCase\_CMF 1.2 software combined with CBCT data allowed to achieve a realistic preview with low X-ray exposure [29].

Marchetti et al. divided virtual surgical planning in four steps (Figure 5): (1) CT data reconstruction; (2) generation of 3D models of hard and soft tissues; (3) various virtual surgical planning and simulation mode; (4) different pre-operative previews of the soft tissues. Skeletal model surgical planning and simulation were performed on the base of the 3D CT; the soft tissue model was made through a physical model. According to clinical options, the software generated a set of simulations and models of the soft tissues; in order to prevent temporomandibular functional problems, an orthodontist assessed the pre-operative plans. From the comparison of simulation and CT surgical outcomes, it was found that the prediction of soft tissue situation had a reliability greater than 91%, with a percentage of error lower than 2 mm for 76%–99%. According to the obtained results, SurgiCase\_CMF was able to give a realistic and accurate preview of the face of the subjects undergoing surgery, but there were significant errors in the lip and chin regions, as in the previous studies [30].

**Figure 5.** Virtual surgical planning proposed by Marchetti et al. [30].

3dMDvultus is a software package in which the rendering function is based on Mass Spring Model. MSM assumes a discretization of a deformable object into n mass points and a set of m connections between each n mass point. In a tetrahedral discretization, for each mesh node a point is allocated, and it is defined a linear spring for all the edges of the mesh. This linear spring follows Hook's law. This model has the advantage of being simple and computationally efficient, but the disadvantage is that the elastic behavior of the model is determined by the spring constant; the value of this constant is an approximation and it has no true bio-mechanical relevance [19].

Khambay et al. used 3dMDvultus to predict the soft tissue changes after surgery. The work included 10 patients who underwent Le Fort I osteotomy. Ten landmarks were used to measure the distances between the pre- and post- operative meshes. Moreover, the percentage of mesh points minor or equal to 2.0 mm were calculated for the full face and for specific anatomical regions. The results demonstrated that the percentage of mesh with errors minor or equal to 2.0 mm for the full face was 85.2%–94.5% and 31.3%–100% for anatomical regions. The range of root mean square (RMS) error was from 2.49 mm to 0.94 mm. The most of mean linear distances computed between the surfaces was equal or less than 0.8 mm, but it increased for the mean absolute distances. From this analysis, it was stated that the use of specific anatomical regions was more significant in clinical practice than the use of full face [31].

As in the previous article, the software considered by Terzic et al. is 3dMDvultus. Three-dimensional photographs of patients' head in natural position (lips and muscles at rest, open eyes, neutral facial expression) of pre- and post-operative were taken; CT and/or CBCT images were imported to the software platform. Pre-operative 3D photograph was fused with pre-operative CT/CBCT images. In the same manner, post-operative bony skull data was fused with pre-operatives by matching bone areas not involved in surgeries. According to actual post-surgical CT, 3D bone segments were reproduced in the pre-operative skull; the software rendering function (mass spring model) was activated. Osteotomy segments were moved to real post-operative position and the soft-tissue rendering generated the textured facial soft tissue prediction. After the fusion of 3D photographs of predicted and real post-operative soft tissue, a horizontal plane was positioned arbitrarily to divide the face in an upper part not involved in surgery, and a lower part including the overlap of simulated and real post-surgery result. Averaged distribution of absolute error showed more discrepancies between predicted and real post-operative outcome in the lower part of the face: in 29.8% of lower halves authors found errors exceeding 3.00 mm. This preliminary study concluded that the software platform had an insufficient accuracy in the forecast of 3D soft-tissue displacements [32].

Ullah et al. carried out a retrospective study in which they focused on the ability of 3dMDvultus to give a prediction of the face appearance of patients undergoing LFI maxillary advancement. The null hypothesis they started from was that the mean difference in absolute distance between the predicted surface generated by the software and the real 3D facial surface measured at eight anatomical regions did not exceed 3 mm. Pre- and post-operative CBCT for each patients were imported into 3dMDvultus and hard tissue and soft tissue were separately segmented and saved in STL (standard triangulation language) format. A template of actual surgical changes was produced using a CAD/CAM software (VRMesh, VirtualGrid, Seattle, WA, USA). Post-operative hard and soft tissue linked together were aligned to the anterior cranial base of pre-operative hard tissue of the same patient. Regions of maxilla and mandible of post-operative models were saved as STL file and imported into 3dMDvultus to generate a template of the actual maxillo-mandibular complex aligned in the same 3D space as the pre-operative image. Using the same software, a LFI osteotomy was simulated for the pre-operative hard tissue. The resulting soft-tissue prediction was exported as an STL file for the analysis. The predicted model and the real surgical result were imported as mesh in VRMesh; the differences between the two models were rendered as a color map, which showed the overall accuracy of the software. After the analysis on the overall surface, the accuracy was evaluated subdividing the mesh into anatomical regions: chin, lower lip, upper lip, nose, right and left nares, right and left paranasal areas. From the analysis of the results, the authors concluded that since the distances between software-generated surfaces and real facial surfaces were lower than 3 mm, 3dMDvultus 3D soft-tissue predictions were acceptable for clinical use in LFI osteotomy. MSM seemed to correctly predict the positions of lip and chin, while improvements should be made for nasal and paranasal regions [33].

In the prospective study by Holzinger et al., 16 patients with open bite dentofacial-dysmorphosis and underwent orthognathic surgery (surgery first) were analyzed. The surgery was planned using conventional sketches as a new software developed by the authors, SOTIRIOS planning software. A conventional pre-surgical planning was carried out using three sides X-rays and, in addition, a CT scan was performed before and 6 months after surgery. A constrained fitting approach using statistical shape modeling (SSM) technique to generate the patient-specific data was used in combination with the anatomical landmarks to compute the patient specific model. A quantitative comparison of the soft tissue prediction and post-operative data was made showing a mean error of 1.46 mm ± 1.53 mm. Authors concluded that the SOTIRIOS planning software was quite accurate and enabled the surgeon to predict the soft tissue outcome [34].

A comparison between Dolphin 3D, ProPlan CMF and an in-house probabilistic finite element (PFEM) [6] was carried out by Knoops et al. [35]. The retrospective study included seven patients who underwent LFI maxillary advancement. From CBCT taken preoperatively, 3D models of bone and soft tissue were constructed. On the base of the advancement and rotation movements measured on the post-operative images, a virtual Le Fort I osteotomy was performed and three different soft tissue predictions were realized with the three softwares. Firstly, the three predictions were compared with the pre-operative model to evaluate differences. After that, the three predictions were compared with the real post-operative outcome, to assess which was the most predictive model. As the patients involved underwent a LFI maxillary advancement, the areas of interest for the comparison with post-operative situation were the upper lip and the paranasal regions. From the outcomes analysis, it resulted that the prediction generated by Dolphin 3D was affected by a general under-prediction of paranasal displacements, while both ProPlan CMF and the PFEM model were characterized by an over-prediction of displacements involving cheilion regions. Average root mean square distances and average percentage of points less than 2 mm between real post-operative images and predictions were computed (Table 3); statistically significant differences resulted from the Friedman test. In particular, the post hoc Wilcoxon signed-rank test proved that RMS for PFEM and RMS for ProPlan had statistically significant lower values with respect to RMS for Dolphin 3D (*p* = 0.016). On the opposite, the differences between ProPlan and PFEM could not be considered statistically significant (*p* = 0.219).

The conclusions of the authors was that ProPlan CMF and PFEM gave soft-tissue predictions with significant accuracies, particularly when using the correct post-operative maxillary position; because of the limitations introduced by a landmark-based algorithm and a sparse architecture, Dolphin 3D resulted to be inaccurate when surgery involved large maxillary advancements, particularly on lateral points.

Nam et al. [36] tried to assess if a 3D virtual surgery could accurately predict the soft tissue outcome of 29 patients who underwent bimaxillary orthognatic surgery (LFI, BSSO and genioplasty) performed by the same surgeon. The predicted outcome was generated using the Simplant Pro program. All simulations were performed by the same operator. The produced simulation of the soft tissue result was superimposed on the post-operative of the patients, and 10 landmarks were designed and positioned on the two models; the simulation error was determined by measuring distances (x, y, z) between the same landmarks on predicted and real images. From this analysis, an accuracy of 52.8% resulted, disagreeing with previous study reporting accuracies around 80%. Most of the errors were found for landmarks of lower lips, mouth corner and chin area (pogonion and menton), with a general larger error in points located on the mandible. Authors suggested as a possible cause of this trend the complexity of the surgery, involving more movements in the mandible rather than in the maxilla; it could determine lower errors in the maxilla but higher inaccuracy in the mandible. Inaccuracies reported for 18 of the 30 soft-tissue landmarks were considered statistically significant (*p* < 0.05), with errors concentrated at mouth corners. Possible explanations addressed by the authors were the ethnic

differences between the data used for developing the software (Caucasian population) and the data used in this study (Asian population); moreover, it was not negligible that errors up to 2 mm could be due to the positioning of soft tissue landmarks. In conclusion, authors stated that 3D software in soft tissue prediction could have great potential, but the accuracy must be improved.

Even if it is out of the scope of the present work, it is important to state that not only the surgical planning plays a crucial role in the optimal process of care of the patients. It is not negligible that maxillofacial surgery, as well as other branches of surgeries, could lead to medical complications such as venous thromboembolism, which can have fatal consequences when moving to the lung (DVT/PE) [37]. Moreover, surgeons must manage patients with rare disorders of hemostasis, which need particular care in the perioperative phase. The works by Simurda et al. [38,39], and Ghadimi et al. [40] have been referred to in order to deepen understanding of this topic.

#### **4. Discussion and Conclusions**

In this systematic review we have summarized twenty-four studies focused on the evaluation of soft-tissue prediction software for maxillofacial surgery. The treated software are commonly used in clinical procedures to have a preview of the patients' outcome after corrective surgery; the most common methodologies on which these software are based are the mass spring model, the finite element model and the mass tensor model, each of them with its advantages and disadvantages. MSM presents an easy architecture and a low memory usage but it has no real biomechanical foundations; FEM is more relevant from a biomechanical point of view, but it has a high computational cost and a rather long simulation time; MTM tries to combine the advantages of MSM and of FEM [18]. In particular, the lip region is the most difficult to predict, since lips do not rely on bony structures; moreover, it is difficult to replicate the same lip position from one acquisition to another. Due to this particularity, lips are supposed to be particularly involved in involuntary displacements caused by surgeries performed on different facial areas. The difficulty in accurately predicting the appearance of this sub region is a crucial point in the development of a prediction method, since it is particularly related to aesthetic self-perception and satisfaction of the patient. With respect to the traditional algorithms based on hard- to soft-tissue ratio, methodologies involving biomechanical models of viscoelasticity of the soft tissue to mimic its elastic deformation seems to be able to give a more realistic simulation. Despite the number of computer programs dedicated to three-dimensional prediction and surgical planning, maxillofacial practitioners often prefer to rely on conventional two-dimensional methods. This fact could be explained with the overall results given by these methodologies, which seem to be only partially relevant and/or reliable in the practice of daily clinical activity. Concerning the works based on Maxilim, Shafi et al. found that this software makes a clinically satisfactory prediction of soft tissue changes, but with an over-prediction of the position of the upper lip [21]; according to Liebregts et al., its accuracy is influenced by the magnitude of maxillary and mandibular movements, by the age of the patient and by the usage of V-Y closure, and these authors further reported a low predictability of lower and upper lips [22]. Moreover, a general tendency was reported towards an under-prediction in the area around the cheek and chin and towards an over-prediction at the mandibular inferior border [24]. It was suggested that Maxilim inaccuracy in the prediction of the lip region might be due to the modelling algorithm, i.e., mass tensor, and might indicate a non-linear response of the upper lip [21]. The study by De Riu et al., although not directly focused on soft tissue prediction, reports that inaccuracies in virtual planning are primarily due to the soft-tissue virtual model [25].

With regard to Dolphin software, Resnick et al. reported that the ability to predict 3D soft tissue changes was limited, with an accuracy that may be acceptable for linear changes but not for lateral facial points [28]. Peterman et al. reported various degrees of accuracy at each soft-tissue landmark in the horizontal and vertical axes; lower lip prediction was the least accurate. A similar result of low prediction accuracy for the lower lip region is reported by Pektas et al. [11]. Conversely, two articles affirmed that the statistical differences between Dolphin soft tissue prediction and the real outcome of the patient did not invalidate the software prediction, as it was close to real surgical results [12,13]. Dolphin 3D imaging uses a landmark-based photographic morphing algorithm that was developed for two-dimensional prediction and then adapted to three-dimensional prediction; this could be an explanation for its imprecisions. Moreover, it requires plotting 79 landmarks on the CT volume (42 for the hard tissue and 37 for the soft tissue); the localization of such a large number of points seems to require too much time for a stable usage in clinical daily practice.

Among the other considered software, SurgiCase\_CMF was analyzed by Bianchi et al. [29] and Marchetti et al. [30]; they reported that the percentage of error was less than 2 mm for the majority of patients, but in both studies an important error in the areas of lips and chin was reported. Khambay et al., Terzic et al. and Ullah et al. used 3dMDvultus to predict soft tissue changes; in general, they reported that its use on specific anatomical regions was more meaningful than on the full face [31–33]. Moreover, they reported an insufficient prediction accuracy particularly for the lower part of the face with errors exceeding 3 mm [32], and for nasal and paranasal regions [33]. Surprisingly, prediction results seemed to be reliable for lips and chin.

The results reported for TIOPS showed significant inaccuracies in the predicted horizontal position of the lower lips and a relatively high variability of the predicted individual surgical outcome. The authors suggest care is needed when presenting the predicted image to patients before the operation [16,41]. The SOTIRIOS planning software and the CASSOS program were found to be quite accurate [34] and able to give a clinically useful prediction of the soft-tissue profile change [17]. OrthoForecast software resulted as being in general an accurate tool for soft-tissue prediction; indeed, more than half of the involved evaluators assessed the predicted images to be very similar or similar to the actual post-operative images [15].

As a result of this review, it seems that Dolphin, CASSOS and OrthoForecast are the most accurate software in soft-tissue prediction. However, their reliability is not sufficient when further surgery is carried out and the complexity of the operation increases; moreover, they provide only partial information on the prediction, due to the lack of the third dimension. From this viewpoint, all articles agree on the potentiality of 3D methodologies. Maxilim produced clinically meaningful results, but its inaccuracy in specific facial areas (particularly in the lip region) led to the conclusion that further steps are needed to give surgeons a tool for predicting facial movements (including involuntary movements) that significantly supports the planning of the intervention in terms of aesthetic success. The comparison between 2D and 3D methodologies shows that the first reach better results in terms of accuracy and surgeon satisfaction. Consequently, the latter are still not able to meet medical needs.

The analysis of these studies has led us to consider some issues that are still pending. Almost all of them report problems of reliability of the software or inaccuracy in specific areas; this shows that an ameliorative method would be necessary to improve prediction results. For this reason, it would be useful to clearly understand the way the existing software compute the tissue displacements. Moreover, it is important to notice that some of the considered software (Dolphin 3D) are not based purely on 3D methods, but they are an adaptation of a 2D algorithm to 3D data; the transfer of 2D data to a 3D representation may cause errors. This aspect affects both hard- and soft-tissue prediction, and it is more evident in soft tissue simulations, due to the complexity of its behavior (involuntary displacements). The tendency of surgeons to prefer conventional methods suggests that, despite the general positive opinion emerging from the studies, the reliability of automatic systems is not sufficient, since it has not significantly increased above the two-dimensional methods commonly used in clinical routine. Indeed, it seems to lack consistency between the results obtained with prediction software and practitioner satisfaction. As previously reported, if the statistical differences in the evaluated measurements do not invalidate surgical predictability software, which parameters can be used to evaluate the precision of the prediction? It may be that the evaluation of automatic output is qualitative more than quantitative. Surgeon opinion on prediction software is that, even if the soft-tissue simulation outputs are plausible with a possible result of the interventions, these do not correspond to the effective post-operative outcome; moreover, this discrepancy is even more evident when the complexity of surgery increases. Neither hard- to soft-tissue ratio-based algorithms nor those based on viscoelastic models seem to

meet surgeon demand. Another aspect is that currently soft tissue prediction seems to be focused on specific intervention effects more than on involuntary displacements; presently, these displacements affect, in an unpredictable manner, the facial outcome of patients. An accurate and reliable prevision method should consider these displacements to give surgeons an effective and helpful tool.

To improve prediction, a clear technical scientific evaluation approach should be established, which seems to be lacking in certain studies. Moreover, it would be better to clarify the level of precision of the predicted outcome needed by surgeons. These aspects are also an incentive for researchers to investigate the correlation between the extent of soft tissue changes with respect to STTs (soft-tissue thickness) and BMI (body mass index).

**Author Contributions:** Conceptualization, F.M. and E.Z.; methodology, J.P.A.S., S.N., F.M.; software, E.C.O. and S.N.; validation, S.N., E.C.O., and J.P.A.S.; formal analysis, F.M.; investigation, S.N.; resources, J.P.A.S. and E.Z.; data curation, E.C.O. and S.N.; writing—original draft preparation, S.N., E.C.O., and J.P.A.S.; writing—review and editing, F.M. and E.Z.; visualization, F.M.; supervision, F.M., E.Z., E.V., G.R.; project administration, E.V. and G.R.; funding acquisition, E.V. and G.R.

**Funding:** This research received no external funding.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## *Article* **The E**ff**ect of Patellar Tendon Release on the Characteristics of Patellofemoral Joint Squat Movement: A Simulation Analysis**

#### **Jianping Wang, Yongqiang Yang, Dong Guo, Shihua Wang, Long Fu and Yu Li \***

School of Mechanical and Power Engineering, Henan Polytechnic University, Jiaozuo 454000, China; wjp@hpu.edu.cn (J.W.); y17864208259@163.com (Y.Y.); 15538949519@163.com (D.G.);

211805020010@home.hpu.edu.cn (S.W.); fl1@gaokowl.com (L.F.)

**\*** Correspondence: liyu@hpu.edu.cn

Received: 7 August 2019; Accepted: 11 October 2019; Published: 14 October 2019

**Abstract:** Objectives: This paper studies the patellar tendon release's effect on the movement characteristics of the artificial patellofemoral joint squat to provide reference data for knee joint surgery. Methods: Firstly, the dynamic finite element model of the human knee joint under squatting was established. Secondly, in the above no-release models, the release of 30% of the attachment area at the upper end, the lower end, or both ends of the patellar tendon were conducted, respectively. Then the simulations of all above four models were conducted. Finally, the results of the simulation were compared and analyzed. Results: The simulation results show that, after releasing the patellar tendon (compared with the no-release simulation's results), the relative flexion, medial-lateral rotation, medial-lateral tilt, and superior-inferior shift of the patella relative to the femur increased; the medial-lateral shift and anterior-posterior shift of the patella relative to the femur decreased. Conclusion: In this paper, the maximum flexion angle of the patella increased after the patellar tendon being released (compared with the no-release model), which indicated that the mobility of knee joint was improved after the patellar tendon release. The simulation data in this paper can provide technical reference for total knee arthroplasty.

**Keywords:** total knee arthroplasty; patellar tendon; patellofemoral joint; squat movement; dynamic finite element analysis

#### **1. Introduction**

With the application of prostheses in total knee arthroplasty (TKA), post-operative complications concomitantly usually occur [1]. One of these complications, often being ignored, is patella baja [2–4]. At the same time, relevant studies [5–9] show that, after TKA, the patellar tendon suffers certain reductions. Other studies [6,10] indicate that the patellar tendon can be shortened, the position of the patella lowered, reducing the knee joint movement and variation. A slightly low patella does not, however, affect the function of the knee joint. Check points during surgery must include patella contracture syndrome [11] and true patella baja [10], which can be treated by release of scar operation, patellar tendon lengthening, and patellar ligament insertion. Patients with pseudo-patella baja and slight joint line up, and those that do not affect knee function, are treated with follow-up observation. However, excessive upward movement of the joint line can result in serious symptoms [12], such as the need for surgery to replace a thin pad [13] and/or patellar ligament release. This paper is for studying the patellar tendon release's effect on the movement characteristics of the artificial patellofemoral joint squat to provide reference data for knee joint surgery.

#### **2. Materials and Methods**

#### *2.1. The Establishment of a Finite Element Model*

In this paper, a healthy male volunteer with a height of 173 cm and a weight of 60 kg was scanned by medical instrument computed tomography (CT) in the range of 10 cm above and below the knee joint center. The correlative parameter was set to 120 kVp and 150 mA, with a scanning interval of 1 mm. The volunteer was scanned with magnetic resonance imaging (MRI). MR (GE Sigma HD1.5T, Boston, MA, USA) was used for proton density weighting (PDWI) and chemical shift fat suppression weighted scanning, and the image coordinates of each layer of different weighted phases were consistent. Flux density is 1.5 T, and slide thickness is 1 mm [14]. Based on the CT and MRI images of the healthy volunteer, a three-dimensional geometric anatomical model including all bone tissues and main soft tissues was established (as shown in Figure 1). A three-dimensional solid model of press fit condylar (PFC) femoral prosthesis (PFC sigma, DePuy orthopaedics, Warsaw, IN, USA.) commonly used in TKA was established. Then the total knee replacement operation was simulated, and the prosthesis was assembled on the three-dimensional geometric anatomical model of the knee joint (as shown in Figure 2a). The model of the knee joint was divided into finite element meshes after TKA. In the finite element model of the knee joint, most of the unit types were hexahedrons (c3d8r, an eight-node linear brick element with reduced integration), and very few were pentahedrons (C3D6, a six-node linear triangular prism element) and tetrahedrons (C3D4, a four-node linear tetrahedron element), with a total unit number of 31,767, and 37,031 nodes (as shown in Figure 2b). The volunteer was provided details of the study and signed an informed consent form.

**Figure 1.** (**a**) Contour tracing of bone from CT; (**b**) MRI image; (**c**) the geometric anatomy model of human knee.

**Figure 2.** (**a**) 3D geometric anatomical model of knee joint after TKA; (**b**) finite element model of knee joint after TKA.

Hypermesh software, version 11.0 (Altair Engineering Corp, MI, USA) was used to simulate the release at ends of the patellar tendon. The width of the patellar tendon at the upper end was 6.84 mm before release (Figure 3a) and 4.788 mm after release (Figure 3b). The height of the patellar tendon at the lower end was 8.09 mm before release (Figure 3c) and 5.633 mm after release (Figure 3d). Figure 3 shows the finite element models before and after release the 30% of the patellar tendon at both ends (clinically, one-third of the central patellar tendon is often transplanted to repair damaged ligaments [15–18]. Therefore, 30% of the patellar tendon was used as an example). Figure 3a represents the model before release of the patellar tendon at the upper end. Figure 3b represents the model released the patellar tendon at the upper end. Figure 3c represents the model before releasing the patellar tendon at the lower end. Figure 3d represents the model releasing the patellar tendon at the lower end. The establishment of the model was based on the model built by Wang Jian-ping, and the validity of the model was verified by in vitro cadaveric experiments [19].

**Figure 3.** The finite element model of patellar tendon: (**a**) the upper before release; (**b**) the upper released; (**c**) the lower before release; (**d**) the lower released.

#### *2.2. The Establishment of a Coordinate System for the Knee Joint's Motion*

The corresponding reference coordinate systems of femur, tibia and patella are established by the following methods: (1) The femur posterior condylar cortical bone was contour-fitted to two laps, with one line through the center of two circles as the X axis of the femur. Parallel with lines of force and through the centerline of the femoral medial condyle represent the Z axis of the femur. Regarding the line, this passes through the center of the femoral medial condyle and perpendicular to the X-axis and Z-axis, as the Y-axis. (2) The intersection point of the Z-axis and the plane of the tibial joint line serves as the origin of the tibial coordinate system, and the femoral coordinate system moves to this point to obtain the tibial coordinate system. (3) The center of the patellar anterior surface moves inward 10 mm, then the center serves as the origin of the patellar coordinate system. The coordinate system of the femur was relocated to the origin, and the patellar coordinate system can be obtained [20]. The bone coordinate systems are shown in Figure 4.

**Figure 4.** Coordinate system of bone. (**a**) femur; (**b**) tibia; (**c**) patella.

#### *2.3. Loading Conditions and the Setting of Material Properties*

For this study, 400 N [19,21] of force was applied to the quadriceps, and the direction of the force was parallel with the femoral shaft, pointing to the starting point of the quadriceps. At the same time, half of the body's weight (0.5BW) was applied along the force line of the knee joint, that is, a force of 300 N perpendicular to the ground was applied to the center of the femoral head. The boundary conditions were defined: the reference point of the ankle corresponding to the rotation center was fully fixed with six degrees of freedom. The material properties of bone tissue and femoral prosthesis were defined as isotropic and linear elastic, respectively. The material properties of polymer polyethylene pad were defined as non-linear elastic-plastic deformability [22,23]. The material properties of soft tissue were defined as nonlinear elastic [24]. The friction coefficient of polymer polyethylene and cobalt chromium molybdenum material was defined as 0.04 [22]. The penalty function with weighted factor was adopted [25,26]. For the femur, polymer polyethylene pad, tibial plateau, patellar prosthesis, and other soft tissues, seven surface contact pairs were defined in the finite element model after TKA. Moreover, nine surface contact pairs were defined in the finite element model of normal human.

#### **3. Results**

The four different knee joint finite element models were imported into the finite element analysis ABAQUS software, version 6.10 (Dassault SIMULIA Corp., Paris, France). Within 0–135 degrees, the data of motion characteristics at different flexion degrees were calculated and analyzed, as shown in Figure 5.

**Figure 5.** A sample of simulation results after TKA. Flexion degrees: (**a**) 0◦; (**b**) 30◦; (**c**) 90◦; (**d**) 135◦.

#### *3.1. Medial-Lateral Shift and Flexion of Patella Along the Inner and Outer Axes*

Figure 6a shows the medial-lateral shift data of the patellofemoral joint within 0–135 degrees of knee flexion. As the figure shows, from 0 to 30 degrees of flexion, the patella relative to the femur shifted outward. After releasing the patellar tendon, the patella shifted outward greatly. The lateral shift of the patellar tendon released at the lower end was larger than that of release at the upper end. The lateral shift was the largest when both ends were released simultaneously. From 40 degrees to 135 degrees of flexion, the patella relative to the femur shifted inward. After releasing the patellar tendon, the patella medial shift was smaller than that of no-release. The medial shift of the patellar tendon released at different ends was basically balanced. After the upper end of the patellar tendon was released, the maximum of the patella medial shift was 2.62 mm at 110 degrees of flexion, which was 11% lower than that of 2.94 mm at 100 degrees of flexion for the no-release model. After the lower end of the patellar tendon was released, the maximum medial shift of the patella was 2.63 mm at 110 degrees of flexion, and decreased 11% by compared with the no-release model. After both ends of the patellar tendon were released, the maximum medial shift of the patella was 2.61 mm at 110 degrees of flexion, and decreased 11% by compared with the no-release model.

Figure 6b shows the changes of patella flexion before and after the release. At the same flexion degree of the knee, the release of the patellar tendon at the upper end and lower end was compared with the no-release model, the flexion of the patella relative to femur was slightly larger. The patella flexion after both ends release (compared with the no-release model), in the same flexion degree of the knee, showed a smaller degree of patella flexion (within 130–135 degrees of flexion, being slightly larger than the no-release model). At 135 degrees of knee flexion, for the three releasing models: the upper end release, lower end release, and both ends release models, the patellar flexion reached maximums of 92, 92.73, and 93.17 degrees, respectively. Compared with the no-release model of 89 degrees, this shows an increase of 3%, 4%, and 5%, respectively.

**Figure 6.** Patellofemoral movement data at different flexion degrees before and after the patellar tendon was released. (**a**) Medial-lateral shift; (**b**) flexion of patella relative to femur; (**c**) superior-inferior shift; (**d**) medial-lateral rotation; (**e**) anterior-posterior shift; (**f**) medial-lateral tilt.

#### *3.2. Superior-Inferior Shift and Medial-Lateral Rotation of the Patella Along Its Upper and Lower Axes*

As Figure 6c shows, during the knee joint's squatting movement, the patella relative to the femur shifted downward. The patella inferior shift of no-release model was slightly smaller than that of the upper end release. For the lower end release simulation, the patella inferior shift was slightly smaller within early 0–20 degrees of flexion and slightly larger after the 20-degree flexion than that of the no-release model. Within all flexion degrees, the patella inferior shift of the both ends release was slightly larger than that of no-release model, other than 20-degree flexion. After releasing the upper end of the patellar tendon, the maximum inferior shift of the patella was 43.80 mm at 110-degree flexion, which increased 6% compared with 41.42 mm at 100-degree flexion of the no-release model. For the lower end release model, the maximum inferior shift of the patella was 43.49 mm at 100-degree flexion, which increased 5% compared with the no-release model. After releasing the patellar tendon at both ends, the maximum inferior shift of the patellar was 43.43 mm at 100-degree flexion, which increased 5% compared with the no-release model.

Figure 6d shows the data of medial-lateral rotation. The data of the patella relative to the femur showed lateral rotation first and then medial rotation. Compared to the no-release model and the upper end model, the knee flexion degree of the initial patella medial rotation was greater for the lower end model and both-ends model. Meanwhile, it is shown that the patella medial-rotation for all three release models as larger than that of the no-release model. At 60 degrees of knee flexion, the patella medial-rotation reached the maximum of 5.41 degrees, 5.70 degrees, and 5.72 degrees, respectively, for the upper-end, lower-end and both-ends models, which increased 1%, 7%, and 7%, respectively, compared with the 5.34 degree value of the no-release model.

#### *3.3. Anterior-Posterior Shift and Medial-Lateral Tilt of Patella Along the Anteroposterior Axis*

Figure 6e shows the patella translated backward relative to the femur. During the knee joint's squatting movement, the patella relative to the femur translated backward. Within all flexion degrees, the patella backward translations of all three release models were smaller than that of no-release model, other than the 30-degree flexion, and the patella posterior shift was the least when both ends were released. At 135 degrees of knee flexion, the patella posterior shift reached the maximum values of 61.82 mm, 62.57 mm, and 61.07 mm, respectively, for the upper-end, lower-end, and both-ends models, which decreased 5%, 4%, and 6% compared with the value 64.97 mm of the no-release model.

Figure 6f shows the data of medial-lateral tilt. Before and after patellar tendon release, within all flexion degrees, the patella showed an inward movement pattern relative to the femur other than 30-degree flexion. At the same degree of knee flexion, after the release of the patella tendon, the medial tilt angle of the patella relative to the femur was larger than that of no release. After 80 degrees of knee flexion, when both ends of the patellar tendon were released at the same time, the patella medial tilt angle was the largest. At 135 degrees of knee flexion, the patella medial tilt angle reached the maximum values of 5.51 degrees, 5.01 degrees, and 5.92 degrees, respectively, for the upper-end, lower-end, and both-ends models, which increased 21%, 10%, and 30%, respectively, compared with the 4.55 degree value of the no-release model.

#### **4. Discussion**

This paper investigates the patellar tendon release's effect on the movement characteristics of the artificial patellofemoral joint squat. The simulation results show that, after the patellar tendon being released, respectively, at the upper end, the lower end, and both ends, being compared with the no-release simulation results, the relative flexion, medial-lateral rotation, medial-lateral tilt, and superior-inferior shift of the patella relative to the femur increased, and the medial-lateral shift and anterior-posterior shift of the patella relative to the femur decreased. Different release models have different effects on the squatting motion characteristics of the patellofemoral joint. In this study, the data of patellofemoral joint squatting after TKA of the no-release model was analyzed and compared with those of all three release models and other related research results of the natural knee joint.

#### *4.1. Medial-Lateral Shift and Flexion of the Patella*

For the medial-lateral shift of the patella relative to the femur the maximum patella medial-shift of all three release models were balanced, and the maximum medial-shift was decreased compared with that of the no-release model (the maximum reduction was about 11% for the both-ends model). The patella shifted outward first and then inward relative to the femur, which is similar to the related studies [27,28].

For the flexion movement of the patella, relative to the femur, previous studies [27,28] showed that the flexion of patella changes linearly with the flexion of the patellofemoral joint. It is consistent with the results of this study. In this study, the maximum flexion angle of the patella increased after the patellar tendon being released (compared with the no-release model, the maximum of the both-ends model increased by 5%), which indicated that the mobility of the knee joint was improved after the patellar tendon release.

#### *4.2. Superior-Inferior Shift and Medial-Lateral Rotation of Patella*

For the superior-inferior shift of the patella, relative to the femur, after the patellar tendon was released in different ends, the patella shifted downward relative to the femur. In the high flexion, the patella was moved up slightly, and the maximum inferior shift (about 43.5 mm at 110 degrees) was increased 5% compared with that of no-release model. The maximum increase occurred in the upper-end release model and the minimum one occurred at both-ends release model. Related studies [27,29] also showed that the patella, with the knee flexion, shifted downward relative to the femur, and the maximum flexion angle of the knee was 110 degrees in their studies (the maximum inferior shift was about 50 mm), which was basically similar to this paper (in this paper, the maximum flexion angle of the knee was 135 degrees, and the maximum inferior shift was 43.8 mm).

For the medial-lateral rotation of the patella, relative to the femur, after the patellar tendon was released at different ends, the motion pattern of the patella relative to the femur was lateral rotation first and then medial rotation, which is consistent with that of the no-release model. For all the above three release models, the maximum medial rotation of the patella increased; among them, the increase of the upper end release was the minimum (1%), and the increase of the both-ends release were the maximum (7%). Heegaard et al. [30] showed that within 0–110 degrees of knee flexion, the patellar rotation pattern changed from medial to lateral rotation, and then lateral rotation within 110–135 degrees of knee flexion. This phenomenon indicates that the patellar motion changes during the high flexion, which might be related to the way of femoral fixation and passive loading.

#### *4.3. Anterior-Posterior Shift and Medial-Lateral Tilt of the Patella*

For the anterior-posterior shift of the patella, relative to the femur, the results of Azmy et al.'s [31] in vitro experiments showed that with the flexion of the knee joint, the patella shifted backward; after 30 degree knee flexion, posterior movement rate of the patella increased, and the maximum posterior movement of the patella reached 32 mm. The results of Baldwin et al. [29], based on the experimental dynamic finite element model, also show that the patella continues to shift backward with the knee flexion, and the maximum posterior movement of the patella is about 30 mm. The motion analysis results of patella anterior-posterior shift in this paper are consistent with those of the above studies; with knee flexion, patella continued to shift backward. Meanwhile, the analysis results of this paper gave the data of the posterior shift of the patella during high knee flexion. At 135 degree knee flexion, the posterior shift of the patella after the upper end, the lower end, and both ends of the patellar tendon were released to reach the maximum values of 61.82 mm, 62.57 mm, and 61.07 mm, respectively, which was 5%, 4%, and 6% lower than the no-release values of 64.97 mm.

For the medial-lateral tilt of the patella relative to the femur with the flexion of the knee joint, the patella relative to the femur exhibits an inward motion pattern before and after release. At the same degree of flexion, the medial tilt angle of the patella relative to the femur was larger than that of the no-release model, the increase of the both ends release was the maximum (30%), and the increase of the lower end release was the minimum (10%). The relevant findings [27,29,31,32] also showed that the patella, with the flexion of knee, tilted inward relative to the femur, which is consistent with the results of this study. At the same time, there are different research results for the medial-lateral tilt motion of the patella [33]. The results of Amis et al. [33] showed that the patella continued to tilt outward with knee flexion, which may be related to the loading method that the quadriceps were divided into six bundles. There are differences in related studies, which may be related to the fact that the patella was less constrained by the femoral articular surface in the direction of patella medial-lateral tilt.

Several limitations of the present study should be noted. In this study, only 30% of the patella was released, and others were not analyzed. The finite element model of the human knee joint after TKA was established in this paper, based on the consideration of all bone tissue and the main soft tissue of the knee; although the transverse patellar ligaments, which play an important role in the movement of the patella, are added, the effect of the joint capsule, hamstring, or gastrocnemius on knee joint activity was not considered. In the subsequent analysis of the knee joint, these tissues could be taken into account to make the model more in line with human physiology.

The anterior cruciate ligament (ACL) has the functions of restricting the forward movement of the tibia and controlling the internal and external valgus of the knee [34–36]. The ACL is extremely vulnerable during exercise and it is difficult to repair itself. If not treated in time, the patient's daily activities will be seriously affected [37]. The posterior cruciate ligament (PCL), aside from providing restraint to posterior tibial translation, is reported to have a considerable role in providing rotational stability to the knee [38]. Injuries to the PCL are far less common than injuries to other knee structures, such as the ACL or the menisci [39]. In addition, it has been found that the clinical outcomes of PCL reconstruction are less satisfactory and less predictable relative to those of ACL reconstruction [39]. The medial patellofemoral ligament (MPFL) is one of the most important restraints to lateral displacement of the patella from 0◦ to 30◦ of flexion, providing up to 60% of lateral patellofemoral stability [40]. In the last decade, the understanding of the pathoanatomic mechanisms involved in lateral patellar dislocation (LPD) has evolved tremendously and reconstruction of the medial patellofemoral ligament (MPFL), as the most important passive restraint against LPD, has evolved to become an established operative procedure [41]. This paper mainly investigates the patellar tendon release's effect on the movement characteristics of the artificial patellofemoral joint squat. However, different knee ligaments play different roles. Based on the models used in this study, we will further study the ligaments' reconstruction and other related issues.

#### **5. Conclusions**

This paper investigates the patellar tendon release's effect on the movement characteristics of the artificial patellofemoral joint squat. The post-TKA patellofemoral joint's motion characteristics in different patellar tendon release pattern were obtained. Compared with the no-release model, the flexion, medial-lateral rotation, medial-lateral tilt, and superior-inferior shift of the patella relative to the femur increased; and the medial-lateral shift and anterior-posterior shift of the patella relative to the femur decreased. In this paper, the maximum flexion angle of the patella increased after the patellar tendon being released (compared with the no-release model), which indicated that the mobility of knee joint was improved after the patellar tendon release. Different release models have different effects on the squatting motion characteristics of patellofemoral joint. In this paper, the simulation results can be used to moderately release the ends of the patellar tendon in the TKA surgery, in order to improve the motility of the patella relative to the femur. The results of this paper can provide a reference for the study of the pathology and rehabilitation of the knee joint, as well as related operations.

**Author Contributions:** Conceptualization: J.W.; data curation: L.F. and Y.L.; investigation: J.W. and L.F.; methodology: J.W. and Y.L.; resources: J.W.; software: J.W., Y.Y., and D.G.; supervision: J.W. and Y.L.; validation: Y.Y., D.G., and S.W.; writing—original draft: J.W.; writing—review and editing: Y.Y., D.G., S.W. and Y.L.

**Funding:** This research was funded by National Natural Science Foundation of China, grant number 31370999.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## *Article* **User Friendliness of a Wearable Visual Behavior Monitor for Cataract and Refractive Surgery**

#### **Bojan Pajic 1,2,3,4,\*, Pavel Zakharov 5, Brigitte Pajic-Eggspuehler <sup>1</sup> and Zeljka Cvejic <sup>2</sup>**


Received: 23 December 2019; Accepted: 17 March 2020; Published: 24 March 2020

**Abstract:** A prospective feasibility study was conducted to determine whether a new wearable device, the Visual Behavior Monitor (VBM), was easy to use and did not present any difficulties with the daily activities of patients. Patients for cataract surgery and refractive lens exchange were randomly selected and screened for inclusion in the study. A total of 129 patients were included in the study as part of a multicenter study. All measurements were performed before surgery. Upon inclusion, patients were trained to wear the device, instructed to wear it for a minimum of 36 h, and were scheduled to return in one week. The VBM measures the distance at which patients' visual activities are performed, the level of illumination, and head translational and rotational movements along the three axes. On the follow-up visit, patients completed a questionnaire about their experience in wearing the device. All patients underwent standard diagnostic testing, with their cataract grade determined by the Lens Opacities Classification System (LOCS) classification. Results indicate that 87% of patients felt comfortable using the wearable device while 8% of patients responded as not feeling comfortable (5% of patients did not respond to the question). In addition, 91% of patients found it easy to attach the wearable to the magnetic clip while 4% of patients did not find it easy, and 5% of patients did not respond. Overall, patients found the device easy to use, with most reporting that the device was not intrusive.

**Keywords:** cataract surgery; refractive surgery; visual behavior monitor

#### **1. Introduction**

Cataract surgery has become the most common operation worldwide. In the United States alone, 3.6 million cataract surgeries are performed each year. This corresponds to a cataract surgery rate (CSR) of 11,000 per 1 million inhabitants. In economically well-developed countries such as the United States, Europe, Australia, or Japan, the CSR is between 4000 and 11,000 [1,2]. Yet, refractive outcomes, and patient satisfaction, following cataract surgery remain a key challenge for ophthalmologists today, particularly when it involves a presbyopia-correcting intraocular lens (IOL) [3].

In an era where we strive to improve on our outcomes by using increasingly sophisticated biometry and IOL power calculations, it is important not to lose sight of the impact of residual or increased refractive errors that may not be apparent using the spherical equivalent or viewing the components of the refractive error in isolation on the patient's vision. We should use the tools that facilitate surgeons to assess and manage the patient's refractive error in its entirety.

Approaching the intended post-operative spherical equivalent, however, often does not achieve spectacle independence. Uncorrected residual spherocylindrical refractive errors appear to have far greater adverse effects on unaided visual acuity than may be evident using a spherical equivalent or the individual sphere and cylinder. If the correction of presbyopia is desired, there are three options available for IOLs in the course of cataract surgery. One is the installation of a monovision where the dominant eye is set for distance and the non-dominant eye for near. A difference of up to 2 D is accepted by the patient. However, before proceeding to surgery, it is essential to try wearing contact lenses to see if the correction is tolerated [4,5]. Most patients reported the surgery met their expectations for decreased dependence on spectacles (93%) [6]. The extended-depth-of-focus-IOL (EDOF-IOL) provides an extended depth-of-field range. Actually, it is a bifocal IOL with a low near-addition of 1.5–2.0 D. This results in excellent distance and intermediate correction without the disadvantages of a multifocal IOL. However, this procedure does not completely eliminate the need for spectacles [5,7]. The third possibility of correction is the implantation of a trifocal IOL. One is located in the distance and receives 50% of the incident light, the second is at 30–40 cm and receives 30% of the light, and the third is at about 60 cm and receives 20% of the light. The strength of trifocal IOLs is that they are glasses-free. The weakness is that fixed distances are given, which requires some getting used to, and that the intermediate part is the weakest. As only a part of the light used for image generation is sharp, there is a loss of contrast, which is noticeable. In addition, undesirable optical side effects in the form of so-called halos or glare occur, which impair vision, especially at twilight and at night [5,8]. There is no ideal correction for all distances, but through objective measurement using the Visual Behavior Monitor (VBM) regarding distance measurement, ambient lighting and head posture will potentially determine a much more accurate setting and support the choice of IOL. Whatever additional technique is used, an important measure of success is the difference between the intended and the actual postoperative refractive outcome. It is necessary to have a method that is sensitive to discrepancies between the intended and the actual outcome.

Hence, to minimize potential side effects and improve patient outcomes, IOLs must be chosen with the visual needs of a patient in mind. One such method to assess a patient's individual visual needs is to determine patient expectations and visual behavior pre-operatively, to improve post-operative outcomes.

Visual behavior monitors are developed through various technological approaches, which should be as precise as possible but also cost-effective. Various studies use active lighting based on infrared LEDs. For example, one study [9] proposes a system that uses 3D vision techniques to estimate and track the 3D line of sight of a person with multiple cameras. The method is based on a simplified eye model and uses the Purkinje images of an infrared light source to determine the eye position. This information is used to estimate the line of sight. In a technological development, a system with active infrared LED illumination and one camera is implemented. Due to the LED illumination, the method can easily find the eyes and the system uses this information to locate the remaining facial features [10]. The authors suggest estimating the local direction of gaze analytically based on the pupil position. Almost all active systems described in the literature have been tested in simulated environments, but not in real ones. A moving person presents new challenges such as variable lighting, changing background, and vibrations that have to be considered in real systems. In a further work [11], an industrial prototype called Copilot is presented whose application tends to take into account the real situation. This system uses an infrared LED illumination to find the eyes. It uses a simple subtraction procedure to find the eyes and calculates only one validated parameter, the percentage eye closure (PERCLOS). This system currently works in low-light conditions [12]. All of these technologies, including ours, tend to collect even more precise data and reflect reality. Further improvements can be expected in the future by combining different technologies. Many technologies are developed in such a way that they are easy to handle and can be spread further. For example, the investigation tools are built into smartphones [13]. This is the way we initially went. However, a further simplification of the operation has become necessary after initial experiences. This led to today's VBM, which does not

require any technical knowledge of the patient. Due to the simplicity, the measurement results are much more precise and reflect reality much better.

Results from a literature review suggest that there is a greater utility of using measurable patient reported outcomes such as functioning, satisfaction, and quality of life as opposed to clinical outcomes when considering the benefits of cataract surgery [1,14].

Psychometric tests in the form of cataract surgery outcome questionnaires have been developed, identifying the limitations of performing daily tasks due to an individual's vision, a trait known as visual functioning [15,16]. Analyzing the visual behavior of cataract patients is extremely important for defining the visual strategy adopted in the patient's care.

McAlinden et al. [16] studied the responsiveness of several questionnaires in measuring patient outcomes. The investigation revealed that "Catquest-9SF," a Rasch modified 9 item questionnaire originally developed in Sweden [17], was the most responsive as measured by an effect size (ES) statistic of 1.45 (95% CI, 1.22–1.67). Although the "Catquest-9SF" includes a measure of visual functioning, a questionnaire known as the "Visual Disability Assessment (VDA)" (ES, 1.09; 95% CI, 0.86–1.32) also includes a measure of mobility of the patient, which may be attractive for patients and clinicians. Finally, the VBM measurement supports the solution finding for a planned surgery.

Morlock et al. [18] developed a Patient-Reported Spectacle Independence Questionnaire (PRSIQ) to assess spectacle independence following cataract surgery. The questionnaire was developed using conventional qualitative and modern measurement theory methods, and evidence was garnered for the use of PRSIQ as a measure of spectacle independence.

While the use of questionnaires and visual assessments in measuring patient outcomes are effective, there are limitations regarding the subjectivity of the results. Current psychometric measures may be subject to systematic error such as recall bias and could cause judgment errors resulting in patient discontentment with surgical outcomes [19,20]. Visual information based on patient behavior can be harnessed using technological tools that perform eye tracking, analyzing visual behavior under various circumstances. Eye tracking technologies can be used to attain and process objective information of a patient's everyday activities, crucial for selecting an appropriate IOL.

Integrating technology that collects objective visual behavior data into pre-operative patient assessments has shown to improve the current understanding of patient behavioral gaze. The data generated could aid clinicians to provide a highly individual evaluation of a patient's visual needs, including distance estimates of the patient's sight, and determine the right course of action to improve post-operative patient outcomes and patient satisfaction.

To date, most approaches have focused on diagnostic testing in a clinical setting. In-clinic measurements cannot fully replicate a patient's daily visual needs. The Visual Behavior Monitor (VBM, Vivior AG, Zurich, Switzerland) is a wearable device to be attached on the spectacles of a patient and perform continuous measurements as he or she goes about their daily activities. The challenge of unassisted monitoring in outpatient settings is ensuring correct use of the monitoring system and reliability of the data. This is especially the case for the senior cataract patients with compromised visual function due to age-related cataract and/or presbyopia. While aiming to develop a compact and non-obtrusive system, it is essential to ensure that the reduced form-factor of the device does not impede the handling of the device by the elderly population. Another challenge is to design the user interfaces of the system in a way that can be correctly interpreted by patients with different levels of exposure to computer and information technologies. Most of the recent progress in the information technologies has been pushed forward generally by the young generations of earlier adopters and, thus, interfaces are shaped for and by this population group, who have enough time and interest to work with new technologies. It is critically important to redesign the user experience from scratch with wearable systems for all population groups.

Inadequately designed systems ignoring the specifics of the population groups might not only be inefficient in delivering required objective information, but also lead to reduced patient satisfaction due to frustration from handling the device. This is aggravated by the patient awareness that the outcome of the surgery might depend on his or her ability to handle the system. With all the above in mind, the system should be comfortable for use by the patients of various age groups, educational backgrounds with different levels of exposure to information technology, and possibly with compromised vision.

The present study was initiated to understand the feasibility of the non-intrusive monitoring of visual behavior and VBM handling by, particularly, a cataract-age population without supervision of a healthcare professional (HCP). The study is designed to evaluate the ease of use, as well as patient perception of wearing the device. This can potentially answer the question of whether patients are willing to use and able to handle the wearable device mounted on the spectacles and are comfortable to use it in order to perform the required measurements.

#### **2. Materials and Methods**

Patients presenting for cataract surgery or refractive lens exchange were randomly selected to take part in this feasibility study to evaluate the VBM. Patients who met the inclusion criteria were recruited consecutively for the study. Upon enrolment, patients completed a modified Cat-Quest 9 questionnaire. For German-speaking sites, the questionnaire was modified to improve readability in German. Patients were trained and fitted with the VBM and instructed to wear the device for at least 36 h. Patients then returned a week later to complete a questionnaire on the user friendliness of the device and schedule surgery. In our study, we used a modified CATQ-9 questionnaire, which refers to the patient's ability to perform certain tasks. The questionnaire consists of 9 questions, with each question being rated on a scale of 1–4—from 1 for "no problems" to 4 for "very big problems." In general, it is asked whether the patient has problems with vision in everyday life and how satisfied the patient is with his or her vision. In particular, the questions ask how well the patient sees on the computer, how well he/she recognizes faces, how well he/she recognizes prices when shopping, how well he/she walks on uneven ground, how he/she watches television with subtitles, and what is his/her favorite pastime and hobbies. This covers a large part of daily life.

The following inclusion criteria are applied: Age over 18 years old, participants considering IOL procedures, participants considering refractive laser procedures, and participants willing and able to participate. As exclusion criteria, we have defined the following points for the study: Unwillingness to participate, manifestation of macular disease based on OCT(Optical Coherence Tomography) /fundus photography, corneal scar or sign of corneal decompensation, standard surgery exclusion criteria, inability to make all postoperative visits, history or manifestation of other pre-existing ocular conditions and comorbidity (amblyopia, monophthalmia), other predisposing sight-threatening ocular conditions (uveitis, diabetic retinopathy, ARMD, macular dystrophy, retinal detachment, neuro-ophthalmic disease, abnormal pupils [deformation] and iris), history of corneal or retinal surgery, evidence of iris or chorioretinal neovascularization, history of anterior segment pathology (aniridia), any scar from corneal endothelium damage (chemical burn, herpetic keratitis, cornea guttata), any medical history contraindicative of standard cataract surgery including those taking medications known to potentially complicate cataract surgery (e.g., α1a-selective alpha blocker), and participating in another clinical trial.

The primary endpoint of the study is testing the design and usability of the wearable and the station.

The Visual Behavior Monitor is a wearable device to continuously monitor visual behavior of the user. During the study, two functionally equivalent versions of the system were used, called Mark II and Mark III. The description below focuses on the latest system, Mark III.

The VBM is fixed on the spectacles frame with help of an adapter. Multiple adapters are provided to patients to allow attachment on all spectacles used by the patient. Patients can easily move the device from the adapter on one spectacle to the adapter on another when changing spectacles. The VBM is configured to start and stop measurements automatically when attached and detached from the adapter, as well as when spectacles are taken off.

After the VBM measurement was carried out, the respondent was asked the following questions: 1. Is it easy to attach the magnetic clip to the spectacles frames? 2. Is it easy is to attach the wearable

to the magnetic clip? 3. Do you always feel comfortable using the wearable? 4. Is it easy to use the wearable?

Regarding the station, the following questions were asked. 5. Is the connection between the station and the wearable always working? 6. Is the information provided on screen always clear to you?

VBM (Figure 1) incorporates the following sensors: Two optical time-of-flight distance sensors directed forward and 30◦ downward, ambient light sensor with separate red, green, and blue channels directed forward, combined ultraviolet and ambient light sensor directed upward, as well as motion and orientation sensors: Accelerometer, gyroscope, and magnetometer. The VBM does not include imaging sensors, such as cameras and location sensors, and does not infringe on the privacy of the wearer and other people. The VBM features an energy-independent real-time clock for reliable timestamping of the measurements.

**Figure 1.** Visual Behavior Monitor (VBM) fixed on the spectacles.

Measurements data are stored in the encrypted format in the device and the progress of data collection is shown on the embedded screen or on the screen of the station in hours and minutes. The onboard memory capacity is more than 1000 h. The VBM is charged overnight and the battery capacity is sufficient for a complete day of continuous measurements. When the device is returned to the clinic, data are transferred from the device to an electronic tablet and then uploaded to the internet server (cloud), decrypted, processed, and visualized via web-interface.

Cloud processing is performed by advanced machine learning algorithms, which are trained to recognize various types of patient's visual activities, such as looking on a desktop computer screen, reading handheld materials, driving a car, and viewing objects at a distance. The environment is further analyzed based on the recognized activities to derive relevant metrics such as viewing distances to the objects of visual activities, such as distance to a computer screen, to a desktop, or to handheld material. The data are further aggregated to create relevant statistical representations.

The VBM is intended to provide behavioral information to the attending HCP and to support ophthalmic surgeons or other healthcare professionals in the planning of cataract and refractive surgery. This is expected to improve participant satisfaction through the objective evaluation of individual vision needs and allowing practitioners to more adequately address those needs. Based on the individual

profile of viewing distances, it is possible to select the optimal IOL solution by overlaying the defocus curve of the IOL, representing the expected visual acuity as a function of defocus (visual distance). Additional measured parameters, such as typical light levels and their spectral content, can be also considered. For example, monofocal IOLs deliver good far vision but cannot provide clear vision in near and intermediate distance zones, while lenses with advanced optics, such as bifocal and trifocal (multifocal) lenses, as well as lenses with an extended depth of focus, are capable of providing sufficient visual acuity and spectacle independence both in the near and intermediate zones, depending on the lens. At the same time, advanced optics can cause halos and glare in low-light conditions. Thus, for a patient with far distance vision needs with a large share of low-light activities, such as night driving, it would be reasonable to recommend a monofocal lens or blended vision solution, while an advanced optics IOL would be preferable for the patient with a range of activities with the viewing distance distributed between all vision zones and requiring often switching between zones.

The information provided by the VBM allows healthcare professionals to objectively assess the needs of the patient and explain the benefits and limitations of the available solutions to the patient. However, the decision should include other factors, like personal preferences of the patient and tolerance to the multifocality of IOLs with advanced optics, as well as other conditions and, thus, should not exclude discussion with the patient. The measuring elements used in the VBM are standard components. Furthermore, the VBM is not intended to replace state-of-the art diagnostic methods in refractive treatment. Nevertheless, it is very important to have hard facts to support the choice of IOL suitability or surgical method.

Patients use the VBM in their normal daily routine, with the only excluded activities being those that involve water exposure of the wearable. Some other activities may be exempt based on a benefit/risk assessment. Those activities are meant to be communicated verbally to the attending HCP during pre-operative interview.

A dedicated tablet is provided to the clinics for measurement data upload on the cloud and data review. The embedded cellular data module of the tablet is used for data upload and report review on the cloud, independent of the clinic's infrastructure.

The clinic has access to a detailed report of all data the VBM collects. A simplified patient report including recommendations is available for the patient.

The healthcare professional has an option to review aggregated metrics of visual behavior of the patient. For example, the viewing distance from all measurement days, collected by the patient, can be displayed as a histogram of the additional optical power (refraction) of the eye required to accommodate the objects of visual activity (Figure 2).

Further, it is possible to relate the viewing distances to the particular head inclinations in order to visualize viewing distances in a two-dimensional plane (Figure 3). Such visualization is helpful to educate patients on the relation of activities and visual distances as visual zones can be linked to the actual physical space.

Additionally, the system allows us to review activities chronologically, as performed by the patient during the selected measurement day. The distribution of viewing distances can be overlaid with other relevant parameters, such as ambient light exposure (Figure 4).

**Figure 2.** Distribution of viewing distance and required additional optical power of the eye (refraction) to observe objects of visual activity. The plot demonstrates domination of the intermediate viewing distance related to the computer work at around −1.75 D (0.57 m) for the particular patient (Patient A). Color coding is used to display the viewing zone, with green representing near viewing distances, blue representing intermediate, and yellow representing far.

This kind of plot allows the demonstration of viewing conditions during various activities to the patient. The time of the day information helps the patient to relate his or her daily routine to the measurements and communicate to the HCP the importance of certain activities. The HCP is also able to relate viewing distance to other parameters. Figure 5 demonstrates the history of ultraviolet light exposure during the study day together with ambient light, which provides additional information about the patient's lifestyle.

The study has been conducted by five sites in three European countries: Germany, Switzerland, and Ireland (Table 1). A total of 129 patients were enrolled between August 2018 and May 2019. Patient age ranged from 49 to 82 (mean: 65). In all cases, the VBM measurement was performed before cataract surgery.

**Figure 4.** Example of the viewing distances and ambient light measurements with VBM during a day (Patient B). The example demonstrates the clear domination of far distances (upper plot), which is a marked difference in visual behavior from Patient A.

**Figure 5.** Example of the ultraviolet light exposure (quantified with UV index) as solid colored bars compared to ambient light level shown as empty bars during a measurement day shown in Figure 4 (Patient B). The level of ultraviolet light indicates that the example day was mostly spent outdoors.

**Table 1.** Vivior feasibility study sites, ethics committees, and enrollment.


All subjects gave their informed consent for inclusion before they participated in the study. The patients were informed in great detail about the content of the work, its usefulness, and its feasibility before the start of the study. The study was conducted in accordance with the Declaration of Helsinki, and the protocol was approved by the appropriate ethics committees for each site (Table 1).

#### **3. Results**

Of the total 178 eyes implanted, 119 eyes were implanted with monofocal intraocular lenses (67%), 39 eyes were implanted with trifocals or multifocals IOLs (22%), 15 eyes were implanted with toric lenses (8%), and the remaining five eyes were implanted with lenses that were not identified by the surgeon (3%). From these 129 patients included in the study, 178 cataract surgeries ultimately resulted, i.e., both eyes were not operated on in all cases, but only those that had an indication for surgery.

The pre-operative cataract scores distribution is shown in Figure 6 and the Cat-Quest 9 Task results are shown in Figure 7. A grade 3 or greater cataract manifested itself in 34% of the eyes, while no cataract and no reported cataract were accounted for in 29% of the eyes. The cataract and presbyopia limitations were also reflected in the pre-operative vision questionnaire with the vast majority of patients reporting some to very great difficulties in mastering everyday life and the majority of the patients being dissatisfied with the current eyesight.

**Figure 6.** Lens opacities classification of cataract study patients.

**Figure 7.** Pre-operative Cat-Quest 9 task results.

We have also investigated whether the difficulties in using the VMB system are related to the patients' visual problems. Based on the patients' answers to the usability questionnaire, the data were statistically evaluated using the parameters of corrected distant visual acuity and the Lens Opacities Classification System (LOCS) score. These results are shown in Table 2.

**Table 2.** Distribution of the preoperative corrected distance visual acuity and Lens Opacities Classification System (LOCS) III score of the patients separated by the answers in the usability questionnaire.


Surprisingly, the groups that had difficulties with the VBM presented the higher visual acuity (lower logMAR values), as well as lower cataract scores (better vision), for all the questions, except for one concerning "information provided on the screen." The study included 15 patients with both eyes cataract-free and a relatively high distance-corrected visual acuity preparing for the clear lens exchange. This group was typically more demanding toward the usability of the system and often appeared in the group answering "no" in the questionnaire in Figure 8. This resulted in the proportion of patients without cataract in the "no" group being between 20% and 40%, while the mean was around 12%. The exception was question 6 where there was only a slightly higher proportion of clear lens patients in the "no" group (16% vs. 11%). Based on this data, we were not able to statistically attribute the difficulties reported in the usability questionnaire to the vision deficiencies of the patients for all questions except question 6, which also reported the highest negative feedback (34% not agreeing that the information provided on the screen is always clear).

**Figure 8.** Results summary of usability questionnaire.

We have also studied the patient satisfaction with the vision results (Figure 9).

**Figure 9.** Post-operative Cat-Quest 9 results showing that surgical results met or exceeded patients' expectations.

A high level of patient satisfaction was reported, with 84% of patients meeting or even exceeding expectations (Figure 9).

#### **4. Discussion**

In today's world, visual acuity is not only becoming more and more important, but the demands on the visual system have increased significantly, whether at work or in leisure time. Various studies have shown that 61% of the patients prefer an explicit decision by the physician regarding the cataract surgery [21,22]. This strong dependency and importance of the physician's decision-making must be compared to the importance of patient involvement in order to maintain high patient satisfaction [23,24]. Previous studies show that the improving patients' perceived level of understanding is an important factor for the high level of patient satisfaction with the results of a cataract surgery [1,25–27]. There is an increasing demand to understand patient needs and expectations when planning cataract and refractive surgeries. It is generally accepted that understanding patient behavior and expectations and including this information in the planning process results in better surgical outcomes. Psychometric tests in the form of cataract surgery outcome questionnaires have been developed, identifying the limitations of performing daily tasks due to an individual's vision, a trait known as visual functioning [15,16]. The use of questionnaires and visual assessments in measuring patient outcomes is effective; however, the subjectivity of these methods is questionable. Nevertheless, questionnaires are very important in addition to objective measurement, because they reflect the important subjective part of the patient's feelings. Its importance and introduction have been used in the literature for a long time where specific questionnaires have been developed for specific questions [23], as we have used it in our study. Current psychometric measures may be subject to systematic error such as recall bias and could cause judgment errors resulting in patient discontent with surgical outcomes [19,20]. Integrating technology, such as the wearable VBM, which collects objective visual behavior data for pre-operative patient assessments, has shown improvement in the current understanding of patient visual behavior. It was very gratifying to note that the application of the VBM was, in the vast majority of cases, problem-free. Thus, the VBM measuring device could be installed with reasonable effort and be reliably used by the patient. In this sense, the goal was achieved to design a VBM in a way that not only works reliably but is also easy to use. Most of the study patients stated that they had severe to very severe visual problems before the surgery. None of the patients who reported very great difficulties in their daily tasks reported troubles in attaching the clip to the spectacles frame or VBM to the clip. The results of the analysis of the patients' visual performance based on their responses to the questionnaire showed that the group that was not satisfied with the user-friendliness had, on average, a higher visual acuity and a lower LOCS score, with the exception of one question related to the presentation of the information on the screen. This effect of the higher reported usability problems of patients with better vision can be attributed to the higher demands of patients with clear natural lenses. The feedback from patients on the system handling was mostly related to the adaptation to a specific spectacle frame, the correct

alignment of the VBM and clip, and the strength of the magnetic attachment. This feedback has been considered in the development of the next generation of the system. The size of the clip has been adapted to make it easier to handle, and the shape of the clip and wearable has been redesigned to visually match each other to make the attachment more intuitive. The clip has a defined arrow shape that points forward to indicate the direction of attachment. The magnetic attachment has been replaced by a purely mechanical attachment to prevent the design from falling off during heavy movement, while at the same time being easy to attach and detach. The reported difficulties in understanding the information on the station's screen and connection problems between the wearable and the station have led to the decision to make the next unit self-sufficient and to omit an additional station for data acquisition, thus eliminating both station-related difficulties.

The objective data generated could aid clinicians in providing a highly individualized evaluation of a patient's visual needs, including estimated viewing distances, and determine the right course of action to improve post-operative patient outcomes and patient satisfaction. The subjective and objective data of patients often coincide. However, the patient's subjective feelings are not always correctly interpreted by them. The objective data collected can close this circle to a higher percentage and thus have potential to increase patient satisfaction. Including objective data of the visual behavior of each individual patient is another important step for a successful customized treatment.

Prior to cataract- and refractive surgery in addition to a comprehensive eye exam in terms of the overall health of the eyes, refraction measurements (amount of nearsightedness, farsightedness, astigmatism), determining the corneal curvature, as well as the length of the eye, it is particularly important to analyze the environmental conditions in the eye and ophthalmological concomitant diseases. The objectively measured ophthalmological findings, as well as the VBM data collected, have the potential to further improve treatment outcomes.

There is no ideal correction for all distances currently available, but the objective measurements of viewing distances, ambient light conditions, and head inclinations with devices such as the VBM will potentially allow better characterization of a patient's lifestyle and support the choice of IOL solutions better matching a specific patient's lifestyle and visual needs. Further studies with extended questions are ongoing.

#### **5. Conclusions**

To test the VBM device usability among patients, a multi-center feasibility study was conducted. Patients found the device easy to use, with most reporting that the device was not intrusive. The patients, even those that were elderly with significant cataracts, were able to comfortably use the VBM. Hence, the results of this feasibility study demonstrated that the concept of wearable monitoring of visual behavior is, in general, accepted by the patients. VBM may be used correctly and reliably by the patient. The difficulties of handling the system by the patients have been considered in the design of the next generation of the VBM.

The objective data collected with the VBM may have the potential to improve treatment outcomes, and the next research phase will seek to objectively assess the impact of the VBM on surgeons' treatment decision making.

**Author Contributions:** B.P. provided the treatment indication, developed the study design, acquired clinical data, and contributed to writing the paper. Z.C. contributed substantially to the methodology and substantially contributed to writing the paper. B.P.-E. performed substantial data curation. P.Z. performed data analysis and contributed to the study design. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Conflicts of Interest:** Bojan Pajic, Brigitte Pajic-Eggspuehler, and Zeljka Cvejic declare no financial interest. It is to declare that Pavel Zakharov is an employee of Vivior AG.

#### **References**


© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## *Article* **Pressure Observer Based Adaptive Dynamic Surface Control of Pneumatic Actuator with Long Transmission Lines**

#### **Deyuan Meng 1,\*, Bo Lu 2, Aimin Li 1, Jiang Yin <sup>1</sup> and Qingyang Li <sup>1</sup>**


**\*** Correspondence: tinydreams@126.com; Tel.: +86-1516-2126-417

Received: 22 July 2019; Accepted: 26 August 2019; Published: 3 September 2019

#### **Featured Application: This research dedicates to develop a pressure observer based nonlinear controller for precise position control of the MRI-compatible pneumatic servo system with long transmission lines.**

**Abstract:** In this paper, the needle insertion motion control of a magnetic resonance imaging (MRI) compatible robot, which is actuated by a pneumatic cylinder with long transmission lines, is considered and a pressure observer based adaptive dynamic surface controller is proposed. The long transmission line is assumed to be an intermediate chamber connected between the control valve and the actuator in series, and a nonlinear first order system model is constructed to characterize the pressure losses and time delay brought by it. Due to the fact that MRI-compatible pressure sensors are not commercially available, a globally stable pressure observer is employed to estimate the chamber pressure. Based on the model of the long transmission line and the pressure observer, an adaptive dynamic surface controller is further designed by using the dynamic surface control technique. Compared to the traditional backstepping design method, the proposed controller can avoid the problem of "explosion of complexity" since the repeated differentiation of virtual controls is no longer required. The stability of the closed-loop system is analytically proven by employing the Lyapunov theory. Extensive experimental results are presented to demonstrate the effectiveness and the performance robustness of the proposed controller.

**Keywords:** pneumatic servo system; long transmission line; pressure observer; dynamic surface control; position tracking

#### **1. Introduction**

The magnetic resonance imaging (MRI) technique is widely used in clinical diagnosis due to its ability to image without the use of ionizing x-rays and superior soft tissue contrast as compared to computed tomography (CT) scanning. Recently, pneumatically actuated MRI-compatible robots, which enable real-time magnetic resonance (MR) image-guided needle placement, are designed for brachytherapy and biopsy by several researchers [1–15]. To fulfill the requirements for MR compatibility of the robotic systems, pneumatic valves are commonly placed outside the scanner room in the aforementioned works. Therefore, long transmission lines between the actuators and valves are used. Since long transmission lines have a significant influence on the pressure dynamics of the pneumatic system and MRI-compatible pressure sensors are not commercially available, precise position control is one of the main technical challenges in robot development.

Indeed, considerable research effort has been devoted to addressing the issue related to long transmission lines. Richer et al. [16] suggested a formula for the time delay and amplitude attenuation between the mass flows at the outlet and inlet of the line, as well as a formula for the pressure drop along the tube. Although this approach was later employed in many other studies (Jiang et al. [14]; Richer et al. [17]), it has been proven improper for long transmission lines (>2 m) [18]. Yang et al. [8] adopted a first order transfer function with time delay to describe the 9-m transmission line dynamics. However, for some unknown reason, the authors omitted the transmission line dynamics in the subsequent controller design. Based on the discretization of a general set of equations from fluid mechanics, two similar distributed models of long pneumatic transmission lines were derived by Li et al. [19] and Krichel et al. [20]. However, this type of model is not suitable for controller design for their high order and low computational speed. Turkseven et al. [21] developed a simplified distributed model to characterize the additional dynamics brought by the long transmission lines.

As previously mentioned, direct pressure measurement is unavailable for the lack of MRI-compatible pressure sensors. However, pressure states are commonly used in controller design for precise position control of a pneumatic actuator. To solve this problem, pressure observers were used in place of a pressure sensor in many studies. Pandian et al. [22] investigated two design methods of an observer to estimate the chamber pressures in a pneumatic cylinder. While a continuous gain observer was used to estimate one of the chamber pressures with the assumption that the other one is measured in the first method, a sliding mode observer was utilized to estimate both chamber pressures in the second method. In this study, the mass flow rate was assumed to be exactly known, which is apparently too restrictive in practice. Wu et al. [23] conducted an analysis on observability and concluded that the pressure states in the pneumatic servo system are not locally observable from the measurement of the output motion within several regions in the state space. Based on the actuator pressure dynamics, various observers are developed, such as the sliding-mode observer (Bigras [24]), energy-based Lyapunov observer (Gulati et al. [25]), and adaptive nonlinear observer (Langjord et al. [26]). It should be noted that the performance of these observers relies on the accuracy of the valve model. Driver et al. [27] developed a pressure estimation algorithm by utilizing the measured actuation force and the hypothesized average air pressure in the actuator. Turkseven et al. [18] presented a pressure observation method by using the measured force and piston displacement. However, the requirement for force sensing in the MRI environment is particularly burdensome because the MRI-compatible force sensor is not available commercially and its development cost is high.

Recently, the backstepping design method has been proven to be a very effective way to develop nonlinear robust controllers for pneumatic servo systems. However, this method has the problem of "explosion of complexity" since the requirement of repeated differentiation of virtual controls. Thus, a practical implementation is difficult. To solve this problem, Swaroop et al. [28] proposed the dynamic surface control (DSC) method, in which the calculation of the virtual control variable's derivative was prevented by introducing a filter at each design step. Since then, the DSC method has been the topic of significant research efforts and a number of excellent theoretical contributions have been made. The applications of DSC can be found in many engineering fields, for example, hydraulic servo systems [29], underwater/autonomous surface vehicles [30], mobile wheeled inverted pendulum [31], pneumatic artificial muscle [32], and servo motor [33].

In this study, the needle insertion motion control of the MRI compatible robot developed in our lab is considered. The robot is actuated by a pneumatic cylinder with long transmission lines. The focus of this paper is dealing with the issue of long transmission lines and realizing a high accuracy control of a pneumatic actuator with a pressure observer. Therefore, the long transmission line is assumed to be an intermediate chamber connected between the control valve and the actuator in series, and a pressure observer based adaptive dynamic surface control is proposed. The main contributions are: (1) The long transmission line dynamic is approximated as a nonlinear first order system, which can estimate the pressure losses and time delay through the long transmission line precisely in real time. (2) A pressure observer, which is proven to be globally stable, is developed to estimate the chamber pressure in a pneumatic actuator with a long transmission line. (3) In contrast to most of the existing nonlinear robust controllers synthesized by the backstepping method, by using the dynamic surface control (DSC) technique, the proposed controller can cope with the problem of "explosion of complexity", since the repeated differentiation of virtual controls is no longer required. The rest of this paper is organized as follows: Section 2 presents the dynamic models and problem statement; Section 3 gives the design and stability proof of the pressure observer based adaptive dynamic surface controller; Section 4 presents the experimental results to demonstrate the performance of the proposed controller; and Section 5 draws the conclusions.

#### **2. Dynamic Models and Problem Formulation**

As shown in Figure 1, a 5-DOF (degree of freedom) MRI compatible robot for abdominal and thoracic punctures was built in our laboratory. The robot's mechanism design was developed such that all motions were decoupled and actuated by pneumatic cylinders. With the help of the robot, the physician could manipulate the needle remotely without moving the patient out of the MRI scanners. A schematic of the pneumatic servo system driving the needle insertion motion is depicted in Figure 2. The pneumatic cylinder (customized product from XMC Corp., Ningbo, China), which had a 150 mm stroke and 10 mm diameter bore, was made of nonmagnetic material to ensure MRI compatibility. While Chamber A was controlled by a proportional directional control valve (FESTO MPYE-5-M5-010-B), a tank was used to maintain the pressure of Chamber B at a specified constant level for safety reasons. The control valve was connected to the cylinder via a 10 m transmission line since it had to be placed away from the MRI scanner. A pressure observer was needed to estimate the pressure of Chamber A, since direct measurement was expensive for the lack of MRI-compatible pressure sensors. The main purpose of this paper was to find a way to deal with the issue of the long transmission line and realize high accuracy control of a pneumatic cylinder with pressure observation. For the purpose of comparison, the pressures at the two ends of the long transmission line, as well as the supply pressure were measured by three pressure sensors (FESTO SDET-22T-D10-G14-I-M12). The piston position was obtained with an optical position sensor (Micronor MR328). The controller was programmed in Simulink on a dSPACE DS1103 control system with a 1 ms sampling period.

**Figure 1.** Virtual prototype of a 5-DOF (degree of freedom) magnetic resonance imaging (MRI) compatible robot for abdominal and thoracic punctures.

**Figure 2.** Schematic of the pneumatic servo system used to drive the needle insertion motion.

As shown in [34], the motion of the piston–rod–needle assembly can be expressed as:

$$
\dot{\mathbf{p}} \cdot m\ddot{\mathbf{x}} = p\_\mathbf{a} A\_\mathbf{a} - p\_\mathbf{b} A\_\mathbf{b} - p\_\mathbf{0} A\_\mathbf{r} - b\dot{\mathbf{x}} - F\_\mathbf{f} S\_\mathbf{f}(\dot{\mathbf{x}}) - d\_\mathbf{f},\tag{1}
$$

where *x*, . *<sup>x</sup>*, and .. *x* are the piston position, velocity, and acceleration, respectively, *m* is the normal mass of the piston–rod–needle assembly, *p*<sup>a</sup> and *p*<sup>b</sup> denotes the absolute pressures of actuator chambers, *p*<sup>0</sup> is the ambience pressure, *A*<sup>a</sup> and *A*<sup>b</sup> are the cross section areas of piston chambers, *A*<sup>r</sup> is the cross section area of the rod, *b* is the total load and cylinder viscous friction coefficient, *F*<sup>f</sup> is the unknown friction coefficient, *S*f( . *x*) is a continuous function, which is always chosen as *S*f( . *x*) = <sup>2</sup> <sup>π</sup> arctan(<sup>900</sup> . *x*), *F*f*S*f( . *x*) is utilized as the smooth approximation for the usual static discontinuous Coulomb friction force, and *d*<sup>f</sup> represents the unmodeled dynamics and external disturbances.

With the assumption that the discharging and charging processes are both isothermal, the pressure dynamic in the actuator chamber can be modeled as [16,35]:

$$
\dot{p}\_{\text{a}} = \frac{RT\_{\text{s}}}{V\_{\text{a}}} \dot{m}\_{\text{a}} - \frac{A\_{\text{a}} p\_{\text{a}}}{V\_{\text{a}}} \dot{\mathbf{x}}, \quad V\_{\text{a}} = V\_{\text{a}0} + A\_{\text{a}} (\frac{L\_{\text{c}}}{2} + \mathbf{x}), \tag{2}
$$

where *<sup>R</sup>* is the gas constant, *Ts* is the ambient temperature, *<sup>V</sup>*<sup>a</sup> is the volume of the cylinder chamber A, . *m*<sup>a</sup> is the mass flow entering or exiting the chamber A, *V*a0 is the dead volume of the cylinder chamber A, and *<sup>L</sup>*<sup>c</sup> is stroke of the actuator. . *m*a can be calculated by [16]:

$$
\dot{m}\_{\text{a}} = (p\_l - p\_{\text{a}}) \frac{p\_{\text{v}} A\_l D\_l}{32 R T\_{\text{s}} L\_l \mu} \,\prime \tag{3}
$$

where *pl* is the average air pressure in the transmission line, *p*<sup>v</sup> is the measured air pressure at the work port of the control valve, *Ll*, *Al*, and *Dl* denote the length, cross section area, and inner diameter of the transmission line, respectively, and μ is the dynamic viscosity of air.

As shown in Figure 2, the long transmission line is assumed to be an intermediate chamber connected between the control valve and the actuator in series, thus, the following equation is formulated to represent the pressure dynamic of the line [16,35]:

$$
\dot{p}\_l = \frac{RT\_\mathbf{s}}{V\_I} (\dot{m}\_\mathbf{v} - \dot{m}\_\mathbf{a}),
\tag{4}
$$

where *Vl* is the volume of the transmission line, and . *m*v is the mass flow rate through the valve, which can be calculated by [35]:

$$\dot{m}\_{\rm V} = A(u) \mathbb{K}\_{\rm q}(p\_{\rm u}, p\_{\rm d}, T\_{\rm u}) = \begin{cases} A(u) \mathbb{C}\_{\rm d} \mathbb{C}\_{1} \frac{p\_{\rm u}}{\sqrt{T\_{\rm u}}} \prime & \frac{p\_{\rm d}}{p\_{\rm u}} \le p\_{\rm r} \\\ A(u) \mathbb{C}\_{\rm d} \mathbb{C}\_{1} \frac{p\_{\rm u}}{\sqrt{T\_{\rm u}}} \sqrt{1 - \left(\frac{p\_{\rm u}}{1 - p\_{\rm r}}\right)^{2}} & \frac{p\_{\rm d}}{p\_{\rm u}} > p\_{\rm r} \end{cases} \tag{5}$$

where *A*(*u*) is the effective valve orifice area, *u* is the control valve's control input, *C*<sup>d</sup> is the discharge coefficient, *C*<sup>1</sup> is a constant with a value of 0.0404, *p*<sup>u</sup> and *p*<sup>d</sup> are the upstream pressure and the downstream pressure, respectively, *T*u is the upstream temperature of air, and *p*r is the critical pressure ratio. Therefore, . *<sup>m</sup>*<sup>v</sup> <sup>=</sup> *<sup>A</sup>*(*u*)*K*q(*p*s, *<sup>p</sup>*v, *<sup>T</sup>*s) when the line is charging, and . *m*<sup>v</sup> = *A*(*u*)*K*q(*p*v, *p*0, *T*s) when the line is discharging.

Choose the state vectors *x* = [*x*1, *x*2, *x*3, *x*4] <sup>T</sup> as *<sup>x</sup>*<sup>1</sup> <sup>=</sup> *<sup>x</sup>*, *<sup>x</sup>*<sup>2</sup> <sup>=</sup> . *x*, *x*<sup>3</sup> = *p*<sup>a</sup> and *x*<sup>4</sup> = *pl*, thus the system dynamics can be expressed in a state-space form as:

$$\begin{cases} \dot{\mathbf{x}}\_1 = \mathbf{x}\_2\\ m\dot{\mathbf{x}}\_2 = A\_\mathbf{a}\mathbf{x}\_3 - A\_\mathbf{b}p\_\mathbf{b} - p\_0 A\_\mathbf{r} - b\mathbf{x}\_2 - A\_\mathbf{f} S\_\mathbf{f}(\mathbf{x}\_2) - d\_\mathbf{f} \\\ \dot{\mathbf{x}}\_3 = \frac{RT\_\mathbf{a}}{V\_\mathbf{a}}\dot{m}\_\mathbf{a} - \frac{A\_\mathbf{a}\mathbf{x}\_3}{V\_\mathbf{a}}\mathbf{x}\_3\\\ \dot{\mathbf{x}}\_4 = \frac{RT\_\mathbf{a}}{V\_l}(\dot{m}\_\mathbf{V} - \dot{m}\_\mathbf{a}) \end{cases} . \tag{6}$$

The control objective is to synthesize a control input *u* for the system Equation (6) such that *x* tracks the desired trajectory *x*<sup>d</sup> with a guaranteed transient and final tracking accuracy. Instead of using the pressure sensor, a globally stable pressure observer is developed to estimate the pressure in chamber A. The parametric uncertainties due to unknown *b* and *A*<sup>f</sup> and unknown external disturbances will be explicitly considered in this paper.

**Assumption 1:** *The desired trajectory xd is at least second-order di*ff*erentiable, and xd,* . *<sup>x</sup>*d*, and* .. *x*<sup>d</sup> *are bounded. Then, there exists a compact set* Ω<sup>0</sup> = *x*d, . *x*d, .. *x*d <sup>T</sup> : *<sup>x</sup>*<sup>2</sup> <sup>d</sup> <sup>+</sup> . *x* 2 <sup>d</sup> <sup>+</sup> .. *x* 2 <sup>d</sup> ≤ *B*<sup>0</sup> *such that x*d, . *x*d, .. *x*d <sup>T</sup> <sup>∈</sup> <sup>Ω</sup>0*, where B*<sup>0</sup> *is a positive constant.*

**Assumption 2:** *The extent of parametric uncertainties and external disturbances can be predicted and given by b*min ≤ *b* ≤ *b*max, *F*fmin ≤ *F*<sup>f</sup> ≤ *F*fmax*, and d*min ≤ *d*<sup>f</sup> ≤ *d*max.

#### **3. Controller Design**

#### *3.1. Pressure Observer*

The following observer is proposed to estimate the pressure of the actuator chamber *A* and the average air pressure in the transmission line:

$$\begin{cases}
\dot{\rho}\_{\mathbf{a}} = \frac{RT\_{\mathbf{a}}}{V}\hat{\dot{m}}\_{\mathbf{a}} - \frac{A\_{\theta}\dot{\mathbf{x}}}{V\_{\mathbf{a}}}\hat{\rho}\_{\mathbf{a}} \\
\dot{\hat{\rho}}\_{l} = \frac{RT\_{\mathbf{a}}}{V\_{l}}(\hat{\dot{m}}\_{\mathbf{v}} - \dot{\hat{m}}\_{\mathbf{a}})
\end{cases},\tag{7}$$

where *<sup>p</sup>*ˆa and *<sup>p</sup>*ˆ*<sup>l</sup>* represent the estimates of *<sup>p</sup>*<sup>a</sup> and *pl*, <sup>ˆ</sup> . *<sup>m</sup>*<sup>a</sup> and <sup>ˆ</sup> . *m*v denote the estimated mass flow rates according to Equation (3) and Equation (5) based on *p*ˆa and *p*ˆ*l*. It should be noted that the proposed observers are closed loop because of the relationship between the estimated mass flow rates and the estimated pressures. Proof of the convergence between the actual pressures and the estimated ones will be given later.

To verify the convergence between the actual pressures and the estimated pressures, the following Lyapunov function candidate was considered:

$$V\_{\rm o} = V\_{\rm o1} + V\_{\rm o2}, V\_{\rm o1} = \frac{1}{2} (\widetilde{p}\_{\rm a} V\_{\rm a})^2,\\ V\_{\rm o2} = \frac{1}{2} (\widetilde{p}\_{\rm l} V\_{\rm l})^2,\tag{8}$$

where *<sup>p</sup>*<sup>a</sup> <sup>=</sup> *<sup>p</sup>*<sup>a</sup> <sup>−</sup> *<sup>p</sup>*ˆa, and *pl* <sup>=</sup> *pl* <sup>−</sup> *<sup>p</sup>*ˆ*l*.

The time derivative of *V*o1 is:

$$\begin{array}{lcl}\dot{V}\_{\text{o1}} &= \left(p\_{\text{a}}V\_{\text{a}} - \hat{p}\_{\text{a}}V\_{\text{a}}\right)\left(\dot{p}\_{\text{a}}V\_{\text{a}} + p\_{\text{a}}\dot{V}\_{\text{a}} - \dot{\hat{p}}\_{\text{a}}V\_{\text{a}} - \hat{p}\_{\text{a}}\dot{V}\_{\text{a}}\right) \\ &= RT\_{\text{s}}V\_{\text{a}}\left(p\_{\text{a}} - \hat{p}\_{\text{a}}\right)\left(\dot{m}\_{\text{a}} - \hat{\dot{m}}\_{\text{a}}\right) \end{array} \tag{9}$$

According the Equation (3), one can obtain:


Therefore, . *V*o1 is negative semi-definite. The time derivative of *V*o2 is:

$$\begin{array}{rcl} \dot{V}\_{\text{o2}} &= (p\_l V\_l - \dot{p}\_l V\_l)(\dot{p}\_l V\_l - \dot{\hat{p}}\_l V\_l) \\ &= R T\_s V\_l (p\_l - \hat{p}\_l)(\dot{m}\_\text{V} - \dot{\hat{m}}\_\text{V}) - R T\_s V\_l (p\_l - \hat{p}\_l)(\dot{m}\_\text{a} - \hat{\hat{m}}\_\text{a}) \end{array} \tag{10}$$

According the Equation (3), one can obtain:


For the mass flow rate relationship given by Equation (5), it has proven in [25] that the term (*pl* <sup>−</sup> *<sup>p</sup>*ˆ*l*)( . *<sup>m</sup>*<sup>v</sup> <sup>−</sup> <sup>ˆ</sup> . *<sup>m</sup>*v) is always non-positive. Thus, . *V*o2 and *V*<sup>o</sup> are negative semi-definite, one can conclude that the proposed pressure observers are globally Lyapunov stable with regard to the errors in the estimated pressures.

#### *3.2. Adaptive Dynamic Surface Controller*

*Step* 1: Differentiate the trajectory tracking error *e*<sup>1</sup> = *x*<sup>1</sup> − *x*<sup>d</sup> with respect to time leads to:

$$
\dot{\mathbf{e}}\_1 = \mathbf{x}\_2 - \dot{\mathbf{x}}\_d.\tag{11}
$$

Consider *x*<sup>2</sup> as the first virtual control input, the virtual control law is designed as follows:

$$
\overline{\mathbf{x}}\_{\text{2d}} = \dot{\mathbf{x}}\_{\text{d}} - k\_1 \mathbf{e}\_{1\prime} \qquad k\_1 > 0,\tag{12}
$$

where *k*<sup>1</sup> is a positive control gain. Following the dynamic surface control theory [28], the signal *x*2d is fed to a low pass filter to obtain a new one *x*2d for the next design step:

$$
\pi\_1 \dot{\mathbf{x}}\_{\text{2d}} + \mathbf{x}\_{\text{2d}} = \overline{\mathbf{x}}\_{\text{2d}}, \qquad \mathbf{x}\_{\text{2d}}(0) = \overline{\mathbf{x}}\_{\text{2d}}(0), \qquad \tau\_1 > 0,\tag{13}
$$

where τ<sup>1</sup> is the filter parameter.

*Step* 2: Define the first surface error as *s*<sup>1</sup> = *x*<sup>2</sup> − *x*2d, its time derivative can be derived as:

$$\dot{s}\_1 = \dot{\mathbf{x}}\_2 - \dot{\mathbf{x}}\_{2d} = \overline{A}\_\mathbf{a} \mathbf{x}\_3 - \frac{1}{m} \Big[ p\_\mathbf{b} A\_\mathbf{b} + p\_0 A\_\mathbf{r} + b \mathbf{x}\_2 + F\_\mathbf{f} \mathbf{S}\_\mathbf{f}(\dot{\mathbf{x}}) + d\_\mathbf{f} \Big] - \frac{1}{\tau\_1} (\overline{\mathbf{x}}\_{2d} - \mathbf{x}\_{2d}),\tag{14}$$

where *A*<sup>a</sup> = *A*a/*m*.

Define *<sup>x</sup>*<sup>3</sup> <sup>=</sup> *<sup>x</sup>*<sup>3</sup> <sup>−</sup> *<sup>x</sup>*ˆ3, and let <sup>ς</sup> <sup>=</sup> *<sup>x</sup>*ˆ3 be the second virtual control input, the virtual control law is designed as follows:

$$\overline{\boldsymbol{\varepsilon}}\_{\rm d} = \frac{1}{\overline{A}\_{\rm a}} \left\{ \frac{1}{m} [p\_{\rm b} \boldsymbol{A}\_{\rm b} + p\_{0} \boldsymbol{A}\_{\rm r} + \hat{\boldsymbol{b}} \mathbf{x}\_{2} + \hat{\boldsymbol{F}}\_{\rm f} \mathbf{S}\_{\rm f}(\dot{\boldsymbol{x}}) + \hat{\boldsymbol{d}}\_{\rm f}] + \frac{1}{\tau\_{1}} (\overline{\mathbf{x}}\_{\rm 2d} - \mathbf{x}\_{2d}) - k\_{2} \mathbf{s}\_{1} - \frac{h\_{1}^{2}}{4 \eta\_{1}} \mathbf{s}\_{1} \right\},\tag{15}$$

where *k*<sup>2</sup> > 0 is a positive control gain, η<sup>1</sup> > 0 is a controller parameter, and *h*<sup>1</sup> is a known function that will be determined later, ˆ *b*, *F*ˆ f,and ˆ *d*<sup>f</sup> are the estimates of *b*, *F*f, and *d*f, respectively. ˆ *b*, *F*ˆ f, and ˆ *d*<sup>f</sup> are updated by: .

$$\dot{\hat{b}} = \text{Proj}\_{\hat{b}}(-\frac{\mathcal{V}\_1}{m}\mathbf{x}\_2\mathbf{s}\_1),\tag{16}$$

$$\dot{\mathcal{F}}\_{\mathbf{f}} = \text{Proj}\_{\mathbf{f}\_{\mathbf{f}}^{\*}}(-\frac{\mathcal{V}\_{2}}{m}S\_{\mathbf{f}}(\mathbf{x}\_{2})\mathbf{s}\_{1}),\tag{17}$$

$$\dot{d}\_{\mathbf{f}} = \text{Proj}\_{\mathcal{d}\_{\mathbf{f}}}(-\frac{\mathcal{Y}^{\mathbf{3}}}{m}\mathbf{s}\_{1}),\tag{18}$$

where γ<sup>1</sup> > 0, γ<sup>2</sup> > 0,and γ<sup>3</sup> > 0 are the observer gains, and the projection mapping is defined as:

$$\text{Proj}\_{\boldsymbol{\xi}}(\cdot) = \begin{cases} 0, & \text{if} \quad \boldsymbol{\xi} = \boldsymbol{\xi}\_{\text{max}} \quad \text{and} \quad \cdot > 0 \\\ 0, & \text{if} \quad \boldsymbol{\xi} = \boldsymbol{\xi}\_{\text{min}} \quad \text{and} \quad \cdot < 0 \\\ \cdot, & \text{otherwise} \end{cases} \tag{19}$$

where ξ is a symbol that can be replaced by *b*, *F*f, and *d*f.

Similarly, the signal ς<sup>d</sup> is fed to a low pass filter to obtain a new one ς<sup>d</sup> for the next design step:

$$
\tau\_2 \dot{\varsigma}\_\mathrm{d} + \varsigma\_\mathrm{d} = \overline{\varsigma}\_\mathrm{d}, \varsigma\_\mathrm{d}(0) = \overline{\varsigma}\_\mathrm{d}(0), \ \tau\_2 > 0,\tag{20}
$$

where τ<sup>2</sup> is the filter parameter.

*Step* 3: Define the second surface error as *s*<sup>2</sup> = ς − ςd, its time derivative can be derived as:

$$\begin{split} \dot{\varsigma}\_{2} &= \dot{\varsigma}\_{\text{r}} - \dot{\varsigma}\_{\text{d}} = \dot{\mathfrak{X}}\_{3} - \dot{\varsigma}\_{\text{d}} = \frac{RT\_{\text{s}}}{V\_{\text{a}}} \dot{\mathfrak{m}}\_{\text{a}} - \frac{A\_{\text{s}} \mathfrak{x}\_{2}}{V\_{\text{a}}} \dot{\mathfrak{x}}\_{3} - \frac{1}{\tau\_{2}} (\overline{\zeta}\_{\text{d}} - \zeta\_{\text{d}}) \\ &= -\frac{RT\_{\text{s}} \mathfrak{a}}{V\_{\text{a}}} \dot{\mathfrak{x}}\_{4} + \left(\frac{RT\_{\text{s}} \mathfrak{a}}{V\_{\text{a}}} - \frac{A\_{\text{s}} \mathfrak{x}\_{2}}{V\_{\text{a}}}\right) \dot{\mathfrak{x}}\_{3} - \frac{1}{\tau\_{2}} (\overline{\zeta}\_{\text{d}} - \zeta\_{\text{d}}) \end{split} \tag{21}$$

where *a* = *<sup>p</sup>*v*AlDl* <sup>32</sup>*RT*s*Ll*<sup>μ</sup> .

Let ζ = −*x*ˆ4 be the third virtual control input, the virtual control law is designed as follows:

$$\overline{\zeta}\_{\rm d} = \frac{V\_{\rm a}}{RT\_{\rm s}a} \left[ \left( \frac{A\_{\rm a} \chi\_2}{V\_{\rm a}} - \frac{RT\_{\rm s}a}{V\_{\rm a}} \right) \mathfrak{k}\_3 + \frac{1}{\tau\_2} (\overline{\zeta}\_{\rm d} - \zeta\_{\rm d}) - k\_3 s\_2 \right] \tag{22}$$

where *k*<sup>3</sup> > 0 is a positive control gain. Similarly, the signal ζ<sup>d</sup> is fed to a low pass filter to obtain a new one ζ<sup>d</sup> for the next design step:

$$
\tau\_3 \dot{\zeta}\_\mathbf{d} + \zeta\_\mathbf{d} = \overline{\zeta}\_\mathbf{d}, \quad \zeta\_\mathbf{d}(0) = \overline{\zeta}\_\mathbf{d}(0), \ \tau\_3 > 0,\tag{23}
$$

where τ<sup>3</sup> is the filter parameter.

*Step* 4: Define the second surface error as *s*<sup>3</sup> = ζ − ζd, its time derivative can be derived as:

$$
\dot{s}\_3 = \dot{\zeta} - \dot{\zeta}\_\mathbf{d} = -\dot{\mathbf{\dot{x}}}\_\mathbf{d} - \dot{\zeta}\_\mathbf{d} = -\frac{RT\_\mathbf{s}}{V\_I} (\hat{\dot{m}}\_\mathbf{v} - \hat{\dot{m}}\_\mathbf{a}) - \frac{1}{\tau\_3} (\overline{\zeta}\_\mathbf{d} - \zeta\_\mathbf{d}).\tag{24}
$$

Similarly, the following control law *<sup>q</sup>*<sup>d</sup> for <sup>ˆ</sup> . *m*v is proposed:

$$q\_{\rm d} = \hat{\dot{m}}\_{\rm a} - \frac{RT\_{\rm s}}{V\_{I}} \left[ \frac{1}{\tau\_{\rm 3}} (\overline{\zeta}\_{\rm d} - \zeta\_{\rm d}) - k\_{\rm 4} s\_{\rm 3} \right] \tag{25}$$

where *k*<sup>4</sup> > 0 is a positive control gain.

Once the *q*<sup>d</sup> is calculated, the desired effective valve orifice area *A*(*u*) can be calculated by:

$$A(\mu) = \begin{cases} \frac{q\_{\rm d}}{K\_{\rm q}(p\_{\rm s}p\_{\rm V}, T\_{\rm s})} & q\_{\rm d} > 0\\ \frac{q\_{\rm d}}{K\_{\rm q}(p\_{\rm V}p\_{\rm V}, T\_{\rm s})} & q\_{\rm d} \le 0 \end{cases}.\tag{26}$$

Thus, the input signal *u* for the proportional-directional control valve could be obtained according to the relation between the input signal and effective valve orifice area.

The proposed pressure observer based adaptive dynamic surface control is illustrated in Figure 3.

**Figure 3.** Block diagram of the control system.

*3.3. Proof of Stability*

**Theorem 1.** *Consider the closed loop pneumatic servo system consisting of the plant Equation (6), the nonlinear control law Equation (25), the pressure observers Equation (7), and the parameter and disturbance estimation algorithm Equations (16)–(18) under the Assumption 1–2. If there exists a set of the feedback gains and the filter constants satisfying* <sup>γ</sup> <sup>=</sup> min *<sup>k</sup>*<sup>1</sup> <sup>−</sup> 1, *<sup>k</sup>*<sup>2</sup> <sup>−</sup> <sup>1</sup> <sup>2</sup> <sup>−</sup> *<sup>A</sup>*a, *<sup>k</sup>*<sup>3</sup> <sup>−</sup> *<sup>A</sup>*<sup>a</sup> <sup>2</sup> , <sup>1</sup> <sup>τ</sup><sup>1</sup> <sup>−</sup> 1, <sup>1</sup> <sup>τ</sup><sup>2</sup> <sup>−</sup> <sup>1</sup> <sup>2</sup> <sup>−</sup> *<sup>A</sup>*<sup>a</sup> 2 > 0 *, then the closed loop system is uniformly and ultimately bounded.*

**Proof.** Consider the following Lyapunov function candidate:

$$\begin{array}{ll} V\_{\mathsf{C}} = V\_{\mathsf{C}\mathsf{S}} + V\_{\mathsf{C}\mathsf{x}} + V\_{\mathsf{O}\_{\mathsf{V}}} \ V\_{\mathsf{C}\mathsf{S}} = \frac{1}{2} \mathbf{c}\_{1}^{2} + \frac{1}{2} \mathbf{s}\_{1}^{2} + \frac{1}{2} \mathbf{s}\_{2}^{2} + \frac{1}{2} \mathbf{s}\_{3}^{2} \ \ V\_{\mathsf{C}\mathsf{x}} = \frac{1}{2} \mathbf{z}\_{1}^{2} + \frac{1}{2} \mathbf{z}\_{2}^{2} + \frac{1}{2} \mathbf{z}\_{3}^{2} \\\ V\_{\mathsf{O}} = \frac{1}{2} \mathbf{y}\_{1}^{-1} \mathbf{\widehat{b}}^{2} + \frac{1}{2} \mathbf{y}\_{2}^{-1} \mathbf{\widehat{F}}\_{\mathsf{f}}^{2} + \frac{1}{2} \mathbf{y}\_{3}^{-1} \mathbf{\widehat{d}}\_{\mathsf{f}}^{2} \end{array} \tag{27}$$

where *z*<sup>1</sup> = *x*2d − *x*2d, *z*<sup>2</sup> = ς<sup>d</sup> − ςd, *z*<sup>3</sup> = ζ<sup>d</sup> − ζd. Substituting Equation (12) into Equation (11) gives:

$$
\dot{\varepsilon}\_1 = s\_1 + z\_1 - k\_1 \varepsilon\_1. \tag{28}
$$

Substituting Equation (15) into Equation (14) and noting *x*ˆ3 = *s*<sup>2</sup> + ς<sup>d</sup> = *s*<sup>2</sup> + *z*<sup>2</sup> + ς<sup>d</sup> yields:

$$\dot{s}\_1 = \overline{A}\_\mathbf{a} (s\_2 + z\_2) - \frac{1}{m} \left[ \overline{b} \mathbf{x}\_2 + \overline{F}\_\mathbf{f} \mathbf{S}\_\mathbf{f} (\dot{\mathbf{x}}) + \overline{d}\_\mathbf{f} \right] + \overline{A}\_\mathbf{a} \overline{x}\_3 - k\_2 s\_1 - \frac{h\_1^2}{4 \eta\_1} s\_1 \tag{29}$$

where *<sup>b</sup>* <sup>=</sup> *<sup>b</sup>* <sup>−</sup> <sup>ˆ</sup> *<sup>b</sup>*, *<sup>F</sup>*<sup>f</sup> <sup>=</sup> *<sup>F</sup>*<sup>f</sup> <sup>−</sup> *<sup>F</sup>*<sup>ˆ</sup> f, *d* <sup>f</sup> <sup>=</sup> *<sup>d</sup>*<sup>f</sup> <sup>−</sup> <sup>ˆ</sup> *d*f. Substituting Equation (22) into Equation (21) and noting *x*ˆ4 = *s*<sup>3</sup> + ζ<sup>d</sup> = *s*<sup>3</sup> + *z*<sup>3</sup> + ζ<sup>d</sup> yields:

$$
\dot{s}\_2 = \frac{RT\_{\text{sf}}a}{V\_{\text{a}}}(s\_3 + z\_3) - k\_3s\_2. \tag{30}
$$

Substituting Equation (25) into Equation (24) gives:

$$
\dot{s}\_3 = -k\_4 s\_3. \tag{31}
$$

Therefore, the time derivative of *V*<sup>1</sup> = *V*cs + *V*<sup>θ</sup> is:

$$\begin{array}{l} \dot{V}\_{1} = & -k\_{1}\mathbf{e}\_{1}^{2} - k\_{2}\mathbf{s}\_{1}^{2} - k\_{3}\mathbf{s}\_{2}^{2} - k\_{4}\mathbf{s}\_{3}^{2} + \mathbf{e}\_{1}\mathbf{s}\_{1} + \mathbf{e}\_{1}\mathbf{z}\_{1} + \overline{A}\_{\mathbf{a}}(\mathbf{s}\_{1}\mathbf{s}\_{2} + \mathbf{s}\_{1}\mathbf{z}\_{2}) + \frac{RT\_{\mathbf{a}}\mathbf{s}}{V\_{\mathbf{a}}}(\mathbf{s}\_{3} + \mathbf{z}\_{3}) \\ & + \boldsymbol{\gamma}\_{1}^{-1}\overline{b}\left[\dot{b} - (-\frac{\boldsymbol{\gamma}\_{1}}{m}\mathbf{x}\_{2}\mathbf{s}\_{1})\right] + \boldsymbol{\gamma}\_{2}^{-1}\overline{F}\_{\mathbf{f}}\left[\dot{\mathbf{f}}\_{\mathbf{f}} - (-\frac{\boldsymbol{\gamma}\_{2}}{m}\mathbf{S}\_{\mathbf{f}}(\mathbf{x}\_{2})\mathbf{s}\_{1})\right] + \boldsymbol{\gamma}\_{3}^{-1}\overline{d}\_{\mathbf{f}}\left[\dot{\hat{d}}\_{\mathbf{f}} - (-\frac{\boldsymbol{\gamma}\_{3}}{m}\mathbf{s}\_{1})\right] \\ & + s\_{1}\left\{-\frac{h\_{1}^{2}}{4\eta\_{1}}\mathbf{s}\_{1} - \frac{1}{m}\left[\overline{b}\mathbf{x}\_{2} + \overline{F}\_{\mathbf{f}}\mathbf{S}\_{\mathbf{f}}(\mathbf{x}) + \overline{d}\_{\mathbf{f}}\right] + \overline{A}\_{\mathbf{a}}\overline{x}\_{3}\right\} \end{array} \tag{32}$$

It has been proved in [36] that with the projection mapping, the following properties hold:

$$\left[\overline{b}\left[\text{Proj}\_{\mathbb{B}}(-\frac{\mathcal{V}\_{1}}{m}\mathbf{x}\_{2}\mathbf{s}\_{1})-(-\frac{\mathcal{V}\_{1}}{m}\mathbf{x}\_{2}\mathbf{s}\_{1})\right]\right]\leq 0.\tag{33a}$$

$$\left| \widetilde{F}\_{\mathbf{f}} \Big| \text{Proj}\_{\mathbb{F}\_{\mathbf{f}}} (-\frac{\mathcal{V}\_{2}}{m} S\_{\mathbf{f}}(\mathbf{x}\_{2}) \mathbf{s}\_{1}) - (-\frac{\mathcal{V}\_{2}}{m} S\_{\mathbf{f}}(\mathbf{x}\_{2}) \mathbf{s}\_{1}) \right| \leq 0. \tag{34b}$$

$$\widetilde{d}\_{\mathbf{f}} \Big[ \mathrm{Proj}\_{\hat{d}\_{\mathbf{f}}} (-\frac{\mathcal{V}3}{m} s\_1) - (-\frac{\mathcal{V}3}{m} s\_1) \Big] \le 0. \tag{33c}$$

According to Assumption 1–2, there exists a known function *h*<sup>1</sup> satisfies:

$$h\_1(t) \ge \frac{1}{m} \left[ |b\_{\text{max}} - b\_{\text{min}}| \| \mathbf{x}\_2 \| + |F\_{\text{fmax}} - F\_{\text{fmin}}| \| \mathbf{S}\_{\text{f}}(\mathbf{x}\_2) \| + |d\_{\text{max}} - d\_{\text{min}}| \right] + \overline{A}\_{\text{a}} \| \rho\_1 \|,\tag{34}$$

where <sup>ρ</sup><sup>1</sup> is the bound of *<sup>x</sup>*<sup>3</sup> as proven in the above step. From the smoothed sliding mode control theory, the following inequality holds.

$$\left\| s\_1 \left\{ -\frac{h\_1^2}{4\eta\_1} s\_1 - \frac{1}{m} \left[ \overline{b} \mathbf{x}\_2 + \widetilde{F}\_\mathbf{f} S\_\mathbf{f}(\dot{\mathbf{x}}) + \widetilde{d\_\mathbf{f}} \right] + \overline{A}\_\mathbf{a} \widetilde{\mathbf{x}}\_3 \right\} \right\| \le \eta\_1. \tag{35}$$

Substituting Equation (34) into Equation (32) yields

$$\dot{V}\_1 \le -k\_1 c\_1^2 - k\_2 s\_1^2 - k\_3 s\_2^2 - k\_4 s\_3^2 + c\_1 s\_1 + c\_1 z\_1 + \overline{A}\_\mathfrak{a} (s\_1 s\_2 + s\_1 z\_2) + \frac{RT\_\mathfrak{s} a}{V\_\mathfrak{a}} (s\mathfrak{z} + z\mathfrak{z}) + \eta\_1. \tag{36}$$

The following inequalities can be obtained by using Young's inequality.

$$\begin{array}{ll} e\_1 s\_1 \le \frac{e\_1^2}{2} + \frac{s\_1^2}{2}, e\_1 z\_1 \le \frac{e\_1^2}{2} + \frac{s\_1^2}{2}, s\_1 s\_2 \le \frac{s\_1^2}{2} + \frac{s\_2^2}{2} \\\ s\_1 z\_2 \le \frac{s\_1^2}{2} + \frac{s\_2^2}{2}, s\_2 s\_3 \le \frac{s\_2^2}{2} + \frac{s\_2^2}{2}, s\_2 z\_3 \le \frac{s\_2^2}{2} + \frac{s\_2^2}{2} \end{array} \tag{37}$$

Substituting Equation (37) into Equation (36) gives:

$$\begin{array}{ll} \dot{V}\_1 \le & -(k\_1 - 1)c\_1^2 - (k\_2 - \frac{1}{2} - \overline{A}\_a)s\_1^2 - (k\_3 - \frac{\overline{A}\_a}{2} - \frac{RT\_a a}{V\_a})s\_2^2 - (k\_4 - \frac{RT\_a a}{2V\_a})s\_3^2 + \\\ b\_2 z\_1^2 + \overline{A}\_a^2 z\_2^2 + \frac{RT\_a a}{2V\_a} z\_3^2 + \eta\_1 \end{array} \tag{38}$$

Following the similar analysis as done in [28], the following inequalities hold.

$$\dot{z}\_1 + \frac{z\_1}{\tau\_1} \le \left| \dot{z}\_1 + \frac{z\_1}{\tau\_1} \right| \le \xi\_1, \quad \dot{z}\_2 + \frac{z\_2}{\tau\_2} \le \left| \dot{z}\_2 + \frac{z\_2}{\tau\_2} \right| \le \xi\_2, \quad \dot{z}\_3 + \frac{z\_3}{\tau\_3} \le \left| \dot{z}\_3 + \frac{z\_3}{\tau\_3} \right| \le \xi\_3,\tag{39}$$

where ξ1(*e*1,*s*1, *z*1, *x*d, . *x*d, .. *x*d), ξ2(*e*1,*s*1,*s*2, *z*1, *z*2, *x*d, . *x*d, .. *x*d), and ξ2(*e*1,*s*1,*s*2,*s*3, *z*1, *z*2, *z*3, *x*d, . *x*d, .. *x*d) are three continuous functions. Due to the fact that the set <sup>Ω</sup><sup>0</sup> is compact in *<sup>R</sup>*<sup>3</sup> and the set Ω<sup>1</sup> = [*e*1,*s*1,*s*2,*s*3, *z*1, *z*2, *z*3] <sup>T</sup> : *e*<sup>2</sup> <sup>1</sup> <sup>+</sup> *<sup>s</sup>*<sup>2</sup> <sup>1</sup> <sup>+</sup> *<sup>s</sup>*<sup>2</sup> <sup>2</sup> <sup>+</sup> *<sup>s</sup>*<sup>2</sup> <sup>3</sup> <sup>+</sup> *<sup>z</sup>*<sup>2</sup> <sup>1</sup> <sup>+</sup> *<sup>z</sup>*<sup>2</sup> <sup>2</sup> <sup>+</sup> *<sup>z</sup>*<sup>2</sup> <sup>3</sup> ≤ 2*B*<sup>1</sup> (where *<sup>B</sup>*<sup>1</sup> <sup>&</sup>gt; 0) is compact in *<sup>R</sup>*7, <sup>Ω</sup><sup>0</sup> <sup>×</sup> <sup>Ω</sup><sup>1</sup> is thus compact in *<sup>R</sup>*10. Therefore, <sup>ξ</sup>1, <sup>ξ</sup>2, and <sup>ξ</sup><sup>3</sup> have maximums on <sup>Ω</sup><sup>0</sup> <sup>×</sup> <sup>Ω</sup>1. Let *<sup>M</sup>*1, *<sup>M</sup>*2, and *<sup>M</sup>*<sup>3</sup> be the maximums of ξ1, ξ2, and ξ<sup>3</sup> on Ω<sup>0</sup> × Ω1, respectively, one can obtain that:

$$z\_1 \dot{z}\_1 + \frac{z\_1^2}{\tau\_1} \le M\_1 |z\_1|, \quad z\_2 \dot{z}\_2 + \frac{z\_2^2}{\tau\_2} \le M\_2 |z\_2|, \quad z\_3 \dot{z}\_3 + \frac{z\_3^2}{\tau\_3} \le M\_3 |z\_3|. \tag{40}$$

With the use of Young's inequality, one can obtain that:

$$z\_1 \dot{z}\_1 \le -\frac{z\_1^2}{\tau\_1} + \frac{z\_1^2}{2} + \frac{M\_1^2}{2}, \quad z\_2 \dot{z}\_2 \le -\frac{z\_2^2}{\tau\_2} + \frac{z\_2^2}{2} + \frac{M\_2^2}{2}, \quad z\_3 \dot{z}\_3 \le -\frac{z\_3^2}{\tau\_3} + \frac{z\_3^2}{2} + \frac{M\_3^2}{2}.\tag{41}$$

Thus, the following inequalities can be obtained.

$$\dot{V}\_{\text{cx}} \le -\frac{1}{\tau\_1} z\_1^2 - \frac{1}{\tau\_2} z\_2^2 - \frac{1}{\tau\_3} z\_3^2 + \frac{1}{2} z\_1^2 + \frac{1}{2} M\_1^2 + \frac{1}{2} z\_2^2 + \frac{1}{2} M\_2^2 + \frac{1}{2} z\_3^2 + \frac{1}{2} M\_3^2. \tag{42}$$

Thus, the following inequalities can be obtained.

$$\begin{array}{ll}\dot{V}\_{1}\leq & -(k\_{1}-1)\varepsilon\_{1}^{2}-(k\_{2}-\frac{1}{2}-\overline{A}\_{\mathrm{a}})\mathbf{s}\_{1}^{2}-(k\_{3}-\frac{\overline{A}\_{\mathrm{a}}}{2}-\frac{RT\_{\mathrm{a}}\mathfrak{a}}{V\_{\mathrm{a}}})\mathbf{s}\_{2}^{2}-(k\_{4}-\frac{RT\_{\mathrm{a}}\mathfrak{a}}{2V\_{\mathrm{a}}})\mathbf{s}\_{3}^{2}-\\ & (\frac{1}{\tau\_{1}}-1)\boldsymbol{z}\_{1}^{2}-(\frac{1}{\tau\_{2}}-\frac{1}{2}-\frac{\overline{A}\_{\mathrm{a}}}{2})\mathbf{z}\_{2}^{2}-(\frac{1}{\tau\_{3}}-\frac{1}{2}-\frac{RT\_{\mathrm{a}}\mathfrak{a}}{2V\_{\mathrm{a}}})\mathbf{z}\_{3}^{2}+\frac{1}{2}M\_{1}^{2}+\frac{1}{2}M\_{2}^{2}+\frac{1}{2}M\_{3}^{2}+\eta\_{1}\end{array} \tag{43}$$

By choosing *k*<sup>1</sup> > 1, *k*<sup>2</sup> > <sup>1</sup> <sup>2</sup> <sup>+</sup> *<sup>A</sup>*a, *<sup>k</sup>*<sup>3</sup> <sup>&</sup>gt; *<sup>A</sup>*<sup>a</sup> <sup>2</sup> <sup>+</sup> *RT*s*<sup>a</sup> <sup>V</sup>*<sup>a</sup> , *<sup>k</sup>*<sup>4</sup> <sup>&</sup>gt; *RT*s*<sup>a</sup>* <sup>2</sup>*V*<sup>a</sup> , <sup>τ</sup><sup>1</sup> <sup>&</sup>lt; 1, <sup>τ</sup><sup>2</sup> <sup>&</sup>lt; <sup>2</sup> *A*a+1 , and τ<sup>3</sup> < <sup>2</sup>*V*<sup>a</sup> *RT*s*a*+*V*<sup>a</sup> , one has:

$$\gamma = \min \left\{ k\_1 - 1, k\_2 - \frac{1}{2} - \overline{A}\_\mathbf{a}, k\_3 - \frac{\overline{A}\_\mathbf{a}}{2} - \frac{RT\_\mathbf{s}a}{V\_\mathbf{a}}, k\_4 - \frac{RT\_\mathbf{s}a}{2V\_\mathbf{a}}, \frac{1}{\tau\_1} - 1, \frac{1}{\tau\_2} - \frac{1}{2} - \frac{\overline{A}\_\mathbf{a}}{2}, \frac{1}{\tau\_3} - \frac{1}{2} - \frac{RT\_\mathbf{s}a}{2V\_\mathbf{a}} \right\} > 0 \tag{44}$$

Then, it can be readily obtained that the following inequality holds.

$$
\dot{V} \le -2\gamma V + \eta\_\prime \tag{45}
$$

where η = <sup>1</sup> 2*M*<sup>2</sup> <sup>1</sup> <sup>+</sup> <sup>1</sup> 2*M*<sup>2</sup> <sup>2</sup> <sup>+</sup> <sup>1</sup> 2*M*<sup>2</sup> <sup>3</sup> + η1. -

Clearly, if γ>η/2*B*1, then . *V* ≤ 0 on set Ω1, thus *V* will remain *V*(*t*) ≤ *B*<sup>1</sup> for all *t*, provided that the initial conditions satisfy *V*(0) ≤ *B*1. Therefore, the errors *e*1,*s*1,*s*2,*s*3, *z*1, *z*2, *z*<sup>3</sup> are bounded. Hence, one can conclude that the closed loop system is uniformly and ultimately bounded.

#### **4. Experimental Results**

In this section, the proposed controller was implemented for the servo control of the pneumatic servo system as shown in Figure 2. The cylinder forward chamber was controlled by the valve through a 10 m long transmission line with 4 mm inside diameter, while the pressure in the return chamber was maintained at about 3 bar by utilizing a tank. A gel phantom was used to perform the following needle insertion experiments. The system physical parameters and the parameters of the controller are given in Table 1.


**Table 1.** Parameters of the system and the proposed controller.

The same control algorithm as the proposed controller, but in which the long transmission line was characterized as a part of the cylinder's dead volume and the chamber pressure was directly measured, was first tested for tracking a sinusoidal trajectory and a smooth square trajectory. As shown in Figure 4, the performance was poor for practical application, which indicates the need for an advanced method to address the issue of long transmission line. The performance of the proposed controller was tested on three types of reference trajectories: A 0.5 Hz sinusoidal signal, a smooth square signal, anda4s periodic signal. Figure 5 shows the response of the system, and the tracking errors are presented in Figure 6. The maximum absolute tracking errors (max*t*{|*e*1|}) were about 2.51 mm, 2.68 mm, and 1.81 mm, while the final steady-state tracking errors (max*t*≥5{|*e*1|}) were about 0.96 mm, 1.05 mm, and 1.39 mm. As seen, the proposed controller could significantly improve the system performance, which indicates the effectiveness of the proposed long transmission line compensation method. However, it should be noted that these needle insertion experiments were performed in a gel phantom, whose resisting force behavior was simpler than real soft tissue. The interaction between needle and soft tissue is very complex and may decrease the needle insertion precision in practical. Since the current research focused on addressing the issue of long transmission line and realizing high accuracy control of a pneumatic actuator with a pressure observer, modeling of needle-tissue interaction and further improving needle insertion precision would be the subject of the next phase of this research.

The process of parameter and disturbance estimations is shown in Figure 7. It can be seen that the estimates converged quickly, and the tracking error improvement could be achieved within several seconds. Comparisons between the observed chamber pressure and the measured one are shown in Figure 8. As can be seen, the estimated chamber pressure was close to its actual value, which indicates the effectiveness of the proposed pressure observer. However, the estimation error during the charging process was slightly bigger than the one during the discharging process. This might be due to the fact that charging process was close to adiabatic and the discharging process was close to isothermal.

**Figure 4.** Tracking response of the same control algorithm as the proposed controller but without using long transmission line compensation.

**Figure 5.** Tracking response of the proposed controller for three different types of trajectories.

**Figure 6.** Tracking errors of the proposed controller for three different types of trajectories.

**Figure 7.** Parameter and disturbance estimations of the proposed controller when the reference is a 0.5 Hz sinusoidal trajectory.

**Figure 8.** Actual and observed chamber pressure for the system tracking three different types of references.

Experiments were also conducted to verify the performance robustness of the proposed controller. To simulate a sudden disturbance acting on the system, a big step signal was added to the output of the position sensor at *t* = 22 s, and removed 10 s later. Figure 9 shows the control accuracy of the proposed controller in this scenario for tracking a sinusoidal trajectory with a frequency of 0.5 Hz and amplitude of 62.5 mm. It can be seen that the system performance did not deteriorate except the transient spikes when the sudden disturbance was added or removed.

**Figure 9.** Tracking error of the proposed controller for 0.5 Hz sinusoidal trajectory with disturbance.

#### **5. Conclusions**

In this paper, the precise motion control of a MRI compatible 1-DOF pneumatic servo system was considered. The long transmission line was characterized as an intermediate chamber connected between the valve and the cylinder in series, and a nonlinear first order system was used to approximate its dynamics. Simultaneously, a globally stable pressure observer was developed to estimate the chamber pressure. Based on the model of the long transmission line and the pressure observer, a pressure observer based adaptive dynamic surface controller was developed and the stability of the closed-loop system was proved via the Lyapunov method. In contrast to most of the existing nonlinear controllers synthesized with the backstepping method, by employing the dynamic surface control technique, the proposed controller could cope with the problem of "explosion of complexity", since the repeated differentiation of virtual controls was no longer required. The experimental results confirmed that the proposed controller was effective and had good performance robustness to sudden disturbances, thus enabling future application in pneumatically actuated MRI-compatible robots. However, precise motion control of pneumatic actuator did not necessarily lead to precise position control needle tip. Modeling of the interaction between needle and soft tissue, and incorporating a more accurate needle insertion force model in the controller design is an essential requirement for practical robot-assisted needle insertion. Thus, these issues will be explored further in the next phase of this research.

**Author Contributions:** D.M. and B.L. proposed the pressure observe based dynamic surface controller for precise position tracking control of the MRI-compatible pneumatic servo system. A.L. designed the experiments. J.Y. and Q.L. performed the experiments and analyzed the data. D.M. and A.L. wrote the paper.

**Funding:** This research was funded by [the Fundamental Research Funds for the Central Universities] grant number [2015XKMS020], [the National Natural Science Foundation of China] grant number [51505474], [the China Postdoctoral Science Foundation ] grant number [2016T90520], and [a Project Funded by the Priority Academic Program Development of Jiangsu Higher Education Institutions]. And the APC was funded by the Fundamental Research Funds for the Central Universities.

**Acknowledgments:** This work was supported by the National Natural Science Foundation of China (Grant No. 51505474), the Fundamental Research Funds for the Central Universities (Grant No. 2015XKMS020), the China Postdoctoral Science Foundation (Grant No. 2016T90520), and a Project Funded by the Priority Academic Program Development of Jiangsu Higher Education Institutions.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## *Article* **Tip Estimation Method in Phantoms for Curved Needle Using 2D Transverse Ultrasound Images**

#### **Zihao Li 1, Shuang Song 1,\*, Li Liu 2,\* and Max Q.-H. Meng <sup>2</sup>**


Received: 8 October 2019; Accepted: 27 November 2019; Published: 5 December 2019

**Abstract:** Flexible needles have been widely used in minimally invasive surgeries, especially in percutaneous interventions. Among the interventions, tip position of the curved needle is very important, since it directly affects the success of the surgeries. In this paper, we present a method to estimate the tip position of a long-curved needle by using 2D transverse ultrasound images from a robotic ultrasound system. Ultrasound is first used to detect the cross section of long-flexible needle. A new imaging approach is proposed based on the selection of numbers of pixels with a higher gray level, which can directly remove the lower gray level to highlight the needle. After that, the needle shape tracking method is proposed by combining the image processing with the Kalman filter by using 3D needle positions, which develop a robust needle tracking procedure from 1 mm to 8 mm scan intervals. Shape reconstruction is then achieved using the curve fitting method. Finally, the needle tip position is estimated based on the curve fitting result. Experimental results showed that the estimation error of tip position is less than 1 mm within 4 mm scan intervals. The advantage of the proposed method is that the shape and tip position can be estimated through scanning the needle's cross sections at intervals along the direction of needle insertion without detecting the tip.

**Keywords:** needle tracking; tip estimation; transverse ultrasound images; percutaneous interventions

#### **1. Introduction**

With the help of the beveled-tip needle, percutaneous interventions and therapies have been widely involved in current clinical procedures, such as brachytherapy [1,2], tissue biopsy [3,4], and drug delivery [5,6]. In intervention procedures, less needle misplacement will lead to a more reliable treatment and a more accurate medical practice. According to the clinical studies [7,8], the needle is easy to be deflected, which will cause needle tip misplacement and may lead to unsafe procedures. Due to the needle-tissue interaction, improper insertion force or physiological motions, such as breathing may cause targets or obstacles to be unstable, which will lead to an unexpected error. To address the challenge, real-time feedback is highly required. Usually, medical imaging devices are used, such as ultrasound (US) [9], computerized tomography (CT) [10,11], or magnetic resonance imaging (MRI) [2,12]. Generally, the image-guided percutaneous interventions are conducted with the use of CT or MRI. However, ultrasound-guided procedures are more attractive due to their advantages such as none ionizing radiations and real-time detection.

Many studies for the guidance of the needle during the insertion operation have been conducted with the help of ultrasound devices, and 2D ultrasound images are quite general to use, especially for the sagittal one (shown in Figure 1). Elif et al. proposed to use circular Hough transform to locate the needle tip accurately, even when the imaging is out-of-plane [13]. Kaya et al. localized the needle axis and estimated the needle tip by using a Gabor Filter in sagittal US images [14]. To execute in real-time, they improved the processing time by applying the bin packing method [15]. Recently, a template-based tracking method with the efficient second-order minimization optimization method has been used to track the needle [16]. In recent studies, more and more novel ideas have been used to locate the needle and evaluate its tip to sagittal US images, such as the use of signal attenuation maps [17], convolution neural networks (CNN) [18], and maximum likelihood estimation sample consensus (MLESAC) method [19]. However, a demerit of using sagittal US images is that out-of-plane bending of the needle cannot be detected. Therefore, the methods applied on sagittal US images are not suitable for the needle which may be deflected by the inevitable factors, especially for long needles.

**Figure 1.** Two methods by using 2D ultrasound for detection.

An alternate solution for this problem is to use a 3D US image, which has been widely studied in recent researches. Yue et al. used a RANSAC method to detect the needle in a 3D US situation and Kalman filter has been used to reduce the error [20]. Chatelain et al. used the particle filtering to track a robot-guided flexible needle by using 3D US [21]. In addition, a convolutional neural network with conventional image processing techniques has also been used to track and detect the needle [22] and a naive Bayesian classifier was used to localize the needle among 3D US volume voxels [23]. However, the large 3D US volumetric dataset would make it difficult to obtain and process the online data.

Due to the above disadvantages, sagittal US images and 3D US volume are not suitable for a long flexible needle. To locate the needle accurately, methods that use transverse US images (shown in Figure 1) have been used successfully in some studies. For example, Vrooijink et al. [24] present a method to track the flexible needle during the insertion into a gelatin tissue by using 2D US images perpendicular to its needle tip. However, the background is pure without noise, which makes it impractical. Waine et al. [25–27] focus on the research about the needle insertion, in permanent prostate brachytherapy (PPB), where needles are typically 200 mm and easily to be deflected, indicating the fact that the rectum limits the movement of US probe. As a result, it is hard to acquire the sagittal images of the curved needle to observe the deflection during needle insertion. This is because the sagittal method has a strong relationship to the movement of the US probe when the needle is out of the view-field of the US images, and this movement maybe deforms the prostate as well as affects the needle target. For a deflected needle, the transverse US image is a better choice for its detection. Compared to the sagittal images, the transverse US images are easy to be acquired when the US probe scanning along the needle, no matter how much the needle is curved.

In this paper, we present a method to track a long-curved needle from the 2D transverse US images and estimate its tip for the guidance of needle insertion. Ultrasound is first used to detect the cross-sections of the long-flexible needle (shown in Figure 2 STEP 1). The needle shape tracking method combined needle detection with Kalman filter develops an accurate location and a robust tracking procedure with scan intervals from 1 mm to 8 mm (shown in Figure 2 STEP 2). Unlike the previous study [26], the 3D needle positions obtained from 2D US images and optical tracking systems have been used in KF for the precise location. The curve fitting method is then used to achieve the shape reconstruction and the needle tip position is estimated based on its length and the curve fitting result (shown in Figure 2 STEP 3). The advantage of the proposed method is that the shape and tip position can be estimated through scanning the needle's cross-sections at intervals along the direction of needle insertion without detecting the tip. Besides, a novel histogram method is introduced to detect the needle in image processing, which can improve the needle localization under the effect of needle comet tail and the poor reflection, despite of the abrupt intensity changes. In addition, a robotic ultrasound system (RUS) [28] is built to evaluate the proposed needle tip estimation method. Results showed that the estimation of the tip position is less than 1 mm with 4 mm scan intervals.

**Figure 2.** The proposed tip estimation method. STEP 1: 2D transverse US images with needle cross-sections are collected by using RUS; STEP 2: Needle cross-sections are detected and tracked in the successive US images; STEP 3: Needle shape is constructed and its tip is estimated.

The rest of this paper is organized as follows. The proposed methods will be introduced in detail in Section 2. Section 3 intends to represent the experimental setup and the results. Finally, the discussion and conclusions are detailed in Section 4.

#### **2. Materials and Methods**

The proposed needle tip estimation method in successive transverse US images can be divided into three stages: needle detection, needle shape tracking, and tip estimation. The processing diagram are shown in the Figure 3. The needle location will be manually selected as an initial region of interest (ROI) by the binary method in the first US image. After that, the prediction of needle position from a Kalman filter can be transformed in the transverse US images. At the same time, the next needle position in this ROI can be found through the histogram method. The KF is then updated for the current precise needle position and prediction of the next position. After all the cross-sections of the needle have been collected,

the needle shape can be fitted by the curve fitting method (part C in Figure 3). Finally, the position of the needle tip can be estimated from the curve fitting result based on the length of the needle. In this study, the ROI is set as a square window with a length of three times larger than the needle diameter and its center represents the needle position. Needle detection (part B in Figure 3) is mainly about image processing, while the Kalman filter is used for needle shape tracking (part A in Figure 3).

**Figure 3.** The pipeline of curved needle tip estimation. As the ROI configuration finished, needle tracking with needle detection begins step by step. Part A shows the needle tracking by KF, while Part B shows the needle detection to locate the needle tip. Part C shows the shape construction and tip estimation.

#### *2.1. Needle Detection*

Needle detection is mainly for identifying the cross-section of the needle in the US images by using the binary method and the histogram method. As the ultrasound is quite sensitive to the metal, the needle-inside area can be brighter than others. A binary method [26] is first used to select ROI in

the first image and then a histogram method is used to locate the needle despite the US shadows and poor reflection.

Binary method intends to strengthen the contrast of brightness to highlight the brighter area to select them. This method is used for the initialization which supposes to find the candidates in the first image. It includes intensity normalization, background reinforcement, and brightness enhancement. The center of the area is the location of the needle. During the experiment, a histogram method is proposed to find the needle accurately. The histogram method contains intensity normalization and background reinforcement. The histogram method tends to find an area of high-intensity pixels, which intends to find the upper face of the needle and then locate the needle based on the diameter. The background reinforcement part can be described as follows:

$$\begin{aligned} \min & \quad & I\_t \\ \text{s.t.} & \quad & \sum\_{I\_t}^{255} n\_{I\_t} \le \mathcal{M}\delta \\ & \quad & I\_t \in [0, 255] \end{aligned} \tag{1}$$

where *Ij* is the gray level of the pixel and *nIj* is the number of pixels which have the gray level of *Ij*. *M* is the size of ROI, and *δ* is the manually selected parameter to limit the bright pixels. In this work, *M* is a square of 45 × 45, and *δ* is set to 0.08 based on empirical tests.

There are two unexpected situations that may affect the position accuracy, namely the comet tail and the poor reflection. The comet tail will affect the size of the needle area and usually lead to a larger area than the actual size (Figure 4a). On the contrary, the poor reflection makes the needle area look much smaller in the image (Figure 4b). Therefore, the accurate location should be intended to eliminate the effect of shadows and poor reflection. In Figure 4, there are two examples which are used the two methods relatively. As the example shown in Figure 4a, the noise (yellow circle in the histogram of ROI) may probably be concerned as the needle while the needle is just concerned about a few pixels (red circle in histogram of ROI). The two methods can both filter the noise and locate the needle correctly. However, in ROI configuration (shown in Figure 3), the histogram method intends to find more candidates than the binary method, since the former focuses on the higher intensity pixels while the latter focuses on the area and intensity. Therefore, the binary method is more feasible in ROI configuration.

However, in the case of poor reflection, the needle displays a little in the image, and the area of the needle is smaller than the expected. Because the needle would reflect as long as the image gain is high enough or the sound power is big enough, it would reveal apparently compared to its surroundings. An example of the histogram of ROI is shown in Figure 4b, the red circle represents the upper surface of the needle. Moreover, the ROI square is darker than the one in Figure 4a, while the settings of the ultrasound are the same in Figure 4. From the figure, the histogram method seems to be more accurate than the binary method. In fact, the error of two methods, in this case, is 0.34 mm with 1.2 mm diameter of the needle. It is not that obvious to judge the accuracy. In this study, we use the histogram method during the experiment.

(**b**) Poor reflection.

**Figure 4.** Two cases may generate the errors: (**a**) the comet tail of needle with binary method process, histogram method process and the histogram of ROI; (**b**) the poor reflection of the needle with histogram method process and the histogram of ROI.

#### *2.2. Needle Shape Tracking*

As indicated in previous researches [20,29,30], Kalman filter has been successfully used for tracking needle in the successive ultrasound images. In this study, the Kalman filter is used to improve the estimation of the needle location in successive frames. As shown in Figure 5, the applied Kalman filter has two processes, prediction and update. The prediction stage intends to locate the needle position previously and set the ROI (red and yellow square in Figure 5) to find the needle precisely with a small window, which is supposed to reduce the computation. The update stage is the result of the needle position after the measurement from the histogram method.

**Figure 5.** The two steps of the Kalman filter. As the next US image is acquired, the previous state (*xprevious*, *yprevious*, *zprevious*, *xprevious*, *yprevious*, *zprevious*) is used to predict the needle position (*xprediction*, *yprediction*, *zprediction*, which then will be transformed to (*Ix*\_*prediction*, *Iy*\_*prediction*) in the image as the center of ROI. The yellow square is the ROI corresponding to (*Ix*\_*prediction*, *Iy*\_*prediction*) and the red one is the update step in KF by using the measurement data from needle detection to locate the needle with its center (*Ix*\_*update*, *Iy*\_*update*) as the needle position. Finally, the measurement state (*xm*, *ym*, *zm*, *xm*, *ym*, *zm*) and the current state (*xc*, *yc*, *zc*, *xc*, *yc*, *zc*) can be obtained.

The state prediction ˆ*ti* intends to represent the prediction state of the transverse needle center position (*x*, *y*, *z*) in the reference frame with the change of the needle position (*x*, *y*, *z*) at sample *i* according to the state t. (*x*, *y*, *z*) are the difference between the previous needle position and the current needle position. *ti* is the result from the previous iteration, which is as follows:

$$\begin{aligned} t\_i &= \begin{bmatrix} x\_i \\ y\_i \\ z\_i \\ \triangle x\_i \\ \triangle y\_i \\ \triangle z\_i \end{bmatrix} \end{aligned} \tag{2}$$

where *x*<sup>1</sup> and *y*<sup>1</sup> are set to be 0, while *z*<sup>1</sup> is equal to the scan interval. The prediction equations are as follows:

$$
\hat{t}\_i = At\_{i-1} \tag{3}
$$

$$
\dot{P}\_i = A P\_{i-1} A^T + Q.\tag{4}
$$

The measurement update equations are as follows:

$$\mathcal{K}\_{l} = \mathring{P}\_{l} \boldsymbol{H}^{T} \left( \boldsymbol{H} \mathring{P}\_{l} \boldsymbol{H}^{T} + \boldsymbol{R} \right)^{-1},\tag{5}$$

$$t\_i = \hat{t}\_i + \mathcal{K}\_i (m\_i - H\hat{t}\_i),\tag{6}$$

$$P\_l = (I - K\_l H)\dot{P}\_{l\prime} \tag{7}$$

where A, H, R, and Q are as follows:

$$A = \begin{bmatrix} 1 & 0 & 0 & 1 & 0 & 0 \\ 0 & 1 & 0 & 0 & 1 & 0 \\ 0 & 0 & 1 & 0 & 0 & 1 \\ 0 & 0 & 0 & 1 & 0 & 0 \\ 0 & 0 & 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 0 & 0 & 1 \end{bmatrix},\tag{8}$$

$$H = \begin{bmatrix} I\_{6 \times 6} \end{bmatrix}, R = Q = 10^{-6} \times \left[ I\_{6 \times 6} \right].\tag{9}$$

A is the state transition matrix, H is the measurement matrix. *P*ˆ *<sup>i</sup>* and *Pi* are the priori and posteriori estimate error covariance, and R and Q are the measurement error covariance and processing error covariance, respectively. *Ki* is the Kalman gain at sample i. *mi* is the measurement state from the needle detection.

The 3D prediction position (*xprediction*, *yprediction*, *zprediction*) is obtained from the previous state. Before needle detection in the current US image, the 3D prediction position should be transformed on the image plane frame as 2D position (*Ix*\_*prediction*, *Iy*\_*prediction*). After the update, the needle position (*Ix*\_*update*, *Iy*\_*update*) in the image will be transformed into the 3D position (*xm*, *ym*, *zm*). Meanwhile, we can get (*xm*, *ym*, *zm*) from:

$$(\triangle x\_m \triangle y\_m, \triangle z\_m) = (x\_m, y\_m, z\_m) - (x\_{previous}, y\_{previous}, z\_{previous}).\tag{10}$$

As a result, the measurement state *mi* in this frame is acquired as(*xm*, *ym*, *zm*, *xm*, *ym*, *zm*). Through the measurement update, the current state can be obtained as (*xc*, *yc*, *zc*, *xc*, *yc*, *zc*) from Equation (6). The relationship between these transformations will be described in the next subsection. In the previous study [26], the data from the image has only two dimensions, lacking the data from the direction along the movement of the probe, which leads to an incomplete location. Moreover, the space information is more capable to locate the needle accurately than the plane information. Therefore, 3D positions have been used in the KF for the precise location. The KF in this work is not only used for filtering, but also for predicting the next needle position in the US image. The ROI for the next iteration is centered around the needle position of the Kalman Filtering prediction, which can help to remove the outliers from the ROI.

#### *2.3. Tip Estimation*

Before tip estimation, 2D points should be transformed into 3D points based on the relationship of each frame. The relationship among the reference frame, probe frame, marker frame, and image frame are specified in Figure 6.

**Figure 6.** The relationship of the frames.

As shown in Figure 6, the image has one plane with 2 axes (axis x and axis y) and axis z is vertical to the image. Moreover, the probe has the same frame with image, except that the probe frame is designed in millimeters and the image frame is set in pixels. Equation (9) implies the transformation from image to reference:

$$Point1 = T\_{marker}^{ref} \times T\_{probe}^{marker} \times T\_{image}^{probe} \times Point2 \tag{11}$$

where *Point*<sup>1</sup> and *Point*<sup>2</sup> are the points on the reference frame and image frame, respectively, *<sup>T</sup>ref marker* is the transformation from the reference frame to the marker frame, *Tmarker probe* is the transformation from the marker frame to the probe frame, *Tprobe image* is the transformation from the probe frame to the image frame. Through this transformation, the needle position in the image can be directly re-defined in the reference

frame for the needle tracking and curve fitting.

The tip estimation has two steps: curve fitting and tip estimation. In this study, the third-order curve line is used for the shape construction where the equations can be written as:

$$f(\mathbf{x}) = \sum\_{k=0}^{3} a\_k \mathbf{x}^k. \tag{12}$$

$$\log(\mathbf{x}) = \sum\_{k=0}^{3} b\_k \mathbf{x}^k,\tag{13}$$

where *f*(*x*) and *g*(*x*) are the equations to fit the line along the x, which is the axis with the same direction of the insertion. (*a*0, *a*1, *a*2, *a*3) and (*b*0, *b*1, *b*2, *b*3) are the free parameters of the needle shape model.

After sample points of the inflected needle have been obtained, the least-square curve fitting method will be used to fit these points as a cubic line. The target functions to fit the cubic line can be defined as follows:

$$F(a\_0, a\_1, a\_2, a\_3) = \min \sum\_{i=1}^n (f(\mathbf{x}\_i) - y\_i)^2 = \min \sum\_{i=1}^n (\sum\_{k=0}^3 a\_k \mathbf{x}\_i^k - y\_i)^2 \tag{14}$$

$$G(b\_0, b\_1, b\_2, b\_3) = \min \sum\_{i=1}^{n} (g(\mathbf{x}\_i) - z\_i)^2 = \min \sum\_{i=1}^{n} (\sum\_{k=0}^{3} b\_k \mathbf{x}\_i^k - z\_i)^2. \tag{15}$$

where *n* is the number of the points (*n* ≥ 4) and (*xi*, *yi*, *zi*) is the position of point *i*. By applying the *l*<sup>2</sup> norm minimization in the two-dimensional Euclidean space, it can be formulated as:

$$\underset{t \in R}{\text{arg min}} \; \|X\mathfrak{a} - Y\|\_{2^\prime} \tag{16}$$

where *a* = [*a*0, *a*1, *a*2, *a*3] , *Y* = [*y*1, *y*2,..., *yn*] and *X* can be written as:

$$X = \begin{bmatrix} 1 & \mathbf{x}\_1 & \mathbf{x}\_1^2 & \mathbf{x}\_1^3 \\ 1 & \mathbf{x}\_2 & \mathbf{x}\_2^2 & \mathbf{x}\_2^3 \\ \vdots & \vdots & \vdots & \vdots \\ 1 & \mathbf{x}\_n & \mathbf{x}\_n^2 & \mathbf{x}\_n^3 \end{bmatrix}. \tag{17}$$

The solution can be estimated as follows:

$$\mathfrak{a} = \left(\boldsymbol{X}^T\boldsymbol{X}\right)^{-1}\boldsymbol{X}^T\boldsymbol{Y}.\tag{18}$$

From this solution, *f*(*x*) and *g*(*x*) can be obtained to construct needle shape. The tip position can then be estimated by the following optimum solutions based on the length of the needle:

$$\begin{aligned} \min & \quad \| \int\_{tail\_x}^{tip\_x} \sqrt{1 + f'(\mathbf{x})^2 + \mathbf{g}'(\mathbf{x})^2} d\mathbf{x} - L \|\_2 \\ \text{s.t.} & \quad tip\_x > tail\_x \end{aligned} \tag{19}$$

where *tailx* is the measured value of the tail position from the optical tracker, *tipx* is the expected value of the tip position in axis *x*, *L* is the length of the needle.

#### **3. Results**

#### *3.1. Experimental Platform Setup*

To verify the proposed tip tracking and shape sensing method, a robotic ultrasound system has been built, which includes a KUKA IIWA robot arm, a Wisonic ultrasound scanner, an NDI optical tracker, an NDI electromagnetic (EM) tracker, and a computer. As shown in Figure 7, the US probe is mounted on the effector of the robot arm by the gripper attached to the passive marker. The phantom or ex-vivo (like chicken breast in Figure 7) is punctured by an 18G beveled-tip needle with 200 mm long, while the needle tip is completely exposed for validation. The diameter of the needle is 15 pixels in the image. The NDI optical tracker is used to localize the marker bound with the probe, while NDI electromagnetic tracker is used to validate the tip position.

Experiments have been taken in a water tank, which provides a liquid environment for the ultrasound. And the needle is placed in water or inserted in the silica gel phantom (shown in Figure 8a), pork and chicken breast. The depth of the ultrasound is set to 4 cm. In this study, the needle is usually detected in 1 to 3 cm from the US probe. During the experiment, the robot arm automatically moves with the ultrasound probe along the direction of the needle insertion without any contact with the tissue (shown in Figure 8b). The whole scan length is at most 160 mm which depends on the scan intervals (shown in Table 1). The scan interval decreases with the increasing collected points.


**Table 1.** Scan lengths with different scan intervals.

**Figure 7.** The devices of the experiment.

Before data collection, the US image and the marker need to be calibrated. After that, the experiment starts after the needle finished puncturing manually. The robot arm is used to move the US probe scanning along the needle. Meanwhile, pose data are collected from the optical tracker and US images from the ultrasound scanner. Finally, the tail of the needle and its tip are measured by the optical and

electromagnetic sensors, respectively, for the curve fitting and the tip validation. The needle is inserted manually, imitating the real situation of percutaneous intervention.

(**a**) Phantom with inserted needle (**b**) The probe scans without any contact

**Figure 8.** The phantom used in the experiment and the movement of the probe.

#### *3.2. Tip Estimation*

Four kinds of platforms have been used in the experiments: water, phantom, chicken, and pork. Each platform was tested several times. Figure 9 shows one test in chicken breast. In this case, the US probe moved along the needle in the chicken breast every 4 mm. The black square point on the left is the needle tail position and the blue line is the estimated needle shape. The green points are the detected needle and considered as the center of the needle, the blue point is the estimated tip position and the red point is measured tip position from EM. The estimation error is 0.69 mm in this test.

**Figure 9.** Experiment in chicken breast with a 4 mm scan interval.

The error of the algorithm is shown in Table 2, which suggests that the errors increase with the increase of scan intervals. Figure 10 shows the results of the experiment on four platforms. The mean errors are all under 0.4 mm with a 1 mm scan interval in the four experiments, while the error is around 1 mm with an 8 mm interval.

**Figure 10.** The experiment with different scan intervals on the four platforms.


**Table 2.** The results of tip estimation (mm).

#### **4. Discussion and Conclusions**

Needle insertion guided by ultrasound images is widely used for percutaneous interventions. However, the needle detection due to its deflection by the inevitable factors is challenging during the needle insertion. Such factors include needle-tissue interaction, improper insertion force, physiological motions, and so on. Automatic needle detection with needle tracking in 2D transverse US images could overcome these limitations and estimate needle tip through the curve-fitting method. The target of this study is to develop a robust needle detection and tracking method with the help of ultrasound images to estimate the needle tip precisely and accurately. We used a histogram method to detect the needle in transverse US images to decrease the effects of comet tail and poor reflection. In subsequent post-processing, the needle was tracked by the Kalman filter tracking method in consecutive US images with the help of the displacement of the probe. A third-order curve fitting method has been used to estimate the needle tip. When the probe is moved by the robot arm, the scanning time is different. We assume that the time when the probe stops to collect the data is the same. The less scan interval we choose, the more collection points we can obtain and the more scanning time it takes. Therefore, the scanning time

mainly depends on the number of scanning points while the accuracy lies in how short the scan interval is. In other words, the accuracy of the tip estimation can be improved by reducing the scan interval to collect more needle positions. However, this will consume more scanning time and reduce the efficiency of tracking. Inversely, fewer collecting positions would cost less time but may lead to a more possibility of the failure of shape construction and a more possibility of large error of the tip estimation. As a result, how to balance precision and scanning time is very important to make the proposed method more efficient. In our experiment, it is found that a 4 mm scan interval has an error less than 1 mm, which is a better choice to satisfy both requirements.

In the proposed method, needle shape tracking has a great contribution to the accurate localization, since the needle can be tracked precisely by Kalman filter through its prediction and update. However, needle shape tracking is heavily dependent on the scan interval, especially for a large curved needle, since Kalman filter is generally well functioned in the lineal system. As a result, if the needle is deflected during insertion, the Kalman filter would make mistakes and wrongly predict the needle position when the scan interval is large. In this study, it is found that Kalman filter would fail if the scan interval is more than 8 mm. This may due to the impact on the prediction of KF with a large deviation. Moreover, the deviation will not be eliminated even with the change of the ROI size. However, the Histogram method showed the accurate and effective detection of the needle, but it relies on the brightness of the image as the needle could not be easily detected where there are plenty of pixels with the highest intensity (which has a max value of 255). However, this condition can be controlled by the setting of the ultrasound scanner to expand the gray level of the image properly.

The proposed method still has its limitations. During the experiment, time is needed for the data collection from the US scanner and optical tracking system, and the movement of the robot arm, which is affected by the scan length and scan interval. However, it is very hard to acquire the whole position of a long needle in one scan for any medical image sensor. Therefore, in the future, it is valuable to find a method to reduce the times of needle detection in order to reduce the time for the tip estimation. Moreover, when the robot arm moves with the US probe precisely, the tissue and needle have a possibility to be deformed by the probe motion. Hence, it is necessary to make the robot arm move smoothly as well as correctly on the surface of tissue. In addition, patient motion is the biggest uncertain problem, which leads to the failure of needle insertion and detection.

In the proposed system, we use the 2D US scanner for the detection of the needle in various kinds of tissue. However, the 3D US scanner can also be used in this system. Although it has volume data and the detection method is different, the tracking method is able to be similar, as we also use the 3D positions for tracking in this study. Moreover, time is also needed for data collection and the movement of the robot arm, especially for the long needle, which is easily out of view-field of US volumes or images. Therefore, we use 2D US images in this study.

In this paper, a method for tracking a long-curved needle from the 2D transverse US images and the tip estimation is represented and demonstrated with RUS. Ultrasound is first used to detect the cross section, with the probe moving along a long-flexible needle. Needle shape tracking method combined needle detection with Kalman filter by using 3D needle positions develops an accurate location and a robust tracking procedure. Needle shape is then constructed by using the curve fitting method and its tip position is estimated based on the former result. A histogram method is introduced to detect the needle in image processing, which can improve the needle localization despite the abrupt intensity changes. This new imaging approach is proposed based on the selection of numbers of pixels with a higher gray level, which can directly remove the lower gray level to highlight the needle. Results of the experiments suggest that the detection of the needle by the histogram method and Kalman Filter has high precision with minimum error 0.13 mm with a 1 mm scan interval in the phantom experiment and maximum error 1.35 mm with a 8 mm scan interval in the pork experiment. With the increase of the scan interval, the mean

error would rise. Moreover, results showed that the estimation of the tip position is less than 1 mm within 4 mm scan intervals. We suggest choosing a 4 mm scan interval to balance the precision and scanning time to maximize efficiency. In the future, we would make the experiments of how long the scan length is the best length to estimate the needle tip. The proposed method would be a great assist to surgeons to locate the needle tip when they perform percutaneous insertion procedures with a long flexible needle, such as prostate brachytherapy.

**Author Contributions:** Conceptualization, M.Q.-H.M.; validation, Z.L.; writing—original draft, Z.L.; writing—review and editing, S.S. and L.L.

**Funding:** This research was funded in part by the National Natural Science Foundation of China, grant number 61803123, and in part by the National Key R&D Program of China, grant number 2018YFB1307700, and in part by the Natural Science Foundation of Guangdong Province, China, grant number 2018A030310565.

**Acknowledgments:** For this kind of study, no formal ethics approval is required by the institutional ethics com. Informed consent was obtained from all individual participants included in the study.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **Abbreviations**

The following abbreviations are used in this manuscript:


#### **References**




c 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

#### *Article*
