wolffd@0: /* C mex version of parzen.m
wolffd@0: [B,B2] = parzen(feat, mu,  Sigma, Nproto); 
wolffd@0: */
wolffd@0: #include "mex.h"
wolffd@0: #include <stdio.h>
wolffd@0: #include <math.h>
wolffd@0: 
wolffd@0: #define PI 3.141592654
wolffd@0: 
wolffd@0: void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]){
wolffd@0:   int    D, M, Q, T, d, m, q, t;
wolffd@0:   double *data, *mu, *SigmaPtr, *N, Sigma;
wolffd@0:   double *B, *dist, *B2, tmp;
wolffd@0:   const int* dim_mu;
wolffd@0:   double const1, const2, sum_m, sum_d, diff;
wolffd@0:   int Dt, DMq, Dm, MQt, Mq;
wolffd@0:   int dims_B2[3];
wolffd@0: 
wolffd@0:   int ndim_mu, i, save_B2;
wolffd@0:   
wolffd@0:   data = mxGetPr(prhs[0]);
wolffd@0:   mu = mxGetPr(prhs[1]);
wolffd@0:   SigmaPtr = mxGetPr(prhs[2]);
wolffd@0:   Sigma = *SigmaPtr;
wolffd@0:   N = mxGetPr(prhs[3]);
wolffd@0:   
wolffd@0:   D = mxGetM(prhs[0]);
wolffd@0:   T = mxGetN(prhs[0]);
wolffd@0: 
wolffd@0:   ndim_mu = mxGetNumberOfDimensions(prhs[1]);
wolffd@0:   dim_mu = mxGetDimensions(prhs[1]);
wolffd@0:   D = dim_mu[0];
wolffd@0:   M = dim_mu[1];
wolffd@0:   /* printf("parzenC: nlhs=%d, D=%d, M=%d, T=%d\n", nlhs, D, M, T);  */
wolffd@0: 
wolffd@0:   /* If mu is mu(d,m,o,p), then [d M Q] = size(mu) in matlab sets Q=o*p,
wolffd@0:      i.e.. the size of all conditioning variabeles */
wolffd@0:   Q = 1;
wolffd@0:   for (i = 2; i < ndim_mu; i++) {
wolffd@0:     /* printf("dim_mu[%d]=%d\n", i, dim_mu[i]); */
wolffd@0:     Q = Q*dim_mu[i];
wolffd@0:   }
wolffd@0: 
wolffd@0:   /* M = max(N) */
wolffd@0:   M = -1000000;
wolffd@0:   for (i=0; i < Q; i++) {
wolffd@0:     /* printf("N[%d]=%d\n", i, (int) N[i]); */
wolffd@0:     if (N[i] > M) {
wolffd@0:       M = (int) N[i];
wolffd@0:     }
wolffd@0:   }
wolffd@0: 
wolffd@0:   /*  printf("parzenC: nlhs=%d, D=%d, Q=%d, M=%d, T=%d\n", nlhs, D, Q, M, T); */
wolffd@0: 
wolffd@0:   plhs[0] = mxCreateDoubleMatrix(Q,T, mxREAL);
wolffd@0:   B = mxGetPr(plhs[0]);
wolffd@0: 
wolffd@0:   if (nlhs >= 2)
wolffd@0:     save_B2 = 1;
wolffd@0:   else
wolffd@0:     save_B2 = 0;
wolffd@0:   
wolffd@0:   if (save_B2) {
wolffd@0:     /* printf("parzenC saving B2\n");  */
wolffd@0:     /*plhs[1] = mxCreateDoubleMatrix(M*Q*T,1, mxREAL);*/
wolffd@0:     dims_B2[0] = M;
wolffd@0:     dims_B2[1] = Q;
wolffd@0:     dims_B2[2] = T;
wolffd@0:     plhs[1] = mxCreateNumericArray(3, dims_B2, mxDOUBLE_CLASS, mxREAL);
wolffd@0:     B2 = mxGetPr(plhs[1]);
wolffd@0:   } else {
wolffd@0:     /* printf("parzenC not saving B2\n"); */
wolffd@0:   }
wolffd@0:   /*
wolffd@0:   plhs[2] = mxCreateDoubleMatrix(M*Q*T,1, mxREAL);
wolffd@0:   dist = mxGetPr(plhs[2]);
wolffd@0:   */
wolffd@0:   const1 = pow(2*PI*Sigma, -D/2.0);
wolffd@0:   const2 = -(1/(2*Sigma));
wolffd@0:  
wolffd@0:   for (t=0; t < T; t++) {
wolffd@0:     /* printf("t=%d!\n",t); */
wolffd@0:     Dt  = D*t;
wolffd@0:     MQt = M*Q*t;
wolffd@0:     for (q=0; q < Q; q++) {
wolffd@0:       sum_m = 0;
wolffd@0:       DMq = D*M*q;
wolffd@0:       Mq = M*q;
wolffd@0: 
wolffd@0:       for (m=0; m < (int)N[q]; m++) {
wolffd@0: 	sum_d = 0;
wolffd@0: 	Dm = D*m;
wolffd@0: 	for (d=0; d < D; d++) {
wolffd@0: 	  /* diff = data(d,t) - mu(d,m,q) */
wolffd@0: 	  /*diff = data[d + D*t] - mu[d + D*m + D*M*q]; */
wolffd@0: 	  diff = data[d + Dt] - mu[d + Dm + DMq];
wolffd@0: 	  sum_d = sum_d + diff*diff;
wolffd@0: 	}
wolffd@0: 	/* dist[m,q,t] = dist[m + M*q + M*Q*t] = dist[m + Mq + MQt] = sum_d */
wolffd@0: 	tmp = const1 * exp(const2*sum_d);
wolffd@0: 	sum_m = sum_m + tmp;
wolffd@0: 	if (save_B2)
wolffd@0: 	  B2[m + Mq + MQt] = tmp;
wolffd@0:       }
wolffd@0: 
wolffd@0:       if (N[q]>0) {
wolffd@0: 	B[q + Q*t] = (1.0/N[q]) *  sum_m;
wolffd@0:       } else {
wolffd@0: 	B[q + Q*t] = 0.0;
wolffd@0:       }
wolffd@0:     }
wolffd@0:   }
wolffd@0: }
wolffd@0: 
wolffd@0: 
wolffd@0: