% ss_edge_worm( filestem, plane_fileroot, k_min, k_max, k_inc, ks, spacing, filter_bank_filestem, rc_Crop, G1_init_thresh ) % ---------------------------------------------------------------------------------------------------------------------------- % % Component of Scale Space Edges, version 1.0. % % Find most scale space edge segments over a range of scales, with a particular search spacing. % Scale space edge following with heuristic annealing to improve edge segment continuity. % % WARNING: It is currently futile to try to read and understand this code. % % Input: % % filestem - Image file stem. % % plane_fileroot - File root of sample scale image planes (e.g. 'Clown_plane_k' % where the full name is like 'Clown_plane_k[#].mat') % % k_min - Minimum k (scale = 2^k) for edge search. Must be integral % and have a sample scale image plane, like 'Clown_plane_k[k_min].mat'. % % k_max - Maximum k (scale = 2^k) for edge search. % % k_inc - Scale resolution (increment). % % ks - Initial search scale power (scale = 2^ks) list. % % spacing - Search grid spacing list, at each initial search scale. % % filter_bank_filestem % - File stem of filter bank, like 'ss_DoG_filters'. % % rc_Crop - Cropping rectangle, relative to sample planes. % [ x_min x_max y_min y_max ] where [ -1 0 0 0] % indicates complete image. % % G1_init_thresh - Normalized gradient threshold for starting an edge segment. % In principle, the edge tracker finds all edges, regardless % of gradient magnitude, so this should be set to 0.0. In % practice, edges that have no points with a gradient above % ~.005 have no physical significance (they are due to noise). % Raising the threshold to this value will reduce the number of % low salience short segments significantly, and speed up the % the search a by a small factor (~10%). % % Output: % % A matlab file with a name like: % % [ filestem '_edges_info ].mat % % containing lists of all derived edge information (points, segments, % and nodes), a segment index map, a gradient map, and a scale map. % % % Mark Dow, last modified March 2002 % University of Oregon % Brain Development Lab / Lewis Center for NeuroImaging % dow@braindev.uoregon.edu % function ss_edge_worm( filestem, plane_fileroot, k_min, k_max, k_inc, ks, spacing, filter_bank_filestem, rc_Crop, G1_init_thresh ); clear global global k_mn global k_ic global k_mx global ss_gray global szFB % Size of filter bank. global szFHW % Filter half width. global im_offsets global filtsize global n_filters global scale global sigma global sz_FullImage global coords_ikn global k_kn global coords global hG0 global hG1x global hG1y global hG1_ikt global hG2_ikt global hG3_ikt global rcCrop global flim_visited global im_visited_seg global scale_jump_thresh global flim_scale global scale_intersect_thresh global seg_n % Hardcoded parameters. min_length = 1; % Don't record edges with length < min_length. scale_jump_thresh = 1.05; % Hysteresis factor for G1 ratio when jumping scale increments. % 1.0 for no hysteresis. G1'/G1 > 1/s_j_t. scale_intersect_thresh =1.2; % G3_thresh = -9999; % Third derivative in gradient direction threshold for valid edge. % -9999 for no filtering based on third derivative. nStopAtSegment = 99999; b_SplitSegments = 1; b_ExtendDanglingEnds = 0; k_mn = k_min; k_mx = k_max; k_ic = k_inc; % Load filter file. load( filter_bank_filestem ); % sigma % coords_ik % k_kn % hG0_ik % hG1_ikt, hG2_ikt, hG3_ikt szFB = size( hG2_ikt ); coords = zeros( 2*szFB(1), 2 ); % TO DO: make this k_max - k_min + 1 dimensional; a set of coords for each integral % scale. % coords_ik = coords_ik.*nk; szFHW = max(max(abs(coords_ik(:,1,1)))); % The filter width is constant for % all filters in this case. coords_ik = permute( coords_ik, [1 3 2] ); coords_ikn = zeros( 2*szFB(1), 2, k_max - k_min + 1 ); for n = 1 : k_max - k_min + 1 coords_ikn( 1:szFB(1), :, n ) = (2^(n-1))*coords_ik(:, :, 1); coords_ikn( szFB(1)+1:2*szFB(1), :, n ) = -(2^(n-1))*coords_ik(:, :, 1); end % Load scale space image: % For each integral scale: for k_n = k_min : k_max if k_n <= k_max % Load a single scale image. fullfilestem = [ plane_fileroot num2str(k_n) ]; load( fullfilestem ); % flim_gray % At first scale image: if k_n == k_min % Get size of full image. sz_FullImage = size( flim_gray ); msg = ['\n\n Starting edge search ...' ]; fprintf( msg ); msg = ['\n\n Size of full (resampled) image: ' num2str(sz_FullImage(2)) ' x ' num2str(sz_FullImage(1)) ]; fprintf( msg ); % Limit cropping rectangle to its intersection with the image. if rc_Crop(1) ~= -1 if rc_Crop(1) <= 0 rc_Crop(1) = 1; end if rc_Crop(3) <= 0 rc_Crop(3) = 1; end if rc_Crop(2) >= sz_FullImage(1) rc_Crop(2) = sz_FullImage(1); end if rc_Crop(4) >= sz_FullImage(2) rc_Crop(4) = sz_FullImage(2); end msg = ['\n\n Finding edges in cropping rectangle:' ]; fprintf( msg ); msg = ['\n x_min ' num2str(rc_Crop(3)) ]; fprintf( msg ); msg = ['\n x_max ' num2str(rc_Crop(4)) ]; fprintf( msg ); msg = ['\n y_min ' num2str(rc_Crop(1)) ]; fprintf( msg ); msg = ['\n y_max ' num2str(rc_Crop(2)) ]; fprintf( msg ); msg = ['\n Size of cropped image: ' num2str( rc_Crop(4) - rc_Crop(3) ) ' x ' num2str( rc_Crop(2) - rc_Crop(1) ) ]; fprintf( msg ); else rc_Crop(1) = 1; rc_Crop(2) = sz_FullImage(1); rc_Crop(3) = 1; rc_Crop(4) = sz_FullImage(2); msg = ['\n\n No cropping.' ]; fprintf( msg ); end sz_Crop = [ ( rc_Crop(2) - rc_Crop(1) + 1 ) ( rc_Crop(4) - rc_Crop(3) + 1 ) ]; flim_gray_cropped = flim_gray( rc_Crop(1):rc_Crop(2), rc_Crop(3):rc_Crop(4) ); % Calculate kMax such that any cropped image dimension is >= 2^(kMax+1): kMax = -1; sz_K = sz_Crop; while ( sz_K(1) > 1 & sz_K(2) > 1 ) sz_K = floor( sz_K./2 ); kMax = kMax + 1; end scale_max = 2^kMax; % If desired maximum scale larger than cropped image permits: if k_max > kMax % Warn user and change maximum scale. msg = ['\n\n Warning: Maximum scale requested is larger than the raw image permits. Maximum scale recalculated.' ]; fprintf( msg ); % Reset maximum scale. k_max = kMax k_mx = k_max; % TO DO: Make sure this new max works with the lower limit. end % Make sure initial search scale power, k_s is less than % maximum scale. bIsksReset = 0; for i = 1:length( ks ) if ks(i) >= k_max if bIsksReset == 0; ks(i) = k_max - .6; bIsksReset = 1; else % To Do: Should remove element from list. ks(i) = k_max; end end end % To Do: Warn user that initial scales are reset, indicate new initial scales. if bIsksReset == 1; msg = ['\n\n Warning: Initial scales reset.' ]; fprintf( msg ); ks end % Initialize scale space image. ss_gray = zeros( sz_FullImage(1), sz_FullImage(2), k_max - k_min + 1 ); end % At first scale image. % Load a single scale plane into scale space. ss_gray( 1:sz_FullImage(1), 1:sz_FullImage(2), k_n - k_min + 1 ) = flim_gray; end end % Each integral scale. clear flim_gray; rcCrop = rc_Crop; flim_visited = zeros(sz_Crop); im_visited_seg = zeros(sz_Crop); flim_gradient = zeros(sz_Crop); flim_scale = zeros(sz_Crop); flim_seg = zeros(sz_Crop); flim_edge_number = zeros(sz_Crop); % Edge point / segment / node tracking data structure. ep_max = floor( ( ( rcCrop(2) - rcCrop(1) )*( rcCrop(4) - rcCrop(3) ) )/5 ); ep_high = 0; seg_max = floor( ep_max/4 ); ep( 1:ep_max, 1:5 ) = 0; seg( 1:seg_max, 1:8 ) = 0; node( 1:2*seg_max, 1:7 ) = 0; ep_buf_max = ( ( rcCrop(2) - rcCrop(1) )+( rcCrop(4) - rcCrop(3) ) )*50; eps( 1:ep_buf_max, 1:5 ) = 0; epf( 1:ep_buf_max, 1:5 ) = 0; ns_termination = 0; % step_state at termination of edge segment %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % For a uniform coverage sampling of k-space, climb to gradient ridge and follow % edge in both directions % edge_list = zeros( [9999 5] ); nth_e = 0; nth_ep = 0; edge_number = -1; seg_n = 0; i_node = 0; eps_n = 0; epf_n = 0; tic nth_row = -1; n_rows = sum( floor( ( ( rcCrop(2)-rcCrop(1) )./spacing ) + 1 ) ); space_i = 0; b_intercept_s = 0; b_intercept_f = 0; for k_0 = ks k_init = k_0; space_i = space_i + 1; % For each nth pixel: for x = rcCrop(1): spacing(space_i) : rcCrop(2) nth_row = nth_row + 1; msg = ['\n\n Done with ' num2str(nth_row) ' of ' num2str(n_rows) ' rows.' ]; fprintf( msg ); msg = ['\n Searching along row ' num2str(x) ' at an initial k of ' num2str(k_0) ]; fprintf( msg ); for y = rc_Crop(3) : spacing(space_i) :rc_Crop(4) if seg_n < nStopAtSegment nth_ep = 0; % Follow gradient at one scale, looking for gradient ridge (G2 zero crossing). [ xr, yr, G1, G2, G3, tha, state ] = SS_GRidgeFindWorm( x, y, x, y, k_0, 0, 1 ); % If gradient ridge is found at original scale: if state == 2 % Look for better ridge at adjacent scales. [ xr, yr, kt, G1, G2, G3, tha, state ] = SS_KRidgeFindWorm( xr, yr, -1, -1, k_0, G1, G2, G3, tha, 0 ); % If best ridge has positive third derivative, % and scale space coordinate has not been marked: if state == 2 ... & G3 > G3_thresh/k_0 ... & G1 > G1_init_thresh/k_0 ... & ( flim_scale( xr - rcCrop(1) + 1, yr - rcCrop(3) + 1 ) == 0 ... | abs( flim_scale( xr - rcCrop(1) + 1, yr - rcCrop(3) + 1 ) - ( kt + 2 ) ) > scale_intersect_thresh ) nth_e = nth_e + 1; edge_number = nth_e; % New segment. seg_n = seg_n + 1; % Store first edge point. eps_n = eps_n + 1; eps( eps_n, 1 ) = xr - rcCrop(1) + 1; eps( eps_n, 2 ) = yr - rcCrop(3) + 1; eps( eps_n, 3 ) = G1; eps( eps_n, 4 ) = kt; eps( eps_n, 5 ) = tha; % Mark edge maps. nth_ep = nth_ep + 1; edge_list( nth_ep, 1 ) = xr - rcCrop(1) + 1; edge_list( nth_ep, 2 ) = yr - rcCrop(3) + 1; edge_list( nth_ep, 3 ) = G1; edge_list( nth_ep, 4 ) = kt + 2; edge_list( nth_ep, 5 ) = flim_scale( xr - rcCrop(1) + 1, yr - rcCrop(3) + 1 ); flim_scale( xr - rcCrop(1) + 1, yr - rcCrop(3) + 1 ) = kt + 2; flim_seg( xr - rcCrop(1) + 1, yr - rcCrop(3) + 1 ) = seg_n; flim_gradient( xr - rcCrop(1) + 1, yr - rcCrop(3) + 1 ) = max(G1, .001); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % Step in +edge direction until no edge found: % % xr_0 = xr; yr_0 = yr; G1_0 = G1; G2_0 = G2; G3_0 = G3; tha_0 = tha; kt_0 = kt; num_edge = 0; b_stepping = 1; % While stepping: while b_stepping == 1 num_edge = num_edge + 1; % Take one step along a ss edge in +edge direction. [ xr, yr, kt, G1, G2, G3, tha, step_state ] = SS_KRidgeStepWorm( xr, yr, kt, G1, G2, G3, tha, 1 ); % If still stepping, % and scale space coordinate has not been marked: if step_state == 2 ... & ( flim_scale( xr - rcCrop(1) + 1, yr - rcCrop(3) + 1 ) == 0 ... | abs( flim_scale( xr - rcCrop(1) + 1, yr - rcCrop(3) + 1 ) - ( kt + 2 ) ) > scale_intersect_thresh ) ... & G3 > G3_thresh/k_0 % Mark new edge point: % Find last step increment. xr0 = eps( eps_n, 1 ); yr0 = eps( eps_n, 2 ); xrp = xr - rcCrop(1) + 1; yrp = yr - rcCrop(3) + 1; dxr = xrp - xr0; dyr = yrp - yr0; % If not an immediate (8) neighbor; if abs(dxr) > 1 | abs(dyr) > 1 % Find coordinates of intermediate points: if abs(dxr) > 0 sign_x = dxr/abs(dxr); else sign_x = 0; end if abs(dyr) > 0 sign_y = dyr/abs(dyr); else sign_y = 0; end if abs(dxr) >= abs(dyr) nn = abs(dxr); rx = 1; ry = abs(dyr/dxr); else nn = abs(dyr); ry = 1; rx = abs(dxr/dyr); end % For each intermediate point: for n_s = 1:nn xrn = xr0 + sign_x*round(n_s*rx); yrn = yr0 + sign_y*round(n_s*ry); % Store point data. eps_n = eps_n + 1; eps( eps_n, 1 ) = xrn; eps( eps_n, 2 ) = yrn; eps( eps_n, 3 ) = G1; eps( eps_n, 4 ) = kt; eps( eps_n, 5 ) = tha; % Mark edge maps. nth_ep = nth_ep + 1; edge_list( nth_ep, 1 ) = xrn; edge_list( nth_ep, 2 ) = yrn; edge_list( nth_ep, 3 ) = G1; edge_list( nth_ep, 4 ) = kt + 2; edge_list( nth_ep, 5 ) = flim_scale( xrn, yrn ); if flim_scale( xrn, yrn ) == 0 flim_scale( xrn, yrn ) = kt + 2; flim_seg( xrn, yrn ) = seg_n; flim_gradient( xrn, yrn ) = max(G1, .001); end end % Else is an 8 neighbor: else % Store point data. eps_n = eps_n + 1; eps( eps_n, 1 ) = xr - rcCrop(1) + 1; eps( eps_n, 2 ) = yr - rcCrop(3) + 1; eps( eps_n, 3 ) = G1; eps( eps_n, 4 ) = kt; eps( eps_n, 5 ) = tha; % Mark edge maps. nth_ep = nth_ep + 1; edge_list( nth_ep, 1 ) = xr - rcCrop(1) + 1; edge_list( nth_ep, 2 ) = yr - rcCrop(3) + 1; edge_list( nth_ep, 3 ) = G1; edge_list( nth_ep, 4 ) = kt + 2; edge_list( nth_ep, 5 ) = flim_scale( xr - rcCrop(1) + 1, yr - rcCrop(3) + 1 ); if flim_scale( xr - rcCrop(1) + 1, yr - rcCrop(3) + 1 ) == 0 ... | abs( flim_scale( xr - rcCrop(1) + 1, yr - rcCrop(3) + 1 ) - ( kt + 2) ) > scale_intersect_thresh flim_scale( xr - rcCrop(1) + 1, yr - rcCrop(3) + 1 ) = kt + 2; flim_seg( xr - rcCrop(1) + 1, yr - rcCrop(3) + 1 ) = seg_n; flim_gradient( xr - rcCrop(1) + 1, yr - rcCrop(3) + 1 ) = max(G1, .001); end end % Else stop stepping, handle result in +edge direction: else b_stepping = 0; % Boundary in forward position: % step_state == -1 % Intersection with an existing edge: % step_state == 3 % Self intersection: % step_state == -4 % Already visited: % step_state == -2 % Couldn't find ridge: % step_state == 1 % Stepping, and G3 above threshold, but intersected a scale-space edge: % step_state == 2 % If intersection with an existing edge: if step_state == 3 | step_state == 2 b_intercept_s = 1; edge_number = nth_e; % Find last step increment. xr0 = eps( eps_n, 1 ); yr0 = eps( eps_n, 2 ); xrp = xr - rcCrop(1) + 1; yrp = yr - rcCrop(3) + 1; dxr = xrp - xr0; dyr = yrp - yr0; % If not an immediate (8) neighbor; if abs(dxr) > 1 | abs(dyr) > 1 % Find coordinates of intermediate points: if abs(dxr) > 0 sign_x = dxr/abs(dxr); else sign_x = 0; end if abs(dyr) > 0 sign_y = dyr/abs(dyr); else sign_y = 0; end if abs(dxr) >= abs(dyr) nn = abs(dxr); rx = 1; ry = abs(dyr/dxr); else nn = abs(dyr); ry = 1; rx = abs(dxr/dyr); end for n_s = 1:nn xrn = xr0 + sign_x*round(n_s*rx); yrn = yr0 + sign_y*round(n_s*ry); % Store point data. eps_n = eps_n + 1; eps( eps_n, 1 ) = xrn; eps( eps_n, 2 ) = yrn; eps( eps_n, 3 ) = G1; eps( eps_n, 4 ) = kt; eps( eps_n, 5 ) = tha; % Mark edge maps. nth_ep = nth_ep + 1; edge_list( nth_ep, 1 ) = xrn; edge_list( nth_ep, 2 ) = yrn; edge_list( nth_ep, 3 ) = G1; edge_list( nth_ep, 4 ) = kt + 2; edge_list( nth_ep, 5 ) = flim_scale( xrn, yrn ); if flim_scale( xrn, yrn ) == 0 flim_scale( xrn, yrn ) = kt + 2; flim_seg( xrn, yrn ) = seg_n; flim_gradient( xrn, yrn ) = max(G1, .001); end end % Else is an 8 neighbor: else % Mark edge point start, segment, and node. % Store point data. eps_n = eps_n + 1; eps( eps_n, 1 ) = xr - rcCrop(1) + 1; eps( eps_n, 2 ) = yr - rcCrop(3) + 1; eps( eps_n, 3 ) = G1; eps( eps_n, 4 ) = kt; eps( eps_n, 5 ) = tha; % Mark edge maps. nth_ep = nth_ep + 1; edge_list( nth_ep, 1 ) = xr; edge_list( nth_ep, 2 ) = yr; edge_list( nth_ep, 3 ) = G1; edge_list( nth_ep, 4 ) = kt + 2; edge_list( nth_ep, 5 ) = flim_scale( xr, yr ); if flim_scale( xr, yr ) == 0 flim_scale( xr, yr ) = kt + 2; flim_seg( xr, yr ) = seg_n; flim_gradient( xr, yr ) = max(G1, .001); end end % Code end node with segment that was intersected. step_state = flim_seg( xr - rcCrop(1) + 1, yr - rcCrop(3) + 1 ) + 10; end ns_termination = step_state; end % stop stepping, handle result in +edge direction end % stepping in +edge direction % % % Step in +edge direction until no edge found: % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % Step in -edge direction until no edge found: % % xr = xr_0; yr = yr_0; G1 = G1_0; G2 = G2_0; G3 = G3_0; tha = tha_0; kt = kt_0; num_edge = 0; b_stepping = 1; % While stepping: while b_stepping == 1 num_edge = num_edge + 1; % Take one step along a ss edge in -edge direction. [ xr, yr, kt, G1, G2, G3, tha, step_state ] ... = SS_KRidgeStepWorm( xr, yr, kt, G1, G2, G3, tha, -1 ); if seg_n == 42 iii =0; end % If best ridge has positive third derivative, % and scale space coordinate has not been marked: if step_state == 2 ... & ( flim_scale( xr - rcCrop(1) + 1, yr - rcCrop(3) + 1 ) == 0 ... | abs( flim_scale( xr - rcCrop(1) + 1, yr - rcCrop(3) + 1 ) - ( kt + 2 ) ) > scale_intersect_thresh ) ... & G3 > G3_thresh/k_0 if (epf_n == 0) xr0 = xr_0; yr0 = yr_0; else xr0 = epf( epf_n, 1 ); yr0 = epf( epf_n, 2 ); end xrp = xr - rcCrop(1) + 1; yrp = yr - rcCrop(3) + 1; dxr = xrp - xr0; dyr = yrp - yr0; if abs(dxr) > 1 | abs(dyr) > 1 if abs(dxr) > 0 sign_x = dxr/abs(dxr); else sign_x = 0; end if abs(dyr) > 0 sign_y = dyr/abs(dyr); else sign_y = 0; end if abs(dxr) >= abs(dyr) nn = abs(dxr); rx = 1; ry = abs(dyr/dxr); else nn = abs(dyr); ry = 1; rx = abs(dxr/dyr); end for n_s = 1:nn xrn = xr0 + sign_x*round(n_s*rx); yrn = yr0 + sign_y*round(n_s*ry); epf_n = epf_n + 1; epf( epf_n, 1 ) = xrn; epf( epf_n, 2 ) = yrn; epf( epf_n, 3 ) = G1; epf( epf_n, 4 ) = kt; epf( epf_n, 5 ) = tha; % Mark edge maps. nth_ep = nth_ep + 1; edge_list( nth_ep, 1 ) = xrn; edge_list( nth_ep, 2 ) = yrn; edge_list( nth_ep, 3 ) = G1; edge_list( nth_ep, 4 ) = kt + 2; edge_list( nth_ep, 5 ) = flim_scale( xrn, yrn ); if flim_scale( xrn, yrn ) == 0 flim_scale( xrn, yrn ) = kt + 2; flim_seg( xrn, yrn ) = seg_n; flim_gradient( xrn, yrn ) = max(G1, .001); end end else % Mark edge point start, segment, and node. epf_n = epf_n + 1; epf( epf_n, 1 ) = xr - rcCrop(1) + 1; epf( epf_n, 2 ) = yr - rcCrop(3) + 1; epf( epf_n, 3 ) = G1; epf( epf_n, 4 ) = kt; epf( epf_n, 5 ) = tha; % Mark edge maps. nth_ep = nth_ep + 1; edge_list( nth_ep, 1 ) = xr - rcCrop(1) + 1; edge_list( nth_ep, 2 ) = yr - rcCrop(3) + 1; edge_list( nth_ep, 3 ) = G1; edge_list( nth_ep, 4 ) = kt + 2; edge_list( nth_ep, 5 ) = flim_scale( xr - rcCrop(1) + 1, yr - rcCrop(3) + 1 ); if flim_scale( xr - rcCrop(1) + 1, yr - rcCrop(3) + 1 ) == 0 flim_scale( xr - rcCrop(1) + 1, yr - rcCrop(3) + 1 ) = kt + 2; flim_seg( xr - rcCrop(1) + 1, yr - rcCrop(3) + 1 ) = seg_n; flim_gradient( xr - rcCrop(1) + 1, yr - rcCrop(3) + 1 ) = max(G1, .001); end end % Else not a new edge point: else %step_state % Stop stepping. b_stepping = 0; % Boundary in forward position: % step_state == -1 % Intersection with an existing edge: % step_state == 3 % Self intersection: % step_state == -4 % Already visited: % step_state == -2 % Couldn't find ridge: % step_state == 1 % G3 threshold not reached: % step_state == 2 % If intersection with an existing edge: if step_state == 3 | step_state == 2 b_intercept_f = 1; edge_number = nth_e; if (epf_n == 0) xr0 = xr_0; yr0 = yr_0; else xr0 = epf( epf_n, 1 ); yr0 = epf( epf_n, 2 ); end xrp = xr - rcCrop(1) + 1; yrp = yr - rcCrop(3) + 1; dxr = xrp - xr0; dyr = yrp - yr0; if abs(dxr) > 1 | abs(dyr) > 1 if abs(dxr) sign_x = dxr/abs(dxr); else sign_x = 0; end if abs(dyr) sign_y = dyr/abs(dyr); else sign_y = 0; end if abs(dxr) >= abs(dyr) nn = abs(dxr); rx = 1; ry = abs(dyr/dxr); else nn = abs(dyr); ry = 1; rx = abs(dxr/dyr); end for n_s = 1:nn xrn = xr0 + sign_x*round(n_s*rx); yrn = yr0 + sign_y*round(n_s*ry); epf_n = epf_n + 1; epf( epf_n, 1 ) = xrn; epf( epf_n, 2 ) = yrn; epf( epf_n, 3 ) = G1; epf( epf_n, 4 ) = kt; epf( epf_n, 5 ) = tha; % Mark edge maps. nth_ep = nth_ep + 1; edge_list( nth_ep, 1 ) = xrn; edge_list( nth_ep, 2 ) = yrn; edge_list( nth_ep, 3 ) = G1; edge_list( nth_ep, 4 ) = kt + 2; edge_list( nth_ep, 5 ) = flim_scale( xrn, yrn ); if flim_scale( xrn, yrn ) == 0 flim_scale( xrn, yrn ) = kt + 2; flim_seg( xrn, yrn ) = seg_n; flim_gradient( xrn, yrn ) = max(G1, .001); end end else % Mark edge point start, segment, and node. epf_n = epf_n + 1; epf( epf_n, 1 ) = xr - rcCrop(1) + 1; epf( epf_n, 2 ) = yr - rcCrop(3) + 1; epf( epf_n, 3 ) = G1; epf( epf_n, 4 ) = kt; epf( epf_n, 5 ) = tha; % Mark edge maps. nth_ep = nth_ep + 1; edge_list( nth_ep, 1 ) = xr; edge_list( nth_ep, 2 ) = yr; edge_list( nth_ep, 3 ) = G1; edge_list( nth_ep, 4 ) = kt + 2; edge_list( nth_ep, 5 ) = flim_scale( xr, yr ); if flim_scale( xr, yr ) == 0 flim_scale( xr, yr ) = kt + 2; flim_seg( xr, yr ) = seg_n; flim_gradient( xr, yr ) = max(G1, .001); end end step_state = flim_seg( xr - rcCrop(1) + 1, yr - rcCrop(3) + 1 ) + 10; end nf_termination = step_state; end % handle result in -edge direction end % stepping in -edge direction % % % Step in -edge direction until no edge found: % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Find length/scale ratio. % if (eps_n + epf_n - 2) > min_length if (eps_n + epf_n - 2) > 0 lk_rat = (eps_n + epf_n - 2)/(( ( sum( 2.^epf( 1:epf_n-1, 4 ) ) + sum( 2.^eps( 1:eps_n-1, 4 ) ) )/(eps_n + epf_n - 2) )); else lk_rat = 0; end % If segment is connected and short (compared to scale), and doesn't connect % distant points on an edge: % To Do: Check for distant connection. if (ns_termination > 10 | ns_termination == -1)... & (nf_termination > 10 | nf_termination == -1)... & (lk_rat < 0) ... & (ns_termination == seg_n+10)... & (nf_termination == seg_n+10) % not self connected % connected, not boundary % Erase segment map between ends. for i = 1:epf_n-1 if flim_seg( epf( i, 1 ), epf( i, 2 ) ) == seg_n flim_seg( epf( i, 1 ), epf( i, 2 ) ) = 0; flim_scale( epf( i, 1 ), epf( i, 2 ) ) = 0; flim_gradient( epf( i, 1 ), epf( i, 2 ) ) = 0; end end for i = 1:eps_n-1 if flim_seg( eps( i, 1 ), eps( i, 2 ) ) == seg_n flim_seg( eps( i, 1 ), eps( i, 2 ) ) = 0; flim_scale( eps( i, 1 ), eps( i, 2 ) ) = 0; flim_gradient( eps( i, 1 ), eps( i, 2 ) ) = 0; end end % Backstep a segment. seg_n = seg_n - 1; % Else write edge point buffers into edge point array, mark segment, mark nodes: else % Write start direction end point. i_node = i_node + 1; seg( seg_n, 1) = ep_high + 1; seg( seg_n, 3) = i_node; node( i_node, 5 ) = node( i_node, 5 ) + 1;% + b_intercept_s; node( i_node, node(i_node,5) ) = seg_n; node( i_node, 6 ) = ns_termination; node( i_node, 7 ) = seg( seg_n, 1); % Write start edge points. first = ep_high + 1; last = ep_high + eps_n; ep( first:last, 1 ) = eps( eps_n:-1:1, 1 ); ep( first:last, 2 ) = eps( eps_n:-1:1, 2 ); ep( first:last, 3 ) = eps( eps_n:-1:1, 3 ); ep( first:last, 4 ) = eps( eps_n:-1:1, 4 ); ep( first:last, 5 ) = eps( eps_n:-1:1, 5 ); % Write finish edge points. first = ep_high + eps_n + 1; last = ep_high + eps_n + epf_n; ep( first:last, 1 ) = epf( 1:epf_n, 1 ); ep( first:last, 2 ) = epf( 1:epf_n, 2 ); ep( first:last, 3 ) = epf( 1:epf_n, 3 ); ep( first:last, 4 ) = epf( 1:epf_n, 4 ); ep( first:last, 5 ) = epf( 1:epf_n, 5 ); % Write finish direction end point. ep_high = ep_high + eps_n + epf_n; i_node = i_node + 1; seg( seg_n, 2) = ep_high; seg( seg_n, 4) = i_node; node( i_node, 5 ) = node( i_node, 5 ) + 1;% + b_intercept_f; node( i_node, node(i_node,5) ) = -1*seg_n; node( i_node, 6 ) = nf_termination; node( i_node, 7 ) = seg( seg_n, 2); % Write to screen edge segment found. % msg = ['\n Edge segment # ' num2str(seg_n) ', length: ' num2str( eps_n + epf_n ) ', ave. gradient: ' num2str( sum(edge_list(1:nth_ep,3))/( eps_n + epf_n ) ) ]; % fprintf( msg ); fprintf('\n Edge segment # '); fprintf( '%3d', seg_n ); fprintf(', length '); fprintf( '%4d', eps_n + epf_n ); fprintf(', gradient '); fprintf( '%0.3f', sum( edge_list(1:nth_ep,3) )/nth_ep ); fprintf(', scale '); fprintf( '%2.2f', 2^( -2 + sum( edge_list(1:nth_ep,4) )/nth_ep ) ); % Write edge to edge map. if nth_ep >= min_length for epn = 1:nth_ep flim_gradient( edge_list(epn,1), edge_list(epn,2) ) = max(edge_list(epn,3), .001); flim_scale( edge_list(epn,1), edge_list(epn,2) ) = edge_list(epn,4); flim_edge_number( edge_list(epn,1), edge_list(epn,2) ) = edge_number; end else for epn = 1:nth_ep flim_scale( edge_list(epn,1), edge_list(epn,2) ) = edge_list(epn,5); end end %%%%%%%%%%%%%%%%%%%% % Split segments at intersecting nodes: if b_SplitSegments == 1 seg_n_new = seg_n; i_SegIntersecting = seg_n; if seg_n == 39 ii = 0; end % If start node is intersecting an edge: if node( seg( i_SegIntersecting,3), 6 ) > 10 i_SegToSplit = node( seg( i_SegIntersecting,3), 6 ) - 10; i_epOfNode = seg( i_SegIntersecting,1); xn = ep( i_epOfNode, 1); yn = ep( i_epOfNode, 2); b_BeforeNewNode = 1; % If the node point is not an endpoint of be split, or it is self % intersecting not at its own endpoint: if ( ~( ep( seg( i_SegToSplit, 1 ), 1 ) == xn & ep( seg( i_SegToSplit, 1 ), 2 ) == yn ) ... & ~( ep( seg( i_SegToSplit, 2 ), 1 ) == xn & ep( seg( i_SegToSplit, 2 ), 2 ) == yn ) ) ... | ( ( i_SegToSplit == i_SegIntersecting ) ... & ~( ep( seg( i_SegToSplit, 1 ), 1 ) == ep( seg( i_SegToSplit, 2 ), 1 ) ) ... & ~( ep( seg( i_SegToSplit, 1 ), 2 ) == ep( seg( i_SegToSplit, 2 ), 2 ) ) ) % For each mid-point along intersected segment, until the new node point % is reached: for i_epOnSegToSplit = seg( i_SegToSplit, 1 ) + 1 : seg( i_SegToSplit, 2 ) - 1 % If this point is the new node: if ep( i_epOnSegToSplit, 1 ) == xn & ep( i_epOnSegToSplit, 2 ) == yn ... & b_BeforeNewNode == 1 seg_n_new = seg_n_new + 1; i_epOldLastPoint = seg( i_SegToSplit, 2 ); seg( i_SegToSplit, 2 ) = i_epOnSegToSplit; % set new last point of split segment i_nodeOldLastNode = seg( i_SegToSplit, 4 ); seg( i_SegToSplit, 4 ) = seg( i_SegIntersecting,3); % set new index to node of last point % Mark new segments at this node. node( seg( i_SegIntersecting, 3 ), node( seg( i_SegIntersecting, 3 ), 5 ) + 1 ) = i_SegToSplit; node( seg( i_SegIntersecting, 3 ), node( seg( i_SegIntersecting, 3 ), 5 ) + 2 ) = -1*seg_n_new; node( seg( i_SegIntersecting, 3 ), 5 ) = node( seg( i_SegIntersecting, 3 ), 5 ) + 2; node( seg( i_SegIntersecting, 3 ), 6 ) = 0; node( seg( i_SegIntersecting, 3 ), 7 ) = i_epOnSegToSplit; b_BeforeNewNode = 0; end end else % To Do: Else if a dangling end, combine % segments into one segment. seg_n end % Make a new segment for the points after the new node. if b_BeforeNewNode == 0 seg( seg_n_new, 1 ) = seg( i_SegToSplit, 2 ); seg( seg_n_new, 2 ) = i_epOldLastPoint; seg( seg_n_new, 3 ) = seg( i_SegIntersecting,3); seg( seg_n_new, 4 ) = i_nodeOldLastNode; if abs( node( i_nodeOldLastNode, 1 ) ) == i_SegToSplit node( i_nodeOldLastNode, 1 ) = seg_n_new; else if abs( node( i_nodeOldLastNode, 2 ) ) == i_SegToSplit node( i_nodeOldLastNode, 2 ) = seg_n_new; else if abs( node( i_nodeOldLastNode, 3 ) ) == i_SegToSplit node( i_nodeOldLastNode, 3 ) = seg_n_new; end end end % For each point after the new node point: for i_epNewSeg = seg( seg_n_new, 1 ) : seg( seg_n_new, 2 ); flim_seg( ep( i_epNewSeg, 1 ), ep( i_epNewSeg, 2 ) ) = seg_n_new; end end end % If end node is intersecting an edge: if node( seg( i_SegIntersecting,4), 6 ) > 10 i_SegToSplit = node( seg( i_SegIntersecting,4), 6 ) - 10; i_epOfNode = seg( i_SegIntersecting, 2 ); xn = ep( i_epOfNode, 1); yn = ep( i_epOfNode, 2); b_BeforeNewNode = 1; % If the node point is not an endpoint of be split, or it is self % intersecting not at its own endpoint: if ( ~( ep( seg( i_SegToSplit, 1 ), 1 ) == xn & ep( seg( i_SegToSplit, 1 ), 2 ) == yn ) ... & ~( ep( seg( i_SegToSplit, 2 ), 1 ) == xn & ep( seg( i_SegToSplit, 2 ), 2 ) == yn ) ) ... | ( ( i_SegToSplit == i_SegIntersecting ) ... & ~( ep( seg( i_SegToSplit, 1 ), 1 ) == ep( seg( i_SegToSplit, 2 ), 1 ) ) ... & ~( ep( seg( i_SegToSplit, 1 ), 2 ) == ep( seg( i_SegToSplit, 2 ), 2 ) ) ) % For each mid-point along intersected segment, until the new node point % is reached: for i_epOnSegToSplit = seg( i_SegToSplit, 1 ) + 1 : seg( i_SegToSplit, 2 ) - 1 % If this point is the new node: if ep( i_epOnSegToSplit, 1 ) == xn & ep( i_epOnSegToSplit, 2 ) == yn ... & b_BeforeNewNode == 1 seg_n_new = seg_n_new + 1; i_epOldLastPoint = seg( i_SegToSplit, 2 ); seg( i_SegToSplit, 2 ) = i_epOnSegToSplit; % set new last point of segment i_nodeOldLastNode = seg( i_SegToSplit, 4 ); seg( i_SegToSplit, 4 ) = seg( i_SegIntersecting,4); % set new index to node of last point % Mark new segments at this node. node( seg( i_SegIntersecting, 4 ), node( seg( i_SegIntersecting, 4 ), 5 ) + 1 ) = i_SegToSplit; node( seg( i_SegIntersecting, 4 ), node( seg( i_SegIntersecting, 4 ), 5 ) + 2 ) = -1*seg_n_new; node( seg( i_SegIntersecting, 4 ), 5 ) = node( seg( i_SegIntersecting, 4 ), 5 ) + 2; node( seg( i_SegIntersecting, 4 ), 6 ) = 0; node( seg( i_SegIntersecting, 4 ), 7 ) = i_epOnSegToSplit; b_BeforeNewNode = 0; end end else % To Do: Else if a dangling end, combine % segments into one segment. seg_n end % Make a new segment for the points after the new node. if b_BeforeNewNode == 0 seg( seg_n_new, 1 ) = seg( i_SegToSplit, 2 ); seg( seg_n_new, 2 ) = i_epOldLastPoint; seg( seg_n_new, 3 ) = seg( i_SegIntersecting, 4 ); seg( seg_n_new, 4 ) = i_nodeOldLastNode; if abs( node( i_nodeOldLastNode, 1 ) ) == i_SegToSplit node( i_nodeOldLastNode, 1 ) = seg_n_new; else if abs( node( i_nodeOldLastNode, 2 ) ) == i_SegToSplit node( i_nodeOldLastNode, 2 ) = seg_n_new; else if abs( node( i_nodeOldLastNode, 3 ) ) == i_SegToSplit node( i_nodeOldLastNode, 3 ) = seg_n_new; end end end % For each point after the new node point: for i_epNewSeg = seg( seg_n_new, 1 ) : seg( seg_n_new, 2 ); flim_seg( ep( i_epNewSeg, 1 ), ep( i_epNewSeg, 2 ) ) = seg_n_new; end end end seg_n = seg_n_new; end % Split segments at intersecting nodes: %%%%%%%%%%%%%%%%%%%% end % Reset. eps_n = 0; epf_n = 0; b_intercept_s = 0; b_intercept_f = 0; end % G3 > 0 (each new edge) end % Each new gradient ridge at original scale end % If seg_n ><= end % Each y end % Each x end toc % % Found all continuous segments % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % Extend dangling ends % if b_ExtendDanglingEnds % For each segment: for i_seg_n = 1:seg_n fprintf( ['\n Edge number: ' num2str(i_seg_n) ] ) % For each end of segment: for i_segment_end = 3:4 % 3 - start direction end (segment array index) % 4 - finish direction end b_StartExtension = 1; if i_segment_end == 3 direction = 1; else direction = -1; end % fprintf( ['\n Direction: ' num2str(direction) ] ) % If an end is dangling (node belongs to only one segment): if node( seg( i_seg_n, i_segment_end), 5 ) == 1 b_stepping = 1; % To Do: Get initial conditions again. xr = ep( seg( i_seg_n, i_segment_end-2 ), 1); yr = ep( seg( i_seg_n, i_segment_end-2 ), 2); G1 = ep( seg( i_seg_n, i_segment_end ), 3); k = ep( seg( i_seg_n, i_segment_end ), 4); G2 = -1; % not currently used as criteria for edge "extension" G3 = -1; tha = ep( seg( i_seg_n, i_segment_end ), 5); %%%%%%%%%%%%%%%%% % % While stepping: % while b_stepping == 1 state = 2; % alive, following edge extension % Find next best edge point: % Find 8-neighbor maximal direction. dx = -1*round( cos(tha+direction*pi/2) ); dy = -1*round( sin(tha+direction*pi/2) ); % % Collapse to 4-neighbor sub-optimal direction % if abs( cos(tha+pi/2) ) >= abs( sin(tha+pi/2) ) % dy = 0; % else % dx = 0; % end % For three forward neighbors, find best gradient: % Get unit step direction: dxc = dx; dyc = dy; if dxc == 0 step_tha = sign(dyc)*pi/2; else step_tha = atan(dyc/dxc); end if dxc < 0 step_tha = step_tha + pi; end step_tha = mod( step_tha + 2*pi, 2*pi ); % Get optimal gradient in center edge direction. [ xc, yc, kc, G1c, thac ] = SS_KGMax( xr + dxc, yr + dyc, k, step_tha, direction ); dG1c = 0; % If step and gradient are consistent with this edge: if ( abs( thac - mod( step_tha + pi/2 , 2*pi ) ) < pi/2 ) dG1c = (G1c - G1)/sqrt( dxc*dxc + dyc*dyc ); end % Check for boundary intersection. if xc < rcCrop(1) ... | xc > rcCrop(2) ... | yc < rcCrop(3) ... | yc > rcCrop(4) state = -1; end % Check for edge intersection. if state == 2 ... & flim_seg( xc - rcCrop(1) + 1, yc - rcCrop(3) + 1 ) ~= 0 ... & flim_seg( xc - rcCrop(1) + 1, yc - rcCrop(3) + 1 ) ~= i_seg_n state = 5; end % Get optimal gradient in counter clock-wise from center direction: if dx ==0 & dy == 1 dxa = -1; dya = 1; end if dx ==1 & dy == 1 dxa = 0; dya = 1; end if dx ==1 & dy == 0 dxa = 1; dya = 1; end if dx ==1 & dy == -1 dxa = 1; dya = 0; end if dx ==0 & dy == -1 dxa = 1; dya = -1; end if dx == -1 & dy == -1 dxa = 0; dya = -1; end if dx == -1 & dy == 0 dxa = -1; dya = -1; end if dx == -1 & dy == 1 dxa = -1; dya = 0; end if dxa == 0 step_tha = sign(dya)*pi/2; else step_tha = atan(dya/dxa); end if dxa < 0 step_tha = step_tha + pi; end step_tha = mod( step_tha + 2*pi, 2*pi ); % Get optimal gradient in ccw from center direction: [ xa, ya, ka, G1a, thaa ] = SS_KGMax( xr + dxa, yr + dya, k, step_tha, direction ); dG1a = 0; % If step and gradient are consistent with this edge: if ( abs( thaa - mod( step_tha + pi/2, 2*pi ) ) < pi/2 ) dG1a = (G1a - G1)/sqrt( dxa*dxa + dya*dya ); end % Check for boundary intersection. if xa < rcCrop(1) ... | xa > rcCrop(2) ... | ya < rcCrop(3) ... | ya > rcCrop(4) state = -1; end % Check for edge intersection. if state == 2 ... & flim_seg( xa - rcCrop(1) + 1, ya - rcCrop(3) + 1 ) ~= 0 ... & flim_seg( xa - rcCrop(1) + 1, ya - rcCrop(3) + 1 ) ~= i_seg_n state = 5; end % Get step clock-wise from center: if dx ==0 & dy == 1 dxb = 1; dyb = 1; end if dx ==1 & dy == 1 dxb = 1; dyb = 0; end if dx ==1 & dy == 0 dxb = 1; dyb = -1; end if dx ==1 & dy == -1 dxb = 0; dyb = -1; end if dx ==0 & dy == -1 dxb = -1; dyb = -1; end if dx == -1 & dy == -1 dxb = -1; dyb = 0; end if dx == -1 & dy == 0 dxb = -1; dyb = 1; end if dx == -1 & dy == 1 dxb = 0; dyb = 1; end if dxb == 0 step_tha = sign(dyb)*pi/2; else step_tha = atan(dyb/dxb); end if dxb < 0 step_tha = step_tha + pi; end step_tha = mod( step_tha + 2*pi, 2*pi ); % Get optimal gradient in cw direction. [ xb, yb, kb, G1b, thab ] = SS_KGMax( xr + dxb, yr + dyb, k, step_tha, direction ); dG1b = 0; % If step and gradient are consistent with this edge: if ( abs( thab - mod( step_tha + pi/2, 2*pi ) ) < pi/2 ) dG1b = (G1b - G1)/sqrt( dxb*dxb + dyb*dyb ); end % Check for boundary intersection. if xb < rcCrop(1) ... | xb > rcCrop(2) ... | yb < rcCrop(3) ... | yb > rcCrop(4) state = -1; end % Check for edge intersection. if state == 2 ... & flim_seg( xb - rcCrop(1) + 1, yb - rcCrop(3) + 1 ) ~= 0 ... & flim_seg( xb - rcCrop(1) + 1, yb - rcCrop(3) + 1 ) ~= i_seg_n state = 5; end % Pick best direction. if dG1c >= dG1a & dG1c >= dG1b xr = xc; yr = yc; kt = kc; G1 = G1c; tha = thac; end if dG1b >= dG1a & dG1b > dG1c xr = xb; yr = yb; kt = kb; G1 = G1b; tha = thab; end if dG1a > dG1c & dG1a > dG1b xr = xa; yr = ya; kt = ka; G1 = G1a; tha = thaa; end % Check for edge intersection. if state == 2 ... & flim_seg( xr - rcCrop(1) + 1, yr - rcCrop(3) + 1 ) ~= 0 state = 3; end % If still stepping: if ( state == 2 | state == 5 ) % Mark edge point start, segment, and node. ep_high = ep_high + 1; if b_StartExtension seg( i_seg_n, 3 + 2*(i_segment_end - 2) ) = ep_high; seg( i_seg_n, 4 + 2*(i_segment_end - 2) ) = ep_high; b_StartExtension = 0; else seg( i_seg_n, 4 + 2*(i_segment_end - 2) ) = ep_high; end ep( ep_high, 1 ) = xr - rcCrop(1) + 1; ep( ep_high, 2 ) = yr - rcCrop(3) + 1; ep( ep_high, 3 ) = G1; ep( ep_high, 4 ) = kt; ep( ep_high, 5 ) = tha; % Mark edge maps. nth_ep = nth_ep + 1; edge_list( nth_ep, 1 ) = xr - rcCrop(1) + 1; edge_list( nth_ep, 2 ) = yr - rcCrop(3) + 1; edge_list( nth_ep, 3 ) = G1; edge_list( nth_ep, 4 ) = kt + 2; edge_list( nth_ep, 5 ) = flim_scale( xr - rcCrop(1) + 1, yr - rcCrop(3) + 1 ); % If valid scale space coordinate has not been marked as an edge, % or self intersection not at gradient of last point: if ( flim_scale( xr - rcCrop(1) + 1, yr - rcCrop(3) + 1 ) == 0 ... | ((flim_seg( xr - rcCrop(1) + 1, yr - rcCrop(3) + 1 ) == i_seg_n) & (flim_gradient( xr - rcCrop(1) + 1, yr - rcCrop(3) + 1) ~= G1))) % Mark map flim_seg( xr - rcCrop(1) + 1, yr - rcCrop(3) + 1 ) = i_seg_n; flim_gradient( xr - rcCrop(1) + 1, yr - rcCrop(3) + 1 ) = G1; flim_scale( xr - rcCrop(1) + 1, yr - rcCrop(3) + 1 ) = kt + 2; % Else if intersecting an existing scale-space edge: else % if abs( flim_scale( xa - rcCrop(1) + 1, ya - rcCrop(3) + 1 ) - ( kt + 2) ) <= scale_intersect_thresh ... % | abs( flim_scale( xb - rcCrop(1) + 1, yb - rcCrop(3) + 1 ) - ( kt + 2) ) <= scale_intersect_thresh ... % | abs( flim_scale( xc - rcCrop(1) + 1, yc - rcCrop(3) + 1 ) - ( kt + 2) ) <= scale_intersect_thresh % Stop stepping. b_stepping = 0; % To Do: Mark new node. % end end if state == 5 % Stop stepping. b_stepping = 0; end % Else intersecting image boundary: else % Stop stepping. b_stepping = 0; % To Do: Mark new node. end end % while stepping % % Stepping. % %%%%%%%%%%%%%% end % each dangling segment end end % each segment end end % each segment end % if b_ExtendDanglingEnds % % % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % % % max_G = max( ep( :, 3 ) ); % fprintf('\n Maximum normalized gradient: '); % fprintf( '%1.2f', max_G ); % fprintf('\n'); % % seg_length(1:seg_n) = 0; % seg_ave_grad(1:seg_n) = 0; % seg_ave_scale(1:seg_n) = 0; % % % For each segment: % for i = 1:seg_n % % % Calculate salience measures. % seg_length(i) = seg(i, 2) - seg(i, 1) + 1; % seg_ave_grad(i) = sum( ep( seg(i, 1):seg(i, 2) , 3 ) )/seg_length(i); % seg_ave_scale(i) = sum( ep( seg(i, 1):seg(i, 2) , 4 ) )/seg_length(i); % % fprintf('\n Salience metrics of edge segment # '); % fprintf( '%3d', i ); % fprintf(', length '); % fprintf( '%3d', seg_length(i) ); % fprintf(' gradient '); % fprintf( '%1.3f', seg_ave_grad(i) ); % fprintf(' ave. scale '); % fprintf( '%1.3f', 2^seg_ave_scale(i) ); % end % % % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Trim segment array to number of segments. seg_temp = seg( 1:seg_n, : ); clear seg; seg = seg_temp; clear seg_temp; %save(filestem, 'ep', 'seg', 'node', 'ep_high', 'seg_n', 'i_node', 'flim_gray_cropped', 'flim_seg', 'flim_gradient', 'flim_scale', 'seg_length', 'seg_ave_grad', 'seg_ave_scale' ) save( [ filestem '_edge_info' ], 'ep', 'seg', 'node', 'ep_high', 'seg_n', 'i_node', 'flim_gray_cropped', 'flim_seg', 'flim_gradient', 'flim_scale' ) % figure % imagesc(flim_visited) % colormap(gray) % figure % imagesc(im_visited_seg) % colormap(cool) % figure % m = size(get(gcf,'colormap'),1); % rand('state',0); % randclr = rand(m,3); % randclr(1,:) = 0; % imagesc(flim_seg) % colormap(randclr) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % Step once along gradient ridge. % function [ xr, yr, kt, G1, G2, G3, tha, step_state ] = SS_KRidgeStepWorm( xr, yr, kt, G1, G2, G3, tha, direction ); global rcCrop global k_mn global k_mx global k_ic global flim_scale global scale_intersect_thresh global seg_n global im_visited_seg step_state = 2; state_n = -999; xr0 = xr; yr0 = yr; d0_max = sqrt( 2*max( 2^kt, 2 ) ); %xy0 = [xr0 yr0] xr_n = 0; yr_n = 0; G1_n = 0; G2_n = 0; G3_n = 0; tha_n = 0; % Get best 8-neighbor edge direction. dx = 0; dy = 0; if abs( G1*sin(tha) ) >= abs( G1*cos(tha) ) dx = direction*sign( sin(tha) ); if abs( cos(tha) ) > .5 dy = direction*-1*sign( cos(tha) ); end else dy = direction*-1*sign( cos(tha) ); if abs( sin(tha) ) > .5 dx = direction*sign( sin(tha) ); end end xr_n = xr + dx; yr_n = yr + dy; % If forward position is not in bounds: if xr_n < rcCrop(1) ... | xr_n > rcCrop(2) ... | yr_n < rcCrop(3) ... | yr_n > rcCrop(4) % Stop stepping. step_state = -1; return; end % If forward position intersects another edge: if flim_scale( xr_n - rcCrop(1) + 1, yr_n - rcCrop(3) + 1 ) ~= 0 ... & abs( flim_scale( xr_n - rcCrop(1) + 1, yr_n - rcCrop(3) + 1 ) - ( kt + 2) ) <= scale_intersect_thresh % Worm terminates with intersect edge code. % To Do: Get the new node point stats again? xr = xr_n; yr = yr_n; step_state = 3; return; end % Look for gradient ridge adjacent to input position: % To Do: Think about waggle k increment. waggle_inc_l = max( 1, round( min(0, kt - k_mn)/k_ic) ); waggle_inc_h = max( 1, round( min(0, k_mx - kt)/k_ic) ); % Waggle low. k_l = max(k_mn, kt - waggle_inc_l*k_ic); [ xr_l, yr_l, G1_l, G2_l, G3_l, tha_l, state_l ] = SS_GRidgeFindWorm( xr_n, yr_n, xr, yr, k_l, 1, 0 ); % [ xr_l, yr_l, G1_l, G2_l, G3_l, tha_l, state_l ] = SS_KRidgeFindWorm( xr_n, yr_n, xr, yr, k_l, G1, G2, G3, tha, 0 ); d0_l = sqrt( (xr_n - xr_l)^2 + (yr_n - yr_l)^2 ); d0_l2 = sqrt( (xr - xr_l)^2 + (yr - yr_l)^2 ); % Waggle high. k_h = min(k_mx - k_ic, kt + waggle_inc_h*k_ic); [ xr_h, yr_h, G1_h, G2_h, G3_h, tha_h, state_h ] = SS_GRidgeFindWorm( xr_n, yr_n, xr, yr, k_h, 1, 0 ); % [ xr_h, yr_h, G1_h, G2_h, G3_h, tha_h, state_h ] = SS_KRidgeFindWorm( xr_n, yr_n, xr, yr, k_h, G1, G2, G3, tha, 0 ); d0_h = sqrt( (xr_n - xr_h)^2 + (yr_n - yr_h)^2 ); d0_h2 = sqrt( (xr - xr_h)^2 + (yr - yr_h)^2 ); % Waggle middle. [ xr_m, yr_m, G1_m, G2_m, G3_m, tha_m, state_m ] = SS_GRidgeFindWorm( xr_n, yr_n, xr, yr, kt, 1, 0 ); % [ xr_m, yr_m, G1_m, G2_m, G3_m, tha_m, state_m ] = SS_KRidgeFindWorm( xr_n, yr_n, xr, yr, kt, G1, G2, G3, tha, 0 ); d0_m = sqrt( (xr_n - xr_m)^2 + (yr_n - yr_m)^2 ); d0_m2 = sqrt( (xr - xr_m)^2 + (yr - yr_m)^2 ); % If any new positions intersect previously mapped edge, not itself: if ( state_l == 3 & d0_l2 > 0 ) ... | ( state_h == 3 & d0_h2 > 0 ) ... | ( state_m == 3 & d0_m2 > 0 ) state_n = 3; if state_l == 3 & d0_l2 > 0 %& im_visited_seg( xr_l - rcCrop(1) + 1, yr_l - rcCrop(3) + 1 ) ~= seg_n xr = xr_l; yr = yr_l; G1 = G1_l; G2 = G2_l; G3 = G3_l; tha = tha_l; end if state_h == 3 & d0_h2 > 0 xr = xr_h; yr = yr_h; G1 = G1_h; G2 = G2_h; G3 = G3_h; tha = tha_h; end % if ( flim_scale( xr_m - rcCrop(1) + 1, yr_m - rcCrop(3) + 1 ) ~= 0 ... % & abs( flim_scale( xr_m - rcCrop(1) + 1, yr_m - rcCrop(3) + 1 ) - ( kt + 2) ) <= scale_intersect_thresh ) if state_m == 3 & d0_m2 > 0 %& im_visited_seg( xr_m - rcCrop(1) + 1, yr_m - rcCrop(3) + 1 ) ~= seg_n xr = xr_m; yr = yr_m; G1 = G1_m; G2 = G2_m; G3 = G3_m; tha = tha_m; end % Else not intersecting previous mapped edge: else state_n = 2; % default to if state_l ~= 2 | ( xr_l == xr0 & yr_l == yr0 ) | d0_l > d0_max | d0_l2 == 0 G1_l = 0; end if state_h ~= 2 | ( xr_h == xr0 & yr_h == yr0 ) | d0_h > d0_max | d0_h2 == 0 G1_h = 0; end if state_m ~= 2 | ( xr_m == xr0 & yr_m == yr0 ) | d0_m > d0_max | d0_m2 == 0 G1_m = 0; end % if ( state_l == 2 & G1_l ~= 0 ) | ( state_m == 2 & G1_m ~= 0 ) | ( state_h == 2 & G1_h ~= 0 ) % if ( state_l == 2 ) | ( state_m == 2 ) | ( state_h == 2 ) if G1_l > G1_m ... & G1_l > G1_h xr_n = xr_l; yr_n = yr_l; G1_n = G1_l; G2_n = G2_l; G3_n = G3_l; tha_n = tha_l; kt = k_l; % state_n = state_l; else if G1_m >= G1_h xr_n = xr_m; yr_n = yr_m; G1_n = G1_m; G2_n = G2_m; G3_n = G3_m; tha_n = tha_m; % state_n = state_m; else if G1_h > G1_m xr_n = xr_h; yr_n = yr_h; G1_n = G1_h; G2_n = G2_h; G3_n = G3_h; tha_n = tha_h; kt = k_h; % state_n = state_h; end end end end % no intersection with current edge % end if G1_n == 0 G1_n = G1; end % If new gradient ridge found: if state_n == 2 % Step to +/- edge position. xr = xr_n; yr = yr_n; G1 = G1_n; G2 = G2_n; G3 = G3_n; tha = tha_n; step_state = 2; % Else: else % Stop stepping. step_state = state_n; end % If stepping (successful at finding an adjacent ridge): if step_state == 2 % Look for better ridge at adjacent scales. % [ xr, yr, kt, G1, G2, G3, tha, state ] = SS_KRidgeFindWorm( xr, yr, xr0, yr0, kt, G1, G2, G3, tha, 0 ); end % if G1_l == 0 & G1_h == 0 & G1_m == 0 % step_state = -4; % end if ( xr == xr0 & yr == yr0 & step_state == 2) xr = xr + dx; yr = yr + dy; % [ xr, yr, kt, G1, G2, G3, tha, step_state ] = SS_KRidgeStepWorm( xr_n, yr_n, kt, G1, G2, G3, tha, direction ); end % If self intersection with same point (loop): if ( xr == xr0 & yr == yr0 & step_state == 2) step_state = -4; end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % Follow gradient ridge through scale, staying at G2 zero crossing. % function [ x, y, kt, G1, G2, G3, tha, state ] = SS_KRidgeFindWorm( x, y, xl, yl, k_0, G1, G2, G3, tha, b_StopIfVisited ); global k_mn global k_mx global k_ic global scale_jump_thresh G1_h = -999; kt = k_0; d0_skip_max = sqrt( 2*max( 2^k_0, 2 ) ); d0_max = sqrt( 1*max( 2^k_0, 2 ) ); % Follow ridge down scale, if rising, until valid ridge ends: kr = k_0; x0 = x; y0 = y; xr = x; yr = y; state = 2; state_h = 2; d0 = 1; % While scale larger than minimum: while state == 2 ... & kr >= k_mn + k_ic % If at integral scale: if abs( kr - 0 - floor(kr + .001) ) < .0001... & kr >= k_mn + k_ic % Check for gradient ridge over integral scale boundary: kr = kr - k_ic; [ xr, yr, G1_h, G2_h, G3_h, tha_h, state_h ] = SS_GRidgeFindWorm( xr, yr, x0, y0, kr, 1, b_StopIfVisited ); d0 = sqrt( (xr - x)^2 + (yr - y)^2 ); % Step to adjacent scale if it is a neighboring ridge. if state_h == 2 ... & d0 <= d0_skip_max ... & ( xr ~= xl & yr ~= yl ) x = xr; y = yr; G1 = G1_h; G2 = G2_h; G3 = G3_h; tha = tha_h; kt = kr; state = state_h; else if state_h == -2 & b_StopIfVisited state = state_h; end kr = k_mn; end end if state == 2 ... & d0 <= d0_skip_max ... & ( xr ~= xl & yr ~= yl )... & kr >= k_mn + k_ic % Test next scale lower: kr = kr - k_ic; [ xr, yr, G1_h, G2_h, G3_h, tha_h, state_h ] = SS_GRidgeFindWorm( xr, yr, x0, y0, kr, 1, b_StopIfVisited ); d0 = sqrt( (xr - x)^2 + (yr - y)^2 ); end if state_h == 2 ... & G1_h >= (1/scale_jump_thresh)*G1 ... & d0 <= d0_max ... & ( xr ~= xl & yr ~= yl )... & kr >= k_mn + k_ic x = xr; y = yr; G1 = G1_h; G2 = G2_h; G3 = G3_h; tha = tha_h; kt = kr; state = state_h; else if state_h == -2 & b_StopIfVisited state = state_h; end kr = k_mn; end end % while kr >= k_mn + k_ic % If didn't step to a lower scale: if kt == k_0 % Follow ridge up scale, if rising, until valid ridge ends: kr = k_0; xr = x; yr = y; state = 2; while kr <= k_mx - 2*k_ic & state == 2; % If just below an integral scale: if abs( kr - (1-k_ic) - floor(kr) ) < .0001 % Jump over integral scale boundary: kr = kr + k_ic; [ xr, yr, G1_h, G2_h, G3_h, tha_h, state_h ] = SS_GRidgeFindWorm( xr, yr, x0, y0, kr, 1, b_StopIfVisited ); d0 = sqrt( (xr - x)^2 + (yr - y)^2 ); % Step to adjacent scale if still a ridge. if state_h == 2 ... & d0 <= d0_skip_max ... & ( xr ~= xl & yr ~= yl ) x = xr; y = yr; G1 = G1_h; G2 = G2_h; G3 = G3_h; tha = tha_h; kt = kr; state = state_h; else if state_h == -2 & b_StopIfVisited state = state_h; end kr = k_mx; % exit loop end end if state == 2 ... & d0 <= d0_skip_max ... & ( xr ~= xl & yr ~= yl ) ... & kr <= k_mn - 2*k_ic % Test next scale higher: kr = kr + k_ic; [ xr, yr, G1_h, G2_h, G3_h, tha_h, state_h ] = SS_GRidgeFindWorm( xr, yr, x0, y0, kr, 1, b_StopIfVisited ); d0 = sqrt( (xr - x)^2 + (yr - y)^2 ); end if state == 2 ... & G1_h >= (1/scale_jump_thresh)*G1 ... & d0 <= d0_max ... & ( xr ~= xl & yr ~= yl ) x = xr; y = yr; G1 = G1_h; G2 = G2_h; G3 = G3_h; tha = tha_h; kt = kr; state = state_h; else if state_h == -2 & b_StopIfVisited state = state_h; end kr = k_mx; % exit loop end end % while kr <= k_mx - 2*k_ic end % If didn't step to a lower scale: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % Follow gradient at one scale, looking for G2 zero crossing. % function [ xr, yr, G1, G2, G3, tha, state ] = SS_GRidgeFindWorm( x, y, x0, y0, k, b_local, b_stop_if_visited ); global ss_gray global szFB global szFHW % Filter half width. global sz_FullImage % Size of full image. global coords_ikn global coords global k_kn global k_mn global im_offsets % List of indices into image, pointing to convolution coordinates. global hG1_ikt global hG2_ikt global hG3_ikt global rcCrop global flim_visited global im_visited_seg global flim_scale global scale_intersect_thresh global seg_n n_steps_max = 99999; if b_local == 1 n_steps_max = max( 2*2^k, 3 ); end nk = 2^floor(k + .0001); fk = k - floor(k + .0001); k_offsets = (floor(k + .0001) - k_mn)*sz_FullImage(1)*sz_FullImage(2); % If ss position has been visited: if b_stop_if_visited ... & abs( flim_visited( x - rcCrop(1) + 1, y - rcCrop(3) + 1 ) - ( k + 2 ) ) <= scale_intersect_thresh ... & im_visited_seg( x - rcCrop(1) + 1, y - rcCrop(3) + 1 ) ~= seg_n % Worm terminates with visited code. xr = x; yr = y; G0 = -1; G1 = -1; G2 = -1; G3 = -1; tha = -1; state = -2; return; end % Initialize worm's state: xr = x; yr = y; n_steps = 0; G3 = -1; state = 2; % Alive. % Calculate gradient components and angle: % Get index into list of fractional k filters available. i_fk = find( round(1000*k_kn) == round(1000*fk) ); if isempty(i_fk) 'Error in SS_G1_RidgeWorm: looking for filter not in this filter bank at ' k_scale = k valid_filter_scale_values = k_kn' return end % Find derivatives at this position: % Get coordinates for filter at this scale. coords = coords_ikn( :, :, floor(k) - k_mn + 1 ); GetImageOffsets( xr, yr ); % Calculate gradient and angle. hG1x(1:szFB(1)) = hG1_ikt( :, i_fk, 1 ); hG1x(szFB(1)+1:2*szFB(1)) = -hG1_ikt( :, i_fk, 1 ); G1x = sum( hG1x'.*ss_gray( k_offsets + im_offsets ) ); hG1y(1:szFB(1)) = hG1_ikt( :, i_fk, 2 ); hG1y(szFB(1)+1:2*szFB(1)) = -hG1_ikt( :, i_fk, 2 ); G1y = sum( hG1y'.*ss_gray( k_offsets + im_offsets ) ); G1 = sqrt( G1x*G1x + G1y*G1y ); if G1 == 0 xr = x; yr = y; G0 = -1; G1 = -1; G2 = -1; G3 = -1; tha = -1; state = 1; return; end if G1x == 0 tha = pi/2; else tha = atan(G1y/G1x); end if G1x < 0 tha = tha + pi; end tha = mod( tha + 2*pi, 2*pi ); n_tha = mod( round( tha*(2*szFB(3)/(2*pi) ) ), 2*szFB(3) ); % Calculate second derivative. hG2(1:szFB(1)) = hG2_ikt( :, i_fk, mod( n_tha, szFB(3) ) + 1 ); hG2(szFB(1)+1:2*szFB(1)) = hG2_ikt( :, i_fk, mod( n_tha, szFB(3) ) + 1 ); G2 = -sum( hG2'.*ss_gray( k_offsets + im_offsets ) ); % Set second derivative sign indicator. G2_s = sign(G2); G1x_s = sign(cos(tha)); G1y_s = sign(sin(tha)); % Mark this position as visited. flim_visited( xr - rcCrop(1) + 1, yr - rcCrop(3) + 1 ) = k + 2; im_visited_seg( xr - rcCrop(1) + 1, yr - rcCrop(3) + 1 ) = seg_n; % Find 8-neighbor maximal direction. fl_dx = G2_s*G1x/G1; fl_dy = G2_s*G1y/G1; dx = round(fl_dx); dy = round(fl_dy); % Collapse to 4-neighbor sub-optimal direction if abs(G1x) >= abs(G1y) dy = 0; else dx = 0; end % Find 4-neighbor maximal direction. % dx = 0; % dy = 0; % if abs(G1x) >= abs(G1y) % % dx = G2_s*sign(G1x); % else % % dy = G2_s*sign(G1y); % end % % Get best 8-neighbor edge direction. % dx = 0; % dy = 0; % if abs( G1*cos(tha) ) >= abs( G1*sin(tha) ) % dx = G2_s*sign( cos(tha) ); % if abs( sin(tha) ) > .5 % dy = G2_s*sign( sin(tha) ); % end % else % dy = G2_s*sign( sin(tha) ); % if abs( cos(tha) ) > .5 % dx = G2_s*sign( cos(tha) ); % end % end % While worm is alive: while state == 2 & n_steps < n_steps_max % Check for valid forward position: % If forward position is not in bounds: if x + dx < rcCrop(1) ... | x + dx > rcCrop(2) ... | y + dy < rcCrop(3) ... | y + dy > rcCrop(4) % Worm terminates with boundary code. state = -1; return; end % If forward position has been visited: if b_stop_if_visited ... & abs( flim_visited( x + dx - rcCrop(1) + 1, y + dy - rcCrop(3) + 1 ) - ( k + 2 ) ) <= scale_intersect_thresh %... % & im_visited_seg( x + dx - rcCrop(1) + 1, y + dy - rcCrop(3) + 1 ) ~= seg_n % Worm terminates with visited code. state = -2; return; end if flim_scale( x + dx - rcCrop(1) + 1, y + dy - rcCrop(3) + 1 ) ~= 0 ... & abs( flim_scale( x + dx - rcCrop(1) + 1, y + dy - rcCrop(3) + 1 ) - ( k + 2) ) <= scale_intersect_thresh % Worm terminates with intersect edge code. state = 3; % return; end % Climb ss gradient: % Find derivatives in forward direction: GetImageOffsets( x + dx, y + dy ); % Calculate gradient and angle: G1x_p = sum( hG1x'.*ss_gray( k_offsets + im_offsets ) ); G1y_p = sum( hG1y'.*ss_gray( k_offsets + im_offsets ) ); G1_p = sqrt( G1x_p*G1x_p + G1y_p*G1y_p ); if G1x_p == 0 tha_p = pi/2; else tha_p = atan(G1y_p/G1x_p); end if G1x_p < 0 tha_p = tha_p + pi; end tha_p = mod( tha_p + 2*pi, 2*pi ); n_tha_p = mod( round( tha_p*(2*szFB(3)/(2*pi) ) ), 2*szFB(3) ); % Calculate second derivative. hG2(1:szFB(1)) = hG2_ikt( :, i_fk, mod( n_tha_p, szFB(3) ) + 1 ); hG2(szFB(1)+1:2*szFB(1)) = hG2_ikt( :, i_fk, mod( n_tha_p, szFB(3) ) + 1 ); G2_p = -sum( hG2'.*ss_gray( k_offsets + im_offsets ) ); % Mark this position as visited. flim_visited( x + dx - rcCrop(1) + 1, y + dy - rcCrop(3) + 1 ) = k + 2; im_visited_seg( x + dx - rcCrop(1) + 1, y + dy - rcCrop(3) + 1 ) = seg_n; % If forward neighbor has same second derivative sign, take a step: if G2_p*G2 > 0 % Assign the new state. n_steps = n_steps + 1; xr = x + dx; yr = y + dy; G1 = G1_p; G1x = G1x_p; G1y = G1y_p; G2 = G2_p; tha = tha_p; % Find 8-neighbor maximal direction. fl_dxn = G2_s*G1x/G1; fl_dyn = G2_s*G1y/G1; if x + round(dx + fl_dxn) == x + dx && y + round(dy + fl_dyn) == y + dy dx = round(dx + 1.42*fl_dxn); dy = round(dy + 1.42*fl_dyn); else dx = round(dx + fl_dxn); dy = round(dy + fl_dyn); end % Collapse to 4-neighbor direction if abs(G1x) >= abs(G1y) dy = 0; else dx = 0; end % Else worm is at second derivative zero crossing. Find optimal adjacent % pixel and terminate worm. else %G3 % Calculate third derivative of this position. hG3(1:szFB(1)) = sign( szFB(3) - n_tha )*hG3_ikt( :, i_fk, mod( n_tha, szFB(3) ) + 1 ); hG3(szFB(1)+1:2*szFB(1)) = -sign( szFB(3) - n_tha )*hG3_ikt( :, i_fk, mod( n_tha, szFB(3) ) + 1 ); G3 = -sum( hG3'.*ss_gray( k_offsets + im_offsets ) ); %G3 % Calculate third derivative of forward position. if n_tha_p ~= n_tha hG3(1:szFB(1)) = sign( szFB(3) - n_tha_p )*hG3_ikt( :, i_fk, mod( n_tha_p, szFB(3) ) + 1 ); hG3(szFB(1)+1:2*szFB(1)) = -sign( szFB(3) - n_tha_p )*hG3_ikt( :, i_fk, mod( n_tha_p, szFB(3) ) + 1 ); end G3_p = -sum( hG3'.*ss_gray( k_offsets + im_offsets ) ); %G3 if abs( G2 - G2_p ) == 0 frac = 1; else frac = abs( G2 )/abs( G2 - G2_p ); end %G3 % Interpolate third derivative at this position G3 = G3 + frac*( G3_p - G3 ); % If forward position's second derivative is closer to zero, take a step: if abs(G2_p) < abs(G2) % Assign the new state. xr = x + dx; yr = y + dy; G1 = G1_p; G1x = G1x_p; G1y = G1y_p; G2 = G2_p; tha = tha_p; end % Normalize first and third derivative. f = 2^(-.5*k); G3 = G3*f; G1 = G1*f; % Worm terminates with success code. state = 2; return; end end % while alive %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % Find neighboring scale-space point of maximum gradient for a single space % position. [ To Do: Valid points must have gradient angle consistent with previous gradient % and step direction to the position.] % function [ x, y, kt, G1, tha ] = SS_KGMax( xr, yr, k, step_tha, direction ); global k_mn global k_mx global k_ic % Scale increment global im_offsets % List of indices into image, pointing to convolution coordinates. global coords_ikn global coords global scale_jump_thresh G1 = 0; % Get coordinates for filter at this scale. coords = coords_ikn( :, :, floor(k) - k_mn + 1 ); % Initialize im_offsets list. GetImageOffsets( xr, yr ); % Get gradient at initial scale: [ G1, tha ] = SS_G1( xr, yr, k ); kt = k; % Get gradient at incrementally lower scale. dk = -1*k_ic; G1_mk = 0; if floor(k + dk) - k_mn >= 0 % Get coordinates for filter at this scale. coords = coords_ikn( :, :, floor(k + dk) - k_mn + 1 ); % Initialize im_offsets list. GetImageOffsets( xr, yr ); [ G1_mk, tha_mk ] = SS_G1( xr, yr, max( k + dk, k_mn+k_ic ) ); % If higher gradient at lower scale, climb down scale to peak gradient. while G1_mk > (1/scale_jump_thresh)*G1 & dk + k_ic > k_mn + k_ic kt = k + dk; G1 = G1_mk; tha = tha_mk; dk = dk - k_ic; if floor(k + dk) - k_mn >= 0 % Get coordinates for filter at this scale. coords = coords_ikn( :, :, floor(k + dk) - k_mn + 1 ); % Initialize im_offsets list. GetImageOffsets( xr, yr ); [ G1_mk, tha_mk ] = SS_G1( xr, yr, max( k + dk, k_mn+k_ic ) ); end end end % Get gradient at incrementally higher scale: dk = k_ic; % Get coordinates for filter at this scale. coords = coords_ikn( :, :, floor(k + dk) - k_mn + 1 ); % Initialize im_offsets list. GetImageOffsets( xr, yr ); [ G1_pk, tha_pk ] = SS_G1( xr, yr, min( k + dk, k_mx-1 ) ); % If higher gradient at higher scale, climb up scale to peak gradient. while G1_pk > (1/scale_jump_thresh)*G1 & dk + k_ic < k_mx - 1 kt = k + dk; G1 = G1_pk; tha = tha_pk; dk = dk + k_ic; [ G1_pk, tha_pk ] = SS_G1( xr, yr, min( k + dk, k_mx-1 ) ); end x = xr; y = yr; return; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % Find gradient at a single scale-space position. % function [ G1, tha ] = SS_G1( xr, yr, k ); global ss_gray global szFB global szFHW % Filter half width. global sz_FullImage % Size of full image. global k_kn global k_mn global im_offsets % List of indices into image, pointing to convolution coordinates. global hG1_ikt global rcCrop fk = k - floor(k + .0001); % fractional part of scale k_offsets = (floor(k + .0001) - k_mn)*sz_FullImage(1)*sz_FullImage(2); % Full image offset for every integral scale above minimum. % Calculate gradient components and angle: % Get index into list of fractional k filters available. i_fk = find( round(1000*k_kn) == round(1000*fk) ); if isempty(i_fk) 'Error in SS_KGMax: looking for filter not in this filter bank at ' k_scale = k valid_filter_scale_values = k_kn' return end % Find derivative at this position: % Calculate gradient and angle. hG1x(1:szFB(1)) = hG1_ikt( :, i_fk, 1 ); hG1x(szFB(1)+1:2*szFB(1)) = -hG1_ikt( :, i_fk, 1 ); G1x = sum( hG1x'.*ss_gray( k_offsets + im_offsets ) ); hG1y(1:szFB(1)) = hG1_ikt( :, i_fk, 2 ); hG1y(szFB(1)+1:2*szFB(1)) = -hG1_ikt( :, i_fk, 2 ); G1y = sum( hG1y'.*ss_gray( k_offsets + im_offsets ) ); G1 = sqrt( G1x*G1x + G1y*G1y ); if G1x == 0 tha = pi/2; else tha = atan(G1y/G1x); end if G1x < 0 tha = tha + pi; end tha = mod( tha + 2*pi, 2*pi ); % Normalize first derivative. f = 2^(-.5*k); G1 = G1*f; return; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % GetImageOffsets( x, y ) % ----------------------- % % Sets a 1-D list of indices into an image, corresponding to a set of % filter coordinates. Folds off image coordinates back onto image. % % Mark Dow, September 2001 function GetImageOffsets( x, y ); global szFB % Size of filter bank. global szFHW % Filter half width. global sz_FullImage % Size of full image. global coords % n x 2 list of zero based filter coordinates. global im_offsets % List of indices into image, pointing to convolution coordinates. c = coords; c(:,1) = c(:,1) + x; c(:,2) = c(:,2) + y; % If filter overlaps an edge of the image, % fold off image coordinates back onto image. if x - szFHW < 1 i_x_low = find( c(:,1) < 1 ); c( i_x_low ) = -c( i_x_low ) + 1; end if x + szFHW > sz_FullImage(1) i_x_high = find( c(:,1) > sz_FullImage(1) ); c( i_x_high ) = -c( i_x_high ) + 2*sz_FullImage(1); end if y - szFHW < 1 i_y_low = find( c(:,2) < 1 ); c( i_y_low + 2*szFB(1) ) = -c( i_y_low + 2*szFB(1) ) + 1; end if y + szFHW > sz_FullImage(2) i_y_high = find( c(:,2) > sz_FullImage(2) ); c( i_y_high + 2*szFB(1) ) = -c( i_y_high + 2*szFB(1) ) + 2*sz_FullImage(2); end im_offsets = c(:,1) + (c(:,2)-1)*sz_FullImage(1);