ccdist.c

Computing 2D Skeleton

/* code for computing absolute distance function of a 2D shape */
/* the positon of the initial front and its position as it advances are
marked as points on grid lines */

/* input: image: piecewise constant /*

/* The boundary of the shape generically crosses gridlines
   between nodes */

/* output: absolute distance function  */


#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <time.h>
#include <fcntl.h>
#include <unistd.h>
#include <sioux.h>
#include <MacTypes.h>
#include <Fonts.h>
#include <math.h>

#define MAX(A,B) (A>B) ? (A) : (B)
#define MIN(A,B) (A<B) ? (A) : (B)
#define ABS(A) (A>=0) ? (A) : (-A)

/* imagesizes:
   ccSchzSchz:  376H x 796W; sample count: 35
   ccSchzNorms: 364H x 784W; sample count: 48
   ccNorms: 324H x 780W; sample count: 20 
   ccKids : 400H x 896W; sample count: 30
   ccDrugs: 324H x 708W; sample count: 14
 */

#define sample_count 14

#define ydim 324 /* image height - dimension for i */
#define xdim 708 /* image width - dimension for j */

#define MaxDist 10000

typedef struct element
{
	float dist;
	int row;
	int col;
} HeapNode;

FILE *infile,*outfile;
char inname[128],outname[128],dir_name[128];


/* row and col are indices Tmat, not in pic */

unsigned char pic[ydim][xdim],outpic[ydim][xdim];


int length, heaplength;
  
int IArray[ydim*xdim],JArray[ydim*xdim];
   /* IArray and JArray store the coordinates of 
   the bdry nodes the shape */

  unsigned char Cmat[ydim][xdim];
  int Bmat[ydim][xdim];
  float dist[ydim][xdim],Tmat[ydim+2][xdim+2];

/* Cmat marks 'far' points as 3, 'close' or 'trial points' as 2
and 'accepted pts' as 1.
far points = points not yet visited,
trial points: points next to the front and outside the front.
accepted pts = points inside the front.  */

/* Bmat stores backpointer from the grid points to their
index in the heap  */

/* Tmat stores d-values (=time when the front arrives at the node) */

  
  HeapNode H[ydim*xdim],TmpNode[1],Head[1];
  /* Note that H,TmpNode and Head are pointers to HeapNode */
/* array for storing the heap. Note that heap index = array index +1 */

  

void calc_dist();
void upheap(HeapNode *,int, HeapNode *);
void downheap(HeapNode *,int);
void locateboundary();
float upwinddist(int,int); 
void calc_dist();
void output_results();
void run();

void calc_dist()
{
  int i,j;
  int Iind, Jind,NIind,NJind,pos;
  /*Initialize */
  
    for(i=0;i<ydim;i++) for(j=0;j<xdim;j++) Cmat[i][j]=3;

    for(i=0;i<ydim+2;i++) for(j=0;j<xdim+2;j++) Tmat[i][j] = MaxDist;
    
    /*find the points neighboring the bdry and
    calculate their distances, store them in Tmat */
    /* bdry points are the grid nodes next to the initial front
    which generically crosses the grid lines between grid nodes */
  
    length = 0;

    locateboundary();
    
    /* Begin advancing the front */
    
    /* first, tag bdry points as accepted points */
    for(i=0;i<length;i++)
    {
      Iind = IArray[i];
      Jind = JArray[i];
      Cmat[Iind-1][Jind-1] = 1; /* point part of the front */
    }

    /* Start of the heap:  neighbors of bdry pts are initial trial points */
      heaplength = 0; /* number of elements on the heap */
      for(i=0;i<length;i++)
      {
      Iind = IArray[i];
      Jind = JArray[i];
      /*West Neighbor*/
      NIind = Iind;
      NJind = Jind-1;
      if(NJind > 0)
	if (Cmat[NIind-1][NJind-1] == 3)
	/*  its d-value is unassigned */
	  {
	    TmpNode->dist = upwinddist(NIind,NJind);
	    TmpNode->row = NIind;
	    TmpNode->col = NJind;
	    Cmat[NIind-1][NJind-1] = 2; /* a new trial point */
	    upheap(H,heaplength,TmpNode);
	    heaplength++;
	  }

      /*North Neighbor*/
      NIind = Iind-1;
      NJind = Jind;
      if(NIind > 0)
	if (Cmat[NIind-1][NJind-1] == 3)
	  {
	    TmpNode->dist = upwinddist(NIind,NJind);
	    TmpNode->row = NIind;
	    TmpNode->col = NJind;
	    Cmat[NIind-1][NJind-1] = 2;
	    upheap(H,heaplength,TmpNode);
	    heaplength++;
	  }


      /*East Neighbor*/
      NIind = Iind;
      NJind = Jind+1;
      if(NJind <= xdim)
      if (Cmat[NIind-1][NJind-1] == 3)
	  {
	    TmpNode->dist = upwinddist(NIind,NJind);
	    TmpNode->row = NIind;
	    TmpNode->col = NJind;
	    Cmat[NIind-1][NJind-1] = 2;
	    upheap(H,heaplength,TmpNode);
	    heaplength++;
	  }

      
      /*South Neighbor*/
      NIind = Iind+1;
      NJind = Jind;
      if(NIind <= ydim)
      if (Cmat[NIind-1][NJind-1] == 3)
	  {
	    TmpNode->dist = upwinddist(NIind,NJind);
	    TmpNode->row = NIind;
	    TmpNode->col = NJind;
	    Cmat[NIind-1][NJind-1] = 2;
	    upheap(H,heaplength,TmpNode);
	    heaplength++;
	  }

    }

   /*march the front forward */
   while(heaplength>0)
   {
     *Head = *H;
     downheap(H,heaplength);
     heaplength--;
     Iind = Head->row;
     Jind = Head->col;
     Cmat[Iind-1][Jind-1] = 1;
     Tmat[Iind][Jind] = Head->dist;

     /*West Neighbor*/
     NIind = Iind;
     NJind = Jind-1;
     if((NJind>0) && (Cmat[NIind-1][NJind-1]>1))
     {
	   TmpNode->dist = upwinddist(NIind,NJind);
	   TmpNode->row = NIind;
	   TmpNode->col = NJind;
	   /*Update the heap*/
	 if (Cmat[NIind-1][NJind-1]==2) /*already a trial point*/
	   {
	     pos = Bmat[NIind-1][NJind-1]; /* heap index of the point */
	     upheap(H,pos-1,TmpNode); /* pos-1 is the array index of the point */
	   }
	 else
	   {
	     Cmat[NIind-1][NJind-1] = 2;
	     upheap(H,heaplength,TmpNode);
	     heaplength++;
	   }  
     }
     /*North Neighbor*/
     NIind = Iind-1;
     NJind = Jind;
     if((NIind>0) && (Cmat[NIind-1][NJind-1]>1))
     {
	   TmpNode->dist = upwinddist(NIind,NJind);
	   TmpNode->row = NIind;
	   TmpNode->col = NJind;
	   
	 if (Cmat[NIind-1][NJind-1]==2) 
	   {
	     pos = Bmat[NIind-1][NJind-1];
	     upheap(H,pos-1,TmpNode);
	   }
	 else
	   {
	     Cmat[NIind-1][NJind-1] = 2;
	     upheap(H,heaplength,TmpNode);
	     heaplength++;
	   }
     }
     
     /*East Neighbor*/
     NIind = Iind;
     NJind = Jind+1;
     if((NJind<=xdim) && (Cmat[NIind-1][NJind-1]>1))
     {
	   TmpNode->dist = upwinddist(NIind,NJind);
	   TmpNode->row = NIind;
	   TmpNode->col = NJind;
	   
	 if (Cmat[NIind-1][NJind-1]==2) 
	   {
	     pos = Bmat[NIind-1][NJind-1];
	     upheap(H,pos-1,TmpNode);
	   }
	 else
	   {
	     Cmat[NIind-1][NJind-1] = 2;
	     upheap(H,heaplength,TmpNode);
	     heaplength++;
	   }
     }
     
     /*South Neighbor*/
     NIind = Iind+1;
     NJind = Jind;
     if((NIind<=ydim) && (Cmat[NIind-1][NJind-1]>1))
     {
	   TmpNode->dist = upwinddist(NIind,NJind);
	   TmpNode->row = NIind;
	   TmpNode->col = NJind;
	   
	 if (Cmat[NIind-1][NJind-1]==2) 
	   {
	     pos = Bmat[NIind-1][NJind-1];
	     upheap(H,pos-1,TmpNode);
	   }
	 else
	   {
	     Cmat[NIind-1][NJind-1] = 2;
	     upheap(H,heaplength,TmpNode);
	     heaplength++;
	   }
     }
   }

   for(i=0;i<ydim;i++) for(j=0;j<xdim;j++)
   {
      if (Tmat[i+1][j+1]<MaxDist)
      dist[i][j] = Tmat[i+1][j+1];
      else printf("(%d %d) unassigned\n",i,j);
   }
}