match.c

A salient-boundary extraction software package based on the paper: S. Wang, T. Kubota, J. M. Siskind

#include "match.h"
#define NDEBUG   /* remove this to turn assertions on */
#include <assert.h>

#define GREEDY_DUAL_CHANGE 
#define SWITCH_LEVEL 32
#define MAX_BAD_PORTION 0.25

#define PRINT_LEVEL 0
#define SHRINKS_SIZE 1000
/* infinity */
#define INTEGER_MIN -999999999
#define INTEGER_MAX  999999999
/* for labels */
#define NIX 0
#define PLUS 1
#define MINUS 2
/* for status */
#define UNMATCHED 0
#define HALVES 1
#define MATCHED 2
/* for x */
#define HALF 2

/*
#define OTHEREND(e,n,list)  ((list) + (e)->nod1 + (e)->nod2 - ((n) - (list)))
*/
#define OTHEREND(e,n,list) (((e)->nod1 + (list)) == (n) ?                \
                             (e)->nod2 + (list) : (e)->nod1 + (list))
#define OTHEREND_INT(e,n) ((e)->nod1 + (e)->nod2 - (n))

#define FIND_SURFACE(n) {                                                \
            while((n)->blossom_parent != -1) {                           \
                (n)=G->nodelist+(n)->blossom_parent;                     \
            }                                                            \
        }

#define SHRINK_SORT_CONST 3
#define BADCOUNT_CONST 10000
#define PNODE(n) (G->nodelist + (n))
#define PEDGE(e) (G->edgelist + (e))

#define INODE(n) ((n) - G->nodelist)
#define IEDGE(e) ((e) - G->edgelist)

typedef struct edge {
    int   slack;
    int   label; /* added by S. Wang */
    char  mark;  
    char  x; 
    int   ptrs[2];
    int   nod1;
    int   nod2;
    int   orig_nod1;
    int   orig_nod2;
} edge;

typedef struct node {
    int   edg_list;
    int   matched_edge;

    /* These are for the augmenting path tree */
    int   child;
    int   sibling;
    int   parent;
    int   parent_edge;

    /* These are for the blossom_tree */
    int   blossom_next;
    int   blossom_parent;
    int   blossom_root;
    int   blossom_next_edge;
    int   penultimate;

    int   pi;
    int   mark;      /* unused pseudo nodes are linked by mark */
    int   tree_root;

    char  status;
    char  label;
    char  hit;

    int   dummy;    /* just to pad to multiple of 8 bytes */
} node;

typedef struct nodeptr {
    int next;
    int surf;
    int delta;
    int dad;
    int sum;
    int status;
#ifdef GREEDY_DUAL_CHANGE
    char tree_processed;
    int gnext;      /* used to pad to multiple of 8 bytes and in reordering */
#endif
} nodeptr;

typedef struct shrink_ary_T {
    int node_i;
    int node_j;
    int node_k;
    int edge;
    int size;
    int next;
} shrink_ary_T;

typedef struct shrink_T {
    shrink_ary_T *ary;
    int length;
} shrink_T;

typedef struct expand_T {
    int *node;
    int length;
} expand_T;

typedef struct stats {
    int expand_count;
    int shrink_count;
    int dualchange_count;
    int dualzero_count;
} stats;

typedef struct graph {
    edge *edgelist;
    node *nodelist;
    int nnodes;
    int nedges;
    int max_nedges;
    int unused;
    int unmatched;
    nodeptr *roots;
    int unmatched_count;
}graph;

typedef struct srkstuff {
    expand_T expands;
    shrink_T shrinks;
    int shrinks_max;
    int expands_max;
    int shrinks_[SHRINKS_SIZE];
} srkstuff;

typedef struct stack {
    int*             ary;
    int              count;
} stack;

#ifdef  CC_PROTOTYPE_ANSI

extern void
    adjust_match (graph *G),
    init_tree (graph *G, node *n),
    clear_tree (graph *G, node *n),
    make_cycle_halves (graph *G, node *node_i, node *node_j, node *node_k,
        edge *e),
    shrink_tree (graph *G, node *newnode),
    label_penultimate (graph *G, node *n, int label),
    fix_matching (graph *G, node *oldnode, node *x),
    fix_tree (graph *G, node *oldnode, node *x, node *y),
    flip_cycle (graph *G, node *node_i),
    augment_path (graph *G, node *node_i, int fractional),
    vereinigung (graph *G, int x, int y),
    set_vereinigung (graph *G, node *n),
    unmatch_edge (graph *G, edge *e),
    dual_match_len (graph *G, int fractional, double *val),
    fix_match (graph *G, int blossom), 
    *CCutil_allocrus (unsigned int size),
    CCutil_freerus (void *p);

extern int 
    perfect_match (graph *G, int *elen),
    build_graph (graph *G, int ncount, int ecount, int *elist, int *elen),
    init (graph *G, srkstuff *srk),
    match_main_frac (graph *G, stats *scount, srkstuff *srk),
    match_main (graph *G, stats *scount, srkstuff *srk, int use_all),
    augment_blossom (graph *G, edge *, int fractional, srkstuff *srk),
    lower_edges (graph *G, node *old),
    expand_blossom (graph *G, node *oldnode, stats *scount),
    lift_edges (graph *G, node *newnode),
    add_to_tree (graph *G, edge *e, node *node_i, node *node_j, int fractional),
    checkout_node (graph *G, node *n, int fractional, srkstuff *srk),
    grow_tree_no_shrink (graph *G, node *n, int fractional, srkstuff *srk),
    grow_tree (graph *G, node *n, int fractional, stats *scount, srkstuff *srk,
               int *found_aug),
    apply_dual_change (graph *G, node *n, int delta),
    find_parity_sum (graph *G, int n),
    parity_correction (graph *G, stats *scount),
    find_single_dual_change (graph *G, node *n),
    find_dual_change_forest (graph *G, node *n),
    make_dual_change_forest (graph *G, stats *scount),
    match_main_more_in_forest (graph *G, stats *scount, srkstuff *srk),
    match_main_forest (graph *G, stats *scount, srkstuff *srk),
    match_main_tree (graph *G, stats *scount, srkstuff *srk),
    make_match (graph *G),
    write_match (graph *G, int *elen, int *elen2, int *esolid, int *elabel_penul, char* match_file),
    label_match (graph *G, int *elen, int *elen2, int *elen_bak, int *esolid, 
                 int *cycle1_min, int *cycle2_min, int *cycle1_orig_min, int *elabel),
    CCutil_readint (FILE *f);

extern node
    *shrink_blossom (graph *G, node *node_i, node *node_j, node *node_k,
        edge *e, stats *scount),
    *common_parent (graph *G, node *node_i, node *node_j, int *size),
    *find_below (graph *G, node *n, int blossom);

#else

extern void
    adjust_match (),
    init_tree (),
    clear_tree (),
    make_cycle_halves (),
    shrink_tree (),
    label_penultimate (),
    fix_matching (),
    fix_tree (),
    flip_cycle (),
    augment_path (),
    vereinigung (),
    set_vereinigung (),
    unmatch_edge (),
    dual_match_len (),
    fix_match (),
    *CCutil_allocrus (),
    CCutil_freerus ();

extern int 
    perfect_match (),
    build_graph (),
    init (),
    match_main_frac (),
    match_main (),
    augment_blossom (),
    lower_edges (),
    expand_blossom (),
    lift_edges (),
    add_to_tree (),
    checkout_node (),
    grow_tree_no_shrink (),
    grow_tree (),
    apply_dual_change (),
    find_parity_sum (),
    parity_correction (),
    find_single_dual_change (),
    find_dual_change_forest (),
    make_dual_change_forest (),
    make_match (),
    match_main_more_in_forest (),
    match_main_forest (),
    match_main_tree (),
    write_match (),
    label_match (),
    CCutil_readint ();

extern node
    *shrink_blossom (),
    *common_parent (),
    *find_below ();

#endif

#if PRINT_LEVEL
static graph *PG = (graph *) NULL;
#endif

#ifdef CC_PROTOTYPE_ANSI
int min_ratio_cycle (int ncount, int ecount, int **elist, int **elen, 
    int **elen2, int **esolid, char *mat_filename, int NumIteration)
#else
int min_ratio_cycle (ncount, ecount, elist, elen, elen2, esolid, 
                     mat_filename, NumIteration)
int ncount, ecount;
int **elist, **elen, **elen2, **esolid;
char *mat_filename;
#endif
{
    graph G;
    int **elabel, **elabel_penul, **elen_bak, **elen2_bak, **elen_bak2;
    int hh, i, end1, end2, cycle1_min, cycle2_min, cycle1_orig_min;
    int iter;
    char current_filename[255];
    int node_inv[ncount], merge_w1[ncount], merge_w2[ncount]; 
    //edge *e;

    hh =0;
    for (i = 0; i < ecount; i++)
      hh = hh + ((*esolid)[i])%2; /* 0: dashed; 1: solid; 2: boundary dashed */

    fprintf (stderr, "# of node = %i\n", ncount);
    fprintf (stderr, "# of solid-edges = %i\n", hh);
    fprintf (stderr, "# of dashed-edges= %i\n", ecount-hh);
  
    elabel = CC_SAFE_MALLOC(1, int *); 
    elabel_penul = CC_SAFE_MALLOC(1, int *);
    elen_bak = CC_SAFE_MALLOC(1, int *);
    elen2_bak = CC_SAFE_MALLOC(1, int *);
    elen_bak2 = CC_SAFE_MALLOC(1, int *);

    *elabel = CC_SAFE_MALLOC(ecount, int);
    if (!(*elabel)) {
        fprintf (stderr, "out of memory in getedgelist\n");
        CC_FREE (*elabel, int);
        return 1;
    }

    *elabel_penul = CC_SAFE_MALLOC(ecount, int);
    if (!(*elabel_penul)) {
        fprintf (stderr, "out of memory in getedgelist\n");
        CC_FREE (*elabel_penul, int);
        return 1;
    }

    *elen_bak = CC_SAFE_MALLOC(ecount, int);
    if (!(*elen_bak)) {
        fprintf (stderr, "out of memory in getedgelist\n");
        CC_FREE (*elen_bak, int);
        return 1;
    }

    *elen2_bak = CC_SAFE_MALLOC(ecount, int);
    if (!(*elen2_bak)) {
        fprintf (stderr, "out of memory in getedgelist\n");
        CC_FREE (*elen2_bak, int);
        return 1;
    }

    *elen_bak2 = CC_SAFE_MALLOC(ecount, int);
    if (!(*elen_bak2)) {
       fprintf (stderr, "out of memory in getedgelist\n");
       CC_FREE (*elen_bak2, int);
       return 1;
    }

    for (iter = 0; iter < NumIteration; iter ++){

        sprintf(current_filename, "%s-%d.cycle", mat_filename, iter);

        for (i = 0; i < ncount; i++) 
	     node_inv[i] = 0;

        for (i = 0; i < ecount; i++) {
          (*elabel)[i] = 0;  /* initialize the edge labels*/
          (*elen_bak)[i] = (*elen)[i]; /* make a copy of edge weight */          
          (*elen2_bak)[i] = (*elen2)[i]; /* make a copy of second edge weight */
          (*elen)[i] = (*elen)[i]*2; /* double all the first weights */          
          (*elen2)[i] = (*elen2)[i]*2; /* double all the second weights */

        }

        for (i = 0; i < ecount; i++) {
	   if ((*esolid)[i] == 1) {
	      end1 = (*elist)[i*2];
              end2 = (*elist)[i*2+1];
              merge_w1[end1] = (*elen_bak)[i]; /* set solid-edge weights to be zero */
              merge_w2[end1] = (*elen2_bak)[i];
              merge_w1[end2] = (*elen_bak)[i];
              merge_w2[end2] = (*elen2_bak)[i];
              (*elen)[i] = 0;
              (*elen2)[i] = 0;
	   } 
        }
       
        for (i = 0; i < ecount; i++) {
	  if ((*esolid)[i] != 1) {
	     end1 = (*elist)[i*2]; /* merge solid weights to the dashed ones */
             end2 = (*elist)[i*2+1];
             (*elen)[i] = (*elen)[i] + merge_w1[end1] + merge_w1[end2];
             (*elen2)[i] = (*elen2)[i] + merge_w2[end1] + merge_w2[end2];
	  } 
        } 


    #if PRINT_LEVEL
        PG = &G;
    #endif

        cycle1_min = 600;
        //cycle1_orig_min = 600;
        cycle1_orig_min = 0; //This is for negative cycles
        cycle2_min=1;

        for (i = 0; i < ecount; i++) 
           (*elen_bak2)[i] = (*elen)[i]; /* make a temp copy of current edge weight */          

        /* linearly transforming the (first) edge weights */
        while (cycle1_min!=0) {
  
            for (i = 0; i < ecount; i++) {
	      (*elabel_penul)[i] = (*elabel)[i]; /* record the min-ratio-cycle in the penultimate step */
              (*elen)[i]=(*elen_bak2)[i] * cycle2_min - cycle1_orig_min * (*elen2)[i] ;
	    }

            G.edgelist = (edge *) NULL;
            G.nodelist = (node *) NULL;
            G.unused = -1;
            G.unmatched = -1;
            G.roots = (nodeptr *) NULL;

            if (build_graph (&G, ncount, ecount, *elist, *elen)) {
                fprintf (stderr, "build_graph failed\n");
                goto CLEANUP;
            }

            /* do perfect matching on solid-dashed graph G */
            if (perfect_match (&G, *elen)){
                fprintf (stderr, "perfect matching failed\n");
                goto CLEANUP;
            }

            /* find the (minimum-cycle-ratio) negative weight cycle with perfect matching */
            if (label_match (&G, *elen, *elen2, *elen_bak2, *esolid, 
                       &cycle1_min, &cycle2_min, &cycle1_orig_min, *elabel)){
                fprintf (stderr, "match labelling failed\n");
                goto CLEANUP;
            }
        }

        for (i = 0; i < ecount; i++) {
            (*elen)[i]=(*elen_bak)[i]; /* recover the edge weight */
            (*elen2)[i]=(*elen2_bak)[i];
	}

        hh =0;
        for (i = 0; i < ecount; i++)
            hh = hh + (*elabel_penul)[i];
        fprintf (stderr, "# edges(nodes) in MRC = %i\n", hh);

        /* if the mrc ratio is zero, go ahead and output the result */
        if (current_filename !=  (char *) NULL) {
            if (write_match (&G, *elen, *elen2, *esolid, *elabel_penul, current_filename)) {
                fprintf (stderr, "write_match failed\n");
                goto CLEANUP;
            }
        }

        for (i = 0; i < ecount; i++) {
          if (((*elabel_penul)[i] == 1) && ((*esolid)[i] != 1)) {
	      node_inv[G.edgelist[i].orig_nod1] = 1;
              node_inv[G.edgelist[i].orig_nod2] = 1;
              // following added for negative cycle with mirror edges
              node_inv[G.edgelist[i].orig_nod1-(2*(G.edgelist[i].orig_nod1 % 2)-1)] = 1;
              node_inv[G.edgelist[i].orig_nod2-(2*(G.edgelist[i].orig_nod2 % 2)-1)] = 1;
              //
	  }
	}

        for (i = 0; i < ecount; i++)
	  if ((*esolid)[i] != 1)
	    if ((node_inv[G.edgelist[i].orig_nod1] == 1) || (node_inv[G.edgelist[i].orig_nod2] == 1)) {
               //(*elen)[i]=600;
               //(*elen2)[i]=1;
               (*elen)[i] = 600; // modified for negative cycle
               (*elen2)[i] = 0;
            }
    }

CLEANUP:

    CC_IFFREE (G.nodelist, node);
    CC_IFFREE (G.edgelist, edge);
    CC_IFFREE (G.roots, nodeptr);
    CC_FREE (*elist, int);
    CC_FREE (*elabel, int); 
    CC_FREE (*elabel_penul, int);
    CC_FREE (*elen, int);
    CC_FREE (*elen_bak, int);
    CC_FREE (elabel, int *);
    CC_FREE (elabel_penul, int *);
    CC_FREE (elen_bak, int *);
    CC_FREE (elen2_bak, int *);
    CC_FREE (elen_bak2, int *);

    return 0;
}

#ifdef CC_PROTOTYPE_ANSI
extern int perfect_match (graph *G, int *elen)
#else
extern int perfect_match (G, elen)
graph *G;
int *elen;
#endif
{
    srkstuff srk;
    stats scount;
    double val;
    int finished = 1;
    int use_all_trees=0;

    srk.expands.node = (int *) NULL;
    srk.shrinks.ary = (shrink_ary_T *) NULL;