# Research on the Bolt Positioning System Based on Multieye Vision Industrial Robots

Authors: Jun Kong; Qiang Li; Yi Yin; Cen Li; Wei Luo
Journal: Shock and Vibration
 (2022)

---

## Abstract

In order to improve the accuracy of bolt positioning for industrial robots, this paper studies the bolt positioning system of industrial robots combined with multieye vision technology. This paper introduces three reconstruction algorithms ART, SIRT, and SART from the perspective of theory and implementation. Moreover, the projection matrix calculated with 0, 1 weighting, length weighting, and linear interpolation weighting is used for three reconstruction algorithms to carry out reconstruction experiments. In addition, this paper combines the actual working conditions of bolt positioning to construct the system of this paper and conducts system simulation research combined with the working conditions of the robot. The research shows that the bolt positioning system based on the multieye vision industrial robot proposed in this paper has a good performance in bolt positioning.


---

## Body

## 1. Introduction

Robots are mechanical devices that perform work automatically. It can either obey human commands, run preprogrammed programs, or act according to principles and programs formulated with artificial intelligence technology. Moreover, its task is to assist or replace human work in production, construction, and other works.Literature [1] developed a bolt-fastening system for automotive engine end covers. The system takes the robot as the core, adopts an intelligent robot, and is equipped with a corresponding tightening mechanism, jacking mechanism, and pressing mechanism to ensure that the production qualification rate reaches more than 99%. At the same time, it can be applied to the tightening of bolts at other types of end caps, with high flexibility. . A bolt replacement robot for replacing the isolation switch on the high-voltage line of the power station is designed. The robot is equipped with a vision sensor, a four-axis mobile platform, and an end effector [2]. Through visual positioning, the bolts that need to be replaced can be found on the isolation plate, and mechanical disassembly and replacement of new bolts can be performed at the end. This design replaces manual work on high-voltage wires using robot operation, improves efficiency, reduces the danger of manual high-altitude work, and improves the automatic maintenance level of the power station [3]. A wind turbine tower connecting bolt inspection robot is designed, which improves the working condition of manual tightening of the connecting bolts of the wind turbine tower. Its structure consists of three parts: a circumferential motion mechanism, an adaptive mechanism, and a three-point clamping mechanism. The mechanism consists of an electric push rod and a clamp, which drives a torque wrench for maintenance [4]. This design solves the maintenance problem of wind turbines and greatly reduces maintenance costs; a bolt-tightening robot for tightening bolts has been developed. The robot is a high-power robot equipped with a 170 F gasoline engine, an electromagnetic clutch, coupling, energy storage body, reversing mechanism, etc. [5], built a control system for the robot, and completed the actual test. This design improves the work efficiency of rail fastening in railway maintenance work. Many excellent results have been achieved in the research of bolt assembly robots, but there are few types of research on light-load vision-guided robots in production lines. At the same time, further research and development are needed for the flexible adaptive structure in the process of helical motion. Research design.Literature [6] proposes a stereo vision ranging system integrated into a humanoid robot. The system uses the image processing library OpenCV and OpenGL to design a corresponding automatic recognition software system according to the requirements of automatic feature recognition and detection of corresponding points and uses OpenGL. Open Graphics Library shows the 3D model obtained from the reconstruction. It was finally used in the actual test of the robot and got good results, which improved the work efficiency. Literature [7] developed and tested a vision-guided grasping system for Phalaenopsis tissue culture seedlings. The system applies a binocular stereo vision algorithm to calculate the 3D coordinates of the grasping point and uses an image processing algorithm that locates the grasping point to determine the appropriate grasping point on the root. Meanwhile, his research team developed and tested a device suitable for gripping Phalaenopsis tissue culture seedlings. Finally, the binocular vision localization algorithm is integrated with the robotic grasper to construct an automatic grasping system. The experimental results show that the automatic grasping system has a success rate of 78.2% in grasping the seedlings in the proper position. The welding seam tracking and feedback technology of the welding robot are researched, and computer vision is used to identify and find the position of the welding crease and then weld the position of the crease. The system is a binocular system based on two CCD cameras. The cameras are installed on opposite sides of the outer hollow shaft to capture images of the welding seam [8]. At the same time, the electromagnetic air valve and two cylinders are used to work with the welding device. The research solves the problem of weld positioning accuracy during mechanical welding and improves the level of welding automation [9].For the production line where the robot base and the RV reducer are connected, a binocular vision guidance scheme for assembly is designed. The scheme uses HALCON for binocular vision processing and camera calibration at the same time. Through median filtering, adaptive K-means segmentation of lab color space, template matching based on image pyramid, subpixel edge detection, and so on, the collected images are contour fitted, and the coordinates of feature points of threaded holes are obtained [10]. Finally, guide the robot to assemble; after completing the algorithm development, use VisualStudio to design control software and finally complete the research and development of the entire assembly system.Literature [11] designed and developed a part recognition and detection system under the binocular camera and used the system to combine with an industrial robot to complete the actual grasping measurement. In this system, the contour of the part edge is identified by the improved Canny algorithm, and the feature point extraction and stereo matching are carried out by using the scale-invariant feature conversion method; the mathematical model of the pose detection system is established by using the stereo vision 3D reconstruction method. The coordinates of the parts on the worktable are obtained, and the parts are grasped through the control of programmed software. This research plays a very important role in the field of automatic loading and unloading of industrial robots [12]. In-depth research on the visual guidance technology of the bolt tightening robot was carried out, and the design of binocular software was completed on the VisualStudio platform; the three parts of the binocular target, image correction, bolt feature point extraction, and pose measurement were completed [13]. In the process of feature extraction, image preprocessing, dynamic threshold segmentation, Minkowski addition and expansion, Minkowski subtraction and erosion, subpixel precision contour, and rammer edge fitting are adopted, and finally, the six corner coordinates of the bolt are extracted to guide the robot to grasp [14].In this paper, multi-eye vision technology is used to study the bolt positioning system of industrial robots to improve the bolt positioning effect of industrial robots.


## 2. Basic Principles of Binocular Vision Imaging

### 2.1. Physical Basis of Binocular Vision Imaging

Robotic vision light is a high-energy electromagnetic wave with a certain energy and penetrating ability, which can penetrate some substances (such as human tissue) that visible light cannot pass through, as shown in Figure1. Generally, visible light has a longer wavelength, and when a photon hits an object, part of it is reflected and most of it is absorbed, while the wavelength of the robot’s vision light is extremely short, and the photon contains high energy. The penetration of the robot’s visual light is related to information such as the equivalent atomic number and density of the irradiated material, and the transmittance of the robot’s visual light is stronger for the material with a lower atomic number and vice versa. The transmission of robot vision light is an important basis for binocular vision imaging.Figure 1
Schematic diagram of the material that transmits light to the robot vision.A beam of robotic vision rays injected into a homogeneous material is considered, as shown in Figure2(a). Considering Beer’s theorem, we have(1)I=I0exp−μΔx⇒μΔx=lnII0.Figure 2
Schematic diagram of visual light passing through. (a) Schematic diagram of robot vision light passing through homogeneous material. (b) Schematic diagram of the robot’s vision of light passing through inhomogeneous materials.
(a)(b)It can be seen from formula (1) that the object with a high μ value causes more attenuation of the robot’s visual photons than the object with a low μ value. For example, the μ of bone is higher than that of soft tissue, indicating that it is more difficult for robotic vision photons to penetrate bone than soft tissue. On the other hand, the μ value of air is almost 0, indicating that the input and output of X rays hardly change on the path through the air.When the material scanned by theX-ray is inhomogeneous, the medium distributed along path 1 can be discretized into several continuous small blocks. When these small pieces are small enough, the medium inside the small pieces can be considered to be homogeneous and have the same attenuation coefficient. The thickness of each discrete small block is assumed to be Δx, and the attenuation coefficients of each discrete small block are μ0,μ1,μ2,μ3,…,μn, respectively, as shown in Figure 2.The ray intensity of the robot vision light after passing through the first small block isI1, and the ray intensity after passing through the second block is I2….The final transmission intensity isIn, then there are(2)I1=I0exp−μ1Δx,(3)I2=I1exp−μ2Δx.Substituting formulas (2) into (3), we get(4)I2=I0exp−μ1+μ2Δx.It continues to accumulate attenuation values along the propagation direction of the robot’s vision light until the final transmission intensityIn of the robot’s vision light when it leaves the illuminated object, as shown in the following equation:(5)In=I0exp−μ1+μ2+⋯+μnΔx.We take the positive exponent of formula (5) and express it in summed form, we get(6)p=∑i=1nμiΔx=lnI0In.p in formula (6) is the projection. If the incident intensity I0 and outgoing intensity In of the X ray are known, a linear formula with μii=1,2,…,n as the last known is obtained according to formula (6). When Δx⟶0, formula (6) can represent the summation of continuous variation, and its integral form is(7)p=lnI0In=∫Lμl.In formula (7), μl is a continuous function of the decay rate with respect to path l. The process of finding the attenuation coefficient function μl by projection is called back projection. If the two-dimensional density function f (x, y) is used to describe the attenuation rate of the two-dimensional plane, then the problem of binocular vision imaging can be expressed as the measured linear integral of an object is given, and the decay rate at each point is calculated to produce two-dimensional density data.


### 2.2. Analytical Reconstruction Algorithm

Two-dimensional Radon transform is a projection transform of straight line integral, and it can be defined in many forms. This paper takes the most commonly used case as an example.We assume that in the plane areaR2, any point (x, y) can be represented by polar coordinates ρ,θ, ρ represents the distance from the line to the origin, θ represents the angle between the straight line l and the positive y-axis, and the function f (x, y) is the image to be reconstructed. We assume that any straight line l∈R2, then the two-dimensional Radon transform of the function f (x, y) is defined as(8)p=∫lfx,ydl=∫lf^r,θdl=∫−∞+∞f^L2+l2,ϕ+tan−1lLdl.The formula of the straight line can be expressed asρ=xcosβ+ysinβ in polar coordinates, then formula (8) can be further expressed as(9)p=∫−∞+∞fx,ydl=∫−∞+∞∫−∞+∞fx,yδxcosβ+ysinβ−ρdxdy.Among them, (x, y) represents the position of the reconstructed pixel in the Cartesian coordinate system, and δx represents the sampling function. This process is the integration of the image in a straight line. According to the previous physical principle of the X-ray, we can abstract the attenuation process of the X-ray into this integration process, and the projection value of the X-ray is the value of the Radon transform at angle β. If the projection is performed at multiple angles, the Radon value of each angle can be obtained, and the two-dimensional image f (x, y) of the original plane can be reconstructed by performing the Radon inverse transformation on these obtained Radon values.The formula for the two-dimensional inverse Radon transform is(10)fx,y=f^r,θ=12π2∫0π∫−∞+∞1rcosθ−β−ρ∂p∂ρdρdβ.Three-dimensional Radon transform is a generalization of two-dimensional, extending the line integral in two-dimensional to an area integral, each area integral corresponds to a point in Radon space. This point is the intersection of the plane with the normal to the plane through the origin. The three-dimensional Radon transform space is composed of all transform values.Figure3 shows the Fourier reconstruction method, which is also the derivation basis of the commonly used filtered back projection (FBP) reconstruction algorithm. As can be seen from the figure, one-dimensional Fourier transform is first performed for each projection data:(11)fx,y=∫−∞+∞Pβte−j2πωtdt.Figure 3
Schematic diagram of the central slice theorem.Pβt in formula (11) is to calculate the projection of f (x, y) on the β direction according to the Radon transform:(12)Pβt=∫−∞+∞∫−∞+∞fx,yδxcosβ+ysinβ−ρdxdy.


### 2.3. Iterative Reconstruction Algorithm

The filtered back projection (FBP) algorithm has certain limitations in practical application. For example, it requires that the projection data must be completely and uniformly distributed, the formula of the filtered back projection is continuous, and the image must be discretized during implementation. In this case, an iterative algorithm is a good choice.The biggest difference between the concept of the iterative reconstruction algorithm and the analytical reconstruction algorithm is that the former discretizes continuous imagesfr,θ. The algorithm divides the entire image area into J=n×n finite number of pixels, which are represented by f^r,θ. Figure 4 shows the process of reconstructing the image after the ray is discretized:Figure 4
Schematic diagram of discretized image projection.Among them,x1,x2,…,x9 represent the corresponding pixel values. It can be seen from the figure that the sum of the rays is(13)p1=w11x1+w12x2+⋯+w19x9,p2=w21x1+w22x2+⋯+w29x9,p3=w31x1+w32x2+⋯+w39x9,⋯p6=w61x1+w62x2+⋯+w69x9.Formula (13) can be expressed in a more compact form:(14)pi=∑j=19wijxji=1,2,…,6.Alternatively, it can be expressed in matrix representation:(15)p=Wx.In formula (15), p=p1,p2,…,p6T,x=x1,x2,…,x9T,W is the matrix of 6×9. Formulas (13) to (15) are derived from the special case of 9 pixels and 6 rays. We then generalize it to the general case and assume the general case of J pixels and I rays. At this time, x is the J-dimensional vector, which is called the image vector, I is the I-dimensional vector, which is called the measurement vector, and W is the I×J matrix, which is called the projection matrix.The task of iterative reconstruction is to findx according to the measured p, and the known projection matrix WW can be determined according to the system geometry, focal spot shape, detector response, and other physical parameters of the binocular vision system). wij in (13 to (15) represents the weighting factor of i-ray to j pixel, and there are many ways to calculate the weighting factor. In the simplest case, the weighting factor is set to 0 or 1 according to whether the ray passes through the pixel, as shown in the following formula:(16)wij=1,I−ray passes throughj−pixel,0,other.The ray can also be regarded as having a certain width, and the width is assumed to beτ (usually taking the pixel width δ). This thick line covers a part of the area of the pixel, and the ratio of the coverage area to the area of the pixel is the weighting factor of the pixel to the projection of the ray. For example, the gray value of the j pixel is xi, the area of the overlapping area between the i ray and the j pixel is Δs, and the ratio of it to the pixel area δ2 is wij=Δs/δ2, which is the weighting factor of the i ray to the j pixel.With the above foundation, the reconstruction problem of the binocular visual image can be transformed into a linear formula-solving process. The most intuitive way is to find the inverse matrixW−1 of matrix W, so as to get(17)x=W−1p.The second solution is to accumulate all the ray values passing through thej pixel to get the j pixel value:(18)xi=∑i=1Iwijpi,j=1,2,…,J.It can be written as a matrix problem as(19)x=WTp.Formula (19) is the form of back-projection reconstruction in the case of discrete pixels, and artifacts are very serious when using this reconstruction formula. However, formula (19) helps us understand the iterative reconstruction algorithm. So far, the main problems encountered by iterative reconstruction may be (1) generally, the number of pixels and the number of rays are extremely large, and it is difficult to directly find W−1. Even if W−1 is stored as a sparse matrix, it still requires a large amount of calculation; (2) in some cases, the number of projections is much smaller than the number of pixels, and the linear formula system may have an infinite number of solutions; (3) in the actual acquisition process, it may be affected by factors such as physical deviation or projection noise, and the reconstruction result cannot be obtained. Therefore, it is necessary to introduce an error value and estimate a set of solutions to make it optimal under a certain optimal criterion. Therefore, formula (15) can be modified as follows:(20)p=Wx+e.Here,e is the error vector, which can be measurement deviation and additional noise, such as detector electronics noise. According to the principle of numerical calculation and optimization, the estimation process can be implemented iteratively, which generally includes the following steps: (1) the image is discrete and initialized; (2) selection of iterative methods; (3) selection of optimal criteria. The iterative methods are as described above, mainly including classical iterations (such as ART, SIRT, and others) and statistical-based iterations (such as EM, MAP, and others). Step 3 also has more options: least squares criterion, maximum uniformity and smoothing criterion, maximum entropy criterion, and Bayesian criterion.


### 2.4. Evaluation Criteria for Reconstruction Algorithms

When using simulated data for testing, since the exact parameters of the simulated data are known, the reconstructed data can be accurately numerically compared with the original data to make an objective evaluation of the reconstruction quality. A commonly used simulation data model is the Shepp. LogaJl standard head model, which is composed of many ellipses with different sizes and densities. Common image numerical evaluation criteria are(1)
The measured valueε of the image similarity coefficient is as follows:(21)ε=∑i=1Nxi−x¯xi∗−x∗¯∑i=1Nxi−x¯2∑i=1Nxi∗−x∗¯21/2.(2)
The normalized RMS distance measurementd is as follows:(22)d=∑i=1Nxi−xi∗2∑i=1Nxi−x¯21/2.(3)
The normalized mean absolute distance measurementr is(23)r=∑i=1Nxi−xi∗∑i=1Nxi.(4)
The image signal-to-noise ratio (SNR) is(24)snr=10×lg∑i=1Nxi2∑i=1Nxi−xi′2.Among them, N is the number of pixels in the reconstructed image;xi is the gray value of the $x_{i}$ pixel in the model image; xi∗ is the gray value of the $x_{i}$ pixel in the reconstructed image; x¯ is the average gray value in the model image; x∗¯ is the average gray value in the reconstructed image.The four measurements above highlight different aspects of image quality. The image similarity coefficientε reflects the similarity between the reconstructed image and the simulated image. The larger the ε is, the more similar the two images are, and when ε is 1, the two images are identical. The normalized RMS distance measurement value d is more sensitive to reflect the error of the local situation, and if there is a large deviation of individual pixels, it will lead to a large d. The normalized mean absolute distance measurement value r is more sensitive to reflect the small error situation of most points. Contrary to d, it emphasizes the importance of more small errors rather than a small number of large errors. The signal-to-noise ratio (SNR) measures the ratio of the image signal to the noise signal, often expressed in decibels.


### 2.5. Calculation of Projection Matrix

The exact projection matrix plays a decisive role in reconstructing an image. Among all the methods for calculating projection matrices, the simplest model is the 0,1 model, as shown in Figure5(a). We assume that the projection matrix entry wij is 1 when the ray pj passes through the pixel, xi, and 0 otherwise. Figure 5(b) shows the projection matrix calculation based on length weighting. The value of the projection matrix entry wij is the value of the length of the ray pj intercepted by the pixel xi. Figure 5(c) shows the projection matrix calculation based on area weighting. In this model, the ray is regarded as having a certain width, and the value of the projection matrix item wij is the ratio of the area covered by the ray pj of the pixel xi to the pixel area. In practical applications, the internal attenuation of objects is continuous, and the attenuation values inside the discrete pixel cells are not completely equal. Therefore, the reconstructed image is only a discretized approximation of the real image, and the reconstructed image may be grainy. In theory, we can reconstruct a continuous description of the original image by using interpolation methods.Figure 5
Three basic projection matrix calculation methods. (a) Projection matrix based on 0, 1. (b) Length-based projection matrix. (c) Area-based projection matrix.
(a)(b)(c)As shown in Figure6, a ray pi is injected into the reconstruction area and is sampled at equal intervals, and the center of the interpolation kernel function hu,v is placed on the sampling point. All reconstructed pixels in the range of the revalued kernel are accumulated and appropriately weighted with the revalued kernel. Figure 6 shows that the sampled value sik at the point Xsik,Ysik is calculated from the adjacent pixel values. The value of sik is calculated by formula (6):(25)sik=∑ihXsik−Xxi,Ysik−Yxj×xi.(26)sik=∑jhXsik−Xxj,Ysik−Yxj×xj.Figure 6
Schematic diagram of the calculation of the sampling valuesik at the point Xsik,Ysik.The value of the projected pixel valuepi corresponding to the ray pi is the accumulation of all the sampling values sik along the ray:(27)pi=∑k∑jhXsik−Xxj,Ysik−Yxj×xj.Formula (26) is a discrete approximation of formula (27):(28)pi=∫∑jhXsik−Xxj,Ysik−Yxj×xjdsi.Formula (27) is rearranged, and we get(29)pi=∑jxj∫hXsik−Xxj,Ysik−Yxjdsi.Formula (28) is shown in Figure 7. Similar to formula (14), a projected pixel value pi is calculated as follows:(30)pi=∑kxj⋅wij.Figure 7
Schematic diagram of interpolation calculation of pixelxj and ray pi.Therefore, the weights are calculated by integrating the value kernel along the ray:(31)wij=∫hXsik−Xxj,Ysik−Yxjdsi.


### 2.6. Algebraic Reconstruction Method (ART)

The algebraic reconstruction algorithm (ART) is proposed to solve the problem of 3D object reconstruction. ART can be written as a linear geometry problem withWx=p. Here, x is an unknown (N × 1-dimensional) column vector to hold all N=n3 voxels in the reconstruction region of size n×n×n. p is an R-dimensional column vector, and $R$ is obtained by multiplying the number of pixels of each projection Rm by the number of projections M in the set of total projected images Pψ, that is, R=Rm×M. Pψ refers to the set of all projected images during one scan. $W$ is a R×N projection matrix, and the element wij in W represents the influence of the voxel vj on the ray ri. Wx=p can be written in the form of a linear system of formulas in formula (13). As mentioned above, it is very difficult to solve this linear system of formulas directly. So, here we introduce Kaczmarz’s method of solving this system of linear formulas.The expression for the update process of the reconstruction regionx is shown in formula (31):(32)xjk+1=xjk+λpi−∑n=1Nwinxnk∑n=1Nwin2wij.λ in formula (31) is a relaxation factor whose value range is in the (0,1] interval. But in general, if the λ value is too close to 0, it will lead to overoptimization. The algorithm calculates the formulas in formula (18) in sequence. After completion, some reconstructed pixels may not necessarily meet the convergence conditions, and the next iteration can be performed in the same way. Figure 8 shows the geometric process of the Kaczmarz method. The two straight lines in the figure can be expressed as two linear formulas, and the process shown in the figure is the solution process.Figure 8
Schematic diagram of the solution process of the Kaczmarz method.


### 2.7. Combined Reconstruction Method (SIRT)

The joint iteration method (SIRT) was proposed by Gilbert shortly after ART was proposed, and it is a parallel computing form of the ART algorithm. In this method, all the pixels of a certain projectionp are calculated first, and then the voxels of the entire reconstruction area are updated. Before the updated value is added to the voxel value, it needs to be weighted and normalized with the weighted value.As shown in formula (19), if w is a nonsingular matrix, its least squares solution is(33)x=WTW−1WTp.Among them,WTp represents the back-projection operation of p. If WTW−1 is viewed as a two-dimensional filter, formula (32) is the aforementioned two-dimensional filtered back-projection. Formula (32) can be transformed into(34)x=WTWWT−1p.In formula (33), WWT−1 is a one-dimensional ρ filter, which filters p. Formula (33) is solved iteratively:(35)x0=WTp,xk+1=xk+λkWTp−WTWxk=xk+λkWTp−Wxk.In formula (34), the projection value p is used as the initial value of the back projection. At the k + 1-th iteration, the algorithm uses the k-th iteration result xk to add the correction value to obtain xk+1, which is the back projection WTp−Wxk of the correction value and the kth-estimated error vector. Thus, the correction value for each voxel is the sum of the error values for all rays passing through that voxel, not just one ray. Therefore, the correction process of SIRT is called a point-by-point correction. This is the biggest difference from the ART algorithm and the fundamental reason why the SIRT algorithm can suppress noise; some random errors are averaged out by the common contribution of all rays passing through the voxel. In order to facilitate the iterative calculation, SIRT can also be rewritten as(36)xjk+1=xjk+λ∑pi∈Pφpi−∑n=1Nwinxnkwij.


### 2.8. Joint Algebraic Reconstruction Method (SART)

The SART algorithm does not correct each pixel (each ray) of the projected image, but it first calculates the projected image of the entire reconstructed area (denoted asPφ at the angle φ). Each pixel value in Pφ contributes to the corrected value of each voxel, and the updated value of each voxel is obtained by accumulating these contributions on each voxel. If the correction terms are simply added, the noise that may exist in the projected image will be added to the reconstructed image to produce artifacts, so weighting needs to be performed when updating. Figure 9 is a schematic diagram of the SART projection matrix calculation.Figure 9
Schematic diagram of SART projection matrix calculation.The way SART updates the reconstruction areax can be expressed in the following form:(37)xjk+1=xjk+∑pi∈Pφλpi−∑n=1Nwinxnk/∑n=1Nwinwij∑pi∈Pφwij.There are two significant differences between formulas (36) and (31) : 1. The correction term for a particular voxel xj is calculated by calculating the adjacent pixels pi in the projected image and weighting the influence of the pixel pi on each voxel by the coefficient wij. 2. Although the ART method guided by Kaczmarz’s method requires the sum of squared weights, the idea of SART regards ART as the inverse process of volume rendering.


## 2.1. Physical Basis of Binocular Vision Imaging

Robotic vision light is a high-energy electromagnetic wave with a certain energy and penetrating ability, which can penetrate some substances (such as human tissue) that visible light cannot pass through, as shown in Figure1. Generally, visible light has a longer wavelength, and when a photon hits an object, part of it is reflected and most of it is absorbed, while the wavelength of the robot’s vision light is extremely short, and the photon contains high energy. The penetration of the robot’s visual light is related to information such as the equivalent atomic number and density of the irradiated material, and the transmittance of the robot’s visual light is stronger for the material with a lower atomic number and vice versa. The transmission of robot vision light is an important basis for binocular vision imaging.Figure 1
Schematic diagram of the material that transmits light to the robot vision.A beam of robotic vision rays injected into a homogeneous material is considered, as shown in Figure2(a). Considering Beer’s theorem, we have(1)I=I0exp−μΔx⇒μΔx=lnII0.Figure 2
Schematic diagram of visual light passing through. (a) Schematic diagram of robot vision light passing through homogeneous material. (b) Schematic diagram of the robot’s vision of light passing through inhomogeneous materials.
(a)(b)It can be seen from formula (1) that the object with a high μ value causes more attenuation of the robot’s visual photons than the object with a low μ value. For example, the μ of bone is higher than that of soft tissue, indicating that it is more difficult for robotic vision photons to penetrate bone than soft tissue. On the other hand, the μ value of air is almost 0, indicating that the input and output of X rays hardly change on the path through the air.When the material scanned by theX-ray is inhomogeneous, the medium distributed along path 1 can be discretized into several continuous small blocks. When these small pieces are small enough, the medium inside the small pieces can be considered to be homogeneous and have the same attenuation coefficient. The thickness of each discrete small block is assumed to be Δx, and the attenuation coefficients of each discrete small block are μ0,μ1,μ2,μ3,…,μn, respectively, as shown in Figure 2.The ray intensity of the robot vision light after passing through the first small block isI1, and the ray intensity after passing through the second block is I2….The final transmission intensity isIn, then there are(2)I1=I0exp−μ1Δx,(3)I2=I1exp−μ2Δx.Substituting formulas (2) into (3), we get(4)I2=I0exp−μ1+μ2Δx.It continues to accumulate attenuation values along the propagation direction of the robot’s vision light until the final transmission intensityIn of the robot’s vision light when it leaves the illuminated object, as shown in the following equation:(5)In=I0exp−μ1+μ2+⋯+μnΔx.We take the positive exponent of formula (5) and express it in summed form, we get(6)p=∑i=1nμiΔx=lnI0In.p in formula (6) is the projection. If the incident intensity I0 and outgoing intensity In of the X ray are known, a linear formula with μii=1,2,…,n as the last known is obtained according to formula (6). When Δx⟶0, formula (6) can represent the summation of continuous variation, and its integral form is(7)p=lnI0In=∫Lμl.In formula (7), μl is a continuous function of the decay rate with respect to path l. The process of finding the attenuation coefficient function μl by projection is called back projection. If the two-dimensional density function f (x, y) is used to describe the attenuation rate of the two-dimensional plane, then the problem of binocular vision imaging can be expressed as the measured linear integral of an object is given, and the decay rate at each point is calculated to produce two-dimensional density data.


## 2.2. Analytical Reconstruction Algorithm

Two-dimensional Radon transform is a projection transform of straight line integral, and it can be defined in many forms. This paper takes the most commonly used case as an example.We assume that in the plane areaR2, any point (x, y) can be represented by polar coordinates ρ,θ, ρ represents the distance from the line to the origin, θ represents the angle between the straight line l and the positive y-axis, and the function f (x, y) is the image to be reconstructed. We assume that any straight line l∈R2, then the two-dimensional Radon transform of the function f (x, y) is defined as(8)p=∫lfx,ydl=∫lf^r,θdl=∫−∞+∞f^L2+l2,ϕ+tan−1lLdl.The formula of the straight line can be expressed asρ=xcosβ+ysinβ in polar coordinates, then formula (8) can be further expressed as(9)p=∫−∞+∞fx,ydl=∫−∞+∞∫−∞+∞fx,yδxcosβ+ysinβ−ρdxdy.Among them, (x, y) represents the position of the reconstructed pixel in the Cartesian coordinate system, and δx represents the sampling function. This process is the integration of the image in a straight line. According to the previous physical principle of the X-ray, we can abstract the attenuation process of the X-ray into this integration process, and the projection value of the X-ray is the value of the Radon transform at angle β. If the projection is performed at multiple angles, the Radon value of each angle can be obtained, and the two-dimensional image f (x, y) of the original plane can be reconstructed by performing the Radon inverse transformation on these obtained Radon values.The formula for the two-dimensional inverse Radon transform is(10)fx,y=f^r,θ=12π2∫0π∫−∞+∞1rcosθ−β−ρ∂p∂ρdρdβ.Three-dimensional Radon transform is a generalization of two-dimensional, extending the line integral in two-dimensional to an area integral, each area integral corresponds to a point in Radon space. This point is the intersection of the plane with the normal to the plane through the origin. The three-dimensional Radon transform space is composed of all transform values.Figure3 shows the Fourier reconstruction method, which is also the derivation basis of the commonly used filtered back projection (FBP) reconstruction algorithm. As can be seen from the figure, one-dimensional Fourier transform is first performed for each projection data:(11)fx,y=∫−∞+∞Pβte−j2πωtdt.Figure 3
Schematic diagram of the central slice theorem.Pβt in formula (11) is to calculate the projection of f (x, y) on the β direction according to the Radon transform:(12)Pβt=∫−∞+∞∫−∞+∞fx,yδxcosβ+ysinβ−ρdxdy.


## 2.3. Iterative Reconstruction Algorithm

The filtered back projection (FBP) algorithm has certain limitations in practical application. For example, it requires that the projection data must be completely and uniformly distributed, the formula of the filtered back projection is continuous, and the image must be discretized during implementation. In this case, an iterative algorithm is a good choice.The biggest difference between the concept of the iterative reconstruction algorithm and the analytical reconstruction algorithm is that the former discretizes continuous imagesfr,θ. The algorithm divides the entire image area into J=n×n finite number of pixels, which are represented by f^r,θ. Figure 4 shows the process of reconstructing the image after the ray is discretized:Figure 4
Schematic diagram of discretized image projection.Among them,x1,x2,…,x9 represent the corresponding pixel values. It can be seen from the figure that the sum of the rays is(13)p1=w11x1+w12x2+⋯+w19x9,p2=w21x1+w22x2+⋯+w29x9,p3=w31x1+w32x2+⋯+w39x9,⋯p6=w61x1+w62x2+⋯+w69x9.Formula (13) can be expressed in a more compact form:(14)pi=∑j=19wijxji=1,2,…,6.Alternatively, it can be expressed in matrix representation:(15)p=Wx.In formula (15), p=p1,p2,…,p6T,x=x1,x2,…,x9T,W is the matrix of 6×9. Formulas (13) to (15) are derived from the special case of 9 pixels and 6 rays. We then generalize it to the general case and assume the general case of J pixels and I rays. At this time, x is the J-dimensional vector, which is called the image vector, I is the I-dimensional vector, which is called the measurement vector, and W is the I×J matrix, which is called the projection matrix.The task of iterative reconstruction is to findx according to the measured p, and the known projection matrix WW can be determined according to the system geometry, focal spot shape, detector response, and other physical parameters of the binocular vision system). wij in (13 to (15) represents the weighting factor of i-ray to j pixel, and there are many ways to calculate the weighting factor. In the simplest case, the weighting factor is set to 0 or 1 according to whether the ray passes through the pixel, as shown in the following formula:(16)wij=1,I−ray passes throughj−pixel,0,other.The ray can also be regarded as having a certain width, and the width is assumed to beτ (usually taking the pixel width δ). This thick line covers a part of the area of the pixel, and the ratio of the coverage area to the area of the pixel is the weighting factor of the pixel to the projection of the ray. For example, the gray value of the j pixel is xi, the area of the overlapping area between the i ray and the j pixel is Δs, and the ratio of it to the pixel area δ2 is wij=Δs/δ2, which is the weighting factor of the i ray to the j pixel.With the above foundation, the reconstruction problem of the binocular visual image can be transformed into a linear formula-solving process. The most intuitive way is to find the inverse matrixW−1 of matrix W, so as to get(17)x=W−1p.The second solution is to accumulate all the ray values passing through thej pixel to get the j pixel value:(18)xi=∑i=1Iwijpi,j=1,2,…,J.It can be written as a matrix problem as(19)x=WTp.Formula (19) is the form of back-projection reconstruction in the case of discrete pixels, and artifacts are very serious when using this reconstruction formula. However, formula (19) helps us understand the iterative reconstruction algorithm. So far, the main problems encountered by iterative reconstruction may be (1) generally, the number of pixels and the number of rays are extremely large, and it is difficult to directly find W−1. Even if W−1 is stored as a sparse matrix, it still requires a large amount of calculation; (2) in some cases, the number of projections is much smaller than the number of pixels, and the linear formula system may have an infinite number of solutions; (3) in the actual acquisition process, it may be affected by factors such as physical deviation or projection noise, and the reconstruction result cannot be obtained. Therefore, it is necessary to introduce an error value and estimate a set of solutions to make it optimal under a certain optimal criterion. Therefore, formula (15) can be modified as follows:(20)p=Wx+e.Here,e is the error vector, which can be measurement deviation and additional noise, such as detector electronics noise. According to the principle of numerical calculation and optimization, the estimation process can be implemented iteratively, which generally includes the following steps: (1) the image is discrete and initialized; (2) selection of iterative methods; (3) selection of optimal criteria. The iterative methods are as described above, mainly including classical iterations (such as ART, SIRT, and others) and statistical-based iterations (such as EM, MAP, and others). Step 3 also has more options: least squares criterion, maximum uniformity and smoothing criterion, maximum entropy criterion, and Bayesian criterion.


## 2.4. Evaluation Criteria for Reconstruction Algorithms

When using simulated data for testing, since the exact parameters of the simulated data are known, the reconstructed data can be accurately numerically compared with the original data to make an objective evaluation of the reconstruction quality. A commonly used simulation data model is the Shepp. LogaJl standard head model, which is composed of many ellipses with different sizes and densities. Common image numerical evaluation criteria are(1)
The measured valueε of the image similarity coefficient is as follows:(21)ε=∑i=1Nxi−x¯xi∗−x∗¯∑i=1Nxi−x¯2∑i=1Nxi∗−x∗¯21/2.(2)
The normalized RMS distance measurementd is as follows:(22)d=∑i=1Nxi−xi∗2∑i=1Nxi−x¯21/2.(3)
The normalized mean absolute distance measurementr is(23)r=∑i=1Nxi−xi∗∑i=1Nxi.(4)
The image signal-to-noise ratio (SNR) is(24)snr=10×lg∑i=1Nxi2∑i=1Nxi−xi′2.Among them, N is the number of pixels in the reconstructed image;xi is the gray value of the $x_{i}$ pixel in the model image; xi∗ is the gray value of the $x_{i}$ pixel in the reconstructed image; x¯ is the average gray value in the model image; x∗¯ is the average gray value in the reconstructed image.The four measurements above highlight different aspects of image quality. The image similarity coefficientε reflects the similarity between the reconstructed image and the simulated image. The larger the ε is, the more similar the two images are, and when ε is 1, the two images are identical. The normalized RMS distance measurement value d is more sensitive to reflect the error of the local situation, and if there is a large deviation of individual pixels, it will lead to a large d. The normalized mean absolute distance measurement value r is more sensitive to reflect the small error situation of most points. Contrary to d, it emphasizes the importance of more small errors rather than a small number of large errors. The signal-to-noise ratio (SNR) measures the ratio of the image signal to the noise signal, often expressed in decibels.


## 2.5. Calculation of Projection Matrix

The exact projection matrix plays a decisive role in reconstructing an image. Among all the methods for calculating projection matrices, the simplest model is the 0,1 model, as shown in Figure5(a). We assume that the projection matrix entry wij is 1 when the ray pj passes through the pixel, xi, and 0 otherwise. Figure 5(b) shows the projection matrix calculation based on length weighting. The value of the projection matrix entry wij is the value of the length of the ray pj intercepted by the pixel xi. Figure 5(c) shows the projection matrix calculation based on area weighting. In this model, the ray is regarded as having a certain width, and the value of the projection matrix item wij is the ratio of the area covered by the ray pj of the pixel xi to the pixel area. In practical applications, the internal attenuation of objects is continuous, and the attenuation values inside the discrete pixel cells are not completely equal. Therefore, the reconstructed image is only a discretized approximation of the real image, and the reconstructed image may be grainy. In theory, we can reconstruct a continuous description of the original image by using interpolation methods.Figure 5
Three basic projection matrix calculation methods. (a) Projection matrix based on 0, 1. (b) Length-based projection matrix. (c) Area-based projection matrix.
(a)(b)(c)As shown in Figure6, a ray pi is injected into the reconstruction area and is sampled at equal intervals, and the center of the interpolation kernel function hu,v is placed on the sampling point. All reconstructed pixels in the range of the revalued kernel are accumulated and appropriately weighted with the revalued kernel. Figure 6 shows that the sampled value sik at the point Xsik,Ysik is calculated from the adjacent pixel values. The value of sik is calculated by formula (6):(25)sik=∑ihXsik−Xxi,Ysik−Yxj×xi.(26)sik=∑jhXsik−Xxj,Ysik−Yxj×xj.Figure 6
Schematic diagram of the calculation of the sampling valuesik at the point Xsik,Ysik.The value of the projected pixel valuepi corresponding to the ray pi is the accumulation of all the sampling values sik along the ray:(27)pi=∑k∑jhXsik−Xxj,Ysik−Yxj×xj.Formula (26) is a discrete approximation of formula (27):(28)pi=∫∑jhXsik−Xxj,Ysik−Yxj×xjdsi.Formula (27) is rearranged, and we get(29)pi=∑jxj∫hXsik−Xxj,Ysik−Yxjdsi.Formula (28) is shown in Figure 7. Similar to formula (14), a projected pixel value pi is calculated as follows:(30)pi=∑kxj⋅wij.Figure 7
Schematic diagram of interpolation calculation of pixelxj and ray pi.Therefore, the weights are calculated by integrating the value kernel along the ray:(31)wij=∫hXsik−Xxj,Ysik−Yxjdsi.


## 2.6. Algebraic Reconstruction Method (ART)

The algebraic reconstruction algorithm (ART) is proposed to solve the problem of 3D object reconstruction. ART can be written as a linear geometry problem withWx=p. Here, x is an unknown (N × 1-dimensional) column vector to hold all N=n3 voxels in the reconstruction region of size n×n×n. p is an R-dimensional column vector, and $R$ is obtained by multiplying the number of pixels of each projection Rm by the number of projections M in the set of total projected images Pψ, that is, R=Rm×M. Pψ refers to the set of all projected images during one scan. $W$ is a R×N projection matrix, and the element wij in W represents the influence of the voxel vj on the ray ri. Wx=p can be written in the form of a linear system of formulas in formula (13). As mentioned above, it is very difficult to solve this linear system of formulas directly. So, here we introduce Kaczmarz’s method of solving this system of linear formulas.The expression for the update process of the reconstruction regionx is shown in formula (31):(32)xjk+1=xjk+λpi−∑n=1Nwinxnk∑n=1Nwin2wij.λ in formula (31) is a relaxation factor whose value range is in the (0,1] interval. But in general, if the λ value is too close to 0, it will lead to overoptimization. The algorithm calculates the formulas in formula (18) in sequence. After completion, some reconstructed pixels may not necessarily meet the convergence conditions, and the next iteration can be performed in the same way. Figure 8 shows the geometric process of the Kaczmarz method. The two straight lines in the figure can be expressed as two linear formulas, and the process shown in the figure is the solution process.Figure 8
Schematic diagram of the solution process of the Kaczmarz method.


## 2.7. Combined Reconstruction Method (SIRT)

The joint iteration method (SIRT) was proposed by Gilbert shortly after ART was proposed, and it is a parallel computing form of the ART algorithm. In this method, all the pixels of a certain projectionp are calculated first, and then the voxels of the entire reconstruction area are updated. Before the updated value is added to the voxel value, it needs to be weighted and normalized with the weighted value.As shown in formula (19), if w is a nonsingular matrix, its least squares solution is(33)x=WTW−1WTp.Among them,WTp represents the back-projection operation of p. If WTW−1 is viewed as a two-dimensional filter, formula (32) is the aforementioned two-dimensional filtered back-projection. Formula (32) can be transformed into(34)x=WTWWT−1p.In formula (33), WWT−1 is a one-dimensional ρ filter, which filters p. Formula (33) is solved iteratively:(35)x0=WTp,xk+1=xk+λkWTp−WTWxk=xk+λkWTp−Wxk.In formula (34), the projection value p is used as the initial value of the back projection. At the k + 1-th iteration, the algorithm uses the k-th iteration result xk to add the correction value to obtain xk+1, which is the back projection WTp−Wxk of the correction value and the kth-estimated error vector. Thus, the correction value for each voxel is the sum of the error values for all rays passing through that voxel, not just one ray. Therefore, the correction process of SIRT is called a point-by-point correction. This is the biggest difference from the ART algorithm and the fundamental reason why the SIRT algorithm can suppress noise; some random errors are averaged out by the common contribution of all rays passing through the voxel. In order to facilitate the iterative calculation, SIRT can also be rewritten as(36)xjk+1=xjk+λ∑pi∈Pφpi−∑n=1Nwinxnkwij.


## 2.8. Joint Algebraic Reconstruction Method (SART)

The SART algorithm does not correct each pixel (each ray) of the projected image, but it first calculates the projected image of the entire reconstructed area (denoted asPφ at the angle φ). Each pixel value in Pφ contributes to the corrected value of each voxel, and the updated value of each voxel is obtained by accumulating these contributions on each voxel. If the correction terms are simply added, the noise that may exist in the projected image will be added to the reconstructed image to produce artifacts, so weighting needs to be performed when updating. Figure 9 is a schematic diagram of the SART projection matrix calculation.Figure 9
Schematic diagram of SART projection matrix calculation.The way SART updates the reconstruction areax can be expressed in the following form:(37)xjk+1=xjk+∑pi∈Pφλpi−∑n=1Nwinxnk/∑n=1Nwinwij∑pi∈Pφwij.There are two significant differences between formulas (36) and (31) : 1. The correction term for a particular voxel xj is calculated by calculating the adjacent pixels pi in the projected image and weighting the influence of the pixel pi on each voxel by the coefficient wij. 2. Although the ART method guided by Kaczmarz’s method requires the sum of squared weights, the idea of SART regards ART as the inverse process of volume rendering.


## 3. Bolt Positioning of Industrial Robots Based on Multieye Vision

The whole scheme is centered on the robot, with the end effector with grasping and prescrewing functions. The control system is responsible for the implementation of trajectory motion, grasping, and prescrewing, and the vision system realizes the positioning of parts. The whole scheme is shown in Figure10.Figure 10
Robot system scheme composition.Common robot types include Cartesian robots, cylindrical coordinate robots, articulated robots, SCARA robots, spherical robots, and so on, as shown in Figure11. The articulated six-degree-of-freedom robot is the most common one and plays an important role in industrial production. The SCARA plane articulated robot has the least interference in space and has the optimal structural solution.Figure 11
Schematic diagram of the workflow of the robot.Based on the above, the effect of the bolt positioning system based on the multivision industrial robot proposed in this paper is verified, and the bolt positioning accuracy is calculated. A total of 36 groups were marked with 1000 bolts in each group, and the test results were obtained as shown in Table1 and Figure 12.Table 1
Statistical table of the accuracy of bolt positioning based on multivision industrial robots.
NumAccuracy (%)198.6297.8397.4499.2599.5697.9797.2899.7997.21097.91198.01299.11398.914100.01597.21697.31799.51898.61998.52098.72199.42298.02398.92497.12598.12697.72797.82898.92997.83097.93199.93299.93398.33498.83598.13698.2Figure 12
Statistical diagram of the accuracy of bolt positioning based on multieye vision industrial robots.From the above research, it can be seen that the bolt positioning system based on the multieye vision industrial robot proposed in this paper has a good performance in bolt positioning.


## 4. Conclusion

Binocular vision measurement technology has been relatively mature and widely used in industrial fields. However, there are few types of research on the application of the bolt assembly problem. At the same time, most of the current research studies on the rotation matching of the bolts use the method of pressing and rotating, as well as autonomously matching within the rotation range. There is a large impact, and the bolts placed in the through holes cannot be grasped. This paper studies the bolt positioning system of an industrial robot combined with multieye vision technology. This research extracts important theoretical support for the bolt grasping process of the bolt assembly robot, which is of great significance for promoting the automation process in the field of industrial assembly.



---

*Source: 7563199-2022-08-24.xml*