디지털 영상 처리 - 노이즈만 있을 때 복원하기(공간 필터링) : 적응 필터 구현

 안녕하세요. 지난 포스팅의 디지털 영상 처리 - 주파수 도메인 필터링에 의한 노이즈 감소에 대해서 알아보았습니다. 그런데 이전에 공간 필터링에 대해서 구현할 때 평균 필터와 순서-통계 필터만 구현하고 적응 필터링은 구현하지 않았습니다. 오늘은 적응 필터링을 구현해보도록 하겠습니다. 적응 필터링의 전체 코드는 아래의 링크를 참고해주시길 바랍니다.

 

skawngus1111/DIP

 

1. 적응적, 지역적 노이즈 감소 필터(Adaptive, Local Noise Reduction Filter)

 

 우측 상단의 영상은 평균은 0이고 분산은 1000인 가우시안 노이즈에 오염된 영상을 보여주고 있습니다. 그리고 그 오른쪽 영상과 아래 영상은 각각 산술 평균 필터와 기하 평균 필터를 적용한 결과입니다. 그리고 왼쪽 아래 영상은 적응적, 지역적 노이즈 감소 필터를 적용한 결과입니다. 실제로 3가지 필터링을 비교해보았을 때 적응 필터를 적용했을 때 가장 좋은 퀄리티의 영상을 복원한 것을 볼 수 있습니다. 

 

function [FilteredImage] = AdaptiveLocalNoiseReductionFilter(image,windowSize)

if size(image, 3) == 3 
    image = rgb2gray(image);
end

kernel = ones([windowSize, windowSize], 'double'); % 커널 정의
n = floor(windowSize/2);

[height, width] = size(image); % 입력 영상의 높이와 너비 추출
FilteredImage = zeros([height, width]); % 입력 영상의 크기에 따라서 스무딩 영상 초기화
image = double(image);
PaddedImage = zeros([height+2*n, width+2*n]);
PaddedImage(n+1:height+n, n+1:width+n) = image;
LocalMean = zeros([height, width]);
LocalVar = zeros([height, width]);

for i=n+1:height+n
    for j=n+1:width+n
        LocalMean(i-n, j-n) = mean(PaddedImage(i-n:i+n, j-n:j+n), 'all');
        LocalVar(i-n, j-n) = var(PaddedImage(i-n:i+n, j-n:j+n), 0, 'all');
    end
end

nvar = sum(LocalVar(:))/(height * width);
LocalVar = max(LocalVar,nvar);   
temp = PaddedImage(n+1:height+n, n+1:width+n);

FilteredImage = temp - nvar./LocalVar .* (temp - LocalMean);
FilteredImage = uint8(FilteredImage);
end

 위에는 제가 구현한 코드입니다. 영상을 패딩하는 부분까지는 동일하지만 여기에서 지역적 영상의 평균과 분산을 계산합니다. 이를 각각 LocalMean과 LocalVar에 저장해줍니다. 그 다음에 전체 영상의 노이즈에 대한 분산을 추정하기 위해서 지역적 영상의 분산의 평균을 계산하여 이를 추정된 노이즈 분산으로 취급합니다. 이때 책에서 $\sigma^{2}_{\eta} > \sigma^{2}_{L}$을 방지하기 위해서 max 함수를 통해 해결합니다. 나머지는 책에서 본 공식에 그대로 대입하여 해결합니다. 

 

2. 적응적 중간값 필터(Adaptive Median Filter)

 다음으로 적응적 중간값 필터의 효과를 보도록 하겠습니다. 왼쪽 영상은 아주 심하게 소금-후추 노이즈에 의해 오염된 영상을 보여주고 있습니다. 그리고 가운데 영상은 $11 \times 11$ 크기의 중간값 필터를 적용한 결과이고 오른쪽 영상은 $S_{\text{max}} = 11$ 그리고 $S_{\text{init}} = 3$로 지정했을 때, 적응적 중간값 필터를 적용했을 때 결과를 보여주고 있습니다. 두 결과를 비교해보면 소금-후추 노이즈는 전체적으로 제거되었지만 중간값 필터를 적용했을 때는 영상 내의 많은 디테일이 손실되었습니다. 그에 반해 적응적 중간값 필터를 적용했을 때는 디테일이 상대적으로 많이 보존되는 것을 볼 수 있습니다. 

 

function [FilteredImage] = AdaptiveMedianFilter(image,StartWindowSize, MaxWindowSize)
if size(image, 3) == 3 
    image = rgb2gray(image);
end

n = floor(StartWindowSize/2);
[height, width] = size(image); % 입력 영상의 높이와 너비 추출
FilteredImage = zeros([height, width]); % 입력 영상의 크기에 따라서 스무딩 영상 초기화
image = double(image);
PaddedImage = zeros([height+2*n, width+2*n]);
PaddedImage(n+1:height+n, n+1:width+n) = image;

CurrentWindowSize = StartWindowSize;

for i=n+1:height+n
    for j=n+1:width+n
        % STEP A-1. Extract Local Image and sort by pixel value
%         LocalImage = PaddedImage(i-n:i+n, j-n:j+n);
        SortByValue = ExtractLocalImage(PaddedImage, CurrentWindowSize, i, j, height, width);
%         SortByValue = sort(reshape(LocalImage, [1, CurrentWindowSize^2]));
        
        % STEP A-2. Get min, max, median of Local Image
        zmin = SortByValue(1);
        zmax = SortByValue(end);
        zmed = median(SortByValue);
        zcur = PaddedImage(i, j);
        
        % STEP A-3. Calculate A1 and A2
        A1 = zmed - zmin;
        A2 = zmed - zmax;
        
        % STEP A-4. Check sign of A1 and A2
        if (A1 > 0) && (A2 < 0) % Enter to Level B
            % STEP B-1. Calculate B1 and B2
            B1 = zcur - zmin;
            B2 = zcur - zmax;
            
            % STEP B-2. Check sign of B1 and B2
            if (B1 > 0) && (B2 < 0)
                FilteredImage(i - n, j - n) = zcur;
            else
                FilteredImage(i - n, j - n) = zmed;
            end
        else % Increase the size of window and repeat Level A until size of window is larger than max size window
            CurrentWindowSize = CurrentWindowSize + 2;
            if CurrentWindowSize > MaxWindowSize
                FilteredImage(i - n, j - n) = zmed;
                CurrentWindowSize = StartWindowSize;
            end
        end
    end
end
FilteredImage = uint8(FilteredImage);

end

function LocalImage = ExtractLocalImage(image, WindowSize, i, j, height, width)
cnt = 0;
LocalImage = [];
for si=-floor(WindowSize/2) : floor(WindowSize/2)
    if (i + si) < 1 || (i + si) > height
        continue
    end
    for sj=-floor(WindowSize/2) : floor(WindowSize/2)
        if (j + sj) < 1 || (j + sj) > width
            continue
        end
        cnt = cnt + 1;
        LocalImage = [LocalImage, image(i + si, j + sj)];
    end
end
LocalImage = sort(LocalImage);
end

 제가 구현한 적응적 중간값 필터 함수는 크게 2가지로 이루어져있습니다. AdaptiveMedianFilter는 전체 적응적 중간값 필터를 적용하는 함수이고 ExtractLocalImage는 윈도우의 크기가 달라질 때마다 지역적 영상을 추출하는 함수입니다.