proj/mag/harp/libs/matlab/mfile-mex/hmi-mask-patch/concomponent.c

#include "mex.h"  /* must appear first */
#include <stdio.h>
#include <string.h>
#include <math.h>
#include "mexhead.h" /* my mex defines */
#include "Doc/concomponent_docstring.h"  /* autogenerated from this file */

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

%concomponent: identify connected components
% 
% [y,nr]=concomponent(x,nbr);
% * identifies connected components of a binary labeling x, storing
% the result as a like-sized image y.  The number of regions is
% optionally returned in nr.
% * On-region components of y are in the range 1..nr.  Inputs x 
% of 0 or NaN are treated as off-region, but passed intact through 
% to output region map as 0 or NaN.
% * Either the 4-pixel neighborhood (N/S/E/W) template can be used,
% or the 8-pixel neighborhood, which includes diagonal pixels.
% Default is 8.
% * This is implemented as a MEX file.
% 
% Inputs:
%   real x(m,n); -- 0, 1, or NaN
%   opt int nbr = 8; -- 4 or 8
% 
% Outputs:
%   int y(m,n); -- 0, NaN, or 1..nr
%   opt int nr;
% 
% See Also: region_bb

% turmon may 2001

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

/* standard boilerplate */
                   
#define NARGIN_MIN  1      /* min number of inputs */
#define NARGIN_MAX  2      /* max number of inputs */
#define NARGOUT_MIN 0      /* min number of output args */
#define NARGOUT_MAX 2      /* max number of output args */

#define ARG_x   0   /* x = data   */
#define ARG_nbr 1   /* nbr = neighbor template */
#define ARG_y   0   /* y = labels */
#define ARG_nr  1   /* region count */

static const char *progname = "concomponent";
#define PROGNAME concomponent
static const char *in_specs[NARGIN_MAX] = {
  "RM",
  "IS"};
static const char *in_names[NARGIN_MAX] = {
  "x",
  "nbr"};
static const char *out_names[NARGOUT_MAX] = {"y", "nr"};

// this is just for brevity when compiled as a C library
#define SHORTNAME ccp

/******************* Utilities ***************/

/* 
 * squeeze out unused region numbers from region map
 * returns final number of regions
 * This means recycling labels that turn out not to be used because
 * they are equivalent to other labels.
 * Algorithm is: scan the region table T, looking for minimal
 * regions (T[r] == r).  Such regions are the prototype for
 * their class.  When one is found, replace all entries in
 * the equivalence table equal to r (including entry r itself)
 * with the current count of the number of minimal regions r_now.
 * This amounts to relabeling all final outputs r with r_now.
 */

static
int
squeeze_regions(int *T,   /* table */
        int Nr)   /* number of valid entries */

{
  int r, rp;     /* region counters */
  int r_now;     /* current region */

  for (r_now = r = 1; r <= Nr; r++)
    if (T[r] == r) {
      /* current region is its own minimal region.
       * renumber it, and all regions mapping to it, 
       * as r_now.
       * Actually, only do the renumbering if
       * it would result in a numbering change.
       */
      if (r_now != r)
    for (rp = r; rp <= Nr; rp++)
      if (T[rp] == r)
        T[rp] = r_now;
      r_now++; /* update current region */
  }
  return r_now-1;
}


/*
 * add the equivalance r1 == r2 to the table T
 * Before this equivalence was added, we knew that (r1,T[r1]) were
 * equivalent, and also that (r2,T[r2]) were equivalent.  The added
 * fact that (r1,r2) are equivalent means that all four states
 * are equivalent.  Their minimal class is the min z of the earlier
 * two minimal classes T[r1], T[r2].  
 * The way to signal the new equivalence is to replace each 
 * instance of r1, r2, T[r1], or T[r2] in the minimal class list 
 * by z.  In some cases, we would know that no instance of one of 
 * these would occur, but it is better to not separate into cases.
 */
static
void
add_equivalence(int *T,
        int Nr,
        int r1,
        int r2)
{
  int z;      /* new minimal region */
  int s1, s2; /* old minimal regions */
  int r;      /* region counter */

  if (r1 == r2) return; /* trivial: nothing to do */
  /* if r1, r2 already have same minimal class listed in T, 
     they are equivalent and there is nothing to update */
  if (T[r1] == T[r2]) return; 
  /* printf("\tadding %4d <-> %4d\t(%4d %4d)\n", r1, r2, T[r1], T[r2]); */
  /* the new minimal region for the equivalence class containing 
     r1 and r2 is the smaller of the previous minimal regions */
  z = T[r1] < T[r2] ? T[r1] : T[r2];
  /* save these old minimal regions */
  s1 = T[r1];
  s2 = T[r2];
  /* replace old instances of the older minimal regions 
   Before this equivalence, we knew that (r1,T[r1]) */
  /* NB: must go up to and including entry #Nr because table starts at 1 
     However, can start at z rather than 1 because no table entry
     lower than z can be numbered r1, r2, s1, or s2, since z is 
     smaller than them all: r < z => T[r] < z => T[r] < r1,r2,s1,s2. */
  for (r = z; r <= Nr; r++)
    if (T[r] == r1 ||
    T[r] == r2 ||
    T[r] == s1 ||
        T[r] == s2) 
      T[r] = z;
}


/* 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 */

/*
 * concomponent: label connected components
 * - We assume x is in the range 1..K, or 0, or NaN.  If 0 or NaN, the 
 * corresponding pixel is excluded from the calculation, and its value
 * is propagated through.  Otherwise, connected components are labeled
 * starting at 1.
 * - The algorithm is to scan through all pixels, putting a region number 
 * in the corresponding entry of y if x is active.  The current region 
 * number is l_new.  If the pixel is isolated (no previously-scanned
 * neighbors were labeled) then a fresh region number is chosen and l_new 
 * is incremented.  Otherwise, a neighbor's region number is used 
 * ("recycling")
 * - Either the 4 or 8 neighbor template may be used.  In case of the
 * 4-nbr template, two neighbors have been examined before the current
 * pixel, so we need to look at those two to recycle labels.  In the 8-
 * nbr template, we look at the four prior pixels.
 * - The rub: there is more than one neighbor, so the neighbors' labels
 * can conflict.  So we see that, in the first scan, we are really 
 * identifying equivalence classes of labels.  We keep a table T to keep
 * track of the equivalences.
 * - The table has one entry per active region-label.  Regions start at 1 
 * so the physical table has one entry more than the number of active
 * labels.  The table starts out at an arbitrary size and is grown 
 * if necessary, this is checked whenever a fresh region number is
 * chosen.
 * - The fundamental property of the table is: the entry for region 
 * r contains the smallest-numbered ("minimal") region equivalent to r.  
 * This makes T[r] unique.  It implies that T[r] <= r.  Equality means
 * that r is itself minimal -- such a region is the prototype for
 * its class since it was the first component uncovered in the scan.
 * Since T[r] is minimal, T[r1] == T[r2] iff r1 and r2 are
 * in the same equivalence class.
 * - T is initialized with T[r] = r, and whenever a neighbor-conflict
 * is discovered, an equivalence is added to T (see below).
 * - After the first scan, a preliminary labeling (correct up to
 * equivalences in T) is determined.  Because we want the final
 * region numbers not to have gaps, we have to squeeze out the
 * "blank" entries in T.  (E.g., if region 2 maps to region 1, then
 * region 2 will not appear in the output and region 3 (say) might
 * need to be renumbered 2.)
 * - Finally, the pixels are scanned in order and their numbers
 * updated to reflect the equivalences: just setting all labels
 * in an equivalence class to the number of the minimal class:
 * y[n][m] = T[y[n][m]].
 */
static
void
concomponent(
         double **x,      /* input image */
         double **y,      /* result: labeling */
         double *n_reg,   /* result: number of regions */
         int n_nbr,       /* number of neighbors */
         int M,           /* dims of image, labeling */
         int N)

{
  void *mxCalloc();
  /* offsets to neighbors that were already scanned in the order we use
     (-1,0)=up; (0,-1)=left; (-1,-1)=left+up; (1,-1)=left+down
     note, by stopping after the first two neighbors, the 4-pixel
     neighbor template (up/left/down/right) can be used painlessly. */
  int delta_m[] = {-1,0,-1,1};  /* previous neighbor offsets */
  int delta_n[] = {0,-1,-1,-1}; /* previous neighbor offsets */
  int m, n;      /* loop counters */
  int nbr;       /* loop counter  */
  int mp, np;    /* position of neighbor pixel */
  int l_new = 0; /* current label counter */
  int l_use;     /* label to use */
  int r;         /* region counter */
  double *box;   /* convenience pointer */
  int *T;        /* minimal region table */
  int r_max = 1024; /* starting point for #regions in table */

  /* initialize table of minimal regions */
  T = (int *) mxCalloc(r_max, sizeof(*T));
  for (r = 0; r < r_max; r++) T[r] = r; /* regions are self-linked */
  /* Step 1:
   * loop over all pixels, finding contiguous regions 
   */
  for (n = 0; n < N; n++) 
    for (m = 0; m < M; m++) {
      /* skip non-object pixels (0 or NaN) */
      if (!ActiveLabel(x[n][m])) {
    y[n][m] = x[n][m]; /* NB: 0/NaN values are copied through */
    continue;
      }
      /* else, are on an object: must label it */
      /* first, search neighbors to recycle labels and find conflicts */
      /* 4-nbr template -> 2 prior neighbors, 8-nbr template -> 4 priors */
      l_use = 0; /* assume un-labeled */
      for (nbr = 0; nbr < n_nbr/2; nbr++) {
    /* neighboring pixel coordinates */
    mp = m + delta_m[nbr];
    np = n + delta_n[nbr];
    /* check the neighbor label */
    if (!ISLEGAL(mp,M) || !ISLEGAL(np,N) || !ActiveLabel(y[np][mp]))
      continue; /* skip bad indexes and unlabeled pixels */
    /* at this point, the neighbor pixel is valid and on-object */
    /* 1: locate an already-labeled neighbor */
    if (!l_use) {
      l_use = y[np][mp]; /* this will be the pixel label */
      continue; /* must continue, to find conflicting labels */
    }
    /* 2: search for label conflicts */
    if (l_use != y[np][mp]) {
      /* pixel label conflicts with another adjacent label:
         enter the conflict as a region equivalence in the table T. 
         This table currently has l_new entries. */
      /* printf("(%d,%d): %d <-> %g\n", n, m, l_use, y[np][mp]); */
      add_equivalence(T, l_new, l_use, (int) y[np][mp]);
      continue;
    }
      }
      /* printf("(%d,%d): <- %d\n", n, m, l_use); */
      /* assign pixel label */
      if (l_use)
    y[n][m] = l_use; /* already found label */
      else {
    y[n][m] = ++l_new; /* create new label */
    /* check if table should be lengthened.  Since the table starts at
       one, and its next entry will be l_new+1, we require that
       T[l_new+1] be valid. */
    if ((l_new + 1) >= r_max) {
      /* must expand region table */
      int *Tp; /* temporary for new table */
      r_max *= 4; /* increase size */
      Tp = (int *) mxCalloc(r_max, sizeof(*T));
      for (r = 0; r < r_max; r++) Tp[r] = r; /* start with T[r] = r */
      /* copy former table to new table */
      memcpy((void *) Tp, (void *) T, (r_max/4) * sizeof(int));
      mxFree(T); /* dump old table */
      T = Tp; /* install new table */
    }
      }
    }
  /* Step 2:
   * squeeze out unused region numbers from region map
   */
  *n_reg = (double) squeeze_regions(T, l_new);
  /* Step 3:
   * loop over all pixels, updating labels 
   */
  for (n = 0; n < N; n++)
    for (m = 0; m < M; m++) {
      if (ActiveLabel(y[n][m]))
    y[n][m] = T[(int) (y[n][m])];
    }
  mxFree(T);
  return;
}
   

/*
 * Gateway routine
 */
#ifdef StaticP  /* undefined under mex */
StaticP
#endif
void 
mexFunction(
        int nlhs, 
        mxArray *plhs[], 
        int nrhs,
        const mxArray *prhs[])
{
  int m,  n;         /* size of x */
  mxArray *plhs_bb;  /* slot to hold bb matrix structure */
  double nr_result;  /* slot for nr number */
  double **x2, **y2; /* 2d indexing */
  int nr_wanted;     /* is the nr output wanted */
  int n_reg;         /* number of regions found */
  int n_nbr;         /* neighbor template */
  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);

  /* 
   * optional args 
   */
  n_nbr = (nrhs <= ARG_nbr) ? 8 : mxGetScalar(prhs[ARG_nbr]);
  if (n_nbr != 4 && n_nbr != 8)
    mexErrMsgTxt((snprintf(errstr, sizeof(errstr),
             
              "%s: Expect neighbor template = 4 or 8, got %d",
              progname, n_nbr), errstr));

  /*
   * create space for output
   */
  /* 1: find how big it must be */
  m  = mxGetM(prhs[ARG_x]);
  n  = mxGetN(prhs[ARG_x]);

  /* 2: make the space */
  /* first, for y */
  /* paradoxically works even if nlhs == 0, in which case strictly
      speaking plhs[0] need not exist (i.e. plhs == NULL in this case)
      MATLAB seems to always have plhs non-null to allow for transmission
      of the result via 'ans' */
  plhs[ARG_y] = mxCreateDoubleMatrix(m, n, mxREAL); 
  /* then, for nr */
  nr_wanted = (nlhs == (ARG_nr+1)); /* if nr output arg supplied */

  /*
   * do the computation
   */
  /* set up ordinary 2-d pointers */
  x2 = mxt_make_matrix2(prhs[ARG_x], -1, -1, 0.0), 
  y2 = mxt_make_matrix2(plhs[ARG_y], -1, -1, 0.0),
  concomponent(x2, y2, 
           /* region number */
           &nr_result, 
           /* sizes */
           n_nbr,
           m, n);
  mxFree(x2);
  mxFree(y2);

  /* plug nr in to its slot if wanted (else, there's no slot!) */
  if (nr_wanted) {
    plhs[ARG_nr] = mxCreateDoubleMatrix(1, 1, mxREAL); 
    *(mxGetPr(plhs[ARG_nr])) = nr_result;
  }
}


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

---