Lib_ContinuumRemoval.R 8.01 KB
Newer Older
jbferet's avatar
jbferet committed
1
# ===============================================================================
Florian de Boissieu's avatar
Florian de Boissieu committed
2
# biodivMapR
jbferet's avatar
jbferet committed
3
4
5
6
7
8
9
10
11
# Lib_ContinuumRemoval.R
# ===============================================================================
# PROGRAMMERS:
# Jean-Baptiste FERET <jb.feret@irstea.fr>
# Copyright 2018/07 Jean-Baptiste FERET
# ===============================================================================
# This Library is dedicated to the computation of the continuum removal
# ===============================================================================

12
13
14
15
16
17
18
19
# prepares data to run multithreaded continuum removal
#
# @param Spectral.Data initial data matrix (nb samples x nb bands)
# @param nbCPU
# @param Spectral information about spectral bands
#
# @return samples from image and updated number of pixels to sampel if necessary
#' @importFrom snow splitRows
De Boissieu Florian's avatar
De Boissieu Florian committed
20
#' @importFrom future plan multiprocess sequential
21
22
23
24
#' @importFrom future.apply future_lapply
Apply.Continuum.Removal <- function(Spectral.Data, Spectral, nbCPU = 1) {
  if (!length(Spectral$WaterVapor) == 0) {
    Spectral.Data <- Spectral.Data[, -Spectral$WaterVapor]
jbferet's avatar
jbferet committed
25
26
27
  }

  # split data to perform continuum removal on into reasonable amount of data
28
29
  nb.Values <- dim(Spectral.Data)[1] * dim(Spectral.Data)[2]
  if (nb.Values > 0) {
jbferet's avatar
jbferet committed
30
31
    # corresponds to ~ 40 Mb data, but CR tends to requires ~ 10 times memory
    # avoids memory crash
32
33
34
    Max.nb.Values <- 2e6
    nb.CR <- ceiling(nb.Values / Max.nb.Values)
    Spectral.Data <- splitRows(Spectral.Data, nb.CR)
jbferet's avatar
jbferet committed
35
36
    # perform multithread continuum removal
    plan(multiprocess, workers = nbCPU) ## Parallelize using four cores
37
38
    Schedule.Per.Thread <- ceiling(nb.CR / nbCPU)
    Spectral.Data.tmp <- future_lapply(Spectral.Data, FUN = ContinuumRemoval, Spectral.Bands = Spectral$Wavelength, future.scheduling = Schedule.Per.Thread)
jbferet's avatar
jbferet committed
39
    plan(sequential)
40
    Spectral.Data <- do.call("rbind", Spectral.Data.tmp)
jbferet's avatar
jbferet committed
41
42
43
44
    rm(Spectral.Data.tmp)
  } else {
    # edit 31-jan-2018
    # resize to delete first and last band as in continuum removal
45
    Spectral.Data <- Spectral.Data[, -c(1, 2)]
jbferet's avatar
jbferet committed
46
47
48
49
50
  }
  gc()
  return(Spectral.Data)
}

51
52
53
54
55
56
57
58
59
60
# Computes continuum removal for matrix shaped data: more efficient than
# processing individual spectra
# the convex hull is based on the computation of the derivative between R at a
# given spectral band and R at the following bands
#
# @param Minit initial data matrix (nb samples x nb bands)
# @param Spectral.Bands information about spectral bands
#
# @return samples from image and updated number of pixels to sampel if necessary
ContinuumRemoval <- function(Minit, Spectral.Bands) {
jbferet's avatar
jbferet committed
61
62

  # Filter and prepare data prior to continuum removal
63
64
65
66
67
68
  CR.data <- Filter.Prior.CR(Minit, Spectral.Bands)
  Minit <- CR.data$Minit
  nb.Bands <- dim(Minit)[2]
  CR.data$Minit <- c()
  Spectral.Bands <- CR.data$Spectral.Bands
  nb.Samples <- CR.data$nb.Samples
jbferet's avatar
jbferet committed
69
  # if samples to be considered
70
  if (nb.Samples > 0) {
jbferet's avatar
jbferet committed
71
72
    # initialization:
    # spectral band corresponding to each element of the data matrix
73
    Lambda <- repmat(matrix(Spectral.Bands, nrow = 1), nb.Samples, 1)
jbferet's avatar
jbferet committed
74
75
    # prepare matrices used to check evolution of the CR process
    # - elements still not processed through continuum removal: initialization to 1
76
    Still.Need.CR <- matrix(1, nrow = nb.Samples, ncol = nb.Bands)
jbferet's avatar
jbferet committed
77
    # - value of the convex hull: initially set to 0
78
    Convex.Hull <- matrix(0, nrow = nb.Samples, ncol = nb.Bands)
jbferet's avatar
jbferet committed
79
80
    # - reflectance value for latest interception with convex hull:
    # initialization to value of the first reflectance measurement
81
    Intercept.Hull <- repmat(matrix(Minit[, 1], ncol = 1), 1, nb.Bands)
jbferet's avatar
jbferet committed
82
    # - spectral band of latest interception
83
    Latest.Intercept <- repmat(matrix(Spectral.Bands[1], ncol = 1), nb.Samples, nb.Bands)
jbferet's avatar
jbferet committed
84
    # number of spectral bands found as intercept
85
    nb.Intercept <- 0
jbferet's avatar
jbferet committed
86
87
    # continues until arbitrary stopping criterion:
    # stops when reach last spectral band (all values before last = 0)
88
89
    while (max(Still.Need.CR[, 1:(nb.Bands - 2)]) == 1 & (nb.Intercept <= (nb.Bands / 2))) {
      nb.Intercept <- nb.Intercept + 1
jbferet's avatar
jbferet committed
90
      # Mstep give the position of the values to be updated
91
92
      Update.Data <- matrix(1, nrow = nb.Samples, ncol = nb.Bands)
      Update.Data[, nb.Bands] <- 0
jbferet's avatar
jbferet committed
93
94
      # initial step: first column set to 0; following steps: all bands below
      # max of the convex hull are set to 0
95
      Update.Data[which((Lambda - Latest.Intercept) < 0)] <- 0
jbferet's avatar
jbferet committed
96
      # compute slope for each coordinate
97
      Slope <- (Minit - Intercept.Hull) / (Lambda - Latest.Intercept) * Still.Need.CR
jbferet's avatar
jbferet committed
98
      # set current spectral band and previous bands to -9999
99
100
      if (!length(which(Still.Need.CR == 0)) == 0) {
        Slope[which(Still.Need.CR == 0)] <- -9999
jbferet's avatar
jbferet committed
101
      }
102
103
      if (!length(which(is.na(Slope))) == 0) {
        Slope[which(is.na(Slope))] <- -9999
jbferet's avatar
jbferet committed
104
105
      }
      # get max index for each row and convert into linear index
106
      Index.Max.Slope <- RowToLinear(max.col(Slope, ties.method = "last"), nb.Samples, nb.Bands)
jbferet's avatar
jbferet committed
107
108
      # !!!! OPTIM: replace repmat with column operation
      # update coordinates of latest intercept
109
      Latest.Intercept <- repmat(matrix(Lambda[Index.Max.Slope], ncol = 1), 1, nb.Bands)
jbferet's avatar
jbferet committed
110
      # update latest intercept
111
      Intercept.Hull <- repmat(matrix(Minit[Index.Max.Slope], ncol = 1), 1, nb.Bands)
jbferet's avatar
jbferet committed
112
      # values corresponding to the domain between the two continuum maxima
113
      Update.Data[which((Lambda - Latest.Intercept) >= 0 | Latest.Intercept == Spectral.Bands[nb.Bands])] <- 0
jbferet's avatar
jbferet committed
114
      # values to eliminate for the next analysis: all spectral bands before latest intercept
115
      Still.Need.CR[which((Lambda - Latest.Intercept) < 0)] <- 0
jbferet's avatar
jbferet committed
116
117
      # the max slope is known, as well as the coordinates of the beginning and ending
      # a matrix now has to be built
118
      Convex.Hull <- Convex.Hull + Update.Data * (Intercept.Hull + sweep((Lambda - Latest.Intercept), 1, Slope[Index.Max.Slope], "*"))
jbferet's avatar
jbferet committed
119
    }
120
121
122
    CR_Results0 <- Minit[, 2:(nb.Bands - 2)] / Convex.Hull[, 2:(nb.Bands - 2)]
    CR_Results <- matrix(0, ncol = (nb.Bands - 3), nrow = nb.Samples)
    CR_Results[CR.data$Samples.To.Keep, ] <- CR_Results0
jbferet's avatar
jbferet committed
123
  } else {
124
    CR_Results <- matrix(0, ncol = (nb.Bands - 3), nrow = nb.Samples)
jbferet's avatar
jbferet committed
125
  }
126
127
  list <- ls()
  rm(list = list[-which(list == "CR_Results")])
jbferet's avatar
jbferet committed
128
129
130
131
132
  rm(list)
  gc()
  return(CR_Results)
}

133
134
135
136
137
138
139
140
141
142
143
144
# Filter data prior to continuum removal:
# - values are expected to be real reflectance values between 0 and 10000
# - negative values may occur, so a +100 value is applied to avoid negative
# - possibly remaining negative values are set to 0
# - constant spectra are eliminated
#
# @param Spectral.Bands
# @param Minit initial data matrix, n rows = n samples, p cols = p spectral bands
#
# @return updated Minit
#' @importFrom matrixStats rowSds
Filter.Prior.CR <- function(Minit, Spectral.Bands) {
jbferet's avatar
jbferet committed
145
146

  # number of samples to be processed
147
  nb.Samples <- nrow(Minit)
jbferet's avatar
jbferet committed
148
  # make sure there is no negative values
149
150
  Minit <- Minit + 100.0
  Minit[which(Minit < 0)] <- 0
jbferet's avatar
jbferet committed
151
  # eliminate invariant spectra
152
153
154
155
156
157
158
159
160
  SD <- rowSds(Minit)
  elim <- which(SD == 0 | is.na(SD))
  keep <- which(!SD == 0 & !is.na(SD))
  nbi0 <- nrow(Minit)
  nbj <- ncol(Minit)
  if (nbj == 1 & nbi0 > 1) {
    Minit <- matrix(Minit, nrow = 1)
    nbi0 <- 1
    nbj <- ncol(Minit)
jbferet's avatar
jbferet committed
161
  }
162
163
164
  Minit <- Minit[keep, ]
  if (length(keep) == 1) {
    Minit <- matrix(Minit, nrow = 1)
jbferet's avatar
jbferet committed
165
166
  }
  # add negative values to the last column and update spectral bands
167
168
169
170
  Minit <- cbind(Minit, matrix(-9999, ncol = 1, nrow = nb.Samples))
  nb.Bands <- ncol(Minit)
  Spectral.Bands[nb.Bands] <- Spectral.Bands[nb.Bands - 1] + 10
  my_list <- list("Minit" = Minit, "Spectral.Bands" = Spectral.Bands, "nb.Samples" = nb.Samples, "Samples.To.Keep" = keep)
jbferet's avatar
jbferet committed
171
172
173
  return(my_list)
}

174
175
176
177
178
179
180
181
182
#
# @param MM
# @param nbi
# @param nbj
#
# @return
RowToLinear <- function(MM, nbi, nbj) {
  adj <- seq(1:nbi)
  MaxCont <- ((MM - 1) * (nbi)) + adj
jbferet's avatar
jbferet committed
183
184
185
  return(MaxCont)
}

186
187
188
189
190
191
192
193
194
195
196
# R equivalent of repmat (matlab)
#
# @param X initial matrix
# @param m nb of replications in row dimension
# @param n nb of replications in column dimension
#
# @return
repmat <- function(X, m, n) {
  mx <- dim(X)[1]
  nx <- dim(X)[2]
  matrix(t(matrix(X, mx, nx * n)), mx * m, nx * n, byrow = T)
jbferet's avatar
jbferet committed
197
}