ccfeatures.c

Computing 2D Skeleton

/* Program:  extracts features                 */

 /*


/* Assumes that image is piecewise constant                 */

/* Assumes that the image boundary is at least (N-1)/2 pixels
   away from the frame boundary                             */


/*                                                          */
/*                                                          */
/* Inputs:  the image and the ABSOLUTE distance function    */
/*                                                          */
/*                                                          */
/* Date:     June, 2003                                    */
/*                                                          */
/*                                                          */




#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>


/* imagesizes:
   ccSchzSchz:  376H x 796W; sample count: 12+23=35
   ccSchzNorms: 364H x 784W; sample count: 21+27=48
   ccNorms: 324H x 780W; sample count: 10+10=20 
   ccKids : 400H x 896W; sample count: 15+15=30
   ccDrugs: 324H x 708W; sample count: 7+7=14
 */
  
 
/* ccSchzSchz */
/*
#define caseno 1 
#define grp_1_count 12
#define grp_2_count 23
#define sample_count 35 
#define ydim 376  
#define xdim 796
*/ 

/* ccSchzNorms */
/*
#define caseno 2
#define grp_1_count 21
#define grp_2_count 27
#define sample_count 48 
#define ydim 364  
#define xdim 784 
*/

/* ccNorms */
/*
#define caseno 3
#define grp_1_count 10
#define grp_2_count 10
#define sample_count 20 
#define ydim 324  
#define xdim 780
*/

/* ccKids */
/*
#define caseno 4
#define grp_1_count 15
#define grp_2_count 15
#define sample_count 30 
#define ydim 400  
#define xdim 896 
*/
 
/* ccDrugs */

#define caseno 5
#define grp_1_count 7
#define grp_2_count 7
#define sample_count 14 
#define ydim 324  
#define xdim 708



#define num_axis_samples 72
#define num_spl_samples 6  /* splenium */
#define num_hk_samples 34 /* genu hook */

#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 NN 40 /* size of local nbhd for finding ctrs  */


#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];
char dir_name[128];

/* Program variables */
unsigned char pic[ydim][xdim];
unsigned char outpic[ydim][xdim];
unsigned char outgr[ydim][xdim];
unsigned char sk[ydim][xdim];
unsigned char ctrs[ydim][xdim];
unsigned char axis[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];
int length[2][xdim*ydim];
float v[ydim][xdim];
float gr[ydim][xdim],grd[ydim][xdim]; /*deriv.of v along skeleton*/
float tol[(N-1)/2],tol_d[(N-1)/2]; /* margin of error in angles */
float cost[2][xdim*ydim];


float axis_distance[1000],spl_distance[100],hk_distance[250];
int Iarr[100],Jarr[100];
int sample_index;
int num_axis_pts; /* for calculating  arc_distances along the axis */
int num_spl_pts;
int num_hk_pts;
float axis_length[sample_count];
float axis_rho[sample_count][num_axis_samples];
float spl_rho[sample_count][num_spl_samples];
float hk_rho[sample_count][num_hk_samples];
float rel_rbn[sample_count],rel_gn[sample_count],rel_spl[sample_count];


int Ictr_spl,Jctr_spl,Ictr_gn,Jctr_gn;
int gn_i,gn_j; /* on the bdry of the disk in the genu */
int spl_i,spl_j; /* on the bdry  of the disk in the splenium */
int hk_i0, hk_j0; /* beginning of the axis of the genu hook  */ 
int hk_i,hk_j; /* the end of genu hook  */
int initial_i,initial_j,final_i,final_j;  

int ibb,inn,jbb,jnn;
int debug_flag;

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 calc_features();
void find_axis();
void clean_up();
void extend_axis_spl();
void extend_axis_gn();
void calc_distance();
void sample_rho();
void calc_CCareas();
void output_pics();
void output_features();
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; /*frame bdry*/
        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++)