1 // This is core/vil/algo/vil_sobel_3x3.cxx
2 #include "vil_sobel_3x3.h"
3 //:
4 // \file
5 // \brief Apply gradient operator to 2D planes of data
6 // \author Tim Cootes
7
8 //: Compute gradients of single plane of 2D data using 3x3 Sobel filters
9 // Computes both i and j gradients of an ni x nj plane of data
10 template <>
11 void
vil_sobel_3x3_1plane(const unsigned char * src,std::ptrdiff_t s_istep,std::ptrdiff_t s_jstep,float * gi,std::ptrdiff_t gi_istep,std::ptrdiff_t gi_jstep,float * gj,std::ptrdiff_t gj_istep,std::ptrdiff_t gj_jstep,unsigned ni,unsigned nj)12 vil_sobel_3x3_1plane(const unsigned char * src,
13 std::ptrdiff_t s_istep,
14 std::ptrdiff_t s_jstep,
15 float * gi,
16 std::ptrdiff_t gi_istep,
17 std::ptrdiff_t gi_jstep,
18 float * gj,
19 std::ptrdiff_t gj_istep,
20 std::ptrdiff_t gj_jstep,
21 unsigned ni,
22 unsigned nj)
23 {
24 const unsigned char * s_data = src;
25 float * gi_data = gi;
26 float * gj_data = gj;
27
28 if (ni == 0 || nj == 0)
29 return;
30 if (ni == 1)
31 {
32 // Zero the elements in the column
33 for (unsigned j = 0; j < nj; ++j)
34 {
35 *gi_data = 0;
36 *gj_data = 0;
37 gi_data += gi_jstep;
38 gj_data += gj_jstep;
39 }
40 return;
41 }
42 if (nj == 1)
43 {
44 // Zero the elements in the column
45 for (unsigned i = 0; i < ni; ++i)
46 {
47 *gi_data = 0;
48 *gj_data = 0;
49 gi_data += gi_istep;
50 gj_data += gj_istep;
51 }
52 return;
53 }
54
55 // Compute relative grid positions
56 // o1 o2 o3
57 // o4 o5
58 // o6 o7 o8
59 const std::ptrdiff_t o1 = s_jstep - s_istep;
60 const std::ptrdiff_t o2 = s_jstep;
61 const std::ptrdiff_t o3 = s_istep + s_jstep;
62 const std::ptrdiff_t o4 = -s_istep;
63 const std::ptrdiff_t o5 = s_istep;
64 const std::ptrdiff_t o6 = -s_istep - s_jstep;
65 const std::ptrdiff_t o7 = -s_jstep;
66 const std::ptrdiff_t o8 = s_istep - s_jstep;
67
68 const unsigned ni1 = ni - 1;
69 const unsigned nj1 = nj - 1;
70
71 s_data += s_istep + s_jstep;
72 gi_data += gi_jstep;
73 gj_data += gj_jstep;
74
75 for (unsigned j = 1; j < nj1; ++j)
76 {
77 const unsigned char * s = s_data;
78 float * pgi = gi_data;
79 float * pgj = gj_data;
80
81 // Zero the first elements in the rows
82 *pgi = 0;
83 pgi += gi_istep;
84 *pgj = 0;
85 pgj += gj_istep;
86
87 for (unsigned i = 1; i < ni1; ++i)
88 {
89 // Compute gradient in i
90 // Note: Multiply each element individually
91 // to ensure conversion to float before addition
92 *pgi = (0.125f * s[o3] + 0.25f * s[o5] + 0.125f * s[o8]) - (0.125f * s[o1] + 0.25f * s[o4] + 0.125f * s[o6]);
93 // Compute gradient in j
94 *pgj = (0.125f * s[o1] + 0.25f * s[o2] + 0.125f * s[o3]) - (0.125f * s[o6] + 0.25f * s[o7] + 0.125f * s[o8]);
95
96 s += s_istep;
97 pgi += gi_istep;
98 pgj += gj_istep;
99 }
100
101 // Zero the last elements in the rows
102 *pgi = 0;
103 *pgj = 0;
104
105 // Move to next row
106 s_data += s_jstep;
107 gi_data += gi_jstep;
108 gj_data += gj_jstep;
109 }
110
111 // Zero the first and last rows
112 for (unsigned i = 0; i < ni; ++i)
113 {
114 *gi = 0;
115 gi += gi_istep;
116 *gj = 0;
117 gj += gj_istep;
118 *gi_data = 0;
119 gi_data += gi_istep;
120 *gj_data = 0;
121 gj_data += gj_istep;
122 }
123 }
124
125 //: Compute gradients of single plane of 2D data using 3x3 Sobel filters
126 // Computes both i and j gradients of an ni x nj plane of data
127 template <>
128 void
vil_sobel_3x3_1plane(const unsigned char * src,std::ptrdiff_t s_istep,std::ptrdiff_t s_jstep,double * gi,std::ptrdiff_t gi_istep,std::ptrdiff_t gi_jstep,double * gj,std::ptrdiff_t gj_istep,std::ptrdiff_t gj_jstep,unsigned ni,unsigned nj)129 vil_sobel_3x3_1plane(const unsigned char * src,
130 std::ptrdiff_t s_istep,
131 std::ptrdiff_t s_jstep,
132 double * gi,
133 std::ptrdiff_t gi_istep,
134 std::ptrdiff_t gi_jstep,
135 double * gj,
136 std::ptrdiff_t gj_istep,
137 std::ptrdiff_t gj_jstep,
138 unsigned ni,
139 unsigned nj)
140 {
141 const unsigned char * s_data = src;
142 double * gi_data = gi;
143 double * gj_data = gj;
144
145 if (ni == 0 || nj == 0)
146 return;
147 if (ni == 1)
148 {
149 // Zero the elements in the column
150 for (unsigned j = 0; j < nj; ++j)
151 {
152 *gi_data = 0;
153 *gj_data = 0;
154 gi_data += gi_jstep;
155 gj_data += gj_jstep;
156 }
157 return;
158 }
159 if (nj == 1)
160 {
161 // Zero the elements in the column
162 for (unsigned i = 0; i < ni; ++i)
163 {
164 *gi_data = 0;
165 *gj_data = 0;
166 gi_data += gi_istep;
167 gj_data += gj_istep;
168 }
169 return;
170 }
171
172 // Compute relative grid positions
173 // o1 o2 o3
174 // o4 o5
175 // o6 o7 o8
176 const std::ptrdiff_t o1 = s_jstep - s_istep;
177 const std::ptrdiff_t o2 = s_jstep;
178 const std::ptrdiff_t o3 = s_istep + s_jstep;
179 const std::ptrdiff_t o4 = -s_istep;
180 const std::ptrdiff_t o5 = s_istep;
181 const std::ptrdiff_t o6 = -s_istep - s_jstep;
182 const std::ptrdiff_t o7 = -s_jstep;
183 const std::ptrdiff_t o8 = s_istep - s_jstep;
184
185 const unsigned ni1 = ni - 1;
186 const unsigned nj1 = nj - 1;
187
188 s_data += s_istep + s_jstep;
189 gi_data += gi_jstep;
190 gj_data += gj_jstep;
191
192 for (unsigned j = 1; j < nj1; ++j)
193 {
194 const unsigned char * s = s_data;
195 double * pgi = gi_data;
196 double * pgj = gj_data;
197
198 // Zero the first elements in the rows
199 *pgi = 0;
200 pgi += gi_istep;
201 *pgj = 0;
202 pgj += gj_istep;
203
204 for (unsigned i = 1; i < ni1; ++i)
205 {
206 // Compute gradient in i
207 // Note: Multiply each element individually
208 // to ensure conversion to double before addition
209 *pgi = (0.125 * s[o3] + 0.25 * s[o5] + 0.125 * s[o8]) - (0.125 * s[o1] + 0.25 * s[o4] + 0.125 * s[o6]);
210 // Compute gradient in j
211 *pgj = (0.125 * s[o1] + 0.25 * s[o2] + 0.125 * s[o3]) - (0.125 * s[o6] + 0.25 * s[o7] + 0.125 * s[o8]);
212
213 s += s_istep;
214 pgi += gi_istep;
215 pgj += gj_istep;
216 }
217
218 // Zero the last elements in the rows
219 *pgi = 0;
220 *pgj = 0;
221
222 // Move to next row
223 s_data += s_jstep;
224 gi_data += gi_jstep;
225 gj_data += gj_jstep;
226 }
227
228 // Zero the first and last rows
229 for (unsigned i = 0; i < ni; ++i)
230 {
231 *gi = 0;
232 gi += gi_istep;
233 *gj = 0;
234 gj += gj_istep;
235 *gi_data = 0;
236 gi_data += gi_istep;
237 *gj_data = 0;
238 gj_data += gj_istep;
239 }
240 }
241
242 //: Compute gradients of single plane of 2D data using 3x3 Sobel filters
243 // Computes both x and j gradients of an nx x nj plane of data
244 template <>
245 void
vil_sobel_3x3_1plane(const float * src,std::ptrdiff_t s_istep,std::ptrdiff_t s_jstep,float * gi,std::ptrdiff_t gi_istep,std::ptrdiff_t gi_jstep,float * gj,std::ptrdiff_t gj_istep,std::ptrdiff_t gj_jstep,unsigned ni,unsigned nj)246 vil_sobel_3x3_1plane(const float * src,
247 std::ptrdiff_t s_istep,
248 std::ptrdiff_t s_jstep,
249 float * gi,
250 std::ptrdiff_t gi_istep,
251 std::ptrdiff_t gi_jstep,
252 float * gj,
253 std::ptrdiff_t gj_istep,
254 std::ptrdiff_t gj_jstep,
255 unsigned ni,
256 unsigned nj)
257 {
258 const float * s_data = src;
259 float * gi_data = gi;
260 float * gj_data = gj;
261
262 if (ni == 0 || nj == 0)
263 return;
264 if (ni == 1)
265 {
266 // Zero the elements in the column
267 for (unsigned j = 0; j < nj; ++j)
268 {
269 *gi_data = 0;
270 *gj_data = 0;
271 gi_data += gi_jstep;
272 gj_data += gj_jstep;
273 }
274 return;
275 }
276 if (nj == 1)
277 {
278 // Zero the elements in the column
279 for (unsigned i = 0; i < ni; ++i)
280 {
281 *gi_data = 0;
282 *gj_data = 0;
283 gi_data += gi_istep;
284 gj_data += gj_istep;
285 }
286 return;
287 }
288
289 // Compute relative grid positions
290 // o1 o2 o3
291 // o4 o5
292 // o6 o7 o8
293 const std::ptrdiff_t o1 = s_jstep - s_istep;
294 const std::ptrdiff_t o2 = s_jstep;
295 const std::ptrdiff_t o3 = s_istep + s_jstep;
296 const std::ptrdiff_t o4 = -s_istep;
297 const std::ptrdiff_t o5 = s_istep;
298 const std::ptrdiff_t o6 = -s_istep - s_jstep;
299 const std::ptrdiff_t o7 = -s_jstep;
300 const std::ptrdiff_t o8 = s_istep - s_jstep;
301
302 const unsigned ni1 = ni - 1;
303 const unsigned nj1 = nj - 1;
304
305 s_data += s_istep + s_jstep;
306 gi_data += gi_jstep;
307 gj_data += gj_jstep;
308
309 for (unsigned j = 1; j < nj1; ++j)
310 {
311 const float * s = s_data;
312 float * pgi = gi_data;
313 float * pgj = gj_data;
314
315 // Zero the first elements in the rows
316 *pgi = 0;
317 pgi += gi_istep;
318 *pgj = 0;
319 pgj += gj_istep;
320
321 for (unsigned i = 1; i < ni1; ++i)
322 {
323 // Compute gradient in i
324 *pgi = 0.125f * (s[o3] + s[o8] - (s[o1] + s[o6])) + 0.25f * (s[o5] - s[o4]);
325 // Compute gradient in j
326 *pgj = 0.125f * (s[o1] + s[o3] - (s[o6] + s[o8])) + 0.25f * (s[o2] - s[o7]);
327
328 s += s_istep;
329 pgi += gi_istep;
330 pgj += gj_istep;
331 }
332
333 // Zero the last elements in the rows
334 *pgi = 0;
335 *pgj = 0;
336
337 // Move to next row
338 s_data += s_jstep;
339 gi_data += gi_jstep;
340 gj_data += gj_jstep;
341 }
342
343 // Zero the first and last rows
344 for (unsigned i = 0; i < ni; ++i)
345 {
346 *gi = 0;
347 gi += gi_istep;
348 *gj = 0;
349 gj += gj_istep;
350 *gi_data = 0;
351 gi_data += gi_istep;
352 *gj_data = 0;
353 gj_data += gj_istep;
354 }
355 }
356
357 //: Compute gradients of single plane of 2D data using 3x3 Sobel filters
358 // Computes both x and j gradients of an nx x nj plane of data
359 template <>
360 void
vil_sobel_3x3_1plane(const double * src,std::ptrdiff_t s_istep,std::ptrdiff_t s_jstep,double * gi,std::ptrdiff_t gi_istep,std::ptrdiff_t gi_jstep,double * gj,std::ptrdiff_t gj_istep,std::ptrdiff_t gj_jstep,unsigned ni,unsigned nj)361 vil_sobel_3x3_1plane(const double * src,
362 std::ptrdiff_t s_istep,
363 std::ptrdiff_t s_jstep,
364 double * gi,
365 std::ptrdiff_t gi_istep,
366 std::ptrdiff_t gi_jstep,
367 double * gj,
368 std::ptrdiff_t gj_istep,
369 std::ptrdiff_t gj_jstep,
370 unsigned ni,
371 unsigned nj)
372 {
373 const double * s_data = src;
374 double * gi_data = gi;
375 double * gj_data = gj;
376
377 if (ni == 0 || nj == 0)
378 return;
379 if (ni == 1)
380 {
381 // Zero the elements in the column
382 for (unsigned j = 0; j < nj; ++j)
383 {
384 *gi_data = 0;
385 *gj_data = 0;
386 gi_data += gi_jstep;
387 gj_data += gj_jstep;
388 }
389 return;
390 }
391 if (nj == 1)
392 {
393 // Zero the elements in the column
394 for (unsigned i = 0; i < ni; ++i)
395 {
396 *gi_data = 0;
397 *gj_data = 0;
398 gi_data += gi_istep;
399 gj_data += gj_istep;
400 }
401 return;
402 }
403
404 // Compute relative grid positions
405 // o1 o2 o3
406 // o4 o5
407 // o6 o7 o8
408 const std::ptrdiff_t o1 = s_jstep - s_istep;
409 const std::ptrdiff_t o2 = s_jstep;
410 const std::ptrdiff_t o3 = s_istep + s_jstep;
411 const std::ptrdiff_t o4 = -s_istep;
412 const std::ptrdiff_t o5 = s_istep;
413 const std::ptrdiff_t o6 = -s_istep - s_jstep;
414 const std::ptrdiff_t o7 = -s_jstep;
415 const std::ptrdiff_t o8 = s_istep - s_jstep;
416
417 const unsigned ni1 = ni - 1;
418 const unsigned nj1 = nj - 1;
419
420 s_data += s_istep + s_jstep;
421 gi_data += gi_jstep;
422 gj_data += gj_jstep;
423
424 for (unsigned j = 1; j < nj1; ++j)
425 {
426 const double * s = s_data;
427 double * pgi = gi_data;
428 double * pgj = gj_data;
429
430 // Zero the first elements in the rows
431 *pgi = 0;
432 pgi += gi_istep;
433 *pgj = 0;
434 pgj += gj_istep;
435
436 for (unsigned i = 1; i < ni1; ++i)
437 {
438 // Compute gradient in i
439 *pgi = 0.125 * (s[o3] + s[o8] - (s[o1] + s[o6])) + 0.25 * (s[o5] - s[o4]);
440 // Compute gradient in j
441 *pgj = 0.125 * (s[o1] + s[o3] - (s[o6] + s[o8])) + 0.25 * (s[o2] - s[o7]);
442
443 s += s_istep;
444 pgi += gi_istep;
445 pgj += gj_istep;
446 }
447
448 // Zero the last elements in the rows
449 *pgi = 0;
450 *pgj = 0;
451
452 // Move to next row
453 s_data += s_jstep;
454 gi_data += gi_jstep;
455 gj_data += gj_jstep;
456 }
457
458 // Zero the first and last rows
459 for (unsigned i = 0; i < ni; ++i)
460 {
461 *gi = 0;
462 gi += gi_istep;
463 *gj = 0;
464 gj += gj_istep;
465 *gi_data = 0;
466 gi_data += gi_istep;
467 *gj_data = 0;
468 gj_data += gj_istep;
469 }
470 }
471