xxxxx.cpp

信息隐藏中用于数字隐写的常用算法:LSB替换LSB匹配,包括随机的和排序的,以及对文件和文件夹进行操作,用CxImage类能快速读取各种格式的图象

	biot = dmatrix(1, phmm->N, 1, T);
	for (i = 1; i <= phmm->N; i++) 
		for (t = 1; t <= T; t++) {
			biot[i][t] = log(phmm->B[i][O[t]]);  //O[t]+1 代替O[t]
		}
 
        /* 1. Initialization  */
 
        for (i = 1; i <= phmm->N; i++) {
                delta[1][i] = phmm->pi[i] + biot[i][1];
                psi[1][i] = 0;
        }
 
        /* 2. Recursion */
 
        for (t = 2; t <= T; t++) {
                for (j = 1; j <= phmm->N; j++) {
                        maxval = -VITHUGE;
                        maxvalind = 1;
                        for (i = 1; i <= phmm->N; i++) {
                                val = delta[t-1][i] + (phmm->A[i][j]);
                                if (val > maxval) {
                                        maxval = val;
                                        maxvalind = i;
                                }
                        }
 
                        delta[t][j] = maxval + biot[j][t]; 
                        psi[t][j] = maxvalind;
 
                }
        }
 
        /* 3. Termination */
 
        *pprob = -VITHUGE;
        q[T] = 1;
        for (i = 1; i <= phmm->N; i++) {
                if (delta[T][i] > *pprob) {
                        *pprob = delta[T][i];
                        q[T] = i;
                }
        }
 
 
	/* 4. Path (state sequence) backtracking */

	for (t = T - 1; t >= 1; t--)
		q[t] = psi[t+1][q[t+1]];

}
 

void BaumWelch(HMM *phmm, int T, int *O, double **alpha, double **beta,
	double **gamma, int *pniter, 
	double *plogprobinit, double *plogprobfinal)
{
	int	i, j;
	int	t, l = 0;
	double theita = 0.0;  //估计嵌入率的特别参数
	double p = 0.0;		  //嵌入率

	double	logprobf, logprobb;
	double	numeratorA, denominatorA;
//	double	numeratorB, denominatorB;

	double ***xi, *scale;
	double delta, deltaprev, logprobprev;

	deltaprev = 10e-70;

	xi = AllocXi(T, phmm->N);
	scale = dvector(1, T);
	
//	Forward(phmm, T, O, alpha, &logprobf);
	ForwardWithScale(phmm, T, O, alpha, scale, &logprobf);
	*plogprobinit = logprobf; /* log P(O |intial model) */
//	Backward(phmm, T, O, beta, &logprobb);
	BackwardWithScale(phmm, T, O, beta, scale, &logprobb);
	ComputeGamma(phmm, T, alpha, beta, gamma);
	ComputeXi(phmm, T, O, alpha, beta, xi);
	logprobprev = logprobf;


	do  {	

		theita = 0.0;

		/* reestimate frequency of state i in time t=1 */
		for (i = 1; i <= phmm->N; i++) 
			phmm->pi[i] = .001 + .999*gamma[1][i];

		/* reestimate transition matrix  and symbol prob in
		   each state */
		for (i = 1; i <= phmm->N; i++) { 
			denominatorA = 0.0;
			for (t = 1; t <= T - 1; t++) 
				denominatorA += gamma[t][i];

			for (j = 1; j <= phmm->N; j++) {  
				numeratorA = 0.0;
				for (t = 1; t <= T - 1; t++) 
					numeratorA += xi[t][i][j];
				phmm->A[i][j] = .001 + 
						.999*numeratorA/denominatorA;
			}

			//需要修改

			/*
			denominatorB = denominatorA + gamma[T][i]; 
			for (k = 1; k <= phmm->M; k++) {
				numeratorB = 0.0;
				for (t = 1; t <= T; t++) {
					if (O[t] == k) 
						numeratorB += gamma[t][i];
				}

				phmm->B[i][k] = .001 +
						.999*numeratorB/denominatorB;
			}
			*/
			for (t = 1; t <= T; t++)
			{
//				theita += xi[t][i][i];
				theita += gamma[t][O[t]+1]; //O[t]+1 代替O[t]
			}
		}

		//对LSB替换方法 估计B, 待验证!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
		theita = theita / T;

		p = 2 * fabs(1 - theita);
		if (p > 1)
			p = 1;
		if (p < 0)
			p = 0;
		for (i = 1; i <= phmm->N; i++)   //又错，待改
		{
			memset(phmm->B[i], 0, phmm->N);
			phmm->B[i][i] = theita;
			if (i % 2 == 1)
				phmm->B[i][i+1] = 1 - theita;
			else 
				phmm->B[i][i-1] = 1 - theita;
		}
		phmm->B[phmm->N][phmm->N -1] = 1 - theita;
		phmm->B[phmm->N][phmm->N] = theita;
		// end of 修改

//	    Forward(phmm, T, O, alpha, &logprobf);
		ForwardWithScale(phmm, T, O, alpha, scale, &logprobf);
//		Backward(phmm, T, O, beta, &logprobb);
		BackwardWithScale(phmm, T, O, beta, scale, &logprobb);
		ComputeGamma(phmm, T, alpha, beta, gamma);
		ComputeXi(phmm, T, O, alpha, beta, xi);

		/* compute difference between log probability of 
		   two iterations */
		delta = logprobf - logprobprev; 
		logprobprev = logprobf;
		l++;
		
	}
	while (delta > DELTA); /* if log probability does not 
                                  change much, exit */ 
 
	*pniter = l;
	*plogprobfinal = logprobf; /* log P(O|estimated model) */
	FreeXi(xi, T, phmm->N);
	free_dvector(scale, 1, T);
}

double BaumWelch2(HMM *phmm, int T, int *O, double **alpha, double **beta,
	double **gamma, int *pniter, 
	double *plogprobinit, double *plogprobfinal)
{
	int	i, j;
	int	t, l = 0;
	double theita = 0.0;  //估计嵌入率的特别参数
	double p = 0.0;		  //嵌入率
	int *q = NULL;        //通过Viterbi估计实际序列q
	int **psi = NULL;     //提供给Viterbi
	double **ddelta = NULL; //
	double logtemp = 0.0;       //

	double	logprobf, logprobb;
	double	numeratorA, denominatorA;
//	double	numeratorB, denominatorB;

	double ***xi, *scale;
	double delta, deltaprev, logprobprev;

	deltaprev = 10e-70;

	xi = AllocXi(T, phmm->N);
	scale = dvector(1, T);
	
//	Forward(phmm, T, O, alpha, &logprobf);
	ForwardWithScale2(phmm, T, O, alpha, scale, &logprobf);
	*plogprobinit = logprobf; /* log P(O |intial model) */
//	Backward(phmm, T, O, beta, &logprobb);
	BackwardWithScale2(phmm, T, O, beta, scale, &logprobb);
	ComputeGamma(phmm, T, alpha, beta, gamma);
	ComputeXi(phmm, T, O, alpha, beta, xi);
	logprobprev = logprobf;

	q = ivector(1, T);
	ddelta = dmatrix(1, T, 1, phmm->N);
	psi = imatrix(1, T, 1, phmm->N);

	ViterbiLog(phmm, T, O, ddelta, psi, q, &logtemp);
//	Viterbi(phmm, T, O, ddelta, psi, q, &logtemp);
	do  
	{	

		theita = 0.0;

		/* reestimate frequency of state i in time t=1 */
		for (i = 1; i <= phmm->N; i++) 
		{
			phmm->pi[i] = .001 + .999*gamma[1][i];
		}

		/* reestimate transition matrix  and symbol prob in
		   each state */
		for (i = 1; i <= phmm->N; i++) 
		{ 
			denominatorA = 0.0;
			for (t = 1; t <= T - 1; t++) 
			{
				denominatorA += gamma[t][i];
			}

			for (j = 1; j <= phmm->N; j++) 
			{  
				numeratorA = 0.0;
				for (t = 1; t <= T - 1; t++) 
				{
					numeratorA += xi[t][i][j];
				}
				phmm->A[i][j] = .001 + .999*numeratorA/(denominatorA + deltaprev);
			}
		}


		//for (t = 1; t <= T; t++)
		//		{
		//			theita += gamma[t][O[t]+1];   // O[t]+1 代替O[t]
		//		}
		
		

		for (t = 1; t <= T; t++)
		{
			if (O[t] == q[t])
				theita = theita + 1;
		}
		
		//对LSB替换方法 估计B, 待验证!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
		
		theita = theita / T;

//		if (theita < 0.5)
//			theita = 0.5;
//		if (theita > 1)
//			theita = 1;
		for (i = 1; i <= phmm->N; i++)   
		{
			memset( &(phmm->B[i][1]), 0, phmm->N * sizeof(double));   // &(phmm->B[i][1]) ，非常重要
			phmm->B[i][i] = theita;
			if (i % 2 == 1)
				phmm->B[i][i+1] = 1 - theita;
			else 
				phmm->B[i][i-1] = 1 - theita;
		}

//	    Forward(phmm, T, O, alpha, &logprobf);
		ForwardWithScale2(phmm, T, O, alpha, scale, &logprobf);
//		Backward(phmm, T, O, beta, &logprobb);
		BackwardWithScale2(phmm, T, O, beta, scale, &logprobb);
		ComputeGamma(phmm, T, alpha, beta, gamma);
		ComputeXi(phmm, T, O, alpha, beta, xi);

		ViterbiLog(phmm, T, O, ddelta, psi, q, &logtemp);
//		Viterbi(phmm, T, O, ddelta, psi, q, &logtemp);

		/* compute difference between log probability of 
		   two iterations */
		delta = logprobf - logprobprev; 
		logprobprev = logprobf;
		l++;
		
	}
	while (delta > DELTA); /* if log probability does not 
                                  change much, exit */ 
 
	*pniter = l;
	*plogprobfinal = logprobf; /* log P(O|estimated model) */
	FreeXi(xi, T, phmm->N);
	free_dvector(scale, 1, T);

	free_ivector(q, 1, T);
	free_imatrix(psi, 1, T, 1, phmm->N);
	free_dmatrix(ddelta, 1, T, 1, phmm->N);

	return theita;
}


void ComputeGamma(HMM *phmm, int T, double **alpha, double **beta, double **gamma)
{

	int 	i, j;
	int	t;
	double	denominator;

	for (t = 1; t <= T; t++) {
		denominator = 0.0;
		for (j = 1; j <= phmm->N; j++) {
			gamma[t][j] = alpha[t][j]*beta[t][j];
			denominator += gamma[t][j];
		}

		for (i = 1; i <= phmm->N; i++) 
			gamma[t][i] = gamma[t][i]/denominator;
	}
}

void ComputeXi(HMM* phmm, int T, int *O, double **alpha, double **beta, 
	double ***xi)
{
	int i, j;
	int t;
	double sum;

	for (t = 1; t <= T - 1; t++)
	{
		sum = 0.0;	
		for (i = 1; i <= phmm->N; i++) 
		{
			for (j = 1; j <= phmm->N; j++)
			{
				xi[t][i][j] = alpha[t][i]*beta[t+1][j]*(phmm->A[i][j])*(phmm->B[j][O[t+1]+1]);   // //O[t+1]+1 代替O[t+1]
				sum += xi[t][i][j];
			}
		}

		for (i = 1; i <= phmm->N; i++) 
			for (j = 1; j <= phmm->N; j++) 
				xi[t][i][j]  /= sum;
	}
}

double *** AllocXi(int T, int N)
{
	int t;
	double ***xi;

	xi = (double ***) malloc(T*sizeof(double **));

	xi --;

	for (t = 1; t <= T; t++)
		xi[t] = dmatrix(1, N, 1, N);
	return xi;
}

void FreeXi(double *** xi, int T, int N)
{
	int t;



	for (t = 1; t <= T; t++)
		free_dmatrix(xi[t], 1, N, 1, N);

	xi ++;
	free(xi);

}


float *vector(int nl, int nh)
//int nl,nh;
{
	float *v;

	v=(float *)calloc((unsigned) (nh-nl+1),sizeof(float));
	memset(v, 0, (nh-nl+1) * sizeof(float));   //自加，赋0值
	return v-nl;
}

int *ivector(int nl, int nh)
//int nl,nh;
{
	int *v;

	v=(int *)calloc((unsigned) (nh-nl+1),sizeof(int));
	memset(v, 0, (nh-nl+1) * sizeof(int));   //自加，赋0值
	return v-nl;
}

double *dvector(int nl, int nh)
//int nl,nh;
{
	double *v;

	v=(double *)calloc((unsigned) (nh-nl+1),sizeof(double));
	memset(v, 0, (nh-nl+1) * sizeof(double));   //自加，赋0值
	return v-nl;
}



float **matrix(int nrl, int nrh, int ncl, int nch)
//int nrl,nrh,ncl,nch;
{
	int i;
	float **m;

	m=(float **) calloc((unsigned) (nrh-nrl+1),sizeof(float*));
	m -= nrl;

	for(i=nrl;i<=nrh;i++) {
		m[i]=(float *) calloc((unsigned) (nch-ncl+1),sizeof(float));
		memset(m[i], 0, (nch-ncl+1) * sizeof(float));   //自加，赋0值
		m[i] -= ncl;
	}
	return m;
}

double **dmatrix(int nrl, int nrh, int ncl, int nch)
//int nrl,nrh,ncl,nch;
{
	int i;
	double **m;

	m=(double **) calloc((unsigned) (nrh-nrl+1),sizeof(double*));
	m -= nrl;

	for(i=nrl;i<=nrh;i++) {
		m[i]=(double *) calloc((unsigned) (nch-ncl+1),sizeof(double));
		memset(m[i], 0, (nch-ncl+1) * sizeof(double));   //自加，赋0值
		m[i] -= ncl;
	}
	return m;
}

int **imatrix(int nrl, int nrh, int ncl, int nch)
//int nrl,nrh,ncl,nch;
{
	int i,**m;

	m=(int **)calloc((unsigned) (nrh-nrl+1),sizeof(int*));
	m -= nrl;

	for(i=nrl;i<=nrh;i++) {
		m[i]=(int *)calloc((unsigned) (nch-ncl+1),sizeof(int));
		memset(m[i], 0, (nch-ncl+1) * sizeof(int));   //自加，赋0值
		m[i] -= ncl;
	}
	return m;
}



float **submatrix(float **a, int oldrl, int oldrh, int oldcl, int oldch, int newrl, int newcl)
//float **a;
//int oldrl,oldrh,oldcl,oldch,newrl,newcl;
{
	int i,j;
	float **m;

	m=(float **) calloc((unsigned) (oldrh-oldrl+1),sizeof(float*));
	m -= newrl;

	for(i=oldrl,j=newrl;i<=oldrh;i++,j++) m[j]=a[i]+oldcl-newcl;

	return m;
}



void free_vector(float *v, int nl, int nh)
//float *v;
//int nl,nh;
{
	free((char*) (v+nl));
}

void free_ivector(int *v, int nl, int nh)
//int *v,nl,nh;
{
	free((char*) (v+nl));
}

void free_dvector(double *v, int nl, int nh)
//double *v;
//int nl,nh;
{
	free((char*) (v+nl));
}



void free_matrix(float **m, int nrl, int nrh, int ncl, int nch)
//float **m;
//int nrl,nrh,ncl,nch;
{
	int i;

	for(i=nrh;i>=nrl;i--) free((char*) (m[i]+ncl));
	free((char*) (m+nrl));
}

void free_dmatrix(double **m, int nrl, int nrh, int ncl, int nch)
//double **m;
//int nrl,nrh,ncl,nch;
{
	int i;

	for(i=nrh;i>=nrl;i--)
		free((char*) (m[i]+ncl));
	free((char*) (m+nrl));
}

void free_imatrix(int **m, int nrl, int nrh, int ncl, int nch)
//int **m;
//int nrl,nrh,ncl,nch;
{
	int i;

	for(i=nrh;i>=nrl;i--) free((char*) (m[i]+ncl));
	free((char*) (m+nrl));
}



void free_submatrix(float **b, int nrl, int nrh, int ncl, int nch)
//float **b;
//int nrl,nrh,ncl,nch;
{
	free((char*) (b+nrl));
}



float **convert_matrix(float *a, int nrl, int nrh, int ncl, int nch)
//float *a;
//int nrl,nrh,ncl,nch;
{
	int i,j,nrow,ncol;
	float **m;

	nrow=nrh-nrl+1;
	ncol=nch-ncl+1;
	m = (float **) calloc((unsigned) (nrow),sizeof(float*));
	m -= nrl;
	for(i=0,j=nrl;i<=nrow-1;i++,j++) m[j]=a+ncol*i-ncl;
	return m;
}



void free_convert_matrix(float **b, int nrl, int nrh, int ncl, int nch)
//float **b;
//int nrl,nrh,ncl,nch;
{
	free((char*) (b+nrl));
}


int  hmmgetseed(void) 
{
	return ((int) getpid());
}

/* 
** hmmsetseed() sets the seed of the random number generator to a
** specific value.
*/
void hmmsetseed(int seed) 
{
	srand(seed);
}

/*
**  hmmgetrand() returns a (double) pseudo random number in the
**  interval [0,1).
*/

double hmmgetrand(void)
{
	return (double) rand()/RAND_MAX;
}