dibapi.cpp

影像融合与融合精度评价源码影像融合与融合精度评价源码

						{
							g = v[j*n+ip];
							h = v[j*n+iq];
							v[j*n+ip] = g - s * (h + g * tau);
							v[j*n+iq] = h + s * (g - h * tau);
						}
					}
				}
			}
        }
        for (ip = 1; ip<n; ip++)
		{
            b[ip] = b[ip] + z[ip];
            d[ip] = b[ip];
            z[ip] = 0.0;
        }
    }
}
void  WINAPI   eigsrt(double * d, double* v, int n,double * eva)
{
	double p;
	int i,j,k;
    for (i = 1; i<n - 1; i++)
	{
        k = i;
        p = d[i];
        for (j = i + 1; j<n; j++)
		{
            if (d[j] >= p)
			{
				k = j;
                p = d[j];
			}
        }
        if (k != i )
		{
            d[k] = d[i];
            d[i] = p;
			p=eva[i];
			eva[i]=eva[k];
			eva[k]=p;
            for (j = 1; j<n; j++)
			{
                p = v[j*n+i];
                v[j*n+i] = v[j*n+k];
                v[j*n+k] = p;
            }
		}
    }
}

BOOL   WINAPI  PCA(double * v,double * Means, int n,int *multi_image,int width,int height)
{
	int ww,hh;
	ww=width*(n-1);
//	ww=width;
	hh=height;
	double *pa;
	BYTE* inter;
	inter=new BYTE [ww*hh];
	pa=new double [width*height];
	if(pa==NULL || inter==NULL)
	{
		AfxMessageBox("Memory no enough!");
		return FALSE;
	}

	long size;
	size=width*height;
	int l,i,j;
    for(l=1;l<n;l++)
	{
		for(i=0;i<size;i++)
		{
			pa[i]=multi_image[(l-1)*size+i]-Means[l];
		}
		for(i=0;i<size;i++)
		{
		   double p;
		   p=0.0;
           for(int k=1;k<n;k++)       //n 可以控制需要的特征向量数
		   {
			   p+=v[k*n+l]*multi_image[(k-1)*size+i];
		   }
		   pa[i]=p;
		}
		
		Norm(pa,size);
		for(i=0;i<height;i++)
			for(j=0;j<width;j++)
			{
				inter[i*ww+j+(l-1)*width]=(BYTE)pa[i*width+j];
			}
	}
    FILE  * stream;
	if( (stream = fopen( "主成分变换数据.img", "wb" )) != NULL)
	{
		for(int k=1;k<n;k++)
		{
		  for(i=0;i<height;i++)
		  for(j=0;j<width;j++)
		  {
			  int temp;
			  temp=inter[(k-1)*width+i*ww+j];
		      fwrite(&temp, 2, 1, stream );
		  }
		}
		fclose( stream );
	}

// 产生新的影像
		LPSTR                   pDib;
		DWORD                   dwTargetBitsSize;
		LPBITMAPINFO            lpbmi=NULL;

		int L_W;
		L_W = (int)ww % sizeof(DWORD);
		if(L_W != 0)
			L_W = sizeof(DWORD) - L_W;

		HDIB hNewDIB=CreateDIB(ww,hh,24);
		if (hNewDIB == 0)	return true;

		pDib = (LPSTR) ::GlobalLock((HGLOBAL) hNewDIB);	
		lpbmi= (LPBITMAPINFO)pDib;
		dwTargetBitsSize=lpbmi->bmiHeader.biHeight*BytesPerLine((LPBYTE)&(lpbmi->bmiHeader));
		lpbmi->bmiHeader.biSizeImage=dwTargetBitsSize;
		DWORD  wBytesPerLine=BytesPerLine((LPBYTE)&(lpbmi->bmiHeader));

		BYTE * image1;
		image1=new BYTE [wBytesPerLine*hh];

		for(i=0;i<hh;i++)
		{
			 for(j=0;j<ww;j++)
			 {
				BYTE im;
				im=inter[i*ww+j];
				LONG lOffset=PIXEL_OFFSET(i,j,wBytesPerLine);
					*(image1+lOffset++)=im;
					*(image1+lOffset++)=im;
					*(image1+lOffset)=im;
			 }
		 }
		memcpy((BYTE*)FindDIBBits(pDib),image1,wBytesPerLine*hh);
		delete [] image1;

		/*影像存储*/
		HDIB hNewDIB1;
		hNewDIB1=ConvertDIBFormat(hNewDIB,8,NULL);
		char Buffer[80];
		wsprintf(Buffer,"主成分变换影像.bmp");
		CFile RFile;
		RFile.Open(Buffer, CFile::modeCreate | CFile::modeWrite);
 		::SaveDIB(hNewDIB1,RFile);
		RFile.Close();
		::GlobalUnlock((HGLOBAL) hNewDIB);
	
	delete [] pa;
	delete [] inter;
	return TRUE;
}
void   WINAPI    Norm(double * buffer, long size)
{
	long i;
	int max=0;
	int min=0;
	for(i=0;i<size;i++)
	{
		if(max<buffer[i]) max=(int)buffer[i];
		if(min>buffer[i]) min=(int)buffer[i];
	}
	double scale=0.0;
	scale=(double)(255-0)/(max-min);
	for(i=0;i<size;i++)
		buffer[i]=(int)(((int)buffer[i]-min)*scale);
	return;
}

void  WINAPI  Means(int * buffer,long size,int sp,double * means)
{
	long i,k;
	for(k=0;k<sp;k++)
	{
		means[k]=0;
		for(i=0;i<size;i++)
			means[k]+=buffer[k*size+i];
		means[k]/=size;
	}
}

void  WINAPI Variance(int * buffer,double * means,long size,int sp,double * variance)
{
	int i,j;
	for(i=0;i<sp;i++)
	{
		for(j=i;j<sp;j++)
		{
			int m;
			double var;
			var=0.0;
			for(m=0;m<size;m++)
			{
				var+=(buffer[i*size+m]-means[i])*(buffer[j*size+m]-means[j]);
			}
			variance[i*sp+j]=variance[j*sp+i]=var/(size-1);
		}
	}
}

void  WINAPI Correlation(int * buffer,long size,int sp,double * cof)
{
	int j,i1,i2;
	for(j=0;j<size;j++)
   	{
    	for(i1=0;i1<sp;i1++)
     	for(i2=0;i2<sp;i2++)
      	cof[i1*sp+i2]+=1.0*buffer[i1*size+j]*buffer[i2*size+j]/size;
   	}
}
//求行列式
double WINAPI  bsdet(double * a,int n)
{
	int i,j,k,is,js,l,u,v;
    double f,det,q,d;
    f=1.0; 
	det=1.0;
    for (k=0; k<=n-2; k++)
    {
		q=0.0;
        for (i=k; i<=n-1; i++)
        for (j=k; j<=n-1; j++)
          { l=i*n+j; d=fabs(a[l]);
	    if (d>q) { q=d; is=i; js=j;}
          }
        if (q+1.0==1.0)
          { det=0.0; return(det);}
        if (is!=k)
          { f=-f;
            for (j=k; j<=n-1; j++)
              { u=k*n+j; v=is*n+j;
                d=a[u]; a[u]=a[v]; a[v]=d;
              }
          }
        if (js!=k)
          { f=-f;
            for (i=k; i<=n-1; i++)
              { u=i*n+js; v=i*n+k;
                d=a[u]; a[u]=a[v]; a[v]=d;
              }
          }
        l=k*n+k;
        det=det*a[l];
        for (i=k+1; i<=n-1; i++)
          { d=a[i*n+k]/a[l];
            for (j=k+1; j<=n-1; j++)
              { u=i*n+j;
                a[u]=a[u]-d*a[k*n+j];
              }
          }
      }
    det=f*det*a[n*n-1];
    return(det);
}

int  WINAPI  brinv(double * a,int n)
  { 
	int *is,*js,i,j,k,l,u,v;
    double d,p;
    //is=malloc(n*sizeof(int));
	is=new int[n];
	js=new int[n];
	if(is == NULL || js == NULL)
	{
		AfxMessageBox("开辟内存出错，退出！",MB_OK | MB_ICONHAND);
		return 0;
	}
    for (k=0; k<=n-1; k++)
      { d=0.0;
        for (i=k; i<=n-1; i++)
        for (j=k; j<=n-1; j++)
          { l=i*n+j; 
		    p=fabs(a[l]);
            if (p>d) 
			{ d=p; is[k]=i; js[k]=j;}
          }
        if (d+1.0==1.0)
          { 
			delete []is;
			delete []js;
			AfxMessageBox("err**not inv\n");
            return(0);
          }
        if (is[k]!=k)
          for (j=0; j<=n-1; j++)
            { u=k*n+j; v=is[k]*n+j;
              p=a[u]; a[u]=a[v]; a[v]=p;
            }
        if (js[k]!=k)
          for (i=0; i<=n-1; i++)
            { u=i*n+k; v=i*n+js[k];
              p=a[u]; a[u]=a[v]; a[v]=p;
            }
        l=k*n+k;
        a[l]=1.0/a[l];
        for (j=0; j<=n-1; j++)
          if (j!=k)
            { u=k*n+j; a[u]=a[u]*a[l];}
        for (i=0; i<=n-1; i++)
          if (i!=k)
            for (j=0; j<=n-1; j++)
              if (j!=k)
                { u=i*n+j;
                  a[u]=a[u]-a[i*n+k]*a[k*n+j];
                }
        for (i=0; i<=n-1; i++)
          if (i!=k)
            { u=i*n+k; a[u]=-a[u]*a[l];}
      }
    for (k=n-1; k>=0; k--)
      { if (js[k]!=k)
          for (j=0; j<=n-1; j++)
            { u=k*n+j; v=js[k]*n+j;
              p=a[u]; a[u]=a[v]; a[v]=p;
            }
        if (is[k]!=k)
          for (i=0; i<=n-1; i++)
            { u=i*n+k; v=i*n+is[k];
              p=a[u]; a[u]=a[v]; a[v]=p;
            }
      }
	delete []is;
	delete []js;
    return(1);
}

void   WINAPI   Matricx_Mul(double a[],double b[],int m,int n,int k,double c[])
 {
 	int i,j,l,u;
	for(i=0;i<=m-1;i++)
	for(j=0;j<=k-1;j++)
		{
		u=i*k+j;
		c[u]=0.0;
		for(l=0;l<=n-1;l++)
			c[u]=c[u]+a[i*n+l]*b[l*k+j];
		}
		return;
}
void   WINAPI   Matricx_tMul(double a[],double b[],int m,int n,int k,double c[])
{
	int i,j,l,u;
	for(i=0;i<n;i++){
		for(j=0;j<k;j++){
			u=i*k+j;c[u]=0.0;
			for(l=0;l<m;l++)
				c[u]=c[u]+a[i+l*n]*b[l*k+j];
			}
		}
		return;

}

void  WINAPI  cal_zz(double * zz,int kindsnum,int sp_no,int * images_buffer,int m_width,int m_height)
{
	if(kindsnum==0)
	{
		AfxMessageBox("Please input the number of the classified type ");
		return;
	}
	
	int j,k;
	//calculate the histogram
/*	for(k=0;k<sp_no;k++)
		for(i=0;i<m_height;i++)
		for(j=0;j<m_width;j++)
		{
			int temp;
			temp=images_buffer[i*m_width+j];
			hism[k*256+temp]++;
		}*/
	static double seg1=0.3333;
	static double seg2=0.6666;
	for(j=0;j<kindsnum;j++)
		for(k=0;k<sp_no;k++)
		zz[j*sp_no+k]=5+j*(250.0/kindsnum);
}
//求逆
void   WINAPI   bssgj(double *a,int n)
{ 
	int i,j,k,m;
    double w,g,*b;
    b=new double [n];
    for (k=0; k<=n-1; k++)
      { w=a[0];
        if (fabs(w)+1.0==1.0)
          { delete [] b; printf("fail\n"); return;}
        m=n-k-1;
        for (i=1; i<=n-1; i++)
          { g=a[i*n]; b[i]=g/w;
            if (i<=m) b[i]=-b[i];
            for (j=1; j<=i; j++)
              a[(i-1)*n+j-1]=a[i*n+j]+g*b[j];
          }
        a[n*n-1]=1.0/w;
        for (i=1; i<=n-1; i++)
          a[(n-1)*n+i-1]=b[i];
      }
    for (i=0; i<=n-2; i++)
    for (j=i+1; j<=n-1; j++)
      a[i*n+j]=a[j*n+i];
    delete [] b;
}
COLORREF   WINAPI   ClassColorIndex(int n)
{
	COLORREF  col;
	switch(n)
	{
		case 1: //红
			{
				col=RGB(255,0,0);
				break;
			}
		case 2: //蓝
			{
				col=RGB(0,0,255);
				break;
			}
		case 3: //黄
			{
				col=RGB(255,255,0);
				break;
			}
		case 4: //紫
			{
				col=RGB(255,0,255);
				break;
			}
		case 5: //青
			{
				col=RGB(0,255,255);
				break;
			}
		case 6:  //绿
			{
				col=RGB(0,255,0);
				break;
			}
		case 7:  //橙
			{
				col=RGB(252,125,1);
				break;
			}
		case 8:  //棕
			{
				col=RGB(117,41,1);
				break;
			}
	    case 9:  //粉
			{
				col=RGB(206,8,248);
				break;
			}
	    case 10:  //褐黄
			{
				col=RGB(219,174,0);
				break;
			}
	    case 11:  //粉
			{
				col=RGB(0,219,185);
				break;
			}
		 case 12:  //青
			{
				col=RGB(52,205,189);
				break;
			}
		 case 13:  //绿黄
			{
				col=RGB(182,233,27);
				break;
			}
		 case 14:  //粉
			{
				col=RGB(0,219,185);
				break;
			}
		 case 15:  //橙黄
			{
				col=RGB(233,197,27);
				break;
			}
        default :  //黑
			{
				col=RGB(0,0,0);
				
			}
	}
	return col;
}
double WINAPI  Distence(double X1[],double X2[],int n)
{
  int i;
  double s;
  for(i=0;i<n;i++)
  {
	 s+=pow((X1[i]-X2[i]),2);
  }
  s=sqrt(s);
  return s;
}
void   WINAPI  BWClustering(int * image,int *Y,int n,int width,int height)
{
    double k=0.5;
	CArray<double,double>  Z;
	int i,j,l,m;
	int * image1;
	image1=new int [n*width*height];
	for(i=0;i<width*height;i++)
		for(l=0;l<n;l++)
		{
			image1[i*n+l]=image[i*n+l];
		}
	for(i=0;i<n;i++)
	{
		Z.Add(image[i]);
		image[i]=-5;
	}
	double dd=0,d1,d2;
	double * X1,*X2;
	X1=new double [n];
	X2=new double [n];
	for(i=1;i<width*height;i++)
	{
		for(l=0;l<n;l++)X1[l]=Z[l];
		for(l=0;l<n;l++)X2[l]=image[i*n+l];
		d1=Distence(X1,X2,n);