avwstats.c

/* {{{ Copyright etc. */

/*  avwstats.c - basic image stats measurements

    Stephen Smith, FMRIB Image Analysis Group

    Copyright (C) 1999-2002 University of Oxford  */

/*  CCOPYRIGHT */

/* }}} */
/* {{{ defines, includes and typedefs */

#include "libss/libss.h"
#include "libss/libavw.h"

/* }}} */
/* {{{ usage */

void usage(void)
{
  printf("\nUsage: avwstats <input> [options]\n\n");
  printf("-l <lthresh> : set lower threshold\n");
  printf("-u <uthresh> : set upper threshold\n");
  printf("-r           : output <robust min intensity> <robust max intensity>\n");
  printf("-R           : output <min intensity> <max intensity>\n");
  printf("-e           : output mean entropy ; mean(-i*ln(i))\n");
  printf("-v           : output <voxels> <volume>\n");
  printf("-V           : output <voxels> <volume> (for nonzero voxels)\n");
  printf("-m           : output mean\n");
  printf("-M           : output mean (for nonzero voxels)\n");
  printf("-s           : output standard deviation\n");
  printf("-S           : output standard deviation (for nonzero voxels)\n");
  printf("-w           : output smallest ROI <xmin> <xsize> <ymin> <ysize> <zmin> <zsize> <tmin> <tsize> containing nonzero voxels\n");
  printf("-x           : output co-ordinates of maximum voxel\n\n");
  printf("Note - thresholds are not inclusive ie lthresh<allowed<uthresh\n\n");
  exit(1);
}

/* }}} */
/* {{{ main */

int main(argc, argv)
  int   argc;
  char  *argv [];
{
  /* {{{ vars */

image_struct im;
int i;

/* }}} */

  if (argc<3)
    usage();

  avw_read(argv[1],&im);

  im.thresh=im.min=im.max=0; /* flag up that the thresholds need finding */

  for (i = 2; i < argc; i++) {
    if (!strcmp(argv[i], "-l"))
      /* {{{ set lower threshold */

{
  i++;
  if (argc<i+1)
    {
      printf("Error: no value given following -l\n");
      usage();
    }

  im.lthresh=atof(argv[i]);
  im.thresh=im.min=im.max=0; /* flag up that the thresholds need re-finding */
}

/* }}} */
    else if (!strcmp(argv[i], "-u"))
      /* {{{ set upper threshold */

{
  i++;
  if (argc<i+1)
    {
      printf("Error: no value given following -u\n");
      usage();
    }

  im.uthresh=atof(argv[i]);
  im.thresh=im.min=im.max=0; /* flag up that the thresholds need re-finding */
}

/* }}} */
    else if (!strcmp(argv[i], "-r"))
      /* {{{ robust range */

{
  if (im.thresh==im.min)
    find_thresholds(&im,0.1);

  printf("%f %f ",(double)im.thresh2,(double)im.thresh98);
}

/* }}} */
    else if (!strcmp(argv[i], "-R"))
      /* {{{ complete range */

{
  if (im.thresh==im.min)
    find_thresholds(&im,0.1);

  printf("%f %f ",(double)im.min,(double)im.max);
}

/* }}} */
    else if (!strcmp(argv[i], "-e"))
      /* {{{ mean entropy */

/* warning - hist uses >=lthresh etc not >lthresh */

#define HISTOGRAM_BINS 1000

{
  int hist[HISTOGRAM_BINS], j, count;
  double entropy=0;

  if (im.thresh==im.min)
    find_thresholds(&im,0.1);

  count=find_histogram(&im,hist,HISTOGRAM_BINS);

  for(j=0;j<HISTOGRAM_BINS;j++)
    entropy -= ((double)hist[j]/count) * log((double)hist[j]/count);

  printf("%f ",entropy/log((double)HISTOGRAM_BINS));
}

/* }}} */
    else if (!strcmp(argv[i], "-v"))
      /* {{{ voxels/volume between thresholds */

{
  int j, count=0;

  for(j=0;j<im.x*im.y*im.z*im.t;j++)
    if ( (im.i[j]>im.lthresh) && (im.i[j]<im.uthresh) )
      count++;

  printf("%d %f ",count,count*im.xv*im.yv*im.zv);
}

/* }}} */
    else if (!strcmp(argv[i], "-V"))
      /* {{{ nonzero voxels/volume between thresholds */

{
  int j, count=0;

  for(j=0;j<im.x*im.y*im.z*im.t;j++)
    if ( (im.i[j]>im.lthresh) && (im.i[j]<im.uthresh) && (im.i[j]!=0) )
      count++;

  printf("%d %f ",count,count*im.xv*im.yv*im.zv);
}

/* }}} */
    else if (!strcmp(argv[i], "-m"))
      /* {{{ mean */

{
  int j, count=0;
  float sum=0;

  for(j=0;j<im.x*im.y*im.z*im.t;j++)
    if ( (im.i[j]>im.lthresh) && (im.i[j]<im.uthresh) )
      {
	sum+=im.i[j];
	count++;
      }

  if (count>0)
    printf("%f ",sum/count);
  else
    printf("0 ");
}

/* }}} */
    else if (!strcmp(argv[i], "-M"))
      /* {{{ nonzero mean */

{
  int j, count=0;
  float sum=0;

  for(j=0;j<im.x*im.y*im.z*im.t;j++)
    if ( (im.i[j]>im.lthresh) && (im.i[j]<im.uthresh) && (im.i[j]!=0) )
      {
	sum+=im.i[j];
	count++;
      }

  if (count>0)
    printf("%f ",sum/count);
  else
    printf("0 ");
}

/* }}} */
    else if (!strcmp(argv[i], "-s"))
      /* {{{ standard deviation */

{
  int j, count=0;
  float sum=0, sumsq=0;

  for(j=0;j<im.x*im.y*im.z*im.t;j++)
    if ( (im.i[j]>im.lthresh) && (im.i[j]<im.uthresh) )
      {
	sum+=im.i[j];
	sumsq+=im.i[j]*im.i[j];
	count++;
      }

  if (count>0)
    {
      float tmpf = (sumsq - sum * sum / count) / count;
      if (tmpf>0) tmpf=sqrt(tmpf);
      else        tmpf=0;
      printf("%f ",tmpf);
    }
  else
    printf("0 ");
}

/* }}} */
    else if (!strcmp(argv[i], "-S"))
      /* {{{ standard deviation */

{
  int j, count=0;
  float sum=0, sumsq=0;

  for(j=0;j<im.x*im.y*im.z*im.t;j++)
    if ( (im.i[j]>im.lthresh) && (im.i[j]<im.uthresh) && (im.i[j]!=0) )
      {
	sum+=im.i[j];
	sumsq+=im.i[j]*im.i[j];
	count++;
      }

  if (count>0)
    {
      float tmpf = (sumsq - sum * sum / count) / count;
      if (tmpf>0) tmpf=sqrt(tmpf);
      else        tmpf=0;
      printf("%f ",tmpf);
    }
  else
    printf("0 ");
}

/* }}} */
    else if (!strcmp(argv[i], "-x"))
      /* {{{ co-ordinates of max voxel */

{
  int x,y,z,t, xm=0,ym=0,zm=0,tm=0;
  float maxi=im.i[0];

  for(t=0;t<im.t;t++) for(z=0;z<im.z;z++) for(y=0;y<im.y;y++) for(x=0;x<im.x;x++)
    if ( (im.i[t*im.x*im.y*im.z+z*im.x*im.y+y*im.x+x]>im.lthresh) &&
	 (im.i[t*im.x*im.y*im.z+z*im.x*im.y+y*im.x+x]<im.uthresh) &&
	 (im.i[t*im.x*im.y*im.z+z*im.x*im.y+y*im.x+x]>maxi) )
      {
	maxi=im.i[t*im.x*im.y*im.z+z*im.x*im.y+y*im.x+x];
	xm=x; ym=y; zm=z; tm=t;
      }

  if (t<2)
    printf("%d %d %d ",xm,ym,zm);
  else
    printf("%d %d %d %d ",xm,ym,zm,tm);
}

/* }}} */
    else if (!strcmp(argv[i], "-w"))
      /* {{{ nonzero ROI */

{
  int x,y,z,t,xmin=im.x-1,xmax=0,ymin=im.y-1,ymax=0,zmin=im.z-1,zmax=0,tmin=im.t-1,tmax=0;

  for(t=0;t<im.t;t++) for(z=0;z<im.z;z++) for(y=0;y<im.y;y++) for(x=0;x<im.x;x++)
    if ( (im.i[t*im.x*im.y*im.z+z*im.x*im.y+y*im.x+x]>im.lthresh) &&
	 (im.i[t*im.x*im.y*im.z+z*im.x*im.y+y*im.x+x]<im.uthresh) &&
	 (im.i[t*im.x*im.y*im.z+z*im.x*im.y+y*im.x+x]!=0) )
      {
	if (x<xmin) xmin=x;
	if (x>xmax) xmax=x;
	if (y<ymin) ymin=y;
	if (y>ymax) ymax=y;
	if (z<zmin) zmin=z;
	if (z>zmax) zmax=z;
	if (t<tmin) tmin=t;
	if (t>tmax) tmax=t;
      }

  printf("%d %d %d %d %d %d %d %d ",xmin,1+xmax-xmin,ymin,1+ymax-ymin,zmin,1+zmax-zmin,tmin,1+tmax-tmin);
}

/* }}} */
    else
      usage();
  }

  printf("\n");
  return(0);
}

/* }}} */