addpath('C:\Documents and Settings\Steve\My Documents\Research\matlab')
addpath('C:\Documents and Settings\Steve\My Documents\Research\optNew\linlsqUtil')

% petFileName = 'C:\Documents and Settings\Steve\My Documents\Research\Human\HeadAndNeck\HN003\MatlabData\cuATSMoptFlipped.am';

% PET = am2geom(petFileName);
TDP = open_presc('C:\Documents and Settings\Steve\My Documents\Research\optNew\linlsqOpt\input\HN003laplace\continuous\PTV70dosePlus.bin');
% TDPfull = full3D(TDP);
% TDPbinary = single(~~TDPfull);
% PETcorrected = PET.data.*TDPbinary;
% PETsparse = sparse3D(PETcorrected);
% PETnorm = PETsparse;
% PETnorm.non_zero_values = PETsparse.non_zero_values./(max(PETsparse.non_zero_values));

presc = TDP;

ind = {};

for k=1:10;
    
    dpresc = (max(presc.non_zero_values)-min(presc.non_zero_values))./k;
                    
    minpresc = min(presc.non_zero_values);

    for n=1:Nstep;
        ind{n} = find(presc.non_zero_values >= (n-1)*dpresc + minpresc ...
            & presc.non_zero_values < n*dpresc + minpresc);
        presc.non_zero_values(ind{n}) = mean(presc.non_zero_values(ind{n}));
    end
    
%     for i=1:k;
%         ind{i} = find(PETnorm.non_zero_values >= (i-1)./k & PETnorm.non_zero_values < i./k);
%         PETnorm.non_zero_values(ind{i}) = double(20*i./k);
%     end

%     PETpresc = PETsparse;
%     PETpresc.non_zero_values = 70.+PETnorm.non_zero_values;
%     PETnorm.non_zero_values = PETsparse.non_zero_values./(max(PETsparse.non_zero_values));
    
    mkdir(['laplace' num2str(k)]);
    cd(['laplace' num2str(k)]);
    write_presc(presc);

    presc = TDP;
end

for k=12:2:20;
    
    dpresc = (max(presc.non_zero_values)-min(presc.non_zero_values))./k;
                    
    minpresc = min(presc.non_zero_values);
    
    for n=1:Nstep;
    ind{n} = find(presc.non_zero_values >= (n-1)*dpresc + minpresc ...
        & presc.non_zero_values < n*dpresc + minpresc);
    presc.non_zero_values(ind{n}) = mean(presc.non_zero_values(ind{n}));
    end
    
%     for i=1:k;
%         ind{i} = find(PETnorm.non_zero_values >= (i-1)./k & PETnorm.non_zero_values < i./k);
%         PETnorm.non_zero_values(ind{i}) = double(20*i./k);
%     end

%     PETpresc = PETsparse;
%     PETpresc.non_zero_values = 70.+PETnorm.non_zero_values;
%     PETnorm.non_zero_values = PETsparse.non_zero_values./(max(PETsparse.non_zero_values));


    mkdir(['laplace' num2str(k)]);
    cd(['laplace' num2str(k)]);
    write_presc(presc);
   
    presc = TDP;
end