CUDA architecture parallel optimization three-dimensional deformation measurement method based on novel correlation function constraint
1. A CUDA architecture parallel optimization stereo deformation measurement method based on novel correlation function constraint is characterized by comprising the following steps:
s01: obtaining a relative pose parameter relation among multiple cameras;
s02: establishing a novel correlation function of speckle three-dimensional registration based on a joint constraint relation between cameras, limiting the search of three-dimensional registration points between image pairs in an area near an epipolar line by the novel correlation function, and performing time sequence registration and three-dimensional registration on sampling scattered spots in images of different sequences;
s03: and carrying out three-dimensional reconstruction on the space coordinates of the scattered spots.
2. The method for measuring the stereo deformation based on the parallel optimization of the CUDA architecture under the constraint of the novel correlation function as claimed in claim 1, wherein the method for obtaining the relative pose parameter relationship among the multiple cameras in the step S01 comprises:
s11: setting a measurement system to be composed of eta +1 cameras, wherein the main camera is xi, and the rest eta cameras are all slave cameras, and establishing a main camera collinearity equation based on a space point P according to a perspective projection principle:
λξ[uξ vξ 1]＝Mξ[Rξ tξ][Xw Yw Zw 1]T (3)
wherein λ isξIs the main camera projection scale coefficient; (u)ξ,vξ) Collecting pixel coordinates of feature points in a picture for a main camera; [ X ]wYw Zw 1]TIs a spatial three-dimensional point homogeneous coordinate; mξIs the internal parameter of the main camera; (R)ξ,tξ) An external reference matrix of a main camera;
s12: respectively solving first-order partial derivatives of the camera external parameters and the three-dimensional coordinate components of the space points through a formula (3) to obtain a collinearity error equation:
wherein xiξIs the deviation value of the main camera collinearity error; b isξ,CξAre respectively (u)ξ,vξ)TTo camera extrinsic parametersAnd the first partial derivative of the three-dimensional coordinate component of the spatial point; omegaξ,κξIs a rotation matrix RξA corresponding euler angle; t is txξ,tyξ,tzξIs an element of the main camera translation matrix, [ delta ]ξc δξp]TThe external parameter correction value and the three-dimensional coordinate correction value of the main camera are the step length of each iteration:
in the formula (5), Δ ωξ,Δκξ,Δtxξ,Δtyξ,ΔtzξAre first order partial derivatives of camera extrinsic parameters, respectivelyAn increment of (d); Δ X, Δ Y, Δ Z are the first partial derivatives of the three-dimensional coordinate components, respectively;
s13: matrix gammaξActual coordinates representing feature points and reprojection coordinates calculated by using collinearity equationThe difference between:
s14: according to the relative external reference relationship between the master camera and the slave camera, the collinearity equation of the eta camera is expressed as:
λη[uη vη 1]T＝Mη[RξηRξ tξη+Rξηtξ][Xw Yw Zw 1]T (7)
in formula (7), ληIs projecting the scale coefficients from the camera; (u)η,vη) The pixel coordinates of the characteristic point of the image under the eta camera are obtained; mηAn internal reference matrix of the eta camera; (R)ξη，tξη) Is the external parameter matrix of the slave camera;
s15: and obtaining the common linear error equation of the eta slave camera by the same method:
in the formula (8), xiηIs a deviation from camera co-linearity error; b isη,CηFirst-order partial derivatives of the external parameters and the spatial point coordinates from the camera image coordinates, respectively; [ delta ] isηc δηp]TRespectively obtaining an external reference correction value and a three-dimensional coordinate correction value of the slave camera; matrix arrayTo representThe difference between the actual coordinates of the feature points and the reprojection coordinates calculated by using a collinearity equation;
s16: let B be ═ Bξ B1 B2 … Bn]T,C＝[Cξ C1 C2 … Cη]TThen Jacobi matrix J ═ B C]Therefore, the normal matrix is represented as:
s17: iterative solution is carried out by using a Levenberg-Marquardt algorithm, and the established normal equation is expressed as follows:
wherein the content of the first and second substances,i is an identity matrix;
and obtaining a multi-camera network external parameter relative relation matrix.
3. The method for measuring the stereo deformation based on the parallel optimization of the CUDA architecture under the constraint of the novel correlation function as claimed in claim 1 or 2, wherein the method for establishing the novel correlation function of the speckle stereo registration in the step S02 comprises:
s21: and performing speckle image registration on the selected camera pair, wherein the speckle search range of the right camera imaging surface is represented as:
Fp(x,e)＝k0x+e(emin≤e≤emax) (11)
in the formula (11), x is the abscissa of an arbitrary point in the left camera object image, and k0Is the slope of the polar line equation, e is the intersection point coordinate of the polar line and the longitudinal axis of the coordinate system in the right camera image; e.g. of the typeminAnd emaxRespectively representing search pairs in the vicinity of epipolar lines in the right camera imageThe lower limit and the upper limit of the sub-pixel of the corresponding speckle point;
s22: floating in the upper and lower boundaries of polar line equation, seeking optimal correlation coefficient peak value, then substituting polar line equation obtained based on multi-camera joint constraint optimization into zero-mean normalized least square distance correlation function, and rewriting the generated novel correlation function expression into x'iAnd e constraint form:
in formula (12), f (x)i,yj) Is the reference image point (x) before deformationi,yj) The gray value of (d); meaning of k, g (x'i,y′j) Is the corresponding same name point (x ') in the deformed target image'i,y′j) Gray value of (f)mIs the average gray value of the sub-area before deformation; gmIs the average gray value of the sub-area after deformation; m is the registration sub-region center point to window boundary distance pixel size.
4. The parallel optimization stereo deformation measurement method based on the novel correlation function constraint CUDA framework as claimed in claim 3, wherein in step S02, parallel acceleration is implemented by using the CUDA framework when time sequence registration and stereo registration are performed on sampling scattered spots in different sequence images, and a Kernel function interface of the novel correlation function is constructed.
5. The CUDA architecture parallel optimization stereo deformation measurement method based on the novel correlation function constraint, according to claim 4, wherein the parallel acceleration method comprises:
s23: before starting a parallel program, importing the sub-pixel matching data of the speckle subarea into a constant memory, and importing the image data into a texture memory; after the speckle registration is finished, the calculation result is transmitted back to the host end; the parallel connection point of the mutual connection between the host end and the equipment end is a Kernel function;
s24: starting a Kernel function to search for the whole pixels and the sub-pixels of the speckle images according to the pre-configured thread blocks and the thread quantity, and simultaneously returning a flag bit for matching the iteratively converged speckle points; after all the thread calculations are completed, data obtained by speckle matching is transmitted back to a CPU memory area, and the sub-pixel searching process of scattered speckles is finished; during data transmission and access, the accumulation operation and access of the memory are carried out in a thread synchronization mode.
6. The method for measuring the parallel optimization stereo deformation of the CUDA architecture based on the novel correlation function constraint, according to claim 4, wherein the Kernel function processing step comprises:
inputting the number of threads, a left camera reference image, a right camera reference image and a target image sequence, defining a registration matrix space and a calibration matrix space;
downloading a pre-calculated data packet, k0,e,emin,emax；
Defining a Kernel function thread index, and initializing parallel computing parameters;
decomposing the novel correlation function, including: creating an inverse matrix space of a Hessian matrix, creating a Hessian matrix space, creating a Jacobi matrix space and creating a shape function parameter space;
iterative solution of speckle stereo registration of novel correlation function to obtain shape function vectorAnd corresponding target map astigmatic spot coordinates (x'i,y′j)；
Blocking the threads, ensuring the synchronization of the threads and releasing the memory space;
and outputting the time sequence registration data and the stereo registration data.
7. The method for measuring parallel optimization stereo deformation of a CUDA framework based on novel correlation function constraints as claimed in claim 5, wherein in step S24, the sub-pixel interpolation is performed in parallel, and 2M +1 threads are created to perform the sub-pixel interpolation parallel calculation on 2M +1 columns of speckle points respectively, wherein 2M +1 is the size of a speckle registration sub-area window.
8. A CUDA architecture parallel optimization stereo deformation measurement system based on novel correlation function constraint is characterized by comprising:
the multi-camera network external reference relative relationship matrix acquisition module acquires relative pose parameter relationships among the multiple cameras;
the correlation function building and registering module is used for building a novel correlation function of speckle three-dimensional registration based on a joint constraint relation between cameras, the novel correlation function limits the search of three-dimensional registration points between image pairs to an area near an epipolar line, and carries out time sequence registration and three-dimensional registration on sampling scattered spots in different sequence images;
and the reconstruction module is used for carrying out three-dimensional reconstruction on the space coordinates of the scattered spots.
9. A CUDA architecture parallel optimization stereo deformation measurement device based on novel correlation function constraints is characterized by comprising the CUDA architecture parallel optimization stereo deformation measurement method based on the novel correlation function constraints, which is claimed in any one of claims 1 to 7.
10. A storage medium, wherein the storage medium stores the CUDA architecture parallel optimization stereo deformation measurement method based on the novel correlation function constraint according to any one of claims 1 to 7.
Background
The three-dimensional deformation information measurement based on the vision measurement mode provides data reference and support for performance analysis, geometric deformation monitoring, bearing capacity evaluation and the like of the measured target. However, in the process of measuring the three-dimensional deformation, the scale of the operation data is large, for example, in order to improve the accuracy of the deformation measurement, a large number of speckle points are usually required to participate in the registration operation, so the scale of the data amount for registration is large, and meanwhile, as the number of frames of the image and the size of the image increase, the amount of the correlation operation in the speckle sub-area and the amount of the interpolation calculation of the sub-pixels both increase greatly, which consumes a lot of time and reduces the efficiency of the three-dimensional deformation measurement. The method has the advantages that the high-precision measurement of the three-dimensional deformation of the test piece based on the vision measurement mode needs to be realized at present, the cost of time and large calculation amount is needed, and the calculation speed is difficult to obtain the cross-order improvement only through the optimization and improvement on the algorithm. Therefore, the theoretical algorithm of deformation measurement needs to be improved, and meanwhile, an optimal parallel coupling operation mechanism needs to be built by combining the characteristics and complexity of the improved algorithm, so that the three-dimensional deformation measurement efficiency based on the visual measurement mode is improved.
Currently, there are two key problems in deformation measurement based on digital image correlation method. The three-dimensional speckle registration problem among different camera pairs is solved, if the image point coordinates corresponding to the maximum correlation coefficient are found by still adopting the traditional correlation function as the judgment standard of the similarity of the speckle image sub-regions among different cameras, the whole target image needs to be searched in the three-dimensional matching process, so that the searching efficiency is reduced, and due to the existence of parallax, the sub-pixel coordinates of the homonymous scattered spots directly obtained through two-dimensional correlation calculation have large deviation, and the reconstruction precision of the three-dimensional coordinates of the space of the subsequent scattered spots is seriously influenced. The second problem is a hardware acceleration mechanism problem related to speckle three-dimensional registration, researches on decomposition and parallel optimization of a speckle subregion registration algorithm based on a CUDA (compute unified device architecture) architecture are still few at present, and the methods generally adopt a hardware combination mode to accelerate a processing layer. For example, by analyzing the parallel operation problem of correlation coefficients in the digital image correlation deformation measurement, the foreign branch researchers adopt a platform mode combining a GPU and an FPGA to solve the problem that speckle matching is limited by operation resources; still other scholars seek a peak distribution of correlation function matches through multiple graphics processors, reducing overall run time. In addition, in some methods, a CUDA socket model is adopted to perform parallel accelerated calculation of speckle matching, but complexity and adaptability of an algorithm are not considered, so that the problem of limitation that memory overflow and the like are easily caused by improper cyclic splitting is easily caused. In China, professional research on the aspect is more scarce, and research shows that a small number of researchers realize accelerated research on speckle time sequence matching in deformation measurement through multithreading, and a certain acceleration effect is achieved. However, the above methods are directed at the acceleration of two-dimensional matching of speckle subregions in a visual measurement mode, do not consider and perform parallel coupling operation of speckle subregion three-dimensional registration between image pairs, and do not provide a specific implementation mode and an optimization scheme of speckle three-dimensional registration and time sequence registration parallel calculation based on a CUDA framework. The invention is achieved accordingly.
Disclosure of Invention
1. Objects of the invention
Aiming at the technical problems that in the three-dimensional deformation information measurement based on the visual measurement mode, the data volume scale of speckle point registration is large, simultaneously, along with the increase of image frame number and image size, the related calculation amount of a speckle subarea and the interpolation calculation amount of sub-pixels can be greatly increased, and the three-dimensional deformation measurement efficiency is greatly reduced in time consumption, the parallel optimization three-dimensional deformation measurement method and system based on the novel related function constraint CUDA framework are provided.
2. The technical scheme adopted by the invention
A CUDA architecture parallel optimization stereo deformation measurement method based on novel correlation function constraint comprises the following steps:
s01: obtaining a relative pose parameter relation among multiple cameras;
s02: establishing a novel correlation function of speckle three-dimensional registration based on a joint constraint relation between cameras, limiting the search of three-dimensional registration points between image pairs in an area near an epipolar line by the novel correlation function, and performing time sequence registration and three-dimensional registration on sampling scattered spots in images of different sequences;
s03: and carrying out three-dimensional reconstruction on the space coordinates of the scattered spots.
In an preferable technical solution, the method for obtaining the relative pose parameter relationship among the multiple cameras in step S01 includes:
s11: setting a measurement system to be composed of eta +1 cameras, wherein the main camera is xi, and the rest eta cameras are all slave cameras, and establishing a main camera collinearity equation based on a space point P according to a perspective projection principle:
λξ[uξ vξ 1]＝Mξ[Rξ tξ][Xw Yw Zw 1]T (3)
wherein λ isξIs the main camera projection scale coefficient; (u)ξ,vξ) Collecting pixel coordinates of feature points in a picture for a main camera; [ X ]w Yw Zw 1]TIs a spatial three-dimensional point homogeneous coordinate; mξIs the internal parameter of the main camera; (R)ξ,tξ) Of a main cameraAn external reference matrix;
s12: respectively solving first-order partial derivatives of the camera external parameters and the three-dimensional coordinate components of the space points through a formula (3) to obtain a collinearity error equation:
wherein xiξIs the deviation value of the main camera collinearity error; b isξ,CξAre respectively (u)ξ,vξ)TTo camera extrinsic parametersAnd the first partial derivative of the three-dimensional coordinate component of the spatial point;is a rotation matrix RξA corresponding euler angle; t is txξ,tyξ,tzξIs an element of the main camera translation matrix, [ delta ]ξc δξp]TThe external parameter correction value and the three-dimensional coordinate correction value of the main camera are the step length of each iteration:
in the formula (5), the reaction mixture is,respectively, the increment of the first-order partial derivative of the external parameter of the camera; Δ X, Δ Y, Δ Z are the first partial derivatives of the three-dimensional coordinate components, respectively;
s13: matrix gammaξActual coordinates representing feature points and reprojection coordinates calculated by using collinearity equationThe difference between:
s14: according to the relative external reference relationship between the master camera and the slave camera, the collinearity equation of the eta camera is expressed as:
λη[uη vη 1]T＝Mη[Rξη Rξ tξη+Rξηtξ][Xw Yw Zw 1]T (7)
in the formula (7), λ is a projection scale coefficient from the camera; (u)η,vη) The pixel coordinates of the characteristic point of the image under the eta camera are obtained; mηAn internal reference matrix of the eta camera; (R)ξη，tξη) Is the external parameter matrix of the slave camera;
s15: and obtaining the common linear error equation of the eta slave camera by the same method:
in the formula (8), xiηIs a deviation from camera co-linearity error; b isη,CηFirst-order partial derivatives of the external parameters and the spatial point coordinates from the camera image coordinates, respectively; [ delta ] isηc δηp]TRespectively obtaining an external reference correction value and a three-dimensional coordinate correction value of the slave camera; matrix gammaηExpressing the difference between the actual coordinates of the characteristic points and the reprojection coordinates calculated by using a collinearity equation;
s16: let B be ═ Bξ B1 B2 … Bη]T,C＝[Cξ C1 C2 … Cη]TThen Jacobi matrix J ═ B C]Therefore, the normal matrix is represented as:
s17: iterative solution is carried out by using a Levenberg-Marquardt algorithm, and the established normal equation is expressed as follows:
wherein γ ═ γξ γ1 γ2 … γη]T(ii) a I is an identity matrix;
and obtaining a multi-camera network external parameter relative relation matrix.
In a preferred embodiment, the method for establishing the novel correlation function of the speckle stereo registration in step S02 includes:
s21: and performing speckle image registration on the selected camera pair, wherein the speckle search range of the right camera imaging surface is represented as:
Fp(x,e)＝k0x+e(emin≤e≤emax) (11)
in the formula (11), x is the abscissa of an arbitrary point in the left camera object image, and k0Is the slope of the polar line equation, e is the intersection point coordinate of the polar line and the longitudinal axis of the coordinate system in the right camera image; e.g. of the typeminAnd emaxRespectively representing the lower limit and the upper limit of searching the sub-pixel of the corresponding scattered spot in the right camera image near the epipolar line;
s22: floating in the upper and lower bounds of an polar line equation, seeking an optimal correlation coefficient peak value, substituting the polar line equation obtained based on multi-camera joint constraint optimization into a zero-mean normalized least square distance correlation function, and rewriting the generated novel correlation function expression into a relation with xiConstraint forms of' and e:
in formula (12), f (x)i,yj) Is the reference image point (x) before deformationi,yj) The gray value of (d); meaning of k, g (x'i,y′j) Is the corresponding same name point (x ') in the deformed target image'i,y′j) Gray value of (f)mIs the average gray value of the sub-area before deformation;gmis the average gray value of the sub-area after deformation; m is the registration sub-region center point to window boundary distance pixel size.
In a preferred technical solution, when performing time sequence registration and stereo registration on sampling scattered spots in different sequence images in step S02, parallel acceleration is implemented based on a CUDA architecture, and a Kernel function interface of a novel correlation function is constructed.
In a preferred technical solution, the parallel acceleration method includes:
s23: before starting a parallel program, importing the sub-pixel matching data of the speckle subarea into a constant memory, and importing the image data into a texture memory; after the speckle registration is finished, the calculation result is transmitted back to the host end; the parallel connection point of the mutual connection between the host end and the equipment end is a Kernel function;
s24: starting a Kernel function to search for the whole pixels and the sub-pixels of the speckle images according to the pre-configured thread blocks and the thread quantity, and simultaneously returning a flag bit for matching the iteratively converged speckle points; after all the thread calculations are completed, data obtained by speckle matching is transmitted back to a CPU memory area, and the sub-pixel searching process of scattered speckles is finished; during data transmission and access, the accumulation operation and access of the memory are carried out in a thread synchronization mode.
In a preferred embodiment, the Kernel function processing step includes:
inputting the number of threads, a left camera reference image, a right camera reference image and a target image sequence, defining a registration matrix space and a calibration matrix space;
downloading a pre-calculated data packet, k0,e,emin,emax；
Defining a Kernel function thread index, and initializing parallel computing parameters;
decomposing the novel correlation function, including: creating an inverse matrix space of a Hessian matrix, creating a Hessian matrix space, creating a Jacobi matrix space and creating a shape function parameter space;
iterative solution of speckle stereo registration of novel correlation function to obtain shape function vectorAnd corresponding target map astigmatic spot coordinates (x'i,y′j)；
Blocking the threads, ensuring the synchronization of the threads and releasing the memory space;
and outputting the time sequence registration data and the stereo registration data.
In a preferred technical solution, in step S24, performing parallel computation on sub-pixel interpolation, and performing sub-pixel interpolation parallel computation on speckle points in 2M +1 columns by creating 2M +1 threads, respectively, where 2M +1 is a size of a speckle registration sub-area window.
The invention also discloses a CUDA architecture parallel optimization three-dimensional deformation measurement system based on novel correlation function constraint, which comprises the following steps:
the multi-camera network external reference relative relationship matrix acquisition module acquires relative pose parameter relationships among the multiple cameras;
the correlation function building and registering module is used for building a novel correlation function of speckle three-dimensional registration based on a joint constraint relation between cameras, the novel correlation function limits the search of three-dimensional registration points between image pairs to an area near an epipolar line, and carries out time sequence registration and three-dimensional registration on sampling scattered spots in different sequence images;
and the reconstruction module is used for carrying out three-dimensional reconstruction on the space coordinates of the scattered spots.
The invention also discloses a CUDA architecture parallel optimization three-dimensional deformation measuring device based on the novel correlation function constraint, which comprises the CUDA architecture parallel optimization three-dimensional deformation measuring method based on the novel correlation function constraint.
The invention also discloses a storage medium, and the storage medium stores the CUDA architecture parallel optimization three-dimensional deformation measurement method based on the novel correlation function constraint.
3. Advantageous effects adopted by the present invention
(1) Aiming at the problem that speckle registration consumes a long time in stereo deformation measurement, a CUDA heterogeneous parallel optimization stereo deformation measurement method based on novel correlation function constraint is provided. Firstly, a novel correlation function is established through the joint constraint relation among multiple cameras, the search of the stereo registration points among image pairs can be limited to the area near the epipolar line instead of the whole image, and therefore the search space is reduced. Then, a GPU parallel mechanism of a CUDA source program is compiled by adopting NVCC, and a parallel operation program design principle based on novel correlation function stereo registration and an optimization scheme of speckle subregion registration are provided. The technology in the invention is coupled with Kernel function to carry out decomposition of novel related function and design of parallel algorithm, and is not limited by heavy load function. The operation performance of the three-dimensional deformation measurement method reaches a better state from two aspects of three-dimensional registration algorithm optimization of the speckle subarea and coupling acceleration with hardware. Experimental results show that the CUDA heterogeneous parallel optimization algorithm based on the novel correlation function constraint respectively achieves the calculation acceleration ratios of speckle time sequence registration and three-dimensional registration of 20.39 times and 17.87 times under the condition of ensuring the measurement precision and stability, and provides technical bedding and reference for real-time measurement and output of three-dimensional deformation of a test piece.
(2) The deformation measurement method based on the visual pattern is improved from the theoretical algorithm, the complexity of the algorithm is reduced, and the accuracy and the stability of the speckle three-dimensional registration are improved, which is a key technology based on image texture matching. Then, considering the coupling factor of the improved algorithm and hardware acceleration, adopting a GPU parallel programming mode of compiling a CUDA source program by NVCC, solving the problem of barrier when the Mex script interacts with different programming languages, simultaneously being not limited by a heavy-load function, and realizing efficient operation of speckle image sub-area time sequence matching and three-dimensional matching in three-dimensional deformation information measurement. In addition, the independence of sub-pixel interpolation calculation is considered, a parallel calculation mode of the sub-pixel row and column interpolation of the speckle image is realized, and the sub-pixel search speed of speckle matching is further improved. Therefore, even if a computer system with a hardware configuration not higher is selected, a good acceleration effect can be still obtained based on the technology of the invention, accurate and stable speckle subregion registration data and a three-dimensional reconstruction cloud picture are obtained, the hardware cost and the development cost of a rapid measurement system based on visual deformation are reduced, and the rapid measurement and result reproduction system has practical engineering application value for realizing the rapid measurement and result reproduction of the three-dimensional deformation.
Drawings
FIG. 1 is a flow chart of a parallel optimization stereo deformation measurement method of a CUDA architecture based on novel correlation function constraints according to the present invention;
FIG. 2 is a schematic diagram of a speckle stereo registration search area optimized based on multi-camera network joint constraint in the embodiment;
FIG. 3 is a schematic block diagram of a CUDA architecture parallel optimization stereo deformation measurement system based on novel correlation function constraints;
FIG. 4 is a flowchart of the image sub-region matching parallel operation of the present invention;
FIG. 5 is a timing matching time-consuming comparison diagram according to the present invention;
FIG. 6 is a comparison diagram of time consumption for stereo matching according to the present invention;
FIG. 7 is a graph of displacement of a scattering spot measured at different times in accordance with the present invention;
FIG. 8 is a Z-direction shift cloud plot when t is 0.8s according to the present invention;
FIG. 9 is a Z-direction shift cloud plot for the present invention at t-1.5 s;
fig. 10 is a Z-direction displacement cloud chart when t is 2.0 s.
Detailed Description
The technical solution in the embodiment of the present invention is clearly and completely described below with reference to the accompanying drawings in the embodiment of the present invention. It is to be understood that the described embodiments are merely exemplary of the invention, and not restrictive of the full scope of the invention. All other embodiments, which can be derived by a person skilled in the art from the embodiments of the present invention without inventive step, are within the scope of the present invention.
The present invention will be described in further detail with reference to the accompanying drawings.
Example 1
As shown in fig. 1, a method for measuring a cubic deformation based on parallel optimization of a CUDA architecture based on a novel correlation function constraint includes the following steps:
s01: obtaining a relative pose parameter relation among multiple cameras;
s02: establishing a novel correlation function of speckle three-dimensional registration based on a joint constraint relation between cameras, limiting the search of three-dimensional registration points between image pairs in an area near an epipolar line by the novel correlation function, and performing time sequence registration and three-dimensional registration on sampling scattered spots in images of different sequences;
s03: and carrying out three-dimensional reconstruction on the space coordinates of the scattered spots.
In a preferred embodiment, the method for obtaining the relative pose parameter relationship among the multiple cameras in step S01 includes:
s11: setting a measurement system to be composed of eta +1 cameras, wherein the main camera is xi, and the rest eta cameras are all slave cameras, and establishing a main camera collinearity equation based on a space point P according to a perspective projection principle:
λξ[uξ vξ 1]＝Mξ[Rξ tξ][Xw Yw Zw 1]T (3)
wherein λ isξIs the main camera projection scale coefficient; (u)ξ,vξ) Collecting pixel coordinates of feature points in a picture for a main camera; [ X ]w Yw Zw 1]TIs a spatial three-dimensional point homogeneous coordinate; mξIs the internal parameter of the main camera; (R)ξ,tξ) An external reference matrix of a main camera;
s12: respectively solving first-order partial derivatives of the camera external parameters and the three-dimensional coordinate components of the space points through a formula (3) to obtain a collinearity error equation:
wherein xiξIs the deviation value of the main camera collinearity error; b isξ,CξAre respectively (u)ξ,vξ)TTo camera extrinsic parametersAnd the first partial derivative of the three-dimensional coordinate component of the spatial point;is a rotation matrix RξA corresponding euler angle; t is txξ,tyξ,tzξIs an element of the main camera translation matrix, [ delta ]ξc δξp]TThe external parameter correction value and the three-dimensional coordinate correction value of the main camera are the step length of each iteration:
in the formula (5), the reaction mixture is,respectively, the increment of the first-order partial derivative of the external parameter of the camera; Δ X, Δ Y, Δ Z are the first partial derivatives of the three-dimensional coordinate components, respectively;
s13: the matrix gamma represents the actual coordinates of the characteristic points and the reprojection coordinates calculated by the collinearity equationThe difference between:
s14: according to the relative external reference relationship between the master camera and the slave camera, the collinearity equation of the eta camera is expressed as:
λη[uη vη 1]T＝Mη[RξηRξ tξη+Rξηtξ][Xw Yw Zw 1]T (7)
in formula (7), ληIs projecting the scale coefficients from the camera; (u)η,vη) The pixel coordinates of the characteristic point of the image under the eta camera are obtained; mηAn internal reference matrix of the eta camera; (R)ξη，tξη) Is the external parameter matrix of the slave camera;
s15: and obtaining the common linear error equation of the eta slave camera by the same method:
in the formula (8), xiηIs a deviation from camera co-linearity error; b isη,CηFirst-order partial derivatives of the external parameters and the spatial point coordinates from the camera image coordinates, respectively; [ delta ] isηc δηp]TRespectively obtaining an external reference correction value and a three-dimensional coordinate correction value of the slave camera; matrix gammaηExpressing the difference between the actual coordinates of the characteristic points and the reprojection coordinates calculated by using a collinearity equation;
s16: let B be ═ Bξ B1 B2 … Bη]T,C＝[Cξ C1 C2 … Cη]TThen Jacobi matrix J ═ B C]Therefore, the normal matrix is represented as:
s17: iterative solution is carried out by using a Levenberg-Marquardt algorithm, and the established normal equation is expressed as follows:
wherein γ ═ γξ γ1 γ2 … γη]T(ii) a I is an identity matrix;
and obtaining a multi-camera network external parameter relative relation matrix.
In a preferred embodiment, the method for establishing the novel correlation function of the speckle stereo registration in step S02 includes:
s21: and performing speckle image registration on the selected camera pair, wherein the speckle search range of the right camera imaging surface is represented as:
Fp(x,e)＝k0x+e(emin≤e≤emax) (11)
in the formula (11), x is the abscissa of an arbitrary point in the left camera object image, and k0Is the slope of the polar line equation, e is the intersection point coordinate of the polar line and the longitudinal axis of the coordinate system in the right camera image; e.g. of the typeminAnd emaxRespectively representing the lower limit and the upper limit of searching the sub-pixel of the corresponding scattered spot in the right camera image near the epipolar line; the defined search area is shown in fig. 2.
S22: in FIG. 2 by [ e ]min,emax]Setting a gray search area, floating in the upper and lower boundaries of an epipolar equation, seeking an optimal correlation coefficient peak value, substituting the epipolar equation obtained based on multi-camera joint constraint optimization into a zero-mean normalized least-squares distance correlation function, and rewriting the generated novel correlation function expression into x'iAnd e constraint form:
in formula (12), f (x)i,yj) Is the reference image point (x) before deformationi,yj) The gray value of (d); meaning of k, g (x'i,y′j) Is the corresponding same name point (x ') in the deformed target image'i,y′j) Gray value of (f)mIs the average gray value of the sub-area before deformation; gmIs the average gray value of the sub-area after deformation; m is the registration sub-region center point to window boundary distance pixel size.
In a preferred embodiment, in step S02, when performing temporal registration and stereo registration on sampled scattered spots in different sequence images, a Kernel function interface based on CUDA architecture is used to implement parallel acceleration and construct a novel correlation function.
One specific implementation of the parallel acceleration method may be that the method includes the following steps:
s23: before starting a parallel program, importing the sub-pixel matching data of the speckle subarea into a constant memory, and importing the image data into a texture memory; after the speckle registration is finished, the calculation result is transmitted back to the host end; the parallel connection point of the mutual connection between the host end and the equipment end is a Kernel function;
s24: starting a Kernel function to search for the whole pixels and the sub-pixels of the speckle images according to the pre-configured thread blocks and the thread quantity, and simultaneously returning a flag bit for matching the iteratively converged speckle points; after all the thread calculations are completed, data obtained by speckle matching is transmitted back to a CPU memory area, and the sub-pixel searching process of scattered speckles is finished; during data transmission and access, the accumulation operation and access of the memory are carried out in a thread synchronization mode.
In a preferred embodiment, the Kernel function processing step includes:
inputting the number of threads, a left camera reference image, a right camera reference image and a target image sequence, defining a registration matrix space and a calibration matrix space;
downloading a pre-calculated data packet, k0,e,emin,emax；
Defining a Kernel function thread index, and initializing parallel computing parameters;
decomposing the novel correlation function, including: creating an inverse matrix space of a Hessian matrix, creating a Hessian matrix space, creating a Jacobi matrix space and creating a shape function parameter space;
iterative solution of speckle stereo registration of novel correlation function to obtain shape function vectorAnd corresponding target map astigmatic spot coordinates (x'i,y′j)；
Blocking the threads, ensuring the synchronization of the threads and releasing the memory space;
and outputting the time sequence registration data and the stereo registration data.
In a preferred embodiment, in step S24, the sub-pixel interpolation is performed in parallel, and 2M +1 threads are created to perform sub-pixel interpolation parallel calculation on speckle points of 2M +1 columns, respectively, where 2M +1 is the size of the speckle registration sub-area window.
As shown in fig. 3, the present invention also discloses a CUDA architecture parallel optimization stereo deformation measurement system based on the novel correlation function constraint, which includes:
the multi-camera network external reference relative relationship matrix acquisition module 10 acquires the relative pose parameter relationship among the multiple cameras;
a correlation function building and registering module 20, which builds a novel correlation function of speckle stereo registration based on the joint constraint relationship between cameras, the novel correlation function limits the search of stereo registration points between image pairs to the area near the epipolar line, and carries out time sequence registration and stereo registration on the sampled scattered spots in different sequence images;
and the reconstruction module 30 is used for carrying out three-dimensional reconstruction on the space coordinates of the scattered spots and fitting the cloud images before and after deformation.
The following describes the specific technical solution of the measurement system of the present invention in further detail through a complete embodiment. The workflow of the measurement system comprises the following steps:
the method comprises the following steps: and establishing a novel correlation function of speckle three-dimensional registration based on the joint constraint relation between the cameras.
Step 11: setting a measurement system to be composed of eta +1 cameras, recording a main camera as xi, and setting the rest eta cameras as slave cameras, and establishing a main camera collinearity equation based on a space point P according to a perspective projection principle:
λξ[uξ vξ 1]＝Mξ[Rξ tξ][Xw Yw Zw 1]T (3)
wherein λ isξIs the main camera projection scale coefficient; (u)ξ,vξ) Collecting pixel coordinates of feature points in a picture for a main camera; [ X ]w Yw Zw 1]TIs a spatial three-dimensional point homogeneous coordinate; mξIs the internal parameter of the main camera; (R)ξ,tξ) Is the external parameter matrix of the main camera.
Step 12: respectively solving first-order partial derivatives of the camera external parameters and the three-dimensional coordinate components of the space points through a formula (3) to obtain a collinearity error equation:
wherein the content of the first and second substances,
furthermore, xiξIs the deviation value of the main camera collinearity error; b isξ,CξAre respectively (u)ξ,vξ)TTo camera extrinsic parametersAnd the first partial derivative of the three-dimensional coordinate component of the spatial point;is a rotation matrix RξA corresponding euler angle; t is txξ,tyξ,tzξAre elements of the primary camera translation matrix. [ delta ] isξc δξp]TThe external parameter correction value and the three-dimensional coordinate correction value of the main camera are the step length of each iteration:
in the formula (5), the reaction mixture is,respectively, the increment of the first-order partial derivative of the external parameter of the camera; Δ X, Δ Y, Δ Z are the first partial derivatives of the three-dimensional coordinate components, respectively.
Step 13: matrix gammaξActual coordinates representing feature points and reprojection coordinates calculated by using collinearity equationThe difference between:
step 14: according to the relative external reference relationship between the master camera and the slave camera, the collinearity equation of the eta camera can be expressed as:
λη[uη vη 1]T＝Mη[RξηRξ tξη+Rξηtξ][Xw Yw Zw 1]T (7)
in formula (7), ληIs projecting the scale coefficients from the camera; (u)η,vη) The pixel coordinates of the characteristic point of the image under the eta camera are obtained; mηAn internal reference matrix of the eta camera; (R)ξη，tξη) Is the external parameter matrix of the slave camera.
Step 15: the collinearity error equation of the eta slave camera can be obtained by the same method:
in the formula (8), xiηIs a deviation from camera co-linearity error; b isη,CηFirst-order partial derivatives of the external parameters and the spatial point coordinates from the camera image coordinates, respectively; [ delta ] isηc δηp]TRespectively obtaining an external reference correction value and a three-dimensional coordinate correction value of the slave camera; matrix gammaηAnd representing the difference between the actual coordinates of the feature points and the reprojected coordinates calculated by using the collinearity equation.
Step 16: let B be ═ Bξ B1 B2 … Bη]T,C＝[Cξ C1 C2 … Cη]TThen Jacobi matrix J ═ B C]The normal matrix can be expressed as:
and step 17: iterative solution is carried out by using a Levenberg-Marquardt algorithm, and the established normal equation can be expressed as follows:
wherein γ ═ γξ γ1 γ2 … γη]T(ii) a I is the identity matrix.
Step 18: after the solution of the equation of the formula (10) is completed, an accurate matrix of the external reference relative relation of the multi-camera network can be obtained, and the three-dimensional registration search area of the homonymous scattered spots among the image pairs imaged by different cameras is limited in the adjacent range of epipolar constraint by combining the epipolar geometry principle. For speckle image registration for a selected camera pair, the speckle search range for the right camera imaging plane can be expressed as:
Fp(x,e)＝k0x+e(emin≤e≤emax) (11)
in the formula (11), x is the abscissa of an arbitrary point in the left camera object image, and k0Is the slope of the polar line equation, e is the intersection point coordinate of the polar line and the longitudinal axis of the coordinate system in the right camera image; e.g. of the typeminAnd emaxRespectively, the lower limit and the upper limit of searching the sub-pixel of the corresponding scattered spot in the right camera image near the epipolar line, and the defined search area is shown in fig. 2.
Step 19: in FIG. 2 by [ e ]min,emax]And setting a yellow search area, floating in the upper and lower boundaries of an epipolar equation, and seeking an optimal correlation coefficient peak value. And then substituting an epipolar line equation obtained based on multi-camera joint constraint optimization into a zero-mean normalized minimum square distance correlation function, and rewriting the generated novel correlation function expression into x'iAnd e constraint form:
in formula (12), f (x)i,yj) Is the reference image point (x) before deformationi,yj) The gray value of (d); meaning of k g (x'i,y′j) Is the corresponding same name point (x ') in the deformed target image'i,y′j) The gray value of (d). f. ofmIs the average gray value of the sub-area before deformation; gmIs the average gray value of the sub-area after deformation; m is the registration sub-region center point to window boundary distance pixel size. The established novel correlation function can reduce the search area of the speckle three-dimensional registration between the camera pairs, improve the search efficiency and improve the overall speed of three-dimensional deformation measurement to a certain extent.
Step two: a parallel implementation and optimization scheme of speckle subregion three-dimensional registration based on novel correlation function constraint by adopting a CUDA framework is designed, and a Kernel function interface of a novel correlation function is constructed at the same time.
Step 21: in three-dimensional deformation measurement, the calculation data of the correlation coefficient occupies most of the total calculation amount, and a parallel acceleration mechanism implemented based on a CUDA (compute unified device architecture) is used for realizing high-speed calculation of scattered spot matching. Speckle parallel matching flow chart based on CUDA (compute unified device architecture) as shown in FIG. 4, the sub-pixel search flow of the speckle image subarea comprises two parts of serial processing and parallel computing.
Step 22: because the data transmission between the CPU end and the video memory module can reduce the overall efficiency of GPU acceleration, the sub-pixel matching data of the speckle subarea is imported into a constant memory and the image data is imported into a texture memory before the parallel program is started; and after the speckle registration is finished, transmitting the calculation result back to the host end.
Step 23: the parallel connection points of the interconnection between the host end and the equipment end are Kernel functions, and the Kernel functions based on the novel correlation function speckle stereo matching are realized as shown in Table 1.
TABLE 1 speckle stereo matching Kernel function implementation based on novel correlation function
Step 24: and finally, starting a Kernel function to search for the integer pixels and the sub-pixels of the speckle images according to the pre-configured thread blocks and the thread quantity, and simultaneously returning the flag bits of the scattered spots matched with iterative convergence. And after all the thread calculations are completed, returning the data obtained by speckle matching to a CPU memory area, and ending the sub-pixel searching process of the scattered speckles. The running time of the parallel computing part is recorded by a timing statement in the program calling the mexw64 script file.
Step 25: parallel calculations are also performed on the sub-pixel interpolation. In the process of matching the speckle subareas of the target image, gray value interpolation operation at the sub-pixel positions among scattered spots is not influenced mutually, and 2M +1 threads are created to perform sub-pixel interpolation parallel calculation on the speckle points in the 2M +1 columns respectively. Where 2M +1 is the size of the speckle registration sub-region window. Thus, in the parallel operation mode, the serial mode (2M +1) can be achieved through 2M +1 times of calculation2The amount of secondary calculation.
Step 26: the parallel program design and the execution flow based on the CUDA architecture have great influence on the operation rate, so that the problem of the sequence of the accumulation operation and the access of the memory is solved in a thread synchronization mode when data transmission and access are involved in the design of the image speckle subregion matching optimization scheme: especially when the parallel program needs to perform frequent local space operation, the reasonable use of texture memory can enhance the overall performance of the program.
In another embodiment, the invention further discloses a CUDA architecture parallel optimization stereo deformation measurement device based on the novel correlation function constraint, which comprises the CUDA architecture parallel optimization stereo deformation measurement method based on the novel correlation function constraint.
In another embodiment, the invention further discloses a storage medium, where the method for measuring the parallel optimized stereo deformation of the CUDA architecture based on the novel correlation function constraint is stored.
The performance test results and analysis are as follows:
the computer used for testing is equipped with an Intel Xeon E5-2620 processor with a main frequency of 2.10GHz and an operating memory of 32GB, and a CUDAToolkit V8.0 is installed. The GPU is an NVIDIAQuadro P2000 series video card, 1024 stream processors are arranged in the GPU, the video memory bit width is 128 bits, and the video memory capacity is 5 GB. Firstly, compiling CUDA file codes and file targets at a terminal or through function calling by an NVCC compiler of an NVIDIAGPU is required to be carried out for subsequent Mex or Mexcuda calling. After the parallel computing code is written in the file with the cu suffix, the parallel computing code is compiled through a Mexcuda instruction to generate a Mex-Function file for subsequent calling to realize GPU parallel computing. A CU file may be created in the Visual Studio project and a CUDA file selected in the Build Customizaiton. In the actual measurement process, the adopted bi-cubic spline interpolation method takes a 7 multiplied by 7 window at the center of an area to be interpolated as a calculation unit, and solves the gray values at the sub-pixel positions of the rows and columns of the speckle image.
The same speckle ROI area is tested respectively through a digital image correlation matching algorithm based on CPU operation and CUDA heterogeneous parallel matching operation, and comparison on the calculation efficiency is carried out. Sequentially selecting 5, 20, 50, 100, 150, 200, 250 and 300 test points to perform time sequence sub-pixel search and stereo sub-pixel search respectively, wherein the matching running time is as shown in fig. 5 and fig. 6.
The results in fig. 5 and fig. 6 show that when the number of the test points participating in the matching is small, the time consumption of the serial matching method based on the CPU and the CUDA heterogeneous parallel operation is not much different, and the parallel acceleration is small. This is because the amount of data participating in the matching operation is small and is not yet in a state of being completely saturated in the CPU calculation capability; when the GPU-side kernel function is started to perform parallel operation, the calculated data needs to be transmitted from the host end to the equipment end, the result needs to be transmitted back to the host end after the operation is completed, and when the data scale is small, the parallel calculation efficiency is not enough to cover the delay of memory access and data transmission, so that the parallel acceleration effect is not obvious. However, as the number of matching points is increased, the time consumption of the conventional serial operation mode is increased more, and the time consumption of the parallel operation is increased more slowly. The proportion of non-computation time consumed on a GPU to the total algorithm execution time is smaller and smaller, and for time-sequence matching time consumption and stereo matching time consumption, the heterogeneous parallel computation speed based on the CUDA is far higher than that of a traditional serial computation mode. And the acceleration effect is more obvious along with the continuous increase of the scale of the matching points. When the number of matching points is 300, the time sequence matching acceleration ratio of the speckle subarea is 20.39 times, the stereo matching acceleration ratio is 17.87 times, and the acceleration effect is obvious.
And then verifying whether the parallel acceleration optimization scheme influences the deformation measurement precision or not, and building a measurement system hardware structure, wherein the deformation of the test piece is realized by external force, the left side of the test piece is rigidly connected to a target simulator sliding block through a bolt, and the right side of the test piece is fixed on a table body, so that the dynamic deformation of the wing surface is simulated through the movement of the target simulator. The sampling frequency of the camera was 5Hz and the exposure time was 7.07 ms. The moving speed of the target simulator is 2cm/s, and the moving time is 2 s.
Speckle points on a certain line are selected at equal intervals to carry out spatial three-dimensional reconstruction, and because the target simulator carries out motion in the horizontal direction, the deformation in the depth direction is calculated in a targeted manner, and a curve in the Z direction is obtained and is shown in FIG. 7. In addition, the total station is used for sampling and verifying the deformation of the test piece, and the obtained calibration data is consistent with the measurement result in fig. 7, so that the precision and the stability of the measurement result can be ensured on the basis of realizing the speed deformation measurement.
For convenience of analysis and reproduction of actual effects, deformation amounts of the speckle pattern in the Z direction at the time when t is 0.8s, t is 1.5s, and t is 2.0s in the dynamic captured image are fitted, and then cloud images of fitting results at different times are obtained as shown in fig. 8, 9, and 10. It can be seen that the displacement cloud in the Z-direction is consistent with the actual deformation.
In summary, the invention aims at the problems that the realization of high-precision deformation measurement needs time and large calculation amount as cost, the calculation speed is difficult to obtain across-order improvement only through optimization and improvement on the algorithm, and the like because the data scale is large in the process of three-dimensional deformation measurement, the operation performance of the three-dimensional deformation measurement method reaches a better state through the two aspects of three-dimensional registration algorithm optimization of speckle subregions and coupling acceleration with hardware, the efficient parallel operation of time sequence matching and three-dimensional matching of the speckle image subregions in the three-dimensional deformation measurement is realized, and the high-precision and quick measurement functions of the three-dimensional deformation based on the visual non-contact measurement mode are further realized.
The above description is only for the preferred embodiment of the present invention, but the scope of the present invention is not limited thereto, and any changes or substitutions that can be easily conceived by those skilled in the art within the technical scope of the present invention are included in the scope of the present invention. Therefore, the protection scope of the present invention shall be subject to the protection scope of the claims.
