wm: Fotopoulos/common.c comparison

comparison Fotopoulos/common.c @ 0:be303a3f5ea8

import

author	Peter Meerwald <pmeerw@cosy.sbg.ac.at>
date	Sun, 12 Aug 2007 13:14:34 +0200
parents
children	cbecc570129d

comparison

equal deleted inserted replaced

--1:000000000000
+:be303a3f5ea8
+#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();
+}

Mercurial > hg > wm

comparison Fotopoulos/common.c @ 0:be303a3f5ea8