Uniformly select multi-shell subsets from a multi-shell scheme
This is a demo to uniformly separate multiple subsets from the multi-shell scheme in Human Connectome Project (HCP) using Mixed Integer Linear Programming (MILP)
Reference:
1. "Single- and Multiple-Shell Uniform Sampling Schemes for Diffusion MRI Using Spherical Codes", Jian Cheng, Dinggang Shen, Pew-Thian Yap, Peter J. Basser, IEEE Transactions on Medical Imaging, 2017.
2. "Designing Single- and Multiple-Shell Sampling Schemes for Diffusion MRI Using Spherical Code", Jian Cheng, Dinggang Shen, Pew-Thian Yap, MICCAI 2014.
Copyright (c) 2013, Jian Cheng (jian.cheng.1983@gmail.com)
Contents
read HCP scheme
grad_hcp_all{1}=ReadDirections('/home/jcheng/mywork/output/sampling/HCP_Q3/grad_b1000.txt');
grad_hcp_all{2}=ReadDirections('/home/jcheng/mywork/output/sampling/HCP_Q3/grad_b2000.txt');
grad_hcp_all{3}=ReadDirections('/home/jcheng/mywork/output/sampling/HCP_Q3/grad_b3000.txt');
VisualizeMultiShellScheme(grad_hcp_all{1}, grad_hcp_all{2}, grad_hcp_all{3});
title(['Original 90x3 HCP scheme']);
extract 30x3 samples from 90x3 scheme
clear params grbParams params.numSamples = [30, 30, 30]; % set a lower bound, it can be 0 or covering radius from an initial guess params.lbCost = [0.33,0.33,0.33,0.15]; % params.lbCost = [0,0,0,0]; params.w=0.5; % grb parameters % MIPFocus 1 grbParams.MIPFocus=1; grbParams.TimeLimit=600; grbParams.OutputFlag=true; % grbParams.DisplayInterval=1; % run [grad,grb, indexMatrix] = OptimalSamplingMultiSubsetsFromDifferentSets(grad_hcp_all, params, grbParams); % the selected indices for each shell index_1=find(indexMatrix{1}==1); index_2=find(indexMatrix{2}==1); index_3=find(indexMatrix{3}==1);
MIPFocus = 1
Academic license - for non-commercial use only
Optimize a model with 2649 rows, 274 columns and 8208 nonzeros
Variable types: 4 continuous, 270 integer (270 binary)
Coefficient statistics:
Matrix range [1e+00, 2e+00]
Objective range [2e-01, 5e-01]
Bounds range [1e-01, 1e+00]
RHS range [3e+00, 3e+01]
Presolve time: 0.01s
Presolved: 2649 rows, 274 columns, 8208 nonzeros
Variable types: 4 continuous, 270 integer (270 binary)
Presolved: 274 rows, 2923 columns, 8482 nonzeros
Presolve removed 274 rows and 2923 columns
Root relaxation: objective 3.882960e-01, 680 iterations, 0.03 seconds
Nodes | Current Node | Objective Bounds | Work
Expl Unexpl | Obj Depth IntInf | Incumbent BestBd Gap | It/Node Time
0 0 0.38830 0 178 - 0.38830 - - 0s
0 0 0.38830 0 184 - 0.38830 - - 0s
0 0 0.38830 0 188 - 0.38830 - - 0s
0 0 0.38830 0 191 - 0.38830 - - 0s
0 0 0.38830 0 191 - 0.38830 - - 0s
0 0 0.38830 0 193 - 0.38830 - - 0s
0 0 0.38830 0 191 - 0.38830 - - 0s
0 0 0.38830 0 200 - 0.38830 - - 0s
0 0 0.38830 0 189 - 0.38830 - - 0s
H 0 0 0.2400000 0.38830 61.8% - 0s
0 0 0.38830 0 183 0.24000 0.38830 61.8% - 0s
H 0 0 0.2428519 0.38830 59.9% - 1s
0 2 0.38830 0 180 0.24285 0.38830 59.9% - 1s
H 176 176 0.2440524 0.38830 59.1% 226 4s
H 230 223 0.2448611 0.38830 58.6% 193 4s
H 234 227 0.2448675 0.38830 58.6% 192 4s
241 224 0.33673 44 179 0.24487 0.38830 58.6% 191 5s
H 267 224 0.2453033 0.38830 58.3% 177 5s
468 417 0.30596 53 148 0.24530 0.38830 58.3% 151 10s
591 482 0.37847 14 247 0.24530 0.38830 58.3% 134 16s
599 487 0.37361 16 237 0.24530 0.38830 58.3% 134 21s
669 531 0.37371 19 237 0.24530 0.38830 58.3% 131 26s
762 634 0.35826 26 208 0.24530 0.38830 58.3% 125 30s
905 731 0.31388 35 173 0.24530 0.38830 58.3% 117 35s
984 781 0.37895 9 248 0.24530 0.38830 58.3% 115 41s
1099 869 0.35599 24 183 0.24530 0.38830 58.3% 110 45s
1142 902 0.35316 17 270 0.24530 0.36921 50.5% 115 50s
1163 917 0.26780 69 270 0.24530 0.36881 50.4% 116 59s
1164 921 0.36801 56 269 0.24530 0.36881 50.3% 119 60s
1708 1203 0.35933 61 252 0.24530 0.36383 48.3% 114 66s
1833 1234 0.34061 67 228 0.24530 0.36279 47.9% 112 71s
1845 1237 0.34041 68 225 0.24530 0.36279 47.9% 112 77s
2287 1394 0.25886 88 68 0.24530 0.36106 47.2% 108 86s
2698 1481 infeasible 78 0.24530 0.35999 46.8% 105 94s
2749 1481 0.34747 66 235 0.24530 0.35896 46.3% 104 95s
3918 1901 0.33343 74 205 0.24530 0.35689 45.5% 99.2 100s
4393 2167 0.33738 71 217 0.24530 0.35608 45.2% 96.6 111s
4650 2293 0.32520 70 211 0.24530 0.35545 44.9% 96.0 115s
4896 2411 0.31401 73 165 0.24530 0.35514 44.8% 94.4 123s
5350 2680 0.32393 73 195 0.24530 0.35468 44.6% 94.3 125s
5590 2821 0.32411 73 193 0.24530 0.35468 44.6% 93.7 133s
6014 3069 infeasible 75 0.24530 0.35413 44.4% 93.0 135s
6285 3214 0.34174 69 228 0.24530 0.35407 44.3% 92.5 143s
6578 3405 0.34298 70 236 0.24530 0.35380 44.2% 91.1 152s
6864 3562 0.34825 64 232 0.24530 0.35363 44.2% 90.4 156s
7168 3708 0.32945 65 207 0.24530 0.35353 44.1% 90.2 166s
7479 3924 0.34595 64 223 0.24530 0.35328 44.0% 89.9 175s
7767 4071 0.34268 68 224 0.24530 0.35311 43.9% 89.5 184s
8096 4245 0.34902 65 235 0.24530 0.35302 43.9% 89.3 187s
8418 4439 0.34654 66 226 0.24530 0.35276 43.8% 88.9 191s
8593 4546 0.34238 67 219 0.24530 0.35265 43.8% 88.7 198s
9445 5061 0.33150 75 209 0.24530 0.35225 43.6% 87.5 209s
9955 5287 0.32819 76 184 0.24530 0.35200 43.5% 86.9 216s
10409 5554 0.34429 65 223 0.24530 0.35183 43.4% 86.7 227s
11149 5937 0.34245 68 218 0.24530 0.35161 43.3% 85.5 232s
11859 6326 0.28950 84 126 0.24530 0.35132 43.2% 84.7 240s
12593 6832 0.31327 78 180 0.24530 0.35108 43.1% 83.9 248s
13277 7168 infeasible 75 0.24530 0.35102 43.1% 83.1 250s
14188 7700 0.31776 74 191 0.24530 0.35079 43.0% 82.5 258s
15043 8181 0.33841 65 221 0.24530 0.35029 42.8% 82.0 270s
17081 9368 0.34399 63 221 0.24530 0.34977 42.6% 80.5 279s
17994 9857 0.33757 66 214 0.24530 0.34958 42.5% 80.2 283s
18922 10385 infeasible 87 0.24530 0.34931 42.4% 80.1 294s
20000 11022 0.31568 75 188 0.24530 0.34912 42.3% 79.4 304s
20701 11383 0.33299 66 216 0.24530 0.34907 42.3% 79.1 315s
21868 12033 0.34407 66 214 0.24530 0.34876 42.2% 78.8 320s
23005 12677 0.34139 67 213 0.24530 0.34848 42.1% 78.4 329s
24216 13333 0.32214 72 186 0.24530 0.34818 41.9% 78.0 343s
H24451 13473 0.2453033 0.34816 41.9% 78.0 343s
24779 13578 0.30193 76 135 0.24530 0.34810 41.9% 77.9 352s
25878 14193 0.33503 72 210 0.24530 0.34796 41.8% 77.4 361s
26805 14679 0.32068 67 210 0.24530 0.34784 41.8% 77.5 368s
27950 15259 0.31320 82 179 0.24530 0.34759 41.7% 77.3 376s
28927 15785 0.34027 65 215 0.24530 0.34741 41.6% 77.0 388s
30478 16604 0.28296 76 129 0.24530 0.34720 41.5% 76.3 394s
31293 16991 0.33997 65 231 0.24530 0.34704 41.5% 76.2 404s
32408 17578 0.33734 68 212 0.24530 0.34686 41.4% 76.1 409s
33753 18309 0.28921 84 125 0.24530 0.34672 41.3% 75.7 415s
34651 18782 0.34376 68 227 0.24530 0.34660 41.3% 75.6 426s
35951 19452 0.34030 71 208 0.24530 0.34648 41.2% 75.3 434s
37304 20199 0.32439 76 199 0.24530 0.34629 41.2% 74.9 437s
38428 20771 0.31272 74 182 0.24530 0.34610 41.1% 74.9 443s
39669 21446 0.32893 70 205 0.24530 0.34604 41.1% 74.7 452s
40774 22081 0.32437 71 214 0.24530 0.34585 41.0% 74.5 458s
41656 22544 0.32818 75 207 0.24530 0.34576 41.0% 74.3 466s
42748 23174 0.27776 79 120 0.24530 0.34564 40.9% 74.2 471s
43797 23693 0.33174 69 226 0.24530 0.34552 40.9% 74.1 476s
44667 24099 0.34206 66 228 0.24530 0.34541 40.8% 74.1 480s
45680 24622 0.32981 68 214 0.24530 0.34535 40.8% 74.0 492s
46358 24908 0.31552 75 173 0.24530 0.34527 40.8% 73.9 502s
47652 25569 0.34197 69 229 0.24530 0.34512 40.7% 73.7 505s
48786 26128 0.30992 70 179 0.24530 0.34505 40.7% 73.5 510s
49739 26630 0.32950 69 199 0.24530 0.34495 40.6% 73.4 518s
50442 26943 0.33310 69 206 0.24530 0.34484 40.6% 73.4 526s
52435 27971 0.31412 78 172 0.24530 0.34472 40.5% 73.2 540s
54821 29298 0.33110 74 210 0.24530 0.34452 40.4% 72.9 545s
56060 29826 0.32992 74 197 0.24530 0.34440 40.4% 72.7 553s
57000 30303 0.34154 69 227 0.24530 0.34434 40.4% 72.7 560s
57559 30521 0.30759 80 178 0.24530 0.34430 40.4% 72.7 566s
58671 31032 0.33593 73 207 0.24530 0.34416 40.3% 72.5 570s
60638 31992 0.33873 71 221 0.24530 0.34399 40.2% 72.4 578s
61817 32621 0.32681 72 196 0.24530 0.34389 40.2% 72.2 582s
63024 33193 0.32974 70 223 0.24530 0.34380 40.2% 72.0 586s
63551 33385 0.32449 76 194 0.24530 0.34376 40.1% 72.0 595s
64676 33956 0.33517 67 212 0.24530 0.34366 40.1% 71.8 600s
Cutting planes:
Clique: 94
MIR: 19
StrongCG: 154
Flow cover: 159
Zero half: 1
Explored 65647 nodes (4712674 simplex iterations) in 600.02 seconds
Thread count was 8 (of 8 available processors)
Solution count 7: 0.245303 0.245303 0.244868 ... 0.24
Time limit reached
Best objective 2.453033102693e-01, best bound 3.436474180708e-01, gap 40.0908%
visualize the result
VisualizeMultiShellScheme(grad{1},grad{2},grad{3});
title({'Estimated scheme.', ['Combined covering radius = ', num2str(CoveringRadius([grad{1};grad{2};grad{3}])*180/pi), ' degree']});
