ccribbons.c

Computing 2D Skeleton

/* Program:  finds ribbon-like parts                        */


/* TRACKS BRANCHES by following ridges of grd  */

 /*  Problem: Does not handle triple points very well       */

/* problem: cannot distinguish between short necks
(# pixel = beta) and short protrusion; the program
 ignores both (see test().) Perhaps, both should be
 filled in as bdry pts because when two parts are
 separated by a very narrow neck, they might as well
 be treated as disjoint. 
 	Alternatively, if the branch can be continued from
 	two ends, it is not a short protrusion. check if
 	the complement of sk1 in a small nbhd of the complete
 	skeleton is connected or not */
 
 /* possible improvements: 
 
 it might be possible to simplify fill_gaps(). At present,
 to bridge a gap in sk1, it tracks along ridges of grd in sk2
 in order to find a path along the gap with minimum cost.
 NOTE: We cannot simply track along sk2
 in case there is a junction along the way. Nor can we
 trace gaps in sk1 by tracking values of v since, if a
 saddle point is missed, v will decrease towards the saddle
 point and then increase.
 
 In prune(),if the branch is too long (here, >10 pixels),
 calculate the actual length by dividing the branch into
 tow segments as done in fill_gaps().
  
 /*


/* Assumes that image is piecewise constant                 */


/*                                                          */
/*                                                          */
/* Inputs:  the image and the ABSOLUTE distance function    */
/*                                                          */
/*                                                          */
/* Date:     August, 2002                                   */
/*                                                          */
/*                                                          */





#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <time.h>
#include <fcntl.h>
#include <unistd.h>
#include <sioux.h>
#include <MacTypes.h>
#include <Fonts.h>
#include <math.h>

#define ydim 364  /* vertical y-dimension: index i  */
#define xdim 784  /* horizontal x-dimension:  index j  */


#define upper_threshold 70
#define middle_threshold 45
#define lower_threshold 10

/*
 (1) upper threshold = degree angle between 0 and 90
                   = minimum angle that a chord must make
                     with the tangent to shape bdry

 (2) middle threshold is used to create points bridging
 gaps in the skeleton obtained by thresholding at
 upper_threshold. If middle threshold is too high, it will
 miss some of the gaps. If it is too low,then it
 finds too many points going beyond the gap
 
 (3) lower_threshold is used to constuct fairly full skeleton 
 which is then used to find the endpts of the branches
 of the skeleon obtained after bridging short gaps created
 by thresholding. 
*/

#define beta 2.01
#define L 20
/* beta is the cost of each connected component of sk1 */
/* eliminates branches with number of pixels<=L if their cost
is beyond the threshold */ 

#define N 15 /* max size for local nbhd in ridge_gr()   */
/* N must be odd */


#define MAX(A,B) ((A)>(B) ? (A) : (B))
#define MIN(A,B) ((A)>(B) ? (B) : (A))
#define ABS(A) ((A)>=0 ? (A) : -(A))

/* file names */
FILE *infile,*outfile;
char buf[64],buf1[64],buf2[64],inname[128];
char d_name[128],outname[128];

/* Program variables */
unsigned char pic[ydim][xdim];
unsigned char outsk[ydim][xdim];
unsigned char outsgmt[ydim][xdim];
unsigned char outgr[ydim][xdim];
unsigned char sk[ydim][xdim];
unsigned char sgmt[ydim][xdim];
unsigned char nbhd_size[ydim][xdim];
unsigned char bdry[ydim][xdim];
unsigned char endpts[ydim][xdim];
unsigned char tmp[ydim][xdim];
int mark[2][ydim][xdim],branchcnt[2];
int Iarray[2][xdim*ydim],Jarray[2][xdim*ydim],arr_len[2];
float distance[xdim*ydim];
int length[2][xdim*ydim];
float cost[2][xdim*ydim];

float v[ydim][xdim];
float gr[ydim][xdim],grd[ydim][xdim];
float tol[(N-1)/2],tol_d[(N-1)/2]; /* margin of error in angles *

/* global variables */
int ibb,inn,jbb,jnn;

void local_grad();
void convert_gr();
void ridge_gr();
void construct_sk();
void initialize_sk();
void ave_cost(int,int);
void confirm_endpts();
void fill_gaps();
void complete_sk();
void prune();
void connect(int,int,int,int,int);
void segment();
void output_results();
void run();

 
   
   void local_grad()
   /* 
    CALCULATES grad USING 4 DIFFERENT WAYS: 
    given coordinates and then rotating it by angle
    whose tangent equals 0.5,1 and 2   */
    
   {
      int i,j,k,i1,j1,i2,j2,i3,j3;
      int count,flag;
      double x,y,q,q1; 
      
      for (i=0;i<ydim;i++) for (j=0;j<xdim;j++)
      {
        if (bdry[i][j]>0) gr[i][j]=1;
        else if (i*j*(ydim-1-i)*(xdim-1-j)==0) gr[i][j]=1;
        else
        {
         ibb = i-1; inn = i+1; jbb = j-1; jnn = j+1;
         flag=0;
         k = pic[i][j];
         i1=i;j1=jbb;
         i2=ibb;j2=jbb;
         i3=ibb;j3=j;
         count=bdry[i1][j1]+bdry[i2][j2]+bdry[i3][j3];
         if (count==3)
         {
           if (pic[i2][j2]==k) gr[i][j]=0.707107;
           else gr[i][j]=1;
           flag = 1;
         }

         
         if (flag==0)
         {
           i1=i3;j1=j3;
           i2=ibb;j2=jnn;
           i3=i;j3=jnn;
           count=bdry[i1][j1]+bdry[i2][j2]+bdry[i3][j3];
           if (count==3)
           {
             if (pic[i2][j2]==k) gr[i][j]=0.707107;
             else gr[i][j]=1;
             flag = 1;
           }
         }
         
         if (flag==0)
         {
           i1=i3;j1=j3;
           i2=inn;j2=jnn;
           i3=inn;j3=j;
           count=bdry[i1][j1]+bdry[i2][j2]+bdry[i3][j3];
           if (count==3)
           {
             if (pic[i2][j2]==k) gr[i][j]=0.707107;
             else gr[i][j]=1;
             flag = 1;
           }
         }
         
         if (flag==0)
         {
           i1=i3;j1=j3;
           i2=inn;j2=jbb;
           i3=i;j3=jbb;
           count=bdry[i1][j1]+bdry[i2][j2]+bdry[i3][j3];
           if (count==3)
           {
             if (pic[i2][j2]==k) gr[i][j]=0.707107;
             else gr[i][j]=1;
             flag = 1;
           }
         }
        
         
         if (flag==0)
         {
           x = v[i][jnn]-v[i][jbb];
           y = v[inn][j]-v[ibb][j];
         
           q = (x*x+y*y)/4;
        
           x = v[inn][jnn]-v[ibb][jbb];
           y = v[inn][jbb]-v[ibb][jnn];
           q1 = (x*x+y*y)/8;
           if (q1<q) q = q1; 
         
           x = (v[inn][jbb]+v[inn][j]-v[ibb][j]-v[ibb][jnn])/2;
           y = (v[inn][jnn]+v[i][jnn]-v[i][jbb]-v[ibb][jbb])/2;
           q1 = (x*x+y*y)/5;
           if (q1<q) q = q1; 
         
           x = (v[inn][jbb]+v[i][jbb]-v[i][jnn]-v[ibb][jnn])/2;
           y = (v[inn][jnn]+v[inn][j]-v[ibb][j]-v[ibb][jbb])/2;
           q1 = (x*x+y*y)/5;
           if (q1<q) q = q1; 
         
           if (q>1) q = 1; 
           if (q>0) gr[i][j] = sqrt(q);
           else gr[i][j] = 0;
         }
         gr[i][j] = MIN(gr[i][j],1);
        }
      }
   }
      
  
   
  
 
 void ridge_gr()
 /* traveses the boundary pixels of NxN nbhds */
 {
    int i,j,m,n,k,kk,p,q;
    int k0,k1,count;
    int h,hh,H;
    int flag,reset_flag,bdry_flag;
    int ii[4*N-4],jj[4*N-4],check[4*N-4];
    float v0,z,d[4*N-4];
    float slope[4*N-4],last,now,next;
    float grx,grmn,accept;

    
    local_grad(); /* makes a rough estimate of ridge_gr */
                       /* tends to be more accurate near jcns */ 
    
    for (i=0;i<ydim;i++) for (j=0;j<xdim;j++) nbhd_size[i][j]=1;
    
    h = (N-1)/2;
    H = 8*h;
    for (p=-h;p<h+1;p++) for (q=-h;q<h+1;q++)
    if ((p+h)*(p-h)*(q+h)*(q-h)==0)
    {
      k = p+q+2*h;
      if (p-q>0) k = H-k;
      z = p*p+q*q;
      d[k] = sqrt((double) z); /* distance from the origin */
      ii[k] = p; jj[k] = q; 
    }
  
   
   
    /* begin scanning */
    for (i=h;i<ydim-h;i++) for (j=h;j<xdim-h;j++)
    if (bdry[i][j]==0)  if (gr[i][j]<0.99) 
    {     
     
     /* adjust the nbhd size to avoid bdry */
     reset_flag = 0;
     do
     {
/* printf("h=%d\n",h);*/
       bdry_flag = 0;
       hh = h+1;
       for (m=i-hh;m<i+hh+1;m++) for (n=j-hh;n<j+hh+1;n++)
          if (bdry[m][n]>0) bdry_flag = 1; 
  /*  if (bdry[m][n]>0)
{printf("bdry pt at (%d %d), v=%f\n",m,n,v[m][n]);flag=1;} */
       if (bdry_flag>0) {h--, reset_flag=1;}
     }
     while ((bdry_flag>0)&&(h>1));
     
     
     
     if (bdry_flag==0) 
     { 
       nbhd_size[i][j] = h;
       if (reset_flag>0) /* recalculate nbhd geometry */
       {
         H = 8*h;
         for (p=-h;p<h+1;p++) for (q=-h;q<h+1;q++)
         if ((p+h)*(p-h)*(q+h)*(q-h)==0)
         {
           k = p+q+2*h;
           if (p-q>0) k = H-k;
           z = p*p+q*q;
           d[k] = sqrt((double) z); /* distance from the origin */
           ii[k] = p; jj[k] = q; 
         }
       }
        
      /* calculate directional gradients */
      /* direction points from the  perimeter to the center (i,j) */ 
       v0 = v[i][j];
       for (m=i-h;m<i+h+1;m++) for (n=j-h;n<j+h+1;n++)
       if ((m-i+h)*(m-i-h)*(n-j+h)*(n-j-h)==0)
       /* scans only the perimeter of the nbhd */
       {
          p = m-i; q = n-j;
          k = p+q+2*h;
          if (p-q>0) k = H-k;
          slope[k] = (v0-v[m][n])/d[k];
       }
       
       /* find local maxima of slopes around the nbhd perimeter */
       for (k=0;k<H;k++) check[k]=0;
       last = slope[H-1];
       now = slope[0];
       z = -2; /* z tracks value of global max */
       for (k=0;k<H;k++)
       {
         next = slope[(k+1)%H];
         if (now-last>0) if (now-next>0) /* local max */
         {
           check[k] = 1; /* mark the direction */
           if (now>z) {z = now; kk = k;} /* current global max */
         }
         last = now; now = next;
       }
 /*  if (i==442) if (j==369)  */
/* if (i>286) if (i<289) if (j==103) */
/* {  printf("\nat (%d %d), nbhd: %d\n",i,j,h);
   for (k=0;k<H;k++) if (check[k]>0)
printf("k:%d (%d %d),gr:%f,chk:%d\n",k,ii[k],jj[k],slope[k],check[k]); 
 }*/
  
    /* calculate cosine of angles between successive maximal rays */

        accept=(h-1)/sqrt((double) h*h+0.25);
 /* Maximal rays should have slope of 1; but due to
 inaccuarcy in slopes, arising because the maximal rays may
 be off the grid line contributing error factor of h/sqrt(h*h+.5*.5)
 and because the point (i,j) may be off the ridge contributing the
 factor (h-1)/h, accept if the slope > (h-1)/sqrt((double) h*h+0.25)
 */
       count = 0; /* number of maximal rays */
       grmn = 1;
       /* grmn tracks minimum value; */
      
       if (z>accept) /* first ray to shape bdry */
       {
        k0 = kk;
        for (k=kk;k<kk+H;k++)
        if ((k%H)==k0)
        {
          k0 = k%H;  /* current ray */
          flag = 0;
          /* find next ray to shape bdry */
          for (q=1;q<H;q++) 
          {
            if (check[(k+q)%H]>0) /* is local max */
            if (slope[(k+q)%H]>accept) 
            {
              flag = 1;
              count++;
              k1 = (k+q)%H; /* next ray */
              break;
            }  
          }
        
          if (flag>0) /* calc ridge gradient */
          {
            grx = (ii[k0]*ii[k1]+jj[k0]*jj[k1])/(d[k0]*d[k1]);
            grx = ABS(1+grx)/2;
            if (grx>0.0001) grx = sqrt((double) grx);
            if (grx<grmn)  grmn=grx;
            
/* if (i==442) if (j==369)
printf("(%d %d) k0 %d k1 %d, kk %d gr %f %f\n",i,j,k0,k1,count,grx,gr[i][j]); */

            k0 = k1; /* new current ray */
                     
          } /* else k0 remains unchanged */
        }
       }

       if (count>=2) if (grmn<gr[i][j]) gr[i][j]=grmn;
       /* if count<2; the point (i,j) is not on the ridge */ 
     }
     if (reset_flag>0)
     {
       h = (N-1)/2;
       H = 8*h;
       for (p=-h;p<h+1;p++) for (q=-h;q<h+1;q++)
       if ((p+h)*(p-h)*(q+h)*(q-h)==0)
       {
         k = p+q+2*h;
         if (p-q>0) k = H-k;
         z = p*p+q*q;
         d[k] = sqrt((double) z); /* distance from the origin */
         ii[k] = p; jj[k] = q;
       }
     }
   }
  
   convert_gr();
/*  p=364; q=294; 
 for (i=p-1;i<p+2;i++) for (j=q-1;j<q+2;j++)
printf("grd(%d %d)=%f\n",i,j,grd[i][j]);  */  
 }
 



 
 void convert_gr()
 /* converts gr to angle in degrees */
 /* removes skeleton points too near the bdry and
 bear small angles */
 {
   int i,j,h;
   float x;
   
   h = (N-1)/2;
   for (j=2;j<=h;j++)
   {
     tol[j] = atanf((float) (0.5/j));
     tol_d[j] =  57.2958*tol[j];
   }
   tol[1] = tol[2];
   tol_d[1]=tol_d[2];
   
   /* convert to angle that the chord makes with
   tangent to shape bdry (radians) */
   for (i=0;i<ydim;i++) for (j=0;j<xdim;j++)
   {
      x = acosf(gr[i][j]);
      
      /* threshold weak edges */
      h = nbhd_size[i][j];
      if (h==1)
      if (x<tol[h]) x = 0;
      
      /* convert to degrees */
      grd[i][j] = 57.2958*x; 
   }
 }
        
   
 
  
 void construct_sk()
 {
    initialize_sk(); 
    fill_gaps();
    complete_sk();
    prune(); 
/*    segment(); */
 }
 
 
 
  
 
 void initialize_sk()