1.9

DESCRIPTION

@@ -8,10 +8,7 @@ Authors@R: c(
     comment = c(ORCID = "0000-0002-1057-3721")),
     person("Pascal", "Mettes", 
              role = "ctb", 
-             comment = "Author of the initial C++ implementation of the SLIC Superpixel algorithm for image data"),
-    person("Charles", "Jekel", 
-             role = "ctb", 
-             comment = "Author of underlying C++ code for dtw"))
+             comment = "Author of the initial C++ implementation of the SLIC Superpixel algorithm for image data"))
 Description: Creates superpixels based on input spatial data. 
   This package works on spatial data with one variable (e.g., continuous raster), many variables (e.g., RGB rasters), and spatial patterns (e.g., areas in categorical rasters).
   It is based on the SLIC algorithm (Achanta et al. (2012) <doi:10.1109/TPAMI.2012.120>), and readapts it to work with arbitrary dissimilarity measures. 

src/distances.cpp

@@ -1,5 +1,4 @@
 #include "distances.h"
-#include "dtw/include/DTW.hpp"
 // using namespace cpp11::literals; // so we can use ""_nm syntax
 
 double get_vals_dist(const std::vector<double>& values1, const std::vector<double>& values2,
@@ -9,7 +8,7 @@ double get_vals_dist(const std::vector<double>& values1, const std::vector<doubl
   } else if (dist_name == "jsd"){
     return jensen_shannon(values1, values2);
   } else if (dist_name == "dtw"){
-    return dtw3(values1, values2);
+    return dtw1d(values1, values2);
   } else if (dist_name == "dtw2d"){
     return dtw2d(values1, values2);
   } else if (dist_name != ""){
@@ -43,50 +42,65 @@ double manhattan(const std::vector<double>& values1, const std::vector<double>&
   return dist;
 }
 
-double dtw3(const std::vector<double>& values1, const std::vector<double>& values2){
-  double p = 2;  // the p-norm to use; 2.0 = euclidean, 1.0 = manhattan
+double dtw1d(const std::vector<double>& values1, const std::vector<double>& values2){
   int len1 = values1.size();
+  int len2 = values2.size();
+  if (len1 == 0 || len2 == 0) {
+    return NAN;
+  }
 
-  std::vector<std::vector<double> > a;
-  std::vector<std::vector<double> > b;
+  std::vector<double> prev(len2);
+  std::vector<double> curr(len2);
 
-  a.reserve(len1);
-  b.reserve(len1);
+  prev[0] = std::fabs(values1[0] - values2[0]);
+  for (int j = 1; j < len2; j++) {
+    prev[j] = prev[j - 1] + std::fabs(values1[0] - values2[j]);
+  }
 
-  for(int i = 0; i < len1; i++){
-    std::vector<double> pair1(2); std::vector<double> pair2(2);
-    pair1[0] = i;
-    pair1[1] = values1[i];
-    pair2[0] = i;
-    pair2[1] = values2[i];
-    a.push_back(pair1); b.push_back(pair2);
+  for (int i = 1; i < len1; i++) {
+    curr[0] = prev[0] + std::fabs(values1[i] - values2[0]);
+    for (int j = 1; j < len2; j++) {
+      double cost = std::fabs(values1[i] - values2[j]);
+      double best = std::fmin(std::fmin(prev[j], curr[j - 1]), prev[j - 1]);
+      curr[j] = cost + best;
+    }
+    prev.swap(curr);
   }
 
-  double scost = DTW::dtw_distance_only(a, b, p);
-  return(scost);
+  return prev[len2 - 1];
 }
 
 double dtw2d(const std::vector<double>& values1, const std::vector<double>& values2){
-  double p = 2;  // the p-norm to use; 2.0 = euclidean, 1.0 = manhattan
   int len1 = values1.size() / 2;
+  int len2 = values2.size() / 2;
+  if (len1 == 0 || len2 == 0) {
+    return NAN;
+  }
 
-  std::vector<std::vector<double> > a;
-  std::vector<std::vector<double> > b;
+  std::vector<double> prev(len2);
+  std::vector<double> curr(len2);
 
-  a.reserve(len1);
-  b.reserve(len1);
+  auto cost = [&](int i, int j) {
+    double dx = values1[i + len1] - values2[j + len2];
+    double dy = values1[i] - values2[j];
+    return std::sqrt(dx * dx + dy * dy);
+  };
 
-  for(int i = 0; i < len1; i++){
-    std::vector<double> pair1(2); std::vector<double> pair2(2);
-    pair1[0] = values1[i + (len1 - 1)];
-    pair1[1] = values1[i];
-    pair2[0] = values2[i + (len1 - 1)];
-    pair2[1] = values2[i];
-    a.push_back(pair1); b.push_back(pair2);
+  prev[0] = cost(0, 0);
+  for (int j = 1; j < len2; j++) {
+    prev[j] = prev[j - 1] + cost(0, j);
+  }
+
+  for (int i = 1; i < len1; i++) {
+    curr[0] = prev[0] + cost(i, 0);
+    for (int j = 1; j < len2; j++) {
+      double best = std::fmin(std::fmin(prev[j], curr[j - 1]), prev[j - 1]);
+      curr[j] = cost(i, j) + best;
+    }
+    prev.swap(curr);
   }
 
-  double scost = DTW::dtw_distance_only(a, b, p);
-  return(scost);
+  return prev[len2 - 1];
 }
 
 

src/distances.h

@@ -11,7 +11,7 @@ double euclidean(const std::vector<double>& values1, const std::vector<double>&
 double manhattan(const std::vector<double>& values1, const std::vector<double>& values2);
 double jensen_shannon(const std::vector<double>& values1, const std::vector<double>& values2);
 double custom_log2(double x);
-double dtw3(const std::vector<double>& values1, const std::vector<double>& values2);
+double dtw1d(const std::vector<double>& values1, const std::vector<double>& values2);
 double dtw2d(const std::vector<double>& values1, const std::vector<double>& values2);
 double custom_distance(const std::vector<double>& values1, const std::vector<double>& values2,
                        const std::string& dist_name);