A three-dimensional digitalization adaptive modeling method for cultural relics

Abstract

The present application 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 three-dimensional digital adaptive modeling method for cultural relics, 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 multiple scale features through feature decomposition; S2, calculating an adaptive filtering threshold for removing noise points according to the scale features, and obtaining a preprocessed point set; the present application can effectively reduce the complexity of algorithm operation through global implicit values, effectively reduce the error rate of point cloud data through axial alignment of point cloud data, reduce the processing time through isosurface, and finally perform smoothing processing through a three-dimensional smoothing operator, so that the carving details of the cultural relic can be better preserved, the three-dimensional model is more accurate, the operation efficiency is improved, and high accuracy is achieved.

Description

Technical Field

[0001] This invention relates to the field of cultural relic modeling technology, and in particular to a three-dimensional digital adaptive modeling method for cultural relics. Background Technology

[0002] Cultural relic modeling generally refers to the process of digitally reconstructing cultural relics in three dimensions using computer technology and related software tools. By collecting information such as the geometric shape, surface texture, and color of the cultural relic, a three-dimensional digital model is constructed to realize the digital protection, display, research, and dissemination of the cultural relic. Existing three-dimensional modeling methods for cultural relics generally adopt image measurement modeling or scanning modeling methods. These methods involve acquiring and extracting features from multiple images, and then using two-dimensional to three-dimensional algorithms to convert planar images into three-dimensional models. Alternatively, they can directly obtain point cloud data of the cultural relic and extract three-dimensional features. However, these methods have several drawbacks, including large data requirements, huge computational load, low efficiency, and low model processing accuracy, making it difficult to accurately construct cultural relic models and resulting in an unsatisfactory user experience. Summary of the Invention

[0003] To address the shortcomings of existing technologies, this invention provides a three-dimensional digital adaptive modeling method for cultural relics, which solves the technical problems of low computational efficiency and low modeling accuracy in existing technologies, thereby achieving the goal of significantly improving computational efficiency and modeling accuracy.

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

[0005] S1. Obtain the 3D point cloud data of the cultural relic and generate a point cloud dataset B. s Based on point cloud dataset B s Determine the neighborhood point set P k And multiple scale features λ are obtained through eigenvalue decomposition. n s ;

[0006] S2, Based on the scale characteristic λ n s Calculate the adaptive filtering threshold T used to remove noise points. y And obtain the preprocessed 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, obtain the surface fitting equation N for calculating the partial derivatives using the surface fitting method. h ;

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

[0009] S5, calculating a feature value sequence (λ max , λ max-1 … λ1, λ0) based on the simplified point set, and determining three axial vectors {u1, u2, u3} through the feature value sequence (λ max , λ max-1 … λ1, λ0);

[0010] S6, constructing a directional bounding box model according to the three axial vectors {u1, u2, u3} and calculating an adaptive Gaussian kernel φ j (D), calculating a global implicit value j based on the adaptive Gaussian kernel φ r (D);

[0011] S7, calculating an isosurface feature value η for constructing a three-dimensional digital model of cultural relics based on the global implicit value s , and performing smoothing processing on the three-dimensional digital model of cultural relics to obtain an optimized three-dimensional model.

[0012] Further, in step S1, the following steps are specifically implemented:

[0013] S11, dividing the three-dimensional point cloud data of cultural relics into a cubic grid with multiple levels according to a level radius r s , and building an octree structure according to the level radius, the expression being:

[0014] r s = 2 s · δ

[0015] wherein r s represents the s-th level radius, and s = 0, 1, 2... s;

[0016] S12, defining a point cloud data set B s = {p i , i = 1, 2... i}, p i representing the i-th data point, and selecting any point p a in the point cloud data set as a center point, obtaining a neighborhood point set P a = {p k , l = 1, 2... l} of the point p l ;

[0017] S13, calculating a covariance matrix value X k based on the neighborhood point set P c , the calculation formula being:

[0018]

[0019] wherein K represents p l the number of points, denotes the mean value of the coordinates of the neighborhood points;

[0020] S14, performing eigenvalue decomposition on the covariance matrix value X c to obtain a plurality of scale features λ n s wherein λ n s denotes the feature of the n-th scale in the s-th hierarchy, and the expression is:

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

[0022] wherein eig denotes performing eigenvalue decomposition on the covariance matrix value X c

[0023] Further, in step S2, the following steps are specifically implemented:

[0024] S21, calculating a multi-scale curvature feature ratio Tz n s according to the plurality of scale features λ b , and the calculation formula is:

[0025]

[0026] wherein λ k s denotes the k-th scale feature in the plurality of scale features λ n s ; and N denotes the number of scale features λ n s .

[0027] S22, calculating a fused scale feature value Rh b according to the multi-scale curvature feature ratio Tz c , and the calculation formula is:

[0028]

[0029] wherein Rh c denotes the c-th fused scale feature value;

[0030] S23, calculating a first statistical quantity μ k according to the neighborhood point set P s , and the calculation formula is:​

[0031]

[0032] Where, μ s This represents the s-th first statistic;

[0033] S24. Based on the first statistic μ s Calculate the second statistic α t The calculation formula is:

[0034]

[0035] Where, α t This represents the t-th second statistic;

[0036] S25. Based on the first statistic μ s Second statistic α t Calculate the adaptive filter threshold T y The calculation formula is:

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

[0038] Among them, T y This represents the y-th adaptive filtering threshold;

[0039] S26. Based on the adaptive filtering threshold T y Noise removal;

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

[0041] like Then point p l These are normal data points;

[0042] S27. Repeat step S26 until the point cloud dataset is detected and a preliminary preprocessed point set is obtained;

[0043] S28. Based on the removal point t T Calculate supplementary point B bc The calculation formula is:

[0044]

[0045] Among them, t T+1 and t T-1 Representing the rejection point t respectively T Two adjacent points;

[0046] S29, Add point B bc Add it to the initial preprocessing point set to obtain the preprocessing point set.

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

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

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

[0050] Where a, b, and c represent the surface fitting equation N, respectively. h The coefficient;

[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 This represents the z-axis coordinate value of the preprocessing point;

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

[0055]

[0056] Among them, Q f This represents the f-th least squares value;

[0057] S34. Fit 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. Based on the surface function fe (x,y,z) compute partial derivatives, including first partial derivative f x1 , second partial derivative f y2 , third partial derivative f xy3 , fourth partial derivative f xx4 and fifth partial derivative f yy5 , the formula is:

[0061]

[0062] wherein, denotes the partial derivative of the surface function f e (x,y,z).

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

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

[0065]

[0066] wherein, denotes the normal vector of the surface function f e (x,y,z);

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

[0068] E g = f x1 · f x1

[0069] F h = f x1 · f y2

[0070] G j = f y2 · f y2

[0071] wherein, E g , F h and G j respectively denote the g, h and j-th curvature values;

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

[0073]

[0074] wherein, L k , M l and N m respectively represent the kth, lth and mth sag values;

[0075] S44, the Gaussian curvature K for screening the surface feature is calculated n , the average curvature H pj and the global curvature

[0076] S45, the steps S41-S44 are repeated to traverse all the points in the preprocessed point set;

[0077] S46, the first curvature threshold Y q and the second curvature threshold Y w are calculated, and the calculation formula is as follows:

[0078]

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

[0080] S47, the screening of the preprocessed point set is completed according to the curvature threshold Y q ;

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

[0082] If K n < Y w and , the point is a reducible point and is eliminated;

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

[0084] S48, the step S47 is repeated to complete the reduction of all the preprocessed point set to obtain a reduced point set and end.

[0085] Further, the calculation formula of the Gaussian curvature K n , the average curvature H pj and the global curvature is as follows:

[0086]

[0087] Where H represents the mean curvature H pj The quantity.

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

[0089] S51. Define the simplified point pj in the simplified point set. t Coordinates are (A) x B y C z Given a quantity of D, calculate all simplified points pj. t The center of mass O zx The coordinates are calculated using the following formula:

[0090]

[0091] Among them, O zx Indicates the coordinates of the x-th centroid;

[0092] S52. Randomly select a simplified point pj t PC as a reference point t Calculate the reference point and other simplified points pj t Mean distance L between points d The calculation formula is:

[0093]

[0094] Among them, L d This represents the distance from the d-th point to the mean.

[0095] S53, Based on the average point distance L d Calculate the weighting coefficient w h The calculation formula is:

[0096]

[0097] Among them, w h This represents the h-th weight coefficient;

[0098] S54, Based on weighting coefficient w h Calculate the covariance matrix Xc used to determine the axial quantities of the directed bounding box. p The calculation formula is:

[0099]

[0100] Where s represents the sum of weight coefficients, and m represents the simplified point pj t Weighted summation;

[0101] S55, according to the covariance matrix Xc p Calculate the characteristic equation and get a plurality of eigenvalues λ E , E = 1, 2... E, the formula is:

[0102]

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

[0104] Wherein, det represents the determinant of Xc p -λ E I, I represents the unit matrix, represents the eigenvector;

[0105] S56, arrange a plurality of eigenvalues λ E in ascending order to get the eigenvalue sequence (λ max , λ max-1 ... λ1, λ0), set the eigenvalues λ max , λ max-1 and λ max-2 as the axis eigenvalues, and set the eigenvectors corresponding to the eigenvalues λ max , λ max-1 and λ max-2 as the three axis vectors {u1, u2, u3} of the oriented bounding box.

[0106] Further, in step S6, the following steps are implemented:

[0107] S61, construct the oriented bounding box model according to the axis vector and the simplified point set, and calculate the alignment point Dq u for aligning the simplified points, the formula is:

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

[0109]

[0110] Wherein, R d represents the transpose value of the dth axis vector;

[0111] S62, define the oriented bounding box space as a network area along the axis vector {u1, u2, u3}, and calculate the division edge length Ψ d according to the point distance mean L b , the formula is:

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

[0113] wherein Ψ b denotes the b-th partition edge length;

[0114] S63, according to the partition edge length Ψ b the network area is divided into a plurality of sub-areas S c , the approval value σ h of each sub-area is calculated

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

[0116] wherein σ h denotes the h-th approval value;

[0117] S64, according to the approval value σ h the adaptive Gaussian kernel φ j (D) is calculated, and the calculation formula is:

[0118]

[0119] wherein φ j (D) denotes the j-th adaptive Gaussian kernel, Dq k and Dq k-1 denote the k-th and k-1-th alignment points, respectively;

[0120] S65, the center point O D and the half-length mean value Bc z of each sub-area are obtained by the model coordinate method, and the global implicit value is calculated The calculation formula is:

[0121]

[0122] wherein, denotes the g-th normal vector.

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

[0124] S71, according to the global implicit value the isosurface feature value η r is calculated, and the calculation formula is:

[0125]

[0126] wherein η r denotes the r-th isosurface feature value;

[0127] S72, according to the isosurface feature value η rA three-dimensional digital model of the cultural relic is constructed by a cubic reconstruction method;

[0128] S73, define an area of each face sheet in the three-dimensional digital model of the cultural relic as A s , and a normal vector as Calculate the vertex normal vector The calculation formula is:

[0129]

[0130] Wherein, P represents the number of face sheets;

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

[0132]

[0133] Wherein, u f represents the f-th axial vector;

[0134] S75, define the vertex d in of the three-dimensional digital model of the cultural relic, and the 1-ring domain vertex as N id Calculate the vertex offset Ld v of the three-dimensional digital model of the cultural relic, and the calculation formula is:

[0135]

[0136] Wherein, I represents the number of 1-ring domain vertices;

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

[0138]

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

[0140] Wherein, Z represents the number of vertex offsets Ld v ;

[0141] S77, perform surface smoothing processing on the three-dimensional digital model of the cultural relic according to the model smoothing value L ph , and obtain an optimized three-dimensional model.

[0142] Through the above technical scheme, the present application provides a three-dimensional digital adaptive modeling method for cultural relics, which at least has the following beneficial effects:

[0143] 1. The application can make the point cloud data more accurate through multi-scale feature processing and adaptive filtering processing of three-dimensional point cloud data, can accelerate three-dimensional data search through octree structure, can reduce the case of feature point misdeletion through curvature threshold method while quantifying the surface bending degree, and the use of adaptive threshold can improve the robustness of the algorithm.

[0144] 2. The application can extract features of the surface based on the fitted surface, determine the reducibility of the point cloud data according to the characteristic value, complete the simplification of the point cloud data, and construct a directed bounding box model to prepare for the subsequent three-dimensional model construction, so that the simplification step of the point cloud data is fast and efficient, the simplified model can reduce the operation steps and the number of data points, and the algorithm operation efficiency is greatly improved while ensuring the accuracy of the vertex.

[0145] 3. The application can effectively reduce the complexity of algorithm operation through global implicit value, effectively reduce the error rate of point cloud data through axial alignment of point cloud data, reduce processing time through isosurface, and finally perform smoothing processing through a three-dimensional smoothing operator, so that the carving details of cultural relics can be better preserved, the three-dimensional model is more accurate, the operation efficiency is improved, and high accuracy is achieved. BRIEF DESCRIPTION OF DRAWINGS

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

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

[0148] In order to make the above-mentioned purposes, features and advantages of the present application more obvious, easy to understand, the present application will be further described in detail below in combination with the drawings and specific embodiments. The realization process of how to apply technical means to solve technical problems and achieve technical effects of the present application can be fully understood and implemented.

[0149] In order to solve the technical problems of low operation efficiency and low modeling accuracy of the prior art, the embodiment proposes a three-dimensional digital adaptive modeling method for cultural relics, which can greatly improve the operation efficiency and improve the modeling accuracy, as shown in Figure 1 The method comprises the following steps:

[0150] S1, obtaining three-dimensional point cloud data of cultural relics to generate point cloud data set B s , determining the neighborhood point set P s based on the point cloud data set B kAnd through the characteristic decomposition of multiple scale features λ n s In the process of acquiring point cloud data, due to the accuracy of the device and the influence of the dim environment light, it is easy to cause the point cloud data to be inaccurate, in order to solve this problem, the specific implementation steps are as follows:

[0151] S11, the three-dimensional point cloud data of cultural relics is divided into a plurality of levels according to the level radius r s The octree structure is built according to the level radius, and the expression is:

[0152] r s =2 s ·δ

[0153] Wherein, r s Indicates the s-th level radius, and s=0,1,2...s; The octree structure is a commonly used three-dimensional data structure, which can improve the processing efficiency of three-dimensional data preprocessing, and has good accuracy, which is not described here.

[0154] S12, define the point cloud data set B s ={p i ,i=1,2...i}, p i Indicates the i-th data point, and any point p a In the point cloud data set is selected as the center point, and the neighborhood point set P a ={p k ,l=1,2...l} of point p l Is obtained; p l Indicates the data point in the neighborhood point set, and the neighborhood point set is the data point within the level radius r s Of the center point.

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

[0156]

[0157] Wherein, K represents the number of p l , Indicates the mean value of the neighborhood point coordinates;

[0158] S14, the covariance matrix value X c Carry out eigenvalue decomposition to obtain a plurality of scale features λ n s , wherein λ n s Indicates the s-th level of the n-th scale feature, and the expression is:

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

[0160] wherein eig represents eigen decomposition of the covariance matrix value X c , through multi-scale feature processing and adaptive filtering processing of three-dimensional point cloud data by algorithm, the point cloud data can be more accurate, through octree structure can accelerate three-dimensional data search, quantize the surface bending degree, at the same time, through the curvature threshold method, reduce the feature point misdeletion situation, and 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 y for removing noise points, and obtain the pretreatment point set; on the basis of step S1, the point cloud data is adaptively filtered, in order to solve this problem, the detailed steps are:

[0162] S21, according to a plurality of scale features λ n s Calculate the multi-scale curvature feature ratio Tz b , the calculation formula is:

[0163]

[0164] wherein λ k s represents a plurality of scale features λ n s The kth scale feature, N represents the number of scale features λ n s ;

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

[0166]

[0167] wherein Rh c represents the cth fusion scale feature value;

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

[0169]

[0170] Where, μ s Let K represent the s-th first statistic, and K represent the point p in the neighborhood set. l Quantity;

[0171] S24. Based on the first statistic μ s Calculate the second statistic α t The calculation formula is:

[0172]

[0173] Where, α t This represents the t-th second statistic;

[0174] S25. Based on the first statistic μ s Second statistic α t Calculate the adaptive filter threshold T y The calculation formula is:

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

[0176] Among them, T y This represents the y-th adaptive filtering threshold;

[0177] S26. Based on the adaptive filtering threshold T y Noise removal;

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

[0179] like Then point p l These are normal data points;

[0180] S27. Repeat step S26 until the point cloud dataset is detected and a preliminary preprocessed point set is obtained;

[0181] S28. Based on the removal point t T Calculate supplementary point B bc The calculation formula is:

[0182]

[0183] Among them, t T+1 and t T-1 Representing the rejection point t respectively T Two adjacent points;

[0184] S29, Add point B bcThe point cloud data is more accurate through multi-scale feature processing and adaptive filtering processing of the three-dimensional point cloud data, the three-dimensional data search is accelerated through the octree structure, the curvature threshold method is used to reduce the feature point deletion error, and the adaptive threshold is used to improve the robustness of the algorithm.

[0185] S3, a three-dimensional rectangular coordinate system is constructed based on any point in the pre-processed point set, and a surface fitting equation N h is obtained according to the point in the three-dimensional rectangular coordinate system h The partial derivative value is calculated; since the point cloud data in the pre-processed point set is large, the step-by-step operation is very complex, which seriously affects the operation efficiency, in order to solve this problem, the detailed implementation steps are as follows:

[0186] S31, a point p d is selected in the pre-processed point set as the origin, a three-dimensional rectangular coordinate system is established, wherein the coordinates of any point p i in the pre-processed point set are (x i , y i , z i ), and a surface fitting equation N h of the pre-processed point set is constructed based on the three-dimensional rectangular coordinate system, and the expression is:

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

[0188] Wherein, a, b and c represent the coefficients of the surface fitting equation N h ;

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

[0190]

[0191] Wherein, H represents the number of pre-processed points, and z i represents the z-axis coordinate value of the pre-processed point;

[0192] S33, the values of the coefficients a, b and c are calculated based on the least square value Q f , and the calculation formula is:

[0193]

[0194] Wherein, Q f represents the fth least square value;

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

[0196]

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

[0198] S35, according to the curved surface function f e (x, y, z) to calculate the partial derivative value, the partial derivative value includes the first partial derivative value f x1 , the second partial derivative value f y2 , the third partial derivative value f xy3 , the fourth partial derivative value f xx4 And the fifth partial derivative value f yy5 , the calculation formula is:

[0199]

[0200]

[0201] Where, Indicates the partial derivative of the curved surface function f e (x, y, z), through the construction of the curved surface fitting equation, the feature of the curved surface can be extracted on the basis of the fitting curved surface, and the reducibility of the point cloud data is determined according to the characteristic value, the point cloud data is simplified, and the oriented bounding box model is constructed, which prepares for the subsequent three-dimensional model construction. The simplification step of point cloud data is fast and efficient, the simplified model can reduce the operation steps and the number of data points, while ensuring the accuracy of the vertex, and greatly improving the operation efficiency of the algorithm.

[0202] S4, calculating the normal vector according to the partial derivative value Based on the normal vector Calculate the first curvature threshold Y q And the second curvature threshold Y w For simplifying the point cloud data, and get the simplified point set; On the basis of step S3, the characteristics of the three-dimensional curved surface need to be further obtained by the partial derivative value, and the point cloud data is simplified based on the characteristics. The specific solution steps are as follows:

[0203] S41, calculating the normal vector according to the first partial derivative value f x1 And the second partial derivative value f y2 In the partial derivative value The calculation formula is:

[0204]

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

[0206] S42. Based on the first partial derivative value f x1 Second partial derivative value 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 These represent the g-th, h-th, and j-th curvature values, respectively.

[0211] S43. Based on the third partial derivative value f xy3 The fourth partial derivative value f xx4 and the fifth partial derivative value 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 These represent the k-th, l-th, and m-th curvature depth values, respectively.

[0214] S44. Calculate the Gaussian curvature K used for screening surface features. n Mean curvature H pj and global curvature The calculation formula is:

[0215]

[0216] Where H represents the mean curvature H pj Quantity;

[0217] S45, repeat steps S41-S44, traverse all points in the pretreatment point set;

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

[0219]

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

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

[0222] If K n ≥ Y w and Then mark as a reserved point;

[0223] If K n < Y w and Then it is a reducible point and is rejected;

[0224] If K n < Y w or Then mark as a reserved point; the reserved point generally indicates that the surface feature is very obvious, which belongs to the point on the boundary or cambered surface. The reducible point is generally the point on the plane or relatively flat circular arc surface. These points do not need to be calculated completely, but only a small amount is needed to complete, so they can be rejected.

[0225] S48, repeat step S47 to complete the reduction of all pretreatment point set to obtain the reduced point set and end. Through the construction of surface fitting equation, the feature of the surface can be extracted on the basis of the fitting surface, and the reducibility of the point cloud data is determined according to the characteristic value, the reduction of the point cloud data is completed, and the oriented bounding box model is constructed, which prepares for the subsequent three-dimensional model construction. The reduction step of point cloud data is fast and efficient, and the reduced model can reduce the operation steps and the number of data points. While ensuring the accuracy of the vertex, the algorithm operation efficiency is greatly improved.

[0226] S5, calculate the point distance mean L d based on the reduced point set, calculate the characteristic value sequence (λ d , λ max ... λ1, λ0) according to the point distance mean L max-1 , and determine three axial vectors; after the reduction of the point cloud data, a directed bounding box model needs to be constructed on this basis, which is helpful to establish a three-dimensional model. In order to solve this problem, the detailed implementation steps are as follows:

[0227] S51, define a simplified point pj in the simplified point set t The coordinates of the centroid O of all simplified points pj are calculated, and the calculation formula is: x y z t zx

[0228]

[0229] where O represents the xth centroid coordinate; zx

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

[0231]

[0232] where L represents the dth average distance; d

[0233] S53, calculate the weight coefficient w according to the average distance L d , and the calculation formula is: h

[0234]

[0235] where w represents the hth weight coefficient; h

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

[0237]

[0238]

[0239] where s represents the sum of the weight coefficients, and m represents the weighted sum of the simplified points pj t

[0240] S55, calculate the characteristic equation according to the covariance matrix Xc p and obtain a plurality of characteristic values λ E , E = 1, 2... E, and 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 The eigenvalue sequence (λ) is obtained by arranging the values ​​in ascending order. max ,λ max-1 ...λ1,λ0), the eigenvalues ​​λ max , λ max-1 and λ max-2 Set as the axis eigenvalue, and with eigenvalue λ max , λ max-1 and λ max-2 The corresponding feature vectors are set as the three axis vectors {u1, u2, u3} of a 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 using 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 adaptive Gaussian kernel φ j (D) Calculate global implicit values Building upon step S5, it is necessary to improve the efficiency of model construction. To address this issue, the specific implementation steps are as follows:

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

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

[0248]

[0249] Among them, R d This represents the transpose of the d-th axis vector; the construction of the directed bounding box model is a commonly used method in 3D imaging, which will not be elaborated here.

[0250] S62. Define the directed bounding box space as the network region along the axis vector {u1,u2,u3}, based on the mean point distance L. d Calculate the side length Ψb , the calculation formula is:

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

[0252] Wherein, Ψ b represents the b-th partition edge length;

[0253] S63, according to the partition edge length Ψ b The network area is divided into a plurality of sub-regions S c , calculate the approval value σ h of each sub-region, the calculation formula is:

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

[0255] Wherein, σ h represents the h-th approval value;

[0256] S64, according to the approval value σ h Adaptive Gaussian kernel φ j (D) is calculated, the calculation formula is:

[0257]

[0258] Wherein, φ j (D) represents the j-th adaptive Gaussian kernel, Dq k and Dq k-1 respectively represent the k and k-1 alignment points;

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

[0260]

[0261] Wherein, represents the g-th normal vector, the global implicit value can effectively reduce the complexity of algorithm operation, the axial alignment point cloud data can effectively reduce the error rate of point cloud data, the processing time is reduced by using the isosurface, and finally the smoothing processing by the three-dimensional smoothing operator can better preserve the sculpture details of cultural relics, make the three-dimensional model more accurate, improve the operation efficiency and achieve higher accuracy.

[0262] S7, according to the global implicit value The isosurface characteristic value η used for constructing the three-dimensional digital model of cultural relics is calculatedr And the cultural relics three-dimensional digital model is smoothed to obtain an optimized three-dimensional model. On the basis of step S6, the model needs to be smoothed to retain more three-dimensional model details. To solve this problem, the detailed solving steps are as follows:

[0263] S71, according to the global implicit value Calculate the characteristic value of the isosurface η r The calculation formula is:

[0264]

[0265] Wherein, η r represents the rth characteristic value of the isosurface;

[0266] S72, according to the characteristic value of the isosurface η r The three-dimensional digital model of cultural relics is constructed by cubic reconstruction method;

[0267] S73, define the area of each face sheet in the three-dimensional digital model of cultural relics as A s , and the normal vector is Calculate the vertex normal vector The calculation formula is:

[0268]

[0269] Wherein, P represents the number of face sheets; The face sheet is a graphic that constitutes the outer surface of the three-dimensional digital model of cultural relics. A three-dimensional graphic can be composed of multiple small graphics by a method similar to calculus. The vertex normal vector at the intersection of the image can be calculated by the area and normal vector of the face sheet.

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

[0271]

[0272] Wherein, u f represents the fth axial vector;

[0273] S75, define the vertex d in of the three-dimensional digital model of cultural relics, the 1-ring domain vertex is N id , calculate the vertex offset Ld v of the three-dimensional digital model of cultural relics, and the calculation formula is:

[0274]

[0275] Wherein, I represents the number of 1-ring domain vertices;

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

[0277]

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

[0279] Wherein, Z represents the number of vertex offset Ld v ;

[0280] S77, according to the model smoothing value L ph The surface smoothing processing is carried out on the three-dimensional digital model of cultural relics, and the optimized three-dimensional model is obtained, the complexity of algorithm operation can be effectively reduced through the global implicit value, the error rate of point cloud data can be effectively reduced through the axial alignment point cloud data, the processing time is reduced through the isosurface, and finally through the smoothing processing of three-dimensional smoothing operator, the sculpture details of cultural relics can be better preserved, the three-dimensional model is more accurate, the operation efficiency can be improved, and high accuracy can be achieved.

[0281] Those skilled in the art can understand that all or part of the steps in the above-mentioned embodiment methods can be completed by programs instructing related hardware, therefore, the present application can adopt a completely hardware embodiment, a completely software embodiment, or an embodiment combining software and hardware aspects. Moreover, the present application can adopt the form of a computer program product implemented on one or more computer usable storage media (including but not limited to disk memory, CD-ROM, optical memory, etc.) containing computer usable program codes.

[0282] Each embodiment in the specification is described in a progressive manner, and each embodiment focuses on the difference from other embodiments. The same or similar parts between each embodiment can be referred to each other. For each of the above embodiments, since it is basically similar to the method embodiment, it is described relatively simply, and the related part can be referred to the part of the method embodiment.

[0283] The above embodiments are described in detail, and the principles and implementation modes of the present application are described by applying specific examples. The above embodiment is only used to help understand the method and core idea of the present application; at the same time, for those skilled in the art, according to the idea of the present application, the specific implementation mode and application range will be changed, and the above description should not be understood as the limitation of the present application.

Claims

1. A method for adaptive three-dimensional digital modeling of cultural relics, characterized in that, The method includes the following steps: S1. Obtain the 3D point cloud data of the cultural relic and generate a point cloud dataset B. s Based on point cloud dataset B s Determine the neighborhood point set P k And multiple scale features λ are obtained through eigenvalue decomposition. n s ; S2, Based on the scale characteristic λ n s Calculate the adaptive filtering threshold T used to remove noise points. y And obtain the preprocessed 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, obtain the surface fitting equation N for calculating the partial derivatives using the surface fitting method. h ; S4. Calculate the normal vector based on the partial derivative values. Based on normal vector Calculate the first curvature threshold Y used to simplify the point cloud data. q Second curvature threshold Y w And obtain a simplified point set, including: Calculate the Gaussian curvature K used for screening surface features. n Mean curvature H pj and global curvature Calculate the first curvature threshold Y q Second curvature threshold Y w The calculation formula is: Among them, Y q H represents the q-th first curvature threshold; H represents the number of preprocessing points. According to the first curvature threshold Y q Complete the screening of the preprocessed point set; If K n ≥Y w And H pj ≥Y q If so, then mark it as a reserved point; If K n <Y w And H pj <Y q Then, the points can be simplified and eliminated; If K n <Y w or H pj <Y q If so, then mark it as a reserved point; S5. Calculate the eigenvalue sequence (λ) based on a 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 adaptive Gaussian kernel φ j (D) Calculate the global implicit value ▽Yn v ; S7, based on the global implicit value ▽Yn v Calculate the isosurface eigenvalues ​​η used to construct a three-dimensional digital model of the cultural relic. r Furthermore, the 3D digital model of the cultural relic is smoothed to obtain an optimized 3D 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 processed according to the layer radius r. s Divide the data into a multi-level cubic mesh and construct an octree structure based on the level radius, as expressed in the following expression: r s =2 s ·d Where, r s Let represent the radius of the s-th level, where s = 0, 1, 2...s; S12, Define point cloud dataset B s ={p i p, i = 1, 2, ..., i i Let p represent the i-th data point, where any point p is randomly selected from the point cloud dataset. a As the center point, obtain point p. a neighborhood point set P k ={p l ,l=1,2...l}; S13, Based on the neighborhood point set P k Calculate the covariance matrix value X c The calculation formula is: Where K represents the neighborhood point p l Quantity, This represents the average coordinates of the neighboring points; S14, change the covariance matrix value X c Eigenvalue decomposition yields multiple scale features λ n s , where λ n s The feature at the nth scale in the s-th level is expressed as: {λ0 s ,λ1 s ,...l n s }=eig(X c ) Where eig represents the expression for the covariance matrix value X c Perform feature decomposition.

3. The adaptive modeling method according to claim 2, 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: Where, λ k s Indicates the feature λ at multiple scales n s The k-th scale feature, N represents the scale feature λ. n s Quantity; S22, Based on the multi-scale curvature characteristic ratio Tz b Calculate the fusion scale eigenvalue Rh c The calculation formula is: Among them, Rh c This represents the c-th fusion scale feature value; S23, Based on the neighborhood point set P k Calculate the first statistic μ s The calculation formula is: Where, μ s This represents the s-th first statistic; S24. Based on the first statistic μ s Calculate the second statistic α t The calculation formula is: Where, α t This represents the t-th second statistic; S25. Based on the first statistic μ s Second statistic α t Calculate the adaptive filter threshold T y The calculation formula is: T y =μ s +1.5·a t ·(1-Rh c ) Among them, T y This represents the y-th adaptive filtering threshold; S26. Based on the adaptive filtering threshold T y Noise removal; like Then point p l Remove noise and mark the removal point t. T ; like Then point p l These are normal data points; S27. Repeat step S26 until the point cloud dataset is detected and a preliminary preprocessed point set is obtained; S28. Based on the removal point t T Calculate supplementary point B bc The calculation formula is: Among them, t T+1 and t T-1 Representing the rejection point t respectively T Two adjacent points; S29, Add point B bc Add it to the initial preprocessing point set to 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 from the preprocessing point set. d Establish a three-dimensional rectangular coordinate system with the origin, and define any point p in the preprocessed point set. i The coordinates are (x i ,y i ,z i ), and based on the three-dimensional rectangular coordinate system, construct the surface fitting equation N for the preprocessed point set. h The expression is: N h =ax 2 +bxy+cy 2 Where a, b, and c represent the surface fitting equation N, respectively. h The coefficient; 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 This represents the z-axis coordinate value of the preprocessed point; S33, Based on least squares value Q f Calculate the values ​​of coefficients a, b, and c respectively, using the following formula: Among them, Q f This represents the f-th least squares value; S34. Fit 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. Based on the surface function f e Calculate the partial derivatives of (x, y, z), including the first partial derivative f. x1 The second partial derivative value f y2 The third partial derivative value f xy3 The fourth partial derivative value f xx4 and the fifth partial derivative value f yy5 The calculation formula is: in, Represents the surface function f e The partial derivatives of (x,y,z).

5. The adaptive modeling method according to claim 4, characterized in that, Step S4 includes: S41. Based on the first partial derivative value f in the partial derivative values... x1 Second partial derivative value f y2 Calculate the normal vector The calculation formula is: in, Represents the surface function f e The normal vector of (x,y,z); S42. Based on the first partial derivative value f x1 Second partial derivative value 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 These represent the g-th, h-th, and j-th curvature values, respectively. S43. Based on the third partial derivative value f xy3 The fourth partial derivative value f xx4 and the fifth partial derivative value 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 These represent the k-th, l-th, and m-th curvature depth values, respectively.

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 quantity.

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 Coordinates are (A) x B y C z Given a quantity of D, calculate all simplified points pj. t The center of mass O zx The coordinates are calculated using the following formula: Among them, O zx Indicates the coordinates of the x-th centroid; S52. Randomly select a simplified point pj t PC as a reference point t Calculate the reference point and other simplified points pj t Mean distance L between points d The calculation formula is: Among them, L d This represents the distance from the d-th point to the mean. S53, Based on the average point distance L d Calculate the weighting coefficient w h The calculation formula is: Among them, w h This represents the h-th weight coefficient; S54, Based on weighting coefficient w h Calculate the covariance matrix Xc used to determine the axial quantities of the directed bounding box. p The calculation formula is: Where s represents the sum of weight coefficients, and m represents the simplified point pj t Weighted summation; S55. Based on 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 The eigenvalue sequence (λ) is obtained by arranging the values ​​in ascending order. max ,λ max-1 ...λ1,λ0), the eigenvalues ​​λ max , λ max-1 and λ max-2 Set as the axis eigenvalue, and with eigenvalue λ max , λ max-1 and λ max-2 The corresponding feature vectors 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 vectors and the set of reduced points, and calculate the alignment point Dq used to align the reduced 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 This represents the transpose of the d-th axis vector; S62. Define the directed bounding box space as the network region along the axis vector {u1,u2,u3}, based on the mean point distance L. d Calculate the side length Ψ b The calculation formula is: P b =1.5·median(L d ) Among them, Ψ b This represents the length of the b-th partition side; S63, Based on the division of side length Ψ b The network region is divided into multiple sub-regions S. c Calculate the approval value σ for each sub-region. h The calculation formula is: s h =median(L d inS c ) Where, σ h This represents the h-th approved value; S64. Based on the approved value σ h Calculate the adaptive Gaussian kernel φ j (D), the calculation formula is: Where, φ j (D) represents the j-th adaptive Gaussian kernel, Dq k and Dq k-1 These represent the k-th and (k-1)-th alignment points, respectively. S65. Obtain the center point O of each sub-region using the model coordinate method. D and half-length mean Bc z Calculate the global implicit value ▽Yn v The calculation formula is: in, Let g represent 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, Based on the global implicit value ▽Yn v Calculate the eigenvalue η of the isosurface r The calculation formula is: or r =0.5·median(|▽Yn v |) Where, η r This represents the eigenvalue of the r-th isosurface; S72. Based on the isosurface characteristic value η r Constructing three-dimensional digital models of cultural relics using the cube reconstruction method; S73. Define the area of ​​each facet in a three-dimensional digital model of an artifact as A. s The normal vector is Calculate vertex normal vectors The calculation formula is: Where P represents the number of facets; S74, Based on vertex normal vectors Calculate the smoothing intensity Qp h The calculation formula is: Among them, u f Represents the f-th axis vector; S75. Define the vertices d of the three-dimensional digital model of the cultural relic. in The vertex of the 1-ring domain is N. id Calculate the vertex offset Ld of the 3D digital model of the cultural relic. v The calculation formula is: Where I represents the number of vertices in a 1-ring neighborhood; S76, Based on vertex offset Ld v Calculate the smoothing value L of the model ph The calculation formula is: W(x) i =diag(Qp h ) Where Z represents the vertex offset Ld v Quantity; S77, Based on the model smoothing value L ph The surface of the three-dimensional digital model of the cultural relic is smoothed to obtain an optimized three-dimensional model.

Images