Elimination of Image Saturation Effects on Multifractal Statistics Using the 2D WTMM Method

Imaging artifacts such as image saturation can restrict the computational analysis of medical images. Multifractal analyses are typically restricted to self-affine, everywhere singular, surfaces. Image saturation regions in these rough surfaces rob them of these core properties, and their exclusion decreases the statistical power of clinical analyses. By adapting the powerful 2D Wavelet Transform Modulus Maxima (WTMM) multifractal method, we developed a strategy where the image can be partitioned according to its localized response to saturated regions. By eliminating the contribution from those saturated regions to the partition function calculations, we show that the estimation of the multifractal statistics can be correctly calculated even with image saturation levels up to 20% (where 20% is the number of saturated pixels over the total number of pixels in the image).

1 Introduction

Instrumentation noise, hardware failures, and inadequate camera or spectrometer calibration can lead to imaging artifacts (Walz-Flannigan et al., 2018; Liu H. et al., 2021). In many circumstances, a quality control procedure eliminates imaging data that do not pass a minimal quality requirement and the samples are re-imaged. However, in certain situations, in particular for imaging data obtained from human samples (e.g., screening mammography), re-imaging is often much more difficult, if not impossible. In other situations, images containing artifacts may still be usable for subjective analyses via visual inspection, but would otherwise be inadequate for objective, computational image analysis pipelines. Therefore, efforts are needed to palliate the deficiencies caused by imaging artifacts on sensitive computational analyses.

Image saturation is a type of distortion where a portion of the acquired image is limited to some maximal pixel value (Figure 1). It can cause problems for computer vision algorithms that assume linearity, unless saturated pixels are identified and handled appropriately (Hasinoff, 2014). Typical approaches to deal with saturated pixels are to either ignore them (eliminate from the analysis) or interpolate their values based on neighboring pixels (Masood et al., 2009). Image saturation also affects calculations in the Fourier domain and needs a strategy to be mitigated (Wetzstein et al., 2010).

FIGURE 1

FIGURE 1. Example of saturation effects on a 360, ×, 360 pixel subregion from a mammogram. (A) Original image. (B) Mesh representation of image (A) as a rough surface. (C) Same image as in (A), but saturated at 20%. (D) Mesh representation of image (C) as a rough surface.

The multifractal analysis of self-affine rough surfaces (Arneodo et al., 2000; Decoster et al., 2000; Arneodo et al., 2003; Wendt et al., 2009) is very valuable in many applications with underlying multiscale non-linear variability, from astrophysics (Khalil et al., 2006; Kestener et al., 2010), geophysics (Roux et al., 2000), texture-based segmentation (Pascal et al., 2018), to biomedical imaging (Kestener et al., 2001; Marin et al., 2017; Gerasimova-Chechkina et al., 2021).

Saturation can easily perturb the multifractal analysis, especially for negative statistical order moments, as discussed in detail below. Indeed, any method that uses either an increment or gradient-based approach to estimate the multifractal signature of rough surfaces (Arneodo et al., 2000; Decoster et al., 2000; Arneodo et al., 2003; Khalil et al., 2006) risks being affected by regions of saturation within the image. Mathematically-speaking, such artifactual images do not possess truly self-affine properties (e.g., everywhere continuous but nowhere differentiable) because these saturation regions are flat and thus infinitely differentiable. Therefore, any estimation of multifractal statistics from such images would be incorrect. To the best of our knowledge, no approach exists in the current literature to rescue the effects of image saturation on multifractal statistics.

The 2D Wavelet Transform Modulus Maxima (WTMM) method has been adapted and applied in three different forms: a multiscale segmentation method (Kestener et al., 2001; Khalil et al., 2007; Roland et al., 2009; Grant et al., 2010; Kestener et al., 2010; McAteer et al., 2010; Batchelder et al., 2014; Marin et al., 2018; Liu J. et al., 2021), a multiscale anisotropy method (Khalil et al., 2006, 2009; Tilbury et al., 2021), and a multifractal formalism (Arneodo et al., 2000; Decoster et al., 2000; Roux et al., 2000; Arneodo et al., 2003; Khalil et al., 2006). The 2D WTMM multifractal method is a multiscale formalism perfectly suited for the analysis of self-affine rough surfaces such as mammograms by identifying density fluctuations and spatial correlations within these surfaces (Batchelder et al., 2014; Plourde et al., 2016; Marin et al., 2017; Gerasimova-Chechkina et al., 2021).

In this manuscript, we demonstrate how the 2D WTMM multifractal method, thanks to the ease with which its space-scale skeleton can be partitioned, allows one to eliminate the effects of saturation on the multifractal statistics. Moreover, we show that the implementation of the strategy used to eliminate these effects does not preclude the analysis of normal (non-saturated) images. In Section 2, a detailed recapitulation of the 2D WTMM multifractal method is presented. The effects of image saturation on the calculation of the multifractal statistics and the rescue method are presented in Section 3, followed by the conclusion and discussion in Section 4.

2 Materials and Methods

2.1 The 2D WTMM Multifractal Method

Let f be a function from ${\mathbb{R}}^2$ into $\mathbb{R}$ and S_h the set of all points x₀ so that the Hölder exponent of f at x₀ is h. The singularity spectrum D(h) of f is the function which associates with any h, the fractal dimension of S_h

$D\left(h\right)=D_F\left\{\mathbf{x}\in {\mathbb{R}}^2\text{\hspace{0.17em},\hspace{0.17em}}h\left(\mathbf{x}\right)=h\right\}.$

The continuous wavelet transform is defined as

${\mathbf{T}}_{\psi }\left[f\right]\left(\mathbf{b},a\right)=\left(\begin{array}{c}\hfill T_{{\psi }_1}\left[f\right]=a^{-2}\int {\psi }_1\left(a^{-1}\left(\mathbf{x}-\mathbf{b}\right)\right)f\left(\mathbf{x}\right)d^2\mathbf{x}\hfill \\ \hfill T_{{\psi }_2}\left[f\right]=a^{-2}\int {\psi }_2\left(a^{-1}\left(\mathbf{x}-\mathbf{b}\right)\right)f\left(\mathbf{x}\right)d^2\mathbf{x}\hfill \end{array}\right)$

where b is the space parameter, $f\in L^2(\mathbb{R})$, a represents scale, x = (x, y), and where the first-order wavelets are defined as (Arneodo et al., 2000)

${\psi }_1\left(x,y\right)=\frac{\partial \mathcal{G}\left(x,y\right)}{\partial x}\text{\hspace{0.17em}and\hspace{0.17em}}{\psi }_2\left(x,y\right)=\frac{\partial \mathcal{G}\left(x,y\right)}{\partial y}$

with

$\mathcal{G}\left(x,y\right)=e^{-\left(x^2+y^2\right)/2}=e^{-|\mathbf{x}{|}^2/2}.$

The wavelet transform can be expressed in terms of its modulus ${\mathcal{M}}_{\psi }[f](\mathbf{b},a)$ and argument ${\mathcal{A}}_{\psi }[f](\mathbf{b},a)$:

${\mathbf{T}}_{\psi }\left[f\right]\left(\mathbf{b},a\right)=\left({\mathcal{M}}_{\psi }\left[f\right]\left(\mathbf{b},a\right),{\mathcal{A}}_{\psi }\left[f\right]\left(\mathbf{b},a\right)\right),(1)$

where

${\mathcal{M}}_{\psi }\left[f\right]\left(\mathbf{b},a\right)=\sqrt{{\left(T_{{\psi }_1}\left[f\right]\left(\mathbf{b},a\right)\right)}^2+{\left(T_{{\psi }_2}\left[f\right]\left(\mathbf{b},a\right)\right)}^2}$

${\mathcal{A}}_{\psi }\left[f\right]\left(\mathbf{b},a\right)=Arg\left(T_{{\psi }_1}\left[f\right]\left(\mathbf{b},a\right)+i{T}_{{\psi }_2}\left[f\right]\left(\mathbf{b},a\right)\right).$

The Wavelet Transform Modulus Maxima (WTMM) are locations b where ${\mathcal{M}}_{\psi }[f](\mathbf{b},a)$ is a local maximum in the angular direction of ${\mathcal{A}}_{\psi }[f](\mathbf{b},a)$ for a given scale a. The WTMM capture the gradient changes in the underlying rough surface. The WTMM are on connected chains, called maxima chains (Arneodo et al., 2000). This process is repeated at every scale. An example for a 2D fractional Brownian motion (fBm) surface (Mandelbrot and Ness, 1968) with Hurst roughness exponent H = 0.7 is provided in Figure 2A, where the black edge detection lines shown in Figures 2B–D are the WTMM chains.

FIGURE 2

FIGURE 2. (A) A 512 × 512-pixel fBm surface with H = 0.7. The maxima chains, i.e., positions where the modulus ${\mathcal{M}}_{\psi }[f](\mathbf{b},a)$ from Eq. 1 is maximum, are shown as black lines while the WTMMM are represented as black dots for scale a = 20 pixels (B), a = 56 pixels (C), and a = 158 pixels (D). Also shown in (D) is the overlay of the maxima chains and WTMMM on the original image smoothed with the Gaussian function $\mathcal{G}$ at scale a = 158 pixels. The WTMMM are then connected across scales to make the WT Skeleton shown in Figure 3.

The WTMM maxima (WTMMM) are defined as the points along the maxima chains where ${\mathcal{M}}_{\psi }[f](\mathbf{b},a)$ is locally maximum (black dots in Figures 2B–D). These WTMMM are connected through scales, a, and form individual maxima lines. The set of all maxima lines is called the Wavelet Transform (WT) space-scale skeleton, as illustrated in Figure 3A.

FIGURE 3

FIGURE 3. (A) WT Skeleton: each vertical line connecting the WTMMM through the scales, by linking the WTMMM closest at the next scale, is known as a maxima line. Each WTMMM point along these maxima lines have a modulus value at each (x, y, log₂(a)) position. (B) Sheaf: One can study the behavior of the modulus as a function of the scale parameter for each maxima line by eliminating the positional information (see Eq. 3). The three vertical purple lines correspond to scales 20, 56, 158 pixels represented in Figures 2B–D. In these two plots, the scale parameter is expressed in units of 7*2^a pixels.

Along a maxima line pointing to the singularity x₀ in the rough surface as a → 0⁺, denoted ${\mathcal{L}}_{{\mathbf{x}}_{\mathbf{0}}}(a)$, the WTMMM follow (Mallat and Hwang, 1992; Arneodo et al., 2000)

${\mathcal{M}}_{\psi }\left[f\right]\left({\mathcal{L}}_{{\mathbf{x}}_{\mathbf{0}}}\left(a\right)\right)\sim a^{h\left({\mathbf{x}}_{\mathbf{0}}\right)},a\to 0^{+}(2)$

where h (x₀) is the Holder roughness exponent. Note that Eq. 2 only holds if the wavelet order (= 1) is greater than the Holder exponent (Mallat and Hwang, 1992; Arneodo et al., 2000), which is a safe assumption to make in this study since all surfaces considered here have roughness exponents less than 1.

This means that h (x₀) can be estimated by considering

$\frac{{\log }_2\left({\mathcal{M}}_{\psi }\left[f\right]\left({\mathcal{L}}_{{\mathbf{x}}_{\mathbf{0}}}\left(a\right)\right)\right)}{{\log }_2\left(a\right)}\sim h\left({\mathbf{x}}_{\mathbf{0}}\right).(3)$

Figure 3B shows such a log-log plot for all of the maxima lines shown in Figure 3A. This is called a sheaf.

2.2 Partition Function and Statistical Order Moments

Let $\mathcal{L}(a)$ be the set of maxima lines at scale a and define the partition functions as (Arneodo et al., 2000, 2003; Khalil et al., 2006):

$Z\left(q,a\right)=\sum _{l\in \mathcal{L}\left(a\right)}{\left(\underset{\left(\mathbf{b},a^{\prime }\right)\in l,a^{\prime }\le a}{\sup }{\mathcal{M}}_{\psi }\left[f\right]\left(\mathbf{b},a^{\prime }\right)\right)}^q(4)$

where q are the statistical order moments. Note that negative q values give more weight to small modulus values while positive q values give more weight to large modulus values. The following power-law relationship allows us to estimate the roughness of a surface (see Arneodo et al. (2003) and references therein):

$Z\left(q,a\right)\sim a^{\tau \left(q\right)},a\to 0^{+}.(5)$

For a monofractal rough surface such as a fBm surface, this τ(q) spectrum is a linear function of q where the slope of τ(q) gives an estimate for the Hurst exponent, H, i.e., (Arneodo et al., 2000):

$\tau \left(q\right)=qH-2.(6)$

However for multifractal rough surfaces, τ(q) is nonlinear which highlights the varying roughness exponents in the underlying surface (Decoster et al., 2000).

A Legendre transform can be applied to the τ(q) spectrum to obtain the D(h) spectrum of singularities (Arneodo et al., 2000, 2003; Khalil et al., 2006):

$D\left(h\right)=\underset{q}{\min }\left(qh-\tau \left(q\right)\right).(7)$

Given the numerical impediments related to the Legendre transform (Arneodo et al., 1995), one can avoid directly performing the Legendre transform by considering h and D(h) as mean quantities defined in a canonical ensemble, i.e., with respect to their Boltzmann weights computed from the WTMMM (Arneodo et al., 2000):

$W_{\psi }\left[f\right]\left(q,l,a\right)=\frac{{\left|{\sup }_{\left(\mathbf{b},a^{\prime }\right)\in l,a^{\prime }\le a}{\mathcal{M}}_{\psi }\left[f\right]\left(\mathbf{b},a^{\prime }\right)\right|}^q}{Z\left(q,a\right)}$

where Z (q, a) was defined in Eq. 4. One can then compute the expectation values:

$h\left(q,a\right)=\sum _{l\in \mathcal{L}\left(a\right)}\ln \left|\underset{\left(\mathbf{b},a^{\prime }\right)\in l,a^{\prime }\le a}{\sup }{\mathcal{M}}_{\psi }\left[f\right]\left(\mathbf{b},a^{\prime }\right)\right|W_{\psi }\left[f\right]\left(q,l,a\right)(8)$

and

$D\left(q,a\right)=\sum _{l\in \mathcal{L}\left(a\right)}{W}_{\psi }\left[f\right]\left(q,l,a\right)\ln W_{\psi }\left[f\right]\left(q,l,a\right).(9)$

This gives the following:

$h\left(q\right)=\frac{d\tau \left(q\right)}{dq}=\underset{a\to 0^{+}}{\lim }\frac{h\left(q,a\right)}{\ln a}(10)$

$D\left(q\right)=\underset{a\to 0^{+}}{\lim }\frac{D\left(q,a\right)}{\ln a},(11)$

from which we obtain the D(h) singularity spectrum. Numerical calculations following the steps outlined above are presented in Figure 4 for sample fBm surfaces with H = 0.1, 0.3, 0.5, 0.7.

FIGURE 4

FIGURE 4. (A) h (a, q) curves (Eq. 8) and (B) D (a, q) curves (Eq. 9) for q values from −3 to 5 for fBm surfaces with H = 0.5. (C) h(q) (Eq. 10) for fBms with H = 0.1, 0.3, 0.5, 0.7. (D) List of q values used. Corresponding (E) D(q) (Eq. 11) and (F) D(h).

Statistical order moments (q values) play a crucial role in the 2D WTMM multifractal method. These values allow one to emphasize different singularity strengths of the underlying surface by weighing the modulus of the wavelet transform along the maxima lines. The quantity and quality of the underlying data will determine the range of the available statistical order moments (for an in-depth discussion, see Khalil et al., 2006).

2.3 Numerical Implementation

We followed previously established numerical calibration studies using the 2D WTMM multifractal method (Arneodo et al., 2000; Decoster et al., 2000; Khalil et al., 2006). For each Hurst value explored, 32 synthetic fBm surfaces 1,024 × 1,024 pixels in size were generated using a Fourier filtering synthesis algorithm (see Arneodo et al. (2000) and references therein). The results presented in this manuscript correspond to the average partition functions over these sets of 32 surfaces.

Saturation was introduced by considering the cumulative density distribution of pixel values for a given image and determining the critical pixel value corresponding to the saturation level desired. Then all pixel values above that critical value were changed to the maximal value.

The calculations of h(q) and D(q) were obtained from the slope of the linear fit of the log-log representation of the h (q, a) (Eq. 10) and D (q, a) (Eq. 11) curves, respectively. In this manuscript, a fixed range of scales was used for fBm surfaces of all H values and all saturation levels: from 17 to 56 pixels.

The numerical calculations described in this section were performed using Xsmurf, a Tcl\Tk software package that runs C-based routines (https://github.com/pkestene/xsmurf).

3 Results

3.1 Effects of Saturation on WTMM Statistics

The apparition of extraneous maxima lines in the sheafs of saturated fBm images are displayed in Figure 5. Compared to the original, non-saturated fBm surfaces, two categories of new maxima lines appear in higher numbers with increasing saturation: 1) maxima lines associated with extremely low modulus values, and 2) other maxima lines with an overall negative slope.

FIGURE 5

FIGURE 5. (A) A sample fBm surface with H = 0.7 and its corresponding sheaf (right). (B) The same fBm surface as in (A) but saturated at 1% and the corresponding sheaf directly below. (C–E) Saturated fBms at (C) 5%, (D) 10%, and (E) 20% and sheafs below. Maxima lines associated with extremely low modulus values, and others with an overall negative slope, appear in higher numbers with increasing saturation. These extraneous maxima lines will interfere with h (a, q) and D (a, q) calculations, which then reduces the number of statistical order moments as shown Figure 6.

Keeping these maxima lines in the skeleton that is fed to the partition function calculation (Eq. 4) affects the behavior of the h (a, q) and D (a, q) curves for a range of q values, as shown in Figures 6A,B. The non-linear behavior of log-log representation of the h (a, q) and D (a, q) curves for several q values precludes the calculation of h(q) (Eq. 10) and D(q) (Eq. 11), as shown in Figures 6C,D.

FIGURE 6

FIGURE 6. (A) h (a, q), (B) D (a, q), (C) h(q), and (D) D(q) for H = 0.7 with 20% saturation. The colors corresponds to the different q values as described in the legend in Figure 4D. Image saturation reduces the range of statistical order moments because of the lack of a good power law fit in the h (a, q) and D (a, q) plots. Thus in the corresponding h(q) and D(q) plots there are less statistics for saturated images (black dots). Similar plots are shown in the Supplementary Material for saturation levels 1, 5, and 10%.

As a consequence, whereas a range of −3 ≤ q ≤ 5 was adequate for the non-saturated fBms (Figure 4), this range becomes limited to $\sim 0<q\le 5$ for the saturated fBms of H = 0.7, regardless of the saturation level (1, 5, 10, 20%).

3.2 The Rescue Method on Saturated Surfaces

We propose an approach to eliminate the two categories of extraneous maxima lines that are responsible for the bad power-law behavior of the h (a, q) and D (a, q) curves. As discussed below, this technique is referred to as the rescue method. In short, it consists of eliminating a subset of maxima lines by applying a first threshold filter on the modulus value at scale a = 1 (referred to below as MF), and a second adjusted slope threshold filter (referred to below as SF).

Let $l\in \mathcal{L}(a)$ be a maxima line ending at scale a_max and define ${\mathcal{M}}_{l(a=a^{\prime })}$ as the value of the modulus for maxima line l at scale a′. The modulus filter (MF) is used to eliminate the maxima lines that have a low modulus value at scale a = 1, i.e., when ${\mathcal{M}}_{l(a=1)}\le MF$.

Next, we define the adjusted slope, m_adj, as:

$m_{\mathit{adj}}=\frac{{\log }_2\left({\mathcal{M}}_{l\left(a=1\right)}\right)-{\log }_2\left({\mathcal{M}}_{l\left(a=a_{max}\right)}\right)}{{\log }_2\left({\mathcal{M}}_{l\left(a=a_{max}\right)}\right)}.$

Therefore, to calculate m_adj, one takes the starting modulus value of a maxima line (at scale a = 1) and subtracts the ending modulus value (at scale a = a_max), all divided by the ending modulus value. For a given slope filter threshold, SF, the subset of maxima lines that satisfy m_adj ≥ SF, are eliminated.

3.2.1 Determination of Numerical Filter Parameters SF and MF

When comparing sheafs from saturated surfaces to sheafs from non-saturated surfaces, we empirically determined that any maxima lines with ${\mathcal{M}}_{l(a=1)}$ of 16 (2⁴) or less should be removed, i.e., log₂ (MF) = 4.

The data mining approach that we took to determine the optimal value for SF is explored in Figure 7 for a fBm surface with H = 0.7 with 20% saturation, where a fixed value of log₂ (MF) = 4 was used, with SF = 0.2, 0.5, 0.8 in the left (Figures 7A,D,G), center (Figures 7B,E,H), and right (Figures 7C,F,I), respectively. Our goal was to explore different values of SF that would remove the maxima lines that caused the saturated sheaf to look different from an unsaturated sheaf. For instance, in Figure 7G (SF = 0.2), the resulting filtered saturated sheaf has maxima lines that are truncated (too many good maxima lines were eliminated) while Figure 7I has many negative sloping, low valued maxima lines compared to the unsaturated case (not enough maxima lines were eliminated). Figure 7H displays what was determined to be the optimal filtering coefficients. Figure 7J is the sheaf corresponding to the unsaturated fBm surface. Note that although only three different values of SF and a single value of MF are explored in Figure 7, we numerically explored a much larger collection of values for these parameters (data not shown).

FIGURE 7

FIGURE 7. Determination of the numerical values for SF and MF for a fBm surface with H = 0.7 and saturation level of 20%. Maxima lines that were eliminated due to the modulus filter $({\mathcal{M}}_{l(a=1)}\le MF)$, are shown in red, while maxima lines that were eliminated due to the adjusted slope filter (m_adj ≥ SF) are shown in green. The remaining maxima lines are shown in black. Here a fixed modulus filter value of log₂ (MF) = 4 is used while the left column has SF = 0.2, the middle column has SF = 0.5, and the right column has SF = 0.8. In (A–C) the complete original sheaf is shown with the eliminated and kept maxima lines, in (D–F) only the eliminated maxima lines are shown, and in (G–I) only the kept maxima lines are shown. As a comparison, the sheaf from the unsaturated surface is shown (J).

A sample fBm surface with H = 0.7 at 20% saturation level is shown in Figures 8A,B. Also shown are the maxima chains, color-coded according to their elimination process (red based on the MF and green based on the SF) in Figure 8C. And Figure 8D shows the corresponding skeleton of eliminated maxima lines.

FIGURE 8

FIGURE 8. (A) A sample fBm surface with H = 0.7. (B) The same surface as in (A) but saturated at 20%. (C) The saturated surface overlaid with maxima chains (at scale 20 pixels) that were kept (black), those that were eliminated by applying the modulus filter log₂ (MF) = 4 (red), and those eliminated by applying the slope filter SF = 0.5 (green). (D) The corresponding extraneous maxima lines that were removed.

3.2.2 Pruned Sheafs

Using the filtering parameters selected through the data mining approach outlined above, (log₂ (MF) = 4 and SF = 0.5) these resulting “pruned” sheafs were then passed to the partition functions calculations (Eq. 4). From these, we get h(a,q) (Eq. 8) and D(a,q) (Eq. 9) shown in Figures 9A,B, respectively, where we can see the lack of stair step behavior compared to their saturated counterparts in Figures 6A,B and Supplementary Figures S1–S3. With these now robust h(a,q) and D(a,q) curves, from which a reliable power law fit can be obtained, we are able to expand the range of statistical order moments for the calculation of h(q) and D(q), as displayed in Figures 9C,D. For reference, the saturated images had a q range of $\sim 0<q\le 5$ while with the rescue method the range was expanded to −2 ≤ q < 5.

FIGURE 9

FIGURE 9. Efficacy of the rescue method. The h (a, q) (A) and D (a, q) (B) plots have better power law fits for a wider range of statistical order moments relative to their saturated counterparts. This expanded range of q values is reflected in the h(q) (C) and D(q) (D) plots. The saturated images had a q range 0.1 < q < 5 while the rescue method had −2 ≤ q ≤ 5.

The low modulus maxima lines (shown in red in Figures 7 and 8) are likely responsible for the step size behavior of the h (a, q) curves for negative q value shown in Figures 6 and Supplementary Figures S1–S3. However, Table 1 shows that only eliminating these maxima lines would be insufficient to eliminate the effects of saturation. Indeed, the number of maxima lines removed depends on the saturation level, which also dictates which of the two filtering processes removes more lines. For example, for the set of fBm surfaces with H = 0.7 and saturation level 1%, using log₂ (MF) = 4 and SF = 0.5, the average percentage of removed maxima lines was 0.4% due to MF and 4.8% due to SF. For saturation level 5%, these averages were 3.9% due to MF and 4.9% due to SF. For saturation level 10%, more maxima lines were removed due to MF (9.3%) than due to SF (4.9%). And for saturation level 20%, the average percentages of removed maxima lines were 21.8% due to MF and 4.6% due to SF. The percentages of removed maxima lines by filtering method for each Hurst value and for each saturation level are listed in Table 1.

TABLE 1

TABLE 1. Percentages of maxima lines removed by filtering processes.

3.3 The Rescue Method on Unsaturated Surfaces

We used the filtering parameters empirically determined above (log₂ (MF) = 4 and SF = 0.5) on the sheaf of a non-saturated fBm surface with H = 0.7. There were only very minor differences between the h(q) and D(q) curves from the unsaturated surfaces with non-filtered sheafs (normal) vs. the unsaturated surfaces with filtered sheafs (treated) (Figure 10).

FIGURE 10

FIGURE 10. The h(q) (A) and D(q) (B) plots show that the rescue method does not fundamentally alter the behavior of the statistics when applied to unsaturated images.

We note that the rescue method will limit the range of q values only if the underlying surface was saturated. This is to be expected as it is due to the removal of maxima lines from the sheaf that is fed to the partition function. Therefore, one should only compare the range of q values of a pruned sheaf to that of a smaller unsaturated image, for which the range of q values will be smaller (see Khalil et al. (2006)).

4 Discussion

The multifractal analysis of rough surfaces inherently assumes that these surfaces are scale-invariant and everywhere singular. These conditions are not met when image saturation regions are present. The rescue method proposed here takes full advantage of the space-scale partitioning capability that is intrinsic to the 2D WTMM multifractal method, which allows for the elimination of a subset of the wavelet transform skeleton corresponding to saturated regions.

The concept of pruning a sheaf was first introduced in Kestener et al. (2001) and then further explored in Batchelder et al. (2014). However, pruning a sheaf to mitigate image saturation effects had yet to be explored. As a contribution to the medical imaging field, data with image artifacts that would have otherwise been rejected could now potentially be rescued and included for multifractal statistical analyses.

Although several applications may benefit from this rescue method, as discussed in the introduction, an immediate benefit is for the sliding-window analysis of mammograms (Marin et al., 2017; Gerasimova-Chechkina et al., 2021). To offer the best possible computational aid to the radiologists in their interpretation of these mammograms, it is critical to extract as much quantitative information as possible from each subregion of each mammogram. A key to this success lies in reducing the number of subregions that have to be rejected, in particular due to image saturation, for which the implementation of this rescue method will be crucial.

4.1 Limitations and Future Explorations

In this study, only positive values of the Hurst exponent were considered for the calibration fBm surfaces (H = 0.1, 0.3, 0.5, 0.7), which is motivated, and justified, by the Hurst values measured in mammogram subregions (Marin et al., 2017; Gerasimova-Chechkina et al., 2021). Despite only using four Hurst values, the effects are similar, which precludes the need to investigate smaller intervals. We have no reason to suspect that exploring other values, for example H = 0.2, 0.4, 0.6, would produce results that require a substantially different approach. However, future efforts should be undertaken to investigate the effects of saturation on surfaces with −1 < H < 0. It is likely that the slope filter (SF) would have to be calibrated. One could also expand the investigation on multifractal surfaces (Decoster et al., 2000; Roux et al., 2000).

The numerical determination of MF and SF was done subjectively. We observe that a constant value of SF = 0.5 worked well for all Hurst values considered in this study, whereas MF is strongly correlated with the Hurst exponent. A future effort could be to implement objective approaches. For example, an automated determination of these two thresholds could be implemented based on the calculation of outliers in the probability density functions of $\log (\mathcal{M})$ and m_adj.

Another improvement to the approach could be to use two thresholds instead of one for SF (e.g., SF_high and SF_low). Maxima lines for which m_adj ≥ SF_high would be eliminated, those for which m_adj < SF_low would be kept, and those in between would require an additional classification procedure. And finally, there may be methods from machine learning to use in order to improve the maxima lines classification, assuming the unsaturated image has scale-invariant properties.

We note that the calculations of the Hurst values reported in this manuscript are underestimates of the theoretical values. This is a known phenomenon for 1D fBm signals (Muzy et al., 1991, 1994; Arneodo et al., 1995) and 2D fBm surfaces (Arneodo et al., 2000). The underestimates reported here are perhaps slightly more accentuated than previous works (Arneodo et al., 2000; Khalil et al., 2006). This could be due to the fact that a fixed universal range of scales was used for every set of fBm surfaces and all saturation levels in this study, as mentioned in Section 2.3.

The rescue method presented in this manuscript could be used on any medical imaging where the 2D WTMM multifractal method is applicable, such as lung CT scans, or histology slides (Khalil et al., 2009). It is also likely that the rescue method could be integrated into the 1D WTMM method, and therefore applicable to 1D signals such as electro-encephalograms (Richard et al., 2015) or thermography time series (Gerasimova et al., 2013, 2014), for which the MF and SF values would need to be fine-tuned for the application.

Data Availability Statement

The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

Author Contributions

Study design: JJ and AK. Image analysis and statistical analysis: JJ and AK. Paper written by JJ and AK.

Funding

Research reported in this publication was supported by the National Cancer Institute of the National Institutes of Health under award number R15CA246335 (approximately $12,000). The content is solely the responsibility of the authors and does not necessarily represent the official views of the National Institutes of Health.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Acknowledgments

We thank David Bradley and Peter Stechlinski as well as the CompuMAINE Lab members for helpful technical discussions.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fphys.2022.921869/full#supplementary-material

Supplementary Figure S1 | (A) h(a, q), (B) D(a, q), (C) h(q), and (D) D(q) for H = 0.7 with 1% saturation. The colors corresponds to the different q values as described in the legend in Figure 4D.

Supplementary Figure S2 | (A) h(a, q), (B) D(a, q), (C) h(q), and (D) D(q) for H = 0.7 with 5% saturation. The colors corresponds to the different q values as described in the legend in Figure 4D.

Supplementary Figure S3 | (A) h(a, q), (B) D(a, q), (C) h(q), and (D) D(q) for H = 0.7 with 10% saturation. The colors corresponds to the different q values as described in the legend in Figure 4D.

References

Arneodo A., Bacry E., Muzy J. F. (1995). The Thermodynamics of Fractals Revisited with Wavelets. Phys. A Stat. Mech. its Appl. 213, 232–275. doi:10.1016/0378-4371(94)00163-N

CrossRef Full Text | Google Scholar

Arnéodo A., Decoster N., Kestener P., Roux S. G. (2003). A Wavelet-Based Method for Multifractal Image Analysis: From Theoretical Concepts to Experimental Applications. Adv. Imaging Electr. Phys. 126, 1–92. doi:10.1016/s1076-5670(03)80014-9

CrossRef Full Text | Google Scholar

Arnéodo A., Decoster N., Roux S. G. (2000). A Wavelet-Based Method for Multifractal Image Analysis. I. Methodology and Test Applications on Isotropic and Anisotropic Random Rough Surfaces. Eur. Phys. J. B 15, 567–600. doi:10.1007/s100510051161

CrossRef Full Text | Google Scholar

Batchelder K. A., Tanenbaum A. B., Albert S., Guimond L., Kestener P., Arneodo A., et al. (2014). Wavelet-based 3d Reconstruction of Microcalcification Clusters from Two Mammographic Views: New Evidence that Fractal Tumors Are Malignant and Euclidean Tumors Are Benign. PloS one 9, e107580. doi:10.1371/journal.pone.0107580

PubMed Abstract | CrossRef Full Text | Google Scholar

Decoster N., Roux S., Arneodo A. (2000). A Wavelet-Based Method for Multifractal Image Analysis. II. Applications to Synthetic Multifractal Rough Surfaces. Eur. Phys. J. B-Condensed Matter Complex Syst. 15, 739–764. doi:10.1007/s100510051179

CrossRef Full Text | Google Scholar

Gerasimova E., Audit B., Roux S. G., Khalil A., Argoul F., Naimark O., et al. (2013). Multifractal Analysis of Dynamic Infrared Imaging of Breast Cancer. EPL 104, 68001. doi:10.1209/0295-5075/104/68001

CrossRef Full Text | Google Scholar

Gerasimova E., Audit B., Roux S. G., Khalil A., Gileva O., Argoul F. o., et al. (2014). Wavelet-based Multifractal Analysis of Dynamic Infrared Thermograms to Assist in Early Breast Cancer Diagnosis. Front. Physiol. 5, 176. doi:10.3389/fphys.2014.00176

PubMed Abstract | CrossRef Full Text | Google Scholar

Gerasimova-Chechkina E., Toner B. C., Batchelder K. A., White B., Freynd G., Antipev I., et al. (2021). Loss of Mammographic Tissue Homeostasis in Invasive Lobular and Ductal Breast Carcinomas vs. Benign Lesions. Front. Physiol. 12. doi:10.3389/fphys.2021.660883

CrossRef Full Text | Google Scholar

Grant J., Verrill C., Coustham V., Arneodo A., Palladino F., Monier K., et al. (2010). Perinuclear Distribution of Heterochromatin in Developing c. elegans Embryos. Chromosome Res. 18, 873–885. doi:10.1007/s10577-010-9175-2

PubMed Abstract | CrossRef Full Text | Google Scholar

Hasinoff S. W. (2014). Saturation (Imaging). Boston, MA: Springer US, 699–701. doi:10.1007/978-0-387-31439-6-483

CrossRef Full Text | Google Scholar

Kestener P., Conlon P. A., Khalil A., Fennell L., McAteer R. T. J., Gallagher P. T., et al. (2010). Characterizing Complexity in Solar Magnetogram Data Using a Wavelet-Based Segmentation Method. ApJ 717, 995–1005. doi:10.1088/0004-637x/717/2/995

CrossRef Full Text | Google Scholar

Kestener P., Lina J., St-Jean P., Arneodo A. (2001). Wavelet-based Multifractal Formalism to Assist in Diagnosis in Digitized Mammograms. Image Anal. Stereol. 20, 169–174.

Google Scholar

Khalil A., Aponte C., Zhang R., Davisson T., Dickey I., Engelman D., et al. (2009). Image Analysis of Soft-Tissue In-Growth and Attachment into Highly Porous Alumina Ceramic Foam Metals. Med. Eng. Phys. 31, 775–783. doi:10.1016/j.medengphy.2009.02.007

PubMed Abstract | CrossRef Full Text | Google Scholar

Khalil A., Grant J. L., Caddle L. B., Atzema E., Mills K. D., Arneodo A. (2007). Chromosome Territories Have a Highly Nonspherical Morphology and Nonrandom Positioning. Chromosome Res. 15, 899–916. doi:10.1007/s10577-007-1172-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Khalil A., Joncas G., Nekka F., Kestener P., Arnéodo A. (2006). Morphological Analysis of H I Features. II. Wavelet‐based Multifractal Formalism. Astrophys. J. Suppl. S 165, 512–550. doi:10.1086/505144

CrossRef Full Text | Google Scholar

Liu H., Cao S., Ling Y., Gan Y. (2021a). Inpainting for Saturation Artifacts in Optical Coherence Tomography Using Dictionary-Based Sparse Representation. IEEE Photonics J. 13, 1–10. doi:10.1109/jphot.2021.3056574

CrossRef Full Text | Google Scholar

Liu J., Enderlin E. M., Marshall H.-P., Khalil A. (2021b). Automated Detection of Marine Glacier Calving Fronts Using the 2-d Wavelet Transform Modulus Maxima Segmentation Method. IEEE Trans. Geosci. Remote Sens. 59, 9047–9056. doi:10.1109/TGRS.2021.3053235

CrossRef Full Text | Google Scholar

Mallat S., Hwang W. L. (1992). Singularity Detection and Processing with Wavelets. IEEE Trans. Inf. Theory 38, 617–643. doi:10.1109/18.119727

CrossRef Full Text | Google Scholar

Mandelbrot B. B., Van Ness J. W. (1968). Fractional Brownian Motions, Fractional Noises and Applications. SIAM Rev. 10, 422–437. doi:10.1137/1010093

CrossRef Full Text | Google Scholar

Marin Z., Batchelder K. A., Toner B. C., Guimond L., Gerasimova-Chechkina E., Harrow A. R., et al. (2017). Mammographic Evidence of Microenvironment Changes in Tumorous Breasts. Med. Phys. 44, 1324–1336. doi:10.1002/mp.12120

PubMed Abstract | CrossRef Full Text | Google Scholar

Marin Z., Wallace J. K., Nadeau J., Khalil A. (2018). Wavelet-based Tracking of Bacteria in Unreconstructed off-axis Holograms. Methods 136, 60–65. doi:10.1016/j.ymeth.2017.09.003

PubMed Abstract | CrossRef Full Text | Google Scholar

Masood S. Z., Zhu J., Tappen M. F. (2009). Automatic Correction of Saturated Regions in Photographs Using Cross-Channel Correlation. Comput. Graph. Forum 28, 1861–1869. doi:10.1111/j.1467-8659.2009.01564.x

CrossRef Full Text | Google Scholar

McAteer R. T. J., Kestener P., Arneodo A., Khalil A. (2010). Automated Detection of Coronal Loops Using a Wavelet Transform Modulus Maxima Method. Sol. Phys. 262, 387–397. doi:10.1007/s11207-010-9530-7

CrossRef Full Text | Google Scholar

Muzy J. F., Bacry E., Arneodo A. (1994). The Multifractal Formalism Revisited with Wavelets. Int. J. Bifurc. Chaos 04, 245–302. doi:10.1142/s0218127494000204

CrossRef Full Text | Google Scholar

Muzy J. F., Bacry E., Arneodo A. (1991). Wavelets and Multifractal Formalism for Singular Signals: Application to Turbulence Data. Phys. Rev. Lett. 67, 3515–3518. doi:10.1103/physrevlett.67.3515

PubMed Abstract | CrossRef Full Text | Google Scholar

Pascal B., Pustelnik N., Abry P., Serres M., Vidal V. (2018). “Joint Estimation of Local Variance and Local Regularity for Texture Segmentation. Application to Multiphase Flow Characterization,” in 2018 25th IEEE International Conference on Image Processing (Athens, Greece: ICIP), 2092–2096. doi:10.1109/ICIP.2018.8451380

CrossRef Full Text | Google Scholar

Plourde S. M., Marin Z., Smith Z. R., Toner B. C., Batchelder K. A., Khalil A. (2016). Computational Growth Model of Breast Microcalcification Clusters in Simulated Mammographic Environments. Comput. Biol. Med. 76, 7–13. doi:10.1016/j.compbiomed.2016.06.020

PubMed Abstract | CrossRef Full Text | Google Scholar

Richard C. D., Tanenbaum A., Audit B., Arneodo A., Khalil A., Frankel W. N. (2015). Swdreader: a Wavelet-Based Algorithm Using Spectral Phase to Characterize Spike-Wave Morphological Variation in Genetic Models of Absence Epilepsy. J. Neurosci. Methods 242, 127–140. doi:10.1016/j.jneumeth.2014.12.016

PubMed Abstract | CrossRef Full Text | Google Scholar

Roland T., Khalil A., Tanenbaum A., Berguiga L., Delichère P., Bonneviot L., et al. (2009). Revisiting the Physical Processes of Vapodeposited Thin Gold Films on Chemically Modified Glass by Atomic Force and Surface Plasmon Microscopies. Surf. Sci. 603, 3307–3320. doi:10.1016/j.susc.2009.09.021

CrossRef Full Text | Google Scholar

Roux S. G., Arnéodo A., Decoster N. (2000). A Wavelet-Based Method for Multifractal Image Analysis. III. Applications to High-Resolution Satellite Images of Cloud Structure. Eur. Phys. J. B 15, 765–786. doi:10.1007/s100510051180

CrossRef Full Text | Google Scholar

Tilbury K., Han X., Brooks P., Khalil A. (2021). Multiscale Anisotropy Analysis of Second-Harmonic Generation Collagen Imaging of Mouse Skin. J. Biomed. Opt. 26, 065002. doi:10.1117/1.jbo.26.6.065002

CrossRef Full Text | Google Scholar

Walz-Flannigan A. I., Brossoit K. J., Magnuson D. J., Schueler B. A. (2018). Pictorial Review of Digital Radiography Artifacts. RadioGraphics 38, 833–846. doi:10.1148/rg.2018170038

PubMed Abstract | CrossRef Full Text | Google Scholar

Wendt H., Roux S. G., Jaffard S., Abry P. (2009). Wavelet Leaders and Bootstrap for Multifractal Analysis of Images. Signal Process. 89, 1100–1114. doi:10.1016/j.sigpro.2008.12.015

CrossRef Full Text | Google Scholar

Wetzstein G., Ihrke I., Heidrich W. (2010). “Sensor Saturation in Fourier Multiplexed Imaging,” in IEEE Computer Society Conference on Computer Vision and Pattern Recognition, 545–552. doi:10.1109/CVPR.2010.5540166

CrossRef Full Text | Google Scholar