recutil.c File Reference

General routines for image reconstruction procedures. More...

#include "libtpcrec.h"

Go to the source code of this file.

Functions
void	recSinTables (int views, float *sinB, float *sinBrot, float rotation)
void	recInterpolateSinogram (float *srcsino, float *newsino, int srcrays, int newrays, int views)
int	bit_rev_int (int x, int n)
void	set_os_set (int os_sets, int *set_seq)
int	recGetStatistics (float *buf, int n, float *osum, float *omin, float *omax, int skip_zero_mins)

Detailed Description

General routines for image reconstruction procedures.

Based on the code written by Sakari Alenius.

Author

Vesa Oikonen

Definition in file recutil.c.

Function Documentation

bit_rev_int()

int bit_rev_int	(	int	x,
		int	n )

Bit-reverse for integers in a list, used by set_os_set().

See also

set_os_set()

Returns

the converted integer.

Parameters

x	Integer to bit-convert.
n	Number of elements.

Definition at line 77 of file recutil.c.

  {
  unsigned char c=(char)x;
  switch(n) {
    case 4:
      c = (c&1)<<1 | (c&2)>>1;
      break;
    case 8:
      c = ((c&1) << 2) | (c&2) | ((c&4) >> 2);
      break;
    case 16:
      c = ((c&1) << 3) | ((c&2) << 1) | ((c&4) >> 1) | ((c&8) >> 3);
      break;
    case 32:
      c = ((c&1)<<4) | ((c&2)<<2) | (c&4) | ((c&8)>>2) | ((c&16)>>4);
      break;
    case 64:
      c = ((c&1)<<5) | ((c&2)<<3) | ((c&4)<<1) | ((c&8)>>1) | ((c&16)>>3) |
          ((c&32)>>5);
      break;
    case 128:
      c = ((c&1)<<6) | ((c&2)<<4) | ((c&4)<<2) | (c&8) | ((c&16)>>2) |
          ((c&32)>>4) | ((c&64)>>6);
      break;
    default: break;
  }
  return((int)c);
}

Referenced by set_os_set().

recGetStatistics()

int recGetStatistics	(	float *	buf,
		int	n,
		float *	osum,
		float *	omin,
		float *	omax,
		int	skip_zero_mins )

Get the sum and minimum and maximum values from a list of floats, optionally ignoring zero values from the minimum.

Returns

Returns the number of non-zero values.

Parameters

buf	Pointer to the float array of length n.
n	Array size.
osum	Sum of array values is returned here; enter NULL if not needed.
omin	Minimum value is returned here; enter NULL if not needed.
omax	Maximum value is returned here; enter NULL if not needed.
skip_zero_mins	Skip zero values when determining the minimum; 1 or 0.

Definition at line 154 of file recutil.c.

  {
  int i, nonzeroes=0;
  float sum, min, max;
 
  if(osum!=NULL) *osum=nanf("");
  if(omin!=NULL) *omin=nanf("");
  if(omax!=NULL) *omax=nanf("");
 
  i=0; while(i<n && !isfinite(buf[i])) i++;
  if(i==n) return(0);
 
  sum=0.0; max=buf[i]; min=nanf(""); 
  if(!skip_zero_mins || min!=0.0) min=buf[i];
  for(; i<n; i++) {
    if(!isfinite(buf[i])) continue;
    sum+=buf[i];
    if(buf[i]>max) max=buf[i];
    if(buf[i]==0.0) {
      if(skip_zero_mins) continue;
    } else {
      nonzeroes++;
    }
    if(isnan(min) || buf[i]<min) min=buf[i];
  }
  if(osum!=NULL) *osum=sum;
  if(omin!=NULL) *omin=min;
  if(omax!=NULL) *omax=max;
 
  return(nonzeroes);
}

Referenced by mrp().

recInterpolateSinogram()

void recInterpolateSinogram	(	float *	srcsino,
		float *	newsino,
		int	srcrays,
		int	newrays,
		int	views )

Interpolate sinogram so that the sinogram bin width is the same as the image pixel width.

Parameters

srcsino	Source sinogram data, size of srcrays*views.
newsino	New interpolated sinogram data calculated here; must be allocated with size newrays*views.
srcrays	Number of rays (bins, columns) in the original sinogram.
newrays	Number of rays (bin, columns) in the new interpolated sinogram.
views	Number of projection views (angles, rows) in both sinograms.

Definition at line 37 of file recutil.c.

  {
  if(srcrays==newrays) {
    /* no interpolation needed; just copy and return */
    for(int i=0; i<views*newrays; i++) newsino[i]=srcsino[i];
    return;
  }
 
  float diff, *sp, *np, oposition, weight;
  np=newsino;
  for(int j=0; j<views; j++) {
    sp=srcsino + j*srcrays;
    for(int i=0; i<newrays; i++) {
      oposition=(float)(i*srcrays)/(float)newrays;
      weight = oposition - (float)(int)oposition;
      sp=srcsino + j*srcrays + (int)oposition;
      diff = *(sp+1) - *sp;
      *np++ = *sp + diff*weight;
    }
  }
  np--; *np=srcsino[srcrays*views - 1];
  return;
}

Referenced by mrp(), and trmrp().

recSinTables()

void recSinTables	(	int	views,
		float *	sinB,
		float *	sinBrot,
		float	rotation )

Pre-compute the sine tables for back-projection.

Parameters

views	Number of views (sinogram rows).
sinB	Array of sine values, length 3*views/2
sinBrot	Array of sine values with rotation, length 3*views/2; enter NULL if not needed.
rotation	Rotation (degrees).

Definition at line 12 of file recutil.c.

  {
  int i, n;
  n=3*views/2;
  for(i=0; i<n; i++) 
    sinB[i]=(float)sin((double)i*(M_PI/(double)views));
  if(sinBrot==NULL) return;
  double rot=M_PI*rotation/180.0;
  for(i=0; i<n; i++) 
    sinBrot[i]=(float)sin((double)i*(M_PI/(double)views) + rot);
}

Referenced by fbp(), imgFBP(), imgMRP(), imgReprojection(), mrp(), reprojection(), and trmrp().

set_os_set()

void set_os_set	(	int	os_sets,
		int *	set_seq )

Make the Ordered Subset process order (bit-reversed sequence).

Parameters

os_sets	Nr of OS sets.
set_seq	Array of length os_sets to be filled here.

Definition at line 113 of file recutil.c.

  {
  if(os_sets==0) {
    return;
  } else if(os_sets==1) {
    set_seq[0]=0;
  } else if(os_sets==2) {
    set_seq[0]=0;
    set_seq[1]=1;
  } else if(os_sets==4) {
    int i=0;
    set_seq[i++]=0; set_seq[i++]=2; set_seq[i++]=1; set_seq[i++]=3;
    return;
  } if(os_sets==8) {
    int i=0;
    set_seq[i++]=0; set_seq[i++]=4; set_seq[i++]=2; set_seq[i++]=6;
    set_seq[i++]=1; set_seq[i++]=5; set_seq[i++]=3; set_seq[i++]=7;
    return;
  } else if(os_sets==16) {
    int i=0;
    set_seq[i++]=0; set_seq[i++]=8;  set_seq[i++]=4; set_seq[i++]=12;
    set_seq[i++]=2; set_seq[i++]=10; set_seq[i++]=6; set_seq[i++]=14;
    set_seq[i++]=1; set_seq[i++]=9;  set_seq[i++]=5; set_seq[i++]=13;
    set_seq[i++]=3; set_seq[i++]=11; set_seq[i++]=7; set_seq[i++]=15;
    return;
  } else {
    set_seq[0]=0; set_seq[1]=1;
    for(int i=1; i<os_sets; i++) set_seq[i]=bit_rev_int(i, os_sets);
  }
}

Referenced by mrp(), and trmrp().