fix some issues of the new sliceshow

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];

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);

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:
 

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,15 +245,15 @@ 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
             if ~exist(optRoast.T2,'file')
                 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
 

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);

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;

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,15 +238,15 @@ 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
 ef_mag = ef_mag.*nan_mask_brain;
 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)