Image Anal Stereol 2005;24:135-143 Original Research Paper DECIPHERING THE FINEST IMPRINT OF GLACIAL EROSION: OBJECTIVE ANALYSIS OF STRIAE PATTERNS ON BEDROCK Piet Stroeven1, Arjen Peter Stroeven2, Jing Hu1, Jean-Louis Chermant3 and Michel Coster3 1Faculty of Civil Engineering & Geosciences, Delft University of Technology, Delft, The Netherlands; 2Department of Physical Geography and Quaternary Geology, Stockholm University, Stockholm, Sweden; 3LARMAUR, FRE CNRS 2717, Université de Rennes 1, Rennes, France e-mail: P.Stroeven@ct.tudelft.nl; arjen.stroeven@natgeo.su.se; J.Hu@citg.tudelft.nl; Jean-louis.chermant@univ-rennes1.fr; Michelcoster@aol.com (Accepted September 19, 2005) ABSTRACT The aim of this study is to compare the efficiency of different mathematical and statistical geometrical methods applied to characterise the orientation distribution of striae on bedrock for deciphering the finest imprint of glacial erosion. The involved methods include automatic image analysis techniques of Fast Fourier Transform (FFT), and the experimental investigations by means of Saltikov’s directed secants analysis (rose of intersection densities), applied to digital and analogue images of the striae pattern, respectively. In addition, the experimental data were compared with the modelling results made on the basis of Underwood’s concept of linear systems in a plane. The experimental and modelling approaches in the framework of stereology yield consistent results. These results reveal that stereological methods allow a reliable and efficient delineation of different families of glacial striae from a complex record imprinted in bedrock. Keywords: Fast Fourier Transform, orientation distribution, rose of intersection densities, striae patterns, Underwood concept. INTRODUCTION Glaciers and ice sheets are natural phenomena that have shaped the landscape of northern Hemisphere continents to a considerable degree (Boulton and Clark, 1990; Kleman et al., 1997; 2005). Glaciers on the one hand are elongated bodies of ice, which typically are 10 to 100 km long; originate in high mountain valleys and of which ice flow direction follows the gradient of the underlying topography (Stroeven, 2004). Today there are 10,000’s of glaciers situated in mountain regions all across the world. Ice sheets, on the other hand, are dome-shaped bodies of ice, which can be several kilometres thick, and 1000’s of kilometres across. These ice sheets, of which there currently are only two, on Greenland and Antarctica, exhibit ice flow directions which are dictated by the ice sheet surface slope (ice flowing from ice spreading centres, ice domes, to their margins following the steepest surface gradient). For these large continent-wide bodies of ice, an ice flow direction at a certain location, therefore, signifies a particular ice surface slope direction which relates to large-scale ice sheet surface configurations of domes, saddles and margins (Kleman and Borgström, 1996). During the last ice age, at the height of which, around 21,000 calendar years ago, there were large ice sheets centred on north America (Laurentide ice sheet, Cordilleran ice sheet), Scandinavia (Fennoscandian ice sheet) and northern Russia (Eurasian ice sheet), a considerable portion of the northern Hemisphere landmass was therefore covered by ice (Denton and Hughes, 1981; Ehlers and Gibbard, 2004). These former ice sheets left a legacy of their existence as they changed the landscape underneath (Sugden and John, 1976). To do so, ice sheets need to be warm-based, that is, the pressure underneath the ice body must be large enough to allow for melting to occur at the given ambient temperature (Paterson, 1994). During warm-based conditions the ice is allowed to slide along its substrate (often facilitated by the presence of a film of water when the substrate is impermeable) or to deform its substrate, conditions which can induce the formation of ice flow-parallel sediment features (landforms). If the ice pressure drops (ice sheet thinning) or ambient ice temperature decreases, conditions may fall below pressure melting conditions and water will refreeze. When the pressure conditions 135 Stroeven P et al: Stereological analysis of striae patterns on bedrock at the ice sheet-substratum interface are such that they are below pressure melting, the ice sheet is said to be cold-based and frozen to its substrate, yielding conditions inapt for either erosion or deposition (sedimentation). When ice sheets are cold-based, they may advance and retreat across an area without the generation of landforms through erosional or depositional processes, save for the imprints generated by meltwater (Hättestrand, 1998). Given the potential paucity of ice flow directional data for some regions formerly covered by cold-based ice sheets, it is important that evidence which is uncovered can be analysed to a high degree of precision and reliability. An area where there is a paucity of glacial geomorphologic data is from northeastern Arctic Canada, presumably because the former Innuitian ice sheet was largely cold based across the islands that form the archipelago. However, Hättestrand and Stroeven (1996) documented the presence of striae, ice sheet erosional imprints, on bedrock outcrops across Baillie-Hamilton Island (BHI), strongly supporting the notion of an extensive and thick ice sheet overrunning the island from different directions at different times. Striae are scratches on bedrock surfaces come about by the grinding of rocks held in the basal layer of the overlying ice. Even at the spatial scale considered, of suites of cross-cutting striae on a single bedrock outcrop (Fig. 1), and although the length of time required to form “individual” striae may be in the 10-1-101 years range, the length of time represented by the cross-cutting striae record remains enigmatic and could relate to small-scale changes in ice flow direction (and reflect changes that have come about in 10-102 years), large-scale changes in ice sheet topology (shifts in the location of ice-spreading centres, i.e,. domes, saddles and ice margins; taking 103-104 years) or result from ice flow during different glaciations that occurred on timescales of 104-105 years (Glasser and Bennett, 2004). Hättestrand and Stroeven (1996) inferred that the striae record reflected a succession of different ice flow directions forced by a massive migration of an ice dome from a central full-glacial location to the north towards a late glacial position to the west. Measurements of ice flow direction on a striated bedrock surface are normally accomplished using a compass. We here want to investigate whether field studies of striae can be accompanied by additional measurements (and verification) on photographs of the striated bedrock surfaces in a laboratory or computerized setting. To that end, a bedrock surface photograph from BHI (Fig. 1) was analysed by means of Fast Fourier Transform (FFT) and directed secants techniques. METHODS Measurements on a sample of the striated bedrock on BHI were used to illustrate an objective and efficient methodology for an assessment of the orientation direction of the (supposed) three families of striae. Not all striae in one family run perfectly parallel due to small-scale changes in ice sheet configuration, as argued earlier. However, as a criterion for the visual determination of striae direction in the field by compass, the angular differences between the three sets of striae are larger than the within-set variability in orientation direction. The interest therefore is in determining the three main directions in an objective and efficient way. To determine an orientation direction, automatic or manual methods can be used (Underwood, 1970; Serra, 1982; Coster and Chermant, 1989). The most often used methods are the rose of directions or of intercepts. Using a grid with an automatic image analyzer, there is a bias due to the digitisation of the frame of measurements: hexagonal or square grid. According to the digitization, the error is more or less important for the investigated directions, in spite of the use of a corrective factor which can be calculated from trigonometry (Chaix and Grillon, 1996). An important new, more elegant, and rapid method to determine orientation directions of elements on images is by application of the FFT technique. Fig. 1. Photograph displaying striae patterns on bedrock from Baillie-Hamilton Island, Arctic Canada. Box indicates the area of analysis in this study (cf. Fig. 2a and Fig. 5). 136 Image Anal Stereol 2005;24:135-143 The FFT technique has been successfully applied in various fields for investigating orientation distribution characteristics of linear systems in plane projections or in space, e.g., short fibres in plastics (Gonzalez et al., 1994) and fibres in concretes (Redon et al., 1998a, b). FFT is used to decompose an image (spatial domain) into its sine and cosine components. The output of the transformation represents the image in the Fourier or frequency domain. In the Fourier domain, each point represents a particular frequency contained in the analysed image. Hence, FFT is used to access the geometric characteristics of an image. Because the image in the Fourier domain is decomposed into its sinusoidal components, it is a straightforward procedure to examine or process certain frequencies, thereby influencing the geometric structure in the image. In the present study, we will restrict the discussion to the two-dimensional Fast Fourier Transform (2D FFT) of a 4-connexity lattices scan of the striated bedrock. There are various forms of the 2D FFT and the one used in this study is the 2D Fast Hartley Transform. The mathematical details are described in Gonzales and Woods (1992). As an alternative to the automatic image analysis by means of FFT, Saltikov’s method of directed secants (Saltikov, 1958; 1974) was applied to an analogue image of the rock for the construction of roses of the number of intersections, or of intersection densities after normalizing by the total grid line length. As previously mentioned, it is well known that such roses offer biased information when derived from digitised images, limiting data generation in an automated set up (Stroeven et al., 2001). By separate application to the three families of striae, three roses of intersection densities are obtained that yield quantitative information on the relative strengths of the signals. Application of Underwood’s concept for a combination of three linear systems allows construction of individual roses of intersection densities, and a total rose of intersection densities, for the striae in the original image (Underwood, 1970). Comparison of this idealized theoretical concept with the experimental data gives insight into representativeness of the information provided by the sample, and on the degree of orientation uniformity of the families of striae. RESULTS FAST FOURIER TRANSFORM The original photograph of striae patterns (Fig. 1) was subjected to certain image process procedures, using an automatic image analyzer Matra Pericolor 3100 (MS2I, Saint-Quentin en Yvelines, France) in order to enhance the contrast between the striae and the background, resulting in a digital image (Fig. 2a) highlighting the linear traces, which is more suitable for the FFT analysis. Fig. 2b is the corresponding image of the Fourier domain. The Fourier Transform produces a complex number valued output image which can be displayed with two parts, i.e., magnitude and phase images. The magnitude image allows identifying the structural patterns in the spatial domain image. In the phase image, the value of each point determines the phase of the corresponding frequency. In general, only the magnitude of the Fourier Transform is displayed in image processing, as it contains most of the information of the geometric structure of the spatial domain image. In this study, we will restrict ourselves to displaying only the magnitude of the Fourier Transform (the Fourier image). (a) (b) Fig. 2. (a) Spatial domain (image) highlighting the striae in Fig. 1 box; (b) the corresponding image in the Fourier domain. 137 Stroeven P et al: Stereological analysis of striae patterns on bedrock In most implementations the Fourier image is shifted in such a way that the DC-value (i.e., the image mean) F(0,0) is displayed in the centre of the Fourier image. The further away from the centre an image point is, the higher is its corresponding frequency. The result (Fig. 2b) reveals that the DC-value is by far the largest component of the image. The Fourier image contains components of all frequencies, but that their magnitudes get smaller for higher frequencies. Hence, low frequencies contain more image information than the higher ones. However, because the dynamic range of the Fourier coefficients (based on the intensity values in the Fourier image) is too large to be displayed in Fig. 2b, all other values appear as black. If the measurements are based on a 24-bits colour image (instead of an 8-bit grey level image) and a well adapted Look-Up Table, the visualisation of the Fourier image will be better legible. The Fourier image displays three dominating striae directions, two passing slightly inclined to the vertical axis through the centre, at 74° and 57° to the horizontal axis, and a more horizontal third direction at 138° with the positive horizontal axis. The vertical and horizontal stripes (lines) in the Fourier image originate from the regular patterns in the background of the original image. The third family of striae is less pronounced than the other two dominant families. As an effective alternative, a length-weighted average orientation direction could be easily established manually on a hand-made copy of the original image (will be discussed in the following section). It should be noted that the major directions in the Fourier image are determined by the orientation distribution of striae in the digital image (i.e., the spatial domain image). For example, the most significant direction at 74° with the horizontal axis (Fig. 2b), reveals that the image intensity in the spatial domain changes the most if we analyse the image in this direction. This implies a preferred orientation of the striae in a direction perpendicular to this, i.e., at 164° with the positive horizontal axis. This most dominant family of striae is shown segmented from the rest of the image in Fig. 3. Fig. 4 shows the same operation performed on the second strongest signal, a family of striae at 147° (57+90) with the positive horizontal axis. In conclusion, Figs. 2-4 clearly reveal that the principal orientations of the two major striae families, i.e., a1 = 164° and a2 = 147° differ from each other by c. 17°. This result is quite compatible with two independent field estimates of 18° (relative to local sun angle) and 24° (when determined graphically using topographic features in the field; cf. Fig. 5 in Hättestrand and Stroeven, 1996). Although the signal strength for the third family is less significant, the Fourier domain still permits a rough estimation of the preferred orientation of striae responsible for this weak signal of a3 = 48° (138-90). Visual inspection of Fig. 2a indicates a larger within-set variation in direction in the third family of striae, when compared to the previous two and dominant families of striae. ROSE OF INTERSECTION DENSITIES In view of the afore-mentioned biases for rose of intersection densities induced by digitisation in automatic image analysis, a magnification of the original striae image has been manually copied on a superimposed transparency, using different colours for the three families of striae (Fig. 5). This allows an unambiguous assessment of the number of intersections for each family of striae with a superimposed rotating grid (at angular steps of 15°), because the unique decision on the geometric features of each striae is maintained for all grid directions. Total grid line length varied between 400 mm and 1000 mm on the 21,147 mm2 image. An optimum grid line length is obtained when the numbers of intersections are of the same order of magnitude as the number of striae. With the line lengths indicated, this is only realized for the two most dominant families. The resulting roses of intersection densities are presented in Fig. 6. Because the scale factor between the original image and the transparency copy has not been accounted for, only relative measures will be considered. Note that the intersection density with a grid attains a maximum in a direction perpendicular to the orientation direction of the striae. Hence, by rotating the roses of intersection densities in Fig. 6 with 90°, true orientation distributions of the families of striae and the relative strengths of the signals are obtained. Table 1 shows that the principal direction of the roses of intersection densities of the three major families of striae thus derived mimics the results obtained through the FFT analysis. Table 1 moreover lists the strengths of the individual orientation directions determined by Saltikov’s directed secants analysis (maximum radius lengths), PL1, PL2 and PL3, varying from 0.223, over 0.063 to 0.030 mm-1 (see Fig. 6d). 138 Image Anal Stereol 2005;24:135-143 (a) (b) Fig. 3. (a) Spatial domain image highlighting the most dominant family of striae; (b) Fourier domain image of the most dominant family of the striae. &¦ '.. ¦¦¦-¦;>¦¦¦*¦'.;? '. ¦. '¦ ' ' ' *' ': -¦"¦¦' ili#S §§P'; || (a) (b) Fig. 4. (a) Spatial domain image highlighting the second dominant family of striae; (b) Fourier domain image of the second dominant family of striae. Fig. 5. Manual copy of the three dominant families of striae in Fig. 2a. Table 1. Strength and preferred orientation of striae direction by stereological techniques. PL1 red PL2 green PL3 black Max. radius of rose 0.223 0.063 0.030 Max. orientation at 75 53 135 Preferred degree (angular step of 15°) 165 143 45 Preferred degree obtained by FFT 164 147 48 Preferred degree (angular step of 60°) according to 160 141 60 Serra’s elliptical model (Serra, 1982) 139 Stroeven P et al: Stereological analysis of striae patterns on bedrock (a) (c) (b) (d) Fig. 6. Rose of intersection densities for (a) red, (b) green and (c) black portions of the striae pattern in Fig. 5, and (d) the superimposition of all the three families of striae. On the basis of an angular increment of every 60° (thus measuring number of intersections only at 0°, 60° and 120°), Serra (1982) proposed an elliptical model to estimate the preferred direction of distribution. According to this model, we find the preferred degree of 160°, 141° and 60°, respectively, for the red, green and black portions of the striae patterns. These results, although based on less labour-intensity for counting, yield similar data as in the case of an angular step of 15°, especially for the two most dominating families (Table 1). This implies that, when orientation distribution of linear systems is relatively uniform (as the red and green portions), a larger angular increment for rose of intersections can be adopted, thus reducing the labour intensity of the measurement procedures. In addition, it should be noted that FFT also allows for construction of the rose of directions which gives information on the preferred orientation of distribution comparable to that obtained by the Saltikov’s directed secants methods. The roses obtained by Saltikov’s method and by FFT can be associated via a certain integration equation. For details, please consult Redon et al. (1998b). DISCUSSION UNDERWOOD’S CONCEPT OF LINEAL SYSTEMS The individual roses of intersection densities (Fig. 6a-c) can be conceived as lineal systems that can be described successively by (Underwood, 1970) 140 Image Anal Stereol 2005;24:135-143 PL 1(0) = PL 1 sin(6-a1) = 0.223 sin(0-165) PL2(#) = PL2sin(#-a2) = 0.063sin(#-143), (1) PL3(6) = PL3 sin(0-a3) = 0.030sin(0-45) in which a1, a2 and a3 are the angles enclosed by the positive horizontal axis and the principle orientation axis of the respective systems (Table 1). The respective maximum radii (PL1, PL2 and PL3; Table 1), derived from Saltikov’s directed secants analysis, represent the signal intensities of the three families. The modelled roses for the three portions are shown in Fig. 7a-c. Total rose of intersection densities for three systems of lineal features in a plane is given by (Underwood, 1970) Pl{9) = Pl1sin0-cc1) + Pl2sin0- + PL3sin(#- a2) + a3) (2) The sum of the three linear systems of striae thus obtained using Equations (1) and (2) is plotted in Fig. 7d (indicated by the filled area). This result is contrasted with the straightforward sum of the experimental data (outline in Fig. 7d). It clearly reveals a satisfactory agreement between the experimental investigation and the calculation. (a) (b) (c) (d) Fig. 7. Modelled rose (Eq. 1, following Underwood’s concept of linear systems) of intersection densities for (a) red, (b) green and (c) black portions of the striae pattern in Fig. 5, and (d) comparison between the total rose of intersection densities of the three linear systems (of striae) obtained by adding the experimental measurements (outline) and by the modelled roses (Eq. 2) (filled area). 141 Stroeven P et al: Stereological analysis of striae patterns on bedrock We have pointed out that one family of striae (black subset in Fig. 5) obviously experiences larger deviations in striae direction than the other two subsets. Hence, we introduce a third measure of orientation direction, where the orientation information of the weakest system is manually determined by measuring the orientation of all individual striae of this family and multiplication by striae length. After summation and normalization using the total length of striae in this family, a weighted-average orientation direction of this system is derived of a3 = 59°. If the calculation of the total rose of intersection densities by Underwood’s concept of lineal systems is repeated, and replacing the value of a3 by 59°, an even smaller deviation from the experimental curve is obtained. The deviation between the two total roses can be quantitatively expressed by: ö = 2Y\PLM) -PlexpM)\2 , (3) where Gi = 0,15,30,...165° The value of S is reduced from 0.006 to 0.003 mm-2 by substituting 59° for 45° in a3. This implies that when the deviation of linear systems within a dataset is relatively large, a similar correction method can be employed in manual approaches. The close correspondence between the experimental and modelled roses of intersection densities for the strongest family of striae points towards a large degree of orientation uniformity (comparing Figs. 6a and 7a). The relative smoothness in the experimental curves reveals the image to be quite representative for this purpose. In contrast, the “noised” character of the experimental rose of intersection densities for the weakest system (Fig. 6c) demonstrates the original image not to be representative for characterising this set of striae due to the significant non-uniformity of orientation direction within this set. This is expected in view of the relative small number and especially small lengths of the striae in this family. This reflects the natural law of statistic concept of heterogeneity. The deviations between the modelled rose of intersection densities and the experimental one are obviously flooded in the scatter resulting from non-representativeness. The reflected representativeness of the intermediate system is coupled with restricted orientation uniformity (visible in the original image and the copy). Of course, a larger sample would have provided for more accurate information on non-uniformity of orientation. The results derived in this study using FFT analyses, and the experimental and particularly the modelling approaches using stereological concepts of directed secants, compare favourably with field attempts to determine striae directions using (i) compass directions relative to the local sun angle and (ii) graphically using compass directions towards topographic features (Hättestrand and Stroeven, 1996). Field-derived values thus derived for (a1 - a2) and («2 - a3) are 18° and 72°, and 24° and 83°, respectively. This contribution compares efficiency and reliability of three mathematical and statistical methods for the purpose of characterising the orientation distribution of glacial striae on a bedrock sample. The 2D Fast Fourier Transform allows for the fastest estimation of the preferred directions of orientation, provided that the signal intensities are strong enough. The stereological approach, by applying Saltikov’s directed secants analysis, provides an equally reliable but less efficient assessment of the two dominant groups of striae patterns. Also here, one set of striae is insufficiently well determined, as illustrated using a third methodology of manual weighed-average orientation direction. Based on the rose of intersection densities, the Underwood concept for a combination of three linear systems in a plane allows for the construction of a total rose of intersection densities from the individual roses. In this study, the theoretical model yields satisfactory agreement with the experimental data and field measurements. Deviation of this idealized concept (modelled rose of intersection densities) from the experimental data gives insight into the representativeness of the information derived from the sample, and on the degree of orientation uniformity of the respective families of striae. The original image is representative for orientation distribution analysis with respect to the most dominant family of striae, which shows a high degree of orientation uniformity. To satisfy the fundamentals of sampling science, a larger bedrock sample is required if non-uniform orientation directions are to be accurately determined, such as that of the least dominant striae direction. In this case, a manual approach of weighted-average orientation direction described in this study can be used to improve the reliability of the modelling concept. ACKNOWLEDGEMENTS The authors are grateful to Ms. Liliane Chermant for her contributions to the measurements of preferred orientations of distribution by FFT techniques. 142 Image Anal Stereol 2005;24:135-143 REFERENCES Boulton GS, Clark CD (1990). A highly mobile Laurentide ice-sheet revealed by sattelite images of glacial lineations. Nature 346:813-7. Chaix JM, Grillon F (1996). On the rose of directions measurements on the discrete grid of an automatic image analyser. J Microsc 184:208-13. Coster M, Chermant JL (1989). Précis d’Analyse d’Images, 2nd édition. Paris: Les Editions du CNRS. Denton GH, Hughes TJ (eds.) (1981). The Last Great Ice Sheets. New York: Wiley-Interscience. Ehlers J, Gibbard PL (2004). Quaternary glaciations: extent and chronology, Part I: Europe. Amsterdam: Elsevier. Glasser NF, Bennett MR (2004). Glacial erosional landforms: origins and significance for palaeoglaciology. Prog Phys Geogr 28:43-75. Gonzales RC, Woods RE (1992). Digital Image Processing. New York: Addison-Wesley Publishing Company. Gonzalez LM, Cambrera FL, Sanchez-Bajo F, Pajares A (1994). Measurement of fibre orientation in short-fibre composites. Acta Metall Mater 42(3): 689-94. Hättestrand C (1998). The glacial geomorphology of central and northern Sweden. Sveriges Geologiska Undersökning 85:1-47. Hättestrand C, Stroeven AP (1996). Field evidence for wet-based ice sheet erosion from the south-central Queen Elizabeth Islands, Northwest Territories, Canada. Arc Alp Res 28:466-74. Kleman J, Borgström I (1996). Reconstruction of palaeo-ice sheets: The use of geomorphological data. Earth Surf Proc Land 21:893-909. Kleman J, Hättestrand C, Borgström I, Stroeven AP (1997). Fennoscandian paleoglaciology reconstructed using a glacial geological inversion model. J Glaciol 43:283-99. Kleman J, Hättestrand C, Stroeven AP, Jansson KN, De Angelis H, Borgström I (2005). Reconstruction of paleo-ice sheets - inversion of their glacial geomorphological record. In: Knight PG, ed. Glaciology and Earth's Changing Environment. London, Blackwell Publishing Ltd. Paterson WSB (1994). The Physics of Glaciers. Oxford: Pergamon Press. Redon C, Chermant L, Chermant JL, Coster M (1998a). Assessment of fibre orientation in reinforced concrete using Fourier image transform. J Microsc 191(3): 258-65. Redon C, Chermant L, Coster M, Chermant JL (1998b). Measurements of orientation in a discretised space by Fourier transform: Automatic investigation of fibre orientation in a reinforced concrete. Acta Stereol 17: 93-8. Saltikov SA (1958). Stereometric Metallography (in Russian), 2nd edition. Moscow: Metallurgizdat. Saltikov SA (1974). Stereometrische Metallographie. Leipzig: VEB Deutsche Verlag für Grundstoffesindustrie. Serra J (1982). Image Analysis and Mathematical Morphology. New York: Academic Press. Stroeven AP (2004). Glacier. In: Goudie AS, ed. Encyclopedia of Geomorphology. Routledge, 454-9. Stroeven P, Stroeven AP, Dalhuisen DH (2001). Image analysis of 'natural' concrete samples by automatic and manual procedures. Cem Concr Compos 23:227-36. Sugden DE, John BS (1976). Glaciers and Landscape: A Geomorphological Approach. London: Edward Arnold, 376. Underwood EE (1970). Quantitative Stereology. Reading (Massachusetts): Addison-Wesley Publishing Co. 143