[go: up one dir, main page]

CN113052825A - Real-time three-dimensional deformation measurement method based on GPU parallel acceleration - Google Patents

Real-time three-dimensional deformation measurement method based on GPU parallel acceleration Download PDF

Info

Publication number
CN113052825A
CN113052825A CN202110333748.3A CN202110333748A CN113052825A CN 113052825 A CN113052825 A CN 113052825A CN 202110333748 A CN202110333748 A CN 202110333748A CN 113052825 A CN113052825 A CN 113052825A
Authority
CN
China
Prior art keywords
interest
point
deformation
dimensional
gpu
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Granted
Application number
CN202110333748.3A
Other languages
Chinese (zh)
Other versions
CN113052825B (en
Inventor
董守斌
林傲宇
蒋震宇
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
South China University of Technology SCUT
Original Assignee
South China University of Technology SCUT
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by South China University of Technology SCUT filed Critical South China University of Technology SCUT
Priority to CN202110333748.3A priority Critical patent/CN113052825B/en
Publication of CN113052825A publication Critical patent/CN113052825A/en
Application granted granted Critical
Publication of CN113052825B publication Critical patent/CN113052825B/en
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Classifications

    • G—PHYSICS
    • G06—COMPUTING OR CALCULATING; COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00—Image analysis
    • G06T7/0002—Inspection of images, e.g. flaw detection
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01B—MEASURING LENGTH, THICKNESS OR SIMILAR LINEAR DIMENSIONS; MEASURING ANGLES; MEASURING AREAS; MEASURING IRREGULARITIES OF SURFACES OR CONTOURS
    • G01B11/00—Measuring arrangements characterised by the use of optical techniques
    • G01B11/16—Measuring arrangements characterised by the use of optical techniques for measuring the deformation in a solid, e.g. optical strain gauge
    • G—PHYSICS
    • G06—COMPUTING OR CALCULATING; COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T1/00—General purpose image data processing
    • G06T1/20—Processor architectures; Processor configuration, e.g. pipelining
    • G—PHYSICS
    • G06—COMPUTING OR CALCULATING; COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T1/00—General purpose image data processing
    • G06T1/60—Memory management
    • G—PHYSICS
    • G06—COMPUTING OR CALCULATING; COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00—Image analysis
    • G06T7/80—Analysis of captured images to determine intrinsic or extrinsic camera parameters, i.e. camera calibration
    • G06T7/85—Stereo camera calibration
    • G—PHYSICS
    • G06—COMPUTING OR CALCULATING; COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00—Indexing scheme for image analysis or image enhancement
    • G06T2207/10—Image acquisition modality
    • G06T2207/10004—Still image; Photographic image
    • G06T2207/10012—Stereo images
    • Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02T—CLIMATE CHANGE MITIGATION TECHNOLOGIES RELATED TO TRANSPORTATION
    • Y02T90/00—Enabling technologies or technologies with a potential or indirect contribution to GHG emissions mitigation

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Quality & Reliability (AREA)
  • Image Processing (AREA)
  • Image Analysis (AREA)

Abstract

本发明公开了一种基于GPU并行加速的实时三维变形测量方法，包括以下步骤：1)通过立体标定获得投影矩阵；2)拍摄物体变形前后的图像；3)选择兴趣区域和兴趣点；4)传输投影矩阵、图像和兴趣点到GPU；5)计算兴趣点变形前的三维坐标；6)计算兴趣点在变形中各个时刻的三维变形；7)将三维变形数据传输回CPU。该方法通过将变形前左相机图像作为所有匹配中的参考图像，从而兴趣点对应的Hessian矩阵等IC‑GN预计算的数据得以复用；基于CUDA异构计算平台开发的GPU加速变形测量程序可以发挥GPU硬件设备的计算性能，针对GPU程序的访存等优化技术使得三维变形测量的计算速度大大提高，满足了实时三维变形测量的需求。The invention discloses a real-time three-dimensional deformation measurement method based on GPU parallel acceleration. Transfer the projection matrix, the image and the interest point to the GPU; 5) Calculate the three-dimensional coordinates of the interest point before deformation; 6) Calculate the three-dimensional deformation of the interest point at each moment in the deformation; 7) Transfer the three-dimensional deformation data back to the CPU. This method uses the left camera image before deformation as the reference image in all matching, so that the IC-GN pre-computed data such as the Hessian matrix corresponding to the interest point can be reused; the GPU-accelerated deformation measurement program developed based on the CUDA heterogeneous computing platform can Taking advantage of the computing performance of GPU hardware devices and optimizing the memory access of GPU programs, the calculation speed of 3D deformation measurement is greatly improved, which meets the needs of real-time 3D deformation measurement.

Description

Real-time three-dimensional deformation measurement method based on GPU parallel acceleration
Technical Field
The invention relates to the technical field of optical measurement, in particular to a real-time three-dimensional deformation measurement method based on GPU parallel acceleration.
Background
In the fields of science and engineering, the three-dimensional digital image correlation method is widely applied to three-dimensional deformation measurement due to the advantages of simple device, non-contact type and the like. The three-dimensional digital image correlation method can measure the three-dimensional appearance and the full-field three-dimensional deformation of a curved surface object, and has very rich application scenes. However, since it directly processes a high-resolution digital image, the amount of calculation thereof is also considerable, with a problem that the calculation takes a long time. With the development of digital image acquisition technology, the image resolution and sampling rate are improved, and the problem is more prominent, so that the application of a three-dimensional digital image correlation method in some real-time monitoring scenes and the like is limited. Secondly, as a measuring method, it is also important to maintain its high accuracy.
In recent years, researchers have made great efforts to improve the computational efficiency of three-dimensional digital image correlation methods. However, the results of these studies are not satisfactory and always give a trade-off between accuracy and efficiency. There are a number of existing schemes that improve computational efficiency by optimizing correlation algorithms, reducing redundant computations, and using multi-threading techniques to speed up programs. However, because the computing capability of the multi-core processor is limited and the adopted computing strategy is simple, the highest computing speed reported at present is only 50000POI/s, and the requirement of real-time processing is far from being met; the speed is achieved on the premise of only using a simple first-order shape function, and the precision is not as good as that of a scheme using a second-order shape function.
Disclosure of Invention
The invention aims to overcome the defects of the prior art and provides a real-time three-dimensional deformation measuring method based on GPU parallel acceleration. In all matching, the left camera image L before deformation is used0As a reference image, the pre-calculation data of the IC-GN algorithm can be repeatedly used, so that a large amount of calculation is saved. The method has the advantages that the algorithm is accelerated by developing a GPU program based on CUDA, the computing capability of hardware is fully exerted, real-time three-dimensional deformation measurement is achieved, the computing speed can exceed 40 frames per second when the number of interest points is about 10000, and the problems that the computing speed of an existing three-dimensional deformation measurement method is low and real-time measurement requirements cannot be met are solved.
In order to achieve the purpose, the technical scheme provided by the invention is as follows: a real-time three-dimensional deformation measurement method based on GPU parallel acceleration comprises the following steps:
1) using Zhangyingyou scaling method to make stereo scaling on left and right fixed cameras to obtain projection matrix M of left and right camerasL、MR；
2) Synchronously shooting left and right images L of the surface of the target object before deformation by using the two cameras0、 R0And the left and right images L at the ith moment in the deformation processi、RiWherein i ═ 1,2,3, …, n;
3) in the image L0Selecting an interest area, and taking a batch of interest points P at equal intervals in the interest areaL0；
4) Copying a projection matrix of the camera, all shot images and interest points to a GPU;
5) on the GPU, for each interest point PL0Firstly, an IC-GN algorithm is used for searching the image R0Corresponding point P onR0Then using PL0、PR0The coordinates of the two points are calculated by a triangulation method to obtain the three-dimensional coordinates P of the interest point before deformationW0；
6) On the GPU, the following operations are carried out on each time i in the deformation process: first, for each point of interest PL0Use the IC-GN algorithm to find it in the image Li、RiCorresponding point P onLiAnd PRiThen using PLi、PRiThe coordinates of the two points are calculated by triangulation to obtain the three-dimensional coordinates P of the interest point at the moment iwiFinally with PWiMinus PW0Obtaining three-dimensional deformation data D of the interest point at the moment iWi；
7) And copying the three-dimensional deformation data of all the interest points at each moment back to the CPU, so as to obtain the three-dimensional deformation of the surface of the object at each moment.
In step 1), the projection matrix has a size of 3 × 4, which represents the relationship between the three-dimensional space coordinates and the two-dimensional coordinates on the camera image.
In step 3), the interest area refers to an area which needs to be measured and is designated by a user, and the interest point interval in the area is also selected according to the measurement requirement.
In step 4), the shot image is firstly stored in a page-locked memory allocated by using a runtime function cudamallocost in the CUDA, and then the data is copied to the GPU by using a runtime function cudaMemcpy in the CUDA; the camera parameters are copied to the constant memory of the GPU, enabling it to be accessed at high speed.
In step 5) and step 6), the IC-GN algorithm is used, wherein the input is a batch of interest points, reference images and target images, and a CUDA thread block is used to process a computation task corresponding to an interest point, including the following:
a. carrying out precalculation: for each interest point, calculating a corresponding Hessian matrix and storing data; using a CUDA thread block to complete the calculation task of each interest point, wherein the Hessian matrix is
Figure RE-GDA0003079967210000031
Figure RE-GDA0003079967210000032
Representing the accumulation of all pixel positions within a reference sub-area, which is a 33 x 33 sub-image centered at the point of interest; ψ is the coordinate of the point of interest, # is the local coordinate in the reference sub-area, # R (ψ + ζ) is the gradient of the reference image,
Figure RE-GDA0003079967210000033
is a Jacobian matrix, and T represents the transposition of the matrix;
b. and (c) for each interest point, setting the coordinate of the interest point as (x, y), and estimating the initial value p of the deformation vector of the interest point by using an image feature assisted method (u, u)x,uy,uxx,uxy,uyy,v,vx,vy,vxx,vxy,vyy) (ii) a Wherein u, v are translation amounts; u. ofx,uy,vx,vyIs the first order gradient component; u. ofxx,uxy,uyy,vxx,vxy,vyyIs the second order gradient component; here, a CUDA thread is used to accomplish a programCalculating the interest points;
c. for each interest point, iteratively updating the corresponding deformation vector p according to the following steps:
c1, calculating the deformation vector increment delta p:
Figure RE-GDA0003079967210000034
wherein H-1The inverse of the Hessian matrix is represented,
Figure RE-GDA0003079967210000035
and
Figure RE-GDA0003079967210000036
normalized coefficients for the reference and target sub-regions respectively,
Figure RE-GDA0003079967210000037
representing the target sub-area after subtraction of the mean value of the gray levels,
Figure RE-GDA0003079967210000038
representing the reference subarea after subtracting the gray mean value, and W (zeta; p) represents a transformation function from the local coordinate of the reference subarea to the local coordinate of the target subarea;
c2, calculating new transformation function W (ζ; p') as W (W) by using Δ p-1(ζ; Δ p); p) wherein W-1(ζ; Δ p) is the inverse of W (ζ; Δ p); then extracting a new deformation vector p 'from the new transformation function W (zeta; p');
c3, updating p, namely making p equal to p';
c4, repeating c1 to c3 until | | | Δ p | | < 0.001, and | | · | | | represents the modular length of the vector;
d. and (3) calculating the coordinates of the interest point at the corresponding point of the target map as (x ', y'), wherein x 'is x + u, and y' is y + v.
In step 6), the IC-GN algorithm used does not need to be pre-calculated, and the data stored after the IC-GN algorithm is pre-calculated in step 5) is directly used.
In step 5) and step 6), triangulation methods are used, including the following:
a. let the coordinate of the interest point on the left view be (x)L,yL) The point coordinate on the right view is (x)R,yR) The corresponding three-dimensional coordinate is (X)W,YW,ZW)；
b. (X) is calculated in the following mannerW,YW,ZW) Each CUDA thread calculates a three-dimensional coordinate;
Figure RE-GDA0003079967210000041
wherein A is a coefficient matrix, and
Figure RE-GDA0003079967210000042
b is a matrix of constant terms, an
Figure RE-GDA0003079967210000043
mL rcRepresenting a projection matrix MLM in the r-th row and c-th columnR rcRepresenting a projection matrix MRWhere r is 1,2,3, and c is 1,2,3, 4.
Compared with the prior art, the invention has the following advantages and beneficial effects:
1. in all matching of three-dimensional deformation measurement, the invention leads the left camera image L before deformation to be0As a reference image, data obtained by IC-GN pre-calculation such as Hessian matrix and the like corresponding to the interest points can be repeatedly utilized in the processing of images at all times in the three-dimensional deformation process, so that a large amount of calculation is saved, the processing speed is accelerated, and the calculation time is shortened.
2. The GPU program is developed based on the CUDA platform, so that development cost can be saved, development difficulty is reduced, meanwhile, the special CUDA program and corresponding program optimization can furthest exert the calculation performance of GPU hardware, and the requirement of measuring three-dimensional deformation in real time is met.
Detailed Description
The present invention will be described in further detail with reference to examples, but the embodiments of the present invention are not limited thereto.
The real-time three-dimensional deformation measurement method based on GPU parallel acceleration provided by the embodiment comprises the following steps:
1) using a Zhangyingyou calibration method to carry out three-dimensional calibration on the left and right fixed cameras to obtain the projection matrixes M of the left and right camerasL、MR。
2) Synchronously shooting left and right images L of the surface of the target object before deformation by using the two cameras0、 R0And the left and right images L at the ith moment in the deformation processi、RiWhere i is 1,2,3, …, n.
3) In the image L0Selecting an interest area, and taking a batch of interest points P at equal intervals in the interest areaL0(ii) a The interest area refers to an area which is specified by a user and needs to be measured, and the interest point interval in the area is also selected according to the measurement requirement.
4) Copying a projection matrix of the camera, all shot images and interest points to a GPU; the method comprises the steps that a shot image is stored in a page locking memory distributed by using a runtime function cudaMallocost in a CUDA (compute unified device architecture), and then data are copied to a GPU (graphics processing unit) by using a runtime function cudaMemcpy function in the CUDA; the projection matrix of the camera is copied to the constant memory of the GPU, allowing it to be accessed at high speed.
5) On the GPU, for each interest point PL0Firstly, an IC-GN algorithm is used for searching the image R0Corresponding point P onR0Then using PL0、PR0The coordinates of the two points are calculated by a triangulation method to obtain the three-dimensional coordinates P of the interest point before deformationW0。
6) On the GPU, the following operations are carried out on each time i in the deformation process: first, for each point of interest PL0Use the IC-GN algorithm to find it in the image Li、RiCorresponding point P onLiAnd PRiThen using PLi、PRiThe coordinates of the two points are calculated by triangulation to obtain the three-dimensional coordinate of the interest point at the time iMark PwiFinally with PWiMinus PW0Obtaining three-dimensional deformation data D of the interest point at the moment iWi(ii) a The IC-GN algorithm used does not need to be pre-calculated, and the data stored after the IC-GN algorithm is pre-calculated in the step 5) is directly utilized.
7) And copying the three-dimensional deformation data of all the interest points at each moment back to the CPU, so as to obtain the three-dimensional deformation of the surface of the object at each moment.
In step 5) and step 6), the IC-GN algorithm is used, wherein the input is a batch of interest points, reference images and target images, and a CUDA thread block is used to process a computation task corresponding to an interest point, including the following:
a. carrying out precalculation: for each interest point, calculating a corresponding Hessian matrix and storing data; using a CUDA thread block to complete the calculation task of each interest point, wherein the Hessian matrix is
Figure RE-GDA0003079967210000061
Figure RE-GDA0003079967210000062
Representing the accumulation of all pixel positions within a reference sub-area, which is a 33 x 33 sub-image centered at the point of interest; ψ is the coordinate of the point of interest, # is the local coordinate in the reference sub-area, # R (ψ + ζ) is the gradient of the reference image,
Figure RE-GDA0003079967210000063
is a Jacobian matrix, and T represents the transposition of the matrix;
b. and (c) for each interest point, setting the coordinate of the interest point as (x, y), and estimating the initial value p of the deformation vector of the interest point by using an image feature assisted method (u, u)x,uy,uxx,uxy,uyy,v,vx,vy,vxx,vxy,vyy) (ii) a Wherein u, v are translation amounts; u. ofx,uy,vx,vyIs the first order gradient component; u. ofxx,uxy,uyy,vxx,vxy,vyyIs twoAn order gradient component; here, a CUDA thread is used to perform a point of interest computing task;
c. for each interest point, iteratively updating the corresponding deformation vector p according to the following steps:
c1, calculating the deformation vector increment delta p:
Figure RE-GDA0003079967210000071
wherein H-1The inverse of the Hessian matrix is represented,
Figure RE-GDA0003079967210000072
and
Figure RE-GDA0003079967210000073
normalized coefficients for the reference and target sub-regions respectively,
Figure RE-GDA0003079967210000074
representing the target sub-area after subtraction of the mean value of the gray levels,
Figure RE-GDA0003079967210000075
representing the reference subarea after subtracting the gray mean value, and W (zeta; p) represents a transformation function from the local coordinate of the reference subarea to the local coordinate of the target subarea;
c2, calculating new transformation function W (ζ; p') as W (W) by using Δ p-1(ζ; Δ p); p) wherein W-1(ζ; Δ p) is the inverse of W (ζ; Δ p); then extracting a new deformation vector p 'from the new transformation function W (zeta; p');
c3, updating p, namely making p equal to p';
c4, repeating c1 to c3 until | | | Δ p | | < 0.001, and | | · | | | represents the modular length of the vector;
d. and (3) calculating the coordinates of the interest point at the corresponding point of the target map as (x ', y'), wherein x 'is x + u, and y' is y + v.
In step 5) and step 6), triangulation methods are used, including the following:
a. setting the point of interest atThe coordinates on the left view are (x)L,yL) The point coordinate on the right view is (x)R,yR) The corresponding three-dimensional coordinate is (X)W,YW,ZW)；
b. (X) is calculated in the following mannerW,YW,ZW) Each CUDA thread calculates a three-dimensional coordinate;
Figure RE-GDA0003079967210000076
wherein A is a coefficient matrix, and
Figure RE-GDA0003079967210000077
b is a matrix of constant terms, an
Figure RE-GDA0003079967210000081
mL rcRepresenting a projection matrix MLM in the r-th row and c-th columnR rcRepresenting a projection matrix MRWhere r is 1,2,3, and c is 1,2,3, 4.
The above embodiments are preferred embodiments of the present invention, but the present invention is not limited to the above embodiments, and any other changes, modifications, substitutions, combinations, and simplifications which do not depart from the spirit and principle of the present invention should be construed as equivalents thereof, and all such changes, modifications, substitutions, combinations, and simplifications are intended to be included in the scope of the present invention.

Claims (7)

1.一种基于GPU并行加速的实时三维变形测量方法，其特征在于，包括以下步骤：1. a real-time three-dimensional deformation measurement method based on GPU parallel acceleration, is characterized in that, comprises the following steps: 1)使用张正友标定法对左右两个固定的相机进行立体标定，获得左、右相机的投影矩阵ML、MR；1) Use Zhang Zhengyou's calibration method to perform stereo calibration on the left and right fixed cameras to obtain the projection matrices ML and MR of the left and right cameras; 2)使用上述两个相机同步拍摄目标物体表面在变形前的左、右图像L0、R0和变形过程中第i个时刻的左、右图像Li、Ri，其中i＝1,2,3,…,n；2) Use the above two cameras to simultaneously capture the left and right images L 0 , R 0 of the surface of the target object before deformation and the left and right images Li and R i at the ith moment during the deformation process, where i=1,2 ,3,…,n; 3)在图像L0上选择兴趣区域，在该区域中等间隔取一批兴趣点PL0；3) Select a region of interest on the image L 0 , and take a batch of interest points P L0 at medium intervals in the region; 4)将相机的投影矩阵、所有拍摄的图像和兴趣点拷贝到GPU上；4) Copy the camera's projection matrix, all captured images and points of interest to the GPU; 5)在GPU上，对每个兴趣点PL0，先使用IC-GN算法寻找其在图像R0上的对应点PR0，然后使用PL0、PR0这两个点的坐标，通过三角测量法计算得到兴趣点在变形前的三维坐标PW0；5) On the GPU, for each point of interest P L0 , first use the IC-GN algorithm to find its corresponding point P R0 on the image R 0 , and then use the coordinates of the two points P L0 and P R0 to measure by triangulation The three-dimensional coordinate P W0 of the interest point before the deformation is obtained by calculating the method; 6)在GPU上，对变形过程中的每个时刻i进行如下操作：首先，对每个兴趣点PL0使用IC-GN算法寻找其在图像Li、Ri上的对应点PLi和PRi，然后使用PLi、PRi这两个点的坐标通过三角测量法计算出兴趣点在时刻i的三维坐标Pwi，最后用PWi减去PW0获得兴趣点在时刻i的三维变形数据DWi；6) On the GPU, perform the following operations on each moment i in the deformation process: First, use the IC-GN algorithm for each interest point P L0 to find its corresponding points P Li and P on the images Li and Ri . Ri , then use the coordinates of the two points P Li and P Ri to calculate the three-dimensional coordinate P wi of the point of interest at time i by triangulation, and finally subtract P W0 from P Wi to obtain the three-dimensional deformation data of the point of interest at time i D Wi ; 7)将所有兴趣点在各个时刻的三维变形数据拷贝回CPU，即可获得物体表面在各个时刻的三维变形。7) Copy the three-dimensional deformation data of all interest points at each moment back to the CPU, and then the three-dimensional deformation of the object surface at each moment can be obtained. 2.根据权利要求1所述的一种基于GPU并行加速的实时三维变形测量方法，其特征在于：在步骤1)中，投影矩阵大小为3×4，表示三维空间坐标和相机图像上二维坐标之间的关系。2. a kind of real-time three-dimensional deformation measurement method based on GPU parallel acceleration according to claim 1, is characterized in that: in step 1), the projection matrix size is 3 × 4, represents two-dimensional coordinates on three-dimensional space coordinates and camera image relationship between coordinates. 3.根据权利要求1所述的一种基于GPU并行加速的实时三维变形测量方法，其特征在于：在步骤3)中，兴趣区域是指用户指定的需要测量的区域，区域中的兴趣点间隔也是根据测量需求选择的。3. a kind of real-time three-dimensional deformation measurement method based on GPU parallel acceleration according to claim 1, is characterized in that: in step 3) in, the area of interest refers to the area that needs to be measured that the user specifies, the interest point interval in the area It is also selected according to the measurement needs. 4.根据权利要求1所述的一种基于GPU并行加速的实时三维变形测量方法，其特征在于：在步骤4)中，拍摄的图像首先被存储在使用CUDA中的运行时函数cudaMallocHost分配的锁页内存中，然后使用CUDA中的运行时函数cudaMemcpy函数，将数据拷贝到GPU；相机的投影矩阵被拷贝到GPU的常量内存，让其能够被高速访问。4. a kind of real-time three-dimensional deformation measurement method based on GPU parallel acceleration according to claim 1, is characterized in that: in step 4), the image of taking is first stored in the lock that the runtime function cudaMallocHost in use CUDA distributes In the page memory, and then use the runtime function cudaMemcpy function in CUDA to copy the data to the GPU; the camera's projection matrix is copied to the constant memory of the GPU, so that it can be accessed at high speed. 5.根据权利要求1所述的一种基于GPU并行加速的实时三维变形测量方法，其特征在于：在步骤5)和步骤6)中，使用的IC-GN算法中，输入是一批兴趣点、参考图像和目标图像，使用一个CUDA线程块来处理一个兴趣点相应的计算任务，包括以下内容：5. a kind of real-time three-dimensional deformation measurement method based on GPU parallel acceleration according to claim 1, is characterized in that: in step 5) and step 6), in the IC-GN algorithm of use, input is a batch of interest points , the reference image and the target image, using a CUDA thread block to process the computation tasks corresponding to a point of interest, including the following: a、进行预计算：对每个兴趣点，计算相应的Hessian矩阵，并存储数据；使用一个CUDA线程块来完成每个兴趣点的计算任务，其中Hessian矩阵为
Figure FDA0002997371300000021
Figure FDA0002997371300000022
表示对参考子区内的所有像素位置的累加，参考子区是以兴趣点为中心的大小为33×33的子图像；ψ是兴趣点的坐标，ζ是参考子区中的局部坐标，
Figure FDA0002997371300000023
是参考图像的梯度，
Figure FDA0002997371300000024
是雅克比矩阵，T表示矩阵的转置；a. Perform pre-calculation: For each interest point, calculate the corresponding Hessian matrix and store the data; use a CUDA thread block to complete the calculation task of each interest point, where the Hessian matrix is
Figure FDA0002997371300000021
Figure FDA0002997371300000022
Represents the accumulation of all pixel positions in the reference sub-region. The reference sub-region is a sub-image with a size of 33 × 33 centered on the point of interest; ψ is the coordinate of the point of interest, ζ is the local coordinate in the reference sub-region,
Figure FDA0002997371300000023
is the gradient of the reference image,
Figure FDA0002997371300000024
is the Jacobian matrix, and T represents the transpose of the matrix; b、对每个兴趣点，设其坐标为(x,y)，使用图像特征辅助的方法估计它的变形矢量初值p＝(u,ux,uy,uxx,uxy,uyy,v,vx,vy,vxx,vxy,vyy)；其中u,v是平移量；ux,uy,vx,vy是一阶梯度分量；uxx,uxy,uyy,vxx,vxy,vyy是二阶梯度分量；这里使用一个CUDA线程来完成一个兴趣点的计算任务；b. For each interest point, set its coordinates as (x, y), and use the image feature-assisted method to estimate the initial value of its deformation vector p=(u, u x , u y , u xx , u xy , u yy ,v,v x ,v y ,v xx ,v xy ,v yy ); where u,v are translations; u x ,u y ,v x ,v y are first-order gradient components; u xx ,u xy , u yy , v xx , v xy , v yy are the second-order gradient components; here a CUDA thread is used to complete the calculation task of a point of interest; c、对每个兴趣点，按以下步骤迭代更新其对应的变形矢量p：c. For each interest point, iteratively update its corresponding deformation vector p according to the following steps: c1、计算变形矢量增量Δp：c1. Calculate the deformation vector increment Δp:
Figure FDA0002997371300000025
Figure FDA0002997371300000025
 其中，H-1表示Hessian矩阵的逆，
Figure FDA0002997371300000031
和
Figure FDA0002997371300000032
分别是参考子区和目标子区的归一化系数，
Figure FDA0002997371300000033
表示减去灰度均值后的目标子区，
Figure FDA0002997371300000034
表示减去灰度均值后的参考子区，W(ζ；p)表示参考子区局部坐标到目标子区局部坐标的变换函数；where H -1 represents the inverse of the Hessian matrix,
Figure FDA0002997371300000031
and
Figure FDA0002997371300000032
are the normalization coefficients of the reference sub-region and the target sub-region, respectively,
Figure FDA0002997371300000033
represents the target sub-area after subtracting the gray mean value,
Figure FDA0002997371300000034
represents the reference sub-region after subtracting the gray mean value, W(ζ; p) represents the transformation function from the local coordinates of the reference sub-region to the local coordinates of the target sub-region; c2、使用Δp计算新的变换函数W(ζ；p')＝W(W-1(ζ；Δp)；p)，其中W-1(ζ；Δp)是W(ζ；Δp)的逆变换；然后从新的变换函数W(ζ；p')中提取新的变形矢量p'；c2. Use Δp to calculate a new transformation function W(ζ; p')=W(W -1 (ζ; Δp); p), where W -1 (ζ; Δp) is the inverse transformation of W(ζ; Δp) ; then extract a new deformation vector p' from the new transformation function W(ζ; p'); c3、更新p，即令p＝p'；c3, update p, that is, let p=p'; c4、重复c1至c3，直到||Δp||＜0.001，||·||表示向量的模长；c4. Repeat c1 to c3 until ||Δp||<0.001, and ||·|| represents the modulo length of the vector; d、计算兴趣点在目标图的对应点坐标为(x',y')，其中x'＝x+u，y'＝y+v。d. Calculate the coordinates of the corresponding point of the interest point in the target image as (x', y'), where x'=x+u, y'=y+v. 6.根据权利要求1所述的一种基于GPU并行加速的实时三维变形测量方法，其特征在于：在步骤6)中，使用的IC-GN算法无须进行预计算，直接利用步骤5)IC-GN算法预计算后所存储的数据。6. a kind of real-time three-dimensional deformation measurement method based on GPU parallel acceleration according to claim 1, is characterized in that: in step 6), the IC-GN algorithm of use need not carry out pre-calculation, directly utilizes step 5) IC-GN The data stored after precomputing the GN algorithm. 7.根据权利要求1所述的一种基于GPU并行加速的实时三维变形测量方法，其特征在于：在步骤5)和步骤6)中，使用的三角测量法，包括以下内容：7. a kind of real-time three-dimensional deformation measurement method based on GPU parallel acceleration according to claim 1, is characterized in that: in step 5) and step 6), the triangulation method of using, comprises the following content: a、设兴趣点在左视图上的坐标为(xL,yL)，在右视图上的点坐标为(xR,yR)，相应的三维坐标为(XW,YW,ZW)；a. Let the coordinates of the point of interest on the left view be (x L , y L ), the coordinates of the point on the right view be (x R , y R ), and the corresponding three-dimensional coordinates are (X W , Y W , Z W ); b、按以下方式计算(XW,YW,ZW)，每个CUDA线程计算一个三维坐标；b. Calculate (X W , Y W , Z W ) as follows, and each CUDA thread calculates a three-dimensional coordinate;
Figure FDA0002997371300000035
Figure FDA0002997371300000035
 其中，A是系数矩阵，且
Figure FDA0002997371300000036
b是常数项矩阵，且
Figure FDA0002997371300000041
mL rc表示投影矩阵ML的在第r行第c列的元素，mR rc表示投影矩阵MR的在第r行第c列的元素，其中r＝1,2,3，c＝1,2,3,4。where A is the coefficient matrix, and
Figure FDA0002997371300000036
b is a constant term matrix, and
Figure FDA0002997371300000041
m L rc represents the element in the rth row and cth column of the projection matrix ML , m R rc represents the element in the rth row and the cth column of the projection matrix MR , where r=1, 2, 3, c=1 ,2,3,4.
CN202110333748.3A 2021-03-29 2021-03-29 A Real-time 3D Deformation Measurement Method Based on GPU Parallel Acceleration Active CN113052825B (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202110333748.3A CN113052825B (en) 2021-03-29 2021-03-29 A Real-time 3D Deformation Measurement Method Based on GPU Parallel Acceleration

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202110333748.3A CN113052825B (en) 2021-03-29 2021-03-29 A Real-time 3D Deformation Measurement Method Based on GPU Parallel Acceleration

Publications (2)

Publication Number Publication Date
CN113052825A true CN113052825A (en) 2021-06-29
CN113052825B CN113052825B (en) 2023-04-21

Family

ID=76516013

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202110333748.3A Active CN113052825B (en) 2021-03-29 2021-03-29 A Real-time 3D Deformation Measurement Method Based on GPU Parallel Acceleration

Country Status (1)

Country Link
CN (1) CN113052825B (en)

Cited By (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN114526682A (en) * 2022-01-13 2022-05-24 华南理工大学 Deformation measurement method based on image feature enhanced digital volume image correlation method
CN115861024A (en) * 2022-11-15 2023-03-28 西安电子科技大学 A method for online real-time detection of circular non-coded markers

Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2000007501A1 (en) * 1998-08-03 2000-02-17 Cardiac Pathways Corporation A dynamically alterable three-dimensional graphical model of a body region
WO2010005973A2 (en) * 2008-07-07 2010-01-14 The Johns Hopkins University Automated surface-based anatomical analysis based on atlas-based segmentation of medical imaging
CN102003946A (en) * 2010-09-02 2011-04-06 北京航空航天大学 High-temperature three-dimensional digital image related measurement system and measurement method
CN104634266A (en) * 2013-11-08 2015-05-20 南京化工职业技术学院 Mechanical sealing end surface deformation measurement system based on binocular vision DIC and measurement method thereof
CN110414573A (en) * 2019-07-09 2019-11-05 北京卫星制造厂有限公司 A kind of Route Dependence digital picture Image Matching accelerated based on GPU

Patent Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2000007501A1 (en) * 1998-08-03 2000-02-17 Cardiac Pathways Corporation A dynamically alterable three-dimensional graphical model of a body region
WO2010005973A2 (en) * 2008-07-07 2010-01-14 The Johns Hopkins University Automated surface-based anatomical analysis based on atlas-based segmentation of medical imaging
CN102003946A (en) * 2010-09-02 2011-04-06 北京航空航天大学 High-temperature three-dimensional digital image related measurement system and measurement method
CN104634266A (en) * 2013-11-08 2015-05-20 南京化工职业技术学院 Mechanical sealing end surface deformation measurement system based on binocular vision DIC and measurement method thereof
CN110414573A (en) * 2019-07-09 2019-11-05 北京卫星制造厂有限公司 A kind of Route Dependence digital picture Image Matching accelerated based on GPU

Non-Patent Citations (3)

* Cited by examiner, † Cited by third party
Title
JIANG ZHENYU 等: "Path-independent digital image correlation with high accuracy，speed and robustness", <<网页>> *
LINGQI ZHANG等: "High accuracy digital image correlation powered by GPU-based parallel computing", <<网页>> *
黄健文: "高精度数字图像相关法的并行算法设计及异构计算实现" *

Cited By (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN114526682A (en) * 2022-01-13 2022-05-24 华南理工大学 Deformation measurement method based on image feature enhanced digital volume image correlation method
CN115861024A (en) * 2022-11-15 2023-03-28 西安电子科技大学 A method for online real-time detection of circular non-coded markers

Also Published As

Publication number Publication date
CN113052825B (en) 2023-04-21

Similar Documents

Publication Publication Date Title
CN115690324A (en) Neural radiation field reconstruction optimization method and device based on point cloud
CN113160294B (en) Image scene depth estimation method, device, terminal device and storage medium
CN106981083B (en) The substep scaling method of Binocular Stereo Vision System camera parameters
JP3651590B2 (en) Method for restoring 3D scene structure and camera motion directly from points, lines and / or from image intensity
CN110070598B (en) Mobile terminal for 3D scanning reconstruction and 3D scanning reconstruction method thereof
Coorg et al. Spherical mosaics with quaternions and dense correlation
WO2022198901A1 (en) Digital speckle correlation rapid implementation method for extracting seed points on basis of grid
CN107230225A (en) The method and apparatus of three-dimensional reconstruction
CN105678757B (en) A kind of ohject displacement measuring method
CN111862299A (en) Human body three-dimensional model construction method, device, robot and storage medium
CN118887346B (en) Gaussian radiation field three-dimensional reconstruction method based on geometric distillation of pre-training optical flow model
CN108564620B (en) A Scene Depth Estimation Method for Light Field Array Cameras
CN115409949B (en) Model training method, perspective image generation method, device, equipment and medium
CN113052825B (en) A Real-time 3D Deformation Measurement Method Based on GPU Parallel Acceleration
CN114509018A (en) Full-field real-time bridge deflection measurement method
CN110375732A (en) Monocular camera pose measurement method based on inertial measurement unit and point line characteristics
WO2018107427A1 (en) Rapid corresponding point matching method and device for phase-mapping assisted three-dimensional imaging system
WO2018133027A1 (en) Grayscale constraint-based method and apparatus for integer-pixel search for three-dimensional digital speckle pattern
CN115601438A (en) External parameter calibration method, device and autonomous mobile equipment
CN118485779A (en) Sparse view angle indoor scene reconstruction method based on mutual information optimization neural radiation field
CN109741389B (en) A Local Stereo Matching Method Based on Region-Based Matching
CN109978957B (en) Binocular System Calibration Method Based on Quantum Behavior Particle Swarm
CN114820810A (en) Analysis method based on Tsai camera plane calibration algorithm
CN110148168A (en) A trinocular camera depth image processing method based on large and small double baselines
CN109300148A (en) Multi-source image registration method based on method collaboration

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant