Wavelet computing improvements

dali/kernels/signal/wavelet/mother_wavelet.cu

@@ -12,158 +12,117 @@
 // See the License for the specific language governing permissions and
 // limitations under the License.
 
-#include <vector>
+#include <cmath>
 #include "dali/kernels/signal/wavelet/mother_wavelet.cuh"
+#include "dali/core/math_util.h"
 
 namespace dali {
 namespace kernels {
 namespace signal {
 
 template <typename T>
-__device__
-T HaarWavelet(T t, T a, T b) {
-  T x = std::pow(2.0, a) - b;
-  if (0.0 <= x && x < 0.5) {
-    return std::pow(2.0, a / 2.0);
+HaarWavelet<T>::HaarWavelet(const std::vector<T> &args) {
+  if (args.size() != 0) {
+    throw new std::invalid_argument("HaarWavelet doesn't accept any arguments.");
   }
-  if (0.5 <= x && x < 1.0) {
-    return -std::pow(2.0, a / 2.0);
-  }
-  return 0.0;
 }
 
 template <typename T>
-__device__
-T DaubechiesWavelet(T t, T a, T b) {
+__device__ T HaarWavelet<T>::operator()(const T &t) const {
+  if (0.0 <= t && t < 0.5) {
+    return 1.0;
+  }
+  if (0.5 <= t && t < 1.0) {
+    return -1.0;
+  }
   return 0.0;
 }
 
-template <typename T>
-__device__
-T SymletWavelet(T t, T a, T b) {
-  return 0.0;
-}
+template class HaarWavelet<float>;
+template class HaarWavelet<double>;
 
 template <typename T>
-__device__
-T CoifletWavelet(T t, T a, T b) {
-  return 0.0;
+MeyerWavelet<T>::MeyerWavelet(const std::vector<T> &args) {
+  if (args.size() != 0) {
+    throw new std::invalid_argument("MeyerWavelet doesn't accept any arguments.");
+  }
 }
 
 template <typename T>
-__device__
-T BiorthogonalWavelet(T t, T a, T b) {
-  return 0.0;
+__device__ T MeyerWavelet<T>::operator()(const T &t) const {
+  T psi1 = (4/(3*M_PI)*t*std::cos((2*M_PI)/3*t)-1/M_PI*std::sin((4*M_PI)/3*t))/(t-16/9*std::pow(t, 3.0));
+  T psi2 = (8/(3*M_PI)*t*std::cos(8*M_PI/3*t)+1/M_PI*std::sin((4*M_PI)/3)*t)/(t-64/9*std::pow(t, 3.0));
+  return psi1 + psi2;
 }
 
-template <typename T>
-__device__
-T MeyerWavelet(T t, T a, T b) {
-  return 0.0;
-}
+template class MeyerWavelet<float>;
+template class MeyerWavelet<double>;
 
 template <typename T>
-__device__
-T GaussianWavelet(T t, T a, T b) {
-  return 0.0;
+MexicanHatWavelet<T>::MexicanHatWavelet(const std::vector<T> &args) {
+  if (args.size() != 1) {
+    throw new std::invalid_argument("MexicanHatWavelet accepts exactly one argument - sigma.");
+  }
+  this->sigma = *args.begin();
 }
 
 template <typename T>
-__device__
-T MexicanHatWavelet(T t, T a, T b) {
-  return 0.0;
+__device__ T MexicanHatWavelet<T>::operator()(const T &t) const {
+  return 2/(std::sqrt(3*sigma)*std::pow(M_PI, 0.25))*(1-std::pow(t/sigma, 2.0))*std::exp(-std::pow(t, 2.0)/(2*std::pow(sigma, 2.0)));
 }
 
+template class MexicanHatWavelet<float>;
+template class MexicanHatWavelet<double>;
+
 template <typename T>
-__device__
-T MorletWavelet(T t, T a, T b) {
-  return 0.0;
+MorletWavelet<T>::MorletWavelet(const std::vector<T> &args) {
+  if (args.size() != 1) {
+    throw new std::invalid_argument("MorletWavelet accepts exactly 1 argument - C.");
+  }
+  this->C = *args.begin();
 }
 
 template <typename T>
-__device__
-T ComplexGaussianWavelet(T t, T a, T b) {
-  return 0.0;
+__device__ T MorletWavelet<T>::operator()(const T &t) const {
+  return C * std::exp(-std::pow(t, 2.0)) * std::cos(5 * t);
 }
 
+template class MorletWavelet<float>;
+template class MorletWavelet<double>;
+
 template <typename T>
-__device__
-T ShannonWavelet(T t, T a, T b) {
-  return 0.0;
+ShannonWavelet<T>::ShannonWavelet(const std::vector<T> &args) {
+  if (args.size() != 0) {
+    throw new std::invalid_argument("ShannonWavelet doesn't accept any arguments.");
+  }
 }
 
 template <typename T>
-__device__
-T FbspWavelet(T t, T a, T b) {
-  return 0.0;
+__device__ T ShannonWavelet<T>::operator()(const T &t) const {
+  return sinc(t - 0.5) - 2 * sinc(2 * t - 1);
 }
 
+template class ShannonWavelet<float>;
+template class ShannonWavelet<double>;
+
 template <typename T>
-__device__
-T ComplexMorletWavelet(T t, T a, T b) {
-  return 0.0;
+FbspWavelet<T>::FbspWavelet(const std::vector<T> &args) {
+  if (args.size() != 0) {
+    throw new std::invalid_argument("FbspWavelet accepts exactly 3 arguments -> m, fb, fc in that order.");
+  }
+  this->m = *args.begin();
+  this->fb = *(args.begin()+1);
+  this->fc = *(args.begin()+2);
 }
 
 template <typename T>
-MotherWavelet<T>::MotherWavelet(const WaveletName& name) {
-  switch(name) {
-    case WaveletName::HAAR:
-    waveletFunc = &HaarWavelet;
-    break;
-
-    case WaveletName::DB:
-    waveletFunc = &DaubechiesWavelet;
-    break;
-
-    case WaveletName::SYM:
-    waveletFunc = &SymletWavelet;
-    break;
-
-    case WaveletName::COIF:
-    waveletFunc = &CoifletWavelet;
-    break;
-
-    case WaveletName::BIOR:
-    waveletFunc = &BiorthogonalWavelet;
-    break;
-
-    case WaveletName::MEY:
-    waveletFunc = &MeyerWavelet;
-    break;
-
-    case WaveletName::GAUS:
-    waveletFunc = &GaussianWavelet;
-    break;
-
-    case WaveletName::MEXH:
-    waveletFunc = &MexicanHatWavelet;
-    break;
-
-    case WaveletName::MORL:
-    waveletFunc = &MorletWavelet;
-    break;
-
-    case WaveletName::CGAU:
-    waveletFunc = &ComplexGaussianWavelet;
-    break;
-
-    case WaveletName::SHAN:
-    waveletFunc = &ShannonWavelet;
-    break;
-
-    case WaveletName::FBSP:
-    waveletFunc = &FbspWavelet;
-    break;
-
-    case WaveletName::CMOR:
-    waveletFunc = &ComplexMorletWavelet;
-    break;
-
-    default:
-    throw new std::invalid_argument("Unknown wavelet name.");
-  }
+__device__ T FbspWavelet<T>::operator()(const T &t) const {
+  return std::sqrt(fb)*std::pow(sinc(t/std::pow(fb, m)), m)*std::exp(2*M_PI*fc*t);
 }
 
+template class FbspWavelet<float>;
+template class FbspWavelet<double>;
+
 }  // namespace signal
 }  // namespace kernel
 }  // namespace dali

dali/kernels/signal/wavelet/mother_wavelet.cuh

@@ -21,40 +21,113 @@
 #include "dali/core/util.h"
 #include "dali/kernels/kernel.h"
 
+#include <vector>
+
 namespace dali {
 namespace kernels {
 namespace signal {
 
-enum class WaveletName {
-  HAAR,
-  DB,
-  SYM,
-  COIF,
-  BIOR,
-  MEY,
-  GAUS,
-  MEXH,
-  MORL,
-  CGAU,
-  SHAN,
-  FBSP,
-  CMOR
+// wavelets are represented by functors
+// they can store any necessary parameters
+// they must overload () operator
+
+template <typename T>
+class HaarWavelet {
+  static_assert(std::is_floating_point<T>::value,
+    "Data type should be floating point");
+ public:
+  HaarWavelet() = default;
+  HaarWavelet(const std::vector<T> &args);
+  ~HaarWavelet() = default;
+
+  __device__ T operator()(const T &t) const;
+};
+
+template <typename T>
+class MeyerWavelet {
+  static_assert(std::is_floating_point<T>::value,
+    "Data type should be floating point");
+ public:
+  MeyerWavelet() = default;
+  MeyerWavelet(const std::vector<T> &args);
+  ~MeyerWavelet() = default;
+
+  __device__ T operator()(const T &t) const;
+};
+
+template <typename T>
+class GaussianWavelet {
+  static_assert(std::is_floating_point<T>::value,
+    "Data type should be floating point");
+ public:
+  GaussianWavelet() = default;
+  GaussianWavelet(const std::vector<T> &args);
+  ~GaussianWavelet() = default;
+
+  __device__ T operator()(const T &t) const;
+
+ private:
+  uint8_t N;
 };
 
 template <typename T>
-class MotherWavelet {
+class MexicanHatWavelet {
   static_assert(std::is_floating_point<T>::value,
     "Data type should be floating point");
+ public:
+  MexicanHatWavelet() = default;
+  MexicanHatWavelet(const std::vector<T> &args);
+  ~MexicanHatWavelet() = default;
+
+  __device__ T operator()(const T &t) const;
+
+ private:
+  T sigma;
+};
+
+template <typename T>
+class MorletWavelet {
+  static_assert(std::is_floating_point<T>::value,
+    "Data type should be floating point");
+ public:
+  MorletWavelet() = default;
+  MorletWavelet(const std::vector<T> &args);
+  ~MorletWavelet() = default;
+
+  __device__ T operator()(const T &t) const;
+
+ private:
+  T C;
+};
 
+template <typename T>
+class ShannonWavelet {
+  static_assert(std::is_floating_point<T>::value,
+    "Data type should be floating point");
  public:
-  MotherWavelet(const WaveletName &name);
-  ~MotherWavelet() = default;
+  ShannonWavelet() = default;
+  ShannonWavelet(const std::vector<T> &args);
+  ~ShannonWavelet() = default;
 
-  __device__ T (*waveletFunc)(T t, T a, T b);
+  __device__ T operator()(const T &t) const;
 };
 
-template class MotherWavelet<float>;
-template class MotherWavelet<double>;
+template <typename T>
+class FbspWavelet {
+  static_assert(std::is_floating_point<T>::value,
+    "Data type should be floating point");
+ public:
+  FbspWavelet() = default;
+  FbspWavelet(const std::vector<T> &args);
+  ~FbspWavelet() = default;
+
+  __device__ T operator()(const T &t) const;
+
+ private:
+  T m;
+  T fb;
+  T fc;
+};
 
 }  // namespace signal
 }  // namespace kernel