radon.c

Go to the documentation of this file.

/*****************************************************************************/
int RADON_TEST;
int RADON_VERBOSE;
/*****************************************************************************/
#include "libtpcrec.h"
/*****************************************************************************/
 
/*****************************************************************************/
/* Initialization and memory handling for radon transform data. */
/*****************************************************************************/
void radonEmpty(RADON *radtra)
{
  if(RADON_VERBOSE){
   printf("RADON: radonEmpty() started. \n");
   fflush(stdout);
  }
  //if radtra is already empty
  if(radtra->status<RADON_STATUS_INITIALIZED){
    if(RADON_VERBOSE){
      printf("RADON: radon object already empty: status %i \n",radtra->status);
      fflush(stdout);
    }
    //return;
  }
  //free the memory occupied by sine 
  free(radtra->sines);
 
  if(RADON_VERBOSE){
    printf("RADON: radonEmpty() finished. \n");
    fflush(stdout);
  }
  return;
}
/*****************************************************************************/
int radonSet(RADON *radtra,int mode,int imgDim,int viewNr,int binNr)
{
  int i;
 
  if( RADON_VERBOSE ){
    printf("RADON: radonSet() started.\n");
    fflush(stdout);
  }
 
  // Set the parameters.
  radtra->mode=mode;
  
  if(imgDim <= 0) return -1;
  else radtra->imgDim=imgDim;
 
  if(viewNr <= 0) return -2;
  else radtra->viewNr=viewNr;
 
  if(binNr <= 0) return -3;
  else radtra->binNr=binNr;
 
  radtra->sampleDist=(float)imgDim/(float)binNr;
 
  // Calculate and set center bin for current geometrics.
  if((binNr%2) != 0){
    radtra->half = (binNr - 1)/2 + 1;
    radtra->centerBin = radtra->half - 1;
  } else {
    radtra->half = binNr/2;
    // In the case binNr is even there is no center bin.
    radtra->centerBin = -1;
  }
 
  /* Set the sine table to contain values of sine to cover the required values 
     of cosine as well. */
  radtra->sines=(float*)calloc(3*viewNr/2,sizeof(float));
  if(radtra->sines==NULL) return -5;
  //Put the values sin(view*pi/viewNr) for view=0:3*viewNr/2 - 1 in the table
  for(i=0; i< 3*viewNr/2; i++) {
    radtra->sines[i]=(float)sin((M_PI/(double)viewNr) * (double)i);
  }
  radtra->status=RADON_STATUS_INITIALIZED;
 
  if( RADON_VERBOSE ){
    printf("RADON: radonSet() done.\n");
    fflush(stdout);
  }
  return 0;
}
/*****************************************************************************/
/* Get functions for Radon data */
/*****************************************************************************/
int radonGetMO(RADON *radtra)
{
  return radtra->mode;
}
/*****************************************************************************/
int radonGetID(RADON *radtra)
{
  return radtra->imgDim;
}
/*****************************************************************************/
int radonGetNV(RADON *radtra)
{
  return radtra->viewNr;
}
/*****************************************************************************/
int radonGetNB(RADON *radtra)
{
  return radtra->binNr;
}
/*****************************************************************************/
float radonGetSD(RADON *radtra)
{
  return radtra->sampleDist;
}
/*****************************************************************************/
int radonGetHI(RADON *radtra)
{
  return radtra->half;
}
/*****************************************************************************/
int radonGetCB(RADON *radtra)
{
  return radtra->centerBin;
}
/*****************************************************************************/
float radonGetSin(RADON *radtra, int nr)
{
  return radtra->sines[nr];
}
/*****************************************************************************/
/* THE ACTUAL TRANSFORM FUNCTIONS */
/*****************************************************************************/
int radonFwdTransform(
  RADON *radtra, int set, int setNr, float *imgdata, float *scndata
) {
  float *imgptr, *scnptr; // pointers for the image and sinogram data
  float *Y, *X, *Xptr, *Yptr; // pointers for the values of LOR in integer points
  double sinus, cosinus, tanus; // sin, cosine and tangent
  double shift_x, shift_y; // shift for LOR
  double scalef; // scaling factor for angle
  int xp, xn, yp, yn, z; //integer points
  int x_left, x_right, y_top, y_bottom; // limits for fast search
  int col, row, view, bin; //counters
  int binNr, viewNr, imgDim, mode;
  int half, center = -1;
  float sdist;
  float dx, dy , loi = 1;  // delta x and y and length of intersection
 
  if( RADON_VERBOSE ){
    printf("RADON: radonFwdTransform() started.\n");
    fflush(stdout);
  }
 
  // Check that the data for this radon transform is initialized.
  if(radtra->status != RADON_STATUS_INITIALIZED) return -1;
 
  // Retrieve the parameters from given radon transform object.
  mode=radonGetMO(radtra);
  imgDim=radonGetID(radtra);
  binNr=radonGetNB(radtra);
  viewNr=radonGetNV(radtra);
  sdist=radonGetSD(radtra); 
  half=radonGetHI(radtra);
  center=radonGetCB(radtra);
 
  // Array for storing values f_x.
  X=(float*)calloc(imgDim+1,sizeof(float));  
  // Array for storing values f_y.
  Y=(float*)calloc(imgDim+1,sizeof(float));  
  if(X==NULL || Y==NULL){
    return -2;
  }
  Xptr=X;
  Yptr=Y;
 
  imgptr=imgdata;
  scnptr=scndata;

  // Find pixel coordinates of the contributing pixels for lines of response
  // belonging to first 1/4th of angles. From these pixel coordinates
  // others can be calculated via symmetries in projection space.
  // N.B. The line of response: a*x+b*y+c
  // => solve y: y = (s - x*cos(theta))/sin(theta)
  //    solve x: x = (s - y*sin(theta))/cos(theta)
 
  for(view=set; view<=viewNr/4; view=view+setNr){
    // Choose the type of the line of response according to view number.
 
    // view=0 -> sin(theta)=0
    if(view==0){
 
      // Length of intersection is 1.
      loi = 1.;
 
      // Choose column according to sample distance for angles 0 and pi/2.
      col = 0;
      for(bin=0; bin<half; bin++){
  col = floor((float)(bin+.5*sdist)*sdist); 
 
  /* Iterate through the entries in the image matrix.
     Calculate raysums for two LORs in the same distance from origin 
           (do halfturn). */
  for(row=0; row<imgDim; row++) {
   
    scnptr[bin] += loi * imgptr[row*imgDim + col];
    if(bin != center)
      scnptr[binNr-bin-1] += loi * imgptr[row*imgDim + (imgDim - 1 - col)];
 
    scnptr[binNr*(viewNr/2) + bin] += 
            loi * imgptr[(imgDim - 1 - col)*imgDim + row];
    if(bin != center)
      scnptr[binNr*(viewNr/2) + (binNr-bin-1)] += 
              loi * imgptr[col*imgDim + row]; 
  }
 
      }
      // End of view==0 (handles angles 0 and pi/2).   
    } else { // angles != 0
 
      // Set sine and cosine for this angle.
      sinus = (double)radonGetSin(radtra,view);
      cosinus = (double)radonGetSin(radtra,viewNr/2 + view);
      tanus = sinus/cosinus;
 
      // Set shift from origin for the first line of response (-n/2,theta).
      // NOTE that image matrix is on cartesian coordinate system where origin
      // is in the middle and shift is in pixels.
      shift_y = -(imgDim/2 -.5*sdist)/sinus;
      shift_x = -(imgDim/2 -.5*sdist)/cosinus;
 
      // Evaluate the function of the first LOR in integer points [-n/2,n/2].
      // NOTE that image matrix is on cartesian coordinate system where origin
      // is in the middle.
      z=-imgDim/2;
      for(col=0; col<imgDim+1; col++){
  Yptr[col]=(float)(shift_y - z/tanus);
    Xptr[col]=(float)(shift_x - z*tanus);
  z++;
      }
 
      // Set shift from the first LOR.
      shift_y = (double)(sdist/sinus);
      shift_x = (double)(sdist/cosinus);
      
      // Set scaling for angle.
      scalef = sinus + cosinus;
 
      // Iterate through half the bins in this view,
      // and determine coordinates of pixels contributing to this LOR.
      // NOTE that shift is added according to 'bin' in every loop.
      // Calculate also the length of intersection.
      // Others are determined via symmetry in projection space.
 
      for(bin=0; bin<half; bin++) {
 
  /* Set the number of intersected pixels to zero. */
        /* np = 0; */
 
  // Limit (x-)indices for fast search.
  // Note that indices are non-negative integers.
  x_left = floor((float)(Xptr[imgDim] + bin*shift_x + imgDim/2));
  if(x_left < 0) x_left = 0;
 
  x_right = floor((float)(Xptr[0] + bin*shift_x + imgDim/2));
  if(x_right <= 0) x_right = 1;
  if(x_right > imgDim) x_right = imgDim - 1; 
 
  /* Iterate through the values in vector Y, in integer points 
           [x_left,x_rigth]. */
  for(z=x_left; z <= x_right; z++) {
       
    xp = z; //positive x-coordinate
    xn = imgDim - 1 - xp; //negative x-coordinate
 
    // Look y from left.  yp=positive y-coordinate
    yp = imgDim/2 - floor(Yptr[xp + 1] + bin*shift_y) - 1;
    yn = imgDim - 1 - yp;
 
    // If the y-value for this x (z) is inside the image grid.
    if(yp < imgDim && yp >= 0 && yn < imgDim && yn >= 0 && xp >= 0 && 
             xp < imgDim && xn < imgDim && xn >= 0)
          {
 
      if(!mode) {       
        loi = 1;
      } else {
        // Compute delta x and y.
        dx = fabs((float)(xp + 1 - (imgDim/2 + Xptr[imgDim - yp] + 
                        bin*shift_x)));       
        dy = fabs((float)(floor(Yptr[xp + 1] + bin*shift_y) + 1 - 
                        (Yptr[xp + 1] + bin*shift_y)));
        if(dx > 1 || dx < 0) dx = 1;
        if(dy > 1 || dy < 0) dy = 1;
        loi = sqrt(dx*dx + dy*dy);
      }
      // Case: theta.
      // Add img(x,y)*k to the raysum of LOR (view,bin)
      scnptr[view*binNr + bin] += loi * imgptr[yp*imgDim + xp];
      if(bin != center)
        // Add img(-x,-y)*k to the raysum of LOR (view,binNr-bin)
        scnptr[view*binNr + binNr - 1 - bin] += 
                loi * imgptr[yn*imgDim + xn];
 
      if(view != viewNr/4) {        
        // Mirror the original LOR on y-axis, i.e. x->-x
        // Case: pi-theta.
        // Add img(-x,y)*k to the raysum of LOR (viewNr-view,bin)
        scnptr[(viewNr - view)*binNr + bin] += 
                loi * imgptr[yp*imgDim + xn];
        // Add img(x,-y)*k to the raysum of LOR (viewNr-view,binNr-bin)
        if(bin != center)
          scnptr[(viewNr-view)*binNr + binNr - 1 - bin] += 
                  loi * imgptr[yn*imgDim + xp];
        
        // Mirror the LOR on line x=y, i.e. x->y.
        // Case: pi/2-theta
        // Add img(-y,-x)*k to the raysum of LOR (viewNr/2-view,bin)
        scnptr[(viewNr/2 - view)*binNr + bin] += 
                loi * imgptr[xn*imgDim + yn];
        // Add img(y,x)*k to the raysum of LOR (viewNr/2-view,binNr-bin)
        if(bin != center)
          scnptr[(viewNr/2 - view)*binNr + binNr - 1 - bin] += 
                  loi * imgptr[xp*imgDim + yp];
      }
 
      // Mirror the LOR on line x=y, and on y-axis i.e x->y and y->-y.
      // Case: pi/2+theta
      // Add img(y,-x)*k to the raysum of LOR (viewNr/2+view,bin)
      scnptr[(viewNr/2 + view)*binNr + bin] +=
              loi * imgptr[xn*imgDim + yp];
      // Add img(y,-x)*k to the raysum of LOR (viewNr/2+view,binNr-bin)
      if(bin != center)
        scnptr[(viewNr/2 + view)*binNr + binNr - 1 - bin] +=
                loi * imgptr[xp*imgDim + yn];       
    }
  }
 
  // Limit (y-)indices for fast search.
  // Note that indices are non-negative integers.
  y_bottom = floor((float)(Yptr[imgDim] + bin*shift_y + imgDim/2));
  if(y_bottom < 0) y_bottom = 0;
  if(y_bottom > imgDim) y_bottom = 0; 
 
  y_top = floor((float)(Yptr[0] + bin*shift_y + imgDim/2));
  if(y_top > imgDim) y_top = imgDim-1;
  if(y_top <= 0) y_top = 1;
 
  /* Iterate through the values in vector X, in integer points 
           [y_bottom,y_top]. */
  for(z=y_top; z >= y_bottom; z--) {
   
    // Look y from this location.
    xp = floor(Xptr[z] + bin*shift_x) + imgDim/2;
    xn = imgDim - 1 - xp;
 
    yp = imgDim - z - 1;
    yn = imgDim - yp - 1;
 
    // If the x-value for this y (z) is inside the image grid.
    if(yp < imgDim && yp >= 0 && yn < imgDim && yn >= 0 && xp >= 0 && 
             xp < imgDim && xn < imgDim && xn >= 0)
          {
     
      if(!mode) {
        loi = 1;
      } else {
        dx = (float)(Xptr[z] + bin*shift_x + imgDim/2 - xp);
        dy = (float)(Yptr[xp] + bin*shift_y - z + imgDim/2); 
        if(dy > 1 || dy < 0){
    dx = dx - (float)(Xptr[z+1] + bin*shift_x + imgDim/2 - xp);
    dy = 1;
        }
        loi = sqrt(dx*dx + dy*dy);
      }
      
      // Case: theta.
      // Add img(x,y)*k to the raysum of LOR (view,bin)
      scnptr[view*binNr + bin] += loi * imgptr[yp*imgDim + xp];
      // Add img(-x,-y)*k to the raysum of LOR (view,binNr-bin)
      if(bin != center)
        scnptr[view*binNr + binNr - 1 - bin] += 
                loi * imgptr[yn*imgDim + xn];
 
      if(view != viewNr/4) {              
        // Mirror the LOR on y-axis, i.e. x->-x
        // Case: pi-theta.
        // Add img(-x,y)*k to the raysum of LOR (viewNr-view,bin)
        scnptr[(viewNr - view)*binNr + bin] += loi * imgptr[yp*imgDim + xn];
        // Add img(x,-y)*k to the raysum of LOR (viewNr-view,binNr-bin)
        if(bin != center)
    scnptr[(viewNr-view)*binNr + binNr - 1 - bin] += 
                  loi * imgptr[yn*imgDim + xp];
        
        // Mirror the LOR on line y=x, i.e. y->x.
        // Case: pi/2 - theta.
        // Add img(y,x)*k to the raysum of LOR (viewNr/2-view,bin)
        scnptr[(viewNr/2 - view)*binNr + bin] +=
                loi * imgptr[xn*imgDim + yn];
        // Add img(-y,-x)*k to the raysum of LOR (viewNr/2-view,binNr-bin)
        if(bin != center)
    scnptr[(viewNr/2 - view)*binNr + binNr - 1 - bin] +=
                  loi * imgptr[xp*imgDim + yp];
      }
 
      // Mirror the LOR on line x=y, and on y-axis i.e x->y and y->-y.
      // Case: pi/2 + theta.
      // Add img(y,-x)*k to the raysum of LOR (viewNr/2+view,bin)
      scnptr[(viewNr/2 + view)*binNr + bin] += loi * imgptr[xn*imgDim + yp];
      // Add img(-y,x)*k to the raysum of LOR (viewNr-view,binNr-bin)
      if(bin != center)
        scnptr[(viewNr/2 + view)*binNr + binNr - 1 - bin] +=
                loi * imgptr[xp*imgDim + yn];
    }
  }
 
  // If mode==0 scale with sin(theta)+cos(theta).
  if(mode==0) {
    // Case: theta.
    scnptr[view*binNr + bin] /= scalef;
    if(bin != center)
      scnptr[view*binNr + binNr - 1 - bin] /= scalef;
 
    if(view != viewNr/4){             
      // Mirror the LOR on y-axis, i.e. x->-x
      // Case: pi-theta.
      scnptr[(viewNr - view)*binNr + bin] /= scalef;
      if(bin != center)
        scnptr[(viewNr-view)*binNr + binNr - 1 - bin] /= scalef;
        
      // Mirror the LOR on line y=x, i.e. y->x.
      // Case: pi/2 - theta.
      scnptr[(viewNr/2 - view)*binNr + bin] /= scalef;
      if(bin != center)
        scnptr[(viewNr/2 - view)*binNr + binNr - 1 - bin] /= scalef;
    }
 
    // Mirror the LOR on line x=y, and on y-axis i.e x->y and y->-y.
    // Case: pi/2 + theta.
    scnptr[(viewNr/2 + view)*binNr + bin] /= scalef;
    if(bin != center)
      scnptr[(viewNr/2 + view)*binNr + binNr - 1 - bin] /= scalef;
  }
      }// END of X loop
    }// End of view>0.
  }// END OF VIEW LOOP
  free(X);
  free(Y);
 
  if(RADON_VERBOSE) {
    printf("RADON: radonFwdTransform() finished.\n");
    fflush(stdout);
  }
 
  return 0;
}// END OF FORWARD (spa -> pro) RADON TRANSFORM
/*****************************************************************************/
int radonFwdTransformEA(
  RADON *radtra, int set, int setNr, float *imgdata, float *scndata)
{
  float *imgptr, *scnptr; // pointers for the image and sinogram data
  float *Y, *X, *Xptr, *Yptr; // pointers for the values of LOR in integer points
  double sinus, cosinus, tanus; // sin, cosine and tangent
  double shift_x, shift_y; // shift for LOR
  int xp, xn, yp, yn, z, xp2; //integer points
  int x_left, x_right, y_top, y_bottom; // limits for fast search
  int col, col1, col2, row, view, bin; //counters
  int binNr, viewNr, imgDim, errors=0;
  int half, center = -1;
  float sdist;
  float a=0,b=0, c=0, d=0, g=0, h=0, A;
  float dx, dy , eps1 = 0, eps2 = 0, eps3 = 0;  // delta x and y and factors.
 
  if( RADON_VERBOSE ){
    printf("RADON: radonFwdTransformEA() started.\n");
    fflush(stdout);
  }
 
  // Check that the data for this radon transform is initialized.
  if(radtra->status != RADON_STATUS_INITIALIZED) return -1;
 
  // Retrieve the parameters from given radon transform object.
  //mode=radonGetMO(radtra); // Should be 2.
  imgDim=radonGetID(radtra);
  binNr=radonGetNB(radtra);
  viewNr=radonGetNV(radtra);
  sdist=radonGetSD(radtra); 
  half=radonGetHI(radtra);
  center=radonGetCB(radtra);
 
  // Array for storing values f_x.
  X=(float*)calloc(imgDim+1,sizeof(float));  
  // Array for storing values f_y.
  Y=(float*)calloc(imgDim+1,sizeof(float));  
  if(X==NULL || Y==NULL){
    return -2;
  }
  Xptr=X;
  Yptr=Y;
 
  imgptr=imgdata;
  scnptr=scndata;

  // Find pixel coordinates of the contributing pixels for tubes of response
  // belonging to first 1/4th of angles. From these pixel coordinates
  // others can be calculated via symmetries in projection space.
  // N.B. The line of response: a*x+b*y+c
  // => solve y: y = (s - x*cos(theta))/sin(theta)
  //    solve x: x = (s - y*sin(theta))/cos(theta)
 
  for(view=set; view<=viewNr/4; view=view+setNr){
    // Choose the type of the line of response according to view number.
 
    // view=0 -> sin(theta)=0
    if(view==0){
 
      // Choose column(s) according to sample distance for angles 0 and pi/2.
      col1 = 0;
      col2 = 0;
      for(bin = 0; bin < half; bin++){
 
  col1 = col2;
  col2 = floor((float)(bin + 1)*sdist); 
 
  // Determine factor epsilon.
  if(col1 == col2){ 
    eps1 = sdist;
    eps2 = 0;
    eps3 = 0;
  }
  if((col2-col1) == 1){
    eps1 = (float)(col2 - (bin)*sdist);
    eps2 = 0;
    eps3 = (float)((bin+1)*sdist - col2);
 
  } 
  // If sdist > pixel size!
  if((col2-col1) > 1){
    eps1 = (float)(col1 + 1 - (bin)*sdist);
    eps2 = 1; // middle pixel.
    eps3 = (float)((bin+1)*sdist - col2);
 
  }
 
  /* Iterate through the entries in the image matrix.
     Calculate raysums for two LORs in the same distance from origin 
           (do halfturn). */
  for(row=0; row<imgDim; row++) {
 
    if(!eps3){   
      scnptr[bin] += eps1 * imgptr[row*imgDim + col1];
      if(bin != center)
        scnptr[binNr-bin-1] += 
                eps1 * imgptr[row*imgDim + (imgDim - 1 - col1)];
 
      scnptr[binNr*(viewNr/2) + bin] +=
              eps1 * imgptr[(imgDim - 1 - col1)*imgDim + row];
      if(bin != center)
        scnptr[binNr*(viewNr/2) + (binNr-bin-1)] += 
                eps1 * imgptr[col1*imgDim + row]; 
    }
 
    if(eps3 && !eps2){
      scnptr[bin] += eps1 * imgptr[row*imgDim + col1] + 
                           eps3 * imgptr[row*imgDim + col2];
      if(bin != center)
        scnptr[binNr-bin-1] += 
                eps1 * imgptr[row*imgDim + (imgDim - 1 - col1)] + 
                eps3*imgptr[row*imgDim+(imgDim-1-col2)];
 
      scnptr[binNr*(viewNr/2) + bin] +=
              eps1*imgptr[(imgDim-1-col1)*imgDim+row] + 
              eps3*imgptr[(imgDim-1-col2)*imgDim+row];
      if(bin != center)
        scnptr[binNr*(viewNr/2) + (binNr-bin-1)] +=
                eps1 * imgptr[col1*imgDim + row] + 
                eps3 * imgptr[col2*imgDim + row] ;
    }
 
    if(eps3 && eps2) {
      for(col = col1; col<=col2; col++) {
        if(col == col1){
           scnptr[bin] += eps1 * imgptr[row*imgDim + col1];
     if(bin != center)
       scnptr[binNr-bin-1] += eps1 * imgptr[row*imgDim + 
                    (imgDim - 1 - col1)];
     scnptr[binNr*(viewNr/2) + bin] += 
                   eps1 * imgptr[(imgDim - 1 - col1)*imgDim + row];
     if(bin != center)
       scnptr[binNr*(viewNr/2) + (binNr-bin-1)] += 
                     eps1 * imgptr[col1*imgDim + row];
        }
        if(col == col2) {
    scnptr[bin] += eps3 * imgptr[row*imgDim + col2];
    if(bin != center)
      scnptr[binNr-bin-1] += eps3 * imgptr[row*imgDim + 
                   (imgDim - 1 - col2)];
    scnptr[binNr*(viewNr/2) + bin] += 
                  eps3 * imgptr[(imgDim - 1 - col2)*imgDim + row];
    if(bin != center)
      scnptr[binNr*(viewNr/2) + (binNr-bin-1)] += 
                    eps3 * imgptr[col2*imgDim + row];
        }
        if(col != col1 && col != col2) {
    scnptr[bin] += eps2 * imgptr[row*imgDim + col];
    if(bin != center)
      scnptr[binNr-bin-1] += 
                    eps2 * imgptr[row*imgDim + (imgDim - 1 - col)];
    scnptr[binNr*(viewNr/2) + bin] += 
                  eps2 * imgptr[(imgDim - 1 - col)*imgDim + row];
    if(bin != center)
      scnptr[binNr*(viewNr/2) + (binNr-bin-1)] += 
                    eps2 * imgptr[col*imgDim + row];
        }
      }
    }
  }
      }
      // End of view==0 (handles angles 0 and pi/2).
    } else {
      
      // Set sine and cosine for this angle.
      sinus = (double)radonGetSin(radtra,view);
      cosinus = (double)radonGetSin(radtra,viewNr/2 + view);
      tanus = sinus/cosinus;
      
      // Set shift from origin for the first line of response (-n/2,theta).
      // NOTE that image matrix is on cartesian coordinate system where origin
      // is in the middle and shift is in pixels.
      shift_y = -(imgDim/2)/sinus;
      shift_x = -(imgDim/2)/cosinus;
      
      // Evaluate the function of the first LOR in integer points [-n/2,n/2].
      // NOTE that image matrix is on cartesian coordinate system where origin
      // is in the middle.
      z=-imgDim/2;
      for(col=0; col<imgDim+1; col++){
  Yptr[col]=(float)(-z/tanus + shift_y);
    Xptr[col]=(float)(-z*tanus + shift_x);
  z++;
      }
      
      // Set shift from the first TOR.
      shift_y = (double)(sdist/sinus);
      shift_x = (double)(sdist/cosinus);      
      
      // Iterate through half the bins in this view,
      // and determine coordinates of pixels contributing to this TOR.
      // NOTE that shift is added according to 'bin' in every loop.
      // Others are determined via symmetry in projection space.
 
      for(bin=0; bin<half; bin++){
 
  // Limit (x-)indices for fast search.
  // Note that indices are non-negative integers.
 
  x_left = floor((float)(Xptr[imgDim] + bin*shift_x + imgDim/2));
  if(x_left < 0) x_left = 0;
 
  x_right = floor((float)(Xptr[0] + bin*shift_x + imgDim/2));
  if(x_right <= 0) x_right = 1;
  if(x_right > imgDim) x_right = imgDim - 1; 
 
  /* Iterate through the values in vector Y, in integer points 
           [x_left,x_rigth]. */
  for(z=x_left; z <= x_right; z++) {
       
    xp = z; //positive x-coordinate
    xn = imgDim - 1 - xp; //negative x-coordinate
 
    // Look y from left.  yp=positive y-coordinate
    yp = imgDim/2 - floor(Yptr[xp + 1] + bin*shift_y) - 1;
    yn = imgDim - 1 - yp;
 
    // If the y-value for this x (z) is inside the image grid.
    if(yp < imgDim && yp >= 0 && yn < imgDim && yn >= 0 && xp >= 0 && 
             xp < imgDim && xn < imgDim && xn >= 0)
          {
 
      /* NOTE that pixels found in this part are always hit from right 
               side. Compute a := |AF| and b := |FB|. */
      a = (float)(xp + 1 - (imgDim/2 + Xptr[imgDim - yp] + bin*shift_x));
      b = (float)(floor(Yptr[xp + 1] + bin*shift_y) + 1 - (Yptr[xp + 1] + 
                 bin*shift_y));
 
      // Calculate the area of lower triangle.
      A = a*b/2;
      // c := |FC|
      c = (float)(xp + 1 - (imgDim/2 + Xptr[imgDim - yp] + 
                (bin+1)*shift_x)); 
      if(c > 0){
        // d := |FD|
        d = (float)(floor(Yptr[xp + 1] + (bin+1)*shift_y) + 1 - 
                  (Yptr[xp + 1] + (bin+1)*shift_y));
        // Subtract the area of upper triangle.
        A = A - c*d/2;
      }
 
      eps1 = A;
      if( (eps1 < 0 || eps1 > 1) && RADON_VERBOSE){
        printf("RADON: Error in factor: eps1=%.5f \n",eps1);
        errors++;
      }
 
      // Case: theta.
      // Add img(x,y)*k to the raysum of TOR (view,bin)
      scnptr[view*binNr + bin] += eps1 * imgptr[yp*imgDim + xp];
      // Add img(-x,-y)*k to the raysum of LOR (view,binNr-bin)
      if(bin != center)
        scnptr[view*binNr + binNr - 1 - bin] += 
                eps1 * imgptr[yn*imgDim + xn];
        
      if(view != viewNr/4){
        // Mirror the original LOR on y-axis, i.e. x->-x
        // Case: pi-theta.
        // Add img(-x,y)*k to the raysum of LOR (viewNr-view,bin)
        scnptr[(viewNr - view)*binNr + bin] += 
                eps1 * imgptr[yp*imgDim + xn];
        // Add img(x,-y)*k to the raysum of LOR (viewNr-view,binNr-bin)
        if(bin != center)
    scnptr[(viewNr-view)*binNr + binNr - 1 - bin] +=
                  eps1 * imgptr[yn*imgDim + xp];
        
        // Mirror the LOR on line x=y, i.e. x->y.
        // Case: pi/2-theta
        // Add img(-y,-x)*k to the raysum of LOR (viewNr/2-view,bin)
        scnptr[(viewNr/2 - view)*binNr + bin] += 
                eps1 * imgptr[xn*imgDim + yn];
        // Add img(y,x)*k to the raysum of LOR (viewNr/2-view,binNr-bin)
        if(bin != center)
    scnptr[(viewNr/2 - view)*binNr + binNr - 1 - bin] += 
                  eps1 * imgptr[xp*imgDim + yp];
      }
 
      // Mirror the LOR on line x=y, and on y-axis i.e x->y and y->-y.
      // Case: pi/2+theta
      // Add img(y,-x)*k to the raysum of LOR (viewNr/2+view,bin)
      scnptr[(viewNr/2 + view)*binNr + bin] += 
              eps1 * imgptr[xn*imgDim + yp];
      // Add img(y,-x)*k to the raysum of LOR (viewNr/2+view,binNr-bin)
      if(bin != center)
        scnptr[(viewNr/2 + view)*binNr + binNr - 1 - bin] +=
                eps1 * imgptr[xp*imgDim + yn];
    }
  }
  
  // Limit (y-)indices for fast search.
  // Note that indices are non-negative integers.
  y_bottom = floor((float)(Yptr[imgDim] + bin*shift_y + imgDim/2));
  if(y_bottom < 0) y_bottom = 0;
  if(y_bottom > imgDim) y_bottom = 0; 
 
  y_top = floor((float)(Yptr[0] + bin*shift_y + imgDim/2));
  if(y_top > imgDim) y_top = imgDim;
  if(y_top <= 0) y_top = 1;
 
  /* Iterate through the values in vector X, in integer points 
           [y_bottom,y_top]. */
  for(z=y_top; z >= y_bottom; z--) {
   
    // Look y from this location.
    xp = floor(Xptr[z] + bin*shift_x) + imgDim/2;
    xn = imgDim - 1 - xp;
 
    yp = imgDim - z - 1;
    yn = imgDim - yp - 1;
 
    // If the x-value for this y (z) is inside the image grid.
    if(yp < imgDim && yp >= 0 && yn < imgDim && yn >= 0 && xp >= 0 && 
             xp < imgDim && xn < imgDim && xn >= 0)
          {
      eps1=eps2=eps3=0;
      dx = (float)(Xptr[z] + bin*shift_x + imgDim/2 - xp);
      dy = (float)(Yptr[xp] + bin*shift_y - z + imgDim/2); 
      if(dy < 1){ // Cases 3,4,5 and 6.
        // 1. PART
        // a := |HA|
        a = dy;
        // b := |HB| (b < 1 for angles in [0,pi/4))
        b = dx;
        // h := height of rectangle R.
        h = a + shift_y;
        if(h > 1){ // Cases 3,4 and 5.
    h = 1;
    g = b + shift_x;
    if(g > 1){ // Cases 3 and 4.
      g = 1;
      xp2 =floor(Xptr[z+1] + (bin+1)*shift_x) + imgDim/2;
      if(xp == xp2){ // Case 4.
        // c := |FC|
        c = (float)(xp + 1 - (imgDim/2 + Xptr[z+1] + 
                        (bin+1)*shift_x));        
        // d := |FD| 
        d = (float)(floor(Yptr[xp + 1] + (bin+1)*shift_y) + 1 - 
                        (Yptr[xp + 1] + (bin+1)*shift_y));
        eps1 = 1 - (a*b + c*d)/2;
        eps2 = 0;
        // Lower triangle on the next pixel (x+1,y).
        eps3 = (1 - d)*(b + shift_x - 1)/2; 
        if((eps1<0 || eps1>1 || eps2<0 || eps2>1 || eps3<0 || eps3>1) 
                       && RADON_VERBOSE) printf(
                      "4: (%i,%i) eps1=%f | (%i,%i) eps2=%f | (%i,%i) eps3=%f\n",
                      xp,yp,eps1,xp+1,yp-1,eps2,xp+1,yp,eps3);
      } else { // Case 3.
        // c=d=0.
        eps1 = 1 - a*b/2;
 
        // Calculate area on pixel in place (xp+1,yp-1).
        dy = (float)(Yptr[xp+1] + (bin+1)*shift_y - (z + 1) + 
                         imgDim/2);
        if(dy < 1){ // Case 11.
          // c := |HC|
          c = dy;
          // d := |HD|
          d = (float)(xp + 2 - (imgDim/2 + Xptr[z+1] + 
                          (bin+1)*shift_x));
          eps2 = c*d/2;
        } else { // Cases 9 and 10.
          dx = (float)(xp + 2 - (imgDim/2 + Xptr[z+1] + 
                           (bin+1)*shift_x));
          if(dx < 1) { // Case 10.
      // g := length of rectangle Q.
      g = dx;
      // c := |CD| (on x-axis).
      c = dx - (float)(xp + 2 - (imgDim/2 + Xptr[z+2] + 
                           (bin+1)*shift_x));
      // Rectangle Q - triangle c*h (h = 1).
      eps2 = g - c/2;
          } else { // Case 9.
      // c := |FC|
      c = (float)(xp + 2 - (imgDim/2 + Xptr[z+2] + 
                            (bin+1)*shift_x));        
      // d := |FD| 
      d = (float)(floor(Yptr[xp + 2] + (bin+1)*shift_y) + 2 - 
                            (Yptr[xp + 2] + (bin+1)*shift_y));
      // Rectangle Q - triangle CFD.
      eps2 = 1 - c*d/2;
          }
        }
 
        // Calculate area on pixel in place (xp+1,yp).
        dx = (float)(xp + 2 - (imgDim/2 + Xptr[z] + (bin+1)*shift_x));
        if(dx < 1){ // Case 10.
          // g := length of rectangle Q.
          g = dx;
          // c := |CD| (on x-axis).
          c = dx - (float)(xp + 2 - (imgDim/2 + Xptr[z+1] + 
                          (bin+1)*shift_x));
          // Rectangle Q - triangle c*h (h = 1).
          eps3 = g - c/2;
        } else { // Case 9.
          // c := |FC|
          c = (float)(xp + 2 - (imgDim/2 + Xptr[z+1] + 
                           (bin+1)*shift_x));       
          // d := |FD| 
          d = (float)(floor(Yptr[xp + 2] + (bin+1)*shift_y) + 1 - 
                           (Yptr[xp + 2] + (bin+1)*shift_y));
          // Rectangle Q - triangle CFD.
          eps3 = 1 - c*d/2;
        }
        if((eps1<0 || eps1>1 || eps2<0 || eps2>1 || eps3<0 || eps3>1) 
                       && RADON_VERBOSE) printf(
                      "3/v%i: (%i,%i) eps1=%f | (%i,%i) eps2=%f | (%i,%i) eps3=%f\n",
                      view,xp,yp,eps1,xp+1,yp-1,eps2,xp+1,yp,eps3);
      }
    } else { // Case 5. (g < 1)
      // c := |DC|.
      c = g - (imgDim/2 + Xptr[z+1] + (bin+1)*shift_x - xp);
      // d := heigth
      d = 1;
      eps1 = g*h - (a*b + c*d)/2;
      eps2 = 0;
      eps3 = 0;
      if((eps1<0 || eps1>1 || eps2<0 || eps2>1 || eps3<0 || eps3>1) 
                     && RADON_VERBOSE) printf(
                   "5: (%i,%i) eps1=%f | (%i,%i) eps2=%f | (%i,%i) eps3=%f\n",
                   xp,yp,eps1,xp+1,yp-1,eps2,xp+1,yp,eps3);
    }
        } else{ // Case 6 (h <= 1). 
    // g := legth of rectangle R
    g = b + shift_x;
    if(g > 1) // Should always be < 1 for angles in [0,pi/4)
      g = 1; 
    eps1 = (g*h - a*b)/2;
    eps2 = 0;
    eps3 = 0;
    if((eps1<0 || eps1>1 || eps2<0 || eps2>1 || eps3<0 || eps3>1) &&
                   RADON_VERBOSE) printf(
                  "6: (%i,%i) eps1=%f | (%i,%i) eps2=%f | (%i,%i) eps3=%f\n",
                  xp,yp,eps1,xp+1,yp-1,eps2,xp+1,yp,eps3);
        }
      } else { // 2. PART
        // Cases 2,7 and 8. (dy >= 1).
        // a := |HA|
        a = 1;
        // b := |HB| (>=1)
        b = dx - (float)(Xptr[z+1] + bin*shift_x + imgDim/2 - xp);
        // h := heigth of rectangle R
        h = 1;
        // g := length of rectangle R
        g = dx + shift_x;
        if(g > 1){ // Cases 2 and 8.
    g = 1 - (float)(Xptr[z+1] + bin*shift_x + imgDim/2 - xp);
                // positive x-coordinate (bin+1)
    xp2 = floor(Xptr[z+1] + (bin+1)*shift_x) + imgDim/2; 
    if(xp == xp2){ // Case 8.
        // c := |FC|
        c = (float)(xp + 1 - (imgDim/2 + Xptr[z+1] + 
                         (bin+1)*shift_x));       
        // d := |FD| 
        d = (float)(floor(Yptr[xp + 1] + (bin+1)*shift_y) + 1 - 
                         (Yptr[xp + 1] + (bin+1)*shift_y));
        eps1 = g*h - (a*b + c*d)/2;
        eps2 = 0;
        // Lower triangle on the next pixel (x+1,y).
        eps3 = (1 - d)*((Xptr[z] + (bin+1)*shift_x + imgDim/2) - 
                            (xp+1))/2; 
        if((eps1<0 || eps1>1 || eps2<0 || eps2>1 || eps3<0 || eps3>1) 
                       && RADON_VERBOSE) printf(
                     "8: (%i,%i) eps1=%f | (%i,%i) eps2=%f | (%i,%i) eps3=%f\n",
                     xp,yp,eps1,xp+1,yp-1,eps2,xp+1,yp,eps3);
      } else{ // Case 2.
        // c=d=0.
        eps1 = g*h - a*b/2;
        /* Pixel in place (xp+1,yp-1) should have been found in 
                       the previous step. */
        eps2 = 0;
        // Calculate area on pixel in place (xp+1,yp).
        dx = (float)((imgDim/2 + Xptr[z] + (bin+1)*shift_x) - (xp+1));
        if(dx < 1){ // Case 10 (trapezium).
          // g := bottom of trapezium Q.
          g = dx;
          // c := top of trapezium Q.
          c = (float)((imgDim/2 + Xptr[z+1] + (bin+1)*shift_x) - 
                                  (xp+1));
          // Area of trapezium Q. Heigth = 1.
          eps3 = (g + c)/2;
          if((eps1<0 || eps1>1 || eps2<0 || eps2>1 || eps3<0 || 
                          eps3>1) && RADON_VERBOSE) printf(
                      "2/10: (%i,%i) eps1=%f | (%i,%i) eps2=%f | (%i,%i) eps3=%f\n",
                      xp,yp,eps1,xp+1,yp-1,eps2,xp+1,yp,eps3);
        } else { // Case 9.
          // c := |FC|
          c = (float)(xp + 2 - (imgDim/2 + Xptr[z+1] + 
                           (bin+1)*shift_x));       
          // d := |FD| 
          d = (float)(floor(Yptr[xp + 2] + (bin+1)*shift_y) + 1 - 
                           (Yptr[xp + 2] + (bin+1)*shift_y));
          // Rectangle Q - triangle CFD.
          eps3 = 1 - c*d/2;
          if((eps1<0 || eps1>1 || eps2<0 || eps2>1 || eps3<0 || 
                          eps3>1) && RADON_VERBOSE) printf(
                      "2/9: (%i,%i) eps1=%f | (%i,%i) eps2=%f | (%i,%i) eps3=%f\n",
                      xp,yp,eps1,xp+1,yp-1,eps2,xp+1,yp,eps3);
        }
      }
        } else { // Case 7. (g < = 1)
    // Area of the parallelogram R.
    eps1 = sdist/cosinus;
    eps2 = 0;
    eps3 = 0;
    if((eps1<0 || eps1>1 || eps2<0 || eps2>1 || eps3<0 || eps3>1) 
                  && RADON_VERBOSE) printf(
                "7: (%i,%i) eps1=%f | (%i,%i) eps2=%f | (%i,%i) eps3=%f\n",
                xp,yp,eps1,xp+1,yp-1,eps2,xp+1,yp,eps3);
        }
      }
      if(!eps2 && !eps3){ // Cases 5,6 and 7.
        // Case: theta.
        // Add img(x,y)*k to the raysum of LOR (view,bin)
        scnptr[view*binNr + bin] += eps1 * imgptr[yp*imgDim + xp];
        // Add img(-x,-y)*k to the raysum of LOR (view,binNr-bin)
        if(bin != center)
    scnptr[view*binNr + binNr - 1 - bin] += 
                  eps1 * imgptr[yn*imgDim + xn];
        if(view != viewNr/4){     
          // Mirror the LOR on y-axis, i.e. x->-x
          // Case: pi-theta.
          // Add img(-x,y)*k to the raysum of LOR (viewNr-view,bin)
          scnptr[(viewNr - view)*binNr + bin] += 
                  eps1 * imgptr[yp*imgDim + xn];
          // Add img(x,-y)*k to the raysum of LOR (viewNr-view,binNr-bin)
          if(bin != center)
      scnptr[(viewNr-view)*binNr + binNr - 1 - bin] += 
                    eps1 * imgptr[yn*imgDim + xp];
        
          // Mirror the LOR on line y=x, i.e. y->x.
          // Case: pi/2 - theta.
          // Add img(y,x)*k to the raysum of LOR (viewNr/2-view,bin).
          scnptr[(viewNr/2 - view)*binNr + bin] += 
                  eps1 * imgptr[xn*imgDim + yn];
          // Add img(-y,-x)*k to the raysum of LOR (viewNr/2-view,binNr-bin)
          if(bin != center)
      scnptr[(viewNr/2 - view)*binNr + binNr - 1 - bin] += 
                    eps1 * imgptr[xp*imgDim + yp];
        }
 
        // Mirror the LOR on line x=y, and on y-axis i.e x->y and y->-y.
        // Case: pi/2 + theta.
        // Add img(y,-x)*k to the raysum of LOR (viewNr/2+view,bin)
        scnptr[(viewNr/2 + view)*binNr + bin] += 
                eps1 * imgptr[xn*imgDim + yp];
        // Add img(-y,x)*k to the raysum of LOR (viewNr-view,binNr-bin)
        if(bin != center)
    scnptr[(viewNr/2 + view)*binNr + binNr - 1 - bin] += 
                  eps1 * imgptr[xp*imgDim + yn];
      } else {
        if(!eps2) { // <=> eps3 != 0 & eps2 = 0 <=> Cases 3,4 and 8.
    if(xp + 1 < imgDim && xn - 1 >= 0){
      // Case: theta.
      // Add img(x,y)*k to the raysum of LOR (view,bin)
      scnptr[view*binNr + bin] +=
                    eps1 * imgptr[yp*imgDim + xp] + 
                    eps3 * imgptr[yp*imgDim + xp+1];
      // Add img(-x,-y)*k to the raysum of LOR (view,binNr-bin)
      if(bin != center)   
        scnptr[view*binNr + binNr - 1 - bin] += 
                      eps1 * imgptr[yn*imgDim + xn] + 
                      eps3 * imgptr[yn*imgDim + xn-1];
      if(view != viewNr/4) {
        // Mirror the LOR on y-axis, i.e. x->-x
        // Case: pi-theta.
        // Add img(-x,y)*k to the raysum of LOR (viewNr-view,bin)
        scnptr[(viewNr - view)*binNr + bin] += 
                      eps1 * imgptr[yp*imgDim + xn] + 
                      eps3 * imgptr[yp*imgDim + xn-1];
        // Add img(x,-y)*k to the raysum of LOR (viewNr-view,binNr-bin)
        if(bin != center)
          scnptr[(viewNr-view)*binNr + binNr - 1 - bin] += 
                        eps1 * imgptr[yn*imgDim + xp] + 
                        eps3 * imgptr[yn*imgDim + xp+1];
      
        // Mirror the LOR on line y=x, i.e. y->x.
        // Case: pi/2 - theta.
        // Add img(y,x)*k to the raysum of LOR (viewNr/2-view,bin).
        scnptr[(viewNr/2 - view)*binNr + bin] += 
                      eps1 * imgptr[xn*imgDim + yn] +
                      eps3 * imgptr[(xn-1)*imgDim + yn];
        /* Add img(-y,-x)*k to the raysum of LOR 
                       (viewNr/2-view,binNr-bin) */
        if(bin != center)
          scnptr[(viewNr/2 - view)*binNr + binNr - 1 - bin] += 
                        eps1 * imgptr[xp*imgDim + yp] +
                        eps3*imgptr[(xp+1)*imgDim+yp];
      }
      // Mirror the LOR on line x=y, and on y-axis i.e x->y and y->-y.
      // Case: pi/2 + theta.
      // Add img(y,-x)*k to the raysum of LOR (viewNr/2+view,bin)
      scnptr[(viewNr/2 + view)*binNr + bin] += 
                    eps1 * imgptr[xn*imgDim + yp] +
                    eps3 * imgptr[(xn-1)*imgDim + yp];
      // Add img(-y,x)*k to the raysum of LOR (viewNr-view,binNr-bin)
      if(bin != center)
        scnptr[(viewNr/2 + view)*binNr + binNr - 1 - bin] +=
                      eps1 * imgptr[xp*imgDim + yn] +
                      eps3*imgptr[(xp+1)*imgDim+yn];
    }
        } else { // <=> eps2!=0 && eps3!=0 <=> Case 3.
    if(xp + 1 < imgDim && xn - 1 >= 0 && yp-1 >= 0 && yn+1 < imgDim) {
      // Case: theta.
      // Add img(x,y)*k to the raysum of LOR (view,bin)
      scnptr[view*binNr + bin] += 
                    eps1 * imgptr[yp*imgDim + xp] + 
                    eps3 * imgptr[yp*imgDim + xp+1] +
        eps2 * imgptr[(yp-1)*imgDim + xp+1];
      // Add img(-x,-y)*k to the raysum of LOR (view,binNr-bin)
      if(bin != center)
        scnptr[view*binNr + binNr - 1 - bin] += 
                      eps1 * imgptr[yn*imgDim + xn] + 
                      eps3 * imgptr[yn*imgDim + xn-1] + 
          eps2 * imgptr[(yn+1)*imgDim + xn-1];
 
      if(view != viewNr/4){     
        // Mirror the LOR on y-axis, i.e. x->-x
        // Case: pi-theta.
        // Add img(-x,y)*k to the raysum of LOR (viewNr-view,bin)
        scnptr[(viewNr - view)*binNr + bin] += 
                      eps1 * imgptr[yp*imgDim + xn] +
                      eps3 * imgptr[yp*imgDim + xn-1] +
          eps2 * imgptr[(yp-1)*imgDim + xn-1]; 
        // Add img(x,-y)*k to the raysum of LOR (viewNr-view,binNr-bin)
        if(bin != center)
          scnptr[(viewNr-view)*binNr + binNr - 1 - bin] +=
                        eps1 * imgptr[yn*imgDim + xp] +
                        eps3 * imgptr[yn*imgDim + xp+1] +
            eps2 * imgptr[(yn+1)*imgDim + xp+1];
      
        // Mirror the LOR on line y=x, i.e. y->x.
        // Case: pi/2 - theta.
        // Add img(y,x)*k to the raysum of LOR (viewNr/2-view,bin)      
        scnptr[(viewNr/2 - view)*binNr + bin] +=
                      eps1 * imgptr[xn*imgDim + yn] +
                      eps3 * imgptr[xn*imgDim + yn-1] +
          eps2 * imgptr[(xn-1)*imgDim + (yn+1)];
        // Add img(-y,-x)*k to the raysum of LOR (viewNr/2-view,binNr-bin)
        if(bin != center)
          scnptr[(viewNr/2 - view)*binNr + binNr - 1 - bin] +=
                        eps1 * imgptr[xp*imgDim + yp] +
                        eps3*imgptr[(xp+1)*imgDim+yp] +
            eps2*imgptr[(xp+1)*imgDim+(yp-1)];
      }
 
      // Mirror the LOR on line x=y, and on y-axis i.e x->y and y->-y.
      // Case: pi/2 + theta.
      // Add img(y,-x)*k to the raysum of LOR (viewNr/2+view,bin)
      scnptr[(viewNr/2 + view)*binNr + bin] += 
                    eps1 * imgptr[xn*imgDim + yp] +
                    eps3 * imgptr[(xn-1)*imgDim + yp] +
        eps2 * imgptr[(xn-1)*imgDim + (yp-1)];
      // Add img(-y,x)*k to the raysum of LOR (viewNr-view,binNr-bin)
      if(bin != center)
        scnptr[(viewNr/2 + view)*binNr + binNr - 1 - bin] +=
              eps1 * imgptr[xp*imgDim + yn] +
              eps3*imgptr[(xp+1)*imgDim+yn] + 
          eps2*imgptr[(xp+1)*imgDim+yn+1];
    }
        }
      } 
    }
  }
      }// END of X loop
    }// End of view>0.
  }// END OF VIEW LOOP
  free(X);
  free(Y);
  if(RADON_VERBOSE) {
    printf("RADON: radonFwdTransformEA() finished with %d errors.\n", errors);
    fflush(stdout);
  }
  return 0;
} // END OF EXACT AREA FORWARD TRANSFORM
/*****************************************************************************/
int radonFwdTransformSA(
  RADON *radtra, int set, int setNr, float *imgdata, float *scndata
) {
  float *imgptr, *scnptr, *imgorigin;
  float sinus, cosinus; 
  float t;   // distance between projection ray and origo.
  int   mode, imgDim, binNr, viewNr, halfImg, centerBin, view;
  int x, y, xright, ybottom;
  float fract, tpow2;
 
  // Retrieve the data from given radon transform object.
  mode=radonGetMO(radtra); // Should be 3 or 4.
  imgDim=radonGetID(radtra);
  binNr=radonGetNB(radtra); // Should be equal to imgDim.
  viewNr=radonGetNV(radtra);
  // Set the center ray and the square of it.
  centerBin=binNr/2;
  tpow2 = centerBin*centerBin;
 
  if(imgDim != binNr) return -1;
 
  // Set half of the image dimension.
  halfImg = imgDim/2;
  
  // Transform one angle at a time.
  for(view=set; view<viewNr; view+=setNr) {
 
    imgorigin = imgdata + imgDim*(halfImg - 1) + halfImg;
 
    sinus = radonGetSin(radtra,view);
    cosinus = radonGetSin(radtra,viewNr/2 + view);
 
    y = halfImg - 2;
    
    if((float)y > centerBin){
      y = (int)(centerBin);
      ybottom = -y;
    } else{
      ybottom = -halfImg + 1;
    }
 
    for(; y >= ybottom; y--) {
      xright = (int)sqrt(/*fabs(*/tpow2 - ((float)y+0.5) * 
                         ((float)y+0.5)/*)*/) + 1;
      if(xright >= halfImg){
  xright = halfImg - 1;
  x = -halfImg;
      } else
  x = -xright;
      
      //      t = centerBin - (float)y * sinus + ((float)(x + 1)) * cosinus;
      t = centerBin + (float)y * sinus + ((float)(x + 1)) * cosinus;
      
      imgptr = imgorigin - y*imgDim + x;
      scnptr = scndata + view*binNr;
      //      for(; x <= xright; x++, t += cosinus){
      for(; x <= xright; x++, t += cosinus){
  if(mode == 3){ // If the linear interpolation is to be utilised.
    fract = t - (float)(int)t;
    *(scnptr+(int)t) += *imgptr * (1.0 - fract);
    *(scnptr+(int)t + 1) += *imgptr++ * fract;
  }
  else // If the nearest neighbour interpolation is to be utilised.
    *(scnptr+(int)(t + 0.5)) += *imgptr++;
      }
    }
  }
 
  return 0;
}// END OF FORWARD TRANSFORM USING THE IMAGE ROTATION APPROACH
/*****************************************************************************/
int radonFwdTransformPRM(
  PRMAT *mat, int set, int setNr, float *imgdata, float *scndata
) {
  float *imgptr, *scnptr; // pointers for the image and sinogram data
  unsigned int imgDim=128, binNr=256, viewNr=192, view, bin, row, col;
  unsigned int half, center;
  int xp, xn, yp, yn;
  float fact;
 
  // Retrieve the data from given projection matrix.
  imgDim=prmatGetID(mat);
  binNr=prmatGetNB(mat); // Should be equal to imgDim.
  viewNr=prmatGetNV(mat);
 
  // Calculate and set the center bin.
  if((binNr%2) != 0){
    half = (binNr - 1)/2 + 1;
    center = half - 1;
  } else{
    half = binNr/2;
    // In the case binNr is even there is no center bin.
    center = -1;
  }
 
  imgptr = imgdata;
  scnptr = scndata;
 
  // Draw sinogram according to projection matrix.
  for(view=set; view<=viewNr/4; view=view+setNr){
    for(bin=0; bin<half; bin++){
      row = view*half + bin;
      for(col=0; col<prmatGetPixels(mat,row); col++){
 
  fact = prmatGetFactor(mat,row,col);
  xp = prmatGetXCoord(mat,row,col);
  xn = imgDim - xp - 1;
 
  yp = prmatGetYCoord(mat,row,col);
  yn = imgDim - yp - 1;
 
  // Case: theta.
  // Add img(x,y)*k to the raysum of LOR (view,bin)
  scnptr[view*binNr + bin] += fact * imgptr[yp*imgDim + xp];
  // Add img(-x,-y)*k to the raysum of LOR (view,binNr-bin)
  if(bin != center)
    scnptr[view*binNr + binNr - 1 - bin] += fact * imgptr[yn*imgDim + xn];
 
  if(view != 0 && view != viewNr/4){              
    // Mirror the LOR on y-axis, i.e. x->-x
    // Case: pi-theta.
    // Add img(-x,y)*k to the raysum of LOR (viewNr-view,bin)
    scnptr[(viewNr - view)*binNr + bin] += fact * imgptr[yp*imgDim + xn];
    // Add img(x,-y)*k to the raysum of LOR (viewNr-view,binNr-bin)
    if(bin != center)
      scnptr[(viewNr-view)*binNr + binNr - 1 - bin] += 
              fact * imgptr[yn*imgDim + xp];
        
    // Mirror the LOR on line y=x, i.e. y->x.
    // Case: pi/2 - theta.
    // Add img(y,x)*k to the raysum of LOR (viewNr/2-view,bin)      
    scnptr[(viewNr/2 - view)*binNr + bin] += fact * imgptr[xn*imgDim + yn];
    // Add img(-y,-x)*k to the raysum of LOR (viewNr/2-view,binNr-bin)
    if(bin != center)
      scnptr[(viewNr/2 - view)*binNr + binNr - 1 - bin] += 
              fact * imgptr[xp*imgDim + yp];
  }
 
  // Mirror the LOR on line x=y, and on y-axis i.e x->y and y->-y.
  // Case: pi/2 + theta.
  // Add img(y,-x)*k to the raysum of LOR (viewNr/2+view,bin)
  scnptr[(viewNr/2 + view)*binNr + bin] += fact * imgptr[xn*imgDim + yp];
  // Add img(-y,x)*k to the raysum of LOR (viewNr-view,binNr-bin)
  if(bin != center)
    scnptr[(viewNr/2 + view)*binNr + binNr - 1 - bin] +=
            fact * imgptr[xp*imgDim + yn];    
      }
    }
  }
 
  return 0;
 
}// END OF FORWARD TRANSFORM WITH A PROJECTION MATRIX
/*****************************************************************************/
// BACK TRANSFORM METHODS
/*****************************************************************************/
int radonBackTransform(
  RADON *radtra, int set, int setNr, float *scndata, float *imgdata
) {
  float *imgptr, *scnptr; // pointers for the image and sinogram data
  float *Y, *X, *Xptr, *Yptr; // pointers for the values of LOR in integer points
  double sinus, cosinus, tanus; // sin, cosine and tangent
  double shift_x, shift_y; // shift for LOR
  int xp, xn, yp, yn, z; //integer points
  int x_left, x_right, y_top, y_bottom; // limits for fast search
  int col, row, view, bin; //counters
  int binNr, viewNr, imgDim;
  int half, center = -1, mode = 0;
  float sdist;
  float dx, dy , loi = 1;  // delta x and y and length of intersection
 
  //Check that the data for this radon transform is initialized
  if(radtra->status != RADON_STATUS_INITIALIZED) return -1;
 
  //Retrieve the data from given radon transform object
  mode=radonGetMO(radtra);
  imgDim=radonGetID(radtra);
  binNr=radonGetNB(radtra);
  viewNr=radonGetNV(radtra);
  sdist=radonGetSD(radtra);
  half=radonGetHI(radtra);
  center=radonGetCB(radtra);

  // Find pixel coordinates of the contributing pixels for lines of response
  // corresponding to first 1/4th of angles. From these pixel coordinates
  // others can be calculated via symmetries in projection space.
  // N.B. The line of response: a*x+b*y+c=cos(view)*x+sin(view)*y+k*bin=0
  // => solve y: y=-x*(cos(view)/sin(view)) - k*bin
  //    solve x: x=-y*(sin(view)/cos(view))) - k*bin
 
  X=(float*)calloc(imgDim+1,sizeof(float));  
  Y=(float*)calloc(imgDim+1,sizeof(float));  
  if(X==NULL || Y==NULL){
    return -2;
  }
  Xptr=X;
  Yptr=Y;
  imgptr=imgdata;
  scnptr=scndata;
  for(view=set; view<viewNr/4; view+=setNr){
    // Choose the type of the line of response according to view number
 
    // Length of intersection is 1.
    loi = 1;
    
    // 0. type: view=0 -> sin(view)=0
    if(view==0){
 
      // Choose the pixels according to sample distance and pixel size,
      // for angles 0 and pi/2.
      for(bin=0; bin<binNr; bin++){
  col=floor((float)(bin+.5*sdist)*sdist); 
  if(col==imgDim) col=imgDim-1;
 
  // Iterate through the entries in the image matrix.
  // Calculate raysums for LORs in the same (absolute) distance from origin.
  for(row=0; row<imgDim; row++){
    imgptr[row*imgDim + col] += loi*scnptr[bin];
    imgptr[(imgDim - 1 - col)*imgDim + row] +=
            loi*scnptr[binNr*(viewNr/2) + bin];
  }
      }
      // End of view==0 (handles angles 0 and pi/2).   
    } else {
      // Set sine and cosine for this angle.
      sinus=(double)radonGetSin(radtra,view);
      cosinus=(double)radonGetSin(radtra,viewNr/2 + view);
      tanus=sinus/cosinus;

      // Set shift from origin for the first line of response (-n/2,theta).
      // NOTE that image matrix is on cartesian coordinate system where origin
      // is in the middle and shift is in pixels.
      shift_y = -(imgDim/2 -.5*sdist)/sinus;
      shift_x = -(imgDim/2 -.5*sdist)/cosinus;
 
      // Evaluate the function of the first LOR in integer points [-n/2,n/2].
      // NOTE that image matrix is on cartesian coordinate system where origin
      // is in the middle.
      z=-imgDim/2;
      for(col=0; col<imgDim+1; col++){
  Yptr[col]=(float)(shift_y - z/tanus);
    Xptr[col]=(float)(shift_x - z*tanus);
  z++;
      }
 
      // Set shift from the first LOR.
      shift_y = (double)(sdist/sinus);
      shift_x = (double)(sdist/cosinus);
            
      // Iterate through half the bins in this view,
      // and determine coordinates of pixels contributing to this LOR.
      // NOTE that shift is added according to 'bin' in every loop.
      // Calculate also the length of intersection.
      // Others are determined via symmetry in projection space.
 
      for(bin=0; bin<half; bin++){
 
  // Limit (x-)indices for fast search.
  // Note that indices are non-negative integers.
 
  x_left = floor((float)((Xptr[imgDim] + bin*shift_x) + imgDim/2));
  if(x_left < 0) x_left = 0;
  x_right = floor((float)((Xptr[0] + bin*shift_x) + imgDim/2));
  if(x_right <= 0) x_right = 1;
  if(x_right > imgDim) x_right = imgDim - 1; 
 
  /* Iterate through the values in vector Y, in integer points 
           [x_left,x_rigth]. */
  for(z=x_left; z <= x_right; z++){
       
    xp = z; //positive x-coordinate
    xn = imgDim - 1 - xp; //negative x-coordinate
 
    // Look y from left.
    yp = imgDim/2 - floor(Yptr[xp + 1] + bin*shift_y) - 1;
    yn = imgDim - 1 - yp;
 
    // If the y-value for this x (z) is inside the image grid.
    if(yp < imgDim && yp >= 0 && yn < imgDim && yn >= 0 && xp >= 0 && 
             xp < imgDim && xn < imgDim && xn >= 0)
          {
      if(!mode){        
        loi = 1;
      } else {
        // Compute delta x and y from 'positive' coordinates. 
        dx = 1 - (float)((Xptr[imgDim - yp] + bin*shift_x) + imgDim/2 - xp);
        dy = 1 - (float)(yp - (imgDim/2 - (Yptr[xp + 1] + bin*shift_y) - 1));
        
        if(dx > 1 || dx < 0) dx = 1;
        if(dy > 1 || dy < 0) dy = 1;
        loi = sqrt(dx*dx + dy*dy);
      }
 
      // Case: theta.
      // Add img(x,y)*k to the raysum of LOR (view,bin)
      imgptr[yp*imgDim + xp] += loi*scnptr[view*binNr + bin];
      // Add img(-x,-y)*k to the raysum of LOR (view,binNr-bin)
      if(bin != center)
        imgptr[yn*imgDim + xn] += 
                loi*scnptr[view*binNr + binNr - 1 - bin];
 
      if(view != viewNr/4){       
        // Mirror the original LOR on y-axis, i.e. x->-x
        // Case: pi-theta.
        // Add img(-x,y)*k to the raysum of LOR (viewNr-view,bin)
        imgptr[yp*imgDim + xn] += loi*scnptr[(viewNr - view)*binNr + bin];
        // Add img(x,-y)*k to the raysum of LOR (viewNr-view,binNr-bin)
        if(bin != center)
    imgptr[yn*imgDim + xp] += 
                  loi*scnptr[(viewNr-view)*binNr + binNr - 1 - bin];
        
        // Mirror the LOR on line x=y, i.e. x->y.
        // Case: pi/2-theta
        // Add img(-y,-x)*k to the raysum of LOR (viewNr/2-view,bin)
        imgptr[xn*imgDim + yn] += loi*scnptr[(viewNr/2 - view)*binNr + bin];
        // Add img(y,x)*k to the raysum of LOR (viewNr/2-view,binNr-bin)
        if(bin != center)
    imgptr[xp*imgDim + yp] += 
                  loi*scnptr[(viewNr/2 - view)*binNr + binNr - 1 - bin];
      }
 
      // Mirror the LOR on line x=y, and on y-axis i.e x->y and y->-y.
      // Case: pi/2+theta
      // Add img(y,-x)*k to the raysum of LOR (viewNr/2+view,bin)
      imgptr[xn*imgDim + yp] += loi*scnptr[(viewNr/2 + view)*binNr + bin];
      // Add img(y,-x)*k to the raysum of LOR (viewNr/2+view,binNr-bin)
      if(bin != center)
        imgptr[xp*imgDim + yn] += 
                loi*scnptr[(viewNr/2 + view)*binNr + binNr - 1 - bin];
        
    }
  }
 
  // Limit (y-)indices for fast search.
  // Note that indices are non-negative integers.
  y_bottom = floor((float)(Yptr[imgDim] + bin*shift_y + imgDim/2));
  if(y_bottom < 0) y_bottom = 0;
  if(y_bottom > imgDim) y_bottom = 0; 
 
  y_top = floor((float)(Yptr[0] + bin*shift_y + imgDim/2));
  if(y_top > imgDim) y_top = imgDim;
  if(y_top <= 0) y_top = 1;
 
  /* Iterate through the values in vector X, in integer points 
           [y_bottom,y_top]. */
  for(z=y_top; z >= y_bottom; z--) {
   
    // Look y from this location.
    xp = floor(Xptr[z] + bin*shift_x) + imgDim/2;
    xn = imgDim - 1 - xp;
 
    yp = imgDim - z - 1;
    yn = imgDim - yp - 1;
 
    // If the x-value for this y (z) is inside the image grid.
    if(yp < imgDim && yp >= 0 && yn < imgDim && yn >= 0 && xp >= 0 && 
             xp < imgDim && xn < imgDim && xn >= 0)
          {
      if(!mode){
        loi = 1;
      } else{
        dx = (float)((Xptr[z] + bin*shift_x) + imgDim/2 - xp);
        dy = (float)(yp - (imgDim/2 - (Yptr[xp] + bin*shift_y) - 1));
        if(dy > 1 || dy < 0) dy = 1;
        loi = sqrt(dx*dx + dy*dy);
      }
      // Case: theta.
      // Add img(x,y)*k to the raysum of LOR (view,bin)
      imgptr[yp*imgDim + xp] += loi*scnptr[view*binNr + bin];
      // Add img(-x,-y)*k to the raysum of LOR (view,binNr-bin)
      if(bin != center)
        imgptr[yn*imgDim + xn] +=
                loi*scnptr[view*binNr + binNr - 1 - bin];
 
      if(view != viewNr/4){             
        // Mirror the LOR on y-axis, i.e. x->-x
        // Case: pi-theta.
        // Add img(-x,y)*k to the raysum of LOR (viewNr-view,bin)
        imgptr[yp*imgDim + xn] += loi*scnptr[(viewNr - view)*binNr + bin];
        // Add img(x,-y)*k to the raysum of LOR (viewNr-view,binNr-bin)
        if(bin != center)
    imgptr[yn*imgDim + xp] +=
                  loi*scnptr[(viewNr-view)*binNr + binNr - 1 - bin];
        
        // Mirror the LOR on line y=x, i.e. y->x.
        // Case: pi/2 - theta.
        // Add img(y,x)*k to the raysum of LOR (viewNr/2-view,bin)      
        imgptr[xn*imgDim + yn] += loi*scnptr[(viewNr/2 - view)*binNr + bin];
        // Add img(-y,-x)*k to the raysum of LOR (viewNr/2-view,binNr-bin)
        if(bin != center)
    imgptr[xp*imgDim + yp] +=
                  loi*scnptr[(viewNr/2 - view)*binNr + binNr - 1 - bin];
      }
 
      // Mirror the LOR on line x=y, and on y-axis i.e x->y and y->-y.
      // Case: pi/2 + theta.
      // Add img(y,-x)*k to the raysum of LOR (viewNr/2+view,bin)
      imgptr[xn*imgDim + yp] += loi*scnptr[(viewNr/2 + view)*binNr + bin];
      // Add img(-y,x)*k to the raysum of LOR (viewNr-view,binNr-bin)
      if(bin != center)
        imgptr[xp*imgDim + yn] += 
                loi*scnptr[(viewNr/2 + view)*binNr + binNr - 1 - bin];
    }
  }
      } // END of X loop
    } // End of view>0.
  } // END OF VIEW LOOP
  free(X);
  free(Y);
  return 0;
} // END OF BACK (pro -> spa) TRANSFORM WITH 0/1 OR LOI-MODEL.
/*****************************************************************************/
int radonBackTransformEA(
  RADON *radtra, int set, int setNr, float *scndata, float *imgdata
) {
  float *imgptr, *scnptr; // pointers for the image and sinogram data
  float *Y, *X, *Xptr, *Yptr; // pointers for the values of LOR in integer points
  double sinus, cosinus, tanus; // sin, cosine and tangent
  double shift_x, shift_y; // shift for LOR
  int xp, xn, yp, yn, z, xp2; //integer points
  int x_left, x_right, y_top, y_bottom; // limits for fast search
  int col, col1, col2, row, view, bin; //counters
  int binNr, viewNr, imgDim, errors=0;
  int half, center = -1;
  float sdist;
  float a=0,b=0, c=0, d=0, g=0, h=0, A;
  float dx, dy , eps1 = 0, eps2 = 0, eps3 = 0;  // delta x and y and factors.
 
  // Check that the data for this radon transform is initialized.
  if(radtra->status != RADON_STATUS_INITIALIZED) return -1;
 
  // Retrieve the data from given radon transform object.
  //mode=radonGetMO(radtra); // Should be 2.
  imgDim=radonGetID(radtra);
  binNr=radonGetNB(radtra);
  viewNr=radonGetNV(radtra);
  sdist=radonGetSD(radtra); 
  half=radonGetHI(radtra);
  center=radonGetCB(radtra);
 
  // Array for storing values f_x.
  X=(float*)calloc(imgDim+1,sizeof(float));  
  // Array for storing values f_y.
  Y=(float*)calloc(imgDim+1,sizeof(float));  
  if(X==NULL || Y==NULL){
    return -2;
  }
  Xptr=X;
  Yptr=Y;
 
  imgptr=imgdata;
  scnptr=scndata;

  // Find pixel coordinates of the contributing pixels for tubes of response
  // belonging to first 1/4th of angles. From these pixel coordinates
  // others can be calculated via symmetries in projection space.
  // N.B. The line of response: a*x+b*y+c
  // => solve y: y=-x/tan(view) + s*(cos(theta)/tan(theta) + sin(theta))
  //    solve x: x=-y*tan(theta) + s*(sin(theta)*tan(theta) + cos(theta))
 
  for(view=set; view<=viewNr/4; view+=setNr){
    // view=0 -> sin(theta)=0
    if(view==0){
 
      // Choose column(s) according to sample distance for angles 0 and pi/2.
      col1 = 0;
      col2 = 0;
      for(bin = 0; bin < half; bin++){
 
  col1 = col2;
  col2 = floor((float)(bin + 1)*sdist); 
 
  // Determine factor epsilon.
  if(col1 == col2){ 
    eps1 = sdist;
    eps2 = 0;
    eps3 = 0;
  }
  if((col2-col1) == 1){
    eps1 = (float)(col2 - (bin)*sdist);
    eps2 = 0;
    eps3 = (float)((bin+1)*sdist - col2);
 
  } 
  // If sdist > pixel size!
  if((col2-col1) > 1){
    eps1 = (float)(col1 + 1 - (bin)*sdist);
    eps2 = 1; // middle pixel.
    eps3 = (float)((bin+1)*sdist - col2);
 
  }
 
  /* Iterate through the entries in the image matrix.
     Calculate raysums for two LORs in the same distance from origin 
           (do halfturn). */
  for(row=0; row<imgDim; row++){
    if(!eps3){   
      imgptr[row*imgDim + col1] += eps1 * scnptr[bin]; 
      if(bin != center)
        imgptr[row*imgDim + (imgDim - 1 - col1)]+= 
                eps1 * scnptr[binNr-bin-1] ;
      imgptr[(imgDim - 1 - col1)*imgDim + row] += 
              eps1 * scnptr[binNr*(viewNr/2) + bin];
      if(bin != center)      
        imgptr[col1*imgDim + row]+= 
                eps1 * scnptr[binNr*(viewNr/2) + (binNr-bin-1)];  
    }
    if(eps3 && !eps2){
      imgptr[row*imgDim + col1] += eps1 * scnptr[bin];
      imgptr[row*imgDim + col2] += eps3 *scnptr[bin];
      if(bin != center){ 
        imgptr[row*imgDim + (imgDim - 1 - col1)] += 
                eps1 * scnptr[binNr-bin-1];
        imgptr[row*imgDim+(imgDim-1-col2)] += 
                eps3*scnptr[binNr-bin-1];
      }
 
      imgptr[(imgDim-1-col1)*imgDim+row] += 
              eps1* scnptr[binNr*(viewNr/2) + bin]; 
      imgptr[(imgDim-1-col2)*imgDim+row] += 
              eps3*scnptr[binNr*(viewNr/2) + bin];
      if(bin != center){        
        imgptr[col1*imgDim + row] += 
                eps1 * scnptr[binNr*(viewNr/2) + (binNr-bin-1)];
        imgptr[col2*imgDim + row] += 
                eps3 *scnptr[binNr*(viewNr/2) + (binNr-bin-1)]; 
      }
    }
    if(eps3 && eps2){
      for(col = col1; col<=col2; col++){
        if(col == col1){          
    imgptr[row*imgDim + col1]+= eps1 * scnptr[bin];
     if(bin != center)       
       imgptr[row*imgDim + (imgDim - 1 - col1)]+= 
                     eps1 * scnptr[binNr-bin-1] ;
     imgptr[(imgDim - 1 - col1)*imgDim + row]+= 
                   eps1 * scnptr[binNr*(viewNr/2) + bin];
     if(bin != center)       
       imgptr[col1*imgDim + row]+= 
                     eps1 * scnptr[binNr*(viewNr/2) + (binNr-bin-1)];
        }
        if(col == col2){    
    imgptr[row*imgDim + col2]+= eps3 * scnptr[bin];
    if(bin != center)     
      imgptr[row*imgDim + (imgDim - 1 - col2)]+=
                    eps3 * scnptr[binNr-bin-1];
    imgptr[(imgDim - 1 - col2)*imgDim + row]+=
                  eps3 * scnptr[binNr*(viewNr/2) + bin];
    if(bin != center)     
      imgptr[col2*imgDim + row]+= 
                    eps3 * scnptr[binNr*(viewNr/2) + (binNr-bin-1)];
        }
        if(col != col1 && col != col2){   
    imgptr[row*imgDim + col]+= eps2 * scnptr[bin] ;
    if(bin != center)     
      imgptr[row*imgDim + (imgDim - 1 - col)]+=
                    eps2 * scnptr[binNr-bin-1];
    imgptr[(imgDim - 1 - col)*imgDim + row]+=
                  eps2 * scnptr[binNr*(viewNr/2) + bin];
    if(bin != center)     
      imgptr[col*imgDim + row]+=
                    eps2 * scnptr[binNr*(viewNr/2) + (binNr-bin-1)];
        }
      }
    }
  }
      }
      // End of view==0 (handles angles 0 and pi/2).   
    } else {
 
      // Set sine and cosine for this angle.
      sinus = (double)radonGetSin(radtra,view);
      cosinus = (double)radonGetSin(radtra,viewNr/2 + view);
      tanus = sinus/cosinus;
 
      // Set shift from origin for the first line of response (-n/2,theta).
      // NOTE that image matrix is on cartesian coordinate system where origin
      // is in the middle and shift is in pixels.
      shift_y = -(imgDim/2)/sinus;
      shift_x = -(imgDim/2)/cosinus;
      
      // Evaluate the function of the first LOR in integer points [-n/2,n/2].
      // NOTE that image matrix is on cartesian coordinate system where origin
      // is in the middle.
      z=-imgDim/2;
      for(col=0; col<imgDim+1; col++){
  Yptr[col]=(float)(-z/tanus + shift_y);
    Xptr[col]=(float)(-z*tanus + shift_x);
  z++;
      }
 
      // Set shift from the first LOR.
      shift_y = (double)(sdist/sinus);
      shift_x = (double)(sdist/cosinus);      
      
      // Iterate through half the bins in this view,
      // and determine coordinates of pixels contributing to this LOR.
      // NOTE that shift is added according to 'bin' in every loop.
      // Calculate also the length of intersection.
      // Others are determined via symmetry in projection space.
 
      for(bin=0; bin<half; bin++){
 
  // Limit (x-)indices for fast search.
  // Note that indices are non-negative integers.
 
  x_left = floor((float)(Xptr[imgDim] + bin*shift_x + imgDim/2));
  if(x_left < 0) x_left = 0;
 
  x_right = floor((float)(Xptr[0] + bin*shift_x + imgDim/2));
  if(x_right <= 0) x_right = 1;
  if(x_right > imgDim) x_right = imgDim - 1; 
 
  /* Iterate through the values in vector Y, in integer points 
           [x_left,x_rigth]. */
  for(z=x_left; z <= x_right; z++) {
       
    xp = z; //positive x-coordinate
    xn = imgDim - 1 - xp; //negative x-coordinate
 
    // Look y from left.
    yp = imgDim/2 - floor(Yptr[xp + 1] + bin*shift_y) - 1;
    yn = imgDim - 1 - yp;
 
    // If the y-value for this x (z) is inside the image grid.
    if(yp < imgDim && yp >= 0 && yn < imgDim && yn >= 0 && xp >= 0 && 
             xp < imgDim && xn < imgDim && xn >= 0)
          {
 
      // NOTE that pixels found in this part are always hit from right side.
      // Compute a := |AF| and b := |FB|.
      a = (float)(xp + 1 - (imgDim/2 + Xptr[imgDim - yp] + bin*shift_x));
      b = (float)(floor(Yptr[xp + 1] + bin*shift_y) + 1 - (Yptr[xp + 1] +
                        bin*shift_y));
 
      // Calculate the area of lower triangle.
      A = a*b/2;
 
      // c := |FC|
      c = (float)(xp + 1 - (imgDim/2 + Xptr[imgDim - yp] +
                        (bin+1)*shift_x)); 
 
      if(c > 0){
        // d := |FD|
        d = (float)(floor(Yptr[xp + 1] + (bin+1)*shift_y) + 1 - 
                          (Yptr[xp + 1] + (bin+1)*shift_y));
        // Subtract the area of upper triangle.
        A = A - c*d/2;
      }
 
      eps1 = A;
      if((eps1 < 0 || eps1 > 1) && RADON_VERBOSE){
        printf("RADON: Error in factor: eps1=%.5f \n",eps1);
        errors++;
      }
 
      // Case: theta.
      // Add img(x,y)*k to the raysum of TOR (view,bin)     
      imgptr[yp*imgDim + xp]+= eps1 * scnptr[view*binNr + bin];
      // Add img(-x,-y)*k to the raysum of LOR (view,binNr-bin)
      if(bin != center)       
        imgptr[yn*imgDim + xn]+= 
                eps1 * scnptr[view*binNr + binNr - 1 - bin];
        
      if(view != viewNr/4){
        // Mirror the original LOR on y-axis, i.e. x->-x
        // Case: pi-theta.
        // Add img(-x,y)*k to the raysum of LOR (viewNr-view,bin)       
        imgptr[yp*imgDim + xn]+= 
                eps1 * scnptr[(viewNr - view)*binNr + bin];
        // Add img(x,-y)*k to the raysum of LOR (viewNr-view,binNr-bin)
        if(bin != center)   
    imgptr[yn*imgDim + xp]+=
                  eps1 * scnptr[(viewNr-view)*binNr + binNr - 1 - bin];
        
        // Mirror the LOR on line x=y, i.e. x->y.
        // Case: pi/2-theta
        // Add img(-y,-x)*k to the raysum of LOR (viewNr/2-view,bin)
        imgptr[xn*imgDim + yn]+=
                eps1 * scnptr[(viewNr/2 - view)*binNr + bin];
        // Add img(y,x)*k to the raysum of LOR (viewNr/2-view,binNr-bin)
        if(bin != center)   
    imgptr[xp*imgDim + yp]+=
                  eps1 * scnptr[(viewNr/2 - view)*binNr + binNr - 1 - bin];
      }
 
      // Mirror the LOR on line x=y, and on y-axis i.e x->y and y->-y.
      // Case: pi/2+theta
      // Add img(y,-x)*k to the raysum of LOR (viewNr/2+view,bin)     
      imgptr[xn*imgDim + yp]+= eps1 * scnptr[(viewNr/2 + view)*binNr + bin];
      // Add img(y,-x)*k to the raysum of LOR (viewNr/2+view,binNr-bin)
      if(bin != center)       
        imgptr[xp*imgDim + yn]+= 
                eps1 * scnptr[(viewNr/2 + view)*binNr + binNr - 1 - bin];
    }
  }
  
  // Limit (y-)indices for fast search.
  // Note that indices are non-negative integers.
  y_bottom = floor((float)(Yptr[imgDim] + bin*shift_y + imgDim/2));
  if(y_bottom < 0) y_bottom = 0;
  if(y_bottom > imgDim) y_bottom = 0; 
 
  y_top = floor((float)(Yptr[0] + bin*shift_y + imgDim/2));
  if(y_top > imgDim) y_top = imgDim;
  if(y_top <= 0) y_top = 1;
 
  // Iterate through the values in vector X, in integer points [y_bottom,y_top].
  for(z=y_top; z >= y_bottom; z--) {
   
    // Look y from this location.
    xp = floor(Xptr[z] + bin*shift_x) + imgDim/2;
    xn = imgDim - 1 - xp;
 
    yp = imgDim - z - 1;
    yn = imgDim - yp - 1;
 
    // If the x-value for this y (z) is inside the image grid.
    if(yp < imgDim && yp >= 0 && yn < imgDim && yn >= 0 && xp >= 0 
             && xp < imgDim && xn < imgDim && xn >= 0)
          {
      eps1=eps2=eps3=0;
     
      dx = (float)(Xptr[z] + bin*shift_x + imgDim/2 - xp);
      dy = (float)(Yptr[xp] + bin*shift_y - z + imgDim/2); 
 
      if(dy < 1){ // Cases 3,4,5 and 6.
        // 1. PART
        // a := |HA|
        a = dy;
        // b := |HB| (b < 1 for angles in [0,pi/4))
        b = dx;
        // h := height of rectangle R.
        h = a + shift_y;
        if(h > 1){ // Cases 3,4 and 5.
    h = 1;
    g = b + shift_x;
    if(g > 1){ // Cases 3 and 4.
      g = 1;
      xp2 =floor(Xptr[z+1] + (bin+1)*shift_x) + imgDim/2;
      if(xp == xp2){ // Case 4.
        // c := |FC|
        c = (float)(xp + 1 - (imgDim/2 + Xptr[z+1] + (bin+1)*shift_x));
        // d := |FD| 
        d = (float)(floor(Yptr[xp + 1] + (bin+1)*shift_y) + 1 - 
                                (Yptr[xp + 1] + (bin+1)*shift_y));
        eps1 = 1 - (a*b + c*d)/2;
        eps2 = 0;
        // Lower triangle on the next pixel (x+1,y).
        eps3 = (1 - d)*(b + shift_x - 1)/2; 
        if((eps1<0 || eps1>1 || eps2<0 || eps2>1 || eps3<0 || eps3>1)
                       && RADON_VERBOSE) printf(
                    "4: (%i,%i) eps1=%f | (%i,%i) eps2=%f | (%i,%i) eps3=%f\n",
                    xp,yp,eps1,xp+1,yp-1,eps2,xp+1,yp,eps3);
      } else{ // Case 3.
        // c=d=0.
        eps1 = 1 - a*b/2;
        // Calculate area on pixel in place (xp+1,yp-1).
        dy = (float)(Yptr[xp+1] + (bin+1)*shift_y - (z + 1) + 
                                 imgDim/2);
        if(dy < 1){ // Case 11.
          // c := |HC|
          c = dy;
          // d := |HD|
          d = (float)(xp + 2 - (imgDim/2 + Xptr[z+1] + 
                                  (bin+1)*shift_x));
          eps2 = c*d/2;
        } else { // Cases 9 and 10.
 
          dx = (float)(xp + 2 - (imgDim/2 + Xptr[z+1] + 
                                  (bin+1)*shift_x));
          if(dx < 1) { // Case 10.
      // g := length of rectangle Q.
      g = dx;
      // c := |CD| (on x-axis).
      c = dx - (float)(xp + 2 - (imgDim/2 + Xptr[z+2] + 
                           (bin+1)*shift_x));
      // Rectangle Q - triangle c*h (h = 1).
      eps2 = g - c/2;
          } else { // Case 9.
      // c := |FC|
      c = (float)(xp + 2 - (imgDim/2 + Xptr[z+2] + 
                                   (bin+1)*shift_x));       
      // d := |FD| 
      d = (float)(floor(Yptr[xp + 2] + (bin+1)*shift_y) + 2 - 
                                    (Yptr[xp + 2] + (bin+1)*shift_y));
      // Rectangle Q - triangle CFD.
      eps2 = 1 - c*d/2;
          }
        }
        // Calculate area on pixel in place (xp+1,yp).
        dx = (float)(xp + 2 - (imgDim/2 + Xptr[z] + (bin+1)*shift_x));
        if(dx < 1){ // Case 10.
          // g := length of rectangle Q.
          g = dx;
          // c := |CD| (on x-axis).
          c = dx - (float)(xp + 2 - (imgDim/2 + Xptr[z+1] + 
                          (bin+1)*shift_x));
          // Rectangle Q - triangle c*h (h = 1).
          eps3 = g - c/2;
        } else { // Case 9.
          // c := |FC|
          c = (float)(xp + 2 - (imgDim/2 + Xptr[z+1] + 
                                  (bin+1)*shift_x));        
          // d := |FD| 
          d = (float)(floor(Yptr[xp + 2] + (bin+1)*shift_y) + 1 - 
                                  (Yptr[xp + 2] + (bin+1)*shift_y));
          // Rectangle Q - triangle CFD.
          eps3 = 1 - c*d/2;
        }
        if((eps1<0 || eps1>1 || eps2<0 || eps2>1 || eps3<0 || eps3>1)
                        && RADON_VERBOSE) printf(
                     "3/v%i: (%i,%i) eps1=%f | (%i,%i) eps2=%f | (%i,%i) eps3=%f\n",
                     view,xp,yp,eps1,xp+1,yp-1,eps2,xp+1,yp,eps3);
      }
    } else{ // Case 5. (g < 1)
      // c := |DC|.
      c = g - (imgDim/2 + Xptr[z+1] + (bin+1)*shift_x - xp);
      // d := heigth
      d = 1;
      eps1 = g*h - (a*b + c*d)/2;
      eps2 = 0;
      eps3 = 0;
      if((eps1<0 || eps1>1 || eps2<0 || eps2>1 || eps3<0 || eps3>1) 
                     && RADON_VERBOSE) printf(
                   "5: (%i,%i) eps1=%f | (%i,%i) eps2=%f | (%i,%i) eps3=%f\n",
                   xp,yp,eps1,xp+1,yp-1,eps2,xp+1,yp,eps3);
    }
        } else { // Case 6 (h <= 1). 
    // g := legth of rectangle R
    g = b + shift_x;
    if(g > 1) // Should always be < 1 for angles in [0,pi/4)
      g = 1; 
    eps1 = (g*h - a*b)/2;
    eps2 = 0;
    eps3 = 0;
    if((eps1<0 || eps1>1 || eps2<0 || eps2>1 || eps3<0 || eps3>1) &&
                   RADON_VERBOSE) printf(
                  "6: (%i,%i) eps1=%f | (%i,%i) eps2=%f | (%i,%i) eps3=%f\n",
                  xp,yp,eps1,xp+1,yp-1,eps2,xp+1,yp,eps3);
        }
      } else { // 2. PART
              // Cases 2,7 and 8. (dy >= 1).
        // a := |HA|
        a = 1;
        // b := |HB| (>=1)
        b = dx - (float)(Xptr[z+1] + bin*shift_x + imgDim/2 - xp);
        // h := heigth of rectangle R
        h = 1;
        // g := length of rectangle R
        g = dx + shift_x;
              if(g > 1){ // Cases 2 and 8.
    g = 1 - (float)(Xptr[z+1] + bin*shift_x + imgDim/2 - xp);
                 // positive x-coordinate (bin+1)
    xp2 = floor(Xptr[z+1] + (bin+1)*shift_x) + imgDim/2;
    if(xp == xp2){ // Case 8.
                  // c := |FC|
      c = (float)(xp + 1 - (imgDim/2 + Xptr[z+1] + (bin+1)*shift_x));
                  // d := |FD| 
      d = (float)(floor(Yptr[xp + 1] + (bin+1)*shift_y) + 1 - 
                              (Yptr[xp + 1] + (bin+1)*shift_y));
 
      eps1 = g*h - (a*b + c*d)/2;
      eps2 = 0;
      // Lower triangle on the next pixel (x+1,y).
      eps3 = (1 - d)*((Xptr[z] + (bin+1)*shift_x + imgDim/2) - 
                         (xp+1))/2; 
      if((eps1<0 || eps1>1 || eps2<0 || eps2>1 || eps3<0 || eps3>1) 
                      && RADON_VERBOSE) printf(
                    "8: (%i,%i) eps1=%f | (%i,%i) eps2=%f | (%i,%i) eps3=%f\n",
                     xp,yp,eps1,xp+1,yp-1,eps2,xp+1,yp,eps3);
    } else { // Case 2.
      // c=d=0.
      eps1 = g*h - a*b/2;
 
      /* Pixel in place (xp+1,yp-1) should have been found in 
                     the previous step. */
      eps2 = 0;
 
      // Calculate area on pixel in place (xp+1,yp).
 
      dx = (float)((imgDim/2 + Xptr[z] + (bin+1)*shift_x) - (xp+1));
      if(dx < 1){ // Case 10 (trapezium).
        // g := bottom of trapezium Q.
        g = dx;
                    // c := top of trapezium Q.
        c = (float)((imgDim/2 + Xptr[z+1] + (bin+1)*shift_x) - 
                                (xp+1));
                    // Area of trapezium Q. Heigth = 1.
        eps3 = (g + c)/2;
        if((eps1<0 || eps1>1 || eps2<0 || eps2>1 || eps3<0 || 
                        eps3>1) && RADON_VERBOSE) printf(
                    "2/10: (%i,%i) eps1=%f | (%i,%i) eps2=%f | (%i,%i) eps3=%f\n",
                      xp,yp,eps1,xp+1,yp-1,eps2,xp+1,yp,eps3);
      } else { // Case 9.
          // c := |FC|
          c = (float)(xp + 2 - (imgDim/2 + Xptr[z+1] + 
                                 (bin+1)*shift_x));       
          // d := |FD| 
          d = (float)(floor(Yptr[xp + 2] + (bin+1)*shift_y) + 1 - 
                                  (Yptr[xp + 2] + (bin+1)*shift_y));
          // Rectangle Q - triangle CFD.
          eps3 = 1 - c*d/2;
          if((eps1<0 || eps1>1 || eps2<0 || eps2>1 || eps3<0 || 
                          eps3>1) && RADON_VERBOSE) printf(
                        "2/9: (%i,%i) eps1=%f | (%i,%i) eps2=%f | (%i,%i) eps3=%f\n",
                        xp,yp,eps1,xp+1,yp-1,eps2,xp+1,yp,eps3);
      }
    }
        } else { // Case 7. (g < = 1)
 
    // Area of the parallelogram R.
    eps1 = sdist/cosinus;
    eps2 = 0;
    eps3 = 0;
    if((eps1<0 || eps1>1 || eps2<0 || eps2>1 || eps3<0 || eps3>1) && 
                   RADON_VERBOSE) printf(
                  "7: (%i,%i) eps1=%f | (%i,%i) eps2=%f | (%i,%i) eps3=%f\n",
                  xp,yp,eps1,xp+1,yp-1,eps2,xp+1,yp,eps3);
        }
      }
      if(!eps2 && !eps3) { // Cases 5,6 and 7.
        // Case: theta.
        // Add img(x,y)*k to the raysum of LOR (view,bin)       
        imgptr[yp*imgDim + xp]+= eps1 * scnptr[view*binNr + bin];
        // Add img(-x,-y)*k to the raysum of LOR (view,binNr-bin)
        if(bin != center)   
    imgptr[yn*imgDim + xn]+= 
                  eps1 * scnptr[view*binNr + binNr - 1 - bin];
 
        if(view != viewNr/4) {      
          // Mirror the LOR on y-axis, i.e. x->-x
          // Case: pi-theta.
          // Add img(-x,y)*k to the raysum of LOR (viewNr-view,bin)
    imgptr[yp*imgDim + xn]+= 
                  eps1 * scnptr[(viewNr - view)*binNr + bin];
          // Add img(x,-y)*k to the raysum of LOR (viewNr-view,binNr-bin)
          if(bin != center)   
      imgptr[yn*imgDim + xp]+= 
                    eps1 * scnptr[(viewNr-view)*binNr + binNr - 1 - bin];
        
          // Mirror the LOR on line y=x, i.e. y->x.
          // Case: pi/2 - theta.
          // Add img(y,x)*k to the raysum of LOR (viewNr/2-view,bin)
          imgptr[xn*imgDim + yn]+= 
                  eps1 * scnptr[(viewNr/2 - view)*binNr + bin];
          // Add img(-y,-x)*k to the raysum of LOR (viewNr/2-view,binNr-bin)
          if(bin != center)   
      imgptr[xp*imgDim + yp]+= 
                    eps1 * scnptr[(viewNr/2 - view)*binNr + binNr - 1 - bin];
        }
 
        // Mirror the LOR on line x=y, and on y-axis i.e x->y and y->-y.
        // Case: pi/2 + theta.
        // Add img(y,-x)*k to the raysum of LOR (viewNr/2+view,bin)
        imgptr[xn*imgDim + yp]+= 
                eps1 * scnptr[(viewNr/2 + view)*binNr + bin];
        // Add img(-y,x)*k to the raysum of LOR (viewNr-view,binNr-bin)
        if(bin != center)
    imgptr[xp*imgDim + yn]+= 
                  eps1 * scnptr[(viewNr/2 + view)*binNr + binNr - 1 - bin];
 
      } else {
        if(!eps2) { // <=> eps3 != 0 & eps2 = 0 <=> Cases 3,4 and 8.
    if(xp + 1 < imgDim && xn - 1 >= 0) {
      // Case: theta.
      // Add img(x,y)*k to the raysum of LOR (view,bin)
      imgptr[yp*imgDim + xp] += eps1 * scnptr[view*binNr + bin];
      imgptr[yp*imgDim + xp+1] += eps3 *scnptr[view*binNr + bin];
      // Add img(-x,-y)*k to the raysum of LOR (view,binNr-bin)
      if(bin != center){    
        imgptr[yn*imgDim + xn] += 
                      eps1 * scnptr[view*binNr + binNr - 1 - bin]; 
        imgptr[yn*imgDim + xn-1] += 
                      eps3 * scnptr[view*binNr + binNr - 1 - bin];
      }     
 
      if(view != viewNr/4) {
        // Mirror the LOR on y-axis, i.e. x->-x
        // Case: pi-theta.
        // Add img(-x,y)*k to the raysum of LOR (viewNr-view,bin)
        imgptr[yp*imgDim + xn]+= 
                      eps1 * scnptr[(viewNr - view)*binNr + bin];
        imgptr[yp*imgDim + xn-1] += 
                      eps3 * scnptr[(viewNr - view)*binNr + bin];
        /* Add img(x,-y)*k to the raysum of LOR 
                       (viewNr-view,binNr-bin) */
        if(bin != center) {
          imgptr[yn*imgDim + xp] += 
                        eps1 * scnptr[(viewNr-view)*binNr + binNr - 1 - bin];  
          imgptr[yn*imgDim + xp+1] += 
                        eps3 *scnptr[(viewNr-view)*binNr + binNr - 1 - bin];
        }
        // Mirror the LOR on line y=x, i.e. y->x.
        // Case: pi/2 - theta.
        // Add img(y,x)*k to the raysum of LOR (viewNr/2-view,bin)
        imgptr[xn*imgDim + yn]+= 
                      eps1 * scnptr[(viewNr/2 - view)*binNr + bin];   
        imgptr[(xn-1)*imgDim + yn] += 
                      eps3 *scnptr[(viewNr/2 - view)*binNr + bin];
        /* Add img(-y,-x)*k to the raysum of LOR 
                       (viewNr/2-view,binNr-bin) */
        if(bin != center) {
          imgptr[xp*imgDim + yp]+= 
                        eps1 * scnptr[(viewNr/2 - view)*binNr + binNr - 1 - bin]; 
          imgptr[(xp+1)*imgDim+yp]+=
                        eps3 * scnptr[(viewNr/2 - view)*binNr + binNr - 1 - bin];
        }
      }
 
      // Mirror the LOR on line x=y, and on y-axis i.e x->y and y->-y.
      // Case: pi/2 + theta.
      // Add img(y,-x)*k to the raysum of LOR (viewNr/2+view,bin)
      imgptr[xn*imgDim + yp]+= 
                    eps1 * scnptr[(viewNr/2 + view)*binNr + bin];  
      imgptr[(xn-1)*imgDim + yp] += 
                    eps3 * scnptr[(viewNr/2 + view)*binNr + bin];
      // Add img(-y,x)*k to the raysum of LOR (viewNr-view,binNr-bin)
      if(bin != center){        
        imgptr[xp*imgDim + yn]+= 
                      eps1 * scnptr[(viewNr/2 + view)*binNr + binNr - 1 - bin];
        imgptr[(xp+1)*imgDim+yn] += 
                      eps3*scnptr[(viewNr/2 + view)*binNr + binNr - 1 - bin];
      }
    }
        } else { // <=> eps2!=0 && eps3!=0 <=> Case 3.
    if(xp+1 < imgDim && xn-1 >= 0 && yp-1 >= 0 && yn+1 < imgDim) {
 
      // Case: theta.
      // Add img(x,y)*k to the raysum of LOR (view,bin)
      imgptr[yp*imgDim + xp] += eps1 * scnptr[view*binNr + bin];
      imgptr[yp*imgDim + xp+1] += eps3 *scnptr[view*binNr + bin]; 
      imgptr[(yp-1)*imgDim + xp+1] += 
                    eps2 * scnptr[view*binNr + bin];
      // Add img(-x,-y)*k to the raysum of LOR (view,binNr-bin)
      if(bin != center){
        imgptr[yn*imgDim + xn] += 
                      eps1 * scnptr[view*binNr + binNr - 1 - bin]; 
        imgptr[yn*imgDim + xn-1] += 
                      eps3 * scnptr[view*binNr + binNr - 1 - bin];
        imgptr[(yn+1)*imgDim + xn-1] += 
                      eps2 * scnptr[view*binNr + binNr - 1 - bin];
      }
      if(view != viewNr/4) {      
        // Mirror the LOR on y-axis, i.e. x->-x
        // Case: pi-theta.
        // Add img(-x,y)*k to the raysum of LOR (viewNr-view,bin)
        imgptr[yp*imgDim + xn] +=
                      eps1 * scnptr[(viewNr - view)*binNr + bin];
        imgptr[yp*imgDim + xn-1] +=
                      eps3 *scnptr[(viewNr - view)*binNr + bin]; 
        imgptr[(yp-1)*imgDim + xn-1]+= 
                      eps2 * scnptr[(viewNr - view)*binNr + bin]; 
        /* Add img(x,-y)*k to the raysum of LOR 
                       (viewNr-view,binNr-bin) */
        if(bin != center) {
          imgptr[yn*imgDim + xp] +=
                        eps1 * scnptr[(viewNr-view)*binNr + binNr - 1 - bin]; 
          imgptr[yn*imgDim + xp+1] +=
                        eps3 * scnptr[(viewNr-view)*binNr + binNr - 1 - bin]; 
          imgptr[(yn+1)*imgDim + xp+1] +=
                        eps2 * scnptr[(viewNr-view)*binNr + binNr - 1 - bin];
        }
 
        // Mirror the LOR on line y=x, i.e. y->x.
        // Case: pi/2 - theta.
        // Add img(y,x)*k to the raysum of LOR (viewNr/2-view,bin)
        imgptr[xn*imgDim + yn]+=
                      eps1 * scnptr[(viewNr/2 - view)*binNr + bin];  
        imgptr[xn*imgDim + yn-1] +=
                      eps3 *scnptr[(viewNr/2 - view)*binNr + bin]; 
        imgptr[(xn-1)*imgDim + (yn+1)] +=
                      eps2 * scnptr[(viewNr/2 - view)*binNr + bin];
        /* Add img(-y,-x)*k to the raysum of LOR 
                      (viewNr/2-view,binNr-bin) */
        if(bin != center) {
          imgptr[xp*imgDim + yp]+= 
                        eps1 * scnptr[(viewNr/2 - view)*binNr + binNr - 1 - bin]; 
          imgptr[(xp+1)*imgDim+yp] +=
                        eps3*scnptr[(viewNr/2 - view)*binNr + binNr - 1 - bin];
          imgptr[(xp+1)*imgDim+(yp-1)] +=
                        eps2*scnptr[(viewNr/2 - view)*binNr + binNr - 1 - bin];
        }
      }
 
      // Mirror the LOR on line x=y, and on y-axis i.e x->y and y->-y.
      // Case: pi/2 + theta.
      // Add img(y,-x)*k to the raysum of LOR (viewNr/2+view,bin).
      imgptr[xn*imgDim + yp] +=
                    eps1 * scnptr[(viewNr/2 + view)*binNr + bin]; 
      imgptr[(xn-1)*imgDim + yp] +=
                    eps3 * scnptr[(viewNr/2 + view)*binNr + bin];
      imgptr[(xn-1)*imgDim + (yp-1)] +=
                    eps2 *scnptr[(viewNr/2 + view)*binNr + bin];
      // Add img(-y,x)*k to the raysum of LOR (viewNr-view,binNr-bin)
      if(bin != center){        
        imgptr[xp*imgDim + yn]+=
                      eps1 * scnptr[(viewNr/2 + view)*binNr + binNr - 1 - bin]; 
        imgptr[(xp+1)*imgDim+yn] +=
                      eps3 * scnptr[(viewNr/2 + view)*binNr + binNr - 1 - bin];
        imgptr[(xp+1)*imgDim+yn+1] +=
                      eps2*scnptr[(viewNr/2 + view)*binNr + binNr - 1 - bin];
      }
    }
        }
      } 
    }
  }
      } // END of X loop
    } // End of view>0.
  } // END OF VIEW LOOP
  free(X);
  free(Y);
 
  if(RADON_VERBOSE) {
    printf("RADON: radonBackTransformEA() finished with %d errors.\n", errors);
    fflush(stdout);
  }
 
  return 0;
 
} // END OF EXACT AREA BACK TRANSFORM.
/*****************************************************************************/
int radonBackTransformSA(
  RADON *radtra, int set, int setNr, float *scndata, float *imgdata
) {
  float *imgptr, *scnptr, *imgorigin;
  float sinus, cosinus; 
  float t;   // distance between projection ray and origo.
  int   mode, imgDim, binNr, viewNr, halfImg, centerBin, view;
  int x, y, xright, ybottom;
  float tpow2;
 
  // Retrieve the data from given radon transform object.
  mode=radonGetMO(radtra); // Should be 3 or 4.
  imgDim=radonGetID(radtra);
  binNr=radonGetNB(radtra); // Should be equal to imgDim.
  viewNr=radonGetNV(radtra);
  // Set the center ray and the square of it.
  centerBin=binNr/2;
  tpow2 = centerBin*centerBin;
 
  if(imgDim != binNr) return -1;
 
  // Set half of the image dimension.
  halfImg = imgDim/2;
  
  // Transform one angle at a time.
  for(view=set; view<viewNr; view+=setNr){
 
    imgorigin = imgdata + imgDim*(halfImg - 1) + halfImg;
 
    sinus = radonGetSin(radtra,view);
    cosinus = radonGetSin(radtra,viewNr/2 + view);
 
    y = halfImg - 2;
    
    if((float)y > centerBin){
      y = (int)(centerBin);
      ybottom = -y;
    } else{
      ybottom = -halfImg + 1;
    }
    for(; y >= ybottom; y--){
      xright = 
        (int)sqrt(/*fabs(*/tpow2 - ((float)y+0.5) * ((float)y+0.5)/*)*/) + 1;
      if(xright >= halfImg){
  xright = halfImg - 1;
  x = -halfImg;
      }
      else
  x = -xright;

      //      t = centerBin - (float)y * sinus + ((float)(x + 1)) * cosinus;
      t = centerBin + (float)y * sinus + ((float)(x + 1)) * cosinus;
      
      imgptr = imgorigin - y*imgDim + x;
      scnptr = scndata + view*binNr;
      for(; x <= xright; x++, t += cosinus){
  if(mode == 3){ // If the linear interpolation is to be utilised.
      *imgptr++ += 
              ( *(scnptr+(int)t) + (*(scnptr+(int)t+1) - *(scnptr+(int)t))
     * (t - (float)(int)t) );
  } else {// If the nearest neighbour interpolation is to be utilised.
    *imgptr++ += *(scnptr+(int)(t + 0.5)); /* (float)(1.0 / (float)views)*/
        }
      }
    }
  }
  return 0;
}
/*****************************************************************************/
int radonBackTransformPRM(
  PRMAT *mat, int set, int setNr, float *scndata, float *imgdata
) {
  float *imgptr, *scnptr; // pointers for the image and sinogram data
  unsigned int imgDim=128, binNr=256, viewNr=192, view, bin, row, col;
  unsigned int half, centerbin;
  int xp, xn, yp, yn;
  float fact;
 
  // Retrieve the data from given projection matrix.
  imgDim=prmatGetID(mat);
  binNr=prmatGetNB(mat); // Should be equal to imgDim.
  viewNr=prmatGetNV(mat);  
 
  // Calculate and set center bin for current geometrics.
  if((binNr%2) != 0){
    half = (binNr - 1)/2 + 1;
    centerbin = half - 1;
  } else {
    half = binNr/2;
    // In the case binNr is even there is no center bin.
    centerbin = -1;
  }
 
  imgptr = imgdata;
  scnptr = scndata;
 
  // Draw sinogram according to projection matrix.
  for(view=set; view<=viewNr/4; view=view+setNr){
    for(bin=0; bin<half; bin++){
      row = view*half + bin;
      for(col=0; col<prmatGetPixels(mat,row); col++){
 
  fact = prmatGetFactor(mat,row,col);
  xp = prmatGetXCoord(mat,row,col);
  xn = imgDim - 1 - xp;
  yp = prmatGetYCoord(mat,row,col);
  yn = imgDim - 1 - yp;
      
  // Add img(x,y)*k to the raysum of LOR (view,bin)
  imgptr[yp*imgDim + xp] += fact * scnptr[view*binNr + bin];
  if(bin != centerbin)
    // Add img(-x,-y)*k to the raysum of LOR (view,binNr-bin)
    imgptr[yn*imgDim + xn] += fact * scnptr[view*binNr + binNr - 1 - bin];
 
  if(view != 0 && view != viewNr/4){        
    // Mirror the original LOR on y-axis, i.e. x->-x
    // Case: pi-theta.
    // Add img(-x,y)*k to the raysum of LOR (viewNr-view,bin)
    imgptr[yp*imgDim + xn] += fact * scnptr[(viewNr - view)*binNr + bin];
    // Add img(x,-y)*k to the raysum of LOR (viewNr-view,binNr-bin)
    if(bin != centerbin)
      imgptr[yn*imgDim + xp] += 
              fact * scnptr[(viewNr-view)*binNr + binNr - 1 - bin];
        
    // Mirror the LOR on line x=y, i.e. x->y.
    // Case: pi/2-theta
    // Add img(-y,-x)*k to the raysum of LOR (viewNr/2-view,bin)
    imgptr[xn*imgDim + yn] += fact * scnptr[(viewNr/2 - view)*binNr + bin];
    // Add img(y,x)*k to the raysum of LOR (viewNr/2-view,binNr-bin)
    if(bin != centerbin)
      imgptr[xp*imgDim + yp] += 
              fact * scnptr[(viewNr/2 - view)*binNr + binNr - 1 - bin];
  }
 
  // Mirror the LOR on line x=y, and on y-axis i.e x->y and y->-y.
  // Case: pi/2+theta
  // Add img(y,-x)*k to the raysum of LOR (viewNr/2+view,bin)
  imgptr[xn*imgDim + yp] += fact * scnptr[(viewNr/2 + view)*binNr + bin];
  // Add img(y,-x)*k to the raysum of LOR (viewNr/2+view,binNr-bin)
  if(bin != centerbin)
    imgptr[xp*imgDim + yn] += 
            fact * scnptr[(viewNr/2 + view)*binNr + binNr - 1 - bin];
      
      } 
    }
  }
 
  return 0;
}
/*****************************************************************************/
/* PROCEDURES FOR SETTING PROJECTION MATRIX IN THE PROJECTION MATRIX 
   DATA STRUCTURE */
/*****************************************************************************/
int radonSetBases(
  RADON *radtra,ELLIPSE *elli, PRMAT *mat
) {
  // temporary pointers for storing the coordinates and factors
  unsigned short int **coords, *factors, **coptr, *facptr; 
  // pointers for the values of LOR in integer points
  float *Y, *X, *Xptr, *Yptr;
  double sinus, cosinus, tanus; // sin, cosine and tangent
  double shift_x, shift_y; // shift for LOR
  int xp, xn, yp, yn, z; //integer points
  int x_left, x_right, y_top, y_bottom; // limits for fast search
  int col, row, view, bin, pix; //counters
  int binNr, viewNr, imgDim, views, rows, mode, half; // properties
  float diam, sdist;
  float scale, sum=0, sqr_sum=0, min=0, max=0;
  float dx, dy , loi = 1;  // delta x and y and length of intersection
 
  // Check that the data for this radon transform is initialized.
  if(radtra->status != RADON_STATUS_INITIALIZED) return -1;
 
  // Retrieve the data from given radon transform object.
  mode=radonGetMO(radtra);
  imgDim=radonGetID(radtra);
  binNr=radonGetNB(radtra);
  viewNr=radonGetNV(radtra);
  sdist=radonGetSD(radtra);
  half=radonGetHI(radtra);
 
  // Set information in the projection matrix data stucture.
  if(binNr == 192) {
    mat->type=PRMAT_TYPE_ECAT931;
  } else if(binNr == 281){
    mat->type=PRMAT_TYPE_GE;
  } else
    mat->type=PRMAT_TYPE_NA;
 
  mat->viewNr=viewNr;
  mat->binNr=binNr;
  mat->imgDim=imgDim;
 
  mat->mode = mode;
  // Allocate memory in PRMAT struct and set fov.
  mat->prmatfov = (int*)calloc(2,sizeof(int));
  mat->prmatfov[0] = ellipseGetMajor(elli);
  mat->prmatfov[1] = ellipseGetMinor(elli);
 
  // rows := number of lines in the base set = (views*half)
  views = viewNr/4 + 1;
  rows = views*half;
 
  mat->factor_sqr_sum=calloc(rows,sizeof(float));
  mat->dimr=rows; 
  //entries in one line
  mat->dime=calloc(rows,sizeof(int)); 
  mat->_factdata=(unsigned short int***)calloc(rows,sizeof(unsigned short int**)); 
  //if not enough memory
  if(!mat->_factdata || !mat->dime) return(4);
  
  // Compute scaling factor, assuming that maximum value is 1.
  scale=65534./1; 
  mat->scaling_factor=1./scale; //is the inverse of scale
 
  /* Allocate (temporary) coordinate and factor arrays to maximum size 
    (= 2*imgDim).
     Coords contains pairs (x,y), coordinates in cartesian coordinate system. */
  coords = (unsigned short int**)calloc(2*imgDim,sizeof(unsigned short int*));
  factors = (unsigned short int*)calloc(2*imgDim,sizeof(unsigned short int));
  if(!coords || !factors) return 5;
 
  for(col=0; col<2*imgDim; col++){
    coords[col] = (unsigned short int*)calloc(2,sizeof(unsigned short int));
    if(!coords[col]) return -1;
  }
 
  coptr = coords;
  facptr = factors;

  // Set diameter of a pixel.
  diam = sqrt(2);
 
  // Array for storing values A_x.
  X=(float*)calloc(imgDim+1,sizeof(float));  
  // Array for storing values A_y.
  Y=(float*)calloc(imgDim+1,sizeof(float));  
  if(X==NULL || Y==NULL){
    return -2;
  }
  Xptr=X;
  Yptr=Y;

  // Find pixel coordinates of the contributing pixels for lines of response
  // belonging to first 1/4th of angles. From these pixel coordinates
  // others can be calculated via symmetries in projection space.
  // N.B. The line of response: a*x+b*y+c
  // => solve y: y=-x/tan(view) + s*(cos(theta)/tan(theta) + sin(theta))
  //    solve x: x=-y*tan(theta) + s*(sin(theta)*tan(theta) + cos(theta))
 
  for(view=0; view<views; view++){
     // Angle theta=0 -> sin(theta)=0.
    if(view==0){
 
      //If the mode switch is on zero do the 0-1 model transform.
      if(mode == 0)
  loi = 1.;
      else
  loi=1/diam;
      
      // Choose column according to sample distance for angle 0.
      col = 0;
      for(bin=0; bin<half; bin++){
  col = floor((float)(bin+.5*sdist)*sdist); 
 
  pix = 0; // length of list pixels (and factors)
  sqr_sum = 0;
 
  // Iterate through the entries in the image matrix.
  // Set factors for pixels intersected by lines in the base set.
  for(row=0; row<imgDim; row++){
     
    /* Check that pixel and pixels hit by symmetrical lines lay inside 
             the FOV.  */
    if(ellipseIsInside(elli,row,col) && ellipseIsInside(elli,row,imgDim - 
             1 - col) && ellipseIsInside(elli,imgDim-1-col,row) && 
             ellipseIsInside(elli,col,row))
          {  
      coptr[pix][0] = (unsigned short int)col;
      coptr[pix][1] = (unsigned short int)row;
      // Convert loi to unsigned short int.
      facptr[pix] = (unsigned short int)(scale*loi);
      // Look for minimal and maximal factor.
      if(min>loi) min=loi;
      if(max<loi) max=loi;
      // Compute square sums in every row.
      sqr_sum += loi*loi; 
      sum += loi;
      pix++;
    }
  }
 
        /* Allocate memory in factor pointer (dynamically according to number 
           of pixels intersected). */
  if(pix){
    mat->_factdata[bin]=
            (unsigned short int**)calloc(pix,sizeof(unsigned short int*)); 
    if(!mat->_factdata[bin]) return -1;
  }
 
  // Allocate leaves.
  for(col=0; col<pix; col++) {
    mat->_factdata[bin][col]=
            (unsigned short int*)calloc(3,sizeof(unsigned short int)); 
    if(!mat->_factdata[bin][col]) return -1;
  }
 
  // Put now values in coordinates and factors array to result pointer.
  mat->fact = mat->_factdata;
  for(col=0; col<pix; col++) {
    mat->fact[bin][col][0] = coptr[col][0]; // x-coodinate
    mat->fact[bin][col][1] = coptr[col][1];  // y-coordinate
    mat->fact[bin][col][2] = facptr[col];  // factor
  }
  
  //Set also the number of pixels belonging to each row and square sums.
  mat->dime[bin]=pix;
  mat->factor_sqr_sum[bin]=sqr_sum;
      
      } // END OF BIN-LOOP FOR ANGLE 0.     
      // End of view=0.
    } else {
 
      // Set sine and cosine for this angle.
      sinus = (double)radonGetSin(radtra,view);
      cosinus = (double)radonGetSin(radtra,viewNr/2 + view);
      tanus = sinus/cosinus;
 
      // Set shift from origin for the first line of response (-n/2,theta).
      // NOTE that image matrix is on cartesian coordinate system where origin
      // is in the middle and shift is in pixels.
      shift_y = -(imgDim/2 -.5*sdist)/sinus;
      shift_x = -(imgDim/2 -.5*sdist)/cosinus;
 
      // Evaluate the function of the first LOR in integer points [-n/2,n/2].
      // NOTE that image matrix is on cartesian coordinate system where origin
      // is in the middle.
      z=-imgDim/2;
      for(col=0; col<imgDim+1; col++){
  Yptr[col]=(float)(shift_y - z/tanus);
    Xptr[col]=(float)(shift_x - z*tanus);
  z++;
      }
 
      // Set shift from the first LOR.
      shift_y = (double)(sdist/sinus);
      shift_x = (double)(sdist/cosinus);
 
      // Iterate through half the bins in this view,
      // and determine coordinates of pixels contributing to this LOR.
      // NOTE that shift is added according to 'bin' in every loop.
      // Calculate also the length of intersection.
      // Others are determined via symmetry in projection space.
 
      for(bin=0; bin<half; bin++){
 
  coptr = coords;
  facptr = factors;
  pix = 0;
  sqr_sum = 0;
  
  // Limit (x-)indices for fast search.
  // Note that indices are non-negative integers.
 
  x_left = floor((float)(Xptr[imgDim] + bin*shift_x + imgDim/2));
  if(x_left < 0) x_left = 0;
 
  x_right = floor((float)(Xptr[0] + bin*shift_x + imgDim/2));
  if(x_right <= 0) x_right = 1;
  if(x_right > imgDim) x_right = imgDim - 1; 
 
  /* Iterate through the values in vector Y, in integer points 
           [x_left,x_rigth]. */
  for(z=x_left; z <= x_right; z++) {
       
    xp = z; //positive x-coordinate
    xn = imgDim - 1 - xp; //negative x-coordinate
 
    // Look y from left.
    yp = imgDim/2 - floor(Yptr[xp + 1] + bin*shift_y) - 1;
    yn = imgDim - 1 - yp;
 
    // If the y-value for this x (z) is inside the image grid.
    if(yp < imgDim && yp >= 0 && yn < imgDim && yn >= 0 && xp >= 0 && 
             xp < imgDim && xn < imgDim && xn >= 0) {
 
      if(!mode) {
        loi = 1;
      } else {
        // Compute delta x and y.
        dx = (float)(xp + 1 - (imgDim/2 + Xptr[imgDim - yp] + bin*shift_x));
        dy = (float)(floor(Yptr[xp + 1] + bin*shift_y) + 1 - 
                          (Yptr[xp + 1] + bin*shift_y));
        if(dx > 1 || dx < 0) dx = 1;
        if(dy > 1 || dy < 0) dy = 1;
        loi = sqrt(dx*dx + dy*dy);
        loi=loi/diam;
      }
 
      /* Check that pixel and pixels hit by symmetrical lines lay inside 
               the FOV.  */
      if(ellipseIsInside(elli,yp,xp) && ellipseIsInside(elli,yn,xn) &&
         ellipseIsInside(elli,yp,xn) && ellipseIsInside(elli,yn,xp) &&
         ellipseIsInside(elli,xn,yn) && ellipseIsInside(elli,xp,yp) &&
         ellipseIsInside(elli,xn,yp) && ellipseIsInside(elli,xp,yn)) {
   
        // Put coordinates and factors in the list.
        coptr[pix][0] = xp;
        coptr[pix][1] = yp;
        facptr[pix] = (unsigned short int)(scale*loi);
       
        // Look for minimal and maximal factor.
        if(min>loi) min=loi;
        if(max<loi) max=loi;
        sqr_sum += loi*loi;             
        sum += loi;
        pix++;             
      }
    }
  }
 
  // Limit (y-)indices for fast search.
  // Note that indices are non-negative integers.
  y_bottom = floor((float)(Yptr[imgDim] + bin*shift_y + imgDim/2));
  if(y_bottom < 0) y_bottom = 0;
  if(y_bottom > imgDim) y_bottom = 0; 
 
  y_top = floor((float)(Yptr[0] + bin*shift_y + imgDim/2));
  if(y_top > imgDim) y_top = imgDim-1;
  if(y_top <= 0) y_top = 1;
 
  /* Iterate through the values in vector X, in integer points 
           [y_bottom,y_top]. */
  for(z=y_top; z >= y_bottom; z--) {
   
    // Look y from this location.
    xp = floor(Xptr[z] + bin*shift_x) + imgDim/2;
    xn = imgDim - 1 - xp;
 
    yp = imgDim - z - 1;
    yn = imgDim - yp - 1;
 
    // If the x-value for this y (z) is inside the image grid.
    if(yp < imgDim && yp >= 0 && yn < imgDim && yn >= 0 && xp >= 0 && 
             xp < imgDim && xn < imgDim && xn >= 0) {
     
      if(!mode){
        loi = 1;
      } else {
        dx = (float)(Xptr[z] + bin*shift_x + imgDim/2 - xp);
        dy = (float)(Yptr[xp] + bin*shift_y - z + imgDim/2); 
        if(dy > 1 || dy < 0) {
    dx = dx - (float)(Xptr[z+1] + bin*shift_x + imgDim/2 - xp);
    dy = 1;
        }
        loi = sqrt(dx*dx + dy*dy)/diam;       
      }
 
      /* Check that pixel and pixels hit by symmetrical lines lay inside 
               the FOV.  */
      if(ellipseIsInside(elli,yp,xp) && ellipseIsInside(elli,yn,xn) &&
         ellipseIsInside(elli,yp,xn) && ellipseIsInside(elli,yn,xp) &&
         ellipseIsInside(elli,xn,yn) && ellipseIsInside(elli,xp,yp) &&
         ellipseIsInside(elli,xn,yp) && ellipseIsInside(elli,xp,yn)){
 
        // Put coordinates and factors in the list.
        coptr[pix][0] = xp;
        coptr[pix][1] = yp;
        facptr[pix] = (unsigned short int)(scale*loi);
     
        // Look for minimal and maximal factor.
        if(min>loi) min=loi;
        if(max<loi) max=loi;
        sqr_sum += loi*loi;          
        sum += loi;
        pix++;        
      }
    }       
  }
  // Allocate memory in result pointer.
  if(pix){
    mat->_factdata[view*half + bin]=
            (unsigned short int**)calloc(pix,sizeof(unsigned short int*)); 
    if(!mat->_factdata[view*half + bin]) return -1;
  }
 
  // Allocate leaves.
  for(col=0; col<pix; col++){
    mat->_factdata[view*half + bin][col]=
            (unsigned short int*)calloc(3,sizeof(unsigned short int)); 
    if(!mat->_factdata[view*half + bin][col]) return -1;
  }
 
  // Put now values in coordinates and factors array to result pointer.
  mat->fact = mat->_factdata;
  for(col=0; col<pix; col++){
    mat->fact[view*half + bin][col][0] = coptr[col][0]; 
    mat->fact[view*half + bin][col][1] = coptr[col][1]; 
    mat->fact[view*half + bin][col][2] = facptr[col];  
  }
 
  // Update also lists mat->dime and mat->factor_sqr_sums
  mat->dime[view*half + bin]=pix;
 
  mat->factor_sqr_sum[view*half + bin]=sqr_sum;
  
      }// END of bin loop
       
    }// End of view>0.
   
  }// END OF VIEW LOOP
 
  free(X);
  free(Y);
  free(coords);
  free(factors);
 
  mat->min = min;
  mat->max = max;
  mat->factor_sum = sum;
  mat->status = PRMAT_STATUS_BS_OCCUPIED;
  return 0;
 
} // END OF SETTING BASE LINES.
/*****************************************************************************/
int radonSetBasesEA(RADON *radtra,ELLIPSE *elli, PRMAT *mat)
{
  // temporary pointers for storing the coordinates and factors
  unsigned short int **coords, *factors, **coptr, *facptr;
  // pointers for the values of LOR in integer points
  float *Y, *X, *Xptr, *Yptr;
  double sinus, cosinus, tanus; // sin, cosine and tangent
  double shift_x, shift_y; // shift for LOR
  int xp, xn, yp, yn, z, xp2; //integer points
  int x_left, x_right, y_top, y_bottom; // limits for fast search
  int col, col1, col2, row, view, bin, pix, errors=0; //counters
  int binNr, viewNr, imgDim, mode,  views, rows, half;
  float sdist;
  float a=0,b=0, c=0, d=0, g=0, h=0, A;
  // delta x and y and factors epsilon.
  float dx, dy , eps1 = 0, eps2 = 0, eps3 = 0;
  float min=0, max=0, sum=0, scale, sqr_sum=0;
 
  // Check that the data for this radon transform is initialized.
  if(radtra->status != RADON_STATUS_INITIALIZED) return -1;
 
  // Retrieve the data from given radon transform object.
  mode=radonGetMO(radtra); // Should be 2.
  imgDim=radonGetID(radtra);
  binNr=radonGetNB(radtra);
  viewNr=radonGetNV(radtra);
  sdist=radonGetSD(radtra);
  half=radonGetHI(radtra);
 
  /*printf("radonSetBases(): m=%i, id=%i, b=%i, v=%i, sd=%f, h=%i.\n",
           mode,imgDim,binNr,viewNr,sdist,half);*/
 
  // Set information in the projection matrix data stucture.
  if(binNr == 192){
    mat->type=PRMAT_TYPE_ECAT931;
  } else if(binNr == 281) {
    mat->type=PRMAT_TYPE_GE;
  } else
    mat->type=PRMAT_TYPE_NA;
 
  mat->viewNr=viewNr;
  mat->binNr=binNr;
  mat->imgDim=imgDim;
 
  mat->mode = mode;
  // Allocate memory in PRMAT struct and set fov.
  mat->prmatfov = (int*)calloc(2,sizeof(int));
  mat->prmatfov[0] = ellipseGetMajor(elli);
  mat->prmatfov[1] = ellipseGetMinor(elli);
 
  // rows := number of lines in the base set = (views*half)
  views = viewNr/4 + 1;
  rows = views*half;
 
  mat->factor_sqr_sum=calloc(rows,sizeof(float));
  mat->dimr=rows; 
  //entries in one line
  mat->dime=calloc(rows,sizeof(int)); 
  mat->_factdata=(unsigned short int***)calloc(rows,sizeof(unsigned short int**)); 
  //if not enough memory
  if(!mat->_factdata || !mat->dime) return(4);
  
  // Compute scaling factor, assuming that maximum value is 1.
  scale=65534./1; 
  mat->scaling_factor=1.0/scale; //is the inverse of scale
 
  /* Allocate (temporary) coordinate and factor arrays to maximum size 
     (= 4*imgDim). Coords contains pairs (x,y), coordinates in cartesian 
      coordinate system. */
  coords = (unsigned short int**)calloc(4*imgDim,sizeof(unsigned short int*));
  factors = (unsigned short int*)calloc(4*imgDim,sizeof(unsigned short int));
  if(!coords || !factors) return 5;
 
  for(col=0; col<4*imgDim; col++){
    coords[col] = (unsigned short int*)calloc(2,sizeof(unsigned short int));
    if(!coords[col]) return -1;
  }
 
  coptr = coords;
  facptr = factors;
 
  // Array for storing values f_x.
  X=(float*)calloc(imgDim+1,sizeof(float));  
  // Array for storing values f_y.
  Y=(float*)calloc(imgDim+1,sizeof(float));  
  if(X==NULL || Y==NULL){
    return -2;
  }
  Xptr=X;
  Yptr=Y;

  // Find pixel coordinates of the contributing pixels for tubes of response
  // belonging to first 1/4th of angles. From these pixel coordinates
  // others can be calculated via symmetries in projection space.
  // N.B. The line of response: a*x+b*y+c
  // => solve y: y=-x/tan(view) + s*(cos(theta)/tan(theta) + sin(theta))
  //    solve x: x=-y*tan(theta) + s*(sin(theta)*tan(theta) + cos(theta))
  for(view=0; view<views; view++){
    // view=0 -> sin(theta)=0
    if(view==0){
 
      // Choose column(s) according to sample distance for angles 0 and pi/2.
      col1 = 0;
      col2 = 0;
      for(bin = 0; bin < half; bin++){
 
  pix = 0; // length of list pixels (and factors)
  sqr_sum = 0;
  
  col1 = col2;
  col2 = floor((float)(bin + 1)*sdist); 
  
  // Determine factor epsilon.
  if(col1 == col2){ 
    eps1 = sdist;
    eps2 = 0;
    eps3 = 0;
  }
  if((col2-col1) == 1){
    eps1 = (float)(col2 - (bin)*sdist);
    eps2 = 0;
    eps3 = (float)((bin+1)*sdist - col2);     
  } 
  // If sdist > pixel size!
  if((col2-col1) > 1){
    eps1 = (float)(col1 + 1 - (bin)*sdist);
    eps2 = 1; // middle pixel.
    eps3 = (float)((bin+1)*sdist - col2);
  }
 
  // Iterate through the entries in the image matrix.
  for(row=0; row<imgDim; row++){
 
    /* Check that pixel and pixels hit by symmetrical lines lay inside 
             the FOV.  */
    if(ellipseIsInside(elli,row,col1) && 
             ellipseIsInside(elli,row,imgDim - 1 - col1) &&
       ellipseIsInside(elli,imgDim-1-col1,row) && 
             ellipseIsInside(elli,col1,row)) { 
 
      if(!eps3){      
        coptr[pix][0] = (unsigned short int)col1;
        coptr[pix][1] = (unsigned short int)row;
        // Convert eps to unsigned short int.
        facptr[pix] = (unsigned short int)(scale*eps1);
        // Look for minimal and maximal factor.
        if(min>eps1) min=eps1;
        if(max<eps1) max=eps1;
        // Compute square sums in every row.
        sqr_sum += eps1*eps1; 
        sum += eps1;
        pix++;
      }
 
      if(eps3 && !eps2){
        coptr[pix][0] = (unsigned short int)col1;
        coptr[pix][1] = (unsigned short int)row;
        // Convert eps to unsigned short int.
        facptr[pix] = (unsigned short int)(scale*eps1);
        // Look for minimal and maximal factor.
        if(min>eps1) min=eps1;
        if(max<eps1) max=eps1;
        // Compute square sums in every row.
        sqr_sum += eps1*eps1; 
        sum += eps1;
        pix++;
 
        coptr[pix][0] = (unsigned short int)col2;
        coptr[pix][1] = (unsigned short int)row;
        // Convert eps to unsigned short int.
        facptr[pix] = (unsigned short int)(scale*eps3);
        // Look for minimal and maximal factor.
        if(min>eps3) min=eps3;
        if(max<eps3) max=eps3;
        // Compute square sums in every row.
        sqr_sum += eps3*eps3; 
        sum += eps3;
        pix++;
      }
 
      if(eps3 && eps2){
        for(col = col1; col<=col2; col++){
 
    if(col == col1){
      coptr[pix][0] = (unsigned short int)col1;
      coptr[pix][1] = (unsigned short int)row;
      // Convert eps to unsigned short int.
      facptr[pix] = (unsigned short int)(scale*eps1);
      // Look for minimal and maximal factor.
      if(min>eps1) min=eps1;
      if(max<eps1) max=eps1;
      // Compute square sums in every row.
      sqr_sum += eps1*eps1; 
      sum += eps1;
      pix++;
    }
    
    if(col == col2){
      coptr[pix][0] = (unsigned short int)col2;
      coptr[pix][1] = (unsigned short int)row;
      // Convert eps to unsigned short int.
      facptr[pix] = (unsigned short int)(scale*eps3);
      // Look for minimal and maximal factor.
      if(min>eps3) min=eps3;
      if(max<eps3) max=eps3;
      // Compute square sums in every row.
      sqr_sum += eps3*eps3; 
      sum += eps3;
      pix++;
    }
 
    if(col != col1 && col != col2){
      coptr[pix][0] = (unsigned short int)col;
      coptr[pix][1] = (unsigned short int)row;
      // Convert eps to unsigned short int.
      facptr[pix] = (unsigned short int)(scale*eps2);
      // Look for minimal and maximal factor.
      if(min>eps2) min=eps2;
      if(max<eps2) max=eps2;
      // Compute square sums in every row.
      sqr_sum += eps2*eps2; 
      sum += eps2;
      pix++;
    }
    }
       }
    }// Pixel is inside the fov
  }// END OF ROW-LOOP
 
        /* Allocate memory in factor pointer (dynamically according to number 
           of pixels intersected). */
  if(pix){
    mat->_factdata[bin]=  // 0 angle
           (unsigned short int**)calloc(pix,sizeof(unsigned short int*));
    if(!mat->_factdata[bin]) return -1;
  }
  // Allocate leaves.
  for(col=0; col<pix; col++){
    mat->_factdata[bin][col]=
           (unsigned short int*)calloc(3,sizeof(unsigned short int)); 
    if(!mat->_factdata[bin][col]) return -1;
  }
 
  // Put now values in coordinates and factors array to result pointer.
  mat->fact = mat->_factdata;
  for(col=0; col<pix; col++){
    mat->fact[bin][col][0] = coptr[col][0]; // x-coodinate
    mat->fact[bin][col][1] = coptr[col][1];  // y-coordinate
    mat->fact[bin][col][2] = facptr[col];  // factor
  }
  
  //Set also the number of pixels belonging to each row and square sums.
  mat->dime[bin]=pix;
  mat->factor_sqr_sum[bin]=sqr_sum;
 
      } // END OF BIN-LOOP FOR ANGLE 0
      // End of view==0 (handles angles 0,pi/4,pi/2,3pi/4).
    } else { // Angles in (0,pi)
 
      // Set sine and cosine for this angle.
      sinus = (double)radonGetSin(radtra,view);
      cosinus = (double)radonGetSin(radtra,viewNr/2 + view);
      tanus = sinus/cosinus;
 
      // Set shift from origin for the first line of response (-n/2,theta).
      // NOTE that image matrix is on cartesian coordinate system where origin
      // is in the middle and shift is in pixels.
      shift_y = -(imgDim/2)/sinus;
      shift_x = -(imgDim/2)/cosinus;
      
      // Evaluate the function of the first LOR in integer points [-n/2,n/2].
      // NOTE that image matrix is on cartesian coordinate system where origin
      // is in the middle.
      z=-imgDim/2;
      for(col=0; col<imgDim+1; col++){
  Yptr[col]=(float)(-z/tanus + shift_y);
    Xptr[col]=(float)(-z*tanus + shift_x);
  z++;
      }
 
      // Set shift from the first LOR.
      shift_y = (double)(sdist/sinus);
      shift_x = (double)(sdist/cosinus);      
 
      // Iterate through half the bins in this view,
      // and determine coordinates of pixels contributing to this LOR.
      // NOTE that shift is added according to 'bin' in every loop.
      // Calculate also the length of intersection.
      // Others are determined via symmetry in projection space.
 
      for(bin=0; bin<half; bin++){
 
  pix = 0;
  sqr_sum = 0;
 
  // Limit (x-)indices for fast search.
  // Note that indices are non-negative integers.
 
  x_left = floor((float)(Xptr[imgDim] + bin*shift_x + imgDim/2));
  if(x_left < 0) x_left = 0;
 
  x_right = floor((float)(Xptr[0] + bin*shift_x + imgDim/2));
  if(x_right <= 0) x_right = 1;
  if(x_right > imgDim) x_right = imgDim - 1; 
 
  /* Iterate through the values in vector Y, in integer points 
           [x_left,x_rigth]. */
  for(z=x_left; z <= x_right; z++) {
       
    xp = z; //positive x-coordinate
    xn = imgDim - 1 - xp; //negative x-coordinate
 
    // Look y from left.
    yp = imgDim/2 - floor(Yptr[xp + 1] + bin*shift_y) - 1;
    yn = imgDim - 1 - yp;
 
    // If the y-value for this x (z) is inside the image grid.
    if(yp < imgDim && yp >= 0 && yn < imgDim && yn >= 0 && xp >= 0 && 
             xp < imgDim && xn < imgDim && xn >= 0) {
 
      /* Check that pixel and pixels hit by symmetrical lines lay inside 
               the FOV. */ 
      if(ellipseIsInside(elli,yp,xp) && ellipseIsInside(elli,yn,xn) &&
         ellipseIsInside(elli,yp,xn) && ellipseIsInside(elli,yn,xp) &&
         ellipseIsInside(elli,xn,yn) && ellipseIsInside(elli,xp,yp) &&
         ellipseIsInside(elli,xn,yp) && ellipseIsInside(elli,xp,yn)){
 
        /* NOTE that pixels found in this part are always hit from right 
                 side. */
        // Compute a := |AF| and b := |FB|.
        a = (float)(xp + 1 - (imgDim/2 + Xptr[imgDim - yp] + bin*shift_x));
        b = (float)(floor(Yptr[xp + 1] + bin*shift_y) + 1 - 
                          (Yptr[xp + 1] + bin*shift_y));
 
        // Calculate the area of lower triangle.
        A = a*b/2;
 
        // c := |FC|
        c = (float)(xp + 1 - (imgDim/2 + Xptr[imgDim - yp] + 
                          (bin+1)*shift_x)); 
 
        if(c > 0) {
 
    // d := |FD|
    d = (float)(floor(Yptr[xp + 1] + (bin+1)*shift_y) + 1 - 
                            (Yptr[xp + 1] + (bin+1)*shift_y));
        
    // Subtract the area of upper triangle.
    A = A - c*d/2;
        }
        
        eps1 = A;
        
        if(eps1 < 0 || eps1 > 1) errors++;
 
        // Put coordinates and factors in the list.
        coptr[pix][0] = xp;
        coptr[pix][1] = yp;
        facptr[pix] = (unsigned short int)(scale*eps1);
        
        // Look for minimal and maximal factor.
        if(min>eps1) min=eps1;
        if(max<eps1) max=eps1;
        sqr_sum += eps1*eps1;
        sum += eps1;
        pix++;
      }// Is inside the FOV.
    }
  }
  
  // Limit (y-)indices for fast search.
  // Note that indices are non-negative integers.
  y_bottom = floor((float)(Yptr[imgDim] + bin*shift_y + imgDim/2));
  if(y_bottom < 0) y_bottom = 0;
  if(y_bottom > imgDim) y_bottom = 0; 
 
  y_top = floor((float)(Yptr[0] + bin*shift_y + imgDim/2));
  if(y_top > imgDim) y_top = imgDim;
  if(y_top <= 0) y_top = 1;
 
  /* Iterate through the values in vector X, in integer points 
           [y_bottom,y_top]. */
  for(z=y_top; z >= y_bottom; z--) {
   
    // Look y from this location.
    xp = floor(Xptr[z] + bin*shift_x) + imgDim/2;
    xn = imgDim - 1 - xp;
 
    yp = imgDim - z - 1;
    yn = imgDim - yp - 1;
 
    // If the x-value for this y (z) is inside the image grid.
    if(yp < imgDim && yp >= 0 && yn < imgDim && yn >= 0 && xp >= 0 && 
             xp < imgDim && xn < imgDim && xn >= 0) {
 
      /* Check that pixel and pixels hit by symmetrical lines lay inside 
               the FOV. */
      if(ellipseIsInside(elli,yp,xp) && ellipseIsInside(elli,yn,xn) &&
         ellipseIsInside(elli,yp,xn) && ellipseIsInside(elli,yn,xp) &&
         ellipseIsInside(elli,xn,yn) && ellipseIsInside(elli,xp,yp) &&
         ellipseIsInside(elli,xn,yp) && ellipseIsInside(elli,xp,yn)){
 
        eps1=eps2=eps3=0;
     
        dx = (float)(Xptr[z] + bin*shift_x + imgDim/2 - xp);
        dy = (float)(Yptr[xp] + bin*shift_y - z + imgDim/2); 
 
        if(dy < 1){ // Cases 3,4,5 and 6.
 
    // 1. PART
    // a := |HA|
    a = dy;
 
    // b := |HB| (b < 1 for angles in [0,pi/4))
    b = dx;
 
    // h := height of rectangle R.
    h = a + shift_y;
 
    if(h > 1){ // Cases 3,4 and 5.
      h = 1;
      g = b + shift_x;
      if(g > 1){ // Cases 3 and 4.
        g = 1;
        xp2 =floor(Xptr[z+1] + (bin+1)*shift_x) + imgDim/2;
 
        if(xp == xp2){ // Case 4.
          // c := |FC|
          c = (float)(xp + 1 - (imgDim/2 + Xptr[z+1] + 
                                 (bin+1)*shift_x));       
          // d := |FD| 
          d = (float)(floor(Yptr[xp + 1] + (bin+1)*shift_y) + 1 - 
                                 (Yptr[xp + 1] + (bin+1)*shift_y));
          eps1 = 1 - (a*b + c*d)/2;
          eps2 = 0;
          // Lower triangle on the next pixel (x+1,y).
          eps3 = (1 - d)*(b + shift_x - 1)/2; 
        } else { // Case 3.
          // c=d=0.
          eps1 = 1 - a*b/2;
          // Calculate area on pixel in place (xp+1,yp-1).
          dy = (float)(Yptr[xp+1] + (bin+1)*shift_y - (z + 1) + 
                                   imgDim/2);
          if(dy < 1){ // Case 11.
      // c := |HC|
      c = dy;
      // d := |HD|
      d = (float)(xp + 2 - (imgDim/2 + Xptr[z+1] + 
                                    (bin+1)*shift_x));
      eps2 = c*d/2;
          } else{ // Cases 9 and 10.
      dx = (float)(xp + 2 - (imgDim/2 + Xptr[z+1] + 
                                    (bin+1)*shift_x));
      if(dx < 1) { // Case 10.
        // g := length of rectangle Q.
        g = dx;
        // c := |CD| (on x-axis).
        c = dx - (float)(xp + 2 - (imgDim/2 + Xptr[z+2] + 
                                           (bin+1)*shift_x));
        // Rectangle Q - triangle c*h (h = 1).
        eps2 = g - c/2;
      } else { // Case 9.
        // c := |FC|
        c = (float)(xp + 2 - (imgDim/2 + Xptr[z+2] + 
                                      (bin+1)*shift_x));        
        // d := |FD| 
        d = (float)(floor(Yptr[xp + 2] + (bin+1)*shift_y) + 
                                      2 - (Yptr[xp + 2] + (bin+1)*shift_y));
        // Rectangle Q - triangle CFD.
        eps2 = 1 - c*d/2;
      }
          }
 
          // Calculate area on pixel in place (xp+1,yp).
          dx = (float)(xp + 2 - (imgDim/2 + Xptr[z] + 
                                  (bin+1)*shift_x));
          if(dx < 1) { // Case 10.
      // g := length of rectangle Q.
      g = dx;
      // c := |CD| (on x-axis).
      c = dx - (float)(xp + 2 - (imgDim/2 + Xptr[z+1] + 
                                        (bin+1)*shift_x));
      // Rectangle Q - triangle c*h (h = 1).
      eps3 = g - c/2;
          } else { // Case 9.
      // c := |FC|
      c = (float)(xp + 2 - (imgDim/2 + Xptr[z+1] + 
                                    (bin+1)*shift_x));        
      // d := |FD| 
      d = (float)(floor(Yptr[xp + 2] + (bin+1)*shift_y) + 
                                    1 - (Yptr[xp + 2] + (bin+1)*shift_y));
      // Rectangle Q - triangle CFD.
      eps3 = 1 - c*d/2;
          }
        }
      } else { // Case 5. (g < 1)
        // c := |DC|.
        c = g - (imgDim/2 + Xptr[z+1] + (bin+1)*shift_x - xp);
        // d := heigth
        d = 1;
 
        eps1 = g*h - (a*b + c*d)/2;
        eps2 = 0;
        eps3 = 0;
      }
    } else { // Case 6 (h <= 1). 
      // g := legth of rectangle R
      g = b + shift_x;
      if(g > 1) // Should always be < 1 for angles in [0,pi/4)
        g = 1; 
 
      eps1 = (g*h - a*b)/2;
      eps2 = 0;
      eps3 = 0;
    }
        } else { // 2. PART. Cases 2,7 and 8. (dy >= 1).
    
    // a := |HA|
    a = 1;
 
    // b := |HB| (>=1)
    b = dx - (float)(Xptr[z+1] + bin*shift_x + imgDim/2 - xp);
 
    // h := heigth of rectangle R
    h = 1;
        
    // g := length of rectangle R
    g = dx + shift_x;
 
    if(g > 1){ // Cases 2 and 8.
      g = 1 - (float)(Xptr[z+1] + bin*shift_x + imgDim/2 - xp);
      xp2 = floor(Xptr[z+1] + (bin+1)*shift_x) + imgDim/2;
      if(xp == xp2){ // Case 8.
        // c := |FC|
        c = (float)(xp + 1 - (imgDim/2 + Xptr[z+1] + 
                                (bin+1)*shift_x));        
        // d := |FD| 
        d = (float)(floor(Yptr[xp + 1] + (bin+1)*shift_y) + 1 - 
                                (Yptr[xp + 1] + (bin+1)*shift_y));
        eps1 = g*h - (a*b + c*d)/2;
        eps2 = 0;
        // Lower triangle on the next pixel (x+1,y).
        eps3 = (1 - d)*((Xptr[z] + (bin+1)*shift_x + imgDim/2) 
                                     - (xp+1))/2; 
      } else { // Case 2.
        // c=d=0.
        eps1 = g*h - a*b/2;
        /* Pixel in place (xp+1,yp-1) should have been found in 
                       the previous step. */
        eps2 = 0;
 
        // Calculate area on pixel in place (xp+1,yp).
        dx = (float)((imgDim/2 + Xptr[z] + (bin+1)*shift_x) - (xp+1));
        if(dx < 1){ // Case 10 (trapezium).
          // g := bottom of trapezium Q.
          g = dx;
          // c := top of trapezium Q.
          c = (float)((imgDim/2 + Xptr[z+1] + (bin+1)*shift_x) - 
                                  (xp+1));
          // Area of trapezium Q. Heigth = 1.
          eps3 = (g + c)/2;
        } else { // Case 9.
          // c := |FC|
          c = (float)(xp + 2 - (imgDim/2 + Xptr[z+1] + 
                                  (bin+1)*shift_x));        
          // d := |FD| 
          d = (float)(floor(Yptr[xp + 2] + (bin+1)*shift_y) + 1 - 
                                 (Yptr[xp + 2] + (bin+1)*shift_y));
          // Rectangle Q - triangle CFD.
          eps3 = 1 - c*d/2;
        }
      }
    } else { // Case 7. (g < = 1)
      // Area of the parallelogram R.
      eps1 = sdist/cosinus;
      eps2 = 0;
      eps3 = 0;
    }
        }
        if(eps1 <= 0 || eps1 > 1) {
    errors++;
    printf("Error: eps1 = %f \n",eps1);
    eps1 = fabs(eps1);
        }
        if(eps2 < 0 || eps2 > 1){
    errors++;
    eps2 = fabs(eps2);
        }
        if(eps3 < 0 || eps3 > 1){
    errors++;
    eps1 = fabs(eps3);
        }
        if(!eps2 && !eps3){ // Cases 5,6 and 7.
    // Put coordinates and factors in the list.
    coptr[pix][0] = xp;
    coptr[pix][1] = yp;
    facptr[pix] = (unsigned short int)(scale*eps1);
    pix++;
    // Look for minimal and maximal factor.
    if(min>eps1) min=eps1;
    if(max<eps1) max=eps1;
    sqr_sum += eps1*eps1;
    sum += eps1;
        } else{
    if(!eps2){ // <=> eps3 != 0 & eps2 = 0 <=> Cases 3,4 and 8.
      if(xp + 1 < imgDim && xn - 1 >= 0){
        // Put coordinates and factors in the list.
        coptr[pix][0] = xp;
        coptr[pix][1] = yp;
        facptr[pix] = (unsigned short int)(scale*eps1);
        
        // Look for minimal and maximal factor.
        if(min>eps1) min=eps1;
        if(max<eps1) max=eps1;
        sqr_sum += eps1*eps1;
        sum += eps1;
        pix++;
 
        // Put coordinates and factors in the list.
        coptr[pix][0] = xp+1;
        coptr[pix][1] = yp;
        facptr[pix] = (unsigned short int)(scale*eps3);
        
        // Look for minimal and maximal factor.
        if(min>eps3) min=eps3;
        if(max<eps3) max=eps3;
        sqr_sum += eps3*eps3;
        sum += eps3;      
        pix++;
      }
    } else{ // <=> eps2!=0 && eps3!=0 <=> Case 3.
      if(xp+1 < imgDim && xn-1 >= 0 && yp-1 >= 0 && yn+1 < imgDim) {
        // Put coordinates and factors in the list.
        coptr[pix][0] = xp;
        coptr[pix][1] = yp;
        facptr[pix] = (unsigned short int)(scale*eps1);
        
        // Look for minimal and maximal factor.
        if(min>eps1) min=eps1;
        if(max<eps1) max=eps1;
        sqr_sum += eps1*eps1;
        sum += eps1;
        pix++;
 
        // Put coordinates and factors in the list.
        coptr[pix][0] = xp+1;
        coptr[pix][1] = yp;
        facptr[pix] = (unsigned short int)(scale*eps3);
        
        // Look for minimal and maximal factor.
        if(min>eps3) min=eps3;
        if(max<eps3) max=eps3;
        sqr_sum += eps3*eps3;
        sum += eps3;
        pix++;
 
        // Put coordinates and factors in the list.
        coptr[pix][0] = xp+1;
        coptr[pix][1] = yp-1;
        facptr[pix] = (unsigned short int)(scale*eps2);
        
        // Look for minimal and maximal factor.
        if(min>eps2) min=eps2;
        if(max<eps2) max=eps2;
        sqr_sum += eps2*eps2;
        sum += eps2;      
        pix++;
      }
    }
        }
      
      }// Pixel is inside the FOV.      
    }
  }
  /* Allocate memory in factor pointer (dynamically according to number 
           of pixels intersected). */
  if(pix){
    mat->_factdata[half*view + bin]=
           (unsigned short int**)calloc(pix,sizeof(unsigned short int*)); 
    if(!mat->_factdata[half*view + bin]) return -1;
  }
 
  // Allocate leaves.
  for(col=0; col<pix; col++){
    mat->_factdata[half*view + bin][col]=
            (unsigned short int*)calloc(3,sizeof(unsigned short int)); 
    if(!mat->_factdata[half*view + bin][col]) return -1;
  }
 
  // Put now values in coordinates and factors array to result pointer.
  mat->fact = mat->_factdata;
  for(col=0; col<pix; col++){
    mat->fact[half*view + bin][col][0] = coptr[col][0]; // x-coodinate
    mat->fact[half*view + bin][col][1] = coptr[col][1];  // y-coordinate
    mat->fact[half*view + bin][col][2] = facptr[col];  // factor
  }
  
  //Set also the number of pixels belonging to each row and square sums.
  mat->dime[half*view + bin]=pix;
  mat->factor_sqr_sum[half*view + bin]=sqr_sum;
 
      }// END OF BIN-LOOP
  
    }// End of view>0.
   
  }// END OF VIEW LOOP
  
  free(X);
  free(Y);
  free(coords);
  free(factors);
 
  mat->min = min;
  mat->max = max;
  mat->factor_sum = sum;
  mat->status = PRMAT_STATUS_BS_OCCUPIED;
 
  if(RADON_VERBOSE) {
    printf("RADON: radonSetBasesEA() finished with %d errors.\n", errors);
    fflush(stdout);
  }
 
  return 0;
}// END OF SETTING BASE LINES (EXACT AREA).
/*****************************************************************************/
int radonSetLUT(RADON *radtra, ELLIPSE *elli, PRMAT *mat)
{
  unsigned int **tmpData, **tmpptr; // coordinates of lines response.
  unsigned int *coords, *iptr;
  unsigned int col, row, pix=0, p=0, lors=0, view, bin; //counters
  unsigned int binNr, viewNr, imgDim, views, half, center;
  int ret=0;
  int xp, xn, yp, yn;
  float diam, sampledist;
 
  //Check that projections are set.
  if(mat->status < PRMAT_STATUS_BS_OCCUPIED) return -1;
 
  //Retrieve data from given radon transform object
  imgDim=radonGetID(radtra);
  binNr=radonGetNB(radtra);
  viewNr=radonGetNV(radtra);
  sampledist=radonGetSD(radtra);
  half=radonGetHI(radtra);
  center=radonGetCB(radtra);
  diam=sqrt(2); 
 
  /* Allocate memory for temporary list where we store coordinates of lines of 
     response. Maximum number of lines (in one angle) hitting a pixel is 
     [diameter of a pixel]/[sample distance]. */
  lors = ceil(diam/sampledist) + 1;
  tmpData = (unsigned int**)calloc(imgDim*imgDim,sizeof(unsigned int*));
  for(p=0; p<imgDim*imgDim; p++){
    tmpData[p]=(unsigned int*)calloc(lors*viewNr,sizeof(unsigned int));
    if(!tmpData[p]){
      fprintf(stderr, "Error: not enough memory.\n");
      return(4); 
    }
  }
  tmpptr = tmpData;
 
  /* Initialize the first entry in the list, to keep the number of entries 
     in each row. And then fill the list from second entry. */
  for(p=0; p<imgDim*imgDim; p++) tmpptr[p][0]=0;
  views = viewNr/4;
  // Analyse the given projection matrix.
  for(view=0; view<views + 1; view++){
    for(bin=0; bin<half; bin++){
      row = view*half + bin;
      for(col=0; col<prmatGetPixels(mat,row); col++){
 
  xp = prmatGetXCoord(mat,row,col);
  xn = imgDim - 1 - xp;
  yp = prmatGetYCoord(mat,row,col);
  yn = imgDim - 1 - yp;
 
  if(xp != 0 || yp != 0){
    if(view == 0){
      /* Put coordinates of coincidence line in the next free place and 
               increase the counter. */
      tmpptr[yp*imgDim + xp][++tmpptr[yp*imgDim + xp][0]] = bin;
      if(bin != center)
        tmpptr[yp*imgDim+(imgDim-1-xp)][++tmpptr[yp*imgDim+(imgDim-1-xp)][0]]=
                binNr-bin-1;
      tmpptr[(imgDim-1-xp)*imgDim+yp][++tmpptr[(imgDim-1-xp)*imgDim+yp][0]]=
              binNr*(viewNr/2) + bin;
      if(bin != center)
        tmpptr[xp*imgDim+yp][++tmpptr[xp*imgDim+yp][0]] = 
                binNr*(viewNr/2) + (binNr-bin-1);
    }
 
    if(view == viewNr/4) {
      
      // Add img(x,y)*k to the raysum of LOR (view,bin)
      tmpptr[yp*imgDim + xp][++tmpptr[yp*imgDim + xp][0]] = 
              (viewNr/4)*binNr + bin;
      // Add img(-x,-y)*k to the raysum of LOR (view,binNr-bin)
      if(bin != center)
        tmpptr[yn*imgDim + xn][++tmpptr[yn*imgDim + xn][0]] = 
               (viewNr/4)*binNr + binNr - 1 - bin;
      
      // Rotate the LOR 90 degrees, i.e. x->-x.
      // Case: pi/2 + theta = 3pi/4.
      // Add img(-y,x)*k to the raysum of LOR (viewNr/2+view,bin)
      tmpptr[xn*imgDim + yp][++tmpptr[xn*imgDim + yp][0]] = 
              (viewNr/2 + viewNr/4)*binNr + bin;
      // Add img(y,-x)*k to the raysum of LOR (viewNr/2+view,binNr-bin)
      if(bin != center)
        tmpptr[xp*imgDim + yn][++tmpptr[xp*imgDim + yn][0]] = 
                (viewNr/2 + viewNr/4)*binNr + binNr - 1 - bin;
    }
 
    if(view != 0 && view != viewNr/4) {
      
      // Add img(x,y)*k to the raysum of LOR (view,bin)
      tmpptr[yp*imgDim + xp][++tmpptr[yp*imgDim + xp][0]] = 
              view*binNr + bin;
      if(bin != center)
        // Add img(-x,-y)*k to the raysum of LOR (view,binNr-bin)
        tmpptr[yn*imgDim + xn][++tmpptr[yn*imgDim + xn][0]] = 
                view*binNr + binNr - 1 - bin;
    
      // Mirror the LOR on line x=y, and on y-axis i.e x->y and y->-y.
      // Case: pi/2+theta
      // Add img(y,-x)*k to the raysum of LOR (viewNr/2+view,bin)
      tmpptr[xn*imgDim + yp][++tmpptr[xn*imgDim + yp][0]] = 
              (viewNr/2 + view)*binNr + bin;
      // Add img(y,-x)*k to the raysum of LOR (viewNr/2+view,binNr-bin)
      if(bin != center)
        tmpptr[xp*imgDim + yn][++tmpptr[xp*imgDim + yn][0]] =
                (viewNr/2 + view)*binNr + binNr - 1 - bin;
          
      // Mirror the original LOR on y-axis, i.e. x->-x
      // Case: pi-theta.
      // Add img(-x,y)*k to the raysum of LOR (viewNr-view,bin)
      tmpptr[yp*imgDim + xn][++tmpptr[yp*imgDim + xn][0]] =
              (viewNr - view)*binNr + bin;
      // Add img(x,-y)*k to the raysum of LOR (viewNr-view,binNr-bin)
      if(bin != center)
        tmpptr[yn*imgDim + xp][++tmpptr[yn*imgDim + xp][0]] =
                (viewNr-view)*binNr + binNr - 1 - bin;
        
      // Mirror the LOR on line x=y, i.e. x->y.
      // Case: pi/2-theta
      // Add img(-y,-x)*k to the raysum of LOR (viewNr/2-view,bin)
      tmpptr[xn*imgDim + yn][++tmpptr[xn*imgDim + yn][0]] =
              (viewNr/2 - view)*binNr + bin;
      // Add img(y,x)*k to the raysum of LOR (viewNr/2-view,binNr-bin)
      if(bin != center)
        tmpptr[xp*imgDim + yp][++tmpptr[xp*imgDim + yp][0]] = 
                (viewNr/2 - view)*binNr + binNr - 1 - bin;  
    }
  }
      }
    }
  }
 
  pix = 0;
  // Get the number of pixels inside the FOV.
  for(row = 0; row < imgDim; row++){
    for(col = 0; col < imgDim; col++){
      if(ellipseIsInside(elli,row,col))
  pix++;
    }
  }
 
  // Analyse the tmp data array.
  p = 0; // index of a pixel inside the FOV
  coords = (unsigned int*)calloc(pix,sizeof(int));
  iptr = coords;
  for(row = 0; row < imgDim; row++){
    for(col = 0; col < imgDim; col++){
      if(ellipseIsInside(elli,row,col)){
  // Put the number of coincidence lines hitting this pixel into the list.
  iptr[p] = tmpptr[row*imgDim + col][0];
  p++; // increase pixel counter
      }
    }
  }
 
  // Allocate memory for the look-up table.
  ret = prmatAllocate(mat,0,pix,coords);
  if(ret){
    free((unsigned int**)tmpData);
    free((int*)coords);
    return(ret);
  }
 
  /* Put pixel coordinates, number of lines and the coordinates of 
     the coincidence lines in the structure. */
  p = 0;
  tmpptr = tmpData;
  iptr = coords;
  for(row = 0; row < imgDim; row++) {
    for(col = 0; col < imgDim; col++) {
      if(ellipseIsInside(elli,row,col)) {
  // Put pixel coordinates in the list.
  mat -> lines[p][0] = row*imgDim + col;
  // Put the number of lines hitting this pixel in the list.
  mat -> lines[p][1] = iptr[p];
  // Put the coordinates of the coincidence lines in the list.
  for(lors = 0; lors < iptr[p]; lors++)
    mat -> lines[p][lors + 2] = tmpptr[row*imgDim + col][lors+1];
 
  p++; // increase pixel counter
      }
    }
  }
 
  printf("\n");
  mat->status = PRMAT_STATUS_LU_OCCUPIED;
 
  free((unsigned int**)tmpData);
  free((int*)coords);
 
  return 0;
}// END OF SETTING LOOK-UP TABLE.
/*****************************************************************************/
int radonSetLORS(RADON *radtra, ELLIPSE *elli, PRMAT *mat)
{
  unsigned int *coords, *iptr;
  unsigned int imgDim=128, binNr=256, viewNr=192, view, bin, pixNr=0;
  unsigned int row, col, rows, half, centerbin, p;
  int ret=0;
  int xp, xn, yp, yn;
  float fact, sqr_sum;
  PRMAT ext_mat; // temporary extented projection matrix
 
  if(elli) {} // to prevent compiler warning about not being used.
 
  printf("radonSetLors() started. \n");
 
  //Retrieve data from given radon transform object
  imgDim=radonGetID(radtra);
  binNr=radonGetNB(radtra);
  viewNr=radonGetNV(radtra);
  
  // Calculate and set center bin for the current geometrics.
  if((binNr%2) != 0){
      half = (binNr - 1)/2 + 1;
      centerbin = half - 1;
  }
  else{
    half = binNr/2;
    // In the case binNr is even there is no center bin.
    centerbin = -1;
  }
 
  // Initiate extented projection matrix.
  prmatInit(&ext_mat);
  
  // Prepare to allocate extented projection matrix.
  rows = viewNr*binNr;
  coords = (unsigned int*)calloc(rows,sizeof(int));
  iptr = coords;
  // Determine the number of hit pixels for EVERY row.
  p = viewNr/4;
  // Extend the projection matrix according to given base set.
  for(view=0; view<p + 1; view++){
    for(bin=0; bin<half; bin++){
      row = view*half + bin;
      pixNr = prmatGetPixels(mat,row);
      // Line in the base set.
      iptr[view*binNr + bin] = pixNr;
      // Symmetrical lines.
      if(bin != centerbin)
  iptr[view*binNr + binNr - 1 - bin] = pixNr;
 
      if(view != 0 && view != viewNr/4){        
  // Mirror the original LOR on y-axis, i.e. x->-x
  // Case: pi-theta.
  iptr[(viewNr - view)*binNr + bin] = pixNr;
  if(bin != centerbin)
    iptr[(viewNr-view)*binNr + binNr - 1 - bin] = pixNr;
        
  // Mirror the LOR on line x=y, i.e. x->y.
  // Case: pi/2-theta
  iptr[(viewNr/2 - view)*binNr + bin] = pixNr;
  if(bin != centerbin)
    iptr[(viewNr/2 - view)*binNr + binNr - 1 - bin] = pixNr;
      }
 
      // Mirror the LOR on line x=y, and on y-axis i.e x->y and y->-y.
      // Case: pi/2+theta
      iptr[(viewNr/2 + view)*binNr + bin] = pixNr;
      if(bin != centerbin)
  iptr[(viewNr/2 + view)*binNr + binNr - 1 - bin] = pixNr;
      
    }
  }
 
  // Allocate the extented projection matrix.
  ret = prmatAllocate(&ext_mat,1,rows,coords);
  if(ret){
    free(coords);
    return(ret);
  }
  printf("Allocation done \n");
    
  // Extend the projection matrix according to given base set.
  for(view=0; view<p + 1; view++){
    for(bin=0; bin<half; bin++){
      row = view*half + bin;
      for(col=0; col<prmatGetPixels(mat,row); col++){
 
  sqr_sum = prmatGetFactorSqrSum(mat,row);
  // Notice that we want the factor in unsigned short int.
  fact = mat->fact[row][col][2];
  xp = prmatGetXCoord(mat,row,col);
  xn = imgDim - 1 - xp;
  yp = prmatGetYCoord(mat,row,col);
  yn = imgDim - 1 - yp;
      
  // Line in the base set.
  ext_mat.fact[view*binNr + bin][col][0] = xp;
  ext_mat.fact[view*binNr + bin][col][1] = yp;
  ext_mat.fact[view*binNr + bin][col][2] = fact;
  ext_mat.factor_sqr_sum[view*binNr + bin] = sqr_sum;
 
  // Then symmetrical lines.
  if(bin != centerbin){
    ext_mat.fact[view*binNr + binNr - 1 - bin][col][0] = xn;
    ext_mat.fact[view*binNr + binNr - 1 - bin][col][1] = yn;
    ext_mat.fact[view*binNr + binNr - 1 - bin][col][2] = fact;
    ext_mat.factor_sqr_sum[view*binNr + binNr - 1 - bin] = sqr_sum;
  }
 
  if(view != 0 && view != viewNr/4){        
    // Mirror the original LOR on y-axis, i.e. x->-x
    // Case: pi-theta.
    ext_mat.fact[(viewNr - view)*binNr + bin][col][0] = xn;
    ext_mat.fact[(viewNr - view)*binNr + bin][col][1] = yp;
    ext_mat.fact[(viewNr - view)*binNr + bin][col][2] = fact;
    ext_mat.factor_sqr_sum[(viewNr - view)*binNr + bin] = sqr_sum;
 
    if(bin != centerbin){
 
      ext_mat.fact[(viewNr - view)*binNr +binNr - 1 - bin][col][0] = xp;
      ext_mat.fact[(viewNr - view)*binNr +binNr - 1 - bin][col][1] = yn;
      ext_mat.fact[(viewNr - view)*binNr +binNr - 1 - bin][col][2] = fact;
      ext_mat.factor_sqr_sum[(viewNr - view)*binNr +binNr - 1 - bin] = 
              sqr_sum;
    }
        
    // Mirror the LOR on line x=y, i.e. x->y.
    // Case: pi/2-theta
    ext_mat.fact[(viewNr/2 - view)*binNr + bin][col][0] = yn;
    ext_mat.fact[(viewNr/2 - view)*binNr + bin][col][1] = xn;
    ext_mat.fact[(viewNr/2 - view)*binNr + bin][col][2] = fact;
    ext_mat.factor_sqr_sum[(viewNr/2 - view)*binNr + bin] = sqr_sum;
 
    if(bin != centerbin){
      ext_mat.fact[(viewNr/2 - view)*binNr +binNr-1-bin][col][0] = yp;
      ext_mat.fact[(viewNr/2 - view)*binNr +binNr-1-bin][col][1] = xp;
      ext_mat.fact[(viewNr/2 - view)*binNr +binNr-1-bin][col][2] = fact;
      ext_mat.factor_sqr_sum[(viewNr/2-view)*binNr+binNr-1-bin] = sqr_sum;
    }     
  }
 
  // Mirror the LOR on line x=y, and on y-axis i.e x->y and y->-y.
  // Case: pi/2+theta
  ext_mat.fact[(viewNr/2 + view)*binNr + bin][col][0] = yp;
  ext_mat.fact[(viewNr/2 + view)*binNr + bin][col][1] = xn;
  ext_mat.fact[(viewNr/2 + view)*binNr + bin][col][2] = fact;
  ext_mat.factor_sqr_sum[(viewNr/2 + view)*binNr + bin] = sqr_sum;
  
  // Add img(y,-x)*k to the raysum of LOR (viewNr/2+view,binNr-bin)
  if(bin != centerbin){
    ext_mat.fact[(viewNr/2 + view)*binNr + binNr-1-bin][col][0] = yn;
    ext_mat.fact[(viewNr/2 + view)*binNr + binNr-1-bin][col][1] = xp;
    ext_mat.fact[(viewNr/2 + view)*binNr + binNr-1-bin][col][2] = fact;
    ext_mat.factor_sqr_sum[(viewNr/2 +view)*binNr +binNr-1-bin] = sqr_sum;
  }
      } 
    }
  }
 
 
  // Empty old data from the given projection matrix.
  free((float*)mat->factor_sqr_sum);
 
  for(row=0; row<mat->dimr; row++){ //every row
    for(col=0; col<mat->dime[row]; col++){ //every factor
      free((unsigned short int*)mat->_factdata[row][col]);
    }
  }
  
  free((int*)mat->dime);
 
  if(mat->dimr>0) free(mat->_factdata); 
  mat->dimr=0; 
 
  // Allocate memory in mat structure for new extented matrix.
  ret = prmatAllocate(mat, 1, rows, coords);
  if(ret){
    free(coords);
    prmatEmpty(&ext_mat);
    return(ret);
  }
  printf("Allocation done \n");
 
  // Copy the matrix from the temporary projection matrix.
  for(row=0; row<prmatGetRows(&ext_mat); row++){
    // Square sums.
    mat->factor_sqr_sum[row] = prmatGetFactorSqrSum(&ext_mat,row);
    for(col=0; col<prmatGetPixels(&ext_mat,row); col++){
      // Coordinates and factors.
      mat->fact[row][col][0] = prmatGetXCoord(&ext_mat,row,col);
      mat->fact[row][col][1] = prmatGetYCoord(&ext_mat,row,col);
      mat->fact[row][col][2] = ext_mat.fact[row][col][2];
    }
  }
 
  // Release the temporary projection matrix.
  prmatEmpty(&ext_mat);
  free(coords);
 
  return 0;
 
 
}//END OF radonSetLORS
/*****************************************************************************/
 
/*****************************************************************************/