1 Mosaic Image Generation
The following paragraphs 1.1 to 1.3 of this supplemental file are a complete reprint of paragraphs 2.2 to 2.4 (including the accompanying figures) from (1), with permission of Deutsche Akademie der Naturforscher Leopoldina – Nationale Akademie der Wissenschaften. No changes were made except for adjusting the numbering of the paragraphs and the figures and reformatting the literature references. The acronym CLSM (confocal laser scanning microscopy) used in those paragraphs is synonymic to CCM (corneal confocal microscopy). Paragraph 1.4 and the complete chapter 2 have been written newly for this supplemental file.
1.1 CLSM Image Registration. The recorded image data are processed in a subsequent offline process. In a first step, every focus image stack is registered completely, i.e. the individual images are transformed into a common threedimensional coordinate system. This step is required, as the eyes move permanently during recording even in case of excellent patient compliance. As a result, successive images are shifted by a certain vector.
Moreover, eye movements also cause characteristic distortion artifacts in the images because of the way the images are acquired sequentially pixel by pixel. Motion of the eye in horizontal direction (nasaltemporal) leads to horizontal shear, motion in vertical direction (inferior-superior) causes vertical compression or strain of image information (Fig. 1). Image registration therefore has to correct both the lateral offset of the individual images of a focus image stack by appropriate translations and the movement-caused distortions.
Supplementary Figure 1
Directly adjacent individual images of a focus image stack in the area of the SNP with strong motion-induced distortion in the second image; top row: Prior to correction of motion artifacts; bottom row: After correction of motion artifacts
In principle, the absolute extent of motion-induced distortions is reduced with shorter recording intervals. Individual image lines (recorded in about 65 μs each) can be assumed to be free of distortions in good approximation. The task of image registration can then be described as determining offset vectors for all image lines, which compensate motion-caused distortions and correctly align the individual images. By the comparison of two directly succeeding individual images, partial areas can be related to each other. Ideally, this might be done for every single image line, which would directly yield the offset values of the image lines of one of the images in relation to the other image. Unfortunately, the image information of a single image line is not sufficient to make a robust correlation. Every image is therefore decomposed first into strip-shaped partial images of a certain height. With the help of the phase correlation function (2), the corresponding image area in the previous image is identified for every partial image (Fig. 2).
Supplementary Figure 2.
Schematic representation of non-linear registration; the strip-shaped partial images of the second image are matched in the best possible way with the previous image using a correlation method.
The registration result of every partial image directly relates to a position change of the eye between the known recording times of the two corresponding image areas. These relations can be translated into a system of linear equations from which the complete motion of the eye during the recording of the focus image stack can be estimated using a mathematical minimization approach (3,4). Taking into account the model for the development of motion artifacts explained above, the correction vectors of all image lines can be derived directly from the calculated motion trajectory of the eye and the motion artifacts can be corrected (Fig. 1).
1.2 SNP Reconstruction. The anterior surface of Bowman's membrane – near which the SNP is located – has a much higher reflectivity than adjacent tissue layers. This boundary layer is visible as a light area in section images along the z-axis of the volume reconstructed from the registered image stack (Fig. 3).
Supplementary Figure 3.
Section in the x-z-plane of a reconstructed volume with ACM deformations; the boundary surface between the basal epithelium and Bowman's membrane is visible as a highly reflective layer.
By tracking this highly reflective layer, starting from an automatically identified seed point, a depth map of the SNP can be calculated (5,6). Using this depth map, a two-dimensional image of the epithelial basal membrane can be reconstructed from the registered image stack. It shows the SNP over the complete area, also in case of ACM deformations (Fig. 4).
Supplementary Figure 4.
Overview of the SNP reconstruction steps; top: Five individual images of a focus stack; bottom left: Depth map calculated from the focus image stack (rendered as relief plot and with color coding of depth values); bottom right: Reconstructed projection image of the SNP.
1.3 Fusion of SNP Images. The image registration and SNP image reconstruction process described in sections 1.1 and 1.2 is carried out for all acquired image stacks of a subject. The intermediate result at this point is a set of twodimensional reconstructed SNP images with the sub-basal nerve structures visible over the entire field of view. The visible area of the SNP in each of these images is usually approximately equal to the size of the field of view of the microscope, but it can be increased or decreased, depending on the local arrangement of the ACM and the eye movements during the scan.
As the microscope is not moved during the examination of a subject and the subjects themselves are instructed to keep their gaze at a fixed point, positional changes of the imaged part of the SNP can be attributed almost exclusively to unavoidable unconscious eye movements. The reconstructed SNP images usually show a certain amount of partially overlapping areas and can therefore be combined in a larger mosaic image in principle. This first requires them to be aligned to each other and after this image registration, they are fused into the resulting mosaic image.
Although motion artifacts have already been corrected during CLSM image registration (see section 1.1), small distortions may still remain in the reconstructed SNP images. Thus, simple rigid approaches will not be sufficient for the image registration task. To compensate the remaining image distortions, we are currently using the same image registration algorithm for the reconstructed SNP images as for the CLSM image stacks.
The fusion of the registered images is finally done by calculating a weighted average of all overlapping image areas. Weights are introduced, because the signal-to-noise ratio is higher in central areas of the images and decreases towards the image borders. The weight of a pixel is then defined by a linear function of its distance to the closest image border, reaching a maximum weight of 1 at a certain distance threshold. All pixels with a border distance exceeding the threshold are equally weighted with 1. Figure 5 shows a mosaic image generated by 61 reconstructed SNP images.
Supplementary Figure 5.
Result of image fusion from 61 focus image stacks; the white frame shows the size of an individual HRT image (400 × 400 μm2) for comparison.
1.4 Implementation Details and Runtime Requirements. The algorithms described in the previous paragraphs were specifically designed for processing continuous CCM (or more accurately HRT) image sequences acquired as described in the Methods section of the paper. They are implemented entirely in C++. The described process is mostly automated, but not entirely. The tracing of the SNP layer (section 1.2) inside the reconstructed volume of a 3D stack sometimes drifts into wrong tissue layers at the edges of the volume. These regions have to be removed manually. The mosaic generation takes about 3 hours of computation time per patient, plus about 20 minutes for the manual removal of erroneous tissue.
2 Image Processing
SNP reconstruction largely eliminated contrast changes based on ACM, but its features may still be visible in the form of gentle, large-scaled brightness variations. Tonal correction, unshading and local adaptive elimination of these variations cleared the mosaic image and enhanced relevant thin image features (nerve fibers) which are depicted as sudden and narrow brightness changes in the image area. The resulting image was then convolved with a series of line-detecting Gabor kernels (7) which additionally enhanced nerve fibers. The quality of this line detection approach considerably depends on the Gabor kernel parameters, which were adjusted to occurring fiber properties and to image properties resulting from prior preprocessing steps.
The filtered image was then over-segmented using a minimum-error thresholding method (8) and subsequently skeletonized. The skeleton pixels were used as seed points in a following adaptive local thresholding process (9), that expands the skeleton based on local grayscale distributions and produces a final segmentation image containing all detected structures. In some areas of the resulting image, a variety of wrongly segmented image components may appear. On the one hand, they are induced by artifacts in the reconstructed mosaic image of the SNP (sudden changes in brightness, blurring, or doubling of image features) and on the other hand, segmentation errors are caused by low contrasting fibers or contrast-changing progressions of fibers. Other structures like Langerhans cells and fibrotic tissue may also be found in the imaged area. The segmented image had to be corrected by removing those structures based on morphological properties (size, elongation). Interrupted nerve fiber segments were then reconnected in a series of sub-steps. Hereby, the segmented nerve fiber network was skeletonized again and fiber end-points were identified. The pixel coordinates of the corresponding fiber sections (medial axes) were interpolated with B-splines (10) which were continued towards the estimated end-point direction and widened with a search cone. The seeking process advances iteratively over a defined maximum reconnection distance for all end-points in parallel and can connect two endings as well as one ending to a fiber in order to create a branch. The thickness of such reconstructed fibers was estimated on the basis of the original fiber sections involved. Finally, a pruning operation was applied in order to identify and eliminate very short branches. After this series of preprocessing steps a masking operation for image areas near the image borders was carried out in order to eliminate possible border artifacts and the resulting image was ready for quantitative morphometrical analysis.
A flow chart of all relevant image processing steps is presented in Fig. 6 and some illustrations depicting the results of the individual steps are shown in Fig. 7.
Supplementary Figure 7.
Illustration of the effects for the individual image processing steps; section of an original image (A), enhanced contrast (B), Gabor-filtering (C), segmentation (D), adaptive local thresholding and nerve fiber reconstruction (E), and graph reconstruction with skeleton (green) and branching points (yellow) (F).
As an explanatory note, no fractal analysis in order to identify primary or secondary branches was conducted and therefore no distinction was made between main fibers and any order of branches. Some recent studies (11) that were based on single CCM images with a size of 0.16 μm2 have performed such a differentiation. While nerve fibers can be classified as primary/secondary fibers mainly based on general fiber orientation in single CCM images, this becomes increasingly difficult when looking at larger sections of the SNP where nerve fiber orientations become less uniform as a radial pattern can be observed. Moreover some of the large-scale reconstructions present a rather uniform and net-like configuration without any discernible hierarchical pattern. Finally the irregular shape of the image area border makes it difficult in some cases to decide whether or not several nerve fiber segments really belong to a single nerve fiber. Because of the above reasons a definition for distinct nerve fibers and thus the calculation of our parameters was chosen that was unambiguously applicable to all possible configurations of the SNP. In the present study, a single nerve fiber is defined as a fiber segment terminated by branching points, end points, and/or image borders.
A large-scale reconstruction with an image area size of 13.92 conventional CCM images, or 2.23 mm2 respectively, took 26:50 minutes to process on a high-end desktop computer running Mathematica 9.0.1 (Wolfram Research Inc., Champaign, IL, USA). The most time (20:08 minutes) was needed for the iterative search cone method that reconnects nerve fibers (535 potentially reconnectable single fiber ends over a maximum reconnection distance of 40 pixels) because this step could not be parallelized. Actual image analysis took 5:29 minutes.
Köhler B, Allgeier S, Stachs O, Winter K, Bretthauer, G. Software-based Imaging and Quantitative Analysis of the Corneal Sub-basal Nerve Plexus. Nova Acta Leopold 2014; 119(401) (in press).
Takita K, Aoki T, Sasaki Y, Higuchi T, Kobayashi K. High-accuracy subpixel image registration based on phase-only correlation. IEICE T Fund Electr 2003; E86-A: 1925–1934.
Allgeier S, Eberle F, Köhler B, Maier S, Stachs O, Guthoff RF, Bretthauer G. Ein neues Mosaikbildverfahren zur großflächigen Darstellung des subbasalen Nervenplexus der Kornea in vivo. Biomed Tech 2010; 55(Suppl. 1): [4 pages].
Allgeier S, Köhler B, Eberle F, Maier S, Stachs O, Zhivov A, Bretthauer G. Elastische Registrierung von in-vivo-CLSM-Aufnahmen der Kornea. In Bildverarbeitung für die Medizin 2011. Handels H, Ehrhardt J, Deserno TM, Meinzer H-P, Tolxdorff T, Eds. Heidelberg, Springer, 2011, p. 149–153.
Allgeier S, Zhivov A, Eberle F, Köhler B, Maier S, Bretthauer G, Guthoff RF, Stachs O. Image reconstruction of the subbasal nerve plexus with in vivo confocal microscopy. Invest Ophthalmol Vis Sci 2011; 52: 5022–5028.
Köhler B, Allgeier S, Eberle F, Maier S, Zhivov A, Guthoff RF, Bretthauer G. Ein neues Bildverarbeitungsverfahren zur zuverlässigen Erfassung des subbasalen Nervenplexus der Kornea in vivo. Biomed Tech 2010; 55(Suppl. 1): [4 pages].
Jain AK, Farrokhnia F. Unsupervised texture segmentation using Gabor filters. Pattern Recogn 1991;24: 1167–1186.
Kittler J, Illingworth J. Minimum error thresholding. Pattern Recogn 1986; 19: 41–47.
Scheibe P, Wüstling P, Voigt C, Hierl T, Braumann UD. Inkrementelle lokal-adaptive Binarisierung zur Unterdrückung von Artefakten in der Knochenfeinsegmentierung. In Bildverarbeitung für die Medizin 2012. Tolxdorff T, Deserno TM, Handels H, Meinzer H-P, Eds. Heidelberg, Springer, 2012, p. 189–194.
Dierckx P. Curve and Surface Fitting with Splines. Oxford, Oxford University Press, 1993.
Petropoulos IN, Alam U, Fadavi H, Asghar O, Green P, Ponirakis G, Marshall A, Boulton AJ, Tavakoli M, Malik RA. Corneal nerve loss detected with corneal confocal microscopy is symmetrical and related to the severity of diabetic polyneuropathy. Diabetes Care 2013; 36: 3646–3651.
The authors thank the staff of the GDS, especially M. Behler and M. Schroers-Teuber, for their excellent technical support.
The GDS was initiated and financed by the German Diabetes Center, which is funded by the German Federal Ministry of Health (Berlin, Germany), the Ministry of Innovation, Science, Research and Technology of the state North Rhine-Westphalia (Düsseldorf, Germany), and grants from the German Federal Ministry of Education and Research (BMBF) to the German Center for Diabetes Research (DZD e.V.). Parts of this work were supported by a grant from the Ministry of Science, Research and the Arts of Baden-Württemberg (AZ 33-7533-7-11.6-9/3/1).
Duality of Interest
No potential conflicts of interest relevant to this article were reported.
D.Z. designed the study. D.Z., N.P., A.Z., S.A., K.W., I.Z., J.B., A.S., S.P., B.K., and O.S. researched data. D.Z. and N.P. wrote the manuscript. R.F.G. and M.R. contributed to discussion. D.Z., N.P., A.Z., S.A., K.W., I.Z., J.B., A.S., S.P., B.K., O.S., R.F.G., and M.R. reviewed and edited the manuscript. D.Z. is the guarantor of this work and, as such, had full access to all of the data in the study and takes responsibility for the integrity of the data and the accuracy of the data analysis.
Diabetes. 2014;63(7):2454-2463. © 2014 American Diabetes Association, Inc.