diff --git a/.github/workflows/build.yml b/.github/workflows/build.yml index 86e2a6f..e161a32 100644 --- a/.github/workflows/build.yml +++ b/.github/workflows/build.yml @@ -15,7 +15,7 @@ jobs: fail-fast: false matrix: os: [ubuntu, macos, windows] - python-version: ["3.9", "3.10", "3.11", "3.12"] + python-version: ["3.10", "3.11", "3.12", "3.13", "3.14"] name: ${{ matrix.os }} - py${{ matrix.python-version }} runs-on: ${{ matrix.os }}-latest defaults: diff --git a/README.md b/README.md index 4a44d7b..91c6e2f 100644 --- a/README.md +++ b/README.md @@ -11,7 +11,7 @@ hippocampal shape and thickness analysis as described in our recent [publication Please see the [documentation pages](https://deep-mi.org/hipsta/dev/index.html) for a general overview, usage information and examples, and output description. Brief usage information is also available [here](hipsta/doc/DOCUMENTATION.md). -Some suggestions for running the script can be found in the [tutorial](TUTORIAL.md). +Some suggestions for running the script can be found in the [tutorial](https://deep-mi.org/hipsta/dev/tutorials/index.html). ## Current status: @@ -75,7 +75,7 @@ exist as an environment variable and point to a valid FreeSurfer installation. or the ASHS software. A custom segmentation is also permissible (some restrictions and settings apply; see [Supported Segmentations](https://github.com/Deep-MI/Hipsta#supported-segmentations)). -3. Python 3.8 or higher including the lapy, numpy, scipy, nibabel, pyvista, and +3. Python 3.10 or higher including the lapy, numpy, scipy, nibabel, pyvista, and pyacvd libraries, among others. See `requirements.txt` for a full list, and use `pip install -r requirements.txt` to install. diff --git a/TUTORIAL.md b/TUTORIAL.md index 5c06139..18d135e 100644 --- a/TUTORIAL.md +++ b/TUTORIAL.md @@ -239,9 +239,9 @@ argument. Since no explicit labels for the hippocampal head are present in ASHS segmentations, it will be necessary to use the `--automask-head` flag and -optionally also the `--automask-tail-margin` argument. If the identification of -the boundary between hippocampal body tail does not work well given the existing -labels, using the `--automask-tail` flag and optionally the `--automask-head-margin` +optionally also the `--automask-head-margin` argument. If the identification of +the boundary between hippocampal body and tail does not work well given the existing +labels, using the `--automask-tail` flag and optionally the `--automask-tail-margin` argument may be advised as well. ### Additional arguments for ASHS and FreeSurfer segmentations @@ -264,10 +264,10 @@ The primary output of the hippocampal shape and thickness analysis are surface files and associated thickness values in the `thickness` folder. The estimated hippocampal thickness values will be stored in csv tables: -| Filenane | Contents | +| Filename | Contents | |----------------------------------------------------|-----------------------------------------------------------| | `lh.grid-segments-x.csv`, `rh.grid-segments-x.csv` | Length estimates in the medial-to-lateral dimension | -| `lh.grid-segments-y.csv`, `rh.grid-segments-y.csv` | Length estimates in the medial-to-lateral dimension | +| `lh.grid-segments-y.csv`, `rh.grid-segments-y.csv` | Length estimates in the posterior-to-anterior dimension | | `lh.grid-segments-z.csv`, `rh.grid-segments-z.csv` | Thickness estimates in the exterior-to-interior dimension | Here, x corresponds to the medial-->lateral dimension, y to the @@ -303,6 +303,57 @@ provided for interior, mid, and exterior surfaces: | `.-surface.-curv.mgh` | Overlays of mean and gaussian curvature estimtes for corresponding surfaces | +## Extensions to parahippocampal and entorhinal cortex + +The Penn ABC-3T ASHS atlas includes several labels beyond the hippocampus, +including parahippocampal, entorhinal, BA35, and BA36 regions. The analysis can +extended to parts of these regions as well. The corresponding processing +variants can be called with the `--lut ashs-penn_abc_3t_t2_ext` and `--lut ashs-penn_abc_3t_t2_ent` +arguments. The former will add (parts of) parahippocampal and entorhinal cortex +as well as BA35 to the analysis, but still restrict the analysis to the +hippocampal body region along the longitudinal axis. +The latter will perform a separate analysis of the entorhinal cortex. Due to +more complex or variable anatomy, these approaches should to some extent be +considered experimental and may require additional QC or manual edits. + +The expected output when using `--lut ashs-penn_abc_3t_t2_ext` is shown in the +following figure; it is simply an extension of the analysis beyond the +subiculum and presubiculum into the cortex. + +![](images/scaled/ctx-01.png) + +A potential issue that might arise in this analysis could be due to large shifts +between segmentations in successive slices, especially for anisotropic voxel +sizes, leading to holes in the surface mesh. Another potential issue is too close +proximity of the subiculum / CA1 and cortical regions so that the surface mesh +will incorrectly create a connection. Manual edits or changing the +spatial filter settings may help resolving these issues (see the following +section on Troubleshooting for details). + +![](images/scaled/ctx-03a.png) +![](images/scaled/ctx-03b.png) +![](images/scaled/ctx-03c.png) +![](images/scaled/ctx-02.png) + +The expected output when using `--lut ashs-penn_abc_3t_t2_ent` is shown below. +Here, we leave out the hippocampal labels and restrict the analysis the +complete entorhinal cortex segmentation. + +![](images/scaled/erc-01a.png) +![](images/scaled/erc-01b.png) +![](images/scaled/erc-01c.png) + +Care should be taken that the entorhinal cortex is cut across its full extent +at its anterior and posterior ends. The left figure shows to conservative +cutting, the middle and right figure more appropriate cuts. Cutting can be +adjusted using the `--cut-range` and `--automask-head-margin` / `--automask-tail-margin` +arguments (see the following section on Troubleshooting for details). + +![](images/scaled/erc-02a.png) +![](images/scaled/erc-02b.png) +![](images/scaled/erc-02c.png) + + ## Troubleshooting ### Logfiles diff --git a/VERSION b/VERSION index bf42b94..2774f85 100644 --- a/VERSION +++ b/VERSION @@ -1 +1 @@ -0.9.1.3 +0.10.0 \ No newline at end of file diff --git a/doc/tutorials/index.md b/doc/tutorials/index.md index ad0b2f1..f663550 100644 --- a/doc/tutorials/index.md +++ b/doc/tutorials/index.md @@ -1,5 +1,3 @@ -======= - ```{include} ../../TUTORIAL.md :relative-images: ``` \ No newline at end of file diff --git a/hipsta/cfg/atlases.py b/hipsta/cfg/atlases.py index 481fb07..15c0f59 100644 --- a/hipsta/cfg/atlases.py +++ b/hipsta/cfg/atlases.py @@ -26,7 +26,7 @@ def get_atlases(lut): LUTINDEX = [226, [233, 235, 237, 239, 245], 234, 236, 238, 240, 240, 242, 244, 246] - LUTDICT = dict(zip(LUTLABEL, LUTINDEX)) + LUTDICT = dict(zip(LUTLABEL, LUTINDEX, strict=True)) HSFLIST = [234, 236, 238, 240, 246] @@ -37,7 +37,7 @@ def get_atlases(lut): LUTINDEX = [226, [233, 235, 237, 239, 245], 234, 236, 238, 240, 240, 242, 244] - LUTDICT = dict(zip(LUTLABEL, LUTINDEX)) + LUTDICT = dict(zip(LUTLABEL, LUTINDEX, strict=True)) HSFLIST = [234, 236, 238, 240] @@ -63,10 +63,62 @@ def get_atlases(lut): LUTINDEX = [1, 2, 4, 3, 3, 5, 8, 8, 9, 10, 11, 12, 20, 5] - LUTDICT = dict(zip(LUTLABEL, LUTINDEX)) + LUTDICT = dict(zip(LUTLABEL, LUTINDEX, strict=True)) HSFLIST = [8, 1, 2, 4] + elif lut == "ashs-penn_abc_3t_t2_ent": + LOGGER.info("Found internal, modified look-up table for the Penn ABC-3T ASHS Atlas for T2-weighted MRI.") + + LUTLABEL = [ + "ca1", + "ca2", + "ca3", + "ca4", + "dg", + "tail_orig", + "subiculum", + "presubiculum", + "entorhinal", + "ba35", + "ba36", + "parahippocampal", + "head", + "tail", + ] + + LUTINDEX = [1, 2, 4, 3, 3, 5, 8, 8, 9, 10, 11, 12, 20, 5] + + LUTDICT = dict(zip(LUTLABEL, LUTINDEX, strict=True)) + + HSFLIST = [9] + + elif lut == "ashs-penn_abc_3t_t2_ext": + LOGGER.info("Found internal, modified look-up table for the Penn ABC-3T ASHS Atlas for T2-weighted MRI.") + + LUTLABEL = [ + "ca1", + "ca2", + "ca3", + "ca4", + "dg", + "tail_orig", + "subiculum", + "presubiculum", + "entorhinal", + "ba35", + "ba36", + "parahippocampal", + "head", + "tail", + ] + + LUTINDEX = [1, 2, 4, 3, 3, 5, 8, 8, 9, 10, 11, 12, 20, 5] + + LUTDICT = dict(zip(LUTLABEL, LUTINDEX, strict=True)) + + HSFLIST = [12, 10, 9, 8, 1, 2, 4] + elif lut == "ashs-umcutrecht_7t": LOGGER.info("Found internal, modified look-up table for ASHS UMC Utrecht 7T atlas.") @@ -86,7 +138,7 @@ def get_atlases(lut): LUTINDEX = [1, 2, 2, 3, 4, 5, 6, 5, 7, 8, 20] - LUTDICT = dict(zip(LUTLABEL, LUTINDEX)) + LUTDICT = dict(zip(LUTLABEL, LUTINDEX, strict=True)) HSFLIST = [2, 3, 4, 6] @@ -105,7 +157,7 @@ def get_atlases(lut): lut = lut[lut[1].str.lower().isin(["presubiculum", "subiculum", "head", "tail", "ca1", "ca2", "ca3", "ca4"])] - LUTDICT = dict(zip(lut[1], lut[0])) + LUTDICT = dict(zip(lut[1], lut[0], strict=True)) HSFLIST = list(lut[1]) @@ -121,6 +173,7 @@ def get_atlases(lut): LUTDICT["bndtail"] = 2260 LUTDICT["bndhead"] = 2320 LUTDICT["bndca4"] = 2420 + LUTDICT["bndba35"] = 100 # return diff --git a/hipsta/cfg/config.py b/hipsta/cfg/config.py index 3f5e8e3..6c30ce4 100644 --- a/hipsta/cfg/config.py +++ b/hipsta/cfg/config.py @@ -1,5 +1,5 @@ """ -This module provides configuration defaults for the hippocampal shape and +This module provides configuration defaults for the hippocampal shape and thickness analysis package. """ @@ -10,6 +10,8 @@ def get_defaults(x): # note: if changing these values, remember to adjust the argparse help # messages as well as the run_hipsta docstring also # options + start_with_edited_labels=False, + start_with_edited_masks=False, no_cleanup=False, no_crop=False, upsample=False, diff --git a/hipsta/computeCubeParam.py b/hipsta/computeCubeParam.py index 2ba962d..7b02992 100755 --- a/hipsta/computeCubeParam.py +++ b/hipsta/computeCubeParam.py @@ -12,8 +12,6 @@ from scipy import sparse as sp from scipy import stats as st -# from sklearn.decomposition import PCA - # ============================================================================== # LOGGING @@ -39,7 +37,7 @@ def getSeam(v4c, t4c, i4c, k4c, v4cBndOpenKeep, t4cBndOpen, anisoLaplEvec): # getSeam subfunction #1 # case 1: one nan in the tetra - def getSeamCase1(v4c, t4c, i4c, k4c, vfuncXEv1, e4cBndOpen, newTetra, newVtcs, newVtcsAdj, newVtcsSgn, i): + def getSeamCase1(v4c, t4c, i4c, k4c, vfuncXEv1, e4cBndOpen, newTetra, newVtcs, newVtcsAdj, i): # do we have one or two neg values? # case 1a: one neg value if np.sum(vfuncXEv1[t4c[i, :]] < 0) == 1: @@ -88,16 +86,7 @@ def getSeamCase1(v4c, t4c, i4c, k4c, vfuncXEv1, e4cBndOpen, newTetra, newVtcs, n # also append i4c and k4c i4c = np.concatenate((i4c, st.mode(i4c[t4c[i, :]])[0])) k4c = np.concatenate((k4c, st.mode(k4c[t4c[i, :]])[0])) - # let's remember the new indices, and whether or not it's a - # 'positive' or 'negative' tetra (in terms of the 2nd - # eigenfunction) newVtcs = np.concatenate((newVtcs, [idxPt1])) - tmpSgn = np.unique(vfuncXEv2[t4c[i, :]][np.logical_not(np.isnan(vfuncXEv2[t4c[i, :]]))]) - if len(tmpSgn) == 1: - newVtcsSgn = np.concatenate((newVtcsSgn, tmpSgn)) - else: - newVtcsSgn = np.concatenate((newVtcsSgn, np.array([0]))) - raise RuntimeError("Inconsistency while creating seam 1, exiting.") else: idxPt1 = newVtcsAdj[t4c[i, idxPos1], t4c[i, idxNeg1]] @@ -113,16 +102,7 @@ def getSeamCase1(v4c, t4c, i4c, k4c, vfuncXEv1, e4cBndOpen, newTetra, newVtcs, n # also append i4c and k4c i4c = np.concatenate((i4c, st.mode(i4c[t4c[i, :]])[0])) k4c = np.concatenate((k4c, st.mode(k4c[t4c[i, :]])[0])) - # let's remember the new indices, and whether or not it's a - # 'positive' or 'negative' tetra (in terms of the 2nd - # eigenfunction) newVtcs = np.concatenate((newVtcs, [idxPt2])) - tmpSgn = np.unique(vfuncXEv2[t4c[i, :]][np.logical_not(np.isnan(vfuncXEv2[t4c[i, :]]))]) - if len(tmpSgn) == 1: - newVtcsSgn = np.concatenate((newVtcsSgn, tmpSgn)) - else: - newVtcsSgn = np.concatenate((newVtcsSgn, np.array([0]))) - raise RuntimeError("Inconsistency while creating seam 2, exiting.") else: idxPt2 = newVtcsAdj[t4c[i, idxPos2], t4c[i, idxNeg1]] @@ -184,16 +164,7 @@ def getSeamCase1(v4c, t4c, i4c, k4c, vfuncXEv1, e4cBndOpen, newTetra, newVtcs, n # also append i4c and k4c i4c = np.concatenate((i4c, st.mode(i4c[t4c[i, :]])[0])) k4c = np.concatenate((k4c, st.mode(k4c[t4c[i, :]])[0])) - # let's remember the new indices, and whether or not it's a - # 'positive' or 'negative' tetra (in terms of the 2nd - # eigenfunction) newVtcs = np.concatenate((newVtcs, [idxPt1])) - tmpSgn = np.unique(vfuncXEv2[t4c[i, :]][np.logical_not(np.isnan(vfuncXEv2[t4c[i, :]]))]) - if len(tmpSgn) == 1: - newVtcsSgn = np.concatenate((newVtcsSgn, tmpSgn)) - else: - newVtcsSgn = np.concatenate((newVtcsSgn, np.array([0]))) - raise RuntimeError("Inconsistency while creating seam 3, exiting.") else: idxPt1 = newVtcsAdj[t4c[i, idxPos1], t4c[i, idxNeg1]] @@ -209,16 +180,7 @@ def getSeamCase1(v4c, t4c, i4c, k4c, vfuncXEv1, e4cBndOpen, newTetra, newVtcs, n # also append i4c and k4c i4c = np.concatenate((i4c, st.mode(i4c[t4c[i, :]])[0])) k4c = np.concatenate((k4c, st.mode(k4c[t4c[i, :]])[0])) - # let's remember the new indices, and whether or not it's a - # 'positive' or 'negative' tetra (in terms of the 2nd - # eigenfunction) newVtcs = np.concatenate((newVtcs, [idxPt2])) - tmpSgn = np.unique(vfuncXEv2[t4c[i, :]][np.logical_not(np.isnan(vfuncXEv2[t4c[i, :]]))]) - if len(tmpSgn) == 1: - newVtcsSgn = np.concatenate((newVtcsSgn, tmpSgn)) - else: - newVtcsSgn = np.concatenate((newVtcsSgn, np.array([0]))) - raise RuntimeError("Inconsistency while creating seam 4, exiting.") else: idxPt2 = newVtcsAdj[t4c[i, idxPos1], t4c[i, idxNeg2]] @@ -237,13 +199,13 @@ def getSeamCase1(v4c, t4c, i4c, k4c, vfuncXEv1, e4cBndOpen, newTetra, newVtcs, n raise RuntimeError("Inconsistency while creating seam (incorrect number of negative values), exiting.") # return - return v4c, t4c, i4c, k4c, e4cBndOpen, newTetra, newVtcs, newVtcsAdj, newVtcsSgn + return v4c, t4c, i4c, k4c, e4cBndOpen, newTetra, newVtcs, newVtcsAdj # -------------------------------------------------------------------------- # getSeam subfunction #2 # case 2: two nans in tetra - def getSeamCase2(v4c, t4c, i4c, k4c, vfuncXEv1, e4cBndOpen, newTetra, newVtcs, newVtcsAdj, newVtcsSgn, i): + def getSeamCase2(v4c, t4c, i4c, k4c, vfuncXEv1, e4cBndOpen, newTetra, newVtcs, newVtcsAdj, i): # get indices idxPos1 = np.where(vfuncXEv1[t4c[i, :]] > 0)[0][0] idxNeg1 = np.where(vfuncXEv1[t4c[i, :]] < 0)[0][0] @@ -279,16 +241,7 @@ def getSeamCase2(v4c, t4c, i4c, k4c, vfuncXEv1, e4cBndOpen, newTetra, newVtcs, n # also append i4c and k4c i4c = np.concatenate((i4c, st.mode(i4c[t4c[i, :]])[0])) k4c = np.concatenate((k4c, st.mode(k4c[t4c[i, :]])[0])) - # let's remember the new indices, and whether or not it's a - # 'positive' or 'negative' tetra (in terms of the 2nd - # eigenfunction) newVtcs = np.concatenate((newVtcs, [idxPt1])) - tmpSgn = np.unique(vfuncXEv2[t4c[i, :]][np.logical_not(np.isnan(vfuncXEv2[t4c[i, :]]))]) - if len(tmpSgn) == 1: - newVtcsSgn = np.concatenate((newVtcsSgn, tmpSgn)) - else: - newVtcsSgn = np.concatenate((newVtcsSgn, np.array([0]))) - raise RuntimeError("Inconsistency while creating seam 5, exiting.") else: idxPt1 = newVtcsAdj[t4c[i, idxPos1], t4c[i, idxNeg1]] @@ -303,14 +256,14 @@ def getSeamCase2(v4c, t4c, i4c, k4c, vfuncXEv1, e4cBndOpen, newTetra, newVtcs, n ) # return - return v4c, t4c, i4c, k4c, e4cBndOpen, newTetra, newVtcs, newVtcsAdj, newVtcsSgn + return v4c, t4c, i4c, k4c, e4cBndOpen, newTetra, newVtcs, newVtcsAdj # -------------------------------------------------------------------------- # getSeam subfunction #3 # case 3: zero nans in the tetra (could happen for boundary tetras # towards two sides) - def getSeamCase3(v4c, t4c, i4c, k4c, vfuncXEv1, e4cBndOpen, newTetra, newVtcs, newVtcsAdj, newVtcsSgn, i): + def getSeamCase3(v4c, t4c, i4c, k4c, vfuncXEv1, e4cBndOpen, newTetra, newVtcs, newVtcsAdj, i): # determine which edges of the current tetra belong to the # boundary surface numEdgeBnd = np.where( @@ -406,8 +359,8 @@ def getSeamCase3(v4c, t4c, i4c, k4c, vfuncXEv1, e4cBndOpen, newTetra, newVtcs, n ), ] ] = np.nan - v4c, t4c, i4c, k4c, e4cBndOpen, newTetra, newVtcs, newVtcsAdj, newVtcsSgn = getSeamCase2( - v4c, t4c, i4c, k4c, tmpvfuncXEv1, e4cBndOpen, newTetra, newVtcs, newVtcsAdj, newVtcsSgn, i + v4c, t4c, i4c, k4c, e4cBndOpen, newTetra, newVtcs, newVtcsAdj = getSeamCase2( + v4c, t4c, i4c, k4c, tmpvfuncXEv1, e4cBndOpen, newTetra, newVtcs, newVtcsAdj, i ) elif ( len(numEdgeBnd) == 3 @@ -432,8 +385,8 @@ def getSeamCase3(v4c, t4c, i4c, k4c, vfuncXEv1, e4cBndOpen, newTetra, newVtcs, n ), ] ] = np.nan - v4c, t4c, i4c, k4c, e4cBndOpen, newTetra, newVtcs, newVtcsAdj, newVtcsSgn = getSeamCase1( - v4c, t4c, i4c, k4c, tmpvfuncXEv1, e4cBndOpen, newTetra, newVtcs, newVtcsAdj, newVtcsSgn, i + v4c, t4c, i4c, k4c, e4cBndOpen, newTetra, newVtcs, newVtcsAdj = getSeamCase1( + v4c, t4c, i4c, k4c, tmpvfuncXEv1, e4cBndOpen, newTetra, newVtcs, newVtcsAdj, i ) elif ( len(numEdgeBnd) == 3 @@ -459,14 +412,14 @@ def getSeamCase3(v4c, t4c, i4c, k4c, vfuncXEv1, e4cBndOpen, newTetra, newVtcs, n ), ] ] = np.nan - v4c, t4c, i4c, k4c, e4cBndOpen, newTetra, newVtcs, newVtcsAdj, newVtcsSgn = getSeamCase1( - v4c, t4c, i4c, k4c, tmpvfuncXEv1, e4cBndOpen, newTetra, newVtcs, newVtcsAdj, newVtcsSgn, i + v4c, t4c, i4c, k4c, e4cBndOpen, newTetra, newVtcs, newVtcsAdj = getSeamCase1( + v4c, t4c, i4c, k4c, tmpvfuncXEv1, e4cBndOpen, newTetra, newVtcs, newVtcsAdj, i ) else: raise RuntimeError("Unknown getSeam() case, exiting.") # return - return v4c, t4c, i4c, k4c, e4cBndOpen, newTetra, newVtcs, newVtcsAdj, newVtcsSgn + return v4c, t4c, i4c, k4c, e4cBndOpen, newTetra, newVtcs, newVtcsAdj # -------------------------------------------------------------------------- # main part of function @@ -474,7 +427,6 @@ def getSeamCase3(v4c, t4c, i4c, k4c, vfuncXEv1, e4cBndOpen, newTetra, newVtcs, n # initialize newVtcs = np.empty(shape=(0), dtype=int) - newVtcsSgn = np.empty(shape=(0), dtype=int) newVtcsAdj = sp.lil_matrix((len(v4c), len(v4c)), dtype=int) newTetra = np.empty(shape=(0, 4), dtype=int) @@ -482,10 +434,6 @@ def getSeamCase3(v4c, t4c, i4c, k4c, vfuncXEv1, e4cBndOpen, newTetra, newVtcs, n vfuncXEv1 = np.zeros(np.shape(v4c)[0]) * np.nan vfuncXEv1[v4cBndOpenKeep] = anisoLaplEvec[:, 1] - # assign second eigenfunction - vfuncXEv2 = np.zeros(np.shape(v4c)[0]) * np.nan - vfuncXEv2[v4cBndOpenKeep] = np.sign(anisoLaplEvec[:, 2]) - # find tetra with zeros: first, transform boundary trias to edges; here, # the important trick is to start working with the boundary trias from # t4cBndOpen, which is to avoid trias that have e.g. 2 points at the upper @@ -521,8 +469,8 @@ def getSeamCase3(v4c, t4c, i4c, k4c, vfuncXEv1, e4cBndOpen, newTetra, newVtcs, n # if there is one nan, we need to make sure that the three vtcs form a surface triangle; # alternatively we check if it is at least a surface edge if np.isin(t4cBndOpen, t4c[i, np.logical_not(np.isnan(vfuncXEv1[t4c[i]]))]).all(axis=1).any(): - v4c, t4c, i4c, k4c, e4cBndOpen, newTetra, newVtcs, newVtcsAdj, newVtcsSgn = getSeamCase1( - v4c, t4c, i4c, k4c, vfuncXEv1, e4cBndOpen, newTetra, newVtcs, newVtcsAdj, newVtcsSgn, i + v4c, t4c, i4c, k4c, e4cBndOpen, newTetra, newVtcs, newVtcsAdj = getSeamCase1( + v4c, t4c, i4c, k4c, vfuncXEv1, e4cBndOpen, newTetra, newVtcs, newVtcsAdj, i ) elif np.isin(e4cBndOpen, t4c[i, np.logical_not(np.isnan(vfuncXEv1[t4c[i]]))]).all(axis=1).any(): # need to temporarily set one vtx to nan @@ -541,8 +489,8 @@ def getSeamCase3(v4c, t4c, i4c, k4c, vfuncXEv1, e4cBndOpen, newTetra, newVtcs, n ) tmpvfuncXEv1 = vfuncXEv1.copy() tmpvfuncXEv1[t4c[i, tmpEdgeNot]] = np.nan - v4c, t4c, i4c, k4c, e4cBndOpen, newTetra, newVtcs, newVtcsAdj, newVtcsSgn = getSeamCase2( - v4c, t4c, i4c, k4c, tmpvfuncXEv1, e4cBndOpen, newTetra, newVtcs, newVtcsAdj, newVtcsSgn, i + v4c, t4c, i4c, k4c, e4cBndOpen, newTetra, newVtcs, newVtcsAdj = getSeamCase2( + v4c, t4c, i4c, k4c, tmpvfuncXEv1, e4cBndOpen, newTetra, newVtcs, newVtcsAdj, i ) else: raise RuntimeError("Vtx " + str(i) + " not in triangle or edge, skipping") @@ -560,14 +508,14 @@ def getSeamCase3(v4c, t4c, i4c, k4c, vfuncXEv1, e4cBndOpen, newTetra, newVtcs, n ) if len(np.where(tmpEdge)) > 1: raise RuntimeError("Found " + str(len(np.where(tmpEdge))) + " edges, exiting.") - v4c, t4c, i4c, k4c, e4cBndOpen, newTetra, newVtcs, newVtcsAdj, newVtcsSgn = getSeamCase2( - v4c, t4c, i4c, k4c, vfuncXEv1, e4cBndOpen, newTetra, newVtcs, newVtcsAdj, newVtcsSgn, i + v4c, t4c, i4c, k4c, e4cBndOpen, newTetra, newVtcs, newVtcsAdj = getSeamCase2( + v4c, t4c, i4c, k4c, vfuncXEv1, e4cBndOpen, newTetra, newVtcs, newVtcsAdj, i ) elif np.sum(np.isnan(vfuncXEv1[t4c[i]])) == 0: # if there are no nans, the checking is done within the getSeamCase3 function - v4c, t4c, i4c, k4c, e4cBndOpen, newTetra, newVtcs, newVtcsAdj, newVtcsSgn = getSeamCase3( - v4c, t4c, i4c, k4c, vfuncXEv1, e4cBndOpen, newTetra, newVtcs, newVtcsAdj, newVtcsSgn, i + v4c, t4c, i4c, k4c, e4cBndOpen, newTetra, newVtcs, newVtcsAdj = getSeamCase3( + v4c, t4c, i4c, k4c, vfuncXEv1, e4cBndOpen, newTetra, newVtcs, newVtcsAdj, i ) else: @@ -581,7 +529,7 @@ def getSeamCase3(v4c, t4c, i4c, k4c, vfuncXEv1, e4cBndOpen, newTetra, newVtcs, n t4c = np.concatenate((t4c, newTetra)) # return - return v4c, t4c, i4c, k4c, newVtcs, newVtcsSgn, newTetra + return v4c, t4c, i4c, k4c, newVtcs, newTetra # ------------------------------------------------------------------------------ @@ -682,10 +630,15 @@ def computeCubeParam(params): # post-process eigenfunction (order, flipping) # decide whether or not to flip anisoLaplEvec[:, 1] and anisoLaplEvec[:, 2]: + + # for the hippocampus: # - anisoLaplEvec[:, 1] should be zero at the extrema of 234 and 240. # - anisoLaplEvec[:, 2] should have extremal values at 234 and 240. - # We therefore check if the zeros of Evec1 are located in 234 or 240/2420 - # rather than 236 or 238. If not, we change order. + # - We therefore check if the zeros of Evec1 are located in 234 or 240/2420 + # rather than 236 or 238. If not, we change order. + + # for the ERC: + # - ... # assign first eigenfunction vfuncXEv1 = np.zeros(np.shape(v4c)[0]) * np.nan @@ -693,141 +646,219 @@ def computeCubeParam(params): # get edges that contain zeros, and turn into list of vertices e4cBndOpen = np.concatenate((t4cBndOpen[:, (0, 1)], t4cBndOpen[:, (0, 2)], t4cBndOpen[:, (1, 2)]), axis=0) e4cBndOpen = np.unique(np.sort(e4cBndOpen, axis=1), axis=1) - # zeroEdgeIdx = e4cBndOpen[np.where(np.abs(np.sum(np.sign(vfuncXEv1[e4cBndOpen]), axis=1)) == 0)[0]] - # decide whether or not to flip anisoLaplEvec[:, 1] (should be inf -> sup) + if params.LUT == "freesurfer" or params.LUT == "freesurfer-no_ml" or params.LUT == "ashs-penn_abc_3t_t2" or params.LUT == "ashs-penn_abc_3t_t2_ext" or params.LUT == "ashs-umcutrecht_7t" or os.path.isfile(params.LUT): # noqa: E501 - # check EV1 for flips; this is done indirectly: we get the locations of - # highest and lowest values in 234 (PrSbc) or 236 (Sbc); highest values - # should have a higher value on the z-axis than lower values. + # if working on the hippocampus: + + # decide whether or not to flip anisoLaplEvec[:, 1] (should be inf -> sup) + + # check EV1 for flips; this is done indirectly: we get the locations of + # highest and lowest values in 234 (PrSbc) or 236 (Sbc); highest values + # should have a higher value on the z-axis than lower values. + # we now convert from the surface coords to ras coords, because we cannot + # assume the surface coord system to be always the same (rightly or not). + # we rely on appropriate transformation matrices, though. + + prsbc_or_sbc = np.logical_or(k4c[hsfList] == params.LUTDICT["presubiculum"], k4c[hsfList] == params.LUTDICT["subiculum"]) # noqa: E501 + + img = nb.load(os.path.join(params.OUTDIR, params.HEMI + ".image.mgz")) + mat = np.linalg.inv(img.header.get_vox2ras_tkr()) + pts = np.concatenate((v4cBndOpenRm, np.ones((v4cBndOpenRm.shape[0], 1))), axis=1) + v4cBndOpenRm_ras = np.matmul(np.matmul(pts, mat.transpose()), img.header.get_vox2ras().transpose()) + + if np.median(v4cBndOpenRm_ras[np.where(np.logical_and(anisoLaplEvec[:, 1] > 0, prsbc_or_sbc))[0], 2]) > np.median(v4cBndOpenRm_ras[np.where(np.logical_and(anisoLaplEvec[:, 1] < 0, prsbc_or_sbc))[0], 2]): # noqa: E501 + pass + elif np.median(v4cBndOpenRm_ras[np.where(np.logical_and(anisoLaplEvec[:, 1] > 0, prsbc_or_sbc))[0], 2]) < np.median(v4cBndOpenRm_ras[np.where(np.logical_and(anisoLaplEvec[:, 1] < 0, prsbc_or_sbc))[0], 2,]): # noqa: E501 + LOGGER.info("Flipping EV1") + anisoLaplEvec[:, 1] = -anisoLaplEvec[:, 1] + else: + io.write_vfunc(os.path.join(cutTetraMeshDir, params.HEMI + ".lapy.aLBO.EV1.psol"), anisoLaplEvec[:, 1]) + io.write_vfunc(os.path.join(cutTetraMeshDir, params.HEMI + ".lapy.aLBO.EV2.psol"), anisoLaplEvec[:, 2]) + nb.freesurfer.save( + nb.freesurfer.MGHImage(dataobj=anisoLaplEvec[:, 1].astype("float32"), affine=None), + filename=os.path.join(cutTetraMeshDir, params.HEMI + ".lapy.aLBO.EV1.mgh"), + ) + nb.freesurfer.save( + nb.freesurfer.MGHImage(dataobj=anisoLaplEvec[:, 2].astype("float32"), affine=None), + filename=os.path.join(cutTetraMeshDir, params.HEMI + ".lapy.aLBO.EV2.mgh"), + ) + raise RuntimeError("Inconsistency detected for EV1, exiting.") + + # decide whether or not to flip anisoLaplEvec[:, 2] (should be 234 -> 240) + + if (np.median(anisoLaplEvec[np.where(prsbc_or_sbc)[0], 2]) < 0): + pass + elif (np.median(anisoLaplEvec[np.where(prsbc_or_sbc)[0], 2]) > 0): + LOGGER.info("Flipping EV2") + anisoLaplEvec[:, 2] = -anisoLaplEvec[:, 2] + else: + io.write_vfunc(os.path.join(cutTetraMeshDir, params.HEMI + ".lapy.aLBO.EV1.psol"), anisoLaplEvec[:, 1]) + io.write_vfunc(os.path.join(cutTetraMeshDir, params.HEMI + ".lapy.aLBO.EV2.psol"), anisoLaplEvec[:, 2]) + nb.freesurfer.save( + nb.freesurfer.MGHImage(dataobj=anisoLaplEvec[:, 1].astype("float32"), affine=None), + filename=os.path.join(cutTetraMeshDir, params.HEMI + ".lapy.aLBO.EV1.mgh"), + ) + nb.freesurfer.save( + nb.freesurfer.MGHImage(dataobj=anisoLaplEvec[:, 2].astype("float32"), affine=None), + filename=os.path.join(cutTetraMeshDir, params.HEMI + ".lapy.aLBO.EV2.mgh"), + ) + raise RuntimeError("Inconsistency detected for EV2, exiting.") + + elif params.LUT == "ashs-penn_abc_3t_t2_ent": + + # if working on the ERC: + + # decide whether or not to flip anisoLaplEvec[:, 1] (should be inf -> sup). + # we now convert from the surface coords to ras coords, because we cannot + # assume the surface coord system to be always the same (rightly or not). + # we rely on appropriate transformation matrices, though. + + erc = k4c[hsfList] == params.LUTDICT["entorhinal"] + + img = nb.load(os.path.join(params.OUTDIR, params.HEMI + ".image.mgz")) + mat = np.linalg.inv(img.header.get_vox2ras_tkr()) + pts = np.concatenate((v4cBndOpenRm, np.ones((v4cBndOpenRm.shape[0], 1))), axis=1) + v4cBndOpenRm_ras = np.matmul(np.matmul(pts, mat.transpose()), img.header.get_vox2ras().transpose()) + + if np.median(v4cBndOpenRm_ras[np.where(np.logical_and(anisoLaplEvec[:, 1] > 0, erc))[0],2]) > np.median(v4cBndOpenRm_ras[np.where(np.logical_and(anisoLaplEvec[:, 1] < 0, erc))[0],2]): # noqa: E501 + pass + elif np.median(v4cBndOpenRm_ras[np.where(np.logical_and(anisoLaplEvec[:, 1] > 0, erc))[0],2]) < np.median(v4cBndOpenRm_ras[np.where(np.logical_and(anisoLaplEvec[:, 1] < 0, erc))[0],2]): # noqa: E501 + LOGGER.info("Flipping EV1") + anisoLaplEvec[:, 1] = -anisoLaplEvec[:, 1] + else: + io.write_vfunc(os.path.join(cutTetraMeshDir, params.HEMI + ".lapy.aLBO.EV1.psol"), anisoLaplEvec[:, 1]) + io.write_vfunc(os.path.join(cutTetraMeshDir, params.HEMI + ".lapy.aLBO.EV2.psol"), anisoLaplEvec[:, 2]) + nb.freesurfer.save( + nb.freesurfer.MGHImage(dataobj=anisoLaplEvec[:, 1].astype("float32"), affine=None), + filename=os.path.join(cutTetraMeshDir, params.HEMI + ".lapy.aLBO.EV1.mgh"), + ) + nb.freesurfer.save( + nb.freesurfer.MGHImage(dataobj=anisoLaplEvec[:, 2].astype("float32"), affine=None), + filename=os.path.join(cutTetraMeshDir, params.HEMI + ".lapy.aLBO.EV2.mgh"), + ) + raise RuntimeError("Inconsistency detected for EV1, exiting.") + + # decide whether or not to flip anisoLaplEvec[:, 2] (should be ... -> BA35) + + any_positive_ba35 = np.where(np.logical_and(anisoLaplEvec[:, 2] > 0, k4c[hsfList] == params.LUTDICT["bndba35"]))[0] # noqa: E501 + any_negative_ba35 = np.where(np.logical_and(anisoLaplEvec[:, 2] < 0, k4c[hsfList] == params.LUTDICT["bndba35"]))[0] # noqa: E501 + + if len(any_negative_ba35) == 0 and len(any_positive_ba35) > 0: + pass + elif len(any_positive_ba35) == 0 and len(any_negative_ba35) > 0: + LOGGER.info("Flipping EV2") + anisoLaplEvec[:, 2] = -anisoLaplEvec[:, 2] + elif len(any_positive_ba35) > len(any_negative_ba35): + logging.warning("Could not determine EV2 orientation, check output.") + pass + elif len(any_positive_ba35) <= len(any_negative_ba35): + LOGGER.info("Flipping EV2") + anisoLaplEvec[:, 2] = -anisoLaplEvec[:, 2] + logging.warning("Could not determine EV2 orientation, check output.") + else: + io.write_vfunc(os.path.join(cutTetraMeshDir, params.HEMI + ".lapy.aLBO.EV1.psol"), anisoLaplEvec[:, 1]) + io.write_vfunc(os.path.join(cutTetraMeshDir, params.HEMI + ".lapy.aLBO.EV2.psol"), anisoLaplEvec[:, 2]) + nb.freesurfer.save( + nb.freesurfer.MGHImage(dataobj=anisoLaplEvec[:, 1].astype("float32"), affine=None), + filename=os.path.join(cutTetraMeshDir, params.HEMI + ".lapy.aLBO.EV1.mgh"), + ) + nb.freesurfer.save( + nb.freesurfer.MGHImage(dataobj=anisoLaplEvec[:, 2].astype("float32"), affine=None), + filename=os.path.join(cutTetraMeshDir, params.HEMI + ".lapy.aLBO.EV2.mgh"), + ) + raise RuntimeError("Inconsistency detected for EV2, exiting.") - if np.median( - v4cBndOpenRm[ - np.where( - np.logical_and( - anisoLaplEvec[:, 1] > 0, - np.logical_or( - k4c[hsfList] == params.LUTDICT["presubiculum"], k4c[hsfList] == params.LUTDICT["subiculum"] - ), - ) - )[0], - 2, - ] - ) > np.median( - v4cBndOpenRm[ - np.where( - np.logical_and( - anisoLaplEvec[:, 1] < 0, - np.logical_or( - k4c[hsfList] == params.LUTDICT["presubiculum"], k4c[hsfList] == params.LUTDICT["subiculum"] - ), - ) - )[0], - 2, - ] - ): - pass - elif np.median( - v4cBndOpenRm[ - np.where( - np.logical_and( - anisoLaplEvec[:, 1] > 0, - np.logical_or( - k4c[hsfList] == params.LUTDICT["presubiculum"], k4c[hsfList] == params.LUTDICT["subiculum"] - ), - ) - )[0], - 2, - ] - ) < np.median( - v4cBndOpenRm[ - np.where( - np.logical_and( - anisoLaplEvec[:, 1] < 0, - np.logical_or( - k4c[hsfList] == params.LUTDICT["presubiculum"], k4c[hsfList] == params.LUTDICT["subiculum"] - ), - ) - )[0], - 2, - ] - ): - LOGGER.info("Flipping EV1") - anisoLaplEvec[:, 1] = -anisoLaplEvec[:, 1] - else: - io.write_vfunc(os.path.join(cutTetraMeshDir, params.HEMI + ".lapy.aLBO.EV1.psol"), anisoLaplEvec[:, 1]) - io.write_vfunc(os.path.join(cutTetraMeshDir, params.HEMI + ".lapy.aLBO.EV2.psol"), anisoLaplEvec[:, 2]) - raise RuntimeError("Inconsistency detected for EV1, exiting.") - - # decide whether or not to flip anisoLaplEvec[:, 2] (should be 234 -> 240) - - if ( - np.median( - anisoLaplEvec[ - np.where( - np.logical_or( - k4c[hsfList] == params.LUTDICT["presubiculum"], k4c[hsfList] == params.LUTDICT["subiculum"] - ) - )[0], - 2, - ] - ) - < 0 - ): - pass - elif ( - np.median( - anisoLaplEvec[ - np.where( - np.logical_or( - k4c[hsfList] == params.LUTDICT["presubiculum"], k4c[hsfList] == params.LUTDICT["subiculum"] - ) - )[0], - 2, - ] - ) - > 0 - ): - LOGGER.info("Flipping EV2") - anisoLaplEvec[:, 2] = -anisoLaplEvec[:, 2] else: - io.write_vfunc(os.path.join(cutTetraMeshDir, params.HEMI + ".lapy.aLBO.EV1.psol"), anisoLaplEvec[:, 1]) - io.write_vfunc(os.path.join(cutTetraMeshDir, params.HEMI + ".lapy.aLBO.EV2.psol"), anisoLaplEvec[:, 2]) - raise RuntimeError("Inconsistency detected for EV2, exiting.") + LOGGER.warning("Warning: could not determine flipping of eigenfunction. As a result, medial/lateral and/or anterior/posterior orientations may be inconsistent. Check your results.") # noqa: E501 io.write_vfunc(os.path.join(cutTetraMeshDir, params.HEMI + ".lapy.aLBO.EV1.psol"), anisoLaplEvec[:, 1]) io.write_vfunc(os.path.join(cutTetraMeshDir, params.HEMI + ".lapy.aLBO.EV2.psol"), anisoLaplEvec[:, 2]) + nb.freesurfer.save( + nb.freesurfer.MGHImage(dataobj=anisoLaplEvec[:, 1].astype("float32"), affine=None), + filename=os.path.join(cutTetraMeshDir, params.HEMI + ".lapy.aLBO.EV1.mgh"), + ) + nb.freesurfer.save( + nb.freesurfer.MGHImage(dataobj=anisoLaplEvec[:, 2].astype("float32"), affine=None), + filename=os.path.join(cutTetraMeshDir, params.HEMI + ".lapy.aLBO.EV2.mgh"), + ) + # -------------------------------------------------------------------------- # define boundary conditions (note that we changed dimensions and directions!) # get seam - v4c, t4c, i4c, k4c, newVtcs, newVtcsSgn, newTetra = getSeam( + v4c, t4c, i4c, k4c, newVtcs, newTetra = getSeam( v4c, t4c, i4c, k4c, v4cBndOpenKeep, t4cBndOpen, anisoLaplEvec ) - # get 2nd principal axis, note that direction may be arbitrary, will - # possibly be switched below; note that this replaces EV2-based newVtcsSgn, - # which we may remove in the future + # x: + + if params.LUT in ( + "freesurfer", + "freesurfer-no_ml", + "ashs-penn_abc_3t_t2", + "ashs-penn_abc_3t_t2_ent", + "ashs-umcutrecht_7t", + ) or os.path.isfile(params.LUT): + + # if working on the hippocampus without ctx or on the ERC: + + # x: medial -> lateral (will initially run along the first axis of the mesh; may be + # flipped during the left / right determination below) + # srxX will have -1 for all newVtcs where dim 0 coordinates are lower than their mean + # snkX will have +1 for all newVtcs where dim 0 coordinates are higher than their mean + + vfuncX = np.zeros(np.shape(v4c)[0]) + srcX = np.zeros(np.shape(v4c)[0]) + srcX[newVtcs[v4c[newVtcs, 0] < np.mean(v4c[newVtcs, 0])]] = -1 + snkX = np.zeros(np.shape(v4c)[0]) + snkX[newVtcs[v4c[newVtcs, 0] >= np.mean(v4c[newVtcs, 0])]] = 1 + vfuncX = srcX + snkX + + # determine left / right; here we map to ras (and don't rely on ras_tkr initially) + + img = nb.load(os.path.join(params.OUTDIR, params.HEMI + ".image.mgz")) + mat = np.linalg.inv(img.header.get_vox2ras_tkr()) + pts = np.concatenate((v4c, np.ones((v4c.shape[0], 1))), axis=1) + ras = np.matmul(np.matmul(pts, mat.transpose()), img.header.get_vox2ras().transpose()) + + if params.HEMI == "lh": + # first make sure that ras coords are indeed in left hemi + if (ras[:, 0] < 0.0).all(): + # lateral ras coords (snk, +1) should be more negative than medial ras coords (src, -1), otherwise flip + if np.mean(ras[np.where(vfuncX > 0)[0], 0]) > np.mean(ras[np.where(vfuncX < 0)[0], 0]): + vfuncX = -vfuncX + else: + raise RuntimeError("Ambiguous hemisphere info, exiting.") + elif params.HEMI == "rh": + # first make sure that ras coords are indeed in right hemi + if (ras[:, 0] > 0.0).all(): + # lateral ras coords (snk, +1) should be more positive than medial ras coords (src, -1), otherwise flip + if np.mean(ras[np.where(vfuncX > 0)[0], 0]) < np.mean(ras[np.where(vfuncX < 0)[0], 0]): + vfuncX = -vfuncX + else: + raise RuntimeError("Ambiguous hemisphere info, exiting.") + + elif params.LUT == "ashs-penn_abc_3t_t2_ext": + + # if working on the hippocampus with ctx: - # pca = PCA(n_components=3) - # pca_scores = pca.fit_transform(v4c) + # x: inf -> sup (phc -> ca3). cannot use medial -> lateral, because phc may bend below hc + # will transform surface coords to ras coords and evaluate dim 2 (3rd dim) - # x: medial -> lateral + img = nb.load(os.path.join(params.OUTDIR, params.HEMI + ".image.mgz")) + mat = np.linalg.inv(img.header.get_vox2ras_tkr()) + pts = np.concatenate((v4c[newVtcs, :], np.ones((v4c[newVtcs, :].shape[0], 1))), axis=1) + ras = np.matmul(np.matmul(pts, mat.transpose()), img.header.get_vox2ras().transpose()) - vfuncX = np.zeros(np.shape(v4c)[0]) - srcX = np.zeros(np.shape(v4c)[0]) - srcX[newVtcs[v4c[newVtcs, 0] < np.mean(v4c[newVtcs, 0])]] = -1 - ###srcX[newVtcs[pca_scores[newVtcs,1]<0]] = -1 - ###srcX[newVtcs[v4c[newVtcs,0]>0]] = 1 - ###srcX[newVtcs[newVtcsSgn>0]] = 1 - snkX = np.zeros(np.shape(v4c)[0]) - snkX[newVtcs[v4c[newVtcs, 0] >= np.mean(v4c[newVtcs, 0])]] = 1 - ###snkX[newVtcs[pca_scores[newVtcs,1]>0]] = 1 - ###snkX[newVtcs[v4c[newVtcs,0]<0]] = -1 - ###snkX[newVtcs[newVtcsSgn<0]] = -1 - vfuncX = srcX + snkX + vfuncX = np.zeros(np.shape(v4c)[0]) + srcX = np.zeros(np.shape(v4c)[0]) + srcX[newVtcs[ras[:, 2] < np.mean(ras[:, 2])]] = -1 + snkX = np.zeros(np.shape(v4c)[0]) + snkX[newVtcs[ras[:, 2] >= np.mean(ras[:, 2])]] = 1 + vfuncX = srcX + snkX # y: Tail (226) -> Head (232) (post --> ant) @@ -835,7 +866,7 @@ def computeCubeParam(params): vfuncY[np.where(i4c == params.LUTDICT["jointtail"])[0]] = -1 vfuncY[np.where(i4c == params.LUTDICT["jointhead"])[0]] = 1 - # z: inf->sup + # z: inf -> sup vfuncZ = np.zeros(np.shape(v4c)[0]) srcZ = np.zeros(np.shape(anisoLaplEvec[:, 1])[0]) @@ -844,28 +875,6 @@ def computeCubeParam(params): snkZ[anisoLaplEvec[:, 1] > 0] = 1 vfuncZ[v4cBndOpenKeep] = srcZ + snkZ - # determine left / right - - img = nb.load(os.path.join(params.OUTDIR, params.HEMI + ".image.mgz")) - mat = np.linalg.inv(img.header.get_vox2ras_tkr()) - pts = np.concatenate((v4c, np.ones((v4c.shape[0], 1))), axis=1) - ras = np.matmul(np.matmul(pts, mat.transpose()), img.header.get_vox2ras().transpose()) - - if params.HEMI == "lh": - if (ras[:, 0] < 0.0).all(): - # lateral (snk, +1) should be more negative than medial (src, -1), otherwise flip - if np.mean(ras[np.where(vfuncX > 0)[0], 0]) > np.mean(ras[np.where(vfuncX < 0)[0], 0]): - vfuncX = -vfuncX - else: - raise RuntimeError("Ambiguous hemisphere info, exiting.") - elif params.HEMI == "rh": - if (ras[:, 0] > 0.0).all(): - # lateral (snk, +1) should be more positive than medial (src, -1), otherwise flip - if np.mean(ras[np.where(vfuncX > 0)[0], 0]) < np.mean(ras[np.where(vfuncX < 0)[0], 0]): - vfuncX = -vfuncX - else: - raise RuntimeError("Ambiguous hemisphere info, exiting.") - # -------------------------------------------------------------------------- # compute ABtetra and Poisson functions @@ -924,6 +933,8 @@ def computeCubeParam(params): # v4cRmBndRm = tetMesh4cRmBnd.v # t4cRmBndRm = tetMesh4cRmBnd.t + tetMesh4cRmBnd.orient_() + TriaMesh.write_vtk(tetMesh4cRmBnd, filename=os.path.join(cutTetraMeshDir, params.HEMI + ".rm.bnd.seam.rm.cut.vtk")) io.write_vfunc( @@ -940,6 +951,19 @@ def computeCubeParam(params): io.write_vfunc(os.path.join(cutTetraMeshDir, params.HEMI + ".poisson1.rm.bnd.seam.rm.cut.psol"), P1[v4cRmBndRmKeep]) io.write_vfunc(os.path.join(cutTetraMeshDir, params.HEMI + ".poisson2.rm.bnd.seam.rm.cut.psol"), P2[v4cRmBndRmKeep]) + nb.freesurfer.save( + nb.freesurfer.MGHImage(dataobj=P0[v4cRmBndRmKeep].astype("float32"), affine=None), + filename=os.path.join(cutTetraMeshDir, params.HEMI + ".poisson0.rm.bnd.seam.rm.cut.mgh"), + ) + nb.freesurfer.save( + nb.freesurfer.MGHImage(dataobj=P1[v4cRmBndRmKeep].astype("float32"), affine=None), + filename=os.path.join(cutTetraMeshDir, params.HEMI + ".poisson1.rm.bnd.seam.rm.cut.mgh"), + ) + nb.freesurfer.save( + nb.freesurfer.MGHImage(dataobj=P2[v4cRmBndRmKeep].astype("float32"), affine=None), + filename=os.path.join(cutTetraMeshDir, params.HEMI + ".poisson2.rm.bnd.seam.rm.cut.mgh"), + ) + # write out cube w4cRm = np.vstack((P0, P1, P2)).transpose() diff --git a/hipsta/computeThickness.py b/hipsta/computeThickness.py index 7edda95..834589e 100755 --- a/hipsta/computeThickness.py +++ b/hipsta/computeThickness.py @@ -275,10 +275,10 @@ def computeThickness(params): dfx = pd.DataFrame(llxLgth) - dz = zip(range(0, dfx.shape[1]), ["z" + str(z) for z in range(0, dfx.shape[1])]) + dz = zip(range(0, dfx.shape[1]), ["z" + str(z) for z in range(0, dfx.shape[1])], strict=True) dfx = dfx.rename(columns=dict(dz)) - dy = zip(range(0, dfx.shape[0]), ["y" + str(z) for z in range(0, dfx.shape[0])]) + dy = zip(range(0, dfx.shape[0]), ["y" + str(z) for z in range(0, dfx.shape[0])], strict=True) dfx = dfx.rename(index=dict(dy)) dfx.to_csv(os.path.join(OUT_DIR, params.HEMI + ".grid-segments-x.csv")) @@ -321,10 +321,10 @@ def computeThickness(params): dfy = pd.DataFrame(llyLgth) - dz = zip(range(0, dfy.shape[1]), ["z" + str(z) for z in range(0, dfy.shape[1])]) + dz = zip(range(0, dfy.shape[1]), ["z" + str(z) for z in range(0, dfy.shape[1])], strict=True) dfy = dfy.rename(columns=dict(dz)) - dx = zip(range(0, dfy.shape[0]), ["x" + str(z) for z in range(0, dfy.shape[0])]) + dx = zip(range(0, dfy.shape[0]), ["x" + str(z) for z in range(0, dfy.shape[0])], strict=True) dfy = dfy.rename(index=dict(dx)) dfy.to_csv(os.path.join(OUT_DIR, params.HEMI + ".grid-segments-y.csv")) @@ -367,10 +367,10 @@ def computeThickness(params): dfz = pd.DataFrame(llzLgth) - dy = zip(range(0, dfz.shape[1]), ["y" + str(z) for z in range(0, dfz.shape[1])]) + dy = zip(range(0, dfz.shape[1]), ["y" + str(z) for z in range(0, dfz.shape[1])], strict=True) dfz = dfz.rename(columns=dict(dy)) - dx = zip(range(0, dfz.shape[0]), ["x" + str(z) for z in range(0, dfz.shape[0])]) + dx = zip(range(0, dfz.shape[0]), ["x" + str(z) for z in range(0, dfz.shape[0])], strict=True) dfz = dfz.rename(index=dict(dx)) dfz.to_csv(os.path.join(OUT_DIR, params.HEMI + ".grid-segments-z.csv")) diff --git a/hipsta/createTetraLabels.py b/hipsta/createTetraLabels.py index 0d4d0ce..a5f11c2 100755 --- a/hipsta/createTetraLabels.py +++ b/hipsta/createTetraLabels.py @@ -96,31 +96,58 @@ def createTetraLabels(params): label=[params.LUTDICT["tail"], params.LUTDICT["head"]], bndlabel=[params.LUTDICT["bndtail"], params.LUTDICT["bndhead"]], ) - _createBoundaryMask( - hbtFile=HBTFILE, - mskFile=MSKFILE, - outFile=os.path.join(params.OUTDIR, "tetra-labels", params.HEMI + "." + "labels-ca4.mgz"), - label=[params.LUTDICT["ca4"]], - bndlabel=[params.LUTDICT["bndca4"]], - ) - # merge bnd and ca4 masks (add ca4 to bnd) - - cmd = ( - os.path.join(os.environ.get("FREESURFER_HOME"), "bin", "mris_calc") - + " " - + "--output " - + os.path.join(params.OUTDIR, "tetra-labels", params.HEMI + "." + "labels-bnd.mgz") - + " " - + os.path.join(params.OUTDIR, "tetra-labels", params.HEMI + "." + "labels-bnd.mgz") - + " " - + "lowerlimit " - + os.path.join(params.OUTDIR, "tetra-labels", params.HEMI + "." + "labels-ca4.mgz") - ) + if params.LUTDICT["ca4"] in params.HSFLIST: + _createBoundaryMask( + hbtFile=HBTFILE, + mskFile=MSKFILE, + outFile=os.path.join(params.OUTDIR, "tetra-labels", params.HEMI + "." + "labels-ca4.mgz"), + label=[params.LUTDICT["ca4"]], + bndlabel=[params.LUTDICT["bndca4"]], + ) + + # merge bnd and ca4 masks (add ca4 to bnd) + cmd = ( + os.path.join(os.environ.get("FREESURFER_HOME"), "bin", "mris_calc") + + " " + + "--output " + + os.path.join(params.OUTDIR, "tetra-labels", params.HEMI + "." + "labels-bnd.mgz") + + " " + + os.path.join(params.OUTDIR, "tetra-labels", params.HEMI + "." + "labels-bnd.mgz") + + " " + + "lowerlimit " + + os.path.join(params.OUTDIR, "tetra-labels", params.HEMI + "." + "labels-ca4.mgz") + ) + + print(cmd) + + subprocess.run(cmd.split(), capture_output=True) + + if params.LUT == "ashs-penn_abc_3t_t2_ent": + _createBoundaryMask( + hbtFile=HBTFILE, + mskFile=MSKFILE, + outFile=os.path.join(params.OUTDIR, "tetra-labels", params.HEMI + "." + "labels-ba35.mgz"), + label=[params.LUTDICT["ba35"]], + bndlabel=[params.LUTDICT["bndba35"]], + ) + + # merge bnd and ba35 masks (add ba35 to bnd) + cmd = ( + os.path.join(os.environ.get("FREESURFER_HOME"), "bin", "mris_calc") + + " " + + "--output " + + os.path.join(params.OUTDIR, "tetra-labels", params.HEMI + "." + "labels-bnd.mgz") + + " " + + os.path.join(params.OUTDIR, "tetra-labels", params.HEMI + "." + "labels-bnd.mgz") + + " " + + "lowerlimit " + + os.path.join(params.OUTDIR, "tetra-labels", params.HEMI + "." + "labels-ba35.mgz") + ) - print(cmd) + print(cmd) - subprocess.run(cmd.split(), capture_output=True) + subprocess.run(cmd.split(), capture_output=True) # load bnd image @@ -163,7 +190,7 @@ def createTetraLabels(params): nb.freesurfer.MGHImage( dataobj=np.array([float(x) for x in vtxLabels])[vRmBndKeep].astype("float32"), affine=None ), - filename=os.path.join(params.OUTDIR, "tetra-labels", params.HEMI + ".rm.bnd.label.mgh"), + filename=os.path.join(params.OUTDIR, "tetra-labels", params.HEMI + ".rm.bnd.tetra.mgh"), ) # return diff --git a/hipsta/cutTetra.py b/hipsta/cutTetra.py index 3fbd10b..625e89e 100755 --- a/hipsta/cutTetra.py +++ b/hipsta/cutTetra.py @@ -346,7 +346,7 @@ def cutTetra(params): ) A = sp.eye(np.shape(v4c)[0]) A = A + sp.coo_matrix( - (np.ones((np.shape(idx)[0])), (idx[:, 0], idx[:, 1])), shape=(np.shape(v4c)[0], np.shape(v4c)[0]) + (np.ones(np.shape(idx)[0]), (idx[:, 0], idx[:, 1])), shape=(np.shape(v4c)[0], np.shape(v4c)[0]) ) A = A + A.transpose() A[A > 0] = 1 diff --git a/hipsta/doc/DOCUMENTATION.md b/hipsta/doc/DOCUMENTATION.md index e72ce8a..4cf5940 100644 --- a/hipsta/doc/DOCUMENTATION.md +++ b/hipsta/doc/DOCUMENTATION.md @@ -8,6 +8,7 @@ thickness analysis. This page provides usage information as well as a technical description of processing steps and outputs. An explanative description can be found in the tutorial. + ## Usage: ### As a command-line script: @@ -31,12 +32,15 @@ Required arguments: --filename Filename of a segmentation file. --hemi Hemisphere. Either 'lh' or 'rh'. - --lut - Look-up table: a text file with numeric and verbal segmentation labels. 'freesurfer', 'ashs-penn_abc_3t_t2', and 'ashs-umcutrecht_7t' are keywords for built-in tables. + --lut + Look-up table: a text file with numeric and verbal segmentation labels. 'freesurfer', 'ashs-penn_abc_3t_t2', 'ashs-penn_abc_3t_t2_ext', 'ashs-penn_abc_3t_t2_ent', and 'ashs-umcutrecht_7t' are keywords for built-in tables. --outputdir Directory where the results will be written. Optional arguments: + --start-with-edited-labels + Start with edited labels. Requires an edited lh.labels.mgz or rh.labels.mgz file in an existing output directory. Existing files will be overwritten. + Start with edited masks. Requires an edited lh.mask.mgz or rh.mask.mgz file in an existing output directory. Existing files will be overwritten. --no-cleanup Do not remove files that may be useful for diagnostic or debugging purposes, but are not necessary otherwise. --no-crop Do not crop image. --upsample Upsample to the smallest voxel edge length. @@ -130,6 +134,7 @@ This script performs the following major processing steps: 10. map subfield values (and other volume-based data, optional) 11. create supplementary files for visualization + ## Outputs: 1. Thickness and curvature estimates @@ -223,14 +228,22 @@ Also the intermediate volume and surface files can be useful for quality control ## Supported segmentations: -- A hippocampal subfields segmentation based on FreeSurfer 7.1.1 or later should +- A hippocampal subfields segmentation based on FreeSurfer 7.1.1 or later should work as-is. Use `--lut freesurfer`. -- Two ASHS atlases are currently supported, the Penn ABC-3T ASHS Atlas for - T2-weighted MRI and the UMC Utrecht 7T atlas. Use `--lut ashs-penn_abc_3t_t2` +- Two ASHS atlases are currently supported, the Penn ABC-3T ASHS atlas for + T2-weighted MRI and the UMC Utrecht 7T atlas. Use `--lut ashs-penn_abc_3t_t2` or `--lut ashs-umcutrecht_7t`. If additional labels for the hippocampal head - (label value: 20) and tail (label value: 5) labels are missing, use the `--automask-head` - and/or `--automask-tail` arguments. It may also make sense to upsample the + (label value: 20) and tail (label value: 5) labels are missing, use the `--automask-head` + and/or `--automask-tail` arguments. It may also make sense to upsample the segmentations to isotropic voxel size using `--upsample`. +- Processing variants of the Penn ABC-3T ASHS atlas can be called with the + `--lut ashs-penn_abc_3t_t2_ext` and `--lut ashs-penn_abc_3t_t2_ent` arguments. + The first one will add (parts of) parahippocampal and entorhinal cortex as + well as BA35 to the analysis, but still restrict the analysis to the + hippocampal body region. The second one will perform a separate analysis of + the entorhinal cortex. Due to more complex or variable anatomy, these approaches + should to some extent be considered experimental and may require additional QC + or manual edits. - If using a custom look-up table (`--lut `), the expected format for the file is: ` `. `R`, `G`, `B`, `A` are numerical values for RGB colors and transparency (alpha) and will be ignored. @@ -256,7 +269,7 @@ Also the intermediate volume and surface files can be useful for quality control or the ASHS software. A custom segmentation is also permissible (some restrictions and settings apply; see [Supported Segmentations](#supported-segmentations)). -3. Python 3.8 or higher including the lapy, numpy, scipy, nibabel, pyvista, and +3. Python 3.10 or higher including the lapy, numpy, scipy, nibabel, pyvista, and pyacvd libraries. See `requirements.txt` for a full list. 4. The gmsh package (version 2.x; http://gmsh.info) must be installed. Can be diff --git a/hipsta/doc/documentation.py b/hipsta/doc/documentation.py index 24a11a8..3ec96c6 100644 --- a/hipsta/doc/documentation.py +++ b/hipsta/doc/documentation.py @@ -1,5 +1,5 @@ """ -This module provides documentation for the hippocampal shape and +This module provides documentation for the hippocampal shape and thickness analysis package. """ diff --git a/hipsta/hipsta.py b/hipsta/hipsta.py index 7366b80..3d3d46e 100644 --- a/hipsta/hipsta.py +++ b/hipsta/hipsta.py @@ -144,10 +144,11 @@ def _parse_arguments(): required.add_argument( "--lut", dest="lut", - help="Look-up table: a text file with numeric and verbal segmentation labels. 'freesurfer', \ - 'ashs-penn_abc_3t_t2', and 'ashs-umcutrecht_7t' are keywords for built-in tables.", + help="Look-up table: a text file with numeric and verbal segmentation labels. \'freesurfer\', \ + 'ashs-penn_abc_3t_t2', 'ashs-penn_abc_3t_t2_ent', 'ashs-penn_abc_3t_t2_ext', and \ + 'ashs-umcutrecht_7t\' are keywords for built-in tables.", default=None, - metavar="", + metavar="", required=False, ) required.add_argument( @@ -162,6 +163,26 @@ def _parse_arguments(): # optional arguments optional = parser.add_argument_group("Optional arguments") + optional.add_argument( + "--start-with-edited-labels", + dest="start_with_edited_labels", + help="Start with edited labels. Requires an edited lh.labels.mgz or rh.labels.mgz \ + file in an existing output directory. Existing files will be overwritten.", + default=get_defaults("start_with_edited_labels"), + action="store_true", + required=False, + ) + + optional.add_argument( + "--start-with-edited-masks", + dest="start_with_edited_masks", + help="Start with edited masks. Requires an edited lh.mask.mgz or rh.mask.mgz \ + file in an existing output directory. Existing files will be overwritten.", + default=get_defaults("start_with_edited_masks"), + action="store_true", + required=False, + ) + optional.add_argument( "--no-cleanup", dest="no_cleanup", @@ -492,6 +513,11 @@ def _evaluate_args(args): class settings: pass + # processing + + settings.START_WITH_EDITED_LABELS = args.start_with_edited_labels + settings.START_WITH_EDITED_MASKS = args.start_with_edited_masks + # clean up settings.CLEANUP = not (args.no_cleanup) @@ -598,11 +624,21 @@ class settings: class params: pass - params.FILENAME = os.path.abspath(args.filename) params.HEMI = args.hemi params.LUT = args.lut params.OUTDIR = args.outputdir + if settings.START_WITH_EDITED_LABELS is False and settings.START_WITH_EDITED_MASKS is False: + params.FILENAME = os.path.abspath(args.filename) + elif settings.START_WITH_EDITED_LABELS is False and settings.START_WITH_EDITED_MASKS is True: + params.FILENAME = os.path.join(params.OUTDIR, params.HEMI + ".mask.mgz") + elif settings.START_WITH_EDITED_LABELS is True and settings.START_WITH_EDITED_MASKS is False: + params.FILENAME = os.path.join(params.OUTDIR, params.HEMI + ".labels.mgz") + elif settings.START_WITH_EDITED_LABELS is True and settings.START_WITH_EDITED_MASKS is True: + raise RuntimeError("Start with edited labels and Start with edited masks cannot both be true.") + else: + raise RuntimeError("Error parsing input arguments (Start with edited labels / Start with edited masks).") + params.LUTDICT = LUTDICT params.HSFLIST = HSFLIST @@ -670,13 +706,19 @@ def _check_params(params): if ( params.LUT != "freesurfer" and params.LUT != "ashs-penn_abc_3t_t2" + and params.LUT != "ashs-penn_abc_3t_t2_ent" + and params.LUT != "ashs-penn_abc_3t_t2_ext" and params.LUT != "ashs-umcutrecht_7t" and not os.path.isfile(params.LUT) ): raise RuntimeError( - "Look-up table can only be 'freesurfer', 'ashs-penn_abc_3t_t2', 'ashs-umcutrecht_7t', \ - or an existing file, but not " - + params.LUT + "Look-up table can only be " + \ + "'freesurfer', " + \ + "'ashs-penn_abc_3t_t2', " + \ + "'ashs-penn_abc_3t_t2_ent', " + \ + "'ashs-penn_abc_3t_t2_ext', " + \ + "'ashs-umcutrecht_7t', " + \ + "or an existing file, but not " + params.LUT ) # check HSFLIST @@ -700,48 +742,51 @@ def _run_analysis(params): # process image (1) - LOGGER.info("Starting convertFormat() ...") - params = convertFormat(params) + if params.internal.START_WITH_EDITED_LABELS is False and params.internal.START_WITH_EDITED_MASKS is False: + LOGGER.info("Starting convertFormat() ...") + params = convertFormat(params) - LOGGER.info("Starting cropImage() ...") - params = cropImage(params) + LOGGER.info("Starting cropImage() ...") + params = cropImage(params) - LOGGER.info("Starting upsampleImage() ...") - params = upsampleImage(params) + LOGGER.info("Starting upsampleImage() ...") + params = upsampleImage(params) - LOGGER.info("Starting copy_image_to_main() ...") - params = copy_image_to_main(params) + LOGGER.info("Starting copy_image_to_main() ...") + params = copy_image_to_main(params) # process labels (2) - LOGGER.info("Starting autoMask() ...") - params = autoMask(params) + if params.internal.START_WITH_EDITED_LABELS is False and params.internal.START_WITH_EDITED_MASKS is False: + LOGGER.info("Starting autoMask() ...") + params = autoMask(params) - LOGGER.info("Starting createLabels() ...") - params = createLabels(params) + LOGGER.info("Starting createLabels() ...") + params = createLabels(params) - LOGGER.info("Starting mergeMolecularLayer() ...") - params = mergeMolecularLayer(params) + LOGGER.info("Starting mergeMolecularLayer() ...") + params = mergeMolecularLayer(params) - LOGGER.info("Starting copy_labels_to_main() ...") - params = copy_labels_to_main(params) + LOGGER.info("Starting copy_labels_to_main() ...") + params = copy_labels_to_main(params) # process mask (3) - LOGGER.info("Starting binarizeMask() ...") - params = binarizeMask(params) + if params.internal.START_WITH_EDITED_MASKS is False: + LOGGER.info("Starting binarizeMask() ...") + params = binarizeMask(params) - LOGGER.info("Starting gaussFilter() ...") - params = gaussFilter(params) + LOGGER.info("Starting gaussFilter() ...") + params = gaussFilter(params) - LOGGER.info("Starting longFilter() ...") - params = longFilter(params) + LOGGER.info("Starting longFilter() ...") + params = longFilter(params) - LOGGER.info("Starting closeMask() ...") - params = closeMask(params) + LOGGER.info("Starting closeMask() ...") + params = closeMask(params) - LOGGER.info("Starting copy_mask_to_main() ...") - params = copy_mask_to_main(params) + LOGGER.info("Starting copy_mask_to_main() ...") + params = copy_mask_to_main(params) # create surface (4) @@ -875,8 +920,9 @@ def run_hipsta(filename, hemi, lut, outputdir, **kwargs): hemi : Hemisphere. Either \'lh\' or \'rh\'. lut : - Look-up table: a text file with numeric and verbal segmentation labels. \'freesurfer\', 'ashs-penn_abc_3t_t2' - and \'ashs-umcutrecht_7t\' are keywords for built-in tables. + Look-up table: a text file with numeric and verbal segmentation labels. \'freesurfer\', + 'ashs-penn_abc_3t_t2', 'ashs-penn_abc_3t_t2_ent', 'ashs-penn_abc_3t_t2_ext', and \'ashs-umcutrecht_7t\' + are keywords for built-in tables. outputdir : Directory where the results will be written. diff --git a/hipsta/processLabels.py b/hipsta/processLabels.py index 1758ae7..67e32bb 100755 --- a/hipsta/processLabels.py +++ b/hipsta/processLabels.py @@ -56,22 +56,22 @@ def autoMask(params): labelCA1 = params.LUTDICT["ca1"] labelCA2 = params.LUTDICT["ca2"] labelCA3 = params.LUTDICT["ca3"] - # if "parahippocampal" in params.LUTDICT.keys(): - # labelPHC = params.LUTDICT["parahippocampal"] - # else: - # labelPHC = None - # if "entorhinal" in params.LUTDICT.keys(): - # labelENT = params.LUTDICT["entorhinal"] - # else: - # labelENT = None - # if "ba35" in params.LUTDICT.keys(): - # labelBA35 = params.LUTDICT["ba35"] - # else: - # labelBA35 = None + if "parahippocampal" in params.LUTDICT.keys(): + labelPHC = params.LUTDICT["parahippocampal"] + else: + labelPHC = None + if "entorhinal" in params.LUTDICT.keys(): + labelENT = params.LUTDICT["entorhinal"] + else: + labelENT = None + if "ba35" in params.LUTDICT.keys(): + labelBA35 = params.LUTDICT["ba35"] + else: + labelBA35 = None # if "ba36" in params.LUTDICT.keys(): # labelBA36 = params.LUTDICT["ba36"] # else: - # labelBA36 = None + # labelBA36 = None # head if params.internal.AUTOMASK_HEAD is True: @@ -80,11 +80,26 @@ def autoMask(params): # remove any potentially existing head labels dat[dat == labelHead] = 0 # - if labelCA2 is not None and labelCA3 is not None: - # find instance of CA2/CA3 - idxHead = np.argwhere(np.logical_or(dat == labelCA2, dat == labelCA3)) + if params.LUT == "ashs-penn_abc_3t_t2_ext": + if labelCA2 is not None and labelCA3 is not None and labelENT is not None and labelBA35 is not None: + # # find instance of CA2/CA3 or ENT or BA35 + # idxHead = np.argwhere(np.logical_or(dat == labelBA35, np.logical_or(dat == labelENT, np.logical_or(dat == labelCA2, dat == labelCA3)))) # noqa: E501 + # find instance of CA2/CA3 + idxHead = np.argwhere(np.logical_or(dat == labelCA2, dat == labelCA3)) + else: + raise RuntimeError("Insufficient label information, exiting.") + elif params.LUT == "ashs-penn_abc_3t_t2_ent": + if labelENT is not None: + # find instance of ENT + idxHead = np.argwhere(dat == labelENT) + else: + raise RuntimeError("Insufficient label information, exiting.") else: - raise RuntimeError("Insufficient label information, exiting.") + if labelCA2 is not None and labelCA3 is not None: + # find instance of CA2/CA3 + idxHead = np.argwhere(np.logical_or(dat == labelCA2, dat == labelCA3)) + else: + raise RuntimeError("Insufficient label information, exiting.") # get min/max row/col/slice if imgDimsAPDir == -1: cutFrom = np.min(idxHead[:, imgDimsAP]) + params.internal.AUTOMASK_HEAD_MARGIN @@ -118,23 +133,45 @@ def autoMask(params): # remove any potentially existing tail labels dat[dat == labelTail] = 0 # - if labelSbc is not None and labelCA1 is not None and labelCA2 is not None and labelCA3 is not None: - ## find instance of Sbc/CA1/CA2/CA3 - # idxTail = np.intersect1d( - # np.argwhere(dat == labelSbc)[:, imgDimsAP], - # np.intersect1d( - # np.argwhere(dat == labelCA1)[:, imgDimsAP], - # np.intersect1d( - # np.argwhere(dat == labelCA2)[:, imgDimsAP], np.argwhere(dat == labelCA3)[:, imgDimsAP] - # ), - # ), - # ) - # find instance of Sbc/CA1 - idxTail = np.intersect1d( - np.argwhere(dat == labelSbc)[:, imgDimsAP], np.argwhere(dat == labelCA1)[:, imgDimsAP] - ) + if params.LUT == "ashs-penn_abc_3t_t2_ext": + if labelPHC is not None and labelSbc is not None and labelCA1 is not None and labelCA2 is not None and labelCA3 is not None: # noqa: E501 + # find instance of Sbc/CA1/CA2/CA3 and PHC + idxTail = np.intersect1d(np.argwhere(dat == labelPHC)[:, imgDimsAP], + np.intersect1d(np.argwhere(dat == labelSbc)[:, imgDimsAP], + np.intersect1d(np.argwhere(dat == labelCA1)[:, imgDimsAP], + np.intersect1d(np.argwhere(dat == labelCA2)[:, imgDimsAP], + np.argwhere(dat == labelCA3)[:, imgDimsAP])))) + # # find instance of Sbc/CA1/CA2/CA3 + # idxTail = np.intersect1d(np.argwhere(dat == labelSbc)[:,imgDimsAP], + # np.intersect1d(np.argwhere(dat == labelCA1)[:,imgDimsAP], + # np.intersect1d(np.argwhere(dat == labelCA2)[:,imgDimsAP], + # np.argwhere(dat == labelCA3)[:,imgDimsAP]))) + else: + raise RuntimeError("Insufficient label information, exiting.") + elif params.LUT == "ashs-penn_abc_3t_t2_ent": + if labelENT is not None: + # find instance of ENT + idxTail = np.argwhere(dat == labelENT)[:, imgDimsAP] + else: + raise RuntimeError("Insufficient label information, exiting.") else: - raise RuntimeError("Insufficient label information, exiting.") + if labelSbc is not None and labelCA1 is not None and labelCA2 is not None and labelCA3 is not None: + ## find instance of Sbc/CA1/CA2/CA3 + # idxTail = np.intersect1d( + # np.argwhere(dat == labelSbc)[:, imgDimsAP], + # np.intersect1d( + # np.argwhere(dat == labelCA1)[:, imgDimsAP], + # np.intersect1d( + # np.argwhere(dat == labelCA2)[:, imgDimsAP], np.argwhere(dat == labelCA3)[:, imgDimsAP] + # ), + # ), + # ) + # find instance of Sbc/CA1 + idxTail = np.intersect1d( + np.argwhere(dat == labelSbc)[:, imgDimsAP], np.argwhere(dat == labelCA1)[:, imgDimsAP] + ) + else: + raise RuntimeError("Insufficient label information, exiting.") # get min/max row/col/slice if imgDimsAPDir == -1: cutFrom = np.max(idxTail) - params.internal.AUTOMASK_TAIL_MARGIN diff --git a/hipsta/removeBoundaryMask.py b/hipsta/removeBoundaryMask.py index ac7bfab..da55fd2 100755 --- a/hipsta/removeBoundaryMask.py +++ b/hipsta/removeBoundaryMask.py @@ -5,6 +5,7 @@ import os +import nibabel as nb import numpy as np from lapy import TetMesh, io @@ -162,6 +163,13 @@ def removeBoundaryMask(params): io.write_vfunc(os.path.join(params.OUTDIR, "tetra-cut", params.HEMI + ".rm.bnd.tetra.psol"), vIdx) + nb.freesurfer.save( + nb.freesurfer.MGHImage( + dataobj=vIdx[vRmBndKeep].astype("float32"), affine=None + ), + filename=os.path.join(params.OUTDIR, "tetra-cut", params.HEMI + ".rm.bnd.tetra.mgh"), + ) + # ------------------------------------------------------------------------- # return diff --git a/hipsta/utils/create_supplementary_files.py b/hipsta/utils/create_supplementary_files.py index 14ec868..1ea3941 100644 --- a/hipsta/utils/create_supplementary_files.py +++ b/hipsta/utils/create_supplementary_files.py @@ -86,34 +86,58 @@ def createSupplementaryFiles(params): lstBnd = [21.5, 22.5, 23.5] + elif LUT == "ashs-penn_abc_3t_t2_ext": + # get labels + + hsfBnd = nb.load(os.path.join(OUT_DIR, HEMI + ".mid-surface.hsf.mgh")) + hsfBnd = np.array(hsfBnd.get_fdata()).flatten() + + # recode to consecutive values + + hsfBnd[hsfBnd == 12] = 31 + hsfBnd[hsfBnd == 10] = 32 + hsfBnd[hsfBnd == 9] = 33 + hsfBnd[hsfBnd == 8] = 34 + hsfBnd[hsfBnd == 1] = 35 + hsfBnd[hsfBnd == 2] = 36 + hsfBnd[hsfBnd == 4] = 37 + + lstBnd = [31.5, 32.5, 33.5, 34.5, 35.5, 36.5] + else: logging.info( - "Creation of subfield boundary overlays is only supported for the freesurfer and ashs look-up tables." + "Creation of subfield boundary overlays is only supported for the freesurfer and (some) ashs look-up tables." # noqa: E501 ) - raise ValueError() - # get mid-surface + if LUT == "freesurfer" or LUT == "ashs-penn_abc_3t_t2" or LUT == "ashs-penn_abc_3t_t2_ext" or LUT == "ashs-umcutrecht_7t": # noqa: E501 - triaMidRm = TriaMesh.read_vtk(os.path.join(OUT_DIR, HEMI + ".mid-surface.vtk")) + # get mid-surface - # + triaMidRm = TriaMesh.read_vtk(os.path.join(OUT_DIR, HEMI + ".mid-surface.vtk")) - vBnd = np.empty((0, 3)) - tBnd = np.empty((0, 3)) + # - for i in lstBnd: - lvlBnd = levelsetsTria(triaMidRm.v, triaMidRm.t, hsfBnd, i) - if len(lvlBnd[0][0]) > 0: - vBnd_lvl = np.array(lvlBnd[0][0]) - tBnd_lvl = np.array(lvlBnd[1][0])[:, [0, 1, 1]] - 1 - tBnd = np.concatenate((tBnd, tBnd_lvl + len(vBnd)), axis=0) - vBnd = np.concatenate((vBnd, vBnd_lvl), axis=0) + vBnd = np.empty((0, 3)) + tBnd = np.empty((0, 3)) - # output + for i in lstBnd: + lvlBnd = levelsetsTria(triaMidRm.v, triaMidRm.t, hsfBnd, i) + if len(lvlBnd[0][0]) > 0: + vBnd_lvl = np.array(lvlBnd[0][0]) + tBnd_lvl = np.array(lvlBnd[1][0])[:, [0, 1, 1]] - 1 + tBnd = np.concatenate((tBnd, tBnd_lvl + len(vBnd)), axis=0) + vBnd = np.concatenate((vBnd, vBnd_lvl), axis=0) + + # output + + if len(vBnd) > 0: + TriaMesh(v=vBnd, t=tBnd.astype("int")).write_vtk( + filename=os.path.join(OUT_DIR, HEMI + ".mid-surface.hsf-bnd.vtk") + ) + else: + logging.info("Could not create boundary overlays.") - if len(vBnd) > 0: - TriaMesh(v=vBnd, t=tBnd.astype("int")).write_vtk( - filename=os.path.join(OUT_DIR, HEMI + ".mid-surface.hsf-bnd.vtk") - ) else: - logging.info("Could not create boundary overlays.") + logging.info( + "Not creating subfield boundary overlays." + ) diff --git a/images/fig1.png b/images/fig1.png deleted file mode 100644 index ba28263..0000000 Binary files a/images/fig1.png and /dev/null differ diff --git a/images/fig2.png b/images/fig2.png deleted file mode 100644 index 5c9f7cc..0000000 Binary files a/images/fig2.png and /dev/null differ diff --git a/images/fig3.png b/images/fig3.png deleted file mode 100644 index 27d7d1c..0000000 Binary files a/images/fig3.png and /dev/null differ diff --git a/images/fig4.png b/images/fig4.png deleted file mode 100644 index 94f441a..0000000 Binary files a/images/fig4.png and /dev/null differ diff --git a/images/image.png b/images/image.png deleted file mode 100644 index 864fb48..0000000 Binary files a/images/image.png and /dev/null differ diff --git a/images/int_mid_ext.png b/images/int_mid_ext.png deleted file mode 100644 index 3e2f7ca..0000000 Binary files a/images/int_mid_ext.png and /dev/null differ diff --git a/images/labels.png b/images/labels.png deleted file mode 100644 index aac9ac1..0000000 Binary files a/images/labels.png and /dev/null differ diff --git a/images/mask.png b/images/mask.png deleted file mode 100644 index f9a60fd..0000000 Binary files a/images/mask.png and /dev/null differ diff --git a/images/mask_w_surf.png b/images/mask_w_surf.png deleted file mode 100644 index 3984e37..0000000 Binary files a/images/mask_w_surf.png and /dev/null differ diff --git a/images/mid.png b/images/mid.png deleted file mode 100644 index 0ba8c23..0000000 Binary files a/images/mid.png and /dev/null differ diff --git a/images/mid_gauss_curv.png b/images/mid_gauss_curv.png deleted file mode 100644 index e549fc2..0000000 Binary files a/images/mid_gauss_curv.png and /dev/null differ diff --git a/images/mid_hsf.png b/images/mid_hsf.png deleted file mode 100644 index 91cd96d..0000000 Binary files a/images/mid_hsf.png and /dev/null differ diff --git a/images/mid_mean_curv.png b/images/mid_mean_curv.png deleted file mode 100644 index 8b743cd..0000000 Binary files a/images/mid_mean_curv.png and /dev/null differ diff --git a/images/mid_thickness.png b/images/mid_thickness.png deleted file mode 100644 index 3c99b6f..0000000 Binary files a/images/mid_thickness.png and /dev/null differ diff --git a/images/mid_tria.png b/images/mid_tria.png deleted file mode 100644 index d3b17f9..0000000 Binary files a/images/mid_tria.png and /dev/null differ diff --git a/images/ov1.png b/images/ov1.png deleted file mode 100644 index 93a5573..0000000 Binary files a/images/ov1.png and /dev/null differ diff --git a/images/ov2.png b/images/ov2.png deleted file mode 100644 index 4e9f6ac..0000000 Binary files a/images/ov2.png and /dev/null differ diff --git a/images/scaled/ctx-01.png b/images/scaled/ctx-01.png new file mode 100644 index 0000000..d06cb61 Binary files /dev/null and b/images/scaled/ctx-01.png differ diff --git a/images/scaled/ctx-02.png b/images/scaled/ctx-02.png new file mode 100644 index 0000000..7835dd4 Binary files /dev/null and b/images/scaled/ctx-02.png differ diff --git a/images/scaled/ctx-03a.png b/images/scaled/ctx-03a.png new file mode 100644 index 0000000..b9cd36c Binary files /dev/null and b/images/scaled/ctx-03a.png differ diff --git a/images/scaled/ctx-03b.png b/images/scaled/ctx-03b.png new file mode 100644 index 0000000..8133ef5 Binary files /dev/null and b/images/scaled/ctx-03b.png differ diff --git a/images/scaled/ctx-03c.png b/images/scaled/ctx-03c.png new file mode 100644 index 0000000..485804f Binary files /dev/null and b/images/scaled/ctx-03c.png differ diff --git a/images/scaled/erc-01a.png b/images/scaled/erc-01a.png new file mode 100644 index 0000000..1f041b9 Binary files /dev/null and b/images/scaled/erc-01a.png differ diff --git a/images/scaled/erc-01b.png b/images/scaled/erc-01b.png new file mode 100644 index 0000000..9ea60a7 Binary files /dev/null and b/images/scaled/erc-01b.png differ diff --git a/images/scaled/erc-01c.png b/images/scaled/erc-01c.png new file mode 100644 index 0000000..f34e511 Binary files /dev/null and b/images/scaled/erc-01c.png differ diff --git a/images/scaled/erc-02a.png b/images/scaled/erc-02a.png new file mode 100644 index 0000000..46dd718 Binary files /dev/null and b/images/scaled/erc-02a.png differ diff --git a/images/scaled/erc-02b.png b/images/scaled/erc-02b.png new file mode 100644 index 0000000..3d9741a Binary files /dev/null and b/images/scaled/erc-02b.png differ diff --git a/images/scaled/erc-02c.png b/images/scaled/erc-02c.png new file mode 100644 index 0000000..cc43805 Binary files /dev/null and b/images/scaled/erc-02c.png differ diff --git a/images/surf.png b/images/surf.png deleted file mode 100644 index bc0c8f1..0000000 Binary files a/images/surf.png and /dev/null differ diff --git a/images/surf_cut.png b/images/surf_cut.png deleted file mode 100644 index 88dfd1b..0000000 Binary files a/images/surf_cut.png and /dev/null differ diff --git a/images/surf_uncut.png b/images/surf_uncut.png deleted file mode 100644 index 2c6489e..0000000 Binary files a/images/surf_uncut.png and /dev/null differ diff --git a/pyproject.toml b/pyproject.toml index 636d565..b54db47 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -6,7 +6,7 @@ build-backend = 'setuptools.build_meta' name = 'hipsta' description = 'A python package for hippocampal shape and thickness analysis' license = {file = 'LICENSE'} -requires-python = '>=3.9' +requires-python = '>=3.10' authors = [ {name = 'Kersten Diers', email = 'kersten.diers@dzne.de'}, {name = 'Martin Reuter', email = 'martin.reuter@dzne.de'} @@ -25,10 +25,11 @@ classifiers = [ 'Operating System :: Unix', 'Operating System :: MacOS', 'Programming Language :: Python :: 3 :: Only', - 'Programming Language :: Python :: 3.9', 'Programming Language :: Python :: 3.10', 'Programming Language :: Python :: 3.11', 'Programming Language :: Python :: 3.12', + 'Programming Language :: Python :: 3.13', + 'Programming Language :: Python :: 3.14', 'Natural Language :: English', 'License :: OSI Approved :: MIT License', 'Intended Audience :: Science/Research', diff --git a/requirements.txt b/requirements.txt index 9be0f1d..03a6cb1 100644 --- a/requirements.txt +++ b/requirements.txt @@ -1,5 +1,5 @@ nibabel>=5.2 -numpy>=1.22.3, <2.0 +numpy>=1.22.3 pandas>=1.3.0 plotly>=5.6.0 pyacvd>=0.2.9