Image Anal Stereol 2012;31:27-42 Original Research Paper
NETWORKS OF NANOPARTICLES IN ORGANIC-INORGANIC COMPOSITES: ALGORITHMIC EXTRACTION AND STATISTICAL ANALYSIS
Ralf Thiedmann1, Aaron Spettl1, Ole Stenzel^*1, Thomas Zeibig1, James C. Hindson2, Zineb Saghi3, Neil C. Greenham2, Paul A. Midgley3 and Volker Schmidt1
institute of Stochastics, Ulm University, Helmholtzstrasse 18, 89069 Ulm, Germany; 2Cavendish Laboratory, University of Cambridge, Thomson Av., Cambridge CB3 OHE, UK; department of Materials Science & Metallurgy, University of Cambridge, Pembroke Street, Cambridge CB2 3QZ, UK
e-mail: ralf.thiedmann@uni-ulm.de, aaron.spettl@uni-ulm.de, ole.stenzel@uni-ulm.de ^corresponding author), thomas.zeibig@uni-ulm.de, jch64@cam.ac.uk, zs256@cam.ac.uk, ncgl 1 @cam.ac.uk, pam33@cam.ac.uk, volker. schmidt @ uni-ulm. de
(Received July 1, 2011; revised December 13, 2011; accepted December 18, 2011) ABSTRACT
The rising global demand in energy and the limited resources in fossil fuels require new technologies in renewable energies like solar cells. Silicon solar cells offer a good efficiency but suffer from high production costs. A promising alternative are polymer solar cells, due to potentially low production costs and high flexibility of the panels. In this paper, the nanostructure of organic-inorganic composites is investigated, which can be used as photoactive layers in hybrid-polymer solar cells. These materials consist of a polymeric (OCiCio-PPV) phase with CdSe nanoparticles embedded therein. On the basis of 3D image data with high spatial resolution, gained by electron tomography, an algorithm is developed to automatically extract the CdSe nanoparticles from grayscale images, where we assume them as spheres. The algorithm is based on a modified version of the Hough transform, where a watershed algorithm is used to separate the image data into basins such that each basin contains exactly one nanoparticle. The algorithm is validated by application to (semi-) synthetic data. After the extraction of the nanoparticles, neighboring ones are connected to form a 3D network that is related to the transport of electrons in polymer solar cells. A detailed statistical analysis of the CdSe network morphology is accomplished, which allows deeper insight into the hopping percolation pathways of electrons.
Keywords: 3D watershed, charge transport, electron tomography, Hough transform, network morphology, nanoparticle system.
INTRODUCTION
The rising global demand in energy and the limited resources in fossil fuels require new technologies in renewable energies like solar cells. Classical silicon solar cells offer a good efficiency but suffer from high production costs. A promising alternative are polymer or hybrid-polymer solar cells due to potentially low production costs and high flexibility of the panels. In this paper, we investigate the nanostructure of organic-inorganic composites that can be used as photoactive layers in hybrid-polymer solar cells. More precisely, blends of CdSe nanoparticles (w 6.5 nm diameter spheres) with poly[2-methoxy-5 -(3' ,7' -dimethyloctyloxy)-1,4-phenylene vinylene] (OCiCio-PPV) are investigated. Solar cells based on polymer:nanoparticle blends, in particular CdSe nanoparticles, have been fabricated for a number of years with a great deal of success (Greenham et a/., 1996; Huynl\etal, 2002; Wang et al,, 2006; Dayal et al,, 2010).
It is clear that there is a close relationship between nanostructure and functionality of this kind of composite materials when using them in hybridpolymer solar cells. Under exposure to light, photons are absorbed in the polymer phase and neutral excited states known as excitons are generated, which diffuse in the polymer phase within a limited life time. If an exciton reaches the surface of a CdSe nanoparticle, it is split up into an electron (negative charge) in the CdSe and a hole (positive charge) in the polymer phase. In addition, CdSe nanoparticles can also absorb photons directly, followed by hole transfer to the polymer. The electrons have to hop from CdSe nanoparticle to nanoparticle to reach the electrode. Since an electron can only hop over short distances between nanoparticles, the morphology of the network of nanoparticles as a whole is of great importance for the performance of the solar cell.
Therefore, it is essential to develop tools in order to extract the ensemble of CdSe nanoparticles
from 3D grayscale images obtained by electron tomography. Connecting neighboring CdSe particles yields a network through which the electrons can hop to reach the electrode. A detailed statistical analysis of the network morphology is accomplished, which allows deeper insight into the hopping percolation pathways of electrons.
The algorithm to extract the ensemble of CdSe nanoparticles is applied to grayscale image data gained by high-angle annular dark-field scanning transmission electron microscopy (HAADF-STEM). Examples of image data for this type of composite material are displayed in Fig. 1, where the upper part shows a thin 2D section and the lower part a 3D cutout of a suitable binarization obtained by global thresholding. In the grayscale image on the upper part of Fig. 1 the gray values are inverted, i.e., the nanoparticles are displayed dark (low gray values) and the background is bright (high gray values).
Fig. 1. CdSe nanoparticles embedded in a OC1C10-PPV phase: Thin 2D slice (top; the scale bar corresponds to 50 nm) and 3D cutout ofbinarized data (bottom; obtained by global thresholding).
Note that binarization and subsequent manual extraction of nanoparticles in 3D grayscale images is accompanied by an enormous effort. Furthermore, small failures in the binarization procedure can lead to quite different results. Hence, in the present paper, a computer algorithm is developed for automatic (i.e. non-interactive) extraction of spherical objects directly from 3D grayscale images. The precise extraction of the CdSe nanoparticles from electron tomography images and the statistical analysis of the morphology of CdSe networks are of significant importance to characterize transport properties of the considered composite material with respect to electron hopping.
In this paper we focus on the case of spherical nanoparticles, so-called nanodots. However, the basic idea of our extraction algorithm also works for nanoparticles with other shapes. For example, for elongated nanoparticles, so-called nanorods, just the maximization step in the Hough transform for the detection of nanorods will be computationally a bit more complex.
The paper is organized as follows. At first, the imaging technique and the resulting 3D data are described, which are the basis of our investigations. Then, a robust algorithm of image processing is developed for fully automatic extraction of nanoparticles from 3D grayscale images. Subsequently, the morphology of the system of extracted spheres representing the CdSe nanoparticles is analyzed with respect to properties that are relevant e.g. to quenching of excitons. Consecutively, network graphs based on the detected spheres are constructed, whose structural properties are closely related to charge transport. The algorithm to extract nanoparticles is validated by detecting nanoparticles from (semi-) synthetic 3D image data, where the true configuration of the nanoparticles is known. Finally, the results are summarized and an outlook to possible future research topics is given.
3D IMAGES OF CDSE-(OC1C10-PPV) COMPOSITES
In this section, the considered organic-inorganic composite is briefly described. Moreover, the imaging technique applied to generate 3D images of this material with high spatial resolution is explained. For further details on materials and imaging, we refer to Hindson et al. (2011).
DESCRIPTION OF THE MATERIAL
The photoactive layer of the considered hybridpolymer solar cell consists of a blend of CdSe nanoparticles (w 6.5 nm diameter spheres with butyl-amine ligands) with poly[2-methoxy-5-(3',7'-dimethyloctyloxy)-l,4-phenylene vinylene] (OC1C10-PPV). The ratio of CdSe nanoparticles to OC1C10-PPV was 6:1 by weight. The photovoltaic device was built using regular architecture. Cleaned substrates were coated with PEDOT:PSS which was annealed under nitrogen for 30 minutes at 230°C. The active layer was spin-coated at 2000 r.p.m. which yields a thickness of 50-70 nm, and annealed at 150°C for 20 minutes. Films for tomography were floated off the substrate onto water, using the PEDOT:PSS as a sacrificial layer, and were then transferred onto TEM grids.
IMAGING TECHNIQUE AND TOMOGRAPHIC RECONSTRUCTION
To investigate the nanomorphology of the ensemble of CdSe nanoparticles embedded in the polymeric phase, we used 3D images obtained by high-angle annular dark-field scanning transmission electron microscopy (HAADF-STEM). In general, electron tomography enables to obtain three dimensional images with a nanometer resolution. Thus, it is an ideal method to analyze the structure of nanoparticle networks in organic photovoltaics. The HAADF-STEM has an advantage in systems where the two components differ significantly in atomic number, Z, since it mainly detects electrons that have undergone Rutherford scattering. The 3D grayscale images were obtained using a FEI Tecnai F20 field emission gun TEM operated at 200 kV. The tilt range was between 120 and 142 degrees in 1 or 2 degree steps, depending on the exact topography of the selected area. The image series was aligned and reconstructed using a SIRT algorithm in Inspect3D software.
DESCRIPTION OF IMAGE DATA AND BASIC NOTATION
The complete 3D data set analyzed in this paper has a size of 917 x 872 x 404 voxels, where each voxel represents (0.4 nm)3. An example of a 2D slice can be seen in the upper part of Fig. 1.
W ' A. i • ^ - / J «Fig. 2. Cutout of a 2D thin section with densely packed particles; the scale bar corresponds to 50 nm.
For the statistical analysis described later on, we considered a cutout of 421 x 421 x 91 voxels, which is a part of the complete data set. Beforehand, in order to obtain a large (cuboid-shaped) cutout with as much nanoparticles as possible, the 3D image has been rotated. The considered cutout contains only that part of the whole data set where the CdSe nanoparticles are rather densely packed (Fig. 2), corresponding to the high-density regions described in Hindson et al, (2011).
The values of all parameters are given in voxels, where we regard the darker phase, i.e., the particles, as the foreground phase, and the brighter phase as the background phase (Fig. 2). Note that for the 8bit grayscale images considered, a gray value of 0 indicates black and a gray value of 255 white.
The algorithm developed in order to automatically extract the CdSe nanoparticles from grayscale images is based on mathematical methods. Therefore, it is unavoidable to introduce a certain minimum of notation.
Let I = {I(x.y.z). (x.y.z) € I)} denote the 3D image, where I(x.y. z) € {()..... 255} is the gray value of the 3D image I at grid point (x. y. z) € D c Z3 and the set D is interpreted as a sampling window. Furthermore, in order to simplify the handling of boundary effects, we consider the following (formal) continuation of the image outside of D. For locations (x.y. z) € Z3 \D not in the domain I) of the image I, the gray values I(x,y,z) are assumed to be such that they do not affect the applied operations of image processing (which will be I(x,y,z) = 0 for (x. y. z) € Z3 \ D in most cases).
ALGORITHM FOR NON-INTERACTIVE EXTRACTION OF SPHERES
In this section, we describe an algorithm to extract single particles from grayscale images like
those shown in Fig. 2, where we assume that each nanoparticle can be seen as a 3D sphere. This assumption is supported by the production technique of the particles. The segmentation algorithm developed is based on the Hough transform (HT) for spherical objects (Duda and Hart, 1972; Ballard, 1981; Burger and Burge, 2008). The HT is a robust algorithm to detect parameterized geometric objects, but suffers from long computer runtimes when the number of parameters which describe the objects is high. However, since spheres can be described quite simply by four parameters, the HT is adequate for extracting such objects. Note that this HT-based approach is similar to the idea of maximum-likelihood estimation in mathematical statistics, where the parameter vector (here: position and radius of the CdSe nano-particle) is chosen such that the identified object is the most plausible {i.e., most likely) one given the 3D image data. Hence, the HT is mainly conceived for extracting one single but not several (densely packed) objects from an image. It can be used to extract an ensemble of objects, but either the number of objects has to be known in advance or has to be determined by extensive computer simulations. In this work, we go one step further and perform a data preprocessing, where we divide the sampling window D into disjoint regions each containing exactly one object to be detected. Subsequently, each region is transformed using the HT-principle in order to identify the object located in this region. This procedure reduces computer runtime of the HT even further and is much more robust. The splitting of the domain D into separated regions (so-called 'basins') is based on an iterative thresholding procedure detecting the initial points of the subsequent 3D watershed segmentation. These initial points are called the 'springs' of the 'basins' in the following (Beucher and Lantuejoul, 1979; Beucher and Meyer, 1993; Baere and Lehmann, 2006).
SMOOTHING OF GRAYSCALE IMAGES
A first step of image pre-processing, which usually precedes statistical image analysis, is a smoothing of the data to reduce noise, i.e., small artifacts resulting from the imaging and reconstruction techniques. Thus, we apply a 3D symmetric Gaussian filter G of size 3x3x3 and standard deviation o = 0.4 to the image I. Note that such a small filter has only negligible influence on the geometric structure of the data we are interested in.
The application of the filter can be seen as a convolution of the image I with the (matrix) filter G = (G(i,j,k),i,j,k € {—1,0,1}). This means that the gray values I'(x,y,z) of the smoothed image I' = {I'{x,y,z), (x.y.z) € D}, where the filter G of size
3x3x3 with centered origin has been applied to, are given by
I'(x,y,z) =
l l l
L L L I(x + i,y + j,z + k)-G(i,j,k) .
i=-1 j=-\k=-1
A 2D slice of the smoothed 3D image, together with the corresponding 2D slice of the original image, can be seen in Fig. 3.
Fig. 3. 2D sections of the original (top) and smoothed (bottom) image.
MORPHOLOGICAL CLOSING
To emphasize and additionally smooth the spherical structures that have to be extracted from the given image, we apply a morphological closing with structuring element of size 5x5x5 and centered at the origin (Serra, 1982). More precisely, the gray values I"{x,y,z) of the dilated image I" = {I"{x,y,z), {x,y, z) € D} are given by
I" (x, v, z) = max I' (x + i,y + j,z + k).
—2<i,j,k<2
The size 5 x 5 x 5 is chosen because it can be warranted that no nanoparticles are lost. Note that for applying the gray level dilation, the values of the image I' outside its domain D are put equal to 0. Analogously, the gray level erosion applied to the image I" is given by
I"'(x,y,z)= min l"(x + i,y+j,z+k),(x,y,z)£D,
where the values of the image I" outside its domain D are put equal to 255. The result of these morphological operators is an enhancement of the spherical shape of the CdSe nanoparticles that shall be detected. Their edges are emphasized and, in addition, their interior is smoothed which improves the edge detection considered below.
Note that the morphological closing has no significant influence on the extracted objects with respect to changes of their sizes. Examples of thin planar sections of the data after application of morphological closing are shown in Fig. 4.
Fig. 4. 2D sections of the original data (top), smoothed image I" after gray level dilation (middle), smoothed image I'" after morphological closing (bottom).
GRADIENT IMAGE
The application of HT in order to extract the nanoparticles first requires the detection of edges of the
objects of interest. One possibility to detect edges in the image I'", i.e., the locations of rapid changes in the gray values, is the consideration of its gradient image. Changes in gray values are mathematically described by partial derivatives, where the locations	are
supposed to form a 'continuous' set. Then, the vector Vl'"(x,y,z) with
VIm(x,y,z) =
id \
I l^.v,;) J
containing the partial derivatives is considered, i.e., VI"'(x,y,z) describes the changes of the gray values of I'" at location (a\ v. z) in the directions of the a*-, y-, and z-axis, respectively. The length \Vf'"(x.y.z)\ of the vector VI"'(x,y,z) describes the (total) strength of changes of the gray values at position (x. y. z), i.e., the longer the vector the stronger the changes in gray values. Since edges correspond to strong changes in gray values, this provides a possibility of edge detection.
Fig. 5. 2D section of gradient magnitude image of I'" (brightened up for demonstration purposes).
A challenge in practical application is the numerical evaluation of the partial derivatives since only a discrete domain D is given. Several solutions to overcome this classical problem have been proposed in literature. We used an algorithm implemented in Avizo, a software of the Visualization Sciences Group Cwww.vsg3d.com), to obtain an approximation of the gradient values. A 2D slice of the 3D gradient image G is shown in Fig. 5, where the magnitudes \Vf'"(x.y.z)\ are plotted as grayscale values, i.e., areas of high variation are displayed as bright regions. Later on in this section we use this gradient image G when considering the Hough transform.
ITERATIVE THRESHOLDING
To determine the initial points of the watershed transformation considered later on in this section, i.e. the 'springs' of the 'basins', we propose the following procedure of an iterative thresholding. The appropriate choice of these initial points is crucial for the subsequent detection of spherical objects since each basin shall contain exactly one object.
The idea of our iterative thresholding is to detect the centers of the objects to be extracted, which are assumed to be the locally darkest parts of the image I'". In particular, we are looking for a set S of initial points for the watershed transformation. To begin with, we put S = 0, where 0 denotes the empty set. Then, we consider the (global) threshold t = 1, which is applied to the image I'", and perform the following steps:
1)	The image I'" is thresholded with threshold value t to get the binarized image I't", i.e.,
T'"(xvA-{ 1 if '"'(-W)^, h 0 if 7///(a. _V)Z)>f)
and the phase of I'", where I"' equals 1, is called foreground. Note that the foreground is growing with increasing value of t, i.e., {foreground of l't"} C {foreground of I',", }.
2)	All connected clusters of the foreground in l't", i.e., the parts of I't" with value 1, are determined. This can be done using e.g. the algorithm of Hoshen and Kopelman (1976). Let {C\.. ■■■C'k(lj} denote the set of separated (foreground) clusters in I"' with volumes \C\\,...,	respectively. If the volume of a cluster, say cluster C\, exceeds a predetermined size, in our case 100 voxels, we check whether this cluster has already contributed to an initial point for the watershed transformation or not. Therefore, we check if one of the points from S = {si,....srii ,} is contained in C\. If no point of the current set S is contained in Cj, i.e., Sj ^ C\ holds for all j <G {1..... /7, i}, then the barycenter of C\ is considered as potential 'spring' of a basin to be constructed. If, in addition, the bary center of C\ is contained in C\, the bary center is added to the set S of initial points. Note that only clusters smaller than 2000 voxels are allowed to generate an initial point. These conditions assure that only (at least roughly) spherically shaped objects are detected by the object-detection algorithm considered later on in this section.
3)	The current threshold value t is increased by one and the procedure described above is performed again, starting with step 1. The procedure is stopped at that value of t if, for the first time, the
complete binary image is said to be foreground, i.e., I't"(x,y,z) = 1 for all (x,y,z) £ D.
The result is a set of initial points 5 = {si,..., sn}, which is used for the watershed transformation. Due to the construction principle of the set S, each point of S is located close to the center of a spherical object that shall be detected later on. Note that the 'critical' size of the cluster, i.e., the minimum size that the cluster is considered to generate a point in S, controls the number of detected objects.
Fig. 6. Schematic ID example of watershed transformation with 3 basins.
WATERSHED TRANSFORMATION
In the previous section, we determined the set S = {si,... ,s„} of initial points in order to construct the watershed transformation of I'". Note that the watershed transformation subdivides the domain D of I'" into pairwise disjoint sets, so-called basins B\,... ,Bn, i.e., BiD Bj = 0 for any i / j and U" i Bt = D, where Si <G Bi for each i <G {1,..., n}.
The intuitive idea of the watershed transformation is that at the initial points si,..., sn, the flooding of I'" starts. Therefore, the gray values of I'" are interpreted as level curves. The boundaries of the basins B\,...,Bn are the locations where different 'valleys' meet, see the schematic illustration displayed in Fig. 6 for the one-dimensional case. More information about the watershed transformation and related algorithms can be found e.g. in Beucher and Lantuejoul (1979); Beucher and Meyer (1993); Roerdink and Meijster (2000); Baere and Lehmann (2006).
In the following sections, we show that the watershed transformation described above leads to a decomposition of the domain D into pairwise disjoint basins B\,...,Bn, where each basin contains exactly one object {i.e. CdSe nanoparticle) to be detected. A 2D slice of the watershed barriers obtained in this way is displayed in Fig. 7. Note that these basins are used to decompose the domain of the gradient image G, considered previously in this section, into disjoint (sub-) regions. They are considered for the detection of spheres using a modified version of the well-known Hough transform. More precisely, for the application
of the HT, the number of spheres to be detected needs to be known in advance. By the watershed transformation, not only the number of spheres is known, but also the area where the HT is looking for every single sphere can be reduced.
Fig. 7. 2D section of watershed barriers.
The watershed transformation detects basins in the whole image, independently of the local contrast. Because each basin is considered as a candidate for the detection of spheres, it is reasonable to ignore basins where the grayscale values do not differ too much. The criterion used to quantify the local contrast in basins is as follows. We compute the standard deviation of the grayscale values si =
(l(W)€5ia'"(^v,z)-m/)2/(|fii|-l))1/2 for each basin Bi, where nii = Yi{x,y,z)eBiI'"(x^y^)/\Bi\ denotes the sample mean of grayscale values in the same basin and \Bi\ the number of voxels of basin B,. If the standard deviation is smaller than a given threshold, then the basin is ignored in the following step. For the current image data a (minimum) standard deviation of 15 has been determined as a suitable value. The remaining basins are denoted by B[,...,B'n,, where n' < n.
HOUGH TRANSFORM FOR THE DETECTION OF SPHERES
The decomposition of the domain D into disjoint basins B\,.... B'n,, where each basin contains exactly one object to be detected, provides optimal preconditions for the Hough transform (HT) (Ballard, 1981; lahne, 2005; Burger and Burge, 2008). The HT is a powerful tool to detect (single) geometric objects from grayscale images provided that the objects to be detected possess a parametric description. In case of a three-dimensional sphere, there are four parameters describing it uniquely: the coordinates x,v,z of the center and the radius r.
The idea of the HT is the consideration of the so-called Hough space, which, in our case, has four
dimensions, representing the parameters a\ _v, z, and r. Note that for the application to image data, only a discrete set of values for a\ _v, z, and r is considered. In its original version, the HT is applied to binary images only, which leads to the following interpretation of the Hough space. Each vector {x,y,z, r) represents one sphere in the Hough space, where the value at (x, v. z. r) in the Hough space provides the number of boundary voxels of the foreground phase of the original binary image contributing to the surface of this sphere, i.e., the number of voxels on the surface/edge of the (discretized) sphere in the original image. The sphere that we are interested in can then be determined just by looking for the global maximum in the Hough space. This procedure is plausible since a sphere is chosen in such a way that the largest number of detected edge points is covered by that sphere.
However, in our case, edge detection is based on the gradient image, which is a grayscale image. In view of this and since any binarization yields inaccuracies, we used a slightly modified version of the HT that can be applied to grayscale images. Instead of counting the number of edge points covered by the surface of a sphere, we sum up their grayscales in the gradient image. This leads to the detection of an object, where binarization artefacts can be avoided. Furthermore, for all voxels outside the currently processed basin, a gradient of zero can be assumed. Therefore structures outside the basin cannot affect the sphere-detection.
It is obvious that using this version of HT, larger spheres are preferred since more edge points can be covered. To reduce this overestimation of large spheres, we applied a rescaling of the values in the Hough space, dividing each value by the square root sjr of the corresponding radius r, before searching for the global maximum. This kind of rescaling has been chosen because then the result of the sphere extraction provides a good optical fit. However, other scaling factors like r 1 or 1/logr are also possible. The rescaling does not affect the detected radii directly, but it allows to prefer smaller / larger radii which are plausible within a certain range. Later on, the extraction algorithm is validated indicating that the rescaling is adequate.
Note that this HT-algorithm is applied to each basin B\,..., B'n, separately, which leads to the detection of n' spherical objects in the whole image. An illustration is given in Fig. 8, where 2D slices of detected spheres have been added to the data.
Fig. 8. Detection of spheres using the Hough transform; arrow indicates Mo strongly overlapping spheres.
POST-PROCESSING OF DETECTED SPHERES
Some of the extracted spheres may strongly overlap; see Fig. 8 for an example of such (almost completely) overlapping spheres indicated by an arrow. A reason for this overlapping is that there are too many basins identified by the iterative thresholding procedure which has been described previously in this section. Also note that the spheres identified by the HT may exceed the boundary of the basin. Due to the varying grayscale values at the centers of the spheres, it happens that the iterative thresholding procedure detects two different initial points which actually belong to one and the same particle. These two initial points then lead to two different basins in the watershed algorithm and, accordingly, two spheres are determined by the Hough transform although it should be only one. From a physical point of view, it is clear that such a strong overlapping cannot be a correct detection of spheres. Hence, we apply the following post-processing of detected spheres to reduce the overlapping.
Fig. 9. Detected spheres after post-processing.
To avoid an over segmentation, we are looking for all pairs of spheres where the volume of overlapping is larger than 5% with respect to the volume of the smaller sphere. Subsequently such a pair of spheres is removed from the list of detected spheres and their basins obtained by the watershed segmentation are recombined. From this newly created basin, a new sphere is extracted and added to the list of detected spheres, where the same techniques are used as described above. This procedure is repeated several times until no pair of spheres is found with an overlapping of 5% or more of the volume of the smaller sphere. In the result, a clearly visible reduction of overlapping effects can be observed (Fig. 9 and Table 1). Note that a small overlapping of the spheres seems reasonable since in reality the nanoparticles have no perfect spherical shape.
Table 1 .Analysis of overlapping spheres (after the first run of sphere-extraction algorithm).
post-processing before after		
total volume		
fraction	0.290	0.241
of spheres		
volume fraction		
covered	0.012	0.001
several times		
The volume fraction of extracted spheres decreases considerably by the post-processing described above (Table 1). However, after the repeated extraction of spheres discussed in the following section, the total volume fraction of all extracted spheres is equal to 0.311.
REPRESENTATION AS MARKED POINT PATTERN
Let B denote the inorganic part (i.e., the CdSe nanoparticles) of the composite material. In the preceding sections we showed that the set B can quite nicely be represented by a system of spheres which can be seen as functional of a marked point pattern {si,n),... ,(sm,rm), in < n', where the points si are the centers of the spheres and the marks r, the corresponding radii. In other words, the set B can be approximated by the union of these spheres, i.e., B & U-Ii b(siJi), where b(s, r)={i'eR3: |s' - s \ < r}.
REPEATED EXTRACTION OF SPHERES
Since, in general, the result of sphere extraction described above does not cover all regions that can be interpreted as nano-sized particles, we applied the sphere-detection algorithm again. The idea is to first remove already detected spheres from the image and then search for spheres in the remaining parts.
% ' v i
- ; jp ^ i
. /
Fig. 10. Removal of detected spheres from I"' (top) and gradient image G (bottom; brightened for demonstration purposes).
The previously detected spheres should not influence the current search for spheres. Therefore, they are removed from the image and the same search algorithm as described above is applied again. To remove the spheres (si,ri),... , (sm,rm) detected in the first run of the algorithm we apply a dilation of the spheres, which can be seen as a 'blowing up' of the union of spheres IJi" i b(si. ''/ )■ More formally, using Minkowski addition, a dilated union of spheres is obtained, which is given by SD = {s + s' : s £ b(o,4), s' € U"ii b(si,ri)} (Serra, 1982). The set SD is then compared to the filtered image I'" and the gradient image G, described in the previous section, where all voxels covered by SD are set to white or black, respectively, i.e., they are set to background (Fig. 10). Note that in this thin 2D slice, due to the dilation of the spheres there are also parts removed from spheres located in neighboring slices.
Fig. 11. Detected spheres of the second run (top), final result of sphere extraction (bottom).
With these modified images, we performed the same procedure of sphere-detection as described in the previous section. An example of the resulting set of spheres is shown in Fig. 11 (top).
The spheres detected in the second run are assumed to be lower-ranked. They are only added to the previously found set of spheres (si, n),..., (sm,rm) if they do not overlap too much wherefore the same criterion (5% overlapping volume) is applied as described in the previous section.
The final set of spheres representing the CdSe nanoparticles is considered as marked point pattern and denoted by (si,r\),...,{sk,rk), where k > in (Fig. 11 (bottom), and Fig. 12). In the following, besides the marked point pattern (s\,ri),... ,(sk,n) itself, we will analyze the union of spheres \Ji=\b{si,ri) representing the system of CdSe nanoparticles with respect to electrical conductivity, where we assume that this system can be described by a 3D graph. The volume fraction of the union \Jl=\b(si,ri) of spheres after repeated extraction is equal to 0.311.
Fig. 12. Final result of sphere extraction in the same slice as in Fig. 1 (top), 3D cutout of detected spheres (bottom).
STATISTICAL ANALYSIS OF EXTRACTED SPHERES
Like in the preceding sections of this paper, for the statistical analysis provided in this and the following section we restricted ourselves to a densely packed part of nanoparticles of the whole data set which can be assumed to be statistically homogeneous, a required condition on the data for the applied analysis tools.
Several structural characteristics of the marked point pattern (si,ri),..., (Sk,rk) of extracted spheres representing the CdSe nanoparticles have been considered. To begin with, we analyze the distribution of particle sizes. Then, we consider some structural characteristics of particle locations, where we analyze the pair correlation function and the nearest-neighbor distance distribution function of sphere centers. We also computed the cumulative distribution function of spherical contact distances to the union of spheres Uf=i b(si,rj). Finally, the subsequent section
deals with three-dimensional graphs formed by these spheres.
PARTICLE SIZES
The sphere representation of the CdSe nanoparticles derived in the previous sections enables us to determine the distribution of particle sizes by considering the distribution of radii r\,..., rk of the extracted spheres. The histogram of these radii is shown in Fig. 13. Note that the variability of radii is rather small. More than 90% of the analyzed spheres possess radii between 5 and 10 voxels, where the mean value is equal to 7.82 voxels. This corresponds to an average sphere diameter of 6.23 nm. This is in reasonable agreement with the mean particle diameter determined experimentally for nanoparticles dispersed on a surface in the absence of any polymer (Hindson et al, 2011).
0.6
0.5
o
0 0.4
3
cr
s 0.3
0
1	0.2
0
0.1
0
0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 radius in voxels
Fig. 13. Histogram of radii.
PARTICLE LOCATIONS
In order to analyze structural properties of the point pattern s\,...,sk of particle locations (i.e., centers of detected spheres) we consider the pair correlation function and the nearest-neighbor distance distribution of s\,...,sk. These characteristics are popular in statistical analysis of geometrically complex point patterns (Diggle, 2003; Illian etal, 2008; Gelfand et al, 2010). Furthermore, we determine the cumulative distribution function of the (voxel-wise) spherical contact distances to the union of spheres |J*=i b(si,rj).
The first diagram in Fig. 14 shows the pair correlation function g{r),r > 0, of sphere centers si,..., Sk, where g{r) is the relative frequency of pairs of sphere centers with distance r from each other. Note that in the case of complete spatial randomness (CSR, i.e. homogeneous Poisson point process), that is, there is no interaction between the points of the considered point pattern, it holds that g{r) = 1. Furthermore, g{r) > 1 indicates clustering of points, whereas g{r) < 1 means repulsion of points. In our case, a repulsion of spheres can be observed, which is about up to
distances of 14 voxels (with a minimum distance of r = 8 between pairs of sphere centers), and, for distances between 16 and 23 voxels, there seems to be a moderate clustering of sphere centers. However, in view of the relative small variability of radii mentioned earlier in this section it is quite clear that the function g(r) shown in the first diagram of Fig. 14, cannot be the pair correlation function of a so-called 'dense packing' considered in mathematical physics. In this case, there should be much larger deviations of the values of g(r) from level 1, having in mind the small variability of radii; see also Fig. 6.11 of Illian et a I. (2008).
Similar structural properties can be concluded from the second diagram in Fig. 14 which shows the cumulative distribution function D(r), r > 0, of nearest-neighbor distances between pairs of sphere centers. In case of CSR it holds that D(r) = 1 — exp(— Xk^1), where A is the mean number of sphere centers per unit volume, d the dimension and Kd the volume of the rf-dimensional unit sphere. Thus, in this case, D(r) is a strictly concave function, whereas the diagram shown in the middle of Fig. 14 first runs along the j-axis, then having a convex part, and being concave only later on (for r > 15). Notice that this type of diagram clearly indicates repulsion of sphere centers for small distances with the same minimum distance r = 8 as indicated by the pair correlation function. We also remark that the distances between pairs of spheres (i.e. CdSe nanoparticles) are of great importance for the transport of electrons by the system of CdSe particles. However, note that the distribution function D(r) does not (yet) take into account the size of the extracted particles.
The third diagram in Fig. 14 shows the cumulative distribution function H(r),r > 0, of the minimum distance from an arbitrary location in the polymer phase, chosen at random, to the union of spheres \Ji=ib(si,ri). Note that H(r) can be determined by spherical dilations of the union of spheres Uf=i b{si, r,) (Serra, 1982). This characteristic also provides important information regarding the performance of polymer solar cells since it is closely related to quenching probabilities of excitons (Oosterhout etal, 2009). Suppose that tq is the (average) diffusion length of excitons. Then, H(tq) is the fraction of those voxels classified as polymer, whose minimum distance to the set Uf=i b(siji), i.e. the CdSe phase, does not exceed the expected diffusion length of excitons. This means that the larger the values H(r) for r <G (0,ro), i.e., the more voxels classified as polymer are closer to the union of CdSe nanoparticles than the expected diffusion length of excitons, the more charges can be generated. As expected for these aggregated regions containing high densities of nanoparticles, typical
distances required for exciton diffusion to an interface are small (less than 2 nm), so exciton dissociation is expected to be efficient (Hindson et al, 2011).
1.6
C 1.4
0
'o 1.2
C
•3 1
C
1	0-8
2	0.6
o
^ 0.4 nJ
^ 0.2 0
f 0.8 Xi n5 X!
O 0.6
CL
Q)
0
0 5 10 15 20 25 distance in voxels
1
f 0.8 X!
o 0.6
0
0 2 4 6 8 10 12 distance in voxels
Fig. 14. Pair correlation function of sphere centers si,...,Sk (top), cumulative distribution function of nearest-neighbor distances between sphere centers (middle) and spherical contact distances to the union of spheres Uf=1	(bottom).
3D GRAPHS OF NANOPARTICLES
In this section we briefly describe how 3D graphs can be defined which connect the detected spheres (si,r\),...,according to given rules. Furthermore, connectivity properties of these graphs are investigated, which are important for the electrical conductivity.
From a qualitative inspection of the images and from the measured device performance (Hindson et a I, 2011) we would expect electron transport to be relatively efficient through these aggregated regions containing high densities of nanoparticles. The analysis below suggests that
distance In voxels
hopping between particles over distances larger than 1 nni is required for this to be the case. It would be interesting in the future to perform sphere extraction and transport analysis on regions of the film with lower nanoparticle densities where the electron transport is expected to be less efficient.
other if the distance dtJ = min{ |.v — y | : x € bis, ./',)■>■ € b(sj,rj)} = max{|sj — sj\ — r, — rj,0} between the corresponding spheres is smaller than some predefined threshold dmax > 0.
For our analysis, we use different values for dmax,
i.e., di!Jx = 2, dmax = 4, and dm}x = 6. For each of these values of dmax, the result is a 3D graph ('V,E), where the set of vertices V is equal to the set (si,...,Sk) of sphere centers. The set of edges E consists of the (undirected) line segments between pairs of points from V such that the distance between the corresponding spheres is smaller than dmax as explained above. Cutouts of such 3D graphs obtained in this way are shown in Fig. 15.
5" 0.2
0 2 4 6 8 10 12 14 16 coordination number
0.3
'Jlk.. .'
0 2 4 6 8 10 12 14 16 coordination number
Fig. 15. Cutouts of 3D graphs with hopping distance dmax = 2 (top), d,,^ = 4 (middle), and d,„ax = 6 (bottom).
MODEL DESCRIPTION
Transport of electrons within the CdSe-(OCiCio-PPV) composite depends on its interior connectivity, i.e., on the existence of percolation pathways towards the electrodes via a network of nanoparticles whose maximal neighboring distance is below the hopping distance of electrons. Since the electrons can hop between the nanoparticles only within a certain maximum distance dmax > 0, we analyze the detected spheres (si, n),..., (s*-, ?>) with respect to their connectivity using a graph representation for the ensemble of these spheres. This means that those pairs Si, sj, i / j, of sphere centers are connected to each
0.3
— o 1
Q)	1
0
0 2 4 6 8 10 12 14 16 coordination number
Fig. 16. Histograms of the coordination number for dmax = 2 (top), dmax = 4 (middle), and dinax = 6 (bottom).
COORDINATION NUMBER
A common characteristic to describe the connectivity of a graph is the distribution of degrees of its vertices, i.e., the number of edges emanating
from each vertex. By some authors, this characteristic is called the coordination number of the graph. The histograms of the coordination numbers, which have been obtained for the three graphs of nanoparticles introduced in earlier in this section, are displayed in Fig. 16. It is clearly visible that the shape of the histograms essentially depends on the choice of the hopping distance dmax. For dmax = 2 the histogram is quite narrow, whereas for dmax = 4 and dmax = 6 it becomes wider and more symmetric. In particular, the mean coordination number increases from 2.2 for dmax = 2 to 4.7 and 7.0 for dmax = 4 and dmax = 6, respectively.
MINIMUM-SPANNING TREE
Another characteristic describing the connectivity of graphs is the so-called minimum spanning tree (MST), which is a popular tool in graph theory (lungnickel, 1999; Diestel, 2005). The idea is to consider the sub-graph with the minimum length but the same connectivity as the original graph, i.e., all nodes that can be connected in the original graph by a sequence of edges are also connected in the MST. Then, as a characteristic to describe connectivity, the relative length £ of the MST is considered, i.e.,
length of the MST length of the original graph '
where we obtain the values £ = 0.83, £ = 0.48, and £ = 0.32 for the graphs with hopping distance dmax = 2, dmax = 4, and dmax = 6, respectively. These values can be interpreted as follows. The graph with an assumed hopping distance of dmax = 2 has no good connectivity at all, i.e., it is hard to find any good percolation pathways. More precisely, for dmax = 2 the graph consists of 119 isolated subgraphs. In contrast, for dmax = 4, the graph consists of 9 subgraphs, where for dmax = 6, the whole graph is connected. Thus, the graph with hopping distance dmax = 6 has very good connectivity properties, i.e., the charge transport to the electrodes should work very well in this case.
GEOMETRIC TORTUOSITY
For describing transport processes in composite materials, the tortuosity of their phases is an important characteristic. It is usually defined as the ratio of the mean effective path length through a material divided by the material thickness. In this paper, we use a geometric approach to describe this kind of property. We consider shortest path lengths instead of effective path lengths. This has the essential advantage that not only a single value is obtained, like the effective path length being a mean value, but a whole
distribution of local geometric tortuosities can be considered, which contains much more information. Additionally, the shortest-path approach can be seen as a purely structural method, i.e., it does not dependent on physical constraints (Decker et al, 1998; Peyrega et al, 2009; Thiedmann et al, 2009).
For electrical conductivity within CdSe-(OCiCio-PPV) composites the Euclidean distances, which electrons have to go, are not important, but rather the number of hops of charges until reaching the electrode. Thus, shortest paths are considered in the way that the number of hops is counted. This approach is quite close to the physical understanding, where a large number of hops causes increased likelihood of losses by recombination. The normalization is done with respect to the material thickness, which is supposed to be the Euclidean distance between the left-hand and right-hand side of the boundary of the sampling window. Thus, the number of hops relative to the material thickness is considered. This quantity is related to the required energy of an electron to move one voxel size closer towards the electrode. Note that in contrast to the standard definition of tortuosity, values smaller than 1 are possible.
The results for maximum hopping distance of d,„ax = 4 and dmax = 6 are displayed in Fig. 17. Note that both histograms are quite narrow. Anyhow, it can be seen that they (slightly) depend on the selected maximum hopping distance dmax.
100
.1 40 20 0
0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 hops per voxel
100
.1 40
il Q)
20 0
0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 hops per voxel
Fig. 17. Histograms of geometric tortuosity for dmax = 4 (top) and dmax = 6 (bottom).
L
Because of the bad connectivity of the graph with maximum hopping distance of dmax = 2, for most starting points there does not exist a path to the end points. Thus, in this case, it is not reasonable to compute the histogram of tortuosity.
the their radii (Fig. 20). In particular, the mean distance between the centers of two matched spheres is equal to 0.6 voxels and the mean of the (absolute) difference between the radii is equal to 0.13 voxels. In comparison to the mean radius of 7.82 voxels these differences are very small.
VALIDATION
To validate our sphere-extraction algorithm, we use the extracted spheres bis,.r(), i 1 ...k, to draw them into an empty 3D image. After adding a Gaussian blur and white noise (Burger and Burge, 2008) we apply the algorithm once again. This approach allows us to compare the detected spheres with the (known) spheres given in the input image.
GENERATION OF SYNTHETIC INPUT IMAGE
We use the centers and radii of the spheres detected from the original image to create a new 3D image with the same dimensions. Similar as in real data, the spheres are drawn in a light gray (gray value 175), the background is a dark gray (gray value 100). To make the extraction of the spheres more difficult and especially more realistic, we apply a Gaussian blur with standard deviation 2.5 and insert white noise subsequently. White noise is modeled by a Gaussian random variable with expectation zero and standard deviation </20, i.e., we add a realization to the gray value of all voxels independently. Clamping ensures that all gray values remain in the interval [0,255] (Fig. 18). This leads to data similar to the real 3D data set, at least by visual inspection.
COMPARISON OF INPUT AND OUTPUT
A visual comparison of the detection results shows that the algorithm works quite well (Fig. 19). To compare the sets of spheres used as input and obtained as output, we match the spheres using a criterion based on their overlapping. For every "input" sphere, the sphere with the highest overlapping from the output is detected. Analogously, for every sphere of the output we detect the input sphere with the highest overlapping. In both cases an overlapping of at least 5% is required. Then, we consider two spheres to be a match only if both have been assigned to each other. A very small minimum percentage is used because we want to match as much spheres as possible, even if they do not match perfectly. Using this approach, about 95% of all spheres are matched. For the pairs of matched input and output spheres, we can compute the distance of their centers and the difference of
Fig. 18. 2D sections of 3D image based on the extracted spheres, once shown as perfect spheres (top) and after Gaussian blur with white noise (bottom).
Fig. 19. Validation of sphere-extraction, once with the expected spheres (top) and once with the image used as input for the algorithm (bottom).
la 0.3
0 i 2 3 4 center distances in voxels
la 0.3
-4-3-2-1 0 1 2 3 4 radius difference in voxels
Fig. 20. Distances between centers (top) and differences between radii (bottom) for all pairs of matched spheres.
In view of these results it is not surprising that also the sample mean and variance of the radii match nearly perfectly (Table 2).
Table 2. Analysis of radii in voxels for validation of the sphere-extraction algorithm.
	input radii	output radii
sample mean	7.82	7.81
sample variance	3.44	3.43
CONCLUSIONS
In the present paper, we have developed and validated an algorithm for the detection of CdSe nanoparticles in CdSe-(OCiCi0-PPV) composites based on 3D grayscale images. The considered cutout of the 3D image was a representative sample of a dense-packed region of CdSe nanoparticles. In a first step, we have filtered the 3D image data and applied a morphological closing to reduce noise and to enhance edges of the objects. The developed iterative thresholding procedure combined with a cluster detection is applied subsequently to locate potential centers of CdSe nanoparticles. For their precise localization, we decomposed the image domain into pairwise disjoint basins, where in each basin exactly one object is detected using a modified Hough transform for sphere extraction. A post-processing of detected particles reduces overlapping artifacts. A
second run of the same procedure is considered where, after removing the previously detected spheres, some additional particles are found. They are then added to get the final sphere representation of the system of CdSe nanoparticles.
The second part of the paper concerns the statistical analysis of the system of detected particles. The spherical particles are interpreted as a marked point pattern, wherefore tools from point-process statistics are applied. In addition, we have connected neighboring CdSe nanoparticles to form 3D graphs, where the fact that electrons can cross small gaps between disjoint nanoparticles by hopping is taken into account. The considered characteristics are mainly focused on describing the connectivity of the graphs, which is strongly related with the transport processes of electrons within the considered organic-inorganic composite.
The results obtained in the present paper can be considered as a first step towards the development of spatial stochastic models for systems of nanoparticles in organic-inorganic composites. By means of such models, non-interactive statistical methods can be established in order to detect and quantify structural advantages/disadvantages of various types of nanoparticle systems in organic-inorganic composites. Moreover, model-based Monte Carlo simulations of virtual materials can be accomplished, where the general intention of such computer experiments is the virtual design of materials with improved physical properties.
ACKNOWLEDGEMENTS
This research has been initiated while Volker Schmidt was visiting the Isaac Newton Institute for Mathematical Sciences, Cambridge, in March/April 20i0. Part of this work has been supported by Deutsche Forschungsgemeinschaft (DFG) under the Priority Programme: 'Elementary Processes of Organic Photovoltaics' (SPP 1355). The imaging work was supported by the Engineering and Physical Sciences Research Council, U.K., and by the European Union under the Framework 6 program under a contract for an Integrated Infrastructure Initiative, Reference 026019 ESTEEM. We thank the anonymous referees for their helpful comments.
REFERENCES
Beare R, Lehmann G (2006). The watershed transform in ITK—discussion and new developments. Insight J (Online) http://hdl.handle.net/1926/202 Ballard DH (1981). Generalizing the Hough transform to detect arbitrary shapes. Pattern Recogn 13(2):111—22.
5
Beucher S, Lantuejoul C (1979). Use of watersheds in contour detection. In: Proc Int Worksh on Image Processing, Real-Time Edge and Motion Detection/Estimation. Rennes, France.
Beucher S, Meyer F (1993). The morphological approach to segmentation: the watershed transformation. In: Dougherty ER, Ed. Mathematical Morphology in Image Processing. New York: Marcel Dekker. 433-81.
Burger W, Burge MJ (2008). Digital Image Processing - An Algorithmic Approach using Java. Berlin: Springer.
Dayal S, Reese MO, Ferguson AJ, Ginley DS, Rumbles G, Kopidakis N (2010). The effect of nanoparticle shape on the photocarrier dynamics and photovoltaic device performance of Poly(3-hexylthiophene):CdSe nanoparticle bulk heterojunction solar cells. Adv Funct Mater 20:2629-35.
Decker L, Jeulin D, Tovena I (1998). 3D morphological analysis of the connectivity of a porous medium. Acta Stereol 17(1):107-112.
Diestel R (2005). Graph Theory. Heidelberg: Springer.
Diggle PJ (2003). Statistical Analysis of Spatial Point Patterns. London: Arnold.
Duda RO, Hart PE (1972). Use of the Hough transformation to detect lines and curves in pictures. Commun ACM 15:11-5.
Gelfand AE, Diggle PJ, Fuentes M, Guttorp P (2010). Handbook of Spatial Statistics. London Boca Raton: Chapman & Hall/CRC.
Greenham NC, Peng X, Alivisatos AP (1996). Charge separation and transport in conjugated-polymer/semiconductor-nanocrystal	composites
studied by photoluminescence quenching and photoconductivity. Phys Rev B 54:17628-37.
Hindson JC, Saghi Z, Hernandez-Garrido JC, Midgley PA, Greenham NC (2011). Morphological study of nanoparticle:polymer solar cells using high-angle annular dark-field electron tomography. Nano Lett
11:904-9.
Hoshen J, Kopelman R (1976). Percolation and cluster distribution. I. Cluster multiple labeling technique and critical concentration algorithm. Phys Rev B 14:3438-45.
Huynh WU, Dittmer JJ, Alivisatos AP (2002). Hybrid nanorod-polymer solar cells. Science 295:2425-7.
Illian J, Penttinen A, Stoyan H, and Stoyan D (2008). Statistical Analysis and Modeling of Spatial Point Patterns. Chichester: J. Wiley & Sons.
Jahne B (2005). Digital Image Processing, 6th revised and extended edition. Berlin: Springer.
Jungnickel D (1999). Graphs, Networks and Algorithms. Berlin: Springer.
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.
Peyrega C, Jeulin D, Delisee C, Malvestio J (2009). 3D morphological modelling of a random fibrous network. Proc 10th Eur Congr Stereol Image Anal, Milan. 12530.
Roerdink JBTM, Meijster A. (2000). The watershed transform: Definitions, algorithms, and parallellization strategies. Fundam Inform 41:187-228.
Serra J (1982). Image Analysis and Mathematical Morphology Vol. I. London: Academic Press.
Thiedmann R, Hartnig C, Manke I, Schmidt V, Lehnert W (2009). Local structural characteristics of pore space in GDL's of PEM fuel cells based on geometric 3D graphs. J Electrochem Soc 156:B 1339-47.
Wang P, Abrusci A, Wong HMP, Svensson M, Andersson MR, Greenham NC (2006). Photoinduced charge transfer and efficient solar energy conversion in a blend of a red polyfluorene copolymer with CdSe nanoparticles. Nano Lett 6:1789-93.