# ECE496Y Final Report: Edge Machine Learning for Detecting Freezing of Gait
in Parkinson’s Patients

Purnoor Ghuman, Tyama Lyall, Usama Mahboob, Alia Aamir
Supervisor: Professor Xilin Liu
The Edward S. Rogers Sr. Department of Electrical and Computer Engineering
University of Toronto

###### Abstract

Parkinson’s disease is a common neurological disease, entailing a multitude of
motor deficiency symptoms. In this project, we developed a device with an
uploaded edge machine learning algorithm that can detect the onset of freezing
of gait symptoms in a Parkinson’s patient. The algorithm achieved an accuracy
of 83.7% in a validation using data from ten patients. The model was deployed
in a microcontroller Arduino Nano 33 BLE Sense Board model and validated in
real-time operation with data streamed to the microcontroller from a computer.

###### Index Terms:

Parkinson’s disease, microcontroller, Freezing of gait, deep learning, CNN

##  I Introduction

Parkinson’s disease (PD) is a neurodegenerative disorder, the second most
common following after Alzheimer’s disease. PD patients experience a multitude
of movement disorders, in particular, one of the common movement disorders is
the Freezing of Gait (FoG), prevalent in the late stages of PD [1]. A
patient’s walking and turning is impeded, thereby increasing instability and
the risk of falls, and as a result potential injuries.

Previous treatment for FoG has included a variety of options; such as medical
treatments and invasive procedures, each with varying results and potential
adverse effects [2]. Studies have also been done on cognitive exercises such
as Rhythmic Auditory Stimulation (RAS). RAS introduces an external auditory
stimulus for patients’ to synchronize their steps, enabling them to initiate
walking with a higher rate of success. This solution requires PD patients to
be aided by medical professionals, but there is yet to be a feasible solution
that detects the onset of FoG to make this form of treatment feasible [3].
Innovative research has gone into wearable devices with acceleration sensors
that detect the onset of FoG via the frequency components of limbs, but had
inconsistent results due to inability to detect FoG with higher accuracy [4].

The above studies give way to our motivation, which is the utilization of
machine learning. Machine learning algorithms can be trained on a set of
patient data and with some optimization, learn the onset of FoG episodes. This
allows the timely delivery of rhythmic auditory and electrical stimulation to
the patient to enable them to overcome difficulty with walking initiation and
turning [4]. The motivation for this project is to implement a non-invasive,
feasible method for detecting the onset of FoG in PD patients. The gap is a
device which can run a machine learning algorithm locally on an embedded
system. This device should accurately predict, based on input data of the
patient’s movements, whether a patient is experiencing FoG or not and output a
signal indicating this to an external source. An algorithm like this which
processes data locally is called an edge machine learning algorithm, resulting
in more convenience for PD patients wearing it to successfully detect the
onset of FoG symptoms [5, 6, 7, 8].

The goal of this project is to develop a microcontroller device with an
uploaded edge machine learning algorithm that can detect the onset of freezing
of gait symptoms in a PD patient from their acceleration data in a timely
manner with limited power resources [9, 7].

Figure 1: Introduction of proposed FoG detection using edge machine learning.

##  II Methods

###  II-A System Level Overview

The final design being implemented by the team consists of two main aspects.
The first part is the training and development of the machine learning model.
For this part, the initial step is the cleaning of the PD patient acceleration
data, shown in Figure 2. Unwanted timestamps, where the patient was not
walking, are removed in this step. Also, only thigh acceleration is considered
in this project because an adequate idea of a patient’s gait can be determined
from their thigh acceleration data.

Figure 2: Illustration of the data processing steps. Figure 3: Simplified
system context diagram.

The cleaned data, consisting of multiple timestamps with their corresponding
acceleration values and classifications, is then segmented into windows of 2
seconds, with each window having its own classification as FoG, denoted as a
1, or No-FoG, denoted as a 2. Next, our team trained a Convolutional Neural
Network (CNN) with the time windows as inputs to classify each 2 second window
as FoG or No-FoG.

In the second part of the project, the trained CNN model is then extracted as
a TensorFlow Lite model and uploaded onto an Arduino BLE Sense 33
microcontroller [7]. The input raw test data is stored on a computer in a CSV
file and written into the board using a script which utilizes the “pandas”
library. A connection between the script and the Arduino port is made to write
data into the board, process a window using the loaded CNN model, and extract
the result, all done in real time.

The system context diagram shown in Figure 3 summarizes the above information
succinctly, showing each working module of the project.

###  II-B Module Level Design

####  II-B1 Data Pre-processing

We took the data from the .csv files and separated them in terms of the type
of acceleration data e.g., thigh vs. ankle acceleration. Used the matplot
library to plot graphs of acceleration values against a time axis. Analyzed
the graphs to determine where and how often the data is useful as there are
segments in time where the patient is not moving so some time intervals are
not useful for our classification task.

We converted patient data files to suitable format (.csv) for later
visualization and manipulation. Removed data not a part of the experiment
(data with annotation ‘0’), determined by reading through provided
documentation. In addition, from data visualization, removed periods of no-
walking for each of the csv files using Microsoft Excel.

The results of data cleanup was data ready for further manipulation and usage
by model. Data visualization enabled further data cleanup through graphing the
data points and analyzing them visually.

####  II-B2 Data Segmentation and Augmentation

We developed a subroutine to divide the continuous acceleration data into two
second windows and return a list of data points, where each element represents
one 2 second window, with the corresponding accelerations across it and a
corresponding list of labels for each window. From the patient data, FoG is
denoted with a 2 and No-Fog is denoted with a 1, but the labels are passed
through a function to change the labels to 0 for No-FoG and 1 for FoG to
correspond with labels typically used in binary classification tasks. To
augment the data and yield more data points, a sliding window approach was
used, where each window has an overlap with the previous of one second.

In order to generate more instances of the minority class, instances of FoG in
our case, oversampling was performed. This involves randomly selecting a
predetermined amount of data points from the minority class to duplicate in
order to improve the model’s ability to learn on this class [10].

The final results of this module is an augmented dataset, ready to be passed
into a model. The dimension of the input data is (8711, 129, 3), where 8711 is
the number of two second windows created, 129 is the number of timestamps in
each window, and 3 is the number of accelerations provided for each timestamp,
one for each dimension x, y, z.

Figure 4: A simplified model diagram.

####  II-B3 Model Implementation

We model design entailed the initial research and outlining of initial model
implementation. Decisions were made here to start building CNN with a simple
2-layer MaxPooling and Convolutional layers model. We also made the decision
to use an optimizer to adjust weights and the learning rate to reduce losses
during training [11]. Concluded using Adaptive Moment Estimation (Adam) as an
optimizer during model training, due to its fast convergence and efficient
process [11, 12].

We researched and utilized Keras library for the implementation of a CNN model
[6, 14]. From the initial model design, obtained additional information
regarding how to structure the model. Initial implementation used in training
included the usage of 2 Convolution layers, a MaxPooling layer, and 2 Dense
layers [13]. Dense layers were utilized, in which we defined an activation
function, to change dimensions of expanded units into a classification output
[14, 15]. The model was then compiled using the previously mentioned ADAM
optimizer and our choice of loss function [11]. The loss function we used was
binary cross entropy loss, utilized for the purpose of determining how well
our predictions are [16, 8].

The model was trained and validated using the dataset splits previously
mentioned, along with the usage of the Keras library fit method. This method
enables model training and validation, with the corresponding ability to
adjust several parameters [17].

The final results of this module were our initial model, which can be found in
Appendix F. We also conducted and received results of its initial training and
validation, which can be found in Appendix G. From observing the validation,
as well as the accuracy corresponding to No-FoG and FoG classes, it was
determined further iteration was required, which is discussed in the
subsequent module.

####  II-B4 Model Optimization

Selected model hyperparameters to tune (number of filters in convolutional
layer, learning rate, epochs, batch size) and used Keras.tuner library and
scikitlearn.GridSearchCV to learn optimal values for each one. See Appendix H
for values over which tuning was conducted on each hyperparameter and
corresponding code [17].

Regularization techniques are used to prevent the model from overfitting to
the training data. Techniques used include the introduction of a dropout layer
in the model structure, which randomly drops some inputs to improve the
model’s generalization [18] and the use of early stopping during training,
which will stop training if the validation loss stops decreasing over many
epochs. In order to improve the model’s accuracy on the minority class, class
weights were calculated based on the ratio of positive to negative examples in
the training data and those weights are used to bias the classifications made
during training.

The final result was a model with hyperparameters tuned using the validation
set to maximize accuracy on both classes and to generalize well to brand new
data.

###  II-C Deployment to Microcontroller

We worked on converting tensorflow model to a tensorflow lite model and
optimized it for memory footprint in order to ensure a small size so it can
fit the memory constraint [18]. In order to optimize we decided to analyze the
final model and merge any common nodes and remove any unused nodes from the
graph. We also decided to freeze the graph by converting the weights into
constants so it would take up less space than having them as variables. We
also decided to quantize the inputs and outputs by decreasing the precision of
the variables so that they would take up less space. This only resulted in a
minimal decrease in accuracy so we decided to keep it.

We set up the Arduino board to upload the model and run continuously while
taking in inputs and calculating output of the model and displaying it on the
serial monitor. We had to set up error logging to ensure any errors caught
during inference are caught and had to account for how much dynamic and static
memory is available on the board to determine what size of the input can be
sent to the board at a time to perform inference and that there was enough
space for the board to store the resulting output classification.

We designed an efficient method to write raw acceleration data into the
microcontroller, which would be then processed by the embedded machine
learning model (Module 2.3.5), and read the produced output from the board.
This is implemented using a Python script, which utilizes the pandas library,
to write data “chunks”, or windows, into the Arduino board [19]. After the
script has written enough data to cover a full window of patient
accelerations, the window is passed into the model for processing, which is
then read by the script, as shown in Appendix L.

Since an Arduino board cannot take multiple rows of input at once, all data
needs to be passed into the board individually, byte-by-byte. This presents a
challenge, as on the Arduino board side, we need to know when to stop
collecting data from the Python script and process. This is done by defining
two states called ‘Read’ and ‘Process’. In the ‘Read’ state, the board reads
the data provided by a computer and populates the appropriate variables, until
the right number of variables have been populated. At that point, the state
switches to ‘Process’ and data is passed into the machine learning model,
after which an output is produced and the state is switched back to ‘Read’.

The implementation of this module results in our raw test data being stored on
a computer and streamed directly to the Arduino board, instead of taking up
space on the board itself.

##  III Results

###  III-A Assessment of Final Design

Overall, assessment of the final design can be broken down into two primary
parts: the binary classification model and the hardware implementation. While
the two aspects are interconnected and work together, their performance can be
assessed separately.

The model has been tuned using a validation set and trained on the training
set, and performs reasonably well on new data, with an overall accuracy of
83.7%. The accuracy on no FoG examples was 84% while the accuracy on FoG
examples was 80%. The accuracy is not as high as the desired threshold
originally set by the team, largely due to the imbalance in the data provided.
Development and training of the model was hindered by the fact that the
dataset contained much more data for No-FoG than data featuring instances of
FoG, making it challenging to develop a model that could recognize FoG
instances with high accuracy. The original plan to use an SVM for
classification was discarded, as it was impossible to train an SVM that was
capable of differentiating between the two classes; it would classify every
point either as one class or the other. The team moved on to researching deep
learning approaches and eventually decided on a CNN due to their weight
sharing capabilities, which lower the number of learned parameters and
therefore the size of the model. Due to space constraints, the final model has
three convolutional layers followed by one fully connected layer with two
outputs, one for each class. Through data standardization, weight biasing, and
various regularization techniques outlined in section 2.3.4, a final model
design was decided on. The final model size was 0.4 MB, which is well below
the memory requirement established by the team.

The hardware implementation involves downloading the trained CNN model to a
header file which can be uploaded onto the Arduino microcontroller. From
there, input data is written to the board using a Python script. Another
program is also written to the board, which can take in the data being written
and segment it continuously into two second windows. Each window is then
passed into the trained model for classification and the output of the model
will be written back to the computer to indicate when FoG is detected. Further
details on how each requirement is tested and verified by the above design is
discussed in the next section.

####  III-A1 Verification of Correct Input into Arduino Board

To confirm that the Arduino board is correctly receiving the acceleration data
being provided from a CSV File, a simple check is performed in the Python
script providing the data. In this check, one acceleration value is provided
as input into the board. The board then takes that acceleration value and
increments it by one and prints it on the monitor. The Python script, after
waiting for something to be printed on the Arduino terminal, consumes the
printed value and prints it on its own terminal. The input value and the
output received from the Arduino board is shown in Figure M3 in Appendix M.
The Python code performing this check can be seen in Figure M1 in Appendix M
and the Arduino code is shown in Figure M2 in Appendix M. This check confirms
that acceleration data being read into the board is received correctly and
that the processing done locally is also correct, as the output value must
always be an increment by one of the input value.

####  III-A2 Verification of Correct Input into Arduino Board

To confirm that the machine learning model is, in fact, uploaded on the board,
we must place the “model.h” file containing the model in the same folder as
the Arduino file uploading it on the board. After the upload happens, it must
be ensured that the upload happened successfully and that the correct number
of bytes are written on to the board ( 400,000 in this case). This can be
confirmed by looking at the output of the upload shown in Figure 3, where a
successful upload is shown and 384952 bytes of memory of the board has been
filled.

####  III-A3 Verification of ML Model Accuracy

To confirm that the model’s accuracy meets our objective, we generated the
trained model’s predictions on the test set, which it had never seen before.
Then, compare the predictions with the true labels and generate a percentage
from that which represents the accuracy.The accuracy on no FoG examples was
84% and the accuracy on FoG examples was 80%, yielding an overall accuracy of
83.7%. To further quantify the model’s performance, we can calculate its
accuracy on each class, and generate a confusion matrix, which will give us
the number of false positives and negatives as seen in Fig. 5.

Figure 5: Confusion matrix metrics and accuracy on each class for final model.

####  III-A4 Verification of Model Memory Requirement

To confirm that the model size will not exceed the 1MB flash memory capability
of the microcontroller, we can write code to generate a file model.h
containing our trained model, and then output the size of that file in bytes.
As seen in Figure 5, the model size is 298,600 bytes, or around 0.3 MB,
meaning it meets the memory requirement.

####  III-A5 Verification of Classification Output

For classification, a 2 second window (input) is sent to the board which uses
the uploaded model to perform inference based on this input and outputs a
classification result of either FoG or No-FoG. This result is then displayed
on the serial monitor of the computer that is connected to the board as the
board itself has no screen to display the result. As it has not been fully
implemented on the board yet the inference and output will be shown as follows
in Fig. 6.

Figure 6: Code to print classification onto serial monitor.

####  III-A6 Verification of Memory Resources Utilized on the Board

Given that the board has a total of 1 megabyte of flash memory, we had to make
sure that the size of our classification sketch was below that. This led to a
limitation on the number of layers in the model, otherwise the size would
increase exponentially. We had to find a good balance between having a complex
model with multiple layers and weights to ensure good accuracy but also making
sure that its size does not exceed the 1 megabyte constraint.

The model uses 384952 bytes of flash and only 21% of the dynamic memory and
40% of storage memory. The rest will be used to store input (2 second windows)
to be sent to the model for inference. Therefore the memory constraint set by
the board’s onboard memory is satisfied.

##  IV Conclusion

After extensive research and multiple iterations, the team was successful in
producing a CNN machine learning model that achieves the accuracy and
specificity set out in the requirements. We were able to meet all the
requirements laid out in our project and therefore accomplish our goal of
correctly predicting the onset of FoG based on the acceleration data provided
in real time that is run locally on the microcontroller.

As can be seen in the verification and testing section, the model is small
enough to be uploaded on the board’s limited memory and is able to run
inference within 5 seconds. Therefore, it is suitable to help patients
identify FoG onset while giving them enough time to be able to react such as
sitting down before it becomes active. The accuracy of the model is high
enough to be utilized by patients without giving inaccurate results.

Throughout the project, we had to make multiple changes in order to achieve an
accurate model. For example, initially we had planned to extract features from
the data manually, but that yielded very low accuracy [19]. We had to iterate
over multiple variations of the model, from using an SVM to a Neural Network
(NN) to finally using an CNN. Through this project we realized that in order
for the model to be compatible with the microcontroller we were restricted to
a few libraries such as Tensorflow. From this project, we learned that there
are many different algorithms that can be used to classify FoG onset, but only
a few finely tuned models were able to achieve all the requirements such as a
high accuracy and a small enough memory footprint so that it can be uploaded
to the microcontroller.

In the future, we can potentially use the acceleration data taken directly
from the microcontroller’s sensors as a patient walks to run inference on data
as it is collected to be used in a real-case scenario [20].

## References

  * [1] T. Bikias, D. Iakovakis, S. Hadjidimitriou, V. Charisis, and L. J. Hadjileontiadis, “Deepfog: an imu-based detection of freezing of gait episodes in parkinson’s disease patients via deep learning,” _Frontiers in Robotics and AI_ , p. 117, 2021.
  * [2] B. Su, C. Smith, and E. Gutierrez Farewik, “Gait phase recognition using deep convolutional neural network with inertial measurement units,” _Biosensors_ , vol. 10, no. 9, p. 109, 2020.
  * [3] X. Liu, M. Zhang, A. G. Richardson, T. H. Lucas, and J. Van der Spiegel, “Design of a closed-loop, bidirectional brain machine interface system with energy efficient neural feature extraction and pid control,” _IEEE transactions on biomedical circuits and systems_ , vol. 11, no. 4, pp. 729–742, 2016.
  * [4] M. Bachlin, M. Plotnik, D. Roggen, I. Maidan, J. M. Hausdorff, N. Giladi, and G. Troster, “Wearable assistant for parkinson’s disease patients with the freezing of gait symptom,” _IEEE Transactions on Information Technology in Biomedicine_ , vol. 14, no. 2, pp. 436–446, 2009.
  * [5] M. Merenda, C. Porcaro, and D. Iero, “Edge machine learning for ai-enabled iot devices: A review,” _Sensors_ , vol. 20, no. 9, p. 2533, 2020.
  * [6] L. Borzì, I. Mazzetta, A. Zampogna, A. Suppa, G. Olmo, and F. Irrera, “Prediction of freezing of gait in parkinson’s disease using wearables and machine learning,” _Sensors_ , vol. 21, no. 2, p. 614, 2021.
  * [7] X. Liu and A. G. Richardson, “Edge deep learning for neural implants: a case study of seizure detection and prediction,” _Journal of Neural Engineering_ , vol. 18, no. 4, p. 046034, 2021.
  * [8] ——, “A system-on-chip for closed-loop optogenetic sleep modulation,” in _2021 43rd Annual International Conference of the IEEE Engineering in Medicine & Biology Society (EMBC)_. IEEE, 2021, pp. 5678–5681.
  * [9] X. Liu, M. Zhang, and J. Van der Spiegel, “A low-power multifunctional cmos sensor node for an electronic façade,” _IEEE Transactions on Circuits and Systems I: Regular Papers_ , vol. 61, no. 9, pp. 2550–2559, 2014\.
  * [10] J. Brownlee, “Random oversampling and undersampling for imbalanced classification,” _Machine Learning Mastery_ , 2020.
  * [11] S. Doshi, “Various optimization algorithms for training neural network,” _Towards data science_ , vol. 13, 2019.
  * [12] X. Liu, H. Zhu, T. Qiu, S. Y. Sritharan, D. Ge, S. Yang, M. Zhang, A. G. Richardson, T. H. Lucas, N. Engheta _et al._ , “A fully integrated sensor-brain–machine interface system for restoring somatosensation,” _IEEE Sensors Journal_ , vol. 21, no. 4, pp. 4764–4775, 2020.
  * [13] E. Allibhai, “Building a convolutional neural network (cnn) in keras,” _Internet draft. Towards Data Science_ , 2018.
  * [14] P. Sharma, “Keras dense layer explained for beginners,” _Retrieved from machinelearningknowledge. ai/website: https://machinelearningknowledge. ai/keras-dense-layer-explained-for-beginners_ , 2020.
  * [15] X. Liu, A. G. Richardson, and J. Van der Spiegel, “An energy-efficient compressed sensing-based encryption scheme for wireless neural recording,” _IEEE Journal on Emerging and Selected Topics in Circuits and Systems_ , vol. 11, no. 2, pp. 405–414, 2021.
  * [16] S. Sharma, S. Sharma, and A. Athaiya, “Activation functions in neural networks,” _towards data science_ , vol. 6, no. 12, pp. 310–316, 2017.
  * [17] S. Lau, “A walkthrough of convolutional neural network—hyperparameter tuning,” _Medium. Retrieved_ , vol. 23, 2019.
  * [18] X. Liu, H. Zhu, M. Zhang, X. Wu, A. G. Richardson, S. Y. Sritharan, D. Ge, Y. Shu, T. H. Lucas, and J. Van der Spiegel, “A fully integrated wireless sensor-brain interface system to restore finger sensation,” in _2017 IEEE International Symposium on Circuits and Systems (ISCAS)_. IEEE, 2017, pp. 1–4.
  * [19] X. Liu, M. Zhang, A. G. Richardson, T. H. Lucas, and J. Van der Spiegel, “A 12-channel bidirectional neural interface chip with integrated channel-level feature extraction and pid controller for closed-loop operation,” in _2015 IEEE Biomedical Circuits and Systems Conference (BioCAS)_. IEEE, 2015, pp. 1–4.
  * [20] X. Liu, M. Zhang, and J. Van der Spiegel, “A low power multi-mode cmos image sensor with integrated on-chip motion detection,” in _2013 IEEE international symposium on circuits and systems (ISCAS)_. IEEE, 2013, pp. 2416–2419.


11institutetext: Graduate School of Informatics, Kyoto University, Kyoto,
Japan. 11email<EMAIL_ADDRESS>
22institutetext: Graduate School of Madicine, Kyoto University, Kyoto, Japan.

# Shape Reconstruction from Thoracoscopic Images using Self-supervised Virtual
Learning

Tomoki Oya 11 Megumi Nakao 22 Tetsuya Matsuda 11

###### Abstract

Intraoperative shape reconstruction of organs from endoscopic camera images is
a complex yet indispensable technique for image-guided surgery. To address the
uncertainty in reconstructing entire shapes from single-viewpoint occluded
images, we propose a framework for generative virtual learning of shape
reconstruction using image translation with common latent variables between
simulated and real images. As it is difficult to prepare sufficient amount of
data to learn the relationship between endoscopic images and organ shapes,
self-supervised virtual learning is performed using simulated images generated
from statistical shape models. However, small differences between virtual and
real images can degrade the estimation performance even if the simulated
images are regarded as equivalent by humans. To address this issue, a
Variational Autoencoder is used to convert real and simulated images into
identical synthetic images. In this study, we targeted the shape
reconstruction of collapsed lungs from thoracoscopic images and confirmed that
virtual learning could improve the similarity between real and simulated
images. Furthermore, shape reconstruction error could be improved by 16.9$\%$.

###### Keywords:

3D shape reconstruction Virtual learning Variational Autoencoder Endoscopic
image.

##  1 Introduction

Endoscopic camera images have been used to reach organs during surgery.
However, it is difficult to grasp the three-dimensional structure of organs
inside the body owing to the high uncertainty caused by organ deformation and
narrow visibility during surgery. Therefore, research is underway to
reconstruct the 3D shape of organs from endoscopic images that are acquired
during surgery. In this paper, we focus on the reconstruction of 3D shapes
using single-viewpoint endoscopic camera images, which is a key technique for
image-guided surgery.

Shape reconstruction of biological organs using endoscopic camera images
reported so far are: A study on the recovery of three-dimensional information
of the entire stomach including texture from endoscopic images of the
stomach[1], and model-based positioning of the liver[2]. The former
reconstructs shape information from video, but only recovers the shape of the
organ surface within the visible range, and the latter calculates the
deformation using parameter optimization, but the computational cost and
stability of the optimization calculation remains an issue. To address these
issues, there have been reports on attempts to predict 3D shapes and
deformations using statistical shape models and graph convolution networks,
which describe the morphological characteristics of biological organs in terms
of low-dimensional parameters[3][4][5][6]. In the field of medical imaging,
data augmentation with simulated data generated from statistical shape models
has been shown to be effective. Tang et al. achieved an improvement in
segmentation accuracy with a data augmentation method using statistical shape
models compared with a method using real images[7]. However, it has been
reported that if excessive simulated images are used for training, the
accuracy decreases in contrast to the increase in training data; hence, the
problem of treating simulated and real images equally still remains.

In this paper, we propose a framework for self-supervised virtual learning for
the shape reconstruction of deformed organs using image translation with
common latent variables between simulated and real images from endoscopic
camera images. This method solves the problem of insufficient training data in
medical imaging by using statistically generated simulated images. In this
paper, we refer to this learning method as “virtual learning.” We use a
Variational Autoencoder (VAE)[8] to extract latent variables that are common
to both simulated and real images, and construct a framework that accepts two
images with different features that humans may regard as equivalent as inputs
and converts them into identical synthetic images. Although some studies[9]
have reported the generation of medical images that represent pathological,
morphological, and anatomical variations of VAEs, there have been no reports
of their application to simulation-based learning, which is the subject of
this study. Since our framework can improve the similarity between simulated
and real images generated from 3D CT images of individual patients, it can
reconstruct the shape of lungs from endoscopic images in thoracoscopic lung
cancer resection with higher accuracy than conventional methods.

##  2 Method

Figure 1: Examples of real and simulated images in endoscopic images, (a) real
image $I_{R}$ extracted from a surgical video, (b) different appearances
generated by changing weights, and (c) simulated image $I_{S}$ corresponding
to the real image

###  2.1 Dataset and Problem Definition

In this study, 128 still images were collected from videos of of thoracoscopic
resections of 9 cases of the right lung and 3 cases of the left lung. The lung
shape $M$ during insufflation was constructed from the CT data taken before
the surgery, and the lung shape $\hat{M}$ during collapse was constructed from
the cone-beam CT data taken during the surgery. The lung shape is a mesh model
consisting of 402 vertices and 796 triangles.

Fig. 1 shows an example of an endoscopic image $I_{R}$ and a simulated image
$I_{S}$ generated from the statistical shape models. In Fig. 1(a), part of the
surface of the upper lobe of the lung and the inner wall of the thoracic
cavity can be seen. Fig. 1(b) shows an example of reproducing different lung
shapes by changing the weights of the principal components of the statistical
shape models as well as the camera parameters. The lung statistical shape
model is a statistical model with low-dimensional parameters describing the
displacement vectors per vertex obtained from the vertex-correlated lung shape
meshes obtained by mesh deformation alignment[10][11] between $M$ and
$\hat{M}$. Owing to the collapse deformation, the lung shape at the time of
air inclusion, which can be obtained before surgery, differs significantly
from the lung shape at the time of collapse during surgery. Therefore, the
simulated image reproduces lungs with various appearances and solves the
problem pertaining to a lack of data, which is a persistent issue in the field
of medical imaging. Fig. 1(c) shows the results of manually adjusting each
parameter to correspond to the surgical scene under the supervision of a
surgeon. Since the lung surface is transparent and the tumor location,
bronchus, and vascular structures can be comprehended, it can be used for
surgical assistance. If the parameters of statistical shape models can be
automatically calculated using the endoscopic images and the two can be
aligned, the shape of organs and anatomical structures can be identified
during surgery.

In endoscopic surgery, several situations are considered for the endoscopic
images and target organs to be captured during an operation. In this study, we
considered the following situations that are commonly observed in most
surgeries.

  * •

The initial shape $M$ of the organ to be operated on is obtained from the
3D-CT measured before the operation.

  * •

A part of the surface of the organ can be observed in the endoscopic image
$I_{R}$, but more than 50 $\%$ of the entire organ is occluded.

  * •

Information on the position of the endoscope camera and the direction of the
eye are not available.

  * •

The organ is deformed during observation, but its deformation is statistically
predictable.

In this study, we aimed to obtain the 3D shape of an organ under surgery
directly from the endoscope image via virtual learning using a simulated image
$I_{S}$ that can be generated from the statistical shape models under the
above conditions, and using image translation with common latent variables
between the endoscope image $I_{R}$ and the simulated image $I_{S}$. The
images used in this study were standardized to 180×120 pixels, and 10,000
simulated images generated for each case were used for training.

###  2.2 Framework

Fig. 2 shows the proposed framework, which consists of a VAE to achieve image
transformation with common latent variables between real and simulated images,
a convolution neural network (CNN) to learn the relationships between
synthetic images and model parameters. In the VAE, the real $I_{R}$ and the
simulated $I_{S}$ images, which have different features but are regarded as
equivalent by humans, are transformed into an identical synthetic image. Then,
in the CNN, we perform self-supervised virtual learning, using the synthetic
image and the initial shape $M$ of the organ as input, and outputting the
deformation mesh and camera parameters. The details of the VAE are described
in Section 2.3, and the details of the CNN are described in Section 2.4.

Figure 2: Proposed Framework. Variational Autoencoder (VAE) trains the common
latent variables in the real image $I_{R}$ and the simulated image $I_{S}$,
while the convolutional neural network (CNN) trains the relationship between
the synthetic image and the model parameters $\theta$.

###  2.3 Image translation

We propose a transformation function $f$ that translates real and simulated
images, which have different detailed features but are regarded as equivalent
by humans, into an identical synthetic image. In this study, the input images
to VAE are the real image $I_{R}$ and the simulated image $I_{S}$, which is an
artificial representation of $I_{R}$. The transformation function $f$ plays
the role of translating $I_{R}$ to $I_{S}$ and reconstructing $I_{S}$ to
$I_{S}$.

VAE consists of an encoder, which converts input data $x$ into latent
variables $z$, and a decoder, which reconstructs the same input data $x$ from
the latent variables $z$. The encoder and decoder in our study consists of
three fully connected layers. Let $\phi$ and $\theta$ be the parameters in the
neural networks of the encoder and decoder, and let $q_{\phi}(z|x)$ and
$p_{\theta}(x|z)$ be the probability distributions of the encoder and decoder,
respectively. $p(x)$ is a Gaussian distribution and the encoder learns the
parameter $\phi$ for $q_{\phi}(z|x)$, which generates parameters $\mu$ and
$\sigma$ of $p(x)$, and the decoder learns the parameter $\theta$ for
$p_{\theta}(x|z)$, which reconstructs $x$ from $z$.

Since $I_{R}$ and $I_{S}$ are compressed into a common latent variable, and
$f(I_{R})$ and $f(I_{S})$ are output from the common latent variable, we
assumed that it would be possible to reconstruct the identical synthetic image
from the common latent variable; therefore, we adopted VAE in this study. Let
$g$ be the transformation function of VAE in the usual learning method[12].
When the real image $I_{R}$ is input into the variable function $g$, the
relationship between the input $I_{R}$ and the output $g(I_{R})$ is
$g(I_{R})\approx I_{R}$. When the simulated image $I_{S}$ is input into the
variable function $g$, the relationship between the input $I_{S}$ and the
output $g(I_{S})$ is $g(I_{S})\approx I_{S}$. By contrast, in this study, when
the real image $I_{R}$ is input into the transformation function $f$ of VAE,
the relationship between the input $I_{R}$ and the output $f(I_{R})$ is
$f(I_{R})\approx I_{S}$. When the simulated image $I_{S}$ is used as input,
the relationship between the input $I_{S}$ and the output $f(I_{S})$ is
$f(I_{S})\approx I_{S}$. By using the transformation function $f$ to transform
two images with different features into the identical synthetic image, we
assumed it would be possible to learn that $I_{R}$ and $I_{S}$ are equivalent.
In addition, synthetic images $f(I_{R})$ and $f(I_{S})$, which are transformed
from $I_{R}$ and $I_{S}$, are highly similar to each other.

In VAE, the error between the input image to the encoder and the image
reconstructed by the decoder must be minimized. In addition, a regularization
term that minimizes the error between the distribution of $q_{\phi}(z|x)$
estimated by the encoder and $p_{\theta}(z)$, which is the prior distribution
of $z$(both Gaussian), is required. Therefore, in VAE training, we introduce a
loss function that is a weighted sum of the reconstruction error shown in the
first term of Eq. (1) and the Kullback-Leibler Divergence (KLD) shown in the
second term of Eq. (1).

$\mathcal{L}=\|x-\hat{x}\|^{2}_{2}+\lambda_{KL}D_{KL}(q(z|x)\parallel p(z))$
(1)

###  2.4 Virtual learning

In this section, we describe virtual learning using the simulated images that
can be generated statistically. In statistical shape models, the collapse
deformations assumed during surgery can be generated from the preoperative CT
models of individual patients before surgery by changing the weights related
to the lung collapse volume. The organ deformations are based on the average
collapsed deformations by learning. The vertex coordinates
$\upsilon_{i}^{{}^{\prime}}$ of the generated $\hat{M}$ can be obtained by the
Eq. (2).

$\upsilon_{i}^{{}^{\prime}}=\upsilon_{i}+\bar{u}+\omega\cdot u_{i}$ (2)

where $\upsilon_{i}$ is the position of the vertex of the lung shape during
aeration, $\bar{u}$ is the average displacement due to collapse deformation,
and $u_{i}$ is the displacement vector obtained from the lung statistical
shape model. $\omega$ is the weight of the displacement, where $\omega<1.0$
represents a squeezed shape compared with the average displacement due to
collapse deformation, and $\omega>1.0$ represents a bulged shape compared with
the average displacement due to collapse deformation. In this study, the
parameters were varied with a width of 0.1 around the average displacement due
to collapse deformation ($\omega=1.0$).

Next, $\hat{M}$ is rendered with the set camera parameters to generate the
simulated image $I_{S}$. The camera parameters that are expected to change
during the surgery are the camera position in three dimensions and the focus
point in two dimensions, which is a total of five dimensions. In this study, a
large number of simulated images IS were generated by changing the parameters
with a width of 15 mm in the depth direction of the camera position, 25 mm in
the other two directions, 15 mm in the horizontal and vertical directions of
the focus point for the right lung case, and 10 mm in the three directions of
the camera position and two directions of the focus point mentioned above for
the left lung case, centering on the camera parameters assumed in advance by
the surgeon. To implement the image transformation described in Section 2.3
accurately, in this study, $I_{S}$ was converted into a binary image of the
lung region and the background, and then into a lung contour image by
morphological transformation[13].

### Reconstruction Error

Since the purpose of this study was to obtain the 3D shape of an organ
adjusted to an endoscope image, it was desirable to minimize not only the
difference between the estimated shape and the target shape but also the
deviation of the position and orientation of the obtained shape. In the
proposed method, Reconstruction Error, defined in Eq. (3), is introduced as a
loss function to calculate the error between mesh vertices after projection.

$\mathcal{L}_{R}=\frac{1}{n}\sum_{i=1}^{n}\|M_{\theta}\upsilon_{i}-M_{\hat{\theta}}\hat{\upsilon_{i}}\|$
(3)

where $M_{\theta}$ and $M_{\hat{\theta}}$ are the target and estimated
projection matrices, respectively, $\upsilon_{i}$ and $\hat{\upsilon_{i}}$ are
the target and estimated mesh vertex points, respectively, and $M$ is a
$4\times 4$ matrix used for the perspective projection transformation in
generating the rendering image of the organ mesh with camera parameters.

### Parameter Loss

Since the space in the thoracic cavity is limited, the motion of the
endoscopic camera is restricted within a certain range. To constrain the range
of variations in the camera parameters, the proposed method introduces the
Parameter Loss defined in Eq. (4) as a loss function to calculate the error
between the estimated value $\hat{\theta}$ and the true value $\theta$ of the
camera parameters and the displacement $\omega$ in the collapse deformation of
the lung.

$\mathcal{L}_{P}=\frac{1}{n}\sum_{i=1}^{n}|\theta_{i}-\hat{\theta}_{i}|^{2}_{2}$
(4)

The loss function $\mathcal{L}$ for the entire virtual learning is defined as
a weighted linear sum of two loss functions as Eq. (5):

$\mathcal{L}=\mathcal{L}_{R}+\lambda_{P}\mathcal{L}_{P}$ (5)

##  3 Experiments

By comparing shape reconstruction methods, the effectiveness of the proposed
framework in organ shape reconstruction using a single endoscope image was
confirmed. The entire network was implemented in Python 3.6.8, Keras, and
TensorFlow GPU libraries. Virtual learning was performed using the Adam
optimizer with the following parameters: batch size = 60, number of epochs =
500, dropout rate = 0.5, and learning rate = $1.0\times 10^{-3}$ . A single
NVIDIA GeForce RTX 2070 GPU was used for training, and the virtual learning
took 1.4 hours.

The implementation of this framework enables us to output 3D shapes aligned to
the endoscope image from the initial shape. In this study, mean absolute error
(MAE), defined as the shape reconstruction error shown in Section 2.4, was
used as an evaluation metric for shape reconstruction accuracy. MAE evaluates
the average error between 3D shape vertices in the camera coordinate system.
Based on the comparison conditions used in the study by Tang et al., this
study compared the results using only real images [7], using simulated images
generated from statistical shape models (virtual learning) [7], and synthetic
images obtained by transforming simulated images using VAE(the proposed
method). As a combination of the loss functions, $\lambda_{KL}=5$,
$\lambda_{P}=0.5$ were adopted based on the results of several sets of
experiments.

Table 1: Shape reconstruction error for the proposed method (Proposed), real image only (Real), and virtual learning (Virtual). The mean $\pm$ standard deviation of MAE [$\mathrm{mm}$]. | Real |  | Virtual |  | Proposed |
---|---|---|---|---|---|---
| 75.7 $\pm$ 33.6 |  | 27.2 $\pm$ 15.5 |  | 22.6 $\pm$ 12.7 |

Cross-validation was performed for all 12 cases except for one test case.
Table 1 shows the evaluation results of the shape reconstruction in 11 cases,
excluding Case 12 owing to the small amount of data. As a result, it can be
confirmed that the shape reconstruction error of this method is significantly
smaller than that of the method based on Tang et al.’s idea (one-way analysis
of variance, ANOVA; $p<0.05$ significance level).

Figure 3: Shape reconstruction results for the proposed method (Proposed),
real image only (Real), and virtual learning (Virtual). The proposed method
can more accurately estimate the tumor position, vascular structure, and
posture.

Fig. 3 shows the visualization of the shape reconstruction results. In Fig.
3(a), the tumor location and vascular structure of the 3D image of the
proposed method are closer to the target image than those of other methods.
The MAE of the method using real images only is 112.6 mm, virtual learning is
18.7 mm, and the proposed method is 9.7 mm. Fig. 3(b) shows that the output
image of the proposed method is closer in appearance to the target image. Fig.
3(c) shows that the pose of the 3D image is closer to the target image than
the other methods; however, there is a shift in the tumor position. Virtual
learning can reduce the vertex-level error. In addition, the proposed image
transformation method enables the constraint of camera parameters to be
applied to the real image, and thus an appropriate pose can be estimated.

##  4 Conclusion

In this paper, we proposed a framework of self-supervised virtual learning for
the problem of organ shape reconstruction from endoscopic camera images.
Virtual learning based on simulated images generated from statistical shape
models was used to cope with the small amount of training data, and the
transformation to a common generated image with common latent variables
between real and simulated images was realized using VAE. However, the
variations in real images are limited, and image transformation results are
varied. In the future, we plan to improve the performance of image
transformation.

## References

  * [1] Widya, A. R., Monno, Y., Okutomi, M., Suzuki, S., Gotoda, T. and Miki, K.: Stomach 3D Reconstruction Using Virtual Chromoendoscopic Images, IEEE Journal of Translational Engineering in Health and Medicine, Vol. 9, pp. 1–11 (2021).
  * [2] Koo, B., Ozgur, E., Roy, B. L., Buc, E. and Bartoli, A.: Deformable Registration of a Preoperative 3D Liver Volume to a Laparoscopy Image Using Contour and Shading Cues, International Conference on Medical Image Computing and Computer Assisted Intervention (MICCAI), pp. 326–334 (2017).
  * [3] Wu, S., Nakao, M., Tokuno, J., Chen-Yoshikawa, T. and Matsuda, T.: Reconstructing 3D Lung Shape from a Single 2D Image during the Deaeration Deformation Process using Model-based Data Augmentation, 2019 IEEE EMBS International Conference on Biomedical Health Informatics (BHI), pp. 1–4 (2019).
  * [4] Nakao, M., Nakamura, M. and Matsuda, T.: Imageto-Graph Convolutional Network for Deformable Shape Reconstruction from a Single Projection Image, International Conference on Medical Image Computing and Computer Assisted Intervention (MICCAI), pp. 259–268 (2021).
  * [5] Modrzejewski, R., Collins, T., Bartoli, A., Hostettler, A. and Marescaux, J.: Soft-Body Registration of Preoperative 3D Models to Intra-operative RGBD Partial Body Scans, International Conference on Medical Image Computing and Computer Assisted Intervention (MICCAI), pp. 39–46 (2018).
  * [6] Ketcha, M. D., De Silva, T., Uneri, A., Jacobson, M. W., Goerres, J., Kleinszig, G., Vogt, S., Wolinsky, J.-P. and Siewerdsen, J. H.: Multi-stage 3D-2D registration for correction of anatomical deformation in image-guided spine surgery, Phys Med Biol, Vol. 62, No. 11, pp. 4604–4622 (2017).
  * [7] Tang, Z., Chen, K., Pan, M., Wang, M. and Song, Z.: An Augmentation Strategy for Medical Image Processing Based on Statistical Shape Model and 3D Thin Plate Spline for Deep Learning, IEEE Access, Vol. 7, pp. 133111-133121 (2019).
  * [8] Kingma, D. P. and Welling, M.: Auto-Encoding Variational Bayes, Porc. 2nd International Conference on Learning Representations (ICLR), (2014).
  * [9] Sengupta, S., Athwale, A., Gulati, T., Zelek, J. and Lakshminarayanan, V.: FunSyn-Net: enhanced residual variational auto-encoder and image-to-image translation network for fundus image synthesis, Porc. SPIE Medical Imaging, 113132M, (2020).
  * [10] Maekawa, H., Nakao, M., Mineura, K., Chen-Yoshikawa, T. F. and Matsuda, T.: Model-based registration for pneumothorax deformation analysis using intraoperative cone-beam CT images, Proc. 42nd Annual International Conference of the IEEE Engineering in Medicine and Biology Society (EMBC), pp. 5818–5821 (2020).
  * [11] Nakao, M., Maekawa, H., Mineura, K., Chen-Yoshikawa, T. F., Date, H. and Matsuda, T.: Kernel-based modeling of pneumothorax deformation using intraoperative cone-beam CT images, Proc. SPIE Medical Imaging, 115980P, (2021).
  * [12] Wan, Z., Zhang, Y. and He, H.: Variational autoencoder based synthetic data generation for imbalanced learning, 2017 IEEE Symposium Series on Computational Intelligence (SSCI), pp. 1-7 (2017).
  * [13] Serra, J.: Introduction to mathematical morphology, Porc. Computer Vision, Graphics, and Image Processing, Vol. 35, No. 3, pp. 283–305 (1986).


# Self-crossing geodesics

Anton Petrunin

##  1\. Introduction

The rubber band on the picture is pulled around a pebble, and it crosses
itself at several points.

The combinatorics of self-crossings can be described by a closed plane curve —
it is the rubber band in a parametrization of the surface with one point
removed. For example, if you could turn the pebble around you would see that
the self-crossings are described by the plane curve on the right diagram.

We assume that the surface of the pebble is strongly convex, smooth, and
frictionless; in this case, the rubber band models a closed geodesic. Suppose
that we are interested in possible patterns of self-crossings; more precisely:

_What are the possible combinatoric types of self-crossings of a closed
geodesic on a strongly convex smooth closed surface?_

Consider the six possible patterns with three double crossings. Configurations
1, 2, 3, and 4 can be realized as mirror-symmetric geodesics on mirror-
symmetric surfaces; the projections on the plane of symmetry are sketched.

Further, we will discuss _forbidden_ configurations; that is, configurations
that cannot appear for closed geodesic. These are configurations 5 and 6.

This question is a good exercise — it could be explained to anyone, but an
answer requires a considerable part of the theory. The reader is welcome to
check that there are no forbidden patterns with less than 3 double crossings
and try the cases with more self-crossings (check [2, Figures 15–17]). By the
way, _is there an algorithm for solving the general case?_ Our question is
closely related to the so-called _flat knot types_ of geodesics; see [1] and
the references therein.

In what follows, we discuss the Gauss–Bonnet formula as well as the
Alexandrov–Toponogov theorem and apply them to forbid configurations 5 and 6.
These theorems are covered in our textbook [3] which I like, altho they are
treated in plenty of other places as well.

## Gauss–Bonnet and no 5

Suppose that $\Delta$ is an $n$-gon with geodesic sides in a surface $\Sigma$.
Recall that by the _Gauss–Bonnet formula_ the sum of the external angles of
$\Delta$ equals

$2{\hskip 0.5pt\cdot\nobreak\hskip 0.5pt}\pi{\hskip 0.5pt\cdot\nobreak\hskip
0.5pt}\chi(\Delta)-\int\limits_{\Delta}K,$

where $\chi(\Delta)$ denotes the _Euler characteristic_ of $\Delta$ and $K$,
the Gauss curvature of $\Sigma$.

Further, we assume that $\Sigma$ is a closed strongly convex surface. In this
case,

  * $\diamond$

$\Sigma$ has strictly positive Gauss curvature;

  * $\diamond$

$\Sigma$ is homeomorphic to the sphere and therefore $\chi(\Sigma)=2$;

  * $\diamond$

$\Delta$ is homeomorphic to the disc and therefore $\chi(\Delta)=1$.

It follows that the sum of the internal angles of $\Delta$ is lager than
$(n-\nobreak 2){\hskip 0.5pt\cdot\nobreak\hskip 0.5pt}\pi$. In particular, if
$\Delta$ is a triangle with angles $\alpha$, $\beta$, and $\gamma$, then

$None$ $\alpha+\beta+\gamma>\pi.$

The Gauss–Bonnet formula can be applied to the whole surface; it implies that
the integral of Gauss curvature along $\Sigma$ is exactly $4{\hskip
0.5pt\cdot\nobreak\hskip 0.5pt}\pi$.

No 5 is forbidden. Suppose there is a geodesic with self-crossings as on the
diagram; it divides the surface $\Sigma$ into one triangle, say $\Delta$, one
hexagon, and three monogons. Denote by $\alpha$, $\beta$, and $\gamma$ the
internal angles of $\Delta$.

Note that three monogons have internal angles $\alpha$, $\beta$, and $\gamma$.
By Gauss–Bonnet, the integral of Gauss curvature along each monogon is
$\pi+\nobreak\alpha$, $\pi+\nobreak\beta$, and $\pi+\nobreak\gamma$
respectively. By $({*})$ the integral of Gauss curvature along the three
monogons exceeds $4{\hskip 0.5pt\cdot\nobreak\hskip 0.5pt}\pi$. But $4{\hskip
0.5pt\cdot\nobreak\hskip 0.5pt}\pi$ is the integral of Gauss curvature along
the _whole_ surface — a contradiction.

## Alexandrov–Toponogov and no 6

Let $\Delta$ be a geodesic triangle with angles $\alpha$, $\beta$, and
$\gamma$ on the surface $\Sigma$. Assume that the sides of $\Delta$ are
length-minimizing among the curves _in_ $\Delta$ with the same endpoints, then
the inequality $({*})$ can be made more exact.

Namely consider the so-called _model triangle_ $\tilde{\Delta}$ of $\Delta$;
that is, $\tilde{\Delta}$ is a plane triangle with equal corresponding sides.
Since the sides are length-minimizing, they satisfy the triangle inequality;
therefore the model triangle is defined.

Denote by $\tilde{\alpha}$, $\tilde{\beta}$ and $\tilde{\gamma}$ the angles of
$\tilde{\Delta}$ respectively. Then

$None$
$\alpha>\tilde{\alpha},\qquad\beta>\tilde{\beta},\qquad\text{and}\qquad\gamma>\tilde{\gamma}.$

Since $\tilde{\alpha}+\tilde{\beta}+\tilde{\gamma}=\pi$, this inequality
implies $({*})$.

The inequality $({*}{*})$ easily follows from the proof of the
_Alexandrov–Toponogov theorem_. The latter implies that $({*}{*})$ holds for
triangles with length-minimizing sides in the _whole_ surface. The proof is
left as an exercise for those familiar with the Alexandrov–Toponogov theorem;
others may simply accept it as true.

No 6 is forbidden. Suppose that such a geodesic $\xi$ exists; assume that its
arcs and angles are labeled as in the leftmost part of the diagram below.
Applying the Gauss–Bonnet formula to the quadrangle and pentagon that $\xi$
cuts from the surface, we get that

$None$ $2{\hskip 0.5pt\cdot\nobreak\hskip 0.5pt}\alpha<\beta+\gamma,\quad
2{\hskip 0.5pt\cdot\nobreak\hskip 0.5pt}\beta+2{\hskip
0.5pt\cdot\nobreak\hskip 0.5pt}\gamma<\pi+\alpha,\quad\text{and
therefore}\quad\alpha<\tfrac{\pi}{3}.$

Consider the part of $\xi$ without the arc labeled by $a$. It cuts from the
surface a pentagon $\Delta$ with sides and angles shown in the middle part of
the diagram.

Let us add additional vertices on the sides of $\Delta$ so that each side
becomes length-minimizing in $\Delta$. Choose a vertex of $\Delta$ and join it
by shortest paths in $\Delta$ to every other vertex. Consider a model triangle
for each triangle in the obtained subdivision of $\Delta$; the model triangles
lie in the plane and we suppose that they share sides as in $\Delta$. By the
comparison inequality $({*}{*})$, the angles of the model triangles do not
exceed the corresponding angles of the original triangle. Therefore, the model
triangles form a convex plane polygon, say $\tilde{\Delta}$, such that

  * $\diamond$

The five angles of $\tilde{\Delta}$ that correspond to the angles of $\Delta$
do not exceed those.

  * $\diamond$

Each side of $\tilde{\Delta}$ equals to the corresponding small side of
$\Delta$.

It remains to show that no convex plane polygon meets these two conditions.

Let us orient the sides $\tilde{\Delta}$ counterclockwise; denote the obtained
vectors by $s_{1},\dots,s_{k}$. Note that the vectors $s_{i}$ point in the
complement of white sectors shown below with angles marked. The sum of the
magnitudes of the vectors in each black sector is also marked (each black
sector corresponds to a side of $\Delta$).

By
$(\smash{\raisebox{-2.15277pt}{\begin{tabular}[]{@{}cc@{}}\lx<EMAIL_ADDRESS>*&*\end{tabular}}})$ we can choose a vector $r$ as on the diagram, so that
$\varphi>\nobreak\tfrac{\pi-\beta}{2}$ and $\psi>\tfrac{\pi-\gamma}{2}$. Note
that for any unit vectors $u$, $v$, $v^{\prime}$, $w$, and $w^{\prime}$ in the
marked black sectors, we have

$\begin{aligned} \langle r,u\rangle<0,\quad\langle r,v\rangle<0,\quad\langle
r,w\rangle&<0,\\\ \langle r,v\rangle+\langle
r,v^{\prime}\rangle<0,\quad\langle r,w\rangle+\langle
r,w^{\prime}\rangle&<0.\end{aligned}\qquad\qquad\qquad\qquad\qquad\qquad\qquad$

It follows that

$\langle r,s_{1}\rangle+\dots+\langle
r,s_{k}\rangle<0.\qquad\qquad\qquad\qquad\qquad\qquad$

But the vectors $s_{i}$ circumambulate $\tilde{\Delta}$; so, the sum has to
vanish — a contradiction.

Acknowledgments. I want to thank Arseniy Akopyan, Maxim Arnold, Serge
Tabachnikov, David Kramer, and Sergio Zamora Barrera for their help. This work
was partially supported by the NSF grant DMS-2005279, and Simons Foundation
grant no 584781.

## References

  * [1] S. Angenent  “Curve shortening and the topology of closed geodesics on surfaces”  In _Ann. of Math. (2)_ 162.3, 2005, pp. 1187–1241  DOI: 10.4007/annals.2005.162.1187
  * [2] V.. Arnold  “Remarks on the enumeration of plane curves”  In _Topology of real algebraic varieties and related topics_ 173, Amer. Math. Soc. Transl. Ser. 2, 1996, pp. 17–32  DOI: 10.1090/trans2/173/03
  * [3] A. Petrunin and S. Zamora Barrera  “What is differential geometry: curves and surfaces”, 2020