function scan_feature(FeatureFN,N_spv,N_sph,FeatureFN2)
% SCAN_FEATURE  Scans through all feature frames, and prints the result
%	scan_feature(FeatureFN,N_spv,N_sph)
%	scan_feature(FeatureFN,N_spv,N_sph,FeatureFN2)
%	
%	scan_feature('features.hips',2,2,'trfeatures.hips');
%	scan_feature('trfeatures.hips',4,1,'features.hips');

if exist(sprintf('%s_000',FeatureFN))
  N_features=1;
  while exist(sprintf('%s_%03d',FeatureFN,N_features))
    N_features=N_features+1;
  end
else
  N_features=hsize(FeatureFN,'Frames');
end
N_sp=N_spv*N_sph;

for n=0:N_features-1
  if exist(sprintf('%s_%03d',FeatureFN,n))
    v=hipsread(sprintf('%s_%03d',FeatureFN,n));
  else
    v=hipsread(FeatureFN,[],[],n);
  end
  if (exist('FeatureFN2')==1) 
    if exist(sprintf('%s_%03d',FeatureFN2,n))
      v2=hipsread(sprintf('%s_%03d',FeatureFN2,n));
    else 
      v2=hipsread(FeatureFN2,[],[],n);
    end
  else
    v2=[min(v(:)) max(v(:))];
  end
  MinMax=[min([v(:); v2(:)]) max([v(:); v2(:)])];
  [M,N]=size(v);
  MN=min([M N]);
  if MN>64
    N_dec=round(MN/64);
    v=v(1:N_dec:M,1:N_dec:N);
  end
  if rem(n,N_sp)==0, clf; end
  subplot(N_spv,N_sph,rem(n,N_sp)+1);
  v(1,1:2)=MinMax; % 			Ensure proper scaling
  gimage(v);
  title(sprintf('Feature %d',n));
  if (rem(n,N_sp)==N_sp-1) | (n==N_features-1)
    print -dps2 scan_feature.ps
    !lpr scan_feature.ps
    !lpq
  end
end