commonproc.cpp

医学图象处理系统

	 
	if( (hFile = _findfirst(FileExtname, &c_file )) != -1L )
	{
		LastFilename = Path + c_file.name;	
		//TRACE("next image filename is %s\n",LastFilename);
		while( _findnext( hFile, &c_file ) == 0 )            
		{
			CString Nowfile = Path + c_file.name;
			if(Nowfile == Filename) OK=TRUE;
			if(OK)
			{
				Filename = LastFilename;
				return TRUE;
			}
			LastFilename = Nowfile;			
		}		 
	}
	_findclose( hFile );  
	return FALSE;
}
void L_TraceRect(CString VarName, CRect rect)
{
	TRACE("%s LeftUp point is (%d ,%d) , RightBottom point is (%d ,%d) w,h = (%d ,%d)\n",\
		   VarName, rect.left,  rect.top,  rect.right, rect.bottom, rect.Width(), rect.Height());
}

void pIntBubble(BYTE *p,int n)
{
	int m,k,j,i;
	BYTE d;
    k=0; m=n-1;
    while (k<m)
      { j=m-1; m=0;
        for (i=k; i<=j; i++)
          if (p[i]>p[i+1])
            { d=p[i]; p[i]=p[i+1]; p[i+1]=d; m=i;}
        j=k+1; k=0;
        for (i=m; i>=j; i--)
          if (p[i-1]>p[i])
            { d=p[i]; p[i]=p[i-1]; p[i-1]=d; k=i;}
      }
}

void IntSplit(BYTE *p,int n,int *m)
{
	int i,j,k,l;
	BYTE t;
    i=0; j=n-1;
    k=(i+j)/2;
    if ((p[i]>=p[j])&&(p[j]>=p[k])) l=j;
    else if ((p[i]>=p[k])&&(p[k]>=p[j])) l=k;
    else l=i;
    t=p[l]; p[l]=p[i];
    while (i!=j)
      { while ((i<j)&&(p[j]>=t)) j=j-1;
        if (i<j)
          { p[i]=p[j]; i=i+1;
            while ((i<j)&&(p[i]<=t)) i=i+1;
            if (i<j)
              { p[j]=p[i]; j=j-1;}
          }
      }
    p[i]=t; *m=i;
    return;
}

void IntQuickSort(BYTE *p,int n)
{
	int m,i0,*i;
	BYTE *s;
    i=&i0;
    if (n>10)
      { IntSplit(p,n,i);
        m=i0;
        IntQuickSort(p,m);

		s=p+(i0+1);
        m=n-(i0+1);
        IntQuickSort(s,m);
      }
    else pIntBubble(p,n);
}

void int_to_uchar(int *r,BYTE *in,int size)
{
	int i,
	max_r=r[0],
    min_r=r[0];
	double fTemp;
	for (i=1; i<size; i++)
	{
	    if ( r[i] > max_r )
			max_r=r[i];
	    else if ( r[i] < min_r )
			min_r=r[i];
	}

	max_r-=min_r;
	if(max_r!=0)
	{
		for (i=0; i<size; i++)
		{
			fTemp=(r[i]-min_r)*255.0/max_r+0.5; 
			in[i] = (BYTE)(fTemp);
		}
	}
}

void float_to_uchar(float *r,BYTE *in,int size)
{
	int i;
	float max_r=r[0];
    float min_r=r[0];
	double fTemp;
	for (i=1; i<size; i++)
	{
	    if ( r[i] > max_r )
			max_r=r[i];
	    else if ( r[i] < min_r )
			min_r=r[i];
	}

	max_r-=min_r;
	if(max_r!=0)
	{
		for (i=0; i<size; i++)
		{
			fTemp=(r[i]-min_r)*255.0/max_r+0.5; 
			in[i] = (BYTE)(fTemp);
		}
	}
}

// VGaussianKernel Generate a 1D Gaussian kernel.
BOOL VGaussianKernel (int ncoeffs, float *coeffs, float s)
{
    float sum, x, b = (float)( -1.0 / (2.0 * s * s) );
    float *p1, *p2, *pend = coeffs + ncoeffs;

    // Setup depends on whether the number of coefficients asked for is odd or even:
    if (ncoeffs & 1) // ncoeffs is a odd number
	{
		p1    = & coeffs[ncoeffs / 2];
		*p1-- = 1.0;
		p2    = p1 + 2;
		x     = 1.0;
		sum   = 0.5;
    } 
	else
	{
		p1  = & coeffs[ncoeffs / 2];
		p2  = p1 + 1;
		x   = 0.5;
		sum = 0.0;
    }

    // Fill the vector with coefficients, then make them sum to 1.0:
    while (p2 < pend) 
	{
		sum += *p1-- = *p2++ = (float)exp (x * x * b);
		x += 1.0;
    }
    sum += sum;
    for (p1 = coeffs; p1 < pend; p1++)
		*p1 /= sum;

    return TRUE;
}
// VGaussianD1Kernel Generate a kernel for the first derivative of a 1D Gaussian.
BOOL VGaussianD1Kernel (int ncoeffs, float *coeffs, float s)
{
    float x, a = (float)( 1.0 / (sqrt (2.0 * PI) * s * s * s) );
    float b = (float)( -1.0 / (2.0 * s * s) );
    float *p1, *p2, *pend = coeffs + ncoeffs;

    // Setup depends on whether the number of coefficients asked for is odd or even:
    if (ncoeffs & 1) 
	{
		p1 =  coeffs + ncoeffs/ 2;
		*p1-- = 0.0;
		p2 = p1 + 2;
		x = 1.0;
    }
	else 
	{
		p1 = & coeffs[ncoeffs / 2];
		p2 = p1 + 1;
		x = 0.5;
    }
    
    while (p2 < pend) // Fill the vector with coefficients:
	{
		/*float tx = x - 0.5f;  // this code is get derivation from integration of small area
		*p2 = 0;
		for(int i=0; i <11; i++,tx += 0.1f)
		{
			*p2 += a * tx * (float)exp (tx * tx * b);
		}
		*p2 /= 10; //*/
		*p2 = a * x * (float)exp (x * x * b);
		*p1-- = -*p2++;
		x += 1.0;
    }
    return TRUE;
}
//  VGaussianSplineKernel
//
//  Generate a 1D Gaussian kernel with its tails splined to zero.
//  The kernel is Gaussian over the domain [-2s,2s].
//  Over [-4s,-2s) and over (2s,4s], it splines to zero.
//  Beyond -4s and 4s it is zero.
//  Thus ncoeffs should be at least 8s + 1 to get a smoothly splined kernel.
// 
//  Based on Mark Nitzberg, David Mumford, Takahiro Shiota, `Filtering,
//  Segmentation and Depth', Lecture Notes in Computer Science 662,
//  Springer-Verlag, 1993.
BOOL VGaussianSplineKernel (int ncoeffs, float *coeffs, float s)
{
    float sum, x, t, b = (float)( -1.0 / (2.0 * s * s) );
    float k2 = (float)( 2.0 * s ) , k4 = (float) ( 4.0 * s );
    float *p1, *p2, *pend = coeffs + ncoeffs;

    // Setup depends on whether the number of coefficients asked for is odd or even:
    if (ncoeffs & 1) 
	{
		p1    = & coeffs[ncoeffs / 2];
		*p1-- = 1.0;
		p2    = p1 + 2;
		x     = 1.0;
		sum   = 0.5;
    }
	else
	{
		p1  = & coeffs[ncoeffs / 2];
		p2  = p1 + 1;
		x   = 0.5;
		sum = 0.0;
    }

    // Fill the vector with coefficients, then make them sum to 1.0:
    while (p2 < pend) 
	{
		if (x <= k2)
			t = (float)exp (x * x * b);
		else if (x <= k4) 
		{
			t  = (float)( 4.0 - x / s );
			t *= (float)( 1.0 / (16.0 * M_E * M_E) * t * t * t );
		} 
		else 
			t = 0.0;
		sum += *p1-- = *p2++ = t;
		x   += 1.0;
    }
    sum += sum;
    for (p1 = coeffs; p1 < pend; p1++)
		*p1 /= sum;

    return TRUE;
}

//***********************************************************
//函数类别： 图像边缘提取
//函数名称：                       
//           Make_Rfilter_Mask
//函数用途：
//           产生R filter[ R(x) = exp(-|x|/λ)/(2λ) ] 
//			 的原始和梯度模板 
//说    明:  利用 Muhittin Gokmen 的文章" λτ-Space Representation 
//			 of Images and Geeralized Edge Detector", PAMI Vol.19.NO.6 JUNE 1997
//原始作者： 陆  宏  伟
//原始日期： 19/12/1999
//***********************************************************
void Make_Rfilter_Mask(int masksize, float lambda,  
					   float *lpDmask, float *lpmask, float maxresponse)
{
    int	  i, hM= masksize/ 2;    
    float gN0 = 0.5f* maxresponse/ lambda;			 // 
    float gD0 = 0.5f* maxresponse/(lambda * lambda); //

    for(lpDmask[hM] = 0.0f, i = -hM; i < 0; i++)
    {
		lpDmask[hM+ i] = -float( gD0* exp(i/ lambda) );
		lpDmask[hM- i] = -lpDmask[hM+ i];
    }
	if( lpmask != NULL )
	{		
		for(lpDmask[hM] = gN0, i = -hM; i < 0; i++)
			lpmask[hM+ i] = lpmask[hM- i] = float( gN0 * exp(i/ lambda) );
    }
}
//***********************************************************
//函数类别： 图像边缘提取
//函数名称：                       
//           Make_Second_Rfilter_Mask
//函数用途：
//           产生2 order R filter[ R(x) = exp(-|x|/sqrt(2)/pow(λ,0.25))
//                    /(2*pow(λ,0.25))*cos(-|x|/sqrt(2)/pow(λ,0.25)- π/4)]
//			 的原始和梯度模板 
//说    明:  利用 Muhittin Gokmen 的文章" λτ-Space Representation 
//			 of Images and Geeralized Edge Detector", PAMI Vol.19.NO.6 JUNE 1997
//原始作者： 陆  宏  伟
//原始日期： 19/12/1999
//***********************************************************
void Make_Second_Rfilter_Mask(int masksize, float lambda, float tau, 
							  float *lpDmask, float *lpmask, float maxresponse)
{
	int   i, hM= masksize/ 2;
    float Tg0, Tg1;
    
    float lambdaq    = float(sqrt(sqrt(lambda)));
    float lambdasqrt = float(sqrt(lambda));
    float cons2      = float(sqrt(2.0)* lambdaq);

    float gN0 = 0.5f* maxresponse/ lambda;
    float gN1 = 0.5f* maxresponse/ lambdaq;
    float gD0 = 0.5f* maxresponse/ (lambda * lambda);
    float gD1 = 0.5f* maxresponse/ lambdasqrt;

	for(lpDmask[hM] = 0.0f, i = -hM; i < 0; i++)
    {
		Tg0 = float( gD0 * exp(i/ lambda) );
		Tg1 = float( gD1 * exp(i/ cons2) * cos( i/ cons2) );
		lpDmask[hM- i] = (1.0f - tau) * Tg0 + tau * Tg1;
		lpDmask[hM+ i] = -lpDmask[hM- i];
    }
	if( lpmask != NULL )
	{		
		for(lpDmask[hM] = gN0, i = -hM; i < 0; i++)
		{
			Tg0 = float( gN0 * exp(i/ lambda) );
			Tg1 = float( gN1 * exp(i/ cons2) * cos(-i/ cons2 - PI/4.0) );
			lpmask[hM+ i] = lpmask[hM- i] = (1.0f - tau) * Tg0 + tau * Tg1;
		}
    }
}
//***********************************************************
//函数类别： 图像边缘提取
//函数名称：                       
//           Make_GED_filter_mask
//函数用途：
//           产生General Edge Detector filter 的原始和梯度模板 
//说    明:  利用 Muhittin Gokmen 的文章" λτ-Space Representation 
//			 of Images and Geeralized Edge Detector", PAMI Vol.19.NO.6 JUNE 1997
//原始作者： 陆  宏  伟
//原始日期： 20/12/1999  今天澳门回归了!!!
//***********************************************************
void Make_GED_filter_mask(int masksize, float lambda, float tau, 
						  float *lpDmask, float *lpmask, float maxresponse)
{
    int   i, hM= masksize/ 2;
    float gN0, gD0, d1, d2, dSqrt, x, quarta;
    float a = lambda * tau;
    float b = lambda * ( 1.0f - tau);
    float discrim = b * b - 4.0f * a;
    float abs_discrim = (discrim > 0.0) ? discrim : -discrim;
    
    if (abs_discrim < 0.01) // case III
	{	
		d1 = (float)sqrt(b / (2.0 * a));
		d2 = 1.0f/ d1;

		gN0 = 0.5f* maxresponse/ b;
		gD0 = 0.5f* d1 * maxresponse / b;

		for(lpDmask[hM] = 0.0f, i = -hM; i < 0; i++)
		{
			lpDmask[hM+ i] = float( gD0 * i * exp(i* d1) );
			lpDmask[hM- i] = -lpDmask[hM+ i];
		}
		if( lpmask != NULL )
		{		
			for(lpDmask[hM] = gN0, i = -hM; i < 0; i++)
				lpmask[hM+ i] = lpmask[hM- i] = float( gN0 * exp(i* d1) * ( d2 - i) );
		}
	}
    else if (fabs(a) < 0.000001) // case V -  membrane 
	{
		d1  = (float)sqrt(b);
		gN0 = maxresponse / d1;
		gD0 = maxresponse / b;

		for(lpDmask[hM] = 0.0f, i = -hM; i < 0; i++)
		{
			lpDmask[hM+ i] = -float( gD0 * exp(i/ d1) );
			lpDmask[hM- i] = -lpDmask[hM+ i];
		}
		if( lpmask != NULL )
		{		
			for(lpDmask[hM] = gN0, i = -hM; i < 0; i++)
				lpmask[hM+ i] = lpmask[hM- i] = float( gN0 * exp(i/ d1) );
		}
	}
    else if (fabs(b) < 0.000001) // case IV -  Plate  
	{
		dSqrt = (float)sqrt(a);
		d2	  = (float)sqrt(dSqrt);
		d1    = 1.0f / (float)( sqrt(2.0) * d2);
		gN0   = maxresponse / (4.0f * dSqrt * d1);
		gD0   = maxresponse / (2.0f * dSqrt);

		for(lpDmask[hM] = 0.0f, x= -d1, i = 1; i <= hM; i++, x -= d1)
		{
			lpDmask[hM- i] = float( gD0 * exp(x)* sin(x) );
			lpDmask[hM+ i] = -lpDmask[hM- i];
		}
		if( lpmask != NULL )
		{		
			for(lpDmask[hM] = gN0, x= -d1, i = 1; i <= hM; i++, x -= d1)
				lpmask[hM+ i] = lpmask[hM- i] = float( gN0 * exp(x)* (cos(x) - sin(x)) );
		}
	} 
    else if (discrim > 0.0 ) // case I
	{
		d1 = (float)sqrt((b + sqrt(discrim))/(2.0 * a));
		d2 = (float)sqrt((b - sqrt(discrim))/(2.0 * a));

		gN0 = 1.0f * maxresponse/(2.0f* a* (d2* d2 - d1* d1));
		gD0 = gN0;
    
		for(lpDmask[hM] = 0.0f, i = -hM; i < 0; i++)
		{
			lpDmask[hM+ i] = -float( gD0 * (exp(i* d2) - exp(i* d1)) );
			lpDmask[hM- i] = -lpDmask[hM+ i];
		}
		if( lpmask != NULL )
		{		
			for(lpDmask[hM] = gN0, i = -hM; i < 0; i++)
				lpmask[hM+ i] = lpmask[hM- i] = float( gN0 * ( exp(i* d1)/ d1 - exp(i* d2)/ d2 ) );
		}
	} 
    else if (discrim < 0.0) // case II
	{
		dSqrt  = (float)sqrt(a);
        quarta = (float)sqrt(dSqrt);
		d1	   = (float)( cos(0.5* atan(sqrt(4.0*a/(b*b) - 1)))) / quarta ;
		d2     = (float)( sin(0.5* atan(sqrt(4.0*a/(b*b) - 1)))) / quarta ;
		gN0    = maxresponse/ (4.0f * dSqrt);
		gD0    = maxresponse/ (4.0f * a * d1 * d2);
    
		for(lpDmask[hM] = 0.0f, i = -hM; i < 0; i++)
		{
			lpDmask[hM+ i] = float( gD0 * exp(i* d2) * sin(i* d1) );
			lpDmask[hM- i] = -lpDmask[hM+ i];
		}
		if( lpmask != NULL )
		{		
			for(lpDmask[hM] = gN0, i = -hM; i < 0; i++)
				lpmask[hM+ i] = lpmask[hM- i] = float( gN0 * exp(i* d2)*( cos(i* d1)/ d2 + sin(i* d1)/ d1 ) );
		}
	} 
}

int GetDriveCount()
{
    int  i, iDriveType, cDrives = 0;
    char szDrive[] = "x:\\\0";

    for (szDrive[0]='C'; szDrive[0] <= 'Z'; szDrive[0]++)
    {
        i = (int) (szDrive[0] - 'C');
        iDriveType = ::GetDriveType( szDrive );

        if (iDriveType == DRIVE_FIXED   ||
            iDriveType == DRIVE_REMOTE  ||
            iDriveType == DRIVE_RAMDISK)
        {
            cDrives++;
        }
    }
    return cDrives;
}