//****************************************************************************
//
// Matlab C routine file:  CC4
//
// Written 4/99 by N. Howe.
//
//****************************************************************************

#include "mex.h"

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

#ifndef MIN
#define MIN(a, b) (((a) < (b)) ? (a) : (b))
// macro for computing minimum
#endif

#ifndef MAX
#define MAX(a, b) (((a) > (b)) ? (a) : (b))
// macro for computing maximum
#endif

//****************************************************************************
//
// Gateway driver to call the calculation from Matlab.
//
// This is the Matlab entry point
//
//****************************************************************************

void 
mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) {
  int i, x, y;
  int nrow, ncol, ncomp, nfinal, id, id1, id2;
  int *equiv, *final, *xmpx, *xmpy;
  double *arr, *comp, *outncomp, *outxmpx, *outxmpy;
	
  // check for proper number of arguments
  if (nrhs != 1) {
    mexErrMsgTxt("Exactly one argument required.");
  }
  if (nlhs > 3) {
    mexErrMsgTxt("Too many output arguments.");
  }

  // check format of arguments
  if (!mxIsDouble(prhs[0])||mxIsComplex(prhs[0])) {
    mexErrMsgTxt("Inputs must be type double.");
  }
  nrow = mxGetM(prhs[0]);
  ncol = mxGetN(prhs[0]);
  arr = mxGetPr(prhs[0]);

  // allocate output space
  plhs[0] = mxCreateDoubleMatrix(nrow, ncol, mxREAL);
  comp = mxGetPr(plhs[0]);

  // process
  ncomp = 0;
  for (y = 0; y <ncol; y++) {
    ncomp++;
    comp[y*nrow] = (double)ncomp;
    for (x = 1; x < nrow; x++) {
      if (arr[y*nrow+x] != arr[y*nrow+x-1]) {
        ncomp++;
      }
      comp[y*nrow+x] = (double)ncomp;
    }
  }
  equiv = (int*)mxMalloc(ncomp*sizeof(int));
  final = (int*)mxMalloc(ncomp*sizeof(int));
  equiv--;	// make 1-indexed
  final--;	// make 1-indexed
  for (i = 1; i <= ncomp; i++) {
    equiv[i] = i;
  }
  for (y = 1; y <ncol; y++) {
    for (x = 0; x < nrow; x++) {
      if (arr[y*nrow+x] == arr[(y-1)*nrow+x]) {
        id1 = (int)comp[y*nrow+x];
        id2 = (int)comp[(y-1)*nrow+x];
        while (equiv[id1] != id1) {
          id1 = equiv[id1];
        }
        while (equiv[id2] != id2) {
          id2 = equiv[id2];
        }
        id = MIN(id1, id2);
        equiv[MAX(id1, id2)] = id;
        id1 = (int)comp[y*nrow+x];
        while (equiv[id1] != id) {
          id2 = id1;
          id1 = equiv[id1];
          equiv[id2] = id;
        }
        id1 = (int)comp[(y-1)*nrow+x];
        while (equiv[id1] != id) {
          id2 = id1;
          id1 = equiv[id1];
          equiv[id2] = id;
        }
      }
    }
  }
  nfinal = 0;
  for (i = 1; i <= ncomp; i++) {
    id = i;
    while (equiv[id] != id) {
      id = equiv[id];
    }
    equiv[i] = id;
    if (id == i) {
      nfinal++;
      final[i] = nfinal;
    }
  }

  // allocate space for example indices if requested
  if (nlhs == 3) {
    plhs[1] = mxCreateDoubleMatrix(1, nfinal, mxREAL);
    outxmpy = mxGetPr(plhs[1]);
    plhs[2] = mxCreateDoubleMatrix(1, nfinal, mxREAL);
    outxmpx = mxGetPr(plhs[2]);
    outxmpx--;  // make 1-indexed
    outxmpy--;  // make 1-indexed
    for (i = 1; i <= nfinal; i++) {
      outxmpx[i] = outxmpy[i] = 0;
    }
  }

  // compose final labels
  for (i = 0; i < nrow*ncol; i++) {
    comp[i] = final[equiv[(int)comp[i]]];
    if ((nlhs > 2)&&(outxmpx[(int)(comp[i])] == 0)) {
      outxmpx[(int)(comp[i])] = i/nrow+1;
      outxmpy[(int)(comp[i])] = i%nrow+1;
    }
  }

  // create output
  if (nlhs == 2) {
    plhs[1] = mxCreateDoubleMatrix(1, 1, mxREAL);
    outncomp = mxGetPr(plhs[1]);
    *outncomp = (double)(nfinal);
  }

  // free stuff
  mxFree(equiv+1);
  mxFree(final+1);
}

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