function [R,options] = Im2RandomImageFeature(Im,options);
%A random image feature is designed based on the RandSeed option, and applied
%to the input image 'Im'  This feature consists of a concatination of filters
%morphological operations, and resizings.



if(nargin < 2), options = [];, end
d.RandSeed = 1;
d.MaxFeatureSize = 1000;
d.RandFeatureTypes = [];
d.caParams = [];

options = ResolveMissingOptions(options,d);

seedrand(options.RandSeed);
if(isempty(options.RandFeatureTypes))
  RandFeatureTypes = ceil(rand(1000,1)*5);
else
  RandFeatureTypes = options.RandFeatureTypes;
end
if(not(iscell(options.caParams)))
  [caParams,RandFeatureTypes] = BuildCAParams(RandFeatureTypes,[size(Im,1),size(Im,2),4],options);
else
  caParams = options.caParams;
end
options.RandFeatureTypes = RandFeatureTypes;
options.caParams = caParams;
Im = FourWayFilter(Im);
R = ExecuteFeature(Im, RandFeatureTypes, caParams);
R = R(:);

function R = ExecuteFeature(Im,RandFeatureTypes, caParams)
for iExec = 1:length(RandFeatureTypes)
  clear R
  switch(RandFeatureTypes(iExec))
    case 1 % maxfilter
      for iOri = 1:4
        R(:,:,iOri) = maxfilter_arbitrarySupport(Im(:,:,iOri), caParams{iExec}.filter);
      end
    case 2 % minfilter
      for iOri = 1:4
        R(:,:,iOri) = maxfilter_arbitrarySupport(Im(:,:,iOri), caParams{iExec}.filter);
      end
    case 3 % regufilter
      for iOri = 1:4
        R(:,:,iOri) = abs(imfilter(Im(:,:,iOri), caParams{iExec}.filter));
      end
    case 4 % generic resize
      for iOri = 1:4
        R(:,:,iOri) = imresize(Im(:,:,iOri), caParams{iExec}.newsize,caParams{iExec}.method);
      end
    case 5 % decmiate
      if(caParams{iExec}.dim == 1)
         R = Im(1:caParams{iExec}.span:end,:,:);
      end
      if(caParams{iExec}.dim == 2)
         R = Im(:, 1:caParams{iExec}.span:end,:);
      end
  end
  Im = R;
end    
options.RandFeatureTypes = RandFeatureTypes;
options.caParams = caParams;



function [caParams,RandFeatureTypes] = BuildCAParams(RandFeatureTypes,sz,options);
caParams = {};
n = 0;
seedrand(options.RandSeed);
newsz = sz;
while(prod(newsz) > options.MaxFeatureSize)
  n = n+1;
  switch(RandFeatureTypes(n))
    case 1 % maxfilter
      caParams{n}.filter = PickRandomMorphFilter(rand(2,1));
    case 2 % minfilt
      caParams{n}.filter = PickRandomMorphFilter(rand(2,1));
    case 3 % filt
      caParams{n}.filter = PickRandomLinFilter(rand(2,1));
    case 4 % generic resize
      caParams{n}.newsize(1:2) = ceil((.25 + rand * .6) * newsz(1:2));
      r = rand > .5;
      if(r)
         caParams{n}.method = 'bilinear';
      else
         caParams{n}.method = 'nearest';
      end
      newsz = [caParams{n}.newsize,4];
    case 5 % decimate
      caParams{n}.dim = ceil(rand * 2);       
      caParams{n}.span = ceil(rand * 8);
      l = length(1:caParams{n}.span:newsz(caParams{n}.dim));
      newsz(caParams{n}.dim) = l;
  end
end
RandFeatureTypes = RandFeatureTypes(1:n);

function f = PickRandomLinFilter(r)
% the first rand val determines the type, deriv or gaus
% the second determines orientation or sigma
if(r(1) < .5)  % gausfilter
   f = GausImg([7 7],[4 4],ceil(r(2) * 10));
else
   f = Gabor(7,3,r(2)*360,1,.25);
end



function f = PickRandomMorphFilter(r)
% morph filters are binary 5 x 5 matricies which are a band in one orientation
f = zeros(5);
o = r(1) * pi;  % first random var deterimines orientation
Th = r(2) * 7;  % second random var deterimines thickness
u = [sin(-o),cos(-o)];
u = u / norm(u);
for x = 1:7 
  for y = 1:7
    x2 = x - 7 / 2 - .5;  
    y2 = y - 7 / 2 - .5;  
    r = [y2, x2];
    d = abs(r * u');
    if(d < Th)
      f(y,x) = 1;
    end
  end
end 

function I = maxfilter_arbitrarySupport(I,Msupport);
%function I = maxfilter_arbitrarySupport(I,Msupport);
%
%support must be a binary matrix.  The centroid of the support will be at
%[floor(size(support,1)/2), floor(size(support,2)/2)]

sz = size(Msupport);
centroid = floor(sz);
fsupport = find(Msupport);
ndxsupport = NDind2sub(sz,fsupport);
nPoints = size(fsupport,1);
ndxsupport = ndxsupport - repmat(centroid, [nPoints,1]);
padamnt = max(ceil(sz/2));
ndxsupport = ndxsupport + padamnt;
ndxsupport = ndxsupport + padamnt;
padI = padimage(I,padamnt,-inf);
BigMat = zeros(nPoints, prod(size(I)));
sI = size(I);
for i = 1:nPoints
  crop = padI( ndxsupport(i,1) +  (1:sI(1))-1,ndxsupport(i,2) + (1:sI(2)) - 1 );
  BigMat(i,:) = crop(:);
end
I = reshape(max(BigMat,[],1),sI);
