prmat.c

Go to the documentation of this file.

/*****************************************************************************/
int PRMAT_TEST;
int PRMAT_VERBOSE;
/*****************************************************************************/
#include "libtpcrec.h"
/*****************************************************************************/
 
/*****************************************************************************/
/* Initialization and memory handling for prmat data      */
/*****************************************************************************/
void prmatInit(PRMAT *mat)
{
  if(PRMAT_VERBOSE) printf("PRMAT:prmatInit() \n"); 
  
  /*init the information of projection matrix*/
  //buffer mat to contain zero
  memset(mat, 0, sizeof(PRMAT)); 
  mat->type=PRMAT_TYPE_NA;
  mat->viewNr=0;
  mat->binNr=0;
  mat->mode=PRMAT_DMODE_01;
  mat->prmatfov= (int*)NULL;
  mat->max=mat->min=0.0;
  mat->scaling_factor=1.0;
  mat->factor_sum=0.0;
  mat->factor_sqr_sum=(float*)NULL;
 
  /*init look-up table*/
  mat->nrp=0;
  mat->lines=mat->_linesdata=(unsigned int**)NULL; //data array
 
  /*init projection data*/
  mat->dimr=0;
  mat->dime=(unsigned int*)NULL;
  mat->fact=mat->_factdata=(unsigned short int***)NULL; //data array
  
  mat->status=PRMAT_STATUS_INITIALIZED;
}
/*****************************************************************************/
void prmatEmpty(PRMAT *mat)
{
  unsigned int ir, ip;
  if(PRMAT_VERBOSE) printf("PRMAT:prmatEmpty() \n");
  
  //if mat is allready empty
  if(mat->status==PRMAT_STATUS_INITIALIZED) return; 
  
  //if look-up table is occupied
  if(mat->status==PRMAT_STATUS_LU_OCCUPIED) { 
    /*free the occupied memory*/
    for(ir=0; ir<mat->nrp; ir++) { //every row
      free((unsigned int*)mat->_linesdata[ir]);
    }
    if(mat->nrp > 0) free(mat->_linesdata);
  }
 
  //if projection matrix is occupied
  if(mat->status==PRMAT_STATUS_PR_OCCUPIED || 
     mat->status==PRMAT_STATUS_BS_OCCUPIED )
  { 
    free((float*)mat->factor_sqr_sum);
    for(ir=0; ir<mat->dimr; ir++){ //every row
      for(ip=0; ip<mat->dime[ir]; ip++){ //every factor
  free((unsigned short int*)mat->_factdata[ir][ip]);
      }
    }
 
    free((int*)mat->dime);
  
    if(mat->dimr>0) free(mat->_factdata); 
    mat->dimr=0; 
  }
  return;
}
/*****************************************************************************/
int prmatAllocate(PRMAT *mat, int set, unsigned int rows, unsigned int *coords)
{
  unsigned int *iptr, ir, ip;
  if(PRMAT_VERBOSE) printf("PRMAT:prmatAllocate(*mat,%d,*coords) \n",rows);
  
  /*check the arguments*/
  if(mat->status==PRMAT_STATUS_UNINITIALIZED) return(1);
  if(rows<1) return(2);
  
  //If set != 0 allocate for projection matrix.
  if(set){
    /*allocate memory*/
    mat->factor_sqr_sum=calloc(rows,sizeof(float));
    mat->dimr=rows; //for dimr
    mat->dime=calloc(rows,sizeof(int)); 
    mat->_factdata=
     (unsigned short int***)calloc(rows,sizeof(unsigned short int**)); 
    //if not enough memory
    if(mat->_factdata==NULL || mat->dime==NULL) return(4);
    //iterate through the rows
    for(ir=0,iptr=coords; ir<rows; ir++){
      //copy the coords vector in the same loop
      mat->dime[ir]=*iptr; 
      //*iptr is the number of (non-zero) factors in this row
      mat->_factdata[ir]=
       (unsigned short int**)calloc(*iptr,sizeof(unsigned short int*)); 
      //if not enough memory
      if(mat->_factdata[ir]==NULL) return(4); 
      //iterate through the factors in this row
      for(ip=0; ip<(*iptr); ip++){ 
  //NOTE that the leaves are coordinates and factor in this order
  mat->_factdata[ir][ip]=
          (unsigned short int*)calloc(3,sizeof(unsigned short int)); 
  if(mat->_factdata[ir][ip]==NULL) return(4);
      }
      iptr++;
    }
    /*set data pointers; these pointers are used to avoid data loss*/
    mat->fact=mat->_factdata;    
  } else {
    //If set was 0 allocate for look-up table.
    //allocate memory for look up table
    mat->nrp=rows; 
    mat->_linesdata=(unsigned int**)calloc(rows,sizeof(unsigned int*)); 
    //if not enough memory
    if(mat->_linesdata==NULL) return(4);
    //iterate through the rows
    for(ir=0,iptr=coords; ir<rows; ir++){
      //iptr is the number of lors belonging to this pixel
      mat->_linesdata[ir]=(unsigned int*)calloc(*iptr + 2,sizeof(unsigned int)); 
      //if not enough memory
      if(mat->_linesdata[ir]==NULL) return(4); 
      iptr++;
    }
    //set data pointers; these pointers are used to avoid data loss
    mat->lines=mat->_linesdata;    
  }
  return(0);
}
/*****************************************************************************/
/* GET procedures for PRMAT datatype. USE ONLY THESE FOR GETTING. */
/*****************************************************************************/
/* GET FUNCTIONS FOR THE PARAMETERS */
/*****************************************************************************/
unsigned int prmatGetNV(PRMAT *mat)
{
  return mat->viewNr;
}
/*****************************************************************************/
unsigned int prmatGetNB(PRMAT *mat)
{
  return mat->binNr;
}
/*****************************************************************************/
unsigned int prmatGetID(PRMAT *mat)
{
  return mat->imgDim;
}
/*****************************************************************************/
/* GET FUNCTIONS FOR ACCESSING LOOK-UP TABLE*/
/*****************************************************************************/
unsigned int prmatGetPIX(PRMAT *mat)
{
  return mat->nrp;
}
/*****************************************************************************/
unsigned int prmatGetPixCoord(PRMAT *mat,int row)
{
  return mat->lines[row][0];
}
/*****************************************************************************/
unsigned int prmatGetRays(PRMAT *mat,int row)
{
  return mat->lines[row][1];
}
/*****************************************************************************/
unsigned int prmatGetBinView(PRMAT *mat, int row, int ind)
{
  return mat->lines[row][ind+2];
}
/*****************************************************************************/
/* GET FUNCTIONS FOR ACCESSING PROJECTION MATRIX BY A LINE OF RESPONSE.*/
/*****************************************************************************/
unsigned int prmatGetRows(PRMAT *mat)
{
  //NOTE no checking on parameters to accelerate iteration
  return(mat->dimr);
}
/*****************************************************************************/
unsigned int prmatGetPixels(PRMAT *mat, int row)
{
  return(mat->dime[row]);
}
/*****************************************************************************/
unsigned int prmatGetCoord(PRMAT *mat, int row, int pix)
{
  return (int)(mat->fact[row][pix][1])*mat->imgDim + (int)(mat->fact[row][pix][0]);
}
/*****************************************************************************/
unsigned int prmatGetXCoord(PRMAT *mat, int row, int pix)
{
  return (int)(mat->fact[row][pix][0]);
}
/*****************************************************************************/
unsigned int prmatGetYCoord(PRMAT *mat, int row, int pix)
{
  return (int)(mat->fact[row][pix][1]);
}
/*****************************************************************************/
float prmatGetFactor(PRMAT *mat, int row, int pix)
{
  return((float)(mat->scaling_factor)*(mat->fact[row][pix][2]));
}
/*****************************************************************************/
float prmatGetMajor(PRMAT *mat)
{
  return(mat->prmatfov[0]);
}
/*****************************************************************************/
float prmatGetMinor(PRMAT *mat)
{
  return(mat->prmatfov[1]);
}
/*****************************************************************************/
float prmatGetMin(PRMAT *mat)
{
  return(mat->min);
}
/*****************************************************************************/
float prmatGetMax(PRMAT *mat) {
  return(mat->max);
}
/*****************************************************************************/
float prmatGetFactorSum(PRMAT *mat)
{
  return(mat->factor_sum);
}
/*****************************************************************************/
float prmatGetFactorSqrSum(PRMAT *mat, int row)
{
  return(mat->factor_sqr_sum[row]);
}
/*****************************************************************************/
/* ARCHIVING METHODS */
/*****************************************************************************/
int prmatSaveMatrix(PRMAT *mat)
{
  char matrixfile[FILENAME_MAX], semi_x_c[2], semi_y_c[2], mod[2];
  int semi_x=mat->prmatfov[0], semi_y=mat->prmatfov[1], mode=mat->mode, itmp;
  unsigned short int x_coordinate, y_coordinate, value, p_coord, lors;
  unsigned int lor, row, fac, ind;
  float ftmp;
  FILE *fp;
  
  //Compose the file name to write
  sprintf(semi_x_c,"%i",semi_x);
  sprintf(semi_y_c,"%i",semi_y);
  //decide the type
  if(mat->type==PRMAT_TYPE_ECAT931) strcpy(matrixfile,"prE");
  else 
    if(mat->type==PRMAT_TYPE_GE) strcpy(matrixfile,"prG");
    else strcpy(matrixfile,"prN");
  //decide the mode 
  sprintf(mod,"%i",mode);
  strcat(matrixfile,mod);
  strcat(matrixfile,semi_x_c);
  strcat(matrixfile,semi_y_c);
  strcat(matrixfile,".prm");
  
  if(access(matrixfile, 0) != -1) {
    if(PRMAT_VERBOSE) printf("PRMAT:File %s already exists.\n", matrixfile);
    return(0);
  }
  
  fp=fopen(matrixfile,"wb+");
  if(fp==NULL) return(1);
 
  // Save status.
  fwrite(&(mat->status),sizeof(int),1,fp);  
 
  // Save type of the scanner.
  fwrite(&(mat->type),sizeof(int),1,fp);  
 
  // Save scanner geometrics.
  fwrite(&(mat->viewNr),sizeof(int),1,fp);  
  fwrite(&(mat->binNr),sizeof(int),1,fp);  
  fwrite(&(mat->imgDim),sizeof(int),1,fp);  
 
  // Save mode.
  fwrite(&(mat->mode),sizeof(int),1,fp);  
 
  // Save FOV.
  fwrite(&(mat->prmatfov[0]),sizeof(int),1,fp);
  fwrite(&(mat->prmatfov[1]),sizeof(int),1,fp); 
 
  // Save min, max, sum and scaling factor.
  fwrite(&(mat->min),sizeof(float),1,fp);
  fwrite(&(mat->max),sizeof(float),1,fp);
  fwrite(&(mat->factor_sum),sizeof(float),1,fp);
  fwrite(&(mat->scaling_factor),sizeof(float),1,fp);  
 
  // If projection matrix is set.
  if(mat->status >= PRMAT_STATUS_BS_OCCUPIED){
 
    fwrite(&(mat->dimr),sizeof(int),1,fp);
    //square sums
    for(row=0;row<prmatGetRows(mat);row++){
      ftmp=prmatGetFactorSqrSum(mat,row);
      fwrite(&ftmp,sizeof(float),1,fp);
    }
    //dime
    for(row=0;row<prmatGetRows(mat);row++){
      itmp=prmatGetPixels(mat,row);
      fwrite(&itmp,sizeof(int),1,fp);
    }
    //the actual data _factdata
    for(row=0;row<prmatGetRows(mat);row++){
      //coordinates and values in every row
      for(fac=0;fac<prmatGetPixels(mat,row);fac++){
  x_coordinate=prmatGetXCoord(mat,row,fac);
  fwrite(&x_coordinate,sizeof(unsigned short int),1,fp);
  y_coordinate=prmatGetYCoord(mat,row,fac);
  fwrite(&y_coordinate,sizeof(unsigned short int),1,fp);
  //NOTE that we wan't to save the values in unsigned short ints
  value=mat->fact[row][fac][2];
  fwrite(&value,sizeof(unsigned short int),1,fp);
      }
    }
  }// END OF SAVING PROJECTION MATRIX.
 
  // If look-up table is set.
  if(mat->status == PRMAT_STATUS_LU_OCCUPIED){
    //save first the statical parts of the structure
    fwrite(&(mat->nrp),sizeof(int),1,fp);
    //save lines: first coordinates of pixels inside the fov,
    for(row=0;row<prmatGetPIX(mat);row++){
      p_coord=prmatGetPixCoord(mat,row);
      fwrite(&p_coord,sizeof(unsigned short int),1,fp);
    }
    // then number of lors hitting a pixel.
    for(row=0; row<prmatGetPIX(mat); row++){
      lors=prmatGetRays(mat,row);
      fwrite(&lors,sizeof(unsigned short int),1,fp);
    }
 
    //save lors
    for(row=0;row<prmatGetPIX(mat);row++){
      for(ind=0; ind<prmatGetRays(mat,row); ind++){
  lor = prmatGetBinView(mat,row,ind);
  fwrite(&lor,sizeof(unsigned int),1,fp);
      }
    }
  }// END OF SAVING LOOK-UP TABLE.
 
  fclose(fp);
  return(0);
}
/*****************************************************************************/
int prmatReadMatrix(char *fname, PRMAT *mat)
{
  unsigned int fov[2], *dimentry=NULL, *nrlor, *iptr, dimtmp, count=0;
  unsigned int row, lor=0, ind, fac;
  unsigned short int  *nrlor_c, *nrlor_l, *usiptr, *usiptr_c;
  unsigned short int x_coordinate=0, y_coordinate=0, value=0;
  float *sqrsum=NULL, *fptr, sqrtmp;
  FILE *fp;
 
  if(PRMAT_VERBOSE) printf("\n PRMAT:prmatReadMatrix(%s,mat) \n",fname);
  if(mat->status==PRMAT_STATUS_UNINITIALIZED) return(1);
  prmatEmpty(mat);
 
  fp=fopen(fname,"rb");
  if(fp==NULL) return(1);
 
  //read first the statical parts of the structure
  if(fread(&(mat->status),sizeof(int),1,fp)==0) {fclose(fp); return(2);}
  if(PRMAT_VERBOSE) printf("PRMAT: status: %i \n",mat->status);
  if(fread(&(mat->type),sizeof(int),1,fp)==0) {fclose(fp); return(2);}
  if(PRMAT_VERBOSE) printf("PRMAT: scanner type: %i \n",mat->type);
  if(fread(&(mat->viewNr),sizeof(int),1,fp)==0) {fclose(fp); return(2);}
  if(fread(&(mat->binNr),sizeof(int),1,fp)==0) {fclose(fp); return(2);}
  if(fread(&(mat->imgDim),sizeof(int),1,fp)==0) {fclose(fp); return(2);}
  if(PRMAT_VERBOSE) 
    printf("PRMAT: scanner geometrics: (%i,%i,%i) \n",
           mat->viewNr,mat->binNr,mat->imgDim);
  if(fread(&(mat->mode),sizeof(int),1,fp)==0) {fclose(fp); return(2);}
  if(PRMAT_VERBOSE) printf("PRMAT: discretisation mode: %i \n",mat->mode);
  if(fread(fov,sizeof(int),1,fp)==0) {fclose(fp); return(2);}
  if(fread((fov+1),sizeof(int),1,fp)==0) {fclose(fp); return(2);}
  mat->prmatfov=calloc(2,sizeof(int)); 
  mat->prmatfov[0]=fov[0];
  mat->prmatfov[1]=fov[1];
  if(PRMAT_VERBOSE) 
    printf("PRMAT: FOV: (%i,%i) \n", mat->prmatfov[0],mat->prmatfov[1]);
  // Read min, max, sum and scaling factor.
  if(fread(&(mat->min),sizeof(float),1,fp)==0) {fclose(fp); return(2);}
  if(PRMAT_VERBOSE) printf("PRMAT: min: %f \n",mat->min);
  if(fread(&(mat->max),sizeof(float),1,fp)==0) {fclose(fp); return(2);}
  if(PRMAT_VERBOSE) printf("PRMAT: max: %f \n",mat->max);
  if(fread(&(mat->factor_sum),sizeof(float),1,fp)==0) {fclose(fp); return(2);}
  if(PRMAT_VERBOSE) printf("PRMAT: sum of all factors: %f \n",mat->factor_sum);
  if(fread(&(mat->scaling_factor),sizeof(float),1,fp)==0) {fclose(fp); return(2);}
  if(PRMAT_VERBOSE) printf("PRMAT: scaling factor: %f \n",mat->scaling_factor);  
 
  // If projection matrix has been occupied.
  if(mat->status >= PRMAT_STATUS_BS_OCCUPIED){
    if(PRMAT_VERBOSE) printf("PRMAT: reading projection matrix \n");  
    if(fread(&(mat->dimr),sizeof(int),1,fp)==0) {fclose(fp); return(2);}
    if(PRMAT_VERBOSE) 
      printf("PRMAT: number of rows in the projection matrix: %i \n",mat->dimr);  
    //read square sums
    sqrsum=calloc(mat->dimr,sizeof(float));  
    for(row=0, fptr=sqrsum;row<(mat->dimr); row++, fptr++) {
      if(fread(&sqrtmp,sizeof(float),1,fp)==0) {fclose(fp); return(2);}
      *fptr=sqrtmp;
    }
    //read dime
    dimentry=(unsigned int*)calloc(mat->dimr,sizeof(int));  
    count = 0;
    for(row=0,iptr=dimentry;row<(mat->dimr);row++,iptr++){
      if(fread(&dimtmp,sizeof(int),1,fp)==0) {fclose(fp); return(2);}
      *iptr=dimtmp;
      count += dimtmp;
    }
 
    //ready to allocate; NOTE that dime is also set in function prmatAllocate()    
    prmatAllocate(mat,1,mat->dimr,dimentry);
    if(PRMAT_VERBOSE)
      printf("PRMAT:prmatAllocate(mat) done in prmatReadMatrix(%s,mat)\n",fname);
 
    //now we can set square sums
    for(row=0, fptr=sqrsum;row<(mat->dimr);row++, fptr++){
      mat->factor_sqr_sum[row]=*fptr;
    }
    //read the actual data _factdata
    for(row=0;row<(mat->dimr);row++){
      //read coordinates and values in every row
      for(fac=0;fac<prmatGetPixels(mat,row);fac++){
  if(fread(&x_coordinate,sizeof(unsigned short int),1,fp)==0) {
          fclose(fp); return(2);}
  mat->fact[row][fac][0]=x_coordinate;
  if(fread(&y_coordinate,sizeof(unsigned short int),1,fp)==0) {
          fclose(fp); return(2);}
  mat->fact[row][fac][1]=y_coordinate;
  if(fread(&value,sizeof(unsigned short int),1,fp)==0) {
          fclose(fp); return(2);}
  mat->fact[row][fac][2]=value;
      }
    }
  }// END OF READING PROJECTIONS
 
  // If look-up table has been occupied.
  if(mat->status == PRMAT_STATUS_LU_OCCUPIED){
    if(PRMAT_VERBOSE) printf("PRMAT: reading look-up table \n");  
    if(fread(&(mat->nrp),sizeof(int),1,fp)==0) {fclose(fp); return(2);}
    if(PRMAT_VERBOSE)
      printf("PRMAT: number of pixels inside the FOV: %i \n",mat->nrp);
    //read nrlor
    nrlor_c=(unsigned short int*)calloc(mat->nrp,sizeof(unsigned short int));  
    nrlor_l=(unsigned short int*)calloc(mat->nrp,sizeof(unsigned short int));  
 
    //first coordinates of pixels.
    usiptr = nrlor_c;
    for(row=0; row<mat->nrp; row++, usiptr++){
      // Read next coordinate from file.
      if(fread(usiptr,sizeof(unsigned short int),1,fp)==0) {fclose(fp); return(2);}
    }
    //then number of lors hitting a pixel
    usiptr = nrlor_l;
    for(row=0; row<mat->nrp; row++, usiptr++){
      // Read next lornr from file.
      if(fread(usiptr,sizeof(unsigned short int),1,fp)==0) {fclose(fp); return(2);}
    }
 
    if(PRMAT_VERBOSE){
      // Get number of lors in lut.
      usiptr = nrlor_l;
      count = 0;
      for(row=0; row<mat->nrp; row++, usiptr++){
  count += *usiptr;
      }
      printf("PRMAT: nr of lors in lut:  %i \n",count);
      //printf("  nr of lors in ave:  %i \n",count/mat->nrp);
    }
 
    // Copy list usi format to integer list.
    nrlor = (unsigned int*)calloc(mat->nrp,sizeof(int));
    usiptr = nrlor_l;
    iptr = nrlor;
    for(row=0; row<mat->nrp; row++, iptr++){
      *iptr = (int)*usiptr++;
    }    
 
    // Now we can allocate memory for look-up table.
    prmatAllocate(mat, 0, mat->nrp, nrlor);
    if(PRMAT_VERBOSE)
      printf("PRMAT:prmatAllocate(mat) done in prmatReadMatrix(%s,mat)\n",fname);
    // And fill mat->lines.
    usiptr_c = nrlor_c;
    usiptr = nrlor_l;
    for(row=0;row<mat->nrp;row++ , usiptr++, usiptr_c++){
      mat->lines[row][0] = *usiptr_c;
      mat->lines[row][1] = *usiptr;
    }
 
    // Read the data _linesdata.
    for(row=0;row<prmatGetPIX(mat);row++){
      for(ind=0; ind<prmatGetRays(mat,row); ind++){
  if(fread(&lor,sizeof(unsigned int),1,fp)==0) {fclose(fp); return(2);}
  mat->lines[row][ind+2] = lor;
      }
    }
 
    free((unsigned short int*)nrlor_l);
    free((unsigned short int*)nrlor_c);
  }// END OF READING LOOK-UP TABLE
 
  if(PRMAT_VERBOSE) printf("PRMAT: prmatReadMatrix() done. \n");
  free((float*)sqrsum); free((int*)dimentry);
 
  fclose(fp);
 
  return(0);
}
/*****************************************************************************/
 
/*****************************************************************************/