GLCM(Gray-Level Co-occurrence Matrix) & Haralick texture features

o GLCM ( Gray-Level Co-occurrence Matrix)

항상 느끼는 것이지만 알고리즘은 공부하고 사용하지 않다보면 잘 잊어버리게되네요.

이것도 예전에 공부했던 내용인데, 다시 보니 기억이....

그래서 결국 이참에 정리를 하게 됩니다.

먼저 GLCM 구하는 방법을 보면 다름과 같습니다.

위의 그림 하나면 설명이 끝납니다.

좌측 메트릭스는 영상이라고 생각해봅시다. 영상의 픽셀 값은 설명을 위해 0~3까지 총 4단계로 구성되죠.

따라서 GLCM 도 4 * 4 메트릭스가 나오는 것입니다.

먼저 가로로 두개씩 보는 것입니다. 예를들어, 좌측 그림에서 가장 상단의 좌측 픽셀 두개는 0, 1 입니다.

이것의 의미는 (0, 1) 이라는 의미, 즉, x = 0, y= 1 이라는 의미입니다. 이 좌표는 GLCM의 위치를

의미합니다.

해당 위치에 +1을 카운트 해주면 됩니다.

GLCM의 좌측 상단 첫번째 값은 2가 됩니다. 왜냐하면 좌측 영상에서 0, 0이 2개가 존재하기 때문이죠.

매우 간단하죠?

여기서 한가지 더 크기에 달라지면 카운트 값이 많아지고 적어지고 하겠죠? 잘 보시면 GLCM의 합은

12 입니다. 왜냐하면 가로로 짝지어서 좌표를 생성하기 떄문에 가로로는 총 3개를 추출할 수 있고,

세로로는 4개를 추출 할 수 있기 때문이죠. 따라서 3*4는 12.   이것은 정규화 하는 과정을 수행합니다.

정규화는 당연히 합이 1이 되도록 만드는 것이므로 모든 GLCM의 값을 12로 나누면 되겠죠.

막상 설명하려니 힘드네요. 참 쉬운데... 어찌 쉽게 타자로 설명할 방~~~법이 없네요.ㅋㅋ

 

※ 아래 예제는 참고만 하시기 바랍니다. 제가 직접 코딩해보니 약간 틀린 곳이 있습니다.

o Haralick Texture Features

사실 지금 설명하는 것은 현재 보고 있는 문서 기준이기 때문에 정확하게

이것이 Haralick Texture Features이다 라고 말할 순 없지만, 어느정도 다른 사람들이

참고하여 공부하는데 도움이 되었으면 합니다.

1. Angular Second Moment (ASM) Feature

 

수식이 어렵게 나와있지만, 매우 간단하다. GLCM의 각 위치의 값을 제곱하여 다 더한 값이 ASM의 값이다.

2. Contrast Feature

 

위 수식의 표현은 참 애매하다. 하지만 알고보면 쉽다는거..

 

 

위의 그림은 GLCM의 예제이다. a)는 입력 영상이고, b)는 GLCM을 구한 것이다. c)는 GLCM을 정규화한 것이다. 아마 잘 보면 b)가 어떻게 나왔는지. 이해가 안될수도 있을 것이다. 왜냐하면 한쪽 방향에 대해서 구한 것이 아니라 두개 방향에 대해 누적하여 생성한 것이기 떄문이다. 세로 양방향에 대해서 누적하면 다음과 같이 나온다. 즉 (y,y+1)과 (y+1,y)의 방향에 대해서 누적한 결과가 b)인 것이다.

위의 예제는 contrast를 구한 예제인데, 보면 수식과는 계산 방법이 조금 다름을 알 수 있다. 결국 풀어서 생각하면 위의 예제와 같다. 왜냐하면 수식대로 n=|i-j| 이기 때문에 i와 j의 차이가 같은 GLCM의 값들을 모두 더한 후 n^2을 곱하여 모두 더한 결과이기 때문이다. Contrast는 대비를 뜻하는데, 이 수식이 왜 대비를 의미하는지는 한번 생각해보기 바란다. 힌트를 주자면 GLCM의 어떤 임의의 위치의 값은 무엇을 의미하는지 잘 생각해보자. 대각성분의 위치에 있는 값들은 평평한 픽셀의 빈도수를 의미하지 않을까?

3. Entropy Feature

 

엔트로피 특징의 수식은 다음과 같다. 엔트로피는 불확실성을 의미한다고 생각하면 쉽다. 예를 들어 동전을 던져 앞이 나올 확률은 0.5, 뒤가 나올 확률도 0.5이다. 이것을 엔트로피 구하는 수식에 대입해보면, 값이 1이 나올 것이다. 결국 불확실성이 크다는 의미이다. 어떤 문제에서 한 쪽으로 치우쳐져 있다면, 엔트로피는 작을 것이다. 이 부분도 한번 생각해보아야 할 문제이다. 엔트로피가 높다는 것은 무슨 의미를 갖는지...

 

4. Variance Feature

 

분산 특징도 크게 어려운 것은 없다. 한가지 의문을 갖자면, 왜 i에 대해서만 구하는가 이다. 또 뮤는 어떻게 구하는지..

위의 예제는 위에서 제시한 GLCM에 대해 mean과 variance의 계산 예제이다. 왜 mean과 variance를 i에 대해서만 구하는지 고민해봐야겟다. 이 값이 의미하는 것이 무엇인지도 함께 고민해보자.

 

5. Correlation Feature

이 특징은 GLCM의 선형적인 상관관계를 보여주는 값이다.

 

 

 

 

Correlation 특징은 위의 수식대로만 풀면 간단하게 구할 수 있다. 이 부분은 특별히 설명할 부분이 없는 듯하다. 한가지 코멘트를 달자면, 좀 특이한 경우인 입력 영상을 만들어서 실제로 구해보면 공부하는데 도움이 되지 않을까 싶다 예를들면 x축 또는 y축에 대칭인 영상이나, 대각선으로 대칭인 영상을 이용하여 GLCM을 구하고 Correlation 특징을 구하는 것도 좋은 공부의 기회가 될 것 같다.

6. Inverse Difference Moment(IDM) Feature

 

이름에서도 알수 있듯이 차이와 반대되는 특징? ㅋㅋ 너무 직역했나, 여튼 의미적으로 이름과 수식이 참 잘 맞아떨어지는 것 같다. i-j의 제곱이 분모에 들어가 있는데, 결국 이것의 의미는 입력 영상의 픽셀 차이가 크다면 특징에 작은 가중치가 적용되고 차이가 작다면, 즉 GLCM의 대각 성분의 경우라면 큰 가중치가 적용되어 특징을 생성한다는 것을 알 수 있다. 의미적으로 보면 특징 이름과 매우 유사하게 입력 영상에서 차이, 즉 대비가 큰 픽셀의 빈도가 많다면 값이 작아지고 대비가 작은 픽셀의 빈도가 많으면 이 특징 값은 커질 것임을 예측할 수 있다.

7. Sum Average Feature

 

이것도 수식 표현이 어려운것 같다. 쉬운 이해를 위해 아래 예제를 보면 쉽게 이해할 듯 하다.

 

위의 예제와 같이 대각선에 있는 GLCM의 합을 더해서 i와 곱하는 형태이다. 수식을 풀어보자면, i는 GLCM에서 x축과 y축의 인덱스의 합이다. 예를들어 (0,3) 이라면 i=0+3 = 3 인 것이다. 따라서 위쪽에 예제로 든 GLCM은 4*4 형태이다 그렇다면 i의 범위는 2(n-1) 이므로 2(4-1) = 6 인 것이다. 여튼 이건 정말 수식이 난해하게 표현된 듯 하다.

8. Sum Variance Feature

 

 

 

9. Sum Entropy Feature

 

 

 

10. Difference Variance Feature 

 

 

F'10 은 아래의 Difference Average로 계산하면 되며, 이것을 이용하여 F10 Difference Variance 를 구한다.

 

11. Difference Entropy Feature

 

 

12. Information Measures of Correlation Feature 1

13. Information Measures of Correlation Feature 2

아래 수식들은 위의 12번과 13번을 구하기 위한 서브 수식을 보여준다.

 

12번 특징과 13번 특징을 구하는 예제이다.  

 

 

 

 

GLCM 구하는 코드

void CoOccurrenceMatrix::CalculateCoocurrenceMatrix ( void )
{
	int x, y;
	int d=1; // distance
	for (y = 0; y < (int ) m_BMPfile -> m_iImageHeight ; y++ )
	{
		for(x= 0; x <(int ) m_BMPfile -> m_iImageWidth ; x++ )
		{
			UpdatePixel ( x, y, x-d, y-d );
			UpdatePixel ( x, y, x, y-d );
			UpdatePixel ( x, y, x+d, y-d );
			UpdatePixel ( x, y, x-d, y );
			UpdatePixel ( x, y, x+d, y );
			UpdatePixel ( x, y, x-d, y+d );
			UpdatePixel ( x, y, x, y+d );
			UpdatePixel ( x, y, x+d, y+d );
		}
	}
	// normalization
	for (int i=0;i<Ng;i++)
		for(int j=0;j<Ng;j++)
			m_COM [i][j] = m_COM [i][j]/ normal ;
}

void inline CoOccurrenceMatrix::UpdatePixel ( int x1 , int y1 , int x2 , int y2)
{
	// Make sure the neighbour pixel exists ( can be e.g. negative ):
	if( x2 < 0 || x2 >= ( int) m_BMPfile -> m_iImageWidth
		|| y2 < 0|| y2 >= ( int) m_BMPfile -> m_iImageHeight )
		return ;
	
	unsigned int pixel , neighbour ;
	pixel = m_BMPfile -> m_ImageData [y1* m_BMPfile -> m_iImageWidth + x1 ];
	neighbour = m_BMPfile -> m_ImageData [y2* m_BMPfile -> m_iImageWidth +x2 ];
	m_COM [ pixel ][ neighbour ] ++ ;

}

13개의 특징을 빠르게 구하는 코드( 3번의 For 루틴으로 모두 구함)

//////////////////////////////////////////////////////////////////////////
//First 
for ( i = 0; i < Ng; i++ )
{
	for( j = 0; j < Ng; j++ )
	{
		m_fEnergy += m_COM [i][j]* m_COM [i][j]; // energy
		m_fContrast += ((i-j)*(i-j))* m_COM [i][j]; // contrast
		m_fHomogeneity += m_COM [i][j ]/(1+(i-j)*(i-j)); // homogeneity
		if( m_COM [i][j] != 0 )
			m_fEntropy += m_COM [i][j]* log( m_COM [i][j]) ;// entropy
		ux= ux + i* m_COM [i][j]; // mean of p_x
		uy= uy + j* m_COM [i][j]; // mean of p_y
		pxy [i+j] += m_COM [i][j]; // p_{x+y}(k)
		if (i >=j) pdxy [i-j] += m_COM [i][j]; // p_{x-y}(k)
		else pdxy [j-i] += m_COM [i][j];
		px[i] += m_COM [i][j]; // p_x
	}
}
//////////////////////////////////////////////////////////////////////////
//Second
for( i = 0; i < Ng; i++ )
{
	for( j = 0; j < Ng; j++ )
	{
		stdDevix += ((i-ux)*(i-ux)* m_COM [i][j]);
		stdDeviy += ((j-uy)*(j-uy)* m_COM [i][j]);
		m_fVariance += (i-ux)*(i-ux)* m_COM [i][j]; // variance
	}

	if( px[i] != 0 ) hx = hx + px[i]* log(px[i]);
	
	DiffAver = DiffAver + i* pdxy [i];
	m_fSumAver += (2*i)*pxy [2*i] ; // sum average
	m_fSumAver += (2*i+1)*pxy [2*i+1] ;
	
	if( pxy [2*i ]!= 0 )
		m_fSumEntr += pxy [2*i]* log(pxy [2*i]); // sum entropy
	if( pxy [2*i +1]!= 0 )
		m_fSumEntr += pxy [2*i +1]* log(pxy [2*i +1]) ;
	if( pdxy [i] != 0 )
		m_fDiffEntr += pdxy [i]* log( pdxy [i]); // different entropy
}

//////////////////////////////////////////////////////////////////////////
//Third
for( i = 0; i < Ng; i++ )
{
	for ( j = 0; j < Ng; j++) 
	{
		m_fCorrelation +=((i-ux)*(j-ux)* m_COM [i][j]) /( stdDevix * stdDevix );
		if (( px[i ]!= 0) &&( px[j ]!= 0))
		{
			hxy1 += m_COM [i][j]* log(px[i]* px[j]);
			hxy2 += px[i]* px[j]* log(px[i]* px[j]);
		}
	}
	m_fSumVari +=(2*i- m_fSumAver ) *(2*i- m_fSumAver )*pxy [2*i];
	m_fSumVari +=(2* i+1- m_fSumAver ) *(2* i+1- m_fSumAver )*pxy [2*i +1];
	m_fDiffVari +=(i- DiffAver )*(i- DiffAver )* pdxy [i ];// different variance
}
hxy1 = -hxy1 ;
hxy2 = -hxy2 ;
m_fInfMeaCor1 = ( m_fEntropy - hxy1 )/hx; // InfMeaCor1
m_fInfMeaCor2 = sqrt (1- exp ( -2*( hxy2 - m_fEntropy ))); // InfMeaCor2