Using qmeshcut to compute an isosurface from an arbitrary grid of points

[ftadel]

Hello,

We have an arbitrary grid of 3D points, one value per point, and would like to use qmeshcut to generate isosurfaces at different levels. The output surface of qmeshcut can be plotted correctly (although I'm questioning the orientation of the faces). But it is not considered as a closed mesh, and it can't be repaired with meshcheckrepair.

Here is some example data: example.zip

And an example script:

load example.mat;
dt = delaunayTriangulation(GridLoc);
[P,v,Faces] = qmeshcut(dt.ConnectivityList, GridLoc, V, 9760);

>> VTA = surfvolume(P, Faces);
Error using surfvolume
open surface is detected, you have to close it first, consider meshcheckrepair() with meshfix option

>> [P, Faces] = meshcheckrepair(P, Faces, 'meshfix')
Fixing asin tolerance to 1.745329e-004 
meshfix C:\Users\franc\AppData\Local\Temp\iso2mesh-Francois\pre_sclean.off 
Cleaning intersections, degeneracies ... 
Saving output mesh to 'C:\Users\franc\AppData\Local\Temp\iso2mesh-Francois\pre_sclean_fixed.off' 

P =
     []

Faces =
     []

Am I doing something wrong here?

[fangq]

@ftadel, try this instead

load example.mat
dt = delaunayTriangulation(GridLoc);
[no,val,fc] = qmeshcut(dt.ConnectivityList, GridLoc, V, 9760);

fc=[fc(:,[1,2,3]); fc(:,[1,3,4])];     % convert quads into triangles
[no1,fc1]=meshcheckrepair(no,fc,'dup');    % remove duplicate nodes
[no1,fc1]=meshcheckrepair(no1,fc1,'meshfix');   % generate a close surface
[no2,el2]=s2m(no1,fc1,1,0.001);    % test if the surface is water-tight

the main changes are: the facedata output of qmeshcut stores a quad per triangle - because when cutting a tetrahedron, the cross-section may have up to 4 nodes. I have to convert it to triangles first before calling meshcheckrepair

also, qmeshcut does not merge nodes, so one has to call meshcheckrepair with dup option to remove duplicated nodes/elements.

on the other side, I assume you could also call convhull to get the surface, or you could use isosurface.

[fangq]

I want to add that if you want to remap your val vector after removing duplicated nodes, you can add this before the meshfix call

[temp, ia]=ismember(no1, no, 'rows');
newval=val(ia);

although in this case, it does not matter because your val is a constant.

[ftadel]

Thanks for your prompt and useful reply!

It works great, except that in the specific use case, we may have separate blobs (reference and recording site for SEEG recordings). We don't know a priori the number of separate objects in the isosurface, it depends on the threshold. If I call meshcheckrepair(...,'meshfix'), it keeps only one connected object. Would you have a solution for dealing with multiple closed objects?

I assume you could also call convhull to get the surface, or you could use isosurface.

The Matlab functions convhull and boundaries give one object only, we'd need to split the different blobs first. The function isosurface can only be executed on a regular grid, with the values V being a volume (3D matrix). For this application, we need to be able to work with arbitrary clouds of points.

[fangq]

@ftadel, instead of using qmeshcut, why not using v2s or vol2surf to extract the surfaces? they are guarantee to be self-intersection free, and are relatively fast.

the only downside is that when there are multiple "islands" in your volume, v2s may miss some objects, it is a known issue of CGAL, see #7, if this is an issue, you can use a workaround , see FAQ 9.

alternatively, you can compile a new cgalsurf binary using this https://github.com/fangq/iso2mesh/pull/41

[ftadel]

@ftadel, instead of using qmeshcut, why not using v2s or vol2surf to extract the surfaces?

In input we may have an arbitrary grid of points with one value per point, not necessarily a volume. For example: the nodes from a tetrahedral mesh that is denser closer to the areas of interest, such as the contacts of the SEEG shafts or holes in the skull. In order to use isosurface, v2s or vol2surf, I would need to interpolate the grid values on a volume first, right?

Are there efficient tools in iso2mesh to do this, i.e. interpolating a grid of 10,000 to 1,000,000 points on a volume? Or what would be your preferred approach in Matlab?

[fangq]

@ftadel, sorry, I kept forgetting, you are creating isosurfaces from unstructured mesh data. interpolation can be quite inefficient depends on the domain.

I think the qmeshcut approach is currently the best approach for what iso2mesh provides. to handle separate/disjointed surfaces, you first call finddisconnsurf(), which returns one connected component per cell element, and then you can handle each surface from the output.

[ftadel]

Thank you for this nice new track. With finddisconnsurf, I can split the different blobs, then call meshcheckrepair('meshfix') on each blob.

I tried executing something like this, in order to compute the isosurface + its volume:

        % Compute isosurface
        [P,v,fc] = qmeshcut(dt.ConnectivityList, GridLoc, V, Thresh);
        % Convert quads into triangles
        Faces = [fc(:,[1,2,3]); fc(:,[1,3,4])];
        % Remove duplicate nodes
        [P,Faces] = meshcheckrepair(P, Faces, 'dup');

  % === THE FOLLOWING IS ONLY TO COMPUTE THE VOLUME "VTA"
        % Separate different blobs (possibly one around each active contact)
        splitFaces = finddisconnsurf(Faces);
        % Process each disconnected element separately
        VTA = 0;
        nodes = cell(2,length(splitFaces));
        for i = 1:length(splitFaces)
            % Generate a closed surface
            [nodes{1,i}, nodes{2,i}] = meshcheckrepair(P, splitFaces{i}, 'meshfix');
            % Test if the surface is water-tight
            % [nodes{1,i}, nodes{2,i}] = s2m(nodes{1,i}, nodes{2,i}, 1, 0.001);
            % Compute its volume
            % VTA = VTA + surfvolume(nodes{1,i}, nodes{2,i});
        end
        % Concatenate all the fixed blobs together
        [P, Faces] = mergesurf(nodes{:});

However, I still can't run surfvolume or s2m on each object, because of Tetgen errors. With this test data: test.zip

K>> surfvolume(node, face)
generating tetrahedral mesh from closed surfaces ...
creating volumetric mesh from a surface mesh ...
Error using surf2mesh
Tetgen command failed

Error in fillsurf (line 20)
[no,el]=surf2mesh(node,face,[],[],1,1);

Error in surfvolume (line 26)
[no,el]=fillsurf(node,face);

And:

K>> s2m(node, face, 1, 0.001)
generating tetrahedral mesh from closed surfaces ...
creating volumetric mesh from a surface mesh ...
Error using surf2mesh
Tetgen command failed

Error in s2m (line 40)
    [node,elem,face]=surf2mesh(v,f,[],[],keepratio,maxvol,regions,holes);

Another question I have is related to the version of Tetgen used by these two functions: it seems that we can't specify 'tetgen1.5', like when calling directly surf2mesh. Is there any technical reason for this? or is it just because it wasn't worth the effort to edit all the functions calling surf2mesh?

I'm sorry, this might not be the best channel to discuss the computation of the volume of this mesh... Maybe I should close this issue and open a new one?

@tmedani Do you have enough information to pick it up from here? Maybe you're much better than I am with these surface manipulations. Now you have it in the Brainstorm interface, you can play with it a bit a try to improve it? Including the computation of the VTA with this new approach? (note that we also need to identify the activated areas that touch the limit of the source space, as qmeshcut does not close the surface at all - if we close then later with meshcheckrepair('meshfix'), it does not follow the surface of the brain...)

[fangq]

@ftadel, the test data works fine on my Linux and windows boxes, see screenshot, are you using a mac?

Another question I have is related to the version of Tetgen used by these two functions: it seems that we can't specify 'tetgen1.5', like when calling directly surf2mesh. Is there any technical reason for this? or is it just because it wasn't worth the effort to edit all the functions calling surf2mesh?

s2m does support the method input https://github.com/fangq/iso2mesh/blob/master/s2m.m#L4

but if you want to change all other dependencies (like fillsurf) of s2m/surf2mesh, right now there is no flag to do that, unless you link tetgen binary to tetgen1.5 inside bin/

Update: I tested the above script on a Mac (imac, intel i5), it also worked.

[fangq]

However, I still can't run surfvolume or s2m on each object, because of Tetgen errors.

does your per surface node list node{1,i} contain the full node list or just the node of the single component?

to avoid using isolated nodes, you can do

[nodecomp, facecomp]= meshcheckrepair(P, splitFaces{i}, 'isolated');

before passing it to meshcheckrepair(..., 'meshfix')

[ftadel]

I reinstalled completely Iso2mesh and that solved my Tetgen issue. I'm not sure what I had done but I guess there was something wrong with the binaries... Sorry I didn't try that before.

I think we have everything we needed to keep working for now. Thanks for all your precious help!