ccfeatures.c

Computing 2D Skeleton

     a = upper_threshold;
     b = lower_threshold;
     
     len=0;
     for (i=0;i<ydim;i++) for (j=0;j<xdim;j++)
     {
       if (sk[i][j]==1)
       {
         Iarray[0][len]=i;
         Jarray[0][len]=j;
         len++;
       }
       else sk[i][j]=0;
     }
       
     arr_len[0]=len;
     ave_cost(1,0);  
    
     
     /* now extend each branch downstream (increasing v)
         when it meets another branch or cannot be extended
         any further, mark the point by setting tmp=1; 
         If the criteria for setting tmp=1 are not met,
         the branch is eliminated in the next step when
          we ascend  */   
     
     for (i=0;i<ydim;i++) for (j=0;j<xdim;j++)
       tmp[i][j]=endpts[i][j]=0;
       

     for (k=0;k<len;k++)
     /* temporary storage */
     {
       Iarray[1][k]=Iarray[0][k];
       Jarray[1][k]=Jarray[0][k];
     }

     while (len>0) 
     {
          m = Iarray[1][0];
          n = Jarray[1][0];
          mk = mark[0][m][n];
          v0 = v[m][n];
          
 /* if (mk==123) */
/*  if (m==362) if (n==7)  
 {
  printf("%d %d v=%f:\n",m,n,v[m][n]);
  for (i=m-1;i<m+2;i++)
  {
  for (j=n-1;j<n+2;j++) 
  printf("(%d %d),sk=%d mk=%d ",i,j,sk[i][j],mark[0][i][j]);
  printf("\n");
  }
}  */       
          flag1 = flag = 0;
          if (sk[m][n]==2)
          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]==mk)
            flag1=1; /* still touching the original branch */
          
          if (sk[m][n]==2)
          for (p=m-1;p<m+2;p++) for (q=n-1;q<n+2;q++)
          if (sk[p][q]>0) if (mark[0][p][q]!=mk)
          /* meets another branch */
          {
            if (sk[p][q]==1) flag=1;
            else if (flag1==0) flag=1;
          }
          
          if (flag>0) tmp[m][n]=1;
          else  /* try to extend the branch */
          {
            vmax = 0;
            for (p=m-1;p<m+2;p++) for (q=n-1;q<n+2;q++)
           /* if (sk[p][q]==0) *//* not on current skeleton */
            if (grd[p][q]>b) if (v[p][q]>vmax)
            {
              vmax = v[p][q];
              ii=p; jj=q;
            }

            if (vmax>0) if (vmax<v0)
            {
               d0 = (ii-m)*(ii-m)+(jj-n)*(jj-n);
               for (p=ii-1;p<ii+2;p++) for (q=jj-1;q<jj+2;q++)
               if ((p-m)*(p-m)+(q-n)*(q-n) - d0 >2)
               /* picks out correct points in the 5x5 nbhd to compare */
             /*  if (sk[p][q]==0) *//* not on current skeleton */
               if (grd[p][q]>b) if (v[p][q]>vmax)
                 vmax = v[p][q];
            }
          
            /* add to the skeleton if (ii,jj) is downstream */
            if (vmax>=v0) if (sk[ii][jj]==0)
            {
              sk[ii][jj] = 2;
              mark[0][ii][jj] = mk;
              Iarray[1][len]=ii; /* insert it in the array */
              Jarray[1][len]=jj;
              len++;
/* if (mk==201) printf("%d %d\n",ii,jj);*/                        
            }
            
            else /* has the branch reached lower_threhshold? */
            if (flag1==0)
            {
              vmax = 0;
              for (p=m-1;p<m+2;p++) for (q=n-1;q<n+2;q++)
              if (sk[p][q]==0) /* not on current skeleton */
              if (v[p][q]>vmax)
              {
                vmax = v[p][q];
                ii=p; jj=q;
              }

              if (vmax>0) if (vmax<v0)
              {
                 d0 = (ii-m)*(ii-m)+(jj-n)*(jj-n);
                 for (p=ii-1;p<ii+2;p++) for (q=jj-1;q<jj+2;q++)
                 if ((p-m)*(p-m)+(q-n)*(q-n) - d0 >2)
                 /* picks out correct points in the 5x5 nbhd to compare */
                 if (sk[p][q]==0) /* not on current skeleton */
                 if (v[p][q]>vmax)
                   vmax = v[p][q];
              }
              if (vmax>=v0) tmp[m][n]=1;
            }
          }

 /* if (m==362) if (n==7) 
  for (i=m-2;i<m+3;i++)
  {
    for (j=n-2;j<n+3;j++) 
    printf("%f ",v[i][j]); printf("\n"); 
  } 
  if (m==362) if (n==7) 
  for (i=m-2;i<m+3;i++)
  {
    for (j=n-2;j<n+3;j++) 
    printf("%f ",grd[i][j]); printf("\n"); 
  } */

/* if (m==362) if (n==7)  
{
  for (i=m-2;i<m+3;i++) for (j=n-2;j<n+3;j++)  
 { printf("(%d %d),sk=%d mk=%d tmp=%d",i,j,sk[i][j],mark[0][i][j],tmp[i][j]);
  printf(" v=%f grd=%f\n",v[i][j],grd[i][j]);}
} */      
      
          /* delete current pt from the array */
          len--;
          if (len>0)
          {
            Iarray[1][0]=Iarray[1][len];
            Jarray[1][0]=Jarray[1][len];
          }
       }
      
/* printf("mark[0]=%d\n",mark[0][212][222]);
for (i=209;i<214;i++) {for (j=211;j<216;j++)
printf("%d ",mark[0][i][j]);printf("\n");} */
       
       
      /*  now ascend from pts with tmp>0 towards the original
       skeleton and mark its endpts.
        this step eliminates extraneous pts of sk2 */
     
       /* initial list of pts on sk2 */
       len = 0;
       for (i=0;i<ydim;i++) for (j=0;j<xdim;j++)
       if (tmp[i][j]>0) if (sk[i][j]==2)
       {
         Iarray[1][len]=i;
         Jarray[1][len]=j;
         len++;
         endpts[i][j]=11;
/* if (mark[0][i][j]==127) printf("%d %d\n",i,j); */
       }        
      
      /* extend the list by ascending along branches */
      /* mark the endpts when sk1 is reached */
       while (len>0)
       {
          m = Iarray[1][0];
          n = Jarray[1][0];
          mk = mark[0][m][n];
          flag=0;

/* if (m==325) if (n==63) */
/* if (mk==201) 
{
  printf("%d %d :\n",m,n);
  for (i=m-1;i<m+2;i++)
  {
  for (j=n-1;j<n+2;j++) 
  printf("(%d %d),sk=%d mk=%d ",i,j,sk[i][j],mark[0][i][j]);
  printf("\n");
  }
} */
  
            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]==mk) if (flag==0)
            {
         /* sk2 branch terminates on meeting sk1 for the first time;
            sk2 branch has higher values of v than sk1 branch */
              endpts[p][q] = 21;
              flag=1;
/* if (mk==201) printf("(%d %d) at (%d %d)\n\n",p,q,m,n); */            
            }
            
            if (flag==0)
            for (p=m-1;p<m+2;p++) for (q=n-1;q<n+2;q++)
            if (mark[0][p][q]==mk) if (tmp[p][q]==0)
            {
              tmp[p][q]=1;
              Iarray[1][len]=p;
              Jarray[1][len]=q;
              len++;
/* if (mk==122) printf("(%d %d) from (%d %d)\n",p,q,m,n); */ 
            }
      
          len--;
          if (len>0)
          {
            Iarray[1][0]=Iarray[1][len];
            Jarray[1][0]=Jarray[1][len];
          }
       }      
       
       for (i=0;i<ydim;i++) for (j=0;j<xdim;j++)
       if (sk[i][j]==2)
       {
         mark[0][i][j]=0;
         if (tmp[i][j]==0) sk[i][j]=0;
       }
      
       for (i=0;i<ydim;i++) for (j=0;j<xdim;j++)
       if (sk[i][j]==2) if ((k=endpts[i][j])>0)
       for (m=i-1;m<i+2;m++) for (n=j-1;n<j+2;n++)
       {
         if (sk[m][n]==1) endpts[m][n]=k;
         else endpts[m][n]=0;
       }
 }
 
 
  
  
  
  void prune()
 /* eliminate the branch of sk1 if its cost
     is beyond the cost threshold */
  {
     int k,m,n;
     int mk,len,len0;
     int i1,j1,i2,j2,kk;
     int Iarr[L],Jarr[L];
     float b,d;
     
     for (mk=1;mk<branchcnt[0]+1;mk++)  
     if (length[0][mk]<=L) if (length[0][mk]>0)
     { 

        len0=length[0][mk];
        len=kk=0;
        for (k=0;k<arr_len[0];k++)
        if (mark[0][m=Iarray[0][k]][n=Jarray[0][k]]==mk)
        {
           Iarr[len]=m;
           Jarr[len]=n;
           len++;
           if (endpts[m][n]>0)
           { i1=m; j1=n; kk=1; }
           if (len>=len0) break;
        }
    
        if (kk>0)
        {
          d=0;
          for (k=0;k<len;k++)
          {
            i2=Iarr[k]; j2=Jarr[k];
            i2=i2-i1;
            if (i2<0) i2 = -i2;
            if (i2!=0) i2++;
            j2=j2-j1;
            if (j2<0) j2 = -j2;
            if (j2!=0) j2++;
            b=i2*i2+j2*j2;
            if (b>d) d = b;
          }
          if (d<1) d=1;
          else d = sqrt((double) d);
        }
        else
        {
/*for (k=0;k<len;k++) tmp[Iarr[k]][Jarr[k]]=1;       
          printf("cannot find endpts: ");
          printf("pt (%d %d) mark=%d\n",m,n,mk);
for (k=0;k<len;k++) printf(" (%d %d)\n",Iarr[k],Jarr[k]);*/
          
          d=0;
          for (k=0;k<len-1;k++)
          {
            i1=Iarr[k]; j1=Jarr[k];
            for (kk=k+1;kk<len;kk++)
            {
              i2=Iarr[kk]; j2=Jarr[kk];
              i2=i2-i1;
              if (i2<0) i2 = -i2;
              if (i2!=0) i2++;
              j2=j2-j1;
              if (j2<0) j2 = -j2;
              if (j2!=0) j2++;
              b=i2*i2+j2*j2;
              if (b>d) d = b;
            }
          }
          if (d<1) d=1;
          else d = sqrt((double) d);
       }
      
       
      if (cost[0][mk]>1-beta/d)
      {
       /*
        kk=7;
        for (k=0;k<len;k++)
        if ((m=nbhd_size[Iarr[k]][Jarr[k]])<kk) kk=m;
 
        if (kk>1) */
        
        {
          length[0][mk]=0;
          for (k=0;k<len;k++) sk[Iarr[k]][Jarr[k]]=2;
        }
      }
    }
  }
  
  
  
  
  void calc_features()
  {
     find_axis();
     calc_distance(); 
     sample_rho(); 
     calc_CCareas();  
  }
  
  
  
  
  void find_axis()
  /* also finds the centers of genu and the splenium */
  /* axis marked 1 between genu ctr and splenium ctr */
  /* extensions of the axis into the genu and the splenium
  are marked as 2  */
  /* By definition, splenium and genu are disks centered
  at their respective centers plus the end regions beyond
  the disks */
  {
     int i,j,m,n,k1,k2,jmin,jmax;
     int p,q,len,flag;
     float x;
     
     for (i=0;i<ydim;i++) for (j=0;j<xdim;j++)
     axis[i][j]=ctrs[i][j]=0;
     
     /* find ctrs of genu and splenium */
    
     jmin = xdim; jmax = 0;
     for (i=0;i<ydim;i++) for (j=0;j<xdim;j++) 
     if (pic[i][j]<200)
     {
       if (j<jmin) jmin=j;
       if (j>jmax) jmax=j;
     }
     /* splenium */
     k1 = jmin+(jmax-jmin)/4;
     x = 0;
     for (i=0;i<ydim;i++) for (j=jmin;j<k1;j++)
     if (pic[i][j]<200)
     if (sk[i][j]>0) if (v[i][j]>x)
     {
       x = v[i][j];
       m = i; n = j;
     }
     ctrs[m][n]=1;
     axis[m][n]=1;
     Ictr_spl=m; Jctr_spl=n;
     
     /* genu */
     k2 = jmin+2*(jmax-jmin)/3;
     x = 0;
     for (i=0;i<ydim;i++) for (j=k2;j<jmax;j++)
     if (pic[i][j]<200)
     if (sk[i][j]>0) if (v[i][j]>x)
     {
       x = v[i][j];
       m = i; n = j;
     }
     ctrs[m][n]=1;
     axis[m][n]=1;
     Ictr_gn=m; Jctr_gn=n;
     
     
     /* find axis in middle portion */
     /* store points in I and J arrays */
     k2 = jmin+(jmax-jmin)/2;
     len = 0;
     for (i=0;i<ydim;i++) for (j=k1;j<k2;j++)
     if (sk[i][j]>0) if (pic[i][j]<200)
     {
        axis[i][j]=1;
        if ((j<k1+2)||(j>k2-2))
        {
          Iarray[0][len]=i;
          Jarray[0][len]=j;
          len++;
        }
     }
     
     /* extend the axis to ctrs */
     while (len>0)
     {
        m = Iarray[0][0];
        n = Jarray[0][0];
        flag = 0;
        
        for (p=m-1;p<m+2;p++) for (q=n-1;q<n+2;q++)
        if (ctrs[p][q]>0) flag = 1;
        
        if (n<k1) if (m>Ictr_spl) flag = 1;
        if (n>k2) if (m>Ictr_gn) flag = 1;
 
        if (flag==0)
        for (p=m-1;p<m+2;p++) for (q=n-1;q<n+2;q++)
        if (sk[p][q]>0) if (axis[p][q]==0)
        {
          axis[p][q]=1;
          Iarray[0][len]=p;
          Jarray[0][len]=q;
          len++;
        }
        len--;
        if (len>0)
        {
          Iarray[0][0]=Iarray[0][len];
          Jarray[0][0]=Jarray[0][len];
        }
     }
     
     /* extend the axis further into splenium */
      flag = 0;  
      m = Ictr_spl;
      n = Jctr_spl;
      while (flag==0)
      {  
        flag = 1;
        for (p=m;p<m+2;p++) for (q=n-1;q<n+2;q++)
        if (sk[p][q]==1) if (axis[p][q]==0)
        {
          axis[p][q]=1;
          m = p;
          n = q;
          flag = 0;
          break;
        }
      }
      
      ctrs[Ictr_spl][Jctr_spl]=0;
      ctrs[m][n]=1;
      Ictr_spl = m;
      Jctr_spl = n;
      
      clean_up();
 
      extend_axis_spl(); 
      extend_axis_gn(); 

  }
  
  
   
  
   void clean_up()
  /* tracks the main axis by marking points on the axis
  (K/2-1) pixels 'apart' */
  {
     int i,j,m,n,nn,ii,jj,p,q;
     int len,flag;
     float di,dj,x;
     int K=8; /* K must be even. If K is too small, tracking
              may pick up protrusion branches; probably 8
              is the minimum value */

     i = Ictr_spl; j=Jctr_spl;
     len=0;
     Iarray[0][len]=i;
     Jarray[0][len]=j;
     len++;
 
       x = 1.001; 
       for (m=i-K/2;m<i+K/2+1;m++) for (n=j-K/2;n<j+K/2+1;n++)
       if ((m-i+K/2)*(m-i-K/2)*(n-j+K/2)*(n-j-K/2)==0)
       if (axis[m][n]==1) if (gr[m][n]<x)
       {
         x = gr[m][n];
         p = m;
         q = n;
       }