1D Haar Wavelets Screencast Series

Vladimir Kulyukin

---

Introduction

Below is a list of my screencasts on 1D Haar wavelets compiled for easy reference. The Java source code is here. This post briefly explains how to use the source code. These two posts (post 1 and post 2) discuss the how 1D Haar wavelets can be used to detect edges.

---

Simple Haar Wavelets: Part 01: Signal Functions, Samples, & Step Functions

---

Simple Haar Wavelets: Part 02: Basic Haar Wavelet Transform

---

Simple Haar Wavelets: Part 03: Shifts & Dilations of the Basic Haar Transform

---

Simple Haar Wavelets: Part 4a: Ordered Fast Haar Wavelet Transform

---

Simple Haar Wavelets: Part 4b: Ordered Fast Haar Wavelet Transform

---

Screencast 05: Simple Haar Wavelets: Part 05: In-Place 1D Fast Haar Wavelet Transform

---

Simple Haar Wavelets: Part 06: Ordered Fast Inverse Haar Wavelet Transform

---

References

Y. Nievergelt. "Wavelets Made Easy", Birkhauser, Boston, USA, 1999.

---

Fast Multiplication of Polynomials with Fast Fourier Transform and Complex Nth Roots of Unity

Vladimir Kulyukin

---

Introduction

The Python source code for this post is here. Two polynomials of degree n can be represented as arrays of n+1 coefficients. For example, 9 - 10x + 7x^2 + 6x^3 can be represented as [9, -10, 7, 6] and -5 + 4x - 2x^3 can be represented as [-5, 4, 0, -2]. The brute-force method to multiply two polynomials is to multiply each coefficient in the first array by each coefficient in the second array and then add them up, which runs in O(n^2).

def mult_poly(A, B):
    lna, lnb = len(A), len(B)
    n = lna + lnb - 1
    rslt = [0 for i in range(n)]
    for i in range(0, lna):
        for j in range(0, lnb):
            rslt[i+j] += A[i]*B[j]
    return rslt

Multiplying two polynomials by hand gives us -45 + 86x - 75x^2 - 20x^3 + 44x^2 - 14x^5 - 12x^6. Let us evaluate it in the Python interpreter.

>>> poly01 = [9, -10, 7, 6]
>>> poly02 = [-5, 4, 0, -2]
>>> mult_poly(poly01, poly02)
[-45, 86, -75, -20, 44, -14, -12]

A faster way to implement two polynomials is possible if we use the Fast Fourier Transform (FFT) to evaluate the two polynomials at the 2n complex 2n-th roots of unity with the FFT, do the pointwise multiplication of the resulting polynomials, and then use the Inverse Fast Fourier Transform (IFFT) to reconstruct (interpolate) the coefficients of the product polynomial. This algorithm is in O(nlogn).

---

Complex N-th Roots of Unity

To implement the O(nlogn) algorithm outlined above, we need a function to compute complex nth roots of unity. Here is an implementation of the formula for the n-th complex root of unity from Chapter 32 in "Introduction to Algorithms" by Cormen et al.  book. The error argument is added to control the complexity of the real and imaginary coefficients. When the real or imaginary parts of a complex number is less than the error, it is set to 0.

def omega_with_error(k, n, error=0):
    if n <= 0:
        return None
    else:
        x = cmath.exp((2 * cmath.pi * 1j * k)/n)
        xr = x.real
        xi = x.imag
        if abs(xr) < error: xr = 0
        if abs(xi) < error: xi = 0
        return complex(xr, xi)
 
Another commonly used version of omega has -1 in the exponent.

def omega2_with_error(k, n, error=0):
    if n <= 0:
        return None
    else:
        x = cmath.exp(-(2 * cmath.pi * 1j * k)/n)
        xr = x.real
        xi = x.imag
        if abs(xr) < error: xr = 0
        if abs(xi) < error: xi = 0
        return complex(xr, xi)

---

Fast Polynomial Multiplication

It is assumed that the two coefficient arrays, A and B, that represent the polynomials have a length equal to an integral power of 2. If they are not, they need to be padded with zeros at the right. This can be done as follows.

def pad_with_zeros(a):
    ln = len(a)
    padded_a = copy.copy(a)
    if ln > 2 and is_pow_of_2(ln):
        return padded_a
    while True:
        padded_a.append(0)
        if is_pow_of_2(len(padded_a)):
            return padded_a

Since the polynomials are evaluated with 2n complexs 2n-th roots of unity, the coefficients must be padded until they are of length 2n, which can be done with pad_with_zeros_until_required_length.

def pad_with_zeros_until_required_length(a, n):
    padded_a = copy.copy(a)
    while True:
        padded_a.append(0)
        if dft.is_pow_of_2(len(padded_a)):
            return padded_a

With these tools in place, we can implement a function convolve_poly(A, B) to multiply (convolve) two polynomials.

def convolve_poly(A, B, wf=dft.omega2, error=0):
    ## 1. ensure that coefficient arrays A and B are of the same length
    lenA, lenB = len(A), len(B)
    if lenA != lenB:
        raise Exception('Coefficient arrays must be of the same length')
    ## 2. pad A and B until their lengths are an integral power of 2
    padded_A = dft.pad_with_zeros(A)
    padded_B = dft.pad_with_zeros(B)
    ## 3. pad A and B more so that their lengths are 2n, where n
    ## is the length of padded A
    padded_A = pad_with_zeros_until_required_length(A, 2*len(padded_A))
    padded_B = pad_with_zeros_until_required_length(B, 2*len(padded_A))
    ## 4. compute dft of padded A and B
    fft_A = dft.recursive_fft(padded_A, wf=wf, error=error)
    fft_B = dft.recursive_fft(padded_B, wf=wf, error=error)
    ## 5. compute pointwise multiplication of A and B
    C = []
    for a, b in zip(fft_A, fft_B):
        C.append(a*b)
    ## 6. compute inverse dft and return the first 2*lenA-2 elements.
    return dft.recursive_inverse_fft(C, wf=wf, error=error)[0:2*lenA-1]

If necessary, we can split the real and imaginary coefficients, which is done with convolve_poly_split(A, B).

def convolve_poly_split(A, B, wf=dft.omega2, error=0):
    lenA, lenB = len(A), len(B)
    if lenA != lenB:
        raise Exception('Coefficient arrays must be of the same length')
    ## 1. pad A and B with zeros until their lengths are an integral
    ## power of 2
    padded_A = dft.pad_with_zeros(A)
    padded_B = dft.pad_with_zeros(B)
    ## 2. pad A and B more so that their length is 2*n, where n is
    ## the length of padded_A.
    padded_A = pad_with_zeros_until_required_length(padded_A, 2*len(padded_A))
    padded_B = pad_with_zeros_until_required_length(padded_B, 2*len(padded_A))
    ## 3. compute dft of padded A and B
    fft_A = dft.recursive_fft(padded_A, wf=wf, error=error)
    fft_B = dft.recursive_fft(padded_B, wf=wf, error=error)
    ## 4. compute pointwise multiplication of A and B
    C = []
    for a, b in zip(fft_A, fft_B):
        C.append(a*b)
    ## 5. compute inverse dft of A*B and split the output into
    ## real and imaginary numbers.
    reals, imags = dft.recursive_inverse_fft_split(C, wf=wf, error=error)
    ## 6. the degree of A and B is lenA-1. the degree of A*B is
    ## 2*(lenA-1) = 2*lenA - 2.
    return reals[0:2*lenA-1], imags[0:2*lenA-1] 

---

Example

Let us multiply two polynomials -10 + x - x^2 + 7x^3 and 3 - 6x + 8x^3. Doing it by hand gives us
-30 + 63x - 9x^2 - 53x^3 - 34x^4 - 8x^5 + 56x^6.

poly03 = [-10, 1, -1, 7]
poly04 = [3, -6, 0, 8]
print 'POLY03 x POLY04'
poly03xpoly04 = convolve_poly_split(poly03, poly04)
print poly03xpoly04
poly03xpoly04 = convolve_poly(poly03, poly04)
print poly03xpoly04
poly03xpoly04 = mult_poly(poly03, poly04)

Below is the output.

POLY03 x POLY04
([-29.999999999999993, 63.0, -9.0, -53.0, -34.00000000000001, -8.0, 56.0], [-7.105427357601002e-15, 5.452154178752075e-15, 1.887379141862766e-15, -9.004867857552575e-15, -1.4210854715202004e-14, 5.452154178752075e-15, 1.942890293094024e-14])
[(-29.999999999999993-7.105427357601002e-15j),(63+5.452154178752075e-15j), (-9+1.887379141862766e-15j),(-53-9.004867857552575e-15j),(-34.00000000000001-1.4210854715202004e-14j), (-8+5.452154178752075e-15j), (56+1.942890293094024e-14j)]

---

Fast Fourier Transform (FFT) and Inverse Fast Fourier Transform in Python

Vladimir Kulyukin

---

Introduction

This post is an implementation of the fast fourier transform (FFT) and its inverse (IFFT) as outlined in Ch. 32 of "Introduction to Algorithms' by Cormen, Leiserson, and Rivest. The Python source is here. There are four top-level functions: recursive_fft(), recursive_inverse_fft(), recursive_fft_split(),  recursive_inverse_fft_split(). Let us define a sample array.

>>> a = [8, 4, 8, 0]

Now we can call recursive_fft to obtain the fast fourier transform of this array. The second argument is a function that computes complex n-th roots of unity. There are four functions for computing complex roots of unity: omega_with_error(), omega2_with_error(), omega(), and omega2(). For example, omega_with_error() implements the formula for the n-th complex root of unity on p. 784 of Cormen's Algorithms book. The only addition is the error argument: when the real or imaginary parts of a complex number is less than error, it is set to 0.

def omega_with_error(k, n, error=0):
    if n <= 0:
        return None
    else:
        x = cmath.exp((2 * cmath.pi * 1j * k)/n)
        xr = x.real
        xi = x.imag
        if abs(xr) < error: xr = 0
        if abs(xi) < error: xi = 0
        return complex(xr, xi)

The function omega2_with_error() is another formula for the n-th complex root of unity. For example, this formula is used in this FFT calculator. Here is how you can obtain the FFT of the input array. The key "wf" can be read as "omega function."

>>> fft_a = recursive_fft(a, wf=omega_with_error, error=0.001)
>>> fft_a
[20, 4j, 12, -4j]

We can also split the output into the real and imaginary components.

>>> real, imag = recursive_fft_split(fft_a, wf=omega_with_error, error=0.001)
>>> real
[20, 0.0, 12, 0.0]
>>> imag
[0, 4.0, 0, -4.0]

To restore the original array, you can call recursive_inverse_fft():

>>> recursive_inverse_fft(ar, wf=omega_with_error, error=0.001)
[(8+0j), (4+0j), (8+0j), 0j]



The output of the inverse FFT can also be split into real and imaginary components:

>>> recursive_inverse_fft_split(ar, wf=omega_with_error, error=0.001)
([8.0, 4.0, 8.0, 0.0], [0.0, 0.0, 0.0, 0.0])

---

Edge Detection with 1D Inplace Fast Haar Wavelet Transform Applied to Image Columns

Vladimir Kulyukin

---

Problem Statement

In a previous post, we used the 1D inplace fast Haar wavelet transform to design a simple technique to detect edges in images by applying the transform to image rows. In this post, we will detect edges by applying the same transform to image columns. We will again use OpenCV to do the coding.  The implementation of the 1D inplace fast Haar wavelet transform remains unchanged (source code is here; look for the function inPlaceFastHaarWaveletTransform).

---

cv::Mat Manipulation

We need to extract columns from cv::Mat structures. Unlike in this post where we applied the Haar transform to image rows, in this post we apply the transform to columns. Edges are detected by looking at the absolute magnitudes of Haar wavelets at specific frequencies (e.g., 1/2^{n-1}, where n is the length of the image column which has to be an integral power of 2 for the transform to work). We will again use grayscale images that can be constructed in OpenCV as follows: 

cv::Mat mat(n, n, CV_8U, cv::Scalar(0));

Grayscale cv::Mats consist of uchar rows and columns.  Hence, a simple function to convert a column of uchars into an array of doubles.

double* ucharColAryToDoubleAry(cv::Mat &img, int colnum, int num_rows) {
  double *colAry = new double[num_rows];
  for(int r = 0; r < num_rows; r++) {
    colAry[r] = (double)img.at<uchar>(r, colnum);
  }
  return colAry;
}

---

Edge Detection

The function haarDetectVerEdges, given below, uses the highest frequency Haar wavelets computed in the current image column to detect edges. The function goes through the image specified in the first argument column by column. It takes the wavelet bandwidth defined by the next two arguments: lower_theta and upper_theta. The fourth parameter, rv, specifies the value on the scale [0, 255] that we want a pixel to contain if that pixel contains an edge. If the value of the wavelet is negative, then the current pixel contains an edge. If the value of the wavelet is positive, the above pixel contains an edge.

void haarDetectVerEdges(cv::Mat &img, int lower_theta, int upper_theta, int rv)
{
  const int nr = img.rows;
  const int nc = img.cols * img.channels();
  double* converted_col_data;
  double hd = 0.0;

  for(int c = 0; c < nc; c++) {
    converted_col_data = ucharColAryToDoubleAry(img, c, nr);
    inPlaceFastHaarWaveletTransform(converted_col_data, nr);
    for(int r = 0; r < nr; r += 2) {
      hd = converted_col_data[r];    
      if ( std::abs(hd) > lower_theta && std::abs(hd) < upper_theta ) {
    if ( hd < 0 ) {
      img.at<uchar>(r, c) = (uchar)rv;
    }
    else if ( hd > 0 && c > 0 ) {
      img.at<uchar>(r-1, c) = (uchar)rv;
    }
      }
    }
    delete [] converted_col_data;
  }
}

---

Tests

Let us write a few tools tests. The functions createGrayscaleSquareMat(), drawCircle(), etc. are given in this post. The first test, testColHaar01(), creates a grayscale 256 x 256 image and draws several ellipses on it. The created image is shown in Figure 1 and is displayed with cv::imshow(). 

void testColHaar01() {
  cv::Mat img = cv::Mat::zeros(256, 256, CV_8U);
  drawEllipse(img, 90);
  drawEllipse(img, 0);
  drawEllipse(img, 45);
  drawEllipse(img, -45);
  cv::imshow("col orig", img);
  haarDetectVerEdges(img, 125, 128, 255);
  cv::imshow("col edge", img);
  cv::waitKey(0);
}

Figure 1. Image with ellipses

When the function haarDetecVerEdges is applied to the image in Figure 1, the resulting image is given in Figure 2.

Figure 2. Edges detected in Figure 1

The second test, testColHaar02() , draws one shallow and one filled rectangle, as shown in Figure 3, and then applies the 1D inplace fast Haar transform to its columns. The result is shown in Figure 4.

 void testColHaar02() {
  cv::Mat img = cv::Mat::zeros(256, 256, CV_8U);
  cv::rectangle(img, cv::Point(10, 10), cv::Point(40, 40), cv::Scalar(255), 1);
  cv::rectangle(img, cv::Point(40, 40), cv::Point(70, 70), cv::Scalar(255), -1);
  cv::imshow("col orig", img);
  haarDetectHorEdges(img, 125, 128, 255);
  cv::imshow("col edge", img);
  cv::waitKey(0);
}

Figure 3. Image with two squares

Figure 4. Edges detected in Figure3

---

Edge Detection with 1D Inplace Fast Haar Wavelet Transform Applied to Image Rows

Vladimir Kulyukin

---

Problem Statement

In this post, we will use the 1D inplace fast Haar wavelet transform to design a simple technique to detect edges in images. We will use OpenCV to do the coding.  We will start with an implementation of the 1D inplace fast Haar wavelet transform.

const double loge_of_2 = log(2);
const double log10_of_2 = log10(2);
double log2(double x) { return log10(x)/log10_of_2; }

void inPlaceFastHaarWaveletTransformWithNumSweeps(double* row_data, int size, int num_sweeps) {
  int I = 1;
  int GAP_SIZE = 2;
  int NUM_SAMPLE_VALS = size;
  if ( num_sweeps < 1 || num_sweeps > size ) return;
  double a = 0;
  double c = 0;
  for(int SWEEP_NUM = 1; SWEEP_NUM <= num_sweeps; SWEEP_NUM++) {
    NUM_SAMPLE_VALS /= 2;
    for(int K = 0; K < NUM_SAMPLE_VALS; K++) {
      a = (row_data[GAP_SIZE * K] + row_data[GAP_SIZE * K + I])/2;
      c = (row_data[GAP_SIZE * K] - row_data[GAP_SIZE * K + I])/2;
      row_data[GAP_SIZE * K] = a;
      row_data[GAP_SIZE * K + I] = c;
    }
    I = GAP_SIZE;
    GAP_SIZE *= 2;
  }
}

void inPlaceFastHaarWaveletTransform(double* row_data, int n) {
  if ( n == 0 || n % 2 != 0 ) return;
  int num_sweeps = log(n)/loge_of_2;
  inPlaceFastHaarWaveletTransformWithNumSweeps(row_data, n, num_sweeps);
}

We can test this code as follows.

void displayArray(double* ary, int size) {
  for(int i = 0; i < size; i++) { std::cout << ary[i] << " "; }
  std::cout << std::endl << std::flush;
}
double row_data[32] = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 255, 255, 255, 255, 255};
inPlaceFastHaarWaveletTransform(row_data, 32);
displayArray(row_data, 32);

This code produces the following output.

39.8438 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -39.8438 0 0 0 0 0 0 0 -79.6875 0 -63.75 -127.5 -95.625 0 0 0

---

cv::Mat Manipulation

The next technical problem we need to solve is how to extract rows from cv::Mat structures. The above Haar transform can then be applied to each row. Edges can be detected by looking at the absolute magnitudes of Haar wavelets at specific frequencies (e.g., 1/2^{n-1}, where n is the length of the image row which has to be an integral power of 2). We will use grayscale images. cv::Mat mat(n, n, CV_8U, cv::Scalar(0));

If you work in RGB images, check OpenCV documentation for appropriate conversions. The cv::Mat structures consist of uchar rows and columns.  Hence, a simple function to convert a row of uchars into an array of doubles.

double* ucharAryToDoubleAry(uchar* row_data, int size) {
  double* rslt = new double[size];
  for(int i = 0; i < size; i++) {
    rslt[i] = (double)row_data[i];
  }
  return rslt;
}

---

Actual Detection of Edges
The function, given below, that detects edges uses the highest frequency Haar wavelets computed in the current image row. The function goes through the image specified in the first argument row by row. It takes the wavelet bandwidth defined by the next two arguments: lower_theta and upper_theta. The fourth parameter, rv, specifies the value on the scale [0, 255] that we want a pixel to contain if that pixel contains an edge. If the value of the wavelet is negative, then the current pixel contains an edge. If the value of the wavelet is positive, the previous pixel contains an edge.

void haarDetectHorEdges(cv::Mat &img, int lower_theta, int upper_theta, int rv)
{
  const int nr = img.rows;
  const int nc = img.cols * img.channels();
  uchar* row_data;
  double* converted_row_data;
  double hd = 0.0;

  for(int r = 0; r < nr; r++) {
    row_data = img.ptr<uchar>(r);
    converted_row_data = ucharAryToDoubleAry(row_data, nc);
    inPlaceFastHaarWaveletTransform(converted_row_data, nc);
    for(int c = 0; c < nc; c += 2) {
      hd = converted_row_data[c];     
      if ( std::abs(hd) > lower_theta && std::abs(hd) < upper_theta ) {
    if ( hd < 0 ) {
      img.at<uchar>(r, c) = (uchar)rv;
    }
    else if ( hd > 0 && c > 0 ) {
      img.at<uchar>(r, c-1) = (uchar)rv;
    }
      }
    }
    delete [] converted_row_data;
  }
}

---

Tests
Now on to a few simple tools in OpenCV that we can use in our tests. One thing to keep in mind is that the length and width of each image must be equal to an integral power of 2. Otherwise, the Haar transform will not work.

cv::Mat createGrayscaleSquareMat(int n) {
  cv::Mat mat(n, n, CV_8U, cv::Scalar(0));
  return mat;
}

cv::Mat createGrayscaleSquareMatWithLine(int n, int y1, int x1, int y2, int x2) {
  cv::Mat mat(n, n, CV_8U, cv::Scalar(0));
  cv::Point p1(y1, x1);
  cv::Point p2(y2, x2);
  cv::line(mat, p1, p2, cv::Scalar(255), 1, 8);
  return mat;
}

void drawCircle(cv::Mat &img, int center_y1, int center_x1, int rad, int thickness) {
  circle(img,
     cv::Point(center_y1, center_x1),
     rad,
     cv::Scalar(255),
     thickness,
     8);
}

void drawEllipse(cv::Mat &img, double angle)
{
  int thickness = 2;
  int lineType = 8;
  const int nr = img.rows;
  cv::ellipse(img,
          cv::Point((int)nr/2, (int)nr/2),
          cv::Size((int)nr/4, (int)nr/16),
          angle,
          0,
          360,
          cv::Scalar(255),
          thickness,
          lineType);
}

The first test, test01(), creates a grayscale 256 x 256 image and draws a line in it. The created image is shown in Figure 1 and is displayed with cv::imshow(). 

void test01() {
  cv::Mat img = createGrayscaleSquareMatWithLine(256, 10, 10, 120, 120);
  cv::imshow("orig", img);
  haarDetectHorEdges(img, 125, 128, 255);
  cv::imshow("edge", img);
  cv::waitKey(0);
}
 

Figure 1. A grayscale image with a line

 When the function haarDetectHorEdges is applied, the resulting image is shown in Figure 2.

Figure 2. Edges detected in Figure 1

The second test, test02(), draws three white circles in a 256 x 256 grayscale image, as shown in Figure 3, and detects edges, as shown in Figure 4.

void test02() {
  cv::Mat img = createGrayscaleSquareMat(256);
  drawCircle(img, 170, 170, 50, 5);
  drawCircle(img, 120, 120, 50, 5);
  drawCircle(img, 170, 120, 50, 5);
  cv::imshow("orig", img);
  haarDetectHorEdges(img, 125, 128, 255);
  cv::imshow("edge", img);
  cv::waitKey(0);
}

 

Figure 3. Three white circles in a grayscale image

Figure 4. Edges detected in Figure 3

The third test, test03(), draws four ellipses circles in a 256 x 256 grayscale image, as shown in Figure 5, and detects edges, as shown in Figure 6.
 
void test03() {
  cv::Mat img = cv::Mat::zeros(256, 256, CV_8U);
  drawEllipse(img, 90);
  drawEllipse(img, 0);
  drawEllipse(img, 45);
  drawEllipse(img, -45);
  cv::imshow("orig", img);
  haarDetectHorEdges(img, 125, 128, 255);
  cv::imshow("edge", img);
  cv::waitKey(0);
}

Figure 5. Four ellipses in a grayscale image

Figure 6. Edges detected in Figure 5

One final note. We have deliberately simplified our images to show how this method works. Grayscale images with a greater variety of pixel values will require more sophisticated thresholds and wavelet frequencies.

---