Three-dimensional digital adaptive modeling method for cultural relics

Abstract

The invention relates to the technical field of cultural relic modeling, solves the technical problems of low operation efficiency and low modeling accuracy in the prior art, and particularly relates to a cultural relic three-dimensional digital adaptive modeling method which comprises the following steps: S1, obtaining three-dimensional point cloud data of a cultural relic to generate a point cloud data set, determining a neighborhood point set based on the point cloud data set and obtaining a plurality of scale features through feature decomposition; according to the method, the complexity of algorithm operation can be effectively reduced through global implicit values, the error rate of the point cloud data can be effectively reduced through axial alignment of the point cloud data, the processing time is shortened by using contour surfaces, and the processing efficiency is improved. And finally, through smoothing processing of a three-dimensional smoothing operator, carving details of the cultural relics can be better reserved, the three-dimensional model is more accurate, the operation efficiency can be improved, and high accuracy can be achieved.

Description

Technical Field

[0001] The present invention relates to the technical field of cultural relic modeling, and in particular to a three-dimensional digital adaptive modeling method for cultural relics. Background Art

[0002] Cultural relic modeling generally refers to the process of using computer technology and related software tools to perform three-dimensional digital reconstruction of cultural relics. By collecting information such as the geometric shape, surface texture, and color of the cultural relics, a three-dimensional digital model of the cultural relics is constructed to achieve digital protection, display, research, and dissemination of the cultural relics. Existing cultural relic three-dimensional modeling methods generally use image measurement modeling or scanning modeling to construct them. By performing image acquisition and feature extraction on multiple images, the two-dimensional to three-dimensional algorithm is used to convert the plane image into a three-dimensional model, or the point cloud data of the cultural relic is directly obtained and extracted based on the three-dimensional features. However, these methods have various defects, including large data requirements, huge amount of calculations, low efficiency, and low model processing accuracy. It is difficult to accurately construct the cultural relic model, and the user experience is unsatisfactory. Summary of the Invention

[0003] In response to the shortcomings of the existing technology, the present invention provides a three-dimensional digital adaptive modeling method for cultural relics, which solves the technical problems of low computing efficiency and low modeling accuracy in the existing technology, and achieves the purpose of greatly improving computing efficiency and enhancing modeling accuracy.

[0004] To solve the above technical problems, the present invention provides the following technical solution: a method for three-dimensional digital adaptive modeling of cultural relics, the method comprising the following steps:

[0005] S1. Obtain the 3D point cloud data of the cultural relics to generate point cloud dataset B s , based on point cloud dataset B s Determine the neighborhood point set P k And obtain multiple scale features λ through feature decomposition n s ;

[0006] S2, according to the scale feature λ n s Calculate the adaptive filter threshold T for removing noise points y , and obtain the preprocessing point set;

[0007] S3. Select any point in the preprocessing point set as the origin to construct a three-dimensional rectangular coordinate system. In the three-dimensional rectangular coordinate system, the surface fitting equation N for calculating the partial derivative is obtained by the surface fitting method. h ;

[0008] S4. Calculate the normal vector based on the partial derivative Based on normal vector Calculate the first curvature threshold Y for simplifying the point cloud data q and the second curvature threshold Y w , and obtain a simplified point set;

[0009] S5, calculate the eigenvalue sequence (λ based on the simplified point set max ,λ max-1 …λ1,λ0), and through the eigenvalue sequence (λ max ,λ max-1 …λ1,λ0) determine the three axis vectors {u1,u2,u3};

[0010] S6. Construct a directed bounding box model based on the three axis vectors {u1, u2, u3} and calculate the adaptive Gaussian kernel φ j (D), based on the adaptive Gaussian kernel φ j (D) Calculate global implicit values

[0011] S7, according to the global implicit value Calculate the isosurface eigenvalue η used to construct the 3D digital model of cultural relics r , and the three-dimensional digital model of the cultural relics is smoothed to obtain an optimized three-dimensional model.

[0012] Furthermore, in step S1, the specific implementation steps are as follows:

[0013] S11, the three-dimensional point cloud data of the cultural relics is based on the level radius r s Divide into cube grids with multiple levels, and build an octree structure based on the level radius. The expression is:

[0014] r s =2 s ·δ

[0015] Among them, r s represents the radius of the sth level, where s=0,1,2...s;

[0016] S12. Define point cloud dataset B s ={p i ,i=1,2...i},p i Represents the i-th data point, and selects a point p in the point cloud dataset a As the center point, get point p a The neighborhood point set P k ={p l ,l=1,2...l};

[0017] S13, according to the neighborhood point set P k Calculate the covariance matrix value X c , the calculation formula is:

[0018]

[0019] Where K represents p l the number of Represents the mean value of the neighborhood point coordinates;

[0020] S14, the covariance matrix value X c Perform eigenvalue decomposition to obtain multiple scale features λ n s , where λ n s Represents the feature of the nth scale in the sth level, and the expression is:

[0021] {λ0 s ,λ1 s ,...λ n s}=eig(X c )

[0022] Among them, eig represents the covariance matrix value X c Perform eigendecomposition.

[0023] Furthermore, in step S2, the specific implementation steps are as follows:

[0024] S21. Based on multiple scale features λ n s Calculate the multi-scale curvature characteristic ratio Tz b , the calculation formula is:

[0025]

[0026] Among them, λ k s Represented in multiple scale features λ n s The kth scale feature, N represents the scale feature λ n s the number of

[0027] S22, according to the multi-scale curvature characteristic ratio Tz b Calculate the fusion scale eigenvalue Rh c , the calculation formula is:

[0028]

[0029] Among them, Rh c represents the cth fusion scale eigenvalue;

[0030] S23, according to the neighborhood point set P k Calculate the first statistic μ s , the calculation formula is:

[0031]

[0032] Among them, μ s represents the sth first statistic;

[0033] S24, according to the first statistic μ s Calculate the second statistic α t , the calculation formula is:

[0034]

[0035] Among them, α t represents the t-th second statistic;

[0036] S25, according to the first statistic μ s and the second statistic α t Calculate the adaptive filtering threshold T y , the calculation formula is:

[0037] T y =μ s +1.5·α t ·(1-Rh c )

[0038] Among them, T y represents the yth adaptive filtering threshold;

[0039] S26, according to the adaptive filtering threshold T y Remove noise;

[0040] like Then point p l Remove the noise and mark the removal point t T ;

[0041] like Then point p l is a normal data point;

[0042] S27, repeat step S26 until the point cloud data set is detected and a preliminary pre-processed point set is obtained;

[0043] S28, according to the elimination point t T Calculate supplement point B bc , the calculation formula is:

[0044]

[0045] Among them, t T+1 and t T-1 Respectively represent and eliminate point t T Two adjacent points;

[0046] S29, add point B bc Add to the preliminary preprocessing point set and obtain the preprocessing point set.

[0047] Furthermore, in step S3, the specific implementation steps are as follows:

[0048] S31. Select a point p in the preprocessing point set d Establish a three-dimensional rectangular coordinate system for the origin and define any point p in the preprocessing point set i The coordinates of (x i ,y i ,z i ), and construct the surface fitting equation N of the preprocessing point set based on the three-dimensional rectangular coordinate system h , the expression is:

[0049] N h =ax 2 +bxy+cy 2

[0050] Among them, a, b and c represent the surface fitting equation N h The coefficient of

[0051] S32, according to the surface fitting equation N h Calculate the least squares value Q f , the calculation formula is:

[0052]

[0053] Where H represents the number of preprocessing points, z i Indicates the z-axis coordinate value of the preprocessing point;

[0054] S33, based on the least squares value Q f Calculate the values ​​of coefficients a, b and c respectively, and the calculation formula is:

[0055]

[0056] Among them, Q f represents the fth least squares value;

[0057] S34, fitting the surface equation N h Convert to surface function f e (x,y,z), the expression is:

[0058]

[0059] Among them, f e (x,y,z) represents the e-th surface function;

[0060] S35, according to the surface function fe (x, y, z) calculates the partial derivative, which includes the first partial derivative f x1 , the second partial derivative f y2 , the third partial derivative f xy3 、The fourth partial derivative f xx4 and the fifth partial derivative f yy5 , the calculation formula is:

[0061]

[0062] in, Represents the surface function f e Partial derivatives of (x,y,z).

[0063] Furthermore, in step S4, the specific implementation steps are as follows:

[0064] S41, according to the first partial derivative f in the partial derivative x1 and the second partial derivative f y2 Calculating the normal vector The calculation formula is:

[0065]

[0066] in, Represents the surface function f e The normal vector of (x,y,z);

[0067] S42, according to the first partial derivative f x1 and the second partial derivative f y2 Calculate the curvature value E g 、F h and G j , the calculation formula is:

[0068] E g =f x1 ·f x1

[0069] F h =f x1 ·f y2

[0070] G j =f y2 ·f y2

[0071] Among them, E g 、F h and G j Represent the g, h and j-th curvature values ​​respectively;

[0072] S43, according to the third partial derivative f xy3 、The fourth partial derivative f xx4and the fifth partial derivative f yy5 Calculate the depth value L k 、M l and N m , the calculation formula is:

[0073]

[0074] Among them, L k 、M l and N m Represent the kth, lth and mth depth values ​​respectively;

[0075] S44. Calculate the Gaussian curvature K used to screen surface features n , mean curvature H pj and global curvature

[0076] S45. Repeat steps S41-S44 to traverse all points in the preprocessing point set;

[0077] S46. Calculate the first curvature threshold Y q and the second curvature threshold Y w , the calculation formula is:

[0078]

[0079] Among them, Y q represents the qth curvature threshold;

[0080] S47, according to the curvature threshold Y q Complete the screening of the preprocessing point set;

[0081] If K n ≥Y w and It is marked as a reserved point;

[0082] If K n <Y w and These are points that can be streamlined and eliminated;

[0083] If K n <Y w or It is marked as a reserved point;

[0084] S48. Repeat step S47 to simplify all pre-processed point sets to obtain a simplified point set and end.

[0085] Furthermore, the Gaussian curvature K n , mean curvature H pj and global curvature The calculation formula is:

[0086]

[0087] Where H represents the mean curvature H pj The number of

[0088] Furthermore, in step S5, the specific implementation steps are as follows:

[0089] S51. Define the streamlined point pj in the streamlined point set t The coordinates are (A x ,B y ,C z ) and the number is D, calculate all the simplified points pj t The center of mass O zx The coordinates are calculated as follows:

[0090]

[0091] Among them, O zx represents the xth centroid coordinate;

[0092] S52, randomly select a simplified point pj t As a reference point pc t , calculate the reference point and other simplified points pj t The mean distance between the points L d , the calculation formula is:

[0093]

[0094] Among them, L d represents the mean value of the dth point distance;

[0095] S53, according to the point distance mean L d Calculate the weight coefficient w h , the calculation formula is:

[0096]

[0097] Among them, w h represents the hth weight coefficient;

[0098] S54, based on the weight coefficient w h Calculate the covariance matrix Xc used to determine the axial volume of the oriented bounding box p , the calculation formula is:

[0099]

[0100] Among them, s represents the sum of weight coefficients, and m represents the simplified point pj t The weighted sum of

[0101] S55, according to the covariance matrix Xc p Calculate the characteristic equation and obtain multiple eigenvalues ​​λ E ,E=1,2...E, the calculation formula is:

[0102]

[0103] det(Xc p -λ E I)=0

[0104] Where det represents Xc p -λ E The determinant of I, where I represents the identity matrix, represents the eigenvector;

[0105] S56, multiple eigenvalues ​​λ E Arrange in ascending order to obtain the eigenvalue sequence (λ max ,λ max-1 ...λ1,λ0), the eigenvalue λ max ,λ max-1 and λ max-2 Set as the axis eigenvalue, and the eigenvalue λ max ,λ max-1 and λ max-2 The corresponding eigenvectors are set as the three axis vectors {u1,u2,u3} of the directed bounding box.

[0106] Furthermore, in step S6, the specific implementation steps are as follows:

[0107] S61. Construct a directed bounding box model based on the axis vector and the simplified point set, and calculate the alignment point Dq used to align the simplified points. u , the calculation formula is:

[0108] Dq u =R d (pj t -O zx )

[0109]

[0110] Among them, R d Represents the transposed value of the d-th axis vector;

[0111] S62, along the axis vector {u1, u2, u3}, the directed bounding box space is defined as the network area, according to the point distance mean L d Calculate the partition side length Ψ b , the calculation formula is:

[0112] Ψ b =1.5·median(L d )

[0113] Among them, b represents the length of the b-th partition side;

[0114] S63, based on the division of side length Ψ b Divide the network area into multiple sub-areas S c , calculate the approval value σ for each sub-region h , the calculation formula is:

[0115] σ h =median(L d inS c )

[0116] Among them, σ h represents the hth approved value;

[0117] S64, according to the approved value σ h Calculate the adaptive Gaussian kernel φ j (D), calculated as:

[0118]

[0119] Among them, φ j (D) represents the jth adaptive Gaussian kernel, Dq k and Dq k-1 denote the kth and k-1th alignment points respectively;

[0120] S65. Obtain the center point O of each sub-region by using the model coordinate method. D and half-length mean Bc z , calculate the global implicit value The calculation formula is:

[0121]

[0122] in, Represents the g-th normal vector.

[0123] Furthermore, in step S7, the specific implementation steps are as follows:

[0124] S71. According to the global implicit value Calculate the isosurface eigenvalue η r , the calculation formula is:

[0125]

[0126] Among them, η r represents the rth isosurface eigenvalue;

[0127] S72, according to the isosurface eigenvalue η rConstruct a three-dimensional digital model of the cultural relic through the cube reconstruction method;

[0128] S73. Define the area of ​​each facet in the 3D digital model of the cultural relic as A. s , the normal vector is Calculating vertex normals The calculation formula is:

[0129]

[0130] Where P represents the number of patches;

[0131] S74, according to the vertex normal vector Calculate the smoothing strength Qp h , the calculation formula is:

[0132]

[0133] Among them, u f represents the f-th axis vector;

[0134] S75. Define the vertex d of the three-dimensional digital model of the cultural relic in , the vertex of 1 ring area is N id , calculate the vertex offset Ld of the three-dimensional digital model of the cultural relic v , the calculation formula is:

[0135]

[0136] Where I represents the number of vertices in the 1-ring domain;

[0137] S76, according to the vertex offset Ld v Calculate the model smoothing value L ph , the calculation formula is:

[0138]

[0139] W(x) i =diag(Qp h )

[0140] Among them, Z represents the vertex offset Ld v the number of

[0141] S77, smoothing value L according to the model ph The surface of the three-dimensional digital model of the cultural relics is smoothed to obtain an optimized three-dimensional model.

[0142] By means of the above technical solution, the present invention provides a method for three-dimensional digital adaptive modeling of cultural relics, which has at least the following beneficial effects:

[0143] 1. The present invention can make point cloud data more accurate through multi-scale feature processing and adaptive filtering of three-dimensional point cloud data. The octree structure can accelerate three-dimensional data search, quantify the degree of surface curvature, and reduce the accidental deletion of feature points through the curvature threshold method. The use of adaptive threshold can also improve the robustness of the algorithm.

[0144] 2. The present invention constructs a surface fitting equation to extract features of the surface based on the fitted surface, and determines the simplification of the point cloud data based on the eigenvalues, thereby completing the simplification of the point cloud data and constructing a directed bounding box model to prepare for the construction of a three-dimensional model in the subsequent steps. The simplification steps of the point cloud data are fast and efficient. The simplified model can reduce the number of calculation steps and data point operations, greatly improving the algorithm's computational efficiency while ensuring the accuracy of the vertices.

[0145] 3. The present invention can effectively reduce the complexity of algorithm operations through global implicit values, effectively reduce the error rate of point cloud data through axial alignment of point cloud data, reduce processing time by using isosurfaces, and finally, through smoothing processing of a three-dimensional smoothing operator, better preserve the carving details of cultural relics, make the three-dimensional model more accurate, and improve computing efficiency while achieving higher accuracy. BRIEF DESCRIPTION OF THE DRAWINGS

[0146] The drawings described herein are used to provide a further understanding of the present application and constitute a part of the present application. The illustrative embodiments of the present application and their descriptions are used to explain the present application and do not constitute an improper limitation on the present application. In the drawings:

[0147] Figure 1 The present invention is a flowchart of a method for three-dimensional digital adaptive modeling of cultural relics. DETAILED DESCRIPTION

[0148] To make the above-mentioned objectives, features, and advantages of the present invention more clearly understood, the present invention is further described below in detail with reference to the accompanying drawings and specific embodiments. This will enable a full understanding of how this application uses technical means to solve technical problems and achieve technical effects, and to implement the invention accordingly.

[0149] Due to the technical problems of low computational efficiency and low modeling accuracy in the existing technology, this embodiment proposes a three-dimensional digital adaptive modeling method for cultural relics, which can greatly improve computational efficiency and enhance modeling accuracy. Figure 1 As shown, the method includes the following steps:

[0150] S1. Obtain the 3D point cloud data of the cultural relics to generate point cloud dataset B s , based on point cloud dataset B s Determine the neighborhood point set P kAnd obtain multiple scale features λ through feature decomposition n s In the process of acquiring point cloud data, due to the accuracy of the equipment and the influence of dim ambient light, it is easy to cause inaccurate point cloud data. In order to solve this problem, the specific implementation steps are as follows:

[0151] S11, the three-dimensional point cloud data of the cultural relics is based on the level radius r s Divide into cube grids with multiple levels, and build an octree structure based on the level radius. The expression is:

[0152] r s =2 s ·δ

[0153] Among them, r s Represents the radius of the s-th level, and s = 0, 1, 2...s; the octree structure is a commonly used three-dimensional data structure that can improve the processing efficiency of three-dimensional data preprocessing and has good accuracy, which will not be described in detail here.

[0154] S12. Define point cloud dataset B s ={p i ,i=1,2...i},p i Represents the i-th data point, and selects a point p in the point cloud dataset a As the center point, get point p a Neighborhood point set P k ={p l ,l=1,2...l};p l Represents the data points in the neighborhood point set, which is the radius r of the center point level s Data points within the range.

[0155] S13, according to the neighborhood point set P k Calculate the covariance matrix value X c , the calculation formula is:

[0156]

[0157] Where K represents p l the number of Represents the mean value of the neighborhood point coordinates;

[0158] S14, the covariance matrix value X c Perform eigenvalue decomposition to obtain multiple scale features λ n s , where λ n s Represents the feature of the nth scale in the sth level, and the expression is:

[0159] {λ0s ,λ1 s ,...λ n s}=eig(X c )

[0160] Among them, eig represents the covariance matrix value X c By performing feature decomposition and performing multi-scale feature processing and adaptive filtering on the three-dimensional point cloud data through the algorithm, the point cloud data can be made more accurate. The octree structure can accelerate the three-dimensional data search, quantify the degree of surface curvature, and reduce the accidental deletion of feature points through the curvature threshold method. The use of adaptive threshold can improve the robustness of the algorithm.

[0161] S2, according to the scale feature λ n s Calculate the adaptive filtering threshold T for removing noise points y , and obtain the preprocessed point set; based on step S1, perform adaptive filtering on the point cloud data. In order to solve this problem, the detailed steps are as follows:

[0162] S21. Based on multiple scale features λ n s Calculate the multi-scale curvature characteristic ratio Tz b , the calculation formula is:

[0163]

[0164] Among them, λ k s Represented in multiple scale features λ n s The kth scale feature, N represents the scale feature λ n s the number of

[0165] S22, according to the multi-scale curvature characteristic ratio Tz b Calculate the fusion scale eigenvalue Rh c , the calculation formula is:

[0166]

[0167] Among them, Rh c represents the cth fusion scale eigenvalue;

[0168] S23, according to the neighborhood point set P k Calculate the first statistic μ s , the calculation formula is:

[0169]

[0170] Among them, μ s represents the sth first statistic, K represents the point p in the neighborhood point set l the number of

[0171] S24, according to the first statistic μ s Calculate the second statistic α t , the calculation formula is:

[0172]

[0173] Among them, α t represents the t-th second statistic;

[0174] S25. According to the first statistic μ s and the second statistic α t Calculate the adaptive filtering threshold T y , the calculation formula is:

[0175] T y =μ s +1.5·α t ·(1-Rh c )

[0176] Among them, T y represents the yth adaptive filtering threshold;

[0177] S26, according to the adaptive filtering threshold T y Remove noise;

[0178] like Then point p l Remove the noise and mark the removal point t T ;

[0179] like Then point p l is a normal data point;

[0180] S27, repeat step S26 until the point cloud data set is detected and a preliminary pre-processed point set is obtained;

[0181] S28, according to the elimination point t T Calculate supplement point B bc , the calculation formula is:

[0182]

[0183] Among them, t T+1 and t T-1 Represent and eliminate point t respectively T Two adjacent points;

[0184] S29, add point B bcIt is added to the preliminary preprocessing point set and the preprocessing point set is obtained. Through multi-scale feature processing and adaptive filtering of three-dimensional point cloud data, the point cloud data can be made more accurate. The octree structure can accelerate the three-dimensional data search, quantify the degree of surface curvature, and reduce the accidental deletion of feature points through the curvature threshold method. The use of adaptive threshold can improve the robustness of the algorithm.

[0185] S3, construct a three-dimensional rectangular coordinate system based on any point in the preprocessing point set, and obtain the surface fitting equation N according to the points in the three-dimensional rectangular coordinate system h , according to the surface fitting equation N h Calculate the partial derivatives; Due to the large amount of point cloud data in the preprocessing point set, the step-by-step calculation steps are very complicated, which seriously affects the efficiency of the calculation. In order to solve this problem, the detailed implementation steps are as follows:

[0186] S31. Select a point p in the preprocessing point set. d As the origin, a three-dimensional rectangular coordinate system is established, where any point p in the preprocessing point set is i The coordinates of (x i ,y i ,z i ), and construct the surface fitting equation N of the preprocessing point set based on the three-dimensional rectangular coordinate system h , the expression is:

[0187] N h =ax 2 +bxy+cy 2

[0188] Among them, a, b and c represent the surface fitting equation N h The coefficient of

[0189] S32, according to the surface fitting equation N h Calculate the least squares value Q f , the calculation formula is:

[0190]

[0191] Where H represents the number of preprocessing points, z i Indicates the z-axis coordinate value of the preprocessing point;

[0192] S33, based on the least squares value Q f Calculate the values ​​of coefficients a, b and c respectively, and the calculation formula is:

[0193]

[0194] Among them, Q f represents the fth least squares value;

[0195] S34, fitting the surface equation N h Convert to surface function f e (x,y,z), the expression is:

[0196]

[0197] Among them, f e (x,y,z) represents the e-th surface function;

[0198] S35, according to the surface function f e (x, y, z) calculates the partial derivative, which includes the first partial derivative f x1 , the second partial derivative f y2 , the third partial derivative f xy3 、The fourth partial derivative f xx4 and the fifth partial derivative f yy5 , the calculation formula is:

[0199]

[0200]

[0201] in, Represents the surface function f e The partial derivatives of (x, y, z) can be used to construct surface fitting equations to extract surface features based on the fitted surface. The simplification of point cloud data can be determined based on the eigenvalues, completing the simplification of point cloud data. At the same time, a directed bounding box model is constructed to prepare for the subsequent three-dimensional model construction. The point cloud data simplification steps are fast and efficient. The streamlined model can reduce the number of calculation steps and data point operations, greatly improving the algorithm's computational efficiency while ensuring vertex accuracy.

[0202] S4. Calculate the normal vector based on the partial derivative Based on normal vector Calculate the first curvature threshold Y for simplifying the point cloud data q and the second curvature threshold Y w , and obtain a simplified point set; based on step S3, it is necessary to further obtain the characteristics of the three-dimensional surface through partial derivatives, and simplify the point cloud data based on this feature. The specific solution steps are as follows:

[0203] S41, according to the first partial derivative f in the partial derivative x1 and the second partial derivative f y2 Calculating the normal vector The calculation formula is:

[0204]

[0205] in, Represents the surface function f e The normal vector of (x,y,z);

[0206] S42, according to the first partial derivative f x1 and the second partial derivative f y2 Calculate the curvature value E g 、F h and G j , the calculation formula is:

[0207] E g =f x1 ·f x1

[0208] F h =f x1 ·f y2

[0209] G j =f y2 ·f y2

[0210] Among them, E g 、F h and G j Represent the g, h and j-th curvature values ​​respectively;

[0211] S43, according to the third partial derivative f xy3 、The fourth partial derivative f xx4 and the fifth partial derivative f yy5 Calculate the depth value L k 、M l and N m , the calculation formula is:

[0212]

[0213] Among them, L k 、M l and N m Represent the kth, lth and mth depth values ​​respectively;

[0214] S44. Calculate the Gaussian curvature K used to screen surface features n , mean curvature H pj and global curvature The calculation formula is:

[0215]

[0216] Where H represents the mean curvature H pj the number of

[0217] S45. Repeat steps S41-S44 to traverse all points in the preprocessing point set;

[0218] S46. Calculate the first curvature threshold Y q and the second curvature threshold Y w , the calculation formula is:

[0219]

[0220] Among them, Y q represents the qth curvature threshold;

[0221] S47, according to the curvature threshold Y q Complete the screening of the preprocessing point set;

[0222] If K n ≥Y w and It is marked as a reserved point;

[0223] If K n <Y w and These are points that can be streamlined and eliminated;

[0224] If K n <Y w or The points are marked as retained points. Retained points generally indicate that the surface features are very obvious, and they are points on the boundary or arc surface. The points that can be simplified are generally points on the plane or relatively flat arc surface. These points do not need to be fully calculated, and only a very small amount is needed to complete them, so they can be eliminated.

[0225] S48. Repeat step S47 to streamline all pre-processed point sets to obtain a streamlined point set and end. By constructing a surface fitting equation, it is possible to extract features of the surface based on the fitted surface, and determine the simplification of the point cloud data based on the eigenvalues, thereby completing the simplification of the point cloud data and constructing a directed bounding box model to prepare for the construction of a three-dimensional model in subsequent steps. The step of streamlining point cloud data is fast and efficient. The streamlined model can reduce the number of calculation steps and data point operations, greatly improving the algorithm's computational efficiency while ensuring vertex accuracy.

[0226] S5. Calculate the mean distance L based on the simplified point set d , according to the point distance mean L d Calculate the eigenvalue sequence (λ max ,λ max-1 ...λ1,λ0), and determine the three axis vectors; after the point cloud data is simplified, it is necessary to build a directed bounding box model on this basis to help build a three-dimensional solid model. To solve this problem, the detailed implementation steps are as follows:

[0227] S51. Define the simplified point pj in the simplified point set t The coordinates are (A x ,B y ,C z ) and the number is D, calculate all the simplified points pj t The center of mass O zx The coordinates are calculated as follows:

[0228]

[0229] Among them, O zx represents the xth centroid coordinate;

[0230] S52, randomly select a simplified point pj t As a reference point pc t , calculate the reference point and other simplified points pj t The mean distance between the points L d , the calculation formula is:

[0231]

[0232] Among them, L d represents the mean value of the dth point distance;

[0233] S53, according to the point distance mean L d Calculate the weight coefficient w h , the calculation formula is:

[0234]

[0235] Among them, w h represents the hth weight coefficient;

[0236] S54, based on the weight coefficient w h Calculate the covariance matrix Xc used to determine the axial volume of the oriented bounding box p , the calculation formula is:

[0237]

[0238]

[0239] Among them, s represents the sum of weight coefficients, and m represents the simplified point pj t The weighted sum of

[0240] S55, according to the covariance matrix Xc p Calculate the characteristic equation and obtain multiple eigenvalues ​​λ E ,E=1,2...E, the calculation formula is:

[0241]

[0242] det(Xc p -λ E I)=0

[0243] Where det represents Xc p -λ E The determinant of I, where I represents the identity matrix, represents the eigenvector;

[0244] S56, multiple eigenvalues ​​λ E Arrange in ascending order to obtain the eigenvalue sequence (λ max ,λ max-1 ...λ1,λ0), the eigenvalue λ max ,λ max-1 and λ max-2 Set as the axis eigenvalue, and the eigenvalue λ max ,λ max-1 and λ max-2 The corresponding eigenvectors are set as the three axis vectors {u1, u2, u3} of the directed bounding box. Here, we need to use the directed bounding box method to construct an accurate 3D model of the cultural relic based on the simplified point set through the three axis vectors of the directed bounding box.

[0245] S6. Construct a directed bounding box model based on the three axis vectors {u1, u2, u3} and calculate the adaptive Gaussian kernel φ j (D), based on the adaptive Gaussian kernel φ j (D) Calculate global implicit values Based on step S5, it is necessary to improve the efficiency of model construction. To solve this problem, the specific implementation steps are as follows:

[0246] S61. Construct a directed bounding box model based on the axis vector and the simplified point set, and calculate the alignment point Dq used to align the simplified points. u , the calculation formula is:

[0247] Dq u =R d (pj t -O zx )

[0248]

[0249] Among them, R d Represents the transposed value of the d-th axis vector; the establishment of a directed bounding box model is a commonly used method in three-dimensional images and will not be described in detail here.

[0250] S62, along the axis vector {u1, u2, u3}, the directed bounding box space is defined as the network area, according to the point distance mean L d Calculate the partition side length Ψb , the calculation formula is:

[0251] Ψ b =1.5·median(L d )

[0252] Among them, b represents the length of the b-th partition side;

[0253] S63, based on the division of side length Ψ b Divide the network area into multiple sub-areas S c , calculate the approval value σ for each sub-region h , the calculation formula is:

[0254] σ h =median(L d inS c )

[0255] Among them, σ h represents the hth approved value;

[0256] S64, according to the approved value σ h Calculate the adaptive Gaussian kernel φ j (D), calculated as:

[0257]

[0258] Among them, φ j (D) represents the jth adaptive Gaussian kernel, Dq k and Dq k-1 denote the kth and k-1th alignment points respectively;

[0259] S65. Obtain the center point O of each sub-region by using the model coordinate method. D and half-length mean Bc z , calculate the global implicit value The calculation formula is:

[0260]

[0261] in, Represents the g-th normal vector. The global implicit value can effectively reduce the complexity of the algorithm operation. The axial alignment of the point cloud data can effectively reduce the error rate of the point cloud data. The use of the isosurface can reduce the processing time. Finally, the smoothing processing of the three-dimensional smoothing operator can better preserve the carving details of the cultural relics, make the three-dimensional model more accurate, improve the calculation efficiency and achieve higher accuracy.

[0262] S7, according to the global implicit value Calculate the isosurface eigenvalue η used to construct the 3D digital model of cultural relicsr , and smooth the three-dimensional digital model of the cultural relics to obtain an optimized three-dimensional model. Based on step S6, the model needs to be smoothed to retain more three-dimensional model details. In order to solve this problem, the detailed solution steps are as follows:

[0263] S71. According to the global implicit value Calculate the isosurface eigenvalue η r , the calculation formula is:

[0264]

[0265] Among them, η r represents the rth isosurface eigenvalue;

[0266] S72, according to the isosurface eigenvalue η r Construct a three-dimensional digital model of the cultural relic through the cube reconstruction method;

[0267] S73. Define the area of ​​each facet in the 3D digital model of the cultural relic as A. s , the normal vector is Calculating vertex normals The calculation formula is:

[0268]

[0269] Here, P represents the number of facets. Facets are the graphics that make up the outer surface of the three-dimensional digital model of the cultural relic. Through a method similar to calculus, multiple tiny graphics can be combined into a three-dimensional solid figure. The area and normal vector of the facets can be used to calculate the vertex normal vector of the image intersection position.

[0270] S74, according to the vertex normal vector Calculate the smoothing strength Qp h , the calculation formula is:

[0271]

[0272] Among them, u f represents the f-th axis vector;

[0273] S75. Define the vertex d of the three-dimensional digital model of the cultural relic in , the vertex of 1 ring area is N id , calculate the vertex offset Ld of the three-dimensional digital model of the cultural relic v , the calculation formula is:

[0274]

[0275] Where I represents the number of vertices in the 1-ring domain;

[0276] S76, according to the vertex offset Ld v Calculate the model smoothing value L ph , the calculation formula is:

[0277]

[0278] W(x) i =diag(Qp h )

[0279] Among them, Z represents the vertex offset Ld v the number of

[0280] S77, smoothing value L according to the model ph The three-dimensional digital model of the cultural relic is smoothed to obtain an optimized three-dimensional model. The global implicit value can effectively reduce the complexity of the algorithm operation, and the axial alignment of the point cloud data can effectively reduce the error rate of the point cloud data. The processing time is reduced by using the isosurface. Finally, the smoothing processing of the three-dimensional smoothing operator can better preserve the carving details of the cultural relic, making the three-dimensional model more accurate, which can not only improve the computing efficiency but also achieve higher accuracy.

[0281] Those skilled in the art will appreciate that all or part of the steps in the above-mentioned embodiment methods can be accomplished by instructing the relevant hardware through a program. Therefore, the present application may take the form of a complete hardware embodiment, a complete software embodiment, or an embodiment combining software and hardware. Furthermore, the present application may take the form of a computer program product implemented on one or more computer-usable storage media (including but not limited to disk storage, CD-ROM, optical storage, etc.) containing computer-usable program code.

[0282] Each embodiment in this specification is described in a progressive manner, with each embodiment focusing on the differences from other embodiments. Reference can be made to the same or similar parts between the embodiments. For the above embodiments, since they are basically similar to the method embodiments, the description is relatively simple. For relevant parts, refer to the partial description of the method embodiments.

[0283] The above embodiments provide a detailed introduction to the present invention. Specific examples are used herein to illustrate the principles and implementation methods of the present invention. The description of the above embodiments is only used to help understand the method of the present invention and its core ideas. At the same time, for those skilled in the art, according to the ideas of the present invention, there may be changes in the specific implementation methods and application scopes. In summary, the contents of this specification should not be understood as limiting the present invention.

Claims

1. A method for adaptive three-dimensional digital modeling of cultural relics, characterized in that: The method comprises the following steps: S1. Obtain the 3D point cloud data of the cultural relics to generate point cloud dataset B s , based on point cloud dataset B s Determine the neighborhood point set P k And obtain multiple scale features λ through feature decomposition n s ; S2, according to the scale feature λ n s Calculate the adaptive filtering threshold T for removing noise points y , and obtain the preprocessing point set; S3. Select any point in the preprocessing point set as the origin to construct a three-dimensional rectangular coordinate system. In the three-dimensional rectangular coordinate system, the surface fitting equation N for calculating the partial derivative is obtained by the surface fitting method. h ; S4. Calculate the normal vector based on the partial derivative Based on normal vector Calculate the first curvature threshold Y for simplifying the point cloud data q and the second curvature threshold Y w , and obtain a simplified point set; S5, calculate the eigenvalue sequence (λ based on the simplified point set max ,λ max-1 ...λ1,λ0), and through the eigenvalue sequence (λ max ,λ max-1 ...λ1,λ0) determine the three axis vectors {u1,u2,u3}; S6. Construct a directed bounding box model based on the three axis vectors {u1, u2, u3} and calculate the adaptive Gaussian kernel φ j (D), based on the adaptive Gaussian kernel φ j (D) Calculate global implicit values S7, according to the global implicit value Calculate the isosurface eigenvalue η used to construct the 3D digital model of cultural relics r , and the three-dimensional digital model of the cultural relics is smoothed to obtain an optimized three-dimensional model.

2. The adaptive modeling method according to claim 1, characterized in that: In step S1, the specific implementation steps are as follows: S11, the three-dimensional point cloud data of the cultural relics is based on the level radius r s Divide into cube grids with multiple levels, and build an octree structure based on the level radius. The expression is: r s =2 s ·d Among them, r s represents the radius of the sth level, where s=0,1,2...s; S12. Define point cloud dataset B s ={p i ,i=1,2...i},p i Represents the i-th data point, and selects a point p in the point cloud dataset a As the center point, get point p a The neighborhood point set P k ={p l ,l=1,2...l}; S13, according to the neighborhood point set P k Calculate the covariance matrix value X c , the calculation formula is: Where K represents p l the number of Represents the mean value of the neighborhood point coordinates; S14, the covariance matrix value X c Perform eigenvalue decomposition to obtain multiple scale features λ n s , where λ n s Represents the feature of the nth scale in the sth level, and the expression is: {λ0 s ,λ1 s ,...l n s }=eig(X c ) Among them, eig represents the covariance matrix value X c Perform eigendecomposition.

3. The adaptive modeling method according to claim 1, characterized in that: In step S2, the specific implementation steps are as follows: S21. Based on multiple scale features λ n s Calculate the multi-scale curvature characteristic ratio Tz b , the calculation formula is: Among them, λ k s Represented in multiple scale features λ n s The kth scale feature, N represents the scale feature λ n s the number of S22, according to the multi-scale curvature characteristic ratio Tz b Calculate the fusion scale eigenvalue Rh c , the calculation formula is: Among them, Rh c represents the cth fusion scale eigenvalue; S23, according to the neighborhood point set P k Calculate the first statistic μ s , the calculation formula is: Among them, μ s represents the sth first statistic; S24, according to the first statistic μ s Calculate the second statistic α t , the calculation formula is: Among them, α t represents the t-th second statistic; S25. According to the first statistic μ s and the second statistic α t Calculate the adaptive filtering threshold T y , the calculation formula is: T y =μ s +1.5·a t ·(1-Rh c ) Among them, T y represents the yth adaptive filtering threshold; S26, according to the adaptive filtering threshold T y Remove noise; like Then point p l Remove the noise and mark the removal point t T ; like Then point p l is a normal data point; S27, repeat step S26 until the point cloud data set is detected and a preliminary pre-processed point set is obtained; S28, according to the elimination point t T Calculate supplement point B bc , the calculation formula is: Among them, t T+1 and t T-1 Respectively represent and eliminate point t T Two adjacent points; S29, add point B bc Add to the preliminary preprocessing point set and obtain the preprocessing point set.

4. The adaptive modeling method according to claim 1, characterized in that: In step S3, the specific implementation steps are as follows: S31. Select a point p in the preprocessing point set d Establish a three-dimensional rectangular coordinate system for the origin and define any point p in the preprocessing point set i The coordinates of (x i ,y i ,z i ), and construct the surface fitting equation N of the preprocessing point set based on the three-dimensional rectangular coordinate system h , the expression is: N h =ax 2 +bxy+cy 2 Among them, a, b and c represent the surface fitting equation N h The coefficient of S32, according to the surface fitting equation N h Calculate the least squares value Q f , the calculation formula is: Where H represents the number of preprocessing points, z i Indicates the z-axis coordinate value of the preprocessing point; S33, based on the least squares value Q f Calculate the values ​​of coefficients a, b and c respectively, and the calculation formula is: Among them, Q f represents the fth least squares value; S34, fitting the surface equation N h Convert to surface function f e (x,y,z), the expression is: Among them, f e (x,y,z) represents the e-th surface function; S35, according to the surface function f e (x, y, z) calculates the partial derivative, which includes the first partial derivative f x1 , the second partial derivative f y2 , the third partial derivative f xy3 、The fourth partial derivative f xx4 and the fifth partial derivative f yy5 , the calculation formula is: in, Represents the surface function f e Partial derivatives of (x,y,z).

5. The adaptive modeling method according to claim 1, characterized in that: In step S4, the specific implementation steps are as follows: S41, according to the first partial derivative f in the partial derivative x1 and the second partial derivative f y2 Calculating the normal vector The calculation formula is: in, Represents the surface function f e The normal vector of (x,y,z); S42, according to the first partial derivative f x1 and the second partial derivative f y2 Calculate the curvature value E g 、F h and G j , the calculation formula is: E g =f x1 ·f x1 F h =f x1 ·f y2 G j =f y2 ·f y2 Among them, E g 、F h and G j Represent the g, h and j-th curvature values ​​respectively; S43, according to the third partial derivative f xy3 、The fourth partial derivative f xx4 and the fifth partial derivative f yy5 Calculate the depth value L k 、M l and N m , the calculation formula is: Among them, L k 、M l and N m Represent the kth, lth and mth depth values ​​respectively; S44. Calculate the Gaussian curvature K used to screen surface features n , mean curvature H pj and global curvature S45. Repeat steps S41-S44 to traverse all points in the preprocessing point set; S46. Calculate the first curvature threshold Y q and the second curvature threshold Y w , the calculation formula is: Among them, Y q represents the qth curvature threshold; S47, according to the curvature threshold Y q Complete the screening of the preprocessing point set; If K n ≥Y w and It is marked as a reserved point; If K n <Y w and These are points that can be streamlined and eliminated; If K n <Y w or It is marked as a reserved point; S48. Repeat step S47 to simplify all pre-processed point sets to obtain a simplified point set and end.

6. The adaptive modeling method according to claim 5, characterized in that: The Gaussian curvature K n , mean curvature H pj and global curvature The calculation formula is: Where H represents the mean curvature H pj The number of 7. The adaptive modeling method according to claim 1, characterized in that: In step S5, the specific implementation steps are as follows: S51. Define the simplified point pj in the simplified point set t The coordinates are (A x ,B y ,C z ) and the number is D, calculate all the simplified points pj t The center of mass O zx The coordinates are calculated as follows: Among them, O zx represents the xth centroid coordinate; S52, randomly select a simplified point pj t As a reference point pc t , calculate the reference point and other simplified points pj t The mean distance between the points L d , the calculation formula is: Among them, L d represents the mean value of the dth point distance; S53, according to the point distance mean L d Calculate the weight coefficient w h , the calculation formula is: Among them, w h represents the hth weight coefficient; S54, based on the weight coefficient w h Calculate the covariance matrix Xc used to determine the axial volume of the oriented bounding box p , the calculation formula is: Among them, s represents the sum of weight coefficients, and m represents the simplified point pj t The weighted sum of S55, according to the covariance matrix Xc p Calculate the characteristic equation and obtain multiple eigenvalues ​​λ E ,E=1,2...E, the calculation formula is: it(Xc p -λ E I)=0 Where det represents Xc p -λ E The determinant of I, where I represents the identity matrix, represents the eigenvector; S56, multiple eigenvalues ​​λ E Arrange in ascending order to obtain the eigenvalue sequence (λ max ,λ max-1 ...λ1,λ0), the eigenvalue λ max ,λ max-1 and λ max-2 Set as the axis eigenvalue, and the eigenvalue λ max ,λ max-1 and λ max-2 The corresponding eigenvectors are set as the three axis vectors {u1,u2,u3} of the directed bounding box.

8. The adaptive modeling method according to claim 1, characterized in that: In step S6, the specific implementation steps are as follows: S61. Construct a directed bounding box model based on the axis vector and the simplified point set, and calculate the alignment point Dq used to align the simplified points. u , the calculation formula is: Dq u =R d (p.j. t -ISLAND zx ) R d =[u1 T ,u2 T ,u3 T ] Among them, R d Represents the transposed value of the d-th axis vector; S62, along the axis vector {u1, u2, u3}, the directed bounding box space is defined as the network area, according to the point distance mean L d Calculate the partition side length Ψ b , the calculation formula is: P b =1.5·median(L d ) Among them, b represents the length of the b-th partition side; S63, based on the division of side length Ψ b Divide the network area into multiple sub-areas S c , calculate the approval value σ for each sub-region h , the calculation formula is: s h =median(L d inS c ) Among them, σ h represents the hth approved value; S64, according to the approved value σ h Calculate the adaptive Gaussian kernel φ j (D), calculated as: Among them, φ j (D) represents the jth adaptive Gaussian kernel, Dq k and Dq k-1 denote the kth and k-1th alignment points respectively; S65. Obtain the center point O of each sub-region by using the model coordinate method. D and half-length mean Bc z , calculate the global implicit value The calculation formula is: in, Represents the g-th normal vector.

9. The adaptive modeling method according to claim 1, characterized in that: In step S7, the specific implementation steps are as follows: S71. According to the global implicit value Calculate the isosurface eigenvalue η r , the calculation formula is: Among them, η r represents the rth isosurface eigenvalue; S72, according to the isosurface eigenvalue η r Construct a three-dimensional digital model of the cultural relic through the cube reconstruction method; S73. Define the area of ​​each facet in the 3D digital model of the cultural relic as A. s , the normal vector is Calculating vertex normals The calculation formula is: Where P represents the number of patches; S74, according to the vertex normal vector Calculate the smoothing strength Qp h , the calculation formula is: Among them, u f represents the f-th axis vector; S75. Define the vertex d of the three-dimensional digital model of the cultural relic in , the vertex of 1 ring area is N id , calculate the vertex offset Ld of the three-dimensional digital model of the cultural relic v , the calculation formula is: Where I represents the number of vertices in the 1-ring domain; S76, according to the vertex offset Ld v Calculate the model smoothing value L ph , the calculation formula is: W(x) i =diag(Qp h ) Among them, Z represents the vertex offset Ld v the number of S77, smoothing value L according to the model ph The surface of the three-dimensional digital model of the cultural relics is smoothed to obtain an optimized three-dimensional model.

Images