proj/mag/harp/libs/matlab/mfile-mex/assignment/assignment.c

#include "mex.h"  /* must appear first */
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include "mexhead.h"
#include "Doc/assignment_docstring.h"  /* autogenerated from this file */
#include "graphtypes.h"

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

%assignment: solve linear sum assignment problem
%
% [pi1,pi2,s]=assignment(rho)
% * Solves linear-sum assignment problem with *integer* costs
% given in rho.
% * Outputs the maximizing permutation pi1, its inverse pi2,
% and the maximizing cost s.
%
% Inputs:    
%  int rho(m,n)
%
% Outputs:
%  int pi1(m)
%  int pi2(n)
%  int s
%
% See Also: 

% turmon apr 2008
% original authors: Benjamin J. Bornstein, Becky Castano

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

/* constants and globals used for error checking */
#define NARGIN_MIN  1      /* min number of inputs */
#define NARGIN_MAX  1      /* max number of inputs */
#define NARGOUT_MIN 1      /* min number of outputs */
#define NARGOUT_MAX 3      /* max number of outputs */

#define ARG_rho     0
#define ARG_pi1     0
#define ARG_pi2     1
#define ARG_s       2

#define PROGNAME assignment
static const char *progname = "assignment";
static const char *in_specs[NARGIN_MAX] = {
  "IM"};
static const char *in_names[NARGIN_MAX] = {
  "rho"};
static const char *out_names[NARGOUT_MAX] = { 
  "pi1",
  "pi2",
  "score" };


static
Graph 
translateMatrix2Graph(int *size, const mxArray *matrix, int *total_weight);


#ifdef StaticP  /* undefined under mex */
StaticP
#endif
void
mexFunction(int nlhs, 
        mxArray *plhs[], 
        int nrhs,
        const mxArray *prhs[])
{
  char errstr[256];
  int i, j;
  int M, N;  /* sizes of input array */
  int score;
  int size;
  int weight;
  int *Mate;
  Graph graph;
  double *pi1;
  double *pi2;
  Edge  p      = NULL;

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


  /* find highest weighted matching */
  graph = translateMatrix2Graph(&size, prhs[ARG_rho], &weight);
  Mate = Weighted_Match(graph, 1, 1);

  /* set up outputs */
  M = mxGetM(prhs[ARG_rho]);
  N = mxGetN(prhs[ARG_rho]);
  plhs[ARG_pi1] = mxCreateDoubleMatrix(1, M, mxREAL);
  plhs[ARG_pi2] = mxCreateDoubleMatrix(1, N, mxREAL);
  plhs[ARG_s]   = mxCreateDoubleMatrix(1, 1, mxREAL);
  pi1 = mxGetPr(plhs[ARG_pi1]);  /* maps 1:M to 1:N or 0 */
  pi2 = mxGetPr(plhs[ARG_pi2]);  /* maps 1:N to 1:M or 0 */

  /* extract outputs from return values */
  size  = Degree(graph, 0);  /* = M + N */
  /* vertices are ordered as 1..N then 1..M in slots 1..N+M,
   * so pi2 appears first, in vertices 1..N */
  for (i = 1; i <= N; i++) {
    pi2[i-1] = Mate[i] ? Mate[i] - N : 0; /* rhs maps N+1..N+M -> 1..M */
#ifdef MEX_DEBUG
    printf("%d -> \t%d\n", i, (int) pi2[i-1]);
#endif
  }
  /* pi1 appears in vertices N+1..N+M */
  for (i = N+1; i <= N+M; i++) {
    pi1[i-(N+1)] = Mate[i] ? Mate[i] : 0; /* identity map 1..N -> 1..N */
#ifdef MEX_DEBUG
    printf("%d -> \t%d\n", i, (int) pi1[i-(N+1)]);
#endif
  }
  /* find the cost of that assignment */
  score = 0;
  for (i = 1; i <= size; i++) {
    p = FirstEdge(graph, i);
    for (j = 1; j <= Degree(graph, i); j++) {
      if (EndPoint(p) == Mate[i]) {
        score += ELabel(p);
#ifdef MEX_DEBUG
    printf("%d -> \t%d\t%d\n", i, Mate[i], ELabel(p));
#endif
      }
      p = NextEdge(p);
    }
  }
  /* return score */
  *mxGetPr(plhs[ARG_s]) = (double) score * 0.5;
}


static
Graph 
translateMatrix2Graph(int *size, const mxArray *matrix, int *total_weight)
{
  int edges    = 0;
  int degree   = 0;
  int vlabel   = 0;
  int elabel   = 0;
  int adj_node = 0;
  int i        = 0;
  int j        = 0;
  int classes1 = 0;
  int classes2 = 0;
  int weight   = 0;
  char   c    = 0;
  double *ptr = (double *) NULL;
  Graph graph = (Graph) NULL;


#ifdef MEX_DEBUG
  mexPrintf("Debug: Entering translateMatrix2Graph().\n");
#endif

  classes1 = mxGetN(matrix);
  classes2 = mxGetM(matrix);
  ptr      = mxGetPr(matrix);

  *size = classes1 + classes2;  /* number of nodes in graph */
  edges = classes1 * classes2;
       
  graph = NewGraph( *size );

  for (i = 1; i <= *size; ++i) {
    /*
    // NOTE: Already done in NewGraph().
    // NOTE: Commented out by Ben Bornstein.
    // NOTE: NLabel(graph, i) = i;
    */

    Xcoord(graph, i) = 0;  
    Ycoord(graph, i) = 0;
  }

  *total_weight = 0;

#ifdef MEX_DEBUG
  mexPrintf("Debug: Will call AddEdge() %d times.\n", edges);
#endif

  for (i = 1; i <= classes1; ++i) {
    for (j = 1; j <= classes2; ++j) {
      weight = (int) *ptr++;
      *total_weight += weight;
      AddEdge(graph, i, j + classes1, weight);
    }
  }

#ifdef MEX_DEBUG
  mexPrintf("Debug: Done calling AddEdge() %d times.\n", edges);
  mexPrintf("Debug: Leaving  translateMatrix2Graph().\n");
#endif

  return graph;
}

---