%Comparison of Geometric Cases
%Last updated: 12/5/2007 (Sukanta Basu)
clear all; close all; clc;

%Works for Matlab Release 2007a - Multi-threading**************************
%setNumberOfComputationalThreads(2);
%Works for Matlab Release 2007b - Multi-threading**************************
maxNumCompThreads(2);
%**************************************************************************
T0 = clock;

%cd C:\Sukanta\Laptop\Works_Ongoing\QPF\Cases\Fake\Geometric;
cd D:\Works_Ongoing\QPF\Cases\Fake\Perturbed;
fw = fopen('Perturbed_FQI_prc75_Th1.dat','w');

Th = 1;
I0 = load('fake000_g240.txt'); I0 = I0*0.254; 
I1 = load('fake001_g240.txt'); I1 = I1*0.254;  
I2 = load('fake002_g240.txt'); I2 = I2*0.254; 
I3 = load('fake003_g240.txt'); I3 = I3*0.254; 
I4 = load('fake004_g240.txt'); I4 = I4*0.254; 
I5 = load('fake005_g240.txt'); I5 = I5*0.254; 
I6 = load('fake006_g240.txt'); I6 = I6*0.254; 
I7 = load('fake007_g240.txt'); I7 = I7*0.254; 

%Set opt = 1 for figures***************************************************
opt = 0;
if opt == 1
    figure(1); imagesc(I0); title('Perturbed_0','FontSize',18); colorbar; caxis([0 20]); print -djpeg Perturbed_0.jpg;
    figure(2); imagesc(I1); title('Perturbed_1','FontSize',18); colorbar; caxis([0 20]); print -djpeg Perturbed_1.jpg;
    figure(3); imagesc(I2); title('Perturbed_2','FontSize',18); colorbar; caxis([0 20]); print -djpeg Perturbed_2.jpg;
    figure(4); imagesc(I3); title('Perturbed_3','FontSize',18); colorbar; caxis([0 20]); print -djpeg Perturbed_3.jpg;
    figure(5); imagesc(I4); title('Perturbed_4','FontSize',18); colorbar; caxis([0 20]); print -djpeg Perturbed_4.jpg;
    figure(6); imagesc(I5); title('Perturbed_5','FontSize',18); colorbar; caxis([0 20]); print -djpeg Perturbed_5.jpg;
    figure(7); imagesc(I6); title('Perturbed_6','FontSize',18); colorbar; caxis([0 20]); print -djpeg Perturbed_6.jpg;
    figure(8); imagesc(I7); title('Perturbed_7','FontSize',18); colorbar; caxis([0 20]); print -djpeg Perturbed_7.jpg;   
    pause;
end
%**************************************************************************
I0 = I0.*(I0>Th);
I1 = I1.*(I1>Th);
I2 = I2.*(I2>Th);
I3 = I3.*(I3>Th);
I4 = I4.*(I4>Th);
I5 = I5.*(I5>Th);
I6 = I6.*(I6>Th);
I7 = I7.*(I7>Th);

Nsurrogates = 10;                  %Number of surrogates
prc         = 0.75;                %Percentile for PHD calculation
HS          = zeros(size(Nsurrogates));

for i = 1:Nsurrogates
    I0x = MakeSurrogate(I0); 
    I0x = I0x.*(I0x>Th);
    opt = 0;
    if opt == 1
        figure(1); subplot 121; imagesc(I0);  axis square;
                   subplot 122; imagesc(I0x); axis square; 
                   pause(1);
    end
    HS(i) = Hausdorff(I0,I0x,prc);
end
HS_Mean = mean(HS); HS_Std = std(HS);
fprintf(fw,'Percentile = %d\t Threshold = %d\t HS_Mean = %f\t  HS_Std = %f\n',prc,Th,HS_Mean,HS_Std);

[H1,Q1,Q1x] = Hausdorff(I0,I1,prc); fprintf(fw,'H1 = %f\t Q1 = %f\t Q1x = %f\n',H1,Q1,Q1x);
[H2,Q2,Q2x] = Hausdorff(I0,I2,prc); fprintf(fw,'H2 = %f\t Q2 = %f\t Q2x = %f\n',H2,Q2,Q2x);
[H3,Q3,Q3x] = Hausdorff(I0,I3,prc); fprintf(fw,'H3 = %f\t Q3 = %f\t Q3x = %f\n',H3,Q3,Q3x);
[H4,Q4,Q4x] = Hausdorff(I0,I4,prc); fprintf(fw,'H4 = %f\t Q4 = %f\t Q4x = %f\n',H4,Q4,Q4x);
[H5,Q5,Q5x] = Hausdorff(I0,I5,prc); fprintf(fw,'H5 = %f\t Q5 = %f\t Q5x = %f\n',H5,Q5,Q5x);
[H6,Q6,Q6x] = Hausdorff(I0,I6,prc); fprintf(fw,'H6 = %f\t Q6 = %f\t Q6x = %f\n',H6,Q6,Q6x);
[H7,Q7,Q7x] = Hausdorff(I0,I7,prc); fprintf(fw,'H7 = %f\t Q7 = %f\t Q7x = %f\n',H7,Q7,Q7x);

FQI1 = (H1/HS_Mean)/Q1;     FQI1x = (H1/HS_Mean)/Q1x;
FQI2 = (H2/HS_Mean)/Q2;     FQI2x = (H2/HS_Mean)/Q2x;
FQI3 = (H3/HS_Mean)/Q3;     FQI3x = (H3/HS_Mean)/Q3x;
FQI4 = (H4/HS_Mean)/Q4;     FQI4x = (H4/HS_Mean)/Q4x;
FQI5 = (H5/HS_Mean)/Q5;     FQI5x = (H5/HS_Mean)/Q5x;
FQI6 = (H6/HS_Mean)/Q6;     FQI6x = (H6/HS_Mean)/Q6x;
FQI7 = (H7/HS_Mean)/Q7;     FQI7x = (H7/HS_Mean)/Q7x;

ET   = etime(clock,T0);

fprintf(fw,'FQI1 = %f\t FQI1x = %f\n',FQI1,FQI1x);
fprintf(fw,'FQI2 = %f\t FQI2x = %f\n',FQI2,FQI2x);
fprintf(fw,'FQI3 = %f\t FQI3x = %f\n',FQI3,FQI3x);
fprintf(fw,'FQI4 = %f\t FQI4x = %f\n',FQI4,FQI4x);
fprintf(fw,'FQI5 = %f\t FQI5x = %f\n',FQI5,FQI5x);
fprintf(fw,'FQI6 = %f\t FQI6x = %f\n',FQI6,FQI6x);
fprintf(fw,'FQI7 = %f\t FQI7x = %f\n',FQI7,FQI7x);
fprintf(fw,'Elapsed Time = %d\n',ET);

fclose all;