High-precision scanner tree ring cross dating rapid analysis method

Abstract

The invention discloses a high-precision scanner tree ring cross dating rapid analysis method, and relates to the technical field of tree ring analysis and chronological determination. The method comprises the steps that a sample section of the tree ring is obtained through a wood slicing machine, the sample section is processed to obtain an observation face of the tree ring, the observation face is scanned to obtain an original digital image, and a high-precision scanner is adopted for scanning. According to the method, high-precision scanning imaging, polar coordinate enhancement processing and automatic annual ring boundary extraction are combined, so that full-process automatic processing of tree annual rings from image acquisition to width sequence generation is realized, manual participation links are greatly reduced, manual interpretation errors and subjective factor interference are avoided, and the tree annual ring width sequence generation efficiency is improved. The objectivity, consistency and repeatability of annual ring identification and measurement results are improved, the overall efficiency of large-batch sample treatment is improved, and the actual requirements of tree annual ring research on high throughput and rapid analysis are met.

Description

TECHNICAL FIELD

[0001] The application belongs to the technical field of tree ring analysis and chronology determination, and particularly relates to a high-precision scanner tree ring cross-dating rapid analysis method. BACKGROUND

[0002] At present, tree rings are used for reconstructing climate change, evaluating forest growth conditions and determining the age of wooden cultural relics, forming a relatively complete tree ring chronology research system. The common process includes wood sample preparation, polishing of the surface of the tree ring, observing the structure of the tree ring under a magnifying glass or a microscope, and recording the width and sequence information of the tree ring in the dating software. Some studies introduce high-resolution scanning equipment to perform area array scanning on the tree ring sample to obtain a digital image reflecting the texture and gray scale changes of the tree ring. Researchers select a measurement cut line or an area on the image, identify the tree ring boundary, and generate a tree ring width sequence or a gray scale sequence. Based on the cross-dating principle, the tree ring sequences of different samples are compared and matched to construct a longer time scale of the chronological sequence, which provides a time reference for climate analysis and cultural relic restoration applications. Related work has begun to focus on the cooperation between high-resolution scanning, tree rings and cross-dating, and has accumulated various measurement software and statistical methods.

[0003] The prior art still has some problems in realizing automatic and high-throughput analysis, such as: the key steps still rely on manual interpretation and manual operation, the number of high-resolution scanned images is large and the resolution is high, researchers need to view the images one by one, repeatedly label the tree ring boundary, manually adjust the cross-dating position of the tree ring sequence, it is difficult to complete the processing task of a large number of samples within a limited time, the manual labeling and comparison process is easily affected by personal experience and fatigue state, the stability of the cross-dating result is insufficient between different batches and different operators, and it is difficult to meet the needs of rapid analysis and batch processing in tree ring research and application.

[0004] Therefore, we provide a high-precision scanner tree ring cross-dating rapid analysis method to solve the above problems. SUMMARY

[0005] In view of the deficiencies of the prior art, the present application provides a high-precision scanner tree ring cross-dating rapid analysis method, which solves the technical problem of how to realize automatic analysis of rapid, stable and high-consistency dating of tree rings by combining high-precision scanning, image enhancement and automatic cross-matching.

[0006] To solve the above technical problems, the present application is realized by the following technical scheme.

[0007] The application discloses a high-precision scanner tree ring cross-dating fast analysis method.

[0008] S2. The original digital image is subjected to position estimation to obtain a tree ring center position, and polar coordinate enhancement processing is conducted on the tree ring center position to obtain a ring enhancement image, wherein the polar coordinate enhancement processing comprises multi-scale band-pass filtering and directional convolution operation.

[0009] S3. Ring boundary extraction is conducted on the ring enhancement image to obtain a tree ring width sequence.

[0010] S4. Standardization pretreatment is conducted on the tree ring width sequence to form standardized ring data.

[0011] S5. Cross-dating matching is conducted on the standardized ring data to form a cross-dating result.

[0012] The application further provides that the processing comprises polishing, polishing and removing surface impurities, the resolution of the high-precision scanner is 1200 dpi, and the scanning area of the high-precision scanner is 200 mm*200 mm, and the color depth is 24 bits.

[0013] The application further provides that the position estimation comprises density distribution estimation and texture direction feature estimation, the density distribution estimation is gray value extraction conducted on the original digital image, the gray value extraction is calculation of pixel density of the original digital image, and a preliminary center position is obtained based on the pixel density.

[0014] The application further provides that the texture direction feature estimation takes the preliminary center position as a reference, extracts a radial texture structure of the original digital image based on the reference, conducts symmetry evaluation on the radial texture structure to obtain a center offset vector, adjusts the preliminary center position based on the center offset vector, and the tree ring center position is obtained through the adjustment.

[0015] The application further provides that the polar coordinate enhancement processing conducts polar coordinate conversion on the original digital image to obtain a polar coordinate image, the polar coordinate conversion takes the tree ring center position as an origin, and converts the original digital image from a rectangular coordinate system into a polar coordinate system.

[0016] The application further provides that the multi-scale band-pass filtering removes low-frequency noise from the polar coordinate image to form a ring feature image, the low-frequency noise removal adopts a Gaussian filter, and the Gaussian filter comprises a smoothing kernel, a standard deviation and a filter scale parameter.

[0017] The application is further configured that the directional convolution operation performs edge detection and directional filtering on the annual ring feature image, the edge detection performs gradient extraction on the annual ring feature image, annual ring edge positions are obtained through the gradient extraction, the gradient extraction includes horizontal pixel gradient extraction and vertical pixel gradient extraction, the gradient extraction adopts a Sobel operator, and the directional filtering performs contrast enhancement on the annual ring feature image based on the annual ring edge positions, and the annual ring enhanced image is obtained based on the contrast enhancement.

[0018] The application is further configured that the annual ring boundary extraction performs global gray scale statistics on the annual ring enhanced image to obtain a gray value threshold, extracts a boundary position of the annual ring enhanced image by setting the gray value threshold, and obtains the tree ring width sequence by performing morphological erosion on the boundary position.

[0019] The application is further configured that the standardization preprocessing adopts a Z-score standardization method, and the standardized annual ring data includes standardized annual ring width data, standardized annual ring gray data and standardized annual ring density data.

[0020] The application is further configured that the cross-dating matching includes the following steps: S51. Similarity calculation is performed on the standardized annual ring data, and a similarity matrix is obtained by calculating the Euclidean distance of each pair of annual ring data points.

[0021] S52. An optimal matching path is obtained by performing path extraction on the similarity matrix.

[0022] S53. Time registration and offset extraction are performed on the optimal matching path, the time registration maps the geometric position of the optimal matching path into a sequence time corresponding relationship, the offset extraction performs time offset adjustment on the standardized annual ring data based on the sequence time corresponding relationship, and the cross-dating result is obtained through the time offset adjustment.

[0023] The application has the following beneficial effects: the application combines high-precision scanning imaging, polar coordinate enhancement processing and automatic annual ring boundary extraction, realizes full-process automatic processing of tree ring from image acquisition to width sequence generation, greatly reduces the manual participation link, avoids manual interpretation errors and subjective factor interference, improves the objectivity, consistency and repeatability of annual ring identification and measurement results, improves the overall efficiency of large quantities of sample processing, and meets the actual needs of high-throughput and rapid analysis of tree ring research.

[0024] The application realizes high-precision automatic alignment and time registration between different sample annual ring sequences by standardizing the pretreatment of the annual ring width sequence and combining the cross-dating matching method based on the similarity matrix and the optimal path extraction, improves the matching accuracy and stability of cross-dating, reduces the uncertainty caused by manual repeated adjustment, makes the dating result more reliable, and provides more stable and reliable technical support for climate change research, forest growth evaluation, wood cultural relic dating and other applications. BRIEF DESCRIPTION OF DRAWINGS

[0025] In order to more clearly illustrate the technical solutions of the embodiments of the application, the drawings needed for the embodiment description will be briefly introduced as follows.

[0026] Fig. 1 The overall method flowchart of the application.

[0027] Fig. 2 The position estimation and polar coordinate enhancement processing flowchart of the application.

[0028] Fig. 3 The cross-dating matching flowchart of the application. DETAILED DESCRIPTION

[0029] The technical solutions in the embodiments of the application will be described below with reference to the drawings in the embodiments of the application. The described embodiments are only some of the embodiments of the application, not all the embodiments.

[0030] Embodiment 1 Please refer to Figs. 1-3 The application is a high-precision scanner tree ring cross-dating rapid analysis method, which comprises the following steps: S1. Obtain the sample section of the tree ring by a wood slicing machine, process the sample section to obtain the observation surface of the tree ring, and scan the observation surface to obtain the original digital image, wherein the scanning is performed by a high-precision scanner. The processing includes polishing, polishing and removing surface impurities. The resolution of the high-precision scanner is 1200 dpi, the scanning area of the high-precision scanner is 200mmx200mm, and the color depth is 24 bits.

[0031] S2. The original digital image is used to estimate the location of the tree ring center. Polar coordinate enhancement processing is then applied to this center to obtain an enhanced tree ring image. This polar coordinate enhancement includes multi-scale bandpass filtering and directional convolution operations. Location estimation involves density distribution estimation and texture direction feature estimation. Density distribution estimation extracts grayscale values ​​from the original digital image, calculating the pixel density to obtain the initial center location. Texture direction feature estimation uses the initial center location as a reference, extracting the radial texture structure from the original digital image. Symmetry evaluation of the radial texture structure yields a center offset vector, which is used to adjust the initial center location, ultimately obtaining the tree ring center location. Polar coordinate enhancement transforms the original digital image into a polar coordinate image, using the tree ring center location as the origin. This transformation converts the original digital image from a Cartesian coordinate system to a polar coordinate system. Multi-scale bandpass filtering removes low-frequency noise from the polar coordinate image to form the tree ring feature image. This low-frequency noise removal uses a Gaussian filter, which includes a smoothing kernel, standard deviation, and filtering scale parameters. Directional convolution is used to perform edge detection and directional filtering on the annual ring feature image. Edge detection extracts gradients from the annual ring feature image to obtain the position of the annual ring edges. Gradient extraction includes horizontal pixel gradient extraction and vertical pixel gradient extraction, and the Sobel operator is used for gradient extraction. Directional filtering enhances the contrast of the annual ring feature image based on the position of the annual ring edges, and the enhanced annual ring image is obtained based on the contrast enhancement.

[0032] S3. Extract tree ring width sequences from the tree ring enhanced image. Tree ring boundary extraction involves performing global grayscale statistics on the tree ring enhanced image to obtain a grayscale threshold. The boundary positions of the tree ring enhanced image are then extracted by setting the grayscale threshold, and morphological erosion is performed on the boundary positions to obtain the tree ring width sequences.

[0033] In the enhanced image of tree rings, the gray-level histogram of the whole image is statistically analyzed. The global threshold T is automatically calculated using the Otsu method. Pixels with gray values ​​greater than T are identified as bright bands of tree rings. Two erosion operations are performed using 3×3 circular structuring elements to remove isolated noise. The radial distance between adjacent boundary points is subtracted along a fixed angle direction to obtain the width of a single tree ring, forming a sequence of tree ring widths arranged in chronological order.

[0034] S4. Standardize the tree ring width sequence to form standardized tree ring data. The standardization preprocessing adopts the Z-score standardization method. The standardized tree ring data includes standardized tree ring width data, standardized tree ring grayscale data, and standardized tree ring density data.

[0035] S5. Perform cross-dating matching on the standardized tree-ring data to generate cross-dating results. Cross-dating matching includes the following steps: S51. Similarity calculation is performed on the standardized ring data, and a similarity matrix is obtained by calculating the Euclidean distance of each pair of ring data points.

[0036] S52. The optimal matching path is obtained by path extraction on the similarity matrix.

[0037] S53. Time registration and offset extraction are performed on the optimal matching path, the time registration maps the geometric position of the optimal matching path to the sequence time correspondence, the offset extraction adjusts the time offset of the standardized ring data based on the sequence time correspondence, and the cross-dating result is obtained by time offset adjustment.

[0038] Example 2 Please refer to Fig. 2 On the basis of Example 1, the center position of the original digital image of the tree ring obtained by high-precision scanning is adaptively estimated, and the polar coordinate conversion, multi-scale band-pass filtering and directional convolution enhancement processing are combined to realize the high-contrast enhancement and reliable expression of the ring texture.

[0039] A tree ring sample image obtained by high-precision scanning is taken as input, the resolution of the original digital image is 1200 dpi, the image size is 4800x4800 pixels, and the image format is 24-bit RGB color image.

[0040] 1. Gray value extraction and pixel density calculation The original RGB color image is grayed, and the gray value of each pixel point is calculated by using the following weighted formula: The gray image I(x, y) with a gray value range of 0-255 is obtained.

[0041] Pixel density statistical analysis is performed on the gray image: a statistical window of 50x50 pixels is taken as a statistical window, the window step is set to 25 pixels, and the whole 4800x4800 image is traversed by sliding, and about 38000 statistical windows are obtained. For each window The gray mean value and the gray variance are calculated. It is found through statistics that the gray variance of multiple windows near the center region of the image is significantly higher than that of the peripheral region. The statistical results of one of the maximum variance windows are: Window center coordinates: (2412, 2389), μ = 131.4, =1458.2.

[0042] Select the window center point as the preliminary center position: .

[0043] 2. Radial texture extraction and symmetry evaluation Take the preliminary center position (2412, 2389) as the reference point to perform radial texture extraction: Angle range: 0°-360°, angle step: 1°, radius range: 1-2100 pixels.

[0044] In each angle direction θ, extract the radial gray curve .

[0045] Then calculate the correlation coefficient for the curve of the symmetric direction pair (θ, θ+180°): wherein, Gr(θ, r) represents the gray value of the rth radial pixel point in the angle direction θ, Gr(θ+180°, r) represents the gray value of the rth radial pixel point in the symmetric direction (θ+180°) opposite to θ, Gr(θ) represents the average gray value of the entire radial gray curve in the direction θ, Gr(θ+180°) represents the average gray value of the entire radial gray curve in the symmetric direction (θ+180°). ρ(θ) represents the radial texture similarity index of the direction θ and its symmetric direction (θ+180°). When ρ(θ) is close to 1, it indicates that the textures of the two directions are highly symmetric, and the center position is accurate. When ρ(θ) is close to 0, it indicates that the texture correlation is weak. When ρ(θ) is close to -1, it indicates that the texture trends are opposite, and the center is significantly offset.

[0046] 180 groups of correlation coefficients are calculated, and the results show that: The maximum correlation coefficient is 0.97, and the minimum correlation coefficient is 0.63, which appears at θ=91° / 271°.

[0047] It is shown that the annual ring symmetry is the worst in the direction of θ=91° / 271°, indicating that the center is offset.

[0048] 3. Center offset vector calculation and center correction In the direction of θ=91°, the extreme point statistics are performed on the corresponding radial gray curve, and the left and right symmetric annual ring positions are compared to obtain the average offset: X-direction offset: -6 pixels, Y-direction offset: +4 pixels.

[0049] The center offset vector is obtained: Δc=(−6,+4) Accordingly, the preliminary center is corrected: The error between the measured and artificial measured real annual ring geometric center is within ±4 pixels, about 0.085 mm.

[0050] 4. Polar coordinate conversion The polar coordinate conversion is performed on the original gray-scale image with the corrected tree ring center (2406, 2393) as the polar coordinate origin: Radial resolution: 1 pixel, angular resolution: 0.5°, radial maximum length: 2200 pixels.

[0051] The coordinate mapping is completed by bilinear interpolation, and the polar coordinate image size is 720x2200.

[0052] Among them, 720 corresponds to the number of angular direction sampling points, and 2200 corresponds to the number of radial pixel points.

[0053] 5. Multi-scale band-pass filtering The polar coordinate image is subjected to multi-scale Gaussian filtering to construct three groups of Gaussian kernels: Table 1: Multi-scale Gaussian filtering data table.

[0054] The filtering process is as follows: The low-frequency background component L(r, θ) is extracted using The high-frequency detail component H(r, θ) is extracted using .

[0055] The band-pass response is constructed: After processing, the low-frequency background change of the polar coordinate image is effectively suppressed, and the annual ring strip texture contrast is increased by about 37%, obtaining the annual ring feature image .

[0056] 6. Directional convolution operation and contrast enhancement The Sobel operator convolution is performed on the annual ring feature image . Horizontal direction operator: Vertical direction operator: The gradient amplitude is obtained: After statistics: The gradient mean , and the gradient standard deviation .

[0057] The threshold is set to: Extract G(r, θ) > 24.5 area as the candidate area of the annual ring edge.

[0058] Directional filtering and contrast enhancement are performed: The gray contrast of the light and dark bands of the annual ring is increased by about 42%, and the annual ring boundary definition is enhanced, obtaining the annual ring enhancement image .

[0059] Example 3 Please refer to Fig. 3 On the basis of Example 1 and Example 2, by performing similarity calculation, optimal matching path search and time offset extraction on the standardized annual ring sequence, automatic cross-dating and time accurate alignment between the sample annual ring sequence and the reference annual ring sequence are realized.

[0060] 1. Data preparation Two groups of annual ring width sequences that have been standardized are selected as input, wherein: The standardized sample annual ring sequence A: A={−0.42, 0.15, 0.38,−0.21, 0.64, 0.91, 0.27,−0.33,−0.58, 0.12}.

[0061] The standardized reference annual ring sequence B: B={0.36,−0.18, 0.59, 0.88, 0.31,−0.26,−0.61, 0.09, 0.41,−0.14, 0.72, 0.95}.

[0062] Among them, the length of sequence A is 10, and the length of sequence B is 12.

[0063] 2. Similarity calculation and similarity matrix construction Take the i-th data point A(i) in sequence A and the j-th data point B(j) in sequence B as a matching unit, and use Euclidean distance to construct similarity measure: For example: D(1,1)=|−0.42−0.36|=0.78.

[0064] D(1,2)=|−0.42+0.18|=0.24.

[0065] D(2,3)=|0.15−0.59|=0.44.

[0066] D(6,10)=|0.91−(−0.14)|=1.05.

[0067] The distance values of all A(i) and B(j) are calculated in the above manner to construct a 10x12 similarity matrix D. The smaller the value in the matrix, the higher the similarity of the two annual ring data points.

[0068] 3. Optimal matching path extraction The dynamic programming method is used to search the path on the similarity matrix D. The path search rules are set as follows: the starting point is D(1, 1), the ending point is D(10, 12), each step is only allowed to move to the right, down or right-down direction, the cumulative cost is the sum of all D(i, j) on the path, and the path with the minimum cumulative cost is selected as the optimal matching path.

[0069] The minimum cumulative cost path is calculated as: (1, 2) to (2, 3) to (3, 4) to (4, 5) to (5, 6) to (6, 7) to (7, 8) to (8, 9) to (9, 10) to (10, 11).

[0070] The cumulative distance of the path is 2.94, which is smaller than the cumulative distance of other candidate paths, and is determined as the optimal matching path.

[0071] 4. Time alignment and offset extraction The coordinate relationship in the optimal matching path is converted into a time corresponding relationship, for example: It can be obtained that: The average offset of each matching point is: It is shown that the sample annual ring sequence A is overall lagged by 1 year on the time axis relative to the reference sequence B.

[0072] The sample sequence A is adjusted in time offset so that A(1) corresponds to B(2). The adjusted annual ring time sequence corresponding relationship is: A(1) corresponds to the year N-9, and A(10) corresponds to the year N.

[0073] The accurate alignment of the sample annual ring sequence and the reference annual ring sequence on the time axis is completed, and the cross-dating result of the sample is obtained.

[0074] The above only describes certain exemplary embodiments of the present application in a descriptive manner, and it is needless to say that the described embodiments can be modified in various ways without departing from the spirit and scope of the present application for those skilled in the art. Therefore, the above figures and descriptions are illustrative in nature and should not be understood as limiting the scope of protection of the claims of the present application.

Claims

1. A high-precision scanner tree-ring cross-dating rapid analysis method, The method comprises the following steps: S1. obtaining a sample section of tree rings by a wood slicing machine, processing the sample section to obtain an observation surface of the tree rings, and scanning the observation surface to obtain an original digital image by using a high-precision scanner; S2. performing position estimation on the original digital image to obtain a tree ring center position, and performing polar coordinate enhancement processing on the tree ring center position to obtain a ring enhancement image, wherein the polar coordinate enhancement processing comprises multi-scale band-pass filtering and directional convolution operation; S3. performing ring boundary extraction on the ring enhancement image to obtain a tree ring width sequence; S4. performing standardization preprocessing on the tree ring width sequence to form standardized ring data; S5. performing cross-dating matching on the standardized ring data to form a cross-dating result.

2. The high precision scanner tree-ring cross-dating rapid analysis method according to claim 1, characterized in that: The processing comprises polishing, polishing and removing surface impurities, the resolution of the high-precision scanner is 1200 dpi, the scanning area of the high-precision scanner is 200 mm x 200 mm, and the color depth is 24 bits.

3. The high precision scanner tree-ring cross-dating rapid analysis method of claim 1, wherein: The position estimation comprises density distribution estimation and texture direction feature estimation, the density distribution estimation performs gray value extraction on the original digital image, the gray value extraction is to calculate the pixel density of the original digital image, and a preliminary center position is obtained based on the pixel density.

4. The high precision scanner tree-ring cross-dating rapid analysis method of claim 3, wherein: The texture direction feature estimation takes the preliminary center position as a reference, extracts a radial texture structure of the original digital image based on the reference, performs symmetry evaluation on the radial texture structure to obtain a center offset vector, adjusts the preliminary center position based on the center offset vector, and obtains the tree ring center position through the adjustment.

5. The high precision scanner tree-ring cross-dating rapid analysis method of claim 1, wherein: The polar coordinate enhancement processing performs polar coordinate conversion on the original digital image to obtain a polar coordinate image, the polar coordinate conversion takes the tree ring center position as the origin, and converts the original digital image from a rectangular coordinate system to a polar coordinate system.

6. The high precision scanner tree-ring cross-dating rapid analysis method of claim 5, wherein: The multi-scale band-pass filtering removes low-frequency noise from the polar coordinate image to form a ring feature image, the low-frequency noise removal adopts a Gaussian filter, and the Gaussian filter comprises a smoothing kernel, a standard deviation and a filter scale parameter.

7. The high precision scanner tree-ring cross-dating rapid analysis method of claim 6, wherein: The directional convolution operation performs edge detection and directional filtering on the ring feature image, the edge detection performs gradient extraction on the ring feature image to obtain a ring edge position, the gradient extraction comprises horizontal pixel gradient extraction and vertical pixel gradient extraction, the gradient extraction adopts a Sobel operator, the directional filtering performs contrast enhancement on the ring feature image based on the ring edge position, and the ring enhancement image is obtained based on the contrast enhancement.

8. The high precision scanner tree-ring cross-dating rapid analysis method of claim 1, wherein: The ring boundary extraction performs global gray scale statistics on the ring enhancement image to obtain a gray value threshold, extracts a boundary position of the ring enhancement image by setting the gray value threshold, and obtains the tree ring width sequence by performing morphological erosion on the boundary position.

9. The high precision scanner tree-ring cross-dating rapid analysis method of claim 1, wherein: The standardization pretreatment adopts a Z-score standardization method, and the standardized ring data includes standardized ring width data, standardized ring gray data and standardized ring density data.

10. The high precision scanner tree-ring cross-dating rapid analysis method of claim 1, wherein: The cross-dating matching includes the following steps: S51. Similarity calculation is performed on the standardized ring data, and a similarity matrix is obtained by calculating the Euclidean distance of each pair of ring data points; S52. An optimal matching path is obtained by performing path extraction on the similarity matrix; S53. Time registration and offset extraction are performed on the optimal matching path, the time registration maps the geometric position of the optimal matching path to a sequence time corresponding relationship, the offset extraction adjusts the time offset of the standardized ring data based on the sequence time corresponding relationship, and the cross-dating result is obtained through the time offset adjustment.

Images