diff --git a/brainCrop.m b/brainCrop.m index 4a856bd..24811c5 100644 --- a/brainCrop.m +++ b/brainCrop.m @@ -4,8 +4,8 @@ % (c) 2024 Gavin Hsu and Andrew Birnbaum function bbox = brainCrop(subj) - [dirname,baseFilename,ext] = fileparts(subj); - seg = [dirname,'\',baseFilename,'_masks',ext]; + [dirname,baseFilename] = fileparts(subj); + seg = [dirname filesep baseFilename,'_masks.nii']; data = load_untouch_nii(seg); wm = data.img == 1; thresh = [0,0,100]; diff --git a/meshByIso2mesh.m b/meshByIso2mesh.m index 83feb0a..614c50e 100644 --- a/meshByIso2mesh.m +++ b/meshByIso2mesh.m @@ -45,7 +45,7 @@ allMaskShow(data.img>0) = numOfTissue + 2; % sliceshow(allMask,[],[],[],'Tissue index','Segmentation. Click anywhere to navigate.') -sliceshow(allMaskShow,[],[],[],'Tissue index','Segmentation. Click anywhere to navigate.',[],mri2mni,[]) +sliceshow(allMaskShow,[],[],[],'Tissue index','Segmentation. Click anywhere to navigate.',[],mri2mni) drawnow % allMask = uint8(allMask); diff --git a/releasingNotes.txt b/releasingNotes.txt index 7b4cbe9..3a725b9 100644 --- a/releasingNotes.txt +++ b/releasingNotes.txt @@ -16,6 +16,8 @@ Added a function to set the headers of the MRI and segmentation masks so that th Added code to handle MRIs with nonzero qform in their header. +Upgraded the sliceshow() function for better visualization of slices in a 3D volume. + ROAST V3.0 (codename: Hell's Kitchen) -- 2019/09/02: diff --git a/reviewRes.m b/reviewRes.m index b7077fc..c6b67c0 100644 --- a/reviewRes.m +++ b/reviewRes.m @@ -71,7 +71,7 @@ function reviewRes(subj,simTag,tissue,fastRender,tarTag) % ROAST by Andrew Birnbaum % April 2024 % -% (c) May 2024 Gavin Hsu +% (c) May 2024, sliceshow improved by Gavin Hsu and Andrew Birnbaum fprintf('\n\n'); disp('=============================================================') @@ -245,7 +245,7 @@ function reviewRes(subj,simTag,tissue,fastRender,tarTag) if ~exist(subjRasRSPD,'file') error(['The subject MRI you provided ' subjRasRSPD ' does not exist. Check if you run through resampling or zero-padding if you tried to do that.']); else - data = load_untouch_nii(subjRasRSPD); sliceshow(data.img,[],'gray',[],[],'MRI: Click anywhere to navigate.',[],mri2mni,[]); drawnow + data = load_untouch_nii(subjRasRSPD); sliceshow(data.img,[],'gray',[],[],'MRI: Click anywhere to navigate.',[],mri2mni); drawnow end if ~isempty(optRoast.T2) %T2 specified @@ -253,7 +253,7 @@ function reviewRes(subj,simTag,tissue,fastRender,tarTag) error(['T2 file ' optRoast.T2 ' does not exist. You used that to run ROAST but maybe later deleted it.']); else data = load_untouch_nii(optRoast.T2); - sliceshow(data.img,[],'gray',[],[],'MRI: T2. Click anywhere to navigate.',[],mri2mni,[]); drawnow + sliceshow(data.img,[],'gray',[],[],'MRI: T2. Click anywhere to navigate.',[],mri2mni); drawnow end end else @@ -300,7 +300,7 @@ function reviewRes(subj,simTag,tissue,fastRender,tarTag) allMaskShow = masks.img; allMaskShow(gel.img>0) = numOfTissue + 1; allMaskShow(elec.img>0) = numOfTissue + 2; - sliceshow(allMaskShow,[],[],[],'Tissue index','Segmentation. Click anywhere to navigate.',[],mri2mni,[]) + sliceshow(allMaskShow,[],[],[],'Tissue index','Segmentation. Click anywhere to navigate.',[],mri2mni) drawnow else @@ -519,7 +519,13 @@ function reviewRes(subj,simTag,tissue,fastRender,tarTag) nan_mask(find(mask)) = 1; cm = colormap(jet(2^11)); cm = [1 1 1;cm]; -bbox = brainCrop(subjRasRSPD); +if strcmp(tissue,'white') | strcmp(tissue,'gray') | strcmp(tissue,'brain') + bbox = brainCrop(subjRasRSPDSeg); + pos = round(mean(bbox)); +else + bbox = []; + pos = []; +end if isRoast @@ -531,13 +537,13 @@ function reviewRes(subj,simTag,tissue,fastRender,tarTag) end figName = ['Voltage in Simulation: ' simTag]; - sliceshow(vol_all.*nan_mask,[],cm,[],'Voltage (mV)',[figName '. Click anywhere to navigate.'],[],mri2mni,bbox); drawnow + sliceshow(vol_all.*nan_mask,pos,cm,[],'Voltage (mV)',[figName '. Click anywhere to navigate.'],[],mri2mni,bbox); drawnow figName = ['Electric field in Simulation: ' simTag]; for i=1:size(ef_all,4), ef_all(:,:,:,i) = ef_all(:,:,:,i).*nan_mask; end ef_mag = ef_mag.*nan_mask; dataShowVal = ef_mag(~isnan(ef_mag(:))); - sliceshow(ef_mag,[],cm,[min(dataShowVal) prctile(dataShowVal,95)],'Electric field (V/m)',[figName '. Click anywhere to navigate.'],ef_all,mri2mni,bbox); drawnow + sliceshow(ef_mag,pos,cm,[min(dataShowVal) prctile(dataShowVal,95)],'Electric field (V/m)',[figName '. Click anywhere to navigate.'],ef_all,mri2mni,bbox); drawnow else diff --git a/sliceshow.m b/sliceshow.m index 5749c84..e9259ca 100644 --- a/sliceshow.m +++ b/sliceshow.m @@ -20,7 +20,7 @@ function sliceshow(img,pos,color,clim,label,figName,vecImg,mri2mni,bbox) % % vecImg is the 3D vectorial field, e.g. an electric field. % -% mri2mni is the mapping from MRI voxel space to the MNI space, so +% mri2mni is the mapping from MRI voxel space to the MNI space. % % bbox is the output of brainCrop(). Use this to display only the brain. % @@ -31,33 +31,14 @@ function sliceshow(img,pos,color,clim,label,figName,vecImg,mri2mni,bbox) % (c) August 2019, Yu (Andy) Huang % (c) 2024, Gavin Hsu and Andrew Birnbaum -if nargin >= 2 && ~isempty(bbox) - minR = bbox(1,1); - maxR = bbox(2,1); - minA = bbox(1,2); - maxA = bbox(2,2); - minS = bbox(1,3); - maxS = bbox(2,3); - img = img(minR:maxR, minA:maxA, minS:maxS); -else - minR = 0; - minA = 0; - minS = 0; -end - if nargin<1 || isempty(img) - error('At least give us a volume to display') + error('At least give us a volume to display'); else [Nx,Ny,Nz] = size(img); - mydata.img = img; end if nargin<2 || isempty(pos) - mydata.pos = round(size(img)/2); - mydata.pos_crop = mydata.pos + [minR,minA,minS]; -else - mydata.pos=pos-[minR,minA,minS]; - mydata.pos_crop = mydata.pos + [minR,minA,minS]; + pos = round(size(img)/2); end if nargin<3 || isempty(color) @@ -89,18 +70,12 @@ function sliceshow(img,pos,color,clim,label,figName,vecImg,mri2mni,bbox) end if nargin<7 || isempty(vecImg) - mydata.vecImg = []; + vecImg = []; else - if nargin >= 2 && ~isempty(bbox) - vecImg = vecImg(minR:maxR, minA:maxA, minS:maxS,:); - end temp = size(vecImg); if any(temp(1:3)~=[Nx,Ny,Nz]) || temp(4)~=3 error('Vector field does not have correct size.'); end - mydata.vecImg = vecImg; - [xi,yi,zi] = ndgrid(1:Nx,1:Ny,1:Nz); - mydata.xi = xi; mydata.yi = yi; mydata.zi = zi; end if nargin<8 || isempty(mri2mni) @@ -110,9 +85,33 @@ function sliceshow(img,pos,color,clim,label,figName,vecImg,mri2mni,bbox) error('Unrecognized format of the voxel-to-MNI mapping.'); end mydata.mri2mni = mri2mni; - mydata.mnipos = round(mydata.mri2mni*[mydata.pos_crop 1]'); end +if nargin<9 || isempty(bbox) + minR = 1; minA = 1; minS = 1; + maxR = Nx; maxA = Ny; maxS = Nz; +else + minR = bbox(1,1); + maxR = bbox(2,1); + minA = bbox(1,2); + maxA = bbox(2,2); + minS = bbox(1,3); + maxS = bbox(2,3); +end + +% adjustment for bbox +img = img(minR:maxR, minA:maxA, minS:maxS); mydata.img = img; +temp = pos-[minR-1,minA-1,minS-1]; +if any(temp<=0), error('Voxel selected falls outside of the bounding box.'); end +mydata.pos = temp; +mydata.pos_crop = mydata.pos + [minR-1,minA-1,minS-1]; +if ~isempty(vecImg), vecImg = vecImg(minR:maxR, minA:maxA, minS:maxS,:); end +mydata.vecImg = vecImg; +[Nx,Ny,Nz] = size(img); +[xi,yi,zi] = ndgrid(1:Nx,1:Ny,1:Nz); +mydata.xi = xi; mydata.yi = yi; mydata.zi = zi; +if ~isempty(mydata.mri2mni), mydata.mnipos = round(mydata.mri2mni*[mydata.pos_crop 1]'); end + whratio = 1.0187; %Width-to-height ratio w = 8.5; %Width in inches % link the callback function to new figure @@ -138,8 +137,8 @@ function myCallback(fh,~,minR, minA, minS) case mydata.h(2); mydata.pos([2 3]) = pos; case mydata.h(3); mydata.pos([1 2]) = pos; end -mydata.pos_crop = mydata.pos + [minR,minA,minS]; -mydata.mnipos = round(mydata.mri2mni*[mydata.pos_crop 1]'); +mydata.pos_crop = mydata.pos + [minR-1,minA-1,minS-1]; +if ~isempty(mydata.mri2mni), mydata.mnipos = round(mydata.mri2mni*[mydata.pos_crop 1]'); end % if new position is valid, update mydata as figure property and display % otherwise do nothing @@ -211,7 +210,7 @@ function showimages(minR,minA,minS,fh) ylim(h,[0,max(size(mydata.img))]) h(4) = nexttile(mytile,4); axis(h(4),'off'); box(h(4),'off'); clim(h(4),mydata.clim); axtoolbar(h(4),{}); -mydata.pos_crop = mydata.pos + [minR,minA,minS]; +mydata.pos_crop = mydata.pos + [minR-1,minA-1,minS-1]; % Arrange user input and coordinates in a UI Panel ip = uipanel(fh,'BackgroundColor','white','AutoResizeChildren','off'); @@ -273,16 +272,16 @@ function updateX(src,event,minR, minA, minS,fh) if ~isempty(str2double(event.Value)) switch src.Tag case 'x' - mydata.pos(1) = round(str2double(event.Value))-minR; + mydata.pos(1) = round(str2double(event.Value))-(minR-1); mydata.pos_crop(1) = str2double(event.Value); case 'y' - mydata.pos(2) = round(str2double(event.Value))-minA; + mydata.pos(2) = round(str2double(event.Value))-(minA-1); mydata.pos_crop(2) = str2double(event.Value); case 'z' - mydata.pos(3) = round(str2double(event.Value))-minS; + mydata.pos(3) = round(str2double(event.Value))-(minS-1); mydata.pos_crop(3) = str2double(event.Value); end - mydata.mnipos = round(mydata.mri2mni*[mydata.pos_crop 1]'); + if ~isempty(mydata.mri2mni), mydata.mnipos = round(mydata.mri2mni*[mydata.pos_crop 1]'); end end if ~sum( mydata.pos<1 | mydata.pos>size(mydata.img) ) set(fh,'UserData',mydata); @@ -303,9 +302,9 @@ function updateXMNI(src,event,minR, minA, minS,fh) mydata.mnipos(3) = round(str2double(event.Value)); end mydata.pos_crop = round(mydata.mri2mni\mydata.mnipos); - mydata.pos(1) = mydata.pos_crop(1)-minR; - mydata.pos(2) = mydata.pos_crop(2)-minA; - mydata.pos(3) = mydata.pos_crop(3)-minS; + mydata.pos(1) = mydata.pos_crop(1)-(minR-1); + mydata.pos(2) = mydata.pos_crop(2)-(minA-1); + mydata.pos(3) = mydata.pos_crop(3)-(minS-1); end if ~sum( mydata.pos<1 | mydata.pos>size(mydata.img) ) set(fh,'UserData',mydata); diff --git a/start_seg.m b/start_seg.m index 1857fb2..cec66ab 100644 --- a/start_seg.m +++ b/start_seg.m @@ -62,7 +62,7 @@ function start_seg(T1,T2,Template,norm) if strcmp(ext,'.hdr'), ref = [dirname filesep t1name '.img']; end t1Data = load_untouch_nii(ref); - sliceshow(t1Data.img,[],'gray',[],[],'MRI: Click anywhere to navigate.',[],[],[]); drawnow + sliceshow(t1Data.img,[],'gray',[],[],'MRI: Click anywhere to navigate.'); drawnow matlabbatch{1}.spm.spatial.preproc.channel.vols = {ref}; % image to be segmented matlabbatch{1}.spm.spatial.preproc.channel.biasreg = 0.001; % P(beta) % 0.0001; @@ -81,7 +81,7 @@ function start_seg(T1,T2,Template,norm) % fprintf('Using %s to segment T1 and T2 images: %s %s\n', Template, ref, ref2); t2Data = load_untouch_nii(ref2); - sliceshow(t2Data.img,[],'gray',[],[],'MRI: T2. Click anywhere to navigate.',[],[],[]); drawnow + sliceshow(t2Data.img,[],'gray',[],[],'MRI: T2. Click anywhere to navigate.'); drawnow matlabbatch{1}.spm.spatial.preproc.channel(2).vols = {ref2}; % the 2nd image aiding segmentation matlabbatch{1}.spm.spatial.preproc.channel(2).biasreg = 0.001; % 0.0001; diff --git a/visualizeRes.m b/visualizeRes.m index 5b4ec12..f568996 100644 --- a/visualizeRes.m +++ b/visualizeRes.m @@ -9,7 +9,7 @@ function visualizeRes(subj,subjRasRSPD,spmOut,segOut,T2,node,elem,face,inCurrent % April 2018 % August 2019 callable by roast_target() % -% (c) May 2024 Gavin Hsu and Andrew Birnbaum +% (c) May 2024, sliceshow improved by Gavin Hsu and Andrew Birnbaum if ndims(varargin{1})==3 isRoast = 1; @@ -35,11 +35,11 @@ function visualizeRes(subj,subjRasRSPD,spmOut,segOut,T2,node,elem,face,inCurrent if showAll if ~strcmp(subjName,'nyhead') disp('showing MRI and segmentations...'); - data = load_untouch_nii(subjRasRSPD); sliceshow(data.img,[],'gray',[],[],'MRI: Click anywhere to navigate.',[],mri2mni,[]); drawnow + data = load_untouch_nii(subjRasRSPD); sliceshow(data.img,[],'gray',[],[],'MRI: Click anywhere to navigate.',[],mri2mni); drawnow if ~isempty(T2) %T2 specified data = load_untouch_nii(T2); - sliceshow(data.img,[],'gray',[],[],'MRI: T2. Click anywhere to navigate.',[],mri2mni,[]); drawnow + sliceshow(data.img,[],'gray',[],[],'MRI: T2. Click anywhere to navigate.',[],mri2mni); drawnow end else disp('NEW YORK HEAD selected, there is NO MRI for it to show.') @@ -64,7 +64,7 @@ function visualizeRes(subj,subjRasRSPD,spmOut,segOut,T2,node,elem,face,inCurrent allMaskShow = masks.img; allMaskShow(gel.img>0) = numOfTissue + 1; allMaskShow(elec.img>0) = numOfTissue + 2; - sliceshow(allMaskShow,[],[],[],'Tissue index','Segmentation. Click anywhere to navigate.',[],mri2mni,[]) + sliceshow(allMaskShow,[],[],[],'Tissue index','Segmentation. Click anywhere to navigate.',[],mri2mni) drawnow end @@ -238,7 +238,7 @@ function visualizeRes(subj,subjRasRSPD,spmOut,segOut,T2,node,elem,face,inCurrent if isRoast figName = ['Voltage in Simulation: ' uniTag]; - sliceshow(vol_all.*nan_mask_brain,[],cm,[],'Voltage (mV)',[figName '. Click anywhere to navigate.'],[],mri2mni,bbox); drawnow + sliceshow(vol_all.*nan_mask_brain,round(mean(bbox)),cm,[],'Voltage (mV)',[figName '. Click anywhere to navigate.'],[],mri2mni,bbox); drawnow end for i=1:size(ef_all,4), ef_all(:,:,:,i) = ef_all(:,:,:,i).*nan_mask_brain; end @@ -246,7 +246,7 @@ function visualizeRes(subj,subjRasRSPD,spmOut,segOut,T2,node,elem,face,inCurrent dataShowVal = ef_mag(~isnan(ef_mag(:))); if isRoast figName = ['Electric field in Simulation: ' uniTag]; - sliceshow(ef_mag,[],cm,[min(dataShowVal) prctile(dataShowVal,95)],'Electric field (V/m)',[figName '. Click anywhere to navigate.'],ef_all,mri2mni,bbox); drawnow + sliceshow(ef_mag,round(mean(bbox)),cm,[min(dataShowVal) prctile(dataShowVal,95)],'Electric field (V/m)',[figName '. Click anywhere to navigate.'],ef_all,mri2mni,bbox); drawnow else for i=1:size(targetCoord,1)