Automatic detection of pith location along norway spruce timber boards on the basis of optical scanning

Knowledge of annual ring width and location of pith in relation to board cross-sections, and how these properties vary in the longitudinal direction of boards, is relevant for many purposes, such as assessment of shape mechanical properties and stability of sawn timber. Hence, the present research aims at developing a novel method and an algorithm, based on data obtained from optical surface scanning, by which the pith location along the length of sawn timber boards can be determined accurately and automatically. The first step of the method is to identify clear wood sections, free of defects along boards. Then time-frequency analysis, using the continuous wavelet transform, is applied to detect the surface annual ring width distribution of the four sides of the selected sections. Finally, the pith location is estimated by comparing annual ring width distributions on the different surfaces, and assuming that annual rings are concentric circles with the pith in the centre. The proposed algorithm was applied to a total sample of 104 Norway spruce boards. Results indicate that optical scanners and the suggested automatic method allow for accurate detection of annual ring width and location of pith along boards. For a sample of boards with the pith located within the cross-section, a mean error of 2.6 mm and 3.2 mm in the depth and thickness direction, respectively, was obtained. For a sample of boards of which 60% with pith located outside the cross-section, a mean discrepancy between automatically and manually determined pith locations of 3.9 mm and 5.8 mm in depth and thickness direction, respectively, was obtained.


Background
Different mechanical and physical properties of clear wood from softwood species can be related to the distance from the pith of the log. For Norway spruce (Picea abies), density, longitudinal modulus of elasticity (MOE), and modulus of rupture (MOR) increase significantly in the radial direction from pith to bark, whereas the longitudinal shrinkage coefficient and annual ring width decrease from pith to bark (Blouin et al. 2007;Ormarsson et al. 1999). Mechanical properties of sawn structural timber depend both on clear wood properties and on occurrence of knots, the latter being oriented in the direction from the pith and outward with larger knot diameter at larger distance from the pith (Kliger et al. 1998;Johansson 2003), meaning that relationships between different properties of sawn timber are not identical to those valid for clear wood. Nonetheless, the pith location and annual ring width, and how these change in the longitudinal direction of a board, are relevant for assessment of stiffness and strength (Hu et al. 2016) and shape stability (Ormarsson et al. 1999) of sawn timber. Knowledge of the pith location is also needed to establish detailed and accurate three-dimensional (3D) models of sawn timber, including geometry of knots and local fibre orientation on the basis of surface scanning, and attempts to develop such models have been made by some researchers, for example Hu et al. (2016) and Lukacevic et al. (2019). Furthermore, other requirements, for example regarding visual appearance of wood products, sometimes exclude the use of pieces with pith visible on surfaces. In other cases, pieces with pith enclosed within the board cross-section should not be used. Thus, it would be of practical value if commercially available scanners used for automated assessment of wood specimens could be used to accurately determine the pith location.
Over the years, some attempts have been made to automatically detect the pith location of timber boards. (Briggert et al. 2016) developed a method to reconstruct the 3D geometry of knots on the basis of data from surface laser scanning of Norway spruce timber boards. The method comprised detection of knot areas visible on the longitudinal surfaces of the board by means of tracheid effect scanning. Since branches of Norway spruce grow from the pith of the log towards the bark, Briggert et al. (2016) further utilised the detected orientation of knots to estimate the pith location along the length direction of the board. However, application of the method required that pith was located outside the cross-section of the assessed board. In addition, to be able to determine which knot surfaces (visible on different board surfaces) are part of the same knot, knowledge of an approximate pith location was needed already at the outset. This information was obtained by examination of the end cross-section at one of the board's ends. Perlin et al. (2018) proposed a method to locate the pith of a wood cross-section by utilising an ultrasonic tomography measurement technique. The method comprised mounting a fixed transmitter transducer and moving the receiver transducer around the cross-section of the specimen to record several readings of ultrasonic pulse velocities (UPVs). The basis for the proposed method was that acoustic waves travel faster in radial direction than in tangential direction. Therefore, by mapping the direction of the highest UPVs values, according to Perlin et al. (2018) the pith can be located at a position where most of these high velocity paths intersect. However, only two test specimens were used to validate the proposed method, a 25 cm diameter circular Eucalyptus grandis specimen and a 20 cm square Aplueialeiocarpa specimen. According to Perlin et al. (2018), the accuracy of the proposed method could be affected by the presence of internal defects within the timber cross-section.
In addition to the above mentioned studies, numerous studies have utilised images of cross-sections of logs generated from computer tomography (CT) X-ray scanning to predict the pith location of logs (Longuetaud et al. 2004;Bhandarkar et al. 1999;Wu and Liew 2000;Andreu and Rinnhofer 2001;Jaeger et al. 1999). Most of these studies (Longuetaud et al. 2004;Bhandarkar et al. 1999;Wu and Liew 2000;Andreu and Rinnhofer 2001) involved detection of growth rings on the cross-sectional CT images of the log slices, with the assumption that the growth rings are concentric circles centred at the pith, and application of Hough transform (HT) to the detected growth rings to estimate the pith location of the log slices. The HT-based method has been applied to both softwood (Longuetaud et al. 2004;Wu and Liew 2000;Andreu and Rinnhofer 2001) and hardwood logs (Bhandarkar et al. 1999). The work presented in Bhandarkar et al. (1999) comprised pith location of specimens obtained from four different hardwood species, namely white ash, red oak, black walnut, and hard maple. The corresponding results were validated using a total of 2176 CT images: 424 from ash, 633 from red oak, 456 from black walnut, and 663 from hard maple. According to Bhandarkar et al. (1999), more than 96% of the determined pith locations were within a 15 mm margin of error and more than 77% of the results were within 7.5 mm margin of error. In addition, Longuetaud et al. (2004) used the HT-based method to predict the pith location of softwood logs taken from 8 different trees. A total of 12,839 CT cross-section images of the logs were used. The results from Longuetaud et al. (2004) showed a 95th percentile of 1.97 mm and 99th percentile of 3.62 mm accuracy. In addition to the HT-based method, Jaeger et al. (1999) proposed a method to locate the pith position on the basis of CT images of Norway spruce logs by using the cross-sectional density map. The work was based on the fact that wood density closer to the pith is lower than what it is farther away from the pith. However, in the work presented by Longuetaud (2005) the accuracy of the HTbased CT method was compared against the density based CT method and it was concluded that the HT-based method is more accurate than the method based on density variation over the cross-section.
In general, results obtained on the basis of CT X-ray scanning (Longuetaud et al. 2004;Bhandarkar et al. 1999;Wu and Liew 2000;Andreu and Rinnhofer 2001;Jaeger et al. 1999) are more accurate than those obtained using the other proposed methods (Briggert et al. 2016;Perlin et al. 2018). However, application of the CT based methods has this far been limited to pith location of logs, where the pith is certainly located within the cross-section. Additionally, the high cost of running a CT scanner makes its application limited. On the other hand, the methods presented by Briggert et al. (2016) and Perlin et al. (2018) were applied only to very few samples, and the pith locations determined in Briggert et al. (2016) were not validated against actual pith locations. Hence, the present research aims at developing a method and an algorithm, based on data obtained from optical surface scanning, by which pith location along the length of sawn timber boards can be determined accurately and automatically.

Purpose, objectives and limitations
The purpose of the present study is to examine the possibility of developing an accurate and robust method and algorithm, solely based on information obtained from optical scanning of longitudinal surfaces, to estimate the pith location of Norway spruce timber boards. The work focuses on two objectives: 1. To detect the annual ring width distribution, which is equivalent to the local annual ring width across each of the four longitudinal timber board surfaces at selected positions along the board, and 2. To determine the pith location, and the average radial distance between annual rings, on the basis of the detected annual ring width distributions across surfaces at the selected positions along the timber board.
The scope of the current study is limited to applications on knot-free clear wood cross-sections of dried and planed Norway spruce timber boards. The fact that drying causes some distortion of the board cross-section is, however, not considered. In this work, annual rings at clear wood sections of boards, assumed to be perfectly rectangular in shape, are assumed to be concentric circles, the pith located at the centre, with a constant radial distance between the circles. However, the pith location for cross-sections that include knots can be estimated using linear interpolation between the pith locations determined at adjacent clear wood sections.

Material and data obtained from scanning
In this study, a total sample of 104 Norway spruce timber boards with nominal dimensions 45 × 145 × 4500 mm 3 was analysed. Regarding the position of the pith, the sample contained both boards with the pith located outside and boards with the pith located inside their cross-sections. All boards included in this study had already been dried to 12% MC and examined, manually and by means of an optical scanner, within a previous research project reported by Olsson and Oscarsson (2017), which facilitated the present study.
Basic data, used to detect annual rings on board surfaces, was obtained using an optical industry wood scanner equipped with LED lights, colour cameras, multi-sensor cameras, and line and dot lasers. In the scanning procedure, the timber board was fed in length-wise direction with constant speed through the scanner/measurement unit by aid of external conveyor belts attached to the scanner. Data delivered by the scanner consists of red, green and blue (RGB) channel images and data of local in-plane fibre direction of all the four sides of the scanned timber board. An approximate pixel size in the RGB images is 0.8 × 0.07 mm 2 (lengthwise × crosswise resolution), and the resolution for the local in-plane fibre direction data is approximately 1 × 4.4 mm 2 (lengthwise × crosswise). The in-plane fibre directions were determined by utilising the so-called tracheid effect, which means that when light from a concentrated source illuminates a wood surface, it propagates more in the direction parallel to the fibres (tracheids) than in perpendicular direction (Soest et al. 1993;Briggert et al. 2018). The lengthwise resolutions of both RGB images and in-plane fibre direction data depend on the speed of the board fed through the scanner, the lower the speed the higher the resolution, whereas the crosswise resolution is independent of speed.

Automatic procedure for calculation of pith location
The procedure and method employed to detect annual ring width distribution and determine the pith location along timber boards consists of three automatic steps.
Step 1: Identify knot-free clear wood sections along boards on the basis of knowledge of local fibre orientation obtained from tracheid effect scanning.
Step 2: Detect local annual ring width on all four sides on the basis of grayscale images of the wood surfaces. This is done for clear wood sections selected among the identified, in Step 1, clear wood sections.
Step 3: Determine the pith location of the selected clear wood sections using a search algorithm by which an error norm is minimised. Here it is assumed that annual rings are concentric circles, centred at the pith, and that the radial distance between adjacent annual rings is constant.
Each of the three steps is presented in detail in separate sub-sections below. Calculations are performed using the software Matlab including the wavelet toolbox (MATLAB 2018).

Identification of clear wood sections
Since the current work is limited to determining the pith location of clear wood sections of boards, the first step of the algorithm is to identify clear wood sections along boards. For this purpose, data of in-plane fibre directions obtained from the surface laser scanning is utilised. Herein, a clear wood section is defined as the centre of a 10 mm long segment, in longitudinal direction, across the four sides within which a maximum of 10% of all the determined in-plane fibre directions have an angle that exceeds 12 o with respect to the longitudinal direction of the board. Figure. 1a shows colour images of the four surfaces of a board and Fig. 1b shows local fibre angles on the surfaces of the same board represented by coloured fields. The lines drawn across the images in Fig. 1a mark the 10 mm long segments in which more than 10 % of the fibre angles exceed 12 o . Consequently, the centre of any 10 mm segment in areas that are not marked by lines drawn across the board images in Fig. 1a represents, according to the definition above, a clear wood section.

Detection of local annual ring width distribution on the board surfaces
Once the clear wood sections along the timber board are identified, the next step is to detect, on these sections, the surface annual ring widths or, more precisely, the distances between points of intersection of annual rings and board surfaces. For this purpose, a grayscale image, which is a single channel image among the three RGB channels obtained from the optical scanning, is used. In the current work, the green channel is used to produce the grayscale images, but similar results can also be achieved by using the other two channels. Each pixel in the grayscale image has a light intensity value between 0 (black) and 1 (white).
To filter out random noise of the raw grayscale images obtained from the scanner, an adaptive pixel wise low pass Wiener filtering (Lim 1990) is applied and the resulting filtered grayscale image is utilised further. In the applied filter, the light intensity of each pixel is adjusted based on the mean and standard deviation of the light intensity of the neighbouring pixels. In the current work, a 3 × 5 (lengthwise × crosswise) pixel neighbourhood is adopted.  Figure 2c shows the corresponding light intensity variation across the same section as indicated in Fig. 2b, and Fig. 2d shows calculated local wavelengths which corresponds to annual ring width distribution over the same part of the section as represented in Fig. 2b-c.
As can be seen by comparing the grayscale image shown in Fig. 2b with the light intensity plot shown in Fig. 2c, the latewood areas, visible on the surface of the board, correspond to local minima of the light intensity plot. Therefore, calculating the distance between the positions of local minima over the board width means calculating the width between annual rings, locally, over the section. By considering the light intensity variation over the section as a time-domain signal, the distance between consecutive annual rings from the grayscale image corresponds to wavelengths of the time-domain signal. In this context, 'time' actually represents the position in the For analysing the frequency content of a time-domain signal fast Fourier transform (FFT) is often used. However, FFT only gives the frequency content of the signal with no information on where, in time domain, a specific frequency occurs. Consequently, FFT is not very suitable for the current application where knowledge of localisation of different frequencies in time domain is needed. Thus, to transform the time domain signal with both time and frequency localisation, the continuous wavelet transform (CWT) (Lilly and Olhede 2012) was applied to obtain results such as those displayed in Fig. 2d. As in FFT, CWT measures the similarity between a signal and an analysing function. In CWT, the analysing function is not a sinusoidal function, which is used for FFT, but a waveform function called wavelet that has a limited duration in time and a mean value of zero. The CWT of a time domain signal x(t), which is the light intensity as a function of the position across a side of the board, can be expressed (Lilly 2017) as where s is a scaling parameter and a translation parameter of the wavelet, . The parameter s represents frequency scaling of the wavelet, giving frequency localisation of the signal, and the parameter represents the shift in time of the wavelet, giving time localisation.
Several different families of wavelets have been suggested, as for example Morlet wavelet, Bump wavelet and generalised Morse wavelet (Lilly and Olhede 2012;Lilly 2017;Jiang and Suter 2017). Each wavelet family has different features and they are suitable for different purposes. For instance, the Morlet wavelet has equal variance in time and frequency, giving similar time and frequency localisation of a signal (Lilly and Olhede 2012;Jiang and Suter 2017). The Bump wavelet, on the other hand, has a narrow variance in frequency and wider variance in time giving better frequency localisation than that of a Morlet wavelet (Jiang and Suter 2017). The generalised Morse wavelet can be designed to have frequency and time localisation capabilities of different wavelet families, such as a narrower variance in frequency similar to the Bump wavelet or equal variance in time and frequency like the Morlet wavelet, A lower value, on the other hand, leads to a better time localisation. For a more detailed explanation of the generalised Morse wavelet the reader is advised to study Lilly and Olhede (2012) and Lilly (2017).
As can be observed from Fig. 2c-d, the considered time-domain signal or the light intensity variation of a single clear wood section may have a considerable frequency and amplitude variability across the board width. Since the generalised Morse wavelet is an exact analytic wavelet, it can handle signals with time-varying frequency and amplitude variability better than the approximately analytic wavelets such as Morlet wavelet (Lilly and Olhede 2008). According to Lilly and Olhede (2008) analytic wavelets are complex-valued wavelets which have zero power at negative frequencies. Hence, due to the nature of the time-domain signal at hand, the generalised Morse wavelet is chosen as a suitable wavelet function for detecting the local annual ring distribution on board surfaces. Herein, parameter values set to = 3 and = 27 are chosen in order to get a minimised Heisenberg area, which can be interpreted as a minimised product of the standard deviations of the wavelet function in time and frequency domain, respectively (Lilly and Olhede 2012; Lilly 2017). Figure 3 shows an example where the CWT method, using generalised Morse wavelet with the specified and values, is applied to a clear wood section indicated by the blue dashed line drawn at the longitudinal position 490 mm in Fig. 3a. Figure 3b shows the corresponding light intensity signal of the selected clear wood section. Figure 3c illustrates the time and frequency localised power distribution obtained after the transformation. The (3) a , = 2 e ∕ power distribution is the value of the wavelet transform, W ( , s) , calculated using Eq. 1. Figure 3d shows the detected wavelength distribution (representing surface annual ring widths) calculated, over the board width, as the inverse of the frequencies indicated by highest power values, see Fig. 3c, over the board surface.

Estimation of pith location and average annual growth rate
When the surface annual ring width distributions visible on the four sides of the selected clear wood section have been detected, the third and final step of the algorithm concerning determination of the pith location and radial distance between annual rings follows. The annual rings visible on the surface of a timber board show where annual rings intersect with the four sides of the board. Thus, the annual ring width distribution visible on the surfaces of the board can be related to the pith location and the radial growth rate of the tree. For known surface annual ring width distributions, as determined in Sect. 3.1.2, an iterative algorithm is employed to calculate coordinates for the pith location (x, y) and the average radial distance between annual rings (r), to determine (x, y, r) such that a defined error norm is minimised. Figure 4, designed to facilitate the description of this step, shows in Fig. 4a, a coordinate system and a photograph of a knot-free board cross-section on which the true pith location (x p , y p ) is indicated with a blue cross mark, and the four sides of the cross-section are marked with blue lines. Figure 4b, shows a coordinate system and a rectangular area drawn with red lines, representing an artificial board cross-section on which equidistant concentric circles, with a distance r 1 and a centre point (x 1 , y 1 ) indicated by a red cross mark, are drawn. Figure 4d shows a rectangular area/board crosssection where (x p , y p ) and (x 1 , y 1 ) are marked by blue and red cross marks, respectively. Figure 4c shows calculated wavelengths over the top surface of the cross-section shown in Fig. 4a-b, drawn with blue and red lines, respectively, and Fig. 4e shows, correspondingly, the calculated wavelengths of the bottom surfaces.
In the first iteration of a procedure to determine the pith location, an arbitrary pith location and an average annual ring width (x 1 , y 1 , r 1 ) are assumed and used to generate an artificial cross-section with growth rings, see Fig. 4b. Then, the distances between consecutive intersection points of the artificially generated growth rings and the four sides of the board cross-section are calculated. The red lines in Fig. 4c and e show such calculated artificial annual ring width distribution plotted for the top and bottom side of the crosssection, respectively. However, as can be seen in Fig. 4c, and Fig. 4e, there is a considerable discrepancy between the actually detected annual ring width distribution, the blue lines, and the artificially calculated distribution, the red lines. This means that the initially assumed parameters, which are the pith location and average annual ring width (x 1 , y 1 , r 1 ) , are very different from the actual true pith location (x p , y p ) and actual average annual ring width (r p ) of the cross-section. The discrepancy between the actually detected and the artificially generated annual ring width distribution is quantified using a cost function, where a high cost (error norm) indicates a large discrepancy between the real and artificial annual ring width distribution and a low cost indicates a small discrepancy. The defined cost function, giving the cost valid for the j th iteration, is calculated as where N is the number of local wavelength data points gathered from the continuous time-frequency results of the CWT, see Fig. 3c, with a constant resolution from the four surfaces of the evaluated clear wood section, d ik is the local wavelength, detected using CWT, in position i (out of the N data points) on surface k, and d ij is the corresponding wavelength calculated on the basis of the artificial annual Fig. 3 a Clear wood section indicated on a grayscale image, b light intensity variation of the section, c the time and frequency localised power distribution obtained after the CWT operation, and d the resulting localised wavelength/annual ring width distribution over the board width rings, determined by the parameter values (x j , y j , r j ) , intersecting with the surface at the same position. When the cost function is calculated for the j th iteration, an optimisation algorithm is employed to identify a position (x j+1 , y j+1 , r j+1 ) giving a lower cost, C j+1 . Eventually, after several iterations, the global minimum point of the cost function, which corresponds to the determined pith location and average annual ring width, within a predefined search domain in the x-and y-direction is identified. Due to its robustness and relatively fast convergence, a simplex optimisation algorithm as presented by Lagarias et al. (1998) is used for the minimisation of the defined cost function. Figure 5b shows a three-dimensional plot of the cost function calculated for a grid of points (the grid employed is indicated on the photograph of the cross-section shown in Fig. 5b) over a predefined search domain in x-and y-direction, for a clear wood section. As can be seen, the calculated minimum point of the cost function coincides closely with the actual pith location of the cross-section which is indicated on the photograph of the cross-section shown in Fig. 5a.

Manual determination of pith location
The estimated pith locations obtained from the automatic procedure described in Sect. 3.1 are validated against manually determined pith locations. The total sample of 104 boards, see Sect. 2, was divided into two samples. Sample one consisted of four boards, in all of which the pith was located within the cross-section. Sample two consisted of 100 boards, of which some with the pith located inside and some with the pith located outside the cross-section. Different methods for manual determination of pith location were used for the two samples. For sample one, pith locations were obtained by cutting the board transversally such that the same clear wood sections as those evaluated in the automatic procedure, see description in Sect. 3.1.1, were evaluated manually by measuring, on the cross-section, the distance from the sides of each board to obtain x-and y-coordinates of the pith with respect to a pre-defined coordinate system. Altogether, 45 clear wood sections were evaluated for the four boards of sample one. Figure 6b shows the predefined coordinate system, an image of an evaluated clear wood section and an example of a manually determined pith location, with x-and y-coordinates indicated on the rulers infolded on two sides of the cross-section. Fig. 4 a A 45 × 145 mm 2 clear wood cross-section with pith location (x p , y p ) , b artificially generated growth rings for an assumed pith location and assumed annual ring width (x 1 , y 1 , r 1 ) , c, e for top and bottom surface, respectively, real detected annual ring width distribution corresponding to a true pith location (x p , y p ) , which is represented by blue lines, and the calculated annual ring width distribution for the artificially generated growth rings corresponding to the parameters (x 1 , y 1 , r 1 ) , represented by red lines, and d a 45 × 145 mm 2 rectangular area with positions (x p , y p ) and (x 1 , y 1 ) marked For sample two, a method presented by Hu et al. (2018) was employed. In this method, the pith location of a clear wood cross-section is determined by manually fitting arcs corresponding to concentric circles, which are printed on a transparent plastic sheet, to the annual rings visible on the end cross-sections of the boards. Figure 6a shows the transparent plastic sheet with several arcs printed on it, and Fig. 6c presents an example where the described method is used to estimate the pith location of a cross-section.
The work on sample two was performed by two technicians, hired for the task in connection with the investigation presented by Olsson and Oscarsson (2017). For this sample, the pith location was manually determined only at the two end sections, i.e., at the root and top end of each of the 100 boards included in the sample, giving in total, 200 evaluated clear wood sections. It is reasonable to assume that the accuracy obtained using this method is at its highest if the pith is located inside the cross-section. In cases where the pith is located far outside the cross-section, and if the annual rings do not coincide perfectly with concentric circles, the accuracy of the manually determined pith locations might be lower.

Result and discussion
As described in Sects. 2 and 3.2, a total of 104 Norway spruce timber boards, divided into two samples, were analysed. Estimation of pith location using the automatic procedure described in Sect. 3 was done for about 10 clear wood sections per board, giving a total of 1045 clear wood sections. For boards in sample one, manual determination of pith locations was done at the very same cross-sections as those evaluated using the automatic procedure. However, for boards in sample two, pith locations were determined manually only at end cross-sections of the boards, giving a total of 200 sections where manually and automatically determined pith locations could be compared. Therefore, to make the best comparison possible, additional clear wood sections located as closely as possible to the board end section were selected for automatic determination of pith location. Thus, the average distance between automatically evaluated sections and board end sections was about 5 mm. In conclusion, the results presented below consist of comparisons of manually and automatically determined pith locations of 45 sections for sample one, and of 200 sections for sample two.
In Fig. 7, the distances between automatically and manually detected pith locations, in x-and y-direction, respectively, of a total of 45 sections of sample one are displayed. Different colours are used for each of the four boards to represent evaluated cross-sections. In Fig. 8, images of two adjacent sides from the scanner and photographs of two selected cross-sections of board no. 1 are shown. Sections along the board where pith locations are estimated are indicated by black cross marks. Red and blue lines, drawn on top of the images from scanning, represent linear interpolations between the automatically and manually determined pith locations, respectively. The photographs of the crosssections (marked Sec 1.1 and Sec 2.2) show sections, out of the evaluated sections for board no. 1, with the largest and smallest discrepancy between automatically and manually determined locations of pith. The largest discrepancy was 8 mm and 8 mm, respectively, in x-and y-direction. The smallest discrepancy was 0 mm in each direction, representing a perfect match in both x-and y-direction.
In Fig. 9, the manually determined pith locations of sample two are displayed. Pith locations of the 100 root end cross-sections are indicated by red dots, and pith locations of the 100 top end cross-sections are indicated by blue dots. Out of the altogether 200 pith locations, around 60% lie outside the cross-section. More than 95% of the determined pith locations were located within ± 15 mm distance from the top surfaces of the boards, which implies that the boards were sawn from logs using a 2X-log cutting pattern. A closer look at the distribution of blue dot marks (Fig. 9) indicates that the technician who manually determined pith location in the top end of the boards worked, in y-direction, with a precision of 5 mm increments. This in turn gives an indication of the precision obtained for manually determined pith location for sample two.
In Table 1, results of a statistical analysis of the same data as displayed in Fig. 7 are presented. This includes mean values, standard deviations, medians, and 80 th , 85 th , 90 th and 95 th percentiles of discrepancies between automatically and manually determined pith locations in x-and y-direction, respectively. Since all the evaluated cross-sections of sample one contained pith, and manual determination of pith location was done by direct measurements of distances from board edges to pith, the discrepancies presented in Table 1 should, for the most part, represent errors of the automatically detected pith locations. Even so, the errors are small, with mean errors of 2.6 mm and 3.2 mm, respectively, in x-and y-direction, and with errors smaller than 6.3 mm and 7.5 mm, respectively, in the two directions for 90% of the evaluated board cross-sections.
The histograms presented in Fig. 10 show the discrepancy between manually and automatically determined pith locations, in x-and y-direction, respectively, of sample two. In Fig. 6 a Plastic sheet used to manually determine the pith location for end cross-sections of boards of sample two, and b manual determination of pith location, applied to sections of boards of sample one, c example of manually determined true pith location Table 2, results from a statistical analysis , which is similar to those presented for sample one in Table 1, of the pith location results of sample two are presented. Not surprisingly, the discrepancies are larger for sample two than for sample one, with mean discrepancies of 3.9 and 5.8 mm, respectively, in x-and y-directions, and discrepancies of up to 9.1 mm and 13.5 mm, respectively, in the two directions for 90 percent of the pith locations. However, for sample two the manually determined pith locations are more uncertain than those for the first sample, and thus the discrepancy between automatically and manually detected pith locations should not be interpreted as errors of the automatic procedure, but rather as upper limits for such errors. The discrepancies in the y-direction are larger than those in the x-direction. This is partly explained by the greater difficulty when manually determining the pith location in the y-direction (in x-direction, the technician could partly rely on estimated symmetry of the annual rings seen on the cross-section). However, it is also likely that the automatic determination of pith is more accurate in x-directions than in y-direction. The width of the investigated timber boards was 145 mm in x-direction and only 45 mm in y-direction, and the pith was always located 'within' the board cross-section with respect to the x-direction, but not with respect to the y-direction. Thus, there is for the x-direction (1) more relevant annual ring width data available, that means more data is available on the top and bottom surface of the board compared to the left and right side, and (2) annual ring width data available on both sides of the actual pith location. In y-direction, on the other hand, small errors of determined annual ring width distributions on surfaces, obtained as described in Sect. 3.1.2, or the fact that real annual rings are not actually concentric circles with equal distance in between, as assumed herein, may lead to significant errors of the estimated pith locations. Of course, when annual rings do not resemble concentric circles and pith is located outside the cross-section, some Regarding computational time, to estimate the pith location along a single board, evaluating about ten clear wood sections along the board, the automatic procedure described in Sect. 3.1 takes a total of around 1.8 s on a common workstation, using computer code in the software Matlab with a mixture of Matlab built-in functions and functions defined by the authors. Out of the total computation time, around 3.5% was spent on the time frequency analysis, see Sect. 3.1.2, and the remaining 96.5% was spent on the iterative calculations to minimise the cost function presented in Sect. 3.1.3. However, in the present research the effort put to minimise the calculation time was limited. Thus, with some more effort and/or implementation in other software, the calculation time of the proposed algorithm would match a typical industrial scanning speed of around 1 second per board.

Conclusion
The aim of the present research was (1) to detect the surface annual ring width distribution and (2) to estimate the pith location for automatically selected clear wood sections along timber boards based on the detected annual ring width distribution. Detection of the surface annual ring width distribution was done by applying a time frequency analysis on the light intensity variation of the four sides of the section. A continuous wavelet transform (CWT) with the generalised Morse wavelet was used for the time frequency analysis. Following the time frequency analysis, a new automatic method to estimate the pith location of timber boards was presented. By application of this method, which is solely based on the information obtained from optical scanning of the four sides of timber boards, it was possible to estimate the pith location accurately even in cases where the pith was located outside the cross-section of the board. For a sample of boards with the pith located within the cross-section, mean errors as small as 2.6 mm and 3.2 mm in the x-and y-direction, respectively, were obtained. For a sample of boards of which 60% had pith located outside the cross-section, the proposed method still delivered a mean discrepancy between the automatically and manually determined pith locations of only 3.9 mm and 5.8 mm in the x-and y-direction, respectively.
The scope of the current study was limited to application to knot-free clear wood cross-sections of planed Norway spruce timber boards. Annual growth rings at clear wood sections were assumed to be concentric circles, with the pith located at the centre, with a constant radial distance between them. For cross-sections that include knots, position of pith would be estimated using linear interpolation between the pith locations determined at adjacent clear wood sections. Further work should focus on eliminating the assumption of constant radial distance between annual rings, on improved accuracy and on reduced calculation time to fully meet industry requirements on speed. Moreover, the potential of utilising knowledge of pith location for improved assessment and grading of wood, as for example for improved accuracy of machine strength grading of sawn timber, should be evaluated.