Form factor (Fourier shape transform) of polygon and polyhedron

The Fourier transform of the indicator function of arbitrary polygons and polyhedra is computed for complex wavevectors. Using the divergence theorem and Stokes' theorem, closed expressions are obtained. Apparent singularities, all removable, are discussed in detail. Loss of precision due to cancellation near the singularities can be avoided by using series expansions.

1. Introduction. The term "form factor" has different meanings in science and in engineering. Here, we are concerned with the form factor of a geometric figure as defined in the physical sciences, namely the Fourier transform of the figure's indicator function, also called the shape transform of the figure.
This form factor has important applications in the emission, detection, and scattering of radiation. Two-dimensional shape transforms are used in the theory of reflector antennas [1]. The shape transform of three-dimensional nanoparticles is used to interpret neutron and x-ray small-angle scattering (SAS, SANS, SAXS) [2]. A particularly rich multitude of mostly polyhedral shapes is observed for particles grown on a substrate [3], as observed by grazing-incidence neutron and x-ray small-angle scattering (GISAS, GISANS, GISAXS) [4]. Thence it comes that possibly the most extensive collection of particle shape transforms published to this date can be found in the documentation [4,5] of a GISAXS simulation software called IsGISAXS. The oldest specimens are the form factors of a cylinder and of a sphere, worked out by Lord Rayleigh in his electromagnetic theory of light scattering [6].
Analytic expressions for form factors of polyhedra can be obtained by straightforward triple integration, using Fubini's theorem. However, except for very simple shapes like a rectangular box, oriented along the coordinate axes, the resulting expressions look more complicated than one might expect for such a plain problem; they do not reflect the symmetry of the geometric figure, and they contain removable singularities that cause two difficulties: Division by zero at the singularities, unless analytic continuations are used, and loss of arithmetic precision due to cancellation of terms near the singularities, unless special precautions are taken.
Deriving and implementing numerically stable algorithms for a considerable collection of Platonic solids, prisms, truncated prisms, pyramids, frusta, and more would be an immense labor. It is preferable to derive the shape transform of arbitrary polygons and polyhedra in generic terms, and deal once and for all with all possible singularities. For polygons, a generic expression for the form factor has been derived decades ago [1]; singularities were not discussed. For polyhedra, a first-order approximation to the form factor is used in a novel GISAXS simulation software, called HipGISAXS [7]. Here, exact expressions for the form factors of polygons and polyhedra will be derived. Proper treatement of all singularities will lead to practicable algorithms for computing at arbitrary wavevectors. All results presented here have already found an application in another novel GISAXS simulation software, called BornAgain [8].
For a Fourier transform in the usual sense, and for most applications, wavevectors are real: q ∈ R 3 . In GISAS, however, the incident and scattered radiation may undergo substantial absorption, which can be modeled by an imaginary part of q. Similarly, an imaginary part of q can be used to described light amplification in a laser medium. Therefore, we admit complex wavevectors q ∈ C 3 .
Polygonal and polyhedral form factors have strong similarities in their mathematical structure, reflecting the analogy of Stokes' theorem and the divergence theorem. This paper strives to bring out these similarities by a purposedly parallel treatment of the two-and the three-dimensional case. It is organized as follows: Section 2 provides some basics for computing the form factor of a planar figure embedded in three-dimensional space. Section 3 deals with the form factor of a polygon; it provides a closed analytic expression and power series coefficients. Section 4 treats the form factor of three-dimensional figures; section 5 specializes to a polyhedron. Section 6 works out simplifications for figures with inversion centers. TODO section 7 As for notation: Vector products (dot products and cross products) are linear in both components; complex conjugation, where needed, is denoted explicitely. And of course 0 0 = 1 [9].

Form factor of a planar figure.
Definition 1 (oriented plane). A plane, given by a normal unit vectorn and a signed distance r ⊥ from the origin, shall be denoted as The orientation of the plane is induced byn as follows: If the ordered triple (b 1 , b 2 ,n) of pairwise orthogonal vectors is a positively oriented (right-handed) base of the R 3 , then the ordered pair (b 1 , b 2 ) is a positively oriented base of E(n, r ⊥ ).
For a loop in the plane E(n, r ⊥ ), the following statements are equivalent: (1) The loop has a positive winding number. (2) The loop runs counterclockwise. (3) The loop fulfills the right-hand-rule with respect ton.
Definition 2 (form factor of a 2d figure in 3d space). The form factor of a twodimensional planar figure Γ ⊂ E(n, r ⊥ ), embedded in three-dimensional space, at a wavevector q ∈ C 3 is where r := r ⊥ + r , r ⊥ := r ⊥n , and r is given by the integration variables.
Definition 3 (vector decomposition induced byn). In a given plane, characterized by a normal vectorn, the subscripts ⊥ and , affixed to a vector v ∈ C 3 , indicate the projection of v onton ( perpendicular to the plane) respectively. The subscript × shall denote an in-plane vector normal to v, This definition is compatible with, and extends, the notation of Definitions 1 and 2 where r was constructed as r ⊥ + r .
If it happens that a given vector v is exactly or almost parallel ton, then the computation of v according to (4) becomes inaccurate, and can result in v having a relatively strong out-of-plane component. To reduce this spurious component, one can iterate Furthermore, if |v |/|v ⊥ | is smaller than the machine epsilon, then it is adequate to let v = 0.
Proposition 4 (factorization of f ). The form factor (2) of a planar figure Γ ⊂ E(n, r ⊥ ) can be factorized as Proof. Go back to the definition (2) of f (q, Γ): The factor f (q , Γ) shall be called the in-plane form factor of Γ.
Proposition 5 (continuity at q = 0). The in-plane form factor of a planar figure Γ ⊂ E(n, r ⊥ ) is continuous at q = 0, and has the limit where Ar(Γ) is the area of Γ.
Proof. The continuity and the first equality in (9) hold because for q → 0 and r ∈ Γ the convergence e iq r → 1 of the integrand of Definition 2 is uniform. The second equality is Definition 2 with q = 0.
Proposition 6 (series expansions). The form factor of a planar figure Γ possesses for any q ∈ C 3 two different absolutely convergent series expansions which coincide for in-plane wave vectors, f n (q ) = φ n (q ). For generic q, Proof. Expand the exponential in the integrand of (2) in q or q to obtain (10a) or (10b). Define the radii (13) a := max j |V j | and b := max j |V j | of a sphere respectively a circle that contain all vertices of Γ. Then proves the absolute convergence of both series. The factorization of f n in (12a) is obtained by straightforward binomial expansion of (q r + q ⊥ r) n .
Remark 7 (choice of origin). It may be advisable to translate the origin before starting a numeric computation of f . This can be done using a standard property of Fourier transforms Two particular choices of the origin have special advantages: An origin at the center of the enclosing circle minimizes a, whereas an origin at the center of gravity lets the expansion coefficient f 1 vanish.
Remark 8 (termination of numeric summation). In general, the integral formula (11a) is no practicable starting point for computing the f m (q , Γ). It allows us, however, to make the following observation: For odd m, but not for even m, it may happen that f m = 0 for some q = 0. In particular, f 1 ≡ 0 for all q if the origin is chosen at the center of gravity. Therefore the termination criterion in a numeric implementation of (10a) or (10b) must only rely on terms with even m.

Form factor of a polygon.
Definition 9 (simple polygonal vertex chain). A simple polygonal vertex chain of size J is a sequence of points V 1 , . . . , V J , made cyclic by the convention V 0 ≡ V J , that lie all in a plane, and when connected by edges V 0 V 1 , . . . , V J−1 , V J , concatenated in this order, form a non-intersecting loop.
The plane E(n, r ⊥ ), given by V 1 , . . . , V J , is oriented such that ∂Γ has a winding number +1 with respect ton.
Throughout the following, we assume that the origin is chosen such that V j = 0 for all j ∈ {1, . . . , J}. For a given polygonal vertex chain, the normal vector can then be computed as and the location of the plane as with arbitrary j ∈ {1, . . . , J}. If vertex coordinates are ill-conditioned, then it may be indicated to use (16) and (17) with a special choice of j, or to average over several j. Proposition 10 (form factor of a polygon). A J-gon Γ ⊂ E(n, r ⊥ ) be given by a simple polygonal vertex chain Proof. We parametrize the edges of the polygon by r j (λ) := R j + E j λ with −1 ≤ λ ≤ +1. Stokes's theorem then takes the form With the choice G := q * × e iqr , this yields With q = 0: which can easily be brought into the form (18). Figure 1 illustrates the so determined form for about the simplest polygon, an equilateral triangle. |f (q)| is plotted as function of q for three different directionsq.
Remark 11 (sinc is accurate). The cardinal sine function is best implemented by literally following this definition. Any available implementation of sin(z) will have full floating-point accuracy for |z| → 0, and so sin(z)/z will have.
Remark 12 (edge-independent terms cancel). Since the vertices of Γ form a closed loop, the sum of the edge vectors vanishes, The sum over j in (18) contains a factor sinc(qE j ) exp(iqR j ). The leading term in a q -expansion of this factor is exp(iq ⊥ r ⊥ ). It can be drawn in front of the sum, which then vanishes per (23). If the next-to-leading terms are relatively small, then the cancellation of the leading term will cause a loss of accuracy in the resulting form factor. Therefore, the analytic result (18) is not practicable for small q , where it must be replaced by the series expansion (10b). Alternatively, if not only q but also q ⊥ is small, then the q-expansion (10a) can be used. Expansion coefficients will be provided below in Proposition 15.
The second panel of Figure 1 demonstrates how this cancellation shows up in a double-precision implementation of (18). For q L close to or below the machine epsilon of 2 · 10 −16 , resulting values for some q are wrong by O((Lq ) −1 ).
Remark 13 (area formula). The case q = 0 is covered by Proposition 5. The area of a polygon can be conveniently computed using the surveyor's formula [10], which is based on a triangular tesselation.
Remark 14 (relation to literature result). A closed expression for the form factor of the polygon is known since long [1, Eq. 6]. In our notation and after a few obvious simplifications, it reads This expression is beautiful for being short and symmetric. However, for each j there are two q planes for which the denominator vanishes. A practical implementation of (25) that takes proper care of these singularities would be more complicated than an algorithm based on our Proposition 10.
Proof. Let us demonstrate the equivalence of (25) with our Proposition 10. Start from Equation (21). In the exponential functions, insert the definitions of E j and R j :  Figure 1, as function of wavenumber q for wavevector directionq = (0, 1, 0). The chain of black circles is computed using the analytic expression (18). The colored curves are computed using the series expansion (10a) with coefficients (28) up to the indicated orders n.
Shuffle indices j + 1 → j, use qE j = q E j , and employ standard vector identities: Insert the definition of q × , employ another vector identity, usenq = 0, and obtain (25).
Proposition 15 (expansion coefficients). The coefficients of the series expansion (10a) of the form factor can be computed as Proof. Expand the functions sinc(qE) and exp(iqR) in (18): Sort by powers of q: For ν = 0, the sum over l yields the constant 1, which cancels per (23) under the sum over j. Therefore the outer sum may as well start at ν = 1. Substitute n = ν − 1 to obtain the series (10a) with coefficients (28). Figure 2 shows that the series expansion works well even beyond the first few oscillations in f (|q|), provided a sufficient number of terms is summed.
Remark 16 (area formula from expansion coefficient). From (11a) we know that f 0 (q, Γ) = Ar(Γ). On the other hand, per Proposition 15, It was not necessary for the above proofs, but may nevertheless be interesting to explicitly show how (31) can be algebraically transformed into the q-independent form (24). Insert the definitions of E j and R j into (31): Multiply out, and shuffle indices (j −1) → j for some terms under the sum: From here, simplify using the same standard vector identities as were employed in Remark 14. Use f 0 (q, Γ) = Ar(Γ) to obtain (24).
Remark 17 (resulting algorithm). Let us summarize sections 2 and 3 by outlining an algorithm for reliably computing the form factor of any polygon Γ, given through an oriented vertex chain. The algorithm involves preselected constants c and c that determine when to use a series expansion instead of the analytic formula. Do not use this algorithm if the symmetry of Γ allows for a simpler computation, as discussed below in section 6.
If the input coordinates cannot be trusted: Test whether all vertices lie in a plane; otherwise terminate with error message. Check whether the coordinate origin lies well inside Γ; otherwise determine a new origin, for instance from the center of gravity of Γ, and apply (15) at the end of the computation.
Compute E j , R j ,n and r ⊥ , and decompose q into q + q ⊥ . If q = 0, then return f (q, Γ) = e iq ⊥ r ⊥ Ar(Γ). If aq < c, then sum the series (10a) with f 0 (q, Γ) = e iq ⊥ r ⊥ Ar(Γ), and other coefficients f n (q, Γ) given by (28). Terminate the summation and return the result if a term for even n, relative the absolute value of the sum acquired so far, falls below the machine epsilon. If aq < c , then sum the series (10b) with coefficients f n (q , Γ) still given by (28). Apply the same termination criterion as in the previous case, and return the result. In the remaining cases, return the result of the direct computation (18).

Form factor of a three-dimensional figure.
Definition 18 (form factor). The form factor of a figure Π ⊂ R 3 at a wavevector q ∈ C 3 is the Fourier transform of its indicator function, given by the integral Proposition 19 (continuity at q = 0). The form factor of a figure Π is continuous at q = 0, and has the limit where Vol(Π) is the volume of Π.
The proof is fully analogous to that of Proposition 5. The later Lord Rayleigh had already noted about the scattering of light by small particles [11]: "The leading term, we have seen, depends only on the volume ; but the same would not be true for those that follow . . . there is no difficulty in proceeding further; but I have not arrived at any results of interest." Proposition 20 (series expansion). The form factor of a figure Π possesses for any q ∈ C 3 the absolutely convergent series expansion Proof. Expand the exponential function in the integrand of (34) to obtain (36). Use the radius a defined in (13) to derive the bound This proves the absolute convergence of the series and thereby the existence and uniqueness of (36).
Remark 21 (termination of numeric summation). The same considerations as in Remark 8 apply: coefficients F n may vanish for odd n; therefore the termination criterion in a numeric implementation of (36) must only rely on terms with even n.

Form factor of a polyhedron.
Proposition 22 (form factor of a polyhedron). An orientable polyhedron Π be given by its K faces Γ k (k = 1, . . . , K). Each face Γ k ∈ E(n k , r ⊥k ) be an J k -gon, given by the simple polygonal vertex chain V k1 , . . . , V kJ k , and withn k pointing towards the outside of Π. Then the form factor for q = 0 is Proof. We first address the case q = 0. For a polyhedron, the divergence theorem takes the form With the choice H := q * e iqr , this yields With the notation (2), this proves (39).
Remark 23 (volume formula). The case q = 0 is covered by Proposition 5. The volume of a polyhedron can be conveniently computed from This formula, which is based on a tetrahedral tesselation. can be retrieved from remote literature [12, Cap. II, § 3, III 171], or easily be derived by applying the divergence theorem (40) with the choice H := r/3.
Remark 24 (face-independent contributions cancel). For a closed surface, the integral of surface normals cancels. For a polyhedron, the integral becomes a sum: For proof, use the divergence theorem (40) with a k-and r-independent H.
In full analogy with Remark 12, this cancellation can lead to a loss of accuracy in (39) if f (q, Γ k )/ Ar(Γ k ) is dominated by k-independent terms. According to the series expansion (10a), this is indeed the case; the leading term is just f 0 (q, Γ k ) = Ar(Γ k ). Therefore, the analytic result (39) is not practicable for small q, where it must be replaced by the series expansion (36). Coefficients, unaffected by cancellation, will be provided below in Proposition 25.
Proposition 25 (expansion coefficients). The coefficients of the series expansion (36) can be computed as Proof. Insert (10a) in (39): Per Proposition 6, the f ν are absolutely convergent. Therefore the two sums can be exchanged. Per Remark 24, the term f 0 (q, Γ k ) = Ar(Γ k ) cancels under the sum over k. Therefore the infinite sum may as well start at ν = 1. Substitute n = ν − 1 to obtain the series (36) with coefficients (44). Figure 3 illustrates the crossover between the domains of analytic computation and series expansion. The implementation, described below in section 7, uses IEEE double-precision arithmetics. The analytic form factor (39) with per-face contributions (18) fails grossly for qa 10 −5 , as was to be expected from Remark 24. On the other side, the series expansion (36) requires more and more terms if qa 10. This leaves a wide intermediate range where both methods agree well enough to allow a match with decent accuracy, as further discussed in section 7. Remark 26 (resulting algorithm). Let us summarize sections 4 and 5 by outlining an algorithm for reliably computing the form factor of any given polyhedron Π, given through a set of K polygonal faces Γ k . The algorithm involves a preselected constant C that determines when to use the series expansion instead of the analytic formula. Do not use this algorithm if the symmetry of Π allows for a simpler computation, as discussed below in section 6.
If the input coordinates cannot be trusted, check whether the origin lies well inside Π. Otherwise determine a new origin, for instance from the center of gravity of Π, and apply (15) at the end of the computation.
If q = 0, then return F (0, Π) = Vol(Π). If aq < C, then sum the series (36), using F 0 = Vol(Π) and other coefficients from (44). Terminate the summation and return the result if a term for even n, relative the absolute value of the sum acquired so far, falls below the machine epsilon. Otherwise compute F (q, Π) according to (39) as a weighted sum over face form factors. To compute these f (q, Γ k ), follow the algorithm described in Remark 17.
6. Symmetric figures. For figures with certain symmetries, the form factor computation can be considerably simplified. Besides the speed benefit, this also improves the accuracy by setting cancelling terms exactly to zero.
Remark 27 (polygon with symmetry S 2 ). If a planar (2J)-gon Γ has a perpendicular twofold symmetry axis (Schoenflies group S 2 ), and thereby an inversion center at point r ⊥ , then its form factor has the symmetry f (q ⊥ + q ) = f (q ⊥ − q ). Its real-space coordinates transform as V j+J = −V j for j = 1, . . . ,J, and thereby R j+J = −R j and E j+J = −E j . This allows the following simplifications: The area can be computed as The form factor for q = 0 can be computed as In contrast to (18), the term under j E j has no constant contribution, but is of order q R j . There is no cancellation for q → 0, and no need to use a series expansion for the accurate computation of f .
Remark 28 (polyhedron with symmetry C i ). If a (2K) has an inversion center at the origin 0 (Schoenflies group C i ), then its face indices can be chosen such that Γ k+K = 0 − Γ k , and the face form factors have the symmetry f (q, Γ k+K ) = f (−q, Γ k ) = f (q, Γ k ) * . The decomposition of q is invariant: q k+K⊥ = q k⊥ , and idem for q k . In contrast, q k+K× = −q k× because it involvesn k+K = −n k . With this, the computation of F (q, Π) can be simplified as follows: The volume (42) can be computed as The analytic form factor (39) for q = 0 can be computed as E j sinc(qE j ) cos(qR j ), obtained from (7) and (18). If Γ has symmetry S 2 they become If q k is small, then thef (q, Γ k ) should be computed from the series expansion From (11a) we see that the f n are even/odd functions of q (but not of q) for even/odd n. Thereby (54)φ n (q, Γ) = 2i sin(q ⊥ r ⊥ )f n (q , Γ) for even n, 2 cos(q ⊥ r ⊥ )f n (q , Γ) for odd n, with f n computed using (28). In the remaining case of small q, the coefficients of the series (36) take the form which offers hardly any advantage over the original form (44).
Remark 29 (prism). For prisms, it is advisable to factorize (34) from the onset. Let a prism Π have a base Γ that determines a normaln. Let the prism extend in normal direction from −h/2 to h/2. Then for all q Per Remark 11, no series expansion is needed for q ⊥ → 0.
7. Numeric tests. A form factor computation for generic polyhedra, based on all the above, has been implemented as part of the GISAS simulation package BornAgain [8,13]. The source code is available under the GNU General Public License. All floating-point numbers, internally and externally, have double precision. Tests were performed on standard PC's with amd64 architecture under Linux, hence in IEEE arithmetics.
The code underwent extensive tests for internal consistency and for compatibility with previous computations. Those reference computations were based on analytical expressions obtained by three successive integrations in Cartesian coordinates; they had been checked against the reference code IsGISAXS [4,5].
The internal tests comprise the symmetry, specialization, and continuity tests. All these tests were performed, as far as applicable, for a suite of particle shapes (pyramidal frusta with 2-, 3-, 4-, and 6-fold symmetry, cuboctahedron, truncated cube, regular dodecahedron and icosahedron), for different wavevector directionsq (along symmetry axis, slightly off such axis, or in completely unsymmetric directions; purely real or with small imaginary parts), and for a logarithmically wide range of magnitudes q. The main result of these tests is the worst-case relative deviation between two form factor computations, Symmetry tests are performed for particle shapes that are invariant under some rotation or reflection R. They yield δ[F (q), F (Rq)] 5 · 10 −10 .
In specialization tests, special parameter sets are chosen for which two otherwise different figures Π 1 , Π 2 coincide. For instance, the rectangular base of a pyramid Π 1 is made a square, so that F (q, Π 1 ) can be compared with the form factor of a square pyramid Π 2 . Or the dihedral angle of a pyramid is set to 90 • so that the pyramid coincides with a prism. These tests yield δ[F (q, Π 1 ), F (q, Π 2 )] 3 · 10 −10 .
Finally, continuity tests search for possible discontinuities due to a change in the computational method. Using a precompiler switch, additional source code lines are activated that tell the test program whether the analytic expression or a series expansion was used to compute a form factor, and in the latter case, at which order the summation was terminated. For given directionq, bisection in q is used to determine where such a change in computation method happens. Then, the form factor slightly below and slightly above this threshold is determined. With η = 8 · 10 −16 chosen as a few times the machine epsilon, the continuity measure is δ cont := δ[F (q(1 − η)), F (q(1 + η))]. According to Remark 17, the switch between series expansion and analytic expression is determined by a parameter c. The optimum value of this parameter is determined empirically so that δ cont is miminimized.