ctlib.c

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

/*
*  Copyright 2006 by OCT lab.
*  Author: Lihui.
*  2006.12.13  
*  Version: 1.0 
*/
/*======CTLib.c=======*/


#include "CTLib.h"



/*
//////////////////////////////////////////////////////////////////////////
//CTart()
void CTart(float *p,float **W,float *F,int K,int S,int M,int N,int Nadj,float RelaxFactor)//
{

	float AthwartP=0.0;//逆投影值
	float AdjFactor=0.0;//象素调整因子
	float sumW_ik=0.0;//射线对应投影矩阵某行元素之和 
	float sumW_ik2=0.0;//射线对应投影矩阵某行内积
	float *AdjF=new float[M*N];//象素调整值

    int i_k=0;//用于判断迭代过程中对应射线序列号
    
    
    

   for (int k=0;k<=Nadj*K*S;k++)//迭代过程
   {
       i_k=mod(k,K*S);
	   sumW_ik=VectorSum(W[i_k],M*N);
      if (sumW_ik!=0)
	  {
          AthwartP=p[i_k]-VectorMulA(W[i_k],F,M*N);
		  sumW_ik2=VectorMulA(W[i_k],W[i_k],M*N);
		  AdjFactor=RelaxFactor*AthwartP/sumW_ik2;
          VectorNumMul(W[i_k],AdjFactor,AdjF,M*N,0);
          VectorSum2(F,AdjF,F,M*N,0);
        
	  }//end if sumW_ik
   
      for(int j=0;j<M*N;j++)//虑除小于0的值
	  {
		  if(F[j]<0)
		  {
			  F[j]=0;
		  }//end if F[j]
	  }//end for j

   }//end for k
   delete []AdjF;
}

//end CTart()
////////////////////////////////////////////////////
*/
//////////////////////////////////////////////////////////////////////////
//CTsirt()
//buffers used by SIRT
float Wt_buffer[NP_MAX*RPV_MAX*RECONMN_MAX*RECONMN_MAX];//投影矩阵的转置Sirt
float P1_buffer[NP_MAX*RPV_MAX];//投影估计向量
float F1_buffer[RECONMN_MAX*RECONMN_MAX];//象素调整值
float F2_buffer[RECONMN_MAX*RECONMN_MAX];//象素调整值

void CTsirt(float *P,float *W,float *F,int K,int S,int M,int N,int Nadj,float RelaxFactor)
{
    int k=0;//loop variable
    int j=0;
       
    float *Wt=&Wt_buffer[0];//投影矩阵的转置
    float *P1=&P1_buffer[0];//投影估计向量
    float *F1=&F1_buffer[0];//象素调整值
    float *F2=&F2_buffer[0];//象素调整值
    
    //float *Wt=(float*)malloc(sizeof(float)*K*S*M*N);//投影矩阵的转置
    //float *P1=(float*)malloc(sizeof(float)*K*S);//投影估计向量
    //float *F1=(float*)malloc(sizeof(float)*M*N);//象素调整值
    //float *F2=(float*)malloc(sizeof(float)*M*N);//象素调整值
    
    MatrixTrans(W,Wt,K*S,M*N);
    MVMul(Wt,P,F,M*N,K*S);

   for (k=0;k<=Nadj;k++)//迭代过程
   {
      MVMul(W,F,P1,K*S,M*N);
      MVMul(Wt,P,F1,M*N,K*S);
      MVMul(Wt,P1,F2,M*N,K*S);
      VectorSum2(F1,F2,F1,M*N,1);
      VectorNumMul(F1,RelaxFactor,F1,M*N,0);
      VectorSum2(F,F1,F,M*N,0);

   
      for(j=0;j<M*N;j++)//虑除小于0的值
	  {
		  if(F[j]<0)
		  {
			  F[j]=0;
		  }//end if F[j]
	  }//end for j

   }//end for k
   
   //free(Wt);
   //free(P1);
   //free(F1);
   //free(F2);


}

//end CTsirt()
////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////
//CTMert()两投影方向的最大熵重建算法
void CTMert(float *P,float *F,int N)
{
    int i=0;
    int j=0;
	for(i=0;i<N;i++)
	{
		for(j=0;j<N;j++)
		{
			F[i*N+(N-j-1)]=4*P[i]*P[N-1+j];
		}
	}
}
//end CTMert()
/////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////
//matrix operation 

////////////////////////////////////////////////////////////////////////////
//VectorDotMul()向量对应元素相乘
void VectorDotMul(float *m1,float *m2,float *mo,int Nc)
{
    int n=0;
	for(n=0;n<Nc;n++)
	{
		mo[n]=m1[n]*m2[n];
	}
}

//////////////////////////////////////////////////////////////////////////
//VectorSum()向量元素求和
float VectorSum(float *m,int Nc)
{
    int n=0;
	float sum=0.0;
	for(n=0;n<Nc;n++)
	{
		sum+=m[n];
	}
	return sum;
}
///////////////////////////////////////////////////////////////////////////
//VectorMulA()向量内积
float VectorMulA(float *m1,float *m2,int Nc)
{
    int n=0;
	float mo=0.0;
	for(n=0;n<Nc;n++)
	{
		mo+=m1[n]*m2[n];
	}
	return mo;
}

////////////////////////////////////////////////////////////////////////
//VectorMulB()向量叉乘
void VectorMulB(float *m1,float *m2,float *mo,int Nc)
{
    int i=0;
    int j=0;
    
	for(i=0;i<Nc;i++)
	{
		for(j=0;j<Nc;j++)
		{
			mo[i*Nc+j]=0.0;
		}
	}

	for(i=0;i<Nc;i++)
	{
		for(j=0;j<Nc;j++)
		{
			mo[i*Nc+j]=m1[i]*m2[j];
		}
	}
}

/////////////////////////////////////////////////////////////////////////
//MatrixMul() M×N矩阵和N×M矩阵相乘得M×M矩阵
void MatrixMul(float *m1,float *m2,float *mo,int Nr,int Nc)
{
    int i=0;
    int j=0;
    int m=0;
    int n=0;
    
	for(i=0;i<Nr;i++)
	{
		for(j=0;j<Nr;j++)
		{
			mo[i*Nr+j]=0.0;
		}
	}

	for(m=0;m<Nr;m++)
	{
		for(n=0;n<Nr;n++)
		{
			for(i=0;i<Nc;i++)
			{
			   mo[m*Nr+n]+=m1[m*Nc+i]*m2[i*Nr+n];
			}//end for i
		}//end for n
	}//end for m
}
////////////////////////////////////////////////////////////////////

/////////////////////////////////////////////////////////////////////////
//MatrixMul() M×N矩阵和N维向量矩阵相乘得M维向量
void MVMul(float *m1,float *m2,float *mo,int Nr,int Nc)
{
    int i=0;
    int m=0;
    
	for(i=0;i<Nr;i++)
	{
		
			mo[i]=0.0;
		
	}

	for(m=0;m<Nr;m++)
	{
		
	
			for(i=0;i<Nc;i++)
			{
			   mo[m]+=m1[m*Nc+i]*m2[i];
			}//end for i
		
	}//end for m
}
////////////////////////////////////////////////////////////////////
//MatrixTrans()矩阵转置
void MatrixTrans(float *m1,float *m2,int M,int N)
{
    int m=0;
    int n=0;
    
	for(m=0;m<M;m++)
	{
		for(n=0;n<N;n++)
		{
			m2[n*M+m]=m1[m*N+n];
		}
	}
}



//////////////////////////////////////////////////////////////////
//VectorSum2()向量对应元素相加减
void VectorSum2(float *m1,float *m2,float *mo,int Nc,int mod)
{
    int n=0;
    
	for(n=0;n<Nc;n++)
	{
		switch(mod)
		{
		case 0:
			mo[n]=m1[n]+m2[n];
			break;
		case 1:
			mo[n]=m1[n]-m2[n];
			break;
		default:
			break;
		}
	}
	
}

////////////////////////////////////////////////////////////////////////
//VectorNumMul()向量对应元素与常数相乘除
void VectorNumMul(float *m1,float num,float *mo,int Nc,int mod)
{
    int n=0;
    
	for(n=0;n<Nc;n++)
	{
		switch(mod)
		{
		case 0:
			mo[n]=m1[n]*num;
			break;
		case 1:
			mo[n]=m1[n]/num;
			break;
		default:
			break;
		}
	}
}

///////////////////////////////////////////////////////////////////////
//mod()  used in CTart() only
int mod(int n1,int n2)
{
	while(n1>=n2)
	{
		n1-=n2;
	}
	return n1;
}

/////////////////////////////////////////////////////////////////////////
//maxv()求向量中最大元素
float maxv(float *Vector,int N)
{
    int n=0;
	float temp=0.0;
	
	for(n=0;n<N;n++)
	{
		if (Vector[n]>temp)
		{
			temp=Vector[n];

		}
	}
	return temp;
}

///////////////////////
/////////////////////////////////////////////////////////////////////////
//maxv2()求二维向量中最大元素
float maxv2(float *Matrix,int M,int N)
{   
    int m=0;
    int n=0;
	float temp=0.0;
	for(m=0;m<M;m++)
		for(n=0;n<N;n++)
	{