\section{Introduction}

Respect for tissue \cite{martin1997objective}, or force sensitivity \cite{goh2012global}, is considered an important skill for performing safe surgery and requires good control of applied forces. Thus, knowledge of the force exerted by a robotic system on the surgical environment is important during robot-assisted minimally invasive surgery (RMIS) to enable safe tissue handling. Force information can be used to provide haptic feedback to the surgeon, automatically and objectively evaluate their force sensitivity for training and credentialing purposes, and inform the decisions and movements of an autonomous agent.

Force information has been difficult to obtain for the above purposes in RMIS because there is no native distal force sensing in commercial RMIS systems. This is due in part to designers needing to meet the strict requirements for biocompatibility and sterilizability of RMIS instruments while ensuring cost-effectiveness \cite{EnayatiReview2016,HadiHosseinabadi2022}. While researchers have explored many approaches to developing force sensors that attempt to address the above constraints, none have gained commercial adoption. Furthermore, many designs contain complex electromechanical components that require specialized knowledge to manufacture, assemble, or integrate. This limits their adoption even in the research community.

To make up for the lack of feasible force sensing options for RMIS tools, researchers aiming to improve force sensitivity have often relied on existing general-purpose commercially available force sensors like those from ATI Industrial Automation (Apex, NC, USA) that are placed in or under the artificial tissue being manipulated. In such a set up, researchers are limited to RMIS studies using only a single end effector \cite{Bahar2020,Chua20201Dynamics} or measuring a single force value for both end-effectors \cite{Brown2017apms,sensubGaleazzi2022}. This approach prevents the study of studying bimanual force-critical tasks such as those shown in Fig.\,\ref{fig:concept} and thus limits applicability to real surgery.

 One approach that circumvents the need for end-effector force sensors is indirect force sensing. This has been explored using physics-based \cite{fontanelli2017modelling,Wang2019dynamic} or neural network models \cite{Yilmaz2020dynamic} of the robot to predict joint torques and using vision-based finite element \cite{Haouchine2018FEA} or deep learning methods \cite{aviles2014recurrent,aviles2016towards,Marban2019rnn,ChuaNeuralForce}. However, these approaches need to be trained or benchmarked against a ground truth. To achieve this, researchers have often used a single environmental force sensor like the ATI sensors noted above. For methods that rely entirely on the robot's internal state, this approach is feasible because each end effector can be trained separately. However, for methods that rely on measuring environmental changes, such as vision-based methods, this approach has limited applicability to bimanual manipulations where internal ground truth force data cannot easily be resolved.

In this work, we present the the design and characterization of a 3-degree-of-freedom (3-DoF) force sensor for RMIS research. By relaxing the strict constraints on size, biocompatibility, and sterilizability, we realize a sensor that is easy to manufacture and integrate into existing hardware while achieving satisfactory performance for a wide variety of research use cases. To capitalize on its manufacturability and modularity, we have open-sourced the design to enable researchers to adapt the sensor for their desired application both within RMIS and for other robotic applications where deploying commercial sensors is difficult due to size, cost, and customizability requirements.

\begin{figure}[!t]
    \centering
    \vspace{1.8mm}
    \includegraphics[width=\linewidth]{fig/concept_design.pdf}
    \caption{Concept renderings of the 3-DoF force sensor design in example use cases requiring bimanual manipulation. (a) View of sensorized patient-side manipulators during suturing. Close-up of sensorized forceps during (b) blunt dissection, (c) running the bowel, and (d) suturing.}
    \label{fig:concept}
\end{figure}

\section{Background}

Previous works have placed sensors at different locations of the RMIS instrument. These include the jaws, wrist, lower shaft, trocar, upper shaft, and the instrument base. These works have also employed various types of sensing technology, with metal strain gauges, capacitive sensors, fiber-Bragg gratings, and infrared (IR) light intensity measurement being the most common technologies employed.

Force sensors located at the jaw have been implemented as both 2-DoF \cite{Hong2012} and 3-DoF sensors \cite{Yu2018} using strain gauges and custom jaw flexures. Custom jaw flexures were also employed in \cite{Kim2018} for 3-DoF force and 2-DoF moment sensing using capacitive elements. Such jaw flexures allow forces to be sensed when any part of the tool tip interacts with the environment. This is important during tasks such as blunt dissection (Fig.\,\ref{fig:concept}b) and running the bowel (Fig.\ref{fig:concept}c) where the tip and back of the tool tip are used for manipulation. This is in contrast to approaches which place sensors on the grasping surface of each jaw and thus only allow forces to be sensed when the environment is grasped \cite{Kim2015,Dai2017}. By locating the sensors at the jaw, grip force can also be computed from the force measurements at each jaw.

Locating force sensors above the articulated tool wrist reduces the electromechanical integration complexity typical of sensors located at the jaw. In \cite{Li2017wrist}, a Stewart platform with strain gauges was used to measure 3-DoF forces at the articulated wrist of a custom RMIS tool. In \cite{Lee2015}, 6-DoF forces and moments were measured at the articulated wrist using capacitive sensors. Torque sensors were also embedded in the drive pulleys and used to both measure grip force and compensate for noise in the wrist force sensors due to drive cable actuation. Sensors have also been located on the lower shaft of RMIS tools with both \cite{Shahzada2016fbgshaft} and \cite{Du2022fbgshaft} using fiber-Bragg gratings to measure 2-DoF lateral forces.

Sensors have also been placed at the interface of the patient's body (at the trocar), or outside the body, for example, on the upper shaft or instrument base. In \cite{Kim2017trocar}, strain gauges were used at the trocar to measure 2-DoF lateral forces, while in \cite{Fontanelli2020trocar}, IR intensity measurement was used instead. At the upper shaft, \cite{hadi2021shaft} used IR intensity to measure 3-DoF forces and moments though for the three directions of force, accuracy metrics for only the two lateral directions of force were reported. At the instrument base and upper shaft, \cite{Novoseltseva2018} used strain gauges to measure 3-DoF forces, with force measurements along the main axis of tool showing poorer accuracy relative to those along the lateral directions.

\section{Methods}

\subsection{Target Design Requirements}

To realize a force sensor that is accessible to the research community, the design of the sensor should be easily manufacturable and integrated into existing RMIS tools. This makes sensors located on the upper shaft of the instrument base particularly suitable \cite{HadiHosseinabadi2022}. However, these designs typically lack accuracy along the main axis of the RMIS tool. Furthermore, they are unable to measure grip force, which can be useful for evaluating surgical skill or for providing feedback to improve tissue manipulation.

RMIS research is often performed on ex-vivo or dry lab tasks. This relaxes the requirements on biocompatibility, sterilizability, and size. Thus, placing sensors at the tool jaws does not require complex jaw designs and can utilize small-size commercial load cells. At the same time, locating the sensor at the jaw reduces measurement noise and enables more accurate measurement of force along the main axis of the tool. Additionally, the jaw sensor placement allows for straightforward grip force measurement.

Based on these considerations, we designed a jaw-mounted sensor that can be customized to suit different RMIS tools and different use cases. Our target use case of tissue manipulation requires that all parts of the jaw be able to sense force, and thus, unlike in \cite{Kim2015,Dai2017}, the sensing elements cannot be solely mounted on the grasping surfaces of each jaw.

Tissue manipulation forces can be up to 3.8\,N in the lateral direction, $-$10.3\,N in the axial compression direction, and 7.8\,N in the axial retraction directions \cite{Toledo1999}. However, there is a need to balance these requirements against the current capabilities of small size commercial load cells. Based on these considerations, our sensor target force ranges are $\pm3\,\text{N}$ in the lateral directions and $\pm5\,\text{N}$ in the axial direction. In terms of accuracy requirements, the average kinesthetic force just noticeable difference for the human hand is 12.5\% \cite{Vicentini2010}, and thus the sensor requires a minimum sensor accuracy of 0.375\,N in all directions for error imperceptibility.

\subsection{Electromechanical Design}
Based on the above design requirements, we designed a 3-degree-of-freedom force sensor located at the tool jaws. As shown in Fig.\,\ref{fig:sensor_mech}a, the sensor comprises five main parts: (1) the base, (2) two load cell arrays, (3) the sensing plate and rod, and (4) the jaw attachment. An optional strain relief bracket (also shown in Fig.\,\ref{fig:sensor_mech}a) can be added to help secure and route wires.

The base is 3D-printed in 6061 aluminum and its geometry can be modified to interface with different RMIS tool jaws. In this paper we present a design that interfaces with the da Vinci Surgical System (Intuitive Surgical, Inc., Sunnyvale, CA, USA) large needle driver jaw using a M2$\times$3 set screw. One of two load cell arrays is placed above the top surface of the base and is electrically isolated using Kapton tape.

\begin{figure}[!t]
    \centering
    \vspace{1.8mm}
    \includegraphics[width=0.75\linewidth]{fig/design_new.pdf}
    \caption[Electromechanical design of the force sensor.]{(a) Exploded view of the force sensor mounted to one jaw of the da Vinci large needle driver. (b) Arrangement of the load cells on the PCB of the load cell array. (c) Fully assembled force sensor.}
    \label{fig:sensor_mech}
\end{figure}

\begin{figure}[!t]
    \centering
    \includegraphics[width=0.75\linewidth]{fig/FBD.pdf}
    \caption[Front and side views of the force sensor with top views of the load cell arrays.]{Front and side views of the force sensor with top views of the load cell arrays. Load cells are numbered 1 through 8 corresponding to Eqns.\,(\ref{eqn:moment_x})\,--\,(\ref{eqn:force_z}).}
    \label{fig:FBD}
\end{figure}

Each load cell array is a 9.5$\times$8.5$\times$1.6\,mm 2-layer FR4 printed circuit board with four HSFPAR003A load cells (Alps Alpine, Japan) soldered along the perimeter (Fig.\,\ref{fig:sensor_mech}b). The load cells measure compression forces of up to 8\,N and rely on a piezoresistive full Wheatstone bridge which allows for good temperature stability. The bridge outputs are amplified using AD623 (Analog Devices, Norwood, MA, USA) instrumentation amplifiers with a gain of 21. This results in a sensor response of 3.063\,$\text{N}\,\text{V}^{-1}$. The amplified analog signals from each sensor were recorded on a PC using an Arduino Mega with serial communication at 125\,Hz.

The sensing plate and rod is machined out of 303 stainless steel for high stiffness. A second load cell array is placed on the top face of the plate in opposition to the first load cell array. The two load cell arrays and the sensing rod and plate are attached to the base using four M1.2$\times$12\,mm screws.

A jaw attachment, which replaces the original tool jaws for grasping, is machined out of 6061 aluminum and is attached to the sensing rod by a M2$\times$3\,mm set screw. The jaw attachment is interchangeable, allowing for researchers to machine different shapes to suit the task they are studying. Here we fabricated a generic shape for tissue retraction and palpation that has a height of 12\,mm. During manufacturing of the load cell arrays, there are small deviations in the heights of each load cell after soldering. Thus, we enabled consistent contact between the sensing plate and the individual load cells on each side by inserting metal shims. The sensors were preloaded up to a maximum of 1.5\,N. The fully assembled sensor is shown in Fig.\,\ref{fig:sensor_mech}c with overall dimensions of 9.5$\times$8.5$\times$23.8mm and a weight of 3.33\,g.

\subsection{Sensing Principle}
\label{sensingprinciple}
The sensing principle relies on the moment balance about the lateral axes of the device, henceforth referred to as the sensor x- and y-axis, and force balance in the main axis, henceforth referred to as the sensor z-axis. Assuming an interaction at the tip of the jaw attachment, and neglecting the contribution of shear forces to the moment balance, the resulting force and moment equations are

\begin{equation}
    M_x =  F_y H - \frac{Lc}{2} (v_2 + v_6- v_4 - v_8) \text{,}
\label{eqn:moment_x}
\end{equation}
\begin{equation}
    M_y =  F_x H - F_z D - \frac{Wc}{2} (v_1 + v_7- v_3 - v_5) \text{, and}
\label{eqn:moment_y}
\end{equation}
\begin{equation}
    F_z =  c (v_1+v_2+v_3+v_4-v_5-v_6-v_7-v_8)\text{,}
\label{eqn:force_z}
\end{equation}
where $H=15.85$\,mm, $D=5.50$\,mm,  $L=3.45$\,mm, and  $W=2.95$\,mm. $c=3.063\,\text{N}\,\text{V}^{-1}$ is the voltage change per unit force, and $v_i$ the voltage output of the $i\textsuperscript{th}$ sensor as labeled in Fig.\,\ref{fig:FBD}. From Eqns. (\ref{eqn:moment_x})\,--\,(\ref{eqn:force_z}), we can express the measured forces as
\begin{equation}
    \begin{bmatrix}
    F_x \\ F_y \\ F_z
    \end{bmatrix}
    = A
    \begin{bmatrix}
    v_1 \\ \vdots \\ v_8
    \end{bmatrix}\text{,}
\label{eqn:force_voltage}
\end{equation}
where $A\in	\mathbb{R}^{8\times3}$ is a sensitivity matrix that maps sensor voltage outputs to forces.

\subsection{Calibration Method}
The  actual value of $A$ can be estimated through linear least-squares fitting to calibration data. To improve the quality of calibration, we also fit a constant offset term for each direction of force. This results in the estimated sensitivity matrix $A^+$ having a dimension of $3\times9$.

To perform the calibration, the sensor was mounted on a Nano17 force sensor (ATI, Apex, NC, USA). The tip of the jaw attachment was affixed to a 3-axis linear stage as shown in Fig.\,\ref{fig:static_cal} and loaded in each Cartesian axis in increments of $0.5\pm0.1\,\,\text{N}$ through the target sensing range of 0 to $\pm3$\,N in the x- and y-directions and 0 to $\pm5$\,N in the z-direction. To reduce calibration errors due to possible hysteretic behavior, data was collected during loading and unloading. The quality of the calibration was evaluated using the root mean square error (RMSE), the normalized root mean square deviation (NRMSD) which is the RMSE normalized by the measurement range of the sensor, the coefficient of determination ($R^2$), and the hysteresis, which is the maximum difference between corresponding measured forces during loading and unloading normalized by the maximum force \cite{figliola2020hysteresis}.

\begin{figure}
    \centering
    \vspace{1.8mm}
    \includegraphics[width=0.6\linewidth]{fig/static_cal.pdf}
    \caption{Setup used for static calibration of a single force sensor.}
    \label{fig:static_cal}
\end{figure}

\begin{figure}[!t]
    \centering
    \includegraphics[width=\linewidth]{fig/gripper_angles.pdf}
    \caption[Gripper joint angle definitions.]{(a) Definitions of angles at the gripper of the sensorized RMIS tool. (b) Angle definitions for the software reported jaw angle $\theta^{'}_{\text{jaw}}$ and the actual jaw angle $\theta_{\text{jaw}}$.}
    \label{fig:gripper_angles}
\end{figure}

\begin{figure*}[!t]
\centering
\vspace{1.8mm}
    \includegraphics[width=\linewidth]{fig/exp_all.pdf}
    \caption{Setups for evaluating performance of force sensors (shown here with strain relief brackets) when mounted on the da Vinci Large Needle Driver tool. (a) The flat tissue manipulation (with inset showing from left to right, retraction, palpation and scraping movements), (b) the cylindrical stem manipulation, and (c) grip force measurement setups.}
    \label{fig:dual_jaw_exp}
\end{figure*}

\subsection{Performance Evaluation}

\subsubsection{Single Jaw Evaluation}
\label{ch4:singlejawevalmethod}
The single jaw evaluation was performed by exerting varying loads on the tip of the jaw attachment in all three Cartesian directions while the sensor was mounted to a Nano17 force sensor using the calibration setup shown in Fig.\,\ref{fig:static_cal} without the tip translation fixture. The sensor accuracy was determined by computing the RMSE and NRMSD of the force measurements and averaging them over three trials. The maximum error over three trials was also calculated.

\subsubsection{Dual Jaw Evaluation}
\label{dualjawevalmethod}
To measure the manipulation forces at the end-effector, each jaw of the da Vinci large needle driver needs to be instrumented with a force sensor. The forces measured in the reference frame of each sensor are first resolved into the reference frame of a da Vinci Research Kit (dVRK) \cite{kazanzides2014dvrk}. This is done using the robot forward kinematic model and the joint position estimates from the motor encoders (6 joints and the gripper angle) to obtain the individual jaw poses in the robot reference frame. The resultant force, $F_r$, is thus
\begin{multline}
    \prescript{0}{}{F_r} =  \prescript{0}{6}{T}(\theta_1,\dots,\theta_6,\theta_G) \prescript{6}{R}{T}(\theta_R) \prescript{R}{}{F} \\ + \prescript{0}{7}{T}(\theta_1,\dots,\theta_5,\theta_7,\theta_G) \prescript{7}{L}{T}(\theta_L) \prescript{L}{}{F}
\label{eqn:resultantforce}
\end{multline}
where $\prescript{j}{i}{T}$ are transformation matrices describing transformations mapping frame $i$ to $j$, with frame 0 being the dVRK origin, frames 6 and 7 describing the left and right tool gripper jaws respectively, and the L and R frames describing the local coordinate frames of each force sensor. The values of $\theta_{0}$ to $\theta_7$ are joint rotation angles, $\theta_R$ and $\theta_L$ are fixed rotation angles, and $\theta_G$ is the angle between the x-axis of frame 5 and the bisector of $\theta_6$ and $\theta_7$ (Fig.\,\ref{fig:gripper_angles}a).

Due to backlash and stretching of the tool actuation tendons, the computed jaw angle $\theta^{'}_\text{jaw}$ during grasping, as derived from the estimates of $\theta_6$ and $\theta_7$ and values of $\theta_R$ and $\theta_L$, is smaller than the actual sensorized tool jaw angle $\theta_\text{jaw}$ (Fig.\,\ref{fig:gripper_angles}b). To ensure that the correct joint angles are used during the pose computation, we define $\theta_\text{jaw}$ such that
\begin{equation}
    \theta_\text{jaw} = (\theta_6 - \theta_R) + (\theta_7 - \theta_L) =
\begin{cases}
    \theta^{'}_\text{jaw} ,& \text{if } \theta^{'}_\text{jaw}>\theta_\text{min}\\
    \theta_\text{min},              & \text{otherwise}
\end{cases}
\label{eqn:gripper jaw}
\end{equation}
where $\theta_6=\theta_7$, and $\theta_\text{min}$ is the minimum jaw angle during grasp.

Because each jaw is instrumented with a force sensor, the grasp force between the two jaws can be obtained. The grasp force was computed by using a two-point grasp model and resolving the forces measured at each sensor into the line of action between the two grasp points. Applying the rules derived by Yoshikawa and Nagai \cite{Yoshikawa1991}, the grasp force for a two-point grasp is \hlight{
\begin{equation}
    F_g = \text{min}(\, |(\prescript{G}{R}{T}\prescript{R}{}{F})\cdot\hat{j}| \, , \, |(\prescript{G}{L}{T}\prescript{L}{}{F})\cdot\hat{j}|)\text{,}
    \label{eqn:grasp}
\end{equation}
where $G$ is denotes the gripper frame of reference as shown in Fig.\,\ref{fig:gripper_angles}a and $\hat{j}=\begin{bmatrix}0&1&0\end{bmatrix}^\intercal$.}

To evaluate the sensor on realistic tissue manipulation tasks, we designed two environments that enable different types of manipulation forces to be exerted by an instrumented RMIS tool mounted on a teleoperated dVRK. The first environment consisted of an artificial silicone tissue (Limbs and Things, Savannah, GA, USA) placed over a sponge with a Nano17 force sensor placed underneath (Fig.\,\ref{fig:dual_jaw_exp}a). In this environment the tool can be teleoperated to perform palpation, scraping, and tissue retraction. However, due to the low friction of the silicone as well as the need to limit grasp forces in software to protect the sensor from damage, the tissue retraction force achievable in this environment were low compared to the sensor's operating range. The second environment consisted of a cylindrical silicone stem mounted on top of a Nano17 force sensor (Fig.\,\ref{fig:dual_jaw_exp}b). This setup thus allowed the teleoperator to exert higher retraction forces on the environment. In these two setups, the ground truth force during the teleoperated interactions can be obtained and used to evaluate the RMSE, RMSD, and maximum error of the resultant force measured from the dual jaw sensors over three trials of each task. However, the grip forces cannot be evaluated.

\begin{table*}[!t]
\centering
\caption{Summary of static calibration results for both sensors.}
\begin{tabular}{ccccccccccccc}
\toprule
\multirow{2}{*}{Sensor} & \multicolumn{3}{c}{RMSE (N)} & \multicolumn{3}{c}{NRMSD (\%)} & \multicolumn{3}{c}{$R^2$} & \multicolumn{3}{c}{Hysteresis (\%)}  \\
                        & x     & y     & z            & x     & y     & z              & x     & y     & z      & x+, x-    & x+, x-    & z+, z-       \\
\cmidrule(r){1-1}\cmidrule(lr){2-4}\cmidrule(lr){5-7}\cmidrule(lr){8-10}\cmidrule(r){11-13}
A                       & 0.023 & 0.056 & 0.044        & 0.388 & 0.876 & 0.438          & 0.999 & 0.996 & 0.999  & 3.96, 3.48 & 3.25, 2.85 & 2.78, 2.37    \\
B                       & 0.032 & 0.048 & 0.045        & 0.531 & 0.814 & 0.445          & 0.999 & 0.997 & 0.999  & 2.13, 2.73 & 2.15, 1.62 & 3.17, 3.36    \\
\bottomrule
\end{tabular}
\label{tbl:static_cal}
\end{table*}

To evaluate grip force, a separate experiment was devised. This involved attaching 3D-printed cantilevers to each interface of the Nano17 force sensor, leaving the sensor ungrounded as shown in Fig.\,\ref{fig:dual_jaw_exp}c. The sensorized RMIS tool was then teleoperated to grasp and release the opposing cantilevers five times each in three different trials. The sensor accuracy was determined by computing the RMSE and normalized root mean square deviation of the force measurements and averaging them over three trials. The maximum error over three trials was also calculated.

\section{Results and Discussion}

\subsection{Static Calibration}

\begin{figure}[!t]
    \centering
    \includegraphics[width=\linewidth]{fig/sensor_response_35.pdf}
    \caption{Load cell (LC) response under loading in each Cartesian direction. Load forces were measured using an ATI Nano17 force sensor within the calibration set up shown in Fig.\,\ref{fig:static_cal}.}
    \label{fig:sensor_response}
\end{figure}

\begin{figure}[!t]
    \centering
    \includegraphics[width=\linewidth]{fig/calibration_performance_35.pdf}
    \caption{Forces measured by the force sensors in each axis when loaded and unloaded independently in each Cartesian direction.}
    \label{fig:calibration_response}
\end{figure}

\begin{table*}[!t]
\centering
\vspace{1.8mm}
\caption{Summary of single jaw evaluation results for both sensors.}
\begin{tabular}{cccccccccc}
\toprule
\multirow{2}{*}{Sensor} & \multicolumn{3}{c}{RMSE (N)} & \multicolumn{3}{c}{NRMSD (\%)} & \multicolumn{3}{c}{Max Error (N)}  \\
                        & x     & y     & z            & x     & y     & z              & x     & y     & z                  \\
\cmidrule(lr){1-1}\cmidrule(lr){2-4}\cmidrule(lr){5-7}\cmidrule(lr){8-10}
A                       & 0.111$\pm$0.016& 0.105$\pm$0.015& 0.064$\pm$0.004       & 1.845$\pm$0.278 & 1.745$\pm$0.246 & 0.635$\pm$0.036          & 0.483 & 0.415 & 0.325              \\
B                       & 0.117$\pm$019& 0.146$\pm$0.013& 0.126$\pm$0.012        & 1.945$\pm$0.316 & 2.43$\pm$0.220  & 1.264$\pm$0.123          & 0.573 & 0.654 & 0.536              \\
\bottomrule
\end{tabular}
\label{table:single_jaw}
\end{table*}

\hlight{As described in Sec. \ref{dualjawevalmethod}, sensing of manipulation and grip forces at the end-effector required two sensors to be fabricated. Thus, the calibration was performed on two sensors, A and B, each corresponding to one jaw. The load cell responses during loading in each Cartesian direction during calibration are shown in Fig.\,\ref{fig:sensor_response} and indicate that when the sensing principle described in Sec. \ref{sensingprinciple} predicted a response  (dashed lines) from a given load cell, there was an appropriate response from that corresponding load cell. Additionally, there were some unexpected responses when no response was predicted due to sensor crosstalk and uneven plate contact that arose from slight errors in manufacturing. Overall, the sensors displayed good linearity over its functional range (Fig.\,\ref{fig:calibration_response}), with the redundant sensing architecture of the sensor mitigating any detrimental effects of crosstalk on sensor performance. The low deviation from unity in both loading and unloading seen in Fig.\,\ref{fig:calibration_response} also indicates that the sensor has low hysteresis despite its non-monolithic design. The results of the calibration procedure are summarized Table \ref{tbl:static_cal}. }

\subsection{Single Jaw Evaluation}

The results of the single jaw evaluation are summarized in Table \ref{table:single_jaw}. The low RMSEs (up to 0.146\,N for sensor B in the y-direction) showed that the sensor performance seen in calibration translated to the good performance in the dynamic loading scenario of the single jaw evaluation. As seen in Fig\,\ref{fig:single_jaw_eval}a, from 25\,s to 40\,s, the sensor could accurately track fast changes in applied force while maintaining the desired accuracy of below 0.375\,N RMSE. The sensor showed high maximum error (up to 0.654\,N for sensor B in the y direction) when loaded in multiple directions with one direction loaded close to max of the calibration range. Thus, the measured versus reference force plot in Fig.\,\ref{fig:single_jaw_eval}b, show some deviation from unity for all three Cartesian directions.

\begin{figure}[!t]
    \centering
    \includegraphics[width=0.9\linewidth]{fig/single_jaw_eval_35B.pdf}
    \caption{Selected results from single jaw evaluation of sensor B. (a) Recorded forces over time. (b) Measured force from the force sensor versus the reference force in each of the three Cartesian directions.}
    \label{fig:single_jaw_eval}
\end{figure}

\subsection{Dual Jaw Evaluation}

For both of the dual jaw evaluation tasks, we used a minimum jaw angle of $\theta_\text{min}=8.4^\circ$. This was required because the gripper does not fully close during grasping and thus the dVRK would report an incorrect gripper pose. This minimum jaw angle was empirically determined to reduce the error in the x-direction force measurements of both tasks with the RMSE of the sensor being below the minimum target of 0.375\,N as summarized in Table \ref{tbl:dual_jaw}. The force measurements from the sensor with respect to the ground truth for selected trials of the flat tissue and cylindrical stem manipulation tasks are shown in Fig.\,\ref{fig:dual_jaw_eval}.\hlight{ The plots show good tracking performance for lateral forces (x and y-direction in Fig.\,\ref{fig:dual_jaw_eval}a and b), as well as for palpation (z-direction in Fig.\,\ref{fig:dual_jaw_eval}a), and tension (z-direction in Fig.\,\ref{fig:dual_jaw_eval}b z-direction). }

\hlight{In our evaluation we identified three main sources of error. First, is the pose uncertainty of the jaws during grasping due to the cable-driven design of the dVRK robot where encoders are not placed directly on each joint of the surgical tool. Second, is the varying point of force application on the sensor, in which its calibration was done with forces applied only to the tip of the jaw attachment (Fig.\,\ref{fig:FBD} and Fig.\,\ref{fig:static_cal}). Third, is the small misalignment between the robot base coordinate frame and the reference force sensor coordinate frame.}

\begin{table*}[!t]
\centering
\vspace{1.8mm}
\caption{Summary of dual jaw evaluation results for the flat tissue and cylindrical stem manipulation tasks.}
\begin{tabular}{cccccccccc}
\toprule
\multirow{2}{*}{Task} & \multicolumn{3}{c}{RMSE (N)} & \multicolumn{3}{c}{NRMSD (\%)} & \multicolumn{3}{c}{Max Error (N)}  \\
                      & x & y & z                    & x     & y     & z              & x     & y     & z                  \\
\cmidrule(lr){1-1}\cmidrule(lr){2-4}\cmidrule(r){5-7}\cmidrule(lr){8-10}
Flat tissue &0.142$\pm$0.020& 0.078$\pm$0.013& 0.097$\pm$0.008 & 2.367$\pm$0.338 & 1.300$\pm$0.225 & 0.980$\pm$0.082          & 0.725 & 0.452 & 0.458              \\%[5pt]
Cylindrical stem &0.089$\pm$0.020& 0.149$\pm$0.023& 0.139$\pm$0.012 & 1.485$\pm$0.327 & 2.481$\pm$0.389 & 1.392$\pm$0.117          & 0.357 & 0.484 & 0.458              \\
\bottomrule
\end{tabular}
\label{tbl:dual_jaw}
\end{table*}

\subsection{Grip Force Evaluation}
\label{gripforceeval}

For the grip force evaluation, the RMSEs and standard deviations over three trials were 0.156$\pm$0.017\,N, and the maximum errors were 0.287\,N. A sample force plot over time using the default minimum jaw angle is shown in Fig.\,\ref{fig:grip_eval}.

The dVRK requires the teleoperator to momentarily close the jaws to trigger teleoperation, this movement causes the tool jaws to snap together, resulting in an impulsive load on the sensors. To prevent damage to the sensors during this movement, we limited the maximum grip force in software. Because of this limit, the highest peak grip force achieved during the evaluation was 1.45\,N and the average peak grip force was 1.35\,N which is below the range of our sensor. \hlight{The limit also prevented us from evaluating the dual jaw performance of the sensor in both manipulation tasks up to the same range used in single jaw evaluation. The sample would slip from grasp before the higher forces were reached. With modification of the underlying dVRK teleoperation code, it would be possible for the sensors to be tested to its lateral force upper range of 3\,N.}

\begin{figure}[!t]
    \centering
    \includegraphics[width=0.9\linewidth]{fig/dual_jaw_plots_35_v2.pdf}
    \caption[Selected results for the dual jaw evaluation tasks.]{Selected results for the (a) flat tissue  and (b) cylindrical stem manipulation tasks.}
    \label{fig:dual_jaw_eval}
\end{figure}

\begin{figure}[!t]
    \centering
    \includegraphics[width=0.9\linewidth]{fig/grip_force_35.pdf}
    \caption[Selected grip force measurement result.]{Selected grip force measurement result. Estimated jaw torque is derived from the motor current measurements at the instrument base as reported through the dVRK API.}
    \label{fig:grip_eval}
\end{figure}

\section{Conclusion}

In this work, we presented a 3-DoF force sensor specifically designed to facilitate RMIS research. The sensor can be used with an existing RMIS tool, is manufacturable using low-volume manufacturing techniques, and has interchangeable jaws that allow for adaptation to different RMIS tasks.

The current design is tolerant to manufacturing variability, leading to robust performance despite their presence. In both single-jaw standalone evaluation and dual-jaw evaluation on an RMIS tool, the sensor met the target accuracy specification of less than 0.375\,N RMSE. In single jaw evaluation, this accuracy was verified within the target sensing range of $\pm3$\,N for the lateral (x and y) directions and $\pm5$ for the axial (z) direction of force. Future work will investigate approaches to increase its robustness to manufacturing variability and enhance the assembly rigidity while maintaining or enhancing the sensing range.

In the dual jaw evaluation, the chief contributor to error was the pose uncertainty of the tool wrist and jaws. The uncertainty was most pronounced during high force manipulations and in grasping. To enable consistent, maximally accurate force measurements in cable-driven RMIS platforms like the dVRK, pose measurement approaches such as those based on stereovision \cite{Li2020Super,Lu2020SuperDeep} or robot-state information \cite{Oneill2018Gripper,kong2018grip,stephens2019gripconditions} will need to be further developed to improve real-time accuracy. Even with the limitation in pose measurement accuracy, the dual jaw sensor meets the human perception-based performance specifications and is a promising tool for enabling RMIS research that requires force information for bimanual tasks. The sensor designs have been made available at \url{https://github.com/enhanced-telerobotics/RMIS_force_sensor/} to allow researchers to manufacture and modify the sensor for use in applications such as providing haptic feedback, performing robotic force control, and collecting bimanual RMIS manipulation datasets that inform data-driven computational methods.

\section{Acknowledgements}

The authors would like to thank Mark Cutkosky, Thomas Daunizeau, Michael Raitor, and Guan Rong Tan for their guidance on electromechanical design and sensor calibration. The surgical robot CAD model is courtesy of Koray Okan.

\addtolength{\textheight}{-6.2cm}

\bibliographystyle{IEEEtran}


\section{Introduction}\label{into}

The present-day value of the cosmic microwave background (CMB) temperature has been accurately measured by COBE-FIRAS, $T_0= 2.7255\pm 0.0006$ K \cite{Fixsen}. If the expansion of the Universe is adiabatic, photon number is conserved, and the CMB spectrum was originally a black-body, the CMB temperature will evolve as $T(z)=T_0(1+z)$. This is a robust prediction of standard cosmology, but it is violated in many contexts, from exotic astrophysical processes to particle physics and cosmological models beyond $\Lambda$CDM \cite{Avgoustidis,Euclid}. Deviations from this behaviour would imply new physics, motivating efforts towards measurements at non-zero redshifts.

There are two paths for such measurements: the thermal Sunyaev-Zel'dovich effect at low redshifts, and spectroscopy of suitable molecular or atomic species at higher redshifts. The first measurement (as opposed to an upper limit) was achieved in 2000 \cite{Srianand}, and their number has been steadily increasing. Recently, the redshift range of these measurements has been extended to $z\sim6.34$ \cite{Riechers}.

These measurements have been used to constrain particular models, but there is no detailed study of their constraining power, as compared to that of other background cosmology data sets at similar redshifts. We present the first such study, showing that CMB temperature measurements are as constraining as Type Ia supernova and Hubble parameter data, and have important synergies with them. We discuss three different cosmological models, each with a different motivation.

\section{Background cosmology data}\label{data}

CMB temperature measurements are effectively background cosmology data. Therefore, our benchmark for their constraining power is two other canonical background cosmology data sets. The first is the Pantheon compilation of Scolnic \textit{et al.} \cite{Scolnic,Riess}. This is a 1048 supernova data set, compressed into 6 correlated measurements of $E^{-1}(z)$ (where $E(z)=H(z)/H_0$ is the dimensionless Hubble parameter) in the redshift range $0.07<z<1.5$. The second is the compilation of 38 Hubble parameter measurements by Farooq \textit{et al.} \cite{Farooq}: this includes both data from cosmic chronometers and from baryon acoustic oscillations. Together, the two sets contain measurements up to redshift $z\sim2.36$, and when using the two in combination we will refer tho this as the cosmological data. The Hubble constant is always analytically marginalized \cite{Anagnostopoulos}.

The CMB temperature constraints are listed in Table \ref{tab1}. At low redshifts they come from the thermal Sunyaev-Zel'dovich (SZ) effect, specifically from 815 Planck clusters in 18 redshift bins \cite{Hurier} and 158 SPT clusters in 12 redshift bins \cite{Saro}. At high redshifts we have spectroscopic measurements at various wavelengths and using several molecular or atomic species. The recent work of \cite{Klimenko} updates some earlier analyses \cite{Srianand0,Noterdaeme1,Noterdaeme2,Krogager}, specifically including the contribution of collisional excitation in the diffuse interstellar medium to the excitation temperature of the tracer species. The highest redshift measurement is at $z\sim6.34$, but must of them are until $z\sim2.5$, so the redshift range of the cosmological and CMB data sets is comparable. In our analysis we do not need the present-day value of the CMB temperature: $T_0$ is analytically marginalized, exactly as is done for $H_0$.

We report constraints for three classes of models, from on a standard statistical likelihood analysis, based on Matlab and Python codes that have been custom-built for this work, but validated against each other and also by comparison with analogous results for the same models in the literature, when these exist.

\begin{table}
\begin{center}
\caption{Measurements of the CMB temperature, with one-sigma uncertainties. The top and bottom parts of the table correspond to SZ and spectroscopic measurements respectively.}
\label{tab1}
\begin{tabular}{| c | c | c |}
\hline
Redshift & $T(z)$ & Reference \\
\hline
$0.037$ & $2.888\pm0.041$ & \cite{Hurier} \\
$0.072$ & $2.931\pm0.020$ & \cite{Hurier} \\
$0.125$ & $3.059\pm0.034$ & \cite{Hurier} \\
$0.129$ & $3.01_{-0.11}^{+0.14}$ & \cite{Saro} \\
$0.171$ & $3.197\pm0.032$ & \cite{Hurier} \\
$0.220$ & $3.288\pm0.035$ & \cite{Hurier} \\
$0.265$ & $3.44_{-0.13}^{+0.16}$ & \cite{Saro} \\
$0.273$ & $3.416\pm0.040$ & \cite{Hurier} \\
$0.332$ & $3.562\pm0.052$ & \cite{Hurier} \\
$0.371$ & $3.53_{-0.14}^{+0.18}$ & \cite{Saro} \\
$0.377$ & $3.717\pm0.065$ & \cite{Hurier} \\
$0.416$ & $3.82_{-0.15}^{+0.19}$ & \cite{Saro} \\
$0.428$ & $3.971\pm0.073$ & \cite{Hurier} \\
$0.447$ & $4.09_{-0.19}^{+0.25}$ & \cite{Saro} \\
$0.471$ & $3.943\pm0.113$ & \cite{Hurier} \\
$0.499$ & $4.16_{-0.20}^{+0.27}$ & \cite{Saro} \\
$0.525$ & $4.380\pm0.120$ & \cite{Hurier} \\
$0.565$ & $4.075\pm0.157$ & \cite{Hurier} \\
$0.590$ & $4.62_{-0.26}^{+0.36}$ & \cite{Saro} \\
$0.619$ & $4.404\pm0.195$ & \cite{Hurier} \\
$0.628$ & $4.45_{-0.23}^{+0.31}$ & \cite{Saro} \\
$0.676$ & $4.779\pm0.279$ & \cite{Hurier} \\
$0.681$ & $4.72_{-0.27}^{+0.39}$ & \cite{Saro} \\
$0.718$ & $4.933\pm0.371$ & \cite{Hurier} \\
$0.742$ & $5.01_{-0.33}^{+0.49}$ & \cite{Saro} \\
$0.783$ & $4.515\pm0.621$ & \cite{Hurier} \\
$0.870$ & $5.356\pm0.617$ & \cite{Hurier} \\
$0.887$ & $4.97_{-0.19}^{+0.24}$ & \cite{Saro} \\
$0.972$ & $5.813\pm1.025$ & \cite{Hurier} \\
$1.022$ & $5.37_{-0.18}^{+0.22}$ & \cite{Saro} \\
\hline
$0.89$ & $5.08\pm0.10$ & \cite{Muller} \\
$1.73$ & $7.9^{+1.7}_{-1.4}$ & \cite{Klimenko} \\
$1.77$ & $6.6^{+1.2}_{-1.1}$ & \cite{Klimenko} \\
$1.78$ & $7.2\pm0.8$ & \cite{Cui} \\
$1.97$ & $7.9\pm1.0$ & \cite{Ge} \\
$2.04$ & $8.6^{+1.9}_{-1.4}$ & \cite{Klimenko}  \\
$2.34$ & $10\pm4$ & \cite{Srianand} \\
$2.42$ & $9.0^{+0.9}_{-0.7}$ & \cite{Klimenko} \\
$2.53$ & $9.8^{+0.7}_{-0.6}$ & \cite{Klimenko} \\
$2.63$ & $10.8^{+1.4}_{-3.3}$ & \cite{Klimenko} \\
$2.69$ & $10.4^{+0.8}_{-0.7}$ & \cite{Klimenko} \\
$3.02$ & $12.1^{+1.7}_{-3.2}$ & \cite{Molaro} \\
$3.09$ & $12.9^{+3.3}_{-4.5}$ & \cite{Klimenko} \\
$3.29$ & $15.2^{+1.0}_{-4.2}$ & \cite{Klimenko} \\
$6.34$ & $23.1^{+7.1}_{-6.7}$ & \cite{Riechers} \\
\hline
\end{tabular}
\end{center}
\end{table}

\section{The Jetzer {\it et al.} model}\label{model1}

This model \cite{Jetzer1,Jetzer2} is based on an earlier idea by Lima \cite{Lima}. It is a decaying dark energy model where photon creation is possible, leading (among others) to a non-standard temperature redshift relation.  We assume flat universes. For sufficiently low redshfts for the radiation to be negligible, the Friedmann equation can be written
\begin{equation}
E^2(z)=\frac{3\Omega_m-m}{3-m}(1+z)^3+\frac{3(1-\Omega_m)}{3-m}(1+z)^m\,.
\end{equation}
The dark energy decay is parametrized by $m$ (with $m=0$ recovering flat $\Lambda$CDM) and it can be thought of as an effective dark energy equation of state $w_{eff}=m/3-1$.

\begin{figure}
\begin{center}
\includegraphics[width=\columnwidth]{figure1a.png}
\includegraphics[width=\columnwidth]{figure1b.png}
\includegraphics[width=\columnwidth]{figure1c.png}
\end{center}
\caption{\label{fig1}Constraints on the Jetzer \textit{et al.} model with $w_r=1/3$. Blue, green, and red curves show constraints from supernova, Hubble parameter and temperature data, while black ones depict the combined constraints. The top panel shows the one, two and three sigma constraints on the $m$--$\Omega_m$ plane, while the others show the posterior marginalized likelihoods for the individual parameters.}
\end{figure}

The model's coupling to radiation must necessarily be weak, as otherwise it would be already ruled out \cite{Chluba}: such a model can only be a small perturbation on the standard $\Lambda$CDM behaviour. Under these assumptions the temperature-redshift relation can be written
\begin{equation}
\frac{T(z)}{T_0}=(1+z)^{3w_r}\left[\frac{(m-3\Omega_m)+m(\Omega_m-1)(1+z)^{m-3}}{(m-3)\Omega_m}\right]^{w_r}\,,
\end{equation}
where we also relax the assumption of a strict radiation fluid, allowing it to have a generic constant equation of state $w_r$, with the standard case being $w_r=1/3$.

An early analysis \cite{Jetzer1,Jetzer2}, which using coeval data analogous to the present one (but kept $T_0$ fixed) finds constraints on $m$ of the order of $\sigma_m\sim0.1$ for the standard value of $w_r$, while if $w_r$ is free to vary the constraint on $m$ is of order $\sigma_m\sim0.2$, with $\sigma_{wr}\sim0.03$. On the other hand, \cite{Noterdaeme2}, with additional temperature measurements, obtained $\sigma_m\sim0.075$. Later \cite{Hurier,Saro} both independently obtained $\sigma_m\sim0.05$, for the standard value of $w_r$. A recent analysis \cite{Lucca} of a broader class of models uses the full the Planck data (but not the Table \ref{tab1} data) and finds, for the specific model we study, $\sigma_m\sim0.07$. On the other hand, \cite{Riechers} obtains uncertainties of$\sigma_m\sim0.05$, with $\sigma_{wr}\sim0.01$, but assumes a matter density fixed to the Planck value.

We restrict ourselves to the data introduced in the previous section. Figure \ref{fig1} shows our constraints, for the standard equation of state for radiation, $w_r=1/3$. We separately show the constraints from supernovas, Hubble parameter and temperature data, to highlight how their combination breaks degeneracies. We see that the three data sets are mutually consistent, and the key role of the temperature data in constraining $m$. Using the cosmology data alone one finds $\Omega_m=0.32\pm0.04$ and a weak constraint $m=0.23\pm0.18$, while if using the temperature data alone the matter density is unconstrained while $m=0.06\pm0.05$. For the combined data sets, we find
\begin{eqnarray}
\Omega_m&=&0.30\pm0.02\\
m&=&0.07\pm0.05\,,
\end{eqnarray}
which compares favourably with previous constraints---recall that we analytically marginalize $T_0$ as well as $H_0$. The combined data set improves the constraint on the matter density by a factor of two.

The above does assume the standard equation of state for radiation, $w_r=1/3$, but this can be relaxed, in which case the one-sigma constraints from the full data are
\begin{eqnarray}
\Omega_m&=&0.32\pm0.03\\
m&=&0.23\pm0.13\\
w_r&=&0.35\pm0.02\,;
\end{eqnarray}
here there are additional degeneracies with weaken the constraint on $m$, while $w_r$ is tightly constrained. Overall, there are no statistically significant deviations from the $\Lambda$CDM behaviour.

\section{The Canuto {\it et al.} model}\label{model2}

The motivation for this model \cite{Canuto1,Canuto2} is the fact that, although the effects of scale invariance vanish in the presence of particles with non-zero rest masses, one may assume that on cosmological scales empty space should still be scale invariant. This assumption ultimately leads to a bimetric theory, with a time-dependent function $\lambda$ playing the role of a scale transformation factor relating the ordinary matter frame to a scale invariant frame. The first of these can be thought of as the atomic (or physical) frame, while the second is a gravitational frame, in which the ordinary Einstein equations would still hold \cite{Canuto2}.

We assume homogeneous and isotropic universes, and a vanishing curvature parameter---previous work shows that relaxing this assumption has no significant impact \cite{Covariant}. We assume a generic power-law function $ \lambda(t)=(t/t_0)^p=x^p$, where $t_0$ is the present age of the universe and we have defined a dimensionless age of the universe, $x$; note that $\lambda(t_0)=1$, so $\Lambda$CDM is recovered for $p=0$. With these assumptions the Friedmann equation is
\begin{equation}
\left(E(z,x)+\frac{p}{2x}\Omega_\lambda\right)^2=\Omega_m(1+z)^{3(1+w)}x^{-p(1+3w)}+\Omega_{\Lambda}x^{2p}\,,
\label{scalef2}
\end{equation}
where $\Omega_\lambda=2/t_0H_0$; note that there is a consistency condition $(1+p\Omega_\lambda/2)^2=\Omega_m+\Omega_{\Lambda}$.

For numerical convenience this can be re-written as
\begin{eqnarray}
E(z,x)&=&\frac{\Omega_\lambda}{2x}\left[-p+\sqrt{N(z,x)}\right]\\
\frac{\Omega_\lambda^2}{4}N(z,x)&=&\Omega_m(1+z)^{3(1+w)}x^{2-p(1+3w)}+\Omega_\Lambda x^{2(1+p)}\,,
\end{eqnarray}
together with
\begin{equation}
\frac{dx}{dz}=- \frac{x}{1+z}\times\frac{1}{\sqrt{N(z,x)}-p}\,,
\end{equation}
with the initial condition $x=1$ at $z=0$. Finally, the temperature-redshift relation is \begin{equation}
T(z)=T_0(1+z)x^{-p/2}\,.
\end{equation}

\begin{figure*}
\begin{center}
\includegraphics[width=\columnwidth]{figure2a.png}
\includegraphics[width=\columnwidth]{figure2b.png}
\end{center}
\caption{\label{fig2}Constraints on the Canuto \textit{et al.} model discussed in the text. Blue, red and black curves show constraints from cosmology, temperature data, and the combined constraints, respectively. The panels show the posterior marginalized likelihoods for $p$ and the matter density.}
\end{figure*}

The model can be studied in two separate contexts. The first assumes $\Omega_\Lambda=0$, making it an alternative to $\Lambda$CDM. However, previous work \cite{Covariant} shows that such a scenario is ruled out: for low-redshift cosmological data one could obtain a statistically reasonable fit for $\Omega_m=0.11\pm0.03$ and a matter equation of state $w_m=0.51\pm0.16$, but such exotic values would be in conflict with higher-redshift data \cite{Planck}. Therefore, in what follows we focus on the second context, where $\Omega_\Lambda\neq0$ and the model is a parametric extension of $\Lambda$CDM.

Figure \ref{fig2} shows our results. We have chosen a uniform prior on $\Omega_\lambda$ corresponding to a present age of the universe from 13.5 to 27 Gyr; the lower limit corresponds to the age of the oldest identified galaxy, GN-z11 \cite{Eleven}, further corroborated by estimates from galaxy clusters \cite{Valcin}. We also assume a standard equation of state for matter $w=0$, and restrict our analysis to $p\le0$, in which case we effectively have a decaying cosmological constant---except if $p=0$, which recovers $\Lambda$CDM. These choices do not significantly impact our results.

From the cosmology data we find the one-sigma constraint $\Omega_m=0.268\pm0.020$ and the two-sigma lower limit $p>-0.18$; the temperature data does not provide a significant constraint on $\Omega_m$ but yields the lower limit $p>-0.16$. For the full data set we find
\begin{eqnarray}
\Omega_m&=&0.272\pm0.018\\
p&>&-0.14\,.
\end{eqnarray}
The addition of the temperature data improves the matter density contraint by $10\%$ and the limit on the parameter $p$ by $24\%$. We also note that in this model there is a mild preference for an age of the universe that is larger than the canonical one, as can be seen in the right panel of Fig. \ref{fig2}.

\section{The fractional cosmology model}\label{model3}

Finally, we consider the recently proposed fractional cosmology model \cite{Fractional}. This is based on mathematical formalism of fractional calculus \cite{Tarasov}; for our purposes, it suffices to note that it is a parametric extension of the usual concepts of differentiation and integration.

The model is envisaged as an alternative to $\Lambda$CDM, in the sense that it is assumed to contain no cosmological constant. Nevertheless, we will in what follows include a cosmological constant in the relevant equations. The Friedmann equation now has the form
\begin{equation}
E^2+\frac{1-\alpha}{H_0t_0} \frac{E}{x}=\Omega_m(1+z)^{3}x^{\alpha-1}+\Omega_\Lambda\,,
\label{fraceq1}
\end{equation}
where $t_0$ and $x$ have the same definitions as the previous section, and $\alpha$ is the fractional calculus parameter: in standard calculus one has $\alpha=1$. Since $E(z=0)=1$, we have $1/(H_0t_0)=(\Omega_m+\Omega_\Lambda-1)/(1-\alpha)$, so can also write
\begin{equation}
E^2+(\Omega_m-1)\frac{E}{x}=\Omega_m(1+z)^3x^{\alpha-1}+\Omega_\Lambda\,. \label{fraceq2}
\end{equation}
From this we see that for $\Omega_m+\Omega_\Lambda<1$ one needs $\alpha>1$; conversely, $\alpha<1$ would require $\Omega_m+\Omega_\Lambda>1$, otherwise one would have a negative age of the universe.

As in the previous section, we can also write this as
\begin{equation}
E(z,x)=\frac{1-\Omega_m}{2x}+ \sqrt{\frac{(1-\Omega_m)^2}{4x^2}+\Omega_m(1+z)^3x^{\alpha-1}+\Omega_\Lambda}\,,
\end{equation}
together with
\begin{equation}
\frac{dx}{dz}=\frac{1-\Omega_m-\Omega_\Lambda}{1-\alpha}\times\frac{1}{(1+z)E(z,x)}\,,
\end{equation}
Finally the temperature redshift relation is
\begin{equation}
    T(z)=T_0(1+z)x^{(\alpha-1)/3}\,;
\end{equation}
these recover the standard $\Lambda$CDM behaviour for $\alpha=1$.

\begin{figure*}
\begin{center}
\includegraphics[width=\columnwidth]{figure3a.png}
\includegraphics[width=\columnwidth]{figure3b.png}
\end{center}
\caption{\label{fig3}Constraints on the fractional model discussed in the text. Blue and red curves show constraints from cosmology and temperature data, respectively. The panels show the posterior marginalized likelihoods for the matter density; the left panel shows the physical region $\alpha>1$, while the right panel shows the $\alpha<1$ region, which is unphysical (implying a negative age of the universe) unless $\Omega_m>1$. Notice that no model parameters simultaneously provide a good fit to the two data sets}
\end{figure*}

A recent analysis \cite{Fractional}, using the Pantheon supernova and cosmic chronometers data and assuming $\Omega_\Lambda=0$, finds $\alpha=2.8\pm0.2$ and $\Omega_m=0.23\pm0.04$ (implying an age of the universe of $t_0\sim33.6$ Gyr). That work does not use baryon acoustic oscillation data, and treats the Hubble constant as an additional parameter, In our analysis, we do include baryon acoustic oscillation data (which is part of the Hubble parameter compilation), as well as CMB temperature data, and analytically marginalize $H_0$ and $T_0$.

Our results, for $\Omega_\Lambda=0$, are in the left panel of Figure \ref{fig3}. We assume generous uniforms prior for the matter density, $\Omega_m\in]0,1[$, and the fractional calculus parameter, $\alpha\in]1,4]$. The cosmological data leads to $\Omega_m=0.36\pm0.02$ and a fractional calculus parameter as large as possible, but this is incompatible with the temperature data, which leads to the two-sigma lower limit $\Omega_m>0.76$. A large matter density is preferred because it ensures $x\sim1$, approximately preserving the temperature-redshift relation and allowing it to fit the existing data. The model is therefore ruled out. Since the best-fit parameters of the two data sets are mutually incompatible, it makes no statistical sense to combine them in a join likelihood.

The mismatch in best-fit values from cosmology and temperature data also applies to the $\alpha<1$ region. This is shown in the right panel of Figure \ref{fig3}, where we also allowed for $\Omega_m>1$ to enable positive ages of the universe. Nevertheless. constraints on the matter density are almost unchanged---$\Omega_m=0.38\pm0.02$ from cosmology data and $\Omega_m=0.9\pm0.1$ from temperature data---and even in this wider parameter space, the model would be ruled out.

Finally, consider the case with $\Omega_\Lambda\neq0$, which would make the model a parametric extension of $\Lambda$CDM. For simplicity we fix the age of the universe to the standard value \cite{Planck}. The Friedmann equation can then be written
\begin{equation}
E(z,x)=\frac{1-\alpha}{2x}+ \sqrt{\frac{(1-\alpha)^2}{4x^2}+\Omega_m(1+z)^3x^{\alpha-1}+\alpha-\Omega_m}\,,
\end{equation}
with
\begin{equation}
\frac{dx}{dz}=-\,\frac{1}{(1+z)E(z,x)}\,,
\end{equation}
and the same temperature-redshift relation. Evidently, $\alpha=1$ corresponds to the $\Lambda$CDM model. The two data sets are now mutually compatible, and the addition of the temperature data improves constraints on the matter density from $\Omega_m=0.39_{-0.08}^{+0.06}$ to $\Omega_m=0.32\pm0.04$, again a significant improvement, while the overall constraint on the fractional calculus parameter is $\alpha=1.048\pm0.043$.

\begin{figure*}
\begin{center}
\includegraphics[width=\columnwidth]{figure4a.png}
\includegraphics[width=\columnwidth]{figure4b.png}
\end{center}
\caption{\label{fig4}Constraints on the fractional cosmology model with $\Omega_\Lambda\neq0$. Blue and red curves show constraints from cosmology and temperature data, while black ones depict the combined constraints. The panels show the posterior marginalized likelihoods for $\alpha$ and the matter density.}
\end{figure*}

\section{Conclusions}\label{concl}

Our analysis shows that recent progress in measuring the temperature of the cosmic microwave background at various non-zero redshifts hasled to a data set which plays a significant role constraining fundamental physics and cosmology paradigms. In particular, we have shown that these measurements have a constraining power that is comparable to that of other background cosmology observables probing similar redshifts, specifically Type Ia supernovae and Hubble parameter measurements.

In next decade, the high-resolution ELT spectrograph \cite{HIRES}, ANDES, will significantly improve the sensitivity of these measurements in the deep matter era, while next-generation ground and possibly space CMB experiments \cite{CORE} will lead to comparable progress in lower-redshift Sunyaev-Zel'dovich effect measurements.Together, these will ensure that temperature measurements will remain a competitive probe of new physics.

More broadly, the fact that the observed temperature-redshift relation is consistent with the canonical prediction is further evidence for the fact that $\Lambda$CDM is a robust paradigm. While it is undoubtedly a phenomenological approximation to a yet unknown more fundamental model, any plausible alternative model must closely reproduce its low-redshift behaviour.

\section*{Acknowledgements}

This work was financed by FEDER---Fundo Europeu de Desenvolvimento Regional funds through the COMPETE 2020---Operational Programme for Competitiveness and Internationalisation (POCI), and by Portuguese funds through FCT - Funda\c c\~ao para a Ci\^encia e a Tecnologia in the framework of the project POCI-01-0145-FEDER-028987 and PTDC/FIS-AST/28987/2017. CJM also acknowledges FCT and POCH/FSE (EC) support through Investigador FCT Contract 2021.01214.CEECIND/CP1658/CT0001. The project that led to this work was started during AstroCamp 2020. AMV was partially supported by Ci\^encia Viva OCJF funds.

\bibliographystyle{model1-num-names}


\section{Introduction}
\IEEEPARstart{H}{andwriting} Recognition has been a challenging task in image processing and computational linguistics. This task is traditionally divided into two main categories: online handwriting recognition and offline handwriting recognition. Because of the linguistic aspect and low-resource status of Vietnamese, recognizing Vietnamese handwriting images is still an arduous task. As collecting handwriting images is problematic, there are not many large-scale and high-standard dataset for researching handwriting recognition in Vietnamese.

Previous works \cite{Shen20161,krishnan2016generating} attempted to synthetically conduct offline handwriting datasets by using handwriting fonts to render images. However, datasets constructed by these approaches is not guaranteed to be natural as (1) for the approach proposed in \cite{krishnan2016generating}, the stable form of computer fonts opposes to the various handwriting styles of human, which obviously makes font-rendered handwriting images have different distribution to realistic handwriting images; and (2) for the approach proposed in \cite{Shen20161}, the colors of strokes and background are not considered carefully, and re-arranging segmented characters to a defined line does not clearly express the unconstrained properties of realistic handwriting images. To this end, we analysed and conducted experiments and then propose a simple but effective method, called Transferring method, to construct an offline handwriting dataset from an online handwriting dataset. Using this method, we introduce the first offline handwriting dataset in Vietnamese in order to motivate research community to explore and conduct more experiments to ultimately create an effective framework to tackle offline handwriting recognition task in Vietnamese.
\section{Challenge from writing characteristics of Vietnamese language}
Like other Western languages, Vietnamese uses the modern Latin script for its writing system. However, there is a clear difference that Vietnamese relies heavily on diacritics, which makes Vietnamese handwriting more complicated in a distinct way. In the Vietnamese alphabet, there are seven letters that always have their attached diacritics: \emph{ă} (breve), \emph{â}, \emph{ê}, \emph{ô} (circumflex), \emph{ơ}, \emph{ư} (horn) and \emph{đ} (bar); and five additional diacritics used to designate tone: grave (as in \emph{à}), acute (as in \emph{á}), hook above (as in \emph{ả}), tilde (as in \emph{ã}), and dot below (as in \emph{ạ}) \cite{truong2015vietnamese}. When being scrawled, some diacritics may be mistakenly read or seen as others, especially for computer vision systems. Specifically, a hook above glyph may look like a horn glyph (\emph{ủ} versus \emph{ư}), an acute glyph may be carelessly written to be like a grave glyph (\emph{á} versus \emph{à}), and a dot below under a letter \emph{i} may make the character \emph{ị} look like an exclamation mark (\emph{ị} versus \emph{!}). The examples above, and many more mistakes caused by carelessly written diacritics, frequently happen in real life and be an obstacle for optical character recognition for Vietnamese handwriting. Those flaws are also what we have encountered in our research, and we will analyze how the models stumble through them in Section 5.
\section{Related works}
Marti et al. \cite{Marti2002TheIA} constructed the IAM dataset by creating forms containing computer-printed texts for volunteers to write down these texts in their own writing style. After collecting these filled forms, line segmentation algorithms and word segmentation algorithms were used to extract lines from handwritten texts and words from extracted handwritten lines. On the other hand, the RIMES dataset \cite{6065550} was constructed based on the content of mails. In both cases, a large number of annotators were called in to annotate the datasets, and the annotation processes were re-corrected many times for the quality assurance.

Motivated by previous works and the shortage of Vietnamese offline handwriting dataset, we propose a novel method, called Transferring method, to construct an offline handwriting dataset in Vietnamese by synthesizing fundamental attributes from different handwriting datasets. Our conducted experiments ensure the resulting dataset totally keep mainly essential properties of an offline handwriting image dataset, which bridge the gap between handwriting images captured in the real world and images constructed synthetically. Consequently, we introduce a novel offline handwriting image dataset in Vietnamese constructed using the Transferring method, named UIT-HWDB.
\section{UIT-HWDB Dataset}
\subsection{Transferring method}
We present in this section the Transferring method to synthetically construct a Vietnamese handwriting image dataset. The novel constructed dataset must have the fully essential properties of a realistic handwriting image dataset. To achieve this requirement, the Transferring method first constructs the handwritten characters into an image, then applies the color variation to make the image more natural. The way of considering these real-world factors is detailed as follows:

\begin{algorithm}
    \caption{Pseudocode for the Transferring method ($U(a, b)$ is the discrete uniform distribution in range $[a, b]$, $Beta(\alpha, \beta)$ is the beta distribution).}\label{alg:cap}
    \begin{algorithmic}
        \Require \\
            \begin{itemize}
                \item $coords = ({x_n, y_n})_{n \in \mathbb{N}}$: the sequence of coordinates of a pen stroke.
                \item $all\_strokes$: the sequence of coordinates sequences of all strokes in the image.
                \item $bg\_dis = ({\alpha_m, \beta_m})_{m \in \mathbb{N}}$: the sequence of $m$ pairs of parameters for approximated beta distributions of background color.
                \item $stroke\_dis = ({\alpha_m, \beta_m})_{m \in \mathbb{N}}$: the sequence of $m$ pairs of parameters for approximated beta distributions of stroke color.
            \end{itemize}

        \Ensure a handwriting image containing characters drawed by the input coordinates.
        \\
        \Procedure{render\_image}{$all\_strokes$}
            \State $image \gets 255$ for all pixels
            \For{$coords$ in $all\_strokes$}
                \For{$ith$ in $[1:length(coords)]$}
                    \State $image[coords_{ith-1}:coords_{ith}] \gets 0$
                \EndFor
            \EndFor
            \State \textbf{return} image
        \EndProcedure
        \\
        \Procedure{render\_color}{image}
            \State $stroke\_ps \gets []$ \Comment{Sequence of stroke color variables}
            \State $bg\_ps \gets []$ \Comment{Sequence of background color variables}
            \For{$ith$ in [1:m]}
                \State $\alpha, \beta \gets stroke\_dis[ith]$
                \State $stroke\_ps.append(X \sim Beta(\alpha, \beta))$
            \EndFor
            \For{$ith$ in [1:m]}
                \State $\alpha, \beta \gets bg\_dis[ith]$
                \State $bg\_ps.append(X \sim Beta(\alpha, \beta)$)
            \EndFor
            \State $ith \sim U(1, m)$ \Comment{Randomly use $ith$ random variable, regarding both stroke and background color}
            \\
            \For{each pixel $p$ in image}
                \If{$p$ is 0} \Comment{Is stroke pixel}
                    \State $p \gets stroke\_ps[ith] \times 255.0$
                \Else \Comment{Is background pixel}
                    \State $p \gets bg\_ps[ith] \times 255.0$
                \EndIf
            \EndFor
        \EndProcedure
        \\

        \State $image \gets RENDER\_IMAGE(all\_strokes)$
        \State $image \gets RENDER\_COLOR(image)$
    \end{algorithmic}
\end{algorithm}

\subsubsection{Human handwriting style}\hspace{\fill} \\
The handwriting style is one of the specific properties that distinguish handwritten text from computer-printed text. Moreover, the variation of handwriting style is the most important characteristic required for a handwriting dataset, both online and offline.

Previous works \cite{krishnan2016generating} automatically used handwriting fonts to construct offline handwriting images, which causes the gap between natural handwriting images and synthetic handwriting images. To bridge this gap, we inherit handwritten characters from an online handwriting dataset, the VNOnDB \cite{nguyen2018database}. This approach guarantees the vast variation in natural writing style of human.

\subsubsection{Color}\hspace{\fill} \\
Recent works \cite{le2019end,le2018recognizing,9335877} constructed handwriting images in Vietnamese by connecting coordinates of pen stroke in the VNOnDB dataset \cite{nguyen2018database}. They kept the ink color as black (0 for pixel value) and the background color as white (255 for pixel value). This approach is intuitively not natural. Ignoring the color factors also means that we unintentionally reduce the complexity of the image compared with an image captured in reality, which causes the gap between research and realistic applications. To address these downsides, we take into account the complexity of the color.

For more details, we analyze the colors from the IAM dataset \cite{Marti2002TheIA}. The IAM dataset was collected manually, containing the real-world variation of both stroke colors and background colors. Therefore, the distributions of color in this dataset are natural enough to bridge the gap between synthetic images and real-world captured images. In the IAM dataset, for each image subset we intend to sample the stroke color values and background color values to get the corresponding color sampling distributions. We find that the color distributions have a finite range from 0 to 255 and have variously flexible shapes and skewnesses. After visualizing the color sampling distributions and taking into account the aforementioned statement, we finally come up with approximating the beta distribution for stroke colors and background colors of each subset. We then calculate and have a list of $\alpha$ and $\beta$ shape parameters for the beta distribution of the stroke colors and background colors. We apply to the strokes and background of handwriting images with random colors that follow the stroke and background color beta distribution which were generated from randomly picked pairs of $\alpha$ and $\beta$ parameters from the list. In this way, we ensure the nature of colors for handwriting images created synthetically.

\subsubsection{Stroke-width}\hspace{\fill} \\
Although most synthetic handwriting datasets are constructed in various ways that resemble images captured in the real world, they still miss the stroke-width variation. But in our work, experiments show that this factor is not essential and can be slightly ignored to keep the Transferring method simple without reducing the complexity of the dataset. We will carefully analyze this statement in Section 5.

\subsubsection{Annotation}\hspace{\fill} \\
Ground truth texts for the dataset are inherited from the ground truth texts of the VNOnDB dataset \cite{nguyen2018database}. However, from VNOnDB, there are some incorrect labels for handwriting images. Therefore, we also re-correct their mislabeled ground truth texts during the image constructing process.

\subsection{UIT-HWDB dataset}
Using the Transferring method with the VNOnDB as the online handwriting dataset and the IAM as the source for extracting colors, we construct the first novel offline handwriting dataset, the UIT-HWDB dataset. In more detail, our dataset has two parts: UIT-HWDB-word (110,745 unconstrained handwritten-word images) and UIT-HWDB-line (7,273 unconstrained handwritten-line images).

\begin{table}[ht]
    \renewcommand{\arraystretch}{1.3}
    \caption{Statistical comparison between handwriting image datasets (Para. stands for Paragraph, * indicates dataset in English, ** indicates dataset in French, $^{+}$ indicates dataset in German, $^{-}$ indicates dataset in Latin, $^{++}$ indicates dataset in Vietnamese).}
    \label{tab:1a}
    \centering
    \begin{tabular}{lccc}
        \hline
        Dataset             & Words   & Lines  & Para. \\ \hline
        RIMES**               & 83,493  & 11,365 & 1,500 \\
        IAM*                & 115,320 & 13,353 & -     \\
        Washinton Database*         & 4,894   & 656    & -     \\
        Parzival Database$^{+}$   & 23,487  & 4,477  & -     \\
        Saint Gall Database$^{-}$ & 11,597  & 1,410  & -     \\
        UIT-HWDB (ours)$^{++}$     & 110,745 & 7,273  & - \\ \hline
    \end{tabular}
\end{table}

We additionally reconstruct the test sets for the UIT-HWDB dataset as the original test sets of the VNOnDB generally contain easy-to-read images, which can lead to over-confidence in the ability of handwriting recognition methods. Specifically, we manually select both easy-to-read and hard-to-read handwriting images for our test set to evaluate the ability of recognition methods on many levels. The dataset is available at \url{https://github.com/hieunghia-pat/UIT-HWDB-dataset} for research purposes.

\subsubsection{Human evaluation}\hspace{\fill} \\
Our Transferring method is a semi-automatic method, which is undoubtedly prone to potential errors while constructing images. To carefully evaluate the errors of our Transferring method, we randomly collect 1,200 images from the UIT-HWDB-word train set and 600 images from the UIT-HWDB-line train set, then invite four people to manually evaluate the images (300 word-level images and 150 line-level images for each person). For the evaluation metric, we use the Cohen's kappa. As our method keeps the handwriting style, background color and stroke color consistent, the main factor that can cause the rendered characters to be abnormal is the choice of stroke width (detailed in section 5). Therefore the guideline for the four people to evaluate the images is to observe whether the rendered characters are deformed because of the improperly applied stroke width (stroke-width error) or the handwriting process of the annotators (non-stroke-width error). Results show the annotators' agreement for non-stroke-width error with Cohen's kappa coefficient $\mathbb{\kappa} = 1$ and $0$ for stroke-width error on the three collected sets. Our Transferring method totally keeps the original shape of characters.
\section{Experiments}
\subsection{Baseline methods}
The early approach for offline handwriting recognition is the combination of Convolutional Neural Network (CNN) \cite{o2015introduction}, Recurrent Neural Network (RNN) \cite{sherstinsky2020fundamentals}, and Connectionist Temporal Classification (CTC) loss \cite{graves2006connectionist}, which we call CTC-based methods. \cite{michael2019evaluating} proposed another combination using attention module \cite{cho2014learning} in their RNN layers. To enhance the ability of attention mechanisms in handwriting recognition methods, \cite{kang2020pay} proposed transformer-based methods, which have inherited the recent success of the transformer architecture in Natural Language Processing (NLP). In this experiment, we use TransformerOCR \cite{feng2020scene} and NRTR \cite{8978180} as the transformer-based methods. For attention-based methods, we implement Attention-based Encoder Decoder (AED) \cite{le2018recognizing} and finally we use CRNN \cite{shi2016end} and a method proposed in \cite{8269951} (BiCRNN for short) as the CTC-based methods.

\subsection{Evaluation metric}
For the Optical Character Recognition task, Damerau–Levenshtein Distance is the widely-used metric for measuring the performance of recognition methods. In this work, we use both two instances of this distance: Character Error Rate (CER) and Word Error Rate (WER).

\subsection{Experiments on the stroke-width factor}

In this experiment, we implemented the baseline methods on the word-level and line-level images.

To carefully analyze the effect of stroke-width factor on the ability of current handwriting image recognition methods, we created two versions for the UIT-HWDB-line dataset and UIT-HWDB-word dataset, named line-v1 and line-v2 (word-v1 and word-v2, respectively). The first version of these two datasets contains handwriting images without any variation in stroke width, while those in the second version contain the variation in the width of handwritten strokes. line-v1 and line-v2 are used to investigate the ability of TransformerOCR (transformer-based method) and BiCRNN (CTC-based method), while word-v1 and word-v2 are used to analyze the ability of NRTR (transformer-based method) and CRNN (CTC-based method). Note that we padded blank pixels to images in all training schemes in order to let all images have the same size and keep their original scale.

To mimic the variation of stroke-width factor, we approximate the width of pen stroke following the assumption: stroke-width is thinner when the writers draw up than when they draw down. Therefore, based on the ratio between the difference of two horizontal coordinates $ |\Delta x|$ and the difference of two vertical coordinates $|\Delta y|$ of two adjacent points, we construct a function $w$ which approximates the stroke width of handwritten characters:
\begin{equation}
    w(\theta) = m \cdot d(\theta)
\end{equation}
\begin{equation}
    d(\theta) = \displaystyle\frac{1}{1+e^{\alpha\theta + \beta}}
\end{equation}
In this experiment, we set $ \alpha = -0.1$, $\beta = 1.13$, $\theta = \arctan{\displaystyle\frac{\Delta y}{\Delta x}}$ and $m$ is the maximum thickness value a stroke can have.

To construct images with constant stroke width, we set for $\theta$ a constant value which makes the function $d(\theta)$ close to 1.

In practice, we observed the proper range for getting $m$ is $[2, 5]$. For the UIT-HWDB dataset, we randomly selected a value for $m \sim U(2, 5)$ where $U(a,b)$ is the discrete uniform distribution in range $[a, b]$  for every image.

\begin{table}
    \renewcommand{\arraystretch}{1.3}
    \caption{Results of TransformerOCR and BiCRNN on two versions of UIT-HWDB-line dataset (the first two rows were trained on line-v1 and last two rows were trained on line-v2).}
    \label{tab:3}
    \centering
    \begin{tabular}{lcccc}
        \hline
                & \multicolumn{2}{c}{line-v1} & \multicolumn{2}{c}{line-v2} \\ \cline{2-5}
                & CER          & WER          & CER          & WER          \\ \hline
        TransformerOCR & \textbf{11.59}        & \textbf{21.38}        & \textbf{11.42}        & \textbf{21.67}        \\
        BiCRNN & 12.06        & 31.84        & 12.28        & 31.49        \\ \hline
        TransformerOCR & 16.88        & 30.61        & 12.18        & \textbf{22.33}        \\
        BiCRNN & \textbf{11.50}        & \textbf{29.78}        & \textbf{10.77}        & 27.54        \\ \hline
    \end{tabular}
\end{table}

The experimental results are shown in Table \ref{tab:3}. According to Table \ref{tab:3}, when we trained TransformerOCR and BiCRNN on line-v1 and then evaluated these two models on both line-v1 and line-v2 test sets, we recognized they had approximately the same results. However, when we trained these models on line-v2 and then evaluated them on the line-v1 test set and line-v2 test set, they performed well on the line-v2 but yielded a few drawbacks on the line-v1. In conclusion, these results showed that with their strong CNN structures, TransformerOCR and BiCRNN, when being trained on version 1, can generalize well the pattern of handwritten text in images, hence they concurrently yielded the same results when predicting on the test set of version 2.

\begin{table}[!ht]
    \renewcommand{\arraystretch}{1.3}
    \caption{Results of NRTR and CRNN on the two versions of the UIT-HWDB-word dataset (the first two rows were trained on word-v1 and last two rows were trained on word-v2).}
    \label{tab:4}
    \centering
    \begin{tabular}{lcccc}
        \hline
                & \multicolumn{2}{c}{word-v1}    & \multicolumn{2}{c}{word-v2}    \\ \cline{2-5}
                & CER           & WER            & CER           & WER            \\ \hline
        NRTR & \textbf{8.28} & 22.00 & \textbf{10.10}         & 25.48          \\
        CRNN & 9.64 & \textbf{20.27} & 10.74         & \textbf{23.26}          \\ \hline
        NRTR & 18.01         & 37.73          & \textbf{8.32} & 21.21 \\
        CRNN & \textbf{16.66}         & \textbf{31.93}          & 9.77 & \textbf{20.62} \\ \hline
    \end{tabular}
\end{table}

Nevertheless, coming from NRTR and CRNN, we obtained another behavior. As Table \ref{tab:4} indicates, these two models only performed well on the version of the dataset that they had been trained. These results show that CRNN and NRTR can catch the specific attributes of each version of the dataset; thus, they are quite sensitive when making a prediction on another version of the images.

Finishing these experiments, we conclude that with current deep learning methods, stroke-width does not participate in forming the complexity of offline handwriting recognition task. Therefore we can ignore this factor to keep our Transferring method simple without lacking the nature of a handwriting image dataset.

\subsection{Experiments on baseline methods}
\begin{table}[!ht]
    \renewcommand{\arraystretch}{1.3}
    \caption{Results of baseline methods on UIT-HWDB-word and UIT-HWDB-line test set.}
    \label{table:2}
    \centering
    \begin{tabular}{l|cc|cc}
        \cline{1-5}
                       & \multicolumn{2}{c|}{UIT-HWDB-word} & \multicolumn{2}{c}{UIT-HWDB-line} \\ \hline
        Model          & CER               & WER                & CER               & WER                \\ \hline
        CRNN           & 9.93            & 20.76            & 31.23  & 100.00  \\
        BiCRNN         & 8.07            & 18.08            & 11.76           & 30.94            \\
        NRTR           & 8.25            & 21.31            & 49.72  & 91.70   \\
        TransformerOCR & \textbf{5.29}            & \textbf{10.37}            & \textbf{11.42}           & \textbf{21.67}            \\
        AED        & 6.70            & 14.68            & 12.53           & 35.31            \\ \hline
    \end{tabular}
\end{table}

We evaluated the baseline methods on UIT-HWDB-word and UIT-HWDB-line containing word-level and line-level handwriting images, respectively. Our results are shown in Table \ref{table:2}.

As shown in Table \ref{table:2}, both NRTR and CRNN obtained bad results on the UIT-HWDB-line dataset. This is because CRNN and NRTR were originally proposed to predict cropped-word scene text images, hence their backbone (CNN structure) were designed to be simple in order to enhance their performance \cite{shi2016end,8978180}. In contrast, Puigcerver et al. \cite{8269951} and Feng et al. \cite{feng2020scene} designed the BiCRNN and the TransformerOCR with every deep CNN structures. Therefore, these two models outperform NRTR and CRNN on recognizing both line-level and word-level handwriting images.

\subsection{Results analysis of the baseline methods}
As the results in our experiments, TransformerOCR achieved the best results among other methods on both the word-level and line-level dataset. Therefore, we focus on analyzing the common mistakes that TransformerOCR suffered.

\begin{figure*}[ht!]
    \centering
    \includegraphics[width=0.75\textwidth]{figure/TransformerOCR_error_2.jpg}

    \caption{Wrong-number prediction because of the low frequency of numbers in training dataset.}
    \label{fig:error_transformer}
\end{figure*}

\begin{figure*}[ht!]
    \centering
    \includegraphics[width=0.75\textwidth]{figure/TransformerOCR_wrong.png}
    \caption{Examples for failed cases that the Transformer can not read scribbled characters.}
    \label{fig:wrong_transformer}
\end{figure*}

After observing the images in the test set and the respective prediction of the TransformerOCR, we conducted the two types of error of the TransformerOCR which are (1) the wrong prediction of number characters and (2) the misunderstanding of scrawled handwritten characters.

For the first mistake (Figure \ref{fig:error_transformer}), as the domain of the UIT-HWDB dataset is daily newspaper, the occurrence of numerical characters is lower than of alphabet characters,
which causes the model to not learn enough to make the correct decision when facing these numerical characters. Statistically, the frequency of numerical characters (including Roman numerical characters) in the UIT-HWDB-word train set is approximately 0.065 and in the UIT-HWDB-word test set is approximately 0.0066 (we calculated these frequencies on the word-level dataset because the line-level is the line-segmented version of paragraph-level, and the word-level is the word-segmented version of the line-level, which indicates the frequency of number characters in the word-level dataset is the frequency of them in the line-level and paragraph-level dataset). But as the low occurrence of number characters on the train set as well as the test set, this type of mistake is typically not the main factor to reduce the performance of the TransformerOCR method.

The last type of error is the misunderstanding of scribbled characters which is the main challenge of recognizing handwriting images as well as the main challenge of our dataset. As depicted in Figure \ref{fig:wrong_transformer}, TransformerOCR totally failed to achieve the acceptable prediction for images containing scrawled handwritten characters. For humans, we find these images are hard-to-read, but this does not mean reading these images is infeasible. The only way to read these words is to infer them
together with meaning of the lines. Therefore we suggest the insight to improve the TransformerOCR or deep learning methods in general for this challenge is finding a way to exploit the meaning of the sentence to make models have better inference.
\section{Conclusion and Future Works}
We have introduced the Transferring method that synthetically forms a novel handwriting image dataset, which is useful for low-resource languages including Vietnamese. In addition, we presented a new benchmark for evaluating handwriting image recognition methods in Vietnamese constructed using this Transferring method. Finally, we have analyzed the confusion of state-of-the-art methods when performing on our dataset. We have identified that the appearance of number characters and badly scrawled text in the already linguistically complicated dataset are the causes behind the shortcomings in the predictions of the models. In the future, we continue to conduct experiments with various augmentation methods to tackle the shortage of number characters in our dataset. Moreover, we will analyze and research the way to effectively combine the handwriting recognition methods with pre-trained language models (like BERTology \cite{rogers2020primer}) to make them infer better when making predictions for scribbled handwriting images.