source: trunk/Softwares/MiBench/src/c/susan-susan.c @ 141

Last change on this file since 141 was 117, checked in by rosiere, 15 years ago

1) Platforms : add new organization for test
2) Load_Store_Unit : add array to count nb_check in store_queue
3) Issue_queue and Core_Glue : rewrite the issue network
4) Special_Register_Unit : add reset value to register CID
5) Softwares : add multicontext test
6) Softwares : add SPECINT
7) Softwares : add MiBench?
7) Read_queue : inhib access for r0
8) Change Core_Glue (network) - dont yet support priority and load balancing scheme

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