imageprocess.cpp

医学图象处理系统

#include "stdafx.h"
#include "ImageLAB.h"
#include "AllImage.h"
#include "ImageProcess.h"
#include "CommonProc.h"
#include "MainFrm.h"
//image process function start from here 

void Image_LinearScalar(BYTE *InputImage, int ImageSize, int a, int b)
{
	int max_r = InputImage[0]; // 图像灰度的最大值
	int min_r = InputImage[0]; // 图像灰度的最小值
	int  iTmpGray;

	for (int i=1; i<ImageSize; i++) // 寻找灰度的最大值及最小值
	{
	    if ( InputImage[i] > max_r )	   max_r = InputImage[i];
		else if ( InputImage[i] < min_r )  min_r = InputImage[i];
	}
	// -----参数的获得-----
	max_r -= min_r;
	if( max_r == 0) return;
	if( a > b ) 
	{
		iTmpGray = a;
		a = b;
		b = iTmpGray;
	}
	double  ratio = (double)(b- a)/ (double)max_r;
	
	// -----建立查找表-----
	BYTE LUT[256];
	
	for(i= min_r; i<256; i++)
	{
		iTmpGray = int( (i- min_r)* ratio + a + 0.5 );
		if( iTmpGray < 0 )
			LUT[i] = 0;
		else if( iTmpGray < 255 )
			LUT[i] = BYTE( iTmpGray );
		else
			LUT[i] = 255;			
	}
	// -----进行线性变换-----
	for(i= 0; i<ImageSize; i++)
	{
		InputImage[i] = LUT[ InputImage[i] ]; 
	}
}

void Image_AutoLevels(BYTE *InputImage, int ImageSize, int a, int b)
{
	int hist[256], iTmpGray;
	BYTE *temp, *end = InputImage + ImageSize;

	memset(hist, 0, 256*sizeof(int) );
	for(temp=InputImage; temp < end; temp++)
	{		
		hist[*temp]++;
	}

	int i, max_r, min_r, Sum, Threshold = int( ImageSize * 0.03 + 0.5 );

	for (Sum= 0, i=0; i<256; i++) 
	{
	    Sum += hist[i];
		if( Sum > Threshold ) 
		{
			min_r = i;
			break;
		}		
	}

	for (Sum= 0, i=255; i>=0; i--) 
	{
	    Sum += hist[i];
		if( Sum > Threshold ) 
		{
			max_r = i;
			break;
		}		
	}
	// -----参数的获得-----
	max_r -= min_r;
	if( max_r == 0) return;
	if( a > b ) 
	{
		iTmpGray = a;
		a = b;
		b = iTmpGray;
	}
	double  ratio = (double)(b- a)/ (double)max_r;
	
	// -----建立查找表-----
	BYTE LUT[256];
	
	for(i= 0; i<256; i++)
	{
		iTmpGray = int( (i- min_r)* ratio + a + 0.5 );
		if( iTmpGray < 0 )
			LUT[i] = 0;
		else if( iTmpGray < 255 )
			LUT[i] = BYTE( iTmpGray );
		else
			LUT[i] = 255;			
	}
	// -----进行线性变换-----
	for(temp=InputImage; temp < end; temp++)
	{		
		*temp = LUT[ *temp ];
	}	
}

void Image_Histogramenhance(BYTE *InputImage,  int ImageSize)
{
	int hist[256];
	BYTE *temp, *end = InputImage + ImageSize;
	memset(hist, 0, 256*sizeof(int) );

	for(temp=InputImage; temp < end; temp++)
	{		
		hist[*temp]++;
	}

	for(int k=1; k < 256; k++)
		hist[k] += hist[k-1];

	for(k=0; k < 256; k++)
		hist[k] = long(hist[k]/(double)ImageSize*255.0+0.5);

	/*----	equalization  ----*/
	for(temp=InputImage; temp < end; temp++)
	{		
		 *temp = BYTE(hist[*temp]);
	}
}
void Image_AutoBinary(BYTE* InputImage , int ImgSize)
{
	if(InputImage==NULL) return;
	
	BYTE *end = InputImage + ImgSize;

	BYTE Threshold, MeanGray, MinGray = InputImage[0];
	double TotalGray = 0.0;
	for(BYTE *temp = InputImage; temp<end; temp++)
	{
		TotalGray += *temp;
		if(*temp < MinGray) MinGray = *temp;
	}
	MeanGray = BYTE(TotalGray/ ImgSize+0.5);
	Threshold = (MeanGray + MinGray)/2;

	TRACE("Threshold = %d ((mean + min)/ 2 method)  \n", Threshold);
	for(temp = InputImage; temp<end; temp++)
	{
		*temp = (*temp>Threshold) ? 255 : 0;
	}
}
//***********************************************************
//函数类别： 二值图像的处理
//函数名称：                       
//           Image_AutoBinary_MaxExpection
//函数用途：
//           对图像进行二值化
//           利用<图像图形学报>99.6
//           "一种基于灰度期望的图像二值化算法" 一文编制程序
//原始作者： 陆  宏  伟
//原始日期： 28/8 /1999
//修改日期:  
//***********************************************************
void Image_AutoBinary_MaxExpection(BYTE* InputImage , // 数据的指针
					   int ImgSize)       // 数据的大小
{
	if(InputImage==NULL) return;
	
	int hist[256];
	BYTE Threshold, *temp, *end = InputImage + ImgSize;
	memset(hist,0,256*sizeof(int));
	
	for(temp=InputImage; temp < end; temp++)
	{		
		hist[*temp]++;
	}

	double TotalGray = 0.0;
	for(int i=0; i<256; i++)
	{
		TotalGray += i * hist[i];
	}
	Threshold = BYTE( TotalGray / ImgSize + 0.5 );

	TRACE("Threshold = %d (MaxExpection method)  \n", Threshold);
	for(temp = InputImage; temp<end; temp++)
	{
		*temp = (*temp>Threshold) ? 255 : 0;
	}
}
//***********************************************************
//函数类别： 二值图像的处理
//函数名称：                       
//           Image_AutoBinary_Qtsu
//函数用途：
//           对图像进行二值化
//           利用<图像图形学报>99.6
//           改进Qtsu算法程序
//原始作者： 陆  宏  伟
//原始日期： 28/08/1999
//修改日期:  11/12/1999
//***********************************************************
void Image_AutoBinary_Qtsu(BYTE* InputImage , // 数据的指针
						   int ImgSize)       // 数据的大小
{
	if(InputImage==NULL) return;
	
	int    i, hist[256], IntegralHist[256], IntegralHistL[256];
	double Sigma[256];
	BYTE   Threshold, *temp, *end = InputImage + ImgSize;
	memset(hist,  0, 256*sizeof(int));
	memset(Sigma, 0, 256*sizeof(double));
	
	for(temp= InputImage; temp < end; temp++)
	{		
		hist[*temp]++;
	}

	memcpy(IntegralHist,  hist, 256*sizeof(int));
	memcpy(IntegralHistL, hist, 256*sizeof(int));

	for(IntegralHistL[0] = 0, i=1; i<256; i++)
	{
		IntegralHist[i] += IntegralHist[i-1];
		IntegralHistL[i] = i* IntegralHistL[i] + IntegralHistL[i-1];		
	}

	int Inverse1, Inverse2;
	double MaxV =0, m02 = (double) IntegralHistL[255] / (double)ImgSize;
	m02 *= m02;

	double ma, mb, pa, pb;
	for(i=1; i<255; i++)
	{		
		Inverse1 = ImgSize - IntegralHist[i];
		Inverse2 = IntegralHistL[255] - IntegralHistL[i];
		if(Inverse1 == 0 || IntegralHist[i] == 0) continue;

		pa = (double)IntegralHist[i]  / ImgSize;
		pb = (double)Inverse1		  / ImgSize;
		ma = (double)IntegralHistL[i] / IntegralHist[i];
		mb = (double)Inverse2		  / Inverse1;
		Sigma[i] = pa* ma* ma + pb* mb* mb - m02;

		if(Sigma[i] > MaxV)
		{
			MaxV = Sigma[i];
			Threshold = BYTE(i);
		}
		//TRACE(" %3d %5d %20.5f \n",i, hist[i], Sigma[i]);
	}

	TRACE("Threshold = %d (Qtsu method)  \n", Threshold);
	for(temp = InputImage; temp<end; temp++)
	{
		*temp = (*temp>Threshold) ? 255 : 0;
	}//*/
}
//***********************************************************
//函数类别： 二值图像的处理
//函数名称：                       
//           Image_AutoBinary_Qtsu
//函数用途：
//           Image_AutoBinary_EntropicCorrelation
//           利用<图像图形学报>99.5
//           "精子运动图像的多目标检测与分割"
//原始作者： 陆  宏  伟
//原始日期： 14/12/1999
//修改日期:  14/12/1999
//***********************************************************
void Image_AutoBinary_EntropicCorrelation(BYTE* InputImage , // 数据的指针
										  int ImgSize)       // 数据的大小
{
	if(InputImage==NULL) return;
	
	int    i, hist[256], ArrarySize =  256*sizeof(double);
	double Fai[256],  Hist1[256], Hist2[256];
	BYTE   Threshold, *temp, *end = InputImage + ImgSize;
	memset(hist, 0, 256*sizeof(int));
	memset(Fai,  0, ArrarySize);
	
	for(temp= InputImage; temp < end; temp++)
	{		
		hist[*temp]++;
	}
	for(i=0; i<256; i++)
	{
		Hist1[i]  = (double)hist[i] / (double)ImgSize;
		Hist2[i] = Hist1[i] * Hist1[i];
	}

	for(i=1; i<256; i++)
	{
		Hist1[i] += Hist1[i-1];
		Hist2[i] += Hist2[i-1];
	}

	double Inverse1, Inverse2;
	double MaxV =0;

	for(i=1; i<255; i++)
	{		

		Inverse1 = 1 - Hist1[i];
		Inverse2 = Hist2[255] - Hist2[i];
		if( Inverse1 < 1e-6 || Hist1[i] < 1e-6 ) continue;
		if( Inverse2 < 1e-6 || Hist2[i] < 1e-6 ) continue;

		Fai[i] = 2*log(Hist1[i] * Inverse1) - log(Hist2[i] * Inverse2);

		if(Fai[i] > MaxV)
		{
			MaxV = Fai[i];
			Threshold = BYTE(i);
		}
		//TRACE(" %3d %5d %20.5f \n",i, hist[i], Fai[i]);
	}

	TRACE("Threshold = %d (EntropicCorrelation method)  \n", Threshold);
	for(temp = InputImage; temp<end; temp++)
	{
		*temp = (*temp>Threshold) ? 255 : 0;
	}//*/
}
void Image_ManualBinary(BYTE* InputImage , int ImgSize,BYTE Threshold,  int Type)
{
	if(InputImage==NULL) return;
	
	BYTE *end = InputImage + ImgSize;
	switch (Type)
	{
		case 0: // >Threshold ? 255 : 0
		{
			for(BYTE *temp = InputImage; temp<end; temp++)
			{
				if(*temp>Threshold)
					*temp = 255;
				else
					*temp = 0;
			}
			break;
		}
		case 1: // >Threshold ? 0 : 255
		{
			for(BYTE *temp = InputImage; temp<end; temp++)
			{
				if(*temp>Threshold)
					*temp = 0;
				else
					*temp = 255;
			}
			break;
		}
		case 2: // >Threshold ? 255: 不变
		{
			for(BYTE *temp = InputImage; temp<end; temp++)
			{
				if(*temp>Threshold)
					*temp = 255;
			}
			break;
		}
		case 3: // >Threshold ? 不变 : 0
		{
			for(BYTE *temp = InputImage; temp<end; temp++)
			{
				if(*temp<Threshold)
					*temp = 0;
			}
			break;
		}
	}
}
void Image_Sobel(BYTE* InputImage , int ImgW , int ImgH)
{
	int	 Imagesize = ImgW*ImgH;
	if (InputImage==NULL || Imagesize==0) return ;

	int i,j;
	BYTE *Image ;
	BYTE *tmpImage = new BYTE[Imagesize];
    BYTE **RowAddress = new BYTE*[ImgH];
  
	//initialize 
	memset(tmpImage, 0, Imagesize);
	RowAddress[0] = InputImage;
	for(i=1; i<ImgH; i++)
	{
		RowAddress[i] = RowAddress[i-1] + ImgW;
	}

	double dx,dy;
	int    nG;
	Image = tmpImage + ImgW;
	for(j=1; j<ImgH-1; j++)
	{		
		for(i=1; i<ImgW-1; i++)
		{
			dx  = RowAddress[j-1][i+1] - RowAddress[j-1][i-1];
			dx += 2*(RowAddress[j][i+1] - RowAddress[j][i-1]);
			dx += RowAddress[j+1][i+1] - RowAddress[j+1][i-1];
			dy  = RowAddress[j+1][i-1] - RowAddress[j-1][i-1];
			dy += 2*(RowAddress[j+1][i] - RowAddress[j-1][i]);
			dy += RowAddress[j+1][i+1] - RowAddress[j-1][i+1];
			nG  = int(sqrt(dx*dx+dy*dy) + 0.5);
			Image[i] = (nG>255) ? 255: nG; 
		}
		Image += ImgW;
	}

	memcpy(InputImage, tmpImage, Imagesize);
	if(  tmpImage!=NULL)  delete []tmpImage;
	if(RowAddress!=NULL)  delete []RowAddress;
}

void Image_EdgeOrention_Sobel(BYTE* InputImage , int ImgW , int ImgH)
{
	int	 Imagesize = ImgW*ImgH;
	if (InputImage==NULL || Imagesize==0) return ;

	int i,j;
	BYTE *Image1, *Image2;
	BYTE *MagntiudeImage = new BYTE[Imagesize];
	BYTE *OrientationImage = new BYTE[Imagesize];
    BYTE **RowAddress = new BYTE*[ImgH];
  	//initialize 
	memset(MagntiudeImage, 0, Imagesize);
	memset(OrientationImage, 0, Imagesize);
	RowAddress[0] = InputImage;
	for(i=1; i<ImgH; i++)
	{
		RowAddress[i] = RowAddress[i-1] + ImgW;
	}

	double dx, dy, SlopeAngle;
	int    nG , dirNum = 16;
	Image1 = MagntiudeImage + ImgW;
	Image2 = OrientationImage + ImgW;
	for(j=1; j<ImgH-1; j++)
	{		
		for(i=1; i<ImgW-1; i++)
		{
			dx  = RowAddress[j-1][i+1] - RowAddress[j-1][i-1];
			dx += 2*(RowAddress[j][i+1] - RowAddress[j][i-1]);
			dx += RowAddress[j+1][i+1] - RowAddress[j+1][i-1];
			dy  = RowAddress[j+1][i-1] - RowAddress[j-1][i-1];
			dy += 2*(RowAddress[j+1][i] - RowAddress[j-1][i]);
			dy += RowAddress[j+1][i+1] - RowAddress[j-1][i+1];
			nG  = int(sqrt(dx*dx+dy*dy) + 0.5);
			Image1[i] = (nG>255) ? 255: nG;	
			if(nG<5||fabs(dx)<0.01)
			{
				Image2[i] = 0;
			}
			else 
			{
				SlopeAngle = atan2(dy, dx);
				if( SlopeAngle<0 ) SlopeAngle += 2*PI;
				BYTE WhichDegree = BYTE( dirNum* SlopeAngle/ PI+ 1.5);
				Image2[i] = WhichDegree<<2;
			}
		}
		Image1 += ImgW;
		Image2 += ImgW;
	}