ccfeatures.c

Computing 2D Skeleton

       Iarray[0][len]=p;
       Jarray[0][len]=q;
       len++;

     flag = 0;
     while (flag==0)
     {
/* printf("i=%d j=%d\n",i,j);
 if (i==236) if (j==133)
 for (m=i-K;m<i+K+1;m++) for (n=j-K;n<j+K+1;n++)
 printf("axis=%d gr=%f\n",axis[m][n],gr[m][n]); */
       for (m=i-K;m<i+K+1;m++) for (n=j-K;n<j+K+1;n++)
       if (m==Ictr_gn) if (n==Jctr_gn)
       {
         ii = m; jj = n;
         flag = 1;
       }
       
       if (flag==0)
       {
         x = 1.001;
         for (m=i-K;m<i+K+1;m++) for (n=j-K;n<j+K+1;n++)
         if ((m-i+K)*(m-i-K)*(n-j+K)*(n-j-K)==0)
         if (axis[m][n]==1) if (gr[m][n]<=x)
         {
           x = gr[m][n];
           ii = m;
           jj = n;
         }
       }
       Iarray[0][len]=ii;
       Jarray[0][len]=jj;
       len++;
       
/* 
 printf("(%d %d) num=%d rho=%f\n",i,j,len,v[i][j]);
*/
       
       for (m=i-K;m<i+K+1;m++) for (n=j-K;n<j+K+1;n++)
       if (axis[m][n]==1)  axis[m][n]=2; /* visited */
       i = p; j = q;
       p = ii; q = jj;
     }     
     
     
     for (i=0;i<ydim;i++) for (j=0;j<xdim;j++) tmp[i][j]=0;
     for (m=0;m<len-1;m++)
     {
       i = Iarray[0][m];
       j = Jarray[0][m];
       ii = Iarray[0][m+1];
       jj = Jarray[0][m+1];
       nn = MAX(ABS(ii-i),ABS(jj-j));
       di = ii-i;
       di /= nn;
       dj = jj-j;
       dj /= nn;

       for (n=0;n<nn;n++)
       {
         ii = i+n*di;
         jj = j+n*dj+0.0001;
         for (p=ii-1;p<ii+2;p++) for (q=jj-1;q<jj+2;q++)
           tmp[p][q] = axis[p][q];
       }
     }
     for (i=0;i<ydim;i++) for (j=0;j<xdim;j++)
     {
        if (tmp[i][j]>0) axis[i][j]=1;
        else axis[i][j] = 0;
     }
  }
  
  
  
  void extend_axis_spl()
     /* extend the axis to the bdry of the splenium */
  {
     int i,j,ii,jj,ib,jb,m,n,p,q;
     float d,x,y,z,ui,uj;
     float dist[201][201];
     
     /* template */
     for (i=0;i<201;i++) for (j=0;j<201;j++)
     {
       d = (i-100)*(i-100)+(j-100)*(j-100);
       dist[i][j]=sqrt((double)d);
     }
     
     /* find a point half way to the bdry */
     i = Ictr_spl;
     j = Jctr_spl;

     x  = v[i][j]/2;
     d = 0;
     for (p=i-x-1;p<i+x+2;p++) for (q=j-x-1;q<j+x+2;q++)
     if (axis[p][q]>0)
     if ((y=dist[p-i+100][q-j+100])<=x)
     if (y>d)
     {
       d = y;
       ii=p;
       jj=q;
     }

     /* find the bdry pt on line from (ii,jj) to (i,j) */
/*  printf("nearby spl_pt: (%d %d) ctr: (%d %d) ",ii,jj,i,j); */ 
     x  = v[i][j];
     d = (i-ii)*(i-ii)+(j-jj)*(j-jj);
     d = sqrt((double)d);
     ui = (i-ii)/d;
     uj = (j-jj)/d;
     z = 0;
     m = i+ui*x;
     n = j+uj*x;
     for (p=m-10;p<m+11;p++) for (q=n-10;q<n+11;q++)
     if (bdry[p][q]>0) if (pic[p][q]<200)
     if ((y=((p-i)*ui+(q-j)*uj)/dist[p-i+100][q-j+100])>z)
     {
       z = y;
       ib=p;
       jb=q;
     }
/*  printf("spl_bdry: (%d %d)\n",ib,jb); */
    
    /* extend the axis */
    x = dist[ib-i+100][jb-j+100];
    for (m=1;m<x;m++)
    {
      p = i+ui*m;
      q = j+uj*m;
      if (bdry[p][q]>0) break;
      else if (pic[p][q]<200)
      {
        axis[p][q]=2;
        spl_i = p;
        spl_j = q;      
      }
    }
  }
  
  
  
  
  void extend_axis_gn()
  {
     int i,j,ib,jb,ih,jh,m,n,p,q;
     int len,flag,hook_flag=0;
     float d,x,y,z;
     float v0,dist[201][201];
     
     /* template */
     for (i=0;i<201;i++) for (j=0;j<201;j++)
     {
       d = (i-100)*(i-100)+(j-100)*(j-100);
       dist[i][j]=sqrt((double)d);
     }
     
     
     /* find the point on the hook axis nearest
     to the center */
    
     /* if all the points on the disk centered at (i,j)
     have sk=0, that means that the hook is almost
     nonexistent. That is the hook bdry is concentric
     with the genu disk. */
/* printf("ctr:(%d %d)\n",Ictr_gn,Jctr_gn); */    
     i = Ictr_gn;
     j = Jctr_gn; 
     x  = v[i][j];
     z = 1;
     d = 0;
     /* find the point on the disk circumference with
     the largest value of gr */
     for (p=i-x-1;p<i+x+2;p++) for (q=j-x-1;q<j+x+2;q++)
     if (axis[p][q]==0) if (q<=j)
     if ((y=dist[p-i+100][q-j+100])<=x)
     if (y>x-1.5) if (gr[p][q]<z)
     {
       flag = 0;
       for (m=p-2;m<p+3;m++) for (n=q-2;n<q+3;n++)
       if (axis[m][n]+bdry[m][n]>0) flag = 1;
       
       if (flag==0)
       {
         z = gr[p][q];
         ih=p;
         jh=q;
/* printf("(%d %d) sk=%d, grd=%f\n",p,q,sk[p][q],grd[p][q]); */
       }
     }
/* printf("gn_ctr (%d %d): nearby pt (%d %d) sk=%d\n",i,j,ih,jh,sk[ih][jh]);*/ 
     
     
     if (sk[ih][jh]==0) /* see if scan along the circumference
                      of the disk missed the skeleton 
                      by cutting along it diagonally */
       for (p=ih-1;p<ih+2;p++) for (q=jh-1;q<jh+2;q++)
         if (sk[p][q]>0) if (sk[ih][jh]!=1)
           sk[ih][jh]=sk[p][q];
     
     if (sk[ih][jh]==0) /* negligible hook; find the farthest bdry pt */
     {
       hook_flag = 1;
       z = x*x;
       for (p=i-100;p<i+101;p++) for (q=j-100;q<j+101;q++)
       if (bdry[p][q]>0) if (p>i) if (q<j)
       for (m=p-1;m<p+2;m++) for (n=q-1;n<q+2;n++)
       if (pic[m][n]<200) if (bdry[m][n]==0)
       if (grd[m][n]>lower_threshold)
       if ((y=(m-i)*(m-i)+(n-j)*(n-j))>z)
       {
         z = y;
         ib = m;
         jb = n;
       }
     }

     if (hook_flag==0)
     /* extend (ih,jh) towards the center */
     {
       m = ih; n = jh;
       axis[m][n]=2;
       v0 = v[m][n];
       x = v[Ictr_gn][Jctr_gn];
       flag=0;
       while (flag==0)
       {
         x = 0;
         for (p=m-1;p<m+2;p++) for (q=n-1;q<n+2;q++)
     /*   if ((axis[p][q]==3)||(ctrs[p][q]>0)) flag = 1; */
         if (axis[p][q]==1)
         {
           flag = 1;
           if (v[p][q]>x)
           {
             hk_i0=p; hk_j0=q;
             x = v[p][q];
           }
         }
         
               
         if (flag==0) /* grow the axis */
         {
           /* try to find a point with sk>0 */
           flag = 1;
           for (p=m-1;p<m+2;p++) for (q=n-1;q<n+2;q++)
           if (axis[p][q]==0)
           if (sk[p][q]>0) if (v[p][q]>=v0-0.0001)
           if (dist[p-Ictr_gn+100][q-Jctr_gn+100]<=x)
           {
             v0 = v[p][q];
             i = p;
             j = q;
             flag = 0;
           }
          
           if (flag>0) /* skeleton pt not found */
           {
            /* find pts on the gray skeleton, but  not on the axis
          and closer to the gn_ctr than (m,n) */
             z = (m-Ictr_gn)*(m-Ictr_gn)+(n-Jctr_gn)*(n-Jctr_gn);
             for (p=m-1;p<m+2;p++) for (q=n-1;q<n+2;q++)
             if (grd[p][q]>lower_threshold) if (axis[p][q]==0)
             if ((y=(p-Ictr_gn)*(p-Ictr_gn)+(q-Jctr_gn)*(q-Jctr_gn))<z)
               axis[p][q]=3;
            
             /* among these points, find one with highest grd */
             z = 0;
             for (p=m-1;p<m+2;p++) for (q=n-1;q<n+2;q++)
             if (axis[p][q]==3)
             {
               axis[p][q]=0;
               if (grd[p][q]>z)
               {
                 z = grd[p][q];
                 i = p;
                j = q;
                 flag = 0;
               }
             }
        
           }
         }
         if (flag==0)
         {
           m = i;
           n = j;
           axis[m][n]=2;
         } 
      }

 
       /* find the farthest bdry pt in the hook */
       /* scan from (ih,jh) outside the disk around gn_ctr */
       x = dist[ih-Ictr_gn+100][jh-Jctr_gn+100];
       jb = jh;
       len = 0;
       Iarray[0][len]=ih;
       Jarray[0][len]=jh;
       len++;

       while (len>0)
       {
         m = Iarray[0][0];
         n = Jarray[0][0];

         for (p=m-1;p<m+2;p++) for (q=n-1;q<n+2;q++)
         if (grd[p][q]>lower_threshold) if (bdry[p][q]==0)
         if (pic[p][q]<200) if (axis[p][q]==0)
         if (dist[p-Ictr_gn+100][q-Jctr_gn+100]>x-1)
         {
         /*  flag=0;
           i=p-Ictr_gn+100;
           j=q-Jctr_gn+100;
           if (i>0) if (i<201) if (j>0) if (j<201)
           if (dist[i][j]<=x-1) flag = 1; 
            
           if (flag==0) */
           Iarray[0][len]=p;
           Jarray[0][len]=q;
           axis[p][q]=3; /* visited the point */
           len++;

           if (q<jb)
           {
             ib=p;
             jb=q;
           }
           else if (q==jb) if (p>ib) ib=p;
         }
         len--;
         if (len>0)
         {
           Iarray[0][0]=Iarray[0][len];
           Jarray[0][0]=Jarray[0][len];
         }
       }
       for (i=0;i<ydim;i++) for (j=0;j<xdim;j++)
          if (axis[i][j]==3) axis[i][j]=0; /* reset visit mark */
     }
/* printf("gn_bd:(%d %d) hook_flag=%d\n",ib,jb,hook_flag); */

    
     /* track back from (ib,jb) towards genu center */
     m = ib; n = jb;
     axis[m][n]=2;
     v0 = v[m][n]; 
     flag=0;
     while (flag==0)
     {
        x =0;
        for (p=m-1;p<m+2;p++) for (q=n-1;q<n+2;q++)
        {
          if (hook_flag==0)
          {
            if (ih==p) if (jh==q) flag = 1;
            if (axis[p][q]==1)
            {
              flag =1;
 printf("*** hook axis does not connect to the main axis extension ***\n");
            }
          }
          else if (axis[p][q]==1) /* no hook and on sk1,find max v */
          {
            flag = 1;
            if (v[p][q]>x)
            {
              hk_i0=p; hk_j0=q;
              x = v[p][q];
            }
          }
        }
        
              
        if (flag==0) /* can grow the axis towards gn_ctr */
        {
          {
            /* try to find a point with increased value of v */
            flag = 1;
            for (p=m-1;p<m+2;p++) for (q=n-1;q<n+2;q++)
            if (grd[p][q]>lower_threshold)
            if (axis[p][q]==0) if (v[p][q]>=v0-0.0001)
            {
              v0 = v[p][q];
              i = p;
              j = q;
              flag = 0;
            }
          }

          if (flag>0) /* no such point found; can happen near bdry */ 
          {
          /* find pts on the gray skeleton, but  not on the axis
          and closer to the gn_ctr than (m,n) */
            z = (m-Ictr_gn)*(m-Ictr_gn)+(n-Jctr_gn)*(n-Jctr_gn);
            for (p=m-1;p<m+2;p++) for (q=n-1;q<n+2;q++)
            if (grd[p][q]>lower_threshold) if (axis[p][q]==0)
            if ((y=(p-Ictr_gn)*(p-Ictr_gn)+(q-Jctr_gn)*(q-Jctr_gn))<z)
            axis[p][q]=3;
          
          /* among these points, find the one with the
          highest value of grd */
            z = 0;
            for (p=m-1;p<m+2;p++) for (q=n-1;q<n+2;q++)
            if (axis[p][q]==3)
            {
              axis[p][q]=0;
              if (grd[p][q]>z)
              {
                z = grd[p][q];
                i = p;
                j = q;
                flag = 0;
              }
            }
          }
/* if (m==243) if (n==531)
 {
    printf("flag=%d\n",flag);
    z = (m-Ictr_gn)*(m-Ictr_gn)+(n-Jctr_gn)*(n-Jctr_gn);
    for (p=m-1;p<m+2;p++) for (q=n-1;q<n+2;q++)
    if (grd[p][q]>lower_threshold) if (axis[p][q]==0)
          if ((y=(p-Ictr_gn)*(p-Ictr_gn)+(q-Jctr_gn)*(q-Jctr_gn))<z)
    {
      printf("y=%f z=%f grd=%f\n",y,z,grd[p][q]);
    }
 } */
 
          if (flag==0)
          {
            m = i;
            n = j;
            axis[m][n]=2;
/* printf("%d %d\n",m,n); */
          } 
       }
    }

    if (hook_flag>0)
    /* find intersection of axis with disk at genu_ctr */
    {
       i = Ictr_gn;
       j = Jctr_gn;
       x  = v[i][j];
       z = dist[ib-i+100][jb-j+100];
       if (z>x)
       for (p=i+1;p<i+x+2;p++) for (q=j-x-1;q<j;q++)
       {
         if (axis[p][q]==2)
         if ((y=dist[p-i+100][q-j+100])<=x)
         if (y>x-1.5)
         {
           ih=p;
           jh=q;
         }
       }
       else
       {
         ih = ib;
         jh = jb;
       }
/* printf("hook_flag=%d; new ih,jh (%d %d)\n",hook_flag,ih,jh); */
    }
     
     gn_i = ih;
     gn_j = jh;
     
     hk_i = ib;
     hk_j = jb;
  }
  
  
 
  
  void calc_distance()
  /* calculates arc-distance at points on the axes
  (K/2-1) pixels 'apart' */
  /* each iteration finds a point distance K away */
  /* next iteration is centered at the previously
  found point at distance K/2 */
  {
     int i,j,m,n,ii,jj,p,q;
     int len,flag;
     float x;
     int K=8; /* K must be even. K=8 is probably the minimum value
                needed for numerical  stability */     
    
     
     /* calculate arc distance along the main axis */
     initial_i = Ictr_spl;
     initial_j = Jctr_spl;
     final_i = Ictr_gn;
     final_j = Jctr_gn; 
   
     for (i=0;i<ydim;i++) for (j=0;j<xdim;j++) tmp[i][j]=0;
        
     i = initial_i; j=initial_j;
     len=0;
     Iarray[0][len]=i;
     Jarray[0][len]=j;
     axis_distance[len]=0;
     len++;
 
       x = 1; 
       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;
       }
       x = (p-i)*(p-i)+(q-j)*(q-j);