function [x,N_tex,r_xx,Components]=cp_image(ImageId,TrainNo,M,N)
% CP_IMAGE  Load image (texture technique comparison experiment)
%	[x,N_tex]=cp_image(ImageId) load the image, N_tex identifies the
%	number of textures in this image.
%	[x,N_tex,r_xx]=cp_image(ImageId,0) load reference image
%	[x,N_tex,r_xx]=cp_image(ImageId,k) load training image number k
%	(starting at 1)

% Initialization
if exist('TrainNo')~=1
  TrainNo=[];
end
x=[];
r_xx=[];
if exist('M')~=1, M=256; end
if exist('N')~=1, N=256; end

% Load/compute data
if isstr(ImageId) 
  % ---------------------------------------------------------------------
  % Load an image identified by its name
  % ---------------------------------------------------------------------
  if strcmp(ImageId,'Nat-5c'), 			% Brodatz
    Components=[77 84 55 53 24];
  elseif strcmp(ImageId,'Nat-16c'),
    Components=[3 4 5 6 9 21 24 29 32 33 54 55 57 68 77 84];
  elseif strcmp(ImageId,'Nat-10'),
    Components=[21 4 36 24 9 37 28 19 29 38];
  elseif strcmp(ImageId,'Nat-5v'), 		% VisTex (MIT)
    Components=[0 17 22 49 56]+2000;
  elseif strcmp(ImageId,'Nat-5v2'),
    Components=[9 16 19 25 33]+2000;
  elseif strcmp(ImageId,'Nat-5v3'),
    Components=[46 7 68 83 9]+2000;
  elseif strcmp(ImageId,'Nat-16v'),
    Components=[7 9 13 14 16 25 33 41 46 51 55 60 62 68 83 94]+2000;
  elseif strcmp(ImageId,'Nat-10v'),
    Components=[25  9 68 33 16 83 46 19 66 94]+2000;
  elseif strcmp(ImageId,'Nat-5m'),		% MeasTex
    Components=[0 5 18 49 36]+3000;
  elseif strcmp(ImageId,'Nat-10m'),
    Components=[52 31 1 56 30 42 54 63 50 12]+3000;
  elseif strcmp(ImageId,'GMRF'),		% Synthetic
    msu_texture_id
    Components=[MSU_GMRF1 MSU_GMRF2 MSU_GMRF3 MSU_GMRF4];
  else
    error(sprintf('Unsupported image identifier parameter "%s"',ImageId));
  end
  N_tex=length(Components);
  
  % FileName=[];
  if isempty(TrainNo) % 		Load entire image
    % FileName=sprintf('%s/datafiler/textures/%s.hips',getenv('HOME'),ImageId);
    if strcmp(ImageId,'Nat-10') ...
	  | strcmp(ImageId,'Nat-10v') ...
	  | strcmp(ImageId,'Nat-10m')
      x=[];
      CompSort=reshape(sort(Components),5,2)';
      for k=1:5
	x=[x [load_texture(CompSort(1,k),128,128,'HistEq','ZeroMean'); ...
	    load_texture(CompSort(2,k),128,128,'HistEq','ZeroMean')]]; 
      end
    else
      [Ref,N_tex]=cp_image(ImageId,0);
      [M,N]=size(Ref);
      x=zeros(M,N);
      for k=min(Ref(:)):max(Ref(:))
	x=x + (Ref==k).*cp_image(Components(k),[],M,N);
      end
    end
  else
    if TrainNo==0 % 			Load reference image
      if strcmp(ImageId,'Nat-10') ...
	    | strcmp(ImageId,'Nat-10v') ...
	    | strcmp(ImageId,'Nat-10m')
	x=[2*ones(128) 5*ones(128) 8*ones(128) 1*ones(128)  4*ones(128) ; ...
	  7*ones(128) 9*ones(128) 3*ones(128) 6*ones(128) 10*ones(128) ];
      else
	FileName=sprintf('%s/datafiler/textures/%s.ref.hips',...
	  getenv('HOME'),ImageId);
	x=hipsread(FileName);
      end
    end
  end
else 
  % ---------------------------------------------------------------------
  % Load an image identified by its components
  % ---------------------------------------------------------------------
  Components=ImageId;
  N_tex=length(Components);
  if isempty(TrainNo) % 		Image component
    for k=1:N_tex
      x=[x load_texture(Components(k),M,N,'HistEq','ZeroMean')]; 
    end
  elseif TrainNo==0 % 		Reference image
    for k=1:N_tex, 
      x=[x k*ones(M,N)]; 
    end
  end
end

% ---------------------------------------------------------------------
% Load training subimages, if requested (must not intersect test region)
% ---------------------------------------------------------------------
if TrainNo>0
  % Compute offsets
  if sum(Components(TrainNo)==[4 29 84]),                    M_o=306; N_o=306;
  elseif sum(Components(TrainNo)==[6 21]),                   M_o=350; N_o=350;
  elseif sum(Components(TrainNo)==[5 9]),                    M_o=350; N_o= 50;
  elseif sum(Components(TrainNo)==[24 38]),                  M_o= 50; N_o=350;
  elseif sum(Components(TrainNo)==[36 37]),                  M_o=256; N_o= 50;
  elseif sum(Components(TrainNo)==[3 32 33 54 55 57 68 77]), M_o= 50; N_o= 50;
  elseif sum(Components(TrainNo)==[1001 1002 1003 1004]),    M_o=0; N_o=0;
  elseif sum(Components(TrainNo)==[53]),                     M_o= 50; N_o=306;
  elseif sum(Components(TrainNo)==[28]),                     M_o= 50; N_o=256;
  else,                                                      M_o=156; N_o=256;
  end
  [x,r_xx]=load_texture(Components(TrainNo),M,N,'Offset',M_o,N_o,...
    'HistEq','ZeroMean');
end
