Hi Nicole,

here's an explanation of the variables (and another code snippit). Let me know if you need more. I have plenty more snippits :).

doug

/*---------------------------------------------------------- MRIxfmCRS2XYZ() - computes the matrix needed to compute the XYZ of the center of a voxel at a given Col, Row, and Slice from the native geometry of the volume

x col y = T * row z slice 1 1

T = [Mdc*D Pxyz0] [0 0 0 1 ]

Mdc = [Vcol Vrow Vslice] V<dim> = the direction cosine pointing from the center of one voxel to the center of an adjacent voxel in the next dim, where dim is either colum, row, or slice. Vcol = [x_r x_a x_s], Vrow = [y_r y_a y_s], Vslice = [z_r z_a z_s]. Vcol can also be described as the vector normal to the plane formed by the rows and slices of a given column (ie, the column normal).

D = diag([colres rowres sliceres]) dimres = the distance between adjacent dim, where colres = mri->xsize, rowres = mri->ysize, and sliceres = mri->zsize.

Pxyz0 = the XYZ location at CRS=0. This number is not part of the mri structure, so it is computed here according to the formula: Pxyz0 = PxyzCenter - Mdc*D*PcrsCenter

PcrsCenter = the col, row, and slice at the center of the volume, = [ ncols/2 nrows/2 nslices/2 ]

PxyzCenter = the X, Y, and Z at the center of the volume and does exist in the header as mri->c_r, mri->c_a, and mri->c_s, respectively.

Note: to compute the matrix with respect to the first voxel being at CRS 1,1,1 instead of 0,0,0, then set base = 1. This is necessary with SPM matrices.

See also: MRIxfmCRS2XYZtkreg, MRItkReg2Native ------------------------------------------------------*/ MATRIX *MRIxfmCRS2XYZ(MRI *mri, int base) { MATRIX *m; MATRIX *Pcrs, *PxyzOffset;

m = MatrixAlloc(4, 4, MATRIX_REAL);

/* direction cosine between columns scaled by distance between colums */ *MATRIX_RELT(m, 1, 1) = mri->x_r * mri->xsize; *MATRIX_RELT(m, 2, 1) = mri->x_a * mri->xsize; *MATRIX_RELT(m, 3, 1) = mri->x_s * mri->xsize;

/* direction cosine between rows scaled by distance between rows */ *MATRIX_RELT(m, 1, 2) = mri->y_r * mri->ysize; *MATRIX_RELT(m, 2, 2) = mri->y_a * mri->ysize; *MATRIX_RELT(m, 3, 2) = mri->y_s * mri->ysize;

/* direction cosine between slices scaled by distance between slices */ *MATRIX_RELT(m, 1, 3) = mri->z_r * mri->zsize; *MATRIX_RELT(m, 2, 3) = mri->z_a * mri->zsize; *MATRIX_RELT(m, 3, 3) = mri->z_s * mri->zsize;

/* Preset the offsets to 0 */ *MATRIX_RELT(m, 1, 4) = 0.0; *MATRIX_RELT(m, 2, 4) = 0.0; *MATRIX_RELT(m, 3, 4) = 0.0;

/* Last row of matrix */ *MATRIX_RELT(m, 4, 1) = 0.0; *MATRIX_RELT(m, 4, 2) = 0.0; *MATRIX_RELT(m, 4, 3) = 0.0; *MATRIX_RELT(m, 4, 4) = 1.0;

/* At this point, m = Mdc * D */

/* Col, Row, Slice at the Center of the Volume */ Pcrs = MatrixAlloc(4, 1, MATRIX_REAL); *MATRIX_RELT(Pcrs, 1, 1) = mri->width/2.0 + base; *MATRIX_RELT(Pcrs, 2, 1) = mri->height/2.0 + base; *MATRIX_RELT(Pcrs, 3, 1) = mri->depth/2.0 + base; *MATRIX_RELT(Pcrs, 4, 1) = 1.0;

/* XYZ offset the first Col, Row, and Slice from Center */ /* PxyzOffset = Mdc*D*PcrsCenter */ PxyzOffset = MatrixMultiply(m,Pcrs,NULL);

/* XYZ at the Center of the Volume is mri->c_r, c_a, c_s */

/* The location of the center of the voxel at CRS = (0,0,0)*/ *MATRIX_RELT(m, 1, 4) = mri->c_r - PxyzOffset->rptr[1][1]; *MATRIX_RELT(m, 2, 4) = mri->c_a - PxyzOffset->rptr[2][1]; *MATRIX_RELT(m, 3, 4) = mri->c_s - PxyzOffset->rptr[3][1];

MatrixFree(&Pcrs); MatrixFree(&PxyzOffset);

return(m); }

Nicole Aucoin wrote:

Hi Doug,

At this point, my main question is: what is the meaning of the normal_ras variables, as saved to the RAS matrix as z_ras? Is it a unit vector that defines the depth direction of the volume?

My goal is to figure out the coordinates in RAS space of the top left corner of the last slice - the explicitly defined points in the bhdr file are all on the first slice.

Nicole

On Thu, 7 Apr 2005, Doug Greve wrote:

...

Hi Nicole,

here's the function we use to read the the .bhdr file. Let me know if this does not answer your questions.

doug

/*------------------------------------------------------*/ int read_bhdr(MRI *mri, FILE *fp) {

char line[STRLEN]; char *l; float tlr=0.; float tla=0.; float tls=0.; /* top left coordinates */ float trr=0.; float tra=0.; float trs=0.; /* top right coordinates */ float brr=0.; float bra=0.; float brs=0.; /* bottom right coordinates */ float xr=0.; float xa=0.; float xs=0.; float yr=0.; float ya=0.; float ys=0.; MATRIX *T, *CRSCenter, *RASCenter;

while(1){ // don't use "while (!feof(fp))"

/* --- read the line --- */
fgets(line, STRLEN, fp);

if (feof(fp))  // wow, it was beyound the end of the file. get out.
  break;

/* --- remove the newline --- */
if(line[strlen(line)-1] == '\n')
  line[strlen(line)-1] = '\0';

/* --- skip the initial spaces --- */
for(l = line;isspace((int)(*l));l++);

/* --- get the varible name and value(s) --- */
if(strlen(l) > 0){
  if(strncmp(l, "cols: ", 6) == 0)
    sscanf(l, "%*s %d", &mri->width);
  else if(strncmp(l, "rows: ", 6) == 0)
    sscanf(l, "%*s %d", &mri->height);
  else if(strncmp(l, "nslices: ", 9) == 0)
    sscanf(l, "%*s %d", &mri->depth);
  else if(strncmp(l, "n_time_points: ", 15) == 0)
    sscanf(l, "%*s %d", &mri->nframes);
  else if(strncmp(l, "slice_thick: ", 13) == 0)
    sscanf(l, "%*s %f", &mri->zsize);
  else if(strncmp(l, "image_te: ", 10) == 0)
    sscanf(l, "%*s %f", &mri->te);
  else if(strncmp(l, "image_tr: ", 10) == 0)
    sscanf(l, "%*s %f", &mri->tr);
  else if(strncmp(l, "image_ti: ", 10) == 0)
    sscanf(l, "%*s %f", &mri->ti);
  else if(strncmp(l, "flip_angle: ", 10) == 0)
    sscanf(l, "%*s %lf", &mri->flip_angle);
  else if(strncmp(l, "top_left_r: ", 12) == 0)
    sscanf(l, "%*s %g", &tlr);
  else if(strncmp(l, "top_left_a: ", 12) == 0)
    sscanf(l, "%*s %g", &tla);
  else if(strncmp(l, "top_left_s: ", 12) == 0)
    sscanf(l, "%*s %g", &tls);
  else if(strncmp(l, "top_right_r: ", 13) == 0)
    sscanf(l, "%*s %g", &trr);
  else if(strncmp(l, "top_right_a: ", 13) == 0)
    sscanf(l, "%*s %g", &tra);
  else if(strncmp(l, "top_right_s: ", 13) == 0)
    sscanf(l, "%*s %g", &trs);
  else if(strncmp(l, "bottom_right_r: ", 16) == 0)
    sscanf(l, "%*s %g", &brr);
  else if(strncmp(l, "bottom_right_a: ", 16) == 0)
    sscanf(l, "%*s %g", &bra);
  else if(strncmp(l, "bottom_right_s: ", 16) == 0)
    sscanf(l, "%*s %g", &brs);
  else if(strncmp(l, "normal_r: ", 10) == 0)
    sscanf(l, "%*s %g", &mri->z_r);
  else if(strncmp(l, "normal_a: ", 10) == 0)
    sscanf(l, "%*s %g", &mri->z_a);
  else if(strncmp(l, "normal_s: ", 10) == 0)
    sscanf(l, "%*s %g", &mri->z_s);
  else { /* --- ignore it --- */ }
}

} /* end while(!feof()) */ // forget to close file fclose(fp);

// vox to ras matrix is // // Xr*Sx Yr*Sy Zr*Sz Tr // Xa*Sx Ya*Sy Za*Sz Ta // Xs*Sx Ys*Sy Zs*Sz Ts // 0 0 0 1 // // Therefore // // trr = Xr*Sx*W + Tr, tlr = Tr // tra = Xa*Sx*W + Ta, tla = Ta // trs = Xs*Sx*W + Ts, tls = Ts // // Thus // // Sx = sqrt ( ((trr-tlr)/W)^2 + ((tra-tla)/W)^2 + ((trs-tls)/W)^2) // since Xr^2 + Xa^2 + Xs^2 = 1 // Xr = (trr-tlr)/(W*Sx) // Xa = (tra-tla)/(W*Sx) // Xs = (trs-tls)/(W*Sx) // // Similar things for others // xr = (trr - tlr) / mri->width; xa = (tra - tla) / mri->width; xs = (trs - tls) / mri->width; mri->xsize = sqrt(xr*xr + xa*xa + xs*xs); if (mri->xsize) // avoid nan { mri->x_r = xr / mri->xsize; mri->x_a = xa / mri->xsize; mri->x_s = xs / mri->xsize; } else // fake values { mri->xsize = 1; mri->x_r = -1; mri->x_a = 0; mri->x_s = 0; } yr = (brr - trr) / mri->height; ya = (bra - tra) / mri->height; ys = (brs - trs) / mri->height; mri->ysize = sqrt(yr*yr + ya*ya + ys*ys); if (mri->ysize) // avoid nan { mri->y_r = yr / mri->ysize; mri->y_a = ya / mri->ysize; mri->y_s = ys / mri->ysize; } else // fake values { mri->ysize = 1; mri->y_r = 0; mri->y_a = 0; mri->y_s = -1; } T = MRIxfmCRS2XYZ(mri,0);

T->rptr[1][4] = tlr; T->rptr[2][4] = tla; T->rptr[3][4] = tls;

//printf("------- read_bhdr: T ---------\n"); //MatrixPrint(stdout,T); //printf("------------------------------\n");

CRSCenter = MatrixAlloc(4,1,MATRIX_REAL); CRSCenter->rptr[1][1] = (mri->width)/2.0; CRSCenter->rptr[2][1] = (mri->height)/2.0; CRSCenter->rptr[3][1] = (mri->depth)/2.0; CRSCenter->rptr[4][1] = 1;

RASCenter = MatrixMultiply(T,CRSCenter,NULL); mri->c_r = RASCenter->rptr[1][1]; mri->c_a = RASCenter->rptr[2][1]; mri->c_s = RASCenter->rptr[3][1];

MatrixFree(&T); MatrixFree(&CRSCenter); MatrixFree(&RASCenter);

#if 0 /* This computation of the center is incorrect because TL is the center of the first voxel (not the edge of the FOV). TR and BR are actually outside of the FOV. */ mri->c_r = (brr + tlr) / 2.0 + (mri->depth - 1) * mri->z_r * mri->zsize / 2.0 ; mri->c_a = (bra + tla) / 2.0 + (mri->depth - 1) * mri->z_a * mri->zsize / 2.0 ; mri->c_s = (brs + tls) / 2.0 + (mri->depth - 1) * mri->z_s * mri->zsize / 2.0 ; #endif

mri->ras_good_flag = 1;

mri->thick = mri->zsize; mri->ps = mri->xsize;

return(NO_ERROR);

} /* end read_bhdr() */

Nicole Aucoin wrote:

...

Hello,

I'm trying to orient bfloat and bshort volumes in 3D Slicer, and to do that, need to calculate the corners of the volumes to calculate an RAS to IJK transform matrix. The .bhdr files that I have contain the top left, top right and bottom right corners, which I'm using as is with the assumption that they are on the first slice, and are in relation to an origin of (0,0,0). I also need a point on the last slice, so I'm using those corners to create vectors that will give me a cross product pointing to a corner on the last slice. I then scale that vector by slice_thick * nslices and add it to the top left corner point. It's almost but not quite right at this point, so I'd like to check my assumptions:

• are the top_left_?, top_right_?, bottom_right_? points on the first

slice?

• is the origin (0,0,0)? if not, how do I find it?

What do the .bhdr file entries normal_? and image_t? denote?

thanks, Nicole