Рассеяние вейвлетного изображения

Для проблем классификации часто полезно отобразить данные в некоторое альтернативное представление, которое отбрасывает неактуальную информацию, сохраняя при этом дискриминационные свойства каждого класса. Вейвлет-рассеяние изображений создает малодисперсионные представления изображений, которые нечувствительны к трансляциям и малым деформациям. Так как трансляции и небольшие деформации в изображении не влияют на принадлежность к классу, коэффициенты преобразования рассеяния обеспечивают функции, из которых можно построить надежные модели классификации.

Вейвлет-рассеяние работает путем каскадирования изображения посредством серии вейвлет-преобразований, нелинейностей и усреднения [1] [3] [5]. Результатом этого глубокого извлечения признаков является то, что изображения в одном классе перемещаются ближе друг к другу в представлении преобразования рассеяния, в то время как изображения, принадлежащие различным классам, перемещаются дальше друг от друга.

KTH-СОВЕТЫ

В этом примере используется общедоступная база данных текстур, база данных изображений KTH-TIPS (Textures in variable Illumination, Pose и Scale) [6]. В этом примере используется набор данных KTH-TIPS в градациях серого. Всего имеется 810 изображений с 10 текстурами и 81 изображение на текстуру. Большинство изображений размером 200 на 200. В этом примере предполагается, что был загружен набор данных в градациях серого KTH-TIPS, который был отменен так, что 10 классов текстур содержатся в отдельных подпапках общей папки. Каждая подпапка называется для класса текстур, который она содержит. Отмена ранжирования загруженного kth_tips_grey_200x200.tar файл достаточен для обеспечения KTH_TIPS папки верхнего уровня и требуемой структуры подпапок.

Используйте imageDatastore для считывания данных. Установите location имущества imageDatastore в папку, содержащую базу данных KTH-TIPS, к которой у вас есть доступ.

location = fullfile(tempdir,'kth_tips_grey_200x200','KTH_TIPS');
Imds = imageDatastore(location,'IncludeSubFolders',true,'FileExtensions','.png','LabelSource','foldernames');

Случайный выбор и визуализация 20 изображений из набора данных.

numImages = 810;
perm = randperm(numImages,20);
for np = 1:20
    subplot(4,5,np)
    im = imread(Imds.Files{perm(np)});
    imagesc(im);    
    colormap gray; axis off;
end

Классификация текстур

В этом примере используется MATLAB™'s параллельная обработка через tall интерфейс массива. Запустите параллельный пул, если он не запущен в данный момент.

if isempty(gcp)
    parpool;
end

Starting parallel pool (parpool) using the 'local' profile ...
Connected to the parallel pool (number of workers: 6).

Для воспроизводимости установите генератор случайных чисел. Перетасуйте файлы набора данных KTH-TIPS и разделите 810 изображений на два случайно выбранных набора, один для обучения и один удерживаемый набор для тестирования. Используйте приблизительно 80% изображений для построения прогностической модели из преобразования рассеяния и остаток для тестирования модели.

rng(100)
Imds = imageDatastore(location,'IncludeSubFolders',true,'FileExtensions','.png','LabelSource','foldernames');
Imds = shuffle(Imds);
[trainImds,testImds] = splitEachLabel(Imds,0.8);

Теперь у нас есть два набора данных. Обучающий набор состоит из 650 изображений, по 65 изображений на текстуру. Тестовый набор состоит из 160 изображений с 16 изображениями на текстуру. Чтобы проверить, подсчитайте метки в каждом наборе данных.

countEachLabel(trainImds)

ans=10×2 table
        Label         Count
    ______________    _____

    aluminium_foil     65  
    brown_bread        65  
    corduroy           65  
    cotton             65  
    cracker            65  
    linen              65  
    orange_peel        65  
    sandpaper          65  
    sponge             65  
    styrofoam          65  

countEachLabel(testImds)

ans=10×2 table
        Label         Count
    ______________    _____

    aluminium_foil     16  
    brown_bread        16  
    corduroy           16  
    cotton             16  
    cracker            16  
    linen              16  
    orange_peel        16  
    sandpaper          16  
    sponge             16  
    styrofoam          16  

Создать tall массивы для изображений с измененным размером.

Ttrain = tall(trainImds);
Ttest = tall(testImds);

Создание структуры рассеяния для входного изображения размером 200 на 200 с помощью InvarianceScale из 150. Гиперпараметр шкалы инвариантности является единственным, который мы задаем в этом примере. Для других гиперпараметров преобразования рассеяния используйте значения по умолчанию.

sn = waveletScattering2('ImageSize',[200 200],'InvarianceScale',150);

Для извлечения функций классификации для каждого обучающего и тестового наборов используйте helperScatImages_mean функция. Код для helperScatImages_mean находится в конце этого примера. helperScatImages_mean изменяет размер изображений до общего размера 200 на 200 и использует структуру рассеяния, sn, для получения матрицы признаков. В этом случае каждая матрица особенности составляет 391 на 7 на 7. Имеются 391 пути рассеяния и каждое изображение коэффициента рассеяния 7 на 7. Наконец, helperScatImages_mean получает среднее по 2-му и 3-му размерам для получения вектора признаков 391 элемента для каждого изображения. Это является значительным сокращением данных с 40 000 элементов до 391.

trainfeatures = cellfun(@(x)helperScatImages_mean(sn,x),Ttrain,'Uni',0);
testfeatures = cellfun(@(x)helperScatImages_mean(sn,x),Ttest,'Uni',0);

Используя tall gather соберите все векторы обучающих и тестовых функций и объедините их в матрицы.

Trainf = gather(trainfeatures);

Evaluating tall expression using the Parallel Pool 'local':
- Pass 1 of 1: Completed in 1 min 39 sec
Evaluation completed in 1 min 39 sec

trainfeatures = cat(2,Trainf{:});
Testf = gather(testfeatures);

Evaluating tall expression using the Parallel Pool 'local':
- Pass 1 of 1: Completed in 23 sec
Evaluation completed in 23 sec

testfeatures = cat(2,Testf{:});

Предыдущий код приводит к двум матрицам с размерами строки 391 и размером столбца, равными количеству изображений в обучающем и тестовом наборах соответственно. Таким образом, каждый столбец является вектором признаков.

Модель и прогноз PCA

В этом примере создается простой классификатор на основе основных компонентов векторов элементов рассеяния для каждого класса. Классификатор реализован в функциях helperPCAModel и helperPCAClassifier. Функция helperPCAModel определяет главные компоненты для каждого класса цифр на основе элементов рассеяния. Код для helperPCAModel находится в конце этого примера. Функция helperPCAClassifier классифицирует задержанные тестовые данные путем нахождения ближайшего совпадения (наилучшая проекция) между основными компонентами каждого вектора признаков теста с обучающим набором и соответствующего назначения класса. Код для helperPCAClassifier находится в конце этого примера.

model = helperPCAModel(trainfeatures,30,trainImds.Labels);
predlabels = helperPCAClassifier(testfeatures,model);

После построения модели и классификации тестового набора определите точность классификации тестового набора.

accuracy = sum(testImds.Labels == predlabels)./numel(testImds.Labels)*100

accuracy = 99.3750

Мы добились 99,375% правильной классификации или 0,625% частоты ошибок для 160 изображений в тестовом наборе. График матрицы путаницы показывает, что наша простая модель неправильно классифицировала одну текстуру.

figure
confusionchart(testImds.Labels,predlabels)

Резюме

В этом примере мы использовали вейвлет-рассеяние изображения, чтобы создать представления текстур с низкой дисперсией для классификации. Используя преобразование рассеяния и простой классификатор главных компонентов, мы достигли 99,375% корректной классификации в задержанном тестовом наборе. Этот результат сравним с современной производительностью базы данных KTH-TIPS [2].

Ссылки

[1] Бруна, J. и С. Маллэт. «Инвариантные сети свертки рассеяния». Транзакции IEEE по анализу шаблонов и машинному интеллекту. т. 35, № 8, 2013, с. 1872-1886.

[2] Хейман, Э., Б. Капуто, М. Фриц и Дж. О. Эклунд. «О значении реальных условий классификации материалов». В «Computer Vision - ECCV 2004» под редакцией Томаса Пайдлы и Йиржи Матаса, 3024: 253-66. Берлин, Гейдельберг: Springer Berlin Heidelberg, 2004. https://doi.org/10.1007/978-3-540-24673-2_21.

[3] Маллат, С. «Инвариантное рассеяние группы». Коммуникации в чистой и прикладной математике. Том 65, номер 10, 2012, стр. 1331-1398.

[4] Пужоль, О. и П. Радева. «Контролируемая классификация текстур для внутрисосудистой характеристики тканей». В «Handbook of Biomedical Image Analysis» под редакцией Джасджита С. Сури, Дэвида Л. Уилсона и Свами Лаксминараян, 57-109. Бостон, Массачусетс: Springer US, 2005. https://doi.org/10.1007/0-306-48606-7_2.

[5] Сифре, Л. и С. Маллат. Конференция IEEE 2013 по компьютерному зрению и распознаванию образов. 2013, стр. 1233-1240. 10.1109/CVPR.2013.163.

[6] Домашняя страница баз данных образов KTH-TIPS. https://www.csc.kth.se/cvap/databases/kth-tips/

Приложение - Вспомогательные функции

helperScatImages_mean

function features = helperScatImages_mean(sf,x)
x = imresize(x,[200 200]);
smat = featureMatrix(sf,x);
features = mean(mean(smat,2),3);
end

helperPCAModel

function model = helperPCAModel(features,M,Labels)
% This function is only to support wavelet image scattering examples in 
% Wavelet Toolbox. It may change or be removed in a future release.
% model = helperPCAModel(features,M,Labels)

% Copyright 2018 MathWorks

% Initialize structure array to hold the affine model
model = struct('Dim',[],'mu',[],'U',[],'Labels',categorical([]),'s',[]);
model.Dim = M;
% Obtain the number of classes
LabelCategories = categories(Labels);
Nclasses = numel(categories(Labels));
for kk = 1:Nclasses
    Class = LabelCategories{kk};
    % Find indices corresponding to each class
    idxClass = Labels == Class;
    % Extract feature vectors for each class
    tmpFeatures = features(:,idxClass);
    % Determine the mean for each class
    model.mu{kk} = mean(tmpFeatures,2);
    [model.U{kk},model.S{kk}] = scatPCA(tmpFeatures);
    if size(model.U{kk},2) > M
        model.U{kk} = model.U{kk}(:,1:M);
        model.S{kk} = model.S{kk}(1:M);
        
    end
    model.Labels(kk) = Class;
end

function [u,s,v] = scatPCA(x,M)
	% Calculate the principal components of x along the second dimension.

	if nargin > 1 && M > 0
		% If M is non-zero, calculate the first M principal components.
	    [u,s,v] = svds(x-sig_mean(x),M);
	    s = abs(diag(s)/sqrt(size(x,2)-1)).^2;
	else
		% Otherwise, calculate all the principal components.
        % Each row is an observation, i.e. the number of scattering paths
        % Each column is a class observation
		[u,d] = eig(cov(x'));
		[s,ind] = sort(diag(d),'descend');
		u = u(:,ind);
	end
end
end

helperPCAClassifier

function labels = helperPCAClassifier(features,model)
% This function is only to support wavelet image scattering examples in 
% Wavelet Toolbox. It may change or be removed in a future release.
% model is a structure array with fields, M, mu, v, and Labels
% features is the matrix of test data which is Ns-by-L, Ns is the number of
% scattering paths and L is the number of test examples. Each column of
% features is a test example.

% Copyright 2018 MathWorks

labelIdx = determineClass(features,model); 
labels = model.Labels(labelIdx); 
% Returns as column vector to agree with imageDatastore Labels
labels = labels(:);


%--------------------------------------------------------------------------
function labelIdx = determineClass(features,model)
% Determine number of classes
Nclasses = numel(model.Labels);
% Initialize error matrix
errMatrix = Inf(Nclasses,size(features,2));
for nc = 1:Nclasses
    % class centroid
    mu = model.mu{nc};
    u = model.U{nc};
    % 1-by-L
    errMatrix(nc,:) = projectionError(features,mu,u);
end
% Determine minimum along class dimension
[~,labelIdx] = min(errMatrix,[],1);   


%--------------------------------------------------------------------------
function totalerr = projectionError(features,mu,u)
    %
    Npc = size(u,2);
    L = size(features,2);
    % Subtract class mean: Ns-by-L minus Ns-by-1
    s = features-mu;
    % 1-by-L
    normSqX = sum(abs(s).^2,1)';
    err = Inf(Npc+1,L);
	err(1,:) = normSqX;
    err(2:end,:) = -abs(u'*s).^2;
    % 1-by-L
    totalerr = sqrt(sum(err,1));
end
end
end