comparison util/Rice Wavelet Toolbox/mirdwt.c @ 78:f69ae88b8be5

added Rice Wavelet Toolbox with my modification, so it can be compiled on newer systems.

+/*
+File Name: mirdwt.c
+Last Modification Date:	%G%	%U%
+Current Version: %M%	%I%
+File Creation Date: Wed Oct 12 08:44:43 1994
+Author: Markus Lang  <[e-mail]>
+
+Copyright: All software, documentation, and related files in this distribution
+           are Copyright (c) 1994  Rice University
+
+Permission is granted for use and non-profit distribution providing that this
+notice be clearly maintained. The right to distribute any portion for profit
+or as part of any commercial product is specifically reserved for the author.
+
+Change History: Fixed code such that the result has the same dimension as the 
+                input for 1D problems. Also, added some standard error checking.
+		Jan Erik Odegard <[e-mail]> Wed Jun 14 1995
+
+*/
+
+#include <math.h>
+/*#include <malloc.h>*/
+#include <stdio.h>
+#include "mex.h"
+#include "matrix.h"
+#if !defined(_WIN32) && !defined(_WIN64)
+#include <inttypes.h>
+#endif
+#define max(A,B) (A > B ? A : B)
+#define min(A,B) (A < B ? A : B)
+#define even(x)  ((x & 1) ? 0 : 1)
+#define isint(x) ((x - floor(x)) > 0.0 ? 0 : 1)
+
+
+void mexFunction(int nlhs,mxArray *plhs[],int nrhs,const mxArray *prhs[])
+
+{
+  double *x, *h,  *yl, *yh, *Lf, *Lr;
+  intptr_t m, n, mh, nh, h_col, h_row, lh, L, i, po2, j;
+  double mtest, ntest;
+
+  /* check for correct # of input variables */
+  if (nrhs>4){
+    mexErrMsgTxt("There are at most 4 input parameters allowed!");
+    return;
+  }
+  if (nrhs<3){
+    mexErrMsgTxt("There are at least 3 input parameters required!");
+    return;
+  }
+  yl = mxGetPr(prhs[0]);
+  n = mxGetN(prhs[0]); 
+  m = mxGetM(prhs[0]); 
+  yh = mxGetPr(prhs[1]);
+  nh = mxGetN(prhs[1]); 
+  mh = mxGetM(prhs[1]); 
+  h = mxGetPr(prhs[2]);
+  h_col = mxGetN(prhs[2]); 
+  h_row = mxGetM(prhs[2]); 
+  if (h_col>h_row)
+    lh = h_col;
+  else  
+    lh = h_row;
+  if (nrhs == 4){
+    L = (intptr_t) *mxGetPr(prhs[3]);
+    if (L < 0)
+      mexErrMsgTxt("The number of levels, L, must be a non-negative integer");
+  }
+  else /* Estimate L */ {
+    i=n;j=0;
+    while (even(i)){
+      i=(i>>1);
+      j++;
+    }
+    L=m;i=0;
+    while (even(L)){
+      L=(L>>1);
+      i++;
+    }
+    if(min(m,n) == 1)
+      L = max(i,j);
+    else
+      L = min(i,j);
+    if (L==0){
+      mexErrMsgTxt("Maximum number of levels is zero; no decomposition can be performed!");
+      return;
+    }
+  }
+  /* check for consistency of rows and columns of yl, yh */
+  if (min(m,n) > 1){
+    if((m != mh) | (3*n*L != nh)){
+      mexErrMsgTxt("Dimensions of first two input matrices not consistent!");
+      return;
+    }
+  }
+  else{
+    if((m != mh) | (n*L != nh)){
+      mexErrMsgTxt("Dimensions of first two input vectors not consistent!");{
+	return;
+      }
+    }
+  }
+  /* Check the ROW dimension of input */
+  if(m > 1){
+    mtest = (double) m/pow(2.0, (double) L);
+    if (!isint(mtest))
+      mexErrMsgTxt("The matrix row dimension must be of size m*2^(L)");
+  }
+  /* Check the COLUMN dimension of input */
+  if(n > 1){
+    ntest = (double) n/pow(2.0, (double) L);
+    if (!isint(ntest))
+      mexErrMsgTxt("The matrix column dimension must be of size n*2^(L)");
+  }
+  plhs[0] = mxCreateDoubleMatrix(m,n,mxREAL);
+  x = mxGetPr(plhs[0]);
+  plhs[1] = mxCreateDoubleMatrix(1,1,mxREAL);
+  Lr = mxGetPr(plhs[1]);
+  *Lr = L;
+  MIRDWT(x, m, n, h, lh, L, yl, yh);
+}
+#define mat(a, i, j) (*(a + (m*(j)+i)))  /* macro for matrix indices */
+
+#ifdef __STDC__
+MIRDWT(double *x, intptr_t m, intptr_t n, double *h, intptr_t lh, intptr_t L,
+       double *yl, double *yh)
+#else
+MIRDWT(x, m, n, h, lh, L, yl, yh)
+double *x, *h, *yl, *yh;
+intptr_t m, n, lh, L;
+#endif
+{
+  double  *g0, *g1, *ydummyll, *ydummylh, *ydummyhl;
+  double *ydummyhh, *xdummyl , *xdummyh, *xh;
+  intptr_t i, j;
+  intptr_t actual_L, actual_m, actual_n, c_o_a, ir, n_c, n_cb, n_c_o, lhm1;
+  intptr_t ic, n_r, n_rb, n_r_o, c_o_a_p2n, sample_f;
+  xh = (double *)(intptr_t)mxCalloc(m*n,sizeof(double));
+  xdummyl = (double *)(intptr_t)mxCalloc(max(m,n),sizeof(double));
+  xdummyh = (double *)(intptr_t)mxCalloc(max(m,n),sizeof(double));
+  ydummyll = (double *)(intptr_t)mxCalloc(max(m,n)+lh-1,sizeof(double));
+  ydummylh = (double *)(intptr_t)mxCalloc(max(m,n)+lh-1,sizeof(double));
+  ydummyhl = (double *)(intptr_t)mxCalloc(max(m,n)+lh-1,sizeof(double));
+  ydummyhh = (double *)(intptr_t)mxCalloc(max(m,n)+lh-1,sizeof(double));
+  g0 = (double *)(intptr_t)mxCalloc(lh,sizeof(double));
+  g1 = (double *)(intptr_t)mxCalloc(lh,sizeof(double));
+  
+  if (n==1){
+    n = m;
+    m = 1;
+  }
+  /* analysis lowpass and highpass */
+  for (i=0; i<lh; i++){
+    g0[i] = h[i]/2;
+    g1[i] = h[lh-i-1]/2;
+  }
+  for (i=1; i<=lh; i+=2)
+    g1[i] = -g1[i];
+  
+  lhm1 = lh - 1;
+  /* 2^L */
+  sample_f = 1;
+  for (i=1; i<L; i++)
+    sample_f = sample_f*2;
+  actual_m = m/sample_f;
+  actual_n = n/sample_f;
+  /* restore yl in x */
+  for (i=0;i<m*n;i++)
+    x[i] = yl[i];
+  
+  /* main loop */
+  for (actual_L=L; actual_L >= 1; actual_L--){
+    /* actual (level dependent) column offset */
+    if (m==1)
+      c_o_a = n*(actual_L-1);
+    else
+      c_o_a = 3*n*(actual_L-1);
+    c_o_a_p2n = c_o_a + 2*n;
+    
+    /* go by columns in case of a 2D signal*/
+    if (m>1){
+      n_rb = m/actual_m;                 /* # of row blocks per column */
+      for (ic=0; ic<n; ic++){            /* loop over column */
+	for (n_r=0; n_r<n_rb; n_r++){    /* loop within one column */
+	  /* store in dummy variables */
+	  ir = -sample_f + n_r;
+	  for (i=0; i<actual_m; i++){    
+	    ir = ir + sample_f;
+	    ydummyll[i+lhm1] = mat(x, ir, ic);  
+	    ydummylh[i+lhm1] = mat(yh, ir, c_o_a+ic);  
+	    ydummyhl[i+lhm1] = mat(yh, ir,c_o_a+n+ic);  
+	    ydummyhh[i+lhm1] = mat(yh, ir, c_o_a_p2n+ic);   
+	  }
+	  /* perform filtering and adding: first LL/LH, then HL/HH */
+	  bpconv(xdummyl, actual_m, g0, g1, lh, ydummyll, ydummylh); 
+	  bpconv(xdummyh, actual_m, g0, g1, lh, ydummyhl, ydummyhh); 
+	  /* store dummy variables in matrices */
+	  ir = -sample_f + n_r;
+	  for (i=0; i<actual_m; i++){    
+	    ir = ir + sample_f;
+	    mat(x, ir, ic) = xdummyl[i];  
+	    mat(xh, ir, ic) = xdummyh[i];  
+	  }
+	}
+      }
+    }
+    
+    /* go by rows */
+    n_cb = n/actual_n;                 /* # of column blocks per row */
+    for (ir=0; ir<m; ir++){            /* loop over rows */
+      for (n_c=0; n_c<n_cb; n_c++){    /* loop within one row */      
+	/* store in dummy variable */
+	ic = -sample_f + n_c;
+	for  (i=0; i<actual_n; i++){    
+	  ic = ic + sample_f;
+	  ydummyll[i+lhm1] = mat(x, ir, ic);  
+	  if (m>1)
+	    ydummyhh[i+lhm1] = mat(xh, ir, ic);  
+	  else
+	    ydummyhh[i+lhm1] = mat(yh, ir, c_o_a+ic);  
+	} 
+	/* perform filtering lowpass/highpass */
+	bpconv(xdummyl, actual_n, g0, g1, lh, ydummyll, ydummyhh); 
+	/* restore dummy variables in matrices */
+	ic = -sample_f + n_c;
+	for (i=0; i<actual_n; i++){    
+	  ic = ic + sample_f;
+	  mat(x, ir, ic) = xdummyl[i];  
+	}
+      }
+    }
+    sample_f = sample_f/2;
+    actual_m = actual_m*2;
+    actual_n = actual_n*2;
+  }
+}
+
+#ifdef __STDC__
+bpconv(double *x_out, intptr_t lx, double *g0, double *g1, intptr_t lh,
+       double *x_inl, double *x_inh)
+#else
+bpconv(x_out, lx, g0, g1, lh, x_inl, x_inh)
+double *x_inl, *x_inh, *g0, *g1, *x_out;
+intptr_t lx, lh;
+#endif
+{
+  intptr_t i, j;
+  double x0;
+ 
+  for (i=lh-2; i > -1; i--){
+    x_inl[i] = x_inl[lx+i];
+    x_inh[i] = x_inh[lx+i];
+  }
+  for (i=0; i<lx; i++){
+    x0 = 0;
+    for (j=0; j<lh; j++)
+      x0 = x0 + x_inl[j+i]*g0[lh-1-j] +
+	x_inh[j+i]*g1[lh-1-j];
+    x_out[i] = x0;
+  }
+}