- update

land-info-lab · Oct 30, 2023 · 31a528c · 31a528c
1 parent 637ba68
commit 31a528c
Show file tree

Hide file tree

Showing 9 changed files with 110 additions and 62 deletions.
diff --git a/R/RcppExports.R b/R/RcppExports.R
@@ -5,8 +5,8 @@ get_depths <- function(px, py, x, y, num) {
     .Call(`_viewscape_get_depths`, px, py, x, y, num)
 }
 
-multiLabel <- function(vpts, dsm_list, max_dis, vpth, h, workers) {
-    .Call(`_viewscape_multiLabel`, vpts, dsm_list, max_dis, vpth, h, workers)
+multiLabel <- function(vpts, dsm, max_dis, vpth, h) {
+    .Call(`_viewscape_multiLabel`, vpts, dsm, max_dis, vpth, h)
 }
 
 visibleLabel <- function(viewpoint, dsm, h, max_dis) {

diff --git a/R/compute_viewshed.R b/R/compute_viewshed.R
@@ -118,38 +118,35 @@ compute_viewshed <- function(dsm,
       return(output)
     }
   }else if (multiviewpoints){
-    # set a new empty vector
-    viewsheds <- c()
+    #viewsheds <- c()
     if (parallel == TRUE){
       if (workers == 0) {
         stop("Please specify the number of CPU cores (workers)")
       }
+      inputs <- split(viewpoints,seq(nrow(viewpoints)))
       if (isTRUE(Sys.info()[1]=="Windows") == FALSE){
-        inputs <- split(viewpoints,seq(nrow(viewpoints)))
         bpparam <- BiocParallel::SnowParam(workers=workers, type="FORK")
-        suppressWarnings(
-          viewsheds <- BiocParallel::bplapply(X = inputs,
-                                              FUN = radius_viewshed,
-                                              dsm = dsm,
-                                              r = r,
-                                              offset = offset_viewpoint,
-                                              BPPARAM = bpparam)
-        )
       }else if (isTRUE(Sys.info()[1]=="Windows") == TRUE){
-        suppressWarnings(
-          viewsheds <- radius_viewshed_m(dsm=dsm,
-                                         r=r,
-                                         viewPts=viewpoints,
-                                         offset=offset_viewpoint,
-                                         workers=workers)
-        )
+        bpparam <- BiocParallel::SnowParam(workers=workers, type="SOCK")
       }
+      suppressWarnings(
+        viewsheds <- BiocParallel::bplapply(X = inputs,
+                                            FUN = radius_viewshed,
+                                            dsm = dsm,
+                                            r = r,
+                                            offset = offset_viewpoint,
+                                            BPPARAM = bpparam)
+      )
     } else {
-      for(i in 1:length(viewpoints[,1])){
-        viewpoint <- c(viewpoints[i,1],viewpoints[i,2])
-        output <- radius_viewshed(dsm, r, viewpoint, offset_viewpoint)
-        viewsheds <- c(viewsheds, output)
-      }
+      # for(i in 1:length(viewpoints[,1])){
+      #   viewpoint <- c(viewpoints[i,1],viewpoints[i,2])
+      #   output <- radius_viewshed(dsm, r, viewpoint, offset_viewpoint)
+      #   viewsheds <- c(viewsheds, output)
+      # }
+      viewsheds <- radius_viewshed_m(dsm=dsm,
+                                     r=r,
+                                     viewPts=viewpoints,
+                                     offset=offset_viewpoint)
     }
     return(viewsheds)
   }

diff --git a/R/utils.R b/R/utils.R
@@ -33,7 +33,7 @@ radius_viewshed <- function(dsm, r, viewPt, offset, offset2 = 0) {
 }
 
 #' @noMd
-radius_viewshed_m <- function(dsm, r, viewPts, offset, offset2 = 0, workers) {
+radius_viewshed_m <- function(dsm, r, viewPts, offset, offset2 = 0) {
   output <- c()
   dsm_list <- list()
   ex <- list()
@@ -58,8 +58,7 @@ radius_viewshed_m <- function(dsm, r, viewPts, offset, offset2 = 0, workers) {
                              dsm=dsm_list,
                              max_dis=distance,
                              vpth=offset,
-                             h=offset2,
-                             workers=workers)
+                             h=offset2)
   for(i in 1:length(z)) {
     out <- new("Viewshed",
                viewpoint = viewPts[i,],

diff --git a/R/zzz.R b/R/zzz.R
@@ -3,6 +3,8 @@
 .onLoad <- function(libname, pkgname) {
   Rcpp::sourceCpp(system.file("extdata/visibleLabel.cpp", package ="viewscape"),
                   env = asNamespace("viewscape"))
+  Rcpp::sourceCpp(system.file("extdata/multiLabel.cpp", package ="viewscape"),
+                  env = asNamespace("viewscape"))
   cat("████████████████████████████████████████████████████████████████████████████████████████████████████\n")
   cat("████████████████████████████████████████████████████████████████████████████████████████████████████\n")
   cat("██████████████████████████████████████████████MM/MMM/MMM/M///M██████████████████████████████████████\n")

diff --git a/inst/extdata/multiLabel.cpp b/inst/extdata/multiLabel.cpp
@@ -0,0 +1,64 @@
+#include <Rcpp.h>
+#include <iostream>
+#include <fstream>
+
+using namespace Rcpp;
+
+
+// [[Rcpp::export]]
+Rcpp::List multiLabel(Rcpp::NumericMatrix &vpts,
+                      Rcpp::List &dsm,
+                      const int max_dis,
+                      const double vpth,
+                      const double h) {
+
+  const int vptnum = vpts.rows();
+  Rcpp::List output(vptnum);
+  for (int i = 0; i < vptnum; i++) {
+    Rcpp::NumericMatrix sub_dsm = dsm[i];
+    int vx = vpts(i,0);
+    int vy = vpts(i,1);
+    const double vz = vpts(i,2) + vpth;
+    Rcpp::NumericVector zl;
+    Rcpp::NumericVector xl;
+    Rcpp::NumericVector yl;
+    const int sub_rows = sub_dsm.rows();
+    const int sub_cols = sub_dsm.cols();
+    int steps;
+    Rcpp::IntegerMatrix visible(sub_rows, sub_cols);
+    for (int j = 0; j < sub_rows; j++) {
+      for (int k = 0; k < sub_cols; k++) {
+        steps = sqrt((vx-k)*(vx-k) + (vy-j)*(vy-j));
+        const double z = sub_dsm(j,k) + h;
+        if (steps <= max_dis) {
+          Rcpp::NumericVector sequence(steps);
+          std::iota(sequence.begin(), sequence.end(), 1);
+          xl = vx + sequence * (k-vx)/steps;
+          yl = vy + sequence * (j-vy)/steps;
+          if(vz<z) {
+            zl = vz + sequence/steps*fabs(vz-z);
+          } else if (vz>z) {
+            zl = vz - sequence/steps*fabs(vz-z);
+          } else if (vz==z) {
+            zl = vz + 0*sequence;
+          }
+          int temp = 1;
+          for (int p = 0; p < steps; p++) {
+            const double d = zl[p] - sub_dsm(yl[p],xl[p]);
+            if (d < 0) {
+              temp = 0;
+              break;
+            }
+          }
+          visible(j,k) = temp;
+        } else {
+          visible(j,k) = 0;
+        }
+      }
+    }
+    output(i) = visible;
+  }
+  return output;
+}
+
+
diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp
@@ -26,18 +26,17 @@ BEGIN_RCPP
 END_RCPP
 }
 // multiLabel
-Rcpp::List multiLabel(Rcpp::NumericMatrix& vpts, Rcpp::List& dsm_list, const int max_dis, const double vpth, const double h, const int workers);
-RcppExport SEXP _viewscape_multiLabel(SEXP vptsSEXP, SEXP dsm_listSEXP, SEXP max_disSEXP, SEXP vpthSEXP, SEXP hSEXP, SEXP workersSEXP) {
+Rcpp::List multiLabel(Rcpp::NumericMatrix& vpts, Rcpp::List& dsm, const int max_dis, const double vpth, const double h);
+RcppExport SEXP _viewscape_multiLabel(SEXP vptsSEXP, SEXP dsmSEXP, SEXP max_disSEXP, SEXP vpthSEXP, SEXP hSEXP) {
 BEGIN_RCPP
     Rcpp::RObject rcpp_result_gen;
     Rcpp::RNGScope rcpp_rngScope_gen;
     Rcpp::traits::input_parameter< Rcpp::NumericMatrix& >::type vpts(vptsSEXP);
-    Rcpp::traits::input_parameter< Rcpp::List& >::type dsm_list(dsm_listSEXP);
+    Rcpp::traits::input_parameter< Rcpp::List& >::type dsm(dsmSEXP);
     Rcpp::traits::input_parameter< const int >::type max_dis(max_disSEXP);
     Rcpp::traits::input_parameter< const double >::type vpth(vpthSEXP);
     Rcpp::traits::input_parameter< const double >::type h(hSEXP);
-    Rcpp::traits::input_parameter< const int >::type workers(workersSEXP);
-    rcpp_result_gen = Rcpp::wrap(multiLabel(vpts, dsm_list, max_dis, vpth, h, workers));
+    rcpp_result_gen = Rcpp::wrap(multiLabel(vpts, dsm, max_dis, vpth, h));
     return rcpp_result_gen;
 END_RCPP
 }
@@ -58,7 +57,7 @@ END_RCPP
 
 static const R_CallMethodDef CallEntries[] = {
     {"_viewscape_get_depths", (DL_FUNC) &_viewscape_get_depths, 5},
-    {"_viewscape_multiLabel", (DL_FUNC) &_viewscape_multiLabel, 6},
+    {"_viewscape_multiLabel", (DL_FUNC) &_viewscape_multiLabel, 5},
     {"_viewscape_visibleLabel", (DL_FUNC) &_viewscape_visibleLabel, 4},
     {NULL, NULL, 0}
 };

diff --git a/src/RcppExports.o b/src/RcppExports.o
diff --git a/src/multiLabel.cpp b/src/multiLabel.cpp
@@ -1,49 +1,40 @@
 #include <Rcpp.h>
 #include <iostream>
 #include <fstream>
+
 using namespace Rcpp;
-#ifdef _WIN32
-  #include <omp.h>  // Include OpenMP for Windows
-#endif
 
 
 // [[Rcpp::export]]
 Rcpp::List multiLabel(Rcpp::NumericMatrix &vpts,
-                      Rcpp::List &dsm_list,
+                      Rcpp::List &dsm,
                       const int max_dis,
                       const double vpth,
-                      const double h,
-                      const int workers) {
+                      const double h) {
 
   const int vptnum = vpts.rows();
   Rcpp::List output(vptnum);
-  #ifdef _WIN32
-  #ifdef _OPENMP
-  #pragma omp parallel for num_threads(workers)
-  #else
-  #pragma omp parallel for num_threads(1)
-  #endif
-  for (unsigned int i = 0; i < vptnum; i++) {
-    const Rcpp::NumericMatrix sub_dsm = dsm_list[i];
+  for (int i = 0; i < vptnum; i++) {
+    Rcpp::NumericMatrix sub_dsm = dsm[i];
     int vx = vpts(i,0);
     int vy = vpts(i,1);
     const double vz = vpts(i,2) + vpth;
-
-    int sub_rows = sub_dsm.nrow();
-    int sub_cols = sub_dsm.ncol();
-
+    Rcpp::NumericVector zl;
+    Rcpp::NumericVector xl;
+    Rcpp::NumericVector yl;
+    const int sub_rows = sub_dsm.rows();
+    const int sub_cols = sub_dsm.cols();
     int steps;
     Rcpp::IntegerMatrix visible(sub_rows, sub_cols);
-    for (unsigned int j = 0; j < sub_rows; j++) {
-      for (unsigned int k = 0; k < sub_cols; k++) {
+    for (int j = 0; j < sub_rows; j++) {
+      for (int k = 0; k < sub_cols; k++) {
         steps = sqrt((vx-k)*(vx-k) + (vy-j)*(vy-j));
         const double z = sub_dsm(j,k) + h;
         if (steps <= max_dis) {
           Rcpp::NumericVector sequence(steps);
           std::iota(sequence.begin(), sequence.end(), 1);
-          Rcpp::NumericVector xl = vx + sequence * (k-vx)/steps;
-          Rcpp::NumericVector yl = vy + sequence * (j-vy)/steps;
-          Rcpp::NumericVector zl;
+          xl = vx + sequence * (k-vx)/steps;
+          yl = vy + sequence * (j-vy)/steps;
           if(vz<z) {
             zl = vz + sequence/steps*fabs(vz-z);
           } else if (vz>z) {
@@ -65,12 +56,8 @@ Rcpp::List multiLabel(Rcpp::NumericMatrix &vpts,
         }
       }
     }
-    #pragma omp critical
-    {
-      output[i] = visible;
-    }
+    output(i) = visible;
   }
-  #endif
   return output;
 }
 

diff --git a/src/viewscape.so b/src/viewscape.so