surf的matlab实现

/* * * Petter Strandmark 2008 * * Uses some code from Andrea Vedaldi's SIFT library * */ #include"mex.h" #include "mexutils.c" #include<stdlib.h> #include<string.h> #include<math.h> typedef struct { int k1 ; int k2 ; double score ; } Pair ; void mexFunction(int nout, mxArray *out[], int nin, const mxArray *in[]) { int K1, K2, ND; double thresh = 0.7 ; int THRESH = 0; enum {MATCHES=0,D} ; mxArray *mxdescr1,*mxdescr2,*mxsign1,*mxsign2; double *descr1,*descr2,*sign1,*sign2; int k; Pair *pairs_begin, *pairs_end, *pairs_iterator; int k1, k2; double* d1; double* M_pt ; double* D_pt = NULL ; if (nin < 2 || nin > 5) { mexErrMsgTxt("2,3,4 or 5 input arguments required"); } for (k=0;k<nin;++k) { if(!mxIsDouble(in[k])) { mexErrMsgTxt("Input must be double"); } } if (nin == 2 || nin == 3) { mxdescr1 = in[0]; mxdescr2 = in[1]; mxsign1 = 0; mxsign2 = 0; } else { mxdescr1 = in[0]; mxdescr2 = in[2]; mxsign1 = in[1]; mxsign2 = in[3]; } if (mxGetNumberOfDimensions(mxdescr1) > 2 || mxGetNumberOfDimensions(mxdescr2) > 2) { mexErrMsgTxt("Descriptors must be two dimensional numeric arrays") ; } K1 = mxGetN(mxdescr1); K2 = mxGetN(mxdescr2); ND = mxGetM(mxdescr1); if(mxGetM(mxdescr1) != ND) { mexErrMsgTxt("The two descriptors must have the same number of rows") ; } if (mxsign1 && (mxGetNumberOfElements(mxsign1) != K1 || mxGetNumberOfElements(mxsign2) != K2)) { mexErrMsgTxt("The sign vector must have the same length as the descriptor"); } descr1 = mxGetPr(mxdescr1); descr2 = mxGetPr(mxdescr2); if (mxsign1) { sign1 = mxGetPr(mxsign1); sign2 = mxGetPr(mxsign2); } if(nin == 3) { THRESH = 2; } else if (nin == 5) { THRESH = 4; } if (THRESH) { if(!uIsRealScalar(in[THRESH])) { mexErrMsgTxt("THRESH should be a real scalar") ; } thresh = *mxGetPr(in[THRESH]) ; } /* ------------------------------------------------------------------ ** Do the job ** --------------------------------------------------------------- */ pairs_begin = (Pair*) mxMalloc(sizeof(Pair) * (K1+K2)) ; pairs_iterator = pairs_begin ; d1 = descr1; for(k1 = 0 ; k1 < K1 ; ++k1, d1 += ND ) { double best = mxGetInf(); double second_best = mxGetInf(); int bestk = -1 ; /* For each point P2[k2] in the second image... */ double* d2 = descr2; for(k2 = 0 ; k2 < K2 ; ++k2, d2 += ND) { if (!mxsign1 || (sign1[k1] == sign2[k2])) { int bin ; double acc = 0 ; for(bin = 0 ; bin < ND ; ++bin) { double delta = d1[bin] - d2[bin]; acc += delta*delta; } /* Filter the best and second best matching point. */ if(acc < best) { second_best = best ; best = acc ; bestk = k2 ; } else if(acc < second_best) { second_best = acc ; } } } /* Lowe's method: accept the match only if unique. */ if( (float) best <= thresh * (float) second_best && bestk != -1) { pairs_iterator->k1 = k1 ; pairs_iterator->k2 = bestk ; pairs_iterator->score = second_best / best ; pairs_iterator++ ; } } /* --------------------------------------------------------------- * Finalize * ------------------------------------------------------------ */ pairs_end = pairs_iterator ; D_pt = NULL ; out[MATCHES] = mxCreateDoubleMatrix (2, pairs_end-pairs_begin, mxREAL) ; M_pt = mxGetPr(out[MATCHES]) ; if(nout > 1) { out[D] = mxCreateDoubleMatrix(1, pairs_end-pairs_begin, mxREAL) ; D_pt = mxGetPr(out[D]) ; } for(pairs_iterator = pairs_begin ; pairs_iterator < pairs_end ; ++pairs_iterator) { *M_pt++ = pairs_iterator->k1 + 1 ; *M_pt++ = pairs_iterator->k2 + 1 ; if(nout > 1) { *D_pt++ = pairs_iterator->score ; } } mxFree(pairs_begin) ; }