n+=*(cp-*p++);      n+=*(cp-*p);      p+=2;      n+=*(cp-*p++);      n+=*(cp-*p++);      n+=*(cp-*p);      p+=x_size-6;      n+=*(cp-*p++);      n+=*(cp-*p++);      n+=*(cp-*p++);      n+=*(cp-*p++);      n+=*(cp-*p++);      n+=*(cp-*p++);      n+=*(cp-*p);      p+=x_size-5;      n+=*(cp-*p++);      n+=*(cp-*p++);      n+=*(cp-*p++);      n+=*(cp-*p++);      n+=*(cp-*p);      p+=x_size-3;      n+=*(cp-*p++);      n+=*(cp-*p++);      n+=*(cp-*p);      if (n<=max_no)        r[i*x_size+j] = max_no - n;    }}/* }}} *//* {{{ susan_principle_small(in,r,bp,max_no,x_size,y_size) */susan_principle_small(in,r,bp,max_no,x_size,y_size)  uchar *in, *bp;  int   *r, max_no, x_size, y_size;{int   i, j, n;uchar *p,*cp;  memset (r,0,x_size * y_size * sizeof(int));  max_no = 730; /* ho hum ;) */  for (i=1;i<y_size-1;i++)    for (j=1;j<x_size-1;j++)    {      n=100;      p=in + (i-1)*x_size + j - 1;      cp=bp + in[i*x_size+j];      n+=*(cp-*p++);      n+=*(cp-*p++);      n+=*(cp-*p);      p+=x_size-2;       n+=*(cp-*p);      p+=2;      n+=*(cp-*p);      p+=x_size-2;      n+=*(cp-*p++);      n+=*(cp-*p++);      n+=*(cp-*p);      if (n<=max_no)        r[i*x_size+j] = max_no - n;    }}/* }}} *//* }}} *//* {{{ smoothing *//* {{{ median(in,i,j,x_size) */uchar median(in,i,j,x_size)  uchar *in;  int   i, j, x_size;{int p[8],k,l,tmp;  p[0]=in[(i-1)*x_size+j-1];  p[1]=in[(i-1)*x_size+j  ];  p[2]=in[(i-1)*x_size+j+1];  p[3]=in[(i  )*x_size+j-1];  p[4]=in[(i  )*x_size+j+1];  p[5]=in[(i+1)*x_size+j-1];  p[6]=in[(i+1)*x_size+j  ];  p[7]=in[(i+1)*x_size+j+1];  for(k=0; k<7; k++)    for(l=0; l<(7-k); l++)      if (p[l]>p[l+1])      {        tmp=p[l]; p[l]=p[l+1]; p[l+1]=tmp;      }  return( (p[3]+p[4]) / 2 );}/* }}} *//* {{{ enlarge(in,tmp_image,x_size,y_size,border) *//* this enlarges "in" so that borders can be dealt with easily */enlarge(in,tmp_image,x_size,y_size,border)  uchar **in;  uchar *tmp_image;  int   *x_size, *y_size, border;{int   i, j;  for(i=0; i<*y_size; i++)   /* copy *in into tmp_image */    memcpy(tmp_image+(i+border)*(*x_size+2*border)+border, *in+i* *x_size, *x_size);  for(i=0; i<border; i++) /* copy top and bottom rows; invert as many as necessary */  {    memcpy(tmp_image+(border-1-i)*(*x_size+2*border)+border,*in+i* *x_size,*x_size);    memcpy(tmp_image+(*y_size+border+i)*(*x_size+2*border)+border,*in+(*y_size-i-1)* *x_size,*x_size);  }  for(i=0; i<border; i++) /* copy left and right columns */    for(j=0; j<*y_size+2*border; j++)    {      tmp_image[j*(*x_size+2*border)+border-1-i]=tmp_image[j*(*x_size+2*border)+border+i];      tmp_image[j*(*x_size+2*border)+ *x_size+border+i]=tmp_image[j*(*x_size+2*border)+ *x_size+border-1-i];    }  *x_size+=2*border;  /* alter image size */  *y_size+=2*border;  *in=tmp_image;      /* repoint in */}/* }}} *//* {{{ void susan_smoothing(three_by_three,in,dt,x_size,y_size,bp) */void susan_smoothing(three_by_three,in,dt,x_size,y_size,bp)  int   three_by_three, x_size, y_size;  uchar *in, *bp;  float dt;{/* {{{ vars */float temp;int   n_max, increment, mask_size,      i,j,x,y,area,brightness,tmp,centre;uchar *ip, *dp, *dpt, *cp, *out=in,      *tmp_image;TOTAL_TYPE total;/* }}} */  /* {{{ setup larger image and border sizes */  if (three_by_three==0)    mask_size = ((int)(1.5 * dt)) + 1;  else    mask_size = 1;  total=0.1; /* test for total's type */  if ( (dt>15) && (total==0) )  {    printf("Distance_thresh (%f) too big for integer arithmetic.\n",dt);    printf("Either reduce it to <=15 or recompile with variable \"total\"\n");    printf("as a float: see top \"defines\" section.\n");    exit(0);  }  if ( (2*mask_size+1>x_size) || (2*mask_size+1>y_size) )  {    printf("Mask size (1.5*distance_thresh+1=%d) too big for image (%dx%d).\n",mask_size,x_size,y_size);    exit(0);  }  tmp_image = (uchar *) malloc( (x_size+mask_size*2) * (y_size+mask_size*2) );  enlarge(&in,tmp_image,&x_size,&y_size,mask_size);/* }}} */  if (three_by_three==0)  {     /* large Gaussian masks */    /* {{{ setup distance lut */  n_max = (mask_size*2) + 1;  increment = x_size - n_max;  dp     = (unsigned char *)malloc(n_max*n_max);  dpt    = dp;  temp   = -(dt*dt);  for(i=-mask_size; i<=mask_size; i++)    for(j=-mask_size; j<=mask_size; j++)    {      x = (int) (100.0 * exp( ((float)((i*i)+(j*j))) / temp ));      *dpt++ = (unsigned char)x;    }/* }}} */    /* {{{ main section */  for (i=mask_size;i<y_size-mask_size;i++)  {    for (j=mask_size;j<x_size-mask_size;j++)    {      area = 0;      total = 0;      dpt = dp;      ip = in + ((i-mask_size)*x_size) + j - mask_size;      centre = in[i*x_size+j];      cp = bp + centre;      for(y=-mask_size; y<=mask_size; y++)      {        for(x=-mask_size; x<=mask_size; x++)	{          brightness = *ip++;          tmp = *dpt++ * *(cp-brightness);          area += tmp;          total += tmp * brightness;        }        ip += increment;      }      tmp = area-10000;      if (tmp==0)        *out++=median(in,i,j,x_size);      else        *out++=((total-(centre*10000))/tmp);    }  }/* }}} */  }  else  {     /* 3x3 constant mask */    /* {{{ main section */  for (i=1;i<y_size-1;i++)  {    for (j=1;j<x_size-1;j++)    {      area = 0;      total = 0;      ip = in + ((i-1)*x_size) + j - 1;      centre = in[i*x_size+j];      cp = bp + centre;      brightness=*ip++; tmp=*(cp-brightness); area += tmp; total += tmp * brightness;      brightness=*ip++; tmp=*(cp-brightness); area += tmp; total += tmp * brightness;      brightness=*ip; tmp=*(cp-brightness); area += tmp; total += tmp * brightness;      ip += x_size-2;      brightness=*ip++; tmp=*(cp-brightness); area += tmp; total += tmp * brightness;      brightness=*ip++; tmp=*(cp-brightness); area += tmp; total += tmp * brightness;      brightness=*ip; tmp=*(cp-brightness); area += tmp; total += tmp * brightness;      ip += x_size-2;      brightness=*ip++; tmp=*(cp-brightness); area += tmp; total += tmp * brightness;      brightness=*ip++; tmp=*(cp-brightness); area += tmp; total += tmp * brightness;      brightness=*ip; tmp=*(cp-brightness); area += tmp; total += tmp * brightness;      tmp = area-100;      if (tmp==0)        *out++=median(in,i,j,x_size);      else        *out++=(total-(centre*100))/tmp;    }  }/* }}} */  }}/* }}} *//* }}} *//* {{{ edges *//* {{{ edge_draw(in,corner_list,drawing_mode) */edge_draw(in,mid,x_size,y_size,drawing_mode)  uchar *in, *mid;  int x_size, y_size, drawing_mode;{int   i;uchar *inp, *midp;  if (drawing_mode==0)  {    /* mark 3x3 white block around each edge point */    midp=mid;    for (i=0; i<x_size*y_size; i++)    {      if (*midp<8)       {        inp = in + (midp - mid) - x_size - 1;        *inp++=255; *inp++=255; *inp=255; inp+=x_size-2;        *inp++=255; *inp++;     *inp=255; inp+=x_size-2;        *inp++=255; *inp++=255; *inp=255;      }      midp++;    }  }  /* now mark 1 black pixel at each edge point */  midp=mid;  for (i=0; i<x_size*y_size; i++)  {    if (*midp<8)       *(in + (midp - mid)) = 0;    midp++;  }}/* }}} *//* {{{ susan_thin(r,mid,x_size,y_size) *//* only one pass is needed as i,j are decremented if necessary to go   back and do bits again */susan_thin(r,mid,x_size,y_size)  uchar *mid;  int   *r, x_size, y_size;{int   l[9], centre, nlinks, npieces,      b01, b12, b21, b10,      p1, p2, p3, p4,      b00, b02, b20, b22,      m, n, a, b, x, y, i, j;uchar *mp;  for (i=4;i<y_size-4;i++)    for (j=4;j<x_size-4;j++)      if (mid[i*x_size+j]<8)      {        centre = r[i*x_size+j];        /* {{{ count number of neighbours */        mp=mid + (i-1)*x_size + j-1;        n = (*mp<8) +            (*(mp+1)<8) +            (*(mp+2)<8) +            (*(mp+x_size)<8) +            (*(mp+x_size+2)<8) +            (*(mp+x_size+x_size)<8) +            (*(mp+x_size+x_size+1)<8) +            (*(mp+x_size+x_size+2)<8);/* }}} */        /* {{{ n==0 no neighbours - remove point */        if (n==0)          mid[i*x_size+j]=100;/* }}} */        /* {{{ n==1 - extend line if I can */        /* extension is only allowed a few times - the value of mid is used to control this */        if ( (n==1) && (mid[i*x_size+j]<6) )        {          /* find maximum neighbour weighted in direction opposite the             neighbour already present. e.g.             have: O O O  weight r by 0 2 3                   X X O              0 0 4                   O O O              0 2 3     */          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];          l[3]=r[(i  )*x_size+j-1]; l[4]=0;                 l[5]=r[(i  )*x_size+j+1];          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];          if (mid[(i-1)*x_size+j-1]<8)        { l[0]=0; l[1]=0; l[3]=0; l[2]*=2;                                                 l[6]*=2; l[5]*=3; l[7]*=3; l[8]*=4; }          else { if (mid[(i-1)*x_size+j]<8)   { l[1]=0; l[0]=0; l[2]=0; l[3]*=2;                                                 l[5]*=2; l[6]*=3; l[8]*=3; l[7]*=4; }          else { if (mid[(i-1)*x_size+j+1]<8) { l[2]=0; l[1]=0; l[5]=0; l[0]*=2;                                                 l[8]*=2; l[3]*=3; l[7]*=3; l[6]*=4; }          else { if (mid[(i)*x_size+j-1]<8)   { l[3]=0; l[0]=0; l[6]=0; l[1]*=2;                                                 l[7]*=2; l[2]*=3; l[8]*=3; l[5]*=4; }          else { if (mid[(i)*x_size+j+1]<8)   { l[5]=0; l[2]=0; l[8]=0; l[1]*=2;                                                 l[7]*=2; l[0]*=3; l[6]*=3; l[3]*=4; }          else { if (mid[(i+1)*x_size+j-1]<8) { l[6]=0; l[3]=0; l[7]=0; l[0]*=2;                                                 l[8]*=2; l[1]*=3; l[5]*=3; l[2]*=4; }          else { if (mid[(i+1)*x_size+j]<8)   { l[7]=0; l[6]=0; l[8]=0; l[3]*=2;                                                 l[5]*=2; l[0]*=3; l[2]*=3; l[1]*=4; }          else { if (mid[(i+1)*x_size+j+1]<8) { l[8]=0; l[5]=0; l[7]=0; l[6]*=2;                                                 l[2]*=2; l[1]*=3; l[3]*=3; l[0]*=4; } }}}}}}}          m=0;     /* find the highest point */          for(y=0; y<3; y++)            for(x=0; x<3; x++)              if (l[y+y+y+x]>m) { m=l[y+y+y+x]; a=y; b=x; }          if (m>0)          {            if (mid[i*x_size+j]<4)              mid[(i+a-1)*x_size+j+b-1] = 4;            else              mid[(i+a-1)*x_size+j+b-1] = mid[i*x_size+j]+1;            if ( (a+a+b) < 3 ) /* need to jump back in image */	    {              i+=a-1;              j+=b-2;              if (i<4) i=4;              if (j<4) j=4;	    }	  }        }/* }}} */        /* {{{ n==2 */        if (n==2)	{          /* put in a bit here to straighten edges */          b00 = mid[(i-1)*x_size+j-1]<8; /* corners of 3x3 */          b02 = mid[(i-1)*x_size+j+1]<8;	  b20 = mid[(i+1)*x_size+j-1]<8;          b22 = mid[(i+1)*x_size+j+1]<8;          if ( ((b00+b02+b20+b22)==2) && ((b00|b22)&(b02|b20)))	  {  /* case: move a point back into line.                e.g. X O X  CAN  become X X X                     O X O              O O O                     O O O              O O O    */            if (b00) 	    {              if (b02) { x=0; y=-1; }              else     { x=-1; y=0; }	    }            else	    {              if (b02) { x=1; y=0; }              else     { x=0; y=1; }	    }            if (((float)r[(i+y)*x_size+j+x]/(float)centre) > 0.7)	    {              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) ) ||                   ( (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) ) )	      {                mid[(i)*x_size+j]=100;                mid[(i+y)*x_size+j+x]=3;  /* no jumping needed */	      }	    }	  }          else          {            b01 = mid[(i-1)*x_size+j  ]<8;            b12 = mid[(i  )*x_size+j+1]<8;            b21 = mid[(i+1)*x_size+j  ]<8;            b10 = mid[(i  )*x_size+j-1]<8;            /* {{{ right angle ends - not currently used */#ifdef IGNORETHIS            if ( (b00&b01)|(b00&b10)|(b02&b01)|(b02&b12)|(b20&b10)|(b20&b21)|(b22&b21)|(b22&b12) )	    { /* case; right angle ends. clean up.                 e.g.; X X O  CAN  become X X O                       O X O              O O O                       O O O              O O O        */              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)&                                    ((b00&((2*r[(i-1)*x_size+j+1])>centre))|(b02&((2*r[(i-1)*x_size+j-1])>centre)))) |                   ((b10)&(mid[(i-1)*x_size+j-2]>7)&(mid[(i)*x_size+j-2]>7)&(mid[(i+1)*x_size+j-2]>7)&                                    ((b00&((2*r[(i+1)*x_size+j-1])>centre))|(b20&((2*r[(i-1)*x_size+j-1])>centre)))) |                   ((b12)&(mid[(i-1)*x_size+j+2]>7)&(mid[(i)*x_size+j+2]>7)&(mid[(i+1)*x_size+j+2]>7)&                                    ((b02&((2*r[(i+1)*x_size+j+1])>centre))|(b22&((2*r[(i-1)*x_size+j+1])>centre)))) |                   ((b21)&(mid[(i+2)*x_size+j-1]>7)&(mid[(i+2)*x_size+j]>7)&(mid[(i+2)*x_size+j+1]>7)&                                    ((b20&((2*r[(i+1)*x_size+j+1])>centre))|(b22&((2*r[(i+1)*x_size+j-1])>centre)))) )	      {                mid[(i)*x_size+j]=100;                if (b10&b20) j-=2;                if (b00|b01|b02) { i--; j-=2; }  	      }	    }#endif/* }}} */            if ( ((b01+b12+b21+b10)==2) && ((b10|b12)&(b01|b21)) &&                 ((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)))|                (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)))) )	    { /* case; clears odd right angles.                 e.g.; O O O  becomes O O O                       X X O          X O O                       O X O          O X O     */              mid[(i)*x_size+j]=100;              i--;               /* jump back */              j-=2;              if (i<4) i=4;              if (j<4) j=4;	    }	  }	}/* }}} */        /* {{{ n>2 the thinning is done here without breaking connectivity */        if (n>2)        {          b01 = mid[(i-1)*x_size+j  ]<8;          b12 = mid[(i  )*x_size+j+1]<8;          b21 = mid[(i+1)*x_size+j  ]<8;          b10 = mid[(i  )*x_size+j-1]<8;          if((b01+b12+b21+b10)>1)          {            b00 = mid[(i-1)*x_size+j-1]<8;            b02 = mid[(i-1)*x_size+j+1]<8;	    b20 = mid[(i+1)*x_size+j-1]<8;	    b22 = mid[(i+1)*x_size+j+1]<8;            p1 = b00 | b01;            p2 = b02 | b12;            p3 = b22 | b21;            p4 = b20 | b10;            if( ((p1 + p2 + p3 + p4) - ((b01 & p2)+(b12 & p3)+(b21 & p4)+(b10 & p1))) < 2)            {              mid[(i)*x_size+j]=100;              i--;              j-=2;              if (i<4) i=4;              if (j<4) j=4;            }          }        }/* }}} */      }}/* }}} *//* {{{ susan_edges(in,r,sf,max_no,out) */