Surface extraction from Colin27 and remesh

[Daniele Ancora]

Dear Dr. Fang,

first of all I thank you very much for the iso2mesh toolbox, I find it very useful and carefully designed. Congrats!

I'm writing here concerning a problem I recently have using the Colin27 model (version 2L). The basic idea of what I would like to do is the following:

1. isolate node, face and elem for each of the 4 regions (Skull, CSF, gray and white matter)
2. merge the node and face of all of the previously isolated regions
3. remesh the model to get another Colin27 mesh

here what I've done:

load MMC_Collins_Atlas_Mesh_Version_2L.mat

%% Surface exterior boundary extraction 

% 1-scalp, 2-CSF, 3-gray matter, 4-white matter

tempelem = elem;

tempelem(:,5) = double(tempelem(:,5) >= 1);         % I set internal layers to the same value of the external

SKALPelem = tempelem(find(tempelem(:,5) == 1),:);   % I find all the elements having the ID I gave before

SKALPface = volface(SKALPelem(:,1:4));              % find the boundary of the element list, i.e. the surface

tempelem = elem;

tempelem(:,5) = double(tempelem(:,5) >= 2);

CSFelem = tempelem(find(tempelem(:,5) == 1),:);

CSFface = volface(CSFelem(:,1:4));

tempelem = elem;

tempelem(:,5) = double(tempelem(:,5) >= 3);

GRAYelem = tempelem(find(tempelem(:,5) == 1),:);

GRAYface = volface(GRAYelem(:,1:4));

tempelem = elem;

tempelem(:,5) = double(tempelem(:,5) >= 4);

WHITEelem = tempelem(find(tempelem(:,5) == 1),:);

WHITEface = volface(WHITEelem(:,1:4));

% here I exclude all the nodes that are not in the surface, i.e. are

% isolated respect the facelist

[SKALPnode,SKALPface] = meshcheckrepair(node,SKALPface,'isolated');

[CSFnode,CSFface] = meshcheckrepair(node,CSFface,'isolated');

[GRAYnode,GRAYface] = meshcheckrepair(node,GRAYface,'isolated');

[WHITEnode,WHITEface] = meshcheckrepair(node,WHITEface,'isolated');

% adding the ID to the head layers

SKALPface(:,4) = 1;

CSFface(:,4) = 2;

GRAYface(:,4) = 3;

WHITEface(:,4) = 4;

figure;

subplot(221), plotmesh(SKALPnode,SKALPface); title('SKALP external boundary');

subplot(222), plotmesh(CSFnode,CSFface); title('CSF external boundary');

subplot(223), plotmesh(GRAYnode,GRAYface); title('GRAY external boundary');

subplot(224), plotmesh(WHITEnode,WHITEface); title('WHITE external boundary');

axis equal;

Until here everything seems running fine, as you can see from the plottings. But then when i try to merge them back:

[newnode,newface] = mergemesh(SKALPnode,SKALPface(:,1:3), CSFnode,CSFface(:,1:3),...

                              GRAYnode,GRAYface(:,1:3), WHITEnode,WHITEface(:,1:3));

[newnode,newface] = meshcheckrepair(newnode,newface,'isolated');

% define internal region for the elements

p1 = [90 50 132];

p2 = [90 50 127];

p3 = [94 94 130];

p4 = [111 99 135];

regions=[p1;p2;p3;p4];

[node2,elem2,face2] = surf2mesh(newnode,newface(:,1:3),[],[],1,100,regions,[]);

i get various warnings such as:

Warning:  Point 39920 is identical with point 60760. 

Warning:  Point 38972 is identical with point 59812. 

...

Warning:  Two subfaces are found duplicated at (60742, 39906, 60918) 

  Subface of facet #122433 is deleted. 

Warning:  Two subfaces are found duplicated at (40079, 60742, 39906) 

  Subface of facet #122434 is deleted. 

...

 

and then tetgen.exe stop working. 

As you replied in the thread:

https://groups.google.com/forum/?hl=en#!msg/iso2mesh-users/MdJNHCBP3h0/bJRVrzjRf8IJ;context-place=searchin/iso2mesh-users/tetgen$20has$20stopped$20working%7Csort:relevance

it may be due to interacting surfaces, although I seems really not from the plots.

Do you have any hints to accomplish this task? I believed it was quite simple to do but I am not stepping forward even trying other ways.

I thank you for your attention, hoping that you can give me some ideas

D.

[Qianqian Fang]

On 08/30/2016 10:55 AM, Daniele Ancora wrote:

Dear Dr. Fang,

first of all I thank you very much for the iso2mesh toolbox, I find it very useful and carefully designed. Congrats!

I'm writing here concerning a problem I recently have using the Colin27 model (version 2L). The basic idea of what I would like to do is the following:

1. isolate node, face and elem for each of the 4 regions (Skull, CSF, gray and white matter)
2. merge the node and face of all of the previously isolated regions
3. remesh the model to get another Colin27 mesh

hi Daniele

you don't actually need to re-extract the surfaces of each region,
simply use the "face" variable in the mesh package.

here is a short script that works on Linux with tetgen 1.4.3:

load MMC_Collins_Atlas_Mesh_Version_2L.mat
plotmesh(node,face,'x>100')

p1 = [90 50 132];
p2 = [90 50 127];
p3 = [94 94 130];
p4 = [111 99 135];
regions=[p1;p2;p3;p4];

[node2,elem2,face2] = s2m(node,face,1,100,regions);

if you want to remove isolated nodes before calling surf2mesh,
you can do so by inserting the below line before surf2mesh.

[node,face]=removeisolatednode(node,face);

actually, if you don't do so, all isolated nodes will be part of the
output mesh (i.e., the new mesh will contain all the nodes in the old volumetric
mesh).

here what I've done:

....

i get various warnings such as:

Warning:  Point 39920 is identical with point 60760. 

Warning:  Point 38972 is identical with point 59812. 

...

Warning:  Two subfaces are found duplicated at (60742, 39906, 60918) 

  Subface of facet #122433 is deleted. 

Warning:  Two subfaces are found duplicated at (40079, 60742, 39906) 

  Subface of facet #122434 is deleted. 

...

 

and then tetgen.exe stop working.

I've seen problems running the script using tetgen 1.4.3 (the default tetgen used
in iso2mesh). I also noticed that tetgen 1.5 is, sometimes, more robust. So if
you encounter complains from tetgen, try to rename the tetgen1.5.* to tetgen.*
under iso2mesh/bin, sometime the meshing can proceed smoothly.

Qianqian

As you replied in the thread:

https://groups.google.com/forum/?hl=en#!msg/iso2mesh-users/MdJNHCBP3h0/bJRVrzjRf8IJ;context-place=searchin/iso2mesh-users/tetgen$20has$20stopped$20working%7Csort:relevance

it may be due to interacting surfaces, although I seems really not from the plots.

Do you have any hints to accomplish this task? I believed it was quite simple to do but I am not stepping forward even trying other ways.

I thank you for your attention, hoping that you can give me some ideas

D.

--
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 post to this group, send email to iso2mes...@googlegroups.com.
Visit this group at https://groups.google.com/group/iso2mesh-users.
For more options, visit https://groups.google.com/d/optout.

The information in this e-mail is intended only for the person to whom it is
addressed. If you believe this e-mail was sent to you in error and the e-mail
contains patient information, please contact the Partners Compliance HelpLine at
http://www.partners.org/complianceline . If the e-mail was sent to you in error
but does not contain patient information, please contact the sender and properly
dispose of the e-mail.

[Daniele Ancora]

Dear Dr. Fang,

Thank you for the very fast feedback. Maybe I haven't explained well. I really want to extract those surface separately in order to deform (each of) them and to do this I am trying to extract them singularly and then remesh them.

I'm doing so because I want just to slightly change the nodes at the boundary of the various layers!

I hope you can suggest me something
Best
D.

[Qianqian Fang]

On 08/30/2016 02:15 PM, Daniele Ancora wrote:
> Dear Dr. Fang,
>
> Thank you for the very fast feedback. Maybe I haven't explained well. I really want to extract those surface separately in order to deform (each of) them and to do this I am trying to extract them singularly and then remesh them.
>
> I'm doing so because I want just to slightly change the nodes at the boundary of the various layers!

the "face" array contains a label column (4th column), each layer
is separated by a unique label.

does that allow you to do what you have planned?

Qianqian

[Daniele Ancora]

Not actually,

I would like to extract the list of the nodes of each boundary together with the face list to have a set of 4 nodes and faces array, then merge them together. This is what I was trying to do in the first lines of the code I have posted there.

D.

[Qianqian Fang]

On 08/30/2016 03:19 PM, Daniele Ancora wrote:

Not actually,

I would like to extract the list of the nodes of each boundary together with the face list to have a set of 4 nodes and faces array, then merge them together. This is what I was trying to do in the first lines of the code I have posted there.

but I thought you could extract the nodes/triangles for each
layer from node/face as well. maybe I did not fully understand
your plan.

is the below script close to what you want to achieve?

load MMC_Collins_Atlas_Mesh_Version_2L.mat
labels=unique(face(:,end))

for i=1:length(labels)
    [no{i},fc{i}]=removeisolatednode(node,face(face(:,end)==labels(i),1:3));
end
% you know have nodes/faces for each layer, you can make changes, then merge
[noall,fcall]=mergemesh(no{1},fc{1},no{2},fc{2},no{3},fc{3},no{4},fc{4});

plotmesh(noall,fcall,'x>100')

p1 = [90 50 132];
p2 = [90 50 127];
p3 = [94 94 130];
p4 = [111 99 135];
regions=[p1;p2;p3;p4];

[node2,elem2,face2] = s2m(noall,fcall,1,100,regions);

this still works without any issue on my Linux box.

D.

On Tuesday, August 30, 2016 at 9:25:07 PM UTC+3, q.fang wrote:

On 08/30/2016 02:15 PM, Daniele Ancora wrote:
> Dear Dr. Fang,
>
> Thank you for the very fast feedback. Maybe I haven't explained well. I really want to extract those surface separately in order to deform (each of) them and to do this I am trying to extract them singularly and then remesh them.
>
> I'm doing so because I want just to slightly change the nodes at the boundary of the various layers!

the "face" array contains a label column (4th column), each layer
is separated by a unique label.

does that allow you to do what you have planned?

Qianqian

>
> I hope you can suggest me something
> Best
> D.
>

[Daniele Ancora]

Oh great! This is beautifully working smoothly!

I figured out that in my method something was wrong in the extraction of the GREYnode and GREYface (# label=3). In fact not merging that mesh didn't make tetgen crash, but still was missing one piece in the puzzle. I Thank you for you very elegant and educative solution!

There is still a problem though, which is the element labeling after s2m. Apparently I cannot manage to obtain this information after the surface2mesh operation, resulting in a mesh which had lost the information of the element ID.

I'm not such an expert with this toolbox so I would be glad if you can help me on this side as well, I am willing to use this model to run simulations in MMC, after long time using MCX (great software both of them)

Best

D.

[Qianqian Fang]

On 08/31/2016 06:45 AM, Daniele Ancora wrote:

Oh great! This is beautifully working smoothly!

I figured out that in my method something was wrong in the extraction of the GREYnode and GREYface (# label=3). In fact not merging that mesh didn't make tetgen crash, but still was missing one piece in the puzzle. I Thank you for you very elegant and educative solution!

There is still a problem though, which is the element labeling after s2m. Apparently I cannot manage to obtain this information after the surface2mesh operation, resulting in a mesh which had lost the information of the element ID.

are you saying you do not have any the labels in elem2 after calling s2m?
how many columns in elem2?

It works on my system though. did you see anything in the printed logs?

a workaround is to rename tetgen1.5.x to tetgen.x. tetgen1.5 labels all
regions automatically. but I hope you can figure out what was wrong
with your current script.  To debug, you may do

plotmesh(node2,face2,'facealpha',0.2,'linestyle','none');
hold on;
plotmesh(regions,'ro');

and see if your region seeds are located inside each segment. if you
still have trouble, please post your full script.

[Daniele Ancora]

Dear Prof,

thank you again, for the answer. although i renamed the tetgen1.5 to tetgen, the results are still wrong. In particular what is happing right now is that there are more regions than they should be:

load MMC_Collins_Atlas_Mesh_Version_2L.mat

labels=unique(face(:,end))

for i=1:length(labels)

    [no{i},fc{i}]=removeisolatednode(node,face(face(:,end)==labels(i),1:3));

end

% you know have nodes/faces for each layer, you can make changes, then merge

[noall,fcall]=mergemesh(no{1},fc{1},no{2},fc{2},no{3},fc{3},no{4},fc{4});

plotmesh(noall,fcall,'x>100')

p1 = [90 50 132];

p2 = [90 50 127];

p3 = [94 94 130];

p4 = [111 99 135];

regions=[p1;p2;p3;p4];

[node2,elem2,face2] = s2m(noall,fcall,1,100,regions);

figure

subplot(121), plotmesh(node2,face2,'facealpha',0.2,'linestyle','none');

subplot(122), plotmesh(node2,elem2,'facealpha',0.2,'linestyle','none');

hold on;

plotmesh(regions,'ro');

the figure is attached.

the weird thing is that if i run now 

labels=unique(elem2(:,end))

labels =

     1

     2

     3

     4

     5

     6

meaning that the regions in the elements are 6 and not 4 like in the faces. this happens also if i do not put 'regions'

[node2,elem2,face2] = s2m(noall,fcall,1,100);

while with the original tetgen I get only one region in the element list, being it with 5 columns where the last one is always 0.

The problem is only with the element list, not with the face, and cannot figure out what is going on!

Best

[Daniele Ancora]

I have not found any solution to this strange behavior. I have tried also from other pc and operating systems but always find the same problems.

In case any other have tried to run my script and have any suggestion I please you to let me know.

Best regards.

D.