ccaxis.c

Computing 2D Skeleton

     /* extend the axis to ctrs */
     ctrs[Ictr1][Jctr1]=ctrs[Ictr2][Jctr2]=1;
     while (len>0)
     {
        m = Iarray[1][0];
        n = Jarray[1][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>Ictr1) flag = 1;
        if (n>k2) if (m>Ictr2) 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[1][len]=p;
          Jarray[1][len]=q;
          len++;
        }
        len--;
        if (len>0)
        {
          Iarray[1][0]=Iarray[1][len];
          Jarray[1][0]=Jarray[1][len];
        }
     }
    
     /* extend the axis further into splenium */
      flag = 0;  
      m = Ictr1;
      n = Jctr1;
      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[Ictr1][Jctr1]=0;
      ctrs[m][n]=1;
      Ictr1 = m;
      Jctr1 = 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 = Ictr1; j=Jctr1;
     len=0;
     Iarray[1][len]=i;
     Jarray[1][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;
       }
       Iarray[1][len]=p;
       Jarray[1][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==Ictr2) if (n==Jctr2)
       {
         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[1][len]=ii;
       Jarray[1][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[1][m];
       j = Jarray[1][m];
       ii = Iarray[1][m+1];
       jj = Jarray[1][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()
  {
     int i,j,ib,jb,m,n,p,q;
     float d,x,y,z,ui,uj;
     float dist[101][101];
     
     /* template */
     for (i=0;i<101;i++) for (j=0;j<101;j++)
     {
       d = (i-50)*(i-50)+(j-50)*(j-50);
       dist[i][j]=sqrt((double)d);
     }
     
     /* extend the axis to the bdry of the splenium */
     
     i = Ictr1;
     j = Jctr1;
     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+50][q-j+50])<=x)
     if (y>d)
     {
       d = y;
       ib=p;
       jb=q;
     }

     /* find the bdry pt on line from (ib,jb) to (i,j)*/
/* printf("b: (%d %d) ctr: (%d %d) ",ib,jb,i,j); */
     x  = v[i][j];
     d = (i-ib)*(i-ib)+(j-jb)*(j-jb);
     d = sqrt((double)d);
     ui = (i-ib)/d;
     uj = (j-jb)/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+50][q-j+50])>z)
     {
       z = y;
       ib=p;
       jb=q;
     }
/*  printf("bdry: (%d %d)\n",ib,jb); */
    /* extend the axis */
    x = dist[ib-i+50][jb-j+50];
    for (m=1;m<x;m++)
    {
      p = i+ui*m;
      q = j+uj*m;
      if (bdry[p][q]>0) break;
      else
      {
        axis[p][q]=2;
        initial_i = p;
        initial_j = q;
      }
    }
/* printf("axis0 %d %d\n",axis_i0,axis_j0); */
  }
  
  
  
  
  
  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 disk. */
     
     i = Ictr2;
     j = Jctr2;
     x  = v[i][j];
     z = 1;
     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 (q<=j)
     if ((y=dist[p-i+100][q-j+100])<=x)
     if (y>x-1) 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("h: (%d %d) sk=%d\n",ih,jh,sk[ih][jh]); */
     
     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[Ictr2][Jctr2];
     flag=0;
     while (flag==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 (flag==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-Ictr2+100][q-Jctr2+100]<=x)
          {
            v0 = v[p][q];
            i = p;
            j = q;
            flag = 0;
          }
          
          if (flag>0)
          {
            z = (m-Ictr2)*(m-Ictr2)+(n-Jctr2)*(n-Jctr2);
            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-Ictr2)*(p-Ictr2)+(q-Jctr2)*(q-Jctr2))<z)
              axis[p][q]=3;
            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 */
     x = v[Ictr2][Jctr2];
     jb = jh;
     len = 0;
     Iarray[1][len]=ih;
     Jarray[1][len]=jh;
     len++;

     while (len>0)
     {
        m = Iarray[1][0];
        n = Jarray[1][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-Ictr2+100][q-Jctr2+100]>x-1)
        {
          flag=0;
          i=p-Ictr2+100;
          j=q-Jctr2+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[1][len]=p;
          Jarray[1][len]=q;
          axis[p][q]=3;
          len++;
          if (q<jb)
          {
            ib=p;
            jb=q;
          }
          else if (q==jb) if (p>ib) ib=p;
        }
        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 (axis[i][j]==3) axis[i][j]=0;
     }
  
/* printf("b:(%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)
     {
        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 (hook_flag>0) if (ctrs[p][q]>0) flag = 1;
        }      
        if (flag==0)
        {
         {
          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)
        {
          z = (m-Ictr2)*(m-Ictr2)+(n-Jctr2)*(n-Jctr2);
          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-Ictr2)*(p-Ictr2)+(q-Jctr2)*(q-Jctr2))<z)
          axis[p][q]=3;
          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;
            }