output labels or the data of "elem(:,5)" from surf2mesh

[tmedani]

Dear Dr @fangq ,
I have already posted this issue in brain2mesh :
fangq/brain2mesh#8 (comment)

I thought that we have solved this issue,
however, we still have it....

This is the scenario :
We use 3 surfaces to generate volume mesh to different regions,
let say in the surface I have 3 regions [1 2 3];

When I use the surf2mesh function, I get the three regions, which is great,
but the labels are not sorted correctly.

Here is my script :
the data files are here :
https://github.com/brainstorm-tools/brainstorm3/tree/master/defaults/anatomy/ICBM152

bstAnatomyPath = '\brainstorm\brainstorm3\defaults\anatomy\ICBM152'
%% 1- Load the surfaces
% Here I use the default subject but it can be any of the surfaces
head  = load(fullfile(bstAnatomyPath,'tess_head.mat'));
inner  = load(fullfile(bstAnatomyPath,'tess_innerskull.mat'));
outer  = load(fullfile(bstAnatomyPath,'tess_outerskull.mat'));
brain  = load(fullfile(bstAnatomyPath,'tess_cortex_pial_low.mat'));

%% 2- Merge the surfaces```

[newnode,newelem]=mergemesh(head.Vertices,head.Faces,...
                                                     outer.Vertices,outer.Faces,...
                                                     inner.Vertices,inner.Faces);

% Find the seed point for each region
center_inner = mean(inner.Vertices);
[~,~,~,~,seedRegion1]=raysurf(center_inner,[0 0 1],inner.Vertices,inner.Faces);
[~,~,~,~,seedRegion2]=raysurf(center_inner,[0 0 1],outer.Vertices,outer.Faces);
[~,~,~,~,seedRegion3]=raysurf(center_inner,[0 0 1],head.Vertices,head.Faces);


%% Generate the volume Mesh
% Generate volume mesh
factor_bst = 1.e-6;
maxvol = 0.1; 
keepratio = 1;
regions = [seedRegion1;seedRegion2;seedRegion3];
% The order is important, the output label will be related to this order,
% which is related to the conductivity value.

[node,elem,face]=surf2mesh(newnode,newelem,...
                                                min(newnode),max(newnode),...
                                                keepratio,maxvol*factor_bst,regions,[]);  

First problem:

elemID = unique(elem(:,5))
 % result : elemID = 0 1 2

figure;
subplot(1,2,2)
col = ['r','y','b'];
elemID = unique(elem(:,5));
for ind = 1 : length(elemID)
plotmesh(node,elem((elem(:,5)==elemID(ind)),:) ,'y>0','facecolor',col(ind));
hold on;
grid on; grid minor;
end
legend({'Inner','Outer','Scalp'})
hold on
coor = mean((inner.Vertices));
%quiver3(coor(1),coor(2),coor(3),0,0,0.1,'LineWidth',2);
subplot(1,2,1) ; % figure
plotmesh(inner.Vertices,inner.Faces ,'y>0','facecolor','r'); hold on;    
plotmesh(outer.Vertices,outer.Faces ,'y>0','facecolor','y'); hold on;
plotmesh(head.Vertices,head.Faces ,'y>0','facecolor','b'); hold on;
legend({'Inner','Outer','Scalp'});
grid on; grid minor;

the plot is here : left are the surfaces, right is the volume

For this specific problem since it's stable we choose to fix it manually and we set this : elem((elem(:,5)==0),5) = 3;

This is working in this specific case, but we want to avoid this intervention.
Do you have any recommendation?

Second problem:
Now let assume that we are using this manual correction,
the second problem occurs when we change the keepratio from 1 to 0.8
The output label are not stable also with this parameter.
Here are new results :

We will really appreciate your help to solve these issues.

Thank you in advance for your help.
Best Regards
Takfarinas

[fangq]

@tmedani, sorry for the delay.

the reason you got 0 labels in some of the regions was because your seeds were not generated properly.

because in this case, 3 of the surfaces are enclosed by the outer layer, therefore, raysurf.m may not return the seed that is specific for a given region.

raysurf.m works by shooting a ray across a closed surface, finding a pair of intersection points, and use the mid-point as the seed. This only works when a closed surface is not enclosed/enclosing another surface - because the mid-point may be inside the surface being enclosed.

to do this properly, you need to use something like below:

%% define seeds along the electrode axis
[t,baryu,baryv,faceidx]=raytrace(orig,v0,allnode,allface);
t=sort(t(faceidx));
t=(t(1:end-1)+t(2:end))*0.5;
seedlen=length(t);
seeds=repmat(orig(:)',seedlen,1)+repmat(v0(:)',seedlen,1).*repmat(t(:)-1,1,3);

where the allnode/allface are the combined surfaces of the 3 layers. The raytrace call returns all intersection points with distance from the ray origin (orig) stored in t. So, the mid-point of every segment between adjacent intersection points defines a seed that belongs to a specific region.

Please give above code a try and let me know if you need me to explain this further.

[tmedani]

Dear @fangq
Thank you for your answer and no worry for the delay ;)

I tried your solution, it works fine, thanks.

However, since I have three surfaces (in this example), I tried just to use the 3 upper points,
in this case, the label 0 is still here.
So I took the whole 6 points (length(t)) and the outputs are perfect.

I changed the last line to :
seeds=repmat(orig(:)',seedlen,1)+repmat(v0(:)',seedlen,1).*repmat(t(:),1,3);

I removed the (-1)

Thank you again for your help.

Kind regards
Takfarinas