Figure 5a,b show the closeup view photographs of the machine vision module under its two actual usage conditions, while their corresponding working principles are illustrated in Figure 6a,b, respectively. As shown in Figures 5a and 6a, when the machine vision module is applied to visually position the drilltip of an original microdrill sample (relative to the end face of the grinding wheel) before the grinding off operation, the microdrill sample is moved to a prescribed position where its drilltip (and its central axis) would align with the outer boundary edge of the grinding wheel end-face (GWEF), i.e., the edge point G, but also keep a very small distance along the Y-direction. Meanwhile, the optical subset is synchronously moved to the front of the telecentric lens, where the X-directional distance between point G and the primary principle point L of the set of the PCA lenses would equal the effective focal length *f* of the PCA lenses (whose nominal value is 50 mm). That is, the ideal condition of GL = *f* = 50 mm would exist. When an object is placed at a Y-Z plane where the effective focal point of the PCA lenses (i.e., point G at that ideal condition) is located, a half-size virtual image of that object would be exactly formed at another Y-Z plane where point V is located, for which the ideal condition of VL = *f*/2 = 25 mm would exist. When the parallel light beam is projected from the backlight illuminator, it would be reflected by the reflection mirror (set at an inclination angle of 45◦ ) and would then reach the telecentric lens and finally be sensed by the CMOS camera. A ray passing through points G, V, and L would reach a reference point R on the reflection mirror and would then be perpendicularly reflected to reach a reference point T on the telecentric lens. The ideal condition of this optical design is that the sum of the three distances VL, LR (denoted by *d*LR), and RT (denoted by *d*RT) would equal the working distance of the telecentric lens (denoted by *dw* with a nominal value of 64.9 mm), that is, the ideal condition of VL + LR + RT = *f*/2 + *d*LR + *d*RT = *dw* = 64.9 mm or LR + RT = *d*LR + *d*RT = 39.9 mm would exist. Therefore, when the ring illuminator is not switched on, a half-size virtual image of the dark shadows of the drilltip of the microdrill sample and of part of the grinding wheel would be sensed and captured by the CMOS camera. Such a virtual image of dark shadows (called a side-viewed image hereinafter) would be captured under a resultant

magnification of 8, because of the optical effect of the PCA lenses, and would further be used to position the drilltip (i.e., to measure the axial distance between the drilltip and the GWEF). Due to a lower magnification, the area of the resultant field of view (FOV) of this optical design would be four times that of the camera-lens subassembly itself; such a larger FOV should be more appropriate for the application of positioning the drilltip. In other words, a larger FOV means that a larger axial distance between the drilltip and the grinding wheel could be maintained when capturing their side-viewed images, so that the interference between the two objects could be more easily avoided. ant magnification of 8, because of the optical effect of the PCA lenses, and would further be used to position the drilltip (i.e., to measure the axial distance between the drilltip and the GWEF). Due to a lower magnification, the area of the resultant field of view (FOV) of this optical design would be four times that of the camera-lens subassembly itself; such a larger FOV should be more appropriate for the application of positioning the drilltip. In other words, a larger FOV means that a larger axial distance between the drilltip and the grinding wheel could be maintained when capturing their side-viewed images, so that the interference between the two objects could be more easily avoided.

condition of VL + LR + RT = f/2 + dLR + dRT = dw = 64.9 mm or LR + RT = dLR + dRT = 39.9 mm would exist. Therefore, when the ring illuminator is not switched on, a half-size virtual image of the dark shadows of the drilltip of the microdrill sample and of part of the grinding wheel would be sensed and captured by the CMOS camera. Such a virtual image of dark shadows (called a side-viewed image hereinafter) would be captured under a result-

Machines 2023, 11, x FOR PEER REVIEW 8 of 34

Figure 5. Closeup views of the machine vision module under its actual usage conditions: (a) the application for positioning the drilltip; (b) the application for measuring the CSGPs. **Figure 5.** Closeup views of the machine vision module under its actual usage conditions: (**a**) the application for positioning the drilltip; (**b**) the application for measuring the CSGPs.

As shown in Figures 5b and 6b, when the machine vision module is applied to visually measure the CSGPs of a ground microdrill sample, the microdrill sample would be moved to another prescribed position where its ground cross-section would be located at the front of the telecentric lens and be focused properly; for which, the axial (Y-directional) distance between the ground cross-section and the telecentric lens would ideally equal the nominal working distance of the telecentric lens (i.e., dw = 64.9 mm). Meanwhile, the optical subset would be synchronously moved to keep a proper distance from the front of the telecentric lens. When the ring illuminator is switched on, part of the diffused light beams projected from it would reach the ground cross-section and would then be reflected into the telecentric lens and finally be sensed by the CMOS camera. Therefore, a bright surface image of the ground cross-section (called an axial-viewed image hereinafter) would be sensed and captured by the CMOS camera under a magnification of 16. Such a high magnification would lead to a tiny FOV and a fine resolution, which should be appropriate for the application of measuring the CSGPs. In addition, the parallel light beam projected from the backlight illuminator would not influence the measurement, since it would not illuminate the ground cross-section and not be reflected into the telecentric lens by the optical subset. Thereby, the backlight illuminator would not need to be switched off under this usage condition. As a result, the setup of the machine vision module, characterized by several conventional machine vision components in combination with the innovative and lower cost optical subset, could enable the drilltip to be visually positioned and the As shown in Figures 5b and 6b, when the machine vision module is applied to visually measure the CSGPs of a ground microdrill sample, the microdrill sample would be moved to another prescribed position where its ground cross-section would be located at the front of the telecentric lens and be focused properly; for which, the axial (Y-directional) distance between the ground cross-section and the telecentric lens would ideally equal the nominal working distance of the telecentric lens (i.e., *dw* = 64.9 mm). Meanwhile, the optical subset would be synchronously moved to keep a proper distance from the front of the telecentric lens. When the ring illuminator is switched on, part of the diffused light beams projected from it would reach the ground cross-section and would then be reflected into the telecentric lens and finally be sensed by the CMOS camera. Therefore, a bright surface image of the ground cross-section (called an axial-viewed image hereinafter) would be sensed and captured by the CMOS camera under a magnification of 16. Such a high magnification would lead to a tiny FOV and a fine resolution, which should be appropriate for the application of measuring the CSGPs. In addition, the parallel light beam projected from the backlight illuminator would not influence the measurement, since it would not illuminate the ground cross-section and not be reflected into the telecentric lens by the optical subset. Thereby, the backlight illuminator would not need to be switched off under this usage condition. As a result, the setup of the machine vision module, characterized by several conventional machine vision components in combination with the innovative and lower cost optical subset, could enable the drilltip to be visually positioned and the CSGPs to be visually measured via captured side- and axial-viewed images, respectively.

CSGPs to be visually measured via captured side- and axial-viewed images, respectively.

Figure 6. Illustrations of the working principles of the machine vision module: (a) the application for positioning the drilltip; (b) the application for measuring the CSGPs. **Figure 6.** Illustrations of the working principles of the machine vision module: (**a**) the application for positioning the drilltip; (**b**) the application for measuring the CSGPs.

#### 2.2. System Integration *2.2. System Integration*

The entire functional block diagram of the developed system is depicted in Figure 7. For accomplishing the system integration, a host personal computer (PC) was prepared to control, manipulate, and monitor the hardware modules of the entire machine. The functional connections between the host PC and the hardware modules were divided into three subsets: the motion control subset, the switching subset, and the image acquisition The entire functional block diagram of the developed system is depicted in Figure 7. For accomplishing the system integration, a host personal computer (PC) was prepared to control, manipulate, and monitor the hardware modules of the entire machine. The functional connections between the host PC and the hardware modules were divided into three subsets: the motion control subset, the switching subset, and the image acquisition

subset. The motion control subset was applied to manipulate the stepping motors of the dual-axis motion module via a National Instruments (NI) PCI-7340 motion control card installed on the motherboard of the host PC. The motion control card was electrically connected to the drivers of the two stepping motors. The stepping motor of the X-axis LMT was driven via an OM VEXTA RKD514L-A stepping motor driver, and that of the Y-axis LMT via an OM VEXTA RKD507-A stepping motor driver, while their angular resolutions were set to 5000 and 2000 step/rev, respectively. Thereby, the X- and Y-axis LMTs could both achieve a fine positioning resolution of 1 µm/step according to the nominal leads of their ballscrew units. The signals of the installed optical linear encoders were fed back to the motion control card for sensing the reference positions of both LMTs as well as for executing closed-loop control, and the extreme positions of both LMTs were detected by Renishaw A-9531-0251 limit switches (magnets) installed in them. subset. The motion control subset was applied to manipulate the stepping motors of the dual-axis motion module via a National Instruments (NI) PCI-7340 motion control card installed on the motherboard of the host PC. The motion control card was electrically connected to the drivers of the two stepping motors. The stepping motor of the X-axis LMT was driven via an OM VEXTA RKD514L-A stepping motor driver, and that of the Y-axis LMT via an OM VEXTA RKD507-A stepping motor driver, while their angular resolutions were set to 5000 and 2000 step/rev, respectively. Thereby, the X- and Y-axis LMTs could both achieve a fine positioning resolution of 1 µm/step according to the nominal leads of their ballscrew units. The signals of the installed optical linear encoders were fed back to the motion control card for sensing the reference positions of both LMTs as well as for executing closed-loop control, and the extreme positions of both LMTs were detected by Renishaw A-9531-0251 limit switches (magnets) installed in them.

Figure 7. Functional block diagram of the improved automation system. **Figure 7.** Functional block diagram of the improved automation system.

The switching subset was applied to switch on/off the induction motor and the ring illuminator via an NI PCI-6221 data acquisition (DAQ) card installed on the motherboard of the host PC. The DAQ card was electrically connected to a relay unit aĴached to the induction motor and the light source regulator of the ring illuminator. (Two customized Fadracer Technology VR-ML1505 type light source regulators were used to power and adjust the ring and backlight illuminators, respectively; that of the backlight illuminator did not need to be connected to and manipulated by the switching subset.) Two independent analog voltage signals were generated and outpuĴed by the DAQ card for triggering/untriggering the relay unit and the light source regulator, respectively, thereby the induction motor and the ring illuminator were accordingly switched on/off. The switching subset was applied to switch on/off the induction motor and the ring illuminator via an NI PCI-6221 data acquisition (DAQ) card installed on the motherboard of the host PC. The DAQ card was electrically connected to a relay unit attached to the induction motor and the light source regulator of the ring illuminator. (Two customized Fadracer Technology VR-ML1505 type light source regulators were used to power and adjust the ring and backlight illuminators, respectively; that of the backlight illuminator did not need to be connected to and manipulated by the switching subset.) Two independent analog voltage signals were generated and outputted by the DAQ card for triggering/untriggering the relay unit and the light source regulator, respectively, thereby the induction motor and the ring illuminator were accordingly switched on/off.

The image acquisition subset was applied to manipulate the CMOS camera for capturing the required digital images and then storing them in the host PC, for which the CMOS camera was electrically connected to the host PC via a universal serial bus (USB) port. Each captured digital image (a side- or axial-viewed image) consisted of a 2592 1944 array of pixels with 8-bit grayscale values ranging from 0 (black) to 255 (white). For a captured side-viewed image, its conversion factor corresponding to the physical dimension was found to be 0.323 µm/pixel by means of spatial calibration for computer vision [29,30]; that for a captured axial-viewed image was found to be 0.165 µm/pixel. Thereby, real-world FOV areas of the captured side- or axial-viewed images were approximately The image acquisition subset was applied to manipulate the CMOS camera for capturing the required digital images and then storing them in the host PC, for which the CMOS camera was electrically connected to the host PC via a universal serial bus (USB) port. Each captured digital image (a side- or axial-viewed image) consisted of a 2592 × 1944 array of pixels with 8-bit grayscale values ranging from 0 (black) to 255 (white). For a captured sideviewed image, its conversion factor corresponding to the physical dimension was found to be 0.323 µm/pixel by means of spatial calibration for computer vision [29,30]; that for a captured axial-viewed image was found to be 0.165 µm/pixel. Thereby, real-world FOV areas of the captured side- or axial-viewed images were approximately 837 × 628 or 428 × 321 µm2 , respectively. Additionally, by means of the sub-pixel localization algorithm [29,30], the-

oretically estimated resolutions of 1/25 pixel for the side- and axial-viewed images (i.e., 0.0129 and 0.0066 µm, respectively) could be achieved. Two captured side-viewed images are shown in Figure 8a,b, in which dark shadows of part of the grinding wheel and of the drilltip of an original or ground microdrill sample with a bright background can be clearly observed. Two captured axial-viewed images are shown in Figure 8c,d, in which ground cross-sections of a UC or ST type microdrill sample (the bright objects) with a dark background can be clearly observed. axial-viewed images (i.e., 0.0129 and 0.0066 µm, respectively) could be achieved. Two captured side-viewed images are shown in Figure 8a,b, in which dark shadows of part of the grinding wheel and of the drilltip of an original or ground microdrill sample with a bright background can be clearly observed. Two captured axial-viewed images are shown in Figure 8c,d, in which ground cross-sections of a UC or ST type microdrill sample (the bright objects) with a dark background can be clearly observed.

tion algorithm [29,30], theoretically estimated resolutions of 1/25 pixel for the side- and

, respectively. Additionally, by means of the sub-pixel localiza-

Machines 2023, 11, x FOR PEER REVIEW 11 of 34

837 628 or 428 321 µm2

 Figure 8. Digital images captured by the machine vision module: (a) a side-viewed image of an **Figure 8.** Digital images captured by the machine vision module: (**a**) a side-viewed image of an original microdrill sample; (**b**) another side-viewed image of a ground microdrill sample; (**c**) an axialviewed image of a ground cross-section; (**d**) an axial-viewed image of another ground cross-section.

original microdrill sample; (b) another side-viewed image of a ground microdrill sample; (c) an ax-

ial-viewed image of a ground cross-section; (d) an axial-viewed image of another ground crosssection. After the three subsets had been successfully set up, the required human-machine interface software and key programs (for motion control, logic control, and image pro-After the three subsets had been successfully set up, the required human-machine interface software and key programs (for motion control, logic control, and image processing) were developed, integrated, and tested in the NI LabVIEW environment, which enabled the constructed hardware system to be applied to execute an AMP that is introduced in Section 3.

#### cessing) were developed, integrated, and tested in the NI LabVIEW environment, which **3. The Automated Measuring Process (AMP)**

enabled the constructed hardware system to be applied to execute an AMP that is introduced in Section 3. 3. The Automated Measuring Process (AMP) In order to establish an AMP that is suitable for measuring the CSWTs and CSODs of both ST and UC type microdrills with the use of the improved automation system, the In order to establish an AMP that is suitable for measuring the CSWTs and CSODs of both ST and UC type microdrills with the use of the improved automation system, the AMP proposed by Chang et al. [28] is extended and improved in this work. As a result, a more comprehensive AMP with specific image processing operations, which primarily consists of five essential sub-processes, is established, and its main flowchart is shown in Figure 9. The main steps in this flowchart are described as follows:

AMP proposed by Chang et al. [28] is extended and improved in this work. As a result, a more comprehensive AMP with specific image processing operations, which primarily consists of five essential sub-processes, is established, and its main flowchart is shown in

As can be seen, the DP, GO, AF, WTM, and ODM sub-processes must be included in the established AMP. Although the DP, GO, and WTM sub-processes are adopted in Chang et al.'s AMP [28], the AF and ODM sub-processes are newly added ones, while the WTM sub-process presented in this work is an improved one with better applicability.

Figure 9. Main flowchart of the established AMP. **Figure 9.** Main flowchart of the established AMP.

#### 3.1. Drilltip Positioning (DP) Sub-Process *3.1. Drilltip Positioning (DP) Sub-Process*

This sub-process is basically identical to that adopted by Chang et al. [28]. When the microdrill sample has been moved to the first prescribed position (during Step 4), it is then moved a small axial distance (whose value was set to 50 µm (≈154.80 pixels in a sideviewed image)) toward the GWEF (via the Y-axis LMT) and a side-viewed image is captured and stored as a temporary file. The two captured images shown in Figure 8a,b cor-This sub-process is basically identical to that adopted by Chang et al. [28]. When the microdrill sample has been moved to the first prescribed position (during Step 4), it is then moved a small axial distance (whose value was set to 50 µm (≈154.80 pixels in a side-viewed image)) toward the GWEF (via the Y-axis LMT) and a side-viewed image is captured and stored as a temporary file. The two captured images shown in Figure 8a,b correspond to the counting index *i* = 1 and *i* > 1, respectively.

respond to the counting index i 1 and i > 1, respectively. Subsequently, edge detection for the captured side-viewed image is performed, as shown in Figure 10a or Figure 10b; a rectangular region of interest (ROI), numbered ROI #1, is assigned to appropriately cover the shadows of the two objects appearing in the FOV. The concept of one-dimensional (1D) edge detection [25,28–34] (described in Appendix A) is applied for performing a bilateral edge detection. To this end, multiple horizontal grayscale profiles (with an equal vertical interval of 1 pixel) are scanned along the axial direction (i.e., the Y-direction) with their search directions being bilaterally outward from Subsequently, edge detection for the captured side-viewed image is performed, as shown in Figure 10a or Figure 10b; a rectangular region of interest (ROI), numbered ROI #1, is assigned to appropriately cover the shadows of the two objects appearing in the FOV. The concept of one-dimensional (1D) edge detection [25,28–34] (described in Appendix A) is applied for performing a bilateral edge detection. To this end, multiple horizontal grayscale profiles (with an equal vertical interval of 1 pixel) are scanned along the axial direction (i.e., the Y-direction) with their search directions being bilaterally outward from the middle of ROI #1. According to the detected locations where sudden changes of grayscale values

occur, two groups of edge points (i.e., the left and right edges indicated in the figures) with a Y-directional sub-pixel resolution can thus be obtained. Furthermore, as shown in Figure 10c or Figure 10d, the rightmost edge point among the left group (denoted by G) and the leftmost edge point among the right group (denoted by D) are extracted to calculate their horizontal distance, *dp*. The visually measured distance *dp* is regarded as the positioning distance between the drilltip (or a certain ground cross-section) and the GWEF at that moment. figures) with a Y-directional sub-pixel resolution can thus be obtained. Furthermore, as shown in Figure 10c or Figure 10d, the rightmost edge point among the left group (denoted by G) and the leftmost edge point among the right group (denoted by D) are extracted to calculate their horizontal distance, dp. The visually measured distance dp is regarded as the positioning distance between the drilltip (or a certain ground cross-section) and the GWEF at that moment.

the middle of ROI #1. According to the detected locations where sudden changes of grayscale values occur, two groups of edge points (i.e., the left and right edges indicated in the

Machines 2023, 11, x FOR PEER REVIEW 14 of 34

Figure 10. Illustrative examples of the image processing procedure for the DP sub-process: (a) edge detection for a side-viewed image; (b) edge detection for another side-viewed image; (c) measuring the positioning distance between the drilltip and the GWEF; (d) measuring the positioning distance between a certain ground cross-section and the GWEF. **Figure 10.** Illustrative examples of the image processing procedure for the DP sub-process: (**a**) edge detection for a side-viewed image; (**b**) edge detection for another side-viewed image; (**c**) measuring the positioning distance between the drilltip and the GWEF; (**d**) measuring the positioning distance between a certain ground cross-section and the GWEF.

The criterion for stopping this sub-process is dp being less than or equal to a prescribed positioning distance (the condition of dp ≤ 500 µm (≈1547.98 pixels in a side-viewed image) was given). If not, the cycle of moving the microdrill sample forward a small axial distance, capturing and storing a side-viewed image, and performing the bilateral edge detection is repeated until the given criterion can be achieved. The actual positioning dis-The criterion for stopping this sub-process is *dp* being less than or equal to a prescribed positioning distance (the condition of *dp* ≤ 500 µm (≈1547.98 pixels in a side-viewed image) was given). If not, the cycle of moving the microdrill sample forward a small axial distance, capturing and storing a side-viewed image, and performing the bilateral edge detection is repeated until the given criterion can be achieved. The actual positioning distance *dp* obtained in the final cycle is stored for further usage (in the GO sub-process).

#### tance dp obtained in the final cycle is stored for further usage (in the GO sub-process). *3.2. Grind-Off (GO) Sub-Process*

3.2. Grind-Off (GO) Sub-Process This sub-process is similar to that adopted by Chang et al. [28]. Two situations corresponding to the counting index i 1 (for grinding off an original microdrill) and i > 1 (for regrinding off a ground microdrill) both exist in this sub-process and their schematic dia-This sub-process is similar to that adopted by Chang et al. [28]. Two situations corresponding to the counting index *i* = 1 (for grinding off an original microdrill) and *i* > 1 (for regrinding off a ground microdrill) both exist in this sub-process and their schematic diagrams are depicted in Figure 11a,b, respectively.

grams are depicted in Figure 11a,b, respectively.

 Figure 11. Schematic diagrams of the GO sub-process: (a) grinding off an original microdrill; (b) **Figure 11.**  For the situation of i 1, the axial distance between the cross-section of the original microdrill sample to be ground off and the GWEF is (dp + lc(1)), as shown in Figure 11a; the actual positioning distance dp has been obtained in the DP sub-process, and lc(1) is the specified axial distance from the drilltip to the reference point O of the 1st specified axial position. The microdrill sample is moved a prescribed transverse distance along the negative X-direction (via the X-axis LMT) and is then moved an axial distance of (dp + lc(1)) along the positive Y-direction (via the Y-axis LMT), so that the drill body to be ground off can be moved to a specific location with a small transverse distance (of several millimeters) from the inner annular portion of the grinding wheel, where the reference point O can be aligned with the GWEF. Subsequently, the microdrill sample is fed transversely along the For the situation of *i* = 1, the axial distance between the cross-section of the original microdrill sample to be ground off and the GWEF is (*dp* + *lc*(1)), as shown in Figure 11a; the actual positioning distance *dp* has been obtained in the DP sub-process, and *lc*(1) is the specified axial distance from the drilltip to the reference point O of the 1st specified axial position. The microdrill sample is moved a prescribed transverse distance along the negative X-direction (via the X-axis LMT) and is then moved an axial distance of (*dp* + *lc*(1)) along the positive Y-direction (via the Y-axis LMT), so that the drill body to be ground off can be moved to a specific location with a small transverse distance (of several millimeters) from the inner annular portion of the grinding wheel, where the reference point O can be aligned with the GWEF. Subsequently, the microdrill sample is fed transversely along the positive X-direction (via the X-axis LMT) to perform the grinding off operation, that is, to roughly grind off the drill body and finish the ground cross-section.

positive X-direction (via the X-axis LMT) to perform the grinding off operation, that is, to roughly grind off the drill body and finish the ground cross-section. For the situation of i > 1, the axial distance between the cross-section of the ground microdrill sample to be reground off and the GWEF is (dp + c(i)), as shown in Figure 11b; the actual positioning distance dp has been obtained in the DP sub-process, and c(i) is an For the situation of *i* > 1, the axial distance between the cross-section of the ground microdrill sample to be reground off and the GWEF is (*dp* + ∆*c*(*i*) ), as shown in Figure 11b; the actual positioning distance *dp* has been obtained in the DP sub-process, and ∆*c*(*i*) is an increment amount calculated by

$$
\Delta\_{\mathcal{L}(i)} = l\_{\mathcal{L}(i)} - l\_{\mathcal{L}(i-1)} \tag{1}
$$

() = () − (ିଵ) (1) which is the axial distance from the previously ground cross-section to the reference point O of the ith specified axial position. Similar to the situation of i 1, the microdrill sample is moved a prescribed transverse distance along the negtive X-direction (via the X-axis LMT) and is then moved an axial distance of (dp + c(1)) along the positive Y-direction (via the Y-axis LMT), so that the remaining drill body to be reground off can also be moved to a specific location with a small transverse distance (of several millimeters) from the inner annular portion of the grinding wheel, where the reference point O can be aligned with the GWEF. Subsequently, the microdrill sample is fed transversely along the positive Xdirection (via the X-axis LMT) to perform the grinding off operation, i.e., to roughly which is the axial distance from the previously ground cross-section to the reference point O of the *i*th specified axial position. Similar to the situation of *i* = 1, the microdrill sample is moved a prescribed transverse distance along the negtive X-direction (via the X-axis LMT) and is then moved an axial distance of (*dp* + ∆*c*(1)) along the positive Y-direction (via the Y-axis LMT), so that the remaining drill body to be reground off can also be moved to a specific location with a small transverse distance (of several millimeters) from the inner annular portion of the grinding wheel, where the reference point O can be aligned with the GWEF. Subsequently, the microdrill sample is fed transversely along the positive X-direction (via the X-axis LMT) to perform the grinding off operation, i.e., to roughly regrind off the remaining drill body and finish the reground cross-section.

regrind off the remaining drill body and finish the reground cross-section. The feed rate for performing the grinding off operations was set to 0.1 mm/s. When the microdrill sample has been fed to a specific position leaving a small transverse distance

increment amount calculated by

(of several millimeters) from the outer boundary of the GWEF, this sub-process is stopped. The ground microdrill sample is then transversely moved to the front of the telecentric lens (during Step 8), as shown by the ends of the moving trajectories of point O shown in Figure 11.

#### *3.3. Auto-Focusing (AF) Sub-Process*

This sub-process is based on the application of two-dimensional discrete Fourier transform (2D-DFT) [29,31–33]. When the microdrill sample has been moved to the second prescribed position (during Step 8), an axial-viewed image of its ground cross-section is then captured and stored as a temporary file. For a digital grayscale image *g*(*x*, *y*) consisting of an *N* × *M* matrix of grayscale values, its 2D-DFT can be computed by [29,31–33]

$$E\_F(\mu, v) = \frac{1}{MN} \sum\_{x=0}^{M-1} \sum\_{y=0}^{N-1} g(x, y) e^{-j2\pi(\frac{\mu x}{M} + \frac{vy}{N})} \tag{2}$$

in which *u* and *v* are the 2D frequency variables and *j* = √ −1. For the captured axialviewed image (with *N* = 1944 and *M* = 2592), its Fourier spectrum can thus be computed and depicted.

In Figure 12a, four axial-viewed images captured from an unfocused condition to a clearly focused condition are presented, while their corresponding centered (shifted) Fourier spectra are depicted in Figure 12b. (Note that the four Fourier spectra are all depicted from *u* = 0 to *u* = 0.35*M* and from *v* = 0 to *v* = 0.35*N*, and the origin (i.e., *u* = 0 and *v* = 0) is located at the center of each 2D spectral plot.) As can be observed, higher spectrum intensities (with brighter colors) are concentrated in the low-frequency domains of the 2D spectral plots, especially the clearly focused one. Therefore, a weighted sum of squares of the spectrum intensities in the low-frequency domain (from *u* = 0 to *u* = 0.25*M* and from *v* = 0 to *v* = 0.25*N*), denoted by *Swl*, can be obtained by

$$\mathcal{S}\_{wl} = \sum\_{\substack{\mu=0 \\ \mu \neq 0}}^{(M/4)-1} \sum\_{v=0}^{(N/4)-1} \left(\mu^2 + v^2\right) \left| E\_F(\mu, v) \right|^2 \tag{3}$$

which can be used to evaluate the focusing condition. By specifying a focusing search range according to the reference position of the Y-axis LMT being increased from Y = Y*f*(b) to Y = Y*f*(f), the value of the weighted sum of squares, *Swl*, is a function of the Y-axial position, i.e., *Swl* = *Swl*(Y). Thereby, a relative focusing index, Γ, can defined as

$$\Gamma(\mathbf{Y}) = \frac{\mathbf{S}\_{\text{wl}}(\mathbf{Y})}{\mathbf{S}\_{\text{wl}(\text{max})}} \text{ for } \mathbf{Y}\_{f(\mathbf{b})} \le \mathbf{Y} \le \mathbf{Y}\_{f(\mathbf{f})} \tag{4}$$

in which *Swl*(max) is the global maximum value of *Swl*(Y) within the specified focusing search range. Figure 13 shows an illustrative example of the variation curve of the relative focusing index Γ with respect to the Y-axial relative position (within a range of 880 µm), for which 89 axial-viewed images were captured sequentially with a positional increment of 10 µm. As can be observed, when the condition of Γ(Y) = 1 (i.e., *Swl*(Y) = *Swl*(max)) is achieved, the axial-viewed image captured at that Y-axial position can be regarded as an optimally focused one. Nevertheless, obtaining the variation curve of Γ(Y) for each ground cross-section is quite time-consuming.

Figure 12. Illustration of the AF sub-process: (a) axial-viewed images captured from an unfocused condition to a clearly focused condition; (b) the corresponding centered Fourier spectra. **Figure 12.** Illustration of the AF sub-process: (**a**) axial-viewed images captured from an unfocused condition to a clearly focused condition; (**b**) the corresponding centered Fourier spectra.

In order to efficiently search for the optimal focusing position within the specified

focusing search range, the bisection method [35] is applied in this sub-process. To this end,

positions (within the search range) are captured and stored, and their yielded values of Swl are calculated and ordered (denoted by Swl(BS,min), Swl(BS,mid), and Swl(BS,max) in ascending order). For evaluating the focusing condition, a relative difference ratio, rf, can be defined

as

In order to efficiently search for the optimal focusing position within the specified focusing search range, the bisection method [35] is applied in this sub-process. To this end, three axial-viewed images corresponding to the backmost, middle, and frontmost Y-axial positions (within the search range) are captured and stored, and their yielded values of *Swl* are calculated and ordered (denoted by *Swl*(BS,min), *Swl*(BS,mid), and *Swl*(BS,max) in ascending order). For evaluating the focusing condition, a relative difference ratio, *rf* , can be defined as Machines 2023, 11, x FOR PEER REVIEW 18 of 34 ௪(BS,max) − ௪(BS,mid)

$$r\_f = \frac{\mathcal{S}\_{\text{wl}(\text{BS,max})} - \mathcal{S}\_{\text{wl}(\text{BS,mid})}}{\mathcal{S}\_{\text{wl}(\text{BS,max})}} \tag{5}$$

(5)

The criterion for stopping this sub-process is *rf* being less than or equal to a prescribed small value (the condition of *rf* ≤ 0.05 was given). If not, the search range is updated (according to the range between the two Y-axial positions where *Swl*(BS,mid) and *Swl*(BS,max) occur, respectively), and the cycle of capturing and storing the three axial-viewed images within the updated search range and then calculating and ordering their yielded values of *Swl* is repeated until the given criterion can be achieved. Among the three axial-viewed images captured in the final cycle, the one yielding a value of *Swl*(BS,max) can meet the condition of Γ(Y) ≈ 1 (i.e., *Swl*(BS,max) ≈ *Swl*(max)) and is stored as an optimally focused one for further usage (in the WTM and ODM sub-processes). small value (the condition of rf ≤ 0.05 was given). If not, the search range is updated (according to the range between the two Y-axial positions where Swl(BS,mid) and Swl(BS,max) occur, respectively), and the cycle of capturing and storing the three axial-viewed images within the updated search range and then calculating and ordering their yielded values of Swl is repeated until the given criterion can be achieved. Among the three axial-viewed images captured in the final cycle, the one yielding a value of Swl(BS,max) can meet the condition of Γ(Y) ≈ 1 (i.e., Swl(BS,max) ≈ Swl(max)) and is stored as an optimally focused one for further usage (in the WTM and ODM sub-processes).

Figure 13. Illustrative example of the variation curve of the relative focusing index Γ with respect to the Y-axial relative position. **Figure 13.** Illustrative example of the variation curve of the relative focusing index Γ with respect to the Y-axial relative position.

#### 3.4. Web Thickness Measuring (WTM) Sub-Process *3.4. Web Thickness Measuring (WTM) Sub-Process*

This sub-process is improved from that adopted by Chang et al. [28]. An illustrative example of this sub-process is presented in Figure 14. As the sequential image processing operations shown in Figure 14a–k are basically identical to those used by Chang et al. [28], their concepts are described briefly. The sequential operations shown in Figure 14l–o regard a specific approach called the best-fit circle (BFC) approach, whose fundamental has This sub-process is improved from that adopted by Chang et al. [28]. An illustrative example of this sub-process is presented in Figure 14. As the sequential image processing operations shown in Figure 14a–k are basically identical to those used by Chang et al. [28], their concepts are described briefly. The sequential operations shown in Figure 14l–o regard a specific approach called the best-fit circle (BFC) approach, whose fundamental has been preliminarily studied by the authors [34] and is explained in this section.

been preliminarily studied by the authors [34] and is explained in this section.

Figure 14a shows an original axial-viewed image obtained at the end of Step 9, and a global coordinate system o-xy is consistently set on the image. This original grayscale image is enhanced via a grayscale transformation function [28] for adjusting its brightness, contrast, and gamma values, so that the ground cross-section can be highlighted relative

temporary file and is transformed into a binary image (consisting of pixels with 1-bit binary values of zeros (black) and ones (white)) via a thresholding approach [29,30]. The transformed binary image, shown in Figure 14c, is dealt with by specific morphological operations [29–31] to fill holes (black pixels) within the cross-section and to eliminate thin protrusions and isolated particles (white pixels) around the cross-section, as shown by the ൝ () () ()

ൡ =

can thus be determined by [38]

⎣ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎡ ଶ

ూೖ

ୀଵ

ూೖ

ୀଵ

ୀଵ

 ూೖ

() = ඩ(

ూೖ

ୀଵ

 ూೖ

 ూೖ

ୀଵ

ୀଵ

ୀଵ

The center and diameter of the fiĴed circle CF(k), denoted by OF(k) and DF(k), respectively,

ൠ = ቊ

ଶ + ()

The root-sum-square (RSS) error of the LS circle-fiĴing for the kth iteration is defined as

 1 ూೖ

−()/2 −()/2

ଶ + () + () + ())

⎦ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎤ ିଵ

⎩ ⎪ ⎪ ⎪ ⎨

⎪ ⎪ ⎪

⎧− (

− (

ూೖ

ୀଵ

ూೖ

ୀଵ

ూೖ

− (

ଶ + ଶ )

ଶ + ଶ )

⎪ ⎪ ⎪ ⎬

⎪ ⎪ ⎪ ⎫

(8)

(11)

ଶ + ଶ )

ୀଵ ⎭

ቋ (9)

ଶ − 4() (10)

ଶ

ూೖ

ୀଵ

ూೖ

ୀଵ

ୀଵ

() = ൜

() = ට()

ଶ +

 ଶ

 ూೖ

> ూ() ూ()

Figure 14. Illustration of the sequential operations of the WTM sub-process. CSWT. As a result, the visually measured CSWT w can be obtained by **Figure 14.** Illustration of the sequential operations of the WTM sub-process.

Figure 14l shows the fiĴed circle CF(1) for the 1st iteration. Obviously, the major parts of theoretical points on contours ФP1 and ФP2 (i.e., all points of {F}1) are enclosed by circle CF(1). That is, circle CF(1) is too large to be a representative common tangent circle for measuring

performing the next circle-fiĴing iteration. For the kth iteration, the group of theoretical

Thereby, the updated group {F}2 is within circle CF(1), as shown in Figure 14m. Subsequently, an updated fiĴed circle CF(2) can be obtained by using Equations (8)–(10) with all points of {F}2 being involved, as shown in Figure 14n. It can be observed that circle CF(2) is smaller than circle CF(1). Expectedly, the following circle-fiĴing iterations would lead to the fiĴed circles being gradually closer to a representative common tangent circle for measuring the CSWT. The iterative procedure is terminated if one of the following conditions is satisfied: (i) the RSS error εF(k) is less than 0.5 pixels (or 0.1 µm) or (ii) the number of all points of {F}k is less than 4. When the iterative procedure is terminated at the kth iteration, as shown in Figure 14o, the fiĴed circle CF(k) is regarded as the resultant best-fit circle (denoted by Ct(BF)), i.e., the representative common tangent circle for measuring the

ଶ +

ଶ − (()

ଶ /4) > 0ൟ (12)

{}ାଵ = {} ∖ ൛(, )ห

points {F}k is thus updated by

Figure 14a shows an original axial-viewed image obtained at the end of Step 9, and a global coordinate system o-xy is consistently set on the image. This original grayscale image is enhanced via a grayscale transformation function [28] for adjusting its brightness, contrast, and gamma values, so that the ground cross-section can be highlighted relative to the background. The enhanced grayscale image, shown in Figure 14b, is stored as a temporary file and is transformed into a binary image (consisting of pixels with 1-bit binary values of zeros (black) and ones (white)) via a thresholding approach [29,30]. The transformed binary image, shown in Figure 14c, is dealt with by specific morphological operations [29–31] to fill holes (black pixels) within the cross-section and to eliminate thin protrusions and isolated particles (white pixels) around the cross-section, as shown by the modified binary image shown in Figure 14d. The exterior contour of the cross-section (formed by a group of white pixel points) in this modified binary image (denoted by Φ) is extracted via a boundary extraction operation [31], as shown in Figure 14e, and the centroid (i.e., the average position of all pixel points) of contour Φ (denoted by O*c*) is calculated, as indicated in Figure 14f. As shown in Figure 14g, a certain pixel point E1 on contour Φ is extracted as a reference point by finding a pixel point having the shortest Euclidean distance among those between all pixel points on contour Φ and centroid O*c*. Thereby, another reference point E2 can be determined by finding a point satisfying the condition of **O***c***E**1 + **O***c***E**2 = **0** (see Figure 14h), which may not exactly coincide with (but is quite close to) contour Φ. Points E1 and E2 are regarded as two reference points for the pair of flute contours of the cross-section (as shown by contour ΦF of Figure 2). Subsequently, two rectangular ROIs centered at points E1 and E2 with their orientations dependent on vector ±**E**1**E**2, numbered ROI #2 and ROI #3, respectively, are assigned to appropriately cover parts of the pair of flute contours, as shown by the operation shown in Figure 14i. In addition, the length of each ROI (vertical to vector ±**E**1**E**2) can be 80 to 95 percent of the Euclidean distance between points E1 and E2, and the width of each ROI (parallel to vector ±**E**1**E**2) is approximately 80 percent of that Euclidean distance. The determined ROI #2 and ROI #3 are then assigned to the enhanced grayscale image (see Figure 14j) for performing unidirectional edge detections within the two ROIs, as shown by the operation shown in Figure 14k. The inward direction parallel to vector ±**E**1**E**2 is defined as the unitary search direction (USD) of each ROI, and the concept of 1D edge detection [25,28–34] (described in Appendix A) is applied for the unidirectional edge detection. To this end, multiple USD-parallel grayscale profiles (with an equal interval of 1 pixel vertical to the USD) are scanned along the USD of each ROI. Two groups (sets) of detected edge points on the pair of flute contours (denoted by {**P**1} and {**P**2}, respectively) with a sub-pixel resolution (along their USDs) can be obtained. Furthermore, two planar curves fitted from {**P**1} and {**P**2}, respectively, can be obtained, as shown by the two fitted flute contours ΦP1 and ΦP2 shown in Figure 14l, while their parametric vector equations can be expressed as

$$\boldsymbol{\Phi}\_{\rm P1} = \boldsymbol{\Phi}\_{\rm P1}(\boldsymbol{w}\_{1}) = \begin{Bmatrix} \boldsymbol{x}\_{\boldsymbol{\Phi}\_{\rm P1}}(\boldsymbol{w}\_{1})\\ \boldsymbol{y}\_{\boldsymbol{\Phi}\_{\rm P1}}(\boldsymbol{w}\_{1}) \end{Bmatrix} = \begin{Bmatrix} \boldsymbol{a}\_{\rm P1x}\boldsymbol{w}\_{1}^{3} + \boldsymbol{b}\_{\rm P1x}\boldsymbol{w}\_{1}^{2} + \boldsymbol{c}\_{\rm P1x}\boldsymbol{w}\_{1} + \boldsymbol{d}\_{\rm P1x} \\ \boldsymbol{a}\_{\rm P1y}\boldsymbol{w}\_{1}^{3} + \boldsymbol{b}\_{\rm P1y}\boldsymbol{w}\_{1}^{2} + \boldsymbol{c}\_{\rm P1y}\boldsymbol{w}\_{1} + \boldsymbol{d}\_{\rm P1y} \end{Bmatrix} \tag{6}$$

$$\Phi\_{\rm P2} = \Phi\_{\rm P2}(w\_2) = \begin{Bmatrix} x\_{\Phi\_{\rm P2}}(w\_2) \\ y\_{\Phi\_{\rm P2}}(w\_2) \end{Bmatrix} = \begin{Bmatrix} a\_{\rm P2x}w\_2^3 + b\_{\rm P2x}w\_2^2 + c\_{\rm P2x}w\_2 + d\_{\rm P2x} \\ a\_{\rm P2y}w\_2^3 + b\_{\rm P2y}w\_2^2 + c\_{\rm P2y}w\_2 + d\_{\rm P2y} \end{Bmatrix} \tag{7}$$

in which *w*1 and *w*2 are two independent variables regarding the order of points for each fitted contour, and the x- and y-components of each fitted contour are represented as cubic polynomials of *w*1 or *w*2 with their coefficients being solved by using the least-squares (LS) polynomial-fitting approach [35].

The BFC approach, illustrated in Figure 14l–o, is based on the application of the LS circle-fitting approach [36–38] in combination with a specific iterative procedure [34]. An initial group (set) of theoretical points on both fitted contours (denoted by {**F**}*k* = {**Φ**P1, **Φ**P2}*k* for *k* = 1) can be calculated via Equations (6) and (7). The coordinate of each theoretical point can be represented as **F***m* = {*x*F*m y*F*m*} T for *m* = 1, 2, . . . , *n*F*k* , and *n*F*k* is the number of

all points of {**F**}*k* . The LS circle-fitting approach is repeatedly executed with all points of {**F**}*k* being involved and updated in every iteration. For the *k*th iteration, three coefficients *a*F(*k*) , *b*F(*k*) , and *c*F(*k*) of the equation of a fitted circle *C*F(*k*) can be solved by [38]

 *a*F(*k*) *b*F(*k*) *c*F(*k*) = *n*F*k* ∑ *m*=1 *x* 2 F*m n*F*k* ∑ *m*=1 *x*F*my*F*m n*F*k* ∑ *m*=1 *x*F*m n*F*k* ∑ *m*=1 *x*F*my*F*m n*F*k* ∑ *m*=1 *y* 2 F*m n*F*k* ∑ *m*=1 *y*F*m n*F*k* ∑ *m*=1 *x*F*m n*F*k* ∑ *m*=1 *y*F*m n*F*k* ∑ *m*=1 1 −1 − *n*F*k* ∑ *m*=1 *x*F*m x* 2 F*m* + *y* 2 F*m* − *n*F*k* ∑ *m*=1 *y*F*m x* 2 F*m* + *y* 2 F*m* − *n*F*k* ∑ *m*=1 *x* 2 F*m* + *y* 2 F*m* (8)

The center and diameter of the fitted circle *C*F(*k*) , denoted by OF(*k*) and *D*F(*k*) , respectively, can thus be determined by [38]

$$\mathbf{O}\_{\mathbf{F}(k)} = \begin{Bmatrix} \mathbf{x}\_{\mathrm{O\!F}(k)} \\ \mathbf{y}\_{\mathrm{O\!F}(k)} \end{Bmatrix} = \begin{Bmatrix} -a\_{\mathrm{F}(k)}/2 \\ -b\_{\mathrm{F}(k)}/2 \end{Bmatrix} \tag{9}$$

$$D\_{\rm F(k)} = \sqrt{a\_{\rm F(k)}^2 + b\_{\rm F(k)}^2 - 4c\_{\rm F(k)}} \tag{10}$$

The root-sum-square (RSS) error of the LS circle-fitting for the *k*th iteration is defined as

$$\varepsilon\_{\rm F(k)} = \sqrt{\sum\_{m=1}^{n\_{\rm Fk}} (\mathbf{x}\_{\rm Fm}^2 + \mathbf{y}\_{\rm Fm}^2 + a\_{\rm F(k)} \mathbf{x}\_{\rm Fm} + b\_{\rm F(k)} \mathbf{y}\_{\rm Fm} + c\_{\rm F(k)})^2} \tag{11}$$

Figure 14l shows the fitted circle *C*F(1) for the 1st iteration. Obviously, the major parts of theoretical points on contours ΦP1 and ΦP2 (i.e., all points of {**F**}1) are enclosed by circle *C*F(1). That is, circle *C*F(1) is too large to be a representative common tangent circle for measuring the CSWT. Therefore, partial points of {**F**}1 that are outside circle *C*F(1) are excluded for performing the next circle-fitting iteration. For the *k*th iteration, the group of theoretical points {**F**}*k* is thus updated by

$$\{\mathbf{F}\}\_{k+1} = \{\mathbf{F}\}\_k \left\{ (\mathbf{x}\_{\text{Fm}}, y\_{\text{Fm}}) \Big| \mathbf{x}\_{\text{Fm}}^2 + y\_{\text{Fm}}^2 - (D\_{\text{F}(k)}^2/4) > 0 \right\} \tag{12}$$

Thereby, the updated group {**F**}2 is within circle *C*F(1), as shown in Figure 14m. Subsequently, an updated fitted circle *C*F(2) can be obtained by using Equations (8)–(10) with all points of {**F**}2 being involved, as shown in Figure 14n. It can be observed that circle *C*F(2) is smaller than circle *C*F(1). Expectedly, the following circle-fitting iterations would lead to the fitted circles being gradually closer to a representative common tangent circle for measuring the CSWT. The iterative procedure is terminated if one of the following conditions is satisfied: (i) the RSS error *ε*F(*k*) is less than 0.5 pixels (or 0.1 µm) or (ii) the number of all points of {**F**}*k* is less than 4. When the iterative procedure is terminated at the *k*th iteration, as shown in Figure 14o, the fitted circle *C*F(*k*) is regarded as the resultant best-fit circle (denoted by *Ct*(BF)), i.e., the representative common tangent circle for measuring the CSWT. As a result, the visually measured CSWT *w* can be obtained by

$$w = D\_{t(\text{BF})} = D\_{\text{F}(k)} \tag{13}$$

in which *Dt*(BF) is the diameter of the best-fit circle *Ct*(BF). In addition, the center of the best-fit circle *Ct*(BF), denoted by O*t*(BF), is equal to OF(*k*) . When the data of obtained *Dt*(BF) and **O***t*(BF) are stored, this sub-process is stopped. The presented WTM sub-process with the BFC approach can be applied to measure the CSWTs of both ST and UC type microdrills.

#### *3.5. Outer Diameter Measuring (ODM) Sub-Process*

This sub-process is based on sequential image processing operations, similar to those in the WTM sub-process. An illustrative example of this sub-process is presented in Figure 15. As the sequential operations shown in Figure 15a–e are essentially identical to those shown in Figure 14a–e, their concepts will not be described again. The sequential operations shown in Figure 15f–h are also based on the BFC approach, and those shown in Figure 15i–l regard a specific approach called the fitted dividing line (FDL) approach; their fundamentals have been preliminarily studied by the authors [34] and are explained in this section.

According to the sequential operations shown in Figure 15a–e, an original axial-viewed image, obtained at the end of Step 9, is dealt with to obtain the exterior contour of the cross-section, Φ, in its modified binary image. Subsequently, as shown by the operations shown in Figure 15f–h, the BFC approach is again applied. That is, an initial group (set) of pixel points on the exterior contour Φ (denoted by {**Q**}*l* for *l* = 1) is extracted, for which the coordinate of each pixel point can be represented as **Q***m* = {*x*R*m y*R*m*} T for *m* = 1, 2, . . . , *n*Q*l* , and *n*Q*l* is the number of all points of {**Q**}*l* . The LS circle-fitting approach is repeatedly executed with all points of {**Q**}*l* being involved and updated in every iteration. For the *l*th iteration, three coefficients *a*Q(*l*) , *b*Q(*l*) , and *c*Q(*l*) of the equation of a fitted circle *C*Q(*l*) can be solved by [38]

 *a*Q(*l*) *b*Q(*l*) *c*Q(*l*) = *n*Q*l* ∑ *m*=1 *x* 2 Q*m n*Q*l* ∑ *m*=1 *x*Q*my*Q*m n*Q*l* ∑ *m*=1 *x*Q*m n*Q*l* ∑ *m*=1 *x*Q*my*Q*m n*Q*l* ∑ *m*=1 *y* 2 Q*m n*Q*l* ∑ *m*=1 *y*Q*m n*Q*l* ∑ *m*=1 *x*Q*m n*Q*l* ∑ *m*=1 *y*Q*m n*Q*l* ∑ *m*=1 1 −1 − *n*Q*l* ∑ *m*=1 *x*Q*m x* 2 Q*m* + *y* 2 Q*m* − *n*Q*l* ∑ *m*=1 *y*Q*m x* 2 Q*m* + *y* 2 Q*m* − *n*Q*l* ∑ *m*=1 *x* 2 Q*m* + *y* 2 Q*m* (14)

The center and diameter of the fitted circle *C*Q(*l*) , denoted by OQ(*l*) and *D*Q(*l*) , respectively, can thus be determined by [38]

$$\mathbf{O}\_{\mathbf{Q}(l)} = \begin{Bmatrix} \mathbf{x}\_{\mathbf{O}\_{\mathbf{Q}}(l)} \\ y\_{\mathbf{O}\_{\mathbf{Q}}(l)} \end{Bmatrix} = \begin{Bmatrix} -a\_{\mathbf{Q}(l)}/2 \\ -b\_{\mathbf{Q}(l)}/2 \end{Bmatrix} \tag{15}$$

$$D\_{\mathbf{Q}(l)} = \sqrt{a\_{\mathbf{Q}(l)}^2 + b\_{\mathbf{Q}(l)}^2 - 4c\_{\mathbf{Q}(l)}} \tag{16}$$

Additionally, the RSS error of the LS circle-fitting for the *l*th iteration is defined as

$$\varepsilon\_{\mathbf{Q}(l)} = \sqrt{\sum\_{m=1}^{n\_{\mathbf{Q}l}} (x\_{\mathbf{Q}m}^2 + y\_{\mathbf{Q}m}^2 + a\_{\mathbf{Q}(l)}x\_{\mathbf{Q}m} + b\_{\mathbf{Q}(l)}y\_{\mathbf{Q}m} + c\_{\mathbf{Q}(l)})^2} \tag{17}$$

The fitted circle *C*Q(1) for the 1st iteration is illustrated in Figure 15f. Obviously, the major parts of pixel points on contour Φ (i.e., all points of {**Q**}1) are enclosed by circle *C*Q(1). In other words, circle *C*Q(1) is too small to overlap the pair of margin contours of the crosssection (as shown by contour ΦM of Figure 2). Therefore, partial points of {**Q**}1 that are inside circle *C*Q(1) are excluded for performing the next circle-fitting iteration. For the *l*th iteration, the group of pixel points {**Q**}*l* is thus updated by

$$\{\mathbf{Q}\}\_{l+1} = \{\mathbf{Q}\}\_l \Big/ \left\{ (x\_{\mathbf{Q}m}, y\_{\mathbf{Q}m}) \Big| x\_{\mathbf{Q}m}^2 + y\_{\mathbf{Q}m}^2 - (D\_{\mathbf{Q}(l)}^2/4) < 0 \right\} \tag{18}$$

Thereby, the updated group {**Q**}2 is out of circle *C*Q(1), and an updated fitted circle *C*Q(2) can be accordingly obtained by using Equations (14)–(16) with all points of {**Q**}2 being involved, as illustrated in Figure 15g. Expectedly, the following circle-fitting iterations would lead to the fitted circles being gradually closer to the pair of margin contours. The iterative procedure is terminated if one of the following conditions is satisfied: (i) the RSS error *ε*Q(*l*) is less than 0.5 pixels (or 0.2 µm) or (ii) the number of all points of {**Q**}*l* is less than 15. When

the iterative procedure is terminated at the *l*th iteration, as shown in Figure 15h, the fitted circle *C*Q(*l*) is regarded as the resultant best-fit circle (denoted by *Cb*(BF)) that overlaps the pair of margin contours. The center and diameter of the best-fit circle *Cb*(BF), denoted by O*b*(BF) and *Db*(BF), can also be obtained. ൝ ୕() ୕() ൡ = ⎣ ⎢ ⎢ ⎢ ⎢ ⎢ ୕୕ ୀଵ ୕ ୀଵ ୕ ୀଵ ୕ ్ ୀଵ ୕ ్ ୀଵ 1 ్ ୀଵ ⎦ ⎥ ⎥ ⎥ ⎥ ⎥ ⎩ ⎪ ⎪ ⎪ ⎨ − ୕(୕ ଶ + ୕ ଶ ) ୀଵ − (୕ ଶ + ୕ ଶ ) ్ ୀଵ ⎭ ⎪ ⎪ ⎪ ⎬ (14)

⎥ ⎥ ⎥ ⎤ ିଵ

⎪ ⎪ ⎪

⎧− ୕(୕

్

ୀଵ

్

ଶ + ୕ ଶ )

⎪ ⎪ ⎪ ⎫

 ୕ ్

ୀଵ

్

= ௧(BF) = () (13)

Machines 2023, 11, x FOR PEER REVIEW 22 of 34

3.5. Outer Diameter Measuring (ODM) Sub-Process

this section.

can be solved by [38]

୕()

⎢ ⎢ ⎢ ⎡ ୕ ଶ ్

ୀଵ

్

in which Dt(BF) is the diameter of the best-fit circle Ct(BF). In addition, the center of the bestfit circle Ct(BF), denoted by Ot(BF), is equal to OF(k). When the data of obtained Dt(BF) and Ot(BF) are stored, this sub-process is stopped. The presented WTM sub-process with the BFC approach can be applied to measure the CSWTs of both ST and UC type microdrills.

This sub-process is based on sequential image processing operations, similar to those in the WTM sub-process. An illustrative example of this sub-process is presented in Figure 15. As the sequential operations shown in Figure 15a–e are essentially identical to those shown in Figure 14a–e, their concepts will not be described again. The sequential operations shown in Figure 15f–h are also based on the BFC approach, and those shown in Figure 15i–l regard a specific approach called the fiĴed dividing line (FDL) approach; their fundamentals have been preliminarily studied by the authors [34] and are explained in

According to the sequential operations shown in Figure 15a–e, an original axialviewed image, obtained at the end of Step 9, is dealt with to obtain the exterior contour of the cross-section, Ф, in its modified binary image. Subsequently, as shown by the operations shown in Figure 15f–h, the BFC approach is again applied. That is, an initial group (set) of pixel points on the exterior contour Ф (denoted by {Q}l for l = 1) is extracted, for which the coordinate of each pixel point can be represented as Qm = {xRm yRm}T for m 1, 2, …, nQl, and nQl is the number of all points of {Q}l. The LS circle-fiĴing approach is repeatedly executed with all points of {Q}l being involved and updated in every iteration. For the lth iteration, three coefficients aQ(l), bQ(l), and cQ(l) of the equation of a fiĴed circle CQ(l)

୕୕

ଶ ్

్

ୀଵ

Figure 15. Illustration of the sequential operations of the ODM sub-process. Additionally, the RSS error of the LS circle-fiĴing for the lth iteration is defined as **Figure 15.** Illustration of the sequential operations of the ODM sub-process.

The center and diameter of the fiĴed circle CQ(l), denoted by OQ(l) and DQ(l), respectively,

ൠ = ቊ

ଶ + ୕()

−୕()/2 −୕()/2

ቋ (15)

ଶ − 4୕() (16)

୕() = ൜

୕() = ට୕()

can thus be determined by [38]

The FDL approach, illustrated in Figure 15i–l, is based on the following procedure. Firstly, pixel points on the pair of margin contours are extracted from the exterior contour Φ. To this end, the Euclidean distances between all points of {**Q**}1 and circular center O*b*(BF), denoted by *lm* for *m* = 1, 2, . . . , *n*Q1, are calculated by

$$l\_m = \left\| \mathbf{Q}\_m - \mathbf{O}\_{b\left(\text{BF}\right)} \right\| = \sqrt{\left(\mathbf{x}\_{\text{Qm}} - \mathbf{x}\_{\text{O}\left(\text{BF}\right)}\right)^2 + \left(y\_{\text{Qm}} - y\_{\text{O}\left(\text{BF}\right)}\right)^2} \tag{19}$$

A criterion to judge whether a pixel point is on the pair of margin contours is

$$\left| l\_m - D\_{b(\text{BF})} / 2 \right| \le \varepsilon\_{\text{M}} \tag{20}$$

in which *ε*M is a specified small constant which can be arbitrarily given as 3 pixels (or 0.5 µm). Thus, a pixel point is regarded as a margin contour point if its calculated distance *lm* is quite close to the radius (*Db*(BF)/2) of the best-fit circle *Cb*(BF). As a result, a group (set) of pixel points on the pair of margin contours (denoted by {**M**}) can be extracted from the exterior contour Φ, as shown in Figure 15i; the coordinate of each pixel point can be represented as **M***m* = {*x*M*m y*M*m*} T for *m* = 1, 2, . . . , *n*M, and *n*M is the number of all points of {**M**}. Secondly, by using the LS line-fitting approach [35] (i.e., the well-known linear regression) with all points of {**M**} being involved, a fitted line (denoted by *L*M) can be obtained via solving its slope *a*M and intercept *b*M by

$$\begin{Bmatrix} \begin{matrix} a\_{\mathsf{M}} \\ b\_{\mathsf{M}} \end{matrix} \end{Bmatrix} = \begin{bmatrix} \sum\_{m=1}^{n\_{\mathsf{M}}} \mathbf{x}\_{\mathsf{M}m}^{2} & \sum\_{m=1}^{n\_{\mathsf{M}}} \mathbf{x}\_{\mathsf{M}m} \\ \sum\_{n\_{\mathsf{M}}}^{n\_{\mathsf{M}}} \mathbf{x}\_{\mathsf{M}m} & \sum\_{m=1}^{n\_{\mathsf{M}}} \mathbf{1} \\ \sum\_{m=1}^{n} \mathbf{x}\_{\mathsf{M}m} & \sum\_{m=1}^{n} \mathbf{1} \end{bmatrix}^{-1} \begin{Bmatrix} \sum\_{m=1}^{n\_{\mathsf{M}}} \mathbf{x}\_{\mathsf{M}m} y\_{\mathsf{M}m} \\ \sum\_{m=1}^{n\_{\mathsf{M}}} y\_{\mathsf{M}m} \end{Bmatrix} \tag{21}$$

or

$$a\_{\mathbf{M}} = \frac{n\_{\mathbf{M}} \sum\_{m=1}^{n\_{\mathbf{M}}} \mathbf{x}\_{\mathbf{M}m} y\_{\mathbf{M}m} - \sum\_{m=1}^{n\_{\mathbf{M}}} \mathbf{x}\_{\mathbf{M}m} \sum\_{m=1}^{n\_{\mathbf{M}}} y\_{\mathbf{M}m}}{n\_{\mathbf{M}} \sum\_{m=1}^{n\_{\mathbf{M}}} \mathbf{x}\_{\mathbf{M}m}^2 - \left(\sum\_{m=1}^{n\_{\mathbf{M}}} \mathbf{x}\_{\mathbf{M}m}\right)^2} \tag{22}$$

$$b\_{\rm M} = \frac{\sum\_{m=1}^{\rm{\rm{\rm{M}}}} \mathbf{x}\_{\rm{\rm{Mm}}}^2 \sum\_{m=1}^{\rm{\rm{\rm{M}}}} y\_{\rm{\rm{Mm}}} - \sum\_{m=1}^{\rm{\rm{\rm{M}}}} \mathbf{x}\_{\rm{\rm{Mm}}} \sum\_{m=1}^{\rm{\rm{\rm{M}}}} \mathbf{x}\_{\rm{\rm{Mm}}} y\_{\rm{\rm{Mm}}}}{n\_{\rm{\rm{M}}} \sum\_{m=1}^{\rm{\rm{\rm{\rm{M}}}}} x\_{\rm{\rm{Mm}}}^2 - \left(\sum\_{m=1}^{\rm{\rm{\rm{M}}}} x\_{\rm{\rm{Mm}}}\right)^2} \tag{23}$$

As shown in Figure 15j, the fitted line *L*M intersects with the pair of margin contours. Thirdly, a line perpendicular to the fitted line *L*M and passing through circular center O*b*(BF), denoted by *L*D, can be obtained via calculating its slope *a*D and intercept *b*D by

$$a\_{\rm D} = -1/a\_{\rm M} \tag{24}$$

$$b\_{\rm D} = y\_{\rm O\_{b}(\rm BF)} - a\_{\rm D} \chi\_{\rm O\_{b}(\rm BF)} = y\_{\rm O\_{b}(\rm BF)} + \chi\_{\rm O\_{b}(\rm BF)} / a\_{\rm M} \tag{25}$$

The obtained line *L*D is called a dividing line, which can be used to divide the group of points {**M**} into two sub-groups {**M**1} and {**M**2}, as indicated in Figure 15k. The criteria for judging whether the *m*th pixel point belongs to sub-group {**M**1} or {**M**2} are

$$\mathbf{M}\_m = \begin{cases} \mathbf{x}\_{\mathbf{M}m} \\ \mathbf{y}\_{\mathbf{M}m} \end{cases} \in \{ \mathbf{M}\_1 \} \text{ if } y\_{\mathbf{M}m} - (a\_{\mathbf{D}} \mathbf{x}\_{\mathbf{M}m} + b\_{\mathbf{D}}) > 0 \tag{26}$$

$$\mathbf{M}\_m = \begin{cases} \mathbf{x}\_{\text{Mm}} \\ \mathbf{y}\_{\text{Mm}} \end{cases} \in \{ \mathbf{M}\_2 \} \text{ if } \mathbf{y}\_{\text{Mm}} - (a\_{\text{D}} \mathbf{x}\_{\text{Mm}} + b\_{\text{D}}) < 0 \tag{27}$$

Finally, the two centroids of {**M**1} and {**M**2} (i.e., the two average positions of all pixel points in individual sub-groups, denoted by E3 and E4, respectively) are calculated, as indicated in Figure 15k. Thereby, two rectangular ROIs centered at points E3 and E4 with their orientations being dependent on vector ±**E**3**E**4, numbered ROI #4 and ROI #5, respectively, are assigned to appropriately cover major parts of the pair of margin contours, as shown by the operation shown in Figure 15l. In addition, the length of each ROI (vertical to vector ±**E**3**E**4) is approximately 90 percent of the chord length of the margin contour, which can be estimated by the Euclidean distance between the first and last pixel points of {**M**1} or {**M**2}, and the width of each ROI (parallel to vector ±**E**3**E**4) can be 50 to 90 percent of the estimated chord length.

Furthermore, the determined ROI #4 and ROI #5 are assigned to the enhanced grayscale image (see Figure 15m) for performing unidirectional edge detections within the two ROIs, as shown by the operation shown in Figure 15n. The inward direction parallel to vector ±**E**3**E**4 is defined as the USD of each ROI, and the concept of 1D edge detection [25,28–34] (described in Appendix A) is applied for the unidirectional edge detection. That is, multiple USD-parallel grayscale profiles (with an equal interval of 1 pixel vertical to the USD) are scanned along the USD of each ROI. Two groups (sets) of edge points on the pair of margin contours (denoted by {**S**1} and {**S**2}, respectively) with a sub-pixel resolution (along their USDs) can be obtained. Subsequently, the LS circle-fitting approach [36–38], as shown in Equations (8)–(11) or Equations (14)–(17), is again executed once with all points of {**S**1} and {**S**2} being involved, and a resultant LS circle (denoted by *Cb*(LS)), i.e., the representative minimum bounding circle for measuring the CSOD, can be solved, as shown in Figure 15o. As a result, the visually measured CSOD *D* can be obtained by

$$D = D\_{b(\text{LS})} \tag{28}$$

in which *Db*(LS) is the diameter of the LS circle *Cb*(LS). In addition, the center of the LS circle *Cb*(LS), denoted by O*b*(LS), is also obtained. When the data of obtained *Db*(LS) and **O***b*(LS) are stored, this sub-process is stopped. The presented ODM sub-process with the BFC and FDL approaches can be applied to measure the CSODs of both ST and UC type microdrills.

In addition, after the execution of both the WTM and ODM sub-processes, the concentricity of the representative common tangent and minimum bounding circles, *Ct*(BF) and *Cb*(LS), of a ground cross-section (denoted by ∆*C*) can be evaluated via calculating the eccentricity of circular centers O*t*(BF) and O*b*(LS), that is,

$$\Delta \mathbf{C} = \left\| \mathbf{O}\_{t(\text{BF})} - \mathbf{O}\_{b(\text{LS})} \right\| \tag{29}$$

The amount of concentricity ∆*C* can be a rule of thumb to evaluate the radial symmetry of the ground cross-section.

### **4. Experimental Results and Discussion**

In order to evaluate the actual effectiveness of the developed system, two series of experiments for measuring the CSWTs and CSODs of certain ST and UC type microdrill samples produced by several companies were conducted. The nominal outer diameter (denoted by *Dn*) and flute length (denoted by *lf* ) of each measured ST type sample were 0.30 and 5.5 mm, respectively. The measured UC type samples were of two specifications: (i) *Dn* = 0.25 mm and *lf* = 3.9 mm and (ii) *Dn* = 0.20 mm and *lf* = 3.0 mm. The nominal overall length and shank diameter of each sample were 38.1 and 3.175 mm (i.e., 1.5 and 0.125 inches), respectively.

#### *4.1. Measuring Uncertainty and Accuracy Tests*

This series of experiments aimed at testing the measuring uncertainty and accuracy that could be achieved by the developed system. A ST type microdrill sample (called sample ST1) and two UC type ones (called samples UC1 (with *Dn* = 0.25 mm) and UC2 (with *Dn* = 0.20 mm)) were measured. For each sample, its drill body was ground off to two specified axial positions of *lc* = 0.15*lf* and 0.45*lf* , while its two ground cross-sections were repeatedly focused, captured, and measured 15 times. For performing the first measurement of each ground cross-section, the established AMP was completely executed (from Step 1 to Step 14) once. For performing the remaining 14 measurements of each ground cross-section, the established AMP was repeatedly executed with Step 4 to Step 7 (i.e., the DP and GO sub-processes) being skipped. As each measured sample was manually removed from and put into the microdrill fixture again after each measurement (according to Step 2 and Step 14), 15 axial-viewed images of each ground cross-section appearing in the FOV with different orientations (angular positions) could be focused and captured for measuring the CSGPs. In addition, for comparison purposes, a Sage Vision SG-1022 measuring microscope (in combination with RVS-250 metrology software), that can achieve a resolution of 0.01 µm and an accuracy of ±3 µm, was manually operated to repeatedly measure the CSGPs of each ground cross-section 15 times, as shown in Figure 16. cross-section, the established AMP was repeatedly executed with Step 4 to Step 7 (i.e., the DP and GO sub-processes) being skipped. As each measured sample was manually removed from and put into the microdrill fixture again after each measurement (according to Step 2 and Step 14), 15 axial-viewed images of each ground cross-section appearing in the FOV with different orientations (angular positions) could be focused and captured for measuring the CSGPs. In addition, for comparison purposes, a Sage Vision SG-1022 measuring microscope (in combination with RVS-250 metrology software), that can achieve a resolution of 0.01 µm and an accuracy of ±3 µm, was manually operated to repeatedly measure the CSGPs of each ground cross-section 15 times, as shown in Figure 16.

The amount of concentricity ∆C can be a rule of thumb to evaluate the radial symmetry of

In order to evaluate the actual effectiveness of the developed system, two series of experiments for measuring the CSWTs and CSODs of certain ST and UC type microdrill samples produced by several companies were conducted. The nominal outer diameter (denoted by Dn) and flute length (denoted by lf) of each measured ST type sample were 0.30 and 5.5 mm, respectively. The measured UC type samples were of two specifications: (i) Dn = 0.25 mm and lf = 3.9 mm and (ii) Dn = 0.20 mm and lf = 3.0 mm. The nominal overall length and shank diameter of each sample were 38.1 and 3.175 mm (i.e., 1.5 and 0.125

This series of experiments aimed at testing the measuring uncertainty and accuracy that could be achieved by the developed system. A ST type microdrill sample (called sample ST1) and two UC type ones (called samples UC1 (with Dn = 0.25 mm) and UC2 (with Dn = 0.20 mm)) were measured. For each sample, its drill body was ground off to two specified axial positions of lc 0.15lf and 0.45lf, while its two ground cross-sections were repeatedly focused, captured, and measured 15 times. For performing the first measurement of each ground cross-section, the established AMP was completely executed (from Step 1 to Step 14) once. For performing the remaining 14 measurements of each ground

Machines 2023, 11, x FOR PEER REVIEW 26 of 34

the ground cross-section.

inches), respectively.

4. Experimental Results and Discussion

4.1. Measuring Uncertainty and Accuracy Tests

Figure 16. Measuring the CSGPs of a microdrill sample by a measuring microscope. **Figure 16.** Measuring the CSGPs of a microdrill sample by a measuring microscope.

The statistical results of this series of experiments are listed in Tables 1–3, in which the CSWT, CSOD, and concentricity measured by the developed automation system are The statistical results of this series of experiments are listed in Tables 1–3, in which the CSWT, CSOD, and concentricity measured by the developed automation system are denoted by *w*(A), *D*(A), and ∆*C*(A), respectively, and those measured by the measuring microscope are denoted by *w*(M), *D*(M), and ∆*C*(M), respectively. For each CSGP, the maximum and minimum values, mean value (denoted by *µ*), standard deviation (denoted by *σ*), uncertainty (denoted by *ũ*), and statistical upper and lower bounds (denoted by *β*(*u*) and *β*(*l*) , respectively) evaluated from its 15 measured values are presented. For each CSWT (*w*(A) or *w*(M)) and CSOD (*D*(A) or *D*(M)), the uncertainty *ũ* is evaluated via the well-known 3*σ*-band approach [39], that is,

$$
\ddot{\mu} = \mathbf{3}\sigma \tag{30}
$$

for which the theoretical probability within the band of ±3*σ* is 99.73% for a Gaussian distribution assumption. Thereby, its statistical upper and lower bounds are

∼

$$\mathcal{B}\_{(u)} = \mu + \tilde{u} = \mu + \mathfrak{G}\sigma \tag{31}$$

$$\mathcal{B}\_{(l)} = \mu - \tilde{\mu} = \mu - \mathbf{3}\sigma \tag{32}$$

The uncertainty range of ±*ũ* (or ±3*σ*) can be regarded as the measuring repeatability for each CSWT and CSOD. However, the 3*σ*-band approach is not suitable for evaluating the uncertainty of each concentricity (∆*C*(A) or ∆*C*(M)), since the resultant statistical lower bound may be a negative value that violates the essence of ∆*C* ≥ 0 based on Equation (29). Hence, a modified approach is used to evaluate the uncertainty of each concentricity. To this end, its 15 measured values are further divided into two data groups by comparing

their individual magnitudes with the mean value *µ*, and the mean values and standard deviations of the group with smaller magnitudes (denoted by *µS* and *σS*) and of that with larger magnitudes (denoted by *µL* and *σL*) can be obtained. The condition of 0 ≤ *µS* < *µL* essentially exists. Therefore, the modified statistical upper and lower bounds for each concentricity are estimated by

$$
\beta\_{(u)} = \mu\_L + 1.6449 \sigma\_L \tag{33}
$$

$$
\beta\_{(l)} = \mu\_S - 1.6449 \sigma\_S \tag{34}
$$

for which the theoretical probability within the band of ±1.6449*σ* is 90% for a Gaussian distribution assumption. Accordingly, the modified uncertainty *ũ* for each concentricity, as a rule of thumb, is defined as

$$
\tilde{\mu} = \beta\_{(\mu)} - \beta\_{(l)} = (\mu\_L - \mu\_S) + 1.6449(\sigma\_L + \sigma\_S) \tag{35}
$$

which is a conservatively estimated band based on the Gaussian distributions of the two divided data groups.

As can be observed in Tables 1–3, the representative statistical values (i.e., the results of *µ* ± *ũ*) of *w*(A) were 118.86 ± 0.538 and 169.11 ± 0.547 µm for sample ST1, 120.12 ± 0.619 and 149.86 ± 0.964 µm for sample UC1, and 125.74 ± 0.950 and 141.68 ± 0.701 µm for sample UC2, respectively; those of *w*(M) were 118.21 ± 1.771 and 169.29 ± 1.646 µm for sample ST1, 120.24 ± 2.196 and 149.01 ± 1.959 µm for sample UC1, and 126.02 ± 1.725 and 141.16 ± 1.829 µm for sample UC2, respectively. Thereby, the developed system could achieve a worst uncertainty range of ±1 µm for measuring the CSWTs, while that achieved by the measuring microscope was ±2.2 µm. It can also be found that, for each cross-section, the absolute difference between the mean values of *w*(A) and *w*(M) was significantly less than 1 µm. By calculating the maximum absolute differences between the two measured extremes of *w*(A) and the two of *w*(M) for each cross-section, the obtained results were 2.11 and 1.35 µm for sample ST1, 2.03 and 2.05 µm for sample UC1, and 1.26 and 1.62 µm for sample UC2. That is, for a conservative evaluation, the relative deviation between the measured CSWTs obtained by the developed system and by the measuring microscope could be less than 2.2 µm. As a result, based on the narrow uncertainty range of ±1 µm and the small relative deviation of 2.2 µm, the developed system could achieve good repeatability and accuracy for measuring the CSWTs.

**Table 1.** Statistical results for the measurements of sample ST1 (*Dn* = 0.30 mm and *l f* = 5.5 mm).

**Table 2.** Statistical results for the measurements of sample UC1 (*Dn* = 0.25 mm and *l f* = 3.9 mm).

**Table 3.** Statistical results for the measurements of sample UC2 (*Dn* = 0.20 mm and *l f* = 3.0 mm).

Furthermore, the representative statistical values of *D*(A) were 294.89 ± 0.584 and 291.94 ± 0.525 µm for sample ST1, 222.67 ± 0.456 and 224.52 ± 0.985 µm for sample UC1, and 201.46 ± 0.679 and 186.93 ± 0.646 µm for sample UC2, respectively; those of *D*(M) were 295.09 ± 2.069 and 291.70 ± 1.694 µm for sample ST1, 222.58 ± 1.958 and 224.62 ± 1.641 µm for sample UC1, and 201.88 ± 1.782 and 186.66 ± 1.744 µm for sample UC2, respectively. Thereby, the developed system could achieve a worst uncertainty range of ±1 µm for measuring the CSODs, while that achieved by the measuring microscope was ±2.1 µm. Additionally, for each cross-section, the absolute difference between the mean values of *D*(A) and *D*(M) was considerably less than 0.5 µm. By calculating the maximum absolute differences between the two measured extremes of *D*(A) and the two of *D*(M) for each cross-section, the obtained results were 1.31 and 1.47 µm for sample ST1, 1.71 and 1.80 µm for sample UC1, and 1.48 and 1.36 µm for sample UC2. In other words, the relative deviation between the measured CSODs obtained by the developed system and by the measuring microscope could be conservatively less than 2 µm. As a result, based on the narrow uncertainty range of ±1 µm and the small relative deviation of 2 µm, the developed system could achieve good repeatability and accuracy for measuring the CSODs.

In addition, all of the six mean values of ∆*C*(A) ranged between 1.13 and 2.69 µm, but those of ∆*C*(M) ranged between 6.91 and 10.25 µm; their corresponding absolute differences notably ranged between 4.22 and 8.62 µm. In addition, all of the six estimated uncertainty bands of ∆*C*(A) ranged between 0.770 and 1.481 µm, however, those of ∆*C*(M) ranged between 14.995 and 19.304 µm. The unreasonably large absolute differences (greater than 4 µm) and uncertainty bands of ∆*C*(M) (greater than 14.5 µm) implied that the measuring microscope was not an adequate device to measure the amount of concentricity. Since the measuring microscope was manually operated to find the representative common tangent and minimum bounding circles of a ground cross-section, the positional stabilities of their

circular centers sensitively depended on the skill and dexterity of the operator, although the uncertainty ranges of their diameters could be nonsensitively within ±2.2 µm. As compared with the used measuring microscope, an estimated uncertainty band less than 1.5 µm indicated that the developed system could achieve a narrower uncertainty band and a more confident measured value for the amount of concentricity.

In summary, from a conservative viewpoint, the developed system could achieve an uncertainty range of ±1.5 µm and an accuracy of 2.5 µm for the measurements of the CSWTs and CSODs of both ST and UC type microdrills. The developed system could also achieve a conservatively estimated uncertainty band of 2 µm, as well as a confident measured value, for the measurement of the concentricity, which could not be achieved by the used measuring microscope.

## *4.2. Applicability Tests for Sampling Inspection*

This series of experiments aimed at evaluating the applicability of the developed system for automatically executing the sampling inspection procedure. An ST type microdrill sample (called sample ST2) and two UC type ones (called samples UC3 (with *Dn* = 0.25 mm) and UC4 (with *Dn* = 0.20 mm)) were measured. For each sample, five axial positions of *lc* = 0.15*lf* , 0.3*lf* , 0.45*lf* , 0.6*lf* , and 0.75*lf* were specified, and their corresponding cross-sections were sequentially ground and measured by executing the established AMP once; the measuring microscope (see Figure 16) was also employed to measure each ground cross-section. To this end, when a sample had been ground off to a specified axial position and measured, the AMP was immediately paused (at the end of Step 12), and the sample was temporarily removed from the microdrill fixture and then manually measured once by the use of the measuring microscope. Subsequently, the sample was thereupon put back into and held by the microdrill fixture, and the paused AMP was resumed for the remaining measuring cycles. Such operations were repeated until all the five axial positions of the drill body had been ground and measured. In addition, the required time for sequentially executing the DP and GO sub-processes (from Step 4 to Step 8) once was approximately 2.5 to 2.75 min, and that for sequentially executing the AF, WTM, and ODM sub-processes (from Step 9 to Step 12) once was approximately 1.5 to 1.75 min. Therefore, the developed system took approximately 4 to 4.5 min to destructively and visually measure the CSGPs at each specified axial position of a microdrill sample.

The statistical results of this series of experiments are listed in Tables 4–6. Based on the uncertainty and accuracy testing results presented in Section 4.1, the values of ∆*C*(M) measured by the measuring microscope are not considered. For all the three samples, the 15 absolute differences between the measured values of *w*(A) and *w*(M) ranged between 0.37 and 1.83 µm, while those between the measured values of *D*(A) and *D*(M) ranged between 0.35 and 1.74 µm. In other words, for each cross-section, the developed system could measure its CSWT and CSOD with a good accuracy of less than 2 µm. The measured values of ∆*C*(A) also reasonably ranged between 0.80 and 3.37 µm. Therefore, the established AMP could be executed successfully by the developed system.

**Table 4.** Statistical results for the measurements of sample ST2 (*Dn* = 0.30 mm and *l f* = 5.5 mm).

a Fitted linear regression model: *w*(A) ≈ 159.700(*lc*/*lf* ) + 95.666 µm (with a norm of residuals of 2.913 µm).

**Table 5.** Statistical results for the measurements of sample UC3 (*Dn* = 0.25 mm and *l f* = 3.9 mm).

a Fitted linear regression model: *w*(A) ≈ 80.490(*lc*/*lf* ) + 108.490 µm (with a norm of residuals of 0.749 µm).

**Table 6.** Statistical results for the measurements of sample UC4 (*Dn* = 0.20 mm and *l f* = 3.0 mm).

a Fitted second-order regression model: *w*(A) ≈ −72.724(*lc*/*lf* ) 2 + 113.384(*lc*/*lf* ) + 100.330 µm (with a norm of residuals of 2.761 µm).

In addition, the measured values of *w*(A) of each sample obviously appeared a progressively increasing trend with the increased axial position, as shown in Figure 17, which can be fitted and formulated via some regression models. By using the LS polynomial-fitting approach [35], a first-order (linear) regression model for a set of measured values of *w*(A) can be represented as [24,25,28]

$$w\_{\rm (A)} \approx a\_1 \left(l\_c / l\_f\right) + a\_0 \tag{36}$$

and a second-order regression model for it can be represented as

$$w\_{\rm (A)} \approx a\_2 \left(l\_c / l\_f\right)^2 + a\_1 \left(l\_c / l\_f\right) + a\_0 \tag{37}$$

As can be observed in Figure 17, the measured values of *w*(A) of samples ST2 and UC3 both had a linearly increasing trend and are fitted by using Equation (36), while those of sample UC4 had a nonlinearly increasing trend and are fitted by using Equation (37). The resultant regression models are listed in the last rows (footnotes) of Tables 4–6, and the norm of residuals of each regression model was significantly less than 3 µm. The well-fitted regression models implied that the helical flutes of these samples were properly machined according to the web taper design of microdrills [5,23–25,28]. The measured values and fitted regression models of the CSWTs could be further used to compare with their corresponding theoretical designs, which should benefit the practical quality control tasks.

Accordingly, the applicability of the developed system for automating the sampling inspection procedure could be verified by the presented experimental results. Therefore, the actual effectiveness of the developed system for carrying out the destructive and visual measurements of the CSWTs and CSODs of both ST and UC type microdrills could be validated.

Figure 17. Measured values and fiĴed curves of the CSWTs of certain microdrill samples. **Figure 17.** Measured values and fitted curves of the CSWTs of certain microdrill samples.

#### 5. Conclusions **5. Conclusions**

validated.

An automation system improved from the one proposed by Chang et al. [28], for carrying out the destructive and visual measurements of the CSGPs of both ST and UC type microdrills, has been presented in this paper. The primary conclusions are as follows: An automation system improved from the one proposed by Chang et al. [28], for carrying out the destructive and visual measurements of the CSGPs of both ST and UC type microdrills, has been presented in this paper. The primary conclusions are as follows:

UC4 had a nonlinearly increasing trend and are fiĴed by using Equation (37). The resultant regression models are listed in the last rows (footnotes) of Tables 4–6, and the norm of residuals of each regression model was significantly less than 3 µm. The well-fiĴed regression models implied that the helical flutes of these samples were properly machined according to the web taper design of microdrills [5,23–25,28]. The measured values and fiĴed regression models of the CSWTs could be further used to compare with their correspond-

Accordingly, the applicability of the developed system for automating the sampling inspection procedure could be verified by the presented experimental results. Therefore, the actual effectiveness of the developed system for carrying out the destructive and visual measurements of the CSWTs and CSODs of both ST and UC type microdrills could be

ing theoretical designs, which should benefit the practical quality control tasks.

Accordingly, the actual effectiveness of the developed system had been validated, since it could actually provide good repeatability and accuracy for carrying out the destructive and visual measurements of the CSGPs of both ST and UC type microdrills. Therefore, the developed system could effectively and comprehensively automate the conventional sampling inspection procedure.

**Author Contributions:** Conceptualization, W.-T.C.; methodology, W.-T.C.; software, W.-T.C. and Y.-Y.L.; validation, W.-T.C. and Y.-Y.L.; formal analysis, W.-T.C.; investigation, W.-T.C. and Y.-Y.L.; resources, W.-T.C.; data curation, W.-T.C. and Y.-Y.L.; writing—original draft preparation, W.-T.C.; visualization, W.-T.C. and Y.-Y.L.; funding acquisition, W.-T.C. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research is funded by National Science and Technology Council of Taiwan under Grant Nos. 101-2218-E-019-003, 110-2221-E-019-059, and 111-2221-E-019-073, and is also funded by National Taiwan Ocean University under Grant Nos. NTOU-RD-AA-2012-104011 and NTOU-RD-AA-2013-104041.

**Data Availability Statement:** Not applicable.

**Conflicts of Interest:** The authors declare no conflict of interest.

### **Appendix A**

The fundamental of the 1D edge detection used in the DP, WTM, and ODM subprocesses, which is based on the concept of 1D grayscale profile scanning [29–33] and has also been adopted by several studies [25,28,34], is described here. By setting a local coordinate system o*r*-x*r*y*r* on an assigned rectangular ROI with its x*r*-axis parallel to one of the two pairs of opposite sides of the ROI (i.e., parallel to the search direction) and its origin o*r* coinciding with a corner of the ROI, an arbitrary 1D grayscale profile within the ROI, denoted by *gp*, can be represented as

$$\mathbf{g}\_p = \mathbf{g}\_p(\mathbf{x}\_{r\prime} y\_{r(\mathbf{S})}) \tag{A1}$$

in which *xr* is a positional variable (that is, *xr* = 0, 1, . . . , *nr*x (in unit of pixels) with *nr*x being the total number of pixels of the x*r*-directional side of the ROI), and *yr*(S) is a specified y*r*-directional position (a constant in unit of pixels). The criterion applied to evaluate whether an edge point (*xr*(E), *yr*(S)) exists is given as follows:

$$\max\_{\mathbf{x}\_{r(\mathcal{E})} \in (0, n\_{\rm rx})} \Delta \mathbf{g}\_p \left( \mathbf{x}\_{r(\mathcal{E})}, y\_{r(\mathcal{S})} \right) \ge \mathbf{C}\_{\rm E} \tag{A2}$$

in which *C*E is a given contrast grayscale value [29], and the criterion value ∆*gp* is defined as

$$\Delta \mathbf{g}\_p \left( \mathbf{x}\_{r(\mathbf{E})}, \mathcal{Y}\_{r(\mathbf{S})} \right) = \frac{\begin{vmatrix} \sum\_{\mathbf{r}\_r = \mathbf{r}\_r(\mathbf{E})}^{\mathbf{r}\_r(\mathbf{E}) - \mathbf{S}\_\mathbf{E}/2} & \mathbf{g}\_p \left( \mathbf{x}\_r \mathbf{y}\_{r(\mathbf{S})} \right) - \sum\_{\mathbf{r}\_r = \mathbf{r}\_r(\mathbf{E}) + \mathbf{S}\_\mathbf{E}/2}^{\mathbf{r}\_r(\mathbf{E}) + (\mathbf{W}\_\mathbf{E} + \mathbf{S}\_\mathbf{E}/2)} \mathbf{g}\_p \left( \mathbf{x}\_r \mathbf{y}\_{r(\mathbf{S})} \right) \right)}{\mathbf{W}\_\mathbf{E}} \\ \text{for} \begin{cases} \begin{aligned} \mathbf{W}\_\mathbf{E} + \mathbf{S}\_\mathbf{E}/2 \le \mathbf{x}\_{r(\mathbf{E})} \le n\_{r\mathbf{x}} - \left( \mathbf{W}\_\mathbf{E} + \mathbf{S}\_\mathbf{E}/2 \right) \\ 2 \mathbf{W}\_\mathbf{E} + \mathbf{S}\_\mathbf{E} < n\_{r\mathbf{x}} \\ 0 < \mathbf{S}\_\mathbf{E} \le \mathbf{W}\_\mathbf{E} \end{cases} \end{cases} \end{cases} \tag{A3}$$

where *W*E and *S*E are given filter width and steepness values (both in pixels) [29], respectively. As referred to Equations (A2) and (A3), an edge point would be located at (*xr*(E), *yr*(S)) when the absolute difference between the two average grayscale values calculated from *xr* = *xr*(E) − (*W*E + *S*E/2) to *xr* = *xr*(E) − *S*E/2 and from *xr* = *xr*(E) + *S*E/2 to *xr* = *xr*(E) + (*W*E + *S*E/2), respectively, is greater than or equal to the given contrast grayscale value *C*E. The location of (*xr*(E), *yr*(S)) can be determined via Equation (A2) to find a local maximum of ∆*gp*, and a fitted 1D grayscale profile, denoted by *gˆp*, can be further obtained as

$$\mathfrak{g}\_p = \mathfrak{g}\_p\left(\mathfrak{x}\_{r\prime} y\_{r(\mathbb{S})}\right) \text{ for } \mathfrak{x}\_{r(\mathbb{E})} - \mathfrak{S}\_{\mathbb{E}}/2 \le \mathfrak{x}\_r \le \mathfrak{x}\_{r(\mathbb{E})} + \mathfrak{S}\_{\mathbb{E}}/2\tag{A4}$$

To this end, curve-fitting approaches [35,40,41] with cubic polynomials or spline curves can be applied. Subsequently, a refined edge point *x* ∗ *r*(E) , *yr*(S) can be determined by numerically solving the following equation:

$$\left(\mathfrak{x}\_{r(\mathbb{E})'}^\* y\_{r(\mathbb{S})}\right) - \mathfrak{C}\_{\mathbb{E}} = \mathbf{0} \tag{A5}$$

The refined edge point *x* ∗ *r*(E) , *yr*(S) can achieve a sub-pixel resolution along the search direction (i.e., the x*r*-direction). As a result, multiple refined edge points with an equal y*r*-directional interval of 1 pixel (within the ROI) can all be detected.

Finally, each refined edge point is represented in a global coordinate system o-xy (that is consistently set on the captured digital image) via the well-known homogeneous transformation matrix [42], that is,

$$
\begin{Bmatrix} \mathbf{x}\_{(\mathrm{E})} \\ \mathbf{y}\_{(\mathrm{E})} \\ 1 \end{Bmatrix} = \begin{bmatrix} \cos \Delta \theta & -\sin \Delta \theta & \Delta x \\ \sin \Delta \theta & \cos \Delta \theta & \Delta y \\ 0 & 0 & 1 \end{bmatrix} \begin{Bmatrix} \mathbf{x}\_{r(\mathrm{E})}^{\*} \\ \mathbf{y}\_{r(\mathrm{S})} \\ 1 \end{Bmatrix} \tag{A6}$$

in which ∆*x* and ∆*y* are the x- and y-components of a vector from origin O to origin O*r* , respectively, and ∆*θ* is the angle of the positive x*r*-axis counterclockwise measured from the positive *x*-axis. The transformed edge coordinate (*x*(E), *y*(E)) is the refined edge point represented in the global coordinate system o-xy. Multiple transformed edge coordinates can therefore be obtained and stored for further usage.

### **References**

**Disclaimer/Publisher's Note:** The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

## *Article* **Improving Industrial Robot Positioning Accuracy to the Microscale Using Machine Learning Method**

**Vytautas Bucinskas , Andrius Dzedzickis \* , Marius Sumanas , Ernestas Sutinys, Sigitas Petkevicius, Jurate Butkiene, Darius Virzonis and Inga Morkvenaite-Vilkonciene**

> Department of Mechatronics, Robotics and Digital Manufacturing, Vilnius Gediminas Technical University, 10223 Vilnius, Lithuania

**\*** Correspondence<EMAIL_ADDRESS>
**Abstract:** Positioning accuracy in robotics is a key issue for the manufacturing process. One of the possible ways to achieve high accuracy is the implementation of machine learning (ML), which allows robots to learn from their own practical experience and find the best way to perform the prescribed operation. Usually, accuracy improvement methods cover the generation of a positioning error map for the whole robot workspace, providing corresponding correction models. However, most practical cases require extremely high positioning accuracy only at a few essential points on the trajectory. This paper provides a methodology for the online deep Q-learning-based approach intended to increase positioning accuracy at key points by analyzing experimentally predetermined robot properties and their impact on overall accuracy. Using the KUKA-YouBot robot as a test system, we perform accuracy measurement experiments in the following three axes: (i) after a long operational break, (ii) using different loads, and (iii) at different speeds. To use this data for ML, the relationships between the robot's operating time from switching on, load, and positioning accuracy are defined. In addition, the gripper vibrations are evaluated when the robot arm moves at various speeds in vertical and horizontal planes. It is found that the robot's degrees of freedom (DOFs) clearances are significantly influenced by operational heat, which affects its static and dynamic accuracy. Implementation of the proposed ML-based compensation method resulted in a positioning error decrease at the trajectory key points by more than 30%.

**Keywords:** online machine learning; deep-q learning; positioning accuracy; industrial robot; vibrations

## **1. Introduction**

Robots are used in various processes, including manufacturing, entertainment, services, and scientific research. To maintain a technical edge and thereby remain competitive, more and more businesses are applying advanced technology and programming solutions to their operational processes [1]. Such a wide application encourages the development of universal robotic systems and requires research of their capabilities and performance characteristics.

In general, industrial robots provide high-level static and dynamic positioning accuracy. Nevertheless, they must be maintained to ensure that they continue to meet the conditions for which they have been programmed and in which they operate [2]. Therefore, for each specific task, it is important to determine the following: (i) positioning accuracy (positioning error between stated and the real position of arm end effector); (ii) repeatability (positioning error between real positions of arm end effector performing repeating movements); (iii) other parameters, which are considered to be unique characteristics of the particular machine [3]. Positioning accuracy depends on a number of actions including, but not restricted to the following: (i) the parameters of the drives guiding the robot's movements; (ii) tolerances in manufacturing parts of the machinery; (iii) tolerances due to the articulation of the robot's chains [4]; (iv) control algorithms [5,6]; (v) dynamic properties of the mobile platform [7,8], and (vi) robot arm properties [9,10]. Each of the factors may

**Citation:** Bucinskas, V.; Dzedzickis, A.; Sumanas, M.; Sutinys, E.; Petkevicius, S.; Butkiene, J.; Virzonis, D.; Morkvenaite-Vilkonciene, I. Improving Industrial Robot Positioning Accuracy to the Microscale Using Machine Learning Method. *Machines* **2022**, *10*, 940. https://doi.org/10.3390/ machines10100940

Academic Editors: Raul D.S.G. Campilho, Francisco J. G. Silva and Dan Zhang

Received: 11 August 2022 Accepted: 13 October 2022 Published: 17 October 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

become important for the accuracy of the robot, depending on its type, lifting capacity, and operational conditions [4]. Due to the complex nature of robots, the increasing positioning accuracy of robotic installations develops challenging tasks. The positioning errors in offline mode can be determined using special algorithms [11] and mathematical models [12–15]. Common positioning accuracy problems in path and trajectory planning can be resolved by improving offline programming software when the code for a robot is generated automatically or by implementing online feedback [16–20]. However, typical compensation methods are limited by the chaotic nature of the positioning errors. Therefore, more adaptive methods, such as machine learning (ML), can increase positioning accuracy by compensating for positioning errors in particular cases [21,22]. Advances in ML with respect to simple error compensation are in the accountancy of feature-chaotic robot error distribution about positioning points. ML ensures respect for robot fluctuation and fits into the prescribed tolerance field. There are a few cases [23–26] of ML approaches compensating for the positioning errors. Chen and Zhang developed an ML kinematicsbased positioning error compensation method for high-precision mechanical machining operations [26]. This positioning error compensation method provided good results, but it contains sophisticated procedures due to its complexity and use of excessive data. More to say, it combines the following: analytical modeling; extended Kalman filtering; spatial interpolation algorithm; an adaptive mesh division algorithm, and an inverse distance weighted interpolation algorithm. To use it at a few trajectory points, this method becomes too costly.

Moreover, ML could be implemented to improve the trajectory-planning process [27–29] or reduce vibrations [30,31].

The main aim of our research is to create a methodology to increase industrial robot positioning accuracy and minimize robot end-of-arm vibrations by applying ML in online mode. This methodology will be used to compensate for positioning errors by shifting the destination point or altering the moving velocity of the end-of-arm reference point before defining the forward kinematic task. The method is focused on destination point approach accuracy in defined arbitrary robot positions. Moreover, our proposed approach suits well for industrial robots with typically closed controllers since position correction is performed externally with respect to the robot controller; thus, modified motion commands to target coordinates are processed in a standard way.

#### **2. Concept of Research**

This paper focused on creating a universal methodology for most types of industrial robots, which can increase positioning accuracy at the robot trajectory endpoints. One of the pillars of ML is mathematical optimization, which involves the numerical computation of parameters for a system designed to make decisions based on unseen data [32–34]. While the ML procedure is based on existing data collection, such a method all the time remains retrospective and reflects errors in previous applications with corresponding load cases. The use of the robot for precise machining or assembly operations requires data about the vibration level and settling time. This will be critically important in addition to static positioning error values and directions. The ML procedure enables a decrease in positioning error values, and end-of-arm vibration influences the resulting accuracy. Nevertheless, the robot workspace mapping procedure, which provides error values and direction at any point of the workspace, all the time remains inaccurate for each particular case.

It is possible to calculate optimal parameters for a given learning problem using currently available data [35]. The collection of data required for ML is a complex process, defined by the robot design, its characteristics, and the aim of ML implementation.

The procedure of implementing our proposed ML method is divided into the following three phases: (i) initial preparation, (ii) positioning task formulation, and (iii) optimization (Figure 1).

zation (Figure 1).

**Figure 1.** Simplified block diagram of positioning accuracy increase procedure by implementing ML method. RF—reference point. **Figure 1.** Simplified block diagram of positioning accuracy increase procedure by implementing ML method. RF—reference point.

It is possible to calculate optimal parameters for a given learning problem using currently available data [35]. The collection of data required for ML is a complex process, defined by the robot design, its characteristics, and the aim of ML implementation.

The procedure of implementing our proposed ML method is divided into the following three phases: (i) initial preparation, (ii) positioning task formulation, and (iii) optimi-

The success of the ML procedure depends on the quantity and quality of data generated during initial preparation. This phase includes experiments for the definition of the main robot's characteristics, such as positioning accuracy, reference point vibration level, statistical evaluation, and analysis of obtained data. The initial preparation procedure should be performed for all robot trajectory points of interest. This procedure must be repeated in cases of mechanical wear, change of tool configuration, or essential variation of environmental conditions. The success of the ML procedure depends on the quantity and quality of data generated during initial preparation. This phase includes experiments for the definition of the main robot's characteristics, such as positioning accuracy, reference point vibration level, statistical evaluation, and analysis of obtained data. The initial preparation procedure should be performed for all robot trajectory points of interest. This procedure must be repeated in cases of mechanical wear, change of tool configuration, or essential variation of environmental conditions.

The obtained data shows the technical conditions of the robot since these data consist of a set of dependencies between positioning errors, vibration level, reference point position, operating time, robot speed, and load. In the general case, these dependencies can be expressed as follows: The obtained data shows the technical conditions of the robot since these data consist of a set of dependencies between positioning errors, vibration level, reference point position, operating time, robot speed, and load. In the general case, these dependencies can be expressed as follows:

$$\begin{cases} \Delta\_{\mathbf{x},\mathbf{y},\mathbf{z}} = F(\mathbf{x},\mathbf{y},\mathbf{z},\mathbf{V},\mathbf{M},\mathbf{t}\_{op})\\ \varepsilon\_{\mathbf{x},\mathbf{y},\mathbf{z}} = F(\mathbf{x},\mathbf{y},\mathbf{z},\mathbf{V},\mathbf{M}) \end{cases} \tag{1}$$

where: ∆*x*,*y*,*z* is the overall positioning error, *εx*,*y*,*z* is the overall level of reference point vibrations, *x y z* are coordinates of robot reference point, *V* is the speed of reference point, *M* is load, *top* is the time from the moment when a robot was switched on.

The second phase of our proposed ML method is the Formulation of the positioning task. In this phase, it is necessary to define the required position and speed of the reference point (*x y z*, *V*) and load (*M*). Moreover, it is necessary to specify the following limitations: acceptable overall positioning error ∆*min*, and overall vibration level of reference point *εmin*. In general, case limitations can be identified as follows:

$$\begin{cases} \Delta\_{\mathbf{x},\mathbf{y},\mathbf{z}} \le \Delta\_{\text{min}}\\ \varepsilon\_{\mathbf{x},\mathbf{y},\mathbf{z}} \le \varepsilon\_{\text{min}} \end{cases} \tag{2}$$

Moreover, to ensure stable operation of the ML algorithm and avoid looping in the algorithm (when the desired accuracy cannot be achieved), it is necessary to include a parameter specifying the maximum number of iterations *kmax*.

The third phase is the Optimization procedure, which contains the cycle during which the optimization algorithm is used to obtain the most suitable compensation parameters in order to achieve the required positioning accuracy and an acceptable level of reference point vibrations.

The ML algorithm runs with each new positioning task; ML selects compensation parameters and predicts the expected positioning error ∆*x,y,z,* and reference point vibration level. This prediction is based on experimentally defined data obtained during the Initial preparation phase. After defining the expected positioning error ∆*x,y,z* and reference point vibration level *εx,y,z*, their values are compared with the limitations. If positioning error and reference point vibration level correspond to the limitations of Equation (2), then modified coordinates of the required reference point position and speed are transferred to the robot controller. In case if results do not correspond to the limitations, the data is saved and checked if the maximum number of iterations *kmax* are reached. The optimization of compensation parameters is performed many times until the limit of iterations is reached. In case if the maximum number of iterations is not achieved, the program selects the best result from all previous iterations.

## *Deep q-Learning Algorithm*

The deep Q-learning algorithm was implemented to realize the mentioned ML procedure. The advantage of this algorithm in comparison to others is emphasized due to its efficiency in similar problems and its relatively simple implementation [36–38]. The deep Q-learning algorithm combines the Q-learning algorithm and a deep artificial neural network. The idea of Q-learning is based on the perception of the environment and state to take respective actions in order to achieve the maximum reward. A neural network enables the algorithm to operate in a much larger environment and optimize calculations procedure by enabling approximation features. It allows the algorithm to observe the pattern in the environment and discover the optimal sequences of actions instead of calculating and evaluating each state and the value of each action at each point in the environment. The same principle was used to optimize the definition of the algorithm loss function. We used the stochastic gradient descent method (SGD) [35] to define the function by the theoretical gradient calculated from randomly selected data points instead of calculating the actual gradient from the entire dataset. Such an approach diminishes the dependency of obtained data from strong dataset influence on functional distribution. Other parameters of the implemented algorithm were selected experimentally by evaluating their impact on positioning accuracy after many trials in a separate study. The configuration of the algorithm used in this study is provided in Table 1.

**Table 1.** Parameters of the implemented algorithm.

The parameters within Table 1 were defined through the experimental running of the algorithm. Process activation through the Softsign procedure corresponds to the robot learning mode, while there are possibilities for other functions. The activation margin (sensitivity to the activation condition) needs redefinition; the ML process converges

correctly and, within 1000 iterations, achieves learning process influence saturation. The implemented optimizer—SGD—corresponds to our purposes and achieves the prescribed result. A number of neurons is chosen according to the desired process parameter resolution, and the number of hidden layers is naturally chosen as one according to the dimensionality of the ML process. The amount of replay memory is defined by the experimental test, and it occurs outside our paper scope as well as the temperature parameter. **Table 1.** Parameters of the implemented algorithm. **Parameter Activation Function Optimizer Hidden Neurons Hidden Layers Replay Memory Temperature**  Value Softsign SGD 37 1 100,000 70 The parameters within Table 1 were defined through the experimental running of the

*Machines* **2022**, *10*, x FOR PEER REVIEW 5 of 21

Furthermore, this paper provides a detailed methodology for experimental research and data collection for ML. algorithm. Process activation through the Softsign procedure corresponds to the robot learning mode, while there are possibilities for other functions. The activation margin

#### **3. Experimental Research** (sensitivity to the activation condition) needs redefinition; the ML process converges cor-

#### *3.1. Experimental Setup* rectly and, within 1000 iterations, achieves learning process influence saturation. The im-

The experiments were performed using the KUKA-YouBot robot (KUKA, Augsburg, Germany) as an industrial robot testbench; it was fixed to a special stable base. The robot's geometric parameters and main characteristics are provided in Table 2 [39]. Positions of the robot's gripper were detected using two USB cameras with a resolution of 1920 × 1080, a checkerboard matrix of 4 × 8 mm size, and a user-defined function, implementing the detect Checkerboard Points procedure in MatLab (Figure 2). According to [40], the use of such a measurement method can ensure an average measurement deviation better than 0.0004 mm. plemented optimizer—SGD—corresponds to our purposes and achieves the prescribed result. A number of neurons is chosen according to the desired process parameter resolution, and the number of hidden layers is naturally chosen as one according to the dimensionality of the ML process. The amount of replay memory is defined by the experimental test, and it occurs outside our paper scope as well as the temperature parameter. Furthermore, this paper provides a detailed methodology for experimental research and data collection for ML. **3. Experimental Research**

**Table 2.** General characteristics KUKA-YoubBot arm [39]. *3.1. Experimental Setup*

**Figure 2.** Experimental setup for measurements of positioning accuracy: (**a**) general view; (**b**) position of the cameras; (**c**) accelerometers mounted on robot gripper. 1—base; 2—robot; 3—USB cameras; 4—based on which are mounted cameras; 5—personal computer; 6—optical checkerboard matrix; 7—holder for cameras; 8—accelerometers. **Figure 2.** Experimental setup for measurements of positioning accuracy: (**a**) general view; (**b**) position of the cameras; (**c**) accelerometers mounted on robot gripper. 1—base; 2—robot; 3—USB cameras; 4—based on which are mounted cameras; 5—personal computer; 6—optical checkerboard matrix; 7—holder for cameras; 8—accelerometers.

The cameras were fixed by a special holder positioned at a 90◦ angle relative to each other to resolve spatial coordinates of the end-of-arm reference point (Figure 2b). The holder with attached cameras was screwed to the same base as the robot. Two special checkerboard patterns (identification marks for cameras) were attached to the gripper (Figure 2b,c). checkerboard patterns (identification marks for cameras) were attached to the gripper (Figure 2b,c). Absolute vibration accelerations of the end-of-arm reference point were measured

The cameras were fixed by a special holder positioned at a 90° angle relative to each other to resolve spatial coordinates of the end-of-arm reference point (Figure 2b). The holder with attached cameras was screwed to the same base as the robot. Two special

**General Information Axes Motion Range, Speed**

Serial kinematics 5 axes Axis 0 (A0) +/−169°, 90°/s Height 655 mm Axis 1 (A1) +90°/−65°, 90°/s Work envelope 0.513 m3 Axis 2 (A2) +146°/−151°, 90°/s Weight 6.3 kg Axis 3 (A3) +/−102°, 90°/s Payload 0.5 kg Axis 4 (A4) +/−167°, 90°/s

*Machines* **2022**, *10*, x FOR PEER REVIEW 6 of 21

**Table 2.** General characteristics KUKA-YoubBot arm [39].

Absolute vibration accelerations of the end-of-arm reference point were measured using industrial accelerometers Ini 603C01 (PCB Piezotronics, Depew, NY, USA) with a measurement range of 0.5~10,000 Hz and an acceleration limit of 51 g. Accelerometers were attached to the gripper in three perpendicular directions, as shown in Figure 2c. Signals from accelerometers were collected using the data acquisition system USB-4432 (National Instruments, Austin, TX, USA). using industrial accelerometers Ini 603C01 (PCB Piezotronics, Depew, NY, USA) with a measurement range of 0.5~10,000 Hz and an acceleration limit of 51 g. Accelerometers were attached to the gripper in three perpendicular directions, as shown in Figure 2c. Signals from accelerometers were collected using the data acquisition system USB-4432 (National Instruments, Austin, TX, USA). Ambient temperature was measured using a digital thermometer MWF-DT-616CT

Ambient temperature was measured using a digital thermometer MWF-DT-616CT (CEM Corporation, Matthews, NC, USA), with a resolution of 0.1 C. Loads were weighted using Silver crest HG01025 (Silvercrest, Corona, CA, USA) scales, with a resolution of 0.001 kg. For the determination of the required warm-up time, the unloaded robot moves at 50% of the max joint speed (which is 90◦/s for all axes) to the trajectory endpoint (Figure 3), rotating joints by the following angles: ϕ0—90◦ , ϕ1—15.2◦ , ϕ2—45.7◦ , ϕ3—44.2◦ , ϕ4—20◦ , and backward. Then the robot returns to the start position, waits for 3 s, and the image of the checkerboard matrix is captured. The experiment trials took place after long operation breaks of 12, 2, and 1 h. In the case of a 2- and 1-h operation break, the warming procedure was applied before the experiment. Firstly, the robot was warmed up using a 30 min warming program and 2- and 1-h operation breaks afterward, correspondingly. (CEM Corporation, Matthews, NC, USA), with a resolution of 0.1 C. Loads were weighted using Silver crest HG01025 (Silvercrest, Corona, CA, USA) scales, with a resolution of 0.001 kg. For the determination of the required warm-up time, the unloaded robot moves at 50% of the max joint speed (which is 90°/s for all axes) to the trajectory endpoint (Figure 3), rotating joints by the following angles: φ0—90°, φ1—15.2°, φ2—45.7°, φ3—44.2°, φ4— 20°, and backward. Then the robot returns to the start position, waits for 3 s, and the image of the checkerboard matrix is captured. The experiment trials took place after long operation breaks of 12, 2, and 1 h. In the case of a 2- and 1-h operation break, the warming procedure was applied before the experiment. Firstly, the robot was warmed up using a 30 min warming program and 2- and 1-h operation breaks afterward, correspondingly.

**Figure 3.** Position of the end-of-arm reference point: (**a**) start (measuring) position; (**b**) endpoint of the trajectory. **Figure 3.** Position of the end-of-arm reference point: (**a**) start (measuring) position; (**b**) endpoint of the trajectory.

The positioning accuracy measurements were performed in the same conditions, but in this case, experiments were carried out with 10%, 50%, and 100% of the maximum joint speed, without load, and loaded by 0.250 kg and 0.500 kg weights. Each measurement was repeated 20 times. The positioning accuracy measurements were performed in the same conditions, but in this case, experiments were carried out with 10%, 50%, and 100% of the maximum joint speed, without load, and loaded by 0.250 kg and 0.500 kg weights. Each measurement was repeated 20 times.

The proposed methodology for the online ML procedure was tested for 2 different robot configurations, imitating the pick-and-place task. In each individual case, using our online ML method, the procedure can be performed at any chosen position in the robot workspace.

### *3.2. Evaluation of the End-of-Arm Reference Point Vibrations*

The dynamic parameters of the robot—the end-of-arm reference point vibration level and settling time when the robot stops at the desired positions—were evaluated using accelerometers (Figure 2c). There are the following three parameters for measuring vibrations: displacement, velocity, and acceleration. The evaluation of vibration acceleration provides the best sensitivity in higher frequencies and suits well for evaluating impacts caused by bearing damage, abnormal gears, and noise. Moreover, this method is recommended when it is necessary to evaluate forces and stresses acting on or rotating parts. Moreover, measurements of acceleration do not require reference points. Therefore, it is more convenient for industrial robots. Moreover, measurements of acceleration do not require reference points. Therefore, it is more convenient for industrial robots. The measurements were carried out at the following three different speeds of the

The proposed methodology for the online ML procedure was tested for 2 different robot configurations, imitating the pick-and-place task. In each individual case, using our online ML method, the procedure can be performed at any chosen position in the robot

The dynamic parameters of the robot—the end-of-arm reference point vibration level and settling time when the robot stops at the desired positions—were evaluated using accelerometers (Figure 2c). There are the following three parameters for measuring vibrations: displacement, velocity, and acceleration. The evaluation of vibration acceleration provides the best sensitivity in higher frequencies and suits well for evaluating impacts caused by bearing damage, abnormal gears, and noise. Moreover, this method is recommended when it is necessary to evaluate forces and stresses acting on or rotating parts.

*Machines* **2022**, *10*, x FOR PEER REVIEW 7 of 21

*3.2. Evaluation of the End-of-Arm Reference Point Vibrations*

The measurements were carried out at the following three different speeds of the joints of the robot: at 10%, 50%, and 100% of the maximum (90◦/s) speed of the joint, using the same robot movement trajectory as shown in Figure 3. The measurement was repeated three times. joints of the robot: at 10%, 50%, and 100% of the maximum (90°/s) speed of the joint, using the same robot movement trajectory as shown in Figure 3. The measurement was repeated three times.

The separate DOF effect on the robot end-of-arm vibrations was evaluated by rotating each robot joint separately. For this purpose, joints I-V at coordinates ϕ0–ϕ4 were rotated clockwise by a 60◦ angle at maximum (90◦/s) speed from the start position Figure 4. The end-of-arm reference point vibrations were measured after the movement of each joint by 60◦ at a maximum (90◦/s) speed. After 10 s, the joints were rotated back to the start position. Measurements were performed when the robot moved in forward and backward directions, and the procedure was repeated three times. The separate DOF effect on the robot end-of-arm vibrations was evaluated by rotating each robot joint separately. For this purpose, joints I-V at coordinates φ0–φ4 were rotated clockwise by a 60° angle at maximum (90°/s) speed from the start position Figure 4. The end-of-arm reference point vibrations were measured after the movement of each joint by 60° at a maximum (90°/s) speed. After 10 s, the joints were rotated back to the start position. Measurements were performed when the robot moved in forward and backward directions, and the procedure was repeated three times.

workspace.

**Figure 4.** The start position of the robot for the measurements of separate joint impact to robot vi-**Figure 4.** The start position of the robot for the measurements of separate joint impact to robot vibrations.

brations. To define the impact of joint I (coordinate φ0, Figure 5) on robot end-of-arm vibrations, joint I was rotated at 60° from the start position, and the vibrations were measured. The vibrations at the following two start positions were evaluated: (i) fully extended po-To define the impact of joint I (coordinate ϕ0, Figure 5) on robot end-of-arm vibrations, joint I was rotated at 60◦ from the start position, and the vibrations were measured. The vibrations at the following two start positions were evaluated: (i) fully extended position (Figure 5a); (ii) compacted position (Figure 5b). *Machines* **2022**, *10*, x FOR PEER REVIEW 8 of 21

**Figure 5.** Start positions of the robot used for the measurements of gripper vibrations: (**a**) maximum peripheral reach on extended position, (**b**) compacted position. In both cases, the reference point was distanced from the *z*-axis origin by 200 mm. **Figure 5.** Start positions of the robot used for the measurements of gripper vibrations: (**a**) maximum peripheral reach on extended position, (**b**) compacted position. In both cases, the reference point was distanced from the *z*-axis origin by 200 mm.

In a fully extended position robot's own weight fully loads joint I. Other joints remain unloaded and stay in a singularity position. The compacted robot position is represented by the situation when joint I is loaded by half of the own weight of other links. In a fully extended position robot's own weight fully loads joint I. Other joints remain unloaded and stay in a singularity position. The compacted robot position is represented by the situation when joint I is loaded by half of the own weight of other links.

All vibration measurements were performed according to the flowchart provided in Figure 6. All vibration measurements were performed according to the flowchart provided in Figure 6.

Signals from three accelerometers were recorded continuously during robot movement between the endpoints of the trajectory. Later acquired data were processed offline (Figure 7), defining transient processes due to robot stops at the trajectory

**Figure 6.** Flowchart of the vibration measurement process.

endpoints.

was distanced from the *z*-axis origin by 200 mm.

Figure 6.

**Figure 6.** Flowchart of the vibration measurement process. **Figure 6.** Flowchart of the vibration measurement process.

Signals from three accelerometers were recorded continuously during robot movement between the endpoints of the trajectory. Later acquired data were processed offline (Figure 7), defining transient processes due to robot stops at the trajectory Signals from three accelerometers were recorded continuously during robot movement between the endpoints of the trajectory. Later acquired data were processed offline (Figure 7), defining transient processes due to robot stops at the trajectory endpoints. *Machines* **2022**, *10*, x FOR PEER REVIEW 9 of 21

**Figure 5.** Start positions of the robot used for the measurements of gripper vibrations: (**a**) maximum peripheral reach on extended position, (**b**) compacted position. In both cases, the reference point

by the situation when joint I is loaded by half of the own weight of other links.

In a fully extended position robot's own weight fully loads joint I. Other joints remain unloaded and stay in a singularity position. The compacted robot position is represented

All vibration measurements were performed according to the flowchart provided in

**Figure 7.** Definition of vibration parameters from experimental data. **Figure 7.** Definition of vibration parameters from experimental data.

Main parameters (end-of-arm reference point maximum peak-to-peak vibration amplitude, settling time) were extracted from raw data, processed statistically, grouped according to testing conditions, and represented in graphical form. Main parameters (end-of-arm reference point maximum peak-to-peak vibration amplitude, settling time) were extracted from raw data, processed statistically, grouped according to testing conditions, and represented in graphical form.

#### *3.3. Calculations 3.3. Calculations*

number.

here:

Overall positioning accuracy from measurements in separate axes was evaluated according to the methodology provided by ISO 9283:1998 [41]. It defines positioning accuracy *APp* as the deviation between a reference position and the mean of attained positions while approaching reference positions from the same direction as follows: Overall positioning accuracy from measurements in separate axes was evaluated according to the methodology provided by ISO 9283:1998 [41]. It defines positioning accuracy *APp* as the deviation between a reference position and the mean of attained positions while approaching reference positions from the same direction as follows:

$$AP\_p = \sqrt{\left(\overline{x} - x\_c\right)^2 + \left(\overline{y} - y\_c\right)^2 + \left(\overline{z} - z\_c\right)^2},\tag{3}$$

= ඥ(̅)ଶ + (ത)ଶ + (̅)ଶ. (4)

∑

+ 3, (6)

ିଵ , (5)

ିଵ . (7)

where are coordinates of the prescribed end-of-arm reference point position, and ̅, ത, ̅ are mean values of the resulting position. where *xc yc zc* are coordinates of the prescribed end-of-arm reference point position, and *x*, *y*, *z* are mean values of the resulting position.

General robot positioning error consists of accumulated values from all errors, and it is an object of ML-obtained compensation. General error from three-axis measurements for the experiments in 1, 2, and 12 h evaluated using the simplified equation, since target General robot positioning error consists of accumulated values from all errors, and it is an object of ML-obtained compensation. General error from three-axis measurements

∑

Robot positioning repeatability *RPl*, also was evaluated using the methodology proposed in ISO 9283:1998. According to it, *RPl* value is the radius of the sphere that defines the closeness of an attained position after *n* attempts to achieve the same position

ୀଵ ; = ඥ(̅− )ଶ + (ത−)ଶ + (̅− )ଶ; = ට∑ (ೕି )̅మ ೕషభ

Robot positioning repeatability was evaluated using the same experimental data as used for the positioning accuracy measurement; values obtained after the 14th minute of

where: are coordinates of the resulting position, and *n* is the measurement cycle

= ̅

ିଵ , ̅= ଵ

̅= ଵ

from the same direction as follows:

∑

̅= ଵ

the experiment were used.

where ̅ ത ̅ are calculated from experimental data as follows:

∑

ିଵ , ത = ଵ

position defined as = 0*,*  = 0 and = 0 as follows:

for the experiments in 1, 2, and 12 h evaluated using the simplified equation, since target position defined as *xc* = 0, *yc* = 0 and *zc* = 0 as follows:

$$AP\_p = \sqrt{(\mathfrak{X})^2 + (\mathfrak{Y})^2 + (\mathfrak{z})^2}. \tag{4}$$

where *x y z* are calculated from experimental data as follows: