imgminmax.c

Go to the documentation of this file.

 
/******************************************************************************/
#include "libtpcimgio.h"
/******************************************************************************/
 
/******************************************************************************/
int imgMax(
  IMG *img,
  float *maxvalue
) {
  if(img->status<IMG_STATUS_OCCUPIED) return(1);
  if(maxvalue==NULL) return(2); else *maxvalue=0.0;
  if(img->dimt<1 || img->dimz<1 || img->dimy<1 || img->dimx<1) return(3);
  float f=nanf("");
  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(!isfinite(f) || img->m[pi][yi][xi][fi]>f) f=img->m[pi][yi][xi][fi];
        }
  *maxvalue=f;
  return(0);
}
/******************************************************************************/
 
/******************************************************************************/
int imgAbsMax(
  IMG *img,
  float *maxvalue
) {
  if(img->status<IMG_STATUS_OCCUPIED) return(1);
  if(maxvalue==NULL) return(2); else *maxvalue=0.0;
  if(img->dimt<1 || img->dimz<1 || img->dimy<1 || img->dimx<1) return(3);
  float f=nanf("");
  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(!isfinite(f) || fabs(img->m[pi][yi][xi][fi])>fabs(f)) f=img->m[pi][yi][xi][fi];
        }
  *maxvalue=f;
  return(0);
}
/******************************************************************************/
 
/******************************************************************************/
int imgRangeMinMax(
  IMG *img,
  IMG_RANGE *r,
  IMG_PIXEL *maxp,
  float *maxv,
  IMG_PIXEL *minp,
  float *minv
) {
  int zi, yi, xi, fi;
  float lmax, lmin;
 
  if(img->status<IMG_STATUS_OCCUPIED) return(1);
  if(img->dimt<1 || img->dimz<1 || img->dimy<1 || img->dimx<1) return(1);
 
  if(r!=NULL) {
    if(r->z1<1 || r->y1<1 || r->x1<1 || r->f1<1) return(2);
    if(r->z2<r->z1 || r->y2<r->y1 || r->x2<r->x1 || r->f2<r->f1) return(3);
    if(r->z2>img->dimz || r->y2>img->dimy || r->x2>img->dimx || r->f2>img->dimt) return(4);
 
    zi=r->z1-1; yi=r->y1-1; xi=r->x1-1; fi=r->f1-1;
    lmax=lmin=nanf("");
    if(maxp!=NULL) {maxp->z=zi+1; maxp->y=yi+1; maxp->x=xi+1; maxp->f=fi+1;}
    if(minp!=NULL) {minp->z=zi+1; minp->y=yi+1; minp->x=xi+1; minp->f=fi+1;}
    for(zi=r->z1-1; zi<r->z2; zi++) {
      for(yi=r->y1-1; yi<r->y2; yi++) {
        for(xi=r->x1-1; xi<r->x2; xi++) {
          for(fi=r->f1-1; fi<r->f2; fi++) {
            if(!isfinite(lmax) || img->m[zi][yi][xi][fi]>lmax) {
              lmax=img->m[zi][yi][xi][fi];
              if(maxp!=NULL && isfinite(img->m[zi][yi][xi][fi])) {
                maxp->z=zi+1; maxp->y=yi+1; maxp->x=xi+1; maxp->f=fi+1;}
            }
            if(!isfinite(lmin) || img->m[zi][yi][xi][fi]<lmin) {
              lmin=img->m[zi][yi][xi][fi];
              if(minp!=NULL && isfinite(img->m[zi][yi][xi][fi])) {
                minp->z=zi+1; minp->y=yi+1; minp->x=xi+1; minp->f=fi+1;}
            }
          }
        }
      }
    }
  } else {
    zi=yi=xi=fi=0; lmax=lmin=nanf("");
    if(maxp!=NULL) {maxp->z=zi+1; maxp->y=yi+1; maxp->x=xi+1; maxp->f=fi+1;}
    if(minp!=NULL) {minp->z=zi+1; minp->y=yi+1; minp->x=xi+1; minp->f=fi+1;}
    for(zi=0; zi<img->dimz; zi++) {
      for(yi=0; yi<img->dimy; yi++) {
        for(xi=0; xi<img->dimx; xi++) {
          for(fi=0; fi<img->dimt; fi++) {
            if(!isfinite(lmax) || img->m[zi][yi][xi][fi]>lmax) {
              lmax=img->m[zi][yi][xi][fi];
              if(maxp!=NULL && isfinite(img->m[zi][yi][xi][fi])) {
                maxp->z=zi+1; maxp->y=yi+1; maxp->x=xi+1; maxp->f=fi+1;}
            }
            if(!isfinite(lmin) || img->m[zi][yi][xi][fi]<lmin) {
              lmin=img->m[zi][yi][xi][fi];
              if(minp!=NULL && isfinite(img->m[zi][yi][xi][fi])) {
                minp->z=zi+1; minp->y=yi+1; minp->x=xi+1; minp->f=fi+1;}
            }
          }
        }
      }
    }
  }
  if(maxv!=NULL) *maxv=lmax;
  if(minv!=NULL) *minv=lmin;
  if(!isfinite(lmax) && (maxp!=NULL || maxv!=NULL)) return(5); 
  if(!isfinite(lmin) && (minp!=NULL || minv!=NULL)) return(5); 
  return(0);
}
/******************************************************************************/
 
/******************************************************************************/
int imgMinMax(
  IMG *img,
  float *minvalue,
  float *maxvalue
) {
  return (imgRangeMinMax(img, NULL, NULL, maxvalue, NULL, minvalue));
}
/******************************************************************************/
 
/******************************************************************************/
int imgFrameMinMax(
  IMG *img, 
  int frame, 
  float *minvalue, 
  float *maxvalue
) {
  if(img->status<IMG_STATUS_OCCUPIED) return(1);
  if(minvalue==NULL || maxvalue==NULL) return(2);
  *minvalue=*maxvalue=0.0; 
  int fi=frame-1;
  if(img->dimt<frame || img->dimz<1 || img->dimy<1 || img->dimx<1) return(3);
  if(frame<1) return(4);
 
  float mi, ma;
  mi=ma=nanf("");
  for(int pi=0; pi<img->dimz; pi++)
    for(int yi=0; yi<img->dimy; yi++)
      for(int xi=0; xi<img->dimx; xi++) {
        if(!isfinite(ma) || img->m[pi][yi][xi][fi]>ma) ma=img->m[pi][yi][xi][fi];
        if(!isfinite(mi) || img->m[pi][yi][xi][fi]<mi) mi=img->m[pi][yi][xi][fi];
      }
  *minvalue=mi; *maxvalue=ma;
  if(!isfinite(mi) || !isfinite(ma)) return(5);
  return(0);
}
/******************************************************************************/
 
/******************************************************************************/
int imgReadMinMax(
  const char *fname, 
  float *fmin, 
  float *fmax
) {
  int fi=0, ret;
  IMG img;
  float frmin, frmax;
 
  if(IMG_TEST) printf("imgReadMinMax(%s, *fmin, *fmax)\n", fname);
  imgInit(&img);
  while((ret=imgReadFrame(fname, fi+1, &img, 0)) == 0) {
    if(imgMinMax(&img, &frmin, &frmax)!=0) {imgEmpty(&img); return STATUS_FAULT;}
    if(fi==0) {
      if(fmin!=NULL) *fmin=frmin;
      if(fmin!=NULL) *fmax=frmax;
    } else {
      if(fmin!=NULL && !(*fmin<=frmin)) *fmin=frmin;
      if(fmax!=NULL && !(*fmax>=frmax)) *fmax=frmax;
    }
    fi++;
  } /* next frame */
  imgEmpty(&img);
  if(ret==STATUS_NOMATRIX && fi>0) return STATUS_OK;
  else return ret;
}
/******************************************************************************/
 
/******************************************************************************/
int imgSmoothMax(
  IMG *img,
  float *maxvalue,
  IMG_PIXEL *p 
) {
  int pi, yi, xi, fi;
  float f, v;
 
  if(img->status<IMG_STATUS_OCCUPIED) return(1);
  if(maxvalue==NULL && p==NULL) return(2);
  if(img->dimt<1 || img->dimz<1 || img->dimy<3 || img->dimx<3) return(3);
  if(maxvalue!=NULL) *maxvalue=0.0;
  if(p!=NULL) p->x=p->y=p->z=p->f=1;
  f=-1.0E20;
  for(pi=0; pi<img->dimz; pi++)
    for(yi=1; yi<img->dimy-1; yi++)
      for(xi=1; xi<img->dimx-1; xi++)
        for(fi=0; fi<img->dimt; fi++) {
          v=img->m[pi][yi-1][xi-1][fi]+
            img->m[pi][yi-1][xi  ][fi]+
            img->m[pi][yi-1][xi+1][fi]+
            img->m[pi][yi  ][xi-1][fi]+
            img->m[pi][yi  ][xi  ][fi]*2.0+
            img->m[pi][yi  ][xi+1][fi]+
            img->m[pi][yi+1][xi-1][fi]+
            img->m[pi][yi+1][xi  ][fi]+
            img->m[pi][yi+1][xi+1][fi];
          v*=0.1;
          if(v>f) {
            f=v; if(p!=NULL) {p->x=xi+1; p->y=yi+1; p->z=pi+1; p->f=fi+1;}}
        }
  if(maxvalue!=NULL) *maxvalue=f;
  return(0);
}
/******************************************************************************/
 
/******************************************************************************/
int imgGetPeak(
  IMG *img,
  float beforeTime,
  IMG_PIXEL *p,
  int verbose
) {
  int zi, yi, xi, fi, mf;
  float f;
 
  if(verbose>0) printf("imgGetPeak(img, %g, p, %d)\n", beforeTime, verbose);
  if(img->status<IMG_STATUS_OCCUPIED) return(1);
  if(p==NULL) return(2);
  if(img->dimt<1 || img->dimz<1 || img->dimy<1 || img->dimx<1) return(3);
  if(beforeTime<img->mid[0]) {
    if(verbose>0) fprintf(stderr, "Error: invalid max search time setting.\n");
    return(4);
  }
  f=nanf(""); mf=img->dimt; p->x=p->y=p->z=p->f=1;
  for(zi=0; zi<img->dimz; zi++) {
    for(yi=0; yi<img->dimy; yi++) {
      for(xi=0; xi<img->dimx; xi++) {
        for(fi=0; fi<img->dimt; fi++) if(img->mid[fi]<=beforeTime) {
          if(!isfinite(img->m[zi][yi][xi][fi])) continue;
          if(isfinite(f) && img->m[zi][yi][xi][fi]<f) // lower
            continue;
          if(isfinite(f) && img->m[zi][yi][xi][fi]==f) { // equal
            // only use this if earlier than in prev max
            if(fi>=mf) continue;
          } 
          f=img->m[zi][yi][xi][fi];
          p->x=xi+1; p->y=yi+1; p->z=zi+1; p->f=fi+1;
          mf=fi;
        }
      }
    }
  }
  if(!isfinite(f)) return(5);
  if(verbose>2) printf("maxval := %g\n", f);
  return(0);
}
/******************************************************************************/
 
/******************************************************************************/
int imgGetMaxTime(
  IMG *img,
  IMG *mimg,
  const int w,
  int verbose
) {
  if(verbose>0) printf("imgGetMaxTime(*img, *mimg, %d)\n", w);
  
  if(img==NULL || img->status!=IMG_STATUS_OCCUPIED) return(1);
  if(mimg==NULL) return(2);
  if(mimg->status==IMG_STATUS_OCCUPIED) imgEmpty(mimg);
 
  /* 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(mimg, img->dimz, img->dimy, img->dimx, 1);
  if(ret) return(ret);
  /* set image header information */
  imgCopyhdr(img, mimg);
  mimg->start[0]=img->start[0]; mimg->end[0]=img->end[img->dimt-1];
  mimg->mid[0]=(mimg->start[0]+mimg->end[0])/2.0;
  mimg->unit=CUNIT_UNKNOWN;
 
  if(w==0) {
    for(int zi=0; zi<img->dimz; zi++) {
      for(int yi=0; yi<img->dimy; yi++) {
        for(int xi=0; xi<img->dimx; xi++) {
          /* Find the frame with max value */
          int ti=0, mi=0; 
          double mv=img->m[zi][yi][xi][ti];
          for(ti=1; ti<img->dimt; ti++) {
            if(img->m[zi][yi][xi][ti]<mv) continue;
            mi=ti; mv=img->m[zi][yi][xi][ti];
          }
          if(mv>0.0) mimg->m[zi][yi][xi][0]=img->mid[mi];
          else mimg->m[zi][yi][xi][0]=0.0;
        }
      }
    }
    return(0);
  }
 
  if(w==1) { // Same as the equation for mean residence time (MRT) for PTACs in pharmacokinetics
    for(int zi=0; zi<img->dimz; zi++) {
      for(int yi=0; yi<img->dimy; yi++) {
        for(int xi=0; xi<img->dimx; xi++) {
          /* Compute the value weighted time */
          double sumw=0.0, sumt=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=1.0;
            sumt+=img->m[zi][yi][xi][ti]*img->mid[ti]*fdur;
            sumw+=img->m[zi][yi][xi][ti]*fdur;
          }
          sumt/=sumw;
          if(sumt>0.0 && sumw>0.0) mimg->m[zi][yi][xi][0]=sumt;
          else mimg->m[zi][yi][xi][0]=0.0;
        }
      }
    }
    return(0);
  }
 
  if(w>1) {
    for(int zi=0; zi<img->dimz; zi++) {
      for(int yi=0; yi<img->dimy; yi++) {
        for(int xi=0; xi<img->dimx; xi++) {
          /* Find the frame with max value */
          int ti=0, mi=0; 
          double mv=img->m[zi][yi][xi][ti];
          for(ti=1; ti<img->dimt; ti++) {
            if(img->m[zi][yi][xi][ti]<mv) continue;
            mi=ti; mv=img->m[zi][yi][xi][ti];
          }
          /* calculate weighted mean from subsequent frames, if possible */
          if(mi<1 || mi>img->dimt-2) continue;
          int i1, i2; i1=mi-1; i2=mi+1; if(i1>0 && i2<img->dimt-1) {i1--; i2++;}
          double sumw=0.0, sumt=0.0;
          for(int i=i1; i<=i2; i++) {
            if(!(img->m[zi][yi][xi][i]>0.0)) continue;
            sumt+=img->m[zi][yi][xi][i]*img->mid[i];
            sumw+=img->m[zi][yi][xi][i];
          }
          sumt/=sumw;
          if(sumt>0.0) mimg->m[zi][yi][xi][0]=sumt;
          else mimg->m[zi][yi][xi][0]=0.0;
        }
      }
    }
    return(0);
  }
 
  return(0);
}
/******************************************************************************/
 
/******************************************************************************/
int imgGetMaxFrame(
  IMG *img,
  IMG *mimg,
  int verbose
) {
  if(verbose>0) printf("imgGetMaxFrame()\n");
  
  if(img==NULL || img->status!=IMG_STATUS_OCCUPIED) return(1);
  if(mimg==NULL) return(2);
  if(mimg->status==IMG_STATUS_OCCUPIED) imgEmpty(mimg);
 
  /* 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(mimg, img->dimz, img->dimy, img->dimx, 1); if(ret) return(ret);
  /* set image header information */
  imgCopyhdr(img, mimg);
  mimg->start[0]=img->start[0]; mimg->end[0]=img->end[img->dimt-1];
  mimg->mid[0]=(mimg->start[0]+mimg->end[0])/2.0;
 
  /* Go through every pixel */
  int ti, zi, xi, yi;
  double mv; int mi;
  for(zi=0; zi<img->dimz; zi++) {
    for(yi=0; yi<img->dimy; yi++) for(xi=0; xi<img->dimx; xi++) {
      ti=mi=0; mv=img->m[zi][yi][xi][ti];
      for(ti=1; ti<img->dimt; ti++) {
        if(img->m[zi][yi][xi][ti]<mv) continue;
        mi=ti; mv=img->m[zi][yi][xi][ti];
      }
      if(mv>1.0E-008) mimg->m[zi][yi][xi][0]=1.0+mi;
      else mimg->m[zi][yi][xi][0]=0.0;
    }
  }
 
  return(0);
}
/******************************************************************************/
 
/******************************************************************************/
int imgAvg(
  IMG *img,
  IMG_RANGE *r,
  float *avg
) {
  int zi, yi, xi, fi;
  long long n=0;
 
  if(img->status<IMG_STATUS_OCCUPIED) return(1);
  if(r!=NULL) {
    if(r->z1<1 || r->y1<1 || r->x1<1 || r->f1<1) return(2);
    if(r->z2<r->z1 || r->y2<r->y1 || r->x2<r->x1 || r->f2<r->f1) return(3);
    if(r->z2>img->dimz || r->y2>img->dimy || r->x2>img->dimx || r->f2>img->dimt) return(4);
  }
  if(avg==NULL) return(5);
 
  *avg=0.0;
  if(r!=NULL) {
    for(zi=r->z1-1; zi<r->z2; zi++) {
      for(yi=r->y1-1; yi<r->y2; yi++) {
        for(xi=r->x1-1; xi<r->x2; xi++) {
          for(fi=r->f1-1; fi<r->f2; fi++) if(!isnan(img->m[zi][yi][xi][fi])) {
            *avg+=img->m[zi][yi][xi][fi]; n++;
          }
        }
      }
    }
  } else {
    for(zi=0; zi<img->dimz; zi++) {
      for(yi=0; yi<img->dimy; yi++) {
        for(xi=0; xi<img->dimx; xi++) {
          for(fi=0; fi<img->dimt; fi++) if(!isnan(img->m[zi][yi][xi][fi])) {
            *avg+=img->m[zi][yi][xi][fi]; n++;
          }
        }
      }
    }
  }
  if(n>0) *avg/=(float)n;
  return(0);
}
/******************************************************************************/
 
/******************************************************************************/
float f_kth_smallest(
  float *data,
  long long int n,
  long long int k
) {
  long long int i, j, l, m;
  float x, s;
 
  l=0; m=n-1;
  while(l<m) {
    x=data[k]; i=l; j=m;
    do {
      while(data[i]<x) i++;
      while(x<data[j]) j--;
      if(i<=j) {s=data[i]; data[i]=data[j]; data[j]=s; i++; j--;}
    } while(i<=j);
    if(j<k) l=i;
    if(k<i) m=j;
  }
  return(data[k]);
}
/******************************************************************************/
 
/******************************************************************************/
float fmedian(
  float *data,
  long long int n
) {
  long long int k;
  float d1, d2;
 
  if(n<1) return(0.0);
  if(n%2) {
    k=(n-1)/2; return(f_kth_smallest(data, n, k));
  } else {
    k=n/2; d1=f_kth_smallest(data, n, k-1); d2=f_kth_smallest(data, n, k);
    return(0.5*(d1+d2));
  }
}
/******************************************************************************/
 
/******************************************************************************/
float fmean(
  float *data,
  long long int n,
  float *sd
) {
  long long int i;
  float sumsqr=0.0, sqrsum=0.0, avg;
 
  if(n<1 || data==NULL) {if(sd!=NULL) *sd=0.0; return(0.0);}
 
  for(i=0; i<n; i++) {sumsqr+=data[i]*data[i]; sqrsum+=data[i];}
  avg=sqrsum/(float)n; if(sd==NULL) return(avg);
  if(n==1) {
    *sd=0.0;
  } else {
    sqrsum*=sqrsum;
    *sd=sqrt( (sumsqr - sqrsum/(float)n) / (float)(n-1) );
  }
  return(avg);
}
/******************************************************************************/
 
/******************************************************************************/
void fMinMaxFin(
  float *data,
  long long int n,
  float *fmin,
  float *fmax
) {
  if(fmin!=NULL) *fmin=nanf("");
  if(fmax!=NULL) *fmax=nanf("");
 
  long long int i;
  for(i=0; i<n; i++) if(isfinite(data[i])) break;
  if(i==n) return; // no finite values found
  float mi, ma;
  mi=ma=data[i++];
  for(; i<n; i++) if(isfinite(data[i])) {
    if(data[i]>ma) ma=data[i];
    else if(data[i]<mi) mi=data[i];
  }
  if(fmin!=NULL) *fmin=mi;
  if(fmax!=NULL) *fmax=ma;
  return;
}
/******************************************************************************/
 
/******************************************************************************/