Skip to content

GitLab

Commit d68ec4e5 authored by jbferet's avatar jbferet
Browse files

change NbIter into nb_partitions

parent 4243b973
Hide whitespace changes
Inline Side-by-side
R/Lib_ImageProcess.R
View file @ d68ec4e5
......@@ -423,12 +423,12 @@ Extract.Samples.From.Image <- function(ImPath, coordPix, MaxRAM = FALSE, Progres
# @param ImPath path for image
# @param HDR path for hdr file
# @param ImPathShade path for shade mask
# @param NbIter number of iterations
# @param nb_partitions number of k-means then averaged
# @param Pix.Per.Iter number of pixels per iteration
#
# @return samples from image and updated number of pixels to sampel if necessary
Get.Random.Subset.From.Image <- function(ImPath, HDR, ImPathShade, NbIter, Pix.Per.Iter) {
nbPix2Sample <- NbIter * Pix.Per.Iter
Get.Random.Subset.From.Image <- function(ImPath, HDR, ImPathShade, nb_partitions, Pix.Per.Iter) {
nbPix2Sample <- nb_partitions * Pix.Per.Iter
# get total number of pixels
nbpix <- as.double(HDR$lines) * as.double(HDR$samples)
# 1- Exclude masked pixels from random subset
......@@ -452,9 +452,9 @@ Get.Random.Subset.From.Image <- function(ImPath, HDR, ImPathShade, NbIter, Pix.P
# if number of pixels to sample superior to number of valid pixels, then adjust iterations
if (NbValidPixels < nbPix2Sample) {
nbPix2Sample <- NbValidPixels
NbIter <- ceiling(NbValidPixels / Pix.Per.Iter)
Pix.Per.Iter <- floor(NbValidPixels / NbIter)
nbPix2Sample <- NbIter * Pix.Per.Iter
nb_partitions <- ceiling(NbValidPixels / Pix.Per.Iter)
Pix.Per.Iter <- floor(NbValidPixels / nb_partitions)
nbPix2Sample <- nb_partitions * Pix.Per.Iter
}
# Select a subset of nbPix2Sample among pixselected
pixselected <- ValidPixels[1:nbPix2Sample]
......
R/Lib_MapAlphaDiversity.R
View file @ d68ec4e5
......@@ -292,9 +292,9 @@ convert_PCA_to_SSD <- function(ReadWrite, Spectral.Species.Path, HDR.SS, HDR.SSD
compute_SSD <- function(Image.Chunk, window_size, nbclusters, MinSun, pcelim, Index.Alpha = "Shannon") {
nbi <- floor(dim(Image.Chunk)[1] / window_size)
nbj <- floor(dim(Image.Chunk)[2] / window_size)
nbIter <- dim(Image.Chunk)[3]
SSDMap <- array(NA, c(nbi, nbj, nbIter * nbclusters))
shannonIter <- FisherAlpha <- SimpsonAlpha <- array(NA, dim = c(nbi, nbj, nbIter))
nb_partitions <- dim(Image.Chunk)[3]
SSDMap <- array(NA, c(nbi, nbj, nb_partitions * nbclusters))
shannonIter <- FisherAlpha <- SimpsonAlpha <- array(NA, dim = c(nbi, nbj, nb_partitions))
PCsun <- matrix(NA, nrow = nbi, ncol = nbj)
# which spectral indices will be computed
......@@ -312,14 +312,14 @@ compute_SSD <- function(Image.Chunk, window_size, nbclusters, MinSun, pcelim, In
lj <- ((jj - 1) * window_size) + 1
uj <- jj * window_size
# put all iterations in a 2D matrix shape
ijit <- t(matrix(Image.Chunk[li:ui, lj:uj, ], ncol = nbIter))
ijit <- t(matrix(Image.Chunk[li:ui, lj:uj, ], ncol = nb_partitions))
# keep non zero values
ijit <- matrix(ijit[, which(!ijit[1, ] == 0)], nrow = nbIter)
ijit <- matrix(ijit[, which(!ijit[1, ] == 0)], nrow = nb_partitions)
nb.Pix.Sunlit <- dim(ijit)[2]
PCsun[ii, jj] <- nb.Pix.Sunlit / window_size**2
if (PCsun[ii, jj] > MinSun) {
# for each iteration
for (it in 1:nbIter) {
for (it in 1:nb_partitions) {
lbk <- (it - 1) * nbclusters
SSD <- as.vector(table(ijit[it, ]))
ClusterID <- sort(unique(ijit[it, ]))
......
R/Lib_MapBetaDiversity.R
View file @ d68ec4e5
......@@ -39,7 +39,7 @@ map_beta_div <- function(Input.Image.File, Output.Dir, window_size,
Output.Dir.SS <- Define.Output.SubDir(Output.Dir, Input.Image.File, TypePCA, "SpectralSpecies")
Output.Dir.BETA <- Define.Output.SubDir(Output.Dir, Input.Image.File, TypePCA, "BETA")
load(file = WS_Save)
Beta <- compute_beta_metrics(Output.Dir.SS, MinSun, Nb.Units.Ordin, NbIter,
Beta <- compute_beta_metrics(Output.Dir.SS, MinSun, Nb.Units.Ordin, nb_partitions,
nbclusters, pcelim, scaling = scaling,
nbCPU = nbCPU, MaxRAM = MaxRAM)
# Create images corresponding to Beta-diversity
......@@ -93,7 +93,7 @@ compute_NMDS <- function(MatBCdist) {
# @param SSD.Path ath for spectral species distribution file
# @param Sample.Sel Samples selected during ordination
# @param coordTotSort coordinates of sunlit spatial units
# @param NbIter number of iterations
# @param nb_partitions number of k-means then averaged
# @param nbclusters number of clusters
# @param pcelim number of CPUs available
# @param nbCPU number of CPUs available
......@@ -102,7 +102,7 @@ compute_NMDS <- function(MatBCdist) {
#' @importFrom snow splitRows
#' @importFrom future plan multiprocess sequential
#' @importFrom future.apply future_lapply
ordination_to_NN <- function(Beta.Ordination.sel, SSD.Path, Sample.Sel, coordTotSort, NbIter, nbclusters, pcelim, nbCPU = FALSE) {
ordination_to_NN <- function(Beta.Ordination.sel, SSD.Path, Sample.Sel, coordTotSort, nb_partitions, nbclusters, pcelim, nbCPU = FALSE) {
nb.Sunlit <- dim(coordTotSort)[1]
# define number of samples to be sampled each time during paralle processing
nb.samples.per.sub <- round(1e7 / dim(Sample.Sel)[1])
......@@ -117,7 +117,7 @@ ordination_to_NN <- function(Beta.Ordination.sel, SSD.Path, Sample.Sel, coordTot
OutPut <- future_lapply(id.sub,
FUN = ordination_parallel, coordTotSort = coordTotSort, SSD.Path = SSD.Path,
Sample.Sel = Sample.Sel, Beta.Ordination.sel = Beta.Ordination.sel, Nb.Units.Ordin = Nb.Units.Ordin,
NbIter = NbIter, nbclusters = nbclusters, pcelim = pcelim, future.scheduling = Schedule.Per.Thread,
nb_partitions = nb_partitions, nbclusters = nbclusters, pcelim = pcelim, future.scheduling = Schedule.Per.Thread,
future.packages = c("vegan", "dissUtils", "R.utils", "tools", "snow", "matlab")
)
plan(sequential)
......@@ -134,12 +134,12 @@ ordination_to_NN <- function(Beta.Ordination.sel, SSD.Path, Sample.Sel, coordTot
# @param Sample.Sel
# @param Beta.Ordination.sel
# @param Nb.Units.Ordin
# @param NbIter
# @param nb_partitions
# @param nbclusters
# @param pcelim
#
# @return list of mean and SD of alpha diversity metrics
ordination_parallel <- function(id.sub, coordTotSort, SSD.Path, Sample.Sel, Beta.Ordination.sel, Nb.Units.Ordin, NbIter, nbclusters, pcelim) {
ordination_parallel <- function(id.sub, coordTotSort, SSD.Path, Sample.Sel, Beta.Ordination.sel, Nb.Units.Ordin, nb_partitions, nbclusters, pcelim) {
# get Spectral species distribution
coordPix <- coordTotSort[id.sub, ]
......@@ -147,7 +147,7 @@ ordination_parallel <- function(id.sub, coordTotSort, SSD.Path, Sample.Sel, Beta
# compute the mean BC dissimilarity sequentially for each iteration
MatBCtmp <- matrix(0, nrow = nrow(id.sub), ncol = Nb.Units.Ordin)
SSDList <- list()
for (i in 1:NbIter) {
for (i in 1:nb_partitions) {
lb <- 1 + (i - 1) * nbclusters
ub <- i * nbclusters
SSDList[[1]] <- SSD.NN[, lb:ub]
......@@ -155,7 +155,7 @@ ordination_parallel <- function(id.sub, coordTotSort, SSD.Path, Sample.Sel, Beta
MatBCtmp0 <- compute_BCdiss(SSDList, pcelim)
MatBCtmp <- MatBCtmp + MatBCtmp0
}
MatBCtmp <- MatBCtmp / NbIter
MatBCtmp <- MatBCtmp / nb_partitions
# get the knn closest neighbors from each kernel
knn <- 3
OutPut <- compute_NN_from_ordination(MatBCtmp, knn, Beta.Ordination.sel)
......@@ -167,7 +167,7 @@ ordination_parallel <- function(id.sub, coordTotSort, SSD.Path, Sample.Sel, Beta
# @param Output.Dir directory where spectral species are stored
# @param MinSun minimum proportion of sunlit pixels required to consider plot
# @param Nb.Units.Ordin maximum number of spatial units to be processed in Ordination
# @param NbIter number of iterations
# @param nb_partitions number of iterations
# @param nbclusters number of clusters defined in k-Means
# @param scaling
# @param nbCPU
......@@ -176,7 +176,7 @@ ordination_parallel <- function(id.sub, coordTotSort, SSD.Path, Sample.Sel, Beta
#
# @return
#' @importFrom labdsv pco
compute_beta_metrics <- function(Output.Dir, MinSun, Nb.Units.Ordin, NbIter, nbclusters, pcelim, scaling = "PCO", nbCPU = FALSE, MaxRAM = FALSE) {
compute_beta_metrics <- function(Output.Dir, MinSun, Nb.Units.Ordin, nb_partitions, nbclusters, pcelim, scaling = "PCO", nbCPU = FALSE, MaxRAM = FALSE) {
# Define path for images to be used
SSD.Path <- paste(Output.Dir, "SpectralSpecies_Distribution", sep = "")
ImPathSunlit <- paste(Output.Dir, "SpectralSpecies_Distribution_Sunlit", sep = "")
......@@ -208,7 +208,7 @@ compute_beta_metrics <- function(Output.Dir, MinSun, Nb.Units.Ordin, NbIter, nbc
MatBC <- matrix(0, ncol = Nb.Units.Ordin, nrow = Nb.Units.Ordin)
SSDList <- list()
BC.from.SSD <- list()
for (i in 1:NbIter) {
for (i in 1:nb_partitions) {
lb <- 1 + (i - 1) * nbclusters
ub <- i * nbclusters
SSDList[[1]] <- Sample.Sel[, lb:ub]
......@@ -216,7 +216,7 @@ compute_beta_metrics <- function(Output.Dir, MinSun, Nb.Units.Ordin, NbIter, nbc
BC.from.SSD <- compute_BCdiss(SSDList, pcelim)
MatBC <- MatBC + BC.from.SSD
}
MatBC <- MatBC / NbIter
MatBC <- MatBC / nb_partitions
# Perform Ordination based on BC dissimilarity matrix
print("perform Ordination on the BC dissimilarity matrix averaged from all iterations")
......@@ -233,7 +233,7 @@ compute_beta_metrics <- function(Output.Dir, MinSun, Nb.Units.Ordin, NbIter, nbc
# Perform nearest neighbor on spatial units excluded from Ordination
print("BC dissimilarity between samples selected for Ordination and remaining")
coordTotSort <- Sunlit.Pixels$coordTotSort
Ordination.est <- ordination_to_NN(Beta.Ordination.sel, SSD.Path, Sample.Sel, coordTotSort, NbIter, nbclusters, pcelim, nbCPU = nbCPU)
Ordination.est <- ordination_to_NN(Beta.Ordination.sel, SSD.Path, Sample.Sel, coordTotSort, nb_partitions, nbclusters, pcelim, nbCPU = nbCPU)
# Reconstuct spatialized beta diversity map from previous computation
Sunlit.HDR <- Get.HDR.Name(ImPathSunlit)
......@@ -445,7 +445,7 @@ get_sunlit_pixels <- function(ImPathSunlit, MinSun) {
# Output.Dir.SS <- Define.Output.SubDir(Output.Dir, Input.Image.File, TypePCA, paste("SpectralSpecies_", nbclusters, sep = ""))
# Output.Dir.BETA <- Define.Output.SubDir(Output.Dir, Input.Image.File, TypePCA, paste("BETA_", nbclusters, sep = ""))
# load(file = WS_Save)
# Beta <- compute_beta_metrics(Output.Dir.SS, MinSun, Nb.Units.Ordin, NbIter, nbclusters, pcelim, scaling = scaling, nbCPU = nbCPU, MaxRAM = MaxRAM)
# Beta <- compute_beta_metrics(Output.Dir.SS, MinSun, Nb.Units.Ordin, nb_partitions, nbclusters, pcelim, scaling = scaling, nbCPU = nbCPU, MaxRAM = MaxRAM)
# # Create images corresponding to Beta-diversity
# print("Write beta diversity maps")
# Index <- paste("BetaDiversity_BCdiss_", scaling, sep = "")
......
R/Lib_MapSpectralSpecies.R
View file @ d68ec4e5
......@@ -63,7 +63,7 @@ map_spectral_species <- function(Input.Image.File, Output.Dir, PCA.Files, TypePC
## 2- PERFORM KMEANS FOR EACH ITERATION & DEFINE SPECTRAL SPECIES ##
print("perform k-means clustering for each subset and define centroids")
# scaling factor subPCA between 0 and 1
Kmeans.info <- init_kmeans(dataPCA, Pix.Per.Iter, NbIter, nbclusters, nbCPU)
Kmeans.info <- init_kmeans(dataPCA, Pix.Per.Iter, nb_partitions, nbclusters, nbCPU)
## 3- ASSIGN SPECTRAL SPECIES TO EACH PIXEL ##
print("apply Kmeans to the whole image and determine spectral species")
apply_kmeans(PCA.Files, PC.Select, ImPathShade, Kmeans.info, Spectral.Species.Path, nbCPU, MaxRAM)
......@@ -71,11 +71,11 @@ map_spectral_species <- function(Input.Image.File, Output.Dir, PCA.Files, TypePC
return()
}
# computes k-means from NbIter subsets taken from dataPCA
# computes k-means from nb_partitions subsets taken from dataPCA
#
# @param dataPCA initial dataset sampled from PCA image
# @param Pix.Per.Iter number of pixels per iteration
# @param NbIter number of iterations
# @param nb_partitions number of k-means then averaged
# @param nbCPU
# @param nbclusters number of clusters used in kmeans
#
......@@ -83,21 +83,21 @@ map_spectral_species <- function(Input.Image.File, Output.Dir, PCA.Files, TypePC
#' @importFrom future plan multiprocess sequential
#' @importFrom future.apply future_lapply
#' @importFrom stats kmeans
init_kmeans <- function(dataPCA, Pix.Per.Iter, NbIter, nbclusters, nbCPU = 1) {
init_kmeans <- function(dataPCA, Pix.Per.Iter, nb_partitions, nbclusters, nbCPU = 1) {
m0 <- apply(dataPCA, 2, function(x) min(x))
M0 <- apply(dataPCA, 2, function(x) max(x))
d0 <- M0 - m0
dataPCA <- Center.Reduce(dataPCA, m0, d0)
# get the dimensions of the images, and the number of subimages to process
dataPCA <- split(as.data.frame(dataPCA), rep(1:NbIter, each = Pix.Per.Iter))
dataPCA <- split(as.data.frame(dataPCA), rep(1:nb_partitions, each = Pix.Per.Iter))
# multiprocess kmeans clustering
plan(multiprocess, workers = nbCPU) ## Parallelize using four cores
Schedule.Per.Thread <- ceiling(length(dataPCA) / nbCPU)
res <- future_lapply(dataPCA, FUN = kmeans, centers = nbclusters, iter.max = 50, nstart = 10,
algorithm = c("Hartigan-Wong"), future.scheduling = Schedule.Per.Thread)
plan(sequential)
Centroids <- list(NbIter)
for (i in (1:NbIter)) {
Centroids <- list(nb_partitions)
for (i in (1:nb_partitions)) {
Centroids[[i]] <- res[[i]]$centers
}
my_list <- list("Centroids" = Centroids, "MinVal" = m0, "MaxVal" = M0, "Range" = d0)
......@@ -130,10 +130,10 @@ apply_kmeans <- function(PCA.Path, PC.Select, ImPathShade, Kmeans.info, Spectral
SeqRead.Shade <- Where.To.Read(HDR.Shade, nbPieces)
# create output file for spectral species assignment
HDR.SS <- HDR.PCA
nbIter <- length(Kmeans.info$Centroids)
HDR.SS$bands <- nbIter
nb_partitions <- length(Kmeans.info$Centroids)
HDR.SS$bands <- nb_partitions
HDR.SS$`data type` <- 1
Iter <- paste('Iter', 1:nbIter, collapse = ", ")
Iter <- paste('Iter', 1:nb_partitions, collapse = ", ")
HDR.SS$`band names` <- Iter
HDR.SS$wavelength <- NULL
HDR.SS$fwhm <- NULL
......@@ -158,8 +158,8 @@ apply_kmeans <- function(PCA.Path, PC.Select, ImPathShade, Kmeans.info, Spectral
Location.RW$lenBin.PCA <- SeqRead.PCA$ReadByte.End[i] - SeqRead.PCA$ReadByte.Start[i] + 1
Location.RW$Byte.Start.Shade <- SeqRead.Shade$ReadByte.Start[i]
Location.RW$lenBin.Shade <- SeqRead.Shade$ReadByte.End[i] - SeqRead.Shade$ReadByte.Start[i] + 1
Location.RW$Byte.Start.SS <- 1 + (SeqRead.Shade$ReadByte.Start[i] - 1) * nbIter
Location.RW$lenBin.SS <- nbIter * (SeqRead.Shade$ReadByte.End[i] - SeqRead.Shade$ReadByte.Start[i]) + 1
Location.RW$Byte.Start.SS <- 1 + (SeqRead.Shade$ReadByte.Start[i] - 1) * nb_partitions
Location.RW$lenBin.SS <- nb_partitions * (SeqRead.Shade$ReadByte.End[i] - SeqRead.Shade$ReadByte.Start[i]) + 1
compute_spectral_species(PCA.Path, ImPathShade, Spectral.Species.Path, Location.RW, PC.Select, Kmeans.info, nbCPU)
}
return()
......@@ -184,7 +184,7 @@ apply_kmeans <- function(PCA.Path, PC.Select, ImPathShade, Kmeans.info, Spectral
compute_spectral_species <- function(PCA.Path, ImPathShade, Spectral.Species.Path, Location.RW, PC.Select, Kmeans.info, nbCPU = 1) {
# characteristics of the centroids computed during preprocessing
nbIter <- length(Kmeans.info$Centroids)
nb_partitions <- length(Kmeans.info$Centroids)
nbCentroids <- nrow(Kmeans.info$Centroids[[1]])
CentroidsArray <- do.call("rbind", Kmeans.info$Centroids)
......@@ -216,7 +216,7 @@ compute_spectral_species <- function(PCA.Path, ImPathShade, Spectral.Species.Pat
SS.Format <- ENVI.Type2Bytes(HDR.SS)
# for each pixel in the subset, compute the nearest cluster for each iteration
Nearest.Cluster <- matrix(0, nrow = Location.RW$nbLines * HDR.PCA$samples, ncol = nbIter)
Nearest.Cluster <- matrix(0, nrow = Location.RW$nbLines * HDR.PCA$samples, ncol = nb_partitions)
# rdist consumes RAM --> divide data into pieces if too big and multiprocess
nbSamples.Per.Rdist <- 25000
if (length(keepShade) > 0) {
......@@ -225,14 +225,14 @@ compute_spectral_species <- function(PCA.Path, ImPathShade, Spectral.Species.Pat
plan(multiprocess, workers = nbCPU) ## Parallelize using four cores
Schedule.Per.Thread <- ceiling(nbSubsets / nbCPU)
res <- future_lapply(PCA.Chunk, FUN = RdistList, CentroidsArray = CentroidsArray, nbClusters = nrow(Kmeans.info$Centroids[[1]]), nbIter = nbIter, future.scheduling = Schedule.Per.Thread)
res <- future_lapply(PCA.Chunk, FUN = RdistList, CentroidsArray = CentroidsArray, nbClusters = nrow(Kmeans.info$Centroids[[1]]), nb_partitions = nb_partitions, future.scheduling = Schedule.Per.Thread)
plan(sequential)
res <- do.call("rbind", res)
Nearest.Cluster[keepShade, ] <- res
rm(res)
gc()
}
Nearest.Cluster <- array(Nearest.Cluster, c(Location.RW$nbLines, HDR.PCA$samples, nbIter))
Nearest.Cluster <- array(Nearest.Cluster, c(Location.RW$nbLines, HDR.PCA$samples, nb_partitions))
# Write spectral species distribution in the output file
fidSS <- file(
description = Spectral.Species.Path, open = "r+b", blocking = TRUE,
......@@ -249,22 +249,22 @@ compute_spectral_species <- function(PCA.Path, ImPathShade, Spectral.Species.Pat
return()
}
# Compute distance between each pixel of input data and each of the nbClusters x nbIter centroids
# Compute distance between each pixel of input data and each of the nbClusters x nb_partitions centroids
#
# @param InputData
# @param CentroidsArray
# @param nbClusters
# @param nbIter
# @param nb_partitions
#
# @return ResDist
#' @importFrom fields rdist
RdistList <- function(InputData, CentroidsArray, nbClusters, nbIter) {
RdistList <- function(InputData, CentroidsArray, nbClusters, nb_partitions) {
# number of pixels in input data
nbPixels <- nrow(InputData)
# compute distance between each pixel and each centroid
cluster.dist <- rdist(InputData, CentroidsArray)
# reshape distance into a matrix: all pixels from iteration 1, then all pixels from it2...
cluster.dist <- matrix(aperm(array(cluster.dist, c(nbPixels, nbClusters, nbIter)), c(1, 3, 2)), nrow = nbPixels * nbIter)
cluster.dist <- matrix(aperm(array(cluster.dist, c(nbPixels, nbClusters, nb_partitions)), c(1, 3, 2)), nrow = nbPixels * nb_partitions)
ResDist <- matrix(max.col(-cluster.dist), nrow = nbPixels)
return(ResDist)
}
R/Lib_PerformPCA.R
View file @ d68ec4e5
......@@ -17,13 +17,13 @@
#' @param TypePCA Type of PCA: "PCA" or "SPCA"
#' @param FilterPCA boolean. If TRUE 2nd filtering based on PCA
#' @param Excluded.WL boolean. Water Vapor Absorption domains (in nanometers). Can also be used to exclude spectific domains
#' @param NbIter numeric. Number of iteration to estimate diversity from the raster.
#' @param nb_partitions numeric. Number of repetitions to estimate diversity from the raster (averaging repetitions).
#' @param nbCPU numeric. Number fo CPUs in use.
#' @param MaxRAM numeric. Maximum size of chunk in GB to limit RAM allocation when reading image file.
#'
#' @return list of paths corresponding to resulting PCA files
#' @export
perform_PCA <- function(ImPath, ImPathShade, Output.Dir, Continuum.Removal = TRUE, TypePCA = "SPCA", FilterPCA = FALSE, Excluded.WL = FALSE, NbIter = 20, nbCPU = 1, MaxRAM = 0.25) {
perform_PCA <- function(ImPath, ImPathShade, Output.Dir, Continuum.Removal = TRUE, TypePCA = "SPCA", FilterPCA = FALSE, Excluded.WL = FALSE, nb_partitions = 20, nbCPU = 1, MaxRAM = 0.25) {
# define the path corresponding to image, mask and output directory
ImNames <- list()
ImNames$Input.Image <- ImPath
......@@ -34,12 +34,12 @@ perform_PCA <- function(ImPath, ImPathShade, Output.Dir, Continuum.Removal = TR
# Extract valid data subset and check validity
print("Extract pixels from the images to perform PCA on a subset")
# define number of pixels to be extracted from the image for each iteration
Pix.Per.Iter <- define_pixels_per_iter(ImNames, NbIter = NbIter)
nb.Pix.To.Sample <- NbIter * Pix.Per.Iter
Pix.Per.Iter <- define_pixels_per_iter(ImNames, nb_partitions = nb_partitions)
nb.Pix.To.Sample <- nb_partitions * Pix.Per.Iter
ImPathHDR <- Get.HDR.Name(ImPath)
HDR <- read.ENVI.header(ImPathHDR)
# extract a random selection of pixels from image
Subset <- Get.Random.Subset.From.Image(ImPath, HDR, ImPathShade, NbIter, Pix.Per.Iter)
Subset <- Get.Random.Subset.From.Image(ImPath, HDR, ImPathShade, nb_partitions, Pix.Per.Iter)
# if needed, apply continuum removal
if (Continuum.Removal == TRUE) {
Subset$DataSubset <- apply_continuum_removal(Subset$DataSubset, Spectral, nbCPU = nbCPU)
......@@ -47,9 +47,9 @@ perform_PCA <- function(ImPath, ImPathShade, Output.Dir, Continuum.Removal = TR
# if number of pixels available inferior number initial sample size
if (Subset$nbPix2Sample < nb.Pix.To.Sample) {
nb.Pix.To.Sample <- Subset$nbPix2Sample
NbIter <- ceiling(nb.Pix.To.Sample / Pix.Per.Iter)
Pix.Per.Iter <- floor(nb.Pix.To.Sample / NbIter)
nb.Pix.To.Sample <- NbIter * Pix.Per.Iter
nb_partitions <- ceiling(nb.Pix.To.Sample / Pix.Per.Iter)
Pix.Per.Iter <- floor(nb.Pix.To.Sample / nb_partitions)
nb.Pix.To.Sample <- nb_partitions * Pix.Per.Iter
}
DataSubset <- Subset$DataSubset
# clean reflectance data from inf and constant values
......@@ -89,7 +89,7 @@ perform_PCA <- function(ImPath, ImPathShade, Output.Dir, Continuum.Removal = TR
ImPathShade <- filter_PCA(ImPath, HDR, ImPathShade, Shade.Update, Spectral, Continuum.Removal, PCA.model, PCsel, TypePCA, nbCPU = nbCPU, MaxRAM = MaxRAM)
## Compute PCA 2 based on the updated shade mask ##
# extract a random selection of pixels from image
Subset <- Get.Random.Subset.From.Image(ImPath, HDR, ImPathShade, NbIter, Pix.Per.Iter)
Subset <- Get.Random.Subset.From.Image(ImPath, HDR, ImPathShade, nb_partitions, Pix.Per.Iter)
# if needed, apply continuum removal
if (Continuum.Removal == TRUE) {
Subset$DataSubset <- apply_continuum_removal(Subset$DataSubset, Spectral, nbCPU = nbCPU)
......@@ -97,9 +97,9 @@ perform_PCA <- function(ImPath, ImPathShade, Output.Dir, Continuum.Removal = TR
# if number of pixels available inferior number initial sample size
if (Subset$nbPix2Sample < nb.Pix.To.Sample) {
nb.Pix.To.Sample <- Subset$nbPix2Sample
NbIter <- ceiling(nb.Pix.To.Sample / Pix.Per.Iter)
Pix.Per.Iter <- floor(nb.Pix.To.Sample / NbIter)
nb.Pix.To.Sample <- NbIter * Pix.Per.Iter
nb_partitions <- ceiling(nb.Pix.To.Sample / Pix.Per.Iter)
Pix.Per.Iter <- floor(nb.Pix.To.Sample / nb_partitions)
nb.Pix.To.Sample <- nb_partitions * Pix.Per.Iter
}
DataSubset <- Subset$DataSubset
# # # assume that 1st data cleaning is enough...
......@@ -130,7 +130,7 @@ perform_PCA <- function(ImPath, ImPathShade, Output.Dir, Continuum.Removal = TR
write_PCA_raster(ImPath, ImPathShade, PCA.Path, PCA.model, Spectral, Nb.PCs, Continuum.Removal, TypePCA, nbCPU, MaxRAM = MaxRAM)
# save workspace for this stage
WS_Save <- paste(Output.Dir2, "PCA_Info.RData", sep = "")
save(PCA.model, Pix.Per.Iter, NbIter, ImPathShade, file = WS_Save)
save(PCA.model, Pix.Per.Iter, nb_partitions, ImPathShade, file = WS_Save)
PCA.Files <- PCA.Path
return(PCA.Files)
}
......@@ -558,10 +558,10 @@ pca <- function(X, type) {
# defines the number of pixels per iteration based on a trade-off between image size and sample size per iteration
#
# @param ImNames Path and name of the images to be processed
# @param NbIter number of iterations peformed to average diversity indices
# @param nb_partitions number of iterations peformed to average diversity indices
#
# @return Pix.Per.Iter number of pixels per iteration
define_pixels_per_iter <- function(ImNames, NbIter = NbIter) {
define_pixels_per_iter <- function(ImNames, nb_partitions = nb_partitions) {
ImPath <- ImNames$Input.Image
ImPathShade <- ImNames$Mask.list
# define dimensions of the image
......@@ -593,7 +593,7 @@ define_pixels_per_iter <- function(ImNames, NbIter = NbIter) {
# trade-off between number of pixels and total pixel size
# maximum number of pixels to be used
Max.Pixel.Per.Iter <- 250000
Max.Pixel.Per.Iter.Size <- NbIter * (Max.Pixel.Per.Iter * as.double(HDR$bands) * Image.Format$Bytes) / (1024^3)
Max.Pixel.Per.Iter.Size <- nb_partitions * (Max.Pixel.Per.Iter * as.double(HDR$bands) * Image.Format$Bytes) / (1024^3)
# maximum amount of data (Gb) to be used
Max.Size.Per.Iter <- 0.3
# amount of data available after first mask
......@@ -606,11 +606,11 @@ define_pixels_per_iter <- function(ImNames, NbIter = NbIter) {
Pix.Per.Iter <- Max.Pixel.Per.Iter
} else if (ImDataGb < Max.Pixel.Per.Iter.Size) {
# define number of pixels corresponding to ImDataGb
Pix.Per.Iter <- floor(Nb.Pixels.Sunlit / NbIter)
Pix.Per.Iter <- floor(Nb.Pixels.Sunlit / nb_partitions)
}
# if size too important, adjust number of pixels to match Max.Size.Per.Iter
} else if (Max.Pixel.Per.Iter.Size >= Max.Size.Per.Iter) {
Pix.Per.Iter <- floor((Max.Size.Per.Iter * (1024^3) / (as.double(HDR$bands) * Image.Format$Bytes)) / NbIter)
Pix.Per.Iter <- floor((Max.Size.Per.Iter * (1024^3) / (as.double(HDR$bands) * Image.Format$Bytes)) / nb_partitions)
}
return(Pix.Per.Iter)
}
......
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