Least-squares diffraction imaging using shaping regularization by anisotropic smoothing

Regularization

We use shaping regularization to constrain the model (Fomel, 2007). We penalize edge diffractions by two shaping operators: thresholding and smoothing by anisotropic diffusion. Iterative application of the first operator forces the model to describe the data with the fewest parameters possible (Daubechies et al., 2004) and therefore denoises and emphasizes edge diffraction signatures in the direction perpendicular to the edge. The second operator emphasizes the signatures of edge diffractions along the edge by enforcing their continuity.

Smoothing based on anisotropic diffusion enhances flow-like structures and completes interrupted lines (Weickert, 1998). The input is an unfiltered image, smoothness of which is increasing while diffusion is proceeding in time (Fehmers and Höcker, 2003). In a discrete form anisotropic diffusion can be formulated as follows:

$$\frac{\mathbf{U}_k-\mathbf{U}_{k-1}}{\Delta t} = - \mathbf{D}^T_{im} \mathbf{V} \mathbf{V}^T \mathbf{D}_{im} \mathbf{U}_k\;,$$ (3)

where $\mathbf{U}_k$ is an image at a diffusion time step $k$ with $\Delta t$ sampling interval, $\mathbf{D}_{im}$ is a matrix operator combining PWDs in inline and crossline directions in the image domain and $\mathbf{V}$ is a matrix of eigenvectors $\mathbf{v}$ , which are parallel to linear features' orientations along the dominant local slopes. Classic anisotropic diffusion partial differential equation (Weickert, 1998) requires three-dimensional structure tensor and gradient. In equation 3, we utilize PWD-based structure tensor instead of its three-dimensional counterpart based on derivatives along the coordinate axes and a combination of PWDs in inline and crossline directions in the image domain instead of a three-dimensional gradient. After terms rearrangement, equation 3 takes the form of a linear least-squares solution with forward modeling operator corresponding to the identity matrix:

$$\mathbf{U}_k = \left(\mathbf{I} + \epsilon^2 \mathbf{D}^T_{im} \mathbf{V} \mathbf{V}^T \mathbf{D}_{im}\right)^{-1} \mathbf{U}_{k-1} ,$$ (4)

where $\epsilon$ corresponds to the time step $\Delta t$ . At diffusion time step $k$ we update the previous image $\mathbf{U}_{k-1}$ with $N$ conjugate gradients iterations until the desired smoothness is achieved in the output image $\mathbf{U}_{k}$ . Total number of diffusion time steps $K$ , number of conjugate gradients iterations $N$ and parameter $\epsilon$ control the smoothness of the final result.

Figure 1 shows anisotropic smoothing operator action on a zigzag pattern contaminated with Gaussian noise. Isotropic smoothing operator (see, e.g., Weickert (1998) for details) action is equal for all the directions and thus results in sharpness loss incurred by smoothing across the edges (Figure 1c). Figure 1d shows that anisotropic smoothing operator (equation 4) preserves edges by smoothing along them and thus results in both S/N ratio enhancement and sharpness of the image.

Flow-like coherent noise patterns in Figure 1d are induced by using ``true'' azimuths of edges across the whole image including both signal and noise regions. The result of combining both regularization operators - thresholding and anisotropic smoothing - is shown in Figure 1f: noise and flow-like artifacts are suppressed in regions with no signal. Thus, the combination of thresholding and anisotropic smoothing is required for proper edge diffraction regularization.

In the following examples, edge orientation estimation is done for each location individually based on a structure tensor.

zigzag0SEG,zigzagSEG,diffuseSEG,anisodiffuseSEG,zigzag-thrSEG,anisodiffuse-thrSEG
Figure 1. Zigzag pattern: (a) with no noise; (b) with Gaussian noise added; (c) corresponds to (b) processed by isotropic smoothing operator (notice smoothing across the edges); (d) corresponds to (b) processed by anisotropic smoothing operator (edges are highlighted and preserved). Flow-like coherent noise patterns in (d) are induced by using ``true'' azimuths of edges across the whole image including both signal and noise regions. Thresholding operator action on (b) is shown in (e). The result of combining both regularization operators - anisotropic smoothing and thresholding - is shown in (f). Noise including flow-like artifacts is attenuated.

Least-squares diffraction imaging using shaping regularization by anisotropic smoothing