Commit 09ebd07a authored by Feret Jean-Baptiste's avatar Feret Jean-Baptiste

- fix bug when no shade mask provided and no radiometric filtering performed

- fix bug when continuum removal applied on pixels with constant values which were not filtered in previous stages
parent 647bed48
......@@ -66,21 +66,22 @@ ContinuumRemoval <- function(Minit, Spectral_Bands) {
CR_data$Minit <- c()
Spectral_Bands <- CR_data$Spectral_Bands
nbSamples <- CR_data$nbSamples
nbSamplesUpDate <- length(CR_data$SamplesToKeep)
# if samples to be considered
if (nbSamples > 0) {
# initialization:
# spectral band corresponding to each element of the data matrix
Lambda <- repmat(matrix(Spectral_Bands, nrow = 1), nbSamples, 1)
Lambda <- repmat(matrix(Spectral_Bands, nrow = 1), nbSamplesUpDate, 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 = nbSamples, ncol = nbBands)
Still.Need.CR <- matrix(1, nrow = nbSamplesUpDate, ncol = nbBands)
# - value of the convex hull: initially set to 0
Convex_Hull <- matrix(0, nrow = nbSamples, ncol = nbBands)
Convex_Hull <- matrix(0, nrow = nbSamplesUpDate, 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, nbBands)
# - spectral band of latest interception
Latest.Intercept <- repmat(matrix(Spectral_Bands[1], ncol = 1), nbSamples, nbBands)
Latest.Intercept <- repmat(matrix(Spectral_Bands[1], ncol = 1), nbSamplesUpDate, nbBands)
# number of spectral bands found as intercept
nb.Intercept <- 0
# continues until arbitrary stopping criterion:
......@@ -88,7 +89,7 @@ ContinuumRemoval <- function(Minit, Spectral_Bands) {
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 = nbSamples, ncol = nbBands)
Update_Data <- matrix(1, nrow = nbSamplesUpDate, 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
......@@ -103,7 +104,7 @@ 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"), nbSamples, nbBands)
Index.Max.Slope <- RowToLinear(max.col(Slope, ties.method = "last"), nbSamplesUpDate, nbBands)
# !!!! OPTIM: replace repmat with column operation
# update coordinates of latest intercept
Latest.Intercept <- repmat(matrix(Lambda[Index.Max.Slope], ncol = 1), 1, nbBands)
......@@ -163,8 +164,9 @@ filter_prior_CR <- function(Minit, Spectral_Bands) {
if (length(keep) == 1) {
Minit <- matrix(Minit, nrow = 1)
}
nbSamplesUpDate <- nrow(Minit)
# add negative values to the last column and update spectral bands
Minit <- cbind(Minit, matrix(-9999, ncol = 1, nrow = nbSamples))
Minit <- cbind(Minit, matrix(-9999, ncol = 1, nrow = nbSamplesUpDate))
nbBands <- ncol(Minit)
Spectral_Bands[nbBands] <- Spectral_Bands[nbBands - 1] + 10
my_list <- list("Minit" = Minit, "Spectral_Bands" = Spectral_Bands, "nbSamples" = nbSamples, "SamplesToKeep" = keep)
......
......@@ -233,7 +233,7 @@ filter_PCA <- function(ImPath, HDR, MaskPath, Shade_Update, Spectral, CR, PCA_mo
nbLines <- SeqRead_Shade$Lines_Per_Chunk[i]
lenBin <- SeqRead_Shade$ReadByte_End[i] - SeqRead_Shade$ReadByte_Start[i] + 1
ImgFormat <- "Shade"
if (!MaskPath == "") {
if ((!MaskPath == FALSE) & (!MaskPath == "")) {
Shade_Chunk <- read_image_subset(MaskPath, HDR_Shade, Byte_Start, lenBin, nbLines, Shade_Format, ImgFormat)
} else {
Shade_Chunk <- ones(nbLines * HDR$samples, 1)
......@@ -306,8 +306,12 @@ filter_PCA <- function(ImPath, HDR, MaskPath, Shade_Update, Spectral, CR, PCA_mo
write_PCA_raster <- function(ImPath, MaskPath, PCA_Path, PCA_model, Spectral, Nb_PCs, CR, TypePCA, nbCPU = 1, MaxRAM = 0.25) {
ImPathHDR <- get_HDR_name(ImPath)
HDR <- read_ENVI_header(ImPathHDR)
ShadeHDR <- get_HDR_name(MaskPath)
HDR_Shade <- read_ENVI_header(ShadeHDR)
if ((!MaskPath == FALSE) & (!MaskPath == "")) {
ShadeHDR <- get_HDR_name(MaskPath)
HDR_Shade <- read_ENVI_header(ShadeHDR)
} else {
HDR_Shade <- FALSE
}
# 1- create hdr and binary files corresponding to PCA file
HDR_PCA <- HDR
HDR_PCA$bands <- Nb_PCs
......@@ -337,7 +341,17 @@ write_PCA_raster <- function(ImPath, MaskPath, PCA_Path, PCA_model, Spectral, Nb
# read image file sequentially
Image_Format <- ENVI_type2bytes(HDR)
PCA_Format <- ENVI_type2bytes(HDR_PCA)
Shade_Format <- ENVI_type2bytes(HDR_Shade)
if (typeof(HDR_Shade) == 'list') {
Shade_Format <- ENVI_type2bytes(HDR_Shade)
} else if (typeof(HDR_Shade) == 'logical'){
Shade_Format <- FALSE
}
# if (!HDR_Shade == FALSE) {
# Shade_Format <- ENVI_type2bytes(HDR_Shade)
# } else {
# Shade_Format <- FALSE
# }
# prepare for sequential processing: SeqRead_Image informs about byte location to read
nbPieces <- split_image(HDR, LimitSizeGb = MaxRAM)
SeqRead_Image <- where_to_read(HDR, nbPieces)
......@@ -357,13 +371,16 @@ write_PCA_raster <- function(ImPath, MaskPath, PCA_Path, PCA_model, Spectral, Nb
Byte_Start <- SeqRead_Shade$ReadByte_Start[i]
nbLines <- SeqRead_Shade$Lines_Per_Chunk[i]
lenBin <- SeqRead_Shade$ReadByte_End[i] - SeqRead_Shade$ReadByte_Start[i] + 1
ImgFormat <- "Shade"
Shade_Chunk <- read_image_subset(MaskPath, HDR_Shade, Byte_Start, lenBin, nbLines, Shade_Format, ImgFormat)
elimShade <- which(Shade_Chunk == 0)
keepShade <- which(Shade_Chunk == 1)
Image_Chunk <- Image_Chunk[keepShade, ]
if (typeof(HDR_Shade) == 'list') {
ImgFormat <- "Shade"
Shade_Chunk <- read_image_subset(MaskPath, HDR_Shade, Byte_Start, lenBin, nbLines, Shade_Format, ImgFormat)
elimShade <- which(Shade_Chunk == 0)
keepShade <- which(Shade_Chunk == 1)
Image_Chunk <- Image_Chunk[keepShade, ]
} else {
keepShade <- matrix(1,ncol = 1,nrow = nrow(Image_Chunk))
}
# apply Continuum removal if needed
if (CR == TRUE) {
Image_Chunk <- apply_continuum_removal(Image_Chunk, Spectral, nbCPU = nbCPU)
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment