Computational Techniques for Automated Analysis of Animal Tissue Histological Images

Faculdade de Engenharia da Universidade do Porto Computational techniques for automated analysis of animal tissue histological images Frederico A. R. B. Junqueira Master Thesis Integrated Masters in Bioengineering Supervisor: Prof. João Manuel R. S. Tavares Associate Professor of the Mechanical Engineering Department, FEUP Co-Supervisor: Augusto Manuel Rodrigues Faustino Associate Professor of the Pathology and Molecular Immunology Department ICBAS June 2015 ii iii Agradecimentos Em primeiro lugar quero agradecer ao Professor João Tavares pela disponibilidade e aconselhamento fornecidos ao longo destes últimos meses. Sem a sua orientação este trabalho não seria possível. Quero agradecer também ao Professor Augusto Faustino pela sua disponibilidade, simpatia e vontade de ajudar. Quero agradecer do fundo do coração à Lia, que me acompanhou a 100% nesta jornada longa e nunca me deixou perder o Norte. Sem o seu carinho, muita paciência e apoio incondicional este trabalho não estaria completo. À minha mãe, ao meu pai, à minha irmã, ao meu irmão e a minha avó pelas gargalhadas, carinho e todo o apoio nesta etapa final. Ao Pedro pelos bons conselhos e ombro amigo sempre que precisei. Ao companheiros da tese Morgana, Jessica e Ricardo por me ouvirem nos momentos de stress. iv v “Everybody has a plan until he gets punched in the face.” Myke Tyson vi vii Resumo O estudo de tecidos celulares fornece uma incontestável fonte de conhecimento e compreensão acerca do corpo humano e do ambiente que o rodeia. Aceder a esta informação é, portanto, crucial para determinar e diagnosticar uma grande variedade de patologias, detetáveis somente ao nível microscópico. A histologia desempenha um papel importante na observação de células e suas características anatómicas, e igualmente para o diagnóstico clínico de patologias involvendo uma anormal conformação celular. Nas imagens histológicas, algoritmos de segmentação semi-automáticos ou automáticos são capazes de separar e identificar estruturas celulares de acordo com as suas diferenças morfológicas. Estes algoritmos de segmentação são a primeira abordagem a sistemas de visão computacional e, no que respeita à histopatologia, o diagnóstico automático de imagens histológicas. Como as amostras histológicas têm uma espessura reduzida, as características volumétricas são quase imperceptíveis, correespondendo a perdas de informação valiosas, principalmente topográficas e volumétricas, críticas para um correcto diagnóstico. Consequentemente, a combinação de algoritmos de segmentação e reconstrução 3D aplicados a datasets de imagens histológicas fornecem uma maior informação acerca da patologia analisada e estruturas microscópicas, destacando regiões anormais. Tendo isto em consideração, o presente trabalho focou-se em desenvolver algoritmo computacional automático capaz de realizar reconstrução 3D de superfícies de tecidos relevantes em secções histológicas 2D. Uma primeira abordagem foi desenvolvida focada em destacar as estruturas relevantes nas secções de tecido. Depois, um estudo feito com base em algoritmos de registo de imagem foi levado a cabo para descobrir qual a metodologia mais indicada para alinhar as secções provenientes de datasets de imagens. Combinando os melhores métodos de processamento e registo de imagem (resultado DICE para caso 1: ; resultado DICE para caso 2: ; resultado DICE para caso 3: ), avaliados em abordagens prévias, em conjunto com um algoritmo de reconstrução 3D foi possível uma representação volumétrica das estruturas de tecidos pertinentes do dataset de imagens alinhadas. viii ix Abstract The study of cellular tissues provides an incontestable source of information and comprehension about the human body and the surrounding environment. Accessing this information is, therefore, crucial to determine and diagnose a wide variety of pathologies detectable only at a microscopic scale. Histology plays an important role in the observation of cells and their anatomical features, and so for clinical diagnosis of all the pathologies involving abnormal cellular conformation. In the histological images, semi-automated or automated segmentation algorithms are able to separate and identify cellular structures according to morphological differences. These segmentation algorithms are the first approach for computational vision systems and, concerning histopathology, the automated diagnose of histological images. Since the histological samples are thin, the volumetric features are almost unnoticeable, corresponding to losses of valuable information, mainly topographical and volumetric data, critical for a correct diagnostic. Hence, the combination of segmentation and 3D reconstruction algorithms applied to histological image datasets provides more information about the analysed pathology and microscopic structures, highlighting abnormal areas. Taking this into consideration, the present work focussed on developing an automatic computational algorithm capable of performing the 3D surface reconstruction of relevant tissue structures of 2D histological slices. A first approach was developed fixed on highlighting the relevant structures from the tissue sections. After that, a study on image registration algorithms was conducted to find the most suited methodology to align the slices from histological image datasets. Combining the top-performing image processing and registration methods (DICE score: 0.9267±0.0337 for Case 1; 0.9367±0.0356 for Case 2; 0.9683±0.0283 for Case 3), evaluated in the previous approaches, with a 3D surface reconstruction algorithm it was possible to calculate a volumetric representation of the pertinent tissue structures from the aligned image dataset. xvi xvii List of Tables Table I. Characterization of the three image datasets studied in this work. ..................... 38 Table II. Table containing the mean DICE scores and the elapsed times for the intensitybased registration implementation, with different types of transformation (rigid, similarity and affine) and models (reference slice, pairwise), on images from both top-performing pre-processing methods (YCbCr color transformation and CLAHEred). ..................................................................................................... 66 Table III. Table containing the mean DICE scores and the elapsed times for the featurebased registration implementation with pairwise model on images with both smoothed and non-smoothed image normalization pre-processing. ......................... 69 Table IV. Table containing the mean DICE scores and the elapsed times for the intensitybased non-rigid registration implementation, with two different algorithms (B-spline and Demon) and models (reference slice, pairwise), on images obtained from CLAHEred pre-processing method. ......................................................................... 71 Table V. Table containing the mean DICE scores and the elapsed times for the intensitybased rigid registration method implementation, on two different image datasets (Case 1 and 2), previously pre-processed by CLAHE-red algorithm. ......................... 73 Table VI. Table containing the mean DICE scores and the elapsed times for the intensitybased rigid registration method implementation, on all the image datasets (Case 1, 2 and 3), considering two different grayscale images (hematoxylin and eosin images), obtained through pre-processing (CLAHE-red and YCbCr transformation) and segmentation (stain deconvolution algorithm). ................................................. 77 xviii xix Abbreviations and Acronyms CNS Central Nervous System CT Computer Tomography ECM Extracellular Matrix H & E Haematoxylin and Eosin HSV Hue – Saturation - Value ITK Insight Segmentation and Registration Toolkit µm Micrometer mm Millimeter MRI Magnetic Resonance Imaging nm Nanometer PCA Principal Component analysis PET Positron Emission Tomography PNS Peripheral Nervous System 3D Three Dimensions 2D Two Dimensions RGB Red – Green - Blue CLAHE Contrast Limited Adaptive Histogram Equalization OD Optical Density xx 1 Chapter 1 Introduction 1.1 - Motivation Histological studies provide an important help in the understanding of some complex pathophysiological processes concerning diseases at the cellular scale. These studies are considered the gold standard for assessing the natural response of a cellular tissue in face of a pathology or therapeutic intervention (Chakravarty, Bedell et al. 2008). To produce histopathology slides, a rather complex protocol must be executed involving a substantial amount of human labor and information processing (Randell, Ruddle et al. 2012). Although in-vivo imaging techniques, such as the MRI and PET, assess anatomical and pathological information without invasive procedures, they require extensive validation when compared to histological ex-vivo examinations (Chakravarty, Bedell et al. 2008). Visual interpretation, the core of most medical diagnostic procedures and the final diagnostic decision for cancer and other diseases, is based on tissue examination. This method requires a long time, intensive manual labor to produce viable results and presents a sampling bias that promotes intraand inter-reviewer discrepancies when analysing histological tissues (Sertel, Kong et al. 2009). Thus, it is clear the need for automated processes concerning morphology diagnostics in medicine, to improve the diagnostic accuracy and provide a fast and reliable second opinion to histopathologists. Automated systems can reduce human factor mistakes and increase the speed of diagnostic processes (Nedzved, Belotserkovsky et al. 2005). The volumetric data analysis from relevant tissue structures visible from 2D histological slices is not often a straightforward process, requiring a great amount of experience from histopathologists (Koshi, Holla et al. 1997). Therefore, three dimensional reconstruction of tissue samples at a microscopic resolution reveals significant potential to improve the study of Introduction 2 disease processes when structural or spatial modifications are involved (D. 1978). The combination of 3D image reconstruction methodologies with staining techniques provides a better understanding on functional information concerning the cellular structures (Roberts, Magee et al. 2012). Hence, there is an urge to develop fully automated approaches for tissue analysis in histological section images, combining the best computational methods to produce an accurate and reliable 3D reconstruction algorithm, enhancing this way the medical study and clinical diagnostic of various diseases. The implementation of these algorithms could provide a better insight into the intricate spatial relations between the studied cell tissues and surrounding tissues. With this in mind, a computational framework was developed in this study, composed by several image processing and registration algorithms and culminating in a 3D tissue reconstruction method. To accomplish this, a previous literature review was performed on the most suited methods to perform each task, as well as histological image notions, crucial to correctly analyse the image datasets tested to validate the developed algorithm, and also to define the relevant tissue structure to be reconstructed. 1.2 - Objectives The present study aimed, primarily, the review and evaluation of currently implemented image processing and registration techniques in the literature, ranging from standard preprocessing methods to complex image registration frameworks, and recently developed algorithms for histological image analysis. With all the concepts and information gathered from the literature review, it was possible to pursue the main goal of this study, which consisted in the development of a completely automatic computational framework, capable of performing accurate histological image alignment and 3D reconstruction of tissues. Therefore, providing detailed volumetric information concerning relevant tissues features, unobtainable through 2D conventional slice analysis. The objectives behind each step developed for the final framework, are explained below. Image pre-processing – first step, developed to accomplish the highest color differentiation between tissues with distinct stains, with several contrast enhancement methods being tested. Image segmentation – aiming to provide the most accurate discrimination between different tissue stains and/or other interesting structures in the histological images, previously pre-processed. Image registration – for the implementation of this step, the great focus was to achieve the most correct slice alignment, thus, mimicking the original disposition of histological tissues in natural conditions, before the tissue preparation procedure. 3 3D surface reconstruction – this final procedure enables the interpolation and surface reconstruction of the previously processed and aligned stack of slices. The main objective with this step is to provide detailed three dimensional insight over commonly studied 2D tissue structures, such as neoplasic tissues. 1.3 - Document Structure The present work is divided in 4 chapters, besides the introduction. A brief description on the contents and subjects addressed in each one is provided below. Chapter 2 - Literature Review: On this chapter, an introduction to histology and its importance as a way to assess cellular responses to pathogens or treatments at a microscopic level is addressed, as well as the standard histological tissue sample preparation, for microscopic observation. This last topic is also discussed, including standard and recent techniques to perform microscopic observation on histological tissue sections. In this chapter, it is also presented a review on both standard and recently developed image processing methods, including pre-processing techniques, segmentation methods and registration algorithms, implemented on stained histological images. Chapter 3 – Methodology: This chapter describes the implemented methodology in this project, approaching first the pre-processing techniques tested, followed by image registration frameworks and culminating in the final workflow. The histological image dataset used to test the developed computational framework is also described in this chapter. Chapter 4 – Results and Discussion: In this chapter, all the results obtained either from pre-processing methods, segmentation or registration methodologies are presented and systematically discussed. Chapter 5 – Conclusion: This last chapter, comprehends the final conclusions about the results, obtained through the performed study, in addition to future development perspectives. 1.4 - Principal Contributions The principal contributions provided by the present work can be subdivided into two domains: the literature review and the developed computational framework. The literature review presented represents an introduction to researchers or developers interested and unfamiliar with histological tissue image processing methods, to fundamental Introduction 4 histology concepts as well as computational techniques, ranging from image pre-processing to image registration algorithms, suited for histological tissue samples. Furthermore, it also presents current state of the art methodologies for image registration, implemented in histological image datasets, which obtained successful slice alignment results. The pre-processing preliminary evaluation on color space transformation methods, best suited to enhance color contrast in stained tissue images, is complete and the resultant image examples presented in this document provide a great insight on these simple algorithms to increase RGB color contrast, not only applicable to histological slices. The CLAHE-red technique was conceived for the present work, and it was proven to be the best pre-processing method to enhance stain contrast in H & E histological images, originating accurate tissue segmentation results, with simple clustering methods (kmeans). The stain deconvolution algorithm, despite being based on previous works on the literature, the computational framework that enabled automatic stain discrimination was conceived and developed in this project. This was accomplished using simple techniques and the results obtained were consistent for most of the tested histology slices. The entire computational framework proposed in this project enables the reconstruction of a three dimensional volume based on real histological tissue structures. It was achieved through the implementation of computationally cheap algorithms, and were obtained highly detailed volumetric representation, not only of differently stained tissues but also, in some cases, the accurate reconstruction of neoplasic tissues, present in the considered histological image datasets. 5 Chapter 2 Literature Review This chapter presents the essential concepts required to understand the topics under study as well as all the research done so far on the subject. Firstly, it will be presented an overview on the histological concepts concerning the current laboratory approach for acquisition of samples, as well as the relevant features of the different types of cellular tissues, since these images from the tissue samples represent the case study of this project. A study on the most suitable computational methods to process and extract information from the histological images is reviewed and analyzed further in this chapter. This literature review culminates in the presentation and analysis of the most accurate 3D reconstruction algorithms for biological images, regarding the future reconstruction of certain relevant portions of cellular tissues in the histological samples. 2.1. Histology Histology is the science that is devoted to study the detailed morphology of cells and tissues concerning the way in which these constitute the different organs in the body, at a microscopic level. The methods implemented by histologists require the study of living cells outside the conditions in which their development is natural, imposing a controlled environment (Junqueira and Carneiro 1987). Histological studies provide an important help in the understanding of some complex pathophysiological processes concerning diseases at the cellular scale. Since these studies are also fundamental to evaluate the performance of new therapies and drug agents, they are considered the gold standard for assessing the natural response of a cellular tissue in face of a Literature Review 12 highlighting most of the significant cellular structures in the tissues. This technique stains the cell nucleus and other acidic structures in blue or violet (haematoxylin) and in pink the cytoplasm and collagen (eosin). Other techniques enable the contrast of more specific tissue structures or organisms by recurring to histochemical reactions, and also, in the case of immunohistochemical stains, to assess the presence or absence of a certain protein. These methods are performed only when the H & E staining fails to accentuate the contrast of the studied structure in the tissue (Junqueira and Carneiro 1987, Randell, Ruddle et al. 2012). Representative images of the visual appearance under the microscope of the referred staining methods applied on cellular tissues are presented in Figure 7. Figure 7. Histology picture of a set of cells lining a duct stained with H & E (on the left), and a histology image stained with immunohistochemical techniques to enhance, in red, the presence of the protein actin in the cells (on the right). Adapted from Leeds University Histology Guide (Michelle Peckham 2003). In the frozen fixation, the tissues to be analysed are rapidly frozen (replacing the chemical fixation and inclusion steps in the previous protocol), and then stained with H & E technique. Despite producing lower quality slides this method acquires slices in a shorter time, ideal when is required a fast examination of the tissue (Randell, Ruddle et al. 2012). 2.1.3. Microscopy and Histological Sample Observation After the preparation of histological samples (more details in section 2.1.2) the microscopic cellular structures present in them are observed under the microscope. In this section, the prevalent types of microscopy implemented to visualize and analyze those tissue slices as well as innovative methods to perform the observation and diagnostic of histological images are addressed. There are two major types of microscopy devices, the light or optical microscopes and the electronic microscopes. The most critical factor concerning a microscope is its resolution power or resolution limit, which is measured by the minimum distance between two particles in the image (Junqueira and Carneiro 1987, Randell, Ruddle et al. 2012). Functional and operational details concerning different microscopes of both groups are explained further. 2.1. Histology 13 The conventional light or optical microscope exhibits images of the stained tissues through illumination, which transverses the sample, generated by a light source. It is composed by both mechanical and optical parts and has a limit resolution of 0.2 µm. The optical part comprises three sets of lenses, namely the condenser, the objectives and ocular lenses. The first condenses the light from the source to the histological sample, the objectives collect the light that crossed the sample and projects an augmented version of the received image, ranging the magnification from 2.5x to 40x, into the ocular lens also contributing for the final magnification in a factor of 10. The final magnification is then, the product of both objective and ocular magnification. However, by convention the ocular magnification factor is not included in image descriptions. Besides the normal light microscope, optical microscopy also comprises other two major types of microscopes, the confocal and the fluorescence microscopes (Junqueira and Carneiro 1987, Randell, Ruddle et al. 2012). Confocal microscopes allow the focusing of thinner sections in the image, avoiding the observation of overlapping planes of the tissue, fact that degrades and reduces the image’s definition. In order to perform this specific focus, the light beam that crosses the histological sample is narrow and the tissue’s image must transverse a small orifice. Consequently, this setup only allows the focussed plane of the original image to reach the detector, blocking all other consecutive planes. Since only a thin section is focussed at a time it is possible the three dimensional (3D) reconstruction by gathering all the planes of the analyzed tissue, through a computational algorithm (application later explored in the following sections) (Junqueira and Carneiro 1987). In fluorescence microscopy, the analyzed samples are lighted by a mercury light source and, by recurring to certain filters the wave-length of the projected light can be regulated. Certain biological structures present in the tissue sample have affinity to fluorescent substances that when excited by the projected light they answer by emitting light in specific wave-length. Through the application of this technique certain biological components exhibit bright colors in the observed image, being highlighted from the surroundings (Junqueira and Carneiro 1987). Electronic microscopy is based on the interaction between electrons and the tissues present in the sample to be analyzed. Considering that light microscopes have a limit resolution of 0.2 µm, electron microscopy represents a more accurate solution, offering a more detailed image of smaller components in the studied tissue with a limit resolution of approximately 3 nm. Nowadays, exist two types of electron microscopes, transmission and scanning electron microscopes (Junqueira and Carneiro 1987). Transmission microscopes possess a resolving power of approximately 3 nm, thus allowing the detailed observation of isolated biomolecules or particles 400 thousand times magnified. For entire tissue samples, the magnifying power is reduced to 120 thousand times, still a high resolution when compared with the optical microscope. The operating mode of this microscope is based on the detour of electrons when in contact with magnetic fields analogous to lens’ Literature Review 14 light reflection in the optical microscope. The electron beam is produced upon heating a tungsten cathode, and, due to a voltage potential between the latter and the anode, the electrons are accelerated and transverse in high speed the microscope tube. In the tube, the beam is condensed through an electromagnetic lens (coils) and interacts with the tissue sample, traversing it and being consecutively amplified by a sequence of magnifying lenses. In the end, the electrons reach a detector (fluorescence plate) and imprint a black and white image of the analyzed sample. The printed grayscale is done according to the amount of electrons that crossed the microscope’s column and so, the tissue sample. Darker spots are electron-dense areas, meaning that more electrons traversed the tissue unaltered, not encountering any structure (Junqueira and Carneiro 1987). Scanning microscopes acquire almost 3D images from the surface of tissues and cells in the analyzed sample. To perform this, the tissue is covered with a metallic coating, and a narrow electron beam is directed to the sample going through the entire surface of the tissue, without traversing it, in opposition to the transmission microscopes. The emitted electrons reflect on the surface and are collected by a detector, amplifying them and, with the intervention of other electronic components, a signal is produced in the form of a black and white image, similar to the transmission microscope. The images produced by this electron microscopy equipment can be consulted in a monitor or stored (Junqueira and Carneiro 1987). Examples of biological images collected from some of the previously referred types of microscopes are depicted in Figure 8. Figure 8. Comparative images acquired from optical microscopy (a) and electron microscopy (b). Image (a) is a light micrograph of Harderian gland from a neonate Alligator mississipiensis stained with Methyl Green-Pyronin Y (bv-blood vessel; hg-Harderian gland and ln-lymphatic node). Adapted from (Rehorek and Smith 2007). Image (b) is a micrograph of a section of mouse liver stained in a saturated solution of uranyl acetate (m-mitochondria and p-highly dense RNA particles). Adapted from (Watson 1958). 2.1. Histology 15 Recent alternatives to the microscopic current approaches have been developed. In order to counter the extensive time dispensed in learning and accustoming to the microscope usage and considering the decrease in use of these devices in medical schools (Randell, Ruddle et al. 2013), these institutions have been using virtual slides (Histopathology slides scanned and stored as digital images), for teaching purposes. These slides allow a greater interaction between the students and the relevant morphological features present in the visualized tissue (Kumar, Velan et al. 2004). In Leeds University, Randell R. and collaborators have developed a virtual reality microscope that consists of a wall-sized high-definition display (Powerwall) capable of rendering gigapixel virtual slides in real time (Figure 9). This system provides a five times greater slide area than conventional microscopes with equivalent magnification, and since it has a wall-size is better suited for group interpretation. This novel approach enables students to cooperatively interpret the displayed images, showing a more interactive apprenticeship. A complementary study, performed by Treanor D. and co-workers, aimed to verify this new solution for virtual slides analysis as a viable replacement for conventional microscopy in the histopathologists’ investigation and diagnostic routine. In fact, the diagnostic made by consulting virtual slices takes 60% longer, mainly due to the considerable amount of time spent to navigate across the entire image in the small display size, provided by common computer monitors and inadequate user interfaces. With this in mind, the aim of this study was to assess if by increasing the display size, using the Powerwall, the diagnostic would reach similar speed when comparing to conventional microscopy. The performed test in this study involved a simple diagnosis, finding small objects in the image, a decision about a lymph node and score a tissue microarray. By using the virtual microscope, histopathologists performed clinical diagnostics and all the other assigned tasks in similar times as when using a conventional microscope (Treanor, Jordan‐Owers et al. 2009, Randell, Hutchins et al. 2012, Randell, Ruddle et al. 2013). Figure 9. Leeds University wall-sized virtual microscope. Adapted from (Randell, Hutchins et al. 2012). In the University of South Carolina School of Medicine, according to (Blake, Lavoie et al. 2003), the transition and implementation of virtual slides and virtual microscopes for teaching Literature Review 16 purposes was performed. The histological slides were scanned and viewed up to a 400x magnification recurring to the MrSID viewer (wavelet-based multiresolution seamless image database, property of LizardTech (Hovanes, Deal et al. 1999)) and the computer as a virtual microscope. The stated approach possesses useful features, including effective microscope and telescope functions providing greater versatility for tissue sample study and increased speed in localizing the structures of interest, when compared to the conventional microscope. In light of the stated, digital pathology promises interesting advantages, both in terms of efficiency and safety considering conventional microscopy procedures (Randell, Ruddle et al. 2012). Potential advantages associated to a digital system reside in the possibility to alert histopathologists about the presence of new slides or cases to be analyzed (similar to the workflow in radiology diagnostics) and provide an easier cooperation between technicians when investigating a particular case. The latter is extremely important in the workflow of specialists, since this digital method allows a faster and safer way to share microscopic visualizations of tissue samples with other specialists, from other labs and also countries, to obtain second opinions, an extremely important procedure to ensure a flawless diagnostic. With the digital procedure, slides can be simultaneously sent to several histopathologists and, since the physical transportation of those slides is inexistent, there is a reduced risk of losing or mixing them, thus avoiding an erroneous diagnostic (Della Mea, Demichelis et al. 2006, Gilbertson, Ho et al. 2006, Nakhleh 2008). Also with the purpose of providing a fast and accurate second opinion to doctors and histopathologists, several computational methods are being developed to process and analyze digital tissue images. Some of those methods are introduced in following sections. 2.2. Image Processing 17 2.2. Image Processing The first approach in order to acquire visual features and information from images, in the particular case of this study, from histological images, involves some computational strategies constituting the image processing procedure. The beginnings of image processing trace back to the middle of the 20th century, when it started to be applied to improve microscope image’s quality, basically through frequency filtering (signal-to-noise ratio, contrast and image restoration methods). Real developments were made since then, and the analog image processing was replaced by digital image processing with the advent of powerful computers capable of applying sophisticated algorithms to large images in an acceptable amount of time (Bonnet 2004). Since the visual interpretation is the core of most medical diagnostic procedures and the final diagnostic decision, for cancer and other diseases, is based on tissue examination, medicine represents a large application field for image processing and analysis algorithms (Bengtsson 2003). However, visual interpretation and evaluation present several weaknesses. For pathologists, it is a time-consuming, cumbersome and tedious process to analyze a large number of tissue samples in practice, thus requiring a long time and intensive manual labor to produce viable results. Besides from this problem, visual evaluations can, in many cases, be subject to unacceptable inter and intra-reviewer discrepancies (20% discrepancy between central and institutional reviewers, as reported by Teot L.A. et al. in (Teot, Sposto et al. 2007)), due to the sampling bias, confirming that it represents an error-prone method (Bengtsson 2003, Kong, Sertel et al. 2009, Sertel, Kong et al. 2009). To overcome the stated weaknesses, established in the currently used visual evaluation process, allied to the fact that digital images are growing in popularity, computational methods resorting to automated image processing and analysis algorithms are being developed (Bengtsson 2003, Kong, Sertel et al. 2009). The automatic processing and analysis of tissue images provides reliable data, accelerates data acquisition process and by allowing digital image management it can replace other evaluation methods, more expensive and impossible to execute (Cisneros, Cordero et al. 2011). Automated systems can exclude human factor mistakes and increase the speed of diagnostic processes. These systems represent an important asset considering that the amount of experienced specialists that conduct a correct histological analysis is reduced or concentrated in big medical centers. Therefore, this leads to an accumulation of cases poorly or misdiagnosed, conducting to incorrect untimely treatments and ultimately resulting in disablement or death. Considering the abovementioned it is clear the need for automated processes concerning morphology diagnostics in medicine, to improve the diagnostic accuracy and compensate the scarce number of specialists (Nedzved, Belotserkovsky et al. 2005). Literature Review 18 The key challenges in histological image computational analysis are automated cellular segmentation and classification in tissue images. Nevertheless, due to the complex nature and variety of histological images, it is difficult to develop automatic segmentation methods applicable to any type of those images (Nedzved, Belotserkovsky et al. 2005, Chomphuwiset, Magee et al. 2011). The general procedure for automatic image processing and analysis can be divided in several steps, starting with the acquisition of digital histology images, which can range from diverse resolutions depending on the application and the size of the biological structure in study, on the histological sample. The following step is the image processing to identify the target tissues or biological structures in question, comprising, as a standard framework, image enhancement, image segmentation, feature extraction and implementation of machine learning algorithms (Caicedo 2009). To perform each of the previously referred stages a wide variety of computational methods can be implemented, according to different purposes (for example, automation of mass screening of histological specimens or quantitative analysis of a significant structure in the tissue) (Bengtsson 2003). The image processing and analysis pipeline that is going to be produced in this work consists of three particular steps, the image pre-processing, image segmentation and 3D reconstruction (in section 2.3). A state of the art on methods for both these steps is presented in more detail in further sections. 2.2.1. Image Pre-Processing Although segmentation is the most important step in image processing and analysis, it is unusual to achieve a consistent and useful segmentation using only a single procedure. In order to obtain a successful segmentation, algorithms typically apply a constructed combination of methods, including a wide variety of preprocessing steps (Beare and Lehmann 2006). To process histological images, an initial preprocessing step must be applied to reduce the computational costs through multi-scale image decomposition (Gonzalez and Woods 2008). This initial process produces low resolution images that can be analyzed to locate interesting structures and allow the implementation of other image processing steps only on those structures’ pixels. The preprocessing step is meant also to restore the images, by reducing image noise, low intensity contrast and intensity inhomogeneities present in the histological data. To perform this, methods such as image smoothing, denoise and enhancement can be applied (He, Long et al. 2010). Image smoothing is commonly performed recurring to spatial filtering methods used to remove image high frequency noise. Image denoising methods are implemented to remove image noise produced in image acquisition and compression processes (Aubert and Kornprobst 2006) and image enhancement techniques favour an increase in contrast between the regions of interest and the background, being the adaptive filters the most commonly employed methods (Gonzalez and Woods 2008). 2.2. Image Processing 19 Considering the amount of manual labor involved in tissue samples preparation (section 2.1.2), this process tends to introduce certain types of artifacts that require proper image preprocessing techniques to be countered. The majority of artifacts found in histological images are based on orientation differences found in the sections mounted in glass slides, variable luminance gradient (depending on the slide region where the tissue observed), nontissue noise produced by dust or bubbles and staining variations (variable tissue thickness and stain concentrations originate color variations in the histological stained structures). Therefore, image pre-processing techniques are applied to deal with acquisition artifacts and defective histology sections (Mosaliganti, Pan et al. 2006). Some techniques are specially designed to deal with histological artifacts present in the digital images, namely the defective section exclusion and principal component analysis (PCA) alignment. In order to improve the 3D reconstruction robustness defective sections have to be identified and removed from the registration process. Since all images are acquired with the same magnification, tissue sizes in consecutive images should not suffer significant variations. Thus, when a large variation is verified it is probably due to broken or defective sections. To eliminate them the relevant structure areas were computed for each image and plotted against section location, using binary masks (masks containing information about the tissue pixel location. Tissue pixels identified and stored as binary masks). Spikes in this plot are potential defective sections (Mosaliganti, Pan et al. 2006). Principal component analysis alignment is used to estimate tissue orientation, according to prior knowledge of typical structure arrangement, concerning the studied tissue. Since tissue orientations are used to initialize registration methods (section 2.3), by using this technique, the likelihood of converging to a more reliable global solution is increased (Mosaliganti, Pan et al. 2006). The staining conditions verified in different histological slices suffer considerable variations (Figure 10). Therefore, after the being digitized, images present considerable color ranges. To normalize color distributions present across the slices, histogram equalization represents a viable solution (Sertel, Catalyurek et al. 2009). Histogram equalization is a wellknown and widely used image enhancement technique, due to its simplicity, high performance in almost all types of images. This technique is performed by remapping of grey-levels in an image based on a probability distribution of the input grey-levels, stretching the dynamic range of the image histogram. Thus, resulting in overall image contrast enhancement. The drawbacks of this method are noticeable in images with high and low mean brightness. The result is a significant change in the image outlook, whereas the purpose was only to enhance the contrast. Histogram equalization is best suited to enhance the edges between different structures, but, in return, reduces local details within those structures, producing over enhancement and saturation artefacts (Kaur, Kaur et al. 2011). The histogram equalization procedure is based on the assumption that the processed image presents uniform image quality in all regions, and therefore, one single grayscale mapping provides similar contrast enhancement throughout all Literature Review 20 these regions. But, when distributions of grayscale intensities are variable according to different regions in the image, the previous assumption is invalid. Facing this, an adaptive histogram equalization technique, capable of determining the mapping for each pixel based on its local grayscale distribution (surrounding pixels), could significantly outperform the standard method. However, when grayscale distribution is highly localized, full histogram equalization might not be desirable to transform very low contrast images. By limiting the contrast allowed through histogram equalization, the possibility of two very close grayscales being mapped to significantly different grayscales, i.e. high slope segments present in the grayscale mapping curve, can be avoided. Combining the contrast limiting approach with the previously mentioned adaptive histogram equalization method, the result is referred to as Contrast Limited Adaptive Histogram Equalization (CLAHE) (Reza 2004). Figure 10. Sequence of histological images from MKI (Mitosis-karyorrhexis index) cells, varying in color information due to staining differences. Adapted from (Sertel, Catalyurek et al. 2009). For medical images, color enhancement represents a valuable tool to aid in visualization, detection and segmentation of specific tissue structures. An effective approach to increase color contrast on an image and maintain its hue (the pure color) is to transform the RGB (R - Red, G - Green, B - Blue) color into HSV (Hue, Saturation, Value) color space, modifying only the saturation (luminance) value in the image’s pixels. The transformation is given by: 𝐻=𝑐𝑜𝑠−1×( 1 2[(𝑅−𝐺)+(𝑅−𝐵)] √(𝑅−𝐺)2+(𝑅−𝐵)(𝐺−𝐵)) (1) 𝑆=1− 3 𝑅+𝐺+𝐵[min⁡(𝑅,𝐺,𝐵)] (2) 𝑉=1 3(𝑅+𝐺+𝐵) (3) The luminance manipulation, through grey-level enhancement processes, without affecting the other two components in HSV color is possible due to the lack of correlation between these components. The usual process starts by performing the color transform to 2.2. Image Processing 21 convert the image in HSV color where the luminance or color saturation can be modified, disregarding the other components. Then, the reverse transform back to RGB color is applied in order to ascertain the effects of the produced modification (Bautista and Yagi 2010). Hukkanen J. and co-.workers implemented a pre-processing method to improve the efficiency of nuclei segmentation (more information in section 2.2.2.) in histological images. The pre-processing method performs the conversion of H & E stained histological images originally in RGB color space into CIE L*a*b color space. The L component, the luminosity component, is then denoted as a grey-level image, which is further processed to obtain the segmentation. The L*a*b color space consists of a luminosity layer “L*”, a chromaticity layer “a*” (indicating the color location in the red-green axis) and a chromaticity layer “b*” (indicating the color location in the blue-yellow axis)(Hukkanen, Hategan et al. 2010). In (Tabesh, Teverovskiy et al. 2007) is presented a study concerning image features for cancer diagnosis and histological grading of prostate images. The features representing color, texture and morphological details were combined in a supervised learning framework. The first stage in this framework involved pre-processing techniques, including background removal and image histogram matching to a reference image. The background was identified and then removed from the analysis through color tissue image transformation, from RGB color space into YCbCr color space (Gonzalez and Woods 2008), and posterior thresholding (section 2.2.2.) of the luminance (Y) component with a global empirically determined threshold. After this, the binary mask containing the tissues of interest is refined via closing and opening operations to fill gaps between tissue structures and remove small artefacts from the image. A convex hull operation (Gonzalez and Woods 2008) is then, applied to ensure the integration of lumens as tissue of interest, avoiding its exclusion from the binary mask. The second pre-processing step, implemented in this study, consisted of an histogram matching between the analysed histological image and a reference image, through the transformation⁡𝐹𝑟−1[𝐹𝑖(𝑥)], where 𝑥 is the pixel value in each of the red, blue and green channels. 𝐹𝑟 and 𝐹𝑖 are, the cumulative distribution function of pixel values for the input and reference images, respectively. Histogram matching is performed to mitigate color variations produced by staining and illumination conditions, which can affect segmentation efficiency. The Stain Deconvolution technique is a pre-processing method, based on color deconvolution, which aims to deconvolve the applied stains on a certain RGB color image, to generate separate images, where each grayscale image shows the distribution of a single stain. This algorithm assumes that the chemical stains implemented to dye the tissue slides follows the Beer-Lambert Law of absorption, which provides a logarithmic relationship between the original RGB color channels and a stain matrix. This complex methodology is presented by both Literature Review 28 based one, and it was capable of reconstructing selected tumor invasion fronts from a considerable extend of histological serial sections (90-500 slices). Through this technique it was achieved an extremely detailed 3D reconstruction of invasion of solid tumors, as visualized in Figure 13. Figure 13. Two views of 3D reconstruction of uterine cervix carcinoma tumor invasion fronts, from different histological specimens. Adapted from (Braumann, Kuska et al. 2005). An extension of the previous study was carried out by Ryk J. and co-workers, where the Large Image Microscope Array (LIMA), a vibratome capable of sectioning tissues to 40 mm thickness slices, is used to section entire organs, particularly lungs, into slice images in order to establish direct correlation between lung pathology and Computer Tomography (CT) images. A camera and a stereomicroscope mounted on the vibratome scan the entire surface area of the tissue, controlled by a custom software, responsible for the entire process automation. The alignment is accomplished through the combination of a custom code and the Insight Segmentation and Registration Toolkit (ITK), resulting in the registration of high magnification and resolution pathology images with the corresponding CT images (de Ryk, Namati et al. 2004). Another technique developed by (Rehorek and Smith 2007) aims for the generation of 3D images from specific microanatomical structures contained in tissue sections. The presented 3D reconstruction technique not only allows the concurrent visualization of multiple structures or tissues but also enables the analysis of spatial topography from the histological sections. This technique involves, foremost, the visual identification and manual delimitation of the region of interest’s position in a 2D section digital image. All the marked digital images were then aligned to previous slices based on fiducial landmarks (reference structures in the image), accomplished by rendering the second of two adjacent tissue sections through manual rotation of the image on top of the other slice. A semiautomatic method was developed by (Kurien, Boyce et al. 2005) to produce threedimensional reconstructions of invasive breast carcinoma, recurring to common laboratory equipment to evaluate the spatial arrangement of parenchymal cells. The tumour studied in this work was stained immunohistochemically to reconstruct two 3D images, one for normal parenchymal cells and the other for malignant. The digital histological sections were acquired using a microscope, a scanner and a camera connected to a computer, and the alignment was 2.3. Image Registration and 3D Reconstruction 29 performed using a semiautomatic method, allowing manual interaction through a graphical interface. The referred method uses cross correlation coefficient as an integrity slice fit measure and an automatic process based on the Fibonacci search algorithm (Ramaprabha, Balaji et al. 2012) to achieve an automatic alignment. Finally, the resultant reconstructed volume was obtained using maximum, minimum point projection and back to front opacity blending. The results obtained through this method were distinct and accurate 3D reconstructed models, contributing to a more comprehensive and explicit visualization (high resolution) of spatial arrangement of normal and malignant parenchymal tissues and their relation to the surrounding tissues (Figure 14). This approach provided an insight on invasive breast carcinoma proliferation, unobtainable through conventional visualization of 2D histological sections. Figure 14. 3D reconstruction of invasive breast carcinoma immunohistochemically stained, illustrating the spatial arrangement of the different parenchymal tissues. Adapted from (Kurien, Boyce et al. 2005). In (Ourselin, Roche et al. 2001) it is emphasized the need of a robust method to perform the alignment of histological sections for 3D reconstruction. In this work, it is proposed an intensity-based method to register the image slices, first by using a block matching strategy (Jain 1981), allowing the computation of local displacements between image slices, and then, a rigid transformation, estimated by those local measures. The entire process is fully automated and integrated into a multi-scale framework in order to improve accuracy and computation times. The results obtained experimentally highlight the capacity of the proposed algorithm to reach sub-pixel accuracy while being able to compensate large displacements between slices. Conventional 3D histopathology is limited by low resolution, time and difficulty with acquiring a large number of images with a microscope, the absence of a fully integrated system for 3D reconstruction (Namati, De Ryk et al. 2007) and, in the case of manually guided 3D reconstruction, the time required for the whole process (Petrie, Flynn et al. 2002). Concerning those limitations, Roberts N. and collaborators developed a 3D histopathology software using Literature Review 30 automated virtual slide scanners to produce high-resolution digital images and 3D reconstructions of the cellular tissue. The software is applicable to any type of tissue, requires minimal manual intervention, once the slides are prepared, and reveals reasonably robustness over a wide variety of data or type of application. The system digitises automatically the virtual slides, which subsequently communicates with the image serving software to align the image and produce the visualization. To perform the alignment it uses a high-resolution registration followed by a multi-level registration method, whereby the user is able to manually select, zoom and re-register the area of interest (Roberts, Magee et al. 2012). The registration algorithm proposed by the author starts by performing a sequential slice-to-slice image based registration, a multi-stage method based on extension of phase correlation (De Castro and Morandi 1987), which consists of a first rigid alignment ignoring scale that serves as input to a non-rigid registration method that divides the input image into equally spaced square patches, individually aligned through the phase correlation. A non-rigid B-spline based transform is estimated through a least squares error minimising method and applied to the whole image to approximate a set of points from each patch. The B-spline is applied at multiple increasing resolutions and B-spline grid sizes. This registration method uses as reference a representative central virtual slide (generally the slice containing the largest portion of the tissue to be studied). The subsequent images are aligned to their neighbours and then concatenated to form a 3D volumetric dataset for the user to visualize, after the completion of rigid and non-rigid transforms. As previously referred in this study, it was also developed a user-interface to promote the interaction and user selection of sub-areas of interest in the image to re-register at higher resolution. This was performed using the explained non-rigid multiple-level resolution method, providing this way sub-cellular accuracy reconstructions. After the 3D volume reconstruction the developed software enables interactive segmentation of the volumetric structures by implementing several techniques, including a manual method called “color example thresholding” in which the user defines a threshold and selects a certain color in the image and the software annotates all the pixels sharing RGB color with the selected one and the considered threshold. Another approach implements region growing algorithms (Efford 2000) to segment spatially connected similar pixels (more information on section 2.2.2), using the user selected points as seeds and the threshold as color similarity criteria. In the end, it was implemented an iso-surfacing using marching cubes (Lorensen and Cline 1987) and mesh decimation (code from VTK (Schroeder, Martin et al. 2003)) to render the segmented volume and allow the visualization of distinct parts of the volume separately (Magee, Treanor et al. 2008, Roberts, Magee et al. 2012). In (Chakravarty, Bedell et al. 2008) is described a method for 3D reconstruction of two dimensional histological sections from mouse brain to create volumetric data, a fundamental step in the analysis of ex-vivo data to validate in-vivo imaging techniques. To achieve the latter, the reconstructed volume is directly mapped into in-vivo anatomical MRI volumes, thus enabling accurate validation of recent imaging technologies, as well as the integration of histological cellular or molecular data in in-vivo structural and functional data. The procedure 2.3. Image Registration and 3D Reconstruction 31 described in this work is an extension of a previous study by Chakravarty et al. (Chakravarty, Bertrand et al. 2006) performed for the reconstruction of human basal ganglia and thalamus, where the histological sections suffered first manual alignment and color adjustment, and then non-linear correction techniques to enhance the image alignment. In (Chakravarty, Bedell et al. 2008) it was used a digital blockface to perform image registration and a volume derived from the average of ten MRI volumes of mice (Lau, Lerch et al. 2008) to serve as volumetric reference. The method stated by the author starts with the alignment of histological sections to the corresponding blockface image. First, a slice-to-slice alignment of the blockface data was executed and a transformation mapping each slice to the next, converging this process to the center slice (reference), was estimated. Then, for the histological dataset it was applied and estimated a linear transformation to map all the slices into the corresponding slice in the blockface reconstructed volume. Both two dimensional transformations were adapted from the linear registration technique proposed in (Collins, Neelin et al. 1994). Posteriorly a nonlinear morphological correction was implemented, consisting of a slice-to-slice warping to correct morphological inconsistencies between slices. To perform the estimation of these transformations the ANIMAL algorithm (an iterative algorithm that estimates a 3D deformation field on a lattice of nodes when a source volume is matched to a target volume) was used (Collins and Evans 1997). Then, in order to counter intensity inhomogeneities between registered slices, possibly produced by staining densities and slice thickness irregularities that can challenge the correct visualization of the images, a nonlinear intensity correction was performed. This procedure first applies the nonlinear transformations estimated in the morphological correction to consecutive series of four slices, to grant that all of them present morphological consistency. After this, each slice was portioned into equally sized square patches, assuming that sufficiently small areas would enable the estimation of a first order polynomial scaling factor to match the joint histograms of each square by using a least-trimmed squares polynomial estimation (Prima, Ayache et al. 2001). Once performed this process for each square patch, a grid containing the scaling factors previously estimated was interpolated to match the resolution of the histological image, and then, the produced interpolation field was multiplied by all the slices to enhance the grey-level consistency. Finally, the histological reconstructed volume was directly warped to the stated MRI template volume reference. The results obtained through this routine demonstrate an effective alignment of the histological volume with the MRI template. These results enhance the fact that the combination of tissue preparation techniques and automated image processing allows for a more comprehensive, multi-modal evaluation of pathology or therapeutic intervention effectiveness in rodent models of CNS disease. The results obtained through the stated technique are presented in Figure 15 (Chakravarty, Bedell et al. 2008). Literature Review 32 Figure 15. Volumetric results from the reconstruction of serial histological slices acquired by Chakravarty M. and collaborators. From left to the right: Blockface reconstructed, original data stacked without any correction, reconstruction model after slice-to-blockface linear alignment, resulting reconstruction from nonlinear morphological correction and reconstruction after nonlinear intensity correction. Adapted from (Chakravarty, Bedell et al. 2008). Johnson J. and co-workers presented a hybrid landmark/intensity-based deformable registration algorithm. This algorithm applies an iterative process by producing accurate correspondences between image structures near landmark locations and elements separate from them by matching corresponding landmarks and image intensities, respectively. First are registered landmarks disregarding intensity and then intensity differences between slices are minimized (Johnson and Christensen 2002). Despite the accurate results, this technique was not applied to reconstruct 3D images from 2D slices and also revealed that when applied to unclassified data based on landmarks does not lead to reliable registrations. Another method for automatic registration of histology sections was presented by Arganda C. and collaborators, consisting in a technique for non-deformable registration recurring to Sobel transforms and segmented contours. The goal of this project was to accurately align tissue sections in volumetric data and also to detect and render relevant structures in 3D. The algorithm proposed accomplished this by finding the best suited rigid body transformation (translation and rotation) of the images being registered (applied globally), through maximization of a matching function based on image correlation (applied locally on specific locations revealed by segmentation methods). A multiresolution pyramidal approach was then implemented, reaching the best registration transformation in increasing resolution stages (Arganda-Carreras, Fernandez-Gonzalez et al. 2004). In (Mosaliganti, Pan et al. 2006) was presented a mutual information based registration approach, having as basis the maximization of mutual information (MI) (Maes, Collignon et al. 1997), an effective similarity measure to register multi-modal images when the image intensities are not linearly correlated. The proposed technique includes four stages, including the transform, metric, optimizer and interpolator phase, in order to register a consecutive slice into a stationary image. The transform stage is modelled, as previously referred approaches, 2.3. Image Registration and 3D Reconstruction 33 to a rigid 2D transform, which allows rotation and translation. To perform the optimizer and interpolation stage were incorporated a regular step gradient optimizer (Maes, Vandermeulen et al. 1999) and a bi-linear interpolation (Maes, Collignon et al. 1997), respectively. The registration was acquired through MI optimization multiresolution strategies, as proved in previous studies that can provide similar robustness to direct registration combined with an increased computational speed when compared to other multiresolution techniques (Studholme, Hill et al. 1996, Maes, Vandermeulen et al. 1999). This approach was performed applying 3-level image pyramids (Figure 16), with image magnifications of 10x, 20x and 50x. Starting with the transform obtained for the lower magnification images are scaled and employed to initialize the next higher magnification, repeating this process until the highest resolution (the improvement in accuracy for resolutions above 50x do not compensate the computational costs). After this, a two-level optimization (Figure 17) was executed to achieve a higher MI in the reconstructed model. This process introduces a restricted translation and rotation around the PCA initialization (described in 2.2.1), altering the converged solutions. This step is repeated until no improvement is verified in the mutual information model. The presented study developed a multi-resolution MI based registration algorithm combined with a novel optimization strategy that allows not only a reduction of manual intervention in the registration process but also a higher chance of obtaining a converged global solution (Mosaliganti, Pan et al. 2006). Figure 16. Schematic representation of the ITK-based registration framework. The transforms pass from a lower to a higher resolution based on the 3-level image pyramids shown. Adapted from (Mosaliganti, Pan et al. 2006). Figure 17. Representation of the regular step gradient descent (A,C) and the two-level optimizer (B,D) on 2 slides from placenta image dataset (PCA was applied initially as pre-processing), used in (Mosaliganti, Pan et al. 2006). Literature Review 34 2.4. Key Issues 35 2.4. Key Issues The present section summarizes the contents addressed in the literature review, highlighting the fundamental concepts to be retained.  Histological techniques are considered the gold standard for assessing the natural response of a cellular tissue in face of a pathology or therapeutic intervention.  The existent fundamental types of tissue are epithelial (covers body surfaces), connective (provides support and connection to all other tissues in the body), muscular (designed for contraction) and nervous (ensures the communication between sensorial cells, CNS and muscles).  In order to visualize histological samples and their tissue structures under the microscope, a previous preparation of the tissue comprising fixation, inclusion and staining techniques must be performed.  Digital histological slices provide numerous advantages when compared to conventional microscopy, namely in the possibility to alert histopathologists about the presence of new slides for analysis and provide easier cooperation between technicians when investigating a particular case (i.e. a second opinion).  Automated computational systems, comprising segmentation and 3D reconstruction algorithms can providing a fast and accurate second opinion to doctors and histopathologists.  Computational systems can also reduce the workload of histopathologists, who perform a rather time consuming and laborious task, subjected to high inter and intra-reader variability.  Image pre-processing methods are applied to the histological images to reduce the visual noise and enhance the contrast between the interest structures and the remainder tissue.  Segmentation is the most relevant process in image processing and analysis systems, enabling the distinction and identification of some particular biological tissue from the remainder components in the input image. Literature Review 36  Three-dimensional reconstruction of tissue samples at a microscopic resolution reveals significant potential to improve the study of disease processes when structural or spatial modifications are involved.  Image registration is considered an optimization process that aligns two or more images, aiming to find the image transform responsible for the best alignment or closest similarity of structures of interest between consecutive images. 37 Chapter 3 Methodology In this chapter, the methodology developed and all the algorithms constituting it will be thoroughly explored, as well as all the image datasets containing animal cell tissue, tested to provide the necessary validation for the proposed framework. Firstly, the experimental database analyzed in this study, cordially supplied by the Pathology Laboratory of the Institute of Biomedical Sciences Abel Salazar, is going to be addressed, including the description of all three different tissue datasets, the properties of the digital image acquisition device and the relevant tissue regions to be processed and highlighted by the proposed methodology. The second section of the methodology presents an exhaustive description and explanation of all the computational methods developed and employed to perform image processing and registration on the dataset, culminating with the final framework implemented integrating the most suited pre-processing, segmentation, registration and 3D reconstruction of the cell tissues. The presented computational framework was implemented in Matlab R2014a® (Inc., Natick, Massachusetts, United States). 3.1. Dataset The three image datasets studied in this project were acquired and prepared, using the standard method detailed in section 2.1.2., in the Pathology laboratory in Institute of Biomedical Sciences Abel Salazar. In order to obtain the digital images an Olympus scanner program was used, creating an image for each slice in the three studied cases, described further in this section. Some features found in images from the datasets that can undermine the image processing and analysis proposed by the algorithm develop in this work, are also addressed in the present section. Methodology 44 image threshold (section 2.2.2.) for the Cb grayscale image is determined using the Otsu algorithm (Otsu 1975), and then applied to the image to obtain a binary mask of the relevant structure in the image (the background intensity becomes 0). The final mask is obtained after small adjustments with morphological operators, namely image closing and opening (inthis order). The effectiveness of all the presented image pre-processing methods was tested with a kmeans algorithm (section 2.2.2.), based on (Mathworks), upgraded to perform colored image segmentation. This clustering algorithm was implemented in this stage due to its simplicity and low computational cost, enabling the image segmentation in multiple classes, or, in this case, cell tissues, to determine the top-performing pre-processing method in stain color discrimination. The color upgrade was achieved by reshaping the incoming RGB image matrix to a row × column × 3 matrix, thus gathering on each column the pixel intensity information of each color channel image. These pixel intensity triplets were used in the segmentation procedure to discriminate classes. The kmeans was applied for 3 and then for 4 classes, in order to separate the image background from two or three different cell tissues. The general scheme of the first workflow implemented is presented in Figure 22. Pre-Processing Color space transformation and manipulation (All the color space were tested) CLAHE Image Normalization + Contrast enhancement + Mask Segmentation Kmeans: for 3/4 classes Figure 22. Schematic model representing the first workflow implemented in the study. Several contrast enhancement pre-processing techniques were applied to the dataset, and their efficiency was tested with the kmeans algorithm, for 3 and 4 classes. 3.2. Workflow implementation 45 3.2.2. Second Approach – Based on Image Registration The knowledge acquired from the pre-processing study, concerning the best methods to provide an accurate segmentation, led to the registration test, where several registration techniques were approached. To determine the most appropriate, two criteria were taken into account: the DICE score (Sørensen 1948, Cheung 2012) and the computational cost. Four registration algorithms were tested in two different registration models, pairwise and reference slice (more details in section 2.3.). The first algorithm implemented is a feature-based registration method and the second an intensity-based algorithm. The feature registration method only enables rotation and translation to the moving images, since the only type of transformation implemented is similarity. On the other hand, the intensity–based model allows other two types of transformation, the rigid and the affine (section 2.3.). The other two algorithms tested, the Demon algorithm and the B-spline algorithm, which are both capable of performing non-rigid registration. All these methods are explained further in this section. The automatic feature-based registration starts by detecting image features in both images, moving and fixed, mainly through the implementation of the Speed-Up Robust Features (SURF) algorithm (Bay, Ess et al. 2008), an algorithm that searches blob features, regions of the image where several properties (for example, brightness and intensity) remain constant. The detected features are then extracted including their location in the image, through the pixels surrounding the interest point or feature blob (Bay, Ess et al. 2008). The detected and extracted feature regions in the moving and fixed images are matched using parallel hierarchical clustering trees (Muja and Lowe 2012), resulting in a pair of indexes from the matched features. The locations of those matching points are also retrieved and the transformation is performed based on the matched points of both fixed and moving image, resorting to M-estimator Sample Consensus (MSAC) algorithm (Torr and Zisserman 2000) to exclude outliers. The inliers of both sets of matching points are mapped and originate a 2D geometric transform object (Hartley and Zisserman 2003). In the last step, the geometric transformation object generated is applied to the moving image, creating this way the registered image. The scale and angle applied (considering that it is a similarity transform) with the transform object to create the final image are recovered to serve as quality control of the whole procedure, allowing this way the removal of completely distorted images, possible outcomes of the registration process. The stated framework was performed based on (MathWorks). Concerning the intensity-based registration three key components have to be previously defined in order to configure the whole registration process, the optimizer, the metric and the Methodology 46 transformation type (section 2.3.). The metric is, in fact, the image similarity metric responsible for evaluating the registration’s accuracy and the optimizer (Regular Step Gradient Optimizer, whose algorithm is described in (Pennec, Cachier et al. 1999)) defines the procedure for minimization or maximization of this similarity metric (Mean Squares metric implementation in (MathWorks)). This registration method is an iterative process that can be performed in three different transformation types (rigid, similarity and affine) and always requires two images, a fixed and a moving image. The whole registration process starts with an internally determined transformation matrix combined with the transformation type specified, determining the image transformation that is going to be applied to the moving image with bilinear interpolation (MathWorks). After the interpolation the metric compares both transformed moving image and fixed image, computing the metric value and then the optimizer checks for ending conditions to stop the registration process. These conditions can be the maximum number of iterations (defined by the optimizer) or a certain metric value threshold. If the stop condition does not verify the optimizer adjusts the transformation matrix to initiate a new registration cycle. The maximum number of iterations and the optimizer step size the registration can be altered to improve the registration, but always with a greater computational cost (MathWorks). The explained workflow is presented in Figure 23Figure 23. Figure 23. Representation of the intensity-based registration framework implemented. Adapted from (MathWorks). The first non-rigid registration methodology tested was the B-Spline Grid, Image and Point Registration developed by (Kroon 2008), based on the algorithm developed in (Rueckert, Sonoda et al. 1999). This is an intensity based registration technique whose algorithm implements a grid of B-spline control points (section 2.3.) that control the transformation of 3.2. Workflow implementation 47 the moving image over the fixed one. It measures the registration error, through squared pixel distance (Vercauteren, Pennec et al. 2009), a similarity criterion based on the information theory and calculates the amount of information existent in the registered moving image about the fixed one (considering a registration procedure). The B-Spline method applies the Fast Limited Memory Optimizer (Kroon 2009), a Quasi-Newton optimizer, to move the control points, in order to achieve the optimal registration between both images with minimal similarity error. The implemented B-Spline method can also performed rigid and affine transformations (section 2.3.). The Demon algorithm is a non-rigid registration technique faster and rather simpler than the B-Spline. This algorithm was first described by (Thirion 1998) and followed by (Wang, Dong et al. 2005), and the methodology followed in the present study was developed by (Kroon 2008). For each pixel a velocity, or movement, is defined by this method, using the intensity differences and gradient information. The velocity matrix is smoothed by a Gaussian filter and iteratively applied to transform the moving image and register it onto the fixed image. The transformation is optimized by a limit memory BFGS optimizer (Liu and Nocedal 1989) in an iterative and multi-resolution way. The Demon algorithm also performs affine registration (section 2.3.). In the second approach were also tested combinations of the previously described algorithms, starting with a rigid registration technique followed by a more accurate and computational demanding non-rigid registration algorithm, as performed in previous studies (Roberts, Magee et al. 2012). In order to align all the images from the datasets available, two different registration models, the reference slice model and the pairwise model, were implemented for all the abovementioned registration methods, starting from the middle slice (in the datasets) since, generally, it is the section with most tissue (Roberts, Magee et al. 2012). Through the reference model the registration procedure is performed considering only as fixed image the middle slice from the image dataset, thus, being all the slices registered to reference section. The pairwise model is performed in a cascade process starting from the center slice and performing registration in pairs of slices (moving image becomes fixed image in the next alignment) in two directions - until the top slice in the first run and the first slice in the second run. The workflow of the registration approach combined with the best suited pre-processing method is presented in Figure 24. Methodology 48 3.2.3. Final Approach – Based on the complete workflow with the 3D reconstruction In the third and final approach, the most successful pre-processing and registration method, or combined methods, for the tested image datasets were connected and the 3D reconstruction final step was added to the workflow. In the pre-processing stage the previously tested methods were subjected to an efficiency test, but to a different segmentation technique, the stain deconvolution method, explained in more detail further in this section. The method implemented for this first step was the CLAHE but with a difference in the algorithm presented in the section 3.2.1., since only the red channel was subjected to the histogram equalization. This was performed due to the high influence of the red channel in both Hematoxylin and Eosin stain, proving, through the segmentation, to be the best stain discriminant method. The image was also subjected to a FIR filter (filters with finite impulse response) convolution to remove line artifacts or deformations on the analyzed slice (section 3.1.). Image Registration Intensitybased registration Type: Rigid and Affine Featurebased registration Type: Similarity Demon algorithm Type: nonrigid B-spline registration Type: nonrigid Combination of two registration models Pre-Processing 1. YCbCr transformation (with normalization) 2. CLAHE-red Figure 24. Schematic model representing the second workflow implemented in the study. In the first step, the most successful pre-processing method was applied to the image dataset, followed by the image registration step, with several methods being tested. 3.2. Workflow implementation 49 The stain deconvolution technique developed for this study was based on the work carried out by (Ruifrok and Johnston 2001, Chan 2013) and also from (Unpublished 2015). The stain deconvolution is a segmentation method that is capable of deconvolve a stained image (for example with H&E) and generate separate images, each one with the distribution of a single stain (Hematoxylin, Eosin and Background, in the case of the tested datasets). This method assumes that the chemicals employed to stain the tissues slides follow the Beer-Lambert Law of absorption (11): 𝐼=⁡𝐼𝑜𝑒−𝑆𝐷 (11) where 𝐼 is the RGB image, 𝐼𝑜 is the intensity of the source, 𝐷 the stain density map and 𝑆 the stain absorption matrix. This equation grants a logarithmic relation between the original RGB color channels and the stain intensities. The implemented algorithm starts by converting the RGB image 𝐼 into the optical density place (OD) – 𝐼𝑂𝐷 - through the previous equation in the form (12): 𝐼𝑂𝐷 =⁡−log⁡(𝐼 𝐼𝑜) (12) After the OD conversion, the image pixels are projected into the Maxwellian chromaticity plane. In the Maxwellian plane the distance between two points is proportional to their chromatic dissemblance (Maxwell and Zaidi 1993), or within this study, the pixels that belong to the same stain are projected closer than the pixels belonging to different stains. To calculate the pair of corresponding Maxwellian coordinates (𝑚1,𝑚2) for each pixel 𝑝, the following equation (13) was implemented: [𝑚1 𝑚2]=⁡[0.7071 −0.7071 0 −0.4082 −0.4082 0.8165]𝑝 ‖𝑝‖1 (13) The described Maxwellian plane property and the pair of Maxwellian coordinates allows the implementation of an unsupervised classifier, for example kmeans, to assign a label for each pixel according to its stain. With the classification performed, the following step is the Stain Matrix estimation. For this study, it is considered the existence of three stains or classes, the Hematoxylin (class 1), Eosin (class 2) and the Background (class 3). The Maxwellian coordinates of the mean value in each of these classes - 𝜇𝑐=(𝜇𝑚1𝑐,𝜇𝑚2𝑐), obtained through the kmeans classification, is converted back to the OD space through the equations (14) exposed below: 𝑆𝑐,𝑟 =0.7071𝜇𝑚1𝑐 −0.4082𝜇𝑚2𝑐 +0.33 𝑆𝑐,𝑔 =−0.7071𝜇𝑚1𝑐 −0.4082𝜇𝑚2𝑐 +0.33 (14) 𝑆𝑐,𝑏 =1−𝑆𝑐,𝑟 −𝑆𝑐,𝑔 where 𝑆𝑐,𝑟 is the red channel value for the stain vector of the class 𝑐, 𝑆𝑐,𝑔 the green channel value and 𝑆𝑐,𝑏 the blue channel counterpart. The stain vector is generated for each class 𝑐 by 𝑆𝑐=[𝑆𝑐,𝑟,𝑆𝑐,𝑔,𝑆𝑐,𝑏]𝑇, and, in turn, the stain matrix by 𝑆=⁡[𝑆1,𝑆2,𝑆3], being 𝑐=1,2,3 the H, E Methodology 50 and B stains, respectively. The stain density map 𝐷 is obtained through the inverse of the stain matrix as demonstrated in equation (15): 𝐷=𝑆−1𝐼𝑂𝐷 (15) The created stain matrix is then applied to the OD image to each pixel to obtain a final three channel image, composed by images from the three different stains, Hematoxylin, Eosin and Background. The stain images are separated and the first two, the relevant stains, suffer a contrast enhancement process. The contrast enhancement method embedded in the stain deconvolution framework performs a mapping of intensity values in grayscale images, so that the output data (image) presents a saturation of 1 % in high and low intensities, thus an increased contrast in intensities (MathWorks). The registration method implemented was the automatic intensity-based registration algorithm (section 3.2.2.) with transform type rigid and performed taking as constant fixed image the reference slice of the image dataset (the central slice), and as moving images all the slices successively, from the first to the last in the dataset, aligning them to the reference and stacking them to form volumetric data. To execute the 3D reconstruction and display were implemented using two different frameworks. The algorithms tested were the Marching Cubes algorithm (Lorensen and Cline 1987) implemented by (Hammer 2011) and the Isosurface framework from (MathWorks), adapted for the image datasets tested in this study. The Marching Cubes algorithm creates polygonal surface representations of isosurfaces (Lorensen and Cline 1987) of 3D scalar fields, such as the slice stack formed after the registration complete process over all the dataset images. The implemented method consists in a vectorized version of the algorithm, which computes a triangulated mesh of the isosurface within a given 3D matrix of scalar values, according to the isosurface value (the constant value or level set from the 3D image represented in the isosurface), and specified in terms of a face and vertex list. The orientation of the triangles generated is selected according to the normal point from higher to lower values (Hammer 2011). This method enables the visualization of the generated 3D surface, with customizable colors and view. The Isosurface framework from (MathWorks) also performs the extraction of the isosurface data from volumetric according to an isosurface value, although with a different process. This method performs first a volumetric interpolation using linear interpolation with (MathWorks), and then connects the points with equal isosurface value, forming a three dimensional surface. This method is combined with a 3D data smoothing technique (MathWorks) to smooth the edges of the reconstructed surface and a 2D polygon filling method (MathWorks) to enable its coloring and visualization. The final workflow schematic representation is shown in Figure 25. 3.2. Workflow implementation 51 3D Reconstruction Marching Cubes Isosurface Pre-Processing CLAHE in red channel Registration Intensity-based rigid registration with reference slice Segmentation Stain deconvolution Reference Slice New image Registered Image is saved Figure 25. Schematic model representing the third and final workflow implemented in the study. The approach starts by computing a loop for each slice in the dataset (minus the reference middle slice) in which the image suffers pre-processing, segmentation through stain deconvolution and registration with the reference slice as fixed image. All the registered images are stored and stacked to form the 3D data dataset, which will, posteriorly, serve as input of both 3D reconstruction frameworks. 52 53 Chapter 4 Results and Discussion In this chapter, all the experimental results obtained with the three approaches described in the methodology (Chapter 3) and their respective discussion are presented. The preprocessing results will be the first to be analyzed, based on color contrast enhancement and segmentation results, followed by registration frameworks comparison and evaluation (second approach). Finally, the final workflow will be discussed, finishing with the 3D reconstruction results. All the selected parameters and decisions performed are properly indicated and justified, throughout the entire chapter. The workflows in this study were tested on a Windows 64bit Intel Core i7-4700HQ CPU at 2.40 GHz, with 8 GB of RAM computer system. 4.1. Pre-processing stage In this section, results from the application of pre-processing methods, addressed in the first approach (section 3.2.1.) will be presented and discussed further in this section. In the end, the best methods or workflows will be determined based on visual interpretation and a segmentation test, performed with kmeans algorithm. Observing all the lesions and their description, in section 3.1.2., found in images from the three datasets, it can be concluded that they present very distinct natures and occur stained by both hematoxylin and eosin stains, in a wide range of color intensities. In light of the previous, the pre-processing techniques were developed focusing the color contrast enhancement and ultimately the distinction between tissues stained with each one of these stains. The resultant images after implementation of the pre-processing techniques, presented throughout this section, applied to a single slice from the second dataset, serving as an illustrative example. The original image from the selected slice (nº 35) is presented in Figure 26. Results and Discussion 60 Figure 34. Resultant images from YCbCr color space transformation applied to the normalized RGB image (Figure 30), with constant Cb factor (1/16) and decreasing luminance multiplying factors - 1/20 (a), 1/40 (b) and 1/80 (c). Analyzing the Figure 34, showing the result obtained when the YCbCr color space transformation is applied to the normalized image, outcome of the normalization method (Figure 30), it is evident the improvement in terms of image quality and sharpness. A general comparison between the three images shows that they are similar in terms of tissue colors among relevant structures, and thereby the diminution factors smaller than 1/20 do not increase the color contrast. The only difference between these images is the green color noise noticed among pixels from the background, which is higher for lower values. Considering this, the best suited diminution value for these images is 1/20x, since it provides a clearer discrimination between different stained tissues. This value, as well as the blue chrominance factor (1/16) were determined through experimentation, i.e. color contrast enhancement optimization. As noticeable in Figure 34 the diminution of both factors results in an overall image darkening (luminance diminution), particularly in tissues stained with hematoxylin, due to higher blue color influence (blue chrominance diminution). The YCbCr color space transformation preceded by an image normalization was the most accomplished method, of the presented in this section, in increasing the color contrast between Hematoxylin stained, in deep blue shades, and Eosin stained tissues, in magenta. The masking method developed for this study (section 3.2.1.) involved a combination of two color space transformation techniques, HSV and YCbCr, with the predetermined parameters (enhancement factors, 2x for HSV and 1/20x for YCbCr). An automatic threshold operation is followed, in order to create a binary image with the lymph node detached from the background. The final steps for this masking method includes two morphological operations, starting with an image closing, with a 4 pixel radius disk shaped structural element, a b c 4.1. Pre-processing stage 61 followed by an image opening performed with the same element, but with 6 pixel radius. These steps were conducted to achieve a unique final mask structure, as the one presented in Figure 35, for the original image (Figure 26). The masking process is important to increase the influence of each different tissue to be segmented in the segmentation procedure. Figure 35. Mask structure, created from the original image to remove the background (a). To corroborate the visual interpretation and selection of the most adequate pre-processing method to enhance color contrast based on the resultant images, in the datasets prepared for this study, a simple segmentation procedure was performed, using kmeans (section 2.2.2.) for each method and considering three different classes (hematoxylin stained tissues, eosin stained tissues and background), with 15 repetitions and using square Euclidean distance to assign the clusters. The segmentation results are shown in Figure 36 as well as the processed images used to acquire them, respectively, masked with the structure on Figure 35. Comparing the masked images originated from the different pre-processing methods and the respective three class segmentation result, represented in Figure 36, it is clearly noticeable that CLAHE and YCbCr transformation were the best performing methods in distinguishing the eosin and hematoxylin stained tissues. Although CLAHE was the top method to distinguish both stains in the segmentation (Figure 36.b), the YCbCr technique proved to be the best method to enhance the contrast between the lesion (neoplasic tissue) and healthy tissues, visible in the respective segmentation result (Figure 36.g). The three class segmentation result correspondent to the HSV enhancement presents mostly over segmented tissues, not being able to accurately find neither eosin nor hematoxylin stained tissues. On the contrary, the L*a*b transformation result presents under segmentation of both stains and an inconsistent pixel label throughout the entire lymph node. Facing the previous results, and in an attempt to improve the segmentation accuracy another test was conducted but considering four classes, dealing this way with tissues presenting lower concentration on eosin stain color, and possibly segment additional tissues, such as the blood vessels. The segmentation results are presented on Figure 37. Through comparison between the results for the four class segmentation and the 3 class counterpart in each image, it can be concluded that was an overall improvement in the tissue discrimination, more evident in the two methods with the worst outcome in the previous test (HSV and L*a*b color space transformation). Results and Discussion 62 The CLAHE-red segmented image presents a higher differentiation in the tissue surrounding blood vessels, maintaining almost unaltered the other pixel’s assigned labels, in relation to Figure 36.b. A similar comparison can be made for both YCbCr segmentation results, were the fourth label was assigned to pixels surrounding the blood vessels (Figure 37.d), although less evident than the CLAHE pixel labelling. Despite a considerable improvement in HSV transformed image segmentation, it does not present a satisfactory outcome, being a large amount of eosin stained tissue still with miss assigned labels. The CIE L*a*b results were clearly Figure 36. Masked images, resultant from the contrast enhancement techniques applied to the original image (left column) – CLAHE-red (a), HSV transformation (c), L*a*b transformation (e) and YCbCr transformation (g). Kmeans results of the respective left image (right column) considering 3 classes – CLAHE-red (b), HSV transformation (d), L*a*b transformation (f) and YCbCr transformation (h). a b c d e f g h 4.1. Pre-processing stage 63 better than the previous segmentation, but also inconsistent, with a considerable amount of pixels belonging to hematoxylin and eosin stains being assigned to a third label. All things considered, and to assess the efficiency of top-performing pre-processing workflows, for stained tissues contrast enhancement, when applied to the other two datasets, the implementation was made and the resultant images are presented in Figure 40, regarding case 1, and Figure 39 for case 3. The selected slices, to serve as example (slice nº35), from case 1 and case 3 are shown in Figure 38. Figure 37. Segmentation results considering four classes, for CLAHE image (a), HSV enhanced image (b), L*a*b transformed image (c) and YCbCr image (d). a b c d Figure 38. Original image (slice nº 35) from case 1 (on the left). Original image (slice nº35) from the third dataset (on the right). Results and Discussion 64 Analyzing the outcomes from the segmentation process for both slices, from case 1 (Figure 39Figure 40.b and Figure 39Figure 40.d) and case 3 (Figure 40.b and Figure 40.d), it can be concluded that both pre-processing workflows produced a satisfactory segmentation, thus achieving an accurate discrimination of both different stains, hematoxylin and eosin, in the tissues. The fact that both automatic pre-processing techniques demonstrated consistency in Figure 39. Masked images, resultant from the contrast enhancement techniques applied to the original image from case 1 – CLAHE-red (a), YCbCr (c). Results from kmeans segmentation performed in the respective left image (b, d). Figure 40. Masked images, resultant from the contrast enhancement techniques applied to the original image from case 3 – CLAHE-red (a) and YCbCr (c). Results from kmeans segmentation performed in the respective preprocessed image, for CLAHE-red is image (b) and for YCbCr is (d). a b c a d c d b 4.1. Pre-processing stage 65 datasets with large morphological and color differences (case 1, case 2 and 3), reinforces their quality, and thus their automatic color contrast enhancing capability. Comparing both methods, it is noticeable a higher differentiation between distinct stains (Hematoxylin and Eosin) with the CLAHE-red approach, being the correspondent tissues accurately separated, even in the collision tumor where both cell lines were clearly identified. On the other hand, the YCbCr workflow produced a tissue differentiation in the collision tumor in case 1 and divided different tissues near the inner neoplasic tissue in the case 3 image. During the pre-processing technique’s development and testing, they presented several intra an inter-case inconsistencies in terms of contrast enhancement efficiency (visible when comparing Figure 36, Figure 40 and Figure 39) and, therefore, in the final segmentation accuracy. These errors are mainly due to staining incongruities in slices from different datasets and even from slices of the same case study, a fact already mentioned in section 3.1.2., presenting, in the last scenario, identical structures with slightly different colors. Since these pre-processing techniques rely on pixel intensity and color properties and have to be automatic, the aforementioned fact is the most probable cause for some errors in the segmentation. Despite some techniques, and respective results, being explored and described to mask blood vessels from the images, namely, CIE L*a*b and normalization techniques, after several trials, intense testing and combination of techniques, the development of an automatic method capable of performing accurate distinction of blood vessels, in all images contained in the three datasets was proven impossible. In some cases, only the exterior blood vessels were removed, in other cases occurred a partial removal on the same. Facing this fact, the regions containing blood vessels in the images were not removed with the developed mask technique, implemented in the final pre-processing workflow. CLAHE-red histogram equalization with background removal (masking process) and the YCbCr color space, also with background removal, preceded by a normalization step, both top performing methods in the segmentation process, were implemented to all slices from the three image datasets, in order to improve image registration, a process exhaustively explored in the next section. Results and Discussion 66 4.2. Registration stage In the present section, the registration methods, or combination of methods, are going to be explored and their efficiency in slice alignment will be compared, based on the DICE score and the computational cost required to compute the entire dataset image registration. The computational cost was measured as the time required to complete the registration of all the slices. The dataset selected to test the registration methods was case 3, containing 100 slices, since it is the complete dataset with the lowest number of images. All the experiments, using different registration algorithms were conducted for reference slice and pairwise registration models (section 3.2.2.). All the images submitted to a registration method suffered a previous grayscale conversion followed by a contrast adjustment with (MathWorks). The first method explored was the automatic - intensity based technique, due to its simplicity and its three types of transformation (rigid, similarity and affine). The best suited pre-processed image to perform registration, among the resultant from the pre-processing methods, described in the previous section, was determined for this first test and used for the other intensity-based methods (Demon and B-spline algorithms). The mean DICE score, average DICE score from the 100 registered slices, measured to assess the registration quality, and the time elapsed during the registration of the entire dataset are presented in Table II, for all the experiments conducted. Table II. Table containing the mean DICE scores and the elapsed times for the intensity-based registration implementation, with different types of transformation (rigid, similarity and affine) and models (reference slice, pairwise), on images from both top-performing pre-processing methods (YCbCr color transformation and CLAHE-red). Firstly, the images resulting from the YCbCr color space transformation (the method was applied to all slices in the dataset) masked and unmasked were tested with rigid Pre-processing method Transformation Intensitybased registration Reference Slice Pairwise YCbCr color space transformation With Mask Rigid DICE: 0.6632±0.1739 Time: 54.808968 s DICE: 0.6425±0.2612 Time: 53.710495 s Without Mask Rigid DICE: 0.9114±0.0288 Time: 66.556130 s DICE: 0.8897±0.0411 Time: 66.166884 s Similarity DICE: 0.9033±0.0642 Time: 81.728404 s DICE: 0.8906±0.0422 Time: 72.233815 s Affine DICE: 0.9054±0.0626 Time: 74.174649 s DICE: 0.8897±0.0432 Time: 67.782040 s CLAHE-red Without Mask Rigid DICE: 0.9683±0.0283 Time: 55.990489 s DICE: 0.9419±0.0385 Time: 47.265646 s Registration stage 67 transformation, to assess if the background removal would improve the registration process, as it did with the segmentation (section 4.1.). Comparing the DICE scores it is clear that without masking process the registration is considerably more accurate (Masked DICE = 0.6632±0.1739; Without mask DICE = 0.9114±0.0288; for reference slice model), although with a higher computational cost. These results prove not only that the masking method implemented is irregular in terms of removed structures, but also that a masked image possesses far lesser information, for an intensity-based registration method, than the same image with background. Considering the discrepancy between the obtained DICE scores, the following experiments were all conducted without masked images, to comprehend which one of the registration types performs the best registration considering the available datasets. Comparing the higher rigid score with both similarity and affine DICE scores it is possible to deduce that when the reference model is performed the rigid transformation slightly outperforms the other two types when the model is the reference slice (rigid DICE = 0.9114±0.0288; similarity DICE = 0.9033±0.0642; affine DICE = 0.9054±0.0626), and, when implementing the pairwise framework the three scores are similar (rigid DICE = 0.8897±0.0411; similarity DICE = 0.8906±0.0422; affine DICE = 0.8897±0.0432), with slightly higher DICE result for the similarity type. After this test, the set of YCbCr images was substituted by the resulting images from the CLAHE-red method to verify which pre-processing technique produces the best images for intensity-based registration. This test was conducted without background removal and to the top scoring registration type – rigid registration. The comparison between the DICE scores indicates that the best pre-processing method, for the considered dataset intensity-based image registration is the CLAHE-red (for reference slice model CLAHE rigid DICE = 0.9683±0.0283, YCbCr rigid DICE = 0.9114±0.0288; for pairwise model CLAHE rigid DICE = 0.9419±0.0385, YCbCr rigid DICE = 0.8897±0.0411). All the tests performed with reference slice model consistently present a higher DICE score when confronted with their pairwise model counterparts, suggesting that, for this case 3 the best model to perform intensity-based registration, and therefore, the entire image dataset alignment is the registration of each image to the reference slice. In the present work the defined reference slice is the middle one (slice nº50 in the considered dataset – case 3), because it is the most representative in terms of tissues, according to (Chakravarty, Bedell et al. 2008)). Analyzing the computing time of all the different registration types the results (Table II), it can be observed that the rigid transformation is the less computationally expensive, followed by the affine and then the similarity transform, with the highest time. These results go according to expected, because the rigid transform consists only in translation and rotation, being simpler than the other two (similarity – translation, rotation and scale; affine – translation, rotation, scale and sheer). Comparing, on the other hand, the times obtained by both registration models, it is visible a consistent increased computational cost in the registration methods performed using reference slice. Results and Discussion 68 Alignment images and pair of fixed and moving image for the most relevant registration experiments are shown in Figure 42, Figure 42 and Figure 43. The chosen moving image is slice nº35, the same as in the pre-processing section (4.1.). Figure 42. Images representing intensity-based image registration pairwise model and rigid transformation type, from YCbCr preprocessing. Image (a) represents the fixed image (slice nº36), (c) the moving image (slice nº35) and (b) the rigid registration overlapping image (Fixed and registered images). b a d c a b c Figure 41. Images representing intensity-based image registration performed with reference slice model and two types of transformation, from YCbCr pre-processing. Image (a) represents the fixed image (slice nº50) and (c) the moving image (slice nº35) for both rigid (b) and affine (d) registration type overlapping image (Fixed and registered images). Registration stage 69 After experimenting the automatic intensity-based registration algorithm, explained in section 3.2.2., an automatic feature-based similarity registration approach was implemented through an algorithm also referred in this section. Based on the top-scoring results obtained with the first method, the feature-based technique was applied to the CLAHE-red processed images, but without any viable results. This was due to the lack of matching points verified in the registration process of most of the slices contained in the dataset, either recurring to reference slice or pairwise model. Taking this into account and that feature-based methods perform image alignment based on features correspondence (between moving and fixed images), such as points, lines and contours, another pre-processing technique was applied to the dataset images, the image normalization (section3.2.1.). This method was implemented due to its capability to detach the lymph node (in this case) from the background, promoting its contour extraction. Aiming also to improve the contour extraction for the feature registration method a gaussian smoothing filter, with kernel dimensions 4x4, was applied to the normalized image. The registration results, DICE mean score and computational time, for smoothed and non-smoothed images are presented in Table III. The reference slice registration model was tested for this set of images, but it was incapable of providing a full dataset registration, and thus, the DICE score was not measured. Table III. Table containing the mean DICE scores and the elapsed times for the featurebased registration implementation with pairwise model on images with both smoothed and non-smoothed image normalization pre-processing. Pre-processing method Feature-based similarity registration (Pairwise model) Image Normalization With smoothing DICE score: 0.9218±0.0481 Time: 30.537524 s Without smoothing DICE score: 0.8412± 0.0999 Time: 27.915660 s a b c Figure 43. Images representing intensity-based image registration performed with reference slice model and rigid transformation type, from CLAHE-red pre-processing. Image (a) represents the fixed image (slice nº50), (c) the moving image (slice nº35) and (b) the rigid registration overlapping image (Fixed and registered images). Results and Discussion 76 Analyzing Table VI it was possible to conclude which of the stain deconvolution outcome images is the most suited to perform image registration. The highest mean DICE scores verified belong to registrations performed with hematoxylin images, in all three datasets, justified by its greater representativeness in terms of high intensity structures, as it can be visualized in Figure 49. In case 2 both images, hematoxylin and eosin, share similar structures, with high pixel intensity, a fact, corroborated by the closest proximity in the DICE scores measured. In terms of computational time, the less demanding are the dataset registrations performed with the eosin images, because of the same reason explained for the DICE score, less amount of pixels to be correlated in the procedure. It is also noticeable that, for case 1 the registration was not only computationally expensive but also inaccurate and inconsistent (high standard deviation values), with mean DICE scores below those verified in for CLAHE-red (Table V), suggesting that the stain deconvolution proved itself incapable to distinguish both tissue stains in a substantial portion of slices, from this dataset. Figure 49. Segmentation results obtained with the stain deconvolution algorithm (first column): case 1-(a), 2-(d) and 3 (g). Hematoxylin image extracted with this algorithm (second column): case 1-(b), 2-(e) and 3-(h). Eosin image extracted (third column): case 1-(c), 2-(f) and 3-(i). h g f e d c b a i Final Results 77 Table VI. Table containing the mean DICE scores and the elapsed times for the intensity-based rigid registration method implementation, on all the image datasets (Case 1, 2 and 3), considering two different grayscale images (hematoxylin and eosin images), obtained through pre-processing (CLAHE-red and YCbCr transformation) and segmentation (stain deconvolution algorithm). Case study Registration Hematoxylin image Eosin image 1 DICE: 0.8387±0.2785 Time: 71 min (approx.) DICE: 9.4781e-04±0.0102 Time: 70 min (approx.) 2 DICE: 0.9345±0.0368 Time: 204.695529 s DICE: 0.9179±0.0610 Time: 155.079913 s 3 DICE: 0.9676±0.0296 Time: 183.110440 s DICE: 0.6152±0.0894 Time: 158.491751 s Considering the previous results the selected set of grayscale images chosen to find the registration transformation were the hematoxylin stain images. The transformation found on each slice registration with the reference was applied to the correspondent eosin image, being aligned with the same translation and rotation than the hematoxylin image. Hereupon, both grayscale image datasets (hematoxylin and eosin) were equally aligned, and this framework was performed for each case study. After registering the complete image dataset to the reference slice, all images were stacked creating volume data, in which the z coordinate represents the slice number. All the image datasets (case 1, 2 and 3) were stacked in two rectangular prisms each, one for the hematoxylin and the other for eosin registered images, with length and width defined by the length and width of the reference slice (because the registered image acquires the size of the fixed image), and with height equal to the number of slices contained in the dataset. These final volumes were smoothed with a box filter (a spatial domain low-pass filter in which each point in the output is the average of the surrounding points in the input structure), with 5x5 Figure 50. Images representing intensity-based image registration performed with reference slice model and rigid transformation type, from hematoxylin images obtained through stain deconvolution (Case 3). Image (a) represents the fixed image (slice nº50), (c) the moving image (slice nº35), (b) the registration overlapping image and (d) the transformed eosin stained slice (nº35), with the registration transformation determined for the hematoxylin images. d c b a Results and Discussion 78 kernel, and served as input for both 3D surface reconstruction methods, the Marching Cubes algorithm and the isosurface framework from (MathWorks) described in section 3.2.3. Using the marching cubes algorithm it was only possible the reconstruction of one stained surface at the time, being represented in Figure 51 the hematoxylin and eosin reconstructions from case 3 image dataset, with isovalue equal to 0.5 for hematoxylin stack and 0.65 for eosin. Both isovalues were determined through experimentation, aiming to reconstruct a complete 3D surface for the considered tissue (stained with hematoxylin or eosin) and at the same time reduce the noise produced by background artifacts with intensities similar to the reconstructed structure, thus, reducing its intelligibility to the observer. Figure 51. 3D surface reconstruction from case 3 registered image dataset, using Marching cubes algorithm. Images (a) and (b) are from hematoxylin image dataset surface reconstruction with different views. Images (c) and (d) are eosin stacked images’ surface reconstruction. Structure (1) represents the neoplasic tissue found in Figure 20 – A; structure (2) represents the neoplasic tissue marked in Figure 20 – B and structure (3) is the blood vessel in Figure 20 – C. a b c d 1 1 2 2 3 3 Final Results 79 Implementing the isosurface framework it was possible to combine both tissue stains stacks’ surface reconstruction in one image, using the patch function (MathWorks), and the result obtained for case 2 and 3 can be consulted in Figure 52. To each tissue stain reconstruction was assigned a color, based on the real color presented in the original image dataset. This was accomplished using an auxiliary algorithm, which, first, separates both stains from the CLAHE-red image based on stain deconvolution segmentation, performs an average pixel intensity measurement in the three color channels of both stain images from the slice and then, after the dataset complete analysis, performs another average calculation for all the slices, obtaining two final RGB triplets to be applied to the respective stained tissue 3D reconstruction. For this method were applied two different isovalues, 0.5 and 0.65 (same as in the previous method) for hematoxylin and eosin stacks, respectively. Observing Figure 51 and Figure 52, it is possible to conclude that both cases 2 and 3 present a consistent slice alignment, visible in the surface uniformity along the Z-axis. This is more evident in the marching cubes algorithm results, suggesting that the 3D interpolation performed by this method provides more accurate results. The success of the entire framework is proven not only by this feature, but also by the similarity, in terms of morphology between what is seen in the original 2D image slices, from both datasets, and the result of the 3D surface reconstruction. Focusing on Figure 51 and in case 3 tissue analysis (section 3.1.2.) the Figure 52. 3D surface reconstruction from case 2 (second line) and 3 (first line) registered image dataset, using isosurfaces framework. Image (a)-top view, (b)-bottom view and (c)-side view, represent the case 3 volumetric data combining both hematoxylin and eosin isosurfaces in three different views. Image (d)-top view, (e)-bottom view and (f)-side view, represent the case 2 volumetric data combining both hematoxylin and eosin isosurfaces in three different views. a f e d c b Results and Discussion 80 resemblances are obvious, therefore, enabling the correlation between some relevant structures identified in Figure 20. The most evident structures reconstructed are the neoplasic tissue stained with hematoxylin (Figure 51.a – 1) and the other two stained with eosin (Figure 51.d – 2), and also the blood vessel (Figure 51.c – 3 and Figure 20 – c). These 3D surface reconstructions are consistent and are able to provide important volumetric information about histological slices contained in the datasets, converting a planar observation in a 3D visualization of the tissues stained with hematoxylin and eosin in an either separate or merged version. Analyzing the surface reconstruction performed for case 2 (Figure 52– d,e,f), it is verified a reduced eosin tissue definition, maybe due to the sparse nature of the tissue on the original dataset that decreases the intelligibility of the volumetric data. The colors expressed in the final surface reconstructions present high affinity with the stain colors verified in the CLAHE-red images in Figure 29, pointing the efficiency of the previously described color algorithm in combination with the isosurface framework. The 3D surface reconstruction framework applied to case 1 presented some computational issues, which compromised the result demonstration. The processing power used to test all the frameworks in this study (Intel Core i7-4700HQ CPU at 2.40 GHz, with 8 GB of RAM on a Windows 64bit system) was not sufficient to load and visualize the volumetric data, mainly due to the bad performance of the stain deconvolution method in discriminating tissue stains in most slices contained in the dataset and the large image dimensions. Taking into account all the previous 3D surface reconstructions, resulting from the final framework developed for this study, it can be concluded that the three-dimensional representation of the datasets was accomplished for the majority of the available datasets, and with interesting and reliable results. Especially for case 3, the reconstruction was able to provide accurate 3D volumetric features from neoplasic tissues in addition to overall morphological characteristics of the lymph node. The most relevant are the three dimensional progression on the Z-axis of the lesion over other tissues (Figure 51), difficult to visually estimate from a set of 2D slices. 4.4. Summary In summary, the top performing pre-processing workflows for hematoxylin and eosin stain contrast enhancement were the CLAHE-red histogram equalization followed by a background removal (masking process) and the YCbCr color space, also with background removal, preceded by a normalization step. After testing several methods, it was concluded that the most efficient framework to perform registration in the considered histological image datasets is the automatic intensitybased registration with rigid transformation type and reference slice model on CLAHE-red preprocessed slices. This was deduced based on both DICE score results and computational cost, being the method with lower times of execution and the second highest DICE scores. The Summary 81 highest score was obtained with B-spline non-rigid algorithm, but, since this method originated distortion artifacts to provide accurate registrations, corrupting valuable tissue information, and the times required to compute entire datasets could reach approximately one hour, the rigid method was considered the most suited for the task. Combining the best methods from both pre-processing and registration sections, the final workflow was completed with the 3D surface reconstruction of both hematoxylin and eosin stained tissues registered stack. The eosin and hematoxylin separate images were acquired from an intermediate segmentation step, the stain deconvolution algorithm, a specially developed method for stained tissues differentiation from histological images. Through the combined implementation of these four main steps an accurate and interesting volumetric representation of both hematoxylin and eosin stained tissues was determined. Both 3D surfaces can be shown together or in separate images. 82 83 Chapter 5 Conclusion Automated techniques for image processing and analysis applied to histological studies constitute an important asset in the comprehension and assessment of relevant aspects, volumetric data and models from numerous diseases, unobtainable through 2D tissue slice analysis. Impelled by the relevance of this studies and the urge to assist the histopathologists in performing medical diagnostic, more accurate and reliable algorithms have been developed, providing a fast and trustworthy second opinion. The present study explores histological concepts, including tissue types and sample preparation as well as successful techniques of image processing and analysis, highlighting current reliable approaches in the field of 3D histological tissue reconstruction. With the previous research in mind and considering the need to produce a 3D reconstruction approach for cellular tissue analysis, in this work it was developed an automated computational framework combining the best algorithms and techniques implemented by the reviewed approaches, focused not only on reconstructing the relevant tissues, but also on their spatial relation with the surrounding structures. After conducting the present study, it is possible to conclude that several artifacts present in histological tissue sections, such as differences in stain colors between slices and/or image datasets, can reduce the consistency of both image processing and registration methods when applied to them. On the other hand, through the implementation of accurate registration techniques on the image datasets, the presence of tissue orientation dissimilarities, due to the manual nature inherent to the sample preparation process, was almost mitigated. The color contrast enhancing pre-processing methods developed and applied to the original digital images effectively provided an improvement in discriminating different stained tissues. This can be inferred either from visual analysis of resulting images or from the stain deconvolution segmentation result’s accuracy, and, consequently high performance in the registration process, with DICE scores of 0.9345±0.0368 for case 2 and 0.9676±0.0296 for case Conclusion 84 3. Case 1 achieved poorer results in general, with a final DICE score of 0.8387±0.2785 in the registration process, justified by an inconsistent stain deconvolution procedure, revealing itself insufficient to accurately separate both stains, in a great portion of images from this dataset. A consistent neoplasic tissue discrimination was not obtained through the implemented algorithm, but, by observing the 3D surface reconstruction of both eosin and hematoxylin (section 4.3.) the presence of this tissue is evident. Despite being specified the lesions in section 3.1.2. as tissues to be targeted by the algorithm, aiming to promote an isolated 3D reconstruction of these neoplasic tissues, the stained tissue separation and posterior 3D reconstruction reveals the interaction between the neoplasia and the surrounding healthy tissues. Therefore, the consistent reconstruction of both stained tissues including the lesion, with accurate volumetric shape, constitutes relevant information for diagnostic purpose. Considering all this, it is concluded that the algorithm developed in this study provides a fast and accurate 3D surface volume, comprising both distinct stains represented independently and showing the three dimensional relation between lesions and the surrounding tissue. The proposed framework obtained a satisfactory performance for two of the three available datasets. 5.1. Future Work Perspectives One future work improvement is the preparation and testing of more histological image datasets to further validate the developed algorithm, and also to assess the algorithm’s accuracy when analyzing tissues without H & E stain. This would also instigate a refinement in the computational techniques, adapting the automatic framework to more sets of images and develop a more universal algorithm, capable of performing detailed tissue 3D reconstruction of any histological section. Another future improvement would be the incorporation of a functional and general lesion detector, able to distinguish special tissues from the other stained healthy ones. The improvement of the stain deconvolution method, implementing supervised classification of the stained tissues, using as training set a portion of the entire image dataset, could lead to more accurate staining discrimination and possibly to define other interesting tissues, such as neoplasic tissues. An interesting future development would be the implementation of the produced computational framework in opencv (Open Source Computer Vision), enabling a faster image processing, for larger image datasets, such as case 1, and the creation of an intuitive user interface for input image dataset manipulation and 3D reconstructed volume observation. 85 References Arganda-Carreras, I., R. Fernandez-Gonzalez and C. Ortiz-de-Solorzano (2004). Automatic registration of serial mammary gland sections. Engineering in Medicine and Biology Society, 2004. IEMBS'04. 26th Annual International Conference of the IEEE, IEEE. Aubert, G. and P. Kornprobst (2006). Mathematical problems in image processing: partial differential equations and the calculus of variations, Springer Science & Business Media. Bautista, P. A. and Y. Yagi (2010). "Improving the visualization and detection of tissue folds in whole slide images through color enhancement." Journal of pathology informatics 1: 1-25. Bay, H., A. Ess, T. Tuytelaars and L. Van Gool (2008). "Speeded-up robust features (SURF)." Computer vision and image understanding 110(3): 346-359. Beare, R. and G. Lehmann (2006). "The watershed transform in ITK-discussion and new developments." Insight J: 1-24. Bengtsson, E. (2003). Computerized cell image analysis: Past, present, and future. Image Analysis, Springer: 395-407. Blake, C. A., H. A. Lavoie and C. F. Millette (2003). "Teaching medical histology at the University of South Carolina School of Medicine: transition to virtual slides and virtual microscopes." The Anatomical Record Part B: The New Anatomist 275(1): 196-206. Bonnet, N. (2004). "Some trends in microscope image processing." Micron 35(8): 635-653. Bookstein, F. L. (1980). "The Measurement of Biological Shape and Shape Change." Syst. Zool 29(1): 102-104. Braumann, U.-D., J.-P. Kuska, J. Einenkel, L.-C. Horn, M. Loffler and M. Hockel (2005). "Threedimensional reconstruction and quantification of cervical carcinoma invasion fronts from histological serial sections." Medical Imaging, IEEE Transactions on 24(10): 1286-1307. Caicedo, J. (2009, 06 November). "Features for Histology Images." Retrieved 03 February, 2015, from http://www.informed.unal.edu.co/jccaicedo/docs / review.pdf. Chakravarty, M. M., B. J. Bedell, S. P. Zehntner, A. C. Evans and D. L. Collins (2008). Threedimensional reconstruction of serial histological mouse brain sections. Biomedical Imaging: From Nano to Macro, 2008. ISBI 2008. 5th IEEE International Symposium on, IEEE. Chakravarty, M. M., G. Bertrand, C. P. Hodge, A. F. Sadikot and D. L. Collins (2006). "The creation of a brain atlas for image guided neurosurgery using serial histological data." Neuroimage 30(2): 359-376. Chan, A. C. S. (2013, 21 July 2013). "Tumor cell extraction by color deconvolution." Retrieved 15 June, 2015, from http://web.hku.hk/~ccsigma/color-deconv/color-deconv.html. References 92 Teot, L. A., R. Sposto, A. Khayat, S. Qualman, G. Reaman and D. Parham (2007). "The problems and promise of central pathology review: development of a standardized procedure for the children's oncology group." Pediatric and Developmental Pathology 10(3): 199-207. Thirion, J.-P. (1998). "Image matching as a diffusion process: an analogy with Maxwell's demons." Medical image analysis 2(3): 243-260. Torr, P. H. and A. Zisserman (2000). "MLESAC: A new robust estimator with application to estimating image geometry." Computer Vision and Image Understanding 78(1): 138-156. Treanor, D., N. Jordan‐Owers, J. Hodrien, J. Wood, P. Quirke and R. A. Ruddle (2009). "Virtual reality Powerwall versus conventional microscope for viewing pathology slides: an experimental comparison." Histopathology 55(3): 294-300. Unpublished, W. (2015). A Discriminative Framework for Stain Deconvolution of Histopathology Images in the Maxwellian Space. Vercauteren, T., X. Pennec, A. Perchant and N. Ayache (2009). "Diffeomorphic demons: Efficient non-parametric image registration." NeuroImage 45(1): S61-S72. Vernon, D. (1991). Machine vision-Automated visual inspection and robot vision, Prentice Hall. Wang, H., L. Dong, J. O'Daniel, R. Mohan, A. S. Garden, K. K. Ang, D. A. Kuban, M. Bonnen, J. Y. Chang and R. Cheung (2005). "Validation of an accelerated'demons' algorithm for deformable image registration in radiation therapy." Physics in medicine and biology 50(12): 2887. Watson, M. L. (1958). "Staining of tissue sections for electron microscopy with heavy metals." The Journal of Biophysical and Biochemical Cytology 4(4): 475-478. Weszka, J. S. (1978). "A survey of threshold selection techniques." Computer Graphics and Image Processing 7(2): 259-265. Xu, R. and D. Wunsch (2005). "Survey of clustering algorithms." Neural Networks, IEEE Transactions on 16(3): 645-678. Zitova, B. and J. Flusser (2003). "Image registration methods: a survey." Image and vision computing 21(11): 977-1000.