function v=cp_energy(y,NonlinId,SmoothId,NormId,f_0,Structure,DoDecFB,FeatureFN,FeatDec)
% CP_ENERGY  Local energy function for texture comparison experiment
%	v=cp_energy(y,NonlinId,SmoothId,NormId,[f_0_v f_0_h]) 
%	computes the local energy of the image y corresponding to y being a
%	subband with center frequencies f_0_v and f_0_h (range 0 to 0.5),
%	alternatively (for the averaging filters) with averaging by an f_0_v
%	by f_0_h mask, or by an f_0_v by f_0_h median filter.
%	
%	v=cp_energy(y,'tanh:0.2',SmoothId,NormId,[f_0_v f_0_h]) 
%	sets alpha=0.2 with the 'tanh' nonlinearity (default 0.25).
%	
%	v=cp_energy(y,NonlinId,SmoothId,NormId,[],alpha,Structure,...
%	  DoDecFB,FeatureFN,FeatDec)
%	computes the local energy for the features stored in the file
%	FeatureFN assuming the specified filter structure.
%	
%	Legal local energy parameters are:
%	NonlinId:	SmoothId:	NormId:
%	- ''		- ''		- ''
%	- 'tanh'	- 'Gauss'	- 'Sqrt'
%	- 'Sqr'  	- 'Average'	- 'Log'
%	- 'Abs'		- 'Median'		

% Default parameter values
if exist('f_0')  ~=1,     f_0=[];       end
if exist('Structure')~=1, Structure=[]; end

if isempty(Structure)
  % Nonlinearity
  if     strcmp(lower(NonlinId),lower(''))
    z=y;
  elseif idcmp(NonlinId,'tanh')
    if length(NonlinId)==4
      alpha=0.25;
    else 
      alpha=str2num(NonlinId(6:length(NonlinId))); 
    end
    z=abs(tanh(alpha*y));
  elseif strcmp(lower(NonlinId),lower('Sqr'))
    z=y.^2;
  elseif strcmp(lower(NonlinId),lower('Abs'))
    z=abs(y);
  else
    error(sprintf('Nonlinearity ID "%s" not supported',NonlinId));
  end

  % Smoothing
  if     strcmp(lower(SmoothId),lower(''))
    v=z;
  elseif strcmp(lower(SmoothId),lower('Gauss'))
    Sigma=1./(2*sqrt(2)*f_0); % 	Coefficient for the Gaussian filters
    v=gausslp(z,Sigma);
  elseif strcmp(lower(SmoothId),lower('Average'))
    M_w=f_0(1);
    N_w=f_0(2);
    v=avgsmth(z,M_w,N_w);
  elseif strcmp(lower(SmoothId),lower('Median'))
    v=medfilt2(z,f_0);
  else
    error(sprintf('Smoothing ID "%s" not supported',SmoothId));
  end

  % Normalization
  if     strcmp(lower(NormId),lower(''))
    v=v;
  elseif strcmp(lower(NormId),lower('Log'))
    v=log(v+1);
  elseif strcmp(lower(NormId),lower('Sqrt'))
    v=sqrt(v);
  else
    error(sprintf('Normalization ID "%s" not supported',NormId));
  end
  
else
  % ----------------------------------------------------------------
  % Structure not empty, i.e. processing output from qmf
  % ----------------------------------------------------------------
  % Paremeter check
  if DoDecFB, 
    if FeatDec>1
      fprintf(2,'WARNING: Will not decimate already decimated output\n');
    end
    FeatDec=1; 
  end
  
  % Nonlinearity
  if     strcmp(lower(NonlinId),lower(''))
    NonlinCmd='cat';
  elseif strcmp(lower(NonlinId),lower('tanh'))
    if DoDecFB
      NonlinCmd=sprintf('tanh -sn %d -alpha %f',Structure,alpha);
    else
      NonlinCmd=sprintf('tanh -alpha %f',alpha);
    end
  elseif strcmp(lower(NonlinId),lower('Sqr'))
    NonlinCmd='powerpix -p 2';
  elseif strcmp(lower(NonlinId),lower('Abs'))
    NonlinCmd='abspix';
  else
    error(sprintf('Nonlinearity ID "%s" not supported',NonlinId));
  end

  % Smoothing
  if     strcmp(lower(SmoothId),lower(''))
    SmoothCmd='cat';
  elseif strcmp(lower(SmoothId),lower('Gauss'))
    if DoDecFB
      SmoothCmd =sprintf('blurr -sn %d -d',Structure);
    else
      SmoothCmd =sprintf('blurr -sn %d',Structure);
    end
  elseif strcmp(lower(SmoothId),lower('Average'))
    error(sprintf('Average filtering not supported for this param. combin.'));
  else
    error(sprintf('Smoothing ID "%s" not supported',SmoothId));
  end

  % Normalization
  if     strcmp(lower(NormId),lower(''))
    NormCmd='cat';
  elseif strcmp(lower(NormId),lower('Log'))
    NormCmd='logimg -o 1';
  elseif strcmp(lower(NormId),lower('Sqrt'))
    NormCmd='powerpix -p 0.5';
  else
    error(sprintf('Normalization ID "%s" not supported',NormId));
  end
  
  % Perform command
  TmpFeatureFN=tempnam;
  eval(sprintf('!%s %s | %s | %s | decimate -N %d > %s',...
    NonlinCmd,FeatureFN,SmoothCmd,NormCmd,FeatDec,TmpFeatureFN));
  
  % Prepare for subsequent programs, i.e. break up into separate files.
  if DoDecFB
    eval(sprintf('!dec2seq -sn %d %s > %s',Structure,TmpFeatureFN,FeatureFN));
    eval(sprintf('!mv %s %s',FeatureFN,TmpFeatureFN));
  end
  N_features=hsize(TmpFeatureFN,'Frames');
  for n=0:N_features-1
    eval(sprintf('!subseq -f %d %d %s > %s_%03d',...
      n,n,TmpFeatureFN,FeatureFN,n));
  end

  % Clean
  eval(sprintf('!rm -f %s %s',FeatureFN,TmpFeatureFN));
end
  

