Parahippocampus Issue

Corrected ventricle region-rowing process in cat_vol_partvol(1639) that caused defects and thickness underestimation in the parahippocampal gyri.
ChristianGaser · Jan 23, 2025 · 5778ffd · 5778ffd
1 parent 8f7f915
commit 5778ffd
Show file tree

Hide file tree

Showing 4 changed files with 103 additions and 51 deletions.
diff --git a/cat_surf_createCS2.m b/cat_surf_createCS2.m
@@ -182,17 +182,18 @@
   clear Yms;
 
   % filling
-  Ymf  = max(Ym,min(1,YMF & ~NS(Ya,opt.LAB.HC) & ~( cat_vol_morph( NS(Ya,opt.LAB.HC),'d',2) & NS(Ya,opt.LAB.TH) ))); 
-  Ymfs = cat_vol_smooth3X(Ymf,1); 
-  Ytmp = cat_vol_morph(YMF,'d',3) & Ymfs>2.3/3;
+  % RD202501: harder boundary compared to before and distance morphometry
+  Ymf  = max(Ym,min(1,YMF & ~NS(Ya,opt.LAB.HC) & ~( cat_vol_morph( NS(Ya,opt.LAB.HC),'dd',2,vx_vol) & NS(Ya,opt.LAB.TH) ))); 
+  Ymfs = cat_vol_smooth3X(Ymf,.5); % RD202501: harder boundary compared to before
+  Ytmp = cat_vol_morph(YMF,'dd',3,vx_vol) & Ymfs>2.1/3;
   Ymf(Ytmp) = max(min(Ym(Ytmp),0),Ymfs(Ytmp)); clear Ytmp Ymfs; 
   Ymf = Ymf * 3;
-  
+
   % removing fine WM structures in the hippocampus area to reduce topological and geometrical defects (added RD20190912)
   % use erode to reduce probability of cutting other gyri
-  HCmask = cat_vol_morph( NS(Ya,opt.LAB.HC) , 'de', 1.5, vx_vol); 
+  HCmask = cat_vol_morph( NS(Ya,opt.LAB.HC) , 'de', 1.5, vx_vol) & ~YMF; % RD202501: open only not filled regions
   Ymf( HCmask ) =  min(2,Ymf( HCmask )); clear HCmask; 
-
+  
   % surface output and evaluation parameter 
   Psurf = struct(); 
   res   = struct('euler_characteristic',nan,'defect_size',nan,'lh',struct(),'rh',struct()); 
@@ -338,7 +339,7 @@
       eval(sprintf('Psurf(si).%s = %s;',surffile{sfi},surffile{sfi})); 
     end
 
-    
+
     %% reduce for object area
     switch opt.surf{si}
       case {'lh'},  Ymfs = Ymf .* (Ya>0) .* ~(NS(Ya,opt.LAB.CB) | NS(Ya,opt.LAB.BV) | NS(Ya,opt.LAB.ON) | NS(Ya,opt.LAB.MB)) .* (mod(Ya,2)==1); Yside = mod(Ya,2)==1; 
@@ -355,7 +356,7 @@
     if ~iscerebellum 
       % RD202107:  Use atlas and Shooting template information to close the
       %            hippocampal gyrus but open the hippocampal region.
-      mask_parahipp = NS(Ya,opt.LAB.PH) | NS(Ya,opt.LAB.HC) | NS(Ya,opt.LAB.VT); 
+      mask_parahipp = opt.close_parahipp & (NS(Ya,opt.LAB.PH) | NS(Ya,opt.LAB.HC) | NS(Ya,opt.LAB.VT));
       if isempty(Ytemplate)
         VT = NS(Ya,opt.LAB.PH) | NS(Ya,opt.LAB.HC); % this does not work but also should not create problems
       else
@@ -1979,7 +1980,7 @@ function saveSurf(CS,P)
 
 
 
-    %% RD202107: closing jof parahippocampal gyri
+    %% RD202107: closing parahippocampal gyri
     %  --------------------------------------------------------------------
     %  We are mostly interested in closing of holes in deep parahippocampal 
     %  regions. Hence, we will limit our operations by an enlargements of

diff --git a/cat_surf_createCS3.m b/cat_surf_createCS3.m
@@ -155,15 +155,16 @@
   clear Yms;
 
   % filling
-  Ymf  = max(Ym,min(1,YMF & ~NS(Ya,opt.LAB.HC) & ~( cat_vol_morph( NS(Ya,opt.LAB.HC),'d',2) & NS(Ya,opt.LAB.TH) ))); 
-  Ymfs = cat_vol_smooth3X(Ymf,1); 
-  Ytmp = cat_vol_morph(YMF,'d',3) & Ymfs>2.3/3;
+  % RD202501: harder boundary compared to before and distance morphometry
+  Ymf  = max(Ym,min(1,YMF & ~NS(Ya,opt.LAB.HC) & ~( cat_vol_morph( NS(Ya,opt.LAB.HC),'dd',2,vx_vol) & NS(Ya,opt.LAB.TH) ))); 
+  Ymfs = cat_vol_smooth3X(Ymf,.5); % RD202501: harder boundary compared to before
+  Ytmp = cat_vol_morph(YMF,'dd',3,vx_vol) & Ymfs>2.1/3;
   Ymf(Ytmp) = max(min(Ym(Ytmp),0),Ymfs(Ytmp)); clear Ytmp Ymfs; 
   Ymf = Ymf * 3;
 
   % removing fine WM structures in the hippocampus area to reduce topological and geometrical defects (added RD20190912)
   % use erode to reduce probability of cutting other gyri
-  HCmask = cat_vol_morph( NS(Ya,opt.LAB.HC) , 'de', 1.5, vx_vol); 
+  HCmask = cat_vol_morph( NS(Ya,opt.LAB.HC) , 'de', 1.5, vx_vol) & ~YMF; % RD202501: open only not filled regions
   Ymf( HCmask ) =  min(2,Ymf( HCmask )); clear HCmask; 
 
   % surface output and evaluation parameter 
@@ -291,7 +292,7 @@
     if ~iscerebellum 
       % RD202107:  Use atlas and Shooting template information to close the
       %            hippocampal gyrus but open the hippocampal region.
-      mask_parahipp = NS(Ya,opt.LAB.PH) | NS(Ya,opt.LAB.HC) | NS(Ya,opt.LAB.VT); 
+      mask_parahipp = opt.close_parahipp & NS(Ya,opt.LAB.PH) | NS(Ya,opt.LAB.HC) | NS(Ya,opt.LAB.VT); 
       if isempty(Ytemplate)
         VT = NS(Ya,opt.LAB.PH) | NS(Ya,opt.LAB.HC); % this does not work but also should not create problems
       else

diff --git a/cat_vol_partvol.m b/cat_vol_partvol.m
@@ -328,7 +328,9 @@
     Ybg = cat_vol_laplace3R(Ybg,Ybg==1.5,0.005)<1.5 & Ym<2.9 & Ym>1.8 & Ydiv>-0.02;
     Ya1(Ybg)=LAB.BG;                                                         % basal ganglia
     Ya1(YA==LAB.TH & Ym>1.9 & Ym<2.85 & Ydiv>-0.1)=LAB.TH;                   % thalamus
-    Ya1(YA==LAB.HC & Ym>1.9 & Ym<2.85 & Ydiv>-0.1)=LAB.HC;                   % hippocampus
+    % RD202501: hippocampus definition was not optimal here as we need a clear parahippocampus for surface reconstruction 
+    Ya1( cat_vol_morph(YA==LAB.HC,'dd',3) & Ym>1.5 & Ym<2.5 & ~(cat_vol_morph(YA==LAB.PH,'dd',2) | ...
+      YA==LAB.TH | cat_vol_morph(YA==LAB.NV,'dd',5,vx_vol))) = LAB.HC;  % hippocampus & ~cat_vol_morph(YA==LAB.PH,'d')
     NBVC = BVCstr > 0 && BVCstr <= 1; % RD202306: new BV correction by default
     if NBVC
       % define some high intensity BVs and the neocortex
@@ -345,12 +347,20 @@
          (Yp0>1.5 & Ym<3.5 & Ybgd>1 & Ybgd<8 & (Ybgd>4 | ...
          Ydiv<-0.02+Ybgd/200))) & Ya1==0)=LAB.CT;                            % cerebrum
     end
-    Ya1((Yp0>2.0 & Ym>2.0) & YA==LAB.CB)=LAB.CB;                             % cerebellum
+    % use region-growing in case of the cerebellum to compensate for fine structure
+    Ycb  = single(YA==LAB.CB) + 2*single(cat_vol_morph( YA==LAB.CT,'de',6)); 
+    Ycb(Yp0<1 & Ym<2) = nan; Ycb(cat_vol_morph( YA~=LAB.CB,'de',12)) = nan; 
+    [Ycb,Yd] = cat_vol_downcut(Ycb,Ym,noise); 
+    Ycb = single(YA==LAB.CB) | cat_vol_smooth3X(Ycb==1 & Yd<100,2)>.7; 
+    % set other regions
+    Ya1((Yp0>2.0 & Ym>2.0  & Ym<3.1) & Ycb)=LAB.CB;                          % cerebellum
     Ya1((Yp0>2.0 & Ym>2.0) & YA==LAB.BS)=LAB.BS;                             % brainstem
     Ya1((Yp0>2.0 & Ym>2.0) & YA==LAB.ON)=LAB.ON;                             % optical nerv
     Ya1((Ya1==0 & Yp0<1.5 & Ym<1.5 & Yp0>1.3 & Ym>1.3) & YA==LAB.BV)=LAB.BV; % low-int VB  (updated due to cerebellar erros RD20190929)
     Ya1((Ya1==0 & Yp0>3.0 & Ym>3.2) & YA==LAB.BV)=LAB.BV;                    % high-int VB (updated due to cerebellar erros RD20190929)
     Ya1((Yp0>2.0 & Ym>2.0) & YA==LAB.MB)=LAB.MB;                             % midbrain
+    Ya1((Yp0>2.0 & Ym>2.0) & cat_vol_smooth3X(YA==LAB.MB,2)>.1)=LAB.MB;      % midbrain
+    Ya1((Yp0>2.0 & Ym>2.0) & YA==LAB.PH & Ya1==0)=LAB.CT;                    % added PH as cortex
     %%
     if NBVC
       % create a cerebellar mask to avoid corrections there
@@ -371,9 +381,10 @@
       [~,Yd] = cat_vol_downcut(Ya0,Ym,noise  ); clear Ya0 
       Ya0 = single(Ya1>0 & Ya1~=LAB.BV); Ya0(Ym<1.7) = nan; 
       [~,Yd2] = cat_vol_downcut(Ya0,Ym,noise * 16); clear Ya0
-      Ya1( ((Yd - Yd2)>10) & Ym>2.4 & YA==LAB.CT & ~Ycb3 & ~Yhc3) = LAB.BV; % add highly distant voxels
+      Ya1( ((Yd - Yd2)>10./YbvA) & Ym>2.4 & YA==LAB.CT & ~Ycb3 & ~Yhc3 & YbvA>.7) = LAB.BV; % add highly distant voxels
       Ya1(Ya1==0 & Ym>2.4 & Yd>70 & Yd2>50 & ~cat_vol_morph(YA>1 & YA~=LAB.HD,'d') & ~Ycb3 & ~Yhc3) = LAB.BV; % mask further BV
       Ya1(Ya1==0 & Ym>2.2 & Yd<20 & ~cat_vol_morph(YA>1 & YA~=LAB.HD,'d') & ~Ycb3) = LAB.CT; % add also save brain voxels
+      Ya1(Ya1==LAB.BV & cat_vol_morph(Ya1==LAB.BV,'l',[inf 16])==0) = 0;
       clear Ycb3 Yhc3
 
       %% modified old lines
@@ -407,6 +418,10 @@
     Ya1((Yp0>1.75 & Ym>1.75 & Yp0<2.5 & Ym<2.5) & Ya1==LAB.MB)=0;            % midbrain correction
     Ya1(Ya1==LAB.CT & ~cat_vol_morph(Ya1==LAB.CT,'do',1.4)) = 0; 
 
+    % correction of structures that should be compact
+    Ya1(Ya1==LAB.BS & ~cat_vol_morph(cat_vol_morph(Ya1==LAB.BS,'o'),'c')) = 0; 
+    Ya1(Ya1==LAB.MB & ~cat_vol_morph(cat_vol_morph(Ya1==LAB.MB,'o'),'c')) = 0; 
+    Ya1(Ya1==LAB.PH) = 0; 
 
     %% Mapping of ventricles:
     %  Ventricle estimation with a previous definition of non ventricle CSF
@@ -415,16 +430,19 @@
     %  because dartle failed for large ventricle. 
     %  It is important to use labopen for each side!
     stime = cat_io_cmd('  Ventricle detection','g5','',verb,stime); dispc=dispc+1;
-    Ynv = cat_vol_morph(cat_vol_morph(~Yb,'d',4,vx_vol) | ...
-          (YA==LAB.CB | YA==LAB.BS),'d',2,vx_vol) | cat_vol_morph(YA==LAB.NV,'e',1,vx_vol);
-    Ynv = Ynv | (~cat_vol_morph((Yp0>0 & Yp0<1.5 & (YA==LAB.VT)),'dd',20,vx_vol) & Yp0<1.5);
-    Ynv = single(Ynv & Ym<2 & ~cat_vol_morph(Yp0<2 & (YA==LAB.VT) & Yg<0.2,'d',4,vx_vol));
+    %% RD202501: added parahypocampus here to avoid overgrowing ventricles and filling issues and defects 
+    Yph = Ym>2 & cat_vol_smooth3X(YA==LAB.PH | YA==LAB.HC,2)>.1; 
+    Ynv = cat_vol_morph(~Yb,'d',4,vx_vol) | cat_vol_morph(YA==LAB.CB,'dd',10,vx_vol) | ...% RD202501: extend CB 
+          cat_vol_morph(YA==LAB.BS,'dd',2,vx_vol) | cat_vol_morph(YA==LAB.NV,'e',1,vx_vol) | cat_vol_morph(YA==LAB.TH & Ym<2.5,'e',1,vx_vol); 
+    Ynv = Ynv | Yph | (~cat_vol_morph((Yp0>0 & Yp0<1.5 & (YA==LAB.VT)),'dd',20,vx_vol) & Yp0<1.5);
+    Ynv = single(Ynv & Ym<2 & ~cat_vol_morph(Yp0<2 & (YA==LAB.VT) & Yg<0.2,'d',4,vx_vol) & ~Yph);
     Ynv = smooth3(round(Ynv))>0.5; 
     % between thamlamus
     Ynv = Ynv | (cat_vol_morph(Ya1==LAB.TH,'c',10) & Yp0<2) | YA==LAB.CB | YA==LAB.BS;
     Ynv = smooth3(Ynv)>0.8;
     Yvt = single(smooth3(Yp0<1.5 & (YA==LAB.VT) & Yg<0.25 & ~Ynv)>0.7); 
     Yvt(Yvt==0 & Ynv)=2; Yvt(Yvt==0 & Ym>1.8)=nan; Yvt(Yvt==0)=1.5;
+    Yvt(Yvt==0 & Ydiv<-0.1) = nan; 
 
     %% subcortical stroke lesions
     if exist('Ylesionmsk','var'), Yvt(Ylesionmsk>0.5) = nan; end
@@ -450,8 +468,10 @@
     %% bottleneck
     Yvt2 = cat_vol_laplace3R(Yvt,Yvt==1.5,0.01); % first growing for large regions
     Yvt(cat_vol_morph(Yvt2<1.4,'o',2) & ~isnan(Yvt) & Yp0<1.5) = 1; 
-    Yvt(cat_vol_morph(Yvt2>1.8,'o',2) & ~isnan(Yvt) & Yp0<1.5) = 2; 
-    Yvt2 = cat_vol_laplace3R(Yvt,Yvt==1.5,0.002);
+    Yvt(cat_vol_morph(Yvt2>1.6,'o',2) & ~isnan(Yvt) & Yp0<1.5) = 2; 
+    Ygx = Ym/3 ./ cat_vol_localstat(Ym/3,Ym>1.125,1,3,1); Yvt(Yvt==1.5 & Ygx<.5)=nan; 
+    Ygx = Ym/3 ./ cat_vol_localstat(Ym/3,Ym>1.25,1,3,1); Yvt(Yvt==1.5 & Ygx<.8)=nan; 
+    Yvt2 = cat_vol_laplace3R(Yvt,Yvt==1.5,0.001);
     % remove small objects
     warning('off','MATLAB:cat_vol_morph:NoObject');
     Yvt = cat_vol_morph(Yvt2<1.5, 'l', [10 0.1]);
@@ -851,13 +871,18 @@
     [tmp0,tmp1,Ya1] = cat_vbdist(Ya1,Yb); clear tmp0 tmp1;
     Ya1(Ybv) = LAB.BV; 
 
-    % consider gyrus parahippocampalis
-    Ya1(YA==LAB.PH) = LAB.PH;
+    % consider gyrus parahippocampalis | cat_vol_morph(Ya1==LAB.HC,'e')
+    % RD202501: add defintion of parahippocampal gyrus (with some extensive 
+    %           processing to really have a whole free thing but we will try to first keep it simple)
+    Yph = cat_vol_morph((YA==LAB.PH ),'dd',1.9) & cat_vol_morph(Ya1==LAB.VT | Ya1==LAB.HC,'dd',4,vx_vol) & Ym>2.125 & Ydiv<0.05; 
+    Ya1(Yph>0) = LAB.PH; clear Yph; % parahippocampus
 
     %% side aligment using laplace to correct for missalignments due to the normalization
     stime = cat_io_cmd('  Side alignment','g5','',verb,stime); dispc=dispc+1;
     YBG  = Ya1==LAB.BG | Ya1==LAB.TH;
-    YMF  = Ya1==LAB.VT | Ya1==LAB.BG | Ya1==LAB.TH | Ya1==LAB.HI | Ya1==LAB.TH; % add the thalamus (RD20190913)
+    YMF  = Ya1==LAB.VT | Ya1==LAB.BG | Ya1==LAB.TH | Ya1==LAB.HI | Ya1==LAB.PH; % add the thalamus (RD20190913)
+    YMF(smooth3(YMF)<.5) = 0; 
+    YMF  = cat_vol_morph(cat_vol_morph(YMF,'c',2),'o'); 
     YMF2 = cat_vol_morph(YMF,'d',2,vx_vol) | Ya1==LAB.CB | Ya1==LAB.BS | Ya1==LAB.MB;
     Ymf  = max(Ym,smooth3(single(YMF2*3))); 
     Ycenter = cat_vol_smooth3X(YS==0,6)<0.9 & cat_vol_smooth3X(YS==1,6)<0.9 & ~YMF2 & Yp0>0 & Ym<3.1 & (Yp0<2.5 | Ya1==LAB.BV);
@@ -868,9 +893,10 @@
 
     %% YMF for FreeSurfer fsaverage
     Ysm  = cat_vol_morph(Ys==2,'d',1.75,vx_vol) & cat_vol_morph(Ys==1,'d',1.75,vx_vol);
-    YMF  = cat_vol_morph(Ya1==LAB.VT | Ya1==LAB.BG | Ya1==LAB.HI | (Ya1==LAB.TH & smooth3(Yp0)>1.25),'c',3,vx_vol) & ~Ysm; % changed thalamus (RD20190913)
-    %YMF  = YMF | (cat_vol_morph(YA==LAB.CT & YBG,'c',6) & ~Ysm); 
-    YMF  = Ym<=2.5  & cat_vol_morph(YMF | Ym>2.3,'c',1) & cat_vol_morph(YMF,'d',2,vx_vol);
+    YMF  = cat_vol_morph(Ya1==LAB.VT | Ya1==LAB.PH | Ya1==LAB.BG | Ya1==LAB.HI | ...
+      (Ya1==LAB.TH & smooth3(Yp0)>1.25),'dc',2,vx_vol) & ~Ysm; % changed thalamus (RD20190913)
+    YMF  = Ym<=2.75  & cat_vol_morph(YMF | Ym>2.3,'c',1) & cat_vol_morph(YMF,'dd',2,vx_vol);
+    %YMF  = cat_vol_morph(YMF,'do',1.5);
     YMF  = smooth3(YMF)>0.5;
     clear Ysm;