Non local means 알고리즘, MATLAB 구현 코드 포함

 

 

 

Non local means algorithm은 노이즈 제거에 강력한 성능을 보여주는 denoising 알고리즘이다.

일반적으로, gaussian smoothing과 같은 노이즈 제거 알고리즘을 많이 사용하는데 이러한 smoothing 방식은 local한 데이터들, 즉 해당 픽셀 주변의 정보들만을 이용한다는데에 그 한계가 있다. 그렇기 때문에 edge등이 소실되기 마련이고, 텍스쳐가 뭉개지는등의 단점이 있다. 이를 통해 노이즈가 제거되기도 하지만, 영상의 디테일 또한 소실되는 문제점이 있다.

 

원본 영상(좌), gaussian filter로 스무딩한 영상(우)

 

 

NL means 알고리즘은 다음과같은 아이디어에서 시작된다.

 

 

해당 픽셀의 gaussian kernel에 해당하는 local 영역의 픽셀값 대신,

영상내에서 해당 픽셀 주변 영역과 비슷한 패턴을 갖는 영역들을 찾아내

찾아낸 영역들에 대해서만 평균을 취하면 텍스쳐또한 살리면서 노이즈 제거가 가능하지 않을까?

 

 

연구는 항상 기발한 아이디어에서 시작되는것 같다. 아이디어만 읽고 생각해봐도 틀림없이 위에서 제안하는 결과가 나올것 같다.  

 

 

논문에서 NL means 알고리즘을 구현하는 방식은 다음과같이 제안한다. 우선, 몇가지 reference 코드들을 검색해봤는데, similarity window 영역에 대한 좌표 설정이 내가 이해하지 못한 다른 방법에 의해 결정되는것같아 제껴두고, 일정부분 참고하며 내가 이해한 방식으로만 코드를 작성해보았다.

NL means algorithm

 

1. denoising을 위한 픽셀에 대해 p(i)라고 하자.

2. NL means에 사용할 gaussian kernel을 준비한다. 이 kernel 사이즈는 Ai, Aj window 사이즈와 동일하게 설정한다.

3. 복원하려는 픽셀 p(i) 주변에 gaussian kernel과 동일한 크기의 주변 영역 Ai를 설정한다. 

4. 픽셀 p(i)주변에서 비슷한 패턴을 가진 영역을 탐색하기 위한 영역을 W라고 할때, W 내에서 영역 Aj를 sliding window 방식으로 탐색하며 다음의 식을 계산한다.

 

weight의 총합은 1

위의 식에서 exponential 항을 먼저 보면, Ai와 Aj의 euclidean distance 의 제곱, 즉 squared euclidean distance를 계산한 뒤, 표준편차가 a인 gaussian kernel을 적용시시키는것을 알 수 있다. h는 상수인 파라미터이다.

 

 

w(i, j)의 모든 j에 대한 합은 1이 되어아 햐므로, Z(i)는 exponential 항을 모두 합한 값으로 구해진다.

 

이렇게 구해진 window는 i 주변 window 내의 픽셀 j유사 확률정도를 나타내는 확률 값으로 볼 수 있다. window를 시각화해보면 다음과 같다.

 

빨간색 점에 해당하는 픽셀, 즉 복원하려는 픽셀 p(i)에 대한 similarity window는 우측과 같이 weight들이 계산된다. weight의 분포는 해당 픽셀주변 영역 Ai와 패턴이 비슷한 영역Aj에 대해 높은 가중치로 계산되어짐을 확인할 수 있다.

이렇게 구해진 weight를 noisy image에 해당하는 픽셀값들에 각각 곱한 뒤 합하여 복원되어진 p(i)의 픽셀값을 구하게 된다.

 

물론 이 sim window area를 넓힐수록 좋은성능을 보일것이 당연하나,, 계산복잡도가 무지막지하게 커진다.

 

 

 

Denoised image의 PSNR이 Noisy image에 비해 확연히 상승함을 확인할 수 있으며,

육안으로도 노이즈가 굉장히 많이 감소하였다.

다만, 과하게 노이즈가 제거된감이 있는데 이는 영상에따라 적절한 Similarity window, Area filter size, 그리고 상수 h값이 존재할텐데, 이를 실험적으로 구해야 한다는게 문제... 이것은 hand-crafted algorithm들이 모두 갖는 한계이기도 하니 어쩔수 없다... 무엇보다 가장 큰 단점은 속도인듯하다. 엄청 느림.

 

MATLAB 코드는 아래에 첨부한다.

function out_img = jy_NL_implementation(in_img,win_size,filt_size,h)
% Size of the image
[ori_r_size, ori_c_size]=size(in_img);

% Window
half_win_size = fix(win_size / 2);
half_filt_size = fix(filt_size / 2);

% Memory for the output
out_img=zeros(ori_r_size,ori_c_size);

% Replicate the boundaries of the input image
in_img_p = padarray(in_img,[half_filt_size half_filt_size],'symmetric');

% Used kernel
kernel = make_kernel(half_filt_size);
kernel = kernel / sum(kernel, 'all');

h=h*h;
for o_r_idx = 1:ori_r_size
    for o_c_idx = 1:ori_c_size
        p_r_idx = o_r_idx + half_filt_size;
        p_c_idx = o_c_idx + half_filt_size;
        
        p_r_idx_nei_st = max(half_filt_size+1, p_r_idx - half_win_size);
        p_r_idx_nei_ed = min(ori_r_size + half_filt_size, p_r_idx + half_win_size);
        p_c_idx_nei_st = max(half_filt_size+1, p_c_idx - half_win_size);
        p_c_idx_nei_ed = min(ori_c_size + half_filt_size, p_c_idx + half_win_size);
        
        N1 = in_img_p(p_r_idx-half_filt_size:p_r_idx+half_filt_size, p_c_idx-half_filt_size:p_c_idx+half_filt_size); 
        
        Z=0;
        NL=0;

        Win_img = in_img_p(p_r_idx_nei_st:p_r_idx_nei_ed,p_c_idx_nei_st:p_c_idx_nei_ed);
        W=zeros([p_r_idx_nei_ed - p_r_idx_nei_st, p_c_idx_nei_ed - p_c_idx_nei_st]);
        for p_r_idx_win = p_r_idx_nei_st:p_r_idx_nei_ed
            for p_c_idx_win = p_c_idx_nei_st:p_c_idx_nei_ed
                
                N2 = in_img_p(p_r_idx_win-half_filt_size:p_r_idx_win+half_filt_size, p_c_idx_win-half_filt_size:p_c_idx_win+half_filt_size);
                
                d = sum(kernel.*(N1-N2).*(N1-N2), 'all');
                w = exp(-d/h);
                W(p_r_idx_win - p_r_idx_nei_st + 1, p_c_idx_win - p_c_idx_nei_st + 1) = w;
                Z = Z + w;                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                
                NL = NL + (w*in_img_p(p_r_idx_win,p_c_idx_win));
            end
        end
%         if (o_r_idx==50 && o_c_idx == 50)
%             figure(1)
%             subplot(1,3,1)
%             imshow(uint8(in_img_p))
%             title('Original noisy image');
%             rectangle('Position',[p_r_idx_nei_st, p_c_idx_nei_st, p_r_idx_nei_ed - p_r_idx_nei_st, p_c_idx_nei_ed - p_c_idx_nei_st],'EdgeColor', 'r', 'LineWidth',3)
%             subplot(1,3,2)
%             imshow(uint8(Win_img))
%             title('Sim window area');
%             subplot(1,3,3)
%             imagesc(W)
%             colorbar
%             title('Sim window weight');
%             disp('a')
%         end
        out_img(o_r_idx, o_c_idx) = NL/Z;
    end
end

out_img = uint8(out_img);
end