📄 stattest.cpp
字号:
#include "stdafx.h"
#include "Afx.h"
#include "math.h"
#include "AliMath.h"
#include "..\resource.h"
#include <time.h>
#include <stdlib.h>
#include <stdio.h>
#include "StatTest.h"
#ifdef _DEBUG
#define new DEBUG_NEW
#undef THIS_FILE
static char THIS_FILE[] = __FILE__;
#endif
#define PRN_FILE "E:\\Documents from PC\\backup for G\\result\\PRN_data.dat"
#define LUMINANCE(r,g,b) (0.299 * r + 0.587 * g + 0.114 * b )
#define WIDTHBYTES(bits) (((bits) + 31) / 32 * 4)
#define FeatureNum 7
#define MAXRUN 30
double frequency_test(BYTE* s,int n)
{
//1 degree of freedom;
int i,n1;
n1=0;
for(i=0;i<n;i++)
n1 += s[i]; //number of '1's;
return ((double)pow(2*n1-n,2))/((double)n);
}
double serial_test(BYTE* s,int n)
{
//2 degree of freedom;
int i, n0, n1, n00, n01, n10, n11;
n1=0;
for(i=0; i<n; i++)
n1 += s[i];
n0 = n-n1;
n00=n01=n10=n11=0;
for(i=0;i<n-1;i++)
{
if(s[i]==0)
{
if(s[i+1]==0)
n00++;
else
n01++;
}
else
{
if(s[i+1]==0)
n10++;
else
n11++;
}
}
return 4.0*((double)(n00*n00+n01*n01+n10*n10+n11*n11))/((double)(n-1))
- 2.0*((double)(n0*n0+n1*n1))/((double)n)-1;
}
double poker_test(BYTE* s, int n, int m)
{
int i;
int temp;
double result;
if (m<=0 || n/m< pow(2,m)*5) //for 512*512 image,normally m<12,define block length;
return -1;
int k = n/m;
int j;
int M = 1<<m;
int* c=new int[M];
memset(c,0,M*sizeof(int));
for(j=0; j<k; j++)
{
temp=0;
for(i=0;i<m;i++)
temp += (s[j*m+i]<<i);
c[temp]+=1;
}
result=0.0;
for(i=0;i<M;i++)
result+= pow(c[i],2);
delete [] c;
return (result*(double)(M)/(double)(k)) - (double)(k);
//2^m-1 degree of freedom; 4095 for lenna
}
/*
double run_test(BYTE* s,int n)
{
int i;
double run_result;
int* b=new int[MAXRUN]; //runs count for '1', blocks;
int* g=new int[MAXRUN]; //runs count for '0', gaps ;
double* e=new double[MAXRUN]; //expected runs count for '0' or '1' ;
//expected number of i-run for n-length sequence
//n(i)=(n-i+3)/2^(i+2);
//for 512*512 sequence, 1-run numbers equal to 32768 approximately
//so MAXRUN for 512*512 sequence will equal to 17,18
memset(b,0,MAXRUN*sizeof(int));
memset(g,0,MAXRUN*sizeof(int));
int k=-1;
for(i=0;i<n-1;i++)
{
if(s[i]^s[i+1]) //s[i] xor s[i+1]
{
if( s[i] & ((i-k)<MAXRUN) )
b[i-k]++;
if(!s[i] & ((i-k)<MAXRUN) )
g[i-k]++;
k=i;
}
}
k=0;
for(i=1;i<MAXRUN;i++)
{
e[i]=((double)(n-i+3))/((double)pow(2,i+2));
if(e[i]<5)
{
k=i-1; //searching for the largest integer for which e[i]>=5
break;
}
}
run_result=0.0;
for(i=1;i<=k;i++)
run_result+=((double)(pow(b[i]-e[i],2)+pow(g[i]-e[i],2)))/((double)e[i]);
delete [] b;
delete [] g;
delete [] e;
return run_result;
//2*k-1 degree of freedom; 4095 for lenna
}
*/
double run_test(BYTE* s, int n)
{
int i, n0, n1;
double mean, vari;
double run_result = 0.0;
n1=0;
for(i=0;i<n;i++)
n1 += s[i];
n0 = n - n1;
mean = 1 + 2*n0*n1/n;
vari = (mean - 1)*(mean - 2)/(n - 1);
//未完
return run_result;
}
DWORD run_test2(BYTE* s, int n, BYTE* mask, int m)
{
int i, j;
DWORD run_result = 0;
BOOL flag = FALSE;
for (j = 0; j < n - m; j++)
{
for (i = 0; i < m; i++)
{
if (s[j + i]^mask[i])
{
flag = FALSE;
break;
}
else
flag = TRUE;
}
if (flag == TRUE)
run_result ++;
}
return run_result;
}
double autocor_test(BYTE* s,int n,int d)
{
if(2*d>n)
return -1;
int i,temp;
temp=0;
for(i=0;i<n-d;i++)
temp += s[i]^s[i+d];
// return (double)(2*temp+d-n)/sqrt((double)(n-d));
double result = 0.0;
result = (double)((temp-(n-d)/2)/sqrt((double)(n-d)/2));
if (result < 0)
result = -result;
return result;
}
double maurer_test(BYTE* s,int n)
{
const int L=8; //(6,16);
const int Q=10*(1<<L);
const int K=n/L-Q;
if( K < 100*(1<<L) )
{
//AfxMessageBox("序列长度太短,不能执行Maurer's Entropy Test.");
return 0;
}
int* T=new int[1<<L];
memset(T,0,(1<<L)*sizeof(int));
int i,j;
int temp;
for(j=0;j<Q;j++)
{
temp=0;
for(i=0; i<L; i++)
temp += (s[j*L+i]<<i);
T[temp]=j;
}
double maurer_result=0.0;
int* A=new int[K];
for(j=Q;j<Q+K;j++)
{
temp=0;
for(i=0;i<L;i++)
temp+= (s[j*L+i]<<i);
A[j-Q]=j-T[temp];
T[temp]=j;
maurer_result+=log(A[j-Q]);
}
delete [] T;
delete [] A;
return maurer_result/(K*log(2.0));
}
double serial_test1(BYTE* s,int n,int L)
{
int i,j;
int temp;
int k = n/L; //the number of blocks of length L
int M = 1<<L;
int* c=new int[M];
memset(c, 0, M*sizeof(int));
for(j=0;j<k;j++)
{
temp=0;
for(i=0;i<L;i++)
temp += (s[j*L+i]<<i);
c[temp]++;
}
double result=0.0;
for(i=0;i<M;i++)
result+=pow(c[i]-((double)(n))/((double)(L*M)),2);
delete [] c;
return result*((double)(L*M))/((double)(n));
//2^L-1 degree of freedom; 4095 for lenna
}
double phai_serial_test(BYTE* s, int n, int m)
{
if( m<=0 )
return 0.0; //(for m=0 or m=-1 return -1;else for m<-1,return a error code);
int i,j;
int temp;
int M = 1<<m;
int* c=new int[M];
memset(c, 0, M*sizeof(int) );
temp=0;
for(i=0;i<m;i++)
temp += (s[i]<<i);
c[temp]++;
for(j=m;j<n;j++)
{
temp >>=1;
temp+= (s[j]<<(m-1));
c[temp]++;
}
double result=0.0;
for(i=0; i<M; i++)
result+=pow(c[i]-((double)(n))/((double)(M)),2);
delete [] c;
return result*((double)(M))/((double)(n));
//2^(m-1) degree of freedom; 4095 for lenna
}
double serial_test2(BYTE* s,int n)
{
double PhaiValue[FeatureNum+1];
int l;
double TestValue;
for(l=0 ; l<FeatureNum+1 ; l++)
PhaiValue[l]=phai_serial_test(s, n, l+1);
TestValue = 0.0;
for(l=0 ; l<FeatureNum ; l++)
TestValue += (PhaiValue[l+1]-PhaiValue[l]) / (1<<l);
return TestValue;
}
void GetLSB(LPBYTE pImg, LPBYTE pLSB, int nCount)
{
int i;
for(i=0; i<nCount; i++)
pLSB[i] = pImg[i] & 1;
}
void Embed(LPBYTE pImg, LPBYTE pData, int nCount1, int nCount2, int seed)
{
if( nCount2==0 || nCount2>nCount1)
return;
int i;
double curPos;
double alpha=(double)(nCount1-nCount2)*2.0 / ((double)(nCount2) * (double)(RAND_MAX));
int nCurPos;
srand((unsigned)time(NULL)+seed);
curPos=0.0;
for( i=0 ; i<nCount2 ; i++ )
{
curPos += (rand()*alpha+1);
if(curPos>nCount1)
break;
nCurPos=(int)(curPos)-1;
pImg[nCurPos] = (pImg[nCurPos]>>1)<<1;
pImg[nCurPos] += pData[i]; //Embedding the data to LSB;
}
}
void Embed1(LPBYTE pImg, LPBYTE pData, int nCount1, int nCount2, int seed)
{
if( nCount2==0 || nCount2>nCount1)
return;
int i;
double curPos;
double alpha=(double)(nCount1-nCount2)*2.0 / ((double)(nCount2) * (double)(RAND_MAX));
int nCurPos;
srand((unsigned)time(NULL)+seed);
curPos=0.0;
for( i=0 ; i<nCount2 ; i++ )
{
curPos += (rand()*alpha+1);
if(curPos>nCount1)
break;
nCurPos=(int)(curPos)-1;
/*
pImg[nCurPos] = (pImg[nCurPos]>>1)<<1;
pImg[nCurPos] += pData[i]; //Embedding the data to LSB;*/
/*
if( (pImg[nCurPos]&1)^pData[i] )
{
if(rand()&1 || pImg[nCurPos]==0)
pImg[nCurPos]++;
else
pImg[nCurPos]--;
}*/
if( (pImg[nCurPos]&1)^pData[i] )
{
if( pImg[nCurPos]!=255)
pImg[nCurPos]++;
else
pImg[nCurPos]--;
}
}
}
int Prep_Sequence(BYTE* s,int n)
{
//根据游程的统计分布规律，去除超长的游程
//对长度为n的序列，最长游程为ln(n)/ln(2)-2
//对512*512序列，最长游程为16
//对1024*1024序列，最长游程为18
//因此可以定义最大游程为20
//序列中超过MAXRUN长度的连续0或者连续1将被去掉,分别代之以单个的0或者1;
BYTE* lpByte=new BYTE[n];
int i;
int k=-1;
int q=0;
int nextbit;
for(i=0; i<n; i++)
{
nextbit = (i==n-1)? (1^s[i]) : s[i+1];
if(s[i]^nextbit)
{
if( (i-k) < MAXRUN)
{
memcpy(lpByte+q, s+k+1, i-k);
q += (i-k);
}
else
lpByte[q++]=s[i];
k=i;
}
}
memcpy(s,lpByte,q);
delete lpByte;
return q;
}
double Chi_Test_His(int* pnHis)
{
double chi=0;
int dnf=0;
for(int i=0; i<128; i++)
if(pnHis[2*i]+pnHis[2*i+1]>5)
{
chi+=pow((double)(pnHis[2*i]-pnHis[2*i+1])/2.0,2)*2.0 / (double)(pnHis[2*i]+pnHis[2*i+1]);
dnf++;
}
return ChiProb(chi,dnf-1);
}
LPBYTE ReadPRNData(int nCount)
{
LPBYTE pLSB = new BYTE[nCount];
LPBYTE pData= new BYTE[nCount];
FILE* fPrn=fopen(PRN_FILE,"rb");
if(fPrn!=NULL)
{
fread( pData , 1 , nCount/8+1, fPrn);
fclose(fPrn);
}
else
{
// AfxMessageBox(IDS_PRN_FILE_ERROR,MB_ICONINFORMATION);
return NULL;
}
int i,j,nIndex;
for(i=0 ; i<nCount/8+1 ; i++)
for(j=0;j<8;j++)
{
nIndex=i*8+j;
if( nIndex > nCount-1) break;
pLSB[nIndex] = (pData[i]>>j) & 1;
}
delete pData;
return pLSB;
}
double* Prog_Test(double (*func)(LPBYTE,int), LPBYTE s, int n, int steps)
{
if(steps <= 0 )
return NULL;
double* result = new double[steps];
int n1;
for(int i=0 ; i<steps; i++)
{
n1 = i*n/steps;
result[i] = (*func)(s, n1);
}
return result;
}
//----The following codes are used for RS Analysis proposed by Jiri Fridrich------------------------------------------------------------------------------------------------
double Smoothness(int* pData,int nCount)
{
double val=0.0;
for(int i=0; i<nCount-1; i++)
val += fabs(pData[i]-pData[i+1]);
return val;
}
void Flip(int* pData,int* pMask, int nCount)
{
for(int i=0; i<nCount; i++)
{
if(pMask[i]==0) continue;
if(pMask[i]==-1) pData[i]++;
if(pData[i] & 1) pData[i]--;
else pData[i]++;
if(pMask[i]==-1) pData[i]--;
}
}
void GetRS(LPBYTE pImg,int* pMask,int nCount,int nMask,int* r, int* r1, int* s, int* s1)
{
//pMask算子应该是非负的,取值为0,1
int nGroup=nCount/nMask;
int i, j;
int* pData = new int[nMask];
int* pData1 = new int[nMask];
int* pMask1 = new int[nMask];
double smooth,smooth1,smooth2;
for(j=0; j<nMask; j++)
pMask1[j] = -pMask[j];
int rCount,sCount,r1Count,s1Count;
rCount=sCount=r1Count=s1Count=0;
for(i=0; i<nGroup; i++)
{
for(j=0; j<nMask; j++)
pData[j]=pImg[i*nMask+j];
memcpy(pData1, pData,nMask*sizeof(int));
smooth = Smoothness(pData, nMask);
Flip(pData, pMask, nMask);
smooth1= Smoothness(pData, nMask);
if(smooth1>smooth)
rCount++;
if(smooth1<smooth)
sCount++;
Flip(pData1, pMask1, nMask);
smooth2= Smoothness(pData1, nMask);
if(smooth2>smooth)
r1Count++;
if(smooth2<smooth)
s1Count++;
}
*r = rCount;
*s = sCount;
*r1 = r1Count;
*s1 = s1Count;
delete pData;
delete pData1;
delete pMask1;
return;
}
//----Image Quality Evaluation------------------------------------------------------------------------------------------------
double* ComputeLuminance(LPBYTE pImg, int nWidth, int nHeight, int nComponents)
{
int i,j,rowstride;
double r,g,b;
rowstride = WIDTHBYTES(nWidth*nComponents*8);
double* pLum = new double[nWidth*nHeight];
for(i=0; i<nHeight; i++)
for(j=0; j<nWidth; j++)
{
if(nComponents==1)
pLum[i*nWidth+j]=(double)(pImg[i*rowstride+j]);
else
{
r = (double)pImg[i*rowstride+j*3+2];
g = (double)pImg[i*rowstride+j*3+1];
b = (double)pImg[i*rowstride+j*3 ];
pLum[i*nWidth+j] = LUMINANCE(r,g,b);
}
}
return pLum;
}
double MaxLuminance(double* pLum, int nWidth, int nHeight)
{
int i,j;
double max = 0.0;
double val;
for(i=0; i<nHeight; i++)
for(j=0; j<nWidth; j++)
{
val = pLum[i*nWidth+j];
if(val > max)
max = val;
}
return (max);
}
double Norm2(double* pLum, int nWidth, int nHeight)
{
int i,j;
double sum = 0.0;
double val;
for(i=0; i<nHeight; i++)
for(j=0; j<nWidth; j++)
{
val = pLum[i*nWidth+j];
sum += (val*val);
}
return (sum);
}
/*----------------------------------------------------------------------------
*/
double DiffNorm2(double* pLum1, double* pLum2, int nWidth, int nHeight)
{
int i,j;
double sum = 0.0;
double val1, val2;
for(i=0; i<nHeight; i++)
for(j=0; j<nWidth; j++)
{
val1 = pLum1[i*nWidth+j];
val2 = pLum2[i*nWidth+j];
sum += pow((val1-val2), 2);
}
return (sum);
}
/*----------------------------------------------------------------------------
* Signal to noise ration (SNR)
*
* N 2
* Sum X
* i=1 i
* SNR = 10 log_10 ----------------
* N 2
* Sum (X - X')
* i=1 i i
* where X_i and and X'_i are the original and modified pixel values and N the
* total number of pixels in the image. SNR is used to measure the signal
* quality.
*/
double SNR(double* pLum1, double* pLum2, int nWidth, int nHeight)
{
return (10 * log10(Norm2(pLum1, nWidth, nHeight) / DiffNorm2(pLum1, pLum2, nWidth, nHeight)));
}
/*----------------------------------------------------------------------------
* Peak Signal to noise ration (PSNR)
*
* 2
* N max X
* i i
* PSNR = 10 log_10 ------------------
* N 2
* Sum (X - X')
* i=1 i i
* where X_i and and X'_i are the original and modified pixel values and N the
* total number of pixels in the image. SNR is used to measure the visual
* quality, even though it is not monotonically related to subjective visual
* quality.
*/
double PSNR(double* pLum1, double* pLum2, int nWidth, int nHeight)
{
double max;
max = 255.0; //MaxLuminance(pLum1, nWidth, nHeight);
double diff = DiffNorm2(pLum1, pLum2, nWidth, nHeight);
if(diff==0)
return 10000.0;
return (20 * log10(max) + 10 * log10(nWidth * nHeight) - 10 * log10(diff));
}
//----Image Quality Evaluation------------------------------------------------------------------------------------------------
int hamming_weight(unsigned char a)
{
int dis=0;
for(int i=0; i<8; i++)
if(a & (1<<i))
dis++;
return dis;
}⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -