[C++/MFC] Median Filter

Median Filter?

중간값 필터링이란, 이미지 픽셀의 중간값을 찾아 필터링하는 과정이다.

출처: https://thebook.io/006796/0302/

이런 이미지 만으로는 제대로 된 분석을 하기 어렵기 때문에
이미지 가공이 필요하다.

대표적인 소금-후추 노이즈의 경우 이미지에 노이즈가 들어갔을 때 명확한 영상처리를 위한 필터링이 필요하다. 이에 따라 여러 알고리즘을 활용해 이미지의 중간값을 찾을 수 있다.

---

Median Filter의 처리 방식

출처: https://www.researchgate.net/figure/A-graphical-depiction-of-the-median-filter-operation_fig1_280925268

Median Filter는 이미지의 픽셀을 filter의 Size에 맞게 등분하여, 0~255 사이의 값들의 배열을 정렬해 중간값을 채택하는 방식이다. 이와같은 방식으로 튀는 픽셀값(노이즈)를 잡을 수 있다.

예를들어 size가 3이라면 3x3 사이즈의 마스크로 중간값을 추출할 것이다.

출처: https://ietresearch.onlinelibrary.wiley.com/doi/10.1049/iet-ipr.2016.0737

튀는 픽셀인 255를 처리하는 과정

출처: https://preventionyun.tistory.com/33

히스토그램을 이용해 구현해보려하는데 그냥 sort만으로도 구현가능하다.

---

1. 이미지 불러오기

// ImageDlg.cpp : 구현 파일
//

#include "stdafx.h"
#include "Image.h"
#include "ImageDlg.h"
#include "afxdialogex.h"
#include <vector>
#include <algorithm>
#include <math.h>

void CImageDlg::OnBnClickedButtonOpen()
{
	// 파일 열어서 bitmap으로 출력
	CFileDialog fDlg(TRUE);
	fDlg.DoModal();
	CFile Rfile;

	if (!Rfile.Open(fDlg.GetPathName(), CFile::modeRead))
	{
		MessageBox(_T("Can't Open Image File!"));
		return;
	}

	memset(m_original, 0, IMAGE_SIZE);

	UINT FileLength = (UINT)Rfile.GetLength();
	unsigned char*ps = new unsigned char[FileLength];

	Rfile.Read(ps, FileLength);
	for (int y = 0; y < IMAGE_HEIGHT; y++)
	{
		for (int x = 0; x < IMAGE_WIDTH; x++)
		{
			m_original[y * IMAGE_WIDTH + x] = ps[y * IMAGE_WIDTH + x];
		}
	}

	delete[] ps;

	if(bmpInfo != nullptr)
	{
		delete bmpInfo;
		bmpInfo = nullptr;
	}

	int rwsize = (((512) + 31) / 32 * 4);
	bmpInfo = (BITMAPINFO*)malloc(sizeof(BITMAPINFO) + IMAGE_HEIGHT * sizeof(RGBQUAD));
	bmpInfo->bmiHeader.biBitCount = 8;
	bmpInfo->bmiHeader.biClrImportant = IMAGE_WIDTH;
	bmpInfo->bmiHeader.biClrUsed = IMAGE_WIDTH;
	bmpInfo->bmiHeader.biCompression = 0;
	bmpInfo->bmiHeader.biPlanes = 1;
	bmpInfo->bmiHeader.biSize = 40;
	bmpInfo->bmiHeader.biSizeImage = rwsize * IMAGE_WIDTH;
	bmpInfo->bmiHeader.biHeight = -IMAGE_HEIGHT;
	bmpInfo->bmiHeader.biWidth = IMAGE_WIDTH;
	bmpInfo->bmiHeader.biXPelsPerMeter = 0;
	bmpInfo->bmiHeader.biYPelsPerMeter = 0;

	for (int i = 0; i < 512; i++) {
		bmpInfo->bmiColors[i].rgbRed = bmpInfo->bmiColors[i].rgbGreen
			= bmpInfo->bmiColors[i].rgbBlue = i;
		bmpInfo->bmiColors[i].rgbReserved = 0;
	}

	StretchDIBits(m_memDC1->GetSafeHdc(),
		0, 0, m_rtViewer1.Width(), m_rtViewer1.Height(), 0, 0, IMAGE_WIDTH, IMAGE_HEIGHT,
		m_original, bmpInfo, DIB_RGB_COLORS, SRCCOPY);

	InvalidateRect(m_rtViewer1, FALSE);
}

2. 메디안 필터로 이미지 처리하기

void CImageDlg::OnBnClickedButtonMedian()
{
	// 이미지가 로드되지 않았을 경우, return
	if (m_original == nullptr)
		return;

	// 결과를 저장할 임시 이미지 생성
	BYTE* temp = new BYTE[IMAGE_SIZE];
	memcpy(temp, m_original, IMAGE_SIZE);

	// 필터 크기 설정
	int filterSize = 7; //필터의 크기는 홀수여야함
	int pad = filterSize / 2; 

	// 히스토그램과 중간값을 저장할 변수 생성
	
	int histogram[256] = { 0, }; // 배열의 모두를 모두 0으로 초기화
	int median = 0;
	/*히스토그램에서는 각 숫자의 빈도수를 나타냄. 이미지의 픽셀이 0~255이므로 256개의 빈도수를 가지게 됨.
	따라서 전체 픽셀 수를 가로 * 세로로 구함. 이후 각 히스토그램의 빈도수를 처음부터 누적해서 더하고, 누적 빈도가 전체 픽셀의
	절반이 넘으면 그 픽셀이 중간값이 됨.*/

	// 모든 픽셀에 대해
	for (int y = pad; y < IMAGE_HEIGHT - pad; y++) { // pad는 필터의 절반 크기로, 이미지의 가장자리를 제외한 부분을 순회함
		for (int x = pad; x < IMAGE_WIDTH - pad; x++) {

			// 히스토그램 초기화
			memset(histogram, 0, sizeof(histogram));

			// Window 내의 값들을 히스토그램에 저장
			for (int i = -pad; i <= pad; i++) { // -pad ~ pad 범위만큼 ( -1 ~ 1 처럼 )
				for (int j = -pad; j <= pad; j++) {
					histogram[m_original[(y + i) * IMAGE_WIDTH + (x + j)]]++;
				}
			}

			// 중간값 계산
			int count = 0;
			for (median = 0; median < 256; median++) {
				count += histogram[median];
				if (count > filterSize * filterSize / 2) break; //count가 전체 픽셀의 절반을 넘음
			}

			// 중간값을 임시 이미지에 저장
			temp[y * IMAGE_WIDTH + x] = median;
		}
	}

	// 결과를 m_filtered에 복사
	memcpy(m_filtered, temp, IMAGE_SIZE);
	delete[] temp;

	if (m_memDC2 == NULL)
	{
		m_memDC2 = new CDC;
		m_memBit2 = new CBitmap;

		CClientDC dc(this); 
		m_memDC2->CreateCompatibleDC(&dc);
		m_memBit2->CreateCompatibleBitmap(&dc, m_rtViewer2.Width(), m_rtViewer2.Height());
		m_memDC2->SelectObject(m_memBit2);
		m_memDC2->SetStretchBltMode(4);
	}

	CClientDC dc(this); 

	StretchDIBits(m_memDC2->GetSafeHdc(),
		0, 0, m_rtViewer2.Width(), m_rtViewer2.Height(), 0, 0, IMAGE_WIDTH, IMAGE_HEIGHT,
		m_filtered, bmpInfo, DIB_RGB_COLORS, SRCCOPY);

	InvalidateRect(m_rtViewer2, FALSE);
}

중간값을 찾아야하기 때문에 Filter Size는 홀수로 설정한다.
Filter의 크기가 클 수록 그림이 좀 더 뭉개지는 것을 볼 수 있다.

---

1) Filter Size = 3

2) Filter Size = 7

히스토그램 외에도 Bubble Sort나 Quick Sort 등 다양한 기존 알고리즘으로도 구현할 수 있지만 히스토그램은 O(n)이고 퀵 정렬은 O(nlogn)이기 때문에 속도가 느리다.

실제로 버블 소트로 구현했을 때는 프리징이 걸리고 10초는 지나야 이미지가 출력됐다.

---

히스토그램을 사용하지 않은 선택 정렬 알고리즘을 통한 메디안 필터는 다음과 같다.
주의사항으로는 필터의 크기가 변경되면 배열 windowValues[]의 값 또한 변경돼야한다는 것이다. 예를들어 7x7 사이즈면 49.. 이런 식으로

void CImageDlg::OnBnClickedButtonMedian()
{
	// 이미지가 로드되지 않았으면 반환
	if (m_original == nullptr)
		return;

	// 결과를 저장할 임시 이미지 생성
	BYTE* temp = new BYTE[IMAGE_SIZE];
	memcpy(temp, m_original, IMAGE_SIZE);

	// 필터 크기 설정
	int filterSize = 3; // 필터 크기는 홀수여야 함
	int pad = filterSize / 2;

	// 각 픽셀에 대해
	for (int y = pad; y < IMAGE_HEIGHT - pad; y++) {
		for (int x = pad; x < IMAGE_WIDTH - pad; x++) {

			// 윈도우 내의 픽셀 값을 저장할 배열 생성
			BYTE windowValues[9]; // 3x3 필터를 가정. 만약 filterSize가 바뀌면 이 부분도 맞춰서 바꿔야함. 예, 5x5 = [25]
			int index = 0;

			// 윈도우 내의 값을 배열에 저장
			for (int i = -pad; i <= pad; i++) {
				for (int j = -pad; j <= pad; j++) {
					windowValues[index++] = m_original[(y + i) * IMAGE_WIDTH + (x + j)];
				}
			}

			// STL에서 제공하는 Sort() 대신에 Selection Sort를 사용하여 중간값을 찾음
			for (int i = 0; i < filterSize * filterSize - 1; i++) {
				int minIndex = i;
				for (int j = i + 1; j < filterSize * filterSize; j++) {
					if (windowValues[j] < windowValues[minIndex]) {
						minIndex = j;
					}
				}
				// Swap
				BYTE tempValue = windowValues[i];
				windowValues[i] = windowValues[minIndex];
				windowValues[minIndex] = tempValue;
			}
			/* Quick 정렬을 사용했으면 다음과 같음. 하지만 선택 정렬이 O(n^2)로 더 빠름
			std::sort(window.begin(), window.end());*/

			// 중간값을 임시 이미지에 저장
			temp[y * IMAGE_WIDTH + x] = windowValues[filterSize * filterSize / 2];
		}
	}

	
			
	// 결과를 m_filtered에 복사
	memcpy(m_filtered, temp, IMAGE_SIZE);
	delete[] temp;

	if (m_memDC2 == NULL)
	{
		m_memDC2 = new CDC;
		m_memBit2 = new CBitmap;

		CClientDC dc(this);
		m_memDC2->CreateCompatibleDC(&dc);
		m_memBit2->CreateCompatibleBitmap(&dc, m_rtViewer2.Width(), m_rtViewer2.Height());
		m_memDC2->SelectObject(m_memBit2);
		m_memDC2->SetStretchBltMode(4);
	}

	CClientDC dc(this);

	StretchDIBits(m_memDC2->GetSafeHdc(),
		0, 0, m_rtViewer2.Width(), m_rtViewer2.Height(), 0, 0, IMAGE_WIDTH, IMAGE_HEIGHT,
		m_filtered, bmpInfo, DIB_RGB_COLORS, SRCCOPY);

	InvalidateRect(m_rtViewer2, FALSE);
}