gcut.cpp

本人对经典Maxflow算法的修改

#include "mex.h"
#include "graph.h"
#include <stdio.h>

void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
{
	
	if (nrhs != 2) {
		mexErrMsgTxt("Expecting two arguments.");	
	}

	const mxArray *edges = prhs[0];
	// begin by verifying that the arguments are ok, since mexErrMsgTxt()
	// does not free the memory allocated with new
	// (actually this is a problem anyway since the function is stopped
	// if the user presses Ctrl-C, so we'll have to deal with it)
 	if (!mxIsDouble(edges)){
		mexErrMsgTxt("EDGES should be of class double.");
	}

	int edges_m = mxGetM(edges); //edge_m是边的数量

	if (mxGetN(edges) != 4 || edges_m == 0) {
		mexErrMsgTxt("EDGES should be a N x 4 non-empty matrix.");
	} 

	int max_v = 0; // the maximum node number seen thus far

	double *edges_p = mxGetPr(edges); //edge_p是指向edges,edges指向prhs[0],所以都是指针

	int i;
	for (i=0; i<edges_m; i++) {
		int e1 = (int)edges_p[i];
		int e2 = (int)edges_p[i+edges_m];
		if (e1<1 || e2<1) {
			mexErrMsgTxt("Vertex numbers must be positive integers.");
		}
		if (e1 > max_v) max_v = e1;
		if (e2 > max_v) max_v = e2;
	}	

	// check the TWEIGHTS matrix
	const mxArray *tweights = prhs[1];
	double *tweights_p = mxGetPr(tweights);
	if (!mxIsDouble(tweights)) {
		mexErrMsgTxt("TWEIGHTS should be of class double.");
	}
        int tweights_m = mxGetM(tweights);
        if (mxGetN(tweights) != 3 || tweights_m == 0) {
                mexErrMsgTxt("TWEIGHTS should be a N x 3 non-empty matrix.");
        }
	for (i=0; i<tweights_m; i++) {
		int v = (int)tweights_p[i];
		if (v<1) {
			mexErrMsgTxt("Vertex numbers must be positive integers.");
		}
		if (v > max_v) max_v = v;
	}
	//max_v是edges 和 tweights 中出现的最大节点的数

	// ok, done checking, now do it
	//printf("Number of vertices: %d\n", max_v);
	//printf("Number of edges: %d\n", edges_m);
	//printf("Number of tweights: %d\n", tweights_m);
	
	//Graph *graph = new Graph();
	typedef Graph<double,double,double> GraphType;
	GraphType *graph = new GraphType(max_v, edges_m);
	               /*estimated # of nodes*/ /*estimated # of edges*/ 
	/*例子
	g -> add_node(); 
	g -> add_node(); 

	g -> add_tweights( 0,    capacities   1, 5 );
	g -> add_tweights( 1,    capacities   2, 6 );
	g -> add_edge( 0, 1,     capacities   3, 4 );

	int flow = g -> maxflow();
	例子*/
	
	// mapping from numbers to vertex ids
	// node_array[i] is the node_id of node number i+1
	GraphType::node_id *node_array = new GraphType::node_id[max_v];
	
	for (i=0; i<max_v; i++) {
		node_array[i] = graph->add_node();
	}

	// add edges
	for (i=0; i<edges_m; i++) {
		int v1 = (int)edges_p[i];
		int v2 = (int)edges_p[i+edges_m];
		double w1 = edges_p[i+2*edges_m];
		double w2 = edges_p[i+3*edges_m];
		graph->add_edge(node_array[v1-1], node_array[v2-1], w1, w2);
		//printf("Added edges (%d, %d, %g, %g)\n", v1, v2, w1, w2);
		//fflush(stdout);
	}

	//printf("Done adding edges.\n");
	
	// add edges to terminal nodes
	for (i=0; i<tweights_m; i++) {
		// ASSUME no vertex occurs twice.
		int v = (int)tweights_p[i];
		double s = tweights_p[i+tweights_m];
		double t = tweights_p[i+2*tweights_m];
		graph->add_tweights(node_array[v-1], s, t);
		//printf("Set tweights (%d, %g, %g)\n", v, s, t);
	}


	// really do it
	/*GraphType::flowtype*/double flow = graph->maxflow();
	
	// put result in output parameters (if there are any)
	// put maxflow in 1st parameter
	if (nlhs > 0) {
		// put labeling in 2nd parameter
		plhs[0] = mxCreateDoubleMatrix(1, max_v, mxREAL);
		double *p = mxGetPr(plhs[0]);
		for (i=0; i<max_v; i++) {
			p[i] = (graph->what_segment(node_array[i]) == GraphType::SOURCE)?0:1;
		}
	}
	if (nlhs > 1) {
		plhs[1] = mxCreateDoubleScalar(flow);
	}
	
	// clean up
	delete graph;
	delete [] node_array;	
}