Topological Image Analysis using Persistence and Min-Max Diagrams

(c) Moo K. Chung

mkchung@wisc.edu

Department of Biostastics and Medical Infomatics
Waisman Laboratory for Brain Imaging and Behavior

Description
July 14, 2010

We will show how to construct persistence diagrms  and min-max diagrams   for 1D (simulated functional signal) and 2D data (cortical thickness). The persistence and min-max diagrams are motivated by topological data analysis where the underlying topology of data and images is quantied.

The codes have been tested under a Macbook pro with 4GB memory and MATLAB 7.5. If you are using the Matlab codes below for your publication, please reference  or . The detailed description of cortical thickness data is given in .

Persistence & Min-Max diagrams for 1D functions
July 14, 2010

The construction for the persistence diagrm for 1D function is based on the iterative pairing and deleation algorithm . Note that the min-max diagram  and the persistence diagram  are identical for 1D functions since there is no saddle points in 1D functions. Let us explain 1D example given in the MICCAI 2009 paper . In the interval [0, 1], we construct a signal s and add noise e to obtain simulated signal Y.

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)

This produces a scatter points so we can't obtain critical values directly. We smooth out the scatter plot using heat kernel smoothing with bandwidth 0.0001 and cosine series expansion upto degree 100  . The resulting smoothed signal is stored in the variable wfs.

k=100; sigma=0.0001;
[wfs, beta]=WFS_COS(Y,x,k,sigma);

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

Then using the iterative pairing and deleation algorithm pairing_1D, we determine the pairing of minimums and maximums.  pairing_1D also plots critical values as white dots (Figure 1).

pairs=pairing_1D(x,wfs);

set(gcf,'Color','w')
``` This produces 3 pairs: ```

>> pairs

pairs =

1.3126    1.6816
0.5948    1.1636
0.2130    0.8886

```The result will be slightly different each time you run the code since the added noise is always changing. The pairing rule is explained in Figure 1.  ```

``` ```Figure 1. The births and deaths of components in sublevel sets. We have critical values a,b,c,d,e,f, where a < b < d <f are minimums and c< e are maximums. At y=a, we have a single component marked by a single gray area. When we increase the level to y=b, we have the birth of a new component in addition to the existing component born at a. At the maximum y=c, the two components merge together to form a single component. Following the pairing rule given in Edelsbrunner (2008), we pair (c,b) and (e,d). Other critical values are paired similarly. See  or .

Cortical thickness data
July 17, 2010

We will use the cortical surface data set used in .  There are total 16 autistic and 11 control subjects. Cortical thickness is computed as the L2 distance between the two surfaces. autism_coord is the gray matter surface for autistic subjects and autism_coordw is the white matter surface for control subjects.

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

We will display cortical thickness of the 1st autistic subject on a unit sphere.

figure_trimesh(sphere,thick_au(1,:), 'rwb')

Since the cortical thickness data is fairly noisy, we will smooth using the weighted-sphearicl harmonic representation with degree 42 and bandwidth 0.001  . The detailed explanation of SPHARMsmooth is given here. The smoothed cortical thickness is given in Figure 2.

directory='/basis/';
SPHARMconstruct(directory,42);

L=42;
sigma=0.001;

coord_co=zeros(11,3,40962);
coord_co(:,1,:)=thick_co;
[coord_sm_co,fourier_coeff_co]= SPHARMsmooth(coord_co,L,sigma);

coord_au=zeros(16,3,40962);
coord_au(:,1,:)=thick_au;
[coord_sm_au,fourier_coeff_au]= SPHARMsmooth(coord_au,L,sigma); Figure 2. Weighted spherical harmonic representation of cortical thickness. Critical values are identified as crosses.

Min-max diagram of a function on a sphere
April 28, 2010

The min-max diagram is constructed by pairing minimums and maximums in a particular fashion . To simplify the problem, we identify how critical values are connected using the Delaunay triangulation. Figure 3. Delaunay triangulation of few minimums and maximums.

Let us map cortical thickness given in Figure 2 onto a plane for better visualization.
fourier_coeff_au.x(:,:,1)is the spherical harmonic coefficients of the 1st autistic subject.

square=SPHARM2square(fourier_coeff_au.x(:,:,1),42,0.001);
[lmax, lmin] = figure_extrema(square); Figure 4. Flatmap represention of cortical thickness given in Figure 2. Critical values are identified as crosses.

The minimum and maximums in the cortical thickness data are paired using the iterative pairing and deleation algorithm . The pairing for all critical values identified in Figure 2 and 4 can be done by

value=
squeeze(coord_au(1,1,:));
pairs= pairing_mesh(sphere, value, L, sigma);

The code requires FINDnbr. Since we are expected to obtain topologically invariant pairing, we constructed the pairing in the spherical mesh.

References
April 22, 2010; October 8, 2010

1. Chung, M.K., Bubenik, P., Kim, P.T. 2009. Persistence diagrams of cortical surface data. Information Processing in Medica Imaging (IPMI). Lecture Notes in Computer Science (LNCS). 5636:386-397.
2. Chung, M.K., Singh, V., Kim, P.T., Dalton, K.M., Davidson, R.J. 2009. Topological characterization of signal in brain images using the min-max diagram. 12th International Conference on Medical Image Computing and Computer Assisted Intervention (MICCAI). Lecture Notes in Computer Science (LNCS). 5762:158-166.
3. Pachauri, D., Hinrichs, C., Chung, M.K., Johnson, S.C., Singh, V. and ADNI. 2010. Cortical surface topology based kernels with application to alzheimer's disease. IEEE Transactions on Medical Imaging. under revision.
4. Chung, M.K., Dalton, K.M., Davidson, R.J. 2008. Tensor-based cortical surface morphometry via weighed spherical harmonic representation. IEEE Transactions on Medical Imaging. 27:1143-1151.
5. Chung, M.K., Dalton, K.M., Shen, L., L., Evans, A.C., Davidson, R.J. 2007. Weighted Fourier series representation and its application to quantifying the amount of gray matter. IEEE Transactions on Medical Imaging 26:566-581.
6. Chung, M.K., Robbins,S., Dalton, K.M., Davidson, Alexander, A.L., R.J., Evans, A.C. 2005. Cortical thickness analysis in autism via heat kernel smoothing. NeuroImage 25:1256-1265