Camera-based method for identification of the layout of a robotic workcell

applied sciences Letter Camera-Based Method for Identification of the Layout of a Robotic Workcell Daniel Huczala * , Petr Ošˇcádal , Tomáš Spurný, Aleš Vysocký , Michal Vocetka and Zdenko Bobovský VSB-TU of Ostrava, Faculty of Mechanical Engineering, Department of Robotics, 17. listopadu 2172/15, Poruba, 708 00 Ostrava, Czech Republic; petr[email protected] (P.O.); tomas.spurny[email protected] (T.S.); [email protected] (A.V.); [email protected] (M.V.); [email protected] (Z.B.) *Correspondence: [email protected] Received: 9 September 2020; Accepted: 20 October 2020; Published: 30 October 2020   Featured Application: A fast and low-cost process for automated identification of positions of workcell components, including robots. Suitable for rapid deployment of robotic applications without a need of previous simulations or CAD modeling. Abstract: In this paper, a new method for the calibration of robotic cell components is presented and demonstrated by identification of an industrial robotic manipulator’s base and end-effector frames in a workplace. It is based on a mathematical approach using a Jacobian matrix. In addition, using the presented method, identification of other kinematic parameters of a robot is possible. The Universal Robot UR3 was later chosen to prove the working principle in both simulations and experiment, with a simple repeatable low-cost solution for such a task—image analysis to detect tag markers. The results showing the accuracy of the system are included and discussed. Keywords: robotic manipulator; identification; calibration; Jacobian; workcell layout detection 1. Introduction For robotic arms there has always been a trade off between the repeatability and absolute accuracy of the measurement of a robot’s positioning in 3D space, as examined by Abderrahim [ 1 ] or by Young [ 2 ]. Many manufacturers of the industrial robot present only the repeatability parameter in their datasheets, when it is way more precise than the absolute positioning. The general problem of robot accuracy is described with experiments by Salmani [3]. The absolute positioning of a robot examines how accurately the robot can move to a position with respect to a frame. To achieve better results, parameter identification and robot calibration are performed. Identification is the process in which a real robot’s kinematic (and possibly dynamic) characteristics are compared with its mathematical model. It includes determination of the error values that are afterwards applied into the control system, which improves the robot’s total pose accuracy using a software solution without the need for adjusting the hardware of the robot. A generally suggested method for robot calibration is the use of a laser tracker. The methodology identifies the error parameters of a robot’s kinematic structure, as is described by Nubiola [ 4 ]. The precision may be even increased, as Wu showed in [ 5 ] when trying to filter errors in measurements and finding optimal measurement configurations. In [6] Nguyen added neural network to compensate for non-geometric errors after the parameter identification was performed. Unfortunately, these solutions are very expensive because of the price of a laser tracker. One may rent a laser tracker if needed, but this is also a time consuming process due to the need to perform precise experiments, measurements and evaluations after every error made during the process that Appl. Sci. 2020,10, 7679; doi:10.3390/app10217679 www.mdpi.com/journal/applsci Appl. Sci. 2020,10, 7679 2 of 14 may lead to incorrect final results. Therefore, the wide deployment of laser trackers is ineffective for many manufacturers. There are other methods of robot calibration that tend to avoid the use of laser tracker. In [ 7 ], Joubair proposed a method using planes of a very preciously made granite cube, but acquisition of such a cube is not easy in general. Filion [ 8 ] or Moller [ 9 ] used additional equipment; in their case it is portable photogrammetry system. In [ 10 ] a new methodology is introduced by Lembono, who suggested to use three flat planes in a robot’s workplace with a 2D laser range finder that intersects the planes, but the simulation was not verified by an experiment. A very different approach was taken by Marie [ 11 ] where the elasto-geometrical calibration method based on finite element theory and fuzzy logic modeling was presented. On the other hand, the very precise results that the methodologies above wanted to achieve are not always necessary, and some nuances in a robot’s kinematic structure that appear during its manufacturing process are acceptable for the users of the robot. The problem they may face is determination of the workplace coordinate system (base frame) in which the robot is deployed and eventually the offset of the tool’s center point when a tool is attached to the robot’s mounting flange, when they need to position it absolutely in a world frame. For such applications the typical way to calibrate more robots is to use point markers attached to every robot, as described by Gan [ 12 ]. However, one important condition is that the robots need to be close together so they can approach each other with the point markers and perform the calibration. Additionally, there are a few optical methods using a camera to improve a robot’s accuracy. Arai [13] used two cameras placed in specific positions to track an LED tag that was mounted on a robot; the method we propose allows us to put the camera in any place, in any orientation that will provide good visibility. In [ 14 ] Motta or Watanabe in [ 15 ] attached a camera to a robot and performed the identification process, but this cannot be used for other robots or to track positions of other components at the workplace at the same time. Van Albada describes in [ 16 ] the process of identification for a single robot. Santolaria presented in [17] the use of on-board measurement sensors mounted on a robot. To avoid these restrictions, we propose a solution based on the OpenCV libraries [ 18 ] for Aruco tag detection by a camera, which adds to the calibration process benefits of simplicity, repeatability and low price. The outcomes may be used in offline robot programming, in reconfigurable robotic workplaces and for tracking of components, with as many tag markers and robots as needed, if the visibility for a camera or multiple cameras is provided. There are methods for 2D camera calibration already presented, and they can be divided into two main approaches. The eye-on-hand calibration, wherein the camera is mounted on the robot and a calibration plate is static, and the eye-on-base method with the calibration marker mounted on the robot with static cameras around [ 19 ]. There are also Robot Operating System (ROS) packages [ 20 , 21 ] providing tools for 2D or 3D camera calibration using these two methods. The ROS is an advanced universal platform that may be difficult for some researchers to be able to utilize. Our approach combines both eye-on-hand and eye-on-base calibration processes, avoids using ROS and can be applied not only to localize the base of a robot, but to also localize other devices or objects in the workplace that are either static or of known kinematic structure (multiple robots) in relation to chosen world frame. 2. Materials and Methods When an image with a tag is obtained, the OpenCV library’s algorithm inserts a coordinate system frame in the tag and can calculate transformation from the camera to the tag. If there are tags placed on all important components of a cell, such as manipulated objects or pallets, the transformation between them may be calculated as well. If there is an industrial robot deployed in a workplace, we can attach an end-effector with Aruco tags to it, perform a trajectory with transformation measurements and using mathematical identification methods calculate the precise position of its base, no matter where it is. Appl. Sci. 2020,10, 7679 3 of 14 2.1. Geometric Model of a Robot For such an identification, a geometric model that is as precise as possible of a robot needs to be determined. The Universal Robot UR3 was chosen for demonstrating the function of the proposed solution. Its geometric model used for all calculations is based on modified Denavit–Hartenberg notation (MDH), described in [ 22 ] by Craig. Our geometric model consists of 9 coordinate systems. The “b” frame is the reference coordinate system (world frame); later in our measurements it is represented by a tag marker placed on a rod. The “0” frame represents base frame of the robot. The frames from “1” to “6” represent the joints; the 6th frame position corresponds to the mounting flange. The “e” frame stands for the tool offset, in this case a measuring point that was focused by the sensor. The scheme of the model is illustrated in Figure 1. The MDH parameters are noted in Table 1; the oistands for offset of the ith joints in this position. Table 1. MDH parameters of the UR3 robot. iαi−1[rad] xi−1[mm] zi[mm] θi[rad] oi[rad] 1 0 0 151.90 q10 2π/2 0 119.85 q2π 3 0 243.65 0 q30 4 0 213.25 −9.45 q40 5−π/2 0 83.35 q50 6π/2 0 81.90 q6π Figure 1. UR3 with coordinate frames according to MDH. Appl. Sci. 2020,10, 7679 4 of 14 For a vector q= [q1 , q2 , q3 , q4 , q5 , q6]T representing the joint variables, a homogeneous transformation matrix Tbe(q) gives the position and orientation P of the UR3’s end-effector tool frame “e” with respect to the base frame “b” of the workplace. P=Tbe(q)(1) Tbe(q) = Ab0A01(q1)A12(q2)A23(q3)A34(q4)A45(q5)A56(q6)A6e(2) According to [ 22 ], matrix Ai−1,i in MDH notation is obtained by multiplying rotation matrix Rx along x axis, translation matrix Tx along x axis, rotation matrix Rz along z axis and translation matrix Tzalong zaxis. Ai−1,i=Rx(αi−1)Tx(xi−1)Rz(θi)Tz(di)(3) G=The geometric model of the UR3 is mathematically expressed by the transformation matrix Tbe(q) noted in Equation (2). Matrix Ab0 is displacement between the reference “b” frame and “0” frame; orientation difference is represented by R0 rotational matrix. Matrix A6e is displacement between the mounting flange and “e” frame of the end-effector. The objective of this study is to determine the 12 parameters of Ab0 matrix, so to find the base frame “0” of a robot in a workplace. To be able to achieve this, it is necessary to identify during the calculations also the displacement of the end-effector ( xe , ye , and ze ); however, the rotational part of the A6e can be freely chosen. The reason is that the transform is static (there is no variable between 6th frame of the robot and end-effector “e” frame). For simplicity, we choose Reas identity matrix. Ab0=      x0 R0y0 z0 0 0 0 1      ;A6e=      100xe 010ye 001ze 0 0 0 1      (4) R0=   r11 r12 r13 r21 r22 r23 r31 r32 r33   (5) In general, geometric models are idealized and very difficult to make comply with real conditions due to manufacturing inaccuracy and environmental conditions. Error values can be estimated and included into the mathematical models though. Finding the relations between theoretical and real models is the crucial task of robot calibration. To find such a relation, geometric parameters of a device have to be identified. Robot identification is a process wherein error parameter values are determined using results of a test measurement. In the following simulation the UR3 robot performed a trajectory as described in Section 3. Obtained data of the end-effector position Pc were compared with the robot’s position P(q)based on qifor each joint. The parameters x0 , y0 , z0 , r11 , r12 , r13 , r21 , r22 , r23 , r31 , r32 , r33 , xe , ye and ze were chosen to be identified. The reason for identification of the end-effector frame is because the Aruco tags may be placed on a low-cost end-effector by 3D printing and the designed CAD model transformation might be different than the real solution. On the other hand, we wanted to make the model as simple as possible, so we avoided MDH parameter identification between particular joints and links of the robot, which would lead to the robot’s calibration process. We assume that this simple identification process may compensate for the small errors between the links. When the transformation matrix Tbe(q) was defined earlier, the position vector of the end-effector P was represented by 4th column of Tbe with respect to the reference frame. If everything is ideal, we can consider this coordinate’s equal to the coordinate’s values calculated using the position sensor, as Equation (6) shows, for a specific q , where T0 be means 1st to 3rd elements of the 4th column of the transformation matrix. Pc=T0 be(q) = P(q)(6) Appl. Sci. 2020,10, 7679 5 of 14 2.2. Identification with the Jacobian Method The most common method for parameter identification is the application of a Jacobian, which is also described, for example, in [ 6 , 16 ]. This iterative method utilizes benefits of the Jacobian matrix that is obtained by partial derivative of the position vector (1st to 3rd elements of 4th column) of Tbe with respect to the parameters in X , the parameters that are going to be identified. Symbolically, the Jacobian is expressed as 3×15 Jimatrix in Equation (7), where pistands for a parameter of X. Ji=      ∂Tbex(q,X) ∂p1... ∂Tbex(q,X) ∂pn ∂Tbey(q,X) ∂p1... ∂Tbey(q,X) ∂pn ∂Tbez(q,X) ∂p1... ∂Tbez(q,X) ∂pn      (7) For every measured point i the Ji is determined. By applying all measured points a 3 n× 15 J matrix is obtained, where nis a number of measured points. J=      J1 J2 ... Jn      (8) As a first step of every iteration a position vector Ym is calculated using values Xj , where j represents the iteration step. For the first iteration, guessed values X0 are used. The q is the vector of joint variables measured by robot’s joint position sensors. Ym=Tbe(q,Xi)(9) The next step is to compare and calculate the difference between the position measured by camera Yc and the previously calculated position Ym , so ∆Y is determined. Yc is n× 3 matrix, where n stands for number of measured points, and there are three measured coordinates x,y,z. ∆Y=Yc−Ym(10) The key equation of this method is Equation (11). When a position changes, the Jacobian matrix changes too; therefore, ∆xcan be observed as the change of the parameters in X. ∆Y=J∆x(11) By using matrix operations in Equations (12) and (13), the values of ∆xare determined. (JTJ)−1JT∆Y= (JTJ)−1JTJ∆x(12) ∆x= (JTJ)−1JT∆Y(13) At the ends of iterations we added the computed values to the Xj+1 . A convergent check follows to decide whether another iteration step is needed. Xj+1=Xj+∆x(14) 3. Simulations and Experiment Two types of simulations and one experiment were performed to verify the proposed method of parameter identification. Simulation A was calculated only with absolute positions of the end-effector coordinates determined by CoppeliaSim software with the built-in UR3 model, as shown in Figure 2. Appl. Sci. 2020,10, 7679 6 of 14 The robot moved along the same path in both simulations and the experiment, consisting of 250 points. The robot stopped at each pose and the measurements were taken. The reason for choosing such a path was to obtain coordinates of the joints that were as different as possible; on the other hand, due to the experiment that was performed with cameras we needed to guarantee the visibility of the Aruco tags, which were used for the simulation B and the experiment. Figure 2. Simulated path of the UR3 in CoppeliaSim, simulation A. 3.1. Simulation A The robot was moved along a defined path with fixed points to stop at. Once it stopped, the joint coordinates and the position vector of the end-effector related to the world base frame were acquired. With these two sets of input values, the identification was made using the methods described in the previous chapter. Results may be seen in Section 3.5. 3.2. Simulation B For simulation B and the experiment, there were cameras and OpenCV libraries [ 18 ] applied for image processing to detect the Aruco tags. Based on the previous research by Oscadal [ 23 ], we used a 3D gridboard with tags, which improves the reliability and accuracy of detection in comparison with basic 2D tags. The gridboard represents a coordinate frame; in our case it was the base frame “b” and the end-effector frame “e” as shown in Figures 3and 4. The OpenCV library algorithm can calculate transformation from a camera to a tag. In real-time measurements, we used Equation (15) thanks to which the position of any camera was not important. Matrix Tcb is the transformation from the camera to the base; matrix Tce is the transformation from the camera to the end-effector. The “c” frame is the camera frame. On the other hand, the position of a camera may be saved for later operations. Tbe =T−1 cb Tce (15) In the simulation B (Figure 3) we deployed the image analysis in CoppeliaSim, with a single camera of resolution 1280 × 720 px. The virtual camera was self-calibrated and the detection parameters of the OpenCV library for finding the Aruco tags were set similarly as in [ 23 ]. The dimensions of the tags were 70 ×70 mm with a 6 ×6 bit matrix. Appl. Sci. 2020,10, 7679 7 of 14 Figure 3. Simulation B with tags and a camera. Figure 4. Experiment A setup; point of view of camera 2. 3.3. Experiment A Three independent cameras were placed around the UR3 robot to observe its trajectory and to calculate transformations in a laboratory during a real experiment. Intel RealSense D435i cameras with 1280 × 720 px resolution were used, and even though they are depth cameras, only the simple 2D RGB pictures were analyzed. The specifications of the cameras are shown in Table 2. The cameras were self-calibrated following the methodology used in [23]. As already said before, the robot went through 250 positions on the path. To reduce the inaccuracy of the detection in reality, 10 camera frames were taken for every position, which gave us 2500 measured points. In total for the three cameras, 7500 pictures were analyzed during the calibration process. Table 2. Specifications of Intel RealSense D435i camera [24]. Parameter Value Resolution 1920×1080 px Frame Rate 30 fps Sensor FOV (H ×V×D) 69.4◦×42.5◦×77◦(±3◦) Dimensions 90 ×25 ×25 mm Connection USB-C 3.1 Gen 1 3.4. Experiment B To observe the impact of placement of the world base frame “b,” another experiment was performed as shown in Figure 5. Please, notice the difference in position and rotation of the base frame. Appl. Sci. 2020,10, 7679 8 of 14 Only one camera (the same type, resolution and calibration) was used in this case with robot following a similar trajectory as in the previous experiment A. The robot went through 200 poses; at each pose, five images were taken and analyzed, so 1000 measurements in total were made. Figure 5. Experiment B setup. 3.5. Results The calibration results of simulations are presented in Table 3. The data calculated based on simulation A (without tag detection, only end-effector position tracking) show high-precision identification with a very small difference in comparison to the expected results (difference ∆ is shown in the brackets). Such an accuracy could be achieved using, for example, a laser tracker in reality. For simulation B, when the tags were detected using simulated camera, the error was higher (maximum 2.59 mm for x0 ), which gave us an idea about how accurate the system might be, so the noise from the environment was lowered, but the camera parameters were kept. In Table 4are the results of the experiment A. The best values were obtained when all the results of the three cameras were combined and analyzed together. The error for the base frame ( x0 , y0 , and z0 ) was maximally 7.61 mm in the z0 direction. We performed other measurements following the same strategy; they all provided similar results, which made it made clear that the position of a camera has an influence on the detection accuracy, which is supported by results of Krajnik’s research [25]. In Table 5are the results of the experiment B. It proves the possibility of placing the base frame freely with respect to the robot. The trajectories are compared in Figures 6–10. Measured cameras/measured path are plotted points that were obtained by camera/end-effector position sensor; robot path is a self-check plot of the end-effector path after the identified parameters of X were applied in the transformation matrix. Points of origin are points where the world base frame was determined for every point on the path. In Figures 11 and 12 we can observe the error values in box plots. Appl. Sci. 2020,10, 7679 9 of 14 Table 3. Results of simulations for the Xvector after identification. Parameter Expected Simulation A Simulation B x0[mm] 431.00 430.98 (∆0.02) 433.59 (∆2.59) y0[mm] 555.00 555.02 (∆0.02) 556.01 (∆1.01) z0[mm] −460.00 −460.07 (∆0.07) −458.41 (∆1.59) r11 [-] 1.00 1.00 (∆0.00) 1.006 (∆0.006) r21 [-] 0.00 0.00 (∆0.00) −0.006 (∆0.006) r31 [-] 0.00 0.00 (∆0.00) 0.001 (∆0.001) r12 [-] 0.00 0.00 (∆0.00) 0.003 (∆0.003) r22 [-] 1.00 1.00 (∆0.00) 1.004 (∆0.004) r32 [-] 0.00 0.00 (∆0.00) −0.000 (∆0.000) r13 [-] 0.00 0.00 (∆0.00) −0.006 (∆0.006) r23 [-] 0.00 0.00 (∆0.00) −0.001 (∆0.001) r33 [-] 1.00 1.00 (∆0.00) 1.000 (∆0.000) xe[mm] 0.00 −0.05 (∆0.05) 0.14 (∆0.14) ye[mm] 0.00 −0.02 (∆0.02) 0.23 (∆0.23) ze[mm] 66.00 65.98 (∆0.02) 64.78 (∆1.22) Table 4. Results of experiment A for the Xvector after identification. Parameter Expected Cam 1 Cam 2 Cam 3 Cameras Combined x0[mm] 431.00 432.53 430.81 426.91 430.05 (∆0.95) y0[mm] 555.00 551.15 561.06 555.45 555.96 (∆0.96) z0[mm] −460.00 −451.33 −449.27 −456.496 −452.39 (∆7.61) r11 [-] 1.00 1.006 1.024 1.013 1.014 (∆0.014) r21 [-] 0.00 −0.018 −0.039 −0.015 −0.024 (∆0.024) r31 [-] 0.00 −0.004 −0.013 0.008 −0.003 (∆0.003) r12 [-] 0.00 0.020 0.014 0.007 0.013 (∆0.013) r22 [-] 1.00 1.008 1.019 1.004 1.011 (∆0.011) r32 [-] 0.00 0.000 −0.009 0.010 0.000 (∆0.000) r13 [-] 0.00 0.010 0.010 0.001 0.007 (∆0.007) r23 [-] 0.00 0.012 −0.009 −0.006 −0.001 (∆0.001) r33 [-] 1.00 0.990 0.969 1.011 0.991 (∆0.009) xe[mm] 0.00 0.18 0.49 0.05 0.24 (∆0.24) ye[mm] 0.00 0.52 0.73 0.68 0.65 (∆0.65) ze[mm] 66.00 64.09 61.00 66.03 63.73 (∆2.27) Table 5. Results of experiment B for the Xvector after identification. Parameter Expected [mm] Experiment B [mm] x0[mm] 611.00 613.88 (∆2.88) y0[mm] 29.00 30.94 (∆1.94) z0[mm] −27.00 −25.12 (∆1.88) r11 [-] - 0.506 r21 [-] - −0.873 r31 [-] - 0.012 r12 [-] - 0.872 r22 [-] - 0.515 r32 [-] - 0.020 r13 [-] - −0.016 r23 [-] - −0.013 r33 [-] - 1.043 xe[mm] 0.00 −0.05 (∆0.05) ye[mm] 0.00 0.69 (∆0.69) ze[mm] 66.00 67.59 (∆1.59)