Image Anal Stereol 2011;30:19-30
Original Research Paper
A MULTISCALE APPROACH TO THE REPRESENTATION OF 3D IMAGES, WITH APPLICATION TO POLYMER SOLAR CELLS
Ralf Thiedmann1, Henrik Hassfeld1, Ole Stenzel1, L. Jan Anton Köster2, Stefan D. Oosterhout2, Svetlana S. van Bavel2, Martijn M. Wienk2, Joachim Loos2, Rene A.J. Janssen2 and Volker Schmidt1
1 Institute of Stochastics, Ulm University, Helmholtzstrasse 18, 89069 Ulm, Germany; 2Chemical Engineering and Chemistry, Molecular Materials and Nanosystems, Eindhoven University of Technology, 5600 MB Eindhoven, The Netherlands
e-mail: ralf.thiedmann@uni-ulm.de, henrik.hassfeld@uni-ulm.de, ole.stenzel@uni-ulm.de, l.j.a.koster@tue.nl, s.d.oosterhout@tue.nl, s.s.v.bavel@tue.nl, m.m.wienk@tue.nl, j.loos@tue.nl, r.a.j.janssen@tue.nl, volker. schmidt@uni-ulm. de
(Accepted February 19, 2011)
ABSTRACT
A multiscale approach to the description of geometrically complex 3D image data is proposed which distinguishes between morphological features on a 'macro-scale' and a 'micro-scale'. Since our method is mainly tailored to nanostructures observed in composite materials consisting of two different phases, an appropriate binarization of grayscale images is required first. Then, a morphological smoothing is applied to extract the structural information from binarized image data on the 'macro-scale'. A stochastic algorithm is developed for the morphologically smoothed images whose goal is to find a suitable representation of the macro-scale structure by unions of overlapping spheres. Such representations can be interpreted as marked point patterns. They lead to an enormous reduction of data and allow the application of well-known tools from point-process theory for their analysis and structural modeling. All those voxels which have been 'misspecified' by the morphological smoothing and subsequent representation by unions of overlapping spheres are interpreted as 'micro-scale' structure. The exemplary data sets considered in this paper are 3D grayscale images of photoactive layers in hybrid solar cells gained by electron tomography. These composite materials consist of two phases: a polymer phase and a zinc oxide phase. The macro-scale structure of the latter is represented by unions of overlapping spheres.
Keywords: adaptive thresholding, discrete skeleton, morphological smoothing, multiscale approach, representation by overlapping spheres.
INTRODUCTION
In this paper, the nanomorphologies of composite materials are considered, being blends of two different (organic and inorganic) solid phases. For instance, such materials are used in the photoactive layer of hybrid polymer-zinc oxide solar cells, where the two solid phases play the role of a polymeric electron donor, consisting of poly(3-hexylthiophene), and an inorganic zinc oxide(ZnO)-electron acceptor.
There is a great advantage of polymer solar cells due to their potentially low production costs. However, the efficiency of polymer solar cells critically depends on the intimacy of mixing of the donor and the acceptor semiconductors used in these devices to create charges as well as on the presence of unhindered percolation pathways in the individual solid phases of the composite material to transport positive and negative charges towards electrodes; see e.g., Yang and Loos (2007). In short, there is a high correlation between efficiency and morphology of the considered
solar cells. It is therefore very important to have tools at one's disposal which are suitable to analyze and model the 3D morphology of these materials quantitatively.
As a data processing step for modeling purposes, a stochastic algorithm is developed to find suitable representations of complex 3D images by systems of overlapping spheres. Besides an enormous reduction of data, this allows the application of tools from point-process theory to analyze and model the considered image data.
More precisely, we show that a morphologically smoothed version of the ZnO phase extracted from electron tomography (ET) images and, in particular, its representation by unions of overlapping spheres can be regarded as a 3D marked point pattern where the centers of spheres are the locations of points and the corresponding radii are their marks. The advantage of this off-grid representation of the macro-scale part of the nanostructure is that it can be modeled relatively easily, independent of any given resolution or grid
structure of the underlying ET images. This is the subject of a forthcoming paper where we develop a stochastic simulation model of the 3D nanostmcture; see Stenzel et al. (2011). This model is based on the multiscale approach to the description of image data with geometrically complex structure considered in the present paper
Such models help to get a better insight into the impact of morphology on the performance of polymer solar cells and, simultaneously, they can be used for virtual scenario analyses, where model-based morphologies of solar cells are simulated to identify polymer solar cells with improved nanostmctures.
Note that the algorithm developed in the present paper can be applied to various kinds of 3D data. For instance, another application of this algorithm is considered in Thiedmann et al. (2010), where the microstmcture of porous materials in electrodes of lithium-ion batteries is analyzed and modeled on the basis of its representation by unions of overlapping spheres. However, we emphasize that in Thiedmann et al. (2010) a stochastic point-process model is developed for materials in Li-ion batteries, which uses the representation by unions of overlapping spheres just as a starting point. But, Thiedmann et al. (2010) does not focus on the representation algorithm itself
Correlation of nanomorphology and functionality in polymer solar cells
The descriptive statistical analysis of ET images for hybrid polymer-ZnO solar cells accomplished in Oosterhout et al. (2009) showed that there can be drastic structural differences between the 3D morphologies of photoactive layers made from identical types of material, but having different thicknesses. Moreover, it turned out that the efficiency of these solar cells significantly depends on the thickness of photoactive layers, sometimes called 'films' for brevity. This close relationship between nanomorphology and functionality of polymer solar cells can be explained as follows. Upon exposure to light, photons are absorbed in the polymer phase and so-called 'excitons', i.e., photoexcited electron-hole pairs, evolve. Excitons are neutral quasi-particles which diffuse inside the polymer phase within a limited lifetime; see Shaw et al. (2008). If an exciton reaches the interface to the ZnO phase, it is split up into a free electron (negative charge) in the ZnO and a hole (positive charge) in the polymer phase. This process is commonly referred to as quenching, because it reduces the intrinsic fluorescent decay of the exciton in the polymer Provided that the electrons in the ZnO phase and the holes in the polymer phase reach the electrodes at the top and bottom
of the photoactive layer, respectively, photocurrent is generated. A schematic illustration of the morphology of photoactive layers in hybrid polymer-ZnO solar cells is shown in Fig. 1, where the electrodes are supposed to be parallel to the x-j-plane. For further information about polymer solar cells and the physical processes therein we refer e.g., to Brabec et al. (2008).
Fig. 1. Schematic layout of a polymer-ZnO thin film solar cell, showing the percolation of photogenerated holes (+) and electrons (—); black denotes ZnO and white polymer
Note that the extent of blending of the two materials has a large impact on the efficiency of these solar cells, because not all excitons are quenched due to their limited lifetimes. Thus, a morphology as displayed in Fig. 1, where both materials are mixed intimately, is desirable since the excitons are more likely to reach the interface and charges can be generated.
Furthermore, the existence of unhindered percolation pathways within both phases (ZnO and polymer) is crucial since the generated charges have to be transported to the electrodes throughout the phases. Because of the electric field between the electrodes, these pathways should be preferably monotonous.
Hence, to obtain solar cells with high efficiency, an intimately mixed morphology with monotonous percolation pathways for both charge carriers is desirable and should be taken into account when producing devices. Spatial stochastic models as considered in Stenzel et al. (2011), which are based on the results derived in the present paper, can be used to identify morphologies with improved efficiency by generating virtual morphologies and investigating the transport processes of electrons and excitons, respectively.
Preprocessing of image data
The stochastic algorithm to obtain an efficient representation by unions of overlapping spheres is applicable to binary images. Therefore an
adequate binarization is required first. Since different (location-dependent) brightnesses occur in the ET images considered in the present paper, we perform an adaptive thresholding procedure proposed by Yanowitz and Bruckstein (1989) and Blayvas et al. (2006).
Furthermore, since the representation of a geometrically complex set by unions of overlapping spheres works best if the boundary of the considered set is not too rough, we additionally apply a morphological smoothing. On the other hand, the functionality of the solar cells heavily depends on the fine structures, which are removed by morphological smoothing. Thus, these 'misspecified' voxels are also included into stochastic simulation models considered in Stenzel et al. (2011). They are interpreted as the 'micro-scale' nanostmcture of the binarized ET images, whereas the morphologically smoothed images represented by unions of overlapping spheres describe the nanostmcture on a 'macro-scale'.
Outline
The paper is organized as follows. First the ET data are described which are considered in this paper Then the binarization method of adaptive thresholding is presented which is followed by the explanation of morphological smoothing, and how it leads to a multi-scale representation of the binarized images. In the same section, a stochastic algorithm is introduced to represent complex 3D structures by unions of overlapping spheres. Finally, a brief summary and some concluding remarks are given.
DATA DESCRIPTION
In this section we briefly describe the 3D image data of hybrid polymer-ZnO solar cells gained by electron tomography (ET); see van Bavel et al. (2009) and Oosterhout etal. (2009). In particular, we consider images for three devices with different thicknesses of the photoactive layers: 57 nm, 100 nm, and 167 nm. These are the examplary data, to which the stochastic algorithm developed in the present paper is applied.
Electron tomography grayscale images
For each of the three different layer thicknesses, the 3D ET images are given as stacks of 2D grayscale images (being parallel to the x-j-plane, say), which are numbered according to their location in z-direction. The sizes of these images in the x-j-plane are 934 x 911 voxels for the 57 nm film, and 942 x 911 voxels for the other two thicknesses. Each voxel represents a cube with side length 0.71 nm.
Fig. 2. Cross-sectional image of the photoactive layer of a complete device (reprinted with permission from Oosterhout et al. (2009), © 2009 Nature Publishing Group).
The stacks of 2D images consist of 80, 184 and 208 slices for the 57 nm, 100 nm, and 167 nm films, respectively. These numbers of slices do not correspond directly to the (physical) thicknesses of the photoactive layers, since some slices at the bottom and top have not been used in the following. One reason for this is the relatively rough surface of films, since the thickness of the photoactive layers (polymer-ZnO) varies locally due to the complex technology of blend manufacturing; see Fig. 2. This excludes some slices near the top. Some slices at the bottom could not be used because of tomography reconstruction artefacts due to the limited number of tilting angles. This explains why the number of slices does not correspond exactly to the layer thickness. Note that the voxel size is the same for all three layer thicknesses.
Fig. 3. 2D images of hybrid polymer-ZnO solar cells with thickness of 57 nm (left), 100 nm (center) and 167nm (right).
Fig. 3 shows representative 2D slices for the three different film thicknesses, where the darker parts of the images represent the ZnO phase due to a
higher electron density compared to polymer. The images displayed in Fig. 3 indicate clear structural differences for the three films. With increasing layer thickness the separated domains of polymer and ZnO are getting finer. In particular, the thinnest film, i.e., the photoactive layer with thickness of 57 nm, features large domains of both polymer and ZnO.
Binarization by global thresholding
To apply the stochastic algorithm developed in the present paper, the 3D grayscale images have to be binarized appropriately. Binarization is necessary since we need to decide which voxels are classified as polymer and which as ZnO. An elementary approach to binarize grayscale images is to use a global threshold: voxels are set to white (polymer) if their grayscale value exceeds a certain threshold, and are otherwise set to black (ZnO).
In Oosterhout et al (2009) two different global thresholds have been considered in order to binarize the grayscale images as displayed in Fig. 3. They have been chosen in such a way that the ZnO phase can be assumed to be a subset of the union of foreground voxels (high threshold) and, vice versa, the union of foreground voxels is contained in the ZnO phase (low threshold); see Fig. 4. Thus, it can be assumed that the actual structure of the ZnO phase lies in between these two bounding sets.
Fig. 4. Binarization with global thresholds. Left: original grayscale image, center: binarization with high threshold, right: binarization with low threshold.
A first, still crude quantitative description of structural differences between the three films can be given by the volume fractions of the binarized 3D images, that is the number of foreground (i.e., black) voxels divided by the number of all voxels of the sampling window. The empirical ZnO volume fractions have been computed for the images binarized by the low and high threshold, respectively. The results are displayed in Table 1. The ZnO volume fractions computed for the 57 nm film are rather different from those of the 100 nm and 167 nm films. On the other hand, the ZnO volume fractions of the 100 nm and 167 nm films coincide which is in accordance with
the images shown in Fig. 3; see also Oosterhout et al. (2009).
Table 1. ZnO volume fraction of photoactive layers.
layer	low	high mean
thickness threshold threshold value
57 nm 100 nm 167 nm
0.098 0.133 0.128
0.172 0.295 0.293
0.135 0.214 0.211
BINARIZATION BY ADAPTIVE THRESHOLDING
It is diflflcult to find a single global threshold to binarize the ET images adequately because of the irregular brightness of these images. An example is given in Fig. 5, with the red circle indicating an underexposed domain.
Instead of considering global thresholds, a method of adaptive thresholding can be used for binarization. This method is based on techniques described in Yanowitz and Bruckstein (1989) and Blayvas et al. (2006), where the main idea is to construct a threshold surface T{x,y) which is location-dependent. If the grayscale value of the image at location {x,y) is larger than T{x,y), the voxel with coordinates {x,y) is set to white, otherwise set to black. In the special case of global thresholding, T(x,y) is a constant, whereas in adaptive thresholding, T {x,y) takes local conditions like overexposure or underexposure into account.
Fig. 5. 2D image of 100 nm data with red circle indicating an underexposed area.
Smoothing of grayscale images
Before constructing an adaptive threshold, the grayscale images need to be smoothed to reduce image noise and to make it easier to separate the objects (ZnO domains) from the background (polymer); see e.g., Jähne (2005). In particular, we applied a Gaussian 3x3x3 filter with variance <7^ = 0.4 given in voxel size. It is quite clear that such a small filter does not
remove too much structure or detail which is important for physical properties of the films.
Note that this smoothing is implemented in 3D, whereas the remaining binarization steps are executed in 2D, i.e., they are applied to every single 2D slice, being parallel to the x-y-plane. This is sufficient for adaptive thresholding because the exposure only depends on the (x, 7)-location, but not on its z-coordinate, which can be seen by considering subsequent slices of the 3D data sets. The reconstruction of the 3D tomographic images imply a high correlation of the gray-values in z-direction, i.e., regions / locations of different brightness do not change rapidly for subsequent 2D slices. In addition, the sizes of the regions with different brightness are on a much larger scale than the structure which is analyzed, i.e., the ZnO phase. Hence, locally seen, the 2D adaptive thresholding can be interpreted as a nearly global threshold.
Local maxima of gradient magnitudes
In the next step, the magnitudes of grayscale gradients are derived for the smoothed images obtained in the previous section. This is done according to a method developed in Unser (1999) and Thevenaz et al. (2000), where a spline is put through the grayscale values on the discrete lattice of voxels. More precisely, a function / : R^ ^ [0,1] with continuous domain is constructed such that f{x,y) =	^ (^-^,/-/),
where summation extends over the lattice of voxels describing the image, /(^,/) denotes the grayscale value of the voxel at position (^,/) in the lattice and ^ is a (separable) cubic B-spline, that is (p(x,y) = ß^ {x) ß^ (7), where
0
ifo< |Ai < 1, if 1 < |Ai <2,
if2< |Ai.
The gradient of / is said to be the vector of its first partial derivatives. In the resulting image of gradient magnitudes, each voxel then takes the value of the magnitude of the gradient, i.e., the length of the gradient vector.
Finally, the local maxima of the gradient image are determined. This is accomplished by comparing each voxel with its eight neighboring voxels. To avoid oversegmentation of the background {i.e., the polymer phase), the set of local gradient maxima is thinned, where the original smoothed image is binarized using a global threshold which is calculated by iteration, following Ridler and Calvard (1978). Then, all those
local maxima which belong to the background of this globally thresholded image, are deleted. The reduced set of local gradient maxima indicates voxels where the grayscale values of the original smoothed image are strongly fluctuating. These voxels are important in order to determine the values of the adaptive threshold surface.
Computing the threshold surface
Let A = {{xi,yi), i=\,...,N} denote the (thinned) set of local gradient maxima with {xj,yj) being the corresponding coordinates on the discrete lattice of voxels. For each point of A, the original grayscale value of the (smoothed) image is used as adaptive threshold, i.e., T{xi,yi) = I{xi,yi), i = \,...,N. For the remaining voxels the adaptive threshold is computed by interpolating the grayscale values of points in A, where the method of inverse-distance weighting is used as described in Yanowitz and Bruckstein (1989). Hence, considering an arbitrary voxel {x,y) ^ A, the adaptive threshold T{x,y) is calculated as a weighted sum of T {xj, yi), i=\,...,N, i.e.,
Tix,y) = —— ^(OiT{xi,yi) .
lui OJi J
Here, the summation runs over all points of A within a sufficiently large neighborhood of {x,y), in our case a square of length 100 voxels centered at {x,y), and the weight (Oj is given by the inverse distance between {x,y) and {xj,yj), i.e., (Oj =
Additionally, the threshold surface obtained in this way is shifted up or down until the resulting binarizations yield the mean volume fraction of ZnO given in Table 1. This threshold surface is then used for binarizing the solar cell data by comparing the grayscale values of the smoothed image 1= [l{x,y)) to the values of the threshold surface T = (T{x, y)), i. e., if I{x,y) > T {x,y), then the voxel with coordinates {x,y) is set to white (polymer), otherwise to black (ZnO).
Examples of binarizing the ET images by adaptive thresholding as explained above, are displayed in Fig. 6.
STOCHASTIC REPRESENTATION
In this section we state a stochastic algorithm which we developed to represent a morphologically smoothed version of the ZnO phase by unions of overlapping spheres. This algorithm has been inspired by the stochastic watershed algorithm described in
Angulo and Jeulin (2007) and Noyel et al (2007) where the locations of seeds/markers, i.e., the initial points, for the watershed transformation are chosen at random. Since the initial points for each single run are chosen at random, the result of one run cannot be ensured to be optimal. Hence, severals runs are performed and all resulting watershed lines are subsequently convoluted using kernel density techniques. The result of the kernel density estimation is an 'intensity map' in which the multiply found watershed lines are emphasized. Hence to come up with the final result, the most dominant watershed lines have to be extracted from the intensity map, wherefore several possibilities are proposed in Faessel and Jeulin (2010). Further developments on this segmentation technique are discussed in Cativa Tolosa et al (2009), where the placement of the markers is more specific to avoid an over-segmentation of the data.
The morphology of the ZnO phase obtained by adaptive thresholding of ET images is rather complex; see Fig. 6. This makes it diflftcult to represent a preferably large part of the ZnO phase by unions of spheres, where we simultaneously try to keep the number of spheres as small as possible; see Fig. 7 for a (schematic) illustration in 2D.
The reason for this is that the binarized ET images described in the previous section contain very fine structural components such as 'thin ZnO branches', i.e., thin ZnO parts connected to larger ZnO domains, 'isolated ZnO particles', i.e., small ZnO particles in the polymer domains, and 'polymeric holes', i.e., small polymeric particles inside the ZnO domains; see Fig. 8. Since it is diflflcult to represent these fine structural components by unions of suflflciently large spheres, we first apply a morphological smoothing of the ET images.
Smoothing of binarized innages
The morphological transformations which we use for smoothing the binarized ET images are 'dilation' and 'erosion' which are defined as follows. Let B C M^ be an arbitrary set and let b{s,r) = {sf e R^ : 15—5^1 < r} denote the three-dimensional sphere with center 5 G M^ and radius r > 0. Then, the dilation of B by b{o,r) is given by the Minkowski sum 50 b{o,r) = {5+5^ : se B,si G ö(o,r)}, where b{o,r) is called 'structuring element'. Furthermore, the erosion of B by b{o,r) is given by the Minkowski difference Bg b{o,r) =	^ : 5 G 5}. In general,
erosion is not the inverse of dilation. Further details regarding morphological transformations can be found e.g., in Jähne (2005), Ohser and Schladitz (2009), Serra (1982; 1988), and Soille (2003).
Fig. 6. 2D images of hybrid polymer-ZnO solar cells; left column: original grayscale images of 57 nm, 100 nm, and 167 nm films; right column: corresponding binarized images using adaptive thresholding.
Fig. 7. Schematic 2D representation of a set by union of circles
Let now B denote the ZnO phase observed in a certain sampling window W, i.e.. Be M/ C M^. One aim of morphological smoothing is to remove small isolated ZnO particles in W. This is accomplished by a so-called 'closing' of the polymer phase which leads to the set (5^0 b{o,r)) 0 b{o,r), where B^ = W\Bdenotes the complement of 5with respect to W. This means that we first apply a dilation
and then a subsequent erosion of the polymer phase B^, considering the same structuring element b{o,r) for both morphological operations. Due to different roughnesses of microstructures in the films, different values (in voxel) of the radius r have been used, that IS r = 2V3 for the 57 nm file, and r = V3 for the 100 nm and 167 nm films. The result of this closing is the elimination of small isolated ZnO particles such that all interior voxels are closer to the polymer phase than r. Additionally, some of the thin ZnO branches are removed as well.
Fig. 8. 2D image of a cutout of the 100 nm file with 'thin ZnO branches' (red), 'isolated ZnO particles' (blue) and 'polymeric holes' (yellow).
Another aim of morphological smoothing is to remove small polymeric holes within the ZnO phase. Therefore, in the third step, we dilate the preprocessed ZnO phase B =W\ {{B^®b[o,r)) Qb[o,r)), usmg the same structuring element -6(0, r) as before. Thus, as a further desirable effect, the boundary dB' of the dilated ZnO phase B' = B ® b{o^r) becomes much smoother than the boundary dB oi the original ZnO phase B; see Fig. 9. This makes it easier to represent the set B' by unions of spheres.
Fig. 9. Original image of 57 nm film split up into structural components at two different scales.
Main steps of the algorithm
We first decompose the set B' into its connected components and apply our stochastic algorithm to each of these components separately. This leads to an essential reduction of runtime. To determine the connected components of B', the algorithm of Hoshen and Kopelman (1976) has been used which works very efiiciently.
Let C d B' hQ an arbitrary connected component of the morphologically transformed ZnO phase B'. The main steps of our algorithm can be described as follows:
1)	Place a point 5 at random in the set C and mark it with the maximum radius Tmax such that the sphere '^(^j^max) is completely contained in C. Then, a candidate for the next center of a sphere is located at random in C, where this point can be accepted or rejected, according to the acceptance-rejection rule described later in this section. This process of 'sphere placement' is continued until a certain desired volume fraction of C is covered.
2)	The procedure described above is repeated several times. Then, an intensity map of locations of accepted centers is computed, using a kernel density estimation, which indicates preferred locations of centers.
3)	The local maxima of this intensity map are used to construct the final representation of C by a union of spheres.
Placement of spheres
As already mentioned above, to begin with, a point 5 is located at random in the set C and marked with the maximum radius Tmax such that the sphere -6(5, rmax) is completely contained in C. This step is implemented by sampling a uniformly distributed random variable on the sub-lattice of those voxels belonging to C. The maximum radius Tmax can efiiciently be determined by computing a distance transformation of C; see Saito and Toriwaki (1994). To avoid trivial spheres, the radius r^ax should be greater than or equal to a/3 voxel sizes, otherwise the point 5 is rejected. Then, a candidate for the next center of a sphere is located at random in C, where the radius is chosen in the same way as before. This point is rejected if the corresponding sphere only covers such parts of C which are already covered by the union of previously accepted spheres. Otherwise, the point is accepted. This process of sphere placement is continued until a certain volume fraction of C, say 99%, is covered by spheres. Alternatively, to avoid infinite loops, the sphere-placement algorithm stops after 1000 trials to
7 S
get a new center accepted. Finally, all those spheres are deleted which are completely covered by other spheres. Note that a matching of the original data by the union of spheres in 100% of the volume fraction is not necessary for the subsequent stochastic modeling using marked point processes, since there also a certain variability occurs.
Intensity map
The representation of C by unions of spheres, as described above, is random. Therefore, it is possible to obtain a more or a less efficient result, where a smaller or a larger number of spheres is needed to cover the set C; see Fig. 7. In order to achieve an efficient representation with smallest possible number of spheres, the sphere-placement procedure described above is repeated one hundred times, say, which appears to be sufficient to stabilize the computation of intensity maps for the considered morphologies. However, note that the required number of repetitions of random placements depends on the considered structures.
Subsequently, a 3D intensity map D= {D{s),s G C} is computed by kernel density estimation. In the present paper, we have used a uniform kernel with the size of one voxel for the solar cell data since it turns out that in this case, it works most efficient. This means that for any location 5 G C, the value D{s) is proportional to the number of accepted centers which are located in a certain (infinitesimal) neighborhood of 5 G C. However, note that for another application of our stochastic algorithm to represent 3D image data of electrodes used in Li-ion batteries by unions of overlapping spheres (Thiedmann et al, 2010), a Gaussian kernel with bandwidth of 1.5 voxel sizes worked best.
The intensity map D describes the spatial distribution of frequencies of locations of accepted centers, i.e., if D{s) > D{sf) for some locations s,sf e C, then it is more likely that a center is accepted at location 5 than at location An example of an intensity map is shown in Fig. 10.
Final representation by union of spheres
The intensity map described in the previous section is used to construct the final representation of the connected component C by a union of spheres. For this purpose, the local maxima of the intensity map are ranked in descending order of their intensity values. According to this ranking we then put points at these local maxima within the set C and associate them with their corresponding radii, starting with an empty configuration.
Fig. 10. Left: morphologically smoothed 2D image of 57 nm file, right: corresponding intensity map (image has been brightened up and then inverted for demonstration purposes).
In other words, the first point is put at the local maximum s\ ^ C with the highest intensity value, and the sphere	with corresponding radius
A,max such that Ö(5i,ritnax ) C C is drawn. Then the next point is put at the local maximum 52 G C of the intensity map with the second highest intensity value, and so on. This procedure is continued until 99%, say, of the volume of C is covered by spheres, or until centers of spheres are assigned to all local maxima of the intensity map. As already mentioned in the section above, there might be centers and corresponding spheres that do not contribute to the final representation of C by a union of spheres, because they are completely covered by other spheres. Therefore this last step of 'cleaning' is necessary to delete those redundant spheres.
Restriction of domain to discrete skeleton
The above developed algorithm to find an optimal representation of 3D objects by unions of overlapping spheres can in some cases be further improved w.r.t. runtime. This goal can be achieved by a restriction of the domain C for the random placement of spheres.
More precisely, we consider the discrete skeleton of the set C, which has been introduced by Lantuejoul (1978), see also Serra (1982). We call a sphere Z/^^^ a maximum sphere in C if Z/^^^ C C and Z/^^^ is tangent to at least two distinct points of the boundary dC oi C. The union of all centers of maximum spheres is the discrete skeleton of C. As stated in Lantuejoul (1978), the union of all maximal spheres equals the set C. Therefore, the idea is to restrict the random sphere placement of the developed algorithm to the discrete skeleton. This reduces the number of rejected spheres, since for the algorithm acting on the unrestricted domain, the averaging performed in the construction of intensity maps prefers locations on this skeleton.
r V	r t
H	H
1  1	\ 1 t I ?
film	number of	number of	mean number
	ZnO voxels	spheres	of voxels
			per sphere
57 nm	5,007,344	44,453	112.6436
100 nm	24,238,789	398,833	60.7743
167 nm	28,816,600	480,819	59.9323
where the mean number of voxels covered per single sphere, i.e., the number of ZnO voxels divided by the number of spheres, is considered showing how many ZnO voxels are represented by a single sphere on average.
ZnO voxels
number of spheres
Fig. 11. 2D images of hybrid polymer-ZnO solar cells; left column: smoothed images of 57 nm, 100 nm, and 167 nm films; right column: corresponding representations by unions of spheres.
Empirically, we found out that this restriction of the search domain can improve the runtime of our algorithm essentially. In particular, for large and 'smooth' sets C, the runtime decreases. Whereas for very fine structures as the considered ZnO phases of polymer solar cells, the runtimes are more or less the same.
Some numerical results
The representation of the morphologically smoothed ZnO phase B' by unions of spheres, which is obtained by the stochastic algorithm stated above, is denoted by B". An illustrating example is given in Fig. 11, where a single 2D slice of the macro-scale component B' of a binarized ET image is shown, together with the corresponding slice of its representation B'' by unions of overlapping spheres.
Note that the representation of the geometrically complex set B^ by unions of spheres leads to an enormous data reduction. This can be seen in Table 2.
mean number of voxels per sphere
57 nm 5,007,344 44,453 112.6436 100 nm 24,238,789 398,833	60.7743
167 nm 28,816,600 480,819	59.9323
But an even more important aspect is the fact that the complex structure of the morphologically smoothed ZnO domain is represented by unions of spheres, which can be interpreted as a marked point pattern. The locations of the spheres are the points and the corresponding radii the marks. This representation allows us to analyse and model the set B" using tools from point-process theory; see Stenzel et al. (2011).
Another example of application of the stochastic algorithm developed in the present paper can be found in Thiedmann et al. (2010) where the micro structure of porous materials used in electrodes of lithium-ion batteries is modeled by marked point processes.
Nanostructure on micro-scale
The very fine structural components such as 'thin ZnO branches', 'isolated ZnO particles', and 'polymeric holes', that are removed by the morphological smoothing (see Fig. 9), are important for the physical properties of the films. Hence a multi-scale representation of the ZnO phase is proposed, where each structural component is considered separately. More precisely, we distinguish between the macro-scale component B" of the binarized ET images, which is obtained by the morphological smoothing of the set B and subsequent representation of B' by unions of spheres, and several micro-scale components, which consist of all those voxels that have been misspecified by these image transformations; see also Fig. 9.
In particular, the set B/\B" = {B\B") U {B" \ B) of misspecified voxels is further decomposed into several sub-components, which can be included separately into stochastic simulation models of the ZnO phase as considered in Stenzel etal. (2011). First, two main tvDCS of misspecifications are distinguished:
outer misspecifications and inner misspecifications. Each ZnO voxel belonging to the set B \ B" and therefore constituted as polymer, is called an outer misspecification. Typically, thin branches and isolated ZnO particles are outer misspecifications. On the other hand, each polymer voxel belonging to B" \ B and constituted as ZnO is called an inner misspecification. Inner misspecifications are fiarther subdivided into boundary misspecifications and interior misspecifications. Polymer voxels (belonging to BF) constituted as ZnO and located near the boundary dB" of the macro-scale component B", are called boundary misspecifications, whereas each inner misspecification which is not a boundary misspecification, is called an interior misspecification. Typically, polymeric holes belong to interior misspecifications.
SUMMARY AND CONCLUSIONS
Fig. 12. Flow chart of image processing steps.
In the present paper, a multiscale approach to represent geometrically complex image data in 3D is developed and applied to images of hybrid polymer-ZnO solar cells gained by electron tomography. This is done with the aim to find an efiicient representation of the morphologically smoothed nanostructure by unions of overlapping spheres.
In the first step, the grayscale ET images are smoothed using a Gaussian filter Then the resulting (smoothed) image data is binarized using an adaptive thresholding procedure due to varying brightnesses of the images. A subsequent morphological smoothing of the binarized images provides the basis for the proposed algorithm in order to find an efiicient representation of the morphologically smoothed ZnO phase by unions of overlapping spheres. A flow chart of these image processing steps is given in Fig. 12.
The great advantage of representing geometrically complex structures by unions of spheres is, on one hand, the enormous data reduction, see Table 2, where the mean number of voxels covered per single sphere is shown. On the other hand, it allows the application of tools from point-process theory to fiarther analyze and model the corresponding data sets, because each sphere can be interpreted as a marked point, where the point itself represents the center of the sphere and its mark the corresponding radius.
Since the fianctionality of the considered solar cells heavily depends on the fine structures removed by morphological smoothing, these 'misspecified' voxels are not 'thrown away'. They are interpreted as the 'micro-scale' nanostructure of the binarized ET images, whereas the morphologically smoothed images, represented by unions of overlapping spheres, describe the nanostructure on a 'macro-scale'. Thus, for the description of the nanomorphology of photoactive layers in hybrid polymer-ZnO solar cells, a multiscale approach is proposed, where the macro-scale as well as the micro-scale can be modeled using tools from the theory of random marked point process; see, e.g., Ilhan et af (2008), Kendall et af (2010), Moller and Waagepetersen (2004), and Stoyan et af (1995). This is the subject of a forthcoming paper, where we develop separate stochastic simulation models for the macro-scale and micro-scale of the nanomorphology of photoactive layers in hybrid polymer-ZnO solar cells; see Stenzel etaf (2011).
Another example of application of the algorithm developed in the present paper is considered in Thiedmann et af (2010), where the microstructure of porous materials used in electrodes of lithium-ion batteries is analyzed and modeled on the basis of its representation by unions of overlapping spheres.
These two examplary applications considered in Stenzel et af (2011) and Thiedmann et af (2010), respectively, show that the stochastic algorithm developed in the present paper is a powerful tool for
finding easy to handle representations of geometrically complex 3D structures.
ACKNOWLEDGEMENTS
This research has been supported by Deutsche Forschungsgemeinschaft (DFG) under the Priority Programme: "Elementary Processes of Organic Photovoltaics" (SPP 1355). Part of this work was done while Volker Schmidt was visiting the Isaac Newton Institute for Mathematical Sciences, Cambridge. The authors are grateful to the anonymous referees for their suggestions and comments on the manuscript.
REFERENCES
Angulo J, Jeulin D (2007). Stochastic watershed segmentation. In: Banon GJF, Barrera J, Braga-Neto UM, Hirata NST, eds. Proc 8th Int Symp Math Morphology (ISMM'2007), Rio de Janeiro, Brazil. 265-76.
van Bavel SS, Sourty E, de With G, Loos J (2009). Three-dimensional nanoscale organization of bulk heterojunction polymer solar cells. Nano Letters 9:50713.
Blayvas I, Bruckstein A, Kimmel R (2006). Efficient computation of adaptive threshold surfaces for image binarization. Pattern Recogn 39:89-101.
Brabec C, Scherf U, Dyakonov V, eds. (2008). Organic Photovoltaics: Materials, Device Physics, and Manufacturing Technologies. Weinheim: Wiley-VCH. Cativa Tolosa S, Blacher S, Denis A, Marajofsky A, Pirard JP, Gommes CJ (2009). Two methods of random seed generation to avoid over-segmentation with stochastic watershed: application to nuclear micrographs. JMicrosc 236:79-86.
Faessel M, Jeulin D (2010) Segmentation of 3D microtomographic images of granular materials with the stochastic watershed. JMicrosc 239:17-31. Hoshen J, Kopelman R (1976). Percolation and cluster distribution. I. Cluster multiple labeling technique and critical concentration algorithm. Phys Rev B 14:343845.
Illian J, Penttinen A, StoyanH, Stoyan D (2008). Statistical Analysis and Modeling of Spatial Point Patterns. Chichester: J Wiley & Sons.
Jahne B (2005). Digital Image Processing. 6th ed. Berlin: Springer.
Kendall WS, Molchanov I, eds. (2010). New Perspectives
in Stochastic Geometry. Berlin: Springer. Lantuejoul C (1978). La Squelettisation et son Application aux Mesures Topologiques des Mosaiques Polycristallines. PhD thesis, Ecole des Mines de Paris.
Mßller J, Waagepetersen RP (2004). Statistical Inference and Simulation for Spatial Point Processes. Boca Raton: Chapman & Hall / CRC.
Noyel G, Angulo J, Jeulin D (2007). Random germs and stochastic watershed for unsupervised multispectral image segmentation. In: Apolloni B, Howlett RJ, Jain L, eds. Knowledge-Based Intelligent Information and Engineering Systems. Lect Notes Comput Sci 4694:1724. Berlin: Springer.
Ohser J, Schladitz K (2009). 3D Images of Materials Structures - Processing and Analysis. Weinheim: Wiley-VCH.
Oosterhout SD, Wienk MM, van Bavel SS, Thiedmann R, Koster LJA, Gilot J, Loos J, Schmidt V, Janssen RAJ (2009). The effect of three-dimensional morphology on the efficiency of hybrid polymer solar cells. Nat Mater 8:818-24.
Ridler TW, Calvard S (1978). Picture thresholding using an iterative selection method. IEEE Trans Syst Man Cyb 8:630-2.
Saito T, Toriwaki JI (1994). New algorithms for euclidean distance transformantion of an ^-dimensional digitized picture with applications. Pattern Recogn 27:1551-65.
Serra J (1982). Image Analysis and Mathematical Morphology. London: Academic Press.
Serra J (1988). Image Analysis and Mathematical Morphology, Vol 2: Theoretical Advances. London: Academic Press.
Shaw PE, Ruseckas A, Samuel IDW (2008). Exciton diffusion measurements in poly(3-hexylthiophene). Adv Mater 20:3516-20.
Soille P (2003). Morphological Image Analysis: Principles and Applications. Berlin: Springer.
Stenzel O, Hassfeld H, Thiedmann R, Koster LJA, Oosterhout SD, van Bavel S, Wienk MM, Loos J, Janssen RAJ, Schmidt V (2011). Spatial modeling of the 3D morphology of hybrid polymer-ZnO solar cells, based on electron tomography data. Ann Appl Stat (to appear).
Stoyan D, Kendall WS, Mecke J (1995). Stochastic Geometry and its Applications. 2nd ed. Chichester: J Wiley & Sons.
Thevenaz P, Blu T, Unser M (2000). Interpolation revisited. IEEE Trans Med Imaging 19:739-58.
Thiedmann R, Stenzel O, Spettl A, Shearing PR, Harris SJ, Brandon NP, Schmidt V (2010). Stochastic simulation model for the 3D morphology of composite materials in Li-ion batteries. Research Report, Ulm University.
Unser M (1999). Splines. A perfect fit for signal and image
processing. IEEE Signal Process Mag 16:22-38.	control. Macromolecules 40:1353-62.
Yang X and Loos J (2007). Toward high-performance Yanowitz SD, Bruckstein AM (1989). A new method for polymer solar cells: The importance of morphology	image segmentation. Comput Vision Graph 46:82-95.