Improved Anisotropic Gaussian Filters

Elongated anisotropic Gaussian filters are used for the orientation estimation of fibers. In cases where computed tomography images are noisy, roughly resolved, and of low contrast, they are the method of choice even if being efficient only in virtual 2D slices. However, minor inaccuracies in the anisotropic Gaussian filters can carry over to the orientation estimation. Therefore, this paper proposes a modified algorithm for 2D anisotropic Gaussian filters and shows that this improves their precision. Applied to synthetic images of fiber bundles, it is more accurate and robust to noise. Finally, the effectiveness of the approach is shown by applying it to real-world images of sheet molding compounds.


INTRODUCTION
Gaussian filters have a wide variety of applications in image processing.Whereas isotropic Gaussian filters, being the foundation of scale space theory (Lindeberg, 1996), can be implemented easily, their anisotropic counterparts are more demanding while being just as interesting (Lampert and Wirjadi, 2006): They give a handle on orientation as well as scale, which makes them the cornerstones of orientation space theory (Faas and van Vliet, 2003).Anisotropic Gaussian filters have been employed for denoising images (Yang et al., 1996;Treece, 2020) as they can be adapted to image structures and, hence, preserve edges.Another application is the estimation of orientations using a filter bank of anisotropic Gaussian filters.For example, local fiber directions can be estimated by the direction of the maximal response of anisotropic Gaussian filters that are aligned along a given set of directions (Robb et al., 2007;Wirjadi et al., 2009).
Most established methods for estimating fiber directions are based on gradients, such as calculating the image's Hessian matrix (Eberly et al., 1994;Frangi et al., 1998;Ohser and Schladitz, 2009) or the second-moment matrix of the image's gradient, called the structure tensor (Haglund, 1992;Weickert, 1999;Krause et al., 2010).In both methods, the local fiber directions result from the eigendecomposition of the respective matrix.Wirjadi et al. (2016) compared them to other methods based on 3D images of single synthetic fibers and identified them as the most accurate.Pinter et al. (2018) further investigated their accuracy on images of multiple fibers, finding the structure tensor to be comparably more robust.
In their comparison, Wirjadi et al. (2016) and Pinter et al. (2018) both identified the Maximal Response (MR) method as robust to noise but suffering from the trade-off between runtime and accuracy in 3D.However, there are use cases where the image quality is too low for methods based on local grayvalue derivatives on the one hand and where due to the production process the fibers are known to be oriented in a 2D subspace anyway.This holds for instance true for µCT images of sheet molding compounds, a material where reinforcing fibers lie within a plane.Then, only 2D images have to be analyzed, similarly to stereological approaches.In 2D, the MR method is less restricted regarding runtime and even outperforms other methods due to its robustness with respect to low image contrast (Schladitz et al., 2016).
The accuracy of the direction estimation clearly depends on the accuracy of the filter responses for the considered directions.Under otherwise perfect conditions, this method's results barely depend on contrast as the filter responses scale with the contrast.Although the response differences are smaller, this does not influence which response is maximal.However, computed tomography images are often affected by noise and other artifacts.For low-contrast images, these effects have a much stronger impact on the detected direction of maximal response due to the small differences in responses for varying angles.In this case, small inaccuracies in the anisotropic Gaussian filter can impair the direction estimation further.In this paper, we will consider the case of low contrast between the foreground, i.e., fibers, and noise, while using a low resolution for the fibers.
Anisotropic Gaussian filters in R 2 can be implemented naïvely by filtering in the directions of the major and the minor axis of the Gaussian's contour ellipses, subsequently.However, Geusebroek et al. (2003) derived a more accurate decomposition, where at least one of two filter directions is aligned with an axis of the image grid.Whereas the naïve implementation may need interpolation for filtering in both directions, Geusebroek et al.'s method requires interpolation for at most one filter direction.Lampert and Wirjadi (2006) generalized these results to R d and provided explicit formulas for R 3 .
Besides implementational inaccuracies of the 1D filters, interpolation introduces spatial inhomogeneity into the filter kernels, as Lam and Shi (2007) have shown.Therefore, they propose a modification of Geusebroek et al.'s method which avoids interpolation altogether at the cost of an additional 1D Gaussian filtering step.However, Lam and Shi's modification limits possible half-axis ratios ω = σ 2 /σ 1 , to ω ≥ 0.4142.The ratio can be lowered to ω ≥ 0.1622 at the cost of aliasing effects.In our setting, far smaller ratios are needed to accurately mimic the fiber shape, e.g., ω = 0.025 in the Application section.
In this paper, we therefore suggest another modification not suffering from this restriction.We propose a modification to Geusebroek's decomposition that halves the number of interpolation steps and show that this modification improves performance.Moreover, we consider cubic instead of linear interpolation, which improves accuracy at the cost of speed.Based on synthetic fiber images, we show that the adapted method results in higher accuracy of the maximal response method.Additionally, we show that it outperforms estimation based on the Hessian matrix or the structure tensor in 2D.Finally, we apply our method to real-world images of sheet molding compounds.

MATERIALS AND METHODS
A natural approach to calculating the anisotropic Gaussian filter in R d is to decompose it into a sequence of multiple Gaussian filters in R, which poses a simpler problem (Lampert and Wirjadi, 2006).The recursive scheme with infinite impulse response by Young et al. (1995;2002), using boundary conditions by Triggs and Sdika (2006) has proven efficient and accurate.For the case of R 2 , Geusebroek et al. (2003) propose a decomposition into filters along the x 1 -axis of the image grid and a filter along another direction that generally does not align with the grid.
Initially, consider an axis-aligned Gaussian kernel with standard deviations σ 1 > σ 2 > 0 centered in the origin, i.e., (2) Its contour lines are axis-aligned ellipses with halfaxis ratio ω = σ 2 /σ 1 .We now rotate the kernel to get g σ 1 ,σ 2 ,θ , whose major half axis points in direction ν = (cos(θ ), sin(θ )) T for θ ∈ [0, π).Formally, where A decomposition of the corresponding filter into one-dimensional filters along the coordinate axes is generally not possible.However, Geusebroek et al. (2003) proved that a decomposition into filters along the x 1 -direction and the direction is indeed possible, namely with the kernel The standard deviations σ x , σ ν * can be computed in terms of the rotation angle θ and the standard deviations σ 1 , σ 2 via Fig. 2 illustrates this decomposition; see Geusebroek et al. (2003) for the detailed derivation.

ALGORITHMS FOR ANISOTROPIC GAUSSIAN FILTERING
Gaussian filters can be implemented naïvely by rotating the image with the same matrix that rotates the filter such that its major and minor axes are aligned with the coordinate axes.Then, the image can be filtered along the coordinate axes using Young et al.'s (2002) recursive Gaussian filter.This corresponds to the filter decomposition along the principal axes, see Eq. 3.This way, more memory is consumed as the image does not fit its previous rectangular structure anymore.Moreover, interpolation steps are necessary for both filter directions (Lampert and Wirjadi, 2006).
More advanced algorithms make use of decomposition along other axes, as in Eq. 7: Lampert and Wirjadi's(2006) geometric algorithm circumvents the interpolation along one axis by considering the filter decomposition as a shear of the coordinate axes with a shear matrix V .Hence, the image is sheared with V before filtering along the coordinate axes.Afterwards, the resulting image is transformed back with V −1 .
In Geusebroek et al.'s (2003) line buffer algorithm (Lampert and Wirjadi, 2006) the image is processed in-place as it filters along the x 1 -axis and the ν * -line, see above.The transformation step necessary for filtering along the ν * -line is the inverse shear used in the geometric algorithm.This transformation using interpolation is necessary every time data is read or written.This can be kept minimal by using image line buffers for the filter history.However, as the recursive Gaussian filter consists of a forward and a backward filter, this yields 2 forward and 2 backward transformation steps, yielding 4 interpolation steps per pixel.
In comparison, the geometric algorithm uses only 2 transformations and, thus, interpolations per pixel, which makes it less error-prone compared to the line buffer algorithm.However, the geometric algorithm needs more memory because the transformed image no longer fits the original rectangular shape.

THE HYBRID ALGORITHM
Our improved scheme combines the advantages of both the geometric and the line buffer algorithm: It filters in x 1 -direction with Young's (2002) recursive Gaussian filter as in the line buffer algorithm.The filter in ν * -direction is modified such that the intermediate transformation steps are omitted: As the forward and backward filter move along the same line, the intermediate transformation steps taken together are the identity.Therefore, the result of the forward filter does not need to be transformed but can be stored in-place.This approach requires 2 interpolation steps per pixel, as in the geometric algorithm, while using as little memory as the line buffer algorithm.The difference to the established algorithms is "smarter bookkeeping".Hereafter, we will call this the hybrid algorithm 1 .
An axis-aligned filter is generally more accurate than a filter that is not axis-aligned since the latter requires interpolation.Therefore, we first filter along the axis and, subsequently, in ν * -direction.
So far, we only discussed a decomposition into filters where one filter direction is aligned with the x 1 -axis.Analogously, a decomposition such that one filter direction is aligned with the x 2 -axis is possible (Lampert and Wirjadi, 2006).This may even be advantageous for 45 • ≤ θ ≤ 135 • : The standard deviation σ x of the filter along the x 1 -axis varies over the rotation angles θ , being largest for θ = 0 • and smallest for θ = 90 • .For the line buffer and the hybrid algorithm, filtering along the x 1 -axis smoothes the image in the same direction, in which the interpolation takes place.This may be less error-prone for stronger smoothing.Hence, we propose to decompose the anisotropic filter with an x 2 -aligned axis for 45 This modification is possible for each of the approaches mentioned above.In the following, we call this the major-axis modification.

THEORETICAL PERFORMANCE ANALYSIS
The runtime of the anisotropic Gaussian filter is constant for each pixel and depends only on the rotation angle and not on the variance.The filtering steps require 12 additions and 13 multiplications.Linear interpolation can be implemented with 1 addition and 2 multiplications per pixel.
We further propose to apply cubic interpolation with natural boundary conditions.In our implementation, each cubic interpolation step takes 8 additions and 14 multiplications per pixel.Therefore, we only combine it with the hybrid algorithm.The total complexities per pixel are listed in Table 1.The runtime of the MR method in total depends on the complexity of the employed anisotropic filter algorithm, the discretization of the direction space, i.e., the number of angles considered, and the image size.The dependency on the latter three is linear, thus we discuss the speed of Gaussian filters, only.

FIBER ORIENTATION ESTIMATION Maximal Response of Anisotropic Gaussian Filters
To filter an image of fibers for directions, imitate the elongated shape of a fiber with the d-dimensional anisotropic Gaussian (function) g θ , see Eq. 1.Its parameters give a handle on the orientation, length, and diameter for the fiber model (Lampert and Wirjadi, 2006).
The filter response (g g g θ * f f f )(x x x) to the image f f f is maximal when θ matches the local fiber direction in the point x x x ∈ Ξ, where Ξ is the fiber system.Therefore, one can find the direction ν that maximizes the filter response for all x x x ∈ Ξ (Wirjadi et al., 2016): ) is estimated by calculating the convolution(g g g θ * f f f )(x x x) for a finite set of directions that covers the space as evenly as possible (Schladitz et al., 2016).Hereafter, we will call this the MR method.The method's accuracy is mainly influenced by the parameter σ 2 , which we will set to r/2 in the following.This is motivated by the 2σ rule for the normal distribution, which says that approximately 95% of the data points are within two standard deviations of the mean (Georgii, 2012).Hence, a correctly aligned filter kernel g θ that is centered within the fiber covers the fiber's thickness with 95% of its weight when σ 2 = r/2, where r is the fiber's radius.Pixels that are further away than 2σ 2 are barely taken into account.This ensures that the filter response is maximal when the 2-dimensional Gaussian filter kernel is aligned with the fiber: a much larger variance might take too many pixels outside of the fiber into account, while a much smaller variance results in filter kernels whose main mass is concentrated in an elliptical region that is thinner than the fibers.In this case, the ellipse might fit inside the fiber for several angles θ which makes it harder to accurately detect the direction for which the maximum is attained.The parameter σ 1 corresponds to the elongation of the fiber.We observed that it has rather little influence on the estimation accuracy.Yet, it should hold σ 1 < ℓ/4, where ℓ is the fiber length, analogous to the argument above, and the Gaussian kernel should still be clearly elongated.

Structure Tensor
The image's gradient describes directions in the image as it is aligned with the fiber's surface normal, which in turn is orthogonal to the fiber direction.Further convolution with a Gaussian kernel allows for robust direction estimation against noise.
For smoothing parameters σ , ρ > 0, let the isotropic Gaussian kernel g σ := g σ ,σ and ∇ f the gradient of the image f .The structure tensor is defined as as denoted in Wirjadi et al. (2016) and proposed by Haglund (1992).Its eigenvector corresponding to the smallest eigenvalue describes the local fiber direction.

Hessian Matrix
The image's Hessian matrix with smoothing parameter σ > 0 describes image curvature, which is minimal along the fiber direction.Therefore, the local fiber direction is determined by the Hessian's eigenvector corresponding to the smallest eigenvalue.Analogous to the structure tensor, we set σ = r following Wirjadi et al. (2016) and Pinter et al. .

RESULTS
In this section, we support the theoretical analysis with experimental results.More precisely, we show that the hybrid algorithm is more accurate than the line buffer algorithm.Employing linear interpolation, the hybrid algorithm is indeed faster than the line buffer algorithm.In the first subsection, we test the performance of anisotropic Gaussian filters as such.In the second subsection, we test performance on synthetic fiber bundles with varying noise contrast to the background.In the third subsection, we apply the algorithms to real-life data.
The following experiments were carried out on an Intel(R) Core(TM) i7-7500U CPU @2.70 GHz with 16 GiB of RAM, using the GNU compiler GCC 9.0 on a 64-bit GNU/Linux operating system.

PERFORMANCE OF ANISOTROPIC GAUSSIAN FILTERS
We test the performance of the anisotropic Gaussian filters for the line buffer algorithm using linear interpolation and for the hybrid algorithm using linear interpolation as well as cubic interpolation with natural boundary conditions.

Accuracy
We reconstruct the anisotropic Gaussian filter kernel by calculating the unit impulse response, i.e., applying the anisotropic Gaussian filter to an image of size N × N with N = 512, in which all pixel values are 0 except one pixel in the image center with pixel value 1.For each algorithm and variance combination considered here, we compute the l 2 -deviation between the reconstructed kernel ĝ ĝ ĝθ and the actual kernel g g g θ as a measure of accuracy, i.e., The mean and maximum deviation for the rotation angles θ = 0 • , 1 • , ..., 179 • are reported in Table 2.
The hybrid algorithm with linear interpolation yields more accurate results than the line buffer algorithm.Cubic interpolation is even more accurate, except for σ 1 = 7.0, σ 2 = 4.0.This is most likely due to ringing artifacts, i.e., oscillations of the interpolation kernel, which is a known problem of cubic interpolation (Lehmann et al., 1999) also known as the Runge phenomenon (Gautschi, 2011), or, more generally, Gibb's phenomenon (Hou and Andrews, 1978).However, cubic interpolation improves the approximations substantially for variance combinations that otherwise yield comparably large errors for linear interpolation.Note that smaller variances go along with larger errors because the Gaussian approximation is less precise there, see (Young and Van Vliet, 1995).
The major-axis modification achieves even higher precision compared to its counterpart without modification.Notably, the hybrid algorithm with linear interpolation and major-axis modification often outperforms the hybrid non-modified algorithm with cubic interpolation.
For elongated Gaussian kernels, the l 2 -error changes considerably over all rotations: It is lowest for small angular deviations from the x-axis, that is, θ = 0 • .Between 50 • to 130 • it is considerably larger, peaking around 90 • .This conforms with our motivation for the major-axis modification.Employing the hybrid algorithm, the deviations shrink significantly, see Fig. 2a.

Throughput
We test the algorithms' data throughput by applying the filter 50 times to Gaussian noise images of sizes N × N with N = 100, 130, ..., 4 990 and calculate the trimmed mean excluding top and bottom 10%.Fig. 2b shows that the hybrid algorithm with linear interpolation is slightly faster than the line buffer algorithm, at least for larger image sizes.Cubic interpolation, however, takes considerably more time.This conforms to the theoretical results above.The major-axis modification further slows down the algorithms.For 45 • ≤ θ ≤ 135 • , the filter in ν * -direction iterates over all image columns, while the image pixels are saved adjacently within a line.Therefore, memory access is more expensive than it is without modification, the more so, the larger the image.This can be circumvented at the cost of memory by saving the image adjacently within a column for 45 For all three algorithms without modification, there are two different speed plateaus in the throughput, see Fig. 2b.As Lampert and Wirjadi (2006) argue, the throughput is dependent on the image size: In our implementations -for the hybrid as well as the line buffer algorithm -we use 4 buffers for the filter history while reading from and writing into the same image.For small image dimensions, these buffers fit into the CPU's L1 data cache.For larger image sizes, the buffer sizes exceed the cache size slowing down the computations.In our test setup the drop at N = 1 690 corresponds with the system's size of the L1 data cache, namely 64 KiB, see Fig. 2b.

EXPERIMENTAL VALIDATION OF THE MR METHOD
The experiments in the previous subsection have shown that anisotropic Gaussian filters are generally more accurately calculated with the hybrid than with the line buffer algorithm, especially with the majoraxis modification.This section will show that these results translate to the accuracy of the MR method.

Setup
We evaluate the MR method on synthetic images with known constant fiber orientation.The design of the images is inspired by our application example.There, bundles of nearly parallel thin fibers form the main building block of the microstructure.The synthetic images shall mimic the fiber system within one such bundle.
Given an image of size 512 × 512 pixels, a width parameter w, and an angle θ , we define an image F θ ,w by setting For each θ = 0 • , 1 • , 2 • , ..., 179 • and w = 1, 2, we generate such an image.These gray-value images represent idealized fiber bundles with known fiber direction and a radius of r = πw/2 pixels.Note that we choose a significantly smaller radius than Wirjadi et al. (2016) andPinter et al. (2018), who use cylinders with a radius of at least 2.5 voxels.Our choice is again inspired by our application example.
As background noise, we generate images B of size 512 × 512 with pixels sampled from the uniform distribution in [0, 1].
To model images of varying contrast between background and fiber, we consider images of the form (1 − c)B + cF θ ,w for c ∈ [0, 1], see Fig. 3.The MR method is applied as described below, which yields a mean absolute angular error MAE w.r.t. the known fiber direction.For comparison, we additionally Fig. 4: Mean angle error for 50 noise images overlayed with synthetic fiber images with direction θ = 0 • , ..., 179 • , and with varying contrast.For each noise image contrast combination, the MAE's maximum over fiber directions was calculated.The mean and standard deviations over 50 noise images are depicted as point symbol with bars for each contrast and algorithm.Note, however, that the standard deviations are small and therefore the bars delimiting the interval are in many cases covered by the symbol for the mean value.
apply the algorithm to images preprocessed by a median filter of size 3 × 3.This is motivated by the fact that smoothing with median filters is a typical preprocessing step for real image data.
The MAE is determined as follows.The synthetic fiber image F θ ,w is binarized with a threshold of 0.75.The resulting image serves as a mask to include only pixels within fiber cores.
Further, we want to make sure that we only estimate the estimating bias and do not confound it with a sampling bias, i.e., lower or higher sampling probability for certain directions.For example, in an image, one can place more fibers and longer fibers in diagonal directions than in the horizontal and vertical direction.Therefore, we only evaluate pixels within a circle around the image's center.In order to avoid boundary effects, the circle radius is set to 206 pixels.
For each realization B, we are interested in the maximum of the MAE over all θ per method.We run the MR method with σ 1 = 20.0,σ 2 = 0.75w≈ r/2, the structure tensor with σ = r, ρ = 6.0 and the Hessian matrix with σ = r as argued in the Methods section.

Results
Fig. 4 shows the mean error for unfiltered images and the standard deviations for 50 noise realizations.Considering the MR method, the hybrid algorithm outperforms the line buffer algorithm, especially for low-contrast cases.The cubic interpolation is more accurate than the linear interpolation, especially for the high contrast, but also for the low-contrast setup.
The hybrid algorithm with linear interpolation and major-axis modification performs nearly as well as it does with cubic interpolation.For cubic interpolation, however, it performs just as well with the major-axis modification as it does without it.The effect of low contrast is reduced for the larger fiber diameter of r = π in comparison to r = π/2.
The structure tensor, however, is affected by low contrast for both r = π/2, π.For r = π/2, it is about as accurate as the MR method using the line buffer algorithm, so it is outperformed by the hybrid algorithm.For r = π, it is always outperformed by the MR method with a hybrid algorithm and modification.Note that for higher contrast, it still performs better or as well as the MR method using the line buffer algorithm, but for the lowest contrast, it is outperformed by all algorithms employed in the MR method.
The Hessian matrix is not able to match the performance of the other methods, see Fig. 7(b).
Applying a median filter to noisy images is a common preprocessing step to get rid of noise while preserving edges.However, the errors are considerably larger than for unfiltered images as direction information is lost by the undirected median filter, see Fig. 7(a).

APPLICATION TO SHEET MOULDING COMPOUNDS
In this section, we apply the MR method to low-contrast image data of sheet molding compound materials.Sheet molding compounds (SMC) are a type of material consisting of stacked layers of fibers.In the automotive industry, SMC are of high interest due to their versatile behavior such as light weight, high stiffness, and strength, which is determined by their fiber direction distribution (Orgéas and Dumont, 2012).Computed tomography imaging of SMC is challenging due to the high fiber volume fraction and the low difference in X-ray absorption of fiber and matrix material.In the following, we will use the MR method both for fiber direction estimation and fiber enhancement, which is useful for segmenting the fibers.

Sheet Molding Compound with Glass Fibers
First, we consider an SMC material with glass fibers, see Fig. 5.The image was taken using the µCT device at the Fraunhofer ITWM, Kaiserslautern, Germany, with a voltage of 120 kV, an integration time of 999 ms, and 1 200 projections/angular steps.
The device uses a Feinfocus FXE-225 X-ray tube and a PerkinElmer detector with 2 048 × 048 pixels (Fraunhofer ITWM, 2022).As specified by the manufacturer, the fibers' diameter is 10 µm.The material was scanned with a pixel spacing of 5 µm, deliberately undersampling the fibers for the sake of imaging representative sample volumes.
Following the arguments from the Methods section, we applied the line buffer and the hybrid algorithm with both linear and cubic interpolation to the sample with σ 1 = 0.5, σ 2 = 20.0.For the structure tensor and the Hessian matrix, we used σ = 1.0 and ρ = 6.0.Additional to the direction estimation, the MR method supplies a fiber enhancement with its maximal response.Combined with Frangi et al.'s (1998) enhancement filtering, we employ it to segment the fiber system.Subsequently, we add a postprocessing step based on Sliseris et al.'s work (2015).For further details see Niedziela et al. (2024).
For these images, we do not have a ground truth, but visually, the masks produced by the MR method all appear quite accurate, see Fig. 5.Note that for very diffuse fiber bundles, the hybrid algorithm using cubic interpolation looks the most accurate.This result is consistent with its higher accuracy for the response to a unit impulse: The higher accuracy very likely translates to the response to fibers, which makes the enhanced fibers even more distinguishable, the basis of our segmentation approach.

Sheet Molding Compound with Carbon Fibers
As a second application, we consider images of the material SMCarbon® 24 CF50-3K by POLYNT Composites Germany GmbH.It consists of carbon fibers with a length of 25 mm within a vinyl ester resin.The fiber diameter is not known directly as it also changes under pressure.
The sample was scanned with the X-ray microscope Xradia 520 by Carl Zeiss Microscopy GmbH (Leibniz-Institut für Verbundwerkstoffe, 2022) with a pixel spacing of 24.93 µm, a voltage of 60 keV, a power of 5 W, and 3 201 projections.The exposure time was 2 s, where 20 single images were taken with an exposure time of 0.1 s and then averaged.For the field-of-view, they used 76 mm × 48 mm.We applied our versions of the MR method with σ 1 = 25.0,σ 2 = 2.0.Following Schladitz et al.
(2016), we binarized the maximal filter response using Niblack's local thresholding (1986) with a window size of w = 4σ 2 , and the threshold 0.6.Further, we excluded components that have a pixel size lower than 100 after eroding the mask with a square of size 2 × 2. For the sake of consistency, we applied the structure tensor and the Hessian matrix with σ = 4.0 and ρ = 6.0.Despite there barely being any contrast within the fiber bundles, the MR algorithm provides a fairly accurate estimation of fiber directions, whereas the gradient-based methods are struggling, see Fig. 6.We compared the histograms of different algorithms for the MR method, see supplementary material.Strikingly, the line buffer and the hybrid algorithm using linear interpolation without modification apparently overestimate the direction 90 • while underestimating the neighboring 89 • and 91 • .This conforms to the error behavior of the kernel as plotted in Fig. 2a.The other algorithms show minor deviations from each other.

DISCUSSION
We proposed an alternative algorithm for elongated anisotropic Gaussian filters in 2D, which improves throughput and accuracy.Employed in a numerical scheme for estimating fiber directions, namely the maximal response of anisotropic Gaussian filters, it improves accuracy, especially for noisy images with low contrast.Moreover, it outperforms established methods using the structure tensor or the Hessian matrix on synthetic images of fibers.
We present two real-world data sets of sheet molding compounds, to which we apply the method successfully.This application also inspired our experimental setup, for which we generated synthetic images containing straight and parallel fiber bundles.Here, our modifications to the algorithm yield improved precision.Note that the method still performs well on visibly bent carbon fibers.Fig. 7: APPENDIX -Mean angle error for 50 noise images overlayed with synthetic fiber images with direction θ = 0 • , ..., 179 • , and with varying contrast.For each noise image contrast combination, the MAE's maximum over fiber directions was calculated.The mean and standard deviations over 50 noise images are depicted as point symbol with bars for each contrast and algorithm.Note, however, that the standard deviations are small and therefore the bars delimiting the interval are in most cases covered by the symbol for the mean value.
As argued and done byKrause et al. (2010) andWirjadi et al. (2016), we set σ = r.This choice was empirically confirmed byPinter et al. (2018).The parameter ρ is less impacted by the size of the underlying structure.Hence, we use a constant ρ = 6.0, as was found optimal byPinter et al. (2018).
l 2 -error in e −3 between the reconstructed and the true Gaussian kernel with σ 1 = 25, σ 2 = 2, over all angles for the line buffer and the hybrid algorithm with linear resp.cubic interpolation. of the hybrid and line buffer algorithms.

Fig. 3 :
Fig. 3: Visualization of the experimental data set for varying contrast c.
(a) Original slice, gray-values are spread for improved visibility.(b) Result using the MR method.
(c) Result using the structure tensor.(d)Result using the Hessian matrix.

Fig. 6 :
Fig. 6: Direction estimation of SMC with carbon fibers.The directions are color-coded, see color wheel at the top-right side.
(a) Experimental results after applying the median filter.results for the Hessian matrix.

Table 1 :
Complexity per pixel for different algorithms with interpolation.

Table 2 :
l 2 -deviation in 10 −3 between the reconstructed and the true Gaussian kernel.Mean over all angles θ , maximal error in brackets.