fbp.c

Go to the documentation of this file.

 
/*****************************************************************************/
#include "libtpcrec.h"
/*****************************************************************************/
 
/*****************************************************************************/
int fbp(
  float *sinogram,
  int rays,
  int views,
  int dim,
  int filter,
  float cutoff,
  float *image
) {
  if(sinogram==NULL || rays<2 || views<2 || dim<2 || image==NULL) return(1);
  if(dim%2) return(2);
  if(cutoff<0.05 || cutoff>0.99) return(2);
 
  int i, n;
 
  /* Set FFT array length */
  int lenFFT=(2*rays)-1;
  if(lenFFT<256) lenFFT=256; else if(lenFFT<512) lenFFT=512;
  else if(lenFFT<1024) lenFFT=1024; else if(lenFFT<2048) lenFFT=2048;
  else if(lenFFT<4096) lenFFT=4096; else if(lenFFT<8192) lenFFT=8192;
  else if(lenFFT<16384) lenFFT=16384; else return(3);
 
  /* Set image pixel values to zero */
  for(i=0; i<dim*dim; i++) image[i]=0.0;
 
  /* Allocate memory and set data pointers */
  float *locData=
        calloc((2*lenFFT + (lenFFT+1) + 5*lenFFT/4 + 3*views/2), sizeof(float));
  if(locData==NULL) return(4);
  float *real, *imag;
  float *sinFFT, *cosFFT, *fbpFilter, *sinBP;
  real=locData; imag=real+lenFFT;
  fbpFilter=imag+lenFFT;
  sinFFT=fbpFilter+(lenFFT+1); 
  cosFFT=sinFFT+lenFFT/4; // sinFFT and cosFFT overlap on purpose
  sinBP=sinFFT+5*lenFFT/4;
 
  /* Pre-compute the sine and cosine tables for FFT */
  {
    double df, dg;
    dg=(double)M_PI*2.0/(double)lenFFT;
    for(i=0, df=0.0; i<5*lenFFT/4; i++, df+=dg) sinFFT[i]=(float)sin(df);
  }
  /* Compute reconstruction filter coefficients */
  n=fbpMakeFilter(sinFFT, cosFFT, cutoff, filter, lenFFT, views, fbpFilter, 0); 
  if(n) {free(locData); return(100+n);}
  /* Compute the sine tables for back-projection */
  recSinTables(views, sinBP, NULL, 0.0);
 
  /* Process each sinogram row (projection, view), two at the same time */
  for(int pi=0; pi<views; pi+=2) {
    /* Copy the data to real and imag arrays */
    for(i=0; i<rays; i++) real[i]=sinogram[pi*rays+i];
    for(i=0; i<rays; i++) imag[i]=sinogram[pi*rays+rays+i];
    for(i=rays; i<lenFFT; i++) real[i]=imag[i]=0.0; // zero padding
    /* Forward FFT */
    fbp_fft_bidir_complex_radix2(real, imag, 1, lenFFT, sinFFT, cosFFT);
    /* Filter */
    for(i=0; i<lenFFT; i++) real[i]*=fbpFilter[i];
    for(i=0; i<lenFFT; i++) imag[i]*=fbpFilter[i];
    /* Inverse FFT */
    fbp_fft_bidir_complex_radix2(real, imag, -1, lenFFT, sinFFT, cosFFT);
    /* Back-projection */
    fbp_back_proj_round(real, image+dim*(dim/2-1)+dim/2, dim, pi, views, rays,
                        0.0, 0.0, 1.0, sinBP, sinBP);
    fbp_back_proj_round(imag, image+dim*(dim/2-1)+dim/2, dim, pi+1, views, rays,
                        0.0, 0.0, 1.0, sinBP, sinBP);
  } // next two projections
 
  free(locData);
  return(0);
}
/*****************************************************************************/
 
/*****************************************************************************/
int imgFBP(
  IMG *scn,
  IMG *img,
  int imgDim,
  float zoom,
  int filter,
  float cutoff,
  float shiftX,
  float shiftY,
  float rotation,
  int verbose
) {
  if(verbose>0)
    printf("imgFBP(scn, img, %d, %g, %d, %g, %g, %g, %g)\n",
           imgDim, zoom, filter, cutoff, shiftX, shiftY, rotation);
 
  int i, j, ret, plane, frame, fullBP;
  int lenFFT, halfDim;
  float *sinB, *sinBrot, *sinFFT, *cosFFT;
  float *oReal, *oImag, *real, *imag, *scnData, *imgData, *imgOrigin;
  float *fbpFilter, bpZoom, bpZoomInv, *fptr, *gptr, *row1, *row2, f;
  float offsX, offsY, pixSize; /* note that pixSize is in cm */
 
 
  /* Check the arguments */
  if(scn->status!=IMG_STATUS_OCCUPIED) return(1);
  if(imgDim<2 || imgDim>4096 || imgDim%2) return(1);
  if(zoom<0.01 || zoom>1000.) return(1);
  if(filter<0 || filter>6) return(1);
  if(cutoff<=0.0 || cutoff>=1.0) return(1);
  if(scn->dimx<=1 || scn->dimx>16384) return(1);
 
  /*
   *  Allocate output image
   */
  if(verbose>1) {printf("allocating memory for the image\n"); fflush(stdout);}
  imgEmpty(img); 
  ret=imgAllocate(img, scn->dimz, imgDim, imgDim, scn->dimt);
  if(ret) return(3);
 
  /* Set image "header" information */
  if(verbose>1) {printf("setting image header\n"); fflush(stdout);}
  img->type=IMG_TYPE_IMAGE;
  img->unit=CUNIT_COUNTS; /* This will convert to ECAT counts/sec */
  img->scanStart=scn->scanStart;
  img->axialFOV=scn->axialFOV;
  img->transaxialFOV=scn->transaxialFOV;
  img->sizez=scn->sizez;
  strcpy(img->studyNr, scn->studyNr);
  img->sampleDistance=scn->sampleDistance;
  if(scn->sampleDistance<=0.0) {
    if(scn->dimx==281) {
      img->sizez=4.25;
      scn->sampleDistance=1.95730;
      scn->axialFOV=150.; scn->transaxialFOV=550.;
    } else {
      img->sizez=6.75;
      scn->sampleDistance=3.12932;
      scn->axialFOV=108.; scn->transaxialFOV=600.829;
    }
  }
  pixSize=scn->sampleDistance*(float)scn->dimx/((float)imgDim*zoom);
  img->sizex=img->sizey=pixSize;
  for(plane=0; plane<scn->dimz; plane++)
    img->planeNumber[plane]=scn->planeNumber[plane];
  strcpy(img->radiopharmaceutical, scn->radiopharmaceutical);
  img->isotopeHalflife=scn->isotopeHalflife;
  img->decayCorrection=IMG_DC_NONCORRECTED;
  for(frame=0; frame<scn->dimt; frame++) {
    img->start[frame]=scn->start[frame]; img->end[frame]=scn->end[frame];
    img->mid[frame]=0.5*(img->start[frame]+img->end[frame]);
  }
  img->isWeight=0;
 
  /*
   *  Preparations for reconstruction
   */
  /* Allocate memory for scan and image data arrays */
  if(verbose>1) {printf("allocating memory for matrices\n"); fflush(stdout);}
  i=scn->dimx*scn->dimy + img->dimx*img->dimy;
  float *matData=calloc(i, sizeof(float));
  if(matData==NULL) {
    imgEmpty(img);
    return(4);
  }
  scnData=matData; imgData=matData+(scn->dimx*scn->dimy);
 
  /* Set FFT array length */
  lenFFT=(2*scn->dimx)-1;
  if(lenFFT<256) lenFFT=256; else if(lenFFT<512) lenFFT=512;
  else if(lenFFT<1024) lenFFT=1024; else if(lenFFT<2048) lenFFT=2048;
  else if(lenFFT<4096) lenFFT=4096; else if(lenFFT<8192) lenFFT=8192;
  else if(lenFFT<16384) lenFFT=16384; else return(1);
  if(verbose>2) 
    printf("rays=%d views=%d lenFFT=%d\n", scn->dimx, scn->dimy, lenFFT);
 
  /* Allocate memory for pre-computed tables */
  i=3*scn->dimy/2 + 3*scn->dimy/2 + 5*lenFFT/4;
  if(verbose>1) printf("allocating memory for pre-computed tables\n");
  float *tabData=calloc(i, sizeof(float));
  if(tabData==NULL) {
    imgEmpty(img); free(matData);
    return(4);
  }
  sinB=tabData; sinBrot=tabData+3*scn->dimy/2;
  sinFFT=sinBrot+3*scn->dimy/2;
 
  /* Pre-compute the sine tables for back-projection (note the rotation!) */
  if(verbose>1) printf("computing sine tables for back-projection\n");
  recSinTables(scn->dimy, sinB, sinBrot, rotation);
 
  /* Pre-compute the sine and cosine tables for FFT */
  if(verbose>1) printf("computing sine and cosine tables for FFT\n");
  {
    double df=0.0, dg;
    dg=(double)M_PI*2.0/(double)lenFFT;
    for(i=0; i<5*lenFFT/4; i++, df+=dg) sinFFT[i]=(float)sin(df);
    cosFFT=sinFFT+lenFFT/4;
  }
 
  /* Allocate memory for reconstruction filter coefficients */
  if(verbose>1) printf("allocating memory for reconstruction filters\n");
  float *filData=calloc((5*lenFFT), sizeof(float));
  if(filData==NULL) {
    imgEmpty(img); free(matData); free(tabData);
    return(4);
  }
  fbpFilter=filData; oReal=filData+lenFFT; oImag=oReal+2*lenFFT;
  real=oReal+lenFFT/2; imag=oImag+lenFFT/2;
 
  /* Compute reconstruction filter coefficients */
  ret=fbpMakeFilter(sinFFT, cosFFT, cutoff, filter, lenFFT, scn->dimy, 
                    fbpFilter, verbose-1);
  if(ret) {
    imgEmpty(img); free(matData); free(tabData); free(filData);
    return(5);
  }
  if(verbose>3) {
    printf("\nfbpFilter:\n");
    for(i=0; i<lenFFT; i++) {
      printf("%g ", fbpFilter[i]); 
      if((i+1)%7==0 || i==lenFFT-1) printf("\n");}
    printf("\n");
  }
 
  /* Set the backprojection zoom and inverse (globals) */
  bpZoom=zoom*(float)imgDim/(float)scn->dimx; 
  bpZoomInv=1.0/bpZoom;
  if(verbose>2) printf("bpZoom=%g bpZoomInv=%g\n", bpZoom, bpZoomInv);
 
  /* Check if image corners need to be processed */
  if(zoom>M_SQRT2) fullBP=1; /* profile overlaps the image; 1.41421356 */
  else fullBP=0; /* do not process image corners */
  if(verbose>1) printf("fullBP := %d\n", fullBP);
 
  /* Initialize variables used by back-projection */
  if(verbose>1) printf("initialize variables for back-projection\n");
  halfDim=imgDim/2; 
  offsX=shiftX/pixSize; offsY=shiftY/pixSize;
  for(i=0; i<3*(scn->dimy)/2; i++) {
    sinB[i]*=bpZoomInv; 
    sinBrot[i]*=bpZoomInv;
  }
  imgOrigin=imgData+imgDim*(halfDim-1)+halfDim;
  if(verbose>2) {
    printf("halfDim=%d offsX=%g offsY=%g\n", halfDim, offsX, offsY);
  }
 
  /*
   *  Reconstruct one matrix at a time
   */
  if(verbose>1) printf("reconstruct one matrix at a time...\n");
  for(plane=0; plane<scn->dimz; plane++) 
    for(frame=0; frame<scn->dimt; frame++) {
 
    if(verbose>3) {
      printf("reconstructing plane %d frame %d\n", 
             scn->planeNumber[plane], frame+1);
      fflush(stdout);
    }
 
    /* Copy scan data into the array */
    for(i=0, fptr=scnData; i<scn->dimy; i++) 
      for(j=0; j<scn->dimx; j++)
        *fptr++=scn->m[plane][i][j][frame];
    if(verbose>8) 
      printf("scn->m[%d][158][116][%d]=%e\n",
             plane, frame, scn->m[plane][64][56][frame]);
 
    /* Initiate image buffer */
    for(i=0, fptr=imgData; i<imgDim*imgDim; i++) *fptr++=0.0;
 
    /* Process each sinogram row (projection), two at the same time */
    for(int view=0; view<scn->dimy; view+=2) {
 
      /* Fill FFT buffers with two rows */
      for(j=0; j<2*lenFFT; j++) oReal[j]=0.0; /* zero padding */
      for(j=0; j<2*lenFFT; j++) oImag[j]=0.0; /* zero padding */
      row1=scnData+view*scn->dimx; 
      row2=row1+scn->dimx;
      fptr=real; gptr=imag;
      for(j=0; j<scn->dimx; j++) {
        *fptr++ = *row1++;
        *gptr++ = *row2++;
      }
 
      /* Forward FFT */
      fptr=real; gptr=imag;
      fbp_fft_bidir_complex_radix2(fptr, gptr, 1, lenFFT, sinFFT, cosFFT);
 
      /* Due to symmetry properties, filtering the combined data together */
      for(j=0, fptr=fbpFilter; j<lenFFT; j++) {
        real[j] *= *fptr;
        imag[j] *= *fptr++;
      }
 
      /* Inverse FFT */
      fptr=real; gptr=imag;
      fbp_fft_bidir_complex_radix2(fptr, gptr, -1, lenFFT, sinFFT, cosFFT);
      fptr=real; gptr=imag;
      if(fullBP) {
        fbp_back_proj(fptr, imgData, imgDim, view, scn->dimy, scn->dimx,
                   offsX, offsY, sinB, sinBrot);
        fbp_back_proj(gptr, imgData, imgDim, view+1, scn->dimy, scn->dimx,
                   offsX, offsY, sinB, sinBrot);
      } else {
        fbp_back_proj_round(fptr, imgOrigin, imgDim, view,
                 scn->dimy, scn->dimx, offsX, offsY, bpZoom, sinB, sinBrot);
        fbp_back_proj_round(gptr, imgOrigin, imgDim, view+1,
                 scn->dimy, scn->dimx, offsX, offsY, bpZoom, sinB, sinBrot);
      }
 
    } /* next projection (view) */
 
    /* Copy the image array to matrix data            */
    /* At the same time, correct for the frame length */
    f=img->end[frame]-img->start[frame]; 
    if(f<1.0) f=1.0; 
    f=1.0/f;
    for(i=0, fptr=imgData; i<img->dimy; i++)
      for(j=0; j<img->dimx; j++)
        img->m[plane][i][j][frame] = (*fptr++)*f;
    if(verbose>8) 
      printf("img->m[%d][64][56][%d]=%e\n",
             plane, frame, img->m[plane][64][56][frame]);
 
  } /* next matrix */
 
  /* Free the memory */
  free(matData); free(tabData); free(filData);
 
  if(verbose>1) printf("imgFBP() done.\n");
  return(0);
}
/*****************************************************************************/
 
/*****************************************************************************/
int fbpMakeFilter(
  float *sinFFT,
  float *cosFFT,
  float cutoff,
  int filter,
  int lenFFT,
  int views,
  float *fbpFilter,
  int verbose
) {
  if(verbose>0) 
    printf("fbpMakeFilter(..., %g, %d, %d, %d, ...)\n", 
           cutoff, filter, lenFFT, views);
  if(sinFFT==NULL || cosFFT==NULL || fbpFilter==NULL) return(1);
 
  int i, istop;
  float *ramp, f;
  float invNr, bp_ifft_scale, rampDC, rampSlope;
 
  /* Check arguments */
  if(lenFFT<2 || views<2) return(1);
 
  /*
   *  Make the ramp to be used in building the reconstruction filter.
   *  Multiplied later by a window to get the actual filter.
   *  Ramp filter:
   *    filt[i] = i/N        0 <= i <= N/2
   *    filt[i] = (N-i)/N    N/2+1 <= i <= N-1
   */
  float *locData=calloc(lenFFT, sizeof(float));
  if(locData==NULL) return(2);
  ramp=locData;
 
  invNr=1.0/(float)lenFFT; 
  bp_ifft_scale=M_PI/(float)(views*lenFFT);
  if(verbose>1) printf("invNr=%g bp_ifft_scale=%g\n", invNr, bp_ifft_scale);
  /* Construct the direct method (original VAX) to calculate the ramp function:    */
  ramp[0]=0.0;
  for(i=1; i<=lenFFT/2; i++) 
    ramp[i]=ramp[lenFFT-i]=(float)i*invNr;
  rampDC=ramp[0]; 
  rampSlope=ramp[lenFFT/2]-rampDC; // height of ramp
  if(verbose>1) printf("rampDC=%g rampSlope=%g\n", rampDC, rampSlope);
  /* FFT filter is made by multiplying a ramp by a window function, */
  /* according to the filter type                                   */
  istop=(int)(cutoff*(float)lenFFT+0.5); 
  if(istop>lenFFT/2) istop=lenFFT/2;
 
  switch(filter) {
    case FBP_FILTER_RAMP: /* Ramp */
      fbpFilter[0]=ramp[0]*bp_ifft_scale;
      for(i=1; i<=istop; i++) 
        fbpFilter[i]=fbpFilter[lenFFT-i]=ramp[i]*bp_ifft_scale;
      for(; i<=lenFFT/2; i++) 
        fbpFilter[i]=fbpFilter[lenFFT-i]=0.0;
      break;
    case FBP_FILTER_HANN: /* Hann */
      fbpFilter[0]=ramp[0]*bp_ifft_scale;
      for(i=1; i<=istop; i++) {
        f=0.5+0.5*cosf(M_PI*invNr*(float)i/cutoff);
        fbpFilter[i]=fbpFilter[lenFFT-i]=ramp[i]*f*bp_ifft_scale;
      }
      for(; i<=lenFFT/2; i++) fbpFilter[i]=fbpFilter[lenFFT-i]=0.0;
      break;
    case FBP_FILTER_HAMM: /* Hamm */
      fbpFilter[0]=ramp[0];
      for(i=1; i<=istop; i++) {
        f=0.54+0.46*cosf(M_PI*invNr*(float)i/cutoff);
        fbpFilter[i]=fbpFilter[lenFFT-i]=ramp[i]*f*bp_ifft_scale;
      }
      for(; i<=lenFFT/2; i++) fbpFilter[i]=fbpFilter[lenFFT-i]=0.0;
      break;
    case FBP_FILTER_PARZEN: /* Parzen */
      fbpFilter[0]=ramp[0];
      for(i=1; i<=istop/2; i++) {
        f=(float)i/(float)istop; f=1.0-6.0*f*f*(1.0-f);
        fbpFilter[i]=fbpFilter[lenFFT-i]=ramp[i]*f*bp_ifft_scale;
      }
      for(; i<=istop; i++) {
        f=1.0-(float)i/(float)istop; f=2.0*f*f*f;
        fbpFilter[i]=fbpFilter[lenFFT-i]=ramp[i]*f*bp_ifft_scale;
      }
      for(; i<=lenFFT/2; i++) fbpFilter[i]=fbpFilter[lenFFT-i]=0.0;
      break;
    case FBP_FILTER_SHEPP: /* Shepp */
      fbpFilter[0]=ramp[0];
      for(i=1; i<=istop; i++) {
        f=0.5*(float)i/(float)istop; f=sinf(f*M_PI)/(f*M_PI);
        fbpFilter[i]=fbpFilter[lenFFT-i]=ramp[i]*f*bp_ifft_scale;
      }
      for(; i<=lenFFT/2; i++) fbpFilter[i]=fbpFilter[lenFFT-i]=0.0;
      break;
    case FBP_FILTER_BUTTER: /* Butter */
      fprintf(stderr, "Warning: Butter filter not implemented, using None.\n");
      /* __attribute__((fallthrough)); */
      /* FALLTHROUGH */
    case FBP_FILTER_NONE: /* None */
      for(i=0; i<lenFFT; i++) fbpFilter[i]=bp_ifft_scale;
      rampSlope=0.0; /* Reset the ramp slope value */
      break;
    default:
      fprintf(stderr, "Error: invalid filter requested.\n");
      free(locData); 
      return(1);
  } /* switch */
 
  free(locData);
  return(0);
}
/*****************************************************************************/
 
/*****************************************************************************/
void fbp_fft_bidir_complex_radix2(
  float *real, 
  float *imag, 
  int direct,
  int n,
  float *sine, 
  float *cosine
) {
  int i, j, m, halfn, mmax, istep, ind;
  float c, s, treal, timag;
 
  halfn=n/2;
 
 
#if(0)
  for(i=j=1; i<=n; i++) {
    if(i<j) {
      treal=real[j-1]; timag=imag[j-1]; 
      real[j-1]=real[i-1]; imag[j-1]=imag[i-1]; 
      real[i-1]=treal; imag[i-1]=timag;
    }
    m=halfn; 
    while(j>m) {j-=m; m++; m/=2;} 
    j+=m;
  }
 
  mmax=1;
  while(mmax<n) {
    istep=2*mmax;
    for(m=1; m<=mmax; m++) {
      ind=(int)((float)(halfn*(m-1))/(float)mmax);
      c=cosine[ind]; s=(direct>=0)?sine[ind]:-sine[ind];
      for(i=m; i<=n; i+=istep) {
        j=i+mmax;
        treal=real[j-1]*c - imag[j-1]*s;
        timag=imag[j-1]*c + real[j-1]*s;
        real[j-1]=real[i-1]-treal;
        imag[j-1]=imag[i-1]-timag;
        real[i-1]+=treal;
        imag[i-1]+=timag;
      }
    }
    mmax=istep;
  }
 
#else
  for(i=j=0; i<n; i++) {
    if(i<j) {
      treal=real[j]; timag=imag[j]; 
      real[j]=real[i]; imag[j]=imag[i]; 
      real[i]=treal; imag[i]=timag;
    }
    m=halfn; 
    while(j>=m) {j-=m; m++; m/=2;} 
    j+=m;
  }
 
  mmax=1;
  while(mmax<n) {
    istep=2*mmax;
    for(m=0; m<mmax; m++) {
      ind=(int)((float)(halfn*m)/(float)mmax);
      c=cosine[ind]; s=(direct>=0)?sine[ind]:-sine[ind];
      for(i=m; i<n; i+=istep) {
        j=i+mmax;
        treal=real[j]*c - imag[j]*s;
        timag=imag[j]*c + real[j]*s;
        real[j]=real[i]-treal;
        imag[j]=imag[i]-timag;
        real[i]+=treal;
        imag[i]+=timag;
      }
    }
    mmax=istep;
  }
#endif
}
/*****************************************************************************/
 
/*****************************************************************************/
void fbp_back_proj(
  float *prj,
  float *imgCorner,
  int imgDim,
  int view,
  int views,
  int rays,
  float offsX,
  float offsY,
  float *sinB,
  float *sinBrot
) {
  float si, co, tleft, sir, cor, rayCenter, dif;
  float t; /* distance between projection ray and origin */
  int x, y, halfDim, leftX, rightTopLimit, yBottom, ti;
 
  float *imgp=imgCorner;
 
  halfDim=imgDim/2;
  leftX=-halfDim; 
  rightTopLimit=halfDim-1;
  yBottom=-halfDim;
  rayCenter=0.5*(float)rays; 
  si=sinB[view]; /* Sine and cosine for x,y-shift */
  co=sinB[view+views/2];
  sir=sinBrot[view]; /* Sine and cosine for backprojection */
  cor=sinBrot[view+views/2];
  x=leftX; y=rightTopLimit;
  tleft= rayCenter - (float)y*sir - offsY*si + ((float)(x+1))*cor + offsX*co;
  for(; y>=yBottom; y--) {
    t=tleft;
    for(x=leftX; x<=rightTopLimit; x++, t+=cor) {
      ti=(int)t; dif=t-(float)ti;
      //ti=(int)(t+0.5); *imgp++ += prj[ti];
      *imgp++ += prj[ti] + (prj[ti+1]-prj[ti])*dif;
    }
    tleft+=sir;
  }
}
/*****************************************************************************/
 
/*****************************************************************************/
void fbp_back_proj_round(
  float *prj,
  float *imgOrigin,
  int imgDim,
  int view,
  int views,
  int rays,
  float offsX,
  float offsY,
  float bpZoom,
  float *sinB,
  float *sinBrot
) {
  if(0) printf("fbp_back_proj_round(..., %d, %d, %d, %d, %g, %g, %g, ...)\n",
               imgDim, view, views, rays, offsX, offsY, bpZoom);
 
  float *imgp, si, co, sir, cor, toffs, rayCenter, tPow2, yph, dif;
  float t; /* distance between projection ray and origo */
  int x, y, halfDim, right_top_limit, xrightround, ybottomround, ti;
 
  /* Sine and cosine for x,y-shift */
  si=sinB[view]*offsY; 
  co=sinB[views/2+view]*offsX;
  /* Sine and cosine for backprojection */
  sir=sinBrot[view]; 
  cor=sinBrot[views/2+view];
  halfDim=imgDim/2; 
  y=halfDim-2;
  rayCenter=0.5*(float)rays;
  if((float)y>(rayCenter*bpZoom)) {
    y=(int)(rayCenter*bpZoom);
    ybottomround=-y;
  } else {
    ybottomround=-halfDim+1;
  }
  toffs=rayCenter-si+co;
  tPow2=rayCenter*rayCenter*bpZoom*bpZoom;
  right_top_limit=halfDim-1;
  for(; y>=ybottomround; y--) {
    yph=0.5+(float)y;
    xrightround=(int)sqrt(tPow2 - yph*yph) + 1;
    if(xrightround>=halfDim) {
      xrightround=right_top_limit; 
      x=-halfDim;
    } else {
      x=-xrightround;
    }
    t=toffs-(float)y*sir+((float)(x+1))*cor;
    imgp=imgOrigin-y*imgDim+x;
    for(; x<=xrightround; x++, t+=cor) {
      ti=(int)t; dif=t-(float)ti;
      *imgp++ += prj[ti] + (prj[ti+1] - prj[ti])*dif;
    }
  }
}
/*****************************************************************************/
 
/*****************************************************************************/