fwa2.m

This package is a free collection of Matlab routines for computing wave atom transforms in one, two

function c = fwa2(x,pat,tp)
% fwa2 - 2D forward wave atom transform
% -----------------
% INPUT
% --
% x is a real N-by-N matrix. N is a power of 2.
% --
% pat specifies the type of frequency partition which satsifies
% parabolic scaling relationship. pat can either be 'p' or 'q'.
% --
% tp is the type of tranform.
% 	'ortho': orthobasis
% 	'directional': real-valued frame with single oscillation direction
% 	'complex': complex-valued frame
% -----------------
% OUTPUT
% --
% c is a cell array which contains the wave atom coefficients. If
% tp=='ortho', then c{j}{m1,m2}(n1,n2) is the coefficient at scale j,
% frequency index (m1,m2) and spatial index (n1,n2). If
% tp=='directional', then c{j,d}{m1,m2}(n1,n2) with d=1,2 are the
% coefficients at scale j, frequency index (m1,m2) and spatial index
% (n1,n2). If tp=='complex', then c{j,d}{m1,m2)(n1,n2) with d=1,2,3,4
% are the coefficients at scale j, frequency index (m1,m2) and spatial
% index (n1,n2).
% -----------------
% Written by Lexing Ying and Laurent Demanet, 2007
  
  if( ismember(tp, {'ortho','directional','complex'})==0 | ismember(pat, {'p','q','u'})==0 )    error('wrong');  end
  
  if(strcmp(tp, 'ortho')==1)
    %---------------------------------------------------------
    N = size(x,1);
    H = N/2;
    lst = freq_pat(H,pat);
    %------------------
    f = fft2(x) / sqrt(prod(size(x)));
    A = N;
    c = cell(length(lst),1);
    %------------------
    for s=1:length(lst)
      nw = length(lst{s});
      c{s} = cell(nw,nw);
      for I=0:nw-1
        for J=0:nw-1
          if(lst{s}(I+1)==0 & lst{s}(J+1)==0)
            c{s}{I+1,J+1} = [];
          else
            B = 2^(s-1);
            D = 2*B;
            Ict = I*B;      Jct = J*B; %starting position in freq
            if(mod(I,2)==0)
              Ifm = Ict-2/3*B;        Ito = Ict+4/3*B;
            else
              Ifm = Ict-1/3*B;        Ito = Ict+5/3*B;
            end
            if(mod(J,2)==0)
              Jfm = Jct-2/3*B;        Jto = Jct+4/3*B;
            else
              Jfm = Jct-1/3*B;        Jto = Jct+5/3*B;
            end
            res = zeros(D,D);
            for id=0:1
              if(id==0)
                Idx = [ceil(Ifm):floor(Ito)];      Icf = kf_rt(Idx/B*pi, I);
              else
                Idx = [ceil(-Ito):floor(-Ifm)];      Icf = kf_lf(Idx/B*pi, I);
              end
              for jd=0:1
                if(jd==0)
                  Jdx = [ceil(Jfm):floor(Jto)];      Jcf = kf_rt(Jdx/B*pi, J);
                else
                  Jdx = [ceil(-Jto):floor(-Jfm)];      Jcf = kf_lf(Jdx/B*pi, J);
                end
                res(mod(Idx,D)+1,mod(Jdx,D)+1) = res(mod(Idx,D)+1,mod(Jdx,D)+1) + conj( Icf.'*Jcf ) .* f(mod(Idx,A)+1,mod(Jdx,A)+1);
              end
            end
            c{s}{I+1,J+1} = ifft2(res) * sqrt(prod(size(res)));
          end
        end
      end
    end
    
  elseif(strcmp(tp, 'directional')==1)
    %---------------------------------------------------------
    N = size(x,1);
    H = N/2;
    lst = freq_pat(H,pat);
    %------------------
    f = fft2(x) / sqrt(prod(size(x)));
    A = N;
    c1 = cell(length(lst),1);
    c2 = cell(length(lst),1);
    %------------------
    for s=1:length(lst)
      nw = length(lst{s});
      c1{s} = cell(nw,nw);
      c2{s} = cell(nw,nw);
      for I=0:nw-1
        for J=0:nw-1
          if(lst{s}(I+1)==0 & lst{s}(J+1)==0)
            c1{s}{I+1,J+1} = [];
            c2{s}{I+1,J+1} = [];
          else
            B = 2^(s-1);
            D = 2*B;
            res = zeros(D,D);
            Ict = I*B;      Jct = J*B; %starting position in freq
            if(mod(I,2)==0)
              Ifm = Ict-2/3*B;        Ito = Ict+4/3*B;
            else
              Ifm = Ict-1/3*B;        Ito = Ict+5/3*B;
            end
            if(mod(J,2)==0)
              Jfm = Jct-2/3*B;        Jto = Jct+4/3*B;
            else
              Jfm = Jct-1/3*B;        Jto = Jct+5/3*B;
            end

            res = zeros(D,D);
            Idx = [ceil(Ifm):floor(Ito)];      Icf = kf_rt(Idx/B*pi, I);
            Jdx = [ceil(Jfm):floor(Jto)];      Jcf = kf_rt(Jdx/B*pi, J);
            res(mod(Idx,D)+1,mod(Jdx,D)+1) = res(mod(Idx,D)+1,mod(Jdx,D)+1) + conj( Icf.'*Jcf ) .* f(mod(Idx,A)+1,mod(Jdx,A)+1);
            Idx = [ceil(-Ito):floor(-Ifm)];      Icf = kf_lf(Idx/B*pi, I);
            Jdx = [ceil(-Jto):floor(-Jfm)];      Jcf = kf_lf(Jdx/B*pi, J);
            res(mod(Idx,D)+1,mod(Jdx,D)+1) = res(mod(Idx,D)+1,mod(Jdx,D)+1) + conj( Icf.'*Jcf ) .* f(mod(Idx,A)+1,mod(Jdx,A)+1);
            c1{s}{I+1,J+1} = ifft2(res) * sqrt(prod(size(res)));
            
            res = zeros(D,D);
            Idx = [ceil(Ifm):floor(Ito)];      Icf = kf_rt(Idx/B*pi, I);
            Jdx = [ceil(-Jto):floor(-Jfm)];      Jcf = kf_lf(Jdx/B*pi, J);
            res(mod(Idx,D)+1,mod(Jdx,D)+1) = res(mod(Idx,D)+1,mod(Jdx,D)+1) + conj( Icf.'*Jcf ) .* f(mod(Idx,A)+1,mod(Jdx,A)+1);
            Idx = [ceil(-Ito):floor(-Ifm)];      Icf = kf_lf(Idx/B*pi, I);
            Jdx = [ceil(Jfm):floor(Jto)];      Jcf = kf_rt(Jdx/B*pi, J);
            res(mod(Idx,D)+1,mod(Jdx,D)+1) = res(mod(Idx,D)+1,mod(Jdx,D)+1) + conj( Icf.'*Jcf ) .* f(mod(Idx,A)+1,mod(Jdx,A)+1);
            c2{s}{I+1,J+1} = ifft2(res) * sqrt(prod(size(res)));
          end
        end
      end
    end
    c = [c1 c2];
    
  elseif(strcmp(tp, 'complex')==1)
    %---------------------------------------------------------
    N = size(x,1);
    H = N/2;
    lst = freq_pat(H,pat);
    %------------------
    f = fft2(x) / sqrt(prod(size(x)));
    A = N;
    c1 = cell(length(lst),1);
    c2 = cell(length(lst),1);
    c3 = cell(length(lst),1);
    c4 = cell(length(lst),1);
    %------------------
    for s=1:length(lst)
      nw = length(lst{s});
      c1{s} = cell(nw,nw);
      c2{s} = cell(nw,nw);
      c3{s} = cell(nw,nw);
      c4{s} = cell(nw,nw);
      for I=0:nw-1
        for J=0:nw-1
          if(lst{s}(I+1)==0 & lst{s}(J+1)==0)
            c1{s}{I+1,J+1} = [];
            c2{s}{I+1,J+1} = [];
            c3{s}{I+1,J+1} = [];
            c4{s}{I+1,J+1} = [];
          else
            B = 2^(s-1);
            D = 2*B;
            res = zeros(D,D);
            Ict = I*B;      Jct = J*B; %starting position in freq
            if(mod(I,2)==0)
              Ifm = Ict-2/3*B;        Ito = Ict+4/3*B;
            else
              Ifm = Ict-1/3*B;        Ito = Ict+5/3*B;
            end
            if(mod(J,2)==0)
              Jfm = Jct-2/3*B;        Jto = Jct+4/3*B;
            else
              Jfm = Jct-1/3*B;        Jto = Jct+5/3*B;
            end

            res = zeros(D,D);
            Idx = [ceil(Ifm):floor(Ito)];      Icf = kf_rt(Idx/B*pi, I);
            Jdx = [ceil(Jfm):floor(Jto)];      Jcf = kf_rt(Jdx/B*pi, J);
            res(mod(Idx,D)+1,mod(Jdx,D)+1) = res(mod(Idx,D)+1,mod(Jdx,D)+1) + conj( Icf.'*Jcf ) .* f(mod(Idx,A)+1,mod(Jdx,A)+1);
            c1{s}{I+1,J+1} = ifft2(res) * sqrt(prod(size(res)));
            
            res = zeros(D,D);
            Idx = [ceil(-Ito):floor(-Ifm)];      Icf = kf_lf(Idx/B*pi, I);
            Jdx = [ceil(-Jto):floor(-Jfm)];      Jcf = kf_lf(Jdx/B*pi, J);
            res(mod(Idx,D)+1,mod(Jdx,D)+1) = res(mod(Idx,D)+1,mod(Jdx,D)+1) + conj( Icf.'*Jcf ) .* f(mod(Idx,A)+1,mod(Jdx,A)+1);
            c2{s}{I+1,J+1} = ifft2(res) * sqrt(prod(size(res)));
            
            res = zeros(D,D);
            Idx = [ceil(Ifm):floor(Ito)];      Icf = kf_rt(Idx/B*pi, I);
            Jdx = [ceil(-Jto):floor(-Jfm)];      Jcf = kf_lf(Jdx/B*pi, J);
            res(mod(Idx,D)+1,mod(Jdx,D)+1) = res(mod(Idx,D)+1,mod(Jdx,D)+1) + conj( Icf.'*Jcf ) .* f(mod(Idx,A)+1,mod(Jdx,A)+1);
            c3{s}{I+1,J+1} = ifft2(res) * sqrt(prod(size(res)));
            
            res = zeros(D,D);
            Idx = [ceil(-Ito):floor(-Ifm)];      Icf = kf_lf(Idx/B*pi, I);
            Jdx = [ceil(Jfm):floor(Jto)];      Jcf = kf_rt(Jdx/B*pi, J);
            res(mod(Idx,D)+1,mod(Jdx,D)+1) = res(mod(Idx,D)+1,mod(Jdx,D)+1) + conj( Icf.'*Jcf ) .* f(mod(Idx,A)+1,mod(Jdx,A)+1);
            c4{s}{I+1,J+1} = ifft2(res) * sqrt(prod(size(res)));
          end
        end
      end
    end
    c = [c1 c2 c3 c4];
  end