www.pudn.com > Gesture[20040824].rar > Background.cpp
// Background.cpp: implementation of the CBackground class.
//
//////////////////////////////////////////////////////////////////////
#include "stdafx.h"
#include "Gesture.h"
#include "Background.h"
#include "Image.h"
#include "math.h"
#ifdef _DEBUG
#undef THIS_FILE
static char THIS_FILE[]=__FILE__;
#define new DEBUG_NEW
#endif
static double threshold[3] = {0.5,0.5,0.5};
#define pi 3.1415926535897
double (*P) [3];
double (*L) [3];
double* S;
double (*P2) [3];
//double* xx;
//////////////////////////////////////////////////////////////////////
// Construction/Destruction
//////////////////////////////////////////////////////////////////////
CBackground::CBackground(CCamAvi *avi)
{
BkgAvi = avi;
Th_diff = 10;
m_BkgImage = new CImage;
m_BkgType = GaussianAndThreshold;
m_bBkgReady = false;
displayflag = false;
testflag = false;
imageshow1 = new CImage;
imageshow2 = new CImage;
imageshow3 = new CImage;
}
CBackground::~CBackground()
{
}
IplImage* CBackground::Thresholding_GrayLevel(IplImage *image)
{
CvSize imsize;
imsize.width = image->width;
imsize.height = image->height;
IplImage *grayimage;
if (image->nChannels>1) {
grayimage = cvCreateImage(imsize, image->depth, 1);
iplColorToGray(image, grayimage);
} else grayimage = image;
IplImage *tempimage = cvCreateImage(imsize, image->depth, 1);
IplImage *diffimage = cvCreateImage(imsize, image->depth, 1);
switch ( m_BkgType ) {
case GaussianAndThreshold :
cvAbsDiff(grayimage, GrayBkg, tempimage);
cvThreshold(tempimage, diffimage, Th_diff, 255, CV_THRESH_BINARY);
break;
default :
return NULL;
}
cvReleaseImage(&tempimage);
if (image->nChannels>1) cvReleaseImage(&grayimage);
return diffimage;
}
void CBackground::ConvertColorBkg2GrayBkg()
{
CvSize imsize;
imsize.width = m_BkgImage->Width();
imsize.height = m_BkgImage->Height();
IplImage *image = m_BkgImage->GetImage();
GrayBkg = cvCreateImage( imsize, image->depth, 1);
if ( image->nChannels > 1 )
iplColorToGray( image, GrayBkg);
else cvCopyImage( image, GrayBkg);
}
void CBackground::ExtractBackgroundFromVideo()
{
CvSize imsize;
int length = BkgAvi->GetAviLength();
IplImage* gray[2000];//怎么用动态表示？
BkgAvi->SetReaderPosition(0);
//读第一帧
imsize.width = BkgAvi->GetFramePointer()->Width();
imsize.height = BkgAvi->GetFramePointer()->Height();
BkgAvi->GetFrameFromAvi();
for (int i = 0;i < length;i++){
BkgAvi->GetFrameFromAvi();
gray[i] = cvCreateImage(imsize,BkgAvi->GetFramePointer()->GetImage()->depth,1);
iplColorToGray(BkgAvi->GetFramePointer()->GetImage(), gray[i]);
BkgAvi->NextFrame();
}
ExtractMoGBackground(gray);
}
void CBackground::ExtractGaussianBackground()
{
}
void CBackground::ExtractMoGBackground(IplImage** gray)
{
int length = BkgAvi->GetAviLength();
int h = BkgAvi->GetFramePointer()->Height();
int w = BkgAvi->GetFramePointer()->Width();
m_MoG = new struct MoGParameters [h*w];
unsigned char *pixelprocess= new unsigned char [length];
for (int a= 0;a < h;a++){
for (int b= 0;b < w;b++){
for (int c = 0;c < length;c++){
//得到一个象素过程的数组，然后操作
pixelprocess[c] = gray[c]->imageData[a*w+b];
}
P = (double (*)[3]) (new double[length*3]);
L = (double (*)[3]) (new double[length*3]);
S = new double [length];
P2 = (double (*)[3]) (new double[length*3]);
EM(pixelprocess,a,b);
}
}
imageshow1->Create(w,h,8);
imageshow2->Create(w,h,8);
imageshow3->Create(w,h,8);
for (int d = 0;d < h*w;d++)
imageshow1->GetImage()->imageData[d] = m_MoG[d].miu[0];
for (int e = 0;e < h*w;e++)
imageshow2->GetImage()->imageData[e] = m_MoG[e].miu[1];
for (int f = 0;f < h*w;f++)
imageshow3->GetImage()->imageData[f] = m_MoG[f].miu[2];
//释放所有指针
for (int j = 0;j < length;j++){
cvReleaseImage(&gray[j]);
}
displayflag = true;
}
void CBackground::EM(unsigned char* pixelprocess,int x_co,int y_co)
{
//初始化
int height = BkgAvi->GetFramePointer()->Height();
int width = BkgAvi->GetFramePointer()->Width();
double ome[3];
double miu[3];
double xig[3];
for (int a = 0;a < 3; a++){
ome[a] = 0.333333333333;
xig[a] = 900;
}
int length = BkgAvi->GetAviLength();
miu[0] = pixelprocess[length/4]; miu[1] = pixelprocess[length/2]; miu[2] = pixelprocess[3*length/4];
//调用EM23
EM23(pixelprocess,ome,miu,xig,x_co,y_co,height,width);
}
void CBackground::EM23(unsigned char* pixelprocess,double* ome,double* miu,double* xig,int x_co,int y_co,int height,int width)
{
double miutemp[3];
double temp[3];
bool flag;
double fenzi[3],fenmu[3];
double T[3];
double T2[3];
bool ok = false;
int max,mid,min;
int length = BkgAvi->GetAviLength();
flag = false;
for (int m = 0;m < 3;m++){
miutemp[m] = 255;
T[m] = 0;
T2[m] = 0;
fenzi[m] = 0;
fenmu[m] = 0;
}
//K-mean算法，求得带入EM的初值
double muu[3],muutemp[3];
muu[0] = 60; muu[1] = 120; muu[2] = 180;
bool kflag = false;
int leibie[256],count1=0,count2=0,count3=0;
double tem1=0,tem2=0,tem3=0,te1=0,te2=0,te3=0;
while(!kflag){
muutemp[0] = muu[0];muutemp[1] = muu[1];muutemp[2] = muu[2];
for (int w=0;w<256;w++){
leibie[w] = KJudge(w,muu[0],muu[1],muu[2]);
}
for (int v=0;v<256;v++){
if (leibie[v]==1){tem1 += v;count1 += 1;}
else{
if (leibie[v]==2){tem2 +=v;count2 += 1;}
else{
if (leibie[v]==3){tem3 +=v;count3 += 1;}
}
}
}
muu[0] = tem1/count1;muu[1] = tem2/count2;muu[2] = tem3/count3;
if (((muu[0]-muutemp[0])<=1)&&((muu[1]-muutemp[1])<=1)&&((muu[2]-muutemp[2])<=1)){
kflag = true;
}
}
for (int u=0;u<256;u++){
if (leibie[u]==1){te1 += pow((u-muu[0]),2);}
else{
if (leibie[u]==2){te2 += pow((u-muu[1]),2);}
else{
if (leibie[u]==3){te3 += pow((u-muu[2]),2);}
}
}
}
xig[0] = te1/(count1-1);xig[1] = te2/(count2-1);xig[2] = te3/(count3-1);
miu[0] = muu[0];miu[1] = muu[1];miu[2] = muu[2];
//如果不小于阈值，就进行操作，否则返回
while (!flag){
for (int f = 0;f < length;f++){
P[f][0] = (exp(-(pixelprocess[f]-miu[0])*(pixelprocess[f]-miu[0])/(2*xig[0])))/((pow(2*pi,0.5))*(pow(xig[0],0.5)));
P[f][1] = (exp(-(pixelprocess[f]-miu[1])*(pixelprocess[f]-miu[1])/(2*xig[1])))/((pow(2*pi,0.5))*(pow(xig[1],0.5)));
P[f][2] = (exp(-(pixelprocess[f]-miu[2])*(pixelprocess[f]-miu[2])/(2*xig[2])))/((pow(2*pi,0.5))*(pow(xig[2],0.5)));
L[f][0] = ome[0]*P[f][0];
L[f][1] = ome[1]*P[f][1];
L[f][2] = ome[2]*P[f][2];
S[f] = L[f][0] + L[f][1] + L[f][2];
P2[f][0] = L[f][0]/S[f];
P2[f][1] = L[f][1]/S[f];
P2[f][2] = L[f][2]/S[f];
T[0] += P2[f][0];
T[1] += P2[f][1];
T[2] += P2[f][2];
T2[0] += P2[f][0]*pixelprocess[f];
T2[1] += P2[f][1]*pixelprocess[f];
T2[2] += P2[f][2]*pixelprocess[f];
}
//第三步
//ome[0] = T[0]/length;
//ome[1] = T[1]/length;
//ome[2] = T[2]/length;
miu[0] = T2[0]/T[0];
miu[1] = T2[1]/T[1];
miu[2] = T2[2]/T[2];
for (int h = 0;h < length;h++){
fenzi[0] += ((pixelprocess[h]-miu[0])*(pixelprocess[h]-miu[0])*P2[h][0]);
fenzi[1] += ((pixelprocess[h]-miu[1])*(pixelprocess[h]-miu[1])*P2[h][1]);
fenzi[2] += ((pixelprocess[h]-miu[2])*(pixelprocess[h]-miu[2])*P2[h][2]);
fenmu[0] += P2[h][0];
fenmu[1] += P2[h][1];
fenmu[2] += P2[h][2];
}
xig[0] = fenzi[0]/fenmu[0];
xig[1] = fenzi[1]/fenmu[1];
xig[2] = fenzi[2]/fenmu[2];
//判断flag
for (int j = 0;j < 3;j++){
temp[j] = abs(miu[j] - miutemp[j]);
}
for (int l = 0;l < 3;l++){
if (temp[l] <= threshold[l]){flag = true;}
else {flag = false; break;}
}
//记录上一次的结果
for (int e = 0;e < 3;e++){
miutemp[e] = miu[e];
}
}
//记录miu[i]，需要归类
if (miu[0]>=miu[1]){
if (miu[1]>=miu[2]){max=0;mid=1;min=2;}
else{
if (miu[0]>=miu[2]){max=0;mid=2;min=1;}
else{max=2;mid=0;min=1;}
}
}
else{
if (miu[0]>=miu[2]){max=1;mid=0;min=2;}
else{
if (miu[1]>=miu[2]){max=1;mid=2;min=0;}
else{max=2;mid=1;min=0;}
}
}
m_MoG[x_co*width+y_co].miu[0] = miu[max];
m_MoG[x_co*width+y_co].miu[1] = miu[mid];
m_MoG[x_co*width+y_co].miu[2] = miu[min];
delete(P);
delete(L);
delete(S);
delete(P2);
//delete(x);
}
int CBackground::KJudge(int w,double x,double y,double z)
{
double tt1,tt2,tt3;
tt1 = abs(w-x);
tt2 = abs(w-y);
tt3 = abs(w-z);
if (tt1 <= tt2){
if (tt2 <= tt3){return 1;}
if (tt1 <= tt3){return 1;}
if (tt3 < tt1){return 3;}
}
else{
if (tt1 <= tt3){return 2;}
if (tt2 <= tt3){return 2;}
if (tt3 < tt2){return 3;}
}
}
int CBackground::Judge(double ranran,double* probb)
{
for (int a = 0; a < 256; a++){
if (ranran == probb[a]) { return a; break; }
if ((ranran > probb[a])&&(ranran < probb[a+1])) { return (a+1); break; }
}
}
void CBackground::Test()
{
//给定mu和xg
double mu[3],xg[3],om[3];
mu[0] = 50; mu[1] = 130; mu[2] = 190;
xg[0] = 200,xg[1] = 300,xg[2] = 400;
om[0] = 0.333333333333,om[1] = 0.333333333333,om[2] = 0.333333333333;
double temp0=0,temp1=0,temp2=0;
double ran=0;
double EMP[12000][3],EML[12000][3],EMP2[12000][3],EMS[12000],T[3],T2[3],miutemp[3],temp[3],fenzi[3],fenmu[3];
bool flag=false;
double ome[3],miu[3],xig[3];
//利用公式生成Gauss分布的一个象素过程
for (int a = 0; a < 256; a++){
temp0 = (om[0]*exp(-pow((a-mu[0]),2)/(2*xg[0])))/(pow(2*pi,0.5)*pow(xg[0],0.5));
temp1 = (om[1]*exp(-pow((a-mu[1]),2)/(2*xg[1])))/(pow(2*pi,0.5)*pow(xg[1],0.5));
temp2 = (om[2]*exp(-pow((a-mu[2]),2)/(2*xg[2])))/(pow(2*pi,0.5)*pow(xg[2],0.5));
gauss[a] = temp0+temp1+temp2;
}
//得到概率分布
probb[0] = gauss[0];
for (int b = 1; b < 256; b++){
probb[b] = probb[b-1] + gauss[b];
}
//利用随机数得到象素过程的每个点值
srand( (unsigned)time( NULL ) );
for (int c = 0; c < 12000; c++){
ran = (rand()*probb[255])/RAND_MAX;//ran是0-1之间
pixpro[c] = Judge(ran,probb);
}
//统计得到0到255的象素值的象素的个数
for (int y=0;y<256;y++){
geshu[y] = 0;
}
for (int x=0;x<12000;x++){
geshu[pixpro[x]] += 1;
}
//K-mean算法，求得带入EM的初值
double muu[3],muutemp[3];
muu[0] = 60; muu[1] = 120; muu[2] = 180;
bool kflag = false;
int leibie[256],count1=0,count2=0,count3=0;
double tem1=0,tem2=0,tem3=0,te1=0,te2=0,te3=0;
while(!kflag){
muutemp[0] = muu[0];muutemp[1] = muu[1];muutemp[2] = muu[2];
for (int w=0;w<256;w++){
leibie[w] = KJudge(w,muu[0],muu[1],muu[2]);
}
for (int v=0;v<256;v++){
if (leibie[v]==1){tem1 += v;count1 += 1;}
else{
if (leibie[v]==2){tem2 +=v;count2 += 1;}
else{
if (leibie[v]==3){tem3 +=v;count3 += 1;}
}
}
}
muu[0] = tem1/count1;muu[1] = tem2/count2;muu[2] = tem3/count3;
if (((muu[0]-muutemp[0])<=1)&&((muu[1]-muutemp[1])<=1)&&((muu[2]-muutemp[2])<=1)){
kflag = true;
}
}
for (int u=0;u<256;u++){
if (leibie[u]==1){te1 += pow((u-muu[0]),2);}
else{
if (leibie[u]==2){te2 += pow((u-muu[1]),2);}
else{
if (leibie[u]==3){te3 += pow((u-muu[2]),2);}
}
}
}
xig[0] = te1/(count1-1);xig[1] = te2/(count2-1);xig[2] = te3/(count3-1);
miu[0] = muu[0];miu[1] = muu[1];miu[2] = muu[2];
//带入EM计算，得到新的mu和xg
for (int z = 0;z < 3; z++){
ome[z] = 0.333333333333;
}
for (int d = 0;d < 3;d++){
T[d] = 0;
T2[d] = 0;
fenzi[d] = 0;
fenmu[d] = 0;
}
while (!flag){
for (int f = 0;f < 12000;f++){
EMP[f][0] = (exp(-pow((pixpro[f]-miu[0]),2)/(2*xig[0])))/((pow(2*pi,0.5))*(pow(xig[0],0.5)));
EMP[f][1] = (exp(-pow((pixpro[f]-miu[1]),2)/(2*xig[1])))/((pow(2*pi,0.5))*(pow(xig[1],0.5)));
EMP[f][2] = (exp(-pow((pixpro[f]-miu[2]),2)/(2*xig[2])))/((pow(2*pi,0.5))*(pow(xig[2],0.5)));
EML[f][0] = ome[0]*EMP[f][0];
EML[f][1] = ome[1]*EMP[f][1];
EML[f][2] = ome[2]*EMP[f][2];
EMS[f] = EML[f][0] + EML[f][1] + EML[f][2];
EMP2[f][0] = EML[f][0]/EMS[f];
EMP2[f][1] = EML[f][1]/EMS[f];
EMP2[f][2] = EML[f][2]/EMS[f];
T[0] += EMP2[f][0];
T[1] += EMP2[f][1];
T[2] += EMP2[f][2];
T2[0] += EMP2[f][0]*pixpro[f];
T2[1] += EMP2[f][1]*pixpro[f];
T2[2] += EMP2[f][2]*pixpro[f];
}
miu[0] = T2[0]/T[0];
miu[1] = T2[1]/T[1];
miu[2] = T2[2]/T[2];
for (int h = 0;h < 12000;h++){
fenzi[0] += (pow((pixpro[h]-miu[0]),2)*EMP2[h][0]);
fenzi[1] += (pow((pixpro[h]-miu[1]),2)*EMP2[h][1]);
fenzi[2] += (pow((pixpro[h]-miu[2]),2)*EMP2[h][2]);
fenmu[0] += EMP2[h][0];
fenmu[1] += EMP2[h][1];
fenmu[2] += EMP2[h][2];
}
xig[0] = fenzi[0]/fenmu[0];
xig[1] = fenzi[1]/fenmu[1];
xig[2] = fenzi[2]/fenmu[2];
for (int j = 0;j < 3;j++){
temp[j] = abs(miu[j] - miutemp[j]);
}
for (int l = 0;l < 3;l++){
if (temp[l] <= threshold[l]){flag = true;}
else {flag = false; break;}
}
//记录上一次的结果
for (int e = 0;e < 3;e++){
miutemp[e] = miu[e];
}
}
EMmu[0] = miu[0];
EMmu[1] = miu[1];
EMmu[2] = miu[2];
for (int aa = 0; aa < 256; aa++){
temp0 = (ome[0]*exp(-pow((aa-miu[0]),2)/(2*xig[0])))/(pow(2*pi,0.5)*pow(xig[0],0.5));
temp1 = (ome[1]*exp(-pow((aa-miu[1]),2)/(2*xig[1])))/(pow(2*pi,0.5)*pow(xig[1],0.5));
temp2 = (ome[2]*exp(-pow((aa-miu[2]),2)/(2*xig[2])))/(pow(2*pi,0.5)*pow(xig[2],0.5));
gauss2[aa] = temp0+temp1+temp2;
}
//置标志，显示出来
testflag = true;
//释放所有资源
}
//找出每个象素点属于哪个类别的
/*for (int g = 0;g < length;g++){
double temp0 = P2[0];
double temp1 = abs(int(pixelprocess[g]-miu[1]));
double temp2 = abs(int(pixelprocess[g]-miu[2]));
//判断最大的，中间的，最小的
if (temp0 >= temp1){
if (temp1 >= temp2) {max=0;mid=1;min=2;}
else{
if (temp0 >= temp2) {max=0;mid=2;min=1;}
else {max=2;mid=0;min=1;}
}
}
else{//temp1>temp0
if (temp0 >= temp2) {max=1;mid=0;min=2;}
else{
if (temp1 >= temp2) {max=1;mid=2;min=0;}
else {max=2;mid=1;min=0;}
}
}
xx[g][0] = max;
xx[g][1] = mid;
xx[g][2] = min;
}*/
/*
if (miu[0]>=miu[1]){
if (miu[1]>=miu[2]){max=0;mid=1;min=2;}
else{
if (miu[0]>=miu[2]){max=0;mid=2;min=1;}
else{max=2;mid=0;min=1;}
}
}
else{
if (miu[0]>=miu[2]){max=1;mid=0;min=2;}
else{
if (miu[1]>=miu[2]){max=1;mid=2;min=0;}
else{max=2;mid=1;min=0;}
}
}
*/