Intraoperative Data-Based Haptic Feedback for Arthroscopic Partial Meniscectomy Punch Simulation

The use of intraoperative tool-tissue interaction data for arthroscopic digital twins has the potential to enhance haptic feedback in surgical simulations. This study demonstrates the implementation of an interpolation-based method, and a Kalman filter-based method for using intraoperative interaction data to provide nonlinear haptic feedback in a partial meniscectomy punch simulation. The interpolation-based method was implemented by interpolation in tool-tissue interaction data collected from a cadaveric meniscus. The Kalman filter-based method takes a state-space formulation of the arthroscopic punch force in the one-dimensional space domain, and performs sensor fusion of the interpolation-based force signal, and a reference finite element force, to form a new haptic signal. These methods were demonstrated using a novel inexpensive haptic system for simulating an arthroscopic partial meniscectomy punch. The software, computer-aided design models, and collected data are provided open-source on GitHub. The face validity of the methods was evaluated in user studies, including six experienced orthopedic surgeons, and five non-professionals. The findings showed that experienced surgeons could distinguish between a nonlinear interpolation-based force signal, and a linear-elastic ideally-plastic finite element-based reference signal in a partial meniscectomy punch simulation. Surgeons preferred the interpolation method as a haptic rendering strategy in this study. The Kalman filter method was not as effective as interpolation in recreating nonlinear haptic feedback for the tested parameters, but demonstrated coupling of haptic force signals from intraoperative data and an idealized finite element method in a partial meniscectomy punch simulation.


I. INTRODUCTION
Surgical simulation of arthroscopic procedures enables safe surgical training without causing harm to the patient [1]. A lower volume of some types of arthroscopic surgeries, as well as a high demand for specialized competence, has led to a gap in providing the training needed to go from student The associate editor coordinating the review of this manuscript and approving it for publication was Lei Wei . to specialist. Arthroscopic partial meniscectomy (APM) is an example of an arthroscopic surgical knee procedure which is now less frequently performed. A recent report by Jacquet et al. [2] stated that among 461 orthopedic surgeons, 75% reported to perform more meniscus repairs and 85% less APM than 5 years ago. This, combined with advancements in simulation technology, has led to a paradigm shift from apprenticeship to simulation-based methods in training of resident doctors [3], [4]. Several arthroscopic knee simulators VOLUME 10, 2022 This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ have been developed for this purpose, as pointed out in recent reviews [5], [6]. These typically include standardized training procedures on generic anatomy. However, patient-specific simulation has been identified as an important future direction for surgical simulation systems, and Ryu et al. [7] found that patient-specific simulators have the ability to develop higher-level competencies, for example for experienced surgeons, outside a clinical setting. In a review from 2019, Overtoom et al. [8] assessed the value of haptic and force feedback in surgical simulators. The findings showed that haptic feedback leads to a shorter learning curve, and that surgeons consider haptic feedback as more important than visual feedback in surgical simulators. In our previous work, we have defined a digital twin concept for arthroscopic knee surgery as ''virtual information that fully describes a patient-specific biomechanical system, such as a joint'' [9]. Here, the digital twin expands on conventional surgery simulators by including patient-specific digital information with real-time calibration of simulations using intraoperative data. Intraoperative data can include intraoperative imaging modalities, such as X-ray, ultrasound or arthroscopic images, as well as live tool-tissue interaction data, as shown in Fig. 1. Three use cases were proposed: (i) surgical skills training for resident doctors, (ii) patient specific preoperative planning, and (iii) a post-operative database of rare surgical cases for advanced training and planning. As such, the digital twin fidelity should increase throughout the three scenarios as complexity increases and surgeons are more experienced. The use of live tool-tissue interaction data has the potential to increase face validity in all three scenarios by improving haptic feedback. We adopt the definition by Vaughan [6] that ''face validity measures how a simulator appears similar to the actual procedure''. In the third use case, we argue that the use of intraoperative tool-tissue interaction data is essential to exploit the full potential of the digital twin. To this end, we define the following research question: How can intraoperative tool-tissue interaction data contribute to better estimation of haptic feedback in arthroscopic surgical simulation?
The concept of recording and modelling mechanical interactions in haptic feedback has previously been referred to as the ''haptic camera'' [10], ''reality-based modeling'' [11], and ''haptography'' [12]. For surgical simulation, Greenish et al. [13] demonstrated the collection of cutting data in anatomical tissues for surgical scissors, rendering the cutting [14], and evaluating the performance. In this study, cutting data were collected by equipping the surgical scissors with strain gauges. Similarly, Callaghan et al. [15] described a setup for collecting cutting data from biological tissue using strain gauges and a potentiometer. More advanced setups have also been described, such as equipping a custom arthroscopic grasper with fiber Bragg sensors [16], and a surgical tool with a proximally located force-torque sensor [17].
Okamura [11] described three methods for active haptic rendering of recorded interaction data in the setting of surgical scissors. These were (i) directly using interpolation in a lookup-table of recorded data, (ii) developing a piecewise linear transfer function with fitted parameters, or (iii) developing an analytical mathematical model of the relevant tool and tissue, such as tool friction and tissue stiffness. The third method was later further developed by Mavash et al. [18]. A drawback of these methods is that they do not account for visualization, such as in mesh-based biomechanical soft tissue deformation and cutting methods.
Currently, the gold standard for patient-specific soft tissue and cutting simulations is to calibrate the material properties for simulations based on the finite element method (FEM) from interaction or experimental data. Although some patient-specific imaging methods, such as shear wave elastography [19] and magnetic resonance elastography [20], have the potential to provide material stiffness data for soft tissues, they have not been sufficiently developed. Lauzeral et al. [21] used model order reduction and statistical shape analysis from computer tomography (CT) images to create a patient specific model of a human liver, and modeled breathing motion. Using experimental data, Bojairami et al. [22] recently extracted tissue properties for probe insertion in real-time using cohesive elements, and Seyfi et al. [23] identified the mechanical properties of the human meniscus using an inverse FEM. Similarly, Wu et al. [24] investigated how live data obtained during a robotic endoscopic procedure could be used to correct imprecise FEM results using machine learning. However, the challenge of relying on finite element (FE) models for modeling haptic interactions experienced in surgery is that it is difficult to accurately model the contact forces involved with topology changes, such as in arthroscopic partial meniscectomy. Peterlik et al. [25] presented the constraint-based multi-rate compliant mechanisms framework based on Lagrange multipliers, which provided more flexibility for modelling contact interactions in haptic simulations, but did not consider topology change. Bui et al. [26] later presented the corotational cut finite element method, which models topology change and associated cutting forces, however, the use of intraoperative data for patient-specific calibration of these parameters has not yet been demonstrated. Schulmann et al. [27] recently studied the effect of discretization on parameter identification for patient-specific applications, and found that ''estimated parameters should not be considered as the true parameter value of the organ or tissue, but instead are model-dependent values''. For a true patient-specific simulation, mechanical properties influenced not only by organ or tissue, but also by variations in age, sex and pathology should be included, and fast and accurate model generation is therefore important. As an alternative approach, we argue that using intraoperative data directly in combination with a generic FE-model can be used to provide patient specific haptic feedback in surgical simulations.
This paper explores two methods for enabling patientspecific nonlinear haptic feedback from intraoperative tooltissue interaction data in arthroscopic partial meniscectomy simulations. The first method is based on interpolation of FIGURE 1. A new illustration of the arthroscopic digital twin concept described in [9]. The left side shows clinical input data throughout pre-, intra-and post-operative stages. The right side shows how the clinical input data relates to the output surgical user scenarios. The focus of this paper is exploring methods for using intraoperative tool-tissue interaction data for haptic feedback in a post-operative virtual surgery database.
intraoperative data, as previously demonstrated by Okamura [11]. The second method is a novel application of the well-known Kalman filter (KF), and works by performing sensor fusion of intraoperative and FEM haptic signal data to form a new haptic signal. The mathematical statespace representation of the KF-based haptic arthroscopic punch force signal is formulated in the one-dimensional space domain, with force and stiffness selected as states.
To the best of our knowledge, the Kalman filter has not been explored for this application. The advantage of these methods lies in enabling patient-specific nonlinear haptic feedback while maintaining an idealized FEM as a baseline for surgical simulations, thus limiting the complexity involved in modelling cutting forces using FEM. The methods are explored by implementing a one degree-offreedom (DOF) open-source haptic system simulating a partial meniscectomy punch, based on tool-tissue interaction data. Implementation of a simple FE-model based on 1D spring elements is included with the purpose of generating an idealized linear-elastic ideally-plastic haptic force signal, and we emphasize that sophisticated implementation of FEM for modelling a partial meniscectomy punch is outside the scope of this paper. The main contributions of our work can be summarized as • demonstrating an interpolation-based and a Kalman filter-based haptic control strategy for using intraoperative data to model a nonlinear partial meniscectomy punch; • implementing an inexpensive partial meniscectomy haptic simulation based on tool-tissue interaction data, and make all software and hardware openly available on GitHub; and • performing user studies to assess how experienced orthopedic surgeons consider interpolation-and Kalman filter-based methods to resemble a true partial meniscectomy punch, and compare these with an idealized linearelastic ideally-plastic haptic signal.
The remainder of this paper is organized as follows. First, we describe how intraoperative data were collected from a cadaveric specimen using an arthroscopic punch equipped with a force-and position sensor. Second, we demonstrate how the haptic control strategies were implemented on an open-source haptic system with dimensions similar to those of the arthroscopic punch. Finally, we present a user study conducted to assess the face validity of the interpolation-and Kalman filter methods.

A. DATA COLLECTION EQUIPMENT
In arthroscopic partial meniscectomy of a traumatic bucket handle tear, the torn part of the meniscus is removed by using an arthroscopic punch. In this study, the tool-tissue interaction forces experienced by the surgeon were investigated by performing an experiment on a cadaveric meniscus. For this experiment, an Acufex 1.5 mm upbiter punch (no. 7207735) was used. A prototype for recording punch interaction forces during partial meniscectomy was designed and implemented, and the setup is shown in Fig. 2. We emphasize that this setup is intended as a low-cost, and open source alternative for capturing the main qualitative haptic features of a partial meniscectomy punch in a controlled lab environment.
The physical properties recorded were position and force over time. A COM-09806 potentiometer was used for position sensing. For force sensing, an FSR07CE force-sensitive resistor (FSR) was selected owing to its low-profile design and low cost. The reported hysteresis for the model was 5 % for 100 actuations of 1 kg [28]. The maximum load capacity was 5 kg. An Arduino MKRZERO microcontroller was selected as the basis for software implementation.
To integrate the FSR, a thumb-sized clip-on assembly was designed such that the FSR could be placed between a base attached to the punch and a button. Thus, the contact area with the force sensitive resistor can be maintained constant. The base and button were 3D-printed using a Formlabs Form 3 3D printer. To integrate the potentiometer, two frames were designed to be attached to the handle and trigger of the punch. The frames were designed such that the potentiometer formed the connection, but did not transfer any forces other than those for rotating the potentiometer. The frames were also 3D-printed using a Formlabs Form 3 3D printer. Computer-aided design (CAD) files, bills of materials, wiring diagrams and software can be downloaded from GitHub [29]. The prototype is shown in Fig. 2. The potentiometer and FSR readings were sampled using the Arduino MKRZERO microcontroller board and imported to MATLAB where nonrelevant data were removed and the potentiometer data were filtered using a low-pass filter with a cutoff frequency of 5 Hz.
The FSR was calibrated using a load cell. FSR-readings were recorded from 0 N to 9.0 N in 1.0 N increments, and the calibration was repeated three times. Regression was performed in Excel to establish a relationship between the FSR-reading and the corresponding force. The following relation was established: where r FSR is the FSR-reading, and f * is the estimated FSR-force in Newtons. The root-mean-squared (RMS) error was calculated to assess accuracy, and was found to be 0.56 N.

B. DATA COLLECTION EXPERIMENT
A human cadaveric knee, acquired from Science Care USA [30] was dissected and prepared for partial meniscectomy. Eight series of ten cuts each, where four were performed slowly and four were performed fast, were conducted on the meniscus using the setup described previously. The knee before and after meniscectomy is shown in Fig. 3. Two series of reference cuts in a sheet of paper, where half were slow and half were fast, were also conducted using the same setup. The sampling frequency was 333 Hz. The data set is available on IEEE DataPort [31].    4 shows the FSR-force plotted against the puncher angle for a series of blank reference punches in a sheet of A4-paper, and a series of partial meniscectomy punches. For the reference punches, there was a smaller variation in the interception angle than that in the partial meniscectomy punches.
The power spectral density (PSD) of the unfiltered position signal for partial meniscectomy was obtained. The PSD showed that most of the frequency components were at the low-end of the spectrum, and typically below 5 Hz. Similar to what was found in [13], this indicates that fidelity must only be attained at low frequencies. It is likely that in this lowfrequency region, the material and tool-tissue nonlinearities dominate, and thus the emphasis should be made on recreating these.
By studying the anatomy of the force-position plot, as shown in Fig. 4, the following four points were defined: (i) Interception point: the position where the punch intercepts with the tissue, and the load starts increasing; (ii) yield point: the force and position where the punch cuts through the surface of the tissue, and the force curve decreases rapidly; (iii) bottom-out point: the position where the punch reaches the end of travel, and the force again starts to increase; and (iv) opening point: the position and force where the punch loses contact with the tissue, and the maximum force is reached. System overview for the partial meniscectomy punch simulation. For rendering force signals using only the interpolation method, the Kalman filter block was bypassed, and the interpolation force signal, f IP was sent directly to the motor controller.

III. HAPTIC ARTHROSCOPIC PUNCH A. SYSTEM SETUP
The motivation for developing a novel open-source haptic arthroscopic punch system is to provide an open-source platform for haptic rendering of arthroscopic procedures involving punches. Most generic commercial haptic devices only provide kinesthetic feedback in three or six DOFs. Some high-end systems, such as the Force Dimension Omega 7, also include pinch haptic feedback in a seventh DOF. However, these devices do not resemble the physical shape of the surgical instruments they attempt to represent. Although many commercial arthroscopic simulation systems have these features, most are protected by proprietary rights, and are thus less suitable for open research. To this end, we present our haptic arthroscopic punch system.
The haptic arthroscopic punch system is based on the impedance control strategy, where the force is computed as a function of the position. An overview of the system is shown in Fig. 5. The main components of the system are the haptic device, control software, and visualization software. The force computation strategies are presented in Section III-B.
A haptic device for rendering kinesthetic feedback for meniscus punching was designed as shown in Fig. 6. The hardware consisted of a BGM2804 brushless direct current (BLDC) gimbal motor, an AMT 102 rotary encoder with 2048 increments, a polulu 4 mm shaft clamp, a 4 mm precision shaft, a 4 mm ball bearing, and two custom designed handles with similar geometry to the Acufex upbiter punch. A high-resistance brushless gimbal motor was selected owing to its low weight, smooth motion, and low internal friction. The model has seven pole pairs. The electronics consisted of a 12V external power supply, an Arduino Due microcontroller (84 MHz), and a Simple FOC shield version 2.0.3 [32]. The Simple FOC shield was configured with BLDC driver pinouts 9 (Pwm A), 5 (Pwm B), 6 (Pwm C) and 8 (enable), and encoder pinouts 3 (A), 2 (B) and 4 (I), by soldering.
The control software was implemented on the Arduino Due microcontroller. The position signal θ from the haptic device was registered using an incremental encoder, and the handle position signal θ m was sent to the virtual environment, where collisions were detected and forces were computed as described later. A force signal f is then sent to the motor controller and converted to a control voltage u c which results in a motor torque τ . For matrix algebra operations, the Basic Linear Algebra library (BLA) was employed [33]. To control the BLDC motor, a field-oriented control (FOC) algorithm, provided by the SimpleFOC-library [32], was employed. The FOC algorithm was run with torque motion control and FOC modulation type space vector pulse-width modulation (PWM).
Visualization software was implemented on a computer using the Processing software environment [34]. The deformed mesh y mesh and tool position y tool are sent from the control software through a serial port. A generic anatomic knee model obtained from magnetic resonance imaging scans and segmented using the method described in [35], was imported. A CAD-model of the Acufex arthroscopic punch was developed and imported. The mechanisms were modeled using rigid-body transformations. The visualization view is shown in Fig. 7.
The setup, including CAD-models, bill of materials, control software, visualization software, intraoperative data and anatomical models are openly available on GitHub [29].

B. FORCE COMPUTATION STRATEGIES 1) INTERPOLATION METHOD
The interpolation method takes a sampled reference cut from intraoperative data, and renders this using interpolation, as described by Okamura [11]. Specifically, the interpolation force at time step k, f IP k , is computed as a function of the position, by interpolation of the recorded intraoperative VOLUME 10, 2022  force, f * and position, θ * , as described in (2).
where f IP k is the force signal transmitted to the haptic device at the k'th time step, θ k is the trigger position of the haptic device, f * k is the recorded intraoperative force with corresponding intraoperative tool position θ * k at a given time step.

2) FINITE ELEMENT METHOD
To represent the finite element force signal, a simple linearelastic ideally-plastic finite element model based on 1D spring elements was implemented. An illustration is shown in Fig. 8. The purpose of the finite element signal was not to model the meniscus punch as accurately as possible, but to generate a generic linear-elastic ideally-plastic reference force signal of the meniscus deformation that had similar stiffness to the intraoperative data. Other more sophisticated FEM-strategies provide a much more accurate modeling of this problem. The linear stiffness was calibrated using intraoperative data. The mesh, described by the node position vector x nodes and element midpoint vector x mid was defined by a start point a, an end point b, a finite number of elements N , and an element length, l e , as shown in (3) and (4).
x mid = a a + Assuming that all the elements have the same stiffness k e and length l e , we obtain the following global stiffness matrix where For the current simulation these parameters were set to a = 0.5, b = 2.0, and k e = 1. Boundary conditions were then applied such that displacements at the first and last nodes were zero. Collision detection is performed using the tool position P tool (x, y), and the tool penetration y tool is calculated. Next, an interpolation vector c is defined from x tool and element index i to the left of the tool position is registered. Then, the nonzero elements of c are calculated as follows: and Assuming that y mesh at position x tool is equal to tool penetration y tool , we obtain the following relation: The tool-tissue interaction force f FEM is calculated from the tool penetration y tool by solving (10) using lower-upper (LU) decomposition.
Finally, the mesh displacement vector y mesh is calculated.

3) KALMAN FILTER METHOD
The Kalman filter method is based on sensor fusion of interpolation and FEM force signals to estimate a new haptic force signal at time step k, f k , as shown in Fig. 9. A mathematical state-space representation of the arthroscopic punch force is formulated in the one-dimensional space domain, and the haptic force, f , and stiffness, ζ , are selected as states, as shown in (12) and (13).
where θ k is the change in handle position per time step. The state space representation leads to the system matrix shown in (14). As such, the a priori state estimate at time step k iŝ As the force signal is sampled by alternating between setting the measurement at time step k, z k , equal to f FEM k or f IP k , the measurement matrix is taken as shown in (16).
The discrete Kalman filter, as shown in [36], was then implemented using the update equations shown in (17) and (18) Here,x k−1 is the a posteriori estimate at time step k −1, P − k is the a priori estimate error covariance at time step k, P k−1 is the a posteriori estimate error covariance at time step k − 1, and Q is the system covariance matrix. The measurement update equations are shown in (19), (20) and (21), where K k is the Kalman gain,x k is the a posteriori state estimate, P k is the a posteriori estimate error covariance at time step k, I is the identity matrix, and the rest as previously defined.
For the Kalman filter method, the process noise covariance matrix, Q, and the measurement noise variances, R FEM and R IP , must be set in advance. In this implementation, they are also assumed to be constant. The process noise covariance matrix Q was determined from the recorded intraoperative stiffness, ζ * , and force data, f * , from a sampled punch shown in Fig. 4. The covariances, σ , were calculated using MATLAB, and resulted in the following covariance matrix: .0275 0.02772 0.2772 6.1024 (22) Furthermore, it was assumed that the measurement variance, R FEM , for the FEM signal was low. For the intraoperative data, the measurement variance, R IP was determined from the force sensor calibration data. The FEM-variance was tuned using system simulations in MATLAB so that the KF-force was about equally biased from the interpolation-and FEM signals.
The measurement variances were set as follows: C. HAPTIC SIGNALS DEMONSTRATION Fig. 10 shows the haptic force signals from the presented system calculated using the interpolation, FEM, and Kalman filter methods. For this demonstration, the FEM produces a linear-elastic force signal until a yield force of 7 N is reached, which is represented by an ideal plastic force. The interpolation method recreates one of the nonlinear forceposition curves shown in Fig. 4. After the bottom-out is reached, the signal is reduced to zero. The yield force was set to 7 N, and the interception point to four degrees for the FEM-signal, which is different than from that for the interpolation signal, to show that the Kalman filter method can account for variations in interception. This shows how the Kalman filter method can recreate nonlinear haptic features from intraoperative data, implemented using interpolation, while being synchronized to a linear FEM-signal. Fig. 11 shows how the Kalman filter haptic force signal can be tuned to be more biased towards the FEM-or interpolation force signal, by selecting different values for the measurement variance, here demonstrated by the interpolation variance, R IP . In addition, the Kalman filter was tuned to be responsive to measurement inputs to model the process nonlinearities, resulting in some minor noise in the estimate. In Fig. 11, the three Kalman filter curves were smoothed out using a two-point moving average filter. This demonstrates how the noise can be mitigated, while still sufficiently modelling the process nonlinearities.

D. COMPUTATIONAL PERFORMANCE ANALYSIS
To assess the computational efficiency of the interpolation and Kalman filter methods, the force computation strategies described in Section III-B were set up and tested individually. The methods were also checked for their sensitivity in refining the FE-mesh. The highest refresh rate of 496 Hz was achieved using only the finite element method. The speed was decreased to 326 Hz, going from five to ten elements, and further to 167 Hz going from ten to 15 elements. When using the interpolation method for haptic rendering, and FEM only for visualization, the maximum refresh rate was 433 Hz. The performance decreased to 294 Hz and 160 Hz when the mesh was refined to ten and 15 elements respectively. The maximum refresh rate achieved with the Kalman filter method was 375 Hz, which was 121 Hz slower than that of FEM only. The refresh rate dropped to 276 Hz going from five to ten elements, and to 147 Hz by going from ten to 15 elements.
The computational performance of the Kalman filter method was lower than the FEM and interpolation methods when the FE-mesh consisted of a small number of elements. However, as the FE-mesh was refined the difference in computational performance was diminished, as shown in Fig. 12. This shows that the computational cost of the FE-calculation was much higher than that of the interpolation-and Kalman filter calculations.

IV. USER STUDIES A. OVERVIEW
The purpose of the user study was to investigate how experienced surgeons conceived the interpolation-and Kalman filter-based force computation strategies with respect to face validity, compared with the idealized linear-elastic ideally-plastic FE haptic signal. Application of the KF-and interpolation methods are mostly relevant for patient-specific simulations where intraoperative tool-tissue interaction data are available, such as in a post-operative virtual surgery database for challenging cases, as described in section I and previously elaborated on in [9]. As such, experienced orthopedic surgeons was selected as the targeted user group.
Two user groups were included in the study. The first cohort consisted of orthopedic surgeons working with arthroscopy. Six participants volunteered (N = 6), of whom five were male and one was female. Their ages ranged from 47 to 76 years old. The average experience with arthroscopic surgery was 18 years, and the minimum reported experience was 10 years. The second cohort consisted of adults with no surgical experience. This cohort was included as a control group. Five participants (N = 5), three males and two females, volunteered. Their ages ranged from 23 to 31 years old. No financial compensation was provided to any of the participants.
The participants were asked to complete a physical reference benchmark punch, as well as two categories of virtual partial meniscectomy punch simulations, including (i) haptic feedback without visualization, and (ii) haptic feedback with visualization. The motivation for this study design was to first isolate the haptic rendering to collect data describing only the kinesthetic sensation, and then include the full simulation. The procedure is described in detail in the following section.

B. SETUP AND PROCEDURE
The experimental setup was as shown in Fig. 13. For the FE force signal, the interception point, linear stiffness, yield force and mesh were set as described in section III-B2. The interpolation force signal was based on one of the cuts shown in Fig. 4, and rendered as shown in Fig. 10. For the KF-method, the process variance was set as shown in (22), and the measurement variances as shown in (23) and (24). The methods were tested with five elements corresponding to 496 Hz, 425 Hz and 375 Hz for the FEM, interpolation and KF methods respectively.
First, the participants were given information on the purpose of the experiment and instructions for using the equipment. Informed consent was obtained. They were then asked to fill out personal and background information in a paper questionnaire. If consent was provided, the session was video recorded. The participants were encouraged to express their thoughts during the testing for qualitative assessment.
Participants were then given a kinesthetic reference benchmark. They were required to perform at least three cuts using an Acufex 1.5 mm upbiter punch (no. 7207735) in a piece of bacon. They were informed that the intention was to replicate a partial meniscectomy punch.
Next, the participants were asked to test the interpolation, FEM and KF haptic rendering strategies without visualization. The three methods were anonymized and labeled ''Alternative 1'', ''Alternative 2'' and ''Alternative 3'', corresponding to the interpolation, FEM, and Kalman filter methods, respectively. The alternatives were presented to the participants in this order. The participants were asked to complete at least three virtual cuts for each alternative. After testing the respective method, the participants were asked to rate on a 5-point Likert scale how well they agreed with a given set of statements corresponding to face validity of the method tested. The statements are presented in Table 1. The statements were chosen to be different for the three methods, with variations in biased and unbiased statements to avoid acquiescence bias. After all three alternatives were tested, the participants were asked to state which method resembled the benchmark cut the most and why.
After this, another round of assessment was performed, but now including the real-time 3D-visualization shown in Fig. 7. The alternatives were presented to the participants in reverse order, namely ''Alternative 3'', ''Alternative 2'' and ''Alternative 1''. As before, participants were required to complete at least three virtual cuts for each alternative, and asked to answer the survey questions after testing the respective method. Finally, the participants were again asked which alternative they preferred and why.

1) QUANTITATIVE
The Likert-scale survey results were grouped into professional and non-professional cohorts. Likert-scale results from survey questions with negative bias were inverted so that a high score indicated good performance, and a low score indicated poor performance. The performance was evaluated with respect to how well the respective method resembled the benchmark cut. Mean values (µ) and variances (σ 2 ) were calculated for the data set. The results are presented in Table 2. The interpolation method scored best with µ = 4.1 for the FIGURE 10. Interpolation, FEM and Kalman filter haptic force signals versus device handle position for a virtual partial meniscectomy punch. The graphs show how the interpolation and Kalman filter methods can recreate nonlinear haptic features from intraoperative data. The highlighted points shows how the Kalman filter method can use an idealized FEM-signal to account for variations in interception. The FEM interception point and yield force have been set differently than for the interpolation signal to illustrate this. Graphic renderings were created using Fusion360. professional group, and µ = 4.1 for the non-professional group. The second best method was the Kalman filter method with a score of µ = 3.2 in the professional group, followed FIGURE 12. Difference in computational performance of the Kalman filter, and interpolation method (including FEM for visualization only), relative to the idealized FE-simulation. As the FE-mesh is refined, the difference converges towards zero as the FE-calculations quickly become computationally expensive.
by the FEM with µ = 2.6. The ranking was the same for both cohorts.
To evaluate the differences in face validity between the interpolation, FEM, and Kalman filter methods, three t-tests assuming unequal variances were performed on the survey data using Excel. Data from the two rounds of testing (without and including visualization) was combined to provide a sufficient number of data points for the t-tests. Bonfer-VOLUME 10, 2022  roni correction was performed, setting α = 0.017. The results are presented in Table 3. For both professional and non-professional cohorts, the tests indicates that the users perceived a significant difference between the interpolation method and reference FEM-signal (p < 0.017). There were no significant difference between the other methods.

2) QUALITATIVE
In the survey, participants were asked directly which method they preferred and why. They were asked first after the first round (no visualization), then once more after the second round (including visualization). Each participant had therefore a total of two votes. The results are presented in Table 4.
In the professional cohort, there were a total of seven votes in favor of the interpolation method, four votes for the Kalman filter method, and no votes for FEM. One of the participants, a male surgeon with 20 years of arthroscopic experience, preferred the interpolation method, and reported that this was because ''it feels similar to a meniscus'' and ''it feels more realistic with a discontinuity than constant resistance''. Two other professional participants also reported similar statements.
In the non-professional cohort, there were five votes for the Kalman filter method, four votes for the interpolation TABLE 2. Survey results with means (µ) and variances (σ 2 ) for the interpolation, finite element (FEM) and Kalman Filter (KF) methods. The score is based on five-point Likert data, where maximum is 5 and minimum is 1.

TABLE 3.
Results from three T-tests with Bonferroni correction (α = 0.017). Tests indicate that participants in both the professional and non-professional groups conceived a significant difference between the interpolation-and FEM-methods. There was no significant difference between the other methods. method, and one vote for FEM. A male participant preferring the interpolation method reported that ''It felt more like the benchmark cut, with a snap at the end of the cut.'' Another participant also reported a similar statement. From the participants preferring the Kalman filter method, one participant reported that ''Discontinuity in resistive force was clear and concise, but some minor improvements might be needed.''

V. DISCUSSION
This study has presented an open-source haptic system for the simulation of a partial meniscectomy punch, and demonstrated implementation of an interpolation and a Kalman filter based method for using intraoperative data to provide nonlinear haptic feedback in arthroscopic surgical simulations. The face validity of the interpolation and Kalman filter methods was evaluated in a user validation study, and compared with haptic feedback from an idealized linear-elastic ideallyplastic FEM reference signal. The Kalman filter method demonstrated coupling of intraoperative data and FEM for haptic feedback in a partial meniscectomy punch simulation, but was outperformed by the interpolation method with respect to face validity.
Observing the results from the quantitative user study, it was evident that the interpolation method was significantly able to capture the nonlinear features in the partial meniscectomy punch, when compared with the idealized reference FEM-signal. There was no significant difference between the Kalman filter-based method and the idealized reference FEM-signal, which means that for the selected parameters, the Kalman filter-based method was not as effective in rendering the nonlinear features as the interpolation method. These findings were supported by participant statements, as well as the results shown in Table 4. In the professional cohort, users explained that they preferred the interpolation method because it most accurately modeled the expected TABLE 4. Results show how many participants preferred each of the respective methods when directly asked about which method resembled the benchmark cut the most. Participants were first asked after the first round (no visualization), and then once more after the second round (including visualization). One participant in the professional user group did not report the preferred method after round two, which is why there are only five responses in the professional user group in the second round. In the professional group, the interpolation method was preferred, followed by the KF-method. In the non-professional group, the KF-method obtained one more vote than the interpolation method.
''pop-through effect'' of the partial meniscectomy punch. Looking at the corresponding signals in Fig. 10, it is clear that this effect was most prominent for the interpolation method, and that the Kalman filter method, with the given set of parameters, did not sufficiently capture this nonlinearity. However, the Kalman filter can easily be tuned to be more based on intraoperative data than on the FEM-signal by selecting other values for the measurement variances, R IP and R FEM , as shown in Fig. 11. In the user study, the measurement variances were tuned so that the KF-force was about equally biased from the interpolation-and FEM signals. It is possible that tuning the Kalman filter to be more biased towards the interpolation signal could have improved the results in the user study. However, because the Kalman filter force signal is dependent on the interpolation force signal, it is not really possible for the Kalman filter method to capture the nonlinear features any better than the interpolation signal, and neither is this the purpose. The advantage of the Kalman filter method is that it also accounts for variations in points of interaction by using the FEM-force when no interpolation force is available, to allow for a meaningful simulation where haptic feedback is synchronized with the visual display, as shown in Fig. 10. This is especially important for post-operative digital twin applications, where haptic feedback from interpolation is more challenging to synchronize with FEM owing to patientspecific variations in anatomy and tissue deformations.
As can be observed from Fig. 10, some minor noise resulted in the Kalman filter force signal because the filter was tuned to be quite responsive to the interpolation and FEM signal inputs. This is necessary to be able to model the nonlinearities of the partial meniscectomy punch process. However, as shown in Fig. 11, the noise can be easily reduced by a filtering process, such as a two-point moving average filter.
Considering the user study design, the motivation was to first isolate the haptic rendering to collect user data describing face validity of the kinesthetic sensation only, and then include visualization to evaluate the face validity of the full simulation experience, and see if the methods performed differently under this condition. There was not enough data to show statistically that including visualization affected face validity, and this should be investigated in future studies. However, as can be observed from Table 4, the distribution for the preferred method was very similar in both cohorts going from 'no visualization' to 'including visualization'. The same trend was seen from the Likert-scale survey data. This indicates that including visualization have not significantly affected the performance of the haptic methods in this study.
Observing the results of the performance analysis described in Section III-D, it is important to note that both FEM and interpolation are necessary for the Kalman filter method. Similarly, the FEM is necessary for the interpolation method to provide a meaningful simulation with both visual and haptic feedback, and has therefore been included in the comparison to only compute the deformed mesh. For a low number of elements in the FE-mesh, the Kalman filter-based, and interpolation methods are slower than FEM only, with a performance of 375 Hz, 425 Hz and 496 Hz respectively. However, the advantage of these methods become clearer as the FE-mesh is refined and the associated matrices increase in size. As shown in Fig. 12, the difference in computational speed between the three methods converges towards zero as the FE-mesh is refined. This shows that FE-calculations are much more computationally expensive, and the KF-and interpolation methods are not affected by the number of elements in the FE-mesh. For a moderately sized 3D tetrahedral mesh, these computations should not make the simulation significantly slower. For reference, Liu et al. [37] recently modeled a human meniscus using 5799 tetrahedral elements for the medial meniscus, and 6646 tetrahedral elements for the lateral meniscus.
The intraoperative data collection and haptic arthroscopic punch systems were designed using low-cost components and open-source software, and are openly available on GitHub [29]. These systems, supplemented by the open data set available at IEEE DataPort [31], provide the sufficient hardware and software for intraoperative data-based haptic rendering of a partial meniscectomy punch.
The haptic arthroscopic punch provide an inexpensive alternative to costly kinestaethic haptic devices with pinch feedback. It should be noted that the desired haptic refresh rate of 1000 Hz was not attained, because the system was implemented on an Arduino Due micro controller (84 MHz) with a single thread. However, the achieved simulation speed was sufficient for the purpose of rendering a partial meniscectomy punch because of the low stiffness and relatively small forces involved. Colgate and Henkel [38] showed that for an impedance-controlled haptic system, passivity is influenced by sampling time T , virtual stiffness K , virtual damping B, and physical damping (or friction b, as shown in (25).
This shows that the maximum achievable virtual stiffness is proportional to the sampling rate. It is important to note that b includes not only the physical damping provided by the haptic device, but also by the human operator. When analyzing the force signals shown in Fig. 4, the maximum stiffness for the partial meniscectomy punch was found to be K = 4.39 N/deg. Further, the maximum handle velocity from the haptic signal was found to beθ = 57 deg/s. Analyzing the the system at the given stiffness and handle velocity, and no virtual damping, the human operator must apply an external damping force of F = KT 2θ = 0.33 N to ensure passivity at 375 Hz refresh rate. This is small, and well within the grasping capabilities of most adults. Compared to the target refresh rate of 1000 Hz, this external force would be 0.13 N, which is smaller, but only by 0.2 N. Conversely, by decreasing the refresh rate to 25 Hz, this external force would be 5.0 N, which is more than 60 % of the virtual force rendered, and would be very noticeable for a human operator. No instability issues were reported during the user studies.
In future work, the interpolation and Kalman filter-based force computation strategies should be explored together with a more sophisticated implementation of FEM, such as the corotational formulation [39], total Lagrangian explicit dynamics (TLED) [40] or reduced-order extended Kalman filter FEM [41], [42], [43]. Novel position tracking systems, such as that recently demonstrated for knee arthroscopy by Ma et al. [44], combined with a 6 DOF force-torque sensor, could provide the basis for in-vivo intraoperative data collection for other arthroscopic procedures. An interesting comparison would be to study the simulation speed and face validity between a real-time finite element strategy optimized for speed, such as TLED, in combination with the use of intraoperative data, and a finite element strategy optimized for modelling topology change, such as the corotational cut FEM [26].

VI. CONCLUDING REMARKS
This paper has presented an open-source haptic system for the simulation of a partial meniscectomy punch, and explored an interpolation-based and a Kalman filter-based method for intraoperative interaction data-based patient-specific haptic feedback in arthroscopic partial meniscectomy punch simulation. It was demonstrated that the Kalman filter-based method could couple finite element and intraoperative interaction data signals, and that the resulting haptic signal can be tuned to be biased towards either of the input signals. The system was evaluated with respect to face validity in a user study that included six orthopedic surgeons and five non-professionals. The user study showed that experienced surgeons could distinguish between a linear-elastic ideally-plastic finite element-based, and a nonlinear interpolation based haptic feedback signal in a partial meniscectomy punch simulation (p < 0.017), and that surgeons preferred the interpolation method as a haptic rendering strategy in this situation. These findings indicate that the Kalman filter method was not as effective as interpolation in rendering the nonlinear haptic feedback given the selected parameters.

ETHICAL CONSIDERATIONS
The use of cadaveric specimens have been considered by the Norwegian Regional Ethics Committee (REK Midt), application 214114. User studies were considered according to the General Data Protection Regulation (GDPR) by the Norwegian Centre for Research Data (NSD), application no. 271881.
MORTEN D. PEDERSEN received the M.Sc. degree in engineering cybernetics and the Ph.D. degree in engineering cybernetics from the Norwegian University of Science and Technology (NTNU), Trondheim, Norway, in 2009 and 2017, respectively.
Since 2018, he has been working as an Associate Professor at the Department of Engineering Cybernetics, NTNU, where he teaches several courses in control systems technology. His research interests include modeling and control of fluid-mechanical systems, pedagogical topics related to the teaching of control theory, biomedical applications of control theory, and the fundamentals of cybernetics.