mdwt,midwt enhancements

HACKING

@@ -79,15 +79,15 @@ addressing so you don't have to worry about row major order and column major
 order. The rwt_printf, mat_offset, offset_row, and offset_col macros will be
 very useful if you need to change any of the code that uses the mat() macro.
 
-To understand the code for the transforms themselves, start with lib/src/dwt.c
+To understand the code for the transforms themselves, start with lib/src/dwt.cc
 which is the best documented of the transforms. The rest of them are written
 and structured in a very similar fashion.
 
 The flow of the code is as follows. One of the transforms is called from MATLAB.
 This invokes one of the wrappers from the mex directory. The function here calls
 rwt_matlab_init in lib/src/init.c which calls other init functions. From here 
 the mex wrapper calls the transform in lib/src. For example, the mdwt function
-for the discrete wavlet transform calls dwt() in the lib/src/dwt.c file. This
+for the discrete wavlet transform calls dwt() in the lib/src/dwt.cc file. This
 function has a few helpers in the same file. It allocates memory necessary for
 the transform in dwt_allocate(), calculates the high and low pass coefficients
 in dwt_coefficients(), performs the convolution in dwt_convolution, and frees

bin/compile.m

@@ -2,20 +2,20 @@
 %
 
 if exist('OCTAVE_VERSION', 'builtin')
-  mkoctfile --mex -v -DOCTAVE_MEX_FILE ../mex/mdwt.c   ../lib/src/dwt.c   ../lib/src/init.c ../lib/src/platform.c -I../lib/inc -o omdwt.mex
-  mkoctfile --mex -v -DOCTAVE_MEX_FILE ../mex/midwt.c  ../lib/src/idwt.c  ../lib/src/init.c ../lib/src/platform.c -I../lib/inc -o omidwt.mex
+  mkoctfile --mex -v -DOCTAVE_MEX_FILE ../mex/mdwt.c   ../lib/src/dwt.cc   ../lib/src/init.c ../lib/src/platform.c -I../lib/inc -o omdwt.mex
+  mkoctfile --mex -v -DOCTAVE_MEX_FILE ../mex/midwt.c  ../lib/src/idwt.cc  ../lib/src/init.c ../lib/src/platform.c -I../lib/inc -o omidwt.mex
   mkoctfile --mex -v -DOCTAVE_MEX_FILE ../mex/mrdwt.c  ../lib/src/rdwt.c  ../lib/src/init.c ../lib/src/platform.c -I../lib/inc -o omrdwt.mex
   mkoctfile --mex -v -DOCTAVE_MEX_FILE ../mex/mirdwt.c ../lib/src/irdwt.c ../lib/src/init.c ../lib/src/platform.c -I../lib/inc -o omirdwt.mex
 else
   x = computer();
   if (x(length(x)-1:length(x)) == '64')
-    mex -v -largeArrayDims ../mex/mdwt.c   ../lib/src/dwt.c   ../lib/src/init.c ../lib/src/platform.c -I../lib/inc -outdir ../bin
-    mex -v -largeArrayDims ../mex/midwt.c  ../lib/src/idwt.c  ../lib/src/init.c ../lib/src/platform.c -I../lib/inc -outdir ../bin
+    mex -v -largeArrayDims ../mex/mdwt.c   ../lib/src/dwt.cc   ../lib/src/init.c ../lib/src/platform.c -I../lib/inc -outdir ../bin
+    mex -v -largeArrayDims ../mex/midwt.c  ../lib/src/idwt.cc  ../lib/src/init.c ../lib/src/platform.c -I../lib/inc -outdir ../bin
     mex -v -largeArrayDims ../mex/mrdwt.c  ../lib/src/rdwt.c  ../lib/src/init.c ../lib/src/platform.c -I../lib/inc -outdir ../bin
     mex -v -largeArrayDims ../mex/mirdwt.c ../lib/src/irdwt.c ../lib/src/init.c ../lib/src/platform.c -I../lib/inc -outdir ../bin
   else
-    mex -v -compatibleArrayDims ../mex/mdwt.c   ../lib/src/dwt.c   ../lib/src/init.c ../lib/src/platform.c -I../lib/inc -outdir ../bin
-    mex -v -compatibleArrayDims ../mex/midwt.c  ../lib/src/idwt.c  ../lib/src/init.c ../lib/src/platform.c -I../lib/inc -outdir ../bin
+    mex -v -compatibleArrayDims ../mex/mdwt.c   ../lib/src/dwt.cc   ../lib/src/init.c ../lib/src/platform.c -I../lib/inc -outdir ../bin
+    mex -v -compatibleArrayDims ../mex/midwt.c  ../lib/src/idwt.cc  ../lib/src/init.c ../lib/src/platform.c -I../lib/inc -outdir ../bin
     mex -v -compatibleArrayDims ../mex/mrdwt.c  ../lib/src/rdwt.c  ../lib/src/init.c ../lib/src/platform.c -I../lib/inc -outdir ../bin
     mex -v -compatibleArrayDims ../mex/mirdwt.c ../lib/src/irdwt.c ../lib/src/init.c ../lib/src/platform.c -I../lib/inc -outdir ../bin
   end

bin/mdwt.m

@@ -1,19 +1,19 @@
 function [y,L] = mdwt(x,h,L)
-%    [y,L] = mdwt(x,h,L);
+%    [y,L] = mdwt(x,h[,L[,transdims]]);
 %
 %    Function computes the discrete wavelet transform y for a 1D or 2D input
 %    signal x using the scaling filter h.
 %
 %    Input:
-%	x : finite length 1D or 2D signal (implicitly periodized)
+%       x : input matrix or nd array ( transform will be done on leading 1 or 2 dims)
 %       h : scaling filter
-%       L : number of levels. In the case of a 1D signal, length(x) must be
-%           divisible by 2^L; in the case of a 2D signal, the row and the
-%           column dimension must be divisible by 2^L. If no argument is
-%           specified, a full DWT is returned for maximal possible L.
+%       L : number of levels. In the case of a 1D transform, size(x,1) must be
+%           divisible by 2^L; for a 2D transform, size(x,2) must also be
+%           divisible by 2^L. The default is the maximal possible L.
+%       transdims: 1 or 2 dimensional transform, default is 2 (if size allows )
 %
 %    Output:
-%       y : the wavelet transform of the signal 
+%       y : the wavelet transform of the signal
 %           (see example to understand the coefficients)
 %       L : number of decomposition levels
 %
@@ -36,7 +36,7 @@
 %
 %    2D Example:
 %
-%       load test_image        
+%       load test_image
 %       h = daubcqf(4,'min');
 %       L = 1;
 %       [y,L] = mdwt(test_image,h,L);
@@ -58,16 +58,21 @@
 %              |   L,H   |  H,H   |
 %              |         |        |
 %              `------------------'
-%       
-%       where 
+%
+%       where
 %            1 : High pass vertically and high pass horizontally
 %            2 : Low pass vertically and high pass horizontally
 %            3 : High pass vertically and low  pass horizontally
-%            4 : Low pass vertically and Low pass horizontally 
+%            4 : Low pass vertically and Low pass horizontally
 %                (scaling coefficients)
 %
+%    4D Tensor with 2D transforms across leading  Example:
 %
-%
+%       x=randn(64,64,3,2);
+%       h=[1,1]*sqrt(.5);
+%       y=mdwt(x,h);
+%       y22 = mdwt(x(:,:,2,2),h);
+%       assert (norm(y(:,:,2,2) - y22 ) < 1e-9)
 %
 %    See also: midwt, mrdwt, mirdwt
 %
@@ -76,7 +81,7 @@
   x = x * 1.0;
   if (exist('L'))
     [y,L] = omdwt(x,h,L);
-  else  
+  else
     [y,L] = omdwt(x,h);
   end
 else

bin/midwt.m

@@ -1,18 +1,17 @@
 function [y,L] = midwt(x,h,L)
-%    [x,L] = midwt(y,h,L);
-% 
-%    Function computes the inverse discrete wavelet transform x for a 1D or
-%    2D input signal y using the scaling filter h.
+%    [x,L] = midwt(y,h[,L[,transdims]])
+%
+%    Compute the inverse discrete wavelet transform x for a 1D or
+%    2D input signal (or tensor of signals) y using the scaling filter h.
 %
 %    Input:
-%	y : finite length 1D or 2D input signal (implicitly periodized)
+%       y : input wavelet domain coefficients (transform is done over leading dims)
 %           (see function mdwt to find the structure of y)
 %       h : scaling filter
-%       L : number of levels. In the case of a 1D signal, length(x) must be
-%           divisible by 2^L; in the case of a 2D signal, the row and the
-%           column dimension must be divisible by 2^L.  If no argument is
-%           specified, a full inverse DWT is returned for maximal possible
-%           L.
+%       L : number of levels. In the case of a 1D transform, size(x,1) must be
+%           divisible by 2^L; for a 2D transform, size(x,2) must also be
+%           divisible by 2^L. The default is the maximal possible L.
+%       transdims: 1 or 2 dimensional transform, default is 2 (if size allows )
 %
 %    Output:
 %       x : periodic reconstructed signal
@@ -36,7 +35,7 @@
 if exist('OCTAVE_VERSION', 'builtin')
   if (exist('L'))
     [y,L] = omidwt(x,h,L);
-  else  
+  else
     [y,L] = omidwt(x,h);
   end
 else

lib/inc/rwt_init.h

@@ -6,19 +6,22 @@
 
 #include "rwt_platform.h"
 
+typedef struct {
+size_t nrows;     /*!< The number of rows in the input matrix. Output matrix will match.  */
+size_t ncols;     /*!< The number of columns in the input matrix. Output matrix will match. */
+size_t nmats;     /*!< The number of matrices (actually just the product of any dimensions past 2) */
+int levels;       /*!< L, the number of levels for the transform. */
+int ncoeff;       /*!< Length of h / the number of scaling coefficients */
+double *scalings; /*!< Wavelet scaling coefficients */
+} rwt_init_params;
+
+typedef enum {NORMAL_DWT, REDUNDANT_DWT, INVERSE_DWT, INVERSE_REDUNDANT_DWT} transform_t;
+
 #if defined(MATLAB_MEX_FILE) || defined(OCTAVE_MEX_FILE)
   #include "mex.h"
   #ifndef OCTAVE_MEX_FILE
     #include "matrix.h"
   #endif
-  typedef struct {
-    size_t nrows;     /*!< The number of rows in the input matrix. Output matrix will match.  */
-    size_t ncols;     /*!< The number of columns in the input matrix. Output matrix will match. */
-    int levels;       /*!< L, the number of levels for the transform. */
-    int ncoeff;       /*!< Length of h / the number of scaling coefficients */
-    double *scalings; /*!< Wavelet scaling coefficients */
-  } rwt_init_params;
-  typedef enum {NORMAL_DWT, REDUNDANT_DWT, INVERSE_DWT, INVERSE_REDUNDANT_DWT} transform_t;
 #endif
 
 #ifdef __cplusplus

lib/inc/rwt_platform.h

@@ -40,11 +40,11 @@
   #define rwt_errormsg(msg) printf("\033[91m%s\033[0m\n", msg);
 #endif
 
-#ifndef max
-  #define max(A,B) (A > B ? A : B)
+#ifndef MAX
+#  define MAX(A,B) (A > B ? A : B)
 #endif
-#ifndef min
-  #define min(A,B) (A < B ? A : B)
+#ifndef MIN
+#  define MIN(A,B) (A < B ? A : B)
 #endif
 #define even(x)  ((x & 1) ? 0 : 1)
 

lib/inc/rwt_transforms.h

@@ -5,6 +5,7 @@
 #define TRANSFORMS_H_
 
 #include <math.h>
+#include "rwt_init.h"
 
 #ifdef __cplusplus
 extern "C" {
@@ -15,8 +16,10 @@ extern "C" {
  *  In all cases it is expected that the output array has already been
  *  allocated prior to calling the transform function.
  */
-void   dwt(double *x, size_t nrows, size_t ncols, double *h, int ncoeff, int levels, double *y);
-void  idwt(double *x, size_t nrows, size_t ncols, double *h, int ncoeff, int levels, double *y);
+void dwt_double(const double *x, double *y,const rwt_init_params * parms);
+void dwt_float(const float *x, float * y,const rwt_init_params * parms);
+void idwt_double(double *x, const double *y,const rwt_init_params * parms);
+void idwt_float(float *x, const float *y,const rwt_init_params * parms);
 void  rdwt(double *x, size_t nrows, size_t ncols, double *h, int ncoeff, int levels, double *yl, double *yh);
 void irdwt(double *x, size_t nrows, size_t ncols, double *h, int ncoeff, int levels, double *yl, double *yh);
 

lib/src/CMakeLists.txt

@@ -1,6 +1,6 @@
 include_directories ("${PROJECT_SOURCE_DIR}/lib/inc")
-add_library(dwt dwt.c)
-add_library(idwt idwt.c)
+add_library(dwt dwt.cc)
+add_library(idwt idwt.cc)
 add_library(irdwt irdwt.c)
 add_library(rdwt rdwt.c)
 add_library(platform platform.c)