Keywords:

List of abbreviations

Acknowledgment

I have the pleasure of modestly including in this section a mention to all the people that made this work possible, and supported me in this academic adventure. Without them, none of this would have been accomplished.

First and mostly to my thesis director, Professor at the ARTEMIS Department of Télécom SudParis, Catalin FETITA. His caring and exhaustive dedication motivated me to achieve a methodological and innovative study. My deep gratitude to him for teaching me how to conduct a professional research, and allowing me to grow both professionally and personally.

I greatly thank Prof. Eric PETIT and Prof. Bertrand THIRION for having accepted reviewing this manuscript and feeding back valuable and meaningful considerations. To Prof. Philippe GRENIER for his participation as President of the jury, for his clarity and problem understanding that always inspired me, as well to Prof. Julie DELON for her thoughtful analysis of this work and her constructive remarks.

I would also like to include a special acknowledge to Prof. Pierre-Yves BRILLET, Radiologist at the Avicenne Hospital in Bobigny for his great predisposition to work together, merging radiology and engineering disciplines. His professionalism to tackle questions during the study helped in achieving better practical results. I thank him for his commitment to transmit (and translate) the key medical concepts to be considered, and to have provided a valuable annotated database, a key element of this work. I would like to mention all ARTEMIS department who received me as a member from day one, from the Director of the department Titus ZAHARIA, to all the associated professors. Special mention to Evelyne TARONI for carefully taking care of all the administrative issues during all these years. To all my fellow Ph.D candidate colleagues (Florence, Sébastien, Dancho, Nicolas) with whom I had the pleasure to work and share ideas, thoughts, and feelings. For enriching my creativity and soul, I am grateful to them.

An exceptional mention goes to my family to whom I owe any achievement and progress in my life. To my dear parents Graciela and Roberto for believing in me and supporting me unconditionally with their love and care. And to my beloved sister and brother, Nati y Mati, who inspire me with their extraordinary aptitudes and kindness, and make me become a better person. Even far away, you are always next to me helping carrying on with my dreams.

No phrase will accomplish to thank my dear love Isabel. Always showing me the best human quality I have ever found anywhere, it is a life gift being accompanied by such a lovely person with such a good influence. Thank you for every moment, and being there also in those difficult ones.

Finally, I would like to thank the people of the Maison de l'Argentine (Cité Universitaire) who facilitated my residence in Paris and made me feel at home: Marcelo Balsells, Julia Bernardi, Agustina Primo, Lorena Gabelloni, Madeleine. To them, for having achieved making the place inclusive and harmonious. And to my friends I made there these years, with whom I have shared wonderful and memorable experiences during this work, and hopefully after. You will be always part of my thoughts and heart. By mentioning Fede Ventosinos, who symbolizes due to his kind heart the friendship of our group, I consider you have all been dearly mentioned.

Résumé

Les pathologies infiltrantes diffuses recensent un large groupe de désordres pulmonaires et nécessitent un suivi régulier en imagerie tomodensitométrique (TDM). Une évaluation quantitative est nécessaire pour établir la progression (régionale) de la maladie et/ou l'impact thérapeutique. Cela implique le développement d'outils automatiques de diagnostic assisté par ordinateur (DAO) pour la segmentation du tissu pathologique dans les images TDM, problème adressé comme classification de texture.

Traditionnellement, une telle classification repose sur une analyse des caractéristiques texturales 2D dans les images TDM axiales selon des critères définis par l'utilisateur. Récemment, des techniques d'intelligence artificielle fondées sur l'apprentissage profond, notamment les réseaux neuronaux convolutionnels (CNN), ont démontré des performances meilleures pour résoudre des tâches visuelles. Toutefois, pour les architectures CNN « classiques » il a été prouvé que les performances étaient moins bonnes en classification de texture par rapport à la reconnaissance d'objets, en raison de la dimensionnalité intrinsèque élevée des données texturales.

Dans ce contexte, ce travail propose un système automatique pour l'analyse quantitative des pathologies infiltrantes diffuses du poumon fondé sur une architecture CNN en cascade (conçue spécialement pour l'analyse de texture) et sur un prétraitement spécifique des données d'entrée par filtrage localement connexe (permettant d'atténuer l'intensité des vaisseaux pulmonaires et d'augmenter ainsi le contraste des régions pathologiques). La classification, s'appliquant à l'ensemble du volume pulmonaire, atteint une précision moyenne de 84% (75.8% pour le tissu normal, 90% pour l'emphysème et la fibrose, 81.5% pour le verre dépoli).

Presentation overview

Infiltrative lung diseases (ILDs) enclose a large group of irreversible lung disorders which require regular follow-up with computed tomography (CT) imaging. ILDs are very rare among lung diseases, and include around 300 different phenotypes. Idiopathic Pulmonary Fibrosis (IPF) is the most spread and responsible of nearly 5% of deaths (per year) from lung disease. The attention to such diseases is nowadays motivated by the advent of new therapies and drugs requiring clinical evaluation. The assessment of ILDs relies nowadays on CT imaging modality because of its high spatial resolution and a good lung tissue contrast provided. This avoids surgical intervention for diagnosis and allows assessing the disease progression over time to guide the therapy.

Although there is a radiological consensus achieved on CT image biomarkers that have to be searched for in diagnosing IPFs, image interpretation remains difficult due to two main issues: multiple sources of variability in CT acquisition altering the image appearance, and the ontological overlap between defined phenotypes or classes. A detailed anatomical description of the targeted pathologies will be presented (chapter 1) to show the significant expertise required to achieve a correct diagnosis. This medical context implies the development of automated computer-aided diagnosis (CAD) tools for pathological lung tissue segmentation. In this way, the pathological regions would be objectively and robustly segmented and fully quantified per CT volume.

The latest state of the art techniques dealing with ILD patterns categorization are thus studied in chapter 2. This problem is addressed as pixel-based texture classification. Traditionally, such classification relies on a two-dimensional analysis of axial CT images by means of handcrafted features. Here will be detailed the three main modules that integrate such a classification model: space transformation, feature representation, and classification. A comparative analysis of the best performing methods, focusing on lung texture categorization completes the presentation. We point out that the use of deep learning techniques, especially Convolutional Neural Networks (CNNs) for visual tasks, have shown great improvements with respect to handcrafted heuristics-based methods, and by consequence, this type of approach becomes the selected reference technique for this work.

Our contribution in the development of a CNN-based CAD system for ILD segmentation and quantification of CT volumes is detailed in chapter 3. An exhaustive description of the proprietary database, built (in collaboration with Avicenne Hospital, Bobigny, France) for the specific task during this research, receives a special attention as one of the key elements of this study. Two preprocessing steps are introduced that need to be applied on the experiments: lung field segmentation, and patch filling procedure (to regularize the subpleural region classification). Afterwards, a hybrid method combining an existing parametric approach and a CNN trained on a reduced dataset is studied as an attempt to achieve higher robustness and specificity.

The CNN showed a good potential in regularizing the classification output, but its usage as unique classifier encountered several issues: necessity to extend the training database, need to reconsider the CNN architecture to target texture instead of object recognition, and the need to take into account the bias introduced by normal high densities of vascular tree and thorax cage. Our work considered these aspects to achieve a performant computer-aided diagnosis system of ILDs based on CNNs. The final classification system adopted is then presented and explained. It consists of a preprocessing module performing a locally connected filtering (LCF) followed by a cascade of two CNNs for a hierarchical classification. The specific preprocessing is applied to the lung images to segment and attenuate the vessel densities while preserving high opacities related to pathologies. This also allows completing the structural knowledge related to the lung mask and airways, by the vascular network.

Locally connected filters (LCF) exploit grayscale topological connectivity of each point on the support of a function 𝑓 to filter out 𝑓 values weakly connected with a pre-defined subset 𝑌 𝑠𝑢𝑝𝑝 𝑓 in the local environment of that point. The topological connectivity is verified by means of a morphological filter, namely the grayscale reconstruction by dilation 𝑅 𝑓  (., 𝑌), since the targeted structures here are the blood vessels, exhibiting a positive contrast with respect to their environment. The shape and size of the reference subset Y will determine the type of the researched connectivity. LCF effect is to attenuate (or suppress) the f values which are not "linked" with their kdistant neighborhood via a high-intensity path. On contrary, when such connection exists, the structures are preserved via the reconstruction operator with a slight "flattening" of the grayscale levels. LCF thus not only presents a denoising property similar to the median filter, but also preserves spatial structures selected by a local connectivity configuration. To remove vascular structures in the 3D image datasets, we take into account the linear shape of the vessels and their cross-sectional size, and set-up a filtering scheme combining 2D LCF applied across planes oriented orthogonal to different directions in space. We chose 9 spatial directions corresponding with the 18-connectivity (excluding symmetric orientations). The vascular structures are detected by adaptive thresholding and directional reconstruction by erosion applied to the difference between 𝑓 (original) and the LCF filtered image. The sheet-like structures are minimally affected by the filter. To remove vascular structures of different sizes, the LCF is applied in a multiresolution scheme using 2 levels of decimation. In addition, the vascular structures are extracted at each level of resolution and combined together prior to a final adaptive thresholding and filtering, which selects the highest confidence vessels. The vascular network thus extracted is added to the lung mask obtained (and completed by the airway structures) to reinforce the train-and-test process of the network.

With respect to the cascade of CNNs, the state of the art T-CNN network architecture is first used to discriminate the less correlated classes, by grouping fibrosis and ground glass in one single class. The latter two will be further differentiated using a shallow network architecture with random filters (referred to as Rand-CNN in the following) inspired from a texture synthesis approach.

The T-CNN network consists of 5 cascaded layers of convolutional filters to compute image features, with 2x2 kernels and LeakyReLU activations, followed by just one average pooling, with size equal to the size of final feature maps and three dense layers that act as the classifier. Unlike of standard CNN architectures adopted in the literature, this T-CNN does not integrate pooling layers operations between the convolutional layers, and also uses the smallest kernel size to help extracting the best texture representation. The T-CNN ends with a final 3-way SoftMax layer for pathological or normal classification. We used 49 patients for training (from scratch) the T-CNN with the purpose of reaching a numerous and representative set required by the CNN to perform a valid generalization. All these cases are not considered for testing. The rest of 10 different DICOM image series of whole thorax scans are used as test cases. From the training set, non-overlapping patches of 32x32 pixels were extracted. We consider those patches falling 80% inside the annotated ROIs. Horizontal flipping and rotations are applied in order to artificially increase the number of samples and avoid over-fitting the neural network. Thus, the training dataset consisted of 24265 image patches: 8978 normal, 8871 fibrosis+ground glass, 6416 emphysema. Similarly, 32x32 pixels patches were extracted from the test dataset as input for the T-CNN, centered to each pixel in the test image to which a class will be assigned with a given probability.

The second module of the cascade, Rand-CNN, tackles the subclassification between ground glass and fibrosis by developing a new texture descriptor. Basically, a texture is characterized by the arrangement of local patterns. It was demonstrated that CNNs are able to extract not only the content of an image but also the style representation, or the texture. The second one, of our interest, was achieved by computing the Gram matrix (correlation between the feature maps given a depth of the CNN). By including the feature correlations of multiple layers, a stationary, multi-scale representation of the input image can be obtained, which captures its texture information but not the global spatial arrangement. It was also shown that even with shallow CNNs with random filters, the style representation could be found for a given input image. These works use such information to synthesize texture images based on examples. In our case, we will use this texture space as a feature descriptor of the pathologies under analysis. Merging these concepts, we use a Random-Multiscale network with seven different filter sizes 𝑛 × 𝑛 with 𝑛 = 3; 5; 7; 11; 15; 23; 29, and 128 feature maps each (896 feature maps in total) to transform the input data. Filter coefficients are obtained from a uniform distribution. Notice that no training is needed in contrast with the T-CNN. Instead, the Gram matrix is computed for the resulting feature maps. The element-wise mean and the standard deviation for the whole set of Gram matrices obtained from T-CNN training patches are computed to obtain two signature matrices (𝐸 𝑘 , and 𝑆 𝑘 respectively, 𝑘 = {1,2}) for each of the two classes to be detected. Finally, the classification is made by calculating the discrepancy δ 𝑘 between the Gram matrix of the test patch (𝐺 𝑋 ) and the signature Gram matrixes of the two classes -fibrosis and ground glass. The minimum value for this computation indicates higher resemblance between textures and dictates the decision making.

The 3D test database not being fully annotated, a fully-quantitative evaluation of the classification result was not possible. A semi-quantitative assessment was instead performed using sparse annotations of the ground truth (10 axial images per scan, evenly distributed). The average values obtained for the four defined classes (normal, emphysema, fibrosis and ground glass) were as follows: sensitivity: 51.6%, specificity: 86.4%, accuracy: 84.6%, and precision: 41.2%. If we analyze the normalized confusion matrix for the whole test database, we notice a good prediction of the proposed method for normal and emphysema, but also a less discriminative power between fibrosis and ground glass patches. Several explanations to this have been advanced. One may come from the effect of LCF filtering which also "flattens" the grayscale variation of the fibrosis regions making them more similar to ground glass (especially for the most confusing situations). Another may come from an effect of fuzzy decision at the interface between fibrosis and normal areas (patches falling on these interface zones will be more affected by a misclassification). Also, low precision for some classes (mainly emphysema), may be explained due to the effect of gravitational gradient which modifies the lung texture appearance in the anterior part of the lung (supine patient position). Finally, a lung border effect is acknowledged: few high intensity pixels from thorax cage may be included in the lung field (due to mask regularization) which bias the decision in subpleural areas (effect amplified by the patch filling procedure).

In addition, the manuscript proposes a discussion about the usefulness of each component of the overall developed system. The benefit of using the cascade of CNNs is first demonstrated with respect to the T-CNN alone for the same classes both qualitatively and quantitatively (an increase of 10% in sensitivity, 3.5% in precision, 1.2% in accuracy, and same behavior for specificity). Then, we have investigated the benefit of the prefiltering module in the proposed framework by assessing the system performance with and without LCF on the same training and validation dataset. Overall, an increase of nearly 10% in sensitivity, 3.5% in specificity and 4.5% in accuracy was assessed when using LCF versus no prefiltering. Finally, we reduced the size of the context patch from 32x32 pixels to 16x16 to study its influence on the classification output. The same behavior was observed for both patch sizes with a better spatial sensibility for the 16x16 pixels patch (as the patch is smaller, it is less biased at interfaces, leading to sharper transition between two different regions).

A conclusion chapter summarizes the contributions of this work and defines future directions for improvement to allow the translation of this research into clinical routine.

Context of the study 1.1 Description of Purpose

This research aims at developing a new robust and reliable solution to be implemented in a computer-aided diagnosis system for quantitative follow-up of pulmonary diseases using volumetric computed tomography.

Among the medical imaging acquisition methods available nowadays, Multi Slice Computed Tomography (MSCT) is currently the most widely available and accurate to characterize and quantify lung pathologies subtypes in a manner that can help determine the progression of specific patterns of disease over time. It is the reference imaging modality for the lung morpho-pathological assessment due to its high spatial resolution and the ability to depict air-filled regions for which other imaging modalities such as magnetic resonance imaging or echography fail.

With this consideration, the proposed solution involves methods for body regions segmentation from whole body Computed Tomography (CT) image series. The aim is to provide a Region of Interest (ROI) partitioning of the dataset according to the anatomical and/or identified pathological regions, mainly interstitial lung diseases (ILDs) and emphysema. Additional tools for computer-aided diagnosis will be considered when applicable, i.e., lung segmentation algorithms. Such ROI partitioning will be exploited to measure longitudinal changes between different temporal series (patient follow-up).

To understand the motivation of this work and the choices that have been made hereafter, we briefly introduce in the following the framework of the study in terms of societal issues, and the actual radiology workflow deployed for diagnosis and follow-up of lung diseases.

Background: social impact of lung diseases

This research focuses on the study of a specific category of non-neoplasic lung pathologies including emphysema and Interstitial Lung Diseases (ILDs), and especially on Idiopathic Interstitial Pneumonias (IIPs) [START_REF]American Thoracic Society/European Respiratory Society international multidisciplinary consensus classification of the idiopathic interstitial pneumonias[END_REF] -a variety of ILDs. It aims at developing a computer-aided diagnosis system of these pathologies exploiting the MSCT imaging technology.

The motivation of such a research comes from the increasing healthcare burden associated with this category of chronic diseases and from an unmet clinical need of efficient diagnosis and follow-up tools. Figure 1 summarizes some basic statistics established in 2013 about ILDs impact on European society [66].

In United States, 3.7 million people have been diagnosed with emphysema, of which 16700 deaths registered, whereas other studies reported that 80.9 per 100000 men and 67.2 per 100000 women suffer from ILDs with 31.5 and 26.1 new cases diagnosed per 100000 persons per year, respectively.

It is important to mention that emphysema and IIPs are irreversible diseases that can be slowed but not stopped, and to date, no efficient treatment is available for most of IIPs. Such diseases affect the tissue and space within and around the alveoli (air sacs) and sometimes also the airways, the blood vessels and the pleura, leading to profound impairment in lung physiology. Unfortunately, if the initiated injury or abnormal repair from injury is not halted, progressive tissue damage can lead to worsening physiologic impairment and even death. In this context, early diagnosis and follow-up of these pathologies are the key issues for increasing the quality of life and the survival chances of the patient.

If the suspicion of lung disease can be clinically confirmed or rejected by performing pulmonary function tests and chest radiography, to date only the computed tomography (CT) investigation allows an accurate diagnosis of the type of pathology (IIPs including more than 200 different disease entities). In addition, the visual feed-back provided by CT is very valuable due to the high spatial resolution of this imaging technique which allows accurately depiction and spatial localization of the disease patterns. In this respect, an overview of the radiology workflow based on CT images investigation is presented in the next section.

Radiology workflow: practical scope and limitations

In this section, the evolving challenges of radiologists in the visualization and interpretation of large complex imaging studies in the present day is examined, in connection with advances in information technologies. It will serve as contextual information to present the study case of this work and its motivation.

Modern medical imaging modalities, such as Computed Tomography (CT), generate large data sets that are difficult and time-consuming to review by using the standard axial section view. The radiologist must navigate through hundreds or thousands of images (one thorax CT acquisition produces around 700-1000 images). Image properties can be adjusted through different display settings for qualitative analysis and post processed to include new auxiliary tools to more easily and accurately extract quantitative information from the image data, including morphologic and physiologic metrics for example.

The image series are viewed on computer workstations taking advantage of the inherent digital characteristics of the data, after a given patient record is requested to the storage center (see Figure 2 for a schematic illustration of the radiology workflow). Both the client workstations and the center of storage are part of the hospital PACS (picture archiving and communication system), a healthcare technology for the storage, retrieval, management, distribution and presentation of medical images. A PACS has four major components: imaging systems, such as our reference technique fur lung examination, computed axial tomography (CAT scan); a secure network for distribution and exchange of patient information; workstations or mobile devices for viewing, processing and interpreting images; archives for storage and retrieval of images and related documentation and reports. Note that images are generally still displayed as axial sections. As human observers, radiologists must search for, detect and interpret targets (task known as the interaction of human observers with complex data sets), tools for upgrading the current workflow should be based on human perceptual and attentional abilities and limitations. As the interpretive process involves multiple complex steps, only the essentials steps based on the description proposed in [START_REF] Andriole | Optimizing analysis, visualization, and navigation of large image data sets: one 5000-section CT scan can ruin your whole day[END_REF] are resumed here:

• Search, detection and localization

When interpreting a diagnostic CT imaging, the radiologist's first task is the detection and localization of any potential abnormalities. Human performance limitations lead to errors in this task,

Image acquisition

Image interpretation

Report Medical decision in the form of false-negative or false-positive judgments. Lesions can be missed, the target may not be salient enough for the observer to recognize it as an abnormality, or it can be consciously reported as an artifact or a normal structure (artifacts and other challenges that arise with CT imaging are discussed later in this chapter). The exponential increase in the number of images associated with each examination make the radiologist's search work more complex and time consuming. There are several key psychophysical factors that influence radiologists' ability to perform the detect-and-locate task, including perceptual ones such as target conspicuity and background clutter, as well as attentional factors that arise in the search process. The signal-to-noise ratio is a standard way to quantify the conspicuity of the target. Interventions that boost the signalto-noise ratio would be expected to improve performance. For example, with cross-sectional images, a key constraint on a reader's ability to detect and locate targets is often the physical similarity of the lesion with nearby normal structures ("anatomic noise" that is generated by distracting or indistinguishable structures on the image). Another way to expand the visual capabilities is to use techniques such as 3D rendering that can help the reader separate objects in space and differentiate normal from similarly appearing abnormal structures. Note that if a target lesion is not immediately detectable at first glance, the task becomes a visual search task, with the implied complexity (selective attention to different parts of the image). Though the searcher benefits from knowing a priori the characteristics of the target, this is not always the case in a medical image detection task (the great variability in visual aspects of pulmonary pathologies, case under study, is discussed in this chapter).

• Change Assessment

Another important task is to establish the amount of change of an abnormality over the time, during a longitudinal follow up. Typically, for Interstitial Lung Diseases, this is done in a qualitative manner with pair wise comparisons of images obtained at different time points to estimate the direction and magnitude of change. This can be a difficult perceptual task, made more difficult if the target lesion is poorly seen, or poorly defined, or exhibits multiple or complex changes over time. Acquisition differences between the images to be compared, traduced in changes in the positioning of the patient, reconstruction kernel, or image level selection (all these sources of variability are described in next sections of this introduction chapter), can make comparison harder.

• Target Characterization

A third visual task is to characterize the nature of the target lesion. This task calls on visual and cognitive resources to correctly match the appearance of the lesion to a specific pathologic entity. In the particular case of pulmonary diseases, it will be seen later in this chapter that there exists a supplementary difficulty derived from the ontological definition of each of the pathologies. Radiologists are also called on to combine information obtained from multiple types of imaging studies performed in the same patient. In this aspect, multimodality image fusion or image registration can clarify the nature of the target lesion.

• Communicating Results

Once a lesion is detected, characterized, and followed, a fourth task for the radiologist is to communicate the results of their image interpretation to the referring physician. This is done by translating visual observations into a description of the finding. This filtering and restatement process of a big amount of data implies a data compression and can introduce ambiguities or errors.

Aim of the work

Now the guidelines of radiology workflow being defined, and thereby also the challenges and limitations, the objective of this work is to investigate how the analysis and visualization can be optimized while navigating large image data sets, in order to achieve safer and better quality care for patients and a more efficient and effective work process for radiologists.

In this regard, tools such as computer-aided detection (CAD), specially based on the new paradigm of Artificial Intelligence (AI), could highlight potentially abnormal areas on a set of images, alerting the radiologist to examine them, thus reducing false-negative rates, interpretation time, and increasing true-positive rates as well as radiologists' confidence in their decisions. Accurate quantitative information in patient follow up is particularly important as input criteria for clinical decisions making. Note two additional constraints for this analysis: (1) any fully manual quantitative analysis is not possible in practice so automated tools to help radiologists processing the overwhelming (and increasing) amount of information contained in modern medical imaging studies should be developed; and (2) the measures should cover the whole organ under study (in this case the entire lung volume), instead of being restricted to few selected regions only.

To proceed with the practical background definition, the next section describes the case study involving CT modality for lung disease examination. First, the general principle of CT acquisition protocol is explained, followed by the discussion of sources of variability that have the greatest impact on image quality. The appearance of various ILD patterns on CT images is presented and discussed afterwards in section 1.3.

Clinical investigation of lung pathologies by MSCT

As a lung disease develops, it changes the structure of the lung parenchyma affecting the elastic properties of the lung, thereby the lung function. These structural changes are important to the patho-physiology of many lung disorders, such as chronic obstructive pulmonary disease (COPD), emphysema, and ILDs [START_REF] Smith | Establishing normal reference values in quantitative computed tomography of emphysema[END_REF]. MSCT can provide reliable measures of lung density, thus differentiating between normal or pathologically-induced attenuation of lung parenchyma. Visual assessment of CT scans provides useful information for identifying and sequentially evaluating the extent of altered lung structures, changes in airway dimensions, vascular distribution, expiratory air trapping, pressure of focal or diffuse disease patterns, etc., and provides a means of objectively characterizing and following these pathologic processes [START_REF] Lynch | CT-definable subtypes of chronic obstructive pulmonary disease: a statement of the Fleischner Society[END_REF].

Understanding the image modality employed for lung investigation is mandatory to apprehend not only the potential but also the limitations of the image cartography provided by this technology. It will help designing the appropriate analysis workflow in relation with the morphological patterns we need to discriminate.

We will briefly present here the principle of the CT without elaborating in detail on the theoretical aspects, but instead focusing on the acquisition parameters and their influence on the image appearance and quality.

MSCT general description

CT and MSCT acquisition principle

The CT imaging was introduced into clinical practice by Godfrey Hounsfield and Allan Mc Cormack in 1972. This technology produces high quality data of the inside of the body in cross sectional views by exploiting the principle of image reconstruction from projections [START_REF] Jain | Image Reconstruction from Projections[END_REF]. These projections are obtained as the attenuation profiles of an X-ray beam traversing the body and recorded by a detector for a multitude of directions (theoretically an infinite number of projections), generally along a cross-section plane (Figure 3). The reconstructed radiographic image thus produces a cross-sectional view of the body which depicts the local X-ray attenuation coefficients of the tissue. The linear attenuation values are expressed in Hounsfield Units (HUs) defined as the relative attenuation with respect to the water:

𝐇𝐔 = 𝟏𝟎𝟎𝟎 𝛍 𝐭 -𝛍 𝐰𝐚𝐭𝐞𝐫 𝛍 𝐰𝐚𝐭𝐞𝐫 , [ 1-1 ]
where 𝜇 𝑡 is the linear attenuation of the tissue of interest, and 𝜇 𝑤𝑎𝑡𝑒𝑟 represents the linear attenuation coefficient of the water. Image reconstruction in CT is the mathematical process that generates images from X-ray projection data acquired at many different angles around the patient. This step has a fundamental impact on image quality. The reconstruction process is very important to ensure that the CT data obtained accurately represents the true attenuation properties of the lung tissue that has been scanned. Two major categories of methods exist, analytical reconstruction and iterative reconstruction. Methods based on weighted filtered backprojection (FBP) are one type of analytical reconstruction that is currently widely used on clinical CT scanners because of their computational efficiency and numerical stability. If 𝑓 denotes the 2D attenuation coefficients across a section plane for which an infinity of projections (attenuation profiles) 𝑔(𝑠, 𝜃) along all possible directions 𝜃 are obtained (though CT scanning), the analytic reconstruction formula of 𝑓(𝑥, 𝑦) from its projections (also called the Radon transform) is given in parallel projection geometry by: 𝑓(𝑥, 𝑦) = ℬℱ -1 {|𝜉|𝐺(𝜉, 𝜃)} , [START_REF] Smith | Establishing normal reference values in quantitative computed tomography of emphysema[END_REF][START_REF] Lynch | CT-definable subtypes of chronic obstructive pulmonary disease: a statement of the Fleischner Society[END_REF] with ℬ -the backpropagation operator,

𝐵 𝑔 = ∫ 𝑔(𝑥𝑐𝑜𝑠𝜃 + 𝑦𝑠𝑖𝑛𝜃, 𝜃)𝑑𝜃 𝜋 0 , 𝑔 -the projection (Radon transform) space, 𝑔(𝑠, 𝜃) = ℛ 𝑓(𝑥, 𝑦) = ∫ ∫ 𝑓(𝑥, 𝑦)𝛿(𝑥𝑐𝑜𝑠𝜃 + 𝑦𝑠𝑖𝑛𝜃 -𝑠)𝑑𝑥𝑑𝑦 ∞ -∞ ∞ -∞ , -∞ < 𝑠 < ∞ , 0 ≤ 𝜃 < 𝜋 , 𝐺 = ℱ 𝑔 -the 1D Fourier transform of 𝑔, 𝐺(𝜉, 𝜃) = ∫ 𝑔(𝑠, 𝜃)𝑒 -𝑗2𝜋𝑠𝜉 𝑑𝑠 ∞ -∞ .
In practice, the theoretical filtering of the Fourier spectrum of projections 𝐺(𝜉, 𝜃) by multiplying with |𝜉| will amplify the noise, thus a weighting using a limited-band filter 𝑊 is applied instead:

f ̂(x, y) = ℬℱ -1 {W(ξ)G(ξ, θ)} . [ 1-3 ]
Reconstruction from attenuation profile projections W(𝛏) is usually called the reconstruction kernel (acting in the Fourier domain) and is specific to the type of image enhancement design which is manufacturer-dependent. Ideally, this kernel needs to have no edge enhancement or image smoothing. In practice, typical qualitative CT images of the lungs use a reconstruction kernel that has edge enhancement to bring out the small details of lung anatomy but at the expense of increased image noise and reduced accuracy of the CT attenuation values. Similarly, a "soft" (low-pass) reconstruction kernel that smooths the image, induces also a decrease in the accuracy of the CT values.

Utility of CT scans in lung investigation

Possibly one of the greatest advantages of CT scanner is that it provides densitometry maps of the lungs and useful information on the type and amount of tissue present (from dense bone to air) since the image scale is linearly correlated with the gravimetric density within the biological range. Furthermore, quantitative CT offers the possibility to obtain rich anatomical/structural information on lung parenchyma, thereby allowing longitudinal studies of chronic lung diseases. The analysis of the linear attenuation values assigned to lung voxels provides the basis for quantitative CT assessment of lung tissue. For example, one simple quantitative CT derived metric could be the percent of lung tissue having a density less than -950 HU on a CT scan acquired at total lung capacity (TLC), which is a marker of emphysema [START_REF] Newell | Development of quantitative CT lung protocols[END_REF]. Another example demonstrated in previous studies indicates that the overall extent of lung fibrosis on MSCT (i.e., combined extent of reticulation and honeycombing change) is a strong independent predictor of mortality in patients with idiopathic pulmonary fibrosis (IPF) [START_REF] Iwasawa | Assessment of prognosis of patients with idiopathic pulmonary fibrosis by computer-aided analysis of CT images[END_REF]. Thus, quantitative analysis of MSCT is important for management of patients with IPF.

Sources of variability in CT acquisition

Despite the great advantages of using CT modality, there are large sources of variations in the CT measurements of lung densitometry such that data from different centers and studies is often not comparable. This may also affect the visual assessment and the (computer-based) automatic detection of pathological regions extent and their characterization. This section will comment on some of the major sources of variation in CT analysis which occur commonly in studies today, to demonstrate the strong variability inherent in CT analysis and its limitations [START_REF] Coxson | Sources of variation in quantitative computed tomography of the lung[END_REF]. In implementing automated quantitative routines, it is critical to identify the invariant visual properties despite these sources of variation, in order to avoid inaccuracies in the CT number measurements.

• CT scanners manufacturers, quality control and calibration Some studies have shown that while many CT scanners from different manufacturers (Philips, Siemens, General Electric, etc.) are comparable, there can be significant differences in lung density measurements between scanners [START_REF] Coxson | Sources of variation in quantitative computed tomography of the lung[END_REF]. Importantly, even within a specific CT scanner there can be changes in the lung density measurement over time, most commonly due to incorrect calibration or changes in the X-ray tube. CT scan is a measuring device, and like all measuring devices it must be properly calibrated before use. Quality control of the CT images must use CT scanner test objects (phantoms) to provide frequent periodic checks on the CT scanner calibration to correct any drift in the CT attenuation values. This calibration process assumes that a patient will be positioned in the center of the CT gantry or isocenter when they are scanned (thus, the importance of the patient position during the test). It appears that as long as scanners are properly calibrated and the images are acquired in a consistent and reproducible manner, the data are useful in assessing disease and disease progression.

• CT reconstruction kernel

As a part of an optimal CT protocol for a particular study, it is important to define an optimal reconstruction kernel. In Figure 4 it is shown that different reconstruction kernels from two different CT scanners result in quite different CT values frequency distribution. Here it can clearly be observed that the extent of emphysema, for example, measured as the percent of voxels of CT values below -950 HU would really depend on the kernel chosen to reconstruct the images [START_REF] Coxson | Sources of variation in quantitative computed tomography of the lung[END_REF].

Figure 4. CT value frequency distribution curves for five different reconstruction kernels from two CT manufacturers from a single subject. The vertical arrow points at the -950 HU threshold cut-off for low attenuation analysis of emphysema. It is shown that the reconstruction algorithm can make significant differences in the "extent of emphysema" measured [START_REF] Coxson | Sources of variation in quantitative computed tomography of the lung[END_REF].

• CT radiation dose

The signal to noise ratio of the CT image is proportional to the radiation dose the patient receives. There is clearly a need to keep the patient radiation dose to a minimum but also ensure adequate signal to noise ratio for the study so that meaningful CT attenuation values are obtained. Thus, the CT radiation dose level used for the evaluation of COPD and ILD is driven by the balance between radiation exposure and image quality. Adequate visual characterization can be achieved with reduced dose CT acquisition techniques, as used for lung cancer screening. However, for pathology grading, excessive image noise due to a reduced CT dose can "simulate" emphysema, particularly at quantitative CT, and may impair segmentation of the airways and quantitative evaluation of airway wall thickness [START_REF] Lynch | CT-definable subtypes of chronic obstructive pulmonary disease: a statement of the Fleischner Society[END_REF] (see Figure 5 for a visual example of the influence of radiation dose over the image fidelity).

For a given CT scanner hardware and software configuration, the X-ray radiation dose is determined by the product of the tube current and total exposure time (mAs) and also by the peak (kilo)voltage, kVp, applied to the X-ray tube. It is important that all subjects/patients are scanned at the same kVp since changing the kVp will change the linear attenuation properties of the matter and so change the measured CT linear attenuation values assigned to the lung voxels independent of any other factor. The larger the patient's body mass index (BMI, defined as a person's weight in kilograms divided by the square of the body height in 𝑚 2 ), the greater the mAs value needs to be selected, in order to ensure that enough X-ray photons hit the X-ray detector (the larger mass of the patient results in greater attenuation of the X-ray photons produced). A andC) can be seen as linear "streaks" or a more "mottled" image while the clinical dose CT scans show a "smoother" image. Panel E shows the histogram of x-ray attenuation values from one of the subjects to illustrate that there is a difference in the extent of emphysema measured using either the threshold or percentile analysis (curve 𝒊 clinical dose CT scan, curve 𝒋 low dose CT scan). Example taken from [START_REF] Coxson | Sources of variation in quantitative computed tomography of the lung[END_REF].

𝑖 𝑗

• Biological or physiological differences Subject-specific factors known to influence lung attenuation include age, gender, body size, socio-economic status, depth of inspiration, smoking, and gas-trapping [START_REF] Smith | Establishing normal reference values in quantitative computed tomography of emphysema[END_REF]. However, the largest variation in lung attenuation is driven by the ability of the patient to inspire or expire to a given lung volume, maintain that lung volume for the duration of the CT scan, 12 seconds or less, and not physically move while the CT scanner is acquiring image data (to avoid motion artifacts). A deeper level of inspiration during the scan results in a lower mean lung density and more low attenuation regions because there is more gas volume compared to tissue volume. Differences in lung volume will result in differences in lung density, and even relatively small changes in lung volume may affect measurements of a given pathology and its extent. Lung density can vary by as much as 80-100 HU from full-inspiration to end-expiration. Consistent measurements are most likely to be achieved at maximum inspiration because the variation in CT lung density is the lowest at full inspiration [START_REF] Lynch | Quantitative computed tomography of diffuse lung disease[END_REF].

To conclude, ideally, there would be one CT scanner and software system using exactly the same parameters specified in the CT scanning protocol to reconstruct images of the lungs on all subjects in a clinical trial to study given pathologies. This brief review showed some of the major sources of variations that can modify the visual perception of CT images in real clinical practice and which should be considered in the image analysis algorithm development for longitudinal studies.

Emphysema and interstitial lung diseases: anatomical description

The image appearance and the salient parenchymal radiological features on CT for each of the lung pathologies we are interested in through this work -emphysema, fibrosis, honeycombing, ground glass, and airspace enlargement with fibrosis (AEF) are discussed in the following. As we mentioned in the introduction, the purpose is to describe and define the phenotypic abnormalities that can be identified on visual and quantitative evaluation of CT images in subjects with COPD or ILD, with the goal of contributing to a personalized approach to the treatment. Although these abnormalities often overlap, identification and quantification of the predominant morphologic findings and their grouping into defined subtypes improve diagnostic accuracy, help optimizing treatment, and provide a framework for data comparison in clinical trials [START_REF] Lynch | CT-definable subtypes of chronic obstructive pulmonary disease: a statement of the Fleischner Society[END_REF].

Here are some well-known categories of lung lesions seen in CT imaging. Notice that this taxonomy is neither complete nor widely accepted at present as there are no standardized classes defined in the literature. However, these CT signs frequently appear in pathologic lungs and play important roles in the diagnosis of the disease.

Emphysema

Pulmonary emphysema is a type of COPD defined morphologically by airspace enlargement with alveolar wall destruction and gradual loss of lung tissue. This kind of obstructive airway disease is characterized by limitation of airflow [START_REF] Smith | Establishing normal reference values in quantitative computed tomography of emphysema[END_REF]. Emphysema lesions, or bullae, are visible in CT images as areas of abnormally low attenuation values close to that of air. Emphysema can be divided into three subtypes, or patterns describing the sizes and the spatial distribution of the lesions: centrilobular emphysema (CLE), defined as multiple small low-attenuation areas; paraseptal emphysema (PSE), defined as multiple low-attenuation areas in a single layer along the pleura often surrounded by interlobular septa that is visible as thin white walls; and panlobular emphysema (PLE), defined as a low-attenuation lung with fewer and smaller pulmonary vessels [START_REF] Sorensen | Quantitative analysis of pulmonary emphysema using local binary patterns[END_REF]. Emphysema lung destruction results in replacement of the normal tissue (which has a typical attenuation of about -850 HU on inspiratory CT) with air-containing spaces, of CT attenuation close to -1000 HU [START_REF] Lynch | Quantitative computed tomography of diffuse lung disease[END_REF].

• Centrilobular Emphysema (CLE)

The primary lesion is produced by dilation and destruction of respiratory bronchioles within a single acinus. The centrilobular lesions are formed by coalescence of several primary lesions. Subsequently, the destruction spreads to the entire lung lobule and fuses many destroyed lobules together to produce a pattern of coalescent destruction that sometimes disintegrates to form large bullous lesions. At CT, CLE is characterized by small well-defined or poorly defined areas of low attenuation surrounded by normal lung. Centrilobular pulmonary arteries or arterioles, which are often seen traversing the hypoattenuated areas, mark the center of each lobule. This is the most common type of smoking related emphysema and is usually upper lung predominant [START_REF] Lynch | CT-definable subtypes of chronic obstructive pulmonary disease: a statement of the Fleischner Society[END_REF]. As CLE becomes more severe, the areas of low attenuation become confluent (Figure 6) and the centrilobular distribution becomes less apparent. In most cases, the areas of low attenuation have no visible walls, otherwise, very thin walls may be seen. Advanced destructive emphysema is manifested as a generalized decrease of attenuation of the lung without focal hypoattenuation and represents an advanced stage of CLE (this pattern may be indistinguishable at CT from the panlobular pattern described below). Confluent emphysema may be differentiated from advanced destructive emphysema by the presence of a preserved rim of normal lung attenuation intervening between areas of lung destruction, and by the absence of lobular hyperexpansion, architectural distortion, or splaying or decreased caliber of vessels. • Panlobular Emphysema (PLE)

PLE specifically refers to diffuse emphysematous destruction across the lobule. In general, the extent and severity of alveolar destruction in PLE is milder than that in CLE, but it affects all of the acini within a lung lobule more or less equally (see the pathology example in Figure 7).

• Paraseptal Emphysema (PSE)

Paraseptal emphysema describes emphysematous lesions caused by selective destruction of the distal acinus, and subsequent reports have used it to describe lesions located near the pleural surface close to the chest wall and in the interlobar fissures. In some cases, multiple destroyed acini coalesce to form striking lesions just under the pleural surface on CT scans. CT shows subpleural areas of low attenuation with a well-defined wall (Figure 8). Rows of PSE may mimic honeycombing, but the size of the cysts is larger than that of honeycombing cysts and architectural distortion and other signs of fibrosis are not present. PSE is commonly associated with marked thickening of the walls of proximal bronchi and bronchioles, suggesting a significant airway inflammatory component. PSE occurs across the entire spectrum of minimal involvement to severe parenchymal obstruction and can be progressive. Bullae (avascular low-attenuation areas above 1 cm in diameter, with a thin but perceptible wall) are found in all types of emphysema but are most commonly associated with PSE [START_REF] Lynch | CT-definable subtypes of chronic obstructive pulmonary disease: a statement of the Fleischner Society[END_REF]. Bullae are often located in the upper lobes in both CLE and PSE but are more evenly distributed in the lungs of patients with advanced destructive emphysema.

As a final remark, we show the complete evolution of emphysema as the destruction of lung parenchyma progresses from normal lung tissue to mild, moderate, and severe disease, with complete deletion of the lung architecture (see Figure 9). As all the stages of emphysema must be recognized by the CAD system, it is obvious how complex the classification problem is. It differs from a standard classification problem because of the inherent natural ordering of both the classes, and also the patterns within a single class. It is indeed very difficult to assign a continuous score of disease severity for each pattern by visual evaluation. Clinician only assigns class labels to patterns [START_REF] Kurugol | Ranking and classification of monotonic emphysema patterns with a multi-class hierarchical approach[END_REF]. Note that this idea will also be applied to ILDs described in next section.

Interstitial lung diseases (ILDs)

Diffuse Lung Diseases is a group of lung diseases that primarily affect the lung parenchyma. This group is referred to by the generic term diffuse parenchymal lung diseases (DPLDs), or interstitial lung diseases (ILDs). They are characterized by the inflammation of the lung interstitium (tissue that surrounds and separates the tiny air sacs, alveoli, of the lung parenchyma). Four types of ILDs, namely idiopathic pulmonary fibrosis, honeycombing, ground-glass opacity, and airspace enlargement with fibrosis (AEF) are resulted by such conditions and are briefly described in the following.

Fibrosis

The most relevant CT signs of fibrotic interstitial lung disease are reticular pattern, traction bronchiectasis, ground glass opacity, and honeycombing. Despite standardized definitions, agreement between observers on the presence and extent of these basic patterns is variable. Reasons for these differing opinions include overlap in the CT appearances of honeycombing, traction bronchiectasis, and paraseptal emphysema. Even among experienced thoracic radiologists, agreement on the presence or absence of honeycombing is moderate. An important secondary feature of interstitial fibrosis is distortion and dilatation of bronchi in the fibrotic lung; agreement between observers on the presence or absence of so-called traction bronchiectasis is better than agreement with respect to honeycombing [START_REF] Hansell | CT staging and monitoring of fibrotic interstitial lung diseases in clinical practice and treatment trials: a position paper from the Fleischner Society[END_REF]. Bronchiectasis is defined at CT as a dilated bronchial lumen relative to the adjacent pulmonary artery, lack of bronchial tapering, or identification of bronchi within 1 cm of the pleural surface, and it is a known indicator of the presence of fibrosis surrounding this region [START_REF] Lynch | CT-definable subtypes of chronic obstructive pulmonary disease: a statement of the Fleischner Society[END_REF].

Honeycombing

Honeycombing is one of the characteristic appearances of pulmonary fibrosis in the endstage and is extensive pulmonary fibrosis with alveolar destruction resulting in a loss of alveolar walls. As the pulmonary fibrosis progresses, honeycombing becomes more prominent. On CT scans, it is defined as clustered cystic airspaces, typically of comparable diameters of the order of 3-10 mm, which are usually subpleural and have well-defined walls. Although the definition of honeycombing seems to be strict, recognition of honeycombing on CT is various among each observer. Air-space consolidation in the presence of pulmonary (paraseptal) emphysema can mimic this appearance.

Ground glass

Ground glass is a nonspecific term referring to the appearance of hazy lung opacity and is relatively easy to recognize because of the presence of the thickened alveolar interstitium, or the presence of the fluid partially filling the alveolar spaces. In histology, ground-glass opacity has been seen in patients with the findings of mild or early interstitial inflammation or infiltration and should be considered as consistent with active inflammation only when there are no superimposed findings of fibrosis such as reticulation, architectural distortion, or traction bronchiectasis. In some diseases, such as idiopathic pulmonary fibrosis (IPF) and sarcoidosis, the appearance of ground-glass opacity correlates with disease activity. These areas of ground-glass opacity have been shown to correspond to regions of active alveolitis and precede irreversible changes such as fibrosis and honeycombing [START_REF] Engeler | Ground-glass opacity of the lung parenchyma: a guide to analysis with high-resolution CT[END_REF]. Sometimes it is difficult to decide which lung density is abnormal, that is, whether lung opacity has increased or decreased. Normal lung tissue may have a ground-glass appearance relative to hypoventilated and hypoperfused emphysematous lung.

Air space enlargement with fibrosis (AEF)

AEF is pathologically characterized by multiple thin-walled cystic lesions (i.e. abnormal airspaces), histologically intensive hyalinized fibrosis and a bronchiolo-centric location. The histologic characteristics of AEF are comparable with usual interstitial pneumonia (UIP), idiopathic pulmonary fibrosis (IPF) and centrilobular emphysema (CLE). However, studies have shown significant differences in many parameters studied between them. AEF appears to represent emphysema with fibrotic walls. It frequently accompanies CLE. Fibrosis in AEF differed from that in UIP and IPF, more hyalinized in character and with little granulation tissue [START_REF] Yamada | Airspace enlargement with fibrosis shows characteristic histology and immunohistology different from usual interstitial pneumonia, nonspecific interstitial pneumonia and centrilobular emphysema[END_REF]. For AEF samples CT images, crosssections of lung show multiple thin-walled cystic lesions of various sizes. A panoramic image of a lung section from an AEF patient also shows multiple cystic lesions with thin walls in the subpleura and deeper parenchyma with tendency of respiratory bronchiolar distribution. The surface of the pleura is quite smooth without thickening or constriction. As mentioned before, the appearance seems to be emphysema, but cystic lesions have fibrous wall and subpleural distribution is different from centrilobular emphysema. Also, the cystic lesion walls of AEF are significantly thicker than that of CLE, and the latter had less inflammation than AEF and no granulation tissue.

To summarize ILDs anatomical description, we present several examples of these categories in Figure 10.

Challenges in Emphysema-ILDs classification

Several difficulties are implied when designing robust CAD algorithms and were progressively enumerated in previous sections showing the source of variations in CT imaging and lung texture analysis. The fundamental clinical problems of how to consistently detect, characterize, and differentiate the various ILDs remain diagnostic challenges. We summarize them in this section:

1. Image artifacts/noise: gravitational gradient from top to bottom of the image (anteroposterior), motion artifacts (see Figure 11 for examples); Standardized acquisition methods are required ideally to enable high precision and scientific evaluation and validation of quantitative measures. However, we have shown that, in practice, CT images are created using many different image acquisition and reconstruction parameters and these parameters vary from vendor to vendor. Due to such variations, it is important to determine how changes in these parameters affect radiomics features (i.e. robustness of features with variations in different parameters). Raw data normalization seems an appropriate step to be included in the workflow to come up with reproducible features. Also, this kind of invariance (namely pixel intensity value translations) should be considered to be added to the overall classification system.

Overlapping CT visual features of different lung classes (section 1.3).

It is important to mention the inherent problem relying on ILD classification task: the great overlap between the classes often makes the task an ontological problem which leads to errors in the classification, i.e., mixing different types of lung patterns, overestimating the pathological areas. Figure 12 shows the resemblance between some of the classes to be detected: normal, emphysema, fibrosis, ground glass, honeycombing, AEF.

Figure 13 synthesizes the overlap between classes (very fuzzy frontiers, as it was already seen in the anatomical descriptions of the classes to be a challenging aspect of this task). In the results section later in this work, we will return to this overlap scheme to show these correlations among real classed samples. Given these variables that have an important impact on visual inference, we consider that the orientation for a proper analysis must include not only image processing information (whether it is based on texture analysis, or intensity signatures, etc.), but also a complementary ontological or structural knowledge.

Detection and classification of lung pathologies: state of the art

ILDs require regular quantitative patient follow-up with CT imaging which implies the development of automated quantitative and reproducible CAD tools for lung texture classification. Manual classification and evaluation of pathology extent is tedious and not reproducible, due to significant interobserver and intraobserver variability among other limitations. This inconsistency in clinical assessment of ILDs is a good reason for searching automation on qualitative and quantitative image analysis.

Classification is the task of assigning objects to one of several predefined categories, i.e., mapping an input attribute set 𝑥 into its class label 𝑦. The input data for a classification task is a collection of samples. Each sample, also known as an instance or example, is characterized by a tuple (𝑥 , 𝑦), where 𝑥 is the attribute set (in our case the raw image data with pixel intensity values) and 𝑦 the class label (also known as category or target). Whether the raw input data set can contain continuous features, the class label, on the other hand, must be a discrete attribute. This is a key characteristic that distinguishes classification from regression, a predictive modeling task in which 𝑦 is a continuous attribute. Classification can also be seen as the task of learning a target function 𝑓 that maps each attribute set 𝑥 to one of the predefined class labels 𝑦. This mapping function is also known as a classification model. A classification model is useful for predictive purposes (predict the class label of unknown records). As shown in Figure 14, a classification model can be treated as a black box that automatically assigns a class label when presented with the attribute set of an unknown record.

Traditional pattern recognition systems are composed of three main modules: a space transformation module, a fixed feature extractor and a general-purpose classifier. This means that input images are represented by features that are the actual input to a classifier which reproduces category labels. The complexity of the problem comes from the input data that has high dimensionality. Note that typical size of a CT image is 512x512 pixels meaning a feature vector of 262144 coordinates that can have values of 12-bit depth. It is indeed a high dimensional problem [START_REF] Hastie | The elements of statistical learning : Data mining, inference, and prediction[END_REF], where the number of features 𝑝 is much larger than the number of observations (𝑝 ≫ 𝑁). So, there is a need to find a better way to represent images than pixel intensity values. The space transformations and the feature extractor convert the input raw image into low-dimensional vectors or short strings of symbols that: 1) can be easily matched or compared and 2) are relatively invariant with respect to transformations or distortions of the input that however do not change its nature. The feature extractor contains most of the prior knowledge and is rather specific to the task. It is also the focus of most of the design effort, because it is often entirely handcrafted.

The classifier, on the other hand, often general purpose and, normally (in the case of machine learning approaches) trainable, predicts the output when given an input vector (that, as we saw, cannot be the raw data but instead a representation of it). One of the main problems with this approach is that the recognition accuracy is largely determined by the ability of the designer to come up with an appropriate set of features. This turns out to be a daunting task which, unfortunately, must be redone for each new problem.

In this chapter we focus on textural-based classification and we describe different approaches that can be found in the literature related to space transformation, feature extraction and classification within the framework of lung texture analysis. As discussed in chapter 1, lung disorders are characterized by specific abnormal findings, mostly texture-like in appearance, and it is the occurrence of several such findings that can point toward a specific diagnosis. For this reason, computer analysis of DPLD is alternatively viewed as a texture analysis problem and generally tackled by proposing textural features [START_REF] Sluimer | Computer analysis of computed tomography scans of the lung: a survey[END_REF]. A texture in digital images is simply an arrangement of local patterns. Such patterns can be described by a given spatial organization of gray levels (e.g., random, periodic). Textures are by definition stationary, so they are agnostic to spatial information. We are focusing here on this aspect, while the knowledge enrichment by including ontological information related to the structural distortion induced in the lung tissue will be addressed later.

Space transformation

The first block building a classification system is the space transformation. It represents a data preprocessing step that may include any kind of mathematical transformation (or projection) of the original data into a new space that may or not signify a reduction of dimensionality. The idea is to extract some second layer information from the data, e.g. interesting frequencies of a signal using Wavelet Transform, or simply prepare it to further analysis. The transformation may improve the process considering only the most relevant representation of the original data, or by filtering noisy artifacts. Other reasons for data transformation and dimensionality reduction are: removing redundancy in data, compression of data sets, obtaining transformed and reduced cases containing only relevant variables that aid to design classifiers with better generalization capabilities, projecting high-dimensional data (preserving intrinsic data topology) onto low-dimensional space to visually show clusters and other relationships within data.

Filter banks and wavelets transform

Filtering techniques have several useful properties for texture classification. First, they provide continuous responses to transient patterns in images, as they are able to characterize the similarity between series of contiguous pixels, which are often carrying important information in medical images. Filtering allows us to seek for specific features in the images (i.e., edge or ridge detection) by modeling the shape of the filters either in the spatial or in the frequency domain.

In particular, Gabor filter banks are used for texture analysis. By observing different frequencies and orientations of images, Gabor filters have been shown to mimic how actual human visual system works. In the spatial domain, a two-dimensional Gabor filter is a Gaussian kernel function modulated by a complex sinusoidal plane wave as defined by [START_REF] Movellan | Tutorial on Gabor Filters[END_REF]:

𝑮(𝒙, 𝒚) = 𝒔(𝒙, 𝒚)𝒘 𝒓 (𝒙, 𝒚) , [ 2-1 ]
where 𝑠(𝑥, 𝑦) is a complex sinusoid, known as the carrier, and 𝑤 𝑟 (𝑥, 𝑦) is a 2D Gaussian-shaped function, known as the envelope. They are defined as follows:

𝑠(𝑥, 𝑦) = exp (𝑗(2𝜋(𝑢 0 𝑥 + 𝑣 0 𝑦) + 𝑃)) , [ 2-2 ]
where (𝑢 0 , 𝑣 0 ) defines the spatial frequency that can be expressed in polar coordinates as magnitude 𝐹 0 and direction 𝜔 0 :

𝑢 0 = 𝐹 0 cos 𝜔 0 , 𝑣 0 = 𝐹 0 sin 𝜔 0 , and 𝑃 the phase of the sinusoid respectively.

𝒘 𝒓 (𝒙, 𝒚) = 𝑲 𝐞𝐱𝐩 (-𝝅(𝒂 𝟐 (𝒙 -𝒙 𝟎 ) 𝒓 𝟐 + 𝒃 𝟐 (𝒚 -𝒚 𝟎 ) 𝒓 𝟐 )) , [START_REF] Lynch | CT-definable subtypes of chronic obstructive pulmonary disease: a statement of the Fleischner Society[END_REF][START_REF] Newell | Development of quantitative CT lung protocols[END_REF] where (𝑥 0 , 𝑦 0 ) denotes the spatial location of the function peak (Gaussian envelope), K scales the magnitude of the Gaussian envelope, a and b are scaling parameters of the Gaussian, and the r subscript stands for a rotation operation such that:

(𝑥 -𝑥 0 ) 𝑟 = (𝑥 -𝑥 0 ) cos 𝜃 + (𝑦 -𝑦 0 ) sin 𝜃 , [START_REF] Lynch | CT-definable subtypes of chronic obstructive pulmonary disease: a statement of the Fleischner Society[END_REF][START_REF] Newell | Development of quantitative CT lung protocols[END_REF][START_REF] Iwasawa | Assessment of prognosis of patients with idiopathic pulmonary fibrosis by computer-aided analysis of CT images[END_REF] (𝑦 -𝑦 0 ) 𝑟 = -(𝑥 -𝑥 0 ) sin 𝜃 + (𝑦 -𝑦 0 ) cos 𝜃 .

CAD systems using HRCT images have been proposed by measuring textures derived from Gabor filtered images [START_REF] Ebrahimian | Gabor texture measure in interstitial lung disease discrimination using high resolution computed tomography thorax images[END_REF]. A set of Gabor filters is shown in the next figure. On the other hand, Wavelet Transforms (WTs) have the desirable property of covering the frequency domain at several scales (Figure 16). Translation-invariant wavelets offer a feature set (Wavelet coefficients) able to fit most of the texture functions, being perfectly complementary to the measures of density using other image features like grey level histograms (GLH) [START_REF] Depeursinge | Near-affine-invariant texture learning for lung tissue analysis using isotropic wavelet frames[END_REF].

Figure 16. Representation of different transformation to a 1D signal [START_REF] Misiti | Wavelet toolbox[END_REF]. Note the benefits of using WT by the multiscale analysis obtained.

Continuous Wavelet transforms is defined as:

𝑊 𝜑 (𝑠, 𝜏) = ∫ 𝑓(𝑥)𝜓 𝑠,𝜏 (𝑥)𝑑𝑥 ∞ -∞ , [ 2-5 ] Directions ( 𝜃 = 𝜔 0 = 0, 𝜋 8 , 2𝜋 8 , … , 7𝜋 8
)

Frequencies (𝐹 0 )
where 𝑊 𝜑 is the result of transforming an input continuous function f(x), 𝜓 𝑠,𝜏 is the wavelet function, 𝑠 is the scale, 𝜏 the translation. The same idea can be extended to 2D (discrete) functions (as images) transformation. The wavelet coefficients measure how closely the constituent wavelets of different scales and positions correlates with each section of the signal [START_REF] Misiti | Wavelet toolbox[END_REF]. Note that transforming the original signal to a set of coefficients can be also seen as a reduction of dimensionality which makes WT useful also for compression tasks. WT can be efficiently implemented as filter banks as shown in the next figure.

Structural approaches

Morphological image processing describes a range of image processing techniques that deal with the shape (or morphology) of features in an image. It is useful for image filtering (shape simplification, enhancing object structure, etc.); image segmentation (watersheds); image measurements (area, perimeter, granulometry); pattern recognition; and texture analysis [67] [68]. Mathematical Morphology includes the theory for the analysis of spatial structures in terms of shapes and topology, and is based on integral geometry, lattice algebra and structuring elements. Using simple operations (erosion, opening, closing, dilation) it is possible to build up several morphological operators: hit-or-miss transform, boundary extraction, region filling, extraction of connected components, thinning/thickening, etc. [START_REF] Shyu | ASSERT: A physician-in-the-loop content-based retrieval system for HRCT image databases[END_REF] [START_REF] Kim | Development of an automatic classification system for differentiation of obstructive lung disease using HRCT[END_REF].

𝐻 𝑥 𝐿 𝑦 𝐻 𝑥 𝐻 𝑦 𝐿 𝑥 𝐻 𝑦 𝑗 = 0 𝑗 = 1
Morphological filters can also be used for image space transformation. A morphological filter is an operator : R n → R satisfying two properties -increasingness and idempotence. Morphological filters can be involved in space transformation operations with respect to either content simplification (noise filtering, contrast enhancement) or content segmentation. We recall here several useful filters, starting with the basic ones -the opening and closing -then increasing the complexity -alternating (sequential) filters, connected filters (grayscale reconstruction-based) to end up with a segmentation-oriented operator (though not a filter) -the watershed (preserving only the idempotence property).

• Opening, closing and alternating filters

Let f : R n → R and B  R n a structuring element (SE). The opening of f by B is defined as:

f o B = (f  B)  B , [ 2-6 ]
where  and  denote the erosion and the dilation operator, respectively [START_REF] Serra | Image Analysis and Mathematical Morphology[END_REF]. The closing of f by B is defined as:

f • B = (f  B)  B.
[ 2-7 ]

Figure 19 illustrate the two filters in a 1D example.

AFoc(f, k) = k(f) = (f o kB) • kB , [ 2-8 ] AFco(f, k) = k(f) = (f • kB) o kB
The alternating sequential filters open-close and close-open are defined by (using the previous notations):

ASFoc(f, k) =  k  k-1  k-2 ….  1 (f) [ 2-9 ] ASFco(f, k) =  k  k-1  k-2 ….  1 (f)
Figure 20 shows an example of results obtained by applying the previous filters to an axial lung CT image using a disk structuring element of radius 3 pixels.

(a) original axial CT, f (b) f o B (c) f • B (d) AFoc(f, 3) = 3(f) (e) AFco(f, 3) = 3(f) (f) ASFoc(f, 3) (g) AFco(f, 3)

• Grayscale reconstruction by dilation/erosion

Let f, g: R n → R and g f  . The grayscale reconstruction by dilation of f by g is defined as [START_REF] Vincent | Morphological grayscale reconstruction in image analysis: applications and efficient algorithms[END_REF]:

) ( ) ( ) ( g g f f  =    , [ 2-10 ] with )) ( ( ) ( ) ( ) ( ) 1 ( ) 1 ( ) ( ) 1 ( g g f H g g n f f n f f - =   =     , [ 2-11 ]
where H is a unit structuring element and

) ( ) 1 ( g f 
is the geodesic elementary dilation of g in f (∧ denotes the infimum operator). Figure 21(a) illustrates the reconstruction by dilation of f by g on a 1D example. Note that the reconstruction by dilation of f by g reconstructs only the "hills" of f marked by the regional maxima of g and imposes plateaus at the level of g maxima. In the particular case of binary data,

) (g f  
becomes the binary set reconstruction of f from g.

Let f, g: R n → R and g f  . The grayscale reconstruction by erosion of f by g is defined as [START_REF] Vincent | Morphological grayscale reconstruction in image analysis: applications and efficient algorithms[END_REF]:

) ( ) ( ) ( g g f f  =    , [ 2-12 ] with )) ( ( ) ( ) ( ) ( ) 1 ( ) 1 ( ) ( ) 1 ( g g f H g g n f f n f f - =   =     , [ 2-13 ] where ) ( ) 1 ( g f 
is the geodesic elementary erosion of g in f (∨ denotes the supremum operator). Figure 21(b) illustrates the reconstruction by erosion of f by g on a 1D example. Note that the reconstruction by erosion of f by g fills in the "basins" of the relief f by creating plateaus at a height not overriding the level of the local minima of g. Thus, the function g performs as a "sink" with respect to the "basins" filling of f.

• Grayscale reconstruction by opening/closing

Let f : R n → R and B  R n a structuring element. Using the previous notations, the grayscale reconstruction by opening of f by B is defined as:

𝑅 𝑜 (𝑓, 𝐵) = 𝜌 𝑓 𝛿 (𝑓 𝐵) . [ 2-14 ]
Similarly, the grayscale reconstruction by closing of f by B is given by 𝑅 𝑐 (𝑓, 𝐵) = 𝜌 𝑓 𝜀 (𝑓 𝐵) . [START_REF] Lynch | CT-definable subtypes of chronic obstructive pulmonary disease: a statement of the Fleischner Society[END_REF][START_REF] Newell | Development of quantitative CT lung protocols[END_REF][START_REF] Iwasawa | Assessment of prognosis of patients with idiopathic pulmonary fibrosis by computer-aided analysis of CT images[END_REF][START_REF] Coxson | Sources of variation in quantitative computed tomography of the lung[END_REF][START_REF] Lynch | Quantitative computed tomography of diffuse lung disease[END_REF][START_REF] Sorensen | Quantitative analysis of pulmonary emphysema using local binary patterns[END_REF][START_REF] Kurugol | Ranking and classification of monotonic emphysema patterns with a multi-class hierarchical approach[END_REF][START_REF] Hansell | CT staging and monitoring of fibrotic interstitial lung diseases in clinical practice and treatment trials: a position paper from the Fleischner Society[END_REF][START_REF] Engeler | Ground-glass opacity of the lung parenchyma: a guide to analysis with high-resolution CT[END_REF][START_REF] Yamada | Airspace enlargement with fibrosis shows characteristic histology and immunohistology different from usual interstitial pneumonia, nonspecific interstitial pneumonia and centrilobular emphysema[END_REF][START_REF] Sluimer | Computer analysis of computed tomography scans of the lung: a survey[END_REF][START_REF] Chien | Detection and classification of interstitial lung diseases and emphysema using a joint morphological-fuzzy approach[END_REF][START_REF] Kim | Development of an automatic classification system for differentiation of obstructive lung disease using HRCT[END_REF][START_REF] Xu | MDCT based 3D texture classification of emphysema and early smoking related lung pathologies[END_REF] Figure 22 illustrates the space transformation operated by the last two connected filters applied on the image in Figure 20 (a) using the same SE disk of radius 3 pixels. Note that the original structures are better preserved here through filtering. The previous space transformations using morphological filters (or combination of them) are useful to strengthen selected characteristics of the textural patterns that might be beneficial for feature extraction and classification. Alternatively, they can be used in (pre-) segmentation tasks allowing to target specific regions in the image. In this respect, they can be combined with the watershed operator (defining the ridges separating catchment basins of regional minima in the image [START_REF] Serra | Image Analysis and Mathematical Morphology[END_REF]. In [START_REF] Chien | Detection and classification of interstitial lung diseases and emphysema using a joint morphological-fuzzy approach[END_REF] a three-dimensional approach based on mathematical morphology and fuzzy logic for quantifying and classifying interstitial lung diseases (ILDs) and emphysema was proposed. The methodology is composed of several stages, mainly: (1) an image multi-resolution decomposition scheme based on the 3D 𝑅 𝑐 (𝑓, 𝐵) filter eq. [START_REF] Lynch | CT-definable subtypes of chronic obstructive pulmonary disease: a statement of the Fleischner Society[END_REF][START_REF] Newell | Development of quantitative CT lung protocols[END_REF][START_REF] Iwasawa | Assessment of prognosis of patients with idiopathic pulmonary fibrosis by computer-aided analysis of CT images[END_REF][START_REF] Coxson | Sources of variation in quantitative computed tomography of the lung[END_REF][START_REF] Lynch | Quantitative computed tomography of diffuse lung disease[END_REF][START_REF] Sorensen | Quantitative analysis of pulmonary emphysema using local binary patterns[END_REF][START_REF] Kurugol | Ranking and classification of monotonic emphysema patterns with a multi-class hierarchical approach[END_REF][START_REF] Hansell | CT staging and monitoring of fibrotic interstitial lung diseases in clinical practice and treatment trials: a position paper from the Fleischner Society[END_REF][START_REF] Engeler | Ground-glass opacity of the lung parenchyma: a guide to analysis with high-resolution CT[END_REF][START_REF] Yamada | Airspace enlargement with fibrosis shows characteristic histology and immunohistology different from usual interstitial pneumonia, nonspecific interstitial pneumonia and centrilobular emphysema[END_REF][START_REF] Sluimer | Computer analysis of computed tomography scans of the lung: a survey[END_REF][START_REF] Chien | Detection and classification of interstitial lung diseases and emphysema using a joint morphological-fuzzy approach[END_REF][START_REF] Kim | Development of an automatic classification system for differentiation of obstructive lung disease using HRCT[END_REF][START_REF] Xu | MDCT based 3D texture classification of emphysema and early smoking related lung pathologies[END_REF] used to detect and analyze the different density patterns of the lung texture; (2) for each pattern in the multi-resolution decomposition, six features are computed, for which fuzzy membership functions define a probability of association with a pathology class; (3) for each pathology class, the probabilities are combined up according to the weight assigned to each membership function and two threshold values are used to decide the final class of the pattern. This approach will be explained in detail later as we return to this subject in one of our proposed joint analysis scheme.

Dimensionality reduction techniques

Dimensionality reduction methods convert the high-dimensional data set 𝑋 = (𝑥 1 , 𝑥 2 , … , 𝑥 𝑁 ) into a lower-dimensional dataset 𝑌 = (𝑦 1 , … , 𝑦 𝑀 ). In this work, we are interested in two-or threedimensional dataset representation 𝑌 that can be displayed in a scatterplot. In the following, we refer to the low-dimensional data representation 𝑌 as a map, and to the low-dimensional representations 𝑦 𝑖 of individual datapoints as map points. The aim of dimensionality reduction is to preserve as much of the significant structure of the high-dimensional data as possible in the lowdimensional map [START_REF] Maaten | Visualizing data using t-SNE[END_REF]. In this section we focus on the two most common techniques: principal component analysis, and t-SNE.

Principal component analysis (PCA)

There are a great number of unsupervised methods in the literature (principal components, multidimensional scaling, self-organizing maps, cluster analysis, and principal curves) that attempt to identify low-dimensional manifolds within the original space that represent high data dimensionality. This provides information about the associations among the variables and whether or not they can be considered as functions of a smaller set of underlying variables. In this section we are covering one of them, PCA, to explain the basic behind these clustering methods.

Principal components determine an optimal linear transformation (projection of the data) of a real-valued 𝑝-dimensional data vector 𝑥 into another 𝑞-dimensional space, where 𝑞 ≤ 𝑝. This means a dimensionality reduction, in terms of the number of variables of the original unsupervised data set. The principal components of the 𝑥 set provide a sequence of best linear approximations to that data. PCA is performed by calculating correlations among elements of the original data vectors and finding the representation retaining the maximum nonredundant and uncorrelated intrinsic information.

Basically, the original data set 𝑥 is used to compute its covariance matrix, its eigenvalues and the corresponding eigenvectors arranged in descending order. The arrangement of subsequent columns of a transformation matrix as the normalized eigenvectors, corresponding to the subsequent largest eigenvalues of the covariance matrix, will result in an optimal linear transformation matrix (discussed below). The elements of transformed vector 𝑦 will be uncorrelated and arranged in descending order according to decreasing information content. This allows for a straightforward reduction of dimensionality, and thus data compression, by discarding the lowest information content elements.

Consider 𝑥 ′ the rank-𝑞 linear model for representing a set of 𝑁 observations 𝑥 𝑖 (𝑥 1 , 𝑥 2 , … , 𝑥 𝑁 ) each belonging to a rank-𝑝:

𝑥 𝑖 ′ = 𝑓(𝜆 𝑖 ) = 𝑢 + 𝑉 𝑞 . 𝜆 𝑖 , [ 2-16 ]
where 𝑢 is a location vector in 𝑅 𝑝 , 𝑉 𝑞 is a 𝑝 × 𝑞 matrix with 𝑞 orthogonal unit vectors as columns, and 𝜆 𝑖 is a 𝑞 vector of parameters. This is the parametric representation of an affine hyperplane of rank 𝑞. Even if the procedure shows a successful cluster separation, for nonlinear problems, a nonlinear projection will respond better separating underlying differences between "classes" [START_REF] Hastie | The elements of statistical learning : Data mining, inference, and prediction[END_REF].

Fitting such a model to the data is normally obtained by least squares minimization of the reconstruction error

min 𝑢,𝑉 𝑞, {𝜆 𝑖 } ∑ ‖𝑥 𝑖 -𝑥 𝑖 ′ ‖ 2 𝑁 𝑖=1 = min 𝑢,𝑉 𝑞, {𝜆 𝑖 } ∑ ‖𝑥 𝑖 -𝑢 -𝑉 𝑞 𝜆 𝑖 ‖ 2 𝑁 𝑖=1
.

[ [START_REF] Lynch | CT-definable subtypes of chronic obstructive pulmonary disease: a statement of the Fleischner Society[END_REF][START_REF] Newell | Development of quantitative CT lung protocols[END_REF][START_REF] Iwasawa | Assessment of prognosis of patients with idiopathic pulmonary fibrosis by computer-aided analysis of CT images[END_REF][START_REF] Coxson | Sources of variation in quantitative computed tomography of the lung[END_REF][START_REF] Lynch | Quantitative computed tomography of diffuse lung disease[END_REF][START_REF] Sorensen | Quantitative analysis of pulmonary emphysema using local binary patterns[END_REF][START_REF] Kurugol | Ranking and classification of monotonic emphysema patterns with a multi-class hierarchical approach[END_REF][START_REF] Hansell | CT staging and monitoring of fibrotic interstitial lung diseases in clinical practice and treatment trials: a position paper from the Fleischner Society[END_REF][START_REF] Engeler | Ground-glass opacity of the lung parenchyma: a guide to analysis with high-resolution CT[END_REF][START_REF] Yamada | Airspace enlargement with fibrosis shows characteristic histology and immunohistology different from usual interstitial pneumonia, nonspecific interstitial pneumonia and centrilobular emphysema[END_REF][START_REF] Sluimer | Computer analysis of computed tomography scans of the lung: a survey[END_REF][START_REF] Chien | Detection and classification of interstitial lung diseases and emphysema using a joint morphological-fuzzy approach[END_REF][START_REF] Kim | Development of an automatic classification system for differentiation of obstructive lung disease using HRCT[END_REF][START_REF] Xu | MDCT based 3D texture classification of emphysema and early smoking related lung pathologies[END_REF][START_REF] Zavaletta | Nonlinear histogram binning for quantitative analysis of lung tissue fibrosis in high-resolution CT data[END_REF][START_REF] Depeursinge | Near-affine-invariant texture learning for lung tissue analysis using isotropic wavelet frames[END_REF] It can be demonstrated [START_REF] Hastie | The elements of statistical learning : Data mining, inference, and prediction[END_REF] that we can partially estimate values for 𝑢 and 𝜆 𝑖 obtaining:

𝑢 ̂= 𝑥̅ , [ 2-18 ] 𝜆 𝑖 = 𝑉 𝑞 𝑇 (𝑥 𝑖 -𝑥̅ ) , [ 2-19 ]
where 𝑥̅ is the mean value of the vector 𝑥. This leaves us to find the orthogonal matrix 𝑉 𝑞 by minimizing:

min 𝑉 𝑞 ∑ ‖(𝑥 𝑖 -𝑥̅ ) -𝑉 𝑞 𝑉 𝑞 𝑇 (𝑥 𝑖 -𝑥̅ )‖ 2 𝑁 𝑖=1 . [ 2-20 ]
It is assumed that mean value 𝑥̅ = 0 (otherwise we simply replace the observations by their centered versions 𝑥 𝑖 ̃= 𝑥 𝑖 -𝑥̅ ). Here we can see that the 𝑝𝑥𝑝 matrix 𝑉 𝑞 𝑉 𝑞 𝑇 is a projection matrix which maps each point 𝑥 𝑖 onto its rank-q reconstruction 𝑉 𝑞 𝑉 𝑞 𝑇 𝑥 𝑖 , the orthogonal projection of 𝑥 𝑖 onto

First principal component

Second principal component the subspace spanned by the columns of 𝑉 𝑞 . The solution can be expressed as follows: stack the (centered) observations into the rows of a 𝑁 × 𝑝 matrix 𝑋, and compute the singular value decomposition of 𝑋 [START_REF] Hastie | The elements of statistical learning : Data mining, inference, and prediction[END_REF]:

𝑋 = 𝑈𝐷𝑉 𝑇 . [ 2-21 ]
Here 𝑈 is an 𝑁 × 𝑝 orthogonal matrix (𝑈 𝑡 𝑈 = 𝐼 𝑝 ) whose columns 𝑈 𝑗 are called the left singular vectors; 𝑉 is a 𝑝 × 𝑝 orthogonal matrix (𝑉 𝑡 𝑉 = 𝐼 𝑝 ) with columns 𝑉 𝑗 the right singular vectors, and 𝐷 is a 𝑝 × 𝑝 diagonal matrix, with diagonal elements 𝑑 1 ≥ 𝑑 2 ≥ ⋯ ≥ 𝑑 𝑝 ≥ 0 known as the singular values. For each rank 𝑞, the solution 𝑉 𝑞 to eq. [START_REF] Lynch | CT-definable subtypes of chronic obstructive pulmonary disease: a statement of the Fleischner Society[END_REF][START_REF] Newell | Development of quantitative CT lung protocols[END_REF][START_REF] Iwasawa | Assessment of prognosis of patients with idiopathic pulmonary fibrosis by computer-aided analysis of CT images[END_REF][START_REF] Coxson | Sources of variation in quantitative computed tomography of the lung[END_REF][START_REF] Lynch | Quantitative computed tomography of diffuse lung disease[END_REF][START_REF] Sorensen | Quantitative analysis of pulmonary emphysema using local binary patterns[END_REF][START_REF] Kurugol | Ranking and classification of monotonic emphysema patterns with a multi-class hierarchical approach[END_REF][START_REF] Hansell | CT staging and monitoring of fibrotic interstitial lung diseases in clinical practice and treatment trials: a position paper from the Fleischner Society[END_REF][START_REF] Engeler | Ground-glass opacity of the lung parenchyma: a guide to analysis with high-resolution CT[END_REF][START_REF] Yamada | Airspace enlargement with fibrosis shows characteristic histology and immunohistology different from usual interstitial pneumonia, nonspecific interstitial pneumonia and centrilobular emphysema[END_REF][START_REF] Sluimer | Computer analysis of computed tomography scans of the lung: a survey[END_REF][START_REF] Chien | Detection and classification of interstitial lung diseases and emphysema using a joint morphological-fuzzy approach[END_REF][START_REF] Kim | Development of an automatic classification system for differentiation of obstructive lung disease using HRCT[END_REF][START_REF] Xu | MDCT based 3D texture classification of emphysema and early smoking related lung pathologies[END_REF][START_REF] Zavaletta | Nonlinear histogram binning for quantitative analysis of lung tissue fibrosis in high-resolution CT data[END_REF][START_REF] Depeursinge | Near-affine-invariant texture learning for lung tissue analysis using isotropic wavelet frames[END_REF][START_REF] Xu | Run-length encoding for volumetric texture[END_REF][START_REF] Bartholmai | Quantitative computed tomography imaging of interstitial lung diseases[END_REF][START_REF] Maldonado | Automated quantification of radiological patterns predicts survival in idiopathic pulmonary fibrosis[END_REF] consists of the first 𝑞 columns of 𝑉. The columns of 𝑈𝐷 are called the principal components of 𝑋. The 𝑁 optimal 𝜆 𝑖 in eq. [START_REF] Lynch | CT-definable subtypes of chronic obstructive pulmonary disease: a statement of the Fleischner Society[END_REF][START_REF] Newell | Development of quantitative CT lung protocols[END_REF][START_REF] Iwasawa | Assessment of prognosis of patients with idiopathic pulmonary fibrosis by computer-aided analysis of CT images[END_REF][START_REF] Coxson | Sources of variation in quantitative computed tomography of the lung[END_REF][START_REF] Lynch | Quantitative computed tomography of diffuse lung disease[END_REF][START_REF] Sorensen | Quantitative analysis of pulmonary emphysema using local binary patterns[END_REF][START_REF] Kurugol | Ranking and classification of monotonic emphysema patterns with a multi-class hierarchical approach[END_REF][START_REF] Hansell | CT staging and monitoring of fibrotic interstitial lung diseases in clinical practice and treatment trials: a position paper from the Fleischner Society[END_REF][START_REF] Engeler | Ground-glass opacity of the lung parenchyma: a guide to analysis with high-resolution CT[END_REF][START_REF] Yamada | Airspace enlargement with fibrosis shows characteristic histology and immunohistology different from usual interstitial pneumonia, nonspecific interstitial pneumonia and centrilobular emphysema[END_REF][START_REF] Sluimer | Computer analysis of computed tomography scans of the lung: a survey[END_REF][START_REF] Chien | Detection and classification of interstitial lung diseases and emphysema using a joint morphological-fuzzy approach[END_REF][START_REF] Kim | Development of an automatic classification system for differentiation of obstructive lung disease using HRCT[END_REF][START_REF] Xu | MDCT based 3D texture classification of emphysema and early smoking related lung pathologies[END_REF][START_REF] Zavaletta | Nonlinear histogram binning for quantitative analysis of lung tissue fibrosis in high-resolution CT data[END_REF][START_REF] Depeursinge | Near-affine-invariant texture learning for lung tissue analysis using isotropic wavelet frames[END_REF][START_REF] Xu | Run-length encoding for volumetric texture[END_REF][START_REF] Bartholmai | Quantitative computed tomography imaging of interstitial lung diseases[END_REF] are given by the first 𝑞 principal components (the 𝑁 rows of the 𝑁 × 𝑞 matrix 𝑈 𝑞 𝐷 𝑞 ) [START_REF] Hastie | The elements of statistical learning : Data mining, inference, and prediction[END_REF].

Principal components have many interesting properties, for example, as they are obtained from the eigenvectors of the covariance matrix, they give directions in which the data have maximal variance. So, the linear combination 𝑋𝑣 1 (𝑣 1 the first eigenvector) has the highest variance among all linear combinations of the features; 𝑋𝑣 2 (𝑣 2 the second eigenvector) has the highest variance among all linear combinations satisfying 𝑣 2 orthogonal to 𝑣 1 , and so on. On the other hand, as it is a linear transformation, it is not optimal for non linear problems. Principal curves generalize the principal component line, providing a smooth one-dimensional curved approximation to a set of data points in 𝑅 𝑝 . A principal surface is more general, providing a curved manifold approximation of dimension 2 or more. Kernel PCA [START_REF] Mika | Kernel PCA and de-noising in feature spaces[END_REF] expand the scope of PCA, mimicking what we would obtain if we were to expand the features by non-linear transformations, and then apply PCA in this transformed feature space [START_REF] Hastie | The elements of statistical learning : Data mining, inference, and prediction[END_REF].

Traditional dimensionality reduction techniques such as PCA are linear techniques that focus on keeping the low-dimensional representations of dissimilar datapoints far apart. For highdimensional data that lies on a (near) low-dimensional, non-linear manifold, it is usually more important to keep the low-dimensional representations of very similar datapoints close together, which is typically not possible with a linear mapping. Next section proposes a different approach to the dimensionality reduction task.

t-Distributed Stochastic Neighbor Embedding (t-SNE)

t-Distributed Stochastic Neighbor Embedding is another technique for dimensionality reduction and is particularly well suited for the visualization of high-dimensional datasets. Contrary to PCA, it is a probabilistic technique. t-SNE minimizes the divergence between two distributions: a distribution that measures pairwise similarities of the input objects and a distribution that measures pairwise similarities of the corresponding low-dimensional points in the embedding. The approach analyses the original data (input of the algorithm) and searches how to best represent it using fewer dimensions by matching both distributions. t-SNE is capable of capturing much of the local structure of the high-dimensional data very well, while also revealing global arrangement such as the presence of clusters at several scales. The method is based on Stochastic Neighbor Embedding (SNE) [START_REF] Hinton | Stochastic Neighbor Embedding[END_REF] that converts the high-dimensional Euclidean distances between datapoints into conditional probabilities that represent similarities, but with two key modifications. The cost function used by t-SNE differs from the one used by SNE in two ways:

• t-SNE uses a symmetrized version of the SNE cost function with simpler gradients, • t-SNE uses a Student t-distribution rather than a Gaussian to compute the similarity between two points in the low-dimensional space.

Let us define a high-dimensional dataset of 𝑁 observations 𝑋 = {𝑥 1 , 𝑥 2 , … , 𝑥 𝑁 }. Let 𝑌 = {𝑦 1 , 𝑦 2 , … , 𝑦 𝑁 } be the converted data (dimensionality reduced into two-or three-dimensional data that can be displayed in a scatterplot).

The similarity of datapoint 𝑥 𝑗 to datapoint 𝑥 𝑖 is the conditional probability, p j/i , that 𝑥 𝑖 would pick 𝑥 𝑗 as its neighbor if neighbors were picked in proportion to their probability density under a Student t-distribution centered at 𝑥 𝑖 . For the low-dimensional counterparts 𝑦 𝑖 and 𝑦 𝑗 of the highdimensional datapoints 𝑥 𝑖 and 𝑥 𝑗 , it is possible to compute a similar conditional probability, which we denote by q j/i . t-SNE minimizes the sum of Kullback-Leibler (cost function 𝐶) divergences between the joint probabilities 𝑝 𝑖𝑗 in the high dimensional space and the joint probabilities 𝑞 𝑖𝑗 in the lowdimensional space, over all datapoints using a gradient descent method [START_REF] Maaten | Visualizing data using t-SNE[END_REF]:

𝐶 = 𝐾𝐿(𝑃‖𝑄) = ∑ ∑ 𝑝 𝑖𝑗 log 𝑝 𝑖𝑗 𝑞 𝑖𝑗 𝑗 𝑖 , [ 2-22 ]
where, we set p ii and q ii to zero (because we are only interested in modeling pairwise similarities). We refer to this type of SNE as symmetric SNE, because it has the property that 𝑝 𝑖𝑗 = 𝑝 𝑗𝑖 and 𝑞 𝑖𝑗 = 𝑞 𝑗𝑖 , ∀𝑖, 𝑗. The values for 𝑞 𝑖𝑗 are obtained by means of a Student t-distribution with one degree of freedom:

𝑞 𝑖𝑗 = (1+‖𝑦 𝑖 -𝑦 𝑗 ‖ 2 ) -1 ∑ (1+‖𝑦 𝑘 -𝑦 𝑙 ‖ 2 ) -1 𝑘≠𝑙 , [ 2-23 ]
and the joint probabilities 𝑝 𝑖𝑗 in the high-dimensional space to be the symmetrized conditional probabilities, that is:

𝑝 𝑖𝑗 = p j/i +p i/j 2𝑁
. This ensures that

∑ 𝑝 𝑖𝑗 > 1 2𝑁 𝑗 , ∀𝑥 𝑖 . [ 2-25 ]
As a result of which, each datapoint 𝑥 𝑖 makes a significant contribution to the cost function. The main advantage of the symmetric version of SNE is the simpler form of its gradient (needed to minimize the cost function 𝐶), which is faster to compute. The t-SNE gradient can be then derivated to [START_REF] Maaten | Visualizing data using t-SNE[END_REF]:

𝛿𝐶 𝛿𝑦 𝑖 = 4 ∑ (𝑝 𝑖𝑗 -𝑞 𝑖𝑗 ) (1 + ‖𝑦 𝑖 -𝑦 𝑗 ‖ 2 ) -1 (𝑦 𝑖 -𝑦 𝑗 ) 𝑗 . [ 2-26 ]
The implementation of t-SNE is computationally quite heavy and therefore there are some limitations to the use of this technique. However, in [START_REF] Van Der Maaten | Accelerating t-SNE using tree-based algorithms[END_REF], an algorithm to substantially accelerate t-SNE tree-based algorithms is proposed (that makes it possible to learn embeddings of data sets with millions of points in practice). This paper develops variants of the Barnes-Hut algorithm and of the dual-tree algorithm that approximate the gradient used for learning t-SNE embeddings in 𝑂(𝑁𝑙𝑜𝑔𝑁).

Feature representation

There have been many different efforts using various methods (and a wide range of features) to robustly characterize the discernible patterns in CT scans over the last several decades using two major perspectives: hand-designed heuristics or automatic learning, both detailed in the following.

Feature space engineering (hand-design)

As we mentioned in the previous section, most part of the literature concentrates on finding "good" textural features for lung texture analysis using similar information -the spatial periodicity and scales contained in the images. Another alternative exploits intensity-based signatures to extract discriminating information. The different useful features are presented and discussed in the following.

Texture-based features

Texture-based features target information related to local distribution of gray level patterns in the image. They include GLCM, RLE, and LBP.

• Gray-level co-occurrence matrices (GLCM)

GLCM is a second-order statistics for characterizing the occurrence of gray-level combinations in pairs of spatially related pixels [START_REF] Xu | MDCT based 3D texture classification of emphysema and early smoking related lung pathologies[END_REF], [START_REF] Zavaletta | Nonlinear histogram binning for quantitative analysis of lung tissue fibrosis in high-resolution CT data[END_REF]. In a 2D image of dimensions 𝑁 × 𝑀, GLCM counts the number of occurrences of a sequence of gray level values of a defined length and along a defined direction given by offsets in the x and y directions: 𝛥𝑥 and 𝛥𝑦. GLCM are not rotation invariant. As an approximation, several GLCM matrices can be computed with different 𝛥𝑥/𝛥𝑦 ratios to span 180 degrees. Another drawback pointed out in [START_REF] Depeursinge | Near-affine-invariant texture learning for lung tissue analysis using isotropic wavelet frames[END_REF] is that the size of the matrices increases exponentially with the number of gray levels used in the image. So analyzing 16 bit images from a CT scanner, with a typical resolution of 512×512 would result impractical because the number of possible cooccurrences is huge, increasing the computational complexity and raising data storing issues. A compromise solution could be to reduce the number of levels used, but this would result in an important loss of relevant information that is unacceptable in our context. This can be handled by compressing the grayscale or by considering an alternative representation in terms of sums and differences. A basic example of the GLCM matrix computation is shown in Figure 27.

• Run length encoding (RLE)

RLE is a second-order statistics to capture the coarseness of a texture in specified directions. A run is a set of consecutive pixels displaying the same gray-level intensity along a specified linear direction [START_REF] Xu | Run-length encoding for volumetric texture[END_REF]. For 2D analysis, an element at location (𝑖, 𝑗) in a run-length matrix represents the number of runs associated with gray-level intensity 𝑖 and length of run j along a specific direction.

Note that, in general, there are 4 directions considered, 0°, 45°, 90°, and 135°. Generally, 11 features are extracted from the run-length matrices: short run emphasis, long run emphasis, low gray-level run emphasis, high gray-level run emphasis, short run low gray-level emphasis, short run high graylevel emphasis, long run low gray-level emphasis, long run high gray-level emphasis, gray-level nonuniformity, run length non-uniformity, and run percentage. For 3D analysis, in the statistics, fine textures with similar gray-level intensities contain more short runs, while coarse textures with significantly different gray-level intensities produce more long runs. Figure 27. An example of GLCM for Δx = 1 and Δy = 0 of an image with 8 gray-levels. An element at location (𝒊, 𝒋) in the cooccurrence matrix specifies the frequency of occurrence of two gray-levels separated by a given distance and a given direction from each other [START_REF] Depeursinge | Near-affine-invariant texture learning for lung tissue analysis using isotropic wavelet frames[END_REF].

• Local binary patterns (LBPs)

The LBP has been found to be a powerful feature for texture classification. In [START_REF] Sorensen | Quantitative analysis of pulmonary emphysema using local binary patterns[END_REF], joint LBP and intensity histograms are used for characterizing regions of interest comprising three classes (normal tissue, centrilobular emphysema, and paraseptal emphysema), with a 95% classification accuracy. Even if the specific application is emphysema quantification, this idea is also applicable to other lung disease patterns. LBPs are obtained by thresholding samples in a local neighborhood with respect to the center pixel intensity and is given by

𝐿𝐵𝑃 (𝐱; 𝑅, 𝑃) = ∑ 𝐻(𝐼(𝑥 𝑝 ) -𝐼(𝐱))2 𝑝 𝑃-1 𝑝=0 , [ 2

Intensity-based signatures

Complementary to the characterization of spatial dependencies, the distributions of the colors and intensities within the image carry out relevant information for texture and object recognition. Densitometric measurements in CT imaging could be resumed by the following parameters calculated over the corresponding histogram distribution 𝑋:

• Mean lung density (MLD): weighted intensity mean value;

• Standard deviation from the histogram curve;

• Skewness: this coefficient is a measure for the degree of symmetry in the variable distribution. In probability theory and statistics, skewness is a measure of the asymmetry of the probability distribution of a real-valued random variable about its mean. The skewness of a random variable X is the third standardized moment 𝛾 1 , defined as [START_REF] Doane | Measuring skewness: a forgotten statistic[END_REF] 𝑆𝑘𝑒𝑤𝑛𝑒𝑠𝑠

[𝑋] = 𝛾 1 = 𝐸 [( 𝑋-𝜇 𝜎 ) 3 ] = 𝜇 3 𝜎 3 = 𝐸[(𝑋-𝜇) 3 ] (𝐸[(𝑋-𝜇) 2 ]) 3 2 ⁄ , [ 2-28 ]
where μ is the mean, σ is the standard deviation, E is the expectation operator, 𝜇 3 is the third central moment.

• Kurtosis: this coefficient is a measure for the degree of tailedness in the variable distribution.

The kurtosis is the fourth standardized moment, defined as [START_REF] Decarlo | On the meaning and use of kurtosis[END_REF]:

𝐾𝑢𝑟𝑡𝑜𝑠𝑖𝑠[𝑋] = 𝛾 2 = 𝐸 [( 𝑋-𝜇 𝜎 ) 4 ] = 𝜇 4 𝜎 4 = 𝐸[(𝑋-𝜇) 4 ] (𝐸[(𝑋-𝜇) 2 ]) 2 , [ 2-29 ]
where again, μ is the mean, σ is the standard deviation, E is the expectation operator, and 𝜇 4 is the fourth central moment. To have a glance of these parameters, the next figure presents several examples. Grayscale values in medical images contain valuable information for the characterization of objects and textures, and are complementary data to texture features [START_REF] Iwasawa | Assessment of prognosis of patients with idiopathic pulmonary fibrosis by computer-aided analysis of CT images[END_REF]. The direct measure is the relative amount of lung parenchyma pixels that have attenuation values between certain thresholds (in a given range of Hounsfield units). Measures based on the attenuation histogram disregard the information related to the morphology of the structures, such as shape and size distribution of bullae in emphysema, for example [START_REF] Sorensen | Quantitative analysis of pulmonary emphysema using local binary patterns[END_REF].

Automatic Feature learning

A potentially more interesting scheme is to rely as much as possible on learning the feature extractor itself, i.e., not to develop an algorithm that solves a problem, but instead learns to solve the problem from a big set of examples/experiences with a process of learning. We explore this alternative approach as it was demonstrated that better pattern recognition systems can be built by relying more on automatic learning and less on hand-designed heuristics [START_REF] Lecunn | Gradient-based learning applied to document recognition[END_REF].

Recently, the use of deep learning techniques (artificial models of human visual system, best system known to solve visual tasks), especially Convolutional Neural Networks (CNNs), has shown great improvements with respect to traditional methods (handcrafted features seen in previous section), emerging as the new state of the art for visual tasks. Deep learning methods have the advantage (over traditional CAD approaches) that the features on which relies the classification are trained directly from the image data. Such methods tend more likely to capture the intrinsic features from the imaging dataset used for training than when using handcrafted rules [START_REF] Lecunn | Gradient-based learning applied to document recognition[END_REF]. CNNs particularly have proven to succeed in image classification tasks [START_REF] Krizhevsky | Imagenet classification with deep convolutional neural networks[END_REF] [23] [START_REF] Roth | Improving Computer-Aided Detection Using Convolutional Neural Networks and Random View Aggregation[END_REF]. CNN features are trained from the data in a fully supervised fashion to learn increasingly complex transformations of the input (corresponding to different levels of abstraction) to capture invariances. Importantly, CNNs pretrained on such large datasets have been shown to contain general-purpose feature extractors, transferrable to many other domains (transfer learning techniques).

There are, naturally, other artificial neural network designs available in the literature briefly described in the following (refer to Figure 30 for simplified schemas):

• Recurrent Neural Networks (RNN)

RNNs exploit sequential information (they are not stateless). In a traditional neural network we assume that all inputs are independent of each other. But in RNNs, neurons are fed information not just from the previous layer but also from themselves from the previous pass (directed cycles in their connection graphs). RNNs are called recurrent because they perform the same task for every element of a sequence, but with the addition of also considering the previous computations. RNNs have "memory" which captures information about what has been calculated so far. RNNs have shown great success in many Natural Language Processing (NLP) tasks. The most commonly used type of RNNs are Long Short Term Memory (LSTM). In general, recurrent networks are a good choice for advancing or completing information. However complicated dynamics can make them difficult to train.

• Autoencoders (AE)

An autoencoder is a neural network that is trained to attempt to copy its input to its output. Internally, it has a hidden layer that describes a code used to represent the input. They compress information automatically, hence the name, due to smaller hidden layers than the input and output layers.

• Restricted Boltzmann Machines (RBMs)

Boltzmann machines (BMs) can be defined as bidirectionally connected networks of stochastic processing units, which can be interpreted as neural network models. A BM can be used to learn important aspects of an unknown probability distribution based on samples from this distribution. In general, this learning process is difficult and time-consuming. However, the learning problem can be simplified by imposing restrictions on the network topology, which leads us to restricted Boltzmann machines.

• Deep Belief Networks (DBNs)

Is the name given to stacked architectures of mostly RBMs or AEs. These networks have been shown to be effectively trainable stack by stack, where each AE or RBM only must learn to encode the previous network. This technique is also known as greedy training, where greedy means making locally optimal solutions to get close to an overall answer. DBNs can be trained through contrastive divergence or back-propagation and learn to represent the data as a probabilistic model, just like regular RBMs or AEs. Once trained or converged to a stable state through unsupervised learning, the model can be used to generate new data. If trained with contrastive divergence, it can classify existing data as the neurons have been taught to look for different features.

• Generative Adversarial Networks (GANs)

GANs are constructed with two networks working together (although often a combination of Feed Forward Network and CNNs), with one tasked to generate content and the other to judge content. The discriminating network receives either training data or generated content from the generative network. How well the discriminating network was able to correctly predict the data source is then used as part of the error for the generating network. This creates a form of competition where the discriminator is getting better at distinguishing real data from generated data and the generator is learning to become less predictable to the discriminator. All these presented architectures can be seen as tools to approximate complex functions, discovering proper features automatically in a hierarchical manner from fine to abstract [START_REF] Zhou | Deep Learning for Medical Image Analysis[END_REF]. However, the advantage of using CNNs is the reduced number of parameters that shape the network due to restrictions in the connection design between layers that also leads to a highly adapted network to the specific visual semantic interpretation task. These conditions will be explained in detail in section 2. Generative Adversarial Network (GAN)

Classifiers

Once features representation of input images is achieved, an algorithm is used to detect and create boundaries among the several classes defined by clusters in the transformed space. Briefly, we enumerate in this section the best-established methods in the literature. Most of these methods are machine learning based. Examples include naive Bayes classifiers, fuzzy logic, support vector machines, neural networks and decision tree classifiers. Most of these methods are adjusted in a fully supervised way, where the underlying parameters are tuned up during a learning process and optimized in order to draw decision boundaries that minimize the errors of categorization of instances with known class labels. For those cases, the model generated by a learning algorithm should both fit the input data well and correctly predict the class labels of records it has never seen before. Therefore, a key objective of the learning algorithm is to build models with good generalization capability; i.e., models that accurately predict the class labels of previously unknown samples.

Bayesian Classifier

The Bayes classifier is based on a probability model and assigns the class which has the maximum estimated posterior probability to the feature vector extracted from the ROI. The posterior probability P(w i /v) of a class w i given a feature vector 𝑣 (that represents the occurrence probability of the category 𝑤 𝑖 when the feature vector 𝑣 has been measured) is determined using Bayes' theorem:

𝑃(𝑤 𝑖 /𝑣) = 𝑃(𝑣/𝑤 𝑖 ).𝑃(𝑤 𝑖 ) 𝑃(𝑣)
, where 𝑃(𝑣/𝑤 𝑖 ) is the probability of meaning 𝑣 given 𝑤 𝑖 , 𝑃(𝑤 𝑖 ) and 𝑃(𝑣) are the probability that the category 𝑤 𝑖 and vector 𝑣 appears in all data, respectively. Figure 31 shows an example representing different 𝑣 values and the corresponding posterior probabilities 𝑃(𝑤 𝑖 /𝑣) of a two-class problem. The feature vector 𝑣 would be assigned to the category 𝑤 𝑖 , if 𝑃(𝑤 𝑖 /𝑣) is maximum over all others; that is

𝑃(𝑤 𝑖 /𝑣) > 𝑃(𝑤 𝑗 /𝑣) for all 𝑗 ≠ 𝑖 . [ 2-31 ]
Bayesian decision theory is a fundamental statistical approach (normally applied to classification problem) based on quantifying the tradeoffs between various classification decisions using probability and the costs that accompany such decisions [START_REF] Duda | Pattern classification[END_REF].

The naive Bayesian classifier is a simplified probabilistic classifier based on Bayesian theory which assumes that attributes are all independent. Thus, under the naive assumption, the classifying function is defined as eq. , but 𝑃(𝑣/𝑤 𝑖 ) should be replaced by

𝑃(𝑣/𝑤 𝑖 ) = ∏ 𝑃(𝑣 𝑗 /𝑤 𝑖 ) 𝑑 𝑗=1
, where 𝑣 𝑗 denotes the 𝑗-th feature, and 𝑑 is the number of features. The simplicity of the method allows good performance with small training sets. Indeed, by building probabilistic models, the naive Bayesian classifier is robust to outliers (i.e. feature vectors that are not representative of the class to which they belong). Moreover, it creates soft decision boundaries, which has the effect to avoid overtraining. However, the arbitrary choice of the distribution model for estimating the probabilities along with the lack of flexibility of the decision boundaries results in limited performance for complex multiclass configurations [START_REF] Farag | Detection and recognition of lung abnormalities using deformable templates[END_REF] [31].

Fuzzy logic

Fuzzy theory provides the membership function operating over the range of real numbers [0, 1] instead of the conventional elements of binary logic, such as 0/1, yes/no, and true/ false. The fuzzy system is a knowledge-based system and consists of the fuzzy IF-THEN rules. The fuzzy IF-THEN rules induced from human experts or based on domain knowledge are IF-THEN statements in which some words are used to characterize a continuous membership function. According to different rule principles, fuzzy systems are divided into three types, the pure fuzzy system, the Takagi-Sugeno-Kang (TSK) fuzzy system [START_REF] Sugeno | Structure identification of fuzzy model[END_REF][START_REF] Takagi | Fuzzy identification of systems and its applications to modeling and control[END_REF], and the fuzzy system with fuzzification and defuzzification [START_REF] Zhang | Atlas-driven lung lobe segmentation in volumetric X-ray CT images[END_REF].

A pure fuzzy system is composed of a fuzzy inference unit combined with a fuzzy rule base, which is the collection of fuzzy IF-THEN rules. The basic architecture of the pure fuzzy system is shown in Figure 32. Suppose the death rate has sharp increase in people with 65-75 years old. In the pure fuzzy system, fuzzy IF-THEN rules are represented in the following form:

IF the age of a person is within 65-75 years old THEN with high death rate

With the fuzzy rule, the fuzzy inference unit maps the fuzzy set from the input space U to the output space V. The main problem of the pure fuzzy system is that its inputs and outputs are fuzzy sets, which are represented by words in natural languages, rather than real-valued variables used in engineering systems. A TSK fuzzy system is proposed to solve the problem of the pure fuzzy system. The basic architecture of the TSK fuzzy system is shown in Figure 33. In the TSK system, the form of fuzzy rules is represented in a mathematical formulation as IF a person is 𝑥 years old, and 𝑥 is within 65-75 THEN the death rate of the person is 𝑦 = 𝑐𝑥 , where 𝑐 is a constant. This change makes inputs and outputs to be real-valued variables and the fuzzy inference unit in the TSK system is a weighted average of the values in the THEN parts of fuzzy rules. However, the TSK system has two problems: the lack of representing human knowledge and the difficulty of applying different principles in fuzzy logic. A fuzzy system with fuzzifier and defuzzifier transforms input real-valued variables into a fuzzy set by using a fuzzifier, and the output of this system is transformed from a fuzzy set into realvalued variables by using a defuzzifier. The basic architecture of the fuzzy system with fuzzifier and defuzzifier is shown in Figure 34. Under this architecture, inputs and outputs are represented by realvalued variables so that it is easier to perform analysis in mathematical functions; on the other hand, the human knowledge can be constructed in engineering applications.

Support Vector Machine (SVM)

The SVM is a useful supervised learning method used for classification and regression in data analysis and pattern recognition because it finds effectively the decision boundaries among the different classes. SVMs have shown to be effective to categorize texture in wavelet feature spaces [START_REF] Depeursinge | Near-affine-invariant texture learning for lung tissue analysis using isotropic wavelet frames[END_REF] and in particular lung tissue [START_REF] Depeursinge | Comparative performance analysis of state-of-the-art classification algorithms applied to lung tissue categorization[END_REF].

Generally, there are many hyperplanes that might separate a set of points into two categories. The SVM selects the best one which is with the maximum margin (distance from the hyperplane to the nearest point on each side) as it is shown in Figure 35 (a)). A set of points is separated into two categories, 𝑉 1 and 𝑉 2 , by the maximum-margin hyperplane. However, most classification tasks require a nonlinear hyperplane (curve) for a precise separation. The SVM maps points into a higher dimension space through an appropriate nonlinear mapping function (the kernel function) which is produced from training data. In this space, the points are rearranged and thus the maximum-margin hyperplane can be found instead of the curve for a full separation.

In summary, SVMs allow training generalizable, nonlinear classifiers in high-dimensional spaces using a small training set. This is enabled through the selection of a subset of vectors (called the support vectors), which well characterize the true boundaries between the classes. However, for a multiclass problem, a modification has to be added to the algorithm (one vs all, etc.).

Artificial Neural Networks

This section describes the alternative computation paradigm to the usual one (based on a programmed instruction sequence) that is the artificial neural networks. These networks are large collections of interacting entities with emergent properties and are based on real neurons in the brain [START_REF] Hertz | Introduction to the Theory of Neural Computation[END_REF] illustrated in the next figure (it is not a strict model, as the real neuron involves many complications, here omitted). An example of layered feed-forward networks (called perceptrons) is shown in Figure 37. These systems can design themselves to solve a given task (hetero-association task that includes classification problems) with relatively little external guidance: a list of training set of correct inputoutput pairs as example. When we apply one of the training inputs to the network we can compare the network output or prediction to the target or correct output, and then change the connection strengths 𝑤 𝑖𝑗 to minimize the difference. Note that the learning rule starts from a general first guess at the weights values (random) and then makes successive improvements. Each layer is a new representation of the input as transformations take place in each step through the network. Without the hidden layers, networks are very limited in the in-out mappings they can model. Note that learning the weights of the hidden layers is equivalent to learning features.

Decision trees

Decision Trees (DTs) are a non-parametric supervised learning method used for classification and regression. It is a model that predicts the value of a target variable by learning simple decision rules inferred from the data features. Tree-based methods (tree-like graph) partition the feature space into a set of rectangles, and then fit a simple model (like a constant) in each one. A decision tree is a flowchart-like structure in which each internal node represents a test (can be seen as an experimental outcome of a random variable) on an attribute, each branch represents the outcome of the test, and each leaf node represents a class label (decision taken after computing all attributes). The paths from root to leaf represent classification rules. In a decision tree, each leaf node is assigned a class label. The non-terminal nodes, which include the root and other internal nodes, contain attribute test conditions to separate records that have different characteristics. Classifying a new test sample is straightforward once a decision tree has been constructed. Starting from the root node, we apply the test condition to the record and follow the appropriate branch based on the outcome of the test. This will lead us either to another internal node, for which a new test condition is applied, or to a leaf node. The class label associated with the leaf node is then assigned to the record [START_REF] Hastie | The elements of statistical learning : Data mining, inference, and prediction[END_REF] [START_REF] Safavian | A survey of decision tree classifier methodology[END_REF].

The training set can be split with recursive steps during the tree-growing process selecting an attribute test condition to divide the records into smaller subsets (Figure 38). This task is done by defining a measure for selecting the best split among several possibilities. These measures are defined in terms of the class distribution of the records before and after splitting, and they are based on the degree of impurity of the child nodes. The smaller the degree of impurity, the more skewed the class distribution.

To define the impurity measures, let 𝑚 denotes a node in the tree, representing a region 𝑅 𝑚 with 𝑁 𝑚 observations. We assign the observation in node 𝑚 to class (𝑚) = arg 𝑚𝑎𝑥 𝑘 𝑝̂𝑚 𝑘 , representing the majority class in node 𝑚. The target classification outcome 𝑌 is considered to take values 1,2,…,K. The proportion of class k observations in node 𝑚, 𝑝̂𝑚 𝑘 , is defined as:

𝑝̂𝑚 𝑘 = 1 𝑁 𝑚 ∑ 𝐼(𝑦 𝑖 = 𝑘)
𝑥 𝑖 ∈𝑅 𝑚 where 𝑥 𝑖 are the inputs, and 𝐼(𝑆) is the indicator of the set 𝑆.

The following impurity measures, Gini index, misclassification error, and cross-entropy (Figure 39) are defined as follows [START_REF] Hastie | The elements of statistical learning : Data mining, inference, and prediction[END_REF]: As mentioned, a nice property of tree models is that they are nonparametric approaches for building classification models, which means that they do not require any prior assumptions regarding the type of probability distributions satisfied by the class and other attributes. In addition, techniques developed for constructing decision trees are computationally inexpensive, making it possible to quickly construct models even when the training set size is very large. Furthermore, once a decision tree has been built, classifying a test record is extremely fast, with a worst-case complexity of 𝑂(ℎ), where h is the maximum depth of the tree. Finally, DTs are considered white box models. If a given situation is observable in a model, the explanation for the condition is easily explained by boolean logic. By contrast, in a black box model (e.g., in an artificial neural network), results may be more difficult to interpret. On the other hand, as practical decision-tree learning algorithms are based on heuristic algorithms such as the greedy algorithm, where locally optimal decisions are made at each node, such algorithms cannot guarantee to return the globally optimal decision tree. This can be mitigated by training multiple trees in an ensemble learner, where the features and samples are randomly sampled with replacement. Another drawback is that decision trees tend to overfit on data with a large number of features. Getting the right ratio of samples to number of features is important, since a tree with few samples in high dimensional space is very likely to overfit.

𝐺𝑖𝑛𝑖 𝑖𝑛𝑑𝑒𝑥 = ∑ 𝑝̂𝑚 𝑘 𝑝̂𝑚 𝑘′ = ∑ 𝑝̂𝑚 𝑘 𝐾 𝑘=1 (1 -𝑝̂𝑚 𝑘 ) 𝑘≠𝑘′ , [ 2-34 ] 𝑀𝑖𝑠𝑐𝑙𝑎𝑠𝑠𝑖𝑓𝑖𝑐𝑎𝑡𝑖𝑜𝑛 𝑒𝑟𝑟𝑜𝑟 = 1 𝑁 𝑚 ∑ 𝐼(𝑦 𝑖 ≠ 𝑘(𝑚)) = 1 -𝑝̂𝑚 𝑘 (𝑚) 𝑖∈𝑅 𝑚 , [ 2
Random forests are a substantial modification of bagging that builds a large collection of decorrelated trees, and then averages them. The idea of bagging is to average many noisy but approximately unbiased models, and hence reduce the variance [START_REF] Hastie | The elements of statistical learning : Data mining, inference, and prediction[END_REF].

Comparative performance analysis of SOA classification systems applied to lung texture categorization

The question that arises is which one of the texture features described characterizes the patterns the best and is the most adaptable to the needs of lung tissue analysis in MSCT imaging. Highly flexible texture modeling is required to catch subtle texture signatures of a given lung tissue pattern. In particular, invariance of the texture descriptors to affine transformations (i.e., translation, rotation, and scale) is desirable to obtain a system that is able to learn any texture appearance independently of localizations, orientations, or sizes. In this section we summarize the best suited methods found in the literature to compare performance of classifiers and goodness of features.

Here, a few of the salient (in terms of performance) methods combining approaches mentioned in §2.1 and §2.3 will be compared and discussed with the purpose to decide the research direction to follow in this work, either toward a handcrafted solution, or to an automatic feature learning.

For the first class of solutions, searching to find the best handcrafted features to represent the ILD patterns, the idea consists of transforming the raw data into a new space where a classifier will be certainly more accurate. Two main works show promising results using different both features and decision making criteria:

• In [START_REF] Depeursinge | Near-affine-invariant texture learning for lung tissue analysis using isotropic wavelet frames[END_REF] a blockwise classification with texture (either discrete wavelet frames, DWF or quincunx wavelet frames, QWF) and gray-level histogram (GLH) near-affine-invariant features was presented. The accuracies obtained by applying an SVM classifier showed its potential for effectively learning from the proposed hybrid feature space. The same author, in [START_REF] Depeursinge | Comparative performance analysis of state-of-the-art classification algorithms applied to lung tissue categorization[END_REF], performed a comparative performance analysis of classifiers applied to the same problem. The work compares five common classifiers with respect to their ability to categorize lung tissue patterns in high-resolution computed tomography (HRCT) images of patients affected with interstitial lung diseases (ILD). The evaluated classifiers are naive Bayes, k-nearest neighbor, J48 decision trees, multilayer perceptron, and support vector machines. SVM reached best values for each metric and allowed a high precision on testing sets.

• In [START_REF] Hansell | CT staging and monitoring of fibrotic interstitial lung diseases in clinical practice and treatment trials: a position paper from the Fleischner Society[END_REF], a fast grey-level histogram-based classification method was proposed. The detection and quantification of pulmonary parenchyma in this work is based on histogram signature mapping techniques trained through expert radiologist consensus on assessment of pathologically confirmed data sets obtained from the Lung Tissue Research Consortium (LTRC) [START_REF] Bartholmai | Quantitative computed tomography imaging of interstitial lung diseases[END_REF]. As pre-processing supervised training phase, multiple small volumes of interest were selected by expert radiologists, from HRCT scans of 14 subjects with proven pathological diagnosis of diffuse pulmonary disease or control subjects without ILD. These samples were then used to determine canonical histogram signatures for each of the classes of visual abnormality with automatic cluster affinity techniques. Finally, those signatures of each of the visual classes were used for the volumetric classification of the HRCT data of test subjects showing that the output characterization and quantification of lung parenchyma on HRCT correlates with visual assessment by expert radiologists, and that quantitative shortterm volumetric longitudinal changes on serial HRCT correlate with IPF mortality (quantitative features of ILD on HRCT may represent an accurate and reproducible biomarker for IPF) [START_REF] Maldonado | Automated quantification of radiological patterns predicts survival in idiopathic pulmonary fibrosis[END_REF].

We can find the second class of solutions in the state of the art that are associated with supervised learning techniques. The idea also consists of transforming the raw data into a new space, but this time the mapping function will be learnt to fit a training dataset. Two main works from this defined category, using artificial neural networks, show best outcome for ILD texture categorization:

• Texture analysis is recognized as being particularly relevant to describe the visual content of medical images that is often not characterized by well-defined objects but has well-defined visual aspects. A particular class of ANNs, the convolutional neural networks, CNNs, may "understand" the high-level content in images (useful for representing the shape of an object): convolutional layers are similar to non-linear filter banks, and the fully connected layers capture their spatial layout. To capture the style in images (or texture information without considering the global arrangement) [START_REF] Anthimopoulos | Lung Pattern Classification for Interstitial Lung Diseases Using a Deep Convolutional Neural Network[END_REF] proposed a CNN designed to extract the low-level textural features of the lung tissue, showing high classification performance for a train-validation-test patch scheme. This was achieved by reducing the size of convolutional kernels (minimum 2x2 pixels) and pooling only at the end of the network with respect to a traditional CNN.

• [START_REF] Van Tulder | Learning features for tissue classification with the classification restricted boltzmann machine[END_REF] shows how restricted Boltzmann machines (RBMs) can be used to learn features that are especially suited for texture-based tissue classification, evaluating its performance on ILD detection problem. High classification accuracy was obtained using these learnt filters. It can be seen that the highest performance is achieved by CNN, which is also supported by the outstanding results achieved by these ANNs in computer vision tasks in the last few years [START_REF] Lecunn | Gradient-based learning applied to document recognition[END_REF] [22] [START_REF] Roth | Deep convolutional networks for pancreas segmentation in CT imaging[END_REF] [24] [START_REF] Tarando | Increasing CAD system efficacy for lung texture analysis using a convolutional network[END_REF] with respect to traditional methods. The human visual system modeling has many features that are desirable in artificial systems, especially dealing with all variations and distortions of the data while working on ILDs automatic classification systems. That is, CNN is robust and fault tolerant (the highly connectivity of the network means that errors in a few terms will probably be inconsequential), flexible (it can adjust to new environments by learning), it can deal with information that is fuzzy, probabilistic, noisy and inconsistent, and it is highly parallel. All these properties lead to the ability to generalize properly to new situations and make CNN an attractive candidate for the task of classification of pathologic lung texture; this is the hypothesis that we investigate further in this work.

CNN also introduces a change in the investigation formalism: instead of having to specify every detail of a calculation, we simply have to compile a training set of representative examples. This is the best way to treat a problem when appropriate rules are very hard to know in advance or inordinately expensive and time consuming. However, the choice of the network architecture is important for the problem to address and will influence the performance achievable.

In the next section, a detailed description of CNN architectures is presented together with the basic theory on which relies the proposed approach.

Focus on CNN-based approaches

While more automatic learning is beneficial, no learning technique can succeed without a minimal amount of prior knowledge about the task. In the case of multilayer artificial neural networks (ANNs), a good way to incorporate knowledge is to fit its architecture to the task. To inject prior knowledge in any ANN we can modify the design of the network using knowledge (about invariance for example): the connectivity (to be local or sparse), the pooling, weight constraints (sharing), etc. In the following, the design of the specific type of ANN, a convolutional neural network (CNN) is described to demonstrate the gains with respect to a fully connected ANN.

CNN design

CNN are biologically-inspired variants of multilayer perceptrons (MLPs). We know the visual cortex contains a complex arrangement of cells. These cells are sensitive to small sub-regions of the visual field, called a receptive field. The sub-regions are tiled to cover the entire visual field. These cells act as local filters over the input space and are well suited to exploit the strong spatially local correlation present in natural images.

• Local connectivity (sparse) across space (local receptive fields)

Images have a strong 2D local structure: variables (or pixels) that are spatially nearby are highly correlated. Local correlations are the reasons for the well-known advantages of extracting and combining local elementary visual features before recognizing spatial objects, because configurations of neighboring variables can be classified into a small number of categories (e.g., edges, corners, etc.). CNNs force the extraction of local features by restricting the receptive fields of hidden units to be local. Each unit in a layer receives inputs from a set of units located in a small neighborhood (receptive fields) in the previous layer. These features are then combined by the subsequent layers in order to detect higher order features. The receptive fields of contiguous units in a feature map are centered on corresponding contiguous units in the previous layer. Therefore, receptive fields of neighboring units overlap (see Figure 40). The architecture thus ensures that the learnt "filters" produce the strongest response to a spatially local input pattern. In other words, the inputs of hidden units in layer m are from a subset of units in layer 𝒎 -𝟏, units that have spatially contiguous receptive fields [START_REF]Theano tutorial on CNN[END_REF].

However, as shown above, stacking many such layers leads to (non-linear) "filters" that become increasingly "global" (i.e. responsive to a larger region of pixel space). For example, the unit in hidden layer 𝑚 + 1 in the figure can encode a non-linear feature of width 5 (in terms of pixel space).

• Sharing weights across space (W replication, replicated local feature detector)

Because same objects may show up on different pixels, a feature (edge, etc.) detector useful in one place in the image is likely to be useful somewhere else. In CNNs, this knowledge can be applied by forcing a set of units, whose receptive fields are located at different places on the image, to have identical weight vectors (and so, each filter is replicated across the entire visual field). These replicated units share the same parameterization (weight vector and bias) and thus search for the same feature in the receptive field (Figure 41). Units in a layer are organized in planes within which all the units share the same set of weights. The set of outputs of the units in such a plane is called a feature map (see section 2.5.2.1 for details of how feature maps are built). Units in a feature map are all constrained to perform the same operation on different parts of the image.

Weight sharing or replication increases learning efficiency by greatly reducing the number of free parameters that needs to be learnt (more efficient, and less data and computation time requiered). The constraints on the model enable CNNs to achieve better generalization on vision problems. We are going to have many feature maps to allow each patch of the image to be represented in several ways. Replicating feature detectors also achieves equivariant activities (neural activation identical for a given representation) and invariant knowledge (invariance in the weights).

Figure 41. Example of weight replication illustrated on 3 hidden units belonging to the same feature map. Weights of the same color are shared, thus, constrained to be identical. Replicating units in this way allows for features to be detected regardless of their position in the visual field [START_REF]Theano tutorial on CNN[END_REF].

• Local translation invariance (pooling)

Distortions or shifts of the input can cause the position of salient features to vary. We achieve small amount of translation invariance in units' activities by pooling together the output of replicated feature detectors. Averaging (or taking the maximum value) for neighboring replicated detectors at each level of the deep net reduces the number of inputs to the next layer allowing us to have more different feature maps (see section 2.5.2.2 for details of this particular layer, and Figure 42 for an illustration). The drawback is that after several levels of this pooling we lose precise information about the position where things are in the image. Going forward in the CNN hierarchy allows getting features that are much more complicated (abstract) and more invariant to position.

CNN layers

CNNs are formed by mainly three different types of layers detailed in the following: convolutional layer, pooling layer, and fully connected layer. Optionally a SoftMax layer allows the network to present a probability of obtaining each class given a test input instead of the class itself.

Convolutional layer

The convolutional layers (normally a deep network contains several of them) are conformed by different feature maps, where each one of them corresponds to a different feature computed at different locations of the input. A feature map is obtained by repeated application of a function across sub-regions of the entire image, in other words, by convolution of the input image with a linear filter (defined by weights w ij ), adding a bias (parameter 𝑏 in eq. ) term and then applying a non-linear function, Rectified Linear Unit (RELU) in this case (Figure 43-Figure 45). The output or the activity 𝑦 of a given unit from a given feature map is a nonlinear function of the total input and it is calculated as follows: -37 ] where RF denotes the receptive field,

𝑧 = 𝑏 + ∑ x i . w ij i∈RF , [ 2
y = RELU(z) = { z if z > 0 0 otherwise [ 2-38 ]
Note that the point-wise derivative for RELU can be easily computed by:

𝑑𝑦 𝑑𝑧 = { 1 𝑖𝑓 𝑧 > 0 0 𝑖𝑓 𝑧 ≤ 0 [ 2-39 ]
RELU is actually not differentiable, but it has subdifferential [0,1]. Any value in that interval can be taken as a subderivative, and can be used during the learning process.

Layer m

Layer m-1 Figure 44. Diagram showing how the output y of a given unit j (that belongs to a given feature map in layer m) is calculated from its inputs 𝐱 𝐢 in layer m-1 and the weights that link them together 𝐰 𝐢𝐣 using the equations for z and y showed above.

Figure 45. Example of image filtering through a convolutional layer. An image patch from a CT slice, shown in yellow, is here analyzed. Each unit, marked as dark dots, activates themselves if the feature searched (determined by its weights) was detected. In this case the receptive field is 3 by 3 pixels. Note that recalling the replication of weights, all units in the same feature map search for the same feature (as they share the same weights).

Pooling layer

This layer is in charge of pooling partitions of its input into a set of non-overlapping rectangles and, for each such sub-region, outputs the maximum or the average value (see Figure 46 for an illustrated example). This pooling is a form of non-linear down-sampling and it is useful mainly for two reasons: by eliminating non-maximal values (or averaging depending the type of pooling done), it reduces computation for upper layers; and it provides a form of translation invariance (reducing the resolution of the feature map and reducing the sensitivity of the output to shifts and distortions).

Imagine cascading a pooling layer with a convolutional layer. There are 8 directions in which one can translate the input image by a single pixel. If pooling is done over a 2x2 region, 3 out of these 8 possible configurations will produce exactly the same output at the convolutional layer. For pooling over a 3x3 window, this jumps to 5/8. Since it provides additional robustness to position, pooling is a way of reducing the dimensionality of intermediate representations.

Fully-connected layer

Normally, last layers of a CNN are globally or fully connected. These layers (also called dense layers) act as the classifier itself. They look for combination of local features extracted by the early layers (which explains why do we need a lot of parameters here) to make the final decision. The activation of each output unit o m is calculated as follows (and represented in Figure 47):

o m = RELU(b m + ∑ x i . w im ) i , [ 2-40 ]
where x i is each of the units in the previous layer, w im the weight between the output unit and that x i , and b m the bias term. Then, again as in the convolutional layer, a non-linear activation function is applied. Note that even the equation to compute o m is similar to the case of 𝑦 (eq. ), in this case there are no constraints to the weights: all x i are considered to compute unit activations of the output layer.

SoftMax output

Optionally, this output layer is set after the fully connected layers in the network. It forces outputs to sum to one, so they can represent probability distributions across discrete mutually exclusive alternatives. Output 𝑦 𝑖 of each unit in this layer is obtained by (see Figure 48 for details):

𝑦 𝑖 = 𝑒 𝑧 𝑖 ∑ 𝑒 𝑧 𝑗 𝐾 𝑗=1 [ 2-41 ]
Figure 48. Group of units of a softmax layer from a neural network. Where, 𝒊: softmax output unit, 𝐲 𝐢 : output of unit 𝒊, 𝐳 𝐢 : total input received by 𝒊 unit (accumulated from the layer below), 𝑲: the total number of units in this layer. Note that the output 𝐲 𝐢 represents the probability of the right answer and its value goes from 0 to 1, and depends on all the inputs.

Neural network learning procedure

Learning in the context of artificial neural network consists of adapting the set of weights W to perform a given task (recognizing objects, speech, etc.). The procedure consists of going through the input patterns µ (where µ represents the index of one possible pattern in the training set) and feed their elements 𝑥 𝑛 𝜇 into the network to compute the predicted output 𝑂 𝜇 and compare it with the target 𝑌 𝜇 (see Figure 49): The aim of learning is to minimize the defined error between 𝑂 and 𝑌 summed over all training cases (a measure of that error simply could be the squared difference between 𝑂 and 𝑌, the actual output, and target output). Note that 𝑂 and 𝑌 are vectors that represent classes. An iterative approach is used to solve this (to generalize to more complex systems with several layers). Convergence is to get closer to a good set of W by properly modifying it (the idea is to get the actual output value closer to the target value). There may be not a perfect answer (no set of W to get the desired answer), but there is a set of 𝑊 that gets the best approximation (minimize the error measure) with the given training cases. We start with random guesses and then we adjust to give a better fit with a learning rule:

𝑂 𝜇 = 𝑌 𝜇 . [ 2-42 ]
𝑤 𝑖𝑗 𝑛𝑒𝑤 = 𝑤 𝑖𝑗 𝑜𝑙𝑑 + ∆𝑤 𝑖𝑗 . [ 2-43 ]

Defining the error measure

We are interested in this section in finding a learning rule that allows us to determine the set of weights by successive improvements from an arbitrary starting point. That is why we shall start by defining the error measure or loss/cost function. In this section we consider a simple perceptron without hidden layers just for demonstration. In the next section we generalize the idea to multilayer networks.

A loss function 𝐸 measures the discrepancy between the correct or desired output 𝑌 𝑖 𝜇 for pattern 𝜇 and the output produced by the system 𝑂 𝑖 𝜇 . As we discussed in previous section, learning

𝑥 𝑛 𝑥 𝑛 𝑂 𝜇 ? 𝑂 𝑖 𝜇
consists in finding the value of W that minimizes the average loss function 𝐸 (average of the errors over a set of labeled examples). Ideally, we would aim to find a good minimum of this global loss function with respect to all the parameters in the system. If the loss function measuring the performance can be made differentiable with respect to the system's tunable parameters 𝑊, then we can find a local minimum of E using gradient-based learning [START_REF] Lecun | A tutorial on energy-based learning[END_REF]. A typical error measure is defined as the squared residual summed over all training cases:

𝐸[𝑊] = 1 2 ∑ (𝑌 𝑖 𝜇 -𝑂 𝑖 𝜇 ) 2 𝑖𝜇
. As the output 𝑂 𝑖 𝜇 of the system depends on the set 𝑊, 𝐸 also does, and therefore, 𝐸 will be smaller the better our 𝑊 is. 𝐸 is normally positive but goes to zero as we approach the solution. Note that this cost function depends only on the weights and the problem patterns. The evolution is only affecting the weights, not the activation of the units themselves. Given the error measure 𝐸[𝑊] we can improve on a set of 𝑊 by sliding downhill on the surface it defines in W space. The learning rule (batch delta rule) for linear neurons (where the output of each neuron is a linear combination of the inputs 𝑂 𝑖 = 𝑅𝐸𝐿𝑈(∑ 𝑤 𝑖𝑘 𝑥 𝑘 )

𝑘

) results in:

∆𝑤 𝑖 = -𝜂 𝜕𝐸 𝜕𝑤 𝑖 = -𝜂 𝜕𝐸 𝜕𝑂 𝜕𝑂 𝜕𝑤 𝑖 = 𝜂 ∑ 𝑥 𝑖 𝜇 (𝑌 𝑖 𝜇 -𝑂 𝑖 𝜇 ) 𝑖𝜇 , [ 2-45 ]
where η is the learning rate parameter and the sign in the first term is introduced to get the error go down in each update. This approach requires continuous valued units with differentiable activation functions. The error surface lies in a space with axes for each weight and another for the error (see Figure 50). For linear units and squared error, 𝐸 is a quadratic bowl (vertical cross sections are parabolas and horizontal cross sections are ellipses). For multilayer non linear nets, the error surface will be much more complicated.

The gradient descent rule produces changes in the weight vectors only in the direction of the pattern vectors x i μ . Thus, any component of the weights orthogonal to the patterns is left unchanged by the learning. Within the pattern subspace, the gradient descent rule necessarily decreases the error if η is small enough, because it takes us in the downhill gradient direction. Thus, with enough iterations we approach to the bottom of the valley arbitrarily closely, from any starting point. Note that any point at the bottom of the valley solves the original problem exactly [START_REF] Hertz | Introduction to the Theory of Neural Computation[END_REF].

(a) (b) Figure 50. (a) Illustration of the shape of the error surface for linear units [START_REF] Hertz | Introduction to the Theory of Neural Computation[END_REF], and (b) the same but only for two input weights w1 and w2. For gradient descent learning, the idea is to do steepest descent on this error surface starting from a random set of W.

The quadratic function is not the only possibility. The cross-entropy loss is useful to train a system for a classification problem with 𝑛 classes, especially if the training data set is actually probabilistic or fuzzy as, for instance, in the association of symptoms with causes in medical diagnosis. The idea of this loss function is that we want to maximize the log probability (given by the SoftMax layer, see eq. ) of getting the answer right. We call this negative log likelihood criterion 𝐿 and it has the following form:

𝐿 = -∑ 𝑡 𝑖 𝑁𝑐𝑙𝑎𝑠𝑠𝑒𝑠 𝑖 log 𝑦 𝑖 [ 2-46 ]
where 𝑡 𝑖 is the target indexed by 𝑖, and 𝑦 𝑖 (equal to a 𝑓(𝑧 𝑖 ) as recalled from eq. [ 2.41 ]) is the activation of the 𝑖-th output unit (the probability of the right answer obtained), and 𝑧 𝑖 the total input accumulated from the layer below for the 𝑖-th unit as explained in section 2.5.2.4. The input (log 𝑦 𝑖 )

given is expected to contain log-probabilities of each class. Obtaining log-probabilities in a neural network is easily achieved by adding a LogSoftMax layer in the last layer of the neural network. Note that the choice of the loss function affects the computation of the gradients.

𝐿 has a very large gradient when target value is 1 and output close to 0. This means that 𝐿 has a very steep derivative when the answer is very wrong. Based on eq. [ 2-41 ], we can compute:

𝜕𝐿 𝜕𝑧 𝑖 = ∑ 𝜕𝐿 𝜕𝑦 𝑗 𝑁𝑐𝑙𝑎𝑠𝑠𝑒𝑠 𝑗 𝜕𝑦 𝑗 𝜕𝑧 𝑖 = 𝜕𝐿 𝜕𝑦 𝑖 𝜕𝑦 𝑖 𝜕𝑧 𝑖 + ∑ 𝜕𝐿 𝜕𝑦 𝑗 𝑗≠𝑖 𝜕𝑦 𝑗 𝜕𝑧 𝑖 [ 2-47 ]
where, again, z i = ∑ x j j w ij is the accumulated input of the 𝑖-th unit in the output layer. Keeping in mind that, from eq. we have: The loss function can be minimized by estimating the impact of small variations of the parameter values on the loss function. This is measured by the gradient of the loss function with respect to the parameters. In the next section we develop this approach to determine the learning rule for multilayer ANN: backpropagation.

𝜕𝑦 𝑖 𝜕𝑧 𝑘 = { 𝑒 𝑧 𝑖 ∑ 𝑒 𝑧 𝑗 𝐾𝑐𝑙𝑎𝑠𝑠𝑒𝑠 𝑗=1 -( 𝑒 𝑧 𝑖 ∑ 𝑒 𝑧 𝑗 𝐾𝑐𝑙𝑎𝑠𝑠𝑒𝑠 𝑗=1 ) 2 if i = k - 𝑒 𝑧 𝑖 𝑒 𝑧 𝑘 (∑ 𝑒 𝑧 𝑗 𝐾𝑐𝑙𝑎𝑠𝑠𝑒𝑠 𝑗=1 ) 2 if i ≠ 𝑘 = { 𝑦 𝑖 (1 -𝑦 𝑖 ) if i = k -𝑦 𝑖 𝑦 𝑘 if i ≠ 𝑘 , [ 2

Backpropagation, learning hidden units connections: learning features

The aim is to find good features without requiring insights into the task, that is, automate loop of designing features for a particular task. Backpropagation is a method to change weights in the feed forward network that are used later to obtain the desired feature representation. We can know how fast the error change as we change a hidden activity, so we could use error derivatives with respect to hidden activities. Each hidden activity can affect many output units (can have many separate effects on the error which must be combined). With backpropagation we compute all delta errors for all hidden units at the same time, starting from 𝑑𝐸/𝑑𝑦 𝑗 to 𝑑𝐸/𝑑𝑊 𝑖𝑗 and 𝑑𝐸/𝑑𝑊 𝑘𝑖 as illustrated in Figure 51.

The procedure starts by the initialization of the weights with small random values. Then we choose an input pattern 𝜇 and apply it to the input layer. We then propagate the signal forward through the network using the respective sum and activation function for each unit until the final outputs have all been calculated. Then, we compute the deltas for the output layers by comparing the actual output with the desired ones for the given pattern 𝜇. Next, we compute the deltas for the preceding layers by propagating the errors backward in the sense indicated in Figure 51(b) until a delta has been calculated for every unit. These deltas will be used for updating all connections. Finally, this procedure has to be repeated for every pattern.

We have already demonstrated how to update W ij connections. To find ∆W ki , the inputhidden connections update, we must differentiate with respect to W ki . Using the chain rule we obtain:

∆𝑤 𝑘𝑖 = -𝜂 𝜕𝐸 𝜕𝑤 𝑘𝑖 = -𝜂 ∑ 𝜕𝐸 𝜕𝑦 𝑖 𝜇 𝜕𝑦 𝑖 𝜇 𝜕𝑤 𝑘𝑖 𝜇 , [ 2-50 ]
where, for example, the term

𝜕𝑦 𝑖 𝜇 𝜕𝑤 𝑘𝑖
can be calculated recalling the dependency between 𝑦 𝑖 and 𝑤 𝑘𝑖 :

𝑦 i 𝜇 = RELU( ∑ 𝑤 𝑘𝑖 𝑥 𝑘 k ) . [ 2-51 ] (a) (b)
Figure 51. A section of a multilayer neural network with only one hidden layer is shown in both figures. Where 𝒋: output unit, 𝒊: hidden unit, 𝒌: input unit, 𝒀𝒊: output of unit 𝒊, 𝒁𝒋: total input received by 𝒋 unit, 𝑾𝒊𝒋: connection between hidden units and output, and 𝑾𝒌𝒊: connections between input units and hidden units. Back-propagate error signal to get derivatives for learning is shown in b).

Ensuring a good generalization: overcoming overfitting

We need to ensure that the 𝑊 found will work well for cases not used for training, this is, generalization issue needs to be considered to convert the backpropagation into a learning procedure (our model will be convincing if it fits a lot of data well). The main problem that appears is the overfitting, the downside of using such powerful models as deep (many hidden layers) ANN. This problem can be illustrated with a curve-fitting perspective as shown in Figure 52 where it can be seen that a curve fitted with too many parameters follows all the small details or noise but is very poor for interpolation and extrapolation. Of course, the same is true for neural networks where too many weights in a network could result in a poor generalization. Training data contains information about the regularities in the mapping input-output, but also has noise: target values may be unreliable, and sampling error can occur (accidental regularities for particular cases that are not real and true regularities that we want to generalize). Several techniques are useful for avoiding overfitting the network:

• Artificial increase data, change data: the idea is that we do not say exactly what to do, just add more cases to be generalized and to come up with the right conclusion, which can be considered as a process of education;

• Weight decay (keep 𝑊 small, makes the model simple) as regularization;

• Early stopping [START_REF] Lecun | A tutorial on energy-based learning[END_REF];

• Model averaging: train lot of different networks and average them together;

• Drop out on fully connected layers: randomly remove hidden units when training the network to regularize weights [START_REF] Srivastava | Dropout: a simple way to prevent neural networks from overfitting[END_REF].

As a final comment, we could see this search of a 𝑊 set as a search for a concept, a set of semantic features (vector of feature activities) that explain similarities between other concepts. Many neurons are used for representing each of the concepts leading to a distributed representation.

CNN model evaluation

Deep Learning models involve developing a software program that learns to solve a task instead of a program that solves the task. Traditional debugging step by step is no longer possible. When the overall system is not working as expected, there are several processing stages that should be checked:

▪ Process input? ▪ Debug model? ▪ Post-processing? (improve outcome)

This section discusses about tools to test the ability of our model (learnt features) to solve the learnt task (model generalization). This is, to evaluate whether a claim or the result of a calculation can possibly be true. The goal is to rule out certain classes of obviously false results and to perform a run-through of the system functionality to ensure that each part of the methodology works roughly as expected.

Monitoring the Deep Learning model implies several aspects presented in the following sections: visualization of features (this include the raw data input), training procedure, and results.

Feature visualization

It is essential to visualize data graphically to understand data properties (looking at the distributions of certain variables or potential correlations between variables), clean the data (detection of outliers as in Figure 53), interpret the results, debug and monitor neural network learning. To understand the need of data visualization and its benefits on avoiding overfitting phenomenon, note that the training error of a model can be reduced by increasing the model complexity. Although the training error for such a complex system may be zero, the test error can be large because the model may accidently fit some of the noise points in the training data (degrading the performance due to a poor generalization to the test examples). Overfitting and underfitting are two "pathologies" that are related to the model complexity. The causes are:

• Presence of noise: a given proportion of the training set is mislabeled. The X attribute test condition in model is spurious because it fits the mislabeled training records, which leads to the misclassification of samples in the test set.

• Lack of representative samples: models that make their classification decisions based on a small number of training examples are also susceptible to overfit. This occurs because of a lack of representative samples in the training data, together with a learning algorithm that continues to refine its model even when few training samples are available. Our case study datasets have a large number of variables. They have a high number of dimensions along which the data is distributed. Visually exploring the data can then become challenging. This can be achieved using techniques known as dimensionality reduction as we have seen in §2.1.3. We recall that the images are all essentially 512-by-512 pixel images and therefore have a total of 262144 dimensions, each holding the value of one specific pixel. It will be seen later in the methodology that in our case, the input image to classify could be a 32-by-32 pixels patch from a CT slice. What we can do is reduce the number of dimensions drastically whilst trying to retain as much of the variation in the information as possible.

We presented t-SNE technique to visualize high-dimensional data by giving each datapoint a location in a two or three-dimensional map. t-SNE is better than existing techniques at creating a single map that reveals structure at many different scales. This is particularly important for highdimensional data that lie on several different, but related, low-dimensional manifolds, such as images of objects from multiple classes seen from multiple viewpoints. t-SNE is capable of capturing much of the local structure of the high-dimensional data very well, while also revealing global structure such as the presence of clusters at several scales. Now that we have (let say) the two resulting dimensions, we can again visualize them by creating a scatter plot of the two dimensions and coloring each sample by its respective label. Figure 54 illustrates the performance of t-SNE dimensionality reduction technique on MNIST data set [START_REF] Maaten | Visualizing data using t-SNE[END_REF] from an image processing domain. In this work, we will use this tool for visualization purposes, to actually see in two dimensions the distribution (or point clouds) of raw data samples, and, on the other hand, samples transformed after CNN convolutional layers (just before the fully connected layers), to check if the transformation performed by the CNN helps disentangling the differences between classes, i.e., if the classes clusters are less correlated between each other.

Model training visualization

The errors committed by a classification model are generally divided into two types: training errors and generalization errors. Training error is the number of misclassification errors committed on training records, whereas generalization error is the expected error of the model on previously unseen records. We insist on the fact that a good classification model must not only fit the training data well, it must also accurately classify records it has never seen before. This is important because a model that fits the training data too well can have a poorer generalization error than a model with a higher training error, as we have demonstrated in previous section (overfitting issue).

Evaluation of the performance of a classification model is based on the counts of records (from any subset: training, validation, test) correctly and incorrectly predicted by the model. These counts can be tabulated in a table known as a confusion matrix 𝐶 (example shown in Figure 55). Each entry 𝐶 𝑖𝑗 in this table denotes the number of records from class 𝑖 (target, or ground truth label) predicted to be of class 𝑗 (prediction, or output of the model). Although a confusion matrix provides the information needed to determine how well a classification model performs, summarizing this information with a single number would make it more convenient to compare the performance of different models (score indexes discussed in section 2.5.4.3). There are two main useful quantities to be monitored during training of a neural network: loss function and accuracy. These values are normally represented with plots where the 𝑦-axis is one of the measures and the 𝑥-axis represents units of epochs, which measure how many times every example has been seen during training in expectation (one epoch means that every example has been seen once). It is preferable to track epochs rather than iterations since the number of iterations depends on the arbitrary setting of batch size.

These two plots get intuitions about the training process and the different hyperparameter settings and how they should be changed to get to a more efficient learning. The first quantity that is useful to track during training is the validation/training loss, as it is evaluated on the individual batches during the forward pass. The second important quantity to track while training a classifier is the validation/training accuracy. This plot can provide valuable insights into the amount of overfitting in your model. In Figure 56 there is an example showing the accuracy and the loss over time.

The gap between the training and validation accuracy indicates the amount of overfitting. In practice it is probably necessary to increase regularization (stronger L2 weight penalty, more dropout, etc.) or collect more data. The other possible case is when the validation accuracy tracks the training accuracy fairly well as in Figure 56. In this case, and if the final accuracy is not high enough, could signify that the model capacity is not high enough (making the model larger by increasing the number of parameters could make the model infer more complex functions from input data). Finally, for image processing tasks, it can be helpful to plot the first-layer kernels (convolutional filters or learnt features detectors) visually (Figure 57).

True label Prediction label

Result visualization

The performance of a model can be expressed in terms of its accuracy and error rate. We evaluate the following parameters for each of the proposed methods on a sparse axial image subset of the test database: sensitivity (True Positive Rate, TPR), specificity (True Negative Rate, TNR), accuracy (ACC) and precision (PREC). These indicators are calculated as follows:

𝑇𝑃𝑅 = [𝑇𝑃] [𝑇𝑃]+[𝐹𝑁] , 𝑇𝑁𝑅 = [𝑇𝑁] [𝑇𝑁]+[𝐹𝑃] , 𝐴𝐶𝐶 = [𝑇𝑃]+[𝑇𝑁] 𝑁 , 𝑃𝑅𝐸𝐶 = [𝑇𝑃] [𝑇𝑃] +[𝐹𝑃] , [ 2-52 ]
where

CNN, an example

This section describes a concrete example of a CNN, the architecture of LeNet-5, the Convolutional Neural Network used in [START_REF] Lecunn | Gradient-based learning applied to document recognition[END_REF] for on-line handwriting recognition. As we can observe in Figure 58, it comprises 7 layers (not counting the input), all of which contain trainable parameters. The network input is a 32x32 array that receives a pixel image of a particular handwritten digit scaled to a standard size. This input feeds forward through hidden layers to the 10 output units. Units that share the same set of weights are organized in planes. In the first hidden layer C1 of LeNet-5 they are organized in six planes, each of which generates a feature map. A unit in a feature map has 25 inputs connected to a 5X5 area in the input: the receptive field of the unit. Each unit has therefore 25 trainable coefficients plus a trainable bias (this is the descriptor to be learnt). As stated earlier, all the units in a feature map share the same set of 25 weights and the same bias, so they detect the same feature at all possible locations on the input. The other feature maps in the layer use different sets of weights and biases, thereby extracting different types of local features. A sequential implementation of a feature map would scan the input image with a single unit that has a local receptive field and store the states of this unit at corresponding locations in the feature map. This operation is equivalent to a convolution, followed by an additive bias and squashing function, hence the name convolutional network. The kernel of the convolution is the set of connection weights used by the units in the feature map. Subsampling layers S are layers with equal number of feature maps and half size than the convolutional layer before them. Each unit in each subsampled feature map is connected to a 2x2 neighborhood in the corresponding feature map in previous layer. The four inputs to a unit in S are added, then multiplied by a trainable coefficient, and added to a trainable bias. The result is passed through a sigmoid function. The 2x2 receptive fields are nonoverlapping, therefore feature maps in S have half the number of rows and columns as feature maps in C Layer.

Since all the weights are learned with back propagation, convolutional networks can be seen as synthesizing their own feature extractor. The weight sharing technique has the interesting effect of reducing the number of free parameters, thereby reducing the capacity of the network and reducing the gap between test error and training error. The network in Figure 58 contains 345308 connections, but only 60000 trainable free parameters because of the weight sharing.

Fully convolutional networks (FCN)

At this point is clear that convolutional networks are a powerful visual models based on hierarchies of features that drive advances in image recognition from coarse to fine inference to make a prediction at a pixel level. The typical use of convolutional networks is on classification tasks, where the output to an image is a single class label. However, in many visual tasks, especially in biomedical image processing, the desired output should include localization, i.e., a class label is supposed to be assigned to each pixel. In this section, we present insights of the so called fully convolutional network useful for semantic segmentation [START_REF] Long | Fully convolutional networks for semantic segmentation[END_REF]. It is an end-to-end (image to image) system that takes input of arbitrary size and produces correspondingly-sized output with both efficient learning and inference.

If we recall from previous section, CNNs basic components are convolution, pooling, and activation functions which operate on local input regions and depend only on relative spatial coordinates. It is the fully connected layers that put conditions on the input size (that shall remain fixed). While a general deep net computes a general nonlinear function, a net with only pooling and convolutional layers (not fully connected ones) computes a nonlinear filter, which is called FCN [START_REF] Long | Fully convolutional networks for semantic segmentation[END_REF]. In this way, a FCN predicts dense outputs from arbitrary-sized inputs. Note that both learning and inference are performed on the whole image at a time by dense feedforward computation and backpropagation. As it can be observed in Figure 59, the main idea is to supplement a usual contracting network by successive layers with upsampling operators layers (deconvolution layers), enabling pixelwise prediction and learning. Adding layers and a spatial loss produces an efficient machine for end-to-end dense learning [START_REF] Long | Fully convolutional networks for semantic segmentation[END_REF].

One important modification of this FCN was introduced in [START_REF] Ronneberger | U-net: Convolutional networks for biomedical image segmentation[END_REF], where the upsampling part now has a large number of feature channels, which allow the network to propagate context information to higher resolution layers. As a consequence, the "expansive path" is more or less symmetric to the "contracting path", and yields a u-shaped architecture, giving the name U-Net to this type of network. The architecture (illustrated in Figure 60) consists of a contracting path (left side) to capture context and a symmetric expanding path (right side) that enables precise localization. The contracting path follows the typical architecture of a convolutional network. It consists of the repeated application of two 3x3 convolutions (unpadded convolutions), each followed by a rectified linear unit (ReLU) and a 2x2 max pooling operation with stride 2 for downsampling. At each downsampling step we double the number of feature channels. Every step in the expansive path consists of an upsampling of the feature map followed by a 2x2 convolution ("up-convolution") that halves the number of feature channels, a concatenation with the correspondingly cropped feature map from the contracting path, and two 3x3 convolutions, each followed by a ReLU. The cropping is necessary due to the loss of border pixels in every convolution. At the final layer a 1x1 convolution is used to map each 64-component feature vector to the desired number of classes. In total, the network has 23 convolutional layers. The energy function is computed by a pixel-wise softmax over the final feature map combined with the cross entropy loss function [START_REF] Ronneberger | U-net: Convolutional networks for biomedical image segmentation[END_REF]. As in [START_REF] Long | Fully convolutional networks for semantic segmentation[END_REF], the network does not have any fully connected layers and only uses the valid part of each convolution, i.e., the segmentation map only contains the pixels for which the full context is available in the input image. This strategy allows the seamless segmentation of arbitrarily large images by an overlap-tile strategy (tiling strategy is necessary to apply the network to large images, since otherwise the resolution would be limited by the GPU memory). The U-Net architecture achieves very good performance on very different biomedical segmentation applications [57] [58].

Contraction Expansion

Style representation using CNNs

It has been demonstrated the effectiveness of CNNs applied to object recognition. We believe that this success is due to a learnt representation of the image that makes object information increasingly explicit through the hierarchical processing. Therefore, along the processing hierarchy of the network (going forward deeper in the net), the input image is transformed into representations that increasingly search for the actual content of the image (more abstract one, or high level content in terms of objects and their arrangement) compared to its detailed pixel values. These feature responses of the network are referred to as the content representation. In [START_REF] Gatys | A neural algorithm of artistic style[END_REF], an artificial system based on a Deep Convolutional Neural Network was presented to create artistic images of high perceptual quality. The authors showed that CNNs can "understand" not only the content, but also the style (not the global information, but the texture) of a given input image. They showed this by synthesizing images that mix simultaneously the content and style representation from two different source images. To obtain a style representation of an input image, the method use a feature space built with the correlations between the different filter responses of the feature maps. A visual example is shown in the next figure.

Style reconstruction

Content reconstruction

Input Feature maps (1) (2) Feature maps correlation

Feature maps content

The results presented in [START_REF] Gatys | A neural algorithm of artistic style[END_REF] were generated on the basis of the VGG-Network. Details of this architecture can be viewed in [START_REF] Simonyan | Very deep convolutional networks for large-scale image recognition[END_REF]. For this particular application, there is no need of fully connected layers. To generate a texture from a given source image, features of different sizes are extracted from it. Next, the correlation of the feature responses is computed to obtain a stationary description of the source image. This description is stored in the Gram matrix 𝐺 𝑖𝑗 𝑙 that is simply the inner product between the feature map 𝑖 and 𝑗 in the same layer 𝑙 and defined as follows:

𝐺 𝑖𝑗 𝑙 = 1 M ∑ F ik l . F jk l M k=1
, where F np l = 𝑛 𝑡ℎ 𝑓𝑒𝑎𝑡𝑢𝑟𝑒 𝑚𝑎𝑝 𝑣𝑎𝑙𝑢𝑒 𝑎𝑡 𝑝𝑜𝑠𝑖𝑡𝑖𝑜𝑛 𝑝, 𝑓𝑜𝑟 𝑙𝑎𝑦𝑒𝑟 𝑙, 𝑀 = ℎ𝑒𝑖𝑔ℎ𝑡 𝑡𝑖𝑚𝑒𝑠 𝑡ℎ𝑒 𝑤𝑖𝑑𝑡ℎ 𝑜𝑓 𝑡ℎ𝑒 𝑓𝑒𝑎𝑡𝑢𝑟𝑒 𝑚𝑎𝑝, 𝑙 = 𝑙-𝑡ℎ 𝑎𝑛𝑎𝑙𝑦𝑧𝑒𝑑 𝑙𝑎𝑦𝑒𝑟.

Finally, a new image with the same stationary description is obtained by performing gradient descent on a random image that has been initialized with white noise (to find another image that matches the correlation-matrix representation of the original one) [START_REF] Gatys | A neural algorithm of artistic style[END_REF]. Results of applying this procedure are shown in the next figure. This approach introduced a new parametric texture model based on the powerful feature spaces of CNNs (previously trained for object recognition). Constraining a spatial summary statistic over feature maps suffices to synthesize high-quality textures. However, in [START_REF] Ustyuzhaninov | Texture Synthesis Using Shallow Convolutional Networks with Random Filters[END_REF] it is demonstrated that such optimized multi-layer feature spaces are not imperative for texture modeling. Instead, much simpler shallow and convolutional networks can serve as the basis for novel texture synthesis algorithms. It is show in [START_REF] Ustyuzhaninov | Texture Synthesis Using Shallow Convolutional Networks with Random Filters[END_REF] They demonstrate that the feature space of random shallow CNNs can serve as a surprisingly good model of natural textures. From all the models presented in the paper, the most accurate in terms of results was a single-layer CNNs named Random-Multiscale. Figure 63 illustrates its architecture. It contains standard rectified linear units and convolutions with stride one, no bias and padding ((𝑓 -1)/2) where 𝑓 is the filter-size, to ensure that the spatial dimension of the output feature maps is the same as the input. It contains eight different filter sizes 𝑓𝑥𝑓𝑥3 with 𝑓 = 3; 5; 7; 11; 15; 23; 37; 55, and 128 feature maps each (1024 feature maps in total). Filters are drawn from a uniform distribution. In [START_REF] Ustyuzhaninov | Texture Synthesis Using Shallow Convolutional Networks with Random Filters[END_REF], Figure 64, it is shown the effectiveness of this architecture (similar with respect to the VGG method [START_REF] Gatys | A neural algorithm of artistic style[END_REF] but without any training needed). Note that the Gram matrix allows discriminating textures whereas this is not possible in the original image space (raw pixel values). 65. Samples synthesized from the two described models for three different reference textures (first column). The second column shows samples from the VGG model from [START_REF] Gatys | A neural algorithm of artistic style[END_REF], and the last column from the multi-scale model (with 1024 feature maps) [START_REF] Ustyuzhaninov | Texture Synthesis Using Shallow Convolutional Networks with Random Filters[END_REF].

Original VGG [START_REF] Gatys | A neural algorithm of artistic style[END_REF] Random multi-scale

Machine learning -based classification of lung diseases: proposed approach

Based on described state of the art concepts used for extracting structural and textural information, it is proposed in this section a complete new machine learning-based framework for lung diseases classification. For this purpose, this chapter starts by explaining the construction of a proprietary annotated ILD database in section 3.1, essential for tuning and testing our machine learning method. In section 3.2, it is presented the specific machine learning approach chosen, based on convolutional neural networks, at two different stages of the research: first, we introduce a mixed approach merging parametric and non-parametric classifiers (used as a benchmark to measure the CNNs capabilities and constraints); second, a fully CNN-based approach is proposed as the final classification system.

Building a representative database

As it was already pointed out, pattern recognition approaches remove the element of human observer variation so that no intra-observer variation is inherent in the use of such approaches. However, the algorithms rely on training by human beings, using an agreed definition of a specific pattern. In this section, a secondary objective of this work is described: building a representative database of ILDs cases to help developing and validating CAD systems. The selection of the histological criteria for patient inclusion, the associated clinical parameters, and the content are detailed in the following.

CT examination is based on the principle of thin sections allowing high longitudinal resolution. An acquisition in fine collimation on the entire thorax during a single apnea thus provides high resolution volumetric data. This allows accurate multiplanar reconstructions and the use of post-processing techniques, in particular minimum intensity projection (mIP) as radiological diagnosis/inclusion criteria. As it was briefly explained in §1.1.2, the first step in computed tomography analysis is the careful visualization of images on reconstructed axial sections, moving up and down interactively within the thoracic volume. This stage often confirms the diagnosis of pulmonary fibrosis, but sometimes it is insufficient to distinguish accurately some phenotypes (honeycombing and traction bronchiectasis) and to evaluate the extension of pathological areas (mainly ground glass).

This study is divided into several work steps:

• radiologist agreement on a suitable protocol;

• constitution of the database by the clinical investigator at the level of the hospital participating in the study; • centralized and multidisciplinary diagnosis by the diagnostic validation committee of the study (including a first radiologist); • reading of CT examinations (using mIP) by a senior radiologist (diagnostic validation of the study) to perform the annotations.

Patient selection criteria

In collaboration with the Service of Radiology of the Avicenne Hospital, Bobigny, France, 6 lung parenchyma disorders were selected (emphysema, fibrosis, honeycombing, AEF, ground glass, and traction bronchiectasis, in addition to normal patterns) to analyze 59 cases of these most frequent ILDs at the Avicenne Hospital. We describe in this section both the inclusion and exclusion criteria for the selection of cases for our study. Inclusion criteria:

• diagnosis of idiopathic interstitial lung disease, • ongoing diagnostic investigation (incident cases),

• compliance with CT protocol,

• patient agreement for study.

Criteria for non-inclusion:

• diagnosis of chronic diffuse interstitial lung disease (CDILD) other than idiopathic interstitial pneumopathies (IIP): secondary PID, sarcoidosis, special PID (langerhans cell histiocytosis, lymphangioleiomyomatosis, alveolar lipoproteinosis, chronic eosinophilic pneumonia),

• histologically confirmed diagnosis of IIP other than IPF: NSIP (non specific interstitial pneumonia), DIP (desquamative interstitial pneumonia), RB-ILD (respiratory bronchiolitis interstitial lung disease), AIP (acute interstitial pneumonia),

• specific diagnosis prior to the study,

• technical impossibility to respect the CT protocol and • lack of patient agreement for the study.

Patients with CDILD follow the usual course of management. Patients with IPF (according to CT or histologic criteria), and patients with unclassifiable idiopathic interstitial lung disease constitute our study population and will be included respectively in groups A, B and C.

Acquisition protocol for CT series

A preliminary meeting of the radiologists involved in the study was conducted to ensure the homogeneity of data acquisitions and reconstructions. The study is based on an evaluation of CT examinations performed routinely. MSCT exams are performed according to a high-resolution protocol on any device allowing at least 64 slices per rotation. Acquisitions of the images were performed on patients in supine position, at maximum inspiration, without injection of contrast medium. The thickness of sections is sub-millimetric with standard irradiation parameters (120 kV, > 1 mA.s / kg depending on the weight of the patient). The reconstruction interval is at least millimetric, overlapping (typically 50%, for example, 0.6 mm slice thickness reconstructed every 0.3 mm) performed with both a high and low spatial resolution algorithm (standard pulmonary and mediastinal convolution kernels). The field of view was adapted to the size of the patient, and the reconstruction matrix is at least 512x512 pixels (larger matrices are also available, namely 768x768 pixels).

Annotation procedure

The purpose of high-quality CT annotation is to establish the ground truth for lung tissue classification as well as to show examples of CT findings related to a studied disease for teaching purposes. Indeed, since the annotations are intended for computerized pattern recognition, the ROIs have to delineate pathologic patterns very precisely to avoid the introduction of noise in the training data. The possibility to visualize and delineate two dimensional ROIs in the entire MSCT volume and to set the window level used for displaying the 16-bit DICOM image series on a computer screen was required for annotation. These specifications were fulfilled by exploiting existing graphical software MedView -developed at ARTEMIS Department of Télécom SudParis (see Figure 66). CT exams were read by one expert radiologist in the field of PID, using anonymized data. The radiologist opens an entire DICOM series and then draws precise ROIs in any layer (at well-defined depths) of the CT volume in the axial view. A PGM file format saves the ROIs corresponding to each annotated slice of the volume as a mask (with pixel values of different indexes, representing each of the possible pathologies).

Note that the lung region will be later on segmented automatically, so the radiologist is concentrating in drawing only target ROIs, not also the precise lung contour. The final ground truth mask will be thus limited to the lung surface. Patterns used to describe lung parenchymal disorders are not standardized among radiological communities (referring to the ontological overlap existing between the classes). Series with blur caused by breathing or heart motion, or containing artifacts were also annotated with a different label to decide later if exclude or not those areas from training-testing any classification algorithm. Note that healthy tissue was delineated in the studied series to provide a wide range of the aspects of normal lung parenchyma. The average time for annotating one case was approximately of 30 minutes, necessary for the radiologist to interpret the MSCT image series and to draw annotations highlighting the important events in the series.

For each patient, the annotations are performed over six characteristic regions of the whole volume that are relevant in the investigation protocol during a follow-up of the patient. This consideration is important to correctly sample all the inter-class variations through different regions of the thorax. These regions are illustrated in the next figure: For a given matrix size, different axial resolution (mm/pixel) may appear, Figure 68. To tackle these different resolutions, we perform a pixel size normalization explained in detail later in the proposed approach. This is one of the mandatory preprocessing steps to be included to ensure similar robustness degree of the features with respect to the protocol variability discussed in the first chapter.

Summary of database content

In this section, the content of the database resulting from the research is quantitatively described and analyzed. A distinction is made for each of the two main stages of the database content through the study: the first one, at the beginning of the research, containing 12 cases; and the final one, constituted of 59 cases. All of them include sparse 2D ROI annotated axial images.

In addition to the proprietary database, annotated patterns from public database [START_REF] Depeursinge | Building a reference multimedia database for interstitial lung diseases[END_REF], a multimedia collection of cases with ILDs at the University Hospitals of Geneva, were also included during the first stage of the research (to increase the samples for each of the classes). This dataset contains high-resolution computed tomography image series with 2D annotated regions of pathological lung tissue along with clinical parameters from patients with pathologically proven diagnoses of ILDs. The library contains 128 patients affected with one of the 13 histological criteria of ILDs, 108 image series with annotated lung tissue patterns as well as a comprehensive set of 99 clinical parameters related to ILDs. As it was mentioned, the annotations of the lung image regions are successively carried out by one experienced radiologist. However, no assessment of the quality of the annotations by measuring inter-agreement measurement was carried out. To illustrate the 2D annotated ROIs used as ground truth for training, validation, or test, some examples are presented in Figure 70. Annotations for only one patient is shown in Figure 71 to illustrate changes in ROIs as we move forward from top to bottom over the already mentioned six characteristic regions of the whole volume relevant during a follow-up of the patient.

(a) (b) (c) (d) (e) (f)
Figure 71 (cont. on next page). 11 CT axial images spaced in the longitudinal direction throughout the lung including the six regions described before, for a given patient in the database.

(g) (h) (i) (j) (k)
Figure 71 (cont.). 11 CT axial images spaced in the longitudinal direction throughout the lung (red: fibrosis, blue: emphysema, light blue: traction bronchiectasis, white: ground glass, violet: honeycombing, brown: AEF, green: normal).

Statistical analysis and visualization of raw data

Before addressing the solution to our classification problem, we pointed out the need of visualizing the raw data to detect outliers and to make first inferences about its distribution. In this section, such analysis is illustrated through two perspectives: an intuitive one, related to histogram measurements, and a second one, using dimensionality reduction techniques, already defined in the state of the art section.

Densitometric and histogram analysis

In this section, some basic densitometric measurements of computed tomography images and their corresponding statistic parameters (kurtosis, skewness, mean lung density, standard deviation) are shown. It is a first simple representation of the raw data attempting to visualize the salient features or characteristics (intra-class variability, relation between the visual aspect and histogram vector, etc.) at different scale levels: patch, slice, volume.

Previous works [START_REF] Depeursinge | Near-affine-invariant texture learning for lung tissue analysis using isotropic wavelet frames[END_REF], [START_REF] Anthimopoulos | Lung Pattern Classification for Interstitial Lung Diseases Using a Deep Convolutional Neural Network[END_REF], [START_REF] Depeursinge | Comparative performance analysis of state-of-the-art classification algorithms applied to lung tissue categorization[END_REF], [START_REF] Van Tulder | Learning features for tissue classification with the classification restricted boltzmann machine[END_REF] have shown that the patch size of 32x32 pixels gets the best reliability of the central pixel context. As we aim to reach a pixel by pixel analysis taking into account the immediate contextual surrounding information, this patch scale is discussed in the following:

• Patch level (32x32 pixels):

Patch extraction (for those partially or totally falling inside the lung) can be performed for all the slices containing sparse 2D annotations (Figure 72 illustrates for one patch, but the procedure is done by scanning the full slice). Decisions about conditions imposed for patch extraction, as overlap between contiguous patches, proportion of the patch inside the annotated ROI, and proportion of the patch inside the lung, are discussed in the next sections. • Slice level:

For the slice densitometric evaluation, the histogram of distribution of the density of each slice within the lung was plotted as shown in Figure 77. Also, skewness, kurtosis, and the mean of that distribution (mean lung density, MLD) were measured. In addition, the percentage of the total volume of tissue that had a density between -250 Hounsfield units (HU) and -600 HU was recorded as the percent of high attenuation area (HAA%).

Two main aspects are derived from observing examples in Figure 77. The first one is that for different pathological conditions, different grey level distributions will be computed (as it has been already observed several times, also for a patch level analysis). The second one is the acknowledgement of the existence of a peak located at -1024 HU of the histogram distribution that concentrates a lot of pixels. This cannot be noticed directly from visualizing the image, but it definitely has an impact for any automated processing aimed, as it is shown by the histogram curve and its associated parameters (namely kurtosis). -1020

• Volume level:

Finally, the same analysis is performed for the whole volume (only pixels inside the lung are considered), with similar findings. The histogram peak at -1024 HU is due to the noise level induced by the convolution kernel and/or to data conversion (note the rescale intercept value for all datasets listed in Table 4 is -1024 HU, irrespective to manufacturer). As this type of noise is present all over the lung field (Figure 79), it risks impacting the classification outcome. Consequently, it will be suppressed in the preprocessing phase using 3D Gaussian smoothing with a kernel support of 5x5x5 pixels (radius = 2) and a standard deviation 𝜎 equal to half pixel size (inducing small impact on textural patterns as shown in Figure 79). It can be observed (the cluster marked in red in Figure 81) that some patches of the subpleural region will be included in the training database, which means patches that contain high density due to the presence of pixels outside the lung. This is a decision taken after analyzing all patch configurations while scanning a whole slice (see orange and blue squares representing each of the two possible configurations in Figure 82). We consider necessary to fully represent all the possibilities for each of the classes. Another way to tackle the border patches is adding the patch filling procedure (detailed in the next section), which will replace high densities outside the lung field. This procedure will remove the classification ambiguity for the subpleural patches (Figure 82, orange). However, for mediastinal regions, the presence of large vessels inside the lung mask will mimic the same configuration as subpleural patches (Figure 82, blue). Consequently, this patch configuration should be maintained in the training dataset. From Figure 80, it is clear the persistent overlap between the classes, and the impossibility to essay any cluster separation using only pixel intensity values. A transformation to a new space that allows disentangling the most important features for each class is thus needed.

As the idea of using histogram distribution is recurrent in the literature, we have also plotted the t-SNE of the same database patch histograms in the next figure. It appears that regarding only the distribution of gray level pixels and not the textural information will be insufficient to build up a robust feature vector. However, the next section essays to demonstrate this fact by combining histogram representation of the raw data with a SVM classifier.

Figure 83. Visualization of the t-SNE applied to the database patch histograms.

State of the art parametric approach

The overlap between the classes seems to persist using histograms as discriminating feature vectors, as it has been shown in previous section describing the raw data. However, searching for confirming the direction of the research towards complex intelligent systems, it is proposed an initial analysis using a parametric approach (based on a state of the art work [START_REF] Depeursinge | Near-affine-invariant texture learning for lung tissue analysis using isotropic wavelet frames[END_REF]). The goal is to study whether the histogram that represents the distribution of the pixel values of each patch is a sufficient feature vector for an accurate classification (since it appears as a reference technique, at least in the medical image processing field). To answer this, the original dataset, composed of 11058 (32x32 pixel) patches with values between -1024 HU and 1024 HU representing 4 classes (normal, emphysema, fibrosis and ground glass) was divided randomly in train (90%) and test (10%) corpus, generating 1000 permutations to train a SVM (implemented using scikit-learn [START_REF] Pedregosa | Scikit-learn: Machine learning in Python[END_REF]). The patches are distributed as: 4000 normal, 641 emphysema, 4000 fibrosis, 2417 ground glass. The pipeline of the classification system is illustrated in the next figure. Each sample is represented as histograms (22 bins) to replicate the conditions used in [START_REF] Depeursinge | Near-affine-invariant texture learning for lung tissue analysis using isotropic wavelet frames[END_REF]. The hyperparameters of the SVM model are fixed by doing a grid search, to the following values: optimal cost of the errors C = 10, RBF kernel, Gamma = 0.01, and class weights according to the proportion of each class with respect to the total number of patches.

The input is randomly divided, while keeping the proportion of the classes, in train and test. Three different sets are analyzed to test the SVM ability of categorization of the target samples: Dataset a. Samples were generated entirely from the original dataset. This case corresponds to the original classification problem.

Dataset b.

The set of patches were kept identical to the original for classes representing diseases, but they were artificially created for normal patches to minimize overlapping of normal and other classes. The 4000 patches of this last class are generated randomly according to a normal distribution with a mean µ= 800 HU and standard deviation σ = 50 HU.

Dataset c.

A set where all classes except emphysema were manipulated in the same way as in (b): normal (µ = 800 HU, σ = 50 HU), fibrosis (µ = 100 HU, σ = 25 HU) and ground glass (µ = 650 HU, σ = 100 HU).

The experiments of SVM-based classification obtained for the three datasets above are illustrated in terms of classification accuracy in Figure 85 and Table 5. It can be observed that, for each dataset, the permutation scores are similar: 45.5% average for the original dataset with a standard deviation of 0.8%; 72.2% for the dataset where the normal class was manipulated, with a deviation of 0.08%; and 99.7% for the dataset where all classes except emphysema were manipulated, with a deviation of 0.1%. The average score for dataset c suggests that the model works correctly: if we manipulate the data by making them easily separable, the SVM reaches an accuracy of almost 100%. Even with the dataset that is partially manipulated (experiment 2, normal class biased) the model manages to separate the data with an accuracy of more than 72%.

From the obtained results, it is seen that the patch histogram is not a sufficient measure for the classification task of pulmonary pathologies, suggesting the need of other features to tackle the overlap between the classes (mainly textural ones that conserve spatial patterns information). It has been studied here the reference technique from the state of the art [START_REF] Depeursinge | Near-affine-invariant texture learning for lung tissue analysis using isotropic wavelet frames[END_REF] as a comparative starting point. The spatial information, not contained in the histogram vector, needs to be exploited further as actual input to the classifier.

Developed approach for ILD detection and quantification

Now that the data has been analyzed and some insights were gained into the complex problem of automated ILD classification, we address two questions that arose while searching for the optimal solution. First, in §3.2.2: would CNN be able to improve the classification of an alternative method, using only a restricted training database? As it has been already stated, at the first stage of this study, the included a small number of cases, which implied a constraint with respect to CNNs training. Secondly, in §3.2.3: can CNN alone obtain better results and how should be the design of the network to achieve this objective? Before explaining the proposed approaches, we describe complementary tools to prepare the data for the classification system that are applied to both methods presented in this chapter.

Data Preprocessing

This section concerns DICOM image processing before feeding it to the automated classifier. It includes auxiliary tools like lung segmentation, and "cleaning" tools like patch filling procedure. These modules are detailed in the following.

Lung segmentation

As we saw from the description of patch extraction, samples on the lung border include high densities from the thorax cage leading to configurations which could mislead the classifier about the presence of genuine fibrosis entities. To avoid this, lung segmentation rises as an auxiliary but mandatory task to be solved allowing the selection of intrapulmonary grayscales only. We have developed a robust procedure [START_REF] Fetita | Robust lung identification in MSCT via controlled flooding and shape constraints: dealing with anatomical and pathological specificity[END_REF] to obtain the lung masks even in difficult morphological configurations.

The underlying idea consists of iteratively segmenting the right and left lungs (contrasted regions), by exploiting the concept of controlled relief flooding summarized in Figure 86, and mathematically described by the grayscale reconstruction by erosion morphological operator. The control of the "flooding" is ensured by the definition of a reference set playing the role of a "sink" and preventing the "filling" of marked "valleys" in the relief (Figure 86). If we consider a schematic representation of the lungs, they can be thought of as connected vessels by means of trachea and stem bronchi (Figure 87 (a)). The consequence is that, if no sink constraint is imposed, the thorax relief flooding will equally affect the lungs (Figure 87 (b)) and their separation could be difficult in general. Indeed, the level of the "liquid" extracted from each lung by image subtraction may be different and a thresholding may generate several connected components which cannot be linked back together by regularization because of the risk to connect the lungs (e.g. configuration of Figure 87 (e), (f)). Instead, if the trachea region is "plugged" and a "sink" created in each lung at a time, the lungs can be extracted and regularized separately (Figure 87 (g), (h)). This idea is implemented in our algorithm, which is summarized in Figure 88 and discussed in the following. All steps are performed in 3D if no otherwise mentioned. They are illustrated at a level of a single axial image (Figure 89) with the major intermediary results also shown in 3D using unshaded volume rendering (Figure 91). Ultimately, a second processing pipeline is developed to detect eventual peripheral pathological dense lung tissue which cannot be marked by the flooding procedure. followed by self-adaptive (Otsu) thresholding [START_REF] Wenqing | The Automatic thresholding of gray-level pictures via twodimensional otsu method [J][END_REF], 3D hole filling and largest component selection (Figure 91 (a)). The maximum and median cross-sectional radius of the trachea (axial plane), 𝑅 𝑚𝑎𝑥 , 𝑅 𝑚𝑒𝑑 , are then computed.

Next, we need to suppress the trachea link between lungs (or equivalently, filling it with high values in the original CT images), and we used a controlled region growing in this respect. Because of a large inter-subject variability, a set of constraints is imposed to prevent a penetration of the growing set into the lung field, which would affect the lung shape. In addition to imposing an intensity constraint (grayscale values within trachea are lower than inside the lung texture), a set of "rough lungs" is built to stop the propagation. For this, a longitudinal (𝑧-axis) linear closing of size 𝑅 𝑚𝑎𝑥 is performed to fill-in the vessels holes (Figure 89 (d)), followed by an opening with a spherical SE of radius 2𝑅 𝑚𝑎𝑥 for trachea suppression (Figure 89 (e)). Finally, the mediastinal space is reduced by a linear closing of size 2𝑅 𝑚𝑎𝑥 in the sagittal direction (Figure 89 (f), Figure 91 (b)). The trachea is removed by avoiding penetration in the "rough" lung set and the two largest 3D connected components are kept with holes filled (Figure 89 (g), Figure 91 (c)).

Note that this partial result would be sufficient if the contrast is similar in both lungs (otherwise, the "flooded spaces" are not similarly selected by the threshold in Figure 89 (c), which imposes the more complex procedure in the following). From Figure 89 (g), an erosion with a spherical SE of radius 𝑅 𝑚𝑎𝑥 is performed and the largest connected component per lung is kept as "sink" reference set (Figure 89 (h), Figure 91 (d)).

The original image relief can now be selectively flooded by alternatively setting the reference set 𝑌 in eq. as one of the lung "sink" set extracted (Figure 89 (h)). The components corresponding to each lung are obtained by adaptive thresholding, as summarized in Figure 89 (i)-(l).

A final regularization of each lung (Figure 89 (j),(l)) consisting of largest 3D component selection and spherical closing of size 𝑅 𝑚𝑒𝑑 provides the segmented pulmonary field, Figure 89 (m)-(o). Note that the proposed procedure succeeds to overcome the presence of dense pathological tissue in the lung if it is not located at lung periphery; otherwise the lung segmentation may be incomplete at the peripheral regions, as shown in Figure 90. For this reason, a second processing pipeline is developed to complete the segmentation, targeting the peripheral dense tissue in the lung.

Detection of peripheral pathological dense lung tissue

The underlying idea here is to exploit the rib cage shape to define the peripheral lung region and add the areas comprised between "enclosed" rib cage, mediastinum and the actual segmented lungs. The algorithm is briefly summarized and illustrated in Figure 92 for case of Figure 90.

Results

The proposed approach was tested on a 40 MSCT dataset including several lung pathologies. The contrasted lung regions were correctly segmented and separated in all cases. The detection of peripheral dense lung tissue also performed well on this database (96% sensitivity, 90% specificity, Figure 93 (a)-(c)), but it remains challenging if contrast agent is used, because of the bias introduced in detecting only the bone structures in the thorax: injected blood vessels inside the lungs are also extracted which will bias the thorax cage segmentation; the consequence is the detected lung area enclosed by the altered thorax will be smaller and some part of the peripheral dense tissue will be missed. Note also that, in few cases, some false positives of dense lung tissue occurred in the top anterior part and intercostal bottom posterior part of the lungs (Figure 93 (d), (e)), mainly due to the ribs curvature change in these areas, which was not addressed specifically in the algorithm. These two drawbacks are subject of future research, however they did not impact our study since all exams were performed without contrast agent injection and the eventual false positives can be easily removed with interactive selection.

Patch filling procedure

In contrast to studies in [START_REF] Li | Medical image classification with convolutional neural network[END_REF][START_REF] Li | Lung image patch classification with automatic feature learning[END_REF] which analyze the classification performance on a patchbased approach by mixing patients from all database, in this work we target a patient-based evaluation. This means performing the classification and analysis on the entire lung field from several axial images. To avoid misclassification of some points near to the lungs edge, a patch filling procedure is proposed and detailed in the following.

As it was mentioned in §3.1.5.2, the subpleural patches in the training set are left intact to ensure that this type of patch is integrated into the class definition. Thus, only the test patches receive a final pre-processing before feeding them to the CNN, such that those patches that are not entirely embedded in the lung mask are completed with the useful content. The filling procedure is performed symmetrically to each axis and recursively for the four quadrants of the patch as illustrated in Figure 94. In this way, we avoid misclassification of some points near to the lung edge where the patch encounters very dense regions that could result in artifacts for the classification task.

𝑥 𝑦

The importance of considering the patch-filling procedure is illustrated in Figure 95 where a lung axial slice was analyzed with the same CNN (discussed in §3.2.2.3), first using raw patches for training and testing and secondly applying the lung masking procedure. We notice that using raw data patches, the CNN overestimates the pathologic regions (fibrosis and ground glass) near the lung border irrespective to the fact that the border patch configurations were included in the training dataset as discussed in §3.1.5.2 (Figure 81).

Hybrid method

In this section, we explore the contribution of a CNN to improve the classification of an alternative parametric method. Considering the reduced annotated database in the first stage of this study, this initial experiment aims to identify the strengths and weakness of CNNs, and to estimate what can be expected from its performance applied to lung texture analysis.

Purpose

We had the opportunity to beneficiate from a previously developed CAD system at ARTEMIS Télécom SudParis exploiting a fully-3D multi-scale morphological decomposition of the lung relief [START_REF] Chien | Detection and classification of interstitial lung diseases and emphysema using a joint morphological-fuzzy approach[END_REF] (further on referred to as CLS in this work), which takes advantage of the three-dimensional data information and was able to show good performance in detecting diseased areas. However, two major drawbacks limited its utilization in clinical routine (Figure 96): first, the system was too sensitive to the presence of small patterns which locally could be interpreted as disease at a given analysis scale, thus providing several false alerts (noise); second, different types of lung patterns were interchanged in the case of complex disease (such as small emphysema sacs misled for fibrosis). To overcome these important drawbacks, we have attempted a combination of the existing CAD system with the classification outcome provided by a CNN (as this combination has been already demonstrated as a powerful approach [START_REF] Roth | Deep convolutional networks for pancreas segmentation in CT imaging[END_REF]), specifically tuned-up, to increase the specificity of the classification and the confidence to diagnosis [START_REF] Tarando | Increasing CAD system efficacy for lung texture analysis using a convolutional network[END_REF]. The advantage of using a deep learning approach would be a better regularization of the classification output (because of a deeper insight into a given pathological class over a large series of samples) where the previous CLS system is extra-sensitive due to the multi-scale response on patient-specific, localized patterns.

The hybrid system using both CLS parametric and CNN non-parametric methods for an integrated classification is shown in the figure below and described in the following.

Train DB (2D patches): Learning features Test DB 3D

Figure 97. Flowchart of the hybrid system proposed. The CNN is pre-trained with the 2D annotated database to set its weights properly for the classification task. Then, test images from the 3D database are fed into both methods. The outcomes are combined in a fusion module to generate the third result: a mixed classification [START_REF] Tarando | Increasing CAD system efficacy for lung texture analysis using a convolutional network[END_REF].

3D classification using parametric techniques

The parametric classifier (CLS) used in this approach exploits a three-dimensional multi-scale morphological decomposition of the lung relief and applies fuzzy-based decision rules to assign a class type to each image voxel at each analysis scale [START_REF] Chien | Detection and classification of interstitial lung diseases and emphysema using a joint morphological-fuzzy approach[END_REF]. The workflow of this approach is briefly described in Figure 98 (for a detailed description the reader should refer to [START_REF] Chien | Detection and classification of interstitial lung diseases and emphysema using a joint morphological-fuzzy approach[END_REF]). The final class of each voxel is decided by integrating over the scales. Note that this parametric approach does not use a training procedure, the decision rules being made using a priori knowledge on the disease textural patterns. This approach has several advantages, namely exploiting the full 3D information of the dataset, but presents also few limitations as discussed in previous section (Figure 96), which require improvements based on additional classification techniques. The training set consists of 32x32 pixels size half-overlapping patches. We consider those patches falling entirely inside the marked ROIs in the ground truth (see §3.1.4). Horizontal flipping and rotations are applied in order to artificially increase the number of samples and avoid over-fitting the neural network [START_REF] Krizhevsky | Imagenet classification with deep convolutional neural networks[END_REF]. In this way, the training dataset includes 5576 image patches: 1776 normal, 1648 fibrosis, 940 ground glass, 1212 emphysema. Similarly, 32x32 pixels patches were extracted from the test dataset as input for the CNN, centered to each pixel in the test image to which a CNN class will be assigned with a given probability. We normalized the lung image patches before feeding them to the network by subtracting the training set mean value from the native pixel values and dividing by the training set standard deviation.

The CNN confusion matrix for the training data, after the learning process, is shown in Figure 100, proving the good adjustment of the network weights to fit training dataset. As a final remark, we show the final 5x5 pixel filters (or features) automatically learned by the first layer of CNN in order to disentangle the patch texture differences between the possible classes (Figure 101). As we were looking for, these kernels are not hand-crafted in contrast to the CLS approach. Note that the input of the CNN are patches from the whole slice, centered on each pixel to be classified.

Merging classification results

The classification results of both parametric and non-parametric approaches are merged in a final decision by calculating a simple average between the probabilities for each class and then selecting the maximum value between all the classes. Result integration in an 8-connected neighborhood around each image pixel is applied to both methods, having the objective in the case of CLS of generating the probability map for a given slice, and in the case of the CNN of smoothing the rough output of the classification.

Results

A qualitative evaluation of the results obtained with the proposed approach is presented in Figure 102. We note that the output provided by the CNN has a "coarser" description of the pathological zones but a more regular spatial distribution of the textural classes. This information is exploited by the parametric classifier to cross-check and to regularize its output.

The 3D test database not being fully annotated, a complete quantitative evaluation of the classification result was not possible. We evaluate the following parameters for each of the three methods on a sparse axial image subset of the test database: sensitivity (True Positive Rate, TPR), specificity (True Negative Rate, TNR), Accuracy (ACC) and Precision (PREC). These indicators are calculated as detailed in section 2.5.4.3.

Overall, a semi-quantitative evaluation using sparse 2D patches showed a sharp 10% increase in true positives detection, no variation in TNR that remains the same, an increase of 6.5% in the precision of the mixed classification compared to the CLS outcome, and a slight increment of 0.4% for the accuracy parameter (these values are obtained on average between all the four classes). The detailed quantitative results are organized from Figure 103 to Figure 106 below. The normalized confusion matrix for each method computed for the whole test database is illustrated in Figure 107. The CNN-CLS system increases the prediction of the correct labels for emphysema and fibrosis with similar behavior for normal and ground glass. However, it is interesting to note (also from Figure 102 to Figure 106) that the CNN alone reaches high prediction for emphysema and fibrosis with a more coherent redistribution of misclassed labels. We note less efficacy in predicting normal tissue and more marked overlap between fibrosis and ground glass classes.

Discussion

This hybrid method showed the potential of Deep Learning techniques. Even with reduced training databases, they are still effective for refining and regularizing the output of a complementary classification approach for lung disease diagnosis. With this experiment, we acquired conscience of the causes that lead to a low accuracy in the automated classification:

• Overlap between classes not solved due to:

➢ Intrinsic problem definition

The highly visual correlated classes (mainly fibrosis and ground glass) remain hard to be separated which implies rethinking the architecture to adapt it to a texture recognition problem (LeNet, as well as any "standard" state of the art CNN, works well for object recognition tasks, but not necessarily for texture classification).

➢ Structural noise due to vascular tree

Normal structures as vessels and bronchi must be cleaned out from the image to be classified as they bias the ROIs with high density values leading to misclassifications. To avoid this kind of confusion, a pre-processing step shall be considered to attenuate the vascular tree presence.

➢ Lack of acquisition standardization

Lack of CT acquisition normalization between heterogeneous patients leads to a strong variability of normal and pathological tissue appearance in addition to their ontological definition.

• Insufficient database

It was verified the great impact of reduced database on CNNs generalization. This requires the extension of the training database and selection of more representative cases obtained with different acquisition characteristics. This is the reason of the efforts put on considerably enlarging the database.

With these issues to be tackled and considering that two classifiers are needed for the hybrid method (that results in a long and complex approach) the next step focuses on developing a CNN approach alone for classifying lung pathologies. The efforts were directed to improve the detection of the "normal" class and reducing the confusion between fibrosis and ground glass classes as noticed in §3.2.2.5.

Cascade of Texture-CNN

The second approach developed in this work focused on an independent CNN architecture designed to overcome the limitations encountered by the previous method. We present new solutions for the critical points discussed in §3.2.2.6, mainly, the introduction of the Locally Connected Filtering (LCF) to reduce the effect of vascular tree on biasing patches and a new CNN based framework to capture the texture properties of raw medical data.

Purpose

According to the performance analysis of the CNN in §3.2.2, we conclude that its architecture requires an updating in order to fit our texture classification problem. In [START_REF] Anthimopoulos | Lung Pattern Classification for Interstitial Lung Diseases Using a Deep Convolutional Neural Network[END_REF] a specific network architecture (further referred to as T-CNN in this work) was designed to capture the low-level textural features of the lung tissue showing high classification performance for a train-validation-test patch scheme. The special regard on adapting the network to a texture recognition problem was the motivation to include it in our workflow. However, a major drawback was observed testing this network (Figure 108): the system fails to learn a good representation to disentangle the differences between two highly correlated classes (namely fibrosis vs. ground glass, already proven to be hardly differentiable, section 3.1.5.1, [START_REF] Zavaletta | High resolution multidetector CT-aided tissue analysis and quantification of lung fibrosis[END_REF], as their histogram signatures are almost equal) which leads to interchanged lung patterns. Also, as it will be described in detail later in this section, the accuracy achieved for the other classes is insufficient when we analyze the result obtained for an entire slice which also shows mixed classification patterns (unsuitable for a quantitative CAD). On the other hand, it was also mentioned the great visual appearance overlap between the classes (normal, emphysema, fibrosis, and ground glass) due to the high anatomical variability interintra-classes and the lack of acquisition standardization (different acquisition techniques and radiation dose for different manufacturers, scanner calibration, scanner models, etc.). As illustrated in Figure 109, it often makes the ILD classification task an ontological problem which leads to errors in the categorization, i.e., mixing different types of lung patterns, overestimating the pathological areas. Finally, it is also mandatory to remove or attenuate texture ambiguity due to vascular network, preserving high densities related to fibrosis or ground glass patterns, and therefore, enhance contrast between normal tissue versus pathological one. Taking into account all these considerations, a new CNN-based CAD system is proposed in the next section.

Overall lung texture classification system

To address the issues discussed previously, a cascade of two classifiers is proposed. The approach consists of a combination of the existing CAD system T-CNN achieving a first classification of less correlated patterns, with a second classification provided by a texture matching method (also CNN-based) for those classes that are highly correlated (namely fibrosis and ground glass). The advantage of using this pipeline is a better differentiation between classes than the original system due to the multi-scale texture features analysis performed by the second classifier, which is only applied to the two poorly separable classes. As previously, because only 2D samples were available in the training database, the convolutional network was designed to perform on a two-dimensional basis at the level of axial images.

In addition, a pre-processing step of the input data is added, issued from locally connected filtering (LCF), to remove or attenuate texture ambiguity due to vascular network, while preserving high densities related to fibrosis or ground glass patterns.

The complete classification workflow is shown below and described in the following. The next sections describe each of the modules in detail.

Locally connected filtering

As we saw in the previous CNN analysis ( §3.2.2.5), the poorer CNN results in detecting normal tissue may be explained by the strong variability of its textural appearance. The confusion with the fibrosis tissue may also come from the presence of high densities corresponding to the vascular network which suggests the need of paying a special attention to the vessel tissue.

The idea implemented in the following consists of introducing a specific filtering able to attenuate the vascular lucencies while preserving the other high densities related to other (pathological) structures. Such filtering exploits the knowledge of vascular morphology (a highly branching tree structure of quasi-cylindrical components of various calibers) and achieves to provide at the same time a segmentation of the vascular network completing the lung mask information. The filtering principle relies on morphological connected filters applied locally and extended to a multiresolution approach to cover the large caliber spectrum of the vessels, as discussed in the following.

Locally connected filters (LCF) exploit grayscale topological connectivity of each point on the support of a function 𝑓 to filter out 𝑓 values weakly connected with a pre-defined subset Y 𝑠𝑢𝑝𝑝 𝑓 in the local environment of that point. The topological connectivity is verified by means of a morphological filter, namely the grayscale reconstruction by dilation Rf  (., Y), or the grayscale reconstruction by erosion Rf  (., Y). Since the targeted structures here are the blood vessels, exhibiting a positive contrast with respect to their environment, the former operator Rf  (., Y) will be used. The shape and size of the reference subset Y will determine the type of the researched connectivity. Reminding the definition of Rf  (., Y) as [START_REF] Vincent | Morphological grayscale reconstruction in image analysis: applications and efficient algorithms[END_REF]:

) ( ) (., ) ( g Y R f f  =   , where f H g g f   = ) ( ) ( ) (1  and )) ( ( ) ( ) ( ) ( ) ( g g n f f n f 1 1 - =    , [ 3-2 ]
with H denoting the unitary structuring element,  the morphological dilation operator, ˄ the infimum operator and

𝑔(𝑥) = { 𝑓(𝑥), 𝑥𝑌 -, 𝑒𝑙𝑠𝑒𝑤ℎ𝑒𝑟𝑒 ,
the LCF by dilation is defined as:

x supp f, LCFf,Y  (x, k) = Rf  (x, Y) | Y Nx (k), [ 3-3 ]
where Nx (k) denotes the k-size spatial neighborhood of x,

Nx (k)={ y supp f | d(x,y)  k}, with d a distance function. [ 3-4 ]
LCF thus provides in each point 𝑥 of the function support the connectivity strength between 𝑓(𝑥) and the positive-contrast structures in a region 𝑌 of the 𝑥 neighborhood. If Y is chosen as the set of points on the border of Nx (k), eq. becomes

x supp f, LCFf  (x, k) = Rf  (x, Nx (k)-Nx (k-1)). [ 3-5 ]
LCFf  (•,k) reconstructs locally the 𝑓 value from a k-distant neighborhood by morphological dilation. Its effect is to attenuate (or suppress) the 𝑓 values which are not "linked" with their kdistant neighborhood via a high-intensity path. On contrary, when such connection exists, the structures are preserved via the reconstruction operator. Figure 111 shows an example of LCFf  (•,k) filtering for two 𝑘 values.

LCFf  thus not only presents a denoising property similar to the median filter, but also preserves spatial structures selected by a local connectivity configuration. Figure 112 illustrates the effect of LCFf  on 2D grayscale and binary images for a spatial neighborhood 𝑘 = 5, demonstrating spatial structures preservation compared with a median filter of the same size. In order to remove vascular structures in the 3D image datasets, we take into account the linear shape of the vessels and their size, and set-up a filtering scheme combining 2D LCFf  applied across planes oriented orthogonal to different directions in space. We chose 9 spatial directions corresponding with the 18-connectivity (excluding symmetric orientations), Figure 113 (a). We define the vascular removal LCF of size k as: To remove vascular structures of different size, the VLCF is applied in a multiresolution scheme using 2 levels of decimation. The input at each resolution level is the filtered image from the previous level, enhanced by a multiplication with a small constant (𝛼 > 1) without exceeding the input (Figure 115). In addition, the vascular structures are detected at each level of resolution and combined prior to a final adaptive thresholding and filtering which select the highest confidence vessels. The vascular structure thus extracted is added to the lung mask obtained in §3.2.1.1 to reinforce the train-and-test process of the network as discussed later. Figure 116 shows the multiresolution filtering result together with the extracted vascular information for the image in Figure 114 (a).

x supp f, VLCFf  (x, k) = inf dC18 { LCFf  (x, kd)},

Classification using trained T-CNN

We use the architecture described in [START_REF] Anthimopoulos | Lung Pattern Classification for Interstitial Lung Diseases Using a Deep Convolutional Neural Network[END_REF] as it is especially designed to find a texture representation of the input (Figure 118). The network consists of 5 cascaded layers of convolutional filters to compute image features, with 2x2 kernels and LeakyReLU activations, followed by just one average pooling, with size equal to the size of final feature maps and three fully connected layers that act as the classifier. In contrast to standard CNN architectures adopted in the literature, this T-CNN does not integrate pooling layers operations between the convolutional layers and uses the smallest kernel size to help extracting the best texture representation. The T-CNN ends with a final 3-way SoftMax layer for pathological or normal classification. Again, DropOut is used to avoid overfitting. Note that this preprocessing step intends to normalize the different parameters associated with different scans described in the introduction. Thus, the training dataset consisted of 24265 image patches: 8978 normal, 8871 fibrosis+ground glass, 6416 emphysema. Similarly, 32x32 pixels patches were extracted from the test dataset as input for the T-CNN, centered to each pixel in the test image to which a T-CNN class will be assigned with a given probability.

Figure 119 shows confusion matrixes for both training and validation set (after the learning is finished). Also, the learning curves, accuracy and loss function, for both sets are shown in Figure 120. Note that after the epoch number 39, the model starts to overfit, as the loss and the accuracy for the validation set starts to increase and decrease respectively. The learning is early stopped at this point, saving the best suited set of parameters. Finally, Figure 121 shows the dimensionality reduction applied to T-CNN-transformed training patches (after all convolutional layers, and before the fully connected ones) to check if the learnt features helped to facilitate the classification task.

Classification using Random-CNN and style representation

This module tries to tackle the subclassification between ground glass and fibrosis developing a new texture descriptor. Basically, a texture is characterized by the arrangement of local patterns. The work in [START_REF] Gatys | A neural algorithm of artistic style[END_REF] demonstrated that CNNs are able to extract not only the content of an image but also the style representation, or the texture. The latter, of our interest, was achieved by computing the Gram matrix (𝐺) between the feature maps given a depth of the CNN. By including the feature correlations of multiple layers, they obtained a stationary, multi-scale representation of the input image, which captures its texture information but not the global spatial arrangement, which is more appropriate at describing textures than objects. In [START_REF] Ustyuzhaninov | Texture Synthesis Using Shallow Convolutional Networks with Random Filters[END_REF] it is shown that even with shallow CNNs with random filter, the style representation could be found for a given input image (see also §2.5.7). Such representation is used in the previous works to synthesize texture images based on examples. In the following, we will use this texture space as a feature description of the pathologies under analysis.

Merging these concepts, we implemented a Random-Multiscale network with seven different filter sizes 𝑓𝑥𝑓 with 𝑓 = 3; 5; 7; 11; 15; 23; 29, and 128 feature maps each (896 feature maps in total, one single convolutional layer) to transform the input data. Filters are obtained from a uniform distribution according to [START_REF] Glorot | Understanding the difficulty of training deep feedforward neural networks[END_REF]. Notice that no training is needed in contrast with §3.2.3.4 where the T-CNN had a training procedure. Instead, pathology signature is computed by storing the Gram matrices of the training patches used in §3.2.3.4. The Gram matrix is computed ( §2.5.7) for the resulting feature maps as follows: The element-wise mean and the standard deviation for the whole set of Gram matrixes obtained from patches used for training in section 3.2.3.4 are computed to obtain two signature matrixes (E k , and S k respectively, k = {1,2}) for each of the two classes to be detected. Finally, the classification is made by calculating the discrepancy δ k between the Gram matrix of the test patch (𝐺 𝑥 ) and the signature Gram matrix of the two classes, fibrosis and ground glass, applying the following formula:

𝐺 𝑖𝑗 = 1 M ∑ F ik . F jk M k=1 , [ 3-8 ]
δ k = ∑ ∑ |Gx ij -E ij k | 1+S ij k j i
, [START_REF] Newell | Development of quantitative CT lung protocols[END_REF][START_REF] Iwasawa | Assessment of prognosis of patients with idiopathic pulmonary fibrosis by computer-aided analysis of CT images[END_REF][START_REF] Coxson | Sources of variation in quantitative computed tomography of the lung[END_REF][START_REF] Lynch | Quantitative computed tomography of diffuse lung disease[END_REF][START_REF] Sorensen | Quantitative analysis of pulmonary emphysema using local binary patterns[END_REF][START_REF] Kurugol | Ranking and classification of monotonic emphysema patterns with a multi-class hierarchical approach[END_REF][START_REF] Hansell | CT staging and monitoring of fibrotic interstitial lung diseases in clinical practice and treatment trials: a position paper from the Fleischner Society[END_REF] where 𝑖, 𝑗 represent the index of the matrix element.

The minimum value for this computation indicates higher resemblance between textures and dictates the decision making.

E fibrosis E ground glass |E fibrosis -E ground glass|

Results

If we analyze the normalized confusion matrix for the whole test database (Figure 125), we notice a good prediction of the proposed method for normal and emphysema, but also a less discriminative power between fibrosis and ground glass patches. This can be explained by the effect of LCF filtering which also "flattens" the grayscale variation of the fibrosis regions making them more similar to ground glass (especially for most confusing situations). Another explanation of the quantitative result for fibrosis vs. ground glass comes from an effect of fuzzy decision at the interface between fibrosis and normal areas (patches falling on these interface zones will be more affected by a misclassification). There is always a transition zone between fibrosis and normal classes depicted as ground glass as it has its appearance, but it belongs to the same fibrotic entity semantically speaking. As in previous section, the 3D test database not being fully annotated, a complete quantitative evaluation of the classification result was not possible. A semi-quantitative evaluation of the performance on the annotated axial image subset of the test database was performed. Results are shown in Figure 126 All the sparse 2D annotated slices for one case (case -03) are shown with their respective automated classification to illustrate the stability of the method with respect to highly correlated slices (which is the case when analyzing slices for the same subject).

Note that the cascade if CNNs also responds well to the great variability between lobes (recall the effect produced by the gravitational gradient in the images), i.e. anterior-posterior robustness, and also to the variability in the appearance between top images (above the aortic arch) and end bottom ones (costophrenic angles). For this case, the confusion matrix is shown in Figure 128, and the qualitative result in Figure 129.

Discussion and comparative results

In this section, we compare the latest method versus modified configurations to discuss about mainly three points: (1) the benefits of using the cascade of CNN-based classifiers to progressively differentiate the classes over the state of the art T-CNN alone; (2) the gain of LCF filtering as a preprocessing step to enhance contrast between pathological and healthy tissue; (3) the impact on the outcome of the patch size parameter for both training and testing.

Benefits of using cascade of CNNs

In this section, we compare the developed cascade of CNNs approach with one of the best performing CNN-based state of the art method, T-CNN alone [START_REF] Anthimopoulos | Lung Pattern Classification for Interstitial Lung Diseases Using a Deep Convolutional Neural Network[END_REF], on the same database.

Figure 130 shows that the T-CNN confusion matrices for the training data are well adjusted for both cases with 3 (normal, emphysema, fibrosis+ground glass) and 4 classes (normal, emphysema, fibrosis, ground glass), after the learning process. However, it can be observed a better generalization for the first case, as it has a higher validation accuracy for each class. Examples of the classification methods can be observed in Figure 131 for a qualitative assessment versus cascade of CNNs showing a better consistency of the latter with respect to the ground truth [START_REF] Tarando | Cascade of convolutional neural networks for lung texture classification: overcoming ontological overlapping[END_REF]. Finally, in Figure 132 we can observe the confusion matrix for the second classifier, the Random-CNN which performs better than the regular T-CNN for the highest correlated classes. A semi-quantitative evaluation using sparse 2D axial slices showed an increase on accuracy of 10% of the proposed approach with respect to the state of the art T-CNN applied on the same 4 classes, no variation in TNR, an increase of 3% in the precision of the cascade approach compared to the original T-CNN classification outcome, and a slight increment of 1.24% for the accuracy parameter (these values are obtained on average between all the four classes). The detailed quantitative results are organized from Figure 133 to Figure 136

LCF filtering effect on classification

It has been discussed in the previous section the benefits of using the cascade of CNNs instead of the T-CNN alone to facilitate the texture recognition task progressively. In this case, we are interested on confirming the advantage of LCF filtering to enhance contrast between classes and helping the classification system. With this objective in mind, results for the cascade of CNNs with and without LCF module are compared and discussed in the following. Examples of the classification methods can be observed in Figure 137 for a quality assessment.

The quantitative comparison was performed once again in terms of sensitivity (True Positive Rate, TPR), specificity (True Negative Rate, TNR), Accuracy (ACC) and Precision (PREC).

The detailed quantitative results are organized from Figure 138

Accuracy -Ground Glass

Overall, a semi-quantitative evaluation using sparse axial images shows an increase for all the parameters of the proposed approach using LCF compared to the original cascade alone classification outcome: 10.25% in true positives detection, 3.68% in TNR, 5.3% in the precision, and 4.59% for the accuracy (these values are obtained on average between all the four classes). If we express the difference in performance between the proposed method (cascade+LCF) and the cascade alone, weighted by the proportion (volume) of each class in the ground truth, we obtain the results in Figure 142:

𝜆 𝑝 = ∆𝑃 [𝐶𝑙𝑎𝑠𝑠 𝑖 ] [𝑔𝑟𝑜𝑢𝑛𝑑 𝑡𝑟𝑢𝑡ℎ]
, [START_REF] Newell | Development of quantitative CT lung protocols[END_REF][START_REF] Iwasawa | Assessment of prognosis of patients with idiopathic pulmonary fibrosis by computer-aided analysis of CT images[END_REF][START_REF] Coxson | Sources of variation in quantitative computed tomography of the lung[END_REF][START_REF] Lynch | Quantitative computed tomography of diffuse lung disease[END_REF][START_REF] Sorensen | Quantitative analysis of pulmonary emphysema using local binary patterns[END_REF][START_REF] Kurugol | Ranking and classification of monotonic emphysema patterns with a multi-class hierarchical approach[END_REF][START_REF] Hansell | CT staging and monitoring of fibrotic interstitial lung diseases in clinical practice and treatment trials: a position paper from the Fleischner Society[END_REF][START_REF] Engeler | Ground-glass opacity of the lung parenchyma: a guide to analysis with high-resolution CT[END_REF][START_REF] Yamada | Airspace enlargement with fibrosis shows characteristic histology and immunohistology different from usual interstitial pneumonia, nonspecific interstitial pneumonia and centrilobular emphysema[END_REF] where 𝑖 = 1, … ,4 (for normal, emphysema, fibrosis, ground glass classes), ∆𝑃 = 𝑃 𝑐𝑎𝑠𝑐𝑎𝑑𝑒+𝐿𝐶𝐹 -𝑃 𝑐𝑎𝑠𝑐𝑎𝑑𝑒 is the difference between each of the four parameters 𝑃 (TPR, TNR, PREC, ACC) for the two methods, and [𝑥] denotes the cardinal of the set 𝑥 (here the number of pixels). When investigating the results case by case, we notice an altered performance in TPR for the fibrosis classification (cases 3, 7 and 10 in particular). If we analyze the normalized confusion matrix for each method computed for the whole test database (Figure 143), we notice a better prediction of the proposed method for normal, emphysema and even ground glass but also a less discriminative power between fibrosis and ground glass patches. This can be explained by the effect of LCF filtering which also "flattens" the grayscale variation of the fibrosis regions making them more similar to ground glass (especially for most confusing situations). Future work will target an improvement of LCF filtering so that the mentioned effect is minimized. Another explanation of the quantitative result for fibrosis vs. ground glass comes from an effect of fuzzy decision at the interface between fibrosis and normal areas (patches falling on these interface zones will be more affected by a misclassification). This can be seen for both methods in Figure 144 illustrating the worst cases of Figure 142 (cases 3, 7, 10): there is always a transition zone between fibrosis and normal classes depicted as ground glass. It could suggest that additional postprocessing might be needed to enforce the confidence of the classification.

It can also be seen from Figure 144 that, qualitatively, the proposed approach provides a more acceptable classification versus CNN cascade despite worse quantitative results for fibrosis (Figure 145). This is mainly due to the major improvement in detecting the normal class, which comforts the utilization of LCF as data preprocessing.

Finally, the precision of the ground truth definition for the quantitative assessment remains an open issue because of the complex task and inter-observer variability. We believe that a more accurate delineation of different classes in the ground truth will have a positive impact on the qualitative results of the proposed method, but this would require additional effort for radiologist experts, which is not straightforward to achieve.

Influence of patch size in CNN classification

We are focusing now on the effect of modifying the path size in the accuracy of the result. This verification is motivated by the systematic response given by the classifier while dealing with the transition between fibrosis and normal tissue. These pixels are normally assigned as ground glass. Figure 146 shows an example of this behavior. Also, note that the size of the patch is directly related to the maximum spatial sensibility of the method. At the same time, the need of having enough context of the analyzed pixel demands patch sizes large enough to contain relevant information at the best possible scale. The most popular choice in the literature is a set of train-validation-test of 32x32 pixels patches. Note that the actual space resolution is given by the pixel size (mm/pixel) times the 32x32 pixels patch size. In our case, as the pixel size is normalized for all the acquisitions (to 0.4 mm/pixel), the analyzed patches have a 12.8x12.8 mm resolution.

Figure 146. Illustration of the systematic response of the network with respect to the border regions between normal and fibrosis tissue (marked with red ellipses). Also, the maximum sensibility (in a spatial sense) can be appreciated, given a 32x32 patch size window used to scan the whole slice.

With the double objective of upgrading the system spatial sensibility and providing a better and smooth transition between two different textures, a reduced size of 16x16 pixels is considered (6.4x6.4 mm). From the training patches used in section 3.2.3.2, 16x16 pixels ROIs where extracted just cropping the contained (centered) patch from the original 32x32 ones.

The T-CNN was properly modified and retrained to this new dataset. Figure 147 Qualitative results for same cases shown in previous section are resumed in the next figure. Figure 150 shows the normalized confusion matrix for the 16x16 based system for the whole test database. We notice a similar prediction of the proposed method for 32x32 patch size, except for emphysema detection that has not performed as well as the previous choice of patch size. Qualitative results shown in Figure 149 seem to indicate a more accurate system in terms of spatial sensibility when using 16x16 patch size with respect to 32x32. As the patch size is smaller, it is less biased in the interface between two different classes making this transition sharper.

However, quantitatively, a slight decrease performance can be observed for all the four parameters, possibly because the system is more sensitive to small patterns (leading to a "noisier" output). Other modifications to the system shall be considered to tackle the transition regions as for example adding patches from annotated ground truth with less than 80% intersection with annotated ROIs. In this way, transition regions between normal and fibrosis are properly included in the training set.

Conclusion and future work

Lung texture classification has been proved to be a very challenging task mainly due to the great visual appearance overlap between the classes, image artifacts, intra-inter class variability, etc. Our analysis and understanding of this particular use case guided us towards a specific type of machine learning-based solution: ANNs (CNN in particular), a priori, known to be robust to this type of "noisy" data.

However, it was not obvious that CNNs could respond properly fitting our high dimension low sample medical data. The first experiment in this work aimed to verify this hypothesis showing that Deep Learning techniques, used even with reduced training databases, are still effective for refining and regularizing the output of a complementary classification approach for lung disease diagnosis [START_REF] Tarando | Increasing CAD system efficacy for lung texture analysis using a convolutional network[END_REF]. We presented a new CAD system for lung texture classification issued from combination between a 3D parametric approach and a CNN-based regularization showing increased performance with respect to the baseline technique.

Two main limitations related to the selected approach arose. First, it was verified the great impact of reduced database on CNNs generalization. We considered, thus, imperative to extend the training database (and hopefully increase the representativeness of each class within the data) to improve the supervised CNN performance which greatly relies on this information. Also, the highly visual correlated classes (mainly fibrosis and ground glass) remained hard to be identified.

A new texture-adapted architecture (instead of any state of the art object recognition network) was proposed to tackle the limitation in differentiating the classes: a CAD system for lung texture classification issued from a cascade of different CNN-based classifications allowing a hierarchical analysis making the system robust [START_REF] Tarando | Cascade of convolutional neural networks for lung texture classification: overcoming ontological overlapping[END_REF]. This choice is today supported by a recent publication [START_REF] Saikat | Deep neural networks for texture classification-A theoretical analysis[END_REF] demonstrating limitations of "classic" CNN architectures when applied to texturebased datasets (inherently higher dimensional compared to handwritten digits or other object recognition datasets), and, thus, suggesting the need of redesigning the network or enriching the system to learn textural features from input data. In this work we presented a classification based on a texture feature space that we considered proper for analyzing lung patches, using the features learnt by Random shallow CNNs. This approach was validated by comparison with the state of the art T-CNN showing a better response to the ontological overlap between the classes, especially fibrosis and ground glass.

It was observed that structural noise due to vascular tree had a great impact on biasing ROIs with high density values, thus, misleading the classifier. Normal structures as vessels and bronchi must be cleaned out from the image to be classified. To allow this, an original image pre-processing framework based on locally connected filtering applied in multiresolution was presented to help improving the learning process and boost the performance of CNN for lung texture classification [START_REF] Tarando | Boosting CNN performance for lung texture classification using connected filtering[END_REF]. By removing the dense vascular network from images used by the CNN for lung classification, locally connected filters provide a better discrimination between different lung patterns and help improving the CNN accuracy. However, it may affect some salient features of fibrosis patterns, making them more difficult to discriminate from ground glass. Future work should consider improving the LCF module to preserve discriminating features of pathological areas while still removing vascular structures.

A similar analysis could be added at the level of intrapulmonary vascular tree by investigating the structural remodeling as marker of emphysema (pruning of distal small caliber vessels, as observed in Figure 156 and Figure 157) or fibrosis (artifactual connection with fibrotic tissue). Finally, the patch scale used for the slice analysis might be also reviewed. The reason is that this type of analysis (by sliding patch-window through the slice) is time consuming due to number of operations to resolve (note that to each pixel inside the lung region corresponds its associated contextual patch of its surrounding). The network must be run separately for each patch, and there exists a lot of redundancy because of the patches overlapping. There is a trade-off between localization accuracy and the use of context (at different scales same patterns may appear visually different). Larger patches reduce spatial localization sensibility, while small patches allow the network to see only small context (as it was demonstrated in the discussion of section 3.2.3.7.3). An alternative solution that may have a good localization and use of context at the same time is the U-Net architecture, discussed in §2.5.6, providing an end-to-end classification.

Figure 1 .

 1 Figure 1. Statistics of ILDs (2013) [66].

Figure 2 .

 2 Figure 2. Schematic representation of the radiology workflow in relation with the PACS.

Figure 3 .

 3 Figure 3. Principle of the multi-slice CT imaging technique.

Figure 5 .

 5 Figure 5. CT scans from two subjects scanned using a low dose CT scan (40 mAs) (A and C) and a clinical CT dose (120 mAs)(B and D). The increased noise in the low dose CT scans (A and C) can be seen as linear "streaks" or a more "mottled" image while the clinical dose CT scans show a "smoother" image. Panel E shows the histogram of x-ray attenuation values from one of the subjects to illustrate that there is a difference in the extent of emphysema measured using either the threshold or percentile analysis (curve 𝒊 clinical dose CT scan, curve 𝒋 low dose CT scan). Example taken from[START_REF] Coxson | Sources of variation in quantitative computed tomography of the lung[END_REF].

Figure 6 .

 6 Figure 6. Examples of CLE are shown: a) mild CLE, scattered centrilobular lucencies, separated by large regions of normal lung, involving an estimated 0.5% of upper lung zone; b) moderate CLE, many well-defined centrilobular lucencies that occupy more than 5% of upper lung zone. Paraseptal emphysema is seen in anteromedial right and left lungs (arrows); c) Confluent CLE, multiple lucencies that span several secondary pulmonary lobules (circled in left lung) but are not associated with extensive hyperexpansion of secondary pulmonary lobules or distortion of pulmonary architecture; d) Advanced destructive emphysema, hyperexpansion of secondary pulmonary lobules with distortion of pulmonary architecture [2].

Figure 7 .

 7 Figure 7. PLE CT scan example shows widespread confluent areas of hyperlucency spanning one or several lobules. Some lobules, outlined by intact interlobular septa, appear hyperexpanded (arrowheads) [2].

Figure 8 .

 8 Figure 8. PSE examples in CT: a) mild PSE CT scan in smoker without COPD shows subpleural foci of low attenuation separated by intact interlobular septa along the mediastinum (arrows), measuring less than 1 cm; b) substantial PSE CT scan in patient with GOLD stage I COPD shows numerous well-demarcated areas of subpleural emphysema along chest wall and mediastinal pleural margins [2].

Figure 9 .

 9 Figure 9. Tissue destruction stages due to emphysema. Two sample patches per class are shown. CT intensity window was set to [-1400;200] HU. Severity increases from (a) to (d) with increasing size of low attenuating regions.

Figure 10 .

 10 Figure 10. Examples of all ILDs classes. Circles in colors indicate approximate pathological regions among healthy lung tissue. Red: fibrosis, blue: emphysema, pink: traction bronchiectasis, yellow: ground glass, violet: honeycombing, brown: AEF.

Figure 11 .

 11 Figure 11. Examples of artifacts that appears in CT imaging. First row examples show blurred images (motion artifacts due to breathing or heart motion -blue ellipses) and the second row shows an example of the gravitational gradient induced in gray scale values seen from anterior to posterior region (note that it could lead to a misclassification of emphysema in the anterior region, for example, or of ground glass in the posterior region where the density becomes higher -red ellipses).

2 . 3 .

 23 Pathological and anatomical variability inter-intra classes (in addition, some patients can have mixed restrictive/fibrotic and destructive/obstructive processes, such as in combined pulmonary fibrosis and emphysema syndrome); Lack of standardization: different image acquisition protocols, quality control, different image reconstruction techniques and radiation dose for different manufacturers, scanner calibration, scanner models, etc. (section 1.2);

.

Figure 12 .

 12 Figure 12. Examples of healthy tissue and some typical ILDs patterns: a) normal, b) emphysema, c) fibrosis, and d) ground glass. First row shows patches less correlated among classes. The examples in the second row are more correlated, showing the fuzzy border to identify these pathologies.

Figure 13 .

 13 Figure 13. Illustration of the overlapping among visual appearance of lung diseases under study. For example, AEF is supposed to exist between emphysema (CLE and PLE) and pulmonary fibrosis (in its different degrees of severity, including honeycombing). The intersections represent the ambiguous cases where the classification decision is not entirely clear.

Figure 14 .

 14 Figure 14. Illustration of classification workflow.

Figure 15 .

 15 Figure 15. Representations of Gabor filter bank for 5 frequencies (vertical axis), and 8 directions (horizontal axis), P=0.

Figure 17 .

 17 Figure 17. Schematic diagram of a one level decomposition (one iteration only) of DWT using filters: LPF (low pass filter), HPF (high pass filter) [78], ↓2 indicates subsampling. An example of standard dyadic discrete wavelet transform is shown in the next figure.

Figure 18 .

 18 Figure 18. Two iterations of a standard dyadic wavelet decomposition (b) of the original CT image (a). At each iteration j, the original image is downsampled by a factor of two to obtain a multiscale representation. Example taken from[START_REF] Depeursinge | Near-affine-invariant texture learning for lung tissue analysis using isotropic wavelet frames[END_REF]. The lowand high-pass decomposition filters combined (according to Figure17) are marked as LH, HL and HH.

Figure 19 .

 19 Figure 19. Illustration of the morphological opening (a) and closing (b) of f by the structuring element B on a 1D example. The alternating filters open-close and close-open are defined by:

Figure 20 .

 20 Figure 20. Example of morphological filtering of an axial lung CT image using a disk structuring element B of 3 pixel radius.

Figure 21 .

 21 Figure 21. Illustration of the morphological reconstruction by (a) dilation and (b) erosion of f by g on a 1D example.

  (a) reconstruction by opening (b) reconstruction by closing Figure 22. Examples of reconstruction by opening and closing using a SE disk of 3 pixel radius, applied on image in Figure 20 (a).

Figure 23 -Figure 23 .

 2323 Figure 23-25 show few examples of possible segmentations which can be obtained using morphological filtering.

Figure 24 .

 24 Figure 24. Example of morphological segmentation of ground glass regions based on 𝝆 𝒇 𝜹 operator.

Figure 25 .

 25 Figure 25. Example of pathological region segmentation (fibrosis) of Figure 20 (a) using morphological filters and watershed. 𝑩 𝒓 denotes the disk SE of radius 𝒓.

Figure 26 Figure 26 .

 2626 Figure 26. Illustration of a linear PCA approximation to a half-sphere circumscribed data: (a) two-dimensional principal component surface fit to the half-sphere data (left image); (b) projection of the data onto the first two principal components.Even if the procedure shows a successful cluster separation, for nonlinear problems, a nonlinear projection will respond better separating underlying differences between "classes"[START_REF] Hastie | The elements of statistical learning : Data mining, inference, and prediction[END_REF].

  -27 ] where I is the input image, 𝐱 is the center pixel, 𝒙 𝑝 = [-𝑅 sin ( + 𝐱 are P local samples taken at a radius 𝑅 around 𝐱, and 𝐻(•) is the Heaviside function. An example is shown in the next figure.

Figure 28 .

 28 Figure 28. Illustration of LBP. (a) The LBP relies on two parameters; the circle radius R and the number of samples P on the circle. (b) Local structure is measured w.r.to a given pixel by placing the center of the circle in the position of that pixel. (c) Samples on the circle are binarized by thresholding with the intensity in the center pixel as threshold value. Black is zero and white is one. The example image shown in (b) has an LBP code of 124. (d) Rotating the example image in (b) 𝟗𝟎° clockwise reduces the LBP code to 31, which is the smallest possible code for this binary pattern. This principle is used to achieve rotation invariance [43].

Figure 29 .

 29 Figure 29. Illustration of example curves and their correspondent parameters: first row (a,b) -MLD and standard deviation (examples taken from [17]); second row (c,d,e) -skewness; third row (f,g,h) -kurtosis values.

Figure 30 .

 30 Figure 30. Illustration of different deep model architectures enumerated in [51].

  5. Deep Convolutional Network (DCN) Recurrent Neural Network (RNN) Auto Encoder (AE) Long / Short Term Memory (LSTM) Boltzmann Machine (BM) Restricted BM (RBM) Deep Belief Network (DBN)

Figure 31 .

 31 Figure 31. An example of a two-class partitioning problem based on the corresponding posterior probability 𝑷(𝒘 𝒊 /𝒗) according to a feature vector 𝒗.

Figure 32 .

 32 Figure 32. Architecture of a pure fuzzy system.

Figure 33 .

 33 Figure 33. Architecture of a Takagi-Sugeno-Kang fuzzy system.

Figure 34 .

 34 Figure 34. Architecture of a fuzzy system with fuzzifier and defuzzifier.

Figure 35 .

 35 Figure 35. An example of SVM division plane for classification between two classes 𝑽 𝟏 and 𝑽 𝟐 : (a) linear separable problem; (b) using kernel function to map the support vectors to a higher space where the samples are linearly separable.

Figure 36 .

 36 Figure 36. Real neuron scheme on the left, and its artificial pair on the right. It consists of a nucleus or the cell body that integrates all the information received by the input dendrites and forwards the transmission down the axon to other neurons if the accumulated potential reaches a threshold.

Figure 37 .

 37 Figure 37. Example of a simple perceptron, where circles represents the neurons or units of the network, and the arrows the connections (all unidirectional).

- 35 ]Figure 38 .

 3538 Figure 38. Example of recursive partitioning with two input variables 𝑿 𝟏 and 𝑿 𝟐 : (a) decision tree with five terminal nodes, 𝑹 𝟏 -𝑹 𝟓 , and four splits; (b) partition of a two-dimensional feature space by applying rules in (a) to do the recursive binary splitting; (c) perspective plot of the prediction surface [47].

Figure 39 .

 39 Figure 39. Comparison among the impurity measures for binary classification problems. The image shows the node impurity measures for two-class classification, as a function of the proportion 𝒑 in class 2. Cross-entropy has been scaled to pass through (0.5, 0.5) [47].

Figure 40 .

 40 Figure 40. CNNs exploit spatially-local correlation by enforcing a local connectivity pattern between neurons of adjacent layers. In other words, the inputs of hidden units in layer m are from a subset of units in layer 𝒎 -𝟏, units that have spatially contiguous receptive fields[START_REF]Theano tutorial on CNN[END_REF].

Figure 42 .

 42 Figure 42. An example of max pooling is shown in the figure. In this case, a 3 size neighboring is pooled from one layer m-1 to the next layer m, taking the maximum value of them.

Figure 43 .

 43 Figure 43. RELU function curve y for a given input z.

Figure 46 .

 46 Figure 46. Example of a pooling operation. The receptive field of each unit is a 2x2 area in the previous layer's corresponding feature map. Each unit computes the average of its four inputs, multiplies it by a trainable coefficient, adds a trainable bias, and passes the result through a sigmoid function. Contiguous units have non overlapping contiguous receptive fields.

Figure 47 .

 47 Figure 47. Example of an output layer fully connected with the previous layer. Note that every unit 𝐱 𝐢 is connected with every unit in the next layer with its respective weight 𝐰 𝐢𝐦 . The output vector formed by the 𝐨 𝐦 represents the actual output of the network, a predicted class in a classification problem.

Figure 49 .

 49 Figure 49. Example of the input transformation through the network to obtain an output that represents a predicted class. For a given pattern μ, 𝒙 𝒏 are its pixel values represented by the grid in the figure. The CNN transforms this input vector and outputs its respective class value. The size of 𝑶 𝝁 is the number of classes to be categorized.

  𝑖 (1 -𝑦 𝑖 ) + ∑ 𝑡 𝑗 𝑦 𝑖 = -𝑡 𝑖 + 𝑦 𝑖 ∑ 𝑡 𝑗 = 𝑦 𝑖 -𝑡 𝑖 𝑗≠𝑖 .

Figure 52 .

 52 Figure 52. Illustration of fitting sample data set with increasing number of free parameters from a to c: (a) underfitted curve, (b) a good fit to noisy data, (c) overfitting of the same data, the fit is perfect on the training set (crosses in black in plots) but is likely to be poor on a test set represented by a circle.

Figure 53 .

 53 Figure 53. Example of how data visualization can help to identify appropriate analytics techniques, and clean the data, in this case, detection of outliers (red circle).

Figure 54 .

 54 Figure 54. Visualizations of 6,000 handwritten digits from the MNIST data set [52].

Figure 55 .

 55 Figure 55. Example of confusion matrix of a 9-class classification problem for a given dataset. Note that it can somehow quantify the correlation between classes (which classes are more difficult to classify).

Figure 56 .

 56 Figure 56. Example of accuracy and loss plots for training/validation sets while training. Note the parallelism: while the loss function decreases for each epoch (the difference between output and target decrease), the accuracy increases correspondingly.

Figure 57 .

 57 Figure 57. Examples of visualized weights for the first layer of a neural network are shown in the next figure: (a) noisy features could imply an issue (overfitting, improperly set learning rate); (b) smooth, clean and diverse features are a good indication that the training is performing well.

  [𝑥] denotes the cardinal of the set 𝑥 (here the number of pixels), 𝑇𝑃 = 𝑡𝑟𝑢𝑒 𝑝𝑜𝑠𝑖𝑡𝑖𝑣𝑒 = 𝐷𝑒𝑡𝑒𝑐𝑡𝑖𝑜𝑛 ∩ 𝐺𝑟𝑜𝑢𝑛𝑑 𝑇𝑟𝑢𝑡ℎ , 𝑇𝑁 = 𝑡𝑟𝑢𝑒 𝑛𝑒𝑔𝑎𝑡𝑖𝑣𝑒 = 𝐷𝑒𝑡𝑒𝑐𝑡𝑖𝑜𝑛 ∩ 𝐺𝑟𝑜𝑢𝑛𝑑 𝑇𝑟𝑢𝑡ℎ , 𝐹𝑁 = 𝑓𝑎𝑙𝑠𝑒 𝑛𝑒𝑔𝑎𝑡𝑖𝑣𝑒 = 𝐺𝑟𝑜𝑢𝑛𝑑 𝑇𝑟𝑢𝑡ℎ -(𝐷𝑒𝑡𝑒𝑐𝑡𝑖𝑜𝑛 ∩ 𝐺𝑟𝑜𝑢𝑛𝑑 𝑇𝑟𝑢𝑡ℎ) , 𝐹𝑃 = 𝑓𝑎𝑙𝑠𝑒 𝑝𝑜𝑠𝑖𝑡𝑖𝑣𝑒 = 𝐷𝑒𝑡𝑒𝑐𝑡𝑖𝑜𝑛 -𝐺𝑟𝑜𝑢𝑛𝑑 𝑇𝑟𝑢𝑡ℎ , 𝑁 = 𝑡𝑜𝑡𝑎𝑙 𝑛𝑢𝑚𝑏𝑒𝑟 𝑜𝑓 𝑝𝑖𝑥𝑒𝑙𝑠 𝑜𝑓 𝑡ℎ𝑒 𝑙𝑢𝑛𝑔 𝑚𝑎𝑠𝑘 .

Figure 58 .

 58 Figure 58. Architecture of LeNet-5 used in [21], a CNN for hand digits recognition. Convolutional layers (or the trainable feature detectors) are labeled as 𝑪, sub-sampling layer as 𝑺 and fully connected layer as 𝑭.

Figure 59 .

 59 Figure 59. Transforming fully connected layers into convolution layers enables a classification net to output a heatmap.Adding layers and a spatial loss produces an efficient machine for end-to-end dense learning[START_REF] Long | Fully convolutional networks for semantic segmentation[END_REF].

Figure 60 .

 60 Figure 60. U-net architecture (example for 32x32 pixels in the lowest resolution). Each blue box corresponds to a multichannel feature map. The number of channels is denoted on top of the box. The x-y-size is provided at the lower left edge of the box. White boxes represent copied feature maps. The arrows denote the different operations [56].

Figure 61 .

 61 Figure 61. In the CNN, a given input image is represented as a set of filtered images (feature maps) at each level in the CNN. It is possible to visualize the information at different levels of the CNN by reconstructing the input image from only the network's responses (feature maps content) in that particular layer. Note that reconstruction from lower layers is almost identical (a,b,c), and in deeper layers of the network, detailed pixel information is lost while the high-level content of the image is preserved (d,e). On top of the original CNN representations a new feature space that captures the style of an input image can be built (feature maps correlations). Style reconstructions of the input (1) are shown from (a) to (e) on top of the CNN architecture on different subsets of CNN layers (note that the information of the global arrangement of the scene was discarded) [25].

Figure 62 .

 62 Figure 62. Images that combine the content of a photograph (A) with the style of several artworks (bottom left corner of each panel). The images were created by finding an image that simultaneously matches the content representation of the photograph and the style representation of the paintings[START_REF] Gatys | A neural algorithm of artistic style[END_REF].

Figure 63 .

 63 Figure 63. Architecture of a single-layer Random-CNN used in [26] for texture synthesis. Each blue cube represents the set of 128 feature maps for each kernel with its respective size 𝒇.

  (a) Raw pixel values (b) VGG [25] (c) Random multi-scale Figure 64. Illustration of the similarity measures between textures for two different texture models as well as between raw pixel values of these textures. First row, examples of patches from ten different textures; second row, normalized Euclidean distances between raw pixel values (a) or Gram matrixes corresponding to the two models (b,c) for all pairs of textures. The matrix element (𝒊, 𝒋) corresponds to the median distance between patches from textures 𝒊 and 𝒋. The values in this figure are shown on a log-scale [26].

  As in the previous method, a new texture is generated from a uniform noise image and iteratively optimized to match the Gram matrix of the reference texture (squared distance between matrixes was the similarity measure function). Comparative examples are shown in the next figure.

Figure

  Figure[START_REF] Jain | Image Reconstruction from Projections[END_REF]. Samples synthesized from the two described models for three different reference textures (first column). The second column shows samples from the VGG model from[START_REF] Gatys | A neural algorithm of artistic style[END_REF], and the last column from the multi-scale model (with 1024 feature maps)[START_REF] Ustyuzhaninov | Texture Synthesis Using Shallow Convolutional Networks with Random Filters[END_REF].

Figure 66 .

 66 Figure 66. A screen shot of the graphical tool for the annotation of image regions.

  a. above the aortic arch, b. aortic arch level, c. carina, d. right inferior pulmonary vein, e. hepatic dome, f. end bottom, costophrenic angles.

Figure 67 .

 67 Six CT axial images spaced in the longitudinal direction throughout the lung and showing the characteristic regions for pathology annotation. Note also the variability in the lung cross-section size, imposing a large FOV for the CT acquisition: a) above the aortic arch, b) aortic arch level, c) carina, d) right inferior pulmonary vein, e) hepatic dome, f) end bottom, costophrenic angles.

Figure 68 .

 68 Figure 68. Axial resolution with respect to the field of view (FOV) for the same reconstruction matrix: (a) FOV capturing both lungs, 0.77 mm/pixel resolution, (b) FOV focused on the right lung, 0.35 mm/pixel resolution.

2 .Figure 69 .

 269 Figure 69. Distribution of the various lung tissue patterns in terms of hand-drawn ROIs and the patches extracted from them for each stage. Patches considered in (a) correspond to the training dataset (5576 in total). For the Stage 1 (b), the total number of ROIs patches is 54233. The high imbalance observed in (b) is tackled with data augmentation of those less numerous classes.

Figure 70 .

 70 Original images (top row) and radiological ROIs ground truth examples (bottom row) from our database. Colors indicate pathological/normal regions in the lung. Red: fibrosis, blue: emphysema, light blue: traction bronchiectasis, violet: honeycombing, brown: AEF, green: normal.

Figure 72 .

 72 Figure 72. Illustration of patch extraction from ground truth ROIs, in this case for a fibrosis region, indicated in pink (green: normal, violet: honeycombing).The next figures show some examples for normal, emphysema, fibrosis and ground glass patches from the database. They illustrate the great intra class variability both visually (in textural appearance) and gray level histogram distribution. It can also be observed the overlap between different classes which share similar phenotypes.

Figure 73 .

 73 Figure 73. Illustration of different normal patches from the database and their respective histograms that shows the great intra class variability (and overlap with other classes). From the first one to the last one, the overlap with other classes increases and the decision making becomes more complex.

Figure 74 .

 74 Figure 74. Illustration of different emphysema patches from the database. Respective histograms are also shown to illustrate the great intra class variability (and overlap with other classes). From the first one to the last one, the overlap with other classes increases.

Figure 75 .

 75 Figure 75. Illustration of different fibrosis patches from the database. Respective histograms are also shown to illustrate the great intra class variability (and overlap with other classes). From the first one to the last one, the overlap with other classes increases and the decision making becomes more complex.

Figure 76 .

 76 Figure 76. Illustration of different ground glass patches from the database. Respective histograms are also shown to illustrate the great intra class variability (and overlap with other classes).

Figure 77 .

 77 Figure 77. Illustration of different slice examples from the database and their densitometric evaluation. Top row -original images, second row -annotated ground truth, third row -histogram distribution starting from -1024 HU (peak remarked by the red ellipses in position -1024 HU), fourth row -same histogram from third row but starting from -1020 HU to avoid the noise peak located at -1024 HU.

Figure 78 .

 78 Figure 78. Illustration of different volume examples (same cases as Figure 77) from the database and their densitometric evaluation. Top row -histogram distribution (peak remarked by the red ellipses), second row -same histogram but starting from -1020 HU to avoid the noise peak located at -1024 HU.

Figure 79 .

 79 Figure 79. Illustration of the noise presented in the data, and the effect of applying 3D Gaussian smoothing to remove it. (a) original image slice; (b) lung field of image in (a) with white points representing the -1024 HU peaks locations; (c) image (a) after applying the Gaussian 3D filtering; (d) both histogram curves for original (red) and filtered (blue) images. The red curve has the peak -1024 HU whereas the blue one does not (the peak is not shown, to compare curves in the same scale).

Figure 81 .

 81 Figure 81. Visualization of the same patch samples organized according to resemblance after Barnes-Hut t-SNE. Note the cluster marked with the red ellipse. It contains patches of the subpleural region, which means patches that contain high density due to the presence of pixels outside the lung.

Figure 82 .

 82 Illustration of the two border configurations while scanning a whole slice (a). Orange squares represent subpleural region patches, and blue ones, those patches biased by vascular high densities. Regions preserved by lung mask segmentation showing high density border in vascular areas (b).

Figure 84 .

 84 Figure 84. Workflow of the parametric approach.

Figure 85 .

 85 Figure 85. Accuracy obtained for each dataset over the different permutations of the dataset (represented on the horizontal axis).

Figure 86 .

 86 Figure 86. Concept of relief flooding controlled by sink definition: (a) original image relief, 𝒇; (b) complete "flooding" of 𝒇; (c) selective "flooding" of 𝒇 with respect to a reference set 𝒀 playing the role of a "sink".

Figure 87 .

 87 Schematic representation of the concept of lung segmentation. Simultaneous lung flooding (b) may lead to the impossibility of separating lungs (either under-segmentation (e) or over-segmentation (f)). Alternative lung flooding (g, h) avoids these issues.

Figure 88 .

 88 Figure 88. Flowchart of the lung segmentation algorithm (contrasted lung regions).

3. 2 . 1 . 1 . 1

 2111 Detection of contrasted (low dense) lung regions by flooding The first step consists of extracting the lung fields linked by trachea using a flooding with sinks placed on volume border (Figure 89 (b)(c)), 𝑔 = 𝑅 𝑓 𝛿 (. , 𝑌) , [ 3-1 ]

  (a) original data, 𝑓 (b) relief flooding, 𝑔 (c) threshold of (b)-(a) (d) z-closing of lung relief (e) spherical opening of (d) (f) liniar closing (rough lungs) (g) trachea-removed lungs (h) "sink" set for each lung (i) flooding of (a) w.r. to "red" sink in (h) (j) threshold of flooded amount in (i) (difference (i)-(a)) (k) flooding of (i) w.r. to "blue" sink in (h) (l) threshold of flooded amount in (k) (difference (k)-(i)) (m) regularization of (j), (l) (n) lung segmentation of (a) (o) 3-D rendering of the result Figure 89. Illustration of the main steps of the 3D segmentation algorithm at a level of an axial image (see text).

Figure 90 .

 90 Figure 90. Illustration of an incomplete lung segmentation: (a) original image slice; (b) lung segmentation of (a).

  (a) lung field from flooding (b) "rough" lungs -red (c) trachea removed (d) reference "sink" set per lung Figure 91. 3D rendering of main partially results of the algorithm (see text).

Figure 92 .

 92 Figure 92. Illustrative flowchart for dense lung tissue retrieval: (a) ribs extraction (150 HU), (b) noise filtering and axial hole filling, (c) longitudinal (z-) linear dilation and closing, (d) enclose rib cage around lungs (morphological closing, spherical SE of variable size), (e) remove lungs detected previously (Figure 89 (m)), (f) add soft tissue filtered shape restricted to mediastinum and abdominal structures, (g) detect lung dense area as remaining region enclosed by the rib cage (green/yellow) not near soft tissue (white) and airways (pink), (h) final segmentation (in green, dense tissue added).

  Figure 93 (continues on next page). Some examples of lung segmentation: blue/red -right/left lung (contrasted areas), green -peripheral opacities. Top: axial view; middle: 3D rendering, front view; bottom: 3D rendering, rear view.

Figure 94 .

 94 Figure94. Illustration of patch filling procedure for one patch example marked with an orange square for both raw data and lung mask in (a). Top row shows the evolution of the raw data for each step of the procedure; and the bottom row, the evolution of the respective lung segmentation patch (where black -outside the lung; white -inside the lung; gray -pixels filled by the procedure in each step). (b)-(e) show the filling procedure performed symmetrically to each axis 𝒙, 𝒚 and both diagonals; and recursively for the four quadrants of the patch in (e) marked in blue squares (as it can be observed respectively from (f) to (i)).

Figure 95 .

 95 Figure 95. Example of patch-filling procedure: (a) original patch (32x32) on the lung border of image (d); (b) corresponding patch region in the lung mask. 0 value (black) means outside the lung; (c) outcome of the patch filling procedure: the part out of the lung mask is filled using the rest of the information in the same patch; (e) and (f) show the classification using the CNN described in §3.2.2.3 without and with patch filling respectively (green=normal, blue=emphysema, pink=fibrosis, white=ground glass), which proves the interest of the masking procedure.

Figure 96 .

 96 Figure 96. Limitations of a previous lung texture classification approach [13]: (a) ground glass/fibrosis, (c) fibrosis; (b), (d) classification results showing "noisy" output and mixed patterns in regions belonging to the same class (green=normal, blue=emphysema, pink=fibrosis, white=ground glass, yellow=airways, red=highly dense patterns or vessels).

  Figure98. Schematic workflow of the parametric 3D approach: (a) input (3D) image; (b) multi-resolution decomposition of the lung "relief" by successive morphological filtering with increasing filter size; (c) the tree graph structure for the hierarchic decomposition of the lung relief from (b); (d) each graph node is classed using a fuzzy-logic approach. The fuzzy membership functions are considered with respect to six features extracted from the pattern associated with the graph node: the gray-level mean, the standard deviation, the external mean, the pattern compactness, the pattern relative volume increase over scales and the neighborhood degree of the node. The x-axis represents the computed feature value and the y-axis represents the corresponding fuzzy analysis value; (e) decision integration over the scales; (f) final result.

Figure 100 .

 100 Figure 100. CNN's Confusion Matrix (training set) of tissue pattern classification.

Figure 101 .

 101 Figure 101. Trained filters examples for the first convolutional layer C1, used to scan the input patch image.

Figure 102 .Figure 103 .

 102103 Figure 102. Examples of mixed classification results: top row -original images, second row -annotated ground truth, third row -CLS classification, fourth row -CNN classification, bottom row -merged CNN-CLS classification (green=normal, blue=emphysema, pink=fibrosis, white=ground glass).

Figure 104 .

 104 Figure 104. Score results considering both original (parametric and non-parametric) and merged classification. The vertical axes represent the True Negative Rate for each class, and the horizontal axes the patient under study. Same behavior (on average) between CNN-CLS and CLS.

Figure 105 .

 105 Figure 105. Score results considering both original (parametric and non-parametric) and merged classification. The vertical axes represent the Precision for each class, and the horizontal axes the patient under study. A 6% increase (on average) is shown for CNN-CLS compared to the original CLS classification.

Figure 106 .

 106 Figure 106. Score results considering both original (parametric and non-parametric) and merged classification. The vertical axes represent the Accuracy for each class, and the horizontal axes the patient under study. A slight 0.4% increase (on average) is shown for CNN-CLS compared to the original CLS classification.

Figure 107 .

 107 Figure 107. Normalized Confusion Matrix for the whole test database for the three methods: (a) state of the art parametric CLS [13]; CNN alone; and (c) proposed hybrid classification merging CLS-CNN.

Figure 108 .

 108 Figure 108. Limitations of the state of the art lung texture classification approach T-CNN [27]: (a) emphysema/fibrosis, (c) fibrosis/ground glass; (b), (d) classification results showing locally mixed patterns (green=normal, blue=emphysema, pink=fibrosis, white=ground glass).

Figure 109 .

 109 Figure 109. Lack of CT acquisition normalization between heterogeneous patient datasets leading to a strong variability of normal tissue appearance (red circles). Idem for the rest of pathologies, leading to a strong variability inter and intra classes not related to their ontological definition.

Figure 110 .

 110 Figure 110. Top: Flowchart of the system proposed. The CNN Cascade is pre-tuned with the 2D annotated database. LCF is applied to both training and test slices, from which patches are extracted. Bottom: flowchart of the cascade system. The T-CNN is pre-trained from scratch with the 2D annotated database pre-filtered by LCF. Then, test images from the database are fed into the system in a patch wise mode. The patches classed as fibrosis+ground glass are separated using Random-CNN to compute a texture signature to improve this subclassification.

Figure 111 .

 111 Figure 111. Illustration of the effect of a locally-connected filter on a ROI of an axial image using different window size. (a) Original axial CT showing fibrosis regions and ROI selection; (b) enlarged ROI (including only lung area) and dual representation as topographical relief. The focus here is on a group of vessels (dotted circle) w.r. to fibrosis reticulations (top); (c) LCF of (b) with a window half-size of 3 pixels. Small vessels are removed with minimal impact on fibrosis area; (d) LCF of (b) with a window half-size of 5 pixels. Larger vessels are removed with the cost of a more pronounced relief "flattening" in the fibrosis area.

  Figure 112. Example of LCFf  filtering (b,e) on noisy grayscale (top) and binary (bottom) images (a,d) versus median filtering (c,f) of same size.

[ 3 - 6 ]

 36 with kd denoting the border subset of the 2D spatial neighborhood of size k orthogonal to the d direction, and C18 the set of orientations for the 18-connectivity (Figure113 (b)).

  Figure 114 (b) illustrates the filtering effect of VLCF at the level of an axial CT image for 𝑘 = 3. The vascular structures are further detected by adaptive thresholding and directional reconstruction by erosion applied to the difference f -VLCFf  , Figure114 (c,d). Note that the sheetlike structures are minimally affected by the filter (in the example, the lung fissures).

  Figure 114. Illustration of VLCF on a lung CT image, for 𝒌 = 𝟑: (b) filtered image, (c) difference, (d) vascular structures up to size 𝒌.

Figure 115 .

 115 Figure 115. Flowchart of the multiresolution VLCF filtering. ˄ and ˅ denote infimum and supremum, while ↓ and ↑ decimation and interpolation, respectively. 𝒇 𝑽𝑳𝑪 is the filtering result and 𝒗 𝒇 vascular regions for all scales.

  Figure 116. Result of vascular removal/detection in Figure 114 (a) using multiresolution VLCF: (a) filtered image, (b) vascular information across all scales, (c) vascular segmentation, (d) lung and vessels mask.

Figure 117 Figure 117

 117117 Figure117illustrates two examples of 3D lung masks including the vascular structures (3D volume rendering). Both filtered images and 3D lung mask (that also contains airway-vessel information) will be used in the classification input.

Figure 118 .

 118 Figure 118. Architecture of the proposed T-CNN described in [27], specially designed for lung texture classification. The set of training patches is used to train this T-CNN for 3-class pathology classification. At this point of the research we count on 49 sets of MDCT of Avicenne database rescaled to match 0.4 mm/pixel spacing[START_REF] Anthimopoulos | Lung Pattern Classification for Interstitial Lung Diseases Using a Deep Convolutional Neural Network[END_REF], from which 32x32 pixels non-overlapping patches were extracted. We consider those patches falling 80% inside the marked ROIs. As usual, horizontal flipping and rotations are applied to artificially increase the training set. Furthermore, the native pixel values for each patch were set in the range [-1000, +1000] HU, and mapped to [0, 1] interval as follows:

Figure 119

 119 Figure 119. T-CNN confusion matrix for (a) training set (24265 patches in total) and (b) validation set (1075 patches in total) after the training.

Figure 120 .

 120 Figure 120. T-CNN training curves: loss and accuracy (mean for all classes expressed in %) vs. number of epochs (an entire pass for all training samples). The learning procedure is stopped at epoch 39 (dotted line).

  Figure 121. t-SNR of training patches transformed in feature vectors after all convolutional layers in the T-CNN. It can be observed that the T-CNN managed to create discriminative features to disentangle the differences between the classes. Note that fibrosis and ground glass are merged in one class. In the next figure, we show the final 2x2 pixel filters automatically learned by the first layer of T-CNN to disentangle the patch texture differences between the possible classes.

Figure 122 .

 122 Figure 122. Trained 2x2 filters examples for the first convolutional layer, used to scan the input patch image.

whereF

  np = 𝑛_𝑡ℎ 𝑓𝑒𝑎𝑡𝑢𝑟𝑒 𝑚𝑎𝑝 𝑣𝑎𝑙𝑢𝑒 𝑎𝑡 𝑝𝑜𝑠𝑖𝑡𝑖𝑜𝑛 𝑝, 𝑀 = ℎ𝑒𝑖𝑔ℎ𝑡 𝑡𝑖𝑚𝑒𝑠 𝑡ℎ𝑒 𝑤𝑖𝑑𝑡ℎ 𝑜𝑓 𝑡ℎ𝑒 𝑓𝑒𝑎𝑡𝑢𝑟𝑒 𝑚𝑎𝑝.

Figure 123

 123 Figure123shows this transformation of the raw data to the new space.

Figure 123 .

 123 Figure 123. Architecture of Random-CNN used for ILD patch texture representation. Each blue cube represents the set of 128 feature maps for each kernel with its respective size 𝒇. The Gram matrix is calculated as the new representation of the input (with texture information).

  Figure 124 illustrates the signature matrix 𝐸 for both classes and the absolute difference between them.

Figure 124 .

 124 Figure 124. Top: signature matrix E visualization for both classes and the absolute difference between them. Bottom: vector representation of the symetrical matrix for fibrosis, and the difference matrix.

Figure 125 .

 125 Figure 125. Normalized confusion matrix for the whole test database for the cascade of CNNs in combination with the LCF module.

  and Figure 127 below.

Figure 126 .Figure 127 .

 126127 Figure 126. Average LCF + cascade of CNNs (for all test cases) score parameters for each class. Averages values: TPR=51.6%, TNR=86.4%, PREC=41.2%, ACC=84.6%.

  Figure 128. Normalized (top) and absolute (bottom) confusion matrix for the entire test case (1606238 analyzed pixels inside the lung in total for the 10 slices) obtained using cascade of CNNs with LCF module.

Figure 129 (

 129 Figure 129 (continues on next page). Examples of classification results for 11 sparse 2D slices for one test case: top row -original images, second row -annotated ground truth, third row -cascade of CNNs method for 4 classes (normal, emphysema, fibrosis, ground glass). Color coding -green: normal, white: ground glass, blue: emphysema, pink: fibrosis, red: high density patterns (vessels), yellow: airways.

Figure 129 (

 129 Figure 129 (cont.). Examples of classification results for 11 sparse 2D slices for one test case.

Figure 130

 130 Figure 130. T-CNN's confusion matrix of tissue pattern classification for the same small validation (945 patches) set for 3 and 4 classes. On the right, we can observe the big overlap between Fibrosis and Ground Glass, and the lack of generalization for those classes, giving a low 10% accuracy for Ground Glass. On the other hand, the 3-class T-CNN keeps the accuracy high over 95% on average.

Figure 132 .

 132 Figure 132. Accuracy for the subclassification done by the Random-CNN for the small validation set used in Figure 128 (for fibrosis and ground glass represents an amount of 571 patches).

  below.

Figure 133 .

 133 Figure 133. Score results considering both state of the art T-CNN and proposed classification. The vertical axes represent the True Positive Rate for each class, and the horizontal axes the patient under study. A 10% increase (on average) in true positives detection is shown for the cascade of CNNs compared to the original T-CNN classification.

Figure 134 .

 134 Figure 134. Score results considering both state of the art T-CNN and proposed classification. The vertical axes represent the True Negative Rate for each class, and the horizontal axes the patient under study. Same behavior (on average) between the cascade of CNNs and T-CNN.

Figure 135 .

 135 Figure 135. Score results considering both state of the art T-CNN and proposed classification. The vertical axes represent the Precision for each class, and the horizontal axes the patient under study. A 3% increase (on average) is shown for the cascade of CNNs compared to the original T-CNN classification.

Figure 136 .

 136 Figure 136. Score results considering both state of the art T-CNN and proposed classification. The vertical axes represent the Accuracy for each class, and the horizontal axes the patient under study. A slight 1.24% increase (on average) is shown for the cascade of CNNs compared to the original T-CNN classification.

  to Figure 141 below.

Figure 138 .

 138 Figure 137. Examples of classification results for three different test cases (a) case -02 (b) case -03 (c) case -06. Top roworiginal images, second row -filtered image, third row -annotated ground truth, fourth row -cascade of CNNs method (over original DICOMs images) for 4 classes (normal, emphysema, fibrosis, ground glass), fifth row -Proposed approach for filtered images and the same classes. Color coding -green: normal, with: ground glass, blue: emphysema, pink: fibrosis, red: high density patterns (vessels), yellow: airways.

Figure 139 .Figure 140 .

 139140 Figure 139. Score results considering both cascade of CNNs with and without the LCF module. The vertical axes represent the True Negative Rate for each class, and the horizontal axes the case under study. A 3.68% increase (on average) is shown when using LCF.

Figure 141 .

 141 Figure 141. Score results considering both cascade of CNNs with and without the LCF module. The vertical axes represent the Accuracy for each class, and the horizontal axes the patient under study. A 4.59% increase (on average) is shown when using LCF.

Figure 142 .

 142 Figure 142. Average performance improvement of the proposed method (cascade+LCF vs. cascade), weighted by the proportion of each class in the lungs.

  Figure 143. Normalized Confusion Matrix for the whole test database for both methods: (a) cascade of CNNs; and (b) proposed classification adding the LCF module.

Figure 145 .

 145 Figure 144. Example of collateral negative effects of LCF applied to test cases: (a) case -03, (b) case -07, (c) case -10. Top row -original images, second row -filtered image, third row -annotated ground truth, fourth row -cascade of CNNs method (over original DICOMs images) for 4 classes (normal, emphysema, fibrosis, ground glass), fifth row -proposed approach for filtered images and the same classes. Same color coding as in Figure 137.

  shows confusion matrixes for both training and validation set (after the learning is finished). The learning procedure is early stopped on epoch 35. Learning curves, accuracy and loss function, for both sets are shown in Figure 148.

Figure 147

 147 Figure 147. T-CNN confusion matrix for (a) training set (24265 patches in total) and (b) validation (1075 patches in total) set after the learning process.

Figure 148 .

 148 Figure 148. 16x16 input T-CNN training curves: loss and accuracy (mean for all classes expressed in %) vs. number of epoch. Training is stopped at epoch 35 (dotted line).

  Figure 149. Examples of 16x16 patch classification results for same three test cases of Figure 137 (a) case -02 (b) case -03 (c) case -06. Top row -original images; second row -annotated ground truth; third row -LCF + cascade of CNNs method, with 32x32 pixels input patches for 4 classes (normal, emphysema, fibrosis, ground glass); fourth row -LCF + cascade of CNNs method, with 16x16 pixels input patches for same 4 classes. Color coding -green: normal, white: ground glass, blue: emphysema, pink: fibrosis, red: high density patterns (vessels), yellow: airways.

AFigure 150 .

 150 Figure 150. Normalized confusion matrix for the whole test database for both methods: (a) state of the art cascade of CNNs [45]; and (b) proposed classification adding the LCF module.

Figure 151 .

 151 Figure 151. Score results considering LCF + cascade of CNNs for 32x32 and 16x16 patch size. The vertical axes represent the True Positive Rate for each class, and the horizontal axes the test case under study. A -4.24% variation (on average) in TPR is shown for the smaller patch compared to the original 32x32.

Figure 152 .

 152 Figure 152. Score results considering LCF + cascade of CNNs for 32x32 and 16x16 patch size. The vertical axes represent the True Negative Rate for each class, and the horizontal axes the case under study. A -0.63% variation (on average) in TPR is shown for the smaller patch compared to the original 32x32.

Figure 153 .

 153 Figure 153. Score results considering LCF + cascade of CNNs for 32x32 and 16x16 patch size. The vertical axes represent the Precision for each class, and the horizontal axes the patient under study. A -3.46% variation (on average) in TPR is shown for the smaller patch compared to the original 32x32.

Figure 154 .

 154 Figure 154. Score results considering LCF + cascade of CNNs for 32x32 and 16x16 patch size. The vertical axes represent the Accuracy for each class, and the horizontal axes the patient under study. A -0.04% variation (on average) in TPR is shown for the smaller patch compared to the original 32x32.

Figure 156 .

 156 Figure 156. Vascular remodeling in emphysema: distal vessel pruning correlates with emphysema severity in COPD [86].

Figure 157 .

 157 Figure 157. Correlation with emphysema score LAA%-950HU [86].

Table 1

 1 summarizes the four methods discussed which represent the best results found in the literature.

	Class [Method]	Target space	Best performance	Advantages	ILD categories
	(1) [Histogram signatures + CVM (Cramer Von Mises measure] Distance) as dissimilarity	3D	Correlation mild to moderate with radiologist classes) [19] volumetric proportion of scoring (quantifies total	Fast, dimensionality glyph reduction, visual	N,E,GG,R, H
	(1) [GLH+Wavelet and SVM]	2D	76,9% accuracy on test patches [17]	Affine invariance, dimensionality reduction	N,E,GG,F, M
				Texture description,
	(2) [Texture-CNN (ANN classifier integrated)]	2D	85,5% accuracy on test patches [27]	dimensionality reduction, automatic	N,E,GG,R, H,M,C
				feature learning
	(2) [Convolutional Restricted Boltzman Machine + SVM]	2D	77% accuracy on test patches [42]	Translation invariance, automatic feature learning	N,E,GG,F, M

Table 1 .

 1 Summary of methods and their results for latest state of the art works for ILD classification. N: normal, F: fibrosis, E: emphysema, GG: ground glass, H: honeycombing, R: reticulation, C: consolidation, M: micronodules.

Table 3 .

 3 Distribution of the lung tissue patterns per patient for the Stage 1 of the database.

			#PATCH 32x32 Normalized 0,4mm/pixel (80% in ROI) NO-Overlap
	Case	AEF	Airways	Emphysema	Fibrosis Ground Glass	Honeycombing	Normal
	1	0		34	160	19	227	1397
	2	0		16	15	3	83	1513
	3	3		0	0	49	0	654
	4	0		0	37	2	1	943
	5	0		20	2	0	0	593
	6	2		1	71	66	6	1449
	7	8		1	13	15	0	1036
	8	60		103	134	178	18	923
	9	0		0	98	30	4	608
	10	3		75	303	175	417	1362
	11	10		14	134	2	398	414
	12	0		0	522	0	8	2875
	13	0		0	48	6	0	844
	14	0		0	9	72	0	609
	15	0		0	200	3	2	2095
	16	0		5	0	3	0	1520
	17	0		0	345	110	0	2404
	18	0		0	26	115	0	874
	19	0		0	163	0	0	2565
	20	0		0	77	4	0	627
	21	0		0	64	70	5	599
	22	15		0	272	1	2	2894
	23	0		0	59	27	1	641
	24	0		0	290	18	0	2062
	25	0		0	0	616	0	2620
	26	0		0	173	37	0	2239
	27	0		13	733	26	0	1848
	28	0		0	50	0	0	1217
	29	0		13	20	16	0	1854
	30	0		128	107	71	0	512
	31	0		0	28	0	0	1477
	32	0		0	33	24	94	524
	33	0		0	7	23	0	1612
	34	0		0	87	0	0	948
	35	0		3	17	4	0	161
	36	0		33	8	0	0	57
	37	0		0	8	3	0	63
	38	0		0	2	18	0	316
	39	0		0	22	6	0	154
	40	0		0	0	3	0	245
	41	0		0	0	5	0	188

Table 4

 4

. Acquisition protocol details for each case in the database (Stage 1).

Table 5 .

 5 Mean and standard deviation for each dataset.

Upsampling

t-SNE visualization

t-SNE is applied to the complete database (stage 1) to illustrate the relationship between the classes. Figure 80 and Figure 81 show the outcome of the process applied to the set of original patches from the whole database (as the normal patches are much more numerous with respect to the other classes, only some of them are considered).

Figure 80. Barnes-Hut t-SNE applied to raw patches from the database.

Classification using CNN

Due to the access to only 2D samples in the available training database, the convolutional network was designed to perform on a two-dimensional basis at the level of axial images.

The chosen architecture (shown in Figure 99) presents a reduced capacity with respect to other deep architectures in the state of the art, meaning a reduced number of parameters to be learned, which helps avoiding the over-fitting problem given the dimension of our training set. Following the standard CNN architectures adopted in the literature, we used two cascaded layers of convolutional filters to compute image features, combined with max-pooling operations and ended with fully-connected layers. Our CNN ends with a final 4-way SoftMax layer for pathological or normal classification. Thus, the outcome is a vector that contains the probability for each patch to be classed as each pathology. DropOut is used to avoid overfitting [START_REF] Krizhevsky | Imagenet classification with deep convolutional neural networks[END_REF]. Torch7 framework is used as it efficiently trains the CNN using GPU acceleration [START_REF] Collobert | Torch7: A matlab-like environment for machine learning[END_REF].

The learning parameters related to this technique were selected based on typical values and modified according to empirical trials. The system is trained to minimize an overall loss function (cross-entropy) which measures the probability of an erroneous answer (Gradient-Based Learning).

Figure 99. Architecture of CNN used for ILD patch classification. Each square represents a feature map, that is, the output for each unit. Neurons extract primitive visual features (edges, corners, etc.) which are then combined in subsequent layers to generate high or more abstract features.

The preliminary dataset (stage 0) for this study includes both the public ILD database detailed in [START_REF] Depeursinge | Building a reference multimedia database for interstitial lung diseases[END_REF] (113 patients with sparse 2D annotated axial images), and a proprietary database of Avicenne Hospital, Bobigny, France consisting of 12 volumetric acquisitions with sparse 2D axial annotations. The CT scanning protocol details are listed below. For the ILD database: slice thickness 1-2 mm; modality HRCT; contrast agent none; axial pixel matrix 512x512; x, y spacing 0.4-1 mm. For Avicenne database: slice thickness 1.25 mm; modality MDCT; contrast agent none; axial pixel matrix 512x512; x, y spacing 0.716 mm.

The relatively large annotated training sets used in this work made possible to train our CNN for a 4-class (normal, fibrosis, ground glass, emphysema) lung tissue classification. We used the ILD database only for training the CNN since it was no compatible with the CLS approach requiring volumetric data as input. We added two cases from Avicenne database to the training set to balance the number of patches for each possible class/pathology, and thus avoid biasing the network. All these cases are not considered for testing. Another set of 10 different DICOM image series of whole thorax (CT scans) from Avicenne database are used as test cases. Note that, as the visual overlap is insurmountable for some cases, the improvements of the automated method and its robustness represents a baseline classification that will serve as a starting point to be refined with medical background. As future work, this classification module will be improved to further differentiate between subclasses (and include classes not taken into account so far) by incorporating into analysis the medical context and complementary image processing.

The main direction which will be explored is to include the presence of structured abnormalities in the assignment of a specific pathology class. In this respect, we shall exploit the bronchial tree segmentation on which the presence of traction bronchiectasis will be detected (as marker of fibrosis, see Figure 155). Based on this, a probability map will be built-up assigning higher probability of fibrosis in the surrounding regions of bronchiectasis.