#include "StdAfx.h"
#include "BmpIO.h"
#include"dft.h"
#include<math.h>
#include <fstream>
#include <iostream>
using namespace std;
double **AllocDouble2D(CSize cs,double initVal)
{
double ** filter=(double **) calloc((unsigned) cs.cx, sizeof(double *));
for (int i = 0; i < cs.cx; i++)
{
filter[i] = (double *) calloc((unsigned) cs.cy, sizeof(double));
if (filter[i] == 0)
return NULL;
for(int j = 0; j < cs.cy ; j++)
{
filter[i][j]=initVal;
}
}
return filter;
}
void FreeDouble2D(double **pFilter, CSize cs)
{
for (int i = 0; i < cs.cx; i++)
free(pFilter[i]);
free(pFilter);
}
double *AllocDouble1D( unsigned int length,double initVal)
{
double * filter=(double *)calloc((unsigned) length, sizeof(double));
return filter;
}
void FreeDouble1D(double *pFilter)
{
free(pFilter);
}
void DFT2D( double ** img,double ** rDout,double **iDout,CSize sz,int shiftFlag )
{
double ** rPara=AllocDouble2D(CSize(sz.cx,sz.cx),0);
double ** iPara=AllocDouble2D(CSize(sz.cx,sz.cx),0);
double w0=2*PI/sz.cx;
double **iTmp=AllocDouble2D(sz,0);
double **rTmp=AllocDouble2D(sz,0);
if(shiftFlag)
{
for(int i=0;i<sz.cx;i++)
for(int j=0; j<sz.cy; j++)
{
if((i+j)%2) img[i][j]=-img[i][j];
}
}
for(int i=0;i<sz.cx;i++)
for(int j=i; j<sz.cx; j++)
{
iPara[i][j]=-sin(w0*i*j);
rPara[i][j]=cos(w0*i*j);
iPara[j][i]=iPara[i][j];
rPara[j][i]=rPara[i][j];
}
for(int k=0;k<sz.cy; k++)
for(int i=0;i<sz.cx;i++)
{
iTmp[i][k]=0;
rTmp[i][k]=0;
for(int j=0; j<sz.cx; j++)
{
rTmp[i][k]+=img[j][k]*rPara[i][j];
iTmp[i][k]+=img[j][k]*iPara[i][j];
}
}
FreeDouble2D(iPara,CSize(sz.cx,sz.cx));
FreeDouble2D(rPara,CSize(sz.cx,sz.cx));
rPara=AllocDouble2D(CSize(sz.cy,sz.cy),0);
iPara=AllocDouble2D(CSize(sz.cy,sz.cy),0);
w0=2*PI/sz.cy;
for(int i=0;i<sz.cy;i++)
for(int j=i; j<sz.cy; j++)
{
iPara[i][j]=-sin(w0*i*j);
rPara[i][j]=cos(w0*i*j);
iPara[j][i]=iPara[i][j];
rPara[j][i]=rPara[i][j];
}
for(int k=0;k<sz.cx; k++)
for(int i=0;i<sz.cy;i++)
{
iDout[k][i]=0;
rDout[k][i]=0;
for(int j=0; j<sz.cy; j++)
{
rDout[k][i]=rDout[k][i]+rTmp[k][j]*rPara[i][j]-iTmp[k][j]*iPara[i][j];
iDout[k][i]=iDout[k][i]+iTmp[k][j]*rPara[i][j]+rTmp[k][j]*iPara[i][j];
}
}
FreeDouble2D(iTmp,sz);
FreeDouble2D(rTmp,sz);
FreeDouble2D(iPara,CSize(sz.cy,sz.cy));
FreeDouble2D(rPara,CSize(sz.cy,sz.cy));
}
double **IDFT2D( double ** rDFT,double **iDFT,CSize sz,int shiftFlag/*=1*/ )
{
double **res=AllocDouble2D(sz,0);
double ** rPara=AllocDouble2D(CSize(sz.cx,sz.cx),0);
double ** iPara=AllocDouble2D(CSize(sz.cx,sz.cx),0);
double w0=2*PI/sz.cx;
double **iTmp=AllocDouble2D(sz,0);
double **rTmp=AllocDouble2D(sz,0);
for(int i=0;i<sz.cx;i++)
for(int j=i; j<sz.cx; j++)
{
iPara[i][j]=sin(w0*i*j);
rPara[i][j]=cos(w0*i*j);
iPara[j][i]=iPara[i][j];
rPara[j][i]=rPara[i][j];
}
for(int k=0;k<sz.cy; k++)
for(int i=0;i<sz.cx;i++)
{
iTmp[i][k]=0;
rTmp[i][k]=0;
for(int j=0; j<sz.cx; j++)
{
iTmp[i][k]+=rDFT[j][k]*iPara[i][j]+iDFT[j][k]*rPara[i][j];
rTmp[i][k]+=rDFT[j][k]*rPara[i][j]-iDFT[j][k]*iPara[i][j];
}
}
FreeDouble2D(iPara,CSize(sz.cx,sz.cx));
FreeDouble2D(rPara,CSize(sz.cx,sz.cx));
rPara=AllocDouble2D(CSize(sz.cy,sz.cy),0);
iPara=AllocDouble2D(CSize(sz.cy,sz.cy),0);
w0=2*PI/sz.cy;
for(int i=0;i<sz.cy;i++)
for(int j=i; j<sz.cy; j++)
{
iPara[i][j]=sin(w0*i*j);
rPara[i][j]=cos(w0*i*j);
iPara[j][i]=iPara[i][j];
rPara[j][i]=rPara[i][j];
}
double iSum=0;
for(int k=0;k<sz.cx; k++)
for(int i=0;i<sz.cy;i++)
{
iDFT[k][i]=0;
rDFT[k][i]=0;
for(int j=0; j<sz.cy; j++)
{
res[k][i]=res[k][i]+rTmp[k][j]*rPara[j][i]-iTmp[k][j]*iPara[j][i];
iSum=iSum+(rTmp[k][j]*iPara[j][i]+iTmp[k][j]*rPara[j][i]);
}
}
FreeDouble2D(iTmp,sz);
FreeDouble2D(rTmp,sz);
FreeDouble2D(iPara,CSize(sz.cy,sz.cy));
FreeDouble2D(rPara,CSize(sz.cy,sz.cy));
for(int i=0;i<sz.cx;i++)
for(int j=0; j<sz.cy; j++)
{
res[i][j]=res[i][j]/(sz.cx*sz.cy);
if(shiftFlag)
{
if((i+j)%2) res[i][j]=-res[i][j];
}
}
return res;
}
double GaussHighPass( int x,int y )
{
double gammaH=1.1,gammaL=0.5;
double D0_2=150*150;
double c=1;
//((x*x+y*y)>20000)?0:1;//
return (gammaH-gammaL)*(1-exp(-c*(1.0*x*x+y*y)/D0_2))+gammaL;
}
int Order(double ** rigion,CSize sz,int d)
{
int m=sz.cx;
int n=sz.cy;
double temp;
int sum=0;
int k;
double *a;
a=AllocDouble1D(m*n,0);
for (int i=0;i<m;i++)
for (int j=0;j<n;j++)
{
a[i*n+j]=rigion[i][j];
}
for(int j=0;j<m*n;j++)
{
k=j;
for (int i=0;i<m*n-j;i++)
if (a[j]>a[k+i])
{
temp=a[j];
a[j]=a[k+i];
a[k+i]=temp;
}
}
for (int i=d/2;i<m*n-d/2;i++)
sum+=a[i];
FreeDouble1D(a);
return sum;
}
int mina(double ** rigion,CSize sz)
{
int m=sz.cx;
int n=sz.cy;
double temp;
int mina=0;
int k;
double *a;
a=AllocDouble1D(m*n,0);
for (int i=0;i<m;i++)
for (int j=0;j<n;j++)
{
a[i*n+j]=rigion[i][j];
}
for(int j=0;j<m*n;j++)
{
k=j;
for (int i=0;i<m*n-j;i++)
if (a[j]>a[k+i])
{
temp=a[j];
a[j]=a[k+i];
a[k+i]=temp;
}
}
mina=a[0];
FreeDouble1D(a);
return mina;
}
int maxa(double ** rigion,CSize sz)
{
int m=sz.cx;
int n=sz.cy;
double temp;
int maxa=0;
int k;
double *a;
a=AllocDouble1D(m*n,0);
for (int i=0;i<m;i++)
for (int j=0;j<n;j++)
{
a[i*n+j]=rigion[i][j];
}
for(int j=0;j<m*n;j++)
{
k=j;
for (int i=0;i<m*n-j;i++)
if (a[j]>a[k+i])
{
temp=a[j];
a[j]=a[k+i];
a[k+i]=temp;
}
}
maxa=a[m*n-1];
FreeDouble1D(a);
return maxa;
}
int meda(double ** rigion,CSize sz)
{
int m=sz.cx;
int n=sz.cy;
double temp;
int meda=0;
int k;
double *a;
a=AllocDouble1D(m*n,0);
for (int i=0;i<m;i++)
for (int j=0;j<n;j++)
{
a[i*n+j]=rigion[i][j];
}
for(int j=0;j<m*n;j++)
{
k=j;
for (int i=0;i<m*n-j;i++)
if (a[j]>a[k+i])
{
temp=a[j];
a[j]=a[k+i];
a[k+i]=temp;
}
}
meda=a[m*n/2];
FreeDouble1D(a);
return meda;
}
int pengding(double ** pIma8,CSize m_sizemybmp,int i,int j,int m,int n)
{
int out=0;
int media=0;
int max=0;
int min=0;
CSize sz(m,n);
double ** rigion;
rigion=AllocDouble2D(sz,0);
int all=0;
for (int k =0 ; k <m; k++)
for (int q =0; q < n; q++)
{
if(i-(m-1)/2<0||j-(n-1)/2<0)
rigion[k][q]=pIma8[i][j];
else
rigion[k][q]=pIma8[i-(m-1)/2+k][j-(n-1)/2+q];
}
media=meda(rigion,sz);
max=maxa(rigion,sz);
min=mina(rigion,sz);
ofstream outfile;
outfile.open("E:\\11.txt",ios::app);
//for (int i =0; i <m; i++)
// for (int j = 0; j <n; j++)
// {
// outfile<<rigion[i][j]<<endl;
// }
outfile<<media<<endl;
outfile<<max<<endl;
outfile<<min<<endl;
outfile.close();
if(max>media&&media>min)
out=1;
else
out=-1;
FreeDouble2D(rigion,sz);
return out;
}