radon.c File Reference

Radon transform. More...

#include "libtpcrec.h"

Go to the source code of this file.

Functions
void	radonEmpty (RADON *radtra)
int	radonSet (RADON *radtra, int mode, int imgDim, int viewNr, int binNr)
int	radonGetMO (RADON *radtra)
int	radonGetID (RADON *radtra)
int	radonGetNV (RADON *radtra)
int	radonGetNB (RADON *radtra)
float	radonGetSD (RADON *radtra)
int	radonGetHI (RADON *radtra)
int	radonGetCB (RADON *radtra)
float	radonGetSin (RADON *radtra, int nr)
int	radonFwdTransform (RADON *radtra, int set, int setNr, float *imgdata, float *scndata)
int	radonFwdTransformEA (RADON *radtra, int set, int setNr, float *imgdata, float *scndata)
int	radonFwdTransformSA (RADON *radtra, int set, int setNr, float *imgdata, float *scndata)
int	radonFwdTransformPRM (PRMAT *mat, int set, int setNr, float *imgdata, float *scndata)
int	radonBackTransform (RADON *radtra, int set, int setNr, float *scndata, float *imgdata)
int	radonBackTransformEA (RADON *radtra, int set, int setNr, float *scndata, float *imgdata)
int	radonBackTransformSA (RADON *radtra, int set, int setNr, float *scndata, float *imgdata)
int	radonBackTransformPRM (PRMAT *mat, int set, int setNr, float *scndata, float *imgdata)
int	radonSetBases (RADON *radtra, ELLIPSE *elli, PRMAT *mat)
int	radonSetBasesEA (RADON *radtra, ELLIPSE *elli, PRMAT *mat)
int	radonSetLUT (RADON *radtra, ELLIPSE *elli, PRMAT *mat)
int	radonSetLORS (RADON *radtra, ELLIPSE *elli, PRMAT *mat)

Variables
int	RADON_TEST
	Drive in test mode if not 0.
int	RADON_VERBOSE
	Drive in verbose mode if not 0.

Detailed Description

Radon transform.

Radon data structure contains the parameters defining a Radon transform. Transformations to and from the Radon domain are implemented in this file. Discretisation of a continuos Radon transform has five (5) different implementations in this file.

Author

Jarkko Johansson

Date

2006-06-16

Definition in file radon.c.

Function Documentation

radonBackTransform()

int radonBackTransform	(	RADON *	radtra,
		int	set,
		int	setNr,
		float *	scndata,
		float *	imgdata )

Transforms the given sinogram in Radon domain to spatial domain. Parameters of the discrete Radon transform are stored in the RADON object. Transform is calculated only in those angles belonging into given subset. Subset contains angles starting from the index 'set' with spacing 'setNr'. Discretisation model utilised in this function is '0/1' or 'length of intersection' according to the given RADON object.

Precondition

scndata contains the projections in v x b lines of response.

Postcondition

scndata is mapped into the cartesian space.

Parameters

radtra	contains an initialized radon transform object.
set	tells which set is to be utilized; 0 < set < setNr-1
setNr	number of subsets
scndata	v x b vector contains the projections.
imgdata	n x n vector for storing the back-projection.

Returns

int 0 if ok

Definition at line 1428 of file radon.c.

  {
  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.

Referenced by radonBackTransform().

radonBackTransformEA()

int radonBackTransformEA	(	RADON *	radtra,
		int	set,
		int	setNr,
		float *	scndata,
		float *	imgdata )

Same as 'radonBackTransform()' but discretisation model in this function is 'exact area'.

Parameters

radtra	contains an initialized radon transform object.
set	tells which set is to be utilized; 0 < set < setNr-1
setNr	number of subsets
scndata	v x b vector contains the projections.
imgdata	n x n vector for storing the back-projection.

Returns

int 0 if ok

Definition at line 1695 of file radon.c.

  {
  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.

Referenced by radonBackTransformEA().

radonBackTransformPRM()

int radonBackTransformPRM	(	PRMAT *	mat,
		int	set,
		int	setNr,
		float *	scndata,
		float *	imgdata )

Transforms the given sinogram in Radon domain to spatial domain. Transform is calculated by multiplying with the transpose of the given projection matrix from the right. Transform is calculated only in those angles belonging into given subset. Subset contains angles starting from the index 'set' with spacing 'setNr'. Discretisation model utilised in this function is '0/1', 'length of intersection' or 'exact area' according to the given projection matrix.

Precondition

scndata contains the projections with v x b lines of response.

Postcondition

scndata is mapped into the cartesian space.

Parameters

mat	pointer to structure where coordinates and values are stored.
set	tells which set is to be utilized; 0 < set < setNr-1
setNr	number of subsets
scndata	v x b vector contains the projections.
imgdata	n x n vector for storing the back-projection.

Returns

int 0 if ok

Definition at line 2507 of file radon.c.

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

Referenced by radonBackTransformPRM().

radonBackTransformSA()

int radonBackTransformSA	(	RADON *	radtra,
		int	set,
		int	setNr,
		float *	scndata,
		float *	imgdata )

Same as 'radonBackTransform()' but discretisation model in this function is 'linear interpolation' or 'nearest neighbour interpolation' according to the given Radon transform object.

Parameters

radtra	contains an initialized radon transform object.
set	tells which set is to be utilized; 0 < set < setNr-1
setNr	number of subsets
scndata	v x b vector contains the projections.
imgdata	n x n vector for storing the back-projection.

Returns

int 0 if ok

Note

first written by Sakari Alenius 1998.

Definition at line 2419 of file radon.c.

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

Referenced by radonBackTransformSA().

radonEmpty()

void radonEmpty

(

RADON *

radtra

)

Frees the memory allocated for radon transform. All data is cleared.

Postcondition

radon transform is emptied.

Parameters

radtra

pointer to transform data to be emptied

Definition at line 27 of file radon.c.

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

Referenced by radonEmpty().

radonFwdTransform()

int radonFwdTransform	(	RADON *	radtra,
		int	set,
		int	setNr,
		float *	imgdata,
		float *	scndata )

Performs the perpendicular model of the radon transform for the given vector in spatial domain. The discretisation mode is chosen according to transform object, discretisation mode is 0 or 1. Transform is performed in those angles belonging in the given subset. The set of coincidence lines is divided into subsets in the following way: let i be the number of subsets, then every ith angle in the base set and the ones symmetrical with it belong to the same subset.

Precondition

imgdata contains pixel values for every pixel in dim*dim image grid

Postcondition

imgdata is mapped into the projection space.

Parameters

radtra	contains an initialized radon transform object.
set	tells which set is to be utilized; 0 < set < setNr-1
setNr	number of subsets
imgdata	contains pixel values for every pixel in dim*dim image grid.
scndata	viewNr*binNr vector for storing the projection

Returns

int 0 if ok

Definition at line 216 of file radon.c.

  {
  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

Referenced by radonFwdTransform().

radonFwdTransformEA()

int radonFwdTransformEA	(	RADON *	radtra,
		int	set,
		int	setNr,
		float *	imgdata,
		float *	scndata )

Same as 'radonFwdTransform()' but discretisation model in this function is 'exact area'.

Parameters

radtra	contains an initialized radon transform object.
set	tells which set is to be utilized; 0 < set < setNr-1
setNr	number of subsets
imgdata	contains pixel values for every pixel in dim*dim image grid.
scndata	viewNr*binNr vector for storing the projection

Returns

int 0 if ok

Definition at line 549 of file radon.c.

{
  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

Referenced by radonFwdTransformEA().

radonFwdTransformPRM()

int radonFwdTransformPRM	(	PRMAT *	mat,
		int	set,
		int	setNr,
		float *	imgdata,
		float *	scndata )

Transforms the given intensity image in spatial domain to Radon domain. Transform is calculated by multiplying with the given projection matrix from the left. Transform is calculated only in those angles belonging into given subset. Subset contains angles starting from the index 'set' with spacing 'setNr'. Discretisation model utilised in this function is '0/1', 'length of intersection' or 'exact area' according to the given projection matrix.

Precondition

imgdata contains pixel values for every pixel in n x n image grid

Postcondition

imgdata is mapped into the projection space.

Parameters

mat	contains an initialised projection matrix
set	tells which set is to be utilized; 0 < set < setNr-1
setNr	number of subsets (spacing between indices)
imgdata	contains pixel values for every pixel in dim*dim image grid.
scndata	viewNr*binNr vector for storing the projection

Returns

int 0 if ok

Definition at line 1327 of file radon.c.

  {
  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

Referenced by radonFwdTransformPRM().

radonFwdTransformSA()

int radonFwdTransformSA	(	RADON *	radtra,
		int	set,
		int	setNr,
		float *	imgdata,
		float *	scndata )

Same as 'radonFwdTransform()' but discretisation model in this function is 'linear interpolation' or 'nearest neighbour interpolation' according to the given Radon transform object.

Parameters

radtra	contains an initialized radon transform object.
set	tells which set is to be utilized; 0 < set < setNr-1
setNr	number of subsets
imgdata	contains pixel values for every pixel in dim*dim image grid.
scndata	viewNr*binNr vector for storing the projection

Returns

int 0 if ok

Note

First written by Sakari Alenius 1998.

Definition at line 1237 of file radon.c.

  {
  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

Referenced by radonFwdTransformSA().

radonGetCB()

int radonGetCB

(

RADON *

radtra

)

Returns the center bin in this radon transform.

Parameters

radtra

radon transform for which the center bin is to be returned.

Returns

int index of the center bin.

Definition at line 180 of file radon.c.

{
  return radtra->centerBin;
}

Referenced by radonBackTransform(), radonBackTransformEA(), radonFwdTransform(), radonFwdTransformEA(), radonGetCB(), and radonSetLUT().

radonGetHI()

int radonGetHI

(

RADON *

radtra

)

Returns the half index of the bins in this radon transform.

Parameters

radtra

radon transform for which the half index is to be returned.

Returns

int half index of the bins.

Definition at line 169 of file radon.c.

{
  return radtra->half;
}

Referenced by radonBackTransform(), radonBackTransformEA(), radonFwdTransform(), radonFwdTransformEA(), radonGetHI(), radonSetBases(), radonSetBasesEA(), and radonSetLUT().

radonGetID()

int radonGetID

(

RADON *

radtra

)

Returns the image dimension in this radon transform.

Parameters

radtra

radon transform for which the image dimension is to be returned.

Returns

int image dimension.

Definition at line 128 of file radon.c.

{
  return radtra->imgDim;
}

Referenced by radonBackTransform(), radonBackTransformEA(), radonBackTransformSA(), radonFwdTransform(), radonFwdTransformEA(), radonFwdTransformSA(), radonGetID(), radonSetBases(), radonSetBasesEA(), radonSetLORS(), and radonSetLUT().

radonGetMO()

int radonGetMO

(

RADON *

radtra

)

Returns the discretization model of this radon transform.

Parameters

radtra

radon transform for which the mode is to be returned.

Returns

int mode

Definition at line 119 of file radon.c.

{
  return radtra->mode;
}

Referenced by radonBackTransform(), radonBackTransformSA(), radonFwdTransform(), radonFwdTransformSA(), radonGetMO(), radonSetBases(), and radonSetBasesEA().

radonGetNB()

int radonGetNB

(

RADON *

radtra

)

Returns the number of bins in this radon transform.

Parameters

radtra

radon transform for which the number of bins is to be returned.

Returns

int number of bins.

Definition at line 148 of file radon.c.

{
  return radtra->binNr;
}

Referenced by radonBackTransform(), radonBackTransformEA(), radonBackTransformSA(), radonFwdTransform(), radonFwdTransformEA(), radonFwdTransformSA(), radonGetNB(), radonSetBases(), radonSetBasesEA(), radonSetLORS(), and radonSetLUT().

radonGetNV()

int radonGetNV

(

RADON *

radtra

)

Returns the number of views in this radon transform.

Parameters

radtra

radon transform for which the number of views is to be returned.

Returns

int number of views.

Definition at line 138 of file radon.c.

{
  return radtra->viewNr;
}

Referenced by radonBackTransform(), radonBackTransformEA(), radonBackTransformSA(), radonFwdTransform(), radonFwdTransformEA(), radonFwdTransformSA(), radonGetNV(), radonSetBases(), radonSetBasesEA(), radonSetLORS(), and radonSetLUT().

radonGetSD()

float radonGetSD

(

RADON *

radtra

)

Returns the sample distance in this radon transform.

Parameters

radtra

radon transform for which the sample distance is to be returned.

Returns

float sample distance.

Definition at line 159 of file radon.c.

{
  return radtra->sampleDist;
}

Referenced by radonBackTransform(), radonBackTransformEA(), radonFwdTransform(), radonFwdTransformEA(), radonGetSD(), radonSetBases(), radonSetBasesEA(), and radonSetLUT().

radonGetSin()

float radonGetSin	(	RADON *	radtra,
		int	nr )

Returns the sine for given angle (index).

Parameters

radtra	radon transform for which the sine is to be returned.
nr	index of the angle to be returned (angle=nr*pi/radonGetNB(RADON)).

Returns

float sin(nr*pi/binNr).

Definition at line 192 of file radon.c.

{
  return radtra->sines[nr];
}

Referenced by radonBackTransform(), radonBackTransformEA(), radonBackTransformSA(), radonFwdTransform(), radonFwdTransformEA(), radonFwdTransformSA(), radonGetSin(), radonSetBases(), and radonSetBasesEA().

radonSet()

int radonSet	(	RADON *	radtra,
		int	mode,
		int	imgDim,
		int	viewNr,
		int	binNr )

Sets the data for the 2-D Radon transform.

Precondition

dim > 0

Postcondition

Parameters for the radon transform are set

Parameters

radtra	radon transform for which the paramaters are to be set
mode	discretisation mode
imgDim	image dimension
viewNr	Number of views (angles)
binNr	Number of bins (distances)

Returns

0 if ok

Definition at line 62 of file radon.c.

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

Referenced by radonSet().

radonSetBases()

int radonSetBases	(	RADON *	radtra,
		ELLIPSE *	elli,
		PRMAT *	mat )

Sets the coordinates and factors for intersected pixels according to given special Radon transform operator for every coincidence line in the BASE set. Base set includes coincidence lines in range [0,pi/4].

Precondition

radtra is a radon transform operator && elli defines a field of view && mat is initialized.

Postcondition

coordinates and factors for pixels contributing to coincidence lines in the base set.

Parameters

radtra	special Radon transform operator.
elli	field of view.
mat	pointer to the datastructure where coordinates and values are to be stored.

Returns

0 if ok.

Definition at line 2606 of file radon.c.

  {
  // 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.

Referenced by radonSetBases().

radonSetBasesEA()

int radonSetBasesEA	(	RADON *	radtra,
		ELLIPSE *	elli,
		PRMAT *	mat )

Same as 'radonSetBases()' but discretisation model in this function is 'exact area'.

Parameters

radtra	special Radon transform operator.
elli	field of view.
mat	pointer to the datastructure where coordinates and values are to be stored.

Returns

0 if ok.

Definition at line 2988 of file radon.c.

{
  // 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).

Referenced by radonSetBasesEA().

radonSetLORS()

int radonSetLORS	(	RADON *	radtra,
		ELLIPSE *	elli,
		PRMAT *	mat )

Sets the coordinates and factors for intersected pixels according to given special Radon transform operator for EVERY line of response.

Precondition

radtra is a radon transform operator && elli defines a field of view && mat is initialized.

Postcondition

coordinates and factors for pixels contributing to all lines of response are set.

Parameters

radtra	special Radon transform operator.
elli	field of view.
mat	pointer to the datastructure where coordinates and values are stored for the base lines.

Returns

0 if ok.

Definition at line 3934 of file radon.c.

{
  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

Referenced by radonSetLORS().

radonSetLUT()

int radonSetLUT	(	RADON *	radtra,
		ELLIPSE *	elli,
		PRMAT *	mat )

Sets a look-up table containing coordinates of lines of response intersecting pixels inside a field of view. Lines of response contributing to a pixel are searched from already set projection matrix.

Precondition

radtra is a radon transform operator && elli defines a field of view && projections are set in structure mat.

Postcondition

coordinates of lines response intersecting a pixel are set in a table.

Parameters

radtra	special Radon transform operator.
elli	field of view.
mat	pointer to the datastructure where coordinates and values are to be stored.

Returns

0 if ok.

Definition at line 3728 of file radon.c.

{
  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.

Referenced by radonSetLUT().

Variable Documentation

RADON_TEST

int RADON_TEST

Drive in test mode if not 0.

Definition at line 12 of file radon.c.

RADON_VERBOSE

int RADON_VERBOSE

Drive in verbose mode if not 0.

Definition at line 13 of file radon.c.

Referenced by radonBackTransformEA(), radonEmpty(), radonFwdTransform(), radonFwdTransformEA(), radonSet(), and radonSetBasesEA().