ccfeatures.c

Computing 2D Skeleton

       if (x!=0) x = sqrt((double)x);
       axis_distance[len]=x;
       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==final_i) if (n==final_j)
       {
         ii = m; jj = n;
         flag = 1;
       }
       
       if (flag==0)
       {
         x = 1;
         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 (tmp[m][n]==0) if (gr[m][n]<=x)
         {
           x = gr[m][n];
           ii = m;
           jj = n;
         }
       }
       
       x = (ii-p)*(ii-p)+(jj-q)*(jj-q);
       if (x!=0) x = sqrt((double)x);
       axis_distance[len]=axis_distance[len-1]+x;
       
       Iarray[0][len]=ii;
       Jarray[0][len]=jj;
       len++;      
 /*
 printf("(%d %d) num=%d len=%f rho=%f\n",i,j,len,axis_distance[len-1],v[i][j]);
 */

/* printf("(%d %d) len=%f\n",ii,jj,axis_distance[len-1]); */ 
   
       for (m=i-K;m<i+K+1;m++) for (n=j-K;n<j+K+1;n++)
       if (axis[m][n]==1)  tmp[m][n]=1; /* visited */
       i = p; j = q;
       p = ii; q = jj;
     }
 /*    
 x = axis_distance[len-1]/axis_length[sample_index];
  printf("axis_length=%f  ",axis_length[sample_index]);
 printf("total len =%f  ratio=%f\n",axis_distance[len-1],x); 
 */   
     num_axis_pts = len;
     axis_length[sample_index] = axis_distance[len-1];
    
     
     /* calculate arc distance along the hook axis */
     initial_i = hk_i0;
     initial_j = hk_j0;
     final_i = hk_i;
     final_j = hk_j; 
   
     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[1][len]=i;
     Jarray[1][len]=j;
     hk_distance[len]=0;
     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]==2) if (gr[m][n]<x)
       {
         x = gr[m][n];
         p = m;
         q = n;
       }
       x = (p-i)*(p-i)+(q-j)*(q-j);
       if (x!=0) x = sqrt((double)x);
       hk_distance[len]=x;
       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==final_i) if (n==final_j)
       {
         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]==2) if (tmp[m][n]==0) if (gr[m][n]<=x)
         {
           x = gr[m][n];
           ii = m;
           jj = n;
         }
       }
       
       x = (ii-p)*(ii-p)+(jj-q)*(jj-q);
       if (x!=0) x = sqrt((double)x);
       hk_distance[len]=hk_distance[len-1]+x;
       
       Iarray[1][len]=ii;
       Jarray[1][len]=jj;
       len++;      
 /*
 printf("(%d %d) num=%d len=%f rho=%f\n",i,j,len,hk_distance[len-1],v[i][j]);
 */

/* printf("(%d %d) len=%f\n",ii,jj,hk_distance[len-1]); */ 
   
       for (m=i-K;m<i+K+1;m++) for (n=j-K;n<j+K+1;n++)
       if (axis[m][n]==2)  tmp[m][n]=1; /* visited */
       i = p; j = q;
       p = ii; q = jj;
     }
 /*    
 x = hk_distance[len-1]/axis_length[sample_index];
  printf("axis_length=%f  ",axis_length[sample_index]);
 printf("total len =%f  ratio=%f\n",hk_distance[len-1],x); 
 */   
     num_hk_pts = len;

 
 
 /* calculate spl distance along the splenium axis */
     initial_i = Ictr_spl;
     initial_j = Jctr_spl;
     final_i = spl_i;
     final_j = spl_j; 
 
     for (i=0;i<ydim;i++) for (j=0;j<xdim;j++) tmp[i][j]=0;
        
     i = initial_i; j=initial_j;
     len=0;
     Iarr[len]=i;
     Jarr[len]=j;
     spl_distance[len]=0;
     len++;
 
       x = 1.001; /*  note: points on the extended axis
                  in the splenium may not have gr<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]==2) if (gr[m][n]<x)
       {
         x = gr[m][n];
         p = m;
         q = n;
       }
       x = (p-i)*(p-i)+(q-j)*(q-j);
       if (x!=0) x = sqrt((double)x);
       spl_distance[len]=x;
       Iarr[len]=p;
       Jarr[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==final_i) if (n==final_j)
       {
         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]==2) if (tmp[m][n]==0) if (gr[m][n]<=x)
         {
           x = gr[m][n];
           ii = m;
           jj = n;
         }
       }
       
       x = (ii-p)*(ii-p)+(jj-q)*(jj-q);
       if (x!=0) x = sqrt((double)x);
       spl_distance[len]=spl_distance[len-1]+x;
       
       Iarr[len]=ii;
       Jarr[len]=jj;
       len++;      
 /*
 printf("(%d %d) num=%d len=%f rho=%f\n",i,j,len,spl_distance[len-1],v[i][j]);
 */

/* printf("(%d %d) len=%f\n",ii,jj,spl_distance[len-1]); */ 
   
       for (m=i-K;m<i+K+1;m++) for (n=j-K;n<j+K+1;n++)
       if (axis[m][n]==2)  tmp[m][n]=1; /* visited */
       i = p; j = q;
       p = ii; q = jj;
     }
 /*    
 x = spl_distance[len-1]/axis_length[sample_index];
  printf("axis_length=%f  ",axis_length[sample_index]);
 printf("total len =%f  ratio=%f\n",spl_distance[len-1],x); 
 */
 
     num_spl_pts = len;    
  }
    
  
  
  


 void sample_rho()
  /* calculates aveage relative rho in each segment  */
  {
     int m,n,i,j,i1,j1,i2,j2,k1,k2;
     int p,q,actual_num_samples;
     float d,d1,d2,len,axis_len,step;
     float dist[201][201],x;
     
     /* 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);
     }  
     
     for (m=0;m<num_axis_samples;m++) axis_rho[sample_index][m]=0;
     for (m=0;m<num_hk_samples;m++) hk_rho[sample_index][m]=0;
     for (m=0;m<num_spl_samples;m++) spl_rho[sample_index][m]=0;
     
     axis_len = axis_distance[num_axis_pts-1];
     
     /* sample the main axis */
     len = axis_distance[num_axis_pts-1];
     step = len/(num_axis_samples - 1);

     k2 = 0;
     
     i = Ictr_spl; j = Jctr_spl;
     x = 0;
     n=0;
     for (p=i-step/2;p<i+step/2+1;p++)
     for (q=j-step/2;q<j+step/2+1;q++)
     if (pic[p][q]<200) if (axis[p][q]==1)
     {
       x += v[p][q];
       n++;
     }
     axis_rho[sample_index][0] = x/(axis_len*n);
     
     for (m=1;m<num_axis_samples-1;m++)
     if ((d=m*step)<=len)
     {
       while (axis_distance[k2]<d)
       {
         if (k2<num_axis_pts-1) k2++;
         else break;
       }
       k1 = k2-1;
       i1 = Iarray[0][k1]; j1 = Jarray[0][k1];
       i2 = Iarray[0][k2]; j2 = Jarray[0][k2];
       d1 = axis_distance[k1];
       d2 = axis_distance[k2];
       i = (i2*(d-d1)+i1*(d2-d))/(d2-d1);
       j = (j2*(d-d1)+j1*(d2-d))/(d2-d1);
       x = 0;
       n=0;
       for (p=i-step/2;p<i+step/2+1;p++)
       for (q=j-step/2;q<j+step/2+1;q++)
       if (pic[p][q]<200) if (axis[p][q]==1)
       {
         x += v[p][q];
         n++;
       }
       if (n>0) axis_rho[sample_index][m] = x/(axis_len*n);
     }
     m = num_axis_samples - 1;
     i = Ictr_gn; j = Jctr_gn;
     x = 0;
     n=0;
     for (p=i-step/2;p<i+step/2+1;p++)
     for (q=j-step/2;q<j+step/2+1;q++)
     if (pic[p][q]<200) if (axis[p][q]==1)
     {
       x += v[p][q];
       n++;
     }
     axis_rho[sample_index][m] = x/(axis_len*n);
     
     
     
     
     /* sample the splenium axis */
     len = spl_distance[num_spl_pts-1];
     actual_num_samples = len/step + 1;      
     if (actual_num_samples>num_spl_samples)
        actual_num_samples = num_spl_samples;
     k2 = 0;    
     i = Ictr_spl; j = Jctr_spl;
     x = v[i][j];
     n=1;
     for (p=i-step/2;p<i+step/2+1;p++)
     for (q=j-step/2;q<j+step/2+1;q++)
     if (pic[p][q]<200) if (axis[p][q]==2)
     {
       x += v[p][q];
       n++;
     }
     spl_rho[sample_index][0] = x/(axis_len*n);

     for (m=1;m<actual_num_samples;m++)
     if ((d=m*step)<=len)
     {
       while (spl_distance[k2]<d)
       {
         if (k2<num_spl_pts-1) k2++;
         else break;
       }
       k1 = k2-1;
       i1 = Iarr[k1]; j1 = Jarr[k1];
       i2 = Iarr[k2]; j2 = Jarr[k2];
       d1 = spl_distance[k1];
       d2 = spl_distance[k2];
       i = (i2*(d-d1)+i1*(d2-d))/(d2-d1);
       j = (j2*(d-d1)+j1*(d2-d))/(d2-d1);
       x = 0;
       n=0;

       for (p=i-step/2;p<i+step/2+1;p++)
       for (q=j-step/2;q<j+step/2+1;q++)
       if (pic[p][q]<200) if (axis[p][q]==2)
       {
         x += v[p][q];
         n++;       
       }    
       if (n>0) spl_rho[sample_index][m] = x/(axis_len*n);
     }
     
     
     
     
     
     /* sample the hk axis */
     len = hk_distance[num_hk_pts-1];
     actual_num_samples = len/step + 1;
     
     if (actual_num_samples>num_hk_samples)
        actual_num_samples = num_hk_samples;
    
     k2 = 0;
     
     i = Ictr_gn; j = Jctr_gn;
     x = v[i][j];
     n=1;
     for (p=i-step/2;p<i+step/2+1;p++)
     for (q=j-step/2;q<j+step/2+1;q++)
     if (pic[p][q]<200) if (axis[p][q]==2)
     {
       x += v[p][q];
       n++;
     }
     hk_rho[sample_index][0] = x/(axis_len*n);

     for (m=1;m<actual_num_samples;m++)
     if ((d=m*step)<=len)
     {
       while (hk_distance[k2]<d)
       {
         if (k2<num_hk_pts-1) k2++;
         else break;
       }
       k1 = k2-1;
       i1 = Iarray[1][k1]; j1 = Jarray[1][k1];
       i2 = Iarray[1][k2]; j2 = Jarray[1][k2];
       d1 = hk_distance[k1];
       d2 = hk_distance[k2];
       i = (i2*(d-d1)+i1*(d2-d))/(d2-d1);
       j = (j2*(d-d1)+j1*(d2-d))/(d2-d1);
       x = 0;
       n=0;
       for (p=i-step/2;p<i+step/2+1;p++)
       for (q=j-step/2;q<j+step/2+1;q++)
       if (pic[p][q]<200) if (axis[p][q]==2)
       {
         x += v[p][q];
         n++;
       }
       if (n>0) hk_rho[sample_index][m] = x/(axis_len*n);
     }
  }
  
  
  
  void calc_CCareas()
  /* calc the residual areas of genu and splenium
  after removing the ribbon */
  {
     int i,j,m,n,R,k,jmin,jmax,imin;
     float xs,xg,d;
     float dist[201][201],area;
     
     /* 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);
     }
     
     /* calculate the total area of the CC */
     /* tmp = 0 inside CC, =255 outside  */
     area = 0;
     jmin = xdim; jmax = 0;
     for (i=0;i<ydim;i++) for (j=0;j<xdim;j++) 
     {
       if (pic[i][j]<200)
       {
         area++;
         tmp[i][j] = 0;
         if (j<jmin) jmin=j;
         if (j>jmax) jmax=j;
       }
       else tmp[i][j]=255;
     }
     
     
     /* splenium */
     
     /* first mark discs centered along the main axis with (tmp=254) 
        in the first quarter */
     k = jmin+(jmax-jmin)/4; 
     for (i=0;i<ydim;i++) for (j=jmin;j<k;j++)
     if (axis[i][j]==1)
     {
       d = v[i][j]+1;
       R = d+0.5;
       for (m=i-R;m<i+R+1;m++) for (n=j-R;n<j+R+1;n++)
       if (tmp[m][n]==0) if (dist[m-i+100][n-j+100]<=d)
       tmp[m][n]=254;
     }
     
 
    xs=0;
    k = Jctr_spl+v[Ictr_spl][Jctr_spl]+1.5;
    for (i=Ictr_spl;i<ydim;i++) for (j=jmin;j<k;j++)
     if (tmp[i][j]==0)
     {
       if (bdry[i][j]>0)
         for (m=i-2;m<i+3;m++) for (n=j-2;n<j+3;n++)
           if (tmp[m][n]==254) tmp[i][j]=255;
       if (tmp[i][j]==0)
       {
         xs++;
         tmp[i][j]=1;
       }  
     }
     imin =Ictr_spl - v[Ictr_spl][Jctr_spl]/2;
     for (i=imin;i<ydim;i++) for (j=jmin;j<Jctr_spl;j++)
     if (tmp[i][j]==0)
     {
       if (bdry[i][j]>0)
         for (m=i-2;m<i+3;m++) for (n=j-2;n<j+3;n++)
           if (tmp[m][n]==254) tmp[i][j]=255;
           
       if (tmp[i][j]==0)
       {
         xs++;
         tmp[i][j]=1;
       } 
     }
     d = v[i=Ictr_spl][j=Jctr_spl];
     for (m=i-d-1.5;m<i+d+2.5;m++) for (n=j-d-1.5;n<j+d+2.5;n++)
     if (pic[m][n]<200) if (tmp[m][n]!=1)
     if (dist[m-i+100][n-j+100]<=d+1) {xs++;tmp[m][n]=1;}

     
     /* genu */
     /* mark pts on disks centered on the main axis (tmp=254) */
     imin=ydim; 
     for (i=0;i<ydim;i++) for (j=hk_j-10;j<jmax+1;j++)
     if (axis[i][j]==1)
     {
       d = v[i][j]+1;
       if (i+d<imin) imin=i+d;
       R = d+0.5;
       for (m=i-R;m<i+R+1;m++) for (n=j-R;n<j