burg.c

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <time.h>

#define szof(a) ((int)sizeof(a))
#define cntof(a) ((int)(sizeof(a) / sizeof((a)[0])))

static void
burg(double *coeffs, int m, const double *x, int N) {
  /* init Ak */
  double Ak[m+1];
  memset(Ak, 0, sizeof(Ak));
  Ak[0] = 1.0;

  /* init f and b */
  double f[N], b[N];
  memcpy(f, x, sizeof(*x) * N);
  memcpy(b, x, sizeof(*x) * N);

  /* setup Dk */
  double Dk = 0.0;
  for (int j = 0; j <= N; j++){
    Dk += 2.0 * f[j] * f[j];
  }
  Dk -= f[0] * f[0] + b[N] * b[N];

  /* burg recursion */
  for (int k = 0; k < m; k++) {
    double mu = 0.0;
    for (int n = 0; n <= N - k - 1; n++ ) {
      mu += f[ n + k + 1 ] * b[ n ];
    }
    mu *= -2.0 / Dk;
    for (int n = 0; n <= (k + 1) / 2; n++) {
      double t1 = Ak[n] + mu * Ak[k + 1 - n];
      double t2 = Ak[k + 1 - n] + mu * Ak[n];
      Ak[n] = t1;
      Ak[k + 1 - n] = t2;
    }
    for (int n = 0; n <= N - k - 1; n++ ) {
      double t1 = f[n + k + 1] + mu * b[n];
      double t2 = b[n] + mu * f[n + k + 1];
      f[n + k + 1] = t1;
      b[n] = t2;
    }
    Dk = (1.0 - mu * mu) * Dk - f[k + 1] * f[k + 1] - b[N - k - 1] * b[N - k - 1];
  }
  for (int i = 1; i <= m; ++i) {
    coeffs[i-1] = Ak[i];
  }
}
extern int
main(void) {
  // CREATE DATA TO APPROXIMATE
  double original[128] = {0};
  for (int i = 0; i < cntof(original); i++) {
    original[i] = cosf(i * 0.01) + 0.75 *cosf(i * 0.03) + 0.5 *cosf(i * 0.05) + 0.25 *cosf(i * 0.11);
  }

  // GET LINEAR PREDICTION COEFFICIENTS
  double coeffs[8] = {0.0};
  burg(coeffs, cntof(coeffs), original, cntof(original));

  // LINEAR PREDICT DATA
  double predicted[128];
  for (int i = 0; i < cntof(predicted); ++i) {
    predicted[i] = original[i];
  }
  for (int i = cntof(coeffs); i < cntof(original); i++ ) {
    predicted[i] = 0.0;
    for (int j = 0; j < cntof(coeffs); j++ ) {
      predicted[i] -= coeffs[j] * original[i - 1 - j];
    }
  }
  // CALCULATE AND DISPLAY ERROR
  double err = 0.0;
  for (int i = cntof(coeffs); i < cntof(predicted); i++ ) {
    printf( "Index: %.2d / Original: %.6f / Predicted: %.6f\n", i, original[ i ], predicted[ i ] );
    double delta = predicted[ i ] - original[ i ];
    err += delta * delta;
  }
  printf("Burg Approximation Error: %f\n", err);
  return 0;
}

gistfile1.txt

https://encode.su/threads/3015-Oodle-Lossless-Image-compression
https://aras-p.info/blog/2023/02/01/Float-Compression-3-Filters/
https://github.com/veluca93/fpnge/blob/926df95bce7bd1affaa0163572ac6f0ae692eb95/fpnge.cc#L619
https://www.investopedia.com/terms/a/autocorrelation.asp
https://ruby0x1.github.io/machinery_blog_archive/post/data-structures-part-2-indices/index.html

http://www.emptyloop.com/technotes/
https://github.com/RhysU/ar/blob/master/collomb2009.cpp
https://github.com/r-lyeh-archived/burg/blob/master/burg.hpp
https://www.researchgate.net/profile/Aleksej-Avramovic/publication/251961072_Gradient_edge_detection_predictor_for_image_lossless_compression/links/0c960539cc73768333000000/Gradient-edge-detection-predictor-for-image-lossless-compression.pdf

lossless_image_compression_filter.c

#include <stdio.h>
#include <stdlib.h>
#include <string.h>

static void
auto_corr(double *r, const unsigned char *a, int n, int lag) {
  double mean = 0.0;
  for (int i = 0; i < n; i++) {
    mean += a[i];
  }
  mean /= n;

  for (int k = 0; k <= lag; ++k) {
    double sk = 0.0;
    for (int i = k; i < n; ++i) {
      sk += (a[i] - mean) * (a[i - k] - mean);
    }
    r[k] = sk / n;
  }
  for (int i = lag; i >= 0; --i) {
    r[i] = r[i] / r[0];
  }
}
static void
levinson_durbin(double *a, const double *r, int p) {
  double sum = 0.0;
  double tmp[p + 1];
  for (int i = 0; i < p; ++i) {
    a[i] = tmp[i] = 0.0;
  }
  a[0] = r[0];
  for (int i = 1; i <= p; ++i) {
    sum = r[i];
    for (int j = 0; j < i; ++j) {
      sum -= tmp[j] * r[i - j - 1];
    }
    tmp[i] = sum / a[i - 1];
    for (int j = 0; j < i; ++j) {
      a[j] -= tmp[i] * a[i - j - 1];
    }
    a[i] = tmp[i];
  }
}
static void
innovations_algo(double* c, const double* rho, int order) {
  for(int i = 0; i <= order; ++i) {
    c[i] = 0.0;
  }
  c[1] = rho[1]/rho[0];
  for(int k = 2; k <= order; ++k) {
    double sum = 0.0;
    for(int j = 1; j <= k-1; ++j) {
      sum += rho[j]*c[j]*c[k-j];
    }
    c[k] = (rho[k] - sum) / rho[0];
  }
}
static int
pred(unsigned char* sig, int n, const int* coef, int order) {
  long long v = 0;
  for (int i = 0; i < order; ++i) {
    v += ((long long)coef[i] * (long long)(sig[n - i - 1] << 16)) >> 16;
  }
  return v >> 16;
}

#define MAX_ORDER 3 // Define the maximum order of the predictor

extern int
main(void) {
   // TODO: Load the image data into a buffer
  unsigned char image_data[] = {253, 252, 251, 245, 234, 250, 251, 242, 232, 222, 212, 252, 254, 253, 252, 253, 254, 253, 251, 253}; // Some sample signal data
  int image_size = sizeof(image_data)/sizeof(image_data[0]);

  // calculate auto correlation
  double ac[MAX_ORDER+1] = {0};
  auto_corr(ac, image_data, image_size, MAX_ORDER);
  for (int i = 0; i <= MAX_ORDER; ++i) {
    printf("%f\n", ac[i]);
  }
  printf("\n");

  // calculate coefficents
  printf("LevinsonDurbin:\n");
  int coef[MAX_ORDER+1] = {0};
  double fcoef[MAX_ORDER+1] = {0};
  levinson_durbin(fcoef, ac, MAX_ORDER);
  for (int i = 0; i <= MAX_ORDER; ++i) {
    coef[i] = fcoef[i] * (1 << 16);
    printf("%d [%f]\n", coef[i], fcoef[i]);
  }
  printf("\n");

  printf("Innovations:\n");
  int coef2[MAX_ORDER+1] = {0};
  double fcoef2[MAX_ORDER+1] = {0};
  innovations_algo(fcoef2, ac, MAX_ORDER);
  for (int i = 0; i <= MAX_ORDER; ++i) {
    coef2[i] = fcoef2[i] * (1 << 16);
    printf("%d [%f]\n", coef2[i], fcoef2[i]);
  }
  printf("\n");

  // encode delta between prediction and reality
  char d[image_size];
  memset(d, 0, sizeof(d));
  for( int i = 0; i < MAX_ORDER; ++i) {
    d[i] = image_data[i];
  }
  for (int i = MAX_ORDER; i < image_size; ++i) {
    unsigned char v = pred(image_data, i, coef, MAX_ORDER);
    d[i] = image_data[i] - v;
    printf("[%d] = %d (%d:%d)\n", i, d[i], v, image_data[i]); 
  }

#if 0
  double img_sum = 0.0f;
  double dt_sum = 0.0f;
  double img_min = 256.0f;
  double img_max = 0.0f;
  double dt_min = 256.0f;
  double dt_max = 0.0f;

  for (int i = 0; i < image_size; ++i) {
    img_sum += abs(image_data[i]);
    img_min = image_data[i] < img_min ? image_data[i] : img_min;
    img_max = image_data[i] > img_max ? image_data[i] : img_max;
    dt_min = d[i] < dt_min ? d[i] : dt_min;
    dt_max = d[i] > dt_max ? d[i] : dt_max;
    dt_sum += d[i];
  }
  printf("img mean: %f[%f, %f, %f]:\nd mean: %f[%f %f %f]\n\n", img_sum / image_size, img_min, img_max, img_max - img_min, dt_sum / image_size, dt_min, dt_max, dt_max - dt_min);
#endif

  // decode orignal image from prediction and reality
  unsigned char r[image_size];
  memcpy(r, d, sizeof(d));
  for (int i = MAX_ORDER; i < image_size; ++i) {
    int v = pred(r, i, coef, MAX_ORDER);
    int av = v < 0 ? -v : v;
    r[i] = d[i] + av;
  }
  for(int i = 0; i < image_size; ++i){
    printf("[%d] = %d\n", i, r[i]);
  }
  return 0;
}

optimize.c

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <time.h>

#define szof(a) ((int)sizeof(a))
#define cntof(a) ((int)(sizeof(a) / sizeof((a)[0])))

#define coeff(x,y,c) ((c)*9+(y)*3+(x))
#define bias(c) (((c)+1)*9-1)
#define pix(x,y,w,c) (((y)*(w)+(x))*4+(c))

static double
pred_pix_chl(const unsigned char *img, int w, int h, int x, int y, int chl, const double *coef) {
  double pred = 0;
  for (int dy = 0; dy < 3; ++dy) {
    for (int dx = 0; dx < 3-(dy==2)*2; ++dx) {
      int nx = x - 2 + dx, ny = y - 2 + dy;
      if (nx >= 0 && nx < w && ny >= 0 && ny < h) {
        pred += coef[coeff(dx,dy,chl)] * img[pix(nx,ny,w,chl)];
      } 
    }
  }
  return pred + coef[bias(chl)];
}
static double
pred_err(const unsigned char *img, int w, int h, int x, int y, int chl, const double *coef) {
  double pred = pred_pix_chl(img, w, h, x, y, chl, coef);
  return (double)img[pix(x,y,w,chl)] - pred;
}
static double
opt_coef_sgd(double *coef, const unsigned char *img, int w, int h, int x, int y, 
        double learn_rate) {
  double tot_err = 0.0;
  for (int dy = 0; dy < 3; ++dy) {
    for (int dx = 0; dx < 3-(dy==2)*2; ++dx) {
      int nx = x - 2 + dx, ny = y - 2 + dy;
      if (nx < 0 || nx >= w ||  ny < 0 || ny >= h) {
        continue;
      }
      for (int c = 0; c < 4; ++c) {
        double err = pred_err(img, w, h, nx, ny, c, coef);
        tot_err += err * err;
        for (int k = 0; k < 8; ++k) {
          coef[c*9+k] -= learn_rate * err * (double)img[pix(nx,ny,w,c)];
        }
        coef[bias(c)] -= learn_rate * err;
      }
    }
  }
  return tot_err;
}
static double
opt_coef_adam(double *coef, const unsigned char *img, int w, int h, int x, int y, 
        double learn_rate, double beta1, double beta2,
        double beta1p, double beta2p, double epsilon) {
  /* ADAM: A METHOD FOR STOCHASTIC OPTIMIZATION: https://arxiv.org/pdf/1412.6980.pdf */
  double m[8*4] = {0.0}; // Initialize momentum estimates
  double v[8*4] = {0.0}; // Initialize variance estimates
  double tot_err = 0.0;
  for (int dy = 0; dy < 3; ++dy) {
    for (int dx = 0; dx < 3-(dy==2)*2; ++dx) {
      int nx = x - 2 + dx, ny = y - 2 + dy;
      if (nx < 0 || nx >= w ||  ny < 0 || ny >= h) {
        continue;
      }
      for (int c = 0; c < 4; ++c) {
        double err = pred_err(img, w, h, nx, ny, c, coef);
        tot_err += err * err;
        // Calculate gradients for each coefficient
        for (int k = 0; k < 8; ++k) {
          double gradient = err * (double)img[pix(nx,ny,w,c)];
          /* update momentum estimate */
          m[c*8+k] = beta1 * m[c*8+k] + (1 - beta1) * gradient;
          /* update variance estimate with bias correction */
          v[c*8+k] = beta2 * v[c*8+k] + (1 - beta2) * gradient * gradient;
          /* Update coefficient using bias-corrected momentum and variance */
          double m_hat = m[c*8+k] / (1 - beta1p);
          double v_hat = v[c*8+k] / (1 - beta2p);
          coef[c*9+k] -= learn_rate * m_hat / (sqrt(v_hat) + epsilon);
        }
        coef[bias(c)] -= learn_rate * err;
      }
    }
  }
  return tot_err;
}
extern int
main(int argc, char **argv) {
  srand(time(NULL));

  int w = argc;
  int h = argc;
  unsigned char *img = argv[0];

  const double learn_rate = 0.001;  /* The step length used when following the negative gradient */
  const double beta1 = 0.9;         /* The exponential decay rate for the 1st moment estimates */
  const double beta2 = 0.99;        /* The exponential decay rate for the 2nd moment estimates */
  const double epsilon = 1e-8;      /* A small floating point value to avoid zero denominator */
  const double tar_err = 0.001;
  const int num_it = 32;
  double alpha = 0.01;

  // Iterate over each pixel
  double coef[8*4+4];
  for (int i =  0; i < cntof(coef); ++i) {
    coef[i] = ((double) rand() / (RAND_MAX)) *  0.01;
  }
  for (int y = 0; y < h; ++y) {
    for (int x = 0; x < w; ++x) {
      float beta1p = beta1, beta2p = beta2;
      for (int i = 0; i < num_it; ++i) {
        // Iterate for a fixed number of times or until convergence
        double err = opt_coef_adam(coef, img, w, h, x, y, learn_rate, beta1, beta2, beta1p, beta2p, epsilon);
        if (err < tar_err) {
          break;
        }
        beta1p *= beta1;
        beta2p *= beta2;
      }
      for (int i = 0; i < num_it; ++i) {
        double err = opt_coef_sgd(coef, img, w, h, x, y, alpha);
        if (err < tar_err) {
          break;
        }
      }
    }
  }
  return 0;
}