diff --git a/include/mri2.h b/include/mri2.h index a918a24ae20..17c50895ca1 100644 --- a/include/mri2.h +++ b/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 diff --git a/mri_edit_wm_with_aseg/mri_edit_wm_with_aseg.cpp b/mri_edit_wm_with_aseg/mri_edit_wm_with_aseg.cpp index c1f0bbce4c6..f6ebb4d005e 100644 --- a/mri_edit_wm_with_aseg/mri_edit_wm_with_aseg.cpp +++ b/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]) ; diff --git a/mri_edit_wm_with_aseg/mri_edit_wm_with_aseg.help.xml b/mri_edit_wm_with_aseg/mri_edit_wm_with_aseg.help.xml index b212e72f183..308bbf36c9c 100644 --- a/mri_edit_wm_with_aseg/mri_edit_wm_with_aseg.help.xml +++ b/mri_edit_wm_with_aseg/mri_edit_wm_with_aseg.help.xml @@ -35,6 +35,10 @@ -fillven + -fix-scm-ha ndil + Remove voxels in amyg, ILV, and parts of hippo + -fix-scm-ha-only aseg.presurf.mgz SCM ndil out.mgz + Standalone: -keep keep edits as found in output volume -keep-in diff --git a/utils/mri2.cpp b/utils/mri2.cpp index 26cd075ba3f..382b288a051 100644 --- a/utils/mri2.cpp +++ b/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); +} +