/********************************************************************
* File: Haar.hh
* Author: Keith Schwarz (htiek@cs.stanford.edu)
*
* Implementation of functions to encode a collection of values
* in the range [0, 1] as Haar basis wavelets. The function takes
* as input a set of values (which must always be a power of two
* in number), then produces a representation of those points as a
* sum of Haar wavelets. A similar function exists to invert the
* encoding.
*
* These functions can be used to do signal compression (by converting
* to a Haar basis, then dropping the lowest-order terms).
*/
#ifndef Haar_Included
#define Haar_Included
#include <vector>
#include <algorithm> // For copy
#include <iterator> // For iterator_traits
/**
* Function: HaarTransform(RandomIterator begin, RandomIterator end,
* RandomIterator out);
* Usage: HaarTransform(v.begin(), v.end(), output.begin());
* ---------------------------------------------------------------------------
* Applies the one-dimensional Haar transform to the elements in the range
* [begin, end), storing the result in the range beginning at out. It is
* assumed that all elements in the range [begin, end) have values in the
* range [0, 1].
*
* Internally, all operations are buffered into a temporary vector, and so it
* is permissible for the input and output ranges to be coincident.
*/
template <typename RandomIterator, typename OutputIterator>
OutputIterator HaarTransform(RandomIterator begin, RandomIterator end,
OutputIterator out);
/**
* Function: HaarTransformInverse(RandomIterator begin, RandomIterator end,
* RandomOutputIterator out);
* Usage: HaarTransformInverse(v.begin(), v.end(), recovered.begin());
* ---------------------------------------------------------------------
* Computes the inverse Haar transform of a range of values [begin, end),
* storing the result in the range beginning with out. Internally, elements
* will be buffered to a vector, so it is permissible for the input and output
* ranges to be coincident.
*/
template <typename RandomIterator, typename OutputIterator>
OutputIterator HaarTransformInverse(RandomIterator begin, RandomIterator end,
OutputIterator out);
/* * * * * Implementation Below This Point * * * * */
template <typename RandomIterator, typename OutputIterator>
OutputIterator HaarTransform(RandomIterator begin, RandomIterator end,
OutputIterator out) {
/* Utility typedefs to make it easier to talk about the underlying type of
* elements being iterated over and the distances between iterators.
*/
typedef typename std::iterator_traits<RandomIterator>::value_type T;
typedef typename std::iterator_traits<RandomIterator>::difference_type diffT;
/* Base case: If there are no elements to transform, we just copy over
* the values from the input range to the output.
*/
if (end - begin <= diffT(1))
return std::copy(begin, end, out);
/* To allow for overlapping input and output regions, we create two temporary
* vectors here, one to hold the smooth components and one to hold the
* detail components. We'll write this back at the end of the function.
* Since each recursive call uses n/2 auxiliary space this way, the total
* space overhead is O(n).
*/
std::vector<T> smooth;
std::vector<T> detail;
/* Fill the buffer in by computing the smooth and detail components. */
for (diffT i = 0; i < (end - begin) / diffT(2); ++i) {
/* Smooth component given by (A + B) / 2 */
smooth.push_back((begin[2 * i] + begin[2 * i + 1]) / 2.0);
/* Detail component given by (A - B + 1) / 2. */
detail.push_back((begin[2 * i] - begin[2 * i + 1]) / 2.0 + T(0.5));
}
/* Next, recursively transform the smooth region with the Haar transform. */
out = HaarTransform(smooth.begin(), smooth.end(), out);
/* Finally, write back the detail region. */
return std::copy(detail.begin(), detail.end(), out);
}
template <typename RandomIterator, typename OutputIterator>
OutputIterator HaarTransformInverse(RandomIterator begin, RandomIterator end,
OutputIterator out) {
/* Utility typedefs to make it easier to talk about the underlying type of
* elements being iterated over and the distances between iterators.
*/
typedef typename std::iterator_traits<RandomIterator>::value_type T;
typedef typename std::iterator_traits<RandomIterator>::difference_type diffT;
/* If the input range is empty. */
if (begin == end)
return out;
/* Maintain a work buffer where the elements that have been expanded out
* so far can live. We will actually allocate two buffers here, since one
* of them will be necessary to hold the elements and one will serve as
* scratch space. Each of them will be large enough to hold the entire
* output range.
*/
std::vector<T> workBuffer(end - begin), outputBuffer(end - begin);
/* Initially, copy over the first element from the input range to the
* work buffer.
*/
workBuffer[0] = *begin;
/* Continuously iterate over a larger and larger prefix, inverting the
* transform. Each iteration doubles the length of the prefix.
*/
for (diffT k = 1; k < end - begin; k *= 2) {
/* Iterate across the first k elements and invert the transformation
* to the output range.
*/
for (diffT i = 0; i < k; ++i) {
/* The S value is what's currently stored in the work buffer, and its
* matching D value can be found in the input range at position i + k.
*
* The formula for inversion is
*
* A = S + (D - 0.5);
* B = S - (D - 0.5);
*/
outputBuffer[2 * i] = workBuffer[i] + (begin[i + k] - T(0.5));
outputBuffer[2 * i + 1] = workBuffer[i] - (begin[i + k] - T(0.5));
}
/* After this process completes, exchange the output buffer and the
* work buffer. Now, the work buffer holds the prefix of S values
* we've unpacked so far.
*/
workBuffer.swap(outputBuffer);
}
/* Finally, once we've expanded everything out, the work buffer holds the
* completely untransformed sequence. We can then copy it to its final
* destination.
*/
return std::copy(workBuffer.begin(), workBuffer.end(), out);
}
#endif