1
Finding Optimal Brain Mappings Using Integer Linear Programming Solvers and Other Methods
Heriberto Flores
Advisor: Unal Sakoglu, Ph.D.
Department of Computer Science
Texas A&M University – Commerce 2
Abstract
Functional magnetic resonance imaging (fMRI) datasets, which are snapshots of the brain taken every few seconds, are four dimensional (4D) datasets, in which a three-dimensional structure contains the structural/spatial brain data and the remaining 1D is the temporal data. Before any analysis, 4D data is converted into 2D, spatial and temporal dimensions. This involves mapping of 3D brain MRI data into a one dimensional array. For almost all of the fMRI data, this reordering has been done using simple linear ordering, which is not optimal in the sense that it leads to a highly disconnected and scattered human brain in 1D, where the structure information is lost. Currently, there is no optimal way of mapping the data from 3D to 1D without the loss of this structural information [6].
The purpose of this project has been to find an optimal way of mapping brain data through the study of two methods. Initially, we analyze the availability and function of integer linear programming solvers as a possible way of finding an optimal solution by transforming brain data into an integer linear programming problem instance. Additionally, we investigated the alternate use of the Hilbert space-filling curve as a method of finding an optimal solution in optimal time. The optimal solution, in this case, is expected to retain most of the structural connection information and provide a better representation of the brain in 1D.
3
Acknowledgments:
This research project would not have been possible without the contribution of a few individuals to whom I would like to say thank you to. First and foremost, I would like to thank my advisor, Dr. Unal Sakoglu, for having allowed me to work under his supervision during the duration of this project. Without his contribution, leadership, and supervision this would have been practically impossible. Secondly, I would like to say thank you to Dr. Abdullah Arslan for his contribution and guidance to the project. Also, I would like to say thank you to Kushal Bohra for helping me get acquainted with the beginnings of this project and answering many of the questions I had. Also, Dr. Mutlu Mete for the access to his equipment at the beginning of the project, without it, it would have taken us longer to begin our trials. Lastly, I would like to thank the Computer Science department and the Honors College program for having allowed me the opportunity to complete a project like this.
4
1. Introduction
1.1 Problem Statement
The objective of the project is to find an optimal mapping solution in order to convert a 3D template MRI brain into a 1D brain without losing any of the spatial information. The initial proposal established the goal of formulating and minimizing an integer linear programming (ILP) instance as a method of finding an optimal mapping solution. The ILP instance created in our theory is based on the theories of the Hamiltonian Path Problem (HPP), and the Traveling Salesman Problem (TSP). HPP and TSP are non-deterministic polynomial-time-hard (NP-hard) problems in combinatorial optimization; in other words, it takes a non-deterministically long time to compute and optimize. Since a typical 3D MRI brain contains more than 100,000 volume elements (voxels), the optimization problem instance can contain hundreds of thousands of variables making it practically impossible to compute a solution without the aid of large computer clusters.
The considerable size of the problem instance has also generated an additional set of limiting factors for which previous research and tests have been conducted in efforts of increasing efficiency. The multidimensional data being utilized for this project consists of a structural variable derived from a location in the brain which is represented by a three-dimensional coordinate set {x,y,z}, as well as a temporal variable representing the signal value {s} of the respective voxel identified by the corresponding coordinate set [6]. The brain dataset used has been divided into 64x64x64, or 262,144 voxels. It is imperative to recognize that some of these voxels, however, are considered “fillers”; that is they do not contain data 5
important to the solution, but are there, instead, just to create a perfect cubed figured necessary for the work. As a result, the amount of voxels involved creates a large number of problem constraints which increases the difficulty of finding an optimal solution. Working with large datasets like the one mentioned here is the main reason for the need of specialized software in this project.
1.2 Purpose of the Project
The eventual objective of the project is to develop an instance of an integer linear programming problem (ILP), and use an ILP solver to find an optimized solution for a large brain dataset. In this work, however, the focus is on developing and testing the theory with a smaller, solvable 3D “toy” brain structure. This is done in order to speed up the process and analyze possible results as the original problem instance could take an incalculable amount of time to solve due to the large number of constraints. By using a “toy” case we can manipulate the data involved, analyze, and better understand whether the results are correlated to our expected outcome. Once an optimal solution is found for the toy case the same approach can be applied to formulate and solve the problem instance for larger datasets. The mapping methods found in this study do not only have to apply to brain mapping, but can be used to map any other three-dimensional sets of data into a single dimension for different types of analysis.
1.3 Research Methodology
The first step in this project involved finding an adequate mixed integer linear programming (MILP) solver that would work with our datasets. A MILP is a solver that is able to solve both mixed integer and integer linear programming problems. The biggest 6
obstacle in this step is finding a solver that is able to handle problems with a large number of variables. There are many solvers which limit the amount of variables that can be used. We selected LP_solve as the solver because it has been deemed the most adequate in previous tests. The solver was tested with simpler problem instances containing in excess of 300,000 variables and provided satisfying results.
LP_solve is a MILP solver [5] which can be run in three different environments. It may be used in raw format through the command prompt, as a library in variety of programming languages, and as a standalone solver with its own IDE. In this instance, LP_solve will be used through the command prompt as this provides the solver with the necessary power to handle large variable problems. The solver is designed to use a text file with the ILP problem instance as the input, and outputs the solution into a separate text file. The solution text file contains the order in which each of the voxels is visited. This solution essentially defines the best fit curve for the dataset.
The main goal has been to initially use the solver to find an optimal solution for multiple toy cases. This is done, in part, because smaller cases provide solutions which can be visually analyzed without the need for special software. The cases are small enough where a solution can be found with a simple look at the data. Once the results returned by the solver are satisfying enough then the ILP problem instance for the main dataset can be introduced. The important step here was to make sure that LP_solve would return optimal solutions for the smaller toy cases. Once the main problem instance is fed into the solver then it becomes a waiting game which can last up to months. Eventually, the goal is for the solver to return the order of the best fit curve for the original data which can then be used to analyze the 3D data in a single dimension. 7
2. Integer Linear Programming
2.1 Designing the Problem Instance
As stated above, LP_solve only needs a simple text file with the ILP problem on it in order to find an optimal solution. The solver is designed to read the problem and take into consideration the constraints created in order to deliver results. One of the most challenging aspects of this problem was transforming the multidimensional MRI dataset into an instance of an ILP problem. The dataset is delivered as a 3D (x,y,z) spatial dimension and a single dimensional variable (s) containing the temporal information, thus the original dataset is in the format of (x,y,z,s). The transformation of this data into the ILP problem instance was done using Java programming.
The program written is designed to create an ILP problem instance based on a modified version of the Hamiltonian Path Problem (HPP) known as the Traveling Salesman Problem (TSP). TSP develops an instance in which every point, or voxel, in the dataset must be visited once and will making a full round trip around the voxels. That is, the solution begins and ends at the same location. The ILP formulation in this case can be better expressed by the following equation:
Minimize: ΣΣ
[6] 8
Subject to: ( ) Σ { } Σ { }
,
This function creates an instance of an ILP problem which can then be used by LP_solve to formulate a solution. Creating the instance in Java is no easy task as the accuracy of the instance created has to be assured before proceeding to attempt to solve the main problem. This has to be done due to the large size of the original instance. The large number of variables and constraints in the original problem instance are too large to be able to look for bugs. The main purpose of the toy case is to test the program’s accuracy and the solver’s effectiveness.
2.2 Programming
Java was used as the main programming language in order to build the ILP problem instance. The goal of the program is to create an accurate ILP problem instance representation of the brain dataset. In order to do this the program must take into account each voxel individually as well as considering neighbors. Figures 1 and 2 show part of the code used to create the problem instance:
9
144 problem.print("min: ");
145
146 for( int i = 0; i <= dataSize; i++){
147 for(int j = 0; j <= dataSize; j++){
148 for( int k = 0; k <= dataSize; k++){
149
150 long mask = flag.getMask(i,j,k);
151
152 int index1 = 0;
153 int tempIndex = 0;
154
155 float signal1 = 0;
156 float tempSignal = 0;
157
158 float tempCoef = 0;
159 double coef = 0;
160
161 boolean there = maskedValues.contains(mask);
162 long addr = mask;
163
164 //index1 = maskedValues.indexOf(mask);
165
166 long temp1 = 0, temp2 = 0, temp3 = 0, temp4 = 0, temp5 = 0, temp6 = 0, temp7 = 0, temp8 = 0;
167
168 //Checking and adding possible path address values to a list for building the problem, will have to build the problem
169 //instance here in order to add the coefficients corresponding to the signal values.
170
171 if(there){
172
173 index1 = address.indexOf(mask);
174
175 signal1 = signalValue.get(index1);
176
177 if(maskedValues.contains(flag.getMask(i+1,j,k))){
178 temp1=addr+1;
179 pathAddress.add(addr+1);
180
181 tempIndex = address.indexOf(flag.getMask(i+1,j,k));
182 tempSignal = signalValue.get(tempIndex);
183
184 tempCoef = (float)Math.pow((signal1 - tempSignal), 2);
185
186
187 problem.printf(" + %f x%d",tempCoef, temp1);
188 }
189
190 if(maskedValues.contains(flag.getMask(i-1,j,k))){
191 temp2 = addr+2;
192 pathAddress.add(addr+2);
193
194 tempIndex = address.indexOf(flag.getMask(i-1,j,k));
195 tempSignal = signalValue.get(tempIndex);
196
197 tempCoef = (float)Math.pow((signal1 - tempSignal), 2);
198
199 problem.printf(" + %f x%d",tempCoef, temp2);
200 }
201
202 if(maskedValues.contains(flag.getMask(i,j+1,k))){
203 temp3 = addr+3;
204 pathAddress.add(addr+3);
205
206 tempIndex = address.indexOf(flag.getMask(i,j+1,k));
207 tempSignal = signalValue.get(tempIndex);
208
Figure 1 10
209 tempCoef = (float)Math.pow((signal1 - tempSignal), 2);
210
211 problem.printf(" + %f x%d",tempCoef, temp3);
212 }
213
214 if(maskedValues.contains(flag.getMask(i,j-1,k))){
215 temp4 = addr+4;
216 pathAddress.add(addr+4);
217
218 tempIndex = address.indexOf(flag.getMask(i,j-1,k));
219 tempSignal = signalValue.get(tempIndex);
220
221 tempCoef = (float)Math.pow((signal1 - tempSignal), 2);
222
223 problem.printf(" + %f x%d",tempCoef, temp4);
224 }
225
226 if(maskedValues.contains(flag.getMask(i,j,k+1))){
227 temp5 = addr+5;
228 pathAddress.add(addr+5);
229
230 tempIndex = address.indexOf(flag.getMask(i,j,k+1));
231 tempSignal = signalValue.get(tempIndex);
232
233 tempCoef = (float)Math.pow((signal1 - tempSignal), 2);
234
235 problem.printf(" + %f x%d",tempCoef, temp5);
236 }
237
238 if(maskedValues.contains(flag.getMask(i,j,k-1))){
239 temp6 = addr+6;
240 pathAddress.add(addr+6);
241
242 tempIndex = address.indexOf(flag.getMask(i,j,k-1));
243 tempSignal = signalValue.get(tempIndex);
244
245 tempCoef = (float)Math.pow((signal1 - tempSignal), 2);
246
247 problem.printf(" + %f x%d",tempCoef, temp6);
248 }
249
250 pathAddress.add(addr+7); // Path to z
251 pathAddress.add(addr+8); // path from z
252
253 temp7 = addr+7;
254 temp8 = addr+8;
255
256 float maskSignal = signalValue.get(address.indexOf(mask));
257
258 float toZero = (float)Math.pow((maskSignal-1), 2);
259 float fromZero = (float)Math.pow((1-maskSignal), 2);
260
261 //String message = "This is not working!";
262
263 problem.print(" + " + toZero + " x" + temp7);
264 problem.print(" + " + fromZero + " x" + temp8);
265 // problem.printf("%s%n", message);
266
267 }
268 else continue;
269 }
270 }
271 }
272
Figure 2 11
In the toy case scenario, there are 27 voxels being used. The program must take into account the uniqueness of each one and the expectations from the problem. Every voxel must be masked into a unique variable that can be used in the program and later traced back. Paths between each neighbor must also be identified and given a unique identifier, which can vary depending on direction. All of these values become the variables and constraints of the ILP problem instance. Once the problem is built it can be fed into the solver.
The solver is designed to return a file with the list of the variables corresponding to the voxel locations as well as a list of the variables corresponding to the paths between voxels. In this case, variables representing voxels are masked as u####, where #### represents a unique identifier assigned to each voxel derived from its (x,y,z) coordinates. The following algorithm was used to mask the voxels location into the unique variable identifier.
( ) + z
After the mask for a certain voxel is identified, a second mask is assigned to the paths between each voxel and its neighbors. These variables are identified in the problem instance as x####, where #### is also a unique identifier pertaining to each path. It is worth noting, however, that the path between two neighbors will eventually have two identifiers depending on which direction the path is being traveled.
When results are returned by the solver, then a simple transformation of the variables to return them into their (x,y,z) format is needed in order to analyze them. The solver provides the order in which each voxel was traversed. Essentially this order is the equivalent of the optimal solution to the problem. This order returned can then be applied to similar datasets in order to analyze them. 12
2.3 Toy Case
The purpose in using toy cases, as stated before,
was to provide a quick way of obtaining results that
could be easily analyzed. The toy cases used in this
project were two small sets of data designed as
scaled down versions of the much larger dataset. The
first toy case used was a 3x3x3 sample. The toy case,
made up of the data in Table 1, was a fictitious brain
made up of only 7 voxels. If stretched out, this
sample brain would resemble a straight line in which
each voxel is only connected to two other neighbors.
This limits the amount of choices that the solver has
while deciding which path to travel and greatly
reduces the number of constraints created in the
problem. This is done in order to test LP_solve’s
ability to traverse the brain in one continuous path.
In order to understand how the toy case is
designed it is important to remember that the data represents two aspects of a voxel. The first
three columns represent the location of each individual voxel while the third column
represents the temporal data from that corresponding voxel. As it might be noted from Table
1, some of the values of S are equivalent to 0. This is because the brain has been padded with
“filler” voxels in order to create a perfect cube. The voxels with a value of 0 are typically
located outside of the brain. The expected outcome from the solver is to find a solution that
X Y Z S
0 0 0 1
1 1 1 0.948134
1 1 2 0
1 1 3 0
1 2 1 0.956256
1 2 2 0.942653
1 2 3 0.625394
1 3 1 0
1 3 2 0
1 3 3 0
2 1 1 0.968796
2 1 2 0.945678
2 1 3 0.929487
2 2 1 0
2 2 2 0.963415
2 2 3 0
2 3 1 0
2 3 2 0
2 3 3 0
3 1 1 0
3 1 2 0
3 1 3 0
3 2 1 0
3 2 2 0.985685
3 2 3 0.998546
3 3 1 0
3 3 2 0
3 3 3 0.954822
3 by 3 by 3 Toy Brain
Table 1: 3 by 3 by 3 Toy Brain Data
13
traverses every voxel for which the value of S does not equal 0. When found, the solution provided by LP_solve can be considered as the optimal solution to the problem.
The following figures depict the results returned by the solver in the toy cases. As observed, the solver found the path between all the voxels. The solver actually began in one corner of the cube, traveled through the empty voxels until it found the brain voxels, and then continued to travel through the brain before finishing up with the “filler” voxels. Figure 1 shows that path that the solver traveled hitting every brain voxel continuously in the 3x3x3 case. Figure 2 shows the same type of results with a toy brain of 4x4x4. The solutions found here help ensure that the solver and the program are working as expected.
Figure 3: 3x3 Toy Brain
14
Figure 4: 4x4 Toy Brain
2.4 Future Problem
The next step in the project is to develop the ILP problem instance for the main dataset. The main obstacle in this step is the sheer size of the problem and the inability to obtain intermediate results. This means that the problem must finish running before any results may be analyzed. It also means that there is no way of knowing whether the solver is actually trying to find a solution or if there is an issue with an instance that large. We have attempted to solve the problem many times; however, the lack of computer power at the beginning of the project and the lack of any type of intermediate results has made it nearly impossible to sit and wait for results. The solver ran continuously for 3 months and no solution was obtained. 15
The lack of a solution does not, however, mean that there is not a solution to the case, as there are factors that may have contributed to the solver not delivering results. First, the solver may be limited on what it can handle. LP_solve is a community driven solver which means that it only improves as users of the program contribute. The solver may actually have encountered a limit with all of the constraints involved. There are more powerful solvers out there; however, commercial solvers licenses can run into thousands of dollars.
Instead of waiting for an optimal solution that may or may not exist, we decided to test a suboptimal way of mapping brain data by using the Hilbert space-filling curve. This method, although previously tested, provides an optimal way of mapping brain data into a single dimension and allows for the analysis of such data.
16
3. Hilbert Space-Filling Curve
3.1 Using Hilbert Curve
As stated earlier, finding an optimal solution using the integer linear programming method proposed here can take exponential time, and, at the same time, no solution has been guaranteed. In order to further visualize large 3D datasets we have elected to include a suboptimal method of mapping. This suboptimal method includes finding the Hilbert space-filling curve and applying it to multiple brain activation datasets in order to classify the data. The Hilbert curve is applied to 16 individual datasets representing two different groups, a control group and a patient group. The goal was to find the suboptimal mappings of the brain data and then use those mappings in order to analyze and classify the data. The Hilbert space-filling curve has proven an optimal method of mapping [4].
Figure 5 shows the brain activation map for a member of the control group while Figure 6 shows a similar map for a member of the patient group. These activation maps were analyzed using mricron, which is software that allows us to make visual observations of different areas of the brain data at multiple angles. We can navigate the data by essentially pin pointing the location which we would like to explore. These maps record activation signals in the grey matter of the brain. The grey matter areas are essentially the portions of the brain involved in sensory reception and other vital functions. These images help visualize the density of grey matter in sections of the brain, the loss of this density can and has been attributed to certain diseases before such as Alzheimer’s [3]. The white parts shows show a higher activation signal than the darker sections. 17
AQI_C:
Figure 5: Activation map of a control subject labeled AQI
18
CHD_P:
Figure 6: Activation map of a patient subject labeled CHD
The computation of the Hilbert space filling curve is widely available online. Matlab has functions which allows us to compute either two dimensional or three dimensional Hilbert curves of certain sizes. In our case, we wanted to compute a Hilbert curve over a three dimensional space containing 64x64x64, or 262,144 voxels. Once the Hilbert curve is found over the 3D 19
Table 2: Patient and Control subjects
space, the curve can be applied datasets of the same dimensions. These datasets can then be analyzed in a single dimension plot in order to understand the data better.
We applied the Hilbert curve to the datasets from Table 2. These datasets contains MRI activation data for 16 different subjects. The 16 subjects are divided into two groups, a control and a patient group. Once the Hilbert curve is applied to the data we then transform the data into a single dimensional array using the order in which the Hilbert curve has placed each voxel. This array can then be plotted and the data analyzed to find significant similarities or differences.
3.2 Finding the Hilbert Curve
In order to find the Hilbert curve for a 64x64x64 cube we used the Matlab programming software. This software environment allows us to visualize, analyze, and essentially build algorithms to work with our data. We used the function hilbert3(n), with a value of n = 6, in order to find a curve fitted for a 643 cube. Since this function gives us the Hilbert curve for an area centralized on (0,0,0) we had to expand the script in order to find the curve traversing the data points matching our datasets. The script used for this data is shown on figures 9, 10, and 11.
However, before the Hilbert curve could be applied to the datasets we had to normalize our datasets into a 64x64x64 dimension. The datasets were not available in a perfect cube, therefore, we had to normalize the data by first reducing the
Brain Data#NameGroup1AQI_CC2CHD_PP3CIF_CC4CML_CC5CTT_PP6CVX_PP7DAY_PP8DFJ_CC9DSN_PP10DUQ_PP11DUT_PP12EWR_CC13FHE_CC14FLK_PP15GWW_CC16HCX_CC20
size of the datasets to a size that would fit into a 64x64x64 cube. Since the data is does not make up a perfect cube any empty spaces had to be filled with 0s. This is done in order to easily compare datasets with each other. Once the datasets are normalized the Hilbert curve can then be applied.
The script used begins by reading in the brain activation maps that have already been normalized. Since the Hilbert curve script gives the coordinates in vector coordinates we have to change them into a format which will work with our data by manipulating the coordinates given. The script also fills in the area outside the brain with zeros. It then computes the average for every certain number of values depending on the bin size and stores this value in an array for the classification.
The Hilbert curve is returned through an array as stated above which includes the coordinates of the points ordered in the order in which they were visited. This array is then used to apply the Hilbert curve to each of the datasets available. This is done by rearranging each dataset in the order of the Hilbert curve. Each of these datasets is then saved into an array containing all of the rearranged for classification. It then creates figures depicting the data so it can be visually analyzed. Figure 7 shows the Hilbert curve used for this project. The curve travels every point in the cube, Figure 8 shows a closer look at how it passes through each point. 21
Figure 7: Hilbert curve traversing a 64 by 64 by 64 cube
Figure 8: Zoomed in version of Hilbert ordering 22
datafolder='C:\Users\Eddie\Documents\Thesis\BrainActivationMatlab'
cd(datafolder);
BinSize = 1000;
% This is the hilbert curve that will be used to classify the brain data,
% the hilber curve traverses a 64 by 64 by 64 cubed.
% The hilbert function used here is predefined.
n = 6;
[x,y,z] = hilbert3(n);
originalHilbertCurve = [x;y;z];
% Names of the brain data files to be used
FileNames = ...
{'wzstat1_AQI_C_in_MNI',
'wzstat1_CHD_P_in_MNI',
'wzstat1_CIF_C_in_MNI',
'wzstat1_CML_C_in_MNI',
'wzstat1_CTT_P_in_MNI',
'wzstat1_CVX_P_in_MNI',
'wzstat1_DAY_P_in_MNI',
'wzstat1_DFJ_C_in_MNI',
'wzstat1_DSN_P_in_MNI',
'wzstat1_DUQ_P_in_MNI',
'wzstat1_DUT_P_in_MNI',
'wzstat1_EWR_C_in_MNI',
'wzstat1_FHE_C_in_MNI',
'wzstat1_FLK_P_in_MNI',
'wzstat1_GWW_C_in_MNI',
'wzstat1_HCX_C_in_MNI'}
%tempFile = FileNames(1);
% First, transform the hilbert curve provided above into usable variables
% L is the length of the array, that is the number of variables in the
% hilbert curbe determined by n
L = (2^n)^3;
maxNum = ((L/2-1)+0.5)*(1/L);
minNum = -((L/2-1)+0.5)*(1/L);
% Transforming the curve
tempHilbertCurve = round (( originalHilbertCurve ./(((L/2-1)+0.5)*(1/L))) * 32 + 32 );
newHilbertCurve = tempHilbertCurve;
x_temp = tempHilbertCurve(1,:);
y_temp = tempHilbertCurve(2,:);
z_temp = tempHilbertCurve(3,:);
%figure,hist(x_temp,L*10);
%figure,hist(y_temp,L*10);
%figure,hist(z_temp,L*10);
Figure 9 23
for i=1 : size(tempHilbertCurve,1),
for j=1 : size(tempHilbertCurve,2),
if tempHilbertCurve(i,j) < 33,
newHilbertCurve(i,j) = newHilbertCurve(i,j)+1;
end
end
end
x_new = newHilbertCurve(1,:);
y_new = newHilbertCurve(2,:);
z_new = newHilbertCurve(3,:);
% This is the plot of the new hilbert curve
figure,plot3(x_new, y_new, z_new);xlabel('x'),ylabel('y'),zlabel('z'),
% Once the values of the hilbert curve are found, the following plots and
% the data in linear form first, then the same data is plotted using the
% hilbert curve above
for iFile = 1:size(FileNames,1),
FileRoot = char(FileNames(iFile));
FileName = char(strcat(FileNames(iFile), '.nii'));
my3DMatrix = spm_read_vols(spm_vol(FileName)); %52,63,52
my3DMatrix(find(isnan(my3DMatrix))) = 0;
[X,Y,Z] = size(my3DMatrix);
my3DMatrix_zeropadded = zeros(64,64,64);
my3DMatrix_zeropadded(1:X,1:Y,1:Z) = my3DMatrix;
FileNameCubedMat = char(strcat(FileNames(iFile), '_64cubed.mat'));
save(FileNameCubedMat, 'my3DMatrix_zeropadded')
my3DMatrix_zeropadded_1D = reshape(my3DMatrix_zeropadded, 64^3,1);
VectorLength = length(my3DMatrix_zeropadded_1D);
DecimateRatio = round(VectorLength/BinSize);
for iBin=1:DecimateRatio-1,
my3DMatrix_zeropadded_1D_binned(iBin)= mean(my3DMatrix_zeropadded_1D((iBin-1)*BinSize+1:iBin*BinSize));
end
FileName1DMat = char(strcat(FileNames(iFile), '_1D.mat'))
save(FileName1DMat, 'my3DMatrix_zeropadded_1D')
FileName1DBinnedMat = char(strcat(FileNames(iFile), '_1D_Binned.mat'))
save(FileName1DBinnedMat, 'my3DMatrix_zeropadded_1D_binned')
size(my3DMatrix_zeropadded_1D_binned)
figure,subplot(211),plot(my3DMatrix_zeropadded_1D), title([FileRoot ' in 1D -linear'])
subplot(212),stem(my3DMatrix_zeropadded_1D_binned,'.'),title([FileRoot ' in 1D -linear, binned'])
saveas(gcf, FileRoot, 'fig')
saveas(gcf, FileRoot, 'png')
% Now create the file with the hilbert curve information for the same
% data, as well as a figure depicting the data using the hilbert curve
% Plotting the hilbert curves
Figure 10 24
Lastly, after the Hilbert curve was applied to the datasets we decided to bin the data into different sizes in order to have a closer observation. Since the data contains over 200,000 thousand data members, the process of binning allows us to group the data into different size groups so that we may look and focus on major areas of interest. There are too many single members and looking at every 100 or 1000 can give us the same representation as a smaller group.
for index = 1:size(newHilbertCurve, 2),
myHilbert1DMap(index) = my3DMatrix_zeropadded( newHilbertCurve(1,index),newHilbertCurve(2,index),newHilbertCurve(3,index));
end
for iBin=1:DecimateRatio-1,
myHilbert1DMap_binned(iBin)= mean( myHilbert1DMap((iBin-1)*BinSize+1:iBin*BinSize));
end
% Saving the variables containing the data
FileNameHilbert1DMat = char(strcat(FileNames(iFile), '_Hilbert1D.mat'))
save(FileNameHilbert1DMat, 'myHilbert1DMap')
FileNameHilbertBinnedMat = char(strcat(FileNames(iFile), '_HilbertBinned.mat'))
save(FileNameHilbertBinnedMat, 'myHilbert1DMap_binned')
% Plotting the hilbert curves
figure,
subplot(211), plot(myHilbert1DMap);title( [FileRoot ' in Hilbert Curve 1D']);
subplot(212), plot(myHilbert1DMap_binned,'.'); title( [FileRoot ' in Hilbert Curve binned']);
saveas(gcf, [FileRoot '_Hilbert'], 'fig')
saveas(gcf, [FileRoot '_HilbertFigure'], 'png')
end
Figure 11 25
4. Results
4.1 Analysis
One important step in analyzing all of the sixteen datasets in this project is to be able to compare them to each other. The goal is to find characteristics between members of the same groups that can be used to characterize other members. Figures 12-15 show the datasets arranged alongside similar group members. The goal here is to try to visually identify any characteristics that only show on all members of the same group, and which are likely absent from some, if not all, members of the other group. Figures 12 and 13 show the activation maps from both, the control group and the patient group, arranged in linear ordering. Figures 14 and 15 show these exact same two groups arrange according to the Hilbert ordering. It is important to remember that these figures are a representation of thousands of values, and, therefore, it is extremely difficult to analyze them visually.
One of the main characteristics to look for in these figures is the shape of the curves between members of the same groups. The shape of these curves is not necessarily always similar, however, they do show the areas where activation signals are being reduced or where there is still areas of high activation. While looking at the figures may not fully reveal the characteristics that we are looking for, they can certainly help by pin pointing the location where such clues may be found. 26
Figure 12: The eight Control datasets in linear ordering
Figure 13: The eight Patient datasets in linear ordering 27
Figure 14: The eight Control datasets in Hilbert ordering
Figure 15: The eight Patient datasets in Hilbert ordering 28
Once the Hilbert curve was applied we began to analyze the data obtained from the reordering. Figure 16 shows the datasets ordered in linear form while Figure 17 shows the same datasets ordered using the Hilbert curve. These figures show the activation intensity of the sixteen datasets in two different orderings. Both of these figures picture all of the datasets lined up along similar imaginary lines. For example, in theory while analyzing Figure 17 you should be able to draw a vertical straight line all the way across the sixteen datasets. This line would be passing through the same location in all subjects therefore allowing us to identify the activation patterns along the same area across all subjects. After taking a closer look we can clearly observe that the datasets in Figure 17 are ordered in a way in which the clusters formed by the Hilbert ordering are lined up with each other across all datasets. This allows us to visually analyze the same area throughout the datasets while utilizing these clusters in order to identify any characteristics that can essentially help with the classification of the data. 29
Figure 16: Datasets in linear form 30
Figure 17: Activation maps datasets arranged using Hilbert curve ordering.
31
After looking at all of the datasets as a group it is often helpful to look at the data individually in order to focus on one dataset at a time. The following figures show the first dataset in the groups, belonging to subject AQI_C. Figure 18 shows the data in linear form both by showing all of the data values, and then, by binning the data every into every 1000 values. From this ordering we can see a general shape of the curve and where the data has positive and negative values, however, we do lack a sense of location. We do not know where most of these spikes are located within a section of the brain. The data is ordered by the way it is read and therefore lacks any location information.
Figure 19, however, uses the suboptimal Hilbert curve ordering. In this figure, it can be observed that the data is grouped into different clusters. Because the Hilbert curve traverses the same path in each dataset we know the approximate location of each cluster. That is the cube template used for the Hilbert ordering is split into 8 smaller sections. This aspect could easily be represented in the one dimensional portrayal of the data by dividing it into 8 equal sections The Hilbert curve traverses each section completely before moving on to the next one. This is important to note in the fact that each of these clusters is located in a certain section of the brain. The exact section, however, depends on where the Hilbert ordering began.
One way to use these clusters of data is by looking at the same clusters of data in different datasets. Because the Hilbert ordering is the same throughout the datasets, the multiple clusters would essentially align through all the datasets. This is helpful in analyzing data visually by observing whether there are any significant characteristics in each cluster, and observe how these characteristics are presented across the multiple datasets. This is a way in which we can essentially identify patterns across datasets from similar groups. The activation patterns in the different clusters of the data can provide significant clues in order to classify the data. 32
Figure 18: AQI dataset in linear form, both including all values and binned
Figure 19: AQI dataset using Hilbert curve ordering, both including all values and binned 33
4.2 Classification
Classification in Matlab was used in order to analyze the data obtained from the Hilbert curve ordering. The goal of classification is to analyze known data and its characteristics in order to find a defining similarity. Essentially, we use these similarities in order to try to determine to which group another similar dataset would belong to. There are many classification algorithms which can be applied already built into Matlab, therefore we decided on the use of the fitcdiscr (fit discriminant analysis) classifier. This classifying algorithm uses current known datasets and their characteristics in order to create a classifying object. We can then use this object in order to classify similar sets of data.
Using fitcdicr is very simple. The function only requires two input values, the predictor values and the classification values. The predictor values are the actual data which will be used to train the classifier. This is the data to which other datasets will be compared and assigned based on similar values. The classification values are the groups to which the predictor data belongs to. In our case, we used half of our datasets as the predictor values. Since we already know to which group each of the datasets belongs to we use the groups, control or patient, as the classification values, note that we only have two classification values. The classifier can classify the data using one of six classifying types:
linear
quadratic
diaglinear
diagquadratic
pseudolinear 34
pseudoquadratic.
After multiple tests, we determined that the most optimal classification types for our type of data would be diaglinear, pseudolinear, and pseudoquadratic. These types returned the best classification when tested on the data.
In order to obtain vast results and compare how the data used as the training group affected classification, we divided the datasets into two equal sets. Datasets 1-8 were used as the training group for one round of tests while groups 9-16 were used for the second round. We also tested the groups using two different bin sizes. By changing the size of the binned members we can determine whether certain data values may have a more specific effect on classification or if the dataset overall could serve as the classifying factor.
The following tables include the results of our classification. Tables 3 and 4 show the results from classifying data from both linear and Hilbert ordering using a bin size of 100. Tables 5 and 6 show the same type of results with a bin size of 1000. Each Table is divided into the two different training groups. Each group was used as the training group for the classification of the other. As we can see from the tables, classification was most effective when using the pseudoquadratic type of classification. In several instances, 6 of the 8 groups were classified correctly, that is a correct classification rate of about 75%. This supports that the datasets belonging to the same group, do indeed, have qualities that can be used in order to classify similar data as members of the same group.
The importance of these results is to show that members of the same groups actually share significant similarities in the activation maps. That is they are producing activation signals that help identify each one as a part of a certain group. 35
Table 3: Using Linear ordering with a bin size of 100
GroupCAVCAVCAVAQI_C0000100CHD_P1111111CIF_C0111111CML_C0111111CTT_P1000011CVX_P1111111DAY_P1111111DFJ_C0100000Correct455466DSN_P1001111DUQ_P1111111DUT_P1111111EWR_C0111111FHE_C0111111FLK_P1001111GWW_C0111111HCX_C0111111Correct224444Diaglinear:Pseudolinear:Pseudoquadratic:36
Table 4: Using Hilbert curve ordering with a bin size of 100
GroupCAVCAVCAVAQI_C0110000CHD_P1111111CIF_C0111111CML_C0111111CTT_P1000011CVX_P1001111DAY_P1111111DFJ_C0001000Correct334566DSN_P1000011DUQ_P1011111DUT_P1111111EWR_C0111111FHE_C0111111FLK_P1111111GWW_C0001011HCX_C0110011Correct344544DiaglinearPseudolinearPseudoquadratic37
Table 5: Linear ordering using a bin size of 1000
GroupCAVCAVCAVAQI_C0000000CHD_P1111111CIF_C0111111CML_C0111111CTT_P1000011CVX_P1111111DAY_P1111111DFJ_C0101000Correct454566DSN_P1000011DUQ_P1111111DUT_P1111111EWR_C0111111FHE_C0000001FLK_P1000011GWW_C0111100HCX_C0000011Correct444465DiaglinearPseudolinearPseudoquadratic38
Table 6: Hilbert curve ordering using bin size of 1000
GroupCAVCAVCAVAQI_C0111100CHD_P1111101CIF_C0111101CML_C0111101CTT_P1000011CVX_P1001101DAY_P1111101DFJ_C0000000Correct334456DSN_P1000011DUQ_P1110011DUT_P1111111EWR_C0111111FHE_C0001100FLK_P1111111GWW_C0001011HCX_C0000011Correct663455DiaglinearPseudolinearPseudoquadratic39
4.3 Conclusions
The activation brain mappings used in this project display areas affected by certain stimulations. Activation maps show signals of the brain due to certain stimulations producing different density signals. The activation values in this case helped identify similarities among members of the same control groups. In our case, we have found that the datasets belonging to the same group do in fact share similar enough characteristics that helps the solver identify them. A way to improve this classification would be to obtain more datasets from different groups, then, use these datasets in the classification process.
The inclusion of the Hilbert space filling curve will hopefully show that there is an optimal solution in brain mapping that is to be found. Classification of the data is essential in helping to determine whether the grey matter data in this project can help identify patterns in the data from similar groups. By finding the optimal solution we can improve classification by knowing exactly at which aspect of the data we are looking at in a single dimension.
4.4 In the future
The next step in finding our optimal solution would be to essentially solve the original ILP problem instance. This was found difficult to do because of the lack of computing resources available. There are powerful solvers out there, such as IBM’s CLPEX, that we could attempt to use, however, they require licenses which can run into the thousands of dollars. As of now, classifying data using suboptimal methods such as the Hilbert curve at least gives us an idea into how mapping the brain matter into a 1D can be essential in detecting abnormalities. 40
Bibliography
[1] “Computational Neuroscience: Constructing fMRI connectivity networks: A
whole brain functional parcellation method for node definition.” Journal of neuroscience methods 228 (2014) : 86.
[2] Coutanche, Mark N., Thompson-Schill, Sharon L., Schultz, Robert T."Multi
voxel pattern analysis of fMRI data predicts clinical symptom severity." NeuroImage 57.1 (2011):113.
[3] Frisoni, G B. Zorsan, A. Sabattoli, F. Beltramello, A. Soininen, H. Laakso, M P. “Detection
of grey matter loss in mild Alzheimer’s disease with vosel based morphometry”. J Neurol Neurosurg Psychiatry 2002;73:657–64
[4] Kontos, D., Megalooikonomou, V., Ghubade, N., Faloutsos, C.. (2003).
“Detecting Discriminative Functional MRI Activation Patterns Using Space Filling Curves”. In proceedings of the 25th Annual Conference of the IEEE EMBS. September 17-21, 2003. Cancun, Mexico.
[5] LP_solve, https://lpsolve.sourceforge.net. Accessed 2014.
[6] Sakoglu, Unal., Arslan, Abdullah. (2014). “In Search of Optimal Space-Filling
Curves for 3-D to 1-D Mapping: Application to 3-D Brain MRI Data”. In Proc. of the 6th International Conference on Biometrics and Computational Biology (BICoB). March 24-26, 2014 Las Vegas, Nevada