function varargout=angularD4radialHaarWT(x,levels,precon,meth,tipe,mr)
% f=angularD4radialHaarWT(x,levels,precon,meth,tipe,mr)
%
% Performs D4 wavelet transform in the angular direction (i.e. in the
% first and second index), and Haar wavelet transform in the radial
% direction (i.e. in the third index).
%
% INPUT:
%
% x            The three-dimensional array
% levels=[n1,n2,n3]
%              The number of levels in the three directions
%              n1=0, n2=0, n3=0 is the identity if precon=[0 0 0]
% precon       Array of length 3 identifying preconditioning  
%      [1 1 1] Precondition (NOT an orthogonal transform) so a sequence of
%              ones and a linearly increasing sequence are mapped to a
%              sequence of zeroes in the wavelet bands 
%      [0 0 0] Don't precondition (YES, this is an orthogonal transform)
%              so a sequence of ones and a linearly increasing sequence are
%              not mapped to all zeroes at the edges in the wavelet bands,
%              i.e. the first and last two coefficients are not zero
% meth  'pixel'     Angular reduction by numerical integration
%       'wavelets'  Angular reduction by wavelet transformation
% tipe         'forward'|'inverse'|'transpose'|'inversetranspose', always of
%              the same meth in both indices
% mr         0 The "original" improper, dimensionally sequential transform
%            1 The "proper" multiresolution transform [default]
%
% OUTPUT:
% 
% f            The size(x) matrix with the transform coefficients; scaling
%              coefficients in front, followed by wavelet coefficients. 
%
% EXAMPLE:
%
% angularD4radialHaarWT('demo1')
% angularD4radialHaarWT('demo2')
%
% Last modified by fjsimons-at-alum.mit.edu, Nov 7, 2011
% Last modified by yanhuay-at-princeton.edu, Nov 16, 2011

if ~isstr(x)
  % Get the levels of the decomposition
  defval('levels',[3 3 3])
  defval('meth', 'pixel')
  defval('tipe','forward')
  defval('mr',1)
  defval('precon',[1 1 1])
  
  % Determine the original size of the input
  [neta,nxi,nrad]=size(x);
  
  % Radial reduction from neta*nxi*129 to neta*nxi*37
  if nrad~=37
    x=radialpart129to37(x);
  end

  % Angular reduction from 513*513*nrad to 128*128*nrad
  if neta~=128 || nxi~=128
    x=angularpart513to128(x,meth);
  end
  
  % Haar wavelet transform in radial direction 
  levelHaar=levels(3);
  preconHaar=precon(3);
  x=radialHaarWT(x,levelHaar,preconHaar,tipe);
  
  % D4 wavelet in angular direction
  levelD4=levels(1:2);
  % If the first reduction is by wavelets, no more need for
  % preconditioning, see IL document (2/7/2009) page 7
  preconD4=[1 1]*strcmp(meth,'pixel');
  % Do the angular transform
  x=angularD4WT(x,levelD4,preconD4,tipe,mr);
  
  % Prepare the output
  varns={x};
  varargout=varns(1:nargout);
  
elseif strcmp(x,'demo1')
  disp('Check whether inverse is inverse of forward:')
  num1=513;
  num2=513;
  numr=129;
  for lev=3:3
    x=randn(num1,num2,numr);
    fp=angularD4radialHaarWT(x,[lev lev lev],[],'pixel','forward');
    xp=angularD4radialHaarWT(fp,[lev lev lev],[],'pixel','inverse');

    fw=angularD4radialHaarWT(x,[lev lev lev],[],'wavelets','forward');
    xw=angularD4radialHaarWT(fw,[lev lev lev],[],'wavelets','inverse');
    %Test the error 
    % need to consider whether we could inverse reduction xp/xw back to size of x  
    if numr~=37
      x37=radialpart129to37(x);
    else
      x37=x;
    end
    if num1~=128 || num2~=128
      xsmallp=angularpart513to128(x37,'pixel');
      xsmallw=angularpart513to128(x37,'wavelets');
    else
      xsmallp=x37;
      xsmallw=x37;
    end
    
    errp=max(max(max(abs(xp-xsmallp))));
    disp(['Levels=' num2str(lev) ' The maximal error for pixel method is ' num2str(errp)]);
    if(errp<1e-6)
      disp('Wavelet transform can be completely constructed')
    else
      disp('Wavelet transform can not be completely constructed')
    end
    disp(' ')
    
    errw=max(max(max(abs(xw-xsmallw))));
    disp(['Levels=' num2str(lev) ' The maximal error for wavelet method is ' num2str(errw)]);
    if(errw<1e-6)
      disp('Wavelet transform can be completely constructed')
    else
      disp('Wavelet transform can not be completely constructed')
    end
    disp(' ')
  end
  
elseif strcmp(x,'demo2') 
  disp('Check the inner product:')
  num1=513;
  num2=513;
  numr=129;
  for lev=3
    % Construct a fine kernel
    kernel=ones(num1,num2,numr);
    % By pixel method
    kwav1=angularD4radialHaarWT(kernel,[lev lev lev],[],'pixel','inversetranspose');
    model=ones(128,128,37); 
    modelwav1=angularD4radialHaarWT(model,[lev lev lev],[],'pixel','forward');
    p1=dotproduct3d(modelwav1,kwav1);
    
    % By wavelet
      % reduction
    kwav2=angularD4radialHaarWT(kernel,[lev lev lev],[],'wavelets','inversetranspose');
    model=ones(512,512,37); 
    model=angularD4WT(model,[2 2],[1 1],'forward');
    model=model(1:128,1:128,:);
    modelwav2=angularD4radialHaarWT(model,[lev lev lev],[],'wavelets','forward');
    p2=dotproduct3d(modelwav2,kwav2); 
    errpw=max(max(max(abs((p1-p2)/p1))));
    disp(['kernel=constant. Constant model. Integral err is:' num2str(errpw)]);
    disp(' ')   
  end
end