Commit 40acdaae authored by De Boissieu Florian's avatar De Boissieu Florian
Browse files

merge with master which had most recent modifications

Merge branch 'master' into CRAN

# Conflicts:
#	DESCRIPTION
#	R/Lib_FilterData.R
#	man/check_data.Rd
#	man/perform_radiometric_filtering.Rd
parents be326cc8 5550ea56
......@@ -3,5 +3,4 @@ doc
.Rproj.user
.Rbuildignore
.Rhistory
RESULTS
02_REDACTION
TODO.md
Package: biodivMapR
Title: biodivMapR: an R package for α- and β-diversity mapping using remotely-sensed images
Version: 0.9.0
Authors@R: c( person("Jean-Baptiste", "Feret", email = "jb.feret@teledetection.fr", role = c("aut", "cre")),
person("Florian", "de Boissieu", email = "florian.deboissieu@irstea.fr", role = c("ctb"), comment = "clean code, format as package"))
Authors@R: c(person(given = "Jean-Baptiste",
family = "Feret",
email = "jb.feret@teledetection.fr",
role = c("aut", "cre"),
comment = c(ORCID = "0000-0002-0151-1334")),
person( given = "Florian",
family = "de Boissieu",
email = "florian.deboissieu@irstea.fr",
role = c("ctb"),
comment = "clean code, format as package"))
Description: Biodiversity estimation from multispectral and hyperspectral remote sensing images.
The methods are based on "Féret, J.-B., Asner, G.P., 2014.
Mapping tropical forest canopy diversity using high-fidelity imaging spectroscopy.
......@@ -21,10 +29,12 @@ Imports:
labdsv,
matlab,
matrixStats,
methods,
raster,
rgdal,
R.utils,
snow,
sp,
stringr,
tools,
vegan,
......
# Generated by roxygen2: do not edit by hand
S3method(split,line)
export(check_data)
export(diversity_from_plots)
export(get_projection)
......@@ -12,17 +11,21 @@ export(perform_PCA)
export(perform_radiometric_filtering)
export(raster2BIL)
export(select_PCA_components)
export(split_line)
import(raster)
import(tools)
importFrom(dissUtils,diss)
importFrom(ecodist,nmds)
importFrom(fields,rdist)
importFrom(future,multiprocess)
importFrom(future,plan)
importFrom(future,sequential)
importFrom(future.apply,future_lapply)
importFrom(labdsv,pco)
importFrom(matlab,ones)
importFrom(matlab,padarray)
importFrom(matrixStats,rowSds)
importFrom(methods,is)
importFrom(raster,projection)
importFrom(raster,raster)
importFrom(raster,shapefile)
......@@ -30,7 +33,13 @@ importFrom(raster,writeRaster)
importFrom(rgdal,readOGR)
importFrom(rgdal,writeOGR)
importFrom(snow,splitRows)
importFrom(sp,spTransform)
importFrom(stats,as.dist)
importFrom(stats,kmeans)
importFrom(stats,prcomp)
importFrom(stats,sd)
importFrom(stringr,str_count)
importFrom(utils,file.edit)
importFrom(utils,find)
importFrom(utils,read.table)
importFrom(vegan,fisher.alpha)
......@@ -12,160 +12,159 @@
#' converts a raster into BIL format as expected by DivMapping codes
#'
#' @param Raster.Path character. Full path for the raster to be converted
#' @param Sensor character. Name of the sensor
#' @param Convert.Integer boolean. Should data be converted into integer ?
#' @param Multiplying.Factor numeric. Multiplying factor (eg convert real reflectance values between 0 and 1 into integer between 0 and 10000).
#' @param Output.Directory character. Path to output directory.
#' @param Multiplying.Factor.Last numeric. Multiplying factor for last band.
#' @param Raster_Path character. Full path for the raster to be converted
#' @param Sensor character. Name of the sensor. a .hdr template for the sensor should be provided in extdata/HDR
#' @param Convert_Integer boolean. Should data be converted into integer ?
#' @param Multiplying_Factor numeric. Multiplying factor (eg convert real reflectance values between 0 and 1 into integer between 0 and 10000).
#' @param Output_Dir character. Path to output directory.
#' @param Multiplying_Factor_Last numeric. Multiplying factor for last band.
#'
#' @return Output.Path path for the image converted into ENVI BIL format
#' @return Output_Path path for the image converted into ENVI BIL format
#' @import raster
#' @import tools
#' @export
raster2BIL <- function(Raster.Path, Sensor = "unknown", Output.Directory = FALSE, Convert.Integer = TRUE, Multiplying.Factor = 1, Multiplying.Factor.Last = 1) {
raster2BIL <- function(Raster_Path, Sensor = "unknown", Output_Dir = FALSE, Convert_Integer = TRUE, Multiplying_Factor = 1, Multiplying_Factor_Last = 1) {
# get directory and file name of original image
Input.File <- basename(Raster.Path)
Input.Dir <- dirname(Raster.Path)
Input_File <- basename(Raster_Path)
Input_Dir <- dirname(Raster_Path)
# define path where data will be stored
if (Output.Directory == FALSE) {
Output.Path <- paste(Input.Dir, "/Converted_DivMAP/", file_path_sans_ext(Input.File), sep = "")
if (Output_Dir == FALSE) {
Output_Path <- paste(Input_Dir, "/Converted_DivMAP/", file_path_sans_ext(Input_File), sep = "")
} else {
dir.create(Output.Directory, showWarnings = FALSE, recursive = TRUE)
Output.Path <- paste(Output.Directory, "/", file_path_sans_ext(Input.File), sep = "")
dir.create(Output_Dir, showWarnings = FALSE, recursive = TRUE)
Output_Path <- paste(Output_Dir, "/", file_path_sans_ext(Input_File), sep = "")
}
message("The converted file will be written in the following location:")
print(Output.Path)
Output.Dir <- dirname(Output.Path)
dir.create(Output.Dir, showWarnings = FALSE, recursive = TRUE)
print(Output_Path)
Output_Dir <- dirname(Output_Path)
dir.create(Output_Dir, showWarnings = FALSE, recursive = TRUE)
# apply multiplying factors
message("reading initial file")
Output.Img <- Multiplying.Factor * brick(Raster.Path)
Last.Band.Name <- Output.Img@data@names[length(Output.Img@data@names)]
Output.Img[[Last.Band.Name]] <- Multiplying.Factor.Last * Output.Img[[Last.Band.Name]]
Output_Img <- Multiplying_Factor * brick(Raster_Path)
Last_Band_Name <- Output_Img@data@names[length(Output_Img@data@names)]
Output_Img[[Last_Band_Name]] <- Multiplying_Factor_Last * Output_Img[[Last_Band_Name]]
# convert into integer
if (Convert.Integer == TRUE) {
Output.Img <- round(Output.Img)
if (Convert_Integer == TRUE) {
Output_Img <- round(Output_Img)
}
# write raster
message("writing converted file")
if (Convert.Integer == TRUE) {
r <- writeRaster(Output.Img, filename = Output.Path, format = "EHdr", overwrite = TRUE, datatype = "INT2S")
if (Convert_Integer == TRUE) {
r <- writeRaster(Output_Img, filename = Output_Path, format = "EHdr", overwrite = TRUE, datatype = "INT2S")
} else {
r <- writeRaster(Output.Img, filename = Output.Path, format = "EHdr", overwrite = TRUE)
r <- writeRaster(Output_Img, filename = Output_Path, format = "EHdr", overwrite = TRUE)
}
hdr(r, format = "ENVI")
# remove unnecessary files
File2Remove <- paste(Output.Path, ".aux.xml", sep = "")
File2Remove2 <- paste(file_path_sans_ext(Output.Path), ".aux.xml", sep = "")
File2Remove <- paste(Output_Path, ".aux.xml", sep = "")
File2Remove2 <- paste(file_path_sans_ext(Output_Path), ".aux.xml", sep = "")
file.remove(File2Remove)
File2Remove <- paste(Output.Path, ".prj", sep = "")
File2Remove2 <- paste(file_path_sans_ext(Output.Path), ".prj", sep = "")
File2Remove <- paste(Output_Path, ".prj", sep = "")
File2Remove2 <- paste(file_path_sans_ext(Output_Path), ".prj", sep = "")
if (file.exists(File2Remove)) {
file.remove(File2Remove)
} else if (file.exists(File2Remove2)) {
file.remove(File2Remove2)
}
File2Remove <- paste(Output.Path, ".sta", sep = "")
File2Remove <- paste(file_path_sans_ext(Output.Path), ".sta", sep = "")
File2Remove <- paste(Output_Path, ".sta", sep = "")
File2Remove <- paste(file_path_sans_ext(Output_Path), ".sta", sep = "")
if (file.exists(File2Remove)) {
file.remove(File2Remove)
} else if (file.exists(File2Remove2)) {
file.remove(File2Remove2)
}
File2Remove <- paste(Output.Path, ".stx", sep = "")
File2Remove2 <- paste(file_path_sans_ext(Output.Path), ".stx", sep = "")
File2Remove <- paste(Output_Path, ".stx", sep = "")
File2Remove2 <- paste(file_path_sans_ext(Output_Path), ".stx", sep = "")
if (file.exists(File2Remove)) {
file.remove(File2Remove)
} else if (file.exists(File2Remove2)) {
file.remove(File2Remove2)
}
File2Rename <- paste(file_path_sans_ext(Output.Path), ".hdr", sep = "")
File2Rename2 <- paste(Output.Path, ".hdr", sep = "")
File2Rename <- paste(file_path_sans_ext(Output_Path), ".hdr", sep = "")
File2Rename2 <- paste(Output_Path, ".hdr", sep = "")
if (file.exists(File2Rename)) {
file.rename(from = File2Rename, to = File2Rename2)
}
# change dot into underscore
Output.Path.US <- file.path(
dirname(Output.Path),
gsub(basename(Output.Path), pattern = "[.]", replacement = "_")
Output_Path_US <- file.path(
dirname(Output_Path),
gsub(basename(Output_Path), pattern = "[.]", replacement = "_")
)
if (!Output.Path.US == Output.Path) {
file.rename(from = Output.Path, to = Output.Path.US)
if (!Output_Path_US == Output_Path) {
file.rename(from = Output_Path, to = Output_Path_US)
}
Output.Path.US.HDR <- paste0(Output.Path.US, ".hdr")
if (!Output.Path.US.HDR == paste0(Output.Path, ".hdr")) {
file.rename(from = paste0(Output.Path, ".hdr"), to = Output.Path.US.HDR)
Output_Path_US_HDR <- paste0(Output_Path_US, ".hdr")
if (!Output_Path_US_HDR == paste0(Output_Path, ".hdr")) {
file.rename(from = paste0(Output_Path, ".hdr"), to = Output_Path_US_HDR)
### UTILITY?? ###
file.rename(from = Output.Path.US.HDR, to = Output.Path.US.HDR)
file.rename(from = Output_Path_US_HDR, to = Output_Path_US_HDR)
}
if (!Sensor == "unknown") {
HDR.Temp.Path <- system.file("extdata", "HDR", paste0(Sensor, ".hdr"), package = "biodivMapR")
if (file.exists(HDR.Temp.Path)) {
HDR_Temp_Path <- system.file("extdata", "HDR", paste0(Sensor, ".hdr"), package = "biodivMapR")
if (file.exists(HDR_Temp_Path)) {
message("reading header template corresponding to the sensor located here:")
print(HDR.Temp.Path)
print(HDR_Temp_Path)
# get raster template corresponding to the sensor
HDR.Template <- read_ENVI_header(HDR.Temp.Path)
HDR_Template <- read_ENVI_header(HDR_Temp_Path)
# get info to update hdr file
# read hdr
HDR.input <- read_ENVI_header(get_HDR_name(Output.Path))
if (!is.null(HDR.Template$wavelength)) {
HDR.input$wavelength <- HDR.Template$wavelength
HDR_input <- read_ENVI_header(get_HDR_name(Output_Path))
if (!is.null(HDR_Template$wavelength)) {
HDR_input$wavelength <- HDR_Template$wavelength
}
if (!is.null(HDR.Template$`sensor type`)) {
HDR.input$`sensor type` <- HDR.Template$`sensor type`
if (!is.null(HDR_Template$`sensor type`)) {
HDR_input$`sensor type` <- HDR_Template$`sensor type`
}
if (!is.null(HDR.Template$`band names`)) {
HDR.input$`band names` <- HDR.Template$`band names`
if (!is.null(HDR_Template$`band names`)) {
HDR_input$`band names` <- HDR_Template$`band names`
}
if (!is.null(HDR.Template$`wavelength units`)) {
HDR.input$`wavelength units` <- HDR.Template$`wavelength units`
if (!is.null(HDR_Template$`wavelength units`)) {
HDR_input$`wavelength units` <- HDR_Template$`wavelength units`
}
# define visual stretch in the VIS domain
HDR.input$`default stretch` <- "0 1000 linear"
HDR_input$`default stretch` <- "0 1000 linear"
# write corresponding hdr file
write_ENVI_header(HDR.input, get_HDR_name(Output.Path))
} else if (!file.exists(HDR.Temp.Path)) {
write_ENVI_header(HDR_input, get_HDR_name(Output_Path))
} else if (!file.exists(HDR_Temp_Path)) {
message("Header template corresponding to the sensor expected to be found here")
print(HDR.Temp.Path)
print(HDR_Temp_Path)
message("please provide this header template in order to write info in HDR file")
print(get_HDR_name(Output.Path))
print(get_HDR_name(Output_Path))
message("or manually add wavelength location in HDR file, if relevant")
}
} else if (Sensor == "unknown") {
message("please make sure that the follozing header file contains information required")
print(get_HDR_name(Output.Path))
print(get_HDR_name(Output_Path))
message("or manually add wavelength location in HDR file, if relevant")
}
return(Output.Path)
return(Output_Path)
}
#' Checks if the data to be processed has the format type expected
#'
#' @param Raster.Path full path for the raster to be converted
#' @param Mask is the raster a mask?
#' @param Raster_Path character. full path for the raster to be converted
#' @param Mask boolean. Set true if the raster is a mask
#'
#' @return
#' @export
check_data <- function(Raster.Path, Mask = FALSE) {
HDR.Path <- get_HDR_name(Raster.Path)
check_data <- function(Raster_Path, Mask = FALSE) {
HDR_Path <- get_HDR_name(Raster_Path)
# check if the hdr file exists
if (file.exists(HDR.Path)) {
HDR <- read_ENVI_header(HDR.Path)
if (file.exists(HDR_Path)) {
HDR <- read_ENVI_header(HDR_Path)
if (Mask == FALSE & (!HDR$interleave == "bil") & (!HDR$interleave == "BIL")) {
message("")
message("*********************************************************")
message("The image format may not compatible with the processing chain")
message("Image format expected:")
......@@ -178,46 +177,55 @@ check_data <- function(Raster.Path, Mask = FALSE) {
message("in order to convert your raster data")
message("or use appropriate software")
message("*********************************************************")
message("")
stop()
} else if (Mask == FALSE & ((HDR$interleave == "bil") | (HDR$interleave == "BIL"))) {
if (HDR$`wavelength units` == "Unknown") {
message("")
message("*********************************************************")
message("Please make sure the wavelengths are in nanometers")
message("if not, stop processing and convert wavelengths in nanometers in HDR file")
message("*********************************************************")
message("")
}
if ((!HDR$`wavelength units` == "Nanometers") & (!HDR$`wavelength units` == "nanometers")) {
message("")
if (is.null(HDR$`wavelength units`)) {
message("*********************************************************")
message("Please make sure the wavelengths are in nanometers")
message("if not, stop processing and convert wavelengths in nanometers in HDR file")
message("Image to process is not multispectral/hyperspectral image ")
message("Format is OK, but make sure Continuum_Removal is set to FALSE")
message("*********************************************************")
message("")
} else {
if (HDR$`wavelength units` == "Unknown") {
message("*********************************************************")
message("IF MULTI / HYPERSPECTRAL DATA: ")
message("Please make sure the wavelengths are in nanometers")
message("if not, stop processing and convert wavelengths in nanometers in HDR file")
message("*********************************************************")
}
if ((!HDR$`wavelength units` == "Nanometers") & (!HDR$`wavelength units` == "nanometers")) {
message("*********************************************************")
message("IF MULTI / HYPERSPECTRAL DATA: ")
message("Please make sure the wavelengths are in nanometers")
message("if not, stop processing and convert wavelengths in nanometers in HDR file")
message("*********************************************************")
}
if (HDR$`wavelength units` == "micrometers") {
message("*********************************************************")
message("Please convert wavelengths in nanometers in HDR file")
message("*********************************************************")
stop()
}
if ((HDR$`wavelength units` == "nanometers") | (HDR$`wavelength units` == "Nanometers")) {
message("*********************************************************")
message(" Format of main raster OK for processing ")
message("*********************************************************")
}
}
if (HDR$`wavelength units` == "micrometers") {
message("")
message("*********************************************************")
message("Please convert wavelengths in nanometers in HDR file")
message("*********************************************************")
message("")
stop()
}
if ((HDR$`wavelength units` == "nanometers") | (HDR$`wavelength units` == "Nanometers")) {
message("")
message("*********************************************************")
message(" All information seem OK for image processing ")
message("*********************************************************")
message("")
}
}
} else if (Mask == TRUE & HDR$bands == 1 & ((HDR$interleave == "bil") | (HDR$interleave == "BIL") | (HDR$interleave == "bsq") | (HDR$interleave == "BSQ"))) {
message("*********************************************************")
message(" Format of mask raster OK for processing ")
message("*********************************************************")
} else if (Mask == TRUE & HDR$bands > 1) {
message("*********************************************************")
message(" Mask raster should contain only one layer ")
message(" Please produce a binary mask with a unique layer ")
message("*********************************************************")
stop()
}
} else {
message("")
message("*********************************************************")
message("The following HDR file was expected, but could not be found:")
print(HDR.Path)
print(HDR_Path)
message("The image format may not compatible with the processing chain")
message("Image format expected:")
message("ENVI hdr file with band interleaved by line (BIL) file format")
......@@ -226,7 +234,6 @@ check_data <- function(Raster.Path, Mask = FALSE) {
print("raster2BIL")
message("in order to convert your raster data")
message("*********************************************************")
message("")
stop()
}
return(invisible())
......
......@@ -11,7 +11,7 @@
# prepares data to run multithreaded continuum removal
#
# @param Spectral.Data initial data matrix (nb samples x nb bands)
# @param Spectral_Data initial data matrix (nb samples x nb bands)
# @param nbCPU
# @param Spectral information about spectral bands
#
......@@ -19,33 +19,33 @@
#' @importFrom snow splitRows
#' @importFrom future plan multiprocess sequential
#' @importFrom future.apply future_lapply
apply_continuum_removal <- function(Spectral.Data, Spectral, nbCPU = 1) {
apply_continuum_removal <- function(Spectral_Data, Spectral, nbCPU = 1) {
if (!length(Spectral$WaterVapor) == 0) {
Spectral.Data <- Spectral.Data[, -Spectral$WaterVapor]
Spectral_Data <- Spectral_Data[, -Spectral$WaterVapor]
}
# split data to perform continuum removal on into reasonable amount of data
nb.Values <- dim(Spectral.Data)[1] * dim(Spectral.Data)[2]
nb.Values <- dim(Spectral_Data)[1] * dim(Spectral_Data)[2]
if (nb.Values > 0) {
# corresponds to ~ 40 Mb data, but CR tends to requires ~ 10 times memory
# avoids memory crash
Max.nb.Values <- 2e6
nb.CR <- ceiling(nb.Values / Max.nb.Values)
Spectral.Data <- splitRows(Spectral.Data, nb.CR)
Spectral_Data <- splitRows(Spectral_Data, nb.CR)
# perform multithread continuum removal
plan(multiprocess, workers = nbCPU) ## Parallelize using four cores
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)
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)
plan(sequential)
Spectral.Data <- do.call("rbind", Spectral.Data.tmp)
rm(Spectral.Data.tmp)
Spectral_Data <- do.call("rbind", Spectral_Data_tmp)
rm(Spectral_Data_tmp)
} else {
# edit 31-jan-2018
# resize to delete first and last band as in continuum removal
Spectral.Data <- Spectral.Data[, -c(1, 2)]
Spectral_Data <- Spectral_Data[, -c(1, 2)]
}
gc()
return(Spectral.Data)
return(Spectral_Data)
}
# Computes continuum removal for matrix shaped data: more efficient than
......@@ -54,47 +54,47 @@ apply_continuum_removal <- function(Spectral.Data, Spectral, nbCPU = 1) {
# 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
# @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) {
ContinuumRemoval <- function(Minit, Spectral_Bands) {
# Filter and prepare data prior to continuum removal
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
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
# if samples to be considered
if (nb.Samples > 0) {
if (nbSamples > 0) {
# initialization:
# spectral band corresponding to each element of the data matrix
Lambda <- repmat(matrix(Spectral.Bands, nrow = 1), nb.Samples, 1)
Lambda <- repmat(matrix(Spectral_Bands, nrow = 1), nbSamples, 1)
# prepare matrices used to check evolution of the CR process
# - elements still not processed through continuum removal: initialization to 1
Still.Need.CR <- matrix(1, nrow = nb.Samples, ncol = nb.Bands)
Still.Need.CR <- matrix(1, nrow = nbSamples, ncol = nbBands)
# - value of the convex hull: initially set to 0
Convex.Hull <- matrix(0, nrow = nb.Samples, ncol = nb.Bands)
Convex_Hull <- matrix(0, nrow = nbSamples, ncol = nbBands)
# - reflectance value for latest interception with convex hull:
# initialization to value of the first reflectance measurement
Intercept.Hull <- repmat(matrix(Minit[, 1], ncol = 1), 1, nb.Bands)
Intercept_Hull <- repmat(matrix(Minit[, 1], ncol = 1), 1, nbBands)
# - spectral band of latest interception
Latest.Intercept <- repmat(matrix(Spectral.Bands[1], ncol = 1), nb.Samples, nb.Bands)
Latest.Intercept <- repmat(matrix(Spectral_Bands[1], ncol = 1), nbSamples, nbBands)
# number of spectral bands found as intercept
nb.Intercept <- 0
# continues until arbitrary stopping criterion:
# stops when reach last spectral band (all values before last = 0)
while (max(Still.Need.CR[, 1:(nb.Bands - 2)]) == 1 & (nb.Intercept <= (nb.Bands / 2))) {
while (max(Still.Need.CR[, 1:(nbBands - 2)]) == 1 & (nb.Intercept <= (nbBands / 2))) {
nb.Intercept <- nb.Intercept + 1
# Mstep give the position of the values to be updated
Update.Data <- matrix(1, nrow = nb.Samples, ncol = nb.Bands)
Update.Data[, nb.Bands] <- 0
Update_Data <- matrix(1, nrow = nbSamples, ncol = nbBands)
Update_Data[, nbBands] <- 0
# initial step: first column set to 0; following steps: all bands below
# max of the convex hull are set to 0
Update.Data[which((Lambda - Latest.Intercept) < 0)] <- 0
Update_Data[which((Lambda - Latest.Intercept) < 0)] <- 0
# compute slope for each coordinate
Slope <- (Minit - Intercept.Hull) / (Lambda - Latest.Intercept) * Still.Need.CR
Slope <- (Minit - Intercept_Hull) / (Lambda - Latest.Intercept) * Still.Need.CR
# set current spectral band and previous bands to -9999
if (!length(which(Still.Need.CR == 0)) == 0) {
Slope[which(Still.Need.CR == 0)] <- -9999
......@@ -103,25 +103,25 @@ ContinuumRemoval <- function(Minit, Spectral.Bands) {
Slope[which(is.na(Slope))] <- -9999
}
# get max index for each row and convert into linear index
Index.Max.Slope <- RowToLinear(max.col(Slope, ties.method = "last"), nb.Samples, nb.Bands)
Index.Max.Slope <- RowToLinear(max.col(Slope, ties.method = "last"), nbSamples, nbBands)
# !!!! OPTIM: replace repmat with column operation
# update coordinates of latest intercept
Latest.Intercept <- repmat(matrix(Lambda[Index.Max.Slope], ncol = 1), 1, nb.Bands)
Latest.Intercept <- repmat(matrix(Lambda[Index.Max.Slope], ncol = 1), 1, nbBands)
# update latest intercept
Intercept.Hull <- repmat(matrix(Minit[Index.Max.Slope], ncol = 1), 1, nb.Bands)
Intercept_Hull <- repmat(matrix(Minit[Index.Max.Slope], ncol = 1), 1, nbBands)
# values corresponding to the domain between the two continuum maxima
Update.Data[which((Lambda - Latest.Intercept) >= 0 | Latest.Intercept == Spectral.Bands[nb.Bands])] <- 0
Update_Data[which((Lambda - Latest.Intercept) >= 0 | Latest.Intercept == Spectral_Bands[nbBands])] <- 0
# values to eliminate for the next analysis: all spectral bands before latest intercept
Still.Need.CR[which((Lambda - Latest.Intercept) < 0)] <- 0
# the max slope is known, as well as the coordinates of the beginning and ending
# a matrix now has to be built
Convex.Hull <- Convex.Hull + Update.Data * (Intercept.Hull + sweep((Lambda - Latest.Intercept), 1, Slope[Index.Max.Slope], "*"))
Convex_Hull <- Convex_Hull + Update_Data * (Intercept_Hull + sweep((Lambda - Latest.Intercept), 1, Slope[Index.Max.Slope], "*"))
}
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
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
} else {
CR_Results <- matrix(0, ncol = (nb.Bands - 3), nrow = nb.Samples)
CR_Results <- matrix(0, ncol = (nbBands - 3), nrow = nbSamples)
}
list <- ls()
rm(list = list[-which(list == "CR_Results")])
......@@ -136,15 +136,15 @@ ContinuumRemoval <- function(Minit, Spectral.Bands) {
# - possibly remaining negative values are set to 0
# - constant spectra are eliminated
#
# @param Spectral.Bands
# @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) {