디지털 영상 처리 - 히스토그램 처리 함수 구현

 안녕하세요. 지난 포스팅의 디지털 영상 처리 - 지역적 히스토그램 처리를 마지막으로 히스토그램과 관련된 이론적인 내용을 모두 끝냈습니다. 오늘은 그 동안 알아본 히스토그램(Histogram), 히스토그램 평활화(Histogram Equalization), 히스토그램 지정(Histogram Matching), 지역 히스토그램 평활화(Local Histogram Equalization)을 구현해보겠습니다. 이미 지난 포스팅에서 자세한 내용은 알아보았기 때문에 간단하게 결과 그림과 코드만 제공해드리도록 하겠습니다. 만약 코드를 사용하고 싶으시면 해당 코드들은 제 아래의 깃허브 링크에 정리해놓았으니 확인해보시길 바랍니다. 

github.com/skawngus1111/DIP

skawngus1111/DIP

 

1. 히스토그램(Histogram)

function [hist] = im2hist(image, minValue, maxValue)
% im2his : 입력 영상을 히스토그램으로 변환하는 함수
% input
%   image : 입력 영상
%   minValue : 히스토그램의 최솟값
%   maxValue : 히스토그램의 최댓값
% output 
%   histogram : 입력 영상에 대한 히스토그램

[height, width] = size(image);
hist = zeros(1, maxValue);

for i=1:height
    for j=1:width
        value = image(i, j);
        if (value+1 >= minValue) && (value+1 <= maxValue)
            hist(value+1) = hist(value+1) + 1;
        end
    end
end
end

2. 히스토그램 평활화(Histogram Equalization)

 다음으로 히스토그램 평활화입니다. 대비를 향상시키기 위해서 영상의 히스토그램을 강제로 균등확률분포를 가지도록 만들어주는 변환이라고 언급하였습니다. 

function [new, T] = histequal(image, minValue, maxValue)
% histequal : 입력 영상에 히스토그램 평활화(Histogram Equalization)을 적용하는 함수
%   input
%       1. image : 입력 영상
%       2. minValue, maxValue : 입력 영상 밝기값의 최솟값과 최댓값
%   output 
%       1. new : 평활화된 영상
%       2. cdf : 입력 영상의 누적확률분포(Cumulative Probability Distribution)

% Basic Algorithm
%     1. 입력 영상의 히스토그램을 얻는다.
%     2. 히스토그램을 정규화한다.(확률화)
%     3. 변환 테이블을 만들기 위해 정규화된 히스토그램으로부터 각 밝기별로 CDF를 계산한다.
%     4. STEP 3을 통해 얻은 변환 테이블을 기반으로 r -> s로의 갱신을 수행한다.

% Toy Example

[height, width] = size(image); % 입력 영상의 크기
totalPixel = width * height; % 전체 픽셀의 개수
hist = im2hist(image, minValue, maxValue) / totalPixel; % 히스토그램 정규화(0 ~ 1)

T = zeros(1, maxValue); % 입력 영상 히스토그램의 CDF를 저장할 변수

T(1) = hist(1);
for i=2:maxValue % 최솟값 ~ 최댓값까지 각 값에 대한 CDF 계산
    T(i) = T(i-1) + hist(i);
end

T = roundn((maxValue-1) * T, 0);

new = zeros([height, width], 'uint8');

for i=1:height
    for j=1:width
        new(i, j) = T(image(i, j) + 1);
    end
end

end

3. 히스토그램 지정(Histogram Matching)

function [MatchedImage, GInverse, G] = histMatching(SourceImage, TargetImage, minValue, maxValue)
%histMatching : sourceImage를 targetImage의 히스토그램으로 매칭시키는 함수
%   input 
%       1. sourceImage : 매칭시킬 source 영상
%       2. targetImage : 매칭시킬 target 히스토그램을 가진 영상
%       3. minValue, maxValue : 입력 영상 밝기값의 최솟값과 최댓값
%   output
%       1. MatchedImage : target 영상을 이용해서 source 영상을 매칭시킨 출력 영상
%       2. GInverse : s -> z로 변환하는 G의 역변환
%       3. G : z -> s로의 변환

% Basic Algorithm
%   1. source 영상을 히스토그램 평활화적용
%   2. target 영상에서 평활화된 source 영상의 히스토그램을 갖기 위한 함수 G 계산
%       G(z) = \sum_{j=1}^{k}p_{z}(z_{j})
%   3. 모호성 문제(ambiguous problem)를 방지하기 위해서 s_{k}의 값 중 가장 가까운 G(z_{q}) 찾기
%   4. STEP 3에서 찾은 s_{k}와 G(z_{k})의 쌍을 이용하여 s_{k} -> z_{q}로 매핑

% STEP 1. source 영상(r)을 히스토그램 평활화 적용(s) : r -> s
[~, T] = histequal(SourceImage, minValue, maxValue);

% STEP 2. target 영상에서 평활화된 source 영상의 히스토그램을 갖기 위한 함수 G 계산
[~, G] = histequal(TargetImage, minValue, maxValue);

% STEP 3. 모호성 문제를 방지하기 위해서 s_{k}의 값 중 가장 가까운 G(z_{q}) 찾기
GInverse = zeros(1, maxValue);

for i=1:maxValue
    diff = abs(T(i) - G);
    [~, idx] = min(diff);
    GInverse(i) = idx - 1;
end


% STEP 4. STEP 3에서 찾은 s_{k}와 G(z_{k})의 쌍을 이용하여 s_{k} -> z_{q}로 매핑
MatchedImage = zeros([size(SourceImage, 1), size(SourceImage, 2)], 'uint8');

for i=1:size(SourceImage, 1)
    for j=1:size(SourceImage, 2)
        MatchedImage(i, j) = GInverse(SourceImage(i, j) + 1);
    end
end

end

4. 지역 히스토그램 평활화(Local Histogram Equalization)

지역 히스토그램 평활화는 윈도우의 크기와 이동 픽셀크기를 얼마나 잡느냐에 따라 그 결과가 달라집니다.

function [EqualizedImage] = LocalHistogramEqualization(image, minValue, maxValue, kernelSize, Stride)
%LocalHistogramEqualization : 지역 히스토그램 평활화를 수행하는 함수
%   input 
%       1. image : 입력 영상
%       2. minValue, maxValue : 입력 영상 밝기값의 최솟값과 최댓값
%       3. kernelSize : 지역 히스토그램을 적용할 윈도우의 크기
%       4. Stride : 픽셀을 몇 칸 건너서 이동할 건 지?
%   output
%       1. MatchedImage : target 영상을 이용해서 source 영상을 매칭시킨 출력 영상
%       2. GInverse : s -> z로 변환하는 G의 역변환
%       3. G : z -> s로의 변환

% Basic Algorithm
%   1. kernel size만큼 윈도우 결정
%   2. 윈도우 내부에서 히스토그램 평활화 적용
%   3. Stride만큼 중심 픽셀을 이동하여 STEP 1 ~ 2 반복

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

[height, width] = size(image);
n = floor(kernelSize/2);
EqualizedImage = zeros([height, width], 'uint8');

for i=1:Stride:height
    if (i <= n) || (i >= (height - n))
        continue
    end
    for j=1:Stride:width
        if (j <= n) || (j >= (width - n))
            continue
        end
        [temp, ~] = histequal(image(i-n:i+n, j-n:j+n), minValue, maxValue);
        EqualizedImage(i-n:i+n, j-n:j+n) = temp; 
    end
end

end