diff --git a/DESCRIPTION b/DESCRIPTION index d8ea10ae..d0ee7f2a 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,7 +1,7 @@ Package: ggforce Type: Package Title: Accelerating 'ggplot2' -Version: 0.4.1.9000 +Version: 0.5.0 Authors@R: c(person(given = "Thomas Lin", family = "Pedersen", @@ -24,7 +24,6 @@ Depends: ggplot2 (>= 3.3.6), R (>= 3.3.0) Imports: - Rcpp (>= 0.12.2), grid, scales, MASS, @@ -41,7 +40,8 @@ Imports: cli, vctrs, systemfonts -LinkingTo: Rcpp, RcppEigen +LinkingTo: + cpp11 RoxygenNote: 7.2.3 Suggests: sessioninfo, @@ -51,7 +51,6 @@ Suggests: units (>= 0.8.0), covr Collate: - 'RcppExports.R' 'aaa.R' 'shape.R' 'arc_bar.R' @@ -64,6 +63,7 @@ Collate: 'bspline_closed.R' 'circle.R' 'concaveman.R' + 'cpp11.R' 'diagonal.R' 'diagonal_wide.R' 'ellipse.R' diff --git a/NAMESPACE b/NAMESPACE index b7e7252c..d648dc93 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -161,7 +161,6 @@ import(ggplot2) import(rlang) import(vctrs) importFrom(MASS,fractions) -importFrom(Rcpp,sourceCpp) importFrom(ggplot2,label_parsed) importFrom(ggplot2,layer) importFrom(grDevices,chull) @@ -241,3 +240,4 @@ importFrom(tweenr,tween_t) importFrom(utils,packageVersion) importFrom(withr,with_seed) useDynLib(ggforce) +useDynLib(ggforce, .registration = TRUE) diff --git a/NEWS.md b/NEWS.md index 310275cf..1ec164a7 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,4 +1,4 @@ -# ggforce (development version) +# ggforce 0.5.0 * Fixed a bug that would cause reordering of data in some geoms (#314) * The concaveman package is no longer a dependency for `geom_mark_hull()` (#308) diff --git a/R/RcppExports.R b/R/RcppExports.R deleted file mode 100644 index d97af41e..00000000 --- a/R/RcppExports.R +++ /dev/null @@ -1,35 +0,0 @@ -# Generated by using Rcpp::compileAttributes() -> do not edit by hand -# Generator token: 10BE3573-1514-4C36-9D1C-5A225CD40393 - -splinePath <- function(x, y, degree, knots, detail, type) { - .Call('_ggforce_splinePath', PACKAGE = 'ggforce', x, y, degree, knots, detail, type) -} - -getSplines <- function(x, y, id, detail, type = "clamped") { - .Call('_ggforce_getSplines', PACKAGE = 'ggforce', x, y, id, detail, type) -} - -bezierPath <- function(x, y, detail) { - .Call('_ggforce_bezierPath', PACKAGE = 'ggforce', x, y, detail) -} - -getBeziers <- function(x, y, id, detail) { - .Call('_ggforce_getBeziers', PACKAGE = 'ggforce', x, y, id, detail) -} - -concaveman_c <- function(p, h, concavity, threshold) { - .Call('_ggforce_concaveman_c', PACKAGE = 'ggforce', p, h, concavity, threshold) -} - -enclose_ellip_points <- function(x, y, id, tol) { - .Call('_ggforce_enclose_ellip_points', PACKAGE = 'ggforce', x, y, id, tol) -} - -enclose_points <- function(x, y, id) { - .Call('_ggforce_enclose_points', PACKAGE = 'ggforce', x, y, id) -} - -points_to_path <- function(pos, path, close) { - .Call('_ggforce_points_to_path', PACKAGE = 'ggforce', pos, path, close) -} - diff --git a/R/cpp11.R b/R/cpp11.R new file mode 100644 index 00000000..4f4a7cb8 --- /dev/null +++ b/R/cpp11.R @@ -0,0 +1,33 @@ +# Generated by cpp11: do not edit by hand + +splinePath <- function(x, y, degree, knots, detail, type) { + .Call(`_ggforce_splinePath`, x, y, degree, knots, detail, type) +} + +getSplines <- function(x, y, id, detail, type) { + .Call(`_ggforce_getSplines`, x, y, id, detail, type) +} + +bezierPath <- function(x, y, detail) { + .Call(`_ggforce_bezierPath`, x, y, detail) +} + +getBeziers <- function(x, y, id, detail) { + .Call(`_ggforce_getBeziers`, x, y, id, detail) +} + +concaveman_c <- function(p, h, concavity, threshold) { + .Call(`_ggforce_concaveman_c`, p, h, concavity, threshold) +} + +enclose_ellip_points <- function(x, y, id, tol) { + .Call(`_ggforce_enclose_ellip_points`, x, y, id, tol) +} + +enclose_points <- function(x, y, id) { + .Call(`_ggforce_enclose_points`, x, y, id) +} + +points_to_path <- function(pos, path, close) { + .Call(`_ggforce_points_to_path`, pos, path, close) +} diff --git a/R/ggforce-package.R b/R/ggforce-package.R index 56d13987..39569a9b 100644 --- a/R/ggforce-package.R +++ b/R/ggforce-package.R @@ -1,6 +1,5 @@ #' @useDynLib ggforce #' @import ggplot2 -#' @importFrom Rcpp sourceCpp #' #' @examples #' rocketData <- data.frame( @@ -57,8 +56,9 @@ "_PACKAGE" ## usethis namespace: start -#' @importFrom lifecycle deprecated #' @import rlang #' @import vctrs +#' @importFrom lifecycle deprecated +#' @useDynLib ggforce, .registration = TRUE ## usethis namespace: end NULL diff --git a/cran-comments.md b/cran-comments.md index e2f052e1..4d3b7aa5 100644 --- a/cran-comments.md +++ b/cran-comments.md @@ -1 +1,10 @@ -Quick release to address the sanitizer issues reported by CRAN +Bigger patch release addressing a range of bugs, as well as removing the +concaveman, Rcpp, and RcppEigen dependency in favor of cpp11 + +## revdepcheck results + +We checked 82 reverse dependencies, comparing R CMD check results across CRAN and dev versions of this package. + + * We saw 0 new problems + * We failed to check 0 packages + diff --git a/man/figures/README-example-1.png b/man/figures/README-example-1.png index 6fb06a88..c4de7686 100644 Binary files a/man/figures/README-example-1.png and b/man/figures/README-example-1.png differ diff --git a/revdep/README.md b/revdep/README.md index 8f229790..52e12c71 100644 --- a/revdep/README.md +++ b/revdep/README.md @@ -1,14 +1,2 @@ # Revdeps -## Failed to check (1) - -|package |version |error |warning |note | -|:-------|:-------|:-----|:-------|:----| -|NA |? | | | | - -## New problems (1) - -|package |version |error |warning |note | -|:----------------------------|:-------|:-----|:-------|:----| -|[cubble](problems.md#cubble) |0.1.1 | |__+1__ | | - diff --git a/revdep/cran.md b/revdep/cran.md index a2c1f748..4e76489a 100644 --- a/revdep/cran.md +++ b/revdep/cran.md @@ -1,15 +1,7 @@ ## revdepcheck results -We checked 66 reverse dependencies (65 from CRAN + 1 from Bioconductor), comparing R CMD check results across CRAN and dev versions of this package. +We checked 82 reverse dependencies, comparing R CMD check results across CRAN and dev versions of this package. - * We saw 1 new problems + * We saw 0 new problems * We failed to check 0 packages -Issues with CRAN packages are summarised below. - -### New problems -(This reports the first line of each new failure) - -* cubble - checking re-building of vignette outputs ... WARNING - diff --git a/revdep/failures.md b/revdep/failures.md index c6abf32b..9a207363 100644 --- a/revdep/failures.md +++ b/revdep/failures.md @@ -1,35 +1 @@ -# NA - -
- -* Version: NA -* GitHub: NA -* Source code: https://github.com/cran/NA -* Number of recursive dependencies: 0 - -Run `cloud_details(, "NA")` for more info - -
- -## Error before installation - -### Devel - -``` - - - - - - -``` -### CRAN - -``` - - - - - - -``` +*Wow, no problems at all. :)* \ No newline at end of file diff --git a/revdep/problems.md b/revdep/problems.md index 72167936..9a207363 100644 --- a/revdep/problems.md +++ b/revdep/problems.md @@ -1,41 +1 @@ -# cubble - -
- -* Version: 0.1.1 -* GitHub: https://github.com/huizezhang-sherry/cubble -* Source code: https://github.com/cran/cubble -* Date/Publication: 2022-06-02 12:30:06 UTC -* Number of recursive dependencies: 136 - -Run `cloud_details(, "cubble")` for more info - -
- -## Newly broken - -* checking re-building of vignette outputs ... WARNING - ``` - Error(s) in re-building vignettes: - --- re-building ‘aggregation.Rmd’ using rmarkdown - `summarise()` has grouped output by 'id'. You can override using the `.groups` - argument. - `summarise()` has grouped output by 'cluster', 'id'. You can override using the - `.groups` argument. - `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")' - Adding missing grouping variables: `id` - Quitting from lines 123-131 (aggregation.Rmd) - Error: processing vignette 'aggregation.Rmd' failed with diagnostics: - ... - The key variable is named differently in the two datasets. Coerce the key to id - to bind them together. - Only bind the common variables from both datasets. - --- finished re-building ‘matching.Rmd’ - - SUMMARY: processing the following file failed: - ‘aggregation.Rmd’ - - Error: Vignette re-building failed. - Execution halted - ``` - +*Wow, no problems at all. :)* \ No newline at end of file diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp deleted file mode 100644 index b6612821..00000000 --- a/src/RcppExports.cpp +++ /dev/null @@ -1,142 +0,0 @@ -// Generated by using Rcpp::compileAttributes() -> do not edit by hand -// Generator token: 10BE3573-1514-4C36-9D1C-5A225CD40393 - -#include -#include - -using namespace Rcpp; - -#ifdef RCPP_USE_GLOBAL_ROSTREAM -Rcpp::Rostream& Rcpp::Rcout = Rcpp::Rcpp_cout_get(); -Rcpp::Rostream& Rcpp::Rcerr = Rcpp::Rcpp_cerr_get(); -#endif - -// splinePath -NumericMatrix splinePath(NumericVector x, NumericVector y, int degree, std::vector knots, int detail, std::string type); -RcppExport SEXP _ggforce_splinePath(SEXP xSEXP, SEXP ySEXP, SEXP degreeSEXP, SEXP knotsSEXP, SEXP detailSEXP, SEXP typeSEXP) { -BEGIN_RCPP - Rcpp::RObject rcpp_result_gen; - Rcpp::RNGScope rcpp_rngScope_gen; - Rcpp::traits::input_parameter< NumericVector >::type x(xSEXP); - Rcpp::traits::input_parameter< NumericVector >::type y(ySEXP); - Rcpp::traits::input_parameter< int >::type degree(degreeSEXP); - Rcpp::traits::input_parameter< std::vector >::type knots(knotsSEXP); - Rcpp::traits::input_parameter< int >::type detail(detailSEXP); - Rcpp::traits::input_parameter< std::string >::type type(typeSEXP); - rcpp_result_gen = Rcpp::wrap(splinePath(x, y, degree, knots, detail, type)); - return rcpp_result_gen; -END_RCPP -} -// getSplines -List getSplines(NumericVector x, NumericVector y, IntegerVector id, int detail, std::string type); -RcppExport SEXP _ggforce_getSplines(SEXP xSEXP, SEXP ySEXP, SEXP idSEXP, SEXP detailSEXP, SEXP typeSEXP) { -BEGIN_RCPP - Rcpp::RObject rcpp_result_gen; - Rcpp::RNGScope rcpp_rngScope_gen; - Rcpp::traits::input_parameter< NumericVector >::type x(xSEXP); - Rcpp::traits::input_parameter< NumericVector >::type y(ySEXP); - Rcpp::traits::input_parameter< IntegerVector >::type id(idSEXP); - Rcpp::traits::input_parameter< int >::type detail(detailSEXP); - Rcpp::traits::input_parameter< std::string >::type type(typeSEXP); - rcpp_result_gen = Rcpp::wrap(getSplines(x, y, id, detail, type)); - return rcpp_result_gen; -END_RCPP -} -// bezierPath -NumericMatrix bezierPath(NumericVector x, NumericVector y, int detail); -RcppExport SEXP _ggforce_bezierPath(SEXP xSEXP, SEXP ySEXP, SEXP detailSEXP) { -BEGIN_RCPP - Rcpp::RObject rcpp_result_gen; - Rcpp::RNGScope rcpp_rngScope_gen; - Rcpp::traits::input_parameter< NumericVector >::type x(xSEXP); - Rcpp::traits::input_parameter< NumericVector >::type y(ySEXP); - Rcpp::traits::input_parameter< int >::type detail(detailSEXP); - rcpp_result_gen = Rcpp::wrap(bezierPath(x, y, detail)); - return rcpp_result_gen; -END_RCPP -} -// getBeziers -List getBeziers(NumericVector x, NumericVector y, IntegerVector id, int detail); -RcppExport SEXP _ggforce_getBeziers(SEXP xSEXP, SEXP ySEXP, SEXP idSEXP, SEXP detailSEXP) { -BEGIN_RCPP - Rcpp::RObject rcpp_result_gen; - Rcpp::RNGScope rcpp_rngScope_gen; - Rcpp::traits::input_parameter< NumericVector >::type x(xSEXP); - Rcpp::traits::input_parameter< NumericVector >::type y(ySEXP); - Rcpp::traits::input_parameter< IntegerVector >::type id(idSEXP); - Rcpp::traits::input_parameter< int >::type detail(detailSEXP); - rcpp_result_gen = Rcpp::wrap(getBeziers(x, y, id, detail)); - return rcpp_result_gen; -END_RCPP -} -// concaveman_c -NumericMatrix concaveman_c(NumericMatrix p, IntegerVector h, double concavity, double threshold); -RcppExport SEXP _ggforce_concaveman_c(SEXP pSEXP, SEXP hSEXP, SEXP concavitySEXP, SEXP thresholdSEXP) { -BEGIN_RCPP - Rcpp::RObject rcpp_result_gen; - Rcpp::RNGScope rcpp_rngScope_gen; - Rcpp::traits::input_parameter< NumericMatrix >::type p(pSEXP); - Rcpp::traits::input_parameter< IntegerVector >::type h(hSEXP); - Rcpp::traits::input_parameter< double >::type concavity(concavitySEXP); - Rcpp::traits::input_parameter< double >::type threshold(thresholdSEXP); - rcpp_result_gen = Rcpp::wrap(concaveman_c(p, h, concavity, threshold)); - return rcpp_result_gen; -END_RCPP -} -// enclose_ellip_points -DataFrame enclose_ellip_points(NumericVector x, NumericVector y, IntegerVector id, double tol); -RcppExport SEXP _ggforce_enclose_ellip_points(SEXP xSEXP, SEXP ySEXP, SEXP idSEXP, SEXP tolSEXP) { -BEGIN_RCPP - Rcpp::RObject rcpp_result_gen; - Rcpp::RNGScope rcpp_rngScope_gen; - Rcpp::traits::input_parameter< NumericVector >::type x(xSEXP); - Rcpp::traits::input_parameter< NumericVector >::type y(ySEXP); - Rcpp::traits::input_parameter< IntegerVector >::type id(idSEXP); - Rcpp::traits::input_parameter< double >::type tol(tolSEXP); - rcpp_result_gen = Rcpp::wrap(enclose_ellip_points(x, y, id, tol)); - return rcpp_result_gen; -END_RCPP -} -// enclose_points -DataFrame enclose_points(NumericVector x, NumericVector y, IntegerVector id); -RcppExport SEXP _ggforce_enclose_points(SEXP xSEXP, SEXP ySEXP, SEXP idSEXP) { -BEGIN_RCPP - Rcpp::RObject rcpp_result_gen; - Rcpp::RNGScope rcpp_rngScope_gen; - Rcpp::traits::input_parameter< NumericVector >::type x(xSEXP); - Rcpp::traits::input_parameter< NumericVector >::type y(ySEXP); - Rcpp::traits::input_parameter< IntegerVector >::type id(idSEXP); - rcpp_result_gen = Rcpp::wrap(enclose_points(x, y, id)); - return rcpp_result_gen; -END_RCPP -} -// points_to_path -List points_to_path(NumericMatrix pos, ListOf path, bool close); -RcppExport SEXP _ggforce_points_to_path(SEXP posSEXP, SEXP pathSEXP, SEXP closeSEXP) { -BEGIN_RCPP - Rcpp::RObject rcpp_result_gen; - Rcpp::RNGScope rcpp_rngScope_gen; - Rcpp::traits::input_parameter< NumericMatrix >::type pos(posSEXP); - Rcpp::traits::input_parameter< ListOf >::type path(pathSEXP); - Rcpp::traits::input_parameter< bool >::type close(closeSEXP); - rcpp_result_gen = Rcpp::wrap(points_to_path(pos, path, close)); - return rcpp_result_gen; -END_RCPP -} - -static const R_CallMethodDef CallEntries[] = { - {"_ggforce_splinePath", (DL_FUNC) &_ggforce_splinePath, 6}, - {"_ggforce_getSplines", (DL_FUNC) &_ggforce_getSplines, 5}, - {"_ggforce_bezierPath", (DL_FUNC) &_ggforce_bezierPath, 3}, - {"_ggforce_getBeziers", (DL_FUNC) &_ggforce_getBeziers, 4}, - {"_ggforce_concaveman_c", (DL_FUNC) &_ggforce_concaveman_c, 4}, - {"_ggforce_enclose_ellip_points", (DL_FUNC) &_ggforce_enclose_ellip_points, 4}, - {"_ggforce_enclose_points", (DL_FUNC) &_ggforce_enclose_points, 3}, - {"_ggforce_points_to_path", (DL_FUNC) &_ggforce_points_to_path, 3}, - {NULL, NULL, 0} -}; - -RcppExport void R_init_ggforce(DllInfo *dll) { - R_registerRoutines(dll, NULL, CallEntries, NULL, NULL); - R_useDynamicSymbols(dll, FALSE); -} diff --git a/src/bSpline.cpp b/src/bSpline.cpp index ae0c7ea3..f1ebe836 100644 --- a/src/bSpline.cpp +++ b/src/bSpline.cpp @@ -1,18 +1,27 @@ -#include #include "deBoor.h" -using namespace Rcpp; -std::vector createKnots(int nControl, int degree) { +#include +#include +#include +#include +#include + +using namespace cpp11::literals; + +#include + +cpp11::writable::doubles createKnots(int nControl, int degree) { int nKnots = nControl + degree + 1; - std::vector knots (nKnots, 0); + cpp11::writable::doubles knots; + knots.reserve(nKnots); for (int i = 0; i < nKnots; i++) { if (i < degree + 1) { - knots[i] = 0; + knots.push_back(0); } else if (i < nKnots - degree) { - knots[i] = knots[i-1] + 1; + knots.push_back(knots[i-1] + 1); } else { - knots[i] = knots[i-1]; + knots.push_back(knots[i-1]); } } return knots; @@ -27,7 +36,7 @@ std::vector createOpenKnots(int nControl, int degree) { } return knots; } -std::vector createControls(NumericVector x, NumericVector y) { +std::vector createControls(const cpp11::doubles& x, const cpp11::doubles& y) { int nControls = x.size(); std::vector controls(nControls, Point()); for (int i = 0; i < nControls; i++) { @@ -35,39 +44,42 @@ std::vector createControls(NumericVector x, NumericVector y) { } return controls; } -// [[Rcpp::export]] -NumericMatrix splinePath(NumericVector x, NumericVector y, int degree, std::vector knots, int detail, std::string type) { +[[cpp11::register]] +cpp11::writable::doubles_matrix<> splinePath(cpp11::doubles x, cpp11::doubles y, + int degree, cpp11::doubles knots, + int detail, cpp11::strings type) { std::vector controls = createControls(x, y); - if (type == "closed") { + std::vector knots_vec(knots.begin(), knots.end()); + if (type[0] == "closed") { controls.push_back(controls[0]); controls.push_back(controls[1]); controls.push_back(controls[2]); } - NumericMatrix res(detail, 2); - double zJump = (knots[knots.size()-1-degree] - knots[degree]); - if (type == "clamped") { + cpp11::writable::doubles_matrix<> res(detail, 2); + double zJump = (knots_vec[knots_vec.size()-1-degree] - knots_vec[degree]); + if (type[0] == "clamped") { zJump /= double(detail-1); } else { zJump /= double(detail); } - double z; - Point point; for (int i = 0; i < detail; i++) { - if (i == detail-1 && type == "clamped") { + Point point; + if (i == detail-1 && type[0] == "clamped") { point = controls[controls.size()-1]; } else { - z = knots[degree] + i * zJump; - int zInt = whichInterval(z, knots); - point = deBoor(degree, degree, zInt, z, knots, controls); + double z = knots_vec[degree] + i * zJump; + int zInt = whichInterval(z, knots_vec); + point = deBoor(degree, degree, zInt, z, knots_vec, controls); } res(i, 0) = point.x; res(i, 1) = point.y; } return res; } -// [[Rcpp::export]] -List getSplines(NumericVector x, NumericVector y, IntegerVector id, - int detail, std::string type = "clamped") { +[[cpp11::register]] +cpp11::writable::list getSplines(cpp11::doubles x, cpp11::doubles y, + cpp11::integers id, int detail, + cpp11::strings type) { std::vector nControls; std::vector pathID; nControls.push_back(1); @@ -81,42 +93,39 @@ List getSplines(NumericVector x, NumericVector y, IntegerVector id, } } int nPaths = nControls.size(); - NumericMatrix paths(nPaths * detail, 2); - IntegerVector pathsID(nPaths * detail); + cpp11::writable::doubles_matrix<> paths(nPaths * detail, 2); + cpp11::writable::integers pathsID(nPaths * detail); int controlsStart = 0; - IntegerVector controlInd; - int pathStart = 0; - IntegerVector pathInd; - IntegerVector::iterator pathIter; - int degree; - std::vector knots; - NumericMatrix path; + R_xlen_t pathStart = 0; for (int i = 0; i < nPaths; i++) { - degree = nControls[i] <= 3 ? nControls[i] - 1 : 3; - if (type == "clamped") { + cpp11::writable::doubles knots; + int degree = nControls[i] <= 3 ? nControls[i] - 1 : 3; + if (type[0] == "clamped") { knots = createKnots(nControls[i], degree); - } else if (type == "open") { + } else if (type[0] == "open") { knots = createOpenKnots(nControls[i], degree); - } else if (type == "closed") { - if (nControls[i] < 3) stop("At least 3 control points must be provided for closed b-splines"); + } else if (type[0] == "closed") { + if (nControls[i] < 3) { + cpp11::stop("At least 3 control points must be provided for closed b-splines"); + } degree = 3; knots = createOpenKnots(nControls[i] + 3, degree); } else { - stop("type must be either \"open\", \"closed\", or \"clamped\""); + cpp11::stop("type must be either \"open\", \"closed\", or \"clamped\""); } - controlInd = Range(controlsStart, controlsStart + nControls[i] - 1); - pathInd = Range(pathStart, pathStart + detail - 1); - path = splinePath(x[controlInd], y[controlInd], degree, knots, detail, type); - int j = 0; - for (pathIter = pathInd.begin(); pathIter != pathInd.end(); pathIter++) { - pathsID[*pathIter] = pathID[i]; - paths(*pathIter, 0) = path(j, 0); - paths(*pathIter, 1) = path(j, 1); - j++; + cpp11::writable::doubles x_tmp(x.begin() + controlsStart, x.begin() + controlsStart + nControls[i]); + cpp11::writable::doubles y_tmp(y.begin() + controlsStart, y.begin() + controlsStart + nControls[i]); + cpp11::doubles_matrix<> path = splinePath(x_tmp, y_tmp, degree, knots, detail, type); + for (R_xlen_t j = 0; j < path.nrow(); ++j) { + pathsID[pathStart + j] = pathID[i]; + paths(pathStart + j, 0) = path(j, 0); + paths(pathStart + j, 1) = path(j, 1); } controlsStart += nControls[i]; - pathStart += detail; + pathStart += path.nrow(); } - return List::create(Named("paths") = paths, - Named("pathID") = pathsID); + return cpp11::writable::list({ + "paths"_nm = paths, + "pathID"_nm = pathsID + }); } diff --git a/src/bezier.cpp b/src/bezier.cpp index 72d87e10..020240ed 100644 --- a/src/bezier.cpp +++ b/src/bezier.cpp @@ -1,14 +1,19 @@ -#include +#include +#include +#include +#include -using namespace Rcpp; +using namespace cpp11::literals; -double Bezier2(double t, NumericVector w) { +#include + +double Bezier2(double t, const cpp11::doubles& w) { double t2 = t * t; double mt = 1-t; double mt2 = mt * mt; return w[0]*mt2 + w[1]*2*mt*t + w[2]*t2; } -double Bezier3(double t, NumericVector w) { +double Bezier3(double t, const cpp11::doubles& w) { double t2 = t * t; double t3 = t2 * t; double mt = 1-t; @@ -16,9 +21,9 @@ double Bezier3(double t, NumericVector w) { double mt3 = mt2 * mt; return w[0]*mt3 + 3*w[1]*mt2*t + 3*w[2]*mt*t2 + w[3]*t3; } -// [[Rcpp::export]] -NumericMatrix bezierPath(NumericVector x, NumericVector y, int detail) { - NumericMatrix res(detail, 2); +[[cpp11::register]] +cpp11::writable::doubles_matrix<> bezierPath(const cpp11::doubles& x, const cpp11::doubles& y, int detail) { + cpp11::writable::doubles_matrix<> res(detail, 2); detail = detail - 1; double step = 1.0/detail; double t; @@ -35,15 +40,15 @@ NumericMatrix bezierPath(NumericVector x, NumericVector y, int detail) { res(i, 1) = Bezier3(t, y); } } else { - stop("Only support for quadratic and cubic beziers"); + cpp11::stop("Only support for quadratic and cubic beziers"); } res(detail, 0) = x[x.size() - 1]; res(detail, 1) = y[y.size() - 1]; return res; } -// [[Rcpp::export]] -List getBeziers(NumericVector x, NumericVector y, IntegerVector id, - int detail) { +[[cpp11::register]] +cpp11::writable::list getBeziers(const cpp11::doubles& x, const cpp11::doubles& y, + const cpp11::integers& id, int detail) { std::vector nControls; std::vector pathID; nControls.push_back(1); @@ -56,29 +61,25 @@ List getBeziers(NumericVector x, NumericVector y, IntegerVector id, pathID.push_back(id[i]); } } - int nPaths = nControls.size(); - NumericMatrix paths(nPaths * detail, 2); - IntegerVector pathsID(nPaths * detail); + size_t nPaths = nControls.size(); + cpp11::writable::doubles_matrix<> paths(nPaths * detail, 2); + cpp11::writable::integers pathsID(nPaths * detail); int controlsStart = 0; - IntegerVector controlInd; - int pathStart = 0; - IntegerVector pathInd; - IntegerVector::iterator pathIter; - NumericMatrix path; - for (int i = 0; i < nPaths; i++) { - controlInd = Range(controlsStart, controlsStart + nControls[i] - 1); - pathInd = Range(pathStart, pathStart + detail - 1); - path = bezierPath(x[controlInd], y[controlInd], detail); - int j = 0; - for (pathIter = pathInd.begin(); pathIter != pathInd.end(); pathIter++) { - pathsID[*pathIter] = pathID[i]; - paths(*pathIter, 0) = path(j, 0); - paths(*pathIter, 1) = path(j, 1); - j++; + R_xlen_t pathStart = 0; + for (size_t i = 0; i < nPaths; i++) { + cpp11::writable::doubles x_tmp(x.begin() + controlsStart, x.begin() + controlsStart + nControls[i]); + cpp11::writable::doubles y_tmp(y.begin() + controlsStart, y.begin() + controlsStart + nControls[i]); + cpp11::doubles_matrix<> path = bezierPath(x_tmp, y_tmp, detail); + for (R_xlen_t j = 0; j < path.nrow(); ++j) { + pathsID[pathStart + j] = pathID[i]; + paths(pathStart + j, 0) = path(j, 0); + paths(pathStart + j, 1) = path(j, 1); } controlsStart += nControls[i]; - pathStart += detail; + pathStart += path.nrow(); } - return List::create(Named("paths") = paths, - Named("pathID") = pathsID); + return cpp11::writable::list({ + "paths"_nm = paths, + "pathID"_nm = pathsID + }); } diff --git a/src/concaveman.cpp b/src/concaveman.cpp index 96230676..4fb41af5 100644 --- a/src/concaveman.cpp +++ b/src/concaveman.cpp @@ -1,10 +1,14 @@ -#include +#include +#include +#include + #include "concaveman.h" -using namespace Rcpp; +#include -// [[Rcpp::export]] -NumericMatrix concaveman_c(NumericMatrix p, IntegerVector h, double concavity, double threshold) { +[[cpp11::register]] +cpp11::writable::doubles_matrix<> concaveman_c(cpp11::doubles_matrix<> p, cpp11::integers h, + double concavity, double threshold) { typedef std::array point_type; std::vector points(p.nrow()); for (auto i = 0; i < p.nrow(); ++i) { @@ -17,7 +21,7 @@ NumericMatrix concaveman_c(NumericMatrix p, IntegerVector h, double concavity, d auto chull = concaveman(points, hull, concavity, threshold); - NumericMatrix res(chull.size(), 2); + cpp11::writable::doubles_matrix<> res(chull.size(), 2); for (size_t i = 0; i < chull.size(); ++i) { res(i, 0) = chull[i][0]; diff --git a/src/concaveman.h b/src/concaveman.h index b59e3b21..2ec0884a 100644 --- a/src/concaveman.h +++ b/src/concaveman.h @@ -7,8 +7,6 @@ #pragma once -#include "robust_predicate/geompred.hpp" - #include #include #include @@ -23,6 +21,8 @@ #include #include +#include "robust_predicate/geompred.hpp" + //#define DEBUG // uncomment to dump debug info to screen //#define DEBUG_2 // uncomment to dump second-level debug info to screen diff --git a/src/cpp11.cpp b/src/cpp11.cpp new file mode 100644 index 00000000..1f772df2 --- /dev/null +++ b/src/cpp11.cpp @@ -0,0 +1,83 @@ +// Generated by cpp11: do not edit by hand +// clang-format off + + +#include "cpp11/declarations.hpp" +#include + +// bSpline.cpp +cpp11::writable::doubles_matrix<> splinePath(cpp11::doubles x, cpp11::doubles y, int degree, cpp11::doubles knots, int detail, cpp11::strings type); +extern "C" SEXP _ggforce_splinePath(SEXP x, SEXP y, SEXP degree, SEXP knots, SEXP detail, SEXP type) { + BEGIN_CPP11 + return cpp11::as_sexp(splinePath(cpp11::as_cpp>(x), cpp11::as_cpp>(y), cpp11::as_cpp>(degree), cpp11::as_cpp>(knots), cpp11::as_cpp>(detail), cpp11::as_cpp>(type))); + END_CPP11 +} +// bSpline.cpp +cpp11::writable::list getSplines(cpp11::doubles x, cpp11::doubles y, cpp11::integers id, int detail, cpp11::strings type); +extern "C" SEXP _ggforce_getSplines(SEXP x, SEXP y, SEXP id, SEXP detail, SEXP type) { + BEGIN_CPP11 + return cpp11::as_sexp(getSplines(cpp11::as_cpp>(x), cpp11::as_cpp>(y), cpp11::as_cpp>(id), cpp11::as_cpp>(detail), cpp11::as_cpp>(type))); + END_CPP11 +} +// bezier.cpp +cpp11::writable::doubles_matrix<> bezierPath(const cpp11::doubles& x, const cpp11::doubles& y, int detail); +extern "C" SEXP _ggforce_bezierPath(SEXP x, SEXP y, SEXP detail) { + BEGIN_CPP11 + return cpp11::as_sexp(bezierPath(cpp11::as_cpp>(x), cpp11::as_cpp>(y), cpp11::as_cpp>(detail))); + END_CPP11 +} +// bezier.cpp +cpp11::writable::list getBeziers(const cpp11::doubles& x, const cpp11::doubles& y, const cpp11::integers& id, int detail); +extern "C" SEXP _ggforce_getBeziers(SEXP x, SEXP y, SEXP id, SEXP detail) { + BEGIN_CPP11 + return cpp11::as_sexp(getBeziers(cpp11::as_cpp>(x), cpp11::as_cpp>(y), cpp11::as_cpp>(id), cpp11::as_cpp>(detail))); + END_CPP11 +} +// concaveman.cpp +cpp11::writable::doubles_matrix<> concaveman_c(cpp11::doubles_matrix<> p, cpp11::integers h, double concavity, double threshold); +extern "C" SEXP _ggforce_concaveman_c(SEXP p, SEXP h, SEXP concavity, SEXP threshold) { + BEGIN_CPP11 + return cpp11::as_sexp(concaveman_c(cpp11::as_cpp>>(p), cpp11::as_cpp>(h), cpp11::as_cpp>(concavity), cpp11::as_cpp>(threshold))); + END_CPP11 +} +// ellipseEnclose.cpp +cpp11::writable::data_frame enclose_ellip_points(cpp11::doubles x, cpp11::doubles y, cpp11::integers id, double tol); +extern "C" SEXP _ggforce_enclose_ellip_points(SEXP x, SEXP y, SEXP id, SEXP tol) { + BEGIN_CPP11 + return cpp11::as_sexp(enclose_ellip_points(cpp11::as_cpp>(x), cpp11::as_cpp>(y), cpp11::as_cpp>(id), cpp11::as_cpp>(tol))); + END_CPP11 +} +// enclose.cpp +cpp11::writable::data_frame enclose_points(cpp11::doubles x, cpp11::doubles y, cpp11::integers id); +extern "C" SEXP _ggforce_enclose_points(SEXP x, SEXP y, SEXP id) { + BEGIN_CPP11 + return cpp11::as_sexp(enclose_points(cpp11::as_cpp>(x), cpp11::as_cpp>(y), cpp11::as_cpp>(id))); + END_CPP11 +} +// pointPath.cpp +cpp11::writable::list points_to_path(cpp11::doubles_matrix<> pos, cpp11::list_of< cpp11::doubles_matrix<> > path, bool close); +extern "C" SEXP _ggforce_points_to_path(SEXP pos, SEXP path, SEXP close) { + BEGIN_CPP11 + return cpp11::as_sexp(points_to_path(cpp11::as_cpp>>(pos), cpp11::as_cpp >>>(path), cpp11::as_cpp>(close))); + END_CPP11 +} + +extern "C" { +static const R_CallMethodDef CallEntries[] = { + {"_ggforce_bezierPath", (DL_FUNC) &_ggforce_bezierPath, 3}, + {"_ggforce_concaveman_c", (DL_FUNC) &_ggforce_concaveman_c, 4}, + {"_ggforce_enclose_ellip_points", (DL_FUNC) &_ggforce_enclose_ellip_points, 4}, + {"_ggforce_enclose_points", (DL_FUNC) &_ggforce_enclose_points, 3}, + {"_ggforce_getBeziers", (DL_FUNC) &_ggforce_getBeziers, 4}, + {"_ggforce_getSplines", (DL_FUNC) &_ggforce_getSplines, 5}, + {"_ggforce_points_to_path", (DL_FUNC) &_ggforce_points_to_path, 3}, + {"_ggforce_splinePath", (DL_FUNC) &_ggforce_splinePath, 6}, + {NULL, NULL, 0} +}; +} + +extern "C" attribute_visible void R_init_ggforce(DllInfo* dll){ + R_registerRoutines(dll, NULL, CallEntries, NULL, NULL); + R_useDynamicSymbols(dll, FALSE); + R_forceSymbols(dll, TRUE); +} diff --git a/src/ellipseEnclose.cpp b/src/ellipseEnclose.cpp index 7dd1c1fa..52a9e590 100644 --- a/src/ellipseEnclose.cpp +++ b/src/ellipseEnclose.cpp @@ -1,8 +1,15 @@ -#include -#include +#include +#include +#include +#include +#include +#include +#include -// [[Rcpp::depends(RcppEigen)]] -using namespace Rcpp; +using namespace cpp11::literals; + +#include +#include struct Ellipse { double x; @@ -12,30 +19,30 @@ struct Ellipse { double rad; }; -bool points_on_line(Eigen::MatrixXd points, Ellipse &enc) { +bool points_on_line(const cpp11::doubles_matrix<>& points, Ellipse &enc) { double xmin, ymin, xmax, ymax; - int n = points.rows(); + R_xlen_t n = points.nrow(); if (n == 1) { - enc.x = points.coeff(0, 0); - enc.y = points.coeff(0, 1); + enc.x = points(0, 0); + enc.y = points(0, 1); enc.a = 0; enc.b = 0; enc.rad = 0; return true; } if (n == 2) { - xmin = points.coeff(0, 0); - xmax = points.coeff(1, 0); - ymin = points.coeff(0, 1); - ymax = points.coeff(1, 1); + xmin = points(0, 0); + xmax = points(1, 0); + ymin = points(0, 1); + ymax = points(1, 1); } else { - double x0 = xmin = xmax = points.coeff(0, 0); - double y0 = ymin = ymax = points.coeff(0, 1); - double xdiff = points.coeff(1, 0) - x0; + double x0 = xmin = xmax = points(0, 0); + double y0 = ymin = ymax = points(0, 1); + double xdiff = points(1, 0) - x0; bool vert = xdiff == 0; - double slope; + double slope = std::numeric_limits::infinity(); if (!vert) { - slope = (points.coeff(1, 1) - y0) / xdiff; + slope = (points(1, 1) - y0) / xdiff; } for (int i = 2; i < n; i++) { xdiff = points(i, 0) - x0; @@ -44,7 +51,7 @@ bool points_on_line(Eigen::MatrixXd points, Ellipse &enc) { ymax = std::max(ymax, points(i, 1)); continue; } - if (slope == (points.coeff(i, 1) - y0) / xdiff) { + if (slope == (points(i, 1) - y0) / xdiff) { if (points(i, 0) < xmin) { xmin = points(i, 0); ymin = points(i, 1); @@ -75,81 +82,222 @@ bool points_on_line(Eigen::MatrixXd points, Ellipse &enc) { return true; } -Ellipse khachiyan(Eigen::MatrixXd points, double tol) { +cpp11::writable::doubles_matrix<> transpose(const cpp11::doubles_matrix<>& x) { + cpp11::writable::doubles_matrix<> trans(x.ncol(), x.nrow()); + for (R_xlen_t j = 0; j < x.ncol(); ++j) for (R_xlen_t i = 0; i < x.nrow(); ++i) { + trans(j, i) = x(i, j); + } + return trans; +} + +cpp11::writable::doubles_matrix<> make_q(const cpp11::doubles_matrix<>& p) { + cpp11::writable::doubles_matrix<> Q(p.nrow() + 1, p.ncol()); + for (R_xlen_t j = 0; j < p.ncol(); ++j) for (R_xlen_t i = 0; i < p.nrow(); ++i) { + Q(i, j) = p(i, j); + } + for (R_xlen_t j = 0; j < p.ncol(); ++j) Q(p.nrow(), j) = 1.0; + + return Q; +} + +// Q * diagonal(u) * Qadj +cpp11::writable::doubles_matrix<> solve_x(const cpp11::doubles_matrix<>& Q, const cpp11::doubles& u, const cpp11::doubles_matrix<>& Qadj) { + cpp11::writable::doubles_matrix<> res(Q.nrow(), Qadj.ncol()); + for (R_xlen_t i = 0; i < Q.nrow(); ++i) for (R_xlen_t j = 0; j < Qadj.ncol(); ++j) { + double cell = 0.0; + for (R_xlen_t k = 0; k < u.size(); ++k) { + cell += Q(i, k) * u[k] * Qadj(k, j); + } + res(i, j) = cell; + } + return res; +} + +// diagonal(Qadj * inverse(X) * Q) +cpp11::writable::doubles solve_m(const cpp11::doubles_matrix<>& Qadj, const cpp11::doubles_matrix<>& X, const cpp11::doubles_matrix<>& Q) { + static auto solve = cpp11::package("base")["solve"]; + cpp11::doubles_matrix<> Xinv = cpp11::as_cpp< cpp11::doubles_matrix<> >(solve(X)); + + cpp11::writable::doubles res; + + for (R_xlen_t i = 0; i < Qadj.nrow(); ++i) { + double sum = 0.0; + for (R_xlen_t k = 0; k < Xinv.ncol(); ++k) { + double cell = 0.0; + for (R_xlen_t j = 0; j < Qadj.ncol(); ++j) { + cell += Qadj(i, j) * Xinv(j, k); + } + sum += cell * Q(k, i); + } + res.push_back(sum); + } + + return res; +} + +cpp11::writable::doubles_matrix<> prod_with_diag(const cpp11::doubles_matrix<>& A, const cpp11::doubles& b) { + if (A.ncol() != b.size()) { + cpp11::stop("A must have the same number of columns as there are elements in b"); + } + cpp11::writable::doubles_matrix<> res(A.nrow(), A.ncol()); + for (R_xlen_t i = 0; i < A.nrow(); ++i) for (R_xlen_t j = 0; j < A.ncol(); ++j) { + res(i, j) = A(i, j) * b[j]; + } + return res; +} + +cpp11::writable::doubles_matrix<> prod(const cpp11::doubles_matrix<>& A, const cpp11::doubles_matrix<>& B) { + if (A.ncol() != B.nrow()) { + cpp11::stop("A must have the same number of columns as there are rows in B"); + } + cpp11::writable::doubles_matrix<> res(A.nrow(), B.ncol()); + for (R_xlen_t i = 0; i < A.nrow(); ++i) for (R_xlen_t j = 0; j < B.ncol(); ++j) { + double cell = 0.0; + for (R_xlen_t k = 0; k < A.ncol(); ++k) { + cell += A(i, k) * B(k, j); + } + res(i, j) = cell; + } + return res; +} +cpp11::writable::doubles_matrix<> prod_with_vec(const cpp11::doubles_matrix<>& A, const cpp11::doubles& b) { + if (A.ncol() != b.size()) { + cpp11::stop("A must have the same number of columns as there are elements in b"); + } + cpp11::writable::doubles_matrix<> res(A.nrow(), 1); + for (R_xlen_t i = 0; i < A.nrow(); ++i) { + double cell = 0.0; + for (R_xlen_t k = 0; k < A.ncol(); ++k) { + cell += A(i, k) * b[k]; + } + res(i, 0) = cell; + } + return res; +} +// ------------------------------------- fourth ------------------------------------------------ +// ----------- third ----------------- +// --------------- first -------------------- - second - +// (1.0/d) * (points * u.asDiagonal() * points.adjoint() - (points * u)*(points * u).adjoint() ).inverse() +cpp11::writable::doubles_matrix<> solve_a(const cpp11::doubles_matrix<>& points, const cpp11::doubles& u, const cpp11::doubles_matrix<>& second) { + static auto solve = cpp11::package("base")["solve"]; + cpp11::writable::doubles_matrix<> first = prod(prod_with_diag(points, u), transpose(points)); + cpp11::doubles_matrix<> third = prod(second, transpose(second)); + + for (R_xlen_t i = 0; i < first.nrow(); ++i) for (R_xlen_t j = 0; j < first.ncol(); ++j) { + first(i, j) -= third(i, j); + } + + cpp11::writable::doubles_matrix<> fourth = cpp11::as_cpp< cpp11::writable::doubles_matrix<> >(solve(first)); + + for (R_xlen_t i = 0; i < fourth.nrow(); ++i) for (R_xlen_t j = 0; j < fourth.ncol(); ++j) { + fourth(i, j) *= 1.0 / points.nrow(); + } + + return fourth; +} + +Ellipse solve_enc(const cpp11::doubles_matrix<>& A) { + static auto svd = cpp11::package("base")["svd"]; + Ellipse enc; + cpp11::list svd_a = cpp11::as_cpp(svd(A)); + + cpp11::doubles d(svd_a["d"]); + cpp11::doubles_matrix<> v(svd_a["v"]); + + enc.a = 1.0/std::sqrt(d[1]); + enc.b = 1.0/std::sqrt(d[0]); + if (v(0, 1) == v(1, 0)) { + enc.rad = std::asin(v(1, 1)); + } else if (v(0, 1) < v(1, 0)) { + enc.rad = std::asin(-v(0, 0)); + } else { + enc.rad = std::asin(v(0, 0)); + } + return enc; +} + +Ellipse khachiyan(const cpp11::doubles_matrix<>& points, double tol) { Ellipse enc; if (points_on_line(points, enc)) { return enc; } - points.adjointInPlace(); - int N = points.cols(); - int d = points.rows(); - Eigen::MatrixXd Q; - Q = points; - Q.conservativeResize(d + 1, Eigen::NoChange); - Q.row(d).setOnes(); - Eigen::MatrixXd Qadj = Q.adjoint(); - double error = 1; - double max, step; - Eigen::MatrixXd X, A, V; - Eigen::VectorXd M, u_tmp, c; - Eigen::VectorXd::Index max_i; - Eigen::VectorXd u(N); - u.fill(1/double(N)); + cpp11::doubles_matrix<> pointsAdj = transpose(points); + + R_xlen_t N = pointsAdj.ncol(); + R_xlen_t d = pointsAdj.nrow(); + cpp11::doubles_matrix<> Q = make_q(pointsAdj); + cpp11::doubles_matrix<> Qadj = transpose(Q); + + cpp11::writable::doubles u(N); + std::fill(u.begin(), u.end(), 1/double(N)); + cpp11::writable::doubles u_tmp(N); + double error = 1; while (error > tol) { - X = Q * u.asDiagonal() * Qadj; - M = (Qadj * X.inverse() * Q).diagonal(); - max = M.maxCoeff(&max_i); - step = (max - d - 1)/((d + 1)*(max - 1)); - u_tmp = (1 - step)*u; - u_tmp[max_i] = u_tmp[max_i] + step; - error = (u_tmp - u).norm(); - u = u_tmp; - } - A = (1.0/d) * (points * u.asDiagonal() * points.adjoint() - (points * u)*(points * u).adjoint() ).inverse(); - c = points * u; - enc.x = c[0]; - enc.y = c[1]; - - Eigen::JacobiSVD svd_A(A, Eigen::ComputeThinV); - enc.a = 1.0/std::sqrt(float(svd_A.singularValues()[1])); - enc.b = 1.0/std::sqrt(float(svd_A.singularValues()[0])); - V = svd_A.matrixV(); - if (V(0, 1) == V(1, 0)) { - enc.rad = std::asin(float(V(1, 1))); - } else if (V(0, 1) < V(1, 0)) { - enc.rad = std::asin(-float(V(0, 0))); - } else { - enc.rad = std::asin(float(V(0, 0))); + cpp11::doubles_matrix<> X = solve_x(Q, u, Qadj); + cpp11::doubles M = solve_m(Qadj, X, Q); + R_xlen_t max_i = 0; + double max = M[0]; + for (R_xlen_t i = 1; i < M.size(); ++i) { + if (M[i] > max) { + max_i = i; + max = M[i]; + } + } + double step = (max - d - 1)/((d + 1)*(max - 1)); + cpp11::writable::doubles u_tmp; + u_tmp.reserve(N); + for (R_xlen_t i = 0; i < u.size(); ++i) { + u_tmp.push_back(u[i] * (1 - step)); + } + u_tmp[max_i] += step; + + error = 0.0; + for (R_xlen_t i = 0; i < u.size(); ++i) { + double du = u_tmp[i] - u[i]; + error += du * du; + u[i] = double(u_tmp[i]); + } + error = std::sqrt(error); } + cpp11::doubles_matrix<> c = prod_with_vec(pointsAdj, u); + cpp11::doubles_matrix<> A = solve_a(pointsAdj, u, c); + enc = solve_enc(A); + enc.x = c(0, 0); + enc.y = c(1, 0); + return enc; } -// [[Rcpp::export]] -DataFrame enclose_ellip_points(NumericVector x, NumericVector y, IntegerVector id, double tol) { +[[cpp11::register]] +cpp11::writable::data_frame enclose_ellip_points(cpp11::doubles x, cpp11::doubles y, cpp11::integers id, double tol) { if (x.size() != y.size() || x.size() != id.size()) { - stop("x, y, and id must have same dimensions"); + cpp11::stop("x, y, and id must have same dimensions"); } - const Eigen::Map X(as< Eigen::Map >(x)); - const Eigen::Map Y(as< Eigen::Map >(y)); - std::vector< double > x0, y0, a, b, rad; std::vector< int > splits; splits.push_back(0); int currentId = id[0]; - int i, size; - for (i = 0; i < id.size(); ++i) { + for (R_xlen_t i = 0; i < id.size(); ++i) { if (id[i] != currentId) { currentId = id[i]; splits.push_back(i); } } splits.push_back(id.size()); - for (i = 0; i < splits.size() - 1; ++i) { - size = splits[i+1] - splits[i]; - Eigen::MatrixXd points(size, 2); - points.col(0) = X.segment(splits[i], size); - points.col(1) = Y.segment(splits[i], size); + + cpp11::writable::doubles x0; + cpp11::writable::doubles y0; + cpp11::writable::doubles a; + cpp11::writable::doubles b; + cpp11::writable::doubles rad; + for (size_t i = 0; i < splits.size() - 1; ++i) { + int size = splits[i+1] - splits[i]; + cpp11::writable::doubles_matrix<> points(size, 2); + for (R_xlen_t j = 0; j < points.nrow(); ++j) { + points(j, 0) = x[splits[i] + j]; + points(j, 1) = y[splits[i] + j]; + } Ellipse center = khachiyan(points, tol); x0.push_back(center.x); y0.push_back(center.y); @@ -157,11 +305,12 @@ DataFrame enclose_ellip_points(NumericVector x, NumericVector y, IntegerVector i b.push_back(center.b); rad.push_back(center.rad); } - return DataFrame::create( - Named("x0") = wrap(x0), - Named("y0") = wrap(y0), - Named("a") = wrap(a), - Named("b") = wrap(b), - Named("angle") = wrap(rad) - ); + + return cpp11::writable::data_frame({ + "x0"_nm = x0, + "y0"_nm = y0, + "a"_nm = a, + "b"_nm = b, + "angle"_nm = rad + }); } diff --git a/src/enclose.cpp b/src/enclose.cpp index 4abf10d5..d9391992 100644 --- a/src/enclose.cpp +++ b/src/enclose.cpp @@ -1,9 +1,11 @@ -#include -using namespace Rcpp; +#include +#include +#include + +using namespace cpp11::literals; + +#include -inline int randWrapper(const int n) { - return std::floor(float(unif_rand()*n)); -} struct Circle { double x; double y; @@ -13,19 +15,19 @@ struct Point { double x; double y; }; -bool equalPoints(Point &p1, Point &p2) { - return std::abs(float(p2.x - p1.x)) < 1e-9 && std::abs(float(p2.y - p1.y)) < 1e-9; +bool equalPoints(const Point &p1, const Point &p2) { + return std::abs(p2.x - p1.x) < 1e-9 && std::abs(p2.y - p1.y) < 1e-9; } -bool perpendicularPoints(Point &p1, Point &p2) { - return std::abs(float(p2.x - p1.x)) < 1e-9 || std::abs(float(p2.y - p1.y)) < 1e-9; +bool perpendicularPoints(const Point &p1, const Point &p2) { + return std::abs(p2.x - p1.x) < 1e-9 || std::abs(p2.y - p1.y) < 1e-9; } -bool horizontalPoints(Point &p1, Point &p2) { - return std::abs(float(p2.y - p1.y)) < 1e-9; +bool horizontalPoints(const Point &p1, const Point &p2) { + return std::abs(p2.y - p1.y) < 1e-9; } -bool verticalPoints(Point &p1, Point &p2) { - return std::abs(float(p2.x - p1.x)) < 1e-9; +bool verticalPoints(const Point &p1, const Point &p2) { + return std::abs(p2.x - p1.x) < 1e-9; } -Circle circleByPoints(Point &p1, Point &p2, Point &p3) { +Circle circleByPoints(const Point &p1, const Point &p2, const Point &p3) { Circle results; double X1, Y1, X2, Y2, A1, A2; X1 = p2.x - p1.x; @@ -34,29 +36,31 @@ Circle circleByPoints(Point &p1, Point &p2, Point &p3) { Y2 = p3.y - p2.y; A1 = Y1/X1; A2 = Y2/X2; - if (std::abs(float(A2 - A1)) < 1e-9) stop("Error in circleByPoints: The 3 points are colinear"); - results.x = ( A1*A2*(p1.y - p3.y) + A2*(p1.x + p2.x) - A1*(p2.x+p3.x) )/( 2* (A2 - A1) ); - results.y = -1*( results.x - (p1.x + p2.x)/2 )/A1 + (p1.y + p2.y)/2; + if (std::abs(A2 - A1) < 1e-9) { + cpp11::stop("Error in circleByPoints: The 3 points are colinear"); + } + results.x = ( A1*A2*(p1.y - p3.y) + A2*(p1.x + p2.x) - A1*(p2.x+p3.x) )/( 2.0 * (A2 - A1) ); + results.y = -1.0*( results.x - (p1.x + p2.x)/2.0 )/A1 + (p1.y + p2.y)/2.0; return results; } -Circle encloseOne(Point &p1) { +Circle encloseOne(const Point &p1) { Circle results; results.x = p1.x; results.y = p1.y; - results.r = 0; + results.r = 0.0; return results; } -Circle encloseTwo(Point &p1, Point &p2) { +Circle encloseTwo(const Point &p1, const Point &p2) { if (equalPoints(p1, p2)) return encloseOne(p1); Circle results; double dx = p2.x - p1.x; double dy = p2.y - p1.y; - results.x = p1.x + dx/2; - results.y = p1.y + dy/2; - results.r = std::sqrt(float(dx*dx + dy*dy))/2; + results.x = p1.x + dx/2.0; + results.y = p1.y + dy/2.0; + results.r = std::sqrt(dx*dx + dy*dy)/2.0; return results; } -Circle encloseThree(Point &p1, Point &p2, Point &p3) { +Circle encloseThree(const Point &p1, const Point &p2, const Point &p3) { if (equalPoints(p1, p2)) return encloseTwo(p1, p3); if (equalPoints(p1, p3)) return encloseTwo(p1, p2); if (equalPoints(p2, p3)) return encloseTwo(p1, p2); @@ -64,8 +68,8 @@ Circle encloseThree(Point &p1, Point &p2, Point &p3) { bool perp23 = perpendicularPoints(p2, p3); bool perp13 = perpendicularPoints(p1, p3); Circle results; - if (perp12 + perp23 + perp13 == 3) { - stop("Error in encloseThree: The 3 points are colinear"); + if (perp12 && perp23 && perp13) { + cpp11::stop("Error in encloseThree: The 3 points are colinear"); } else if (perp12 + perp23 + perp13 == 2) { if (perp12) { if (horizontalPoints(p1, p2)) { @@ -108,7 +112,7 @@ Circle encloseDefault(std::vector points) { case 1: return encloseOne(points[0]); case 2: return encloseTwo(points[0], points[1]); case 3: return encloseThree(points[0], points[1], points[2]); - default: stop("Error in encloseDefault - expecting less than 4 points"); + default: cpp11::stop("Error in encloseDefault - expecting less than 4 points"); } } bool inCircle(Circle c, Point p) { @@ -156,10 +160,9 @@ std::vector extendPerimeter(std::vector perimeter, Point p) { } } } - stop("Error in extendPerimeter: Could not enclose points"); + cpp11::stop("Error in extendPerimeter: Could not enclose points"); } Circle enclosePoints(std::vector points) { - //std::random_shuffle(points.begin(), points.end(), randWrapper); std::vector::iterator it = points.begin(); Circle center = {0.0, 0.0, 0.0}; std::vector perimeter; @@ -175,18 +178,19 @@ Circle enclosePoints(std::vector points) { return center; } -// [[Rcpp::export]] -DataFrame enclose_points(NumericVector x, NumericVector y, IntegerVector id) { +[[cpp11::register]] +cpp11::writable::data_frame enclose_points(cpp11::doubles x, cpp11::doubles y, cpp11::integers id) { if (x.size() != y.size() || x.size() != id.size()) { - stop("x, y, and id must have same dimensions"); + cpp11::stop("x, y, and id must have same dimensions"); } - std::vector< double > x0, y0, r; + cpp11::writable::doubles x0; + cpp11::writable::doubles y0; + cpp11::writable::doubles r; std::vector< std::vector > all_points; std::vector points; all_points.push_back(points); int currentId = id[0]; - int i; - for (i = 0; i < id.size(); ++i) { + for (R_xlen_t i = 0; i < id.size(); ++i) { Point p_tmp = {x[i], y[i]}; if (id[i] != currentId) { currentId = id[i]; @@ -195,15 +199,15 @@ DataFrame enclose_points(NumericVector x, NumericVector y, IntegerVector id) { } all_points.back().push_back(p_tmp); } - for (i = 0; i < all_points.size(); ++i) { + for (size_t i = 0; i < all_points.size(); ++i) { Circle center = enclosePoints(all_points[i]); x0.push_back(center.x); y0.push_back(center.y); r.push_back(center.r); } - return DataFrame::create( - Named("x0") = wrap(x0), - Named("y0") = wrap(y0), - Named("r") = wrap(r) - ); + return cpp11::writable::data_frame({ + "x0"_nm = x0, + "y0"_nm = y0, + "r"_nm = r + }); } diff --git a/src/pointPath.cpp b/src/pointPath.cpp index fbcbc75f..dd6f7afc 100644 --- a/src/pointPath.cpp +++ b/src/pointPath.cpp @@ -1,7 +1,12 @@ -#include -#include +#include +#include +#include +#include + +using namespace cpp11::literals; -using namespace Rcpp; +#include +#include double distSquared(std::pair p, std::pair p1) { double x = p1.first - p.first; @@ -21,22 +26,18 @@ std::pair projection(std::pair a, std::pair path, std::vector &res, bool closed_poly) { - int i, j, k; - double dist, shortest_dist = -1; - std::pair point, a, b, close, closest; - point.first = x; - point.second = y; - for (i = 0; i < path.size(); ++i) { - for (j = 0; j < path[i].nrow(); ++j) { +void dist_to_path(double x, double y, const cpp11::list_of< cpp11::doubles_matrix<> >& path, std::vector &res, bool closed_poly) { + double shortest_dist = -1; + std::pair closest; + std::pair point(x, y); + for (R_xlen_t i = 0; i < path.size(); ++i) { + for (R_xlen_t j = 0; j < path[i].nrow(); ++j) { if (j == path[i].nrow() && !closed_poly) break; - a.first = path[i](j, 0); - a.second = path[i](j, 1); - k = j == path[i].nrow() - 1 ? 0 : j + 1; - b.first = path[i](k, 0); - b.second = path[i](k, 1); - close = projection(a, b, point, true); - dist = std::sqrt(distSquared(point, close)); + std::pair a(path[i](j, 0), path[i](j, 1)); + R_xlen_t k = j == path[i].nrow() - 1 ? 0 : j + 1; + std::pair b(path[i](k, 0), path[i](k, 1)); + std::pair close = projection(a, b, point, true); + double dist = std::sqrt(distSquared(point, close)); if (shortest_dist < 0 || dist < shortest_dist) { shortest_dist = dist; closest = close; @@ -49,19 +50,19 @@ void dist_to_path(double x, double y, ListOf path, std::vector path, bool close) { +[[cpp11::register]] +cpp11::writable::list points_to_path(cpp11::doubles_matrix<> pos, cpp11::list_of< cpp11::doubles_matrix<> > path, bool close) { std::vector res_container; - NumericMatrix proj(pos.nrow(), 2); - NumericVector dist(pos.nrow()); - for (int i = 0; i < pos.nrow(); ++i) { + cpp11::writable::doubles_matrix<> proj(pos.nrow(), 2); + cpp11::writable::doubles dist(pos.nrow()); + for (R_xlen_t i = 0; i < pos.nrow(); ++i) { dist_to_path(pos(i, 0), pos(i, 1), path, res_container, close); proj(i, 0) = res_container[0]; proj(i, 1) = res_container[1]; dist[i] = res_container[2]; } - return List::create( - _["projection"] = proj, - _["distance"] = dist - ); + return cpp11::writable::list({ + "projection"_nm = proj, + "distance"_nm = dist + }); }