image.c

Go to the documentation of this file.

 
/*****************************************************************************/
#include "tpcclibConfig.h"
/*****************************************************************************/
#include "tpcimage.h"
/*****************************************************************************/
 
/*****************************************************************************/
static const char *img_modality[] = {
  "Unknown",  // IMG_MODALITY_UNKNOWN
  "PET",      // IMG_MODALITY_PET
  "SPECT",    // IMG_MODALITY_SPECT
  "CT",       // IMG_MODALITY_CT
  "MRI",      // IMG_MODALITY_MRI
0};
char *imgModalityDescr(
  imgmodality c
) {
  if(c<IMG_MODALITY_UNKNOWN || c>=IMG_MODALITY_LAST) return NULL;
  return (char*)img_modality[c];
}
/*****************************************************************************/
 
/*****************************************************************************/
static const char *img_content[] = {
  "Unknown",     // IMG_CONTENT_UNKNOWN
  "image",       // IMG_CONTENT_IMAGE
  "raw",         // IMG_CONTENT_RAW
  "attenuation", // IMG_CONTENT_ATTN
  "polarmap",    // IMG_CONTENT_POLARMAP
0};
char *imgContentDescr(
  imgcontent c
) {
  if(c<IMG_CONTENT_UNKNOWN || c>=IMG_CONTENT_LAST) return NULL;
  return (char*)img_content[c];
}
/*****************************************************************************/
 
/*****************************************************************************/
void imgInit(
  IMG *img
) {
  if(img==NULL) return;
  img->dimt=img->dimx=img->dimy=img->dimz=0;
  iftInit(&img->ih);
  iftInit(&img->oh);
  img->content=IMG_CONTENT_UNKNOWN;
  img->modality=IMG_MODALITY_UNKNOWN;
  img->format=img->oformat=IMG_FORMAT_UNKNOWN;
 
  for(int i=0; i<20; i++) img->scanStart[i]=(char)0;
  for(int i=0; i<MAX_STUDYNR_LEN+1; i++) img->studyNr[i]=(char)0;
 
 
  img->sizex=img->sizey=img->sizez=0.0;
  img->gapx=img->gapy=img->gapz=0.0;
  img->xform[0]=img->xform[1]=0;
  for(int i=0; i<6; i++) img->quatern[i]=0.0;
  for(int i=0; i<12; i++) img->srow[i]=0.0;
  for(int i=0; i<6; i++) img->iop[i]=0.0;
  for(int i=0; i<3; i++) img->ipp[i]=0.0;
  for(int i=0; i<12; i++) img->mt[i]=0.0;
 
  img->isot=ISOTOPE_UNKNOWN;
  img->decayCorrection=DECAY_UNKNOWN;
 
  img->weighting=WEIGHTING_UNKNOWN;
  img->weight=img->prompts=img->randoms=NULL;
 
  img->cunit=img->tunit=UNIT_UNKNOWN;
 
  img->x1=img->x2=img->x=NULL;
  img->m=(float****)0; img->p=0;
  img->_z=(float****)0; img->_y=(float***)0; img->_x=(float**)0; img->_t=(float*)0;
}
/*****************************************************************************/
 
/*****************************************************************************/
void imgFree(
  IMG *img
) {
  if(img==NULL) return;
  iftFree(&img->ih);
  iftFree(&img->oh);
  free(img->weight); free(img->prompts); free(img->randoms);
  free(img->x1); free(img->x2); free(img->x);
  free(img->_t); free(img->_z); free(img->_x); free(img->_y);
  imgInit(img);
}
/*****************************************************************************/
 
/*****************************************************************************/
int imgAllocate(
  IMG *img,
  const unsigned int dimz,
  const unsigned int dimy,
  const unsigned int dimx,
  const unsigned int dimt,
  TPCSTATUS *status
) {
  int verbose=0; if(status!=NULL) verbose=status->verbose;
  if(verbose>0) {
    printf("%s(img, %u, %u, %u, %u)\n", __func__, dimz, dimy, dimx, dimt); fflush(stdout);}
 
  if(img==NULL || dimz<1 || dimy<1 || dimx<1 || dimt<1) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_FAIL);
    return(TPCERROR_FAIL);
  }
 
  /* Remove any previous contents */
  imgFree(img);
 
  /* Memory for frame times */
  int ret=0; {img->x1=calloc(dimt, sizeof(float)); if(img->x1==NULL) ret++;}
  if(!ret) {img->x2=calloc(dimt, sizeof(float)); if(img->x2==NULL) ret++;}
  if(!ret) {img->x=calloc(dimt, sizeof(float)); if(img->x==NULL) ret++;}
  if(ret) {
    imgFree(img);
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_OUT_OF_MEMORY);
    return(TPCERROR_OUT_OF_MEMORY);
  }
 
  /* Memory for weights */
  ret=0;   {img->weight=calloc(dimt, sizeof(float)); if(img->weight==NULL) ret++;}
  if(!ret) {img->prompts=calloc(dimt, sizeof(float)); if(img->prompts==NULL) ret++;}
  if(!ret) {img->randoms=calloc(dimt, sizeof(float)); if(img->randoms==NULL) ret++;}
  if(ret) {
    imgFree(img);
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_OUT_OF_MEMORY);
    return(TPCERROR_OUT_OF_MEMORY);
  }
 
  /* Memory for pixel matrix */
  ret=0;   {img->_t=(float*)malloc((size_t)dimz*dimy*dimx*dimt*sizeof(float)); if(img->_t==NULL) ret++;}
  if(!ret) {img->_x=(float**)malloc((size_t)dimz*dimy*dimx*sizeof(float*)); if(img->_x==NULL) ret++;}
  if(!ret) {img->_y=(float***)malloc((size_t)dimz*dimy*sizeof(float**)); if(img->_y==NULL) ret++;}
  if(!ret) {img->_z=(float****)malloc((size_t)dimz*sizeof(float***)); if(img->_z==NULL) ret++;}
  if(ret) {
    imgFree(img);
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_OUT_OF_MEMORY);
    return(TPCERROR_OUT_OF_MEMORY);
  }
  /* Set matrix data pointers */
  {
    float ***yptr, **xptr, *tptr;
    yptr=img->_y; xptr=img->_x; tptr=img->_t;
    for(unsigned int zi=0; zi<dimz; zi++) {
      img->_z[zi]=yptr;
      for(unsigned int yi=0; yi<dimy; yi++) {
        *yptr++=xptr;
        for(unsigned int xi=0; xi<dimx; xi++) {
          *xptr++=tptr; tptr+=dimt;
        }
      }
    }
  }
  img->m=img->_z;
  img->p=img->_t;
  /* Set all pixel values to zero */
  for(unsigned int zi=0; zi<dimz; zi++)
    for(unsigned int yi=0; yi<dimy; yi++)
      for(unsigned int xi=0; xi<dimx; xi++)
        for(unsigned int ti=0; ti<dimt; ti++)
          img->m[zi][yi][xi][ti]=(float)0.0;
  /* Set matrix dimensions */
  img->dimz=dimz; img->dimy=dimy; img->dimx=dimx; img->dimt=dimt;
 
  statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_OK);
  return(TPCERROR_OK);
}
/*****************************************************************************/
 
/*****************************************************************************/
int imgHasData(
  IMG *img
) {
  if(img==NULL) return(0);
  if(img->dimz<1 || img->dimy<1 || img->dimx<1 || img->dimz<1) return(0);
  if(img->m==NULL || img->p==NULL) return(0);
  if(img->x1==NULL || img->x2==NULL) return(0);
  return(1);
}
/******************************************************************************/
 
/******************************************************************************/
int imgHasTimes(
  IMG *img
) {
  if(imgHasData(img)==0) return(0);
  for(int i=0; i<img->dimt; i++) if(img->x2[i]>0.00000001) return(1);
  return(0);
}
/******************************************************************************/
 
/******************************************************************************/
int imgHasCounts(
  IMG *img
) {
  if(img==NULL || img->dimt<1) return(0);
  int p=0, r=0;
  /* If image has only one frame, then any value > 0 is fine */
  if(img->dimt==1) {
    if(img->prompts[0]>0.000000001) p=1;
    if(img->randoms[0]>0.000000001) r=2;
    return(p+r);
  }
  /* When more frames, we check that frames have different count level */
  for(unsigned short int i=1; i<img->dimt; i++) {
    if(fabsf(img->prompts[i]-img->prompts[i-1])>0.001) p=1;
    if(fabsf(img->randoms[i]-img->randoms[i-1])>0.001) r=2;
    if((p+r)>2) break;
  }
  return(p+r);
}
/******************************************************************************/
 
/******************************************************************************/
int imgHasWeights(
  IMG *img
) {
  if(img==NULL || img->dimt<1) return(0);
  if(img->weighting==WEIGHTING_UNKNOWN || img->weighting==WEIGHTING_OFF) return(0);
  if(img->weighting<WEIGHTING_LAST) return(1);
  return(0);
}
/*****************************************************************************/
 
/*****************************************************************************/
void imgContents(
  IMG *img,
  FILE *fp
) {
  if(fp==NULL) return;
  if(img==NULL) {fprintf(fp, "IMG pointer is NULL\n"); fflush(fp); return;}
 
  fprintf(fp, "modality: %s\n", imgModalityDescr(img->modality));
  fprintf(fp, "content: %s\n", imgContentDescr(img->content));
  fprintf(fp, "format: %s\n", imgFormatDescr(img->format));
  if(img->oformat!=IMG_FORMAT_UNKNOWN)
    fprintf(fp, "output_format: %s\n", imgFormatDescr(img->oformat));
 
  fprintf(fp, "hasData := %d\n", imgHasData(img));
  fprintf(fp, "hasTimes := %d\n", imgHasTimes(img));
  fprintf(fp, "hasCounts := %d\n", imgHasCounts(img));
 
  if(img->scanStart[0]) fprintf(fp, "scanStart: %s\n", img->scanStart);
 
  fprintf(fp, "dimensions: %u x %u x %u x %u\n", img->dimx, img->dimy, img->dimz, img->dimt);
  if(img->sizex!=0.0 || img->sizey!=0.0 || img->sizez!=0.0)
    fprintf(fp, "voxel_size: %g x %g x %g\n", img->sizex, img->sizey, img->sizez);
  if(img->gapx!=0.0 || img->gapy!=0.0 || img->gapz!=0.0)
    fprintf(fp, "voxel_gaps: %g x %g x %g\n", img->gapx, img->gapy, img->gapz);
 
  fprintf(fp, "qform := %d\n", img->xform[0]);
  fprintf(fp, "sform := %d\n", img->xform[1]);
  if(floatNonzeroes(img->quatern, 6)>0) {
    fprintf(fp, "quatern_b := %g\n", img->quatern[0]);
    fprintf(fp, "quatern_c := %g\n", img->quatern[1]);
    fprintf(fp, "quatern_d := %g\n", img->quatern[2]);
    fprintf(fp, "quatern_x_shift := %g\n", img->quatern[3]);
    fprintf(fp, "quatern_y_shift := %g\n", img->quatern[4]);
    fprintf(fp, "quatern_z_shift := %g\n", img->quatern[5]);
  }
  if(floatNonzeroes(img->srow, 12)>0) {
    for(int i=0; i<4; i++) fprintf(fp, "srow_x[%d] := %g\n", i, img->srow[i]);
    for(int i=0; i<4; i++) fprintf(fp, "srow_y[%d] := %g\n", i, img->srow[4+i]);
    for(int i=0; i<4; i++) fprintf(fp, "srow_z[%d] := %g\n", i, img->srow[8+i]);
  }
  if(floatNonzeroes(img->iop, 6)>0) {
    fprintf(fp, "image_orientation_patient := %g", img->iop[0]);
    for(int i=1; i<6; i++) fprintf(fp, ", %g", img->iop[i]);
    fprintf(fp, "\n");
  }
  if(floatNonzeroes(img->ipp, 3)>0) {
    fprintf(fp, "image_position_patient := %g", img->ipp[0]);
    for(int i=1; i<3; i++) fprintf(fp, ", %g", img->ipp[i]);
    fprintf(fp, "\n");
  }
  if(floatNonzeroes(img->mt, 3)>0) {
    fprintf(fp, "image_matrix_transformation_parameters := %g", img->mt[0]);
    for(int i=1; i<12; i++) fprintf(fp, ", %g", img->mt[i]);
    fprintf(fp, "\n");
  }
 
  if(img->isot!=ISOTOPE_UNKNOWN) fprintf(fp, "isotope: %s\n", isotopeName(img->isot));
  if(img->decayCorrection==DECAY_CORRECTED) fprintf(fp, "physical_decay: corrected\n");
  else if(img->decayCorrection==DECAY_NOTCORRECTED) fprintf(fp, "physical_decay: not corrected\n");
 
  if(imgHasWeights(img)) fprintf(fp, "weighting: no\n");
  else fprintf(fp, "weighting: yes\n");
 
  if(img->cunit!=UNIT_UNKNOWN) fprintf(fp, "pixel_unit := %s\n", unitName(img->cunit));
  if(img->tunit!=UNIT_UNKNOWN) fprintf(fp, "time_unit := %s\n", unitName(img->tunit));
 
  fprintf(fp, "original header: %d items\n", img->ih.keyNr);
 
  fflush(fp);
  return;
}
/*****************************************************************************/
 
/*****************************************************************************/
unsigned long long imgNaNs(
  IMG *img,
  int fix
) {
  if(img==NULL) return(0);
  unsigned long long i, n=0, pxlNr=img->dimt*img->dimz*img->dimy*img->dimx;
  for(i=0; i<pxlNr; i++) {
    if(!isfinite(img->p[i])) {
      n++;
      if(fix!=0) img->p[i]=0.0;
    }
  }
  return(n);
}
/*****************************************************************************/
 
/*****************************************************************************/
int imgMinMax(
  IMG *img,
  float *minvalue,
  float *maxvalue
) {
  if(img==NULL || !imgHasData(img)) return(TPCERROR_FAIL);
  float f1, f2;
  f1=f2=nanf("");
  if(minvalue!=NULL) *minvalue=f1;
  if(maxvalue!=NULL) *maxvalue=f2;
  unsigned long long i, n=img->dimt*img->dimz*img->dimy*img->dimx;
  // search first valid number
  for(i=0; i<n; i++) if(isfinite(img->p[i])) break;
  if(i==n) return(TPCERROR_MISSING_DATA);
  f1=f2=img->p[i++];
  for(; i<n; i++) {
    if(img->p[i]>f2) f2=img->p[i];
    else if(img->p[i]<f1) f1=img->p[i];
  }
  if(minvalue!=NULL) *minvalue=f1;
  if(maxvalue!=NULL) *maxvalue=f2;
  return(TPCERROR_OK);
}
/*****************************************************************************/
 
/*****************************************************************************/
int imgXRange(
  IMG *img,
  double *xmin,
  double *xmax
) {
  if(xmin!=NULL) *xmin=nan("");
  if(xmax!=NULL) *xmax=nan("");
  if(img==NULL || !imgHasData(img)) return(TPCERROR_FAIL);
  float f1=nanf(""), f2=nanf("");
  for(int i=0; i<img->dimt; i++) {
    if(isfinite(img->x1[i])) {
      if(isnan(f1) || img->x1[i]<f1) f1=img->x1[i];
    }
    if(isfinite(img->x2[i])) {
      if(isnan(f2) || img->x2[i]>f2) f2=img->x2[i];
    }
  }
  if(xmin!=NULL) *xmin=f1;
  if(xmax!=NULL) *xmax=f2;
  if(isnan(f1) || isnan(f2)) return(TPCERROR_FAIL);
  return(TPCERROR_OK);
}
/*****************************************************************************/
 
/*****************************************************************************/