reprojection.c

Go to the documentation of this file.

 
/*****************************************************************************/
#include "libtpcrec.h"
/*****************************************************************************/
 
/*****************************************************************************/
int imgReprojection(
  IMG *img,
  IMG *scn,
  int verbose
) {
  if(verbose>0) printf("imgReprojection()\n");
 
  int i, j, ret;
  int plane, frame, binNr, viewNr;
  float *scnData, *imgData, *fptr, f;
  float *sinB, bpZoom, bpZoomInv;
  
 
  /* Check the arguments */
  if(img==NULL || scn==NULL || img->status!=IMG_STATUS_OCCUPIED) return(1);
  /* Clear any previous contents in sinogram data structure */
  imgEmpty(scn);
 
  /* Determine the sinogram dimensions */
  if(img->sizez<3.0) { /* HR+ */
    scn->dimx=binNr=288; scn->dimy=viewNr=144;
    scn->sampleDistance=2.25; /* bin size (mm) */
    scn->axialFOV=155.2; scn->transaxialFOV=583.;
  } else if(img->sizez<5.0) { /* GE Advance */
    scn->dimx=binNr=281; scn->dimy=viewNr=336;
    scn->sampleDistance=1.95730; /* bin size (mm) */
    scn->axialFOV=150.; scn->transaxialFOV=550.;
  } else { /* 931 */
    scn->dimx=binNr=192; scn->dimy=viewNr=256;
    scn->sampleDistance=3.12932; /* bin size (mm) */
    scn->axialFOV=108.; scn->transaxialFOV=600.829;
  }
    
  /*
   *  Allocate output sinogram
   */
  if(verbose>1) printf("allocating memory for sinogram\n");
  ret=imgAllocate(scn, img->dimz, scn->dimy, scn->dimx, img->dimt);
  if(ret) return(3);
  
  /* Set sinogram "header" information */
  scn->type=IMG_TYPE_RAW;
  scn->_fileFormat=img->_fileFormat;
  scn->unit=CUNIT_COUNTS;
  scn->scanStart=img->scanStart;
  scn->sizez=img->sizez;
  strcpy(scn->studyNr, img->studyNr);
  for(plane=0; plane<img->dimz; plane++)
    scn->planeNumber[plane]=img->planeNumber[plane];
  strcpy(scn->radiopharmaceutical, img->radiopharmaceutical);
  scn->isotopeHalflife=img->isotopeHalflife;
  scn->decayCorrection=IMG_DC_UNKNOWN;
  for(frame=0; frame<img->dimt; frame++) {
    scn->start[frame]=img->start[frame]; scn->end[frame]=img->end[frame];
    scn->mid[frame]=0.5*(scn->start[frame]+scn->end[frame]);
  }
  scn->isWeight=0;
 
  /*
   *  Preparations for reprojection
   */
  /* Allocate memory for scan and image data arrays */
  if(verbose>1) printf("allocating memory for matrix data\n");
  scnData=(float*)calloc(scn->dimx*scn->dimy, sizeof(float));
  imgData=(float*)calloc(img->dimx*img->dimy, sizeof(float));
  if(scnData==NULL || imgData==NULL) return(4);
  /* Pre-compute the sine tables for back-projection */
  sinB=(float*)calloc((3*viewNr/2), sizeof(float));
  if(sinB==NULL) {free((char*)scnData); free((char*)imgData); return(4);}
  recSinTables(viewNr, sinB, NULL, 0.0);
  /* Set the backprojection zoom and its inverse */
  bpZoom=img->zoom*(float)img->dimx/(float)binNr; bpZoomInv=1.0/bpZoom;
  if(verbose>1) printf("bpZoom=%g bpZoomInv=%g\n", bpZoom, bpZoomInv);
  /* Initialize variables used by back-projection */
  for(i=0; i<3*viewNr/2; i++) sinB[i]*=bpZoomInv;
 
 
  /*
   *  Reprojection of one matrix at a time
   */
  if(verbose>1) {printf("reprojection of matrices\n"); fflush(stdout);}
  for(plane=0; plane<img->dimz; plane++) for(frame=0; frame<img->dimt; frame++) {
    if(verbose>2) {
      printf("  plane %d frame %d\n", img->planeNumber[plane], frame+1);
      fflush(stdout);
    }
 
    /* Copy image data into the array */
    /* At the same time, correct for the frame length */
    f=img->end[frame]-img->start[frame]; if(f<1.0) f=1.0;
    /* if GE Advance, then multiply by voxel volume, */
    /* because it is corrected again in image calibration */
    if(binNr==281 && img->sizex>0 && img->sizey>0 && img->sizez>0)
      f*=0.001*img->sizex*img->sizey*img->sizez;
    for(i=0, fptr=imgData; i<img->dimy; i++) for(j=0; j<img->dimx; j++)
      *fptr++ = f*img->m[plane][i][j][frame];
 
    /* Initiate sinogram buffer to zero */
    memset(scnData, 0, scn->dimx*scn->dimy*sizeof(float));
 
    /* Reproject on one angle (view) at a time */
    for(i=0; i<viewNr; i++) {
      viewReprojection(imgData, scnData+(i*binNr), i, 
          img->dimx, viewNr, binNr, sinB, sinB, 0.0, 0.0, bpZoom);
      /* smoothing by mean */
      if(scn->dimx>=img->dimx)
        reprojectionAvg5(scnData+(i*binNr), binNr);
      else
        reprojectionAvg3(scnData+(i*binNr), binNr);
    }
    
    /* Copy the sinogram data to sinogram matrix */
    /* Correct for zoom */
    f=bpZoomInv*bpZoomInv;
    for(i=0, fptr=scnData; i<viewNr; i++) for(j=0; j<binNr; j++)
      scn->m[plane][i][j][frame] = *fptr++ * f;
 
  }
  if(verbose>1) {printf("  reprojection done.\n"); fflush(stdout);}
  free(imgData); free(scnData); free(sinB);
  
  return(0);
}
/*****************************************************************************/
 
/*****************************************************************************/
int reprojection(
  float *image,
  int dim,
  int rays,
  int views,
  float bpzoom,
  float *sinogram,
  int verbose
) {
  if(verbose>0) {
    printf("reprojection(%d, %d, %d, %g)\n", dim, rays, views, bpzoom);
    fflush(stdout);
  }
  if(image==NULL || rays<2 || views<2 || dim<2 || sinogram==NULL) return(1);
  if(dim%2) return(2);
  if(bpzoom<0.1) return(3);
 
  int i;
 
  /* Initiate sinogram data to zero */
  for(i=0; i<rays*views; i++) sinogram[i]=0.0;
 
  /* Set the backprojection zoom inverse */
  float bpzoomInv;
  bpzoomInv=1.0/bpzoom;
 
  /* Pre-compute the sine tables for back-projection */
  float sinB[3*views/2];
  recSinTables(views, sinB, NULL, 0.0);
  for(i=0; i<3*views/2; i++) sinB[i]*=bpzoomInv;
 
  /* Reproject on one angle (view) at a time */
  for(i=0; i<views; i++) {
    if(verbose>1) {printf("  reprojecting angle %d\n", 1+i); fflush(stdout);}
    viewReprojection(image, sinogram+(i*rays), i, 
                     dim, views, rays, sinB, sinB, 0.0, 0.0, bpzoom);
 
    if(dim>=rays)
      reprojectionAvg5(sinogram+(i*rays), rays);
    else
      reprojectionAvg3(sinogram+(i*rays), rays);
 
/*
    if(i==0 || i==views/4 || i==views/2 || i==3*views/4)
      reprojectionAvg5(sinogram+(i*rays), rays);
*/
  }
 
  return(0);
}
/*****************************************************************************/
 
/*****************************************************************************/
void viewReprojection(
  float *idata,
  float *sdata,
  int view,
  int dim, 
  int viewNr,
  int rayNr,
  float *sinB,
  float *sinBrot,
  float offsX,
  float offsY,
  float bpZoom
) {
  float *iOrigin, sir, cor, si, co, tpow2, t, *imgp, fract, rayCenter, toffs;
  int x, y, yBottom, xRight, halfdim, ti;
 
  if(view>=viewNr) {
    fprintf(stderr, "Error in viewReprojection(): view=%d >= %d\n", view, viewNr);
    fflush(stderr);
    exit(1);
  }
 
  /* Set pointer to image origin (middle of the image data) */
  halfdim=dim/2;
  iOrigin= idata + dim*(halfdim-1) + halfdim;
  /* Set sine and cosine for x,y-shift */
  si=sinB[view]*offsY; 
  co=sinB[viewNr/2+view]*offsX;
  /* Set sin and cos for reprojection */
  sir=sinBrot[view]; 
  cor=sinBrot[viewNr/2+view];
  /* Set rotation point */
  y=halfdim-2;
  rayCenter=(float)rayNr/2.0;
  if((float)y>rayCenter*bpZoom) {
    y=(int)(rayCenter*bpZoom); 
    yBottom=-y;
  } else {
    yBottom=-halfdim+1;
  }
  toffs=rayCenter-si+co;
  tpow2=rayCenter*rayCenter*bpZoom*bpZoom;
  for(; y>=yBottom; y--) {
    xRight=(int)sqrt(tpow2-((float)y+0.5)*((float)y+0.5)) + 1;
    if(xRight>=halfdim) {
      xRight=halfdim-1; 
      x=-halfdim;
    } else {
      x=-xRight;
    }
    /* distance between projection ray and origo */
    t=toffs - (float)y*sir + ((float)(x+1))*cor;
    imgp=iOrigin-y*dim+x;
    for(; x<=xRight; x++, t+=cor, imgp++) {
      ti=(int)t;
      fract=t-(float)ti;
      sdata[ti] += *imgp * (1.0-fract);
      sdata[ti+1] += *imgp * fract;
    }
  }
}
/*****************************************************************************/
 
/*****************************************************************************/
void reprojectionAvg3(
  float *data,
  int n
) {
  int i;
  float tmp, prev=0.0, w;
 
  w=1.0/3.0;
  for(i=1; i<n-1; i++) {
    tmp=(data[i-1] + data[i] + data[i+1])*w;
    data[i-1]=prev;
    prev=tmp;
  }
  data[i-1]=prev; data[i]=0.0;
}
/*****************************************************************************/
 
/*****************************************************************************/
void reprojectionAvg5(
  float *data,
  int n
) {
  int i;
  float tmp, prev2, prev, w;
 
  prev=prev2=data[2]; w=0.2;
  for(i=2; i<n-2; i++) {
    tmp=(data[i-2] + data[i-1] + data[i] + data[i+1] + data[i+2])*w;
    data[i-2]=prev2;
    prev2=prev;
    prev=tmp;
  }
  data[i-2]=prev2; data[i-1]=prev;
}
/*****************************************************************************/
 
/*****************************************************************************/
void reprojectionMed3(
  float *data,
  int n
) {
  int i;
  float me, prev=0.0, a, b, c;
 
  for(i=1; i<n-1; i++) {
    a=data[i-1]; b=data[i]; c=data[i+1];
    if(a>=b && a<=c) me=a; 
    else if(b>=a && b<=c) me=b;
    else me=c;
    data[i-1]=prev;
    prev=me;
  }
  data[i-1]=prev; data[i]=0.0;
}
/*****************************************************************************/
 
/*****************************************************************************/
void viewBackprojection(
  float *prj,
  float *idata,
  int dim, 
  int view,
  int viewNr,
  int rayNr,
  float *sinB,
  float *sinBrot,
  float offsX,
  float offsY,
  float bpZoom
) {
  int halfdim=dim/2;
  float *iOrigin=idata+dim*(halfdim-1)+halfdim; 
  float si, co, sir, cor;
  /* Set sin and cos for x,y-shift */
  si=sinB[view]*offsY; 
  co=sinB[viewNr/2+view]*offsX;
  /* Set sin and cos for backprojection */
  sir=sinBrot[view];
  cor=sinBrot[viewNr/2+view];
 
  /* Set rotation point */
  float rayCenter, toffs, tpow2;
  int x, y, yBottom, xRight;
  rayCenter=(float)rayNr/2.0;
  y=halfdim-2;
  if((float)y>rayCenter*bpZoom) {
    y=(int)(rayCenter*bpZoom); 
    yBottom=-y;
  } else {
    yBottom=-halfdim+1;
  }
  toffs=rayCenter-si+co;
 
  float *iptr, fract, t;
  int ti;
  tpow2=rayCenter*rayCenter*bpZoom*bpZoom;
  for(; y>=yBottom; y--) {
    xRight=(int)sqrt(tpow2-((float)y+0.5)*((float)y+0.5)) + 1;
    if(xRight>=halfdim) {
      xRight=halfdim-1; 
      x=-halfdim;
    } else {
      x=-xRight;
    }
    /* distance between projection ray and origo */
    t=toffs - (float)y*sir + ((float)(x+1))*cor;
    iptr=iOrigin-y*dim+x;
    for(; x<=xRight; x++, t+=cor, iptr++) {
      ti=(int)t;
      fract=t-(float)ti;
      *iptr+= prj[ti] + (prj[ti+1] - prj[ti])*fract;
    }
  }
}
/*****************************************************************************/
 
/*****************************************************************************/