bathmath.cpp

数字摄影测量相关算法编程,在VC环境下实现影像特征点提取.

#include "stdafx.h"
#include "basemath.h"
#include <windows.h>
#include <windowsx.h>
#include <math.h>
#include <stdlib.h>

/*************************************************************************

  Gauss_zy(double a[],double b[],int n,double x[])

   高斯主元消去法-答解线性方程

   Returns:  1 if successful.
			 0 if an error occurs.

  Comments:  a[n*n]----------------Coefficient Matricx
  			 b[n*1]----------------Constant Matricx
			 n	   ----------------The Dimension numbers of the Equation 
			 x[n*1]----------------The result of Equation

*************************************************************************/
 BOOL Gauss_zy(double a[],double b[],int n,double x[])
 {
 int l,k,i,j,is,p,q,*js;
 double d,t;

  js=(int *)GlobalAllocPtr(GHND,n*sizeof(int));

  l=1;
  for(k=0;k<=n-2;k++)
  {d=0.0;
   for(i=k;i<=n-1;i++)
    for(j=k;j<=n-1;j++)
      {t=fabs(a[i*n+j]);
       if(t>d) {d=t;js[k]=j;is=i;}
      }
   if(d+1.0==1.0) l=0;
   else
    {if(js[k]!=k)
      for(i=0;i<=n-1;i++)
        {p=i*n+k;q=i*n+js[k];
         t=a[p];a[p]=a[q];a[q]=t;
        }
      if(is!=k)
        {for(j=k;j<=n-1;j++)
           {p=k*n+j;q=is*n+j;
            t=a[p];a[p]=a[q];a[q]=t;
           }
          t=b[k];b[k]=b[is];b[is]=t;
        } 
    }
    if(l==0)
      {
      //free(js);
      GlobalFreePtr(js);
	  return FALSE;
      }
    d=a[k*n+k];
    for(j=k+1;j<=n-1;j++)
      {p=k*n+j;a[p]=a[p]/d;}
    b[k]=b[k]/d;
    for(i=k+1;i<=n-1;i++)
      {for(j=k+1;j<=n-1;j++)
         {p=i*n+j;
	  a[p]=a[p]-a[i*n+k]*a[k*n+j];
         }
         b[i]=b[i]-a[i*n+k]*b[k];
      }
    }
   d=a[(n-1)*n+n-1];
   if(fabs(d)+1.0==1.0)
     {
     //free(js);
     GlobalFreePtr(js);
	 return FALSE;
     }
   x[n-1]=b[n-1]/d;
   for(i=n-2;i>=0;i--)
     {t=0.0;
      for(j=i+1;j<=n-1;j++)
        t=t+a[i*n+j]*x[j];
      x[i]=b[i]-t;
     }
   js[n-1]=n-1;
   for(k=n-1;k>=0;k--)
     if(js[k]!=k)
       {t=x[k];x[k]=x[js[k]];x[js[k]]=t;}
   //free(js);
   GlobalFreePtr(js);
   return TRUE;
 }
float average(BYTE* A,short M,short i0,short j0,short N)
{
	float ave=0.f;
	for(int i=0;i<N;i++)
	{
		for(int j=0;j<N;j++)
		{
			ave+=(float)(A[(i+i0)*M+j+j0])/(N*N);
		}
	}
	return ave;
}
float sigma(BYTE* A,short N,short i0,short j0,float aveA,short M)
{
	float sg=0.f;
	for(int i=0;i<M;i++)
	{
		for(int j=0;j<M;j++)
		{
			sg+=(float)(A[(i+i0)*N+j+j0]*A[(i+i0)*N+j+j0])/(M*M);
		}
	}
	return (sg-aveA*aveA);
	
}

float sigma_AB(BYTE* A,short i0,short j0,float aveA,short N,		 // 协方差
					BYTE* B,float aveB,short M)
{
	float sgAB=0.f;
	for(int i=0;i<M;i++)
	{
		for(int j=0;j<M;j++)
		{
			sgAB+=(float)(A[(i+i0)*N+(j+j0)]*B[i*M+j])/(M*M);
		}
	}
	return (sgAB-aveA*aveB);

}
float  singlead(BYTE* A,BYTE* B,short N1,short N2,short i0,short j0)  //平均绝对差
{
	float sumad=0.0;
	for(int i=0;i<N1;i++)
	{
		for(int j=0;j<N1;j++)
		{
			sumad+=fabs(B[(i+i0)*N2+(j+j0)]-A[i*N1+j]);
		}
	}
    return sumad/(N1*N1);
}
/*
char* int_boi(long xx)//十进制数转化为二进制数
{
   int i=0;
   long sum=0;
   long rest;
   do {	   
	   rest=xx%2;
	   xx=xx/2;	
	   sum+=rest*pow(10,i);	 
	   i++;
   }while (xx!=0);
   return sum;
}*/