📄 corrcoff.cpp
字号:
// CorrCoff.cpp: implementation of the CCorrCoff class.
//
//////////////////////////////////////////////////////////////////////
#include "CorrCoff.h"
#include "math.h"
//////////////////////////////////////////////////////////////////////
// Construction/Destruction
//////////////////////////////////////////////////////////////////////
CCorrCoff::CCorrCoff()
{
}
CCorrCoff::CCorrCoff(int x,int y,int nx,int ny)
{
m_nSBoxX=x;
m_nSBoxY=y;
m_nBBoxX=nx;
m_nBBoxY=ny;
m_fLowDbz=10;
m_fLowPercengtage=0.80f;
}
CCorrCoff::~CCorrCoff()
{
}
void CCorrCoff::Compute_Coff(int Z11[][1000],int Z22[][1000], int Xm/*所取空间j的格点数*/,int Yn/*所取空间i的格点数*/,float fT)
{
int i,j; //所选资料的格点索引
// int m=0,n=0; //记录盒子的位置
FILE *file;
file=fopen("f:\\wave\\Coff.dat","wt");
float nX1,nY1; ////////
m_fT=fT;
for(i=0;i<=Yn;i=i+m_nSBoxY) ///i表示y方向，j表示x方向
{
for(j=0;j<=Xm;j=j+m_nSBoxX)
{
nX1=-1;
nY1=-1;
Compute_Corr(Z11,Z22,j,i,Xm,Yn,nX1,nY1,file);//i,j格点的起始位置xm 513
// printf("%d %d\n",nY1,nX1);
////////////////////////////////////////////////////////////////
//计算每个格点上U、V风速
if(nX1==-1||nY1==-1) ////////
{
m_U[i/m_nSBoxX][j/m_nSBoxY]=0;
m_V[i/m_nSBoxX][j/m_nSBoxY]=0;
continue;
}
m_U[i/m_nSBoxY][j/m_nSBoxX]=( nX1-(j+m_nSBoxX/2))*0.25f*1000/m_fT;
// float dd=m_U[i/m_nSBoxX][j/m_nSBoxY];
m_V[i/m_nSBoxY][j/m_nSBoxX]=( nY1-(i+m_nSBoxY/2))*0.25f*1000/m_fT;
///////调试用////////////////////////////////////////////////
/*
if( m_U[i/m_nSBoxX][j/m_nSBoxY]==0||
m_V[i/m_nSBoxX][j/m_nSBoxY]==0)
{
int ppp=0;
}*/
// n=n+1;
}
// m=m+1;
}
fclose(file);
}
//////////////////////////////////////////////////////////////////////////////////////////////////////
////搜索方向
void CCorrCoff::Compute_Corr(int Z11[][1000], int Z22[][1000],int x,int y,int Xm,int Yn,float&nX0,float&nY0,FILE *file)
{
int m;
// int ii=(m_nBBoxX-m_nSBoxX)/(2*m_nSBoxX); //在x方向要搜索的次数
// int jj=(m_nBBoxY-m_nSBoxY)/(2*m_nSBoxY); //在y方向要搜索的次数,
//大盒子与小盒子有同样的中心点(x+m_nSBoxX/2,y+m_nSBoxY/2),(x,y)为小盒子左上角坐标
float ss=0.0f;
int mm1,nn1,mm2,nn2;
short *temp1=new short[1000];
short *temp2=new short[1000];
/////////////////对边界进行处理///////////////////////////////////////////////////////
// mm1=x-ii*m_nSBoxX;
mm1=(x+m_nSBoxX/2)-m_nBBoxX/2; ///这两句等价
// mm1=x
if(mm1<0)mm1=0;
// mm2=x+(ii+1)*m_nSBoxX;
mm2=(x+m_nSBoxX/2)+m_nBBoxX/2;
if(mm2>Xm)mm2=Xm;
// nn1=y-jj*m_nSBoxY;
nn1=(y+m_nSBoxY/2)-m_nBBoxY/2;
if(nn1<0)nn1=0;
// nn2=y+(jj+1)*m_nSBoxY;
nn2=(y+m_nSBoxY/2)+m_nBBoxY/2;
if(nn2>Yn)nn2=Yn;
float s1=0.0f;
//////////计算每一个盒子的相关系数////////////////////////////////////////////////////////
int l,k;
int N=0;
for(l=0;l<m_nSBoxX;l++)
{
for(k=0;k<m_nSBoxY;k++)
{
if((Z11[x+l][y+k])>=m_fLowDbz) ///当该地方有反照率值的时候，记下该值
{
s1=s1+(Z11[x+l][y+k]);
N=N+1;
}
}
}
if(N>=m_fLowPercengtage*m_nSBoxX*m_nSBoxY) ///如果所选点数超过了一定值就选出
{
s1=s1/N;
///////计算时间2的平均值////////////////////////////////////////////////////////
int ll,kk;
float rho=0;
// 下面开始搜寻///////////////////////////////////////////////
for(ll=mm1;ll<mm2;ll++)//=ll+m_nSBoxX) //i,j为空间格点数
{
for(kk=nn1;kk<nn2;kk++)//=kk+m_nSBoxY)
{
float ss;
int l,k;
int M=0;
N=0;
float s2=0.0f;
int MM=0;
temp1[MM]=0;
temp2[MM]=0;
///////////
for(l=0;l<m_nSBoxX;l++)
{
for(k=0;k<m_nSBoxY;k++)
{
if((Z11[x+l][y+k])>=m_fLowDbz&&(Z22[ll+l][kk+k])>=m_fLowDbz)
{
temp1[MM]=(Z11[x+l][y+k]);
temp2[MM]=(Z22[ll+l][kk+k]);
MM++;
}
if(Z22[ll+l][kk+k]>=m_fLowDbz)
{
s2=s2+(Z22[ll+l][kk+k]);
M=M+1;
}
}
}
if(M>=m_fLowPercengtage*m_nSBoxX*m_nSBoxY)
{
s2=s2/M;
float s3=0.0f;
float s4=0.0f;
float s5=0.0f;
////////////这个判据并没有起作用？/////////////////////////////////////
if(ll>200&&kk>200)
{
int mmm=0;
}
////////////////////////////////////////////////////////////////////////
if(MM<m_fLowPercengtage*m_nSBoxX*m_nSBoxY)
continue;
for(m=0;m<MM;m++)
{
s3=s3+(temp1[m]-s1)*(temp2[m]-s2);
s4=s4+(temp1[m]-s1)*(temp1[m]-s1);
s5=s5+(temp2[m]-s2)*(temp2[m]-s2);
}
if(s4*s5<=0)
ss=0;
else
ss=s3/float(sqrt((s4*s5)));
if(ss>rho) ////////取最大相关系数，并记下该格点的位置
{
rho=ss;
nX0=ll+m_nSBoxX/2.0f;
nY0=kk+m_nSBoxY/2.0f;
// fprintf(file,"%f %f \n",nX0,nY0);
/// int iii=0;
}
}
}
}
fprintf(file,"%f ",rho);
fprintf(file,"\n");
}
delete[]temp1;
delete[]temp2;
}
void CCorrCoff::Intervence(unsigned char Z1[][1000],unsigned char Z2[][1000],unsigned char Z3[][1000], float theta[],FILE *file)
{
// File *file;
if(file==NULL)
return;
float dx=0.25f; ////x方向的格距
float dy=0.25f; ////y方向的格距
int i,j,k;
float xx; ////x方向的距离，以雷达站点为中心
float yy; ////y方向的距离，以雷达站点为中心
float RR; ////空间某点距离雷达的斜距
float zz; ////空间某点的高度，也就是要插值的高度层
float Alpha; ////方位角
float m_theta; ////仰角
float Pi=3.14159f;
float deg=Pi/180;
// int index;
int Ze; /////最后得到的各点上的值
int m_nIndex=-1;
for(j=512;j>=0;j--) /////512*512的空间距离，从该区域的左上角开始记录数据
{
yy=(j-256)*dy; /////每一行上的y值
for(i=0;i<=512;i++)
{
//////将格点上直角坐标的位置转化的球坐标下
zz=3.0f;
xx=(i-256)*dx; /////x的值
RR=float(sqrt(xx*xx+yy*yy+zz*zz));
Alpha=float(atan(yy/xx)); /////方位
if(xx>=0)Alpha=90-Alpha/deg;
else
Alpha=270-Alpha/deg;
zz=zz-float(sqrt(RR*RR))/17000; ////标准大气折射处理
m_theta=float(asin(zz/RR))/deg; ////仰角
int jj=(int)(RR/0.25f); ////订正到格点
int ii=(int)(Alpha);
for(k=0;k<3;k++) ///判断所要插值的在哪一层
{
float mm=theta[0];
if(m_theta<theta[0])m_nIndex=0;
else if(m_theta>theta[k])m_nIndex=k+1;
}
float Ze1=0.0f,Ze2=0.0f;
switch(m_nIndex)
{
case 0:
Ze=int(((Z1[ii][jj+1]*(Alpha-ii)+Z1[ii+1][jj+1]*(ii+1-Alpha))*((jj+1)*0.25f-RR)
+(Z1[ii][jj]*(Alpha-ii)+Z1[ii+1][jj]*((ii+1)-Alpha))*(RR-jj*0.25f))/0.25f);
break;
case 1:
Ze1=((Z1[ii][jj+1]*(Alpha-ii)+Z1[ii+1][jj+1]*((ii+1)-Alpha))*((jj+1)*0.25f-RR)
+(Z1[ii][jj]*(Alpha-ii)+Z1[ii+1][jj]*((ii+1)-Alpha))*(RR-jj*0.25f))/(0.25f);
Ze2=((Z2[ii][jj+1]*(Alpha-ii)+Z2[ii+1][jj+1]*((ii+1)-Alpha))*((jj+1)*0.25f-RR)
+(Z2[ii][jj]*(Alpha-ii)+Z2[ii+1][jj]*((ii+1)-Alpha))*(RR-jj*0.25f))/(0.25f);
Ze=int((Ze1*(m_theta-theta[0])+Ze2*(theta[1]-m_theta))/(theta[1]-theta[0]));
break;
case 2:
Ze1=((Z2[ii][jj+1]*(Alpha-ii)+Z2[ii+1][jj+1]*((ii+1)-Alpha))*((jj+1)*0.25f-RR)
+(Z2[ii][jj]*(Alpha-ii)+Z2[ii+1][jj]*((ii+1)-Alpha))*(RR-jj*0.25f))/(0.25f);
Ze2=((Z3[ii][jj+1]*(Alpha-ii)+Z3[ii+1][jj+1]*((ii+1)-Alpha))*((jj+1)*0.25f-RR)
+(Z3[ii][jj]*(Alpha-ii)+Z3[ii+1][jj]*((ii+1)-Alpha))*(RR-jj*0.25f))/(0.25f);
Ze=int((Ze1*(m_theta-theta[1])+Ze2*(theta[2]-m_theta))/(theta[2]-theta[1]));
case 3:
Ze=int(((Z3[ii][jj+1]*(Alpha-ii)+Z3[ii+1][jj+1]*((ii+1)-Alpha))*((jj+1)*0.25f-RR)
+(Z3[ii][jj]*(Alpha-ii)+Z3[ii+1][jj]*((ii+1)-Alpha))*(RR-jj*0.25f))/(0.25f));
break;
default:
break;
}
fprintf(file,"%2d ",Ze);
}
fprintf(file,"\n");
}
}
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -