proj/mag/ident/libs/mexfunctions/mrf_segment_wts.c

#include "mex.h"  /* must appear first */
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <memory.h>
#include <limits.h>
#include "mexhead.h" /* my mex defines */
#include "rng.h"     /* random numbers */
#include "Doc/mrf_segment_wts_docstring.h"  /* autogenerated from this file */


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

%mrf_segment_wts: segment with a discrete Markov random field
%
% [yp,post] = mrf_segment_wts(iter,T,beta,alpha,dist,y,lprob1,...,lprobK)
% * Performs iter full sweeps (using temperature schedule T) of
% Gibbs sampling on an input labeling y, with entries in 1..K, to produce
% an output labeling yp.  Posterior probability is optionally returned.
% Posterior is actually an energy function, which is correct with respect
% to changes in y, but does not have the correct scale factors which
% vary with beta, alpha, and dist.
% * Conditional distributions of pixel data x given class y are 
% calculated externally, and given by lprob1...lprobK, as log-
% probabilities.  
% * Per-class biases are given by alpha, which can be empty, indicating
% no bias.  Otherwise, the interpretation is alpha(k) is the prior 
% log-probability of seeing class k.
% * The smoothness parameter of the MRF (Potts model) is beta.  If beta
% is a matrix, beta(k,l) is the smoothness "reward" given to a site of
% class k for having a neighbor of class l.  The scalar beta thus 
% corresponds to a diagonal matrix with repeated beta entries.
% (beta and alpha agree with definitions in the Besag paper below.)
% * Pixel-pixel distances are given by dist, where dist(nu,m,n) gives
% the distance between pixel (i,j) and its neighbor number nu, looking
% up or left.  For the 3x3 neighborhood, pixel s=(m,n) has 8 neighbors s',
% and distances to 4 of them, where s'<s, at offsets:
%   (-1,-1),(0,-1),(1,-1),(-1,0),
% are given in that order, in dist(i,j,:).  The other neighbors 
% have s'>s, and the corresponding distances are listed in the 
% symmetric entries of dist.  If dist=[], it is taken to be
% everywhere 1, thus removing the direction-sensitive smoothing
% (but still smoothing).
% * Classes are 1..K, but labels NaN and 0 are not updated or counted
% as neighbors.  Any NaN in a log-probability forces a NaN in the
% output class.
% * A `clock' at each pixel may speed computation.  This recognizes that,
% if the neighbors of a pixel do not change, the Gibbs sampler is rolling
% a stationary die once per iteration to determine the pixel label.  
% The waiting time until another label change is then geometric,
% and this waiting time can be sampled once to short-circuit a series
% of die rolls.  See Ripley, below.
% * The speedup of this method is greatest when few labels change.  If
% supplied, iter(2) is the threshold (in [0,1]) of #changed/#labels
% for a switchover from the ordinary method to the clock method.  
% (Only one switch is permitted in the iteration sequence.)  
% If iter(2) = 0, the ordinary coin-flip method is used throughout; 
% if =1, the clock method is used throughout.
% * Annealing is permitted through the `temperature' parameter T. 
% Initial temperature is T(1), reduced by a factor of T(2) at each 
% iteration. The default T(2) = 1 suppresses annealing.
% Clocks are reset whenever the cumulative drop in temperature 
% reaches a factor of T(3).  (Roughly, then, clocks are reset every 
% log(T(3))/log(T(2)) iterations.)  If T(4) is present and equals
% zero, additional iterations are done at zero temperature, at the
% end of the normal annealing schedule, until the labels reach a
% fixed point.  (This is equivalent to Besag's ICM.)
% * If iter[1] has a fractional part, that part is multiplied by 2**31 
% and rounded to set the desired random number seed; if not, a 
% pseudorandom seed is generated.  This is provided to allow repeatability.
% * This is implemented as a MEX file.
% 
% Inputs:
%   int iter[1] or [2] = [0 0.05];
%   real T[0] or [1] or [2] or [3] or [4] = [1 1 0.8 0];
%   real beta[1] or [K,K];
%   real alpha[K] or [0] = [];
%   real dist[nbr,m,n] or [0];
%   int y[m,n];  -- 1 <= y <= K, or 0 or NaN
%   real lprob1[m,n];
%   ...
%   real lprobK[m,n];
% 
% Outputs:
%   int yp[m,n];  -- 1 <= y <= K, or 0 or NaN
%   opt real post;
% 
% See Also:  mrf_segment, makemrfdiscwts
% Geman and Geman, Stochastic relaxation, Gibbs distributions, and
% the Bayesian restoration of images, IEEE PAMI Nov. 1984
% J. Besag, On the statistical analysis of dirty pictures, JRSSB, 1986
% B. D. Ripley, Statistical Inference for Spatial Processes,
% Cambridge U., 1988, p. 99

% turmon sep/oct 2006, weighted distance
% turmon 5 march 1999, modified to streamline for batch operations
% MJT 18 Mar 1996

***************************************************************/
/* Updated by -mex2pymex.py-ver1- on Wed Sep 23 16:55:52 2009 */

/**************************************************************
*** WEIGHTED DISTANCE IMPLEMENTATION NOTES
THIS IS AN EXPERIMENTAL VERSION

07/nov/2006

Still to-do: check the kmax and clock variables eg for a small example
  as mentioned below.
Tests on solar images seem OK, but would like to quantitatively test
  the rate of pixel class flips.


04/nov/2006

Still need to check kmax and clock.

Also, one query: I thought the "most probable class and clock > it_left"
property would imply the number of clocks to manage should be decreasing,
in fact, in small examples, should be decreasing from say 16 to say 1 or 2
active pixels.  Does this actually happen?  I couldn't find evidence.

Perhaps something as simple as a count of active clocks?


03/nov/2006

checked kmax/0temp against earlier version: works
calloc seems OK.

todo 
return clock/kmax to matlab for checks


10/24
Suggest an implementation change:
- instead of resetting all the clocks when temperature ratchets down, 
  don't reset large clocks as temperature drops

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

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

WEIGHTED DISTANCE IMPLEMENTATION METHOD

type Count becomes a quasi-count (counts weighted neighbor classes)
3d array count3 records quasi-counts for each site, each class

When s=(n,m) is altered from past_y->this_y, counts for each neighbor (i,j),
(i,j) in {-1,0,1}^2 \ {0,0}, are NOW updated according to:

  count3[past_y][n+i][m+j]--; 
  count3[this_y][n+i][m+j]++;

This changed to:

  dist = wtd_distance((m,n),(m+j,n+i));
  count3[past_y][n+i][m+j] -= dist; 
  count3[this_y][n+i][m+j] += dist;

and we define wtd_distance(s,s') like:

  if (s<s') {
    nu = offset(i,j)    // 0, 1, 2, or 3
    return dist(nu,m,n,nu) // lookup in array
  } else {
    nu = offset(-i,-j)  // 0, 1, 2, or 3
    return dist(nu,m+j,n+i)
  }

A check is done [routine: check_counts()] verifying that the neighbor
quasi-counts, which we update as we go along, are still correct at
the end of the iteration.  They are real-valued, so the comparison
allows small differences.

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

typedef float Count; /* neighbors' per-class quasi-counts */
typedef unsigned int Clock; /* clock for expiration of label */
typedef unsigned char Class;  /* class label */
#define ClassMax 8        /* maximum number of classes (can be increased) */
#define CLOCK_MAX UINT_MAX   /* max integer held by Clock */
#define NBHD_WIDTH 1 /* width of neighborhood around one pixel */
#define NBHD_BLOCK (NBHD_WIDTH*2+1) /* total block size (odd) */
#define NBHD_NUM_2 ((NBHD_BLOCK*NBHD_BLOCK-1)/2) /* number of nbrs / 2 */

#define NARGIN_MIN_const 7 /* constant to do arithmetic on below */
#define NARGIN_MIN  (NARGIN_MIN_const)    /* min number of inputs (one class) */
#define NARGIN_MAX  (NARGIN_MIN_const+ClassMax-1)      /* max number of inputs */
#define NARGOUT_MIN 0      /* min number of output args */
#define NARGOUT_MAX 2      /* max number of output args */
/*#undef NARGIN_MIN_const /* don't need this any more */

#define ARG_iters 0
#define ARG_T     1
#define ARG_beta  2
#define ARG_alpha 3
#define ARG_dist  4
#define ARG_y     5
#define ARG_lprob 6

#define ARG_yp    0
#define ARG_post  1

static const char *progname = "mrf_segment_wts";
#define PROGNAME mrf_segment_wts
static const char *in_specs[NARGIN_MAX] = {
  "RS|RV(2)",
  "RS|RV(2)|RV(3)|RV(4)",
  "RS(1)|RM",
  "RV",
  "RA", /* dist */
  "RM", /* y */
  "RA|RM", /* lprob1 */
  "RM", /* lprob2 */
  "RM", /* lprob3 */
  "RM", /* lprob4 */
  "RM", /* lprob5 */
  "RM", /* lprob6 */
  "RM", /* lprob7 */
  "RM"};/* lprob8 */
static const char *in_names[NARGIN_MAX] = {
  "iters",
  "T",
  "beta",
  "alpha",
  "dist",
  "y",
  "lprob1",
  "lprob2",
  "lprob3",
  "lprob4",
  "lprob5",
  "lprob6",
  "lprob7",
  "lprob8"};
static const char *out_names[NARGOUT_MAX] = {
  "yp", "post"};

// defined just for brevity in a generated #include file
#define SHORTNAME msw

/*
 * Default arguments (used in mexFunction)
 */
#define Iters_count 2
#define Default_iters_count 2
static double Default_iters[Default_iters_count] = { 0.0, 0.05 };
#define T_count 4
#define Default_T_count 4
static double Default_T[Default_T_count]     = { 1.0, 1.0, 0.8, 0.0 };


/*********************************************************************
 *  UTILITIES: Labels
 *********************************************************************/

/* simple utilities */ 
#define ISLEGAL(l,L) (((l)>=0)&&((l)<L)) /* for indexing: 0 <= l < L */
#define ActiveLabel(y) ((!isnan(y)) && ((y) > 0)) /* not NaN and not 0 */


/*
 * set label to inactive if any log-probability is NaN
 */
static
void
find_inactive_labels(
                 double **y,
                     double ***logp,
             int M,
                 int N,
                 int K)
{
  int m, n, k; /* counters */
  const double nan = mxt_getnand();

  for (n = 0; n < N; n++)
    for (m = 0; m < M; m++) {
      for (k = 0; k < K; k++)
    if (isnan(logp[k][n][m])) { 
      y[n][m] = nan; /* ie, inactive */
      break; /* go to next pixel */
    } 
    } /* for m, n */
}


/*
 * check_labels: Are classes y all integers in 1..K, or 0, or NaN?
 */
static
int
check_labels(
         int K, /* correct number of classes excluding 0 */
         double *y,
         int N)
{
  int n;

  for (n = 0; n < N; n++)
    if (!isnan(y[n]))
      if (y[n] < 0 || y[n] > K || floor(y[n]) != y[n])
    return(0);
  return(1);
}


/*********************************************************************
 *  UTILITIES: Counts, aka Quasi-counts
 *********************************************************************/

/*
 * sites2dist: given two sites, return their distance d(s,s').
 * The distance array is a flat list generically indexed by (nu,m,n), 
 *   0 <= nu < NBHD_NUM_2  (varies fastest)
 *   0 <= m < M            (next fastest)
 *   0 <= n < N            (varies slowest)
 * Sites are:
 *   s = (m,n)
 *   s'= (m+j,n+i)
 * and the distance array is size NBHD_NUM_2xMxN. For example, for 3x3
 * blocks centered about each site, sites have 8 neighbors, 
 * NBHD_NUM_2 = 4 and NBHD_BLOCK = 3.  For 5x5 blocks, sites have
 * 24 neighbors, NBHD_NUM_2 = 12, and NBHD_BLOCK = 5.
 * 
 * The trick is finding an offset into this array.  If the array
 * is null, the return value is always 1, so the distances
 * are isotropic (but not zero) in this case.
 * 
 * The neighbor at (j,i) corresponds to an offset
 *   nu = j + NBHD_BLOCK*i,
 * where j (corresponding to m) varies faster than i 
 * (corresponding to n).  
 * This is where the neighbor indexing is critical.
 * Eg, for del=3, we have 3x3 neighborhoods around the central
 * point [s], and del2half = 4 prior neighbors that are ordered 
 * as shown below.  The subsequent neighbors s' marked by [-] are 
 * not stored with s because s is a prior neighbor of each such
 * s'.  Storage is done this way so that every patch of distances
 * (ie, the numbered squares below) is still in the same indexing
 * order as the original images.  So, for example, if we have all 
 * 9 distances, a simple reshape(dists, [3 3]) will put them in
 * an immediately recognizable format.
 *     n->
 *  m  0 3 -
 *  |  1 s -
 *  V  2 - - 
 * In the general del case, neighbors are indexed by 
 * nu = 0...NBHD_NUM_2-1, and nu=0 corresponds to the "farthest"
 * pixel (lexicographically) and nu=max corresponds to the site
 * closest to s.
 * 
 * 1. If s' < s  in lexicographic order, then the distance is stored
 * in the NBHD_NUM_2-length vector associated with s.  The offset, nu,
 * into the vector is found from s-s' = (-j,-i), which is a pair of
 * nonnegative integers, since i < 0 or j < 0 in this case.
 * 2. If s' > s, then the distance is stored in the vector associated
 * with s'.  The offset is found in the same way, from s'-s=(j,i),
 * which is again a pair of nonnegative integers.
 * This could be implemented as a macro, but doing so seems to offer
 * no speed advantage.
 */
static
double
sites2dist(double *dist,
       int M,
       int N,
       int m,
       int n,
       int j,
       int i)

{
  /* the nu below is a precursor to the nu in the narrative above */
  double nu = j + NBHD_BLOCK*i; /* ie, delta_m + BLOCK*delta_n */
  int offset;

  if (!dist) return 1.0; /* null pointer */
  if (nu < 0) /* case 1 above */
    /* s' < s: return the (j,i) neighbor of s=(m,n) */
    /* nu ranges from -NBHD_NUM_2(max offset)...-1(nearest), 
     * we have to "flip" it by doing NBHD_NUM_2 + nu to get
     * the actual offset */
    offset = (n*M+m)*NBHD_NUM_2 + (NBHD_NUM_2 + nu);
  else
    /* s' > s: return the (j,i) neighbor of s'=(m+j,n+i) */
    /* nu ranges from +1(nearest)...NBHD_NUM_2(max offset) */
    offset = ((n+i)*M+(m+j))*NBHD_NUM_2 + (NBHD_NUM_2 - nu);
  /* range check is disabled unless debug flag is on */
  mxAssert(offset >= 0 && offset < M*N*NBHD_NUM_2, "sites2dist range error");
  return dist[offset];
}


/*
 * Initialize neighborhood quasi-counts based on label structure
 */
static
void
init_counts(
        Count ***count3,
        double **y2,   /* in 1..K, or 0 */
        double *dist,  /* can be NULL */
        int K,
        int M, 
        int N)

{
  int m, n; /* counters */
  int i, j; /* counters */
  int k; /* counters */

  /* Computation */
  for (n = 0; n < N; n++)
    for (m = 0; m < M; m++) {
      /* Find quasi-counts, of each label, over neighbors of site (m,n) */
      for (k = 0; k < K; k++)
    count3[k][n][m] = 0; /* reset all quasi-counts */
      if (!ActiveLabel(y2[n][m])) continue; /* just skip finding its nbrs */
      for (i = -NBHD_WIDTH; i <= NBHD_WIDTH; i++)
    if (ISLEGAL(n+i,N))
      for (j = -NBHD_WIDTH; j <= NBHD_WIDTH; j++) 
        if (ISLEGAL(m+j,M) && ((i != 0) || (j != 0))) {
          if (ActiveLabel(y2[n+i][m+j]))
        count3[(int) (y2[n+i][m+j]-1)][n][m] += 
          sites2dist(dist,M,N,m,n,j,i); ; /* y2 in 1..K */
        } 
    } /* for n, m */
}


/*
 * Does current quasi-count array match the current labeling?
 * Return value: number of mismatches, or (-1) for calloc failure
 */
static
int
check_counts(
         int K, /* correct number of classes excluding 0 */
         Count ***count3,
         double **y2,
         double *dist,   /* can be NULL */
         int M,
         int N)
{
  int m, n; /* counters */
  int i, j; /* counters */
  int k; /* counters */
  Count *count1; /* neighbor-counts for one pixel */
  int retval = 0; /* assume OK */
   
  count1 = (Count *) calloc(K, sizeof(Count));
  if (!count1) return -1; /* trouble */
  /* Computation */
  for (n = 0; n < N; n++)
    for (m = 0; m < M; m++) {
      /* Find counts of neighbors of each type */
      for (k = 0; k < K; k++)
    count1[k] = 0; /* reset all quasi-count accumulators */
      if (!ActiveLabel(y2[n][m])) continue; /* just skip finding its nbrs */
      for (i = -NBHD_WIDTH; i <= NBHD_WIDTH; i++)
    if (ISLEGAL(n+i,N))
      for (j = -NBHD_WIDTH; j <= NBHD_WIDTH; j++) 
        if (ISLEGAL(m+j,M) && ((i != 0) || (j != 0))) {
          if (ActiveLabel(y2[n+i][m+j]))
        count1[(int) (y2[n+i][m+j]-1)] += 
          sites2dist(dist,M,N,m,n,j,i); /* y2 in 1..K */
        } 
      /* Check count3 against count1 (float eps = 1.2e-7) */
      for (k = 0; k < K; k++)
    if (fabs(count1[k] - count3[k][n][m]) > 
        1e-4*max(1.0,fabs(count3[k][n][m]))) {
      /*
      mexPrintf("(m,n,k) = (%d,%d,%d) count1=%g, count3=%g del = %g\n",
            m, n, k, count1[k], count3[k][n][m], 
            count1[k] - count3[k][n][m]); */
      retval++;
    }
    } /* for n, m */
  free(count1);
  return(retval);
}


/*********************************************************************
 *  UTILITIES: Clocks
 *********************************************************************/

/*
 * include routines for clock generation
 */

#include "mrf_newclock.h"


/*
 * Reset the clock associated with each image site.  This is
 * done when temperature drops, which makes all probabilities
 * stale.
 *
 * This is nontrivial. We compare iter_remain to the current clock.  
 * There are three cases:
 *  clock <= iter_remain: reset clock.  The temperature drop has
 *     caused the probabilities to change, and the clock value
 *     is stale.
 *  clock > iter_remain and label is NOT the highest-probability
 *     class: reset clock.  The temperature drop will enhance the
 *     probability of the highest-probability class relative to the
 *     others, with the likely consequence of making the clock
 *     (if figured again) smaller than it currently is.
 *  clock > iter_remain and label IS highest-probability: leave 
 *     clock alone.  The temperature drop will cause the probability
 *     of the current class y (because it is highest) to be even
 *     larger relative to the other classes.  Thus its clock,
 *     if chosen again, would be larger than it now is.  Since it
 *     will expire after updates cease, the class will never change.
 * Because so many labels are dominated by one class, the last case
 * is quite common.  This is true with real data, but also true of
 * situations where the class balance is quite delicate.
 *
 * It is critical that kmax, the most probable class, be maintained 
 * outside this routine.  Both kmax and y are in the range 1..K, with
 * 0 representing "stale."
 * 
 */
static
void
reset_clocks(
         Clock **clock2, /* clocks */
         Class **kmax2,  /* class of maximum probability */
         double **y2,    /* current class */
         int iter_remain,  /* iterations remaining */
         int M, 
         int N)
{
  int m, n; /* counters */

  for (n = 0; n < N; n++)
    for (m = 0; m < M; m++) {
      /* reset clock only when it is now smaller than the number
     of iters left, or its class is not the most likely.
     Otherwise, decreasing temperatures imply iterations will
     run out before its clock does.
     Note, both kmax2 and y2 are in the range 1..K. kmax2 = 0
     indicates that kmax is stale.  Such a value will never equal y2.
      */
      if (!ActiveLabel(y2[n][m])) continue; /* clock irrelevant */
      if (clock2[n][m] <= iter_remain || kmax2[n][m] != y2[n][m])
    clock2[n][m] = 0;
    } /* for n, m */
}

/*
 * initialize both clocks and kmax by setting to zero, indicating
 * stale values for that site.
 */
static
void
init_clocks(
        Clock **clock2, /* clocks */
        Class **kmax2,  /* class of maximum probability */
        int M, 
        int N)
{
  int m, n; /* counters */

  for (n = 0; n < N; n++)
    for (m = 0; m < M; m++)
      clock2[n][m] = kmax2[n][m] = 0; /* for initialization */
}


/*********************************************************************
 * ROUTINES FOR LABELING
 * The first (per pixel MAP) is non-iterative, but all others are 
 * iterative and need an outer loop to invoke them repeatedly
 *********************************************************************/

/*
 * Initialize all labels using pixel-wise MAP estimate
 * Note, we only need the per-class probabilities and alpha,
 * there is no smoothness term.
 */
static
void
per_pixel_map_estimate(
                   double **y,
                       double ***logp,
               double *alpha,
                   int M,
                   int N,
                   int K)
{
   int m, n, k; /* counters */
   double one_prob; /* holds one probability */
   double winner; /* highest prob so far */
   const double nan = mxt_getnand();

   for (n = 0; n < N; n++) 
     for (m = 0; m < M; m++) {
       y[n][m] = 0.0; /* "undefined" to begin with */
       for (k = 0; k < K; k++) {
     one_prob = logp[k][n][m] + alpha[k]; /* data term, plus label prior */
     if (isnan(one_prob)) { 
       /* nan -> invalid pixel; quit looping */
       y[n][m] = nan; 
       break;
     } else if (y[n][m] == 0.0 || one_prob > winner) {
       /* valid, and no winner yet, or new winner */
       winner = one_prob;
       y[n][m] = k+1; /* stays in range 1..K */
     }
       } /* for k */
       /* end loop: y[n][m] != 0 */
     } /* for m, n */
}

/*
 * find_posterior: finds the posterior probability of the current
 * labeling y2 using the stored (and assumed valid) quasi-counts
 * in count3.  Returns this value as a double.
 */
static
double
find_posterior(
           double **y2, /* in 1..K, or 0 or NaN */
           Count  ***count3,
           double ***lprob,
           int K,
           double **beta2,
           double *alpha,
           int M, 
           int N)

{
   int m, n; /* counters */
   int l, k; /* class label */
   double post;      /* accumulate log-posterior */
   double mrf_smooth; /* accumulates smoothness term of MRF */

   post = 0.0;
   for (n = 0; n < N; n++)
     for (m = 0; m < M; m++) {
       /* skip inactive pixels */
       if (!ActiveLabel(y2[n][m])) continue; 
       /* Find posterior of the class label y[n][m] */
       k = y2[n][m]-1; /* in 0..K-1 */
       /* find smoothness term if label #l is eventually chosen */
       for (mrf_smooth = l = 0; l < K; l++)
     mrf_smooth += beta2[l][k] * count3[l][n][m]; /* wtd sum of counts */
       /* three terms: class prior, local smoothness, data */
       post += alpha[k] + mrf_smooth + lprob[k][n][m];
     }
   return post;
}


/*
 * Computational routine: iterate_0temp
 * Performs one image sweep using stored neighbor-counts,
 * at "zero temperature" -- no sampling is used, so at each
 * iteration, each pixel is updated with the most likely label.
 * This corresponds to Besag's ICM algorithm.
 *
 * Returns fraction of active pixels that were changed:
 * this is used to determine when a stationary point has been
 * reached.
 */
static
double
iterate_0temp(
          double **y2, /* in 1..K, or 0 or NaN */
          Count  ***count3,
          Class  **kmax2, 
          double ***lprob,
          int K,
          double **beta2,
          double *alpha,
          double *dist,  /* can be NULL */
          int M, 
          int N)

{
   int m, n; /* counters */
   int i, j; /* counters */
   int k, l; /* counters */
   double lprob1; /* log-probability */
   double lprob1_max; /* maximum log-probability */
   const double lprob_max_init = -1E64; /* anything < log(minfloat) */
   double mrf_smooth; /* accumulates smoothness term of MRF */
   Count d1; /* the current quasi-count */
   int this_y, past_y; /* integers for class labels */
   int num_pel_active = 0, num_pel_flip = 0; /* # active or modified pels */

   for (n = 0; n < N; n++)
     for (m = 0; m < M; m++) {
       /* skip inactive pixels */
       if (!ActiveLabel(y2[n][m])) continue; 
       num_pel_active++; /* count active pixels */
       /* if this label is already most probable, continue */
       if (kmax2[n][m] == y2[n][m]) continue;
       /* Find log-probabilities, and the largest log-probability */
       this_y = 0; /* anything in 0..K-1 is OK, it is overwritten in loop */
       for (lprob1_max = lprob_max_init, k = 0; k < K; k++) {
         /* find smoothness term if label #l is eventually chosen */
     for (mrf_smooth = l = 0; l < K; l++)
       mrf_smooth += beta2[l][k] * count3[l][n][m]; /* wtd sum of counts */
         /* three terms: class prior, local smoothness, data */
     lprob1 = alpha[k] + mrf_smooth + lprob[k][n][m];
     if (lprob1 > lprob1_max) {
       /* update winning class */
       lprob1_max = lprob1;
       this_y = k; /* new label candidate, in 0..K-1 */
     }
       }
       /* Update the label y */
       kmax2[n][m] = this_y+1; /* record most probable label in 1..K */
       past_y = y2[n][m]-1; /* old label. NB! in 0..K-1 */
       if (this_y == past_y) 
     continue; /* no need to update neighbor counts, just go on */
       y2[n][m] = this_y+1; /* new label as stored, for class in 1..K */
       num_pel_flip++; /* label has changed */
       /* Update costs in neighborhood for the change in label y */
       /* printf("[%d,%d] replaces %d -> %d\n", m, n, past_y, this_y); */
       for (i = -NBHD_WIDTH; i <= NBHD_WIDTH; i++)
     if (ISLEGAL(n+i,N))
       for (j = -NBHD_WIDTH; j <= NBHD_WIDTH; j++) 
         if (ISLEGAL(m+j,M))
           if ((i != 0) || (j != 0)) {  /* skip the current site */
         /* update the neighborhood quasi-counts */
         d1 = sites2dist(dist,M,N,m,n,j,i);
         count3[past_y][n+i][m+j] -= d1; 
         count3[this_y][n+i][m+j] += d1;
         if (d1 > 0)
           kmax2[n+i][m+j] = 0;  /* old max-prob is stale */
           }
     }
#ifdef DEBUG
   mexPrintf("0T: Pixels changed = %f = %d/%d\n", 
              (double) num_pel_flip/num_pel_active, 
              (int) num_pel_flip, (int) num_pel_active); 
#endif
   return((double) num_pel_flip/num_pel_active); /* efficiency measure */
}


/*
 * Computational routine: iterate_count
 * Performs one image sweep using stored neighbor-counts
 * more efficient when many pixels are changing.
 *
 * Returns fraction of active pixels that were changed:
 * this is a measure of computational burden.
 */
static
double
iterate_count(
          double **y2, /* in 1..K, or 0 or NaN */
          Count  ***count3,
          double ***lprob,
          int K,
          double T,
          double **beta2,
          double *alpha,
          double *dist,  /* can be NULL */
          int M, 
          int N)

{
  int m, n; /* counters */
  int i, j; /* counters */
  int k, l; /* counters */
  double *prob1, *lprob1; /* probabilities and their logs */
  double lprob1_max; /* maximum log-probability */
  const double lprob_max_init = -1E64; /* anything < log(minfloat) */
  const double Tinv = 1.0/T; /* prefer to multiply */
  double mrf_smooth; /* accumulates smoothness term of MRF */
  Count d1; /* the current quasi-count */
  double sum; /* normalization of probability mass function */
  double uni, rng_uniform(); /* random numbers */
  int this_y, past_y; /* integers for class labels */
  int num_pel_active = 0, num_pel_flip = 0; /* # active or modified pels */

  /* Set up per-color bins to hold probabilities */
  prob1  = (double *) calloc(K, sizeof(double));
  lprob1 = (double *) calloc(K, sizeof(double));
  mxAssert(prob1 && lprob1, "mrf_segment_wts: failed calloc(count/lprob)");  
  /* Computation */
  for (n = 0; n < N; n++)
    for (m = 0; m < M; m++) {
      /* skip inactive pixels */
      if (!ActiveLabel(y2[n][m])) continue; 
      num_pel_active++; /* count active pixels */
      /* Find log-probabilities, and the largest log-probability */
      for (lprob1_max = lprob_max_init, k = 0; k < K; k++) {
    /* find smoothness term if label #l is eventually chosen */
    for (mrf_smooth = l = 0; l < K; l++)
      mrf_smooth += beta2[l][k] * count3[l][n][m]; /* wtd sum of counts */
    /* three terms: class prior, local smoothness, data */
    lprob1[k] = Tinv * (alpha[k] + mrf_smooth + lprob[k][n][m]);
    if (lprob1[k] > lprob1_max) lprob1_max = lprob1[k]; 
      }
      /* compute probability normalization, scaled by exp(-lprob1_max) */
      for (sum = k = 0; k < K; k++)
    sum += (prob1[k] = exp(lprob1[k] - lprob1_max)); /* rescale to max */
      /* Note: 0 <= prob1[k] <= 1 and 1 <= sum <= K.  prob1/sum is the pmf
       * of the new class; prob1 itself is un-normalized. */
      /* Draw from the probability */
      uni = rng_uniform()*sum; /* rng_uniform on [0,sum) */
      for (k = 0; k < K; k++)
    if ((uni -= prob1[k]) <= 0)
      break;
      /* Update the label y */
      past_y = y2[n][m]-1; /* old label. NB! in 0..K-1 */
      y2[n][m] = k+1; /* new label as stored, for class in 1..K */
      this_y = y2[n][m]-1; /* new label. NB! {this,past}_y both in 0..K-1 */
      if (past_y == this_y) /* our state unchanged */
    continue; /* no need to update neighbor counts, just go on */
      num_pel_flip++; /* label has changed */
      /* Update costs in neighborhood for the change in label y */
      /* printf("[%d,%d] replaces %d -> %d\n", m, n, past_y, this_y); */
      for (i = -NBHD_WIDTH; i <= NBHD_WIDTH; i++)
    if (ISLEGAL(n+i,N))
      for (j = -NBHD_WIDTH; j <= NBHD_WIDTH; j++) 
        if (ISLEGAL(m+j,M)) 
          if ((i != 0) || (j != 0)) {  /* skip the current site */
        /* update the neighborhood quasi-counts */
        d1 = sites2dist(dist,M,N,m,n,j,i);
        count3[past_y][n+i][m+j] -= d1; 
        count3[this_y][n+i][m+j] += d1;
          }
    }
  /* printf("Efficiency = %f = %d/%d\n", 
     (double) num_pel_flip/num_pel_active, 
     (int) num_pel_flip, (int) num_pel_active); */
  /* tidy up */
  free(lprob1);
  free(prob1);
  return((double) num_pel_flip/num_pel_active); /* efficiency measure */
}


/*
 * Core computational routine: iterate_clock
 * Performs one image sweep using running clocks and neighbor-quasi-counts.
 * It is most efficient when few pixels are changing.
 *
 * Like the other routines here, we keep track of the quasi-count array, 
 * count3.  This contains weighted distances, cumulative over all 
 * neighbors of a given site, to neighbors of class 1,...,K.
 * We also keep track of kmax, the maximum-probability class number,
 * i.e. argmax_{k=1..K} Pr(y(s) == k | y(s'), s' ~ s)
 * This is in the range 1...K, or 0 if it is unknown.  The latter is
 * true, for instance, if any neighbor's class assignment has 
 * changed, which will make the previously computed probabilities
 * stale.
 *
 * Returns fraction of active pixels whose clock needed to be reset:
 * a measure of computational burden.  This is *not* the same as the 
 * fraction of pixels whose label changed, because sometimes clocks
 * need to be recomputed when neighbor classes change.
 */
static
double
iterate_clock(
          double **y2,    /* in 1..K, or 0 or NaN */
          Clock  **clock2, 
          Class  **kmax2, 
          Count  ***count3,
          double ***lprob,
          int K,
          double T,
          int it_left,   /* iterations remaining, for clock generation */
          double **beta2,
          double *alpha,
          double *dist,  /* can be NULL */
          int M, 
          int N)

{
  int m, n; /* counters */
  int i, j; /* counters */
  int k, l; /* counters */
  double *prob1, *lprob1; /* prob's and log-probs, for one pel */
  double lprob1_max; /* maximum log-probability at this pel */
  const double lprob_max_init = -1E64; /* anything < log(minfloat) */
  const double Tinv = 1.0/T; /* prefer to multiply */
  double mrf_smooth; /* accumulates smoothness term of MRF */
  Count d1; /* the current quasi-count */
  double sum; /* normalization of probability mass function */
  double uni, rng_uniform(); /* random numbers */
  int this_y, past_y; /* integers for class labels */
  int curr_class;     /* class label, or sentinel value */
  int num_pel_active = 0, num_pel_flip = 0; /* # active or modified pels */
  int num_pel_expire = 0; /* # pels with expired clocks */

  /* Set up per-color bins to hold probabilities */
  prob1  = (double *) calloc(K, sizeof(double));
  lprob1 = (double *) calloc(K, sizeof(double));
  mxAssert(prob1 && lprob1, "mrf_segment_wts: failed calloc(clock/lprob)");  
  /* Computation */
  for (n = 0; n < N; n++)
    for (m = 0; m < M; m++) {
      /* skip inactive pixels */
      if (!ActiveLabel(y2[n][m])) continue;
      num_pel_active++; /* count active pixels */
      /* At each pel, there are three cases                                 
       *                                                                    
       * DESCRIPTION     STATE           ACTION                             
       * Clock Ticking   clock > 1       decrement clock                    
       * Clock Expired   clock = 1       draw different state               
       *                 curr_class =    choose new clock                   
       *                    <old state>  reset neighbors' clocks            
       * Clock Reset     clock = 0       draw a state (needn't be different)
       *                 curr_class = -1 choose new clock     
       *                                 reset neighbors if state changed 
       */
      /* Case 1 above: If time remains, just decrement the clock */
      if (clock2[n][m] > 1) {
    clock2[n][m]--;
    continue; 
      } 
      num_pel_flip++; /* need to sample, although label may not change */
      if (clock2[n][m] == 1) num_pel_expire++; /* label must change */
      /* ...otherwise, must re-evaluate probabilities and find a new state */
      /* curr_class holds
     (1) class of this pel (0...K-1) if its clock expired
     (2) the sentinel value (-1) if its clock was reset */
      curr_class = (clock2[n][m] == 0) ? -1 : (y2[n][m]-1);
      /* Find log-probabilities, and the largest log-probability */
      for (lprob1_max = lprob_max_init, k = 0; k < K; k++) {
    /* find smoothness term if label #l is eventually chosen */
    for (mrf_smooth = l = 0; l < K; l++)
      mrf_smooth += beta2[l][k] * count3[l][n][m]; /* wtd sum of counts */
    /* three terms: class prior, local smoothness, data */
    lprob1[k] = Tinv * (alpha[k] + mrf_smooth + lprob[k][n][m]);
    if (lprob1[k] > lprob1_max) {
      lprob1_max = lprob1[k];
      kmax2[n][m] = k+1; /* succeeds for at least one k, in range 1..K */
    }
      }
      /* compute probability normalization, scaled by exp(-lprob1_max) */
      for (sum = k = 0; k < K; k++) {
    prob1[k] = exp(lprob1[k] - lprob1_max); /* rescale to max */
    if (k != curr_class) /* skips summing curr_class if clock expired */
      sum += prob1[k];
      }
      /* Note: 0 <= prob1[k] <= 1; if curr_class == -1 then 1 <= sum <= K but
       * if curr_class >= 0, only know 0 <= sum <= K.  If curr_class == -1,
       * prob1/sum is the pmf of the new class; in this case prob1 itself 
       * is un-normalized.  If curr_class >= 0, the last sentence is true
       * for prob1, excluding its curr_class entry, which is understood to
       * have zero probability. */
      /* Draw from the probability */
      uni = rng_uniform()*sum; /* rng_uniform on [0,sum) */
      for (k = 0; k < K; k++) {
    if (k == curr_class) /* skips curr_class if clock expired, as above */
      continue;
    if ((uni -= prob1[k]) <= 0)
      break;
      }
      /* end loop: note k != curr_class for sure */
      /* Update the label y */
      past_y = y2[n][m]-1; /* old label. NB! in 0..K-1 */
      y2[n][m] = k+1; /* new label as stored, for class in 1..K */
      this_y = y2[n][m]-1; /* new label. NB! {this,past}_y both in 0..K-1 */
      /* Update clock */
      if (past_y == this_y) { /* clock reset && our state unchanged */
    /* must choose new clock because nbr changed => cost changed */
    clock2[n][m] = (Clock) new_clock(prob1[this_y]/sum, it_left);
    /* however, need not reset neighbors' clocks: can continue */
    continue;
      }
      /* Update costs in neighborhood for the change in label y */
      /* printf("[%d,%d] replaces %d -> %d\n", m, n, past_y, this_y); */
      for (i = -NBHD_WIDTH; i <= NBHD_WIDTH; i++)
    if (ISLEGAL(n+i,N))
      for (j = -NBHD_WIDTH; j <= NBHD_WIDTH; j++) 
        if (ISLEGAL(m+j,M)) 
          if ((i != 0) || (j != 0)) { /* skip the current site */
        /* update the neighborhood quasi-counts */
        d1 = sites2dist(dist,M,N,m,n,j,i);
        count3[past_y][n+i][m+j] -= d1; /* -> range 0..K-1 */
        count3[this_y][n+i][m+j] += d1; /* -> range 0..K-1 */
        /* and reset the neighbors' clocks and max-k's */
        if (d1 > 0) {
          clock2[n+i][m+j] = 0;  /* old clock is invalid */
          kmax2[n+i][m+j]  = 0;  /* old max-prob is stale */
        }
          }
      /* Choose new clock for the current site (choose nbr clocks as we go) */
      if (curr_class == -1) 
    /* were not avoiding any class => sum as computed contains all probs */
    clock2[n][m] = (Clock) new_clock(prob1[this_y]/sum, it_left);
      else
    /* were avoiding curr_class => sum excludes prob[curr_class] */
    clock2[n][m] = (Clock) new_clock(prob1[this_y]/(sum+prob1[curr_class]),
                     it_left);
    }
#ifdef DEBUG
  printf("Efficiency = %f = %d/%d, %d expired\n", 
     (double) num_pel_flip/num_pel_active, 
     (int) num_pel_flip, (int) num_pel_active, num_pel_expire);
#endif
  /* tidy up */
  free(lprob1);
  free(prob1);
  return((double) num_pel_flip/num_pel_active); /* efficiency measure */
}



/*********************************************************************
 * OUTER LOOP
 * 
 * Controls and initializes the iteration process
 *
 * Progresses through initialization then through various stages
 * of update which consist of sweeps through all pixels in the
 * image using the stored quasi-counts or clocks.
 *********************************************************************/
  
/* stages determining which routine is used for labeling updates */
#define ITER_INIT  0   /* dummy state for startup */
#define ITER_COUNT 1   /* performing "count" updates */
#define ITER_CLOCK 2   /* performing "clock" updates */
#define ITER_0TEMP 3   /* considering final T=0 updates */
#define ITER_END   4   /* exit state */

static
control(double  **yp2,    /* result: labeling */
    double  *postprob,/* result: posterior probability, NULL ok */
    Clock   **clock2, /* place for per-label timers */
    Class   **kmax2,  /* place for per-pixel most-likely-class */
    Count  ***count3, /* place for neighbor count */
    int K,            /* number of classes */
    double *iters,    /* iteration control */
    double *T,        /* temperature control */
    double **beta2,   /* smoothness */
    double *alpha,    /* per-class bias */
    double *dist,     /* inter-site distances */
    double ***lprob3, /* log-probabilities of data */
    int M,
    int N)            /* size of input and result */

{
   double flipfrac=0;      /* how many pixels flipped in this iteration */
   int stage = ITER_INIT;  /* just initializing */
   double iter_flipthresh; /* threshold controlling method using flipfrac */
   int iter, iter_max;     /* counter and target number of iters */
   double T_iter;          /* current temperature (never increases) */
   double T_rat;           /* temperature reduction ratio */
   double T_step;          /* temperature stepsize */
   double T_thresh = -1;   /* thresholded temperature ("clock" input) */
   double T_0temp;         /* zero if doing 0-temperature iterations */
   int mismatch;           /* number of mismatches in quasi-count check */
   int count3_inited = 0;  /* have we initialized count3 array? */

   /* assume RNG is already seeded */
   /* set up symbolic names for inputs */
   iter_max  = iters[0]; /* target iter number */
   iter_flipthresh = iters[1]; /* threshold on flipfrac for method change */
   T_iter = T[0];        /* current temperature */
   T_rat  = T[1];        /* temperature drop per iter (multiplicative) */
   T_step = T[2];        /* temp change until clock reset */
   T_0temp= T[3];        /* T=0 iterations, or not */
   /* done initializing: set stage up */
   if (iter_max > 0) {
     /* some gibbs sampler iters will be done: start at COUNT or CLOCK */
     if (iter_flipthresh < 1) 
       stage = ITER_COUNT; /* typical: start with iterate_count */
     else
       stage = ITER_CLOCK; /* start with iterate_clock */
   } else {
     /* no g.s. iters: either 0 temperature iters, or none */
     if (T_0temp == 0)
       stage = ITER_0TEMP; /* start with iterate_0temp */
     else
       stage = ITER_END;   /* no iters at all */
   }
   /* initialize auxiliary variables if needed */   
   if (stage < ITER_END || postprob != NULL) {
     /* initialize neighbor count variables */   
     init_counts(count3, yp2, dist, K, M, N);
     count3_inited = 1; /* it's OK to check it later for consistency */
     /* reset the per-site clocks and max-probability class */
     init_clocks(clock2, kmax2, M, N);
   }
   /* Loop over several sweeps, using one method for each sweep */
   for (iter = 0; stage < ITER_END; iter++) {
     /* if/else below guarantees iter < iter_max, else stage is 0TEMP */
     if (iter < iter_max) {
       T_iter *= T_rat; /* ordinary annealing: lower the present temp'ture */
     } else {
       stage = ITER_0TEMP; /* check for T=0 updates, and possibly leave */
       T_iter = 0;         /* unused at this stage, but maintain it */
     }
#ifdef DEBUG
     mexPrintf("I = %03d  S%d (%.6f)  T = %0.3f\n", 
           iter, stage, flipfrac, T_iter);
#endif
     /* perform one sweep, depending on selected method */
     if (stage == ITER_COUNT) {
       /* perform a "count" update */
       flipfrac = iterate_count(yp2, count3, lprob3, 
                K, T_iter, beta2, alpha, dist, M, N);
       if (isnan(flipfrac)) flipfrac = 0; /* 0/0 case -> 0 */
       /* switch update method if few labels change */
       if (flipfrac < iter_flipthresh) stage = ITER_CLOCK;
     } else if (stage == ITER_CLOCK) {
       /* refresh on first clock pass, or if temp has changed too much */
       if (T_thresh < 0 || T_iter < T_thresh*T_step) {
     /* printf("resetting clocks\n"); */
         reset_clocks(clock2, kmax2, yp2, iter_max-iter, M, N);
         T_thresh = T_iter; /* reset temperature threshold also */
       }
       /* perform a "clock" update */
       flipfrac = iterate_clock(yp2, clock2, kmax2, count3, lprob3, 
                K, T_thresh, iter_max-iter,
                beta2, alpha, dist, M, N);
       if (isnan(flipfrac)) flipfrac = 0; /* 0/0 case -> 0 */
     } else {
       /* check to perform a zero-temperature update */
       if (T_0temp != 0) {
     stage = ITER_END; /* no T=0 updates, break out of the loop */
     continue;
       }
       flipfrac = iterate_0temp(yp2, count3, kmax2, lprob3, 
                K, beta2, alpha, dist, M, N);
       if (isnan(flipfrac)) flipfrac = 0; /* 0/0 case -> 0 */
       if (flipfrac == 0)
     stage = ITER_END;  /* no more changes -> exit */
     } /* if stage, else */
   } /* for iter */

   if (postprob != NULL)
     /* posterior was desired: easy to find since we have count3 */
     *postprob = find_posterior(yp2, count3, lprob3, K, beta2, alpha, M, N);
   /* consistency check for robustness */
   if (count3_inited)
     if ((mismatch = check_counts(K, count3, yp2, dist, M, N)) != 0) {
       char errstr[200];
       mexErrMsgTxt((snprintf(errstr, sizeof(errstr),
                  "%s: Internal error: "
                  "%d of %d counts inconsistent", 
                  progname, mismatch, K*M*N), errstr));
     }
}


/*********************************************************************
 * GATEWAY ROUTINE
 *********************************************************************/
#ifdef StaticP  /* undefined under mex */
StaticP
#endif
void 
mexFunction(
        int nlhs, 
        mxArray *plhs[], 
        int nrhs,
        const mxArray *prhs[])
{
  int M, N, K;      /* hold numeric parameters */
  int i, k;         /* counters */
  int lprob_input_3d; /* was log-prob input one 3d matrix or not */
  double *postprob; /* posterior probability, or NULL if unwanted */
  double ***lprob3; /* matrix of data log-probabilities (3d) */
  double **yp2;     /* matrix of labels (2d) */
  Count  ***count3; /* matrix of neighbor quasi-counts (3d) */
  Count    *count;  /* pointer to data referred to by above ptrs */
  Clock   **clock2; /* matrix of label clocks (2d) */
  Clock    *clock;  /* pointer to data referred to by above ptrs */
  Class   **kmax2;  /* matrix of most probable classes (2d) */
  Class    *kmax;   /* pointer to data referred to by above ptrs */
  double *iters;    /* iteration input */
  int rng_seed;     /* the rng seed that was used */
  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 args
   */
  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);

  /* find how big the inputs are */
  lprob_input_3d = (nrhs - ARG_lprob == 1); // single 3d logprob input?
  if (!lprob_input_3d) {
    K = nrhs - ARG_lprob; /* number of classes = number of spare inputs */
    M = mxGetM(prhs[ARG_lprob]);
    N = mxGetN(prhs[ARG_lprob]);
  } else {
    if (mxGetNumberOfDimensions(prhs[ARG_lprob]) != 3)
      mexErrMsgTxt((snprintf(errstr, sizeof(errstr),
             "%s: A single lprob must be 3-D",
                progname), errstr));
    M = mxGetDimensions(prhs[ARG_lprob])[0];
    N = mxGetDimensions(prhs[ARG_lprob])[1];
    K = mxGetDimensions(prhs[ARG_lprob])[2];
  }
  /* crude check -- will give a more informative message than later ones */
  if (K > ClassMax || K < 2)
    mexErrMsgTxt((snprintf(errstr, sizeof(errstr),
               "%s: Expect %d <= K <= %d",
               progname, 2, ClassMax), errstr));
   
  start_sizechecking();
  /* beta */
  sizeinit(prhs[ARG_beta]);
  sizeisMN(1, 1); /* beta = [1,1] */
  sizeisMN(K, K); /* beta = [K,K] */
  sizecheck_msg(progname, in_names, ARG_beta);
  /* alpha */
  sizeinit(prhs[ARG_alpha]);
  sizeisMN(0, 0); /* alpha = [] */
  sizeisMN(1, K); /* alpha = 1xK or Kx1 */
  sizecheck_msg(progname, in_names, ARG_alpha);
  /* dist */
  sizeinit(prhs[ARG_dist]);
  sizeisMN(0, 0); /* dist = [] (2d) */
  sizeis3(0, 0, 0); /* dist = [] (3d) */
  sizeisMN(NBHD_NUM_2*M, N); /* allow flattened dist for now */
  sizeis3(NBHD_NUM_2, M, N); /* dist = NBR x M x N */
  sizecheck_msg(progname, in_names, ARG_dist);
  /* check sizes of log-probability matrices */
  if (!lprob_input_3d) /* i.e., if supplied individually */
    for (k = 1; k < K; k++) {
      sizeinit(prhs[ARG_lprob]);
      sizeagree(prhs[ARG_lprob+k]); /* size(cost) == size(x) */
      sizecheck_msg(progname, in_names, ARG_lprob);
    }
  /* check size of y */
  sizeinit(prhs[ARG_y]);
  sizeisMN(0, 0); /* y = [] */
  sizeisMN(M, N); /* y = [M,N] */
  sizecheck_msg(progname, in_names, ARG_y);

  /*
   * create space for output
   */

  /* make the space */
  plhs[ARG_yp] = mxCreateDoubleMatrix(M, N, mxREAL);
  mxAssert(plhs[ARG_yp], "mrf_segment_wts: failed mxCreate (output y)");
  setrange(plhs[ARG_yp], 0.0, 0.0); /* disables scaling */
  yp2 = mxt_make_matrix2(plhs[ARG_yp], -1, -1, 0.0); /* 2d indexing into yp */
  if (nlhs > ARG_post) {
    plhs[ARG_post] = mxCreateDoubleMatrix(1, 1, mxREAL);
    mxAssert(plhs[ARG_post], "mrf_segment_wts: failed mxCreate (output post)");
    setrange(plhs[ARG_post], 0.0, 0.0); /* disables scaling */
    postprob = mxGetPr(plhs[ARG_post]);
  } else {
    postprob = NULL; /* sentinel for not wanted */
  }

  /*
   * Establish index arrays for nonstandard matrices
   */
  /* 1: set up each of K probability maps for 2d indexing */
  if (lprob_input_3d) {
    /* this is freed with one call to mxFree */
    lprob3 = mxt_matrix3_index(mxGetPr(prhs[ARG_lprob]), M, N, K);
    mxAssert(lprob3, "mrf_segment_wts: failed calloc (lprob3)");
  } else {
    lprob3 = (double ***) calloc(K, sizeof(double **)); /* submatrix ptrs */
    mxAssert(lprob3, "mrf_segment_wts: failed calloc (lprob3)");
    for (k = 0; k < K; k++) {
      lprob3[k] = mxt_make_matrix2(prhs[ARG_lprob + k], -1, -1, 0.0);
      mxAssert(lprob3[k], "mrf_segment_wts: failed calloc (lprob3[k])");
    }
  }
  /* 2: set up count as a 3d matrix */
  /*    (cannot use mxt_make_matrix2, which is for double) */
  count  = (Count *)   calloc(K*M*N, sizeof(Count)); /* all counts */
  count3 = (Count ***) calloc(K, sizeof(Count **));  /* ptrs to submat's */
  mxAssert(count && count3, "mrf_segment_wts: failed calloc (count)");  
  for (k = 0; k < K; k++) {
    count3[k] = (Count **) calloc(N, sizeof(Count *)); /* column indexes */
    mxAssert(count3[k], "mrf_segment_wts: failed calloc (count3[k])");  
    for (i = 0; i < N; i++) 
      count3[k][i] = count + k*(M*N) + i*M; /* set up col ptrs */
  }
  /* 3: set up clock as a 2d matrix (again, not of doubles) */
  clock  = (Clock *)  calloc(N*M, sizeof(Clock)); /* all clocks */
  clock2 = (Clock **) calloc(N,   sizeof(Clock *)); /* 2d indexes */
  mxAssert(clock && clock2, "mrf_segment_wts: failed calloc (clock)");  
  for (i = 0; i < N; i++)
    clock2[i] = clock + i*M; /* set up column pointers */

  /* 4: set up kmax as a 2d matrix (again, not of doubles) */
  kmax  = (Class *)  calloc(N*M, sizeof(Class)); /* all clocks */
  kmax2 = (Class **) calloc(N,   sizeof(Class *)); /* 2d indexes */
  mxAssert(kmax && kmax2, "mrf_segment_wts: failed calloc (clock)");  
  for (i = 0; i < N; i++)
    kmax2[i] = kmax + i*M; /* set up column pointers */

  /*
   * Initialize yp from y or lprob
   */
  /* propagate y -> yp, or initialize yp from lprob */
  if (mxGetN(prhs[ARG_y]) == 0) {
    /* initialize yp from log-probabilities */
    per_pixel_map_estimate(yp2, lprob3, 
               mxt_make_vector(prhs[ARG_alpha], K, 
                       (double *) NULL, 0), 
               M, N, K);
  } else {
    /* first, make sure IC is valid */
    if (!check_labels(K, mxGetPr(prhs[ARG_y]), M*N))
      mexErrMsgTxt((snprintf(errstr, sizeof(errstr),
                 "%s: Classes are 0, NaN, or 1..%d", 
                 progname, K), errstr));
    /* then copy IC to new labeling area */
    memcpy((void *) &(yp2[0][0]), 
       (void *) mxGetPr(prhs[ARG_y]), 
       M * N * sizeof(double));
    /* propagate any NaN's in prob-maps through to the labeling */
    find_inactive_labels(yp2, lprob3, M, N, K);
  }

  /*
   * Set up rng seed
   */
  iters = mxt_make_vector(prhs[ARG_iters], 
              Iters_count, Default_iters, Default_iters_count);
  /* note, if below is 0, no seed was given */
  rng_seed = (int) rint(2147483648.0 * fmod(iters[0], 1.0)); 
  /* initialize the RNG in one of two ways */
  if (rng_seed == 0)
    rng_seed = rng_init_state(); /* varies from run to run */
  else
    rng_init_state_seeded(rng_seed); /* deterministic */

  /*
   * do the computation
   */
  /* enter main loop */
  control(yp2,
      postprob, /* posterior output, if desired */
      clock2,   /* pointer structure set up above */
      kmax2,    /* pointer structure set up above */
      count3,   /* pointer structure set up above */
      K,
      iters,    /* found above */
      mxt_make_vector(prhs[ARG_T], T_count, Default_T, Default_T_count),
      mxt_make_matrix1(prhs[ARG_beta], K, K, 1.0),
      mxt_make_vector(prhs[ARG_alpha], K, (double *) NULL, 0), 
      mxGetNumberOfElements(prhs[ARG_dist]) == 0 ? 
        ((double * ) NULL) :
        mxGetPr(prhs[ARG_dist]),
      lprob3,   /* pointer structure set up above */
      M, N);
  /* ensure classes are 1..K, 0, or NaN */
  if (!check_labels(K, &(yp2[0][0]), M*N))
    mexErrMsgTxt((snprintf(errstr, sizeof(errstr),
               "%s: Internal error: labels out of bounds",
               progname), errstr));
  /* free the temporary space up */
  mxFree(yp2); /* pointer vector allocated with mxt_mm2 */
  free(kmax);
  free(kmax2);
  free(clock);
  free(clock2);
  free(count);
  for (k = 0; k < K; k++) 
    free(count3[k]); /* allocated with calloc */
  free(count3);
  if (lprob_input_3d) {
    mxFree(lprob3); /* allocated with mxt_matrix3_index */
  } else {
    for (k = 0; k < K; k++) 
      mxFree(lprob3[k]); /* allocated with mxt_mm2 */
    free(lprob3);  /* allocated with calloc */
  }
}



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

---