--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/Fotopoulos/common.c	Sun Aug 12 13:14:34 2007 +0200
@@ -0,0 +1,571 @@
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <math.h>
+#include <float.h>
+#include "pgm.h"
+#include "common.h"
+
+#define IA      16807
+#define IM      2147483647
+#define AM      (1.0/IM)
+#define IQ      127773
+#define IR      2836
+#define MASK    123459876
+
+static int NN = 0;
+static int m = 0;
+static double two_over_N = 0;
+static double root2_over_rootN = 0;
+static double *C = NULL;
+
+static gray maxval;
+static int format;
+
+void open_image(FILE *in, int *width, int *height)
+{
+
+  pgm_readpgminit(in, width, height, &maxval, &format);
+}
+
+void load_image(int **im, FILE *in, int width, int height)
+{
+  int col, row;
+  gray *rowbuf;
+
+  rowbuf = malloc(sizeof(gray) * width);
+
+  for (row = 0; row < height; row++) {
+    pgm_readpgmrow(in, rowbuf, width, maxval, format);
+    for (col = 0; col < width; col++)
+      im[row][col] = rowbuf[col];
+  }
+
+  free(rowbuf);
+}
+
+void save_image(int **im, FILE *out, int width, int height)
+{
+  int col, row;
+  gray *rowbuf;
+
+  pgm_writepgminit(out, width, height, 255, 0);
+
+  rowbuf = malloc(sizeof(gray) * width);
+
+  for (row = 0; row < height; row++) {
+    for (col = 0; col < width; col++)
+      rowbuf[col] = im[row][col];
+    pgm_writepgmrow(out, rowbuf, width, 255, 0);
+  }
+
+  free(rowbuf);
+}
+
+int ** imatrix(int nrows, int ncols)
+{
+  int **m;
+  int i, j;
+  m = (int **) malloc (nrows * sizeof(int *));
+  for (i = 0; i < nrows; i++) {
+    m[i] = (int *) malloc (ncols * sizeof(int));
+    if (!m[i]) fprintf(stderr, "\nIt's not working");
+  }
+  for (i = 0; i < nrows; i++)
+    for (j = 0; j < ncols; j++)
+      m[i][j] = 0;
+  return m;
+}
+
+void freematrix(int **I, int rows)
+{
+  int k;
+  for (k = 0; k < rows; k++)
+    free (I[k]);
+}
+
+float ran0(long int *idum)
+{
+  long int k;
+  float ans;
+  *idum ^= MASK;
+  k = (*idum) / IQ;
+  *idum = IA * (*idum - k * IQ) - IR * k;
+  if (*idum < 0) *idum += IM;
+  ans = AM * (*idum);
+  *idum ^= MASK;
+  return ans;
+}
+
+float gasdev(long int *idum)
+{
+
+  float v1;
+  v1 = (float) sqrt( -2.0 * log(ran0(idum))) * cos(2 * PI * ran0(idum));
+  return v1;
+}
+
+
+void put_image_from_int_2_double(int **i, double *f, int N)
+{
+  int l, j, k;
+  k = 0;
+  for (l = 0; l < N; l++)
+    for (j = 0; j < N; j++)
+      f[k++] = (double) i[l][j];
+}
+
+void put_image_from_double_2_int(double *f, int **i, int N)
+{
+  int l, j, k;
+  k = 0;
+  for (l = 0; l < N; l++)
+    for (j = 0; j < N; j++)
+      i[l][j] = (int) f[k++];
+}
+
+
+void bitrev(double *f, int len)
+{
+  int i, j, m, halflen;
+  double temp;
+
+  if (len <= 2) return ;  /* No action necessary if n=1 or n=2 */
+  halflen = len >> 1;
+  j = 1;
+  for (i = 1; i <= len; i++) {
+    if (i < j) {
+      temp = f[j - 1];
+      f[j - 1] = f[i - 1];
+      f[i - 1] = temp;
+    }
+    m = halflen;
+    while (j > m) {
+      j = j - m;
+      m = (m + 1) >> 1;
+    }
+    j = j + m;
+  }
+}
+
+void inv_sums(double *f)
+{
+  int ii, stepsize, stage, curptr, nthreads, thread, step, nsteps;
+
+  for (stage = 1; stage <= m - 1; stage++) {
+    nthreads = 1 << (stage - 1);
+    stepsize = nthreads << 1;
+    nsteps = (1 << (m - stage)) - 1;
+    for (thread = 1; thread <= nthreads; thread++) {
+      curptr = NN - thread;
+      for (step = 1; step <= nsteps; step++) {
+        f[curptr] += f[curptr - stepsize];
+        curptr -= stepsize;
+      }
+    }
+  }
+}
+
+//--------------------------------------------------------
+//Foreign code - FCT from Bath Univerity
+//--------------------------------------------------------
+void rarrwrt(double f[], int n)
+{
+  int i;
+
+  for (i = 0; i <= n - 1; i++) {
+    fprintf(stderr, "%4d : %f\n", i, f[i]);
+  }
+}
+
+/* fast DCT based on IEEE signal proc, 1992 #8, yugoslavian authors. */
+
+void fwd_sums(double *f)
+{
+  int ii, stepsize, stage, curptr, nthreads, thread, step, nsteps;
+
+  for (stage = m - 1; stage >= 1; stage--) {
+    nthreads = 1 << (stage - 1);
+    stepsize = nthreads << 1;
+    nsteps = (1 << (m - stage)) - 1;
+    for (thread = 1; thread <= nthreads; thread++) {
+      curptr = nthreads + thread - 1;
+      for (step = 1; step <= nsteps; step++) {
+        f[curptr] += f[curptr + stepsize];
+        curptr += stepsize;
+      }
+    }
+  }
+}
+
+void scramble(double *f, int len)
+{
+  double temp;
+  int i, ii1, ii2, halflen, qtrlen;
+
+  halflen = len >> 1;
+  qtrlen = halflen >> 1;
+  bitrev(f, len);
+  bitrev(&f[0], halflen);
+  bitrev(&f[halflen], halflen);
+  ii1 = len - 1;
+  ii2 = halflen;
+  for (i = 0; i <= qtrlen - 1; i++) {
+    temp = f[ii1];
+    f[ii1] = f[ii2];
+    f[ii2] = temp;
+    ii1--;
+    ii2++;
+  }
+}
+
+void unscramble(double *f, int len)
+{
+  double temp;
+  int i, ii1, ii2, halflen, qtrlen;
+
+  halflen = len >> 1;
+  qtrlen = halflen >> 1;
+  ii1 = len - 1;
+  ii2 = halflen;
+  for (i = 0; i <= qtrlen - 1; i++) {
+    temp = f[ii1];
+    f[ii1] = f[ii2];
+    f[ii2] = temp;
+    ii1--;
+    ii2++;
+  }
+  bitrev(&f[0], halflen);
+  bitrev(&f[halflen], halflen);
+  bitrev(f, len);
+}
+
+void initcosarray(int length)
+{
+  int i, group, base, item, nitems, halfN;
+  double factor;
+
+  m = -1;
+  do {
+    m++;
+    NN = 1 << m;
+    if (NN > length) {
+      fprintf(stderr, "ERROR in FCT-- length %d not a power of 2\n", length);
+      exit(1);
+    }
+  } while (NN < length);
+  if (C != NULL) free(C);
+  C = (double *)calloc(NN, sizeof(double));
+  if (C == NULL) {
+    fprintf(stderr, "Unable to allocate C array\n");
+    exit(1);
+  }
+  halfN = NN / 2;
+  two_over_N = 2.0 / (double)NN;
+  root2_over_rootN = sqrt(2.0 / (double)NN);
+  for (i = 0; i <= halfN - 1; i++) C[halfN + i] = 4 * i + 1;
+  for (group = 1; group <= m - 1; group++) {
+    base = 1 << (group - 1);
+    nitems = base;
+    factor = 1.0 * (1 << (m - group));
+    for (item = 1; item <= nitems; item++) C[base + item - 1] = factor * C[halfN + item - 1];
+  }
+
+  //printf("before taking cos, C array =\n"); rarrwrt(C,N);
+  for (i = 1; i <= NN - 1; i++) C[i] = 1.0 / (2.0 * cos(C[i] * PI / (2.0 * NN)));
+  //printf("After taking cos, Carray = \n"); rarrwrt(C,N);
+}
+
+
+void inv_butterflies(double *f)
+{
+  int stage, ii1, ii2, butterfly, ngroups, group, wingspan, increment, baseptr;
+  double Cfac, T;
+
+  for (stage = 1; stage <= m; stage++) {
+    ngroups = 1 << (m - stage);
+    wingspan = 1 << (stage - 1);
+    increment = wingspan << 1;
+    for (butterfly = 1; butterfly <= wingspan; butterfly++) {
+      Cfac = C[wingspan + butterfly - 1];
+      baseptr = 0;
+      for (group = 1; group <= ngroups; group++) {
+        ii1 = baseptr + butterfly - 1;
+        ii2 = ii1 + wingspan;
+        T = Cfac * f[ii2];
+        f[ii2] = f[ii1]-T;
+        f[ii1] = f[ii1] + T;
+        baseptr += increment;
+      }
+    }
+  }
+}
+
+void fwd_butterflies(double *f)
+{
+  int stage, ii1, ii2, butterfly, ngroups, group, wingspan, increment, baseptr;
+  double Cfac, T;
+
+  for (stage = m; stage >= 1; stage--) {
+    ngroups = 1 << (m - stage);
+    wingspan = 1 << (stage - 1);
+    increment = wingspan << 1;
+    for (butterfly = 1; butterfly <= wingspan; butterfly++) {
+      Cfac = C[wingspan + butterfly - 1];
+      baseptr = 0;
+      for (group = 1; group <= ngroups; group++) {
+        ii1 = baseptr + butterfly - 1;
+        ii2 = ii1 + wingspan;
+        T = f[ii2];
+        f[ii2] = Cfac * (f[ii1]-T);
+        f[ii1] = f[ii1] + T;
+        baseptr += increment;
+      }
+    }
+  }
+}
+
+void ifct_noscale(double *f, int length)
+{
+  if (length != NN) initcosarray(length);
+  f[0] *= INVROOT2;
+  inv_sums(f);
+  bitrev(f, NN);
+  inv_butterflies(f);
+  unscramble(f, NN);
+}
+
+void fct_noscale(double *f, int length)
+{
+  if (length != NN) initcosarray(length);
+  scramble(f, NN);
+  fwd_butterflies(f);
+  bitrev(f, NN);
+  fwd_sums(f);
+  f[0] *= INVROOT2;
+}
+
+void ifct_defn_scaling(double *f, int length)
+{
+  ifct_noscale(f, length);
+}
+
+void fct_defn_scaling(double *f, int length)
+{
+  int i;
+
+  fct_noscale(f, length);
+  for (i = 0; i <= NN - 1; i++) f[i] *= two_over_N;
+}
+
+void ifct(double *f, int length)
+{
+  /* CALL THIS FOR INVERSE 1D DCT DON-MONRO PREFERRED SCALING */
+  int i;
+
+  if (length != NN) initcosarray(length);   /* BGS patch June 1997 */
+  for (i = 0; i <= NN - 1; i++) f[i] *= root2_over_rootN;
+  ifct_noscale(f, length);
+}
+
+void fct(double *f, int length)
+{
+  /* CALL THIS FOR FORWARD 1D DCT DON-MONRO PREFERRED SCALING */
+  int i;
+
+  fct_noscale(f, length);
+  for (i = 0; i <= NN - 1; i++) f[i] *= root2_over_rootN;
+}
+
+/****************************************************************
+    2D FAST DCT SECTION
+****************************************************************/
+
+static double *g = NULL;
+static double two_over_sqrtncolsnrows = 0.0;
+static int ncolsvalue = 0;
+static int nrowsvalue = 0;
+
+void initfct2d(int nrows, int ncols)
+{
+  if ((nrows <= 0) || (ncols < 0)) {
+    fprintf(stderr, "FCT2D -- ncols=%d or nrows=%d is <=0\n", nrows, ncols);
+    exit(1);
+  }
+  if (g != NULL) free(g);
+  g = (double *)calloc(nrows, sizeof(double));
+  if (g == NULL) {
+    fprintf(stderr, "FCT2D -- Unable to allocate g array\n");
+    exit(1);
+  }
+  ncolsvalue = ncols;
+  nrowsvalue = nrows;
+  two_over_sqrtncolsnrows = 2.0 / sqrt(ncols * 1.0 * nrows);
+}
+
+void fct2d(double f[], int nrows, int ncols)
+/* CALL THIS FOR FORWARD 2d DCT DON-MONRO PREFERRED SCALING */
+{
+  int u, v;
+
+  if ((ncols != ncolsvalue) || (nrows != nrowsvalue)){
+    initfct2d(nrows, ncols);
+  }
+  for (u = 0; u <= nrows - 1; u++){
+    fct_noscale(&f[u*ncols], ncols);
+  }
+  for (v = 0; v <= ncols - 1; v++){
+    for (u = 0; u <= nrows - 1; u++) {
+      g[u] = f[u * ncols + v];
+    }
+    fct_noscale(g, nrows);
+    for (u = 0; u <= nrows - 1; u++) {
+      f[u*ncols + v] = g[u] * two_over_sqrtncolsnrows;
+    }
+  }
+}
+
+void ifct2d(double f[], int nrows, int ncols)
+/* CALL THIS FOR INVERSE 2d DCT DON-MONRO PREFERRED SCALING */
+{
+  int u, v;
+
+  if ((ncols != ncolsvalue) || (nrows != nrowsvalue)){
+    initfct2d(nrows, ncols);
+  }
+  for (u = 0; u <= nrows - 1; u++){
+    ifct_noscale(&f[u*ncols], ncols);
+  }
+  for (v = 0; v <= ncols - 1; v++){
+    for (u = 0; u <= nrows - 1; u++) {
+      g[u] = f[u * ncols + v];
+    }
+    ifct_noscale(g, nrows);
+    for (u = 0; u <= nrows - 1; u++) {
+      f[u*ncols + v] = g[u] * two_over_sqrtncolsnrows;
+    }
+  }
+}
+
+void matmul(double **a, double **b, double **r, int N)
+{
+  int i, j, k;
+
+  for (i = 0; i < N; i++)
+    for (j = 0; j < N; j++) {
+      r[i][j] = 0.0;
+      for (k = 0; k < N; k++)
+        r[i][j] += a[i][k] * b[k][j];
+    }
+}
+
+void hartley(double **in, double **out, int N)
+{
+  int k, n;
+  double **h;
+
+  h = dmatrix(N, N);
+  //Building up the transformation matrix
+  for (k = 0; k < N; k++)
+    for (n = 0; n < N; n++)
+      h[k][n] = (cos(2 * PI * k * n / N) + sin(2 * PI * k * n / N)) / sqrt(N);
+
+  // Now we have to multiply the input with the transformation matrix
+  matmul(h, in, out, N);
+  freematrix_d(h, N);
+  return ;
+}
+
+double ** dmatrix(int nrows, int ncols)
+{
+  double **m;
+  int i, j;
+  m = (double **) malloc (nrows * sizeof(double *));
+  for (i = 0; i < nrows; i++) {
+    m[i] = (double *) malloc (ncols * sizeof(double));
+    if (!m[i]) printf("\nIt's not working");
+  }
+  for (i = 0; i < nrows; i++)
+    for (j = 0; j < ncols; j++)
+      m[i][j] = 0.0;
+  return m;
+}
+
+void freematrix_d(double **I, int rows)
+{
+  int k;
+  for (k = 0; k < rows; k++)
+    free (I[k]);
+}
+
+void matrix_i2d(int **i, double **d, int N)
+{
+  int x, y;
+  for (x = 0; x < N; x++)
+    for (y = 0; y < N; y++)
+      d[y][x] = i[y][x];
+}
+
+void matrix_d2i(double **d, int **i, int N)
+{
+  int x, y;
+  for (x = 0; x < N; x++)
+    for (y = 0; y < N; y++)
+      i[y][x] = d[y][x];
+}
+
+double * dvector(long int N)
+{
+  double *m;
+  m = (double *) malloc (N * sizeof(double));
+  if (!m) printf("\nIt's not working");
+  return m;
+}
+
+void put_matrix_2_vector(double **i, double *f, int N)
+{
+  int l, j, k;
+  k = 0;
+  for (l = 0; l < N; l++)
+    for (j = 0; j < N; j++)
+      f[k++] = i[l][j];
+}
+
+void put_vector_2_matrix(double *f, double **i, int N)
+{
+  int l, j, k;
+  k = 0;
+  for (l = 0; l < N; l++)
+    for (j = 0; j < N; j++)
+      i[l][j] = f[k++];
+}
+
+void set_in_binary() {
+#if defined(EMX)
+  _fsetmode(in, "b");
+#elif defined(MINGW)
+  setmode(STDIN_FILENO, O_BINARY);
+#endif
+}
+
+void set_out_binary() {
+#if defined(EMX)  
+  _fsetmode(out, "b");
+#elif defined(MINGW)
+  setmode(STDOUT_FILENO, O_BINARY);
+#endif
+}
+ 
+void wm_init2() {
+  set_out_binary();
+}
+
+void wm_init1() {
+  set_in_binary();
+}
+
+void wm_init() {
+  set_in_binary();
+  set_out_binary();
+}
+