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