Detection of symmetry using a crystallographic image processing algorithm
DD ETECTION OF SYMMETRY USING A CRYSTALLOGRAPHICIMAGE PROCESSING ALGORITHM . A P
REPRINT
Paul Plachinda
Department of PhysicsPortland State UniversityPortland, OR 97201 [email protected]
August 22, 2019 A BSTRACT
This article presents an automated method to quantify and detect symmetry elements in 2D patterns bymeans of image processing. Escher’s woodcuts, a widely recognized didactic tool for crystallographiceducation of students, were used to demonstrate this approach. We also discuss peculiarities in thedetection of black and white symmetry, color symmetry, and detection of the ”hidden” and ”broken”symmetry elements by means of the phase origin map approach. K eywords Symmetry Detection · Symmetry Enforcement · Phase Origin Map
The legacy of the Dutch graphic artist M.C. Escher (1898 - 1972) consists of more than 500 woodcuts, many of whichexhibit exceptionally peculiar symmetric patterns, which combine outstanding artistic beauty with a carefully elaboratedsymmetry. These two unique points have made them a valuable tool for teaching the foundations of symmetry inmultiple elementary textbooks on crystallography.[1, 2] Besides tessellations of the plane Escher’s drawings demonstratemore complex concepts, such as black and white, and color symmetry. The founders of these concepts, Shubnikov andBelov[3, 4] also used Escher’s drawings to illustrate their pioneering ideas. The physical application of these conceptsis magnetic symmetry, where the black and white groups correspond to spin up/spin down states of the same atom inthe unit cell. Color symmetry is used to describe more complex spin arrangements.Teaching symmetry to students can be sometimes complicated since the mathematical representation of it as thought atthe university level sometimes contradicts student’s naive perception. Several papers (e.g.[5] and [6]) have previouslyaddressed comprehension difficulties and provided simple models to understand cornerstone concepts of crystallographyand symmetry.Visual perception of symmetry by humans can be illustrated by several examples:[7] people tend to value more and payattention to high-symmetry elements and perfectly repeating patterns, and to ignore features violating such. Violationof translational symmetry is considered exceptionally more obvious and disturbing rather than violation of rotational ormirror symmetry elements. It is interesting to observe, that traditional Japanese aesthetics (Wabi-sabi) values asymmetrymore than symmetric patterns associated with the Greek ideals of beauty and perfection in the West.[8]From the author’s personal teaching experience, students tend to envision symmetry in two different ways: an orthodoxway - if there is a single feature violating proposed symmetry, the overall symmetry is dictated by this feature; and arelaxed way - if there is enough repeating motifs of the same symmetry, then the image possesses the symmetry of anundistorted motif. Real crystals never possess perfect symmetry as described by their mathematical model due to point All M.C. Escher works c (cid:13) a r X i v : . [ phy s i c s . pop - ph ] A ug PREPRINT - A
UGUST
22, 2019and line defects, grain boundaries, magnetic domains, etc. This paper is aimed to elucidate the concepts of symmetryquantification, refinement, and detection of the "hidden" and "broken" symmetry elements.For more detailed consideration we chose four examples: three Escher’s plane tessellations (Figures 1a,5,9) and one socalled "impossible figure" (Figure 2a). We will discuss the ability of the algorithm developed by Zou and Hovmöller[9]to assist in determination of correct plane group via so called residuals (Examples 4.1 and 4.2). Also we discussusefulness of the Phase Origin map (POM) to reveal "broken" (Example 4.2) and "hidden" (Example 4.3 and 4.4)symmetry elements, often associated with colorization of plane tessellations.
The essence of Crystallographic image processing (CIP) and its the basic steps outlined in the quote of Sir AaronKlug, a Nobel Prize winner and pioneer of the CIP[10]: "The essence of image processing of this type is that it is atwo-step procedure after the first image has been obtained. First the Fourier transform of the raw image is produced.Next, Fourier coefficients are manipulated, or otherwise corrected, and then transformed back again to reproduce thereconstructed image.” Applicability of the CIP to high resolution TEM (HRTEM) images to facilitate crystal structuredetermination has been extensively developed in the research group of Sven Hovmöller. (See e.g. [9]) A great amountof real live examples is use of the CIP techniques in application to crystal structure determination from HRTEM canbe found in [11]. First application of this technique to images other than HRTEM, namely to images coming from ascanning tunnelling microscope (STM), was done by P.Moeck.[12] Also it had been pointed out that the CIP techniquescan aid in removal of scanning artifacts.[13]Here we make the definition of the CIP agnostic of the method how the images were obtained and whether they depict acrystal structure or not: all images that contain any intrinsic symmetry, including just translational symmetry, can beprocessed "crystalographically". To demostrated this concept alongside with some numerical peculiarities of the CIPmethods we will apply them to Escher’s drawings.This family of the CIP methods is based on analysis of a symmetric pattern in the Fourier (inverse) space. Since theimage periodic, its Fourier space representation is build up by a discrete set of complex-valued Fourier components(FC). The Fourier transformation has direct relation with electron or X-ray diffraction - methods widely used for crystalstructure determination. Electrons (or X-rays) diffract on certain crystallographic planes with Miller indexes ( hk ) , or ( hkl ) in the 3D space, however we limit our scope here to 2D drawings/structures. Therefore, each peak in diffractionpicture is labeled by the indexes of the corresponding Miller plane. Similarly each peak in the reciprocal (or Fourier)space is also assigned a pair ( hk ) indicating its diffraction "order". The FCs are also commonly referred in X-ray andelectron crystallography as "structure factors" or "reflections", here, however, we refrain from using this term since itmostly applies to experimental methods of structure determination.A plane (or space groups) possess a certain number of symmetry elements linking together the atomic positions in thereal space. Correspondingly this symmetry is preserved upon the Fourier transformation, and now the magnitudes andthe phases of certain FCs in the reciprocal space are not random anymore, but have to match with those of another FCconnected to the first one by a symmetry opration. Further details are discussed in section 3.3. The CIP algorithm is based on certain relations between the phases ( φ ( hk ) ) and the amplitudes ( F ( hk ) ◦ or ◦ with their actualphases, and separately compares phases of low order FCs (i.e. for which h and k are small, as defined by the authors ofthe program) with their symmetry related mates. R-factor is then employed as a measure of goodness of fit, definedas R = (cid:80) | φ obs − φ calc | (cid:80) | φ obs | . Here φ obs ( hk ) is the experimentally observed phase, and φ calc ( hk ) - the theoretical phaseobtained from a proposed structural model. Thus, two different R-factors are available in ALLSPACE: "vs. theoretical",i.e. where φ calc is calculated based solely on the plane group information, and "vs. other spots" where the programattempts to link pair of symmetry related and calculate mutual phase discrepancy. Further details about the ALLSPACEalgorithm can be found in the web page. 2 PREPRINT - A
UGUST
22, 2019Another approach proposed by Zou and Hovmöller,[9] is an algorithm for the evaluation of the so called "symmetrizedphase", defined as φ sym ( hk ) = arctan (cid:32) (cid:80) j s j w j sin( φ jobs ( hk )) (cid:80) j s j w j cos( φ jobs ( hk )) (cid:33) + (cid:40) ◦ if (cid:80) j s j w j cos( φ jobs ( hk )) > ◦ if (cid:80) j s j w j cos( φ jobs ( hk )) < (1)The summation in Eq. (1) runs over all symmetry-related FCs defined by the desired symmetry group; w j is a weightingfactor, often set to the amplitude of the corresponding FC, s j = 1 if the phases φ sym ( hk ) and φ jobs ( hk ) should be equalfrom the plane group relations and s j = − if they should differ by ◦ . If a FC is not related to any other FC by thesymmetry (except by the Friedel’s law), then: φ sym ( hk ) = φ obs ( hk ) . Derivation of this formula and a computationallyefficient version of it is given in the Appendix.Optimal phasing is determined as a minimum of the "phase residual" functional φ Res [ φ obs ( hk )] = (cid:88) hk F hk | φ obs ( hk ) − φ sym ( hk ) | (cid:30) (cid:88) hk F hk (2)for each plane symmetry group. Here φ sym ( hk ) is the "symmetrized phase" which fulfills the symmetry relations andrestrictions. All information on plane (wallpaper) groups is summarized in the many-volume International Tables forCrystallography.[14, 15] A reader is cardinally referred to this compendium of crystallographic information for spaceand plane group notation, and a comprehensive list of symmetry operators in each group.Explicit relations between the FCs for each plane group are tabulated in e.g. Refs. [16] and[17]. φ sym ( hk ) is the same for all symmetry related FCs. Centrosymmetric plane symmetry groups( p , p mm, p gm, p gg, c mm, p mm, p gm, p mm ) impose further restrictions on possible values of symmetrizedphases limiting them to ◦ or ◦ only. The lower the residuals φ Res are for a given plane symmetry group, the higheris the likelihood that this group is the right one for the pattern under investigation.Similarly an "Amplitude residual functional" can be introduced: F Res [ F obs ( hk )] = (cid:88) hk | F obs ( hk ) − F sym ( hk ) | (cid:30) (cid:88) hk F obs ( hk ) (3)This functional is identical to the R-work factor, widely used in crystallography and structure determination. However,this functional is less important due to the fact the Fourier transform (FT) contains more information in the phase spacethan in the amplitude space.[18] E.g. if we "tamper" the Fourier transform of an image and modify it such way thatthe amplitudes of all FCs are set the same, but the phases are preserved from the original image, the inverse FT willstill yield a recognizable image. If, in turn, the amplitudes are left intact, but all phases are set the same, then thereconstructed image is completely scrambled because it is the phase that determines which FC are present or not, andthe amplitudes only determine the contribution of each FC.Another quantifier to differentiate between plane groups is the ratio F o /F e calculated as the total amplitude of allsystematically absent FCs (that were nevertheless observed) to the total amplitude of all other FCs.[9] Systematicallyabsent FCs are the FCs whose F e amplitudes vanish due to presence of glide line symmetry elements in the planegroup ( p g , p g , c m , c m , p mg , p gm , p gg , c mm , p gm ).[14] (A glide reflection is the composition of areflection in a line and a translation along that line.)We use full plane group symbols instead of the reduced ones (i.e. p g vs. pg ) to emphasize the orientation of asymmetry element with respect to the image as it is was drawn. An F o /F e ratio is thus another quantifier for "degreeof matching" of the image under inspection to a certain to a plane group in addition to the amplitude and the phaseresiduals. If for an experimental image the F o /F e ratio is larger than zero for a given plane symmetry group, the imagedoes not perfectly fit the group potentially due to the presence of hidden and broken symmetry elements — see below.These formulas are implemented in the CRISP program.[19] CRISP is a closed-source program with no possibilityto interact with the code for further examination of the algorithm, therefore we have developed our own softwareimplementation of the aforementioned scheme in order to be able to explore the symmetry deeper.[20, 21, 22] The image processing procedure begins with a Fast Fourier Transform (FFT) of the input image. Information is collectedexclusively from the Friedel-independent (half-plane) part of the reciprocal space. Since any real image is neither truly3
PREPRINT - A
UGUST
22, 2019periodic nor infinite, edge effects such as peak streaking may occur. This causes intensity of the peaks to decrease andmakes exact detection of their position more challenging. To reduce this effect, the image is multiplied by a windowfunction before the FFT. 2D Hann window, a common windowing technique in the FFT, has successfully eliminatedpeak streaking. Since some plane groups possess systematic absences the amplitudes of some peaks will be zero, whichcan adversely affect the indexing routine later, therefore we add one more auxiliary operation - an auto cross correlationof the power spectrum (squared amplitude part of the FFT). As the result of this routine, all peaks that previouslyvanished in the power spectrum are now revealed in their correct positions. Once all of these corrections are applied,the program performs a peak search on the power spectrum. It should be noted that though the peak search is performedon the auto-correlated power spectrum, intensities of found peaks are collected from the original power spectrum.
As the result of the peak search the program obtain number pairs { x i ; y i } , called "image coordinates" in pixels, whichdepend on the size of the image (in pixels). Corresponding pairs in the Fourier space { h i ; k i } , are referred to as"reciprocal lattice coordinates", and depend on the basis vectors assigned to the reciprocal lattice. Within a known basisset ( (cid:126)U , (cid:126)V ) image coordinates and lattice coordinates are mutually linked by the relation: (cid:18) h i k i (cid:19) = (cid:18) U V U V (cid:19) − (cid:18) x i y i (cid:19) (4)Here U , U are the components of the (cid:126)U basis vector and V , V are the components of (cid:126)V , both expressed in pixels.This approach follows somewhat the one described in Ref. [23], however this time several more refinement stepswere introduced several more to determine the basis vectors. The choice of basis is done by screening through thelow-resolution peaks with weighting factors applied to make rectangular/hexagonal lattices more preferable rather thanoblique. The basis, which indexes a sufficient number of low-resolution peaks, is further refined by minimizing thevalue of functional L using a least squares algorithm (LS). L (cid:104) (cid:126)U , (cid:126)V (cid:105) = (cid:88) i (cid:12)(cid:12)(cid:12)(cid:12)(cid:18) U V U V (cid:19) (cid:20) h i k i (cid:21) − (cid:18) x i y i (cid:19)(cid:12)(cid:12)(cid:12)(cid:12) ,where []-represents rounding to the nearest integer. (5)As long as the basis is refined and all peaks are indexed, the program can compute "lattice parameters" ( a ∗ , b ∗ ) as a ∗ = (cid:112) U + U , b ∗ = (cid:112) V + V , a ∗ b ∗ cos γ ∗ = U V + U V . These values can be further refined via a similarLS routine, as a minimum of the functional Q : Q [ a ∗ , b ∗ ] = (cid:80) i (cid:0)(cid:0) x i + y i (cid:1) − (cid:0) h i a ∗ + k i b ∗ + 2 h i k i a ∗ b ∗ cos γ ∗ (cid:1)(cid:1) . The value of γ ∗ is not refined, since in allgroups except p and p , the indexing of which is trivial, it is constrained either to ◦ or to ◦ . The "latticeparameters" are arbitrary numbers, which depend on the processed image’s length calibration and resolution. Thereforeexact numbers do not have any real physical sense. However if the ratio between a ∗ /b ∗ approaches unity and the anglebetween both vectors is significantly close to ◦ or ◦ , the lattice is likely either hexagonal or square respectively.The program stores lattice coordinates and FCs ( h, k, F, φ ) for each indexed peak for further processing. Related FCs are fully determined by the plain group’s symmetry operators. In vol. B of the "International tables" (SeeRef. [15], Table A1.4.4.1. Crystallographic space groups in reciprocal space, p150) all 3D relations are listed. Bydropping the l index in the space groups corresponding to the appropriate plain group (e.g. P ba −→ p gg ) one canobtain phase relations in the plane groups. Friedel-pairs should be dropped too, since they do not carry any additionalinformation. As one can see from table A1.4.4.1 in Ref. [15], certain groups, like p m , c m , p mm , c mm havethe same phase relations. Thus, groups p m / c m and p mm / c mm can be distinguished by systematical absencesdue to centering only. Systematical absences contribute to φ res due to the weighting factors F hk . In the p group thereare no relations, restrictions only, i.e. each phase should be restricted to ◦ / ◦ . Once all related pairs are determined,Eq. (1) is applied to calculate a symmetrized phase for each group of related FCs. An important feature first implemented in CRISP is the "origin search and refine procedure". It utilizes the phase originmap (POM), which is being built by plotting the values of φ Res for the user-selected plane group as a function of φ x and φ y , related to the arguments of the φ Res [ φ obs ] functional as (See Ref. [15], section 1.4.2.5): φ (cid:48) obs ( hk ) = φ obs ( hk ) + ( φ x h + φ y k ) (6)4 PREPRINT - A
UGUST
22, 2019The POM is computed within the range { φ x , φ y } = − ◦ ... ◦ . Then CRISP looks for a minimum on the POMand shifts the phase of each FC ( hk ) by ( φ min x h + φ min y k ) . Usability of the POM goes beyond a simple minimumsearch. Cell metrics can be obtained from the POM as well: its apexes and centering may be recognized by strong,well-localized minima. The symmetry of the POM should be in an ideal case the same as the symmetry of the patternanalyzed. Certain features of the POM will also represent hidden or broken symmetry elements. (See below in theexamples section) Groups with glide lines
Some groups ( p g , p g , p mg , p gm , p gg , p gm ) exhibit non-trivial phase relations,involving phase shift by a "glide factor" k × ◦ ( p g , p mg ) or ( h + k ) × ◦ ( p gg , p gm ). These factorsmake estimation of the phase error using just functional (2) and formula (1) inapplicable for groups with glide lines.Therefore, we redefine s hk in formula (1) as follows: s hk = (cid:26) − if "glide factor" is odd +1 if "glide factor" is even (cid:18) since sincos ( φ + k × ◦ ) = ( − k × sincos ( φ ) (cid:19) (7)Similar correction factors should be applied to functional (2) by multiplying φ sym by the s -factor: φ Res [ φ obs ( hk )] = (cid:88) hk F hk | φ obs ( hk ) − s hk φ sym ( hk ) | (cid:30) (cid:88) hk F hk (8)These corrections guarantee an artifact-free phase origin map with minima corresponding exactly to their real positionas determined by the plane group. After shifting phases to these minima as required by Eq. (6) the reconstructed imageis expected to be free from symmetry artifacts too. "Symmetry enforcement" means the following: For a given plane group, only the Fourier components from thefundamental domain are selected and then the amplitudes and phases from this domain are multiplied by thecorresponding symmetry operators onto the whole plane. From this an image is synthesized as a Fourier sum: I ( x, y ) = (cid:80) h,k F hk exp( iφ hk ) exp( i ( hx + ky )) . The resultant image is called "enforced" in a certain group."Symmetry Refinement" means: A new set of FCs is built, based upon the input data but with symmetry-related FCsset to their respective mean values: phases restricted to either the symmetrized phase or to the nearest ◦ or +180 ◦ incentrosymmetric groups, phases are shifted according to the POM minimum, or any desired offset. This new set ofcoefficients is then used to build a new real-space image using the Fourier sum. Fourier sum is a continuous function,that can be used to synthesize image of any pixel size. Enforcement and/or refinement can be performed for any of theplane groups, regardless of the residuals. The residuals are just a rank of the plane groups in terms of their likelihood ofmatching the original symmetry. Enforced and refined images should be compared with the original ones to ensure thatenforcement and refinement was conducted in the correct plane group. As an example of automatic symmetry detection using phase residuals we will consider Escher’s woodcut "Angel-Devil". The pattern as follows from the overlaid unit cell (Figure 1a), belongs to the p mm group. This imagedemonstrates only one concept: the ability of the program to detect the most probable symmetry by analyzing theresiduals. Hereinafter "probable group" is assumed in the meaning of the groups having the lowest amplitude and phaseresiduals. There are other criteria that can go hand in hand with Fourier methods, e.g. Akaike criteria.[24] This image ischosen because angels and devils are considerably different geometrically and thus no "hidden" symmetry elements canconnect them. After performing FFT of the original image (Figure 1a), auto cross correlation, and peak search Figure1b is obtained. The unit cell in the reciprocal space is marked by the blue and green vectors (online only); green is thex-axis translation vector and blue is the y-axis.Only the Friedel independent part of the space is used for the peak search. The size of the red circles represents theintensity of the actual Fourier transformation, whereas the white dots correspond to the auto cross correlated image.(Note how FC the previously invisible due to the g -glide line are now revealed). This pattern can be indexed in a squarecell with lattice parameters a ∗ = 4 . , b ∗ = 3 . , γ ∗ = 90 ◦ . Slight deviation from the exact relation a ∗ = b ∗ isdue to accumulation of numerical errors and limited number of repeating elements. Green circles indicate that a peak isindexed in the given basis vector set. The peak search procedure had collected 8264 FCs, within the range of h max < ,5 PREPRINT - A
UGUST
22, 2019 | k max | < . On the next step the program searches related FCs, using the relations as extracted from the table inRef. [15] and computes residuals for a user-selected plane group. (a) (b) Figure 1: (a) "Angel-Devil" (No.45), M.C.Escher, 1941, with a p gm unit cell superimposed. (b) Indexed FFT powerspectrum of Figure 1a (central part)Amplitude and phase residuals for different plane groups are summarized in Table 1.Table 1: Amplitude and Phase residuals for Figure 1a Plane Group Amplitude Residuals Phase Residuals F o /F e p2 15.0 24.2 -p1m1 15.7 31.0 -p11m 15.7 31.0 - p g p g p gg p p gm p gm group as well anyin any other italicized group from Table 1 we observe that the original image doesn’t change. Indeed, all italicizedgroups are the subgroups of p gm . The generating group ( p gm ) has slightly higher residuals than all its subgroups.This is because the supergroup also has a higher number of FCs mutually connected by symmetry operators. It is clearlyseen that groups like p mm and its subgroups have the residuals higher.Therefore, formally speaking, every italicized group is a group of symmetry of the image in Figure 1a, but p gm hasthe highest symmetry of all possible groups. Therefore, from the analysis of the residuals and refined images we canconclude about the plane symmetry of the image under inspection. The image refined in the p gm group has less noiseand image processing artifacts due to averaging multiple FCs by symmetry operators.The image in Figure 1a has two elements of motif that are clearly distinct from each other. Therefore it does not containany "hidden" or "broken" symmetry elements. Now we will consider an image with such elements present.6 PREPRINT - A
UGUST
22, 2019
In this paragraph we will consider an Escher’s inspired hexagonal tessellation on the plane. This is one of socalled "Impossible figures".[26] The pattern as follows from Figure 2a, belongs to the p m group. This image isconsiderably more complex due to presence of elements closely resembling each other but yet not connected by theexact ("obvious") symmetry elements. This pattern can be indexed in slightly distorted hexagonal cell with latticeparameters a ∗ = 14 . , b ∗ = 14 . , γ ∗ = 61 . ◦ . Green circles indicate that a peak is indexed in the given basisvector set. The diameter of each red circle indicates initial, not enhanced, amplitude. The peak search procedure hadcollected 1680 peaks, within the range | h max | < , | k max | < (a) (b) Figure 2: (a)"Impossible figure" a hexagonal pattern with a p m group unit cell superimposed. (b) FFT powerspectrum of Figure 2a. FC from the lower semi-plane only are shownOn the next step the program searches for related FCs, using the relations extracted from the table in Ref. [15] andcomputes residuals for all trigonal/hexagonal ( p , p m , p m , p , p mm ) and p groups.The phase and amplitude residuals are listed in the table in Figure 3. The p m group has significantly lower phaseresiduals than all of the other groups. Another group with low phase residuals, p is a subgroup of p m . The reasongroup p has slightly higher residuals, possessing exactly the same symmetry is that there are more symmetry relationstaken into account by formula (1). The image, synthesized using symmetrized phases and p m group symmetryconstrains repeats the input image.The picture in the right side of Figure 3 is the Phase Origin Map (POM). It shows φ Res as a function of { φ x , φ y } . Thevalues of (cid:8) φ min x , φ min y (cid:9) corresponding to the lowest φ Res can be converted into the position of unit cell origin. In Figure3 converted position of the minimum is indicated by a circle. A strong line in the POM corresponds to the m -line,whereas two dimmer lines intersecting the strong line are, as follows from the unit cell geometry, the g -lines. Thismap shows mutual arrangement of the symmetry elements within the unit cell. The POM with symmetry elements,superimposed on top of it is shown in Figure 2a.Black lines in Figure 4 correspond to a "might-have-been" symmetry element – a vertical mirror line. Only a fewfeatures in Figure 2a violate this symmetry element, making it "broken" - this is exactly what makes this image"impossible". The "broken" symmetry group would be p m . This group also has relatively small residuals. Groupswith six fold axes exhibit significantly higher residuals, thus they are not "broken", but rather incorrect. Dark lines inthe POM correspond to the rectangular setting of the trigonal lattice.We would like to emphasize importance of the POM for further crystallographic analysis: By careful reviewing ofthe POM one can reveal broken or missing symmetry elements. Or symmetry elements which reduce the possible7 PREPRINT - A
UGUST
22, 2019Figure 3: (Program Screen shot) Figure 2a refined in the p m group. Positions of the 3-fold axes on the POM aremarked with arrows Figure 4: POM for the p m group(apparent) high symmetry to lower by exclusion of some symmetry operators due to coloring, defects, twinning, grainboundaries, magnetization, etc. This example is devoted to the problem of detection and quantification of black and white symmetry by means of theCIP. The concept of black and white symmetry expands the traditional concept of symmetry by adding a color inversionoperator (’) to other purely geometrical operators.[3] Black-and-white groups possess a grey group (i.e. a group whereblack and white are averaged to grey) as their super group. Strictly speaking, grey groups are not defined for patterns8
PREPRINT - A
UGUST
22, 2019with more than two colors; however, here we will use the concept of "colorblind" images, i.e. images where all thecolors are ignored and boundaries are enhanced using the Sobel Filter.As an example of the ability of the algorithm to distinguish the black and white symmetry, we will consider Escher’s"Lizards" (Figure 5). By inspection we can determine only one pure geometrical symmetry operator acting onthis image: glide lines (dashed red) running though the "cheeks" of the white and the black lizards. Since no othergeometrical elements appear, the plane group seems to be pg y ( p g ) polar.Figure 5: "Lizards" (No.124), M. C. Escher, 1965. g y -lines are shown in red, g (cid:48) x - in blue, the unit cell of the b/w group p (cid:48) g (cid:48) g is shown as well.Initial image processing is done the same way as before. Because of use of the autocorrelation during the peak searchroutine, even peaks corresponding to FCs forbidden by the g-line boost their intensity (e.g. (0,1)). This significantly aidsthe indexing procedure. The pattern can be indexed in the rectangular system with the lattice parameters a ∗ = 2 . , b ∗ = 3 . .On the next step we reconstruct the image and compute phase and amplitude residuals. Since the lizards have very sharpboundaries, FFT produces high number of Fourier components, only a fraction (562) of which was collected. Group p g ( pg y ) has phase residuals significantly lower than other groups and the reconstructed image perfectly matches theoriginal one. Note the much higher residuals for plane groups p and p gg . The picture in the right side of Figure 6 isagain the POM. Since group p g is polar, i.e. has no fixed origin along the y direction, the phase residuals are almostthe same along the blue lines. In reality, they differ a little bit due to presence of black-and-white symmetry. Carefulexamination of the POM in Figure 6 reveals additional highly symmetric features such as yellow and dark red lines,and sharp yellow points. As it was shown before this suggests presence of "hidden" symmetry elements.In this case, the "hidden elements" are the elements that change the color of the lizards along with some geometricalmanipulations on them. Indeed if we refine in a super group of p g , namely p gg we obtain a "gray-averaged" image.(Figure 7a). I.e. plane group p gg is the grey super group for the lizards pattern. In Figure 7a one can still, recognizethe lizards’ pattern. If black-white symmetry is taken into consideration, one immediately identifies (cid:48) axis and g (cid:48) x glidelines (depicted in blue in Figure 5). Combination of these two symmetry operators with a gray g y glide line yields theShubnikov (black and white) group of the image in Figure 5: p (cid:48) gg (cid:48) .We can further investigate symmetry of this image by applying the Sobel filter to Figure 5. The Sobel filter enhancesboundaries and effectively eliminates uniform features (such as colorization), thus yielding a "silhouette" image, wherethe colors are ignored (See Figure 7b). In this case the Fedorov (purely geometrical) plane group p gg is not "broken"anymore and symmetry refinement in this group yields a sharp image. (Figure 8)Again, the values of { φ x , φ y } corresponding to the lowest φ Res (Figure 6 blue lines, Figure 8 green and yellow lines)can be converted into the position of the unit cell origin(s). In Figure 7a converted positions of these minima are plottedas red rings. They hull the fundamental domain of the pattern. The minima in the POM represent arrangement of thesymmetry elements within the cell. In Figure 8 vertical lines are more pronounced than the horizontal, which indicatesthat the g y glide-line is less broken than g x , the reason for that is the situation with black and white lizards’ eyes (see9 PREPRINT - A
UGUST
22, 2019Figure 6: (Program Screen shot) Figure 5 refined in the p g group (a) (b) Figure 7: (a) "Gray image" refined in p gg group. (b) Lizards silhouettes (colorblind image)Figure 7b). While there is a slight difference in the lizards’ eyes in Figure 7b, all eyes are averaged out after refinementin the p gg group. See Figure 8. Since only a minor feature of the picture violates the symmetry of the colorblind p gg group the residuals for p gg are now closer to the "proper" symmetry group – p g . Therefore, comparing those two,we can estimate "how much" symmetry is broken by comparing the residuals between "proper" and colorblind groups. A slightly more complex example is given by the "Reptiles" image (Figure 9). Unlike "Lizards", "Reptiles" have threedistinct shades of gray. Therefore, we cannot define the grey group for them anymore, however, we can still obtain acolorblind image for symmetry estimation. The hexagonal pattern here is easy to recognize by a naked eye. The unit10
PREPRINT - A
UGUST
22, 2019Figure 8: (Program Screen shot) Figure 7b refined in p gg cell overlaid with the image depicts all symmetry features of the picture. Thus, by inspection, the picture possesses p colorblind symmetry and p "proper" symmetry, if colors are not ignored. Therefore the full 3-colored (Belov) planegroup for this image is p | p . We follow the notations in Ref. [4] for three-color plane groups G col = G | G (cid:48) , where G is the Fedorov plane group and the subgroup G (cid:48) ⊃ G of index 3 contains operations that keep the first color fixed.Since both Shubnikov (b/w) and Belov (color) plane groups can be represented in terms of the layer groups, where shiftalong the z axis corresponds to color change, the isomorphic 3D layer group for the lizard pattern is P .Figure 9: "Reptiles", M. C. Escher, 1943, Fragment, with a p | p unit cell superimposed. Symmetry elements depictedin orange are of the proper symmetry, in blue – color-blind.Initial image processing was done in the same way as before. While the first order peaks in the FFT power spectrumstill exhibit a hexagonal arrangement, the hexagonal pattern yields to a 2-fold symmetry for the second order FC. (SeeFigure 10a) 11 PREPRINT - A
UGUST
22, 2019On the next step we reconstruct the image and compute phase and amplitude residuals. Groups p ( φ Res = 20 . ), p ( φ Res = 18 . ), p ( φ Res = 22 . ) have phase residuals significantly lower than other groups with consistent cellmetrics: p m ( φ Res = 27 . ), p m ( φ Res = 25 . ), p mm ( φ Res = 39 . ). Only refinement in the p group givesthe correct reconstructed image, therefore we conclude that the Fedorov group for this image is indeed just p . However,the fact that the residuals for p and p are very close to p , it suggests the presence of "hidden" symmetry, which inthis case is color symmetry.In the POM for p the minima, indicated by the bright spots, correspond non broken 2-axes, and the dark dots correspondto the hidden symmetry elements: 6- and 3-axes. The higher is the order of a symmetry element the more relationsbetween FCs it implies. Therefore 6-axes will appear in the POM darker (i.e. higher residuals) than e.g. 2-axes, sincefor an origin shift to a six-fold axis a pair of { φ x , φ y } should obey more restrictions. Thus we conclude that the darkdots on the POM in Figure 10b are most likely 6-axes. This provides an insight onto the highest symmetry grouppossible, which is now hidden due to colorization. By applying the Sobel filter and extracting the silhouettes of thelizards, we can see that the image refined in the p group now clearly resembles the initial one. (a) (b) Figure 10: (a) Indexed Fourier Transformation of Figure 9. (b) POM for p refinement overlaid with p unit cell (thecell is slightly shifted to demonstrate features on the POM)From Figure 11 one can see that indeed, after colorization is removed, nothing else breaks the highest possiblesymmetry: p , and the residuals for groups p and p (subgroup of p ) are lower than they used to be in the coloredcase.An important factor for practical application of the aforementioned techniques of symmetry quantification and refinementis image resolution. If the number of collected FC is quite small (i.e. image resolution is too low), some FC have nomates, and some phase comparisons cannot be performed fully. Usually for high symmetry groups, where the numberof related FCs is more than two, only 10-20% of all FCs have all required mates. Some comparisons may be totallydegenerated; in this case there is no difference in the phase statistics between some high- and low- symmetry groups.Therefore, a supergroup will normally have residuals higher than its subgroups. This is indeed the case in this examplewhere the p group yields lower residuals than the correct p group. We analyzed the applicability of CIP to symmetry analysis of highly symmetric 2D images. Phase residuals allow us todeduce the most probable plane group, low phase residuals in the groups that don’t yield exactly same reconstructedimage after refinement most likely reflect the fact that there are some "broken" or "hidden" symmetry elements. In thiscase a careful analysis of the phase origin map can often reveal possible candidates for both geometrical and black and12
PREPRINT - A
UGUST
22, 2019Figure 11: Symmetrization of a uniformed silhouette image of Reptileswhite supergroups. An edge detection algorithm, like the Sobel filter can provide some insight about the color-blindsupergroup.The four examples discussed above illustrate three important features of the CIP for symmetry determination: simpleidentification of the most probable plane group, detection of "broken" symmetry elements on an example of amonochrome image, detection of "hidden" and "broken" symmetry elements associated with colorization of the initialpattern. 13
PREPRINT - A
UGUST
22, 2019
Appendix: Derivation of a numerically efficient version of Eq. (1)
Original paper by Zou and Hovmöller[9] does not contain derivation of Eq. (1). However, it can be obtained byconsidering strict geometrical relations in the 4-D ( h, k, F, φ ) Fourier space. Suppose we have a set of symmetry relatedFCs as determined by a given plane group: (1,1), (-1,1), (-1,-1), (1,-1). By the rules of symmetry, their amplitudes andphases are related as well. The amplitudes of all FCs under consideration should be equal, and phase increases by ◦ as we go in the counterclockwise direction from the (1,1) FC (See Figure 12a). (0,0) (1,1)(1,-1)(1,1)(-1,-1) (a) F sym (cid:1) sym (cid:1) (cid:2) (cid:1) (cid:3) (cid:1) (cid:4) (cid:1) (cid:5) F F F F ReIm ReIm (b)
Figure 12: (a) Schematic representation of a set of symmetry related FCs. (b) Graphical representation of thesymmetrization procedureThus, their complex-valued FCs Φ hk = F hk exp( iφ hk ) are mutually linked as: Φ = Φ (1 , = F · e iφ Φ = Φ ( − , = F · e i ( φ + π/ = i Φ Φ = Φ ( − , − = F · e i ( φ +2 π/ = − Φ Φ = Φ (1 , − = F · e i ( φ +3 π/ = − i Φ (9)The fore factors (1 , i, − , − i ) appear due to the symmetry and do not bear any additional structural information.However, if the symmetry is slightly broken, i.e. each amplitude and phase slightly deviates from its ideal value, theequations above transform as: Φ = Φ (1 , = F · e iφ Φ = Φ ( − , = F · e i ( φ + π/ = iF · e iφ Φ = Φ ( − , − = F · e i ( φ +2 π/ = − F · e iφ Φ = Φ (1 , − = F · e i ( φ +3 π/ = − iF · e iφ (10)The amplitudes ( F , F , F , F ) and phases ( φ , φ , φ , φ ) do not differ significantly, thus, ignoring the symmetrycaused fore factors, we can represent the complex-valued FC as a bundle of vectors of almost the same length andpointing in almost the same direction: Figure 12b.The symmetrization procedure replaces this bundle by a single complex number (vector), defined as the complex averageof the bundle: Φ sym = (cid:80) ni Φ i (cid:14) n . The following equations give the amplitude and the phase of the symmetrizedstructural factor: F sym = (cid:80) ni | F i ( h, k ) | n ; φ sym = arctan (cid:18) (cid:80) ni F i sin φ i ( h, k ) (cid:80) ni F i cos φ i ( h, k ) (cid:19) (11)These formulas coincide with Eq. (1), with the difference, that the weighting factor w i is now explicitly expressed interms of the amplitude. Also, the origin of the s factor is now explicitly revealed as the presence of glide lines changesthe fore-factor in Eqs. (9), (10) by rotating the corresponding vector either ◦ or ◦ degrees, and thus it has to berotated back for correct averaging. Complex FCs used for image reconstruction are now given as: Φ = Φ sym ; Φ = − Φ sym ; Φ = i Φ sym ; Φ = − i Φ sym ; (12)From the numerical point of view, formula (1) is not very efficient, since it requires serial calculation of n + 1 trigonometric functions, which is a time-consuming routine. As it is shown here, Eq. (1) computes the argument of a sum14 PREPRINT - A
UGUST
22, 2019of many complex numbers, it is worth keeping coefficients of the Fourier transformation for the image in the complex-valued form (i.e. not to separate phase and amplitude). Then the Eq. (1) can be rewritten as: φ sym = arg ( (cid:80) i s i Φ i ) ,which uses only one trigonometric function instead. Using the latter formula instead of (1) yields two to four times(depends on the number of relations) boost in the performance.15 PREPRINT - A
UGUST
22, 2019
References [1] D. Schwarzenbach.
Crystallography . John Wiley, 1996.[2] R.O. Gould and W. Massa.
Crystal Structure Determination . Springer Berlin Heidelberg, 2013.[3] A.V. Shubnikov and V.A. Koptsik.
Symmetry in Science and Art . Plenum Press, 1974.[4] A.V. Shubnikov and N.V. Belov.
Colored Symmetry . Macmillan, 1964.[5] Allen Nussbaum. The mystery of the fifteenth bravais lattice.
American Journal of Physics , 68(10):950–954,2000.[6] J. Amezcua-López and A. E. Cordero-Borboa. A mechanical lattice aid for crystallography teaching.
The PhysicsTeacher , 26(6):384–389, 1988.[7] Klaus Landwehr. Visual discrimination of the 17 plane symmetry groups.
Symmetry , 3(2):207 – 219, 01 2011.[8] L. Koren.
Wabi-sabi for Artists, Designers, Poets & Philosophers . Wabi-sabi for Artists, Designers, Poets &Philosophers. Imperfect Publishing, 2008.[9] Xiaodong Zou and Sven Hovmöller.
Structure Determination from HREM by Crystallographic Image Processing ,pages 275–300. Springer Netherlands, Dordrecht, 2006.[10] Aaron Klug. From virus structure to chromatin: X-ray diffraction to three-dimensional electron microscopy.
Annual Review of Biochemistry , 79(1):1–35, 2010. PMID: 20192760.[11] Xiaodong Zou, Sven Hovmöller, and Peter Oleynikov. Electron crystallography: Electron microscopy and electrondiffraction.
Electron Crystallography: Electron Microscopy and Electron Diffraction , pages 1–344, 01 2012.[12] P. Moeck. Systems for assessing and enhancing the performance of scanning probe microscopes by quantifyingand enforcing symmetries and periodicities in two dimensions, jun " 5" 2012. US Patent 8,196,216.[13] Jack C. Straton, Bill Moon, Taylor T. Bilyeu, and Peter Moeck. Removal of multiple-tip artifacts from scanningtunneling microscope images by crystallographic averaging.
Advanced Structural and Chemical Imaging , 1(1):14,Nov 2015.[14] Theo Hahn, editor.
International Tables for Crystallography, Volume A: Space Group Symmetry . InternationalTables for Crystallography. Kluwer Academic Publishers, Dordrecht, Boston, London, 5. revised edition edition,2002.[15] U. Shmueli, editor.
International Tables for Crystallography, Volume B: Reciprocal Space . International Tablesfor Crystallography. Springer Netherlands, 2008.[16] X. Zou.
Electron Crystallography of Inorganic Structures: Theory and Practice . Chemical Communications.Stockholm University, 1997.[17] Bart Verberck. Symmetry-adapted fourier series for the wallpaper groups.
Symmetry , 4(3):379–426, 2012.[18] A. V. Oppenheim and J. S. Lim. The importance of phase in signals.
Proceedings of the IEEE , 69(5):529–541,May 1981.[19] Sven Hovmöller. Crisp: crystallographic image processing on a personal computer.
Ultramicroscopy , 41(1):121 –135, 1992.[20] P Moeck, J Straton, B Moon, P Plachinda, M Hietschold, H Glowatzki, N Koch, and JP Rabe. Crystallographicprocessing of scanning probe microscopy images on the example of 2d periodic molecular monolayer arrays.
Microscopy and Microanalysis , 16(S2):456–457, 2010.[21] P. Moeck, P. Plachinda, B. Moon, J. Straton, S. Rouvimov, M. Toader, M. Abdel-Hafiez, and M. Hietschold.Quantifying and enforcing two-dimensional symmetries in scanning probe microscopy images of periodic objects.In , pages 1–2, Dec 2009.[22] Pavel Plachinda, Bill Moon, and Peter Moeck. Crystallographic image processing software for scanning probemicroscopists. APS March Meeting, 2010.[23] Xiangyan Zeng, Bryant Gipson, Zi Yan Zheng, Ludovic Renault, and Henning Stahlberg. Automatic latticedetermination for two-dimensional crystal images.
Journal of Structural Biology , 160(3):353 – 361, 2007. ElectronCrystallography of Membrane Proteins.[24] Peter Moeck. Towards generalized noise-level dependent crystallographic symmetry classifications of more orless periodic crystal patterns.
Symmetry , 10(5), 2018.[25] V. Gramlich and H. Wondratschek.
Graphs for klassengleiche subgroups , pages 415–425. Springer Netherlands,Dordrecht, 2004.[26] Bruno Ernst and Hans de Rijk.