proj/mag/harp/libs/matlab/mfile-mex/standalone/cstat.c

#include "mex.h"  /* must appear first */
#include <stdio.h>
#include <math.h>
#include <string.h>
#include <limits.h> /* INT_MAX */
#include "mexhead.h"
#include "Doc/cstat_docstring.h"  /* autogenerated from this file */

/**************************************************************

%cstat: compute class statistics
%
% [mu,sigma,ct] = cstat(image, label, K, mode)
% * For each `label' value from 1 to K, the mean, standard 
% deviation, and count of the corresponding image pixels are found.  
% NaNs in `image' are skipped.
% * Usable labels are integers from 1..K.  Labels <1, >= K+1, or NaN
% are skipped.  Fractional labels are rounded down.
% * K is the number of classes.  If omitted or given as empty, 
% K is taken to be the maximum value of the labeling.
% * If mode contains 'no-sigma', sigma is not computed.  If mode 
% contains 'no-mu', mu is not computed.  'no-mu' implies 'no-sigma'.
% * If image is given as empty, it is equivalent to 'no-mu', so only
% counts are computed.
% * If mode contains 'unscaled', mu is not normalized by ct, and sigma
% is not computed.
% * When un-computed output arguments must be returned, as in
%   [mu,sigma,ct] = cstat(image, label, [], 'no-mu')
% they are returned as empty.
%
% Inputs:    
%  real image[m,n] or []
%  int  label[m,n]
%  opt  int K = max(label(:))
%  opt  string mode = ''
%
% Outputs:
%  real mu[1,K] or []
%  real sigma[1,K] or []
%  int  ct[1,K]
%

% turmon 7 Apr 2008 enhanced, cleaned up, better error detection
% turmon mar 2011 added modes

***************************************************************/
/* Updated by -mex2pymex.py-ver1- on Fri Sep 25 00:00:48 2009 */

/* constants and globals used for error checking */
#define NARGIN_MIN  2
#define NARGIN_MAX  4
#define NARGOUT_MIN 0
#define NARGOUT_MAX 3

#define ARG_IMAGE  0
#define ARG_LABEL  1
#define ARG_K      2
#define ARG_MODE   3

#define ARG_MU     0
#define ARG_SIGMA  1
#define ARG_CT     2

static const char *progname = "cstat";
#define PROGNAME cstat
static const char *in_specs[NARGIN_MAX] = {
  "RM",
  "IM",
  "IS",
  "SV"};
static const char *in_names[NARGIN_MAX] = {
  "image",
  "label",
  "K",
  "mode"};
static const char *out_names[NARGOUT_MAX] = {
  "mu",
  "sigma",
  "ct"};

/*
 * extracts various mode parameters from a free-form string
 *
 * returns 0 if the string was OK, an explanatory tag otherwise
 */
static
char *
extract_mode(char *mode, int *Munscaled, int *Mnomu, int *Mnosigma) 
{
  char *word;
  char *sep = ", \t";

  if (!mode) return "null"; /* should not happen */
  for (word = strtok(mode, sep); word; word = strtok(NULL, sep)) {
    /* unscaled */
    if      (strcasecmp(word, "unscaled") == 0) { *Munscaled = *Mnosigma = 1; }
    /* sesw, or transposed */
    else if (strcasecmp(word, "no-mu"   ) == 0) { *Mnomu = *Mnosigma = 1;}
    else if (strcasecmp(word, "no-sigma") == 0) { *Mnosigma = 1; }
    /* unrecognized word: not OK */
    else return word; 
  }
  return 0; /* OK */
}

/* compute the minimum and maximum intensity */

static
void
computeminmax(double *input, mwSize n, double *minP, double *maxP)
{
  mwSize i;  
  double min, max;

  // signals not-initialized-yet
  min = max = mxt_getnand(); 
  // skip to first non-nan (ok for n == 0)
  for (i = 0; i < n && isnan(min); i++)
    min = max = input[i];
  // now, ordinary scan from here
  for (; i < n; i++) {
    // nans will fail these if conditions
    if (input[i] < min) min = input[i];
    if (input[i] > max) max = input[i];
  }   
  *minP = min;
  *maxP = max;
}

// take double in 1..K, round down, convert to int in 0..K-1
#define CLASS2INDEX(x) (((int) floor(x)) - 1)

/* 
 * compute the number of valid pixels in each class 
 * (unused)
 */

static
void
computesize(double *image, double *class, const mwSize n,
        double *size, const int nclasses)
{
  mwSize i;
  int index;
  
  /* initialize */
  for (index = 0; index < nclasses; index++)
     size[index] = 0;
  /* accumulate */
  for (i = 0; i < n; i++) {  
    /* NB, check class as well as image */
    if (isnan(class[i])) continue;
    if (image && isnan(image[i])) continue;
    /* Get and check class value */
    index = CLASS2INDEX(class[i]);
    if (index >= 0 && index < nclasses)
      size[index]++;
  }   
}

/* 
 * compute the per-class mean and counts
 *
 * if "image" is NULL, mean is not referenced
 */
static
void
computemean(double *image, 
        double *class, 
        mwSize n,
        double *size,          // set by this routine
        double *mean,          // set by this routine
        int nclasses,
        int unscaled)
{
  mwSize i;
  int index;

  /* initialize */
  for (index = 0; index < nclasses; index++) {
    size[index] = 0;
    if (image) 
      mean[index] = 0;
  }
  /* accumulate */
  for (i = 0; i < n; i++) {
    /* check class and image */
    if (isnan(class[i])) continue;
    if (image && isnan(image[i])) continue;
    /* Get and check class value */
    index = CLASS2INDEX(class[i]);
    if (index >= 0 && index < nclasses) {
      size[index] += 1.0;
      if (image)
    mean[index] += image[i]; // allow null image
    }
  }
  /* normalize */
  if (image && !unscaled)
    for (index = 0; index < nclasses; index++) 
      mean[index] /= size[index]; /* 0/0 goes to NaN */
}

/* per-class standard deviation.
 *
 * have left (mu,sigma) as two-pass for robustness.
 * note that image must be non-null in this routine
 */

static
void
computestd(double *image, 
       double *class, 
       mwSize n, 
       double *mean, 
       double *size, 
       double *std, 
       int nclasses)
{
  mwSize i;
  int index;
  
  /* initialize */
  for (index = 0; index < nclasses; index++)
     std[index] = 0;
  /* accumulate sum of squares */
  for (i = 0; i < n; i++) {  
    /* Check class as well as image */
    if (isnan(class[i]) || isnan(image[i])) continue;
    /* Get and check class value */
    index = CLASS2INDEX(class[i]);
    if (index >= 0 && index < nclasses) {
      std[index] += (image[i] - mean[index]) * (image[i] - mean[index]);
    }
  }   
  /* normalize */
  for (index = 0; index < nclasses; index++) {
    if (size[index] <= 1)
      std[index] = mxt_getnand(); /* not enough data */
    else
      std[index] = sqrt(std[index]/(size[index]-1));
  }   
}

/* gateway routine */

#ifdef StaticP  /* undefined under mex */
StaticP
#endif
void
mexFunction(int nlhs, 
        mxArray *plhs[], 
        int nrhs,
        const mxArray *prhs[])
{
  int nclasses;
  mwSize npix;
  double junk, nclasses_dbl;
  char *mode, *word;           // label for converted mode string and substring
  int empty_image;             // image = [] ?
  int mode_unscaled, mode_nomu, mode_nosigma;
  mxArray *arg_ct, *arg_mu, *arg_sigma;
  char errstr[256];

  /* Hook for introspection (function signature, docstring) */
  if (nrhs < 0) { 
    plhs[0] = mxt_PackSignature((mxt_Signature) (-nrhs), 
                NARGIN_MIN, NARGIN_MAX, 
                NARGOUT_MIN, NARGOUT_MAX, 
                in_names, in_specs, out_names, docstring);
    return;
  }
  /* check number of arguments */
  if ((nrhs < NARGIN_MIN) || (nrhs > NARGIN_MAX))
    mexErrMsgTxt((snprintf(errstr, sizeof(errstr),
               "%s: Expect %d <= input args <= %d",
               progname, NARGIN_MIN, NARGIN_MAX), errstr));
  if ((nlhs < NARGOUT_MIN) || (nlhs > NARGOUT_MAX))
    mexErrMsgTxt((snprintf(errstr, sizeof(errstr),
               "%s: Expect %d <= output args <= %d",
               progname, NARGOUT_MIN, NARGOUT_MAX), errstr));
  mexargparse(nrhs, prhs, in_names, in_specs, NULL, progname);
  start_sizechecking();
  if (mxGetNumberOfElements(prhs[ARG_IMAGE]) > 0) {
    /* if size(image) > 0, size(image) == size(label) */
    sizeinit(prhs[ARG_IMAGE]);
    sizeagree(prhs[ARG_LABEL]); 
    sizecheck_msg(progname, in_names, ARG_IMAGE);
  }

  npix = mxGetNumberOfElements(prhs[ARG_LABEL]);

  /* determine the number of classes */
  if (nrhs <= ARG_K || mxGetNumberOfElements(prhs[ARG_K]) == 0) {
    /* get from labeling */
    computeminmax(mxGetPr(prhs[ARG_LABEL]), npix, &junk, &nclasses_dbl);
  } else {
    nclasses_dbl = mxGetScalar(prhs[ARG_K]);
  }
  /* check # classes */
  nclasses = (int) nclasses_dbl;
  if (isnan(nclasses_dbl) || (nclasses_dbl != (double) nclasses) || (nclasses < 0))
    mexErrMsgTxt((snprintf(errstr, sizeof(errstr),
               "%s: bad supplied/inferred K: %g",
               progname, nclasses_dbl), errstr));

  /* 
   * `mode' input
   */
  empty_image = (mxGetNumberOfElements(prhs[ARG_IMAGE]) == 0);
  mode_unscaled = 0;
  mode_nomu     = empty_image;
  mode_nosigma  = empty_image;
  if (nrhs > ARG_MODE) {
    // OK not to specify it
    if ((mode = mxArrayToString(prhs[ARG_MODE])) == NULL)
      mexErrMsgTxt((snprintf(errstr, sizeof(errstr),
                 "%s: bad mode (non-string?). "
                 "Could not convert mode arg to string.", 
                 progname), errstr));
    if ((word = extract_mode(mode, &mode_unscaled, &mode_nomu, &mode_nosigma)) != NULL)
      mexErrMsgTxt((snprintf(errstr, sizeof(errstr),
                 "%s: bad word <%s> in mode <%.80s>", 
                 progname, word, /* need a fresh copy of the mode */
                 mxArrayToString(prhs[ARG_MODE])), errstr));
    mxFree(mode); // done with the mode string
  }
  // unscaled -> mean is invalid -> sigma will be nonsense
  // no mean computed -> cannot get sigma
  // note that empty_image already implies nomu and nosigma
  if (mode_unscaled || mode_nomu)
    mode_nosigma = 1;

  /* allocate memory for results */
  arg_ct = mxCreateDoubleMatrix(1, nclasses, mxREAL); // always need to keep these
  if (mode_nomu)
    arg_mu = mxCreateDoubleMatrix(0, 0, mxREAL); // return it empty
  else
    arg_mu = mxCreateDoubleMatrix(1, nclasses, mxREAL);
  if (mode_nosigma)
    arg_sigma = mxCreateDoubleMatrix(0, 0, mxREAL); // return it empty
  else
    arg_sigma = mxCreateDoubleMatrix(1, nclasses, mxREAL);
  if (!arg_ct || !arg_mu || !arg_sigma) 
    mexErrMsgTxt((snprintf(errstr, sizeof(errstr),
               "%s: Could not allocate space for results",
               progname), errstr));

  /* class means */
  computemean(mode_nomu ? NULL : mxGetPr(prhs[ARG_IMAGE]), 
          mxGetPr(prhs[ARG_LABEL]), 
          npix,
          mxGetPr(arg_ct),
          mxGetPr(arg_mu),
          nclasses,
          mode_unscaled);
  if (!mode_nosigma) {
    /* class standard deviation */
    computestd(mxGetPr(prhs[ARG_IMAGE]), 
           mxGetPr(prhs[ARG_LABEL]), 
           npix,
           mxGetPr(arg_mu),
           mxGetPr(arg_ct),
           mxGetPr(arg_sigma),
           nclasses); 
  }

  /*
   * Return results
   */
  // mu -- as output 1, it is always returned
  plhs[ARG_MU] = arg_mu;
  // sigma
  if (nlhs > ARG_SIGMA)
    plhs[ARG_SIGMA] = arg_sigma;
  else
    mxDestroyArray(arg_sigma);
  // ct
  if (nlhs > ARG_CT)
    plhs[ARG_CT] = arg_ct;
  else
    mxDestroyArray(arg_ct);
}



/* Hook for generic tail matter */
#ifdef MEX2C_TAIL_HOOK
#include "mex2c_tail.h"
#endif

---