comparison util/Rice Wavelet Toolbox/mirdwt_r.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:	06/14/95	16:22:45
+Current Version: MIRDWT.c	2.4
+File Creation Date: Wed Oct 12 08:44:43 1994
+Author: Markus Lang  <[e-mail]>
+
+Copyright (c) 2000 RICE UNIVERSITY. All rights reserved.
+Created by Markus Lang, Department of ECE, Rice University. 
+
+This software is distributed and licensed to you on a non-exclusive 
+basis, free-of-charge. Redistribution and use in source and binary forms, 
+with or without modification, are permitted provided that the following 
+conditions are met:
+
+1. Redistribution of source code must retain the above copyright notice, 
+   this list of conditions and the following disclaimer.
+2. Redistribution in binary form must reproduce the above copyright notice, 
+   this list of conditions and the following disclaimer in the 
+   documentation and/or other materials provided with the distribution.
+3. All advertising materials mentioning features or use of this software 
+   must display the following acknowledgment: This product includes 
+   software developed by Rice University, Houston, Texas and its contributors.
+4. Neither the name of the University nor the names of its contributors 
+   may be used to endorse or promote products derived from this software 
+   without specific prior written permission.
+
+THIS SOFTWARE IS PROVIDED BY WILLIAM MARSH RICE UNIVERSITY, HOUSTON, TEXAS, 
+AND CONTRIBUTORS AS IS AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, 
+BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS 
+FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL RICE UNIVERSITY 
+OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, 
+EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, 
+PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; 
+OR BUSINESS INTERRUPTIONS) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, 
+WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR 
+OTHERWISE), PRODUCT LIABILITY, OR OTHERWISE ARISING IN ANY WAY OUT OF THE 
+USE OF THIS SOFTWARE,  EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+
+For information on commercial licenses, contact Rice University's Office of 
+Technology Transfer at [e-mail] or [phone]
+
+Change History: Fixed the code such that 1D vectors passed to it can be in
+                either passed as a row or column vector. Also took care of 
+		the code such that it will compile with both under standard
+		C compilers as well as for ANSI C compilers
+		Jan Erik Odegard <[e-mail]> Wed Jun 14 1995
+
+                Fix minor bug to allow maximum number of levels
+
+MATLAB description:
+%function x = mirdwt(yl,yh,h,L);
+% 
+% function computes the inverse redundant discrete wavelet transform y for a
+% 1D or  2D input signal. redundant means here that the subsampling after
+% each stage of the forward transform has been omitted. yl contains the
+% lowpass and yl the highpass components as computed, e.g., by mrdwt. In
+% case of a 2D signal the ordering in yh is [lh hl hh lh hl ... ] (first
+% letter refers to row, second to column filtering).  
+%
+%    Input:
+%       yl   : lowpass component
+%       yh   : highpass components
+%       h    : scaling filter
+%       L    : number of levels. in case of a 1D signal length(yl) must be
+%              divisible by 2^L; in case of a 2D signal the row and the
+%              column dimension must be divisible by 2^L.
+%   
+%    Output:
+%	x    : finite length 1D or 2D signal
+%
+% see also: mdwt, midwt, mrdwt
+
+*/
+#include <math.h>
+#include <stdio.h>
+#include <inttypes.h>
+
+#define max(a, b) ((a) > (b) ? (a) : (b))
+#define mat(a, i, j) (*(a + (m*(j)+i)))  /* macro for matrix indices */
+
+#ifdef __STDC__
+MIRDWT(double *x, uintptr_t m, uintptr_t n, double *h, uintptr_t lh, uintptr_t L,
+       double *yl, double *yh)
+#else
+MIRDWT(x, m, n, h, lh, L, yl, yh)
+double *x, *h, *yl, *yh;
+uintptr_t m, n, lh, L;
+#endif
+{
+  double  *g0, *g1, *ydummyll, *ydummylh, *ydummyhl;
+  double *ydummyhh, *xdummyl , *xdummyh, *xh;
+  long i, j;
+  uintptr_t actual_L, actual_m, actual_n, c_o_a, ir, n_c, n_cb, n_c_o, lhm1;
+  uintptr_t ic, n_r, n_rb, n_r_o, c_o_a_p2n, sample_f;
+  xh = (double *)(uintptr_t)mxCalloc(m*n,sizeof(double));
+  xdummyl = (double *)(uintptr_t)mxCalloc(max(m,n),sizeof(double));
+  xdummyh = (double *)(uintptr_t)mxCalloc(max(m,n),sizeof(double));
+  ydummyll = (double *)(uintptr_t)mxCalloc(max(m,n)+lh-1,sizeof(double));
+  ydummylh = (double *)(uintptr_t)mxCalloc(max(m,n)+lh-1,sizeof(double));
+  ydummyhl = (double *)(uintptr_t)mxCalloc(max(m,n)+lh-1,sizeof(double));
+  ydummyhh = (double *)(uintptr_t)mxCalloc(max(m,n)+lh-1,sizeof(double));
+  g0 = (double *)(uintptr_t)mxCalloc(lh,sizeof(double));
+  g1 = (double *)(uintptr_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, uintptr_t lx, double *g0, double *g1, uintptr_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;
+uintptr_t lx, lh;
+#endif
+{
+  uintptr_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;
+  }
+}