{'fraction_non_alphanumeric': 0.030138841855739926, 'fraction_numerical': 0.006137825939722316, 'mean_word_length': 4.813238188976378, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Robot Calligraphy using Pseudospectral Optimal Control in Conjunction with a Novel Dynamic Brush Model

Chinese calligraphy is a unique art form with great artistic value but difficult to master. In this paper, we formulate the calligraphy writing problem as a trajectory optimization problem, and propose an improved virtual brush model for simulating the real writing process. Our approach is inspired by pseudospectral optimal control in that we parameterize the actuator trajectory for each stroke as a Chebyshev polynomial. The proposed dynamic virtual brush model plays a key role in formulating the objective function to be optimized. Our approach shows excellent performance in drawing aesthetically pleasing characters, and does so much more efficiently than previous work, opening up the possibility to achieve real-time closed-loop control.


I. INTRODUCTION
Making robots write beautiful calligraphy is difficult, as learning and mastering this art form takes humans years of practice. Chinese characters are complex and a calligraphy brush is difficult to manipulate properly. Hence, making a robot achieve comparable results is a worthwhile endeavour, both to expand robot's capabilities into art as well as push our understanding on how to manipulate deformable brushes.
Most relevant research in this area adopts either a learningbased method or an approach based on trajectory optimization. The former includes learning from demonstration [1], [2], or learning from visual feedback given a reference image [3]. By using learning, one can eliminate the difficulty of modeling the behavior of a real calligraphy brush. However, learning methods have a large training cost and may not generalize well to previously unseen characters. On the other hand, trajectory optimization methods do not face these problems: they simulate the writing behavior of an actual brush, and then search for an optimal trajectory for the robot to execute [4], [5]). The difficulty there is that most simulated brush models [5], [6] do not account for the complex ways a brush deforms during the writing process. Being able to capture the complexity of a deformable brush has an important influence on the final performance.
In this paper, we propose a novel virtual brush model intended to capture the dynamics of an actual calligraphy brush, which is an improved version based on an earlier virtual brush model by Kwok et al. [4], which we use as a baseline. We improve on this model by explicitly modeling the brush dynamics, while retaining a much simpler structure compared to the elaborate models in other work [7], [8]. We  To do so, our second contribution is an efficient trajectory optimization-based method based on pseudospectral optimal control [9], [10] that achieves fully automatic writing of Chinese characters given a character's unicode. Pseudospectral methods are based on Legendre or Chebyshev polynomials, which are excellent at approximating and representing continuous trajectories and controls. Our method converges quickly and efficiently even when using many control points.
Finally, we exploit the existence of vector-based character databases such that, given a character's unicode, we immediately have access to the individual character strokes as well as the stroke order. We optimize for each stroke separately, and to obtain even faster convergence we initialize the trajectory for each stroke by a stroke skeleton which is also available in these databases.
Our proposed continuous, nonlinear optimization framework differs from previous work [4], [5], [6] which uses heuristic optimization methods known for slow convergence rates and a high computational cost. Hence, our method can potentially be used in a closed-loop control system, which we intend to pursue in future work. The proposed end-toend system displays excellent writing performance, although there is definitely room for growth.

II. RELATED WORK
There is much work in the robotics community that aims to create art using robots, such as painting and drawing [11], [12], [13], sculpture [14], graffiti [15], etc. Below we focus our discussion on robot calligraphy, most algorithms using a calligraphy brush can be categorized as either learning or trajectory optimization-based.

A. Calligraphy robots using learning-based methods
Examples of simple learning-based methods include Sun et al. 's learning from demonstration [16] [1], Mueller et al. 's trial-and-learn iterative learning method [3]. Some more advanced learning algorithms such as RNN [17], generative adversarial networks [18], deep reinforcement learning [19], and local and global learning models [2] were also explored. However, these methods require many iterations of training to achieve good performance, and have difficulty generalizing to new and/or complicated characters.

B. Virtual brush models
Trajectory optimization-based algorithms mainly rely on virtual brush models, which can be divided into two categories: physics-based models and data-driven models.
Physics-based virtual brush models strive to simulate the physical dynamics of a real brush from experimental observation [20], [21], [22], [23] or physical laws [8], [24]. Strassmann [25] proposes an initial design featuring four basic parameters of a hairy brush based on bristles. Wong et al. [26] propose to use a cone to represent the bundle of the brush and use the cross-section of the cone, an ellipse, to represent the footprint of the brush. Xu et al. [7] propose a detailed virtual brush model for use in synthetic imagery, with good results in creating realistic-looking simulations. However, obtaining and fitting good parameters to the models above is difficult. As such, we propose a virtual brush with an easy structure to fit and implement, geared towards real-time trajectory optimization.
Data-driven virtual brush models are created by measuring and recording actual brush footprints on writing surfaces. Kwok et al. propose a simple virtual brush which draws droplet shapes whose size proportional to the writing height [4]. In their later work [6], they use a camera placed below the writing plane to collect footprints during the writing process. Lam et al. [5] define their writing mark as a polygon connected by eight points and fit their position parameters with the collected footprints. With an eye for efficiency, Baxter et al. [27] build a deformation table to speed up the associated computation while still being able to simulate complex effects.

C. Stroke extraction
Stroke extraction involves separating a character into its comprising strokes and is difficult to do with good accuracy when analyzing only the pixels of an image. There are three main categories of stroke extraction methods: skeletonbased [28], [29], region-based [30], and contour-based [31], [32], [33]. Most of these methods are complex and may not achieve good results as compared to the ground-truth stroke segmentation. As such, we propose using vector-based character database, which provides a quick and accurate way to extract strokes, as well as stroke order.

D. Optimization methods
As mentioned above, Kwok et al. [4], [6] propose using genetic algorithm for direct character-level optimization guided by their brush models. However, their reported computation times are long, and the strokes have to be manually separated. Lam et al. [5] minimize the width difference of the strokes between reference images and a simulated image written by the virtual brush as it moves along the middle axis of each stroke. However, their method is sensitive to small variations in stroke images, and so the results suffer from a loss of smoothness. One of our main contributions below is the use of methods from optimal-control to achieve the desired stroke trajectories.

III. METHODOLOGY OVERVIEW
An overview of our methodology is given in Fig. 2 on the next page. Given a character's unicode we find the corresponding SVG image in a vector-based database, which already comes segmented into the different strokes. For each stroke, we compute initial trajectories from the stroke skeletons, and then optimize from there until our simulated image matches the reference images rendered from the SVG file. At the core of the simulation is a virtual brush model, which we discuss at length in the section below, followed by a discussion of our optimization framework.

IV. TWO VIRTUAL BRUSH MODELS
There are two virtual brush models that we have experimented with in order to simulate the calligraphy writing process: a simple virtual brush and a dynamic virtual brush. Both generate good results, but the dynamic virtual brush is more sophisticated and better at closing the "sim2real" gap.

A. Simple virtual brush model
The simple virtual brush model is similar to Kwok et al. 's work [4] and is used as a baseline comparison in our results. Given the height z of the brush, this model simply fills a circular patch with simulated ink, where the radius of the circle is proportional to z. Note that Kwok et al. also proposed using more complex templates (such as droplets), but we get that behavior for free by modeling the dynamics in our model below.

B. Dynamic virtual brush model
The dynamic virtual brush model has two components: a drawing component, and a dynamic update component. The drawing component describes how the brush leaves a mark on paper depending on its parameters. The updating component then describes how the brush parameters are updated due to deformations when executing an open loop trajectory [x(t), y(t), z(t)] T . Note that we align the task coordinate frame for the robot + brush such that the z = 0 plane corresponds to the plane of the paper.  1) Drawing Component: The dynamic virtual brush model has four state parameters governing its drawing footprint: width w, drag d, offset o, and orientation θ. As shown in Fig. 3, width w and drag d define the size of a brush mark, the offset o simulates the deviation of the brush mark from the center of the vertical brush handle due to bending of the brush hairs when applying pressure, and the orientation θ describes the direction of the writing mark. Together with the end-effector coordinates X .
= (x, y, z), these four parameters define the 7-dimensional brush state X .
Given the state X our model then computes the "root" location of the brush mark as where π(X) is the projection of X in the z = 0 plane and V transforms the scalar orientation θ into a three-dimensional unit vector in that same plane. The shape of the brush mark is parameterized by a quadratic curve which is uniquely determined by the width w and the drag d. Some more specialized brush behavior such as intentional hair-splitting is reserved from more artistic calligraphy styles and are not modeled here.
2) Dynamic Update Component: We adopt a simple discrete state dynamical model described at a high level by where U i = (v xi , v yi , v zi ) is the control at the i th time instance. The dynamics model is given below and was formulated based on careful observation and experimentation.
Above Width(z i+1 ) and Drag(z i+1 ) respectively compute a new value for width and drag, given the updated end-effector height z i+1 . The constants K w and K d model a small amount of inertia in those parameters -we have used 0.02 in our experiments -modeling the fact that the deformations happen gradually and steadily. The orientation and offset updates are more complicated given that they are governed by friction, modeling the fact that the part of the brush in contact with the paper resists movement. In particular, the function Offset(z i+1 ) computes a value for the offset o when friction is overcome, and Offset'(X root i , X i+1 ) computes the value of o that makes the brush stay put. The min operator ensures that the brush stays put until it "snaps" to the smaller value Offset(z i+1 ) when the end-effector X i+1 moves too far out. The orientation θ i+1 is accordingly updates, as X root i+1 is a function of o i+1 . 3) System Identification: The width, drag, and offset functions above were experimentally determined in a separate calibration phase. We measured the width, drag, and offset of writing marks left by a real calligraphy brush, then fit a polynomial relationship between the parameters, width w, drag d, and offset o, and the change in height of the brush, z. We fit the relationships between the parameters width w, drag d, and offset d w.r.t. z by collecting brush footprint data at varying writing heights, shown in Figure 4. The figure shows that width and drag behave almost linearly, but that the offset o has a local maximum. We observe that moving the brush in one direction will cause the tip of the brush to gradually move towards the opposite direction.

A. Review of Pseudospectral Optimal Control
Below we formulate the calligraphy writing process as a trajectory optimization problem, and in particular we adopt some of the machinery from pseudospectral optimal control (PSOC) methods. As such, we briefly review pseudospectral optimal control methods in this section, closely following the exposition by Fahroo et al. [10].
A simplified version of an optimal control problem can be stated in terms of a cost function C and system dynamics f where the objective is to find the optimal control sequence u(t) that minimizes the cost function C. Above, x(t) represents the system's state trajectory. The basic idea in PSOC is to approximate the control trajectory u(t) and the state trajectory x(t) by a polynomial curve with unknown parameters, thereby transforming the original problem into a nonlinear programming problem. To this end, pseudospectral methods choose a specific set of points from the curve for interpolation. For example, in the case of Chebyshev pseudospectral methods, the interpolation points used are the Chebyshev-Gauss-Lobatto (CGL) points [34]: To recover the state x(t) at any arbitrary time t one can use barycentric interpolation, i.e. where We can now express the original dynamics equation with an approximation, where the objective function (10) can be discretized if necessary. The original problem is thus transformed to minimizing the cost C with respect to the two coefficient vectors X, U : representing the values of the state and controls, respectively, at the CGL points.

B. PSOC for Calligraphy
Below we apply these methods to trajectory optimization for open-loop control trajectories of a robot end-effector, with the goal of faithfully reproducing Chinese characters. PSOC methods are generally used for collocated optimal control where the system dynamics are enforced through specialized components of the cost function. However, because our controls are simply the trajectory velocity they can be calculated exactly from the end-effector trajectory, and hence we have no separate control parameters U .
The optimization for a character is decomposed into a series of trajectory optimization problems corresponding to individual strokes of the character. This greatly simplifies the process and is also more computationally efficient. In this, we are helped by the existence of vector-based character databases in which characters are stored as their decomposed individual strokes, as described next. To obtain reference images for each stroke, as well as initialize the nonlinear optimization, we exploit the existence of vector-based character databases. Vector-based images are advantageous because they store individual strokes as Bezier curves, which are themselves continuous polynomial curves defined by a set of control points. This makes the stroke extraction a trivial task. In our case, the database we choose is a Scalable Vector Graphics (SVG) database from MakeMeHanzi [35], and an example of extracted strokes can be seen in Fig. 5.

D. Stroke trajectory representation
The stroke trajectories will be used as the open-loop control trajectories for the robot to draw strokes, and we represent them as three-dimensional trajectories of the endeffector. Each of the x, y, z components are separately represented as 1-dimensional Chebyshev polynomial curves. The trajectories X with three dimensions x, y, z are expanded by interpolating the values at the Chebyshev-Gauss-Lobatto (CGL) points given by Eq. 12 and Eq. 14: Hence, the decision variables are the combination of the three sets of CGL points in the x, y, z dimensions:

E. Optimization for the stroke trajectories
The objective function C s for stroke s minimizes the sumsquared pixel difference between the image V (X s ) produced by simulating the drawing process, and a reference image I s created from the SVG representation. The optimal trajectory parameters X * s for stroke s are obtained as = arg min Above, V (·) is a function that represents the virtual brush simulation, taking a stroke trajectory with pseudospectral parameterization X s and drawing a stroke image according to the given trajectory. R 1 (X s ) is a regularization term to make sure the final trajectory has a tendency to be lift up when a stroke is finished. The optimization is conditioned on the initial brush state X init when the brush starts writing the stroke. In particular, if θ (d−1) last is the final orientation after writing the previous stroke it is constrained to be the initial direction parameter for the current stroke. Finally, we scale the z i to the same magnitude as x i and y i for the convenience of optimization.
In practice we try all polynomial orders between 3 and 8 for each stroke and pick the one with the smallest error, in order to adapt to strokes of varying length and complexity. Finally, a dense and even sampling of the polynomial is performed in V to execute a continuous stroke in simulation.
We use the Levenberg-Marquardt (LM) algorithm, a second-order trust-region method, to minimize Eq. 18. In short, LM switches between a gradient-based search and a second-order Gauss-Newton update by controlling a damping factor λ: The Jacobian matrix J is computed using numerical differentiation at each iteration. We use the the GTSAM library [36] to perform the optimization. GTSAM was originally created to solve simultaneous localization and mapping problems but has been used in many different contexts since, including motion planning [37], [38].

F. Trajectory Initialization Estimates
We use the skeleton of a stroke as our initialization for the x and y coordinates of our trajectory, while the initial z coordinates are set to a fixed value. From our observations, when people write calligraphy, they generally make the brush approximately follow the skeleton of the stroke while varying the height of the brush.
Because we start from a vector-based representation, extracting an initial trajectory for the individual strokes is much simplified. Even when starting from images, there are many good image-based skeleton extraction algorithms, e.g. the Chordal Axis Transform [39], [40], an example result of which is shown in Fig. 6a. In our case, we use the "animation path" provided by the database as the initial estimation for the 2D x i and y i sequence for simplicity. CGL points are sampled on each skeleton path to obtain the pseudospectral representation. Generating an easy estimate for z is not intuitive, and so we just initialize it with a constant sequence. Given the initial trajectory, we can simulate the image formation process using both the simple virtual brush and the dynamic virtual brush, as illustrated in Fig. 6b and Fig. 6c respectively. The position of the written mark from the dynamic virtual brush is different from its given trajectory, which is ignored by most previous research.

G. Sim2Real Analysis and Improvements
It is known that writing brushes need to periodically dip ink to be able to write. To avoid prediction error to accumulate as more strokes, especially complicated ones, are written, we propose to dip ink and restore the brush to its standard state each time after a stroke is written. That means, the θ (d−1) last in (18) becomes a constant value representing the initial states. After the brush writes some complicated strokes or performs some unexpected movement, the tip may become messy. At that time, predicting the state of the brush accurately and finding feasible trajectory to control is difficult and unreliable. As such, we carefully design control algorithms to help the brush dip ink. Given a circle ink stone, the brush is pushed down heavily at first to make the tip flat, and then it slowly moves to the edge of ink stone in different directions with gradually smaller extent. After that, the tip will generally restore the initial state with small error.

VI. RESULTS
Below we present the results of our approach, including photographs of characters that have been drawn by a Franka robot, using inverse kinematics provided by MoveIt! [41]. The characters are written using relatively slow end-effector velocity to prevent excessive jerk and vibrations.   7 shows a comparison between the written results generated by the simple virtual brush model and by the dynamic virtual brush model. The simple virtual brush model is able to generate good touching and brush outlines for some simple strokes (like dot) because of the good linear fitting model in the Sim2Real section, but not all (like the long hook). However, it can be seen that the dynamic virtual brush shows obvious advantages in terms of preserving the relative positions of strokes, since it models the deformation of the real brush. Periodically restoring the brush to its initial state helps the simple brush to behave better, but it still worse than the dynamic virtual brush model.
In Fig. 8, we show more extensive results for three different characters. Both the simulated and written images before and after optimization are shown for easy comparison. From the figure, we can see that the optimization achieves good performance for simulated images. However, because of a "sim2real" gap in the virtual brush model, the proposed method still has a ways to go in terms of approaching the smoothness and definition of detail displayed in the reference images, rendered from the vector-based character database.
A detailed error analysis shows that the optimization achieves almost perfect results in simulation, but that there still is a considerable sim2real gap. Fig. 9 shows the differences between reference images and both simulated and written characters, respectively before and after optimization. As in Figure 7, green and red pixels represent positive and negative differences. Column (c) shows that the optimization can achieve near-perfect results, where the only errors remaining stem from anti-aliasing, as well as a few small errors near the stroke endings. However, we see larger errors when these optimized results are executed on the robot (column e), even though the characters themselves look visually pleasing.  In this paper, we presented a trajectory optimization method to make robots write calligraphy, searching for openloop control trajectories by minimizing a simulation against a reference image. The proposed dynamic virtual brush simulation yields good results in simulation and in turn produces reasonable open-loop control trajectories for a real robot. However, from the results it is clear that the sim2real gap has not been fully closed, and it is an open question whether a closed loop strategy will ever yield master-level calligraphy. Hence, in future work we plan to investigate on both better brush models and the possibility of feedback control around optimized trajectories.