3D Plots of Toroidal Spirals in Octave

Vladimir Kulyukin

---

Problem
  

A toroidal spiral is a space curve whose x, y, and z coordinates are computed with the following parametric equations: 

x = (C + sin(Kt))cos(t); y = (C + sin(Kt))sin(t); z = cos(20t). They are called toroidal because they lie on a torus, as shown in Figure 1. In this post, we will briefly tackle the problem of 3D plotting toroidal spirals in Octave and investigate the impact of the C and K parameters.

 Figure 1. 3D plot of a Toroidal Spiral

---

Solution Outline

The Octave code below shows how the plot in Figure 1 can be generated in Octave. It also shows a generic approach and allows us to play with the C and K parameters.

t = 0:0.01:2*pi;
K = 20;
C = 2;

x2 = arrayfun(@(x) (C + sin(K*x))*cos(x), t);
y2 = arrayfun(@(x) (C + sin(K*x))*sin(x), t);
z2 = arrayfun(@(x) cos(K*x), t);

figure;
plot3(x2, y2, z2);
xlabel('x=(2+sin(20t))cos(t)');
ylabel('y=(2+sin(20t))sin(t)');
zlabel('z=cos(20t)');
view([30, 40]);

Figure 2 shows another 3D plot  with C=4 and K=20.

 

Figure 2. 3D Plot of Another Toroidal Curve

---

Octave Source

t = 0:0.01:2*pi;
K = 20;

C = 0;
x0 = arrayfun(@(x) (C + sin(K*x))*cos(x), t);
y0 = arrayfun(@(x) (C + sin(K*x))*sin(x), t);
z0 = arrayfun(@(x) cos(K*x), t);

figure;
plot3(x0, y0, z0);
xlabel('x=(0+sin(20t))cos(t)');
ylabel('y=(0+sin(20t))sin(t)');
zlabel('z=cos(20t)');
view([30, 40]);


C = 1;
x1 = arrayfun(@(x) (C + sin(K*x))*cos(x), t);
y1 = arrayfun(@(x) (C + sin(K*x))*sin(x), t);
z1 = arrayfun(@(x) cos(K*x), t);

figure;
plot3(x1, y1, z1);
xlabel('x=(1+sin(20t))cos(t)');
ylabel('y=(1+sin(20t))sin(t)');
zlabel('z=cos(20t)');
view([30, 40]);


C = 2;
x2 = arrayfun(@(x) (C + sin(K*x))*cos(x), t);
y2 = arrayfun(@(x) (C + sin(K*x))*sin(x), t);
z2 = arrayfun(@(x) cos(K*x), t);

figure;
plot3(x2, y2, z2);
xlabel('x=(2+sin(20t))cos(t)');
ylabel('y=(2+sin(20t))sin(t)');
zlabel('z=cos(20t)');
view([30, 40]);

C = 3;
x3 = arrayfun(@(x) (C + sin(K*x))*cos(x), t);
y3 = arrayfun(@(x) (C + sin(K*x))*sin(x), t);
z3 = arrayfun(@(x) cos(K*x), t);

figure;
plot3(x3, y3, z3);
xlabel('x=(3+sin(20t))cos(t)');
ylabel('y=(3+sin(20t))sin(t)');
zlabel('z=cos(20t)');
view([30, 40]);

C = 4;
x4 = arrayfun(@(x) (C + sin(K*x))*cos(x), t);
y4 = arrayfun(@(x) (C + sin(K*x))*sin(x), t);
z4 = arrayfun(@(x) cos(K*x), t);

figure;
plot3(x4, y4, z4);
xlabel('x=(4+sin(20t))cos(t)');
ylabel('y=(4+sin(20t))sin(t)');
zlabel('z=cos(20t)');
view([30, 40]);

C = 5;
x5 = arrayfun(@(x) (CONST + sin(K*x))*cos(x), t);
y5 = arrayfun(@(x) (CONST + sin(K*x))*sin(x), t);
z5 = arrayfun(@(x) cos(K*x), t);

figure;
plot3(x5, y5, z5);
xlabel('x=(5+sin(20t))cos(t)');
ylabel('y=(5+sin(20t))sin(t)');
zlabel('z=cos(20t)');
view([30, 40]);

---

Integral Approximation with Infinite Limits and Regular Partitions

Vladimir Kulyukin

---

Problem
  

There is a theorem in integral calculus that states that if a function f is integrable on [a, b], then its definite integral on [a, b] is equal to an infinite limit of the sum of f evaluated at regular points obtained with ever decreasing regular partitions. In other words:

Let us investigate how we can use this theorem to programmatically evaluate definite integrals.

---

Function Abstraction

We need to use functions as objects. Toward that result, we define a one-argument Function class that all other subclasses will override. Source code is here.

public class Function {
   
    public Function() {}
    public double v(double x) { return 0; }
    public double[] generateRangeInterval(double[] domain_values) {
        return null;
    }
   
}

The idea of regular partition can now be encapsulated as a subclass of the Function class (RegularPartition.java):

public class RegularPartitionSum extends Function {
   
    private double mA = 0;
    private double mB = 0;
    private Function mF = null;
   
    public RegularPartitionSum(Function f, double a, double b) {
        mA = a;
        mB = b;
        mF = f;
    }
   
    @Override
    public double v(double n) {
        final double step = (mB - mA)/n;
        final int in = (int)n;
        double sum = 0;
        for(int i = 0; i <= in; i++) {
            sum += mF.v(mA + i*step);
        }
        return step*sum;
    }
   
}

Given a function, we should be able to approximate its limit at infinity to any desirable number of steps, which we abstract into a FunctionLimit class (FunctionLimit.java).

public class FunctionLimit {

public static double limitAtInfinity(Function f, double start, double step, int num_steps) {
        double x = start;
        double y = 0;
        while ( num_steps >= 0 ) {
            y = f.v(x);
            x += step;
            num_steps--;
        }
        return y;
}

We proceed by implementing Riemann's sum (RiemannSum.java) with the method that uses FunctionLimit to estimate a definite integral with regular partitions and infinite limits. We compute the limit of the regular partition sums for a specific number of steps.

public class RiemannSum {
    
    public static double infiniteLimitOfRegularPartition(Function f, double a, double b, int num_steps) {
        Function rpF = new RegularPartitionSum(f, a, b);
        return FunctionLimit.limitAtInfinity(rpF, 1, 1, num_steps);
    }
}

---

Sample Functions

Let us define several functions and use the above implementation to evaluate their definite integrals. All of them come from James Stewart's "Calculus: Early Transendentals," 3rd edition. The name of the function class specifies the number of example and the page.

public class F_example_07_p343 extends Function {
   
    public F_example_07_p343() {}
   
    @Override
    public double v(double x) {
        return Math.sqrt(x);
    }
}

The function example_07_p343(), defined below, when placed in the main method, outputs "RS = 4.668915495400732," which is pretty close to the manual evaluation of the definite intergal.
 
public static void example_07_p343() {
        Function f = new F_example_07_p343();
        System.out.println("RS  = " + RiemannSum.infiniteLimitOfRegularPartition(f, 1, 4, 2000));
}

Here is another sample function:

public class F_ex09_p344 extends Function {
   
    public F_ex09_p344() {}
   
    @Override
    public double v(double x) {
        return x*x*x;
    }
}

We can evaluate this function in this method:

public static void ex_09_p344() {
        System.out.println("RS = " + RiemannSum.infiniteLimitOfRegularPartition(new F_ex09_p344(), 0, 5, 1000));
}

The output value is 156.56234375015595. More sample functions are implemented and tested here.

  

---

Midpoint Rule

The accuracy of our evaluations can be compared with the mid-point rule, which is a more standard method of evaluating definite integrals. Add these methods to RiemannSum and use them to evaluate the sample functions.

public static double midPointArea(Function f, double from, double upto, double step) {
        double area = 0;
        if ( from >= upto ) return 0;
        double currP = from;
        while ( currP <= upto - step ) {
            area += f.v(currP+step/2.0);
            currP += step;
        }
        return step*area;
}

public static double midPointRule(Function f, double a, double b, int n) {
        if ( b < a ) throw new IllegalArgumentException("midPointRule: b < a");
        if ( a == b ) return 0;
        final double step = (b - a)/n;
        return RiemannSum.midPointArea(f, a, b, step);
}
  
 

---

Programmatic Construction of Multi-Scale 2D Matrices for Direct and Inverse 1D Haar Wavelet Transform

Vladimir Kulyukin

---

Problem
  

Given a 1D signal of length 2^n, how can we construct the 2^n x 2^n matrices for 1D direct and inverse Haar wavelet transforms?

---

Direct Matrix

The 2^n by 2^n matrix for 2D direct Haar wavelet transform is obtained by computing 1D Haar transforms of the 2^n basis vectors. Each such vector becomes the corresponding column of the matrix. In other words, the directly transformed 0th base vector becomes the 0th column of the direct matrix, the 1st transformed base vector becomes the 1st column in the direct matrix and so on.  For example, if we want to compute the forward matrix for any signal of length 2^3, we need to convert the following 8 basis vectors:

double[] bv0 = {1, 0, 0, 0, 0, 0, 0, 0};
double[] bv1 = {0, 1, 0, 0, 0, 0, 0, 0};
double[] bv2 = {0, 0, 1, 0, 0, 0, 0, 0};
double[] bv3 = {0, 0, 0, 1, 0, 0, 0, 0};
double[] bv4 = {0, 0, 0, 0, 1, 0, 0, 0};
double[] bv5 = {0, 0, 0, 0, 0, 1, 0, 0};
double[] bv6 = {0, 0, 0, 0, 0, 0, 1, 0};
double[] bv7 = {0, 0, 0, 0, 0, 0, 0, 1};

The 1D Haar transforms of these vectors can be computed by using OneDHaar.orderedFastHaarWaveletTransform(double[] v). Let us quickly define two methods to display the results.

public static void forward_basis_vector(double[] bv) {
        OneDHaar.orderedFastHaarWaveletTransform(bv);
        OneDHaar.displaySample(bv);
}

public static void test_forward_basis_vectors() {
        double[] bv0 = {1, 0, 0, 0, 0, 0, 0, 0};
        double[] bv1 = {0, 1, 0, 0, 0, 0, 0, 0};
        double[] bv2 = {0, 0, 1, 0, 0, 0, 0, 0};
        double[] bv3 = {0, 0, 0, 1, 0, 0, 0, 0};
        double[] bv4 = {0, 0, 0, 0, 1, 0, 0, 0};
        double[] bv5 = {0, 0, 0, 0, 0, 1, 0, 0};
        double[] bv6 = {0, 0, 0, 0, 0, 0, 1, 0};
        double[] bv7 = {0, 0, 0, 0, 0, 0, 0, 1};
       
        forward_basis_vector(bv0);
        forward_basis_vector(bv1);
        forward_basis_vector(bv2);
        forward_basis_vector(bv3);
        forward_basis_vector(bv4);
        forward_basis_vector(bv5);
        forward_basis_vector(bv6);
        forward_basis_vector(bv7);
}

The output of test_forward_basis_vectors() is as follows:

0.125 0.125 0.25 0.0 0.5 0.0 0.0 0.0 
0.125 0.125 0.25 0.0 -0.5 0.0 0.0 0.0 
0.125 0.125 -0.25 0.0 0.0 0.5 0.0 0.0 
0.125 0.125 -0.25 0.0 0.0 -0.5 0.0 0.0 
0.125 -0.125 0.0 0.25 0.0 0.0 0.5 0.0 
0.125 -0.125 0.0 0.25 0.0 0.0 -0.5 0.0 
0.125 -0.125 0.0 -0.25 0.0 0.0 0.0 0.5 
0.125 -0.125 0.0 -0.25 0.0 0.0 0.0 -0.5

Now we can compute the direct or forward 2D Haar matrix for any n. 

public static double[][] computeForwardHaarTransformMatrix(int n) {
        int size = (int) Math.pow(2, n);
        double[] base_vector = new double[size];
        double[][] fhw = new double[size][size];
        for(int col_num = 0; col_num < size; col_num++) {
            for(int i = 0; i < size; i++) {
                if ( i == col_num )
                    base_vector[i] = 1;
                else
                    base_vector[i] = 0;
            }
            OneDHaar.orderedFastHaarWaveletTransform(base_vector);
            for(int row_num = 0; row_num < size; row_num++) {
                fhw[row_num][col_num] = base_vector[row_num];
            }
        }
       
        return fhw;
}

We can test our implementation with testForwardHaarTransformMatrix():

public static void testForwardHaarTransformMatrix() {
        double[][] fhtm3 = OneDHaar.computeForwardHaarTransformMatrix(3);
        double[][] fhtm4 = OneDHaar.computeForwardHaarTransformMatrix(4);
        TwoDHaar.displaySample(fhtm3, 8, 0);
        double[] sample00 = {1, 1, 1, 1, 1, 1, 1, 1};
        double[] sample01 = {1, 1, 1, 1, 1, 1, 1, 1};
        System.out.print("Input:\t"); OneDHaar.displaySample(sample00);
        OneDHaar.orderedFastHaarWaveletTransform(sample00);
        System.out.print("FHWT:\t"); OneDHaar.displaySample(sample00);
        double[] fwd_sample01 = OneDHaar.applyHaarTransformMatrix(fhtm3, sample01);
        System.out.print("MHWT:\t"); OneDHaar.displaySample(fwd_sample01);
        System.out.println();
       
        double[] sample10 = {1, 1, 1, 1, 0, 0, 0, 0};
        double[] sample11 = {1, 1, 1, 1, 0, 0, 0, 0};
        System.out.print("Input:\t"); OneDHaar.displaySample(sample10);
        OneDHaar.orderedFastHaarWaveletTransform(sample10);
        System.out.print("FHWT:\t"); OneDHaar.displaySample(sample10);
        double[] fwd_sample11 = OneDHaar.applyHaarTransformMatrix(fhtm3, sample11);
        System.out.print("MHWT:\t"); OneDHaar.displaySample(fwd_sample11);
       
        double[] sample20 = {12.0, 16.0, 27.0, 32.8, 33.5, 33.5, 39.0, 39.0,
                             40.0, 41.3, 41.3, 42.0, 43.0, 45.0, 35.5, 49.0};
        double[] sample21 = {12.0, 16.0, 27.0, 32.8, 33.5, 33.5, 39.0, 39.0,
                             40.0, 41.3, 41.3, 42.0, 43.0, 45.0, 35.5, 49.0};
        System.out.print("Input:\t"); OneDHaar.displaySample(sample20);
        OneDHaar.orderedFastHaarWaveletTransform(sample20);
        System.out.print("FHWT:\t"); OneDHaar.displaySample(sample20);
        double[] fwd_sample21 = OneDHaar.applyHaarTransformMatrix(fhtm4, sample21);
        System.out.print("MHWT:\t"); OneDHaar.displaySample(fwd_sample21);
 }

---

Inverse Matrix

The 2^n by 2^n matrix for the inverse 2D Haar wavelet transform is obtained by computing the inverse 1D Haar transforms of the 2^n basis vectors. Each such vector becomes the corresponding column of the inverse  matrix. In other words, the inversely transformed 0th base vector becomes the 0th column vector, the  1st inversely transformed base vector becomes the 1st column in the direct matrix and so on.  For example, if we want to compute the inverse matrix for any signal of length 2^3, we need to compute the inverse transforms the following 8 basis vectors:

double[] bv0 = {1, 0, 0, 0, 0, 0, 0, 0};
double[] bv1 = {0, 1, 0, 0, 0, 0, 0, 0};
double[] bv2 = {0, 0, 1, 0, 0, 0, 0, 0};
double[] bv3 = {0, 0, 0, 1, 0, 0, 0, 0};
double[] bv4 = {0, 0, 0, 0, 1, 0, 0, 0};
double[] bv5 = {0, 0, 0, 0, 0, 1, 0, 0};
double[] bv6 = {0, 0, 0, 0, 0, 0, 1, 0};
double[] bv7 = {0, 0, 0, 0, 0, 0, 0, 1};

The inverse 1D Haar transforms of these vectors can be computed by using OneDHaar.orderedFastInverseHaarWaveletTransform(double[] v). Let us quickly define two methods to display the results.

public static void inverse_basis_vector(double[] bv) {
        OneDHaar.orderedFastInverseHaarWaveletTransform(bv);
        OneDHaar.displaySample(bv);
 }
    
public static void test_inverse_basis_vectors() {
        double[] bv0 = {1, 0, 0, 0, 0, 0, 0, 0};
        double[] bv1 = {0, 1, 0, 0, 0, 0, 0, 0};
        double[] bv2 = {0, 0, 1, 0, 0, 0, 0, 0};
        double[] bv3 = {0, 0, 0, 1, 0, 0, 0, 0};
        double[] bv4 = {0, 0, 0, 0, 1, 0, 0, 0};
        double[] bv5 = {0, 0, 0, 0, 0, 1, 0, 0};
        double[] bv6 = {0, 0, 0, 0, 0, 0, 1, 0};
        double[] bv7 = {0, 0, 0, 0, 0, 0, 0, 1};
      
        inverse_basis_vector(bv0);
        inverse_basis_vector(bv1);
        inverse_basis_vector(bv2);
        inverse_basis_vector(bv3);
        inverse_basis_vector(bv4);
        inverse_basis_vector(bv5);
        inverse_basis_vector(bv6);
        inverse_basis_vector(bv7);
}

The output of running test_inverse_basis_vectors() is given below.

1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0
1.0 1.0 1.0 1.0 -1.0 -1.0 -1.0 -1.0 
1.0 1.0 -1.0 -1.0 0.0 0.0 0.0 0.0 
0.0 0.0 0.0 0.0 1.0 1.0 -1.0 -1.0 
1.0 -1.0 0.0 0.0 0.0 0.0 0.0 0.0 
0.0 0.0 1.0 -1.0 0.0 0.0 0.0 0.0 
0.0 0.0 0.0 0.0 1.0 -1.0 0.0 0.0  
0.0 0.0 0.0 0.0 0.0 0.0 1.0 -1.0 

We can now compute the inverse 2D Haar matrix for any n:

 public static double[][] computeInverseHaarTransformMatrix(int n) {
        int size = (int) Math.pow(2, n);
        double[] base_vector = new double[size];
        double[][] ihw = new double[size][size];
        for(int col_num = 0; col_num < size; col_num++) {
            for(int i = 0; i < size; i++) {
                if ( i == col_num )
                    base_vector[i] = 1;
                else
                    base_vector[i] = 0;
            }
            OneDHaar.orderedFastInverseHaarWaveletTransform(base_vector);
            for(int row_num = 0; row_num < size; row_num++) {
                ihw[row_num][col_num] = base_vector[row_num];
            }
        }
       
        return ihw;
}

We can test the computation of the inverse 2D matrix as follows:

public static void testInverseHaarTransformMatrix() {
        double[][] ihtm3 = OneDHaar.computeInverseHaarTransformMatrix(3);
        double[][] ihtm4 = OneDHaar.computeInverseHaarTransformMatrix(4);
        TwoDHaar.displaySample(ihtm3, 8, 0);
       
        double[] sample00 = {1, 1, 1, 1, 1, 1, 1, 1};
        double[] sample01 = {1, 1, 1, 1, 1, 1, 1, 1};
        OneDHaar.orderedFastHaarWaveletTransform(sample00);
        OneDHaar.orderedFastHaarWaveletTransform(sample01);
       
        System.out.print("Input:\t"); OneDHaar.displaySample(sample00);
        OneDHaar.orderedFastInverseHaarWaveletTransform(sample00);
        System.out.print("IFHWT:\t"); OneDHaar.displaySample(sample00);
        double[] isample01 = OneDHaar.applyHaarTransformMatrix(ihtm3, sample01);
        System.out.print("IMHWT:\t"); OneDHaar.displaySample(isample01);
        System.out.println();
       
        double[] sample10 = {1, 1, 1, 1, 0, 0, 0, 0};
        double[] sample11 = {1, 1, 1, 1, 0, 0, 0, 0};
        OneDHaar.orderedFastHaarWaveletTransform(sample10);
        OneDHaar.orderedFastHaarWaveletTransform(sample11);
        System.out.print("Input:\t"); OneDHaar.displaySample(sample10);
        OneDHaar.orderedFastInverseHaarWaveletTransform(sample10);
        System.out.print("IFHWT:\t"); OneDHaar.displaySample(sample10);
        double[] isample11 = OneDHaar.applyHaarTransformMatrix(ihtm3, sample11);
        System.out.print("IMHWT:\t"); OneDHaar.displaySample(isample11);
       
        double[] sample20 = {12.0, 16.0, 27.0, 32.8, 33.5, 33.5, 39.0, 39.0,
                             40.0, 41.3, 41.3, 42.0, 43.0, 45.0, 35.5, 49.0};
        double[] sample21 = {12.0, 16.0, 27.0, 32.8, 33.5, 33.5, 39.0, 39.0,
                             40.0, 41.3, 41.3, 42.0, 43.0, 45.0, 35.5, 49.0};
        OneDHaar.orderedFastHaarWaveletTransform(sample20);
        OneDHaar.orderedFastHaarWaveletTransform(sample21);
        System.out.print("Input:\t"); OneDHaar.displaySample(sample20);
        OneDHaar.orderedFastInverseHaarWaveletTransform(sample20);
        System.out.print("IFHWT:\t"); OneDHaar.displaySample(sample20);
        double[] isample21 = OneDHaar.applyHaarTransformMatrix(ihtm4, sample21);
        System.out.print("IMHWT:\t"); OneDHaar.displaySample(isample21);
  }
  

---

Matrix Multiplication

Here is how the direct or inverse matrices can be applied to the transform or invert the 1D signals:

public static double[] applyHaarTransformMatrix(double[][] htm, double[] v) {
        int num_rows = htm.length;
        if ( num_rows < 1 ) return null;
        int num_cols = htm[0].length;
        if ( num_cols < 1 ) return null;
        if ( num_rows != num_cols ) return null;
        if ( num_rows != v.length ) return null;
        double[] inversed_v = new double[num_cols];
        double dot_product = 0;
        for(int row = 0; row < num_rows; row++) {
            dot_product = 0;
            for(int col = 0; col < num_cols; col++) {
                dot_product += htm[row][col]*v[col];
            }
            inversed_v[row] = dot_product;
        }
        return inversed_v;
    }
  

---

2D FFT: Binarization of 2D FFT Spectra of Image Waves of Different Frequencies and Orientations

Vladimir Kulyukin

---

Problem
  

This is a list of posts on the binarization of the 2D FFT frequency spectra of wave images of different frequencies and orientations. I could not put all images into a single post and had to break it into four different posts each of which contains an example of how to binarize the 2D FFT of wave images.

---

List of Blog Posts

1. Binarization of the 2D FFT frequency spectrum of an image wave moving north-east with a medium frequency
2. Binarization of the 2D FFT frequency spectrum of an image wave moving south-east with a greater frequency
3. Binarization of the 2D FFT frequency spectrum of an image wave moving south-east with an even greater frequency
4. Binarization of the 2D FFT frequency spectrum of a image wave moving east with a high frequency
   

---