Two-dimensional illustration

A sample "aggregate"shape to be analyzed is shown in Fig. 1, which was drawn by hand using a drawing program and then saved as a digital file. This file was converted to an ASCII file of ones and twos (one for background, two for object) using a shareware image processing program. The shape is then represented in the computer as a collection of pixels or list of numbers, located at specific sites in the plane. For a real particle, with an actual size, there will be a certain real length associated with the side length of each square pixel.

Figure 1: A 2-D model aggregate shape. The angle $\theta$ used in the Fourier analysis is defined. The origin of the coordinate axes is at the center of mass of the object.

The center of mass of the object is found using the following equation:

$$\vec{r}_{CM} = \frac{\sum_i \vec{r}_i}{\sum_i 1}$$ (1)

where the sum over i is a sum over all pixels of the object (pixels with label "2"), the position vector is defined as the vector from the (arbitrary) origin to the center of the pixel of interest, and the symbol in the denominator is the numeral one. The denominator is then simply the area of the object in terms of the number of pixels contained in the object.

The surface of the object is found next, numerically, at a given number of angles by extending a line segment from the center of mass to a point that just crosses the digital surface of the object, which is taken to be the edge of the pixel found at the surface in the direction chosen. The lengths of these line segments are denoted R_n. The angles, between 0 and 2$2\pi$, are chosen according to a Gaussian quadrature scheme [29,30] of the order desired, and are illustrated in Fig. 1. This method of finding the surface will cause the numerical surface to contain the small-scale "digital roughness" of the pixellated surface. It is conceivable that some sort of interpolation or smoothing operation could be done at this point, to get rid of the digital roughness to some extent. This has not been done in this paper, for 2-D or for 3-D.

One must note here that this method of finding the surface will only work for a certain class of shapes. Defining this class precisely is left for the 3-D section. For now, it can be simply said that when extending a line segment from the center of mass to the surface, one must not intersect the surface twice. Figure 2 shows a particle that has two features that would invalidate this procedure: an "overhanging" feature, and internal porosity that is resolved on the length scale of the pixel size.

Figure 2: A 2-D model aggregate shape that illustrates the features that are not allowed in the mathematical analysis, including resolved internal porosity, and an "overhanging" region, giving a non-single-valued surface, as measured by directed segments from the center of mass.

The Gaussian quadrature is a method for doing integrals numerically [29]. If a function f is to be integrated from -1 to 1, and a Gaussian quadrature of order N is chosen, then

$$\int_{-1}^1 f(x) dx \approx \sum_{n=1}^N f_n w_n$$ (2)

where the function is evaluated at the points x_n, -1 < x_n < 1, of the Gaussian quadrature, and w_n are the weights of the Gaussian quadrature of that order. In this paper, a 120 point Gaussian quadrature was generally used. If the integration limits are instead a and b, then one uses the linear transformation

$$x_n^{/} = \frac{1}{2}(b-a) x_n + \frac{1}{2}(b+a)$$ (3)

The function is evaluated at these new points, and the sum in eq.(2) is used to evaluate the integral with the same weights. If a = 0 and b = 2/$b=2\pi$, the limits for $\theta$ in the surface problem, then this transformation gives the values of $\theta_n$_n used, for a given value of 0 n $0 \le n \le N$ N.

Now the function R() = ($R(\theta)={(\theta_n, R_n)}$_n , R_n ) is the length from the center of mass to the surface at an angle $\theta_n$_n, where $\theta_n$_n is measured counterclockwise from the x-axis. This function can be constructed to be periodic in $\theta$, by letting _{N + i} = _i + 2$\theta_{N+i}=\theta_i + 2\pi$, and so can be expanded in terms of cosine and sine functions:

$$R(\theta) = \sum_{j=0}^{\infty} [a_j \cos (j\theta) + b_j \sin (j\theta)].$$ (4)

Note that the value of b₀ can be taken to be zero, as sin(0) equals zero identically. Also note that eq. (4) could be written in terms of a complex exponential, which combines the sine and cosine terms.

The values of a_j and b_j are given by the integrals

$$a_j = \frac{1}{2 \pi}\int_0^{2\pi} R(\theta) \cos (j\theta) d\theta$$

$$b_j = \frac{1}{2 \pi}\int_0^{2\pi} R(\theta) \sin (j\theta) d\theta$$ (5)

These integrals are straightforward to do, since the values of $\theta_n$_n were chosen according to the Gaussian quadrature procedure, and so are "ready-made" to do the integrals.

Figure 3 provides a plot of the values of a_j and b_j vs. jfound for the shape shown in Fig. 1. Note that, while there is fluctuation, the values of the coefficients tend to decrease as the value of j increases. By j = 20 the coefficients have become small enough to be negligible. After this point, higher order values of j are just trying to match the pixel-to-pixel digital structure of the exterior of the object, which is probably not an important or realistic part of the real shape.

Figure 3: The cos and sin coefficients as a function of the index j for the shape shown in Fig. 1. The value of a₀ 200 is off the vertical scale so that the smaller coefficients could be seen.

Figure 4 shows the computed shape of the object for different numbers of Fourier coefficients used, i.e.,

$$R(\theta) \approx \sum_{j=0}^{N} [a_j \cos (j\theta) + b_j \sin (j\theta)].$$ (6)

Note that the a₀coefficient is a measure of the average circular size of the object, as can be seen also in eq. (5), and the other coefficients (j $j \geq 1$ 1) can be thought of as perturbing the average circular surface to match the correct shape.

Figure 4: The shape of the analyzed image of Fig. 1 for different numbers of cos and sin coefficients used to re-create the shape. The numeric labels in the figure are equivalent to the j values in the previous figure.

Figure 5 shows the computed area of the shape, using the expansion, as a function of how many expansion coefficients were used. The equation for the area A of a given shape, using N coefficients, is [31]:

$$A = \frac{1}{2} \int_0^{2\pi} [r(\theta)]^2 d\theta$$

$$= \pi a_0^2 + \frac{\pi}{2} \sum_{j=1}^N (a_j^2 + b_j^2)$$ (7)

Note that by N = 10 or so, the area has reached its asymptotic value, which is less than the digital value by only 0.1 %. The digital area is just the number of the pixels times the area of each pixel. In this model case, the area of each pixel is taken to be unity. The slight disagreement reflects the fact that the surface of the object was only sampled at a finite number of points. Choosing surface points according to a higher order Gaussian quadrature would cause this error to be smaller, as long as the numerical accuracy of the computer was not exceeded. However, one does not want to totally reproduce the digital details of the pixellated surface, as presumably these details are not real, but are only an artifact of using square pixels to represent a true curved surface. So there is an upper limit on how fine one wants to resolve the surface.

Figure 5: The computed area of the shape shown in Fig. 2
as a function of the index j.

The technique described for 2-D has been: image acquisition, numerical surface determination, and then Fourier analysis. In 3-D, the image is (1) generated using x-ray tomography, (2) acquired numerically, and then (3) analyzed with spherical harmonic functions, which are the 3-D equivalent of the Fourier series. In the case of the asteroid mentioned earlier [23], the surface was numerically sampled using laser range finder measurements made on an orbiting satellite. The x-ray tomography and particle image acquisition step for the procedure in 3-D is described next.