THESE/codes/graphe/Ncut_9/quadedgep.m

   1 % function [x,y,gx,gy,par,threshold,mag,mage,g,FIe,FIo,mago] = quadedgep(I,par,threshold);\r
   2 % Input:\r
   3 %    I = image\r
   4 %    par = vector for 4 parameters\r
   5 %      [number of filter orientations, number of scales, filter size, elongation]\r
   6 %      To use default values, put 0.\r
   7 %    threshold = threshold on edge strength\r
   8 % Output:\r
   9 %    [x,y,gx,gy] = locations and gradients of an ordered list of edgels\r
  10 %       x,y could be horizontal or vertical or 45 between pixel sites\r
  11 %       but it is guaranteed that there [floor(y) + (floor(x)-1)*nr] \r
  12 %       is ordered and unique.  In other words, each edgel has a unique pixel id.\r
  13 %    par = actual par used\r
  14 %    threshold = actual threshold used\r
  15 %    mag = edge magnitude\r
  16 %    mage = phase map\r
  17 %    g = gradient map at each pixel\r
  18 %    [FIe,FIo] = odd and even filter outputs\r
  19 %    mago = odd filter output of optimum orientation\r
  20 %\r
  21 % Stella X. Yu, 2001\r
  22 \r
  23 \r
  24 \r
  25 function [x,y,gx,gy,par,threshold,mag,mage,g,FIe,FIo,mago] = quadedgep(I,par,threshold);\r
  26 \r
  27 if nargin<3 | isempty(threshold),\r
  28     threshold = 0.2;\r
  29 end\r
  30 \r
  31 [r,c] = size(I);\r
  32 def_par = [8,1,20,3];\r
  33 \r
  34 % take care of parameters, any missing value is substituted by a default value\r
  35 if nargin<2 | isempty(par),\r
  36    par = def_par;\r
  37 end\r
  38 par(end+1:4)=0;\r
  39 par = par(:);\r
  40 j = (par>0);\r
  41 have_value = [ j, 1-j ];\r
  42 j = 1; n_filter = have_value(j,:) * [par(j); def_par(j)];\r
  43 j = 2; n_scale  = have_value(j,:) * [par(j); def_par(j)];\r
  44 j = 3; winsz    = have_value(j,:) * [par(j); def_par(j)];\r
  45 j = 4; enlong   = have_value(j,:) * [par(j); def_par(j)];\r
  46 \r
  47 % always make filter size an odd number so that the results will not be skewed\r
  48 j = winsz/2;\r
  49 if not(j > fix(j) + 0.1),\r
  50     winsz = winsz + 1;\r
  51 end\r
  52 \r
  53 % filter the image with quadrature filters\r
  54 FBo = make_filterbank_odd2(n_filter,n_scale,winsz,enlong);\r
  55 FBe = make_filterbank_even2(n_filter,n_scale,winsz,enlong);\r
  56 n = ceil(winsz/2);\r
  57 f = [fliplr(I(:,2:n+1)), I, fliplr(I(:,c-n:c-1))];\r
  58 f = [flipud(f(2:n+1,:)); f; flipud(f(r-n:r-1,:))];\r
  59 FIo = fft_filt_2(f,FBo,1); \r
  60 FIo = FIo(n+[1:r],n+[1:c],:);\r
  61 FIe = fft_filt_2(f,FBe,1);\r
  62 FIe = FIe(n+[1:r],n+[1:c],:);\r
  63 \r
  64 % compute the orientation energy and recover a smooth edge map\r
  65 % pick up the maximum energy across scale and orientation\r
  66 % even filter's output: as it is the second derivative, zero cross localize the edge\r
  67 % odd filter's output: orientation\r
  68 mag = sqrt(sum(FIo.^2,3)+sum(FIe.^2,3));\r
  69 mag_a = sqrt(FIo.^2+FIe.^2);\r
  70 [tmp,max_id] = max(mag_a,[],3);\r
  71 base_size = r * c;\r
  72 id = [1:base_size]';\r
  73 mage = reshape(FIe(id+(max_id(:)-1)*base_size),[r,c]);\r
  74 mage = (mage>0) - (mage<0);\r
  75 \r
  76 ori_incr=pi/n_filter; % to convert jshi's coords to conventional image xy\r
  77 ori_offset=ori_incr/2;\r
  78 theta = ori_offset+([1:n_filter]-1)*ori_incr; % orientation detectors\r
  79 % [gx,gy] are image gradient in image xy coords, winner take all\r
  80 mago = reshape(FIo(id+(max_id(:)-1)*base_size),[r,c]);\r
  81 ori = theta(max_id);\r
  82 ori = ori .* (mago>0) + (ori + pi).*(mago<0);\r
  83 gy = mag .* cos(ori);\r
  84 gx = -mag .* sin(ori);\r
  85 g = cat(3,gx,gy);\r
  86 \r
  87 % phase map: edges are where the phase changes\r
  88 mag_th = max(mag(:)) * threshold;\r
  89 eg = (mag>mag_th);\r
  90 h = eg & [(mage(2:r,:) ~= mage(1:r-1,:)); zeros(1,c)];\r
  91 v = eg & [(mage(:,2:c) ~= mage(:,1:c-1)), zeros(r,1)];\r
  92 [y,x] = find(h | v);\r
  93 k = y + (x-1) * r;\r
  94 h = h(k);\r
  95 v = v(k);\r
  96 y = y + h * 0.5; % i\r
  97 x = x + v * 0.5; % j\r
  98 t = h + v * r;\r
  99 gx = g(k) + g(k+t);\r
 100 k = k + (r * c);\r
 101 gy = g(k) + g(k+t);\r
 102 \r
 103 \r