Parameter-sparse modification of Fourier methods to analyse the shape of closed contours with application to otolith outlines

Additional keywords: cod otoliths, Greenland halibut otoliths, modified Fourier method, shape analysis.

The popularity of shape analysis in terms of analysis of two-dimensional (2-D) contours has increased in marine research with application, for example, in stock discrimination, in which the shape of otolith outlines from individual fish in different regions is compared, as reported for Micromesistius poutassou (Keating et al. 2014), Gadus morhua (Jónsdóttir et al. 2006; Stransky et al. 2008) and Scomberomorus cavalla (De Vries et al. 2002). The broader issue of morphometric outlines has been reviewed by Christoph Stransky (2014). A general challenge of shape analysis methods of 2-D contours is to describe important shape features with as few descriptors as possible. Why this is so can be illustrated by the frequently applied Fisher’s linear discrimination method (see Johnson and Wichern 2007), where a one-dimensional (1-D) discrimination function is constructed as a linear combination of all the descriptors involved. Because these are generally to be estimated, it is known that the uncertainty increases with the number of descriptors, making the results unduly uncertain if too many descriptors are chosen.

To approximate otolith outlines with a limited number of descriptors, Fourier methods have proved to be appropriate. A natural explanation for this is that if we imagine tracing the contour several times at the same speed, we can obtain a periodic signal with a perfectly smooth join between each loop. Fourier methods have been developed for the analysis of periodic phenomena.

However, for very localised features (landmarks), wavelet analysis can be a more appropriate approach than Fourier methods (see Sadighzadeh et al. 2014). The recently released R-package shapeR (see http://cran.r-project.org), used for analyses of otolith contours, contains options for both Fourier and wavelet analysis (Libungan and Pálsson 2015).

An extensively used method to transform otolith outlines is the Fourier method using elliptical Fourier descriptors (EFDs) introduced by Kuhl and Giardina (1982). Briefly, the EFD method consists of superimposing ellipses with increasing cyclic frequencies and allows for different distances between succeeding contour points. The continuous outline is approximated with straight lines between succeeding sampled contour points. Relatively few frequency components are needed to obtain a good approximation of complex 2-D outlines. However, the method requires four descriptors (parameters) for each ellipse and a major drawback of the method is that these descriptors are not independent of each other, as shown by Haines and Crampton (2000). Another Fourier method (Haines and Crampton 2000) uses the tangent angle to the contour as a basic variable instead of Cartesian coordinates as used by EFDs. The continuous contour is created by a smooth interpolation between succeeding sampled contour points and a set of contour points is created with equidistant spacing along the contour. Only two Fourier coefficients are needed for each frequency component and these are independent of each other. A third method (the ‘MIRR’ method) was recently reported by Reig-Bolaño et al. (2010) and transforms the 2-D contour to a 1-D signal by mirroring the lower part towards the right. Then, equidistant x-values are chosen, and the corresponding contour points, y, are found by interpolation between the closest contour points. Finally, fast Fourier transform (FFT) is applied to the y-values. This method also benefits from the fact that only two parameters are needed for each frequency component and that these parameters are independent of each other. This method was originally named ‘partial reflection’ by Reig-Bolaño et al. (2010), but was recently denoted ‘MIRR’ (Harbitz and Albert 2015) to emphasise the mirroring aspect.

All the Fourier methods described above provide a poor approximation of a pure ellipse based on only one frequency component, in particular EFD and MIRR. The aim of this paper is to explain why this is so and to provide a simple, general and applicable modification that provides a nearly perfect approximation of a pure ellipse based only on one frequency component. In addition, it is shown that the modified method appears to need fewer parameters than the unmodified methods to approximate the large-scale (low-frequency) shapes of real otolith contours. Briefly, the reason why the original EFD method does not perform well in the case of a pure ellipse is that a constant tracing speed around the contour is chosen. This is not in accordance with a common Fourier description of a pure ellipse, and the modified method applies an appropriate dynamic tracing speed. In the case of the MIRR method, a pure sinusoidal function in one dimension has no vertical parts, in contrast with an ellipse, and the modified method compensates for this by choosing particular horizontal x-values. The methods are used to discriminate between samples of Greenland halibut stocks from southern Greenland and Norwegian waters, as well as between Norwegian coastal cod (NCC) and north-east Arctic cod (NEAC) stocks, where discrimination between the stocks has been documented on the basis of EFDs of otolith contours (Stransky et al. 2008; Harbitz and Albert 2015).

To avoid ambiguities with the MIRR method, the original contours are replaced with lasso contours. The lasso contour in this context was introduced by Harbitz and Albert (2015) because the contour can be virtually seen as the result of the tightening of a rope around the outline, thus providing a non-concave shape. In addition, only one contour point is hit by a radial pointing in any direction from any point inside the contour. By calculating the radial contour points based on radials with equidistant angles from the otolith centroid, it is easy to calculate an average (standardised) contour of several otoliths.

The methods outlined were applied to a sample of Greenland halibut, consisting of 83 otoliths from southern Greenland waters in 2007, and a sample of Greenland halibut from Norwegian waters, consisting of 828 otoliths from the same year. In addition, the methods were applied to samples of NCC (367 otoliths from the northern part of the Norwegian coast) and NEAC (243 otoliths from the south-western part of the Barents Sea area). Stock discrimination results based on classical EFD methods are reported for the comparisons of Greenland halibut (Harbitz and Albert 2015) and cod (Stransky et al. 2008) samples. In both cases, discrimination analyses were also made on the basis of genetics, giving clear indications of differences between the stocks.

In order to avoid false discrimination because of differences in fish length and sex distributions for the Greenland halibut samples, the previous analyses calculated the average of repeated random samples of 83 otoliths from the larger Norwegian water sample in such a way that the same sex–length distribution as for the Greenland sample was obtained each time. For the cod samples, a more approximate approach was chosen by restricting the length of the fish to 30–70 cm.

Let the contour be described by a set of straight lines between N succeeding contour points (x₁,y₁),…,(x_N,y_N) with x_i = x(t_i) and y_i = y(t_i), t_i > t_i-1 where t denotes time as the contour is traced in an anti-clockwise direction. Let s denote accumulated distance travelled along the contour. The following increment notations can be introduced:

for i = 1,…,N with x₀ ≡ x_N and y₀ ≡ y_N. In their construction of EFDs, Kuhl and Giardina (1982) assume a constant speed V = Δs/Δt = 1 along the piecewise linear contour, (i.e. Δs = Δt). Thus, the horizontal and vertical speeds (Δx/Δt and Δy/Δt respectively) are constant between succeeding contour points. An approximate contour (x_k(t),y_k(t)) based on the first k frequency components and 4k + 2 EFDs is now given as follows:

where the centroid coordinates a₀ and c₀ can be calculated as:

and a_n, b_n, c_n and d_n are the EFDs for frequency components n, n = 1,2,…,k:

with and where T = t_N is the perimeter.

Note that the expressions for a₀ and c₀ above are much simpler than the corresponding expressions given by Kuhl and Giardina (1982) and later repeated by Lestrel (1989) and Safaee-Rad et al. (1992). For k = 1, a ‘best’ fitting ellipse is obtained. However, it turns out that even for a pure ellipse the first harmonic in general gives a rather poor approximation (see Fig. 1a). To understand why, we consider the standard way to describe a pure ellipse mathematically.

Fig. 1.  Comparison between a pure ellipse (shaded grey) with a minor/major axis ratio of 0.5 and the approximation based on one frequency component (1FC) and the (a) elliptical Fourier descriptors (EFDs) and (b) ‘MIRR’ methods. The grey and black contours are based on the unmodified and modified methods respectively.

A pure ellipse with semi-major and -minor axes lengths equal to a and d respectively can be parameterised as follows:

where .

Per definition, only one harmonic component is needed to describe this ellipse. Based on Eqn 5 the speed V = ds/dt is:

which is clearly a function of t and not constant, except for a circle when a = d and the square root expression is equal to 1. In the case of the pure ellipse, the concept of a constant tracing speed along the contour appears to be in conflict with the aim of reproducing well a contour with as few harmonics as possible.

The main idea behind improving the EFD method is to choose a tracing speed along the contour that is proportional to V(t) given by Eqn 6 based on the ellipse that gives the best fit to the otolith contour. Although V(t) now varies along the contour, it is reasonable to assume that the speed is nearly constant between two succeeding contour points. Thus, we get:

where V_i is the ellipse speed at contour sample point i and Δs_i is the distance between sample points i and i + 1. Because Δt_i is changed, t_i and T must be changed accordingly:

with i = 1,…,N, and T = t_N. In general, the contour has a best-fitted ellipse with a rotation angle Ψ with a horizontal axis that is generally different from zero. To calculate the appropriate value for V_j, the rotation angle has to be taken into account. This is outlined in Appendix 1 to avoid too much mathematics in the main text.

The redefinitions of Δt_i and t_i by Eqns 7 and 8 respectively, along with the new value T = t_N can now be put directly into Eqn 4 to recalculate the EFD. In this manner, a new best-fitting ellipse is defined by the new values of a₁, b₁, c₁ and d₁, and the procedure can be repeated until convergence. A faster alternative, where no iterations are needed, is to determine a best ellipse outside the EFD framework, for example defined as the ellipse with the same second-moments of area (area moments of inertia) as the region enclosed by the contour. This approach is used, for example, by the image processing toolbox in Matlab (Mathworks Inc., Natick, MA, USA).

In order to compare the shape of individual contours from two different groups (e.g. for stock discrimination purposes), the contour has to be standardised with regard to orientation, size, initial contour point and contour tracing direction. How this is done within the EFD framework is described in more detail in Appendix 1.

The concept of MIRR, denoted ‘partial reflection’ by Reig-Bolaño et al. (2010), is illustrated in Fig. 2. First the otolith contour is rotated so that the best-fitting ellipse is oriented with the major axis horizontally. It is quite obvious that this method applied to a pure ellipse will provide a contour that is vertical in both horizontal ends as well as at the midway point. However, with equidistant x-values, the first Fourier component of the y-values will provide a pure sinusoidal function, which by no means is vertical for any x. So, intuitively, equidistant x-values provide a poor approach to an ellipse, with only one frequency component (Fig. 1). However, we can choose the x-values in such a (non-equidistant) way that the first harmonic reproduces the ellipse virtually perfectly. The way to choose these x-values is as follows:

where x_min and x_max are the minimum and maximum x-values of the original (rotated) contour respectively and N is the even number of x-values. The corresponding vertical contour coordinates y₁,…,y_N will not match perfectly with the original contour points, and must be found by interpolation.

Fig. 2.  A Greenland halibut otolith aligned along the major axis of a best-fitted ellipse. The thick curve is the one-dimensional representation of the non-concave lasso contour, where the lower part is mirrored at a vertical line at the maximum x-value of the contour. Note that using the lasso contour, ambiguities such as those shown at x₀ are avoided.

The result in Eqn 9 is based on the following consideration. For equidistant sampling points in the x-direction, covering exactly one period of a sinusoidal function y = d·sin(x), the FFT of the sampled y-values will only have a contribution at the lowest harmonic; all higher harmonics will disappear. From the right expression for the ellipse in Eqn 5, we see how y is expressed as a function of x in this case of an ellipse. By equating the sinusoidal and ellipse expressions for y, we see how the x-values are to be chosen in order to provide a FFT of the y-values with contribution from the first harmonic only, in the case of a pure ellipse. Note that the x-values in Eqn 9 are independent of the ratio b/a between the minor and major half-axes of the ellipse. In order to standardise the coordinates for discrimination and classification analyses, the x-values are first translated to have a mean value equal to zero. Owing to the symmetry seen from Eqn 9, this is simply obtained by subtracting the mean of these x-values from each of the x-values. The y-values are translated to have zero mean at the midway point between the minimum and maximum y-values, but this translation has no influence on the MIRR descriptors involved in a discrimination analysis. After translation, the x- and y-values are standardised with regard to size by dividing with the square root of the original otolith contour (pixel) area. Note that the number (N) of x-values is arbitrary, so any spatial resolution can be obtained. An example is shown in Fig. 3.

Fig. 3.  Mirrored contour of a Greenland halibut otolith with (a) equal spacing (unmodified ‘MIRR’ method) and (b) the spacing applied with the modified MIRR method. In practice, far more x-values are chosen; a modest number is used here for illustration purposes.

For non-convex contours, ambiguities may occur. One way to get around this is to construct a lasso contour, as shown in Fig. 2. The main algorithmic concept to construct the lasso contour is to choose point by point in such a manner that the line between two succeeding points creates an angle with the horizontal that does not decrease when moving in an anti-clockwise direction. For each point determined, the next point is the first one that gives the smallest change in angle.

An approximation of the contour based on k frequency components is obtained by taking the inverse of the FFT applied to the y-values, where all FFT components for frequencies larger than k are removed. For the ellipse, only the first frequency component is maintained.

To quantify the goodness of fit between the original contour (possibly smoothed) and a fitted contour based on k frequency components, we assume straight lines between succeeding contour points. For the approximate contour we can freely choose the number of contour points (N_k) and a conservative choice would be N_k = 2m with a value of m so that N_k > N. Then, the minimum distance from each original contour point to the (continuous) approximate contour is calculated along with the area (A) enclosed by the original contour. Finally, the average (D_fit) of all these N distances divided by the square root of the area (A) of the original contour is taken as the goodness of fit measure:

where (x_j,y_j) is contour point j on the original contour and (x_kj,y_kj) is a point on the approximate contour based on k frequency components. For each (x_j,y_j) point, the corresponding point (x_kj,y_kj) that is closest to (x_j,y_j) is found. Note that the minimisation is relatively quick because it is easily done using simple algebra, which enables an explicit solution where no numerical loop is needed. If time is an issue, we can easily explore the effect of decimating the number of contour points as well as the original contour.

As a discrimination tool, Fisher’s linear discrimination method combined with cross validation (leave one otolith out at a time) was applied (Johnson and Wichern 2007), as reported previously for discrimination analyses of the Greenland halibut samples (Harbitz and Albert 2015) and cod samples (Stransky et al. 2008). In addition to the use of the modified descriptors, the present analyses were run based on the lasso contour, whereas the previously published cod results were based on raw contours.

EFDs for the pure ellipse

Interestingly, the results given in Table 1 of a = 54.852 and b = 32.191 (original EFD) are very close to the results a = 54.915 and b = 31.963 found by Safaee-Rad et al. (1992) on a similar elliptic test case with true values a = 60 and b = 30. However, the author does not agree with their statement that the substantial error is ‘…expected due to quantisation errors’. The modified method provides virtually perfect results, even with a non-smoothed pixel contour, clearly confirming that it is the concept of a constant tracing speed that is the main error explanation in this case. And in this pure mathematical example, quantisation errors are negligible.

Table 1.  Test Case 1: elliptical Fourier descriptors (EFDs) applied to a pure mathematical ellipse
Input: a = 60, b = 30, ψ = 0, x_i = a cos θ, y_i = b sin θ, θ_i = Δθ,2Δθ,…,2π, Δθ = 2π/4096

The results above were obtained within the EFD framework (i.e. new ellipse parameters were calculated for each iteration). Using the pure ellipse parameters as input to the EFD equations to find the appropriate velocities to apply, the results a = 59.999982 and b = 29.999997 were obtained (i.e. virtually perfect results without iterations).

For the original EFD in Table 2 we get very similar erroneous results for a and b as in Test Case 1. Again, superior results are obtained with the modified method for a and b, as illustrated in Fig. 4c, despite the fact that now obvious quantisation effects are present. Again, a similarly good result was obtained by applying the second-moment approach to find the best-fitting ellipse and thus avoiding the iteration process.

Table 2.  Test Case 2: elliptical Fourier descriptors (EFDs) applied to an ellipse contour extracted from a black-and-white image
The black and white pixel image of an ellipse with a = 60, b = 30 and ψ = 45° is shown in Fig. 4a, along with the extracted contour (x₁,y₁),…,(x_N,y_N) with n = 256, based on intensity thresholding

Fig. 4.  (a) A black-and-white image of an ellipse with a minor/major axis ratio of 0.5 rotated 45° along with the detected contour. (b, c) Reproduction of the ellipse based on elliptical Fourier descriptors (EFD) and one frequency component. The dashed line indicates the classic (unmodified) method, whereas the dotted line (virtually a perfect fit) shows results with the modified EFD method.

For the MIRR method, there is no inherent method to find a best rotation angle of the original contour. The ellipse with the same second-moments of area as the original contour is defined as the best-fitting ellipse in this study, and the contour is rotated with an angle equal to the angle between the major axis of the best-fitting ellipse and the horizontal x-axis. With the same data as in Test Case 2 (pixel contour from a black-and-white image), the rotation angle was found to be exactly 45°. Based on one frequency component and equidistant x-values, the MIRR method provided a minor/major axis ratio equal to 0.5696 (i.e. a bit closer to the true ratio 0.5 compared with the ratio of 0.5869 obtained by the (unmodified) EFD method). Applying the modified MIRR with 2048 non-equidistant x-values yielded a ratio 0.5005 and the approximate contour was very close to a perfect ellipse, as in the modified EFD case (see Fig. 1).

The goodness-of-fit deviation measure (D_fit) between approximate and true lasso contours was calculated for all the 83 Greenland halibut contours as well as the 367 NCC contours for the number of frequency components (k) equal to 1, 2,…,10, 15 and 20. Note that in these runs, the best-fitted ellipse is based on the second-moment criterion and not the inherent criterion defined within the EFD framework. The results are shown in Fig. 5. As is clearly seen from Fig. 5a–d, the modified method provides considerably better fit for the smallest frequencies. In particular, an apparent difference is seen for the MIRR method (Fig. 5a, b). For the EFD method, the difference seems to disappear for k > 2 in the Greenland halibut case.

Fig. 5.  Goodness-of-fit results for (a, c, e) Norwegian coastal cod (NCC) and (b, d, f) Greenland halibut otoliths. The greatest improvement is seen for the modified ‘MIRR’ method applied to the Greenland halibut otoliths (b). k, number of frequency components; EFD, elliptical Fourier descriptor. See text for further details.

The MIRR and EFD methods are compared in Fig. 5e, f. As is seen, the two methods perform rather equally. Note, however, that the number of coefficients in the case of the EFD method is 4k, whereas it is 2k in the case of the MIRR method. In addition, the MIRR coefficients are independent of each other, whereas the EFDs are not (Haines and Crampton 2000).

An example of an approximate lasso contour fit to a Greenland halibut lasso contour is shown in Fig. 6. The ‘ringing’ phenomenon clearly seen for the unmodified MIRR is a typical behaviour in Fourier approximation to time series in one dimension, when abrupt changes in the y-values are present. This noise phenomenon seems to disappear with the modified method.

Fig. 6.  Comparison of original lasso contour of Greenland halibut otolith and approximation by (a) unmodified and (b) modified ‘MIRR’ methods. Note that the ‘ringing’ phenomena and sharp corners at the minimum and maximum x-values for the unmodified MIRR have completely disappeared for the modified MIRR. k, number of frequency components.

No apparent different discrimination results were found for the Greenland halibut and cod samples compared with the reported results based on the classical methods either in terms of estimated probabilities of correct classification of a single otolith or estimated standard deviations of these probabilities. This issue is covered in more detail in the Discussion.

Calculations were performed using Matlab R2015a for Macintosh (Mathworks Inc.) with a Macintosh OS X Yosemite Version 10.10.1 processor (Apple Inc., Cupertino, CA, USA) and 8-GB RAM. The time taken to obtain the goodness of fit measure D_fit was typically ~1 s per otolith with a conservative choice of the number of contour points (e.g. 4096). However, a negligibly different result was obtained with a decimation of a factor of 10 with regard to the original as well as the approximate contour. Thus, in this case, between 0.01 and 0.02 s was sufficient per otolith.

To compare the efficacy of the modified versus the original method in the discrimination analyses, a contour must be available for both EFD and MIRR methods and it must be standardised in terms of rotation. Commands in Matlab were used here to extract the contour. A typical cod example was a 572 × 720-pixel image with a 1636-pixel contour. To extract the contour took ~0.02 s. Another 0.03 s was needed by MIRR but not by the classical EFD method to provide the elliptical parameters for the best-fitted ellipse based on second area moments (i.e. the major and minor axes lengths, the rotation angle and the centroid). To convert from a raw pixel contour to the lasso contour, another 0.06 s was typically needed for a 3000-pixel cod otolith contour and 0.04 s for a 7000-pixel Greenland halibut contour. To provide EFDs with 50 frequency components, ~0.03 s was needed. Thus ~0.3 s was needed if the modified method required 10 iterations. If the best-fit ellipse parameters are available, no iterations are needed and the difference between the classical and modified methods was negligible with regard to EFD computing time. For the calculation of the MIRR descriptors, the whole NCC sample was tested and ~0.02 s was required based on 2048 x-values, with negligible difference between the classical and modified approaches.

The discrimination program further required 0.07 s to compare two groups where each group consisted of 83 individuals and there were 40 descriptors for each individual. To conclude, only a small fraction of a second was needed per otolith from contour extraction to the final discrimination. Thus, extensive analyses with replicates based on bootstrapping of individual otoliths are fully feasible, even over a couple of hours.

For a pure ellipse, the modified EFD and MIRR methods reproduce the ellipse virtually perfectly, whereas the original methods perform rather poorly, with only one frequency component. The examples in the present study are a pure mathematical ellipse, with the horizontal axis along the major ellipse axis, and a black-and-white image where the ellipse major axis makes a 45° angle with the horizontal axis. In both cases, the true ratio between the minor and major axes was 0.5. In the case of the image, a modest contour length of 256 was chosen to enhance a possible pixel noise effect. In an actual study, these exercises can be easily reproduced with more appropriate choices of axes lengths, rotation angles and perimeter lengths.

The ellipse case studies demonstrated a rather fast convergence rate with the modified EFD method and indicated that an upper threshold based on the absolute value of the relative change between two succeeding iterations could be used as an automatic termination criterion. However, it remains to be demonstrated that convergence will be obtained in a general case. If convergence problems occur, or the iterations take too long, it is straightforward enough to apply an approach outside the EFD framework to fit the best ellipse and then apply the EFD to calculate the higher-order predictors.

When applied to the actual cod and Greenland halibut otolith lasso contours, the most apparent improvement of the approximate lasso contour with the modified methods was seen for the first few frequency components. However, even for 20 frequency components the modified methods performed slightly better than the unmodified methods, and for a larger number of frequency components there is a negligible difference between the classical and modified methods. So, for the purpose of reproducing approximate contours, there seems to be little to lose by using the modified methods, also when taking the time required into account.

An explanation as to why the discrimination results for the modified methods applied to the Greenland halibut and cod stock samples did not show any apparent difference from the unmodified EFDs could be that the magnitude and variance of the descriptors decrease rather rapidly with frequency, and there are rather strong negative correlations involved. As a consequence, the general rule of an increased uncertainty with an increased number of estimated descriptors may not be true in this case. However, the lack of improved discrimination is not necessarily an indication that the modified methods do not represent any improvement in such analyses. Having several methods that give the same results enhances the reliability, and different results indicate that the analyses should be interpreted more carefully. In fact, good discrimination can be obtained by a method that gives an apparently poor approximation of the contour because of an indirect and not easily interpretable effect where stock differences remain despite a poor contour reproduction. Good and false discrimination results can also be easily obtained because of the pitfalls described by Harbitz and Albert (2015). An interesting challenge is also to investigate whether a combination of predictors from different methods may improve discrimination.

The analyses of actual otolith contours demonstrated herein are limited to lasso contours with non-concave shapes, one reason being to avoid ambiguities with the MIRR method. Although there is a rather strong (pixel noise) smoothing effect in the lasso contour, the abrupt direction changes may introduce unwanted effects, so to study the effect of smoothing is an interesting subject even for lasso contours. The lasso method also effectively masks effects like ripples, which, in principle, may be important features for stock discrimination purposes, for example. So, investigations of how the modified method performs on the original contour are also interesting subjects for the future.

It should be noted that documented discrimination results based on more than 10 frequency components are rare and too many predictors will easily provide too much noise and numerical challenges. Interestingly, the same large score probabilities were obtained in the case of the cod study based on the modified EFDs and the lasso contour as for the previous studies based on classical EFDs and the raw contour. This indicates that the shape difference is on a rather large scale, as demonstrated by the approximate contour in Fig. 7, as well as the comparison between the average shapes in Fig. 8. In case of important small-scale features, the use of EFD is not necessarily an appropriate approach. The lasso contour could, in principle, also be used to study small-scale effects, but by indirect techniques, such as using the difference between the length of the raw contour perimeter and the lasso contour perimeter as a supplemental predictor.

Fig. 7.  An example of a Norwegian coastal cod otolith with ripples, and a reproduction of the contour (black line) by inverse elliptical Fourier descriptors (EFD) based on 10 frequency components.

Fig. 8.  Comparison of the average contour shape of whole Norwegian coastal cod and north-east Arctic cod otoliths.

There is a popular third Fourier method mentioned in the Introduction that is based on tangent angles to the contour at equidistant points along a smooth interpolation between successive sampled contour points (Haines and Crampton 2000). This method also does not perform well on a pure ellipse, based on one frequency component, although it performs far better than the methods analysed herein. The tangent angle method can also rather easily be modified for a pure ellipse, but it is much harder to find a general method that is also applicable to a general outline, so here is a new challenge for the future.

The data material and program codes needed to reproduce the results and produce new ones are available from the author (except for some codes in Matlab’s image processing toolbox). All codes are written in Matlab, and the lasso contour script is also written in R.