Mercurial > hg > smallbox
changeset 78:f69ae88b8be5
added Rice Wavelet Toolbox with my modification, so it can be compiled on newer systems.
line wrap: on
line diff
--- a/SMALLboxSetup.m Wed Mar 23 17:08:55 2011 +0000 +++ b/SMALLboxSetup.m Fri Mar 25 15:27:33 2011 +0000 @@ -103,18 +103,9 @@ 'will not work.\n\n']); end - cd([Sparco_path, FS, 'sparco-1.2', FS, 'tools' ,FS, 'rwt']) - fprintf('Compiling the Rice Wavelet Toolbox MEX interfaces...'); - try - if exist('mdwt' ,'file')~=3, mex mdwt.c mdwt_r.c; end - if exist('midwt' ,'file')~=3, mex midwt.c midwt_r.c; end - if exist('mrdwt' ,'file')~=3, mex mrdwt.c mrdwt_r.c; end - if exist('mirdwt','file')~=3, mex mirdwt.c mirdwt_r.c; end - fprintf('SPARCO Installation Successful!\n'); - catch - warning('Could not compile Rice Wavelet Toolbox MEX interfaces.'); - end + cd(SMALL_path); + fprintf('SPARCO Installation Successful!\n'); else fprintf('\n ******************************************************************'); fprintf('\n\n SPARCO and Rice Wavelet Toolbox are already installed'); @@ -354,28 +345,39 @@ fprintf('\n\n matlab_midi (http://www.kenschutte.com/midi/) is already installed'); end -cd([pwd,FS,'util',FS,'ksvd utils']); +cd([SMALL_path,FS,'util',FS,'ksvd utils']); make cd(SMALL_path); %% - +cd([SMALL_path, FS, 'util', FS, 'Rice Wavelet Toolbox']) + fprintf('Compiling the Rice Wavelet Toolbox MEX interfaces...'); + try + if exist('mdwt' ,'file')~=3, mex mdwt.c; end + if exist('midwt' ,'file')~=3, mex midwt.c; end + if exist('mrdwt' ,'file')~=3, mex mrdwt.c; end + if exist('mirdwt','file')~=3, mex mirdwt.c; end + fprintf('Rice Wavelet Toolbox Installation Successful!\n\n'); + catch + warning('Could not compile Rice Wavelet Toolbox MEX interfaces.\n'); + end +cd(SMALL_path); fprintf('\n ******************************************************************'); fprintf('\n\n Initialising SMALLbox Examples Setup'); - % Need to do a bit of temporary housekeeping first. - cd(SMALL_path); - try - cd(['Problems',FS,'private']); - if exist('addtocols' ,'file')~=3, - fprintf('\n Compiling MEX interfaces for SMALL examples \n'); - make; - end - fprintf('\n SMALLbox Problems Installation Successful! \n'); - catch - fprintf('\n SMALLbox Problems Installation Failed \n'); - end - cd(SMALL_path); +% % Need to do a bit of temporary housekeeping first. +% cd(SMALL_path); +% try +% cd(['Problems',FS,'private']); +% if exist('addtocols' ,'file')~=3, +% fprintf('\n Compiling MEX interfaces for SMALL examples \n'); +% make; +% end +% fprintf('\n SMALLbox Problems Installation Successful! \n'); +% catch +% fprintf('\n SMALLbox Problems Installation Failed \n'); +% end +% cd(SMALL_path);
--- a/examples/Image Denoising/SMALL_ImgDenoise_DL_test_KSVDvsRLSDLA.m Wed Mar 23 17:08:55 2011 +0000 +++ b/examples/Image Denoising/SMALL_ImgDenoise_DL_test_KSVDvsRLSDLA.m Fri Mar 25 15:27:33 2011 +0000 @@ -30,6 +30,7 @@ [pathstr1, name, ext, versn] = fileparts(which('SMALLboxSetup.m')); cd([pathstr1,FS,'data',FS,'images']); load('test_image.mat'); +cd(TMPpath); % [filename,pathname] = uigetfile({'*.png;'},'Select a file containin pre-calculated notes'); % [pathstr, name, ext, versn] = fileparts(filename); % test_image = imread(filename); @@ -41,9 +42,9 @@ % Defining Image Denoising Problem as Dictionary Learning % Problem. As an input we set the number of training patches. for noise_ind=1:1 -for im_num=4:4 +for im_num=2:2 SMALL.Problem = generateImageDenoiseProblem(test_image(im_num).i, 40000, '',256, noise_level(noise_ind)); -SMALL.Problem.name=im_num; +SMALL.Problem.name=int2str(im_num); results(noise_ind,im_num).noisy_psnr=SMALL.Problem.noisy_psnr;
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/util/Rice Wavelet Toolbox/Contents.m Fri Mar 25 15:27:33 2011 +0000 @@ -0,0 +1,18 @@ +% Rice Wavelet Toolbox +% Version 2.4 Dec 2002 +% +% Wavelet Transforms. +% mdwt - Discrete orthogonal wavelet transform using Mallat alg. (1D and 2D) +% midwt - Inverse discrete orthogonal wavelet transform +% mrdwt - Undecimated (redundant) discrete wavelet transform (1D and 2D) +% mirdwt - Inverse undecimated discrete wavelet transform +% daubcqf - Daubecheis filter coefficients +% +% Wavelet Domain Processing. +% denoise - Denoise signal and images by thresholding wavelet coefficients +% HardTh - Hard thresholding +% SoftTh - Soft thresholding +% +% Other. +% makesig - Create Donoho-Johnstone test signals +% compile - compile the Rice Wavelet toolbox
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/util/Rice Wavelet Toolbox/HardTh.m Fri Mar 25 15:27:33 2011 +0000 @@ -0,0 +1,71 @@ +function x = HardTh(y,thld) +% x = HardTh(y,thld); +% +% HARDTH hard thresholds the input signal y with the threshold value +% thld. +% +% Input: +% y : 1D or 2D signal to be thresholded +% thld : threshold value +% +% Output: +% x : Hard thresholded output (x = (abs(y)>thld).*y) +% +% HERE'S AN EASY WAY TO RUN THE EXAMPLES: +% Cut-and-paste the example you want to run to a new file +% called ex.m, for example. Delete out the % at the beginning +% of each line in ex.m (Can use search-and-replace in your editor +% to replace it with a space). Type 'ex' in matlab and hit return. +% +% +% Example: +% y = makesig('WernerSorrows',8); +% thld = 1; +% x = HardTh(y,thld) +% x = 1.5545 5.3175 0 1.6956 -1.2678 0 1.7332 0 +% +% See also: SoftTh +% + +%File Name: HardTh.m +%Last Modification Date: 8/15/95 17:49:37 +%Current Version: HardTh.m 2.4 +%File Creation Date: Mon Jan 31 09:42:50 1994 +%Author: Haitao Guo <harry@jazz.rice.edu> +% +%Copyright (c) 2000 RICE UNIVERSITY. All rights reserved. +%Created by Haitao Guo, 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 techtran@rice.edu or (713) 348-6173 + +x = (abs(y) > thld).*y;
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/util/Rice Wavelet Toolbox/INSTALL Fri Mar 25 15:27:33 2011 +0000 @@ -0,0 +1,26 @@ +In order to install this distribution of Rice Wavelet Tools version 2.4 +released - <Dec 1 2002> + +1. Properly set up your system to create MEX-files. Please refer to the + "Matlab Application Program Guide" to properly set up of your matlab + and C-compiler to be able to compile C-mex files on your system. + All reference documentations are available on the MathWorks web page: + www.mathworks.com + +2. Make a toolbox directory and uncompress/extract all the files. + For example, in the unix environment, + + gunzip rwt.tar.gz + tar xvf rwt.tar + +3. Run MATLAB and change to the temporary directory containing the files. + +4. Compile the toolbox by executing the Matlab command + + compile + +5. Add the toolbox directory to your Matlab path. + +6. For further instructions, please refer to the README file. + +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/util/Rice Wavelet Toolbox/LICENSE Fri Mar 25 15:27:33 2011 +0000 @@ -0,0 +1,38 @@ +This "rice-wlet-tools", version 2.4 +Released - <Dec 1 2002> + +CONDITIONS FOR USE: +Copyright (c) 2000 RICE UNIVERSITY. All rights reserved. + +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 techtran@rice.edu or (713) 348-6173 +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/util/Rice Wavelet Toolbox/README Fri Mar 25 15:27:33 2011 +0000 @@ -0,0 +1,38 @@ +Rice Wavelet Tools version 2.4 +Released - <Dec 1 2002> + +INSTALLATION: +To install this distribution of Rice Wavelet Tools see the INSTALL file. + +SOURCE: + www.dsp.rice.edu/software/RWT + +EMAIL: +For bug reports and questions, send email to wlet-tools@rice.edu + +CONDITIONS FOR USE: +See the LICENSE file + +TOOLBOX FUNCTIONS: + + Wavelet Transforms + mdwt - Discrete orthogonal wavelet transform using the Mallat algorithm (1D and 2D) + midwt - Inverse discrete orthogonal wavelet transform + mrdwt - Undecimated (redundant) discrete wavelet transform (1D and 2D) + mirdwt - Inverse undecimated discrete wavelet transform + daubcqf - Daubechies filter coefficients + + Wavelet Domain Processing + denoise - Denoise signals and images by thresholding wavelet coefficients + HardTh - Hard thresholding + SoftTh - Soft thresholding + + Other + makesig - Create Donoho-Johnstone test signals + compile - Compile the Rice Wavelet Toolbox + +Functions omitted in this version of toolbox can be found in +version 2.01 at www.dsp.rice.edu/software/RWT/RWT-2.01.tar.Z + + +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/util/Rice Wavelet Toolbox/README.sparco Fri Mar 25 15:27:33 2011 +0000 @@ -0,0 +1,21 @@ +What are these funny ".mhelp" files?! +------------------------------------- + +Some users of Matlab 7.0.4 have reported that the Rice Wavelet Toolbox +MEX interfaces aren't being found. Instead, the .m files (which only +contain the help text) are taking precedence, and function calls to +these routines cause an error. + +A simple workaround is to move these .m help files so that Matlab +doesn't have to choose between them and the MEX versions. We've +renamed these files: + +mdwt.m to mdwt.mhelp +midwt.m to midwt.mhelpm +mirdwt.m to mirdwt.mhelpm +mrdwt.dll to mrdwt.mhelpm + +-- The Sparco development team + (Michael P. Friedlander and Ewout van den Berg) +$Id: README.sparco 830 2008-03-25 22:28:04Z mpf $ +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/util/Rice Wavelet Toolbox/SoftTh.m Fri Mar 25 15:27:33 2011 +0000 @@ -0,0 +1,79 @@ +function x = SoftTh(y,thld) +% x = SoftTh(y,thld); +% +% SOFTTH soft thresholds the input signal y with the threshold value +% thld. +% +% Input: +% y : 1D or 2D signal to be thresholded +% thld : Threshold value +% +% Output: +% x : Soft thresholded output (x = sign(y)(|y|-thld)_+) +% +% HERE'S AN EASY WAY TO RUN THE EXAMPLES: +% Cut-and-paste the example you want to run to a new file +% called ex.m, for example. Delete out the % at the beginning +% of each line in ex.m (Can use search-and-replace in your editor +% to replace it with a space). Type 'ex' in matlab and hit return. +% +% +% Example: +% y = makesig('Doppler',8); +% thld = 0.2; +% x = SoftTh(y,thld) +% x = 0 0 0 -0.0703 0 0.2001 0.0483 0 +% +% See also: HardTh +% +% Reference: +% "De-noising via Soft-Thresholding" Tech. Rept. Statistics, +% Stanford, 1992. D.L. Donoho. +% + +%File Name: SoftTh.m +%Last Modification Date: 8/15/95 17:49:48 +%Current Version: SoftTh.m 2.4 +%File Creation Date: Mon Mar 7 10:38:45 1994 +%Author: Haitao Guo <harry@jazz.rice.edu> +% +%Copyright (c) 2000 RICE UNIVERSITY. All rights reserved. +%Created by Haitao Guo, 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 techtran@rice.edu or (713) 348-6173 +% +%Change History: +% + +x = abs(y); +x = sign(y).*(x >= thld).*(x - thld);
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/util/Rice Wavelet Toolbox/compile.m Fri Mar 25 15:27:33 2011 +0000 @@ -0,0 +1,42 @@ +% COMPILE compiles the c files and generates mex files. +% + +%Copyright (c) 2000 RICE UNIVERSITY. All rights reserved. +%Created by Haitao Guo, 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 techtran@rice.edu or (713) 348-6173 + +mex mdwt.c +mex midwt.c +mex mrdwt.c +mex mirdwt.c
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/util/Rice Wavelet Toolbox/daubcqf.m Fri Mar 25 15:27:33 2011 +0000 @@ -0,0 +1,106 @@ +function [h_0,h_1] = daubcqf(N,TYPE) +% [h_0,h_1] = daubcqf(N,TYPE); +% +% Function computes the Daubechies' scaling and wavelet filters +% (normalized to sqrt(2)). +% +% Input: +% N : Length of filter (must be even) +% TYPE : Optional parameter that distinguishes the minimum phase, +% maximum phase and mid-phase solutions ('min', 'max', or +% 'mid'). If no argument is specified, the minimum phase +% solution is used. +% +% Output: +% h_0 : Minimal phase Daubechies' scaling filter +% h_1 : Minimal phase Daubechies' wavelet filter +% +% Example: +% N = 4; +% TYPE = 'min'; +% [h_0,h_1] = daubcqf(N,TYPE) +% h_0 = 0.4830 0.8365 0.2241 -0.1294 +% h_1 = 0.1294 0.2241 -0.8365 0.4830 +% +% Reference: "Orthonormal Bases of Compactly Supported Wavelets", +% CPAM, Oct.89 +% + +%File Name: daubcqf.m +%Last Modification Date: 01/02/96 15:12:57 +%Current Version: daubcqf.m 2.4 +%File Creation Date: 10/10/88 +%Author: Ramesh Gopinath <ramesh@dsp.rice.edu> +% +%Copyright (c) 2000 RICE UNIVERSITY. All rights reserved. +%Created by Ramesh Gopinath, 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 techtran@rice.edu or (713) 348-6173 + +if(nargin < 2), + TYPE = 'min'; +end; +if(rem(N,2) ~= 0), + error('No Daubechies filter exists for ODD length'); +end; +K = N/2; +a = 1; +p = 1; +q = 1; +h_0 = [1 1]; +for j = 1:K-1, + a = -a * 0.25 * (j + K - 1)/j; + h_0 = [0 h_0] + [h_0 0]; + p = [0 -p] + [p 0]; + p = [0 -p] + [p 0]; + q = [0 q 0] + a*p; +end; +q = sort(roots(q)); +qt = q(1:K-1); +if TYPE=='mid', + if rem(K,2)==1, + qt = q([1:4:N-2 2:4:N-2]); + else + qt = q([1 4:4:K-1 5:4:K-1 N-3:-4:K N-4:-4:K]); + end; +end; +h_0 = conv(h_0,real(poly(qt))); +h_0 = sqrt(2)*h_0/sum(h_0); %Normalize to sqrt(2); +if(TYPE=='max'), + h_0 = fliplr(h_0); +end; +if(abs(sum(h_0 .^ 2))-1 > 1e-4) + error('Numerically unstable for this value of "N".'); +end; +h_1 = rot90(h_0,2); +h_1(1:2:N)=-h_1(1:2:N);
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/util/Rice Wavelet Toolbox/denoise.m Fri Mar 25 15:27:33 2011 +0000 @@ -0,0 +1,245 @@ +function [xd,xn,option] = denoise(x,h,type,option) +% [xd,xn,option] = denoise(x,h,type,option); +% +% DENOISE is a generic program for wavelet based denoising. +% The program will denoise the signal x using the 2-band wavelet +% system described by the filter h using either the traditional +% discrete wavelet transform (DWT) or the linear shift invariant +% discrete wavelet transform (also known as the undecimated DWT +% (UDWT)). +% +% Input: +% x : 1D or 2D signal to be denoised +% h : Scaling filter to be applied +% type : Type of transform (Default: type = 0) +% 0 --> Discrete wavelet transform (DWT) +% 1 --> Undecimated DWT (UDWT) +% option : Default settings is marked with '*': +% *type = 0 --> option = [0 3.0 0 0 0 0] +% type = 1 --> option = [0 3.6 0 1 0 0] +% option(1) : Whether to threshold low-pass part +% 0 --> Don't threshold low pass component +% 1 --> Threshold low pass component +% option(2) : Threshold multiplier, c. The threshold is +% computed as: +% thld = c*MAD(noise_estimate)). +% The default values are: +% c = 3.0 for the DWT based denoising +% c = 3.6 for the UDWT based denoising +% option(3) : Type of variance estimator +% 0 --> MAD (mean absolute deviation) +% 1 --> STD (classical numerical std estimate) +% option(4) : Type of thresholding +% 0 --> Soft thresholding +% 1 --> Hard thresholding +% option(5) : Number of levels, L, in wavelet decomposition. By +% setting this to the default value '0' a maximal +% decomposition is used. +% option(6) : Actual threshold to use (setting this to +% anything but 0 will mean that option(3) +% is ignored) +% +% Output: +% xd : Estimate of noise free signal +% xn : The estimated noise signal (x-xd) +% option : A vector of actual parameters used by the +% program. The vector is configured the same way as +% the input option vector with one added element +% option(7) = type. +% +% HERE'S AN EASY WAY TO RUN THE EXAMPLES: +% Cut-and-paste the example you want to run to a new file +% called ex.m, for example. Delete out the % at the beginning +% of each line in ex.m (Can use search-and-replace in your editor +% to replace it with a space). Type 'ex' in matlab and hit return. +% +% Example 1: +% h = daubcqf(6); [s,N] = makesig('Doppler'); n = randn(1,N); +% x = s + n/10; % (approximately 10dB SNR) +% figure;plot(x);hold on;plot(s,'r'); +% +% %Denoise x with the default method based on the DWT +% [xd,xn,opt1] = denoise(x,h); +% figure;plot(xd);hold on;plot(s,'r'); +% +% %Denoise x using the undecimated (LSI) wavelet transform +% [yd,yn,opt2] = denoise(x,h,1); +% figure;plot(yd);hold on;plot(s,'r'); +% +% Example 2: (on an image) +% h = daubcqf(6); load lena; +% noisyLena = lena + 25 * randn(size(lena)); +% figure; colormap(gray); imagesc(lena); title('Original Image'); +% figure; colormap(gray); imagesc(noisyLena); title('Noisy Image'); +% Denoise lena with the default method based on the DWT +% [denoisedLena,xn,opt1] = denoise(noisyLena,h); +% figure; colormap(gray); imagesc(denoisedLena); title('denoised Image'); +% +% +% See also: mdwt, midwt, mrdwt, mirdwt, SoftTh, HardTh, setopt +% + +%File Name: denoise.m +%Last Modification Date: 04/15/97 10:44:28 +%Current Version: denoise.m 2.4 +%File Creation Date: Mon Feb 20 08:33:15 1995 +%Author: Jan Erik Odegard <odegard@ece.rice.edu> +% +%Copyright (c) 2000 RICE UNIVERSITY. All rights reserved. +%Created by Jan Erik Odegard, 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 techtran@rice.edu or (713) 348-6173 +% +%Change History: Fixed output of function and an error in the computation +% of the threshold for redundant denoising. +% <Jan Erik Odegard> <Mon Jul 31, 1995> +% +% This code is composed of several of our old codes for +% wavelet based denoising. In an effort to make the mess +% more manageable we decided to create on code that would +% handle all the various wavelet based denoising methods. +% However, only time will show (as we discover new and +% improved forms of denoising) if we can succeed in our goals. +% <Jan Erik Odegard> <Thu May 11, 1995> +% + +if(nargin < 2) + error('You need to provide at least 2 inputs: x and h'); +end; +if(nargin < 3), + type = 0; + option = []; +elseif(nargin < 4) + option = []; +end; +if(isempty(type)), + type = 0; +end; +if(type == 0), + default_opt = [0 3.0 0 0 0 0]; +elseif(type == 1), + default_opt = [0 3.6 0 1 0 0]; +else, + error(['Unknown denoising method',10,... + 'If it is any good we need to have a serious talk :-)']); +end; +option = setopt(option,default_opt); +[mx,nx] = size(x); +dim = min(mx,nx); +if(dim == 1), + n = max(mx,nx); +else, + n = dim; +end; +if(option(5) == 0), + L = floor(log2(n)); +else + L = option(5); +end; +if(type == 0), % Denoising by DWT + xd = mdwt(x,h,L); + if (option(6) == 0), + tmp = xd(floor(mx/2)+1:mx,floor(nx/2)+1:nx); + if(option(3) == 0), + thld = option(2)*median(abs(tmp(:)))/.67; + elseif(option(3) == 1), + thld = option(2)*std(tmp(:)); + else + error('Unknown threshold estimator, Use either MAD or STD'); + end; + else, + thld = option(6); + end; + if(dim == 1) + ix = 1:n/(2^L); + ykeep = xd(ix); + else + ix = 1:mx/(2^L); + jx = 1:nx/(2^L); + ykeep = xd(ix,jx); + end; + if(option(4) == 0), + xd = SoftTh(xd,thld); + elseif(option(4) == 1), + xd = HardTh(xd,thld); + else, + error('Unknown threshold rule. Use either Soft (0) or Hard (1)'); + end; + if (option(1) == 0), + if(dim == 1), + xd(ix) = ykeep; + else, + xd(ix,jx) = ykeep; + end; + end; + xd = midwt(xd,h,L); +elseif(type == 1), % Denoising by UDWT + [xl,xh] = mrdwt(x,h,L); + if(dim == 1), + c_offset = 1; + else, + c_offset = 2*nx + 1; + end; + if (option(6) == 0), + tmp = xh(:,c_offset:c_offset+nx-1); + if(option(3) == 0), + thld = option(2)*median(abs(tmp(:)))/.67; + elseif(option(3) == 1), + thld = option(2)*std(tmp(:)); + else + error('Unknown threshold estimator, Use either MAD or STD'); + end; + else, + thld = option(6); + end; + if(option(4) == 0), + xh = SoftTh(xh,thld); + if(option(1) == 1), + xl = SoftTh(xl,thld); + end; + elseif(option(4) == 1), + xh = HardTh(xh,thld); + if(option(1) == 1), + xl = HardTh(xl,thld); + end; + else, + error('Unknown threshold rule. Use either Soft (0) or Hard (1)'); + end; + xd = mirdwt(xl,xh,h,L); +else, % Denoising by unknown method + error(['Unknown denoising method',10,... + 'If it is any good we need to have a serious talk :-)']); +end; +option(6) = thld; +option(7) = type; +xn = x - xd;
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/util/Rice Wavelet Toolbox/makesig.m Fri Mar 25 15:27:33 2011 +0000 @@ -0,0 +1,211 @@ +function [x,N] = makesig(SigName,N) +% [x,N] = makesig(SigName,N) Creates artificial test signal identical to the +% standard test signals proposed and used by D. Donoho and I. Johnstone +% in WaveLab (- a matlab toolbox developed by Donoho et al. the statistics +% department at Stanford University). +% +% Input: SigName - Name of the desired signal (Default 'all') +% 'AllSig' (Returns a matrix with all the signals) +% 'HeaviSine' +% 'Bumps' +% 'Blocks' +% 'Doppler' +% 'Ramp' +% 'Cusp' +% 'Sing' +% 'HiSine' +% 'LoSine' +% 'LinChirp' +% 'TwoChirp' +% 'QuadChirp' +% 'MishMash' +% 'Werner Sorrows' (Heisenberg) +% 'Leopold' (Kronecker) +% N - Length in samples of the desired signal (Default 512) +% +% Output: x - vector/matrix of test signals +% N - length of signal returned +% +% See also: +% +% References: +% WaveLab can be accessed at +% www_url: http://playfair.stanford.edu/~wavelab/ +% Also see various articles by D.L. Donoho et al. at +% web_url: http://playfair.stanford.edu/ + +%File Name: makesig.m +%Last Modification Date: 08/30/95 15:52:03 +%Current Version: makesig.m 2.4 +%File Creation Date: Thu Jun 8 10:31:11 1995 +%Author: Jan Erik Odegard <odegard@ece.rice.edu> +% +%Copyright (c) 2000 RICE UNIVERSITY. All rights reserved. +%Created by Jan Erik Odegard, 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 techtran@rice.edu or (713) 348-6173 +% +%Change History: This m-file is a copy of the code provided with WaveLab +% customized to be consistent with RWT. +% Jan Erik Odegard <odegard@ece.rice.edu> Thu Jun 8 1995 +% + +if(nargin < 1) + SigName = 'AllSig'; + N = 512; +elseif(nargin == 1) + N = 512; +end; +t = (1:N) ./N; +x = []; +y = []; +if(strcmp(SigName,'HeaviSine') | strcmp(SigName,'AllSig')), + y = 4.*sin(4*pi.*t); + y = y - sign(t - .3) - sign(.72 - t); +end; +x = [x;y]; +y = []; +if(strcmp(SigName,'Bumps') | strcmp(SigName,'AllSig')), + pos = [ .1 .13 .15 .23 .25 .40 .44 .65 .76 .78 .81]; + hgt = [ 4 5 3 4 5 4.2 2.1 4.3 3.1 5.1 4.2]; + wth = [.005 .005 .006 .01 .01 .03 .01 .01 .005 .008 .005]; + y = zeros(size(t)); + for j =1:length(pos) + y = y + hgt(j)./( 1 + abs((t - pos(j))./wth(j))).^4; + end +end; +x = [x;y]; +y = []; +if(strcmp(SigName,'Blocks') | strcmp(SigName,'AllSig')), + pos = [ .1 .13 .15 .23 .25 .40 .44 .65 .76 .78 .81]; + hgt = [4 (-5) 3 (-4) 5 (-4.2) 2.1 4.3 (-3.1) 2.1 (-4.2)]; + y = zeros(size(t)); + for j=1:length(pos) + y = y + (1 + sign(t-pos(j))).*(hgt(j)/2) ; + end +end; +x = [x;y]; +y = []; +if(strcmp(SigName,'Doppler') | strcmp(SigName,'AllSig')), + y = sqrt(t.*(1-t)).*sin((2*pi*1.05) ./(t+.05)); +end; +x = [x;y]; +y = []; +if(strcmp(SigName,'Ramp') | strcmp(SigName,'AllSig')), + y = t - (t >= .37); +end; +x = [x;y]; +y = []; +if(strcmp(SigName,'Cusp') | strcmp(SigName,'AllSig')), + y = sqrt(abs(t - .37)); +end; +x = [x;y]; +y = []; +if(strcmp(SigName,'Sing') | strcmp(SigName,'AllSig')), + k = floor(N * .37); + y = 1 ./abs(t - (k+.5)/N); +end; +x = [x;y]; +y = []; +if(strcmp(SigName,'HiSine') | strcmp(SigName,'AllSig')), + y = sin( pi * (N * .6902) .* t); +end; +x = [x;y]; +y = []; +if(strcmp(SigName,'LoSine') | strcmp(SigName,'AllSig')), + y = sin( pi * (N * .3333) .* t); +end; +x = [x;y]; +y = []; +if(strcmp(SigName,'LinChirp') | strcmp(SigName,'AllSig')), + y = sin(pi .* t .* ((N .* .125) .* t)); +end; +x = [x;y]; +y = []; +if(strcmp(SigName,'TwoChirp') | strcmp(SigName,'AllSig')), + y = sin(pi .* t .* (N .* t)) + sin((pi/3) .* t .* (N .* t)); +end; +x = [x;y]; +y = []; +if(strcmp(SigName,'QuadChirp') | strcmp(SigName,'AllSig')), + y = sin( (pi/3) .* t .* (N .* t.^2)); +end; +x = [x;y]; +y = []; +if(strcmp(SigName,'MishMash') | strcmp(SigName,'AllSig')), + % QuadChirp + LinChirp + HiSine + y = sin( (pi/3) .* t .* (N .* t.^2)) ; + y = y + sin( pi * (N * .6902) .* t); + y = y + sin(pi .* t .* (N .* .125 .* t)); +end; +x = [x;y]; +y = []; +if(strcmp(SigName,'WernerSorrows') | strcmp(SigName,'AllSig')), + y = sin( pi .* t .* (N/2 .* t.^2)) ; + y = y + sin( pi * (N * .6902) .* t); + y = y + sin(pi .* t .* (N .* t)); + pos = [ .1 .13 .15 .23 .25 .40 .44 .65 .76 .78 .81]; + hgt = [ 4 5 3 4 5 4.2 2.1 4.3 3.1 5.1 4.2]; + wth = [.005 .005 .006 .01 .01 .03 .01 .01 .005 .008 .005]; + for j =1:length(pos) + y = y + hgt(j)./( 1 + abs((t - pos(j))./wth(j))).^4; + end +end; +x = [x;y]; +y = []; +if(strcmp(SigName,'Leopold') | strcmp(SigName,'AllSig')), + y = (t == floor(.37 * N)/N); % Kronecker +end; +x = [x;y]; +y = []; + +% disp(sprintf('MakeSignal: I don*t recognize << %s>>',SigName)) +% disp('Allowable SigNames are:') +% disp('AllSig'), +% disp('HeaviSine'), +% disp('Bumps'), +% disp('Blocks'), +% disp('Doppler'), +% disp('Ramp'), +% disp('Cusp'), +% disp('Crease'), +% disp('Sing'), +% disp('HiSine'), +% disp('LoSine'), +% disp('LinChirp'), +% disp('TwoChirp'), +% disp('QuadChirp'), +% disp('MishMash'), +% disp('WernerSorrows'), +% disp('Leopold'), +%end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/util/Rice Wavelet Toolbox/mdwt.c Fri Mar 25 15:27:33 2011 +0000 @@ -0,0 +1,215 @@ +/* +File Name: mdwt.c +Last Modification Date: 06/14/95 12:56:43 +Current Version: mdwt.c 1.5 +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 + +MATLAB gateway for MDWT.c, discrete wavelet transform +*/ +#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) +#define mat(a, i, j) (*(a + (m*(j)+i))) /* macro for matrix indices */ + +void mexFunction(const int nlhs,mxArray *plhs[],const int nrhs,const mxArray *prhs[]) + +{ + double *x, *h, *y, *Lf, *Lr; + intptr_t m, n, h_col, h_row, lh, L, i, po2, j; + double mtest, ntest; + + /* check for correct # of input variables */ + if (nrhs>3){ + mexErrMsgTxt("There are at most 3 input parameters allowed!"); + return; + } + if (nrhs<2){ + mexErrMsgTxt("There are at least 2 input parameters required!"); + return; + } + x = mxGetPr(prhs[0]); + n = mxGetN(prhs[0]); + m = mxGetM(prhs[0]); + h = mxGetPr(prhs[1]); + h_col = mxGetN(prhs[1]); + h_row = mxGetM(prhs[1]); + if (h_col>h_row) + lh = h_col; + else + lh = h_row; + if (nrhs == 3){ + L = (intptr_t) *mxGetPr(prhs[2]); + 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 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); + y = mxGetPr(plhs[0]); + plhs[1] = mxCreateDoubleMatrix(1,1,mxREAL); + Lr = mxGetPr(plhs[1]); + *Lr = L; + MDWT(x, m, n, h, lh, L, y); +} +#ifdef __STDC__ +MDWT(double *x, intptr_t m, intptr_t n, double *h, intptr_t lh, intptr_t L, double *y) +#else +MDWT(x, m, n, h, lh, L, y) +double *x, *h, *y; +intptr_t m, n, lh, L; +#endif +{ + double *h0, *h1, *ydummyl, *ydummyh, *xdummy; + int *prob; + intptr_t i, j; + intptr_t actual_L, actual_m, actual_n, r_o_a, c_o_a, ir, ic, lhm1; + + xdummy = (double *)mxCalloc(max(m,n)+lh-1,sizeof(double)); + ydummyl =(double *) (intptr_t)mxCalloc(max(m,n),sizeof(double)); + ydummyh = (double *)(intptr_t)mxCalloc(max(m,n),sizeof(double)); + h0 =(double *)(intptr_t)mxCalloc(lh,sizeof(double)); + h1 = (double *)(intptr_t)mxCalloc(lh,sizeof(double)); + + + /* analysis lowpass and highpass */ + if (n==1){ + n = m; + m = 1; + } + for (i=0; i<lh; i++){ + h0[i] = h[lh-i-1]; + h1[i] =h[i]; + } + for (i=0; i<lh; i+=2) + h1[i] = -h1[i]; + + lhm1 = lh - 1; + actual_m = 2*m; + actual_n = 2*n; + + /* main loop */ + for (actual_L=1; actual_L <= L; actual_L++){ + if (m==1) + actual_m = 1; + else{ + actual_m = actual_m/2; + r_o_a = actual_m/2; + } + actual_n = actual_n/2; + c_o_a = actual_n/2; + + /* go by rows */ + for (ir=0; ir<actual_m; ir++){ /* loop over rows */ + /* store in dummy variable */ + for (i=0; i<actual_n; i++) + if (actual_L==1) + xdummy[i] = mat(x, ir, i); + else + xdummy[i] = mat(y, ir, i); + /* perform filtering lowpass and highpass*/ + fpsconv(xdummy, actual_n, h0, h1, lhm1, ydummyl, ydummyh); + /* restore dummy variables in matrices */ + ic = c_o_a; + for (i=0; i<c_o_a; i++){ + mat(y, ir, i) = ydummyl[i]; + mat(y, ir, ic++) = ydummyh[i]; + } + } + + /* go by columns in case of a 2D signal*/ + if (m>1){ + for (ic=0; ic<actual_n; ic++){ /* loop over column */ + /* store in dummy variables */ + for (i=0; i<actual_m; i++) + xdummy[i] = mat(y, i, ic); + /* perform filtering lowpass and highpass*/ + fpsconv(xdummy, actual_m, h0, h1, lhm1, ydummyl, ydummyh); + /* restore dummy variables in matrix */ + ir = r_o_a; + for (i=0; i<r_o_a; i++){ + mat(y, i, ic) = ydummyl[i]; + mat(y, ir++, ic) = ydummyh[i]; + } + } + } + } +} + +#ifdef __STDC__ +fpsconv(double *x_in, intptr_t lx, double *h0, double *h1, intptr_t lhm1, + double *x_outl, double *x_outh) +#else +fpsconv(x_in, lx, h0, h1, lhm1, x_outl, x_outh) +double *x_in, *h0, *h1, *x_outl, *x_outh; +intptr_t lx, lhm1; +#endif + +{ + intptr_t i, j, ind; + double x0, x1; + + for (i=lx; i < lx+lhm1; i++) + x_in[i] = *(x_in+(i-lx)); + ind = 0; + for (i=0; i<(lx); i+=2){ + x0 = 0; + x1 = 0; + for (j=0; j<=lhm1; j++){ + x0 = x0 + x_in[i+j]*h0[lhm1-j]; + x1 = x1 + x_in[i+j]*h1[lhm1-j]; + } + x_outl[ind] = x0; + x_outh[ind++] = x1; + } +}
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/util/Rice Wavelet Toolbox/mdwt.mhelp Fri Mar 25 15:27:33 2011 +0000 @@ -0,0 +1,132 @@ +function [y,L] = mdwt(x,h,L); +% [y,L] = mdwt(x,h,L); +% +% 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) +% 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. +% +% Output: +% y : the wavelet transform of the signal +% (see example to understand the coefficients) +% L : number of decomposition levels +% +% 1D Example: +% x = makesig('LinChirp',8); +% h = daubcqf(4,'min'); +% L = 2; +% [y,L] = mdwt(x,h,L) +% +% 1D Example's output and explanation: +% +% y = [1.1097 0.8767 0.8204 -0.5201 -0.0339 0.1001 0.2201 -0.1401] +% L = 2 +% +% The coefficients in output y are arranged as follows +% +% y(1) and y(2) : Scaling coefficients (lowest frequency) +% y(3) and y(4) : Band pass wavelet coefficients +% y(5) to y(8) : Finest scale wavelet coefficients (highest frequency) +% +% 2D Example: +% +% load test_image +% h = daubcqf(4,'min'); +% L = 1; +% [y,L] = mdwt(test_image,h,L); +% +% 2D Example's output and explanation: +% +% The coefficients in y are arranged as follows. +% +% .------------------. +% | | | +% | 4 | 2 | +% | | | +% | L,L | H,L | +% | | | +% -------------------- +% | | | +% | 3 | 1 | +% | | | +% | L,H | H,H | +% | | | +% `------------------' +% +% 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 +% (scaling coefficients) +% +% +% +% +% See also: midwt, mrdwt, mirdwt +% + +%File Name: mdwt.m +%Last Modification Date: 08/07/95 15:13:25 +%Current Version: mdwt.m 2.4 +%File Creation Date: Wed Oct 19 10:51:58 1994 +%Author: Markus Lang <lang@jazz.rice.edu> +% +%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 techtran@rice.edu or (713) 348-6173 +% +%Change History: +% +%Modification #1 +%Mon Aug 7 11:42:11 CDT 1995 +%Rebecca Hindman <hindman@ece.rice.edu> +%Added L to function line so that it can be displayed as an output +% +%Change History: +% +%Modification #1 +%Thu Mar 2 13:07:11 CDT 2000 +%Ramesh Neelamani<neelsh@ece.rice.edu> +%Revamped the help file +% + + +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/util/Rice Wavelet Toolbox/mdwt_r.c Fri Mar 25 15:27:33 2011 +0000 @@ -0,0 +1,181 @@ +/* +File Name: MDWT.c +Last Modification Date: 06/14/95 13:15:44 +Current Version: MDWT.c 2.4 +File Creation Date: Wed Oct 19 10:51:58 1994 +Author: Markus Lang <lang@jazz.rice.edu> + +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 techtran@rice.edu or (713) 348-6173 + +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 <odegard@ece.rice.edu> Wed Jun 14 1995 + +%y = mdwt(x,h,L); +% +% function computes the discrete wavelet transform y for a 1D or 2D input +% signal x. +% +% Input: +% x : finite length 1D or 2D signal (implicitely periodized) +% h : scaling filter +% L : number of levels. in case of a 1D signal length(x) must be +% divisible by 2^L; in case of a 2D signal the row and the +% column dimension must be divisible by 2^L. +% +% see also: midwt, mrdwt, mirdwt +*/ + +#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__ +MDWT(double *x, uintptr_t m, uintptr_t n, double *h, uintptr_t lh, uintptr_t L, double *y) +#else +MDWT(x, m, n, h, lh, L, y) +double *x, *h, *y; +uintptr_t m, n, lh, L; +#endif +{ + double *h0, *h1, *ydummyl, *ydummyh, *xdummy; + int *prob; + uintptr_t i, j; + uintptr_t actual_L, actual_m, actual_n, r_o_a, c_o_a, ir, ic, lhm1; + + xdummy = (double *)(uintptr_t)mxCalloc(max(m,n)+lh-1,sizeof(double)); + ydummyl =(double *) (uintptr_t)mxCalloc(max(m,n),sizeof(double)); + ydummyh = (double *)(uintptr_t)mxCalloc(max(m,n),sizeof(double)); + h0 =(double *)(uintptr_t)mxCalloc(lh,sizeof(double)); + h1 = (double *)(uintptr_t)mxCalloc(lh,sizeof(double)); + + + /* analysis lowpass and highpass */ + if (n==1){ + n = m; + m = 1; + } + for (i=0; i<lh; i++){ + h0[i] = h[lh-i-1]; + h1[i] =h[i]; + } + for (i=0; i<lh; i+=2) + h1[i] = -h1[i]; + + lhm1 = lh - 1; + actual_m = 2*m; + actual_n = 2*n; + + /* main loop */ + for (actual_L=1; actual_L <= L; actual_L++){ + if (m==1) + actual_m = 1; + else{ + actual_m = actual_m/2; + r_o_a = actual_m/2; + } + actual_n = actual_n/2; + c_o_a = actual_n/2; + + /* go by rows */ + for (ir=0; ir<actual_m; ir++){ /* loop over rows */ + /* store in dummy variable */ + for (i=0; i<actual_n; i++) + if (actual_L==1) + xdummy[i] = mat(x, ir, i); + else + xdummy[i] = mat(y, ir, i); + /* perform filtering lowpass and highpass*/ + fpsconv(xdummy, actual_n, h0, h1, lhm1, ydummyl, ydummyh); + /* restore dummy variables in matrices */ + ic = c_o_a; + for (i=0; i<c_o_a; i++){ + mat(y, ir, i) = ydummyl[i]; + mat(y, ir, ic++) = ydummyh[i]; + } + } + + /* go by columns in case of a 2D signal*/ + if (m>1){ + for (ic=0; ic<actual_n; ic++){ /* loop over column */ + /* store in dummy variables */ + for (i=0; i<actual_m; i++) + xdummy[i] = mat(y, i, ic); + /* perform filtering lowpass and highpass*/ + fpsconv(xdummy, actual_m, h0, h1, lhm1, ydummyl, ydummyh); + /* restore dummy variables in matrix */ + ir = r_o_a; + for (i=0; i<r_o_a; i++){ + mat(y, i, ic) = ydummyl[i]; + mat(y, ir++, ic) = ydummyh[i]; + } + } + } + } +} + +#ifdef __STDC__ +fpsconv(double *x_in, uintptr_t lx, double *h0, double *h1, uintptr_t lhm1, + double *x_outl, double *x_outh) +#else +fpsconv(x_in, lx, h0, h1, lhm1, x_outl, x_outh) +double *x_in, *h0, *h1, *x_outl, *x_outh; +uintptr_t lx, lhm1; +#endif + +{ + uintptr_t i, j, ind; + double x0, x1; + + for (i=lx; i < lx+lhm1; i++) + x_in[i] = *(x_in+(i-lx)); + ind = 0; + for (i=0; i<(lx); i+=2){ + x0 = 0; + x1 = 0; + for (j=0; j<=lhm1; j++){ + x0 = x0 + x_in[i+j]*h0[lhm1-j]; + x1 = x1 + x_in[i+j]*h1[lhm1-j]; + } + x_outl[ind] = x0; + x_outh[ind++] = x1; + } +}
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/util/Rice Wavelet Toolbox/midwt.c Fri Mar 25 15:27:33 2011 +0000 @@ -0,0 +1,226 @@ +/* +File Name: midwt.c +Last Modification Date: 06/14/95 12:55:58 +Current Version: midwt.c 1.4 +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, *y, *Lf, *Lr; + intptr_t m, n, h_col, h_row, lh, L, i, po2, j; + double mtest, ntest; + + /* check for correct # of input variables */ + if (nrhs>3){ + mexErrMsgTxt("There are at most 3 input parameters allowed!"); + return; + } + if (nrhs<2){ + mexErrMsgTxt("There are at least 2 input parameters required!"); + return; + } + y = mxGetPr(prhs[0]); + n = mxGetN(prhs[0]); + m = mxGetM(prhs[0]); + h = mxGetPr(prhs[1]); + h_col = mxGetN(prhs[1]); + h_row = mxGetM(prhs[1]); + if (h_col>h_row) + lh = h_col; + else + lh = h_row; + if (nrhs == 3){ + L = (intptr_t) *mxGetPr(prhs[2]); + 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 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; + MIDWT(x, m, n, h, lh, L, y); +} + +#define mat(a, i, j) (*(a + (m*(j)+i))) /* macro for matrix indices */ + + +#ifdef __STDC__ +MIDWT(double *x, intptr_t m, intptr_t n, double *h, intptr_t lh, intptr_t L, double *y) +#else +MIDWT(x, m, n, h, lh, L, y) +double *x, *h, *y; +intptr_t m, n, lh, L; +#endif +{ + double *g0, *g1, *ydummyl, *ydummyh, *xdummy; + intptr_t i, j; + intptr_t actual_L, actual_m, actual_n, r_o_a, c_o_a, ir, ic, lhm1, lhhm1, sample_f; + xdummy = (double *)mxCalloc(max(m,n),sizeof(double)); + ydummyl = (double *)mxCalloc(max(m,n)+lh/2-1,sizeof(double)); + ydummyh = (double *)(intptr_t)mxCalloc(max(m,n)+lh/2-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; + } + /* synthesis lowpass and highpass */ + for (i=0; i<lh; i++){ + g0[i] = h[i]; + g1[i] = h[lh-i-1]; + } + for (i=1; i<=lh; i+=2) + g1[i] = -g1[i]; + + lhm1 = lh - 1; + lhhm1 = lh/2 - 1; + /* 2^L */ + sample_f = 1; + for (i=1; i<L; i++) + sample_f = sample_f*2; + + if (m>1) + actual_m = m/sample_f; + else + actual_m = 1; + actual_n = n/sample_f; + + for (i=0; i<(m*n); i++) + x[i] = y[i]; + + /* main loop */ + for (actual_L=L; actual_L >= 1; actual_L--){ + r_o_a = actual_m/2; + c_o_a = actual_n/2; + + /* go by columns in case of a 2D signal*/ + if (m>1){ + for (ic=0; ic<actual_n; ic++){ /* loop over column */ + /* store in dummy variables */ + ir = r_o_a; + for (i=0; i<r_o_a; i++){ + ydummyl[i+lhhm1] = mat(x, i, ic); + ydummyh[i+lhhm1] = mat(x, ir++, ic); + } + /* perform filtering lowpass and highpass*/ + bpsconv(xdummy, r_o_a, g0, g1, lhm1, lhhm1, ydummyl, ydummyh); + /* restore dummy variables in matrix */ + for (i=0; i<actual_m; i++) + mat(x, i, ic) = xdummy[i]; + } + } + /* go by rows */ + for (ir=0; ir<actual_m; ir++){ /* loop over rows */ + /* store in dummy variable */ + ic = c_o_a; + for (i=0; i<c_o_a; i++){ + ydummyl[i+lhhm1] = mat(x, ir, i); + ydummyh[i+lhhm1] = mat(x, ir, ic++); + } + /* perform filtering lowpass and highpass*/ + bpsconv(xdummy, c_o_a, g0, g1, lhm1, lhhm1, ydummyl, ydummyh); + /* restore dummy variables in matrices */ + for (i=0; i<actual_n; i++) + mat(x, ir, i) = xdummy[i]; + } + if (m==1) + actual_m = 1; + else + actual_m = actual_m*2; + actual_n = actual_n*2; + } +} + +#ifdef __STDC__ +bpsconv(double *x_out, intptr_t lx, double *g0, double *g1, intptr_t lhm1, + intptr_t lhhm1, double *x_inl, double *x_inh) +#else +bpsconv(x_out, lx, g0, g1, lhm1, lhhm1, x_inl, x_inh) +double *x_inl, *x_inh, *g0, *g1, *x_out; +intptr_t lx, lhm1, lhhm1; +#endif +{ + intptr_t i, j, ind, tj; + double x0, x1; + + for (i=lhhm1-1; i > -1; i--){ + x_inl[i] = x_inl[lx+i]; + x_inh[i] = x_inh[lx+i]; + } + ind = 0; + for (i=0; i<(lx); i++){ + x0 = 0; + x1 = 0; + tj = -2; + for (j=0; j<=lhhm1; j++){ + tj+=2; + x0 = x0 + x_inl[i+j]*g0[lhm1-1-tj] + x_inh[i+j]*g1[lhm1-1-tj] ; + x1 = x1 + x_inl[i+j]*g0[lhm1-tj] + x_inh[i+j]*g1[lhm1-tj] ; + } + x_out[ind++] = x0; + x_out[ind++] = x1; + } +}
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/util/Rice Wavelet Toolbox/midwt.mhelp Fri Mar 25 15:27:33 2011 +0000 @@ -0,0 +1,93 @@ +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. +% +% Input: +% y : finite length 1D or 2D input signal (implicitly periodized) +% (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. +% +% Output: +% x : periodic reconstructed signal +% L : number of decomposition levels +% +% 1D Example: +% xin = makesig('LinChirp',8); +% h = daubcqf(4,'min'); +% L = 1; +% [y,L] = mdwt(xin,h,L); +% [x,L] = midwt(y,h,L) +% +% 1D Example's output: +% +% x = 0.0491 0.1951 0.4276 0.7071 0.9415 0.9808 0.6716 0.0000 +% L = 1 +% +% See also: mdwt, mrdwt, mirdwt +% + +% +% +%File Name: midwt.m +%Last Modification Date: 08/07/95 15:13:52 +%Current Version: midwt.m 2.4 +%File Creation Date: Wed Oct 19 10:51:58 1994 +%Author: Markus Lang <lang@jazz.rice.edu> +% +%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 techtran@rice.edu or (713) 348-6173 +% +%Change History: +% +%Modification #1 +%Mon Aug 7 11:52:33 CDT 1995 +%Rebecca Hindman <hindman@ece.rice.edu> +%Added L to function line so that it can be displayed as an output +% +%Thu Mar 2 13:07:11 CDT 2000 +%Ramesh Neelamani<neelsh@ece.rice.edu> +%Revamped the help file +% + + + +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/util/Rice Wavelet Toolbox/midwt_r.c Fri Mar 25 15:27:33 2011 +0000 @@ -0,0 +1,193 @@ +/* +File Name: MIDWT.c +Last Modification Date: 06/14/95 13:01:15 +Current Version: MIDWT.c 2.4 +File Creation Date: Wed Oct 12 08:44:43 1994 +Author: Markus Lang <lang@jazz.rice.edu> + +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 techtran@rice.edu or (713) 348-6173 + +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 <odegard@ece.rice.edu> Wed Jun 14 1995 + + Fix minor bug to allow maximum number of levels + +decription of the matlab call: +%y = midwt(x,h,L); +% +% function computes the inverse discrete wavelet transform y for a 1D or 2D +% input signal x. +% +% Input: +% x : finite length 1D or 2D input signal (implicitely periodized) +% h : scaling filter +% L : number of levels. in case of a 1D signal length(x) must be +% divisible by 2^L; in case of a 2D signal the row and the +% column dimension must be divisible by 2^L. +% +% see also: mdwt, mrdwt, mirdwt + +*/ + +#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__ +MIDWT(double *x, uintptr_t m, uintptr_t n, double *h, uintptr_t lh, uintptr_t L, double *y) +#else +MIDWT(x, m, n, h, lh, L, y) +double *x, *h, *y; +uintptr_t m, n, lh, L; +#endif +{ + double *g0, *g1, *ydummyl, *ydummyh, *xdummy; + uintptr_t i, j; + uintptr_t actual_L, actual_m, actual_n, r_o_a, c_o_a, ir, ic, lhm1, lhhm1, sample_f; + xdummy = (double *)(uintptr_t)mxCalloc(max(m,n),sizeof(double)); + ydummyl = (double *)(uintptr_t)mxCalloc(max(m,n)+lh/2-1,sizeof(double)); + ydummyh = (double *)(uintptr_t)mxCalloc(max(m,n)+lh/2-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; + } + /* synthesis lowpass and highpass */ + for (i=0; i<lh; i++){ + g0[i] = h[i]; + g1[i] = h[lh-i-1]; + } + for (i=1; i<=lh; i+=2) + g1[i] = -g1[i]; + + lhm1 = lh - 1; + lhhm1 = lh/2 - 1; + /* 2^L */ + sample_f = 1; + for (i=1; i<L; i++) + sample_f = sample_f*2; + + if (m>1) + actual_m = m/sample_f; + else + actual_m = 1; + actual_n = n/sample_f; + + for (i=0; i<(m*n); i++) + x[i] = y[i]; + + /* main loop */ + for (actual_L=L; actual_L >= 1; actual_L--){ + r_o_a = actual_m/2; + c_o_a = actual_n/2; + + /* go by columns in case of a 2D signal*/ + if (m>1){ + for (ic=0; ic<actual_n; ic++){ /* loop over column */ + /* store in dummy variables */ + ir = r_o_a; + for (i=0; i<r_o_a; i++){ + ydummyl[i+lhhm1] = mat(x, i, ic); + ydummyh[i+lhhm1] = mat(x, ir++, ic); + } + /* perform filtering lowpass and highpass*/ + bpsconv(xdummy, r_o_a, g0, g1, lhm1, lhhm1, ydummyl, ydummyh); + /* restore dummy variables in matrix */ + for (i=0; i<actual_m; i++) + mat(x, i, ic) = xdummy[i]; + } + } + /* go by rows */ + for (ir=0; ir<actual_m; ir++){ /* loop over rows */ + /* store in dummy variable */ + ic = c_o_a; + for (i=0; i<c_o_a; i++){ + ydummyl[i+lhhm1] = mat(x, ir, i); + ydummyh[i+lhhm1] = mat(x, ir, ic++); + } + /* perform filtering lowpass and highpass*/ + bpsconv(xdummy, c_o_a, g0, g1, lhm1, lhhm1, ydummyl, ydummyh); + /* restore dummy variables in matrices */ + for (i=0; i<actual_n; i++) + mat(x, ir, i) = xdummy[i]; + } + if (m==1) + actual_m = 1; + else + actual_m = actual_m*2; + actual_n = actual_n*2; + } +} + +#ifdef __STDC__ +bpsconv(double *x_out, uintptr_t lx, double *g0, double *g1, uintptr_t lhm1, + uintptr_t lhhm1, double *x_inl, double *x_inh) +#else +bpsconv(x_out, lx, g0, g1, lhm1, lhhm1, x_inl, x_inh) +double *x_inl, *x_inh, *g0, *g1, *x_out; +uintptr_t lx, lhm1, lhhm1; +#endif +{ + uintptr_t i, j, ind, tj; + double x0, x1; + + for (i=lhhm1-1; i > -1; i--){ + x_inl[i] = x_inl[lx+i]; + x_inh[i] = x_inh[lx+i]; + } + ind = 0; + for (i=0; i<(lx); i++){ + x0 = 0; + x1 = 0; + tj = -2; + for (j=0; j<=lhhm1; j++){ + tj+=2; + x0 = x0 + x_inl[i+j]*g0[lhm1-1-tj] + x_inh[i+j]*g1[lhm1-1-tj] ; + x1 = x1 + x_inl[i+j]*g0[lhm1-tj] + x_inh[i+j]*g1[lhm1-tj] ; + } + x_out[ind++] = x0; + x_out[ind++] = x1; + } +}
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/util/Rice Wavelet Toolbox/mirdwt.c Fri Mar 25 15:27:33 2011 +0000 @@ -0,0 +1,261 @@ +/* +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; + } +}
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/util/Rice Wavelet Toolbox/mirdwt.mhelp Fri Mar 25 15:27:33 2011 +0000 @@ -0,0 +1,138 @@ +function [x,L] = mirdwt(yl,yh,h,L); +% function [x,L] = mirdwt(yl,yh,h,L); +% +% Function computes the inverse redundant discrete wavelet +% transform x for a 1D or 2D input signal. (Redundant means here +% that the sub-sampling 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 the 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 the case of a 1D signal, +% length(yl) 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. +% +% Output: +% x : finite length 1D or 2D signal +% L : number of levels +% +% HERE'S AN EASY WAY TO RUN THE EXAMPLES: +% Cut-and-paste the example you want to run to a new file +% called ex.m, for example. Delete out the % at the beginning +% of each line in ex.m (Can use search-and-replace in your editor +% to replace it with a space). Type 'ex' in matlab and hit return. +% +% +% Example 1: +% xin = makesig('Leopold',8); +% h = daubcqf(4,'min'); +% L = 1; +% [yl,yh,L] = mrdwt(xin,h,L); +% [x,L] = mirdwt(yl,yh,h,L) +% x = 0.0000 1.0000 0.0000 -0.0000 0 0 0 -0.0000 +% L = 1 +% +% Example 2: +% load lena; +% h = daubcqf(4,'min'); +% L = 2; +% [ll_lev2,yh,L] = mrdwt(lena,h,L); % lena is a 256x256 matrix +% N = 256; +% lh_lev1 = yh(:,1:N); +% hl_lev1 = yh(:,N+1:2*N); +% hh_lev1 = yh(:,2*N+1:3*N); +% lh_lev2 = yh(:,3*N+1:4*N); +% hl_lev2 = yh(:,4*N+1:5*N); +% hh_lev2 = yh(:,5*N+1:6*N); +% figure; colormap(gray); imagesc(lena); title('Original Image'); +% figure; colormap(gray); imagesc(ll_lev2); title('LL Level 2'); +% figure; colormap(gray); imagesc(hh_lev2); title('HH Level 2'); +% figure; colormap(gray); imagesc(hl_lev2); title('HL Level 2'); +% figure; colormap(gray); imagesc(lh_lev2); title('LH Level 2'); +% figure; colormap(gray); imagesc(hh_lev1); title('HH Level 1'); +% figure; colormap(gray); imagesc(hl_lev2); title('HL Level 1'); +% figure; colormap(gray); imagesc(lh_lev2); title('LH Level 1'); +% [lena_Hat,L] = mirdwt(ll_lev2,yh,h,L); +% figure; colormap(gray); imagesc(lena_Hat); +% title('Reconstructed Image'); +% +% See also: mdwt, midwt, mrdwt +% +% Warning! min(size(yl))/2^L should be greater than length(h) +% + +%File Name: mirdwt.m +%Last Modification Date: 08/07/95 15:14:21 +%Current Version: mirdwt.m 2.4 +%File Creation Date: Wed Oct 19 10:51:58 1994 +%Author: Markus Lang <lang@jazz.rice.edu> +% +%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 techtran@rice.edu or (713) 348-6173 +% +%Change History: +% +%Modification #1 +%Mon Aug 7 15:09:51 CDT 1995 +%Rebecca Hindman <hindman@ece.rice.edu> +%Added L to function line so that it can be displayed as an output +% +%Modification #2 +%Thursday Mar 2 2000 +% Added Example 2 +% Felix Fernandes <felixf@rice.edu> + + + + + + + + + + + + + + + +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/util/Rice Wavelet Toolbox/mirdwt_r.c Fri Mar 25 15:27:33 2011 +0000 @@ -0,0 +1,219 @@ +/* +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 <lang@jazz.rice.edu> + +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 techtran@rice.edu or (713) 348-6173 + +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 <odegard@ece.rice.edu> 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; + } +}
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/util/Rice Wavelet Toolbox/mrdwt.c Fri Mar 25 15:27:33 2011 +0000 @@ -0,0 +1,242 @@ +/* +File Name: mrdwt.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, h_col, h_row, lh, L, i, po2, j; + double mtest, ntest; + + /* check for correct # of input variables */ + if (nrhs>3){ + mexErrMsgTxt("There are at most 3 input parameters allowed!"); + return; + } + if (nrhs<2){ + mexErrMsgTxt("There are at least 2 input parameters required!"); + return; + } + x = mxGetPr(prhs[0]); + n = mxGetN(prhs[0]); + m = mxGetM(prhs[0]); + h = mxGetPr(prhs[1]); + h_col = mxGetN(prhs[1]); + h_row = mxGetM(prhs[1]); + if (h_col>h_row) + lh = h_col; + else + lh = h_row; + if (nrhs == 3){ + L = (intptr_t) *mxGetPr(prhs[2]); + 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 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); + yl = mxGetPr(plhs[0]); + if (min(m,n) == 1) + plhs[1] = mxCreateDoubleMatrix(m,L*n,mxREAL); + else + plhs[1] = mxCreateDoubleMatrix(m,3*L*n,mxREAL); + yh = mxGetPr(plhs[1]); + plhs[2] = mxCreateDoubleMatrix(1,1,mxREAL); + Lr = mxGetPr(plhs[2]); + *Lr = L; + MRDWT(x, m, n, h, lh, L, yl, yh); +} +#define mat(a, i, j) (*(a + (m*(j)+i))) + + +#ifdef __STDC__ +MRDWT(double *x, intptr_t m, intptr_t n, double *h, intptr_t lh, intptr_t L, + double *yl, double *yh) +#else +MRDWT(x, m, n, h, lh, L, yl, yh) +double *x, *h, *yl, *yh; +intptr_t m, n, lh, L; +#endif +{ + double *tmp; + double *h0, *h1, *ydummyll, *ydummylh, *ydummyhl; + double *ydummyhh, *xdummyl , *xdummyh; + long i, j; + intptr_t actual_L, actual_m, actual_n, c_o_a, ir, n_c, n_cb, n_c_o; + intptr_t ic, n_r, n_rb, n_r_o, c_o_a_p2n, sample_f; + xdummyl = (double *)(intptr_t)mxCalloc(max(m,n)+lh-1,sizeof(double)); + xdummyh = (double *)(intptr_t)mxCalloc(max(m,n)+lh-1,sizeof(double)); + ydummyll = (double *)(intptr_t)mxCalloc(max(m,n),sizeof(double)); + ydummylh = (double *)(intptr_t)mxCalloc(max(m,n),sizeof(double)); + ydummyhl = (double *)(intptr_t)mxCalloc(max(m,n),sizeof(double)); + ydummyhh = (double *)(intptr_t)mxCalloc(max(m,n),sizeof(double)); + h0 = (double *)(intptr_t)mxCalloc(lh,sizeof(double)); + h1 = (double *)(intptr_t)mxCalloc(lh,sizeof(double)); + + if (n==1){ + n = m; + m = 1; + } + /* analysis lowpass and highpass */ + for (i=0; i<lh; i++){ + h0[i] = h[lh-i-1]; + h1[i] =h[i]; + } + for (i=0; i<lh; i+=2) + h1[i] = -h1[i]; + + actual_m = 2*m; + actual_n = 2*n; + for (i=0; i<m*n; i++) + yl[i] = x[i]; + + /* main loop */ + sample_f = 1; + for (actual_L=1; actual_L <= L; actual_L++){ + actual_m = actual_m/2; + actual_n = actual_n/2; + /* 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 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; + xdummyl[i] = mat(yl, ir, ic); + } + /* perform filtering lowpass/highpass */ + fpconv(xdummyl, actual_n, h0, h1, 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(yl, ir, ic) = ydummyll[i]; + mat(yh, ir, c_o_a+ic) = ydummyhh[i]; + } + } + } + + /* 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; + xdummyl[i] = mat(yl, ir, ic); + xdummyh[i] = mat(yh, ir,c_o_a+ic); + } + /* perform filtering: first LL/LH, then HL/HH */ + fpconv(xdummyl, actual_m, h0, h1, lh, ydummyll, ydummylh); + fpconv(xdummyh, actual_m, h0, h1, lh, ydummyhl, ydummyhh); + /* restore dummy variables in matrices */ + ir = -sample_f + n_r; + for (i=0; i<actual_m; i++){ + ir = ir + sample_f; + mat(yl, ir, ic) = ydummyll[i]; + mat(yh, ir, c_o_a+ic) = ydummylh[i]; + mat(yh, ir,c_o_a+n+ic) = ydummyhl[i]; + mat(yh, ir, c_o_a_p2n+ic) = ydummyhh[i]; + } + } + } + } + sample_f = sample_f*2; + } +} + +#ifdef __STDC__ +fpconv(double *x_in, intptr_t lx, double *h0, double *h1, intptr_t lh, + double *x_outl, double *x_outh) +#else +fpconv(x_in, lx, h0, h1, lh, x_outl, x_outh) +double *x_in, *h0, *h1, *x_outl, *x_outh; +intptr_t lx, lh; +#endif +{ + intptr_t i, j; + double x0, x1; + + for (i=lx; i < lx+lh-1; i++) + x_in[i] = x_in[i-lx]; + for (i=0; i<lx; i++){ + x0 = 0; + x1 = 0; + for (j=0; j<lh; j++){ + x0 = x0 + x_in[j+i]*h0[lh-1-j]; + x1 = x1 + x_in[j+i]*h1[lh-1-j]; + } + x_outl[i] = x0; + x_outh[i] = x1; + } +}
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/util/Rice Wavelet Toolbox/mrdwt.mhelp Fri Mar 25 15:27:33 2011 +0000 @@ -0,0 +1,120 @@ +function [yl,yh,L] = mrdwt(x,h,L); +% [yl,yh,L] = mrdwt(x,h,L); +% +% Function computes the redundant discrete wavelet transform y +% for a 1D or 2D input signal. (Redundant means here that the +% sub-sampling after each stage is omitted.) yl contains the +% lowpass and yh the highpass components. In the 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: +% x : finite length 1D or 2D signal (implicitly periodized) +% h : scaling filter +% L : number of levels. In the case of a 1D +% 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. +% +% Output: +% yl : lowpass component +% yh : highpass components +% L : number of levels +% +% HERE'S AN EASY WAY TO RUN THE EXAMPLES: +% Cut-and-paste the example you want to run to a new file +% called ex.m, for example. Delete out the % at the beginning +% of each line in ex.m (Can use search-and-replace in your editor +% to replace it with a space). Type 'ex' in matlab and hit return. +% +% +% Example 1:: +% x = makesig('Leopold',8); +% h = daubcqf(4,'min'); +% L = 1; +% [yl,yh,L] = mrdwt(x,h,L) +% yl = 0.8365 0.4830 0 0 0 0 -0.1294 0.2241 +% yh = -0.2241 -0.1294 0 0 0 0 -0.4830 0.8365 +% L = 1 +% Example 2: +% load lena; +% h = daubcqf(4,'min'); +% L = 2; +% [ll_lev2,yh,L] = mrdwt(lena,h,L); % lena is a 256x256 matrix +% N = 256; +% lh_lev1 = yh(:,1:N); +% hl_lev1 = yh(:,N+1:2*N); +% hh_lev1 = yh(:,2*N+1:3*N); +% lh_lev2 = yh(:,3*N+1:4*N); +% hl_lev2 = yh(:,4*N+1:5*N); +% hh_lev2 = yh(:,5*N+1:6*N); +% figure; colormap(gray); imagesc(lena); title('Original Image'); +% figure; colormap(gray); imagesc(ll_lev2); title('LL Level 2'); +% figure; colormap(gray); imagesc(hh_lev2); title('HH Level 2'); +% figure; colormap(gray); imagesc(hl_lev2); title('HL Level 2'); +% figure; colormap(gray); imagesc(lh_lev2); title('LH Level 2'); +% figure; colormap(gray); imagesc(hh_lev1); title('HH Level 1'); +% figure; colormap(gray); imagesc(hl_lev2); title('HL Level 1'); +% figure; colormap(gray); imagesc(lh_lev2); title('LH Level 1'); +% +% See also: mdwt, midwt, mirdwt +% +% Warning! min(size(x))/2^L should be greater than length(h) +% + +%File Name: mrdwt.m +%Last Modification Date: 08/07/95 15:14:47 +%Current Version: mrdwt.m 2.4 +%File Creation Date: Wed Oct 19 10:51:58 1994 +%Author: Markus Lang <lang@jazz.rice.edu> +% +%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 techtran@rice.edu or (713) 348-6173 +% +%Change History: +% +%Modification #1 +%Mon Aug 7 14:26:26 CDT 1995 +%Rebecca Hindman <hindman@ece.rice.edu> +%Added L to function line so that it can be displayed as an output +% +%Thursday Mar 2 2000 +% Added Example 2 +% Felix Fernandes <felixf@rice.edu> + +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/util/Rice Wavelet Toolbox/mrdwt_r.c Fri Mar 25 15:27:33 2011 +0000 @@ -0,0 +1,213 @@ +/* +File Name: MRDWT.c +Last Modification Date: 09/21/95 15:42:59 +Current Version: MRDWT.c 2.4 +File Creation Date: Wed Oct 12 08:44:43 1994 +Author: Markus Lang <lang@jazz.rice.edu> + +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 techtran@rice.edu or (713) 348-6173 + +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 <odegard@ece.rice.edu> Wed Jun 14 1995 + +MATLAB description: +%[yl,yh] = mrdwt(x,h,L); +% +% function computes the redundant discrete wavelet transform y for a 1D or +% 2D input signal . redundant means here that the subsampling after each +% stage is omitted. yl contains the lowpass and yl the highpass +% components. 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: +% x : finite length 1D or 2D signal (implicitely periodized) +% h : scaling filter +% L : number of levels. in case of a 1D signal length(x) 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: +% yl : lowpass component +% yh : highpass components +% +% see also: mdwt, midwt, mirdwt + + +*/ + +#include <math.h> +#include <stdio.h> +#include <inttypes.h> + +/*#define mat(a, i, j) (a[m*(j)+i]) */ +#define max(a, b) ((a) > (b) ? (a) : (b)) +#define mat(a, i, j) (*(a + (m*(j)+i))) + + +#ifdef __STDC__ +MRDWT(double *x, uintptr_t m, uintptr_t n, double *h, uintptr_t lh, uintptr_t L, + double *yl, double *yh) +#else +MRDWT(x, m, n, h, lh, L, yl, yh) +double *x, *h, *yl, *yh; +uintptr_t m, n, lh, L; +#endif +{ + double *tmp; + double *h0, *h1, *ydummyll, *ydummylh, *ydummyhl; + double *ydummyhh, *xdummyl , *xdummyh; + long i, j; + uintptr_t actual_L, actual_m, actual_n, c_o_a, ir, n_c, n_cb, n_c_o; + uintptr_t ic, n_r, n_rb, n_r_o, c_o_a_p2n, sample_f; + xdummyl = (double *)(uintptr_t)mxCalloc(max(m,n)+lh-1,sizeof(double)); + xdummyh = (double *)(uintptr_t)mxCalloc(max(m,n)+lh-1,sizeof(double)); + ydummyll = (double *)(uintptr_t)mxCalloc(max(m,n),sizeof(double)); + ydummylh = (double *)(uintptr_t)mxCalloc(max(m,n),sizeof(double)); + ydummyhl = (double *)(uintptr_t)mxCalloc(max(m,n),sizeof(double)); + ydummyhh = (double *)(uintptr_t)mxCalloc(max(m,n),sizeof(double)); + h0 = (double *)(uintptr_t)mxCalloc(lh,sizeof(double)); + h1 = (double *)(uintptr_t)mxCalloc(lh,sizeof(double)); + + if (n==1){ + n = m; + m = 1; + } + /* analysis lowpass and highpass */ + for (i=0; i<lh; i++){ + h0[i] = h[lh-i-1]; + h1[i] =h[i]; + } + for (i=0; i<lh; i+=2) + h1[i] = -h1[i]; + + actual_m = 2*m; + actual_n = 2*n; + for (i=0; i<m*n; i++) + yl[i] = x[i]; + + /* main loop */ + sample_f = 1; + for (actual_L=1; actual_L <= L; actual_L++){ + actual_m = actual_m/2; + actual_n = actual_n/2; + /* 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 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; + xdummyl[i] = mat(yl, ir, ic); + } + /* perform filtering lowpass/highpass */ + fpconv(xdummyl, actual_n, h0, h1, 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(yl, ir, ic) = ydummyll[i]; + mat(yh, ir, c_o_a+ic) = ydummyhh[i]; + } + } + } + + /* 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; + xdummyl[i] = mat(yl, ir, ic); + xdummyh[i] = mat(yh, ir,c_o_a+ic); + } + /* perform filtering: first LL/LH, then HL/HH */ + fpconv(xdummyl, actual_m, h0, h1, lh, ydummyll, ydummylh); + fpconv(xdummyh, actual_m, h0, h1, lh, ydummyhl, ydummyhh); + /* restore dummy variables in matrices */ + ir = -sample_f + n_r; + for (i=0; i<actual_m; i++){ + ir = ir + sample_f; + mat(yl, ir, ic) = ydummyll[i]; + mat(yh, ir, c_o_a+ic) = ydummylh[i]; + mat(yh, ir,c_o_a+n+ic) = ydummyhl[i]; + mat(yh, ir, c_o_a_p2n+ic) = ydummyhh[i]; + } + } + } + } + sample_f = sample_f*2; + } +} + +#ifdef __STDC__ +fpconv(double *x_in, uintptr_t lx, double *h0, double *h1, uintptr_t lh, + double *x_outl, double *x_outh) +#else +fpconv(x_in, lx, h0, h1, lh, x_outl, x_outh) +double *x_in, *h0, *h1, *x_outl, *x_outh; +uintptr_t lx, lh; +#endif +{ + uintptr_t i, j; + double x0, x1; + + for (i=lx; i < lx+lh-1; i++) + x_in[i] = x_in[i-lx]; + for (i=0; i<lx; i++){ + x0 = 0; + x1 = 0; + for (j=0; j<lh; j++){ + x0 = x0 + x_in[j+i]*h0[lh-1-j]; + x1 = x1 + x_in[j+i]*h1[lh-1-j]; + } + x_outl[i] = x0; + x_outh[i] = x1; + } +}
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/util/Rice Wavelet Toolbox/setopt.m Fri Mar 25 15:27:33 2011 +0000 @@ -0,0 +1,70 @@ +function option = setopt(opt_par,default); +% option = setopt(opt_par,default); +% +% SETOPT can modify a default option vector with user specified options. +% +% Input: +% opt_par : Users desired option vector +% default : Program default option vector +% +% Output: +% option : New option vector +% +% Example: +% opt_par = [1 2 3 4]; +% default = [1 1 1 1]; +% option = setopt(opt_par,default) +% option = 1 2 3 4 +% + +%File Name: setopt.m +%Last Modification Date: 10/26/95 13:06:18 +%Current Version: setopt.m 2.4 +%File Creation Date: Thu Nov 11 10:40:22 1993 +%Author: Jan Erik Odegard <odegard@ece.rice.edu> +% +%Copyright (c) 2000 RICE UNIVERSITY. All rights reserved. +%Created by Jan Erik Odegard, 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 techtran@rice.edu or (713) 348-6173 +% +%Change History: +% + +if (nargin < 2) + error('You need two inputs'); +end; +len = length(opt_par); +option = zeros(size(default)); +option(1:len) = opt_par(1:len); +option = option + (option == 0).*default;
--- a/util/SMALL_ImgDeNoiseResult.m Wed Mar 23 17:08:55 2011 +0000 +++ b/util/SMALL_ImgDeNoiseResult.m Fri Mar 25 15:27:33 2011 +0000 @@ -21,17 +21,17 @@ im=SMALL.Problem.Original; imnoise=SMALL.Problem.Noisy; -subplot(2, m, 1); imshow(im/maxval); +subplot(2, m, 1); imagesc(im/maxval);colormap(gray);axis off; axis image; % Set aspect ratio to obtain square pixels title('Original image'); -subplot(2,m,m+1); imshow(imnoise/maxval); +subplot(2,m,m+1); imagesc(imnoise/maxval);axis off; axis image; title(sprintf('Noisy image, PSNR = %.2fdB', SMALL.Problem.noisy_psnr )); for i=2:m - subplot(2, m, i); imagesc(SMALL.solver(i-1).reconstructed.Image/maxval);colormap(gray); + subplot(2, m, i); imagesc(SMALL.solver(i-1).reconstructed.Image/maxval);axis off; axis image; title(sprintf('%s Denoised image, PSNR: %.2f dB in %.2f s',... - SMALL.DL(i-1).name, SMALL.solver(i-1).reconstructed.psnr, SMALL.solver(i-1).time )); + SMALL.DL(i-1).name, SMALL.solver(i-1).reconstructed.psnr, SMALL.solver(i-1).time ),'Interpreter','none'); if strcmpi(SMALL.DL(i-1).name,'ksvds') D = kron(SMALL.Problem.basedict{2},SMALL.Problem.basedict{1})*SMALL.DL(i-1).D; else @@ -40,8 +40,8 @@ dictimg = showdict(D,SMALL.Problem.blocksize,... round(sqrt(size(D,2))),round(sqrt(size(D,2))),'lines','highcontrast'); - subplot(2,m,m+i);imagesc(dictimg);colormap(gray); + subplot(2,m,m+i);imagesc(dictimg);axis off; axis image; title(sprintf('%s dictionary in %.2f s',... - SMALL.DL(i-1).name, SMALL.DL(i-1).time)); + SMALL.DL(i-1).name, SMALL.DL(i-1).time),'Interpreter','none'); end \ No newline at end of file
--- a/util/SMALL_plot.m Wed Mar 23 17:08:55 2011 +0000 +++ b/util/SMALL_plot.m Fri Mar 25 15:27:33 2011 +0000 @@ -20,7 +20,7 @@ for i =1:m subplot(m,n, (i-1)*n+1); plot(1:length(SMALL.solver(i).solution), SMALL.solver(i).solution, 'b') - title(sprintf('%s(%s) in %.2f s', SMALL.solver(i).name, SMALL.solver(i).param,SMALL.solver(i).time)) + title(sprintf('%s(%s) in %.2f s', SMALL.solver(i).name, SMALL.solver(i).param,SMALL.solver(i).time),'Interpreter','none') xlabel('Coefficient') % Plot reconstructed signal against original