function GenerateHoughClusters(list_of_matches) % This function takes list of matches as input and computes outliers / % moving objects / background; exports to a file. % First, some constants % Label for moving objects MOVING_OBJECT = 1; % Label for outliers OUTLIER = -1; % Label for background BACKGROUND = 0; % A point must be at most this far away to be part of a cluster NEIGHBORHOOD_THRESHOLD = 15; % Determines whether to graph outliers or not (red) SHOW_OUTLIERS = 0; % It is possible that the plane is tracking the moving object. % In this case the tiepoints will be so close to each other, % the angles won't make any sense. We can detect this by looking % at the line lengths and if they are too smal, we set the angle to 0 LINE_LENGTH_THRESHOLD = 5; %pixels % columns location of image locations in matches locs = [2,3,7,8]; % **************************************************************** % This is the key data structure containing the whole scene % The column format is the following: % x_ref y_ref x_tar y_tar angle distance label scene = zeros(0,8); prevousScene = scene; GlobalScene = scene; % **************************************************************** % Structure of scene refXIndex = 1; refYIndex = 2; tarXIndex = 3; tarYIndex = 4; angleIndex = 5; distanceIndex = 6; labelIndex = 7; timeIndex = 8; % Create the figure and label axes %HoughFigure = figure('Position', [100 100 1024 768],'PaperPosition',[0 0 30 20]); %xlabel('Line angle'); %ylabel('Line length'); %zlabel('Time'); %title('Hough space'); %hold on; % Open the matches list hMap = fopen(list_of_matches, 'r'); fprintf('Reading matches map file\n'); time = 0; label = 1; changeLabel = 0; % Read the list of match files while(~feof(hMap)) fprintf('File # = %d\n', time); % Save the scene so that we can match the next frame with the % clusters from the current frame previousScene = scene; scene = zeros(0,8); matchFile = fgetl(hMap); if (exist(matchFile)==0) fprintf('File # %d does not exist\n', time); time = time +1;pre continue; end % Read in matches matches = load(matchFile); matches = RemoveDuplicates(matches); numMatches= size(matches, 1); % Calculate angles and distances, label all points as "background" (0) for j = 1:numMatches ptRef = [matches(j, locs(1)), matches(j, locs(2))]; ptTar = [matches(j, locs(3)), matches(j, locs(4))]; distance = Distance(ptRef, ptTar); if (distance < LINE_LENGTH_THRESHOLD) angle = 0; else angle = Angle(ptRef, ptTar); end %if (time == 0) % label everything with "0" for now - as background % label = label + 1; % changeLabel = label; %else % changeLabel = FindLabel(ptRef, previousScene); % if (changeLabel == -1) % label = label + 1; % changeLabel = label; % end %end %scene = vertcat(scene, [ptRef, ptTar, angle, distance, changeLabel, time]); scene = vertcat(scene, [ptRef, ptTar, angle, distance, BACKGROUND, time]); end % calculating angles and distances % Calculate median distance and slope (hypothetical center of background cluster) %scene = sortrows(scene, angleIndex); %medianAngle = median(scene(:, angleIndex)); %scene = sortrows(scene, distanceIndex); %medianDistance = median(scene(:, distanceIndex)); % Try to find all Hough clusters, starting with the first point in % the list % Each point could be its own cluster possibleCluster = 1; % Number of points left in the scene (will change) numMatchesLeft = numMatches; while (possibleCluster <= numMatchesLeft) % For now, the first point in the cluster is its own neighbor neighbors = scene(possibleCluster, :); % The first point in the cluster is the seed of the cluster seed = scene(possibleCluster, angleIndex:distanceIndex); % Remove the seed point from the scene scene(possibleCluster,:) = []; %### % Find the neighbors of the seed, recursively - find N neighbors % near the seed within the circle defined by NEIGBORHOOD_THRESHOLD, % then for each of those neighbors find its neighbors within their % individual circles % % As the neighbors are found, they are removed from scene [neighbors scene] = FindNeighbors(neighbors, scene, seed, NEIGHBORHOOD_THRESHOLD); % What is the total number of matches now? (could be numMatches-1) numMatchesLeft = size(scene,1); % Store all the neighbors in a separate data structure - now a full % cluster. clusterContainer = neighbors; % Check how many neighbors were found clusterCount = size(clusterContainer, 1); % ################################################################# % Now we have a cluster of points to plot - here we determine % whether it's the background, a moving object, or an outlier % ################################################################# % If there are more than 2 points but less than 20 percent of % total, it's a moving object, color it BLUE if (clusterCount > 2 && clusterCount < round(0.2*numMatches)) clusterContainer(:, labelIndex) = MOVING_OBJECT * ones(clusterCount, 1); %set(0,'CurrentFigure', HoughFigure); %for i = 1:clusterCount %color = ones(1,3) * clusterContainer(i, labelIndex); %plot3(clusterContainer(i, angleIndex), clusterContainer(i, distanceIndex), ones(1,1)*time, 'Color', color, 'Marker', '.', 'MarkerSize', 5, 'LineStyle', 'none'); %end %plot3(clusterContainer(:, angleIndex), clusterContainer(:, distanceIndex), ones(clusterCount,1)*time, 'Color', color, 'Marker', '.', 'MarkerSize', 5, 'LineStyle', 'none'); else % If there are less than 3 percent, these are % outliers/crossovers, they are RED if (clusterCount < round(0.03*numMatches)) % But FIRST, check the previous frame, if it exists, % whether there are any neighbors around these points % in the previous frame if (time ~= 0) % clusterContainer contains the 'outliers' % previousScene is the previous frame % for each 'outlier' point in clusterContainer we % need to look for possible neighbors in the % previous frame - if at least one exists, then the % neighbor in the previous frame and the current % 'outlier' are not outliers but moving objects % instead numOutliers = size(clusterContainer, 1); for outlier = 1:numOutliers outlierSeed = clusterContainer(outlier, :); outlierNeighbors = zeros(0,8); % since the network is modified (neighbors are % removed), we need to save it and use a fresh % copy for each outlier - it's possible two % outliers may share the same neighbor in the % previous frame outlierNetwork = previousScene; % find neighbors in the previous frame [outlierNeighbors outlierNetwork] = FindNeighbors(outlierNeighbors, outlierNetwork, outlierSeed, NEIGHBORHOOD_THRESHOLD); % how many neighbors for the outlier were found? numOutlierNeighborsInPrevFrame = size(outlierNeighbors,1); if (numOutlierNeighborsInPrevFrame > 0) % If neighbors were found, current % 'outlier' is not an outlier clusterContainer(outlier, labelIndex) = MOVING_OBJECT;% * ones(clusterCount, 1); % The neighbors found in the previous scene are not outliers either for numPointsToEdit = 1:numOutlierNeighborsInPrevFrame % Remove previous frame from global scene % GlobalScene(find(GlobalScene(:, timeIndex) = time),:) = []; GlobalScene(:, timeIndex == time) = []; % Change found neighbors to MOVING_OBJECT outlierNeighbors(:, labelIndex) = MOVING_OBJECT * ones(numOutlierNeighborsInPrevFrame, 1); % Put neighbors and previous frame back into global scene GlobalScene = vertcat(GlobalScene, outlierNeighbors); GlobalScene = vertcat(GlobalScene, outlierNetwork); clusterContainer(outlier, labelIndex) = OUTLIER;% * ones(clusterCount, 1); end % for else clusterContainer(outlier, labelIndex) = OUTLIER;% * ones(clusterCount, 1); end end %for outlier %color = 'red'; %if (SHOW_OUTLIERS) %set(0,'CurrentFigure', HoughFigure); %plot3(clusterContainer(:, angleIndex), clusterContainer(:, distanceIndex), ones(clusterCount,1)*time, 'Color', color, 'Marker', '.', 'MarkerSize', 5, 'LineStyle', 'none'); %end %else % If there are more than 20% of the points in the cluster, % then the cluster is background, the points are GREEN %color = 'green'; %set(0,'CurrentFigure', HoughFigure); %plot3(clusterContainer(:, angleIndex), clusterContainer(:, distanceIndex), ones(clusterCount,1)*time, 'Color', color, 'Marker', '.', 'MarkerSize', 5, 'LineStyle', 'none'); end % if time != 0 end % if (outlier) end % else GlobalScene = vertcat(GlobalScene, clusterContainer); end % while possible cluster >=numMatches time = time + 1; %pause; end % end reading in matches files save GlobalScene.mat GlobalScene % These functions are defined within showMatches, not outside %--------------------------------------------------------------------------------------------- % Calculates angle between two points function angle = Angle(ptA, ptB) %if (ptA(2)-ptB(2) == 0) % angle = -1; % return %end angle = atan2((ptA(2)-ptB(2)),(ptA(1)-ptB(1))); angle = 180/pi * angle; end %--------------------------------------------------------------------------------------------- % Calculates distance between two points % ptB could be a vector, then distance is a vector of distances function distance = Distance(ptA, ptB) % if pointB is a vector [rows cols] = size(ptB); if (rows > 1) for i = 1:rows distance(i) = sqrt(sum((ptA - ptB(i,:)).^2)); end else distance = sqrt(sum((ptA - ptB).^2)); end end %--------------------------------------------------------------------------------------------- % Export current figure to file function ExportFigure(matchesFigure, HoughFigure, matchPath) % extract filename [pathstr, name, ext, versn] = fileparts(matchPath); % construct new filename with new extension matchesimg=fullfile('matches_', name, '.jpg'); % remove backslash from name - matlab 'feature' matchesimg(matchesimg=='\') = []; matchesimg(matchesimg=='/') = []; % construct new filename with new extension houghimg=fullfile('hough_', name, '.jpg'); % remove backslash from name - matlab 'feature' houghimg(houghimg=='\') = []; houghimg(houghimg=='/') = []; % export to jpeg set(0,'CurrentFigure', matchesFigure); print('-djpeg', '-r300', matchesimg); set(o,'CurrentFigure', HoughFigure); print('-djpeg', '-r300', houghimg); %print('-dtiff', 'SIFT.tif') end %--------------------------------------------------------------------------------------------- %----------------------------------------------------------------------------------------------------------------------------------- % The recursive function which, given a seed point, recursively finds all neighbors and their neighbors and so on % in the neighborhood of the original seed within a threshold % Each of the neighbors found around the seed becomes a new seed for % the recursion % neighbors: as input, 1x8 vector containing the original seed point so that it doesn't get lost % network: points we're looking at % seed: where are we starting to look % distance: threshold around the seed to find first neighbors % returns neighbors and modified network (not containing the neighbors) function [neighbors network] = FindNeighbors(neighbors, network, seed, distance) if (size(network, 1) == 0) % neighbors = []; return end % Compute distances from the seed to each point in the network distances = abs(Distance(seed, network(:, angleIndex:distanceIndex))); % Select the indices of the network points which are within the % circle neighborIndices = find(distances < distance); % Termination step. Function terminates when there are no neighbors. numNeighbors = size(neighborIndices, 2); if (numNeighbors == 0) %neighbors = []; return end % Add newly found points to the neighbors matrix neighbors = vertcat(neighbors, network(neighborIndices,:)); % Remove newly found points from the network network(neighborIndices,:) = []; % Recursive step: for each of the newly found neighbors, find neighbors within the threshold % Function terminates when there are no neighbors %recursion = 1; %while( recursion < numNeighbors ) for recursion = 1:numNeighbors newNeighbors = []; newSeed = neighbors(recursion, angleIndex:distanceIndex); %newNeighbors = neighbors(recursion, :); [newNeighbors network]= FindNeighbors(newNeighbors, network, newSeed, distance); if (size(newNeighbors,2) > 0) neighbors = vertcat(neighbors, newNeighbors); end end end %----------------------------------------------------------------------------------------------------------------------------------- % Removes duplicate SIFT matches function frame = RemoveDuplicates(frame) [numPts cols] = size(frame); i = 1; j = 1; while (i < numPts ) j = i+1; while (j < numPts) if ((frame(i,1) == frame(j,1) && frame(i,2) == frame(j,2)) || (frame(i,3) == frame(j,3) && frame(i,4) == frame(j,4))) frame(i,:) = []; numPts = numPts - 1; end j = j + 1; end i = i + 1; end return end %----------------------------------------------------------------------------------------------------------------------------------- function label = FindLabel(ptRef, previousScene) distances = Distance(ptRef, previousScene(:, tarXIndex:tarYIndex)); [rows cols] = size(distances); distances = reshape(distances, [cols rows]); distances = horzcat(distances, previousScene(:, labelIndex)); distances = sortrows(distances, 1); if (distances(1,1) < NEIGHBORHOOD_THRESHOLD) label = distances(1, 2); else label = -1; end end %----------------------------------------------------------------------------------------------------------------------------------- end % end GenerateHoughClusters function