recutil.c

Go to the documentation of this file.

/*****************************************************************************/
#include "libtpcrec.h"
/*****************************************************************************/
 
/*****************************************************************************/
void recSinTables(
  int views, 
  float *sinB, 
  float *sinBrot,
  float rotation
) {
  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);
}
/*****************************************************************************/
 
/*****************************************************************************/
void recInterpolateSinogram(
  float *srcsino,
  float *newsino,
  int srcrays,
  int newrays,
  int views
) {
  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;
}
/*****************************************************************************/
 
/*****************************************************************************/
int bit_rev_int(
  int x,
  int n
) {
  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);
}
/*****************************************************************************/
 
/*****************************************************************************/
void set_os_set(
  int os_sets,
  int *set_seq
) {
  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);
  }
}
/*****************************************************************************/
 
/*****************************************************************************/
int recGetStatistics(
  float *buf,
  int n,
  float *osum,
  float *omin,
  float *omax,
  int skip_zero_mins
) {
  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);
}
/*****************************************************************************/
 
/*****************************************************************************/