ccfeatures.c

Computing 2D Skeleton

   {
     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()
 {
    /* the first two steps identify sk1. sk2 found
    in the first two steps is discarded in step 3 */
    initialize_sk(); /* step 1 */
    fill_gaps();     /* step 2 */
    complete_sk();   /* step 3 */
    prune();         /* step 4 */
 }
 
 
 
  
 
 void initialize_sk()
 {
     int i,j,m,n;
     int mk,mk1,mk2,len0,len1;
     float a,b;

     a = upper_threshold;
     b= middle_threshold;
     /* first try to eliminate short branches by following
     ridges of grd */
     len0=len1=0;
     for (i=0;i<ydim;i++) for (j=0;j<xdim;j++)
     {
       if (grd[i][j]>=(a-tol_d[nbhd_size[i][j]])) 
     /*  if (grd[i][j]>=a) */
       {
         sk[i][j]=1;
         Iarray[0][len0]=i;
         Jarray[0][len0]=j;
         len0++;
       }
       else if (grd[i][j]>b) 
       {
         sk[i][j]=2;
         Iarray[1][len1]=i;
         Jarray[1][len1]=j;
         len1++;
       }
       else sk[i][j]=0;
     }
     
     arr_len[0]=len0;
     ave_cost(1,0);
     
     /* fill in gaps */
     /* if in a 5x5 nbhd of a point above
     middle threshold has two points of sk1,
     add it to sk2 */
     for (i=0;i<ydim;i++) for (j=0;j<xdim;j++)
     if (sk[i][j]==0)
     if (grd[i][j]>(b-tol_d[nbhd_size[i][j]]))
     {
       mk1=mk2=0;
       for (m=i-2;m<i+3;m++) for (n=j-2;n<j+3;n++)
       if (sk[m][n]==1)
       {
         mk = mark[0][m][n];
         if (mk1==0) mk1=mk;
         else if (mk!=mk1) mk2=mk;    
       }
       if (mk2>0)
       {
         sk[i][j]=2;
         Iarray[1][len1]=i;
         Jarray[1][len1]=j;
         len1++;
       }
     } 
     arr_len[1]=len1;
     ave_cost(2,1);
 }
 
 
 
   void ave_cost(int tp, int endpt_flag)
   /* numbers the branches of sk of type tp,
   calculates their average costs */
  /* lists endpts of all the branches if
  the input endpt_flag > 0 */
  {
  
     int i,j,m,n,p,q,k,tpp;
     int mm,nn,ind,len;
     int count,kk,flag,flag1,flag2;
     float a,b,c;
   
     a = upper_threshold;
     ind = tp-1;
     if (tp==1) tpp=2;
     else tpp=1;
     
     for (i=0;i<ydim;i++) for (j=0;j<xdim;j++)
     {
       tmp[i][j]=mark[ind][i][j]=0;
       if (endpt_flag>0) endpts[i][j]=0;
     }
     
     
     /* number each branch of sk, from 1 to branchcnt */
     count = 0;
     for (k=0;k<arr_len[ind];k++)
     {
       m = Iarray[ind][k];
       n = Jarray[ind][k];
       if (mark[ind][m][n]==0)
       {
         count++;
         mark[ind][m][n]=kk=count;       
         len=1;
         b = MAX(0,upper_threshold-tol_d[nbhd_size[m][n]]); 
         c = gr[m][n]/cosf(b/57.2958);     
         flag = 0;

         while (flag==0)
         {
            tmp[m][n]=1; /* do not visit it again */
            
/* if (kk==266) if (tp==1) */
/* if (kk==11){
 for (p=m-1;p<m+2;p++) for (q=n-1;q<n+2;q++)
 if (sk[p][q]>0)
 printf("at (%d %d) sk(%d %d)=%d\n",m,n,p,q,sk[p][q]);
 printf("\n");} */

   
             /* if the function "ave_cost" is called
              with endpt_flag>0, check if (m,n) is an endpt */
 
              if (endpt_flag>0)
              {
                flag2=0;            
                for (p=m-1;p<m+2;p++) for (q=n-1;q<n+2;q++)
                if (sk[p][q]==tpp) 
                   flag2=1;                            
                if (flag2>0)
                {
                  endpts[m][n]=1;
                }
              } 
            
            
            /* find contiguous unmarked points if any */
            flag1=1;
            for (p=m-1;p<m+2;p++) for (q=n-1;q<n+2;q++)
            if (sk[p][q]==tp) if (mark[ind][p][q]==0)
            {
              mark[ind][p][q] = kk;
              len++;
              b = MAX(0,upper_threshold-tol_d[nbhd_size[p][q]]); 
              c += gr[m][n]/cosf(b/57.2958); 
              mm = p; nn = q;
              flag1 = 0;
            }
            
            if (flag1==0)
            {
              m = mm; n = nn; /* next grw pt */
            }
            
            else
            /* look for an unvisited marked pt on the branch */
            {
              for (mm=0;mm<arr_len[ind];mm++)
              {
                p = Iarray[ind][mm];
                q = Jarray[ind][mm];
                if (tmp[p][q]==0) if (mark[ind][p][q]==kk)
                {
                   flag1 = 0;
                   m = p; n = q;
                   break;
                }
              }
            }

            flag = flag1;
         }
         length[ind][kk]=len;
         cost[ind][kk] = c/len;    
       }
     }
     
     branchcnt[ind] = count;
     if (endpt_flag>0) confirm_endpts(); 
  }
  
  
  
  void confirm_endpts()
  /* if the branch extends from the endpt upstream
  and downstream, it cannot be an endpt */
  {
     int m,n,p,q,mk,mm,nn;
     int flag;
     float vmax,vmin;
  
     for (m=0;m<ydim;m++) for (n=0;n<xdim;n++)
     if (endpts[m][n]>0)
     {
       mk=mark[1][m][n];
       vmin=vmax=v[m][n];
       flag = 1;
       
       for (p=m-1;p<m+2;p++) for (q=n-1;q<n+2;q++)
       if (mark[1][p][q]==mk) 
       if (v[p][q]>vmax)
       {
         vmax = v[p][q];
         mm=p; nn=q;
         flag = 0;
       }       
       if (flag==0)
       {
         flag=1;
         for (p=mm-1;p<mm+2;p++) for (q=nn-1;q<nn+2;q++)
         if (mark[1][p][q]==mk) 
         if (v[p][q]>vmax)
         {
           vmax = v[p][q];
           flag = 0;
         }
       }
       
       if (flag==0)
       {
         flag=1;
         for (p=m-1;p<m+2;p++) for (q=n-1;q<n+2;q++)
         if (mark[1][p][q]==mk) 
         if (v[p][q]<vmin)
         {
           vmin = v[p][q];
           mm=p; nn=q;
           flag = 0;
         }
         if (flag==0)
         {
           flag=1;
           for (p=mm-1;p<mm+2;p++) for (q=nn-1;q<nn+2;q++)
           if (mark[1][p][q]==mk) 
           if (v[p][q]<vmin)
           {
             vmin = v[p][q];
             flag = 0;
           }
         }
       }
    
       if (flag==0)
       /* the sk2 branch goes past all sk1 branches */
       /* check if the sk1 branches are too short and
       should be eliminated */
       {
          mk=0;
          for (p=m-1;p<m+2;p++) for (q=n-1;q<n+2;q++)
          if (sk[p][q]==1) if ((mk=mark[0][p][q])>0)
          {
            if (length[0][mk]>L) flag=1;
            else if (cost[0][mk]<1-beta/length[0][mk]) flag = 1;
              /* next to a significant sk1 branch */
          }
       }
       endpts[m][n]=flag;
    }
  }

 
 
 
 void fill_gaps()
 /* try to extend sk1 by filling gaps. 
 Tracks along ridges of grd within sk2 in order to
 find a path along the gap with minimum cost. */
 {
     int i,j,m,n,p,q,kk,mm,nn;
     int i1,i2,j1,j2;
     int mk,mk0,mk1,mk2;
     int count,len;
     int flag,terminal_flag,ignore_flag;
     int Iarr[L],Jarr[L];
     float a,b,c,d;
   
     for (i=0;i<ydim;i++) for (j=0;j<xdim;j++) tmp[i][j]=0;
     
     count=arr_len[1]=0;
     /* test at each endpt of sk2 whether we can extend */ 
     for (i=0;i<ydim;i++) for (j=0;j<xdim;j++)
     if (endpts[i][j]>0) if ((kk=mark[1][i][j])!=0)
     {
        len=flag=terminal_flag=ignore_flag=0;
        Iarr[len]=m=i;
        Jarr[len]=n=j;
        len++;
        tmp[m][n]=1;
        b = MAX(0,upper_threshold-tol_d[nbhd_size[m][n]]); 
        c = gr[m][n]/cosf(b/57.2958); 
       
        /* test if this is possibly an isolated pt of sk2 */
         mk0=mk1=mk2=0;
         for (p=m-1;p<m+2;p++) for (q=n-1;q<n+2;q++)
         if ((mk0=mark[0][p][q])!=0) /* nbhr of sk1 */
         {   
           if (mk1==0) mk1=mk0;
           else if (mk0!=mk1)
           mk2=mk0; /* near at least 2 branches of sk1 */ 
         }
         if (mk1==0)
         printf("nearby pt on sk1 not found at (%d %d)\n",m,n);
         else if (mk2>0) terminal_flag=1; /* May be an isolated sk2 pt */
         else mk0 = mk1;
         /* nbhr of a single  branch of sk1 which is marked mk0 */
         
/*if (kk==mark[1][364][294]) 
{
printf("mk0=%d mk1=%d mk2=%d",mk0,mk1,mk2); 
printf(" term=%d ign=%d\n",terminal_flag,ignore_flag);
} */       

        while ((flag+terminal_flag+ignore_flag==0)&&(len<L))
        {
/* if (kk==422)
{
 printf("k=%d m=%d n=%d \n",k,m,n);
 for (p=m-1;p<m+2;p++) for (q=n-1;q<n+2;q++)
  {
  printf("(%d %d) mk0=%d,",p,q,mark[0][p][q]);
  printf(" len=%d, sk=%d grd=%f\n",len,sk[p][q],grd[p][q]);
 }
} */
       
 /* test if the branch merges into a previosly acquired gap-filling
    piece of  sk1. (mark=0 at the gap-filling pieces of sk1 */
       for (p=m-1;p<m+2;p++) for (q=n-1;q<n+2;q++)
       if (sk[p][q]==1) if (mark[0][p][q]==0)
          ignore_flag=1; 
       
     /* test if the branch has reached another branch of sk1 */
       if (ignore_flag==0) 
       for (p=m-1;p<m+2;p++) for (q=n-1;q<n+2;q++)
       if ((mk=mark[0][p][q])!=0) if (mk!=mk0)
       /* reached another branch of sk1 */
       {
         for (i1=p-1;i1<p+2;i1++) for (j1=q-1;j1<q+2;j1++)
         if (endpts[i1][j1]>0) /* reached another endpt */
         terminal_flag=1;
       }

/* if (k==508)  if (terminal_flag>0) 
for (p=m-1;p<m+2;p++) for (q=n-1;q<n+2;q++) 
printf("mk0=%d, mk0(%d %d)=%d\n",mk0,p,q,mark[0][p][q]); */

       if (terminal_flag+ignore_flag==0) 
       /* extend the branch by descent along grd ridge */
       {
         
        /* by shallowest descent along grd from (m,n) along sk1 */
         a = 0;
         for (p=m-1;p<m+2;p++) for (q=n-1;q<n+2;q++)
         if (mark[1][p][q]==kk)
         if ((b=grd[p][q])>middle_threshold)
         if (tmp[p][q]==0)
         if (b>a)
         {
           mm=p; nn=q;
           a = b;
         }
         if (a>middle_threshold)
         {
           Iarr[len]=m=mm; Jarr[len]=n=nn;
           len++;
           tmp[m][n]=1;
           b = MAX(0,upper_threshold-tol_d[nbhd_size[p][q]]); 
           c += gr[m][n]/cosf(b/57.2958); 
/* if (kk==508) printf("%d %d len=%d\n",m,n,len); */              
         }
         else flag=1; /* cannot extend any further */
       }
/*if (kk==225)
printf("(%d %d) flag=%d a=%f len=%d\n",mm,nn,flag,a,len); */
    
       }
/* if (kk==508)
{
printf("(%d %d) len=%d ",Iarr[0],Jarr[0],len);
printf(" ignore=%d term %d\n",ignore_flag,terminal_flag);
} */         
 
      for (p=0;p<len;p++) tmp[Iarr[p]][Jarr[p]]=0;
      if (terminal_flag>0) if (ignore_flag==0)
      /* join the branch to sk1 if the cost is below threshold */
      {
        if (len==1) d=1;
        else  if (len<10) 
        {
          i1=Iarr[0]; j1=Jarr[0];
          i2=Iarr[len-1]; j2=Jarr[len-1];
          i2=i2-i1;
          if (i2<0) i2 = -i2;
          if (i2!=0) i2++;
          j2=j2-j1;
          if (j2<0) j2 = -j2;
          if (j2!=0) j2++;
          d=i2*i2+j2*j2;
          if (d<1) d=1;
          else d = sqrt((double) d);
        }
        else 
        /* divide the branch into 2 segments to measure
        the length more accurately */ 
        {
          int i3,j3;
          float y,z,d0,d1,d2,dx,dy;
          
          i1=Iarr[0]; j1=Jarr[0];
          i2=Iarr[len-1]; j2=Jarr[len-1];
          dx = j2-j1; dy=i2-i1; d0=i1*dx-j1*dy;
          d=0;
          if ((y=dx*dx+dy*dy) != 0)
          for (p=0;p<len;p++)
          {
            m=Iarr[p]; n=Jarr[p];
            z = m*dx-n*dy-d0;
            z *= z;
            z /= y;
            if (z>d)
            {
              d = z;
              i3 = m; j3 = n;
            }
          }
          i2=i2-i3;
          if (i2<0) i2 = -i2;
          if (i2!=0) i2++;
          j2=j2-j3;
          if (j2<0) j2 = -j2;
          if (j2!=0) j2++;
          d2=i2*i2+j2*j2;
          if (d2>0) d2 = sqrt((double) d2);
          
          i1=i1-i3;
          if (i1<0) i1 = -i1;
          if (i1!=0) i1++;
          j1=j1-j3;
          if (j1<0) j1 = -j1;
          if (j1!=0) j1++;
          d1=i1*i1+j1*j1;
          if (d1>0) d1 = sqrt((double) d1); ;
          
          d = d1+d2;
        }
        
        c = c/len;   
     
/* if (kk==138) printf("len=%d, d=%f c=%f c0=%f\n",len,d,c,1+beta/d); */
        
        if (c<1+beta/d)
        for (p=0;p<len;p++) 
        {       
          sk[m=Iarr[p]][n=Jarr[p]]=1;
/* if (kk==138) printf("(%d %d)\n",m,n);*/ 
          Iarray[0][arr_len[0]]=m;
          Jarray[0][arr_len[0]]=n;
          arr_len[0]++;
          mark[1][m][n]=0;
        }
      }

/* if (kk==194) printf("sk=%d \n",sk[210][391]);*/
    }
 }
 
 
 
 
 
 void complete_sk()
 /* marks endpts of branches of sk2;
 endpt = 11 at the end with highest v and
 endpt = 21 at the end with lowest v  */
 {
     int i,j,m,n,p,q,k,len;
     int ii,jj,d0;
     int mk,flag,flag1;
     float v0,vmax,a,b;