Linköping University | Department of Biomedical Engineering Master’s thesis, 30 hp | Biomedical Engineering Spring 2016 | LiTH-IMT/ERASMUS-R-13/43-SE
Processing of Optical
Coherence Tomography Images
- Filtering and Segmentation of
Pathological Thyroid Tissue
Uppsat
sens titel på engelska kan placeras här
– inklusive undertitel i samma storlek
Daniela Koller
Författare 2:s Namn
Supervisor: Neda Haj-Hosseini Examiner: Göran Salerud
Linköping University SE-581 83 Linköping 013-28 10 00, www.liu.se
Abstract
In the human body, the main function of the healthy thyroid gland is the regulation of the metabolism and hormone production. Included in the thyroid are organized structured and uniformly shaped follicles ranging from 50 −500 µm in diameter. Pathologies lead to morphological changes of these follicles, affecting the density and size, but can also lead to an absence. In this study optical coherence tomography (OCT) was used to examine pathological thyroid tissue by extracting structural information of the follicles from image segmentation. However, OCT images usually include a high amount of speckle noise which affects the segmentation outcome. Due to that, the OCT images need to be improved. The aim of this thesis was to investigate the appropriate filtering methods to enhance the images and thus improve the segmentation outcome.
The images of pathological thyroid tissues with a size of 0.5 −1 cm were scanned by a spectral domain OCT system (Telesto II, Thorlabs GmbH, Germany) using a center wavelength of 1300 nm. The obtained 2D and 3D images were saved as .oct file as well as implemented and visualized in a MATLAB graphical user interface (GUI) for further processing. For image improvement, four filtering enhancement methods were applied to the 2D images such as the enhanced resolution imaging (ERI), adaptive Wiener filter, discrete wavelet transform (DWT) and multi-frame wavelet transform (WT). The processed images were further converted to grayscale and binary images for intensity-based segmentation. The output of all methods were compared and evaluated using signal-to-noise ratio (SNR), contrast-to-noise ratio (CNR), enhanced number of looks (ENL), edge profiles and outcome of the segmented images.
It was demonstrated that the complex DWT (cDWT) with a higher threshold and the multi-frame WT using the haar wavelet showed enhanced results over the other filtering methods. The computed SNR could be increased up to 52% and the ENL value up to 4802%, applying the multi-frame WT, while the CNR could be increased up to 106% for cDWT. The lowest obtained gradient was equal to an intensity decrease of −61% and −68% for multi-frame WT and cDWT, respectively. The filtering method could increase the smoothness of the image while the edge sharpness could be kept. The segmentation could detect both small and large follicles. ERI did not show any improvement in the segmentation but could enhance the structural detail of the image. Larger neighbourhoods of the adaptive Wiener filter showed a highly blurred image and led to merged follicles in the image segmentation.
The wavelet filters DWT and multi-frame WT gave most satisfying results since high and low frequencies were divided into subbands, where individual information on vertical, horizontal and diagonal edges was stored. Applied cDWT had an even higher amount of subbands, so that more information on signal and speckle noise could be specified. Due to this fact, it was possible to achieve a decreased noise level while edge sharpness were maintained. Using a multi-frame image an increased SNR was obtained, as the intensity information stayed constant over the individual frames while the noise information changed.
Wavelet based filtering showed higher improved results in comparison to the adaptive Wiener filter or the ERI in the 2D domain. By applying filtering methods in higher dimensions such as 3D or even 4D, better results in noise reduction are expected. Improved settings for the individual filtering methods as well as enhancement in segmentation are part of the future work.
Acknowledgements
First and foremost, I would like to thank my family and friends from Vienna, supporting me even over a far distance. I would also like to thank Malte Behr for his clarifying and time-saving support in some mathematical issues when my mind reached an impasse. Credit deserves also Olivier Cros and Anders Eklund for their support in medical image filtering.
Furthermore, I want to thank my supervisor Neda Haj-Hosseini and my examiner Göran Salerud. I would also like to thank the Department of Surgery at Linköping University Hospital for providing the tissue samples. I am also very thankful to both, Fiona Simpson and my college Martha Velasco Santascoy for proof-reading my thesis.
Especially the very last weeks were highly work-intensive and demanding. Ankica Babic kept up my cogitation even after long working hours with endless chocolate supply, encouraging conversations and her happy mood. Martha, too, proofed her exceptional ability to cheer me up in every situation with her sunny personality. Thanks for your support!
The time I spend in Linköping was unforgettable and I will always keep it in good memory. I am looking forward to coming back to Sweden but also visiting the home countries of my international friends I had the pleasure to meet here.
You’re braver than you believe, and stronger than you seem, and smarter than you think. – W.t.P.
Life is what happens while you are busy making other plans.
Table of Content
1 Introduction 1
1.1 Medical Imaging Technologies . . . 1
1.1.1 OCT versus Microscope Slides . . . 2
1.1.2 OCT versus Micro-CT System . . . 3
1.2 Optical Coherence Tomography (OCT) . . . 3
1.2.1 Imaging Principles . . . 3
1.2.2 Detection Principles . . . 5
1.2.3 Resolution and Field of View . . . 6
1.2.4 Artifacts and Noise . . . 7
1.3 Medical Application . . . 10
1.3.1 Medical Use of OCT . . . 10
1.3.2 Thyroid Gland . . . 10
1.4 Optical Properties . . . 12
1.4.1 Refractive Index . . . 12
1.4.2 Absorption Coefficient . . . 12
1.4.3 Scattering Coefficient and Anisotropy . . . 12
1.4.4 Reduced Scattering Coefficient . . . 13
1.4.5 Tissue Optical Properties for a 1300-1350nm OCT Laser Source . . . 14
1.5 Image Processing . . . 16
1.5.1 Speckle and Speckle Noise . . . 16
1.5.2 Filtering and Denoising Methods . . . 16
1.5.3 Enhanced Resolution Imaging . . . 18
1.5.4 Image Segmentation . . . 19
2 Materials and Methods 21 2.1 The OCT System . . . 21
2.1.1 Hardware . . . 21
2.1.2 Software and Settings . . . 21
2.2 Tissue Specimen . . . 22
2.3 Image Processing . . . 23
2.3.2 Image Filtering and Enhancement Methods . . . 25
2.3.3 Image Filtering and Enhancement Evaluation . . . 26
2.3.4 Image Segmentation . . . 29
3 Results 31 3.1 Data Implementation and Visualization . . . 31
3.2 Image Filtering and Enhancement Methods . . . 33
3.2.1 Enhanced Resolution Imaging . . . 33
3.2.2 Adaptive Wiener Filter . . . 35
3.2.3 Discrete Wavelet Transform . . . 37
3.2.4 Multi-Frame Wavelet Transform . . . 40
3.2.5 Summary of Evaluation Parameters . . . 43
3.3 Image Segmentation . . . 44
3.3.1 Original Image . . . 44
3.3.2 Enhanced Resolution Imaging . . . 44
3.3.3 Adaptive Wiener Filter . . . 45
3.3.4 Discrete Wavelet Transform . . . 45
3.3.5 Multi-Frame Wavelet Transform . . . 46
4 Discussion 49 4.1 Data Implementation and Visualization . . . 49
4.2 Image Filtering and Segmentation Methods . . . 49
4.2.1 Enhanced Resolution Imaging . . . 50
4.2.2 Adaptive Wiener Filter . . . 50
4.2.3 Discrete Wavelet Transform . . . 51
4.2.4 Multi-Frame Wavelet Transform . . . 52
4.2.5 Evaluation Parameters . . . 52 4.3 Future Work . . . 53 5 Conclusion 55 Bibliography 57 List of Figures 63 List of Tables 65 Acronyms 67
Chapter
1
Introduction
The thyroid gland is an essential organ in the human body as it regulates the metabolism and cell activity. Especially the follicles are responsible for regulating the hormone production, help to protect against diseases or to repair damaged tissue. In diseased tissue, the ordinary function of the thyroid gland is affected and more differences to healthy thyroid tissue can be obtained. Until now it is known that the follicles within the tissue are highly affected by a disease and show morphological changes in shape, size, number and density. Optical coherence tomography (OCT) can be used in addition to other examination methods to analyze pathologies. As OCT provides clear images of tissue structures, it can be a valuable tool for studying the follicular morphology of the thyroid and its respective pathology.
Image segmentation should provide more information on follicular properties. So, properties such as size, diameter, perimeter or density can be calculated and in a perfect case related to a specific disease. An optimum combination of image filtering and segmentation should give the opportunity to extract morphological parameters of follicles in the thyroid. Information on the follicles can be used for intraoperative decision making and for a quick postoperative histology analysis prior to the standard histology examination.
1.1
Medical Imaging Technologies
In the health care sector medical imaging techniques, such as X-ray computer tomography (CT), magnetic resonance imaging (MRI), ultrasound (US), and nuclear imaging modalities as positron emission tomog-raphy (PET) and single-photon emission CT, are essential to improve diagnosis and individual treatment of patients. Even though a high light penetration can be reached into tissue, the low spatial resolution can limit the applications for detailed microscopic information. To overcome these restrictions optical imaging methods such as fluorescence, conventional and confocal microscopy or multi-photon imaging can be used. They have a spatial resolution range of several micrometers [1, 2].
The resolution and penetration depth of OCT, illustrated in Figure 1.1, lies between optical microscopy methods and US imaging [1]. Microscopy generates images by transmission or optical reflection dif-ferences [2] having a resolution of 1 µm [1]. For US the penetration depth is greatly increased, so that the absolutely highest resolution of approximately 150 µm [1] can be achieved by measuring intensity differences of back-scattered sound waves within the tissue [2].
As shown in Table 1.1, different types of imaging systems and technologies can be used for medical applications to get more information on the tissues properties facilitating diagnosis of tissue samples [3].
Figure 1.1 – Comparison of penetration depth and resolution of commonly used medical imaging technologies referred to Popescu et al. [1].
Angiography and fluorescence microscopy can provide higher resolution and smaller fields than OCT, but it needs systemic labelling, which leads to limitations in longitudinal studies [3].
Table 1.1 – Examples of intravital imaging systems used in preclinical cancer research [3]
Imaging System Type of Wave Penetration Depth Resolution
OCT Near-Infrared Light 1 −2 mm 1 −10 µm
US Sound Centimeters 50 µm
Fluorescence Microscopy Visible or Near-Infrared Light 300 −800 µm 1 −5 µm
MRI Radio-Frequency Wave No Limit 10 −100 µm
CT X-Ray No Limit 50 µm
PET γ-Ray No Limit 1 −2 mm
1.1.1 OCT versus Microscope Slides
A pilot study by Fine et al. [4] included the comparison of tissue samples which were scanned with spectral-domain OCT (SD-OCT) and the examination with the microscope of corresponding areas. Both methods showed similar results for fallopian tube samples, so that soft tissue, fat, and blood vessels could be differentiated. Examining breast tissue, areas of inflammation and tumor fragments could be differentiated from normal fat tissue in microscopic and SD-OCT images. Even though some difficulties occurred correlating images that were scanned with SD-OCT with those created by the microscope, tissue structures could be differentiated in both cases. Thus, three-dimensional (3D) SD-OCT images can be used for diagnosing pathologies so that obtaining microscope slides can eventually be bypassed [4].
1.1.2 OCT versus Micro-CT System
Sun et al. [5] compared 3D X-ray images created with a commercial micro-computed tomography (micro-CT) system with those created by a SD-OCT system. It could be seen that microcalcifications and density variations within the tissue were prominent on the X-ray images, in order to differentiate fat and tumor tissue types. The OCT system provided information on the tissue’s morphology and it was also possible to identify all different tissue types. After image processing using different kinds of image filters and kernel sizes, the border of the tumor could be distinguished clearly in the OCT image, while the grainy micro-CT image made a distinction more difficult [5].
1.2
Optical Coherence Tomography (OCT)
OCT is an optical imaging technique based on light interference performing 3D cross-sectional imaging by adjacent depth-scans and can be used in a non-invasive way. With an imaging depth of 2 −3 mm and a spatial resolution of up to 1 µm, OCT is used to image microstructural biological tissue in medical fields such as ophthalmology, dermatology, gastrointestinal endoscopy, interventional cardiology, dentistry, gynaecology, and laryngology [1].
The initial application of OCT was on the retina, as it is transparent and has low scattering properties. Tissue optical properties, such as scattering and absorption, lead to optical attenuation which limits the imaging depth for most biologic tissues. OCT imaging principles show similarities to US, but it uses light
instead of sound. In this case, the propagation velocity at approximately 3 × 108m/s is around a million
times faster. The distances within the tissue are determined by the measurement of the echo time delay of back-scattered and back-reflected light, which can be measured using low-coherence interferometry. In this way, optical properties of micro-structural tissues can be obtained showing morphological and cellular features [2].
1.2.1 Imaging Principles
Spatial information can be obtained by the echo time delay ∆T , where z reflects the travelled distance and v is the velocity of the light wave:
∆T = z
v (1.1)
A generic OCT system measures structures with a resolution of 10 µm corresponding to a time
resolution of approximately 30 × 10−15s. As direct electronic detection is not possible, interferometry or
correlation techniques have to be used [2]. An optical fiber-based Michelson interferometer, as shown in Figure 1.2, is used to explain the basic concept of a generic OCT system. Light from a low-coherence source is evenly split by a fiber coupler and redirected into a reference arm and a sample arm. The use of low-coherence light leads to the observation of the interferometric signal over only a limited depth in the sample. The reference arm uses a reference mirror to reflect light back into the interference system, where a specific reference delay is created. At the sample arm, the light beam is focused on the sample, where light influenced by the material and its refractive indices is scattered back into the optical scanning system. At the fiber coupler, both reflected light beams are recombined, so that a specific interference
pattern is generated. This intensity pattern is recorded by a photodetector at the fiber output of the interferometer [1, 2, 6].
By using a low-coherence light source, interference of both reflected light beams can only occur when the traveled path length matches the coherence length of reflected light [1, 2]. To be able to measure the intensity of back-scattered light and the echo time delay from the sample arm, the interference output has to be detected and demodulated, while the path length of the reference arm has to be measured [2]. When the distance along the light propagation of the reference mirror changes, the back-scattered light forms an interference pattern with corresponding depths of light reflected from the sample. Therefore, light intensities and their dependency on the depth in the sample can be obtained [1].
Figure 1.2 – Schematic setup of a generic Michelson interferometer based on a low time-coherence light source as referred to [6].
An A-scan or depth-scan can be obtained by recording the signals at the detector for a complete movement cycle of the reference mirror which represents the reflectivity profile over the sample’s depth. Multiple adjacent A-scans compose a two-dimensional (2D) cross-sectional or B-scan. A 3D OCT dataset can be provided by obtaining different subsequent B-scans [1, 6].
To create an axial 2D image profile that gives information on the sample depths, the intensity of B-scans following one another in transversal direction is measured. As the axial resolution is independent of the focusing conditions of the light beam, high resolution can only be obtained by the coherence length of the light source. The detected interference signal is equivalent to its electric-field auto-correlation, where the coherence length represents its spatial width. Results that are equivalent to the envelope of the electric field auto-correlation can be obtained by calculating the Fourier transform (FT) of the
power spectrum [2]. The axial resolution ∆z, which is equivalent to the coherence length lc, is inversely
proportional to the power spectrum width [2] and the spectral bandwidth of the source. If the spectral distribution of the light source is Gaussian, it is defined by [1, 2, 6]:
∆z = lc=
2 ln(2) π
λ20
The mean or center wavelength is denoted by λ0 and ∆λ represents its bandwidth or spectral width [1, 7]. To achieve a high axial resolution, optical sources with a broad bandwidth can be used [1, 2].
For an OCT image, the minimum focused spot size of the light beam determines the transverse resolution ∆x, which is independent from the axial resolution [1, 2]:
∆x = 4λ0 π f d = 1.27 λ0 N A (1.3)
In this equation, d refers to the spot size of the light beam projected on the objective lens, and its focal length at the sample arm is given by f . A high transversal resolution can be achieved by focusing the light beam to a minimum spot size, which shows an inverse proportionality to numerical aperture (NA) [1, 2]. ∆x also shows a proportionality to the confocal parameter b which is also defined as depth of focus [2]:
b = π∆x
2
2λ (1.4)
At the focal point, a bigger diameter of the light beam can be provided using a low NA, where the depth of field increases. In the opposite scenario, a high NA leads to a small focused diameter; therefore, it shows an excellent transversal resolution. A low NA is used in most commercial OCT systems to achieve a depth range of some millimeters, leading to a transversal resolution of around 20 −25 µm [1].
To calculate the signal-to-noise ratio (SNR), standard techniques can be used, where the quantum efficiency of the detector is represented by η, the detected power is given by P , hv is the photon energy, and N EB represents the detection noise’s equivalent bandwidth. The optical power has to be increased if a high imaging resolution or acquisition speed is used to keep a given SNR [2].
SN R = 10 log ηP
2hvN EB (1.5)
Light sources with short pulse lasers which have an extremely high power output and short coherence length are used to create high-speed and high-resolution images for research applications [2].
1.2.2 Detection Principles
Mainly used OCT technologies can be divided into time-domain OCT (TD-OCT) and Fourier-domain OCT (FD-OCT), where the latter one can further be distinguished into swept-source OCT (SS-OCT) and
SD-OCT [1, 6]. The photocurrent ID is captured by a detector including the ’DC term’, a pathlength
independent offset, the ’cross-correlation term’ and the ’auto-correlation term’. The cross-correlation component is desired for OCT imaging as they show a proportionality to the reflectivity of the sample. The auto-correlation component represents the interference between different sample reflectors referring to artifacts in a common OCT system design‘[6].
Time-Domain Low Coherence Interferometry
TD-OCT uses a low-coherence, broadband and continuous-wave [6] light source. It is based on a scanning reference delay, which can be measured due to the known path length back-reflected light travels within the reference arm [1]. Using this method, an A-scan can be obtained due to the fast variation of an optical delay line [8]. Depth and structural information can be obtained by the interference pattern that is created
when back-reflected light from the reference and sample arm travelled an equal optical distance. Light echoes can be detected sequentially when the reference mirror is moved stepwise [1].
Fourier-Domain Low Coherence Interferometry
FD-OCT systems have a static reference arm, where the reference mirror is not moving. This method can achieve a higher data acquisition speed in comparison to TD-OCT systems. Two primary methods evolved from the FD-OCT system, SD-OCT and SS-OCT. SD-OCT uses a charge coupled device (CCD) line detector as all depth components are obtained simultaneously. To extract the intensity information as a function of depth, FT is performed on the interference pattern that was detected by the camera. In SS-OCT, a broadband source is used to sequentially scan over the bandwidth of the source. Therefore, a point detector is used to subsequently detect the back-scattered light from the sample. An A-scan is obtained by performing FT on the detected signal [1].
1.2.3 Resolution and Field of View
Figure 1.3 – Schematic setup of a generic sample arm in an OCT system describing both axial and lateral resolution as well as field of view (FOV), related to [6].
An OCT system, shown in Figure 1.3, can be seen as reflection-mode scanning confocal microscope, where the optical fiber is used as the pinhole aperture for light illumination and collection. Considering an ideal reflection of a confocal microscope, the intensity of a point reflector in the focal plane in lateral position is given by [6]:
I(v) = 2J1(v)
v
4
(1.6)
where J1(v) is a first order Bessel function, and v is the normalized lateral range parameter given by
v = 2πx sin(α)/λ0. In an imaging system, the smallest distance were two points in space can still be
distinguished is given by the spacial resolution. The full width of an OCT system of the lateral resolution δx and axial resolution δz at half the maximum power can be defined by [6]:
δx = 0.37 λ0
sin(α) = 0.37
λ0
N A (1.7)
where λ0defines the center wavelength of the light source, α represents half of the objective’s angular
optical aperture and NA can be assumed from N A = sin(α). The lateral FOV is given by [6]:
where f is the frequency used and θmaxdefines the maximum one-sided scan angle. Describing the response of an ideal confocal microscope, the detected intensity from a planar reflector is given by [6]:
I(u) = sin(v/2)
u/2
2
(1.9)
where u defines the normalized axial range parameter given by u = 8πz sin2(α/2)/λ0with z being
the distance to the optical axis. Due to that, the axial resolution is unaffected by the imaging depth as the distinguishable points are parallel to the light beam. The lateral resolution on the other hand gives information about the distance, where two points perpendicular to the light beam can be distinguished. The most commonly used OCT applications in research and clinics use an objective with a low NA, so that
the lateral resolution δx matches approximately the axial resolution δz = lc, see Equation (1.2), which is
defined by the interferometer. This allows imaging with nearly isotropic resolution, where the confocal gate length limiting the axial range shows a much higher value than the lateral resolution. The axial FOV is given by [6]: F OVaxial = 0.565λ sin2α 2 = 0.565λ sin2hsin−12(N A)i (1.10)
where λ represents the used wavelength.
1.2.4 Artifacts and Noise
Artifacts from the ’DC term’ and ’auto-correlation term’ in FD-OCT can be removed by the use of phase modulators at the reference arm and a sequentially acquired phase-shifted spectral interferogram with increasing acquisition time. For optical detection, shot noise is the fundamental limitation, where additional noise sources have to be considered as well [6].
Time-Domain OCT
The SNR for a TD-OCT system with the assumption of a single sample reflector and neglection of the ’auto-correlation term’ is defined by [6]:
SN RT DOCT = hIDi2 T DOCT σT DOCT2 = ρST DOCTRS 2eBT DOCT (1.11) hIDi2
T DOCT is the mean square of the detected photocurrent, where the pathlength of the reference
arm zRand the sample arm zSis equal. In the equation, σ2T DOCT represents the noise variance, ρ the
responsivity of the detector, ST DOCTRSis the power reflected from the sample, e the electric charge and
BT DOCT defines the detection bandwidth. It can be assumed that the amount of reflected light of the
reference is much higher than the one reflected from the sample so that the mean detector photocurrent is dominated by the power of the reference arm [6].
Fourier-Domain OCT
A spectrometer based FD-OCT shows a big advantage in SNR over TD-OCT with an inherent connection between spectrometer-based and swept-source systems. Assuming a single sample reflector and no autocorrelation terms, the general SNR is given by [6]:
SN RF DOCT = hiDi2F DOCT σ2 F DOCT = ρSF DOCT[km]RS 4eBF DOCT M (1.12)
where hiDi2F DOCT is the mean-square peak signal power and σF DOCT2 is the noise variance. The
noise variance is uncorrelated in each spectral channel and in the inverse discrete Fourier summation the variances add incoherently by assuming that the power reflectivity from the reference is much bigger than
from the sample. SF DOCT[km] is the power reflected from the sample that corresponds to the detector’s
spectral channel m, giving a higher peak signal power than the interference patterns which are coherently added. In FD-OCT, interference is sensed in a longer coherence length in each detection channel as in
TD-OCT with a single detection channel. BF DOCT defines the detection bandwidth [6].
For a better comparison, Table 1.2 was created, including following assumptions: equal A-scan length, acquisition time, instantaneous sample arm power, and a rectangular shaped source spectrum. The illumination power in each spectral channel is equal to the total illumination power in a TD-OCT system as
only one channel at a time is illuminated in an SS-OCT system, and therefore SSSOCT[km] = ST DOCT.
An SD-OCT system illuminates and detects all spectral channels simultaneously; therefore, the allowable
power per spectral channel has to be decreased with a factor M , so that SSDOCT[km] = ST DOCTM .
In SS-OCT, the detection bandwidth is limited by the analog-to-digital sampling frequency, giving
BSSOCT = BT DOCT. In SD-OCT the signals of each channel are integrated, so that the detection
bandwidth BSDOCT = BT DOCTM . In both FD-OCT methods the whole depth range is sampled all the time,
so that the SNR improvement is increased by a factor M , and decreased by a factor of 2 due to the fact that positive and negative sample displacements to the reference position are generated. In this way [6]:
SN RSDOCT = SN RSSOCT + ρST DOCTRS 4eBT DOCTM = SN RT DOCT M 2 (1.13)
T able 1.2 – Comparison of SNR expressions between TD-OCT , SD-OCT and SS-OCT , which is normalized to the instantaneous po wer ST D O C T of the sample arm and the detection bandwidth BT D O C T used in TD-OCT [6 ] TD-OCT SS-OCT SD-OCT Mean-Squar e P eak Signal P o wer hiD i 2 ρ 2S 2 TD O C T 2 [R R RS ] ρ 2S 2 TD O C T 4 [R R RS ]M 2 ρ 2S 2 TD O C T 4 [R R RS ] Noise V ariance σ 2 ρS T D O C T RR BT D O C T ρS T D O C T RR BT D O C T M ρS F D O C T [km ]R R BF D O C T M Signal to Noise Ratio S N R = hiD i 2 σ 2 ρS T D O C T RS 2 eB T D O C T ρS T D O C T RS 2 eB T D O C T M 2 ρS T D O C T RS 2 eB T D O C T
1.3
Medical Application
Pathological tissue can be defined as abnormal behaviour of healthy cells leading to cancerous tumor growth [9]. Benign tumor cells are non-cancerous and do not spread over the whole body as they are limited in growth. Malignant tumor cells are cancerous and grow beyond control, in a way that they can spread over the whole body [10]. Given that around 98% of all tumors are of epithelial origin, diagnosis of early cancer formation is essential, especially in intraepithelial stages [11].
1.3.1 Medical Use of OCT
Optical properties provided by OCT can be visualized by scattering and reflection from particles due to different refractive indices within the sample [5]. A big advantage can be seen in the non-invasive use, which can support biopsies and help guiding therapies [12]. When considering particularly pathological tissue, OCT is an important system to detect and diagnose cancerous tissue properties and blood flow [12]. As the penetration depth is restricted to approximately 2 mm, endoscopic OCT can be used to image internal sites to get more information on disease progression. OCT imaging is also used: in dermatology for diagnosing and examining skin tumors and inflammatory skin diseases in vivo [13]; in cardiology for imaging atherosclerotic plaque, coronary stents or noncoronary cardiac applications [14]; in gynecology [15]; in urology [16]; in gastroenterology [14]; and in laryngology [11]. For examining larger tumor volumes, US or CT is more suitable, even though OCT can provide a greater contrast [3].
1.3.2 Thyroid Gland Healthy Thyroid Gland
Healthy thyroid consists of two lobes giving it a butterfly shape and has a longitudinal diameter of 40 −60 mm and an antero-posterior diameter of 13 −18 mm [17]. The thyroid is surrounded by an approximately 150 µm thick capsule [18]. The main task of the human thyroid gland is the production of two main body hormones Triiodothyronine (T3) and Thyroxine (T4). The cells have the unique functionality to absorb iodine, which is used to produce the hormones. The balance of T3 and T4 is kept by the hypothalamus and pituitary in the brain by the release of thyroid stimulating hormones. The metabolism and the activity of cells are regulated by the thyroid gland, where an imbalance can lead to hyperthyroidism (overproduction of T3 and T4) or hypothyroidism (underproduction of T3 and T4) [19].
Normal sized follicles show a well organized structure within the thyroid gland, with a round or oval shape ranging from 50 −500 µm [20]. They can be seen as the functional units which contain homogeneous colloid that is surrounded by thyroid epithelial cells (thyrocytes) [21]. The connective tissue builds a framework binding together different structures. It can also store fat, separate the follicles and carry blood, nerves and lymph vessels [22]. The interchangeable amount of colloid, activity and number of cells surrounding the follicle define the follicle’s volumetric size [21], which also provides information about the gland’s activity [22]. Smaller follicles are centrally situated while the larger ones can be found in the periphery of the thyroid.
Follicle properties such as amount, inner follicular surface area and colloidal volume can provide information on synthesis, secretion and storage of thyroid hormones, as well as on functional and morphological levels and their changes [21]. By examining normal and pathological thyroid samples, Lee et al. [23] showed that OCT can overcome the limitations of US and provide information on differences
between normal and tumor thyroid follicles [23]. A mean colloid volume of 7.02 × 105µm3 for human thyroid glands could be calculated by Kot et al. [22].
Pathology of Thyroid Gland
In case of a malfunctioning thyroid, different benign and malignant diseases can be distinguished [19]. Around 95% to 99% of all thyroid nodules [19] and the majority of all diffuse thyroid diseases [24] are benign, leaving a small percentage of malignancies which are classified by the spread into adjacent or distant tissue parts [19]. Pathologies can further be specified by hormone production, increased thyroid growth, nodules or changes in follicular properties and thyroid cancer [19]. An abnormal production of the hormones T3 and T4 can be specified by hypothyroidism or hyperthyroidism, where epithelial cells can increase to an abnormal size. An inflammation of the thyroid is called thyroiditis and can cause abnormal functionality.
Goiter is an enlargement of the thyroid gland, which appears as a bulge in the neck and can be either toxic or non-toxic [19]. Zhou et al. [20] examined benign thyroid diseases and could observe follicular changes compared to healthy tissue. The multinodular colloid goiter disease showed variable shapes and sizes of the atrophic follicles, as well as a reduction in colloid and increase in cellularity [20]. Thyroid enlargement can be seen in many thyroid diseases [21], where the incidence of malignancy of the gland is usually below 5% [19]. Hashimoto’s thyroiditis leads to enlarged and distorted follicles as well as to an increase of interfollicular cellular density [20].
Graves’ disease is an autoimmune disorder, where the thyroid is attacked by proteins causing hyper-thyroidism [19]. The follicles can show diffuse volume growth (hypertrophy) and an increase of tissue (hyperplasia), which can grow into adjacent muscle cells. A lack of fibrovascular cores can be obtained in tall follicular cells with papillae, but normal follicle cells can rarely be found in small clusters in adjacent lymph cells. At a very late stage of Graves’ disease, a decrease in the amount (follicular atrophy), an increase in density (nodularity) of follicles or an excess of fibrous tissue (fibrosis) can occur [25].
Only 1% of all cancer types are thyroidal malignancies which can be specified into lymphoma, papillary, follicular, anaplastic, and medullary cell carcinoma [19]. A benign tumor such as follicular adenoma shows variations in the size of the follicle having a diameter of 40 −800 µm [20]. Papillary tumors are the most common form of thyroid cancer, occurring in 60-70% of patients, and are mostly defined by chaotic hypervascularity. The solid and homogeneous nodules in follicular cell carcinoma show a tick irregular capsule and emerge around 5-15% [19]. Normally sized and structured follicles are absent and replaced by complex papillae. A follicular variant is characterized by a homogeneous microfollicular pattern of the nodular thyroid. The microfollicles are densely packed and have a size of around 50 µm [20].
Anaplastic thyroid cancer is defined by irregular boundaries and occurs around 5-10%, while medullary cell carcinoma arises around 3-5% showing irregular margin and vascularity. Medullary carcinoma can have a very variable histomorphologic appearance, where normal follicles are absent and dense fibrosis surround tumor cells [20]. Lymphoma is quite rare with an occurrence of around 4% and is mostly defined by a heterogeneous parenchyma [24].
1.4
Optical Properties
The specification of tissue optical properties is necessary for design and interpretation of optical measure-ments, so that measurements can be interpreted and therapies can be planned correctly.
1.4.1 Refractive Index
When light hits another medium, the amount of rays reflected and refracted depends on the dimensionless refractive index n(λ) of the specific medium. The refractive index is wavelength dependent and can be seen as function of the wavelength λ. The absolute refractive index is defined by [26, 27]:
n(λ) = λ0 λm = c/f v/f = c v (1.14)
where the wavelength of light is given by λ0 in vacuum and λmin the medium. The wave frequency
is defined by f , the velocity of light in vacuum by c and in the medium by v. The refraction of light can be described by Snell’s law [26]:
n1sin(θ1) = n2sin(θ2) (1.15)
where the refractive index n1and the incident angle θ1indicate the parameters of the medium from
which the light comes, and the refractive index n2and the incident angle θ2describe the parameters of the
medium, where the light ray travels into [26].
1.4.2 Absorption Coefficient
The intensity of light travelling in a medium decreases as part of the incident light is absorbed over the path length. Applying the Beer-Lambert law, the transmission T is defined by [28]:
T = I
I0
= 10CL= exp(−µaL) (1.16)
where I0is the initial light intensity, is the molar extinction coefficient, C the concentration, L the
pathlength and µathe absorption coefficient defined by µa= 2.303C. Depending on the wavelength,
absorption is normally dominated by the chromophores water and blood, where other absorbers as fat, melanin, bilirubin and beta-carotene also have to be considered in the tissue (see Figure 1.4). For the region around 1300 nm the main absorber is water; however, this absorption is quite low [28].
1.4.3 Scattering Coefficient and Anisotropy
For thin tissue samples (around 100 nm), the scattering coefficient can be measured by a collimated
transmission measurement Tcon the scale of one mean free path (mf p = 1/µs), so that µs= − ln(Tc)/L,
where L specifies the thickness. Anisotropy is difficult to measure as a direct measurement of the scattering function p(θ) involves angular light scattering, where θ describes the scattering deflection angle. By
using µsof collimated transmission measurement and the more robust µ
0
sof diffuse light measurement,
anisotropy can be calculated by g = 1 − µ0s/µs. Having a value of 1, light would scatter completely
Figure 1.4 – Absorption coefficients spectrum of the main chromophores in human tissue showing protein, melanin, oxy- and deoxy-hemoglobin, collagen and water at different wavelengths, based on [29].
reflectance R depends on the attenuation µt(cm−1) and absolute reflected light ρ which were used to
specify the scattering coefficient and anisotropy [28].
Tissue optical properties include the absorption coefficient µa(cm−1), the scattering coefficient
µs(cm−1), and the scattering function p(θ, ψ)(sr−1), where ψ describes the scattering azimuthal angle
(see Figure 1.5), while the refractive index is defined by n. During thin tissue transmission microscopy or confocal reflectance microscopy, including OCT, only few or a single scattering event occurs, where the scattering function is appropriate. For thick tissue, anisotropy g = hcos(θ)i characterizes the scattering direction [28].
Figure 1.5 – After a photon is scattered different optical angles arise due to the change of direction, shown in blue. The scattering deflection angle θ is shown in green while the scattering azimuthal angle ψ is denoted in orange.
1.4.4 Reduced Scattering Coefficient
The reduced scattering coefficient µ0sis sufficient for characterizing skin, bone, brain, breast as well as
other fibrous, soft and fatty tissue [28]. Diffusion of light along a distance using one big step having
isotropic scattering is described by the reduced mean free path mf p0 = 1/µ0s[cm]. Equivalently, it can
angle θ, where a great amount of scattering events have to occur before an absorption event, i.e. µa µ
0
s (Figure 1.6). The reduced scattering coefficient is calculated by [30, 31]:
µ0s= µs(1 − g)[cm−1] (1.17)
Figure 1.6 – The reduced mean free path mpf0is shown in blue with one big step and the mean free path mpf with ten small steps in black [31].
1.4.5 Tissue Optical Properties for a 1300-1350nm OCT Laser Source
Bashkatov et al. [32] studied the optical properties of cortical bone samples in vitro in a spectral range of 800 −2000 nm. Their results for the measured absorption coefficient in the spectral range of 800 −1000 nm agree with the results of other authors ([33–37]) [32]. Results of measured refractive indices for human tissue specimen performed by Tearney et al. [38] at λ = 1300nm are shown in Table 1.3.
Table 1.3 – Obtained values for the refractive index of different human tissue specimens at λ = 1300nm by Tearney et al. [38].
Tissue Measured Refractive Index n
Epidermis (in vivo) 1.34
Heart Muscle (in vitro) 1.382
Dermis (in vitro) 1.40
Mesenteric Adipose (in vitro) 1.467
Stratum Corneum (in vivo) 1.51
Changes in the refractive index at the boundaries and a multi-layer structure of the sample are taken into account in iterative methods using a theoretical model for light propagation. Estimating optical coefficients can lead to errors when analyzing optical parameters in different experiments. Factors such as pathology of the tissue, geometry of irradiation, interface of refractive indices and angular resolution of the photodetector have to be considered. Obtained values in Table 1.4 are collected from different authors presented in [39]. The data output strongly depends on the used measurement technique, preparation of specimen, systematic errors and instrumentation noise [39].
T able 1.4 – T issue optical properties: absorption coef ficient µa , reduced scattering coef ficient µ 0 ,s and anisotrop y g of human tissue samples using a w av elength λ of 1300 − 1350 nm [32 , 39 , 40 ]. T issue W a v elength Anisotr opy Absor ption Coefficient Scattering Coefficient Refer ence λ [nm ] g µa [cm − 1 ] µ 0 [1/s cm ] Cortical Bone (in vitro) 1325 (fix) 0.9 0.6 15 [32 ] Stomach W all Mucosa 1300 0.84 2 7 [40 ] Caucasian Skin 1300 -0.41 14.7 [39 ] Caucasian Skin 1300 -1.76 6.6 [39 ] Caucasian Skin 1300 -1.77 10.69 [39 ] Epidermis 1300 -0.71 25.7 [39 ] Dermis 1300 -1.19 16.2 [39 ] Subcutaneous Adipose T issue 1300 -0.89 7.81 [39 ] Subcutaneous Adipose T issue 1300 -1.05 15.8 [39 ] Sclera 1300 -1.58 16.47 [39 ] White Brain Matter 1300 -2.32 24.35 [39 ] Gray Brain Matter 1300 -1.75 4.33 [39 ] T umor (Glioma) 1300 -2.62 8.88 [39 ] Cranial Bone 1300 -0.49 14.92 [39 ] Maxillary Sinuses Mucous Membrane 1300 -0.67 3.90 [39 ] Stomach W all Mucose 1300 0.833 1.76 6.67 [39 ] Gallbladder T issue 1350 -6.95 4.56 [39 ] Gallbladder Bile 1350 -3.99 2.94 [39 ] Aorta 1300 -4.59 9.47 [39 ]
1.5
Image Processing
1.5.1 Speckle and Speckle Noise
Speckle is caused by a limited bandwidth and is a basic property of images and signals for detection systems using a narrow bandwidth. It mainly occurs due to the influence and change of spatial coherence caused by multiple back-scattering and random delays from multiple forward scattering. Therefore, parameters such as the coherence light properties of the source, multiple scattering, the light beam phase aberration and the detector’s aperture can influence the speckle. In OCT images, these changes in shape of the wavefront lead to destructive and constructive interference in localized regions, speckle [41].
Speckle noise is the pure distortion preventing an ideal image [42], as it leads to reduced contrast and in highly scattering structures makes it difficult to distinguish boundaries [41]. Multiple scattered photons lead to an interference generating speckle noise which shows a grainy appearance and low-intensity features [42, 43]. If all back-scattered light in an OCT image show constructive interference, speckle noise will vanish resulting in a highly improved image [41]. For structural enhancement and image denoising in an ordinary OCT image speckle noise suppression is needed [43]. The speckle pattern does not depend on depth [41] and is not correlated to optical wavelengths, positions, or angles. Therefore, speckle can only be reduced in proportion to the amount of acquired data sets but cannot be removed from generic OCT images [43].
1.5.2 Filtering and Denoising Methods
Speckle suppression methods can be differed into digital denoising algorithms and frame averaging methods [42]. The latter technique can further be divided into single-frame and multi-frame averaging. A single-frame averaging technique utilizes denoising methods such as Wiener filter or wavelets which often includes assumptions on an a priori model for noise and signal. Using multi-frames as denoising technique, repeated B-scans from the same unique position are required to increase the SNR [44]. In the OCT system parameters such as the frequency and incident angle of the laser beam, or the recorded angle of back-scattered light can be changed to decrease speckle noise in an even higher level. This means, that multiple frames from different positions of the sample are generated [42].
For the usage on OCT images, multiple digital denoising methods have been adapted including median filters as well as simple averaging, bayesian estimations, diffusion filters, complex diffusion or wavelet thresholding [42]. The adaptive Wiener filter is one of the most commonly used low-pass filters, where successful reduction of speckle noise could be reported for OCT images [45]. Enhancement in image contrast can be achieved by filtering techniques including rotating kernel transform which however leads to edge blurring. Techniques based on wavelet filtering show good noise reduction while maintaining image sharpness [43]. Wavelet denoising algorithms can include wavelet decompositions as the complex dual-tree discrete wavelet transform (WT) (WT), adaptive wavelet filters, or curvelets transformations [42]. It has been confirmed that filtering methods based on WT are highly suited for reduction of noise for improvisation and achieve enhanced images [42, 43, 46, 47].
Discrete Wavelet Transform
Applying discrete WT (WT) leads to a 2D decomposition of the image into four different subband images, low-low (LL), high-low (HL), low-high (LH) and high-high (HH). The subbands were down-sampled,
where the even indexed columns and rows are kept [48]. The wavelet decomposition generates one approximation coefficient and three detail coefficients including horizontal, vertical and diagonal edge information. The WT is limited, as it is not shift-invariant and shows problems with orientation distinction in multiple dimensions. These limitations can be overcome by cDWT, that is motion selective and the six subbands at each scale are spatially oriented in distinct angles [49].
Stationary Wavelet Transform
The transformation of a multi-level 2D stationary WT (WT) is similar to DWT, where the subbands are equal to the size of the input image as no down-sampling is performed [48]. Adler et al. [43] adapted and enhanced a wavelet thresholding method by Chang et al. [47] for spatial OCT image properties. The wavelet-based, spatially adaptive speckle reduction algorithm was demonstrated on ophthalmologic TD-OCT and FD-OCT images [43].
The Gibbs ringing artifact makes the image reconstruction more difficult, as it appears as spurious signal close to sharp edges. These artifacts could be removed by the application of an undecimated WT resulting in improved denoising. For each level the vertical, horizontal and diagonal details were stored in individual subbands giving the specific edge information of the image. All subbands however included an evenly spread information of speckle noise. The wavelet decomposition gives detail coefficients which can be estimated into signal (speckle) or speckle noise information as true edge features show spatial and multilevel clustering. Adapting the threshold which is applied on the coefficients, the noise can be decreased while the sharpness of the edges is kept [43].
For the image quality evaluation, different regions of interest (ROIs) were chosen over the image, where the average contrast-to-noise ratio (CNR), equivalent number of looks (ENL), SNR as well as the sharpness were measured. ENL is the measure of smoothness of homogeneous regions which is calculated over defined ROIs by measuring the smoothness of the image. The contrast between the background and image feature is determined by the CNR. The image quality metrics are defined by [43, 45]
SN R = 10 log 10mean(I 2 lin) σ2 lin (1.18) CN R = 1 A A X a=1 µa− µb q σ2 a+ σ2b (1.19) EN L = 1 H H X h=1 µ2h σ2h (1.20)
where Ilindescribes the linear image intensity and σlin2 the variance of noise in Ilin. The mean and
variance of the same background area are described by µband σbwhile µaand σaas well as µhand σh
describe the mean and variance of all (A) and only homogeneous (H) ROIs, respectively.
Adler et al. [43] and Ozcan et al. [45] showed that wavelet denoising could improve the SNR, CNR and ENL, while the sharpness of the image slightly degraded compared to the original or Gaussian filtered image [43, 45].
Multi-Frame Wavelet Transform
A multi-frame denoising method was proposed by Mayer et al. [42], where processing of the single frames was done before averaging, and for wavelet denoising it was possible to differ between structural and noise information in the tissue images. For the wavelet transform different wavelets such as ’haar’, ’db8’ or ’sym8’ could be chosen for SWT and ’dualTree’ for cDWT. The single steps include scaling of the single intensity frames with logarithmic transformation followed by wavelet decomposition, weight computation, coefficient weighting, averaging and wavelet reconstruction to obtain the final image. Weighting is performed by either only significance or correlation weight or a combination of both. Significance weight is used for a local noise estimation on the detail coefficients while the presence of structural information in the image frame is provided using the correlation weight on the approximation coefficients. For the image reconstruction, an inverse WT is performed on the averaged coefficients. At least two input frames were required for this filtering method, where an increased amount enhanced the denoised output image [42]. Evaluating the multi-frame wavelet denoising method both the SNR gain as well as the edge integrity was calculated and averaged over different ROIs using [42]
SN Rgain=
σ(Na)
σ(Nf)
− 1 (1.21)
where the noise is estimated by the difference of the gold-standard image Fg which was created by
averaging 455 recorded frames obtained from 35 positions and the original image by σ(Na) ≈ Fa− Fg
and σ(Nf) ≈ Ff − Fg. Faand Ff refer to the averaged and filtered image, respectively [42].
Other Filtering Methods
For denoising volumetric SD-OCT images, Fang et al. [44] presented a new method, multiscale sparsity based tomographic denoising. Using this algorithm, a non-uniform scanning pattern is applied, where few B-scans are captured slowly at a high SNR and the others fast at a nominal SNR. The high-SNR scans are used to denoise the adjacent low-SNR images, as they are expected to show similar noise pattern and texture. A volumetric ocular scan was generated with a densely sampled large FOV having a low-SNR and other scans were taken which were sparsely sampled but had a high-SNR. Low-SNR images were denoised by applying a dictionary from a neighboring averaged high-SNR image. Equally sized patches were extracted from the original, up- and down-sampled image to create a multi-scale structural dictionary which was used to estimate denoised patches to denoise the SD-OCT image [44].
1.5.3 Enhanced Resolution Imaging
The super-resolution technology is a method, where a high-resolution image can be recovered using a single or multiple low-resolution images. The signal-processing technique increases the pixel number for the reconstruction. A modified back projection super-resolution method on ophthalmologic OCT images was applied by Uji et al. [50] to enhance the image quality [50]. They applied the bicubic interpolation algorithm to fourfold magnify a low-resolution image, where the new pixels are calculated from the surrounding ones creating a smoother and enlarged image. After reducing the magnified image again by a factor of 4, the difference between the processed and original image was constructed by subtraction of the minified image. To generate the enhanced-resolution image (ERI) image, the processed image was
converted from 32 to 8 bits, magnified via bicubic interpolation and added to the magnified interpolated original image [50].
The parameters detailed variance (DV) and background variance (BV) were used to compare the quality of the original and processed image. The image was separated into two intensity regions, where DV represented the high intensity containing information about details and edges and BV corresponded to the lower intensity region including the background. The DV/BV ratio increased for ERI processed images as enhanced edges led to a higher DV value, whereas the BV value only showed little changes. Additionally, mean opinion score was used for a subjective comparison by experienced specialists between the enlarged original image without interpolation, the interpolated image as well as the generated ERI image. It was shown that ERI could achieve a higher grading than the interpolated image. The results pointed out that ERI is a simple and low computation method which obtains image enhancement and gives better results as pure magnification methods with or without bicubic interpolation. It is also mentioned that the main concept of this post-processing method does not include an increased image resolution, but an image enhancement [50].
1.5.4 Image Segmentation
Three different algorithms to segment OCT and X-ray images were developed and compared by Sun et al. [5] using MATLAB. Intensity-based segmentation was better suited for X-ray images, as the contrast between adipose and tumor region is greater than for OCT images. Spatial frequency filtering is based on a band-pass filter in the Fourier domain and was optimized by adapting the center frequency as well as the bandwidth. The idea behind this filtering method was the increased density in the tumor tissue with different spatial frequency contents in relation to the adipose parts. The texture-based algorithm gave the best results for the segmentation in the OCT image, as the texture and entropy of adipose and tumor tissue show different properties. By adapting the correlation distance and used kernel size, the tumor region could be clearly separated from other tissue parts. Limitations of all codes occur for automated segmentation as universal settings for different images is very difficult to obtain, so that individual optimization provided the best results [5].
In addition to the use in medical applications, entropy segmentation was used by Satorres Martinez et al. [51] on transparent objects with an irregular surface to detect defects. Nguyen et al. [52] used entropy calculations to detect cracks in images of a pavement surface. By enhancing entropy based segmentation from a simple bi-level threshold to a multi-level threshold, Arora et al. [53] could achieve better results for sonar image segmentation also due to a better pattern recognition. Using maximum entropy combined with multilevel thresholding, Horng [54] could achieve good results in pixel extraction. This demonstrated that entropy based segmentation methods provide good performance in different fields, as it is based on the measurement of variance and randomness within the image [55].
Chapter
2
Materials and Methods
2.1
The OCT System
2.1.1 Hardware
All optical images were taken with the Telesto II (Thorlabs GmbH, Germany) OCT system out of the Telesto Series SD-OCT Imaging Systems. The system included a high resolution base unit with a center wavelength of 1300 nm, a bandwidth of 100 nm and had a maximum imaging depth of 3.5 mm. The axial resolution in air and water was 5.5 nm and 4.2 µm respectively, while the lateral resolution was 13 µm. The OCT system had a FOV of 10 mm×10 mm. The laser was classified according to DIN EN 60825-1:2008 as a laser product class 1M [56].
The optical layout of the used OCT system is illustrated in Figure 2.1, where light exited the fiber from the low coherence light source, was collimated by the collimator and directed into the beam splitter. There, light was split and directed towards the reference arm and the sample arm. In the reference arm, the intensity could be adjusted via a variable aperture before it passed a dispersion compensation plate and was sent back to the beam splitter by an adjustable prism. The reflected light was combined with the collimated light from the source and moved into the sample arm. After passing two galvanometers which allow an axial and transversal scan, the light beam was focused on the sample. The light was back-reflected and back-scattered into the fiber and directed into the detector, where it was processed and forwarded to a computer for further processing [57].
Imaging artifacts such as distortion could occur due to refractive media. The image displayed differences between the optical path length of the reference and sample arm which were calculated from the refractive index and the physical path length. The measured imaging depth in a sample changed over the whole scan as the loss of imaging depth depended on the properties of the material. For samples with a variable surface interface diffraction, shadowing effects or multiple measured structures can occur [57].
2.1.2 Software and Settings
The ThorImage OCT Software Package included with the system was used to acquire, visualize and manage the OCT scans. The software run on a Microsoft Windows 7 operating system and featured 2D as well as 3D data acquisition and processing. Additional processing techniques, such as Doppler Mode or Speckle Variance Mode, were included but not used in this application [58].
Figure 2.1 – Setup of the used SD-OCT system: A light source sends light from a FC/APC (fiber channel / angled physical contact) fiber to a collimator and further to a beam splitter. There, light is split into the reference arm and the sample arm. In the reference arm it is possible to manually adjust the reference intensity with a variable aperture before light passes a dispersion compensation plate and is re-reflected by an adjustable prism. Light which was sent directly to the sample arm is redirected by two galvanometer mirrors and will be focused by an objective lens into the sample. The back-reflected and back-scattered light from the sample is sent back into the fiber and further to a CCD sensor.
3D scans were created using the parameter range as listed in Table 2.1. Preprocessed noise reduction could be achieved by averaging eight B-scans for the image generation. The images themselves were stored as .oct files and exported to another computer for further processing.
Table 2.1 – Used parameters for 3D scans, where X refers to the amount of A-scans in a B-scan, Y gives the number of B-scans in a volume and Z is the amount of depth pixels in an A-scan.
X-Axis Y-Axis Z-Axis
Sizes 100-500 100-200 422 [pixel]
FOV 0.6-2.5 0.6-2 1.5 [mm]
Spacing 5 5 3.55 [µm]
2.2
Tissue Specimen
For image processing and examination, different types of pathological human thyroid samples were used. The Department of Surgery at Linköping University Hospital provided the specimen from patients undergoing a thyroid surgery, which were examined using the OCT system. The regional ethic committee approved the study and patients signed a written consent (Dnr 2014/452-32 and Dnr 2015/463-32). After excision, 4% formalin was used for fixation and the specimens were stored in a fridge to prevent fast autolysis and inhibit the decay of the tissue samples to preserve it for future studies. To scan the 0.5 −1 cm
sized tissue samples, the specimen was placed in a culture dish and covered by a thin saran foil to prevent specimens from oil contamination for an increased refractive index matching. In this way, the imaging resolution could be increased together with a spacer as the difference between the refractive index of oil and tissue is smaller than for air and tissue. In this project only specimen from one patient was used which was scanned at different regions.
2.3
Image Processing
Even though 3D images were obtained by the OCT system, image processing was only performed on 2D images. Therefore, individual B-scans were extracted from the volume. For the performed filtering method multi-frame WT, based on Mayer et al. [42], a multiple frame image was required. Due to that, ten 2D images of the thyroid sample were obtained from the same position, based on the assumption that the image intensity stays consistent while noise changes.
2.3.1 Data Implementation and Visualization
For better usage of the OCT data files, some MATLAB scripts were provided by the system manufacturer for reading 2D images. MATLAB R2016a was used for further processing of the saved .oct files. The provided MATLAB scripts were adapted, so that also volumetric scans could be extracted. For better illustration and a more user-friendly application, the obtained images were further used in a graphical user interface (GUI), where additional perspectives, enhancement and segmentation possibilities could be chosen.
Provided 2D Data Read-Out
The Main script included various MATLAB functions to extract and plotted the obtained .oct intensity data, showed the area of which the image was extracted, and additionally plotted the 2D image in a rotating figure. OCTFileOpen inspected the current path for the existence of the temporary folder OCTData and createed a new subfolder after removing all previous contents and subfolders. Following that, the .oct intensity file was unzipped and extracted into the new folder while all image properties and settings were saved in an additional .xml variable. By running OCTFileGetProperty, specific dataset properties as ’AcquisitionMode’, ’RefractiveIndex’, ’Comment’, ’Study’ and ’ExperimentNumber’ were displayed in the Command Window of MATLAB. The intensity values were loaded by OCTFileGetIntensity which redirected to OCTFileGetRealData. This function extracted information about the image dimensions as length, height and depth but opened, read and reshaped the intensity data in a way that only a 2D matrix was returned from a 3D image including the first B-scan.
Further MATLAB .m files such as OCTFileGetColoredData, OCTFileGetChirp,
OCTFileGetNrRaw-Dataand OCTFileGetRawData were used to read the video image, chirp vector and spectral raw data,
respectively, but were not used for additional processing. OCTFileClose was used to delete the temporary folder with all its contents.
3D Data Extraction and Implementation
To be able to extract not only 2D, but also 3D intensity data directly, OCTFileGetRealData was adapted. There, volumetric data was opened, read and reshaped including the whole data set by increasing the
dimension to 3D. The obtained volume was in form of stacked 2D scans.
GUI for Enhanced Visualization
A GUI was developed in MATLAB to examine all measured and extracted .oct files. By running the main function GUI_3Dimaging, the current working folder was opened, showing all available raw data files that could be selected. Subsequently, the GUI figure was opened and the currently empty plots and values got filled after selecting one of the previously shown .oct data files. The callback function buttonUpdate read the selected file and compared it with the currently shown one, which was empty when the GUI was just loaded. Following that, the newly selected and currently opened file were compared. The existence of a previously saved variable ’handleMY’ in the folder with the current file’s name was checked. If there was a difference between the files and the .mat file did not exist, all values in the GUI figure were set so zero and the green field, showing the currently opened file, was assigned with the name of the newly opened file.
The already provided and adapted MATLAB files OCTFileOpen, OCTFileGetProperty and
OCTFi-leGetIntensitywere called to extract information about the selected OCT data file. If the selected file
was computed for the first time, a new folder with the used file’s name was created. To decrease the computation time, OCTFileOpen was adapted in a way that already processed and stored data were loaded from the individual file folder into the MATLAB workspace without further format conversion or calculation. The extracted image intensity was normalized and the values defining spacing and size of the image were used to define the maximum value in the MATLAB GUI figure. Depending on the dimension range of the data, the file was handled in a different way.
A 3D data file had to be flipped in order to plot it the right way in the 3D slice plot which was not necessary for a 2D data file. Additionally, a meshgrid was created to show the real dimensions of the image on the x-, y- and z- coordinates. As initial value the center 2D B-scan slice was chosen in the plot. After all values were assigned and the images were plotted, the extracted image data was saved in the folder of the selected file. In this way, the processing time for further usage could be reduced.
The MATLAB GUI figure included the possibility to draw a rectangle to select more specific regions in the 2D image. Therefore, the callback function buttonRectangle was called, where an initial rectangle was plotted close to the edges of the 2D image. The user could specify the size and position of the rectangle with the mouse and could fixate it with a double click. If the drawn rectangle exceeded the range of the image, the borders were redefined to the next closest edge. In this way, error messages could be decreased and the application was orientated in a more user-friendly way.
The dimensions of the defined rectangle were extracted by clicking of the ’Cut Image’ button, where the callback function buttonCutImage was called. The obtained values were saved and used to cut the image in the defined range. The maximum size of the image was updated and in case of a 3D image, the slice plot was updated via creating a new meshgrid.
After defining the favoured 2D image, different enhancement and filtering methods such as ’Enhanced Resolution Imaging’, ’Adaptive Wiener Filter’, ’Discrete Wavelet Transform’, and ’Multi-Frame Wavelet Filter’ could be chosen, running the individual developed functions. Edge detection and segmentation was performed automatically within those functions. Additionally, the user could choose individual specifications for the threshold level and a region to delete with a maximum amount of pixels P. Each obtained filtered image as well as results of the edge detection and segmentation was automatically plotted
in a new figure and saved as .fig and .png file in the currently used directory.
2.3.2 Image Filtering and Enhancement Methods
To compare and evaluate the changes of processed 2D OCT images, different types of filters and enhance-ment methods were tested for noise reduction. The user could define nine different ROIs in the selected pathological thyroid image, where the order defining the specific regions was important. The first ROI was drawn over a background area, followed by defining five homogeneous areas and the last three ROIs should include an edge region. The user had to follow these restrictions for correct calculations of the metrics while the size and exact position of the rectangle could be freely chosen.
After fixating the drawn ROI, the color changed to inform about the specified region. A green rectangle indicated the background area while red and blue rectangles defined homogeneous regions and areas including an edge, respectively. These areas were further used for all 2D filtering methods including ERI, the adaptive Wiener filter, the DWT and the multi-frame WT.
Enhanced Resolution Imaging
Based on Uji et al. [50], ERI was applied to enhance the OCT image quality. First, the 2D image specified in the GUI was defined as original image. To create the interpolated image, resizing was applied, where the image was increased fourfold using bicubic interpolation. Following that, the magnified image was minimized by the same factor and subtracted from the original image. After converting the obtained image into an 8 bit integer image, another bicubic interpolation was applied and the output was added to the previously created interpolated image. The resulting enhanced image was of the same size as the
interpolated imageand the fourfold of the original image.
Adaptive Wiener Filter
For the adaptive Wiener filter, three different local neighborhoods of each pixel were chosen with the sizes 3x3, 5x5 and 9x9. The original image was converted into a gray-scale image using the MATLAB function mat2gray. The output image was filtered using the function wiener2 with the specification of the previously defined neighbourhoods.
Discrete Wavelet Transform
The wavelet-based image denoising method examined was based on a single-level 2D DWT [43, 49]. The original image was filtered by applying a 2D separable discrete WT (WT), real dual-tree discrete WT (WT) and cDWT. The code used for these methods was taken from [49]. The decomposition levels were individually defined by applying nearly symmetric low-pass and high-pass filters for orthogonal two-channel perfect reconstruction filter bank.
To be able to apply the different DWT methods, the image was resized, so that the X- and Y-Axis had
a size of 2x. Therefore, the image was increased to the next bigger suitable values, where new rows and
columns were filled up with zeros. Different thresholds for the WT including T = [0.5 − 5] with a step size of s = 0.5, were implemented.
Multi-Frame Wavelet Transform
The multi-frame wavelet WT based on Mayer et al. [42] was another tested filtering method. Multiple OCT frames (n=10) were taken from the same position and a denoised image was calculated by using different wavelets as basis. Based on the results in [42], only the combination of the computed significance weight and the correlation weight were applied for two SWT wavelet functions, haar and db8, as they showed increased results over the individual use. The applied method differed from Mayer et al. as the parameters other than the wavelet were kept constant, see Table 2.2.
Table 2.2 – Parameters used for multi-frame filtering method based on Mayer et al. [42], where the parameters k and p control the amount of noise reduction within the significance and correlation weight, respectively.
Parameters Values
k 1
p 2
Window Size (for median filtering) 2
Theta ’theta2’
Wavelet Function ’haar’ / ’db8’
Decomposition Level 3
Weight Mode Significance Weight & Correlation Weight
The variables k and p were used to control the amount of noise reduction for the significance and correlation weight, respectively. After assigning the variables, the image size was adapted if it was not a multiple of eight by filling up empty rows and lines with zeros. By the use of haar and db8 as wavelet functions for SWT the MATLAB function swt2 was applied, returning the detail and approximation coefficients of the image. Following that, the significance and correlation weight were individually calculated and combined afterwards. For averaging, the mean approximation coefficients were calculated and the output image was reset to its original size. For these image processing steps, framework code provided by Mayer et al. [42] was used.
2.3.3 Image Filtering and Enhancement Evaluation
For evaluating all individual enhanced output images, nine different ROIs were specified in the original image, where the size of the rectangles could be chosen individually. The first ROI specified the noisy background region (marked in green), the next five defined ROIs were drawn over homogeneous regions (marked in red) and the last three rectangles were supposed to cover a region including an edge (marked in blue). Following that, the SNR, see Equation (1.18), CNR, see Equation (1.19) and ENL, see Equation (1.20) were calculated over the defined ROIs for the original and filtered images. Image segmentation was an additional way to evaluate the different filter methods.
Edge Integrity
For evaluation of the changes over the edges of the selected image, an edge detection method, based on Ahmad et al. [59], was applied. The used method was divided into three steps including edge selection, profile extraction and curve comparison.
