Mercurial > hg > smallbox
view util/Rice Wavelet Toolbox/mirdwt.c @ 212:1d134a1b6f95 luisf_dev
Propagate the change from 'mailhe' to 'opt' in two_step_DL to the examples
author | bmailhe |
---|---|
date | Wed, 21 Mar 2012 18:13:02 +0000 |
parents | f69ae88b8be5 |
children |
line wrap: on
line source
/* 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 <lang@jazz.rice.edu> 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 <odegard@ece.rice.edu> 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; } }