dtw1mex.cpp

来自「一个关于数据聚类和模式识别的程序,在生物化学,化学中因该都可以用到.希望对大家有」· C++ 代码 · 共 101 行

// How to compile:
// MATLAB 7.1: mex dtw1mex.cpp d:/users/jang/c/lib/dtw.cpp -output dtw1mex.dll
// Others: mex dtw1mex.cpp d:/users/jang/c/lib/dtw.cpp

#include <string.h>
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <limits.h>		// For using mxGetInf()
#include "mex.h"
#include "dtw.h"

/* Input Arguments */
#define VEC1		prhs[0]		/*      test input, idx = i */
#define VEC2		prhs[1]		/* reference input, idx = j */
#define BEGINCORNER	prhs[2]
#define ENDCORNER	prhs[3]
#define LOWERBOUND	prhs[4]

/* Output Arguments */
#define	MINDISTANCE	plhs[0]
#define DTWPATH 	plhs[1]
#define DTWTABLE 	plhs[2]

void mexFunction(
	int nlhs,	mxArray *plhs[],
	int nrhs, const mxArray *prhs[])
{
	int i, j, size1, size2, beginCorner, endCorner, featureDim, featureDim2, dtwPathLen, needBackTrack, *dtwPath;
	double *vec1, *vec2, **dtwTable, lowerBound;
	POSITION **prevPos;

	// Error checking
	if ((nrhs<2)||(nrhs>5)||(nlhs>3)) {
		char message[200];
		strcpy(message, mexFunctionName());
		strcat(message, " requires 2~5 input arguments and 0~3 output arguments.\n");
		strcat(message, "Usage: [minDistance, dtwPath, dtwTable] = ");
		strcat(message, mexFunctionName());
		strcat(message, "(vec1, vec2, beginCorner, endCorner, lowerBound)");
		mexErrMsgTxt(message);
	}
	featureDim  = mxGetM(VEC1);
	featureDim2 = mxGetM(VEC2);
	if (featureDim!=featureDim2)
		mexErrMsgTxt("Feature dimensions of input vectors mismatch!");
	size1 = mxGetN(VEC1);
	size2 = mxGetN(VEC2);
	if ((size1==0) || (size2==0))
		mexErrMsgTxt("One or two of the input vectors is empty!");

	// Get input arguments (with possible default values)
	vec1 = mxGetPr(VEC1);
	vec2 = mxGetPr(VEC2);
	beginCorner=1;	if (nrhs>=3) beginCorner = (int)(mxGetPr(BEGINCORNER)[0]);
	endCorner=1;	if (nrhs>=4) endCorner = (int)(mxGetPr(ENDCORNER)[0]);
	lowerBound=INF;	if (nrhs>=5) lowerBound=mxGetPr(LOWERBOUND)[0];

	// Create dtwTable & dtwPath
	dtwPathLen=0;	
	dtwTable=(double **)malloc(size1*sizeof(double *));
	for (i=size1-1;i>=0;i--)
		dtwTable[i]=(double *)malloc(size2*sizeof(double));
	needBackTrack=0;
	if (nlhs>1){		// Need to back track to find the optimum path
		needBackTrack=1;
		dtwPath=(int *)malloc(2*(size1+size2)*sizeof(int));	// Max. length of DTW length is (size1+size2);
		prevPos=(POSITION **)malloc(size1*sizeof(POSITION *));
		for (i=size1-1;i>=0;i--)
			prevPos[i]=(POSITION *)malloc(size2*sizeof(POSITION));
	}

	// Create output argument	
	MINDISTANCE = mxCreateDoubleMatrix(1, 1, mxREAL);
	mxGetPr(MINDISTANCE)[0] = dtw1totalDist(vec1, size1, vec2, size2, featureDim, beginCorner, endCorner, dtwTable, needBackTrack, dtwPath, &dtwPathLen, prevPos, lowerBound);
	if (nlhs>1){		// Create output arguments
		DTWPATH = mxCreateDoubleMatrix(2, dtwPathLen, mxREAL);
		for (i=0; i<2*dtwPathLen; i++)
			mxGetPr(DTWPATH)[i] = dtwPath[i]+1;	// Plus one for MATLAB indexing
	}
	if (nlhs>2){
		DTWTABLE = mxCreateDoubleMatrix(size1, size2, mxREAL);
		for (i=0; i<size1; i++)
			for (j=0; j<size2; j++)
				mxGetPr(DTWTABLE)[j*size1+i]=dtwTable[i][j];
		for (i=0; i<size1*size2; i++)
			if (mxGetPr(DTWTABLE)[i]>=INF)
				mxGetPr(DTWTABLE)[i]=mxGetInf();
	}

	// Free memory
	for (i=0; i<size1;i++)
		free(dtwTable[i]);
	free(dtwTable);
	if (nlhs>1){
		for (i=0; i<size1;i++)
			free(prevPos[i]);
		free(prevPos);
		free(dtwPath);
	}
}