summaryrefslogtreecommitdiffstats
path: root/SD-VBS/common/toolbox/MultiNcut/quadedgep2.m
diff options
context:
space:
mode:
Diffstat (limited to 'SD-VBS/common/toolbox/MultiNcut/quadedgep2.m')
-rwxr-xr-xSD-VBS/common/toolbox/MultiNcut/quadedgep2.m188
1 files changed, 188 insertions, 0 deletions
diff --git a/SD-VBS/common/toolbox/MultiNcut/quadedgep2.m b/SD-VBS/common/toolbox/MultiNcut/quadedgep2.m
new file mode 100755
index 0000000..5041377
--- /dev/null
+++ b/SD-VBS/common/toolbox/MultiNcut/quadedgep2.m
@@ -0,0 +1,188 @@
1% function [xs,ys,gx,gy,par,threshold,mag,mage,g,FIe,FIo,mago] = quadedgep(I,par,threshold);
2% Input:
3% I = image
4% par = vector for 4 parameters
5% [number of filter orientations, number of scale, filter size, elongation]
6% To use default values, put 0.
7% threshold = threshold on edge strength
8% Output:
9% [x,y,gx,gy] = locations and gradients of an ordered list of edgels
10% x,y could be horizontal or vertical or 45 between pixel sites
11% but it is guaranteed that there [floor(y) + (floor(x)-1)*nr]
12% is ordered and unique. In other words, each edgel has a unique pixel id.
13% par = actual par used
14% threshold = actual threshold used
15% mag = edge magnitude
16% mage = phase map
17% g = gradient map at each pixel
18% [FIe,FIo] = odd and even filter outputs
19% mago = odd filter output of optimum orientation
20
21% Stella X. Yu, 2001
22
23% This is the multi scale version of the filtering
24% For the moment the parameters are defined by default. See line 34
25
26
27function [x,y,gx,gy,par,threshold,mag_s,mage,g,FIe,FIo,mago] = quadedgep2(I,par,data,threshold);
28
29
30if nargin<4 | isempty(threshold),
31 threshold = 0.1;
32end
33
34[r,c] = size(I);
35def_par = [4,30,3];
36
37display_on = 1;
38
39% take care of parameters, any missing value is substituted by a default value
40if nargin<2 | isempty(par),
41 par = def_par;
42end
43% par(end+1:4)=0;
44% par = par(:);
45% j = (par>0);
46% have_value = [ j, 1-j ];
47% j = 1; n_filter = have_value(j,:) * [par(j); def_par(j)];
48% j = 2; n_scale = have_value(j,:) * [par(j); def_par(j)];
49% j = 3; winsz = have_value(j,:) * [par(j); def_par(j)];
50% j = 4; enlong = have_value(j,:) * [par(j); def_par(j)];
51
52n_ori = par(1); %if it doesn't work, par<-def_par
53
54winsz = par(2);
55enlong = par(3);
56
57% always make filter size an odd number so that the results will not be skewed
58j = winsz/2;
59if not(j > fix(j) + 0.1),
60 winsz = winsz + 1;
61end
62
63% filter the image with quadrature filters
64if (isempty(data.W.scales))
65 error ('no scales entered');
66end
67
68n_scale=length(data.W.scales);
69filter_scales=data.W.scales;
70%
71% if strcmp(data.dataWpp.mode,'multiscale')
72% n_scale=size((data.dataWpp.scales),2);
73% filter_scales=data.dataWpp.scales;
74% else
75% filter_scales=data.dataWpp.scales(1);
76% n_scale=1;
77% end
78% if n_scale>0&&strcmp(data.dataWpp.mode,'multiscale')
79% if (~isempty(data.dataWpp.scales))
80% filter_scales=data.dataWpp.scales;
81% else
82% filter_scales=zeros(1,n_scale);
83%
84% for i=1:n_scale,
85% filter_scales(i)=(sqrt(2))^(i-1);
86% end
87% data.dataWpp.scales=filter_scales;
88% end
89% else filter_scale=1;
90% data.dataWpp.scales=filter_scales;
91% end
92%
93% %%%%%%% juste pour que ca tourne
94% if ~strcmp(data.dataWpp.mode,'multiscale')
95% filter_scales=data.dataWpp.scales(1);
96% n_scale=4;
97% end
98% %%%%%%%%%%%%
99
100FBo = make_filterbank_odd2(n_ori,filter_scales,winsz,enlong);
101FBe = make_filterbank_even2(n_ori,filter_scales,winsz,enlong);
102n = ceil(winsz/2);
103f = [fliplr(I(:,2:n+1)), I, fliplr(I(:,c-n:c-1))];
104f = [flipud(f(2:n+1,:)); f; flipud(f(r-n:r-1,:))];
105FIo = fft_filt_2(f,FBo,1);
106FIo = FIo(n+[1:r],n+[1:c],:);
107FIe = fft_filt_2(f,FBe,1);
108FIe = FIe(n+[1:r],n+[1:c],:);
109
110% compute the orientation energy and recover a smooth edge map
111% pick up the maximum energy across scale and orientation
112% even filter's output: as it is the second derivative, zero cross localize the edge
113% odd filter's output: orientation
114
115[nr,nc,nb] = size(FIe);
116
117FIe = reshape(FIe, nr,nc,n_ori,length(filter_scales));
118FIo = reshape(FIo, nr,nc,n_ori,length(filter_scales));
119
120mag_s = zeros(nr,nc,n_scale);
121mag_a = zeros(nr,nc,n_ori);
122
123mage = zeros(nr,nc,n_scale);
124mago = zeros(nr,nc,n_scale);
125mage = zeros(nr,nc,n_scale);
126mago = zeros(nr,nc,n_scale);
127
128
129
130for i = 1:n_scale,
131 mag_s(:,:,i) = sqrt(sum(FIo(:,:,:,i).^2,3)+sum(FIe(:,:,:,i).^2,3));
132 mag_a = sqrt(FIo(:,:,:,i).^2+FIe(:,:,:,i).^2);
133 [tmp,max_id] = max(mag_a,[],3);
134
135 base_size = nr * nc;
136 id = [1:base_size]';
137 mage(:,:,i) = reshape(FIe(id+(max_id(:)-1)*base_size+(i-1)*base_size*n_ori),[nr,nc]);
138 mago(:,:,i) = reshape(FIo(id+(max_id(:)-1)*base_size+(i-1)*base_size*n_ori),[nr,nc]);
139
140 mage(:,:,i) = (mage(:,:,i)>0) - (mage(:,:,i)<0);
141
142 if display_on,
143 ori_incr=pi/n_ori; % to convert jshi's coords to conventional image xy
144 ori_offset=ori_incr/2;
145 theta = ori_offset+([1:n_ori]-1)*ori_incr; % orientation detectors
146 % [gx,gy] are image gradient in image xy coords, winner take all
147
148 ori = theta(max_id);
149 ori = ori .* (mago(:,:,i)>0) + (ori + pi).*(mago(:,:,i)<0);
150 gy{i} = mag_s(:,:,i) .* cos(ori);
151 gx{i} = -mag_s(:,:,i) .* sin(ori);
152 g = cat(3,gx{i},gy{i});
153
154 % phase map: edges are where the phase changes
155 mag_th = max(max(mag_s(:,:,i))) * threshold;
156 eg = (mag_s(:,:,i)>mag_th);
157 h = eg & [(mage(2:r,:,i) ~= mage(1:r-1,:,i)); zeros(1,nc)];
158 v = eg & [(mage(:,2:c,i) ~= mage(:,1:c-1,i)), zeros(nr,1)];
159 [y{i},x{i}] = find(h | v);
160 k = y{i} + (x{i}-1) * nr;
161 h = h(k);
162 v = v(k);
163 y{i} = y{i} + h * 0.5; % i
164 x{i} = x{i} + v * 0.5; % j
165 t = h + v * nr;
166 gx{i} = g(k) + g(k+t);
167 k = k + (nr * nc);
168 gy{i} = g(k) + g(k+t);
169
170% figure(1);
171% clf;
172% imagesc(I);colormap(gray);
173% hold on;
174% quiver(x,y,gx,gy); hold off;
175% title(sprintf('scale = %d, press return',i));
176
177 % pause;
178 0;
179else
180 x = [];
181 y = [];
182 gx = [];
183 gy =[];
184 g= [];
185 end
186end
187
188