ccctrs.c

Computing 2D Skeleton

  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 find_ctrs()
  {
     int i,j,m,n;
     int clr,mxflag,mnflag,endflag;
     int k=2,kk=3;
     float x,y;
     
     for (i=0;i<ydim;i++) for (j=0;j<xdim;j++) ctrs[i][j]=0;
     
     for (i=0;i<ydim;i++) for (j=0;j<xdim;j++) 
     if (sk[i][j]>0)
     {
        clr = pic[i][j];
       
        endflag=0;
        for (m=i-k;m<i+k+1;m++) for (n=j-k;n<j+k+1;n++)
        if (pic[m][n]==clr)  if (endpts[m][n]>0) endflag=1;
        
        x = v[i][j];        
        mxflag=1;
        if ((endflag==0)&&(sk[i][j]==1)) mnflag=1;
        else mnflag=0;
        
        for (m=i-kk;m<i+kk+1;m++) for (n=j-kk;n<j+kk+1;n++)
        if (pic[m][n]==clr) if (sk[m][n]>0)
        {
           y = v[m][n];
          if (y>x) mxflag=0;
          if (mnflag>0) if (sk[m][n]==1) if (y<x) mnflag=0;
        }
        
        if (mxflag>0) if (confirm_ctr(i,j)>0)
        for (m=i-2;m<i+3;m++) for (n=j-2;n<j+3;n++) ctrs[m][n]=1;
        
        if (mnflag>0) if (confirm_saddlept(i,j)>0)
        for (m=i-3;m<i+4;m++) for (n=j-3;n<j+4;n++) ctrs[m][n]=2;
     }
  }
 
  
  
  
  int confirm_ctr(int i,int j)
  {
     int m,n,clr,mxflag;
     float x;
     
     mxflag=1;
     clr=pic[i][j];
     x = v[i][j];
     for (m=i-NN;m<i+NN+1;m++) for (n=j-NN;n<j+NN+1;n++)
     if (pic[m][n]==clr) if (sk[m][n]>0)
     if (v[m][n]>x) mxflag=0;
     return(mxflag);
  }
  
  
  
  
  int confirm_saddlept(int i,int j)
  {
    int m,n,mm,nn,clr,initial,mnflag,flag;
    int k=2;
    float x;

    initial = flag = mnflag = 0;
    clr = pic[i][j];
    
    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))
    if (pic[m][n]==clr) if (sk[m][n]==1)
    {
       if (initial==0) {mm=m; nn=n; initial=1;}
       else if ((m-i)*(mm-i)+(n-j)*(nn-j)<0) /* obtuse angle */
       flag = 1; /* not an endpoint */
    }
    
    if (flag>0)
    {
      mnflag=1;
      x = v[i][j];
      for (m=i-NN;m<i+NN+1;m++) for (n=j-NN;n<j+NN+1;n++)
      if (pic[m][n]==clr) if (sk[m][n]==1)
      if (v[m][n]<x) mnflag=0; 
    }  
    return(mnflag);
  }
  
  

 
  void connect(m,n,mm,nn,color)
  int m,n,mm,nn,color;
  /* Draw straight line conncting (m,n) to (mm,nn) */
 {       
    int i,j,imin,imax,jmin,jmax;
    float x,d,d0,dx,dy;

    
    imin=MIN(m,mm); imax=MAX(m,mm);
    jmin=MIN(n,nn); jmax=MAX(n,nn);
    
 /*   printf("imin=%d imax=%d jmin=%d jmax=%d",imin,imax,jmin,jmax); */
       
    dx = nn-n; dy = mm-m; d0 = m*dx-n*dy;
    if ((x=dx*dx+dy*dy) != 0)  
    for (i=imin;i<imax+1;i++) for (j=jmin;j<jmax+1;j++)
    {      
       d = i*dx-j*dy - d0;
       d *= d;
       d /= x; /* sq-distance from the st. line */       
       /* mark the point if it is within 1 pixel from the st. line */
        
       if (d<1) if (pic[i][j]==color)
       {             
          sgmt[i][j]=1;
/* if (i==182) if ((j==171)||(j==179))
 {
printf("sgmt pt (%d %d) connecting (%d %d) (%d %d)\n",i,j,m,n,mm,nn);
} */             
       }       
    }
      
 }
 
  
  
  void segment()
  /* segments protusions */
  {
    int i,j,m,n,k0,k,kk,k1,p,q;
    int t1,t2,i1,j1,i2,j2;
    int mk,h,H,flag;
    int ii[4*N-4],jj[4*N-4],check[4*N-4];
    float v0,z,zm,d[4*N-4];
    float slope[4*N-4],last,now,next;
    float grx,grmn,accept;

    for (i=0;i<ydim;i++) for (j=0;j<xdim;j++)
    if (endpts[i][j]==21) if (length[0][mk=mark[0][i][j]]>0)
    {     
      h = nbhd_size[i][j];     
/* if (i==171) if (j==300) */
/*printf("(%d %d), mk=%d, h=%d\n",i,j,mk,h);*/

      if (h>1)
      {
        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) */
      accept=(h-1)/sqrt((double) h*h+0.25); 
       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]; 
       }
       
       grmn = 1;
       
       /* 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 (mk==109)*/
/* if (i==358) if (j==112) */
/*  { printf("\n at (%d %d), nbhd: %d\n",i,j,h);
   for (k=0;k<H;k++) if (check[k]>0)
printf("k:%d (%d %d),slope:%f,chk:%d\n",k,ii[k],jj[k],slope[k],check[k]);} */
       
       flag = 0;
       if (z>accept) /* first ray to shape bdry */
       {
        k0 = kk;
        for (k=kk;k<kk+H;k++)
        if ((k%H)==k0)
        {
          /* find next ray to shape bdry */
          for (q=1;q<H;q++) 
          {
            flag = 0;
            if (check[(k+q)%H]>0) /* is local max */
            if (slope[(k+q)%H]>accept)  
            {
              flag = 1;
              k1 = (k+q)%H; /* next ray */
              break;
            }
          }
          
          if (flag==0) if ((k%H)==kk)
          /* no slope above accept; just find the next max slope */
          {
            now=0;
            for (q=1;q<H;q++)