彭汉川版计算向量概率分布函数---estpa解读

彭汉川版mRMR特征选择算法网址：https://www.mathworks.com/matlabcentral/fileexchange/14608-mrmr-feature-selection-using-mutual-information-computation

estpa函数是为了计算单个向量a的概率分布，其大概思想是：

当a=[1 2 1 2 1]时，其概率分布为 p(x=1)=3/5=0.6, p(x=2)=2/5=0.4，所得P=[0.6 0.4]

当a=[1 5 1 2 1]时，其概率分布为p(x=1)=3/5=0.6, p(x=2)=1/5=0.2, p(x=3)=0/5=0, p(x=4)=0/5=0, p(x=5)=1/5=0.2, 所得P=[0.6 0.2 0 0 0.2]

//=========================================================
//
//This is a prog in the MutualInfo 0.9 package written by
// Hanchuan Peng.
//
//Disclaimer: The author of program is Hanchuan Peng
//      at <penghanchuan@yahoo.com> and <phc@cbmv.jhu.edu>.
//
//The CopyRight is reserved by the author.
//
//Last modification: April/19/2002
//
//========================================================
//
// estpa.cpp
// Calculate the histgram/probability of one vector/image
// By Hanchuan Peng
// Dec/2000
// April/8/2002

#include "miinclude.h"

//return the number of states
template <class T> void copyvecdata(T * srcdata, long len, int * desdata, int& nstate, int &minn, int& maxx);

template <class T> void copyvecdata(T * srcdata, long len, int * desdata, int& nstate, int &minn, int& maxx)
{

  //输入输出变量定义：
  //srcdata:输入向量
  //len：输入向量的长度
  //desdata:输出向量地址
  //nstate:最大值-最小值+1
  //例如输入矩阵为a=[1 2 1 2 1]
  if(!srcdata || !desdata)
  {
    printf("NULL points in copyvecdata()!\n");
    return;
  } 

  long i;

  //note: originally I added 0.5 before rounding, however seems the negative numbers and 
  //      positive numbers are all rounded towarded 0; hence int(-1+0.5)=0 and int(1+0.5)=1;
  //      This is unwanted because I need the above to be -1 and 1.
  // for this reason I just round with 0.5 adjustment for positive and negative differently

  //copy data
  //int minn,maxx;
  if (srcdata[0]>0)
    maxx = minn = int(srcdata[0]+0.5);
  else
    maxx = minn = int(srcdata[0]-0.5);

  int tmp;
  double tmp1;
  for (i=0;i<len;i++)
  {
    tmp1 = double(srcdata[i]);//i=0时，tmp1=1，tmp=2；i=1时，tmp1=2，tmp=3；以此类推
    tmp = (tmp1>0)?(int)(tmp1+0.5):(int)(tmp1-0.5);//round to integers将原始值+/-0.5再取整
    minn = (minn<tmp)?minn:tmp;
    maxx = (maxx>tmp)?maxx:tmp;
    desdata[i] = tmp;
  }//得到minn=2，maxx=3，desdata[]=[2 3 2 3 2]

  //make the vector data begin from 0 (i.e. 1st state)
  for (i=0;i<len;i++)
  {
    desdata[i] -= minn;
  }//desdata[]=[0 1 0 1 0]

  //return the #state
  nstate = (maxx-minn+1);//nstate=3-2+1=2

  return;
}


void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])//matlab调用主函数
{
  // check the arguments检查输入参数是否正确
  //其中nlhs为输出参数个数
  //plhs[]为输出参数数组
  //nrhs为输入参数个数
  //prhs[]为输入参数数组

  if(nrhs != 1 && nrhs !=2 && nrhs!=3)//输入参数不是1/2/3时报错
    mexErrMsgTxt("Usage [marginprob,statelist,cumsumlist] = progname(vector1, maxstatenum, b_returnprob). \n(Both vectors can be images). Max range handled: INT type of the OS");
  if(nlhs > 3)//输出参数>3时报错
    mexErrMsgTxt("Too many output argument <marginalprob_list, statelist, cumsumlist>.");

  if (!mxIsInt8(prhs[0]) && !mxIsUint8(prhs[0]) && !mxIsDouble(prhs[0]) ) //输入参数数据类型不对
    mexErrMsgTxt("The first input argument must be types of INT8 or UINT8 or DOUBLE.");

  //get and check size information

  long i,j;
  //输入数据处理
  void *img1 = (void *)mxGetData(prhs[0]);//获取非数字矩阵的数据值
  long len1 = mxGetNumberOfElements(prhs[0]);//获取数据的数量
  mxClassID type1 = mxGetClassID(prhs[0]);//获取数据的类型

  if (!img1 || !len1)
    mexErrMsgTxt("The input vector is invalid.");

  int b_findstatenum = 1;
  int nstate1 = 0;
  if (nrhs>=2)
  {
    b_findstatenum = 0;
    long MaxGrayLevel = (long) mxGetScalar(prhs[1]);
    nstate1 = MaxGrayLevel;
    if (MaxGrayLevel<=1)
    {
      printf("The argument #state is invalid. This program will decide #state itself.\n");
      b_findstatenum = 1;
    }
  }

  int b_returnprob = 1;
  if (nrhs>=3)
  {
    b_returnprob = (mxGetScalar(prhs[2])!=0);
  }

  //copy data into new INT type array (hence quantization) and then reange them begin from 0 (i.e. state1)

  int * vec1 = new int[len1];
  int nrealstate1=0, minn,maxx;
  switch(type1)
  {
    case mxINT8_CLASS: copyvecdata((char *)img1,len1,vec1,nrealstate1,minn,maxx); break;
    //调用前面的数字取整函数
    case mxUINT8_CLASS: copyvecdata((unsigned char *)img1,len1,vec1,nrealstate1,minn,maxx); break;
    case mxDOUBLE_CLASS: copyvecdata((double *)img1,len1,vec1,nrealstate1,minn,maxx); break;
  }//沿用前面的例子vec1=[0 1 0 1 0],nrealstate1=2

  //update the #state when necessary
  if (nstate1<nrealstate1)
  {
    nstate1 = nrealstate1;
    //printf("First vector #state = %i\n",nrealstate1);
  }

  //generate the marginal-distribution list
  //生成编辑分布

  plhs[0] = mxCreateDoubleMatrix(nstate1,1,mxREAL);
  double *ha = (double *) mxGetPr(plhs[0]);

  for (i=0; i<nstate1;i++)
  {
    ha[i] = 0;
  }//所有值清零
  
  for (i=0;i<len1;i++)
  {
    ha[vec1[i]] += 1;
    //ha[vec1[0]]=ha[0]=1;ha[vec1[1]]=ha[1]=1;循环之后ha[0]=3,ha[1]=2
  }
  
  //return the probabilities, otherwise return count numbers
  if(b_returnprob)
  {
    for (i=0; i<nstate1;i++)
    {
      ha[i] /= len1;
      //ha[0]=3/5=0.6;ha[1]=2/5=0.4
    }
  }

  //return more information
  if (nlhs>=2)  //the second return value is the state list
  {
    plhs[1] = mxCreateDoubleMatrix(nstate1,1,mxREAL);
    double * pstate = (double *)mxGetPr(plhs[1]);
    for (i=0;i<nstate1;i++)
      pstate[i] = minn + i;
  }

  if (nlhs>=3) //the third return value is the cumsum prob list
  {
    plhs[2] = mxCreateDoubleMatrix(nstate1,1,mxREAL);
    double * pcumsum = (double *)mxGetPr(plhs[2]);
    pcumsum[0] = ha[0];
    for (i=1;i<nstate1;i++)
      pcumsum[i] = pcumsum[i-1] + ha[i];
  }

  //free memory
  if (vec1) delete []vec1;


  return;
}