disp('MATLAB Script File:  cofm1.m') ;
%
%   For the Cornea Study, this does manipulations with groups of
%   feature matrices (copied from cocd4.m)
%
%   Currently tuned for finding specific outliers.
%
format compact ;
format short ;


if0l1 = 1 * 10 + 1 ;


igroup = 9 ;  %  1 - Normal       (Central Image = "All K")
              %  2 - Ker Suspect  (Central Image = "All K")
              %  3 - Kerataconic  (Central Image = "All K")
              %  4 - Pre          (Central Image = "All P")
              %  5 - 1 Mo.        (Central Image = "All P")
              %  6 - 1 Yr.        (Central Image = "All P")
              %  7 - Pre - 1 Mo.  (Central Image not inter.)
              %  8 - Pre - 1 Yr.  (Central Image not inter.)
              %  9 - 1 Mo - 1 Yr  (Central Image not inter.)

ictype = 0 ;   %  0 - no centering
               %  1 - center about "central image" for this group
               %  2 - center about mean for this group



if igroup == 1 ;
  namestr = '\matlab\steve\cornea\fm\fnnor.asc' ;
  imagstr = 'Norm' ;
  icent = 119 ;    %  cofmap1.m index for Centering Feature Matrix
              %  119 - Mean All K
  centstr = '- all K' ;
  nincol = 1 ;
elseif igroup == 2 ;
  namestr = '\matlab\steve\cornea\fm\fnksus.asc' ;
  imagstr = 'Ksus' ;
  icent = 119 ;    %  cofmap1.m index for Centering Feature Matrix
              %  119 - Mean All K
  centstr = '- all K' ;
  nincol = 1 ;
elseif igroup == 3 ;
  namestr = '\matlab\steve\cornea\fm\fnker.asc' ;
  imagstr = 'Ker' ;
  icent = 119 ;    %  cofmap1.m index for Centering Feature Matrix
              %  119 - Mean All K
  centstr = '- all K' ;
  nincol = 1 ;
elseif igroup == 4 ;
  namestr = '\matlab\steve\cornea\fm\fnpre.asc' ;
  imagstr = 'Pre' ;
  icent = 122 ;    %  cofmap1.m index for Centering Feature Matrix
              %  122 - Mean all P
  centstr = '- all P' ;
  nincol = 1 ;
elseif igroup == 5 ;
  namestr = '\matlab\steve\cornea\fm\fn1mo.asc' ;
  imagstr = '1 Mo.' ;
  icent = 122 ;    %  cofmap1.m index for Centering Feature Matrix
              %  122 - Mean all P
  centstr = '- all P' ;
  nincol = 1 ;
elseif igroup == 6 ;
  namestr = '\matlab\steve\cornea\fm\fn1yr.asc' ;
  imagstr = '1 Yr.' ;
  icent = 122 ;    %  cofmap1.m index for Centering Feature Matrix
              %  122 - Mean all P
  centstr = '- all P' ;
  nincol = 1 ;
elseif igroup == 7 ;
  namestr = '\matlab\steve\cornea\fm\fnpre1mo.asc' ;
  imagstr = 'Pre - 1Mo' ;
  icent = 0 ;   %  no need to center differences
  centstr = '' ;
  nincol = 2 ;
elseif igroup == 8 ;
  namestr = '\matlab\steve\cornea\fm\fnpre1yr.asc' ;
  imagstr = 'Pre - 1Yr' ;
  icent = 0 ;   %  no need to center differences
  centstr = '' ;
  nincol = 2 ;
elseif igroup == 9 ;
  namestr = '\matlab\steve\cornea\fm\fn1mo1yr.asc' ;
  imagstr = '1Mo - 1Yr' ;
  icent = 0 ;   %  no need to center differences
  centstr = '' ;
  nincol = 2 ;
end ;


nr = 10 ;
nc = 7 ;



disp('Loading feature matrices for this group') ;
fid = fopen(namestr,'r') ;
if nincol == 1 ;
  mfindx = fscanf(fid,'%4d') 
elseif nincol == 2 ;
  mfindx = fscanf(fid,'%4d %4d') ;
  mfindx = reshape(mfindx,2,length(mfindx)/2)' ;
end ;

%  Used this for checking
%mfindx
%pause ;


prestr = '\matlab\steve\cornea\fm\fm42' ;
mfa = [] ;
for i = 1:size(mfindx,1) ;

  if mfindx(i,1) < 10 ;
    eval(['load ' prestr '000' num2str(mfindx(i,1))]) ;
  elseif mfindx(i,1) < 100 ;
    eval(['load ' prestr '00' num2str(mfindx(i,1))]) ;
  elseif mfindx(i,1) < 1000 ;
    eval(['load ' prestr '0' num2str(mfindx(i,1))]) ;
  else ;
    eval(['load ' prestr num2str(mfindx(i,1))]) ;
  end ;
          %  This creates a number of variables, especially fm,
          %  a 19 x 7 "feature matrix"

  if nincol == 2 ;    %  Then need to subtract second feature matrix
    fm1 = fm ;
    if mfindx(i,2) < 10 ;
      eval(['load ' prestr '000' num2str(mfindx(i,2))]) ;
    elseif mfindx(i,2) < 100 ;
      eval(['load ' prestr '00' num2str(mfindx(i,2))]) ;
    elseif mfindx(i,2) < 1000 ;
      eval(['load ' prestr '0' num2str(mfindx(i,2))]) ;
    else ;
      eval(['load ' prestr num2str(mfindx(i,2))]) ;
    end ;
          %  This creates a number of variables, especially fm,
          %  a 19 x 7 "feature matrix"
    fm = fm1 - fm ;
  end ;

  %  Convert c-s to 10 x 7 amplitude - phase form, 
  %        using block from cofmtst1.m
  mamp = abs(fm(1,:)) ;
  mpha = atan2(0,fm(1,:)) * 180 / pi ;
  for i = 1:9 ;
      vamp = sqrt(fm(2*i,:).^2 + fm(2*i+1,:).^2) ;
    mamp = [mamp; vamp] ;
      vpha = atan2(fm(2*i+1,:),fm(2*i,:)) * 180 / pi ;
        vind = find(vpha < 0) ;
      vpha(vind) = vpha(vind) + 360 ;
    mpha = [mpha; vpha] ;
  end ;

  %  Feature matrix form for output
  fmamp = fm(1,:) ;   %  To get signs right
  fmamp = [fmamp; mamp(2:10,:)] ;  %  Use amplitude for rest
  mfa = [mfa fmamp(:)] ;
          %  put matrix into a column vector and add to right
end ;



if  ictype > 0  &  nincol == 1  ;
  disp('Do Centering') ;
  if ictype == 1 ;   %  Then do "central image" centering
    if icent == 119 ;
    %  imagstr = 'Mean All K Study' ;
              %  All 3 groups in Keratoconus study
              %          - called kstud by NKT
              %          - called ker by JTZ
      fmistr = 'mn42ker.mat' ;
    elseif icent == 122 ;
    %  imagstr = 'Mean All P Study' ;
              %  All 4 groups in PRK study
              %          - called pstud by NKT
              %          - called pos by JTZ
              %                   nor by JTZ are not used (same as pre).
      fmistr = 'mn42pos.mat' ;
    end ;
    ['load \matlab\steve\cornea\fm\' fmistr]
    eval(['load \matlab\steve\cornea\fm\' fmistr]) ;
    mamp = mnamp ;
    mpha = mnpha ;
    %  Feature matrix form for output
    fmamp = mamp(1,:) .* cos(mpha(1,:) * pi / 180) ;
              %  Here Jinting used names mnamp & mnpha
    fmamp = [fmamp; mamp(2:10,:)] ;  %  Use amplitude for rest
    
    mfa = mfa - vec2mat(fmamp(:),size(mfindx,1)) ;
  elseif ictype == 2 ;  %  Then Center at Mean

    mfa = mfa - vec2mat(mean(mfa')',size(mfindx,1)) ;

    centstr = '- mean' ;
  end ;
else ;  %  No centering is done
  centstr = '' ;
end ;



d = size(mfa,1) ;
n = size(mfa,2) ;

disp(['joint feature matrix has size ' num2str(d) ' rows ']) ;
disp(['and ' num2str(n) ' cols']) ;


fclose('all') ;

disp('F0 L1 coeffs are:') ;
vcf0l1 = mfa(if0l1,:)' 


disp(['For ' imagstr ':']) ;
disp(['number of cases is: ' num2str(size(mfa,2))])
if igroup == 7 ;
  [cmax imax] = max(vcf0l1) ;
  disp(['max pre file number is: ' num2str(mfindx(imax,1))]) ;
  disp(['max 1mo file number is: ' num2str(mfindx(imax,2))]) ;
  disp(['    max coeff is: ' num2str(cmax)]) ;
  disp(['    max index is: ' num2str(imax)]) ;
elseif igroup == 9 ;
  [cmin imin] = min(vcf0l1) ;
  disp(['min 1mo file number is: ' num2str(mfindx(imin,1))]) ;
  disp(['min 1yr file number is: ' num2str(mfindx(imin,2))]) ;
  disp(['    min coeff is: ' num2str(cmin)]) ;
  disp(['    min index is: ' num2str(imin)]) ;
end ;






