Total variation 최소화로 이미지 노이즈 제거하기

이미지 노이즈 제거(denoising)는 컴퓨터 비전의 주요 분야 중 하나로, 노이즈가 섞인 관측으로부터 깨끗한 이미지를 복원하는 것을 목표로 한다. Rudin, Osher, Fatemi(ROF)가 소개한 total variation 최소화(total variation minimization)는 이미지 그래디언트와 데이터 충실도(data fidelity)에 페널티를 부여함으로써 이 문제를 해결한다. 이 글은 ROF 모델의 수학적 유도를 탐구하고 Chambolle-Pock Primal-Dual 알고리즘으로 이를 구현한다.

위 인터랙티브 데모는 total variation 최소화로 노이즈를 제거하는 과정을 보여준다. 보기 모드에서 원시(primal) 변수나 쌍대(dual) 변수 중 하나를 표시하도록 선택할 수 있다. 원시 변수는 노이즈가 제거된 이미지 자체를 나타낸다. 반면 쌍대 변수는 빨강과 초록 채널에 인코딩된 2D 벡터를 시각화하는데, 각각 x 방향과 y 방향에 대응하며 두 성분 모두 [-1, 1]에서 [0, 1]로 선형 스케일링되고 파랑 채널은 0.5로 고정된다. 파라미터 값을 조정하여 최적화 과정이 어떻게 전개되는지 관찰할 수도 있다.

이미지 노이즈 제거

이미지 노이즈 제거 시나리오에서, 관측된 이미지 \(f\)가 평균이 0인 가산 가우시안 노이즈 \(\eta\)로 손상되었다고 가정한다.

\[f = u + \eta\]

\(u\)를 복원하기 위해, 서로 상충하는 두 목표의 균형을 맞추는 최적화 문제를 푼다.

• 충실도(Fidelity): 복원된 이미지 \(u\)는 관측 \(f\)에 가까워야 한다.

• 정규화(Regularization): 복원된 이미지 \(u\)는 매끄러워야 한다.

위 두 항을 더하여, 최소화할 일반적인 목적 함수는

\[E(u) = \mathcal{R}(u) + \frac{1}{2\lambda} \|u - f\|_2^2\]

이며, 여기서 \(\mathcal{R}(u)\)는 정규화 항이고 \(\lambda\)는 두 항 사이의 트레이드오프를 제어한다. 두 항을 동시에 고려함으로써, 복원된 이미지는 관측값을 어느 정도 따르면서, 매끄러운 표면을 가지도록 강제된다.

Total variation 최소화

\(\mathcal{R}(u)\)는 출력의 특성을 결정한다. 특히 total variation(TV) 정규화 항은 \(L_1\) 노름으로 정의된다. \(L_1\) 노름은 \(L_2\) 노름보다 큰 이상치(이미지 경계선)에 덜 민감하므로, 평탄한 영역의 노이즈 변동을 0으로 억제하면서 날카로운 값 변화를 상대적으로 허용한다.

따라서 TV 노이즈 제거 모델은 다음과 같이 정의된다.

\[\min_{u} \int_{\Omega} \|\nabla u\| dx + \frac{1}{2\lambda} \|u - f\|_2^2\]

여기서 \(\Omega\)는 이미지 영역이다.

Primal-Dual 방법

그런데 \(L_1\) 노름 \(\|x\|\)은 원점 \(x=0\)에서 미분 불가능하다. 이를 효율적으로 풀기 위해, 쌍대 노름(dual norm)을 사용하여 최소화 문제를 안장점(saddle-point) 문제로 변환한다.

쌍대 공식화

\(L_1\) 노름의 성질을 이용해 TV 항을 다음과 같이 다시 쓸 수 있다.

\[\|\nabla u\|_1 = \max_{\|p\|_\infty \le 1} \langle p, \nabla u \rangle\]

여기서 \(p\)는 쌍대 변수다. 그러면 최적화 문제는 다음과 같이 된다.

\[\min_{u} \max_{\|p\|_\infty \le 1} {E(u, p)} = \min_{u} \max_{\|p\|_\infty \le 1} \left( \langle p, \nabla u \rangle + \frac{1}{2\lambda} \|u - f\|_2^2 \right)\]

Chambolle-Pock 알고리즘

Chambolle-Pock 알고리즘은 변수 \(p\)와 \(u\)를 반복적으로 갱신한다.

• dual variable update (\(p\)에 대한 경사 상승)

먼저 목적 함수를 \(p\)에 대해 최소화한다.

\[p^{n+1}=\arg\max_{p \in P}​ \left( {\langle p,\nabla \bar{u}^n \rangle + \frac{1}{2\lambda} \|u - f\|_2^2 − \frac{1}{2\sigma}\|p−p^n\|_2^2} \right)\]

\(-\frac{1}{2\sigma} \|p - p^n\|_2^2\) 항은 안정성을 보장하는 근접(proximal) 페널티다. 근접 페널티의 마이너스 부호는 최적화 문제의 오목성(concavity)을 보존한다. 도함수를 0으로 놓으면 다음을 얻는다.

\[\nabla \bar{u}^n - (p^{n+1}−p^n) = 0\]

그리고 \(p^{n+1}\)을 \(L_\infty\) 공의 표면에 사영한다.

\[p^{n+1} = \text{proj}_P (p^n + \sigma \nabla \bar{u}^n)\]

이 사영은 \(L_1\) 노름으로 정규화하여 쌍대 제약 \(P = \{p: \|p\|_\infty \le 1\}\)을 보장한다.

• primal variable update (\(u\)에 대한 경사 하강)

수반 연산자(adjoint operator) 성질 \(\langle p, \nabla u \rangle = \langle \nabla^* p, u \rangle = \langle -\text{div } p, u \rangle\)를 이용해, 목적 함수를 \(u\)에 대해 최소화한다.

\[u^{n+1} = \arg\min_{u} \left( \langle -\text{div } p^{n+1}, u \rangle + \frac{1}{2\lambda} \|u - f\|_2^2 + \frac{1}{2\tau} \|u - u^n\|_2^2 \right)\]

\(\frac{1}{2\tau} \|u - u^n\|_2^2\) 항은 안정성을 보장하는 근접 페널티다. 도함수를 0으로 놓으면 다음을 얻는다.

\[-\text{div } p^{n+1} + \frac{1}{\lambda}(u^{n+1} - f) + \frac{1}{\tau} (u^{n+1} - u^n) = 0\]

그리고

\[u^{n+1} = \frac{u^n + \tau \text{div } p^{n+1} + \frac{\tau}{\lambda} f}{1 + \frac{\tau}{\lambda}}.\]

이다.

• 과이완(Over-relaxation)

추가로 수렴 속도를 높이기 위해, 아래와 같이 보조 변수를 두어 primal 변수를 업데이트할 수 있다.

\[\bar{u}^{n+1} = u^{n+1} + \theta(u^{n+1} - u^n)\]

보통 \(\theta=1\)을 사용한다.

예제 코드

아래는 위에서 유도한 공식들을 이용해 이미지의 노이즈를 제거하는 코드이다.

function u = tv_denoise(f, lambda, iter)
    % f: Noisy input, lambda: Regularization weight, iter: Iterations
    [rows, cols] = size(f);
    f = double(f);
    
    % Parameters for convergence (L^2 <= 8)
    L = sqrt(8);
    tau = 0.01;
    sigma = 1/(L^2 * tau);
    theta = 1;
    
    % Initialize variables
    u = f;
    u_bar = f;
    p = zeros(rows, cols, 2);
    
    for i = 1:iter
        % Dual Update
        ux = [diff(u_bar, 1, 2), zeros(rows, 1)];
        uy = [diff(u_bar, 1, 1); zeros(1, cols)];
        p(:,:,1) = p(:,:,1) + sigma * ux;
        p(:,:,2) = p(:,:,2) + sigma * uy;
        
        % Magnitude for projection
        L_infty = max(abs(p(:,:,1)), abs(p(:,:,2)));
        p(:,:,1) = p(:,:,1) ./ max(1, L_infty);
        p(:,:,2) = p(:,:,2) ./ max(1, L_infty);
        
        % Primal Update
        u_old = u;
        div_p = [p(:,1,1), diff(p(:,:,1), 1, 2)] + ...
                [p(1,:,2); diff(p(:,:,2), 1, 1)];
        
        u = (u + tau * div_p + (tau/lambda) * f) / (1 + tau/lambda);
        
        % Relaxation
        u_bar = u + theta * (u - u_old);
    end
end

다른 정규화 항

• \(L_2\) 노름 \(\|\nabla u\|_2^2\): 고주파수를 전역적으로 페널티하여 등방성 확산(isotropic diffusion)을 일으킨다. 제곱 노름이 작은 그래디언트보다 큰 그래디언트를 더 페널티하므로 에지가 흐려진다. \(L_2\) 노름은 모든 값에서 미분 가능하므로, 최적화 문제를 푸는 데 쌍대 변수를 사용할 필요가 없다.

• 후버(Huber) 노름 \(h_{\epsilon} (\nabla u)\): \(L_1\)과 \(L_2\) 노름의 혼합으로, 작은 값에는 \(L_2\) 노름을, 큰 값에는 \(L_1\) 노름을 사용한다. \(L_1\) 정규화 항의 계단 효과(staircase effect)를 막으면서도 \(L_2\) 정규화 항보다 큰 그래디언트에 덜 민감하다.

정리하며

Primal-Dual 방법을 통한 TV 최소화는 이미지 노이즈 제거에 잘 활용되는 방법이다. 단순한 필터와 달리, 그래디언트의 \(L_1\) 노름을 이용해 에지의 존재를 명시적으로 모델링한다. 계단 효과를 일으킬 수 있지만, 노이즈가 섞인 이미지의 구조적 무결성을 잘 보존할 수 있다.