%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\Geometric;
fw = fopen('Geometric_FQI_prc75.dat','w');

I0 = load('geom000_g240.txt'); I0 = I0*0.254;
I1 = load('geom001_g240.txt'); I1 = I1*0.254;
I2 = load('geom002_g240.txt'); I2 = I2*0.254;
I3 = load('geom003_g240.txt'); I3 = I3*0.254;
I4 = load('geom004_g240.txt'); I4 = I4*0.254;
I5 = load('geom005_g240.txt'); I5 = I5*0.254;

%Set opt = 1 for figures***************************************************
opt = 0;
if opt == 1
    figure(1); imagesc(I0); title('Geometric_0','FontSize',18); colorbar; print -djpeg Geometric_0.jpg;
    figure(2); imagesc(I1); title('Geometric_1','FontSize',18); colorbar; print -djpeg Geometric_1.jpg;
    figure(3); imagesc(I2); title('Geometric_2','FontSize',18); colorbar; print -djpeg Geometric_2.jpg;
    figure(4); imagesc(I3); title('Geometric_3','FontSize',18); colorbar; print -djpeg Geometric_3.jpg;
    figure(5); imagesc(I4); title('Geometric_4','FontSize',18); colorbar; print -djpeg Geometric_4.jpg;
    figure(6); imagesc(I5); title('Geometric_5','FontSize',18); colorbar; print -djpeg Geometric_5.jpg;
    pause;
end
%**************************************************************************

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>0);
    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 HS_Mean = %f\t  HS_Std = %f\n',prc,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);

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;

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,'Elapsed Time = %d\n',ET);

fclose all;