% Statistical and Computational Methods in Brain Image Analysis
% Chapter 14. Persistent Homology
%
% (C) Moo K. Chung 
% University of Wisconsin-Madison
% mkchung@wisc.edu

%--------------------------------------------------------

% Generating scatter points by clicking mouse
figure; 

axis on
hold on
xy = [];
n = 0;
but = 1;
while but == 1
    [xi,yi, but] = ginput(1);
    plot(xi,yi,'ko')
    n = n+1;
    xy(:,n) = [xi;yi];
end

% Rips complex on a key shaped object
endp = mouse2points(500,500, 'toy-key.tif');
imshow('toy-key.tif')
set(gcf,'Color','white','InvertHardcopy','off');

[adj, prob, vertices] = points2graph(endp,50, 'rips');
figure_graph_color(adj,vertices, ones(size(adj,1),1));


%Rips filtration construction

figure;
subplot(2,3,1);
imshow('toy-key.tif')
set(gcf,'Color','white','InvertHardcopy','off');

[adj, prob, vertices] = points2graph(endp,0, 'rips');
figure_graph_color(adj,vertices, ones(size(adj,1),1));
str= ['0 mm']
title(str)
    

for i=1:5
    subplot(2,3,i+1);
    imshow('toy-key.tif')
    set(gcf,'Color','white','InvertHardcopy','off');
    
    [adj, prob, vertices] = points2graph(endp,20*i, 'rips');
    figure_graph_color(adj,vertices, ones(size(adj,1),1));
    str= [num2str(20*i), 'mm']
    title(str)
end;

%--------------------------------------------------------
% Example: 1D function simultation
x=[0:0.002:1]';
s= x + 7*(x - 0.5).^2 + cos(8*pi*x)/2;
e=normrnd(0,0.2,length(x),1);
Y=s+e;
figure;
plot(x,Y,'ko','MarkerEdgeColor',[0.5 0.5 0.5], 'MarkerFaceColor',[0.7 0.7 0.7], 'MarkerSize',4)


% Weighted-cosine series representation
k=100; sigma=0.0001;
[wfs, beta]=WFS_COS(Y,x,k,sigma);

hold on;
plot(x,wfs,'k','LineWidth',5);

% Peristent diagram
pairs=pairing_1D(x,wfs);
set(gcf,'Color','w')
figure; plot(pairs(1,:),pairs(2,:),'.k')

%------------------------------------------------
% Case study:  autism cortical thickness

load AUTISM.coordinates.mat

thick_au=squeeze(sqrt(sum((autism_coord-autism_coordw).^2,2)));

load unitsphere.mat
figure_trimesh(sphere,thick_au(1,:)')
caxis([2 6]); colormap('hot')


% Weighted-SPHARM with 42 degrees and bandwidth 0.001
L=42;
sigma=0.001;

surf.faces=sphere.faces;
surf.vertices=zeros(40962,3);
surf.vertices(:,1)=thick_au(1,:)';

[surf_smooth, fourier]=SPHARMsmooth2(surf,sphere,L, sigma);
figure; figure_trimesh(sphere,surf_smooth.vertices(:,1))
caxis([2 6]); colormap('hot')

% Flattening
square=SPHARM2square(fourier.x,42,0.001);
[lmax, lmin] = figure_extrema(square); colormap('hot')

% Pairing
value=surf_smooth.vertices(:,1);
pairs= pairing_mesh(sphere, value, L, sigma);
figure; plot(pairs(:,1),pairs(:,2),'ko','MarkerEdgeColor',...
    'k', 'MarkerFaceColor',[0.7 0.7 0.7], 'MarkerSize',5)
xlim([2.5 5]); ylim([2.5 5])
print_pdf('PD')




