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