ccdist.c

Computing 2D Skeleton

void locateboundary() 
{
  /* Locate the pts neighboring the actual bdry.
  Calculate their distance from the bdry. */
 
  int NIind, NJind,i,j;
  int  VP,VW,VN,VE,VS,SW,SN,SE,SS;
  float d;
  
  for(i=2;i<=ydim-1;i++) for(j=2;j<=xdim-1;j++)
  /* indices (i,j) belong to array Tmat, not pic */
  {
	VP = pic[i-1][j-1];
	/*West Neighbor*/
    NIind = i;
	NJind = j-1;
	VW = pic[NIind-1][NJind-1];
	/*North Neighbor*/
	NIind = i-1;
	NJind = j;
	VN = pic[NIind-1][NJind-1];
	/*East Neighbor*/
	NIind = i;
	NJind = j+1;
	VE = pic[NIind-1][NJind-1];
	/*South Neighbor*/
	NIind = i+1;
	NJind = j;
	VS = pic[NIind-1][NJind-1];

	SW = SN = SE = SS = 0;
    if(VW!=VP) SW = 1;
    if(VN!=VP) SN = 1;
	if(VE!=VP) SE = 1;
	if(VS!=VP) SS = 1;
	
	if (SW+SE+SN+SS==0) d = 0;
	else
	{
	  IArray[length] = i;
	  JArray[length] = j;
	  length++;
	  if ((SW+SE==2)||(SN+SS==2)||(SW+SE+SN+SS==1))
		d = 0.5;
	  else d = 0.35355;
	  Tmat[i][j] = d;
	}
  }
}

float upwinddist(int Iind,int Jind)
{
  /*Calculate the value at (Iind,Jind) using the upwind scheme*/
  float v1,v2,d;
  double w1,w2;

  v1 = MIN(Tmat[Iind][Jind-1],Tmat[Iind][Jind+1]);
  v2 = MIN(Tmat[Iind-1][Jind],Tmat[Iind+1][Jind]);
  w1 = MAX(v1,v2);
  w2 = MIN(v1,v2);
  if(w1-w2>=1) d = w2+1;
  else d = (w1+w2+sqrt(2-(w1-w2)*(w1-w2)))/2;

  return d;			
}



void upheap(HeapNode *init, int length,HeapNode *newnode)
{
    /* init points to the root node. length = heap array index of the hole.
    heap index numbering starts at 1. (i.e. the numbering of tree nodes.)
    The index of the heap array = heap index - 1 */
    
    /* For a heap structure, two objects are needed to describe it:
    the array of HeapNodes and the length of the array. 
    length indicates where the upheap should begin, which means that at 
	length+1 the heap array has hole and the upheap should begin there.

    The change of the number of nodes in the heap is maintained separately.
    We should not depend on either upheap or downheap to update it. So that
    the two function can be used flexibly.
	*/
	int flag,index,parentindex,Iind,Jind;

	index = length+1; /* heap index of the hole */

    flag = 1;

	do{
		parentindex =  (int)floor(index/2);
		if(parentindex > 0)
		{
			if(newnode->dist<(init+parentindex-1)->dist)
			/* move the parent down */
			{
				*(init+index-1) = *(init+parentindex-1);
				Iind = (init+parentindex-1)->row;
				Jind = (init+parentindex-1)->col;
				Bmat[Iind-1][Jind-1] = index; /* update backpointer */
				index = parentindex; /* new location of the hole */
			}
			else /* insert the newnode */
			{
				*(init+index-1) = *newnode;
				Iind = newnode->row;
				Jind = newnode->col;
				Bmat[Iind-1][Jind-1] = index;
				flag = 0;
			}
		}
		else /* arrived at the root node */
		{
			*(init+index-1) = *newnode;
			Iind = newnode->row;
			Jind = newnode->col;
			Bmat[Iind-1][Jind-1] = index;
			flag = 0;
		}

	} while(flag == 1);

}

void downheap(HeapNode *init,int length)
/* downheap is called after removing the root node. */
{
	int  index,Iind,Jind;
     
	index=1; /* heap index of the hole */

	while(length>=2*index) /* not a leaf node */
	{
		if(length>=(2*index+1)) /* both children present */
		{
			if((init+2*index-1)->dist <= (init+2*index)->dist)
			/* left child smaller than right child, promote left */
			{
				*(init+index-1) = *(init+2*index-1);
				Iind = (init+2*index-1)->row;
				Jind = (init+2*index-1)->col;
				Bmat[Iind-1][Jind-1] = index; /* update backpointer */
				index = 2*index; /* new location of the hole */
			}
			else /* promote right child */
			{
				*(init+index-1) = *(init+2*index);
				Iind = (init+2*index)->row;
				Jind = (init+2*index)->col;
				Bmat[Iind-1][Jind-1] = index;
				index = 2*index+1;
			}
		}
		else /* only one child present, always on the left */
		{
				*(init+index-1) = *(init+2*index-1);
				Iind = (init+2*index-1)->row;
				Jind = (init+2*index-1)->col;
				Bmat[Iind-1][Jind-1] = index;
				index = 2*index;
		}
	}

	if (index < length)
	/* hole at a leaf which is not the last leaf.
	   insert the last leaf node here and then migrate upward  */
	   upheap(init,index-1,init+length-1);
}


 void output_results()
 {
    /*
    int i,j;
    float x,mx=0,mn=MaxDist+1,scale;

    for (i=0;i<ydim;i++) for (j=0;j<xdim;j++)
    {
       x = dist[i][j];
       mx = MAX(mx,x);
       mn = MIN(mn,x);
    }
    x = mx-mn;
    if (x>1) scale=255/x;
    else scale = 1;
    for (i=0;i<ydim;i++) for (j=0;j<xdim;j++)
      outpic[i][j]=(int) ((dist[i][j]-mn)*scale);

     printf("Writing %s\n",outname);
     if ((outfile = fopen(outname,"w+")) == NULL)
     {
	   printf("Cannot create %s\n",outname);
	   exit(1);
     }
     fwrite(outpic,sizeof(unsigned char),xdim*ydim,outfile);
     fclose(outfile); 
     */
    
     
     printf("Writing %s\n",outname);
     if ((outfile = fopen(outname,"w+")) == NULL)
     {
	   printf("Cannot create %s\n",outname);
	   exit(1);
     }
     fwrite(dist,sizeof(float),xdim*ydim,outfile); 
     fclose(outfile); 
 } 
 



void run(void)
{
   int m,n;
   char str[100][128];
   char sstr[100][128];
   time_t systime;
   	   
   SIOUXSettings.asktosaveonclose = 0;
   systime = time(NULL);
   puts(ctime(&systime));
   
  /* FILE NAMES */
 /* Strings str are input files; sstr are output files */
 
       
   /* SCHIZOPHRENICS */
  /*
   strcpy(dir_name,"::CorpusCallosum:ccSchzSchz:");

   strcpy(str[0],"::CorpusCallosum:ccSchzSchz:1030_2_128.buchar");
   strcpy(str[1],"::CorpusCallosum:ccSchzSchz:1033_2_128.buchar");
   strcpy(str[2],"::CorpusCallosum:ccSchzSchz:1035_2_128.buchar");
   strcpy(str[3],"::CorpusCallosum:ccSchzSchz:1036_2_128.buchar");
   strcpy(str[4],"::CorpusCallosum:ccSchzSchz:1037_2_128.buchar");
   strcpy(str[5],"::CorpusCallosum:ccSchzSchz:1038_4_128.buchar");
   strcpy(str[6],"::CorpusCallosum:ccSchzSchz:468_6_128.buchar");
   strcpy(str[7],"::CorpusCallosum:ccSchzSchz:474_4_128.buchar");
   strcpy(str[8],"::CorpusCallosum:ccSchzSchz:476_6_128.buchar");
   strcpy(str[9],"::CorpusCallosum:ccSchzSchz:479_6_128.buchar");
   strcpy(str[10],"::CorpusCallosum:ccSchzSchz:491_4_128.buchar");
   strcpy(str[11],"::CorpusCallosum:ccSchzSchz:510_4_128.buchar");
   strcpy(str[12],"::CorpusCallosum:ccSchzSchz:511_4_128.buchar");
   strcpy(str[13],"::CorpusCallosum:ccSchzSchz:519_4_128.buchar");
   strcpy(str[14],"::CorpusCallosum:ccSchzSchz:578_4_128.buchar");
   strcpy(str[15],"::CorpusCallosum:ccSchzSchz:579_4_128.buchar");
   strcpy(str[16],"::CorpusCallosum:ccSchzSchz:580_4_128.buchar");
   strcpy(str[17],"::CorpusCallosum:ccSchzSchz:597_4_128.buchar");
   strcpy(str[18],"::CorpusCallosum:ccSchzSchz:614_4_128.buchar");
   strcpy(str[19],"::CorpusCallosum:ccSchzSchz:615_4_128.buchar");
   strcpy(str[20],"::CorpusCallosum:ccSchzSchz:631_4_128.buchar");
   strcpy(str[21],"::CorpusCallosum:ccSchzSchz:635_5_128.buchar");
   strcpy(str[22],"::CorpusCallosum:ccSchzSchz:647_4_128.buchar");
   strcpy(str[23],"::CorpusCallosum:ccSchzSchz:672_4_128.buchar");
   strcpy(str[24],"::CorpusCallosum:ccSchzSchz:677_4_128.buchar");
   strcpy(str[25],"::CorpusCallosum:ccSchzSchz:686_4_128.buchar");
   strcpy(str[26],"::CorpusCallosum:ccSchzSchz:689_4_128.buchar");
   strcpy(str[27],"::CorpusCallosum:ccSchzSchz:690_4_128.buchar");
   strcpy(str[28],"::CorpusCallosum:ccSchzSchz:772_4_128.buchar");
   strcpy(str[29],"::CorpusCallosum:ccSchzSchz:789_4_128.buchar");
   strcpy(str[30],"::CorpusCallosum:ccSchzSchz:790_4_128.buchar");
   strcpy(str[31],"::CorpusCallosum:ccSchzSchz:791_4_128.buchar");
   strcpy(str[32],"::CorpusCallosum:ccSchzSchz:792_4_128.buchar");
   strcpy(str[33],"::CorpusCallosum:ccSchzSchz:793_4_128.buchar");
   strcpy(str[34],"::CorpusCallosum:ccSchzSchz:853_4_128.buchar");
   
   strcpy(sstr[0],"::CorpusCallosum:ccSchzSchz:1030");
   strcpy(sstr[1],"::CorpusCallosum:ccSchzSchz:1033");
   strcpy(sstr[2],"::CorpusCallosum:ccSchzSchz:1035");
   strcpy(sstr[3],"::CorpusCallosum:ccSchzSchz:1036");
   strcpy(sstr[4],"::CorpusCallosum:ccSchzSchz:1037");
   strcpy(sstr[5],"::CorpusCallosum:ccSchzSchz:1038");
   strcpy(sstr[6],"::CorpusCallosum:ccSchzSchz:468");
   strcpy(sstr[7],"::CorpusCallosum:ccSchzSchz:474");
   strcpy(sstr[8],"::CorpusCallosum:ccSchzSchz:476");
   strcpy(sstr[9],"::CorpusCallosum:ccSchzSchz:479");
   strcpy(sstr[10],"::CorpusCallosum:ccSchzSchz:491");
   strcpy(sstr[11],"::CorpusCallosum:ccSchzSchz:510");
   strcpy(sstr[12],"::CorpusCallosum:ccSchzSchz:511");
   strcpy(sstr[13],"::CorpusCallosum:ccSchzSchz:519");
   strcpy(sstr[14],"::CorpusCallosum:ccSchzSchz:578");
   strcpy(sstr[15],"::CorpusCallosum:ccSchzSchz:579");
   strcpy(sstr[16],"::CorpusCallosum:ccSchzSchz:580");
   strcpy(sstr[17],"::CorpusCallosum:ccSchzSchz:597");
   strcpy(sstr[18],"::CorpusCallosum:ccSchzSchz:614");
   strcpy(sstr[19],"::CorpusCallosum:ccSchzSchz:615");
   strcpy(sstr[20],"::CorpusCallosum:ccSchzSchz:631");
   strcpy(sstr[21],"::CorpusCallosum:ccSchzSchz:635");
   strcpy(sstr[22],"::CorpusCallosum:ccSchzSchz:647");
   strcpy(sstr[23],"::CorpusCallosum:ccSchzSchz:672");
   strcpy(sstr[24],"::CorpusCallosum:ccSchzSchz:677");
   strcpy(sstr[25],"::CorpusCallosum:ccSchzSchz:686");
   strcpy(sstr[26],"::CorpusCallosum:ccSchzSchz:689");
   strcpy(sstr[27],"::CorpusCallosum:ccSchzSchz:690");
   strcpy(sstr[28],"::CorpusCallosum:ccSchzSchz:772");
   strcpy(sstr[29],"::CorpusCallosum:ccSchzSchz:789");
   strcpy(sstr[30],"::CorpusCallosum:ccSchzSchz:790");
   strcpy(sstr[31],"::CorpusCallosum:ccSchzSchz:791");
   strcpy(sstr[32],"::CorpusCallosum:ccSchzSchz:792");
   strcpy(sstr[33],"::CorpusCallosum:ccSchzSchz:793");
   strcpy(sstr[34],"::CorpusCallosum:ccSchzSchz:853");