showit       = 1;
k_const      = 2^(1/3) ; # usually 2 or sqrt(2)
%raw_filename = "c[v028.jpg-v029.jpg].u32.rgb.raw";
lp=LOADPATH;
LOADPATH=['/home/mann/comparasolve-1.3:' lp ':'];

global C a c s k;

%showit       = 1;
%k_const      = (2); % usually 2 or sqrt(2)
%raw_filename = "comparagram.u32.rgb.raw";

global C;
%fin =fopen(raw_filename,"r");
%fout=fopen([raw_filename ".out.estacs"],"w");
chan=0;
%while (feof(fin)==0)
  % read 256x256 block from raw file
%  [C count]=fread(fin, 256*256, "unsigned int", 0);
%  [C count]=fread(fin, 256*256, "unsigned int", 0); % GET GREEN ONLY
%  if (count<256*256) break; endif;
%  chan++;
%  C=reshape(C, 256, 256);
%  C=C;%'; %long story... im doing this all wrong i know

  load /home/mann/comparunroll/Jsum03.mat
  C = outputmat;
  % store the original data, compute max ridge D
  Csave=C;
  D=C;
  mc=max(C);for i=1:size(D,1); D(i,:)=D(i,:)-mc; end;
  [I J]=find(D==0);
  D=zeros(256,256);
  for i=1:size(I,1); D(I(i),J(i))=C(I(i),J(i)); end;

  % now normalize weights by row/column
  Csave=kruithoff(Csave);
  D=kruithoff(D)*10000;
%  %D(2     :256-64,256-64:256-5) = D(2     :256-64,256-64:256-5) * 64;  % top edge
%  D(256-64:256-5 ,2     :253  ) = D(256-64:256-5 ,2     :253  ) * 64;  % right edge
%  D(256-64:256-5 ,256-64:256-5) = D(256-64:256-5 ,256-64:256-5) / 64;  % dont do twice
%  %D(2:253,2:128) = D(2:253,2:128) * 16;       % bottom swath
%  D(6:250,6:250) = D(6:250,6:250) * 32;     % entire central region
%  D(250:256,:) = D(250:256,:)./16;             % rightmost column
%  %D(195,256)=max(max(D))*200000;

  % set starting point (a=1 or -1, c=1, k=k_const, s=1)

  % subsequent channels start from last channel's parameter values
  % should only do so if last pass converged
  % also if some didn't and some did, redo non-converged cases 
  % starting from convergent ones
%  if (chan==1)
%    [a c]=trial_ac; 
%    a=sqrt(e^(a)); 
%    c=sqrt(c);  
%    k=sqrt(k_const); 
%    s=sqrt(2);
%  end%if
  
%  if(showit); bfplotmat(C, [0], log(a*a), (c)^2,(k)^2, s*s); endif;
%  x=[a c s]';

  % solve for ridge parameters
  C=D;
%  x=bfgs_hess('comparasolve_f_k', 'comparasolve_d2f_k', x);
%  a=x(1); c=x(2);  s=x(3); %k=x(3); s=x(4);
%  if(showit); bfplotmat(C, [0], log(a*a), (c)^2,(k)^2, s*s); endif;
  imshow(C);

  % solve original data parameters
  C=Csave+D; %PLUS RIDGE!!
  %x=bfgs_hess('comparasolve_f_k', 'comparasolve_d2f_k', x);
  %a=x(1); c=x(2);  s=x(3); %k=x(3); s=x(4);
  %if (showit); bfplotmat(Csave, [0], log(a*a), (c)^2,(k)^2, s*s); endif;
  imshow(C);

%endwhile
%fclose(fin);

