Rémi Decelle
email: remi.decelle@loria.fr

Phuc Ngo

Isabelle Debled-Rennesson

Frédéric Mothe

Fleur Longuetaud

Pith Estimation on Tree Log End Images

Keywords: Agent-based method, Local orientation, Hough transform

In this paper, we present an algorithm for pith estimation from digital images of wood cross-sections. The method is based on a probabilistic approach, namely ant colony optimization (ACO). After introducing the approach, we describe the implementation and the reproduction of the method linking to an online demonstration. Results show that the approach performs as well as state-of-the-art methods. The estimated pith is below 5mm from the ground truth. It is a fast method that could be used in real-time environment. This paper also gives the details about the intern parameter choice and shows how to use the C++ source code for testing, as well as provides limit cases of the proposed method and future improvements.

Introduction

The centre of the annual rings, also called pith, is one of the most important feature to be detected since it can be related to wood quality [START_REF] Akachuku | Variations in pith eccentricity and ring width in teak (tectona grandis lf)[END_REF][START_REF] Fallah | Investigating the out-ofroundness and pith-off-centre in stems of three broadleaved species in hyrcanian forests[END_REF][START_REF] Rune | Basal sweep and compression wood in young scots pine trees[END_REF] and it allows to extract other features on log-end image [START_REF] Entacher | Towards an automated generation of tree ring profiles from ct-images[END_REF][START_REF] Fabijanska | Towards automatic tree rings detection in images of scanned wood samples[END_REF][START_REF] Hanning | Determining the average annual ring width on the front side of lumber[END_REF][START_REF] Nordmark | Models of knots and log geometry of young pinus sylvestris sawlogs extracted from computed tomographic images[END_REF] such as annual rings, ring widths, knots, heartwood and sapwood. In the literature, several methods have been proposed for pith detection on log cross-sections. Most of them [START_REF] Bhandarkar | A system for detection of internal log defects by computer analysis of axial ct images[END_REF][START_REF] Boukadida | Pithextract: A robust algorithm for pith detection in computer tomography images of wood-application to 125 logs from 17 tree species[END_REF][START_REF] Gazo | A fast pith detection for computed tomography scanned hardwood logs[END_REF][START_REF] Longuetaud | Automatic knot detection and measurements from x-ray ct images of wood: A review and validation of an improved algorithm on softwood samples[END_REF][START_REF] Longuetaud | Automatic detection of pith on ct images of spruce logs[END_REF][START_REF] Wei | On the use of x-ray computed tomography for determining wood properties: a review[END_REF] have been developed for X-ray computed tomographic (CT) images. The techniques based on CT images allow an efficient and robust detection of external and internal characteristics of tree logs, including the pith. However, the CT scanners are very expensive, and not every laboratory or wood-industry sites (e.g., sawmills) can acquire such a device.

Recently, there have been some efforts to develop pith detection methods on RGB images of cross sections from tree logs [START_REF] Kurdthongmee | A comparative study of the effectiveness of using popular dnn object detection algorithms for pith detection in cross-sectional images of parawood[END_REF][START_REF] Kurdthongmee | A fast algorithm to approximate the pith location of rubberwood timber from a normal camera image[END_REF][START_REF] Norell | Estimation of pith position in untreated log ends in sawmill environments[END_REF][START_REF] Schraml | Pith estimation on rough log end images using local fourier spectrum analysis[END_REF]. Contrary to CT images, RGB images exhibit disturbances like sawing marks, dirt or ambient light variations, which make the detection more challenging (see Fig. 1). On the other side, the acquisition of such images can be done with low-cost and more accessible devices (e.g., smartphone camera, industrial camera, . . . ) and could be used everywhere, from the forest on the harvester, to the sawmill stocking area, or at the road side for wood sells. Furthermore, the current camera technologies provide images of quality and high resolution, and this allows us to extract the wood quality features, including pith, on such images using the image processing methods.

To the best of our knowledge, four works [START_REF] Kurdthongmee | A comparative study of the effectiveness of using popular dnn object detection algorithms for pith detection in cross-sectional images of parawood[END_REF][START_REF] Kurdthongmee | A fast algorithm to approximate the pith location of rubberwood timber from a normal camera image[END_REF][START_REF] Norell | Estimation of pith position in untreated log ends in sawmill environments[END_REF][START_REF] Schraml | Pith estimation on rough log end images using local fourier spectrum analysis[END_REF] have been published in the context of pith detection on digital images of rough, untreated log ends. Except in [START_REF] Kurdthongmee | A comparative study of the effectiveness of using popular dnn object detection algorithms for pith detection in cross-sectional images of parawood[END_REF] which uses a deep neural network (DNN), the others rely on tree ring analysis and use image processing tools for the detection. More precisely, it is assumed that most tree rings close to the pith approximate a circular shape, and thus the normal directions of these rings would point towards the pith. Based on this idea, the pith detection is generally processed in three steps:

1. Estimate normal directions from tree ring local orientations. 2. Accumulate the normals in an accumulation space. [START_REF] Boukadida | Pithextract: A robust algorithm for pith detection in computer tomography images of wood-application to 125 logs from 17 tree species[END_REF]. Extract the pith from the accumulation space.

Generally, Hough transform [START_REF] Duda | Use of the hough transformation to detect lines and curves in pictures[END_REF] is used as accumulator of normal directions, then the pith is identified at the point with maximum accumulation, or the barycenter of points above a given accumulation threshold. In other words, the pith detection methods differ at the local orientation estimation step. In [START_REF] Norell | Estimation of pith position in untreated log ends in sawmill environments[END_REF], Norell and Borgefors presented two detection methods using two different techniques to estimate the normal directions: the quadrature filters and the Laplacian pyramids. The proposed methods are robust to disturbances; e.g., rot, dirt or snow. However, a prior segmentation of the log end is needed before the detection. Later, Schraml and Uhl [START_REF] Schraml | Pith estimation on rough log end images using local fourier spectrum analysis[END_REF] proposed to compute the local orientations with Fourier spectrum analysis. The approach was fast, robust, and accurate in estimating the pith position. It, however, requires some preconditions about the cross-section size and its location in the image for initializing the computation. Recently, Kurdthongmee et al. [START_REF] Kurdthongmee | A fast algorithm to approximate the pith location of rubberwood timber from a normal camera image[END_REF] used histogram of oriented gradient (HOG) to estimate normals of the tree rings. As stated in the paper, the algorithm provides only an approximation of the pith location, and needs more treatments to identify it exactly.

In this paper, we propose a general method to estimate pith location on digital images taken in realistic environments; e.g., in the sawmill, forest, logyard, or on the road. The raw images are directly processed without any prior segmentation nor knowledge of log end visual appearance, shape or location. In particular, the proposed method must not only provide accurate pith estimation, but also be efficient in computation time to be used in real-time applications. To this end, we consider a smooth gradient based method to compute local orientations; this method was originally used for fingerprint images [START_REF] Hong | Fingerprint image enhancement: algorithm and performance evaluation[END_REF]. Then, a probabilistic approach, based on Ant Colony Optimization (ACO) [START_REF] Dorigo | Ant system: optimization by a colony of cooperating agents[END_REF], is performed to accumulate the normals of tree rings in a robust way. Finally, the pith is located at the barycenter of points above a given accumulation threshold.

The proposed method is described in the following section (Section 2), with the details of the local orientation computation and the ACO algorithm. Section 3 gives the description of the source code and its usage. The experimental results and parameter discussions are addressed in Section 4, followed by the conclusion.

Algorithm for pith detection

The proposed algorithm to estimate pith location on digital images is composed of four steps. Firstly, a pre-processing is applied on input image to remove sawing marks visible on log ends. In case of high-resolution images, a resizing step can be applied to reduce the computation time. Secondly, we compute local orientations for pixels of the pre-processed image. Then, the ACO algorithm is used to accumulate the normals, and finally extract the pith from this accumulator. For an accurate pith estimation, the ACO algorithm and the pith extraction are performed twice: the first to coarsely estimate the pith region, and the second for the precise pith location.

Pre-processing image

Pith estimation methods based on ring analysis strongly depend on local orientation estimations. In this paper, we work with raw images in which there might be sawing marks on rough log ends. The presence of sawing marks perturbs the orientation estimations. To reduce errors induced by sawing marks and also computation time, we perform this pre-processing step. Firstly, the input image is converted into gray-scale, then down-sampled with a factor s using bi-linear interpolation. Secondly, we remove sawing marks using the method proposed in [START_REF] Norell | Automatic counting of annual rings on Pinus sylvestris end faces in sawmill industry[END_REF] which is based on Fast Fourier Transform (FFT).

Typically, sawing marks are straight lines being parallel or in fan-shape and not always evenly spaced. This repetitive pattern suggests that filtering in the Fourier domain is suitable to reduce them. Indeed, in the Fourier spectrum, they correspond to the line passing through the centre with a direction perpendicular to them. In other words, the energy level will be high along this line. Therefore, by reducing this energy and converting filtered spectrum back to spatial domain, sawing marks can be removed or at least reduced. More precisely, we first compute a Fast Fourier Transform (FFT) and filter the Fourier spectrum with a band-pass filter, also remove both horizontal and vertical lines. Then, we threshold it with a value λ to filter high value energy points corresponding to sawing marks. A line fitting, using principal component analysis (PCA), is applied on the obtained points to find the direction of the line. This line is further convoluted with a Gaussian filter of standard deviation σ, and pixel-wise multiplied with the original Fourier spectrum to reduce energy along the line. The result is transformed back into a spatial image using the inverse FFT. The process is illustrated in Fig. 2 (a-e).

Local orientation

After removing sawing marks, we now compute normal directions for blocks of pixels in the image using a smooth gradient based method. The method was used in [START_REF] Hong | Fingerprint image enhancement: algorithm and performance evaluation[END_REF] to assess the local orientation in fingerprint images. It is a least mean square orientation estimation in a local area, namely a window of size w. More precisely, the gradient

∇(u, v) = [∇ x (u, v), ∇ y (u, v)] t is estimated for each block centered at pixel (u, v) as follows. ∇ x (u, v) = i=v+ w 2 i=u-w 2 j=v+ w 2 j=u-w 2 2δ x (i, j)δ y (i, j) (1)
∇ y (u, v) = i=v+ w 2 i=u-w 2 j=v+ w 2 j=u-w 2 δ 2 x (i, j)δ 2 y (i, j) (2)
where δ x (i, j) and δ y (i, j) are the derivatives with respect to x and y of the pixel (i, j). The derivatives are estimated by a Sobel operator [START_REF] Sobel | An isotropic 3x3 image gradient operator[END_REF]. Then, the local orientation of the block centered at (u, v) is computed as:

θ(u, v) = 1 2 tan -1 ∇ y (u, v) ∇ x (u, v) (3)
Fig. 2 (f) shows an example of local orientations estimated by this method.

Ant colony optimization

We now describe the process of accumulating normals using ACO which is an algorithm inspired by the behavior of ant species. Ants deposit pheromones that help other ants of the colony to make the best choice in their goal. Ant system [START_REF] Dorigo | Ant system: optimization by a colony of cooperating agents[END_REF] was the first ACO algorithm, it is applied for solving different combinatorial optimization problems; e.g., traveling salesman problem (TSP), quadratic assignment problem (QAP) and the job-shop scheduling problem (JSP). Since then, a large number of ACO algorithms have been developed to address various problems like edge detection [START_REF] Nezamabadi-Pour | Edge detection using ant algorithms[END_REF][START_REF] Tian | An ant colony optimization algorithm for image edge detection[END_REF]. It has, to our knowledge, never been used as accumulator of local orientations.

The main idea of the proposed method is that a certain number of ants are uniformly placed on the rough log-end image. They can freely move on the image and use normal values as pheromones. Their final goal is the pith. Each ant iteratively lays down pheromones as it moves towards the pith, and they all move towards the pith area where there is a high quantity of pheromones. The process is summarized in Alg. 1. Hereafter, we describe in details the different steps of the proposed ACO algorithm.

Initialization (Line 1 and 2 in Alg. 1) Let K × K be the number of ants and π (t) be the pheromone matrix at iteration t. At the beginning, K ants are placed in an uniform grid on the pre-processed image, and the pheromone matrix π (0) is initialized with random values drawn from a normal distribution.

Computation of probabilistic transition matrix (Line 5 to 9 in Alg. 1) An ant can move randomly with a probability that evolves during the optimization.

(u, v) */ foreach (u, v) ∈ η t k do η t k (u, v) = 1 √ (u-a) 2 +(v-b) 2 +1
/* Compute the probabilistic transition matrix ρ of the ant to move towards a position (u, v)

*/ foreach (u, v) ∈ τ t k do τ t k (u, v) = (π t (u,v)) α (η t k (u,v)) β i,j (π t (i,j)) α (η t k (i,j)) β /*
a d = 1 5 t+1 k=t-3 d k if a d < ε then Break Return π (N )
At iteration t, the probability for an ant k, currently at position (a, b), to move to the position (u, v) is defined by:

τ t k (u, v) = π t (u, v) α η t k (u, v) β i,j π t (i, j) α η t k (i, j) β ( 4
)
where τ t (u, v) is the amount of pheromone at (u, v), η t k (u, v) is the desirability of the k th ant to move towards (u, v), and equal to the inverse distance from (u, v) to (a, b):

η t k (u, v) = 1 (u -a) 2 + (v -b) 2 + 1 (5)
The desirability can be seen as a weighting of pheromone matrix. It aims to ensure ants having a higher probability to move towards local maxima and not towards the global one. α and β are respectively parameters to control the influence of τ t k (u, v) and η t k (u, v). The ratio α β allows to modify the behavior of ants; a high ratio value leads ants to move more quickly to the pheromone peaks, while a low value leads ants to continue to explore areas in image.

Pheromone deposit (Line 12 to 16 in Alg. 1) Each ant is the centre of an image block cluster. The cluster consists of n × n blocks, and each block has a size of ω × ω pixels. For each block cluster, the median value of local orientations (see Section 2.2) is considered as the block orientation. Then a line is drawn according to the orientation and passing through (u, v). All elements of the deposited pheromone matrix ρ t k along the line are incremented by 1. Indeed, depositing pheromones along the whole line allows to include lines intersections which could not happened if pheromone deposit is locally done. Fig. 3 illustrates this step of pheromone deposit of an ant.

Updating the pheromone matrix (Line 17 in Alg. 1) The pheromone matrix is updated once all ants have moved:

π t+1 = (1 -γ)π t + K k=1 ρ t k (6)
where ρ t k is the deposited pheromone matrix of the k th ant at iteration t, and γ is the rate of pheromone evaporation; the higher γ is, the faster pheromones are removed.

The process is repeated maximum N times (Line 3 in Alg. 1). In order to reduce computational time, the pheromone matrix π t is resized by a factor m comparing to the pre-processed image I. In other words, if

I is of size H × W , then π t is of size H m × W m .
Fig. 3: Left: Image with the ants in yellow crosses, the pith in red, a cluster of 3 × 3 blocks of size 8 × 8 pixels and normals (according to the cluster) in green lines. Right: Normal accumulation by pheromone matrix with the considered ant in yellow cross, the 9 green lines corresponding to the normal directions of the block clusters and the pith estimation in red.

Pith extraction

To extract the pith position from the pheromone matrix π, we take the barycentre of all the pixels above κ * max(τ ) in π. Indeed, taking the maximum value of accumulation is less robust than the barycentre of the highest values. Note that the higher κ is, the more sensitive to small variations the pith estimation is. Furthermore, we introduce an early stopping criterion in Alg. 1. At each iteration t, we estimate the pith location and compute the distance between the current and the last estimation. Instead of running for N iterations, the algorithm could stop as soon as the average of the last five distances falls below a threshold ε.

Code sources

Download and installation

The proposed method is implemented in Matlab 2019b and C++ using the open source library OpenCV3 (OPEN Computer Vision). Both implementations are available at the github repository: https://gitlab.com/Ryukhaan/treetrace/-/tree/master/pith

The installation is done with a cmake4 procedure (see README.md 5 ). In the following, we focus on the C++ implementation.

Description and usage

The repository has four packages:

• aco computes the Ant Colony Optimization algorithm for one image;

• normals computes normal accumulations using Bresenham lines;

• orientation computes local orientations for one image;

• ui manages the display (pheromones, ants position on image, and so on).

Once the installation is done, the executable file is in the build directory and named AntColonyPith. More details about the options are given in the command line helper.

.

/AntColonyPith --help

The options can be provided in two ways:

• using command line with usual options,

• providing a JSON file with all parameters (an example of JSON file, namely parameters.json 6 , is provided within the repository). • Output: Two files are created. The first one consists of the detected pith position in CSV format. The second one is an image of the input image with the detected pith denoted by a cross.

Experimental results

Experiments on real images

We experiment our algorithm on two datasets obtained from Douglas fir trees: Besle consists of 65 images and BBF consists of 40 images (see Fig. 4 for some examples). RGB images are converting into grayscale using the usual weighted method. Both datasets include the raw log ends taken in the forest or log yard, the images contain different disturbances such as sawing marks, dirt and light variations. Some visual results are shown in Fig. 5 by running the proposed method with the default parameters. The values of default parameters are given in the file parameters.json 6 . Further experiments about computation time and algorithm convergence are presented in the next sections. Note that we take the average value overall experiments on images in both datasets.

Accuracy of the method

For our implementation, we process twice the described method Section 2. RBG image are converted into grayscale (with the function imread and the option IMREAD GRAYSCALE from Opencv). The first run is to coarsely estimate the pith while the second run is for a precise pith estimation. For the first run, we split the image into 4 × 4 sub-images to manage the sawing marks removal (see Section 2.1). After retrieving the first pith estimation, this latter is converted back to coordinate of the original image. We select a sub-image of size 512 × 512 pixels centered on it and process again the algorithm (including the preprocessing without subdivision). Ground truths were done by two operators. Each operator independently, for each image, pointed the pith. The truth is the average of these two measures.

To determine the parameters' values, we manually minimized over the whole BBF dataset the sum of distances between ground truths and results. Then, -K = 16: the number of ants K × K; -α = 2.0: the control of the pheromone influence in (4); -β = 1.0: the control of the heuristic influence in (4); -γ = 0.07: the evaporation rate in [START_REF] Entacher | Towards an automated generation of tree ring profiles from ct-images[END_REF]; -m = 5 for the first run then m = 2 for the second one: how many pixels an element in the matrix τ stands for; -n = 3: the size of the blocks cluster (see Section 2.3); -ω = 8: the size in pixels of a block (see Section 2.3); -κ = 0.8: threshold to the barycentre (see Section 2.4) -ε = 2 for the first run then ε = 0.5: the thresh to early stop the algorithm.

-N = 50: the maximum number of iterations;

We have compared our results with [START_REF] Kurdthongmee | A comparative study of the effectiveness of using popular dnn object detection algorithms for pith detection in cross-sectional images of parawood[END_REF], [START_REF] Kurdthongmee | A fast algorithm to approximate the pith location of rubberwood timber from a normal camera image[END_REF] and [START_REF] Schraml | Pith estimation on rough log end images using local fourier spectrum analysis[END_REF] on our datasets. For the algorithm of Kurdthongmee et al., [START_REF] Kurdthongmee | A fast algorithm to approximate the pith location of rubberwood timber from a normal camera image[END_REF], we get optimized parameters with a Table 2: Average, standard deviation, minimum and maximum between ground truths and estimated piths by our method and methods of [START_REF] Schraml | Pith estimation on rough log end images using local fourier spectrum analysis[END_REF][START_REF] Kurdthongmee | A fast algorithm to approximate the pith location of rubberwood timber from a normal camera image[END_REF][START_REF] Kurdthongmee | A comparative study of the effectiveness of using popular dnn object detection algorithms for pith detection in cross-sectional images of parawood[END_REF] methods (in mm) and average time to proceed one image (in ms).

Besle

Mean subregion of size 24 × 24 pixels and a quantization factor of 12. We also used optimized parameters of our dataset for Schraml and Uhl algorithm [START_REF] Schraml | Pith estimation on rough log end images using local fourier spectrum analysis[END_REF]. For the comparison with the DNN [START_REF] Kurdthongmee | A comparative study of the effectiveness of using popular dnn object detection algorithms for pith detection in cross-sectional images of parawood[END_REF], we have done a twofold cross-validation. For each imageset, half of images have been used for the training and the other half for the validation. Two models were trained for each imageset by inverting the training and the validation sets. Ground truths consist of a square of 300 × 300 with the pith position at the center. A data augmentation have been processed on-the-flight (i.e. each time each image was transformed before passing through the DNN). The DNN hyperparameters were the same as [START_REF] Kurdthongmee | A comparative study of the effectiveness of using popular dnn object detection algorithms for pith detection in cross-sectional images of parawood[END_REF], only input size have been modified which is 576 × 432 for both imageset (the ground truth is resized according to that). The DNN returns a box with a probability of finding a pith in it. The predicted pith is the center of the box with the highest probability. To compare each method, we have aggregated all predictions from trained models (which gives us predictions for all images).

Table 2 presents a statistical analysis of the three algorithms on our datasets. The deep learning method is the fastest but drawbacks are the learning time and the creation of dataset with ground truths. Our method is, in average, 5 times faster than [START_REF] Schraml | Pith estimation on rough log end images using local fourier spectrum analysis[END_REF] and can be easily parallelized. We can observe that both our method and [START_REF] Schraml | Pith estimation on rough log end images using local fourier spectrum analysis[END_REF] are more accurate than [START_REF] Kurdthongmee | A comparative study of the effectiveness of using popular dnn object detection algorithms for pith detection in cross-sectional images of parawood[END_REF][START_REF] Kurdthongmee | A fast algorithm to approximate the pith location of rubberwood timber from a normal camera image[END_REF]. The results [START_REF] Kurdthongmee | A comparative study of the effectiveness of using popular dnn object detection algorithms for pith detection in cross-sectional images of parawood[END_REF] are worse on BBF imageset, this may be due to the small number of images in it. Fig. 6 presents boxplots for our method, [START_REF] Schraml | Pith estimation on rough log end images using local fourier spectrum analysis[END_REF] and [START_REF] Kurdthongmee | A comparative study of the effectiveness of using popular dnn object detection algorithms for pith detection in cross-sectional images of parawood[END_REF] to better illustrate the differences between them. We excluded [START_REF] Kurdthongmee | A fast algorithm to approximate the pith location of rubberwood timber from a normal camera image[END_REF] since the results are less accurate than the three others. For Besle imageset, our method is a little less accurate than [START_REF] Schraml | Pith estimation on rough log end images using local fourier spectrum analysis[END_REF]. [START_REF] Kurdthongmee | A comparative study of the effectiveness of using popular dnn object detection algorithms for pith detection in cross-sectional images of parawood[END_REF] is even a little less accurate and presents one outlier (7.61 mm). Its first and third quartiles are higher than our method and [START_REF] Schraml | Pith estimation on rough log end images using local fourier spectrum analysis[END_REF]. For BBF imageset, [START_REF] Schraml | Pith estimation on rough log end images using local fourier spectrum analysis[END_REF] has one outlier (7.6 mm) and [START_REF] Kurdthongmee | A comparative study of the effectiveness of using popular dnn object detection algorithms for pith detection in cross-sectional images of parawood[END_REF] has four outliers above 10 mm.

Effect of parameter changes on computation time

Hereafter, the experiments are carried out to analyse the effects of parameter changes. First, we focus on computation time then on precision and convergence. Let estimate the time complexity for one iteration in Alg. 1. Let I be the input image of size H × W . According to Alg. 1, there are K 2 ants, and they can freely move on I. Let now estimate the time complexity for an ant, namely the k th ant. Firstly, the desirability matrix η and the probabilistic transition matrix ρ associated to the k th ant are computed in O( HW m 2 ). Secondly, for the pheromone deposit, we must recall that each ant is the centre of an image block cluster which consists of n×n blocks and each block is of size ω ×ω. Therefore, to estimate the block orientation, we compute n 2 times the median of an array of ω 2 pixels. This operation is done by sorting the array and costs O(n 2 ω 2 log ω). We also compute n 2 times the deposited pheromone matrix ρ t k along the directional line l. In the worst case, the length of l is equal to ( H m ) 2 + ( W m ) 2 . In other words, the pheromone matrix update is done in O( n 2 m √ H 2 + W 2 ). Finally, the ant's position is updated in O(1). Once each ant has moved, the pheromones matrix π t is updated in O( HW m 2 ). Therefore, the total time complexity for one iteration is:

O K 2 n 2 ω 2 log ω + n 2 m H 2 + W 2 + 1 m 2 HW + 1 + 1 m 2 HW (7)
We can simplify this equation by keeping only the main input parameters, which are K, H and W . The total time complexity is therefore: We now validate this theoretical time complexity by the experiments. Firstly, it can be seen in Eq. ( 7) that the method is quadratic with respect to the size of the input image I. This is confirmed by Fig. 7, we have the computation time for one iteration according to I.

O K 2 H 2 + W 2 + HW + 1 (8)
From Eq. [START_REF] Fabijanska | Towards automatic tree rings detection in images of scanned wood samples[END_REF], the method has a linear time complexity with respect to the number of ants K 2 . Indeed, Fig. 8 (a) shows the computation time according to K. The number of ants quadratically increases, and thus the computation time.

Still in Eq. [START_REF] Fabijanska | Towards automatic tree rings detection in images of scanned wood samples[END_REF], the computation time decreases as m increases. Fig. 8 (b) shows the computation time according to m. It can be seen that having a value of m higher than 1 is really computationally helpful. Indeed, the matrix π is widely used during the process, from the pheromones deposit to the pheromones updates (reducing the size of η, τ and π). Let now look at two parameters: the block size ω × ω and the number of clusters n × n around the ant. First, the computation time according to the block size is shown in Fig. 8 (c). As the number of block is at least one, the computation time does not start at 0. For a block size of 3 × 3 it takes in average 243ms, while for a block size of 11 × 11 it is 259ms. Note that the computation time depends on the length of the line l used to update ρ (which depends on local orientations). For n, it should be quadratic. As n should be at least one, the computation time does not start at 0. Fig. 8 (d) shows the computation time according to n, and it is nearly quadratic. Again, this is due to the length of l.

Effects of parameter changes on the convergence

We now analyse the influence of parameters regarding the convergence of the algorithm. More precisely, the algorithm converges if the Euclidean distance between two successive estimations does not vary more than one pixel. We compute the average variations at each iteration over the whole set of images. We also fit a curve ae bx + c for each change in the value of parameters. It is done by First, we focus on the number of ants K × K. Fig. 9 shows the convergence speed according to K. One can see that the more ants are, the fastest the algorithm converges. With only four ants, the algorithm fluctuates between some positions. However, a high number of ants does not speed up convergence but makes the convergence point more stable (the parameter c is lower with a high value of K).

Let us now focus on the block's size ω. Fig. 11 shows the convergence speed according to ω. It seems that large block speeds up convergence but not as sharply as K. It is observed that an increasing in ω seems to slow down the convergence. This is due to an increase in the deposited pheromones. Indeed, the higher ω is, there more pheromones in the wrong places are. Therefore, it requires more iterations to remove those pheromones. Let us now consider the number of blocks n. Fig. 10 shows the convergence speed according to n. Contrary to intuition, a large number of blocks does not lead to an important acceleration of convergence, but it slows the algorithm down. This could be due to a less accurate local orientation with larger blocks.

Finally, we look at m. As a reminder, the higher m is, the smaller τ is. Fig. 12 shows the convergence speed according to m. It appears that m slightly speeds up the convergence. Raising the value of m causes a slightly decelera-Fig. 11: Variation for both axis between the pith estimation at time t and at time t + 1 according to the block size n. tion the convergence. It could be explained by the fact that with a small τ the small variations in the local orientation are not considered when ants deposits pheromones.

Limit cases

Our algorithm relies mainly on tree ring analysis. In other words, if the tree rings are not well presented in the input log-end image, then the detection result could be inaccurate. Fig. 13 shows some examples in which tree ring analysis is difficult and leads to a wrong pith estimation by using the proposed method. These two examples are not from both datasets. Indeed, our algorithm works very well on introduced datasets (no outliers). In both examples, it can be observed that tree rings are barely visible that makes difficult their analysis. Indeed, the pre-processing step may remove information about tree rings, which leads our algorithm to an inaccurate pith estimation. Furthermore, there are many other disturbances on such images; for instance, the sawing marks are not straight lines as we assumed in Section 2.1, or the presence of log-tree cracks.

Image credits

All images used in this paper are from the French National Research Agency, in the framework of the project TreeTrace, ANR-17-CE10-0016. Some samples (images in Fig. 5 and Fig. 13) are available for testing on the github repository 7 .

Conclusion

In this paper, we presented a probabilistic method for detecting pith position on digital images of rough, untreated log ends. More precisely, the proposed method is based on ant colony optimization (ACO) to robustly accumulate the normals of ring tree, then the pith location is extracted from this accumulation space as barycenter of points above a threshold. The experiments demonstrated that the proposed method provides not only an accuracy pith estimation (a distance of less than 5 mm from the ground truths), but also is efficient in computation time and it could be used in real-time applications. In addition to the implementation of the method, an online demonstration is available for testing at: https://ngophuc.github.io/ACO_PithDetection_IPOLDemo It could be noticed that the algorithm has many parameters. Though, they are set with default values and allow a good performance on tested images. Generally, the algorithm provides a very accurate pith estimation. A study on the role and effect of the different parameters is addressed in the paper for a better understanding of the parameters on the presented method. Based on this study, a perspective is to reduce the number of parameters and further provide an automatic approach to determine the best parameters adapted to a given image or a set of images of same characteristic. Furthermore, for reducing computation time of the algorithm, parallelism should be considered in future work.

Fig. 1 :

 1 Fig. 1: Examples of image for pith detection: (a-b) CT images, (c-e) digital images captured in log-yard. Digital images taken in realistic environments may contain disturbances: (c) light condition, (d) soiling and dirt, (e) sawing marks.

  Detected line (e) Sawing mark removal (f) Local orientation

Fig. 2 :

 2 Fig. 2: Removing sawing marks and estimation of local orientations: (a) input grayscale image, (b) FFT of (a), (c) threshold (b) with λ = 0.875, (d) detected line after convoluted with a Gaussian of σ = 6, (e) image after removing sawing marks, and (f) estimated local orientations (red lines) on (e).

Algorithm 1 :

 1 Ant colony optimization for normal accumulation Input: The estimated local orientations I θ The number of ants K × K The maximum number of iterations N The number of block clusters around an ant n × n The size of each block cluster ω × ω pixels Output: The pheromone matrix π Variables : η t k : The desirability matrix of the k th ant at iteration t τ t k : The probabilistic transition matrix of the k th ant at iteration t ρ t k : The deposited pheromone matrix of the k th ant at iteration t π t : The pheromone matrix at iteration t 1 Initialize the positions of the K ants 2 Initialize the pheromone matrix π 0 3 for t = 1 . . . N do 4 for k = 1 . . . K 2 do 5 Let (a, b) be the position of the k th ant /* Compute the deposited pheromone matrix ρ t of the k th ant */ 6 Let B be the n × n block clusters of size ω × ω centered at (u, v) in I θ 7 foreach h ∈ B do 8 Let (ox, oy) be the median value of local orientations of h 9 Let l be the line of orientation (ox, oy) passing through (u, v) Increase ρ t k by 1 along l /* Compute the desirability matrix η of the ant to move towards a position

  Move the ant according to the probabilistic transition matrix τ t k */ Let (x, y) be the position of the maximum probability in τ t k Move the ant to the position (x, y) /* Update the pheromone matrix π t after all K ants moved */ π t+1 = (1γ)π t + K k=1 ρ t k /* Early-stopping criteria */ /* Compute the new pith position */ The current pith position pt+1 is estimated according to π t+1 /* Compute the distance between the current and the last pith position */ dt+1 = pt+1pt 2 /* Compute the average of the last five distances */

•

  Input: The image to be processed; • Command Line: To run the program from the CODESOURCES/build ./AntColonyPith --input=path to image [list of parameters] ./AntColonyPith --input path to image --parameters path to parameters.json For instance, to run the program on harvest.jpeg with default parameters ./AntColonyPith --input ../../samples/harvest.jpeg To run the program on harvest.jpeg with 10×10 ants, α = 1.0 and without animation ./AntColonyPith --input ../../samples/harvest.jpeg --ant=10 --alpha=1.0 --animated=false or ./AntColonyPith --input ../../samples/harvest.jpeg -n 10 -a 1.0 --animated=false To run the program on harvest.jpeg with parameters in parameters.json ./AntColonyPith --input ../../samples/harvest.jpeg --parameters ../AntColonyPith/parameters.json

Fig. 4 :

 4 Fig. 4: Examples from Besle (the first two rows) and BBF (the last row) datasets.

Fig. 5 :

 5 Fig. 5: Pith position (black cross) detected by the proposed method on raw log-end images using default parameters.

Table 1 :

 1 Pre-processing parameters for both steps. H is the height of the Fourier spectrum. values on Besle dataset. Tables 1 shows the parameters obtained for the preprocessing. Parameters for the ACO-based algorithms are set as follows (the values are the same for both phase unless otherwise indicated):

Fig. 6 :

 6 Fig.6: Boxplots of distance between ground truths and pith estimation (in mm) for[START_REF] Schraml | Pith estimation on rough log end images using local fourier spectrum analysis[END_REF][START_REF] Kurdthongmee | A comparative study of the effectiveness of using popular dnn object detection algorithms for pith detection in cross-sectional images of parawood[END_REF] and our method.

Fig. 7 :

 7 Fig. 7: Computation time for one iteration according to the size of input image (in ms).

Fig. 8 :

 8 Fig. 8: Computation time for one iteration according to (a) the number of ants K, (b) the size of τ , (c) the block size ω and (d) the cluster size n.

Fig. 9 :

 9 Fig.9: Variation in both axis between the pith estimation at time t and at time t + 1 according to the number K × K of ants.

Fig. 10 :

 10 Fig. 10: Variation in both axis between the pith estimation at time t and at time t + 1 according to the cluster size ω.

Fig. 12 :

 12 Fig. 12: Variation in both axis between the pith estimation at time t and at time t + 1 according to the parameter m.

Fig. 13 :

 13 Fig. 13: Examples of wrong pith estimation. Left: The tree rings are in low resolution and the sawing marks are not straight lines. Right: the sawing marks are straight lines but tree rings are in low resolution and really small.

.26 1.32 0.44 4.63 1745

			StDev Min Max Time (ms)
	[21]	2.29 0.98 0.39 4.96	8344
	[13]	25.06 21.23 2.15 92.44	667
	[12]	2.88	1.67 0.87 7.61	138
	Our method 2.34	1.02 0.46 5.04	1611
	BBF	Mean StDev Min Max Time (ms)
	[21]	2.39	1.48 0.49 7.59	8660
	[13]	38.38 38.91 3.14 232.74	721
	[12]	12.69 53.55 0.50 341.92	186
	Our method 2

https://opencv.org/

http://www.cmake.org

https://gitlab.com/Ryukhaan/treetrace/-/blob/master/README.md

https://gitlab.com/Ryukhaan/treetrace/-/blob/master/pith/c++/ AntColonyPith/parameters.json

Acknowledgment

This research was made possible by support from the French National Research Agency, in the framework of the project TreeTrace, ANR-17-CE10-0016.


Florian Delconte
email: florian.delconte@loria.fr

Phuc Ngo

Isabelle Debled-Rennesson

Bertrand Kerautret

Van-Tho Nguyen

Thiery Constant

Tree Defect Segmentation using Geometric Features and CNN

Keywords: Wood surface defects, Defect segmentation, Relief map, LIDAR, Centerline, U-Net

Estimating the quality of standing trees or roundwood after felling is a crucial step in forest production trading. The on-going revolution in the forest sector resulting from the use of 3D sensors can also contribute to this step. Among them the terrestrial lidar scanning is a reference descriptive method offering the possibility to segment defects. In this paper, we propose a new reproducible method allowing to automatically segment the defects. It is based on the construction of a relief map inspired from a previous strategy and combining with a convolutional neural network to improve the resulting segmentation quality. The proposed method outperforms the previous results and the source code is publicly available with an online demonstration allowing to test the defect detection without any software installation.

Introduction

In the domain of biological image processing, the wood structures are often exploited to address various objectives, for instance, species identification [START_REF] Carpentier | Tree species identification from bark images using convolutional neural networks[END_REF], wood quality estimation [START_REF] He | Application of deep convolutional neural network on feature extraction and detection of wood defects[END_REF], tree microhabitats identification [START_REF] Rehush | Identifying tree-related microhabitats in TLS point clouds using machine learning[END_REF], tracability [START_REF] Schraml | Matching score models for hyperspectral range analysis to improve wood log traceability by fingerprint methods[END_REF], or plant growing analysis [START_REF] Gélard | 3D leaf tracking for plant growth monitoring[END_REF]. These various applications rely on different image acquisition modalities such as classical 2D bitmap images (including hyperspectral images), 3D point cloud (from multi-view stereo 3D or LiDAR scan) or 3D volumetric images (medical X-Ray CT scanner [START_REF] Krähenbühl | Knot segmentation in 3D CT images of wet wood[END_REF] or ultra sound [START_REF] Gilbert | Use of Sonic Tomography to Detect and Quantify Wood Decay in Living Trees[END_REF]).

The aim of this work is to detect defects located on the trunk surface of living tree (see teasing Figure 1). Various types of defects are identified by biology experts (also called singularity) depending on their origins and their development stages (burls, branch scar, picot, . . . ). Figure 2 illustrates samples of defects on beech and oak species. The detection of such a structure is a key point for the value determination and the optimization of the transformation taking knotiness or aesthetics into consideration. The defect detection on living tree is not an easy task in the image processing domains, since each type of singularities presents numerous geometric variations both inter or intra species. Figure 2 (a,b) illustrates a same defect type on a same species but presenting a very different geometric shape. In this work, the defect detection is addressed by using 3D scan of trunk (as illustrated in Figure 1). Schütt et al. can be considered as the pioneers to exploit the terrestrial LiDAR data for tree defect detection [START_REF] Schütt | Approaches for recognition of wood quality of standing trees based on terrestrial laserscanner data. Laser-scanners for forest and landscape assessment[END_REF]. The authors proposed to localise singularity by combining 3D terrestrial LiDAR with 2D images. After using a cylindrical coordinate transformation, a neural network is trained and used to extract the singularity areas. The method was promising, however the extraction process needs a potential interactive correction and no details are given to reproduce the method nor the result quality measures. Thomas et al.

proposed later an automatic method to detect severe surface defects by using a 2D circle fitting algorithm [START_REF] Thomas | Automated detection of severe surface defects on barked hardwood logs[END_REF]. The method was then improved with a parallel implementation in order to reduce execution time [START_REF] Thomas | Using parallel computing methods to improve log surface defect detection methods[END_REF]. One of the limitation of such an approach, was the minimum detectable defect size was 12.7 centimeters and with a relief higher than 1.27 centimeters. Answering to the previous limitation, Kretschmer et al. [START_REF] Kretschmer | A new approach to assessing tree stem quality characteristics using terrestrial laser scans[END_REF] added more geometry in the singularity detection by using a cylinder fitting based approach. Using a tree reconstruction method [START_REF] Pfeifer | Automatic reconstruction of single trees from terrestrial laser scanner data[END_REF], a series of cylinders is fitted according to the wood main axis and each 3D points is associated to its cylinder part. Such an association allows to generate a distance map that is used to extract manually the defects. The proposed strategy is not automatic, however it allows to detect smaller defects with near 4.3 centimeters for the minimal size and 2 centimeters for the relief height. Observing that the wood trunk does not always fit perfectly a cylinder, Nguyen et al. [START_REF] Nguyen | Segmentation of defects on log surface from terrestrial lidar data[END_REF] proposed an automatic patch based method that allows to better follow the trunk geometry. The main algorithm relies on the recovering of the trunk centerline [START_REF] Kerautret | Centerline detection on partial mesh scans by confidence vote in accumulation map[END_REF] allowing to avoid the cylinder fitting step of previous works.

The methods described in the previous part are designed specifically for the wood defect estimation but it could be interesting to mention other general approaches that are exploited for surface-crack detection in an industrial context. For instance, Tabernik et al. [START_REF] Tabernik | Segmentation-based deep-learning approach for surface-defect detection[END_REF] proposed a segmentation-based deep-learning method to detect surface anomaly. Even if the considered images differ from the tree defect context, their strategy could be interesting to adapt since the proposed architecture does not request training with numerous images. In other context of the train industry, the defect on rail surface were analysed through [START_REF] Faghih-Roohi | Deep convolutional neural networks for detection of rail surface defects[END_REF]. Like the previous works, their approach was based on the deep learning and can detect various defects like weld, squat or joint. Finally, we can refer to another application of micro cup surface inspection from a confocal laser microscopy images that exploits neural network to detect defects on very textured images [START_REF] Weimer | Learning defect classifiers for textured surfaces using neural networks and statistical feature representations[END_REF].

Following the previous strategy introduced by Nguyen et al. [START_REF] Nguyen | Segmentation of defects on log surface from terrestrial lidar data[END_REF], we propose a new method based on the construction of a new relief map image combined with a convolutional neural network (CNN) to precisely segment tree defects. The main overview of the method is given on the Figure 3. The relief map is constructed from the input mesh (upper part of Figure 3) and the convolutional neural network is exploited to segment defect area (lower part of Figure 3) that can be visualized on the original input mesh.The first step of the new approach is described in the following section (Section 2) with the overview of the centerline detection, followed by relief map image construction. Associated to this representation, the segmentation process based on the U-Net architecture is introduced in Section 3. The experiment part presenting the main results and reproducibility links are addressed in Section 4 before concluding.

Geometric tools

The proposed method relies on three geometric tools: (i) the centerline of the wood log, (ii) the reference and delta distances and (iii) the relief map. The first Fig. 3: Pipeline of the proposed method. two tools are introduced in [START_REF] Kerautret | Centerline detection on partial mesh scans by confidence vote in accumulation map[END_REF][START_REF] Nguyen | Segmentation of defects on log surface from terrestrial lidar data[END_REF] and summarized in following sections. The third tool, the relief map, is defined in Section 2.3. It permits to represent the input 3D points with a 2D map characterising the relief of the points, relatively to the centerline of the wood log and from a fitted tangent plane.

Centerline of the wood log

In [START_REF] Kerautret | Centerline detection on partial mesh scans by confidence vote in accumulation map[END_REF] a method is presented to extract the centerline of 3D shapes using solely partial mesh scans of the shapes. The centerline is a polyline with several small segments (see Figure 4 (b,c)). It is obtained by constructing an accumulation map from input faces and normal vectors (see Figure 4 (a)) and by filtering it with a confidence vote. Since the method inputs are only a set of faces, the centerline can also be recovered both from full and partial mesh (see Figure 4 (b,c)). The details of the algorithm are available in the associated reference [START_REF] Kerautret | On the implementation of centerline extraction based on confidence vote in accumulation map[END_REF] and on the GitHub repository: https://github.com/kerautret/CDCVAM Due to the non constant diameter of wood logs, a process of optimization must be done to obtain a smooth centerline of the wood logs. In [START_REF] Nguyen | Segmentation of defects on log surface from terrestrial lidar data[END_REF], the authors used a smoothing process based on cubic spline. Note that the implementation details and reproductive evaluation can be found in the complementary work [START_REF] Nguyen | Algorithms and implementation for segmenting tree log surface defects[END_REF] with the GitHub repository: https://github.com/vanthonguyen/treelogdefectsegmentation

Segmentation of defects on the mesh of partial trunk

Thanks to the ability to compute the centerline even in the mesh built from only one scan of the trunk. As shown in Figure 8, the performance of the segmentation algorithm on these meshes was not far from the mesh of the entire trunk, which was scanned from 4 di erent points of view. The weak point when using just one scan is that the spatial resolution of the border of such a scan is lower than in the center area and there are more borders leading to fewer points in the corresponding cylindrical sectors. Thus, the information about the defects was less precise compared to 4 scans (see Figure 8 (a) and (b) for example). The Table 1 also reflect this finding as the precision is very high (0.910) while the recall is lower (0.659). The multiscan can also benefit from the reduction of occlusions (Figure 8 (g) and (h)). While in the one scan mesh, there are several large missing zones, on the 4 scans mesh, these zones are compensated by points from other 11

Reference and delta distances

In order to easily access the neighborhood of each point on the wood log surface, we work in cylindrical coordinates. A local coordinate system (C i , u i , v i , w i ) is defined for each segment S i of the centerline. A point P (x, y, z) in Cartesian coordinates corresponds to the cylindrical coordinates (r P , θ P , z P ) with:

r P is the distance between P and P , the projection of P on the segment S i of the centerline. z P is the height of P along the centerline.

θ P is the angle formed between the segment P P and the axis v i of local coordinate system associated to S i .

For more details of transformation in cylindrical coordinate, we refer the readers to [START_REF] Nguyen | Segmentation of defects on log surface from terrestrial lidar data[END_REF]. To correctly detect the local relief variation around each point P of the wood log, a rectangular neighborhood is studied, named patch P P , it is proportional to the size and circumference of the wood log (see Figure 5 (a) and [START_REF] Nguyen | Segmentation of defects on log surface from terrestrial lidar data[END_REF] for details). P P characterizes the shape of the log around the point P . The central straight line fitting the points of P P is calculated by a RANSAC based linear regression. The reference distance of the point P, rP , is the distance from this straight line to the centerline (see Figure 5 (b)). The difference between r P and rP represents the relief of the tree at a point. It is called the delta distance: δ P = r P -rP (see Figure 5 (b)). We use the delta distance in the next section to generate the relief maps.

Relief map

A relief map is a 2D representation of the tree mesh. It is obtained by firstly discretizing the cylindrical point space and by secondly, completing the missing information with a multi-resolution analysis. This map is used to segment the defects. It must also allow a reverse operation, i.e., compute from pixels of the map the corresponding 3D points of the mesh.

Cylindrical space discretisation

The relief map represents the unfolding of the wood log. The width of the map is the circumference of the trunk, i.e., 2π * r m with r m the average radius of the trunk. The height of the map is the height of the trunk, obtained by subtracting the z component of the point having the maximum height and the one of minimum height. Each point of the tree mesh is associated to a cell of the relief map. We then calculate a value to represent all the points of a cell. The chosen value is the maximum value of the delta distances of the points associated to the cell. An illustration is given in Figure 6.

The two maps are generated from the same input mesh. On the left, intensity of the pixels is calculated from the distance to the centerline. On the right, intensity of the pixels is calculated from the delta distances. With the map obtained by the distance to the centerline, we can observe the artefacts, the yellow and red traces, due to the non-cylindrical tree, while using the delta distance, these traces disappear and the defects become more visible. Figure 7, on the first line, provides several relief maps deduced from this process. Multi-resolution analysis It is possible that some cells of the relief map do not contain a value because no point is associated with it. To handle this, we propose a multi-resolution analysis to improve the obtained relief map. In the proposed process, for every empty cell, we reduce the resolutions by a factor 1 2 n , with n ∈ N, until the cell contains at least one point. It should be noticed that if n is too high, we may lost information of delta distance for the defect detection. Therefore, during the multi-resolution analysis, we fix the limite of n to four. Figure 8 illustrates the multi-resolution analysis to obtain a value in an empty cell. We consider in this example an array T with a resolution of 10 × 10. Some cells contain black dots, corresponding to the points previously discretized. In others, there is none, like the cell located in [START_REF] Carpentier | Tree species identification from bark images using convolutional neural networks[END_REF][START_REF] Carpentier | Tree species identification from bark images using convolutional neural networks[END_REF]. The multi-resolution analysis is illustrated by the colors of the borders of the table. Respectively, the black, red and blue borders correspond to resolutions reduced by a factor of 1 1 , 1 2 1 , 1 2 2 . We look for the information in the cell (2, 2), then (1, 1) of the red discretization, then (0, 0) of the blue discretization, etc ... The value in the (2, 2) cell is then the maximum of the delta distances of the points in the (0, 0) cell of the blue discretization. To represent the discretization in the form of an image we associate a gray intensity to the delta distance. This intensity is distributed on a fixed scale: 1cm equals ten gray levels starting from -5. The colored relief maps in this article are obtained by applying a color scheme from blue to red. The bottom row in Figure 7 shows the improvement brought to the relief maps with the multi-resolution analysis. We can see in Figure 9, results of the multi-resolution analysis centered on a branch scar type defect.

Segmentation with U-Net architecture

Hereafter, we process the detection of defects on tree barks using the previously obtained relief map. Note that the 3D problem of defect detection on tree bark surfaces becomes a 2D problem of relief map segmentation. More precisely, it is a binary-image-classification in which each pixel of the relief map will be classified as defect or not. We can observe in Figure 7, that the defects on tree barks may have arbitrary size, shape and orientation. Furthermore, the roughness of the tree bark, the variability of the defects on the same species and between the different species make the detection task difficult to automate by conventional segmentation algorithms.

Over the past few years, the deep-learning methods are becoming common and successful for the segmentation task with remarkable performance improvements. Indeed, they often achieve the highest accuracy rates on popular segmentation benchmarks comparing to the classical computer-vision approaches. Furthermore, the deep-learning algorithms can be adapted to different problems as they can learn the hidden high-level features from the image directly, and have capacity to represent and recognize the complex structures.

In this paper, we use a deep learning-based segmentation method, namely U-Net [START_REF] Ronneberger | U-net: Convolutional networks for biomedical image segmentation[END_REF], to detect tree bark defects with the relief maps as input. It should be mentioned that in the context of surface-singularity detection, several neural network architectures have been proposed (see Section 1). As stated in [START_REF] Ronneberger | U-net: Convolutional networks for biomedical image segmentation[END_REF], U-Net enables the model to be trained using a small number of samples, and to create a precise pixel-wise mask of interest objects in the images. Thus, it is a well-suited architecture for our segmentation problem. In the following, we describe a modified version of the original U-Net [START_REF] Ronneberger | U-net: Convolutional networks for biomedical image segmentation[END_REF] for detecting tree bark defects, and the generation of training data from the relief maps.

Segmentation network

The U-Net was first introduced in [START_REF] Ronneberger | U-net: Convolutional networks for biomedical image segmentation[END_REF] as a fully convolutional network (FCN) architecture for biomedical image segmentation, and it was designed for a precise pixel-wise segmentation. Recently, many variants of U-Net architecture have been proposed to address the segmentation of medical images, satellite images such as U-Net++ [START_REF] Zhou | Unet++: A nested u-net architecture for medical image segmentation[END_REF], KU-Net [START_REF] Wagner | Using the u-net convolutional network to map forest types and disturbance in the atlantic rainforest with very high resolution images[END_REF], TernausNet [START_REF] Iglovikov | Ternausnet: U-net with VGG11 encoder pre-trained on imagenet for image segmentation[END_REF], . . . U-Net is well-known for its performance to be trained with very few training images.

U-Net is an auto-encoder architecture. The encoder stage takes the input images and extracts features from objects in the image, then condenses them into smaller layers. These features are propagated in the decoder stage to produce a segmentation.

The original U-Net proposes to use convolution layers followed by max pooling for down-sampling in encoding part, while the decoding part consists of Fig. 10: The architecture for tree bark defect detection based on U-Net [START_REF] Ronneberger | U-net: Convolutional networks for biomedical image segmentation[END_REF]. up-sampling followed by concatenation with the corresponding layer of the encoding part. Each layer is followed by the rectify linear unit (ReLu) activation function, except the last one. This final layer is a 1x1 convolution followed by a pixel-wise soft-max over the final feature map. The cross entropy loss function is used to update weights of the network. The soft-max function redistributes the weights of the final layer of the network in the interval of [0, 1] modeling a probability distribution over predicted output classes. In total the network has 23 convolutional layers. More details of U-Net can be found in [START_REF] Ronneberger | U-net: Convolutional networks for biomedical image segmentation[END_REF].

In this paper, to address our problem of detecting tree bark defects which is defined as two class segmentation, we made several changes to the original U-Net. In order to reduce the over-fitting of the considered neural network, we apply a regularization technique, called drop-out. More precisely, we add two dropout layers in the encoder and decoder, with probability 0.5, to randomly drop some of the connections between layers. In addition, due to the dying ReLU problem [START_REF] Lu | Dying relu and initialization: Theory and numerical examples[END_REF] -i.e., the ReLU neurons become inactive and only output 0 for any input-the Leaky ReLu activation function is employed instead of ReLU from the original architecture. Finally, in the last layer, we use a Sigmoid activation function instead of soft-max function to ensure the output pixel values range between 0 and 1. For the training, we use input image of size 320x320 pixels. The proposed network architecture is illustrated in Figure 10.

Training data

The training dataset for U-Net framework is built from the relief maps generated from tree bark surfaces by the process described in Section 2. It should be mentioned that we only have 25 annotated meshes -i.e., 25 relief maps with ground-truths-of tree bark surfaces for the learning process (for more details, see Section 4.1). Due to this limited number of samples, two strategies have been adopted to augment the existing data while keeping the significant characteristics of learning objects which are tree bark defects. First and foremost, we split each relief map into patches of same size. For this, we carry out two types of cutting (see Figure 11). The first one aims to obtain samples centered on defects. More precisely, we perform the splitting of relief map according to the barycenter of the connected component associated to the defect in the annotated image. Note that if a defect is close to the border, a translation is applied to obtain a patch containing the defect and being included in the map. The second one collects samples that do not contain any defect so that the network can learn tree bark without defect. Some samples of extracted patches with and without defects are given in Figure 12. It should be mentioned that the size of the generated patches is limited by the width and height of our relief maps. As observed in Figure 7, the relief maps may have different sizes because of the discretization being made with respect to the circumference and height of the tree bark (see Section 2). For our framework, the patches are of size 320x320 pixels which is the largest size that could be extracted from the relief maps. The obtained images are then randomly separated into two subsets with a ratio of 7:3 for the training and validation of the network.

After this splitting process, different transformation techniques have also been used on the obtained patches for data augmentation. In particular, we consider the operations: rotation, vertical and horizontal flip, zoom and deletion of rectangular area randomly [START_REF] Devries | Improved regularization of convolutional neural networks with cutout[END_REF]. Note that this data augmentation is performed on the fly, i.e., during the training time.

Experiments

Dataset

We have two datasets for our experiments: INRAE1a and INRAE1b. IN-RAE1a contains 10 trunks of different species: beech (1), birch (1), elm (1), fir (2), red oak (2), wild cherry (2) and service wood [START_REF] Boulch | Convpoint: Continuous convolutions for point cloud processing[END_REF]. INRAE1b contains 15 meshes including alder (4), aspen (4), beech (1), birch (2), horn beam (1), lime (1), red oak [START_REF] Carpentier | Tree species identification from bark images using convolutional neural networks[END_REF]. The first dataset was used in [START_REF] Nguyen | Segmentation of defects on log surface from terrestrial lidar data[END_REF], experiments were carried out on INRAE1a to compare performance and robustness of our method with [START_REF] Kretschmer | A new approach to assessing tree stem quality characteristics using terrestrial laser scans[END_REF] and [START_REF] Nguyen | Segmentation of defects on log surface from terrestrial lidar data[END_REF] on different tree species. The relief maps of INRAE1a are given in Figure 14. The second dataset is used for the training and illustrated in Figure 13.

Both datasets have the ground-truths being made by hand-labeling defects. The ground-truth is given as a set of point indices associated to the defect. These indices are then used to generate the annotated maps for training the network.

Network training

The training process was first performed on the relief maps generated from 15 meshes of INRAE1b. We used the parameters recommended in [START_REF] Nguyen | Algorithms and implementation for segmenting tree log surface defects[END_REF] for computing these maps. After splitting the obtained relief maps into patches (see Section 3.2), we divide the patches of each map into two subsets: 70% for training and 30% for validation. This subdivision allows to have the same bark variability on training and validation. To summarize, the whole dataset has 265 images of size 320x320 pixels, it is partitioned into 204 images for training and 61 images for validation.

The training process is done on GPU (geforce RTX 2080Ti with 12Go RAM). The modified U-Net is implemented using Tensorflow 2. including rotation, vertical and horizontal flip, zoom and deletion of rectangular area randomly [START_REF] Devries | Improved regularization of convolutional neural networks with cutout[END_REF]. We used the Adam optimiser [START_REF] Kingma | Adam: A method for stochastic optimization[END_REF], and set the learning rate at 0.0001, the two parameters β 1 = 0.9 and β 2 = 0.99 (default values in Tensorflow 2.2). We trained our network for 40 epochs, each epoch comprised 63 steps with 10 images per batch. About the parameters of dropout rate δ and Leaky ReLu activation α, several values have been tested, and we come out with δ = 0.5 and α = 0.01 for the smallest loss function (binary cross entropy) on the validation. It should be mentioned that our training is quite fast, it takes about 14 seconds per epoch. In other words, the proposed architecture allows a high-quality segmentation and very fast training -about 10 minutes for the whole training process-with very small dataset. In particular, the prediction takes, on average, 451 miliseconds per map.

Experimental results

The first experiments were performed on INRAE1a. More precisely, the relief maps were generated for the 10 meshes of the dataset, then predicted by our network which is previously trained on INRAE1b. The output prediction is a gray-level image. We threshold this image at 0.5 to obtain a binary image in which the white pixels indicate the tree defect and black is not. The results are given in Figure 14.

To evaluate and compare the methods, we used the classic metrics: precision, recall and F measure (F1). For a fair comparison, we performed the evaluation measures on the mesh points, as done in [START_REF] Nguyen | Segmentation of defects on log surface from terrestrial lidar data[END_REF], but not the predicted maps. As described in Section 2.3, from pixel positions, we can easily retrieve the mesh point indices associated to the pixel, and identify those points for localizing the defects on the mesh. Table 1 shows the obtained results. Generally, the proposed method outperforms both cylindrical-based [START_REF] Kretschmer | A new approach to assessing tree stem quality characteristics using terrestrial laser scans[END_REF] and patch-based [START_REF] Nguyen | Segmentation of defects on log surface from terrestrial lidar data[END_REF] methods. We improve the detection performance, about 41% and 8% better in F1 measure comparing to [START_REF] Kretschmer | A new approach to assessing tree stem quality characteristics using terrestrial laser scans[END_REF] and [START_REF] Nguyen | Segmentation of defects on log surface from terrestrial lidar data[END_REF], respectively. Figure 15 shows a visual comparison on meshes of the results obtained by the proposed method and the patch-based method [START_REF] Nguyen | Segmentation of defects on log surface from terrestrial lidar data[END_REF]. We can see in Figure 15 (a) that our method tends to produce fewer false positives, but sometimes miss small defects. This may due to the fact that the network has not been trained with samples containing small defects. Figure 15 (b) is an example where the detection by our method covers better the form of defects than patch-based method [START_REF] Nguyen | Segmentation of defects on log surface from terrestrial lidar data[END_REF].

We carried out a second experiment to demonstrate the generalization of the network. Using the same parameters, we trained our CNN on mixed data of INRAE1a and INRAE1b, i.e., 25 meshes in total. We generated 5 different folds, each of which contains 20 meshes for training and 5 for validation. Table 2 INRAE1a

Patch method [START_REF] Nguyen | Segmentation of defects on log surface from terrestrial lidar data[END_REF] Cylinder method [START_REF] Kretschmer | A new approach to assessing tree stem quality characteristics using terrestrial laser scans[END_REF] Our method prec recall F1 prec recall F1 prec recall F1 summarizes the distribution of data and the results obtained. For each fold in Table 2, the meshes indicated in Mesh id correspond to test data, and the others are used for training. In this way, we ensure to test and compare our method with the others on all available data. It can be observed that, over the 5 folds, we generally obtain the best F1 measure, and almost better on the precision comparing to [START_REF] Kretschmer | A new approach to assessing tree stem quality characteristics using terrestrial laser scans[END_REF] and [START_REF] Nguyen | Segmentation of defects on log surface from terrestrial lidar data[END_REF]. The worst scoring result by the proposed method is Beech3 in the fold 5, we are at 0.336 for F1. Though, this score is comparable to the best score of 0.493 obtained by [START_REF] Nguyen | Segmentation of defects on log surface from terrestrial lidar data[END_REF].

Note that the measured results of the methods [START_REF] Nguyen | Segmentation of defects on log surface from terrestrial lidar data[END_REF] and [START_REF] Kretschmer | A new approach to assessing tree stem quality characteristics using terrestrial laser scans[END_REF] for the dataset INRAE1a (see Table 2) are from [START_REF] Nguyen | Segmentation of defects on log surface from terrestrial lidar data[END_REF], while the results on the new dataset in Table 2 are obtained by using the source code from the GitHub repository with the recommended parameters described in [START_REF] Nguyen | Algorithms and implementation for segmenting tree log surface defects[END_REF][START_REF] Nguyen | Segmentation of defects on log surface from terrestrial lidar data[END_REF].

An online demonstration for testing the proposed method is available at: https://kerautret.github.io/TLDDC/ 5 Source code to reproduce results

Global view

The source code to reproduce the results presented in the article, including relief map and segmentation, is available at the GitHub repository: https://github.com/FlorianDelconte/TLDDC The repository is composed of different file and directories :

-The directory Centerline contains the centerline code.

-The directory examples contains INRAE1e and INRAE1b meshes, each mesh is accompanied by two files indicating the ground-truth location of

Our method Patch method [START_REF] Nguyen | Segmentation of defects on log surface from terrestrial lidar data[END_REF] Our method Patch method [START_REF] Nguyen | Segmentation of defects on log surface from terrestrial lidar data[END_REF] (a) Redoak2 (b) WildServiceTree

Our method Patch method [START_REF] Nguyen | Segmentation of defects on log surface from terrestrial lidar data[END_REF] Our method Patch method [START_REF] Nguyen | Segmentation of defects on log surface from terrestrial lidar data[END_REF] (c) WildCherry1 (d) WildCherry2

Fig. 15: Comparison on mesh. Yellow is true positive, red is false negative and green is false positive compared to the ground-truth.

defects using mesh faces and mesh points, suffixed by -groundtruth.id and -groundtruth-points.id respectively. -The directory models contains the trained models in the paper: five files of .hdf5 extension for the corresponding k-fold, and one named KFoldAssociation for a relation file between mesh example and the k-fold. -The directories mesures and run contain python and bash scripts to directly reproduce the results of this article. -The code for generating the relief maps is found in UnrolledMap.h, Un-rolledMap.cpp, DefectSegmentationUnroll.h and DefectSegmentationUnroll.cpp.

Fold number Mesh id Patch method [START_REF] Nguyen | Segmentation of defects on log surface from terrestrial lidar data[END_REF] Cylinder method [START_REF] Kretschmer | A new approach to assessing tree stem quality characteristics using terrestrial laser scans[END_REF] Our method prec recall F1 prec recall F1 prec recall F1 Table 2: Comparison results on mixed dataset. Overall row is computed from the sum of T P ,T N ,F P and F N on all the tested meshes.

Installation

For compilation process, the program required this libraries to be installed:

- Instructions for installing on ubuntu 20.04 and debian 10 have been tested and are detailed on GitHub. Once the dependencies are installed and the sources downloaded, the code is built using the following commands: cd TLDDC mkdir build cd build cmake .. -DDGtal DIR=/path/to/DGtalSourceBuild make Two executable files are generated in the build directory:

segunroll allows to generate relief maps from meshes.

-segToMesh allows to project the segmentation of the defects of the relief map towards the mesh.

Usage

To generate the relief map, run the following command from build:

./segunroll -i InputMesh

[-h] [-n] [CenterlineParameters] [ReliefMapParameters]
With -InputMesh is the path to a trunk mesh.

--h is the option for the command line helper.

--n is the option allows to invert the normals of the faces of the meshes 5 .

-CenterlineParameters contains the parameters of centerline computation (--accRadius, --trackStep, --binWidth, --patchWidth, --patchHeight, --voxelSize).

They are set by default with the recommended values in [START_REF] Nguyen | Algorithms and implementation for segmenting tree log surface defects[END_REF]. -ReliefMapParameters contains the parameters for the relief map (--decreaseFactor, --grayscaleOrigin, --intensityPerCm). They are also set with default values.

The following files are created after executing the command:

centerline.off : the generated centerline.

discretisation.txt: the discretization map.

output.pgm: the generated relief map.

To segment the bark tree defects, run the following command from build: -PathToModel contains the model file, one of the five in the models directory.

-Treshold contains the threshold ([0; 255]) to apply on the network prediction.

The following files are created after executing the command:

-outputSEG.pgm: the prediction result by the network.

-outputSEGTRESH.pgm: the segmented image after thresholding.

To project the segmented result on the input mesh, run the following command from build:

./segToMesh -i InputMesh

With InputMesh contains the path to the same mesh used to generate the relief map. The following files are created after executing the command:

output-defect.id : the file containing the id of the points of the mesh belonging to a defect (to compare with groundTruth). -Poutputdefect.off : the output mesh with the segmented defects in green.

To execute these three scripts in succession, run the command from run:

./deep-segmentation.sh PathToModel InputMesh Treshold

To reproduce the measurements in Table 1, run the following command from mesures:

./testINRAE1A.sh This command fill the results.tex file which contains the performance measure presented in this article. To reproduce the Table 2 table, run the follwing command from mesures:

./testK folds.sh After executing this command, five files are created: resultsN.tex (with N = 1...5) containing the performance measure corresponding to the lines in Table 2.

Conclusion

From the difficult problem of tree defect detection, new contributions is proposed in this work with first a new relief map able to locally adapt itself on global shape of the trunk. Such an adaption is important since the trunk geometry may appear with significant variations which make wrong the segmentation of the defect. The second contribution, is to propose a segmentation process based on the U-Net architecture allowing to outperform the previous works. The results, source code and dataset are all available on a git repository allowing the reproduction of the results together with an online demonstration.

In future works, we plan to address the defect classification. Such features will be interesting to get a finer estimation of the wood quality. Other perspectives consist in investigating the 3D point cloud semantic segmentation such as PointNet [START_REF] Qi | Pointnet: Deep learning on point sets for 3d classification and segmentation[END_REF], ConvPoint [START_REF] Boulch | Convpoint: Continuous convolutions for point cloud processing[END_REF], PointCNN [START_REF] Li | Pointcnn: Convolution on xtransformed points[END_REF], . . .

Fig. 1 :

 1 Fig. 1: Overview of the proposed method: input LiDAR 3d points (a) and its reconstructed mesh (b) are used to construct the relief map (c) which is exploited in U-Net (d). The defects are segmented and compared to ground truth (e).

Fig. 2 :

 2 Fig. 2: Examples of defects: scar and picot. Defect areas are highlighted in red.

  (a) Centerline (in red) (b) Branches are highlighted in green (c) The rest of the branches and smaller defects were detected by the patch based method

Figure 5 :

 5 Figure 5: Segmentation results of the proposed method on an European beech, which has a curve trunk and many branches.

10 (

 10 a) accumulation (c) centerline result (red) on full mesh

  (a) Centerline (in red (b) Branches are highlighted in green (b) The rest of the branches and smaller defects were detected by the patch based method

Figure 6 :

 6 Figure 6: Segmentation results of the proposed method on a partial trunk of an oak.

  With all points of the patch (b) With a subsample of the patch

Figure 7 :

 7 Figure 7: Line fitting to the profile of the distances to the centerline of a patch.

Fig. 4 :

 4 Fig. 4: Illustration of the main idea of the centerline extraction algorithm: (a) accumulation step from surface faces f k and f j in the direction of their normal vectors ( -→ n j and -→ n k ); (b) example of tracking step from the 3D accumulation values. Images (c) and (d) show the centerline extraction respectively on full and partial mesh.

  Rectangular patch PP centered on P (b) Delta distance

Fig. 5 :

 5 Fig. 5: (a) A patch in blue, associated to the red point, is used to compute the reference distance of this point. (b) Computation of the reference distance for the red point. See [18] for more details.

Fig. 6 :

 6 Fig. 6: Difference between the distance map (a) and the relief map (b).

Fig. 7 :

 7 Fig. 7: Examples of relief maps. Red color for the stronger reliefs. Blue color for the lower reliefs. The upper row is not improved with a multi-resolution analysis.

Fig. 8 :

 8 Fig. 8: Illustration of the multi-resolution analysis on an example of size 10 × 10. Ajouter description

Fig. 11 :

 11 Fig. 11: Illustration of extracting patches of size 320x320 pixels from the relief and the annotated maps.

  (a) Extracted patches with defects from relief maps (b) Extracted patches with defects from annotated maps (c) Extracted patches without defect from relief maps

Fig. 12 :

 12 Fig. 12: Some samples from training data.

Fig. 13 :

 13 Fig. 13: Illustration of input 3D points of the mixed test.

Fig. 14 :

 14 Fig. 14: Results on INRAE1a. First row: relief maps, second row: ground-truths, Third row: predictions by our network.

-

  python3 ../run/predict.py InputReliefMap PathToModel ThresholdWith InputReliefMap contains the path to the relief map.

  2 [26] andKeras [9]. During the training, data augmentation was applied randomly to the input images,

	WildServiceTree Alder1	Alder2	Alder3	Alder4	Birch2	Birch	Birch4
	Hornbeam4	Red oak2 Redoak3 Redoak4 Redoak1	Beech	Beech3	Linden

Table 1 :

 1 Comparison results: Overall row is computed from the sum of T P ,T N ,F P and F N on all the tested meshes.

	Fir1	0.747 0.769 0.757 0.137 0.937 0.238	0.746 0.857 0.797
	Fir2	0.673 0.775 0.719 0.353 0.452 0.395	0.792 0.801 0.795
	WildCherry1	0.696 0.765 0.728 0.683 0.512 0.584	0.757 0.881 0.813
	WildCherry2	0.846 0.711 0.771 0.661 0.822 0.732	0.799 0.955 0.870
	Redoak1	0.749 0.742 0.744 0.479 0.444 0.459	0.866 0.696 0.770
	Redoak2	0.428 0.833 0.564 0.061 0.400 0.104	0.730 0.428 0.538
	Beech	0.670 0.604 0.634 0.360 0.289 0.320	0.863 0.591 0.701
	Birch	0.733 0.756 0.744 0.607 0.421 0.496	0.774 0.726 0.748
	Elm	0.694 0.755 0.721 0.494 0.309 0.378	0.881 0.642 0.741
	WildServiceTree 0.247 0.741 0.370 0.057 0.463 0.100	0.856 0.504 0.633
	Overall	0.685 0.740 0.710 0.289 0.563 0.380	0.793 0.789 0.790

This research was made possible by support from the French National Research Agency, in the framework of the project WoodSeer, ANR-19-CE10-011.