reprojection.c File Reference

Image reprojection. More...

#include "libtpcrec.h"

Go to the source code of this file.

Functions
int	imgReprojection (IMG *img, IMG *scn, int verbose)
int	reprojection (float *image, int dim, int rays, int views, float bpzoom, float *sinogram, int verbose)
void	viewReprojection (float *idata, float *sdata, int view, int dim, int viewNr, int rayNr, float *sinB, float *sinBrot, float offsX, float offsY, float bpZoom)
void	reprojectionAvg3 (float *data, int n)
void	reprojectionAvg5 (float *data, int n)
void	reprojectionMed3 (float *data, int n)
void	viewBackprojection (float *prj, float *idata, int dim, int view, int viewNr, int rayNr, float *sinB, float *sinBrot, float offsX, float offsY, float bpZoom)

Detailed Description

Image reprojection.

Based on the codes written by Sakari Alenius for Sun UNIX workstations.

Author

Vesa Oikonen

Definition in file reprojection.c.

Function Documentation

imgReprojection()

int imgReprojection	(	IMG *	img,
		IMG *	scn,
		int	verbose )

Image reprojection to 2D sinogram.

Returns

Returns 0 if ok.

Parameters

img	Pointer to IMG struct containing the input image data. img->sizez is used to determine the scanner specific parameters including sinogram dimensions.
scn	Pointer to initiated IMG struct in which the sinogram will be written.
verbose	Verbose level; if zero, then nothing is printed to stderr or stdout

Definition at line 16 of file reprojection.c.

  {
  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);
}

reprojection()

int reprojection	(	float *	image,
		int	dim,
		int	rays,
		int	views,
		float	bpzoom,
		float *	sinogram,
		int	verbose )

Reprojection of one 2D image matrix to sinogram when data is provided as arrays of floats.

See also

imgReprojection, fbp,

Returns

Returns 0 if ok.

Parameters

image	Pointer to float array containing dim*dim image pixel values.
dim	Image x and y dimensions; must be an even number.
rays	Nr of rays (bins or columns) in sinogram data.
views	Nr of views (rows) in sinogram data; usually larger or equal to the image dimensions.
bpzoom	Backprojection zoom factor; imagezoom*dim/rays
sinogram	Pointer to pre-allocated sinogram data; size must be at least rays*views.
verbose	Verbose level; if zero, then nothing is printed to stderr or stdout

Definition at line 154 of file reprojection.c.

  {
  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);
}

Referenced by atnMake().

reprojectionAvg3()

void reprojectionAvg3	(	float *	data,
		int	n )

Average over three samples for image reprojection (1D 3-point mean).

Parameters

data	Pointer to data array
n	Array length

Definition at line 293 of file reprojection.c.

  {
  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;
}

Referenced by imgReprojection(), and reprojection().

reprojectionAvg5()

void reprojectionAvg5	(	float *	data,
		int	n )

Average over five samples for image reprojection (1D 5-point mean).

Parameters

data	Pointer to data array
n	Array length

Definition at line 314 of file reprojection.c.

  {
  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;
}

Referenced by imgReprojection(), and reprojection().

reprojectionMed3()

void reprojectionMed3	(	float *	data,
		int	n )

Median over three samples for image reprojection (1D 3-point median).

Parameters

data	Pointer to data array
n	Array length

Definition at line 336 of file reprojection.c.

  {
  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;
}

viewBackprojection()

void viewBackprojection	(	float *	prj,
		float *	idata,
		int	dim,
		int	view,
		int	viewNr,
		int	rayNr,
		float *	sinB,
		float *	sinBrot,
		float	offsX,
		float	offsY,
		float	bpZoom )

Back-projection of one angle (view)

Parameters

prj	Pointer to source projection data.
idata	Pointer to output image data of size dim*dim (upper left corner).
dim	Image dimension; must be an even number
view	Projection view (in order to get correct sine from tables)
viewNr	Number of projection views (sinogram rows) (to get correct sine from tables)
rayNr	Number of rays (bins, columns)
sinB	Pre-computed sine table for reprojection.
sinBrot	Pre-computed sine table for reprojection with rotation.
offsX	x offset in pixels; shift_x/pixsize.
offsY	y offset in pixels; shift_x/pixsize.
bpZoom	Zoom

Definition at line 359 of file reprojection.c.

  {
  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;
    }
  }
}

Referenced by mrp(), and trmrp().

viewReprojection()

void viewReprojection	(	float *	idata,
		float *	sdata,
		int	view,
		int	dim,
		int	viewNr,
		int	rayNr,
		float *	sinB,
		float *	sinBrot,
		float	offsX,
		float	offsY,
		float	bpZoom )

Reprojection of one angle (view)

Parameters

idata	Pointer to source input image data of size dim*dim
sdata	Pointer to preallocated output sinogram projection view (angle, row)
view	Projection view (in order to get correct sine from tables)
dim	Image dimension; must be an even number
viewNr	Number of projection views (sinogram rows) (to get correct sine from tables)
rayNr	Number of rays (bins, columns)
sinB	Pre-computed sine table for reprojection.
sinBrot	Pre-computed sine table for reprojection with rotation.
offsX	x offset in pixels; shift_x/pixsize.
offsY	y offset in pixels; shift_x/pixsize.
bpZoom	Zoom

Definition at line 216 of file reprojection.c.

  {
  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;
    }
  }
}

Referenced by imgReprojection(), mrp(), reprojection(), and trmrp().