
% Moo K. Chung mkchung@wisc.edu
% October 2010
%
% Lecture 7.
%
%---------------------------
% persistence diagram

x=[0:0.002:1]';
s=10*(x - 0.5).^2 + cos(7*pi*x)/2 +x
%s= x
%s=zeros(length(x),1);
e=normrnd(0,0.2,length(x),1);
Y=s+e;

figure;
plot(x,Y,'ro', 'MarkerEdgeColor',[0.5 0.5 0.5], 'MarkerFaceColor',[0.5 0.5 0.5], 'MarkerSize',5)
hold on;

k=100; sigma=0.001
% 1D heat kernel smoothing
[wfs, beta]=WFS_COS(Y,x,k,sigma);

plot(x,wfs,'r','LineWidth',2);
set(gcf,'Color','w') 


pairs=pairing_1D(x,wfs)
figure;
plot(pairs(:,1), pairs(:,2),'o', 'MarkerEdgeColor',[0 0 1], 'MarkerFaceColor',[.2 .2 1], 'MarkerSize',10)
set(gcf,'Color','w') 

%---------------------------------
% Bonferroni correction simulation

Y=normrnd(0,1,100,100); 
imagesc(Y); colorbar; colormap('bone')
axis off
set(gcf,'Color','w') 
print -r300 -dtiff N01-1

norminv(0.95,0,1)
% 1.65

1-normcdf(1.64,0,1)
% 0.051

[Yl, Yh] = threshold_image(Y, 1.64);
figure; imagesc(Yh);colormap('bone'); colorbar; 
axis off
set(gcf,'Color','w') 
print -r300 -dtiff N01-2

n=100*100;
size(find(reshape(Y,n,1)>=1.64),1)


ans/n
%0.051

% Bonferroni correction
norminv(1-0.05/10000,0,1)
% 4.42

%--------------------------------------------------
% computing Euler characteristic in 2D binary object

% incorrect way

I=imread('toy-key.tif');
figure;imagesc(I); colorbar
eul = bweuler(I)
% 1


% correct way
I=imread('toy-key.tif');
I= imresize(I,.1, 'nearest');
figure;imagesc(I); colorbar
I=(max(max(I))-I);
I=I/max(max(I));
figure;imagesc(I); colorbar

figure;imshow(I)
axis on;

eul = bweuler(I)
% 0

%-------------------------------
% Monte-Carlo Simulation Method

maxY=[];
n=10000
for i=1:n
  Y=normrnd(0,1,100,100);
  maxY(i)=max(max(Y));
end

hist(maxY,100);

% cumulative distribution function
sortY=sort(maxY);
figure;plot(sortY,[1:n]/n,'-')

% threshold corresponding to alpha=0.05 level
% need to find h such that 0.05*n numbers are above h.
alpha=0.05
sortY(n*(1-alpha))

% 4.40
% compare this to the Bonferroni thresholding
% there is not much difference. why?


%-----------------------------------
% isotropic kernel smoothing

Y=normrnd(0,1,100,100); 
imagesc(Y); colorbar; colormap('bone')
axis off
set(gcf,'Color','w')

smooth = gaussblur(Y,5);
figure;imagesc(smooth);colormap('bone');colorbar;
axis off
set(gcf,'Color','w')
print -r300 -dtiff exc1
%-----------------------------
% Euler characteristic Method

eul=[];
i=1;
for h=-2:0.1:2

    h=0.8
    [Yl, Yh] = threshold_image(smooth, h);
    %figure; imagesc(Yh);colormap('bone'); colorbar;
    %axis off
    %set(gcf,'Color','w')
    %print -r300 -dtiff N01-2

    Yh=reshape(Yh,100*100,1);
    excursion=zeros(100*100,1);
    excursion(find(Yh>h))=1;
    excursion=reshape(excursion, 100, 100);
    figure; imagesc(excursion);colormap('bone');
    axis off
    set(gcf,'Color','w')
    print -r300 -dtiff exc3

    eul(i) = bweuler(excursion);
    i=i+1;
end;
figure;plot(-2:0.1:2, eul,'-k')

%------------------------------
% Mean of Euler characteristic

eul=[];
i=1;
for h=-2:0.05:2
    for k=1:100
        Y=normrnd(0,1,100,100);
        smooth = gaussblur(Y,5);
        [Yl, Yh] = threshold_image(smooth, h);
        Yh=reshape(Yh,100*100,1);
        excursion=zeros(100*100,1);
        excursion(find(Yh>h))=1;
        excursion=reshape(excursion, 100, 100);
        eul(i,k) = bweuler(excursion);
    end;
    i=i+1;
end;

figure;
for k=1:50
hold on; line(-2:0.05:2, eul(:,k),'Color', [0.7 0.7 0.7],'LineWidth',1)
end;
set(gcf,'Color','w')
hold on; plot(-2:0.05:2, mean(eul,2),'-k','LineWidth',2)

print -r300 -dtiff ECplot


%-------------------------------
% 3D thresholding simulation
% with signal


signal=zeros(20,20,20);
signal(1:10,1:10,1:10)=0.5;
%noise =  normrnd(0.5, 1, 20,20,20);
%save noise.oct.2010.mat noise
load noise.oct.2010.mat
data=signal+noise;
data = smooth3(data,'gaussian',5);

for h=0.1:0.3:2
   
    figure; 
    p1 = patch(isosurface(data,h), ...
        'FaceColor',[0.4 0.7 0.8],'EdgeColor','none');
         
    alpha(0.93)

    p2 = patch(isocaps(data,h), ...
        'FaceColor','interp','EdgeColor','none');
    %caxis([0 2])
    isonormals(data,p1)
    view(3); axis vis3d tight
    camlight; lighting phong

    axis off
    set(gcf,'Color','w')
    %colormapeditor

    print('-dtiff', '-r300',strcat('3DEC', num2str(h),'.tif'))

end;





