Check regions consistency

[oesteban]

Regions should be disjoint by definition. However, Freesurfer's surfaces are not disjoint, we need to face this situation somehow:

• We can use directly Freesurfer's surfaces, then check the overlap in binary high-resolution derived masks and with some hierarchical criterium fix the overlapping areas. Then we have to go back to mesh. Problem: we are losing the original features (Freesurfer's surfaces are very precise, topologically correct)
• Start with high-resolution binary masks since the beginning for specifying the prior and compute the mesh inside -> anyways, we are always losing freesurfer's properties.

This issue is related to the problem of separated hemispheres.

[oesteban]

@zosso said:

Every mesh separates exactly two regions. And also, this is what freesurfer gives you etc. Also, even if the segmentation model is a priori region based (I'll come back to this in a second), through the theory of shape gradients (to be re-included in the methods section), you get flow problems uniquely concerned with the regions' boundaries, i.e. the very contours I was talking about. (this stems from Green's identities, if you wish).

Well, every mesh separates exactly two regions holds if and only if:

1. The mesh is the surface of a closed unique region (what is not true in freesurfer output, given that we have one region -the white matter- splitted through 2 contours -the hemispheres-).
2. The meshes do not intersect. I don't know the CSF meshes we can get from freesurfer, but in the best of the cases -no intersections at all, because somebody cared about fixing overlaps-, we can have paths where 3 meshes meet. This means that a mesh can limit with 2 or more regions, AND a vertex could belong to 2 or more surfaces -e.g. if 3 meshes meet, then it's 3 surfaces-.

Regarding 1, we have 2 options:

• investigate if we can get the wh and pial surfaces with the hemispheres merged, or
• see if the nodes where they merge are marked somehow, so we can easily reproduce the merge assuming small errors in the normal computation in areas close to the intersection and skip shape gradients computation in these marked vertices. Please, think that the gradients in the corpus callosum are always trying to grow until the other hemisphere is covered. Yes, the other contour makes exactly the opposite, but I don't think it grows with exactly the same strength (specially if we are estimating parameters on each roi).

Regarding 2, if a mesh is surrounded by more than 1 region, then I need to keep track of which region is outside every vertex. This is not difficult to compute, but we need to voxelize the regions (I keep discussion around this after your comments).

THe only moment where you actually need the regions rather than the contours, is for the evaluation of the region statistics (mean and covariance). Here, you need to make sure that the point samples cover only representative voxels of your region, so yes you would want to exclude any csf when sampling the interior of the WM-GM surface. I should have written somewhere in the very first draft something like this: you define your point samples for each region in the original T1 domain (where it is rather easy to get accurate anatomical segmentation within the white matter for example). Also, you can easily think of an appropriate subsampling (in the sense of density of point samples) in order to reduce the computation time at each iteration of statistics-update. Indeed, these point samples in T1 space are mapped to Diffusion space through the current estimate of the deformation field u; there, you interpolate the features for each such mapped sample point and use them to estimate your region's statistics.

This is already implemented and it works. The only else thing that I can do is to establish that the way I add priors to the method is hierarchical (e. g. the layers in Photoshop). What I add first/last is on top and the pixels included in the top region will be removed from lower regions (if they also cover the top region).

Long story short: for each contour mesh, you simply need two pointers to lists of point samples representative for the regions on either side of the contour. One such list of point samples can then be shared by different contours (for example, the GM samples are used as inside region of the pial surface, and outside region of the WM-GM surface...)

I think that, once solved the first problem (surfaces limiting with more than 1 other region) the hierarchical solution is the easiest and we can get rid of having long lists of surrounding samples per vertex. I only would need to keep an identifier of the outer region on each node. And, as we do not intend to change topology (nor the surfaces are allowed to slide one respect another), that should be enough.

What do you think?

[oesteban]

Ok, I've checked with pial and white surfaces and it seems to work out: pial is correctly emptied of white.

Now, I should proceed with the second part: associate the corresponding outer region with each node and not with the inner region (as it is now implemented).

Still, it could be an issue having surfaces split in hemispheres, given the extracted parameters:

Region 0:
    Inside region
        Mean = [0.38439, 0.00082918]
        Covariance Matrix 
                0.0294945 -1.78823e-05
                -1.78823e-05 1.5602e-07

    Outside region
        Mean = [0.0372162, 0.000159651]
        (Inverse) Covariance Matrix 
                0.013507 2.56427e-05
                2.56427e-05 2.09818e-07

Region 1:
    Inside region
        Mean = [0.390598, 0.000800097]
        Covariance Matrix 
                0.0286 -1.24119e-05
                -1.24119e-05 9.87993e-08

    Outside region
        Mean = [0.0371323, 0.000160705]
        (Inverse) Covariance Matrix 
                0.013444 2.57165e-05
                2.57165e-05 2.12753e-07

Region 0 is lh.white and region 1 is rh.white. Uneven gradients could happen.

I pospone this issue to the next milestone, and I open a new one regarding the new approach (outer surface attached to each vertex).

[oesteban]

Example of the regions generated by lh.white, rh.white and lh.pial with the new implementation: