Tekkotsu Homepage
Demos
Overview
Downloads
Dev. Resources
Reference
Credits

susan.cc

Go to the documentation of this file.
00001 //-*-c++-*-
00002 
00003 /* Functions taken from the SUSAN edge detector package by Dr. Stephen Smith.
00004    The original file can be found at susan2l.c
00005    Ported to Tekkotsu/C++ by Matt Carson and Dave Touretzky, and
00006    distributed here by the kind permission of Dr. Smith */
00007 
00008 #include <string.h>
00009 #include <stdlib.h>
00010 #include <cmath>
00011 
00012 #include "susan.h"
00013 
00014 /* {{{ Copyright etc. */
00015 
00016 /**********************************************************************\
00017 
00018   SUSAN Version 2l by Stephen Smith
00019   Oxford Centre for Functional Magnetic Resonance Imaging of the Brain,
00020   Department of Clinical Neurology, Oxford University, Oxford, UK
00021   (Previously in Computer Vision and Image Processing Group - now
00022   Computer Vision and Electro Optics Group - DERA Chertsey, UK)
00023   Email:    steve@fmrib.ox.ac.uk
00024   WWW:      http://www.fmrib.ox.ac.uk/~steve
00025 
00026   (C) Crown Copyright (1995-1999), Defence Evaluation and Research Agency,
00027   Farnborough, Hampshire, GU14 6TD, UK
00028   DERA WWW site:
00029   http://www.dera.gov.uk/
00030   DERA Computer Vision and Electro Optics Group WWW site:
00031   http://www.dera.gov.uk/imageprocessing/dera/group_home.html
00032   DERA Computer Vision and Electro Optics Group point of contact:
00033   Dr. John Savage, jtsavage@dera.gov.uk, +44 1344 633203
00034 
00035   A UK patent has been granted: "Method for digitally processing
00036   images to determine the position of edges and/or corners therein for
00037   guidance of unmanned vehicle", UK Patent 2272285. Proprietor:
00038   Secretary of State for Defence, UK. 15 January 1997
00039 
00040   This code is issued for research purposes only and remains the
00041   property of the UK Secretary of State for Defence. This code must
00042   not be passed on without this header information being kept
00043   intact. This code must not be sold.
00044 
00045 \**********************************************************************/
00046 
00047 /* }}} */
00048 /* {{{ Readme First */
00049 
00050 /**********************************************************************\
00051 
00052   SUSAN Version 2l
00053   SUSAN = Smallest Univalue Segment Assimilating Nucleus
00054 
00055   Email:    steve@fmrib.ox.ac.uk
00056   WWW:      http://www.fmrib.ox.ac.uk/~steve
00057 
00058   Related paper:
00059   @article{Smith97,
00060         author = "Smith, S.M. and Brady, J.M.",
00061         title = "{SUSAN} - A New Approach to Low Level Image Processing",
00062         journal = "Int. Journal of Computer Vision",
00063         pages = "45--78",
00064         volume = "23",
00065         number = "1",
00066         month = "May",
00067         year = 1997}
00068 
00069   To be registered for automatic (bug) updates of SUSAN, send an email.
00070 
00071   Compile with:
00072   gcc -O4 -o susan susan2l.c -lm
00073 
00074   See following section for different machine information. Please
00075   report any bugs (and fixes). There are a few optional changes that
00076   can be made in the "defines" section which follows shortly.
00077 
00078   Usage: type "susan" to get usage. Only PGM format files can be input
00079   and output. Utilities such as the netpbm package and XV can be used
00080   to convert to and from other formats. Any size of image can be
00081   processed.
00082 
00083   This code is written using an emacs folding mode, making moving
00084   around the different sections very easy. This is why there are
00085   various marks within comments and why comments are indented.
00086 
00087 
00088   SUSAN QUICK:
00089 
00090   This version of the SUSAN corner finder does not do all the
00091   false-corner suppression and thus is faster and produced some false
00092   positives, particularly on strong edges. However, because there are
00093   less stages involving thresholds etc., the corners that are
00094   correctly reported are usually more stable than those reported with
00095   the full algorithm. Thus I recommend at least TRYING this algorithm
00096   for applications where stability is important, e.g., tracking.
00097 
00098   THRESHOLDS:
00099 
00100   There are two thresholds which can be set at run-time. These are the
00101   brightness threshold (t) and the distance threshold (d).
00102 
00103   SPATIAL CONTROL: d
00104 
00105   In SUSAN smoothing d controls the size of the Gaussian mask; its
00106   default is 4.0. Increasing d gives more smoothing. In edge finding,
00107   a fixed flat mask is used, either 37 pixels arranged in a "circle"
00108   (default), or a 3 by 3 mask which gives finer detail. In corner
00109   finding, only the larger 37 pixel mask is used; d is not
00110   variable. In smoothing, the flat 3 by 3 mask can be used instead of
00111   a larger Gaussian mask; this gives low smoothing and fast operation.
00112 
00113   BRIGHTNESS CONTROL: t
00114 
00115   In all three algorithms, t can be varied (default=20); this is the
00116   main threshold to be varied. It determines the maximum difference in
00117   greylevels between two pixels which allows them to be considered
00118   part of the same "region" in the image. Thus it can be reduced to
00119   give more edges or corners, i.e. to be more sensitive, and vice
00120   versa. In smoothing, reducing t gives less smoothing, and vice
00121   versa. Set t=10 for the test image available from the SUSAN web
00122   page.
00123 
00124   ITERATIONS:
00125 
00126   With SUSAN smoothing, more smoothing can also be obtained by
00127   iterating the algorithm several times. This has a different effect
00128   from varying d or t.
00129 
00130   FIXED MASKS:
00131 
00132   37 pixel mask:    ooo       3 by 3 mask:  ooo
00133                    ooooo                    ooo
00134                   ooooooo                   ooo
00135                   ooooooo
00136                   ooooooo
00137                    ooooo
00138                     ooo
00139 
00140   CORNER ATTRIBUTES dx, dy and I
00141   (Only read this if you are interested in the C implementation or in
00142   using corner attributes, e.g., for corner matching)
00143 
00144   Corners reported in the corner list have attributes associated with
00145   them as well as positions. This is useful, for example, when
00146   attempting to match corners from one image to another, as these
00147   attributes can often be fairly unchanged between images. The
00148   attributes are dx, dy and I. I is the value of image brightness at
00149   the position of the corner. In the case of susan_corners_quick, dx
00150   and dy are the first order derivatives (differentials) of the image
00151   brightness in the x and y directions respectively, at the position
00152   of the corner. In the case of normal susan corner finding, dx and dy
00153   are scaled versions of the position of the centre of gravity of the
00154   USAN with respect to the centre pixel (nucleus).
00155 
00156   BRIGHTNESS FUNCTION LUT IMPLEMENTATION:
00157   (Only read this if you are interested in the C implementation)
00158 
00159   The SUSAN brightness function is implemented as a LUT
00160   (Look-Up-Table) for speed. The resulting pointer-based code is a
00161   little hard to follow, so here is a brief explanation. In
00162   setup_brightness_lut() the LUT is setup. This mallocs enough space
00163   for *bp and then repositions the pointer to the centre of the
00164   malloced space. The SUSAN function e^-(x^6) or e^-(x^2) is
00165   calculated and converted to a uchar in the range 0-100, for all
00166   possible image brightness differences (including negative
00167   ones). Thus bp[23] is the output for a brightness difference of 23
00168   greylevels. In the SUSAN algorithms this LUT is used as follows:
00169 
00170   p=in + (i-3)*x_size + j - 1;
00171   p points to the first image pixel in the circular mask surrounding
00172   point (x,y).
00173 
00174   cp=bp + in[i*x_size+j];
00175   cp points to a position in the LUT corresponding to the brightness
00176   of the centre pixel (x,y).
00177 
00178   now for every pixel within the mask surrounding (x,y),
00179   n+=*(cp-*p++);
00180   the brightness difference function is found by moving the cp pointer
00181   down by an amount equal to the value of the pixel pointed to by p,
00182   thus subtracting the two brightness values and performing the
00183   exponential function. This value is added to n, the running USAN
00184   area.
00185 
00186   in SUSAN smoothing, the variable height mask is implemented by
00187   multiplying the above by the moving mask pointer, reset for each new
00188   centre pixel.
00189   tmp = *dpt++ * *(cp-brightness);
00190 
00191 \**********************************************************************/
00192 
00193 /* }}} */
00194 /* {{{ Machine Information */
00195 
00196 /**********************************************************************\
00197 
00198   Success has been reported with the following:
00199 
00200   MACHINE  OS         COMPILER
00201 
00202   Sun      4.1.4      bundled C, gcc
00203 
00204   Next
00205 
00206   SGI      IRIX       SGI cc
00207 
00208   DEC      Unix V3.2+ 
00209 
00210   IBM RISC AIX        gcc
00211 
00212   PC                  Borland 5.0
00213 
00214   PC       Linux      gcc-2.6.3
00215 
00216   PC       Win32      Visual C++ 4.0 (Console Application)
00217 
00218   PC       Win95      Visual C++ 5.0 (Console Application)
00219                       Thanks to Niu Yongsheng <niuysbit@163.net>:
00220                       Use the FOPENB option below
00221 
00222   PC       DOS        djgpp gnu C
00223                       Thanks to Mark Pettovello <mpettove@umdsun2.umd.umich.edu>:
00224                       Use the FOPENB option below
00225 
00226   HP       HP-UX      bundled cc
00227                       Thanks to Brian Dixon <briand@hpcvsgen.cv.hp.com>:
00228                       in ksh:
00229                       export CCOPTS="-Aa -D_HPUX_SOURCE | -lM"
00230                       cc -O3 -o susan susan2l.c
00231 
00232 \**********************************************************************/
00233 
00234 /* }}} */
00235 /* {{{ History */
00236 
00237 /**********************************************************************\
00238 
00239   SUSAN Version 2l, 12/2/99
00240   Changed GNUDOS option to FOPENB.
00241   (Thanks to Niu Yongsheng <niuysbit@163.net>.)
00242   Took out redundant "sq=sq/2;".
00243 
00244   SUSAN Version 2k, 19/8/98:
00245   In corner finding:
00246   Changed if(yy<sq) {...} else if(xx<sq) {...} to
00247           if(yy<xx) {...} else {...}
00248   (Thanks to adq@cim.mcgill.edu - Alain Domercq.)
00249 
00250   SUSAN Version 2j, 22/10/97:
00251   Fixed (mask_size>x_size) etc. tests in smoothing.
00252   Added a couple of free() calls for cgx and cgy.
00253   (Thanks to geoffb@ucs.ed.ac.uk - Geoff Browitt.)
00254 
00255   SUSAN Version 2i, 21/7/97:
00256   Added information about corner attributes.
00257 
00258   SUSAN Version 2h, 16/12/96:
00259   Added principle (initial enhancement) option.
00260 
00261   SUSAN Version 2g, 2/7/96:
00262   Minor superficial changes to code.
00263 
00264   SUSAN Version 2f, 16/1/96:
00265   Added GNUDOS option (now called FOPENB; see options below).
00266 
00267   SUSAN Version 2e, 9/1/96:
00268   Added -b option.
00269   Fixed 1 pixel horizontal offset error for drawing edges.
00270 
00271   SUSAN Version 2d, 27/11/95:
00272   Fixed loading of certain PGM files in get_image (again!)
00273 
00274   SUSAN Version 2c, 22/11/95:
00275   Fixed loading of certain PGM files in get_image.
00276   (Thanks to qu@San-Jose.ate.slb.com - Gongyuan Qu.)
00277 
00278   SUSAN Version 2b, 9/11/95:
00279   removed "z==" error in edges routines.
00280 
00281   SUSAN Version 2a, 6/11/95:
00282   Removed a few unnecessary variable declarations.
00283   Added different machine information.
00284   Changed "header" in get_image to char.
00285 
00286   SUSAN Version 2, 1/11/95: first combined version able to take any
00287   image sizes.
00288 
00289   SUSAN "Versions 1", circa 1992: the various SUSAN algorithms were
00290   developed during my doctorate within different programs and for
00291   fixed image sizes. The algorithms themselves are virtually unaltered
00292   between "versions 1" and the combined program, version 2.
00293 
00294 \**********************************************************************/
00295 
00296 namespace DualCoding {
00297 
00298 /* }}} */
00299 /* {{{ setup_brightness_lut(bp,thresh,form) */
00300 
00301 /*void setup_brightness_lut(bp,thresh,form)
00302   uchar **bp;
00303   int   thresh, form;*/
00304 void setup_brightness_lut(uchar **bp, int thresh, int form)
00305 {
00306 int   k;
00307 float temp;
00308 
00309   *bp=(unsigned char *)malloc(516);
00310   *bp=*bp+258;
00311 
00312   for(k=-256;k<257;k++)
00313   {
00314     temp=((float)k)/((float)thresh);
00315     temp=temp*temp;
00316     if (form==6)
00317       temp=temp*temp*temp;
00318     temp=100*std::exp(-temp);
00319     *(*bp+k)= (uchar)temp;
00320   }
00321 }
00322 
00323 
00324 /* }}} */
00325 /* {{{ susan principle */
00326 
00327 /* {{{ susan_principle(in,r,bp,max_no,x_size,y_size) */
00328 
00329 /*void susan_principle(in,r,bp,max_no,x_size,y_size)
00330   uchar *in, *bp;
00331   int   *r, max_no, x_size, y_size;*/
00332 void susan_principle(uchar *in, int *r, uchar *bp, int max_no,
00333          int x_size, int y_size)
00334 {
00335 int   i, j, n;
00336 uchar *p,*cp;
00337 
00338   memset (r,0,x_size * y_size * sizeof(int));
00339 
00340 
00341   for (i=3;i<y_size-3;i++)
00342     for (j=3;j<x_size-3;j++)
00343     {
00344       n=100;
00345       p=in + (i-3)*x_size + j - 1;
00346       cp=bp + in[i*x_size+j];
00347 
00348       n+=*(cp-*p++);
00349       n+=*(cp-*p++);
00350       n+=*(cp-*p);
00351       p+=x_size-3; 
00352 
00353       n+=*(cp-*p++);
00354       n+=*(cp-*p++);
00355       n+=*(cp-*p++);
00356       n+=*(cp-*p++);
00357       n+=*(cp-*p);
00358       p+=x_size-5;
00359 
00360       n+=*(cp-*p++);
00361       n+=*(cp-*p++);
00362       n+=*(cp-*p++);
00363       n+=*(cp-*p++);
00364       n+=*(cp-*p++);
00365       n+=*(cp-*p++);
00366       n+=*(cp-*p);
00367       p+=x_size-6;
00368 
00369       n+=*(cp-*p++);
00370       n+=*(cp-*p++);
00371       n+=*(cp-*p);
00372       p+=2;
00373       n+=*(cp-*p++);
00374       n+=*(cp-*p++);
00375       n+=*(cp-*p);
00376       p+=x_size-6;
00377 
00378       n+=*(cp-*p++);
00379       n+=*(cp-*p++);
00380       n+=*(cp-*p++);
00381       n+=*(cp-*p++);
00382       n+=*(cp-*p++);
00383       n+=*(cp-*p++);
00384       n+=*(cp-*p);
00385       p+=x_size-5;
00386 
00387       n+=*(cp-*p++);
00388       n+=*(cp-*p++);
00389       n+=*(cp-*p++);
00390       n+=*(cp-*p++);
00391       n+=*(cp-*p);
00392       p+=x_size-3;
00393 
00394       n+=*(cp-*p++);
00395       n+=*(cp-*p++);
00396       n+=*(cp-*p);
00397 
00398       if (n<=max_no)
00399         r[i*x_size+j] = max_no - n;
00400     }
00401 }
00402 
00403 /* }}} */
00404 /* {{{ susan_thin(r,mid,x_size,y_size) */
00405 
00406 /* only one pass is needed as i,j are decremented if necessary to go
00407    back and do bits again */
00408 
00409 /*susan_thin(r,mid,x_size,y_size)
00410   uchar *mid;
00411   int   *r, x_size, y_size;*/
00412 void susan_thin(int *r, uchar *mid, int x_size, int y_size)
00413 {
00414   int   l[9], centre, //nlinks, npieces,
00415       b01, b12, b21, b10,
00416       p1, p2, p3, p4,
00417       b00, b02, b20, b22,
00418       m, n, a=0, b=0, x, y, i, j;
00419 uchar *mp;
00420 
00421   for (i=4;i<y_size-4;i++)
00422     for (j=4;j<x_size-4;j++)
00423       if (mid[i*x_size+j]<8)
00424       {
00425         centre = r[i*x_size+j];
00426         /* {{{ count number of neighbours */
00427 
00428         mp=mid + (i-1)*x_size + j-1;
00429 
00430         n = (*mp<8) +
00431             (*(mp+1)<8) +
00432             (*(mp+2)<8) +
00433             (*(mp+x_size)<8) +
00434             (*(mp+x_size+2)<8) +
00435             (*(mp+x_size+x_size)<8) +
00436             (*(mp+x_size+x_size+1)<8) +
00437             (*(mp+x_size+x_size+2)<8);
00438 
00439 /* }}} */
00440         /* {{{ n==0 no neighbours - remove point */
00441 
00442         if (n==0)
00443           mid[i*x_size+j]=100;
00444 
00445 /* }}} */
00446         /* {{{ n==1 - extend line if I can */
00447 
00448         /* extension is only allowed a few times - the value of mid is used to control this */
00449 
00450         if ( (n==1) && (mid[i*x_size+j]<6) )
00451         {
00452           /* find maximum neighbour weighted in direction opposite the
00453              neighbour already present. e.g.
00454              have: O O O  weight r by 0 2 3
00455                    X X O              0 0 4
00456                    O O O              0 2 3     */
00457 
00458           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];
00459           l[3]=r[(i  )*x_size+j-1]; l[4]=0;                 l[5]=r[(i  )*x_size+j+1];
00460           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];
00461 
00462           if (mid[(i-1)*x_size+j-1]<8)        { l[0]=0; l[1]=0; l[3]=0; l[2]*=2; 
00463                                                 l[6]*=2; l[5]*=3; l[7]*=3; l[8]*=4; }
00464           else { if (mid[(i-1)*x_size+j]<8)   { l[1]=0; l[0]=0; l[2]=0; l[3]*=2; 
00465                                                 l[5]*=2; l[6]*=3; l[8]*=3; l[7]*=4; }
00466           else { if (mid[(i-1)*x_size+j+1]<8) { l[2]=0; l[1]=0; l[5]=0; l[0]*=2; 
00467                                                 l[8]*=2; l[3]*=3; l[7]*=3; l[6]*=4; }
00468           else { if (mid[(i)*x_size+j-1]<8)   { l[3]=0; l[0]=0; l[6]=0; l[1]*=2; 
00469                                                 l[7]*=2; l[2]*=3; l[8]*=3; l[5]*=4; }
00470           else { if (mid[(i)*x_size+j+1]<8)   { l[5]=0; l[2]=0; l[8]=0; l[1]*=2; 
00471                                                 l[7]*=2; l[0]*=3; l[6]*=3; l[3]*=4; }
00472           else { if (mid[(i+1)*x_size+j-1]<8) { l[6]=0; l[3]=0; l[7]=0; l[0]*=2; 
00473                                                 l[8]*=2; l[1]*=3; l[5]*=3; l[2]*=4; }
00474           else { if (mid[(i+1)*x_size+j]<8)   { l[7]=0; l[6]=0; l[8]=0; l[3]*=2; 
00475                                                 l[5]*=2; l[0]*=3; l[2]*=3; l[1]*=4; }
00476           else { if (mid[(i+1)*x_size+j+1]<8) { l[8]=0; l[5]=0; l[7]=0; l[6]*=2; 
00477                                                 l[2]*=2; l[1]*=3; l[3]*=3; l[0]*=4; } }}}}}}}
00478 
00479           m=0;     /* find the highest point */
00480           for(y=0; y<3; y++)
00481             for(x=0; x<3; x++)
00482               if (l[y+y+y+x]>m) { m=l[y+y+y+x]; a=y; b=x; }
00483 
00484           if (m>0)
00485           {
00486             if (mid[i*x_size+j]<4)
00487               mid[(i+a-1)*x_size+j+b-1] = 4;
00488             else
00489               mid[(i+a-1)*x_size+j+b-1] = mid[i*x_size+j]+1;
00490             if ( (a+a+b) < 3 ) /* need to jump back in image */
00491       {
00492               i+=a-1;
00493               j+=b-2;
00494               if (i<4) i=4;
00495               if (j<4) j=4;
00496       }
00497     }
00498         }
00499 
00500 /* }}} */
00501         /* {{{ n==2 */
00502 
00503         if (n==2)
00504   {
00505           /* put in a bit here to straighten edges */
00506           b00 = mid[(i-1)*x_size+j-1]<8; /* corners of 3x3 */
00507           b02 = mid[(i-1)*x_size+j+1]<8;
00508     b20 = mid[(i+1)*x_size+j-1]<8;
00509           b22 = mid[(i+1)*x_size+j+1]<8;
00510           if ( ((b00+b02+b20+b22)==2) && ((b00|b22)&(b02|b20)))
00511     {  /* case: move a point back into line.
00512                 e.g. X O X  CAN  become X X X
00513                      O X O              O O O
00514                      O O O              O O O    */
00515             if (b00) 
00516       {
00517               if (b02) { x=0; y=-1; }
00518               else     { x=-1; y=0; }
00519       }
00520             else
00521       {
00522               if (b02) { x=1; y=0; }
00523               else     { x=0; y=1; }
00524       }
00525             if (((float)r[(i+y)*x_size+j+x]/(float)centre) > 0.7)
00526       {
00527               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) ) ||
00528                    ( (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) ) )
00529         {
00530                 mid[(i)*x_size+j]=100;
00531                 mid[(i+y)*x_size+j+x]=3;  /* no jumping needed */
00532         }
00533       }
00534     }
00535           else
00536           {
00537             b01 = mid[(i-1)*x_size+j  ]<8;
00538             b12 = mid[(i  )*x_size+j+1]<8;
00539             b21 = mid[(i+1)*x_size+j  ]<8;
00540             b10 = mid[(i  )*x_size+j-1]<8;
00541             /* {{{ right angle ends - not currently used */
00542 
00543 #ifdef IGNORETHIS
00544             if ( (b00&b01)|(b00&b10)|(b02&b01)|(b02&b12)|(b20&b10)|(b20&b21)|(b22&b21)|(b22&b12) )
00545       { /* case; right angle ends. clean up.
00546                  e.g.; X X O  CAN  become X X O
00547                        O X O              O O O
00548                        O O O              O O O        */
00549               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)&
00550                                     ((b00&((2*r[(i-1)*x_size+j+1])>centre))|(b02&((2*r[(i-1)*x_size+j-1])>centre)))) |
00551                    ((b10)&(mid[(i-1)*x_size+j-2]>7)&(mid[(i)*x_size+j-2]>7)&(mid[(i+1)*x_size+j-2]>7)&
00552                                     ((b00&((2*r[(i+1)*x_size+j-1])>centre))|(b20&((2*r[(i-1)*x_size+j-1])>centre)))) |
00553                    ((b12)&(mid[(i-1)*x_size+j+2]>7)&(mid[(i)*x_size+j+2]>7)&(mid[(i+1)*x_size+j+2]>7)&
00554                                     ((b02&((2*r[(i+1)*x_size+j+1])>centre))|(b22&((2*r[(i-1)*x_size+j+1])>centre)))) |
00555                    ((b21)&(mid[(i+2)*x_size+j-1]>7)&(mid[(i+2)*x_size+j]>7)&(mid[(i+2)*x_size+j+1]>7)&
00556                                     ((b20&((2*r[(i+1)*x_size+j+1])>centre))|(b22&((2*r[(i+1)*x_size+j-1])>centre)))) )
00557         {
00558                 mid[(i)*x_size+j]=100;
00559                 if (b10&b20) j-=2;
00560                 if (b00|b01|b02) { i--; j-=2; }
00561           }
00562       }
00563 #endif
00564 
00565 /* }}} */
00566             if ( ((b01+b12+b21+b10)==2) && ((b10|b12)&(b01|b21)) &&
00567                  ((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)))|
00568                 (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)))) )
00569       { /* case; clears odd right angles.
00570                  e.g.; O O O  becomes O O O
00571                        X X O          X O O
00572                        O X O          O X O     */
00573               mid[(i)*x_size+j]=100;
00574               i--;               /* jump back */
00575               j-=2;
00576               if (i<4) i=4;
00577               if (j<4) j=4;
00578       }
00579     }
00580   }
00581 
00582 /* }}} */
00583         /* {{{ n>2 the thinning is done here without breaking connectivity */
00584 
00585         if (n>2)
00586         {
00587           b01 = mid[(i-1)*x_size+j  ]<8;
00588           b12 = mid[(i  )*x_size+j+1]<8;
00589           b21 = mid[(i+1)*x_size+j  ]<8;
00590           b10 = mid[(i  )*x_size+j-1]<8;
00591           if((b01+b12+b21+b10)>1)
00592           {
00593             b00 = mid[(i-1)*x_size+j-1]<8;
00594             b02 = mid[(i-1)*x_size+j+1]<8;
00595       b20 = mid[(i+1)*x_size+j-1]<8;
00596       b22 = mid[(i+1)*x_size+j+1]<8;
00597             p1 = b00 | b01;
00598             p2 = b02 | b12;
00599             p3 = b22 | b21;
00600             p4 = b20 | b10;
00601 
00602             if( ((p1 + p2 + p3 + p4) - ((b01 & p2)+(b12 & p3)+(b21 & p4)+(b10 & p1))) < 2)
00603             {
00604               mid[(i)*x_size+j]=100;
00605               i--;
00606               j-=2;
00607               if (i<4) i=4;
00608               if (j<4) j=4;
00609             }
00610           }
00611         }
00612 /* }}} */
00613       }
00614 }
00615 
00616 
00617 
00618 /* }}} */
00619 /* {{{ susan_edges(in,r,sf,max_no,out) */
00620 
00621 /*susan_edges(in,r,mid,bp,max_no,x_size,y_size)
00622   uchar *in, *bp, *mid;
00623   int   *r, max_no, x_size, y_size;*/
00624 void susan_edges_internal(uchar *in, int *r, uchar *mid, uchar *bp, 
00625       int max_no, int x_size, int y_size)
00626 {
00627 float z;
00628 int   do_symmetry, i, j, m, n, a, b, x, y, w;
00629 uchar c,*p,*cp;
00630 
00631   memset (r,0,x_size * y_size * sizeof(int));
00632 
00633   for (i=3;i<y_size-3;i++)
00634     for (j=3;j<x_size-3;j++)
00635     {
00636       n=100;
00637       p=in + (i-3)*x_size + j - 1;
00638       cp=bp + in[i*x_size+j];
00639 
00640       n+=*(cp-*p++);
00641       n+=*(cp-*p++);
00642       n+=*(cp-*p);
00643       p+=x_size-3; 
00644 
00645       n+=*(cp-*p++);
00646       n+=*(cp-*p++);
00647       n+=*(cp-*p++);
00648       n+=*(cp-*p++);
00649       n+=*(cp-*p);
00650       p+=x_size-5;
00651 
00652       n+=*(cp-*p++);
00653       n+=*(cp-*p++);
00654       n+=*(cp-*p++);
00655       n+=*(cp-*p++);
00656       n+=*(cp-*p++);
00657       n+=*(cp-*p++);
00658       n+=*(cp-*p);
00659       p+=x_size-6;
00660 
00661       n+=*(cp-*p++);
00662       n+=*(cp-*p++);
00663       n+=*(cp-*p);
00664       p+=2;
00665       n+=*(cp-*p++);
00666       n+=*(cp-*p++);
00667       n+=*(cp-*p);
00668       p+=x_size-6;
00669 
00670       n+=*(cp-*p++);
00671       n+=*(cp-*p++);
00672       n+=*(cp-*p++);
00673       n+=*(cp-*p++);
00674       n+=*(cp-*p++);
00675       n+=*(cp-*p++);
00676       n+=*(cp-*p);
00677       p+=x_size-5;
00678 
00679       n+=*(cp-*p++);
00680       n+=*(cp-*p++);
00681       n+=*(cp-*p++);
00682       n+=*(cp-*p++);
00683       n+=*(cp-*p);
00684       p+=x_size-3;
00685 
00686       n+=*(cp-*p++);
00687       n+=*(cp-*p++);
00688       n+=*(cp-*p);
00689 
00690       if (n<=max_no)
00691         r[i*x_size+j] = max_no - n;
00692     }
00693 
00694   for (i=4;i<y_size-4;i++)
00695     for (j=4;j<x_size-4;j++)
00696     {
00697       if (r[i*x_size+j]>0)
00698       {
00699         m=r[i*x_size+j];
00700         n=max_no - m;
00701         cp=bp + in[i*x_size+j];
00702 
00703         if (n>600)
00704         {
00705           p=in + (i-3)*x_size + j - 1;
00706           x=0;y=0;
00707 
00708           c=*(cp-*p++);x-=c;y-=3*c;
00709           c=*(cp-*p++);y-=3*c;
00710           c=*(cp-*p);x+=c;y-=3*c;
00711           p+=x_size-3; 
00712     
00713           c=*(cp-*p++);x-=2*c;y-=2*c;
00714           c=*(cp-*p++);x-=c;y-=2*c;
00715           c=*(cp-*p++);y-=2*c;
00716           c=*(cp-*p++);x+=c;y-=2*c;
00717           c=*(cp-*p);x+=2*c;y-=2*c;
00718           p+=x_size-5;
00719     
00720           c=*(cp-*p++);x-=3*c;y-=c;
00721           c=*(cp-*p++);x-=2*c;y-=c;
00722           c=*(cp-*p++);x-=c;y-=c;
00723           c=*(cp-*p++);y-=c;
00724           c=*(cp-*p++);x+=c;y-=c;
00725           c=*(cp-*p++);x+=2*c;y-=c;
00726           c=*(cp-*p);x+=3*c;y-=c;
00727           p+=x_size-6;
00728 
00729           c=*(cp-*p++);x-=3*c;
00730           c=*(cp-*p++);x-=2*c;
00731           c=*(cp-*p);x-=c;
00732           p+=2;
00733           c=*(cp-*p++);x+=c;
00734           c=*(cp-*p++);x+=2*c;
00735           c=*(cp-*p);x+=3*c;
00736           p+=x_size-6;
00737     
00738           c=*(cp-*p++);x-=3*c;y+=c;
00739           c=*(cp-*p++);x-=2*c;y+=c;
00740           c=*(cp-*p++);x-=c;y+=c;
00741           c=*(cp-*p++);y+=c;
00742           c=*(cp-*p++);x+=c;y+=c;
00743           c=*(cp-*p++);x+=2*c;y+=c;
00744           c=*(cp-*p);x+=3*c;y+=c;
00745           p+=x_size-5;
00746 
00747           c=*(cp-*p++);x-=2*c;y+=2*c;
00748           c=*(cp-*p++);x-=c;y+=2*c;
00749           c=*(cp-*p++);y+=2*c;
00750           c=*(cp-*p++);x+=c;y+=2*c;
00751           c=*(cp-*p);x+=2*c;y+=2*c;
00752           p+=x_size-3;
00753 
00754           c=*(cp-*p++);x-=c;y+=3*c;
00755           c=*(cp-*p++);y+=3*c;
00756           c=*(cp-*p);x+=c;y+=3*c;
00757 
00758           z = std::sqrt((float)((x*x) + (y*y)));
00759           if (z > (0.9*(float)n)) /* 0.5 */
00760     {
00761             do_symmetry=0;
00762             if (x==0)
00763               z=1000000.0f;
00764             else
00765               z=((float)y) / ((float)x);
00766             if (z < 0) { z=-z; w=-1; }
00767             else w=1;
00768             if (z < 0.5) { /* vert_edge */ a=0; b=1; }
00769             else { if (z > 2.0) { /* hor_edge */ a=1; b=0; }
00770             else { /* diag_edge */ if (w>0) { a=1; b=1; }
00771                                    else { a=-1; b=1; }}}
00772             if ( (m > r[(i+a)*x_size+j+b]) && (m >= r[(i-a)*x_size+j-b]) &&
00773                  (m > r[(i+(2*a))*x_size+j+(2*b)]) && (m >= r[(i-(2*a))*x_size+j-(2*b)]) )
00774               mid[i*x_size+j] = 1;
00775           }
00776           else
00777             do_symmetry=1;
00778         }
00779         else 
00780           do_symmetry=1;
00781 
00782         if (do_symmetry==1)
00783   { 
00784           p=in + (i-3)*x_size + j - 1;
00785           x=0; y=0; w=0;
00786 
00787           /*   |      \
00788                y  -x-  w
00789                |        \   */
00790 
00791           c=*(cp-*p++);x+=c;y+=9*c;w+=3*c;
00792           c=*(cp-*p++);y+=9*c;
00793           c=*(cp-*p);x+=c;y+=9*c;w-=3*c;
00794           p+=x_size-3; 
00795   
00796           c=*(cp-*p++);x+=4*c;y+=4*c;w+=4*c;
00797           c=*(cp-*p++);x+=c;y+=4*c;w+=2*c;
00798           c=*(cp-*p++);y+=4*c;
00799           c=*(cp-*p++);x+=c;y+=4*c;w-=2*c;
00800           c=*(cp-*p);x+=4*c;y+=4*c;w-=4*c;
00801           p+=x_size-5;
00802     
00803           c=*(cp-*p++);x+=9*c;y+=c;w+=3*c;
00804           c=*(cp-*p++);x+=4*c;y+=c;w+=2*c;
00805           c=*(cp-*p++);x+=c;y+=c;w+=c;
00806           c=*(cp-*p++);y+=c;
00807           c=*(cp-*p++);x+=c;y+=c;w-=c;
00808           c=*(cp-*p++);x+=4*c;y+=c;w-=2*c;
00809           c=*(cp-*p);x+=9*c;y+=c;w-=3*c;
00810           p+=x_size-6;
00811 
00812           c=*(cp-*p++);x+=9*c;
00813           c=*(cp-*p++);x+=4*c;
00814           c=*(cp-*p);x+=c;
00815           p+=2;
00816           c=*(cp-*p++);x+=c;
00817           c=*(cp-*p++);x+=4*c;
00818           c=*(cp-*p);x+=9*c;
00819           p+=x_size-6;
00820     
00821           c=*(cp-*p++);x+=9*c;y+=c;w-=3*c;
00822           c=*(cp-*p++);x+=4*c;y+=c;w-=2*c;
00823           c=*(cp-*p++);x+=c;y+=c;w-=c;
00824           c=*(cp-*p++);y+=c;
00825           c=*(cp-*p++);x+=c;y+=c;w+=c;
00826           c=*(cp-*p++);x+=4*c;y+=c;w+=2*c;
00827           c=*(cp-*p);x+=9*c;y+=c;w+=3*c;
00828           p+=x_size-5;
00829  
00830           c=*(cp-*p++);x+=4*c;y+=4*c;w-=4*c;
00831           c=*(cp-*p++);x+=c;y+=4*c;w-=2*c;
00832           c=*(cp-*p++);y+=4*c;
00833           c=*(cp-*p++);x+=c;y+=4*c;w+=2*c;
00834           c=*(cp-*p);x+=4*c;y+=4*c;w+=4*c;
00835           p+=x_size-3;
00836 
00837           c=*(cp-*p++);x+=c;y+=9*c;w-=3*c;
00838           c=*(cp-*p++);y+=9*c;
00839           c=*(cp-*p);x+=c;y+=9*c;w+=3*c;
00840 
00841           if (y==0)
00842             z = 1000000.0f;
00843           else
00844             z = ((float)x) / ((float)y);
00845           if (z < 0.5) { /* vertical */ a=0; b=1; }
00846           else { if (z > 2.0) { /* horizontal */ a=1; b=0; }
00847           else { /* diagonal */ if (w>0) { a=-1; b=1; }
00848                                 else { a=1; b=1; }}}
00849           if ( (m > r[(i+a)*x_size+j+b]) && (m >= r[(i-a)*x_size+j-b]) &&
00850                (m > r[(i+(2*a))*x_size+j+(2*b)]) && (m >= r[(i-(2*a))*x_size+j-(2*b)]) )
00851             mid[i*x_size+j] = 2;  
00852         }
00853       }
00854     }
00855 }
00856 
00857 
00858 /* {{{ edge_draw(in,corner_list,drawing_mode) */
00859 
00860 /*edge_draw(in,mid,x_size,y_size,drawing_mode)
00861   uchar *in, *mid;
00862   int x_size, y_size, drawing_mode;*/
00863 void edge_draw(uchar *in, uchar *mid, int x_size, int y_size, int drawing_mode)
00864 {
00865 int   i;
00866 uchar *inp, *midp;
00867 
00868   if (drawing_mode==0)
00869   {
00870     /* mark 3x3 white block around each edge point */
00871     midp=mid;
00872     for (i=0; i<x_size*y_size; i++)
00873     {
00874       if (*midp<8) 
00875       {
00876         inp = in + (midp - mid) - x_size - 1;
00877         *inp++=255; *inp++=255; *inp=255; inp+=x_size-2;
00878         *inp++=255; ++inp;      *inp=255; inp+=x_size-2;
00879         *inp++=255; *inp++=255; *inp=255;
00880       }
00881       midp++;
00882     }
00883   }
00884 
00885   /* now mark 1 black pixel at each edge point */
00886   midp=mid;
00887   for (i=0; i<x_size*y_size; i++)
00888   {
00889     if (*midp<8) 
00890       //*(in + (midp - mid)) = 0;
00891       *(in + (midp - mid)) = 1;
00892     midp++;
00893   }
00894 }
00895 
00896 } // namespace

DualCoding 5.1CVS
Generated Mon May 9 04:56:28 2016 by Doxygen 1.6.3