imgarithm.c

Go to the documentation of this file.

 
/*****************************************************************************/
#include "libtpcimgp.h"
/*****************************************************************************/
 
/*****************************************************************************/
int imgArithm(
  IMG *img1,
  IMG *img2,
  char operation,
  float ulimit,
  int verbose
) {
  if(verbose>0) printf("%s(img1, img2, '%c', %g, %d)\n", __func__, operation, ulimit, verbose);
 
  /* Check the arguments */
  if(img1->status!=IMG_STATUS_OCCUPIED || img2->status!=IMG_STATUS_OCCUPIED) {
    if(verbose>0) fprintf(stderr, "Invalid image status.\n");
    return(1);
  }
  int ret=0;
  if(img1->dimx!=img2->dimx || img1->dimy!=img2->dimy) ret=1;
  /* Check the plane numbers */
  if(img1->dimz!=img2->dimz) ret=2;
  if(ret==0) for(int pi=0; pi<img1->dimz; pi++)
    if(img1->planeNumber[pi]!=img2->planeNumber[pi]) ret=2;
  /* Check the frame numbers */
  if(img1->dimt!=img2->dimt) ret=3;
  if(ret>0) {
    if(verbose>0) fprintf(stderr, "Image dimensions do not match.\n");
    return(ret);
  }
 
  /* Operate */
  switch(operation) {
    case '+':
      for(int pi=0; pi<img1->dimz; pi++) for(int yi=0; yi<img1->dimy; yi++)
        for(int xi=0; xi<img1->dimx; xi++) for(int fi=0; fi<img1->dimt; fi++)
          img1->m[pi][yi][xi][fi]+=img2->m[pi][yi][xi][fi];
      break;
    case '-':
      for(int pi=0; pi<img1->dimz; pi++) for(int yi=0; yi<img1->dimy; yi++)
        for(int xi=0; xi<img1->dimx; xi++) for(int fi=0; fi<img1->dimt; fi++)
          img1->m[pi][yi][xi][fi]-=img2->m[pi][yi][xi][fi];
      break;
    case '/':
    case ':':
      for(int pi=0; pi<img1->dimz; pi++) for(int yi=0; yi<img1->dimy; yi++)
        for(int xi=0; xi<img1->dimx; xi++) for(int fi=0; fi<img1->dimt; fi++) {
          if(fabs(img2->m[pi][yi][xi][fi])>1.0E-5)
            img1->m[pi][yi][xi][fi]/=img2->m[pi][yi][xi][fi];
          else img1->m[pi][yi][xi][fi]=0.0;
        }
      break;
    case '*':
    case 'x':
    case '.':
      for(int pi=0; pi<img1->dimz; pi++) for(int yi=0; yi<img1->dimy; yi++)
        for(int xi=0; xi<img1->dimx; xi++) for(int fi=0; fi<img1->dimt; fi++)
          img1->m[pi][yi][xi][fi]*=img2->m[pi][yi][xi][fi];
      break;
    default:
      if(verbose>0) fprintf(stderr, "Invalid operation.\n");
      return(10);
  }
 
  /* Check for limits */
  if(ulimit>0.0) {
    for(int pi=0; pi<img1->dimz; pi++) for(int yi=0; yi<img1->dimy; yi++)
      for(int xi=0; xi<img1->dimx; xi++) for(int fi=0; fi<img1->dimt; fi++)
        if(img1->m[pi][yi][xi][fi]>ulimit) img1->m[pi][yi][xi][fi]=ulimit;
  }
 
  return(0);
}
/*****************************************************************************/
 
/*****************************************************************************/
int imgArithmConst(
  IMG *img,
  float operand,
  char operation,
  float ulimit,
  int verbose
) {
  if(verbose>0) printf("%s(img, %g, '%c', %g, %d)\n", __func__, operand, operation, ulimit, verbose);
 
  /* Check the arguments */
  if(img->status!=IMG_STATUS_OCCUPIED) {
    if(verbose>0) fprintf(stderr, "Invalid image status.\n");
    return(1);
  }
  switch(operation) {
    case '+':
      for(int pi=0; pi<img->dimz; pi++) for(int yi=0; yi<img->dimy; yi++)
        for(int xi=0; xi<img->dimx; xi++) for(int fi=0; fi<img->dimt; fi++)
          img->m[pi][yi][xi][fi]+=operand;
      break;
    case '-':
      for(int pi=0; pi<img->dimz; pi++) for(int yi=0; yi<img->dimy; yi++)
        for(int xi=0; xi<img->dimx; xi++) for(int fi=0; fi<img->dimt; fi++)
          img->m[pi][yi][xi][fi]-=operand;
      break;
    case '/':
    case ':':
      if(fabs(operand)<1.0e-100) return(2);
      for(int pi=0; pi<img->dimz; pi++) for(int yi=0; yi<img->dimy; yi++)
        for(int xi=0; xi<img->dimx; xi++) for(int fi=0; fi<img->dimt; fi++)
          img->m[pi][yi][xi][fi]/=operand;
      break;
    case '*':
    case 'x':
    case '.':
      for(int pi=0; pi<img->dimz; pi++) for(int yi=0; yi<img->dimy; yi++)
        for(int xi=0; xi<img->dimx; xi++) for(int fi=0; fi<img->dimt; fi++)
          img->m[pi][yi][xi][fi]*=operand;
      break;
    default:
      if(verbose>0) fprintf(stderr, "Invalid operation.\n");
      return(10);
  }
 
  /* Check for limits */
  if(ulimit>0.0) {
    for(int pi=0; pi<img->dimz; pi++) for(int yi=0; yi<img->dimy; yi++)
      for(int xi=0; xi<img->dimx; xi++) for(int fi=0; fi<img->dimt; fi++)
        if(img->m[pi][yi][xi][fi]>ulimit) img->m[pi][yi][xi][fi]=ulimit;
  }
 
  return(0);
}
/*****************************************************************************/
 
/*****************************************************************************/
int imgArithmFrame(
  IMG *img1,
  IMG *img2,
  char operation,
  float ulimit,
  int verbose
) {
  if(verbose>0) printf("%s(img1, img2, '%c', %g, %d)\n", __func__, operation, ulimit, verbose);
 
  /* Check the arguments */
  if(img1->status!=IMG_STATUS_OCCUPIED || img2->status!=IMG_STATUS_OCCUPIED) {
    if(verbose>0) fprintf(stderr, "Invalid image status.\n");
    return(1);
  }
  int ret=0;
  if(img1->dimx!=img2->dimx || img1->dimy!=img2->dimy) ret=1;
  /* Check the plane numbers */
  if(img1->dimz!=img2->dimz) ret=2;
  if(ret==0) for(int pi=0; pi<img1->dimz; pi++)
    if(img1->planeNumber[pi]!=img2->planeNumber[pi]) ret=2;
  if(ret>0) {
    if(verbose>0) fprintf(stderr, "Image dimensions do not match.\n");
    return(ret);
  }
 
  /* Operate */
  switch(operation) {
    case '+':
      for(int pi=0; pi<img1->dimz; pi++) for(int yi=0; yi<img1->dimy; yi++)
        for(int xi=0; xi<img1->dimx; xi++) for(int fi=0; fi<img1->dimt; fi++)
          img1->m[pi][yi][xi][fi]+=img2->m[pi][yi][xi][0];
      break;
    case '-':
      for(int pi=0; pi<img1->dimz; pi++) for(int yi=0; yi<img1->dimy; yi++)
        for(int xi=0; xi<img1->dimx; xi++) for(int fi=0; fi<img1->dimt; fi++)
          img1->m[pi][yi][xi][fi]-=img2->m[pi][yi][xi][0];
      break;
    case '/':
    case ':':
      for(int pi=0; pi<img1->dimz; pi++) for(int yi=0; yi<img1->dimy; yi++)
        for(int xi=0; xi<img1->dimx; xi++) {
          if(fabs(img2->m[pi][yi][xi][0])>1.0E-8)
            for(int fi=0; fi<img1->dimt; fi++)
              img1->m[pi][yi][xi][fi]/=img2->m[pi][yi][xi][0];
          else for(int fi=0; fi<img1->dimt; fi++) img1->m[pi][yi][xi][fi]=0.0;
        }
      break;
    case '*':
    case 'x':
    case '.':
      for(int pi=0; pi<img1->dimz; pi++) for(int yi=0; yi<img1->dimy; yi++)
        for(int xi=0; xi<img1->dimx; xi++) for(int fi=0; fi<img1->dimt; fi++)
          img1->m[pi][yi][xi][fi]*=img2->m[pi][yi][xi][0];
      break;
    default:
      if(verbose>0) fprintf(stderr, "Invalid operation.\n");
      return(10);
  }
 
  /* Check for limits */
  if(ulimit>0.0) {
    for(int pi=0; pi<img1->dimz; pi++) for(int yi=0; yi<img1->dimy; yi++)
      for(int xi=0; xi<img1->dimx; xi++) for(int fi=0; fi<img1->dimt; fi++)
        if(img1->m[pi][yi][xi][fi]>ulimit) img1->m[pi][yi][xi][fi]=ulimit;
  }
 
  return(0);
}
/*****************************************************************************/
 
/*****************************************************************************/
int imgLn(
  IMG *img
) {
  if(img->status<IMG_STATUS_OCCUPIED) return(1);
  if(img->dimt<1 || img->dimz<1 || img->dimy<1 || img->dimx<1) return(2);
  for(int pi=0; pi<img->dimz; pi++)
    for(int yi=0; yi<img->dimy; yi++)
      for(int xi=0; xi<img->dimx; xi++)
        for(int fi=0; fi<img->dimt; fi++) {
          if(img->m[pi][yi][xi][fi]<=0.0) img->m[pi][yi][xi][fi]=0.0;
          else img->m[pi][yi][xi][fi]=log(img->m[pi][yi][xi][fi]);
        }
  return(0);
}
/*****************************************************************************/
 
/*****************************************************************************/
int imgLog10(
  IMG *img
) {
  if(img->status<IMG_STATUS_OCCUPIED) return(1);
  if(img->dimt<1 || img->dimz<1 || img->dimy<1 || img->dimx<1) return(2);
  for(int pi=0; pi<img->dimz; pi++)
    for(int yi=0; yi<img->dimy; yi++)
      for(int xi=0; xi<img->dimx; xi++)
        for(int fi=0; fi<img->dimt; fi++) {
          if(img->m[pi][yi][xi][fi]<=0.0) img->m[pi][yi][xi][fi]=0.0;
          else img->m[pi][yi][xi][fi]=log10(img->m[pi][yi][xi][fi]);
        }
  return(0);
}
/*****************************************************************************/
 
/*****************************************************************************/
int imgAbs(
  IMG *img
) {
  if(img->status<IMG_STATUS_OCCUPIED) return(1);
  if(img->dimt<1 || img->dimz<1 || img->dimy<1 || img->dimx<1) return(2);
  for(int pi=0; pi<img->dimz; pi++)
    for(int yi=0; yi<img->dimy; yi++)
      for(int xi=0; xi<img->dimx; xi++)
        for(int fi=0; fi<img->dimt; fi++) 
          img->m[pi][yi][xi][fi]=fabsf(img->m[pi][yi][xi][fi]);
  return(0);
}
/*****************************************************************************/
 
/*****************************************************************************/
int imgInv(
  IMG *img
) {
  if(img->status<IMG_STATUS_OCCUPIED) return(1);
  if(img->dimt<1 || img->dimz<1 || img->dimy<1 || img->dimx<1) return(2);
  for(int pi=0; pi<img->dimz; pi++)
    for(int yi=0; yi<img->dimy; yi++)
      for(int xi=0; xi<img->dimx; xi++)
        for(int fi=0; fi<img->dimt; fi++) {
          if(img->m[pi][yi][xi][fi]<=0.0) img->m[pi][yi][xi][fi]=0.0;
          else {
            img->m[pi][yi][xi][fi]=1.0/img->m[pi][yi][xi][fi];
            if(!isfinite(img->m[pi][yi][xi][fi])) img->m[pi][yi][xi][fi]=0.0;
          }
        }
  return(0);
}
/*****************************************************************************/
 
/*****************************************************************************/
int imgFrameIntegral(
  IMG *img,
  int first,
  int last,
  IMG *iimg,
  int verbose
) {
  if(verbose>0) printf("%s(img, %d, %d, iimg, %d)\n", __func__, first, last, verbose);
 
  int zi, yi, xi, fi, ret, times_exist;
  float fstart, fend, dur, x, y, k;
 
  /* Check the arguments */
  if(img==NULL || iimg==NULL || first<0 || first>last) return(1);
  if(img->status!=IMG_STATUS_OCCUPIED) return(2);
  times_exist=imgExistentTimes(img);
  fstart=img->start[first]; fend=img->end[last];
  if(verbose>1) printf("  time_range := %g - %g\n", fstart, fend);
 
  /* Allocate memory for the integral */
  imgEmpty(iimg);
  ret=imgAllocateWithHeader(iimg, img->dimz, img->dimy, img->dimx, 1, img);
  if(ret) {
    if(verbose>0) fprintf(stderr, "Error %d in allocating memory for sum image.\n", ret);
    imgEmpty(iimg); return(3);
  }
 
  /* Set header information */
  iimg->start[0]=fstart; iimg->end[0]=fend;
  iimg->mid[0]=0.5*(iimg->start[0]+iimg->end[0]);
  if(img->type==IMG_TYPE_RAW) iimg->decayCorrection=IMG_DC_NONCORRECTED;
  else iimg->decayCorrection=IMG_DC_CORRECTED;
  iimg->decayCorrFactor[0]=0.0;
  /* Change also the unit, if possible */
  if(img->type==IMG_TYPE_IMAGE) {
    if(img->unit==CUNIT_KBQ_PER_ML) iimg->unit=CUNIT_SEC_KBQ_PER_ML;
  }
 
  /* Integrate the first frame */
  fi=first;
  if(times_exist) {
    dur=img->end[fi]-img->start[fi]; if(dur<0.0) {imgEmpty(iimg); return(4);}
  } else dur=1.0;
  for(zi=0; zi<img->dimz; zi++)
    for(yi=0; yi<img->dimy; yi++)
      for(xi=0; xi<img->dimx; xi++)
        iimg->m[zi][yi][xi][0]=dur*img->m[zi][yi][xi][fi];
  /* Add integrals of the following frames */
  for(fi=first+1; fi<=last; fi++) {
    if(times_exist) {
      dur=img->end[fi]-img->start[fi]; if(dur<0.0) {imgEmpty(iimg); return(4);}
    } else dur=1.0;
    /* Add the integral of this frame to the previous integral */
    for(zi=0; zi<img->dimz; zi++)
      for(yi=0; yi<img->dimy; yi++)
        for(xi=0; xi<img->dimx; xi++)
          iimg->m[zi][yi][xi][0]+=dur*img->m[zi][yi][xi][fi];
    if(times_exist) {
      /* Check whether frames are contiguous */
      dur=img->start[fi]-img->end[fi-1]; if(dur<=1.0E-10) continue;
      /* When not, calculate the integral between frames */
      x=0.5*(img->start[fi]+img->end[fi-1]);
      for(zi=0; zi<img->dimz; zi++)
        for(yi=0; yi<img->dimy; yi++)
          for(xi=0; xi<img->dimx; xi++) {
            k=(img->m[zi][yi][xi][fi]-img->m[zi][yi][xi][fi-1])
              /(img->mid[fi]-img->mid[fi-1]);
            y=img->m[zi][yi][xi][fi-1]+k*(x-img->mid[fi-1]);
            iimg->m[zi][yi][xi][0]+=dur*y;
          }
    }    
  } /* next frame */
 
  /* Set frame times */
  iimg->start[0]=img->start[first];
  iimg->end[0]=img->end[last];
  iimg->mid[0]=0.5*(iimg->start[0]+iimg->end[0]);
 
  return(0);
}
/*****************************************************************************/
 
/*****************************************************************************/
int imgRawCountsPerTime(
  IMG *img,
  int operation
) {
  int fi, zi, xi, yi, notime=0;
  float f;
 
  /* Check the input data */
  if(operation!=0 && operation!=1) return(1);
  if(img==NULL) return(2);
  if(img->status!=IMG_STATUS_OCCUPIED) return(3);
  if(img->type!=IMG_TYPE_RAW) return(0);
  for(fi=0; fi<img->dimt; fi++) {
    f=img->end[fi]-img->start[fi]; if(f<=1.0E-12) notime++;
  }
  if(notime>0) {
    if(img->dimt>1) return(4);
    else return(0); /* just one frame; might be parametric or transmission sinogram,
                       therefore no error. */
  }
  for(fi=0; fi<img->dimt; fi++) {
    f=img->end[fi]-img->start[fi]; if(operation==1) f=1.0/f;
    for(zi=0; zi<img->dimz; zi++)
      for(yi=0; yi<img->dimy; yi++)
        for(xi=0; xi<img->dimx; xi++)
          img->m[zi][yi][xi][fi]*=f;
  }
  return(0);
}
/*****************************************************************************/
 
/*****************************************************************************/
int imgConvertUnit(
  IMG *img,
  char *unit
) {
  int ret, new_unit=0;
  float conversion_factor=1.0;
 
  if(img==NULL || unit==NULL) return(1);
  if(img->unit==CUNIT_UNKNOWN) return(1);
  /* Try to identify requested unit */
  new_unit=imgUnitId(unit); if(new_unit<0) return(1-new_unit);
  /* Check if unit needs no conversion */
  if(img->unit==new_unit) return(0);
  /* Get conversion factor */
  if(img->unit==CUNIT_KBQ_PER_ML && new_unit==CUNIT_BQ_PER_ML)
    conversion_factor=1000.0;
  else if(img->unit==CUNIT_BQ_PER_ML && new_unit==CUNIT_KBQ_PER_ML)
    conversion_factor=0.001;
  else if(img->unit==CUNIT_KBQ_PER_ML && new_unit==CUNIT_NCI_PER_ML)
    conversion_factor=27.027;
  else if(img->unit==CUNIT_NCI_PER_ML && new_unit==CUNIT_KBQ_PER_ML)
    conversion_factor=0.037;
  else if(img->unit==CUNIT_NCI_PER_ML && new_unit==CUNIT_BQ_PER_ML)
    conversion_factor=37.0;
  else if(img->unit==CUNIT_KBQ_PER_ML && new_unit==CUNIT_MBQ_PER_ML)
    conversion_factor=0.001;
  else if(img->unit==CUNIT_MBQ_PER_ML && new_unit==CUNIT_KBQ_PER_ML)
    conversion_factor=1000.0;
  else if(img->unit==CUNIT_PER_SEC && new_unit==CUNIT_PER_MIN)
    conversion_factor=60.0;
  else if(img->unit==CUNIT_PER_MIN && new_unit==CUNIT_PER_SEC)
    conversion_factor=1.0/60.0;
  else if(img->unit==CUNIT_ML_PER_ML && new_unit==CUNIT_ML_PER_DL)
    conversion_factor=0.01;
  else if(img->unit==CUNIT_ML_PER_DL && new_unit==CUNIT_ML_PER_ML)
    conversion_factor=100.0;
  else
    return(10);
 
  /* Convert pixel values */
  ret=imgArithmConst(img, conversion_factor, '*', FLT_MAX, 0);
  if(ret) return(10+ret);
 
  /* Set the unit id */
  img->unit=new_unit;
 
  return(0);
}
/*****************************************************************************/
 
/*****************************************************************************/
int imgAUMC(
  IMG *img,
  IMG *oimg,
  int verbose
) {
  if(verbose>0) printf("%s(*img, *oimg)\n", __func__);
  
  if(img==NULL || img->status!=IMG_STATUS_OCCUPIED) return(1);
  if(img->dimt<2) {
    if(verbose>1) fprintf(stderr, "Error: dynamic image required.\n");
    return(2);
  }
  if(!imgExistentTimes(img)) {
    if(verbose>1) fprintf(stderr, "Error: unknown frame times.\n");
    return(2);
  }
  if(oimg==NULL) return(3);
  if(oimg->status==IMG_STATUS_OCCUPIED) imgEmpty(oimg);
 
  if(imgFramesCheck(img, verbose-1)>0) {
    if(verbose>1) fprintf(stderr, "Error: gaps or overlap between time frames.\n");
    return(4);
  }
 
  /* Allocate memory for one frame */
  int ret;
  if(verbose>1) printf("allocating memory for %dx%dx%d pixels\n", img->dimz, img->dimy, img->dimx);
  ret=imgAllocate(oimg, img->dimz, img->dimy, img->dimx, 1);
  if(ret) return(ret);
  /* set image header information */
  imgCopyhdr(img, oimg);
  oimg->start[0]=img->start[0]; oimg->end[0]=img->end[img->dimt-1];
  oimg->mid[0]=(oimg->start[0]+oimg->end[0])/2.0;
  oimg->unit=CUNIT_UNKNOWN;
 
  /* Calculate sum of (frame time * pixel value) */
  for(int zi=0; zi<img->dimz; zi++) {
    for(int yi=0; yi<img->dimy; yi++) {
      for(int xi=0; xi<img->dimx; xi++) {
        oimg->m[zi][yi][xi][0]=0.0;
        for(int ti=0; ti<img->dimt; ti++) {
          if(isnan(img->m[zi][yi][xi][ti])) continue;
          float fdur=img->end[ti]-img->start[ti]; if(!(fdur>0.0)) fdur=0.0;
          oimg->m[zi][yi][xi][0]+=img->m[zi][yi][xi][ti]*img->mid[ti]*fdur;
        }
      }
    }
  }
 
  return(0);
}
/*****************************************************************************/
 
/*****************************************************************************/
int imgMRT(
  IMG *img,
  IMG *oimg,
  int verbose
) {
  if(verbose>0) printf("%s(*img, *oimg)\n", __func__);
  
  if(img==NULL || img->status!=IMG_STATUS_OCCUPIED) return(1);
  if(img->dimt<2) {
    if(verbose>1) fprintf(stderr, "Error: dynamic image required.\n");
    return(2);
  }
  if(!imgExistentTimes(img)) {
    if(verbose>1) fprintf(stderr, "Error: unknown frame times.\n");
    return(2);
  }
  if(oimg==NULL) return(3);
  if(oimg->status==IMG_STATUS_OCCUPIED) imgEmpty(oimg);
 
  if(imgFramesCheck(img, verbose-1)>0) {
    if(verbose>1) fprintf(stderr, "Error: gaps or overlap between time frames.\n");
    return(4);
  }
 
  /* Allocate memory for one frame */
  int ret;
  if(verbose>1) printf("allocating memory for %dx%dx%d pixels\n", img->dimz, img->dimy, img->dimx);
  ret=imgAllocate(oimg, img->dimz, img->dimy, img->dimx, 1);
  if(ret) return(ret);
  /* set image header information */
  imgCopyhdr(img, oimg);
  oimg->start[0]=img->start[0]; oimg->end[0]=img->end[img->dimt-1];
  oimg->mid[0]=(oimg->start[0]+oimg->end[0])/2.0;
  oimg->unit=CUNIT_UNKNOWN;
 
  /* Calculate AUMC, AUC, and MRT */
  for(int zi=0; zi<img->dimz; zi++) {
    for(int yi=0; yi<img->dimy; yi++) {
      for(int xi=0; xi<img->dimx; xi++) {
        float aumc=0.0, auc=0.0;
        for(int ti=0; ti<img->dimt; ti++) {
          if(isnan(img->m[zi][yi][xi][ti])) continue;
          float fdur=img->end[ti]-img->start[ti]; if(!(fdur>0.0)) fdur=0.0;
          aumc+=img->m[zi][yi][xi][ti]*img->mid[ti]*fdur;
          auc+=img->m[zi][yi][xi][ti]*fdur;
        }
        float mrt=aumc/auc;
        if(isfinite(mrt)) oimg->m[zi][yi][xi][0]=mrt; else oimg->m[zi][yi][xi][0]=0.0;
      }
    }
  }
 
  return(0);
}
/*****************************************************************************/
 
/*****************************************************************************/