Skip to content

Commit

Permalink
Improves sits_bayes_postprocess function. Window can be of any size.
Browse files Browse the repository at this point in the history
  • Loading branch information
@rolfsimoes
rolfsimoes committed Aug 13, 2018
1 parent f5f1a01 commit 4dce4f0
Show file tree
Hide file tree
Showing 4 changed files with 31 additions and 119 deletions.
4 changes: 2 additions & 2 deletions R/RcppExports.R
View file
Expand Up @@ -13,8 +13,8 @@ preprocess_data <- function(data, minimum_value, scale_factor) {
.Call(`_sits_preprocess_data`, data, minimum_value, scale_factor)
}

smooth_estimator_class <- function(data, noise) {
.Call(`_sits_smooth_estimator_class`, data, noise)
smooth_estimator_class <- function(data, window, noise) {
.Call(`_sits_smooth_estimator_class`, data, window, noise)
}

scale_data <- function(data, scale_factor, adj_val = 0.0) {
Expand Down
9 changes: 7 additions & 2 deletions R/sits_bayes_postprocess.R
View file
Expand Up @@ -7,12 +7,16 @@
#' and a noise value, to do a bayesian smoothing process.
#'
#' @param raster_class output raster coverage
#' @param window (matrix) neighborhood window to compute bayesian smooth.
#' The central element index (i, j) is given by
#' i = floor(nrows(window)/2)+1 and j = floor(ncols(window)/2)+1.
#' Elements '0' are excluded from window.
#' @param noise bayesian smoothing parameter
#' @param file file to save the post processed raster
#' @return raster_layers.tb tibble with metadata about the output RasterLayer objects
#'
#' @export
sits_bayes_postprocess <- function(raster_class, noise, file) {
sits_bayes_postprocess <- function(raster_class, window = matrix(1, nrow = 3, ncol = 3, byrow = TRUE), noise = 100, file) {

# allocate matrix of each class
smooth_values <- matrix(NA,
Expand Down Expand Up @@ -43,7 +47,8 @@ sits_bayes_postprocess <- function(raster_class, noise, file) {
values[values > 9999] <- 9999

# apply bayes smooth estimator
smooth_values[, band] <- smooth_estimator_class(values, noise)
new_values <- smooth_estimator_class(data = values, window = window, noise = noise)
smooth_values[, band] <- new_values
}

# generate an output raster layer based on the input first layer
Expand Down
9 changes: 5 additions & 4 deletions src/RcppExports.cpp
View file
Expand Up @@ -45,14 +45,15 @@ BEGIN_RCPP
END_RCPP
}
// smooth_estimator_class
NumericVector smooth_estimator_class(const NumericMatrix& data, const double& noise);
RcppExport SEXP _sits_smooth_estimator_class(SEXP dataSEXP, SEXP noiseSEXP) {
NumericVector smooth_estimator_class(const NumericMatrix& data, const IntegerMatrix& window, const double& noise);
RcppExport SEXP _sits_smooth_estimator_class(SEXP dataSEXP, SEXP windowSEXP, SEXP noiseSEXP) {
BEGIN_RCPP
Rcpp::RObject rcpp_result_gen;
Rcpp::RNGScope rcpp_rngScope_gen;
Rcpp::traits::input_parameter< const NumericMatrix& >::type data(dataSEXP);
Rcpp::traits::input_parameter< const IntegerMatrix& >::type window(windowSEXP);
Rcpp::traits::input_parameter< const double& >::type noise(noiseSEXP);
rcpp_result_gen = Rcpp::wrap(smooth_estimator_class(data, noise));
rcpp_result_gen = Rcpp::wrap(smooth_estimator_class(data, window, noise));
return rcpp_result_gen;
END_RCPP
}
Expand Down Expand Up @@ -86,7 +87,7 @@ static const R_CallMethodDef CallEntries[] = {
{"_sits_apply_transition_matrix", (DL_FUNC) &_sits_apply_transition_matrix, 3},
{"_sits_normalize_data", (DL_FUNC) &_sits_normalize_data, 3},
{"_sits_preprocess_data", (DL_FUNC) &_sits_preprocess_data, 3},
{"_sits_smooth_estimator_class", (DL_FUNC) &_sits_smooth_estimator_class, 2},
{"_sits_smooth_estimator_class", (DL_FUNC) &_sits_smooth_estimator_class, 3},
{"_sits_scale_data", (DL_FUNC) &_sits_scale_data, 3},
{"_sits_scale_matrix_integer", (DL_FUNC) &_sits_scale_matrix_integer, 2},
{NULL, NULL, 0}
Expand Down
128 changes: 17 additions & 111 deletions src/prob_prior_neigh.cpp
View file
@@ -1,114 +1,32 @@
#include <Rcpp.h>
using namespace Rcpp;

NumericVector build_inside_neigh(const NumericMatrix& data,
const int& i, const int& j) {
NumericMatrix neigh(3, 3);
NumericVector build_neigh(const NumericMatrix& data,
const IntegerMatrix& window,
const int& i,
const int& j) {

for (int k = 0; k < 3; ++k) {
for (int l = 0; l < 3; ++l) {
neigh(k, l) = data(i + k - 1, j + l - 1);
}
}

return wrap(neigh);
}

NumericVector build_vborder_neigh(const NumericMatrix& data,
const int& i, const int& j) {
int nrows = data.nrow();
int ncols = data.ncol();
NumericMatrix neigh(3, 2);

// check left border
if (j == 0) {
for (int k = 0; k < 3; ++k) {
for (int l = 0; l < 2; ++l) {
neigh(k, l) = data(i + k - 1, j + l);
}
}
// check right border
} else {
for (int k = 0; k < 3; ++k) {
for (int l = 0; l < 2; ++l) {
neigh(k, l) = data(i + k - 1, j + l - 1);
}
}
}

return wrap(neigh);
}

NumericVector build_hborder_neigh(const NumericMatrix& data,
const int& i, const int& j) {
int nrows = data.nrow();
int ncols = data.ncol();
NumericMatrix neigh(2, 3);

// check upper border
if (i == 0) {
for (int k = 0; k < 2; ++k) {
for (int l = 0; l < 3; ++l) {
neigh(k, l) = data(i + k, j + l - 1);
}
}
// check lower border
} else {
for (int k = 0; k < 2; ++k) {
for (int l = 0; l < 3; ++l) {
neigh(k, l) = data(i + k - 1, j + l - 1);
}
}
}

return wrap(neigh);
}
NumericVector neigh;

NumericVector build_corner_neigh(const NumericMatrix& data,
const int& i, const int& j) {
int nrows = data.nrow();
int ncols = data.ncol();
NumericMatrix neigh(2, 2);
for (int k = 0; k < window.rows(); ++k) {
for (int l = 0; l < window.cols(); ++l) {
int data_i = i + k - window.rows() / 2, data_j = j + l - window.cols() / 2;

// check upper-left corner
if (i == 0 && j == 0) {
for (int k = 0; k < 2; ++k) {
for (int l = 0; l < 2; ++l) {
neigh(k, l) = data(i + k, j + l);
}
}
// check lower-left
} else if (i == nrows - 1 && j == 0) {
for (int k = 0; k < 2; ++k) {
for (int l = 0; l < 2; ++l) {
neigh(k, l) = data(i + k - 1, j + l);
}
}
// check upper-right
} else if (i == 0 && j == ncols - 1) {
for (int k = 0; k < 2; ++k) {
for (int l = 0; l < 2; ++l) {
neigh(k, l) = data(i + k, j + l - 1);
}
}
// check lower-right
} else if (i == nrows - 1 && j == ncols - 1) {
for (int k = 0; k < 2; ++k) {
for (int l = 0; l < 2; ++l) {
neigh(k, l) = data(i + k - 1, j + l - 1);
if (data_i >= 0 && data_j >= 0 &&
data_i < data.nrow() && data_j < data.ncol() && window(k, l) > 0) {
neigh.push_back(data(data_i, data_j) * window(k, l));
}
}
}

return wrap(neigh);
return neigh;
}

double smooth_estimator_pixel(const NumericVector& neigh,
const double& p,
double smooth_estimator_pixel(const double& p,
const NumericVector& neigh,
const double& noise) {

NumericVector log_neigh(neigh);
log_neigh = log(neigh / (10000 - neigh));
NumericVector log_neigh = log10(neigh / (10000 - neigh));
double x = log(p / (10000 - p));
double v = var(log_neigh);
double w1 = v / (noise + v);
Expand All @@ -120,30 +38,18 @@ double smooth_estimator_pixel(const NumericVector& neigh,

// [[Rcpp::export]]
NumericVector smooth_estimator_class(const NumericMatrix& data,
const IntegerMatrix& window,
const double& noise) {

int nrows = data.nrow();
int ncols = data.ncol();

NumericVector result(nrows * ncols);
NumericVector neigh;

int k = 0;
for (int i = 0; i < nrows; ++i) {
for (int j = 0; j < ncols; ++j) {
if (i > 0 && i < nrows - 1 && j > 0 && j < ncols - 1) {
neigh = build_inside_neigh(data, i, j);
} else if ((i == 0 && j == 0) ||
(i == nrows - 1 && j == 0) ||
(i == 0 && j == ncols - 1) ||
(i == nrows - 1 && j == ncols - 1)) {
neigh = build_corner_neigh(data, i, j);
} else if (i == 0 || i == nrows - 1) {
neigh = build_hborder_neigh(data, i, j);
} else if (j == 0 || j == ncols - 1) {
neigh = build_vborder_neigh(data, i, j);
}
result(k++) = smooth_estimator_pixel(neigh, data(i, j), noise);
result(k++) = smooth_estimator_pixel(data(i, j), build_neigh(data, window, i, j), noise);
}
}

Expand Down

0 comments on commit 4dce4f0

Please sign in to comment.