Difficulty getting mesh to work

[Daniel Marshall]

Dear Dr. Fang,

I am glad that I found this toolbox, it looks like exactly what we need! I am working on 3D meshing of segmented microCT data which will be imported into comsol for some electrical simulations.

I am running into some issues which I hope you could provide insight on. Some of these may be simple as I am still getting familiarized with the toolbox

1. When specifying opt for vol2mesh, I have 2 sets of parameters, one for the whole nerve volume, and one volume for the fascicles running through the nerve. I would like to set different meshing parameters, for each (fascicles need a finer mesh). However, I need to perform an affine transformation to both meshes, so I specify opt as follows:

opt(1).radbound=4;

opt(2).radbound=2;

opt(1).side='lower';

opt(2).side='lower';

opt(1).A=diag([1 1 5]); 

opt(2).A=diag([1 1 5]); 

opt(1).B=diag([0 0 0]); 

opt(2).B=diag([0 0 0]);

opt(1).dofix = 1

opt(2).dofix = 1

I recieve the following error:

Arrays have incompatible sizes for this operation.

Error in vol2surf (line 195)

            v0=(opt(i).A*v0'+repmat(opt(i).B(:),1,size(v0,1)))';

Error in vol2mesh (line 65)

        [no,el,regions,holes]=vol2surf(img,ix,iy,iz,opt,dofix,'cgalsurf');

Related documentation

However, when I specify a single opt parameter set:

opt=struct('radbound',3,'A',diag([1 1 5]), 'B', [0 0 0]);

and run the same code:

[node,elem,face]=vol2mesh(test,1:size(test,1),1:size(test,2),1:size(test,3),opt,20,1);

It functions without errors. Am I doing something wrong with my A and B transforms?

2. As I mentioned above, the mesh I am trying to create has two regions. I was interested in trying to create surfaces instead of volumes so I could remesh them to remove anisotropies introduced by the Affine transformation. When I tried doing this, vol2surf creates a surface mesh just fine, but running remeshsurf on the 2 region surface results in the inner region disappearing completely (only the outer nerve surface remains).

3. After the above issue I tried to mesh and remesh the surfaces individually, then recombine via mergesurf. When I do so, the output has no region_id column for the elements, though the documentation says there should be one. This results in my two separate surfaces becoming one region. This is also relevant because I am interested in adding some primitives (i.e. meshacylinder) but it seems merging them via mergesurf or surfboolean results in no region_id output, unless I am doing something wrong.

Thanks so much!

-Daniel

4. 

[Qianqian Fang]

hi Daniel, 

using a struct array for opt in vol2mesh won't work for all "method" options. if you use cgalmesh as the method, the correct way to set label-based density is to use the maxvolume input, see this demo script

https://github.com/fangq/iso2mesh/blob/master/sample/demo_label_sizing.m

if you want to extract surface only, you use vol2surf, instead of vol2mesh

do not use mergemesh to combine meshes unless you know these meshes are not spatially intersecting and are completely disconnected. mergemesh only concatenate meshes, and do not do intersection/boolean operations.

for surfaces, you should call mergesurf() instead (which calls surfboolean - it may not work for all inputs, such as the case I mentioned in this earlier email: https://groups.google.com/g/iso2mesh-users/c/iAXbzRxokXc/m/zK5Ryt22AQAJ)

[Daniel Marshall]

Hi Dr. Fang,

Thank you for your reply!

For the struct based opt, the input that is problematic is the linear transformation input. From the vol2surf documentation:

opt(1,2,..).{A,B}: linear transformation for each surface
My interpretation is that the input should be as follows:

opt(1).A=diag([1 1 5]); 

opt(2).A=diag([1 1 5]); 

opt(1).B=diag([0 0 0]); 

opt(2).B=diag([0 0 0]);

This results in an error, is my syntax wrong?

I do want to create mesh eventually, I was using vol2surf in order to create a 2 region surface, which I then wanted to remesh and then mesh into a volumetric mesh.

i.e. vol2surf>remeshsurf>surf2mesh

The problem is my second region (which is interior to the first) disappears during the remesh operation

Is there a way to have the output of mergesurf retain region id's? According to the documentation the output should be

% newnode: the node coordinates after merging, dimension (nn,3)

% newelem: tetrahedral element or surfaces after merging (nn,4) or (nhn,5)

but my newelem output is only (nn,3) dimensions, with no region id.

Thanks!

[Daniel Marshall]

Hi Dr. Fang,

I wanted to ask if you had any more insight on why I get the error with the opt specifications. Am I doing something wrong with specifying both opt sets when using cgal2surf?

Thanks!

-Daniel

[Fang, Qianqian]

On 7/22/21 12:01 PM, Daniel Marshall wrote:

Hi Dr. Fang,

I wanted to ask if you had any more insight on why I get the error with the opt specifications. Am I doing something wrong with specifying both opt sets when using cgal2surf?

Thanks!

-Daniel

On Friday, July 16, 2021 at 3:21:57 PM UTC-4 Daniel Marshall wrote:

Hi Dr. Fang,

Thank you for your reply!

For the struct based opt, the input that is problematic is the linear transformation input. From the vol2surf documentation:

opt(1,2,..).{A,B}: linear transformation for each surface
My interpretation is that the input should be as follows:

opt(1).A=diag([1 1 5]); 

opt(2).A=diag([1 1 5]); 

opt(1).B=diag([0 0 0]); 

opt(2).B=diag([0 0 0]);

This results in an error, is my syntax wrong?

opt.B is supposed to be a 1x3 or 3x1 vector, not a diag matrix

https://github.com/fangq/iso2mesh/blob/master/vol2surf.m#L195

I do want to create mesh eventually, I was using vol2surf in order to create a 2 region surface, which I then wanted to remesh and then mesh into a volumetric mesh.

i.e. vol2surf>remeshsurf>surf2mesh

The problem is my second region (which is interior to the first) disappears during the remesh operation

if the surface of your 2nd region does not intersect with the outer surface, you can directly combine them using mergemesh; if they do intersect, you can merge them use surfboolean or mergesurf. I don't see a need to use remeshsurf.

remeshsurf is a computationally expensive step to fix invalid mesh that usually can't be fixed by meshcheckrepair('meshfix'). It works by rasterizing the surfaces into a voxelated domain, then call cgalsurf to re-extract the surface. unfortunately it can only handle a single region at a time.

if yo have to do this - you can call surf2vol for each of the surface component to produce rasterized volumes over the same grid, and then adding/combing them together to form a multi-label volume, then call v2m('cgalmesh') to do the final mesh generation (or v2s('cgalmesh')->s2m('tetgen1.5')).

Is there a way to have the output of mergesurf retain region id's? According to the documentation the output should be

% newnode: the node coordinates after merging, dimension (nn,3)

% newelem: tetrahedral element or surfaces after merging (nn,4) or (nhn,5)

but my newelem output is only (nn,3) dimensions, with no region id.

again, if the regions (surface/tet) are not overlapping, mergemesh should retain the label, if not, please send me an example

--
You received this message because you are subscribed to the Google Groups "iso2mesh-users" group.
To unsubscribe from this group and stop receiving emails from it, send an email to iso2mesh-user...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/iso2mesh-users/70984ca9-7b90-4fbe-b0c4-2afca0a2e1e4n%40googlegroups.com.

[Daniel Marshall]

Hi Dr. Fang,

Thank you for the advice, I understand what I was doing wrong now! I did have one more question, I managed to get a mesh for my 2 regions with vol2mesh using:

opt=struct('radbound',3,'A',diag([1 1 5]), 'B', [0 0 0]);

but 

opt=struct('radbound',5,'A',diag([1 1 5]), 'B', [0 0 0]);

or 

opt(1).radbound=4;

opt(2).radbound=2;

opt(1).side='lower';

opt(2).side='lower';

opt(1).A=diag([1 1 5]); 

opt(2).A=diag([1 1 5]); 

opt(1).B=[0 0 0]; 

opt(2).B=[0 0 0];

opt(1).dofix = 1

opt(2).dofix = 1

results in the following error:

extracting surfaces from a volume ...

region 1 centroid : 226.500000 162.666667 1.333333

region 2 centroid : 229.000000 187.666667 2.333333

processing threshold level 1...

Surface Mesh Extraction Utility (Based on CGAL 5.0) 

(modified for iso2mesh by Qianqian Fang) 

http://iso2mesh.sf.net 

 

RNG seed 1648335518 

set initial mesh size to 50 

Final number of points: 5549 

processing threshold level 2...

Surface Mesh Extraction Utility (Based on CGAL 5.0) 

(modified for iso2mesh by Qianqian Fang) 

http://iso2mesh.sf.net 

 

RNG seed 1648335518 

set initial mesh size to 50 

Final number of points: 17073 

surface mesh generation is complete

generating tetrahedral mesh from closed surfaces ...

creating volumetric mesh from a surface mesh ...

Error using surf2mesh (line 117)

Tetgen command failed:

Opening C:\Users\dpm42\AppData\Local\Temp\iso2mesh-dpm42\post_vmesh.poly.

Constructing Delaunay tetrahedralization.

Delaunay seconds:  0.191

Creating surface mesh.

Recovering boundaries.

Error:  Invalid PLC! Two subfaces intersect.

  1st (#10582): (77922, 45597, 1342)

  2nd (#39401): (21036, 77925, 80848)

Program stopped.

Error in vol2mesh (line 75)

[node,elem,face]=surf2mesh(no,el,[],[],1,maxvol,regions,holes);

Is there a common cause of this issue?

Thanks!

-Daniel

[Fang, Qianqian]

On 7/26/21 9:15 AM, Daniel Marshall wrote:

Hi Dr. Fang,

Thank you for the advice, I understand what I was doing wrong now! I did have one more question, I managed to get a mesh for my 2 regions with vol2mesh using:

opt=struct('radbound',3,'A',diag([1 1 5]), 'B', [0 0 0]);

but 

opt=struct('radbound',5,'A',diag([1 1 5]), 'B', [0 0 0]);

hi Daniel,

the v2m call with 'cgalsurf' (default) method only works well for a single closed surface.

for multiple surfaces/levelsets, they have to be well separated (can be one-enclosing another, but there must be enough gap between them).

the error below simply shows that your extracted multiple-shell surfaces are self-intersecting (inner/outer surfaces are in contact)

if you know that the two surfaces are supposed to be touching each other, the best way to deal with this is to use

v2s->surface decoupling or boolean ->s2m

or

volume dilation -> v2m

see the example of the 2nd path here

https://github.com/fangq/Colin27BrainMesh/blob/master/colinsmesh_v2.m#L70-L73

where the thickenbinvol grows the inner volume by specified voxels and ensures that the outer layer is well separated with the inner mask.

I often choose the 1st path, the decoupling can be done via surfboolean('decouple').

Qianqian

To view this discussion on the web visit https://groups.google.com/d/msgid/iso2mesh-users/8790cb79-53bf-43a9-9eeb-c9bf34af9998n%40googlegroups.com.

[Daniel Marshall]

Hi Dr. Fang, 

Thank you for the advice, the surfaces should not be touching and is an error in the input, I will clean the input accordingly.

Thanks!

-Daniel

[Daniel Marshall]

Hi Dr. Fang,

I am back with another question which I hope will be a simple fix. I am getting the below error when using bwislands on a 3d volumetric image, but I do not receive the error when passing a 2d image.

Thanks!

Error using imfill>check_locations (line 283)

Function imfill expected its second input argument, LOCATIONS, to have either 1 or NDIMS(IM) columns.

Error in imfill>parse_inputs (line 259)

    locations = check_locations(locations, size(IM));

Error in imfill (line 122)

[I,locations,conn,do_fillholes] = parse_inputs(args{:});

Error in bwislands (line 27)

imgnew=imfill(img,[I,J]);

[Daniel Marshall]

I am also sometimes getting this error when trying to construct a mesh:

Error using surf2mesh (line 117)

Tetgen command failed:

Assertion failed: neightet.tet != dummytet, file tetgen.cxx, line 21703

Opening C:\Users\dpm42\AppData\Local\Temp\iso2mesh-dpm42\post_vmesh.poly.

Constructing Delaunay tetrahedralization.

Delaunay seconds:  0.125

Creating surface mesh.

Recovering boundaries.

Error in vol2mesh (line 75)

[node,elem,face]=surf2mesh(no,el,[],[],1,maxvol,regions,holes);

 

[Qianqian Fang]

On 7/29/21 4:17 PM, Daniel Marshall wrote:

Hi Dr. Fang,

I am back with another question which I hope will be a simple fix. I am getting the below error when using bwislands on a 3d volumetric image, but I do not receive the error when passing a 2d image.

should be fixed with this patch, give it a try

https://github.com/fangq/iso2mesh/commit/8c35d0d707f5c01b86087d4a657efd3a02f4ff39

To view this discussion on the web visit https://groups.google.com/d/msgid/iso2mesh-users/3b7f168f-4230-4cd1-bb17-d381ef1483d5n%40googlegroups.com.

[Qianqian Fang]

On 7/29/21 4:56 PM, Daniel Marshall wrote:

I am also sometimes getting this error when trying to construct a mesh:

tetgen usually gives this error when the input has self-intersection or invalid surface input (with open-edge etc).

please double check, otherwise, upload the surface input somewhere and I am happy to take a look

To view this discussion on the web visit https://groups.google.com/d/msgid/iso2mesh-users/a8c9fa23-8cb9-4445-9cf2-e2c9e760d4afn%40googlegroups.com.

[Daniel Marshall]

Hi Dr. Fang,

I checked the volumes and they seem ok to the best of my knowledge. Is there an email or other channel I can send the file to you privately?

Thanks!

-Daniel

[Fang, Qianqian]

On 7/30/21 4:38 PM, Daniel Marshall wrote:

Hi Dr. Fang,

I checked the volumes and they seem ok to the best of my knowledge. Is there an email or other channel I can send the file to you privately?

you can upload it to a file sharing site, like google drive or dropbox, and send me a direct link; if you don't have access to those, you can also try file.io: https://file.io/

To view this discussion on the web visit https://groups.google.com/d/msgid/iso2mesh-users/ab5fe1b2-12c5-4b7e-9e8c-dcd5734db80cn%40googlegroups.com.