Chapter 5: Image processing

Here you can find the code listings in R language from the corresponding chapter of the book.

You need to load Rssa, ssabook, lattice, latticeExtra, plyr, fma to run these examples.

Contents

  1. Chapter 5: Image processing
    1. Fragments 5.1.1 (‘Mars’: Input), 5.1.2 (‘Mars’: Decomposition with svd.method = "svd"), 5.1.3 (‘Mars’: Decomposition with svd.method = "eigen")
      1. Produced output
    2. Fragments 5.1.4 (‘Mars’: Decomposition) and 5.1.5 (‘Mars’: Reconstruction)
      1. Produced output
    3. Fragment 5.1.6 (‘Mars’: Decomposition)
      1. Produced output
    4. Fragment 5.1.7 (‘Mars’: Reconstruction)
      1. Produced output
    5. Fragments 5.2.1 (Auxiliary plot of 2D image) and 5.2.2 (‘Mars’: Mask specification and decomposition)
      1. Produced output
    6. Fragment 5.2.3 (‘Mars’: Reconstruction)
      1. Produced output
    7. Fragments 5.2.4 (‘Mars’: Identification) and 5.2.5 (‘Mars’: Improvement by DerivSSA)
      1. Produced output
    8. Fragment 5.3.1 (‘Mars’: Parameter estimation with 2D-ESPRIT)
      1. Produced output
    9. Fragment 5.3.2 (‘Mars’: Parameter estimation with Shaped 2D-ESPRIT)
      1. Produced output
    10. Fragment 5.4.1 (Mars: magnified reconstructions by 2D-SSA and ShSSA)
      1. Produced output
    11. Fragment 5.4.2 (‘Brecon Beacons’: Decomposition)
      1. Produced output
    12. Fragment 5.4.3 (‘Brecon Beacons’: Reconstruction)
      1. Produced output
    13. Fragments 5.4.4 (2D-SSA: centered DFT) and 5.4.5 (2D-SSA: DFT of cumulative reconstructions)
      1. Produced output
    14. Fragments 5.4.6 (‘Kruppel’: Analysis of data given on a cylinder)
      1. Produced output
    15. Fragments 5.4.7 (‘Monet’: decomposition by multivariate 2D-SSA)
      1. Produced output

Fragments 5.1.1 (‘Mars’: Input), 5.1.2 (‘Mars’: Decomposition with svd.method = "svd"), 5.1.3 (‘Mars’: Decomposition with svd.method = "eigen")

data("Mars", package = "Rssa")
# ssa(Mars, kind = "2d-ssa", L = c(25, 25), svd.method = "svd")
print(system.time(ssa(Mars, kind = "2d-ssa", L = c(25, 25),
                      svd.method = "eigen")))

Produced output

   user  system elapsed 
  5.706   0.663   6.877

Fragments 5.1.4 (‘Mars’: Decomposition) and 5.1.5 (‘Mars’: Reconstruction)

print(system.time(s.Mars.25 <- ssa(Mars, kind = "2d-ssa", 
                                   L = c(25, 25))))
r.Mars.25 <-
    reconstruct(s.Mars.25,
                groups = list(Noise = c(12, 13, 15, 16)))
plot(r.Mars.25, cuts = 255, layout = c(3, 1))

Produced output

Reconstruction ‘Mars’: Separated periodic noise, (Lx, Ly)=(25, 25).

   user  system elapsed 
  0.736   0.014   0.866

Fragment 5.1.6 (‘Mars’: Decomposition)

plot(s.Mars.25, type = "vectors", idx = 1:20,
     cuts = 255, layout = c(10, 2),
     plot.contrib = FALSE)
plot(wcor(s.Mars.25, groups = 1:30),
     scales = list(at = c(10, 20, 30)))

Produced output

Eigenarrays ‘Mars’: Eigenarrays, (Lx, Ly)=(25, 25).

Wcor ‘Mars’: w-Correlations (Lx, Ly)=(25, 25).

Fragment 5.1.7 (‘Mars’: Reconstruction)

print(system.time(s.Mars.160.80 <-
  ssa(Mars, kind = "2d-ssa", L = c(160, 80))))
r.Mars.160.80.groups <- list(Noise = c(36, 37, 42, 43))
r.Mars.160.80 <- reconstruct(s.Mars.160.80, 
                             groups = r.Mars.160.80.groups)
plot(r.Mars.160.80, cuts = 255, layout = c(3, 1))

Produced output

Reconstruction ‘Mars’: Reconstruction (Lx, Ly)=(160, 80).

   user  system elapsed 
  0.828   0.032   0.959

Fragments 5.2.1 (Auxiliary plot of 2D image) and 5.2.2 (‘Mars’: Mask specification and decomposition)

plot2d <- function(x) {
  regions <- list(col = colorRampPalette(grey(c(0, 1))));
  levelplot(t(x[seq(nrow(x), 1, -1), ]), aspect = "iso",
            par.settings = list(regions = regions), 
            colorkey = FALSE, 
            scales = list(draw = FALSE, relation = "same"), 
            xlab = "", ylab = "")
}
mask.Mars.0 <- (Mars != 0)
mask.Mars.1 <- (Mars != 255)
Mars[!mask.Mars.0] <- NA
print(system.time(s.Mars.shaped <-
  ssa(Mars, kind = "2d-ssa",
      mask = mask.Mars.1, wmask = circle(15))))
mask.Mars.res <- (s.Mars.shaped$weights > 0)
plot2d(mask.Mars.0)
plot2d(mask.Mars.1)
plot2d(mask.Mars.res)

Produced output

Mask Mask Mask
'Mars' masks specification. Left: specified by NA, center: the parameter mask, right: resulting mask. White and black colors correspond to TRUE and FALSE respectively.
   user  system elapsed 
  0.946   0.026   1.171

Fragment 5.2.3 (‘Mars’: Reconstruction)

r.Mars.shaped.groups <- list(Noise = c(7, 8, 9, 10))
r.Mars.shaped <- reconstruct(s.Mars.shaped,
                    groups = r.Mars.shaped.groups)
plot(r.Mars.shaped, cuts = 255, layout = c(3, 1), 
     fill.color = "yellow")

Produced output

Reconstruction ‘Mars’: Reconstruction, ShSSA.

Fragments 5.2.4 (‘Mars’: Identification) and 5.2.5 (‘Mars’: Improvement by DerivSSA)

plot(s.Mars.shaped, type = "vectors", idx = 1:30, 
     fill.color = "yellow", cuts = 255, layout = c(10, 3), 
     plot.contrib = FALSE)
plot(wcor(s.Mars.shaped, groups = 1:30),
     scales = list(at = c(10, 20, 30)))
plot(s.Mars.shaped)
s.Mars.shaped.deriv <-
  fossa(s.Mars.shaped, nested.groups = list(1:30))
plot(s.Mars.shaped.deriv, type = "vectors", idx = 1:30, 
     fill.color = "yellow", cuts = 255, layout = c(10, 3), 
     plot.contrib = FALSE)
plot(wcor(s.Mars.shaped.deriv, groups = 1:30), 
     scales = list(at = c(10, 20, 30)))
plot(s.Mars.shaped.deriv)

Produced output

Eigenarrays ‘Mars’: Eigenarrays, ShSSA.

Eigenarrays ‘Mars’: Eigenarrays, ShSSA; improvement by DerivSSA.

Wcor Contributions
'Mars': w-Correlations, ShSSA: initial (left) and after DerivSSA.
Wcor Contributions
'Mars': Contribution of elementary components, ShSSA: initial (left) and after DerivSSA.

Fragment 5.3.1 (‘Mars’: Parameter estimation with 2D-ESPRIT)

pe.Mars.160.80 <- parestimate(s.Mars.160.80, 
                              groups = r.Mars.160.80.groups)
print(pe.Mars.160.80)
print(pe.Mars.160.80[[1]])
print(pe.Mars.160.80[[2]])
plot(pe.Mars.160.80, col = c(11, 12, 13, 14))
plot(s.Mars.160.80, type = "vectors", 
     idx = r.Mars.160.80.groups$Noise,
     cuts = 255, layout = c(4, 1), plot.contrib = FALSE)

Produced output

Roots ‘Mars’: Parameter estimation with 2D-ESPRIT, (Lx, Ly)=(160, 80).

Eigenarrays ‘Mars’: Eigenarrays corresponding to the periodic noise, (Lx, Ly)=(160, 80).

       x: period       rate |        y: period       rate
          -5.000  -0.000169 |           10.003  -0.000111
           5.000  -0.000169 |          -10.003  -0.000111
           9.995   0.000175 |            4.999  -0.000093
          -9.995   0.000175 |           -4.999  -0.000093
   period     rate   |    Mod     Arg  |     Re        Im
   -5.000  -0.000169 |  0.99983  -1.26 |  0.30906  -0.95087
    5.000  -0.000169 |  0.99983   1.26 |  0.30906   0.95087
    9.995   0.000175 |  1.00017   0.63 |  0.80897   0.58814
   -9.995   0.000175 |  1.00017  -0.63 |  0.80897  -0.58814
   period     rate   |    Mod     Arg  |     Re        Im
   10.003  -0.000111 |  0.99989   0.63 |  0.80905   0.58755
  -10.003  -0.000111 |  0.99989  -0.63 |  0.80905  -0.58755
    4.999  -0.000093 |  0.99991   1.26 |  0.30879   0.95103
   -4.999  -0.000093 |  0.99991  -1.26 |  0.30879  -0.95103

Fragment 5.3.2 (‘Mars’: Parameter estimation with Shaped 2D-ESPRIT)

pe.Mars.shaped <- parestimate(s.Mars.shaped, 
                              r.Mars.shaped.groups)
print(pe.Mars.shaped)

Produced output

       x: period       rate |        y: period       rate
          -5.007  -0.001403 |           10.015  -0.000408
           5.007  -0.001403 |          -10.015  -0.000408
          10.008   0.000350 |            5.004  -0.001084
         -10.008   0.000350 |           -5.004  -0.001084

Fragment 5.4.1 (Mars: magnified reconstructions by 2D-SSA and ShSSA)

Mars.sh <- r.Mars.shaped$Noise
Mars.rect.sh <- Mars.rect <- r.Mars.25$Noise
Mars.rect.sh[is.na(Mars.sh)] <- NA
library("latticeExtra")
p.part.rect <- plot2d(Mars.rect[60:110, 200:250]) +
                 layer(panel.fill(col = "yellow", 
                                  alpha = 0.2), 
                       under = FALSE) +
                 plot2d(Mars.rect.sh[60:110, 200:250])
p.part.shaped <- plot2d(r.Mars.shaped[[1]][60:110, 200:250]) +
                   layer(panel.fill(col = "yellow"), 
                         under = TRUE)
plot(c(Rectangular = p.part.rect, Shaped = p.part.shaped))

Produced output

Comparison ‘Mars’: comparison of texture reconstructions by 2D-SSA and ShSSA.

Fragment 5.4.2 (‘Brecon Beacons’: Decomposition)

data("brecon", package = "ssabook")
s.brecon <- ssa(brecon, kind = "2d-ssa", L = c(8, 8),
                svd.method = "eigen")
plot(s.brecon, type = "vectors", idx = 1:32,
     cuts = 255, layout = c(8, 4), plot.contrib = FALSE)

Produced output

Eigenarrays ‘Brecon Beacons’: 8 x 8 windows, eigenarrays.

Fragment 5.4.3 (‘Brecon Beacons’: Reconstruction)

r.brecon <- reconstruct(s.brecon, 
                        groups = list(1:3, 4:8, 9:17))
plot(r.brecon, cuts = 255, layout = c(5, 1), 
     par.strip.text = list(cex = 0.75), 
     col = topo.colors(1000))
plot(r.brecon, cuts = 255, layout = c(5, 1), 
     par.strip.text = list(cex = 0.75), 
     type = "cumsum", at = "free", 
     col = topo.colors(1000))

Produced output

Reconstruction ‘Brecon Beacons’: 8 x 8 window, reconstructions (𝕏̃k).

Reconstruction ‘Brecon Beacons’: 8 x 8 window, cumulative reconstructions (𝕐̃k).

Fragments 5.4.4 (2D-SSA: centered DFT) and 5.4.5 (2D-SSA: DFT of cumulative reconstructions)

centered.mod.fft <- function(x) {
  N <- dim(x)
  shift.exp <- exp(2i*pi * floor(N/2) / N)
  shift1 <- shift.exp[1]^(0:(N[1] - 1))
  shift2 <- shift.exp[2]^(0:(N[2] - 1))
  Mod(t(mvfft(t(mvfft(outer(shift1, shift2) * x)))))
}
plot2d(centered.mod.fft(brecon - r.brecon$F1))
plot2d(centered.mod.fft(brecon - r.brecon$F1 - r.brecon$F2))
plot2d(centered.mod.fft(brecon - r.brecon$F1 - 
                          r.brecon$F2 - r.brecon$F3))

Produced output

Reconstruction Reconstruction Reconstruction
'Brecon Beacons': 8 x 8 window, absolute values of the DFT of 𝕏 - 𝕐̃k, k=1,2,3.

Fragments 5.4.6 (‘Kruppel’: Analysis of data given on a cylinder)

data("kruppel", package = "ssabook")
mask <- matrix(0, ncol = 200, nrow = 200)
mask[, 40:160] <- 1
s <- ssa(kruppel, L = c(20, 20), mask = mask,  
         circular = c(TRUE, FALSE))
plot(s, type = "vectors", vectors = "factor", 
     idx = 1:12, aspect = 1.4, col = topo.colors(1000))
rec6 <- reconstruct(s, list("C1-6" = 1:6))
p <- plot(rec6, aspect = 1.4, col = topo.colors(1000))
p$condlevels[[1]] <- 
  sub("Residuals", "Resd.", p$condlevels[[1]], fixed = TRUE)
plot(p)

Produced output

Factor vectors ‘Kruppel’: Factor vectors.

Reconstruction ‘Kruppel’: Original image, reconstruction and residual.

Fragments 5.4.7 (‘Monet’: decomposition by multivariate 2D-SSA)

data("monet", package = "ssabook")
plot(c(0, 1), c(0, 1), type = "n", xlab = "", ylab = "", 
     axes = FALSE)
rasterImage(monet, 0, 0, 1, 1, interpolate = FALSE)
ss <- ssa(monet, L = c(30, 30, 1))
plot(ss, type = "vectors", idx = 1:20, slice = list(k = 1), 
     plot.contrib = FALSE)
rec <- reconstruct(ss, groups = list(smooth = 1:6))
rec$smooth[rec$smooth > 1] <- 1
rec$smooth[rec$smooth < 0] <- 0
plot(c(0, 1), c(0, 1), type = "n", xlab = "", ylab = "", 
     axes = FALSE)
rasterImage(rec$smooth, 0, 0, 1, 1, 
            interpolate = FALSE)
p1 <- plot(rec, slice = list(k = 1), main = "Red", 
           col = c("#000000", "#FF7070"), layout = c(3, 1))
p2 <- plot(rec, slice = list(k = 2), main = "Green", 
           col = c("#000000", "#70FF70"), layout = c(3, 1))
p3 <- plot(rec, slice = list(k = 3), main = "Blue",
           col = c("#000000", "#7070FF"), layout = c(3, 1))
plot(p1, split = c(1, 1, 1, 3), more = TRUE)
plot(p2, split = c(1, 2, 1, 3), more = TRUE)
plot(p3, split = c(1, 3, 1, 3), more = FALSE)

Produced output

Original Smooth
'Monet': Initial (left) and smoothed (right) images.

Eigenarrays ‘Monet’: Eigenarrays.

Channels ‘Monet’: Three channels of the reconstructed and residual images.