function [SO, IO] = process_filesIV(index_file, casefilename)

% [SpecFeatures, ImFeatures] = process_filesII(index_filename, case_filename)
% returns matrices of observations (each row a new observation)

%initialize variables etc.

% spectral range of interest.
lambda_min = 400;
lambda_max = 2500;
l = lambda_min:lambda_max;

% VISIBLE correlation feature templates V1 - V11
M = [450 500 550 600 650 700 750 800 850 900 950];
m = M;
W = [50  50  50  50  50  50  50  50  50  50  50];
for i = 1:length(M)
   L = (M(i)-W(i)):(M(i)+W(i));
   Vtmin(i) = min(L);
   VT{i} = feature_template(L, M(i), W(i))';
end

% NEAR IR correlation feature templates NIR1 - NIR5
M = [1050 1100 1150 1200 1250];
m = [m M];
W = [50   50   50   50   50];
for i = 1:length(M)
   L = (M(i)-W(i)):(M(i)+W(i));
   NIRtmin(i) = min(L);
   NIRT{i} = feature_template(L, M(i), W(i))';
end

% IR correlation feature templats MIR1 - MIR9
M = [1500 1550 1600 1650 1700 1750 2100 2200 2300];
m = [m M];
W = [50   50   50   50   50   50   100  100  100];
for i = 1:length(M)
   L = (M(i)-W(i)):(M(i)+W(i));
   IRtmin(i) = min(L);
   IRT{i} = feature_template(L, M(i), W(i))';
end

% ice/snow correlation feature templates ICE
M = [765 950 1135 1250];
m = [m M];
W = [15  30  35   10];
for i = 1:length(M)
   L = (M(i)-W(i)):(M(i)+W(i));
   ICEtmin(i) = min(L);
   ICET{i} = feature_template(L, M(i), W(i))';
end

T = [VT NIRT IRT ICET];  % templates, m is the template centers
ltmin = [Vtmin NIRtmin IRtmin ICEtmin];

% spectral filters
load rgb.mat
MaxIntensity = sum(sum([R G B]));
% R, G, B, l_rgb

% start reading and processing files;
fid = OpenDataFile(index_file, 3);
L = readline(fid);
Err = 0;
i = 1;

% open CASE file for output.
f1 = 'IDnum  ';
f2 = 'sample_id  ';
f3 = 'basic_class  ';
f4 = 'classification ';
f5 = 'formation  ';

% spectral features
f6 = 'V1 V2 V3 V4 V5 V6 V7 V8 V9 V10 V11 ';
f7 = 'NIR1 NIR2 NIR3 NIR4 NIR5 ';
f8 = 'IR1 IR2 IR3 IR4 IR5 IR6 IR7 IR8 IR9 '
f9 = 'ICE1 ICE2 ICE3 ICE4  ';
f10 = 'Red  Green  Intensity  ';

% image features
f11 = 'CameraRed  CameraGreen  CameraIntensity';

casefile = OpenCaseFile(casefilename, [f1 f2 f3 f4 f5 f6 f7 f8 f9 f10 f11]);

while (L > 0)
   S = ParseDataFileLine(L)
   
   % write first fields of S
   fprintf(casefile,'%d\t', S.idnum);
   fprintf(casefile,'%s\t', S.sample_id,S.basic_class,S.classification, S.formation);
   
   % read spectrum
   if (S.specfile(1) ~= '*')
      [spectrum,lambda] = ReadSpectrum(S.specfile);
      
   
   	if (spectrum ~= -1)  
      	% spectrum read successfully!
      
      	imin = find(lambda == lambda_min);
         imax = find(lambda == lambda_max);
         for i = 1:length(ltmin)
            itmin(i) = find(lambda == ltmin(i));
         end
         
         figure(1);
         PlotSpectra(l,spectrum(imin:imax),S.specfile);
         
         % compute spectral properties
         % extract features features S1, S2, ....
            SF = SpectralCorrelationII(spectrum, T, itmin);
            hold on
			   plot(m,SF,'k*');

            % extract R, G and Intensity
            fmin = find(lambda == l_rgb(1));
            F = FilterCoeffs(spectrum, [R G B], fmin);
            Intensity = sum(F);
            
            r = F(1)/Intensity;
            g = F(2)/Intensity;
            b = F(3)/Intensity;
            Intensity = Intensity/MaxIntensity;
            
            bar([550 600 650],[0 0 r],'r');
 			   bar([550 600 650],[0 g 0],'g');
 			   bar([550 600 650],[b 0 0],'b');
            hold off
            
            % write to case file
            fprintf(casefile,'%g  ',SF, r, g, Intensity);  
            % update output matrix
            SO = [SO; SF r g Intensity];
                                              
         end
      else
         fprintf(casefile,'*  *  *  *  *  *  *  *  *  *  *  *  *  *  ');
      end  % spectrum processing
   
   
   % read image
   if (S.segimagefile(1) ~= '*')
      segimage = ReadImage(S.segimagefile);
%      figure(2)
%      image(segimage);
      
      % compute RGB parameters & write to case file
      [cr, cg, cb] = SegImageRGB(segimage,[255 255 255]);
      Cintensity = cr + cg + cb;
      cr = cr / Cintensity
      cb = cb / Cintensity
      cg = cg / Cintensity;
      Cint = Cintensity/765
      
      % update casefile
      fprintf(casefile,'%g  ', cr, cg,  Cint);
      
   else
      fprintf(casefile,'*  *  *');
   end
         
   drawnow
   fprintf(casefile,'\n');
   L = readline(fid);
end

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