function FeatureSubset = GetC1FeatureSubset(WhichBands, c1params, originalimagesize,options)
%function FeatureSubset = GetC1FeatureSubset(WhichBands, c1params, originalimagesize,options)
%
%The c1 feature vector contains elements from many bands.  This function returns the indicies of the features
%which lie within the indicated bands.
%WhichBands is a row vector containing a list of those bands to be included e.g. [1 5 4]

d.rot = [90 -45 0 45];
d.c1ScaleSS = [1:2:18];
d.RF_siz    = [7:2:39];
d.c1SpaceSS = [8:2:22];
d.div = [4:-.05:3.2];
d.includeborder = 1;
d.includelaplacianfilter = 0;
d.includedoubleorientations = 0;
d.includegrayscalefilter = 0;
d.smallerscalesonly = 0;
c1params = ResolveMissingOptions(c1params,d);

if(nargin < 4)
  options = [];
end
Defoptions.includeOrientations = ones(1,length(c1params.rot));
Defoptions.includeLaplacian = 0;
Defoptions.includeGS = 0;
options = ResolveMissingOptions(options, Defoptions);
[fSiz,filters,c1OL,numSimpleFilters] = init_gabor_alpha(c1params.rot, c1params.RF_siz, c1params.div,c1params);

z = zeros(originalimagesize);
[c1] = C1(z, filters, fSiz, c1params.c1SpaceSS, c1params.c1ScaleSS, c1OL,0);
FeatureSubset = [];
includeFilters = [];
includeFilters = options.includeOrientations;
if(c1params.includelaplacianfilter)
   includeFilters = [includeFilters,options.includeLaplacian];
end
if(c1params.includegrayscalefilter)
   includeFilters = [includeFilters,options.includeGS];
end
for i = 1:length(c1)
  if (intersect(WhichBands,i))
     c1{i} = 0*c1{i}+1;
     for j = 1:length(includeFilters)
       c1{i}(:,:,j) = includeFilters(j);
     end
     FeatureSubset = [FeatureSubset; c1{i}(:)];
  else
     FeatureSubset = [FeatureSubset; 0*c1{i}(:)];
  end
end
FeatureSubset = find(FeatureSubset);
