-diffusiontensors

This command produces tensor fields from noisy input images. These encode smoothing geometries for the corrupted images, each field having the same dimensions as the input image. The tensors comprising each field are in a one-to-one correspondence with the pixels of the input image.

The command implements the "anisotropic smoothing" operation usually attributed to -smooth, one that strives to recover corrupted images (see the example on the left). Typically, -smooth calls the internal implementation of -diffusiontensors when it has not otherwise been given a particular smoothing geometry.

A script author usually calls -diffusiontensors when she wishes to modify smoothing geometries in particular ways. See “Tensors for the Tonsorially Challenged” in the Cookbook for an example.

Belying a close relationship, the parameter list for -diffusiontensors echoes many of -smooth's. Indeed, -smooth essentially passes such parameters unchanged to its internal copy of -diffusiontensors. The command has the following format:

$gmic -diffusiontensors sharpness>=0,0<=anisotropy<=1, alpha[%]>=0,sigma[%]>=0,is_sqrt={0|1}

  1. sharpness: >= 0.0 Optional and defaults to 0.7. Increasing values diminish smoothing and preserves small scale detail. Settings in excess of 1.0 run the risk of suppressing smoothing too much, preserving noise as well as detail.

  2. anisotropy: range: 1.0 ≥ a ≥ 0.0. Optional and defaults to 0.3. Induces greater smoothing along edge contours at the expense of cross-gradient smoothing. Smaller anisotropy values induce a shift toward unidirectional diffusion; the effect tends to be more like conventional blurring. Anisotropy values near unity often causes "false edges" to be detected creating sinuous, hairlike artifacts.

  3. alpha: Pre-structure blur variance. Gaussian variance σ0 ≥ 0. Optional and defaults to 0.6. The command initially calls -structuretensors to locate edges in images. Before doing so, -diffusiontensors applies a conventional isotropic Gaussian blur to the operand image, which, to a degree, assists -structuretensors with edge detection and finding edge orientation. Alpha sets the Gaussian variance, σ, for that operation. As it is based on isotropic blurring, too high a value of alpha can obliterate edges as well as noise.

  4. sigma: Post-structure blur filter. Gaussian variance σ1 ≥ 0. Optional and defaults to 1.1. Following its call to -structuretensors, -diffusiontensors blurs the tensor field it obtains from that command; this parameter sets the variance, σ. For particularly noisy images, the tensor field produced by -structuretensors can have spurious entries; blurring tends to reduce the effects of such outliers. Larger values of sigma tend to favor large-scale laminar smoothing over localized 'eddy' smoothing.

  5. is_sqrt: Boolean flag {0|1} Modifies the computation of the smoothing intensities in the gradient and contour directions. See "Smoothing Geometry," below, for further details.

Smoothing Geometry

For each pixel in a noisy image, -diffusiontensors calculates a tensor embodying a "smoothing geometry" for that pixel. This tensor encodes two vectors, one directed along the gradient, the largest intensity change manifest in the locale of a pixel, and which is likely to be an edge, and the other along the contour of the (suspected) edge. These vectors are orthogonal to each other.

In the presence of noise, and with increased anisotropy settings, the smoothing vector along the contour is stronger than across the gradient. The operation tends to strengthen edges because diffusion proceeds more rapidly in directions parallel to edges, smearing noise pixels without smearing edge gradients.

The process is adaptive. Where intensity differences are great in a particular locale, indicating edges much stronger than the ambient noise, smoothing vectors are correspondingly small, an inverse proportional relationship. Strong edges do not require as much enhancement as degraded ones, and unnecessary smoothing would require more time and possibly introduce artifacts.

-diffusiontensors calculates smoothing geometries from two sources: guidance parameters from the command line, in particular, sharpness, anisotropy, and is_sqrt, and a survey performed by -structuretensors. The survey identifies "edge-like" structures from differences in image intensity in the locality of each pixel.

-structuretensors itself produces a tensor field. The tensors comprising that field characterize the rapidity of change in image intensity along the cardinal directions (width and height) in the neighborhood of associated pixels. Pixels which straddle marked changes in intensity give rise to structure tensors with eigenvalues that differ greatly in value, their associated eigenvectors oriented along and transverse to the gradient.

-diffusiontensors sets out to estimate corresponding smoothing tensors for each pixel, applying sharpness, anisotropy, is_root and the corresponding eigenvalues and eigenvectors computed from the structure tensors. At each pixel, smoothing tensors have the same orientation as the corresponding structure tensor but have different magnitudes, respectively s+ for smoothing strength in the gradient direction and s‒ for smoothing strength in the contour direction.

These new magnitudes have inverse power relationships with the corresponding eigenvalues of the structure tensors. Smoothing is very small where structure tensor eigenvalues exhibit large differences. This phenomenon corresponds to well-defined edges in the source image and these are not regions where one needs a great deal of edge enhancement. As the absolute sizes and differences of structure tensor eigenvalues decline, gradient and contour smoothing magnitudes increase. This is natural, as the smaller eigenvalues correspond to less distinct features which smoothing might profitably enhance.

We could regard s+ and s‒ as the axes of an elliptic diffusion kernel; differences in their magnitudes induces ellipses of greater eccentricity and a correspondingly greater directional effect in the diffusion pattern. The command line parameters which set the overall smoothing magnitude and effect a differentiation between s+ and s‒ are, respectively, sharpness and anisotropy.

Determining Smoothing Strength

Given are:

  1. λ+, the gradient eigenvalue, derived for each pixel in the source image by -structuretensors. Its corresponding eigenvector points along the steepest gradient.
  2. λ‒, the contour eigenvalue, corresponding to the contour vector, oriented at right angles to the gradient vector and pointing along the edge.
  3. sharpness, from the command line, which scales exponents p1 and p2 (see relationships following) and, correspondingly, uniformly reduces smoothing values s+ and s‒.
  4. anisotropy, from the command line. As it ranges over the semi-closed interval [0,…1), this factor increases the inequality p1 ≤ p2 and which, from the “Contour and Gradient Smoothing” relationships following, induces s+ ≤ s‒. That is, as anisotropy increases, smoothing along the gradient, s+, decreases in relation to smoothing along the contour, s‒. The eccentricity of diffusion kernels increase and smoothing becomes increasingly biased in the direction set by contour eigenvectors.
  5. is_sqrt, when set, sharpness is scaled by one half. Since this value sets p1, and p1 in turn scales p2, both these exponents are scaled by one half, and, through the exponential law, induces square root operations on the calculations of s+ and s‒. As with an inverse relationship, the halving of exponents p1 and p2 strengthen s+ and s‒. The effect of setting this flag to one roughly amounts to “super-smoothing.”

With these given parameters, the following relationships calculate smoothness values s+ and s‒:

Exponents
p 1 = { sharpness 2 ; is_sqrt > 0 sharpness ; is_sqrt = 0 } p_1 = left lbrace{ binom{ ital "sharpness" over 2; ital "is_sqrt" > 0}{ ital "sharpness"; ital "is_sqrt" = 0} }right rbrace p 2 = p 1 1 − anisotropy p_2 = p_1 over { 1 - ital "anisotropy" }

Contour and Gradient Smoothing
s - = 1 ( 1 + λ + + λ - ) p 1 s_"-" = 1 over {(1 + %lambda_"+" + %lambda_"-")^{p_1} } s + = 1 ( 1 + λ + + λ - ) p 2 s_"+" = 1 over {(1 + %lambda_"+" + %lambda_"-")^{p_2} }

A Diffusiontensors Pipeline

To help develop a concrete sense of how -diffusiontensors operates, we trace histories of four pixels through an image reconstruction pipeline mediated by this command and -smooth. Apart from the action of -smooth, this narrative is internal to the -diffusiontensors command. Some readers may wish to refer to the command's definition in gmic_stdlib.gmic.

Original picture. We mark four reference pixels. From top to bottom:

  1. (155, 51) Touch point of a small, upward pointing green arrow and the largest arrow

  2. (97, 80) A high contrast transition point along the upper vertical edge of the largest arrow

  3. (40,176) A low contrast transition point along lower edge of the leftmost green arrow

  4. (120,230) A point in a nearly constant field, nowhere near any transition

The same picture, now corrupted with noise at a variety of scales. We apply the -plasma command with iterations through two scales, dependent perturbations (alpha) of seventeen and scale independent perturbations (beta) of fourteen. We overlay a Gaussian noise with a variance of thirty.

Our trial -diffusiontensor parameters are:

  • sharpness: 0.7
  • anisotropy: 0.7
  • alpha: 4.0
  • sigma: 1.0
  • is_sqrt: 1

-diffusiontensors calls -structuretensors to locate the strength and orientation of image features.

  1. -diffusiontensors first applies isotropic blurring with a Gaussian variance of 4 (alpha). To a degree, this filtering affects ambient noise more than it does image features, but a too large variance will erode image features as well.

  2. Following alpha filtering, -diffusiontensors invokes -structuretensors. That command generates a tensor field, the coefficients of the individual tensors are computed from differences between pixels and their immediate left and right, top and bottom neighbors.

  3. The data comprising the tensor field produced by -structuretensors is not itself free from noise; some values will be outliers. Consequently, -diffusiontensors blurs the field produced by -structuretensors, employing a Gaussian variance of 1.0 (sigma). In effect, this blurring averages such outliers with the prevailing values of the neighborhood.

    As with alpha, one can go too far with sigma, as too large a variance can make indistinct smaller features that -structuretensors has detected.

Here are the tensors which -structuretensors calculated for our reference points, reflecting a local averaging from the second blurring operation (sigma):

  1. G [ 155,51 ] = [ 75.02135 62.97602 62.97602 65.44759 ] size 9 { bold {G}_[155,51] = left [ matrix{ 75.02135 # 62.97602 ## 62.97602 # 65.44759 } right ] }
  2. G [ 97,80 ] = [ 381.901 4.90887 4.90887 1.75710 ] size 9 { bold {G}_[97,80] = left [ matrix{ 381.901 # 4.90887 ## 4.90887 # 1.75710 } right ] }
  3. G [ 40,176 ] = [ 6.70915 − 3.57508 − 3.57508 55.02620 ] size 9 { bold {G}_[40,176] = left [ matrix{ 6.70915 # -3.57508 ## -3.57508 # 55.02620 } right ] }
  4. G [ 120,230 ] = [ 6.20982 − 1.00163 − 1.00163 6.77233 ] size 9 { bold {G}_[120,230] = left [ matrix{ 6.20982 # -1.00163 ## -1.00163 # 6.77233 } right ] }

While a compact representation, tensors do not lend themselves to easy visualization. Here, we solve for their eigenvalues and eigenvectors and overlay the data onto their structure tensor field locales. In each case, the gradient vector is the larger of the two and is directed along the most rapid change in intensity; the companion contour vector is barely visible at this scale. As noted in the text, the magnitude of the difference between λ+ and λ\u2500 signal edges, with the contour vector running along it and the gradient vector crossing it.

    1. (155, 51): λ– = 7.0768, θ- = 47.1734°; λ+ = 133.3922, θ+ = – 42.8266°

      The gradient vector (λ+  ∡ θ+) is about at right angles to the edge of the large arrow, oriented at an incline of about 42° in the original image. The contour vector is about five degrees more vertical than the edge; such variances are not surprising in vicinity of the touch point, where noise and blurring have altered the original contours.

    2. (97, 80): λ–. = 1.6937, θ– = 89.2603°; λ+ = 381.9644, θ+ = – 0.7397°

      In spite of noise and blurring, | λ+ – λ– | indicates a much stronger edge in this locale; The edge of the arrow is still clearly visible through the noise and the strength of the edge is reflected by the magnitude of | λ+ \u2500 λ\u2500 |. The gradient vector is only a fraction of a degree off from the perpendicular of the arrow edge, exactly 0° in the original.

    3. (40,176): λ–. = 6.4406, θ– = -4.2089°; λ+ = 55.2893, θ+ = 85.7911°

      Compared to the previous sample, the contrast between the leftmost olive green arrow and the cyan background is much lower, and the smaller magnitudes of the eigenvalues reflect this. Though small, the gradient vector is only a fraction of a degree off from the edge perpendicular.

    4. (120,230) λ–. = 5.4507 θ– = -37.1577 °; λ+ = 7.5314, θ+ = 52.8423°

      This sample, far from any features in the original, exhibits only very small gradient and contour vectors, just two percent of the largest vector at pixel (97, 80). While the data indicates an edge oriented around 142°, the small magnitude of | λ+ \u2500 λ\u2500 | suggests this is a spurious reading.

Given structure tensors and command parameters, we compute the smoothing magnitudes for the four points, harnessing the relationships that are given in the previous section. These magnitudes scale unit vectors with the same orientation as the input structure tensors and form smoothing vectors.

The diagrams on the left depict the smoothing vectors at each of the reference points and are drawn at a larger scale than the previous diagram, approximately at the ratio of 2,421:1.

These vectors instruct -smooth to diffuse the associated image pixel in an anisotropic way. In the first case, at pixel (155,51), diffusion will proceed at the relative magnitude of 0.03123 along an orientation of 47.1734°. The amplitude parameter given to -smooth scales these magnitudes to their working values. The overall tensor field embodies similar instructions for every pixel in the corrupted image, and, in aggregate constitute a smoothing geometry.

The smoothing vectors for the four reference pixels are:

  1. (155,51): s–. = 0.03123, θ– = 47.1734°, s+ = 0.000009589, θ+ = – 42.8266°

  2. (97,80): s–. = 0.0155, θ– = 89.2603°, s+ = 0.00006394, θ+ = – 0.7397°

  3. (40,176): s–. = 0.05517, θ– = -4.2089°, s+ = 0.000009589, θ+ = 85.7911°

  4. (120, 230): s–. = 0.1578, θ– = 52.8423°, s+ = 0.002123, θ+ = 142.8423°

-diffusiontensors calls -eigen2tensors. This command recomposes tensor fields from image pairs. The pair consists of a two channel image of eigenvalues and a two channel image of unit vectors; it produces a single tensor field encoding a smoothing geometry. This is passed as a parameter to -smooth.

In the present case, the two channel eigenvalue image contains the smoothing values; the unit vectors have been carried forward unchanged from the tensor field produced by -structuretensors. This operation gives rise to diffusion tensors for each of the reference points:

  1. T [ 155,51 ] = [ 0.009152 − 0.01521 − 0.01521 0.0253360 ] size 9 { {bold T}_[ 155,51 ] = left [ matrix {0.009152 # -0.01521 ## -0.01521 # 0.0253360 } right ] }
  2. T [ 97,80 ] = [ 0.00001078 − 0.0003918 − 0.0003918 0.0156 ] size 9 { {bold T}_[ 97,80 ] = left [ matrix {0.00001078 # -0.0003918 ## -0.0003918 # 0.0156} right ] }
  3. T [ 40,176 ] = [ 0.0556 0.005292 0.005292 0.0005718 ] size 9 { {bold T}_[ 40,176 ] = left [ matrix { 0.0556 # 0.005292 ## 0.005292 # 0.0005718 } right ] }
  4. T [ 120,230 ] = [ 0.07325 0.07467 0.07467 0.07992 ] size 9 { {bold T}_[ 120,230 ] = left [ matrix {0.07325 # 0.07467 ## 0.07467 # 0.07992} right ] }

The -smooth command blurs images through diffusion under the control of a given tensor field. In a sense, this command places over each pixel in the image an elliptic diffusion kernel, taken from the tensor field and which has been tailored to the pixel's immediate environment by -diffusiontensors. This gives rise to three broad behaviors:

  1. Should there be a strong gradient in the environment, the kernel is very eccentric but also not particularly large. It is oriented along the contour lines of the edge. Over the course of the command's run, the color value of the pixel will “bleed out” into neighboring pixels, the diffusion rate and pattern set by the values of s+ and s‒. Since, due to the strength of the local gradient, these values are not large, the diffusion will be very small. In short, the command adapts to strong edges and leaves them more or less alone.
  2. If the gradient is distinct enough to be distinguished from noise, but only just barely, the kernel is still markedly eccentric and oriented along the contour lines of the edge, but the values of s+ and s‒ are much larger. Over the course of the command's run, the color value of the central pixel will bleed out at faster rates than the previous case, and anisotropically, diffusing most strongly in the direction of θ‒ more or less collinear to the tightest contour lines in the neighborhood of the pixel. This tends to smear noise along, but not across, edges. In this case, the command adapts to weak edges by selectively smoothing along the edge.
  3. If there are no strong gradients in the neighborhood – a low contrast image region – the values λ+, λ─ tend to zero, the denominators of the contour and smoothing gradients tend to one, as do s+ and s‒. This gives rise to less elliptic diffusion kernels. Without a strong directional bias, these kind of kernels behave more like conventional, symmetric blurring kernels. Over the course of the command's run, the color value of the central pixel bleeds more or less uniformly in all directions. In the absence of edges, this is just the kind of smoothing needed; areas with uniform color and luminosity blur in an omnidirectional way, spreading irregularities in all directions.

In practice, we place the -smooth command in a -repeat 5 …-done loop, so that we may apply the diffusion map to a number of intermediary stages. For each stage we use:

  1. amplitude: 800
  2. spatial discretization (dl): 0.8 (default)
  3. angular discretization (da): 30° (default)
  4. diffusion precision: 2 (default)
  5. interpolation: nearest neighbor (default)
  6. Gaussian fast appoximation: enabled (default)

Neither -diffusiontensors nor -smooth can conjure image information wholly lost to noise. We harnessed -plasma to corrupt the image, which blurs as well as perturbs, so, from a spectral point of view, part of the image was destroyed by a low pass filter. Gone irretrievably are sharp corners and edges.

Within those limitations, -diffusiontensors recovered almost all information still manifest in the image. We could filter out most of the salt-and-pepper racket and recover areas of constant tone. Our reconstructed image looks like a somewhat blurred version of the original, the variance, one speculates, not too different from the blurring component in -plasma.

Comparing a magnified bit of the output with the pristine original illustrates some of the inner workings of -diffusiontensor and -smooth. The edges of the smoothed image seem to have shadowy, thin lines, quite unlike the softer regions in the interiors of shapes. These are footprints of very elliptical smoothing kernels, diffusing their pixels along thin, edge hugging paths. If the image hadn't been so terribly shot up with noise, these artifacts would be not nearly so manifest. Had we pushed the anisotropy parameter near its limits of one, these artifacts would start to look a little bit like curling hair strands; often these look quite lovely.

Garry Osgood