cbp.c

关于生物光学程序的代码

#include <math.h> 
#include <stdio.h> 
#include <string.h> 
#include "io_mat.h"
#include "A_comp.h"
#include "fft.h"
#include "allocate.h"


#ifdef __linux__
#include <getopt.h>
#endif
extern int optind;
extern char *optarg;



#define PI_1        3.14159265358979323846
/*#define WRITE_FILTER*/



int VERBOSE=1; 


struct params{
  char mode; /* 't' = transmission, 'e' = emission, */
             /*'E' = emission with external attenuation file */
  char *infile;
  char *outfile;
  char attnfile[256]; /* 'E'-mode: filename of forward projected */
                      /*transm. sinogram for attenuation correction */

  int use_ramp_filter;

  double freq_cutoff; /* filter cutoff between 0 and 1 */
  
  int resol; 
  /* for -r flag, resol=-1 is default, means geom.n_x=geom.n_y=sino.n_t
     if -r flag is given, the resol contains resolution of recon in
     each direction (only square pix supported by Acomp) */

};


/* Routines for filtering sinogram */
void filter_sino(float **Y,struct sino_info *sino,struct params *cmdline);
double* GHamming_Filter(int n,struct params *cmdline,struct sino_info *sino);
double* ramp_Filter(int n,struct params *cmdline,struct sino_info *sino);
void conv(double *out, double *in1, int len1, double *in2, int len2);
void normalize_filter_to_sum1(double *filt, int len);
double hrl(int m, struct sino_info *sino);
double dhrl(double m, struct sino_info *sino);

/* routines for backproj */
void backproject(float **image,float **Y,struct sino_info *sino,struct geom_info *geom);
double interpolate(double findex, float **Y, int angle_index);

/* Special routines for transmission case */
void transmission_transformY(float **Y,float **dosage,struct sino_info *sino);

/* Special routines for emission case */
void attn_transformY(float **Y,float **attn,struct sino_info *sino);

/* Mean-value-cor*/
void mean_value_correct(float **image, double pmass, struct geom_info *geom);
double calc_pmass(float **Y, struct sino_info *sino, struct geom_info *geom);

/* Input/Output */
void read_mat_file(struct params *cmdline,float ***counts,float ***dosage,float ***attn,struct sino_info *sino);
int write_image(char *fname,char *varname,float **image,int rows, int cols);
void parse_cmdline(struct params *cmdline,int argc, char **argv);
void cmdline_error();
double snc(double n);
void cut_freq(double *filter,double del,int n);

/* General */
float min(float **img,int N, int M);
float max(float **img,int N, int M);


int main(int argc, char **argv)
{
  float **Y;           /* sinogram */
  float **image;       /* output image */
  float **dosage;      /* only for transmission case */
  float **attn;        /* for 'E' case,  attenuation sinogram, e.g. forw.proj. trn recon*/
  struct geom_info geom; /* contains reconstruction geometry */
  struct sino_info sino; /* contains sinogram info, def. in Acomp.h */ 
  struct params cmdline; /* contains command-line options, filenames */
  double pmass;
  int i,j;
 

  /* parse command-line */
  parse_cmdline(&cmdline,argc, argv);
  
 
  /* read input */
  read_mat_file(&cmdline,&Y,&dosage,&attn,&sino);
  printf("t_0=%f delta_theta=%f  delta_t=%f\n",sino.t_0,sino.delta_theta,sino.delta_t);

  /* If in transmission mode, take Y=log(dosage/Y) */
  if (cmdline.mode=='t')
    transmission_transformY(Y,dosage,&sino);

  /* If in transmission mode, take Y=Y./attn */
  if (cmdline.mode=='E')
    attn_transformY(Y,attn,&sino);

  
  /* set size and parameters of reconstruction */
  /* image geometry */
  
  /* if -r option is not given, then resol==-1, in that case
     set geom.n_y=geom.n_x=sino.n_t;
  */
  if (cmdline.resol == -1)
    {
      geom.n_y=geom.n_x=sino.n_t;
      geom.delta=sino.delta_t;
      geom.x_0=-0.5*geom.delta*(double)(geom.n_x-1);
      geom.y_0=-0.5*geom.delta*(double)(geom.n_y-1);
    }
  else
    {
      geom.n_y=geom.n_x=cmdline.resol;
      geom.x_0 = geom.y_0 = sino.t_0;
      geom.delta = -2.0*sino.t_0/((double)(geom.n_y-1));
      /* note: assume that reconstruction starts at outer limits
	 of parallel projections. 
      */
    }
  


  /* Calculate pmass of input sinogram for mean */
  /* value normalization later                  */
  pmass=calc_pmass(Y, &sino, &geom);
  
  /* allocate and init output image */
  image = (float **)get_img(geom.n_x, geom.n_y, sizeof(float));
  for(i=0; i<geom.n_y; i++)
    for(j=0; j<geom.n_x; j++) 
      image[i][j] = 0.0;
  
  /* Filter Sinogram */
  filter_sino(Y,&sino,&cmdline);

  /* Backprojection */
  backproject(image,Y,&sino,&geom);

  /* mean-value correction */
  mean_value_correct(image, pmass, &geom);

  /* Write output image as .mat file */
  write_image(cmdline.outfile,"img",image,geom.n_y, geom.n_x);
  
  /* Free output-image and sinogram */
  free_img((void **)image);
  free_img((void **)Y);
  if (cmdline.mode=='t')
    free_img((void **)dosage);

  return(0);
   
}



void filter_sino(float **Y,struct sino_info *sino,struct params *cmdline)
{
  double *data,*filter,*ans;
  int n,i,j;

#ifdef WRITE_FILTER
  FILE *fp;
#endif

  printf("Filtering...\n");

  /* Calculate filter-length */
  n = (int) ceil(log(sino->n_t)/log(2.0));
  n = pow(2,n);

  /* Calculate Filter  of length (2*n) */
  if (cmdline->use_ramp_filter)
    filter=ramp_Filter(n,cmdline,sino);
  else
    filter=GHamming_Filter(n,cmdline,sino);

  /* Allocate space for data */
  ans = (double *)get_spc(4*n, sizeof(double));
  data = (double *)get_spc(2*n, sizeof(double));

  /* Pad data with zeros */
  for(i=sino->n_t;i<2*n;i++)
    data[i] = 0.0;
 

  /* Filter sinogram Y for each projection angle */
  for(i=0;i<sino->n_theta;i++) 
    {
      for(j=0;j<sino->n_t;j++)
	data[j] = (double)Y[i][j];
  
      convlv(data-1,2*n,filter-1,ans-1);
      
      /* insert into Y */
      for(j=0;j<sino->n_t;j++)
	Y[i][j] = (float)ans[j]; 
      
    }



#ifdef WRITE_FILTER
  /* Write filter response to disk for checking */
  fp=fopen("filter","wa");
  for (i=0;i<2*n;i++)
    fprintf(fp,"%g ",filter[i]);
  fclose(fp);
#endif
  
  free((void *)(ans));
  free((void *)(data));
  free((void *)(filter));
} 




/* This reconstruction filter GHamming_Filter has frequency response

      H(omega) = H_ramp(omega) * (0.5 + 0.5cos(pi*omega/omega_c))

   where H_ramp is the ideal ramp reconstruction filter, 
   and omega_c is the cutoff frequency. 
   The filter input parameter specifies the ratio omega/omega_c .
*/
double* GHamming_Filter(int n,struct params *cmdline,struct sino_info *sino)
{
  double *filter;
  int i;
  double del;


  /* allocate filter */
  filter = (double *)get_spc(2*n, sizeof(double));
  
  del = cmdline->freq_cutoff;
  printf("cos delay=%f, so cutoff is %f percent.\n",del,100.0*del); 
  
  /* Calculate impulse response in time domain, 
     hrl is ideal ramp filter response 1/(4*delta_t)*(2*delta(n)-sinc^2(n/2)),
     dhrl is basically same as hrl but w/ double arithm. input, 
          and uses 1/(4*delta_t)*(2*sinc(n)-sinc^2(n/2))
  */
  for(i=-n+1;i<=n;i++)
    filter[i+n-1] = 0.5*(hrl(i,sino) + 0.5*(dhrl(i-1.0/del,sino) + dhrl(i+1.0/del,sino) )); 
  
  /* So far, filter response is only accurate for omega<=omega_c */
  /* convolve with sinc() to chop off frequency domain portion
     that is over cutoff frequency */
  cut_freq(filter,del,n);
  

  /* Resort filter : before, filter had phase delay,*/
  /* i.e. ideal h[0] was at filter[len-1], now for fft-filtering */
  /* resort to fft-symmetry to get lin. phase */
 
  for (i=0;i<=n;i++)
    filter[i]=filter[i+n-1];
  
  for (i=n+1;i<2*n;i++)
    filter[i]=filter[2*n-i];
  

  return(filter);
}


/* This routine implements the ideal ramp filter but with 
   a cutoff frequency as specified on the commandline. 
   So
   H(omega) = H_ramp(omega) * rect(omega/omega_c) 

   where rect(x)=1 for |x|<1
   and H_ramp is the ideal highpass.
*/
double* ramp_Filter(int n,struct params *cmdline,struct sino_info *sino)
{
  double *filter;
  int i;
  double del;

  /* allocate filter */
  filter = (double *)get_spc(2*n, sizeof(double));
  
  del = cmdline->freq_cutoff;
  printf("cos delay=%f, so cutoff is %f percent.\n",del,100.0*del); 
  
  /* Calculate impulse response in time domain, 
     hrl is ideal ramp filter response 1/(4*delta_t)*(2*delta(n)-sinc^2(n/2)),
   */
  for(i=-n+1;i<=n;i++)
    filter[i+n-1] = hrl(i,sino); 
  
  /* So far, filter response is only accurate for omega<=omega_c */
  /* convolve with sinc() to chop off frequency domain portion
     that is over cutoff frequency */
  cut_freq(filter,del,n);
  

  /* Resort filter : before, filter had phase delay,*/
  /* i.e. ideal h[0] was at filter[len-1], now for fft-filtering */
  /* resort to fft-symmetry to get lin. phase */
 
  for (i=0;i<=n;i++)
    filter[i]=filter[i+n-1];
  
  for (i=n+1;i<2*n;i++)
    filter[i]=filter[2*n-i];
  
  return(filter);

}




/* convolution , used in compute_filter to */
/* convolve Gaussian*HammW with ideal recon-filter*/

void conv(double *out, double *in1, int len1, double *in2, int len2)
{
  int i,n;

  for (n=0;n<len1+len2-1;n++)
    {
      out[n]=0.0;
      for (i=0;i<len1;i++)
	if ((n-i>=0) && (n-i<len2))
	  out[n] += (in1[i]*in2[n-i]);
    }
}



void normalize_filter_to_sum1(double *filt, int len)
{
  int i;
  double sum=0.0;
 
  for (i=0;i<len;i++)
    sum+=(filt[i]);
 
  for (i=0;i<len;i++)
    filt[i]/=sum;
}




/* hrl returns 1/(4*delta_t)*(   2*delta(m)- (sinc(m/2))^2)  */
/* this is the ideal reconstruction filter */

double hrl(int m, struct sino_info *sino)
{ 
    float temp;

    if (m == 0)
        return(1/(4.0*sino->delta_t));
    else
    {
        temp = PI_1*m/2;
        temp = sin(temp)/temp;
        temp =  - temp*temp;
        return( temp/(4.0*sino->delta_t) );
    }
}

double dhrl(double m, struct sino_info *sino)
{ 
    double temp,temp2,omeg;

    if (fabs(m)<1e-8)
        return(1.0/(4.0*sino->delta_t));
    else
    {
        omeg = PI_1*m/2.0;
        temp = sin(omeg)/omeg;
        temp2 = 2*sin(omeg*2)/(omeg*2);
        temp = temp2 -temp*temp;
        return( temp/(4.0*sino->delta_t) );
    }
}







void backproject(float **image,float **Y,struct sino_info *sino,struct geom_info *geom)
{
  double x=0.0,y=0.0;
  double theta,t,dat,findex;
  int i,j,k;
  double *cs,*sn; 
       /* arrays to precompute cos(theta_k), sin(theta_k) */
  
  float *im_pt; 
       /* pointer to image[i][j]  Assume image is */
       /* allocated as 1-D array w/ pointers on each row */

  printf("Backprojecting ... \n");


  /* Precompute sin and cos for all projection angles */
 
  cs=(double*)malloc(sizeof(double)*sino->n_theta);
  sn=(double*)malloc(sizeof(double)*sino->n_theta);
  theta=sino->theta_0;
  for (k=0;k<sino->n_theta;k++)
    {
     sn[k]=sin(theta);
     cs[k]=cos(theta);
     theta+=sino->delta_theta;
    }



  /* Start Backprojection */
  im_pt=&image[0][0];
  y=geom->y_0;
  for (i=0;i<geom->n_y;i++)
    {
      x=geom->x_0;
      for (j=0;j<geom->n_x;j++)
	{
	  for (k=0;k<sino->n_theta;k++)
	    {
	      t=x*cs[k]-y*sn[k];
	      findex=(t-sino->t_0)/sino->delta_t;
	      
	      /* Check whether location is out of */
              /* projection area for this angle */
	      /* If so,  set image[i][j] to zero, */
              /* break k loop and go to next image location */
	      if ((findex<0) || (findex > sino->n_t - 1 ))
		{
		  *im_pt=0.0;
		  break;
		}
	     
	      dat=interpolate(findex, Y, k); 
                  /* bi-linear interpolation between par. proj */