Deep Learning to Segment Crossing Points in X-rays of Old Canvases

Proyecto Fin de Carrera Ingeniería de Telecomunicación Formato de Publicación de la Escuela Técnica Superior de Ingeniería Autor: F. Javier Payán Somet Tutor: Juan José Murillo Fuentes Dep. Teoría de la Señal y Comunicaciones Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2013 Trabajo Fin de Máster Ingeniería de Telecomunicación Deep Learning to Segment Crossing Points in X-rays of Old Canvases Autor: Antonio Delgado Bejarano Tutor: Juan José Murillo Fuentes Dpto. Teoría de la Señal y Comunicaciones Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2022 Trabajo Fin de Máster Ingeniería de Telecomunicación Deep Learning to Segment Crossing Points in X-rays of Old Canvases Autor: Antonio Delgado Bejarano Tutor: Juan José Murillo Fuentes Catedrático de Universidad Dpto. Teoría de la Señal y Comunicaciones Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2022 Trabajo Fin de Máster: Deep Learning to Segment Crossing Points in X-rays of Old Canvases Autor: Antonio Delgado Bejarano Tutor: Juan José Murillo Fuentes El tribunal nombrado para juzgar el trabajo arriba indicado, compuesto por los siguientes profesores: Presidente: Vocal/es: Secretario: acuerdan otorgarle la calificación de: El Secretario del Tribunal Fecha: Acknowledgements Th is document closes my stage as Master’s student at the ETSI of University of Seville. I would like to dedicate a few lines to show my gratitude to everyone who have helped me and have shared with me this incredible formative stage. I would personally like to thank my tutor Juan José for his dedication and invaluable help in all the tasks in which I have been lucky to work with him. Thank you for making me feel so comfortable working with you, for always taking into account my ideas and for everything I have learned from you in all the time we have been collaborating. Making this project would not have been possible without your help. In addition, I also thank Museo Nacional del Prado for the images provided to make this project. I would also like to thank my family and friends for their company, advice and support during these two years. The good and bad times are always better with you. Of course, I have to mention my classmates. We have shared two difficult and stressful courses in which we have always helped each other. There are many memories of coffee times, talks and meetings that I will always remember with great affection, I wish you all the best in the future. To conclude, I would like to thank the students and professors who have fought for the new format of the master’s degree that concerns us. I have been fortunate to study this new curriculum, with really interesting and current subjects. Antonio Delgado Bejarano Sevilla, 2022 I Abstract Th is document presents our research and results in the field of study of X-rays of paints. Specifically, we will present a novel method to characterize old canvases based on the estimation of density of threads in a square of one cm side. The presented method uses a Deep Learning model to segment crossing points between vertical and horizontal threads in X-ray plates of plain weaves. Then, the density of vertical and horizontal threads are calculated from the segmentation obtained with the Deep Learning model. Accordingly, throughout this document we will be talking about different aspects. Once we had described the motivation of this work in detail, we will talk about previous approaches to estimation of density of threads in which Deep Learning tools weren’t involved. Later, we will discuss about the Deep Learning model employed in full detail. The design of the architecture will be described, as well as all the variants and tools used to optimize said architecture. Likewise, all the results of the different training sessions carried out will be presented in order to select the best architectures based on different metrics. Once the most interesting models have been selected, they will be used to obtain the segmentation of the crossing points that will serve to measure the distance between threads. We will develop a series of algorithms that allow the segmented image to be binarized and, later, the count can be obtained. All these variants will be tested and compared to offer a global vision of the performance of each one. Finally, we will present an algorithm to apply everything previously developed to full paintings, so that we can obtain their density maps. A series of examples of the density maps obtained for different plates will be presented, as well as the possible matchs found between them. We will close the document with some conclusions and defining a series of future lines of research in this field. III 1 Introduction Try to learn something about everything and all about something. Thomas Huxley In this document we present our research and results in the field of study of paintings using Artificial Intelligence. In particular, we will present a novel method to characterize old canvases based on the estimation of thread densities in square patches of one cm side. The method we will present uses a Deep Learning model 1 to segment crossing points between vertical and horizontal threads in X-rays of plain weaves. Then, the vertical and horizontal thread densities will be calculated from the segmentation obtained with the Deep Learning model. Accordingly, throughout this document we will be talking about different aspects. Once we have described in detail the motivation for doing this work and some fundaments for the study of paintings, we will talk about previous approaches to thread density estimation in which Deep Learning tools weren’t involved. Later, we will discuss about the Deep Learning model employed in full detail. Once the crossing points segmentation was obtained, we will talk about how to estimate the thread density per cm. Different approaches will be developed and compared in order to determine the method that exhibits the best performance according to various metrics. The segmentation obtained using Deep Learning and the code developed to estimate the thread density allow us to get the vertical and horizontal thread density in a single crop of one cm side. Since we are interested in estimating the vertical and horizontal thread density in whole paintings, the last part of this document will describe how we have implemented a solution for this tasks making use of all methods previously described. In this way, we are able to obtain some graphs in a whole painting based on the thread density in both vertical and horizontal directions. These graphs will be named density maps. Finally, it must be taken into account that obtaining the density maps can take several hours when big paintings are processed, even using powerful machines. Since the process of obtaining the density maps of whole paintings can become actually slow, an alternative solution have been developed to take advantage of multiprocessing capabilities. Also, the performance of different machines have been studied in order to provide an order of magnitude of the processing time according to painting size and the capabilities of the machine. A summary of the tasks performed in this work is presented below: • Describe a series of concepts about fabrics, looms or canvases and the use of X-rays in the study of old paintings. • Review the scientific literature to learn about the previous works developed in the field of characterization of X-ray plates of old paintings. • Develop various Deep Learning models to segment crossing points in images of one cm side. Both the image database and the corresponding labels were provided to us by the tutor. • Compare different methods to estimate thread density from the crossing point segmentation obtained previously. • Create an algorithm to obtain density maps of whole paintings making use of some libraries provided by the tutor. 1All the code have been developed using Python. TensorFlow and Keras libraries have been used for Deep Learning tasks. 1 2Chapter 1. Introduction • Enhance said density maps algorithm using multiprocessing, testing the performance in various cases. • Generate density maps from several paintings and compare the generated density maps of different paintings to find matches between them. •Mention some future lines of research in this area. 1.1 Motivation The final goal of this work is generating various graphs per painting that allow conservators to study and compare different artworks. Among all of those graphs, we can highlight two figures: one of them represents the vertical density map and the other represents the horizontal density map. Each density map describes with a color the number of threads per cm for each region of one square centimeter of the painting. The interest in generating density maps of different paintings for analyze and study them is not new. To the best of our knowledge, until now different tools have been developed, highlighting those based on signal processing (performing a frequency analysis of the canvas) and others based on machine learning. In section 2 we will describe in detail all the previous approaches to the estimation of density maps in artworks. However, we have verified that these solutions have numerous limitations that prevent us from using those tools in different scenarios or in appropriate standards of time and ease of use. All this motivated the interest in applying tools provided by Deep Learning to try to obtain a novel method that performs the estimation of density maps quickly and accurately. In the rest of this section we will discuss about why density maps are relevant for the study of a painting. For this, we will present some key points about the fabric of the paintings, the use of X-rays in the study of said paintings and what information can conservators obtain from both vertical and horizontal density maps and the rest of graphs generated. 1.1.1 Fabric The fabric of a painting offers us a lot of useful information to study that artwork. We can see the fabric as a fingerprint of the painting, because the fabric allows the painting to be characterized based on the information obtained from the canvas. Among all the information that can be obtained from the canvas of a painting, we can highlight the following features [1]: •The type of fabric: taffeta, twill or satin •The fabric material: cotton, linen... •The weight of the fabric •The number of threads per centimeter in both vertical and horizontal orientation •The patterns of threads that can be found in the fabric •The deviations of the threads with respect to the horizontal and vertical axis All these features allow conservators to obtain a lot of relevant information about the painting. For example, the analysis of the canvas allow conservators to determine if the fabric remains intact and keeps its original size or if, on the contrary, it has undergone modifications or some pieces are missing. Also, the study of the canvas is essential to dating the artwork and it can help to determine its authorship by comparing it with other artworks with similar fabrics. It is therefore remarkable the information that the study of the canvas provides to conservators when they have to attribute hitherto anonymous artworks, or locating them in time with more precision. Thanks to the characteristics of the canvas, a painting can be located in a specific place and time and, from there, an attribution to an author can be made through these data or through a comparison with the canvas of other paintings. The study of fabrics also allows to make complete studies of authors by knowing the type of canvas they used to choose for their paintings [2]. 1.1.2 Plain weaves A canvas is a fabric made of a certain material that is used as a support for pictorial artworks. Throughout history, cotton and linen have been the most widely used materials to create canvases. Canvases has become the most used support thanks to its resistance to cold and humidity and thanks to how light and easy to transport they are. The use of canvas became popular from the Quattrocento (although there are some earlier examples), since Italian painters had high-quality Venetian fabrics that led to the change from painting on boards or walls to painting on canvas. 1.1 Motivation 3 Linen fabric is very strong and tear resistant, with a wide range of textures. However, it is expensive and can present tensioning problems with humidity. Linen is widely used for oil painting. On the other hand, cotton is more resistant to humidity changes and it is less expensive than linen. However, cotton absorbs more paint and it is more difficult to paint on it. To form the canvas, the fabric is usually adhered to a wooden chassis reinforced in the center, called wood stretcher. Thus, the fabric is taut enough to paint over on it. To prevent the fibers of the canvas from deteriorating due to contact with the oils, a primer layer is applied before painting. This prime layer is formed by a series of chemical products (glycerin, zinc oxide...) that form a smooth and clear surface on which to apply the oil. In fact, painters traditionally applied layers and layers of polished lead primer so that the surface did not have a cloth-like appearance. In this way, painters achieved almost photographic results [3]. Linen fabric [4] Cotton fabric [5] Wood stretcher [6] Figure 1.1 Types of fabrics and frame. On the other hand, regarding the elaboration of the fabrics, we can differentiate mainly three ways of binding the threads: taffeta, twill and satin. Taffeta fabrics are predominant, since they are present in a large majority of the canvases used in artworks [7]. Thanks to taffeta great strength and resistance, it is also used in the textile industry and in upholstery [8]. In addition, in this project we have only worked with paintings made with taffeta fabrics, so we will focus on describing this type of fabric although we have also mentioned twill and satin. The plain weave (taffetea) is characterized by an intertwining of the vertical threads with the horizontal ones. To make the taffeta fabric, a series of parallel threads are arranged on a loom from back to front. These threads are tightened before starting to weave and are called the warp. The separation between the warp threads follows a deterministic pattern that is determined by their placement on the loom. On the other hand, another thread is passed orthogonally to the warp from one side to the other. This other thread, which is the thread to weave, is called the weft. Since the weft follows a manual process, we can expect a certain variability that can be modeled as a Gaussian process. The spacing between the weft threads depends on how the warp threads are tensioned. This crossing of warp and weft threads is what forms the balanced taffeta fabric [9]. Both the warp and the weft are usually perfectly visible, although there are weaves where the threads of one direction can hide the threads of the other direction. Figure 1.2 Plain weave - Taffeta [10]. 4Chapter 1. Introduction When studying a canvas, both the type of textile fiber used and the relationship between the weft and the warp give us a lot of information. Because of this, the count of threads per centimeter in both directions (vertical and horizontal) is used as a fundamental characteristic in the study of canvases [1]. We are interested in knowing (and that is the final goal of this work) the number of vertical and horizontal threads per square centimeter in each fragment of the canvas of that area. When making comparatives between different paintings, data on thread densities can be very useful. For example, we know that all canvases made from the same warp roll will have the same thread densities in the warp direction. Logically, counting threads in squares one centimeter on a side is a very complex and tedious task to do manually. In fact, we can face situations where different people observe a different thread number, or that some threads are not complete on the edges and it is not known whether to count that edge thread or not. In addition, when working with old canvases we have to take into account the variability of densities in a canvas as a result of the use of very old looms and the passage of time. Figure 1.3 Parts of an old loom [11]. As we mentioned before, there are more study parameters apart from the density of threads per centimeter, such as the deviation suffered by the threads with respect to the horizontal and vertical axes. This deviation occurs because the canvas is fixed to the frame by nails, so we can study the deviation to get information about the canvas. The deviation of the threads and the position of the nails can be useful for conservators to study the integrity of the canvas or to compare productions between different artesanal workshops. Finally, it remains for us to talk about how to count the threads of the canvas, since this is not a trivial matter. The front of the canvas (where the image is painted) is covered with paint and primer. For this reason, distinguishing the threads is virtually impossible, even with highly detailed photographs. We could take a photograph of the back of the painting to study the canvas, however this is often not possible because new pieces of cloth are usually added to the back to reinforce the canvas over the years, making it impossible to photograph the original fabric. For all these reasons, the methodology currently adopted (and with which we have worked) consists of analyzing X-ray plates of the pictures [12]. 1.1.3 X-rays As we have mentioned, the use of X-rays to study the fabrics of the paintings is fundamental, since it allows to access to images of the threads that make up the canvas in order to count them. This is possible thanks to the primer that was applied to the canvas to paint over, since as we have commented, said primer has elements such as lead that are opaque to X-rays, allowing the radiographic image to be formed. However, the X-ray plates superimpose all the structures. So, we will see in the same image the threads, the paint, the cracks, 1.1 Motivation 5 the wood stretcher, the damage, the differences in primer between areas of the same painting, etc. All these elements will be treated as noise, since we are only interested in being able to count threads. Finally, we can summarize this chapter in a serie of key ideas: • The fabric of a painting gives us a lot of useful information to study, compare and preserve that painting. • The main feature that we can extract from the fabric of a painting is the thread density per squared cm in both vertical and horizontal directions. •X-rays plates of paintings are used to be able to count the threads of a fabric. Adan y Eva (Rubens, 1628) [13] X-ray plate of Adan y Eva [38] Figure 1.4 Painting and X-ray plate of Adan y Eva by Rubens. Figure 1.5 Vertical density map of Adan y Eva by Rubens. 2 State of art The more you read, the more things you will know. The more you know, the further you will go. Dr Seuss In this chapter we will present the review of the literature carried out to know and compare the solutions provided by other authors when studying and characterizing paintings. Special attention will be focused to those works related to obtaining and analyzing the thread densities of the fabric. Besides, other works that are not so related to the objective of this project will be mentioned, as they can provide certain ideas for future and innovative lines of research. Among all the solutions that we will mention, we will highlight those based on frequency studies, robust but with some limitations, as well as new approaches through Machine Learning and Deep Learning, since our work is based on the design and use of a Deep Learning architecture to obtain the density maps of the paintings. 2.1 Fabric frequency analysis In the previous chapter we have talked about the relevance of studying the canvas when it comes to restoring, analyzing, or preserving the paintings. In this context, the great amount of information that thread density maps can offer us to date and compare artworks has become clear. In the literature we find different works focused on obtaining the density maps. In this first section we will focus on several works that are based on a frequency analysis of the fabric, obtaining good results in many cases. 2.1.1 Using Fourier Transform We find in [14] the first work where a theoretical framework is defined in order to model a fabric through a frequency analysis based on the Fourier Transform (FT). In said work it is commented that we can know a series of parameters that describe a fabric from its Fourier Transform, that is, from the FT of the X-ray plate. A few years later, D.H. Johnson (one of the leading researchers in this field) and other authors publish two papers where they put into practice the aforementioned theoretical framework, thus applying FT to study canvases. In [15], the first of the two works that we mention, the way of obtaining the thread counting maps is briefly presented thanks to the study of the frequency spectrum. This spectrum is obtained from images of 1 cm side cropped from X-ray plates of the fabric. Likewise, a method to compare the obtained density maps and match them is presented, so that it can be concluded if two canvases were obtained from the same loom, with all the consequences that this information has when dating or attributing works. In [16], published three years later, the previous work is extended providing much more information and specifying for various types of fabrics, although the core is the same: obtaining density maps from the frequency spectrum and the subsequent comparison to be able to match them. 7 8Chapter 2. State of art Bearing in mind that we are dealing with an X-ray plate of a painting, the intensity of each pixel depends on the chemical and structural composition of that area of the image, that is, the amount of paint and primer (thickness of the layers), the chemical composition, if that region is overlapping the wood stretcher, if there are nails nearby, etc. The authors also describe how the fabric can be seen thanks to the lead primer. This compound stays in the gaps between the threads and it is seen as a light gray on the plate, thus allowing the structure or pattern of the fabric to be seen (except when the image is saturated, there you can’t see the underlying structures like the thread pattern). The greater the amount of paint and primer, the more X-ray energy is absorbed and therefore the clearer that region (gap between threads) appears on the plate. The following image shows the behavior described in a plain weave fabric: Figure 2.1 X-ray of P1114 Ribera (MNP) [34]. The reason why the authors resort to frequency studies of the X-ray plate is that the threads form a repeating pattern on the plate, as can be seen in the image above. As explained in [16], this repetition pattern can be described as a quasi-periodic function in both vertical and horizontal dimensions. For all these reasons, in the paper it is proposed to use the discrete FT in two dimensions (2D-DFT) of 1 cm square patches to obtain the thread counting maps. In addition, they propose using an overlap of 0.5 cm between two adjacent patches to improve the resolution and to get a better appearance of the resulting maps. Therefore, in order to count the threads automatically, 2D DFT of images with a size of 1cm is used. [16] can be consulted to learn more details of the calculation and the expressions of the FT of a patch. However, we will expose the aspects related to the interpretation of the spectrum, since to obtain the spectrum we will resort to Python libraries that do the calculation autonomously. Let’s see an example of a spectrum obtained from a 1cm cut of a plate of our collection: Figure 2.2 Spectrum of a 1cm patch from a canvas. 2.1 Fabric frequency analysis 9 We can appreciate in the previous image a symmetry in quadrants 2 and 4, as well as in quadrants 1 and 3. This is due to the input image takes real values (the value of the pixels will be in the range 0-1 or 0-255), so the FT will be conjugate symmetric according to the properties of the transform [26]. That is why in the previous representation, where we see the magnitude of the spectrum, two pairs of identical quadrants appear. Another aspect to take into account is how to obtain the thread count from the representation of the spectrum. The location of the maximums that appear on the vertical and horizontal axes give us the value of said count, since the peaks of the magnitude of the FT are due to the periods found in the image. In this case, those periods correspond to the frequency of repetition of the threads in both directions. The maximum of the vertical axis will give us the number of horizontal threads and the maximum of the horizontal axis will give us the number of vertical threads per centimeter (because the image from which we obtained the spectrum has a side of 1cm). This is the simplest approximation, where maxima are located on the axes themselves. However, it may happen that in the image to be analyzed the threads are inclined, normally due to the cusping phenomenon. This phenomenon occurs when placing the fabric in the wood stretcher (fixing it with nails) before applying the primer, which causes the threads to move under tension and the usual pattern is distorted. In [16], rather than treating this as a problem, it is presented as an opportunity to extract more information from the fabric, since the spectrum allows us to extract the angle of inclination of each patch of the fabric. Therefore, to take into account the possible rotation when searching for the maxima on the axes that give us the thread count, it is proposed to perform the search in a wedge-shaped area. This is equivalent to performing a bandpass filtering of the spectrum, removing all the information from the image that is not related to the repetition pattern of the threads and their inclination. The paint, the nails, the flaws, or the wood stretcher are treated as interferences that are not desired in the frequency analysis of the quasi-periodic functions. This also eliminates the low-frequency components associated with the gray tones of the image. Definitely, to find both thread densities and the angle of rotation (if it exists), all you have to do is define the parameters that give us the wedge region to study: the minor radius, the major radius and the opening angle that define that circular sector or wedge. Each fabric requires a specific study region, so these parameters are set through a user interface in the work presented in [16]. A special case can happens if there is more than one maximum in one of the axes. The authors solve this event by choosing the maximum that gives the closest count to the count of contiguous patches. Once the count of each 1cm region is obtained, density maps are generated. This maps represent the number of threads in each region with different colors. As already mentioned, the comparison of these maps is what allows us to know if two paintings were obtained from the same roll. Finally, note that [16] also presents a method for comparing density maps beyond visual appearance, although we will not describe it because it is not within the scope of our work (we are only interested in obtaining density maps). I do think it is interesting to comment that, as described in [16], it is necessary to know if the warp corresponds to the vertical or horizontal orientation of the canvas, since the placement of the fabric on the wood stretcher could have been done in one orientation or another. It may be the case that two frames cut from the same roll have density maps that coincide when one of them is rotated 90º, precisely for this reason (we will see some examples throughout the work). 2.1.2 Power Spectral Density Another interesting work is the one presented in [17], whose authors are researchers and professors at University of Seville and also have extensive experience in this field of study. In said work, the use of FT is recognized as a tool to obtain the density maps of a patch (as we saw in previous works). However, the authors describe some limitations found when using FT to obtain density maps as well as the usefulness that these maps can have when making comparisons between paintings. Mainly, it is argued that the frequency study of fabrics based on FT does not work well when it is applied to paintings whose fabrics are damaged. In fact, it is not uncommon to find old paintings in which a piece of the original canvas is missing or irregularities have appeared as a result of the passage of time or other vicissitudes. In addition, they find that FT also does not work well if the primer is thicker than normal, since it does not allow the pattern of the threads to be distinguished. In summary, the pattern of repetition of the threads must be about distinguishable to some degree for the FT to offer the expected results. 16 Chapter 3. Design of a Deep Learning Model for Thread Counting Figure 3.1 X-ray detail of P001692 Rubens (MNP) [38]. On the other hand, the other scenario where FT has been found to work poorly is in some fabrics where the threads (usually in the horizontal direction) are not well perceived. Instead, a small nodule can simply be seen where the vertical threads cross the horizontal ones. This usually happens because the threads in one direction are quite tight, causing that threads in perpendicular direction can’t be seen. In this case, since the threads of one of the directions cannot be seen clearly, the FT is not capable of analyzing the repetition pattern in said direction. This is the case that has caused the most problems so far, since the results obtained with Aracne [18] in terms of thread count were totally useless. This is a usual situation in paintings by Velazquez. We present an example of a patch taken from the radiographic plate of Prince Baltasar Carlos on Horseback by Velázquez (reference P001180 in MNP) [39]: Figure 3.2 X-ray detail of P001180 Velázquez (MNP) [39]. These problematic cases are the ones that motivate to look for a new solution away from the frequency analysis used until now. For this reason, we will test a Deep Learning alternative with a view to improving the results obtained with frequency analysis. The segmentation of the crossing points would allow, in principle, to solve the problem found in scenarios such as those described in Figures 3.1 and 3.2 and it must be useful for the rest of the cases. 3.2 Deep Learning basis 17 3.2 Deep Learning basis Before starting to describe both the dataset used and the model designed, I would like to present some necessary notions about neural networks, since we will refer to some specific aspects throughout the chapter. The core of Deep Learning are artificial neural networks, known by their acronym ANN. These networks base their operation on the same idea that we find in biological neurons and their interconnection. ANN are tremendously powerful and usable in a huge spectrum of use cases: from automatic translators that follow natural language processing models to recommendation systems for streaming platforms, for example. Since the 90s of the 20th century, ANN use has spread at an enormous speed, thanks to improvements in training algorithms and the increase in computing power of machines, especially thanks to the appearance of powerful GPUs. Their main problem is that they require huge databases to be trained correctly in order to make accurate predictions [40]. As we have said, the basic processing unit that forms a neural network is called neuron due to the biological simile. Neurons performs a simple operation: they receive a series of inputs and provide as output the weighted sum of these inputs. To determine with what intensity each input affects the weighted sum there are some parameters called weights ( w ), which must be adjusted during training. There is also another parameter, called bias ( b ), which represents a constant value that is added to the product of the weights by the inputs [43]. Figure 3.3 Neuron [43]. The output of the neuron, y , therefore corresponds to a linear weighting of the inputs. However, to make the network capable of finding complex patterns (which usually involve non-linearities), we will have to transform the output of the neurons in some way so that they are non-linear, since otherwise we would not achieve anything by grouping many neurons (because several successive linear combinations can be described as a single linear combination). This problem is solved using the activation function, which receives the output of a neuron as input and modifies it according to a specific non-linearity. Among the most outstanding we find the activation functions ReLU, sigmoid or hyperbolic tangent [40]. Once the output of each one of the neurons is non-linearized, it is possible to chain many neurons (or layers of neurons) to make the network distinguish very complex patterns to classify, group data or predict new values [43]. Bias has been omitted to simplify the representation of Figure 3.4, but following the scheme of Figure 3.3 the term b would be added to the product of the weights by the input and, subsequently, all that sum would go through the activation function: Figure 3.4 Activation function [43]. 18 Chapter 3. Design of a Deep Learning Model for Thread Counting ANNs are networks formed by many of these neurons. We can place many neurons in a layer, or many layers with many neurons forming deep networks. It is important to mention that all neurons in a layer are connected to all neurons in the next layer. Because of this, these networks are known as fully connected or feed-forward networks [40]. The input of a neural network will have as many neurons as input data dimensions. For example, if the input is a vector of 100 elements, the first layer of the network must have 100 neurons. The number of output neurons will be given by the application for which the network was built, highlighting in this case the classification or regression networks. In regression networks there is a single output neuron, which provides the value of the regression. In classification networks there are as many output neurons as classes to be classified, activating one of the neurons above the rest depending on what class the input is. Figure 3.5 Feed-forward Network [43]. Despite all the benefits discussed about fully connected ANNs, they have a limitation when working with images. To handle images as input data, images are first converted into a one-dimensional vector with as many inputs as there are pixels in the image. As we have commented, each input data corresponds to a neuron of the first layer, so when working with images there will be an input neuron for each pixel of the image. Later, when more layers with more neurons are added, since they are all connected to each other, the number of parameters would increase considerably when working with images of a certain size. On the other hand, having converted the image into a vector, all the spatial information of the image has been lost. This causes that we do not have access to many spatial features present in an image which can be extracted from the location of each pixel. For example, the relationship of each pixel with its neighbors or the structures that groups of pixels defines can give us a lot of useful information [44]. For all this, convolutional neural networks (CNN) were developed and are the ones used when working with images. Consequently, the type of network that we are going to design for our problem will be a CNN. This type of networks redefines the approach to be able to extract spatial features from images, one of the limitations of fully connected ANNs. Now, the hidden layers are not formed by neurons as we saw before, but by convolutional filters (also called kernels) [40]. These filters are convolved over an array of pixels and they are shifted across the image so that all pixels in the image have convolved with the filters. By convolving the kernels with the group of pixels, some features from that group of pixels are extracted (lines, shapes...). So, when passing the image through the filter we have a map of the image’s features (one map for each filter). Each input of these kernels takes a value that weights the convolution. These weights will be trained in the CNNs (analogous to the parameters of the neurons in fully connected ANNs). As in fully connected ANNs, the result of the kernel convolution is passed through an activation function, the ReLU function being the most typical in CNNs. In addition to convolution, CNNs have another distinctive element: subsampling [40]. Only a fraction of all the features collected by the filters will be preserved (those that contain the most relevant information). To choose the characteristics, a selection criterion is taken: normally the maximum of a group of characteristics is the criterion that is usually used, although there are others such as the average value of that group of characteristics. Subsampling allows to reduce the size of feature maps, which means that the number of network parameters does not grow exponentially out of control as it does when using fully connected ANNs to images. This process of convolution, activation and subsampling is repeated in all the layers, managing to transform an image into a feature vector with dimensions that depend on the design of the architecture. This vector does contain spatial information of the image, encoding said information in fewer entries. In addition, CNNs layers have a very well defined hierarchy: the first layers are 3.2 Deep Learning basis 19 in charge of looking for simple spatial information: straight lines, curves, simple shapes... while the deeper layers detect much more complex shapes or patterns that are based on the characteristics of the first layers. Therefore, the deeper a CNN is, it will be able to detect, a priori, more complex patterns of information within the image. Once features have been extracted from an image, they can serve as input to a fully connected ANN for regression or classification. This structure of CNN to extract features + feed-forward network for classification or regression is very common [43]. Figure 3.6 CNN + Fully-connected ANN [43]. Figure 3.7 Subsampling 2x2 (MaxPool) in CNN [43]. Apart from the types of networks mentioned, some more complex architectures can be found, such as recurrent neural networks (RNN) [50] or transformers and attention mechanisms [51], although we will not describe them because they are beyond the scope of this document. It remains for us to mention how is the learning of a neural network. We are working in the supervised learning paradigm, that is, we have the labels (or ground truth) of the data that we use to train the network. Therefore, network output during training is compared with said reference truth. This comparison serves to measure the error between the expected output and that provided by the network. This error is called cost function or loss funcion, and there are several types depending on the task that we want the network to perform: minimum squared error (MSE), Kullback-Leibler KL divergence, Binary-Crossentropy, etc [40]. The value of cost function will depend on the output value of the network, which depends on the value that the parameters of the different neurons have taken. Therefore, once the error has been calculated for certain values of the parameters, we are interested in knowing how to change these parameters so that the error decreases. This is known as the gradient descent method [40]. The cost function will usually be a function that is neither concave nor convex, so a priori it will not have a global minimum but there will be some local minima. The strategy consists of calculating the partial derivative of the cost function with respect to each of the parameters, thus obtaining the gradient vector of said function. The gradient gives us the direction of maximum growth of said function, so we will move in the opposite direction to the gradient to go in the direction of maximum decrease. Thus, the value of the cost function decreases and, therefore, the error between the output of the network and the ground truth also decreases. Initially, the network parameters take random values (unless they are explicitly initialized to a certain value). The cost function can be represented as a multidimensional surface (as many dimensions as there are parameters in the network), so that by varying the value of the parameters we move along said error 20 Chapter 3. Design of a Deep Learning Model for Thread Counting surface. The gradient makes us move locally along said surface looking for areas where the function (the error) is smaller. In the following figure we can see a representation of how gradient descent works: the surface represents the cost function and the black arrow indicates the zone of maximum descent (given by the opposite direction to the gradient). As we move to a lower area of the surface, the parameters θ1 and θ2 (on the axes that form the plane under the surface) update their value [43]. Figure 3.8 Performance of gradient descent on a surface representing the error [43]. Once the gradient is calculated, the new values of the parameters, ˆw , are obtained by subtracting the value of the calculated gradient from the previous value of the parameters, w. This new set of network parameters will cause the network output to be different, therefore the value of the cost function will be different. This process is iterative and it is repeated until we have reached a local minimum. This process is also influenced by the learning rate [40], η , a hyperparameter of the network. The learning rate weights how much the gradient affects the update of the parameters in each iteration, that is, how much progress is made in each step. ˆw=w−η·∇f(3.1) η is a very important and influential hyperparameter in the training of a model, since if the learning rate is very high the learning will be practically random, jumping from one area of the surface to another without converging. On the other hand, if the learning rate is very low, it is possible to get stuck in local minima of little relevance. This gradient descent method has many more complex and sophisticated variants that optimize the algorithm. The most used are stochastic gradient descent (SGD), Adam, Adagrad or RMSProp [40]. In our work we have chosen the Adam optimizer, since it tends to be of fast convergence, computationally efficient and with little memory consumption, usable for networks with many parameters or input data. In addition, Adam is cataloged as the best optimizer for practically all cases in [45]. It only remains to mention how to calculate those partial derivatives with respect to each parameter of the network, necessary to build the gradient. It is not a viable option to do it by brute force, that is, by calculating each and every one of said partial derivatives, since a network has many neurons, each neuron several parameters and many possible paths (remember that all neurons are connected to all neurons of the next layer). Instead of calculating all this, the backpropagation algorithm [43] is used. The idea of backpropagation is to assign part of the error to each parameter. To analyze this distribution of responsibilities in the error, it is better to follow a back-to-front strategy, that is, from the last to the first layers. This allows backpropagation of the error and also allows working efficiently, reusing the derivatives already calculated [43]. Backpropagation algorithm begins by calculating the derivative of the cost function with respect to the parameters of the neurons in the last layer. To calculate the derivative, you have to see the path that connects each parameter (of each neuron) with the final cost. In the last layer this is relatively simple: first find the weighted sum that the neuron outputs ( z ) and then pass this weighted sum through the activation function ( a ). The output value of the activation function of the neurons of the last layer is precisely the output value of the 3.2 Deep Learning basis 21 network. Therefore, comparing the value of the activation of the neurons of the last layer with the ground truth, the value of the cost function is obtained. Therefore, following the chain rule, the derivative of the cost function (C) with respect to the parameters (w) of the neurons of the last (L) layer can be decomposed as: • The derivative of the cost function with respect to activation, that is, the derivative of the cost function itself (MSE for example). • The derivative of each activation function with respect to the weighted sum of each neuron, that is, the derivative of the activation function itself (sigmoid, hyperbolic tangent...). • The derivative of the weighted sum with respect to the inputs received by the neurons of the last layer, ∂C ∂ωL=∂C ∂aL·∂aL ∂zL·∂zL ∂ωL(3.2) In this way we have already calculated the partial derivatives that allow us to see how the cost function varies with respect to the parameters of the neurons in the last ( L ) layer. If we propagate this error to the previous layer ( L−1 ), the derivative of the cost function with respect to all parameters can be broken down by the chain rule into: • The derivative of the cost function with respect to the activation of the last layer, that is, the derivative of the cost function itself (MSE for example). • The derivative of each activation function with respect to the weighted sum of each neuron in the last layer, that is, the derivative of the activation function itself (sigmoid, hyperbolic tangent...) of the neurons in the last layer. • The derivative of each activation function with respect to the weighted sum of each neuron of the penultimate layer. • The derivative of the weighted sum with respect to the input received by the neurons of the penultimate layer, which corresponds to the value of the output (activation) of the neurons of the previous layer (L−2). ∂C ∂ωL−1=∂C ∂aL·∂aL ∂zL·∂zL ∂aL−1·∂aL−1 ∂zL−1·∂zL−1 ∂ωL−1(3.3) As we walk the network backwards, the number of partial derivatives to be calculated grows ostensibly. However, thanks to the backpropagation algorithm, most of these partial derivatives are already calculated. As we see in the equation of the L−1 layer, the only new thing to be calculated is the second term, that is, how the weighted sum of the L layer varies as a function of the output of the L−1 layer. In this way we move the error from one layer to the previous one, since the only layer that has access to the cost function in its equations is the last one. This is repeated iteratively until the first of the layers, so that we have all the partial derivatives we need just calculated [43]. The training scheme of a network would therefore be as follows: the data that enters the network passes through all the neurons, which will show a certain output according to the value of their parameters (forward propagation). Afterwards, the output of the last layer will be compared with the ground truth, measuring the error according to a cost function. Finally, the error is propagated backwards in order to update all the network parameters so that the error between the network output and the ground truth decreases [43]. 22 Chapter 3. Design of a Deep Learning Model for Thread Counting Figure 3.9 Training scheme of a network [43]. 3.3 Provided dataset Before describing the model, we are going to dedicate a subsection to talk about the images provided to train it. Although the generation of the images and their labeling for training are not part of this work, I consider it relevant to briefly comment on how the dataset with which we are going to work is structured. The set of images has been extracted using 41 X-ray plates from different paintings provided by the Museo Nacional del Prado. The paintings that we can find in the dataset are very different, representing several centuries (from the 17th to the 19th), different authors and many different types of fabrics. Some of the most outstanding authors present in the set of images are Diego Velázquez, Pedro Pablo Rubens, Francisco de Goya, El Greco, Herman van Swanevelt, Miguel de Pret, Luis Paret, Bartolomé Esteban Murillo, Claudio de Lorena, Jean Lemaire, Nicolas Poussin, Gaspard Dughet, Miguel Cabrera, Jan Both or José de Ribera. Numerous patches were obtained from all these radiographic plates, of which 260 labeled patches were provided. All of them were preprocessed before labeling to improve their visualization. Furthermore, they all are 300x300 pixels size and everyone has a resolution of 200 pixels per centimeter. So, each patch has physical dimensions of 1.5cm on each side. The resolution of each plate was adjusted during preprocessing, since originally each radiograph was digitized with a specific resolution. Patch extracted from a painting Labeling of that patch Figure 3.10 A patch and its label. 3.3 Provided dataset 23 Also highlight the heterogeneity of the dataset, since we have images of paintings by many different authors and of different types of pain weaves canvas. In addition, the thread densities per centimeter that we found in the dataset range from 6 to 24 threads per centimeter. This shows that we are going to work with a wide variety of cases, which may be representative of the cases that we would find in a study of all the paintings in a museum with taffeta weave (in terms of thread densities). On the other hand, the dataset contains plates that are better or worse preserved: some with a lot of noise (cracks, damage...) and others with a very good appearance, in which the threads are easily distinguished. For all these reasons, the dataset used is very rich in the number of different paintings (more than 40), in the number of authors represented, in the quality of the images and in the range of thread densities. Figure 3.11 Example of the variety present in the dataset. 3.3.1 Working with the dataset Once the provided dataset was described, it was necessary to adapt some aspects to be able to use the dataset in the training of the Deep Learning model that we will propose. More specifically, there are two common tasks in Deep Learning that are usually made to train a model: the first is separating the available dataset, so that some images are used to train and validate the model and others to check how well it works making predictions with other data unpublished for it [41]. The other common task is to implement a Data Augmentation function which will perform a series of transformations on the images to increase the number of available samples in order to obtain better results [40]. Below we will describe in more detail how we have dealt with these two aspects. Separation of the data set into subsets In Deep Learning, it is common to work with different subsets when training and evaluating the performance of a model. The most common strategy is to divide the dataset into three subsets: a training subset, a validation subset, and a test subset. The usual recommendations that are usually followed argue that it is best to allocate 70% of the data for training, 20% of the data for validation and 10% of the data for testing [40]. However, these proportions are only indicative and may vary from case to case to the programmer’s liking. Actually, it is not even mandatory to have all three sets, although it is recommended. Going a little deeper, the training set (the largest) contains the data that the model will use to fix the parameters of the neurons. This error is not representative of the performance of the model, since the parameters have been adjusted to have a low error for these specific data, so the error in training has a bias that means that the error in real situations is underestimated. [41]. In fact, if the model is badly chosen or if the model is trained for too many iterations, the training error can tend to zero, since during learning the network can overfit its parameters for that specific data set obtaining an error that is practically non-existent. This happens because the network has enough parameters (degrees of freedom) to fix them so that the training error is very close to zero, but this does not mean that the model works well, since the model will not be able to generalize to other data for which it has not been trained. This behavior is called overfitting, and it is something to avoid [47]. 24 Chapter 3. Design of a Deep Learning Model for Thread Counting On the other hand, the validation set is used to adjust the hyperparameters of the model, such as the learning rate, the optimizer, the number of epochs or the batch size. Both the learning rate, η , and the optimizer have already been described previously. For its part, the batch size indicates how many training samples enter the network simultaneously [40]. Normally not all the training data is passed at once, as this would be computationally very expensive. Instead, groups of data are passed consecutively, so that for each group forward propagation is performed, the error is calculated and backpropagation is performed to update the parameters [46]. The number of epochs indicates how many times all the data will pass through the network. The validation data has not been used for training, so the model is impartial in providing an output when it receives this data as input. However, the error in validation is also underestimated, since the hyperparameters of the model have been chosen based on these validation data [41]. The test set is used to obtain a representative error of the model’s performance. Test data is only used to see how the model performs with out-of-sample data, that is, data that is completely new and that the model has not seen before. Therefore, the test error is unbiased and describes how the model will behave in a real production situation. In the case of our work, the division of the data into the three subsets has not been trivial. Due to the great heterogeneity of our data, it was necessary to distribute the images carefully in the three sets, so that in all the sets we find the different densities of threads per centimeter well represented, as well as images with better or worse quality. Also, since we have multiple labeled patches from each painting, care has been taken that all patches from each painting are in the same set. This avoids the case that there are patches from the same painting in the test set and in the training set, since if this happens we would not be being entirely rigorous when evaluating the performance of the network. Data augmentation We have also implemented a Data Augmentation function to increase the number of available samples. To get more samples, the available images are artificially modified, generating variants of each training and validation image in this case. These modifications must be realistic, not random, since the idea is to reproduce variations that can appear as real circumstances. The most typical are rotations, displacements, orientation changes, lighting modifications, vertical or horizontal flips, etc. In this way, many more can be drawn from each sample, considerably increasing the size of the dataset. We have implemented a Data Augmentation function with which we obtain the following samples: • Since the images have a size of 300x300 pixels and a resolution of 200 pixels per cm, four patches will be taken from the four corners with a size of 200x200 pixels, each of these patches representing a piece of fabric with a side of 1cm (for that resolution). This generates 4 samples for each image. Original 300x300 pixel image Getting the 4 corner 200x200 pixel patches Figure 3.12 Obtaining 4 patches from the corners of the original image. • The original 300x300 pixel image is flipped from right to left, then obtaining the four 200x200 pixel patches from each corner. This generates another four samples. 3.3 Provided dataset 25 Figure 3.13 300x300 pixel image flipped from right to left. • The original 300x300 pixel image is flipped top to bottom, then obtaining the four 200x200 pixel patches in each corner. This again generates four new samples. Figure 3.14 300x300 pixel image flipped top to bottom. • From each of the three 300x300 pixel images (the original, the one flipped from right to left and the one flipped from top to bottom) two centered patches of size 200x200 pixels are extracted. The first one covers the region defined by the pixels [50:250, 50:250] and the second covers the region defined by the pixels [65:265, 65:265]. These patches (two per image) are randomly rotated between 2ºand 7º with respect to the vertical, in both directions. Therefore, in total we add 12 new samples obtained from rotating these central patches (four for each of the 300x300 pixel images). • Additionally, the patches in the region [50:250, 50:250] are also rotated with a random angle between 8ºand 12º, in both directions with respect to the vertical. Half the images are obtained than in the previous case because in this case a more aggressive rotation is applied. We want this aggresssive rotation to be less represented in the dataset. By this way 6 new samples are generated (two for each 300x300 pixel image) 32 Chapter 3. Design of a Deep Learning Model for Thread Counting Generative Adversarial Networks (GANs) are also mentioned [61]. This type of network consists of two parts: a generating network and a discriminating network. Broadly speaking, the generating network provides samples similar to the real ones, with the discriminating network being responsible for classifying each generated sample as true or false. GANs can be useful in segmentation tasks in order to increase the number of samples available to train a segmentation network. For example, a GAN can be used to generate new synthesized images and their segmentation masks. Also comment on the existence of other techniques such as: • Models using dilated convolutions [59]. These convolutions are obtained by leaving gaps to increase the receptive field of the convolution. In this way, larger areas are covered without increasing the number of network parameters. • Regional CNNs (R-CNN) [60], which use a region proposal network (RPN) that provides the estimated bounding boxes of the candidate objects to be segmented. • Current models based on attention [51]. Highlight the Vision Transformers [58], tremendously powerful but that need huge amounts of images to be trained. 3.4.2 Development of an architecture for segmentation from scratch Among all the architectures and options described above, we decided to develop our own version of a U-Net for the segmentation of the crossing points. This architecture usually works well with datasets that are not too large (as long as they are accompanied by good Data Augmentation). In addition, there is a lot of code available on the internet that makes it easy to develop our own version of the U-Net with the Keras and TensorFlow. Specifically, we use the code in [62] as a starting point. Highlight from said code the definition of the function conv2dblock. This function defines two consecutive convolutional layers, so that the output of the first layer is the input of the other convolutional layer. Both the number of channels or filters and the size of the kernel can be passed as parameters. In addition, it offers the possibility of applying batch normalization after both convolutions. The batch normalization method [40] is an extra layer that serves to standardize the data after passing through a layer or series of layers, making the training much more stable and reducing the number of epochs needed. The architecture follows the pattern we described for the U-Net in the previous section: conv2dblock + MaxPooling2D on the contracting path and Conv2DTranspose [63], concatenations, and conv2dblock on the expansive path. In addition, after each Pooling operation (MaxPool2D) and each oversampling operation (Conv2DTranspose) a Dropout layer can be applied. Dropout layers [64] act as regularizers to avoid overfitting. The idea is to deactivate a series of neurons randomly, thus preventing the network from overfitting. This way of making the training noisier is comparable to training different architectures in parallel. The percentage of neurons that are deactivated is a hyperparameter of the network that has to be set. With all this we can now define our own U-Net architecture. We start by defining a five layers U-Net. All convolutions use kernels of size 3x3 and are of type "same" [65], that is, a certain padding is applied when doing the convolution so that the output is the same size as the input. The number of filters in the first layer is passed as a parameter when instantiating the model and the number of filters is multiplied by two from layer to layer as we go down the network. First, we set the number of filters in the first layer to 16, doubling in subsequent layers. The Dropout value on each layer will be 0.2 initially. All layers use the ReLU activation function except the last layer. The activation function chosen for the last layer was a sigmoid, so that each segmentation pixel takes a value between 0 and 1. As input we will have tensors of size (N,200,200,1), where N is the number of images that pass through the network simultaneously (batch size). Figure 3.23 Activation functions used [66]. 3.4 Creating a neural network 33 Figure 3.24 conv2dblock [62]. In the following diagram we can see the architecture of the network and how the size of the feature maps varies from input to output. Note that in order to concatenate in the first layer of the expansive path, zero padding had to be done to go from size 24x24 to 25x25 at the first step in the expansive path: Figure 3.25 Architecture used with conv2dblocks modules in blue horizontal arrows, with the number of filters used, ni , indicated on top. Up and down steps involve halving reductions of sizes downwards, and duplicate them upwards, with zero padding if needed. Upcoming features are concatenated with the features at the same level in the encoder. The size of all feature maps is indicated below them. At the output a 1×1convolutional filter with sigmoid function is used. U-Net training For the architecture training we will follow these guidelines: • 5 different trains will be made. Each of them will have two random characteristics that will cause some variability between them: random initialization of the network weights and random choice (within the fixed ranges) of the angles used for the Data Augmentation rotations. • Two Keras callbacks will be used: Model Checkpoint and Early Stopping [67]. Callbacks allow to perform actions during network training. The first one will be configured to save the weights with lowest validation loss. The second one will be configured to stop training if the validation loss have not improved in the last 5 epochs. • We will use Adam as the optimizer, a learning rate of 1e-3, a batch size of 32, and a maximum number of 100 epochs. The shuffle option will be activated so that the images are shuffled when training. In this way, each batch that passes through the network will contain images from different plates. This also helps to prevent overfitting. • For the training we will use the binary cross-entropy function as loss function [68] and the accuracy as metric to evaluate the model [70]. 34 Chapter 3. Design of a Deep Learning Model for Thread Counting Binary cross-entropy function is governed by the following equation: BCE =−1 N N ∑ i=1 yi·log(p(yi))+(1−yi)·log(1−p(yi)) (3.4) In the above equation, yrepresents the value of the label (1 for the crossing point pixels and 0 for the rest). p(y) is the probability that a pixel is a crossing point. It will take a value between 0 and 1, closer to one the more certain it is that that pixel belongs to a crossing point. The deep learning model will predict the probability that each pixel is a crossing point and the value of the loss function is calculated by checking the label. Let’s take the example of a pixel for which the network output predicts with a very high probability that it is a crossing point and taking the label of that pixel value 1. In that case the loss function increases by value - log(p(y)). The logarithm of a number close to 1 will be a negative value but close to zero, so when the value of the logarithm changes sign, the loss function will increase very little. On the other hand, if the probability of being a crossing point is very small, we will have that the logarithm of p(y) will be a high negative number, so when the sign is changed, the loss function will increase a lot. This behavior is analogous for pixels in which the label takes a value of zero. In this way, this loss function takes a high value for bad predictions and a low value for accurate predictions. This is computed and averaged for all the pixels. Figure 3.26 Cross-entropy or log loss [69]. As for the accuracy metric, its formulation is as follows: Accuracy =Number of correct predictions Total of predictions (3.5) The total of predictions is the total of the 200x200 pixels of each image. This metric compares the binary value corresponding to the output that the network has provided in each pixel with the value of that pixel in the label. This means that it binarizes the output image, since the network output takes values between [0-1] and it must be binarized. The threshold that the function uses internally is 0.5, half the range of pixels values. With all the above, we perform the first training of the model with the described hyperparameters and callbacks. We will call this first model getUnet, and it has about two million parameters. All the trainings have been performed on an Intel Xeon E5-2630 v4 processor with 40 CPU Cores and two GPU NVIDIA Tesla P100 16GB. We present below the results obtained in the 5 training sessions: Table 3.1 getUnet training results. getUnet training results Train Number Epochs Train Loss Train Accuracy Valid Loss Valid Accuracy Test Loss Test Accuracy 1 2 0.1847 0.9204 0.2184 0.9079 0.2362 0.9015 2 6 0.1655 0.9273 0.2171 0.9074 0.2317 0.9035 3 6 0.1657 0.9273 0.2050 0.9123 0.2212 0.9058 4 11 0.1552 0.9315 0.2083 0.9134 0.2247 0.9043 5 6 0.1656 0.9273 0.2144 0.9090 0.2245 0.9043 As we can see in the previous table, in terms of losses (both in validation and in test) the best is the third training. The third training is also the best in terms of test accuracy, although the best in terms of validation accuracy is the fourth training. The best weights from each train were saved. Emphasize that only the results 3.4 Creating a neural network 35 of loss and accuracy in validation will be taken into account to select the best model, although we also present the test results for completeness. We can see below the learning curves for training and validation in trainings 3 and 4: Learning curve for the third training Learning curve for the fourth training Figure 3.27 Learning curves of the two best training sessions. It can be seen how in a few epochs the best result is achieved. However, the training is noisy, oscillating around the optimum. To have a more stable training a lower learning rate can be set, so the learning would be slow and sustained. See the results with a learning rate of 1e−4: Table 3.2 getUnet training results, LR 1e−4. getUnet LR 1e−4training results Train Number Epochs Train Loss Train Accuracy Valid Loss Valid Accuracy Test Loss Test Accuracy 1 18 0.1632 0.9282 0.2099 0.9091 0.2369 0.8974 2 7 0.2181 0.9128 0.2375 0.9017 0.2564 0.8917 3 18 0.1703 0.9255 0.2232 0.9049 0.2434 0.8950 4 9 0.1799 0.9229 0.2225 0.9081 0.2319 0.9041 5 16 0.1718 0.9250 0.2187 0.9080 0.2349 0.9017 In this case, the best training sessions are the first and fourth, both in terms of validation and test losses and in terms of accuracy. However, we see that the values achieved are worse than in the previous training sessions. We can also see the learning curves: Learning curve of first training Learning curve of fourth training Figure 3.28 Learning curves of the two best trainings with LR 1e−4. 36 Chapter 3. Design of a Deep Learning Model for Thread Counting As we can see, the training sessions are slower and more stable. We think that better results are achieved in the previous case because the training sessions are noisier. Precisely this oscillation around the optimum causes a lower minimum to be reached in validation than when performing a more stable training. It is also interesting to see the results obtained for some validation images because comparing the results for validation images can be useful to choose the best model. We present below the results for five randomly chosen validation images for the best training sessions. In the figures presented below we can see the following sequence by rows: original image, segmentations and label: Validation Images Segmentation obtained with Training 3, LR 1e−3 Segmentation obtained with Training 4, LR 1e−3 Segmentation obtained with Training 1, LR 1e−4 Segmentation obtained with Training 4, LR 1e−4 Labels Figure 3.29 Segmentation of validation images: . In the group of images we see that training 4 with LR 1e−3 seems to offer better results than training 3 with the same LR. This can be seen especially in the third of the images, where the segmentation obtained with the weights of training 4 completely detect the horizontal thread of the lower part. They both also present problems in the segmentation of the first image, where both have problems with the thick horizontal thread of the lower part. Still, the results are generally good. If we also include the trainings with LR 1e−4 in the comparison, we could conclude that training 4 is somewhat better, especially in the 3.4 Creating a neural network 37 first and third image. In the first, training 4 better detects the thick thread of the lower part, adding fewer spurious points. In the third, it detects the central vertical thread somewhat better. In addition, the difference in thickness of the segmented points between the two best trainings of the model with LR 1e−4 is also striking. At this point we must take into account the purpose we are pursuing, which is none other than measuring the distance between the crossing points. To measure the distance between the segmented points it would be necessary to binarize the segmented image to have a binary image and then measure the distance in some way. For our purpose, the best model will not be the one that performs the most identical segmentation to the label, but the one that best locates the crossing points. As we can see, the label points have a certain thickness and the loss function penalizes if the segmentation is not identical, that is, if segmented points with the same thickness are not obtained (since the error is measured for each pixel). At this point in the project, the validation loss and validation accuracy metrics will help us determine the best models, but later the one that offers the best thread count will be taken as the best model. In any case, the best models in terms of losses and accuracy seem to be the ones that would offer a better count, since they are the ones that segment the vast majority of points. In view of these results, we decided to try another modification to the model. The kernel size of the first and last layers changed from 3x3 to 7x7. This change was motivated by the difference in the size of the threads between the images in the dataset. We think that increasing the size of the kernel in the first layer allows the network to see a more global context in the first feature extraction. Once again, we perform another round of model training with the described hyperparameters and callbacks, keeping the learning rate at 1e−3 , which is the rate with which we have obtained the best results so far. We will call this new model getUnet6, and it has about two million two hundred thousand parameters. We present below the results obtained in the 5 training sessions: Table 3.3 getUnet6 training results. getUnet6 training results Train Number Epochs Train Loss Train Accuracy Valid Loss Valid Accuracy Test Loss Test Accuracy 1 8 0.1621 0.9286 0.2034 0.9140 0.2121 0.9099 2 7 0.1628 0.9284 0.2085 0.9095 0.2262 0.9024 3 9 0.1611 0.9291 0.2062 0.9119 0.2217 0.9035 4 12 0.1640 0.9279 0.2118 0.9120 0.2501 0.8990 5 21 0.1521 0.9328 0.2062 0.9100 0.2271 0.9003 As we can see in the table, the results are similar to those obtained with the original model. Although in general terms the results are somewhat better, the variation with respect to the previous case is not very notable. The first training improves previous results in both loss and accuracy (in validation and test respectively) compared to the fourth training of the getUnet model, which was the best up to now. Let’s look at the learning curve for the first training of the getUnet6 model: Curva aprendizaje primer entrenamiento Figure 3.30 Learning curve of the best model training getUnet6 with LR 1e−3. 38 Chapter 3. Design of a Deep Learning Model for Thread Counting And let’s also see the segmentation of the test images that we saw earlier with the first training of the getUnet6 model. The segmentation obtained with training 4 and LR 1e−3 of the getUnet model, which is the best so far, is also presented in order to be able to compare the segmentations visually in a comfortable way: Validation Images Segmentation obtained with the U-Net getUnet6 Best segmentation obtained with the U-Net getUnet Label Figure 3.31 Segmentation of validation images: Training 1 with LR 1e−3getUnet6. There are not several differences. In the first image, the getUnet6 model better resolves the thick horizontal thread area, besides to not adding spurious points in the central area of the image as the getUnet model does. However, in the third image, the getUnet6 model fails in the bottom horizontal thread, which is fully detected by the getUnet model. This area is really complex, since even appreciating the image it is not clear if there is a thread there or not. In the rest of the images we can see how the getUnet6 model tends to segment the crossing points with a more uniform shape and with a higher gray level, that is, with less uncertainty. In short, it does not seem clear that neither is expressly better than the other, although both offer quite acceptable results and good values for the validation loss and validation accuracy metrics. Tune some hyperparameters When designing an architecture from scratch, it is worth asking if the architecture designed will be optimal for the problem. There are a series of hyperparameters that shape the architecture and govern the training that can be modified: the number of layers, the number of convolutional filters in each layer, the learning rate or the dropout value are some examples. In order to optimize these values, there are tools such as the Keras Tuner package [71]. Its use is quite simple and intuitive: it only requires defining the search ranges of the hyperparameters to be tuned and selecting one of the search algorithms. To define the search ranges of the hyperparameters, you must define a Model-building function that inherits from the HyperModel class. Using the hp object we can define the ranges of the hyperparameters that we want to tune. In our case we have specified the following variations: 3.4 Creating a neural network 39 •Dropout between 0 and 0.5 • The number of filters on the model layers. The number of filters in the first layer is defined as an even integer between 6 and 24. Remember that, due to the way we have designed the architecture, this number doubles as it goes down the network. •Learning rate value of 1e−3or 1e−4. •Apply or not batch normalization in all layers. • The kernel size of the first and second layers (and their corresponding in the expansive path). Both kernels can take values of 3x3, 5x5 and 7x7 (each independently). • The number of layers in the model. Three possibilities are contemplated: do not add any, add one more layer or add two more layers with respect to the getUnet model. If a layer is added on the contracting path, its corresponding layer is added on the expansive path to maintain the U-shape of the architecture. Once the ranges of the hyperparameters have been defined, the search algorithm must be initialized. These algorithms provide an efficient way to test different combinations of all the hyperparameters that we want to tune. Logically, testing all possible combinations (grid search) would be so computationally expensive, and would also involve an enormous amount of time. Keras Tuner offers three search algorithms [72]: •Random search, which tests different combinations of hyperparameters in a totally random way. • BayesianOptimization, which starts by choosing the combinations randomly but, as the search completes, it chooses the next combination of hyperparameters according to the previous experience of combining them. • Hyperband, which performs some initial random training with few epochs to determine the best candidate values for the hyperparameters. This process is repeated iteratively to tune the search. We have used the BayesianOptimization algorithm with the following configuration: •The objective is to minimize losses in validation. • A maximum of 30 trials will be made. Each trial is a different combination of the hyperparameters. In addition, for each trial, 3 executions will be carried out. That is, the training will be repeated three times for each combination of hyperparameters. • An early stopping callback is established. So, in case of not improving the validation loss in three epochs, the execution is stopped and the next execution of that trial starts. Keras Tuner generates a .json file for each trial. In this file we can see the chosen combination for hyperparameters and the values of losses and accuracy that it reaches, both in training and in validation. The average time for each trial has been about 45 minutes. Keras Tuner provided the following results: •0.5 dropout rate. •18 filters on the first layer. This number is duplicated on each new layer. •1e−3Learning rate. • 5x5 convolutional filters in the first layer, intermediate value between the models getUnet and getUnet6. •3x3 convolutional filters in the second layer, as in getUnet and getUnet6 models. •Apply batch normalization on all layers. •Without adding any new layer. Regarding the getUnet6 model, we see that the number of filters increases from 16 to 18, the dropout increases from 0.15 to 0.5 and the size of the convolutional filters of the first layer changes from 7x7 to 5x5. Note the high dropout rate suggested by the tuner, since this dropout value implies randomly deactivating 50% of the neurons during training to avoid overfitting. With all these values we repeat a new session of 5 trainings of this new model that we will call getUnet KT: 40 Chapter 3. Design of a Deep Learning Model for Thread Counting Table 3.4 getUnetKT training results. getUnetKT training results Train Number Epochs Train Loss Train Accuracy Valid Loss Valid Accuracy Test Loss Test Accuracy 1 11 0.1982 0.9147 0.2167 0.9071 0.2449 0.8957 2 13 0.1651 0.9267 0.2076 0.9094 0.2285 0.9006 3 14 0.1597 0.9300 0.2063 0.9128 0.2120 0.9089 4 10 0.1614 0.9288 0.2071 0.9130 0.2133 0.9068 5 9 0.1983 0.9146 0.2304 0.9017 0.2445 0.8973 The disparity in results between some training sessions and others is striking. This is due to the high value of dropout. This dropout rate causes the 50% of neurons to be deactivated. Due to this, the resulting model sometimes works very well, but sometimes it works poorly. Specifically, we can see how trainings 1 and 5 are not good at all, while trainings 3 and 4 offer very good results. Even so, the results are similar to what we got with the getUnet6 model. As already discussed, it will be in the later chapter where we will choose the best model, since we are not so interested in having the best pixel-by-pixel segmentation as in having the best estimate of thread densities. Custom output layer We also tested one last interesting variant of the model. As we know, the last layer has a sigmoid activation function, which provides as output the gray level of each pixel, that is, a value between [0-1]. We can think that if said sigmoid function were narrowed, the number of pixels receiving an intermediate gray level would be reduced, forcing the output of the network to be more binary. However, if the sigmoid is too much narrowed it will become a step function that couldn’t be used for neural network training, since it has a discontinuity at 0 and is therefore not differentiable, so the backpropagation algorithm would not work. Sigmoid narrowing can be simulated if the parameters that govern the shape of the sigmoid were two more weights of the network, that is, two parameters that the network could adjust, as stated in [74]. These two parameters govern where the sigmoid is centered (defaults to 0) and how steep it is. By creating this new custom activation function (called Custom Sigmoid, CS) we could make the learning a bit easier. Now, the model has freedom to shape the last activation function in the way that gets the best results. We present below the results of adding this modification to the model getUnet KT. This model will be called getUnet KT CS: Table 3.5 getUnetKT CS training results. getUnetKT Custom Sigmoid training results Train Number Epochs Train Loss Train Accuracy Valid Loss Valid Accuracy Test Loss Test Accuracy 1 15 0.1601 0.9296 0.2009 0.9139 0.2199 0.9069 2 16 0.1591 0.9300 0.2071 0.9133 0.2223 0.9070 3 15 0.1596 0.9294 0.2086 0.9122 0.2272 0.9048 4 8 0.1717 0.9262 0.2059 0.9113 0.2346 0.8993 5 13 0.1627 0.9291 0.2059 0.9119 0.2169 0.9058 In general, the results are quite good, especially the first and last training. Lastly, let’s see the segmentation obtained with the best training of the getUnet KT and getUnet KT CS models: 3.4 Creating a neural network 41 Validation Images Segmentation obtained with the third training of getUnet KT Best segmentation obtained with the first training of getUnet KT CS Figure 3.32 Segmentation of validation images with getUnet KT and getUnet KT CS. No major differences are observed in either case with respect to the previous results for this little set of five validation images. Both seem to resolve the thick thread in the first image relatively well, as well as the horizontal thread in the lower part of the third image. In short, we have several variants with similar performances in the absence of being tested for their true purpose: the estimation of the thread densities per cm. All the weights of all the trainings have been saved to be able to be loaded later. So, will be able to evaluate the performance of each model when estimating the thread densities per cm. Other improvements tested Finally, briefly comment on a couple of techniques to improve model training, although without noteworthy novelties. We tested the learning rate schedule and weight decay techniques. The learning rate schedule [75] allows adjusting the learning rate in each epoch, making the value of the learning rate vary throughout the training as desired. Until now, a constant learning rate had been applied in all training epochs. It was tried to select a learning rate decay, that is, make the learning rate decrease as the training progresses. In this way, we favor a quick initial adjustment of the weights and, later, a more leisurely and stable learning. So, the parameters of the network gradually refine their value in a controlled manner. We tried to retrain the best model we had so far with the following learning rate decay: LRi+1=LRi 1 + decay rate ·epochi (3.6) Adecay rate of 1/5 was fixed, although we did not get better results than we had previously. On the other hand, it was tried to add the regularization method known as weight decay [76]. This method consists of adding the value of the weights squared (because the weights can be negative) to the loss function and multiplied by a constant called wd that weights how much the value of the weights affects the loss function. This constant also prevents the model from setting all weights to zero to make the loss function decrease. Loss =Binary crossentropy +wd∗(W2)(3.7) 48 Chapter 4. From segmentation to thread counting If we superimpose the located points to the original image, we can see what points we are getting: Figure 4.9 Location of points using local maxima. The result is quite satisfactory for this presented image and for others. From here on, we will consider both thresholding methods (scaled Otsu’s method and the local maxima filter-based method) to perform the following tests, so that we can see which one offers better results for the validation and test sets. 4.2.2 Spatial count First of all, we will implement a spatial count algorithm that allows us to measure the distance in pixels between a crossing point and the rest. Therefore, the first step is to locate the neighborhood crossing points for any crossing point in the image. To perform this task we resort to the KDTree decision tree algorithm [83]. This algorithm is a binary decision tree, so that the nodes represent hyperplanes that separate data points. They are widely used for searching for nearby neighbors. Once the tree is built taking as input the matrix containing the coordinates of the centroids, the query method is used to obtain the nearest neighbors. In our case we will search for the 9 closest neighbors to each centroid. This value is sufficient to always find a neighbor. We have programmed our spatial counting function so that if there are less than 9 centroids, an error exception is thrown. From the search for the 9 nearest neighbors we obtain two matrices: the first is the distance matrix, of size Mx9 (where M is the number of detected centroids), which contains in each row the distance from each centroid to each of the 9 nearest neighbours. The second matrix obtained is the matrix of indices, also of size Mx9, which contains in each row the indices of the nearest neighbors of the crossing point represented in that row. Since we already have the distance matrix, it only remains to see what distances must be taken into account and perform the calculations. We will therefore iterate through each row of these matrices to compute the distance from each crossing point to its neighbors. Since we have in the index matrix the coordinates of each neighbor, we start by defining a series of complex vectors that go from the crossing point that is being analyzed to each of the neighbors found. This vector allows us to obtain the angle between the crossing point in question and each of the neighbors. This angle between a crossing point and its neighbors will help to discriminate the following cases, taking into account the way in which the function scipy measurements center of mass returns the coordinates of the center of mass: if angle =90 then Neighbors are to the right of the crossing point end if if angle =−90 then Neighbors are to the left of the crossing point end if if angle =0then Neighbors are below the crossing point end if if angle =±180 then Neighbors are above the crossing point end if 4.2 Measure the distance between segmented crossing points 49 A tolerance of 25ºis left with respect to these angle values in case they are not completely exact or in case there is a certain inclination of the thread pattern. All that remains is to measure the distance to the nearest neighbor of all those found in each of the four cases described. Thus, we would measure the distance to the nearest neighbors to the left and right, that is, the distance between vertical threads, as well as the distance to the nearest neighbors above and below, that is, the distance between horizontal threads. All these distances are already calculated in the corresponding matrix, we just have to access the correct positions which are those distances that represent the distance of each crossing point to the selected neighbors. This is easy because we also have the index array. The four distances measured by each crossing point as well as the angles that each crossing point makes with each of the neighbors of interest are stored in four new arrays (two for distances and angles of vertical threads and the other two for the horizontal ones). If no neighbor is found in any of the four directions, the corresponding distance and angle will not be added. This can happen at edges, where there will always be at least one direction in which no neighbors exist. We can summarize all this graphically in the following figure: Figure 4.10 Graphical representation of the spatial count method. This is repeated for all rows of the index matrix, that is, for all localized crossing points. Once we have completed the process, we have a matrix with all the distances of interest and another with all the angles of interest. Finally, we also add the possibility of eliminating spurious distances that could be the product of some poorly segmented points. That is, if all the distances are in a certain range but there is some value much above or much below, it might be convenient to eliminate them. The option of passing an upper percentile value and a lower percentile value as a parameter to the function is present, so that the distances above or below the percentiles would be eliminated. Finally, the total distance is calculated as the average of all the distances that have remained in each of the two directions. Likewise, the average tilt of the image is calculated by means of the average of all the angles in each of the two directions. As for the distance, once the average distance in pixels has been obtained, the average distance between threads/cm is obtained as follows: Vertical count (threads cm ) = 200 (pix cm )/Mean vertical distance (pix)(4.1) Horizontal count (threads cm ) = 200 (pix cm )/Mean horizontal distance (pix)(4.2) Therefore, the described spatial count function receives as input an image (and percentile values if desired) and returns the density of vertical and horizontal threads per cm, as well as the tilt of the vertical and horizontal threads. 4.2.3 Frequency count Another option that we could consider is to carry out a frequency analysis of the segmentation. The idea is exactly the same as the one described in previous sections: calculate the 2D-DFT of the image and obtain the maximum of the spectrum on the vertical and horizontal axes. As we know, these maximums represent the value of the vertical and horizontal count. Although we have previously seen how FT applied directly to the cut patches of the plate has several limitations, we believe that its application to the segmented image may be interesting, since a priori they will be much less noisy images in which the repetition patterns will be more clear. Even so, FT may still not be able to deal well with cases where the distance between threads is 50 Chapter 4. From segmentation to thread counting irregular. The algorithm to obtain the thread density from the FT was provided to us. It could be described in the following steps: •The spectrum is calculated by the two-dimensional DFT of Nfft points using Python functions. • A priori values that will define the search area in the spectrum are set: a minimum count of 4 threads/cm, a maximum count of 25 threads/cm and a maximum angle of inclination of the image of 25º. All our images fit within these parameters. This step is necessary, as described in [16]. •The above values are converted to DFT points as follows: minFFT (DFT points) = 4(threads cm ) 200 pix cm ·N f ft (4.3) maxFFT (DFT points) = 25(threads cm ) 200 pix cm ·N f ft (4.4) • Rectangles that shapes the maximum search area are defined. The size of these rectangles will be given by the prior values set (once converted to DFT points). The minimum and maximum count describe how long the rectangle will be, while the maximum angle defines how wide it will be. Let us remember that the spectrum is conjugate symmetric due to image real values, so it is enough to define a rectangle between the origin and the positive X axis and another between the origin and the positive Y axis. • The maximum found in each rectangle is found and the Euclidean distance from the origin to the maximum of each axis is measured. The maximum on the X axis will give us the vertical distance, while the maximum on the Y axis will give us the horizontal distance. •The above DFT-points distances are converted to threads per cm following the equations above. We will test both the frequency counting method and the spatial counting method to see what results we get with each of them. 4.3 Measure the counting error with respect to the labels With all these ingredients we can already test how good our method is when it comes to estimating the thread densities per cm. To do this, we will compare the density of vertical and horizontal threads per cm that we obtain by applying all the steps described (segmentation using the U-Net + binarization + obtaining counts through spatial or frequency analysis) with the real density of vertical and horizontal threads per cm. The real distance between threads will be obtained by applying the spatial counting method directly to the labels and with percentiles (100-0), that is, without discarding any distance, since all the crossing points are correctly segmented in the labels and there are no spurious distances. We have created a function that measures, for each image, the difference between the real count and the one obtained with our method according to two different metrics: • The root mean square error (RMSE) of a set with n images will be measured first. This will give us an absolute error value, that is, in how many threads per cm we are making a mistake. RMSE =v u u t n ∑ i=1 (|valueapprox −valuereal)2 n(4.5) • The percentage error of a set with n images will also be measured, that is, what percentage of error we are making. To measure this error, the absolute value of the difference between the actual count and the estimated will be divided by the actual value. This metric is more interesting because it allows relativizing the error made based on the thread density of each image. It is not the same to count a thread more or less in an image with many threads or with few threads. %error = n ∑ i=1 |valueapprox −valuereal| |valuereal| n·100 (4.6) The function will generate a graph with the percentage error for each image in the set of images that we pass to it as input. Also, it will return the average RMSE and percentage error of all input images. We will 4.3 Measure the counting error with respect to the labels 51 use both the set and the validation set to measure the error and see how the whole process works. Although test results are included for completeness, the choice of the best model will be made based on the validation results. Comment at this point that the test set is quite hard. Of the 5 paintings present in the test set, we highlight the following characteristics: • Plate P001180 [39], by Velázquez, is one of those in which FT does not work, as we saw in Figure 3.2. It has a medium-high thread density and the threads in one direction are only noticeable as small nodules on the crossing points. It’s a low noise plate and there are some similar images in validation. • Plates P07905 [36] and P07906 [37], by Manuel de Pret, are quite complex plates. They have a remarkable noise level and problematic areas, threads that appear and disappear and very distinct and irregular threads. They have a medium density with some areas of slightly higher densities. • Plates P998 [84] and P999 [85], by Murillo, are also complicated and medium count plates. It has threads of different thicknesses, some of them being quite thin, and there are few images with these characteristics in training and validation. They also have a lot of noise, cracks and damage. There is no low densities plates in test set, which if we remember the distribution of the training data, were the predominant ones in training and validation sets. In addition, the qualities of the images are not the best, although not all of them are bad. In total, we have 104 test images, as discussed previously. To perform some checks with validation and test images, we must take into account the different models that we have to evaluate, as well as the different variants that we can add in the procedure. We can summarize it as follows: • We have good results with the getUnet model with LR 1e−3 , with the getUnet6 model, with the getUnet KT model, and with the getUnet CS model. We will test all of these, discarding the getUnet model with LR 1e−4for having poorer loss and accuracy results. • Different forms of thresholding will be tested: Otsu’s method scaling the threshold by 1.2 (as described in the corresponding section) or local maxima filtering method. • Spatial counting and frequency counting can be used. In principle, percentiles will not be applied to eliminate spurious distances in the spatial count, since we want to see how all the developed algorithms work together. If we eliminate these distances we would be masking the problem. Later, it will be possible to debate whether to incorporate these percentiles to filter the distances or not. Therefore, we have four cases: binarize with Otsu and use spatial counting, binarize with Otsu and use frequency counting, binarize with the local maxima filter and use spatial counting, or binarize with the local maxima filter and use frequency counting. In the following figures we can see the errors made in validation and test sets in the four cases mentioned. The weights obtained in the 5 training sessions of each model have been taken into account: 52 Chapter 4. From segmentation to thread counting Validation RMSE Test RSME Validation percentage error Test percentage error Figure 4.11 Box diagram of error when Otsu binarization and spatial count are used. In these images we see the representation of the RMSE and the percentage error in the form of a boxplot [86]. These graphs allow to study the data with ease. For each model, a box is presented, whose interpretation is as follows: • Each box represents the interquartile range (IQR), that is, in the area of the box are the data between the first quartile and the third quartile of the entire distribution. • The central black line represents the median or second quartile. If it is located in the center of the box, then the median matches the mean and mode of the distribution. •The upper and lower whiskers represent the limits of beginning the outliers and their value is limited to 150% of the IQR. •Diamonds indicate outliers, that is, distribution values that exceed the IQR by more than 150%. In this first case we study the data obtained when using Otsu’s method to binarize together with spatial counting to measure the distance between crossing points. It highlights the great variability that exists in the getUnet KT model. This was to be expected, since we commented that using such a high dropout rate there would be trainings where the error would be low and others where the error would be high. This behavior is perfectly described by the wide box, that is, the IQR is very wide because the data is very disparate. The good performance of the getUnet KT CS model also stands out, since it is the one with the lowest mean error of all and exhibits a behavior with few variations between one training and another. If we look at which one offers the least error, the getUnet KT CS model is the best in RMSE, both in test and validation. Also, getUnet KT CS is the best in percentage error in validation, while the best in percentage error in test is the getUnet KT model. As for the choice of the best model we look at the validation data, it is clear that getUnet KT CS is the one that offers the best results in terms of counting, although the getUnet KT model is close and does not should be discarded. The best of all training sessions is the third training of the model getUnet KT CS with an RMSE of 0.141 and a percentage error of 1.17% in validation, followed by the third training of the model getUnet KT CS with a RMSE of 0.168 and a percentage error of 1.37%. In test, the best performance 4.3 Measure the counting error with respect to the labels 53 is provided by the third training of the getUnet KT model, with an RMSE of 0.191 and a percentage error of 1.60%. The results using the scaled Otsu method and spatial counting are quite satisfactory, obtaining really low values. See now what happens when Local MaxFilter binarization and spatial count are used: Validation RMSE Test RSME Validation percentage error Test percentage error Figure 4.12 Box diagram of error when Local MaxFilter binarization and spatial count are used. Again, the variability of the getUnet KT model is noticeable, although in this case there seems to exhibit less variability. The results are similar: the best model according to validation results are getUnet KT CS and getUnet KT, with the getUnet and getUnet6 models being close. In terms of numerical results, no notable differences seems to appear. The lowest validation error is provided by the third training of the getUnet KT CS model, with an RMSE of 0.147 and a percentage error of 1.25%. On the other hand, the lowest test error is provided by the third training of the getUnet KT model, with an RMSE of 0.189 and a percentage error of 1.60%. These results are practically identical to those obtained in the previous case, however there is a noticeable time difference between using the Otsu’s method or filtering local maxima. Local maxima filter is much slower in execution, so in view of the results we will consider that the Otsu’s method is the one that provides the most balanced results between performance and computational cost. Let’s see now the results when Otsu binarization and frequencial count are used: 54 Chapter 4. From segmentation to thread counting Validation RMSE Test RSME Validation percentage error Test percentage error Figure 4.13 Box diagram of error when Otsu binarization and frequencial count are used. When we include the frequency count in the comparison, we see that the results are drastically worse. It has been found that this method does not work well at all when it is applied to binarized segmentation of some ot our test images. The reasons are similar to what we described when we mentioned the situations in which FT did not work well. Most validation and test images are complex and have some noise, which sometimes causes the segmentation to be somewhat irregular. In addition, many of the images have irregular thread patterns that provokes segmentations with irregular distances between threads, something that, as we say, makes FT not working well. For all this, we can conclude that the spatial counting method works better than the frequency counting for our validation and test image sets. However, we will use both methods to compare them over a wider range when we generate the density maps of full paintings. At this point we can conclude that the best combination for the validation set is to use the Otsu’s scaled method to binarize segmentation and spatial counting to measure the distance between crossing points. Besides, we have checked that the third training of both getUnet KT CS and getUnet KT (in that order) offers the best validation results. Otsu binarization plus spatial count and the two training session mencioned for both getUnet KT CS and getUnet KT are also the best for the test set, although validation results are the ones having into account to choose the best model. In the following figure you can see the representation of the error made in the vertical and horizontal count for some test and validation images using the weights of the third training of the model getUnet KT: 4.3 Measure the counting error with respect to the labels 55 Validation Error Test Error Figure 4.14 Error in each image using 3rd training of getUnet KT model. We highlight the few errors above 10% in the test images, with all errors being above 10% caused by the horizontal count. In validation, errors are also generally low. There are some cases where the validation error is above 10%. This may be due to the fact that in validation we have all the Data augmentation samples, so there are samples with extra rotation that usually present difficulties. We can represent the samples in which more error occurs. It is represented in the following figures, from left to right: 56 Chapter 4. From segmentation to thread counting •The input image to be segmented. • The label corresponding to that image. The vertical and horizontal count values obtained by applying the spatial count method to said label are indicated. •The segmented image that the network provides as output. •The binarized image using the scaling Otsu’s method that we have described above. The vertical and horizontal count values obtained by applying the spatial counting method to the binary image resulting are indicated. Let’s see the test samples in which there is more than 8% error to interpret what happens: Figure 4.15 Test images where error bigger than 8% occurs using getUnet KT model. We can see that the images with more than 8% error are, in general, quite complex images. All of them have a high level of noise, cracks and damage. In some, like the first or the last, it is difficult even with the naked eye to determine how many threads there are in each direction. In addition, we see how the thread patterns are irregular, being especially conflicting in the last image, where we see thick threads and very 4.3 Measure the counting error with respect to the labels 57 fine threads that even seem to get lost. The segmentation of the network in these conditions is not entirely satisfactory, although it does a good job considering the scenario in which we operate. You can see how noisy and very complex areas are included as crossing points with an intermediate gray level, while there are crossing points that have not been detected. Let’s also look at the validation samples where there is more than 12% error to see what happens: Figure 4.16 Validation images where error bigger than 12% occurs using getUnet KT model. In this case it is evident that the first image has a considerable rotation, something that can difficult the segmentation. In fact, we can see how numerous false points appear in the lower right corner, possibly due to both the rotation and the presence of cracks in the image. The other two images are again complex. They have a lot of noise and there are very thick threads that break the pattern. The noise causes some threads to seem as only one, as we see in the central area of the second image. This results in a segmentation where numerous points appear between the true crossing points that cause the count to go up, especially in the vertical direction. Now we represent the error made in the vertical and horizontal count for the test and validation images using the weights of the third training of the model getUnet KT CS: 64 Chapter 5. Generation of density maps Figure 5.1 Plate P0998 shown with ImageJ. In red, the region of interest. ImageJ makes it easy to target the pixel values that define the region of interest. Specifically, we need the row and column value of the pixel in the upper left corner as well as the row and column value of the pixel in the lower right corner. Once we have the path of the plate to read it and the information of the region of interest, we proceed to load the plate. To manage the reading and preprocessing of the plates, we were provided with a series of already programmed Python scripts. In them, object-oriented programming is used to manage the complete plate as an object, which has different attributes and methods that makes easy operating with it. When loading the plate with these libraries we can access to some attributes obtained from the metadata of the .ti f file of the plate, highlighting the following: •Number of rows in the image, that is, the height of the painting in pixels. •Number of columns in the image, that is, the width of the painting in pixels. •The resolution with which the plate was digitized, in pixels/cm. With these three values we can operate to obtain information about the physical dimensions of the painting. Specifically, we can obtain the height and width of the plate in cm, dividing the size in pixels by the resolution. Besides, the width and height of the region of interest in pixels are obtained. To obtain the height of the region of interest (number of rows), the difference between the value of the row of the final pixel (lower right corner) and the value of the row of the initial pixel (upper left corner) is obtained. Likewise, to obtain the width of the region of interest (number of columns), the difference between the value of the column of the final pixel and the value of the column of the initial pixel is obtained. 5.1.2 Getting values to move across the plate Once the image has been read and the region of interest has been defined, it is necessary to set and calculate another series of values that will help us to move across the entire painting and obtain the count. On the one hand, we have the values and attributes that will serve to preprocess the extracted patches and on the other the values and attributes that will serve to calculate the number of necessary windows in which we will subdivide the image. Preprocessing values The patches used to train the network that were provided to us had been preprocessed. That is, they applied preprocessing to improve the patches appearance and so facilitate labeling and learning of the network when these patches where cropped from the plates. For this reason, we were also provided with the libraries already programmed to carry out this preprocessing in the different windows in which we will subdivide the image. 5.1 Function to generate density maps 65 The preprocessing consists of a local mean filter, a local standard deviation filter and a tail clipping of the histogram of the image gray levels. The local mean filter consists of calculating the mean of the pixels in a window of a certain size, subtracting the mean value obtained from all the pixels in that window. The local standard deviation filter consists of calculating the value of the deviation of the pixels in a window of a certain size, dividing the value of all the pixels in the window by the calculated deviation. The tail clipping of the histogram is performed to discard its tails, so that if the image is subsequently scaled, it avoids very high or very low ranges to be assign to few values. The preprocessing windows were set to 27 pixels, as this is the value they previously used to extract the training patches. In addition, a margin is taken around the patch to be preprocessed to improve its result when performing the preprocessing. The margin must be at least half of the window. Taking into account that we perform a double filtering with 27 pixel windows and that for each filtering the margin must be at least half, the total margin for both filters must be at least twice the margin of each of the filters, this is, 27 pixels. We have set the value of the margin to 30 pixels. This value is sufficient and it is the same value that the ones used to obtain the training patches. In addition, a resampling is performed if the resolution of the image is not adequate. Let us remember that all the patches that were provided were preprocessed and their resolution adjusted to 200 pixels/cm, therefore all the windows in which we divide the image must also be preprocessed patches in the same way and with a resolution of 200 pixels/cm . For this, the image attribute targetResolution is accessed and set to 200 pixels/cm. On the other hand, it is necessary to define the division into windows of the plate to count each region separately. To obtain the density of threads per cm, we need each window in which we subdivide the image to have physical dimensions of 1cm x 1cm. Therefore, the size of the window in pixels will be given by the scanning resolution of the plate. In addition, we can add some overlap when scrolling the window to improve the final result and obtain smoother density maps. To control this overlap we use the shift fraction variable, which measures how much the window size advances in each iteration. The shift will be a fraction between 0 and 1, so a shift of 1 indicates that there is no overlap because the window advances its position by as many pixels as its size. This value must be chosen according to the dimensions of the plate studied: if the plate is small, a low shift can be used to obtain density maps with high resolution, while if the plate is very large, a high shift is more recommended to evit the process to be so lengthy. The shift fraction values used in each processed plate have been agreed with the tutor of the work. Once the shift fraction is fixed and the size of the window in pixels is obtained from the resolution, the number of pixels of overlap is obtained, that is, how many pixels the window advances in each iteration in rows or in columns: Shi ft (pixels) = Size window (pix)·Shi ft f raction With all this we can now subdivide the region of interest of the plate into windows. The fundamental idea is to obtain the number of rows and columns that the density maps will have, so that each pixel of the density maps physically represents a square window with a side of 1 cm. To obtain the number of rows, we divide the height of the region of interest (in pixels) by the number of pixels that the window advances in each iteration (shift in pixels). Likewise, to obtain the number of columns, we divide the width of the region of interest (in pixels) by the number of pixels that the window advances in each iteration (shift in pixels). The window advance value is the same for rows and columns. The windows, like the training patches, are squares with physical dimensions of 1 cm on each side. In both divisions the lower integer is taken in order to have an integer number of rows and columns. Taking the lower integer ensures that we do not go outside the region of interest. Once the number of rows and columns of the density maps is obtained, the corresponding empty matrices with those dimensions are constructed. The following image shows a graphical description of this division into windows. The Murillo P998 plate, which has a resolution of 200 pixels/cm, has been taken as a reference. Also, a shift of 1/2 has been set, that is, 50% of the window size. Therefore, each window with a side of 1cm has a size of 200x200 pixels. Each yellow or blue square represents a window. As we can see, they overlap 50% (100 pixels) due to the set shift. The difference in colors is only to allow the overlap of two consecutive windows to be differentiated in the example image. This overlap occurs both by rows and by columns, as exemplified in the image: 66 Chapter 5. Generation of density maps Figure 5.2 Division in windows of plate P998 by Murillo. Counting across the painting Now we divide the plate into 1cm side windows using a double f or loop that move acroos as many rows and columns the density maps have. Specifically, for a given row all columns will be iterated. For each row and column index of the density maps given by the double loop, the equivalent location on the plate is obtained. To do this we multiply the row and column index of the loop by the value of the shift. However, it should also be noted that the origin of the density map coordinates is located in the upper left corner of the region of interest. Therefore, to obtain the real index in row and column of the plate, it is necessary to add to the row and column of the initial pixel of the region of interest the row and column indicated by the double loop once transformed to plate pixel coordinates. This row and column value represents the top left corner of each window to process. The following operations are performed for each window: • A cut is made to obtain a patch of 1cm on each side, that is, with as many pixels as the resolution of the image originally. Those pixels will be counted down and to the right of the reference pixel (described in the previous paragraph) that represents the upper left corner. • If the resolution is not 200 pixels/cm, a resampling is performed to adjust the cropped patch to that resolution. Therefore, at this point we have a patch of 1cm side and 200x200 pixels. This matches the patches the U-Net was trained on. •The cropped patch is preprocessed as previously described. • The cropped patch is normalized to be in the range [0-1]. For this, we first subtract the minimum value of the image from all the pixels and, later, we divide all the pixels by the maximum value of the image. •We use some of the previously trained U-Net models to obtain the segmentation of the patch. • The segmentation is binarized using Otsu’s method and scaling the threshold by 1.2, as described in the previous chapter. • The spatial counting method without eliminate any distance is applied to obtain the number of vertical and horizontal threads in that patch, that is, the density of threads per cm in the vertical and horizontal direction. • Additionally, the thread density per cm will also be calculated using the frequency analysis method described. Specifically, it will be used to obtain the count from the segmentation and also to obtain the count by directly applying the FT to the original image once it has been cropped and preprocessed. In this way we can see the density maps obtained in the tree cases for comparison. • The result of the vertical and horizontal count obtained by each of the three ways (direct frequency analysis of the image, frequency analysis of the image segmentation or spatial count of the image segmentation) is stored in the input of the corresponding matrices to generate the density maps. In 5.2 Algorithm improvements 67 this way we will have the value of the number of threads per cm for each patch of the 1cm side square. In addition to the number of threads per cm, we also obtain the angle of deviation of the threads in that image from the spatial counting method. This other metric is also of interest as described at the beginning of this document. • In order to be able to follow the process and check the performance of all the algorithms involved, the obtained patch, the segmentation and the thresholding of the segmentation are sometimes represented. We also print the value of the count obtained with the three methods. Once the double loop has finished we can represent the density maps as an image, offering a simple interpretation of them through color maps. In addition, we can also generate some histograms from the density maps to check the distribution of the number of threads per cm in each of the directions. 5.2 Algorithm improvements In this section we will describe two variants of the previous algorithm that were developed to improve the computation time needed to obtain the density maps of a painting. Obtaining the density maps can take a long time, depending on the size of the painting and the capabilities of the machine where it is running on. At the end of this section we will compare the original algorithm and the two improvements proposed in this section in terms of execution time. 5.2.1 Build tensors to take advantage of GPU capabilities The first improvement consists of building tensors to be able to obtain the network prediction of several images simultaneously. Until now, each patch was passed to the network to obtain its corresponding segmentation individually. However, we could group several patches to obtain the segmentation of all of them simultaneously. If we grouped N patches, the network would receive as input a tensor of (N,200,200,1), that is, N images of 200x200 pixels in gray scale. To implement this we build row tensors, that is, for each row we create a tensor with all the extracted patches in all the columns of that row. Once the row tensor is built, it is passed to the network to obtain a tensor with the segmentation of all the patches of each row, saving a lot of time when operating with tensors instead of with individual images. Subsequently, each image of the row tensor and the segmentation tensor of that row is accessed to perform the count as we did until now. 5.2.2 Use multiprocessing libraries The second improvement, based on the previous one, consisted of using multiprocessing libraries to take full advantage of the capabilities of the machine where the algorithm is executed. Specifically, we use the multiprocessing [88] library and the h5py [89] library for file management to save memory. The idea of this version of the algorithm is to process several rows of the plate simultaneously, specifically as many rows as there are cores in the processor of the computer where it is executed. Now, patches of several rows are cropped and preprocessed at the same time, saving a considerable amount of time. To do this, groups of as many rows as the machine has processors are taken simultaneously. Each time a group of rows is preprocessed, the preprocessed row tensors are saved in an .hdf5 file. Each tensor has a particular identifier that allows the information to be retrieved later. This is done because each preprocessed row tensor can waste a lot of memory, so not every machine could have all the preprocessed tensors stored in memory. hdf5 files allow you to write and read from memory very quickly and efficiently, being perfect for this task. Once the entire plate has been preprocessed, we will read the tensors from the .hdf5 file and pass them through the network to obtain their segmentation. The simplest way would be to read the tensors one by one and obtain the segmentation of each one. However, we could group these tensors to build a mega-tensor. The size of the mega-tensor would be limited by the amount of memory that the GPU in question has, as we are looking to build the largest possible tensor that allows our machine’s ability to get the most out of the GPU. For this, the first thing is to obtain information about the free memory that the GPU has. Using the nvidia smi library we can get the amount of free memory. Likewise, we can consult in Python the size in memory that each 200x200 pixel image occupies. By dividing the free memory by the size of each image, we would obtain the maximum number of images that we could pass to the network at the same time. In practice, we 68 Chapter 5. Generation of density maps have multiplied this value by 80% to never cause an overflow in the GPU memory. What we do then is read several row tensors simultaneously, taking care that the memory waste of the samples of all the tensors does not exceed the threshold defined by the GPU memory. These tensors are joined in a mega-tensor that is passed to the network to obtain the segmentation of all those images. The segmentation tensors are also saved in another .hdf5 file to be able to read them later. This is repeated until there are no image tensors left to read. Between each group of tensors the GPU memory is cleaned to avoid overflow. We have verified that this is necessary because if not TensorFlow throws a memory exception, since if tensors are not removed from the memory of the GPU they remain stored and accumulated. In case the machine does not have a GPU, the network prediction will be much slower, since it must be done using the CPU. If a machine does not have GPU, each row tensor will be read separately and its prediction will be obtained. Finally, we also vary the way in which the segmented images are analyzed in order to parallelize the operations as much as possible. Now we will read several image tensors and segmentation tensors at the same time (as many as the machine has CPU cores). Through a loop, both tensors are recovered, thresholding the segmentation of each image and obtaining the count value through frequency and spatial analysis and performing the frequency analysis directly from the corresponding patch. The advantage is that you will be doing this on several pairs of images (original images and segmentation) from tensors at the same time. The thread count and angle estimation results are stored in matrices to be represented later. With this variant of the algorithm we get the most out of both the CPU and the GPU of each machine. A CPU-only variant of this algorithm has also been developed, that is, it uses the CPU for network prediction as well. This variant is intended to be used on machines that do not have a GPU, although the process slows down a lot. However, we think it is necessary to be tested, since we are interested in characterizing the analysis of plates in terms of computation time. To do this, we will compare how long it would take to generate the density maps for different paintings with different sizes and on different machines, some with CPU and others with CPU+GPU. That is why we will present two comparisons. The first will serve to measure the computation time for different plates using each of the three versions of the algorithm presented. The best version will be selected and another comparison of the execution time will be made for different paintings on different machines with CPU and GPU characteristics of all kinds. 5.2.3 Comparison of the three variants of the algorithm First of all, we execute the three variants of the algorithm on the machine that we have been using to train the network. Recall that this machine has an Intel Xeon E5-2630 v4 processor with 40 CPU Cores and two NVIDIA Tesla P100 16GB graphics cards as GPU. It will be tested on the following paintings: • P07905 by Manuel de Pret [36]. It has dimensions of 6095 rows by 7848 columns and a resolution of 200 pixels/cm. A shift fraction of 1/10 is applied. • P001180 by Velázquez [39]. It has dimensions of 44600 rows by 7100 columns and a resolution of 200 pixels/cm. 1/4 is applied. • P999 by Murillo [85]. It has dimensions of 5687 rows by 7290 columns and a resolution of 200 pixels/cm. A shift fraction of 1/10 is applied. • P007701 by Paret [90]. It has dimensions of 8060 rows by 6395 columns and a resolution of 200 pixels/cm. A shift fraction of 1/10 is applied. • P001171 by Velázquez [91]. It has dimensions of 34468 rows by 44507 columns and a resolution of 149 pixels/cm. A shift fraction of 1/2 is applied. • P1114 by Ribera [34]. It has dimensions of 17815 rows by 23622 columns and a resolution of 118 pixels/cm. A shift fraction of 1/2 is applied. The first three plates have samples in the test set. The last three are previously unreleased and were not part of any of the network’s training sets. As we can see, there are small paintings that have been assigned a low shift, like the first one, while others are huge and have been assigned a higher shift fraction, like the penultimate one. This is interesting to see the necessary processing times in each case. In the following table we collect the comparison of execution times of the complete program with each of the three variants of the algorithm presented. The first variant, where the row tensors were built, will be called the Tensor version. The second variant, in which multiprocessing libraries are used, will be called the Multiprocessing version: 5.2 Algorithm improvements 69 Table 5.1 Execution time of each version of the algorithm. Execution time of each version of the algorithm Plate Original algorithm Tensor version Multiprocessing version P07905 7h 6min 5h 52min 28min P001180 9h 48min 8h 1min 37min P999 7h 11min 5h 54min 29min P007701 11h 1min 9h 22min 46min P001171 1day 1h 20min 21h 2min 1h 49min P1114 14h 34min 12h 12min 58min As we can see, the execution time of the original algorithm is quite high, even on a powerful machine like the one used. Using the Tensor variant reduces processing time by 16% on average across all paintings. However, we see that by using the multiprocessing libraries, the computation time is reduced by more than 90%. On the tested machine, with 40 CPU cores, 40 rows of the plate are processed simultaneously. Also, by building the mega-tensors we are maximizing the use of GPU capabilities. Although it is true that resources are better used in terms of GPU, the main reason for the drastic decrease in execution times is in the total and parallel use of the CPU. Cropping and preprocessing tasks are quite expensive for CPU, so the time is notoriously reduced by parallelizing them. Next, we will try to characterize the computation time required by the paintings mentioned before when they are processed on machines with different capabilities. They will be tested on both powerful machines and laptops, with and without using GPUs. We try to give an idea of the time that may be required to process plates of different sizes (such as those previously selected) depending on the available machine. Note that we will always use the multiprocessing algorithm to improve times as much as possible. The machines used are the following: •PC Sedna: Intel Xeon E5-2630 v4 40 cores CPU + 16GB GPU Tesla P100. •PC Sedna without GPU - Intel Xeon E5-2630 v4 - 40 cores CPU (without using GPU). •PC Azken - INTEL i9 24 cores CPU + 48 GB GPU NVIDIA RTX A600. •PC Azken without GPU - INTEL i9 24 cores CPU (without using GPU). •MacBook Pro 2021 - 10 cores CPU + 16 cores GPU (16GB). •MacBook Pro Intel i5 2016 - 4 cores CPU (without GPU). We collect in the following table the execution times for each one of the paintings in the machines presented: Table 5.2 Execution times of the multiprocessing algorithm on different machines. Execution times of the multiprocessing algorithm on different machines Computer P07905 P001180 P999 P007701 P001171 P1114 Sedna GPU 28min 37min 29min 46min 1h 49min 58min Sedna not GPU 1h 30min 1h 57min 1h 34min 2h 20min 6h 21min 2h 42min Azken GPU 23min 33min 23min 39min 1h 30min 46min Azken not GPU 47min 1h 3min 48min 1h 12min 2h 41min 1h 24min MacBookPro 2021 M1Pro (GPU) 42min 55min 42min 1h 3min 2h 21min 1h 16min MacBookPro 2016 Intel (not GPU) 4h 56min 6h 35min 5h 7min 7h 34min 15h 24min 8h 25min We can see how the differences are notable depending on the machine on which it is executed. The conclusions that we can draw are the following: • At full capacity, the Azken machine is the one that offers the best results. Not surprisingly, it has the most powerful CPU and GPU in the entire comparison. • Although the difference between using or not using the GPU can involve up to three times the time, we see that the main difference is in the use of a powerful CPU. It is true that the task of obtaining the segmentation predictions is greatly accelerated when using a GPU. However, the hardest task in the entire pipeline is the preprocessing of the patches. Preprocessing consumes a lot of CPU resources, so 70 Chapter 5. Generation of density maps CPU capabilites are the bottleneck more than GPU. The Sedna and Azken machines have very powerful CPUs which mean that, even without using the GPUs to obtain the neural network predictions, the times are quite acceptable. • Counting of large plates can be got even on laptops, though they need to have good hardware and preferably a GPU that can run TensorFlow. This is the case of the MacBookPro with M1 Pro chip. With this chip it is possible to use the laptop’s GPU with Keras and TensorFlow libraries. If we add to this a fairly powerful CPU (for a laptop) to which we are taking full advantage, result are also good. • With MacBookPro with Intel i5 times increase exponentially. This laptop does not have a GPU and its quad-core processor is quite limited for tasks of this demand. It would not be advisable to work with machines of this type because they take almost five hours for small plates such as P07905 plate. For all the results presented, we conclude that the proposed method to obtain the density maps of a painting requires certain computational capabilities to be executed in logical times. The algorithm has been optimized as much as possible, implementing two improvements that have significantly reduced the time it takes to obtain the density maps by more than 90%. The main requirement is a powerful CPU, since the most demanding calculations require high computing capacities to be parallelized. In addition, it is highly recommended to have a GPU compatible with Keras and TensorFlow to reduce the prediction time of the neural network as much as possible. 5.3 Density maps obtained in some experiments In this last section we will see the results obtained when applying the proposed algorithm to complete plates. The density maps obtained for a series of paintings of the provided collection will be presented below. This section will also help us to compare the results that would be obtained in the following three scenarios: • When performing a frequential analysis of the plate. This was the most used method until now and it is the one that has been used in softwares such as Aracne. •When performing a frequential analysis of the segmentation obtained with the neural network. •When performing a spatial analysis of the segmentation obtained with the neural network. Besides, the density maps will be obtained using the getUnet KT model and the weights from its third training as well as using the getUnet KT CS model and the weights from its third training. We will use the following paintings to obtain their density maps: • Two bunches of grapes with a fly (reference P07905) by Manuel de Pret [36]. It has dimensions of 6095 rows by 7848 columns and a resolution of 200 pixels/cm. A shift fraction of 1/10 is applied. There are samples drawn from this paintings in the test set. X-ray plate Example patch Figure 5.3 P07905 plate and example patch. 5.3 Density maps obtained in some experiments 71 • Two bunches of grapes (reference P07906) by Manuel de Pret [37]. It has dimensions of 6192 rows by 7968 columns and a resolution of 200 pixels/cm. A shift fraction of 1/10 is applied. There are samples drawn from this painting in the test set. X-ray plate Example patch Figure 5.4 P07906 plate and example patch. • Prince Baltasar Carlos, on Horseback (reference P001180) by Velázquez [39]. It has dimensions of 44600 rows by 7100 columns and a resolution of 200 pixels/cm. A shift fraction of 1/4 is applied. No samples of this painting have been used so far. X-ray plate Example patch Figure 5.5 P001180 plate and example patch. 72 Chapter 5. Generation of density maps • The farewell of the prodigal son (reference P998) by Murillo [84]. It has dimensions of 5633 rows by 7086 columns and a resolution of 200 pixels/cm. A shift fraction of 1/10 is applied. There are samples drawn from this painting in the test set. X-ray plate Example patch Figure 5.6 P998 plate and example patch. • The dissipation of the prodigal son (reference P999) by Murillo [85]. It has dimensions of 5687 rows by 7290 columns and a resolution of 200 pixels/cm. A shift fraction of 1/10 is applied. There are samples drawn from this painting in the test set. X-ray plate Example patch Figure 5.7 P999 plate and example patch. 5.3 Density maps obtained in some experiments 73 • Adam and Eve (reference P001692) by Rubens [38]. It has dimensions of 34035 rows by 25874 columns and a resolution of 138 pixels/cm. A shift fraction of 1/2 is applied. There are samples drawn from this painting in the training set. X-ray plate Example patch Figure 5.8 P001692 plate and example patch. • The Abduction of Europa (reference P001693) by Rubens [92]. It has dimensions of 35552 rows by 38583 columns and a resolution of 188 pixels/cm. A shift fraction of 1/2 is applied. There are samples drawn from this painting in the training set. X-ray plate Example patch Figure 5.9 P001693 plate and example patch. 80 Chapter 5. Generation of density maps count and the spatial count offer good results. This first patch is a complex image, with noise, threads of different widths and not completely equidistant, something that can make it difficult for the FT to obtain the count directly from the cropped patch spectrum. In the second patch all three methods work fine. Regarding the models comparison, both models seems to do a good job, although getUnet KT provides a more compact and stable segmentation, since all points have the same shape. 5.3.2 P07906 Results Now we present the density maps obtained for the P07906 plate. We start with the density maps obtained by performing the frequency analysis of the plate: Vertical density map Horizontal density map Figure 5.22 P07906 density maps performing a frequential analysis of the plate. Next, let’s see those density maps when performing the frequential analysis of the segmentation: Vertical density map Horizontal density map Figure 5.23 P07906 density maps performing a frequential analysis of the segmentation with getUnet KT. 5.3 Density maps obtained in some experiments 81 Vertical density map Horizontal density map Figure 5.24 P07906 density maps performing a frequential analysis of the segmentation with getUnet KT CS. Third, we also present the density maps by performing a spatial analysis of the segmentation: Vertical density map Horizontal density map Figure 5.25 P07906 density maps performing a spatial analysis of the segmentation with getUnet KT. Vertical density map Horizontal density map Figure 5.26 P07906 density maps performing a spatial analysis of the segmentation with getUnet KT CS. 82 Chapter 5. Generation of density maps We can also represent the inclination of the threads in each area of the plate, since this value has been obtained when performing the spatial count: Vertical angle map Horizontal angle map Figure 5.27 P07906 Angle deviation in degrees with getUnet KT /getUnet KT CS. Lastly, let’s see the comparison of the vertical and horizontal count histograms with each of the three methods: FT Horizontal counting FT Segmentation Horizontal counting Spatial Horizontal counting FT Vertical counting FT Segmentation Vertical counting Spatial Vertical counting Figure 5.28 P07906 Counting histograms with getUnet KT. 5.3 Density maps obtained in some experiments 83 FT Horizontal counting FT Segmentation Horizontal counting Spatial Horizontal counting FT Vertical counting FT Segmentation Vertical counting Spatial Vertical counting Figure 5.29 P07906 Counting histograms with getUnet KT CS. For this plate we see that the previous conclusions are repeated. The frequency analysis of the plate generates numerous intense blue zones at the top of the vertical density map and on the sides of the horizontal density map (see Figure 5.22). This indicates that the frequency analysis is not working as expected. The other two methods improve this behavior, especially the spatial count which, again, is the one that offers the best results (see Figures 5.25 and 5.26). In this case, we can appreciate some more noticeable difference between using one model or another. The getUnet KT model performs better than the getUnet KT CS model in some areas in the upper left region of the vertical density map. The tilt maps indicate that in this case there is some more noticeable rotation, especially in the left part of the horizontal density map and in some bands of the vertical density map. Finally, in the histograms we find the same behavior of the previous case. Let’s also see some patch examples: Example with getUnet KT model Example with getUnet KT CS model Figure 5.30 P07906 examples. 84 Chapter 5. Generation of density maps We see that the images are very similar, since the canvas of the plates are similar. In the first patch the frequency analysis of the plate is quite wrong in the vertical count, while between the frequency and spatial count of the segmentation there are almost two threads of difference. In the segmentation there is a vertical thread that has not been completely segmented, although in view of the patch it is not clear that this thread exists because there are areas where it seems to exist and areas where it does not. In fact, if you try to count the number of vertical threads in the patch there appear to be 14 horizontal threads, an intermediate result between both results. In the second patch the spatial count gives a small error, since it counts part of a vertical thread less. The two frequency counts seems to work a little bit better. As we can see, the pattern of the threads is quite homogeneous, a situation where the FT works correctly as we expect, Also, the FT is able to infer the pattern from such noisy segmentation. 5.3.3 P001180 Results Next we present the density maps obtained for the P001180 plate. Let us remember that this plate is one of those where the application of FT to the images directly extracted from the plate does not work well at all. As for P07905 and P07906 plates, some samples from this plate are present in the test set. Let’s see the density maps obtained by performing the frequency analysis of the plate: Vertical density map Horizontal density map Figure 5.31 P001180 density maps performing a frequential analysis of the plate. 5.3 Density maps obtained in some experiments 85 Next, let’s see those density maps when performing the frequential analysis of the segmentation. The first two maps (from left to right) belong to getUnet KT model, while the two last maps belong to getUnetKT CS model: Vertical density map Horizontal density map Vertical density map Horizontal density map Figure 5.32 P001180 density maps performing a frequential analysis of the segmentation. 86 Chapter 5. Generation of density maps Third, we also present the density maps by performing a spatial analysis of the segmentation. The first two maps (from left to right) belong to getUnet KT model, while the two last maps belong to getUnetKT CS model: Vertical density map Horizontal density map Vertical density map Horizontal density map Figure 5.33 P001180 density maps performing a spatial analysis of the segmentation. 5.3 Density maps obtained in some experiments 87 See now the inclination of the threads in each area of the plate, since this value has been obtained when performing the spatial count The deviation in degrees map obtained is practically the same for both model, so we only attach a figure that would be valid for both models: Vertical angle map Horizontal angle map Figure 5.34 P001180 Angle deviation in degrees with getUnet KT /getUnet KT CS. 88 Chapter 5. Generation of density maps Lastly, let’s see the comparison of the vertical and horizontal count histograms with each of the three methods: FT Horizontal counting FT Segmentation Horizontal counting Spatial Horizontal counting FT Vertical counting FT Segmentation Vertical counting Spatial Vertical counting Figure 5.35 P001180 Counting histograms with getUnet KT. FT Horizontal counting FT Segmentation Horizontal counting Spatial Horizontal counting FT Vertical counting FT Segmentation Vertical counting Spatial Vertical counting Figure 5.36 P001180 Counting histograms with getUnet KT CS. The fabric’s frequency analysis works quite poorly, as expected for the type of canvas it is (see Figure 5.31). The result improves a lot when applying the frequency count to the segmentation, as can be seen in Figure 5.32. FT works badly because the threads in one of the directions are barely visible, only perceived 5.3 Density maps obtained in some experiments 89 as small nodules at the crossing points. This is logically solved by applying FT to the segmentation. Finally, spatial counting (see Figure 5.33) further improves the result, offering very smooth density maps, in which the painting’s fingerprint can be perfectly perceived. The improvement obtained in the spatial count may be due to the fact that the segmentation sometimes presents irregular distance between threads, affecting the performance of the FT. From the horizontal map of the spatial count we can deduce that in the upper area there is a kind of seam across the width, which has a somewhat lower count value (it can be seen in yellow and green in the figure 5.33). This is undoubtedly related to the results of the tilt maps, in which we see that there seems to be some alteration in that area, since the rest of the tilt maps are practically neutral. For their part, the histograms show the poor result obtained with the counting based on the FT of the plate, as well as the good results with the other two methods. See now some examples for this plate: Example with getUnet KT model Example with getUnet KT CS model Figure 5.37 P001180 examples. In the first patch we see that the method based on the FT of the patch fails in the vertical count, that is precisely where the pattern is not appreciated. It’s perceptible that the vertical threads are only seen as a widening at the crossing points. The methods based on the segmentation do work really well. In the second patch we see that there is a series of vertical threads very close together in the center of the image, while there are thicker ones in other areas of the image. This means that the distance between threads is not very homogeneous, so neither FT of the patch and the frequency analysis of the segmentation work well. Spatial counting does provide satisfactory results in both cases. In the density maps obtained for this plate there are no great differences between those obtained with the getUnet KT model and those obtained with the getUnet KT CS model, although if we look at the spatial count density maps we would say that the result obtained using the getUnet KT CS model looks somewhat better, especially in the horizontal density map, since there seems to be more continuity in the regions with more intense red color. 96 Chapter 5. Generation of density maps Lastly, let’s see the comparison of the vertical and horizontal count histograms with each of the three methods: FT Horizontal counting FT Segmentation Horizontal counting Spatial Horizontal counting FT Vertical counting FT Segmentation Vertical counting Spatial Vertical counting Figure 5.53 P999 Counting histograms with getUnet KT. FT Horizontal counting FT Segmentation Horizontal counting Spatial Horizontal counting FT Vertical counting FT Segmentation Vertical counting Spatial Vertical counting Figure 5.54 P999 Counting histograms with getUnet KT CS. The canvas of this plate is very similar to that of P998, and this can be seen reliably represented in the density maps. The results on the maps are practically identical, with the same problems and advantages as before. The first method still provides blue areas where it is not able to get the count well, while the second 5.3 Density maps obtained in some experiments 97 method combined with the getUnet KT CS model seems to get the best results by bringing out the most intense red in some areas, although spatial counting also gives very good results. The horizontal inclination maps does exhibits some differences, since before there was a notable inclination of the horizontal threads in some areas and now there is not any inclination of the threads in the entire plate. In the histograms, the bimodal character seen previously appears again. Finally, let’s see some examples of specific patches: Example with getUnet KT model Example with getUnet KT CS model Figure 5.55 P999 examples. In the first patch, all three methods work quite well, which is expected since the pattern is generally well defined. In the second patch, notable differences appear in the horizontal count. It is due to the presence of a thick horizontal thread that breaks the homogeneity of the pattern, which means that only the spatial count works well. Emphasize that in some patches from both Murillo plates which has a lot of noise, the frequency count of the patch can be useful over the spatial count, since the FT is able to see the underlying structures of the image and obtain the count (something that the spatial count cannot do if the segmentation fails due to noise). 98 Chapter 5. Generation of density maps 5.3.6 P001692 Results We present the density maps obtained for the P001692 plate. Some patches from this plate were used in the training set. Let’s see the density maps obtained by performing the frequency analysis of the plate: Vertical density map Horizontal density map Figure 5.56 P001692 density maps performing a frequential analysis of the plate. Next, let’s see those density maps when performing the frequential analysis of the segmentation: Vertical density map Horizontal density map Figure 5.57 P001692 density maps performing a frequential analysis of the segmentation with getUnet KT. 5.3 Density maps obtained in some experiments 99 Vertical density map Horizontal density map Figure 5.58 P001692 density maps performing a frequential analysis of the segmentation with getUnet KT CS. Third, we also present the density maps by performing a spatial analysis of the segmentation: Vertical density map Horizontal density map Figure 5.59 P001692 density maps performing a spatial analysis of the segmentation with getUnet KT. 100 Chapter 5. Generation of density maps Vertical density map Horizontal density map Figure 5.60 P001692 density maps performing a spatial analysis of the segmentation with getUnet KT CS. We can also represent the inclination of the threads in each area of the plate, since this value has been obtained when performing the spatial count: Vertical angle map Horizontal angle map Figure 5.61 P001692 Angle deviation in degrees with getUnet KT /getUnet KT CS. 5.3 Density maps obtained in some experiments 101 Lastly, let’s see the comparison of the vertical and horizontal count histograms with each of the three methods: FT Horizontal counting FT Segmentation Horizontal counting Spatial Horizontal counting FT Vertical counting FT Segmentation Vertical counting Spatial Vertical counting Figure 5.62 P001692 Counting histograms with getUnet KT. FT Horizontal counting FT Segmentation Horizontal counting Spatial Horizontal counting FT Vertical counting FT Segmentation Vertical counting Spatial Vertical counting Figure 5.63 P001692 Counting histograms with getUnet KT CS. This plate is particular, since it presents in some areas the highest counts seen so far. In the density maps, it stands out that the application of the FT directly on the patches again produces intense blue areas, that is, areas where it does not count well (see Figure 5.56). This is mitigated in part by applying the FT on the 102 Chapter 5. Generation of density maps segmentation, since we can see how the density maps are much cleaner compared to the previous ones (see Figure 5.57). It’s true that there are blue areas, especially on the horizontal map, that do not correspond to errors but to counts of that density. The blue areas that we can attribute as errors are those that are randomly distributed throughout the plate, not those that form well-defined regions. The aforementioned improvement is even more noticeable when viewing the density maps obtained with the spatial count, where the different areas with different counts can be clearly appreciated (see Figures 5.59). The central area of the density maps is striking, where there seems to be some kind of vertical seam that unites two different fabrics. This partly explains the two orange and red areas that we see in the central part of the vertical density map. When the seam is there, the surrounding threads tighten and come together more, which causes that in each square centimeter analyzed the count increases. On the other hand, the difference in the density maps obtained with each of the U-Net models is notable on this plate. The getUnet KT CS model offers a clearer result, offering more uniform counts in practically the entire plate. In the tilt maps we see how there is a certain angle in the edge areas, possibly due to the nails. Finally, the histograms again show that a to the left when using the first method, since those unwanted blue regions that appear when using the FT directly on the patches causes those very low counts to be reflected in the histogram. The rest of the results and the histograms are the usual ones, with few differences between the second and third methods (beyond the resolution of the histograms, as usual). Remember that there were some patches extracted from this plate in the training set, so it was expected that the density maps would be good how it is the case. Let’s see some examples of this painting: Example with getUnet KT model Example with getUnet KT CS model Figure 5.64 P001692 examples. In the first image we see how all the methods offer good results. It could be expected that the FT of the patch would not offer a good count due to the presence of threads of different thicknesses, although surprisingly this has not been the case. In the second image we see that the spatial count is somewhat below the real result, which is obtained by the other two ways. In this case the thread pattern is very well defined and it was expected that the frequency methods would work well. 5.3 Density maps obtained in some experiments 103 5.3.7 P001693 Results Next we present the density maps obtained for the P01693 plate, which has a similar canvas to the P001692 and also has some samples in the training set. Let’s see the density maps obtained by performing the frequency analysis of the plate: Vertical density map Horizontal density map Figure 5.65 P001693 density maps performing a frequential analysis of the plate. Next, let’s see those density maps when performing the frequential analysis of the segmentation: Vertical density map Horizontal density map Figure 5.66 P001693 density maps performing a frequential analysis of the segmentation with getUnet KT. 104 Chapter 5. Generation of density maps Vertical density map Horizontal density map Figure 5.67 P001693 density maps performing a frequential analysis of the segmentation with getUnet KT CS. Third, we also present the density maps by performing a spatial analysis of the segmentation: Vertical density map Horizontal density map Figure 5.68 P001693 density maps performing a spatial analysis of the segmentation with getUnet KT. 5.3 Density maps obtained in some experiments 105 Vertical density map Horizontal density map Figure 5.69 P001693 density maps performing a spatial analysis of the segmentation with getUnet KT CS. We can also represent the inclination of the threads in each area of the plate, since this value has been obtained when performing the spatial count: Vertical angle map Horizontal angle map Figure 5.70 P001693 Angle deviation map with getUnet KT /getUnet CS. 112 Chapter 5. Generation of density maps 5.3.9 P49435-1 Results Next we present the density maps obtained for the P49435-1 plate. Some patches from this plate are in the training set. We start with the density maps obtained by performing the frequency analysis of the plate: Vertical density map Horizontal density map Figure 5.83 P49435-1 density maps performing a frequential analysis of the plate. Next, let’s see those density maps when performing the frequential analysis of the segmentation: Vertical density map Horizontal density map Figure 5.84 P49435-1 density maps performing a frequential analysis of the segmentation with getUnet KT. 5.3 Density maps obtained in some experiments 113 Vertical density map Horizontal density map Figure 5.85 P49435-1 density maps performing a frequential analysis of the segmentation with getUnet KT CS. Third, we also present the density maps by performing a spatial analysis of the segmentation: Vertical density map Horizontal density map Figure 5.86 P49435-1 density maps performing a spatial analysis of the segmentation with getUnet KT. 114 Chapter 5. Generation of density maps Vertical density map Horizontal density map Figure 5.87 P49435-1 density maps performing a spatial analysis of the segmentation with getUnet KT CS. We can also represent the inclination of the threads in each area of the plate, since this value has been obtained when performing the spatial count: Vertical angle map Horizontal angle map Figure 5.88 P49435-1 Angle estimation with getUnet KT /getUnet KT CS. 5.3 Density maps obtained in some experiments 115 Lastly, let’s see the comparison of the vertical and horizontal count histograms with each of the three methods: FT Horizontal counting FT Segmentation Horizontal counting Spatial Horizontal counting FT Vertical counting FT Segmentation Vertical counting Spatial Vertical counting Figure 5.89 P49435-1 Counting histograms with getUnet KT. FT Horizontal counting FT Segmentation Horizontal counting Spatial Horizontal counting FT Vertical counting FT Segmentation Vertical counting Spatial Vertical counting Figure 5.90 P49435-1 Counting histograms with getUnet KT CS. In the results of this plate we can see an habitual pattern. With the first method numerous spurious zones appear due to the problems that the FT sometimes has in order to obtain a valid count (see Figure 5.83). This behavior is reduced in part by applying frequency analysis to the segmentation (see Figures 5.84 and 5.85), and further reduced by applying spatial counting to the segmentation (see Figures 5.86 and 5.87). On the other hand, although it is true that the spatial count offers the cleanest results, we can see in the maps 116 Chapter 5. Generation of density maps generated with the frequency count of the segmentation how it seems that this method better counts the areas with high density in the horizontal density map (with the getUnet KT CS model is used). From the histograms we can again highlight the tail that appears in the histogram of the first method, product of the blue spurious areas mentioned. Let’s see some examples for this painting: Example with getUnet KT model Example with getUnet KT CS model Figure 5.91 P49435-1 examples with getUnet KT. In the first patch we can see how the FT of the patch offers an erroneous count in the horizontal direction, since it counts an extra thread, although in general the results with the three methods are similar and acceptable. Recall that there are samples of this plate in the training set. In the second patch we see how the spatial count is the one that offers the correct result in the horizontal direction. In this case, the thread pattern is not completely homogeneous, which can lead to errors in the frequency methods. 5.3 Density maps obtained in some experiments 117 5.3.10 P49435-2 Results Next we present the density maps obtained for the P49435-2 plate, which is so similar respect to the last analyzed plate and there also are some patches from this plate in the training set. We start with the density maps obtained by performing the frequency analysis of the plate: Vertical density map Horizontal density map Figure 5.92 P49435-2 density maps performing a frequential analysis of the plate. Next, let’s see those density maps when performing the frequential analysis of the segmentation: Vertical density map Horizontal density map Figure 5.93 P49435-2 density maps performing a frequential analysis of the segmentation with getUnet KT. 118 Chapter 5. Generation of density maps Vertical density map Horizontal density map Figure 5.94 P49435-2 density maps performing a frequential analysis of the segmentation with getUnet KT CS. Third, we also present the density maps by performing a spatial analysis of the segmentation: Vertical density map Horizontal density map Figure 5.95 P49435-2 density maps performing a spatial analysis of the segmentation with getUnet KT. 5.3 Density maps obtained in some experiments 119 Vertical density map Horizontal density map Figure 5.96 P49435-2 density maps performing a spatial analysis of the segmentation with getUnet KT CS. We can also represent the inclination of the threads in each area of the plate, since this value has been obtained when performing the spatial count: Vertical angle map Horizontal angle map Figure 5.97 P49435-2 Angle deviation in degrees with getUnet KT /getUnet KT CS. 120 Chapter 5. Generation of density maps Lastly, let’s see the comparison of the vertical and horizontal count histograms with each of the three methods: FT Horizontal counting FT Segmentation Horizontal counting Spatial Horizontal counting FT Vertical counting FT Segmentation Vertical counting Spatial Vertical counting Figure 5.98 P49435-2 Counting histograms with getUnet KT. FT Horizontal counting FT Segmentation Horizontal counting Spatial Horizontal counting FT Vertical counting FT Segmentation Vertical counting Spatial Vertical counting Figure 5.99 P49435-2 Counting histograms with getUnet KT CS. The results obtained for this plate are quite similar to those obtained for the previous plate, since both are very similar. Everything said then is verified here as well, although perhaps the presence of lower counts in all the density maps and histograms is striking. Let’s see some examples for this painting: 5.4 Match between density maps 121 Example with getUnet KT model Example with getUnet KT CS model Figure 5.100 P49435-2 examples with getUnet KT. In the first patch we see how the three methods work well despite the high level of image noise. In the second patch it is apreciable that frequential methods do not work well, since the gaps between the threads are variable in size and this makes spatial counting improving the other two methods. 5.4 Match between density maps To finish this chapter, let’s see if there is a match between the density maps of the similar paintings analyzed. We will use the density maps obtained by spatial counting of the segmentation provided by the getUnet KT model. It must be taken into account that the coincidences between plates can only occur in the warp, which is the one that presents a deterministic behavior. If we remember what was mentioned in the introduction to this document, the warp threads follow a pattern determined by how they were placed on the loom and how they were tightened. Therefore, all pieces of fabric cut from the same roll will have matching warp threads. However, this does not happen in the weft, since in that case the pattern is totally random (it can usually be easily modeled as a Gaussian). Knowing these details allows us to deduce that there can only be a match in the density maps in one direction, precisely the one that corresponds to the warp. Although the warp is usually associated with the vertical threads, it is possible that when fixing the fabric to the wood stretcher it has been rotated, so in that case the match would occur between the horizontal threads. 5.4.1 P07906 and P07906 We have not found any combination of the density maps of these two tables that provides the match sought. 128 Chapter 7. References [15] Johnson, D. H., Sun, L., Johnson, C. R. & Hendriks, E. (2010). “Matching canvas weave patterns from processing X-ray images of master paintings”, IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP). Dalas ,TX, USA, 19 March 2010. IEEE, pp. 958,961 [16] Johnson, D. H., Johnson Jr, C. R. & Erdmann, R. G. (2013). “Weave analysis of paintings on canvas from radiographs”. Signal Processing, 93, pp. 527–540 [17] Simois, F. J. & Murillo-Fuentes, J. J. (2018). “On the power spectral density applied to the analysis of old canvases”. Signal Processing, 143, pp. 253 – 268 [18] Universidad de Sevilla (2019). Desarrollan en la ETSI una herramienta para el estudio de lienzos en el Prado. Avaliable in: https://www.us.es/actualidad-de-la-us/desarrollan-en-la-etsi-una-herramienta-para-el- estudio-de-lienzos-en-el-prado [Consulted 08-19-2022]. [19] Maaten, L. & Erdmann, R. G. (2015). “Automatic thread-level canvas analysis: A machine-learning approach to analyzing the canvas of paintings”. IEEE Signal Process. Mag., 32(4), pp. 38-45 [20] Sizyakin, R., Cornelis, B., Meeus, L., Dubois, H., Martens, M., Voronin, V. & Pizurica, A. (2020). “Crack detection in paintings using convolutional neural networks”. IEEE Access, 8, pp. 74535 – 74552 [21] Roberto, J., Ortego, D., Davis & B. (2020). “Toward the automatic retrieval and annotation of outsider art images: A preliminary statement”, AI4HI. Marseille, 11-16 May 2020. ERLA, pp. 16-22 [22] Wikiart (2022). Arte del día. Avaliable in: https://www.wikiart.org/es [Consulted 08-17-2022]. [23] Nemade, R. Nitsure, A. , Hirve, P. & Mane, S.B. (2017). “Detection of forgery in art paintings using machine learning”. International Journal of Innovative Research in Science, Engineering and Technology, 6(5), pp. 8681-8692 [24] Pu, W., Sober, B., Daly, N., Higgitt, C., Daubechies, I. & Rodrigues, M., (2020). “A connected autoencoders based approach for image separation with side information: With applications to art investigation”, IEEE Int. Conf. on Acoustics, Speech and Signal Process. (ICASSP). 4-8 May 2020, Barcelona. IEEE, pp. 2213–2217 [25] Zou, Z., Zhao, P. & Zhao, X.(2021). “Virtual restoration of the colored paintings on weathered beams in the forbidden city using multiple deep learning algorithms”. Advanced Engineering Informatics, 50, pp. 1-15 [26] Bonafonte, A. (2008). Transformada de Fourier. Avaliable in: https://upcommons.upc.edu/bitstream/handle/ 2117/191346/tema2.transffourierv29may2009-2742.pdf [Consulted 08-22-2022]. [27] Tyagi, M. (2021). HOG (Histogram of Oriented Gradients): An overview. Available in: https://towardsdatascience.com /hog-histogram-of-oriented-gradients-67ecd887675f [Consulted 08-22-2022]. [28] MullOverThing (2021). What is an isotropic gaussian?. Avaliable in: https://mulloverthing.com/whatis-an-isotropic-gaussian/ [Consulted 08-22-2022]. [29] Kaggle (2021). ResNet-18. Avaliable in: https://www.kaggle.com/datasets/pytorch/resnet18 [Consulted 08-24-2022]. [30] Education Ecosystem (2018). Understanding K-means Clustering in ML. Avaliable in: https:// towardsdatascience.com/ understanding-k-means-clustering-in-machine-learning-6a6e67336aa1 [Consulted 08-24-2022]. [31] Na8 (2018). Clasificar con K-NN in Python. Avaliable in: https://www.aprendemachinelearning.com/ clasificar-con-k-nearest-neighbor-ejemplo-en-python/ [Consulted 08-24-2022]. [32] Papers With Code (2017). VGG. Avaliable in: https://paperswithcode.com/method/vgg [Consulted 08-24-2022]. 129 [33] Papers with code (2021). Exponential Linear Unit. Avaliable in: https://paperswithcode.com/method/elu [Consulted 08-24-2022]. [34] Ribera, J. (1632). Ixion, Museo Nacional del Prado (P1114). [35] Ribera, J. (1632). Ticio, Museo Nacional del Prado (P1113). [36] Pret, M. (1630-1644). Two bunches of grapes with a fly, Museo Nacional del Prado (P07905). [37] Pret, M. (1630-1644). Two bunches of grapes, Museo Nacional del Prado (P07906). [38] Rubens, P.P. (1630-1644). Adam and Eve, Museo Nacional del Prado (P001692). [39] Velázquez, D. (1634-1635). Prince Baltasar Carlos on Horseback, Museo Nacional del Prado (P001180). [40] Geron, A. (2019). Hands-on machine learning with Scikit-Learn, Keras and TensorFlow: concepts, tools, and techniques to build intelligent systems. O’Reilly. (2nd ed.). [41] Cruces, S. (2022). “Resampling techniques for model selection and predictive evaluation”. Universidad de Sevilla. [42] Samarth, A. (2021). How to split data into three sets (train, validation, and test) And why?. Avaliable in: https://towardsdatascience.com/how-to-split-data-into-three-sets-train-validation-and-test-and-why- e50d22d3e54c [Consulted 08-24-2022]. [43] DotCSV. (2020). Aprende qué son las redes neuronales. [YouTube]. 12 November 2020. Avaliable in: https://www.youtube.com/watch?v=MRIv2IwFTPglist=PL-Ogd76BhmcB9OjPucsnc2-piEE96jJDQ [Consulted 08-26-2022]. [44] Aravind, P. (2020). CNN vs. RNN vs. ANN – Analyzing 3 Types of Neural Networks in Deep Learning. Avaliable in: https://www.analyticsvidhya.com/blog/2020/02/cnn-vs-rnn-vs-mlp-analyzing-3-types-of-neural- networks-in-deep-learning/ [Consulted 08-26-2022]. [45] Doshi, S. (2019). Various Optimization Algorithms For Training Neural Network. Avaliable in: https://towardsdatascience.com/optimizers-for-training-neural-network-59450d71caf6 [Consulted 08-26- 2022]. [46] Shah, T. (2017). About Train, Validation and Test Sets in ML. Avaliable in: https://towardsdatascience.com /train-validation-and-test-sets-72cb40cba9e7 [Consulted 08-26-2022]. [47] Grupo Atico34. (2022). Overfitting. Qué es, causas, consecuencias y cómo solucionarlo. Avaliable in: https://protecciondatos-lopd.com/empresas/overfitting/ [Consulted 08-26-2022]. [48] Minaee, S., Boykov, Y., Porikli, F., Plaza, A., Kehtarnavaz, N. Demetri Terzopoulos. (2022). “Image Segmentation Using Deep Learning: A survey”. IEEE TRANSACTIONS ON PATTERN ANALYSIS AND MACHINE INTELLIGENCE, 44(7), pp. 3523-3542 [49] Lamba, H. (2019). Understanding semantic segmentation with UNET. Avaliable in: https://towardsdatascience.com /understanding-semantic-segmentation-with-unet-6be4f42d4b47 [Consulted 08-26-2022] [50] Sherstinsky, A. (2020). “Fundamentals of Recurrent Neural Network (RNN) and Long Short-Term Memory (LSTM) Network”. Physica D: Nonlinear Phenomena, 404, pp. 1-43 [51] Vaswani, A., Shazeer, N., Parmar, N., Uszkoreit, J., Jones, L., Gomez, N., Kaiser, L. Polosukhin, I. (2017). “Attention is all you need”, NIPS’17: Proceedings of the 31st International Conference on Neural Information Processing Systems, pp. 6000-6010 [52] Simonyan, K. Zisserman, A. (2014). “Very deep convolutional networks for large-scale image recognition”. ICLR 2015, pp. 1409- 1423 130 Chapter 7. References [53] Szegedy, C., Liu, W., Jia, Y., Servant, P., Reed, S., Anguelov, D., Ethan, E., Vanhoucke, V. Rabinovich, A. (2015), "Going deeper with convolutions”. IEEE Conference on Computer Vision and Pattern Recognition (CVPR), pp. 1-9 [54] Krizhevsky, A., Sutskever I., Hinton, G. (2012). “ImageNet classification with deep convolutional neural networks”, Proceedings of the 25th International Conference on Neural Information Processing Systems, Lake Tahoe, Nevada, USA - Volume 1, pp. 1097–1105 [55] He, k., Zhang, X., Ren, S. Sun, J. (2016). “Deep Residual Learning for Image Recognition”, IEEE Conference on Computer Vision and Pattern Recognition (CVPR). Las Vegas, NV, USA, 27 June 2016. IEEE, pp. 770-778 [56] Ronneberger, O., Fischer, P. Brox. (2015). “U-Net: Convolutional Networks for Biomedical Image Segmentation”, Medical Image Computing and Computer-Assisted Intervention - MICCAI 2015. Munich, 5-9 October 2015. Springer International, pp. 241-241 [57] Kolhar, S. Jagtap, J. (2021). “Convolutional neural network based encoder-decoder architectures for semantic segmentation of plants”. Ecological Informatics, 64, pp. 1-12 [58] Kolesnikov, A., Dosovitskiy, A. Weissenborn, D., Heigold, G., Uszkoreit, J., Beyer, L., Minderer, M., Dehghani, M., Houlsby, N., Gelly, S., Unterthiner, T., Zhai, X. (2021). “An Image is Worth 16x16 Words: Transformers for Image Recognition at Scale”. ICLR 2021 [59] Yu, F. Koltun, V. (2016). “MULTI-SCALE CONTEXT AGGREGATION BY DILATED CONVOLUTIONS”. ICLR 2016 [60] Ren, S., He, K., Girshick, R. Sun, J. (2017). “Faster R-CNN: Towards Real-Time Object Detection with Region Proposal Networks”. In: IEEE Transactions on Pattern Analysis and Machine Intelligence. 1 June 2017. IEEE, 39(6), pp. 1137-1149. DOI: 10.1109/TPAMI.2016.2577031. [61] Goodfellow, I., Pouget-Abadie, J., Mirza, M., Xu, B., Warde-Farley, D., Ozair, S., Bengio, Y. (2014). “Generative adversarial nets”. In: Advances in neural information processing systems, pp. 2672–2680 [62] hlamba28 (2019). UNET-TGS. Avaliable in: https://github.com/hlamba28/UNET-TGS/blob/master/ TGS%20UNET.ipynb [Consulted 08-16-2022]. [63] Tensorflow (2022). tf.keras.layers.Conv2DTranspose. Avaliable in: https://www.tensorflow.org /api docs/python/tf/keras/layers/Conv2DTranspose [Consulted 08-16-2022]. [64] Brownlee, J. (2018). A Gentle Introduction to Dropout for Regularizing Deep Neural Networks. Avaliable in: https://machinelearningmastery.com/dropout-for-regularizing-deep-neural-networks/ [Consulted 08-18- 2022]. [65] Artificial Intelligence AWS (2021). What is the difference between same convolution and full convolution in terms of feature map size?. Avaliable in: https://ai.stackexchange.com/questions/27779/what-is-the- difference-between-same-convolution-and-full-convolution-in-terms-of [Consulted 08-16-2022]. [66] Ferro, C., Mayorga, C. N. García, A. (2020). Llenado de series de datos de 2014 a 2019 de PM2.5 por medio de una red neuronal y una regresión lineal. Thesis. Universidad Sergio Arboleda Madrid. DOI: 10.13140/RG.2.2.35092.53126. [67] Keras (2020). Callbacks API. Avaliable in: https://keras.io/api/callbacks/ [Consulted 08-17-2022]. [68] Godoi, D. (2018). Understanding binary cross-entropy / log loss: a visual explanation. Avaliable in: https://towardsdatascience.com/understanding-binary-cross-entropy-log-loss-a-visual-explanation- a3ac6025181a [Consulted 08-30-2022]. [69] ML Cheatsheet (2017). Loss Function. Avaliable in: https://ml-cheatsheet.readthedocs.io/en/latest/ loss-functions.html [Consulted 08-30-2022]. 131 [70] Keras (2022). Accuracy. Avaliable in: https://keras.io/api/metrics/accuracy metrics/ [Consulted 08-30- 2022]. [71] Keras (2022). Tune hyperparameters in your custom training loop. Avaliable in: https://keras.io/guides/ keras-tuner/custom-tuner/ [Consulted 08-30-2022]. [72] Wadekar, S. (2021). Hyperparameter Tuning in Keras. Avaliable in: https://medium.com/swlh/hyperparameter-tuning-in-keras-tensorflow-2-with-keras-tuner-randomsearch- hyperband-3e212647778f [Consulted 08-30-2022]. [73] Wolfram Mathworld (2022). Sigmoid function. Avaliable in: https://mathworld.wolfram.com/SigmoidFunction.html [Consulted 08-30-2022]. [74] zeka0 (2017). Custom layer: trainable threshold. Avaliable in: https://github.com/keras-team/keras/issues/ 6926 [Consulted 09-01-2022]. [75] Katherine, L. (2022). How to Choose a Learning Rate Scheduler for Neural Networks. Avaliable in: https://neptune.ai/blog/how-to-choose-a-learning-rate-scheduler [Consulted 09-01-2022]. [76] Vasani, D. (2019). This thing called Weigth Decay. Avaliable in: https://towardsdatascience.com/thisthing-called-weight-decay-a7cd4bcfccab [Consulted 09-01-2022]. [77] Otsu, N. (1979). “A Threshold Selection Method from Gray-Level Histogram”. IEEE Transactions on System Man Cybernetics, 9, pp. 62-66. DOI: http://dx.doi.org/10.1109/TSMC.1979.4310076. [78] OpenCV (2022). Miscellaneous Image Transformations. Avaliable in: https://docs.opencv.org/4.x/d7/d1b/ group imgproc misc.html-ggaa42a3e6ef26247da787bf34030ed772cad0c5199ae8637a6b195062fea4789fa9 [Consulted 09-01-2022]. [79] Massieu, W. (2020). CENTROS DE MASA, DE GRAVEDAD Y CENTROIDES. Avaliable in: https://www.ipn.mx/ assets/files/cecyt11/docs/Practicas/Fisica/Fisica1/practica-7.pdf [Consulted 09-01-2022]. [80] Scipy docs (2022). Scipy ndimage label. Avaliable in: https://docs.scipy.org/doc/scipy/reference/ generated/scipy.ndimage.label.html [Consulted 09-01-2022]. [81] Scipy docs (2022). Scipy ndimage measurements center of mass. Avaliable in: https://docs.scipy.org/doc/ scipy-0.14.0/reference/generated/scipy.ndimage.measurements.center-of-mass.html [Consulted 09-01-2022]. [82] StackOverflow (2018). Get coordinates of local maxima in 2D array above certain value. Avaliable in: https://stackoverflow.com/questions/9111711/get-coordinates-of-local-maxima-in-2d-array-above-certain- value [Consulted 09-01-2022]. [83] Scipy docs (2022). Scipy spatial KDTREE. Avaliable in:: https://docs.scipy.org/doc /scipy/reference/- generated/scipy.spatial.KDTree.html [Consulted 09-01-2022]. [84] Murillo, B. E. (1660-1665). The farewell of the prodigal son, Museo Nacional del Prado (P998). [85] Murillo, B. E. (1660-1665). The dissipation of the prodigal son, Museo Nacional del Prado (P999). [86] Montes, D. (2018). Diagrama BoxPlot. Avaliable in: https://www.pgconocimiento.com/diagramaboxplot/ [Consulted 09-01-2022]. [87] National Institute of Health (2021). ImageJ. Avaliable in: https://imagej.nih.gov/ij/ [Consulted 08-18- 2022]. [88] Python Docs (2022). Multiprocessing. Avaliable in: https://docs.python.org/3/library/multiprocessing.html [Consulted 17-08-2022]. [89] Python Docs (2022). HDF5 for Python. Avaliable in: https://docs.h5py.org/en/stable/ . [Consulted 08-21-2022]. 132 Chapter 7. References [90] Paret, L. (1777). Self-portrait in the studio, Museo Nacional del Prado (P007701). [91] Velázquez, D. (1630). The forge of Vulcan, Museo Nacional del Prado (P001171). [92] Rubens, P.P. (1630-1644). The abduction of Europa, Museo Nacional del Prado (P001693). [93] Cabrera, M. (1751). Presentation of the Virgin in the Temple, Museo Nacional del Prado (P49435-1). [94] Cabrera, M. (1751). The Annunciation, Museo Nacional del Prado (P49435-2). List of Figures 1.1 Types of fabrics and frame 3 1.2 Plain weave - Taffeta [10] 3 1.3 Parts of an old loom [11] 4 1.4 Painting and X-ray plate of Adan y Eva by Rubens 5 1.5 Vertical density map of Adan y Eva by Rubens 5 2.1 X-ray of P1114 Ribera (MNP) [34] 8 2.2 Spectrum of a 1cm patch from a canvas 8 2.3 Scenario where FT works well 10 2.4 Scenario where FT doesn’t work well 10 2.5 PSD and triangle for thread counting [17] 11 2.6 Aracne interface [18] 12 3.1 X-ray detail of P001692 Rubens (MNP) [38] 16 3.2 X-ray detail of P001180 Velázquez (MNP) [39] 16 3.3 Neuron [43] 17 3.4 Activation function [43] 17 3.5 Feed-forward Network [43] 18 3.6 CNN + Fully-connected ANN [43] 19 3.7 Subsampling 2x2 (MaxPool) in CNN [43] 19 3.8 Performance of gradient descent on a surface representing the error [43] 20 3.9 Training scheme of a network [43] 22 3.10 A patch and its label 22 3.11 Example of the variety present in the dataset 23 3.12 Obtaining 4 patches from the corners of the original image 24 3.13 300x300 pixel image flipped from right to left 25 3.14 300x300 pixel image flipped top to bottom 25 3.15 Obtaining the two central patches that will be rotated later 26 3.16 Transposed of the three 300x300 pixel images 26 3.17 Distribution of thread densities per cm in the three subsets 27 3.18 Difference between all the samples obtained from the same image using the updated Data Augmentation 28 3.19 Difference between all the samples extracted from an image having all of them rotated 90º 28 3.20 U-Net architecture [49] 30 3.21 Variation of dimensions in a U-Net [49] 31 3.22 SegNet architecture [57] 31 3.23 Activation functions used [66] 32 3.24 conv2dblock [62] 33 133 134 List of Figures 3.25 Architecture used with conv2dblocks modules in blue horizontal arrows, with the number of filters used, ni , indicated on top. Up and down steps involve halving reductions of sizes downwards, and duplicate them upwards, with zero padding if needed. Upcoming features are concatenated with the features at the same level in the encoder. The size of all feature maps is indicated below them. At the output a 1×1convolutional filter with sigmoid function is used 33 3.26 Cross-entropy or log loss [69] 34 3.27 Learning curves of the two best training sessions 35 3.28 Learning curves of the two best trainings with LR 1e−435 3.29 Segmentation of validation images: 36 3.30 Learning curve of the best model training getUnet6 with LR 1e−337 3.31 Segmentation of validation images: Training 1 with LR 1e−3getUnet6 38 3.32 Segmentation of validation images with getUnet KT and getUnet KT CS 41 4.1 Sigmoid activation function [73] 43 4.2 Example of segmentation provided by the U-Net 44 4.3 Test image, binarized segmentation and found regions 46 4.4 Test image and grayscale segmentation 46 4.5 Output of the local maxima filter 46 4.6 Otsu thresholding 47 4.7 Binary image where the maximum filter points coincide with those of the original image 47 4.8 Image with the located points 47 4.9 Location of points using local maxima 48 4.10 Graphical representation of the spatial count method 49 4.11 Box diagram of error when Otsu binarization and spatial count are used 52 4.12 Box diagram of error when Local MaxFilter binarization and spatial count are used 53 4.13 Box diagram of error when Otsu binarization and frequencial count are used 54 4.14 Error in each image using 3rd training of getUnet KT model 55 4.15 Test images where error bigger than 8% occurs using getUnet KT model 56 4.16 Validation images where error bigger than 12% occurs using getUnet KT model 57 4.17 Error in each image using 3rd training of getUnet KT Threshold model 58 4.18 Test images where error bigger than 8% occurs using getUnet KT CS model 59 4.19 Validation images where error bigger than 12% occurs using getUnet KT CS model 60 5.1 Plate P0998 shown with ImageJ. In red, the region of interest 64 5.2 Division in windows of plate P998 by Murillo 66 5.3 P07905 plate and example patch 70 5.4 P07906 plate and example patch 71 5.5 P001180 plate and example patch 71 5.6 P998 plate and example patch 72 5.7 P999 plate and example patch 72 5.8 P001692 plate and example patch 73 5.9 P001693 plate and example patch 73 5.10 P1114 plate and example patch 74 5.11 P49435-1 plate and example patch 74 5.12 P49435-2 plate and example patch 75 5.13 P07905 density maps performing a frequential analysis of the plate 76 5.14 P07905 density maps performing a frequential analysis of the segmentation with getUnet KT 76 5.15 P07905 density maps performing a frequential analysis of the segmentation with getUnet KT CS 76 5.16 P07905 density maps performing a spatial analysis of the segmentation with getUnet KT 77 5.17 P07905 density maps performing a spatial analysis of the segmentation with getUnet KT CS 77 5.18 P07905 Angle deviation in degrees with getUnet KT /getUnet CS 77 5.19 P07905 Counting histograms with getUnet KT 78 5.20 P07905 Counting histograms with getUnet KT CS 78 5.21 P07905 examples 79 5.22 P07906 density maps performing a frequential analysis of the plate 80 5.23 P07906 density maps performing a frequential analysis of the segmentation with getUnet KT 80 5.24 P07906 density maps performing a frequential analysis of the segmentation with getUnet KT CS 81 List of Figures 135 5.25 P07906 density maps performing a spatial analysis of the segmentation with getUnet KT 81 5.26 P07906 density maps performing a spatial analysis of the segmentation with getUnet KT CS 81 5.27 P07906 Angle deviation in degrees with getUnet KT /getUnet KT CS 82 5.28 P07906 Counting histograms with getUnet KT 82 5.29 P07906 Counting histograms with getUnet KT CS 83 5.30 P07906 examples 83 5.31 P001180 density maps performing a frequential analysis of the plate 84 5.32 P001180 density maps performing a frequential analysis of the segmentation 85 5.33 P001180 density maps performing a spatial analysis of the segmentation 86 5.34 P001180 Angle deviation in degrees with getUnet KT /getUnet KT CS 87 5.35 P001180 Counting histograms with getUnet KT 88 5.36 P001180 Counting histograms with getUnet KT CS 88 5.37 P001180 examples 89 5.38 P998 density maps performing a frequential analysis of the plate 90 5.39 P998 density maps performing a frequential analysis of the segmentation with getUnet KT 90 5.40 P998 density maps performing a frequential analysis of the segmentation with getUnet KT CS 90 5.41 P998 density maps performing a spatial analysis of the segmentation with getUnet KT 91 5.42 P998 density maps performing a spatial analysis of the segmentation with getUnet KT CS 91 5.43 P998 Angle deviation in degrees with getUnet KT /getUnetKT CS 91 5.44 P998 Counting histograms with getUnet KT 92 5.45 P998 Counting histograms with getUnet KT CS 92 5.46 P998 examples 93 5.47 P999 density maps performing a frequential analysis of the plate 94 5.48 P999 density maps performing a frequential analysis of the segmentation with getUnet KT 94 5.49 P999 density maps performing a frequential analysis of the segmentation with getUnet KT CS 94 5.50 P999 density maps performing a spatial analysis of the segmentation with getUnet KT 95 5.51 P999 density maps performing a spatial analysis of the segmentation with getUnet KT CS 95 5.52 P999 Angle deviation in degrees with getUnet KT /getUnet KT CS 95 5.53 P999 Counting histograms with getUnet KT 96 5.54 P999 Counting histograms with getUnet KT CS 96 5.55 P999 examples 97 5.56 P001692 density maps performing a frequential analysis of the plate 98 5.57 P001692 density maps performing a frequential analysis of the segmentation with getUnet KT 98 5.58 P001692 density maps performing a frequential analysis of the segmentation with getUnet KT CS 99 5.59 P001692 density maps performing a spatial analysis of the segmentation with getUnet KT 99 5.60 P001692 density maps performing a spatial analysis of the segmentation with getUnet KT CS 100 5.61 P001692 Angle deviation in degrees with getUnet KT /getUnet KT CS 100 5.62 P001692 Counting histograms with getUnet KT 101 5.63 P001692 Counting histograms with getUnet KT CS 101 5.64 P001692 examples 102 5.65 P001693 density maps performing a frequential analysis of the plate 103 5.66 P001693 density maps performing a frequential analysis of the segmentation with getUnet KT 103 5.67 P001693 density maps performing a frequential analysis of the segmentation with getUnet KT CS 104 5.68 P001693 density maps performing a spatial analysis of the segmentation with getUnet KT 104 5.69 P001693 density maps performing a spatial analysis of the segmentation with getUnet KT CS 105 5.70 P001693 Angle deviation map with getUnet KT /getUnet CS 105 5.71 P001693 Counting histograms with getUnet KT 106 5.72 P001693 Counting histograms with getUnet KT CS 106 5.73 P001693 examples 107 5.74 P1114 density maps performing a frequential analysis of the plate 108 5.75 P1114 density maps performing a frequential analysis of the segmentation with getUnet KT 108 5.76 P1114 density maps performing a frequential analysis of the segmentation with getUnet KT CS 108 5.77 P1114 density maps performing a spatial analysis of the segmentation with getUnet KT 109 5.78 P1114 density maps performing a spatial analysis of the segmentation with getUnet KT CS 109 5.79 P1114 Angle estimation with getUnet KT /getUnet CS 109 5.80 P1114 Counting histograms with getUnet KT 110 5.81 P1114 Counting histograms with getUnet KT CS 110 136 List of Figures 5.82 P1114 Examples 111 5.83 P49435-1 density maps performing a frequential analysis of the plate 112 5.84 P49435-1 density maps performing a frequential analysis of the segmentation with getUnet KT 112 5.85 P49435-1 density maps performing a frequential analysis of the segmentation with getUnet KT CS 113 5.86 P49435-1 density maps performing a spatial analysis of the segmentation with getUnet KT 113 5.87 P49435-1 density maps performing a spatial analysis of the segmentation with getUnet KT CS 114 5.88 P49435-1 Angle estimation with getUnet KT /getUnet KT CS 114 5.89 P49435-1 Counting histograms with getUnet KT 115 5.90 P49435-1 Counting histograms with getUnet KT CS 115 5.91 P49435-1 examples with getUnet KT 116 5.92 P49435-2 density maps performing a frequential analysis of the plate 117 5.93 P49435-2 density maps performing a frequential analysis of the segmentation with getUnet KT 117 5.94 P49435-2 density maps performing a frequential analysis of the segmentation with getUnet KT CS 118 5.95 P49435-2 density maps performing a spatial analysis of the segmentation with getUnet KT 118 5.96 P49435-2 density maps performing a spatial analysis of the segmentation with getUnet KT CS 119 5.97 P49435-2 Angle deviation in degrees with getUnet KT /getUnet KT CS 119 5.98 P49435-2 Counting histograms with getUnet KT 120 5.99 P49435-2 Counting histograms with getUnet KT CS 120 5.100 P49435-2 examples with getUnet KT 121 5.101 P998 and P999 match 122 5.102 P001692 and P001693 vertical match 123 5.103 P49435-1 and P49435-2 vertical match 124 List of Tables 3.1 getUnet training results 34 3.2 getUnet training results, LR 1e−435 3.3 getUnet6 training results 37 3.4 getUnetKT training results 40 3.5 getUnetKT CS training results 40 5.1 Execution time of each version of the algorithm 69 5.2 Execution times of the multiprocessing algorithm on different machines 69 137