wavelets.cc File Reference

#include <catch2/catch_all.hpp>
#include <memory>
#include <random>
#include "sopt/types.h"
#include "sopt/wavelets/direct.h"
#include "sopt/wavelets/indirect.h"
#include "sopt/wavelets/wavelet_data.h"
#include "sopt/wavelets/wavelets.h"

Include dependency graph for wavelets.cc:

Go to the source code of this file.

Typedefs
using	t_iVector = sopt::Array< sopt::t_uint >

Functions
t_iVector	even (t_iVector const &x)
t_iVector	odd (t_iVector const &x)
template<typename T >
Eigen::Array< typename T::Scalar, T::RowsAtCompileTime, T::ColsAtCompileTime >	upsample (Eigen::ArrayBase< T > const &input)
sopt::t_int	random_integer (sopt::t_int min, sopt::t_int max)
t_iVector	random_ivector (sopt::t_int size, sopt::t_int min, sopt::t_int max)
template<typename T0 >
void	check_round_trip (Eigen::ArrayBase< T0 > const &input_, sopt::t_uint db, sopt::t_uint nlevels=1)
	TEST_CASE ("wavelet data")
	TEST_CASE ("Wavelet transform innards with integer data", "[wavelet]")
	TEST_CASE ("1D wavelet transform with floating point data", "[wavelet]")
	TEST_CASE ("1D wavelet transform with complex data", "[wavelet]")
	TEST_CASE ("2D wavelet transform with real data", "[wavelet]")
	TEST_CASE ("Functor implementation", "[wavelet]")
	TEST_CASE ("Automatic input resizing", "[wavelet]")
	TEST_CASE ("Dirac wavelets")

Typedef Documentation

t_iVector

using t_iVector = sopt::Array<sopt::t_uint>

Definition at line 11 of file wavelets.cc.

Function Documentation

check_round_trip()

template<typename T0 >

void check_round_trip	(	Eigen::ArrayBase< T0 > const &	input_,
		sopt::t_uint	db,
		sopt::t_uint	nlevels = 1
	)

Definition at line 49 of file wavelets.cc.

                                              {
  auto const input = input_.eval();
  auto const &dbwave = sopt::wavelets::daubechies_data(db);
  auto const transform = sopt::wavelets::direct_transform(input, nlevels, dbwave);
  auto const actual = sopt::wavelets::indirect_transform(transform, nlevels, dbwave);
  CAPTURE(actual);
  CAPTURE(input);
  CAPTURE(transform);
  CHECK(input.isApprox(actual, 1e-14));
  CHECK(not transform.isApprox(sopt::wavelets::direct_transform(input, nlevels - 1, dbwave), 1e-4));
}

References sopt::wavelets::daubechies_data(), sopt::wavelets::direct_transform(), and sopt::wavelets::indirect_transform().

Referenced by TEST_CASE().

even()

t_iVector even

(

t_iVector const &

x

)

Definition at line 12 of file wavelets.cc.

                                   {
  t_iVector result((x.size() + 1) / 2);
  for (t_iVector::Index i(0); i < x.size(); i += 2) result(i / 2) = x(i);
  return result;
};

Referenced by TEST_CASE().

odd()

t_iVector odd

(

t_iVector const &

x

)

Definition at line 17 of file wavelets.cc.

                                  {
  t_iVector result(x.size() / 2);
  for (t_iVector::Index i(1); i < x.size(); i += 2) result(i / 2) = x(i);
  return result;
};

Referenced by TEST_CASE().

random_integer()

sopt::t_int random_integer	(	sopt::t_int	min,
		sopt::t_int	max
	)

Definition at line 34 of file wavelets.cc.

                                                       {
  extern std::unique_ptr<std::mt19937_64> mersenne;
  std::uniform_int_distribution<sopt::t_int> uniform_dist(min, max);
  return uniform_dist(*mersenne);
};

References mersenne().

Referenced by TEST_CASE().

random_ivector()

t_iVector random_ivector	(	sopt::t_int	size,
		sopt::t_int	min,
		sopt::t_int	max
	)

Definition at line 39 of file wavelets.cc.

                                                                       {
  extern std::unique_ptr<std::mt19937_64> mersenne;
  t_iVector result(size);
  std::uniform_int_distribution<sopt::t_int> uniform_dist(min, max);
  for (t_iVector::Index i(0); i < result.size(); ++i) result(i) = uniform_dist(*mersenne);
  return result;
};

References mersenne().

Referenced by TEST_CASE().

TEST_CASE() [1/8]

TEST_CASE	(	"1D wavelet transform with complex data"	,
		""	[wavelet]
	)

Definition at line 244 of file wavelets.cc.

                                                                 {
  using namespace sopt;
  using namespace sopt::wavelets;
  SECTION("Round-trip test for complex data") {
    auto input = Array<t_complex>::Random(random_integer(2, 100) * 2).eval();
    auto const &dbwave = daubechies_data(random_integer(1, 38));
    auto const actual = indirect_transform(direct_transform(input, 1, dbwave), 1, dbwave);
    CHECK(input.isApprox(actual, 1e-14));
    CHECK(not input.isApprox(direct_transform(input, 1, dbwave), 1e-4));
  }
}

References sopt::wavelets::daubechies_data(), sopt::wavelets::direct_transform(), sopt::wavelets::indirect_transform(), and random_integer().

TEST_CASE() [2/8]

TEST_CASE	(	"1D wavelet transform with floating point data"	,
		""	[wavelet]
	)

Definition at line 189 of file wavelets.cc.

                                                                        {
  using namespace sopt;
  using namespace sopt::wavelets;
 
  Image<> const data = Image<>::Random(16, 16);
  auto const &wavelet = daubechies_data(4);
 
  // Condition on input fixture data
  REQUIRE((data.rows() % 2 == 0 and (data.cols() == 1 or data.cols() % 2 == 0)));
 
  SECTION("Direct transform == two downsample + convolution") {
    auto const actual = direct_transform(data.row(0), 1, wavelet);
    Array<> high(data.cols() / 2);
    Array<> low(data.cols() / 2);
    down_convolve(high, data.row(0), wavelet.direct_filter.high);
    down_convolve(low, data.row(0), wavelet.direct_filter.low);
    CHECK(low.transpose().isApprox(actual.head(data.row(0).size() / 2)));
    CHECK(high.transpose().isApprox(actual.tail(data.row(0).size() / 2)));
  }
 
  SECTION("Indirect transform == two upsample + convolution") {
    auto const actual = indirect_transform(data.row(0).transpose(), 1, wavelet);
    auto const low = upsample(data.row(0).transpose().head(data.rows() / 2));
    auto const high = upsample(data.row(0).transpose().tail(data.rows() / 2));
    auto expected = copy(data.row(0).transpose());
    convolve_sum(expected, low, wavelet.direct_filter.low.reverse(), high,
                 wavelet.direct_filter.high.reverse());
    CAPTURE(expected.transpose());
    CAPTURE(actual.transpose());
    CHECK(expected.isApprox(actual));
  }
 
  SECTION("Round-trip test for single level") {
    for (t_int i(0); i < 20; ++i) {
      check_round_trip(Array<>::Random(random_integer(2, 100) * 2), random_integer(1, 38), 1);
    }
  }
 
  SECTION("Round-trip test for two levels") {
    check_round_trip(Array<>::Random(8), 1, 2);
    check_round_trip(Array<>::Random(8), 2, 2);
    check_round_trip(Array<>::Random(16), 4, 2);
    check_round_trip(Array<>::Random(52), 10, 2);
  }
 
  t_uint constexpr nlevels = 5;
  SECTION("Round-trip test for multiple levels") {
    for (t_int i(0); i < 10; ++i) {
      auto const n = random_integer(2, nlevels);
      check_round_trip(Array<>::Random(random_integer(2, 100) * (1u << n)), random_integer(1, 38),
                       n);
    }
  }
}

References check_round_trip(), sopt::wavelets::daubechies_data(), sopt::wavelets::direct_transform(), sopt::wavelets::indirect_transform(), n, random_integer(), and upsample().

TEST_CASE() [3/8]

TEST_CASE	(	"2D wavelet transform with real data"	,
		""	[wavelet]
	)

Definition at line 256 of file wavelets.cc.

                                                              {
  using namespace sopt;
  using namespace sopt::wavelets;
  SECTION("Single level round-trip test for square matrix") {
    auto N = random_integer(2, 100) * 2;
    check_round_trip(Image<>::Random(N, N), random_integer(1, 38), 1);
  }
  SECTION("Single level round-trip test for non-square matrix") {
    auto Nx = random_integer(2, 5) * 2;
    auto Ny = Nx + 5 * 2;
    check_round_trip(Image<>::Random(Nx, Ny), random_integer(1, 38), 1);
  }
  SECTION("Round-trip test for multiple levels") {
    for (t_int i(0); i < 10; ++i) {
      auto const n = random_integer(2, 5);
      auto const Nx = random_integer(2, 5) * (1u << n);
      auto const Ny = random_integer(2, 5) * (1u << n);
      check_round_trip(Image<>::Random(Nx, Ny), random_integer(1, 38), n);
    }
  }
}

References check_round_trip(), n, N, and random_integer().

TEST_CASE() [4/8]

TEST_CASE	(	"Automatic input resizing"	,
		""	[wavelet]
	)

Definition at line 295 of file wavelets.cc.

                                                   {
  using namespace sopt;
  auto const wavelet = wavelets::factory("DB3", 4);
  auto const input = Image<t_complex>::Random(256, 128).eval();
  Image<t_complex> output(1, 1);
  wavelet.direct(output, input);
  CHECK(output.rows() == input.rows());
  CHECK(output.cols() == input.cols());
 
  output.resize(1, 1);
  wavelet.indirect(input, output);
  CHECK(output.rows() == input.rows());
  CHECK(output.cols() == input.cols());
}

References sopt::wavelets::factory().

TEST_CASE() [5/8]

TEST_CASE

(

"Dirac wavelets"

)

Definition at line 310 of file wavelets.cc.

                            {
  using namespace sopt;
  auto const wavelet = wavelets::factory("Dirac");
  Image<t_complex> const input = Image<t_complex>::Random(256, 128);
  Image<t_complex> output(1, 1);
 
  wavelet.direct(output, input);
  CHECK(output.isApprox(input));
 
  output = Image<t_complex>::Zero(1, 1);
  wavelet.indirect(input, output);
  CHECK(output.isApprox(input));
}

References sopt::wavelets::factory().

TEST_CASE() [6/8]

TEST_CASE	(	"Functor implementation"	,
		""	[wavelet]
	)

Definition at line 278 of file wavelets.cc.

                                                 {
  using namespace sopt;
  auto const wavelet = wavelets::factory("DB3", 4);
  auto const input = Image<t_complex>::Random(256, 128).eval();
  SECTION("Normal instances") {
    auto const transform = wavelet.direct(input);
    CHECK(transform.isApprox(wavelets::direct_transform(input, wavelet.levels(), wavelet)));
    CHECK(input.isApprox(wavelet.indirect(transform)));
  }
  SECTION("Expression instances") {
    Image<t_complex> output(2, input.cols());
    wavelet.direct(output.row(0).transpose(), input.row(0).transpose());
    wavelet.indirect(output.row(0).transpose(), output.row(1).transpose());
    CHECK(input.row(0).isApprox(output.row(1)));
  }
}

References sopt::wavelets::direct_transform(), and sopt::wavelets::factory().

TEST_CASE() [7/8]

TEST_CASE

(

"wavelet data"

)

Definition at line 62 of file wavelets.cc.

                          {
  for (sopt::t_int num = 1; num < 100; num++) {
    if (num < 39)
      REQUIRE(sopt::wavelets::daubechies_data(num).coefficients.size() == 2 * num);
    else
      REQUIRE_THROWS(sopt::wavelets::daubechies_data(num));
  }
}

References sopt::wavelets::daubechies_data().

TEST_CASE() [8/8]

TEST_CASE	(	"Wavelet transform innards with integer data"	,
		""	[wavelet]
	)

Definition at line 71 of file wavelets.cc.

                                                                      {
  using namespace sopt::wavelets;
 
  t_iVector small(3);
  small << 1, 2, 3;
  t_iVector large(6);
  large << 4, 5, 6, 7, 8, 9;
 
  SECTION("Periodic scalar product") {
    // no wrapping
    CHECK(periodic_scalar_product(large, small, 0) == 1 * 4 + 2 * 5 + 3 * 6);
    CHECK(periodic_scalar_product(large, small, 1) == 1 * 5 + 2 * 6 + 3 * 7);
    CHECK(periodic_scalar_product(large, small, 3) == 1 * 7 + 2 * 8 + 3 * 9);
 
    // with wrapping
    CHECK(periodic_scalar_product(large, small, 4) == 1 * 8 + 2 * 9 + 3 * 4);
    // with wrapping and expression
    CHECK(periodic_scalar_product(large, small.reverse(), 4) == 3 * 8 + 2 * 9 + 1 * 4);
    // wrapping works with offset as well
    CHECK(periodic_scalar_product(large, small, 4 + large.size()) == 1 * 8 + 2 * 9 + 3 * 4);
    CHECK(periodic_scalar_product(large, small, 4 - 3 * large.size()) == 1 * 8 + 2 * 9 + 3 * 4);
 
    // signal smaller than filter
    CHECK(periodic_scalar_product(small, large.head(4), 1) == 4 * 2 + 5 * 3 + 6 * 1 + 7 * 2);
  }
 
  SECTION("Convolve") {
    t_iVector result(large.size());
 
    convolve(result, large, small);
 
    CHECK(result(0) == 1 * 4 + 2 * 5 + 3 * 6);
    CHECK(result(1) == 1 * 5 + 2 * 6 + 3 * 7);
    CHECK(result(3) == 1 * 7 + 2 * 8 + 3 * 9);
    CHECK(result(4) == 1 * 8 + 2 * 9 + 3 * 4);
  }
 
  SECTION("Convolve and sum") {
    t_iVector result(large.size());
    t_iVector noOffset(large.size());
 
    // Check that if high pass is zero, then this is an offseted convolution
    convolve_sum(result, large, small, large, 0 * small);
    convolve(noOffset, large, small);
    CHECK(result(small.size() - 1) == noOffset(0));
    CHECK(result(0) == noOffset(result.size() - small.size() + 1));
 
    // Check same for low pass
    convolve_sum(result, large, 0 * small, large, small);
    CHECK(result(small.size() - 1) == noOffset(0));
    CHECK(result(0) == noOffset(result.size() - small.size() + 1));
 
    // Check symmetry relationships
    auto const trial = [&small, &large](int a, int b, int c, int d) {
      t_iVector result(large.size());
      convolve_sum(result, a * large, b * small, c * large, d * small);
      return result;
    };
 
    // should all be ok as long as arguments sum: (a * b) + (c * d) == (a' * b') + (c' * d')
    CHECK((trial(0, 1, 3, 1) == trial(0, 1, 1, 3)).all());
    CHECK((trial(5, 1, 3, 1) == trial(3, 1, 5, 1)).all());
    CHECK((trial(1, 5, 3, 1) == trial(3, 1, 5, 1)).all());
    CHECK((trial(1, 3, 5, 1) == trial(3, 1, 5, 1)).all());
    CHECK((trial(1, 3, 1, 5) == trial(3, 1, 5, 1)).all());
    CHECK((trial(1, 0, 4, 2) == trial(3, 1, 5, 1)).all());
    CHECK((trial(1, -1, 1, 1) == trial(0, 1, 0, 1)).all());
    CHECK((trial(4, -3, 2, 6) == trial(0, 1, 0, 1)).all());
  }
 
  SECTION("Convolve and Down-sample simultaneously") {
    t_iVector expected(large.size());
    convolve(expected, large, small);
    t_iVector actual(large.size() / 2);
    down_convolve(actual, large, small);
    for (size_t i(0); i < static_cast<size_t>(actual.size()); ++i)
      CHECK(expected(i * 2) == actual(i));
  }
 
  SECTION("Convolve output to expression") {
    t_iVector actual(large.size() * 2);
    t_iVector expected(large.size());
    convolve(actual.head(large.size()), large, small);
    convolve(expected, large, small);
    CHECK((actual.head(large.size()) == expected).all());
  }
 
  SECTION("Copy does copy") {
    auto result = copy(large);
    CHECK(large.data() != result.data());
 
    auto actual = copy(large.head(3));
    CHECK(large.data() != actual.data());
    CHECK(large.data() == large.head(3).data());
  }
 
  SECTION("Convolve, Sum and Up-sample simultaneously") {
    for (sopt::t_int i(0); i < 100; ++i) {
      auto const Ncoeffs = random_integer(2, 10) * 2;
      auto const Nfilters = random_integer(2, 5);
      auto const Nhead = Ncoeffs / 2;
      auto const Ntail = Ncoeffs - Nhead;
 
      auto const coeffs = random_ivector(Ncoeffs, -10, 10);
      auto const low = random_ivector(Nfilters, -10, 10);
      auto const high = random_ivector(Nfilters, -10, 10);
 
      t_iVector actual(Ncoeffs);
      t_iVector expected(Ncoeffs);
      // does all in go, more complicated but compuationally less intensive
      up_convolve_sum(actual, coeffs, even(low), odd(low), even(high), odd(high));
      // first up-samples, then does convolve: conceptually simpler but does unnecessary operations
      convolve_sum(expected, upsample(coeffs.head(Nhead)), low, upsample(coeffs.tail(Ntail)), high);
      CHECK((actual == expected).all());
    }
  }
}

References a, b, even(), odd(), random_integer(), random_ivector(), and upsample().

upsample()

template<typename T >

Eigen::Array<typename T::Scalar, T::RowsAtCompileTime, T::ColsAtCompileTime> upsample

(

Eigen::ArrayBase< T > const &

input

)

Definition at line 23 of file wavelets.cc.

                                    {
  using Matrix = Eigen::Array<typename T::Scalar, T::RowsAtCompileTime, T::ColsAtCompileTime>;
  Matrix result(input.size() * 2);
  for (t_iVector::Index i(0); i < input.size(); ++i) {
    result(2 * i) = input(i);
    result(2 * i + 1) = 0;
  }
  return result;
};

Referenced by TEST_CASE().