External Email - Use Caution

Hi Doug,

thanks a lot for the clarification!

Tim

On August 16, 2019 at 6:19 PM "Greve, Douglas N.,Ph.D." DGREVE@mgh.harvard.edu wrote:

The algorithm in v6 is not as you describe. The function below is based on Anderson Winkler's (cc'ed) work.

/*!   \fn MRISvolumeTH3(MRIS *w, MRIS *p, MRI *vol)   \brief Compute vertex-wise volume based on dividing each obliquely truncated   trilateral pyramid into three tetrahedra. Based on Anderson M. Winkler's   srf2vol matlab script.   \params w is white surface, p is pial. Output is an MRI struct of size nvertices  */ MRI *MRISvolumeTH3(MRIS *w, MRIS *p, MRI *vol, MRI *mask, double *totvol) {   MRI *facevol;   int nthface, v1, v2, v3;   int vno, nfacenbrs, nthfacenbr, nbrfaceno;   MATRIX *Ap, *Bp, *Cp, *Aw, *Bw, *Cw;   MATRIX *C1 = NULL, *C2 = NULL, *C3 = NULL;   double T1, T2, T3, volsum, dval, vtxvolsum;   VERTEX *vtx;

if (!vol) vol = MRIalloc(w->nvertices, 1, 1, MRI_FLOAT);

Ap = MatrixAlloc(3, 1, MATRIX_REAL);   Bp = MatrixAlloc(3, 1, MATRIX_REAL);   Cp = MatrixAlloc(3, 1, MATRIX_REAL);   Aw = MatrixAlloc(3, 1, MATRIX_REAL);   Bw = MatrixAlloc(3, 1, MATRIX_REAL);   Cw = MatrixAlloc(3, 1, MATRIX_REAL);

facevol = MRIalloc(w->nfaces, 1, 1, MRI_FLOAT);   volsum = 0.0;   for (nthface = 0; nthface < w->nfaces; nthface++) {     MRIsetVoxVal(facevol, nthface, 0, 0, 0, 0.0);     if (w->faces[nthface].ripflag) continue;

v1 = w->faces[nthface].v[0];     v2 = w->faces[nthface].v[1];     v3 = w->faces[nthface].v[2];

Ap->rptr[1][1] = p->vertices[v1].x;     Ap->rptr[2][1] = p->vertices[v1].y;     Ap->rptr[3][1] = p->vertices[v1].z;     Bp->rptr[1][1] = p->vertices[v2].x;     Bp->rptr[2][1] = p->vertices[v2].y;     Bp->rptr[3][1] = p->vertices[v2].z;     Cp->rptr[1][1] = p->vertices[v3].x;     Cp->rptr[2][1] = p->vertices[v3].y;     Cp->rptr[3][1] = p->vertices[v3].z;

Aw->rptr[1][1] = w->vertices[v1].x;     Aw->rptr[2][1] = w->vertices[v1].y;     Aw->rptr[3][1] = w->vertices[v1].z;     Bw->rptr[1][1] = w->vertices[v2].x;     Bw->rptr[2][1] = w->vertices[v2].y;     Bw->rptr[3][1] = w->vertices[v2].z;     Cw->rptr[1][1] = w->vertices[v3].x;     Cw->rptr[2][1] = w->vertices[v3].y;     Cw->rptr[3][1] = w->vertices[v3].z;

// As the Ap is the common vertex for all three, it can be used as the origin.     Bp = MatrixSubtract(Bp, Ap, Bp);     Cp = MatrixSubtract(Cp, Ap, Cp);     Aw = MatrixSubtract(Aw, Ap, Aw);     Bw = MatrixSubtract(Bw, Ap, Bw);     Cw = MatrixSubtract(Cw, Ap, Cw);

// The next lines compute the volume for each, using a scalar triple product:     C1 = VectorCrossProduct(Bw, Cw, C1);     T1 = fabs(V3_DOT(Aw, C1));

C2 = VectorCrossProduct(Cp, Bw, C2);     T2 = fabs(V3_DOT(Bp, C2));

C3 = VectorCrossProduct(Cw, Bw, C3);     T3 = fabs(V3_DOT(Cp, C3));

dval = (T1 + T2 + T3) / 6.0;     MRIsetVoxVal(facevol, nthface, 0, 0, 0, dval);     volsum += dval;   }   printf("Total face volume %g\n", volsum);

vtxvolsum = 0;   for (vno = 0; vno < w->nvertices; vno++) {     vtx = &w->vertices[vno];     if (vtx->ripflag) continue;     if (mask && MRIgetVoxVal(mask, vno, 0, 0, 0) < 0.5) continue;     nfacenbrs = w->vertices_topology[vno].num;     volsum = 0.0;     for (nthfacenbr = 0; nthfacenbr < nfacenbrs; nthfacenbr++) {       nbrfaceno = w->vertices_topology[vno].f[nthfacenbr];       volsum += (MRIgetVoxVal(facevol, nbrfaceno, 0, 0, 0) / 3.0);       // divide by 3 because each face participates in 3 vertices     }     MRIsetVoxVal(vol, vno, 0, 0, 0, volsum);     vtxvolsum += volsum;   }   printf("Total vertex volume %g (mask=%d)\n", vtxvolsum, mask == NULL);   *totvol = vtxvolsum;

MRIfree(&facevol);   MatrixFree(&Ap);   MatrixFree(&Bp);   MatrixFree(&Cp);   MatrixFree(&Aw);   MatrixFree(&Bw);   MatrixFree(&Cw);   MatrixFree(&C1);   MatrixFree(&C2);   MatrixFree(&C3);

return (vol); }

On 8/15/2019 4:15 AM, Tim Schäfer wrote:

...

     External Email - Use Caution

Dear FreeSurfer experts,

can anybody refer me to a publication that explains in detail how vertex-wise cortical volume (the data in <subject>/surf/?h.volume) is computed in FreeSurfer 6?

My understanding is the following:

• area is computed as 1/3 the area of all triangles meeting at that vertex.
• thickness is the shortest distance from the vertex (at white surface) to any vertex on pial, then the other way around from that pial vertex, divided by two.
• volume is area * thickness

Is that true, and where can I read up on it?

Thanks,

Tim

-- Dr. Tim Schäfer Postdoc Computational Neuroimaging Department of Child and Adolescent Psychiatry, Psychosomatics and Psychotherapy University Hospital Frankfurt, Goethe University Frankfurt am Main, Germany

---

Freesurfer mailing list Freesurfer@nmr.mgh.harvard.edu https://mail.nmr.mgh.harvard.edu/mailman/listinfo/freesurfer

---

Freesurfer mailing list Freesurfer@nmr.mgh.harvard.edu https://mail.nmr.mgh.harvard.edu/mailman/listinfo/freesurfer

-- Dr. Tim Schäfer Postdoc Computational Neuroimaging Department of Child and Adolescent Psychiatry, Psychosomatics and Psychotherapy University Hospital Frankfurt, Goethe University Frankfurt am Main, Germany