kadjacency.c

经典dijkstra算法程序很完整能够正确运行

/*
*
* 
* Returns the Sparse adjacency matrix of the K neighbours of xi.
*
* Inputs:
* -------
*
*   X      - Data (d x N)
*   K      - Number of neighbourg
*
*   
* Outputs:
* -------
*
*   A      - Sparse adjacency matrix (N x N) with KN non-zeros elements


s                 = 1;
e                 = 3;
d                 = 2;
N                 = 1000;
K                 = 5;
X                 = rand(d , N);
A                 = Kadjacency(X , K);
[path , pathcost] = dijkstra(A , s , e);
gplot(A , X');hold on,plot(X(1 , :) , X(2 , :) , 'k+',X(1 , path) , X(2 , path) , 'r' , 'markersize' , 5 , 'linewidth' , 2), hold off



mex -g -output Kadjacency.dll Kadjacency.c

mex -f mexopts_intelamd.bat -output Kadjacency.dll Kadjacency.c


 Author : S閎astien PARIS : [e-mail]
 -------  Date : 11/01/2007

 Reference ""


*/

#include <math.h>
#include <limits.h>
#include "mex.h"

#define MAX_INF 1.7977e+307

void qsindex( double * , int * , int , int  ); 


void mexFunction(int nlhs, mxArray** plhs ,  int nrhs, const mxArray** prhs)
{
    
    double *X;
	
	int K;
	
	double *A;
	
	
	int d , N;
	
	int i , j, l, k , nzmax;
	
	double *dist;
		
	int *jc , *ir; 
			
	int id , jd , max_inx , ind;
	
	double temp, adist, max_dist , val;
	
	
    if (nrhs != 2)
	{
        
		mxErrMsgTxt("Too few arguments");
		
	}
    
	/* Input 1 */
	
	
    X        = mxGetPr(prhs[0]);
	
    d        = mxGetM(prhs[0]);
	
    N        = mxGetN(prhs[0]);
	
	
	/* Input 2 */
	
    K        = (int) mxGetScalar(prhs[1]);
	
    if (K < 0 || K > N)
	{
        mxErrMsgTxt("K must be < N");
	}
	
	
	
    dist      = mxMalloc(N*sizeof(double));

	nzmax     = K*N;
	
	
	
	
    plhs[0]   = mxCreateSparse(N , N , nzmax , mxREAL);
	
    A         = mxGetPr(plhs[0]);  // (KN x 1) vector
	
    ir        = mxGetIr(plhs[0]);  // (KN x 1) index vector
    
	jc        = mxGetJc(plhs[0]);  // (N+1 x 1) index vector
	
	/* Main Call */
	
	
	k         = 0;
	
	for( j = 0 ; j < N ; j++ ) 
	{
		
		jd            = j*d;
		
		jc[j]         = k;
		
		
		for( i = 0 ; i < N ; i++ ) 
		{
			
			id        = i*d;
			
			adist     = 0.0;
			
			for( l = 0 ; l < d ; l++ ) 
			{
				temp     = (X[l + id] - X[l + jd]);
				
				adist   += (temp*temp); 
			}
			
			if (adist != 0.0)
			{
				
				dist[i]   = adist;
				
			}
			
			else
			{
				
				dist[i]   = MAX_INF;
			}
			
		}
		
		
		for( l = 0 ; l < K ; l++) 
		{
			
			max_dist = MAX_INF;
			
			max_inx  = 0;
			
			
			for( i = 0 ; i < N ; i++ ) 
			{
				if( max_dist > dist[i] ) 
				{
					max_inx  = i;
					
					max_dist = dist[i];
				}
			}
			
			
			ir[k]           = max_inx;
			
			A[k]            = max_dist;
			
			dist[ max_inx ] = MAX_INF;
			
			k++;
			
		}
		
	}
	
	jc[N] = k;
	
	mxFree(dist);
	
}

/* ----------------------------------------------------------------------------- */

void qsindex( double *array, int *index , int left, int right ) 
{
	
	double pivot;	// pivot element.
	
	int ipivot , holex;	// hole index.
	
	int i;
	
	holex          = left + ( right - left ) /2;

	pivot          = array[holex];		     // get pivot from middle of array.

	ipivot         = index[holex];
	
	array[holex]   = array[left];              // move "hole" to beginning of

	index[holex]   = index[left];              // move "hole" to beginning of
	
	holex          = left;			             // range we are sorting.
	
	for ( i = left + 1 ; i <= right ; i++ ) 
	{
		if ( array[ i ] <= pivot ) 
		{
			array[ holex ] = array[ i ];

			index[ holex ] = index[ i ];


			array[ i ]     = array[ ++holex ];

			index[ i ]     = index[ holex ];

		}
	}
	
	if ( holex - left > 1 ) 
	{
	
		qsindex( array, index , left, holex - 1 );
	
	}
	if ( right - holex > 1 ) 
	{
	
		qsindex( array, index , holex + 1, right );
	
	}
	
	array[holex] = pivot;

	index[holex] = ipivot;
	
}