gusucode.com > vision工具箱matlab源码程序 > vision/cvalgMatchFeatures.m
function [indexPairs, matchMetric] = cvalgMatchFeatures(features1in, ... features2in, metric, ... matchPercentage, method, ... maxRatioThreshold, isPrenormalized, ... uniqueMatches, isLegacyMethod) % Main algorithm used by the matchFeatures function. See matchFeatures % for more details. % Copyright 2011 The MathWorks, Inc. %#codegen % Determine output class if (isa(features1in, 'double')) outputClass = 'double'; else outputClass = 'single'; end if isempty(features1in) || isempty(features2in) indexPairs = zeros(2, 0, 'uint32'); matchMetric = zeros(1, 0, outputClass); return; end % cast feature data to expected output class [features1, features2] = castFeatures(features1in, features2in, metric,... method, outputClass); % normalize features using L2-norm if ~isPrenormalized && ~strcmpi(metric, 'hamming') [features1, features2] = normalizeFeatures(features1, features2, method, metric); end % Convert match threshold percent to a numeric threshold matchThreshold = percentToLevel(matchPercentage, size(features1, 1), ... metric, outputClass); % Find matches based on selected method if isLegacyMethod [indexPairs, matchMetric] = findMatchesLegacy(features1, features2, ... metric, method, maxRatioThreshold, matchThreshold, outputClass); elseif strcmpi(method, 'approximate') [indexPairs, matchMetric] = findMatchesApproximate(features1,features2,... metric, maxRatioThreshold, matchThreshold, uniqueMatches, outputClass); else % exhaustive [indexPairs, matchMetric] = findMatchesExhaustive(features1,features2, ... metric, maxRatioThreshold, matchThreshold, uniqueMatches, outputClass); end %========================================================================== % Use Approximate nearest neighbor search to find matches. %========================================================================== function [indexPairs, matchMetric] = findMatchesApproximate(features1,features2,... metric, maxRatioThreshold, matchThreshold, uniqueMatches, outputClass) N2 = cast(size(features2, 2), 'uint32'); [indexPairs, matchMetric] = findApproximateNearestNeighbors(features1,... features2, metric, outputClass); [indexPairs, matchMetric] = removeWeakMatches(indexPairs, ... matchMetric, matchThreshold, metric); [indexPairs, matchMetric] = removeAmbiguousMatches(indexPairs, ... matchMetric, maxRatioThreshold, N2, metric); if uniqueMatches % perform backward match to remove non-unique matches. % avoid search for the same feature more than once. idx = unique(indexPairs(2,:)); reverse_index_pairs = findApproximateNearestNeighbors(... features2(:,idx), features1, metric, outputClass); f2ToF1Matches = zeros(1,N2,'uint32'); f2ToF1Matches(idx) = reverse_index_pairs(2,:); symmetric_pairs = f2ToF1Matches(indexPairs(2,:)) == indexPairs(1,:); else % branch required for codegen symmetric_pairs = true(1,size(indexPairs,2)); end indexPairs = indexPairs(:,symmetric_pairs); matchMetric = matchMetric(1,symmetric_pairs); %========================================================================== % Find approximate nearest neighbors. Return indices to the matching % features and the distance metrics for the 2 nearest neighbors (for the % ratio test). %========================================================================== function [indexPairs, metrics] = findApproximateNearestNeighbors(... features1, features2, metric, outputClass) if isSimMode [indexPairs, metrics] = ocvFlannBasedMatching(... features1, features2, metric); else coder.internal.errorIf(~coder.internal.isTargetMATLABHost,... 'vision:matchFeatures:codegenApproxHostOnly'); [indexPairs, metrics] = ... vision.internal.buildable.matchFeaturesApproxNN.findApproximateNearestNeighbors(... features1, features2, metric); end % type of metric output should match feature input type metrics = cast(metrics, outputClass); % convert to 1-based indexing & cast to uint32 indexPairs = bsxfun(@plus, uint32(indexPairs(1,:)), uint32(1)); indexPairs = vertcat(uint32(1:size(indexPairs,2)), indexPairs(1,:)); %========================================================================== % Find matches using an exhaustive search. %========================================================================== function [indexPairs, matchMetric] = findMatchesExhaustive(features1, features2,... metric, maxRatioThreshold, matchThreshold, uniqueMatches, outputClass) % SCORES is an N1-by-N2 correspondence metric matrix where the rows % correspond to the feature vectors in FEATURES1, and the columns % correspond to the feature vectors in FEATURES2. N1 = uint32(size(features1,2)); N2 = uint32(size(features2,2)); scores = exhaustiveDistanceMetrics(features1, features2, N1, N2, outputClass, metric); [indexPairs, matchMetric] = findNearestNeighbors(scores, metric); [indexPairs, matchMetric] = removeWeakMatches(indexPairs, ... matchMetric, matchThreshold, metric); [indexPairs, matchMetric] = removeAmbiguousMatches(indexPairs, ... matchMetric, maxRatioThreshold, N2, metric); if uniqueMatches uniqueIndices = findUniqueIndices(scores, metric, indexPairs); else % branch required for codegen uniqueIndices = true(1,size(indexPairs,2)); end indexPairs = indexPairs(:, uniqueIndices); matchMetric = matchMetric(1, uniqueIndices); %========================================================================== % Find nearest neighbors using an exhaustive search. Return indices to the % matching features and the distance metrics for the 2 nearest neighbors % (for the ratio test). %========================================================================== function [indexPairs, topTwoMetrics] = findNearestNeighbors(scores, metric) if strcmp(metric, 'normxcorr') [topTwoMetrics, topTwoIndices] = vision.internal.partialSort(scores, 2, 'descend'); else [topTwoMetrics, topTwoIndices] = vision.internal.partialSort(scores, 2, 'ascend'); end indexPairs = vertcat(uint32(1:size(scores,1)), topTwoIndices(1,:)); %========================================================================== % Compute distance metrics %========================================================================== function scores = exhaustiveDistanceMetrics(features1, features2, ... N1, N2, outputClass, metric) % Generate correspondence metric matrix switch metric case 'sad' % Generate correspondence matrix using Sum of Absolute Differences scores = metricSAD(features1, features2, N1, N2, outputClass); case 'normxcorr' % Generate correspondence matrix using Normalized Cross-correlation scores = metricNormXCorr(features1, features2); case 'ssd' % Generate correspondence matrix using Sum of Squared Differences scores = metricSSD(features1, features2, N1, N2, outputClass); otherwise % 'hamming' % Generate correspondence matrix using Sum of Squared Differences scores = metricHamming(features1, features2, N1, N2, outputClass); end %========================================================================== % Find matches using legacy method modes. %========================================================================== function [indexPairs, matchMetric] = findMatchesLegacy(features1, features2, ... metric, method, maxRatioThreshold, matchThreshold, outputClass) N1 = uint32(size(features1,2)); N2 = uint32(size(features2,2)); switch method case 'nearestneighborsymmetric' scores = exhaustiveDistanceMetrics(features1, features2, N1, N2, ... outputClass, metric); % keep this for backward compatibility [indexPairs, matchMetric] = findMatchesNN(scores, metric, ... matchThreshold); case 'threshold' % keep this for backward compatibility scores = exhaustiveDistanceMetrics(features1, features2, N1, N2, ... outputClass, metric); [indexPairs, matchMetric] = findMatchesThreshold(scores, ... metric, matchThreshold); case 'nearestneighbor_old' % keep this for backward compatibility scores = exhaustiveDistanceMetrics(features1, features2, N1, N2,... outputClass, metric); [indexPairs, matchMetric] = findMatchesNN_old(scores, N1, N2, ... metric); [indexPairs, matchMetric] = removeWeakMatches(indexPairs, ... matchMetric, matchThreshold, metric); case 'nearestneighborratio' if N2 > 1 [indexPairs, matchMetric] = findMatchesExhaustive(... features1, features2, metric, maxRatioThreshold, ... matchThreshold, false, outputClass); else % If FEATURES2 contains only 1 feature, we cannot use ratio. % Use NearestNeighborSymmetric instead, resulting in a single % match scores = exhaustiveDistanceMetrics(features1, features2, ... N1, N2, outputClass, metric); [indexPairs, matchMetric] = findMatchesNN(scores, ... metric, matchThreshold); end end %========================================================================== % Remove ambiguous matches using the nearest neighbor ratio test. %========================================================================== function [indexPairs, matchMetric] = removeAmbiguousMatches(indexPairs, ... matchMetric, maxRatio, N2, metric) if N2 > 1 % apply ratio test only if there are more than 1 feature vectors unambiguousIndices = findUnambiguousMatches(matchMetric, maxRatio, metric); else unambiguousIndices = true(1,size(matchMetric,2)); end indexPairs = indexPairs(:, unambiguousIndices); matchMetric = matchMetric(1,unambiguousIndices); %========================================================================== % Enforce uniqueness by applying a bi-directional match constraint %========================================================================== function uniqueIndices = findUniqueIndices(scores, metric, ... indexPairs) if strcmpi(metric,'normxcorr') [~, idx] = max(scores(:,indexPairs(2,:))); else [~, idx] = min(scores(:,indexPairs(2,:))); end uniqueIndices = idx == indexPairs(1,:); %========================================================================== % Find matches using Nearest-Neighbor strategy %========================================================================== function [indexPairs, matchMetric] = findMatchesNN(scores, ... metric, matchThreshold) % Find the maximum(minimum) entry in scores. % Make it a match. % Eliminate the corresponding row and the column. % Repeat. nRows = size(scores, 1); nCols = size(scores, 2); nMatches = min(nRows, nCols); indexPairs = zeros([2, nMatches], 'uint32'); matchMetric = zeros(1, nMatches, 'like', scores); useMax = strcmp(metric, 'normxcorr'); for i = 1:nMatches if useMax [matchMetric(i), ind] = max(scores(:)); [r, c] = ind2sub(size(scores), ind); else [matchMetric(i), ind] = min(scores(:)); [r, c] = ind2sub(size(scores), ind); end indexPairs(:, i) = [r, c]; if useMax scores(r, :) = -inf('like',scores); scores(:, c) = -inf('like',scores); else scores(r, :) = inf('like',scores); scores(:, c) = inf('like',scores); end end [indexPairs, matchMetric] = removeWeakMatches(indexPairs, ... matchMetric, matchThreshold, metric); %========================================================================== % Normalize features to be unit vectors %========================================================================== function [features1, features2] = normalizeFeatures(features1, features2, ... method, metric) % move this to parsing where we map old method values. % normalize the features if strcmp(method, 'nearestneighbor_old') && ... strcmp(metric, 'normxcorr') % for backward compatibility, subtract the mean from features f1Mean = mean(features1); features1 = bsxfun(@minus, features1, f1Mean); f2Mean = mean(features2); features2 = bsxfun(@minus, features2, f2Mean); end % Convert feature vectors to unit vectors features1 = normalizeX(features1); features2 = normalizeX(features2); %========================================================================== % Convert match threshold percent to an numeric threshold %========================================================================== function matchThreshold = percentToLevel(matchPercentage, ... vector_length, metric, outputClass) matchPercentage = cast(matchPercentage, outputClass); vector_length = cast(vector_length, outputClass); if (strcmp(metric, 'normxcorr')) matchThreshold = cast(0.01, outputClass)*(cast(100, outputClass) ... - matchPercentage); else if (strcmp(metric, 'sad')) max_val = cast(2, outputClass)*sqrt(vector_length); elseif (strcmp(metric, 'ssd')) max_val = cast(4, outputClass); else % 'hamming' % the below value assumes that binary features are stored % in 8-bit buckets, which is correct for binaryFeatures class max_val = cast(8*vector_length, outputClass); end matchThreshold = (matchPercentage*cast(0.01, outputClass))*max_val; if strcmp(metric, 'hamming') % Round up since we are dealing with whole bits matchThreshold = round(matchThreshold); end end %========================================================================== % Cast features to expected output class. %========================================================================== function [features1, features2] = castFeatures(features1in, features2in,... metric, method, outputClass) if ~strcmp(metric, 'hamming') if strcmpi(method, 'approximate') % approximate search requires single features1 = cast(features1in, 'single'); features2 = cast(features2in, 'single'); else features1 = cast(features1in, outputClass); features2 = cast(features2in, outputClass); end else % do not cast binary feature data features1 = features1in; features2 = features2in; end %========================================================================== % Find unambiguous matches using David Lowe's disambiguation strategy %========================================================================== function unambiguousIndices = findUnambiguousMatches(topTwoScores, maxRatioThreshold, metric) if strcmpi(metric, 'normxcorr') % If the metric is 'normxcorr', then the scores are cosines % of the angles between the feature vectors. % The ratio of the angles is an approximation of the ratio of % euclidean distances. See David Lowe's demo code. topTwoScores(topTwoScores > 1) = 1; % prevent complex angles topTwoScores = acos(topTwoScores); end % handle division by effective zero zeroInds = topTwoScores(2, :) < cast(1e-6, 'like', topTwoScores); topTwoScores(:, zeroInds) = 1; ratios = topTwoScores(1, :) ./ topTwoScores(2, :); unambiguousIndices = ratios <= maxRatioThreshold; %========================================================================== % Find matches using a bidirectional greedy strategy %========================================================================== function [indexPairs, matchMetric] = findMatchesThreshold(scores, ... metric, matchThreshold) if strcmp(metric, 'normxcorr') inds = find(scores >= matchThreshold); else inds = find(scores <= matchThreshold); end matchMetric = scores(inds)'; [rowInds, colInds] = ind2sub(size(scores), inds); indexPairs = cast([rowInds, colInds]', 'uint32'); %========================================================================== % Find matches using the old nearest neighbor strategy for compatibility %========================================================================== function [indexPairs, matchMetric] = findMatchesNN_old(scores, N1, N2,... metric) % SCORES is an N1-by-N2 correspondence metric matrix where the rows % correspond to the feature vectors in FEATURES1, and the columns % correspond to the feature vectors in FEATURES2. % For each feature vector in FEATURES1 find the best match in FEATURES2. % This is simply the entry along each row with the minimum or maximum % metric value, depending on the metric. row_index_pairs = [(1:N1); zeros(1, N1)]; if (strcmp(metric, 'normxcorr')) [row_scores, row_index_pairs(2, :)] = max(scores, [], 2); else [row_scores, row_index_pairs(2, :)] = min(scores, [], 2); end % For each feature vector in FEATURES2 find the best match in FEATURES1. % This is simply the entry down each col with the minimum or maximum metric % value, depending on the metric. col_index_pairs = [zeros(1, N2); (1:N2)]; if (strcmp(metric, 'normxcorr')) [col_scores, col_index_pairs(1, :)] = max(scores, [], 1); else [col_scores, col_index_pairs(1, :)] = min(scores, [], 1); end % Concatenate both row and column lists M = [row_scores', col_scores]; indexPairs = [row_index_pairs, col_index_pairs]; % Remove duplicate entries in the matches list [trimmed_index_pairs, I, ~] = unique(indexPairs', 'rows'); matchMetric = M(I); indexPairs = trimmed_index_pairs'; %========================================================================== % Generate correspondence metric matrix using Sum of Absolute Differences %========================================================================== function scores = metricSAD(features1, features2, N1, N2, outputClass) % Need to obtain feature vector length to perform explicit row indexing % needed for code generation of variable sized inputs if isSimMode() % call optimized builtin function scores = zeros(N1, N2, outputClass); scores(:,:) = visionSADMetric(features1,features2); else if coder.internal.isTargetMATLABHost scores = vision.internal.buildable.ComputeMetricBuildable.ComputeMetric_core(features1,features2, 'sad', N1, N2, outputClass); else % for portable C code generation vector_length = size(features1, 1); scores = zeros(N1, N2, outputClass); for c = 1:N2 for r = 1:N1 scores(r, c) = sum(abs(features1(1:vector_length, r) - ... features2(1:vector_length, c))); end end end end %========================================================================== % Generate correspondence metric matrix using Sum of Squared Differences %========================================================================== function scores = metricSSD(features1, features2, N1, N2, outputClass) % Need to obtain feature vector length to perform explicit row indexing % needed for code generation of variable sized inputs if isSimMode() % call optimized builtin function scores = zeros(N1, N2, outputClass); scores(:,:) = visionSSDMetric(features1,features2); else if coder.internal.isTargetMATLABHost scores = vision.internal.buildable.ComputeMetricBuildable.ComputeMetric_core(features1,features2, 'ssd', N1, N2, outputClass); else % for portable C code generation vector_length = size(features1, 1); scores = zeros(N1, N2, outputClass); for c = 1:N2 for r = 1:N1 scores(r, c) = sum((features1(1:vector_length, r) - ... features2(1:vector_length, c)).^2); end end end end %========================================================================== % Generate correspondence metric matrix based on Hamming distance %========================================================================== function scores = metricHamming(features1, features2, N1, N2, outputClass) persistent lookupTable; % lookup table for counting bits % Need to obtain feature vector length to perform explicit row indexing % needed for code generation of variable sized inputs if isSimMode() % call optimized builtin function scores = zeros(N1, N2, outputClass); scores(:,:) = visionHammingMetric(features1,features2); else if coder.internal.isTargetMATLABHost scores = vision.internal.buildable.ComputeMetricBuildable.ComputeMetric_core(features1,features2, 'hamming', N1, N2, outputClass); else % for portable C code generation vector_length = size(features1, 1); scores = zeros(N1, N2, outputClass); if isempty(lookupTable) lookupTable = zeros(256, 1, outputClass); for i = 0:255 lookupTable(i+1) = sum(dec2bin(i)-'0'); end end for c = 1:N2 for r = 1:N1 temp = bitxor(features1(1:vector_length, r),... features2(1:vector_length, c)); idx = double(temp) + 1; % cast needed to avoid integer math scores(r,c) = sum(lookupTable(idx)); end end end end %========================================================================== % Generate correspondence metric matrix using Normalized Cross-Correlation %========================================================================== function scores = metricNormXCorr(features1, features2) scores = features1' * features2; %========================================================================== % Remove weak matches %========================================================================== function [indices, matchMetric] = removeWeakMatches(indices, ... matchMetric, matchThreshold, metric) if (strcmp(metric, 'normxcorr')) inds = matchMetric(1,:) >= matchThreshold; else inds = matchMetric(1,:) <= matchThreshold; end indices = indices(:, inds); matchMetric = matchMetric(:, inds); %========================================================================== % Normalize the columns in X to have unit norm. %========================================================================== function X = normalizeX(X) Xnorm = sqrt(sum(X.^2, 1)); X = bsxfun(@rdivide, X, Xnorm); % Set effective zero length vectors to zero X(:, (Xnorm <= eps(single(1))) ) = 0; %========================================================================== function flag = isSimMode() flag = isempty(coder.target);