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
-- Nicole Aucoin nicole@bwh.harvard.edu Surgical Planning Laboratory, Brigham and Women's Hospital Rm 245, 1249 Boylston St., Boston, MA 02215 http://www.slicer.org _______________________________________________ Freesurfer mailing list Freesurfer@nmr.mgh.harvard.edu https://mail.nmr.mgh.harvard.edu/mailman/listinfo/freesurfer