Applying vector median filter on RGB image based on matlab

前言：

最近想看看矢量中值濾波（Vector median filter， VMF）在GRB圖像上的濾波效果，意外的是找了一大圈卻發現網上沒有現成的code，因此經過matab親自實現了一個，須要學習的朋友能夠拿過去用。本文的核心是VMF的matlab實現，最後經過在RGB圖像上應用舉例說明。

 

VMF的數學表達：

 含有N個矢量的集合{C₁，C₂，...C_N}，它的VMF結果以下所示：

其中，C_VM1表示距離全部其餘向量的距離和最小的那個向量。而距離能夠本身定義，經常使用的歐氏距離，曼哈頓距離等等。

matlab code：

%   向量中值濾波,根據歐式距離進行濾波，將RGB圖像中每一個像素位置的三個顏色做爲一個向量總體處理。

% pdist函數解釋
% Pairwise distance between pairs of objects
% Syntax
% D = pdist(X)
% D = pdist(X,distance)
% Description
%
% D = pdist(X)
% 計算 X 中各對行向量的相互距離(X是一個m-by-n的矩陣). 這裏 D 要特別注意，
% D 是一個長爲m(m–1)/2的行向量.能夠這樣理解 D 的生成：首先生成一個 X 的距離方陣，
% 因爲該方陣是對稱的，且對角線上的元素爲0，因此取此方陣的下三角元素，按照Matlab中矩陣的按列存儲原則，
% 此下三角各元素的索引排列即爲(2,1), (3,1), ..., (m,1), (3,2), ..., (m,2), ..., (m,m–1).
% 能夠用命令 squareform(D) 將此行向量轉換爲原距離方陣.(squareform函數是專門幹這事的，其逆變換是也是squareform。)
%
% D = pdist(X,distance) 使用指定的距離.

vectors_set = rand(5,3);%向量集合中的每一行是一個向量，一共含有5個向量

dist = pdist(vectors_set,'euclidean');%使用歐式距離計算向量之間的距離
dist = squareform(dist);%還原回矩陣的形式
dist = sum(dist,2);%求出每一個向量到其餘全部向量的距離和
indx = find(dist==min(dist));%找到距離和最小的向量的index
median_vec = vectors_set(indx,:) %根據index找出中值向量

 

RGB圖像上的應用：

%   向量中值濾波,根據歐式距離進行濾波，能夠用於RGB圖像的去噪。

clc;    % Clear command window.
clear;    % Delete all variables.
close all;    % Close all figure windows except those created by imtool.
imtool close all;    % Close all figure windows created by imtool.
workspace;    % Make sure the workspace panel is showing.
fontSize = 15;
% Read in a standard MATLAB color demo image.
folder = fullfile(matlabroot, '\toolbox\images\imdemos');
baseFileName = 'peppers.png';
% Get the full filename, with path prepended.
fullFileName = fullfile(folder, baseFileName);
if ~exist(fullFileName, 'file')
    % Didn't find it there.  Check the search path for it.
    fullFileName = baseFileName; % No path this time.
    if ~exist(fullFileName, 'file')
        % Still didn't find it.  Alert user.
        errorMessage = sprintf('Error: %s does not exist.', fullFileName);
        uiwait(warndlg(errorMessage));
        return;
    end
end
rgbImage = imread(fullFileName);
subplot(1, 3, 1);
imshow(rgbImage);
title('Original color Image', 'FontSize', fontSize);
subplot(1, 3, 2);
noisyRGB = imnoise(rgbImage,'salt & pepper', 0.05);
imshow(noisyRGB);
title('Image with Salt and Pepper Noise', 'FontSize', fontSize);
subplot(1, 3, 3);
resortedImage = noisyRGB;

[h,w,bands] = size(rgbImage);%彩色圖像的大小
win_radius = 1;%濾波窗口半徑大小
win_size = (2*win_radius+1).^2;%濾波窗口的大小
for i = 1 + win_radius : h - 2
    for j = 1 + win_radius : w-2
        vectors_set = reshape(rgbImage(i-win_radius:i+win_radius,j-win_radius:j+win_radius,:),win_size,1,bands);%將窗口內的向量進行變形
        vectors_set = reshape(permute(vectors_set,[2,3,1]),bands,[])';%將向量矩陣變爲每一行爲一個向量的模式
        [n,d] = size(vectors_set);%n是向量個數，d是向量維度
        % pdist
        % Pairwise distance between pairs of objects
        % Syntax
        % D = pdist(X)
        % D = pdist(X,distance)
        % Description
        %
        % D = pdist(X)
        % 計算 X 中各對行向量的相互距離(X是一個m-by-n的矩陣). 這裏 D 要特別注意，
        % D 是一個長爲m(m–1)/2的行向量.能夠這樣理解 D 的生成：首先生成一個 X 的距離方陣，
        % 因爲該方陣是對稱的，且對角線上的元素爲0，因此取此方陣的下三角元素，按照Matlab中矩陣的按列存儲原則，
        % 此下三角各元素的索引排列即爲(2,1), (3,1), ..., (m,1), (3,2), ..., (m,2), ..., (m,m–1).
        % 能夠用命令 squareform(D) 將此行向量轉換爲原距離方陣.(squareform函數是專門幹這事的，其逆變換是也是squareform。)
        %
        % D = pdist(X,distance) 使用指定的距離.distance能夠取下面圓括號中的值，用紅色標出！
        dist = pdist(vectors_set,'euclidean');
        dist = squareform(dist);
        dist = sum(dist,2);
        indx = find(dist==min(dist));
        median_vec = vectors_set(indx,:);        
        resortedImage(i,j,:)= median_vec(1,:);       
    end
end

imshow(resortedImage);
title('Restored Image', 'FontSize', fontSize);

 

實驗結果：

 

總結：算法的實現並不複雜，可是運算時間較長，能夠進一步優化。

 

另外：在針對RGB圖像的矢量濾波上面，有一種作法是color plane by color plane的作法，雖然已經不算是vector median filter-based的方法了，可是效果還好，能夠借鑑

 

matlab code:

%
% 一個plane一個plane的單獨處理

clc;    % Clear command window.
clear;    % Delete all variables.
close all;    % Close all figure windows except those created by imtool.
imtool close all;    % Close all figure windows created by imtool.
workspace;    % Make sure the workspace panel is showing.
fontSize = 15;
% Read in a standard MATLAB color demo image.
folder = fullfile(matlabroot, '\toolbox\images\imdemos');
baseFileName = 'peppers.png';
% Get the full filename, with path prepended.
fullFileName = fullfile(folder, baseFileName);
if ~exist(fullFileName, 'file')
    % Didn't find it there.  Check the search path for it.
    fullFileName = baseFileName; % No path this time.
    if ~exist(fullFileName, 'file')
        % Still didn't find it.  Alert user.
        errorMessage = sprintf('Error: %s does not exist.', fullFileName);
        uiwait(warndlg(errorMessage));
        return;
    end
end
rgbImage = imread(fullFileName);
% Get the dimensions of the image.  numberOfColorBands should be = 3.
[rows columns numberOfColorBands] = size(rgbImage);
% Display the original color image.
subplot(3, 4, 1);
imshow(rgbImage);
title('Original color Image', 'FontSize', fontSize);
% Enlarge figure to full screen.
set(gcf, 'Position', get(0,'Screensize')); 
% Extract the individual red, green, and blue color channels.
redChannel = rgbImage(:, :, 1);
greenChannel = rgbImage(:, :, 2);
blueChannel = rgbImage(:, :, 3);
% Display the individual red, green, and blue color channels.
subplot(3, 4, 2);
imshow(redChannel);
title('Red Channel', 'FontSize', fontSize);
subplot(3, 4, 3);
imshow(greenChannel);
title('Green Channel', 'FontSize', fontSize);
subplot(3, 4, 4);
imshow(blueChannel);
title('Blue Channel', 'FontSize', fontSize);
% Generate a noisy image.  This has salt and pepper noise independently on
% each color channel so the noise may be colored.
noisyRGB = imnoise(rgbImage,'salt & pepper', 0.05);
subplot(3, 4, 5);
imshow(noisyRGB);
title('Image with Salt and Pepper Noise', 'FontSize', fontSize);
% Extract the individual red, green, and blue color channels.
redChannel = noisyRGB(:, :, 1);
greenChannel = noisyRGB(:, :, 2);
blueChannel = noisyRGB(:, :, 3);
% Display the noisy channel images.
subplot(3, 4, 6);
imshow(redChannel);
title('Noisy Red Channel', 'FontSize', fontSize);
subplot(3, 4, 7);
imshow(greenChannel);
title('Noisy Green Channel', 'FontSize', fontSize);
subplot(3, 4, 8);
imshow(blueChannel);
title('Noisy Blue Channel', 'FontSize', fontSize);
% Median Filter the channels:
redMF = medfilt2(redChannel, [3 3]);
greenMF = medfilt2(greenChannel, [3 3]);
blueMF = medfilt2(blueChannel, [3 3]);
% Find the noise in the red.
noiseImage = (redChannel == 0 | redChannel == 255);
% Get rid of the noise in the red by replacing with median.
noiseFreeRed = redChannel;
noiseFreeRed(noiseImage) = redMF(noiseImage);
% Find the noise in the green.
noiseImage = (greenChannel == 0 | greenChannel == 255);
% Get rid of the noise in the green by replacing with median.
noiseFreeGreen = greenChannel;
noiseFreeGreen(noiseImage) = greenMF(noiseImage);
% Find the noise in the blue.
noiseImage = (blueChannel == 0 | blueChannel == 255);
% Get rid of the noise in the blue by replacing with median.
noiseFreeBlue = blueChannel;
noiseFreeBlue(noiseImage) = blueMF(noiseImage);
% Reconstruct the noise free RGB image
rgbFixed = cat(3, noiseFreeRed, noiseFreeGreen, noiseFreeBlue);
subplot(3, 4, 9);
imshow(rgbFixed);
title('Restored Image', 'FontSize', fontSize);

實驗結果：