bet.c

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

      sprintf(filename,"%s_debug.img",argv[2]);
      ofp=fopen(filename,"wb");
    }
  else
    {
      sprintf(filename,"%s_debug.img",argv[2]);
      ofp=fopen(filename,"ab");
    }

  memcpy(tmpimage,in,sizeof(FDT)*x_size*y_size*z_size);
  im.i=tmpimage;
  draw_surface(&im,hist_max,thresh2,v,pc);
  im.i=in;

  x=x_size*0.52;
  for(z=0;z<z_size;z++)
    for(y=0;y<y_size;y++)
      fwrite(&tmpimage[z*y_size*x_size+y*x_size+x],sizeof(FDT),1,ofp);

  free(tmpimage);
  fclose(ofp);
}

#endif

/* }}} */
  /* {{{ debug movement */

#ifdef DEBUG_MOVEMENT

if ( (iters==50) || (iters%100==0) )
{
  for(mm=0, i=0; i<pc; i++)
    mm += sqrt ( (v[i].x-v[i].xnew)*(v[i].x-v[i].xnew) + (v[i].y-v[i].ynew)*(v[i].y-v[i].ynew) + (v[i].z-v[i].znew)*(v[i].z-v[i].znew) );

  mm /= pc;

  printf("mean movement=%f\n",mm);
}

#endif

/* }}} */
  /* {{{ update tessellation */

for(i=0; i<pc; i++)
{
  v[i].x = v[i].xnew;
  v[i].y = v[i].ynew;
  v[i].z = v[i].znew;
}

/* }}} */
}

/* }}} */
      /* {{{ test surface for non-spherecity */

if (pass<10)
{
  int j, l;
  double intersection=0;

  for(i=0; i<pc; i++) /* loop round all points */
      for(j=0; j<pc; j++) /* inner loop round all points */
	{
	  int do_it=1;
	  
	  if (j==i) /* other point is same as current one - don't use */
	    do_it=0;
      
	  for(l=0; v[i].n[l]>-1; l++) /* other point is connected to current one - don't use */
	    if (j==v[i].n[l])
	      do_it=0;

	  if ( (do_it) &&
	       ( (v[i].x-v[j].x)*(v[i].x-v[j].x) + (v[i].y-v[j].y)*(v[i].y-v[j].y) + (v[i].z-v[j].z)*(v[i].z-v[j].z) < ml*ml ) )
	    {
	      double dist = sqrt ( (v[i].x-v[j].x)*(v[i].x-v[j].x) +
				   (v[i].y-v[j].y)*(v[i].y-v[j].y) +
				   (v[i].z-v[j].z)*(v[i].z-v[j].z) ),
		distorig = sqrt ( (v[i].xorig-v[j].xorig)*(v[i].xorig-v[j].xorig) +
				  (v[i].yorig-v[j].yorig)*(v[i].yorig-v[j].yorig) +
				  (v[i].zorig-v[j].zorig)*(v[i].zorig-v[j].zorig) );

	      tmpf = (distorig/ml0) - (dist/ml);    /* orig distance (in units of mean length) - current distance */
	      tmpf *= tmpf;                         /* weight higher values more highly */
	      /*printf("self-intersection value %f at: %f %f %f\n",tmpf,v[i].x/im.xv,v[i].y/im.yv,v[i].z/im.zv);*/
	      intersection += tmpf;
	    }
	}

  if (verbose) printf("self-intersection total = %.1f (threshold=%.1f)\n",intersection,(double)SELF_INTERSECTION_THRESHOLD);

  if (intersection>SELF_INTERSECTION_THRESHOLD)
    {
      if (verbose) printf("thus will rerun with higher smoothness constraint\n");
      pass++;
    }
  else
    pass=0;
}
else
  pass=0;

/* }}} */
    }

  /* {{{ write brain image and mask */

if ( (output_mask) || (output_brain) )
{
  mask = (FDT *) malloc(sizeof(FDT)*x_size*y_size*z_size);
  im.i=mask;
  fill_surface(&im,v,pc);

  if (apply_thresholding)
    for(i=0; i<z_size*y_size*x_size; i++)
      if (in[i]<threshold)
	mask[i]=0;
}

if (output_mask)
{
  im.min=0;
  im.max=1;
  sprintf(filename,"%s_mask",argv[2]);
  avw_write(filename,im);
}

if (output_brain)
{
  FDT *tmpin = malloc(sizeof(FDT)*x_size*y_size*z_size*t_size);
  int goesneg=0;

  im.min=thresh2;
  im.max=thresh98;

  if ( (im.dtmin<0) && ((double)hist_min<0) )
    goesneg=1;

  for(i=0; i<z_size*y_size*x_size*t_size; i++)
    {
      if ( goesneg )
	{
	  if (mask[i%(z_size*y_size*x_size)]<0.5)
	    tmpin[i]=0;
	  else
	    tmpin[i]=(FDT)( ((double)(MIN(MAX(in[i],thresh2),thresh98))-(double)thresh2) * (im.dtmax-1) / ((double)thresh98-(double)thresh2) + 1 );

	  im.min=0;
	  im.max=(FDT)im.dtmax;
	}
      else
	tmpin[i]=mask[i%(z_size*y_size*x_size)]*in[i];
    }

  im.i=tmpin;
  im.t=t_size;
  avw_write(argv[2],im);
  im.t=1;
  free(tmpin);
}

if ( (output_mask) || (output_brain) )
  free(mask);

/* }}} */
  /* {{{ output cost */

if (output_cost)
{
  im.i=raw;
  im.min=raw[0];
  im.max=raw[0];

  for(i=0; i<z_size*y_size*x_size; i++)
    {
      if (raw[i]>im.max) im.max=raw[i];
      if (raw[i]<im.min) im.min=raw[i];
    }

  sprintf(filename,"%s_cost",argv[2]);
  avw_write(filename,im);
  free(raw);
}

/* }}} */
  /* {{{ find skull and output */

if (output_skull)
{
  FDT *skull = (FDT *) malloc(sizeof(FDT)*x_size*y_size*z_size);

  im.i=skull;
  memset(skull,(unsigned char)0,sizeof(FDT)*x_size*y_size*z_size);
  draw_surface(&im,1,1,v,pc); /* tell skull searching where to start */

  /* {{{ create skull image */

im.i=in;

for(z=0; z<z_size; z++)
  for(y=0; y<y_size; y++)
    for(x=0; x<x_size; x++)
      if (IA(skull,x,y,z)==1)
        {
	  FDT val, val2, minval, maxval;
	  double nx, ny, nz, d, d_max=0, d_min, grad, maxgrad, X, Y, Z, maxJ, lastJ;
	  int xx, yy, zz, j=0;

	  /* {{{ zero this point */

IA(skull,x,y,z)=0;

/* }}} */
	  /* {{{ find nearest node and setup normal */

{
  double mind=10000000;

  for(i=0; i<pc; i++)
  { 
    tmpf = (x*im.xv-v[i].x)*(x*im.xv-v[i].x) +
      (y*im.yv-v[i].y)*(y*im.yv-v[i].y) +
      (z*im.zv-v[i].z)*(z*im.zv-v[i].z);
    if (tmpf<mind)
      {
	j=i;
	mind=tmpf;
      }
  }

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

/* }}} */
	  /* {{{ find minval, d_max and maxval up to SKULL_SEARCH distance */

maxval=threshold;
minval=IA(in,x,y,z);

	  for(d=0; d<SKULL_SEARCH; d+=scale*0.5)
	    {
	      xx=x+FTOI(d*nx/im.xv); yy=y+FTOI(d*ny/im.yv); zz=z+FTOI(d*nz/im.zv);

	      if ( (xx>=0) && (xx<x_size) && (yy>=0) && (yy<y_size) && (zz>=0) && (zz<z_size) )
		{
		  val=IA(in,xx,yy,zz);

		  if (val>maxval)
		    {
		      maxval=val;
		      d_max=d;
		    }

		  if (val<minval)
		      minval=val;
		}
	    }

/* }}} */
	  if (maxval>threshold) /* can we see an intensity peak on the far side of the skull? */
	    {
	      /* {{{ find furthest out point that has small intensity */

maxgrad=1;
d_min=SKULL_START;
maxJ =-1000000;
lastJ=-2000000;

	  for(d=SKULL_START; d<d_max; d+=scale*0.5)
	    {
	      xx=x+FTOI(d*nx/im.xv); yy=y+FTOI(d*ny/im.yv); zz=z+FTOI(d*nz/im.zv);

	      if ( (xx>=0) && (xx<x_size) && (yy>=0) && (yy<y_size) && (zz>=0) && (zz<z_size) )
		{
		  tmpf=d/30 - IA(in,xx,yy,zz)/((double)(thresh98-thresh2));

		  /*		  if ( (tmpf>maxJ) && (lastJ+1.0/30>tmpf*0.5) )*/
		  if (tmpf>maxJ)
		  {
		    maxJ=tmpf;
		    d_min=d;
		  }

		  lastJ=tmpf;
		}
	    }

/* }}} */
	      /* {{{ find _first_ max gradient out from d_min */

maxgrad=0;

X=x+d_min*nx/im.xv; Y=y+d_min*ny/im.yv; Z=z+d_min*nz/im.zv;

if ( (X>0) && (X<x_size-1) && (Y>0) && (Y<y_size-1) && (Z>0) && (Z<z_size-1) )
{
  val2 = TLI(im,X,Y,Z);

	  for(d=d_min+scale; d<d_max; d+=scale*0.5)
	    {
	      X=x+d*nx/im.xv; Y=y+d*ny/im.yv; Z=z+d*nz/im.zv;

	      if ( (X>0) && (X<x_size-1) && (Y>0) && (Y<y_size-1) && (Z>0) && (Z<z_size-1) )
		{
		  val = TLI(im,X,Y,Z);
		  /*printf("%d %d %d   %f %f %f   %d %d\n",x,y,z,X,Y,Z,(int)val,(int)val2);*/

		  grad=val-val2;
		  val2=val;		 

		  if (grad>0) /* this so that we don't do anything if we're still in the same voxel */
		    {
		      if (grad > maxgrad)
			{
			  maxgrad=grad;
			  d_min=d;
			}
		      else
			d=d_max;
		    }
		}
	    }
}

/*printf("\n");*/

/* }}} */
	      /* {{{ mark this point as skull */

if (maxgrad>0)
{
  xx=x+FTOI(d_min*nx/im.xv);
  yy=y+FTOI(d_min*ny/im.yv);
  zz=z+FTOI(d_min*nz/im.zv);

  if ( (xx>=0) && (xx<x_size) && (yy>=0) && (yy<y_size) && (zz>=0) && (zz<z_size) )
    {
      if (code_skull)
	IA(skull,xx,yy,zz)=(FDT)(d_min*100.0);
      else
	IA(skull,xx,yy,zz)=100;
    }
}

/* }}} */
	    }
	}

im.i=skull;

/* }}} */
  /* {{{ output */

im.min=0;
im.max=100;
sprintf(filename,"%s_skull",argv[2]);
avw_write(filename,im);

/* }}} */

  free(skull);
}

/* }}} */
  /* {{{ output overlay (corrupts input image) */

if (output_overlay)
{
  im.min=thresh2;
  im.max=hist_max;
  im.i=in;
  draw_surface(&im,im.max,im.max,v,pc);
  sprintf(filename,"%s_overlay",argv[2]);
  avw_write(filename,im);
}

/* }}} */
  /* {{{ output xtopol surface (corrupts tessellation) */

  /* {{{ output xtopol surface (corrupts tessellation) */

if (output_xtopol)
{
  int xtopol_pc=pc;
  double ax=0, ay=0, az=0, tmpf=0;

#ifdef DEBUG_NORMALS
  xtopol_pc *= 2;
#endif

  for(i=0; i<xtopol_pc; i++)
    {
      ax += v[i].x;
      ay += v[i].y;
      az += v[i].z;
    }
  ax/=(double)xtopol_pc; ay/=(double)xtopol_pc; az/=(double)xtopol_pc; 

  for(i=0; i<xtopol_pc; i++)
    tmpf += sqrt( (v[i].x-ax)*(v[i].x-ax) + (v[i].y-ay)*(v[i].y-ay) + (v[i].z-az)*(v[i].z-az) );
  tmpf/=(double)xtopol_pc;

  for(i=0; i<xtopol_pc; i++)
    {
      v[i].x=(v[i].x-ax)/tmpf;;
      v[i].y=(v[i].y-ay)/tmpf;;
      v[i].z=(v[i].z-az)/tmpf;;
    }

  xtopol_output(v,xtopol_pc,argv[2]);
}

/* }}} */
  /* {{{ COMMENT output xtopol surface (x>0) */

#ifdef FoldingComment

if (output_xtopol)
{
  int xtopol_pc=pc, j, k, l, *mapping;
  double ax=0, ay=0, az=0, tmpf=0;

#ifdef DEBUG_NORMALS
  xtopol_pc *= 2;
#endif

  mapping = (int*)malloc(sizeof(int)*xtopol_pc);

  for(i=0; i<xtopol_pc; i++)
    {
      ax += v[i].x;
      ay += v[i].y;
      az += v[i].z;
    }
  ax/=(double)xtopol_pc; ay/=(double)xtopol_pc; az/=(double)xtopol_pc; 

  for(i=0; i<xtopol_pc; i++)
    tmpf += sqrt( (v[i].x-ax)*(v[i].x-ax) + (v[i].y-ay)*(v[i].y-ay) + (v[i].z-az)*(v[i].z-az) );
  tmpf/=(double)xtopol_pc;

  j=0;
  for(i=0; i<xtopol_pc; i++)
    {
      if (v[i].x-ax>0)
	{
	  mapping[i]=j;
	  v[j].x=v[i].x;
	  v[j].y=v[i].y;
	  v[j].z=v[i].z;
	  l=0;
	  for(k=0; v[i].n[k]!=-1; k++)
	    {
	      if (v[v[i].n[k]].x-ax>0)
		{
		  v[j].n[l]=v[i].n[k];
		  l++;
		}
	    }
	  v[j].n[l]=-1;
	  j++;
	}
    }

  for(i=0; i<j; i++)
    {
      v[i].x=(v[i].x-ax)/tmpf;
      v[i].y=(v[i].y-ay)/tmpf;
      v[i].z=(v[i].z-az)/tmpf;
      for(k=0; v[i].n[k]!=-1; k++)
	v[i].n[k] = mapping[v[i].n[k]];
    }

  xtopol_output(v,j,argv[2]);
}

#endif

/* }}} */
  /* {{{ COMMENT output xtopol surface (y>0) */

#ifdef FoldingComment

if (output_xtopol)
{
  int xtopol_pc=pc, j, k, l, *mapping;
  double ax=0, ay=0, az=0, tmpf=0;

#ifdef DEBUG_NORMALS
  xtopol_pc *= 2;
#endif

  mapping = (int*)malloc(sizeof(int)*xtopol_pc);

  for(i=0; i<xtopol_pc; i++)
    {
      ax += v[i].x;
      ay += v[i].y;
      az += v[i].z;
    }
  ax/=(double)xtopol_pc; ay/=(double)xtopol_pc; az/=(double)xtopol_pc; 

  for(i=0; i<xtopol_pc; i++)
    tmpf += sqrt( (v[i].x-ax)*(v[i].x-ax) + (v[i].y-ay)*(v[i].y-ay) + (v[i].z-az)*(v[i].z-az) );
  tmpf/=(double)xtopol_pc;

  j=0;
  for(i=0; i<xtopol_pc; i++)
    {
      if (v[i].y-ay>0)
	{
	  mapping[i]=j;
	  v[j].x=v[i].x;
	  v[j].y=v[i].y;
	  v[j].z=v[i].z;
	  l=0;
	  for(k=0; v[i].n[k]!=-1; k++)
	    {
	      if (v[v[i].n[k]].y-ay>0)
		{
		  v[j].n[l]=v[i].n[k];
		  l++;
		}
	    }
	  v[j].n[l]=-1;
	  j++;
	}
    }

  for(i=0; i<j; i++)
    {
      v[i].x=(v[i].x-ax)/tmpf;
      v[i].y=(v[i].y-ay)/tmpf;
      v[i].z=(v[i].z-az)/tmpf;
      for(k=0; v[i].n[k]!=-1; k++)
	v[i].n[k] = mapping[v[i].n[k]];
    }

  xtopol_output(v,j,argv[2]);
}

#endif

/* }}} */
  /* {{{ COMMENT output xtopol surface (z>0) */

#ifdef FoldingComment

if (output_xtopol)
{
  int xtopol_pc=pc, j, k, l, *mapping;
  double ax=0, ay=0, az=0, tmpf=0;

#ifdef DEBUG_NORMALS
  xtopol_pc *= 2;
#endif

  mapping = (int*)malloc(sizeof(int)*xtopol_pc);

  for(i=0; i<xtopol_pc; i++)
    {
      ax += v[i].x;
      ay += v[i].y;
      az += v[i].z;
    }
  ax/=(double)xtopol_pc; ay/=(double)xtopol_pc; az/=(double)xtopol_pc; 

  for(i=0; i<xtopol_pc; i++)
    tmpf += sqrt( (v[i].x-ax)*(v[i].x-ax) + (v[i].y-ay)*(v[i].y-ay) + (v[i].z-az)*(v[i].z-az) );
  tmpf/=(double)xtopol_pc;

  j=0;
  for(i=0; i<xtopol_pc; i++)
    {
      if (v[i].z-az>0)
	{
	  mapping[i]=j;
	  v[j].x=v[i].x;
	  v[j].y=v[i].y;
	  v[j].z=v[i].z;
	  l=0;
	  for(k=0; v[i].n[k]!=-1; k++)
	    {
	      if (v[v[i].n[k]].z-az>0)
		{
		  v[j].n[l]=v[i].n[k];
		  l++;
		}
	    }
	  v[j].n[l]=-1;
	  j++;
	}
    }

  for(i=0; i<j; i++)
    {
      v[i].x=(v[i].x-ax)/tmpf;
      v[i].y=(v[i].y-ay)/tmpf;
      v[i].z=(v[i].z-az)/tmpf;
      for(k=0; v[i].n[k]!=-1; k++)
	v[i].n[k] = mapping[v[i].n[k]];
    }

  xtopol_output(v,j,argv[2]);
}

#endif

/* }}} */

/* }}} */

  if (verbose) printf("\n");
  return(0);
}

/* }}} */