function SO = ProcessAsterFiles(index_file, casefilename)

% SO = ProcessAsterFiles(index_file, casefilename)
%
% Process data in index_file to generate a Netica casefile.  Index file is a text file,
% each column representing a field, each row a new observation.  This is exclusively for 
% the Aster spectral library.

% Liam Pedersen, NASA ARC/CMU 1999

%**************************
%initialize variables etc.
%**************************

%************************************************************
% spectral feature extraction global data

% spectral range of interest.
lambda_min = 400;      
lambda_max = 2500;      % ASDI FieldSpec FR spectrometer
%lambda_max = 1000;       % OceanOptics S2000 spectrometer

% create spectral templates T and define wavelength range.
l = lambda_min:lambda_max;
%[T, ltmin, m, SpecLabels] = OOIS2000_SpecFeatureTemplates(1);
[T, ltmin, m, SpecLabels]  = ASDIFR_SpecFeatureTemplates(1);

% spectral filters
load rgb.mat % loads vectors R, G, B and l_rgb.  RGB filter fns + wavelength range.
SpectralFeatureLabels = [SpecLabels 'Red Green Intensity '];

% Number of spectral features
NSpecFeatures = length(ltmin) + 4;

% casefile header labels
label1 = 'IDnum  Name  Type  Class  SubClass  ';
CaseFileLabels = [label1 SpectralFeatureLabels];

%*****************************************************************
%*****************************************************************

% open CASE file for output.
casefile = OpenCaseFile(casefilename, CaseFileLabels);

% Open INDEX file and process each line
Err = 0; i = 1;
fid = OpenDataFile(index_file, 3);
L = readline(fid);

while (L > 0)
   S = ParseDataFileLineII(L);
   
   % read in Aster spectra file
   
   % write first fields of S
   fprintf(casefile,'%d\t', S.idnum);
   
   %********************************
   % Process SPECTRUM (if possible).
   
   if (S.specfile(1) ~= '*')
      % load spectra and extract feature vector: (local function)
      [spectrum,lambda,Name,Type,Class,Subclass] = ReadAsterSpectra(S.specfile);
      if (spectrum ~= -1)  % spectrum read successfully!
         SFV = SpectralFeatureVector(spectrum,lambda,T,ltmin, m, R,G,B,l_rgb,lambda_min,lambda_max);
         disp(Name);
         % write to case file
         fprintf(casefile,'"%s"\t', Name);
         fprintf(casefile,'%s\t', Type, Class, Subclass);

         fprintf(casefile,'%g  ',SFV);  
      
         % update output matrix
         SO = [SO; S.idnum SFV];
      else
         disp('Unable to read spectral file');
         S.specfile
         for i = 1:(NSpecFeatures+3)
            fprintf(casefile,'*  ');
         end
      end
   else
      for i = 1:(NSpecFeatures+3)
         fprintf(casefile,'*  ');
      end
   end  % spectrum processing
      
   drawnow
   fprintf(casefile,'\n');
   L = readline(fid);
end

% close I/O files
fclose(fid);
fclose(casefile);


%**********************************************************************
%**********************************************************************


%*************************************************************
% Process SPECTRUM

function SFV = SpectralFeatureVector(spectrum, lambda, Templates, LT, m, R, G, B, l_rgb, lmin, lmax)

% SFV = SpectralFeatureVector(specfile, Templates, LT, m, R, G, B, l_rgb)
%		spectrum = spectra vector
%     lambda   = wavelength values for each entry in spectrum (1nm spacing assumed)
%		Templates = cell array of feature templates (1 nm spacing)
%		LT = array of starting wavelengths (nm) of feature templates
%		R, G, B, l_rgb are R,G,B filter fns, along with wavelength vector l_rgb (1nm spacing)
%		lmin, lmax = wavelength range (nm).

% Liam Pedersen, NASA ARC/CMU, August 1998

% set for plotting of spectra and spectral features.
plotspecfeatures = 1;

% get indices of wavelength limits in spectrum vector
if (lmin < lambda(1)) lmin = lambda(1); end
if (lmax > max(lambda)) lmax = max(lambda); end
imin = find(lambda == lmin);
imax = find(lambda == lmax);

for i = 1:length(LT)
   itmin(i) = find(lambda == LT(i));
end

% compute spectral template correlations:
SF = SpectralCorrelationII(spectrum, Templates, itmin);

% compute RGI coefficients
dl = lambda(1) - l_rgb(1);
if (dl <= 0)
   fmin = find(lambda == l_rgb(1));
   MaxIntensity = sum(sum([R G B]));
   F = FilterCoeffs(spectrum, [R G B], fmin);
else
   MaxIntensity = sum(sum([R(dl:length(R)) G(dl:length(G)) B(dl:length(B))]));
   F = FilterCoeffs(spectrum, [R(dl:length(R)) G(dl:length(G)) B(dl:length(B))], 1);
end

Intensity = sum(F);
r = F(1)/Intensity;
g = F(2)/Intensity;
b = F(3)/Intensity;
Intensity = Intensity/MaxIntensity;
   
SFV = [SF r g Intensity];

if (plotspecfeatures)
   figure(1);
   clf;
   PlotSpectralFeatures(lmin:lmax, spectrum(imin:imax), [SF r g b], m, '');
   do = 'PlotSpectralFeatures(lmin:lmax, spectrum(imin:imax), [SF r g b], m, '''');';
%   eval(do,'disp(''failed to draw'')');
   drawnow
end
   

