Automated quantification of lung structures from optical coherence tomography images Alex M. Pagnozzi,1,* Rodney W. Kirk,1 Brendan F. Kennedy,1 David D. Sampson,1,2 and Robert A. McLaughlin1 1 2
Optical + Biomedical Engineering Laboratory, School of Electrical, Electronic and Computer Engineering, The University of Western Australia, 35 Stirling Highway, Crawley, Western Australia 6009, Australia Centre for Microscopy, Characterisation and Analysis, The University of Western Australia, 35 Stirling Highway, Crawley, Western Australia, 6009, Australia * [email protected]
Abstract: Characterization of the size of lung structures can aid in the assessment of a range of respiratory diseases. In this paper, we present a fully automated segmentation and quantification algorithm for the delineation of large numbers of lung structures in optical coherence tomography images, and the characterization of their size using the stereological measure of median chord length. We demonstrate this algorithm on scans acquired with OCT needle probes in fresh, ex vivo tissues from two healthy animal models: pig and rat. Automatically computed estimates of lung structure size were validated against manual measures. In addition, we present 3D visualizations of the lung structures using the segmentation calculated for each data set. This method has the potential to provide an in vivo indicator of structural remodeling caused by a range of respiratory diseases, including chronic obstructive pulmonary disease and pulmonary fibrosis. ©2013 Optical Society of America OCIS codes: (110.4500) Optical coherence tomography; (100.0100) Image processing; (100.2960) Image analysis.
References and links 1. 2. 3. 4.
G. Turato, R. Zuin, and M. Saetta, “Pathogenesis and pathology of COPD,” Respiration 68(2), 117–128 (2001). T. C. Kravis, A. Ahmed, T. E. Brown, J. D. Fulmer, and R. G. Crystal, “Pathogenic mechanisms in pulmonary fibrosis: collagen-induced migration inhibition factor production and cytotoxicity mediated by lymphocytes,” J. Clin. Invest. 58(5), 1223–1232 (1976). E. P. Trulock, L. B. Edwards, D. O. Taylor, M. M. Boucek, B. M. Keck, and M. I. Hertz, “Registry of the International Society for Heart and Lung Transplantation: twenty-second official adult lung and heart-lung transplant report--2005,” J. Heart Lung Transplant. 24(8), 956–967 (2005). R. Lozano, M. Naghavi, K. Foreman, S. Lim, K. Shibuya, V. Aboyans, J. Abraham, T. Adair, R. Aggarwal, S. Y. Ahn, M. Alvarado, H. R. Anderson, L. M. Anderson, K. G. Andrews, C. Atkinson, L. M. Baddour, S. BarkerCollo, D. H. Bartels, M. L. Bell, E. J. Benjamin, D. Bennett, K. Bhalla, B. Bikbov, A. Bin Abdulhak, G. Birbeck, F. Blyth, I. Bolliger, S. Boufous, C. Bucello, M. Burch, P. Burney, J. Carapetis, H. Chen, D. Chou, S. S. Chugh, L. E. Coffeng, S. D. Colan, S. Colquhoun, K. E. Colson, J. Condon, M. D. Connor, L. T. Cooper, M. Corriere, M. Cortinovis, K. C. de Vaccaro, W. Couser, B. C. Cowie, M. H. Criqui, M. Cross, K. C. Dabhadkar, N. Dahodwala, D. De Leo, L. Degenhardt, A. Delossantos, J. Denenberg, D. C. Des Jarlais, S. D. Dharmaratne, E. R. Dorsey, T. Driscoll, H. Duber, B. Ebel, P. J. Erwin, P. Espindola, M. Ezzati, V. Feigin, A. D. Flaxman, M. H. Forouzanfar, F. G. Fowkes, R. Franklin, M. Fransen, M. K. Freeman, S. E. Gabriel, E. Gakidou, F. Gaspari, R. F. Gillum, D. Gonzalez-Medina, Y. A. Halasa, D. Haring, J. E. Harrison, R. Havmoeller, R. J. Hay, B. Hoen, P. J. Hotez, D. Hoy, K. H. Jacobsen, S. L. James, R. Jasrasaria, S. Jayaraman, N. Johns, G. Karthikeyan, N. Kassebaum, A. Keren, J. P. Khoo, L. M. Knowlton, O. Kobusingye, A. Koranteng, R. Krishnamurthi, M. Lipnick, S. E. Lipshultz, S. L. Ohno, J. Mabweijano, M. F. MacIntyre, L. Mallinger, L. March, G. B. Marks, R. Marks, A. Matsumori, R. Matzopoulos, B. M. Mayosi, J. H. McAnulty, M. M. McDermott, J. McGrath, G. A. Mensah, T. R. Merriman, C. Michaud, M. Miller, T. R. Miller, C. Mock, A. O. Mocumbi, A. A. Mokdad, A. Moran, K. Mulholland, M. N. Nair, L. Naldi, K. M. Narayan, K. Nasseri, P. Norman, M. O’Donnell, S. B. Omer, K. Ortblad, R. Osborne, D. Ozgediz, B. Pahari, J. D. Pandian, A. P. Rivero, R. P. Padilla, F. Perez-Ruiz, N. Perico, D. Phillips, K. Pierce, C. A. Pope 3rd, E. Porrini, F. Pourmalek, M. Raju, D. Ranganathan, J. T. Rehm, D. B. Rein, G. Remuzzi, F. P. Rivara, T. Roberts, F. R. De León, L. C. Rosenfeld, L. Rushton, R. L. Sacco, J. A. Salomon, U. Sampson, E. Sanman, D. C. Schwebel, M. Segui-Gomez, D. S. Shepard, D. Singh, J. Singleton, K. Sliwa, E. Smith, A. Steer, J. A. Taylor, B. Thomas, I. M. Tleyjeh, J. A. Towbin, T. Truelsen, E. A. Undurraga,
#193941 - $15.00 USD
Received 22 Jul 2013; revised 24 Sep 2013; accepted 26 Sep 2013; published 9 Oct 2013
(C) 2013 OSA 1 November 2013 | Vol. 4, No. 11 | DOI:10.1364/BOE.4.002383 | BIOMEDICAL OPTICS EXPRESS 2383
5. 6. 7. 8. 9. 10. 11. 12. 13. 14. 15. 16. 17. 18. 19. 20. 21. 22. 23. 24. 25. 26. 27. 28. 29. 30.
N. Venketasubramanian, L. Vijayakumar, T. Vos, G. R. Wagner, M. Wang, W. Wang, K. Watt, M. A. Weinstock, R. Weintraub, J. D. Wilkinson, A. D. Woolf, S. Wulf, P. H. Yeh, P. Yip, A. Zabetian, Z. J. Zheng, A. D. Lopez, C. J. Murray, M. A. AlMazroa, and Z. A. Memish, “Global and regional mortality from 235 causes of death for 20 age groups in 1990 and 2010: a systematic analysis for the global burden of disease study 2010,” Lancet 380(9859), 2095–2128 (2012). H. Parameswaran, A. Majumdar, S. Ito, A. M. Alencar, and B. Suki, “Quantitative characterization of airspace enlargement in emphysema,” J. Appl. Physiol. 100(1), 186–193 (2006). E. R. Weibel, C. C. W. Hsia, and M. Ochs, “How much is there really? Why stereology is essential in lung morphometry,” J. Appl. Physiol. 102(1), 459–467 (2006). L. Knudsen, E. R. Weibel, H. J. G. Gundersen, F. V. Weinstein, and M. Ochs, “Assessment of air space size characteristics by intercept (chord) measurement: an accurate and efficient stereological approach,” J. Appl. Physiol. 108(2), 412–421 (2010). D. S. Faffe and W. A. Zin, “Lung parenchymal mechanics in health and disease,” Physiol. Rev. 89(3), 759–775 (2009). L. P. Hariri, M. B. Applegate, M. Mino-Kenudson, E. J. Mark, B. D. Medoff, A. D. Luster, B. E. Bouma, G. J. Tearney, and M. J. Suter, “Volumetric optical frequency domain imaging of pulmonary pathology with precise correlation to histopathology,” Chest 143(1), 64–74 (2013). S. Meissner, L. Knels, C. Schnabel, T. Koch, and E. Koch, “Three-dimensional Fourier domain optical coherence tomography in vivo imaging of alveolar tissue in the intact thorax using the parietal pleura as a window,” J. Biomed. Opt. 15(1), 016030 (2010). M. Mertens, A. Tabuchi, S. Meissner, A. Krueger, K. Schirrmann, U. Kertzscher, A. R. Pries, A. S. Slutsky, E. Koch, and W. M. Kuebler, “Alveolar dynamics in acute lung injury: heterogeneous distension rather than cyclic opening and collapse,” Crit. Care Med. 37(9), 2604–2611 (2009). J. Bickenbach, R. Dembinski, M. Czaplik, S. Meissner, A. Tabuchi, M. Mertens, L. Knels, W. Schroeder, P. Pelosi, E. Koch, W. M. Kuebler, R. Rossaint, and R. Kuhlen, “Comparison of two in vivo microscopy techniques to visualize alveolar mechanics,” J. Clin. Monit. Comput. 23(5), 323–332 (2009). R. A. McLaughlin, X. Yang, B. C. Quirk, D. Lorenser, R. W. Kirk, P. B. Noble, and D. D. Sampson, “Static and dynamic imaging of alveoli using optical coherence tomography needle probes,” J. Appl. Physiol. 113(6), 967– 974 (2012). B. C. Quirk, R. A. McLaughlin, A. Curatolo, R. W. Kirk, P. B. Noble, and D. D. Sampson, “In situ imaging of lung alveoli with an optical coherence tomography needle probe,” J. Biomed. Opt. 16(3), 036009 (2011). D. Cabrera Fernández, H. M. Salinas, and C. A. Puliafito, “Automated detection of retinal layer structures on optical coherence tomography images,” Opt. Express 13(25), 10200–10216 (2005). M. K. Garvin, M. D. Abràmoff, X. Wu, S. R. Russell, T. L. Burns, and M. Sonka, “Automated 3-D intraretinal layer segmentation of macular spectral-domain optical coherence tomography images,” IEEE Trans. Med. Imaging 28(9), 1436–1447 (2009). A. Mishra, A. Wong, K. Bizheva, and D. A. Clausi, “Intra-retinal layer segmentation in optical coherence tomography images,” Opt. Express 17(26), 23719–23728 (2009). A. Yazdanpanah, G. Hamarneh, B. Smith, and M. Sarunic, “Intra-retinal layer segmentation in optical coherence tomography using an active contour approach,” in MICCAI 2009, G-Z Yang, D. J. Hawkes, D. Rueckert, A. Noble, C. J. Taylor, ed. (Springer, Heidelberg, 2009), pp. 649–656. R. J. Zawadzki, A. R. Fuller, D. F. Wiley, B. Hamann, S. S. Choi, and J. S. Werner, “Adaptation of a support vector machine algorithm for segmentation and visualization of retinal structures in volumetric optical coherence tomography data sets,” J. Biomed. Opt. 12(4), 041206 (2007). S. Osher and J. A. Sethian, “Fronts propagating with curvature-dependent speed: algorithms based on HamiltonJacobi formulations,” J. Comput. Phys. 79(1), 12–49 (1988). A. Curatolo, B. F. Kennedy, D. D. Sampson, and T. R. Hillman, “Speckle in optical coherence tomography,” in Advanced Biophotonics: Tissue Optical Sectioning, V. V. Tuchin and R. K. Wang, eds. (Taylor & Francis, 2013), pp. 183–249. J. M. Schmitt, S. H. Xiang, and K. M. Yung, “Speckle in optical coherence tomography,” J. Biomed. Opt. 4(1), 95–105 (1999). A. E. Desjardins, B. J. Vakoc, W. Y. Oh, S. M. R. Motaghiannezam, G. J. Tearney, and B. E. Bouma, “Angleresolved optical coherence tomography with sequential angular selectivity for speckle reduction,” Opt. Express 15(10), 6200–6209 (2007). T. R. Hillman, S. G. Adie, V. Seemann, J. J. Armstrong, S. L. Jacques, and D. D. Sampson, “Correlation of static speckle with sample properties in optical coherence tomography,” Opt. Lett. 31(2), 190–192 (2006). B. F. Kennedy, T. R. Hillman, A. Curatolo, and D. D. Sampson, “Speckle reduction in optical coherence tomography by strain compounding,” Opt. Lett. 35(14), 2445–2447 (2010). B. F. Kennedy, A. Curatolo, T. R. Hillman, C. M. Saunders, and D. D. Sampson, “Speckle reduction in optical coherence tomography images using tissue viscoelasticity,” J. Biomed. Opt. 16(2), 020506 (2011). D. C. Adler, T. H. Ko, and J. G. Fujimoto, “Speckle reduction in optical coherence tomography images by use of a spatially adaptive wavelet filter,” Opt. Lett. 29(24), 2878–2880 (2004). J. Rogowska and M. E. Brezinski, “Image processing techniques for noise removal, enhancement and segmentation of cartilage OCT images,” Phys. Med. Biol. 47(4), 641–655 (2002). P. Perona and J. Malik, “Scale-space and edge detection using anisotropic diffusion,” IEEE Trans. Pattern Anal. 12(7), 629–639 (1990). D. C. Fernandez, “Delineating fluid-filled region boundaries in optical coherence tomography images of the retina,” IEEE Trans. Med. Imaging 24(8), 929–945 (2005).
#193941 - $15.00 USD
Received 22 Jul 2013; revised 24 Sep 2013; accepted 26 Sep 2013; published 9 Oct 2013
(C) 2013 OSA 1 November 2013 | Vol. 4, No. 11 | DOI:10.1364/BOE.4.002383 | BIOMEDICAL OPTICS EXPRESS 2384
31. P. Puvanathasan and K. Bizheva, “Interval type-II fuzzy anisotropic diffusion algorithm for speckle noise reduction in optical coherence tomography images,” Opt. Express 17(2), 733–746 (2009). 32. D. Lopes, “Anisotropic diffusion (Perona & Malik),” http://www.mathworks.com.au/matlabcentral/fileexchange/14995-anisotropic-diffusion-perona-malik. 33. C. Li, “Level set for image segmentation,” http://www.mathworks.com.au/matlabcentral/fileexchange/12711level-set-for-image-segmentation. 34. C. Li, R. Huang, Z. Ding, J. C. Gatenby, D. N. Metaxas, and J. C. Gore, “A level set method for image segmentation in the presence of intensity inhomogeneities with application to MRI,” IEEE Trans. Image Process. 20(7), 2007–2016 (2011). 35. J. Serra, Image Analysis and Mathematical Morphology (Academic Press, 1982). 36. H. Lum and W. Mitzner, “A species comparison of alveolar size and surface forces,” J. Appl. Physiol. 62(5), 1865–1871 (1987). 37. D. Lorenser, X. Yang, R. W. Kirk, B. C. Quirk, R. A. McLaughlin, and D. D. Sampson, “Ultrathin side-viewing needle probe for optical coherence tomography,” Opt. Lett. 36(19), 3894–3896 (2011). 38. L. Scolaro, D. Lorenser, R. A. McLaughlin, B. C. Quirk, R. W. Kirk, and D. D. Sampson, “High-sensitivity anastigmatic imaging needle for optical coherence tomography,” Opt. Lett. 37(24), 5247–5249 (2012). 39. A. E. Lefohn, J. E. Cates, and R. T. Whitaker, “Interactive, GPU-based level sets for 3D brain tumor segmentation,” in MICCAI 2003, R. E. Ellis, T. M. Peters, ed. (Springer, Heidelberg, 2003), pp. 564–572. 40. R. T. Whitaker, “A level-set approach to 3D reconstruction from range data,” Int. J. Comput. Vis. 29(3), 203– 231 (1998). 41. M. Kass, A. Witkin, and D. Terzopoulos, “Snakes: active contour models,” Int. J. Comput. Vis. 1(4), 321–331 (1988). 42. C. Li, C. Xu, C. Gui, and M. D. Fox, “Level set evolution without re-initialization: a new variational formulation,” in Proceedings of IEEE Conference on Computer Vision and Pattern Recognition (IEEE, 2005), pp. 430–436. 43. M. Ochs, J. R. Nyengaard, A. Jung, L. Knudsen, M. Voigt, T. Wahlers, J. Richter, and H. J. G. Gundersen, “The number of alveoli in the human lung,” Am. J. Respir. Crit. Care Med. 169(1), 120–124 (2004). 44. R. A. McLaughlin, L. Scolaro, P. Robbins, C. Saunders, S. L. Jacques, and D. D. Sampson, “Parametric imaging of cancer with optical coherence tomography,” J. Biomed. Opt. 15(4), 046029 (2010). 45. L. Scolaro, R. A. McLaughlin, B. R. Klyen, B. A. Wood, P. D. Robbins, C. M. Saunders, S. L. Jacques, and D. D. Sampson, “Parametric imaging of the local attenuation coefficient in human axillary lymph nodes assessed using optical coherence tomography,” Biomed. Opt. Express 3(2), 366–379 (2012). 46. L. Chin, X. Yang, R. A. McLaughlin, P. B. Noble, and D. D. Sampson, “En face parametric imaging of tissue birefringence using polarization-sensitive optical coherence tomography,” J. Biomed. Opt. 18(6), 066005 (2013). 47. J. P. Williamson, R. A. McLaughlin, W. J. Noffsinger, A. L. James, V. A. Baker, A. Curatolo, J. J. Armstrong, A. Regli, K. L. Shepherd, G. B. Marks, D. D. Sampson, D. R. Hillman, and P. R. Eastwood, “Elastic properties of the central airways in obstructive lung diseases measured using anatomical optical coherence tomography,” Am. J. Respir. Crit. Care Med. 183(5), 612–619 (2011). 48. H. Dayman, “Mechanics of airflow in health and in emphysema,” J. Clin. Invest. 30(11), 1175–1190 (1951). 49. S. Nava and F. Rubini, “Lung and chest wall mechanics in ventilated patients with end stage idiopathic pulmonary fibrosis,” Thorax 54(5), 390–395 (1999).
1. Introduction The primary function of the lungs is the diffusion of oxygen into and the release of carbon dioxide from the bloodstream. Air-filled sacs in the lung, called alveoli, facilitate this exchange of gases. Many lung diseases affect the structure of these alveolar air spaces, leading to impairment of lung function. For instance, chronic obstructive pulmonary disease (COPD) is characterized by the progressive destruction of alveolar walls; specifically, a degradation of the elastin. This reduces the recoil pressure of the lung parenchymal tissue [1], which includes the alveolar walls, blood vessels, and bronchi; causing the alveoli to become permanently dilated. In contrast, pulmonary fibrosis is characterized by the alteration of parenchymal collagen and infiltration of chronic inflammatory cells into the alveolar walls [2], causing stiffening and, thus, impeding alveolar dilation. COPD and pulmonary fibrosis combined account for over three quarters of all lung transplants [3]. Furthermore, COPD alone is estimated to be the third leading cause of death worldwide [4]. Characterization of disease-related changes in lung parenchyma, reflected in the size of lung structures, may allow assessment of disease progression [5]. Such sizing could assist in the assessment and development of new treatments for these pathologies. The current gold standard for evaluating the size of such structures is histological analysis of resected lung tissue. Stereological techniques [6] are used to compute measures that parameterize the lung tissue. One widely adopted measure is the linear intercept (chord) length [7], which is the mean-free distance in the air spaces. Note that this is not a measure of alveolar size, but
#193941 - $15.00 USD
Received 22 Jul 2013; revised 24 Sep 2013; accepted 26 Sep 2013; published 9 Oct 2013
(C) 2013 OSA 1 November 2013 | Vol. 4, No. 11 | DOI:10.1364/BOE.4.002383 | BIOMEDICAL OPTICS EXPRESS 2385
measures the mean distance across structural units comprising alveoli and alveolar ducts. However, the invasiveness of obtaining histological sections precludes longitudinal assessment and limits its use in live models. Furthermore, histological analysis cannot be used to assess the dynamic and viscoelastic properties of the lung, which are altered in lung disease [8]. Optical coherence tomography (OCT) is an optical imaging modality that has been demonstrated to image lung parenchyma. Imaging has been performed ex vivo and in vivo with endoscopic-based techniques on humans [9], as well as in vivo measurements via a thoracic window in rabbit and mouse models [10,11]. This latter technique involves removal of the skin and muscle over the rib cage, allowing access with an external OCT probe, and has been validated using confocal microscopy [12]. However, OCT has a limited penetration depth in tissue of 1-2 mm, and so has been restricted to imaging peripheral tissue [9–11]. Recent work has demonstrated the ability to image structures such as alveoli situated more deeply in the lung using an OCT needle probe [13,14]. This is achieved by encasing the miniaturized imaging optics within a hypodermic needle. Such techniques offer potential for greater coverage of the lung parenchyma, although the insertion of the needle will result in some tissue trauma. Whilst this method has been shown to be capable of visualizing large numbers of alveoli and alveolar ducts, the automated quantification of lung structures has not been addressed. The aim of the present work is to develop a method capable of automatically identifying parenchymal tissue visible in an OCT scan, and from this compute the OCT median chord length to provide an indicator of lung status. An automated method is required because manual sizing of large numbers of lung structures in a 3D OCT volume is infeasible. The accuracy of such automated quantification methods rely on the robustness of techniques to identify and delineate structures, a process referred to as image segmentation. Image segmentation techniques have been used previously in OCT to delineate tissue structures, typically in the eye [15–19]. Such techniques have the potential to identify lung parenchyma, enabling automatic characterization of these structures. Meissner et al. [10] demonstrated that the segmentation and quantification of alveoli is feasible, although the method proposed was not automated and was restricted to a small number of alveoli. In this paper, we present a fully automated segmentation and quantification algorithm capable of delineating large numbers of lung structures and computing estimates of average size. We demonstrate the algorithm on three-dimensional OCT data obtained with a needle probe from two healthy animal models: pig and rat, and validate the automated quantification against manual measurements. In addition, we present the first published 3D volume renderings of segmented lung parenchyma acquired with an OCT needle probe. 2. Methods Our proposed automated segmentation and quantification algorithm is organized into several stages. The first stage reduces speckle noise in the OCT scans to aid in segmentation. The second stage performs a preliminary delineation of structures through use of the Level Set method [20]. The third stage refines the segmentation obtained from the Level Set method, removing erroneously segmented structures. In a final stage of processing, the segmented OCT data is characterized using stereological techniques to compute the median size of the lung air spaces. In the following sections, we describe these stages in detail. 2.1 Noise reduction In OCT, speckle arises from the coherent summation of optical wavefields backscattered from sub-resolution scatterers in the sample arm of the interferometer [21,22]. In the context of this paper, it may be considered as a form of noise that compromises the ability of automated methods to identify structural features and boundaries in OCT images. Many hardware-based methods [22–26], and software-based methods [27,28] have been proposed to perform noise reduction whilst minimizing the effect on image sharpness.
#193941 - $15.00 USD
Received 22 Jul 2013; revised 24 Sep 2013; accepted 26 Sep 2013; published 9 Oct 2013
(C) 2013 OSA 1 November 2013 | Vol. 4, No. 11 | DOI:10.1364/BOE.4.002383 | BIOMEDICAL OPTICS EXPRESS 2386
In our algorithm, anisotropic diffusion [29] was implemented to reduce speckle whilst preserving the definition of edges. Anisotropic diffusion has been demonstrated previously in the removal of speckle from retinal OCT scans [30,31]. The present implementation [32] performs diffusion with a 2D kernel that adapts the degree of smoothing, based on the location and direction of intensity gradients in the image. This kernel promotes smoothing in directions with weak intensity gradients (areas of homogeneous tissue) and minimizes smoothing across strong intensity gradients (edges). The anisotropic diffusion algorithm is formulated as a diffusive process, as shown in Eq. (1): ∇I ( x, y , t ) 2 ∂ I ( x, y , t ) = ∇. exp − ∇I ( x , y , t ) , ∂t κ
(1)
where t is the algorithm iteration number; I(x,y,t) is the OCT radial B-scan intensity at coordinate (x,y) during iteration t; and κ is a positive constant chosen as a scaling factor. The negative exponential in this equation controls the rate of change of intensities at each iteration as a function of the strength of nearby edges features. For a given OCT volume, 2D anisotropic diffusion was performed within each B-scan. Subsequently, additional smoothing was performed between adjacent B-scans using weighted averaging. This operation averaged the intensity value at every individual pixel in the B-scan with the intensity of the corresponding pixels in both adjacent B-scans. This minimizes the effect of signal loss that is apparent on one B-scan but not present on adjacent scans. The loss of signal in these instances is the result of strong attenuation of the OCT signal, which can confound the delineation of lung structures. 2.2 Automatic delineation of lung structures The automatic delineation of lung structures from the OCT image was performed using a Level Set method first introduced by Osher and Sethian [20]. Using a standard approach for Level Sets, we represent the walls of the lung structures present in a single B-scan as a subset of an implicitly defined, higher dimensional function, referred to as a level set function. In our application, the walls visible in a B-scan correspond to a 2D plane (i.e., a level set) intersecting this 3D function. The 2D Level Set algorithm adopted in this work, described in detail in [33,34], does not require that we analytically define the implicit 3D function. Instead, it is initialized with a particular level set (corresponding to an initial segmentation of the B-scan), and then iteratively refined to minimize an error function. The error function is calculated from the difference between intensity values in the anisotropically smoothed OCT data and the implicitly defined 3D level set function; and also as a function of a measurement of smoothness of the shape of the segmentation outline. The error function F is shown in Eq. (2): F (φ , c, b) = ε (φ , c, b) + vL(φ ) + μ R (φ ),
(2)
where φ is the 3D level set function; c is a set of constant values which specifies a piecewise constant model of the OCT radial B-scan, approximating the intensity values within each contour with a constant value; and b is a smoothly varying field representing scanner-specific variation in B-scan intensities (i.e., a bias field). The function ε(φ,c,b) is used as the data term in the level set formulation, dependent upon differences between the OCT intensity values and the level set function. L(φ) and R(φ) are regularization terms, which serve to smooth the contour, and smooth the level set function, respectively. The constant terms ν and µ are weights to adjust the relative impact of the regularization terms. Full details of this Level Set implementation are given in [34]. A strength of the Level Set method is the robustness with which it will converge to an optimal segmentation from a wide range of initial level set configurations. In our implementation, the level set for each B-scan was initialized as a three by three grid pattern of square shapes, with dimensions approximately one tenth the length of the image. After #193941 - $15.00 USD
Received 22 Jul 2013; revised 24 Sep 2013; accepted 26 Sep 2013; published 9 Oct 2013
(C) 2013 OSA 1 November 2013 | Vol. 4, No. 11 | DOI:10.1364/BOE.4.002383 | BIOMEDICAL OPTICS EXPRESS 2387
convergence of the Level Set algorithm, a 3D level set function is computed. The delineation of lung structures is obtained from the zero level set of this function (i.e., the set of all points for which the level set function equals zero). Figure 1 shows representative results on a radial B-scan of rat lung. Figure 1(a) shows a radial B-scan acquired with the OCT needle probe. The zero level set is illustrated in Fig. 1(b), visualized as red contours super-imposed on the smoothed OCT B-scan. The interior of each red contour shape corresponds to an air space. Structures in each B-scan are segmented independently of adjacent B-scans. Performing the algorithm in 2D on a sequence of B-scans, rather than in a 3D implementation, both significantly reduces computation and better tolerates the high variation that can occur between B-scans. Note that the needle hole appears as a central black disc of no signal in Fig. 1(a). The extent of this ‘no signal’ region is extended in Fig. 1(b) as an artifact of the anisotropic diffusion performed to reduce the effects of speckle. However, we found this artifact to have minimal impact on the subsequent quantification described in Section 2.4. 2.3 Refinement of the segmented structures The results of the level set segmentation are refined through the use of morphological opening; by the application of edge-based consistency checks; and the inclusion of physiological a priori restrictions on the upper and lower limits of air space dimensions for a given animal model. A 2D morphological opening operation [35] is first performed on each B-scan of the segmented data, to separate adjacent structures that have been erroneously merged into a single structure. Morphological opening is a standard image processing technique, which consists of applying a binary erosion to the foreground pixels, followed by a binary dilation. This has the effect of separating adjacent structures whilst minimizing any artifactual reduction in the size estimate. Subsequently, an edge-based consistency check was applied to each separate segmented structure. Segmentation of our experimental results identified that the Level Set algorithm would sometimes delineate a structure even where the edge was poorly defined. This is particularly evident at locations far from the scanning probe, where attenuation of the light beam yields a poor signal-to-noise ratio and variations in the intensity gradients were often artifacts of noise rather than indicative of the boundary of underlying structures. Examples of this artifact are illustrated in Fig. 1(c), which shows a binary thresholding of the level set function, with air spaces denoted by white. Several additional erroneously segmented regions are labeled as noisy structures n1 – n9. These poorly delineated artifacts were filtered from the final segmentation by calculating the intensity gradient of the OCT radial B-scan at each delineating pixel, and computing the proportion of edge pixels above an empirically chosen edge strength. If more than half of the edge pixels for a segmented shape were below the threshold value, then the shape would be discarded as an artifact. The result of this was to only retain lung structures with strongly delineated OCT boundaries. Finally, a physiology-based consistency check was applied to specify acceptable bounds for the sizes of air spaces (spanning from alveolus to bronchiole). The alveolar size of many species has been well documented, and so we set our thresholds to range from approximately one quarter, to three times, the alveolar diameter estimated by Lum and Mitzner [36] for both pig and rodent (using mouse in lieu of rat). As there is a broad range of structure sizes, we intentionally set a wide range of acceptable sizes to include potential disease related changes whilst excluding gross segmentation errors. Structures with dimensions outside of this range were discarded. An example of a structure outside the bounds of acceptable size is illustrated in Fig. 1(c), labeled as incorrect structure n10. This artifactual feature was most likely caused by a slight tear in the tissue due to needle probe precession. Figure 1(d) shows the refined results after the filtering algorithm, illustrating the removal of the identified artifacts from the segmentation. After segmentation of each individual B-scan, the algorithm stacked the segmented B-scans together, forming a 3D volumetric segmentation of the lung parenchyma.
#193941 - $15.00 USD
Received 22 Jul 2013; revised 24 Sep 2013; accepted 26 Sep 2013; published 9 Oct 2013
(C) 2013 OSA 1 November 2013 | Vol. 4, No. 11 | DOI:10.1364/BOE.4.002383 | BIOMEDICAL OPTICS EXPRESS 2388
Fig. 1. (a) Radial B-scan image acquired from one rotation of the OCT needle probe. (b) Zero level set (red contour corresponding to the edges of structures) superimposed on the smoothed radial OCT B-scan. (c) Segmented image obtained by binarizing the level set function and performing morphological opening. This image contains both noisy structures (n1-n9) and an incorrect structure (n10). (d) Segmented image after post-processing with noisy and incorrect structures removed by edge-based consistency checking (refer to text).
2.4 Quantification To quantify the resulting segmentation, we adapted stereological methods [6] for determining the size of air spaces. These are standard methods typically applied to a set of 2D histological slices, and adapted here to quantify individual 2D radial scans of the segmented OCT data. We utilized chord length measurements [7] to quantify segmented OCT data. Beginning with a randomly selected point located on the segmented parenchyma, we defined an arbitrarily orientated line segment extending across the adjacent air space, and calculated the distance to the next intersection with parenchyma. Note that the resulting chord often crosses from one alveolus, through an alveolar duct, and terminates in an opposing alveolus. While individual chords may span a single alveolus or a bronchiole, we found the median chord length to be more indicative of the entire acinar air space (the space including the respiratory bronchiole and adjoining alveoli), as noted in [7]. Segmentation of the OCT scans allowed this measurement to be fully automated, and a population of chord lengths was computed over each data volume. The median and standard deviation of the chord lengths were then calculated, providing a characterization of the
#193941 - $15.00 USD
Received 22 Jul 2013; revised 24 Sep 2013; accepted 26 Sep 2013; published 9 Oct 2013
(C) 2013 OSA 1 November 2013 | Vol. 4, No. 11 | DOI:10.1364/BOE.4.002383 | BIOMEDICAL OPTICS EXPRESS 2389
geometry of the air spaces. As the distribution of chord lengths is not symmetric, we report the median instead of the mean. 3. Experiment 3.1 Tissue preparation All animal experiments conformed to institutional ethics approvals. Tissue was obtained from animals sacrificed as part of a different research project approved by The University of Western Australia Animal Ethics Committee, and utilized under institutional tissue sharing protocols. Two animals were used in this research: an adult rat, and an adult pig, and imaged according to the protocol described in [13]. The animals were euthanized on the day of the experiment, and dissection was performed immediately post-mortem. Both animals were exsanguinated prior to removal of lung lobes. Scans were acquired between 1 and 4 hours after dissection. The lungs were kept moist during this period with Krebs solution to prevent dehydration. Both the rat and pig lungs were allowed to deflate passively for 1 to 2 hours after excision. After tracheal cannulation, the lungs were inflated by a syringe filled with isotonic saline to improve image quality and minimize imaging artifacts that arise from the refractive index mismatch between the air-filled alveoli and the tissue of the parenchyma, as described in [13]. 3.2 Optical coherence tomography imaging OCT images were acquired using a spectral-domain OCT scanner with a central wavelength of 836 nm and a spectral bandwidth of 50 nm, giving a measured axial resolution of 9.3 µm in air. For the experiments performed here, the OCT scanner used an exposure time of 13.4 µs and was interfaced to an OCT needle probe, encased within a 30-gauge needle (outer diameter 310 µm). The design of this side-facing OCT needle probe has been detailed previously [37], but is briefly outlined here. A length of single-mode fiber was interfaced to the sample arm of the OCT scanner, and terminated with miniaturized focusing optics. The focusing optics comprise lengths of no-core fiber and graded-index (GRIN) fiber to produce a converging light beam and are terminated with an angle-polished, gold-coated length of nocore fiber to redirect the light beam at right angles to the fiber axis. The fully fused assembly is mounted in a 30-gauge hypodermic needle using adhesive. A small imaging window was etched in the side of the needle shaft, through which light is emitted and collected. The average full width at half maximum (FWHM) transverse resolution of the probe is 12 µm in tissue, over a range of 600 µm. The needle probe was mounted onto a stepper motor and counter-rotated over a range of 360 degrees at a frequency of 1 Hz to acquire 2D radial Bscans. Consecutive radial B-scans were acquired by translating the needle in steps of 5 µm, over a pullback distance of 2 mm. This produced a 3D cylindrical volumetric data set with an imaging radius of approximately 600 µm. 3.3 Data processing and visualization The segmentation algorithm was implemented using MATLAB (Mathworks, Natick, MA, USA) run on a PC running Windows 7 64-bit, with 16 GB of RAM and a quad-core Intel i5 processor. We developed custom software to visualize the 3D segmented data. This software was written in C++ and used the Visualization Toolkit (VTK) (Kitware, NY, USA). A Marching Cubes algorithm was used to create a 3D surface from the segmented data. This rendered surface and the OCT data were then graphically overlaid. The opacity of the OCT data was modified to allow simultaneous visualization of both the segmented data and the underlying OCT data volume, allowing a visual assessment of accuracy of the segmentation algorithm. A sample of 10,000 chord lengths was calculated over each data volume, and the median and standard deviation of these measurements were computed. To validate these measurements, we additionally performed 1,000 manual measurements, where a trained
#193941 - $15.00 USD
Received 22 Jul 2013; revised 24 Sep 2013; accepted 26 Sep 2013; published 9 Oct 2013
(C) 2013 OSA 1 November 2013 | Vol. 4, No. 11 | DOI:10.1364/BOE.4.002383 | BIOMEDICAL OPTICS EXPRESS 2390
human operator selected the beginning and end of each chord, and the median and standard deviation of these were compared against the automated measurements. 4. Results Figures 2 and 3 show representative results from the segmentation algorithm, demonstrated on the rat and pig, respectively. Figure 2(a) shows a representative radial B-scan of the rat lung. As in Fig. 1(a), the needle hole appears as a black void at the center of the image, surrounded by alveoli and bronchioles. Figure 2(b) shows the corresponding segmented radial scan of Fig. 2(a), graphically representing the segmented structures as white regions. This image demonstrates a range of segmented structures, including individual alveoli and acini (i.e., structural units comprising the respiratory bronchioles, alveolar ducts and alveolar sacs). A 3D rendering of the segmentation results is shown in Fig. 2(c) (Media 1), overlaid on the original OCT data.
Fig. 2. Segmentation results obtained from the OCT image of rat lung. (a) Radial OCT image acquired from one rotation of the needle probe. (b) Segmented radial image corresponding to the OCT image in (a). (c) Clipped portion of the 3D segmented data overlaid on the partially transparent OCT data (Media 1). Green structures correspond to structures segmented by the algorithm. Faded green structures indicate the segmented structures lying underneath the OCT data.
Figure 3(a) shows a representative radial scan acquired from the pig lung, and the corresponding segmentation in Fig. 3(b). The 3D rendered visualization of this data set is shown in Fig. 3(c) (Media 2). These visualizations show, to our knowledge, the first 3D volume renderings of lung structures acquired with an OCT needle probe. In these
#193941 - $15.00 USD
Received 22 Jul 2013; revised 24 Sep 2013; accepted 26 Sep 2013; published 9 Oct 2013
(C) 2013 OSA 1 November 2013 | Vol. 4, No. 11 | DOI:10.1364/BOE.4.002383 | BIOMEDICAL OPTICS EXPRESS 2391
visualizations, the structures segmented by the algorithm are represented as green surfaces. Only a clipped portion of the segmented data is shown to improve clarity of the visualization. In addition to the segmentation of structures, Figs. 2 and 3 also illustrate the removal of some structures introduced as an artifact of OCT needle imaging. Some tearing of tissue by the needle probe is apparent above the central needle hole in Fig. 2(a) and to the right of the hole in Fig. 3(a). These regions, which do not correspond to the natural air spaces of the lung, were filtered out by the algorithm and, therefore, do not appear in the segmented data.
Fig. 3. Segmentation results obtained from the OCT image of pig lung. (a) Radial OCT image. (b) Segmented radial image corresponding to the OCT image in (a). (c) A 3D visualization of the segmented data overlaid on the partially transparent OCT data (Media 2). Green structures correspond to air spaces segmented by the algorithm.
Table 1 illustrates the median and the standard deviation of air space sizes, measured both automatically (10,000 measurements) from the segmented data and manually (1,000 measurements) from the raw OCT data, for both the pig and rat data volumes. The automated median chord length calculated from the segmented data was 7% lower than the manual estimate for the pig lung, and 6% lower for the rat lung. These values lie within the standard deviation of the manual measurements. The standard deviation of the automatic chord length measurements differed by less than 1 µm from the standard deviation obtained from the manual measurements for both the pig and rat lung.
#193941 - $15.00 USD
Received 22 Jul 2013; revised 24 Sep 2013; accepted 26 Sep 2013; published 9 Oct 2013
(C) 2013 OSA 1 November 2013 | Vol. 4, No. 11 | DOI:10.1364/BOE.4.002383 | BIOMEDICAL OPTICS EXPRESS 2392
Table 1. Comparison of the median and standard deviation about the median of air space sizes obtained from automated and manual methods for both the pig and rat data. Median and standard deviation of air space size quantified in pig lung [median ± standard deviation] (µm) Automatic Manual % Quantification Quantification Difference 75 ± 5 80 ± 5 7
Median and standard deviation of air space size quantified in rat lung [median ± standard deviation] (µm) Automatic Manual % Quantification Quantification Difference 62 ± 5 65 ± 5 6
5. Discussion The capacity to quantify the dimensions of lung structures, either singly or multiply at different pressure levels, may enable new protocols for assessing lung diseases that affect the structure of the lung. The results presented in this paper highlight the potential for an algorithm to segment and quantify lung structures in an automated way, demonstrated in pig and rat models. The ability to segment structures with a wide range of shapes is an advantage of the Level Set method, and its representation of this complex morphology through the use of a level set function. The accompanying animations demonstrate that the automatic segmentation method achieved a useful segmentation of the OCT data. These results have been quantitatively validated against manual stereological estimates of air space size. The quantification of median size obtained from the segmentation was ~7% lower than the manual measurements for both pig and rat data. This difference can be attributed to regions of low signal-to-noise ratio in the OCT image, which were filtered out to improve robustness of the automated segmentation, but utilized in the manual measurements. However, this error translates to a disparity in chord length of ≤ 5 µm between the segmented data and the OCT data measurements, which is less than one transverse resolution element of the OCT system (12 µm). In addition, this difference is less than the observed physiological changes in the size of airway structures observed under several disease conditions. For example, COPD has been found to double the mean chord length [7] and increase the weighted mean alveolar diameter by approximately 80% [5] in comparison to healthy lung. This suggests that the segmentation is a sufficiently accurate characterization of the OCT data to distinguish healthy lung from various disease processes, and highlights the potential for diagnostically relevant quantifications to be made from this segmentation. We note that some pathological changes could present additional challenges to automated segmentation of OCT images, such as the thinning of the alveoli walls that is observed in COPD. However, extremely high sensitivity OCT imaging probes, such as that demonstrated in [38], should offer an improved ability to detect such features. The segmentation algorithm was observed to sometimes produce erroneous segmentations at locations with a low OCT signal-to-noise ratio (SNR). This is illustrated in Fig. 4, showing an erroneous segmentation of part of a bronchiole wall. Figure 4(a) shows the radial B-scan after anisotropic diffusion to reduce the effects of speckle. The zero level set, delineating the parenchymal structures, is superimposed on the image as a red contour. The bronchiole is highlighted by a yellow circle. The OCT signal from the far wall of the bronchiole has suffered significant attenuation, making the outline hard to distinguish, and resulting in the segmentation highlighted in Fig. 4(b). A small number of these structures are not removed during the post-processing described in Section 2.3, and remain in the final segmentation. Conversely, in a small number of cases, low OCT signal was found to obscure parts of the parenchyma between adjacent air spaces. This resulted in a segmentation that aggregated these air spaces. An example is shown in Fig. 4(c), where the parenchymal wall between an acinus and an area of void created by the needle probe is obscured, resulting in a segmentation which combines these into a single, connected air space. Subsequent filtering of abnormally large air spaces, typically indicative of such segmentation errors, resulted in removal of the acinus, seen in Fig. 4(d). Both of these classes of error could be reduced with the implementation of a 3D Level Set algorithm, which would combine information from
#193941 - $15.00 USD
Received 22 Jul 2013; revised 24 Sep 2013; accepted 26 Sep 2013; published 9 Oct 2013
(C) 2013 OSA 1 November 2013 | Vol. 4, No. 11 | DOI:10.1364/BOE.4.002383 | BIOMEDICAL OPTICS EXPRESS 2393
adjacent B-scans. Extending the Level Set algorithm to 3D is feasible and has previously been demonstrated [39,40], but also results in a significantly increased computational load for the algorithm.
Fig. 4. Two examples illustrating errors in segmentation of lung structures due to low OCT SNR (a) Segmented radial image after anisotropic diffusion. Red contour shows zero level set. Yellow circle highlights bronchiole with far wall poorly delineated. (b) Resulting segmentation omitting far wall of bronchiole. (c) Yellow circle highlights an acinus apparently joined to an area of void caused by the needle probe. The white arrow indicates the location at which the lung parenchyma is obscured by poor OCT SNR (d) Resulting segmentation with acinus omitted.
The Level Set method was chosen for our application as it has several advantages over alternative image segmentation techniques, such as parametric active contours [41]. Firstly, the level set function representation of structure means that complex morphology, such as the intricate network of alveolar sacs in the lung, can be handled in a robust manner [34]. This level set representation can also handle changes to this complex topology, such as the bifurcation or merging of airways, which is common in the lung [34,42]. Furthermore, the Level Set method does not require parameterization, which simplifies numerical computation. An advantage of this particular implementation of level sets is that it exhibits robustness of segmentation in the presence of some inhomogeneity of signal intensity. Using a technique originally designed for inhomogeneity in MRI images, the algorithm addresses this issue by generating a bias field and updating its estimate of the bias field with each iteration, along with the level set function, until it resembles the underlying intensity inhomogeneity in the image data.
#193941 - $15.00 USD
Received 22 Jul 2013; revised 24 Sep 2013; accepted 26 Sep 2013; published 9 Oct 2013
(C) 2013 OSA 1 November 2013 | Vol. 4, No. 11 | DOI:10.1364/BOE.4.002383 | BIOMEDICAL OPTICS EXPRESS 2394
Additional stereological measures may be computed from segmented data, such as the volume fraction of parenchymal tissue [6], total alveolar surface area [6] and total alveoli number [43]. The choice of measure depends both on the pathological variations caused by the disease, and the physiological question being posed. In addition, parametric imaging techniques provide an alternate approach to quantification of the underlying structure of the tissue. Such techniques have been demonstrated to quantify the attenuation coefficient of tissue [44,45], or automatically quantify the degree of birefringence [46]. These approaches rely on segmentation of tissue types, as both the attenuation coefficient and birefringence will vary between different tissues. We suggest that automated segmentation techniques, such as that proposed in this paper, represent one approach to improving the quantification accuracy of these parametric imaging techniques. The method reported here has potential in the assessment of lung compliance, by computing measurements on lung tissue acquired at different levels of expansion. Dynamic OCT scans of lungs undergoing simulated respiration have previously been demonstrated in [13], producing a sequence of lung scans at different stages of alveolar expansion and recruitment. Alternatively, Williamson et al. [47] demonstrated the acquisition of a series of OCT scans at predefined lung pressures. Whilst this latter study quantified changes in bronchial dimensions rather than alveolar spaces, it demonstrates the practicality of acquiring OCT scans at predefined and well-controlled pressures to characterize compliance of lung parenchyma in a range of diseases. Quantification of lung structure dimensions at each pressure, and changes across pressures, could provide an indicator of changes in tissue compliance due to tissue remodeling, which is significant in diseases such as emphysema and pulmonary fibrosis [48,49]. 6. Conclusions This paper has presented a fully automated segmentation and quantification algorithm capable of delineating large numbers of lung structures from an OCT volume. The algorithm was demonstrated on two animals: pig and rat. In both animals, the majority of lung structures visible in the OCT data were accurately segmented. Quantification of air space size was validated against manual measurements, and the automated estimates of median chord length were within 5 µm of manual measurements. This paper also presents the first published 3D volume renderings of lung structures acquired with an OCT needle probe. Our segmentation method may enable creation of new protocols for characterizing lung disease through objective and automated measurement of the air spaces, with the potential to assist in the assessment of new treatments for these diseases. Acknowledgments Robert A. McLaughlin is supported by a fellowship from Cancer Council Western Australia and Brendan F. Kennedy is supported by a Raine Priming Grant from the Raine Medical Research Foundation. This research has been supported in part by funding from the National Breast Cancer Foundation, Australia. We wish to acknowledge A/Prof. Peter Noble and Prof. J. Pillow, Centre for Neonatal Research and Education, UWA, for provision of the lung tissue used in this research. We also wish to acknowledge our associate Bryden Quirk, for his assistance in construction of the OCT needle probe used in these experiments.
#193941 - $15.00 USD
Received 22 Jul 2013; revised 24 Sep 2013; accepted 26 Sep 2013; published 9 Oct 2013
(C) 2013 OSA 1 November 2013 | Vol. 4, No. 11 | DOI:10.1364/BOE.4.002383 | BIOMEDICAL OPTICS EXPRESS 2395
Optical + Biomedical Engineering Laboratory, School of Electrical, Electronic and Computer Engineering, The University of Western Australia, 35 Stirling Highway, Crawley, Western Australia 6009, Australia Centre for Microscopy, Characterisation and Analysis, The University of Western Australia, 35 Stirling Highway, Crawley, Western Australia, 6009, Australia * [email protected]
Abstract: Characterization of the size of lung structures can aid in the assessment of a range of respiratory diseases. In this paper, we present a fully automated segmentation and quantification algorithm for the delineation of large numbers of lung structures in optical coherence tomography images, and the characterization of their size using the stereological measure of median chord length. We demonstrate this algorithm on scans acquired with OCT needle probes in fresh, ex vivo tissues from two healthy animal models: pig and rat. Automatically computed estimates of lung structure size were validated against manual measures. In addition, we present 3D visualizations of the lung structures using the segmentation calculated for each data set. This method has the potential to provide an in vivo indicator of structural remodeling caused by a range of respiratory diseases, including chronic obstructive pulmonary disease and pulmonary fibrosis. ©2013 Optical Society of America OCIS codes: (110.4500) Optical coherence tomography; (100.0100) Image processing; (100.2960) Image analysis.
References and links 1. 2. 3. 4.
G. Turato, R. Zuin, and M. Saetta, “Pathogenesis and pathology of COPD,” Respiration 68(2), 117–128 (2001). T. C. Kravis, A. Ahmed, T. E. Brown, J. D. Fulmer, and R. G. Crystal, “Pathogenic mechanisms in pulmonary fibrosis: collagen-induced migration inhibition factor production and cytotoxicity mediated by lymphocytes,” J. Clin. Invest. 58(5), 1223–1232 (1976). E. P. Trulock, L. B. Edwards, D. O. Taylor, M. M. Boucek, B. M. Keck, and M. I. Hertz, “Registry of the International Society for Heart and Lung Transplantation: twenty-second official adult lung and heart-lung transplant report--2005,” J. Heart Lung Transplant. 24(8), 956–967 (2005). R. Lozano, M. Naghavi, K. Foreman, S. Lim, K. Shibuya, V. Aboyans, J. Abraham, T. Adair, R. Aggarwal, S. Y. Ahn, M. Alvarado, H. R. Anderson, L. M. Anderson, K. G. Andrews, C. Atkinson, L. M. Baddour, S. BarkerCollo, D. H. Bartels, M. L. Bell, E. J. Benjamin, D. Bennett, K. Bhalla, B. Bikbov, A. Bin Abdulhak, G. Birbeck, F. Blyth, I. Bolliger, S. Boufous, C. Bucello, M. Burch, P. Burney, J. Carapetis, H. Chen, D. Chou, S. S. Chugh, L. E. Coffeng, S. D. Colan, S. Colquhoun, K. E. Colson, J. Condon, M. D. Connor, L. T. Cooper, M. Corriere, M. Cortinovis, K. C. de Vaccaro, W. Couser, B. C. Cowie, M. H. Criqui, M. Cross, K. C. Dabhadkar, N. Dahodwala, D. De Leo, L. Degenhardt, A. Delossantos, J. Denenberg, D. C. Des Jarlais, S. D. Dharmaratne, E. R. Dorsey, T. Driscoll, H. Duber, B. Ebel, P. J. Erwin, P. Espindola, M. Ezzati, V. Feigin, A. D. Flaxman, M. H. Forouzanfar, F. G. Fowkes, R. Franklin, M. Fransen, M. K. Freeman, S. E. Gabriel, E. Gakidou, F. Gaspari, R. F. Gillum, D. Gonzalez-Medina, Y. A. Halasa, D. Haring, J. E. Harrison, R. Havmoeller, R. J. Hay, B. Hoen, P. J. Hotez, D. Hoy, K. H. Jacobsen, S. L. James, R. Jasrasaria, S. Jayaraman, N. Johns, G. Karthikeyan, N. Kassebaum, A. Keren, J. P. Khoo, L. M. Knowlton, O. Kobusingye, A. Koranteng, R. Krishnamurthi, M. Lipnick, S. E. Lipshultz, S. L. Ohno, J. Mabweijano, M. F. MacIntyre, L. Mallinger, L. March, G. B. Marks, R. Marks, A. Matsumori, R. Matzopoulos, B. M. Mayosi, J. H. McAnulty, M. M. McDermott, J. McGrath, G. A. Mensah, T. R. Merriman, C. Michaud, M. Miller, T. R. Miller, C. Mock, A. O. Mocumbi, A. A. Mokdad, A. Moran, K. Mulholland, M. N. Nair, L. Naldi, K. M. Narayan, K. Nasseri, P. Norman, M. O’Donnell, S. B. Omer, K. Ortblad, R. Osborne, D. Ozgediz, B. Pahari, J. D. Pandian, A. P. Rivero, R. P. Padilla, F. Perez-Ruiz, N. Perico, D. Phillips, K. Pierce, C. A. Pope 3rd, E. Porrini, F. Pourmalek, M. Raju, D. Ranganathan, J. T. Rehm, D. B. Rein, G. Remuzzi, F. P. Rivara, T. Roberts, F. R. De León, L. C. Rosenfeld, L. Rushton, R. L. Sacco, J. A. Salomon, U. Sampson, E. Sanman, D. C. Schwebel, M. Segui-Gomez, D. S. Shepard, D. Singh, J. Singleton, K. Sliwa, E. Smith, A. Steer, J. A. Taylor, B. Thomas, I. M. Tleyjeh, J. A. Towbin, T. Truelsen, E. A. Undurraga,
#193941 - $15.00 USD
Received 22 Jul 2013; revised 24 Sep 2013; accepted 26 Sep 2013; published 9 Oct 2013
(C) 2013 OSA 1 November 2013 | Vol. 4, No. 11 | DOI:10.1364/BOE.4.002383 | BIOMEDICAL OPTICS EXPRESS 2383
5. 6. 7. 8. 9. 10. 11. 12. 13. 14. 15. 16. 17. 18. 19. 20. 21. 22. 23. 24. 25. 26. 27. 28. 29. 30.
N. Venketasubramanian, L. Vijayakumar, T. Vos, G. R. Wagner, M. Wang, W. Wang, K. Watt, M. A. Weinstock, R. Weintraub, J. D. Wilkinson, A. D. Woolf, S. Wulf, P. H. Yeh, P. Yip, A. Zabetian, Z. J. Zheng, A. D. Lopez, C. J. Murray, M. A. AlMazroa, and Z. A. Memish, “Global and regional mortality from 235 causes of death for 20 age groups in 1990 and 2010: a systematic analysis for the global burden of disease study 2010,” Lancet 380(9859), 2095–2128 (2012). H. Parameswaran, A. Majumdar, S. Ito, A. M. Alencar, and B. Suki, “Quantitative characterization of airspace enlargement in emphysema,” J. Appl. Physiol. 100(1), 186–193 (2006). E. R. Weibel, C. C. W. Hsia, and M. Ochs, “How much is there really? Why stereology is essential in lung morphometry,” J. Appl. Physiol. 102(1), 459–467 (2006). L. Knudsen, E. R. Weibel, H. J. G. Gundersen, F. V. Weinstein, and M. Ochs, “Assessment of air space size characteristics by intercept (chord) measurement: an accurate and efficient stereological approach,” J. Appl. Physiol. 108(2), 412–421 (2010). D. S. Faffe and W. A. Zin, “Lung parenchymal mechanics in health and disease,” Physiol. Rev. 89(3), 759–775 (2009). L. P. Hariri, M. B. Applegate, M. Mino-Kenudson, E. J. Mark, B. D. Medoff, A. D. Luster, B. E. Bouma, G. J. Tearney, and M. J. Suter, “Volumetric optical frequency domain imaging of pulmonary pathology with precise correlation to histopathology,” Chest 143(1), 64–74 (2013). S. Meissner, L. Knels, C. Schnabel, T. Koch, and E. Koch, “Three-dimensional Fourier domain optical coherence tomography in vivo imaging of alveolar tissue in the intact thorax using the parietal pleura as a window,” J. Biomed. Opt. 15(1), 016030 (2010). M. Mertens, A. Tabuchi, S. Meissner, A. Krueger, K. Schirrmann, U. Kertzscher, A. R. Pries, A. S. Slutsky, E. Koch, and W. M. Kuebler, “Alveolar dynamics in acute lung injury: heterogeneous distension rather than cyclic opening and collapse,” Crit. Care Med. 37(9), 2604–2611 (2009). J. Bickenbach, R. Dembinski, M. Czaplik, S. Meissner, A. Tabuchi, M. Mertens, L. Knels, W. Schroeder, P. Pelosi, E. Koch, W. M. Kuebler, R. Rossaint, and R. Kuhlen, “Comparison of two in vivo microscopy techniques to visualize alveolar mechanics,” J. Clin. Monit. Comput. 23(5), 323–332 (2009). R. A. McLaughlin, X. Yang, B. C. Quirk, D. Lorenser, R. W. Kirk, P. B. Noble, and D. D. Sampson, “Static and dynamic imaging of alveoli using optical coherence tomography needle probes,” J. Appl. Physiol. 113(6), 967– 974 (2012). B. C. Quirk, R. A. McLaughlin, A. Curatolo, R. W. Kirk, P. B. Noble, and D. D. Sampson, “In situ imaging of lung alveoli with an optical coherence tomography needle probe,” J. Biomed. Opt. 16(3), 036009 (2011). D. Cabrera Fernández, H. M. Salinas, and C. A. Puliafito, “Automated detection of retinal layer structures on optical coherence tomography images,” Opt. Express 13(25), 10200–10216 (2005). M. K. Garvin, M. D. Abràmoff, X. Wu, S. R. Russell, T. L. Burns, and M. Sonka, “Automated 3-D intraretinal layer segmentation of macular spectral-domain optical coherence tomography images,” IEEE Trans. Med. Imaging 28(9), 1436–1447 (2009). A. Mishra, A. Wong, K. Bizheva, and D. A. Clausi, “Intra-retinal layer segmentation in optical coherence tomography images,” Opt. Express 17(26), 23719–23728 (2009). A. Yazdanpanah, G. Hamarneh, B. Smith, and M. Sarunic, “Intra-retinal layer segmentation in optical coherence tomography using an active contour approach,” in MICCAI 2009, G-Z Yang, D. J. Hawkes, D. Rueckert, A. Noble, C. J. Taylor, ed. (Springer, Heidelberg, 2009), pp. 649–656. R. J. Zawadzki, A. R. Fuller, D. F. Wiley, B. Hamann, S. S. Choi, and J. S. Werner, “Adaptation of a support vector machine algorithm for segmentation and visualization of retinal structures in volumetric optical coherence tomography data sets,” J. Biomed. Opt. 12(4), 041206 (2007). S. Osher and J. A. Sethian, “Fronts propagating with curvature-dependent speed: algorithms based on HamiltonJacobi formulations,” J. Comput. Phys. 79(1), 12–49 (1988). A. Curatolo, B. F. Kennedy, D. D. Sampson, and T. R. Hillman, “Speckle in optical coherence tomography,” in Advanced Biophotonics: Tissue Optical Sectioning, V. V. Tuchin and R. K. Wang, eds. (Taylor & Francis, 2013), pp. 183–249. J. M. Schmitt, S. H. Xiang, and K. M. Yung, “Speckle in optical coherence tomography,” J. Biomed. Opt. 4(1), 95–105 (1999). A. E. Desjardins, B. J. Vakoc, W. Y. Oh, S. M. R. Motaghiannezam, G. J. Tearney, and B. E. Bouma, “Angleresolved optical coherence tomography with sequential angular selectivity for speckle reduction,” Opt. Express 15(10), 6200–6209 (2007). T. R. Hillman, S. G. Adie, V. Seemann, J. J. Armstrong, S. L. Jacques, and D. D. Sampson, “Correlation of static speckle with sample properties in optical coherence tomography,” Opt. Lett. 31(2), 190–192 (2006). B. F. Kennedy, T. R. Hillman, A. Curatolo, and D. D. Sampson, “Speckle reduction in optical coherence tomography by strain compounding,” Opt. Lett. 35(14), 2445–2447 (2010). B. F. Kennedy, A. Curatolo, T. R. Hillman, C. M. Saunders, and D. D. Sampson, “Speckle reduction in optical coherence tomography images using tissue viscoelasticity,” J. Biomed. Opt. 16(2), 020506 (2011). D. C. Adler, T. H. Ko, and J. G. Fujimoto, “Speckle reduction in optical coherence tomography images by use of a spatially adaptive wavelet filter,” Opt. Lett. 29(24), 2878–2880 (2004). J. Rogowska and M. E. Brezinski, “Image processing techniques for noise removal, enhancement and segmentation of cartilage OCT images,” Phys. Med. Biol. 47(4), 641–655 (2002). P. Perona and J. Malik, “Scale-space and edge detection using anisotropic diffusion,” IEEE Trans. Pattern Anal. 12(7), 629–639 (1990). D. C. Fernandez, “Delineating fluid-filled region boundaries in optical coherence tomography images of the retina,” IEEE Trans. Med. Imaging 24(8), 929–945 (2005).
#193941 - $15.00 USD
Received 22 Jul 2013; revised 24 Sep 2013; accepted 26 Sep 2013; published 9 Oct 2013
(C) 2013 OSA 1 November 2013 | Vol. 4, No. 11 | DOI:10.1364/BOE.4.002383 | BIOMEDICAL OPTICS EXPRESS 2384
31. P. Puvanathasan and K. Bizheva, “Interval type-II fuzzy anisotropic diffusion algorithm for speckle noise reduction in optical coherence tomography images,” Opt. Express 17(2), 733–746 (2009). 32. D. Lopes, “Anisotropic diffusion (Perona & Malik),” http://www.mathworks.com.au/matlabcentral/fileexchange/14995-anisotropic-diffusion-perona-malik. 33. C. Li, “Level set for image segmentation,” http://www.mathworks.com.au/matlabcentral/fileexchange/12711level-set-for-image-segmentation. 34. C. Li, R. Huang, Z. Ding, J. C. Gatenby, D. N. Metaxas, and J. C. Gore, “A level set method for image segmentation in the presence of intensity inhomogeneities with application to MRI,” IEEE Trans. Image Process. 20(7), 2007–2016 (2011). 35. J. Serra, Image Analysis and Mathematical Morphology (Academic Press, 1982). 36. H. Lum and W. Mitzner, “A species comparison of alveolar size and surface forces,” J. Appl. Physiol. 62(5), 1865–1871 (1987). 37. D. Lorenser, X. Yang, R. W. Kirk, B. C. Quirk, R. A. McLaughlin, and D. D. Sampson, “Ultrathin side-viewing needle probe for optical coherence tomography,” Opt. Lett. 36(19), 3894–3896 (2011). 38. L. Scolaro, D. Lorenser, R. A. McLaughlin, B. C. Quirk, R. W. Kirk, and D. D. Sampson, “High-sensitivity anastigmatic imaging needle for optical coherence tomography,” Opt. Lett. 37(24), 5247–5249 (2012). 39. A. E. Lefohn, J. E. Cates, and R. T. Whitaker, “Interactive, GPU-based level sets for 3D brain tumor segmentation,” in MICCAI 2003, R. E. Ellis, T. M. Peters, ed. (Springer, Heidelberg, 2003), pp. 564–572. 40. R. T. Whitaker, “A level-set approach to 3D reconstruction from range data,” Int. J. Comput. Vis. 29(3), 203– 231 (1998). 41. M. Kass, A. Witkin, and D. Terzopoulos, “Snakes: active contour models,” Int. J. Comput. Vis. 1(4), 321–331 (1988). 42. C. Li, C. Xu, C. Gui, and M. D. Fox, “Level set evolution without re-initialization: a new variational formulation,” in Proceedings of IEEE Conference on Computer Vision and Pattern Recognition (IEEE, 2005), pp. 430–436. 43. M. Ochs, J. R. Nyengaard, A. Jung, L. Knudsen, M. Voigt, T. Wahlers, J. Richter, and H. J. G. Gundersen, “The number of alveoli in the human lung,” Am. J. Respir. Crit. Care Med. 169(1), 120–124 (2004). 44. R. A. McLaughlin, L. Scolaro, P. Robbins, C. Saunders, S. L. Jacques, and D. D. Sampson, “Parametric imaging of cancer with optical coherence tomography,” J. Biomed. Opt. 15(4), 046029 (2010). 45. L. Scolaro, R. A. McLaughlin, B. R. Klyen, B. A. Wood, P. D. Robbins, C. M. Saunders, S. L. Jacques, and D. D. Sampson, “Parametric imaging of the local attenuation coefficient in human axillary lymph nodes assessed using optical coherence tomography,” Biomed. Opt. Express 3(2), 366–379 (2012). 46. L. Chin, X. Yang, R. A. McLaughlin, P. B. Noble, and D. D. Sampson, “En face parametric imaging of tissue birefringence using polarization-sensitive optical coherence tomography,” J. Biomed. Opt. 18(6), 066005 (2013). 47. J. P. Williamson, R. A. McLaughlin, W. J. Noffsinger, A. L. James, V. A. Baker, A. Curatolo, J. J. Armstrong, A. Regli, K. L. Shepherd, G. B. Marks, D. D. Sampson, D. R. Hillman, and P. R. Eastwood, “Elastic properties of the central airways in obstructive lung diseases measured using anatomical optical coherence tomography,” Am. J. Respir. Crit. Care Med. 183(5), 612–619 (2011). 48. H. Dayman, “Mechanics of airflow in health and in emphysema,” J. Clin. Invest. 30(11), 1175–1190 (1951). 49. S. Nava and F. Rubini, “Lung and chest wall mechanics in ventilated patients with end stage idiopathic pulmonary fibrosis,” Thorax 54(5), 390–395 (1999).
1. Introduction The primary function of the lungs is the diffusion of oxygen into and the release of carbon dioxide from the bloodstream. Air-filled sacs in the lung, called alveoli, facilitate this exchange of gases. Many lung diseases affect the structure of these alveolar air spaces, leading to impairment of lung function. For instance, chronic obstructive pulmonary disease (COPD) is characterized by the progressive destruction of alveolar walls; specifically, a degradation of the elastin. This reduces the recoil pressure of the lung parenchymal tissue [1], which includes the alveolar walls, blood vessels, and bronchi; causing the alveoli to become permanently dilated. In contrast, pulmonary fibrosis is characterized by the alteration of parenchymal collagen and infiltration of chronic inflammatory cells into the alveolar walls [2], causing stiffening and, thus, impeding alveolar dilation. COPD and pulmonary fibrosis combined account for over three quarters of all lung transplants [3]. Furthermore, COPD alone is estimated to be the third leading cause of death worldwide [4]. Characterization of disease-related changes in lung parenchyma, reflected in the size of lung structures, may allow assessment of disease progression [5]. Such sizing could assist in the assessment and development of new treatments for these pathologies. The current gold standard for evaluating the size of such structures is histological analysis of resected lung tissue. Stereological techniques [6] are used to compute measures that parameterize the lung tissue. One widely adopted measure is the linear intercept (chord) length [7], which is the mean-free distance in the air spaces. Note that this is not a measure of alveolar size, but
#193941 - $15.00 USD
Received 22 Jul 2013; revised 24 Sep 2013; accepted 26 Sep 2013; published 9 Oct 2013
(C) 2013 OSA 1 November 2013 | Vol. 4, No. 11 | DOI:10.1364/BOE.4.002383 | BIOMEDICAL OPTICS EXPRESS 2385
measures the mean distance across structural units comprising alveoli and alveolar ducts. However, the invasiveness of obtaining histological sections precludes longitudinal assessment and limits its use in live models. Furthermore, histological analysis cannot be used to assess the dynamic and viscoelastic properties of the lung, which are altered in lung disease [8]. Optical coherence tomography (OCT) is an optical imaging modality that has been demonstrated to image lung parenchyma. Imaging has been performed ex vivo and in vivo with endoscopic-based techniques on humans [9], as well as in vivo measurements via a thoracic window in rabbit and mouse models [10,11]. This latter technique involves removal of the skin and muscle over the rib cage, allowing access with an external OCT probe, and has been validated using confocal microscopy [12]. However, OCT has a limited penetration depth in tissue of 1-2 mm, and so has been restricted to imaging peripheral tissue [9–11]. Recent work has demonstrated the ability to image structures such as alveoli situated more deeply in the lung using an OCT needle probe [13,14]. This is achieved by encasing the miniaturized imaging optics within a hypodermic needle. Such techniques offer potential for greater coverage of the lung parenchyma, although the insertion of the needle will result in some tissue trauma. Whilst this method has been shown to be capable of visualizing large numbers of alveoli and alveolar ducts, the automated quantification of lung structures has not been addressed. The aim of the present work is to develop a method capable of automatically identifying parenchymal tissue visible in an OCT scan, and from this compute the OCT median chord length to provide an indicator of lung status. An automated method is required because manual sizing of large numbers of lung structures in a 3D OCT volume is infeasible. The accuracy of such automated quantification methods rely on the robustness of techniques to identify and delineate structures, a process referred to as image segmentation. Image segmentation techniques have been used previously in OCT to delineate tissue structures, typically in the eye [15–19]. Such techniques have the potential to identify lung parenchyma, enabling automatic characterization of these structures. Meissner et al. [10] demonstrated that the segmentation and quantification of alveoli is feasible, although the method proposed was not automated and was restricted to a small number of alveoli. In this paper, we present a fully automated segmentation and quantification algorithm capable of delineating large numbers of lung structures and computing estimates of average size. We demonstrate the algorithm on three-dimensional OCT data obtained with a needle probe from two healthy animal models: pig and rat, and validate the automated quantification against manual measurements. In addition, we present the first published 3D volume renderings of segmented lung parenchyma acquired with an OCT needle probe. 2. Methods Our proposed automated segmentation and quantification algorithm is organized into several stages. The first stage reduces speckle noise in the OCT scans to aid in segmentation. The second stage performs a preliminary delineation of structures through use of the Level Set method [20]. The third stage refines the segmentation obtained from the Level Set method, removing erroneously segmented structures. In a final stage of processing, the segmented OCT data is characterized using stereological techniques to compute the median size of the lung air spaces. In the following sections, we describe these stages in detail. 2.1 Noise reduction In OCT, speckle arises from the coherent summation of optical wavefields backscattered from sub-resolution scatterers in the sample arm of the interferometer [21,22]. In the context of this paper, it may be considered as a form of noise that compromises the ability of automated methods to identify structural features and boundaries in OCT images. Many hardware-based methods [22–26], and software-based methods [27,28] have been proposed to perform noise reduction whilst minimizing the effect on image sharpness.
#193941 - $15.00 USD
Received 22 Jul 2013; revised 24 Sep 2013; accepted 26 Sep 2013; published 9 Oct 2013
(C) 2013 OSA 1 November 2013 | Vol. 4, No. 11 | DOI:10.1364/BOE.4.002383 | BIOMEDICAL OPTICS EXPRESS 2386
In our algorithm, anisotropic diffusion [29] was implemented to reduce speckle whilst preserving the definition of edges. Anisotropic diffusion has been demonstrated previously in the removal of speckle from retinal OCT scans [30,31]. The present implementation [32] performs diffusion with a 2D kernel that adapts the degree of smoothing, based on the location and direction of intensity gradients in the image. This kernel promotes smoothing in directions with weak intensity gradients (areas of homogeneous tissue) and minimizes smoothing across strong intensity gradients (edges). The anisotropic diffusion algorithm is formulated as a diffusive process, as shown in Eq. (1): ∇I ( x, y , t ) 2 ∂ I ( x, y , t ) = ∇. exp − ∇I ( x , y , t ) , ∂t κ
(1)
where t is the algorithm iteration number; I(x,y,t) is the OCT radial B-scan intensity at coordinate (x,y) during iteration t; and κ is a positive constant chosen as a scaling factor. The negative exponential in this equation controls the rate of change of intensities at each iteration as a function of the strength of nearby edges features. For a given OCT volume, 2D anisotropic diffusion was performed within each B-scan. Subsequently, additional smoothing was performed between adjacent B-scans using weighted averaging. This operation averaged the intensity value at every individual pixel in the B-scan with the intensity of the corresponding pixels in both adjacent B-scans. This minimizes the effect of signal loss that is apparent on one B-scan but not present on adjacent scans. The loss of signal in these instances is the result of strong attenuation of the OCT signal, which can confound the delineation of lung structures. 2.2 Automatic delineation of lung structures The automatic delineation of lung structures from the OCT image was performed using a Level Set method first introduced by Osher and Sethian [20]. Using a standard approach for Level Sets, we represent the walls of the lung structures present in a single B-scan as a subset of an implicitly defined, higher dimensional function, referred to as a level set function. In our application, the walls visible in a B-scan correspond to a 2D plane (i.e., a level set) intersecting this 3D function. The 2D Level Set algorithm adopted in this work, described in detail in [33,34], does not require that we analytically define the implicit 3D function. Instead, it is initialized with a particular level set (corresponding to an initial segmentation of the B-scan), and then iteratively refined to minimize an error function. The error function is calculated from the difference between intensity values in the anisotropically smoothed OCT data and the implicitly defined 3D level set function; and also as a function of a measurement of smoothness of the shape of the segmentation outline. The error function F is shown in Eq. (2): F (φ , c, b) = ε (φ , c, b) + vL(φ ) + μ R (φ ),
(2)
where φ is the 3D level set function; c is a set of constant values which specifies a piecewise constant model of the OCT radial B-scan, approximating the intensity values within each contour with a constant value; and b is a smoothly varying field representing scanner-specific variation in B-scan intensities (i.e., a bias field). The function ε(φ,c,b) is used as the data term in the level set formulation, dependent upon differences between the OCT intensity values and the level set function. L(φ) and R(φ) are regularization terms, which serve to smooth the contour, and smooth the level set function, respectively. The constant terms ν and µ are weights to adjust the relative impact of the regularization terms. Full details of this Level Set implementation are given in [34]. A strength of the Level Set method is the robustness with which it will converge to an optimal segmentation from a wide range of initial level set configurations. In our implementation, the level set for each B-scan was initialized as a three by three grid pattern of square shapes, with dimensions approximately one tenth the length of the image. After #193941 - $15.00 USD
Received 22 Jul 2013; revised 24 Sep 2013; accepted 26 Sep 2013; published 9 Oct 2013
(C) 2013 OSA 1 November 2013 | Vol. 4, No. 11 | DOI:10.1364/BOE.4.002383 | BIOMEDICAL OPTICS EXPRESS 2387
convergence of the Level Set algorithm, a 3D level set function is computed. The delineation of lung structures is obtained from the zero level set of this function (i.e., the set of all points for which the level set function equals zero). Figure 1 shows representative results on a radial B-scan of rat lung. Figure 1(a) shows a radial B-scan acquired with the OCT needle probe. The zero level set is illustrated in Fig. 1(b), visualized as red contours super-imposed on the smoothed OCT B-scan. The interior of each red contour shape corresponds to an air space. Structures in each B-scan are segmented independently of adjacent B-scans. Performing the algorithm in 2D on a sequence of B-scans, rather than in a 3D implementation, both significantly reduces computation and better tolerates the high variation that can occur between B-scans. Note that the needle hole appears as a central black disc of no signal in Fig. 1(a). The extent of this ‘no signal’ region is extended in Fig. 1(b) as an artifact of the anisotropic diffusion performed to reduce the effects of speckle. However, we found this artifact to have minimal impact on the subsequent quantification described in Section 2.4. 2.3 Refinement of the segmented structures The results of the level set segmentation are refined through the use of morphological opening; by the application of edge-based consistency checks; and the inclusion of physiological a priori restrictions on the upper and lower limits of air space dimensions for a given animal model. A 2D morphological opening operation [35] is first performed on each B-scan of the segmented data, to separate adjacent structures that have been erroneously merged into a single structure. Morphological opening is a standard image processing technique, which consists of applying a binary erosion to the foreground pixels, followed by a binary dilation. This has the effect of separating adjacent structures whilst minimizing any artifactual reduction in the size estimate. Subsequently, an edge-based consistency check was applied to each separate segmented structure. Segmentation of our experimental results identified that the Level Set algorithm would sometimes delineate a structure even where the edge was poorly defined. This is particularly evident at locations far from the scanning probe, where attenuation of the light beam yields a poor signal-to-noise ratio and variations in the intensity gradients were often artifacts of noise rather than indicative of the boundary of underlying structures. Examples of this artifact are illustrated in Fig. 1(c), which shows a binary thresholding of the level set function, with air spaces denoted by white. Several additional erroneously segmented regions are labeled as noisy structures n1 – n9. These poorly delineated artifacts were filtered from the final segmentation by calculating the intensity gradient of the OCT radial B-scan at each delineating pixel, and computing the proportion of edge pixels above an empirically chosen edge strength. If more than half of the edge pixels for a segmented shape were below the threshold value, then the shape would be discarded as an artifact. The result of this was to only retain lung structures with strongly delineated OCT boundaries. Finally, a physiology-based consistency check was applied to specify acceptable bounds for the sizes of air spaces (spanning from alveolus to bronchiole). The alveolar size of many species has been well documented, and so we set our thresholds to range from approximately one quarter, to three times, the alveolar diameter estimated by Lum and Mitzner [36] for both pig and rodent (using mouse in lieu of rat). As there is a broad range of structure sizes, we intentionally set a wide range of acceptable sizes to include potential disease related changes whilst excluding gross segmentation errors. Structures with dimensions outside of this range were discarded. An example of a structure outside the bounds of acceptable size is illustrated in Fig. 1(c), labeled as incorrect structure n10. This artifactual feature was most likely caused by a slight tear in the tissue due to needle probe precession. Figure 1(d) shows the refined results after the filtering algorithm, illustrating the removal of the identified artifacts from the segmentation. After segmentation of each individual B-scan, the algorithm stacked the segmented B-scans together, forming a 3D volumetric segmentation of the lung parenchyma.
#193941 - $15.00 USD
Received 22 Jul 2013; revised 24 Sep 2013; accepted 26 Sep 2013; published 9 Oct 2013
(C) 2013 OSA 1 November 2013 | Vol. 4, No. 11 | DOI:10.1364/BOE.4.002383 | BIOMEDICAL OPTICS EXPRESS 2388
Fig. 1. (a) Radial B-scan image acquired from one rotation of the OCT needle probe. (b) Zero level set (red contour corresponding to the edges of structures) superimposed on the smoothed radial OCT B-scan. (c) Segmented image obtained by binarizing the level set function and performing morphological opening. This image contains both noisy structures (n1-n9) and an incorrect structure (n10). (d) Segmented image after post-processing with noisy and incorrect structures removed by edge-based consistency checking (refer to text).
2.4 Quantification To quantify the resulting segmentation, we adapted stereological methods [6] for determining the size of air spaces. These are standard methods typically applied to a set of 2D histological slices, and adapted here to quantify individual 2D radial scans of the segmented OCT data. We utilized chord length measurements [7] to quantify segmented OCT data. Beginning with a randomly selected point located on the segmented parenchyma, we defined an arbitrarily orientated line segment extending across the adjacent air space, and calculated the distance to the next intersection with parenchyma. Note that the resulting chord often crosses from one alveolus, through an alveolar duct, and terminates in an opposing alveolus. While individual chords may span a single alveolus or a bronchiole, we found the median chord length to be more indicative of the entire acinar air space (the space including the respiratory bronchiole and adjoining alveoli), as noted in [7]. Segmentation of the OCT scans allowed this measurement to be fully automated, and a population of chord lengths was computed over each data volume. The median and standard deviation of the chord lengths were then calculated, providing a characterization of the
#193941 - $15.00 USD
Received 22 Jul 2013; revised 24 Sep 2013; accepted 26 Sep 2013; published 9 Oct 2013
(C) 2013 OSA 1 November 2013 | Vol. 4, No. 11 | DOI:10.1364/BOE.4.002383 | BIOMEDICAL OPTICS EXPRESS 2389
geometry of the air spaces. As the distribution of chord lengths is not symmetric, we report the median instead of the mean. 3. Experiment 3.1 Tissue preparation All animal experiments conformed to institutional ethics approvals. Tissue was obtained from animals sacrificed as part of a different research project approved by The University of Western Australia Animal Ethics Committee, and utilized under institutional tissue sharing protocols. Two animals were used in this research: an adult rat, and an adult pig, and imaged according to the protocol described in [13]. The animals were euthanized on the day of the experiment, and dissection was performed immediately post-mortem. Both animals were exsanguinated prior to removal of lung lobes. Scans were acquired between 1 and 4 hours after dissection. The lungs were kept moist during this period with Krebs solution to prevent dehydration. Both the rat and pig lungs were allowed to deflate passively for 1 to 2 hours after excision. After tracheal cannulation, the lungs were inflated by a syringe filled with isotonic saline to improve image quality and minimize imaging artifacts that arise from the refractive index mismatch between the air-filled alveoli and the tissue of the parenchyma, as described in [13]. 3.2 Optical coherence tomography imaging OCT images were acquired using a spectral-domain OCT scanner with a central wavelength of 836 nm and a spectral bandwidth of 50 nm, giving a measured axial resolution of 9.3 µm in air. For the experiments performed here, the OCT scanner used an exposure time of 13.4 µs and was interfaced to an OCT needle probe, encased within a 30-gauge needle (outer diameter 310 µm). The design of this side-facing OCT needle probe has been detailed previously [37], but is briefly outlined here. A length of single-mode fiber was interfaced to the sample arm of the OCT scanner, and terminated with miniaturized focusing optics. The focusing optics comprise lengths of no-core fiber and graded-index (GRIN) fiber to produce a converging light beam and are terminated with an angle-polished, gold-coated length of nocore fiber to redirect the light beam at right angles to the fiber axis. The fully fused assembly is mounted in a 30-gauge hypodermic needle using adhesive. A small imaging window was etched in the side of the needle shaft, through which light is emitted and collected. The average full width at half maximum (FWHM) transverse resolution of the probe is 12 µm in tissue, over a range of 600 µm. The needle probe was mounted onto a stepper motor and counter-rotated over a range of 360 degrees at a frequency of 1 Hz to acquire 2D radial Bscans. Consecutive radial B-scans were acquired by translating the needle in steps of 5 µm, over a pullback distance of 2 mm. This produced a 3D cylindrical volumetric data set with an imaging radius of approximately 600 µm. 3.3 Data processing and visualization The segmentation algorithm was implemented using MATLAB (Mathworks, Natick, MA, USA) run on a PC running Windows 7 64-bit, with 16 GB of RAM and a quad-core Intel i5 processor. We developed custom software to visualize the 3D segmented data. This software was written in C++ and used the Visualization Toolkit (VTK) (Kitware, NY, USA). A Marching Cubes algorithm was used to create a 3D surface from the segmented data. This rendered surface and the OCT data were then graphically overlaid. The opacity of the OCT data was modified to allow simultaneous visualization of both the segmented data and the underlying OCT data volume, allowing a visual assessment of accuracy of the segmentation algorithm. A sample of 10,000 chord lengths was calculated over each data volume, and the median and standard deviation of these measurements were computed. To validate these measurements, we additionally performed 1,000 manual measurements, where a trained
#193941 - $15.00 USD
Received 22 Jul 2013; revised 24 Sep 2013; accepted 26 Sep 2013; published 9 Oct 2013
(C) 2013 OSA 1 November 2013 | Vol. 4, No. 11 | DOI:10.1364/BOE.4.002383 | BIOMEDICAL OPTICS EXPRESS 2390
human operator selected the beginning and end of each chord, and the median and standard deviation of these were compared against the automated measurements. 4. Results Figures 2 and 3 show representative results from the segmentation algorithm, demonstrated on the rat and pig, respectively. Figure 2(a) shows a representative radial B-scan of the rat lung. As in Fig. 1(a), the needle hole appears as a black void at the center of the image, surrounded by alveoli and bronchioles. Figure 2(b) shows the corresponding segmented radial scan of Fig. 2(a), graphically representing the segmented structures as white regions. This image demonstrates a range of segmented structures, including individual alveoli and acini (i.e., structural units comprising the respiratory bronchioles, alveolar ducts and alveolar sacs). A 3D rendering of the segmentation results is shown in Fig. 2(c) (Media 1), overlaid on the original OCT data.
Fig. 2. Segmentation results obtained from the OCT image of rat lung. (a) Radial OCT image acquired from one rotation of the needle probe. (b) Segmented radial image corresponding to the OCT image in (a). (c) Clipped portion of the 3D segmented data overlaid on the partially transparent OCT data (Media 1). Green structures correspond to structures segmented by the algorithm. Faded green structures indicate the segmented structures lying underneath the OCT data.
Figure 3(a) shows a representative radial scan acquired from the pig lung, and the corresponding segmentation in Fig. 3(b). The 3D rendered visualization of this data set is shown in Fig. 3(c) (Media 2). These visualizations show, to our knowledge, the first 3D volume renderings of lung structures acquired with an OCT needle probe. In these
#193941 - $15.00 USD
Received 22 Jul 2013; revised 24 Sep 2013; accepted 26 Sep 2013; published 9 Oct 2013
(C) 2013 OSA 1 November 2013 | Vol. 4, No. 11 | DOI:10.1364/BOE.4.002383 | BIOMEDICAL OPTICS EXPRESS 2391
visualizations, the structures segmented by the algorithm are represented as green surfaces. Only a clipped portion of the segmented data is shown to improve clarity of the visualization. In addition to the segmentation of structures, Figs. 2 and 3 also illustrate the removal of some structures introduced as an artifact of OCT needle imaging. Some tearing of tissue by the needle probe is apparent above the central needle hole in Fig. 2(a) and to the right of the hole in Fig. 3(a). These regions, which do not correspond to the natural air spaces of the lung, were filtered out by the algorithm and, therefore, do not appear in the segmented data.
Fig. 3. Segmentation results obtained from the OCT image of pig lung. (a) Radial OCT image. (b) Segmented radial image corresponding to the OCT image in (a). (c) A 3D visualization of the segmented data overlaid on the partially transparent OCT data (Media 2). Green structures correspond to air spaces segmented by the algorithm.
Table 1 illustrates the median and the standard deviation of air space sizes, measured both automatically (10,000 measurements) from the segmented data and manually (1,000 measurements) from the raw OCT data, for both the pig and rat data volumes. The automated median chord length calculated from the segmented data was 7% lower than the manual estimate for the pig lung, and 6% lower for the rat lung. These values lie within the standard deviation of the manual measurements. The standard deviation of the automatic chord length measurements differed by less than 1 µm from the standard deviation obtained from the manual measurements for both the pig and rat lung.
#193941 - $15.00 USD
Received 22 Jul 2013; revised 24 Sep 2013; accepted 26 Sep 2013; published 9 Oct 2013
(C) 2013 OSA 1 November 2013 | Vol. 4, No. 11 | DOI:10.1364/BOE.4.002383 | BIOMEDICAL OPTICS EXPRESS 2392
Table 1. Comparison of the median and standard deviation about the median of air space sizes obtained from automated and manual methods for both the pig and rat data. Median and standard deviation of air space size quantified in pig lung [median ± standard deviation] (µm) Automatic Manual % Quantification Quantification Difference 75 ± 5 80 ± 5 7
Median and standard deviation of air space size quantified in rat lung [median ± standard deviation] (µm) Automatic Manual % Quantification Quantification Difference 62 ± 5 65 ± 5 6
5. Discussion The capacity to quantify the dimensions of lung structures, either singly or multiply at different pressure levels, may enable new protocols for assessing lung diseases that affect the structure of the lung. The results presented in this paper highlight the potential for an algorithm to segment and quantify lung structures in an automated way, demonstrated in pig and rat models. The ability to segment structures with a wide range of shapes is an advantage of the Level Set method, and its representation of this complex morphology through the use of a level set function. The accompanying animations demonstrate that the automatic segmentation method achieved a useful segmentation of the OCT data. These results have been quantitatively validated against manual stereological estimates of air space size. The quantification of median size obtained from the segmentation was ~7% lower than the manual measurements for both pig and rat data. This difference can be attributed to regions of low signal-to-noise ratio in the OCT image, which were filtered out to improve robustness of the automated segmentation, but utilized in the manual measurements. However, this error translates to a disparity in chord length of ≤ 5 µm between the segmented data and the OCT data measurements, which is less than one transverse resolution element of the OCT system (12 µm). In addition, this difference is less than the observed physiological changes in the size of airway structures observed under several disease conditions. For example, COPD has been found to double the mean chord length [7] and increase the weighted mean alveolar diameter by approximately 80% [5] in comparison to healthy lung. This suggests that the segmentation is a sufficiently accurate characterization of the OCT data to distinguish healthy lung from various disease processes, and highlights the potential for diagnostically relevant quantifications to be made from this segmentation. We note that some pathological changes could present additional challenges to automated segmentation of OCT images, such as the thinning of the alveoli walls that is observed in COPD. However, extremely high sensitivity OCT imaging probes, such as that demonstrated in [38], should offer an improved ability to detect such features. The segmentation algorithm was observed to sometimes produce erroneous segmentations at locations with a low OCT signal-to-noise ratio (SNR). This is illustrated in Fig. 4, showing an erroneous segmentation of part of a bronchiole wall. Figure 4(a) shows the radial B-scan after anisotropic diffusion to reduce the effects of speckle. The zero level set, delineating the parenchymal structures, is superimposed on the image as a red contour. The bronchiole is highlighted by a yellow circle. The OCT signal from the far wall of the bronchiole has suffered significant attenuation, making the outline hard to distinguish, and resulting in the segmentation highlighted in Fig. 4(b). A small number of these structures are not removed during the post-processing described in Section 2.3, and remain in the final segmentation. Conversely, in a small number of cases, low OCT signal was found to obscure parts of the parenchyma between adjacent air spaces. This resulted in a segmentation that aggregated these air spaces. An example is shown in Fig. 4(c), where the parenchymal wall between an acinus and an area of void created by the needle probe is obscured, resulting in a segmentation which combines these into a single, connected air space. Subsequent filtering of abnormally large air spaces, typically indicative of such segmentation errors, resulted in removal of the acinus, seen in Fig. 4(d). Both of these classes of error could be reduced with the implementation of a 3D Level Set algorithm, which would combine information from
#193941 - $15.00 USD
Received 22 Jul 2013; revised 24 Sep 2013; accepted 26 Sep 2013; published 9 Oct 2013
(C) 2013 OSA 1 November 2013 | Vol. 4, No. 11 | DOI:10.1364/BOE.4.002383 | BIOMEDICAL OPTICS EXPRESS 2393
adjacent B-scans. Extending the Level Set algorithm to 3D is feasible and has previously been demonstrated [39,40], but also results in a significantly increased computational load for the algorithm.
Fig. 4. Two examples illustrating errors in segmentation of lung structures due to low OCT SNR (a) Segmented radial image after anisotropic diffusion. Red contour shows zero level set. Yellow circle highlights bronchiole with far wall poorly delineated. (b) Resulting segmentation omitting far wall of bronchiole. (c) Yellow circle highlights an acinus apparently joined to an area of void caused by the needle probe. The white arrow indicates the location at which the lung parenchyma is obscured by poor OCT SNR (d) Resulting segmentation with acinus omitted.
The Level Set method was chosen for our application as it has several advantages over alternative image segmentation techniques, such as parametric active contours [41]. Firstly, the level set function representation of structure means that complex morphology, such as the intricate network of alveolar sacs in the lung, can be handled in a robust manner [34]. This level set representation can also handle changes to this complex topology, such as the bifurcation or merging of airways, which is common in the lung [34,42]. Furthermore, the Level Set method does not require parameterization, which simplifies numerical computation. An advantage of this particular implementation of level sets is that it exhibits robustness of segmentation in the presence of some inhomogeneity of signal intensity. Using a technique originally designed for inhomogeneity in MRI images, the algorithm addresses this issue by generating a bias field and updating its estimate of the bias field with each iteration, along with the level set function, until it resembles the underlying intensity inhomogeneity in the image data.
#193941 - $15.00 USD
Received 22 Jul 2013; revised 24 Sep 2013; accepted 26 Sep 2013; published 9 Oct 2013
(C) 2013 OSA 1 November 2013 | Vol. 4, No. 11 | DOI:10.1364/BOE.4.002383 | BIOMEDICAL OPTICS EXPRESS 2394
Additional stereological measures may be computed from segmented data, such as the volume fraction of parenchymal tissue [6], total alveolar surface area [6] and total alveoli number [43]. The choice of measure depends both on the pathological variations caused by the disease, and the physiological question being posed. In addition, parametric imaging techniques provide an alternate approach to quantification of the underlying structure of the tissue. Such techniques have been demonstrated to quantify the attenuation coefficient of tissue [44,45], or automatically quantify the degree of birefringence [46]. These approaches rely on segmentation of tissue types, as both the attenuation coefficient and birefringence will vary between different tissues. We suggest that automated segmentation techniques, such as that proposed in this paper, represent one approach to improving the quantification accuracy of these parametric imaging techniques. The method reported here has potential in the assessment of lung compliance, by computing measurements on lung tissue acquired at different levels of expansion. Dynamic OCT scans of lungs undergoing simulated respiration have previously been demonstrated in [13], producing a sequence of lung scans at different stages of alveolar expansion and recruitment. Alternatively, Williamson et al. [47] demonstrated the acquisition of a series of OCT scans at predefined lung pressures. Whilst this latter study quantified changes in bronchial dimensions rather than alveolar spaces, it demonstrates the practicality of acquiring OCT scans at predefined and well-controlled pressures to characterize compliance of lung parenchyma in a range of diseases. Quantification of lung structure dimensions at each pressure, and changes across pressures, could provide an indicator of changes in tissue compliance due to tissue remodeling, which is significant in diseases such as emphysema and pulmonary fibrosis [48,49]. 6. Conclusions This paper has presented a fully automated segmentation and quantification algorithm capable of delineating large numbers of lung structures from an OCT volume. The algorithm was demonstrated on two animals: pig and rat. In both animals, the majority of lung structures visible in the OCT data were accurately segmented. Quantification of air space size was validated against manual measurements, and the automated estimates of median chord length were within 5 µm of manual measurements. This paper also presents the first published 3D volume renderings of lung structures acquired with an OCT needle probe. Our segmentation method may enable creation of new protocols for characterizing lung disease through objective and automated measurement of the air spaces, with the potential to assist in the assessment of new treatments for these diseases. Acknowledgments Robert A. McLaughlin is supported by a fellowship from Cancer Council Western Australia and Brendan F. Kennedy is supported by a Raine Priming Grant from the Raine Medical Research Foundation. This research has been supported in part by funding from the National Breast Cancer Foundation, Australia. We wish to acknowledge A/Prof. Peter Noble and Prof. J. Pillow, Centre for Neonatal Research and Education, UWA, for provision of the lung tissue used in this research. We also wish to acknowledge our associate Bryden Quirk, for his assistance in construction of the OCT needle probe used in these experiments.
#193941 - $15.00 USD
Received 22 Jul 2013; revised 24 Sep 2013; accepted 26 Sep 2013; published 9 Oct 2013
(C) 2013 OSA 1 November 2013 | Vol. 4, No. 11 | DOI:10.1364/BOE.4.002383 | BIOMEDICAL OPTICS EXPRESS 2395
