xxxxx.cpp

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

// xxxxx.cpp: implementation of the xxxxx class.
//
//////////////////////////////////////////////////////////////////////

#include "stdafx.h"
#include "LSBAnalysis.h"
#include "xxxxx.h"
#include <stdio.h> 
#include <math.h>
#include <malloc.h>
#include <stdlib.h>
#include "process.h"

#define DELTA 0.001 
#define VITHUGE  100000000000.0

#ifdef _DEBUG
#undef THIS_FILE
static char THIS_FILE[]=__FILE__;
#define new DEBUG_NEW
#endif


void ReadHMM(FILE *fp, HMM *phmm)
{
	int i, j, k;

	fscanf(fp, "M= %d\n", &(phmm->M)); 

	fscanf(fp, "N= %d\n", &(phmm->N)); 

	fscanf(fp, "A:\n");
	phmm->A = (double **) dmatrix(1, phmm->N, 1, phmm->N);
	for (i = 1; i <= phmm->N; i++) { 
		for (j = 1; j <= phmm->N; j++) {
			fscanf(fp, "%lf", &(phmm->A[i][j])); 
		}
		fscanf(fp,"\n");
	}

	fscanf(fp, "B:\n");
	phmm->B = (double **) dmatrix(1, phmm->N, 1, phmm->M);
	for (j = 1; j <= phmm->N; j++) { 
		for (k = 1; k <= phmm->M; k++) {
			fscanf(fp, "%lf", &(phmm->B[j][k])); 
		}
		fscanf(fp,"\n");
	}

	fscanf(fp, "pi:\n");
	phmm->pi = (double *) dvector(1, phmm->N);
	for (i = 1; i <= phmm->N; i++) 
		fscanf(fp, "%lf", &(phmm->pi[i])); 

}

void FreeHMM(HMM *phmm)
{
	free_dmatrix(phmm->A, 1, phmm->N, 1, phmm->N);
	free_dmatrix(phmm->B, 1, phmm->N, 1, phmm->M);
	free_dvector(phmm->pi, 1, phmm->N);
}

/*
** InitHMM() This function initializes matrices A, B and vector pi with
**	random values. Not doing so can result in the BaumWelch behaving
**	quite weirdly.
*/ 

void InitHMM(HMM *phmm, int N, int M, int seed)
{
/*
	int i, j;
	double sum;

	hmmsetseed(seed);	

    phmm->M = M;
	
	phmm->N = N;
	
	phmm->A = (double **) dmatrix(1, phmm->N, 1, phmm->N);

	for (i = 1; i <= phmm->N; i++) 
	{
		sum = 0.0;
		for (j = 1; j <= phmm->N; j++) 
		{
			phmm->A[i][j] = hmmgetrand(); 
			sum += phmm->A[i][j];
		}
		for (j = 1; j <= phmm->N; j++) 
			phmm->A[i][j] /= sum;
	}

	phmm->B = (double **) dmatrix(1, phmm->N, 1, phmm->M);
	
	for (j = 1; j <= phmm->N; j++) 
	{
		sum = 0.0;	
		for (k = 1; k <= phmm->M; k++)
		{
			phmm->B[j][k] = hmmgetrand();
			sum += phmm->B[j][k];
		}
		for (k = 1; k <= phmm->M; k++) 
			phmm->B[j][k] /= sum;
	}
 
	phmm->pi = (double *) dvector(1, phmm->N);
	sum = 0.0;
	for (i = 1; i <= phmm->N; i++) 
	{
		phmm->pi[i] = hmmgetrand(); 
		sum += phmm->pi[i];
	}
	for (i = 1; i <= phmm->N; i++) 
		phmm->pi[i] /= sum;

*/

	int i, j;
	double sum;

	hmmsetseed(seed);	

    phmm->M = M;
	
	phmm->N = N;
	
	phmm->A = (double **) dmatrix(1, phmm->N, 1, phmm->N);

	for (i = 1; i <= phmm->N; i++)              //随机产生A概率分布
	{
		sum = 0.0;
		for (j = 1; j <= phmm->N; j++) 
		{
			phmm->A[i][j] = hmmgetrand(); 
			sum += phmm->A[i][j];
		}
		for (j = 1; j <= phmm->N; j++) 
			phmm->A[i][j] /= sum;
	}
	/*

	for (i = 1; i <= phmm->N; i++)              //均匀产生A概率分布
	{
		for (j = 1; j <= phmm->N; j++) 
		{
			phmm->A[i][j] = 1.0 / (phmm->N * phmm->N);
		}
	}
	*/


    phmm->B = (double **) dmatrix(1, phmm->N, 1, phmm->M);

	double theita = hmmgetrand();
//	double theita = 0.75;
	for (j = 1; j <= phmm->N; j++)   
	{
		phmm->B[j][j] = theita;
		if (j % 2 == 1)
			phmm->B[j][j+1] = 1 - theita;
		else 
			phmm->B[j][j-1] = 1 - theita;
	}

	phmm->pi = (double *) dvector(1, phmm->N);
	for (i = 1; i <= phmm->N; i++)
	{
		phmm->pi[i] = 1.0 / phmm->N;                  //均匀分布初始概率
	}
	/*
	sum = 0.0;                                        //随机分布初始概率
	for (i = 1; i <= phmm->N; i++) 
	{
		phmm->pi[i] = hmmgetrand(); 
		sum += phmm->pi[i];
	}
	for (i = 1; i <= phmm->N; i++) 
		phmm->pi[i] /= sum;
	*/
}


void InitHMM2(BYTE *pImg, int dwWidth, int dwHeight, HMM *phmm, int N, int M, int seed)
{
	int i, j;
	double sum;
	double eps = 1.0e-30;

	/* initialize random number generator */

	//A 初始化
	hmmsetseed(seed);
	phmm->M = M;
	phmm->N = N;
	phmm->A = (double **) dmatrix(1, phmm->N, 1, phmm->N);

	BYTE temp1, temp2;
	for (i = 1; i <= dwHeight; i++)
	{
		for (j = 1; j < dwWidth; j++)
		{
			temp1 = pImg[(i - 1) * dwWidth + (j - 1)] + 1;
			temp2 = pImg[(i - 1) * dwWidth + j] + 1;
			phmm->A[temp1][temp2] ++;
		}
	}

    for (i = 1; i <= phmm->N; i++) 
	{
		sum = 0.0;
		for (j = 1; j <= phmm->N; j++)
		{
			sum += phmm->A[i][j];
		}
		for (j = 1; j <= phmm->N; j++) 
			 phmm->A[i][j] /= (sum + eps);
	}

	//B 初始化
 
    phmm->B = (double **) dmatrix(1, phmm->N, 1, phmm->M);

	double theita = hmmgetrand();
	for (j = 1; j <= phmm->N; j++)   
	{
		memset( &(phmm->B[j][1]), 0, phmm->N);
		phmm->B[j][j] = theita;
		if (j % 2 == 1)
			phmm->B[j][j+1] = 1 - theita;
		else 
			phmm->B[j][j-1] = 1 - theita;
	}
//	phmm->B[phmm->N][phmm->N -1] = 1 - theita;
//	phmm->B[phmm->N][phmm->N] = theita;
	


    phmm->pi = (double *) dvector(1, phmm->N);
//	sum = 0.0;
		//	for (i = 1; i <= phmm->N; i++)
		//	{
		//		phmm->pi[i] = hmmgetrand(); 
		//		sum += phmm->pi[i];
		//	}
		//	for (i = 1; i <= phmm->N; i++) 
		//		phmm->pi[i] /= sum;
		
//	for (i = 1; i <= phmm->N; i++)
	//	{
	//		phmm->pi[i] = 1.0 / phmm->N;                  //均匀分布初始概率
	//	}
	for (i = 1; i <= phmm->N; i++) 
		phmm->pi[i] = 0;
	if(pImg[0]% 2 == 1)
	{
		phmm->pi[pImg[0]+1] = theita;
		phmm->pi[pImg[0]] = 1 - theita;
	}
	else
	{
		phmm->pi[pImg[0]+1] = theita;
		phmm->pi[pImg[0]+2] = 1 - theita;
	}
	
}


void CopyHMM(HMM *phmm1, HMM *phmm2)
{
        int i, j, k;
 
        phmm2->M = phmm1->M;

 
        phmm2->N = phmm1->N;
 
        phmm2->A = (double **) dmatrix(1, phmm2->N, 1, phmm2->N);
 
        for (i = 1; i <= phmm2->N; i++)
                for (j = 1; j <= phmm2->N; j++)
                        phmm2->A[i][j] = phmm1->A[i][j];
 
        phmm2->B = (double **) dmatrix(1, phmm2->N, 1, phmm2->M);
        for (j = 1; j <= phmm2->N; j++)
                for (k = 1; k <= phmm2->M; k++)
                        phmm2->B[j][k] = phmm1->B[j][k];
 
        phmm2->pi = (double *) dvector(1, phmm2->N);
        for (i = 1; i <= phmm2->N; i++)
                phmm2->pi[i] = phmm1->pi[i]; 
 
}

void PrintHMM(FILE *fp, HMM *phmm)
{
        int i, j, k;

	fprintf(fp, "M= %d\n", phmm->M); 
	fprintf(fp, "N= %d\n", phmm->N); 
 
	fprintf(fp, "A:\n");
        for (i = 1; i <= phmm->N; i++) {
                for (j = 1; j <= phmm->N; j++) {
                        fprintf(fp, "%f ", phmm->A[i][j] );
		}
		fprintf(fp, "\n");
	}
 
	fprintf(fp, "B:\n");
        for (j = 1; j <= phmm->N; j++) {
                for (k = 1; k <= phmm->M; k++){
                        fprintf(fp, "%f ", phmm->B[j][k]);
		}
		fprintf(fp, "\n");
	}
 
	fprintf(fp, "pi:\n");
        for (i = 1; i <= phmm->N; i++) {
		fprintf(fp, "%f ", phmm->pi[i]);
	}
	fprintf(fp, "\n\n");
}


void Forward(HMM *phmm, int T, int *O, double **alpha, double *pprob)
{
        int     i, j;   /* state indices */
        int     t;      /* time index */
 
        double sum;     /* partial sum */
 
        /* 1. Initialization */
 
        for (i = 1; i <= phmm->N; i++)
                alpha[1][i] = phmm->pi[i]* phmm->B[i][O[1]+1];  //O[1]+1 代替O[1]
 
        /* 2. Induction */
 
        for (t = 1; t < T; t++) {
                for (j = 1; j <= phmm->N; j++) {
                        sum = 0.0;
                        for (i = 1; i <= phmm->N; i++)
                                sum += alpha[t][i]* (phmm->A[i][j]);
 
                        alpha[t+1][j] = sum*(phmm->B[j][O[t+1]+1]);  // //O[t+1]+1 代替O[t+1]
                }
        }
 
        /* 3. Termination */
        *pprob = 0.0;
        for (i = 1; i <= phmm->N; i++)
                *pprob += alpha[T][i];
 
}

void ForwardWithScale(HMM *phmm, int T, int *O, double **alpha, double *scale, double *pprob)
/*  pprob is the LOG probability */
{
	int	i, j; 	/* state indices */
	int	t;	/* time index */

	double sum;	/* partial sum */

	/* 1. Initialization */

	scale[1] = 0.0;	
	for (i = 1; i <= phmm->N; i++) {
		alpha[1][i] = phmm->pi[i]* (phmm->B[i][O[1]+1]);  //O[1]+1 代替O[1]
		scale[1] += alpha[1][i];
	}
	for (i = 1; i <= phmm->N; i++) 
		alpha[1][i] /= scale[1]; 
	
	/* 2. Induction */

	for (t = 1; t <= T - 1; t++) {
		scale[t+1] = 0.0;
		for (j = 1; j <= phmm->N; j++) {
			sum = 0.0;
			for (i = 1; i <= phmm->N; i++) 
				sum += alpha[t][i]* (phmm->A[i][j]); 

			alpha[t+1][j] = sum*(phmm->B[j][O[t+1]+1]);   // //O[t+1]+1 代替O[t+1]
			scale[t+1] += alpha[t+1][j];
		}
		for (j = 1; j <= phmm->N; j++) 
			alpha[t+1][j] /= scale[t+1]; 
	}

	/* 3. Termination */
	*pprob = 0.0;

	for (t = 1; t <= T; t++)
		*pprob += log(scale[t]);
	
}


void Backward(HMM *phmm, int T, int *O, double **beta, double *pprob)
{
        int     i, j;   /* state indices */
        int     t;      /* time index */
        double sum;
 
 
        /* 1. Initialization */
 
        for (i = 1; i <= phmm->N; i++)
                beta[T][i] = 1.0;
 
        /* 2. Induction */
 
        for (t = T - 1; t >= 1; t--) {
                for (i = 1; i <= phmm->N; i++) {
                        sum = 0.0;
                        for (j = 1; j <= phmm->N; j++)
                                sum += phmm->A[i][j] *
                                        (phmm->B[j][O[t+1]+1])*beta[t+1][j];  // //O[t+1]+1 代替O[t+1]
                        beta[t][i] = sum;
 
                }
        }
 
        /* 3. Termination */
        *pprob = 0.0;
        for (i = 1; i <= phmm->N; i++)
                *pprob += beta[1][i];
 
}

void BackwardWithScale(HMM *phmm, int T, int *O, double **beta, double *scale, double *pprob)
{
        int     i, j;   /* state indices */
        int     t;      /* time index */
	double sum;
 
 
        /* 1. Initialization */
 
        for (i = 1; i <= phmm->N; i++)
                beta[T][i] = 1.0/scale[T]; 
 
        /* 2. Induction */
 
        for (t = T - 1; t >= 1; t--) {
                for (i = 1; i <= phmm->N; i++) {
			sum = 0.0;
                        for (j = 1; j <= phmm->N; j++)
                        	sum += phmm->A[i][j] * 
					(phmm->B[j][O[t+1]+1])*beta[t+1][j];   //O[t+1]+1 代替O[t+1]
                        beta[t][i] = sum/scale[t];
 
                }
        }
 
}

void ForwardWithScale2(HMM *phmm, int T, int *O, double **alpha, double *scale, double *pprob)
/*  pprob is the LOG probability */
{
	int	i, j; 	/* state indices */
	int	t;	/* time index */

	double sum;	/* partial sum */
	double eps = 10e-70;
	/* 1. Initialization */

	scale[1] = 0.0;	
	for (i = 1; i <= phmm->N; i++) {
		alpha[1][i] = phmm->pi[i]* (phmm->B[i][O[1]+1]);   //O[1]+1 代替O[1]
		scale[1] += alpha[1][i];
	}
	for (i = 1; i <= phmm->N; i++) 
		alpha[1][i] /= (scale[1] + eps); 
	
	/* 2. Induction */
	
	for (t = 1; t <= T - 1; t++) {
		scale[t+1] = 0.0;
		for (j = 1; j <= phmm->N; j++) {
			sum = 0.0;
			for (i = 1; i <= phmm->N; i++) 
				sum += alpha[t][i]* (phmm->A[i][j]); 

			alpha[t+1][j] = sum*(phmm->B[j][O[t+1]+1]);   // //O[t+1]+1 代替O[t+1]
			scale[t+1] += alpha[t+1][j];
		}
		for (j = 1; j <= phmm->N; j++) 
			alpha[t+1][j] /= (scale[t+1] + eps); 
	}

	/* 3. Termination */
	*pprob = 0.0;

	for (t = 1; t <= T; t++)
		*pprob += log(scale[t] + eps);
	
}


void BackwardWithScale2(HMM *phmm, int T, int *O, double **beta, double *scale, double *pprob)
{
    int     i, j;   /* state indices */
    int     t;      /* time index */
	double sum;
	double eps = 10e-70;
 
    /* 1. Initialization */
 
    for (i = 1; i <= phmm->N; i++)
            beta[T][i] = 1.0 / (scale[T] + eps); 

    /* 2. Induction */

    for (t = T - 1; t >= 1; t--)
	{
		for (i = 1; i <= phmm->N; i++) 
		{
			sum = 0.0;
            for (j = 1; j <= phmm->N; j++)
				sum += phmm->A[i][j] * (phmm->B[j][O[t+1]+1]) * beta[t+1][j];  // //O[t+1]+1 代替O[t+1]
			beta[t][i] = sum / (scale[t] + eps);  //使用和alpha同样的尺度因子尺度化
		}
    }
}


void Viterbi(HMM *phmm, int T, int *O, double **delta, int **psi, int *q, double *pprob)
{
	int 	i, j;	/* state indices */
	int  	t;	/* time index */	

	int	maxvalind;
	double	maxval, val;

	/* 1. Initialization  */
	
	for (i = 1; i <= phmm->N; i++) {
		delta[1][i] = phmm->pi[i] * (phmm->B[i][O[1]+1]);  //O[1]+1 代替O[1]
		psi[1][i] = 0;
	}	

	/* 2. Recursion */
	
	for (t = 2; t <= T; t++) {
		for (j = 1; j <= phmm->N; j++) {
			maxval = 0.0;
			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*(phmm->B[j][O[t]+1]);  //O[t]+1 代替O[t]
			psi[t][j] = maxvalind; 

		}
	}

	/* 3. Termination */

	*pprob = 0.0;
	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 ViterbiLog(HMM *phmm, int T, int *O, double **delta, int **psi, int *q, double *pprob)
{
        int     i, j;   /* state indices */
        int     t;      /* time index */
 
        int     maxvalind;
        double  maxval, val;
	double  **biot;

	/* 0. Preprocessing */

	for (i = 1; i <= phmm->N; i++) 
		phmm->pi[i] = log(phmm->pi[i]);
	for (i = 1; i <= phmm->N; i++) 
		for (j = 1; j <= phmm->N; j++) {
			phmm->A[i][j] = log(phmm->A[i][j]);
		}