solution for thegap in splitbrain.m

[aarongeller]

splitbrain.m often chooses the wrong mesh set, causing flipped brains and requiring manipulation of the variable thegap to force the right choice. I added a function orientation_good that solves this. My version of splitbrain supports slicing for coronal, axial, sagittal and "oblique coronal" views and is pasted below.

function split = splitbrain(cortex, orientation, b, m)
% Splits the brain and exports a mesh that is a subset of the original
% mesh. Requires the function splitFV (Matlab Exchange)
%
%
% INPUTS:
% cortex - the original mesh. a structure containing tri and vert 
% 
% orientation - slice orientation ('c', 'a', 's', or 'oc'), with
% slices defined as follows, for axes defined medial/lateral = X
% axis, anterior/posterior = Y axis, inferior/superior = Z axis:
%
% c (coronal): defined by a line in the XY plane, angle near n*pi
%              for n in integers including zero
% a (axial): defined by a line in the XZ plane
% s (sagittal): defined by a line in the XY plane, angle near
%               n*pi/2 for n in integers excluding zero
% oc (oblique coronal): defined by a line in the YZ plane
%
% b- intercept
%
% m- slope
%
% OUTPUTS:
% split - a new mesh split at the specified point. A structure containing
% tri and vert

%     (This is a subfunction created as a part of) Omni-planar and surface
%     casting of epileptiform activity (OPSCEA) (UC Case Number SF2020-281)
%     jointly created by Dr. Jon Kleen, Ben Speidel, Dr. Robert Knowlton,
%     and Dr. Edward Chang is licensed for non-commercial research use at
%     no cost by the Regents of the University of California under CC
%     BY-NC-SA 4.0 (https://creativecommons.org/licenses/by-nc-sa/4.0/).
%     Please contact innovation@ucsf.edu if you are interested in using
%     OPSCEA for commercial purposes.

%     The following copyright notice and citation is to be included in any
%     publication, material or media wherein all or a part of Licensed
%     Material is contained, “Certain materials incorporated herein are
%     Copyright © 2016 The Regents of the University of California
%     (REGENTS). All Rights Reserved.

%     Please cite the following paper in your publications if you have used
%     our software in your research, as well as any relevant toolboxes used
%     herein as appropriate (img_pipe, FreeSurfer): Kleen JK, Speidel B,
%     Baud MO, Rao VR, Ammanuel SG, Hamilton LS, Chang EF, Knowlton RC.
%     Accuracy of omni-planar and surface casting of epileptiform activity
%     for intracranial seizure localization. In press at Epilepsia.”

%mind the gap (in mm): prevents triangles that span the gap from inducing unwanted mesh
thegap = max([abs(5.*m) 10]);

% NOTE; occasionally produces unwanted side effects for certain
% slices where a physical gap is inserted between the slice and cortex.
% Adjusting thegap can help for individual cases but will need a more 
% unified solution in a future iteration of OPSCEA software.

if strcmp(orientation, 'c')
    if (min(cortex.cortex.vert(:,2)) < b) && (b < max(cortex.cortex.vert(:,2)))
        [idx, ~] = sort(find(abs((cortex.cortex.vert(:,2) - (m.*cortex.cortex.vert(:,1) + b + thegap))) <= thegap)); 
        % get indices of verts with y between (mx + b) and (mx + b + 2*thegap)
    else 
        % This means intercept lies outside A-P extent of the mesh
        display('Weird case in splitbrain.m...');
        [idx, ~] = sort(find(abs((cortex.cortex.vert(:,2) - (m.*cortex.cortex.vert(:,1) + b + thegap))+2*thegap)<=thegap));
    end

elseif strcmp(orientation, 'a')
    if (min(cortex.cortex.vert(:,3)) < b) && (b < max(cortex.cortex.vert(:,3)))
        [idx, ~] = sort(find(abs((cortex.cortex.vert(:,3) - (m.*cortex.cortex.vert(:,1) + b - thegap))) <= thegap)); 
        % get indices of verts with z between (mx + b - 2*thegap) and (mx + b)
    else 
        % This means intercept lies outside S-I extent of the mesh
        display('Weird case in splitbrain.m...');
        [idx, ~] = sort(find(abs((cortex.cortex.vert(:,3)-(m.*cortex.cortex.vert(:,1) + b + thegap))+2*thegap)<=thegap));
    end

elseif strcmp(orientation, 's')
    [idx, ~] = sort(find(abs((cortex.cortex.vert(:,2) - (m.*cortex.cortex.vert(:,1) + b - thegap))) <= thegap)); 
    % get indices of verts with y between (mx + b - 2*thegap) and (mx + b)

elseif strcmp(orientation, 'oc')
    [idx, ~] = sort(find(abs((cortex.cortex.vert(:,3) - (m.*cortex.cortex.vert(:,2) + b + thegap))) <= thegap)); 
    % get indices of verts with z between (my + b + 2*thegap) and (my + b)
end

if isempty(idx)
    split.vert = cortex.cortex.vert;
    split.tri = cortex.cortex.tri;
else
    mesh.tri = delete_verts(cortex.cortex.tri, idx);        
    mesh.vert = cortex.cortex.vert;

    disp('Generating partial mesh for slice view...')
    FVout = splitFV(mesh.tri,mesh.vert);
    fv.vert = FVout(1).vertices;
    fv.tri = FVout(1).faces;

    if ~orientation_good(fv.vert, m, b, orientation)
        display('splitFV set #1 has wrong orientation, choosing set #2...');
        fv.vert = FVout(2).vertices;
        fv.tri = FVout(2).faces;
    end

    split = fv;
end

function tri = delete_verts(tri, idx)
while(~isempty(intersect(tri(:,2), idx)))
    [~,ia2,~] = intersect(tri(:,2), idx); %looks for the tris that match the indices of the verts on the line
    tri(ia2,:) = []; %remove the tris from the mesh dataset that have already been counted
end

function status = orientation_good(verts, m, b, orientation)
% check that the chunk we're choosing is on the right (correct) side
% of the line.

status = 0;
centroid = mean(verts);

if strcmp(orientation, 'c') && centroid(2) < m*centroid(1) + b
    % for coronal cut want to choose the part behind the plane (we're
    % looking back) so we want centroid below the line
    status = 1;
elseif strcmp(orientation, 'a') && centroid(3) > m*centroid(1) + b
    % for axial cut we want to choose the part above the plane
    % (we're looking up) so we want centroid above the line
    status = 1;
elseif strcmp(orientation, 's') && centroid(1) > (centroid(2) - b)/m
    % for sagittal cut we want the part to the right of the plane
    % (we're looking to the right)
    status = 1;
elseif strcmp(orientation, 'oc') && centroid(2) < (centroid(3) - b)/m
    % for oblique coronal cut we want the part to the right of the
    % plane when viewed from the side (i.e. the posterior part)
    status = 1;
end

[raphaelchristin]

Hi Aaron, thanks again for your contribution. Your changes are now incorporated in the main branch of the repository. Feel free to open another issue if you find another bug!

[aarongeller]

cool, thanks Raphael!

[raphaelchristin]

Hi @aarongeller, I am reopening this issue because it seems that your fix for the orientation issue you encountered caused a bug on our end: #10 explains the situation. I will be reverting back to our old splitbrain.m to avoid the issue, but I will also open a new branch to investigate what in your fix caused the issue. The hope is to find the issue and bring back your version once we can fix the bug.

To kickstart this process, could you provide a bit more detail on the issue that prompted you to write this fix? Does it only happen when you use views other than the 'c' view? We do not use anything else than coronal views and are thinking of removing the orientation parameter altogether in a future version. However, if this is something that you use, I am happy to work toward finding a solution so that we retain the functionality.

[raphaelchristin]

The branch I will use to work on retaining your fix is multiple-orientations-splitbrain-test.

[aarongeller]

Hi Raphael-

The issue with brain flipping was occurring with any orientation, including 'c'. It was caused because the mesh partition provided by splitFV was not in the expected order. It was surmountable by adjusting the variable "thegap" and I believe that is was the original coder was referring to with the comment in the code regarding that variable. Obviously if my version is causing issues for you reverting is reasonable.

[raphaelchristin]

Yes I understand your issue. I am new to the codebase so sometimes it takes me some time to understand the code completely, I hope you'll forgive me. I have reverted to our old version and it seems to have fixed the issue for now, although I would like to investigate and find a way to integrate your changes.

[raphaelchristin]

I'll be looking into it over the next week/longer if necessary. Like I mentioned previously, the aim is to bring back your version in the future.

[raphaelchristin]

Sorry if this causes you trouble, I hope we can fix this quickly. Thanks again for all your contributions!

[aarongeller]

Can you send me a dataset which exhibits issue 10? If I get a chance I can take a look at how it runs on my end.

[raphaelchristin]

I'll look into that!

[raphaelchristin]

Linking #15 for documentation, this issue becomes a priority.

[raphaelchristin]

Found the bug! We used the wrong electrodefile for patients with high density grids. Will fix the patients and merge this fix again.

[aarongeller]

cool!

[raphaelchristin]

I hope to be able to fix #18 before bringing this back. It does not seem as though this would cause an issue in new places (except for one particular patient on our end), but I still hope to figure everything out before merging.

[raphaelchristin]

Hi Aaron, just to clarify the issue you were facing: you were encountering some brains that were flipped like in the following image?

If I understand correctly, your fix chooses the opposite side of the mesh in those cases to make sure that the slice faces the viewer. The current problem, and the only thing that stops me from bringing back your changes, is that on our end it seems that the side of the brain that is facing the viewer is not made transparent. Therefore, it hides from view the side of the brain we are interested about. Here is an example of this:

I am wondering if you made any other changes to the codebase, for this fix or maybe even before. I am especially curious about changes to the OPSCEAsurfslice.m file.

Once again, thanks a lot for all your contributions. Hopefully we can resolve this soon.

[aarongeller]

Hi Raphael-

Actually, I think my use of "flip" was different. It appears that the brain for ID in the top figure is "flipped" in the sense that it's not oriented optimally for the viewer. What I was seeing pertained just to slicing as done by OPSCEAsurfslice for depth electrodes, and involved display of the wrong half of the brain. An example for your plot would be that the anterior half of the brain would be shown for AD instead of the posterior one, and therefore the depth electrode would be blocked from view.

Now, if the issue you're referring to is how the PID electrode is displayed in the bottom figure is splitting the brain with residual left hemisipheric mesh, that's something I handle explicitly by defining a sagittal view ('s' mode) in splitbrain.

I have made several modifications both to OPSCEA.m and OPSCEAsurfslice.m and you're welcome to have my versions, if you'd like. My version of OPSCEA does not use the sagittal view, but I forked my own version of OPSCEA just for structural visualization that does. You're welcome to any of this code.

Best,

Aaron

[raphaelchristin]

Hi Aaron, update on this issue: I was able to find a fix! I moved your orientation_good() function to OPSCEAsurfslice() and used it to rotate the brain in the cases where needed. The current version fixes the issues I have encountered on my end, and I believe that it would fix yours as well, although I cannot test. If something is not addressed, feel free to open an issue and we will work to fix any problem that comes up. Once again, thanks for your help and patience!

[raphaelchristin]

Hi Raphael-

Actually, I think my use of "flip" was different. It appears that the brain for ID in the top figure is "flipped" in the sense that it's not oriented optimally for the viewer. What I was seeing pertained just to slicing as done by OPSCEAsurfslice for depth electrodes, and involved display of the wrong half of the brain. An example for your plot would be that the anterior half of the brain would be shown for AD instead of the posterior one, and therefore the depth electrode would be blocked from view.

Now, if the issue you're referring to is how the PID electrode is displayed in the bottom figure is splitting the brain with residual left hemisipheric mesh, that's something I handle explicitly by defining a sagittal view ('s' mode) in splitbrain.

I have made several modifications both to OPSCEA.m and OPSCEAsurfslice.m and you're welcome to have my versions, if you'd like. My version of OPSCEA does not use the sagittal view, but I forked my own version of OPSCEA just for structural visualization that does. You're welcome to any of this code.

Best,

Aaron

I don't think that that will be necessary for now, but thanks very much for the offer! I will keep it in mind. I have been thinking that one day it would be good to merge everything we both have into one cohesive codebase that we both us. If that would be of interest to you we can discuss that in the future. Thanks again for contributing!

[aarongeller]

sounds good, happy to contribute.