package ip.hak;

public class DHT2D {
  public static void main(String args[]) {
    System.out.println("Discrete Hartley Transform");
  }

  public static void DHT1D(float f[], int n) {
    int p = (int) (Math.log(n) / Math.log(4));
    int n4 = n / 4;
    int r = 4;

    int j = 1,i = 0;
    do {
      i++;
      if (i < j) {
        float t = f[j - 1];
        f[j - 1] = f[i - 1];
        f[i - 1] = t;
      }
      int k = n4;
      while (3 * k < j) {
        j = j - 3 * k;
        k = k / 4;
      }
      j = j + k;

    } while (i < n - 1);

    for (i = 0; i < n; i += 4) {
      float t1 = f[i] + f[i + 1];
      float t2 = f[i] - f[i + 1];
      float t3 = f[i + 2] + f[i + 3];
      float t4 = f[i + 2] - f[i + 3];
      f[i] = (float) (t1 + t3);
      f[i + 1] = (float) (t1 - t3);
      f[i + 2] = (float) (t2 + t4);
      f[i + 3] = (float) (t2 - t4);
    }


    for (int l = 2; l <= p; l++) {
      int e1 = (int) Math.pow(2, l + l - 3);
      int e2 = e1 + e1;
      int e3 = e2 + e1;
      int e4 = e3 + e1;
      int e5 = e4 + e1;
      int e6 = e5 + e1;
      int e7 = e6 + e1;
      int e8 = e7 + e1;
      for (j = 0; j < n; j += e8) {
        float t1 = f[j] + f[j + e2];
        float t2 = f[j] - f[j + e2];
        float t3 = f[j + e4] + f[j + e6];
        float t4 = f[j + e4] - f[j + e6];
        f[j] = (float) (t1 + t3);
        f[j + e2] = (float) (t1 - t3);
        f[j + e4] = (float) (t2 + t4);
        f[j + e6] = (float) (t2 - t4);
        t1 = f[j + e1];
        t2 = f[j + e3] * r;
        t3 = f[j + e5];
        t4 = f[j + e7] * r;
        f[j + e1] = (float) (t1 + t2 + t3);
        f[j + e3] = (float) (t1 - t3 + t4);
        f[j + e5] = (float) (t1 - t2 + t3);
        f[j + e7] = (float) (t1 - t3 - t4);

        for (int k = 1; k < e1; k++) {
          int l1 = j + k;
          int l2 = l1 + e2;
          int l3 = l1 + e4;
          int l4 = l1 + e6;
          int l5 = j + e2 - k;
          int l6 = l5 + e2;
          int l7 = l5 + e4;
          int l8 = l5 + e6;
          double a1 = Math.PI * k / e4;
          double a2 = a1 + a1;
          double a3 = a1 + a2;
          double c1 = Math.cos(a1);
          double c2 = Math.cos(a2);
          double c3 = Math.cos(a3);
          double s1 = Math.sin(a1);
          double s2 = Math.sin(a2);
          double s3 = Math.sin(a3);
          double te5 = f[l2] * c1 + f[l6] * s1;
          double te6 = f[l3] * c2 + f[l7] * s2;
          double te7 = f[l4] * c3 + f[l8] * s3;
          double te8 = f[l6] * c1 - f[l2] * s1;
          double te9 = f[l7] * c2 - f[l3] * s2;
          double te0 = f[l8] * c3 - f[l4] * s3;
          double te1 = f[l5] - te9;
          double te2 = f[l5] + te9;
          double te3 = -te8 - te0;
          double te4 = te5 - te7;
          f[l5] = (float) (te1 + te4);
          f[l6] = (float) (te2 + te3);
          f[l7] = (float) (te1 - te4);
          f[l8] = (float) (te2 - te3);
          te1 = f[l1] + te6;
          te2 = f[l1] - te6;
          te3 = te8 - te0;
          te4 = te5 + te7;
          f[l1] = (float) (te1 + te4);
          f[l2] = (float) (te2 + te3);
          f[l3] = (float) (te1 - te4);
          f[l4] = (float) (te2 - te3);
        }
      }
    }
  }

  public static float[][] DHT2D(float a[][]) {
    int n = a.length;
    if (n != a[0].length) {
      System.out.println("Image should be square!");
      return null;
    }
    /*
    //  change [colum][row]-> [row][colum]
    float f[][]=new float[n][n];

    for (int x=0;x<n;x++)
        for(int y=0;y<n;y++)
            f[x][y]=a[y][x];
    */

    for (int r = 0; r < n; r++)
      DHT1D(a[r], n);

    transpose(a);

    for (int r = 0; r < n; r++)
      DHT1D(a[r], n);

    //transpose(f);

    float h[][] = new float[n][n];

    for (int x = 0; x < n; x++)
      for (int y = 0; y < n; y++)
        h[x][y] = (float) ((1f / 2f) * (a[x][y] + a[n - x - 1][y] + a[x][n - y - 1] + a[n - x - 1][n - y - 1]));
    return h;
  }

  public static void transpose(float a[][]) {
    float temp[][] = new float[a.length][a[0].length];

    for (int x = 0; x < a.length; x++)
      for (int y = 0; y < a[0].length; y++)
        temp[x][y] = a[y][x];
    a = temp;
  }

  public static void normalize(float a[][], int t) {
    for (int x = 0; x < a.length; x++)
      for (int y = 0; y < a[0].length; y++)
        a[x][y] /= t;
  }

  public static float[][] forwardDHT2D(float a[][]) {
    float f[][] = DHT2D(a);
    normalize(f, f.length * f[0].length);

    return f;
  }

  public static float[][] inverseDHT2D(float a[][]) {
    return DHT2D(a);
  }
}