Lib_ContinuumRemoval.R 9.15 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
# prepares data to run multithreaded continuum removal
#
14
# @param Spectral_Data initial data matrix (nb samples x nb bands)
15 16 17 18 19
# @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
#' @importFrom future.apply future_lapply
22
apply_continuum_removal <- function(Spectral_Data, Spectral, nbCPU = 1) {
23
  if (!length(Spectral$WaterVapor) == 0) {
24
    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
  nb.Values <- dim(Spectral_Data)[1] * dim(Spectral_Data)[2]
29
  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
    Max.nb.Values <- 2e6
    nb.CR <- ceiling(nb.Values / Max.nb.Values)
34
    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
    Schedule_Per_Thread <- ceiling(nb.CR / nbCPU)
38
    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 41
    Spectral_Data <- do.call("rbind", Spectral_Data_tmp)
    rm(Spectral_Data_tmp)
jbferet's avatar
jbferet committed
42 43 44
  } 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
  }
  gc()
48
  return(Spectral_Data)
jbferet's avatar
jbferet committed
49 50
}

51 52 53 54 55 56
# 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)
57
# @param Spectral_Bands information about spectral bands
58 59
#
# @return samples from image and updated number of pixels to sampel if necessary
60
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
  nbBands <- dim(Minit)[2]
  CR_data$Minit <- c()
  Spectral_Bands <- CR_data$Spectral_Bands
  nbSamples <- CR_data$nbSamples
69
  nbSamplesUpDate <- length(CR_data$SamplesToKeep)
jbferet's avatar
jbferet committed
70
  # if samples to be considered
71
  if (nbSamples > 0) {
jbferet's avatar
jbferet committed
72 73
    # initialization:
    # spectral band corresponding to each element of the data matrix
74
    Lambda <- repmat(matrix(Spectral_Bands, nrow = 1), nbSamplesUpDate, 1)
jbferet's avatar
jbferet committed
75 76
    # prepare matrices used to check evolution of the CR process
    # - elements still not processed through continuum removal: initialization to 1
77
    Still.Need.CR <- matrix(1, nrow = nbSamplesUpDate, ncol = nbBands)
jbferet's avatar
jbferet committed
78
    # - value of the convex hull: initially set to 0
79
    Convex_Hull <- matrix(0, nrow = nbSamplesUpDate, ncol = nbBands)
jbferet's avatar
jbferet committed
80 81
    # - reflectance value for latest interception with convex hull:
    # initialization to value of the first reflectance measurement
82
    Intercept_Hull <- repmat(matrix(Minit[, 1], ncol = 1), 1, nbBands)
jbferet's avatar
jbferet committed
83
    # - spectral band of latest interception
84
    Latest.Intercept <- repmat(matrix(Spectral_Bands[1], ncol = 1), nbSamplesUpDate, nbBands)
jbferet's avatar
jbferet committed
85
    # number of spectral bands found as intercept
86
    nb.Intercept <- 0
jbferet's avatar
jbferet committed
87 88
    # continues until arbitrary stopping criterion:
    # stops when reach last spectral band (all values before last = 0)
89 90
    # while (max(Still.Need.CR[, 1:(nbBands - 2)]) == 1 & (nb.Intercept <= (nbBands / 2))) {
    while (max(Still.Need.CR[, 1:(nbBands - 2)]) == 1) {
91
      nb.Intercept <- nb.Intercept + 1
92 93 94 95 96 97 98 99 100 101
      # identify samples still needing continuum removal
      Sel <- which(Still.Need.CR[,(nbBands-2)]==1)
      # update variables to process samples needing CR only
      nbSamplesUpDate_tmp <- length(Sel)
      Lambda_tmp <- matrix(Lambda[Sel,],nrow = nbSamplesUpDate_tmp)
      Minit_tmp <- matrix(Minit[Sel,],nrow = nbSamplesUpDate_tmp)
      Latest.Intercept_tmp <- matrix(Latest.Intercept[Sel,],nrow = nbSamplesUpDate_tmp)
      Still.Need.CR_tmp <-matrix(Still.Need.CR[Sel,],nrow = nbSamplesUpDate_tmp)
      Convex_Hull_tmp <-matrix(Convex_Hull[Sel,],nrow = nbSamplesUpDate_tmp)
      Intercept_Hull_tmp <-matrix(Intercept_Hull[Sel,],nrow = nbSamplesUpDate_tmp)
jbferet's avatar
jbferet committed
102
      # Mstep give the position of the values to be updated
103
      Update_Data <- matrix(1, nrow = nbSamplesUpDate_tmp, ncol = nbBands)
104
      Update_Data[, nbBands] <- 0
jbferet's avatar
jbferet committed
105 106
      # initial step: first column set to 0; following steps: all bands below
      # max of the convex hull are set to 0
107
      Update_Data[which((Lambda_tmp - Latest.Intercept_tmp) < 0)] <- 0
jbferet's avatar
jbferet committed
108
      # compute slope for each coordinate
109
      Slope <- (Minit_tmp - Intercept_Hull_tmp) / (Lambda_tmp - Latest.Intercept_tmp) * Still.Need.CR_tmp
jbferet's avatar
jbferet committed
110
      # set current spectral band and previous bands to -9999
111 112
      if (!length(which(Still.Need.CR_tmp == 0)) == 0) {
        Slope[which(Still.Need.CR_tmp == 0)] <- -9999
jbferet's avatar
jbferet committed
113
      }
114 115
      if (!length(which(is.na(Slope))) == 0) {
        Slope[which(is.na(Slope))] <- -9999
jbferet's avatar
jbferet committed
116 117
      }
      # get max index for each row and convert into linear index
118
      Index.Max.Slope <- RowToLinear(max.col(Slope, ties.method = "last"), nbSamplesUpDate_tmp, nbBands)
jbferet's avatar
jbferet committed
119 120
      # !!!! OPTIM: replace repmat with column operation
      # update coordinates of latest intercept
121
      Latest.Intercept_tmp <- repmat(matrix(Lambda_tmp[Index.Max.Slope], ncol = 1), 1, nbBands)
jbferet's avatar
jbferet committed
122
      # update latest intercept
123
      Intercept_Hull_tmp <- repmat(matrix(Minit_tmp[Index.Max.Slope], ncol = 1), 1, nbBands)
jbferet's avatar
jbferet committed
124
      # values corresponding to the domain between the two continuum maxima
125
      Update_Data[which((Lambda_tmp - Latest.Intercept_tmp) >= 0 | Latest.Intercept_tmp == Spectral_Bands[nbBands])] <- 0
jbferet's avatar
jbferet committed
126
      # values to eliminate for the next analysis: all spectral bands before latest intercept
127
      Still.Need.CR_tmp[which((Lambda_tmp - Latest.Intercept_tmp) < 0)] <- 0
jbferet's avatar
jbferet committed
128 129
      # the max slope is known, as well as the coordinates of the beginning and ending
      # a matrix now has to be built
130 131 132 133 134 135 136
      Convex_Hull_tmp <- Convex_Hull_tmp + Update_Data * (Intercept_Hull_tmp + sweep((Lambda_tmp - Latest.Intercept_tmp), 1, Slope[Index.Max.Slope], "*"))
      # update variables
      Convex_Hull[Sel,] <- Convex_Hull_tmp
      Still.Need.CR[Sel,] <- Still.Need.CR_tmp
      Lambda[Sel,] <- Lambda_tmp
      Latest.Intercept[Sel,] <- Latest.Intercept_tmp
      Intercept_Hull[Sel,] <- Intercept_Hull_tmp
jbferet's avatar
jbferet committed
137
    }
138 139 140
    CR_Results0 <- Minit[, 2:(nbBands - 2)] / Convex_Hull[, 2:(nbBands - 2)]
    CR_Results <- matrix(0, ncol = (nbBands - 3), nrow = nbSamples)
    CR_Results[CR_data$SamplesToKeep, ] <- CR_Results0
jbferet's avatar
jbferet committed
141
  } else {
142
    CR_Results <- matrix(0, ncol = (nbBands - 3), nrow = nbSamples)
jbferet's avatar
jbferet committed
143
  }
144 145
  list <- ls()
  rm(list = list[-which(list == "CR_Results")])
jbferet's avatar
jbferet committed
146 147 148 149 150
  rm(list)
  gc()
  return(CR_Results)
}

151 152 153 154 155 156
# 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
#
157
# @param Spectral_Bands
158 159 160 161
# @param Minit initial data matrix, n rows = n samples, p cols = p spectral bands
#
# @return updated Minit
#' @importFrom matrixStats rowSds
162
filter_prior_CR <- function(Minit, Spectral_Bands) {
jbferet's avatar
jbferet committed
163 164

  # number of samples to be processed
165
  nbSamples <- nrow(Minit)
jbferet's avatar
jbferet committed
166
  # make sure there is no negative values
167 168
  Minit <- Minit + 100.0
  Minit[which(Minit < 0)] <- 0
jbferet's avatar
jbferet committed
169
  # eliminate invariant spectra
170 171 172 173 174 175 176 177 178
  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
179
  }
180 181 182
  Minit <- Minit[keep, ]
  if (length(keep) == 1) {
    Minit <- matrix(Minit, nrow = 1)
jbferet's avatar
jbferet committed
183
  }
184
  nbSamplesUpDate <- nrow(Minit)
jbferet's avatar
jbferet committed
185
  # add negative values to the last column and update spectral bands
186
  Minit <- cbind(Minit, matrix(-9999, ncol = 1, nrow = nbSamplesUpDate))
187 188 189
  nbBands <- ncol(Minit)
  Spectral_Bands[nbBands] <- Spectral_Bands[nbBands - 1] + 10
  my_list <- list("Minit" = Minit, "Spectral_Bands" = Spectral_Bands, "nbSamples" = nbSamples, "SamplesToKeep" = keep)
jbferet's avatar
jbferet committed
190 191 192
  return(my_list)
}

193 194 195 196 197 198 199 200 201
#
# @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
202 203 204
  return(MaxCont)
}

205 206 207 208 209 210 211 212 213 214 215
# 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
216
}