[영상신호처리] Filtering 알고리즘 공부

션·2022년 2월 5일
영상신호처리
0

영상신호처리

목록 보기
1/5

2021 - 1 영상신호처리 수업 과제 수행

  • FFT 를 이용한 노이즈 분석 및 소거과정
  • Filtering 에 의한 영상 복원 (mean filter)
  • Filtering에 의한 영상 복원 (median filter)

1. 노이즈 제거를 위한 전처리 과정 간단 소개

  • 어떤 이미지의 노이즈 발견을 위해, 그 이미지의 스펙트럼을 관찰하는 것은 필수적인 과정이다. 그래서 나는 먼저 2D-FFT 함수를 이미지에 적용해 각 픽셀에서의 스펙트럼을 구해냈다.

  • 스펙트럼 이미지의 시각화를 위해, 스펙트럼의 Real part와 Imag part 값을 제곱해 square root를 적용시켜주고, log scale로 scaling 시키는 Compute_Spectrum 함수까지 작성했다.

  • 하지만 아무리 log scale로 Specturm 값을 scale 했더라도, 실제 이미지 출력에 쓰이는 0~255까지의 값보다는 훨씬 큰 값이기 때문에, 이를 scale 하기 위한 Dnormalize 함수를 이용하였다. 스펙트럼 이미지 내에서 Max 값과 Min 값을 찾은다음 그를 기준으로 모든 픽셀의 값을 0~255 사이로 스케일 해주는 함수이다.

  • 또한 MoveCenter라는 중심이동을 시켜주는 별도의 함수를 만들어 사용하였다. 1사분면과 4사분면의 이미지를 바꾸고 2 사분면과 3사분면의 이미지를 바꿔주는 함수인데, 이렇게 하면 중심에 가장 낮은 스펙트럼 값이 위치하게 되므로, 스펙트럼 이미지의 규칙성을 쉽게 발견할 수 있다는 장점이 있었다. 

void MoveCenter(BYTE **img1, BYTE **img2, int width, int height) {
  int x, y;
  int cx = (int)( width * 0.5 );
  int cy = (int)( height * 0.5);

  for (y=0; y<cy; y++) {
  	for (x=0; x<cx; x++) {
  		img2[y][x] = img1[cy + y][cx + x];
  		img2[cy + y][cx + x] = img1[y][x];

  		img2[y][cx+x] = img1[cy+y][x];
  		img2[cy+y][x] = img1[y][cx+x];
		}
	}
}
  • 내가 작성한 MoveCenter 함수의 특징
    Real2, Imag2 의 스펙트럼 이미지를 보기위해 가공된 img1의 이미지를, cx, cy라는 영상의 중암점을 기준으로 1, 4 사분면 끼리 바꾸고 2, 3 사분면끼리 바꾸어 img2에 넣은 형태이다. 고로

스펙트럼 Real2, Imag2 에 바로 대응하는 이미지 = img1 변수
img1을 보기 좋은 형태로 가공한 이미지 = img2 변수

이다.

  • 스펙트럼 영역에서 노이즈 제거 후 이미지를 복구할 때는 2D-FFT의 역과정인 2D-IFFT 함수를 작성해서 사용하였다. 이 함수를 사용하는 과정에서 많은 오류가 발생했었는데, 해결하는 과정은 뒷부분의 disscusion에서 좀 더 자세히 다뤄보고자 한다.

2. Filtering에 의한 영상복원 1 – MoonN.256

  • 위와 같은 전처리 과정을 마치고 난 후, MoonN.256 이미지에 어떤 형태의 노이즈가 있나 확인해 볼 수 있었다. MoonN.256 이미지는 중심 이동된 스펙트럼 상에서 특정한 원을 그리는 노이즈들이 존재하는 형태였다. 이를 제거하기 위해서는 먼저 1)중심에서부터 노이즈까지의 평균적인 거리를 계산하고, 그를 통해 2) 노이즈에서부터 걸러내고 싶은 Band width를 계산해야겠다는 계획을 세웠다.

다음과 같이 노이즈가 스펙트럼 상에서 원형으로 존재하는 것을 알 수 있으며, 영상 정보를 통해 알아본 대각선 방향 노이즈 까지의 거리는 대략 원점으로부터 102.55 정도 였다.

또한, 중심으로부터 가로, 세로 방향에 놓인 노이즈 까지의 거리는 97 이었기에, 102.55와 97의 중간값은 대략 99.415이다.

그래서 Noiselen = 99.415 로 설정하고, Widthlen = 5 정도로 설정한다면 충분히 노이즈를 걸러낼 수 있을 것이다.

  • 문제 해결을 위해 작성한 subroutine program
void Filter1(double **Real2, double **Imag2, int width, int height, double Noiselen,double Widthlen)
{

int cx = (int)(width * 0.5); // 원점에서부터 중간 x 지점까지를 나타내는 변수
int cy = (int)(height * 0.5); // 원점에서부터 중간 y 지점까지를 나타내는 변수

int x, y; // 이미지의 픽셀을 조회하기 위한 변수

for(y = 0; y< cy; y++) {
  for(x = 0; x< cx; x++) { 
  // for 문에서 x=0부터 x<cx 까지로, y= 0부터 y < cy 까지로 지정해, 1 사분면에 대해서만 노이즈와의 거리를 계산하기로 결정했다.

	double d = sqrt((double)(x*x + y*y));
  // d라는 변수는 원점 (0,0) 에서부터, 현재 픽셀 (x, y) 까지의 거리로 설정했다.

	if(abs(d - Noiselen)<=Widthlen) { 
  // d - Noiselen 의 절댓값이 제거하고 싶은 bandwidth의 폭 범위 내에 존재할 경우에는, 현재 픽셀이 노이즈 범위 내에 있는 픽셀이라고 판단할 것이다. 이 때 현재 픽셀의 Real 값, imag 성분을 전부 0으로 만들도록 코드를 작성했다.

// 그런데 현재 for 문이 1사분면의 픽셀에 대해서만 돌고 있으므로, 픽셀의 대칭 특성을 이용해 2,3,4 분면의 지점에서도 값을 제거하는 과정을 진행해 줄 것이다.

  //if문 조건 만족시 제 1사분면에 대해 스펙트럼 성분을 제거해주는 코드
  Real2[y][x] = 0;
  Imag2[y][x] = 0;

  //if문 조건 만족시 제 2사분면에 대해 스펙트럼 성분을 제거해주는 코드
  Real2[y][width - x - 1] = 0;
  Imag2[y][width - x - 1] = 0;


  //if문 조건 만족시 제 3사분면에 대해 스펙트럼 성분을 제거해주는 코드
  Real2[height - y - 1][x] = 0;
  Imag2[height - y - 1][x] = 0;


  //if문 조건 만족시 제 4사분면에 대해 스펙트럼 성분을 제거해주는 코드
  Real2[height - y - 1][width - x - 1] = 0;
  Imag2[height - y - 1][width - x - 1] = 0;

	}
}
}
}
  • 결과 구현 모습

다음과 같이, ring 의 형태로 band reject filtering이 된 것을 확인할 수 있다.

  • discussion
    IFFT-2D 함수 구현시의 오류: Real part와 imag part의 중요성
    초반에 IFFT 함수를 적용만 하면, 자꾸 이미지가 반전되어서 출력되는 문제점을 발견하였다.

원인을 탐색하다보니 2D – IFFT 함수 구현을 잘못했기 때문이라는 것을 알게 되었다. ifft 함수를 적용해야할 Imag 파트 부분에, 계속 Real part 값만 집어 넣고 있었다. 

[처음 잘못 작성한 IFFT_2D의 일부코드]

for (y=0; y<height; y++)
{
for(x=0; x<width; x++)
{

tmp[x].real = Real[y][x];
tmp[x].imag = Imag[y][x];
}
ifft(tmp,LN);
for(x=0; x<width; x++)
{
Real2[y][x] = tmp[x].real;
Imag2[y][x] = tmp[x].real;
}
}
for(x=0; x<width; x++)
{
for(y=0; y<height; y++)
{

tmp[y].real = Real2[y][x];
tmp[y].imag = Imag2[y][x];
}

ifft(tmp, LN);

for (y=0; y<height; y++){
Real2[y][x] = tmp[y].real;
Imag2[y][x] = tmp[y].real;

}

}

Imag part에는 이미지의 위상정보가 들어와야 하는데, 그러한 정보가 들어오지 못하니 이미지가 반전된 형태로 출력하는 것이라고 추측할 수 있었다.

이미지의 스펙트럼 중, imaginary 성분이라는 것이 무엇인지 정확히 와닿지 않았는데, 이번 실수를 계기로 이미지의 위상성분은 이미지의 회전에 관여한다는 사실을 확인할 수 있어 유익했다. 다음에 코드를 작성할 때도, 이런 위상정보를 누락하는 실수를 저지르지는 않았는지 면밀히 확인해야겠다는 결심도 하게 되었다.

  • for 문 간소화를 통한 계산 속도 절감
    제거할 노이즈를 선택하는 과정에서 for문을 썼었는데, 처음에는 이미지의 모든 픽셀에 대해 for문을 적용해야한다는 생각을 했었다. 
double d = sqrt((double)(x*x + y*y));

for (y=0; y<cy; y++)
{
	for (x=0; x<cx; x++)
{
	if(abs(d - Noiselen)<=Widthlen)
{
	Real2[y][x] = 0;
	Imag2[y][x] = 0;
}
}
}
d = sqrt((double)((width-x-1)*(width-x-1) + y*y));
for (y=0; y<cy; y++)
{
	for (x=cx; x<width; x++)
{
	if(abs(d - Noiselen)<=Widthlen)
{
	Real2[y][x] = 0;
	Imag2[y][x] = 0;
}
}
d = sqrt((double)(x*x + (height - y - 1)*(height - y - 1)));
for (y=cy; y<height; y++)
{
	for (x=0; x<cx; x++)
{
	if(abs(d - Noiselen)<=Widthlen)
}
	Real2[y][x] = 0;
	Imag2[y][x] = 0;
}
}
d = sqrt((double)((width-x-1)*(width-x-1) + (height - y - 1)*(height - y - 1)));
for (y=cy; y<height; y++)
{
	for (x=cx; x<width; x++)
{
	if(abs(d - Noiselen)<=Widthlen)
{
	Real2[y][x] = 0;
	Imag2[y][x] = 0;
}
}
// 이런 방법으로, 총 4번의 for문을 돌렸다
// 개선한 코드의 모습, for 문을 1번만 돌려서 처리할 수 있음을 확인할 수 있다.

for(y = 0; y< cy; y++) {
for(x = 0; x< cx; x++) { 
Real2[y][x] = 0;
Imag2[y][x] = 0;

Real2[y][width - x - 1] = 0;
Imag2[y][width - x - 1] = 0;

Real2[height - y - 1][x] = 0;
Imag2[height - y - 1][x] = 0;

Real2[height - y - 1][width - x - 1] = 0;
Imag2[height - y - 1][width - x - 1] = 0;
}
}

노이즈가 중심으로부터 같은 거리에 떨어져 있다는 특성을 이용하면, 모든 픽셀을 조회할 필요가 없었다. (0,0) 에서 영상의 중심까지의 픽셀 (즉, 1사분면) 을 조회한 다음 노이즈를 걸러주고, 점대칭의 특성을 이용해 나머지 2,3,4 사분면에 대해서 노이즈를 제거해주면, for문의 사용을 대폭 줄일 수 있고, 연산속도도 절감할 수 있다.

FFT에 비해 DFT 과정에서 시간이 오래 소모된 이유를 생각해보면, 중첩된 for loop문을 오랜 시간 돌기 때문이었다. 고로 이미지 처리 속도를 높이기 위해서는 for loop 문의 과정을 최대한 간소화 해야한다. 위의 코드를 짜는 과정에서도 마찬가지다. 전체 픽셀을 조회하는 대신 점대칭의 특성을 이용해 일부 픽셀만 조회하는 방법으로 for문을 대폭 간소화할 수 있었고, 이를 통해 계산속도도 절감할 수 있다는 결과를 얻었다. 다음 과제에서도 최대한 for문의 사용을 최소화 하는 방향으로 코드를 짜야겠다는 결심을 하게 되었다.

2. Filtering에 의한 영상복원 2 – HouseN.256

  • 과정
    HouseN 이미지는 위와 다르게, 노이즈가 중심을 기점으로 복소 대칭으로 등장하는 형태였다. 고로 1사분면의 노이즈를 제거하면 4사분면의 노이즈도 처리가 가능하며, 2사분면의 노이즈를 제거하면 3사분면의 노이즈까지 제거할 수 있다. 이러한 특성을 파악하고 다음과 같은 Filter 2 함수의 프로토타입을 작성해보았다.
void Filter2(double **Real2, double **Imag2, int width, int height, int quadrant, int Twidth, int Theight)

몇가지 새로운 argument 를 소개 하자면,
∎ int quadrant : 노이즈가 위치한 사분면을 지시하는 변수이다. quadrant = 1 이면 1, 4 사분면에 위치한 노이즈를 제거하고, quadrant = 2 이면 2, 3 사분면에 위치한 노이즈를 제거한다.
∎ int Twidth : Target width의 줄임말로, 중심으로부터 노이즈까지 떨어진 x 좌표 픽셀값을 의미한다.
∎ int Theight : Target Height의 줄임말로, 중심으로부터 노이즈까지 떨어진 y 좌표 픽셀값을 의미한다.

위와 같은 방법을 이용한다면, 1)어떤 노이즈가 어떤 사분면에 있는지 알고 2) 중심으로부터 얼마 만큼의 좌표에 떨어져 있는지만 안다면 복소 대칭의 형태로 등장하는 모든 노이즈를 제거할 수 있다는 이점이 있다.

  • 문제 해결을 위해 작성한 subroutine program
void Filter2(double **Real2, double **Imag2, int width, int height, int quadrant, int Twidth, int Theight)
{
int cx = (int)(width * 0.5);
int cy = (int)(height * 0.5);

int x, y;

//아래의 코드는 1,4 사분면에 대한 노이즈를 처리하는 과정이다.
if(quadrant == 1)
{
for(y = Theight-5; y< Theight + 6; y++)
{ 
for(x = Twidth - 5; x< Twidth + 6; x++)
{ 

// Noise Target 좌표를 기점으로, ±5 픽셀만큼 떨어진 부분들의 스펙트럼 성분들을 0으로 바꿔주었다.
Real2[y][x] = 0;
Imag2[y][x] = 0;

// 그 후 1사분면에 대칭인 4사분면에 대해서도 노이즈를 제거해주었다.
Real2[height - y - 1][width - x - 1] = 0;
Imag2[height - y - 1][width - x - 1] = 0;

}

}
}

//아래의 코드는 2, 3 사분면에 대한 노이즈를 처리하는 과정이다.
else if(quadrant == 2)
{
for(y = Theight-5; y < Theight + 6; y++)
{
for(x = width - (Twidth + 5); x< width - (Twidth - 6); x++)
{
Real2[y][x] = 0;
Imag2[y][x] = 0;

Real2[height - y - 1][width - x - 1] = 0;
Imag2[height - y - 1][width - x - 1] = 0;

}
}
}
}

영상 정보를 통해, HouseN.256 이미지의 1,4사분면 노이즈는 중심으로부터 56x24 픽셀만큼 떨어져 있음을 알 수 있었고, 2,3 사분면 노이즈는 중심으로부터 38x12 픽셀만큼 떨어져 있음을 알 수 있었다. 이 값들을 Filter2의 Twidth, Theight 인자로 사용할 것이다.

-결과 구현 모습

다음과 같이 점대칭 형태의 노이즈가 잘 제거된 것을 확인할 수 있다.

3. Mean-filtering과 Median-filtering – ApeN.256

  • Mean-filtering 과정
    Mean filter는 NxN window를 이미지에 적용시켜, window 내에 들어온 픽셀값들을 모두 평균하여 바꿔주는 필터링 기법이다. (본 과제에서는 3x3 window를 이용했다.) 이를 구현하기 위해서 원본 이미지를 img1에 담고, img1의 평균 픽셀 값을 img2 라는 새로운 변수에 담는 코드를 작성하기로 결정했다.

  • 문제 해결을 위해 작성한 subroutine program

void MeanFilter(BYTE **img1, BYTE **img2, int width, int height)
{
int x,y,u,v;
//x, y 는 img1의 전체 픽셀을 조회하기 위한 변수들이고, u,v는 3x3 window를 조회하기 위한 변수들이다.

int sum;
//sum 이라는 변수는 3x3 픽셀의 값을 전부 더한 후 평균하기 위한 목적의 변수이다.

sum = 0;


for(y = 1; y < height - 1; y++){
for(x = 1; x<width -1; x++)
{ // 위와 같이, 3x3 window를 움직이기 위한 x 픽셀의 범위는 1부터 width-2 까지이고, y 픽셀의 범위는 1 부터 height-2 까지가 되어야 한다.

for(v = y-1; v<y+2; v++)
{
for(u = x-1; u < x+2; u++)
{
//window는 x,y 좌표를 중심으로, (x-1 ~ x+1), (y-1 ~y+1) 까지 총 9개의 픽셀의 값을 조회한다.
sum += img1[v][u]; 
//그후 그 값들을 sum에 더한다.

}

}

sum = (int)(sum / 9);
//sum에 더해진 9개의 값을 평균 내준 다음, 정수형 변환을 해준다. 

for(v = y-1; v<y+2; v++)
{
for(u = x-1; u < x+2; u++) {
img2[v][u] = sum; 
//그 후 img2의 3x3 만큼의 픽셀들에 img1의 3x3 픽셀들의 평균값으로 채워서 필터링해준다.

}
}
sum = 0;
//sum은 반복문 내에서 다시 사용되어야하기 때문에, 한차례의 반복문이 끝난 후에는 꼭 그 값을 0으로 초기화한다.
}
}

}

  • 결과 구현 모습

  • discussion
    mean filter의 구현은 픽셀들을 평균값 내는 것이기 때문에 구현 자체는 어렵지 않았다. 하지만 한 가지 궁금한 부분이 있었다. 바로 window의 이동 간격을 다르게 움직이면 어떻게 이미지가 달라지는가 였다. 위의 이미지는 window를 이동하는 간격을 1로 한 것인데, 3으로 이동 간격을 움직이면 어떻게 될지 궁금하여 한번 실행해보았다.

window를 움직이는 간격을 3으로 키워보니, 간격이 1일 때보다는 resolution이 좋지 않다는 것을 알 수 있다. 그 대신 이미지 처리에 드는 시간은 훨씬 적어졌음을 확인할 수 있다. window 이동크기가 커지면서, 연산과정도 당연히 더 적어지기 때문이다. 이를 통해 높은 노이즈 제거성능보다 빠른 노이즈 제거 속도를 요하는 작업에선, window 이동폭의 크기를 크게 하여 필터링을 적용해도 좋을 것 같다는 생각이 들었다.

  • Median-filtering 과정
    Median filter는 NxN window를 이미지에 적용시켜, window 내에 들어온 픽셀값들을 정렬해 그 중간값을 3x3 의 픽셀값으로 채택하는 필터이다. 이를 구현하기 위해서 1)원본 이미지를 img1에 담고, 2) img1의 3x3 픽셀들의 값을 정렬해 중간값을 찾고, 3) 그 값을 img2 라는 새로운 변수에 담는 코드를 작성하기로 결정했다.

  • 문제 해결을 위해 작성한 subroutine program

void MedianFilter(BYTE **img1, BYTE **img2, int width, int height)
{
BYTE *tmp = (BYTE *)malloc(sizeof(BYTE) * 9);
//window에 들어온 픽셀을 정렬하기 위해, tmp라는 1x9의 배열을 동적할당을 이용해 만들어 주었다.

int x,y,u,v; 
// x,y는 img1의 픽셀을 조회하기 위한 변수, u,v는 3x3 window의 픽셀을 조회하기 위한 변수이다.

int i, k; // i와 k는 tmp 배열에 버블 정렬 기법을 이용하기 위한 변수들이다. 
int temp, d; // d는 tmp 배열의 index를 지정하기 위한 변수이고( tmp[d] 처럼 ), temp는 tmp[d]에 저장된 값을 담기 위한 임시 변수이다.

d=0; //tmp 배열의 시작 index는 0 이므로, d = 0으로 초기화한다.

for(y = 1; y < height - 1; y++) 
{
for(x = 1; x<width -1; x++) 
{
for(v = y-1; v<y+2; v++) 
{
for(u = x-1; u < x+2; u++) 
{
tmp[d] = img1[v][u]; 
// tmp 배열에 3x3 window의 픽셀값을 대입해준다. 
d++;
// 한 개의 값을 대입한 후에는 다음 픽셀을 위해 index d의 값을 올려준다.

}
}

d=0; // tmp 배열에 데이터가 모두 담긴 상태이므로, 다음 반복문을 위해 0으로 초기화해준다.

//tmp 배열의 버블 정렬을 위한 코드
for (k=0; k<9; k++)
{
for (i=0; i<8; i++)
{
if(tmp[i] > tmp[i+1])
{
temp = tmp[i];
tmp[i] = tmp[i+1];
tmp[i+1] = temp;
//배열의 첫번째 요소부터 마지막 요소까지에 대해 적용한다. 만약 tmp[i] > tmp[i+1] 이라면 두 위치의 값을 swap 하는 코드이다.
}

}

}

//위의 버블정렬이 끝나고 나면 tmp는 작은값부터 차례로 정렬된 형태일 것이다. 이 행렬의 중간값인 tmp[4]를 취한다.
for(v = y-1; v<y+2; v++) 
{
for(u = x-1; u < x+2; u++) 
{
img2[v][u] = tmp[4]; 
//img2의 3x3 영역에 대해, tmp의 중간값인 tmp[4] 로 전부 채워준다.

}
}
}
}

free(tmp);
//tmp는 동적할당된 공간이므로, 꼭 해제를 시켜준다.
}
  • 결과 구현 모습


  • discussion
    위 코드에서, 픽셀의 중간 값을 찾기 위해 픽셀값을 정렬하는 기법으로 버블정렬을 사용하였다. 버블 정렬이란 서로 인접한 두 원소를 검사하여 정렬하는 알고리즘으로, 인접한 두 원소를 비교해 크기가 순서대로 되어있지 않으면 교환하는 방식이다. 코드상으로 구현하기 간단하다는 장점은 있으나 모든 원소끼리 비교해야 하므로 시간이 오래 소요된다는 단점이 있었다. 그래서 mean filter보다 median filter가 연산 속도가 더 많이 소요된다는 사실 또한 확인할 수 있었다.

    ApeN2.256 이미지에 각각 mean filter와 median filter의 적용을 통해, salt 형태의 노이즈가 존재할 때는 mean filter 보다 median filter가 더 유용하다는 것을 확인할 수 있었다.
profile
션
._.
다음 포스트

[영상신호처리] Disparity Map 개념 이해

0개의 댓글