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;

% correlation feature templates
M = [450 500 550 600 650 700 750 800 850 900 950];
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));
   tmin(i) = min(L);
   T(:,i) = feature_template(L, M(i), W(i))';
end

% 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 = 'S1  S2 S3 S4 S5 S6 S7 S8 S9 S10 S11  ';
f7 = 'Red  Green  Intensity  ';

% image features
f8 = 'CameraRed  CameraGreen  CameraIntensity';

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

% create pdf (assumes all output values normalized to be in [-1,1] or [0,1] range
x = -1:0.05:1;
dx = [diff(x) 0];
pdf = zeros(17, length(x));

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(tmin)
            itmin(i) = find(lambda == tmin(i));
         end
         
         figure(1);
         PlotSpectra(l,spectrum(imin:imax),S.specfile);
         
         % compute spectral properties
         % extract features features S1, S2, ....
            SF = SpectralCorrelation(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];
            
            % update pdf's
            for i = 1:length(SF)
               ParzenWindow(i,:) = gaussian(x,SF(i),0.1);
            end
            ParzenWindow(12,:) = gaussian(x,r,0.1);
            ParzenWindow(13,:) = gaussian(x,g,0.1);
            ParzenWindow(14,:) = gaussian(x,Intensity,0.1);
            for i = 1:size(ParzenWindow,1)
               I = sum(dx.*ParzenWindow(i,:));
               pdf(i,:) = pdf(i,:) + ParzenWindow(i,:)/I;
            end
                        
         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);
      
      % update pdf's
      ParzenWindow(15,:) = gaussian(x,cr,0.1);
      ParzenWindow(16,:) = gaussian(x,cg,0.1);
      ParzenWindow(17,:) = gaussian(x,Cint,0.1);
      for i = 15:17
         I = sum(dx.*ParzenWindow(i,:));
         pdf(i,:) = pdf(i,:) + ParzenWindow(i,:)/I;
      end

   else
      fprintf(casefile,'*  *  *');
   end
         
   drawnow
   fprintf(casefile,'\n');
   L = readline(fid);
end

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

% print histograms
I = sum(pdf')';
pdf = pdf ./ (I * ones(1,size(pdf,2)));
figure(3)
plot(x, pdf);