own dwt

Author: nath@W740SU

cshift.m

@@ -0,0 +1,28 @@
+function shifted = cshift(orig, shift, dir)
+
+% 2D Circular Shift
+% 
+% USAGE:
+%    y = cshift2D(x, m)
+% INPUT:
+%    x - M by N array
+%    m - amount of shift
+% OUTPUT:
+%    y - matrix x will be shifed by m samples down
+%
+% WAVELET SOFTWARE AT POLYTECHNIC UNIVERSITY, BROOKLYN, NY
+% http://taco.poly.edu/WaveletSoftware/
+
+[nr, nc] = size(orig);
+if dir == 1
+    n = 0:nr-1;
+    n = mod(n-shift, nr);
+    shifted = orig(n+1,:);
+elseif dir == 2
+    n = 0:nc-1;
+    n = mod(n-shift, nc);
+    shifted = orig(:,n+1);
+end
+end
+
+

main.m

@@ -1,7 +1,7 @@
 clear all; close all;
 
 MIN_THRESH = 4;
-LEVEL = 5;
+LEVEL = 6;
 LOGRESCALE = 0;
 IMAGE = 'castle.png';
 OIMAGE = 'castle';

singledec.m

@@ -1,16 +1,46 @@
 function [A,H,V,D] = singledec(img,LPF,HPF)
-    if size(LPF,1) > 1
-        LPF = LPF';
-    end
-    if size(HPF,1) > 1
-        HPF = HPF';
-    end
-    n = size(img,3);
-    for i =1:n
-        [a_, h_, v_, d_] = dwt2(img(:,:,i), LPF, HPF, 'mode','per');
-        A(:,:,i) = a_;
-        H(:,:,i) = h_;
-        V(:,:,i) = v_;
-        D(:,:,i) = d_;
-    end
+if size(LPF,1) > 1
+    LPF = LPF';
+end
+if size(HPF,1) > 1
+    HPF = HPF';
+end
+L = length(LPF)/2;
+[nr, nc, nl] = size(img);
+for layer = 1:nl
+    imglayer = cshift(img(:,:,layer), -L, 2);
+
+    convresult = conv2(imglayer, LPF);
+    hzlo = convresult(:,1:2:end);
+    hzlo(:, 1:L) = hzlo(:, 1:L) + hzlo(:,(1:L)+nc/2);
+    hzlo = hzlo(:,1:nc/2);
+
+    convresult = conv2(imglayer, HPF);
+    hzhi = convresult(:,1:2:end);
+    hzhi(:, 1:L) = hzhi(:, 1:L) + hzhi(:,(1:L)+nc/2);
+    hzhi = hzhi(:,1:nc/2);
+
+    hzlo = cshift(hzlo, -L, 1);
+    hzhi = cshift(hzhi, -L, 1);
+
+    convresult = conv2(hzlo, LPF');
+    alayer = convresult(1:2:end,:);
+    alayer(1:L, :) = alayer(1:L, :) + alayer((1:L)+nr/2, :);
+    A(:,:,layer) = alayer(1:nr/2,:);
+
+    convresult = conv2(hzhi, LPF');
+    hlayer = convresult(1:2:end,:);
+    hlayer(1:L, :) = hlayer(1:L, :) + hlayer((1:L)+nr/2, :);
+    H(:,:,layer) = hlayer(1:nr/2,:);
+
+    convresult = conv2(hzlo, HPF');
+    vlayer = convresult(1:2:end,:);
+    vlayer(1:L, :) = vlayer(1:L, :) + vlayer((1:L)+nr/2, :);
+    V(:,:,layer) = vlayer(1:nr/2,:);
+
+    convresult = conv2(hzhi, HPF');
+    dlayer = convresult(1:2:end,:);
+    dlayer(1:L, :) = dlayer(1:L, :) + dlayer((1:L)+nr/2, :);
+    D(:,:,layer) = dlayer(1:nr/2,:);
+end
 end

singlerec.m

@@ -5,8 +5,38 @@
     if size(HPF,1) > 1
         HPF = HPF';
     end
-    o = size(A,3);
-    for i =1:o
-        img(:,:,i) = idwt2(A(:,:,i), H(:,:,i), V(:,:,i), D(:,:,i), LPF, HPF,'mode','per');
-    end
+    [nr, nc, nl] = size(A);
+for layer = 1:nl
+    aup = zeros(nr, 2*nc);
+    aup(:,1:2:end) = A(:,:,layer);
+    hup = zeros(nr, 2*nc);
+    hup(:,1:2:end) = H(:,:,layer);
+    vtlo = conv2(aup,LPF )+conv2(hup,HPF);
+    L = length(LPF);
+    vtlo(:, 1:L-2) = vtlo(:,1:L-2) + vtlo(:, 2*nc+(1:L-2));
+    vtlo = vtlo(:,1:2*nc);
+    vtlo = cshift(vtlo, 1-L/2, 2);
+
+    vup = zeros(nr, 2*nc);
+    vup(:,1:2:end) = V(:,:,layer);
+    dup = zeros(nr, 2*nc);
+    dup(:,1:2:end) = D(:,:,layer);
+    vthi = conv2(vup,LPF)+conv2(dup,HPF);
+    vthi(:, 1:L-2) = vthi(:,1:L-2) + vthi(:, 2*nc+(1:L-2));
+    vthi = vthi(:,1:2*nc);
+    vthi = cshift(vthi, 1-L/2, 2);
+
+    vtloup = zeros(2*nr,2*nc);
+    vtloup(1:2:end,:) = vtlo;
+    vthiup = zeros(2*nr,2*nc);
+    vthiup(1:2:end,:) = vthi;
+
+    layerout = conv2(vtloup,LPF')+conv2(vthiup,HPF');
+    layerout(1:L-2,:) = layerout(1:L-2,:) + layerout(2*nr+(1:L-2), :);
+    layerout = layerout(1:2*nr, :);
+    layerout = cshift(layerout, 1-L/2, 1);
+
+    img(:,:,layer) = layerout;
+end
+
 end