Сопоставление внешней формы объекта с использованием нормализованной взаимной корреляции

Хорошо, давайте попробуем... Это решение пытается использовать существующую реализацию normxcorr2 и изменить ее, чтобы решить вашу проблему.

Формула нормализованной взаимной корреляции:

В этом случае вы хотите изменить границы интеграции для каждого окна. Это, в свою очередь, влияет как на стандартные отклонения, так и на саму корреляцию. Разберемся с этим в несколько шагов:

Шаг 1. Получите правильную корреляцию

Мы можем сделать это, изменив изображение шаблона:

template_fix = template;
mean_template_mask = mean(template(mask == 1));
template_fix(mask == 0) = mean_template_mask;
result = normxcorr2(template_fix, query)

Обратите внимание, что, внося это изменение, мы делаем среднее значение шаблона равным среднему значению шаблона внутри маски. таким образом, все пиксели шаблона за пределами маски не участвуют в интегрировании, поскольку они равны среднему значению.

Шаг 2. Исправьте стандартный шаблон

size_mask = sum(mask(:));
size_template = numel(template);
std_template = std2(template);
std_template_masked = sqrt(sum((template_fix(:) - mean_template_mask).^2)/size_mask);
result = result * (std_template/std_template_masked);

Шаг 3. Исправьте стандартный запрос

sum_filt = ones(size(template));
std_query = filter2(query.^2, sum_filt) - filter2(query, sum_filt).^2/size_template;
std_query = sqrt(std_query/size_template);

std_query_mask = filter2(query.^2, mask) - filter2(query, mask).^2/size_mask;    
std_query_mask = sqrt(std_query_mask/size_mask);

result = result .* std_query ./ std_query_mask;

Мой Matlab не отвечает, поэтому у меня не было возможности проверить его на практике. Если я не пропустил некоторые ошибки, он должен быть математически эквивалентен.

Это решение выполняет некоторые дополнительные свертки, но не обрабатывает перекрывающиеся пиксели более одного раза.

Если вы используете одно и то же изображение шаблона несколько раз, вы можете реорганизовать шаги 1 и 2, чтобы выполнить их только один раз для предварительной обработки. Хотя оба не должны быть вычислительно дорогими.

Другой подход: прямой

Вот другой, простой подход, в котором не используется исходная функция normxcorr2. Этот код можно легко оптимизировать для использования памяти в ущерб удобочитаемости.

% q for query, t for template and mask for mask
% shape = 'full' or 'same' or 'valid'

t_mask = t .* mask;
n      = numel(mask);
tq_sum = filter2(t_mask,q, shape);

q_sum  = filter2(mask, q, shape);    
q_mean = q_sum/n;
t_sum  = sum(t_mask(:));
t_mean = t_sum/n;

res1 = tq_sum - t_mean*q_sum - t_sum*q_mean + t_mean*q_mean*n;

t_std = sqrt((sum(t_mask(:).^2) - sum(t_mask(:)).^2/n)/n);
q_std = sqrt((filter2(mask, q.^2, shape) - q_sum.^2/n)/n);

res = res1 ./ (n * t_std * q_std)

Это производная от предыдущего сообщения, где я дал ответ здесь: шаблон Matlab, соответствующий только 0 (или 1) в матрице

Однако в этом решении использовались циклы for, что весьма неэффективно. Таким образом, мы можем использовать im2col, bsxfun и col2im, чтобы сделать это быстрее. im2col берет перекрывающиеся области изображения и помещает их в отдельные столбцы. По сути, это берет скользящие окна вашего изображения, как и при любой пространственной фильтрации изображения, собирает все пиксели в скользящем окне и помещает каждое окно в виде отдельных столбцов.

Предположим, что размер вашего шаблона равен M x N, а размер вашего изображения, в котором вы хотите выполнить поиск, равен R x C, и предположим, что ваш шаблон изображения называется imTemplate, а изображение, которое вы хотите найти, имеет значение imSearch, мы можем выполнить следующую настройку. Предположим также, что оба изображения являются бинарными.

[M, N] = size(imTemplate);
[R, C] = size(imSearch);

%// Cast to double for precision
imTemplate = im2double(imTemplate);
imSearch = im2double(imSearch);

neigh = im2col(imSearch, [M, N]);
templateCol = imTemplate(:); %// Ensures we place template into single column

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

imInvert = ~imTemplate;
imInvertNoBorder = imclearborder(imInvert, 8); %// Search 8-pixel neighbourhood

Мы будем использовать это, чтобы выяснить, какие пиксели мы собираемся удалить из поиска. Это можно сделать:

rowsToRemove = imInvertNoBorder(:) == 1;

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

neigh(rowsToRemove,:) = [];

Теперь мы можем вычислить NCC по всем этим столбцам. Если вы помните, NCC между двумя сигналами выглядит следующим образом:

_{(источник: www.jot.fm)}

Таким образом, нам нужно вычесть среднее значение из каждой окрестности, а также нам нужно вычесть среднее значение из каждого из столбцов. Затем мы вычисляем формулу, как показано выше. Мы можем легко сделать это векторизованным в MATLAB следующим образом:

neighMeanSubtract = bsxfun(@minus, neigh, mean(neigh));
templateMeanSubtract = templateCol - mean(templateCol);

Мы можем вычислить числитель NCC для каждой окрестности (до суммирования) следующим образом:

numerator = bsxfun(@times, neighMeanSubtract, templateMeanSubtract);

Теперь все, что нам нужно сделать, это суммировать все столбцы, и это даст нам наш окончательный числитель:

sumNumerator = sum(numerator);

Знаменатель можно вычислить так:

denominator1 = sqrt(sum(neighMeanSubtract.*neighMeanSubtract));
denominator2 = sqrt(sum(templateMeanSubtract.*templateMeanSubtract));
sumDenominator = denominator1 .* denominator2;

Наконец, наш NCC можно вычислить следующим образом:

NCC = sumNumerator ./ sumDenominator;

Вы заметите, что это одна строка значений. Каждая строка соответствует выходу, определенному в окрестности. Таким образом, нам также нужно преобразовать это обратно в матрицу, поэтому вы можете использовать col2im:

finalOutput = col2im(NCC, [M, N], [R, C]);

Приведенный выше оператор возьмет перекрывающиеся окрестности M x N, определенные в NCC, и изменит их форму, чтобы она стала матрицей R x C. Иногда вы получаете сообщение об ошибке делите на ноль, особенно если окно поиска соседей однородно. Таким образом, вы получите NaN чисел. Предполагается, что области без изменений не имеют корреляции при обработке изображений, поэтому давайте обнулим эти местоположения:

finalOutput(isnan(finalOutput)) = 0;

Если вы хотите найти место с наибольшей корреляцией, просто выполните:

[rowNCC, colNCC] = find(finalOutput == max(finalOutput(:)));

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

maxCoeff = max(abs(finalOutput(:)));
[rowNCC, colNCC] = find(abs(finalOutput) == maxCoeff);

Для вашего удовольствия от копирования и вставки вот весь код:

function [rowNCC, colNCC] = testCorr(imTemplate, imSearch)
    [M, N] = size(imTemplate);
    [R, C] = size(imSearch);

    %// Cast to double for precision
    imTemplate = im2double(imTemplate);
    imSearch = im2double(imSearch);

    neigh = im2col(imSearch, [M, N]);
    templateCol = imTemplate(:); %// Ensures we place template into single column

    imInvert = ~imTemplate;
    imInvertNoBorder = imclearborder(imInvert, 8); %// Search 8-pixel neighbourhood
    rowsToRemove = imInvertNoBorder(:) == 1;
    neigh(rowsToRemove,:) = [];

    neighMeanSubtract = bsxfun(@minus, neigh, mean(neigh));
    templateMeanSubtract = templateCol - mean(templateCol);

    numerator = bsxfun(@times, neighMeanSubtract, templateMeanSubtract);
    sumNumerator = sum(numerator);

    denominator1 = sqrt(sum(neighMeanSubtract.*neighMeanSubtract));
    denominator2 = sqrt(sum(templateMeanSubtract.*templateMeanSubtract));
    sumDenominator = denominator1 .* denominator2;

    NCC = sumNumerator ./ sumDenominator;    

    finalOutput = col2im(NCC, [M, N], [R, C]);
    finalOutput(isnan(finalOutput)) = 0;

    maxCoeff = max(abs(finalOutput(:)));
    [rowNCC, colNCC] = find(abs(finalOutput) == maxCoeff);
end

Удачи!