A method for online real-time detection of circular non-coding landmarks

Abstract

The present invention discloses a circular non-coding mark point online real-time detection method, which is based on CPU-GPU heterogeneous computing and includes: initialization of the CPU and GPU, real-time acquisition of a frame of image, calculation of integral graph and integral square graph on the CPU and asynchronous transmission of image data, image binarization on the GPU, marking and analysis of connected domains on the GPU, edge extraction on the GPU, least squares ellipse fitting on the GPU, asynchronous transmission of detection results from the GPU to the CPU, and judgment on whether to stop image acquisition. If so, the entire detection process is terminated; otherwise, real-time image acquisition and real-time image detection are repeated. The present invention combines the advantages of CPU and GPU, places some steps suitable for serial processing on the CPU, places steps suitable for parallel processing on the GPU, improves the detection speed of circular non-coding mark points, realizes real-time online detection of high-resolution images containing a large number of mark points, and can further improve the detection speed.

Description

Technical Field

[0001] The present invention belongs to the technical field of 3D vision technology, and relates to a circular non-coding landmark point detection method, and further relates to an online real-time detection method for circular non-coding landmark points based on CPU-GPU heterogeneous computing. Background Art

[0002] In some high-precision measurement applications of 3D vision technology, such as real-time 3D motion capture and real-time 3D reconstruction, artificial markers are usually pasted on the surface of the object to be measured, and then the 3D coordinate information of the artificial markers is detected and reconstructed to achieve real-time measurement of the object's motion, deformation, trajectory and posture. Therefore, marker detection technology is the basis for realizing high-precision 3D vision technology. Artificial markers come in various forms and can be divided into square markers, cross markers, circular markers, etc. according to their shapes. Among them, circular non-coded markers have the characteristics of affine invariance, translation invariance, and scale invariance, which are easy to detect and identify, and easy to process and manufacture, and are widely used in the field of 3D vision technology.

[0003] Chinese invention patent ZL 200710133269.7 proposes a circular marker detection method. The method first preprocesses the image to obtain the continuous closed boundary area in the edge map, and then uses the continuous closed boundary area as the feature of the marker, and designs a corresponding feature detection method to detect possible markers in the image. The method has a complex process and the algorithm is not robust, so it is not suitable for the recognition and positioning of circular non-encoded markers in real-time 3D vision technology. The journal article "Sub-pixel Positioning of Circular Markers and Its Application" (Infrared and Laser Engineering, 2008, 37 (S1): 47-50) combines Canny edge detection, target recognition, sub-pixel edge detection and other technologies to propose a sub-pixel positioning method for circular markers in complex backgrounds. Although this method can achieve higher detection accuracy, the algorithm runs slowly and cannot achieve the purpose of real-time detection.

[0004] A new circular landmark point recognition and location algorithm is proposed in the journal article "A fast extraction algorithm for circular landmark points" (Optical Technology, 2013, 39(05): 393-397). This method first quickly obtains the rough center and rough edge points of the ellipse, then performs sub-pixel processing on the edge points, removes noise points based on curvature filtering, and then fits the ellipse. According to the distance relationship between the point and the ellipse, a sub-pixel edge point optimization is performed, and the optimized edge points are used to perform least squares fitting of the ellipse to obtain more accurate center coordinates. This algorithm can only achieve fast detection of a few landmark points (12 landmark points) on low-resolution images (800*600), and cannot achieve the goal of real-time detection for high-resolution images (more than one million) containing more landmark points. An automatic, high-precision circular landmark point recognition and location method is proposed in the journal article "Automatic Recognition and Location of Circular Landmark Points in High-Speed ​​Images" (Journal of Surveying and Mapping Science and Technology, 2020, 37(04): 380-387). The method first uses HOG (Histogram of Oriented Gradient) and SVM (Support Vector Machine) to train the existing data to generate a feature model related to the circular marker, and automatically extracts the area of ​​the image containing the circular marker; secondly, the Canny operator and improved Zernike detection are used to extract the sub-pixel edge information of the circular marker; finally, the least squares fitting method is used to obtain the high-precision center position of the circular marker. Since the Canny operator and improved Zernike are used to extract the integer pixel edge and sub-pixel edge of the marker, the entire detection process is time-consuming and cannot achieve online real-time detection.

[0005] In summary, the existing circular non-coded marker detection algorithms are designed and run on the CPU side. The general process is to first use various methods to obtain the integer pixel edge of the marker; secondly, obtain the sub-pixel edge of the marker based on the integer pixel edge interpolation; finally, use the least squares fitting method to obtain the center coordinates of the circular marker. Since multiple steps such as integer pixel edge extraction, sub-pixel edge extraction and ellipse fitting need to be performed in sequence, the entire process is relatively time-consuming and can only achieve fast detection of low-resolution images containing a few markers. For high-resolution images containing more markers, it is difficult to meet the actual application requirements of online real-time detection on the CPU side. Summary of the invention

[0006] In order to overcome the shortcomings of the above-mentioned prior art, the purpose of the present invention is to provide an online real-time detection method for circular non-coded marker points. The method combines the advantages of both CPU and GPU processors, places some steps suitable for serial processing on the CPU side, and places steps suitable for parallel processing on the GPU side. It overcomes the defect of the prior art that the detection speed drops sharply with the increase of image resolution or the number of marker points, greatly improves the detection speed of circular non-coded marker points, and realizes online real-time detection of high-resolution images containing a large number of marker points.

[0007] In order to achieve the above object, the technical solution adopted by the present invention is:

[0008] The first step is to initialize the CPU side. Allocate the required page-locked memory on the CPU side to store the acquired image data and the calculation results on the CPU side.

[0009] The second step is to initialize the GPU. Allocate the required memory area in the GPU global memory to store the image data transmitted from the CPU and various calculation results on the GPU. Declare two CUDA streams, one of which is a transmission stream for asynchronous data transmission; the other is a processing stream for various image processing on the GPU. Declare CUDA events to determine whether the tasks in the stream have been completed. Declare the required global variables.

[0010] The third step is to capture a frame of image in real time. The camera is controlled to capture a frame of circular non-encoded marker image in real time and store it in the locked page memory allocated on the CPU side in the first step.

[0011] Step 4: CPU-side integral graph and integral square graph calculation and image data asynchronous transmission. The integral graph and integral square graph calculation are implemented on the CPU side, and the circular non-coded marker point image collected in the third step is asynchronously transmitted from the CPU side to the GPU side in the transmission stream declared in the second step, and the calculated integral graph and integral square graph are asynchronously transmitted from the CPU side to the GPU side.

[0012] Step 5: Binarization of the GPU-side image: Based on the integral image and integral square image calculated in step 4, the circular non-encoded marker image is binarized in the processing flow declared in step 2.

[0013] Step 6: Connected domain labeling and analysis on the GPU side: In the processing flow declared in step 2, the connected domains of the binary image obtained in step 5 are labeled and analyzed.

[0014] Step 7: GPU-side edge extraction: In the processing flow declared in step 2, the whole pixel edges and sub-pixel edges in the connected domain marked in step 6 are extracted.

[0015] Step 8: GPU-side least squares ellipse fitting: In the processing flow declared in step 2, the sub-pixel edge extracted in step 7 is subjected to least squares ellipse fitting to obtain the center coordinates of the circular non-coding marker point.

[0016] In the ninth step, the detection result is asynchronously transmitted from the GPU to the CPU. In the transmission stream declared in the second step, the center coordinates of the marker points obtained in the eighth step are asynchronously transmitted from the GPU to the CPU to complete the real-time detection of the one-frame image collected in the third step.

[0017] The tenth step is to determine whether to stop image acquisition. If so, the entire detection process is terminated; otherwise, steps 3 to 9 are repeated to acquire and detect images in real time.

[0018] Compared with the prior art, the present invention has the following beneficial effects:

[0019] (1) The present invention combines the advantages of both CPU and GPU processors, places some steps suitable for serial processing on the CPU side, and places steps suitable for parallel processing on the GPU side, which greatly improves the detection speed of circular non-coding markers.

[0020] (2) The method of the present invention overcomes the defect of the prior art that the detection speed drops sharply as the image resolution or the number of landmarks increases, and can realize real-time online detection of high-resolution images containing a large number of landmarks. The algorithm has higher versatility and faster speed.

[0021] (3) During continuous detection, the detection speed can be further improved because a large amount of time occupied by the initialization of the CPU and GPU is saved. BRIEF DESCRIPTION OF THE DRAWINGS

[0022] Figure 1 Flow chart of the specific operation steps of the present invention.

[0023] Figure 2 One frame of image collected for the present invention contains 65 circular non-encoded marker points.

[0024] Figure 3 for Figure 2 Detection results of circular non-encoding landmarks in . DETAILED DESCRIPTION

[0025] The technical solution of the present invention will be described in detail below with reference to the accompanying drawings.

[0026] The present invention proposes an online real-time detection method for circular non-coding landmarks based on CPU-GPU heterogeneous computing. Figure 1 As shown, the following steps are included:

[0027] The first step is CPU initialization: allocate the required page-locked memory on the CPU to store the acquired image data and the calculation results on the CPU.

[0028] In the embodiment of the present invention, three blocks of locked page memory of size M×N are allocated on the CPU side, which are used for subsequent storage of the collected circular non-encoded marker image, integral image and integral square image, respectively. M and N represent the width and height of the image, respectively.

[0029] The second step is GPU initialization: allocate the required memory area in the GPU global memory to store the image data transmitted from the CPU and various calculation results on the GPU. Declare two CUDA streams, one of which is a transmission stream for asynchronous data transmission; the other is a processing stream for various image processing on the GPU. Declare CUDA events to determine whether the tasks in the stream have been completed. Declare the required global variables.

[0030] In an embodiment of the present invention, a memory area of ​​size M×N×1 is allocated in the global memory of the GPU and texture attributes are bound to it for subsequent storage of grayscale images uploaded from the CPU. A memory area of ​​size M×N×1 is allocated in the global memory of the GPU for subsequent storage of binary images. Two memory areas of size M×N×4 are allocated in the global memory of the GPU for subsequent storage of integral images and integral square images of int type. A memory area of ​​size M×N×4×4 is allocated in the global memory of the GPU for subsequent storage of the coordinates of the circumscribed rectangle of the connected domain of int4 type. A memory area of ​​size M×N×8×2×(2rmax+1) is allocated in the global memory of the GPU. 2 / (2rmin+1) 2 A memory area of ​​size M×N×4 is allocated in the GPU global memory for the subsequent storage of double2-type sub-pixel edge coordinates; where rmin and rmax are the minimum radius and maximum radius of the circular non-coding marker points pre-set by the user. A memory area of ​​size M×N×4 is allocated in the GPU global memory for the subsequent storage of the number of edge pixels in the connected domain of type int. A memory area of ​​size M×N×8×5 is allocated in the GPU global memory for the subsequent storage of double-type ellipse parameters; a memory area of ​​size M×N×8×2 is allocated in the GPU global memory for the subsequent storage of double2-type marker center coordinates. Where M and N represent the width and height of the image, respectively.

[0031] Of the two declared CUDA streams, one is a transmission stream, which is used for asynchronous transmission of data between the CPU and GPU. The other is a processing stream, which is used for various image processing on the GPU.

[0032] Two CUDA events are declared. One event is used to query whether a certain transmission task in the transmission stream has been completed; the other event is used to query whether a certain task in the processing stream has been completed.

[0033] Declare an int-type label matrix lable on the GPU side and allocate a M×N×4 memory area for it in the GPU global memory; declare an int-type root-finding matrix root on the GPU side and allocate a M×N×4 memory area for it in the GPU global memory.

[0034] Declare a coefficient matrix M of double type on the GPU side, and allocate a memory area of ​​M×N×8×25 for it in the GPU global memory; declare a matrix b of double type on the GPU side, and allocate a memory area of ​​M×N×8×5 for it in the GPU global memory.

[0035] On the GPU side, declare an int type global variable rectNums to record the number of circumscribed rectangles of the connected domain, and set the initial value to 0. On the GPU side, declare a bool type global variable isChanged to record whether there is a label change operation, and set the initial value to false.

[0036] The third step is to capture a frame of image in real time: control the camera to capture a frame of circular non-coded marker image in real time, and store the image in a locked page memory of size M×N allocated on the CPU side in the first step. It should be noted that the image captured by the camera must be an 8-bit grayscale image. Figure 2 The figure shows a frame of image containing 65 circular non-coded marker points captured in real time by a Basler industrial camera with a resolution of 5 million (2448*2048).

[0037] The fourth step is the calculation of the integral image and the integral square image on the CPU side and the asynchronous transmission of image data: the calculation of the integral image and the integral square image is realized on the CPU side, and the circular non-coded marker point image collected in the third step is asynchronously transmitted from the CPU side to the GPU side in the transmission stream declared in the second step, and the calculated integral image and integral square image are asynchronously transmitted from the CPU side to the GPU side.

[0038] In an embodiment of the present invention, the integral map is calculated on the CPU side and stored in a page-locked memory of size M×N allocated in the second step; at the same time, the circular non-coding marker point image collected in the third step is asynchronously uploaded to the GPU global memory bound to the texture attribute allocated in the second step using the transmission stream declared in the second step. After the integral map is calculated, the integral square map is calculated on the CPU side and stored in a page-locked memory of size M×N allocated in the second step; at the same time, the calculated integral map is asynchronously uploaded to a GPU global memory of size M×N×4 allocated in the second step using the transmission stream declared in the second step. After the integral square map is calculated, the integral square map is asynchronously uploaded to a GPU global memory of size M×N×4 allocated in the second step using the transmission stream declared in the second step. In addition, the event declared in the second step is inserted into the transmission stream to query whether each step of the upload is completed.

[0039] In the present invention, the integral image refers to an image of size M×N, in which the value of the pixel with coordinates (i, j) is equal to the sum of the grayscale values ​​of all the pixels in the rectangular area formed by the two pixels with coordinates (0, 0) and (i, j) in the image collected in the third step. Where i∈[0,M),j∈[0,N).

[0040] The integral square image refers to an image of size M×N, in which the value of the pixel with coordinates (i, j) is equal to the sum of the squares of the grayscale values ​​of all the pixels in the rectangular area formed by the two pixels with coordinates (0, 0) and (i, j) in the image collected in the third step. Where i∈[0,M),j∈[0,N).

[0041] Step 5: GPU-side image binarization: Based on the integral image and integral square image calculated in step 4, the circular non-encoded marker image is binarized in the processing flow declared in step 2.

[0042] In the embodiment of the present invention, in the fifth step, the status of the event inserted in the transmission stream in the fourth step is first queried. If the event has been completed, the circular non-encoded marker point image uploaded to the GPU global memory in the fourth step is binarized in the processing flow declared in the second step. The specific implementation steps of the GPU-side image binarization are as follows:

[0043] 1) Thread grid division. The CUDA thread grid consists of thread blocks, and thread blocks consist of threads. Since the image is two-dimensional data information, a two-dimensional grid and a two-dimensional thread block are allocated. According to CUDA programming rules, the x dimension (blockDim.x) of a two-dimensional thread block is generally set to a multiple of 32. Since there can be a maximum of 1024 threads in a block, the y dimension (blockDim.y) of a two-dimensional thread block is generally set to an even number between 1 and 32. The specific setting of the y dimension (blockDim.y) of the two-dimensional thread block depends on the usage of shared resources in the streaming multiprocessor (SM); if the shared resources in the streaming multiprocessor (SM), such as registers and shared memory, are used more, the y dimension (blockDim.y) of the two-dimensional thread block is set to a smaller even number or 1.

[0044] After the two-dimensional thread block is set, set the two-dimensional grid. The two-dimensional grid calculation formulas are as follows:

[0045]

[0046] Among them, gridDim.x represents the x dimension of the two-dimensional grid, gridDim.y represents the y dimension of the two-dimensional grid; blockDim.x represents the x dimension of the two-dimensional thread block, blockDim.y represents the y dimension of the two-dimensional thread block, and M and N represent the width and height of the image respectively.

[0047] After the thread grid is divided, each thread in the two-dimensional thread grid is responsible for the binarization processing task of a corresponding pixel in the image. The correspondence between the pixels in the image and the threads in the thread grid is as follows:

[0048]

[0049] Among them, i and j represent the coordinate components of the pixel with coordinates (i,j), i∈[0,M), j∈[0,N); threadIdx represents the thread index in the two-dimensional thread block, threadIdx.x and threadIdx.y are the index components of the x dimension and y dimension respectively, threadIdx.x∈[1,blockDim.x], threadIdx.y∈[1,blockDim.y]; blockIdx represents the thread block index in the two-dimensional thread grid, blockIdx.x and blockIdx.y are the index components of the x dimension and y dimension respectively, blockIdx.x∈[1,gridDim.x], threadIdx.y∈[1,gridDim.y].

[0050] 2) Threshold value calculation. A thread in the two-dimensional thread grid divided in step 1) is responsible for calculating the threshold value of a pixel in the image. The specific calculation method is as follows:

[0051]

[0052] Among them, T(i,j) represents the threshold of the pixel with coordinates (i,j); m(i,j) represents the mean grayscale value of the pixels in the square area with the pixel with coordinates (i,j) as the center and a side length of 2r+1, r represents half the side length of the square area, and is generally taken as r∈[10,20]; s(i,j) represents the standard deviation of the grayscale values ​​of the pixels in the square area; k is a constant, and is generally taken as k∈[0.2,0.5].

[0053] The mean m(i,j) is calculated as follows:

[0054] m(i,j)=[I′(i+r,j+r)+I′(ir,jr)-I′(i+r,jr)-I′(ir,j+r)] / (2r+1) 2

[0055] Where I′(i,j) represents the value of the pixel with coordinates (i,j) in the integral image.

[0056] The standard deviation s(i,j) is calculated as follows:

[0057]

[0058] Where I″(i,j) represents the value of the pixel with coordinates (i,j) in the integral square graph.

[0059] 3) Binarization. Based on the threshold calculated in step 2), each pixel is binarized:

[0060]

[0061] Among them, I b Represents a binary image, I b (i, j) represents the value of the pixel with coordinates (i, j) in the binary image; I(i, j) represents the grayscale value of the pixel with coordinates (i, j) in the circular non-encoded marker image uploaded to the GPU global memory in the fourth step.

[0062] Step 6: GPU-side connected domain labeling and analysis: In the processing flow declared in step 2, the connected domain labeling and analysis are performed on the binary image obtained in step 5. The specific implementation steps are as follows:

[0063] 1) Thread grid division. Divide the two-dimensional thread grid according to the method described in step 5. After the grid division is completed, each thread in the thread grid is responsible for the connected domain labeling operation of a corresponding pixel in the image.

[0064] 2) Initialize the label matrix label and the root matrix root. If the binary image I obtained in the fifth step b If the value of the pixel with coordinate (i,j) is 255, the element value label(i,j) at the corresponding position of the label matrix label declared in the second step is assigned to j×N+i; otherwise, the element value label(i,j) at the corresponding position is assigned to -1. The values ​​of all elements of the root matrix root declared in the second step are assigned to -1. Where i∈[0,M), j∈[0,N), M and N represent the width and height of the image respectively.

[0065] 3) Update the root-finding matrix. For the element with coordinates (i, j) in the label matrix, if its value label(i, j) is not -1, then access the nine elements in the 3×3 square area centered on the element, and assign the non-negative minimum value of the nine element values ​​to the element root(i, j) with coordinates (i, j) in the root-finding matrix root.

[0066] 4) Update the label matrix. Perform root-finding operations on all elements in the root-finding matrix root and update the label matrix. The specific implementation steps are: for the element with coordinates (i, j) in the root-finding matrix root, if its value root(i, j) is a non-negative number, calculate the new coordinates according to the element value and access the element in the root-finding matrix root corresponding to the new coordinates. The new coordinates corresponding to the non-negative element value root(i, j) are (root(i, j) % N, root(i, j) / N). If the element value at the new coordinate position in the root-finding matrix root is non-negative, repeat the aforementioned root-finding operation until the calculated new coordinates are the same as the coordinates calculated last time, then end the root-finding operation, assign the element value in the root-finding matrix root corresponding to the new coordinates calculated last time to the element with coordinates (i, j) in the label matrix lable, and assign the Boolean global variable isChanged to true (true), indicating that there is a change in the label during the root-finding operation.

[0067] 5) If the Boolean global variable isChanged in step 4) is true, set the Boolean global variable isChanged to false, and repeat steps 3) and 4). If the Boolean global variable isChanged in step 4) is false, the connected domain is marked.

[0068] 6) Connected domain analysis. Based on CUDA atomic operations, calculate the minimum xmin, ymin and maximum xmax, ymax of the coordinates of each element in the connected domain marked in step 5) in the x and y dimensions; if the conditions are met: xmax-xmin>=2×rmin+1, ymax-ymin>=2×rmin+1, xmax-xmin<=2×rmax+1, ymax-ymin<=2×rmax+1, then store these four values ​​(xmin, ymin, xmax, ymax) in the global memory of type int4 allocated in the second step, and add 1 to the global variable rectNums. Among them, xmin, ymin, xmax, ymax are the coordinates of the circumscribed rectangle of the connected domain.

[0069] Step 7: GPU-side edge extraction: extract the whole pixel edges and sub-pixel edges in the connected domain marked in step 6 in the processing flow declared in step 2. The specific implementation steps are as follows:

[0070] 1) Divide the thread grid. Divide the two-dimensional thread grid according to the method described in step 5. The x-dimension and y-dimension of the two-dimensional thread block are set to blockDim.x=32, blockDim.y=32, and the x-dimension and y-dimension of the two-dimensional grid are set to gridDim.x=rectNums, gridDim.y=1, so that one thread block is responsible for the edge extraction task of one connected domain.

[0071] 2) Integer pixel edge extraction in the connected domain. According to the coordinates of the connected domain circumscribed rectangle obtained in the sixth step, determine the area where integer pixel edge extraction needs to be performed in the label matrix label; for an element in this area with coordinates (i, j) and a non-negative label value, if the label values ​​of its four neighboring elements are all non-negative, then the element with coordinates (i, j) in the root matrix root is assigned a value of -1. After all elements in the connected domain are processed, the coordinates corresponding to the non-negative elements in the root matrix root are the coordinates of the integer pixel edge points. The number of integer pixel edge points in the connected domain is recorded and stored in the global memory of the int type allocated in the second step. Among them, i∈[0,M), j∈[0,N), M and N represent the width and height of the image respectively.

[0072] 3) Sub-pixel edge extraction in connected domains. According to the coordinates of the connected domain circumscribed rectangle obtained in step 6, determine the area where sub-pixel edge extraction needs to be performed in the root matrix root; if the element label value with coordinates (i, j) in this area is non-negative, then the sub-pixel coordinates of the pixel with coordinates (i, j) are calculated using the spatial moment method in the circular non-encoded marker image uploaded to the GPU global memory in step 4, and stored in the global memory allocated in step 2.

[0073] Step 8: GPU-side least squares ellipse fitting: In the processing flow declared in step 2, perform least squares ellipse fitting on the sub-pixel edge extracted in step 7 to obtain the center coordinates of the circular non-coding marker point. The specific implementation steps are as follows:

[0074] 1) Divide the thread grid. Divide the two-dimensional thread grid according to the method described in step 5. The x-dimension and y-dimension of the two-dimensional thread block are set to blockDim.x=256, blockDim.y=1, and the x-dimension and y-dimension of the two-dimensional grid are set to gridDim.x=rectNums, gridDim.y=1, so that one thread is responsible for the ellipse fitting task in a connected domain.

[0075] 2) Calculate the coefficient matrix M and matrix b. The general form of the ellipse equation is:

[0076] x 2 +Axy+By 2 +Cx+Dy+E=0

[0077] Among them, (A, B, C, D, E) are the five parameters of the general equation of the ellipse.

[0078] For a sub-pixel edge point (x k ,y k ), 1≤k≤n, n represents the number of edge pixels in the connected domain, and the objective function F fitted according to the principle of least squares method is:

[0079]

[0080] In order to minimize the objective function F, the partial derivatives of F with respect to each parameter must be 0, that is:

[0081]

[0082] Then we can get the linear equation M×[ABCDE] T =b. Therefore, the coefficient matrix M and matrix b can be derived as:

[0083]

[0084] It can be seen from the above formula that each element of the coefficient matrix M and matrix b is the cumulative sum of the coordinates of n edge points after the product operation. Therefore, the threads in the thread block are first processed by CUDA expansion technology to process the product operation of the edge point coordinates, and then the parallel reduction algorithm is used to solve the cumulative sum of n threads to obtain the coefficient matrix M and matrix b, and store them in the global memory allocated in the second step.

[0085] 3) Solve the parameters of the general equation of the ellipse. Based on the coefficient matrix M and matrix b calculated in step 2), the parameters of the general equation of the ellipse can be calculated:

[0086] [ABCDE] T =(M T M) -1 M T b

[0087] First, the coefficient matrix M is transposed using the column-by-column read and row-by-row storage technique, and then the positive definite symmetric matrix M is transformed using the Gaussian column pivot elimination method. T M is inverted, and finally matrix multiplication is performed to obtain the ellipse parameters, which are stored in the global memory allocated in the second step.

[0088] 4) Solve the center coordinates of the circular non-coded marker point. Calculate the center coordinates (x c ,y c ), the specific formula is as follows:

[0089]

[0090] The calculated center coordinates of the marker points are stored in the global memory allocated in step 2. It should be noted that the event declared in step 2 is inserted into the processing flow to query whether the solution and storage of the center coordinates of the marker points are completed.

[0091] In the ninth step, the detection results are asynchronously transmitted from the GPU to the CPU: the center coordinates of the marker points obtained in the eighth step are asynchronously transmitted from the GPU to the CPU in the transmission stream declared in the second step, completing the real-time detection of a frame of image collected in the third step.

[0092] Query the status of the event inserted in the eighth step processing stream:

[0093] If the event is completed, the center coordinates of the marker point obtained in the eighth step are asynchronously transmitted from the GPU end to the CPU end in the transmission stream declared in the second step, completing the real-time detection of a frame of image collected in the third step.

[0094] Figure 3 The method of the present invention is shown in Figure 2 As a result of the detection, all 65 circular non-encoding markers in the figure are correctly detected (the center of the marker is marked with a cross). The time consumed by the entire detection process (from the first step to the ninth step) is 8ms.

[0095] Table 1 shows the comparison of detection time when the image contains different numbers of landmarks (the image resolution is 2448*2048).

[0096] Table 1 Detection efficiency of images with a resolution of 2448×2048 containing different numbers of landmarks

[0097]

[0098] It can be seen from the table that the CPU-GPU heterogeneous processing method proposed in the present invention has a greatly improved processing speed compared to the CPU side. The more landmarks contained in the image, the greater the improvement, which can be as much as 101 times (when the image contains 96 landmarks).

[0099] The CPU model used in the above processing is Intel 10900K, and the GPU model is RTX3090.

[0100] The tenth step is to determine whether to stop image acquisition. If so, the entire detection process is terminated; otherwise, steps 3 to 9 are repeated to acquire and detect images in real time.

[0101] Since initialization takes a lot of time, when continuous detection is performed, starting from the second image, the detection time is reduced from 8ms to 5ms for an image with a resolution of 2448×2048 (containing 65 circular non-encoded landmarks).

[0102] Therefore, it can be seen that since the present invention combines the advantages of CPU and GPU, some steps suitable for serial processing are placed on the CPU side for processing, and steps suitable for parallel processing are placed on the GPU side for processing, thereby overcoming the defect of the prior art that the detection speed drops sharply with the increase of image resolution or the number of marker points, greatly improving the detection speed of circular non-coding marker points, and realizing real-time online detection of high-resolution images containing a large number of marker points. In addition, during continuous detection, the detection speed can be further improved because a large amount of time occupied by the initialization of the CPU and GPU sides is saved.

Claims

1. A circular non-coding landmark online real-time detection method based on CPU-GPU heterogeneous computing, characterized in that: The steps include: The first step is CPU initialization: allocate the required page-locked memory on the CPU to store the acquired image data and the calculation results on the CPU; The second step is GPU initialization: allocate the required memory area in the GPU global memory to store the image data transmitted from the CPU and the calculation results of the GPU; declare two CUDA streams, one of which is the transmission stream for asynchronous data transmission; the other is the processing stream for the GPU to process the image; declare CUDA events to determine whether the tasks in the stream have been completed; declare the required global variables; The third step is to collect a frame of image in real time: control the camera to collect a frame of circular non-encoded marker image in real time, and store it in the locked page memory allocated on the CPU side in the first step; Step 4: CPU-side integral graph and integral square graph calculation and image data asynchronous transmission: The integral graph and integral square graph calculation are implemented on the CPU side, and the circular non-coding marker point image collected in the third step is asynchronously transmitted from the CPU side to the GPU side in the transmission stream declared in the second step, and the calculated integral graph and integral square graph are asynchronously transmitted from the CPU side to the GPU side; Step 5: GPU-side image binarization: Based on the integral image and integral square image calculated in step 4, the circular non-encoded marker image is binarized in the processing flow declared in step 2. Step 6: GPU-side connected domain marking and analysis: In the processing flow declared in step 2, the connected domain marking and analysis are performed on the binary image obtained in step 5; Step 7: GPU-side edge extraction: extract the whole pixel edges and sub-pixel edges in the connected domain marked in step 6 in the processing flow declared in step 2; Step 8: GPU-side least squares ellipse fitting: In the processing flow declared in step 2, the least squares ellipse fitting is performed on the sub-pixel edge extracted in step 7 to obtain the center coordinates of the circular non-coding marker point; Step 9: asynchronously transmit the detection result from the GPU to the CPU: asynchronously transmit the center coordinates of the marker point obtained in step 8 from the GPU to the CPU in the transmission stream declared in step 2, and complete the real-time detection of a frame of image collected in step 3; Step 10: determine whether to stop image acquisition. If so, end the entire detection process; Otherwise, repeat the third step to the ninth step to collect and detect images in real time.

2. According to claim 1, the online real-time detection method of circular non-coding markers is characterized in that: The method of calculating the integral graph and the integral square graph and transmitting the image data asynchronously on the CPU side in the fourth step is: The integral image is calculated on the CPU side and stored in a locked page memory of size M×N allocated in the second step; at the same time, the circular non-encoded marker image collected in the third step is asynchronously uploaded to the GPU global memory bound to the texture attribute allocated in the second step using the transmission stream declared in the second step; After the integral graph is calculated, the integral square graph is calculated on the CPU side and stored in a locked page memory of size M×N allocated in the second step; at the same time, the calculated integral graph is asynchronously uploaded to a GPU global memory of size M×N×4 allocated in the second step using the transport stream declared in the second step; After the integral square graph is calculated, the integral square graph is asynchronously uploaded to a GPU global memory of size M×N×4 allocated in the second step using the transmission stream declared in the second step; and the event declared in the second step is inserted into the transmission stream to query whether each step of uploading is completed; Where M and N represent the width and height of the image respectively.

3. The online real-time detection method for circular non-coding marker points according to claim 1 is characterized in that: The fifth step of GPU-side image binarization comprises the following steps: 1) Thread meshing The CUDA thread grid consists of thread blocks, which consist of threads. The x dimension of a two-dimensional thread block is set to a multiple of 32, and the y dimension is set to an even number between 1 and 32. Set up a two-dimensional grid. The two dimensions of the two-dimensional grid are calculated as follows: Among them, gridDim.x represents the x dimension of the two-dimensional grid, gridDim.y represents the y dimension of the two-dimensional grid; blockDim.x represents the x dimension of the two-dimensional thread block, blockDim.y represents the y dimension of the two-dimensional thread block; M and N represent the width and height of the image respectively; After the thread grid is divided, each thread in the two-dimensional thread grid is responsible for the binarization processing task of a corresponding pixel in the image. The corresponding relationship between the pixels in the image and the threads in the thread grid is as follows: Wherein, i and j represent the coordinate components of the pixel with coordinates (i, j), i∈[0,M), j∈[0,N); threadIdx represents the thread index in the two-dimensional thread block, threadIdx.x and threadIdx.y are the index components of the x dimension and y dimension respectively, threadIdx.x∈[1,blockDim.x], threadIdx.y∈[1,blockDim.y]; blockIdx represents the thread block index in the two-dimensional thread grid, blockIdx.x and blockIdx.y are the index components of the x dimension and y dimension respectively, blockIdx.x∈[1,gridDim.x], threadIdx.y∈[1,gridDim.y]; 2) Threshold solution One thread in the divided two-dimensional thread grid is responsible for solving the threshold value of a pixel corresponding to it in the image. The solution method is as follows: Where T(i,j) represents the threshold of the pixel with coordinates (i,j); m(i,j) represents the mean grayscale value of the pixel in the square area with the pixel with coordinates (i,j) as the center and a side length of 2r+1, r represents half of the side length of the square area, and its value is r∈[10,20]; s(i,j) represents the standard deviation of the grayscale value of the pixel in the square area; k is a constant, and its value is k∈[0.2,0.5]; The calculation formula of mean m(i,j) is as follows: m(i,j)=[I′(i+r,j+r)+I′(i-r,j-r)-I′(i+r,j-r)-I′(i-r,j+r)] / (2r+1) 2 Where I′(i,j) represents the value of the pixel with coordinates (i,j) in the integral image; The calculation formula of standard deviation s(i,j) is as follows: Where I″(i,j) represents the value of the pixel with coordinates (i,j) in the integral square map; 3) Binarization Based on the calculated threshold, each pixel is binarized: Among them, I b Represents a binary image, I b (i, j) represents the value of the pixel with coordinates (i, j) in the binary image; I(i, j) represents the grayscale value of the pixel with coordinates (i, j) in the circular non-encoded marker image uploaded to the GPU global memory in the fourth step.

4. According to claim 3, the online real-time detection method of circular non-coding markers is characterized in that: The implementation steps of the sixth step of GPU-side connected domain marking and analysis are as follows: 1) Thread meshing Divide the two-dimensional thread grid according to the method described in step 5. After the grid division is completed, each thread in the thread grid is responsible for the connected domain labeling task of a corresponding pixel in the image; 2) Initialize label matrix label and root matrix root If the binary image I obtained in the fifth step b If the value of the pixel with coordinates (i, j) is 255, the element value label(i, j) at the corresponding position of the label matrix label declared in the second step is assigned to j×N+i; otherwise, the element value label(i, j) at the corresponding position is assigned to -1; the values ​​of all elements of the root matrix root declared in the second step are assigned to -1; 3) Update the root-finding matrix For the element with coordinates (i, j) in the label matrix, if its value label(i, j) is not -1, then access the nine elements in the 3×3 square area centered on the element, and assign the non-negative minimum value of the nine element values ​​to the element root(i, j) with coordinates (i, j) in the root matrix root; 4) Update the label matrix Perform root-finding operation on all elements in the root-finding matrix root and update the label matrix. The steps are as follows: for the element with coordinates (i, j) in the root-finding matrix root, if its value root(i, j) is a non-negative number, calculate the new coordinate according to the element value and access the element in the root-finding matrix root corresponding to the new coordinate. The new coordinate corresponding to the non-negative element value root(i, j) is (root(i, j) % N, root(i, j) / N); if the element value at the new coordinate position in the root-finding matrix root is non-negative, repeat the root-finding operation until the calculated new coordinate is the same as the coordinate calculated last time, then end the root-finding operation, assign the element value in the root-finding matrix root corresponding to the new coordinate calculated last time to the element with coordinates (i, j) in the label matrix lable, and assign the Boolean global variable isChanged to true, indicating that there is a label change in the root-finding operation; 5) If the Boolean global variable isChanged in 4) is true, set the Boolean global variable isChanged to false, and repeat 3) and 4). If the Boolean global variable isChanged in 4) is false, the connected domain marking is completed; 6) Connected domain analysis Based on CUDA atomic operations, calculate the minimum value xmin, ymin and the maximum value xmax, ymax of the coordinates of the elements in each marked connected domain in the x dimension and y dimension; if the conditions are met: xmax-xmin>=2×rmin+1, ymax-ymin>=2×rmin+1, xmax-xmin<=2×rmax+1, ymax-ymin<=2×rmax+1, then store these four values ​​(xmin, ymin, xmax, ymax) in the global memory of int4 type allocated in the second step, and add 1 to the global variable rectNums, where xmin, ymin, xmax, ymax are the coordinates of the circumscribed rectangle of the connected domain.

5. The online real-time detection method for circular non-coding marker points according to claim 4 is characterized in that: The implementation steps of the seventh step of GPU-side edge extraction are as follows: 1) Divide the thread grid Divide the two-dimensional thread grid according to the method described in step 5, wherein the x dimension and y dimension of the two-dimensional thread block are set to blockDim.x=32, blockDim.y=32, and the x dimension and y dimension of the two-dimensional grid are set to gridDim.x=rectNums, gridDim.y=1, so that one thread block is responsible for the edge extraction task of one connected domain; 2) Integer pixel edge extraction in connected domain According to the coordinates of the connected domain circumscribed rectangle obtained in the sixth step, determine the area where integer pixel edge extraction needs to be performed in the label matrix label; for an element with coordinates (i, j) and non-negative label value in this area, if the label values ​​of its four neighboring elements are all non-negative, assign the element with coordinates (i, j) in the root matrix root to -1; after all elements in the connected domain are processed, the coordinates corresponding to the elements with non-negative element values ​​in the root matrix root are the coordinates of the integer pixel edge points, and the number of integer pixel edge points in the connected domain is recorded and stored in the global memory of type int allocated in the second step; 3) Sub-pixel edge extraction in connected domains According to the coordinates of the circumscribed rectangle of the connected domain obtained in the sixth step, determine the area where sub-pixel edge extraction needs to be performed in the root-finding matrix root; if the element label value with coordinates (i, j) in this area is non-negative, then the sub-pixel coordinates of the pixel with coordinates (i, j) in the circular non-encoded marker point image uploaded to the GPU global memory in the fourth step are calculated using the spatial moment method and stored in the global memory allocated in the second step.

6. The online real-time detection method for circular non-coding marker points according to claim 5 is characterized in that: The specific implementation steps of the eighth step of least squares ellipse fitting on the GPU are as follows: 1) Divide the thread grid Divide the two-dimensional thread grid according to the method described in step 5, wherein the x dimension and y dimension of the two-dimensional thread block are set to blockDim.x=256, blockDim.y=1, and the x dimension and y dimension of the two-dimensional grid are set to gridDim.x=rectNums, gridDim.y=1, so that one thread is responsible for the ellipse fitting task in one connected domain; 2) Calculate the coefficient matrix M and matrix b For a sub-pixel edge point (x k ,y k ), 1≤k≤n, n represents the number of edge pixels in the connected domain, then the coefficient matrix M and matrix b are calculated as follows: The threads in the thread block are first used to process the product operation of the edge point coordinates using the CUDA expansion technology, and then the parallel reduction algorithm is used to solve the cumulative sum of the n threads to obtain the coefficient matrix M and the matrix b, and store them in the global memory allocated in the second step; each element of the coefficient matrix M and the matrix b is the cumulative sum of the coordinates of the n edge points after the product operation; 3) Solve the parameters of the general equation of the ellipse Based on the calculated coefficient matrix M and matrix b, the parameters of the general equation of the ellipse are calculated: [A B C D E] T =(M T M) -1 M T b Among them, (A, B, C, D, E) are the five parameters of the general equation of the ellipse; When solving, we first use the column-by-column read-by-row storage technique to transpose the coefficient matrix M, and then use the Gaussian column pivot elimination method to transform the positive definite symmetric matrix M. T M is inverted, and finally matrix multiplication is performed to obtain the parameters of the general equation of the ellipse, which are stored in the global memory allocated in the second step; 4) Solve the center coordinates of the circular non-coding mark point The center coordinates (x c ,y c ), the formula is as follows: The calculated center coordinates of the marker points are stored in the global memory allocated in the second step, and the event declared in the second step is inserted into the processing flow to query whether the solution and storage of the center coordinates of the marker points are completed.

Images