summaryrefslogtreecommitdiffstats
path: root/all_pairs/source/susan/susan.c
diff options
context:
space:
mode:
Diffstat (limited to 'all_pairs/source/susan/susan.c')
-rw-r--r--all_pairs/source/susan/susan.c2014
1 files changed, 2014 insertions, 0 deletions
diff --git a/all_pairs/source/susan/susan.c b/all_pairs/source/susan/susan.c
new file mode 100644
index 0000000..fb990b1
--- /dev/null
+++ b/all_pairs/source/susan/susan.c
@@ -0,0 +1,2014 @@
1/**********************************************************************\
2
3 SUSAN Version 2l by Stephen Smith
4 Oxford Centre for Functional Magnetic Resonance Imaging of the Brain,
5 Department of Clinical Neurology, Oxford University, Oxford, UK
6 (Previously in Computer Vision and Image Processing Group - now
7 Computer Vision and Electro Optics Group - DERA Chertsey, UK)
8 Email: steve@fmrib.ox.ac.uk
9 WWW: http://www.fmrib.ox.ac.uk/~steve
10
11 (C) Crown Copyright (1995-1999), Defence Evaluation and Research Agency,
12 Farnborough, Hampshire, GU14 6TD, UK
13 DERA WWW site:
14 http://www.dera.gov.uk/
15 DERA Computer Vision and Electro Optics Group WWW site:
16 http://www.dera.gov.uk/imageprocessing/dera/group_home.html
17 DERA Computer Vision and Electro Optics Group point of contact:
18 Dr. John Savage, jtsavage@dera.gov.uk, +44 1344 633203
19
20 A UK patent has been granted: "Method for digitally processing
21 images to determine the position of edges and/or corners therein for
22 guidance of unmanned vehicle", UK Patent 2272285. Proprietor:
23 Secretary of State for Defence, UK. 15 January 1997
24
25 This code is issued for research purposes only and remains the
26 property of the UK Secretary of State for Defence. This code must
27 not be passed on without this header information being kept
28 intact. This code must not be sold.
29
30\**********************************************************************/
31
32/**********************************************************************\
33
34 SUSAN Version 2l
35 SUSAN = Smallest Univalue Segment Assimilating Nucleus
36
37 Email: steve@fmrib.ox.ac.uk
38 WWW: http://www.fmrib.ox.ac.uk/~steve
39
40 Related paper:
41 @article{Smith97,
42 author = "Smith, S.M. and Brady, J.M.",
43 title = "{SUSAN} - A New Approach to Low Level Image Processing",
44 journal = "Int. Journal of Computer Vision",
45 pages = "45--78",
46 volume = "23",
47 number = "1",
48 month = "May",
49 year = 1997}
50
51 To be registered for automatic (bug) updates of SUSAN, send an email.
52
53 Compile with:
54 gcc -O4 -o susan susan2l.c -lm
55
56 See following section for different machine information. Please
57 report any bugs (and fixes). There are a few optional changes that
58 can be made in the "defines" section which follows shortly.
59
60 Usage: type "susan" to get usage. Only PGM format files can be input
61 and output. Utilities such as the netpbm package and XV can be used
62 to convert to and from other formats. Any size of image can be
63 processed.
64
65 This code is written using an emacs folding mode, making moving
66 around the different sections very easy. This is why there are
67 various marks within comments and why comments are indented.
68
69
70 SUSAN QUICK:
71
72 This version of the SUSAN corner finder does not do all the
73 false-corner suppression and thus is faster and produced some false
74 positives, particularly on strong edges. However, because there are
75 less stages involving thresholds etc., the corners that are
76 correctly reported are usually more stable than those reported with
77 the full algorithm. Thus I recommend at least TRYING this algorithm
78 for applications where stability is important, e.g., tracking.
79
80 THRESHOLDS:
81
82 There are two thresholds which can be set at run-time. These are the
83 brightness threshold (t) and the distance threshold (d).
84
85 SPATIAL CONTROL: d
86
87 In SUSAN smoothing d controls the size of the Gaussian mask; its
88 default is 4.0. Increasing d gives more smoothing. In edge finding,
89 a fixed flat mask is used, either 37 pixels arranged in a "circle"
90 (default), or a 3 by 3 mask which gives finer detail. In corner
91 finding, only the larger 37 pixel mask is used; d is not
92 variable. In smoothing, the flat 3 by 3 mask can be used instead of
93 a larger Gaussian mask; this gives low smoothing and fast operation.
94
95 BRIGHTNESS CONTROL: t
96
97 In all three algorithms, t can be varied (default=20); this is the
98 main threshold to be varied. It determines the maximum difference in
99 greylevels between two pixels which allows them to be considered
100 part of the same "region" in the image. Thus it can be reduced to
101 give more edges or corners, i.e. to be more sensitive, and vice
102 versa. In smoothing, reducing t gives less smoothing, and vice
103 versa. Set t=10 for the test image available from the SUSAN web
104 page.
105
106 ITERATIONS:
107
108 With SUSAN smoothing, more smoothing can also be obtained by
109 iterating the algorithm several times. This has a different effect
110 from varying d or t.
111
112 FIXED MASKS:
113
114 37 pixel mask: ooo 3 by 3 mask: ooo
115 ooooo ooo
116 ooooooo ooo
117 ooooooo
118 ooooooo
119 ooooo
120 ooo
121
122 CORNER ATTRIBUTES dx, dy and I
123 (Only read this if you are interested in the C implementation or in
124 using corner attributes, e.g., for corner matching)
125
126 Corners reported in the corner list have attributes associated with
127 them as well as positions. This is useful, for example, when
128 attempting to match corners from one image to another, as these
129 attributes can often be fairly unchanged between images. The
130 attributes are dx, dy and I. I is the value of image brightness at
131 the position of the corner. In the case of susan_corners_quick, dx
132 and dy are the first order derivatives (differentials) of the image
133 brightness in the x and y directions respectively, at the position
134 of the corner. In the case of normal susan corner finding, dx and dy
135 are scaled versions of the position of the centre of gravity of the
136 USAN with respect to the centre pixel (nucleus).
137
138 BRIGHTNESS FUNCTION LUT IMPLEMENTATION:
139 (Only read this if you are interested in the C implementation)
140
141 The SUSAN brightness function is implemented as a LUT
142 (Look-Up-Table) for speed. The resulting pointer-based code is a
143 little hard to follow, so here is a brief explanation. In
144 setup_brightness_lut() the LUT is setup. This mallocs enough space
145 for *bp and then repositions the pointer to the centre of the
146 malloced space. The SUSAN function e^-(x^6) or e^-(x^2) is
147 calculated and converted to a uchar in the range 0-100, for all
148 possible image brightness differences (including negative
149 ones). Thus bp[23] is the output for a brightness difference of 23
150 greylevels. In the SUSAN algorithms this LUT is used as follows:
151
152 p=in + (i-3)*x_size + j - 1;
153 p points to the first image pixel in the circular mask surrounding
154 point (x,y).
155
156 cp=bp + in[i*x_size+j];
157 cp points to a position in the LUT corresponding to the brightness
158 of the centre pixel (x,y).
159
160 now for every pixel within the mask surrounding (x,y),
161 n+=*(cp-*p++);
162 the brightness difference function is found by moving the cp pointer
163 down by an amount equal to the value of the pixel pointed to by p,
164 thus subtracting the two brightness values and performing the
165 exponential function. This value is added to n, the running USAN
166 area.
167
168 in SUSAN smoothing, the variable height mask is implemented by
169 multiplying the above by the moving mask pointer, reset for each new
170 centre pixel.
171 tmp = *dpt++ * *(cp-brightness);
172
173\**********************************************************************/
174
175/**********************************************************************\
176
177 Success has been reported with the following:
178
179 MACHINE OS COMPILER
180
181 Sun 4.1.4 bundled C, gcc
182
183 Next
184
185 SGI IRIX SGI cc
186
187 DEC Unix V3.2+
188
189 IBM RISC AIX gcc
190
191 PC Borland 5.0
192
193 PC Linux gcc-2.6.3
194
195 PC Win32 Visual C++ 4.0 (Console Application)
196
197 PC Win95 Visual C++ 5.0 (Console Application)
198 Thanks to Niu Yongsheng <niuysbit@163.net>:
199 Use the FOPENB option below
200
201 PC DOS djgpp gnu C
202 Thanks to Mark Pettovello <mpettove@umdsun2.umd.umich.edu>:
203 Use the FOPENB option below
204
205 HP HP-UX bundled cc
206 Thanks to Brian Dixon <briand@hpcvsgen.cv.hp.com>:
207 in ksh:
208 export CCOPTS="-Aa -D_HPUX_SOURCE | -lM"
209 cc -O3 -o susan susan2l.c
210
211\**********************************************************************/
212
213/**********************************************************************\
214
215 SUSAN Version 2l, 12/2/99
216 Changed GNUDOS option to FOPENB.
217 (Thanks to Niu Yongsheng <niuysbit@163.net>.)
218 Took out redundant "sq=sq/2;".
219
220 SUSAN Version 2k, 19/8/98:
221 In corner finding:
222 Changed if(yy<sq) {...} else if(xx<sq) {...} to
223 if(yy<xx) {...} else {...}
224 (Thanks to adq@cim.mcgill.edu - Alain Domercq.)
225
226 SUSAN Version 2j, 22/10/97:
227 Fixed (mask_size>x_size) etc. tests in smoothing.
228 Added a couple of free() calls for cgx and cgy.
229 (Thanks to geoffb@ucs.ed.ac.uk - Geoff Browitt.)
230
231 SUSAN Version 2i, 21/7/97:
232 Added information about corner attributes.
233
234 SUSAN Version 2h, 16/12/96:
235 Added principle (initial enhancement) option.
236
237 SUSAN Version 2g, 2/7/96:
238 Minor superficial changes to code.
239
240 SUSAN Version 2f, 16/1/96:
241 Added GNUDOS option (now called FOPENB; see options below).
242
243 SUSAN Version 2e, 9/1/96:
244 Added -b option.
245 Fixed 1 pixel horizontal offset error for drawing edges.
246
247 SUSAN Version 2d, 27/11/95:
248 Fixed loading of certain PGM files in get_image (again!)
249
250 SUSAN Version 2c, 22/11/95:
251 Fixed loading of certain PGM files in get_image.
252 (Thanks to qu@San-Jose.ate.slb.com - Gongyuan Qu.)
253
254 SUSAN Version 2b, 9/11/95:
255 removed "z==" error in edges routines.
256
257 SUSAN Version 2a, 6/11/95:
258 Removed a few unnecessary variable declarations.
259 Added different machine information.
260 Changed "header" in get_image to char.
261
262 SUSAN Version 2, 1/11/95: first combined version able to take any
263 image sizes.
264
265 SUSAN "Versions 1", circa 1992: the various SUSAN algorithms were
266 developed during my doctorate within different programs and for
267 fixed image sizes. The algorithms themselves are virtually unaltered
268 between "versions 1" and the combined program, version 2.
269
270\**********************************************************************/
271
272#include "../extra.h"
273#include "wcclibm.h"
274#include "wccfile.h"
275#include "wccmalloc.h"
276
277#define EXP_A 184
278#define EXP_C 16249
279
280float susan_expf(float y)
281{
282 union
283 {
284 float d;
285 struct
286 {
287#ifdef LITTLE_ENDIAN
288 short j, i;
289#else
290 short i, j;
291#endif
292 } n;
293 } eco;
294 eco.n.i = EXP_A*(y) + (EXP_C);
295 eco.n.j = 0;
296 return eco.d;
297}
298
299float susan_sqrtf(float val)
300{
301 float x = val/10;
302
303 float dx;
304
305 float diff;
306 float min_tol = 0.00001f;
307
308 int i, flag;
309
310
311 flag = 0;
312 if (val == 0 ) x = 0;
313 else {
314 for (i=1;i<20;i++)
315 {
316 if (!flag) {
317 dx = (val - (x*x)) / (2.0f * x);
318 x = x + dx;
319 diff = val - (x*x);
320 if (fabs(diff) <= min_tol) flag = 1;
321 }
322 }
323 }
324 return (x);
325}
326
327/* ********** Optional settings */
328
329typedef int TOTAL_TYPE;
330
331#define SEVEN_SUPP /* size for non-max corner suppression; SEVEN_SUPP or FIVE_SUPP */
332#define MAX_CORNERS 15000 /* max corners per frame */
333
334#define FTOI(a) ( (a) < 0 ? ((int)(a-0.5)) : ((int)(a+0.5)) )
335#define abs(a) ( (a) < 0 ? -a : a )
336typedef unsigned char uchar;
337typedef struct {int x,y,info, dx, dy, I;} corner_rep;
338typedef corner_rep CORNER_LIST[MAX_CORNERS];
339
340extern char susan_input[7292];
341struct wccFILE susan_file;
342float susan_dt;
343int susan_bt;
344int susan_principle_conf;
345int susan_thin_post_proc;
346int susan_three_by_three;
347int susan_drawing_mode;
348int susan_susan_quick;
349int susan_max_no_corners;
350int susan_max_no_edges;
351
352int susan_getint( struct wccFILE *fd )
353{
354 int c, i;
355 char dummy[10000];
356
357 c = susan_wccfgetc(fd);
358 _Pragma( "loopbound min 1 max 4" )
359 while (1) { /* find next integer */
360 if (c=='#') /* if we're at a comment, read to end of line */
361 susan_wccfgets(dummy,9000,fd);
362 if (c==EOFX) {
363 /* "Image is not binary PGM." */
364 }
365 if (c>='0' && c<='9')
366 break; /* found what we were looking for */
367 c = susan_wccfgetc(fd);
368 }
369
370 /* we're at the start of a number, continue until we hit a non-number */
371 i = 0;
372 _Pragma( "loopbound min 2 max 3" )
373 while (1) {
374 i = (i*10) + (c - '0');
375 c = susan_wccfgetc(fd);
376 if (c==EOFX) return (i);
377 if (c<'0' || c>'9') break;
378 }
379
380 return (i);
381}
382
383
384void susan_get_image( struct wccFILE *fd,
385 unsigned char **in, int *x_size, int *y_size )
386{
387 char header [100];
388
389 susan_wccfseek( fd, 0, WCCSEEK_SET );
390
391 /* {{{ read header */
392
393 header[0]=susan_wccfgetc( fd );
394 header[1]=susan_wccfgetc( fd );
395 if(!(header[0]=='P' && header[1]=='5')) {
396 /* "Image does not have binary PGM header." */
397 }
398
399 *x_size = susan_getint(fd);
400 *y_size = susan_getint(fd);
401 // dummy read
402 susan_getint(fd);
403
404/* }}} */
405
406 *in = (uchar *) susan_wccmalloc(*x_size * *y_size);
407
408 if (susan_wccfread(*in,1,*x_size * *y_size,fd) == 0) {
409 /* "Image is wrong size.\n" */
410 }
411}
412
413
414void susan_put_image( int x_size, int y_size )
415{
416 int i;
417 _Pragma( "loopbound min 7220 max 7220" )
418 for ( i = 0; i < x_size*y_size; i++ ) {
419 ;
420 }
421}
422
423
424void susan_int_to_uchar( char *r, uchar *in, int size )
425{
426 int i, max_r=r[0], min_r=r[0];
427
428 _Pragma( "loopbound min 0 max 0" )
429 for (i=0; i<size; i++) {
430 if ( r[i] > max_r )
431 max_r=r[i];
432 if ( r[i] < min_r )
433 min_r=r[i];
434 }
435
436 if (max_r == min_r) {
437 /* Should not occur in TACLeBench. */
438 _Pragma( "loopbound min 0 max 0" )
439 for (i=0; i<size; i++) {
440 in[i] = (uchar)(0);
441 }
442
443 return;
444 }
445 max_r-=min_r;
446
447 _Pragma( "loopbound min 0 max 0" )
448 for (i=0; i<size; i++) {
449 in[i] = (uchar)((int)((int)(r[i]-min_r)*255)/max_r);
450 }
451}
452
453
454void susan_setup_brightness_lut( uchar **bp, int thresh, int form )
455{
456 int k;
457 float temp;
458
459 *bp=(unsigned char *)susan_wccmalloc(516);
460 *bp=*bp+258;
461
462 _Pragma( "loopbound min 513 max 513" )
463 for(k=-256;k<257;k++) {
464 temp=((float)k)/((float)thresh);
465 temp=temp*temp;
466 if (form==6)
467 temp=temp*temp*temp;
468 temp=100.0*susan_expf(-temp);
469 *(*bp+k)= (uchar)temp;
470 }
471}
472
473
474void susan_principle( uchar *in, char *r, uchar *bp,
475 int max_no, int x_size, int y_size )
476{
477 int i, j, n;
478 uchar *p,*cp;
479
480 susan_wccmemset(r,0,x_size * y_size * sizeof(int));
481
482 _Pragma( "loopbound min 0 max 0" )
483 for (i=3;i<y_size-3;i++) {
484 _Pragma( "loopbound min 0 max 0" )
485 for (j=3;j<x_size-3;j++) {
486 n=100;
487 p=in + (i-3)*x_size + j - 1;
488 cp=bp + in[i*x_size+j];
489
490 n+=*(cp-*p++);
491 n+=*(cp-*p++);
492 n+=*(cp-*p);
493 p+=x_size-3;
494
495 n+=*(cp-*p++);
496 n+=*(cp-*p++);
497 n+=*(cp-*p++);
498 n+=*(cp-*p++);
499 n+=*(cp-*p);
500 p+=x_size-5;
501
502 n+=*(cp-*p++);
503 n+=*(cp-*p++);
504 n+=*(cp-*p++);
505 n+=*(cp-*p++);
506 n+=*(cp-*p++);
507 n+=*(cp-*p++);
508 n+=*(cp-*p);
509 p+=x_size-6;
510
511 n+=*(cp-*p++);
512 n+=*(cp-*p++);
513 n+=*(cp-*p);
514 p+=2;
515 n+=*(cp-*p++);
516 n+=*(cp-*p++);
517 n+=*(cp-*p);
518 p+=x_size-6;
519
520 n+=*(cp-*p++);
521 n+=*(cp-*p++);
522 n+=*(cp-*p++);
523 n+=*(cp-*p++);
524 n+=*(cp-*p++);
525 n+=*(cp-*p++);
526 n+=*(cp-*p);
527 p+=x_size-5;
528
529 n+=*(cp-*p++);
530 n+=*(cp-*p++);
531 n+=*(cp-*p++);
532 n+=*(cp-*p++);
533 n+=*(cp-*p);
534 p+=x_size-3;
535
536 n+=*(cp-*p++);
537 n+=*(cp-*p++);
538 n+=*(cp-*p);
539
540 if (n<=max_no)
541 r[i*x_size+j] = max_no - n;
542 }
543 }
544}
545
546
547void susan_principle_small( uchar *in, char *r, uchar *bp,
548 int max_no, int x_size, int y_size )
549{
550 int i, j, n;
551 uchar *p,*cp;
552
553 susan_wccmemset(r,0,x_size * y_size * sizeof(int));
554
555 _Pragma( "loopbound min 0 max 0" )
556 for (i=1;i<y_size-1;i++) {
557 _Pragma( "loopbound min 0 max 0" )
558 for (j=1;j<x_size-1;j++) {
559 n=100;
560 p=in + (i-1)*x_size + j - 1;
561 cp=bp + in[i*x_size+j];
562
563 n+=*(cp-*p++);
564 n+=*(cp-*p++);
565 n+=*(cp-*p);
566 p+=x_size-2;
567
568 n+=*(cp-*p);
569 p+=2;
570 n+=*(cp-*p);
571 p+=x_size-2;
572
573 n+=*(cp-*p++);
574 n+=*(cp-*p++);
575 n+=*(cp-*p);
576
577 if (n<=max_no)
578 r[i*x_size+j] = max_no - n;
579 }
580 }
581}
582
583
584uchar susan_median( uchar *in, int i, int j, int x_size )
585{
586 int p[8],k,l,tmp;
587
588 p[0]=in[(i-1)*x_size+j-1];
589 p[1]=in[(i-1)*x_size+j ];
590 p[2]=in[(i-1)*x_size+j+1];
591 p[3]=in[(i )*x_size+j-1];
592 p[4]=in[(i )*x_size+j+1];
593 p[5]=in[(i+1)*x_size+j-1];
594 p[6]=in[(i+1)*x_size+j ];
595 p[7]=in[(i+1)*x_size+j+1];
596
597 _Pragma( "loopbound min 0 max 0" )
598 for(k=0; k<7; k++) {
599 _Pragma( "loopbound min 0 max 0" )
600 for(l=0; l<(7-k); l++) {
601 if (p[l]>p[l+1]) {
602 tmp=p[l];
603 p[l]=p[l+1];
604 p[l+1]=tmp;
605 }
606 }
607 }
608
609 return( (p[3]+p[4]) / 2 );
610}
611
612
613/* this enlarges "in" so that borders can be dealt with easily */
614void susan_enlarge( uchar **in, uchar *tmp_image,
615 int *x_size, int *y_size, int border )
616{
617 int i, j;
618
619 _Pragma( "loopbound min 95 max 95" )
620 for(i=0; i<*y_size; i++) { /* copy *in into tmp_image */
621 susan_wccmemcpy( tmp_image+(i+border)*(*x_size+2*border)+border,
622 *in+i* *x_size, *x_size);
623 }
624
625 _Pragma( "loopbound min 7 max 7" )
626 for(i=0; i<border; i++) { /* copy top and bottom rows; invert as many as necessary */
627 susan_wccmemcpy( tmp_image+(border-1-i)*(*x_size+2*border)+border,
628 *in+i* *x_size,*x_size);
629 susan_wccmemcpy( tmp_image+(*y_size+border+i)*(*x_size+2*border)+border,
630 *in+(*y_size-i-1)* *x_size,*x_size);
631 }
632
633 _Pragma( "loopbound min 7 max 7" )
634 for(i=0; i<border; i++) { /* copy left and right columns */
635 _Pragma( "loopbound min 109 max 109" )
636 for(j=0; j<*y_size+2*border; j++) {
637 tmp_image[j*(*x_size+2*border)+border-1-i]=tmp_image[j*(*x_size+2*border)+border+i];
638 tmp_image[j*(*x_size+2*border)+ *x_size+border+i]=tmp_image[j*(*x_size+2*border)+ *x_size+border-1-i];
639 }
640 }
641
642 *x_size+=2*border; /* alter image size */
643 *y_size+=2*border;
644 *in=tmp_image; /* repoint in */
645}
646
647
648void susan_smoothing( int three_by_three, uchar *in, float dt,
649 int x_size, int y_size, uchar *bp )
650{
651 float temp;
652 int n_max, increment, mask_size,
653 i,j,x,y,area,brightness,tmp,centre;
654 uchar *ip, *dp, *dpt, *cp, *out=in,
655 *tmp_image;
656 TOTAL_TYPE total;
657
658 /* {{{ setup larger image and border sizes */
659
660 if (three_by_three==0)
661 mask_size = ((int)(1.5 * dt)) + 1;
662 else
663 mask_size = 1;
664
665 if ( dt>15 ) {
666 /* "Distance_thresh too big for integer arithmetic." */
667 // Either reduce it to <=15 or recompile with variable "total"
668 // as a float: see top "defines" section.
669 }
670
671 if ( (2*mask_size+1>x_size) || (2*mask_size+1>y_size) ) {
672 /* "Mask size too big for image." */
673 }
674
675 tmp_image = (uchar *)susan_wccmalloc( (x_size+mask_size*2) * (y_size+mask_size*2) );
676 susan_enlarge(&in,tmp_image,&x_size,&y_size,mask_size);
677
678 if (three_by_three==0) {
679 /* large Gaussian masks */
680 /* {{{ setup distance lut */
681
682 n_max = (mask_size*2) + 1;
683
684 increment = x_size - n_max;
685
686 dp = (unsigned char *)susan_wccmalloc(n_max*n_max);
687 dpt = dp;
688 temp = -(dt*dt);
689
690 _Pragma( "loopbound min 15 max 15" )
691 for(i=-mask_size; i<=mask_size; i++) {
692 _Pragma( "loopbound min 15 max 15" )
693 for(j=-mask_size; j<=mask_size; j++) {
694 x = (int) (100.0 * susan_expf( ((float)((i*i)+(j*j))) / temp ));
695 *dpt++ = (unsigned char)x;
696 }
697 }
698
699 /* {{{ main section */
700
701 _Pragma( "loopbound min 95 max 95" )
702 for (i=mask_size;i<y_size-mask_size;i++) {
703 _Pragma( "loopbound min 76 max 76" )
704 for (j=mask_size;j<x_size-mask_size;j++) {
705 area = 0;
706 total = 0;
707 dpt = dp;
708 ip = in + ((i-mask_size)*x_size) + j - mask_size;
709 centre = in[i*x_size+j];
710 cp = bp + centre;
711 _Pragma( "loopbound min 15 max 15" )
712 for(y=-mask_size; y<=mask_size; y++) {
713 _Pragma( "loopbound min 15 max 15" )
714 for(x=-mask_size; x<=mask_size; x++) {
715 brightness = *ip++;
716 tmp = *dpt++ * *(cp-brightness);
717 area += tmp;
718 total += tmp * brightness;
719 }
720 ip += increment;
721 }
722 tmp = area-10000;
723 if (tmp==0)
724 *out++=susan_median(in,i,j,x_size);
725 else
726 *out++=((total-(centre*10000))/tmp);
727 }
728 }
729 } else {
730 /* 3x3 constant mask */
731
732 /* {{{ main section */
733
734 _Pragma( "loopbound min 15 max 15" )
735 for (i=1;i<y_size-1;i++) {
736 _Pragma( "loopbound min 15 max 15" )
737 for (j=1;j<x_size-1;j++) {
738 area = 0;
739 total = 0;
740 ip = in + ((i-1)*x_size) + j - 1;
741 centre = in[i*x_size+j];
742 cp = bp + centre;
743
744 brightness=*ip++; tmp=*(cp-brightness); area += tmp; total += tmp * brightness;
745 brightness=*ip++; tmp=*(cp-brightness); area += tmp; total += tmp * brightness;
746 brightness=*ip; tmp=*(cp-brightness); area += tmp; total += tmp * brightness;
747 ip += x_size-2;
748 brightness=*ip++; tmp=*(cp-brightness); area += tmp; total += tmp * brightness;
749 brightness=*ip++; tmp=*(cp-brightness); area += tmp; total += tmp * brightness;
750 brightness=*ip; tmp=*(cp-brightness); area += tmp; total += tmp * brightness;
751 ip += x_size-2;
752 brightness=*ip++; tmp=*(cp-brightness); area += tmp; total += tmp * brightness;
753 brightness=*ip++; tmp=*(cp-brightness); area += tmp; total += tmp * brightness;
754 brightness=*ip; tmp=*(cp-brightness); area += tmp; total += tmp * brightness;
755
756 tmp = area-100;
757 if (tmp==0)
758 *out++=susan_median(in,i,j,x_size);
759 else
760 *out++=(total-(centre*100))/tmp;
761 }
762 }
763 }
764}
765
766
767void susan_edge_draw( uchar *in, uchar *mid,
768 int x_size, int y_size, int drawing_mode )
769{
770 int i;
771 uchar *inp, *midp;
772
773 if (drawing_mode==0) {
774 /* mark 3x3 white block around each edge point */
775 midp=mid;
776 _Pragma( "loopbound min 7220 max 7220" )
777 for (i=0; i<x_size*y_size; i++) {
778 if (*midp<8) {
779 inp = in + (midp - mid) - x_size - 1;
780 *inp++=255; *inp++=255; *inp=255; inp+=x_size-2;
781 *inp++=255; inp++; *inp=255; inp+=x_size-2;
782 *inp++=255; *inp++=255; *inp=255;
783 }
784 midp++;
785 }
786 }
787
788 /* now mark 1 black pixel at each edge point */
789 midp=mid;
790 _Pragma( "loopbound min 7220 max 7220" )
791 for (i=0; i<x_size*y_size; i++) {
792 if (*midp<8)
793 *(in + (midp - mid)) = 0;
794 midp++;
795 }
796}
797
798
799void susan_thin( char *r, uchar *mid, int x_size, int y_size )
800{
801 int l[9], centre,
802 b01, b12, b21, b10,
803 p1, p2, p3, p4,
804 b00, b02, b20, b22,
805 m, n, a, b, x, y, i, j;
806 uchar *mp;
807
808 _Pragma( "loopbound min 87 max 87" )
809 for (i=4;i<y_size-4;i++) {
810 _Pragma( "loopbound min 68 max 68" )
811 for (j=4;j<x_size-4;j++) {
812 if (mid[i*x_size+j]<8) {
813 centre = r[i*x_size+j];
814 /* {{{ count number of neighbours */
815
816 mp=mid + (i-1)*x_size + j-1;
817
818 n = (*mp<8) +
819 (*(mp+1)<8) +
820 (*(mp+2)<8) +
821 (*(mp+x_size)<8) +
822 (*(mp+x_size+2)<8) +
823 (*(mp+x_size+x_size)<8) +
824 (*(mp+x_size+x_size+1)<8) +
825 (*(mp+x_size+x_size+2)<8);
826
827 /* {{{ n==0 no neighbours - remove point */
828
829 if (n==0)
830 mid[i*x_size+j]=100;
831
832 /* {{{ n==1 - extend line if I can */
833
834 /* extension is only allowed a few times - the value of mid is used to control this */
835
836 if ( (n==1) && (mid[i*x_size+j]<6) ) {
837 /* find maximum neighbour weighted in direction opposite the
838 neighbour already present. e.g.
839 have: O O O weight r by 0 2 3
840 X X O 0 0 4
841 O O O 0 2 3 */
842
843 l[0]=r[(i-1)*x_size+j-1]; l[1]=r[(i-1)*x_size+j]; l[2]=r[(i-1)*x_size+j+1];
844 l[3]=r[(i )*x_size+j-1]; l[4]=0; l[5]=r[(i )*x_size+j+1];
845 l[6]=r[(i+1)*x_size+j-1]; l[7]=r[(i+1)*x_size+j]; l[8]=r[(i+1)*x_size+j+1];
846
847 if (mid[(i-1)*x_size+j-1]<8) { l[0]=0; l[1]=0; l[3]=0; l[2]*=2;
848 l[6]*=2; l[5]*=3; l[7]*=3; l[8]*=4; }
849 else { if (mid[(i-1)*x_size+j]<8) { l[1]=0; l[0]=0; l[2]=0; l[3]*=2;
850 l[5]*=2; l[6]*=3; l[8]*=3; l[7]*=4; }
851 else { if (mid[(i-1)*x_size+j+1]<8) { l[2]=0; l[1]=0; l[5]=0; l[0]*=2;
852 l[8]*=2; l[3]*=3; l[7]*=3; l[6]*=4; }
853 else { if (mid[(i)*x_size+j-1]<8) { l[3]=0; l[0]=0; l[6]=0; l[1]*=2;
854 l[7]*=2; l[2]*=3; l[8]*=3; l[5]*=4; }
855 else { if (mid[(i)*x_size+j+1]<8) { l[5]=0; l[2]=0; l[8]=0; l[1]*=2;
856 l[7]*=2; l[0]*=3; l[6]*=3; l[3]*=4; }
857 else { if (mid[(i+1)*x_size+j-1]<8) { l[6]=0; l[3]=0; l[7]=0; l[0]*=2;
858 l[8]*=2; l[1]*=3; l[5]*=3; l[2]*=4; }
859 else { if (mid[(i+1)*x_size+j]<8) { l[7]=0; l[6]=0; l[8]=0; l[3]*=2;
860 l[5]*=2; l[0]*=3; l[2]*=3; l[1]*=4; }
861 else { if (mid[(i+1)*x_size+j+1]<8) { l[8]=0; l[5]=0; l[7]=0; l[6]*=2;
862 l[2]*=2; l[1]*=3; l[3]*=3; l[0]*=4; } }}}}}}}
863
864 m=0; /* find the highest point */
865 _Pragma( "loopbound min 3 max 3" )
866 for(y=0; y<3; y++)
867 _Pragma( "loopbound min 3 max 3" )
868 for(x=0; x<3; x++)
869 if (l[y+y+y+x]>m) { m=l[y+y+y+x]; a=y; b=x; }
870
871 if (m>0) {
872 if (mid[i*x_size+j]<4)
873 mid[(i+a-1)*x_size+j+b-1] = 4;
874 else
875 mid[(i+a-1)*x_size+j+b-1] = mid[i*x_size+j]+1;
876 if ( (a+a+b) < 3 ) { /* need to jump back in image */
877 i+=a-1;
878 j+=b-2;
879 if (i<4) i=4;
880 if (j<4) j=4;
881 }
882 }
883 }
884
885 /* {{{ n==2 */
886
887 if (n==2) {
888 /* put in a bit here to straighten edges */
889 b00 = mid[(i-1)*x_size+j-1]<8; /* corners of 3x3 */
890 b02 = mid[(i-1)*x_size+j+1]<8;
891 b20 = mid[(i+1)*x_size+j-1]<8;
892 b22 = mid[(i+1)*x_size+j+1]<8;
893 if ( ((b00+b02+b20+b22)==2) && ((b00|b22)&(b02|b20))) {
894 /* case: move a point back into line.
895 e.g. X O X CAN become X X X
896 O X O O O O
897 O O O O O O */
898 if (b00) {
899 if (b02) { x=0; y=-1; }
900 else { x=-1; y=0; }
901 } else {
902 if (b02) { x=1; y=0; }
903 else { x=0; y=1; }
904 }
905 if (((float)r[(i+y)*x_size+j+x]/(float)centre) > 0.7) {
906 if ( ( (x==0) && (mid[(i+(2*y))*x_size+j]>7) && (mid[(i+(2*y))*x_size+j-1]>7) && (mid[(i+(2*y))*x_size+j+1]>7) ) ||
907 ( (y==0) && (mid[(i)*x_size+j+(2*x)]>7) && (mid[(i+1)*x_size+j+(2*x)]>7) && (mid[(i-1)*x_size+j+(2*x)]>7) ) ) {
908 mid[(i)*x_size+j]=100;
909 mid[(i+y)*x_size+j+x]=3; /* no jumping needed */
910 }
911 }
912 } else {
913 b01 = mid[(i-1)*x_size+j ]<8;
914 b12 = mid[(i )*x_size+j+1]<8;
915 b21 = mid[(i+1)*x_size+j ]<8;
916 b10 = mid[(i )*x_size+j-1]<8;
917 /* {{{ right angle ends - not currently used */
918
919#ifdef IGNORETHIS
920 if ( (b00&b01)|(b00&b10)|(b02&b01)|(b02&b12)|(b20&b10)|(b20&b21)|(b22&b21)|(b22&b12) ) {
921 /* case; right angle ends. clean up.
922 e.g.; X X O CAN become X X O
923 O X O O O O
924 O O O O O O */
925 if ( ((b01)&(mid[(i-2)*x_size+j-1]>7)&(mid[(i-2)*x_size+j]>7)&(mid[(i-2)*x_size+j+1]>7)&
926 ((b00&((2*r[(i-1)*x_size+j+1])>centre))|(b02&((2*r[(i-1)*x_size+j-1])>centre)))) |
927 ((b10)&(mid[(i-1)*x_size+j-2]>7)&(mid[(i)*x_size+j-2]>7)&(mid[(i+1)*x_size+j-2]>7)&
928 ((b00&((2*r[(i+1)*x_size+j-1])>centre))|(b20&((2*r[(i-1)*x_size+j-1])>centre)))) |
929 ((b12)&(mid[(i-1)*x_size+j+2]>7)&(mid[(i)*x_size+j+2]>7)&(mid[(i+1)*x_size+j+2]>7)&
930 ((b02&((2*r[(i+1)*x_size+j+1])>centre))|(b22&((2*r[(i-1)*x_size+j+1])>centre)))) |
931 ((b21)&(mid[(i+2)*x_size+j-1]>7)&(mid[(i+2)*x_size+j]>7)&(mid[(i+2)*x_size+j+1]>7)&
932 ((b20&((2*r[(i+1)*x_size+j+1])>centre))|(b22&((2*r[(i+1)*x_size+j-1])>centre)))) ) {
933 mid[(i)*x_size+j]=100;
934 if (b10&b20) j-=2;
935 if (b00|b01|b02) { i--; j-=2; }
936 }
937 }
938#endif
939
940 if ( ((b01+b12+b21+b10)==2) && ((b10|b12)&(b01|b21)) &&
941 ((b01&((mid[(i-2)*x_size+j-1]<8)|(mid[(i-2)*x_size+j+1]<8)))|(b10&((mid[(i-1)*x_size+j-2]<8)|(mid[(i+1)*x_size+j-2]<8)))|
942 (b12&((mid[(i-1)*x_size+j+2]<8)|(mid[(i+1)*x_size+j+2]<8)))|(b21&((mid[(i+2)*x_size+j-1]<8)|(mid[(i+2)*x_size+j+1]<8)))) ) {
943 /* case; clears odd right angles.
944 e.g.; O O O becomes O O O
945 X X O X O O
946 O X O O X O */
947 mid[(i)*x_size+j]=100;
948 i--; /* jump back */
949 j-=2;
950 if (i<4) i=4;
951 if (j<4) j=4;
952 }
953 }
954 }
955
956 /* {{{ n>2 the thinning is done here without breaking connectivity */
957
958 if (n>2) {
959 b01 = mid[(i-1)*x_size+j ]<8;
960 b12 = mid[(i )*x_size+j+1]<8;
961 b21 = mid[(i+1)*x_size+j ]<8;
962 b10 = mid[(i )*x_size+j-1]<8;
963 if((b01+b12+b21+b10)>1) {
964 b00 = mid[(i-1)*x_size+j-1]<8;
965 b02 = mid[(i-1)*x_size+j+1]<8;
966 b20 = mid[(i+1)*x_size+j-1]<8;
967 b22 = mid[(i+1)*x_size+j+1]<8;
968 p1 = b00 | b01;
969 p2 = b02 | b12;
970 p3 = b22 | b21;
971 p4 = b20 | b10;
972
973 if( ((p1 + p2 + p3 + p4) - ((b01 & p2)+(b12 & p3)+(b21 & p4)+(b10 & p1))) < 2) {
974 mid[(i)*x_size+j]=100;
975 i--;
976 j-=2;
977 if (i<4) i=4;
978 if (j<4) j=4;
979 }
980 }
981 }
982 }
983 }
984 }
985}
986
987
988void susan_edges( uchar *in, char *r, uchar *mid, uchar *bp,
989 int max_no, int x_size, int y_size )
990{
991 float z;
992 int do_symmetry, i, j, m, n, a, b, x, y, w;
993 uchar c,*p,*cp;
994
995 susan_wccmemset(r,0,x_size * y_size);
996
997 _Pragma( "loopbound min 89 max 89" )
998 for (i=3;i<y_size-3;i++) {
999 _Pragma( "loopbound min 70 max 70" )
1000 for (j=3;j<x_size-3;j++) {
1001 n=100;
1002 p=in + (i-3)*x_size + j - 1;
1003 cp=bp + in[i*x_size+j];
1004
1005 n+=*(cp-*p++);
1006 n+=*(cp-*p++);
1007 n+=*(cp-*p);
1008 p+=x_size-3;
1009
1010 n+=*(cp-*p++);
1011 n+=*(cp-*p++);
1012 n+=*(cp-*p++);
1013 n+=*(cp-*p++);
1014 n+=*(cp-*p);
1015 p+=x_size-5;
1016
1017 n+=*(cp-*p++);
1018 n+=*(cp-*p++);
1019 n+=*(cp-*p++);
1020 n+=*(cp-*p++);
1021 n+=*(cp-*p++);
1022 n+=*(cp-*p++);
1023 n+=*(cp-*p);
1024 p+=x_size-6;
1025
1026 n+=*(cp-*p++);
1027 n+=*(cp-*p++);
1028 n+=*(cp-*p);
1029 p+=2;
1030 n+=*(cp-*p++);
1031 n+=*(cp-*p++);
1032 n+=*(cp-*p);
1033 p+=x_size-6;
1034
1035 n+=*(cp-*p++);
1036 n+=*(cp-*p++);
1037 n+=*(cp-*p++);
1038 n+=*(cp-*p++);
1039 n+=*(cp-*p++);
1040 n+=*(cp-*p++);
1041 n+=*(cp-*p);
1042 p+=x_size-5;
1043
1044 n+=*(cp-*p++);
1045 n+=*(cp-*p++);
1046 n+=*(cp-*p++);
1047 n+=*(cp-*p++);
1048 n+=*(cp-*p);
1049 p+=x_size-3;
1050
1051 n+=*(cp-*p++);
1052 n+=*(cp-*p++);
1053 n+=*(cp-*p);
1054
1055 if (n<=max_no)
1056 r[i*x_size+j] = max_no - n;
1057 }
1058 }
1059
1060 _Pragma( "loopbound min 87 max 87" )
1061 for (i=4;i<y_size-4;i++) {
1062 _Pragma( "loopbound min 68 max 68" )
1063 for (j=4;j<x_size-4;j++) {
1064 if (r[i*x_size+j]>0) {
1065 m=r[i*x_size+j];
1066 n=max_no - m;
1067 cp=bp + in[i*x_size+j];
1068
1069 if (n>600) {
1070 p=in + (i-3)*x_size + j - 1;
1071 x=0;y=0;
1072
1073 c=*(cp-*p++);x-=c;y-=3*c;
1074 c=*(cp-*p++);y-=3*c;
1075 c=*(cp-*p);x+=c;y-=3*c;
1076 p+=x_size-3;
1077
1078 c=*(cp-*p++);x-=2*c;y-=2*c;
1079 c=*(cp-*p++);x-=c;y-=2*c;
1080 c=*(cp-*p++);y-=2*c;
1081 c=*(cp-*p++);x+=c;y-=2*c;
1082 c=*(cp-*p);x+=2*c;y-=2*c;
1083 p+=x_size-5;
1084
1085 c=*(cp-*p++);x-=3*c;y-=c;
1086 c=*(cp-*p++);x-=2*c;y-=c;
1087 c=*(cp-*p++);x-=c;y-=c;
1088 c=*(cp-*p++);y-=c;
1089 c=*(cp-*p++);x+=c;y-=c;
1090 c=*(cp-*p++);x+=2*c;y-=c;
1091 c=*(cp-*p);x+=3*c;y-=c;
1092 p+=x_size-6;
1093
1094 c=*(cp-*p++);x-=3*c;
1095 c=*(cp-*p++);x-=2*c;
1096 c=*(cp-*p);x-=c;
1097 p+=2;
1098 c=*(cp-*p++);x+=c;
1099 c=*(cp-*p++);x+=2*c;
1100 c=*(cp-*p);x+=3*c;
1101 p+=x_size-6;
1102
1103 c=*(cp-*p++);x-=3*c;y+=c;
1104 c=*(cp-*p++);x-=2*c;y+=c;
1105 c=*(cp-*p++);x-=c;y+=c;
1106 c=*(cp-*p++);y+=c;
1107 c=*(cp-*p++);x+=c;y+=c;
1108 c=*(cp-*p++);x+=2*c;y+=c;
1109 c=*(cp-*p);x+=3*c;y+=c;
1110 p+=x_size-5;
1111
1112 c=*(cp-*p++);x-=2*c;y+=2*c;
1113 c=*(cp-*p++);x-=c;y+=2*c;
1114 c=*(cp-*p++);y+=2*c;
1115 c=*(cp-*p++);x+=c;y+=2*c;
1116 c=*(cp-*p);x+=2*c;y+=2*c;
1117 p+=x_size-3;
1118
1119 c=*(cp-*p++);x-=c;y+=3*c;
1120 c=*(cp-*p++);y+=3*c;
1121 c=*(cp-*p);x+=c;y+=3*c;
1122
1123 z = susan_sqrtf((float)((x*x) + (y*y)));
1124 if (z > (0.9*(float)n)) { /* 0.5 */
1125 do_symmetry=0;
1126 if (x==0)
1127 z=1000000.0;
1128 else
1129 z=((float)y) / ((float)x);
1130 if (z < 0) { z=-z; w=-1; }
1131 else w=1;
1132 if (z < 0.5) { /* vert_edge */ a=0; b=1; }
1133 else { if (z > 2.0) { /* hor_edge */ a=1; b=0; }
1134 else { /* diag_edge */ if (w>0) { a=1; b=1; }
1135 else { a=-1; b=1; }}}
1136 if ( (m > r[(i+a)*x_size+j+b]) && (m >= r[(i-a)*x_size+j-b]) &&
1137 (m > r[(i+(2*a))*x_size+j+(2*b)]) && (m >= r[(i-(2*a))*x_size+j-(2*b)]) )
1138 mid[i*x_size+j] = 1;
1139 } else {
1140 do_symmetry=1;
1141 }
1142 } else {
1143 do_symmetry=1;
1144 }
1145
1146 if (do_symmetry==1) {
1147 p=in + (i-3)*x_size + j - 1;
1148 x=0; y=0; w=0;
1149
1150 /* | \
1151 y -x- w
1152 | \ */
1153
1154 c=*(cp-*p++);x+=c;y+=9*c;w+=3*c;
1155 c=*(cp-*p++);y+=9*c;
1156 c=*(cp-*p);x+=c;y+=9*c;w-=3*c;
1157 p+=x_size-3;
1158
1159 c=*(cp-*p++);x+=4*c;y+=4*c;w+=4*c;
1160 c=*(cp-*p++);x+=c;y+=4*c;w+=2*c;
1161 c=*(cp-*p++);y+=4*c;
1162 c=*(cp-*p++);x+=c;y+=4*c;w-=2*c;
1163 c=*(cp-*p);x+=4*c;y+=4*c;w-=4*c;
1164 p+=x_size-5;
1165
1166 c=*(cp-*p++);x+=9*c;y+=c;w+=3*c;
1167 c=*(cp-*p++);x+=4*c;y+=c;w+=2*c;
1168 c=*(cp-*p++);x+=c;y+=c;w+=c;
1169 c=*(cp-*p++);y+=c;
1170 c=*(cp-*p++);x+=c;y+=c;w-=c;
1171 c=*(cp-*p++);x+=4*c;y+=c;w-=2*c;
1172 c=*(cp-*p);x+=9*c;y+=c;w-=3*c;
1173 p+=x_size-6;
1174
1175 c=*(cp-*p++);x+=9*c;
1176 c=*(cp-*p++);x+=4*c;
1177 c=*(cp-*p);x+=c;
1178 p+=2;
1179 c=*(cp-*p++);x+=c;
1180 c=*(cp-*p++);x+=4*c;
1181 c=*(cp-*p);x+=9*c;
1182 p+=x_size-6;
1183
1184 c=*(cp-*p++);x+=9*c;y+=c;w-=3*c;
1185 c=*(cp-*p++);x+=4*c;y+=c;w-=2*c;
1186 c=*(cp-*p++);x+=c;y+=c;w-=c;
1187 c=*(cp-*p++);y+=c;
1188 c=*(cp-*p++);x+=c;y+=c;w+=c;
1189 c=*(cp-*p++);x+=4*c;y+=c;w+=2*c;
1190 c=*(cp-*p);x+=9*c;y+=c;w+=3*c;
1191 p+=x_size-5;
1192
1193 c=*(cp-*p++);x+=4*c;y+=4*c;w-=4*c;
1194 c=*(cp-*p++);x+=c;y+=4*c;w-=2*c;
1195 c=*(cp-*p++);y+=4*c;
1196 c=*(cp-*p++);x+=c;y+=4*c;w+=2*c;
1197 c=*(cp-*p);x+=4*c;y+=4*c;w+=4*c;
1198 p+=x_size-3;
1199
1200 c=*(cp-*p++);x+=c;y+=9*c;w-=3*c;
1201 c=*(cp-*p++);y+=9*c;
1202 c=*(cp-*p);x+=c;y+=9*c;w+=3*c;
1203
1204 if (y==0)
1205 z = 1000000.0;
1206 else
1207 z = ((float)x) / ((float)y);
1208 if (z < 0.5) { /* vertical */ a=0; b=1; }
1209 else { if (z > 2.0) { /* horizontal */ a=1; b=0; }
1210 else { /* diagonal */ if (w>0) { a=-1; b=1; }
1211 else { a=1; b=1; }}}
1212 if ( (m > r[(i+a)*x_size+j+b]) && (m >= r[(i-a)*x_size+j-b]) &&
1213 (m > r[(i+(2*a))*x_size+j+(2*b)]) && (m >= r[(i-(2*a))*x_size+j-(2*b)]) )
1214 mid[i*x_size+j] = 2;
1215 }
1216 }
1217 }
1218 }
1219}
1220
1221
1222void susan_edges_small( uchar *in, char *r, uchar *mid, uchar *bp,
1223 int max_no, int x_size, int y_size )
1224{
1225 float z;
1226 int do_symmetry, i, j, m, n, a, b, x, y, w;
1227 uchar c,*p,*cp;
1228
1229 susan_wccmemset(r,0,x_size * y_size);
1230
1231 _Pragma( "loopbound min 0 max 0" )
1232 for (i=1;i<y_size-1;i++) {
1233 _Pragma( "loopbound min 0 max 0" )
1234 for (j=1;j<x_size-1;j++) {
1235 n=100;
1236 p=in + (i-1)*x_size + j - 1;
1237 cp=bp + in[i*x_size+j];
1238
1239 n+=*(cp-*p++);
1240 n+=*(cp-*p++);
1241 n+=*(cp-*p);
1242 p+=x_size-2;
1243
1244 n+=*(cp-*p);
1245 p+=2;
1246 n+=*(cp-*p);
1247 p+=x_size-2;
1248
1249 n+=*(cp-*p++);
1250 n+=*(cp-*p++);
1251 n+=*(cp-*p);
1252
1253 if (n<=max_no) {
1254 r[i*x_size+j] = max_no - n;
1255 }
1256 }
1257 }
1258
1259 _Pragma( "loopbound min 0 max 0" )
1260 for (i=2;i<y_size-2;i++) {
1261 _Pragma( "loopbound min 0 max 0" )
1262 for (j=2;j<x_size-2;j++) {
1263 if (r[i*x_size+j]>0) {
1264 m=r[i*x_size+j];
1265 n=max_no - m;
1266 cp=bp + in[i*x_size+j];
1267
1268 if (n>250) {
1269 p=in + (i-1)*x_size + j - 1;
1270 x=0;y=0;
1271
1272 c=*(cp-*p++);x-=c;y-=c;
1273 c=*(cp-*p++);y-=c;
1274 c=*(cp-*p);x+=c;y-=c;
1275 p+=x_size-2;
1276
1277 c=*(cp-*p);x-=c;
1278 p+=2;
1279 c=*(cp-*p);x+=c;
1280 p+=x_size-2;
1281
1282 c=*(cp-*p++);x-=c;y+=c;
1283 c=*(cp-*p++);y+=c;
1284 c=*(cp-*p);x+=c;y+=c;
1285
1286 z = susan_sqrtf((float)((x*x) + (y*y)));
1287 if (z > (0.4*(float)n)) { /* 0.6 */
1288 do_symmetry=0;
1289 if (x==0)
1290 z=1000000.0;
1291 else
1292 z=((float)y) / ((float)x);
1293 if (z < 0) { z=-z; w=-1; }
1294 else w=1;
1295 if (z < 0.5) { /* vert_edge */ a=0; b=1; }
1296 else { if (z > 2.0) { /* hor_edge */ a=1; b=0; }
1297 else { /* diag_edge */ if (w>0) { a=1; b=1; }
1298 else { a=-1; b=1; }
1299 }
1300 }
1301 if ( (m > r[(i+a)*x_size+j+b]) && (m >= r[(i-a)*x_size+j-b]) ) {
1302 mid[i*x_size+j] = 1;
1303 }
1304 } else {
1305 do_symmetry=1;
1306 }
1307 } else {
1308 do_symmetry=1;
1309 }
1310
1311 if (do_symmetry==1) {
1312 p=in + (i-1)*x_size + j - 1;
1313 x=0; y=0; w=0;
1314
1315 /* | \
1316 y -x- w
1317 | \ */
1318
1319 c=*(cp-*p++);x+=c;y+=c;w+=c;
1320 c=*(cp-*p++);y+=c;
1321 c=*(cp-*p);x+=c;y+=c;w-=c;
1322 p+=x_size-2;
1323
1324 c=*(cp-*p);x+=c;
1325 p+=2;
1326 c=*(cp-*p);x+=c;
1327 p+=x_size-2;
1328
1329 c=*(cp-*p++);x+=c;y+=c;w-=c;
1330 c=*(cp-*p++);y+=c;
1331 c=*(cp-*p);x+=c;y+=c;w+=c;
1332
1333 if (y==0)
1334 z = 1000000.0;
1335 else
1336 z = ((float)x) / ((float)y);
1337 if (z < 0.5) { /* vertical */ a=0; b=1; }
1338 else { if (z > 2.0) { /* horizontal */ a=1; b=0; }
1339 else { /* diagonal */ if (w>0) { a=-1; b=1; }
1340 else { a=1; b=1; }
1341 }
1342 }
1343 if ( (m > r[(i+a)*x_size+j+b]) && (m >= r[(i-a)*x_size+j-b]) ) {
1344 mid[i*x_size+j] = 2;
1345 }
1346 }
1347 }
1348 }
1349 }
1350}
1351
1352
1353void susan_corner_draw( uchar *in, CORNER_LIST corner_list,
1354 int x_size, int drawing_mode)
1355{
1356 uchar *p;
1357 int n=0;
1358
1359 _Pragma( "loopbound min 0 max 0" )
1360 while(corner_list[n].info != 7) {
1361 if (drawing_mode==0) {
1362 p = in + (corner_list[n].y-1)*x_size + corner_list[n].x - 1;
1363 *p++=255; *p++=255; *p=255; p+=x_size-2;
1364 *p++=255; *p++=0; *p=255; p+=x_size-2;
1365 *p++=255; *p++=255; *p=255;
1366 n++;
1367 } else {
1368 p = in + corner_list[n].y*x_size + corner_list[n].x;
1369 *p=0;
1370 n++;
1371 }
1372 }
1373}
1374
1375
1376void susan_corners( uchar *in, char *r, uchar *bp,
1377 int max_no, CORNER_LIST corner_list, int x_size, int y_size)
1378{
1379 int n,x,y,sq,xx,yy,
1380 i,j;
1381 float divide;
1382 uchar c,*p,*cp;
1383 char *cgx,*cgy;
1384
1385 susan_wccmemset(r,0,x_size * y_size);
1386
1387 cgx=(char *)susan_wccmalloc(x_size*y_size);
1388 cgy=(char *)susan_wccmalloc(x_size*y_size);
1389
1390 _Pragma( "loopbound min 85 max 85" )
1391 for (i=5;i<y_size-5;i++) {
1392 _Pragma( "loopbound min 66 max 66" )
1393 for (j=5;j<x_size-5;j++) {
1394 n=100;
1395 p=in + (i-3)*x_size + j - 1;
1396 cp=bp + in[i*x_size+j];
1397
1398 n+=*(cp-*p++);
1399 n+=*(cp-*p++);
1400 n+=*(cp-*p);
1401 p+=x_size-3;
1402
1403 n+=*(cp-*p++);
1404 n+=*(cp-*p++);
1405 n+=*(cp-*p++);
1406 n+=*(cp-*p++);
1407 n+=*(cp-*p);
1408 p+=x_size-5;
1409
1410 n+=*(cp-*p++);
1411 n+=*(cp-*p++);
1412 n+=*(cp-*p++);
1413 n+=*(cp-*p++);
1414 n+=*(cp-*p++);
1415 n+=*(cp-*p++);
1416 n+=*(cp-*p);
1417 p+=x_size-6;
1418
1419 n+=*(cp-*p++);
1420 n+=*(cp-*p++);
1421 n+=*(cp-*p);
1422 if (n<max_no){ /* do this test early and often ONLY to save wasted computation */
1423 p+=2;
1424 n+=*(cp-*p++);
1425 if (n<max_no){
1426 n+=*(cp-*p++);
1427 if (n<max_no){
1428 n+=*(cp-*p);
1429 if (n<max_no){
1430 p+=x_size-6;
1431
1432 n+=*(cp-*p++);
1433 if (n<max_no){
1434 n+=*(cp-*p++);
1435 if (n<max_no){
1436 n+=*(cp-*p++);
1437 if (n<max_no){
1438 n+=*(cp-*p++);
1439 if (n<max_no){
1440 n+=*(cp-*p++);
1441 if (n<max_no){
1442 n+=*(cp-*p++);
1443 if (n<max_no){
1444 n+=*(cp-*p);
1445 if (n<max_no){
1446 p+=x_size-5;
1447
1448 n+=*(cp-*p++);
1449 if (n<max_no){
1450 n+=*(cp-*p++);
1451 if (n<max_no){
1452 n+=*(cp-*p++);
1453 if (n<max_no){
1454 n+=*(cp-*p++);
1455 if (n<max_no){
1456 n+=*(cp-*p);
1457 if (n<max_no){
1458 p+=x_size-3;
1459
1460 n+=*(cp-*p++);
1461 if (n<max_no){
1462 n+=*(cp-*p++);
1463 if (n<max_no){
1464 n+=*(cp-*p);
1465
1466 if (n<max_no) {
1467 x=0;y=0;
1468 p=in + (i-3)*x_size + j - 1;
1469
1470 c=*(cp-*p++);x-=c;y-=3*c;
1471 c=*(cp-*p++);y-=3*c;
1472 c=*(cp-*p);x+=c;y-=3*c;
1473 p+=x_size-3;
1474
1475 c=*(cp-*p++);x-=2*c;y-=2*c;
1476 c=*(cp-*p++);x-=c;y-=2*c;
1477 c=*(cp-*p++);y-=2*c;
1478 c=*(cp-*p++);x+=c;y-=2*c;
1479 c=*(cp-*p);x+=2*c;y-=2*c;
1480 p+=x_size-5;
1481
1482 c=*(cp-*p++);x-=3*c;y-=c;
1483 c=*(cp-*p++);x-=2*c;y-=c;
1484 c=*(cp-*p++);x-=c;y-=c;
1485 c=*(cp-*p++);y-=c;
1486 c=*(cp-*p++);x+=c;y-=c;
1487 c=*(cp-*p++);x+=2*c;y-=c;
1488 c=*(cp-*p);x+=3*c;y-=c;
1489 p+=x_size-6;
1490
1491 c=*(cp-*p++);x-=3*c;
1492 c=*(cp-*p++);x-=2*c;
1493 c=*(cp-*p);x-=c;
1494 p+=2;
1495 c=*(cp-*p++);x+=c;
1496 c=*(cp-*p++);x+=2*c;
1497 c=*(cp-*p);x+=3*c;
1498 p+=x_size-6;
1499
1500 c=*(cp-*p++);x-=3*c;y+=c;
1501 c=*(cp-*p++);x-=2*c;y+=c;
1502 c=*(cp-*p++);x-=c;y+=c;
1503 c=*(cp-*p++);y+=c;
1504 c=*(cp-*p++);x+=c;y+=c;
1505 c=*(cp-*p++);x+=2*c;y+=c;
1506 c=*(cp-*p);x+=3*c;y+=c;
1507 p+=x_size-5;
1508
1509 c=*(cp-*p++);x-=2*c;y+=2*c;
1510 c=*(cp-*p++);x-=c;y+=2*c;
1511 c=*(cp-*p++);y+=2*c;
1512 c=*(cp-*p++);x+=c;y+=2*c;
1513 c=*(cp-*p);x+=2*c;y+=2*c;
1514 p+=x_size-3;
1515
1516 c=*(cp-*p++);x-=c;y+=3*c;
1517 c=*(cp-*p++);y+=3*c;
1518 c=*(cp-*p);x+=c;y+=3*c;
1519
1520 xx=x*x;
1521 yy=y*y;
1522 sq=xx+yy;
1523 if ( sq > ((n*n)/2) ) {
1524 if(yy<xx) {
1525 divide=(float)y/(float)abs(x);
1526 sq=abs(x)/x;
1527 sq=*(cp-in[(i+FTOI(divide))*x_size+j+sq]) +
1528 *(cp-in[(i+FTOI(2*divide))*x_size+j+2*sq]) +
1529 *(cp-in[(i+FTOI(3*divide))*x_size+j+3*sq]);
1530 } else {
1531 divide=(float)x/(float)abs(y);
1532 sq=abs(y)/y;
1533 sq=*(cp-in[(i+sq)*x_size+j+FTOI(divide)]) +
1534 *(cp-in[(i+2*sq)*x_size+j+FTOI(2*divide)]) +
1535 *(cp-in[(i+3*sq)*x_size+j+FTOI(3*divide)]);
1536 }
1537
1538 if(sq>290){
1539 r[i*x_size+j] = max_no-n;
1540 cgx[i*x_size+j] = (51*x)/n;
1541 cgy[i*x_size+j] = (51*y)/n;
1542 }
1543 }
1544 }
1545 }}}}}}}}}}}}}}}}}}
1546 }
1547 }
1548
1549 /* to locate the local maxima */
1550 n=0;
1551 _Pragma( "loopbound min 85 max 85" )
1552 for (i=5;i<y_size-5;i++) {
1553 _Pragma( "loopbound min 66 max 66" )
1554 for (j=5;j<x_size-5;j++) {
1555 x = r[i*x_size+j];
1556 if (x>0) {
1557 /* 5x5 mask */
1558 #ifdef FIVE_SUPP
1559 if ((x>r[(i-1)*x_size+j+2]) &&
1560 (x>r[(i )*x_size+j+1]) &&
1561 (x>r[(i )*x_size+j+2]) &&
1562 (x>r[(i+1)*x_size+j-1]) &&
1563 (x>r[(i+1)*x_size+j ]) &&
1564 (x>r[(i+1)*x_size+j+1]) &&
1565 (x>r[(i+1)*x_size+j+2]) &&
1566 (x>r[(i+2)*x_size+j-2]) &&
1567 (x>r[(i+2)*x_size+j-1]) &&
1568 (x>r[(i+2)*x_size+j ]) &&
1569 (x>r[(i+2)*x_size+j+1]) &&
1570 (x>r[(i+2)*x_size+j+2]) &&
1571 (x>=r[(i-2)*x_size+j-2]) &&
1572 (x>=r[(i-2)*x_size+j-1]) &&
1573 (x>=r[(i-2)*x_size+j ]) &&
1574 (x>=r[(i-2)*x_size+j+1]) &&
1575 (x>=r[(i-2)*x_size+j+2]) &&
1576 (x>=r[(i-1)*x_size+j-2]) &&
1577 (x>=r[(i-1)*x_size+j-1]) &&
1578 (x>=r[(i-1)*x_size+j ]) &&
1579 (x>=r[(i-1)*x_size+j+1]) &&
1580 (x>=r[(i )*x_size+j-2]) &&
1581 (x>=r[(i )*x_size+j-1]) &&
1582 (x>=r[(i+1)*x_size+j-2]) )
1583 #endif
1584 #ifdef SEVEN_SUPP
1585 if ((x>r[(i-3)*x_size+j-3]) &&
1586 (x>r[(i-3)*x_size+j-2]) &&
1587 (x>r[(i-3)*x_size+j-1]) &&
1588 (x>r[(i-3)*x_size+j ]) &&
1589 (x>r[(i-3)*x_size+j+1]) &&
1590 (x>r[(i-3)*x_size+j+2]) &&
1591 (x>r[(i-3)*x_size+j+3]) &&
1592
1593 (x>r[(i-2)*x_size+j-3]) &&
1594 (x>r[(i-2)*x_size+j-2]) &&
1595 (x>r[(i-2)*x_size+j-1]) &&
1596 (x>r[(i-2)*x_size+j ]) &&
1597 (x>r[(i-2)*x_size+j+1]) &&
1598 (x>r[(i-2)*x_size+j+2]) &&
1599 (x>r[(i-2)*x_size+j+3]) &&
1600
1601 (x>r[(i-1)*x_size+j-3]) &&
1602 (x>r[(i-1)*x_size+j-2]) &&
1603 (x>r[(i-1)*x_size+j-1]) &&
1604 (x>r[(i-1)*x_size+j ]) &&
1605 (x>r[(i-1)*x_size+j+1]) &&
1606 (x>r[(i-1)*x_size+j+2]) &&
1607 (x>r[(i-1)*x_size+j+3]) &&
1608
1609 (x>r[(i)*x_size+j-3]) &&
1610 (x>r[(i)*x_size+j-2]) &&
1611 (x>r[(i)*x_size+j-1]) &&
1612 (x>=r[(i)*x_size+j+1]) &&
1613 (x>=r[(i)*x_size+j+2]) &&
1614 (x>=r[(i)*x_size+j+3]) &&
1615
1616 (x>=r[(i+1)*x_size+j-3]) &&
1617 (x>=r[(i+1)*x_size+j-2]) &&
1618 (x>=r[(i+1)*x_size+j-1]) &&
1619 (x>=r[(i+1)*x_size+j ]) &&
1620 (x>=r[(i+1)*x_size+j+1]) &&
1621 (x>=r[(i+1)*x_size+j+2]) &&
1622 (x>=r[(i+1)*x_size+j+3]) &&
1623
1624 (x>=r[(i+2)*x_size+j-3]) &&
1625 (x>=r[(i+2)*x_size+j-2]) &&
1626 (x>=r[(i+2)*x_size+j-1]) &&
1627 (x>=r[(i+2)*x_size+j ]) &&
1628 (x>=r[(i+2)*x_size+j+1]) &&
1629 (x>=r[(i+2)*x_size+j+2]) &&
1630 (x>=r[(i+2)*x_size+j+3]) &&
1631
1632 (x>=r[(i+3)*x_size+j-3]) &&
1633 (x>=r[(i+3)*x_size+j-2]) &&
1634 (x>=r[(i+3)*x_size+j-1]) &&
1635 (x>=r[(i+3)*x_size+j ]) &&
1636 (x>=r[(i+3)*x_size+j+1]) &&
1637 (x>=r[(i+3)*x_size+j+2]) &&
1638 (x>=r[(i+3)*x_size+j+3]) )
1639 #endif
1640 {
1641 corner_list[n].info=0;
1642 corner_list[n].x=j;
1643 corner_list[n].y=i;
1644 corner_list[n].dx=cgx[i*x_size+j];
1645 corner_list[n].dy=cgy[i*x_size+j];
1646 corner_list[n].I=in[i*x_size+j];
1647 n++;
1648 if(n==MAX_CORNERS){
1649 /* "Too many corners." */
1650 }
1651 }
1652 }
1653 }
1654 }
1655 corner_list[n].info=7;
1656}
1657
1658
1659void susan_corners_quick( uchar *in, char *r, uchar *bp,
1660 int max_no, CORNER_LIST corner_list, int x_size, int y_size )
1661{
1662 int n,x,y,i,j;
1663 uchar *p,*cp;
1664
1665 susan_wccmemset(r,0,x_size * y_size);
1666
1667 _Pragma( "loopbound min 0 max 0" )
1668 for (i=7;i<y_size-7;i++) {
1669 _Pragma( "loopbound min 0 max 0" )
1670 for (j=7;j<x_size-7;j++) {
1671 n=100;
1672 p=in + (i-3)*x_size + j - 1;
1673 cp=bp + in[i*x_size+j];
1674
1675 n+=*(cp-*p++);
1676 n+=*(cp-*p++);
1677 n+=*(cp-*p);
1678 p+=x_size-3;
1679
1680 n+=*(cp-*p++);
1681 n+=*(cp-*p++);
1682 n+=*(cp-*p++);
1683 n+=*(cp-*p++);
1684 n+=*(cp-*p);
1685 p+=x_size-5;
1686
1687 n+=*(cp-*p++);
1688 n+=*(cp-*p++);
1689 n+=*(cp-*p++);
1690 n+=*(cp-*p++);
1691 n+=*(cp-*p++);
1692 n+=*(cp-*p++);
1693 n+=*(cp-*p);
1694 p+=x_size-6;
1695
1696 n+=*(cp-*p++);
1697 n+=*(cp-*p++);
1698 n+=*(cp-*p);
1699 if (n<max_no){
1700 p+=2;
1701 n+=*(cp-*p++);
1702 if (n<max_no){
1703 n+=*(cp-*p++);
1704 if (n<max_no){
1705 n+=*(cp-*p);
1706 if (n<max_no){
1707 p+=x_size-6;
1708
1709 n+=*(cp-*p++);
1710 if (n<max_no){
1711 n+=*(cp-*p++);
1712 if (n<max_no){
1713 n+=*(cp-*p++);
1714 if (n<max_no){
1715 n+=*(cp-*p++);
1716 if (n<max_no){
1717 n+=*(cp-*p++);
1718 if (n<max_no){
1719 n+=*(cp-*p++);
1720 if (n<max_no){
1721 n+=*(cp-*p);
1722 if (n<max_no){
1723 p+=x_size-5;
1724
1725 n+=*(cp-*p++);
1726 if (n<max_no){
1727 n+=*(cp-*p++);
1728 if (n<max_no){
1729 n+=*(cp-*p++);
1730 if (n<max_no){
1731 n+=*(cp-*p++);
1732 if (n<max_no){
1733 n+=*(cp-*p);
1734 if (n<max_no){
1735 p+=x_size-3;
1736
1737 n+=*(cp-*p++);
1738 if (n<max_no){
1739 n+=*(cp-*p++);
1740 if (n<max_no){
1741 n+=*(cp-*p);
1742
1743 if (n<max_no) {
1744 r[i*x_size+j] = max_no-n;
1745 }
1746 }}}}}}}}}}}}}}}}}}
1747 }
1748 }
1749
1750 /* to locate the local maxima */
1751 n=0;
1752 _Pragma( "loopbound min 0 max 0" )
1753 for (i=7;i<y_size-7;i++) {
1754 _Pragma( "loopbound min 0 max 0" )
1755 for (j=7;j<x_size-7;j++) {
1756 x = r[i*x_size+j];
1757 if (x>0) {
1758 /* 5x5 mask */
1759 #ifdef FIVE_SUPP
1760 if ((x>r[(i-1)*x_size+j+2]) &&
1761 (x>r[(i )*x_size+j+1]) &&
1762 (x>r[(i )*x_size+j+2]) &&
1763 (x>r[(i+1)*x_size+j-1]) &&
1764 (x>r[(i+1)*x_size+j ]) &&
1765 (x>r[(i+1)*x_size+j+1]) &&
1766 (x>r[(i+1)*x_size+j+2]) &&
1767 (x>r[(i+2)*x_size+j-2]) &&
1768 (x>r[(i+2)*x_size+j-1]) &&
1769 (x>r[(i+2)*x_size+j ]) &&
1770 (x>r[(i+2)*x_size+j+1]) &&
1771 (x>r[(i+2)*x_size+j+2]) &&
1772 (x>=r[(i-2)*x_size+j-2]) &&
1773 (x>=r[(i-2)*x_size+j-1]) &&
1774 (x>=r[(i-2)*x_size+j ]) &&
1775 (x>=r[(i-2)*x_size+j+1]) &&
1776 (x>=r[(i-2)*x_size+j+2]) &&
1777 (x>=r[(i-1)*x_size+j-2]) &&
1778 (x>=r[(i-1)*x_size+j-1]) &&
1779 (x>=r[(i-1)*x_size+j ]) &&
1780 (x>=r[(i-1)*x_size+j+1]) &&
1781 (x>=r[(i )*x_size+j-2]) &&
1782 (x>=r[(i )*x_size+j-1]) &&
1783 (x>=r[(i+1)*x_size+j-2]) )
1784 #endif
1785 #ifdef SEVEN_SUPP
1786 if ((x>r[(i-3)*x_size+j-3]) &&
1787 (x>r[(i-3)*x_size+j-2]) &&
1788 (x>r[(i-3)*x_size+j-1]) &&
1789 (x>r[(i-3)*x_size+j ]) &&
1790 (x>r[(i-3)*x_size+j+1]) &&
1791 (x>r[(i-3)*x_size+j+2]) &&
1792 (x>r[(i-3)*x_size+j+3]) &&
1793
1794 (x>r[(i-2)*x_size+j-3]) &&
1795 (x>r[(i-2)*x_size+j-2]) &&
1796 (x>r[(i-2)*x_size+j-1]) &&
1797 (x>r[(i-2)*x_size+j ]) &&
1798 (x>r[(i-2)*x_size+j+1]) &&
1799 (x>r[(i-2)*x_size+j+2]) &&
1800 (x>r[(i-2)*x_size+j+3]) &&
1801
1802 (x>r[(i-1)*x_size+j-3]) &&
1803 (x>r[(i-1)*x_size+j-2]) &&
1804 (x>r[(i-1)*x_size+j-1]) &&
1805 (x>r[(i-1)*x_size+j ]) &&
1806 (x>r[(i-1)*x_size+j+1]) &&
1807 (x>r[(i-1)*x_size+j+2]) &&
1808 (x>r[(i-1)*x_size+j+3]) &&
1809
1810 (x>r[(i)*x_size+j-3]) &&
1811 (x>r[(i)*x_size+j-2]) &&
1812 (x>r[(i)*x_size+j-1]) &&
1813 (x>=r[(i)*x_size+j+1]) &&
1814 (x>=r[(i)*x_size+j+2]) &&
1815 (x>=r[(i)*x_size+j+3]) &&
1816
1817 (x>=r[(i+1)*x_size+j-3]) &&
1818 (x>=r[(i+1)*x_size+j-2]) &&
1819 (x>=r[(i+1)*x_size+j-1]) &&
1820 (x>=r[(i+1)*x_size+j ]) &&
1821 (x>=r[(i+1)*x_size+j+1]) &&
1822 (x>=r[(i+1)*x_size+j+2]) &&
1823 (x>=r[(i+1)*x_size+j+3]) &&
1824
1825 (x>=r[(i+2)*x_size+j-3]) &&
1826 (x>=r[(i+2)*x_size+j-2]) &&
1827 (x>=r[(i+2)*x_size+j-1]) &&
1828 (x>=r[(i+2)*x_size+j ]) &&
1829 (x>=r[(i+2)*x_size+j+1]) &&
1830 (x>=r[(i+2)*x_size+j+2]) &&
1831 (x>=r[(i+2)*x_size+j+3]) &&
1832
1833 (x>=r[(i+3)*x_size+j-3]) &&
1834 (x>=r[(i+3)*x_size+j-2]) &&
1835 (x>=r[(i+3)*x_size+j-1]) &&
1836 (x>=r[(i+3)*x_size+j ]) &&
1837 (x>=r[(i+3)*x_size+j+1]) &&
1838 (x>=r[(i+3)*x_size+j+2]) &&
1839 (x>=r[(i+3)*x_size+j+3]) )
1840 #endif
1841 {
1842 corner_list[n].info=0;
1843 corner_list[n].x=j;
1844 corner_list[n].y=i;
1845 x = in[(i-2)*x_size+j-2] + in[(i-2)*x_size+j-1] + in[(i-2)*x_size+j] + in[(i-2)*x_size+j+1] + in[(i-2)*x_size+j+2] +
1846 in[(i-1)*x_size+j-2] + in[(i-1)*x_size+j-1] + in[(i-1)*x_size+j] + in[(i-1)*x_size+j+1] + in[(i-1)*x_size+j+2] +
1847 in[(i )*x_size+j-2] + in[(i )*x_size+j-1] + in[(i )*x_size+j] + in[(i )*x_size+j+1] + in[(i )*x_size+j+2] +
1848 in[(i+1)*x_size+j-2] + in[(i+1)*x_size+j-1] + in[(i+1)*x_size+j] + in[(i+1)*x_size+j+1] + in[(i+1)*x_size+j+2] +
1849 in[(i+2)*x_size+j-2] + in[(i+2)*x_size+j-1] + in[(i+2)*x_size+j] + in[(i+2)*x_size+j+1] + in[(i+2)*x_size+j+2];
1850
1851 corner_list[n].I=x/25;
1852 /*corner_list[n].I=in[i*x_size+j];*/
1853 x = in[(i-2)*x_size+j+2] + in[(i-1)*x_size+j+2] + in[(i)*x_size+j+2] + in[(i+1)*x_size+j+2] + in[(i+2)*x_size+j+2] -
1854 (in[(i-2)*x_size+j-2] + in[(i-1)*x_size+j-2] + in[(i)*x_size+j-2] + in[(i+1)*x_size+j-2] + in[(i+2)*x_size+j-2]);
1855 x += x + in[(i-2)*x_size+j+1] + in[(i-1)*x_size+j+1] + in[(i)*x_size+j+1] + in[(i+1)*x_size+j+1] + in[(i+2)*x_size+j+1] -
1856 (in[(i-2)*x_size+j-1] + in[(i-1)*x_size+j-1] + in[(i)*x_size+j-1] + in[(i+1)*x_size+j-1] + in[(i+2)*x_size+j-1]);
1857
1858 y = in[(i+2)*x_size+j-2] + in[(i+2)*x_size+j-1] + in[(i+2)*x_size+j] + in[(i+2)*x_size+j+1] + in[(i+2)*x_size+j+2] -
1859 (in[(i-2)*x_size+j-2] + in[(i-2)*x_size+j-1] + in[(i-2)*x_size+j] + in[(i-2)*x_size+j+1] + in[(i-2)*x_size+j+2]);
1860 y += y + in[(i+1)*x_size+j-2] + in[(i+1)*x_size+j-1] + in[(i+1)*x_size+j] + in[(i+1)*x_size+j+1] + in[(i+1)*x_size+j+2] -
1861 (in[(i-1)*x_size+j-2] + in[(i-1)*x_size+j-1] + in[(i-1)*x_size+j] + in[(i-1)*x_size+j+1] + in[(i-1)*x_size+j+2]);
1862 corner_list[n].dx=x/15;
1863 corner_list[n].dy=y/15;
1864 n++;
1865 if(n==MAX_CORNERS){
1866 /* "Too many corners." */
1867 }
1868 }
1869 }
1870 }
1871 }
1872 corner_list[n].info=7;
1873}
1874
1875
1876void susan_call_susan( struct wccFILE *inputFile, int mode )
1877{
1878 uchar *in, *bp, *mid;
1879 int x_size, y_size;
1880 char *r;
1881 CORNER_LIST corner_list;
1882
1883 susan_get_image( inputFile, &in, &x_size, &y_size );
1884
1885 if (susan_dt<0) susan_three_by_three=1;
1886 if ( (susan_principle_conf==1) && (mode==0) )
1887 mode=1;
1888
1889 switch (mode) {
1890 case 0:
1891 /* {{{ smoothing */
1892
1893 susan_setup_brightness_lut(&bp,susan_bt,2);
1894 susan_smoothing(susan_three_by_three,in,susan_dt,x_size,y_size,bp);
1895
1896 break;
1897 case 1:
1898 /* {{{ edges */
1899
1900 r = (char *) susan_wccmalloc(x_size * y_size);
1901 susan_setup_brightness_lut(&bp,susan_bt,6);
1902
1903 if (susan_principle_conf) {
1904 if (susan_three_by_three)
1905 susan_principle_small(in,r,bp,susan_max_no_edges,x_size,y_size);
1906 else
1907 susan_principle(in,r,bp,susan_max_no_edges,x_size,y_size);
1908 susan_int_to_uchar(r,in,x_size*y_size);
1909 } else {
1910 mid = (uchar *)susan_wccmalloc(x_size*y_size);
1911 susan_wccmemset(mid,100,x_size * y_size); /* note not set to zero */
1912
1913 if (susan_three_by_three)
1914 susan_edges_small(in,r,mid,bp,susan_max_no_edges,x_size,y_size);
1915 else
1916 susan_edges(in,r,mid,bp,susan_max_no_edges,x_size,y_size);
1917 if(susan_thin_post_proc)
1918 susan_thin(r,mid,x_size,y_size);
1919 susan_edge_draw(in,mid,x_size,y_size,susan_drawing_mode);
1920 }
1921
1922 break;
1923 case 2:
1924 /* {{{ corners */
1925
1926 r = (char *) susan_wccmalloc(x_size * y_size);
1927 susan_setup_brightness_lut(&bp,susan_bt,6);
1928
1929 if (susan_principle_conf) {
1930 susan_principle(in,r,bp,susan_max_no_corners,x_size,y_size);
1931 susan_int_to_uchar(r,in,x_size*y_size);
1932 } else {
1933 if(susan_susan_quick)
1934 susan_corners_quick(in,r,bp,susan_max_no_corners,corner_list,x_size,y_size);
1935 else
1936 susan_corners(in,r,bp,susan_max_no_corners,corner_list,x_size,y_size);
1937 susan_corner_draw(in,corner_list,x_size,susan_drawing_mode);
1938 }
1939
1940 break;
1941 }
1942
1943 susan_put_image(x_size,y_size);
1944}
1945
1946void susan_init( void )
1947{
1948 volatile int a = 0;
1949 susan_file.data = susan_input;
1950 susan_file.size = 7292;
1951 susan_file.size += a;
1952 susan_file.cur_pos = 0;
1953 susan_file.cur_pos += a;
1954
1955 susan_dt = 4.0;
1956 susan_dt += a;
1957 susan_bt = 20;
1958 susan_bt += a;
1959 susan_principle_conf = 0;
1960 susan_principle_conf += a;
1961 susan_thin_post_proc = 1;
1962 susan_thin_post_proc += a;
1963 susan_three_by_three = 0;
1964 susan_three_by_three += a;
1965 susan_drawing_mode = 0;
1966 susan_drawing_mode += a;
1967 susan_susan_quick = 0;
1968 susan_susan_quick += a;
1969 susan_max_no_corners = 50;
1970 susan_max_no_corners += a;
1971 susan_max_no_edges = 50;
1972 susan_max_no_edges += a;
1973
1974 // mode=0; /* Smoothing mode */
1975 // mode=1; /* Edges mode */
1976 // mode=2; /* Corners mode */
1977
1978 // principle=1; /* Output initial enhancement image only; corners or edges mode (default is edges mode) */
1979 // thin_post_proc=0; /* No post-processing on the binary edge map (runs much faster); edges mode */
1980 // drawing_mode=1; /* Mark corners/edges with single black points instead of black with white border; corners or edges mode */
1981 // three_by_three=1; /* Use flat 3x3 mask, edges or smoothing mode */
1982 // susan_quick=1; /* Use faster (and usually stabler) corner mode; edge-like corner suppression not carried out; corners mode */
1983 // dt=10.0; /* Distance threshold, smoothing mode, (default=4) */
1984 // bt=50; /* Brightness threshold, all modes, (default=20) */
1985}
1986
1987void susan_main( void )
1988{
1989 susan_call_susan( &susan_file, 0 );
1990 susan_wccfreeall();
1991 susan_call_susan( &susan_file, 1 );
1992 susan_wccfreeall();
1993 susan_call_susan( &susan_file, 2 );
1994 susan_wccfreeall();
1995}
1996
1997int susan_return( void )
1998{
1999 return 0;
2000}
2001
2002int main( int argc, char **argv )
2003{
2004 SET_UP
2005 for (jobsComplete=-1; jobsComplete<maxJobs; jobsComplete++){
2006 START_LOOP
2007 susan_init();
2008 susan_main();
2009 STOP_LOOP
2010 }
2011 WRITE_TO_FILE
2012
2013 return susan_return();
2014}
generated by cgit v1.2.2 (git 2.25.0) at 2025-08-30 22:23:19 -0400