条件随机场在二维图像中的应用

#include <math.h> #include "mex.h" /* Sum-Product Belief Propagation on General Graphs * * Usage: * [nodeBel niter nodeMsgs] = BP_General_C(edgePot,nodePot,maximize,maxIter, optTol,starEdge_V,starEdge_E); * * Input: * nodePot(n,k) - Potential at node n for state k * edgePot(k1,k2,e) - Potential on edge e for states k1 and k2 * maxmize - Should be 0, since this code only supports sum-product * maxIter - Maximum number of iterations * opTol - Optimality Tolerance * starEdge_V - Vertex vector in star edge representation * starEdge_E - Edge vector in star edge representation * * Output: * nodeBel(n,k) - Belief at node i for state k * niter - Number of iterations * nodeMsgs(e,k) - Message for state k at edge e * */ int DEBUG = 0; /* Function Declarations */ void checkInput(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) { if (nrhs!=7) mexErrMsgTxt("BP_General_C requires SEVEN Inputs"); if (mxIsChar(prhs[0])||mxIsClass(prhs[0], "sparse")||!mxIsComplex(prhs[0])||mxIsChar(prhs[1])||mxIsClass(prhs[1], "sparse")||!mxIsComplex(prhs[1])) mexErrMsgTxt("Inputs must be COMPLEX, full, and nonstring"); if (mxGetNumberOfDimensions(prhs[0])!=3) mexErrMsgTxt("Edge evidence must be a three dimensional array"); if (mxGetNumberOfDimensions(prhs[1])!=2) mexErrMsgTxt("local evidence must be a two dimensional array"); if (!mxIsClass(prhs[5],"int32")||!mxIsClass(prhs[6],"int32")) mexErrMsgTxt("edge numberings must be int32"); } void printVectorInt(int* X,int nElem) { int i; printf("< "); for(i = 0; i < nElem; i++) { printf("%d ",X[i]);} printf(">\n"); } void printVectorDouble(double* X,int nElem) { int i; printf("< "); for(i = 0; i < nElem; i++) { printf("%f ",X[i]);} printf(">\n"); } void printVectorDoubleComplex(double* X,double* XI,int nElem) { int i; printf("< "); for(i = 0; i < nElem; i++) { printf("%f+i%f ",X[i],XI[i]);} printf(">\n"); } void printMatrixDouble(double* X,int nRows, int nCols) { int i,j; for(i = 0; i < nRows; i++) { printf("< "); for(j = 0; j < nCols; j++) { printf("%f ",X[i+nRows*j]);} printf(">\n");} } void printMatrixDoubleComplex(double* X,double* XI,int nRows, int nCols) { int i,j; for(i = 0; i < nRows; i++) { printf("< "); for(j = 0; j < nCols; j++) { printf("%f+%fi ",X[i+nRows*j],XI[i+nRows*j]);} printf(">\n");} } void print2DArrayDouble(double** X,int nRows,int nCols) { int i,j; for(i=0;i < nRows;i++) { printf("< "); for(j=0;j<nCols;j++) { printf("%f ",X[i][j]); } printf(">\n"); } } void print2DArrayDoubleComplex(double** X,double** XI,int nRows,int nCols) { int i,j; for(i=0;i < nRows;i++) { printf("< "); for(j=0;j<nCols;j++) { printf("%f+%fi ",X[i][j],XI[i][j]); } printf(">\n"); } } void print3DMatrixDouble(double*X,int nRows,int nCols,int height) { int i,j,k; for(k = 0; k < height; k++) { printf("(%d,:,:) = \n",k); for(i = 0; i < nRows; i++) { printf("< "); for(j = 0; j < nCols; j++) { printf("%f ",X[i+nRows*(j+nCols*k)]); } printf(">\n"); } } } int edgeNum(int i,int j,int* starV,int* starE) { int E; if(DEBUG)printf("Looking for edge between %d and %d\n",i,j); for(E=starV[i];E<starV[i+1];E++) { if(DEBUG)printf("Checking %d: Value = %d vs. %d\n",E,starE[E],j); if(j==starE[E]) return E; } mexErrMsgTxt("Could not find Edge Number!"); return 0; } double absoluteDif(double n1,double n2) { if (n1 > n2) return n1-n2; else return n2-n1; } /* return the real part of the multiplication (x+yi)*(u+vi) */ double multiplyR(double x, double y, double u, double v) { return x*u - y*v; } /* return the complex part of the multiplication (x+yi)*(u+vi) */ double multiplyI(double x, double y, double u, double v) { return x*v + y*u; } /* return the real part of the division (x+yi)/(u+vi) */ double divideR(double x, double y, double u, double v) { return (x*u + y*v)/(u*u + v*v); } /* return the complex part of the division (x+yi)/(u+vi) */ double divideI(double x, double y, double u, double v) { /*printf("[x,y,u,v] = [%.3f, %.3f, %.3f, %.3f], result = %.3f\n",x,y,u,v,(-x*v+y*u)/(u*u+v*v));*/ return (-x*v+y*u)/(u*u + v*v); } /* return the absolute value (squared) of (x+yi) */ double absSquared(double x, double y) { return x*x+y*y; } void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) { /* Declare Variables */ int i,j,k,ii,jj,iii,nbrsInd,nbrsInd2,E,iter,n,n2,E2,k_i,k_j, nNodes,nEdges,nStates,maximize,maxIter,converged,nNbrs,nNbrs2, lhs1_dims[2],lhs2_dims[2],lhs3_dims[2],lhs4_dims[3], *starV,*starE, *nbrs,*nbrs2; /* Real Coefficent Variables */ double optTol,sum,realPart,temp2,temp3,temp4, *nodePot,*edgePot, *nodeBel,*nIter,*msgs,*edgeBel, *pot_ij,*temp,*newm, **prod_of_msgs,**old_bel,**new_bel,**old_msg,**new_msg; /* Complex Coefficient Variables */ double sumI,temp2I,temp3I,temp4I, *nodePotI,*edgePotI, *nodeBelI,*msgsI,*edgeBelI, *pot_ijI,*tempI,*newmI, **prod_of_msgsI,**old_belI,**new_belI,**old_msgI,**new_msgI; /* Check Input */ checkInput(nlhs,plhs,nrhs,prhs); /* Get Input Pointers */ nodePot = mxGetPr(prhs[1]); edgePot = mxGetPr(prhs[0]); maximize = mxGetPr(prhs[2])[0]; maxIter = mxGetPr(prhs[3])[0]; optTol = mxGetPr(prhs[4])[0]; starV = mxGetPr(prhs[5]); starE = mxGetPr(prhs[6]); nNodes = mxGetDimensions(prhs[1])[0]; nEdges = mxGetDimensions(prhs[0])[2]; nStates = mxGetDimensions(prhs[1])[1]; nodePotI = mxGetPi(prhs[1]); edgePotI = mxGetPi(prhs[0]); /* Avoid indexing confusion by decrementing the Star arrays */ for(i = 0; i <= nNodes;i++) starV[i]=starV[i]-1; for(j = 0; j < nEdges;j++) starE[j]=starE[j]-1; /* Set-up Output Arrays */ lhs1_dims[0] = nNodes; lhs1_dims[1] = nStates; lhs2_dims[0] = 1; lhs2_dims[1] = 1; lhs3_dims[0] = nEdges; lhs3_dims[1] = nStates; lhs4_dims[0] = nStates; lhs4_dims[1] = nStates; lhs4_dims[2] = nEdges; plhs[0] = mxCreateNumericArray(2,lhs1_dims,mxDOUBLE_CLASS,mxCOMPLEX); plhs[1] = mxCreateNumericArray(2,lhs2_dims,mxDOUBLE_CLASS,mxREAL); plhs[2] = mxCreateNumericArray(2,lhs3_dims,mxDOUBLE_CLASS,mxCOMPLEX); plhs[3] = mxCreateNumericArray(3,lhs4_dims,mxDOUBLE_CLASS,mxCOMPLEX); nodeBel = mxGetPr(plhs[0]); nIter = mxGetPr(plhs[1]); msgs = mxGetPr(plhs[2]); nodeBelI = mxGetPi(plhs[0]); msgsI = mxGetPi(plhs[2]); edgeBel = mxGetPr(plhs[3]); edgeBelI = mxGetPi(plhs[3]); /* Allocate Memory for Auxiliary Arrays */ prod_of_msgs = mxCalloc(nNodes,sizeof(double*)); old_bel = mxCalloc(nNodes,sizeof(double*)); new_bel = mxCalloc(nNodes,sizeof(double*)); old_msg = mxCalloc(nEdges,sizeof(double*)); new_msg = mxCalloc(nEdges,sizeof(double*)); newm = mxCalloc(nStates,sizeof(double)); pot_ij = mxCalloc(nStates*nStates,sizeof(double)); temp = mxCalloc(nStates,sizeof(double)); prod_of_msgsI = mxCalloc(nNodes,sizeof(double*)); old_belI = mxCalloc(nNodes,sizeof(double*)); new_belI = mxCalloc(nNodes