cbp.c

关于生物光学程序的代码

	      
	      *im_pt += dat;
	    }
	  
	  /* Clip to >= zero, normalize */
	  if (*im_pt <= 0.0)
	    *im_pt = 0.0;
	  else
	    *im_pt *= (PI_1/(float)sino->n_theta);

	  im_pt++;
	  x += geom->delta;
	}
      y += geom->delta;
    }
  
  free((void*)cs);
  free((void*)sn);
}



/* Interpolate is used in backproject to */
/* perform bi-linear interploation between parallel projections */

double interpolate(double findex, float **Y, int angle_index)
{
  int index,ce,fl;
  double out;
  
  /* check whether findex is integer */
  index=(int)findex;
  if ((findex-(double)index>-1E-10) && (findex-(double)index<1E-10))
    out=Y[angle_index][index];
  else /* if not integer, then bilinear interpolation between parallel projections */
    {
      ce=(int)(findex+1.0); /* ceiling */
      fl=(int)findex; /* floor  */

      out=(((double)ce)-findex)*Y[angle_index][fl] + (findex-(double)fl)*Y[angle_index][ce];
    }
 
  return(out);

}


void transmission_transformY(float **Y,float **dosage,struct sino_info *sino)
{
  int i,j;
  float temp;

  for(i=0;i<sino->n_theta;i++)
    for(j=0;j<sino->n_t;j++)
      {
	temp = Y[i][j];
	if (temp <= 0.0) /* I'm not sure whether this */
                         /* is good..., CHECK         */
	  temp = 1.0;
	if (dosage[i][j]<=0.0)
	  dosage[i][j]=1.0;
	Y[i][j] = log( dosage[i][j]/temp );
	if (Y[i][j] < 0.0)
	  Y[i][j] = 0.0;
      }
}


void attn_transformY(float **Y,float **attn,struct sino_info *sino)
{
  int i,j;
  float temp;

  printf("performing attn correction.\n");

  for(i=0;i<sino->n_theta;i++)
    for(j=0;j<sino->n_t;j++)
      {
	temp = attn[i][j];
	if (temp <= 0.0) /* I'm not sure whether this */
                         /* is good..., CHECK         */
	  {
	    temp = 1.0;
	    fprintf(stderr,"Warning: attn sinogram has zero. Set to 1.0.\n");
	  }

	Y[i][j] /= temp;

      }
}



void mean_value_correct(float **image, double pmass, struct geom_info *geom)
{
  double imass=0.0;
  float  cor_factor,mx;
  int    i,j;
  
  for (i=0;i<geom->n_y;i++)
    for (j=0;j<geom->n_x;j++)
      imass += (double)image[i][j];
  imass /= ((double)(geom->n_y*geom->n_x)); 
  
  cor_factor=(float)(pmass/imass);
  
  printf("pmass:%f imass:%f cor_factor:%f\n",pmass,imass,cor_factor);

  for (i=0;i<geom->n_y;i++)
    for (j=0;j<geom->n_x;j++)
      image[i][j] *= cor_factor;

  /* Calculate Mean just to printf */
  mx = 0.0;
  for (i=0;i<geom->n_y;i++)
    for (j=0;j<geom->n_x;j++)
      mx += image[i][j];

  mx /= (float)(geom->n_y*geom->n_x);

 
  printf("The mean of the CBP image is %f\n",mx);
  printf("The min,max of the CBP image are %.3f,%.3f\n",min(image,geom->n_y,geom->n_x),max(image,geom->n_y,geom->n_x));
} 




/* Calculates sinogram mean and re-normalizes it     */
/* for the reconstruction geometry in geom           */
/* The output is later used in mean_value_correction */

double calc_pmass(float **Y, struct sino_info *sino, struct geom_info *geom)
{
  double pmass;
  double meanY=0.0;
  int i,j;
  
  for (i=0;i<sino->n_theta;i++)
    for (j=0;j<sino->n_t;j++)
      meanY+=(double)Y[i][j];

  meanY=meanY/((double)(sino->n_theta*sino->n_t));

  pmass = meanY*((double)sino->n_t)*sino->delta_t/(geom->delta*geom->delta*(double)(geom->n_y*geom->n_x)); 

  return(pmass);
}




/* routine to read sinogram data */
void read_mat_file(struct params *cmdline,float ***counts,float ***dosage,float ***attn,struct sino_info *sino)
{
  MATFile *fin;
  float **tmp_counts,**tmp_dosage=NULL,**tmp_attn=NULL;
  int ntheta,nr;


  fin = matOpen(cmdline->infile, "r");
  if (fin == NULL) 
    {
      fprintf(stderr,"Cannot read MatLab file %s \n",cmdline->infile);
      exit(-1);
    }

  /* Routine  read_mat_float_image is in ansi_lib/io_mat.c */
  tmp_counts = read_mat_float_image(&ntheta,&nr,"counts",fin);
  if (tmp_counts == NULL) 
    {
      fprintf(stderr,"Cannot read counts matrix \n");
      exit(-1);
    }


  /* If in transmission-mode, read dosage matrix */
  if (cmdline->mode=='t')
    {
      tmp_dosage = read_mat_float_image(&ntheta,&nr,"dosage",fin);
      if (tmp_dosage == NULL) 
	{
	  fprintf(stderr,"Cannot read dosage \n");
	  exit(-1);
	}
    }


  /* Read n_theta, n_t, delta_theta, delta_t, */
  /* theta_0, t_0,  dec_a=1, dec_t=1          */
  /* Routine in ansi_lib/Acomp.c              */
  read_mat_sino_info(fin, sino);
  matClose(fin);

  
  /* Check if dimensions of counts matrix agree */
  /* with sino->n_theta, sin->n_t from variables in .mat file*/
  if (ntheta!=sino->n_theta)
    {
      fprintf(stderr,"Error in input file: Matrix has %d columns, but n_theta=%d.\n",ntheta,sino->n_theta);
      exit(-1);
    }
  if (nr!=sino->n_t)
    {
      fprintf(stderr,"Error in input file: Matrix has %d rows, but n_t=%d.\n",nr,sino->n_t);
      exit(-1);
    } 

  /* If in 'E' emission mode, load attn-file */
  if (cmdline->mode == 'E')
    {
      /* open file */
      fin = matOpen(cmdline->attnfile, "r");
      if (fin == NULL) 
	{
	  fprintf(stderr,"Cannot read MatLab file %s \n",cmdline->attnfile);
	  exit(-1);
	}
      
      /* read attn-sino*/
      tmp_attn = read_mat_float_image(&ntheta,&nr,"trn",fin);
      if (tmp_attn == NULL) 
	{
	  fprintf(stderr,"Cannot read trn matrix of file: %s \n",cmdline->attnfile);
	  exit(-1);
	}
      matClose(fin);

      /* check size */
      if (ntheta!=sino->n_theta)
	{
	  fprintf(stderr,"Error in input file: Matrix has %d columns, but n_theta=%d.\n",ntheta,sino->n_theta);
	  exit(-1);
	}
      if (nr!=sino->n_t)
	{
	  fprintf(stderr,"Error in input file: Matrix has %d rows, but n_t=%d.\n",nr,sino->n_t);
	  exit(-1);
	} 
    }


  /* Set output pointers */
  *counts=tmp_counts;
  if (cmdline->mode=='t')
    *dosage=tmp_dosage;
  if (cmdline->mode=='E')
    *attn=tmp_attn;

}




/* Writes float image as a .mat file */
int write_image(char *fname,char *varname,float **image,int rows, int cols)
{
  MATFile *fout;
  int status;

  fout = matOpen(fname, "w");
  if (fout == NULL) 
    {
      fprintf(stderr,"Cannot write to MatLab file %s\n",fname);
      exit(-1);
    }
  
  status=write_mat_float_image_cor(image,rows,cols,varname,fout);
  matClose(fout);
  return(status);
}  



void parse_cmdline(struct params *cmdline,int argc, char **argv)
{
  int i,ct;
  int trans_FLAG=0;
  int attn_FLAG=0;

  /* Set defaults */
  cmdline->mode = 'e';
  cmdline->use_ramp_filter = 0;
  cmdline->freq_cutoff = 0.8;
  cmdline->resol=-1;

  /* read cmdline options */
  for (;;)
    { 
      i= getopt(argc, argv, "r:f:tE:R");
      if (i==EOF)
        break;
      else
        {
          if (i == '?')
            cmdline_error();
	  if (i=='r') 
	    cmdline->resol=atoi(optarg);
	  if (i=='f')
	    cmdline->freq_cutoff = atof(optarg);
	  if (i=='t')
	    {
	      cmdline->mode = 't';
	      trans_FLAG=1;
	    }
	  if (i=='E')
	    {
	      cmdline->mode = 'E';
	      attn_FLAG=1;
	      strcpy(cmdline->attnfile,optarg);
	    }
	  if (i=='R')
	    cmdline->use_ramp_filter=1;
   	  
	}
    }
  
  /* check cmdline options */
  ct=optind;
  if ((argc-optind < 2) || (argc-optind > 2))
    {
      fprintf(stderr,"\nError: Incorrect number of arguments.\n\a\n");
      cmdline_error();
    }

  if ((attn_FLAG) && (trans_FLAG))
    {
      fprintf(stderr,"\nError: Must not specify both -E and -t options.\n\a\n");
      cmdline_error();
    }



  cmdline->infile=argv[ct];
  ct++;
  cmdline->outfile=argv[ct];
  ct++;


  if (cmdline->freq_cutoff > 1.0)
    {
     fprintf(stderr,"Error: frequency cutoff=%f, must not be larger than 1. In this implementation, a cutoff >1 does NOT transition to a ramp filter.\n\a\n",cmdline->freq_cutoff);
     exit(-1);
    }



  /* Print options */
  printf("\n");
  if (cmdline->mode == 'e')
    printf("mode = emission, no attenuation file\n");
  
  if (cmdline->mode == 'E')
    printf("mode = emission, attenuation file = %s\n",cmdline->attnfile);
  
  if (cmdline->mode == 't')
    printf("mode = transmission\n");
  
  if (cmdline->use_ramp_filter)
    printf("filter = ramp\n");
  else
    printf("filter = generalized Hamming\n");

  if (cmdline->resol == -1)
    printf("reconstruction resolution = sino.n_t (default)\n");
  else
    printf("reconstruction resolution = %d\n",cmdline->resol);

  
  printf("filter cutoff = %f\n",cmdline->freq_cutoff);
  printf("input file = %s\n",cmdline->infile);
  printf("output file = %s\n",cmdline->outfile);
  printf("\n");

}



void cmdline_error()
{
  fprintf(stderr,"\nUsage: cbp [-r <resolution_in_pixels>] [-R] [-f <cutoff_freq>] [-t] [-E <attenuation-file>] <input_file.mat> <output-file.mat>\n\n");
  fprintf(stderr,"-r: number of pixels in recon in each dir. (default=sino.n_t) \n");
  fprintf(stderr,"-R: Use ramp filter instead of default Gen. Hamming filter \n");
  fprintf(stderr,"-f: Filter cutoff frequency between 0.0 and 1.0, default=0.8\n");
  fprintf(stderr,"-t: Transmission mode (default=emission w/o attenuation file)\n");
  fprintf(stderr,"-E: Allows to specify attenuation file for emission mode\n");
  exit(-1);
}




float min(float **img,int N, int M)
{
  int i,j;
  float min;
  
  min = img[0][0];
  for (i=0;i<N;i++)
    for(j=0;j<M;j++)
      if ( img[i][j] < min )
	min = img[i][j];
  
  return(min);
}


float max(float **img,int N, int M)
{
  int i,j;
  float max;
  
  max = img[0][0];
  
  for (i=0;i<N;i++)
    for(j=0;j<M;j++)
      if (img[i][j] > max)
	max = img[i][j];
  
  return(max);
}




void cut_freq(double *filter,double del,int n)
{
  double flt2[1024];
  int i,k;
  
  for(i=-n+1;i<=n;i++)
    {
      flt2[i+n-1]=0.0;
      for (k=-n+1;k<=n;k++)
	flt2[i+n-1] += (filter[k+n-1] * snc(del*(double)(i-k))*del);
    } 
  for(i=-n+1;i<=n;i++)
    filter[i+n-1] = flt2[i+n-1];
  
}

double snc(double n)
{
  double out,x;

  if (fabs(n)<1e-9)
    out=1.0;
  else
    {
      x=PI_1*n;
      out=sin(x) / x;
    }
  return(out);


}