Skip to content

Instantly share code, notes, and snippets.

@tim-salabim
Last active October 28, 2017 11:10
Show Gist options
  • Save tim-salabim/561afd611fea2acdd14b4b6acce1bc8e to your computer and use it in GitHub Desktop.
Save tim-salabim/561afd611fea2acdd14b4b6acce1bc8e to your computer and use it in GitHub Desktop.
perp points in the making
#include <Rcpp.h>
#include <math.h>
using namespace Rcpp;
// helper functions to convert degrees to radians and vice versa =============
// [[Rcpp::export]]
double deg(double radians) {
return 180 * radians / M_PI;
}
// [[Rcpp::export]]
double rad(double degrees) {
return degrees * M_PI / 180;
}
// ===========================================================================
// convert delta x and delta y to direction and length (or u & v to wd & ws)
// ===========================================================================
// scalar operation to convert delta x and delta y to direction
// [[Rcpp::export]]
double xy2dir(double x, double y) {
double w = deg(atan2(y, x));
if (w <= 0) {
w = w + 360;
}
double dir = 90 - w;
if (w >= 90) {
dir = dir + 360;
}
return dir;
}
// vector version of xy2dir
// [[Rcpp::export]]
NumericVector calcDir(NumericVector x, NumericVector y) {
int x_len = x.size();
// int y_len = y.size();
NumericVector res(x_len); // init result vector
for(int i = 0; i < x_len; i++) {
res[i] = xy2dir(x[i], y[i]);
}
return res;
}
// scalar operation to convert delta x and delta y to length
// [[Rcpp::export]]
double xy2len(double x, double y) {
return sqrt(x*x + y*y);
}
// vector version of xy2len
// [[Rcpp::export]]
NumericVector calcLen(NumericVector x, NumericVector y) {
int x_len = x.size();
// int y_len = y.size();
NumericVector res(x_len); // init result vector
for(int i = 0; i < x_len; i++) {
res[i] = xy2len(x[i], y[i]);
}
return res;
}
// ===========================================================================
// convert direction and length to delta x and delta y (or wd & ws to u & v)
// ===========================================================================
// scalar opration to convert direction and length to delta x
// [[Rcpp::export]]
double dirlen2x(double dir, double len) {
return len * sin(rad(dir));
}
// vector version of dirlen2x
// [[Rcpp::export]]
NumericVector calcDeltaX(NumericVector dir, NumericVector len) {
int res_len = dir.size();
NumericVector res(res_len);
for (int i = 0; i < res_len; i++) {
res[i] = dirlen2x(dir[i], len[i]);
}
return(res);
}
// scalar opration to convert direction and length to delta x
// [[Rcpp::export]]
double dirlen2y(double dir, double len) {
return len * cos(rad(dir));
}
// vector version of dirlen2y
// [[Rcpp::export]]
NumericVector calcDeltaY(NumericVector dir, NumericVector len) {
int res_len = dir.size();
NumericVector res(res_len);
for (int i = 0; i < res_len; i++) {
res[i] = dirlen2y(dir[i], len[i]);
}
return(res);
}
// ===========================================================================
// perp point calculation
// ===========================================================================
// delta x on the right side (to be added to midpoint)
// [[Rcpp::export]]
double perpRX(double dir,
double len) {
double right = dir + 90;
return dirlen2x(right, len);
}
// delta y on the right side (to be added to midpoint)
// [[Rcpp::export]]
double perpRY(double dir,
double len) {
double left = dir + 90;
return dirlen2y(left, len);
}
// delta x on the left side (to be added to midpoint)
// [[Rcpp::export]]
double perpLX(double dir,
double len) {
double right = dir + 270;
return dirlen2x(right, len);
}
// delta y on the left side (to be added to midpoint)
// [[Rcpp::export]]
double perpLY(double dir,
double len) {
double left = dir + 270;
return dirlen2y(left, len);
}
// calculate perp points
// [[Rcpp::export]]
NumericMatrix calcPerpPointR(NumericVector x_offset, NumericVector y_offset, double dist) {
int nrow = x_offset.size();
NumericMatrix mat(nrow, 2);
for (int i = 0; i < nrow; i++) {
double dir = xy2dir(x_offset[i], y_offset[i]);
double len = xy2len(x_offset[i], y_offset[i]);
double x1_offset = dirlen2x(dir, len * 1) + perpRX(dir, dist);
double y1_offset = dirlen2y(dir, len * 1) + perpRY(dir, dist);
mat(i, 0) = x1_offset;
mat(i, 1) = y1_offset;
}
return mat;
}
// calculate perp points
// [[Rcpp::export]]
NumericMatrix calcPerpPointL(NumericVector x_offset, NumericVector y_offset, double dist) {
int nrow = x_offset.size();
NumericMatrix mat(nrow, 2);
for (int i = 0; i < nrow; i++) {
double dir = xy2dir(x_offset[i], y_offset[i]);
double len = xy2len(x_offset[i], y_offset[i]);
double x1_offset = dirlen2x(dir, len * 1) + perpLX(dir, dist);
double y1_offset = dirlen2y(dir, len * 1) + perpLY(dir, dist);
mat(i, 0) = x1_offset;
mat(i, 1) = y1_offset;
}
return mat;
}
// rcpp diff version (from http://gallery.rcpp.org/articles/sugar-diff/)
// [[Rcpp::export]]
NumericVector diff_cpp(NumericVector x){
return diff(x);
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment