Journal of Microscopy, Vol. 258, Issue 3 2015, pp. 223–232
doi: 10.1111/jmi.12236
Received 10 September 2014; accepted 3 February 2015
Seamless stitching of tile scan microscope images F . B . L E G E S S E ∗ , †, O . C H E R N A V S K A I A ∗ , †, S . H E U K E †, T . B O C K L I T Z ∗ , †, T . M E Y E R ∗ , J . P O P P ∗ , † & R . H E I N T Z M A N N ∗ , †, ‡ ∗ Leibniz Institute of Photonic Technology (IPHT) Jena e.V., Jena, Germany
†Institute of Physical Chemistry and Abbe Center of Photonics, Friedrich-Schiller University Jena, Helmholtzweg, Jena, Germany ‡King’s College London, Randall Division of Cell and Molecular Biophysics, NHH, Guy’s Campus, London SE1 1UL, U.K.
Key words. Image blending, image compositing, image stitching.
Summary
Introduction
For diagnostic purposes, optical imaging techniques need to obtain high-resolution images of extended biological specimens in reasonable time. The field of view of an objective lens, however, is often smaller than the sample size. To image the whole sample, laser scanning microscopes acquire tile scans that are stitched into larger mosaics. The appearance of such image mosaics is affected by visible edge artefacts that arise from various optical aberrations which manifest in grey level jumps across tile boundaries. In this contribution, a technique for stitching tiles into a seamless mosaic is presented. The stitching algorithm operates by equilibrating neighbouring edges and forcing the brightness at corners to a common value. The corrected image mosaics appear to be free from stitching artefacts and are, therefore, suited for further image analysis procedures. The contribution presents a novel method to seamlessly stitch tiles captured by a laser scanning microscope into a large mosaic. The motivation for the work is the failure of currently existing methods for stitching nonlinear, multimodal images captured by our microscopic setups. Our method eliminates the visible edge artefacts that appear between neighbouring tiles by taking into account the overall illumination differences among tiles in such mosaics. The algorithm first corrects the nonuniform brightness that exists within each of the tiles. It then compensates for grey level differences across tile boundaries by equilibrating neighbouring edges and forcing the brightness at the corners to a common value. After these artefacts have been removed further image analysis procedures can be applied on the microscopic images. Even though the solution presented here is tailored for the aforementioned specific case, it could be easily adapted to other contexts where image tiles are assembled into mosaics such as in astronomical or satellite photos.
In medical diagnostic applications, optical imaging techniques are required to allow fast, noninvasive, image acquisition combined with high resolution and deep tissue penetration to detect endogenous disease markers as well as morphological tissue alterations. In the diagnostic routine, large tissue areas, with sizes ranging from millimetres up to several centimetres, have to be imaged. Extended biomedical samples, however, do not fit into the field of view of a microscope for the required magnification and numerical aperture. Thus, ex vivo images of large tissue sections are measured stepwise using motorized stages that translate samples in lateral directions. Series of single tiles are obtained and rearranged into a large image mosaic for the whole specimen thereby creating a large overview image with subcellular resolution. Frequently, the visual or automated analysis of these mosaics is compromised by various artefacts. In particular, the uneven illumination and collection efficiency as well as scattering caused by the investigated tissue result in visible brightness variation within each tile of the mosaics. Uneven brightness becomes more pronounced if larger tile sizes, which allow faster scanning, are selected. These are particularly strong in nonlinearoptical-microscopy techniques such as coherent anti-Stokes Raman scattering (CARS) (Cheng et al., 2002; Nan et al., 2006), which proved great diagnostic potential (Meyer et al., 2013a) and is the modality used for capturing the test images in this work. The nonlinearity and two-colour nature of the CARS process increases the susceptibility of the images to laser intensity fluctuations and chromatic aberrations such that the observed brightness changes vary significantly between tiles. Correcting these image artefacts simplifies any automated or manual image interpretation and therefore improves not only the image appearance, but its diagnostic value.
Correspondence to: Thomas Bocklitz. Tel: +49-3641/9-48328; fax: +49-3641/948302; e-mail: [email protected]
C 2015 The Authors C 2015 Royal Microscopical Society Journal of Microscopy
224
F.B. LEGESSE ET AL.
Fig. 1. Illustration of stitching variables. (A) Nonoverlapping tiles A and B representing images in the first dataset. (B) Overlapping tiles C and D representing images of the second dataset with O corresponding to the width of the overlap. In the final stitched image only half of the overlap area on each tile is used. S is the width of the grey-shaded portions of the tiles used as a basis for correction in the algorithm.
The problem of image mosaicking in the context of microscopic tile-scans can be divided into two parts as follows (Chen & Klette, 1999): 1. Image registration: is concerned with aligning matching tiles and performing tile transformations thereby determining the structure of the mosaic. 2. Image merging: refers to the process of seamlessly stitching the tiles together into a single mosaic. Some approaches further split this problem into (Brown & Lowe, 2007): (a) Nonuniform tile brightness correction: adjusting the brightness of the tiles in order to remove uneven illumination or emission collection artefacts. (b) Image blending: modifying the tile brightness so that there is a smooth transition between adjacent tiles. A number of effective approaches have been proposed for image registration. Among these are a robust method based on invariant local features (Brown & Lowe, 2007), a solution optimizing a global image quality that builds on accurate pairwise image registration (Th´evenaz & Unser, 2007) and a method that uses the weighted least squares method on image shifts calculated by cross correlation (Faas et al., 2012). However, for microscopic tile-scan acquisitions there is an adequate meta-data to determine the tile-configurations.
Furthermore, in overlapping tile acquisitions the transformation is reasonably assumed to be only translation, which is true for most microscopic acquisitions using motorized sample stages (Preibisch et al., 2009). This translational offset and in turn overlap width obtained from the meta-data can be further refined by the Phase Correlation Method (Kuglin, 1975). Therefore, the scope of the work presented here is restricted to the problem of image merging. For correcting nonuniform tile brightness during image merging, a gain compensation procedure that minimizes an error function defined as a sum of gain normalized intensity errors for all overlapping pixels has been described (Brown & Lowe, 2007; Faas et al., 2012). The downside of this approach is a requirement for a sufficient overlap to exist between the tiles. One of the most noted and simple image blending methods is based on spatially varying weighting or alpha blending. The problem of optimum seam-width selection posed by this algorithm is effectively alleviated by Multi-band (Burt & Adelson, 1983; Brown & Lowe, 2007) and Gradient domain (P´erez et al., 2003; Agarwala et al., 2004) blending schemes. These methods mostly manipulate grey levels of pixels found only in the vicinity of the boundary between neighbouring tiles. These result in a smooth transition between neighbouring tiles along the boundary but does not account for large, overall illumination difference between tiles that may lead to C 2015 The Authors C 2015 Royal Microscopical Society, 258, 223–232 Journal of Microscopy
STITCHING OF TILE SCAN IMAGES
225
Fig. 2. Adjustment of boundary curve values at corner points to corner averages. (A) Mean intensity curve H (y). The red dots indicate connection points values H (O /2) and H (L − O /2). The vertical dotted lines le f t and ri ght indicate the required correction at both ends to reach corner averages C T1 and C T2 . (B) Correcting function U (y) to subtract. The intensity values of the line at connection points correspond to the required adjustment(S(O /2) = le f t and S(L − O /2) = ri ght ). (C) Modified mean intensity curve Hm (y) = H (y) − U (y).
sharp intensity changes at the edges of the seam (see ‘Results and discussion’ section) (Uyttendaele et al., 2001). Here, we present a novel blending algorithm that removes mosaicking artefacts and corrects the whole image by computing weights that account for both the distance from tile boundaries and the brightness changes between the neighbouring tiles. The method is developed to optimally merge a rectangular grid of tiles of the same size obtained using laser scanning nonlinear microscopes equipped with a motorized stage that moves the sample laterally. Our method works for both nonoverlapping and overlapping tiles although constant vertical and horizontal overlap widths throughout the grid are currently required. The algorithm is implemented in MATLAB (The MathWorks, MA, USA) using the DIPimage toolbox (Quantitative Imaging Group, TU Delft, The Netherlands).
dataset consists of images captured by a Zeiss LSM 510 Meta microscope that feature no overlap area between adjacent tiles. The second dataset consists of CARS images that were acquired by a custom-built laser scanning microscope previously described in Meyer et al. (2013b). In contrast to the first set, these images feature a significant amount(∼30%) of overlap between neighbouring tiles. The motorized stage for moving microscopic slides has an accuracy of ±3μm as stated by the manufacturer (M¨arzh¨auser Wetzlar, Wetzlar, Germany). The size of the tiles are 2048 × 2048 and 4096 × 4096 pixels corresponding to sample areas of approximately 450μm × 450μm and 1mm × 1mm for the first and second dataset, respectively. In order to improve the signal-to-noise ratio (SNR) in the tiles, a 2 × 2 pixel binning was applied. Correcting nonuniform brightness
Methodology Image dataset The stitching algorithm was tested on two sets of CARS mosaic images, i.e. images with and without mutual overlap between neighbouring tiles. The CARS images are acquired by laser scanning microscopy of large biopsies taken from head and neck squamous cell carcinoma (HNSCC). The first C 2015 The Authors C 2015 Royal Microscopical Society, 258, 223–232 Journal of Microscopy
The removal of the uneven illumination and emission collection is based on a multiplicative correction algorithm. In this regard, each tile is multiplied with the inverse of a map reflecting deviation of a pixel value from the mean of an average tile brightness image that is constructed from foreground objects. Thus, the first step in the correction is the subtraction of the background. The applied background or offset value is provided as a predefined argument or determined as a mean value
226
F.B. LEGESSE ET AL.
Table 1. Quality metrics of image mosaics. The numbers indicate a decrease in quality of the mosaics with subsequent re-stitching. Image Mosaic1 Mosaic2 Mosaic3
MSE
PSNR
0.0017 0.0026 0.0040
27.80 25.82 24.00
Table 2. Quality metrics of image Mosaics stitched from modified tiles. The lower quality of the mosaics as compared to the values in Table 1 is due to the artificially introduced smooth intensity pattern. These mosaics display an improvement in quality over the mosaic before the application of the stitching algorithm, which has MSE of 0.0251 and PSNR of 15.19. Image Mosaic1 Mosaic2 Mosaic3
MSE
PSNR
0.0199 0.0204 0.0207
16.12 16.58 15.67
from a dark sample area selected by a user. This background value is subtracted from each of the tiles and the resulting images are compared with a certain threshold value to separate the foreground objects from the background. Stated mathematically, for a set of n tiles T = {t1 (x, y), t2 (x, y), . . . , tn (x, y)} and an estimated background value b, the set of background subtracted tiles G = {g 1 (x, y), g 2 (x, y), . . . , g n (x, y)} is given by: g i (x, y) = ti (x, y) − b,
where v is the chosen threshold value. The arithmetic mean of all these foreground areas yields an average brightness image. This average image, m(x, y), is computed from F using the expression: n f i (x, y) m(x, y) = n i =1 , (3) i =1 H [ f i (x, y)] H [a] =
1 0
if a > 0, . otherwise
(4)
Curves that reveal the general illumination pattern in this image are constructed by averaging pixel brightness values along the x and y axis. For tiles of size M × N pixels, the x − axis and y − axis means are given by the expressions: mx (x) =
N 1 m(x, y), N y=1
m y (y) =
Blending No Blending Linear Blending Average MosaicJ ICE Our Algorithm without an uneven illumination correction Our Complete Algorithm
M 1 m(x, y). (5) M x=1
MSE
PSNR
0.0048 0.0050 0.0053 0.0051 0.0040 0.0048
22.43 22.28 22.12 22.17 22.56 22.78
0.0039
24.13
Herein the assumption that no value in m(x, y) is zero was made. Smoothing these curves with a Gaussian filter (kernelsize σ = 8px) and dividing them by an average brightness value for the overall image, results in two vectors that reflect relative deviations from the mean along each axis. These deviation functions are given by the expressions:
sx (x) =
a 1 w(t, σ )mx (x + t), mavg t=−a
s y (y) =
(1)
where i = {1, . . . , n}. The corresponding set of foreground objects F = { f 1 (x, y), f 2 (x, y), . . . , f n (x, y)} is computed as: g (x, y) if g i (x, y) > v, (2) f i (x, y) = i 0 otherwise
for
Table 3. Quality metrics for mosaics stitched from overlapping subtiles. As before, a reasonable comparison of the metrics can only be made within this table because a different smooth pattern is applied for each dataset corresponding to each of the tables.
a 1 w(t, σ )m y (y + t), mavg t=−a
(6)
where w(t, σ ) is a one-dimensional Gaussian filter of size k = 2a + 1 and standard deviation σ , and mavg denotes the average brightness value for the entire image given as: ⎞ M N 1 1 1 = ⎝ mx (x) + m y (y)⎠ . 2 M x=0 N y=0 ⎛
mavg
(7)
The deviation map is computed by forming the outer product of the row x-deviation vector with the column y-deviation vector as formulated in Eq. (8): sxy (x, y) = sx (x) × s y (y).
(8)
The inverse of this deviation map results a correction map, c xy (x, y) = 1/sxy (x, y), that is multiplied with each of the background subtracted tiles to correct for the major part of the uneven illumination and emission collection. The stack of corrected images GC = {G 1 (x, y), G 2 (x, y), . . . , G n (x, y)} then becomes: G i (x, y) = c xy (x, y) × g i (x, y).
(9)
C 2015 The Authors C 2015 Royal Microscopical Society, 258, 223–232 Journal of Microscopy
STITCHING OF TILE SCAN IMAGES
227
The averaging and the subsequent smoothing with a Gaussian filter for a specific background corrected tile A(x, y) = G i (x, y) are given as: O S
+2 a 2 1 At (y) = w(t, σ ) S A(x, y + t) , t=−a
Al (x) =
a
w(t, σ )
1 S
t=−a
Ab (y) =
a
w(t, σ )
1 S
t=−a
Ar (x) =
a
w(t, σ )
t=−a
Fig. 3. Correction map construction. The value of I P is computed using inverse distance weighting (IDW) of I T ,I B ,I L , and I R which are intensity values on the top ( Atm (y)), bottom (Abm (y)), left (Alm (x)), and right (Ar m (x)) correction curves closest to the point. Table 4. Parameters of the images used in the stitching example in Figure 7.
Binned tile size Number of tiles Horizontal overlap Vertical overlap Stitch width (S) Stitched tile size IDW exponent p Kernel Size σ
Figures 7(A)–(C)
Figures 7(G)–(I)
2048 × 2048 px 7×8 634 px 610 px 100 px 1414 × 1438 px 6 30
1024 × 1024 px 9×7 0 px 0 px 5 px 1024 × 1024 px 6 30
Removing tile artefacts A visible edge appears between adjacent tiles in large image mosaics due to grey level differences across tile boundaries. The algorithm presented below modifies the tile overall intensity variations in order to provide a smooth grey level transition between neighbouring tiles. The meanings of the variables that are used in the subsequent explanation of the algorithm are introduced in Figure 1 displaying the stitching of two horizontally neighbouring tiles without (Fig. 1 A) and with overlap (Fig. 1 B). The grey-shaded segments of width S in Figure 1(A) are used by the stitching algorithm to compute mean curves along the edges. Displayed in Figure 1(B), L and O represent the width of the tiles and the overlap, respectively. In this case, the grey-shaded segments are located symmetrically around the centre of the overlap. The values of these and other parameters used in this section are detailed in Table 4. C 2015 The Authors C 2015 Royal Microscopical Society, 258, 223–232 Journal of Microscopy
x= O2 − 2S
1 S
+2
O 2
S
A(x + t, y) ,
y= O2 − 2S L − O2 + 2S
A(x, y + t) ,
x=L − O2 − 2S L − O2 + 2S
(10)
A(x + t, y) ,
y=L − O2 − 2S
where At (y), Ab (y), Al (x), and Ar (x) denote mean intensity curves located at the top, bottom, left and right edges of the tile, respectively. w(t, σ ) is a one-dimensional Gaussian filter of size k = 2a + 1 and standard deviation σ . Next, the corresponding mean curves between neighbouring edges are mutually averaged. The horizontal boundary curve for vertically adjacent tiles A and B and the vertical boundary curve for horizontally adjacent tiles A and C are given, respectively, by the expressions: H2 (y) = V2 (x) =
Ab (y)+Bt (y) , 2 Ar (x)+C l (x) . 2
(11)
Intensity values of four adjacent averaged boundary curves at corner points are averaged to define each corner target value. All individual boundary curves are then forced to share these target values at the connection points. The modification of boundary curves to adjust neighbouring connection point values to corner averages is illustrated in Figure 2. In Figure 2(A) the values of the curve H (y) at connection points y = O /2 and y = L − O /2 have to be changed to computed corner average values C T1 and C T2 , respectively. This is achieved by constructing a straight line, U (y), which features the values U (O /2) = le f t = H (O /2) − C T1 and U (L − O /2) = ri ght = H (L − O /2) − C T2 as illustrated in Figure 2(B). Subtracting U (y) from H (y) results in a new curve Hm (y), shown in Figure 2(C), with Hm (O /2) = C T1 and Hm (L − O /2) = C T2 . The remaining part of the curve is changed linearly with distance from the connection point. Dividing these modified curves by the original curves yields correction curves for each edge of tile A: Atm (y) =
H1m (y) , At (y)
Abm (y) =
H2m (y) , Ab (y)
Alm (x) =
V1m (x) , Al (x)
Ar m (x) =
V2m (x) . Ar (x)
(12)
A correction map for the whole tile is constructed from the correction curves by utilization of an Inverse Distance Weighting (IDW) scheme. IDW is an interpolation technique that assigns a value to an unknown point using the weighted average
228
F.B. LEGESSE ET AL.
Fig. 4. Blending of two vertically adjacent overlapping tiles with different methods. (A) Stitched image with no blending applied and half the overlap area from each tile removed. (B) Averaging the pixel grey-levels in the overlap area produces sharp brightness changes at the borders of the overlap. (C) The mosaic stitched with alpha blending has a smooth grey-level transition at the boundary but the overall brightness difference between the tiles is still apparent. (D) Mosaic stitched with the algorithm used in this paper is free from the stated artefacts.
of known values. The weights in the scheme are functions of the inverse of the distance from the known points to the point in question. Following Shepard (1968), these weights are power functions of the inverse distance. A schematic of the construction of a correction map is displayed in Figure 3. For the scenario shown in Figure 3, the sample points to be used are values on the edge correction curves at coordinates determined by the point in question(Iq ) and the IDW translates into: Ii Wi i , (13) Iq = Wi
1. A single CARS image of size 2048 × 2048 pixels, with intensity rescaled from 0.0 to 1.0, acquired using our custom-built LSM is cut into 16 nonoverlapping tiles of size 512 × 512 pixels. These tiles are then stitched back into a single image using our algorithm [Mosaic1]. 2. The reconstructed image in experiment 1 (Mosaic1) is then re-cut at the previous locations into tiles of the same size and then re-stitched using our method [Mosaic2]. 3. A portion of the reconstructed image in experiment 1 (Mosaic1), with size 2045 × 2045 pixels, is cut into 25 tiles of size 409 × 409 pixels and re-stitched with our algorithm [Mosaic3].
where Wi = 1/(d i ) p for i = T , B, L , and, R; d i is the shortest distance from each edge and p is a positive real number. Each of the tiles is multiplied by its unique correction map. If there is an overlap, the stitching process is completed by removing half of the overlap area from each tile.
The mosaics resulting from each of the experiments are then compared to the original image using image mosaicking error metrics defined in (Bevilacqua et al., 2010). Specifically, the metrics used here are the Mean Square Error (MSE) and the Peak Signal to Noise Ratio (PSNR) given by the expressions:
Results and discussion The algorithm for removing the tile artefacts accounts for an overall nonuniform brightness of tiles in addition to the grey level jumps at the boundary between the tiles. This advantage of the algorithm over other blending mechanisms is illustrated in Figure 4 which shows the stitching of two vertically adjacent, overlapping tiles. For objectively measuring the quality of the stitching using our algorithm the following set of experiments were performed.
MS E =
x
P S N R = 10 log10
y (R I (x,y)−RM (x,y))
N
2
,
2 N(ma x(RM )−mi n(RM )) 2 , x y (R I (x,y)−RM (x,y))
(14)
where RI (x, y) and RM (x, y) are corresponding pixels values of a region of interest R in the original and mosaic image, respectively, and N is the total number of pixels in R. The entire image was chosen as a region of interest R for our calculations. The square root of the MSE gives the average difference in grey value per pixel between the original and the mosaic image. The PSNR on the other hand gives the ratio of the maximum signal value to the MSE in the logarithmic C 2015 The Authors C 2015 Royal Microscopical Society, 258, 223–232 Journal of Microscopy
STITCHING OF TILE SCAN IMAGES
229
Fig. 5. Stitching of tiles modified to simulate real image acquisition scenarios. (A) Original image. (B) A mosaic assembled from tiles with Gaussian intensity distribution to be used for introducing vignetting artefacts. (C) Image assembled from the modified tiles. (D) The modified image stitched with the algorithm presented in this paper showing the removal of the artefact.
decibel scale. Thus, the better the quality of the stitching, the lower the MSE, the higher the PSNR. Before the computation of the metrics the total intensity of the mosaics are normalized to the total intensity of the original image. The results of these experiments is summarized in Table 1. The apparent differences between the original image and the mosaics are due to (1) modification of the original tile values by the correction map introduced to bring the different average tile brightness to a single mean value, and (2) further changes in the tile values during smooth blending the tiles by forcing the brightness values at the corners to an average value. The quality of the mosaicking also decreases, albeit at a decreasing rate, with more re-stitching as more multiplicative factors are introduced into the operations. However, it should be noted that by omitting the first precorrection step, the main algorithm for removing tile artefacts will not modify the forecast perfect brightness for overlapping tiles. The correction steps of the stitching algorithm become more relevant when there is an actual uneven brightness in the C 2015 The Authors C 2015 Royal Microscopical Society, 258, 223–232 Journal of Microscopy
tiles and an intensity jump across their boundaries. In order to replicate such image acquisition scenarios, each tile was first multiplied with smooth brightness patterns before the re-stitching in experiment 1. For instance, a smooth image of size 512 × 512 pixels having a Gaussian brightness distribution with 1 at the centre and 0.5 at corners is used to simulate vignetting. Figure 5(B) shows a mosaic of such smooth brightness patterns to be applied on the original image in Figure 5(A). Figure 5(C) displays an image assembled from the modified tiles with no blending algorithm applied. Our algorithm removes the introduced artefact as illustrated in Figure 5(D). Table 2 presents the corresponding error metrics computed for these images for the previously explained experiments. For these experiments, the total intensity of stitched image is normalized to the total intensity of the image with no blending. In order to put these results into perspective, the MSE and PSNR of the mosaic before the application of the blending algorithm as compared to the original image are 0.0251 and 15.19, respectively. An improvement in image quality can be
230
F.B. LEGESSE ET AL.
Fig. 6. An image assembled from overlapping tiles using different algorithms. (A) The image with no blending applied with half the overlap area from either side discarded features intensity jumps across tile boundaries. (B) Stitching with Microsoft’s Image Composite Editor (ICE) removes most of the tiling artefacts. The image stitched with the Grid/Collection stitching plugin in ImageJ with the blending options Average (C) and Linear Blend (D) show only a small improvement over the one with no blending. (E) ‘MosaicJ’ leads to similar results as the ‘Linear Blend’ option. (F) The image stitched with our algorithm has apparently significant improvement in quality.
attested visually from Figure 5 and quantitatively (lower MSE and higher PSNR) from Table 2. To compare our algorithm to other competing methods the following experiment was performed. Because the methods available only work on overlapping tiles, the original image in Figure 5(A) is first cut into 16 overlapping tiles of size 587 × 587 pixels. These tiles feature an overlap of 100 pixels in both the horizontal and vertical direction. As before, the sub-tiles obtained are multiplied with the smooth Gaussian brightness. The comparative stitching methods tested are. 1. The Grid/Collection stitching (Preibisch et al., 2009) plugin in ImageJ (Schneider et al., 2012) using the blending options: ‘Average’ and ‘Linear Blending’. 2. The MosaicJ (Faas et al., 2012) plugin in ImageJ. 3. Microsoft’s Image Composite Editor (ICE; Microsoft Research). The total intensity of each mosaic image is normalized to the total intensity of the stitched image with no blending obtained by just removing half the overlap area from each tile.
Figure 6 shows images stitched from overlapping tiles using these methods and our algorithm. Table 3 presents the corresponding error metrics for these mosaics as compared to the original image. As can be inferred from the results our algorithm gives lower MSE value and higher PSNR values than the other methods indicating better quality or similarity with the original images. To demonstrate the general performance of the algorithm on real tile acquisitions, Figure 7 displays the stitching of example images from two datasets. Figures 7(A)–(F) display a mosaic with overlapping tiles and Figures 7(G)–(I) display a mosaic with nonoverlapping tiles. The mosaics are presented before any correction (A, G), after correction for the uneven illumination and signal collection (B, H) and after application of our edge compensation algorithm (C, I). For the overlapping tiles the mosaics stitched using the Grid/Collection stitching (D), MosaicJ (E) and Microsoft’s ICE (F) are also presented.Important stitching parameters are summarized in Table 4. The performance of the preprocessing, and the removal of tile artefacts depend on several conditions, which – if unfulfilled – may lead to new artefacts. (1) The sample C 2015 The Authors C 2015 Royal Microscopical Society, 258, 223–232 Journal of Microscopy
STITCHING OF TILE SCAN IMAGES
231
Fig. 7. Stitching example for tiles with and without overlap. (A) A mosaic of 7 × 8 overlapping tiles. The correction of uneven illumination and signal correction in each of the tiles gives (B). (C) The stitched image with the edge artefacts between the tiles removed. For comparison purposes, also displayed are mosaics stitched with Grid/Collection stitching (D), MosaicJ (E), and ICE (F). (G) A mosaic of 9 × 7 nonoverlapping tiles. Most of the artefacts are removed by the correction of uneven brightness in (H). The stitching algorithm removes the residual visible edge artefacts (I).
possesses a reasonably homogeneous structure and a large number of foreground tiles were acquired. Such homogenous structures feature sample peaks that are not far-removed from the average of multiple tiles at the specific pixel position (I peak /Iavg < 4). Sample inhomogeneities effect in particular the calculation of the average brightness image. Thus, averaging over a large numbers of foreground tiles will ensure best correction of the uneven signal collection. However, C 2015 The Authors C 2015 Royal Microscopical Society, 258, 223–232 Journal of Microscopy
even for large numbers of foreground tiles, comparably dark foreground images with a localized increased intensity, such as CARS images with lipid droplets (bright) among connective tissue (comparatively dark), lead to undesired protrusion within the average brightness tile. Such outlier artefacts can be avoided by the computation of the median instead of the average of all foreground tiles. (2) The sample displays no periodicity in the range of the tile size. Otherwise, these sample
232
F.B. LEGESSE ET AL.
structure contributions are compensated besides the removal of image artefacts during the averaging process. Any loss of sample periodicities, however, can be circumvented by selection of an appropriate tile size. (3) All images should be acquired utilizing a stable experimental setup, i.e. no power fluctuation or beam drift should occur. Such instabilities affect both the correction of uneven signal collection and the removal of tile artefacts as exemplified in Figures 7(A)–(C). Here, the excitation power was oscillating considerably with periodicities above single tile but below the full image acquisition time. As a result of the power fluctuations the corrected image appears shadowed, but displays a seamless transition with decreased overall brightness difference between neighbouring tiles. (4) The overlap regions of neighbouring tiles cover the same area or share reasonably similar scatterer or emitter distribution. As evident from Figures 7(G)–(I) for nonoverlapping tiles, artefacts will arise if the two corresponding grey stitching areas cover regions of clearly different scatterer or emitter concentrations densities. Thus, the correction width S is required to be short for nonoverlapping images and oversampling is beneficial for the stitching result. For images that are in good agreement with all four conditions excellent result are obtained. The majority of the vignetting as well as the visible edge artefact between neighbouring tiles is removed. The morphology of the sample is readily retrieved and the image intensities strongly correlate to the concentration of scatters or emitters under investigation. Any following automated image analysis aiming for texture features or intensity based operations such as co-localization of different image modalities will benefit from the preceding correction. If, however, any of the conditions are not satisfied within reasonable approximation, the corrected images have to be interpreted with caution. In particular, the image intensity is not in strict dependence to the analyte density and quantification is, in general, not possible or may lead to wrong conclusions. Nevertheless, a visual examination of the pure sample’s morphology is simplified as evident by comparison of Figures 7(A) and (C). Although the algorithm is described here for application on images captured by nonlinear microscopy techniques, it can easily be adopted for stitching scan tile images acquired by other imaging methods such as confocal microscopy and in a general context for assembling panoramic views from multiple astronomical or satellite photos. Acknowledgements This work was financially supported by the European Union ¨ Regionale Entwicklung (EFRE)’, via ‘Europ¨aischer Fonds fur ¨ ¨ Bildung, Wissenschaft und the ‘Thuringer Ministerium fur Kultur’ (TMBWK , projects: B578-06001, 14.90 HWP and B714-07037), the European network of excellence p4l (photonics4life), the German Federal Ministry for Science and Education (BMBF) MediCARS (FKZ: 13N10774) and Fiber Health Probe (FKZ:13N12525) and the Abbe School of
Photonics scholarship, granted in cooperation with partners from the German Optics & Photonics industry, the Federal Ministry of Education and Research and the County of Thuringia. The authors are also grateful to the partners who provided the samples for the image analysis: Franziska Kluschke, ¨ Hans-Joachim R¨owert-Huber and Jurgen Lademann (Charit´e - Universit¨atsmedizin Berlin, Department of Dermatology, Venerology and Allergology). References Agarwala, A., Dontcheva, M., Agrawala, M., Drucker, S., Colburn, A., Curless, B., Salesin, D. & Cohen, M. (2004). Interactive digital photomontage. ACM T. Graphic. 23, 294–302. Bevilacqua, A., Gherardi, A. & Piccinini, F. (2010). Quantitative quality assessment of microscopic image mosaicing. WASET 47, 283–286. Brown, M. & Lowe, D.G. (2007). Automatic panoramic image stitching using invariant features. Int. J. Comput. Vision 74, 59–73. Burt, P.J. & Adelson, E.H. (1983). A multiresolution spline with application to image mosaics. ACM T. Graphic. 2, 217–236. Chen, C.Y. & Klette, R. (1999). Image stitching - comparisons and new techniques. Lect. Notes Comput. Sc. 1689, 615–622. Cheng, J., Jia, Y., Zheng, G. & Xie, X. (2002). Laser-scanning coherent antiStokes Raman scattering microscopy and applications to cell biology. Biophys. J. 83, 502–509. Faas, F.G., Avramut, M.C., van den Berg, B.M., Mommaas, A.M., Koster, A.J. & Ravelli, R.B. (2012). Virtual nanoscopy: generation of ultra-large high resolution electron microscopy maps. J. Cell Biol. 198, 457–469. Kuglin, C. (1975). The phase correlation image alignment method. In Proc. Int. Conf. Cybernetics and Society, Sept. 1975, 163–165. Meyer, T., Baumgartl, M., Gottschall, T., Pascher, T., Wuttig, A., ¨ Matth¨aus, C., Romeike, B.F., Brehm, B.R., Limpert, J., Tunnermann, A., Dietzek, B., Schmitt, M. & Popp, J. (2013a). A compact microscope setup for multimodal nonlinear imaging in clinics and its application to disease diagnostics. Analyst 138, 4048–4057. Meyer, T., Chemnitz, M., Baumgartl, M., Gottschall, T., Pascher, T., ¨ Matth¨aus, C., Romeike, B.F., Brehm, B.R., Limpert, J., Tunnermann, A., Schmitt, M., Dietzek, B. & Popp, J. (2013b). Expanding multimodal microscopy by high spectral resolution coherent anti-Stokes Raman scattering imaging for clinical disease diagnostics. Anal. Chem. 85, 6703–6715. Nan, X., Potma, E. & Xie, X. (2006). Nonperturbative chemical imaging of organelle transport in living cells with coherent anti-Stokes Raman scattering microscopy. Biophys. J. 91, 728–735. P´erez, P., Gangnet, M. & Blake, A. (2003). Poisson image editing. ACM T. Graphic. 22, 313–318. Preibisch, S., Saalfeld, S. & Tomancak, P. (2009). Globally optimal stitching of tiled 3d microscopic image acquisitions. Bioinformatics 25, 1463– 1465. Schneider, C., Rasband, W. and Eliceiri, K. (2012). NIH Image to ImageJ: 25 years of image analysis. Nat. Methods 9, 671–675. Shepard, D. (1968). A two-dimensional interpolation function for irregularly-spaced data. ACM National Conference 1968, 517–524. Th´evenaz, P. & Unser, M. (2007). User-friendly semiautomated assembly of accurate image mosaics in microscopy. Microsc. Res. Tech. 70, 135– 146. Uyttendaele, M., Eden, A. & Skeliski, R. (2001). Eliminating ghosting and exposure artifacts in image mosaics. Proc. CVPR IEEE, 2, 509–516. C 2015 The Authors C 2015 Royal Microscopical Society, 258, 223–232 Journal of Microscopy
doi: 10.1111/jmi.12236
Received 10 September 2014; accepted 3 February 2015
Seamless stitching of tile scan microscope images F . B . L E G E S S E ∗ , †, O . C H E R N A V S K A I A ∗ , †, S . H E U K E †, T . B O C K L I T Z ∗ , †, T . M E Y E R ∗ , J . P O P P ∗ , † & R . H E I N T Z M A N N ∗ , †, ‡ ∗ Leibniz Institute of Photonic Technology (IPHT) Jena e.V., Jena, Germany
†Institute of Physical Chemistry and Abbe Center of Photonics, Friedrich-Schiller University Jena, Helmholtzweg, Jena, Germany ‡King’s College London, Randall Division of Cell and Molecular Biophysics, NHH, Guy’s Campus, London SE1 1UL, U.K.
Key words. Image blending, image compositing, image stitching.
Summary
Introduction
For diagnostic purposes, optical imaging techniques need to obtain high-resolution images of extended biological specimens in reasonable time. The field of view of an objective lens, however, is often smaller than the sample size. To image the whole sample, laser scanning microscopes acquire tile scans that are stitched into larger mosaics. The appearance of such image mosaics is affected by visible edge artefacts that arise from various optical aberrations which manifest in grey level jumps across tile boundaries. In this contribution, a technique for stitching tiles into a seamless mosaic is presented. The stitching algorithm operates by equilibrating neighbouring edges and forcing the brightness at corners to a common value. The corrected image mosaics appear to be free from stitching artefacts and are, therefore, suited for further image analysis procedures. The contribution presents a novel method to seamlessly stitch tiles captured by a laser scanning microscope into a large mosaic. The motivation for the work is the failure of currently existing methods for stitching nonlinear, multimodal images captured by our microscopic setups. Our method eliminates the visible edge artefacts that appear between neighbouring tiles by taking into account the overall illumination differences among tiles in such mosaics. The algorithm first corrects the nonuniform brightness that exists within each of the tiles. It then compensates for grey level differences across tile boundaries by equilibrating neighbouring edges and forcing the brightness at the corners to a common value. After these artefacts have been removed further image analysis procedures can be applied on the microscopic images. Even though the solution presented here is tailored for the aforementioned specific case, it could be easily adapted to other contexts where image tiles are assembled into mosaics such as in astronomical or satellite photos.
In medical diagnostic applications, optical imaging techniques are required to allow fast, noninvasive, image acquisition combined with high resolution and deep tissue penetration to detect endogenous disease markers as well as morphological tissue alterations. In the diagnostic routine, large tissue areas, with sizes ranging from millimetres up to several centimetres, have to be imaged. Extended biomedical samples, however, do not fit into the field of view of a microscope for the required magnification and numerical aperture. Thus, ex vivo images of large tissue sections are measured stepwise using motorized stages that translate samples in lateral directions. Series of single tiles are obtained and rearranged into a large image mosaic for the whole specimen thereby creating a large overview image with subcellular resolution. Frequently, the visual or automated analysis of these mosaics is compromised by various artefacts. In particular, the uneven illumination and collection efficiency as well as scattering caused by the investigated tissue result in visible brightness variation within each tile of the mosaics. Uneven brightness becomes more pronounced if larger tile sizes, which allow faster scanning, are selected. These are particularly strong in nonlinearoptical-microscopy techniques such as coherent anti-Stokes Raman scattering (CARS) (Cheng et al., 2002; Nan et al., 2006), which proved great diagnostic potential (Meyer et al., 2013a) and is the modality used for capturing the test images in this work. The nonlinearity and two-colour nature of the CARS process increases the susceptibility of the images to laser intensity fluctuations and chromatic aberrations such that the observed brightness changes vary significantly between tiles. Correcting these image artefacts simplifies any automated or manual image interpretation and therefore improves not only the image appearance, but its diagnostic value.
Correspondence to: Thomas Bocklitz. Tel: +49-3641/9-48328; fax: +49-3641/948302; e-mail: [email protected]
C 2015 The Authors C 2015 Royal Microscopical Society Journal of Microscopy
224
F.B. LEGESSE ET AL.
Fig. 1. Illustration of stitching variables. (A) Nonoverlapping tiles A and B representing images in the first dataset. (B) Overlapping tiles C and D representing images of the second dataset with O corresponding to the width of the overlap. In the final stitched image only half of the overlap area on each tile is used. S is the width of the grey-shaded portions of the tiles used as a basis for correction in the algorithm.
The problem of image mosaicking in the context of microscopic tile-scans can be divided into two parts as follows (Chen & Klette, 1999): 1. Image registration: is concerned with aligning matching tiles and performing tile transformations thereby determining the structure of the mosaic. 2. Image merging: refers to the process of seamlessly stitching the tiles together into a single mosaic. Some approaches further split this problem into (Brown & Lowe, 2007): (a) Nonuniform tile brightness correction: adjusting the brightness of the tiles in order to remove uneven illumination or emission collection artefacts. (b) Image blending: modifying the tile brightness so that there is a smooth transition between adjacent tiles. A number of effective approaches have been proposed for image registration. Among these are a robust method based on invariant local features (Brown & Lowe, 2007), a solution optimizing a global image quality that builds on accurate pairwise image registration (Th´evenaz & Unser, 2007) and a method that uses the weighted least squares method on image shifts calculated by cross correlation (Faas et al., 2012). However, for microscopic tile-scan acquisitions there is an adequate meta-data to determine the tile-configurations.
Furthermore, in overlapping tile acquisitions the transformation is reasonably assumed to be only translation, which is true for most microscopic acquisitions using motorized sample stages (Preibisch et al., 2009). This translational offset and in turn overlap width obtained from the meta-data can be further refined by the Phase Correlation Method (Kuglin, 1975). Therefore, the scope of the work presented here is restricted to the problem of image merging. For correcting nonuniform tile brightness during image merging, a gain compensation procedure that minimizes an error function defined as a sum of gain normalized intensity errors for all overlapping pixels has been described (Brown & Lowe, 2007; Faas et al., 2012). The downside of this approach is a requirement for a sufficient overlap to exist between the tiles. One of the most noted and simple image blending methods is based on spatially varying weighting or alpha blending. The problem of optimum seam-width selection posed by this algorithm is effectively alleviated by Multi-band (Burt & Adelson, 1983; Brown & Lowe, 2007) and Gradient domain (P´erez et al., 2003; Agarwala et al., 2004) blending schemes. These methods mostly manipulate grey levels of pixels found only in the vicinity of the boundary between neighbouring tiles. These result in a smooth transition between neighbouring tiles along the boundary but does not account for large, overall illumination difference between tiles that may lead to C 2015 The Authors C 2015 Royal Microscopical Society, 258, 223–232 Journal of Microscopy
STITCHING OF TILE SCAN IMAGES
225
Fig. 2. Adjustment of boundary curve values at corner points to corner averages. (A) Mean intensity curve H (y). The red dots indicate connection points values H (O /2) and H (L − O /2). The vertical dotted lines le f t and ri ght indicate the required correction at both ends to reach corner averages C T1 and C T2 . (B) Correcting function U (y) to subtract. The intensity values of the line at connection points correspond to the required adjustment(S(O /2) = le f t and S(L − O /2) = ri ght ). (C) Modified mean intensity curve Hm (y) = H (y) − U (y).
sharp intensity changes at the edges of the seam (see ‘Results and discussion’ section) (Uyttendaele et al., 2001). Here, we present a novel blending algorithm that removes mosaicking artefacts and corrects the whole image by computing weights that account for both the distance from tile boundaries and the brightness changes between the neighbouring tiles. The method is developed to optimally merge a rectangular grid of tiles of the same size obtained using laser scanning nonlinear microscopes equipped with a motorized stage that moves the sample laterally. Our method works for both nonoverlapping and overlapping tiles although constant vertical and horizontal overlap widths throughout the grid are currently required. The algorithm is implemented in MATLAB (The MathWorks, MA, USA) using the DIPimage toolbox (Quantitative Imaging Group, TU Delft, The Netherlands).
dataset consists of images captured by a Zeiss LSM 510 Meta microscope that feature no overlap area between adjacent tiles. The second dataset consists of CARS images that were acquired by a custom-built laser scanning microscope previously described in Meyer et al. (2013b). In contrast to the first set, these images feature a significant amount(∼30%) of overlap between neighbouring tiles. The motorized stage for moving microscopic slides has an accuracy of ±3μm as stated by the manufacturer (M¨arzh¨auser Wetzlar, Wetzlar, Germany). The size of the tiles are 2048 × 2048 and 4096 × 4096 pixels corresponding to sample areas of approximately 450μm × 450μm and 1mm × 1mm for the first and second dataset, respectively. In order to improve the signal-to-noise ratio (SNR) in the tiles, a 2 × 2 pixel binning was applied. Correcting nonuniform brightness
Methodology Image dataset The stitching algorithm was tested on two sets of CARS mosaic images, i.e. images with and without mutual overlap between neighbouring tiles. The CARS images are acquired by laser scanning microscopy of large biopsies taken from head and neck squamous cell carcinoma (HNSCC). The first C 2015 The Authors C 2015 Royal Microscopical Society, 258, 223–232 Journal of Microscopy
The removal of the uneven illumination and emission collection is based on a multiplicative correction algorithm. In this regard, each tile is multiplied with the inverse of a map reflecting deviation of a pixel value from the mean of an average tile brightness image that is constructed from foreground objects. Thus, the first step in the correction is the subtraction of the background. The applied background or offset value is provided as a predefined argument or determined as a mean value
226
F.B. LEGESSE ET AL.
Table 1. Quality metrics of image mosaics. The numbers indicate a decrease in quality of the mosaics with subsequent re-stitching. Image Mosaic1 Mosaic2 Mosaic3
MSE
PSNR
0.0017 0.0026 0.0040
27.80 25.82 24.00
Table 2. Quality metrics of image Mosaics stitched from modified tiles. The lower quality of the mosaics as compared to the values in Table 1 is due to the artificially introduced smooth intensity pattern. These mosaics display an improvement in quality over the mosaic before the application of the stitching algorithm, which has MSE of 0.0251 and PSNR of 15.19. Image Mosaic1 Mosaic2 Mosaic3
MSE
PSNR
0.0199 0.0204 0.0207
16.12 16.58 15.67
from a dark sample area selected by a user. This background value is subtracted from each of the tiles and the resulting images are compared with a certain threshold value to separate the foreground objects from the background. Stated mathematically, for a set of n tiles T = {t1 (x, y), t2 (x, y), . . . , tn (x, y)} and an estimated background value b, the set of background subtracted tiles G = {g 1 (x, y), g 2 (x, y), . . . , g n (x, y)} is given by: g i (x, y) = ti (x, y) − b,
where v is the chosen threshold value. The arithmetic mean of all these foreground areas yields an average brightness image. This average image, m(x, y), is computed from F using the expression: n f i (x, y) m(x, y) = n i =1 , (3) i =1 H [ f i (x, y)] H [a] =
1 0
if a > 0, . otherwise
(4)
Curves that reveal the general illumination pattern in this image are constructed by averaging pixel brightness values along the x and y axis. For tiles of size M × N pixels, the x − axis and y − axis means are given by the expressions: mx (x) =
N 1 m(x, y), N y=1
m y (y) =
Blending No Blending Linear Blending Average MosaicJ ICE Our Algorithm without an uneven illumination correction Our Complete Algorithm
M 1 m(x, y). (5) M x=1
MSE
PSNR
0.0048 0.0050 0.0053 0.0051 0.0040 0.0048
22.43 22.28 22.12 22.17 22.56 22.78
0.0039
24.13
Herein the assumption that no value in m(x, y) is zero was made. Smoothing these curves with a Gaussian filter (kernelsize σ = 8px) and dividing them by an average brightness value for the overall image, results in two vectors that reflect relative deviations from the mean along each axis. These deviation functions are given by the expressions:
sx (x) =
a 1 w(t, σ )mx (x + t), mavg t=−a
s y (y) =
(1)
where i = {1, . . . , n}. The corresponding set of foreground objects F = { f 1 (x, y), f 2 (x, y), . . . , f n (x, y)} is computed as: g (x, y) if g i (x, y) > v, (2) f i (x, y) = i 0 otherwise
for
Table 3. Quality metrics for mosaics stitched from overlapping subtiles. As before, a reasonable comparison of the metrics can only be made within this table because a different smooth pattern is applied for each dataset corresponding to each of the tables.
a 1 w(t, σ )m y (y + t), mavg t=−a
(6)
where w(t, σ ) is a one-dimensional Gaussian filter of size k = 2a + 1 and standard deviation σ , and mavg denotes the average brightness value for the entire image given as: ⎞ M N 1 1 1 = ⎝ mx (x) + m y (y)⎠ . 2 M x=0 N y=0 ⎛
mavg
(7)
The deviation map is computed by forming the outer product of the row x-deviation vector with the column y-deviation vector as formulated in Eq. (8): sxy (x, y) = sx (x) × s y (y).
(8)
The inverse of this deviation map results a correction map, c xy (x, y) = 1/sxy (x, y), that is multiplied with each of the background subtracted tiles to correct for the major part of the uneven illumination and emission collection. The stack of corrected images GC = {G 1 (x, y), G 2 (x, y), . . . , G n (x, y)} then becomes: G i (x, y) = c xy (x, y) × g i (x, y).
(9)
C 2015 The Authors C 2015 Royal Microscopical Society, 258, 223–232 Journal of Microscopy
STITCHING OF TILE SCAN IMAGES
227
The averaging and the subsequent smoothing with a Gaussian filter for a specific background corrected tile A(x, y) = G i (x, y) are given as: O S
+2 a 2 1 At (y) = w(t, σ ) S A(x, y + t) , t=−a
Al (x) =
a
w(t, σ )
1 S
t=−a
Ab (y) =
a
w(t, σ )
1 S
t=−a
Ar (x) =
a
w(t, σ )
t=−a
Fig. 3. Correction map construction. The value of I P is computed using inverse distance weighting (IDW) of I T ,I B ,I L , and I R which are intensity values on the top ( Atm (y)), bottom (Abm (y)), left (Alm (x)), and right (Ar m (x)) correction curves closest to the point. Table 4. Parameters of the images used in the stitching example in Figure 7.
Binned tile size Number of tiles Horizontal overlap Vertical overlap Stitch width (S) Stitched tile size IDW exponent p Kernel Size σ
Figures 7(A)–(C)
Figures 7(G)–(I)
2048 × 2048 px 7×8 634 px 610 px 100 px 1414 × 1438 px 6 30
1024 × 1024 px 9×7 0 px 0 px 5 px 1024 × 1024 px 6 30
Removing tile artefacts A visible edge appears between adjacent tiles in large image mosaics due to grey level differences across tile boundaries. The algorithm presented below modifies the tile overall intensity variations in order to provide a smooth grey level transition between neighbouring tiles. The meanings of the variables that are used in the subsequent explanation of the algorithm are introduced in Figure 1 displaying the stitching of two horizontally neighbouring tiles without (Fig. 1 A) and with overlap (Fig. 1 B). The grey-shaded segments of width S in Figure 1(A) are used by the stitching algorithm to compute mean curves along the edges. Displayed in Figure 1(B), L and O represent the width of the tiles and the overlap, respectively. In this case, the grey-shaded segments are located symmetrically around the centre of the overlap. The values of these and other parameters used in this section are detailed in Table 4. C 2015 The Authors C 2015 Royal Microscopical Society, 258, 223–232 Journal of Microscopy
x= O2 − 2S
1 S
+2
O 2
S
A(x + t, y) ,
y= O2 − 2S L − O2 + 2S
A(x, y + t) ,
x=L − O2 − 2S L − O2 + 2S
(10)
A(x + t, y) ,
y=L − O2 − 2S
where At (y), Ab (y), Al (x), and Ar (x) denote mean intensity curves located at the top, bottom, left and right edges of the tile, respectively. w(t, σ ) is a one-dimensional Gaussian filter of size k = 2a + 1 and standard deviation σ . Next, the corresponding mean curves between neighbouring edges are mutually averaged. The horizontal boundary curve for vertically adjacent tiles A and B and the vertical boundary curve for horizontally adjacent tiles A and C are given, respectively, by the expressions: H2 (y) = V2 (x) =
Ab (y)+Bt (y) , 2 Ar (x)+C l (x) . 2
(11)
Intensity values of four adjacent averaged boundary curves at corner points are averaged to define each corner target value. All individual boundary curves are then forced to share these target values at the connection points. The modification of boundary curves to adjust neighbouring connection point values to corner averages is illustrated in Figure 2. In Figure 2(A) the values of the curve H (y) at connection points y = O /2 and y = L − O /2 have to be changed to computed corner average values C T1 and C T2 , respectively. This is achieved by constructing a straight line, U (y), which features the values U (O /2) = le f t = H (O /2) − C T1 and U (L − O /2) = ri ght = H (L − O /2) − C T2 as illustrated in Figure 2(B). Subtracting U (y) from H (y) results in a new curve Hm (y), shown in Figure 2(C), with Hm (O /2) = C T1 and Hm (L − O /2) = C T2 . The remaining part of the curve is changed linearly with distance from the connection point. Dividing these modified curves by the original curves yields correction curves for each edge of tile A: Atm (y) =
H1m (y) , At (y)
Abm (y) =
H2m (y) , Ab (y)
Alm (x) =
V1m (x) , Al (x)
Ar m (x) =
V2m (x) . Ar (x)
(12)
A correction map for the whole tile is constructed from the correction curves by utilization of an Inverse Distance Weighting (IDW) scheme. IDW is an interpolation technique that assigns a value to an unknown point using the weighted average
228
F.B. LEGESSE ET AL.
Fig. 4. Blending of two vertically adjacent overlapping tiles with different methods. (A) Stitched image with no blending applied and half the overlap area from each tile removed. (B) Averaging the pixel grey-levels in the overlap area produces sharp brightness changes at the borders of the overlap. (C) The mosaic stitched with alpha blending has a smooth grey-level transition at the boundary but the overall brightness difference between the tiles is still apparent. (D) Mosaic stitched with the algorithm used in this paper is free from the stated artefacts.
of known values. The weights in the scheme are functions of the inverse of the distance from the known points to the point in question. Following Shepard (1968), these weights are power functions of the inverse distance. A schematic of the construction of a correction map is displayed in Figure 3. For the scenario shown in Figure 3, the sample points to be used are values on the edge correction curves at coordinates determined by the point in question(Iq ) and the IDW translates into: Ii Wi i , (13) Iq = Wi
1. A single CARS image of size 2048 × 2048 pixels, with intensity rescaled from 0.0 to 1.0, acquired using our custom-built LSM is cut into 16 nonoverlapping tiles of size 512 × 512 pixels. These tiles are then stitched back into a single image using our algorithm [Mosaic1]. 2. The reconstructed image in experiment 1 (Mosaic1) is then re-cut at the previous locations into tiles of the same size and then re-stitched using our method [Mosaic2]. 3. A portion of the reconstructed image in experiment 1 (Mosaic1), with size 2045 × 2045 pixels, is cut into 25 tiles of size 409 × 409 pixels and re-stitched with our algorithm [Mosaic3].
where Wi = 1/(d i ) p for i = T , B, L , and, R; d i is the shortest distance from each edge and p is a positive real number. Each of the tiles is multiplied by its unique correction map. If there is an overlap, the stitching process is completed by removing half of the overlap area from each tile.
The mosaics resulting from each of the experiments are then compared to the original image using image mosaicking error metrics defined in (Bevilacqua et al., 2010). Specifically, the metrics used here are the Mean Square Error (MSE) and the Peak Signal to Noise Ratio (PSNR) given by the expressions:
Results and discussion The algorithm for removing the tile artefacts accounts for an overall nonuniform brightness of tiles in addition to the grey level jumps at the boundary between the tiles. This advantage of the algorithm over other blending mechanisms is illustrated in Figure 4 which shows the stitching of two vertically adjacent, overlapping tiles. For objectively measuring the quality of the stitching using our algorithm the following set of experiments were performed.
MS E =
x
P S N R = 10 log10
y (R I (x,y)−RM (x,y))
N
2
,
2 N(ma x(RM )−mi n(RM )) 2 , x y (R I (x,y)−RM (x,y))
(14)
where RI (x, y) and RM (x, y) are corresponding pixels values of a region of interest R in the original and mosaic image, respectively, and N is the total number of pixels in R. The entire image was chosen as a region of interest R for our calculations. The square root of the MSE gives the average difference in grey value per pixel between the original and the mosaic image. The PSNR on the other hand gives the ratio of the maximum signal value to the MSE in the logarithmic C 2015 The Authors C 2015 Royal Microscopical Society, 258, 223–232 Journal of Microscopy
STITCHING OF TILE SCAN IMAGES
229
Fig. 5. Stitching of tiles modified to simulate real image acquisition scenarios. (A) Original image. (B) A mosaic assembled from tiles with Gaussian intensity distribution to be used for introducing vignetting artefacts. (C) Image assembled from the modified tiles. (D) The modified image stitched with the algorithm presented in this paper showing the removal of the artefact.
decibel scale. Thus, the better the quality of the stitching, the lower the MSE, the higher the PSNR. Before the computation of the metrics the total intensity of the mosaics are normalized to the total intensity of the original image. The results of these experiments is summarized in Table 1. The apparent differences between the original image and the mosaics are due to (1) modification of the original tile values by the correction map introduced to bring the different average tile brightness to a single mean value, and (2) further changes in the tile values during smooth blending the tiles by forcing the brightness values at the corners to an average value. The quality of the mosaicking also decreases, albeit at a decreasing rate, with more re-stitching as more multiplicative factors are introduced into the operations. However, it should be noted that by omitting the first precorrection step, the main algorithm for removing tile artefacts will not modify the forecast perfect brightness for overlapping tiles. The correction steps of the stitching algorithm become more relevant when there is an actual uneven brightness in the C 2015 The Authors C 2015 Royal Microscopical Society, 258, 223–232 Journal of Microscopy
tiles and an intensity jump across their boundaries. In order to replicate such image acquisition scenarios, each tile was first multiplied with smooth brightness patterns before the re-stitching in experiment 1. For instance, a smooth image of size 512 × 512 pixels having a Gaussian brightness distribution with 1 at the centre and 0.5 at corners is used to simulate vignetting. Figure 5(B) shows a mosaic of such smooth brightness patterns to be applied on the original image in Figure 5(A). Figure 5(C) displays an image assembled from the modified tiles with no blending algorithm applied. Our algorithm removes the introduced artefact as illustrated in Figure 5(D). Table 2 presents the corresponding error metrics computed for these images for the previously explained experiments. For these experiments, the total intensity of stitched image is normalized to the total intensity of the image with no blending. In order to put these results into perspective, the MSE and PSNR of the mosaic before the application of the blending algorithm as compared to the original image are 0.0251 and 15.19, respectively. An improvement in image quality can be
230
F.B. LEGESSE ET AL.
Fig. 6. An image assembled from overlapping tiles using different algorithms. (A) The image with no blending applied with half the overlap area from either side discarded features intensity jumps across tile boundaries. (B) Stitching with Microsoft’s Image Composite Editor (ICE) removes most of the tiling artefacts. The image stitched with the Grid/Collection stitching plugin in ImageJ with the blending options Average (C) and Linear Blend (D) show only a small improvement over the one with no blending. (E) ‘MosaicJ’ leads to similar results as the ‘Linear Blend’ option. (F) The image stitched with our algorithm has apparently significant improvement in quality.
attested visually from Figure 5 and quantitatively (lower MSE and higher PSNR) from Table 2. To compare our algorithm to other competing methods the following experiment was performed. Because the methods available only work on overlapping tiles, the original image in Figure 5(A) is first cut into 16 overlapping tiles of size 587 × 587 pixels. These tiles feature an overlap of 100 pixels in both the horizontal and vertical direction. As before, the sub-tiles obtained are multiplied with the smooth Gaussian brightness. The comparative stitching methods tested are. 1. The Grid/Collection stitching (Preibisch et al., 2009) plugin in ImageJ (Schneider et al., 2012) using the blending options: ‘Average’ and ‘Linear Blending’. 2. The MosaicJ (Faas et al., 2012) plugin in ImageJ. 3. Microsoft’s Image Composite Editor (ICE; Microsoft Research). The total intensity of each mosaic image is normalized to the total intensity of the stitched image with no blending obtained by just removing half the overlap area from each tile.
Figure 6 shows images stitched from overlapping tiles using these methods and our algorithm. Table 3 presents the corresponding error metrics for these mosaics as compared to the original image. As can be inferred from the results our algorithm gives lower MSE value and higher PSNR values than the other methods indicating better quality or similarity with the original images. To demonstrate the general performance of the algorithm on real tile acquisitions, Figure 7 displays the stitching of example images from two datasets. Figures 7(A)–(F) display a mosaic with overlapping tiles and Figures 7(G)–(I) display a mosaic with nonoverlapping tiles. The mosaics are presented before any correction (A, G), after correction for the uneven illumination and signal collection (B, H) and after application of our edge compensation algorithm (C, I). For the overlapping tiles the mosaics stitched using the Grid/Collection stitching (D), MosaicJ (E) and Microsoft’s ICE (F) are also presented.Important stitching parameters are summarized in Table 4. The performance of the preprocessing, and the removal of tile artefacts depend on several conditions, which – if unfulfilled – may lead to new artefacts. (1) The sample C 2015 The Authors C 2015 Royal Microscopical Society, 258, 223–232 Journal of Microscopy
STITCHING OF TILE SCAN IMAGES
231
Fig. 7. Stitching example for tiles with and without overlap. (A) A mosaic of 7 × 8 overlapping tiles. The correction of uneven illumination and signal correction in each of the tiles gives (B). (C) The stitched image with the edge artefacts between the tiles removed. For comparison purposes, also displayed are mosaics stitched with Grid/Collection stitching (D), MosaicJ (E), and ICE (F). (G) A mosaic of 9 × 7 nonoverlapping tiles. Most of the artefacts are removed by the correction of uneven brightness in (H). The stitching algorithm removes the residual visible edge artefacts (I).
possesses a reasonably homogeneous structure and a large number of foreground tiles were acquired. Such homogenous structures feature sample peaks that are not far-removed from the average of multiple tiles at the specific pixel position (I peak /Iavg < 4). Sample inhomogeneities effect in particular the calculation of the average brightness image. Thus, averaging over a large numbers of foreground tiles will ensure best correction of the uneven signal collection. However, C 2015 The Authors C 2015 Royal Microscopical Society, 258, 223–232 Journal of Microscopy
even for large numbers of foreground tiles, comparably dark foreground images with a localized increased intensity, such as CARS images with lipid droplets (bright) among connective tissue (comparatively dark), lead to undesired protrusion within the average brightness tile. Such outlier artefacts can be avoided by the computation of the median instead of the average of all foreground tiles. (2) The sample displays no periodicity in the range of the tile size. Otherwise, these sample
232
F.B. LEGESSE ET AL.
structure contributions are compensated besides the removal of image artefacts during the averaging process. Any loss of sample periodicities, however, can be circumvented by selection of an appropriate tile size. (3) All images should be acquired utilizing a stable experimental setup, i.e. no power fluctuation or beam drift should occur. Such instabilities affect both the correction of uneven signal collection and the removal of tile artefacts as exemplified in Figures 7(A)–(C). Here, the excitation power was oscillating considerably with periodicities above single tile but below the full image acquisition time. As a result of the power fluctuations the corrected image appears shadowed, but displays a seamless transition with decreased overall brightness difference between neighbouring tiles. (4) The overlap regions of neighbouring tiles cover the same area or share reasonably similar scatterer or emitter distribution. As evident from Figures 7(G)–(I) for nonoverlapping tiles, artefacts will arise if the two corresponding grey stitching areas cover regions of clearly different scatterer or emitter concentrations densities. Thus, the correction width S is required to be short for nonoverlapping images and oversampling is beneficial for the stitching result. For images that are in good agreement with all four conditions excellent result are obtained. The majority of the vignetting as well as the visible edge artefact between neighbouring tiles is removed. The morphology of the sample is readily retrieved and the image intensities strongly correlate to the concentration of scatters or emitters under investigation. Any following automated image analysis aiming for texture features or intensity based operations such as co-localization of different image modalities will benefit from the preceding correction. If, however, any of the conditions are not satisfied within reasonable approximation, the corrected images have to be interpreted with caution. In particular, the image intensity is not in strict dependence to the analyte density and quantification is, in general, not possible or may lead to wrong conclusions. Nevertheless, a visual examination of the pure sample’s morphology is simplified as evident by comparison of Figures 7(A) and (C). Although the algorithm is described here for application on images captured by nonlinear microscopy techniques, it can easily be adopted for stitching scan tile images acquired by other imaging methods such as confocal microscopy and in a general context for assembling panoramic views from multiple astronomical or satellite photos. Acknowledgements This work was financially supported by the European Union ¨ Regionale Entwicklung (EFRE)’, via ‘Europ¨aischer Fonds fur ¨ ¨ Bildung, Wissenschaft und the ‘Thuringer Ministerium fur Kultur’ (TMBWK , projects: B578-06001, 14.90 HWP and B714-07037), the European network of excellence p4l (photonics4life), the German Federal Ministry for Science and Education (BMBF) MediCARS (FKZ: 13N10774) and Fiber Health Probe (FKZ:13N12525) and the Abbe School of
Photonics scholarship, granted in cooperation with partners from the German Optics & Photonics industry, the Federal Ministry of Education and Research and the County of Thuringia. The authors are also grateful to the partners who provided the samples for the image analysis: Franziska Kluschke, ¨ Hans-Joachim R¨owert-Huber and Jurgen Lademann (Charit´e - Universit¨atsmedizin Berlin, Department of Dermatology, Venerology and Allergology). References Agarwala, A., Dontcheva, M., Agrawala, M., Drucker, S., Colburn, A., Curless, B., Salesin, D. & Cohen, M. (2004). Interactive digital photomontage. ACM T. Graphic. 23, 294–302. Bevilacqua, A., Gherardi, A. & Piccinini, F. (2010). Quantitative quality assessment of microscopic image mosaicing. WASET 47, 283–286. Brown, M. & Lowe, D.G. (2007). Automatic panoramic image stitching using invariant features. Int. J. Comput. Vision 74, 59–73. Burt, P.J. & Adelson, E.H. (1983). A multiresolution spline with application to image mosaics. ACM T. Graphic. 2, 217–236. Chen, C.Y. & Klette, R. (1999). Image stitching - comparisons and new techniques. Lect. Notes Comput. Sc. 1689, 615–622. Cheng, J., Jia, Y., Zheng, G. & Xie, X. (2002). Laser-scanning coherent antiStokes Raman scattering microscopy and applications to cell biology. Biophys. J. 83, 502–509. Faas, F.G., Avramut, M.C., van den Berg, B.M., Mommaas, A.M., Koster, A.J. & Ravelli, R.B. (2012). Virtual nanoscopy: generation of ultra-large high resolution electron microscopy maps. J. Cell Biol. 198, 457–469. Kuglin, C. (1975). The phase correlation image alignment method. In Proc. Int. Conf. Cybernetics and Society, Sept. 1975, 163–165. Meyer, T., Baumgartl, M., Gottschall, T., Pascher, T., Wuttig, A., ¨ Matth¨aus, C., Romeike, B.F., Brehm, B.R., Limpert, J., Tunnermann, A., Dietzek, B., Schmitt, M. & Popp, J. (2013a). A compact microscope setup for multimodal nonlinear imaging in clinics and its application to disease diagnostics. Analyst 138, 4048–4057. Meyer, T., Chemnitz, M., Baumgartl, M., Gottschall, T., Pascher, T., ¨ Matth¨aus, C., Romeike, B.F., Brehm, B.R., Limpert, J., Tunnermann, A., Schmitt, M., Dietzek, B. & Popp, J. (2013b). Expanding multimodal microscopy by high spectral resolution coherent anti-Stokes Raman scattering imaging for clinical disease diagnostics. Anal. Chem. 85, 6703–6715. Nan, X., Potma, E. & Xie, X. (2006). Nonperturbative chemical imaging of organelle transport in living cells with coherent anti-Stokes Raman scattering microscopy. Biophys. J. 91, 728–735. P´erez, P., Gangnet, M. & Blake, A. (2003). Poisson image editing. ACM T. Graphic. 22, 313–318. Preibisch, S., Saalfeld, S. & Tomancak, P. (2009). Globally optimal stitching of tiled 3d microscopic image acquisitions. Bioinformatics 25, 1463– 1465. Schneider, C., Rasband, W. and Eliceiri, K. (2012). NIH Image to ImageJ: 25 years of image analysis. Nat. Methods 9, 671–675. Shepard, D. (1968). A two-dimensional interpolation function for irregularly-spaced data. ACM National Conference 1968, 517–524. Th´evenaz, P. & Unser, M. (2007). User-friendly semiautomated assembly of accurate image mosaics in microscopy. Microsc. Res. Tech. 70, 135– 146. Uyttendaele, M., Eden, A. & Skeliski, R. (2001). Eliminating ghosting and exposure artifacts in image mosaics. Proc. CVPR IEEE, 2, 509–516. C 2015 The Authors C 2015 Royal Microscopical Society, 258, 223–232 Journal of Microscopy
