Area of a label

[SherazKhan]

I am wondering, anyone aware of a function, that can estimate label area ?
In Matlab I used to do this as follows:

`function area = triangle_area(faces,verts)

r12 = verts(faces(:,1),:);
r13 = verts(faces(:,3),:) - r12;
r12 = verts(faces(:,2),:) - r12;
area = sqrt(sum(cross(r12',r13').^2,1))/2;
return

function c = cross(a,b)
c = [a(2,:).*b(3,:)-a(3,:).*b(2,:)
a(3,:).*b(1,:)-a(1,:).*b(3,:)
a(1,:).*b(2,:)-a(2,:).*b(1,:)];
return`

`function [PFV,PF]=convert_patch(P,FV,F)
% function [PFV]=convert_patch(P,FV)
% -------------------------------------------------------------
% Given a list of connected nodes (P), gives the structure Faces Vertices of the
% patch
% --------------------------------------------------------------
%
% INPUTS :
% - P : indices of connected nodes
% - FV : structure with fields faces, vertices ,i.e. corresponding to a
% tesselation
% - F : scalar field defined on nodes of FV
% OUTPUTS :
% - PFV : subtesselation of FV based on the nodes P with fields faces,
% vertices
% - PF : scalar field defined on the resulting submesh
% -------------------------------------------------------------------

PFV.vertices=FV.vertices(P,:);
PF=zeros(length(P),size(F,2));
T=[];
nT=1;
for ii=1:size(P,1)
[ind_i,ind_j]=find(FV.faces==P(ii));
for k=1:length(ind_i)
if ismember(ind_i(k),T)
else
convtri=convert(FV.faces(ind_i(k),:),P);
if length(convtri)<=2
else
PFV.faces(nT,:)=convtri;
T=[ind_i(k) T];
nT=nT+1;
end
end
end
PF(ii,:)=F(P(ii),:);
end

function [convtri]=convert(tri,P)

ind1=find(P==tri(1));
ind2=find(P==tri(2));
ind3=find(P==tri(3));

convtri=[ind1 ind2 ind3];`

`% Need freesurfer matlab toolbox
addpath /freesurfer/matlab/

% here provide the surface file you want
[vertices_l, faces_l] = freesurfer_read_surf('lh.pial');
[vertices_r, faces_r] = freesurfer_read_surf('rh.pial');

faces_l=faces_l+1;
faces_r=faces_r+1;

% Computing area in mm2 of the whole brain
area_l =triangle_area(faces_l,vertices_l);

area_r =triangle_area(faces_r,vertices_r);

label_file='xx-lh.label';
labverts = read_label('',label_file);% freesurfer function
labverts = 1+squeeze(labverts(:,1));

FV.faces=faces_l;
FV.vertices=vertices_l;

[PFV,PF]=convert_patch(labverts,FV,zeros(size(vertices_l,1)),1);

area_l_label =triangle_area(PFV.faces,PFV.vertices);`

[larsoner]

I don't know of any existing function, but what you propose (sum up areas of each angle) I'd expect to be fairly quick to compute.

Such a function could go in nibabel or MNE

[SherazKhan]

Just finished coding it in python, now comparing results with mris_anatomical_stats :-)

[agramfort]

and then sending a PR :)

[SherazKhan]

@agramfort @larsoner

Following is the sample script which can go in mn/label.py and mne/tests/test_label.py

Unfortunately, result from my code and freesurfer "mris_anatomical_stats" do not match

I wrote small python wrapper ( compute_area(label_fname, subject)) over freesurfer in the attached code.

There must be some correction factor
https://gist.github.com/SherazKhan/9cbb5520d4efb39bc90f5b175a60a3c5

[SherazKhan]

@agramfort @larsoner any clues, on this discrepancy between estimated areas ?

[agramfort]

no idea how significant is the difference?

[SherazKhan]

Always underestimated by roughly 10-15%, however it depend upon on the surface label is defined, I am currently testing it on white surface.

[larsoner]

There is probably a better algorithm based on the curvature. From some light googling I found this:

https://www.researchgate.net/post/How_can_you_measure_the_surface_area_of_a_mesh_either_directly_or_indirectly_with_computer_graphics

Specifically:

If the mesh represent a plain shell and the triangles are large with respect to the curvature, a better easy approximation is to calculate the average radius of curvature for each triangle and use spherical geometry (area of spherical triangle).

You could also post to the Freesurfer list to get an algorithm description.

[SherazKhan]

I post it on the freesurfer mailing list, but summing area of triangles, suppose to be more accurate, as triangles wrt curvature on high resolution freesurfer surface.

[larsoner]

suppose to be more accurate, as triangles wrt curvature on high resolution freesurfer surface

We know that the triangles are interior to the true smooth, convex surface, so we should be able to leverage our curvature estimates to improve the estimate, even if the algorithm I posted does a poor job.

In any case, we can quantify the error easily enough for a sphere. Let's start with a number of points (163842,) that is probably close to the surface:

def triangle_area(rr, tris):
    r12 = rr[tris[:, 0], :]
    r13 = rr[tris[:, 2], :] - r12
    r12 = rr[tris[:, 1], :] - r12
    cross = np.cross(r12, r13)
    return np.sum(np.sqrt(np.sum(cross * cross, axis=1))) / 2.


surf = mne.surface._get_ico_surface(7)
surf['rr'] /= np.linalg.norm(surf['rr'], axis=-1, keepdims=True)
expected = 4 * np.pi
print(100 * triangle_area(surf['rr'], surf['tris']) / expected)

We get 99.9981, so we're only off by ~0.002%. The white surface isn't anywhere near as smooth as a sphere, but even ico-1 (substituting for the 7 above), which has 42 points, gives us 92.83, so < 8% error. If we drop down to ico-0, which has 12 points, then we get 76.19, so ~24% error. I doubt we're accumulating that much error based on the curvature, the surface would have to be very curved.

Have you tried calculating the surface area of the entire surface to see if that matches? If FreeSurfer includes all triangles that contain any label vertex (as opposed to the triangles that contain only label vertices), that would change the number of triangles used.

[larsoner]

Looking back at the thing I linked above, this seems doable:

calculate the average radius of curvature for each triangle and use spherical geometry (area of spherical triangle)

But it looks like FreeSurfer uses the same def as we do:

https://github.com/freesurfer/freesurfer/blob/e34ae4559d26a971ad42f5739d28e84d38659759/utils/MRIScomputeTriangleProperties_extracted.h#L130

So I'm not sure why there would be a mismatch