ctlib.c

本源码经过上机调试,是CT算法在TI的CCS下编程 可以在DSP硬件和软件仿真条件下运行,同时对CT算法在ARM,MIPS,PC,FPGA等上实现都有借鉴意义.搞CT等重建算法的人值得一看

		if (Matrix[m*N+n]>temp)
		{
			temp=Matrix[m*N+n];
		}
	}
	return temp;
}
///////////////////////
/////////////////////////////////////////////////////////////////////////
//minv()求向量中最小元素
float minv(float *Vector,int N)
{
    int n=0;
	float temp=Vector[0];
	
	for(n=0;n<N;n++)
	{
		if (Vector[n]<temp)
		{
			temp=Vector[n];
		}
	}
	return temp;
}

///////////////////////
//absv()向量元素取绝对值
void absv(float *Vector,int N)
{
    int n=0;
    
	for (n=0;n<N;n++)
	{
		Vector[n]=fabs(Vector[n]);
	}
}

//end matrix operation 
////////////////////////////////////////////////////




/////////////////////////////////////////////////////////////////////////////
//投影矩阵CTw()
//////////////////////////////////////////////////////////////

//////////////////////////////////////////////////////////////
//ctw()
void CTw(float *W,int K,int S,int M,int N)//
{
 
    int J=M*N;
	int i=0;
	int j=0;
	int k=0;
	int s=0;
	int n=0;
	int m=0;
	
	float lx=1.0/(M-1.0),ly=1.0/(N-1.0),d=1.0/(S-1.0),a=PI/K;



    for (k=0;k<K;k++)
	{
		for (s=0;s<S;s++)
		{
            float t=-1.0/2.0*(1.0-sqrt(2.0)*sin(-k*a+PI/4.0))+s*d;


			for (n=0;n<N;n++)
			{
				for(m=0;m<M;m++)
				{
					i=k*S+s;
					j=n*M+m;
					if ((fabs(tan(k*a))>=0.0) && (fabs(tan(k*a))<=1.0))//((fabs(tan(k*a))>=0.0) & (fabs(tan(k*a))<=1.0))
					{
					W[i*J+j]=lx*fabs(sec(k*a))*sinc((t*sec(k*a)+m*lx*tan(k*a)-n*ly)/ly);
					}
					else if((fabs(tan(k*a))>1.0) && ((k*a)<((PI/2)-1e-6)||(k*a)>((PI/2)+1e-6)))//((fabs(tan(k*a))>1.0) & ((k*a)!=(PI/2))) 
					{
						W[i*J+j]=ly*fabs(sec(k*a))/fabs(tan(k*a))*sinc((t*sec(k*a)+m*lx*tan(k*a)-n*ly)/(lx*tan(k*a)));
					}
                    else if ((k*a)>((PI/2)-1e-6)&&(k*a)<((PI/2)+1e-6)) //((k*a)==(PI/2))
					{
						W[i*J+j]=ly*sinc((t+m*lx)/lx);
					}
				}//end for m
			}//end for n
		}//end for s
	}//end for k	
	   		   
}


/////////////////////////////////////////
//sinc()
float sinc(float x)
{
	return sin(x*PI)/(x*PI);
}

///////////////////////////////////////////////
//sec()
float sec(float x)
{

	return 1.0/cos(x);
}

//end 投影矩阵CTw()
/////////////////////////////////////////////////////////////////////////////


/*

////////////////////////////////////////////////////////////////////////
//测试函数Agauss(x,y)四峰非对称高斯函数，可模拟二维四峰非对称温度场分布good
//在X和Y的绝对值小于0.5的范围之内
float Agauss(float x,float y)
{
	 float xi[]={0.15, 0.15, -0.15 ,-0.15};//%四峰的ｘ坐标
	 float yi[]={0.15,-0.15, 0.15 ,-0.15};//%四峰的y坐标
	 float ai[]={1, 0.4 ,0.6, 0.8};//%四峰的幅度
     float ag=0.0;
     if(x>=-0.5&&x<=0.5&&y>=-0.5&&y<=0.5)
	 {
       for(int n=0;n<4;n++)
	   {
          ag=ag+ai[n]*exp(-4.0*log(2)/pow(0.2,2.0)*pow((x-xi[n]),2.0)
			  -4*log(2.0)/pow(0.2,2.0)*pow((y-yi[n]),2.0));
	   }
	 }
    else
    ag=0;
	return ag;
}//end Agauss

////////////////////////////////////////////////////////////////////////

//////////////////////////////////////////////////////////////////////////
//origF()
void origF(float *Fo,int M,int N)
{
	float x0=-0.5;
	float y0=-0.5;
	float d=1.0/M;
	for (int n=0;n<N;n++)
	{
		for (int m=0;m<M;m++)
		{
			Fo[n*M+m]=Agauss(x0+d/2.0+m*d,y0+d/2.0+n*d);
           
		}//取象素中点的值作为该象素值
	}
}
//end orgF()
////////////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////
//CTsimpson()求投影值
//p--投影值，K--方向数，S--每方向投影数
void CTsimpson(float *p,int K,int S)//
{
    //da间隔角度，xra xrb  x轴起止点，yra yrb  y轴起止点，y轴间隔，h积分步长，p投影向量
	//int K=8;
	//int S=26;
	int L=51;
	float da=PI/K;
	float xra=-sqrt(2)/2.0;
	float xrb=sqrt(2)/2.0;
	float yra=-0.5;
	float yrb=0.5;
	float d=(yrb-yra)/(S-1);
	float h=(xrb-xra)/(L-1)/2;

	float a=0.0;
	float yr=0.0;
	float v[2][2];
	float xr=0.0;
	float f1=0.0;
	float f2=0.0;
	float xy[2];
	float fa=0.0;
	float fb=0.0;

  for(int k=0;k<=K-1;k++)
  {
    for(int s=0;s<=S-1;s++)
	{
        a=k*da;  //%确定射线方向
        yr=yra+s*d;  //%确定射线沿着yr轴的位置

        v[0][0]=cos(a);  // %旋转向量使积分线以da大小旋转
        v[0][1]=-sin(a);
		v[1][0]=sin(a);
		v[1][1]=cos(a);

        f1=0.0;
	    f2=0.0;

        for(int  l=0;l<=L-1;l++)
		{
            xr=xra+(2*l+1)*h;
            
			xy[0]=v[0][0]*xr+v[0][1]*yr;
            xy[1]=v[1][0]*xr+v[1][1]*yr;

            f1=f1+Agauss(xy[0],xy[1]);
        }//end for l1

        for(l=1;l<=L-1;l++)
		{
            xr=xra+2*l*h;

			xy[0]=v[0][0]*xr+v[0][1]*yr;
            xy[1]=v[1][0]*xr+v[1][1]*yr;

            f2=f2+Agauss(xy[0],xy[1]);
        }//end for l2

        xy[0]=v[0][0]*xra+v[0][1]*yr;
        xy[1]=v[1][0]*xra+v[1][1]*yr;
        fa=Agauss(xy[0],xy[1]);

        xy[0]=v[0][0]*xrb+v[0][1]*yr;
        xy[1]=v[1][0]*xrb+v[1][1]*yr;
        fb=Agauss(xy[0],xy[1]);

        p[k*S+s]=h/3.0*(fa+4.0*f1+2.0*f2+fb);
    }//end for s
  }//end for k
}//end simpson

//end CTsimpson()求投影值
//////////////////////////////////////////////////////////////////////////

////////////////////////////////////////////////////////////////////////////
//error analysis

///////////////////////
//Eaver()平均相对误差
float Eaver(float *Fo,float *Fr,int J)
{
	float error=0.0;
	float *Ftemp=new float[J];
    VectorSum2(Fo,Fr,Ftemp,J,1);
	absv(Ftemp,J);
	error=VectorSum(Ftemp,J)/(maxv(Fo,J)*J);
	delete []Ftemp;
	return error;
}

///////////////////////
//Emax()最大相对误差
float Emax(float *Fo,float *Fr,int J)
{
	float error=0.0;
	float *Ftemp=new float[J];
    VectorSum2(Fo,Fr,Ftemp,J,1);
	absv(Ftemp,J);
	error=maxv(Ftemp,J)/maxv(Fo,J);
	delete []Ftemp;
	return error;
}

///////////////////////
//Esqrt()均方根误差
float Esqrt(float *Fo,float *Fr,int J)
{
	float error=0.0;
	float *Ftemp=new float[J];
    VectorSum2(Fo,Fr,Ftemp,J,1);
	error=sqrt(VectorMulA(Ftemp,Ftemp,J)/VectorMulA(Fo,Fo,J));
	delete []Ftemp;
	return error;
}

//end error analysis
//////////////////////////////////////////////////////////////////

////////////////////////////////////
//单色测温TOneColor()
void TOneColor(float **F,float Wavelength,double Parameter,int H,int J)
{
   double c1=3.7418e-16;
   double c2=1.4388e-2;
   double wavelength=Wavelength*1e-9;

   for(int h=0;h<H;h++)
   {
	for (int j=0;j<J;j++)
	{
		if( F[h][j]<=0)
		{
			F[h][j]=300;
		}
		else
		{
            F[h][j]=c2/wavelength/(logf(c1/PI/(wavelength*wavelength*wavelength*wavelength*wavelength))-logf(F[h][j])+Parameter);
		}
	}
   }
}
//end TOneColor()
/////////////////////////////////////////////////////////////

/////////////////////////////////////////////
//双波长测温TTwoColor()
void TTwoColor(float **F,float **FWav1,float **FWav2,double Parameter2,int H,int J)
{
   double c1=3.7418e-16;
   double c2=1.4388e-2;
   double WaveLength1=670e-9;
   double WaveLength2=780e-9;
   double tempF2=0.0;
   float Fsum=0.0;

   for(int h=0;h<H;h++)
   {
    tempF2=maxv(FWav2[h],J);
   Fsum=VectorSum(FWav2[h],J);
   Fsum/=J;
	for (int j=0;j<J;j++)
	{
//		if(FWav1[h][j]<0.1||FWav2[h][j]<0.2*tempF2|| FWav1[h][j]* FWav2[h][j]<0.05*tempF2*tempF2)  //    FWav1[h][j]*FWav2[h][j]<=8000
		if( FWav1[h][j]<0.1||FWav2[h][j]<1.7*Fsum|| FWav1[h][j]* FWav2[h][j]<1.7*Fsum*Fsum)

		{
			F[h][j]=300;
		}
		else
			if( FWav1[h][j]/FWav2[h][j]>2.5)
			{
				F[h][j]=2.5;
			}
			else
		{
            F[h][j]=c2*(1/WaveLength2-1/WaveLength1)/(logf(FWav1[h][j]/FWav2[h][j])-5*logf(WaveLength2/WaveLength1)+Parameter2);
//F[h][j]=FWav1[h][j]/FWav2[h][j];
		}

//F[h][j]=FWav2[h][j];
	}
   }
}
//end TTwoColor()
//////////////////////////////////////////////////////////

*/