FixSCMHA. #NF. Allows editing of subcortical mass (ie, wm.mgz, filled…

…, etc) by removing amyg and inf lat vent and part of hippocampus.

include/mri2.h

@@ -147,4 +147,26 @@ int QuadEulerCharChangeCheckReorder(MRI *mri, const char *testname, int decExpec
 MRI *MRIfindBrightNonWM(MRI *mri_T1, MRI *mri_wm);
 MRI *MRIzconcat(MRI *mri1, MRI *mri2, int nskip, MRI *out);
 
+/*!
+  \fn class FixSubCortMassHA
+  \brief This class "fixes" the SCM by removing stray voxels in
+  hippocampus, amyg, and inferior lateral ventricle.  The SCM can be
+  the wm.seg.mgz, wm.mgz, or filled.mgz.  The need for this stems from
+  the fact that there is a lot of WM in the hippo (eg, alveus, fimria)
+  and these often get segmented into the SCM where they create defects
+  or make the surface distorted.  All the voxels in amyg and ILV are
+  removed. Most of the voxels in hippo are removed except for those
+  within nDilate of cortex, wm, and background. These are the voxels
+  most likely to be in the wm of entorhinal cortex.
+ */
+class FixSubCortMassHA {
+public:
+  MRI *subcorticalmass=NULL; // wm.seg, wm, or filled
+  MRI *aseg=NULL;  // aseg.presurf
+  MRI *mask=NULL;  // a temporary mask
+  int nDilate = 1;
+  int FixSCM(void);
+};
+
+
 #endif

mri_edit_wm_with_aseg/mri_edit_wm_with_aseg.cpp

@@ -85,6 +85,8 @@ static int keep_edits_input = 0 ;
 static void usage_exit(int code) ;
 
 static int fcd = 0 ;
+int FixSCMHA = 0;
+int FixSCMHANdil = 0;
 
 int
 main(int argc, char *argv[])
@@ -157,6 +159,14 @@ main(int argc, char *argv[])
   remove_paths_to_cortex(mri_wm, mri_T1, mri_aseg) ;
   edit_segmentation(mri_wm, mri_T1, mri_aseg) ;
   spackle_wm_superior_to_mtl(mri_wm, mri_T1, mri_aseg) ;
+  if(FixSCMHA){
+    FixSubCortMassHA fscmha;
+    fscmha.aseg = mri_aseg;
+    fscmha.subcorticalmass = mri_wm;
+    fscmha.nDilate = FixSCMHANdil;
+    fscmha.FixSCM();
+    MRIfree(&fscmha.mask);
+  }
   if (keep_edits)
   {
     MRI *mri_old ;
@@ -254,6 +264,28 @@ get_option(int argc, char *argv[])
     keep_edits_input = 1 ;
     fprintf(stderr, "preserving editing changes in input volume...\n");
   }
+  else if (!stricmp(option, "fix-scm-ha"))
+  {
+    FixSCMHA = 1 ;
+    FixSCMHANdil = atoi(argv[2]) ; // usually set to 1
+    printf("FixSCM HA %d\n",FixSCMHANdil);
+    nargs = 1;
+  }
+  else if (!stricmp(option, "fix-scm-ha-only"))
+  {
+    // mri_edit_wm_with_aseg aseg.presurf.mgz SCM.mgz ndil out.mgz
+    // SCM = wm.seg.mgz, wm.mgz, filled.mgz, etc; ndil usually = 1
+    FixSubCortMassHA fscmha;
+    fscmha.aseg = MRIread(argv[2]);
+    if(fscmha.aseg==NULL) exit(1);
+    fscmha.subcorticalmass = MRIread(argv[3]);
+    if(fscmha.subcorticalmass==NULL) exit(1);
+    fscmha.nDilate = atoi(argv[4]) ;    
+    fscmha.FixSCM();
+    MRIfree(&fscmha.mask);
+    int err = MRIwrite(fscmha.subcorticalmass,argv[5]);
+    exit(err);
+  }
   else if (!stricmp(option, "debug_voxel"))
   {
     Gx = atoi(argv[2]) ;

mri_edit_wm_with_aseg/mri_edit_wm_with_aseg.help.xml

@@ -35,6 +35,10 @@
     </required-flagged>
     <optional-flagged>
       <argument>-fillven</argument>
+      <argument>-fix-scm-ha ndil</argument>
+      <explanation>Remove voxels in amyg, ILV, and parts of hippo</explanation>
+      <argument>-fix-scm-ha-only aseg.presurf.mgz SCM ndil out.mgz</argument>
+      <explanation>Standalone: </explanation>
       <argument>-keep</argument>
       <explanation>keep edits as found in output volume</explanation>
       <argument>-keep-in</argument>

utils/mri2.cpp

@@ -5931,3 +5931,60 @@ MRI *MRIzconcat(MRI *mri1, MRI *mri2, int nskip, MRI *out)
 
   return(out);
 }
+
+// See class definition for docs
+int FixSubCortMassHA::FixSCM(void)
+{
+  int c;
+
+  mask = MRIalloc(aseg->width,aseg->height,aseg->depth,MRI_INT);
+  MRIcopyHeader(aseg, mask);
+  MRIcopyPulseParameters(aseg, mask);
+
+  // create a mask with WM, cortex, and background. 
+  // Importantly, hippocampus is not in the mask
+#ifdef HAVE_OPENMP
+#pragma omp parallel for 
+#endif
+  for(c=0; c < aseg->width; c++){
+    int r,s;
+    for(r=0; r < aseg->height; r++){
+      for(s=0; s < aseg->depth; s++){
+	int segid = MRIgetVoxVal(aseg,c,r,s,0);
+	if(segid == 0 || segid == 2 || segid == 3 || segid == 41 || segid == 42){
+	  MRIsetVoxVal(mask,c,r,s,0,1);
+	}
+      }
+    }
+  }
+
+  // Dilate the mask. This will dilate the mask into hippocampus
+  printf("Dilating %d voxels in 3d\n",nDilate);
+  for(int n=0; n < nDilate; n++) MRIdilate(mask,mask);
+
+  // Now, remove amyg and ILV and any hippo that is not in the mask
+#ifdef HAVE_OPENMP
+#pragma omp parallel for 
+#endif
+  for(c=0; c < aseg->width; c++){
+    int r,s;
+    for(r=0; r < aseg->height; r++){
+      for(s=0; s < aseg->depth; s++){
+	int segid = MRIgetVoxVal(aseg,c,r,s,0);
+	if(IS_AMYGDALA(segid) || IS_INF_LAT_VENT(segid)){
+	  // If in amyg or ventricle, just zero the input regardless of mask
+	  MRIsetVoxVal(subcorticalmass,c,r,s,0, 0);
+	  continue;
+	}
+	if(MRIgetVoxVal(mask,c,r,s,0)>0.5) continue;
+	if(IS_HIPPO(segid) && MRIgetVoxVal(mask,c,r,s,0)<0.5){
+	  // If in hipp but not in the mask, zero
+	  MRIsetVoxVal(subcorticalmass,c,r,s,0, 0);
+	  continue;
+	}
+      }
+    }
+  }
+  return(0);
+}
+