imageprocess.cpp

医学图象处理系统

		}
	}

	if(  tmpImage!=NULL)  delete []tmpImage;
	if(S_RowAddress!=NULL)  delete []S_RowAddress;
	if(D_RowAddress!=NULL)  delete []D_RowAddress;
}
void Image_NonRMiddleFilter(BYTE* InputImage , int ImgW , int ImgH ,int Masksize)
{
	int	 Imagesize = ImgW*ImgH;
	if (InputImage==NULL || Imagesize==0) return ;
	
	int  i,j;
	int  Area = Masksize*Masksize; 
	int  HalfArea = Area/2;
	int	 r=Masksize/2;  
	BYTE *tmpImage = new BYTE[Imagesize] ,  *Image ;
    BYTE **RowAddress = new BYTE*[ImgH];
    BYTE *blk = new BYTE[Area]; 
	//initialize 
	memcpy(tmpImage, InputImage, Imagesize);

	RowAddress[0] = tmpImage;
	for(i=1 ; i<ImgH ;i++)
	{
		RowAddress[i] = RowAddress[i-1] + ImgW;
	}

	CProgressDlg * MsgHandler = new CProgressDlg();
    MsgHandler->Create(NULL);
    MsgHandler->SetStatus("Nonrecursive Median Filter Process");

	int dx,dy,k;;
	Image = InputImage + r*ImgW;
	for(j=r; j<ImgH-r; j++)
	{		
		for(i=r; i<ImgW-r; i++)
		{
			k=0;
			for(dy=-r; dy<=r; dy++)
				for(dx=-r; dx<=r; dx++)
					blk[k++] =  RowAddress[j+dy][i+dx]; 
			IntQuickSort(blk, Area);
			Image[i] = blk[HalfArea];				 
		}
		Image += ImgW;
		MsgHandler->SetPos(100 * j /ImgH );
	}
	if(       blk!=NULL)  delete []blk;
	if(  tmpImage!=NULL)  delete []tmpImage;
	if(RowAddress!=NULL)  delete RowAddress;
	MsgHandler->DestroyWindow();
    delete  MsgHandler;
}
void Image_RMiddleFilter(BYTE* InputImage , int ImgW , int ImgH ,int Masksize)
{
	int	 Imagesize = ImgW*ImgH;
	if (InputImage==NULL || Imagesize==0) return ;

	int  i,j;
	int  Area = Masksize*Masksize; 
	int  HalfArea = Area/2;
	int	 r=Masksize/2;  
    BYTE **RowAddress = new BYTE*[ImgH];
    BYTE *blk = new BYTE[Area]; 
	//initialize 

	RowAddress[0] = InputImage;
	for(i=1 ; i<ImgH ;i++)
	{
		RowAddress[i] = RowAddress[i-1] + ImgW;
	}

	CProgressDlg * MsgHandler = new CProgressDlg();
    MsgHandler->Create(NULL);
    MsgHandler->SetStatus("Recursive Median Filter Process");

	int dx,dy,k;;
	for(j=r; j<ImgH-r; j++)
	{		
		for(i=r; i<ImgW-r; i++)
		{
			k=0;
			for(dy=-r; dy<=r; dy++)
				for(dx=-r; dx<=r; dx++)
					blk[k++] =  RowAddress[j+dy][i+dx]; 
			IntQuickSort(blk, Area);
			RowAddress[j][i] = blk[HalfArea];				 
		}
		MsgHandler->SetPos(100 * j /ImgH );
	}
	if(       blk!=NULL)  delete []blk;
	if(RowAddress!=NULL)  delete RowAddress;
	MsgHandler->DestroyWindow();
    delete  MsgHandler;
}
//The fast median theory is that remove left colomn and add right colomn
void Image_FastMedian (BYTE* InputImage , int ImgW , int ImgH ,int Masksize)
{
	int	 Imagesize = ImgW*ImgH;
	if (InputImage==NULL || Imagesize==0) return ;

	int  Area = Masksize*Masksize; 
	int  level = Area/2;
	int	 r=Masksize/2;
	BYTE *tmpImage = new BYTE[Imagesize] ,  *Image ;
    BYTE **RowAddress = new BYTE*[ImgH];
	memcpy(tmpImage, InputImage, Imagesize);
	RowAddress[0] = tmpImage;
	for(int i=1 ; i<ImgH ;i++)
	{
		RowAddress[i] = RowAddress[i-1] + ImgW;
	}

	CProgressDlg * MsgHandler = new CProgressDlg();
    MsgHandler->Create(NULL);
    MsgHandler->SetStatus("Fast Median Filter Process");

    int classf[256];
	int ix, iy, p, fx, fy;
    int fax, fex, rank, count;

	int iax = r;
	int iex = ImgW - r;
	int iay = r;
	int iey = ImgH - r;
	int fay = 0;
	int fey = Masksize ;

	Image = InputImage + iay*ImgW;
	for (iy = iay; iy < iey; iy++)
	{
		fax = count = rank = 0;
		fex = Masksize;
		memset(classf,0,256*sizeof(int));
		
		for (fy = fay; fy < fey; fy++)        
		{
			for (fx = fax; fx < fex; fx++)
			{
				p = RowAddress[fy][fx];
				classf[p]++;
			}
		}
		while(count<level)
		{
			count+=classf[rank++];
		}
		rank--;

		for (ix = iax; ix < iex ; ix++)      
		{
			Image[ix] = (BYTE)rank;              
			for (fy = fay; fy < fey; fy++)    
			{
				p = RowAddress[fy][fex + 1];
				classf[p]++; 	             // add a pixel         
				if (p <= rank)   count++;
				p = RowAddress[fy][fax];
				classf[p]--;                 // remove a pixel 
				if (p <= rank)   count--;
			}	  
			while(count  >= level)
		  		count -= classf[rank--];
			while(count < level)
				count += classf[++rank];
			fax++;
			fex++;
		}
		Image[iex] = (BYTE)rank;
		fay++;
		fey++;
		Image += ImgW;
		MsgHandler->SetPos(100 * iy /ImgH );
	}
	if(  tmpImage!=NULL)  delete []tmpImage;
	if(RowAddress!=NULL)  delete []RowAddress;
	MsgHandler->DestroyWindow();
    delete  MsgHandler;
}
//***********************************************************
//函数类别： 图像的边缘提取
//函数名称：  
//           Image_SUSAN_Edge_Normal
//函数用途：
//           图像的边缘提取；
//
//参数说明：
//           CImage   *s_img  图像的对象指针
//           int       bt     阈值
//           BYTE   *mid    边缘信息的输出的指针
//说    明：
//           利用SUSAN边缘检测的原理编制的
//原始作者： 不 是 我
//原始日期： 3/20/1999 
//***********************************************************
void Image_SUSAN_Edge_Normal(BYTE* InputImage , int ImgW , int ImgH ,int bt)
{
	int	 Imagesize=ImgW*ImgH;
	if (InputImage==NULL || Imagesize==0) return ;
	
	int    do_symmetry, i, j, m, n, a, b, x, y, w,max_no=1850;
	int    Width=ImgW, Height=ImgH;
    int    *r  =new int[Imagesize];
	BYTE   *mid = new BYTE[Imagesize];
	BYTE   *in=InputImage;
	BYTE   c,*p,*cp,*bp,bbp[516];
    double temp,z; 

    bbp[0]=bbp[1]=bbp[515]=bbp[514]=0;
    bp=bbp+258;
   
    for(int k=-256;k<257;k++)
	{
        temp=((float)k)/((float)bt);
        temp=temp*temp;
		temp=temp*temp*temp;
        temp=100.0*exp(-temp)+0.5;
        bbp[k+258] = (BYTE)(temp);
	}
    memset(mid,100, Width*Height);
    memset(r,  0,   Width*Height * sizeof(int));

    for(i=3;i<Height-3;i++)
    {
		for(j=3;j<Width-3;j++)
		{
			n=100;
			p=in + (i-3)*Width + j - 1;
			cp=bp + in[i*Width+j];

			n+=*(cp-*p++);
			n+=*(cp-*p++);
			n+=*(cp-*p);
			p+=Width-3; 

			n+=*(cp-*p++);
			n+=*(cp-*p++);
			n+=*(cp-*p++);
			n+=*(cp-*p++);
			n+=*(cp-*p);
			p+=Width-5;

			n+=*(cp-*p++);
			n+=*(cp-*p++);
			n+=*(cp-*p++);
			n+=*(cp-*p++);
			n+=*(cp-*p++);
			n+=*(cp-*p++);
			n+=*(cp-*p);
			p+=Width-6;

			n+=*(cp-*p++);
			n+=*(cp-*p++);
			n+=*(cp-*p);
			p+=2;
			n+=*(cp-*p++);
			n+=*(cp-*p++);
			n+=*(cp-*p);
			p+=Width-6;

			n+=*(cp-*p++);
			n+=*(cp-*p++);
			n+=*(cp-*p++);
			n+=*(cp-*p++);
			n+=*(cp-*p++);
			n+=*(cp-*p++);
			n+=*(cp-*p);
			p+=Width-5;

			n+=*(cp-*p++);
			n+=*(cp-*p++);
			n+=*(cp-*p++);
			n+=*(cp-*p++);
			n+=*(cp-*p);
			p+=Width-3;

			n+=*(cp-*p++);
			n+=*(cp-*p++);
			n+=*(cp-*p);

			if (n<=max_no)
			   r[i*Width+j] = max_no - n;
	   }
	}

	for(i=4;i<Height-4;i++)
	{
		for(j=4;j<Width-4;j++)
		{
			if (r[i*Width+j]>0)
			{
				m=r[i*Width+j];
				n=max_no - m;
				cp=bp + in[i*Width+j];

				if (n>600)
				{
					p=in + (i-3)*Width + j - 1;
					x=0;y=0;
 
   					c=*(cp-*p++);x-=c;y-=3*c;
					c=*(cp-*p++);y-=3*c;
					c=*(cp-*p);x+=c;y-=3*c;
					p+=Width-3; 
    
					c=*(cp-*p++);x-=2*c;y-=2*c;
				    c=*(cp-*p++);x-=c;y-=2*c;
				    c=*(cp-*p++);y-=2*c;
				    c=*(cp-*p++);x+=c;y-=2*c;
				    c=*(cp-*p);x+=2*c;y-=2*c;
				    p+=Width-5;
    
			        c=*(cp-*p++);x-=3*c;y-=c;
			        c=*(cp-*p++);x-=2*c;y-=c;
			        c=*(cp-*p++);x-=c;y-=c;
			        c=*(cp-*p++);y-=c;
					c=*(cp-*p++);x+=c;y-=c;
					c=*(cp-*p++);x+=2*c;y-=c;
					c=*(cp-*p);x+=3*c;y-=c;
					p+=Width-6;

					c=*(cp-*p++);x-=3*c;
					c=*(cp-*p++);x-=2*c;
					c=*(cp-*p);x-=c;
					p+=2;
					c=*(cp-*p++);x+=c;
					c=*(cp-*p++);x+=2*c;
					c=*(cp-*p);x+=3*c;
					p+=Width-6;
    
					c=*(cp-*p++);x-=3*c;y+=c;
					c=*(cp-*p++);x-=2*c;y+=c;
					c=*(cp-*p++);x-=c;y+=c;
					c=*(cp-*p++);y+=c;
					c=*(cp-*p++);x+=c;y+=c;
					c=*(cp-*p++);x+=2*c;y+=c;
					c=*(cp-*p);x+=3*c;y+=c;
					p+=Width-5;

					c=*(cp-*p++);x-=2*c;y+=2*c;
					c=*(cp-*p++);x-=c;y+=2*c;
					c=*(cp-*p++);y+=2*c;
					c=*(cp-*p++);x+=c;y+=2*c;
					c=*(cp-*p);x+=2*c;y+=2*c;
					p+=Width-3;

					c=*(cp-*p++);x-=c;y+=3*c;
					c=*(cp-*p++);y+=3*c;
					c=*(cp-*p);x+=c;y+=3*c;

			        z = sqrt((float)((x*x) + (y*y)));
			        if (z > (0.9*(float)n)) /* 0.5 */
					{
				        do_symmetry=0;
                        if (x==0)    z=1000000.0;
					        else     z=((float)y) / ((float)x);
						if (z < 0) { z=-z; w=-1; }
						    else     w=1;
						if (z < 0.5) { /* vert_edge */ a=0; b=1; }
						else 
						{
							if (z > 2.0) 
							{ /* hor_edge */ 
								a=1; b=0; 
							}
						    else 
							{ /* diag_edge */ 
								if (w>0) { a=1; b=1; }
								else { a=-1; b=1; }
							}
						}
						if( (m > r[(i+a)*Width+j+b]) && 
							(m >= r[(i-a)*Width+j-b]) &&
							(m > r[(i+(2*a))*Width+j+(2*b)]) 
							&& (m >= r[(i-(2*a))*Width+j-(2*b)]) )
						    mid[i*Width+j] = 1;
					}
                    else
						do_symmetry=1;
				}
                else
					do_symmetry=1;

				if (do_symmetry==1)
				{
					p=in + (i-3)*Width + j - 1;
				    x=0; y=0; w=0;

				    /*   |      \
					     y  -x-  w
					     |        \   */

				    c=*(cp-*p++);x+=c;y+=9*c;w+=3*c;
				    c=*(cp-*p++);y+=9*c;
				    c=*(cp-*p);x+=c;y+=9*c;w-=3*c;
				    p+=Width-3; 
  
				    c=*(cp-*p++);x+=4*c;y+=4*c;w+=4*c;
				    c=*(cp-*p++);x+=c;y+=4*c;w+=2*c;
				    c=*(cp-*p++);y+=4*c;
				    c=*(cp-*p++);x+=c;y+=4*c;w-=2*c;
				    c=*(cp-*p);x+=4*c;y+=4*c;w-=4*c;
				    p+=Width-5;
    
				    c=*(cp-*p++);x+=9*c;y+=c;w+=3*c;
				    c=*(cp-*p++);x+=4*c;y+=c;w+=2*c;
				    c=*(cp-*p++);x+=c;y+=c;w+=c;
				    c=*(cp-*p++);y+=c;
				    c=*(cp-*p++);x+=c;y+=c;w-=c;
				    c=*(cp-*p++);x+=4*c;y+=c;w-=2*c;
				    c=*(cp-*p);x+=9*c;y+=c;w-=3*c;
				    p+=Width-6;

				    c=*(cp-*p++);x+=9*c;
				    c=*(cp-*p++);x+=4*c;
				    c=*(cp-*p);x+=c;
				    p+=2;
				    c=*(cp-*p++);x+=c;
				    c=*(cp-*p++);x+=4*c;
				    c=*(cp-*p);x+=9*c;
				    p+=Width-6;
    
				    c=*(cp-*p++);x+=9*c;y+=c;w-=3*c;
				    c=*(cp-*p++);x+=4*c;y+=c;w-=2*c;
				    c=*(cp-*p++);x+=c;y+=c;w-=c;
				    c=*(cp-*p++);y+=c;
				    c=*(cp-*p++);x+=c;y+=c;w+=c;
				    c=*(cp-*p++);x+=4*c;y+=c;w+=2*c;
				    c=*(cp-*p);x+=9*c;y+=c;w+=3*c;
				    p+=Width-5;
 
				    c=*(cp-*p++);x+=4*c;y+=4*c;w-=4*c;
				    c=*(cp-*p++);x+=c;y+=4*c;w-=2*c;
				    c=*(cp-*p++);y+=4*c;
				    c=*(cp-*p++);x+=c;y+=4*c;w+=2*c;
				    c=*(cp-*p);x+=4*c;y+=4*c;w+=4*c;
				    p+=Width-3;

				    c=*(cp-*p++);x+=c;y+=9*c;w-=3*c;
				    c=*(cp-*p++);y+=9*c;
				    c=*(cp-*p);x+=c;y+=9*c;w+=3*c;

				    if (y==0)
					    z = 1000000.0;
				    else
					    z = ((float)x) / ((float)y);
				    if (z < 0.5) { /* vertical */ a=0; b=1; }
				    else { if (z > 2.0) { /* horizontal */ a=1; b=0; }
				    else { /* diagonal */ if (w>0) { a=-1; b=1; }
										else { a=1; b=1; }}}
				    if ( (m > r[(i+a)*Width+j+b]) && 
						 (m >= r[(i-a)*Width+j-b]) &&
					     (m > r[(i+(2*a))*Width+j+(2*b)]) && 
						 (m >= r[(i-(2*a))*Width+j-(2*b)]) )
					mid[i*Width+j] = 2;	
				}
			}
        }
	} 
	Image_SUSAN_thin(r,mid,Width,Height);
	for(i=0;i<Imagesize;i++)
	{
		if(mid[i]<10)
		{
			BYTE *p = InputImage + i - Width - 2;
			*p++ = 255;  *p++ = 255;  *p++ = 255;  p+= Width - 2;
			*p++ = 255;  *p++ = 255;  *p++ = 255;  p+= Width - 2;
			*p++ = 255;  *p++ = 255;  *p++ = 255; 
		}
	}//*/
	for(i=0;i<Imagesize;i++)
	{
		if(mid[i]<10)