meschach_utils.c

JPEG2000实现的源码

/**************************************************************************
**
** Copyright (C) 1993 David E. Steward & Zbigniew Leyk, all rights reserved.
**
**			     Meschach Library
** 
** This Meschach Library is provided "as is" without any express 
** or implied warranty of any kind with respect to this software. 
** In particular the authors shall not be liable for any direct, 
** indirect, special, incidental or consequential damages arising 
** in any way from use of the software.
** 
** Everyone is granted permission to copy, modify and redistribute this
** Meschach Library, provided:
**  1.  All copies contain this copyright notice.
**  2.  All modified copies shall carry a notice stating who
**      made the last modification and the date of such modification.
**  3.  No charge is made for this software or works derived from it.  
**      This clause shall not be construed as constraining other software
**      distributed on the same medium as this software, nor is a
**      distribution fee considered a charge.
**
***************************************************************************/

/*****************************************************************************/
/* File: "meschach_utils.c"                                                  */
/* Description: Eigenvector and eigenvalue analysis code.  Modified from     */
/*              original distribution to be stand-alone, compact, and to     */
/*              conform to JPEG2000 VM coding policies.                      */
/* Author: Austin Lan                                                        */
/* Affiliation: Eastman Kodak Company                                        */
/* Version: VM8.5                                                            */
/* Last Revised: 11 September, 2000                                          */
/*****************************************************************************/

#include <local_services.h>
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>
#include <float.h>
#include <ifc.h>

#include "meschach_utils.h"

#define MEM_KEY "MESCHACH"
#define FALSE 0
#define TRUE 1
#define	MACHEPS	DBL_EPSILON /* From float.h; else use 2.22045e-16 */
#define	SQRT2 1.4142135623730949
#define	sgn(x) ( (x) >= 0 ? 1 : -1 )
#ifndef max
  #define max(a,b) ((a) > (b) ? (a) : (b))
#endif
#ifndef min
  #define min(a,b) ((a) > (b) ? (b) : (a))
#endif
#define	MEM_COPY(from,to,size)	memmove((to),(from),(size))
#define	m_copy(in,out) _m_copy(in,out,0,0)
#define	v_copy(in,out) _v_copy(in,out,0)
#define set_col(mat,col,vec) _set_col(mat,col,vec,0)

/* allocate one object of given type */
#define	NEW(type) ((type *)local_malloc(MEM_KEY, sizeof(type)))

/* allocate num objects of given type */
#define	NEW_A(num, type) ((type *)local_malloc(MEM_KEY, (num) * sizeof(type)))

 /* re-allocate arry to have num objects of the given type */
#define	RENEW(var,num,type) \
    ((var)=(type *)((var) ? \
		    local_realloc((void *)(var), ((num) * sizeof(type))) : \
		    local_malloc(MEM_KEY, (num) * sizeof(type))))

/* A[i][j] <- val */
#define	m_set_val(A,i,j,val) ((A)->me[i][j]  = (val) )

/* x[i] <- val */
#define	v_set_val(x,i,val) ((x)->ve[i]  = (val))

/* returns A[i][j] */
#define	m_entry(A,i,j) ((A)->me[i][j])

/* returns x[i] */
#define	v_entry(x,i) ((x)->ve[i])

/*******************************************************************************/
/*******************************************************************************/
/*******************************************************************************/

/* __zero__ -- zeros an array of floating point numbers */
static void __zero__(register double *dp, register int len)
{
   memset((void *)dp, 0, (len * sizeof(double)));
}

/* v_zero -- zero the vector x */
static VEC *v_zero(VEC *x)
{
   if ( x == NULL )
      local_error("v_zero");
   
   __zero__(x->ve,x->dim);
   
   return x;
}

/* v_get -- gets a VEC of dimension 'dim'
   -- Note: initialized to zero */
VEC *v_get(int size)
{
   VEC *vector;
   
   if (size < 0)
      local_error("v_get");
   
   if ((vector=NEW(VEC)) == (VEC *)NULL )
      local_error("v_get");
   memset((void *)vector, 0, sizeof(VEC));
   
   vector->dim = vector->max_dim = size;
   if ((vector->ve=NEW_A(size,double)) == (double *)NULL )
   {
      local_free((void *)vector);
      local_error("v_get");
   }
   memset((void *)vector->ve, 0, (size * sizeof(double)));
   
   return (vector);
}

/* v_resize -- returns the vector x with dim new_dim
   -- x is set to the zero vector */
static VEC *v_resize(VEC *x, int new_dim)
{
   if (new_dim < 0)
      local_error("v_resize");

   if ( ! x )
      return v_get(new_dim);

   /* nothing is changed */
   if (new_dim == x->dim)
      return x;

   if ( x->max_dim == 0 )	/* assume that it's from sub_vec */
      return v_get(new_dim);
   
   if ( new_dim > x->max_dim )
   {
      x->ve = RENEW(x->ve,new_dim,double);
      if ( ! x->ve )
	 local_error("v_resize");
      memset((void *)x->ve, 0, (new_dim * sizeof(double)));
      x->max_dim = new_dim;
   }
   
   if ( new_dim > x->dim )
      __zero__(&(x->ve[x->dim]),new_dim - x->dim);
   x->dim = new_dim;
   
   return x;
}

/* _v_copy -- copies vector into new area */
static VEC *_v_copy(VEC *in, VEC *out, int i0)
{
   if ( in==NULL )
      local_error("_v_copy");
   if ( in==out )
      return (out);
   if ( out==NULL || out->dim < in->dim )
      out = v_resize(out,in->dim);

   MEM_COPY(&(in->ve[i0]),&(out->ve[i0]),(in->dim - i0)*sizeof(double));

   return (out);
}

/* v_free -- returns VEC & asoociated memory back to memory heap */
int v_free(VEC *vec)
{
   if ( vec==(VEC *)NULL || (int)(vec->dim) < 0 )
      /* don't trust it */
      return (-1);
   
   if ( vec->ve == (double *)NULL ) {
      local_free((void *)vec);
   }
   else
   {
      local_free((void *)vec->ve);
      local_free((void *)vec);
   }
   
   return (0);
}

MAT *m_get(int m, int n)
{
   MAT *matrix;
   int i;
   
   if (m < 0 || n < 0)
      local_error("m_get");   
   
   if ((matrix=NEW(MAT)) == (MAT *)NULL )
      local_error("m_get");
   memset((void *)matrix, 0, sizeof(MAT));

   matrix->m = m; matrix->n = matrix->max_n = n;
   matrix->max_m = m; matrix->max_size = m*n;
   if ((matrix->base = NEW_A(m*n,double)) == (double *)NULL )
   {
      local_free((void *)matrix);
      local_error("m_get");
   }
   memset((void *)matrix->base, 0, (m*n*sizeof(double)));
   if ((matrix->me = (double **)local_malloc(MEM_KEY, (m*sizeof(double *)))) == 
       (double **)NULL )
   {	
      local_free((void *)matrix->base);	
      local_free((void *)matrix);
      local_error("m_get");
   }
   
   /* set up pointers */
   for ( i=0; i<m; i++ )
      matrix->me[i] = &(matrix->base[i*n]);
   
   return (matrix);
}

/* m_resize -- returns the matrix A of size new_m x new_n; A is zeroed
   -- if A == NULL on entry then the effect is equivalent to m_get() */
static MAT *m_resize(MAT *A, int new_m, int new_n)
{
   int i;
   int new_max_m, new_max_n, new_size, old_m, old_n;
   
   if (new_m < 0 || new_n < 0)
      local_error("m_resize");

   if ( ! A )
      return m_get(new_m,new_n);

   /* nothing was changed */
   if (new_m == A->m && new_n == A->n)
      return A;

   old_m = A->m;	old_n = A->n;
   if ( new_m > A->max_m )
   {	/* re-allocate A->me */
      A->me = RENEW(A->me,new_m,double *);
      if ( ! A->me )
	 local_error("m_resize");
      memset((void *)A->me, 0, (new_m * sizeof(double *)));
   }
   new_max_m = max(new_m,A->max_m);
   new_max_n = max(new_n,A->max_n);
   
   new_size = new_max_m*new_max_n;
   if ( new_size > A->max_size )
   {	/* re-allocate A->base */
      A->base = RENEW(A->base,new_size,double);
      if ( ! A->base )
	 local_error("m_resize");
      memset((void *)A->base, 0, (new_size * sizeof(double)));
      A->max_size = new_size;
   }
   
   /* now set up A->me[i] */
   for ( i = 0; i < new_m; i++ )
      A->me[i] = &(A->base[i*new_n]);
   
   /* now shift data in matrix */
   if ( old_n > new_n )
   {
      for ( i = 1; i < min(old_m,new_m); i++ )
	 MEM_COPY((char *)&(A->base[i*old_n]),
		  (char *)&(A->base[i*new_n]),
		  sizeof(double)*new_n);
   }
   else if ( old_n < new_n )
   {
      for ( i = (int)(min(old_m,new_m))-1; i > 0; i-- )
      {   /* copy & then zero extra space */
	 MEM_COPY((char *)&(A->base[i*old_n]),
		  (char *)&(A->base[i*new_n]),
		  sizeof(double)*old_n);
	 __zero__(&(A->base[i*new_n+old_n]),(new_n-old_n));
      }
      __zero__(&(A->base[old_n]),(new_n-old_n));
      A->max_n = new_n;
   }
   /* zero out the new rows.. */
   for ( i = old_m; i < new_m; i++ )
      __zero__(&(A->base[i*new_n]),new_n);
   
   A->max_m = new_max_m;
   A->max_n = new_max_n;
   A->max_size = A->max_m*A->max_n;
   A->m = new_m;	A->n = new_n;
   
   return A;
}

/* _m_copy -- copies matrix into new area */
static MAT *_m_copy(MAT *in, MAT *out, int i0, int j0)
{
   int i;

   if ( in==NULL )
      local_error("_m_copy");
   if ( in==out )
      return (out);
   if ( out==NULL || out->m < in->m || out->n < in->n )
      out = m_resize(out,in->m,in->n);

   for ( i=i0; i < in->m; i++ )
      MEM_COPY(&(in->me[i][j0]),&(out->me[i][j0]),
	       (in->n - j0)*sizeof(double));
   
   return (out);
}

/* m_free -- returns MAT & asoociated memory back to memory heap */
int m_free(MAT *mat)
{
   if ( mat==(MAT *)NULL || (int)(mat->m) < 0 ||
	(int)(mat->n) < 0 )
     return (-1);
   
   if ( mat->base != (double *)NULL ) {
      local_free((void *)(mat->base));
   }
   if ( mat->me != (double **)NULL ) {
      local_free((void *)(mat->me));
   }
   local_free((void *)mat);
   
   return (0);
}

/* get_col -- gets a specified column of a matrix and retruns it as a vector */
static VEC *get_col(MAT *mat, int col, VEC *vec)
{
   int i;
   
   if ( mat==(MAT *)NULL )
      local_error("get_col");
   if ( col >= mat->n )
      local_error("get_col");
   if ( vec==(VEC *)NULL || vec->dim<mat->m )
      vec = v_resize(vec,mat->m);
   
   for ( i=0; i<mat->m; i++ )
      vec->ve[i] = mat->me[i][col];
   
   return (vec);
}

/* __ip__ -- inner product */
static double __ip__(register double *dp1, register double *dp2, int len)
{
   register int	i;
   register double sum;

   sum = 0.0;
   for ( i = 0; i < len; i++ )
      sum  += dp1[i]*dp2[i];
    
   return sum;
}

/* _in_prod -- inner product of two vectors from i0 downwards */
static double _in_prod(VEC *a, VEC *b, int i0)
{
   int limit;

   if ( a==(VEC *)NULL || b==(VEC *)NULL )
      local_error("_in_prod");
   limit = min(a->dim,b->dim);
   if ( i0 > limit )
      local_error("_in_prod");

   return __ip__(&(a->ve[i0]),&(b->ve[i0]),(int)(limit-i0));
}

/* __mltadd__ -- scalar multiply and add c.f. v_mltadd() */
static void __mltadd__(register double *dp1, register double *dp2, 
		       register double s, register int len)
{