A three-dimensional matching and digital restoration method for cultural relic fragments

Abstract

The application discloses a kind of three-dimensional matching and digital repair method of cultural relic fragments, it is related to cultural relic repair technical field, the method comprises the following steps: obtaining the point cloud data of each fragment belonging to the same cultural relic;Using Poisson surface reconstruction method, the triangular mesh model of each fragment is constructed;The corrosion edge in the triangular mesh model of each fragment is repaired;The triangular mesh model of each fragment is registered after edge repair is completed using multi-level matching algorithm, and the optimal splicing position and attitude parameter of each fragment are obtained;According to the optimal splicing position and attitude parameter of each fragment, the triangular mesh model of each fragment after edge repair is completed is spliced, and the initial cultural relic digital model is obtained;The missing area of initial cultural relic digital model is repaired, and the repaired cultural relic digital model is obtained.The application not only improves the accuracy and efficiency of cultural relic repair, but also reduces the risk of manual operation, while ensuring the traceability and reversibility of repair steps.

Description

TECHNICAL FIELD

[0001] The present application relates to the technical field of cultural relic restoration, in particular to a three-dimensional matching and digital restoration method for cultural relic fragments. BACKGROUND

[0002] In the long river of human civilization, archaeological discoveries have always been an important way to explore history and understand ancient societies. The Sanxingdui site, as one of the most remarkable archaeological discoveries of the 20th century, not only reveals the unique cultural features of ancient Shu, but also provides a wealth of bronze artifacts. These artifacts not only have high artistic value, but also carry rich historical information, which is invaluable for studying ancient civilization, cultural heritage, and historical development.

[0003] However, after thousands of years of erosion, most of these artifacts are broken and severely corroded. This not only affects the aesthetic and research value of the artifacts, but also poses great difficulties for their protection and restoration.

[0004] In the field of cultural relic restoration, traditional manual comparison and bonding methods are not only inefficient, but also risk damaging the artifacts, and the restoration results are difficult to ensure consistency. At the same time, although existing digital technologies provide new ideas for cultural relic restoration, high costs, inconvenient equipment, and insufficient precision limit their widespread application in the field of cultural relic restoration. In addition, the large number of bronze fragments, irregular shapes, varying degrees of corrosion, and blurred edges make the restoration work even more difficult. Therefore, how to ensure the historical authenticity and integrity of the artifacts while maintaining the efficiency of the restoration has become a technical problem that needs to be solved in the field of cultural relic restoration. SUMMARY

[0005] The purpose of the present application is to provide a three-dimensional matching and digital restoration method for cultural relic fragments, which can handle different types of corrosion and various specifications of fragments, and achieve high-precision, high-efficiency, and non-damaging cultural relic restoration.

[0006] To achieve the above-mentioned purpose, the present application provides a three-dimensional matching and digital restoration method for cultural relic fragments, comprising:

[0007] Obtaining point cloud data of each fragment belonging to the same cultural relic;

[0008] Based on the point cloud data of each fragment, a Poisson surface reconstruction method is used to construct a triangular mesh model of each fragment;

[0009] Repairing the corroded edges in the triangular mesh model of each fragment to obtain a completed edge-repaired triangular mesh model of each fragment;

[0010] The multi-level matching algorithm is used to register the edge-repaired triangular mesh models of the fragments, so as to obtain the optimal splicing position and optimal attitude parameter of each fragment.

[0011] The edge-repaired triangular mesh models of the fragments are spliced according to the optimal splicing position and optimal attitude parameter of each fragment, so as to obtain an initial digital relic model.

[0012] The missing area of the initial digital relic model is repaired, so as to obtain a repaired digital relic model.

[0013] According to the specific embodiments provided in the application, the application has the following technical effects:

[0014] The application provides a three-dimensional matching and digital repair method for fragments of relics. By obtaining point cloud data of each fragment belonging to the same relic, the application solves the problem that it is difficult to accurately obtain three-dimensional information of fragments of relics in traditional repair methods, and realizes comprehensive digital capture of morphological characteristics of fragments of relics. By constructing a triangular mesh model of each fragment based on point cloud data using a Poisson surface reconstruction method, and repairing eroded edges in the triangular mesh model, the application solves the problems of three-dimensional model construction and edge integrity recovery of fragments of relics, and realizes accurate construction and complete presentation of a three-dimensional model of fragments of relics, thereby providing more accurate boundary conditions for subsequent fragment matching and splicing work. By using a multi-level matching algorithm to register the edge-repaired triangular mesh models, the application solves the problem of accurate matching and splicing of fragments of relics in a three-dimensional space, and realizes high-precision alignment and attitude adjustment between fragments. According to the registration result, the triangular mesh models of the fragments are spliced, so as to obtain an initial digital relic model. This step solves the problem of digital splicing of fragments of relics, and preliminarily realizes digital reconstruction of the whole relic. Finally, by repairing the missing area of the initial digital relic model, the application solves the problem of integrity recovery of the digital relic model, and obtains a repaired digital relic model, thereby realizing high-precision digital reproduction and protection of the original appearance of the relic. It can be seen that the application can not only realize high-precision digital repair of relics, but also significantly improve the efficiency of repair. More importantly, the whole repair process is completely based on a digital model, and does not cause any physical damage to fragments of relics. BRIEF DESCRIPTION OF DRAWINGS

[0015] In order to more clearly illustrate the technical solutions in the embodiments of the application or the prior art, the following will briefly introduce the drawings needed in the embodiments. Obviously, the drawings in the following description are only some embodiments of the application, and for those skilled in the art, other drawings can also be obtained without creative labor.

[0016] Figure 1A flowchart of a three-dimensional matching and digital restoration method of cultural relic fragments is provided for an embodiment of the present application.

[0017] Figure 2 A detailed flowchart of a three-dimensional matching and digital restoration method of cultural relic fragments is provided for an embodiment of the present application. DETAILED DESCRIPTION

[0018] The technical solutions in the embodiments of the present application will be described clearly and completely below with reference to the drawings in the embodiments of the present application. Obviously, the described embodiments are only a part of the embodiments of the present application, rather than all the embodiments. Based on the embodiments in the present application, all other embodiments obtained by those skilled in the art without creative work fall within the scope of protection of the present application.

[0019] In order to make the above-mentioned purposes, features and advantages of the present application more obvious and easy to understand, the present application will be further described in detail below with reference to the drawings and specific embodiments.

[0020] As shown in Figure 1 and Figure 2 , the present application provides a three-dimensional matching and digital restoration method of cultural relic fragments, as shown in the following steps 101 to 106. Among them:

[0021] Step 101, obtaining point cloud data of each fragment belonging to the same cultural relic.

[0022] In an exemplary embodiment, the fragment point cloud data is obtained by a three-dimensional scanning device (3D CaMegaPCP-400). The 3D CaMegaPCP-400 three-dimensional scanner is used to quickly scan each fragment of the same cultural relic, and the scanning frequency is 0.25 seconds. The fragment surface data is obtained by the laser beam emitted by the device irradiating the fragment surface. There are 640 laser points distributed on each laser line, and 480 laser lines are emitted at a time. The collected point cloud data contains more than 300,000 points, including XYZ(three-dimensional coordinates) and RGB(color) information. The scanning accuracy is set to 0.1 millimeter. The modeling accuracy is 0.1 millimeter.

[0023] In another exemplary embodiment, the following technologies can also be used:

[0024] Structured light scanning technology, obtaining three-dimensional information by projecting structured light patterns;

[0025] Photogrammetry technology, reconstructing a three-dimensional model using multi-angle photos.

[0026] Step 102, based on the point cloud data of each fragment, a Poisson surface reconstruction method is used to construct a triangular mesh model of each fragment.

[0027] In an exemplary embodiment, after obtaining the point cloud data of each fragment in the acquisition step 101, the point cloud data is preprocessed, specifically including:

[0028] Noise point removal. Statistical Outlier Removal (SOR) method is used for noise point removal. For each vertex in the point cloud data, the average distance between it and its neighboring points is calculated:

[0029]

[0030] Wherein, μ d represents the average distance between the current vertex and its neighboring points, N' represents the number of neighboring points, d n' represents the distance from the current vertex to the nth' neighboring point.

[0031] Then the standard deviation σ d is calculated according to the following formula:

[0032]

[0033] If the average distance of the current vertex to its neighboring points is out of the threshold range, i.e. μ d < μ d - α' · σ d or μ d > μ d - α' · σ d , the current vertex is marked as a noise point and deleted. Wherein α' is the standard deviation multiple threshold, and α' = 2 in the embodiment.

[0034] Coordinate system unification. The point cloud data in the local coordinate system is converted to the global coordinate system through the rigid transformation matrix T1:

[0035]

[0036] Wherein, R is a 3 × 3 rotation matrix, t is a 3 × 1 translation vector, x1, y1 and z1 represent the horizontal coordinate, vertical coordinate and vertical coordinate of the point cloud data respectively. The optimal transformation matrix is solved by Iterated Closest Points (ICP) algorithm:

[0037]

[0038] Wherein, T * represents the optimal transformation matrix of the point cloud data in the local coordinate system, p b' represents the point cloud in the local coordinate system, q b' represents the corresponding point cloud of p b' in the global coordinate system, and B' represents the total number of effective corresponding point pairs participating in registration.

[0039] Merging of overlapping region point cloud. For the point cloud data of the overlapping region, a voxelization method is used for downsampling merging. The space is divided into a cubic grid with edge length r, and the pixel points falling in the same grid are averaged:

[0040]

[0041] where M' is the number of points in the grid, and in the present embodiment, the voxel size r = 0.5 mm, and p m' represents the point cloud data falling in the same voxel (cubic grid), including the average of XYZ coordinate values (spatial position) and the average of RGB color values (color information), represented by p merged to represent the point cloud data value in each voxel.

[0042] Conversion of grid model. Since the point cloud data is a discrete point set, in order to obtain a continuous surface description and facilitate subsequent texture mapping and other operations, it is converted into a triangular mesh form.

[0043] In this embodiment, the Poisson surface reconstruction method is used to convert the point cloud data of each fragment into a triangular mesh model. First, the normal vector of each vertex is estimated:

[0044]

[0045] where n p represents the normal vector of the current vertex, N p is the neighborhood point set of the current vertex, |N p | is the number of neighborhood points, p is the current vertex, and q is a point in the neighborhood of p, q next represents the next point of q in the neighborhood point set of p. The value of each vertex is p merged , and the normal vector n p of each vertex is obtained. The vector field The exponential function x2 (the target of solving, which defines the object surface, the isosurface with a value of 0, i.e., the surface composed of all points with x2 = 0, where x2 > 0 inside the object; x2 < 0 outside the object; x2 = 0 on the object surface) can be calculated according to the Poisson equation.

[0046] Then the Poisson equation is solved:

[0047]

[0048] where, is the Laplacian operator, represents the vector field p composed of the normal vector n divergence. That is, an isosurface with value 0 can be obtained. Then the zero isosurface is triangulated using the Marching Cubes algorithm: first, the space is divided into a regular cubic grid; second, for each cube, the scalar values (indicator function values) of the 8 vertices are checked to determine whether the isosurface intersects the cube; third, the intersection points are calculated by interpolation on the intersecting edges; and fourth, the intersection points are connected to form triangles, thereby obtaining a triangular mesh model.

[0049] Specifically, the first step is space division and scalar value calculation. First, the space region containing the point cloud needs to be divided into a regular cubic grid, and each cube has 8 vertices. It should be noted that these vertices are not the original point cloud data points, but the vertices of the regular grid. At each grid vertex, we need to calculate the value of the indicator function x2, which is obtained by solving the Poisson equation. When the vertex is inside the object, the value of x2 is greater than 0; when the vertex is outside the object, the value of x2 is less than 0.

[0050] The second step is to determine the intersection of the isosurface. Check each edge of the cube. If the scalar values of the two endpoints of an edge are positive and negative, respectively, it means that the isosurface (i.e., the x2 = 0 surface) intersects the edge. This determination process determines the position of the object surface in space.

[0051] The third step is intersection point calculation. For those edges that are determined to intersect the isosurface, the specific intersection point position of the isosurface is calculated by linear interpolation. Suppose the two endpoints of the edge are v1 and v2, and their corresponding scalar values are s1 and s2, then the intersection point position can be calculated by the formula v1 + t(v2-v1), where t = s1 / (s1-s2). These calculated intersection points are new generated points, which are not the original point cloud data points.

[0052] The fourth step is the generation of triangular mesh. According to the distribution of intersection points in the cube, use the predefined triangular connection mode (usually there is a lookup table of 256 cases) to connect the intersection points into triangles, forming local surface patches. All the triangular patches of the cubes together form the complete object surface.

[0053] In the specific application to debris, the original point cloud data is mainly used in three aspects: calculating the normal vector field, guiding the solution of the Poisson equation, and providing the approximate position information of the surface. The vertices of the final mesh model come from the intersection points calculated by the Marching Cubes (MC) algorithm, which are optimized to better represent the object surface, rather than directly using the original point cloud data points as mesh vertices.

[0054] The above method has multiple advantages: it can well handle noise and incomplete data; it can generate a continuous and smooth surface; it can naturally handle complex topological structures; and it can also fill small data gaps. Overall, the final triangular mesh model is a continuous surface reconstructed based on point cloud data, rather than directly using point cloud data points as mesh vertices, which can achieve better surface reconstruction results.

[0055] In this implementation, the triangular mesh can represent a continuous surface with fragments. Since there may be small holes in the corroded bronze ware, for holes in the triangular mesh model with an area less than a threshold S, a moving least squares method is used to fit a surface to fill the holes:

[0056]

[0057] where w m” (x') is a weight function:

[0058] w m” (x') = exp(-||x'-x' m” || 2 / h 2 ).

[0059] where m" represents the number of holes with an area less than a threshold S, h is a smoothing parameter, x' is a spatial position point in the region to be repaired, f(x') is the value of the fitted surface function at position x', and x' m” are known points on the boundary of the hole. This process uses information from surrounding known points to estimate the surface shape of the hole region through weighted averaging. The final output is a completed continuous surface patch that seamlessly integrates into the original mesh model.

[0060] In this implementation, the Quadric Error Metric (QEM) method is used for mesh simplification. For each edge, its cost is calculated:

[0061]

[0062] where F e is the set of faces adjacent to edge e, and K f is the quadric error matrix of face f(x'). The edges are gradually contracted from small to large cost until the target number of faces is reached.

[0063] Through the above steps, noise points are removed, coordinates are unified, local coordinate systems are converted to a unified coordinate system and overlapping area point cloud merging is processed, mesh model conversion, mesh hole repair, and mesh model simplification are performed. High-precision three-dimensional data is obtained to provide basic data support for subsequent processing.

[0064] Step 103, repairing the corroded edges in the triangular mesh model of each fragment to obtain the triangular mesh model of each fragment with completed edge repair.

[0065] In some embodiments, the edges of the corroded bronze fragment may be unclear and need to be repaired. In order to better match the fragments later, step 103 repairs the corroded edges in the triangular mesh model of each fragment to obtain the triangular mesh model of each fragment with completed edge repair, which specifically includes steps 201-203:

[0066] Step 201, performing edge curvature analysis on the triangular mesh model of the fragment to obtain the curvature analysis result of each vertex in the triangular mesh model.

[0067] Step 202, determining the type of each vertex in the triangular mesh model based on the edge score of each vertex in the triangular mesh model and the adaptive threshold, the type being a complete edge point, a corroded edge point and a transition area point.

[0068] Step 203, repairing the vertices with the type of corroded edge point in the triangular mesh model in a surface extrapolation manner to obtain the triangular mesh model of the fragment with completed edge repair.

[0069] In this implementation, the curvature analysis result in step 201 at least includes: average curvature and maximum principal curvature.

[0070] For each triangular facet T with p as a vertex, the area of the Voronoi region belonging to p in the facet is calculated: if the triangle is an obtuse triangle: the area is taken as half of the area of the triangle; if the triangle is an acute triangle: the Voronoi region is surrounded by the perpendicular bisector of the edge.

[0071] The calculation formula of the average curvature is:

[0072]

[0073] wherein, A Mixed (p) represents the Mixed Voronoi area of the current vertex, p represents the current vertex, p1 and p2 represent the first and second adjacent vertices of the current vertex respectively, T represents a triangle with p, p1 and p2 as vertices, S T represents the area of the triangle, a and a" represent the angles with the first and second adjacent vertices as vertices respectively, H(p) represents the average curvature of the current vertex, N(p) represents the vertex set in the preset neighborhood of the current vertex p, a j and β jrespectively represent an angle with the current vertex as a vertex and an angle with the second adjacent vertex as a vertex, p j represents the three-dimensional coordinates of the jth vertex in the neighborhood.

[0074] The calculation formula of the maximum principal curvature is:

[0075]

[0076] wherein K(p) represents the Gaussian curvature of the current vertex, θ i represents the internal angle of the adjacent triangular facet at the current vertex, n is the number of triangular facets adjacent to the current vertex, and k1(p) represents the maximum principal curvature of the current vertex.

[0077] wherein in the implementation of this embodiment, the calculation formula of the edge score in step 202 is:

[0078] EdgeScore(p)=w1·ΔH(p)+w2·||k1(p)||+w3·D(p).

[0079] wherein EdgeScore(p) represents the edge score of the current vertex, w1, w2 and w3 are respectively a first weight coefficient, a second weight coefficient and a third weight coefficient, ΔH(p) represents the average curvature change rate, and D(p) represents the distance from the current vertex to the nearest non-edge point.

[0080] The calculation formula of the adaptive threshold is:

[0081] Thresh local (p)=OTSU(EdgeScore(N k (p)));

[0082] wherein Thresh local (p) represents the adaptive threshold of the current vertex, N k (p) represents a near-neighbor region formed by k nearest neighbors of the current vertex as a center, k=30, and OTSU represents the maximum inter-class variance method.

[0083] The type of each vertex in the triangular mesh model is determined according to the following conditions:

[0084] If EdgeScore(p)<0.8×Thresh local (p), the current vertex is a complete edge point.

[0085] If EdgeScore(p)>1.2×Thresh local (p), the current vertex is an eroded edge point.

[0086] If 0.8×Thresh local (p) and 1.2×Thresh (p), the current vertex is a partial edge point.(p) < EdgeScore(p) < 1.2 x Thresh local (p), the current vertex is a transition region point.

[0087] Wherein, the step 203 specifically comprises:

[0088] For each eroded edge point, a quadratic surface is fitted in its k-nearest neighbor region (k is set to 20 by default):

[0089] f(x, y) = ax 2 + by 2 + cxy + dx + ey + f.

[0090] The coefficients are solved by using the weighted least squares method:

[0091] w r' = exp(-||p r' -p|| 2 / h 2 ).

[0092] Wherein, p r' represents the coordinates of the nearest neighbor points in the k-nearest neighbor region of the current eroded edge point, h is a smoothing parameter, which is 1 / 3 of the radius of the nearest neighbor region, and w r' represents the weight coefficient, which is calculated in advance according to the distance and determines the importance of each data point in solving the surface coefficients, and is used to determine the influence of each nearest neighbor point on fitting, so as to help solve the to-be-solved surface coefficients a, b, c, d, e, and f by using the weighted least squares method.

[0093] Specifically, we need to minimize the weighted sum of squared errors to solve the coefficients. The error function can be expressed as:

[0094]

[0095] Wherein, x r' , y r' and z r' represent the horizontal coordinate, vertical coordinate and vertical coordinate of p r' respectively. In order to solve the surface coefficients a, b, c, d, e, f, the error function F error needs to be minimized. This can be achieved by taking the partial derivative of F error with respect to each surface coefficient and setting the partial derivative to zero. This will result in a linear equation system, which can be solved by solving the equation system to find the optimal value of the surface coefficients. Specifically, the partial derivatives of F error with respect to a, b, c, d, e, f are taken, and these partial derivatives are set to zero:

[0096]

[0097] A linear equation set containing six equations is obtained, and the values of the surface coefficients a, b, c, d, e and f can be found by solving the equation set.

[0098] The formula for repairing the vertex of the type of corrosion edge point in the triangular mesh model is:

[0099]

[0100] wherein p smooth represents the reconstructed corrosion edge point after repair, p new represents the reconstructed corrosion edge point, p new =f(x+δ x ,y+δ y ), x and y represent the horizontal coordinate of the corrosion edge point and the vertical coordinate of the corrosion edge point respectively, δ x and δ y represent the extrapolation step length of the horizontal coordinate of the corrosion edge point and the extrapolation step length of the vertical coordinate of the corrosion edge point respectively, N(p new ) represents the set of neighboring points in the preset neighborhood of the reconstructed corrosion edge point, f() represents the fitted quadratic surface of the current reconstructed corrosion edge point, p i' represents the i'th vertex in N(p new ), G s represents the spatial domain Gaussian kernel, G r represents the value domain Gaussian kernel, G s (||p new -p i' ||) represents the spatial distance weight of the reconstructed corrosion edge point and the neighboring point, G r (|H(p new )-H(p i' )|) represents the average curvature weight of the reconstructed corrosion edge point and the neighboring point.

[0101] The spatial domain Gaussian kernel G s is used to measure the spatial distance relationship between the reconstructed corrosion edge point and its neighboring points, and the weight is larger when the distance is closer. The spatial position relationship of the reconstructed corrosion edge point is considered through the spatial domain Gaussian kernel, and the local geometric shape is maintained in combination with the fitted quadratic surface f(), so that the geometric features of the reconstructed corrosion edge point are repaired.

[0102] The value domain Gaussian kernel G r is used to measure the average curvature difference between the reconstructed corrosion edge point and its neighboring points, and the weight is larger when the curvature difference is smaller. The local texture change is maintained through the average curvature difference, so that the texture features of the reconstructed corrosion edge point are repaired.

[0103] Specifically, for each corrosion edge point, firstly, surface extrapolation is performed: taking the current corrosion edge point as the center, k' = 20 nearest neighbor points are selected, a quadratic surface f() is fitted using these points, and extrapolation is performed in the x and y directions to obtain new points p new .

[0104] Then, weighted smoothing based on the bilateral filtering idea is performed:

[0105] In the preset neighborhood of p new , a vertex set N(p new ) is collected, a spatial domain weight Gs and a value domain weight Gr are applied to each neighborhood point, and a final repair point Psmooth is obtained through weighted averaging.

[0106] Wherein, Gs(x) = exp(-x 2 / 2σs 2 ), Gr(x) = exp(-x 2 / 2σr 2 ), σs is a spatial domain standard deviation, and σr is a value domain standard deviation.

[0107] In step 104, a multi-level matching algorithm is used to register the triangular mesh models of each fragment with the completed edge repair to obtain the optimal splicing position and optimal attitude parameter of each fragment.

[0108] In some embodiments, through multi-feature fusion and multi-level matching strategy, accurate automatic matching of fragments is realized, and the consistency of the overall splicing is ensured through global optimization. In step 103, a multi-level matching algorithm is used to register the triangular mesh models of each fragment with the completed edge repair to obtain the optimal splicing position and optimal attitude parameter of each fragment, which specifically includes steps 301-304:

[0109] In step 301, feature extraction is performed based on the triangular mesh models of each fragment with the completed edge repair.

[0110] In step 302, the weighted similarity between any two fragments is calculated according to the features extracted from each fragment to obtain a first matching screening result.

[0111] In step 303, a nonlinear optimization objective function is constructed based on the first matching screening result, and the nonlinear optimization objective function is iteratively optimized until a preset matching precision is reached to obtain a second matching screening result.

[0112] In step 304, a global consistency optimization algorithm is used to adjust the splicing position and attitude parameter of each fragment in the second matching screening result to obtain the optimal splicing position and optimal attitude parameter of each fragment.

[0113] In the embodiment, the feature extraction in step 301 comprises at least geometric feature extraction and texture feature extraction.

[0114] The geometric feature extraction comprises constructing a curvature distribution descriptor of each vertex according to curvature analysis structure, and extracting a broken edge curve feature of each patch based on the completed edge-repaired triangular mesh model of each patch.

[0115] The expression of the curvature distribution descriptor is as follows:

[0116] C dist (p) = {k1(q), k2(q) | q ∈ N r (p)}.

[0117]

[0118] Wherein, C dist (p) represents the curvature distribution descriptor of the current vertex, N r (p) represents the neighborhood point set with the current vertex as the center and r1 as the radius, r1 is equal to 2 times the average curvature, k2() represents the second principal curvature value, q represents the vertex in the neighborhood point set.

[0119] The expression of the broken edge curve feature is as follows:

[0120] E feat (s) = [κ(s), τ(s), n(s), θ(s)].

[0121] Wherein, E feat (s) represents the broken edge curve feature of the current patch, s is the broken edge parameter of the current patch, κ(s) is the broken edge curvature of the current patch, τ(s) is the broken edge torsion of the current patch, n(s) is the broken edge surface normal vector of the current patch, and θ(s) is the broken surface direction angle of the current patch.

[0122] The texture feature extraction comprises calculating a broken surface normal map feature of each patch based on the completed edge-repaired triangular mesh model of each patch, and constructing a local height field descriptor according to a preset shape field of each vertex.

[0123] The calculation formula of the broken surface normal map feature is as follows:

[0124]

[0125] Wherein, N v map represents the broken surface normal map feature of the current patch, (u, v) represents the parametric plane coordinates of the current patch, W(u, v) represents the vertex region of the current patch, n qrepresents a vertex normal vector within the vertex region of the current fragment.

[0126] The expression of the local height field descriptor is:

[0127] H desc (p) = {h(q) | q e N d (p)}.

[0128] where H desc (p) represents the local height field descriptor of the current vertex, N d (p) represents a circular neighborhood with the current vertex as the center and a diameter of d, d takes 5 times of the average edge length of the mesh, q represents a vertex within the circular neighborhood, and h(q) represents the relative height value of the qth vertex within the circular neighborhood relative to the current vertex.

[0129] where implementing this embodiment, step 302 specifically includes:

[0130] The weighted similarity of any two fragments f and f' is calculated according to the following formula:

[0131]

[0132] where S ff' represents the weighted similarity of any two fragments, F k” represents the k"th feature, k" = 1, 2, 3, 4 respectively correspond to the curvature distribution descriptor, the broken edge curve feature, the broken surface normal map feature, and the local height field descriptor, w k” represents the weight of the k"th feature, and sim() is a similarity measurement function representing any two fragments.

[0133] If the weighted similarity of any two fragments is greater than a preset similarity, the current any two fragments are added to the feature matching point pair set.

[0134] The fragments of the feature matching point pair set are subjected to first matching screening using the following formula:

[0135]

[0136] where T init represents the first matching screening result, R represents a rotation matrix, t represents a translation vector, M represents the feature matching point pair set, Rxf represents the rotated coordinates of the fragment f, and (f, f') represents the fragment pair in the feature matching point pair set.

[0137] where implementing this embodiment, the construction formula of the nonlinear optimization objective function in step 303 is:

[0138] E(T init ) = E geo(T init )+λ1E norm (T init )+λ2E curv (T init )。

[0139] wherein, E(T init ) represents a value of a non-linear optimization objective function constructed based on the first matching screening result, i.e., the second matching screening result E, Egeo represents a geometric distance item, Enorm represents a normal consistency item, Ecurv represents a curvature continuity item, and λ1 and λ2 represent a first weight coefficient and a second weight coefficient respectively.

[0140] wherein, in the implementation, the step 304 specifically comprises:

[0141] constructing a matching graph G(V, E) according to the second matching screening result; wherein, V is a fragment node set, and E is a matching edge set.

[0142] finding a minimum weight edge set connecting all fragments using the following formula:

[0143] w(e) = -log(S ff' ).

[0144] MST = argmin∑ e∈E w(e).

[0145] wherein, w(e) represents a weight of a matching edge, and MST represents a minimum weight edge set connecting all fragments.

[0146] constructing an initial global transformation matrix for each fragment based on the minimum weight edge set, and adjusting a stitching position and a pose parameter of each fragment in the second matching screening result using the following formula until a preset iteration number is reached, to obtain an optimal stitching position and an optimal pose parameter of each fragment:

[0147]

[0148] wherein, E global represents a global error, (j, j') represents a fragment pair in the second matching screening result, T j and T j' represent global transformation matrices of the fragments j and j' respectively, and T jj' represents a relative transformation matrix of the fragment pair (j, j').

[0149] stitching the triangular mesh models of each fragment with completed edge repair according to the optimal stitching position and the optimal pose parameter of each fragment, to obtain an initial digital model of cultural relics.

[0150] Step 106, repairing the missing area of the initial cultural relic digital model to obtain a repaired cultural relic digital model.

[0151] In some embodiments, the reconstruction reference is provided by the feature database, the accurate reconstruction of the geometric shape is realized, the natural transition and continuity of the texture are ensured, the complete digital model conforming to the historical features is generated, step 106 repairs the missing area of the initial cultural relic digital model to obtain a repaired cultural relic digital model, and specifically includes steps 401-403:

[0152] Step 401, determining the missing area in the initial cultural relic digital model according to the missing boundary curve.

[0153] Step 402, using a B-spline method to perform surface reconstruction on the missing area to obtain a reconstructed area.

[0154] Step 403, synthesizing the texture of the reconstructed area based on the texture feature library to obtain a repaired cultural relic digital model.

[0155] In this implementation, step 401 further includes the following:

[0156] Extracting the boundary points of the missing area:

[0157] Wherein, C boundary represents the boundary point set of the missing area, represents the boundary of the missing area, p b” is the boundary point of the missing area.

[0158] The constraint condition of the missing area boundary is calculated by the following formula:

[0159] B constraint ={n(s),κ(s),θ(s)}.

[0160] Wherein, B constraint represents the boundary constraint condition set, s represents the current boundary point in the missing area boundary point set, n(s) represents the normal vector of the current boundary point, κ(s) represents the boundary curvature of the current boundary point, and θ(s) represents the tangent angle of the current boundary point.

[0161] In this implementation, step 402 specifically includes the following:

[0162] Step 501, based on the matched fragments, extracting the following features for the construction of the geometric feature database:

[0163] Global shape feature. Calculate the overall shape descriptor:

[0164] F global ={V',S',R',P'}.

[0165] where V' is the volume proportion feature, S' is the surface area distribution, R' is the radial distribution function, and P' is the principal direction feature.

[0166] Local feature template. Extract local feature template from complete region:

[0167] F local (p) = {H(p), K(p), D(p)}.

[0168] where H(p) is the mean curvature field, K(p) is the Gaussian curvature field, and D(p) is the shape diameter function.

[0169] Step 502, establish texture feature library based on known region:

[0170] T pattern = {T color , T normal , T rough}.

[0171] where T color is the color distribution pattern, T normal is the normal distribution pattern, and T rough is the roughness distribution pattern.

[0172] Step 503, construct initial control point grid to obtain reconstructed region by the following formula:

[0173]

[0174] where S(u', v') represents the reconstructed region, (u', v') represents the parametric plane coordinates of the control points in the reconstructed region, P i”,j” is the control point in the reconstructed region, N i”,k' (u') and N j”,l (v') represent the first and second B-spline basis functions, respectively, and k' and l' represent the spline order of the first and second B-spline basis functions, respectively.

[0175] Step 504, boundary optimization of the reconstructed region by the following formula:

[0176] E surface = E smooth + λ3E boundary + λ4E feature .

[0177] wherein, Esmooth represents a smooth term of the reconstruction region, Eboundary represents a boundary constraint term of the reconstruction region, Efeature represents a feature preserving term of the reconstruction region, and λ3 and λ4 represent the third weight coefficient and the fourth weight coefficient respectively.

[0178] In another exemplary embodiment, step 402 can also obtain the reconstruction region by using the symmetry feature of the cultural relic; or reconstruct based on a similar cultural relic template.

[0179] In the implementation, step 403 specifically includes:

[0180] The reconstruction region is parameterized by using the following formula:

[0181]

[0182] wherein, E' represents a vertex set of the reconstruction region, i" and j" represent the horizontal coordinate and the vertical coordinate of the current vertex in the reconstruction region respectively, u i” and u j” represent the parametric domain horizontal coordinate and the parametric domain vertical coordinate of the current vertex in the reconstruction region respectively, d i”j” represents the geodesic distance of the current vertex in the reconstruction region, and w i”j” represents the edge weight of the current vertex in the reconstruction region.

[0183] The parameterized reconstruction region is matched with the texture element according to the following formula:

[0184]

[0185] wherein, p' represents the current vertex in the reconstruction region, T new represents a texture similar to the boundary texture of the reconstruction region in the texture feature library, T pattern represents the texture feature library, T boundary represents the boundary texture of the reconstruction region, T represents a texture in the texture feature library, and D() represents a texture similarity measurement function; T new is obtained by matching the similarity between the texture feature library and the boundary texture, and is used as an ideal texture state that the internal vertex of the reconstruction region should reach.

[0186] The texture of the reconstruction region is obtained by minimizing the texture transition error by using the following formula, and the repaired digital model of the cultural relic is obtained:

[0187]

[0188] wherein, E texture represents the texture transition error, Ω represents the reconstruction region, p" represents the current vertex in the boundary region of the reconstruction region, represents the boundary region of the reconstruction region, and Tknown T represents the texture of the non-missing area, and μ represents the boundary weight coefficient of the reconstructed area.

[0189] The original area on the cultural relic fragment that is not damaged contains complete original texture information T known T is used as a reference for the synthesis of the texture of the reconstructed area. known T is used as a reference for the synthesis of the texture of the reconstructed area.

[0190] The whole process can be understood as follows: first, the reconstructed area is parameterized, a unified texture coordinate system is established, then the texture sample that best matches the boundary texture of the reconstructed area is selected from the texture feature library, and finally the continuity of the texture inside the reconstructed area and the smooth transition of the known area boundary are realized by minimizing the texture transition error, which can ensure that the repaired texture not only maintains the original features, but also realizes the natural transition.

[0191] In summary, the method of the present application is designed by modularization, and each functional module is independently encapsulated and equipped with a standardized data interface. Not only can each module work efficiently and collaboratively, but also greatly facilitates the maintenance and upgrading of the system. When a module needs to be updated or optimized, it can be done individually without the need for large-scale changes to the entire system, thereby reducing maintenance costs and improving the flexibility and scalability of the system. The use of data pipeline architecture makes the processing flow clear and efficient, and the processing process is traceable. Under this architecture, each link from input to output is clearly defined and transmitted through a standardized data interface. This not only ensures the accuracy and integrity of the data, but also enables the processing process to be effectively monitored and recorded, providing strong support for subsequent data analysis and problem troubleshooting. By integrating intelligent processing technology and with the support of deep learning algorithms, adaptive parameter adjustment and multi-feature fusion matching are achieved. This enables the system to automatically adjust processing parameters according to different situations, thereby optimizing processing results. At the same time, the application of multi-feature fusion matching technology also improves the accuracy and robustness of matching, providing more reliable technical support for the three-dimensional matching and digital restoration of cultural relic fragments.

[0192] In another exemplary embodiment, in addition to modular design, cloud computing architecture can be used to place computing tasks in the cloud or use distributed systems for multi-node parallel processing.

[0193] The application also provides an application scene of the three-dimensional matching and digital restoration method of cultural relic fragments. Specifically, the three-dimensional matching and digital restoration method of cultural relic fragments provided by the embodiment can be applied in a comprehensive scene of cultural relic protection and digital restoration. The comprehensive scene includes a preliminary processing link, a three-dimensional digital restoration link, and an application and display link of restoration results. The cultural relic fragments first enter the preliminary processing link, and after being cleaned, classified, and preliminarily recorded, enter the three-dimensional digital restoration link. In this link, the three-dimensional matching and digital restoration method of cultural relic fragments provided by the embodiment uses technical means such as high-precision scanning, three-dimensional modeling, edge restoration, and multi-level matching algorithms to accurately reconstruct and digitally restore the cultural relic fragments. After the restoration is completed, the digital model of the cultural relic enters the application and display link of restoration results, which is used for cultural heritage protection, research, education, and public display, and the like. The three-dimensional matching and digital restoration method of cultural relic fragments provided by the embodiment is a key technology in the three-dimensional digital restoration link. Specifically, in this link, efficient and accurate restoration of cultural relic fragments can be achieved, which not only improves the efficiency and quality of cultural relic restoration, but also provides strong technical support for subsequent cultural relic restoration, display, and dissemination.

[0194] The application has the following technical effects:

[0195] By using a high-precision three-dimensional scanning device to obtain original data, more than 300,000 point cloud data are collected each time, edge detection and restoration are directly performed in three-dimensional space, information loss caused by 2D-3D conversion is avoided, an edge feature extraction method based on curvature analysis can accurately capture local geometric details, a matching strategy of multi-feature fusion is used, geometric and texture features are comprehensively considered, and matching accuracy is improved. The restoration accuracy of the application reaches more than 95%, and the scanning accuracy can reach 0.1 mm.

[0196] By using an automatic three-dimensional scanning system, complete data can be obtained in only 0.25 seconds per scan; by using an improved feature extraction algorithm, rapid edge detection and restoration are realized; according to a multi-level matching strategy, the convergence speed is accelerated through gradual optimization from coarse registration to fine registration; and based on an intelligent reconstruction method of a feature database, manual intervention is reduced, and the restoration speed is significantly improved. The average processing time of a single fragment is shortened to 30 minutes by the application, and the efficiency is improved by 80% compared with the traditional manual method.

[0197] The application uses laser three-dimensional scanning technology, does not need to touch the surface of cultural relics, performs fragment matching and restoration verification in a virtual environment, establishes a digital cultural relic feature database, reduces repeated operations on real objects, uses a completely non-contact digital processing method, and avoids secondary damage to cultural relics.

[0198] The application preserves original feature information in the edge repairing process through a complete three-dimensional data acquisition and storage mechanism, records all pose transformation parameters in the matching process, and realizes digital recording of the whole repairing process through scientific management and storage of the feature database, thereby ensuring traceability and reversibility of the repairing process.

[0199] The application can process different types of corrosion and various specifications of fragments through a direct feature analysis method in three-dimensional space, a matching strategy based on multi-feature fusion, improved robustness of matching, and an intelligent reconstruction method based on a database that can handle various missing conditions.

[0200] The technical features of the above embodiments can be combined in any manner, and to make the description concise, not all possible combinations of the technical features in the above embodiments are described, but as long as the combinations of the technical features do not exist contradictions, they should be considered as the scope of the present application.

[0201] The principles and implementation modes of the application are described herein by using specific examples, and the above embodiment descriptions are only used to help understand the method of the application and its core idea; meanwhile, for those skilled in the art, the specific implementation modes and application ranges will be changed according to the idea of the application. In conclusion, the content of the present application should not be understood as a limitation.

Claims

1. A method for three-dimensional matching and digital restoration of fragments of cultural heritage objects, characterized by, The three-dimensional matching and digital restoration method of the cultural relic fragments comprises the following steps: Obtain point cloud data of each fragment belonging to the same cultural relic; Based on the point cloud data of each fragment, a triangular mesh model of each fragment is constructed by using a Poisson surface reconstruction method; Repair the corroded edges in the triangular mesh model of each fragment to obtain a completed edge-repaired triangular mesh model of each fragment, specifically including: performing edge curvature analysis on the triangular mesh model of the fragment to obtain a curvature analysis result of each vertex in the triangular mesh model; determining the type of each vertex in the triangular mesh model based on the edge score of each vertex in the triangular mesh model and an adaptive threshold, the type being a complete edge point, a corroded edge point and a transition area point; repairing the vertex with the type of the corroded edge point in the triangular mesh model by using a surface extrapolation method to obtain the completed edge-repaired triangular mesh model of the fragment; the curvature analysis result at least includes: average curvature and maximum principal curvature; Align the completed edge-repaired triangular mesh models of each fragment by using a multi-level matching algorithm to obtain the optimal splicing position and the optimal attitude parameter of each fragment; Splice the completed edge-repaired triangular mesh models of each fragment according to the optimal splicing position and the optimal attitude parameter of each fragment to obtain an initial cultural relic digital model; Repair the missing area of the initial cultural relic digital model to obtain a repaired cultural relic digital model.

2. The method for three-dimensional matching and digital restoration of fragments of an artifact according to claim 1, characterized in that, The calculation formula of the average curvature is: wherein A Mixed (p) represents the mixed Thiessen polygon area of the current vertex, p represents the current vertex, p1 and p2 represent the first and second adjacent vertices of the current vertex respectively, T represents a triangle with p, p1 and p2 as vertices, S T represents the area of the triangle, a and a' represent the angles with the first and second adjacent vertices as vertices respectively, H(p) represents the average curvature of the current vertex, N(p) represents the vertex set within the preset neighborhood of the current vertex p, a j and β j represent the angles with the current vertex and the second adjacent vertex as vertices respectively, p j represents the three-dimensional coordinates of the jth vertex within the neighborhood; The calculation formula of the maximum principal curvature is: where K(p) denotes the Gaussian curvature of the current vertex, θ i denotes the internal angle of the neighboring triangle at the current vertex, n is the number of neighboring triangles of the current vertex, and k1(p) denotes the maximum principal curvature of the current vertex.

3. The method for three-dimensional matching and digital restoration of fragments of an artifact according to claim 2, characterized in that, The calculation formula of the edge score is: EdgeScore(p)=w1·ΔH(p)+w2·||k1(p)||+w3·D(p); Wherein, EdgeScore(p) represents the edge score of the current vertex, w1, w2 and w3 are respectively the first weight coefficient, the second weight coefficient and the third weight coefficient, ΔH(p) represents the average curvature change rate, D(p) represents the distance from the current vertex to the nearest non-edge point; The calculation formula of the adaptive threshold is: Thresh local (p) = OTSU(EdgeScore(N k (p))); wherein Thresh local (p) denotes the adaptive threshold value of the current vertex, N k (p) denotes the neighborhood region formed by k nearest neighbors of the current vertex, k = 30, OTSU denotes the maximum inter-class variance method; The type of each vertex in the triangular mesh model is determined according to the following conditions: If EdgeScore(p) < 0.8 x Thresh local (p), then the current vertex is a complete edge point. If EdgeScore(p) > 1.2 x Thresh local (p), then the current vertex is an eroded edge point. 0.8 x Thresh local (p) < EdgeScore(p) < 1.2 x Thresh local (p), then the current vertex is a transition region point.

4. The method for three-dimensional matching and digital restoration of fragments of cultural heritage objects according to claim 1, characterized in that, The formula for repairing the vertex with the type of the corroded edge point in the triangular mesh model is: wherein p smooth represents a reconstructed corrosion edge point of which repair is completed, p new represents a reconstructed corrosion edge point, p new = f(x+ δ x , y+ δ y ), x and y respectively represent the horizontal coordinate of the corrosion edge point and the vertical coordinate of the corrosion edge point, δ x and δ y respectively represent the extrapolation step length of the horizontal coordinate of the corrosion edge point and the extrapolation step length of the vertical coordinate of the corrosion edge point N(p new ) represents a set of near neighbor points within a preset neighborhood of the reconstructed corrosion edge point, f() represents a fitting quadratic surface of the current reconstructed corrosion edge point, p i' represents the i'th vertex in N(p new ), G s represents a spatial domain Gaussian kernel, G r represents a value domain Gaussian kernel, G s (||p new -p i' ||) represents a spatial distance weight of the reconstructed corrosion edge point and the near neighbor points, G r (|H(p new )-H(p i' )|) represents an average curvature weight of the reconstructed corrosion edge point and the near neighbor points.

5. The method for three-dimensional matching and digital restoration of fragments of cultural heritage objects according to claim 2, characterized in that, Align the completed edge-repaired triangular mesh models of each fragment by using a multi-level matching algorithm to obtain the optimal splicing position and the optimal attitude parameter of each fragment, specifically including: Extract features based on the completed edge-repaired triangular mesh models of each fragment; Calculate the weighted similarity between any two fragments according to the extracted features of each fragment to obtain a first matching screening result; Construct a nonlinear optimization objective function based on the first matching screening result, and iteratively optimize the nonlinear optimization objective function until a preset matching accuracy is reached to obtain a second matching screening result; Adjust the splicing position and the attitude parameter of each fragment in the second matching screening result by using a global consistency optimization algorithm to obtain the optimal splicing position and the optimal attitude parameter of each fragment.

6. The method for three-dimensional matching and digital restoration of fragments of cultural heritage objects according to claim 5, characterized in that, The feature extraction at least includes: geometric feature extraction and texture feature extraction; The geometric feature extraction is: constructing a curvature distribution descriptor of each vertex according to a curvature analysis result, and extracting a fracture edge curve feature of each fragment based on the completed edge-repaired triangular mesh model of each fragment; an expression of the curvature distribution descriptor is: C dist (p) = {k1(q), k2(q) | q e N r (p)}; where C dist (p) denotes the curvature distribution descriptor of the current vertex, N r (p) denotes the set of neighboring points centered at the current vertex with radius r equal to 2 times the average curvature, k2() denotes the second principal curvature value, q denotes a vertex within the set of neighboring points; An expression of the fracture edge curve feature is: E feat (s) = [K(s), T(s), n(s), Q(s)]; wherein E feat (s) is a fracture edge curve feature of the current fragment, s is a fracture edge parameter of the current fragment, K(s) is a fracture edge curvature of the current fragment, T(s) is a fracture edge torsion of the current fragment, n(s) is a fracture edge surface normal vector of the current fragment, and Q(s) is a fracture surface direction angle of the current fragment. The texture feature extraction is: calculating a fracture surface normal map feature of each fragment based on the completed edge-repaired triangular mesh model of each fragment, and constructing a local height field descriptor according to a preset shape neighborhood of each vertex; A calculation formula of the fracture surface normal map feature is: wherein, represents a normal map feature of the break plane of the current patch, (u, v) represents the parametric plane coordinates of the current patch, W(u, v) represents the vertex region of the current patch, n q represents a vertex normal vector within the vertex region of the current patch; An expression of the local height field descriptor is: H desc (p) = {h(q) | q e N d (p)}; where H desc (p) denotes the local height field descriptor of the current vertex, N d (p) denotes the circular neighborhood with the current vertex as the center and the diameter d, d takes 5 times the average edge length of the mesh, q denotes the vertex in the circular neighborhood, h(q) denotes the relative height value of the qth vertex in the circular neighborhood relative to the current vertex.

7. The method for three-dimensional matching and digital restoration of fragments of cultural heritage objects according to claim 5, characterized in that, According to the features extracted from each fragment, a weighted similarity between any two fragments is calculated, and a first matching screening result is obtained, specifically including: The weighted similarity between any two fragments f and f' is calculated according to the following formula: where S ff' denotes the weighted similarity of any two fragments, F k” denotes the k”-th feature, k” = 1, 2, 3, 4 corresponding to curvature distribution descriptor, fracture edge curve feature, fracture surface normal map feature and local height field descriptor respectively, w k” denotes the weight of the k”-th feature, sim() is a similarity measure function of any two fragments; If the weighted similarity between the current any two fragments is greater than a preset similarity, the current any two fragments are added to a feature matching point pair set; The fragments of the feature matching point pair set are screened by the following formula: where T init denotes the first matching result, R denotes a rotation matrix, t denotes a translation vector, M denotes a set of feature matching point pairs, Rxf denotes a rotated coordinate of fragment f, and (f, f') denotes a pair of fragments in the set of feature matching point pairs.

8. The method for three-dimensional matching and digital restoration of fragments of cultural heritage objects according to claim 1, characterized in that, The missing area of the initial cultural relic digital model is repaired, and a repaired cultural relic digital model is obtained, specifically including: The missing area in the initial cultural relic digital model is determined according to a missing boundary curve; A surface reconstruction is performed on the missing area by using a B-spline method, and a reconstructed area is obtained; The texture of the reconstructed area is synthesized based on a texture feature library, and the repaired cultural relic digital model is obtained.

9. The method for three-dimensional matching and digital restoration of fragments of an artifact according to claim 8, characterized in that, The texture of the reconstructed area is synthesized based on a texture feature library, and the repaired cultural relic digital model is obtained, specifically including: The following formula is used for parameterization of the reconstructed area: wherein E' represents a set of vertices of the reconstruction region, i" and j" represent horizontal and vertical coordinates of a current vertex in the reconstruction region, respectively, u i” and u j” represent parametric horizontal and vertical coordinates of the current vertex in the reconstruction region, respectively, d i”j” represents geodesic distance of the current vertex in the reconstruction region, w i”j” represents edge weight of the current vertex in the reconstruction region. The following formula is used for matching texture elements of the parameterized reconstructed area: where p' denotes the current vertex in the reconstruction region, T new denotes a texture similar to the boundary texture of the reconstruction region in the texture feature library, T pattern denotes the texture feature library, T boundary denotes the boundary texture of the reconstruction region, T denotes a texture in the texture feature library, and D() denotes a texture similarity measurement function. The following formula is used for minimizing a texture transition error, and the texture of the reconstructed area is obtained, and the repaired cultural relic digital model is obtained. wherein E texture denotes a texture transition error, Ω denotes a reconstruction region, p" denotes a current vertex in a boundary region of the reconstruction region, denotes a boundary region of the reconstruction region, T known denotes a texture of a non-missing region, μ denotes a boundary weight coefficient of the reconstruction region.

Images