Depth Profile / Array

[wsnyder4]

Dear BrainVISA experts,

I am hoping to play around with the depth profile for each sulcus, getting the full list of values for the hull junction points geodesic distances to the bottom points. My end goal is to get these points for all folds with a given sulcus label, and compute the standard deviation or MAD for this depth profile array.

Is there any way to get an array for all depth values along a sulcus in python? I am wondering if the function in learn_orientation_prior.py is of use here:

def compute_depthmap(subject, graph, skel, write_data=False):
    fat = aims.FoldGraphAttributes(skel, graph)
    fat.prepareDepthMap()
    depth = fat.getDepth()
    depth_factor = fat.getDepthfactor()
    sX, sY, sZ, sT = depth.getSize()[:4]
    converter = aims.Converter_Volume_S16_Volume_FLOAT()
    depth2 = aims.Volume_FLOAT(sX, sY, sZ, sT)
    converter.convert(depth, depth2)
    depth = depth2
    depth[:] /= depth_factor
    if write_data:
        aims.write(depth, '%s_depth.ima' % subject)
    return depth

Thank you again for your assistance! Any point in the right direction is much appreciated.

All the best,
Will

[denisri]

Hi Will,
Yes this is the way we compute depths maps for the sulci. Using this function will provide the geodesic depths from the brain hull to the bottom of sulci. Then you can get the depth values on the points you need (like bottom voxels for a given sulcus, or each sulcus) to do the stats.
Denis

[wsnyder4]

Thank you Denis for the quick help! One additional quick question: does BV / Morphologist preprocess this array any further before calculating mean depth? I have been able to replicate max depth but not mean depth measures (output by the morphometry process), finding that my mean depth from this array is always lower than expected. I notice the array has many values of 1, but even excluding these still leads to an underestimate of mean depth relative to what the morphometry process outputs.

Thank you again,
Will

[wsnyder4]

For reference, my code is below, using the test case of "S.C.sylvian._left" to compare my mean depth to the mean depth from the morphometry process.

from soma import aims
from soma import aims, aimsalgo
import numpy as np
test_graph = aims.read("left_sulcal_graph.arg")
test_skel = aims.read("left_sulcal_skeleton.nii")

def compute_depthmap( graph, skel):
    fat = aims.FoldGraphAttributes(skel, graph)
    fat.prepareDepthMap()
    depth = fat.getDepth()
    depth_factor = fat.getDepthfactor()
    converter = aims.Converter_Volume_S16_Volume_FLOAT()
    depth2 = aims.Volume_FLOAT(depth.dimX(), depth.dimY(), depth.dimZ(), depth.dimT())
    converter.convert(depth, depth2)
    depth = depth2
    depth[:] /= depth_factor
    return depth

depth = compute_depthmap(test_graph, test_skel)
depth_array = depth.arraydata()


vals = []
for v in test_graph.vertices():
    if 'label' in v:
        if v['label'] == "S.C.sylvian._left":
            bottom_map = v['aims_bottom']
            for i, p in enumerate(bottom_map[0].keys()):
                p2 = (0, p[2], p[1], p[0])
                sulcus_depth = depth_array[p2]
                vals.append(sulcus_depth)

np.max(vals)
np.mean(vals)

[denisri]

At the first glance I don't see what's different from what we do. We actually record depth values for every bottom voxel, as you are doing. Thus it would need to test things throughly to understand. In the C++ source code of FoldGraphAttributes, I read:

  float			vvol, nvvol = 0, tmpf1;

// [...]

  for( iv=_graph.begin(); iv!=ev; ++iv )
    {
      v = *iv;
// [...]

      rc_ptr<BucketMap<Void> >	ssb, btb, otb;
      v->getProperty( "aims_bottom", btb );
      v->getProperty( "aims_other", otb );
      if( v->getProperty( "aims_ss", ssb ) && !ssb->empty() )
        {

// [...]

          if( btb && !btb->empty() )
            {
              np += btb->begin()->second.size();
              v->setProperty( "bottom_point_number", 
                              (int) btb->begin()->second.size() );
// [...]
              // mean depth
              tmpf1 = 0;
              for( ib=btb->begin()->second.begin(), 
                     eb=btb->begin()->second.end(); ib!=eb; ++ib )
                tmpf1 += _depth( ib->first );
              if( btb->begin()->second.size() != 0 )
                {
                  tmpf1 /= btb->begin()->second.size() * depthfactor;
                  v->setProperty( "mean_depth", tmpf1 );
                }

// [...]
        }
    }

which looks pretty similar to what you are doing in python. Here depthfactor is applied in the division, whereas for you it is applied previously in compute_depthmap, and in FoldGraphAttributes this is computed for every vertex where you are doing it by label. To produce the morphometry measurements, we combine vertex-wise measurements for all vertices with the same label, by re-weighting them by their bottom size (bottom_point_number), so it should be the same result at the end.
There must be a little thing somewhere that I don't see right now...

[wsnyder4]

Thank you for the extra info here -- I tried applying the depthfactor in the same way as the C++ code but still got my underestimated mean depth value. Additionally, S.C.sylvian._left had only one fold in my graph, so I didn't need to weight by bottom_point_number.

I do notice that I have 17 as my bottom_point_number despite there being 64 points in the aims_bottom buckets -- but perhaps I'm misunderstanding their relation.

I am wondering if I could store the values directly from the C++ code. Is there any way I could manipulate this from my own BrainVISA install to obtain the verified values used to compute mean depth?

[denisri]

That's weird. When I run your code on an example for me, I actually get the same values. Are you sure your graph is up-to-date ? Haven't you recomputed the skeleton without recomputing the graph ? Or can you try using the brainvisa process "morphologist/sulci/graphmanipulations/compute cortical graph upgrade from old" ?

[wsnyder4]

I should also note, I am using BrainVISA version 5.0.3, and the graphs were labeled with the CNN method and rewritten over unlabeled graphs. About half of my graphs have version 5.0 and the other half are 3.1 and I am not sure why any particular graph received 5.0 or 3.1
I have a bit of a custom adaptation of the brainvisa pipeline to process my data (mostly by specifying my own input and output names), so it could very likely be an issue on my side! I will keep looking in and see if using an updated BrainVISA or using the standard pipeline with the standard databasing fixes my issues.

[denisri]

Uh oh, you're right. We have found a big problem...
When a graph is written, and if it already exists at the output location on disk, the default is to save only elements marked as "modified". This is an optimization to save time when, for instance, manually labeling sulci without modifying buckets. But here the process actually modifies things and doesn't mark them as modified. So an older labeled graph will stay partly around and mess things up.
brainvisa/morpho-deepsulci#8
I'm fixing this, but the binary distributions in use will suffer from this problem...

[wsnyder4]

In the meantime, this would mean that rerunning the pipeline but saving a labeled graph to a separate location would yield accurate data then, correct?

I very much appreciate you looking into this!

[denisri]

It's fixed but we will need to build a new binary package. We can't do that before september, so you will need to erase any previous CNN labeled graph before re-running it, or install your distrib read-write (using the bv tool without arguments) and replace the file with the patch found here:
https://brainvisa.info/download/updates/5.0.4/python/deepsulci/sulci_labeling/capsul/labeling.py

[wsnyder4]

This will work just well! I appreciate the quick turnaround with this patch, and will place it in my read-write distribution.

Have a great holiday!

[denisri]

I am wondering if I could store the values directly from the C++ code. Is there any way I could manipulate this from my own BrainVISA install to obtain the verified values used to compute mean depth?

What exactly do you want to do ? Modify the C++ code to store individual depth values for each bottom voxel ? That would need to modify the C++ source code and recompile it. It is possible as a developer install, but not on a binary "user" distribution. For this, look at https://dev.brainvisa.info, but I don't think this is necessary at this time.

Note that I'm on vacation this evening, for a month, so I won't be able to help you during the next month...