bet.c

mri_toolbox是一个工具用来MRI. 来自于SourceForge, 我上传这个软件,希望能结识对医疗软件感兴趣的兄弟.

/* {{{ Copyright etc. */

/*  bet.c - Brain Extraction Tool

    Stephen Smith, FMRIB Image Analysis Group

    Copyright (C) 1999-2003 University of Oxford  */

/*  CCOPYRIGHT */

/* }}} */
/* {{{ defines, includes and typedefs */

#include "libss/libss.h"
#include "libss/libavw.h"
#include "libss/libtessa.h"

void usage(void);

/* }}} */
/* {{{ usage */

void usage(void)
{
  printf("\nBET (Brain Extraction Tool) v1.2 - FMRIB Analysis Group, Oxford\n\n");
  printf("Usage: bet <input fileroot> <output fileroot> [options]\n\n");

  printf("-o : generate brain surface outline overlaid onto original image\n");
  printf("-m : generate binary brain mask\n");
  printf("-s : generate approximate skull image\n");
  printf("-n : don't generate segmented brain image output\n");
  printf("-f <fractional_threshold> : fractional intensity threshold (0->1); default=0.5; smaller values give larger brain outline estimates\n");
  printf("-g <fractional_threshold> : vertical gradient in fractional intensity threshold (-1->1); default=0; positive values give larger brain outline at bottom, smaller at top\n");
  printf("-t : apply thresholding to segmented brain image and mask\n");
  printf("-c <x y z> : co-ordinates (mm not voxels) for centre of initial brain surface sphere\n");
  printf("-r <r>: head radius (mm not voxels); initial surface sphere is set to half of this\n");
  printf("-v : verbose text output\n\n");

  /*printf("-X : generate xtopol format output; file extensions=coo,dat\n");
  printf("-C : generate approximate image cost function output; fileroot=<output filroot>_cost\n");
  printf("-S : as above but \"colour code\" brain-skull distance\n");*/

  exit(1);
}

/* }}} */
/* {{{ main */

#define TESSELATE_ORDER             5       /* 5 */
#define ITERATIONS                  1000    /* 1000 */
#define BRAIN_THRESHOLD_DEFAULT     0.5     /* 0.5 */
#define COST_SEARCH                 7       /* 7 mm */
#define NORMAL_MAX_UPDATE_FRACTION  0.5     /* 0.5 */
#define LAMBDA_TANGENT              0.5     /* 0.5 */
#define LAMBDA_FIT                  0.1     /* 0.1 */
#define RADIUSMAX                   10.0    /* 10.0 mm */
#define RADIUSMIN                   3.33    /* 3.33 mm */
#define SELF_INTERSECTION_THRESHOLD 4000    /* 4000 */
#define SKULL_SEARCH                30      /* 30 mm */
#define SKULL_START                 -3      /* -3 mm */
#define MIN_FOV                     20      /* 20 mm */

/*#define DEBUG_NORMALS*/
/*#define DEBUG_MOVEMENT*/
/*#define DEBUG_EVOLVE*/

int main(argc, argv)
  int   argc;
  char  *argv [];
{
  /* {{{ vars */

FDT *in, *mask=NULL, *raw=NULL, threshold, thresh2, thresh98,
  hist_min=0, hist_max=0, medianval;
int x_size, y_size, z_size, t_size, x, y, z, i, pc=0, iters, pass=1,
  output_brain=1, output_xtopol=0, output_cost=0, output_mask=0,
  output_overlay=0, output_skull=0, apply_thresholding=0, code_skull=0, verbose=0;
double cgx=-1000000, cgy, cgz, radius=-1000000, scale, ml=0, ml0=0, tmpf,
  brain_threshold=BRAIN_THRESHOLD_DEFAULT, gradthresh=0, incfactor=0,
  rE = 0.5 * (1/RADIUSMIN + 1/RADIUSMAX), rF = 6 / (1/RADIUSMIN - 1/RADIUSMAX);
char filename[1000];
image_struct im;
points_struc *v;

/* }}} */

  /* {{{ process arguments */

if (argc<3)
     usage();

avw_read(argv[1],&im);

in=im.i;
x_size=im.x; y_size=im.y; z_size=im.z;
t_size=MAX(im.t,1); im.t=1;
scale=MIN(im.xv,MIN(im.yv,im.zv)); /* minimum voxel side length */

for (i = 3; i < argc; i++) {
  if (!strcmp(argv[i], "-X"))
    output_xtopol=1;
  else if (!strcmp(argv[i], "-C"))
    output_cost=1;
  else if (!strcmp(argv[i], "-s"))
    output_skull=1;
  else if (!strcmp(argv[i], "-S"))
    /* {{{ colour-coded skull image */

    {
      output_skull=1;
      code_skull=1;
    }

/* }}} */
  else if (!strcmp(argv[i], "-m"))
    output_mask=1;
  else if (!strcmp(argv[i], "-o"))
    output_overlay=1;
  else if (!strcmp(argv[i], "-n"))
    output_brain=0;
  else if (!strcmp(argv[i], "-t"))
    apply_thresholding=1;
  else if (!strcmp(argv[i], "-v"))
    verbose=1;
  else if (!strcmp(argv[i], "-f"))
    /* {{{ fractional brain threshold */

    {
      i++;
      if (argc<i+1) /* option following f hasn't been given */
      {
	printf("Error: no value given following -f\n");
	usage();
      }
      brain_threshold=atof(argv[i]);
      if ( (brain_threshold<=0) || (brain_threshold>=1) )
      {
	printf("Error: value following -f must lie between 0 and 1\n");
	usage();
      }
    }

/* }}} */
  else if (!strcmp(argv[i], "-g"))
    /* {{{ gradient fractional brain threshold */

    {
      i++;
      if (argc<i+1) /* option following g hasn't been given */
      {
	printf("Error: no value given following -g\n");
	usage();
      }
      gradthresh=atof(argv[i]);
      if ( (gradthresh<-1) || (gradthresh>1) )
      {
	printf("Error: value following -g must lie between -1 and 1\n");
	usage();
      }
    }

/* }}} */
  else if (!strcmp(argv[i], "-c"))
    /* {{{ initial sphere centre */

    {
      i++;
      if (argc<i+3) /* options following -c hasn't been given */
      {
	printf("Error: need three values after -c\n");
	usage();
      }
      cgx=atof(argv[i]);
      i++;
      cgy=atof(argv[i]);
      i++;
      cgz=atof(argv[i]);
    }

/* }}} */
  else if (!strcmp(argv[i], "-r"))
    /* {{{ initial sphere radius */

    {
      i++;
      if (argc<i+1) /* option following -r hasn't been given */
      {
	printf("Error: no value given following -r\n");
	usage();
      }
      radius=atof(argv[i]);
      if (radius<0)
      {
	printf("Error: value following -r must be greater than 0\n");
	usage();
      }
    }

/* }}} */
  else
    usage();
}

brain_threshold=pow(brain_threshold,0.275);

if ( !output_xtopol && !output_cost && !output_skull && !output_mask && !output_overlay && !output_brain )
{
  printf("No outputs requested!\n");
  usage();
}

/* }}} */
  /* {{{ image preprocessing */

/* {{{ setup for single-slice data if present */

if (im.x*im.xv<MIN_FOV) im.xv=200;
if (im.y*im.yv<MIN_FOV) im.yv=200;
if (im.z*im.zv<MIN_FOV) im.zv=200;

/* }}} */

if (verbose) printf("\nBET (Brain Extraction Tool) v1.2 - FMRIB Analysis Group, Oxford\n\n");

im.min=im.max=0;
find_thresholds(&im,0.1);
hist_min=im.min; hist_max=im.max; thresh2=im.thresh2; thresh98=im.thresh98; threshold=im.thresh;
if (verbose) printf("hist_min=%.2f thresh2=%.2f thresh=%.2f thresh98=%.2f hist_max=%.2f\n",
       im.min,(double)thresh2,(double)threshold,(double)thresh98,im.max);
if (verbose) printf("THRESHOLD %.1f\n",(double)threshold);

if (cgx<-999999)
{
  c_of_g (im,&cgx,&cgy,&cgz);
  cgx*=im.xv; cgy*=im.yv; cgz*=im.zv;
}
if (verbose) printf("CofG (%.1f,%.1f,%.1f) mm\n",cgx,cgy,cgz);

if (radius<-999999)
     radius = find_radius (im,im.xv*im.yv*im.zv);
if (verbose) printf("RADIUS %.1f\n",radius);

if ( im.thresh98 - im.thresh2 > 1.5 )
{
  FDT *tmpimage = (FDT *) malloc(sizeof(FDT)*x_size*y_size*z_size);

  i=0;
  for(z=0; z<z_size; z++)
    for(y=0; y<y_size; y++)
      for(x=0; x<x_size; x++)
	{
	  FDT tmp=IA(in,x,y,z);
	  
	  if ( (tmp>thresh2) && (tmp<thresh98) &&
	       ( (x*im.xv-cgx)*(x*im.xv-cgx) + (y*im.yv-cgy)*(y*im.yv-cgy) + (z*im.zv-cgz)*(z*im.zv-cgz) < radius*radius ) )
	    tmpimage[i++]=tmp;
	}
  medianval = median(0.5,tmpimage,i);
  if (verbose)   printf("MEDIANVAL %.1f\n",(double)medianval);
  free(tmpimage);
} else {
  medianval=im.thresh98;
  if (verbose)   printf("MEDIANVAL %.1f\n",(double)medianval);
}

if (output_cost) /* prepare cost function image for writing into */
{
  raw = (FDT *) malloc(sizeof(FDT)*x_size*y_size*z_size);
  memset(raw,(unsigned char)0,sizeof(FDT)*x_size*y_size*z_size);
}

/* }}} */

  while (pass>0)
    {
      /* {{{ initialize tessellation */

if (pass==1)
{
  object *old = &ico;                /* start with icosohedral tessellation */
  
  tessa((int)TESSELATE_ORDER,&old);  /* create tessellated triangles; level 5 gives 2562 points */
  pc=points_list(old,&v);            /* convert triangles to vertex list */
  free(old);                         /* don't need triangular tessellation any more */

  /*printf("VERTICES %d\n",pc);*/

  /* measure initial spacing for use later in self-intersection calculations */
  ml0 = sqrt( (v[0].xorig-v[v[0].n[0]].xorig)*(v[0].xorig-v[v[0].n[0]].xorig) +
	      (v[0].yorig-v[v[0].n[0]].yorig)*(v[0].yorig-v[v[0].n[0]].yorig) +
	      (v[0].zorig-v[v[0].n[0]].zorig)*(v[0].zorig-v[v[0].n[0]].zorig) );
  /*printf("ml0=%f\n",ml0);*/

#ifdef DEBUG_NORMALS
  v = (points_struc *)realloc((void *)v,sizeof(points_struc)*2*(pc+10)); /* make space for storing normals */
#endif

}

/* scale vertex positions for this image, set surface small and allow to grow */
for(i=0; i<pc; i++)
{
  v[i].x = v[i].xorig * radius * 0.5 + cgx;
  v[i].y = v[i].yorig * radius * 0.5 + cgy;
  v[i].z = v[i].zorig * radius * 0.5 + cgz;
}

/* }}} */
      /* {{{ find brain surface */

for(iters=0; iters<ITERATIONS; iters++)
{
  /* {{{ find local surface normals */

for(i=0; i<pc; i++)
{
  double nx, ny, nz, tmpf;
  int k, l;

  nx=ny=nz=0.0;

  for(k=0; v[i].n[k]>-1; k++); /* find number of connections */

  for(l=0; l<k; l++) /* for each pair of consecutive neighbours form a vector product to get normal */
    {
      double adx = v[v[i].n[l]].x - v[i].x,
	ady = v[v[i].n[l]].y - v[i].y,
	adz = v[v[i].n[l]].z - v[i].z,
	bdx = v[v[i].n[(l+1)%k]].x - v[i].x,
	bdy = v[v[i].n[(l+1)%k]].y - v[i].y,
	bdz = v[v[i].n[(l+1)%k]].z - v[i].z;
  
      nx += ady*bdz - adz*bdy;
      ny += adz*bdx - adx*bdz;
      nz += adx*bdy - ady*bdx;
    }

  /* make the normal vector of length 1 */
  tmpf = sqrt(nx*nx+ny*ny+nz*nz);
  nx/=(double)tmpf; ny/=(double)tmpf; nz/=(double)tmpf;

  /* {{{ debug normals */

#ifdef DEBUG_NORMALS

if (iters==ITERATIONS-1) /* final iteration */
{
     v[pc+i].x=v[i].x+nx*20.0;
     v[pc+i].y=v[i].y+ny*20.0;
     v[pc+i].z=v[i].z+nz*20.0;
     v[pc+i].n[0]=i;
     v[pc+i].n[1]=-1;
}

#endif

/* }}} */

  v[i].nx=nx;
  v[i].ny=ny;
  v[i].nz=nz;
}

/* }}} */
  /* {{{ find mean connection length every now and then */

if ( (iters==50) || (iters%100==0) ) /* add the 50 as the rate of change is highest at start; thus do 0,50,100,200,..... */
{
  int l;

  ml=0;

  for(i=0; i<pc; i++)
    {
      double mml=0;
      
      for(l=0; v[i].n[l]>-1; l++)
	mml += sqrt( (v[i].x-v[v[i].n[l]].x)*(v[i].x-v[v[i].n[l]].x) +
		     (v[i].y-v[v[i].n[l]].y)*(v[i].y-v[v[i].n[l]].y) +
		     (v[i].z-v[v[i].n[l]].z)*(v[i].z-v[v[i].n[l]].z) );

      ml += mml/l;
    }

  ml /= pc;
  /*printf("iteration=%d, ml=%f\n",iters,ml);*/
}

/* }}} */
  /* {{{ increased smoothing for pass>1 */

if (pass>1)
{
  incfactor=pow((double)10.0,(double)pass);

  if (iters>ITERATIONS*0.75)
    incfactor=4.0*(1.0-((double)iters)/ITERATIONS)*(incfactor-1.0) + 1.0;
}

/* }}} */

  for(i=0; i<pc; i++) /* calculate tessellation update */
    {
      /* {{{ variables, and setup k and normal */

FDT lmin, lmax;
int d, k, l;
double nx=v[i].nx, ny=v[i].ny, nz=v[i].nz, sx, sy, sz, fit, sn, stx, sty, stz;

for(k=0; v[i].n[k]>-1; k++); /* find number of connections */

/* }}} */
      /* {{{ average position of neighbours: smoothing vector */

/* s is vector from current vertex to the mean position of its neighbours */

sx=sy=sz=0;

for(l=0; l<k; l++)
{
  sx += v[v[i].n[l]].x;
  sy += v[v[i].n[l]].y;
  sz += v[v[i].n[l]].z;
}

sx = sx/k - v[i].x;
sy = sy/k - v[i].y;
sz = sz/k - v[i].z;

/* part of s normal to surface, sn = n * (s.n)
   part of s tangent to surface, st = s - sn */

sn = sx*nx + sy*ny + sz*nz; /* this is just the s.n part - will multiply by n later */

stx = sx - nx * sn;
sty = sy - ny * sn;
stz = sz - nz * sn;

/* }}} */
      /* {{{ COMMENT ORIG find intensity-based part of cost function - local max */

#ifdef FoldingComment

{
  FDT lnew;
  int all_inside_image=1;
  double lthresh, local_brain_threshold;

  fit=0;
  lmin=hist_max;
  lmax=thresh98/2+thresh2/2;

  tmpf = brain_threshold + gradthresh * ( (v[i].z-cgz) / radius );
  local_brain_threshold = MIN( 1.0 , MAX( tmpf , 0.0 ) );  

  d=1;   /* start d at 1 not 0 so that boundary is just outside brain not on the actual edge */
  x=FTOI((v[i].x-((double)d)*nx)/im.xv); y=FTOI((v[i].y-((double)d)*ny)/im.yv); z=FTOI((v[i].z-((double)d)*nz)/im.zv);
  if ( (x>=0) && (x<x_size) && (y>=0) && (y<y_size) && (z>=0) && (z<z_size) )
    {
      lnew=IA(in,x,y,z);
      lmin = MIN(lmin,lnew);
      lmax = MAX(lmax,lnew);
    }
  else all_inside_image=0;

  d=COST_SEARCH;
  x=FTOI((v[i].x-((double)d)*nx)/im.xv); y=FTOI((v[i].y-((double)d)*ny)/im.yv); z=FTOI((v[i].z-((double)d)*nz)/im.zv);
  if ( (x>=0) && (x<x_size) && (y>=0) && (y<y_size) && (z>=0) && (z<z_size) )
    {
      lnew=IA(in,x,y,z);
      lmin = MIN(lmin,lnew);
      /*lmax = MAX(lmax,lnew);*/
    }
  else all_inside_image=0;

  if (all_inside_image)
    {
      for(d=2;d<COST_SEARCH;d++)
	{
	  x=FTOI((v[i].x-((double)d)*nx)/im.xv); y=FTOI((v[i].y-((double)d)*ny)/im.yv); z=FTOI((v[i].z-((double)d)*nz)/im.zv);
	  lnew=IA(in,x,y,z);
	  lmin = MIN(lmin,lnew);
	  if (d<COST_SEARCH/2) /* only look relatively locally for maximum intensity */
	    lmax = MAX(lmax,lnew);
	}

      lmin = MAX(lmin,thresh2);   /* so that extreme values don't screw this up */
      lmax = MIN(lmax,thresh98);

      lthresh = (lmax - thresh2)*local_brain_threshold + thresh2;
      tmpf = lmin - lthresh;
      fit = tmpf / ((lmax - thresh2)*0.5); /* scale range to around -1:1 */
  
      if (output_cost)
	{
	  x=(v[i].x/im.xv); y=(v[i].y/im.yv); z=(v[i].z/im.zv);
	  if ( (x>=0) && (x<x_size) && (y>=0) && (y<y_size) &&(z>=0) && (z<z_size) )
	    IA(raw,x,y,z)=thresh98/2 + thresh2/2 + (0.5*fit*((double)thresh98 - (double)thresh2));
	}  
    }
}

#endif

/* }}} */
      /* {{{ find intensity-based part of cost function - local max */

{
  FDT lnew;
  int all_inside_image=1;
  double lthresh, local_brain_threshold=brain_threshold;

  fit=0;
  lmin=medianval;
  lmax=threshold;

  /* {{{ change local threshold if gradient threshold used */

  if (gradthresh!=0)
    {
      tmpf = brain_threshold + gradthresh * ( (v[i].z-cgz) / radius );
      local_brain_threshold = MIN( 1.0 , MAX( tmpf , 0.0 ) );
    }

/* }}} */

  d=1;   /* start d at 1 not 0 so that boundary is just outside brain not on the actual edge */
  x=FTOI((v[i].x-((double)d)*nx)/im.xv); y=FTOI((v[i].y-((double)d)*ny)/im.yv); z=FTOI((v[i].z-((double)d)*nz)/im.zv);
  if ( (x>=0) && (x<x_size) && (y>=0) && (y<y_size) && (z>=0) && (z<z_size) )
    {
      lnew=IA(in,x,y,z);
      lmin = MIN(lmin,lnew);
      lmax = MAX(lmax,lnew);
    }
  else all_inside_image=0;

  d=COST_SEARCH;
  x=FTOI((v[i].x-((double)d)*nx)/im.xv); y=FTOI((v[i].y-((double)d)*ny)/im.yv); z=FTOI((v[i].z-((double)d)*nz)/im.zv);
  if ( (x>=0) && (x<x_size) && (y>=0) && (y<y_size) && (z>=0) && (z<z_size) )
    {
      lnew=IA(in,x,y,z);
      lmin = MIN(lmin,lnew);
    }
  else all_inside_image=0;

  if (all_inside_image)
    {
      for(d=2;d<COST_SEARCH;d++)
	{
	  x=FTOI((v[i].x-((double)d)*nx)/im.xv); y=FTOI((v[i].y-((double)d)*ny)/im.yv); z=FTOI((v[i].z-((double)d)*nz)/im.zv);
	  lnew=IA(in,x,y,z);
	  lmin = MIN(lmin,lnew);
	  if (d<COST_SEARCH/2) /* only look relatively locally for maximum intensity */
	    lmax = MAX(lmax,lnew);
	}

      lmin = MAX(lmin,thresh2);
      lmax = MIN(lmax,medianval);

      lthresh = (lmax - thresh2)*local_brain_threshold + thresh2;

      tmpf = lmin - lthresh;
      if (lmax - thresh2>0)
	fit = tmpf / ((lmax - thresh2)*0.5); /* scale range to around -1:1 */
      else
	fit = tmpf / 0.5;
  
      if (output_cost)
	{
	  x=(v[i].x/im.xv); y=(v[i].y/im.yv); z=(v[i].z/im.zv);
	  if ( (x>=0) && (x<x_size) && (y>=0) && (y<y_size) &&(z>=0) && (z<z_size) )
	    IA(raw,x,y,z)=thresh98/2 + thresh2/2 + (0.5*fit*((double)thresh98 - (double)thresh2));
	}  
    }
}

/* }}} */
      /* {{{ estimate the update */

/* normal component of smoothing */
tmpf = 2 * ABS(sn) / (ml*ml);              /* 1/r ; r=local radius of curvature */
tmpf = 0.5 * ( 1 + tanh((tmpf-rE)*rF) );   /* f(1/r) ; this makes big r have low correction and vice versa */

if ( (pass>1) && (sn > 0) )  /* if need increased smoothing, to avoid self-intersection; */
{                            /* only apply to concave curvature (as seen from outside surface) */
  tmpf *= incfactor;
  tmpf = MIN(tmpf,1);
}

tmpf = sn*tmpf;                                              /* multiply normal component by correction factor */

tmpf += NORMAL_MAX_UPDATE_FRACTION * ml * LAMBDA_FIT * fit ; /* combine normal smoothing and image fit */

v[i].xnew = v[i].x + stx*LAMBDA_TANGENT + tmpf*nx;
v[i].ynew = v[i].y + sty*LAMBDA_TANGENT + tmpf*ny;
v[i].znew = v[i].z + stz*LAMBDA_TANGENT + tmpf*nz;

/* }}} */
    }

  /* {{{ debug evolve: output single slice to show how shape evolves */

#ifdef DEBUG_EVOLVE

#define SPACING 5

if (iters%SPACING==0)
{
  image_struct imd=im;
  FILE *ofp;
  FDT *tmpimage=(FDT *)calloc(x_size*y_size*z_size,sizeof(FDT));

  if (iters==0)
    {
      imd.x=im.y; imd.y=im.z; imd.z=ITERATIONS/SPACING; imd.t=1;
      imd.xv0=im.yv; imd.yv0=im.zv; imd.zv0=1;
      imd.min=thresh2; imd.max=thresh98;
      imd.i=tmpimage;

      sprintf(filename,"%s_debug",argv[2]);
      avw_write(filename,imd);