0001 function head_join = vb_job_head_3shell(proj_root, parm)
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020
0021
0022
0023
0024
0025
0026
0027
0028
0029
0030
0031
0032
0033
0034
0035
0036
0037
0038
0039
0040
0041
0042
0043
0044
0045
0046
0047
0048
0049
0050
0051
0052
0053
0054
0055
0056
0057
0058
0059
0060
0061
0062
0063
0064
0065
0066
0067
0068
0069
0070
0071
0072
0073
0074
0075
0076
0077
0078
0079
0080
0081
0082
0083
0084
0085
0086
0087
0088
0089
0090
0091
0092
0093
0094 proj_root = vb_rm_trailing_slash(proj_root);
0095
0096
0097 head_parm = vb_set_head_parm;
0098
0099
0100 if ~isfield(parm,'Nsurf'),
0101 Nsurf = head_parm.Nsurf;
0102 else
0103 Nsurf = parm.Nsurf;
0104 end;
0105
0106
0107 if ~isfield(parm,'Nvertex'),
0108 switch Nsurf
0109 case 1
0110 Nvertex = [6000];
0111 case 3
0112 Nvertex = head_parm.Nvertex;
0113 end
0114 else
0115 Nvertex = parm.Nvertex;
0116 end;
0117 if length(Nvertex)==1, Nvertex = repmat(Nvertex,[1,Nsurf]); end;
0118
0119
0120 if ~isfield(parm,'Radius_fs'), parm.Radius_fs = head_parm.Radius_fs ;end;
0121 if ~isfield(parm,'Radius_csf'), parm.Radius_csf = head_parm.Radius_csf;end;
0122 if ~isfield(parm,'Radius_scalp'),parm.Radius_scalp = head_parm.Radius_scalp;end
0123 if ~isfield(parm,'Radius_skull'),parm.Radius_skull = head_parm.Radius_skull;end
0124 if ~isfield(parm,'Skull_ratio'), parm.Skull_ratio = head_parm.Skull_ratio ;end
0125
0126 if ~isfield(parm,'Radius_final'),parm.Radius_final = head_parm.Radius_final;end
0127 if ~isfield(parm,'Radius'), parm.Radius = head_parm.Radius ;end;
0128 if ~isfield(parm,'Radius_fill'), parm.Radius_fill = head_parm.Radius_fill;end;
0129 if ~isfield(parm,'Radius_gray'), parm.Radius_gray = head_parm.Radius_gray;end;
0130 if ~isfield(parm,'Glevel'), parm.Glevel = head_parm.Glevel ;end;
0131
0132 if ~isfield(parm,'Nloop'), parm.Nloop = head_parm.Nloop ;end;
0133 if ~isfield(parm,'Nlast_csf'), parm.Nlast_csf = head_parm.Nlast_csf ;end;
0134 if ~isfield(parm,'Nlast_scalp'), parm.Nlast_scalp = head_parm.Nlast_scalp;end;
0135 if ~isfield(parm,'Nlast_skull'), parm.Nlast_skull = head_parm.Nlast_skull;end;
0136
0137 if ~isfield(parm,'Scalp_mode'), parm.Scalp_mode = head_parm.Scalp_mode;end;
0138 if ~isfield(parm,'vstep'), parm.vstep = head_parm.vstep; end;
0139
0140 step = parm.vstep;
0141
0142
0143 if ~isfield(parm,'plot_mode'),
0144 plot_mode = 0;
0145 else
0146 plot_mode = parm.plot_mode;
0147 end;
0148
0149
0150 analyzefile = parm.analyze_file;
0151 if isfield(parm, 'brain_file')
0152 brain_file = [proj_root '/' parm.brain_file];
0153 brain_parm = parm;
0154 brain_parm.brain_file = brain_file;
0155 else
0156 brain_file = [];
0157 end
0158
0159
0160 [headfile, basename] = change_name(parm.head_file);
0161
0162 headmodel{1} = change_name([basename '_CSF' ], Nvertex(1));
0163 headfile = [proj_root '/' headfile];
0164 headsave{1} = [proj_root '/' headmodel{1}];
0165
0166 if Nsurf == 3
0167
0168 headmodel{2} = change_name([basename '_Skull'], Nvertex(2));
0169 headmodel{3} = change_name([basename '_Scalp'], Nvertex(3));
0170 head_join = change_name([basename '_3shell'], sum(Nvertex));
0171 head_join = [proj_root '/' head_join];
0172 headsave{2} = [proj_root '/' headmodel{2}];
0173 headsave{3} = [proj_root '/' headmodel{3}];
0174 else
0175
0176 head_join = headsave{1};
0177 end
0178
0179
0180 [Vdim, Vsize] = analyze_hdr_read(analyzefile);
0181
0182 [V,F,Vs,Fs] = load_face(parm, Vdim, Vsize);
0183
0184
0185 if plot_mode > 0
0186 debug_plot(V, F, analyzefile);
0187 end
0188
0189
0190
0191
0192 parm.Nvertex = Nvertex(1);
0193 parm.Nlast = parm.Nlast_csf;
0194
0195 DW = 20;
0196
0197
0198
0199 if isfield(parm,'freesurf_file') && ~isempty(parm.freesurf_file)
0200
0201
0202 if ~exist(headfile,'file')
0203 [Vhead,Fhead] = vb_fs_load_surface(parm.freesurf_file);
0204 [Fhead,Vhead,XXhead] = vb_out_normal_surf(Fhead,Vhead);
0205 Vhead = Vhead * 0.001;
0206
0207 disp(['--- Save CSF model file: ']);
0208 disp([headfile]);
0209 vb_fsave(headfile,'Vhead','Fhead','XXhead');
0210 else
0211 load(headfile,'Vhead','Fhead','XXhead');
0212 end
0213
0214
0215 if plot_mode > 0
0216 debug_plot(Vhead, Fhead, analyzefile);
0217 end
0218
0219 if ~isempty(brain_file)
0220
0221 B = vb_make_head_surf_model(brain_parm);
0222 Dim = size(B);
0223
0224 Bcsf = zeros(Dim + 2*DW);
0225 Bcsf((1:Dim(1))+DW, (1:Dim(2))+DW, (1:Dim(3))+DW) = B;
0226 Dim = Dim + 2*DW;
0227 else
0228 Dim = ceil((Vdim .* Vsize)/step) + 2*DW;
0229
0230 end
0231
0232
0233 Vhead = vb_spm_right_to_analyze_right_mm(Vhead, Vdim, Vsize);
0234 Vhead = Vhead + step*DW;
0235
0236
0237 B = vb_surf_to_filled_mask(Vhead, Fhead, Dim, step);
0238
0239
0240 B = vb_morphology_operation(B, parm.Radius_fs, step);
0241
0242
0243 if ~isempty(brain_file)
0244
0245
0246
0247 if vb_matlab_version > 6,
0248 B = int8(B) + int8(Bcsf);
0249 else
0250 B = double(B) + double(Bcsf);
0251 end
0252 end
0253
0254
0255 if ~isempty(Vs)
0256 Vs = Vs + step*DW;
0257 Bface = vb_face_to_mask(Vs, Fs, Dim, step);
0258 Radius = [- parm.Radius_scalp, - parm.Radius_skull];
0259 Bface = vb_morphology_operation(Bface, Radius , step);
0260
0261
0262 if vb_matlab_version > 6,
0263 B = int8(B) .* int8(Bface);
0264 else
0265 B = double(B) .* double(Bface);
0266 end
0267 end
0268
0269
0270
0271 [Vhead, Fhead, XXhead, B] = vb_mask_to_surf_expand(B, parm);
0272 Vhead = Vhead - step*DW;
0273
0274 Vhead = vb_analyze_right_mm_to_spm_right(Vhead, Vdim, Vsize);
0275
0276 disp(['--- Save CSF model file: ']);
0277 disp([headsave{1}]);
0278 vb_fsave(headsave{1},'Vhead','Fhead','XXhead');
0279
0280 if Nsurf==1
0281 vb_plot_slice_head(headsave{1}, analyzefile, brain_file);
0282 end
0283 elseif isfield(parm,'curry_file') && ~isempty(parm.curry_file)
0284
0285
0286 if ~exist(headfile,'file')
0287 vb_job_curry2vbmeg(parm.curry_file, headfile, analyzefile);
0288 end
0289
0290
0291 if Nsurf==3, parm.plot_mode = 0; end;
0292
0293 vb_job_smooth_head(headfile, headsave{1}, parm);
0294
0295 elseif isfield(parm,'gray_file') && ~isempty(parm.gray_file)
0296
0297
0298 [B] = vb_make_head_surf_model(brain_parm);
0299
0300
0301 Dim = size(B);
0302
0303
0304 if ~isempty(Vs)
0305 Bface = vb_face_to_mask(Vs, Fs, Dim, step);
0306 Radius = [- parm.Radius_scalp, - parm.Radius_skull];
0307 Bface = vb_morphology_operation(Bface, Radius , step);
0308
0309
0310 if vb_matlab_version > 6,
0311 B = int8(B) .* int8(Bface);
0312 else
0313 B = double(B) .* double(Bface);
0314 end
0315 end
0316
0317
0318
0319 [Vhead, Fhead, XXhead, B] = vb_mask_to_surf_expand(B, parm);
0320 Vhead = vb_analyze_right_mm_to_spm_right(Vhead, Vdim, Vsize);
0321
0322 if isempty(XXhead), return; end;
0323
0324 disp(['--- Save CSF model file: ']);
0325 disp([headsave{1}]);
0326 vb_fsave(headsave{1},'Vhead','Fhead','XXhead');
0327
0328 if Nsurf==1,
0329 vb_plot_slice_head(headsave{1}, analyzefile, brain_file);
0330 end
0331 else
0332 error('No CSF input file');
0333 end
0334
0335 if Nsurf==1,
0336
0337 project_file_saving(parm);
0338 return;
0339 end;
0340
0341 if isempty(V), error('Scalp file is needed for 3 shell model'); end
0342
0343
0344 load(headsave{1},'Vhead','Fhead')
0345
0346
0347 if plot_mode > 0
0348 debug_plot(Vhead, Fhead, analyzefile);
0349 if plot_mode == 1, return; end;
0350 end
0351
0352
0353
0354
0355 [Rmax,Rmin] = vb_distance_brain_face(Vhead,V);
0356
0357 fprintf('Max distance from face to CSF = %5.1f\n',Rmax)
0358 fprintf('Min distance from face to CSF = %5.1f\n',Rmin)
0359
0360 Rstep = round(Rmax/2);
0361 Radius = [Rstep ceil(Rmax + 2 - Rstep)];
0362 fprintf('Expansion Radius from CSF to Scalp = %3.0f\n',sum(Radius))
0363
0364
0365
0366 [B ,Bcsf, Bface] = ...
0367 vb_intersect_face_csf(Vhead,Fhead,V,F,Radius,step,Vdim,Vsize,DW);
0368
0369
0370 if plot_mode > 1
0371 zindx = [40:20:200];
0372 zindx = ceil(zindx/step);
0373 vb_plot_slice( B, [], zindx, 1);
0374 vb_plot_slice( Bcsf, [], zindx, 1);
0375 vb_plot_slice( Bface, [], zindx, 1);
0376 if plot_mode == 2, return; end;
0377 end
0378
0379 parm.Nlast = parm.Nlast_scalp;
0380
0381 if parm.Scalp_mode == 0
0382
0383 parm.Nvertex = Nvertex(3);
0384 parm.Radius = parm.Radius_final;
0385
0386 [Vhead, Fhead, XXhead, B] = vb_mask_to_surf_expand(B, parm);
0387
0388 Vhead = Vhead - step*DW;
0389 Vhead = vb_analyze_right_mm_to_spm_right(Vhead,Vdim,Vsize);
0390
0391 disp(['--- Save Scalp model file: ']);
0392 disp([headsave{3}]);
0393 vb_fsave(headsave{3},'Vhead','Fhead','XXhead');
0394 else
0395
0396 parm.Radius = parm.Radius_final;
0397 parm.Nvertex = Nvertex(3);
0398
0399 [Vhead, Fhead, XXhead, Bface] = vb_mask_to_surf_expand(Bface, parm);
0400
0401 Vhead = Vhead - step*DW;
0402 Vhead = vb_analyze_right_mm_to_spm_right(Vhead,Vdim,Vsize);
0403
0404 disp(['--- Save Face model file: ']);
0405 disp([headsave{3}]);
0406 vb_fsave(headsave{3},'Vhead','Fhead','XXhead');
0407
0408 switch parm.Scalp_mode
0409 case 1
0410
0411 parm.Radius = parm.Radius_final;
0412 B = vb_morphology_operation(B, parm.Radius, step);
0413 case 2
0414
0415 B = Bface;
0416 end
0417 end
0418
0419 clear Bface
0420
0421
0422 if plot_mode > 2
0423 debug_plot(Vhead, Fhead, analyzefile);
0424 if plot_mode == 3, return; end;
0425 end
0426
0427
0428
0429
0430 if isfield(parm,'fs_skull_file') && ~isempty(parm.fs_skull_file)
0431
0432
0433
0434
0435 load(headsave{1},'Vhead','Fhead')
0436 [V,F] = vb_fs_load_surface(parm.fs_skull_file);
0437 V = V * 0.001;
0438
0439
0440 Radius(2) = Radius(2) - parm.Radius_scalp;
0441
0442
0443 [B ,Btmp, Bskl] = ...
0444 vb_intersect_face_csf(Vhead,Fhead,V,F,Radius,step,Vdim,Vsize,DW);
0445
0446 clear Btmp
0447 switch parm.Scalp_mode
0448 case {0,1}
0449
0450 parm.Radius = parm.Radius_final;
0451 B = vb_morphology_operation(B, parm.Radius, step);
0452 case 2
0453
0454 B = Bskl;
0455 end
0456 else
0457
0458 fprintf('Max width of Scalp = %5.1f\n',parm.Radius_scalp)
0459
0460 B = vb_morphology_operation(B, - parm.Radius_scalp , step);
0461
0462 end
0463
0464 Rskull = Rmin * parm.Skull_ratio;
0465 fprintf('Min width of Skull = %5.1f\n',Rskull)
0466
0467 Bcsf = vb_morphology_operation(Bcsf, Rskull, step);
0468
0469
0470
0471 if vb_matlab_version > 6,
0472 B = int8(B) + int8(Bcsf);
0473 ix = find( B(:) > 0 );
0474 B = zeros(size(B),'int8');
0475 B(ix) = 1;
0476 else
0477 B = double(B) + double(Bcsf);
0478 ix = find( B(:) > 0 );
0479 B = zeros(size(B));
0480 B(ix) = 1;
0481 end
0482
0483
0484
0485
0486 parm.Radius = parm.Radius_final;
0487 parm.Nvertex = Nvertex(2);
0488 parm.Nlast = parm.Nlast_skull;
0489
0490 [Vhead, Fhead, XXhead, B] = vb_mask_to_surf_expand(B, parm);
0491
0492 Vhead = Vhead - step*DW;
0493 Vhead = vb_analyze_right_mm_to_spm_right(Vhead,Vdim,Vsize);
0494
0495 disp(['--- Save Skull model file: ']);
0496 disp([headsave{2}]);
0497 vb_fsave(headsave{2},'Vhead','Fhead','XXhead');
0498
0499
0500 vb_head_join_files(headmodel, proj_root, head_join);
0501
0502
0503
0504
0505
0506 vb_plot_slice_head(headsave, analyzefile, brain_file);
0507
0508
0509
0510
0511 load(headsave{3},'Vhead','Fhead');
0512
0513 if exist('surf_face','var') && isfield(surf_face,'V_reduce')
0514 V = surf_face.V_reduce;
0515 F = surf_face.F_reduce;
0516 end
0517
0518 figure
0519 subplot(1,2,1)
0520 vb_plot_surf(V,F,[],[],1,1);
0521 hold on
0522 vb_plot_surf(Vhead,Fhead,'none','k');
0523 alpha 0.8;
0524
0525 title('Scalp surface on face')
0526
0527 view([ -210 30]);
0528
0529 subplot(1,2,2)
0530 vb_plot_surf(Vhead,Fhead,[],[],1,1);
0531
0532 title('Scalp surface')
0533
0534 view([ -210 30]);
0535
0536 project_file_saving(parm);
0537
0538 return
0539
0540 function project_file_saving(parm)
0541
0542
0543
0544 proj_file = get_project_filename;
0545 if isempty(proj_file)
0546 return;
0547 end
0548 project_file_mgr('load', proj_file);
0549 project_file_mgr('add', 'head_parm', parm);
0550
0551
0552 return
0553
0554 function [fnew, basename] = change_name(fold, N);
0555
0556
0557
0558
0559 fext = '.head.mat';
0560 ix_ext = findstr(fold, fext);
0561
0562 if ~isempty(ix_ext),
0563 basename = fold(1:ix_ext-1);
0564 else
0565 basename = fold;
0566 end
0567
0568 if nargin == 2
0569 fnew = sprintf('%s_%d%s', basename, N, fext);
0570 else
0571 fnew = sprintf('%s%s', basename, fext);
0572 end
0573
0574 return
0575
0576 function [V,F,Vs,Fs] = load_face(parm, Vdim, Vsize);
0577
0578
0579 if isfield(parm,'fs_scalp_file') && ~isempty(parm.fs_scalp_file)
0580 [V,F] = vb_fs_load_surface(parm.fs_scalp_file);
0581 V = V * 0.001;
0582 elseif isfield(parm,'face_file') && ~isempty(parm.face_file)
0583 load(parm.face_file,'surf_face');
0584 V = surf_face.V;
0585 F = surf_face.F;
0586 else
0587 V = []; F = [];
0588 Vs = []; Fs = [];
0589 return
0590 end
0591
0592
0593 [F,V] = vb_close_surf(F,V);
0594
0595
0596 Vs = vb_spm_right_to_analyze_right_mm(V, Vdim, Vsize);
0597 Fs = F;
0598
0599
0600 zmin = 2;
0601 Vs(:,3) = max(Vs(:,3),zmin);
0602
0603 return
0604
0605 function debug_plot(V, F, analyzefile)
0606 [B, Vdim, Vsize] = vb_load_analyze_to_right(analyzefile);
0607 Vana = vb_spm_right_to_analyze_right( V, Vdim, Vsize);
0608
0609 zindx = 0.1*[2:2:8];
0610 xindx = 0.1*[2:2:8];
0611
0612 Vmin = min(Vana);
0613 Vmax = max(Vana);
0614 zindx = ceil( zindx*(Vmax(3) - Vmin(3)) + Vmin(3));
0615 xindx = ceil( xindx*(Vmax(1) - Vmin(1)) + Vmin(1));
0616 zindx = zindx( zindx > 0 & zindx <= Vdim(3) );
0617 xindx = xindx( xindx > 0 & xindx <= Vdim(1) );
0618
0619 vb_plot_slice_surf(B, Vana, F, zindx, 'z', [2 2]);
0620 vb_plot_slice_surf(B, Vana, F, xindx, 'x', [2 2]);
0621
0622 return
0623
0624 ix = find( Vs(:,3) >= zmin);
0625 Fs = vb_patch_select2(ix,F,size(Vs,1));