% EDGELINK - Link edge points in an image into lists % % Usage: [edgelist edgeim, etype] = edgelink(im, minlength, location) % % **Warning** 'minlength' is ignored at the moment because 'cleanedgelist' % has some bugs and can be memory hungry % % Arguments: im - Binary edge image, it is assumed that edges % have been thinned (or are nearly thin). % minlength - Optional minimum edge length of interest, defaults % to 1 if omitted or specified as []. Ignored at the % moment. % location - Optional complex valued image holding subpixel % locations of edge points. For any pixel the % real part holds the subpixel row coordinate of % that edge point and the imaginary part holds % the column coordinate. See NONMAXSUP. If % this argument is supplied the edgelists will % be formed from the subpixel coordinates, % otherwise the the integer pixel coordinates of % points in 'im' are used. % % Returns: edgelist - a cell array of edge lists in row,column coords in % the form % { [r1 c1 [r1 c1 etc } % r2 c2 ... % ... % rN cN] ....] % % edgeim - Image with pixels labeled with edge number. % Note that junctions in the labeled edge image will be % labeled with the edge number of the last edge that was % tracked through it. Note that this image also includes % edges that do not meet the minimum length specification. % If you want to see just the edges that meet the % specification you should pass the edgelist to % DRAWEDGELIST. % % etype - Array of values, one for each edge segment indicating % its type % 0 - Start free, end free % 1 - Start free, end junction % 2 - Start junction, end free (should not happen) % 3 - Start junction, end junction % 4 - Loop % % This function links edge points together into lists of coordinate pairs. % Where an edge junction is encountered the list is terminated and a separate % list is generated for each of the branches. % % Note I am not sure if using bwmorph's 'thin' or 'skel' is best for % preprocessing the edge image prior to edgelinking. The main issue is the % treatment of junctions. Skel can result in an image where multiple adjacent % junctions are produced (maybe this is more a problem with my junction % detection code). Thin, on the other hand, can produce different output when % you rotate an image by 90 degrees. On balance I think using 'thin' is better. % Note, however, the input image should be 'nearly thin' otherwise the thinning % operation could shorten the ends of structures. Skeletonisation and thinning % is surprisingly awkward. % % See also: DRAWEDGELIST, LINESEG, MAXLINEDEV, CLEANEDGELIST, % FINDENDSJUNCTIONS, FILLEDGEGAPS % Copyright (c) 1996-2017 Peter Kovesi % Centre for Exploration Targeting % The University of Western Australia % peter.kovesi at uwa edu au % % Permission is hereby granted, free of charge, to any person obtaining a copy % of this software and associated documentation files (the "Software"), to deal % in the Software without restriction, subject to the following conditions: % % The above copyright notice and this permission notice shall be included in % all copies or substantial portions of the Software. % % The Software is provided "as is", without warranty of any kind. % February 2001 - Original version % September 2004 - Revised to allow subpixel edge data to be used % November 2006 - Changed so that edgelists start and stop at every junction % January 2007 - Trackedge modified to discard isolated pixels and the % problems they cause (thanks to Jeff Copeland) % January 2007 - Fixed so that closed loops are closed! % May 2013 - Completely redesigned with a new linking strategy that % hopefully handles adjacent junctions correctly. It runs % about twice as fast too. % November 2017 - Distance bugfix in trackedge function [edgelist, edgeim, etype] = edgelink(im, minlength, location) % Set up some global variables to avoid passing (and copying) of arguments, % this improves speed. global EDGEIM; global ROWS; global COLS; global JUNCT; if ~exist('minlength','var') || isempty(minlength), minlength = 0; end EDGEIM = im ~= 0; % Make sure image is binary. EDGEIM = bwmorph(EDGEIM,'clean'); % Remove isolated pixels % Make sure edges are thinned. Use 'thin' rather than 'skel', see % comments in header. EDGEIM = bwmorph(EDGEIM,'thin',Inf); [ROWS, COLS] = size(EDGEIM); % Find endings and junctions in edge data [RJ, CJ, re, ce] = findendsjunctions(EDGEIM); Njunct = length(RJ); Nends = length(re); % Create a sparse matrix to mark junction locations. This makes junction % testing much faster. A value of 1 indicates a junction, a value of 2 % indicates we have visited the junction. JUNCT = spalloc(ROWS,COLS, Njunct); for n = 1:Njunct JUNCT(RJ(n),CJ(n)) = 1; end % ? Think about using labels >= 2 so that EDGEIM can be uint16, say. % EDGEIM = double(EDGEIM); % Cast to double to allow the use of -ve labelings edgeNo = 0; % Summary of strategy: % 1) From every end point track until we encounter an end point or % junction. As we track points along an edge image pixels are labeled with % the -ve of their edge No. % 2) From every junction track out on any edges that have not been % labeled yet. % 3) Scan through the image looking for any unlabeled pixels. These % correspond to isolated loops that have no junctions. %% 1) Form tracks from each unlabeled endpoint until we encounter another % endpoint or junction. for n = 1:Nends if EDGEIM(re(n),ce(n)) == 1 % Endpoint is unlabeled edgeNo = edgeNo + 1; [edgelist{edgeNo} endType] = trackedge(re(n), ce(n), edgeNo); etype(edgeNo) = endType; end end %% 2) Handle junctions. % Junctions are awkward when they are adjacent to other junctions. We % start by looking at all the neighbours of a junction. % If there is an adjacent junction we first create a 2-element edgetrack % that links the two junctions together. We then look to see if there are % any non-junction edge pixels that are adjacent to both junctions. We then % test to see which of the two junctions is closest to this common pixel and % initiate an edge track from the closest of the two junctions through this % pixel. When we do this we set the 'avoidJunction' flag in the call to % trackedge so that the edge track does not immediately loop back and % terminate on the other adjacent junction. % Having checked all the common neighbours of both junctions we then % track out on any remaining untracked neighbours of the junction for j = 1:Njunct if JUNCT(RJ(j),CJ(j)) ~= 2; % We have not visited this junction JUNCT(RJ(j),CJ(j)) = 2; % Call availablepixels with edgeNo = 0 so that we get a list of % available neighbouring pixels that can be linked to and a list of % all neighbouring pixels that are also junctions. [ra, ca, rj, cj] = availablepixels(RJ(j), CJ(j), 0); for k = 1:length(rj) % For all adjacent junctions... % Create a 2-element edgetrack to each adjacent junction edgeNo = edgeNo + 1; edgelist{edgeNo} = [RJ(j) CJ(j); rj(k) cj(k)]; etype(edgeNo) = 3; % Edge segment is junction-junction EDGEIM(RJ(j), CJ(j)) = -edgeNo; EDGEIM(rj(k), cj(k)) = -edgeNo; % Check if the adjacent junction has some untracked pixels that % are also adjacent to the initial junction. Thus we need to % get available pixels adjacent to junction (rj(k) cj(k)) [rak, cak] = availablepixels(rj(k), cj(k)); % If both junctions have untracked neighbours that need checking... if ~isempty(ra) && ~isempty(rak) % Find untracked neighbours common to both junctions. commonrc = intersect([ra ca], [rak cak], 'rows'); for n = 1:size(commonrc, 1); % If one of the junctions j or k is closer to this common % neighbour use that as the start of the edge track and the % common neighbour as the 2nd element. When we call % trackedge we set the avoidJunction flag to prevent the % track immediately connecting back to the other junction. distj = norm(commonrc(n,:) - [RJ(j) CJ(j)]); distk = norm(commonrc(n,:) - [rj(k) cj(k)]); edgeNo = edgeNo + 1; if distj < distk edgelist{edgeNo} = trackedge(RJ(j), CJ(j), edgeNo, ... commonrc(n,1), commonrc(n,2), 1); else edgelist{edgeNo} = trackedge(rj(k), cj(k), edgeNo, ... commonrc(n,1), commonrc(n,2), 1); end etype(edgeNo) = 3; % Edge segment is junction-junction end end % Track any remaining unlabeled pixels adjacent to this junction k for m = 1:length(rak) if EDGEIM(rak(m), cak(m)) == 1 edgeNo = edgeNo + 1; edgelist{edgeNo} = trackedge(rj(k), cj(k), edgeNo, rak(m), cak(m)); etype(edgeNo) = 3; % Edge segment is junction-junction end end % Mark that we have visited junction (rj(k) cj(k)) JUNCT(rj(k), cj(k)) = 2; end % for all adjacent junctions % Finally track any remaining unlabeled pixels adjacent to original junction j for m = 1:length(ra) if EDGEIM(ra(m), ca(m)) == 1 edgeNo = edgeNo + 1; edgelist{edgeNo} = trackedge(RJ(j), CJ(j), edgeNo, ra(m), ca(m)); etype(edgeNo) = 3; % Edge segment is junction-junction end end end % If we have not visited this junction end % For each junction %% 3) Scan through the image looking for any unlabeled pixels. These % should correspond to isolated loops that have no junctions or endpoints. for ru = 1:ROWS for cu = 1:COLS if EDGEIM(ru,cu) == 1 % We have an unlabeled edge edgeNo = edgeNo + 1; [edgelist{edgeNo} endType] = trackedge(ru, cu, edgeNo); etype(edgeNo) = endType; end end end edgeim = -EDGEIM; % Finally negate image to make edge encodings +ve. % Eliminate isolated edges and spurs that are below the minimum length % ** DISABLED for the time being ** % if nargin >= 2 && ~isempty(minlength) % edgelist = cleanedgelist2(edgelist, minlength); % end % If subpixel edge locations are supplied upgrade the integer precision % edgelists that were constructed with data from 'location'. if nargin == 3 for I = 1:length(edgelist) ind = sub2ind(size(im),edgelist{I}(:,1),edgelist{I}(:,2)); edgelist{I}(:,1) = real(location(ind))'; edgelist{I}(:,2) = imag(location(ind))'; end end clear global EDGEIM; clear global ROWS; clear global COLS; clear global JUNCT; %---------------------------------------------------------------------- % TRACKEDGE % % Function to track all the edge points starting from an end point or junction. % As it tracks it stores the coords of the edge points in an array and labels the % pixels in the edge image with the -ve of their edge number. This continues % until no more connected points are found, or a junction point is encountered. % % Usage: edgepoints = trackedge(rstart, cstart, edgeNo, r2, c2, avoidJunction) % % Arguments: rstart, cstart - Row and column No of starting point. % edgeNo - The current edge number. % r2, c2 - Optional row and column coords of 2nd point. % avoidJunction - Optional flag indicating that (r2,c2) % should not be immediately connected to a % junction (if possible). % % Returns: edgepoints - Nx2 array of row and col values for % each edge point. % endType - 0 for a free end % 1 for a junction % 5 for a loop function [edgepoints endType] = trackedge(rstart, cstart, edgeNo, r2, c2, avoidJunction) global EDGEIM; global JUNCT; if ~exist('avoidJunction', 'var'), avoidJunction = 0; end edgepoints = [rstart cstart]; % Start a new list for this edge. EDGEIM(rstart,cstart) = -edgeNo; % Edge points in the image are % encoded by -ve of their edgeNo. preferredDirection = 0; % Flag indicating we have/not a % preferred direction. % If the second point has been supplied add it to the track and set the % path direction if exist('r2', 'var') && exist('c2', 'var') edgepoints = [edgepoints r2 c2 ]; EDGEIM(r2, c2) = -edgeNo; % Initialise direction vector of path and set the current point on % the path dirn = unitvector([r2-rstart c2-cstart]); r = r2; c = c2; preferredDirection = 1; else dirn = [0 0]; r = rstart; c = cstart; end % Find all the pixels we could link to [ra, ca, rj, cj] = availablepixels(r, c, edgeNo); while ~isempty(ra) || ~isempty(rj) % First see if we can link to a junction. Choose the junction that % results in a move that is as close as possible to dirn. If we have no % preferred direction, and there is a choice, link to the closest % junction % We enter this block: % IF there are junction points and we are not trying to avoid a junction % OR there are junction points and no non-junction points, ie we have % to enter it even if we are trying to avoid a junction if (~isempty(rj) && ~avoidJunction) || (~isempty(rj) && isempty(ra)) % If we have a prefered direction choose the junction that results % in a move that is as close as possible to dirn. if preferredDirection dotp = -inf; for n = 1:length(rj) dirna = unitvector([rj(n)-r cj(n)-c]); dp = dirn*dirna'; if dp > dotp dotp = dp; rbest = rj(n); cbest = cj(n); dirnbest = dirna; end end % Otherwise if we have no established direction, we should pick a % 4-connected junction if possible as it will be closest. This only % affects tracks of length 1 (Why do I worry about this...?!). else distbest = inf; for n = 1:length(rj) dist = sum(abs([rj(n)-r; cj(n)-c])); if dist < distbest rbest = rj(n); cbest = cj(n); distbest = dist; dirnbest = unitvector([rj(n)-r cj(n)-c]); end end preferredDirection = 1; end % If there were no junctions to link to choose the available % non-junction pixel that results in a move that is as close as possible % to dirn else dotp = -inf; for n = 1:length(ra) dirna = unitvector([ra(n)-r ca(n)-c]); dp = dirn*dirna'; if dp > dotp dotp = dp; rbest = ra(n); cbest = ca(n); dirnbest = dirna; end end avoidJunction = 0; % Clear the avoidJunction flag if it had been set end % Append the best pixel to the edgelist and update the direction and EDGEIM r = rbest; c = cbest; edgepoints = [edgepoints r c ]; dirn = dirnbest; EDGEIM(r, c) = -edgeNo; % If this point is a junction exit here if JUNCT(r, c); endType = 1; % Mark end as being a junction return; else % Get the next set of available pixels to link. [ra, ca, rj, cj] = availablepixels(r, c, edgeNo); end end % If we get here we are at an endpoint or our sequence of pixels form a % loop. If it is a loop the edgelist should have start and end points % matched to form a loop. If the number of points in the list is four or % more (the minimum number that could form a loop), and the endpoints are % within a pixel of each other, append a copy of the first point to the end % to complete the loop endType = 0; % Mark end as being free, unless it is reset below if length(edgepoints) >= 4 if abs(edgepoints(1,1) - edgepoints(end,1)) <= 1 && ... abs(edgepoints(1,2) - edgepoints(end,2)) <= 1 edgepoints = [edgepoints edgepoints(1,:)]; endType = 5; % Mark end as being a loop end end %---------------------------------------------------------------------- % AVAILABLEPIXELS % % Find all the pixels that could be linked to point r, c % % Arguments: rp, cp - Row, col coordinates of pixel of interest. % edgeNo - The edge number of the edge we are seeking to % track. If not supplied its value defaults to 0 % resulting in all adjacent junctions being returned, % (see note below) % % Returns: ra, ca - Row and column coordinates of available non-junction % pixels. % rj, cj - Row and column coordinates of available junction % pixels. % % A pixel is avalable for linking if it is: % 1) Adjacent, that is it is 8-connected. % 2) Its value is 1 indicating it has not already been assigned to an edge % 3) or it is a junction that has not been labeled -edgeNo indicating we have % not already assigned it to the current edge being tracked. If edgeNo is % 0 all adjacent junctions will be returned function [ra, ca, rj, cj] = availablepixels(rp, cp, edgeNo) global EDGEIM; global JUNCT; global ROWS; global COLS; % If edgeNo not supplied set to 0 to allow all adjacent junctions to be returned if ~exist('edgeNo', 'var'), edgeNo = 0; end ra = []; ca = []; rj = []; cj = []; % row and column offsets for the eight neighbours of a point roff = [-1 0 1 1 1 0 -1 -1]; coff = [-1 -1 -1 0 1 1 1 0]; r = rp+roff; c = cp+coff; % Find indices of arrays of r and c that are within the image bounds ind = find((r>=1 & r<=ROWS) & (c>=1 & c<=COLS)); % A pixel is avalable for linking if its value is 1 or it is a junction % that has not been labeled -edgeNo for i = ind if EDGEIM(r(i),c(i)) == 1 && ~JUNCT(r(i), c(i)); ra = [ra; r(i)]; ca = [ca; c(i)]; elseif (EDGEIM(r(i),c(i)) ~= -edgeNo) && JUNCT(r(i), c(i)); rj = [rj; r(i)]; cj = [cj; c(i)]; end end %--------------------------------------------------------------------- % UNITVECTOR Normalises a vector to unit magnitude % function nv = unitvector(v) nv = v./sqrt(v(:)'*v(:));